CC BY
48
№ 2(26) 2010
В. В. Андреев, М. И. Семёнов
Программное приложение для решения задач оптимальной параметрической идентификации динамических моделей: применение для прогнозирования динамики социально-экономической системы США
В рамках данной работы разработано программное приложение для решения задач оптимальной параметрической идентификации моделей. Использован алгоритм Левенберга-Марквард-та в модификации Флетчера, дающий лучшие результаты для решения переопределенных нелинейных систем уравнений по сравнению со стандартными методами. Задачи подобного типа возникают при настройке модели на реальныеданные.
Библиотека Optimization Toolbox пакета MATLAB представляет собой очень мощный аппарат для решения широкого круга задач оптимизации. Кроме того, основной пакет MATLAB предоставляет средства для целевой оптимизации, например, оператор «обратный слэш» для решения системы линейных уравнений или функция fminsearch для нахождения локального минимума нелинейных функций. Однако в случае, когда задача описана нелинейной системой уравнений, и требуется минимизация методом наименьших квадратов, применение функции fminsearch является неэффективным. Функция Isqnonlin, входящая в библиотеку Optimization Toolbox, тоже является недостаточно гибкой.
Разработанное программное приложение применено для исследования математической модели, описывающей динамику социально-экономической системы на примере США. Ранее аналогичные модели были использованы для анализа социально-экономической системы России и показали свою достаточную адекватность.
Нелинейные задачи идентификации моделей
Постановка задачи
Пусть задан вектор-столбец у экспериментальных данных, элементы которого у( зависят от переменной х., т. е.
У (Х1У
У (хп \
Цель состоит в том, чтобы получить регрессию у(х) на некоторую функцию f(x,с), где с — вектор неизвестных параметров идентифицируемой модели. Введем вектор-столбец
f (х,, с)
,f (, с),
Разность
f - у
(1)
46
№ 2(26) 2010
представляет собой вектор-столбец так называемой невязки
'г (с)л
КГ (. СЬ
который должен быть нулевым при условии, что модель f(x,с) является физически совершенной, а экспериментальные данные
у. = у(х), /' = 1.....л, получены без какой-либо
погрешности измерений. Очевидно, в реальной ситуации этого не наблюдается, и задача заключается в минимизации вектора невязки г. Для решения обычно применяется метод наименьших квадратов [1]. Он применяется для поиска такого вектора параметров модели с, для которого сумма квадратов невязок
S = г7 г
становится минимальной. Необходимые условия минимума при этом записываются так:
dS( с)/ дс/ == 2 J 7 г = 2ь = 0, j = 1,..., m. (2)
Здесь
'Э г (xvc)/dcj
дг (хп ,c)/dCj ,
Запишем систему уравнений (2) в векторном виде
'aS(c)/Эс1 dS( с)/дст
= 2J7 г = 2ъ = 0.
Здесь
гдг (х1, с)/Э с1 ••• Э г (х1, с)/дс„
д г (хп ,с)/дс1 ••• Эг (х„ ,с)/Эс, матрица Якоби; u = J7г.
m у
Очевидно, что и, г — функции неизвестных параметров с. В общем случае вектор с не может быть получен в конечном виде, и его необходимо рассчитать в процессе итерации. Пусть с на к + 1 шаге итераций имеет простейшую форму
~(к+1) _ ~(к)
= с1К > + Лс1К >.
(4)
Пусть вектор невязки г непрерывен по параметрам с. Тогда представим вектор невязки на к +1 шаге итераций в виде ряда Тейлора:
г(к+ 1) _ -(к)
= Гп' +
т ^г(к)
Ас,1к) +... = Э с w у
/
= г(Ч + Jw Ас(Ч +....
(5)
Здесь — якобиан, соответствующий к шагу итераций.
Умножим уравнение (5) на матрицу ) слева:
(и7)(Чг№+1> = (и7)(Ч+ (и7)(ЧАси +.... (6) Введем обозначение
= (J7 f J
Тогда формула (6) в линейном приближении примет вид:
А1к )Лс(к) > f(k+1)
(7)
(3)
Уравнение (7) является основой для целого ряда методов оптимальной идентификации параметров нелинейных моделей [1].
Метод Ньютона-Рафсона
В методе предполагается, что норма вектора невязки падает так быстро, что вторым членом в формуле (7) можно пренебречь. Следовательно, получим выражение
д (к) Дс№) =
(к)
(8)
которое может быть использовано для получения величин Ас14. Затем по формуле (4)
47
№ 2(26) 2010
находим значения с№+1) для следующего шага итераций. Имеется несколько модификаций метода, например, добавляют в левую часть уравнения (8) член аДс(Ч, где а < 1.
Метод Левенберга-Марквардта
В этом методе [2, 3] делается искусственное предположение, что второй член в формуле (7) может быть аппроксимирован с помощью
МДс^(9)
Здесь — масштабирующий параметр; й — подходящая диагональная матрица весов. Ее часто выбирают равной единичной матрице I или й = сйад (А(Ч ). Следовательно, с учетом выражения (9) уравнение (7) преобразуется к виду
А1к)Ас1к) + Х1к]ОАс1к) = -ь1к) . (10)
Идея принадлежит Левенбергу [2]. Стратегию изменений улучшил Марквардт [3]. Он ввел вместо единичной матрицы I диагональную матрицу весов й = сйад (А(Ч ). Кроме того, Марквардт внес новый фактор Я, отражающий согласование прогнозируемой суммы квадратов с реальной на текущем шаге итерации. Если значение Я лежит в заданных пределах (Я(, Яг), то параметры итерации не изменятся. В противном случае последуют изменения параметра X и вектора и. Значение X уменьшается вдвое, если Я > Я(. Затем, если значение X становится ниже критического уровня Хс, данный параметр будет обнулен. В результате следующий шаг процедуры итерации будет совпадать с алгоритмом Ньютона-Рафсона (8). Если Я < Я(, то вектор и устанавливается так, чтобы его компоненты находились на интервале 2 < и( < 0. Если при этом X = О, то последуют изменения параметров Хс и X. Затем Флетчер значительно улучшил стратегию Марквардта по адаптации X и сделал ее законченной [4].
Уравнение (10) может быть преобразовано к виду
(А(Ч +^)0)Дс"( ' =-г)(Л). (11)
Чем выше масштабирующий параметр X, тем ближе результат к устойчивому решению в направлении наискорейшего спуска. При X = 0 метод совпадает с алгоритмом Ньютона-Рафсона (8), который является менее стабильным, и результат может расходиться.
В рассмотренном алгоритме стратегия поиска основана на сопоставлении прогнозов решения для следующей итерации. Если прогноз близок к реальности, то величина X может быть снижена. Если прогноз плохой, то параметр X должен стать выше для того, чтобы стабилизировать процесс нахождения решения.
Метод Левенберга-Марквардта в модификации Флетчера для решения нелинейной задачи идентификации моделей методом наименьших квадратов
Принцип действия алгоритма [1, 5, 6] показан на блок-схеме (рис. 1).
Пусть имеется общая переопределенная система нелинейных алгебраических уравнений (1). Ее решение, оптимальное в смысле наименьших квадратов, ищется путем минимизации эвклидовой нормы или суммы квадратов невязок. Необходимые условия минимума имеют вид уравнений (2) или (3). Вектор и = и7г должен стремиться к нулевому в точке с*, соответствующей оптимальной идентификации параметров модели. Параметры с после к-о\л итерации ищутся по формуле (4). Процедура вычислений повторяется до тех пор, пока не будет выполнено неравенство
|5(* +1) _ 5У0| <е> (12)
Здесь — сумма квадратов невязок после к-ой итерации.
№ 2(26) 2010
Рис. 1. Алгоритм модифицированного метода Марквардта
Описание системы для решения задач оптимальной параметрической идентификации моделей
Программное приложение для решения задач оптимальной параметрической идентификации динамических моделей разработано с использованием средств пакета МАТ!_АВ. Порядок работы системы показан
на рис. 2, где СНУ — система нелинейных уравнений, СКО — среднеквадратичное отклонение.
Стартовое окно системы представлено на рис. 3. С его помощью осуществляется переход между тремя основными блоками системы. Рассмотрим их назначение.
1. В данном блоке происходит подготовка экспериментальных или статистиче-
49
№ 2(26) 2010
Загрузка экспериментальных или статистических данных
Подготовка данных
Выделение трендов с помощью вейвлет-преобразования Масштабирование данных по значению Изменение количества отсчетов
§
I
§
I *
I
I
ч
I
I £
§ I
И
I
§
48
0
1
I
1=
0 »
1
1
§
&
I
I
си
I
I
о §
Рис. 2. Порядок работы системы параметрической идентификации моделей
ских данных для последующих расчетов (рис. 4).
Желательно, чтобы все данные были сглажены для выявления основных тенденций в их временном поведении, а также приведены к их реальным физическим значениям размерностей. Это осуществляется с помощью вейвлет-преобразования. Ранее авторами в работе [7] было описано приложение, раз-
работанное с применением пакета МАТ1_АВ и предназначенное для обработки статистических данных, включая и выявление трендов. В итоге получаем ожидаемые входные переменные математической модели.
2. Здесь осуществляется задание математической модели в виде символьной системы уравнений, например, системы дифференциальных уравнений (рис. 5).
50
№ 2(26) 2010
Рис. 3. Стартовое окно системы оптимального определения параметров модели по экспериментальным или статистическим данным
Общий вид модели запишем так:
и параметров модели, которые также задаются в данном блоке приложения. Затем осуществляется численный расчет левой части системы (в случае модели, заданной системой дифференциальных уравнений — численное дифференцирование). В результате получим вектор-столбец у. Подбор начального приближения параметров модели в некотором заданном диапазоне осуществляется циклически, исходя из требования минимизации суммы квадратов невязок в численно рассчитанной левой части и ее же, полученной по расчету правой части системы с имеющимися исходными данными выходных переменных и набору параметров, т.е. вектор невязки г равен:
Р = О (...).
Р - у = г.
Здесь многоточие условно обозначает В итоге получаем значения начального всю совокупность переменных состояния приближения параметров модели.
Рис. 4. Окно подготовки экспериментальныхданных
51
№ 2(26) 2010
Рис. 5. Окно задания модели и ее параметров
3. В этом блоке происходит подбор параметров системы (рис. 6). По начальному набору значений параметров модели по алгоритму Левенберга-Марквардта в модификации Флетчера производится минимизация суммы квадратов невязок в исходных данных и данных, рассчитанных по модели. В итоге находим значения параметров модели, при которых она наиболее точно описывает статистические или экспериментальные данные.
Применение приложения для параметрической идентификации модели социально-экономической динамики США
Важное место при принятии управленческих решений занимает умение предвидеть возможные пути развития различных процессов, протекающих в социально-эконо-
мической системе. В предыдущих работах [8-13] на основе математического моделирования была проанализирована ситуация в России в 90-е годы 20 века и в начале 21 века. В этих работах по результатам проведенных исследований был сделан вывод о том, что очередной кризис социально-экономической системы страны может наступить в 2009/2010 годах. Сейчас, когда контуры мирового финансового кризиса довольно четко обозначились, можно утверждать, что разработанные математические модели достаточно корректны. Кроме того, предложенные модели были использованы также для исследования прошедших уже этапов эволюции России (СССР). В частности, они правильно позволили смоделировать распад СССР в 1991 году [9] и дефолт в России 1998 года [8, 10,12,13].
В данной работе на основе математической модели, элементы которой взаимодей-
№ 2(26) 2010
о
3
р и
I
Рис. 6. Окно подбора параметров модели
ствуют в соответствии со схемой «хищник-жертва», представленной на рис. 7, исследована динамика социально-экономической системы США. Здесь стрелка исходит от хищника и заканчивается на жертве. В качестве основных факторов, определяющих динамику социально-экономической системы, здесь выбраны элементы: Х1, обозначающий доходы бюджета; Х2 — валовой внутренний продукт (ВВП); Х3 — расходы на финансирование науки и образования; Х4 — личные потребительские расходы; Х5- дефицит бюджета.
В соответствии с взаимодействиями между «хищниками» и «жертвами» на рис. 7 математическую модель запишем так:
бх^ _ бх
= "Р^Л + Рг*з ~ Рз~ РдХ2<
Рис. 7. Взаимодействия элементов социально-экономической системы США
аТ
= У1 -Угхгхз "Уз*з>
Л
^ = 51х1х4 - 2.
* (1 + 54
= ф1 х^5 +ф2х2х5 -Фзх5 + ф,
(13)
53
№ 2(26) 2010
На рис. 8-12 представлены данные статистики по индексам экономики США с 2000 по 2008 год включительно [14-16], используемые для идентификации параметров математической модели (13). Применим для этой модели разработанное программное приложение. Процесс подготовки данных статистики показан на рис. 4. Сначала из файла загружаются данные статистики, и указывается временной диапазон — с 2000 по 2008 год. Затем они сглаживаются с помощью вейвлета Добеши-8 с опреде-> ■ ленным уровнем приближения [7]. Заданием [а новых значений максимумов и минимумов g кривых добиваемся их соответствия реаль-5 ным масштабам на рис. 8-12. Выбирается SS необходимое количество отсчетов для опи-S сания кривых с учетом быстродействия сис-| темы при расчетах. Полученные результаты
4 сохраняются для дальнейшей работы с эти-I ми данными.
| На втором этапе задаются модель и ее § параметры. При этом правая часть системы
5 дифференциальных уравнений (13) в сим-Л вольном виде заносится в таблицу (рис. 5). | Далее производится численный расчет ле-| вой части системы уравнений (13) по имею-
6 щимся данным статистики. Задав началь-| ный диапазон изменения параметров от g -100 до 100, а также количество циклов
подбора равным 50, рассчитаются началь-г§ ные значения параметров модели. При этом g для решения системы обыкновенных диф-§ ференциальных уравнений (13) использо-г? ваны решатели ode, имеющиеся в пакете 8 MATLAB. Подбор осуществляется путем | начального приближения расчетной и по-| лученной после подбора коэффициентов левых частей модели (13), как показано на Е| рис. 5. В результате получаем начальные | значения параметров модели для дальней-^ шего подбора.
§ На последнем этапе подбора (рис. 6), & указав количество итераций, визуально на-§ блюдаем на экране монитора процесс под-| бора коэффициентов модели в соответст-& вии с алгоритмом Левенберга-Марквардта & в модификации Флетчера.
В результате были получены следующие значения параметров модели (13): а1 =0,0484, а2 = 3,0127-10-4, а3 =0,0555, а4 = 2,476-10-®, а5 =-0,1368, р1 = 0,001, р2 = 0,0061, р3 = 0,0013, р4 = — 0,771 4, у, = -3,2645-10-®, у2= 8,1158-10-4, у3 = -0,0738, 51 = 4,102-10-6 , 62 = 0,0 6 69, 64 =-1,5 4 09, Ф1 = - 1,521-Ю"4, ф2 = -0,0889, ср3 = -0,3905, Ф4= 16,2482. Начальные условия для системы дифференциальных уравнений (13) были заданы близкими к соответствующим статистическим данным на конец 2000 года.
По полученным значениям параметров модели построено решение системы дифференциальных уравнений (13) на интересующем интервале времени с 2000 по конец 2010 года. Результаты прогнозирования различных показателей экономики США на конец 2010 года, полученные по статистическим данным с 2000 по 2008 год включительно, представлены на рис. 13-17 (сплошные кривые). На указанных рисунках пунктирные кривые дублируют статистические данные, приведенные на соответствующих рис. 8-12. «Хвосты» пунктирных кривых после 2008 года на рис. 13-17 соответствуют статистическим данным 2009 года. Видно, что решения модели (13) с 2000 по 2008 год достаточно хорошо описывают тенденции динамики соответствующих статистических данных за тот же период времени.
На рис. 13 видно, что на конец 2009 года доходы бюджета США остались практически на уровне конца 2008 года. В то же время, если бы сохранились тенденции динамики эволюции социально-экономической системы, имевшие место с 2000 по 2008 год, доходы бюджета на конец 2009 года были бы примерно в 1,5 раза выше (сплошная кривая). Из рисунка следует, что при таких тенденциях прогнозируемая на конец 2010 года величина доходов бюджета многократно возрастает. Это означает, что по прогнозу к концу 2010 года должен наступить так называемый «режим с обострением», т.е. кризис социально-экономической системы страны. Если не принять своевременные оптимальные управленческие решения, этот
№ 2(26) 2010
16001-'-'-'-
2000 2002 2004 2006 2008
2000 2002 2004 2006 2008
Рис. 8. Доходы бюджета в млрд долл.
Рис. 9. Валовой внутренний продукт (ВВП) в процентах к предыдущему году
350
300
250
200
150
100
2000 2002 2004 2006 2008
2000 2002 2004 2006 2008
Рис. 10. Расходы на науку и образование в млрд долл.
Рис. 11. Личные потребительские расходы в процентах к предыдущему году
-200
2000 2002 2004 2006 2008
Рис.12. Дефицит бюджета в млрд долл.
кризис будет усугубляться и приведет к тя- личных потребительских расходов граждан
желым последствиям. США за истекший 2009 год не изменилась,
Из рис. 14 и 16 следует, что общая тен- и прогнозируемые сплошные кривые доста-
денция эволюции динамики ВВП, а также точно хорошо согласуются с ними.
55
№ 2(26) 2010
10000 ......................, 8
9000 Данные статистики |
' — Прогнозируемыеданные 1 6
8000 2009 год 1
4
о
§ 7000 ч о 2
6000 0
-.о
X" 5000 X™ -2
4000 -4
3000
2000 -6
2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 Время, год
Рис. 13. Доходы бюджета
§
I §
I *
I
I
ч
I
I £
§ I
И
I §
48
0
1
I 1=
0 »
1
1 §
&
I [
<ъ
I &
о §
2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 Время, год
Рис. 14. Валовой внутренний продукт
• : А ч
• 1 " 1, 1»1 , • ,
1 » « ч •' «; % / « ■ I • » •» • ■ 1. »
■ 1 1« »• \! . ; ч и 1 .
\ - ч | *
■ Д энные статистики \ ^
* 2009 год ^
2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 Время, год
Рис. 15. Расходы на науку и образование
2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 Время, год
Рис. 16. Личные потребительские расходы граждан
1000
Данные статистики — Прогнозируемыеданные 2009 год
2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 Время, год
Рис. 17. Дефицит бюджета
Интересно сравнение статистических данных и прогнозной кривой расходов на науку и образование на конец 2010 года (рис. 15). Прогноз предсказывал (сплошная кривая) снижение указанных расходов практически
до уровня 2001 года. В реальности произошло существенное увеличение указанных расходов (пунктирная кривая). Это свидетельствует о понимании руководством США важности развития новых высоких технологий и подго-
56
№ 2(26) 2010
товки высококвалифицированных специалистов для вывода страны из кризиса.
В то же время в 2009 году произошло резкое увеличение дефицита бюджета (рис. 17). Тенденция, сложившаяся с 2000 по конец 2008 года (сплошная кривая), здесь резко расходится с реальными статистическими данными 2009 года.
Таким образом, проанализировав ситуацию, нельзя сказать, что руководство США предприняло все необходимые меры для вывода страны из кризиса. Ситуация в целом остается достаточно трудной. Без принятия в 2010 году соответствующих решений в 2011 году социально- экономическая система США может оказаться в еще более сложной ситуации.
Заключение
Разработанная система идентификации нелинейных систем способна достаточно быстро и удобно осуществить подбор параметров заданной математической модели, при которых она наиболее точно описывает данные статистики или экспериментальные данные. Это ускоряет процесс параметрической идентификации сложных нелинейных систем.
Описок литературы
1. Гдал Ф., Мюррей У., Райт М. Практическая оптимизация. М.: Мир. 1985. 509 с.
2. Levenberg К. A Method for the Solution of Certain Problems in Least Squares// Quart. Appl. Math. 1944. Vol. 2. Pp 164-168.
3. Marquardt D. An Algorithm for Least-Squares Estimation of Nonlinear Parameters// SIAM J. Appl. Math. 1963. Vol. 11.Pp 431-441.
4. Fletcher Ft. A Modified Marquardt Subroutine for Nonlinear Least Squares// Rpt. AERE-R 6799, Harwell. 1971.
5. Balda M. LMFsolve: Levenberg-Marquardt-Fletch-er's algorithm for nonlinear least squares problem// MathWorks, MATLAB Central, File Exchange, ld = 1606. 2007.
6. More J. J. The Levenberg-Marquardt Algorithm: Implementation and Theory. In book: Numerical
Analysis. Lecture Notes in Mathematics 630/ed. .Si
G. A. Watson: Springer Verlag. 1977. Pp 105- |
116. ^
7. Андреев В. В., Семёнов М. И. Создание и ис- Sä
пользование системы обработки и анализа дан- <§
<ц
ных с применением пакета Matlab// Прикладная информатика. 2008. №2 (14). С. 85-92.
8. Андреев В. В., Карпова О. В. Математическое ^ моделирование социально- экономических процессов в России конца XX и начала XXI веков// Нелинейный мир. 2007. Т. 5. № 12. С. 773777.
9. Карпова О. В., Андреев В. В. Моделирование динамики одной социально- экономической системы на основе модели типа «хищник- жертва» // Математика. Компьютер. Образование: Сб. научных трудов. Том 1/ Под ред. Г. Ю. Ризниченко. — М.-Ижевск: НИЦ «Регулярная и хаотическая динамика». 2007. С. 194-202.
10. Карпова О. В., Андреев В. В. Исследование социально-экономической динамики России на основе модели типа «хищник — жертва» // Математика. Компьютер. Образование: Сб. трудов XV международной конференции. Том 1/ Под общей редакцией Г. Ю. Ризниченко. — Ижевск: Научно-издательский центр «Регулярная и хаотическая динамика». 2008. С. 231-236.
11. Андреев В. В., Карпова О. В. Попытка построения математической модели социально- экономической системы: исследование на примере Чувашской Республики// Вестник Чувашского университета. Гуманитарные науки. 2008. №1. С. 385-390.
12. Андреев В. В., Васильева Е. А. Математическое моделирование и исследование динамики социально- экономической системы России// Известия РАЕН. Дифференциальные уравнения. 2009. №14. С. 25-38.
13 .Андреев В. В., Ярмулина О. О. Математическое моделирование динамики социально-экономической системы (на примере России)// Нелинейный мир. 2009. Т. 7. №6. С. 464-474.
14. http://www.cbo.gov — U. S. Congressional Budget Office.
15. http://www.bea.gov — U. S. Department of Commerce, Bureau of EconomicAnalysis.
16. http://www.gpo.gov/fdsys — U. S. Government Printing Office, Federal Digital System.
