Ассимиляция данных в имитационном моделировании экологических процессов методом минимизации корректирующих возмущений Текст научной статьи по специальности «Математика»

2017 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА Т. 13. Вып. 3

ПРИКЛАДНАЯ МАТЕМАТИКА. ИНФОРМАТИКА. ПРОЦЕССЫ УПРАВЛЕНИЯ

ПРОЦЕССЫ УПРАВЛЕНИЯ

УДК 539.3

А. Г. Топаж1, Е. П. Митрофанов2

АССИМИЛЯЦИЯ ДАННЫХ В ИМИТАЦИОННОМ МОДЕЛИРОВАНИИ ЭКОЛОГИЧЕСКИХ ПРОЦЕССОВ МЕТОДОМ МИНИМИЗАЦИИ КОРРЕКТИРУЮЩИХ ВОЗМУЩЕНИЙ

1 ООО «Бюро Гиперборея», Российская Федерация, 193312, Санкт-Петербург, ул. Подвойского, 40

2 Агрофизический научно-исследовательский институт, Российская Федерация, 195220, Санкт-Петербург, Гражданский пр., 14

Оперативная коррекция текущего вектора состояния динамической модели данными прямых или косвенных измерений в режиме реального времени представляет собой известную задачу теории автоматического управления и называется проблемой ассимиляции или усвоения данных. Математические модели в экологии требуют специфических подходов к этой проблеме по сравнению с традиционными методами, применяемыми в гидрометеорологии. Основное отличие состоит в том, что в этом случае наличествуют, как правило, лишь редкие наблюдения, причем над ограниченным числом косвенных характеристик. В работе предлагается оригинальный подход к ассимиляции данных для задач подобного рода. Он заключается в том, что в исходную систему вводится дополнительное слагаемое, отражающее роль неучитываемых в идеальной модели случайных внешних воздействий. Далее находится такая форма этих возмущений, которая минимизирует взвешенную сумму их интегральной мощности и нормы отклонений наблюдаемых и теоретических характеристик в точках измерения. При этом выбор весового коэффициента позволяет отразить сравнительную степень доверия теоретической модели или фактическим измерениям. Таким образом, формальная постановка задачи записывается в виде проблемы оптимального управления. Приведены примеры применения предложенного метода для нескольких простых моделей. В задаче о равномерном движении материальной точки получено аналитическое решение, дающее, тем не менее, нетривиальные результаты. Задача об ассимиляции данных в модели Лотки—Вольтерра решена численно. Показано, что в этом случае предложенный метод «минимального возмущения» приводит к наименее «травматичной» коррекции эталонной модели по сравнению с известными альтернативными подходами — адаптивной идентификацией параметров или подгонкой начального состояния. Библиогр. 7 назв. Ил. 2.

Ключевые слова: ассимиляция данных, минимизация корректирующих возмущений, экология, имитационное моделирование.

Топаж Александр Григорьевич — доктор технических наук; alex.topaj@gmail.com Митрофанов Евгений Павлович — аспирант; mjeka@bk.ru

Topaj Alexander Grigorievich — doctor of technical sciences; alex.topaj@gmail.com Mitrofanov Evgeniy Pavlovich — postgraduate student; mjeka@bk.ru

© Санкт-Петербургский государственный университет, 2017

A. G. Topaj1, E. P. Mitrofanov2

ASSIMILATION OF DATA IN THE IMITATIVE MODELING OF ENVIRONMENTAL PROCESSES BY THE METHOD OF MINIMIZING CORRECTIVE PERTURBATIONS

1 Bureau Hyperborea, 40, Podvoisky ul., St. Petersburg, 193312, Russian Federation

2 Agrophysical Research Institute, 14, Grazhdanskiy pr., St. Petersburg, 195220, Russian Federation

The on-line correction of the vector of a dynamic model state variable by direct or indirect measurements is a well-known problem of the theory of automatic control. It is called the problem of data assimilation. Ecological models require specific approaches to this problem in comparison to traditional methods used in hydrometeorology. The main difference is that there are only rare observations on a limited number of indirect characteristics. The paper presents an original approach to data assimilation for such cases. The main idea is inclusion an additional term to the initial system, which describes the random external influences that are not counted in the ideal model. Further, one can find a form of these perturbations that minimizes the weighted sum of their integrated power and the norm of deviations of the observed and theoretical characteristics at moments of measurement. Then the choice of weighting factor allows us for comparative confidence to be reflected between the theoretical model and actual measurements. Thus, the formal statement of the problem is written down as a problem of optimal control. Examples of application of the proposed method for several simple models are given. In the problem of the uniform motion of a material point, an analytical solution can be obtained which, nonetheless, yields rather interesting results. The problem of data assimilation for the Lotka—Volterra model has been solved numerically. It is shown that the proposed method of "minimal perturbation" leads to the least "traumatic" correction of the ideal model in comparison with the available alternative approaches; i. e., adaptive parameter identification or backdating adjustment of the initial state. Refs 7. Figs 2.

