CC BY
25
Вычислительные технологии
Том 22, Специальный выпуск 1, 2017
Моделирование кинетики химических реакций с использованием схем вариационного усвоения данных
А. А. Гришина1'2'*, А. В. Пененко1'2
1 Новосибирский государственный университет, Россия
2Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия
*Контактный e-mail: a.a.grishina17@gmail.com
Работа посвящена сравнению алгоритмов вариационного усвоения данных 3DVar в неявной модификации и метода 4DVar со слабыми ограничениями, примененных к задаче химической кинетики, описываемой системой уравнений Роберт-сона. Для этой цели сначала была решена прямая задача, и в ее численное решение внесен шум. Полученные данные о концентрации одного из веществ использованы в качестве "данных измерений" для алгоритма усвоения данных. Определение значений концентраций всех учитываемых в модели веществ и протекающих процессов с течением времени с использованием алгоритмов типа 3DVar и 4DVar проанализировано как теоретически, так и численно. Оказалось, что из-за специфики выбранной постановки более простой алгоритм 3DVar в неявной модификации показал лучшие результаты, чем стандартный алгоритм 4DVar.
Ключевые слова: вариационное усвоение данных, 3DVar, 4DVar, химическая кинетика, система Робертсона, Quasi Steady-State Approximation.
Введение
Динамические задачи совместного использования математических моделей и данных измерений для восстановления неизмеряемых характеристик исследуемого процесса относятся к классу задач усвоения данных. Усвоение данных получило широкое применение в прогностических моделях метеорологии и океанологии, где изучаются разномасштабные быстро изменяющиеся во времени процессы. В [1] можно проследить краткую историю развития методов усвоения данных в задачах прогнозирования основных параметров погоды.
В рассматриваемых задачах реконструкция функций состояния, описывающих поведение процесса, осуществляется путем его прогнозирования на основе априорных предположений, порождающих математическую модель, и последующей коррекции получаемого прогноза решения с учетом поступающих в реальном времени данных измерений. Так как в реальных приложениях данные поступают непрерывно, в задаче усвоения данных приходится решать последовательность обратных задач. В каждой такой задаче реализуются этапы получения численного решения в каждый момент времени: получаемое в процессе математического моделирования численное решение называется
© ИВТ СО РАН, 2017
фоновым прогнозом (далее просто прогнозом), поступающие данные именуются данными измерений, а окончательное решение называют анализом, причем последнее является оптимальным в некотором смысле решением текущей обратной задачи из последовательности. Данная терминология в английском эквиваленте применена, например, в [2].
Особенностью задач усвоения данных является необходимость учитывать как погрешности технологии моделирования, так и погрешности поступающих данных. Возникающие задачи минимизации соответствующих целевых функционалов стимулировали появление различных подходов к решению задач усвоения данных. Это методы взвешенных наименьших квадратов Гаусса; на их основе — метод фильтрации Калма-на и совокупность из, по меньшей мере, семи его ансамблевых модификаций [3, 4]; вариационные методы на основе принципа Эйлера-Лагранжа в формулировках с жесткими и слабыми ограничениями для учета моделей процессов, априорных сведений и доступных данных мониторинга [5]; методы релаксации Ньютона для прямого учета в моделях процессов невязок между рассчитываемыми и измеряемыми значениями функций состояния. Сравнение и краткое изложение стохастического и вариационного подходов представлено, например, в [6]. Существенным продвижением было использование методов сопряженных уравнений в задачах усвоения данных [7-9]. Ранее методы сопряженных уравнений использовались при решении обратных задач, что было отражено в [10].
Основные этапы в развитии теоретических основ и практических применений перечисленных методов усвоения, их сравнительный анализ с выявлением приоритетов и обсуждение перспектив научного поиска представлены в недавних тематических обзорах [2, 4, 11, 12], содержащих около тысячи ссылок.
В настоящей работе рассматриваются два метода: "неявный" аналог ЗВУаг [13-15] и 4ВУаг — методы вариационного усвоения данных в формулировках со слабыми ограничениями [11, 13, 16]. Здесь основной объект исследования — задачи химической кинетики в режимах прогноза и решения обратных задач по заданным целевым критериям. По содержанию статья представляет собой развитие работ [14-16]. Как и в [14], в состав модели процессов вводятся новые искомые функции, обозначающие в технике моделирования неопределенность.
Рассматриваемый "неявный" аналог ЗВУаг отличается от общепринятого метода тем, что вместо некоторой нормы невязки между анализом и прогнозом в целевой функционал включается норма управляющей функции, которая определяет величину этой невязки. Главное отличие методов ЗВУАИ, и 4ВУАИ, заключается в формировании "окна усвоения", содержащего различное количество данных измерений, участвующих в выражении целевого функционала: в методе ЗВУаг используются поступающие непосредственно в момент вычисления данные измерений, тогда как в 4ВУаг учитываются также данные измерений, поступившие в систему ранее.
Мы добавили в функционалы для схем ЗВУаг и 4ВУаг параметры регуляризации, которые обычно не применяются в задачах усвоения данных (по крайней мере, с апостериорными правилами их выбора), зато активно используются в теории регуляризации А.Н. Тихонова [17, 18]. Проблема выбора параметра регуляризации с позиций теории оптимального управления изучалась, например, в [19].
Работа алгоритмов усвоения данных демонстрируется на примере модели Роберт-сона [20] на выбранном конечном интервале. За "точное" решение было принято численное решение прямой задачи химической кинетики с известными начальными данными,
полученное по схеме QSSA [21, 22], а для решения задачи усвоения данных использовался линейно неявный метод Эйлера [20].
1. Постановка задачи усвоения данных
Изложим постановку задачи усвоения данных в химической кинетике. В силу специфики задач химии исследуемые процессы с большим разбросом скоростей реакций порождают жесткие системы с ограниченным интервалом предсказуемости процесса. Это накладывает ограничения на величину шага вычислений. А также в задачах химической кинетики применимо понятие источника веществ, которое будет использоваться в технологии построения схем вариационного усвоения данных.
Рассмотрим временной интервал [Ьа, Ьь]. Предположим, что исследуемый процесс описывается функцией состояния у(£) Е Кга, где п — конечное число элементов вектора состояния. Ее изменение во времени описывается задачей Коши
^ = Б(у(1)), 1Е [1а,и], (1)
у(*«) = V, (2)
где Б : Кга ^ Кга — оператор модели; V Е Кга — начальные данные. Прямая задача для этой системы формулируется так: по известным Б, V найти у(1), Ь Е \Ъа, Ьь]. Пусть в моменты времени ¿т поступают данные измерений:
Нт(у(1т)) = Фт = Ф™ + 6Фт, т = 0,...,М, Ьа <Р < и, (З)
где Фт — вектор данных измерений, поступающих в моменты ¿т; Ф^ — вектор точных измеряемых данных; дФт — вектор погрешностей данных измерений, поступающих в моменты 1т; Нт : Мга ^ Кга(т) — оператор измерений, устанавливающий соответствие между функцией состояния и данными измерений в момент их поступления. Размерность п(т) образа оператора измерений определяется количеством поступивших в "окно усвоения" измерений. Определим обратную задачу: по известным Б, Фт и некоторой информации о 5Фт для всех т = 0,... , М требуется найти у Ь Е \Ъа, 4].
В задаче усвоения данных требуется определить значение функции состояния у( ) на всем рассматриваемом промежутке [Ьа, Ь5] в ситуации, когда в реальном времени, т. е. в процессе вычисления решения задачи усвоения, поступают данные измерений Фт. Формально определим задачу усвоения данных как последовательность обратных задач, проиндексированных "текущим" моментом времени , в которой требуется определить у(Ь), Ь Е \Ъа, и], на основе Б, Фт и некоторой информации о 5Фт для всех т = 0,..., М (и), где М ^ *) = тах{т : 1т < и}.
2. Вариационные алгоритмы усвоения данных 2.1. Неявная модификация метода ЭЮуаг
Приведем этапы получения решения задачи усвоения с использованием метода типа ЗВУаг для постановки (1)-(З), в которой модифицируем уравнение (1) путем введения функции неопределенности ^(Ь) в правую часть уравнения (1):
(у(1)
(И
Б(у(1))+ ф), 1Е [Iа, и]. (4)
В химической кинетике, когда функция состояния y(t) — это вектор неизвестных концентраций веществ, функция неопределенности 'q(t) имеет смысл "искомых" источников веществ. Ее можно трактовать также как корректирующую добавку для последующей минимизации погрешностей модели.
Введем равномерную сетку ut = {Р, j = 0,... , N} и получим связь между анализом на новом шаге по времени ya(Р) и анализом на предыдущем шаге ya(Р-1), которая будет являться условием для минимизации целевого функционала. Мы воспользуемся линейно неявным методом Эйлера. Для его конкретизации перейдем от сеточных функций и операторов к функциям и операторам непрерывного аргумента. В предположении о том, что анализ ya(t), непрерывно зависящий от времени, удовлетворяет модифицированной модели (4):
dy(1 = S y (t))+ V(t), (5)
и непрерывно дифференцируем в окрестности значения анализа на предыдущем шаге по времени, разложим первое слагаемое S(ya(t)) в точке t3 по формуле Тейлора относительно предыдущего шага по времени Р-1:
dy§1 = S(ya(Р-1)) + J(ya(P-1)) {ya(P) - ya(P-1)) +
+V(P) + o{\\{ya(P) - ya(P-1))\\) . (6)
Здесь J(ya(P-1)) = —(y(-— — матрица Якоби правой части оператора модели S
"ja
(ya(t3-1)).
Далее аппроксимируем производную анализа по времени разностью назад с первым порядком, отбрасываем остаточный член в формуле (6). Заметим, что производную функции состояния по времени в левой части уравнения (4) можно аппроксимировать с помощью другой схемы, в таком случае уравнения далее хоть и изменятся, но это не повлияет на суть метода. В нашем случае получаем выражение для анализа и операторов, определенных для сеточного аргумента. Здесь и далее т = Р — Р-1 — шаг равномерной сетки ut.
y3 _ y3-1
Ya y = S(yi-1) + J {yi-1) {yi — yi-1) + v1. (7)
Тогда при заданном rf анализ y3a на новом шаге по времени выражается через предыдущее значение y3a-1 :
yi = yi-1 + тА1-1 (S (yi-1) + rj3) . (8)
Здесь и далее — единичная матрица размерности, равной числу компонент искомой функции состояния; А3-1 = (I — rJ (y3a-1)) . Таким образом, получаем линейно неявный метод Эйлера.
Неявная модификация метода 3DVar состоит в поиске условного минимума выпуклого функционала J3DVar на каждом новом шаге j по времени. Условием минимизации служит выражение для анализа, полученное ранее. С введением функции неопределенности целевой функционал принимает вид функционала Тихонова [17, 18]:
JsDVar(yi) = \\н3yi — W Iii-! + aslltf lß-1 -+ min, (9)
где а3 — параметр регуляризации; Н3, Ф] — оператор измерений и вектор данных измерений, введенные ранее. Они полагаются нулевыми, если в момент Р данных измерений нет. Норма вычисляется в пространстве Кга, а обратные к К, Q матрицы играют роль весовых матриц в скалярном произведении: (х, ъ)д-1 = x*Q-1z, х, ъ Е Кга, (х, ъ) = (х, ъ) 1. Матрица К имеет смысл матрицы ковариаций погрешностей измерений. Она описывает зависимость погрешностей результатов работы измерительных приборов друг от друга, если она известна, и задается диагональной (возможно, единичной), иначе. Q имеет смысл матрицы ковариаций погрешностей задания компонентов функции источников. Если соответствующие матрицы ковариаций будут известны заранее или вычислены, то их можно использовать при задании весовых матриц, что может благоприятно повлиять на качество решения — эти вопросы предполагается рассмотреть в дальнейшем. В проведенных численных экспериментах все весовые матрицы взяты единичными, а балансировка решения между фоновым прогнозом и точным воспроизведением данных измерений происходит за счет выбора параметра регуляризации.
Минимум функционала 3,ъвУаг отыскивается в процессе решения условной экстремальной задачи методом сопряженных множителей Лагранжа. Стационарная точка функции Лагранжа
Сзвъаг(V3, уа, у*) = \Нуа - Ф' 111-1 +а3\\Г? 11^-1 -
-м - у Г1 - гА1-1 (5 (уГ1) + ^ - , у**) (10)
есть условный минимум целевого функционала (9) при условии (8). Здесь у^ — сопряженный множитель Лагранжа.
Для нахождения стационарной точки функционала производные Фреше приравняем функции Лагранжа по ее переменным к нулю и получим систему
= 2 ( Н)* К-1 (Ну{ - Ф>) - у = 0,
ду3а
™УаГ =2 а:,Q-1 г? + т(А1-1)*у{ =0, (11)
д'цэ
9С^аГ = у{ - у{-1 - гАГ 1 (Б (у{-1) + *) = 0.
ду*
Тогда выражение для анализа на шаге ] имеет вид
(I+ а2А1-^-1 А-1)* (н*-* К-1 Н') у{ = (12)
т2
= уГ1 + тА^в {у-1- + -А>-^-1(А>-1У Н)*К-1Ф>.
а
В начальный момент времени можно положить у° = у | исходя из априорного предположения о значении функции состояния при Ь = ¿а.
Изложим традиционный метод ЗВУаг. На каждом шаге ] = 1,... , N по времени метод предполагает реализацию двух этапов:
1. Вычисление фонового прогноза у3р на новом шаге по времени ] по имеющейся математической модели при известном вычисленном значении анализа на предыдущем шаге ( - 1) . В нашем случае на данном этапе предлагается отыскивать решение прямой задачи (1), (2) с помощью линейно неявного метода Эйлера (по аналогии с (8)):
у1 = и у-1) = уГ1 + тА^Б (уГ1) , (1З)
где U является оператором решения прямой задачи в формулировке без функции неопределенности ^(t).
2. Поиск минимума целевого функционала J^Var, отличного от J3DVar из (9):
Xovar(yi) = \Нyi - & III-! + ы - у%-1 mn (14)
yi
где В-1 — весовая матрица в скалярном произведении при невязке между фоновым прогнозом и анализом. Она имеет смысл матрицы ковариаций погрешности модели получения фонового прогноза y3p по отношению к точному решению на шаге j. Так как точное решение предполагается неизвестным, вычисление матриц ковариаций становится отдельной задачей. В статье, как и все весовые матрицы, примем ее единичной.
Заметим, что в "традиционном" методе 3DVar осуществляется поиск безусловного минимума целевого функционала J3,DVar. Здесь в явном виде присутствует шаг прогноза, на котором учитывается связь фонового прогноза y3p с анализом y3a-1 посредством математической модели, а фоновый прогноз y3p связан с анализом y3a по построению целевого функционала. В предлагаемой неявной модификации метода 3DVar в уравнение модели добавляется функция неопределенности ^(t), которая присутствует в выражении целевого функционала, а анализ y3a на новом шаге по времени связан со значением анализа y3a-1 на предыдущем шаге по времени при помощи оператора U решения прямой задачи с включенной функцией неопределенности. Последнее является условием для минимизации целевого функционала J'sDVar. Неявная модификация 3DVar эквивалентна "традиционному" методу 3DVar, если весовая матрица В в "традиционном" 3DVar будет вычислена по формуле
в = -A3-1Q-1(A3-1) *. (15)
аз
Чтобы это показать, приравняем производную Фреше по y3a к нулю:
2 (Н3)*R-1 (Н3y3a - + 2В-1 (y3a - y3p) = (16)
разделим обе части уравнения на 2 и выразим y3p из (13):
(I + В (Н3)*R-1H3) y3a = = y3a-1 + rA3a~lS (y3a-1) + В (Н3)*R-1&. (17) Сравнивая (12) и (17), получаем (15).
2.2. Метод 4БУвг
В отличие от методов типа ЗВУаг, метод 4ВУаг подразумевает вычисление как фонового прогноза, так и анализа не в одном моменте времени, а на некотором временном интервале, образующем "окно усвоения".
В выражение для целевого функционала включаются данные измерений, поступающие в моменты }, к = к°,... ,кш, 1° < < 1к™ = 1м(и\ где т — ширина "окна усвоения", от которой зависит, какое количество данных измерений будет использовано при минимизации целевого функционала; М(1 *) — последнее ожидаемое измерение из (3), данные о котором придут в процессе счета. При фиксированном значении ширины "окна" т моменты времени внутри "окна" обозначим V : < V < 1к™, ] = 0,... ,т.
При получении новых данных измерений Фк^+1 в процессе вычисления "окно усвоения" перемещается, и в результате анализ, вычисленный до момента Ьк1, становится решением задачи усвоения в моменты Р : Ько < Р < Ьк1; начальный учитываемый момент поступления данных измерений обновляется: Ько = Ьк1, последний — добавляется к "окну": 1к™ = Ьк™+1. Все соответствующие значения анализа, фонового прогноза и другие величины, зависящие от шага по времени, изменяются соответственно. Таким образом, сдвиг "окна усвоения" происходит с перекрытием.
Следует отметить, что как при использовании метода типа ЗВУаг, так и для 4ВУаг анализ необходимо найти не только в точках поступления данных измерений, но и в промежутках времени между этими точками. Введение "окна усвоения" позволит решить задачу минимизации для упомянутых промежутков. В случае метода типа ЗВУаг решение изменяется только в точках измерений, а между ними сохраняет ту же динамику, что наблюдается у фонового прогноза.
Приведем этапы поиска анализа методом 4ВУаг при фиксированном "окне усвоения" ширины т до его перемещения, используя обозначения из подразд. 2.1.
1. Вычисление фонового прогноза во всех точках, принадлежащих "окну усвоения". Так как "окна усвоения" перекрываются, то для части нового "окна усвоения", которая входила в предыдущее "окно усвоения" (обозначим эти шаги индексами ] = 1,... ,тр) в качестве фонового прогноза выбирается анализ, полученный на предыдущем шаге, а значения фонового прогноза на оставшейся части "окна усвоения" вычисляются с помощью решения соответствующей прямой задачи:
ур = и у-1) = у>р-1 + тА3,-^ у-1) , + 1,...,т, (18)
где АГ = (I -т. (уГ) )-1 ;3 (ур-1) ^ .
В предположении о том, что фоновый прогноз ур недостаточно точен, вводится вариация функции состояния бу3, связывающая анализ у3а с фоновым прогнозом ур:
у а = ур + йу3, 3 = 0,..., т. (19)
Получим зависимость между вариацией функции состояния на последующем шаге по времени и значением вариации на предыдущем шаге, для чего снова рассмотрим дифференциальную постановку задачи.
Пусть фоновый прогноз ур^) удовлетворяет рассмотренному ранее уравнению (4) при "ц(Р) = 0:
,-1 т, 1) дБ(ур(Р-1))
р(V
~РГ = Б (ур
Б(ур(1)), 1Е [1а, и]. (20)
Анализ и фоновый прогноз связаны соотношением (19), тогда
Б(ур(^ + бу(г)) + ф), 1е [га,и]. (21)
(1(ур(1) + 8у(1))
<к =Б (ур(
В случае если правая часть Б нелинейна, а каждая компонента вектор-функции состояния непрерывно дифференцируема в окрестности фонового прогноза, разложим правую часть (21) по формуле Тейлора при фиксированном так же, как делали на этапе 1 метода типа ЗВУаг:
Б (ур(г) + бу(г)) = Б(ур(1)) + 3(ур(1)) • 5у(Р) + о(\\бу(т . (22)
Вычитая из (21) выражение (20) с учетом (22), получаем Му(Ь)
3(ур(1)) ■ ¿у(*) + ^) + 0(||*у(*)||) . (23)
Далее аппроксимируем левую часть (23) с первым порядком по времени; отбросим остаточный член о(||бу(£)||) , все переменные в правой части (23) возьмем при Ь = Р :
¿у.3 _ ¿уЗ-1
—-— = J(yl) ■ 8у3 + V3, 3 = 1,...,т. (24)
2. Для вычисления анализа на текущем "окне усвоения" необходимо найти минимум целевого функционала оУаг. Он зависит от вариации функции состояния в начальный момент и функции неопределенности во всех последующих точках из "окна усвоения".
Уаг (6у°, V1, . ., V™ ) = = £ 01 Н3(у3а) - *3ЦЪ-1 +^4Цг13Ц1-1 ] + ъЦбуОЦЪ-! ШШ . (25)
3 = 1 7 У1'
Здесь а4 — параметр регуляризации, наличие которого в соответствии с теорией Тихонова [17, 18] позволяет уменьшать либо невязку между анализом и фоновым прогнозом при его увеличении, либо невязку между анализом и данными измерений в случае его уменьшения. Функционал записан унифицированно для всех моментов времени из "окна усвоения", при этом, если в момент ] данные измерений не поступали, оператор Н3 и вектор измерений Ф3 полагаются нулевыми.
Минимум функционала отыскивался, как и в алгоритме типа 3ВУаг, методом сопряженных множителей Лагранжа. Приведем выражение функции Лагранжа для функционала (25) с ограничениями (24), учитывая (19):
¿4 В.аг (6у°, V1, . ., V™ , у* ,..., у?) = £ 01 Н(у3р + бу3) - ФЦЪ-г + ] +
=1
11бу°111-1 - £(6у3 - А3(бу3-1 + т3), у* ), (26)
3 = 1
где у — сопряженные множители Лагранжа. Функции минимума функционалов (10), (26) в задачах усвоения данных определяются как 3ВУаг/4ВУаг-анализ исследуемой системы в соответствующих областях определения моделей процессов и данных измерений в заданных схемах усвоения.
В настоящей работе рассматривается итерационный поиск минимума функции Лаг-ранжа. К началу проведения итераций на всем "окне усвоения" вычисляется фоновый прогноз у3р, ] = 0,... ,т, по формуле (18). При работе с "окном усвоения" в первый раз фоновый прогноз у3 задается равным функции
состояния у| 1=1а из априорных предположений или начальным данным измерений Ф°.
В рамках каждой итерации считаются известными значения бу°, г/3, ] = 1,... ,т, вычисленные на предыдущей итерации или взятые из начального приближения, которое положим нулевым. В первую очередь решается прямая задача: отыскиваются численные значения анализа во всех точках "окна усвоения" по формулам (19), (24).
Отметим, что при выполнении условия (24) производная Фреше функции Лагранжа по сопряженным множителям у1, ] = 1,... ,щ, принимает нулевое значение.
Затем решается сопряженная задача, состоящая в поиске значений сопряженных множителей Лагранжа. Для их получения продифференцируем функцию Лагранжа по переменным бу3, ] = 1,... ,щ, входящим в ее выражение, и приравняем результат к нулю:
= -у* + №У?*+1 + 2Н*П-1(Ну + ду3) - V) = 0, (27)
ды птт* тз-1(и(
д 6уь
у: + 2Н*К-1(Ну + ду™) - V™) = 0. (28)
Сопряженная задача решается в обратном времени: согласно (28) вычисляется у™, затем из (27) рекурсивно выражаются все значения у1, ] = (щ - 1),..., 1.
Затем с использованием градиентного метода происходит поиск ду°, г/3, ] = 1,... , уменьшающих значение функционала (25), в направлении, противоположном производным по этим переменным от функции Лагранжа:
= 2а4В -1ду° + (А1Уу*, (29)
= 2а,Я-1г13 + (Ар)* у**, ¿ = 1,...,ь>. (30)
В качестве реализации градиентного алгоритма был использован метод сопряженных градиентов из свободно распространяемой библиотеки ОБЬ [23].
После выполнения критерия остановки итераций (ортогональность текущего направления поиска минимума и градиента функции Лагранжа (29), (30)) "окно усвоения" перемещается. До получения следующих данных измерений алгоритм работает с целью вычисления фонового прогноза согласно (18).
3. Усвоение данных для системы Робертсона
3.1. Система Робертсона. Задача химической кинетики
Для получения численных результатов применения алгоритмов усвоения данных рассмотрим задачу химической кинетики, описываемую системой Робертсона [20]. Эта система возникла при описании химических реакций, протекающих в клетках некоторых растений [24]. В системе три субстанции участвуют в трех химических реакциях:
( 1() -к1У1® + Ьт ®Уг®, (31)
(И
Ь У1(г) - к2 у2(г) у3(г) - кзУ22(г), (32)
(У2(1) _ „. ^ „. т „. N 7, „. 2
к3 У22(1), 1е [Iа; 4] = [0;1], (33)
2,
с начальными условиями и значениями параметров:
У1(0) = 1, У2 (0) = 0, у3(0) = 0, (34)
к1 = 0.04, к2 = 104, к3 = 3 • 107. (35)
Здесь функция состояния y(t) представлена трехкомпонентным вектором неизвестных концентраций уг(Р),г = 1, 2, 3, реагирующих веществ; кг — константы скоростей химических реакций. Примем константы скоростей реакций и время обезразмеренными, как это сделано в [20]. Так как наибольшее отношение констант скоростей реакций имеет порядок 109, система считается жесткой.
Для проверки алгоритмов усвоения смоделируем данные измерений. Для этого решим прямую задачу, состоящую в определении концентраций реагентов в системе Ро-бертсона (31) на всем промежутке времени при известных начальных данных (34), а затем аддитивно введем "шум" в полученное численное решение.
Введем равномерную сетку на рассматриваемом интервале [0; 1] с шагом т = 10-4 и получим численное решение прямой задачи с помощью метода квазистационарных аппроксимаций QSSA (от англ. Quasi-Steady-State Approximation) [21, 22]. Этот метод одним из первых стал использоваться для решения задач химической кинетики. Алгоритм QSSA применяется к системам ОДУ, правую часть которых можно представить в виде суммы продукции Р и деструкции D:
dyjt) dt
S(y(t)) = Р(y(t)) - D(y(t))y(t), t E [ta; tb].
(36)
Если концентрация вещества увеличивается исходя из стехиометрических уравнений, то соответствующие слагаемые попадают в продукцию, иначе — в деструкцию. Для системы Робертсона они имеют следующий вид:
D(y(t)) • y(t)
'h 0 0
00
h уз&) + hy2(t) 0 00
Р (y(t))
'h У2СО Уз^У h y\(t) h yl (t)
fyi(t)\
I y^(t) I Ww 'my
Р (t) Mt),
fD\ (t) уЩ D ( ) ( )I D(t) V3(t)J
(37)
(38)
Метод QSSA для вычисления значений концентраций веществ на новом шаге по времени при известных значениях концентраций на предыдущем шаге имеет вид
J+i
У3г + Рг Г,
рЦD,
exp (-Dl т)у3 +
если Dlr < £m если D3r > £m
1 — exp (—D\ т)
D
Р r,
иначе.
(39)
104};
Здесь у3 — концентрация вещества г Е {1, 2, 3} в момент времени Р, ] Е {1,... Р3 = Рг(у3); О3 = Ог(у3); ега;п — малое число, вводится для проверки, что вычисления не содержат деления на близкое к машинному нулю число; етах — ограничивающая константа, при превышении которой правая часть будет выражать стационарное значение соответствующего уг. Численная схема (39) является монотонной и дает неотрицательное решение при неотрицательных входных данных.
Численное решение прямой задачи методом QSSA примем за "точное решение" задачи усвоения; в качестве "данных измерений" возьмем решение прямой задачи и внесем в значения концентрации одного из веществ — "шума" — по формуле
=11 +
100%
е) ул^),
т = 1,
М,
(40)
где Фт — построенные "данные измерений"; у1 — компонента имеющегося "точного решения ; 4 — нормально распределенная случайная величина с нулевым средним и единичной дисперсией; р — величина "шума" в процентах, примем р = 5%. В качестве начальных данных использовались первое "измерение" Ф1 для первого вещества и нулевые значения концентраций остальных веществ.
3.2. Результаты численных экпериментов
Для системы Робертсона поставим задачу усвоения данных следующим образом: необходимо определить значения концентраций всех реагентов на всем промежутке времени по поступающим "данным измерений" об одном веществе. В реальных условиях это типичная ситуация, когда из всей совокупности искомых функций состояния система мониторинга может обеспечить измерения только отдельных ее компонент, а технология моделирования использует их для восстановления эволюции исследуемого процесса.
Исходя из начальных значений концентраций в постановке прямой задачи для системы Робертсона, для начала всех трех реакций необходимо, чтобы концентрация первого вещества была положительной, поэтому целесообразно выбрать для усвоения именно его. Вследствие монотонности численной схемы (39) при начальных концентрациях (34) с добавлением шума (40) концентрация вещества 1 будет неотрицательной.
Реализации методов 3ВУаг и 4ВУаг были проверены на отсутствие ошибок путем усвоения "точного" решения и вычисления фонового прогноза по модели получения "точного" решения. Особый интерес представляют тесты для исходной постановки, когда в качестве основной модели вычисления в методах усвоения был выбран линейно неявный метод Эйлера. В тестах при увеличении параметра регуляризации а невязки между анализом и данными измерений возрастали, а между анализом и фоновым прогнозом уменьшались. При уменьшении параметра регуляризации возникал противоположный эффект: невязки между анализом и данными измерений уменьшались, между анализом и фоновым прогнозом — возрастали, что соответствует результатам общей теории регуляризации А.Н. Тихонова [17, 18].
На рисунке приведены графики численного решения для первого, второго и третьего вещества задачи усвоения данных каждым из методов: фоновый прогноз, получаемый при использовании 4ВУаг; усваиваемые "данные измерений" и "точное" решение.
Для метода типа 3ВУаг выбран параметр регуляризации а3 = 10-8, так как этот алгоритм позволял получить более точный анализ по отношению к "точному" решению, когда невязки между "данными" измерений и анализом уменьшались при уменьшении параметра регуляризации. Для метода 4ВУаг были выбраны "окно усвоения" с пятью измерениями и параметр регуляризации а4 = 103. Данный алгоритм позволял получить более точное решение при уменьшении невязок между анализом и фоновым прогнозом. В проведенных численных экспериментах метод 4ВУаг позволил получить более сглаженное решение, нежели алгоритм типа 3ВУаг. Применение 4ВУаг позволило изменить не только значение в моменты усвоения данных, но и скорректировать решение до и после появления новых измерений.
Можно сделать вывод о том, что для малокомпонентной системы Робертсона модификация алгоритма типа 3ВУаг (более простая по сравнению с 4ВУаг) с нужным параметром регуляризации оказалась более эффективной в созданных условиях. В литературе, когда методы усвоения данных применяются к более сложным системам для решения прогностических задач, метод 4ВУаг часто показывает лучшие результаты,
Концентрация вещества у\
= 0.00002 и
4DVar results
--4DVar analysis
— forecast
■ ■ exact solution
• • measurements
0.4 0.S
Time
3DVar resulls
Концентрация вещества y2
- 3DVar analysis
■ ■ exact solution
4DVar analysis forecast exact solution
0.4 0.0
Time
O.OOQOOl— 0.0
Концентрация вещества уз
0.035
0.4 0.Û
Time
4DVar results
- 4DVar analysis
— forecast
— exact solution
0.4 0.Û
Time
Графики численного решения задачи при использовании методов 3DVar (слева) и 4DVar (справа): анализ (analysis) "неявного" 3DVar (аз = 10-8) , 4DVar (a4 = 103), фоновый прогноз (forecast — для 4DVar), "данные измерений" (measurements), "точное" решение (exact solution)
чем метод 3ВУаг. Применение общепринятого метода 3ВУаг, описанного в конце под-разд. 2.1, действительно дало худшие результаты в среднем по трем веществам, чем метод 4ВУаг. Для выяснения причин обнаруженных эффектов при использовании модификации метода 3ВУаг проанализируем модель с позиций теории обратных и некорректных задач. Суммируя уравнения (31), (32) и добавляя уравнение (33), получим
Ш + кзУ22(г) = - уг(1), у3(г) = к3у22(г),
У2(0) = У2, УЗ(0) = ¥3.
Если значения у\(Ь) заданы с погрешностью, то задача определения концентрации второго вещества эквивалентна решению уравнения Рикатти, правая часть которого получается как решение задачи численного дифференцирования. При условии заданных ^2,^3 решение задачи об определении концентраций у2(Ь), у3^) по заданной компоненте функции состояния у\(Ь) единственное, однако выбранная нами постановка является некорректной. Поэтому слишком точное приближение зашумленных данных, осуществляемое алгоритмом 4ВУаг, пагубно сказывается на точности решения.
Остается вопрос о проявлении высокой эффективности алгоритма "неявного" 3ВУаг в тех же условиях. Анализируя график первого вещества алгоритмом типа 3ВУаг, можно отметить, что производная анализа (наклон графика) между моментами усвоения соответствует производной точного решения, а так как ( ) определяется именно производной у\(Ь), то закономерно, что результат определения его концентрации будет лучше. С другой стороны, 4ВУаг при точном приближении данных, в отличие от 3ВУаг, увеличивает гладкость анализа. Как следствие, он существенно изменяет производную У\(¿), что приводит к существенным ошибкам в компоненте решения, соответствующей второму веществу. Вероятно, ситуацию можно исправить, если вместе с измерениями значений концентрации представлять также значения производной по времени. Однако с практической точки зрения получить такие данные измерений достаточно сложно.
Для решения подобных затруднений в теории регуляризации Тихонова предлагается выбирать параметр регуляризации в зависимости от уровня шума. Предполагаем рассмотреть этот вопрос в дальнейшем.
Выводы
Рассмотрены два алгоритма вариационного усвоения данных: метод типа 3ВУаг и алгоритм 4ВУаг. Оба алгоритма предполагают получение анализа — решения серии обратных задач, в совокупности составляющих задачу усвоения данных, — путем минимизации некоторого целевого функционала. В целевой функционал метода типа 3ВУаг входят нормы функций неопределенности, а также нормы невязок между анализом и данными измерений, взятыми только в момент усвоения данных. Целевой функционал метода 4ВУаг содержит аналогичные нормы функций неопределенности и разность между данными измерений и искомым анализом на фиксированном интервале времени — "окне усвоения". Это позволяет изменять траектории решения в фазовом пространстве искомых функций состояния не только в моменты поступления данных измерений, но и в промежутках между ними.
Работа алгоритмов продемонстрирована на примере решения задачи химической кинетики, описываемой системой Робертсона, обладающей свойством жесткости. Задача
усвоения данных состояла в поиске концентраций реагентов на рассматриваемом временном интервале при неизвестных начальных данных и "зашумленных данных измерений", смоделированных на основе численного решения прямой задачи методом QSSA. Последнее было принято за "точное" решение. В основе алгоритмов усвоения данных для получения анализа в методе типа 3DVar и фонового прогноза в 4DVar на новом шаге по времени при известных значениях функции состояния на предыдущем шаге использовалась иная модель — линейно неявный метод Эйлера. "Измеряемым" веществом выбран первый реагент. Значения концентраций второго и третьего реагентов удалось точнее определить с помощью метода типа 3DVar. Для "измеряемого" вещества погрешность численного решения задачи усвоения по отношению к "точному" решению при использовании метода 4DVar оказалось возможным уменьшить, но в проведенных экспериментах это негативно сказалось на точности восстановления концентраций неизме-ряемых веществ. Теоретическое исследование корректности алгоритмов и результатов численных экспериментов показало, что обнаруженные эффекты являются следствием свойств рассматриваемой задачи усвоения данных для модели Робертсона.
Благодарности. Работа выполняется при поддержке программы Президиума РАН и грантов 1.33П и 11.2П/1.3-3; поддержана проектами MK-8214.2016.1 и РФФИ 14-0100125.
Авторы благодарны профессору В.В. Пененко и к.ф.-м.н. Е.А. Цветовой за ценные обсуждения постановок задач и полученных результатов. Выражаем особую благодарность рецензентам за ценные замечания.
Список литературы / References
[1] Ипатова В.М., Шутяев В.П. Алгоритмы и задачи ассимиляции данных для моделей динамики атмосферы и океана. Научно-образовательный курс. М.: Мир, 2013. 29 с. Ipatova, V.M., Shutyaev, V.P. Data assimilation algorithms and problems for atmospheric dynamics and ocean models. Moscow: Mir, 2013: 29 p. Available at: https://mipt.ru/education/chair/mathematics/upload/99f/algsaasimilation.pdf (In Russ.)
[2] Navon, I.M. Data assimilation for numerical weather prediction: A Review. Berlin: Springer, 2009. P. 21-65.
[3] Kalman, R.E. A new approach to linear filtering and prediction problems // Transactions of the ASME J. of Basic Engineering. Ser. D. 1960. Vol. 82. P. 35-45.
[4] Bocquet, M., Elbern, H., Eskes, H., Hirtl, M., Zabkar, R., Carmichael, G.R., Flemming, F., Inness, A., Pagowski, M., Perez Camano, J.L., Saide, P.E., San Jose, R., Sofiev, M., Vira, J., Baklanov, A., Carnevale, C., Grell, G., Seigneur, C.
Data assimilation in atmospheric chemistry models: current status and future prospects for coupled chemistry meteorology models // Atmos. Chem. Phys. 2015. Vol. 15. P. 5325-5358.
[5] Sasaki, Y. A fundamental study of the numerical prediction based on the variational principle //J. Meteor. Soc. Japan. 1955. Vol. 33. P. 30-43.
[6] Красюк Т.В. Усвоение данных: конкуренция методов и проблема усвоения спутниковых наблюдений // Тр. ГМЦ РФ. 2012. Вып. 348. C. 43-60.
Krasyuk, T.V. Data assimilation: competition among methods and the problem of assimilation for satellite assimilation // Proc. of the Hydrometeorological Research Centre of the Russian Federation. 2012. Vol. 348. P. 43-60. (In Russ.)
[7] Пененко В.В., Образцов Н.Н. Вариационный метод согласования полей метеорологических элементов // Метеорология и гидрология. 1976. №. 11. C. 3-16.
Penenko, V.V., Obraztsov, N.N. A variational initialization method for the fields of meteorological elements // Meteorologiya i Gidrologiya. 1976. No. 11. P. 3-16. (In Russ.)
[8] Marchuk, G.I., Penenko, V.V. Application of perturbation theory to problems of simulation of atmospheric processes // Monsoon Dynamics (Joint IUTAM/IUGG Intern. Symp. on Monsoon Dynamics, Delhi, 1977) / J. Lighthill, R.P. Pearce (Eds). Cambridge Univ. Pres, 1981. P. 639-655.
[9] Le Dimet, F.-X., Talagrand, O. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects // Tellus A. 1986. Vol. 38A. P. 97-110.
[10] Марчук Г.И. О постановке некоторых обратных задач // Докл. АН СССР. 1964. Т. 156, № 3. С. 503-506.
Marchuk, G.I. Formulation of some inverse problems // Dokl. AN SSSR. 1994. Vol. 156, No. 3. P. 503-506. (In Russ.)
[11] Singh, K., Sandu, A., Bowman, K., Parrington, M., Jones, D.B., Lee, M. Ozone data assimilation with GEOS-Chem: a comparison between 3D-Var, 4D-Var, and suboptimal Kalman filter approaches // Atmospheric Chemistry and Physics Discussions. 2011. Vol. 11. P. 22247-22300.
[12] Seamless prediction of the earth system: from minutes to months / G. Brunet, S. Jones, P.M. Ruti (Eds). Geneva: WMO, 2015. 471 p.
Available at: http://library.wmo.int/pmb_ged/wmo_1156_en.pdf
[13] Пененко В.В. Вариационные методы усвоения данных и обратные задачи для изучения атмосферы, океана и окружающей среды // Сиб. журн. вычисл. математики. 2009. Т. 12, № 4. C. 421-434.
Penenko, V.V. Variational methods of data assimilation and inverse problems for studying the atmosphere, ocean, and environment // Sib. Zh. Vychislitel'noi Matematiki. 2009. Vol. 12, No. 4. P. 421-434. (In Russ.)
[14] Пененко А.В., Пененко В.В. Прямой метод вариационного усвоения данных для моделей конвекции-диффузии на основе схемы расщепления // Вычисл. технологии. 2014. Т. 19, № 4. C. 69-83.
Penenko, A.V., Penenko, V.V. Direct data assimilation method for convection-diffusion models based on splitting scheme // Comput. Technologies. 2014. Vol. 19, No. 4. P. 69-83. (In Russ.)
[15] Пененко А.В., Пененко В.В. Последовательные алгоритмы усвоения данных в моделях мониторинга качества атмосферы на базе вариационного принципа со слабыми ограничениями // Сиб. журн. вычисл. математики. 2016. Т. 19, № 4. C. 401-418. Penenko, A.V., Penenko, V.V., Tsvetova E.A. Sequential data assimilation algorithms in air quality monitoring models based on a weak-constraint variational principle // Numerical Analysis and Applications. 2016. Vol. 9, No. 4. P. 312-325.
[16] Elbern, H., Strunk, A., Schmidt, H., Talagrand, O. Emission rate and chemical state estimation by 4-dimensional variational inversion // Atmos. Chem. Physics. 2007. Vol. 7. P. 3749-3769.
[17] Тихонов А.Н. О решении некорректно поставленных задач и методе регуляризации // Докл. АН СССР. 1963. Т. 151, № 3. С. 501-504.
Tikhonov, A.N. About ill-posed problems solving and regularization method // Dokl. AN SSSR. 1963. Vol. 151, No. 3 P. 501-504. (In Russ.)
[18] Численные методы решения некорректных задач / А.Н. Тихонов, А.В. Гончарский, В.В. Степанов, А.Г. Ягола. М.: Наука, 1990. 230 c.
Numerical methods for solving ill-posed problems / A.N. Tikhonov, A.V. Goncharskiy, V.V. Stepanov, A.G. Yagola. Moscow: Nauka, 1990. 230 p. (In Russ.)
[19] Агошков В.И. Методы оптимального управления и сопряженных уравнений в задачах математической физики. М.: ИВМ РАН, 2003. 256 с.
Agoshkov, V.I. Optimal control and adjoint equation methods in problems of mathematical physics. Moscow: IVM RAN, 2003. 256 p. (In Russ.)
[20] Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999. 685 с.
Hayrer, E., Wanner, G. Solving ordinary differential equations. Stiff and differential-algebraic problems. Moscow: Mir, 1999. 685 p. (In Russ.)
[21] Pantea, C., Gupta, A., Rawlings, J.B., Craciun, G. The QSSA in chemical kinetics: As taught and as practiced / N. Jonoska, M. Saito (Eds). Discrete and Topological Models in Molecular Biology. Natural Computing Series. Berlin; Heidelberg: Springer-Verlag, 2014. 480 p.
[22] Jay, L.O., Sandu, A., Potra, F.A., Carmichael, G.R. Improved QSSA methods for atmospheric chemistry integration // SIAM J. Sci. Comput. 1997. No. 18(1). P. 182-202.
[23] GNU scientific library reference manual edition 2.2.1, for GSL Version 2.2.1 Available at: https://www.gnu.org/software/gsl/manual/html_node/ (accessed 25.10.2016).
[24] Сабинин Д.А. Физиология развития растений. М.: Изд-во АН СССР, 1963. 15 c. Sabinin, D.A. Physiology of plants growth. Moscow: Izd-vo AN SSSR, 1963. 15 p. (In Russ.)
Поступила в 'редакцию 13 февраля 2017 г.
Chemical kinetics modelling with variational data assimilation schemes
grishina, anastasiia A.1'2'*, penenko, alexey V.1,2 1 Novosibirsk State University, Novosibirsk, 630090, Russia
2Institute of Computational Mathematics and Mathematical Geophysics SB RAS, 630090, Russia * Corresponding author: Grishina, Anastasiia A., e-mail: a.a.grishina17@gmail.com
The paper addresses a comparison of the implicitly modified variational data assimilation 3DVar algorithm with the weak-constrained method 4DVar. Both methods are applied to a chemical kinetics problem based on the Robertson system. To model the measurement data for the Robertson system a direct problem is solved. Noise is accounted in the numerical solution and the real-time data assimilation of one reagent is carried out. The numerical solution for the direct problem is obtained with the help of quasi-steady state approximation, while the direct problem step in data assimilation algorithms is performed using linearly implicit Euler method.
In the paper an implicit analogue of 3DVar is considered. It differs from the conventional method so that the penalty functional contains the norm of control function instead of the norm of the discrepancy between analysis and forecast. The norm of control function determines this discrepancy. The principal distinction between 3DVar and 4DVar methods lies in the formation of an "assimilation window" which contains different number of measurements as a part of cost functional expression, in the 3DVar
© ICT SB RAS, 2017
method the measurements incoming immediately in the moment of computation are solely used, whereas in the 4DVar method the measurements, which have come earlier, are also taken into account. Identifying concentrations of all species and ongoing processes in time with the help of modified 3DVar algorithm and 4DVar are analyzed both theoretically and numerically. It was shown that the 3DVar algorithm demonstrated better performance for the considered conditions.
Keywords: variational data assimilation, 3DVar, 4DVar, chemical kinetics, Robertson system, Quasi Steady-State Approximation.
Acknowledgements. This research was partly supported by RAS Presidium program and grants I.33P and II.2P/I.3-3 and by projects MK-8214.2016.1 and RFBR 14-0100125.
Authors are thankful to Professor V.V. Penenko and to Dr. E.A. Tsvetova for valuable discussions on the results.
Received 13 February 2017
