УДК 004.942:338.27:504.062
ПРЕИМУЩЕСТВА ДИФФЕРЕНЦИАЛЬНЫХ МОДЕЛЕЙ В ЭКОЛОГО-ЭКОНОМИЧЕСКОМ МОДЕЛИРОВАНИИ
А.В. Затонский
Березниковский филиал Пермского национального исследовательского политехнического университета
E-mail: [email protected]
Показано, что использование эколого-экономических моделей в форме обыкновенных дифференциальных уравнений в задачах прогнозирования развития управляемых систем, предпочтительнее, чем полиномиальных моделей и временных трендов, в смысле адекватности и уменьшения возможного разброса прогнозов, необходимых для принятия решений.
Ключевые слова:
Принятие решений, эколого-экономическая модель, обыкновенные дифференциальные уравнения.
Key words:
Decision support, ecological-economic model, ordinary differential equation.
В практике экологического, социального и экономического моделирования часто используются модели динамики вида
Я *0), г(>), 0 = а0 + £ ах 0) + £ ^ (0,
, к
где 6 (/)={х1(/),х2(/),_| - вектор факторов,
6- вектор возмущений, у(-) - реакция исследуемого объекта;
или ,у( х(1), г(1),1) = а0 + П аЛ (1) + П7;(!),
7
либо, для функций одного аргумента -
Я х(1), г(1), Г) = а + £ а,.х(1)' + £&,.г(1)7,
, к
I
либо модели временны х рядов в форме 7(1) = £ а/.
/=0
Поиск в Интернет позволяет обнаружить авторефераты, в которых используются подобные модели: [1. С. 29], [2. С. 27] [3. С. 10 и 16] и т. п. Назначение моделей, обычно - исследование истории и, на ее основании, прогнозирование свойств объекта в зависимости от принятых решений х(^.
Такие модели можно упрощенно понять, например, так: если вкладывать в предприятия (отрасль) инвестиции по графику x1(t), то на выходе получим чистый дисконтированный доход (или другой показатель экономической эффективности) У^О, ^(0) с учетом спроса на продукцию (возмущающего воздействия) ^(0. Дальше обычно речь идет об идентификации ^ и Ьр об учете обратных связей, выраженных некоторой функцией F(y), а точнее
у (х(?X г О1)) = ао + £ а,.*,. О1) +
+£7,(0 - ЯхХ(0, <0)) (1)
к
и т. п. При этом молчаливо принимается предположение, что существует прямая связь между факторами и значением реакции, а единственный динамический элемент в модели - чистое запаздывание (например, в моделях вида _у (х(1)) = а0 + £ ах, (1 -А1)).
,
Однако подобное предположение не всегда близко к реальности. Например, удобряя поле по определенным правилам, можно получить рост урожая (и дальнейшие экономические или социальные бонусы). То есть достоверно, из множества наблюдений, известно, что внесенное количество удобрений x1 ускоряет рост урожая в каком-то диапазоне вноса удобрений:
ду(х,г,1) ^ а + а[х(1), V!: 0 < х(1) < хтах, а > 0, д!
а снижение количества осадков в определенных условиях снижает скорость роста:
ду(х^ « а + Ь г(1), V! : гтт < г(1) < г_
Ь1 > 0, V!: г(1) < 0.
Для сложных систем, особенно учитывающих естественные процессы в природе, массовая идентификация коэффициентов связи между у(Р) и xi(t) без убедительного доказательства их взаимной независимости приводит к порождению «попугайских моделей» [4], адекватно интерполирующих прошлое, но не способных к прогнозу будущего -что, собственно, и требуется при построении моделей поддержки принятия решений.
К ним же относятся попытки экстраполировать у(^ вперед по данным временных рядов (трендов), особенно с учетом ошибок или ненаблюдаемых внешних возмущений, что проиллюстрировано ниже.
Возникает вопрос - что же идентифицировать при построении динамической экономико-математической модели: связь между фактором и реакцией или связь между фактором и динамикой изменения реакции под воздействием этого фактора? В теории автоматического управления, как известно, фактор (или динамика его изменения) линеаризуется, а затем исследуется его влияние на динамику поведения объекта. Подобные подходы к экономико-математическим системам также разработаны очень давно. Например, в [5. С. 99] формулируется модель экологического равновесия
ду(х У, О = у - у2 - х(?^у^ 0 < х(?^ ^ (2)
дt
соответствующей по форме (1), от которой недалеко как до доходности (определяемой здесь квотой вылова x(í)), так и до катастроф в развитии популяции, что, собственно, и рассматривается далее в книге.
Таким образом, на уровне общенаучных рассуждений получается вывод, что в эколого-эконо-мических моделях лучше использовать в качестве основы дифференциальные, а не алгебраические уравнения динамики системы.
Попробуем проверить вывод на нескольких примерах. Добавим в ур. (2) возмущающее воздействие - например, сезонное (периодическое) влияние погодных условий на воспроизведение популяции вида z(t)=sin(b2t) и получим
ду( х у, о=у- у2 - х(1) у- ^1 ).
д1
Построим в МаЛАВ простую модель, положив рост квоты вылова в виде линейной зависимости x(t)=ao+a1t (единица измерения времени - год). В [5] аналитически доказано, что при использовании варианта модели
система теряет устойчивость при х>0,25 (в отсутствие обратной связи, связывающей вылов с текущим значением популяции). Проверим это, включив Manual Switch в нижнее положение, задав Gain=0 (рис. 1), чтобы исключить пока возмущающие воздействия, ai=0 и задавая последовательно ao=0,22, ao=0,25 и ao=0,28. При этом начальное
Рис. 2. Изменение по годам (ось абсцисс) динамики численности популяции (ось ординат) в зависимости отдоли вылова X (ряды данных)
значение (Initial Condition) в блоке интегратора « Y(t)» установим в 1 (полная популяция).
Блок MinMax в модели предназначен, чтобы не допускать падения популяции ниже 0 с последующим аварийным остановом расчетов.
Действительно, при превышении x>0,25 наблюдается катастрофическое снижение популяции в конце десятилетнего периода (рис. 2).
Так как теоретический результат совпадает с расчетным, будем считать, что построенная модель адекватно отражает (2). Включим Manual Switch в верхнее положение (задавая, таким образом, долю от размера популяции) и установим Gain в произвольным образом выбранное значение 0,05; зададим a0=0,2 и a1=0,025, период изменения возмущения z(t) установим равный 2^ (один год). Получим зависимость популяции от времени, показанную на рис. 3 сплошной линией. Попробуем спрогнозировать развитие популяции поданным 1-6 годов, проведя регрессионный анализ. Используем для этого все доступные в MS Excel модели. Полу-
чим уравнения динамики вида y(t) = ^ ait1: i=0
1. y(f)=-0,045481f+0,927281; (3)
2. y(i)=0,008226t2-0,093550t + 0,969766; (4)
3. y(i)=-0,002270f+0,028255f--0,138835t+ 0,987380; (5)
4. y(i)=0,000724t4-0,010835i3+0,060137t2- -0,177522t+0,995103; (6)
5. Xi)=0,930131e-°'055976' (7)
Из рис. 3 очевидно, что экстраполяция в данном случае получается неудачной. Здесь и во всех следующих графиках по оси я отложено время в годах, по оси у - доля популяции от начального значения у(0)=1.
Будем искать решение задачи в виде
т
у(1) = О) +£ а,х(1), полагая x(t)=0,2+0,025t из-
,=1
вестным (так как решение о доле вылова, принятое или планируемое лицом, принимающим решения, идолжно быть известным). Для этого подготовим вспомогательную таблицу и произведем поиск решения наименьшего квадратичного отклонения
/ / \\2 I { т \ \
а:
a,X(t)'
, где J*(t):te{t,j - полу-
ченные при помощи модели (рис. 1) «экспериментальные» значения, при ^=1,4на интервале времени с1по6 года. «Спрогнозируем» развитие ситуации при принятом решении x(t)=0,2+0,025t на период с 7 по 10 год. Получим, вне зависимости от степени полинома, неудовлетворительные по качеству прогнозы (рис. 4).
Модели, по которым произведен расчет на рис. 3 и 4, получились вполне «попугайские» (по терминологии К.С. Лосева): они неплохо интерполируют исходные данные (на интервале до 6 года включительно), но качество экстраполяции оставляет желать лучшего.
Перейдем к моделированию воздействия факторов на динамику объекта. Подберем методом наименьших квадратов коэффициенты уравнений
—Y(t)
...1
— 2 ---3 --4 — 5
I, лет
Рис. 3. Неудачные попытки экстраполяции численности популяции зависимостями У(1): исходные данные: 1) ур. (3); 2) ур. (4); 3) ур. (5); 4) ур. (6); 5) ур. (7)
Рис. 4. Неудачные попытки экстраполяции численности популяции зависимостями: 1) у()=1,291121-1,81922х^); 2)у()=2,245758-9,01517х() + 13,17732ха)2; 3) у()=5,046407-40,798х(1)+131,2891х(1)2-143,871х(г)3;
4) уа)=5,036115-40,6991х()+131,0587х()2-144,083х а)3 +0,87505х()
dy(t) dt
= an +
£a,x(t У + £bj.y(t)j =
j=i
= ao + £ a, (0,2 + 0,025t)' + £ bj.v(t) j
i=l j= 1
для m=1,4, n=1,4 и начального условия y (0)=1. Развернем уравнения:
Ф<0
dt
■ = a0 + a1 (0,2 + 0,025t) + b1 _y(t)
)
dt
= a0 + a (0,2 + 0,0251) +
Ф<0
dt
+a2 (0,2 + 0,0251)2 + b _y(t)
= a0 + a1 (0,2 + 0,0251) + b1 _y(t) + b2 j(t)2
)
dt
= a0 + a1 (0,2 + 0,0251) +
j(t + At) « y(t) + At
a0 +£ a,x(t У + £ j(t У
i=1 j=1
t < 10.
3. Вычисление суммы квадратов отклонений ¿=£(у(^)-у*и))2, tФ6.
4. Оптимизация методом покоординатного спуска с переменным шагом А={10-1,10-2,10-13,10-4,10-5}, уменьшающимся каждый раз, когда при предыдущем значении шага получен локальный минимум.
В результате подбора коэффициентов а0, а, Ь путем решения задачи £^шт, получили следующий набор решений (рис. 5, табл. 1).
Таблица 1. Оценка качества экстраполяционных свойств моделей на основе ОДУ
+a2 (0,2 + 0,0251)2 + b _y(t) + b2 .y(t)2
и т. д.
Используем следующий набор подпрограмм.
1. Линейная одномерная интерполяция для вычисления x(t):tg{t;|
2. Численное решение обыкновенного дифференциального уравнения (ОДУ) методом Эйлера с небольшим шагом:
( m n ^
m n Итера- ций S y(10) y (10) - y*(10) inn p o/n
y (10)
0 0 17 0,11249 0,33933 43,24
0 1 8441 0,0003887 0,69345 15,98
1 0 13917 0,0090426 1,09493 83,14
1 1 70654 0,0000426 0,65809 10,07
1 2 66418 0,0000426 0,65809 10,07
2 1 58456 0,000452 0,65880 10,18
2 2 53054 0,000452 0,65880 10,18
Данные решения позволяют сделать достаточно адекватный прогноз развития системы (рис. 5). Важно даже не то, что его относительная погрешность меньше, чем раньше, а именно адекватность: система развивается примерно так, как получено в результате моделирования, тогда как на рис. 3 и 4 многие экстраполяции неадекватны. Кроме того, в последнем случае существует простой и понятный критерий выбора прогноза ¿, а в двух предыдущих случаях выбрать порядок интер-
Рис. 5. Результаты экстраполяции на основе решения ОДУ без учета взаимного влияния состояния системы и управления:
1) m=0, n=1; 2) m=1, n=1
Рис. 6. Результаты экстраполяции на основе решения ОДУ с учетом взаимного влияния состояния системы и управления: 1) п=1, т=1; 2) п=2, т=2; 3) п=1, т=2; 4) точное решение
полирующего полинома, оценивая сумму квадратов отклонений, не удается.
Рассмотрим модель, учитывающую взаимное влияние управления и состояния системы вида
т ”
,у-х(/Уу(/)’. Внесем соответствующие
& ,=0 ; = 0
изменения в программу. Так, блок решения ОДУ примет вид
for (double t=0; t<tk; t+=dt) // tk - конечное время, 0...10 {
tmp=0; x=xt (t);
for (i=0; i<=im; i++) for (j=0; j<=in; j++)
tmp += ak [i][j]*ipow (x, i)*ipow (y, j); y += dt*tmp;
if (у<0) у=0;
},
где 1ро^ - функция возведения в целочисленную степень, %1 - функция одномерной интерполяции.
Произведя численные эксперименты, получим результаты прогнозов на основании рассмотренных моделей (табл. 2).
Таблица 2. Оценка качества экстраполяционных свойств моделей на основе ОДУ сучетом взаимного влияния управления и состояния системы
m n S у(10) С с с О1 ^У 1 О1
У (10)
0 0 0,11249 0,33933 43,24
1 1 0,0000369 0,65583 9,69
1 2 0,0000031 0,66706 11,57
2 2 0,0000131 0,66906 11,91
СПИСОК ЛИТЕРАТУРЫ
1. Дзюба С.А. Модели управления подсистемами предприятия в сфере среднего бизнеса и их инструментальное обеспечение: автореф. дис. ... д-ра экон. наук. - Иркутск, 2011. - 46 с.
2. Мицек Е.Б. Эконометрическое моделирование инвестиций в основной капитал экономики России: автореф. дис. ... д-ра экон. наук. - Екатеринбург, 2011. - 49 с.
Интересно, что программа не находит в качестве оптимального верное решение a0i= 1, а02=ац=—1 при всех остальных a=0. Квадратичная ошибка в этом случае 6=0,035783, хотя погрешность прогноза у(10)=0,57121 составляет всего 4,46 %. Адекватность моделей сохраняется во всех случаях (рис. 6).
Заключение
Нельзя считать, что описанные неудачные попытки доказывают невозможность удачных аппроксимаций и экстраполяций временны?х трендов в эколого-экономическом моделировании с использованием традиционных и широко распространенных методов. Однако это хорошая иллюстрация того, что использование в качестве основы моделей обыкновенных дифференциальных уравнений может привести к качественному росту прогнозов и, следовательно, принимаемых решений.
3. Миролюбова А.А. Методология моделирования инвестиционного процесса в реальном секторе экономики региона: автореф. дис. ... д-ра экон. наук. - Иваново, 2012. - 33 с.
4. Лосев К.С. Мифы и заблуждения в экологии. - М.: Научный мир, 2010. - 224 с.
5. Арнольд В.А. Теория катастроф. - М.: Наука, 1990. - 128 с.
Поступила 26.03.2012 г.