CC BY
76
УДК 519.622.2; 681.5.015 М. Ю. Хавинсон, М. П. Кулаков
ПОДХОД К ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ СИСТЕМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
Рассматривается подход к оценке параметров систем обыкновенных дифференциальных уравнений с постоянными коэффициентами, позволяющий точечно оценить коэффициенты без нахождения аналитического решения системы. Приведены примеры реализации описываемого подхода в среде MathCAD и результаты теста на погрешность.
Ключевые слова: параметрическая идентификация систем, обыкновенные дифференциальные уравнения, точечная оценка коэффициентов.
Michael Yu. Khavinson, Matvey P. Kulakov AN APPROACH TO PARAMETRIC IDENTIFICATION OF ORDINARY DIFFERENTIAL EQUATIONS
(Institute for Complex Analysis of Regional Problems Far Eastern Branch of the Russian Academy of Sciences, Birobidzhan)
A new approach to parametric identification of ordinary differential equations with constant coefficients is considered in this article. The aррroach allows to estimate the coefficients without finding the analytical solutions of the system. The examples of taking this aррroach implemented in system MathCAD and results of the aptitude test are given.
Key words: parametric identification of systems, ordinary differential equations, point estimation of coefficients.
Введение
Важной задачей моделирования динамики систем различного характера является оценка коэффициентов моделей. В идентификации систем нет универсального метода, позволяющего оценить параметры любой модели [1—4, 9—12], поэтому каждый исследователь самостоятельно выбирает соответствующий его задаче математический и программный инструментарий. Большинство модельных уравнений не имеет аналитического решения, и общий метод идентификации систем должен включать,
Хавинсон Михаил Юрьевич — научный сотрудник (Институт комплексного анализа региональных проблем Дальневосточного отделения Российской академии наук, г. Биробиджан), e-mail: [email protected]
Кулаков Матвей Павлович — младший научный сотрудник (Институт комплексного анализа региональных проблем Дальневосточного отделения Российской академии наук, г. Биробиджан), e-mail: [email protected]
Работа частично поддержана грантами РГНФ проект № 13-12-79001 а (р), ДВО РАН № 13-Ш-В-10-001 и РФФИ № 11-01-98512-р_восток_а
© Хавинсон М.Ю., Кулаков М.П., 2013
113
во-первых, нахождение численного решения, во-вторых, оценку параметров моделей.
В настоящей статье рассматривается комбинированный алгоритм, позволяющий совместить нахождение приближенного решения системы дифференциальных уравнений и точечную оценку ее параметров по реальным данным, заданным в табличном виде, т. е. в виде дискретной функции на конечном временном интервале. Данный подход к оценке параметров системы дифференциальных уравнений позволяет идентифицировать систему, не находя и не аппроксимируя ее аналитическое решение. Предложенный алгоритм является комбинацией метода Рун-ге-Кутты 4 порядка с переменным шагом и градиентных методов оптимизации (в частности, Левенберга-Марквардта). Выбор указанных численных методов обусловлен их относительно высокой вычислительной скоростью в сочетании с эффективностью. Подобный подход к решению задачи параметрической идентификации в общей формулировке описан в работах В. А. Охорзина [7, 8], Б. П. Безручко и др. [3]. В настоящей работе развита идея этих подходов применительно к системам обыкновенных дифференциальных уравнений, произведена программная реализация и проведены тесты на погрешность.
В предлагаемом подходе функционал невязки сложно зависит от численного решения системы обыкновенных дифференциальных уравнений и реальных данных. При его минимизации из некоторой стартовой догадки или начальных значениях оцениваемых параметров производится численное решение системы, вычисляется невязка модельной и реальной динамики, определяется дальнейшее направление спуска и новое приближение параметров. На каждом из последующих шагов спуска процедура повторяется, пока итерационный процесс не сойдется, а значение целевой функции не достигнет некоторого минимума. Более того, для получения максимально правдоподобных и содержательно значимых оценок в алгоритме предлагается рассматривать множество стартовых догадок с многократным поиском минимумов, которых для нелинейных моделей может быть несколько. Для каждого из них вычисляются дополнительные показатели качества полученных оценок.
Алгоритм параметрической идентификации систем
Структура предложенного алгоритма позволяет использовать различные приближенные схемы решения дифференциальных уравнений и методы оптимизации. Общая схема подхода представлена на рис. 1.
Алгоритм состоит из 3 блоков. Блок 1 включает в себя задание системы обыкновенных дифференциальных уравнений вида:
х = Г| х, А|,
где X — вектор производных по времени t гладких определенных в пространстве Я" функций Х1, %ъ---, XI,..., хп, А — матрица действительных
114
параметров системы; Г х, А| е К" — вектор-функция правых частей системы дифференциальных уравнений.
Рис. 1. Блок-схема алгоритма оценки коэффициентов систем дифференциальных уравнений
Вторая часть блока 1 состоит из задания матриц табличных функций X,-=| (Х^,- (Х2), ••• (!'„), (/=1,2,..., Ы) — реальных данных объемом N измерений, а так же начальных условий х1 (¿0) системы и начальной догадки о значениях коэффициентов А0. В моделях начальными условиями, как правило, служат точки табличных функций, например,
первая точка реальных данных хI (^)=(Х,)0 (¿=1,2,___, и). Для получения
наиболее достоверных оценок параметров начальные условия могут быть оценены совместно с параметрами. Начальную догадку о значениях коэффициентов модели подбирает сам исследователь, исходя из содержательного анализа модели. В блоке 2 производится численное решение системы с учетом заданных начальных условий и начальных приближений параметров системы. В блоке 3 расчетные точки модели сравниваются с реальными данными. Блок 3 представляет собой оптимизацию функционала суммы квадратов отклонений модельных и фак-
115
тических значений. При удовлетворительном результате заданные параметры остаются, при неудовлетворительном меняется начальная догадка и происходит повтор блока 2. Качество оценки определяется критериями оптимизационного метода (достижение функционалом глобального минимума), а также дополнительными показателями качества приближения модельных данных к фактическим (минимальная ошибка аппроксимации и максимальный коэффициент корреляции).
Более подробная запись алгоритма для системы автономных дифференциальных уравнений выглядит следующим образом. Пусть математическая модель некоторого процесса имеет вид:
„ |
~~Г = Л аЪ — ^ а1■ ХЪ — ^ Хп
ш
~ХГ~ Ь1 ■ Ь2■ ■ • • ■ Ьк ■ Х1, Х2 ■ ■ • • ■ ХП.
Ш '
^Хп г |
—n = fn\W1■ ■ Х2 ■ ■ • • ■ Хп1
где Х2,---,, Хп — фазовые переменные, t — переменная времени, аъ, а2,..., ар, Ьъ, Ъъ--; Ьк,..., шъ, Ш2,..., Шз — параметры системы, /ъ(аъ аъ-.., а, Х1, Х2,...,, Хп), _/2(Ьъ, Ь2,., Ьк, Хъ, Х2,.,, Хп), /п (Ш1, Ш2,., Шз, Хъ, Х2,.,, Хп) — правые части системы уравнений. Для удобства записи определим, что ]>к>..>Б, тогда матрица А начальных приближений параметров системы имеет вид:
A =
b1 b2 b3
Vwi
w9
w,
0
0
и2 .....* " ••• " у
Согласно схеме Рунге-Кутты 4 порядка [2] имеем:
¿+1 =: xi' ,+Нки + 2К12 + 2К13 + ки\ /6,
2 ¿+1 =; х2| . + И К21 + 2.К22 + 2К23 + К241 /6
где
ХП \ / + 1 = I х„: г + h Кп! +2Кп2 + 2Кп3 + Кп41 /6,
Kn=hfl\al,a2,...aJ,xl,x2,...,x„\, Ки =hfl\al,a2,...aJ,xl+0.5Ku,x2+0.5K2l,...,xn+0.5Knl\, Ки = hf1\a1,a2,...aJ,x1+0.5K12,x2+0.5K22,...,x„+0.5K„2\, Ки =hflla1,a2,...aJ,x1 + 0.5 Кп,х2 + 0.5 К21,...,хп + 0.5Кп1\.
a
a
a
a
a
a
1
2
3
4
5
b
0
k
116
Положив в данную схему начальные догадки о значениях параметров Ао, начальные значения фазовых переменных (х;)о=х; (¿о)=(Х;)о (¿=1,2,..., п) и шаг интегрирования й=1 /г, превосходящий шаг реальных данных в г раз, несложно получить матрицу 8; = | (5"]••• 0!?и);1 (¿=1,2,..., Ыг) численного решения исходной системы, где (Б;0 I — матрица-столбец найденных численно значений Хгф, (Б2) ; — значений Х2(£), ..., (Бп) ; — значений Хп (¿) íо+h, ¿о+2й,...). Пусть определены табличные функции реальных
данных объемом N в виде матриц-столбцов (Х1) ¿, (Х2) ¿, ..., (Хп) ^ с шагом, по времени равным 1 (¿=1,2,., Тогда функционал суммы квадратов отклонений фактических и модельных значений, т. е. невязка, имеет вид:
5д1а1,Ь1,...,\г1,...,а,,Ьк,...\гх' = Хх1 г-I¿V г,г12 + //2£!| Х2| г-I 52] г,г.12 +
'■Г ' I £ I-г
¡=1 г=1
/-1
ГДе //, //2...../¡п — весовые коэффициенты, N — длина ряда данных.
Последний этап алгоритма — минимизация функционала Бц (а1, Ьу..., ш1,., а¡, Ьк,..., шБ) по параметрам системы. В представленном алгоритме (рис. 1) может быть использован любой метод оптимизации.
Реализация параметрической идентификации систем дифференциальных уравнений в MathCAD
При реализации описанного алгоритма в МаШСЛО, равно как и в другой системе компьютерной алгебры, есть ряд особенностей. В качестве примера рассмотрим идентификацию параметров системы, представленной в работе [10]:
сЖх .
— = Ь - ах, Л
- = Сх- Еу, [Ж
где х, у — фазовые переменные, I — переменная времени, Ь, а, С, Е — параметры системы.
Описанный алгоритм можно начать с введения реальных данных. В блоке 1 иногда удобно записать значения табличной функции в виде транспонирования матрицы-строки, также есть возможность считывать данные из файла (например, с расширением или *.хк) с помощью встроенных команд. На рис. 2 матрицы Х1 и Х2 — реальные данные.
Параметры в блоке 1 удобно вводить в уравнения системы, как переменные, производные которых равны нулю. На рисунке 2 а в матрице-столбце Б (¿, х) скаляры хо и Х1 — фазовые переменные, Х2, хэ,..., Х5 — параметры системы. Начальные условия (Хо и Х1 на рис. 2 а) и стартовая догадка (а, Ь, С и Е на рис. 2 а) задаются в виде равенств, в случае многомерной системы с большим числом параметров их удобнее записать в виде матрицы.
117
XI : (4.3 4.6 4,7 4,9 5 5-21 XI (0.485 0.447 0.476 0.509 0.482 0.525 )Т ( Х2-ХЗ-Х0 ^
о о о
к О
хО:- XIи а 0.03 xl Х2|) С : 0.О4
Ь 0.09 Ii0.08
N - rows(Xl) N = 6
i := 0.. N - I г := 10
Ry<x,p) := Rtadap[[slack(x,p),0,N - 1 ,(N - l)r,D]
Sq(x0,xl,a,b,C,E)
раг«-{а b С El ini1<— (xO xl/ Tuitgc*— RytfiniLpar) K-l
Sql<-^ (XI-mnseir ij2
¡ = 0 N-l
(X2i-iwigei.,l2}3
¡ = 0 Sq I + Sq2
Sq(xO,xl,a,b,C,E) ^ 1.1.389 üiven
Sq(xO,xl,a,b,C,E) = 0
p .= M i г io: г! а, b,( . 1 ;
f0.917059343127297 ^ 0.155238340329286 0.869754219370492 ^8.656381938396079,
Sqfxio.Xlo.po.pbM.P.;) - 9.502« 10"
P =
ä)
ü)
Рис. 2. а) реализация блока 1, б) реализация блока 2 и построение невязки, в) минимизация невязки в MathCAD
В блоке 2 встроенная команда КкайарЬ (решать заданное уравнение методом Рунге-Кутты 4 порядка с переменным шагом) записывается как функция ЯуЬ (х, р) от начальной точки интегрирования х и догадки р (рис. 2 б). В приведенном примере шаг интегрирования равен Н=1/г=0,1. На рис. 2 б также представлен вариант составления функционала невязки Бц (х0, у0, а, Ь, С, Е) в среде МаШСАБ.
При нахождении минимума невязки желательно выводить на экран значения невязки до и после процедуры оптимизации (вывод значений до процедуры оптимизации позволит лучше подобрать начальные приближения, после — оценить суммарную погрешность) (рис. 2 в).
Для расчета качества приближения модели и реальных данных можно использовать среднюю ошибку аппроксимации (оценивает степень отклонения расчетных значений от фактических) и коэффициент корреляции (оценивает, насколько модель «улавливает» динамику изменения табличной функции) (рис. 3). Для построения графиков удобно брать большее число точек, чем длина ряда данных, в то время как для вычисления качества приближения необходимо учитывать только те модельные значения, которые по времени соответствуют реальным данным (для этого введена дополнительная функция Яу^(х, у) с шагом интегрирования й=1).
При решении задачи параметрической идентификации для систем дифференциальных уравнений описанным выше способом необходимо учитывать шаг интегрирования, количество точек табличных функций и начальные приближения оцениваемых коэффициентов. Несмотря на то, что алгоритм работает на основе переменного шага интегрирования (при этом решения выводятся через заданный равномерный шаг), необходимо самостоятельно задавать в программе оптимальный шаг для каждой конкретной системы: слишком мелкая сетка может существенно замедлить расчеты, слишком крупная — некорректно описать характер динамических режимов (особенно при наличии колебаний).
118
Рис. 3. Оценка приближения табличной функции и расчетных данных в MathCAD
Длина ряда данных, по которому оцениваются коэффициенты модели, по возможности, должна быть максимально длинной. Как правило, динамика фактической величины зашумлена. В случае рядов данных, число точек которых исчисляется сотнями, можно использовать возможности фильтрования, отбрасывания шума, выделения колебаний и т. д., но довольно многие экономические показатели состоят из 10-50 точек, что делает невозможным фильтрацию шума. В этом случае необходимо искать некоторый пучок решений, соответствующий разбросу данных, а затем с помощью экспертного анализа определять одно или несколько подходящих решений. Очевидно, что в этом случае, при наличии некоторой тенденции, чем длиннее ряд, тем более узким спектром кривых его можно приблизить.
Еще одним очень важным аспектом в применении описанного подхода параметрической идентификации является подбор начальных приближений коэффициентов или стартовой догадки. Как уже было отмечено, относительные небольшие ряды данных и наличие шума может привести к необходимости поиска пучка решений, т. е. поиска разных наборов коэффициентов. Иными словами, минимизируемый в приведенном алгоритме функционал может иметь несколько содержа-
119
тельно равноценных локальных минимумов, которые необходимо находить, варьируя начальные приближения параметров. Эти начальные приближения должен определять эксперт, учитывая особенности моделируемой ситуации и структуры модели.
Для оценки погрешности алгоритма рассмотрены три дифференциальных уравнения, решения которых можно записать аналитически:
х = Ь-ах, (1)
х = гх — жг2 (2)
х + 2кх + са2х = 0 (3)
Выбранная совокупность уравнений позволяет продемонстрировать работу подхода к параметрической идентификации линейных, нелинейных и имеющих периодические решения уравнений. Для каждого уравнения при конкретных значениях параметров задается дискретный ряд решений, полученных из аналитического решения. Затем, варьируя «опции» алгоритма, по представленному ряду данных решается задача параметрической идентификации.
В расчете погрешности описанного алгоритма оценки коэффициентов использован метод двойного пересчета. Перебирался шаг интегрирования й=0.1 и 0.05, количества точек Ы, по которым выполнялась оценка параметров, также варьировались начальных приближений оцениваемых параметров. Рассмотрены комбинации стартовых догадок двух параметров каждого уравнения, соответствующих одной десятой и десятикратному заданному значению каждого параметра (табл. 1, 2 и 3).
Таблица 1
0,1 я 10я
Aa, % Ab, % Acp, % Aa, % Ab, % Acp, %
N=5, h=0,1 0,716 0,426 0,005182 10,274 5,700 0,015737
0,1b
10b 0,352 0,550 0,007535 1,343 2,410 0,004531
N=5, h=0,05 0,716 0,426 0,005182 10,277 5,697 0,015741
0,1b
10b 0,563 0,359 0,007632 2,410 1,343 0,004531
N=10, h=0,1 0,468 0,291 0,006047 4,919 2,859 0,025725
0,1b
10b 1,926 1,321 0,058828 0,175 0,107 0,001705
N=10, h=0,05 0,468 0,291 0,006048 4,919 2,859 0,025726
0,1b
10b 1,922 1,318 0,058749 0,175 0,107 0,001704
N=15, h=0,1 0,06618 0,05588 0,005971 0,643 0,410 0,010300
0,1b
10b 2,262 1,388 0,023747 0,133 0,076 0,002004
N=15, h=0,05 0,0662 0,0559 0,005973 0,634 0,410 0,0104
0,1b
10b 2,258 1,385 0,0237 0,133 0,076 0,002004
120
Для определения качества аппроксимации оценок коэффициентов использован показатель ошибки аппроксимации:
к-к
Ak =
•100%
где к — заданный коэффициент и к — его оценка. Для вычисления относительной ошибки приближения рядов данных применен показатель средней ошибки аппроксимации:
х-х
1 N
X
100%
где X _ матрица-столбец реального ряда данных, — матрица-столбец расчетных данных, N — длина заданного ряда данных.
В табл. 1 — 3 приведены оценки погрешностей предложенного алгоритма параметрической идентификации на примере уравнений (1) — (3) при вариации к, N и стартовых догадок.
Таблица 2
0,1r 10r
Aa, % Ab, % Acp, % Aa, % Ab, % Acp, %
N=5, h=0,1 0,061 0,492 0,015584 4,889 53,166 0,62250
0,1s
10s 0,417 4,617 0,052580 4,237 43,659 0,54581
N=5, h=0,05 0,061 0,492 0,015609 4,889 53,166 0,62249
0,1s
10s 0,417 4,617 0,052580 4,237 43,662 0,54585
N=10, h=0,1 0,139 0,620 0,065336 7,406 29,423 3,790910
0,1s
10s 86,439 550,301 38,19000 0,377 1,603 0,177503
N=10, h=0,05 0,139 0,620 0,065336 7,406 29,422 3,791464
0,1s
10s 86,439 550,300 38,19017 0,377 1,603 0,177503
N=20, h=0,1 0,026 0,0109 0,029774 4,998 2,575 5,400658
0,1s
10s 0,160 0,260 0,106737 0,0410 0,0764 0,027951
N=20, h=0,05 0,026 0,0109 0,029774 4,998 2,575 5,400658
0,1s
10s 0,160 0,260 0,106737 0,0410 0,0764 0,027989
Оценки погрешности алгоритма для уравнений (1) — (3), представленные в табл. 1 — 3, демонстрируют довольно адекватную идентификацию параметров при разных начальных приближениях, независимо от шага и количества точек. Наибольшие отклонения найденных оценок параметров от заданных наблюдаются при малом количестве точек и начальных приближениях 10я и 0.1b; 10s и 0.1г, 10ю и 0.1fc. Несложно заметить, что начальные приближения свободного члена уравнения (1)
121
слабее сказываются на невязке, чем коэффициент при фазовой переменной (очевидно, чем больше порядок переменной, тем существеннее вклад начальных приближений ее в невязку). Кроме того, важно подобрать правильное соотношение начальных приближений, включая знаки. В практическом аспекте это задача лежит в плоскости предварительного анализа моделируемой ситуации и аналитического исследования модели. В целом, с помощью предложенного подхода искомые оценки коэффициентов найдены с довольно хорошей точностью (десятые доли процента отклонения от заданной величины), что обеспечило высокое качество аппроксимации решений (порядок точности в среднем десятые доли процента).
Таблица 3
0,1k 10k
Aa, % Ab, % Acp, % Aa, % Ab, % Acp, %
N=10, h=0,1 0,454 0,0015 0,0237474 1,441 0,0070 0,0348854
0,1 ю
10® 4,651 0,0116 0,2204293 3,076 0,0104 0,0768792
N=10, h=0,05 0,454 0,0015 0,0237325 1,441 0,0070 0,0348854
0,1 ю
10® 0,855 0,0315 0,2261174 3,067 0,0109 0,0748630
N=20, h=0,1 0,662 0,0040 0,0437047 0,481 0,0122 0,1075936
0,1®
10® 1,673 0,0690 0,6519197 0,126 0,0199 0,1784524
N=20, h=0,05 0,369 0,0177 0,1567043 0,481 0,0122 0,1075936
0,1 ю
10® 0,455 0,0131 0,1268999 1,945 0,0135 0,1344047
N=30, h=0,l 0,647 0,00301 0,2722964 0,264 0,00201 0,2436575
0,1 ю
10® 2,699 0,01290 1,1628419 0,214 0,00101 2,6990000
N=30, h=0,05 0,647 0,00300 0,2722965 0,100 0,00414 0,4215779
0,1®
10® 2,470 0,00250 0,6397427 1,353 0,00246 0,2031063
Описанный в работе алгоритм с незначительными модификациями ранее был успешно применен для оценки параметров моделей различных социально-экономических процессов, построенных с использованием систем нелинейных дифференциальных уравнений с постоянными коэффициентами [5, 6, 13]. Вместе с тем он может оказаться полезным в прикладных научных исследованиях других областей науки, а также в образовательном процессе студентов высших учебных заведений.
Список литературы
1. Амосов А. А, Дубинский Ю. Л., Копченова Н. В. Вычислительные методы для инженеров. М.: Высшая школа, 1994. 544 с.
2. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М.: Бином. Лаборатория знаний, 2003. 632 с.
122
3. Безручко Б. П., Смирнов Д. А. Математическое моделирование и хаотические временные ряды. Саратов: ГосУНЦ «Колледж», 2005. 320 с.
4. Калиткин Н. Н. Численные методы. М.: Наука, 1978. 512 с.
5. Курилова Е. В., Кулаков М. П., Хавинсон М. Ю., Фрисман Е. Я. Моделирование динамики добычи минеральных ресурсов в регионе: эконофизический подход / / Информатика и системы управления. 2012 г. № 4(34). С. 3 — 13.
6. Курилова Е. В., Кулаков М. П., Хавинсон М. Ю., Фрисман Е. Я.. Построение и исследование обобщенной модели динамики макрокомпонент регионального развития (на примере Еврейской автономной области) / / Региональные проблемы. 2009. № 1(14). С. 5— 10.
7. Охорзин В. А. Компьютерное моделирование в системе MathCAD. М.: Финансы и статистика, 2006. 144 с.
8. Охорзин В. А. Прикладная математика в системе MathCAD. СПб.: Изд-во «Лань», 2008. 352 с.
9. Самарский А. А. Введение в численные методы. М.: Наука, 1987. 269 с.
10. Фрисман Е. Я., Хавинсон М. Ю., Аносова С. В., Фишман Б. Е., Петров Г. И. Системная динамика регионального развития: подходы к моделированию блока экономики (на примере Еврейской автономной области) / / Пространственная экономика. № 3(11). 2007. С.134—146.
11. Gill P. E., Murray W., Wright M. H. Practical optimization. London: Academic Press, 1981. 496 p.
12. Kahaner D., Moler C., Nash S. Numerical Methods and Software. Prentice-Hall International, Inc., 1989. 571 p.
13. Khavinson M. Yu., Kulakov M. P., Mishchuk S. N. Prediction of Foreign Labor Migration Dynamics at the Regional Level / / Studies on Russian Economic Development. 2013. No.2(24). P. 170—178.
•Jc -Jc -Jc
123