Keywords: assimilation of data, minimization of corrective disturbances, ecology, simulation modeling.

Введение. Динамическое имитационное моделирование и компьютерный эксперимент все глубже проникают в практику научного исследования не только в технических, но и в естественных дисциплинах. При этом возникает весьма непростой вопрос, касающийся взаимодействия результатов численного и данных натурного экспериментов. В качестве примера можно рассмотреть ситуацию с динамическими моделями в агроэкологии, описывающими продукционный процесс культурных растений на сельскохозяйственном поле. За более чем 30-летнюю историю этого направления качественный скачок претерпели не только компьютерная техника, но и технологические средства измерения. С появлением совершенно новых, дистанционных и даже космических средств мониторинга оказалось, что многие характеристики состояния системы «почва—растение—атмосфера» сегодня допускают прямое оперативное измерение с практически произвольно задаваемым уровнем временного и пространственного разрешения. Причем для этого зачастую не требуется никаких специальных усилий или дорогостоящих измерений, достаточно иметь доступ к частично или полностью открытым, оперативно обновляющимся базам данных дистанционного спутникового зондирования. Условно говоря, листовой индекс, коэффициент проективного покрытия или даже степень обеспеченности растений азотным питанием, прямо или косвенно определяемые из вегетационных оптических индексов, в принципе сегодня могут считаться такими же внешними входными данными для модели, как метеорологические параметры — температура или осадки [1]. В существующих моделях эти величины вычисляются на каждом шаге расчета, т. е. представляют собой

компоненты вектора внутренних переменных состояния моделируемой системы. Таким образом, возникает объективная проблема применения оперативно поступающих экспериментальных данных для корректировки модели в режиме «runtime», т. е. непосредственно в ходе расчета. Решение указанного вопроса отнюдь не так очевидно, как кажется на первый взгляд.

Так, в последнее время активно используется подход, связанный с корректировкой величины листового индекса посева в моделях продуктивности сельскохозяйственных растений по данным спутниковых снимков дистанционного зондирования Земли [2]. Не вызывает сомнения, что листовой индекс действительно может быть более-менее адекватно оценен по тем или иным индексам спектрального анализа изображения (NDVI и т. п.). Однако явная коррекция в реальной динамической модели агроэкосистемы только этого показателя для обеспечения достоверности прогнозных расчетов в течение оставшегося сезона вегетации совершенно недостаточна. Как правило, листовой индекс является не ключевой, а вычисляемой переменной состояния и пересчитывается на каждом шаге из биомассы побега. Тогда его коррекция на текущем шаге позволит, возможно, лишь более адекватно смоделировать фотосинтез на следующем шаге, после чего модель вновь вернется к исходному состоянию. Если же воспользоваться обратной зависимостью и пересчитать по листовому индексу биомассу побега, то возникает вопрос — а что делать с неоднозначно связанными с нею на всех шагах расчета другими динамическими переменными: биомассой корня, физиологическим временем, влажностью почвы и др.? Таким образом, примитивная подмена модельного значения измеренным может сработать только для очень простых малопараметрических моделей регрессионного характера, но вряд ли допустима для комплексных динамических моделей агроэкосистем.

Вообще говоря, вопрос оперативной коррекции текущего вектора состояния динамической модели прямыми или косвенными измерениями представляет собой известную задачу теории автоматического управления и называется проблемой ассимиляции или усвоения данных [3]. Математическая теория ассимиляции данных получила наибольшее развитие в приложении к задачам синоптики [4]. Безусловно, динамическое моделирование в экологии, в частности моделирование продукционного процесса растений, имеет ряд специфических особенностей по сравнению с прогнозом погоды. Прежде всего это гораздо больший разброс в семантике составляющих вектора состояния моделируемого объекта. В задачах метеорологического прогнозирования большая размерность вектора состояния определяется в первую очередь пространственным разрешением расчетной сетки, а набор отслеживаемых физических величин сравнительно невелик (температура и давление). В то же время модель агроэкоси-стемы может быть даже точечной (одномерной), однако число значимых переменных состояния у нее довольно велико (биомассы различных органов растения, концентрации питательных веществ и влаги в почве и растении и т. д.). Далее, несмотря на появление новых средств измерения, с математической точки зрения большинство экологических моделей все еще представляют собой ненаблюдаемые системы управления — прямой или косвенной оценке по имеющимся наблюдениям доступна лишь малая часть переменных состояния. И, наконец, существенной особенностью является временная разреженность поступления данных измерений по сравнению с шагом модели. Если в метеорологических приложениях уровень развития информационного обеспечения современных моделей глобальной циркуляции позволяет организовать поступление свежих сведений оперативных измерений чуть ли не на каждом шаге, автоматизация модельных расчетов в экологии и агроэкологии еще не достигла

