CC BY
11
УДК 517.977, 519.688
DOI: 10.14529/ cmsel90402
АЛГОРИТМ ДИНАМИЧЕСКОЙ РЕКОНСТРУКЦИИ ВХОДОВ СТОХАСТИЧЕСКОГО ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ: НАСТРОЙКА ПАРАМЕТРОВ И ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ*
© 2019 Л.А. Мельникова, В.Л. Розенберг
Институт математики и механики им. Н.Н. Красовского УрО РАН (620990 Екатеринбург, ул. С. Ковалевской, д. 16)
E-mail: [email protected], [email protected] Поступила в редакцию: 16.04.2019
Задача реконструкции неизвестных входов стохастического дифференциального уравнения исследуется с позиций подхода теории динамического обращения. Рассматривается постановка, в которой одновременное восстановление возмущений в детерминированном и стохастическом членах уравнения проводится на основе дискретной информации о некотором количестве реализаций случайного процесса. Задача сводится к обратной задаче для системы обыкновенных дифференциальных уравнений, которым удовлетворяют математическое ожидание и ковариационная матрица исходного процесса. Разработан программно-ориентированный алгоритм решения, основанный на конструкциях теории позиционного управления с моделью; получена оценка его точности относительно количества доступных измерению реализаций. Предложена программная процедура настройки параметров алгоритма для получения наилучшего результата аппроксимации различных возмущений, удовлетворяющих априорным ограничениям, в конкретной динамической системе. Искомые зависимости параметров алгоритма от количества измеряемых реализаций определяются эмпирически через решение специальной экстремальной задачи, в которой минимизируется отклонение выхода алгоритма от тестовой функции. Для оптимизации времяемкого процесса адаптации алгоритма к системе, предполагающего моделирование большого числа независимых траекторий стохастического уравнения, используется распараллеливание вычислений. Приведен модельный пример, иллюстрирующий предложенные конструкции. Рассмотрена система, упрощенно описывающая популяционную динамику двух взаимодействующих видов. Представлены результаты расчетов и характеристики эффективности распараллеливания.
Ключевые слова: стохастическое дифференциальное уравнение, динамическая реконструкция,
управляемая модель, настройка параметров, распараллеливание вычислений.
ОБРАЗЕЦ ЦИТИРОВАНИЯ
Мельникова Л.А., Розенберг В.Л. Алгоритм динамической реконструкции входов стохастического дифференциального уравнения: настройка параметров и численные
эксперименты // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2019. Т. 8, № 4. С. 15-29. DOI: 10.14529/cmsel90402.
Введение
Проблематика реконструкции неизвестных характеристик управляемых систем в режиме реального времени на основе неточной и/или неполной информации о фазовом состоянии возникает во многих теоретических и практических исследованиях и по этой причине привлекает пристальное внимание специалистов в различных областях. Задачи динамического восстановления входов могут рассматриваться как специальный случай задач теории идентификации; как правило, они являются некорректными (в силу
‘Статья рекомендована к публикации программным комитетом научной конференции «Параллельные вычислительные технологии (ПаВТ) 2019».
неединственности искомого параметра, как следствия неточности измерений, и в силу разрывности обратного оператора) и требуют применения регуляризирующих процедур. В отличие от достаточно большого количества работ, в которых изучаются апостериорные постановки обратных задач, а разрешающие алгоритмы могут использовать всю историю измерения выхода системы, мы ориентируемся на так называемый метод динамического обращения, предложенный и развитый в работах А.В. Кряжимского, Ю.С. Осипова и их коллег [1—3]. Он основан на сочетании принципов теории позиционного управления [4] и идей теории некорректных задач [5]. Задача реконструкции заменяется задачей позиционного управления вспомогательной динамической системой (моделью) с использованием обратных связей. Фактически необходимо подобрать, на основании поступающей информации о фазовом состоянии, такой закон модельного управления, чтобы движение модели отслеживало динамику реальной системы. При выполнении определенных условий оказывается, что в этом случае модельное управление аппроксимирует (в смысле определенной метрики) неизвестное входное воздействие. Метод динамического обращения был реализован для различных систем, описываемых обыкновенными дифференциальными уравнениями (ОДУ) [1, 2], дифференциально-функциональными уравнениями [6],
уравнениями и вариационными неравенствами с распределенными параметрами [3, 7] и др. Были также разработаны алгоритмы восстановления, функционирующие в условиях неполной информации о фазовом состоянии динамической системы [2, 8].
Настоящая работа продолжает исследования (в рамках указанной теории) задач реконструкции для линейного/квазилинейного стохастического дифференциального уравнения (СДУ) [9-11]. Предполагается, что в системе действуют два неизвестных неслучайных возмущения: одно входит в детерминированный член уравнения и
влияет на математическое ожидание искомого процесса, второе входит в интеграл Ито и характеризует амплитуду случайных помех. Рассматриваются условия, при выполнении которых возможно одновременное восстановление возмущений на основе дискретной по времени информации (в том числе неполной) о некотором, достаточно большом, количестве реализаций фазового вектора. Показывается, что задача для СДУ сводится к обратной задаче для системы ОДУ, описывающей динамику математического ожидания и ковариационной матрицы исходного процесса. Разработана конечношаговая программно-ориентированная разрешающая процедура, базирующаяся на конструкциях теории позиционного управления с моделью, адаптированных ранее к обратным задачам в различных постановках для систем ОДУ [1, 2, 8]; получена оценка качества восстановления в зависимости от количества доступных измерению реализаций [11]. Существенную трудность в применении алгоритмов реконструкции указанного типа представляет их адаптация к конкретной динамической системе. Дело в том, что до сих пор не существует универсальной процедуры подбора модельных параметров даже в случае задач для ОДУ; сложности такого рода «всплывали» во многих численных примерах [2]. Новизна данной статьи состоит в рассмотрении этого прикладного аспекта решения задачи реконструкции и разработке эмпирической процедуры для автоматизации процесса калибровки модели. В рассматриваемой задаче настройка параметров алгоритма предполагает многочисленные прогонки вычислительного модуля и моделирование большого числа траекторий СДУ. Для оптимизации этого процесса используется распараллеливание времяемких вычислительных процедур.
Статья организована следующим образом. В разделе 1 приводится формальная постановка задачи динамической реконструкции входов квазилинейной системы СДУ. Раздел 2 посвящен описанию процедуры сведения задачи для системы СДУ к задаче для системы ОДУ. Конструктивный алгоритм решения последней и результаты, характеризующие его сходимость, представлены в разделе 3. Прикладные аспекты, связанные с настройкой параметров алгоритма на определенную систему и с распараллеливанием этого процесса, обсуждаются в разделе 4. В разделе 5 приводится иллюстративный пример с результатами реконструкции возмущений в модельной системе СДУ. В заключении резюмируются полученные результаты и намечаются направления дальнейших исследований.
1. Постановка задачи реконструкции
Простые линеаризованные модели могут быть полезны при изучении динамики рыночных цен при воздействии случайных факторов, изменения численности многовидовой биологической популяции в стохастической среде или процессов хаотического движения однотипных частиц. В этом контексте рассматривается квазилинейная система СДУ специального вида с диффузией, зависящей от фазового состояния:
dx(t, и) = (A(t)x(t, и) + B(t)ui(t) + f(t)) dt + U2(t)x(t, и) d£(t, ш), ж(0, и) = xq. (1)
Здесь t £ Т = [0, я?], х = (xi,X2, • • •, хп) £ Мп — вектор-столбец, £ £ М; жо — известный детерминированный или случайный (нормально распределенный) вектор начальных условий; си £ fl, (fi,F,P) — вероятностное пространство; £(t,oj) — стандартный скалярный винеровский процесс (т. е. выходящий из нуля процесс с нулевым математическим ожиданием и дисперсией, равной t); A(t) = {oy(i)}, B(t) = {hj(t)} и f(t) = {/*(£)} — непрерывные матричные функции размерности п х п, п х г и п х 1 соответственно. На систему действуют два внешних неслучайных возмущения: вектор u\(t) =
(un(t),ui2(t), ■ ■ ■ ,uir(t)) G Мг и диагональная матрица = {«2i(i)j «22(i)5 • • • ,U2n(t)} £
Mnxn. Воздействие щ входит в детерминированную компоненту и влияет на математическое ожидание искомого процесса. Поскольку U2xd£ = (u2iXid^,U22X2d^,... ,U2nxnd^), то можно считать, что вектор U2 = («215 U22,.U2n) регулирует амплитуду случайных помех. Полагаем, что векторы щ и U2 принимают значения из априори известных выпуклых компактов Sui и Su2, щ(-) £ L2(T;Mr), и2(-) £ L2(T]Rn).
Под решением системы (1) понимаем случайный процесс, удовлетворяющий при любом t с вероятностью 1 соответствующему интегральному тождеству, содержащему в правой части стохастический интеграл Ито. Известно, что сделанные предположения обеспечивают существование единственного решения, являющегося нормальным марковским процессом с непрерывными реализациями [12].
Рассматриваемая задача формулируется следующим образом. В дискретные, достаточно частые, моменты времени т* £ Т, т* = id, ó = d/l, i £ [0 : l], поступает информация о некотором количестве N реализаций случайного процесса х(тг). Считаем, что I = l(N) и существуют оценки mf математического ожидания процесса m(t) = Mx(t) и Df ковариационной матрицы D(t) = M(x(t) — m(t))(x(t) — m{t))' (штрих означает транспонирование) такие, что выполняется соотношение
(2)
причем h(N), g(N) —> 0 при N —> оо. Символ || • || обозначает евклидову норму в соответствующем пространстве. Статистические процедуры построения стандартных оценок mf и Df [13] допускают модификации, обеспечивающие выполнение (2).
Задача состоит в построении алгоритма динамической реконструкции неизвестных возмущений щ(-) и дгО), определяющих случайный процесс ж(-), по дискретной информации о его реализациях, причем вероятность сколь угодно малого отклонения приближений от искомых входов в метрике соответственно пространств L2(T;Mr) и L2(T;Mn) должна быть близка к 1 при достаточно большом N и специальным образом согласованном с N шаге временной дискретизации 5 = S(N) = $/Z(JV).
Квазилинейность системы (1) позволяет, в соответствии с методом моментов [14], заменить задачу для СДУ задачей для системы ОДУ, которой удовлетворяют математическое ожидание и ковариационная матрица искомого процесса. После чего возможно организовать процедуру динамической реконструкции возмущений в детерминированном и стохастическом членах правой части в ситуации, когда доступны одновременные измерения достаточно большого количества траекторий (например, движения однотипных частиц).
2. Редукция задачи
Следуя [9—11], заменим задачу восстановления для квазилинейного СДУ задачей для системы ОДУ. Нам понадобится
Лемма. Стандартные оценки mf математического ожидания т(ъ) и Df ковариационной матрицы D(Ti), построенные по N (N > 1) реализациям
х1 (т^, х2(т^,..., xN(т*) случайных величин х(тг) по следующим правилам [13]:
N
(3)
г=1
м
1
1
N
mf = ^7^ЖГ(Тг),
dn =
N
(xr(n) - mf)(xr(n) - m.
N - 1 ^
Г— 1
Ny
i ) 5
(4)
обеспечивают выполнение свойства (2).
Отметим, что оценку (4) с целью оптимизации вычислений (в частности, для их распараллеливания) удобно переписать в виде:
Df =
1 N
хЪУ
Г= 1
ХГ(п)(хГ(7Д)'-
N
N-1'
-т,-
N/N\f
I (
\т\
(5)
Доказательство Леммы приведено в [11], где получены следующие параметры оценки (2):
6(N) = C1/ATin{Q’Q(1/2+3£)}, l(N) = d/5(N),
h(N) = C2/N1/2~e, g(N) = C3/lVmin{1"“’( (6)
Здесь 0<е<1/2,0<а<1 — параметры, позволяющие регулировать порядки малости в аппроксимационной и вероятностной частях оценки (2). Очевидно, S(N) -А 0, а также h(N), g(N) 0, как требуется в соотношении (2), при любых допустимых значениях ей а. Через Ci обозначаем вспомогательные константы, которые зависят от структуры системы и от заданных ограничений на возмущения (но не от самих функций, подлежащих реконструкции) и могут быть выписаны явно.
Кратко опишем процедуру получения системы ОДУ, которой удовлетворяют математическое ожидание и ковариационная матрица исходного процесса (подробное изложение приведено в [10, 11]). Введем обозначения: то = Мхо, Dq = M(xq— то){хо — то)'. В силу квазилинейности исходной системы и равенства нулю математического ожидания интеграла Ито величина m(t) не зависит от u2(t); ее динамика задается уравнением
m(t) = A(t)m(t) + B(t)ui(t) + f(t), m £ Mn, m(0) = mo-
Ковариационная матрица D(t) явно не зависит от rai(t); в соответствии со схемой метода моментов [14], получаем следующее матричное уравнение:
D{t) = A(t)D(t) + D(t)A'(t) + U2{t){D{t) + D £ Rnxn, D{0) = D0.
Его удобно переписать в виде более традиционного для рассматриваемых задач векторного уравнения, размерность которого, с учетом симметричности матрицы D(t), равна nd = (га2 + га)/2. Вводится вектор d(t) = {ds(t)}, s £ [1 : гад], состоящий из последовательно записанных и пронумерованных элементов матрицы D(t), взятых построчно, начиная с элемента, расположенного на главной диагонали; его координаты находятся по элементам матрицы D(t) = {dij(t)}, i,j £ [1 : га]:
ds(t) = dij(t), i < j, s = (n — i/2)(i — 1) + j.
Затем матричное уравнение записывается в виде
d(t) = A(t)d(t) + B(d(t),m(t))u3(t), d(to) = do.
Здесь матрица A(t) : T —$■ Rn<lXnd и диагональная матрица В(d(t), m(t)) : T —> RAdXTld могут быть выписаны явно [11], вектор щ(Ь) = {«3S(i)}, s £ [1 : rid] получается из u2(t)
по формулам щ8 = u2iii2j, г < j, s = (га — г/2)(г — 1) + j, он принимает значения
из некоторого выпуклого компакта SU3 £ Mnd и принадлежит пространству L2(T]Wld). Начальное состояние do получается из Do, а измерения df — из Df.
Таким образом, имеем систему уравнений размерности га + гад:
m{t) = A(t)m(t) + B(t)ui(t) + f(t), rra(0) = m0, (7)
d(t) = A(t)d(t) + B(d(t),m(t))u3(t), d(0) = do. (8)
Переформулируем исходную задачу восстановления. Исходя из (2), считаем, что по ходу развития процесса в дискретные моменты времени п £ Т поступает неточная информация mf, df о фазовом состоянии системы (7), (8) такая, что выполняется соотношение
р( max {\\mf - т(тг)\\ ,\\df - d(Ti)\\} < h(N)) = 1 - g(N), (9)
\ie[l:Z(iV)] /
где h(N), g(N) —> 0 при N —> oo. Следует указать алгоритм динамической реконструкции неизвестных возмущений rai(-) и щ(-) по информации (9), причем вероятность сколь угодно малого отклонения приближений от искомых входов в метрике соответственно пространств L2(T;Rr) и L2(T; Rnd) должна быть близка к 1 при достаточно большом N и специальным образом согласованном с N шаге выполнения измерений ó = S(N) = d/fN). Отметим, что при дополнительных предположениях восстанавливается и функция и2{-) с помощью координат вектора газ(-), равных га^, i £ [1 : га].
Система (7), (8) нелинейна по фазовым переменным, но линейна по управлению. Сформулированная для нее задача соответствует рассмотренной в [1, 2]. В настоящей работе модифицируется конечношаговый разрешающий алгоритм, предложенный ранее для ОДУ, и показывается, что имеет место конструктивное согласование его параметров с количеством доступных измерению реализаций исходного случайного процесса.
3. Алгоритм реконструкции возмущений
Опишем поэтапно искомый алгоритм для системы (7), (8). В начальный момент то = О фиксируется значение N, определяются величины lN = l(N), hN = h(N) и gN = g(N) (см. (6)) и строится равномерное разбиение промежутка Т с шагом 6N = d/lN: Ti £ Т, г* = idN, i £ [0 : lN]. Вводится дискретная модель, динамика и начальное состояние которой определяются соотношениями
wm(Ti+1) = Wm{n) + (+ B(Ti)vii + f{n))SN, wm{0) = m0, (10)
Wd{n+i) = wd(n) + (A(Ti)df + B{df ,mf)vzi)5N, геД0) = d0. (11)
Здесь i £ [0 : lN — 1], wm, wd — модельные переменные, аппроксимирующие фазовое состояние системы (7), (8), гД — управляющие воздействия, вычисляемые в момент г* по правилам, которые конкретизируются ниже.
Работа алгоритма разбивается на lN однотипных шагов. На г-м шаге, который выполняется на интервале (тц Tj+i], исходными данными для вычислений служат оценки mf, df и сформированное к этому моменту состояние модели wd(Ti). Модельные
управления определяются следующим образом: гД и представляют собой единственные решения экстремальных задач
Vu = argmin |2(гст(т*) - mf, В(тДщ) + аf ||щ||2 : щ £ 5Д j , (12)
v3i = argmin {2(лДт;) - df, B{df, mf )v3) + ad |ф3||2 : v3 £ 5из} , (13)
где (•, •) — скалярное произведение в соответствующем евклидовом пространстве, af, аf — параметры регуляризации. После вычисления управлений (12) и (13) пересчитывается согласно (10) и (11) состояние модели wm(Ti+i), wd(Ti+i). Процесс заканчивается в конечный момент времени $.
Пусть 77* = t/*(m(-), d(-)) — множество всех возмущений и(-) = (щ(-),и3(-)) £ L2(Т; Mr+nd), порождающих пару (m(-), d(-)). Минимальный по норме пространства L2 (Т; Mr+nd) элемент множества Н* обозначим через «*(•), его существование и единственность обусловлены выпуклостью 77* и строгой выпуклостью нормы в L2 (Т; Mr+nd), см. [1, 2]. Функция uN(t) = (vf(7), vf(t)), vf (7) = vf, vf(7) = vf, 7 £ [т;,тт), i £ [0 : lN -1] определяется из соотношений (12), (13). Приведем основной результат [11] о сходимости алгоритма.
Теорема. Пусть выполняются условия согласования параметров
WV
aN а?
6N + hN 6N + hN
a
N
a
N
0 при N —> oo.
(14)
Тогда последовательность {иД-)} компактна в L2(Т; Mr+nd) и имеет место сходимость
Р (ДД‘) -«*(•) ||L2(T;ir+"d) "А 0) -А 1 при N—>00.
(15)
В предположении об ограниченности вариации реальных возмущений справедлива следующая оценка точности алгоритма относительно количества реализаций процесса, доступных измерению:
где К i, К 2 и К3 — некоторые положительные константы.
Заметим, что константы в (16), (17) зависят от параметров исходной системы, от априорных структурных ограничений, в частности, на вариацию реальных возмущений, но не зависят от самих восстанавливаемых функций. Выбор показателей степени величины 1/N в аппроксимационной и вероятностной частях оценки типа (16) допускает определенный произвол в диапазоне изменения е и а из соотношений (6). В теореме критерием является равенство показателей, его обеспечивают соотношения (17). Очевидно, возможны другие варианты. В случае, когда вектор щ = («21, «22, ■ ■ ■, Щп) единственен, а его координаты неотрицательны (ограничение представляется естественным для характеристики помехи с нулевым средним, ниже считаем эти условия выполненными), алгоритм (10) (14) восстанавливает и саму функцию U2{-) в смысле метрики пространства 42(Т;МП) с оценкой точности вида (16).
4. Настройка параметров алгоритма
Количество N доступных измерению траекторий исходного СДУ (1) определяет как погрешность hN входных оценок (2), (9), так и параметры алгоритма 6N, а!^ и обеспечивающие сходимость (15), (16). Соотношения не прописаны явно, как это необходимо для расчетов, а задаются в виде порядка сходимости при 1/N —)• Ос точностью до констант, для которых оценки сверху могут быть выписаны явно, но сами значения нуждаются в уточнении. Отметим подтвержденную экспериментально чувствительность алгоритма к ограничениям на функции, формирующие исходную систему (1), на множества допустимых значений возмущений Su\ и Su2, на вариацию возмущений и др. Необходимо получить соотношения (17), обеспечивающие успешную работу алгоритма по восстановлению в конкретной системе вида (1) любых возмущений, удовлетворяющих априорным ограничениям (в этом суть калибровки алгоритма). На данном этапе полагаем, что «оптимальные» зависимости параметров определяются эмпирически через решение следующей экстремальной задачи:
(^ A3*)=argminj £ PN\KM3(-)-uĄ-)\\L2{Ty. (Аь А2, К3) G SK 1 . (18)
I NeiN J
Здесь L2(T) = L2(T; Mr+nd), к (-) — выход алгоритма, зависящий от констант
из (17); IN — множество рассматриваемых значений N, /3N — весовые коэффициенты, Xlive/™ PN = 1; &к — некоторое множество допустимых значений вектора (К 1, К2, К3), как правило, параллелепипед в положительном октанте. Фиксируется некоторая реальная функция «*(•), вводится равномерная сетка по К1, К2, К3 и полным перебором по ней
2/13
(16)
где D\ и D2 — некоторые положительные константы.
Как показано в [П], для получения оценки (16) следует положить
hN = 1/!V6/13, SN = Ai/IV6/13, a% = A2/JV4/13, a$ = A3/1V4/13, (17)
решается задача (18), что требует многочисленных прогонок вычислительного модуля с моделированием большого числа независимых траекторий СДУ (1) для формирования оценок (2) для различных N. Эксперименты подтверждают применимость решения (18) для восстановления других реальных возмущений при сохранении заданных ограничений. Для получения динамики уравнения (1) используется метод Эйлера с заменой винеровского процесса последовательностью случайных импульсов (аппроксимационная схема детально описана в [15]); его среднеквадратический порядок точности для квазилинейного СДУ равен 0(Ss), где — шаг моделирования. Этот шаг должен быть мал относительно шага óN, с которым поступает информация о траекториях и, соответственно, работает процедура восстановления неизвестных функций u\(t) и u2(t). Таким образом, в процедуре подбора параметров алгоритма для конкретной системы задействованы две временные сетки, причем моделирование с более мелким шагом имеет место только на стадии настройки алгоритма, при решении реальной задачи входная информация должна поступать из «внешнего мира». Отметим, что вероятностный характер оценки сходимости (16) предполагает усреднение выхода алгоритма при фиксированном наборе параметров по серии запусков.
Вычислительные эксперименты показали, что описанная процедура, подразумевая множество прогонов на сетке по нескольким параметрам алгоритма, требует значительных ресурсов при расчете на последовательной машине, однако допускает естественное (и эффективное) распараллеливание. Оно базируется на стандартной схеме «мастер-рабочий» [16] с единым загрузочным MPI-модулем, почти равномерным распределением вычислительной нагрузки по всем ядрам и отсутствием обменов между рабочими модулями. Моделирование большого числа независимых траекторий СДУ с мелким временным шагом является единственной времяемкой частью процедуры, поэтому оно распределяется поровну между всеми ядрами. В момент поступления входной информации об оценках (2), (3), (5) все модули проводят расчет одинаковых по размеру порций xr(Ti) и усредняют как сами величины, так и их квадраты, в соответствии с формулами (3), (5). Затем рабочие отправляют свой вклад в оценки мастеру, который добавляет свою порцию и, усредняя полученные суммы по количеству задействованных ядер (легитимность операции обеспечивается одинаковостью порций), окончательно формирует (3) и (5). Все вычисления в соответствии с алгоритмом реконструкции (10)—(13) выполняются мастером. Блок-схема основного фрагмента параллельной процедуры приведена на рис. 1.
5. Модельный пример
Рассмотрим квазилинейную систему, описывающую случайный процесс, близкий к известному средне-возвратному процессу Орнштейна—Уленбека [12]:
dxi(t) = (~2xi(t) + 0.1жг(£) + un(t) + ui2(t) + 1 )dt + U2i(t)xi(t) d£(t),
dx2(t) = (O.lxi(t) - x2(t) + 2ui2(t) +2)dt + u22(t)x2(t) d£(t),
t G T = [0,1], жДО) = 1,x2(0) = 1. (19)
Здесь u(t) = (ui(t), u2(t)), ui(t) = («П(i), «12^)) и u2(t) = («21 (i), «22(i)) — неизвестные возмущения, подлежащие восстановлению, uu(t) G [—2,2], U2i(t) £ [0,2], i = 1,2;
£(i) — стандартный скалярный винеровский процесс. Уравнения типа (19) используются, в частности, в некоторых простейших моделях, описывающих динамику численности
Рис. 1. Схема основного фрагмента параллельной процедуры. Операции, выполняемые только мастером, помечены символом «М», только рабочими — символом «W» (остальные операции выполняются на всех ядрах)
относительно устойчивых популяций биологических особей [12]. В таком случае величина x(t) = (xi(t), Х2(t)) представляет текущую численность (в условных единицах) двух взаимодействующих видов; структура детерминированной части уравнения определяет некоторые численности (например, средние величины за прошедший интервал), стремление вернуться к которым «подсознательно» присутствует у популяций. Такому возврату могут препятствовать соседи и воздействие внешних факторов через возмущение u\(t), влияющее на математическое ожидание процесса, и возмущение u2(t), характеризующее амплитуду случайных шумов. Измерению (в дискретные моменты времени) доступны реализации x(t), соответствующие различным колониям. Запишем для (19) систему ОДУ, соответствующую (7), (8):
mi(t) = -2mi(t) + 0.1m2(t) + иц (t) + u\2(t) + 1, m2(t) = O.lmi(t) - m2(t) + 2ui2(t) + 2, di(t) = -Ad\{t) + 0.2 d2(t) + {di(t) + ml(t))u3i(t), d2(t) = 0.1di(t) - 3d2(t) + 0.1d3(t) + {d2(t) + mi{t)m2{t))u32{t), d3(t) = 0.2d2(t) - 2d3(t) + (d3(t) + ml(t))u33(t),
®i(0) = 1, ®2(0) = 1, di(0) = 0, d2(0) = 0, d3(0) = 0. (20)
Здесь u31 = iĄ-l, u32 = u21u22, u33 = iĄ2.
Целью вычислительных экспериментов являлась адаптация алгоритма к конкретной динамической системе (19)—(20), т.е. получение универсального набора зависимостей его параметров óN, аи а^ от N через К3, К2, К3 (см. (17)) для восстановления любых
возмущений, удовлетворяющих заданным ограничениям. В первой серии для решения (18) были выбраны функции
u\(t) = (sin 10t, 1), U2(t) = (л/t, 10-\/t3/3 — t2/2 + 3t/16)
и параметры Ki G [0.1,5], i = 1,2,3, шаг сетки no Ki равный 0.05, IN = {103,106}, /3N = {0.5,0.5}. Были получены оптимальные значения констант (в смысле (18)) в соотношениях (17) для системы (19) и принятых ограничений на возмущения: К* = 0.6, К% = 0.35, К% = 3.5. Результаты восстановления функций u\(t) и полученные для
найденных Ki и различных N, в принципе согласуются с основным утверждением статьи, подтверждая сходимость вида (15): чем больше N, тем меньше погрешность реконструкции, см. рис. 2, 3. Затем полученные соотношения параметров были протестированы на
Рис. 2. Восстановление возмущений г*ц, «12 (верхний ряд) и «21) «22 (нижний ряд): реальные функции изображены пунктирной линией, результаты их восстановления — сплошной; горизонтальная ось соответствует времени t, вертикальная — значениям возмущений. Параметры: N = 103, SN = 0.025, Ss = SN/102, = 0.04, а^ = 0.4,
погрешность ||«дг(-) — «*(•)!!l2(t) = 0.868
других возмущениях, удовлетворяющих заданным геометрическим ограничениям. На рис. 4 приведен пример восстановления «г(£) = (л/2(1 — 2|t — 0.51), 1 + 0.5sin30t) со значениями параметров алгоритма из варианта, представленного на рис. 3. Точность реконструкции сравнима с результатами, полученными для тестового (в смысле (18)) возмущения. Отметим, что вторая координата возмущения «2 восстанавливается хуже первой, что объясняется большей амплитудой шума во втором уравнении системы (19).
Для иллюстрации показателей, характеризующих качество распараллеливания, был выбран вариант из рис. 3 с симуляцией 128000 траекторий уравнения (19). Моделирование проводилось в Институте математики и механики им. Н.Н.Красовского УрО РАН на гибридном суперкомпьютере «Уран», который состоит из 196 вычислительных узлов, оснащенных процессорами Intel Хеоп и графическими ускорителями NVIDIA Tesla, и имеет пиковую производительность более 215 Тфлопс. Отметим, что в общей сложности
Рис. 3. Восстановление возмущений иц, «12 (верхний ряд) и «21, «22 (нижний ряд): реальные функции изображены пунктирной линией, результаты их восстановления — сплошной; горизонтальная ось соответствует времени t, вертикальная — значениям возмущений. Параметры: N = 106, óN = 0.001, Ss = SN /102, = 0.005, а^ = 0.05,
погрешность Ци^(-) — «*(•)||ц2(т) = 0.118
пользователям доступно 1940 вычислительных ядер CPU, 314 платы GPU и 4 Тб оперативной памяти. В наших экспериментах использовался компактный раздел кластера, называемый Apollo (16 узлов по два 18-и ядерных процессора Intel Xeon CPU Е5-2697 с частотой 2.30 ГГц, оперативная память 256 Гб, локальный жесткий диск 1 Тб); его узлы объединены высокоскоростной сетью Infiniband нового поколения со скоростью передачи данных 100 Г бит/с.
Рис. 4. Восстановление альтернативных возмущений «21 и «22: реальные функции
изображены пунктирной линией, результаты их восстановления — сплошной; горизонтальная ось соответствует времени t, вертикальная — значениям возмущений. Параметры: N = 106, SN = 0.001, Ss = SN/102, ct^ = 0.005, a^ = 0.05, погрешность \\uN(-) -«*(•)||l2(t) = 0.169
Результаты тестирования представлены на рис. 5, где приведены графики зависимостей ускорения и эффективности параллельного алгоритма от количества ядер. Напомним, что ускорение Sp и эффективность Ер вычисляются, соответственно, по формулам: Sp = Т\/Тр, Ер = Sp/p, Тр — время выполнения программы (для одного набора параметров) на
кластере для р ядер, Т\ — время работы последовательного алгоритма. Для «идеального» параллельного алгоритма ускорение равно числу ядер; в этом случае имеем единичную эффективность.
Рис. 5. Ускорение Sp (слева) и эффективность Ер (справа) в зависимости от количества ядер: реальные характеристики изображены сплошной линией, идеальные величины — пунктирной; горизонтальная ось соответствует количеству ядер р, вертикальная — значениям Sp и Ер
Из рис. 5 следует, что эффективность распараллеливания достаточно высока, причем с ростом числа задействованных процессоров она не падает ниже разумного уровня (например, прир = 128 имеем Sp = 103.5 и Ер = 0.81). Это важно, поскольку для настройки алгоритма на конкретную систему и его тестирования (в вероятностном смысле) требуется достаточно большое количество запусков программы.
Заключение
В работе рассмотрена обратная задача для квазилинейной стохастической системы с диффузией, зависящей от фазового состояния. Постановка предполагает динамическую реконструкцию двух неизвестных неслучайных возмущений, входящих в детерминированную и стохастическую части системы, на основе точных дискретных измерений некоторого количества реализаций искомого случайного процесса. Задача для СДУ заменяется обратной задачей для системы двух ОДУ, которым удовлетворяют математическое ожидание и ковариационная матрица исходного процесса. Для конструктивного алгоритма восстановления получена оценка точности относительно количества доступных измерению реализаций. Ключевым результатом статьи следует считать разработку процедуры адаптации алгоритма к конкретной динамической системе и ограничениям на ее структурные элементы; это обусловило необходимость прогонки вычислительного модуля на сетке параметров и моделирования большого числа траекторий СДУ для получения статистических оценок моментов. С целью экономии вычислительных ресурсов введено естественное и достаточно эффективное распараллеливание основных времяемких процедур.
В качестве перспективных направлений исследований отметим создание теоретического обоснования процедуры автоматизации настройки алгоритма через поиск экстремума некоторого критерия качества реконструкции на множестве допустимых параметров и апробацию методики на других алгоритмах восстановления входов динамических систем.
Литература
1. Кряжимский А.В., Осипов Ю.С. О моделировании управления в динамической системе // Изв. АН СССР. Техн. кибернетика. 1983. № 2. С. 51-60.
2. Осипов Ю.С., Кряжимский А.В., Максимов В.И. Методы динамического восстановления входов управляемых систем. Екатеринбург: Изд-во УрО РАН, 2011. 291 с.
3. Максимов В.И. Задачи динамического восстановления входов бесконечномерных систем. Екатеринбург: Изд-во УрО РАН, 2000. 305 с.
4. Красовский Н.Н., Субботин А.И. Позиционные дифференциальные игры. М.: Наука, 1984. 456 с.
5. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука, 1978. 142 с.
6. Близорукова М.С., Максимов В.И. Об одном алгоритме реконструкции траектории и управления в системе с запаздыванием // Труды Института математики и механики УрО РАН. 2012. Т. 18, № 1. С. 109-122. DOI: 10.1134/S0081543813020065.
7. Кадиев А.М., Максимов В.И. О реконструкции управлений в параболическом уравнении // Дифференц. уравнения. 2007. Т. 43, № 11. С. 1545-1552.
8. Кряжимский А.В., Осипов Ю.С. Об устойчивом позиционном восстановлении управления по измерениям части координат / / Некоторые задачи управления и устойчивости. 1989. С. 33-47.
9. Розенберг В.Л. Задача динамического восстановления неизвестной функции в линейном стохастическом дифференциальном уравнении // Автом. и телемех. 2007. № 11. С. 76-87. DOI: 10.1134/S0005117907110069.
10. Розенберг В.Л. Восстановление амплитуды случайной помехи в линейном стохастическом уравнении по измерениям части координат / / Журнал вычисл. математики и мат. физики. 2016. Т. 56, № 3. С. 377-386. DOI: 10.7868/S0044466916030169.
11. Rozenberg V.L. Dynamical Input Reconstruction Problem for a Quasi-linear Stochastic System // IFAC-PapersOnLine. 2018. Vol. 51, no. 32. P. 727-732. DOI: 10.1016/j.ifacol.2018.11.460.
12. Оксендаль Б. Стохастические дифференциальные уравнения. Введение в теорию и приложения. М.: Мир, 2003. 408 с.
13. Королюк В.С., Портенко Н.И., Скороход А.В., Турбин А.Ф. Справочник по теории вероятностей и математической статистике. М.: Наука, 1985. 640 с.
14. Черноусько Ф.Л., Колмановский В.Б. Оптимальное управление при случайных возмущениях. М.: Наука, 1978. 272 с.
15. Милыптейн Г.Н. Численное интегрирование стохастических дифференциальных уравнений. Свердловск: Изд-во УрГУ, 1988. 224 с.
16. Гергель В.П. Теория и практика параллельных вычислений. М.: БИНОМ, 2007. 424 с.
Мельникова Лидия Антоновна, математик 1-й кат., Институт математики и механики
им. Н.Н. Красовского УрО РАН (Екатеринбург, Российская Федерация)
Розенберг Валерий Львович, к.ф.-м.н., с.и.с., Институт математики и механики
им. Н.Н. Красовского УрО РАН (Екатеринбург, Российская Федерация)
DOI: 10.14529/cmsel90402
ALGORITHM OF DYNAMICAL INPUT RECONSTRUCTION FOR A STOCHASTIC DIFFERENTIAL EQUATION: TUNING OF PARAMETERS AND NUMERICAL EXPERIMENTS
© 2019 L.A. Melnikova, V.L. Rozenberg
Krasovskii Institute of Mathematics and Mechanics,
Ural Branch of the Russian Academy of Sciences (S. Kovalevskoi str. 16, Yekaterinburg, 620990 Russia)
E-mail: [email protected], [email protected] Received: 16.04.2019
The problem of reconstructing unknown inputs in a stochastic differential equation is investigated by means of the approach of the theory of dynamic inversion. The statement when the simultaneous reconstruction of disturbances in the deterministic and stochastic terms of the equation is performed from the discrete information on a number of realizations of the stochastic process is considered. The problem is reduced to an inverse problem for ordinary differential equations describing the mathematical expectation and covariance matrix of the process. A software-oriented solving algorithm based on constructions of the theory of positional control with a model is designed; an estimate for its convergence rate with respect to the number of measurable realizations is obtained. A program procedure for the automatic tuning of the algorithm’s parameters in order to have the best approximation results for different disturbances satisfying a priori constraints in a specific dynamical system is proposed. Desired dependencies of the algorithm’s parameters on the number of measured realizations are determined empirically via solving a specific extremal problem, where the deviation of the algorithm’s output from some test function is minimized. To optimize the time-taking adaptation process assuming the simulation of a large number of independent trajectories of the stochastic process, the parallelization of calculations is applied. A model example illustrating the method proposed is given. A system approximately describing the population dynamics of two interacting biological species is considered. Calculation results and parallelization efficiency characteristics are presented.
Keywords: stochastic differential equation, dynamical reconstruction, controlled model, parameter tuning, parallelization of calculations.
FOR CITATION
Melnikova L.A., Rozenberg V.L. Algorithm of Dynamical Input Reconstruction for a Stochastic Differential Equation: Tuning of Parameters and Numerical Experiments. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering. 2019. vol. 8, no. 4. pp. 15-29. (in Russian) DOI: 10.14529/cmsel90402.
This paper is distributed under the terms of the Creative Commons Attribution-Non Commercial 3.0 License which permits non-commercial use, reproduction and distribution of the work without further permission provided the original work is properly cited.
References
1. Kryazhimskii A.V., Osipov Yu.S. Modelling of a Control in a Dynamic System. Engrg. Cybernetics. 1984. vol. 21, no. 2, pp. 38-47.
2. Osipov Yu.S., Kryazhimskii A.V., Maksimov V.I. Dynamic Recovery Methods for Inputs of Control Systems. Yekaterinburg, Publishing of the UB of RAS, 2011. 291 p. (in Russian)
3. Maksimov V.L Dynamical Inverse Problems of Distributed Systems. VSP, Utrecht-Boston, 2002. 269 p.
4. Krasovskii N.N., Subbotin А.I. Game-Theoretical Control Problems. Springer, New York, 1988. 517 p.
5. Tikhonov A.N., Arsenin V.Ya. Solutions of Ill-Posed Problems. Wiley, New York, 1981. 258 p.
6. Blizorukova M.S., Maksimov V.I. On a Reconstruction Algorithm for the Trajectory and Control in a Delay System. Proc. Steklov Inst. Math. 2013. vol. 280, no. 1, pp. S66-S79. DOI: 10.1134/S0081543813020065.
7. Kadiyev A. M., Maksimov V.I. On the Reconstruction of Controls in a Parabolic Equation. Differential Equations. 2007. vol. 43, no. 11, pp. 1585-1593. DOI: 10.1134/S0012266107110134.
8. Kryazhimskii A.V., Osipov Yu.S. On a Stable Positional Recovery of Control from Measurements of a Part of Coordinates. Some Problems of Control and Stability. Sverdlovsk, Publishing of the UB of Academy of Sciences of the USSR, 1989. pp. 33-47. (in Russian)
9. Rozenberg V.L. Dynamic Restoration of the Unknown Function in the Linear Stochastic Differential Equation. Autom. Remote Control. 2007. vol. 68, no. 11, pp. 1959-1969. DOI: 10.1134/S0005117907110069.
10. Rozenberg V.L. Reconstruction of Random-Disturbance Amplitude in Linear Stochastic Equations from Measurements of Some of the Coordinates. Comput. Math. Math. Phys. 2016. vol. 56, no. 3, pp. 367-375. DOI: 10.1134/S0965542516030143.
11. Rozenberg V.L. Dynamical Input Reconstruction Problem for a Quasi-linear Stochastic System. IFAC-PapersOnLine. 2018. vol. 51, no. 32. pp. 727-732. DOI: 10.1016/j.ifacol.2018.11.460.
12. Oksendal B. Stochastic Differential Equations: an Introduction with Applications. Springer, Berlin, 1985. 408 p.
13. Korolyuk V.S., Portenko N.I., Skorokhod A.V., Turbin A.F. Handbook on Probability Theory and Mathematical Statistics. Moscow, Nauka, 1985. 640 p. (in Russian)
14. Chernous’ko F.L., Kolmanovskii V.B. Optimal Control under Random Perturbation. Moscow, Nauka, 1978. 272 p. (in Russian)
15. Milshtein G.N. Numerical Integration of Stochastic Differential Equations. Sverdlovsk, Publishing of the Ural State University, 1988. 224 p. (in Russian)
16. Gergel’ V.P. Theory and Practice of Parallel Computing. Moscow, Binom, 2007. 424 p. (in Russian)