подобного совершенства. Таким образом, математические модели в экологии зачастую требуют специфических подходов к оперативной коррекции модели данными измерений. В ряде работ рассматриваются альтернативные методы осуществления этой процедуры в практических задачах мониторинга и предсказания урожайности сельскохозяйственных посевов с использованием динамических моделей продуктивности [5]. В статье предлагается оригинальный подход к ассимиляции данных, основанный на методах решения вариационных задач в теории автоматического управления. Описываются примеры его применения для двух очень простых моделей, проводится сравнение с известными решениями.

Постановка задачи ассимиляции данных. Рассмотрим постановку задачи ассимиляции данных в общем виде. Есть динамическая модель в дискретном времени

Хк+1 = I (Хк , Щ, ) (1)

с вектором состояния х. С помощью (1) осуществляется рекурентный пересчет этого вектора в реальном времени с учетом управления и и неконтролируемых воздействий w. И в какой-то момент времени к = К нам приходит некоторое измерение величины, зависящей от х, т. е. имеем вектор у к. При этом правило вычисления у от вектора состояния х задано априори:

Ук = у(хк),

а размерность вектора у существенно меньше размерности вектора х, т. е., как правило, обратного однозначного отображения у в х не существует. Ясно, что в реальности ожидаемое для идеального случая равенство

Ук = у(хк)

не выполняется. Причин этого может быть несколько: ошибки измерений, неверная интерпретация правила вычисления у от х или, наконец, неточность самой модели. Причем источников такой неточности, в свою очередь, может быть много — неверное задание начального состояния хо, неточности во входных данных щ, Wk, к = 1,..., К, неопределенность величин параметров модели или влияние вообще неучтенных в модели эффектов и процессов. Однако модель такова, какова она есть, и, несмотря на текущее рассогласование с данными измерений, необходимо продолжать расчеты по ней далее для моментов времени к = К +1,..., N, возможно каким-то образом скорректировав текущий вектор состояния для лучшего удовлетворения данных наблюдений. Заметим, что речь в таком случае идет именно об оперативной корректировке самого вектора состояния модели, а не ее параметров (адаптивная идентификация) или закона управления (адаптивное управление). То есть сама модель остается неизменной, а мы просто подменяем текущую величину Хк на некое специальным образом скорректированное значение Хк = Х(Хк, ук, ук) и продолжаем расчет с данного состояния по той же модели по крайней мере до получения очередного измерения, требующего корректировки.

Проблема состоит в том, каким должен быть выбран алгоритм корректировки. Ясно, что он напрямую зависит от того, какой причиной исследователь объясняет себе полученное рассогласование. Если подозрение падает на неверное задание в модели тех или иных постоянных параметров, требуются их «переидентификация» и пересчет всей предыдущей траектории с уточненными значениями. Этот подход очень близок проблеме адаптивной параметрической идентификации. Предположение о неточности задания начальных условий влечет за собой решение той или иной

обратной задачи для их модификации. Примеры использования обоих подходов, равно как и их разнообразные комбинации в моделях природных и биологических систем, изложены в [5]. Ниже рассматривается метод, основанный на априорном предположении о том, что используемая модель структурно неполна для адекватного описания имитируемого объекта.

Метод минимизации корректирующих возмущений. В качестве возможного теоретически обоснованного и логически непротиворечивого метода решения задачи оперативной коррекции произвольной динамической модели предлагается следующий принципиальный подход.

Запишем уравнение динамики поведения моделируемой системы в виде

Xfc+i = f (xk, Uk, Wk) + ik, (2)

который практически полностью совпадает с исходной теоретической моделью исследуемого объекта за исключением добавления в правую часть вектора ik, отражающего интегральное влияние неучитываемых в опорной модели факторов и эффектов на каждом шаге рекуррентного пересчета вектора состояния. При этом, исходя из физики процесса, некоторые из составляющих вектора ik могут быть положены тождественно равными нулю (тем самым мы берем на себя смелость утверждать, что некоторые уравнения, составляющие эволюционный оператор f, описывают реальность абсолютно точно). Полагаем, что именно неучтенное влияние этих неизвестных и неопределенных факторов привело к тому, что очередные измерения оказались несогласованы с ожидаемыми величинами. Понятно, что количество вариантов возмущений, которые могли бы привести к полученному рассогласованию, огромно. Однако мы склонны доверять рассматриваемой модели и считаем, что в общем и целом она адекватно описывает реальность. Тогда можно поставить и в идеале решить следующую задачу: а какими должны быть минимальные по мощности неучтенные и неизвестные нам возмущения ik , при которых траектория динамики системы может породить вектор состояния в момент K, полностью или в большой мере соответствующий проведенному в этот момент измерению у к ? Разрешив такую задачу, т. е. найдя подобные «минимально обеспечивающие» возмущения ik, можно просто проинтегрировать измененную модель (2) на интервале времени от 0 до K и получить некое теоретически обоснованное состояние Хк, которое естественно принять в качестве опорной точки дальнейшего расчета модели в режиме реального времени. Таким образом, формальная постановка задачи может быть записана в виде проблемы оптимального управления:

Хк = Xk(ik) | k = K, xk+i = f (xk, Uk, Wk) + i(k), i(k) = arg nun J(x,i) = arg nun((1 - y)• || у(хк) - У к || +Y• || ik ||), (3)

где коэффициент 1 — y естественно интерпретировать как относительный показатель доверия измерениям по отношению к модели. При 7 ^ 0 решение будем искать в классе возмущений, обеспечивающих точное совпадение модели с измерениями в опорный момент времени, а выбор конкретного возмущения из такого класса, согласно критерию малости интегральной мощности, будет достигаться в пределе, подобно тому, как это происходит при регуляризации некорректных задач методом сингулярного разложения по малому параметру [6].

Примеры применения метода. Приведем примеры использования предложенного подхода для двух моделей динамики систем, записанных в непрерывном времени.

Пример 1. Задача о равномерном движении материальной точки.

Рассматривается простейшая модель движения материальной точки в предположении отсутствия внешних сил:

Х1 = Х2, Х1 |4=о = 0, (4)

¿2 = 0, Х2 |4=о = V.

Закон движения, определяемый этой тривиальной моделью, представляет собой равномерное прямолинейное перемещение точки с постоянной скоростью V, и тогда естественно ожидать, что в момент времени Ь = Т координата точки будет равна ¿1 (Т) = V ■ Т. Однако полученные измерения координаты, основанные, например, на данных телеметрии, показывают, что она составляет величину X. При этом явные измерения скорости, т. е. величины ¿2(Т), недоступны. Ставится задача оперативной коррекции вектора состояния модели для производства дальнейших расчетов в реальном времени при Ь > Т. Основной вопрос заключается в том, какое согласованное опорное значение в момент времени Т следует приписать переменной состояния ¿2, если переменная ¿1 будет скорректирована из условия полного доверия измерениям, т. е. ей будет приписана величина X. Рассмотрим решение такой задачи согласно предложенному алгоритму. Запишем измененную модель с учетом неопределенных внешних возмущений. Ясно, что для описываемого случая нужно добавить дополнительное слагаемое только во второе уравнение системы (4), так как первое уравнение представляет собой строгое тождество. Добавка же аддитивной составляющей в уравнение для ускорения в данном случае имеет прозрачную физическую аналогию — по нашему мнению, рассогласование результатов модели и эксперимента было вызвано действием неизвестной внешней силы, а не, например, неточностью задания начальных условий. Полагаем, что и координата, и скорость точки в момент начала отсчета были измерены точно, и им можно доверять. Таким образом, для модели с учетом возмущений имеем

¿1 = ¿2, ¿1 |4=о = 0,

(5)

¿2 = и, ¿2 |£=0 = V.

В полном соответствии с базовым принципом будем искать динамику возмущения как решение вариационной задачи, т. е. из условия минимизации функционала

Т

-1 = (1 - 7) ■ (¿1(Т) - XX)2 + 7 ■ Т3 ■ [ и2{Ь)йг

(здесь множитель

Т3

во втором слагаемом функционала добавлен как масштабный коэффициент для приведения обоих слагаемых к одной размерности). Поставленная задача синтеза оптимального управления может быть решена аналитически с помощью метода максимума Понтрягина.

Гамильтониан вариационной задачи записывается следующим образом:

Н = А1 ■ ¿2 + Л2 ■ и + 7 ■ Т3 ■ и2,

откуда сразу вытекает выражение для управления, доставляющего ему минимальное значение:

- Л2 и°Р* - 2 • 7 • Т3 •

Система дифференциальных уравнений для сопряженных переменных принимает вид

• дН

М = —т— =0,

дxl /РЛ

• ОН л (6)

\<2 — — —- — —А1

дХ2

о граничными условиями трансверсальности на правом конце интервала управления:

А1(Т) = 2-(1 - 7) ■ (¿1 (Т) - X),

Л2 (Т) =0. 1 )

Решение системы уравнений для сопряженных переменных (6) и (7) тривиально, и тогда получаем

А1(Ь) = 2 ■ (1 - 7) ■ (¿1 (Т) - X), А2(Ь) = 2 ■ (1 - 7) ■ (¿1 (Т) - X) ■ (Т - Ь), (8)

г^) = ■ (Х1(Т) - X) ■ (Т -1).

Таким образом, оптимальный закон управления представляет собой динамику возмущений в виде силы, линейно убывающей от времени и достигающей нуля в точке получения измерений Т. Подставив выражение для возмущения (8) в исходные уравнения для переменных состояния (5), можно аналитически разрешить систему, получив выражения для динамики переменных состояния как явные функции от времени:

6 ■ 7 2 ■ 7 Т

1-7 (Х1(т)-х).(т-г)3, (9)

6 ■ 7 ■ Т3

[,\ тг 1 -7 Х\{Т) -X 1 -7 , (ГГЛ у» № .42 - --^--2 • 7 • Г3 ' (-г'1(Т) - А ) • (Т - *) .

Записав первое уравнение системы (9) для момента времени Ь = Т, разрешаем выведенное алгебраическое выражение относительно неизвестной ¿1 (Т) :

Итак, «оптимально скорректированная» характеристика для координаты представляет собой средневзвешенную величину измеренного и модельного значений, где веса определяются введенным показателем 7. Для случая абсолютного доверия измерениям (7 = 0) она, естественно, оказывается равной XX.

Наконец, аналогичным образом для скорректированного значения скорости находим, что

т = 1/ - - 1~Л У-Т-Х= 7-(7-1) 3-(1-7) X

~2( ' 2 2 • (7 + 1) Т 2 • (2 • 7 + 1) 2 • (2 • 7 + 1) Т'

332 Вестник СПбГУ. Прикладная математика. Информатика... 2017. Т. 13. Вып. 3

Для случая полного доверия измерениям (7 = 0) окончательно имеем

Таким образом, наилучшим образом «согласованное» с полученными измерениями координаты значение скорости в точке измерения в предположении о минимальном отклонении реальной динамики системы от «эталонной» модели на всем промежутке интегрирования оказывается достаточно просто вычисляемым из начальной и средней скоростей движения. Численное исследование этого решения в среде имитационного моделирования системной динамики АпуЬс^с показывает, что для рассматриваемого случая линейно-убывающая от времени функция возмущения действительно имеет минимальную интегральную квадратичную норму по сравнению с другими тестируемыми формами возмущений (разовый импульс, сила постоянной интенсивности и т. д.), обеспечивающими переход системы из заданных условий в момент времени 0 в точку с измеренными координатами в момент времени Т.

Пример 2. Динамика двухвидового взаимодействия. Запишем классическую систему уравнений Лотки—Вольтерра, описывающую взаимодействие двух антагонистичных видов (хозяин—паразит, хищник—жертва и т. д.):

где а,1 и о,2 — суммарная приспособленность (коэффициент размножения) первого и второго видов; 61 и 62 — коэффициенты влияния межвидового взаимодействия на их динамику. Все параметры считаются положительными, так что знаки при соответствующих членах в системе уравнений (10) показывают, что в отсутствие взаимодействующего вида численность первой популяции имеет возрастающий характер, а второй — убывающий. В свою очередь, взаимодействие позитивно влияет на второй вид и отрицательно на первый.

В рассматриваемом случае будем интерпретировать переменные состояния Х1 и Х2 как суммарную биомассу растительности и численность единственного вида растительноядных животных в обсуждаемом ареале. Тогда соответствующая задача ассимиляции данных может быть поставлена следующим образом. Пусть известны начальные значения переменных состояния. Тогда система уравнений (10) может быть проинтегрирована численно, порождая, как известно, при определенных значениях параметров временную динамику мощностей популяций в виде сдвинутых по фазе друг относительно друга периодических колебаний. В момент времени Т с достаточной точностью становится известной фактическая биомасса растительности Х1 (оцениваемая, например, при помощи анализа космических снимков выбранной площади). В то же время численность животных прямому измерению не поддается. Требуется, исходя из отклонения измеренной и модельной величин Х1, оценить поправку к модельной величине Х2(Т), чтобы продолжить интегрирование модели с нового состояния. Применим к поставленной задаче предложенную выше методику, выбрав в качестве корректирующих малых возмущений, «ломающих» исходную модель, неизвестные процессы миграции растительноядных животных за пределы рассматриваемой области. Тогда измененная система будет записана в виде

Х1 = (01 — 61 • Х2) • Х1, Х2 = ( — 02 + 62 • Х1) • Х2,

(10)

Х1 = (01 — 61 • Х2) • Х1,

Х2 = ( — 02 + 62 • Х1) • Х2 + и,

(11)

а вариационная постановка задачи ассимиляции данных будет выглядеть следующим образом: найти такой закон изменения во времени интенсивности миграционных процессов и(Ь), который, обладая минимальной интегральной нормой, обеспечивал бы максимальную близость величин Х1 и Х1(Т).

В принципе, поставленная задача может быть сформулирована в общем виде (3) как стандартная проблема оптимального управления. Однако, в отличие от рассмотренного примера 1, ее оказывается невозможно разрешить аналитически, а характер получаемых уравнений таков, что даже нахождение численного решения для конкретных параметров может вызвать значительные затруднения. Потому в дальнейшем ограничимся рассмотрением вырожденного случая полного доверия измерениям. При этом критерий оптимальности вырождается в требование минимизации интегральной нормы управления

т

7 = ^ и2(г)аь, (12)

о

а к стандартным начальным условиям динамической системы (11) добавляется добавочное условие на конце интервала интегрирования для переменной Х1 , т. е. совокупность всех граничных условий записывается как

Х1(0) = Х10, Х2(0) = Х20, Х1(Т )= Х1. (13)

Гамильтониан вариационной задачи (11)—(13)

Н = А1 • (01 — 61 • Х2) • Х1 + Л2 • (—02 + 62 • Х1) • Х2 + Л2 • и + и2,

откуда оптимальное управление выражается как ир = — А2/2. Уравнения для сопряженных переменных имеют вид

А1 = —Л1 • (01 — 61 • Х2) — А2 • 62 • Х2,

• (14)

АЛ2 = А1 • 61 • Х1 — А2 • (62 • Х1 — 02),

а единственное дополнительное граничное условие трансверсальности (соответствующее уравнению для Х2 со свободным правым концом)

А2(Т ) = 0. (15)

Соотношения (11) с подставленным туда выражением для ир и (14) с граничными условиями (13), (15) представляют собой типичную краевую задачу для четырех переменных. Она может быть решена численно методом «пристрелки», а именно итерационным перебором неизвестных начальных значений для сопряженных переменных, где на каждом шаге осуществляются интегрирование получившейся задачи Коши и вычисление невязки значений на правом конце для заданных величин Х1 и А2 . В результате сразу имеем как явное выражение для закона управления, так и «скорректированную» динамику переменных состояния, обеспечивающую точное соответствие измеренной величины Х1 ее модельному значению в момент измерения. Для проведения соответствующих расчетов были построены модели как исходной системы (10), так и расширенной уравнениями для сопряженных переменных динамической системы (11)—(15) в среде имитационного моделирования AnyLogic. При этом решение краевой задачи методом «пристрелки» осуществлялось при помощи реализованной в этой среде специальной надстройки для параметрической идентификации исследуемых моделей.

В качестве опорного набора параметров для исследования были выбраны следующие значения: а1 = 0.15; а2 = 0.2; Ъ1 = 0.1; Ъ2 = 0.2; х10 = 2; х20 = 1; Т = 50. Ожидаемая величина динамической переменной Х1 в конце интеграла моделирования для «эталонной» системы должна быть равной 0.585. Измеренное значение, по которому осуществлялась коррекция модели, полагалось равным Х1 =2. Выборочные результаты решения поставленной задачи ассимиляции данных методом минимизации случайного возмущения приведены на рис. 1, где представлены динамика переменных состояния для исходной (без учета миграции) и скорректированной (с ми-

Time (conventional units)

Рис. 1. Временной ход переменных состояния для «эталонной» и «скорректированной» моделей двухвидового взаимодействия (A) и динамика минимального по мощности корректирующего возмущения (Б)

нимальной по мощности миграцией) моделей, а также временной ход сопряженной переменной \2 и определяемый ею закон оптимального управления. Нетрудно заметить, что в рамках предложенного подхода действительно оказывается возможным обеспечить достаточно существенное изменение конечного значения измеряемой переменной малым по норме корректирующим управлением (имеющим в данном случае смысл неучитываемых процессов эмиграции и иммиграции).

Следующим этапом исследования было сравнение различных методов ассимиляции данных для рассматриваемой простейшей модели двухвидового взаимодействия.

В частности, предложенный метод минимизации корректирующего возмущения сравнивался с альтернативными известными подходами — адаптивной идентификацией параметров и коррекцией начального состояния. В первом случае равенство модельной и фактической величин х\ (Т) обеспечивается калибровкой того или иного постоянного параметра исходной модели. Изучались все возможные случаи, т. е. параметризация производилась последовательно для каждого из параметров системы (10) — а\,а2, Ь2. Во втором случае «попадание в цель» обусловливается выбором адекватных начальных условий. Здесь, исходя из семантики задачи, варьированию подвергалась только величина х2о, поскольку именно невозможность точного измерения численности популяции первичных консументов служит в рассматриваемом случае причиной неполной наблюдаемости системы управления. Все вычисления также выполнены в соответствующих надстройках (экспериментах специального характера) среды имитационного моделирования AnyLogic.

Результаты проведенного сравнения подходов в терминах кривых роста отдельно для условной биомассы растительности и численности растительноядных представлены на рис. 2 (шкала по вертикальной оси обрезана из соображений выбора адекватного масштаба для удобства восприятия). Нетрудно заметить, что для

15 20 25 30 35

Time (conventional units)

Рис. 2. Кривые динамики переменных состояния xi (А) и x2 (Б) для разных способов коррекции модели

рассматриваемого случая метод «минимального возмущения» приводит к наиболее «щадящей» коррекции модели. Во-первых, соответствующие кривые динамики переменных состояния оказываются наиболее близки к «эталонной» модели, особенно на начальном периоде интегрирования. Во-вторых, отклонение конечной величины ненаблюдаемой переменной х2(Т ), которая будет принята за опорное начальное значение на следующих этапах расчета, также в таком случае наименьшее среди всех вариантов.

Заключение. Предлагаемый метод ассимиляции данных оперативных измерений в экологическом моделировании, как представляется, имеет ряд преимуществ по сравнению с имеющимися алгоритмами. Это не только наименее «травматичный» способ коррекции эталонной модели — в конце концов полученные при исследовании примера 2 результаты могут оказаться нехарактерными и невоспроизводимыми для другого набора исходных параметров или других динамических систем. Здесь, безусловно, открывается широкое поле для дальнейших исследований. Но в его пользу можно привести и соображения «идеологического» характера. В частности, важными позитивными особенностями рассматриваемого подхода являются независимость и полная изолированность расчетов для каждого нового поступления данных оперативных измерений. Действительно, для каждого нового периода интегрирования между измерениями при последовательном расчете модели заново решается задача минимизации корректирующего возмущения в пределах именно этого и только этого периода. Причем ни параметры модели, ни граничные значения переменных состояния во временных точках получения измерений не меняются — общая траектория системы остается гладкой, в то время как, например, использование аналогичного подхода при адаптивной идентификации либо приводит к тому, что порождающая модель имеет нестационарные, ступенчато меняющиеся от времени параметры, либо требует при поступлении новых данных полного пересчета всей предыдущей траектории, в том числе и с учетом всех уже проведенных ранее измерений.

Конечно, предлагаемый метод гораздо более трудоемок, чем формализованные и стандартизованные подходы адаптивной параметрической идентификации, фильтрации Калмана или решения обратной задачи определения начального состояния. В частности, он содержит проблему синтеза оптимального управления минимальной мощности, обеспечивающего максимальную степень соответствия измерений модельным расчетам. Ясно, что для сложных, плохо формализованных имитационных моделей реальных природных объектов о нахождении аналитического решения здесь речь не идет, а использование численных методов может потребовать значительных вычислительных и временных ресурсов. Однако мощности современных компьютеров и применение продвинутых средств автоматизации компьютерного эксперимента с динамическими имитационными моделями [7] гарантируют практическую реализуемость подхода. Огромное значение приобретает выбор адекватной семантики и формализации корректирующих возмущений — иными словами, четкое выделение автором модели ее потенциальных «тонких и чувствительных мест», в которых формальное введение неучитываемого в модели произвольного внешнего фактора позволит максимально «ненарушающим» способом осуществить оперативную корректировку непосредственно в процессе расчета. В любом случае не вызывает сомнений, что проблема ассимиляции данных в экологическом моделировании должна решаться не дилетантски, а на серьезном профессиональном уровне, т. е. с привлечением современного и продвинутого математического аппарата.

Литература

1. Брыксин В. М., Евтюшкин А. В., Хворова Л. А. Разработка программного комплекса оценки урожайности зерновых культур с использованием математического моделирования и данных дистанционного зондирования Земли // Материалы Всерос. конференции (с международным участием) «Математические модели и информационные технологии в сельскохозяйственной биологии: итоги и перспективы», 14-15 октября 2010 г. СПб.: АФИ, 2010. С. 76-80.

2. Клещенко А. Д., Найдина Т. А. Использование данных дистанционного зондирования для моделирования физиологических процессов растений в динамических моделях прогнозирования урожая // Современные проблемы дистанционного зондирования Земли из космоса. 2011. Т. 8, № 1. C. 170-178.

3. Шутяев В. П. Операторы управления и итерационные алгоритмы в задачах вариационного усвоения данных. М.: Наука, 2001. 240 с.

4. Kalnay E. Atmospheric modeling, data assimilation and predictability. Cambridge: Cambridge University Press, 2003. 341 p.

5. Olioso A., Inoue Y., Ortega-Farias S., Brisson N. Future directions for advanced evapotranspiration modeling: Assimilation of remote sensing data into crop simulation models and SVAT models // Irrigation and Drainage Systems. 2005. Vol. 19. P. 377-412.

6. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. М.: Наука, 1974. 223 с.

7. Medvedev SS., Topaj A. Crop simulation model registrator and polyvariant analysis // IFIP Advances in Information and Communication Technology (AICT). 2011. Vol. 359. P. 295-301.

Для цитирования: Топаж А. Г., Митрофанов Е. П. Ассимиляция данных в имитационном моделировании экологических процессов методом минимизации корректирующих возмущений // Вестник Санкт-Петербургского университета. Прикладная математика. Информатика. Процессы управления. 2017. Т. 13. Вып. 3. С. 326-338. DOI: 10.21638/11701/spbu10.2017.309

References

1. Bryksin V. M., Evtiushkin A. V., Khvorova L. A. Razrabotka programmnogo kompleksa otsenki urozhainosti zernovykh kul'tur s ispol'zovaniem matematicheskogo modelirovaniia i dannykh distantsionnogo zondirovaniia Zemli [Development of a software system for assessing crop yields using mathematical modeling and remote sensing data from the Earth]. Proceeding of Russian Conference "Mathematical models and information technologies in agricultural biology: results and prospects". October 14-15, 2010. Saint Petersburg, ARI Press, 2010, pp. 76-80. (In Russian)

2. Kleshchenko A. D., Naidina T. A. Ispol'zovanie dannykh distantsionnogo zondirovaniia dlia modelirovaniia fiziologicheskikh protsessov rastenii v dinamicheskikh modeliakh prognozirovaniia urozhaia [The use of remote sensing data to model the physiological processes of plants in dynamic crop forecasting models]. Sovremennyie problemy distantsionnogo zondirovaniia Zemli iz kosmosa [Modern problems of remote sensing of the Earth from space], 2011, vol. 8, no. 1, pp. 170-178. (In Russian)

3. Shutiaev V. P. Operatory upravleniia i iteratsionnye algoritmy v zadachakh variatsionnogo usvoeniia dannykh [Operators of control and iterative algorithms in problems of variational data assimilation]. Moscow, Nauka Publ., 2001, 240 p. (In Russian)

4. Kalnay E. Atmospheric Modeling, Data Assimilation and Predictability. Cambridge, Cambridge University Press, 2003, 341 p.

5. Olioso A., Inoue Y., Ortega-Farias S., Brisson N. Future directions for advanced evapotranspiration modeling: Assimilation of remote sensing data into crop simulation models and SVAT models. Irrigation and Drainage Systems, 2005, vol. 19, pp. 377-412.

6. Tikhonov A. N., Arsenin V. Ia. Metody resheniia nekorrektnykh zadach [Methods for solving ill-posed problems]. Moscow, Nauka Publ., 1974, 223 p. (In Russian)

7. Medvedev S., Topaj A. Crop simulation model registrator and polyvariant analysis. IFIP Advances in Information and Communication Technology (AICT), 2011, vol. 359, pp. 295-301.

For citation: Topaj A. G., Mitrofanov E. P. Assimilation of data in the imitative modeling of environmental processes by the method of minimizing corrective perturbations. Vestnik of Saint Petersburg University. Applied Mathematics. Computer Science. Control Processes, 2017, vol. 13, iss. 3, pp. 326-338. DOI: 10.21638/11701/spbu10.2017.309

Статья рекомендована к печати проф. Л. А. Петросяном. Статья поступила в редакцию 27 апреля 2017 г. Статья принята к печати 8 июня 2017 г.