\documentclass{cmi}
\usepackage{upgreek}
\usepackage{cite}

\begin{document}


\classify{УДК 517.977, 519.688}

\renewcommand*{\thefootnote}{\fnsymbol{footnote}} %to make a number in end of title delete this line


\title{АЛГОРИТМ ДИНАМИЧЕСКОЙ РЕКОНСТРУКЦИИ ВХОДОВ
СТОХАСТИЧЕСКОГО ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ:
НАСТРОЙКА ПАРАМЕТРОВ И ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ\footnote{Статья рекомендована к публикации программным
комитетом научной конференции <<Параллельные вычислительные технологии (ПаВТ) 2019>>.}%
}

\author{\copyright~2019 Л.А.~Мельникова, В.Л.~Розенберг}

%\author{User1 \\ \href{mailto:me@somewhere.com}{me@somewhere.com}
%	\and User2 \\ email %\href{mailto:someone@somewhere.com}{someone@somewhere.com} }


\address{Институт математики и механики им. Н.Н.~Красовского УрО РАН\\
	(620990 Екатеринбург, ул.~С.Ковалевской, д.~16)}

\email{E-mail: \href{mailto:meln@imm.uran.ru}{meln@imm.uran.ru},
\href{mailto:rozen@imm.uran.ru}{rozen@imm.uran.ru}}


\received{Поступила в редакцию: 16.04.2019}

\maketitle{}


\renewcommand*{\thefootnote}{\arabic{footnote}}
\setcounter{footnote}{0}
\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}

% оформление аннотации
\begin{abstract}%
Задача реконструкции неизвестных входов
стохастического дифференциального уравнения
исследуется с позиций подхода теории динамического обращения.
Рассматривается постановка, в которой одновременное восстановление
возмущений в детерминированном и стохастическом членах уравнения
проводится на основе дискретной информации
о некотором количестве реализаций случайного процесса.
Задача сводится к обратной задаче для системы обыкновенных дифференциальных уравнений, которым удовлетворяют
математическое ожидание и ковариационная матрица исходного процесса.
Разработан программно-ориентированный алгоритм решения, основанный
на конструкциях теории позиционного управления с моделью; получена оценка его
точности относительно количества доступных измерению реализаций.
Предложена программная процедура настройки параметров алгоритма для получения наилучшего
результата аппроксимации различных возмущений, удовлетворяющих априорным ограничениям,
в конкретной динамической системе.
Искомые зависимости параметров алгоритма от количества измеряемых реализаций определяются
эмпирически через решение специальной экстремальной задачи, в которой минимизируется
отклонение выхода алгоритма от тестовой функции.
Для оптимизации времяемкого процесса адаптации алгоритма к системе,
предполагающего моделирование большого числа независимых траекторий стохастического уравнения,
используется распараллеливание вычислений.
Приведен модельный пример, иллюстрирующий предложенные конструкции.
Рассмотрена система, упрощенно описывающая популяционную динамику двух
взаимодействующих видов. Представлены результаты расчетов и характеристики
эффективности распараллеливания.

\keywords{стохастическое дифференциальное уравнение, динамическая реконструкция,
управляемая модель, настройка параметров, распараллеливание вычислений.}
\end{abstract}
\flushleft{\textbf{ОБРАЗЕЦ ЦИТИРОВАНИЯ}}
\justifying
\begin{citationplace}
Мельникова~Л.А., Розенберг~В.Л. Алгоритм динамической реконструкции входов стохастического дифференциального уравнения:
настройка параметров и численные эксперименты // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2019. Т.~X, №~Y. С.~Z1--Z2.
DOI:~10.14529/cmseXXXXXX.
\end{citationplace}

%Верхние колонтитулы
\markboth{Л.А.~Мельникова, В.Л.~Розенберг}{Алгоритм динамической реконструкции входов\ldots}
% Заголовки разделов формируются при помощи команд  \section{}, \subsection{}, \subsubsection{}

\section*{Введение}
Проблематика реконструкции неизвестных характеристик управляемых систем
в режиме реального времени на основе неточной и/или неполной информации о фазовом состоянии
возникает во многих теоретических и практических исследованиях и по этой причине
привлекает пристальное внимание специалистов в различных областях.
Задачи динамического восстановления входов могут рассматриваться как специальный случай
задач теории идентификации; как правило, они являются некорректными (в силу неединственности искомого
параметра, как следствия неточности измерений, и в силу разрывности обратного оператора)
и требуют применения регуляризирующих процедур.
В отличие от достаточно большого количества работ, в которых изучаются апостериорные постановки обратных задач,
а разрешающие алгоритмы могут использовать всю историю измерения выхода системы,
мы ориентируемся на так называемый метод динамического обращения,
предложенный и развитый  в работах А.В.~Кряжимского, Ю.С.~Осипова
и их коллег \cite{kryazh,oskr2011,maxim2000}.
Он основан на сочетании принципов теории позиционного управления
\cite{kras} и идей теории некорректных задач \cite{tikh}.
Задача реконструкции заменяется задачей позиционного управления
вспомогательной динамической системой (моделью) с использованием обратных связей.
Фактически необходимо подобрать, на основании поступающей информации о фазовом состоянии,
такой закон модельного управления, чтобы движение модели отслеживало
динамику реальной системы. При выполнении определенных условий оказывается,
что в этом случае модельное управление аппроксимирует (в смысле определенной метрики)
неизвестное входное воздействие.
Метод динамического обращения был реализован для различных
систем, описываемых  обыкновенными дифференциальными
уравнениями (ОДУ) \cite{kryazh,oskr2011}, дифференциально-функциональными уравнениями \cite{maxbliz},
уравнениями и вариационными неравенствами с распределенными параметрами \cite{maxim2000, maxim2007} и др.
Были также разработаны алгоритмы восстановления, функционирующие в условиях неполной информации
о фазовом состоянии динамической системы \cite{oskr2011,oskr1989}.

Настоящая работа продолжает исследования (в рамках указанной теории)
задач реконструкции для линейного/квазилинейного стохастического дифференциального уравнения (СДУ) \cite{rozen2007, rozen2016, rozen2018}.
Предполагается, что в системе действуют два неизвестных неслучайных возмущения: одно входит в детерминированный член уравнения
и влияет на математическое ожидание искомого процесса, второе входит в интеграл Ито и характеризует амплитуду случайных помех.
Рассматриваются условия, при выполнении которых возможно одновременное восстановление
возмущений на основе дискретной по времени информации (в том числе неполной)
о некотором, достаточно большом, количестве реализаций фазового вектора.
Показывается, что задача для СДУ сводится к обратной задаче для системы ОДУ, описывающей динамику
математического ожидания и ковариационной матрицы исходного процесса.
Разработана конечношаговая программно-ориентированная разрешающая процедура, базирующаяся
на конструкциях теории позиционного управления с моделью, адаптированных ранее к обратным задачам
в различных постановках для систем ОДУ \cite{kryazh,oskr2011,oskr1989};
получена оценка качества восстановления в зависимости от количества доступных измерению реализаций \cite{rozen2018en}.
Существенную трудность в применении алгоритмов реконструкции указанного типа
представляет их адаптация к конкретной динамической системе.
Дело в том, что до сих пор не существует универсальной процедуры подбора модельных параметров
даже в случае задач для ОДУ; сложности такого рода <<всплывали>> во многих численных примерах \cite{oskr2011}.
Новизна данной статьи состоит в рассмотрении этого прикладного аспекта
решения задачи реконструкции и разработке эмпирической процедуры для автоматизации процесса калибровки модели.
В рассматриваемой задаче настройка параметров алгоритма предполагает многочисленные прогонки
вычислительного модуля и моделирование большого числа траекторий СДУ.
Для оптимизации этого процесса используется распараллеливание времяемких вычислительных процедур.

Статья организована следующим образом. В разделе~1 приводится формальная постановка
задачи динамической реконструкции входов квазилинейной системы СДУ.
Раздел~2 посвящен описанию процедуры сведения задачи для системы СДУ
к задаче для системы ОДУ. Конструктивный алгоритм решения последней и
результаты, характеризующие его сходимость, представлены в разделе~3.
Прикладные аспекты, связанные с настройкой параметров алгоритма на определенную
систему и с распараллеливанием этого процесса, обсуждаются в разделе~4.
В разделе~5 приводится иллюстративный пример с результатами реконструкции
возмущений в модельной системе СДУ. В заключении резюмируются полученные результаты
и намечаются направления дальнейших исследований.

\setcounter{equation}{0}
\section{Постановка задачи реконструкции}
Простые линеаризованные модели могут быть полезны при изучении
динамики рыночных цен при воздействии случайных факторов,
изменения численности многовидовой биологической популяции в стохастической среде или
процессов хаотического движения однотипных частиц. В этом контексте
рассматривается квазилинейная система СДУ специального вида с диффузией, зависящей от фазового состояния:
\begin{equation}\label{2_1}
dx(t,\omega )=(A(t)x(t,\omega) + B(t)u_1(t) + f(t))\,dt + U_2(t)x(t,\omega)\, d\xi (t,\omega), \quad x(0,\omega)=x_{0}.
\end{equation}

Здесь $t\in T=[0,\vartheta]$,\
$x = (x_1, x_2, \ldots, x_n) \in \mathbb{R}^{n}$~--- вектор-столбец,\
$\xi \in \mathbb{R}$;\
$x_{0}$ --- известный детерминированный или случайный (нормально распределенный) вектор начальных условий;\
$\omega \in \Omega$,\ $(\Omega, F, P)$ --- вероятностное пространство;\
$\xi(t,\omega )$ --- стандартный скалярный винеровский процесс
(т.\,е. выходящий из нуля процесс с нулевым математическим ожиданием и
дисперсией, равной $t$);\
$A(t) = \{a_{ij}(t)\}$, $B(t)=\{b_{ij}(t)\}$ и $f(t)=\{f_{i}(t)\}$~---  непрерывные матричные функции размерности
$n\times n$, $n\times r$ и $n\times 1$ соответственно.
На систему действуют два внешних неслучайных возмущения:
вектор $u_1(t) = (u_{11}(t), u_{12}(t), \ldots, u_{1r}(t)) \in \mathbb{R}^{r}$ и
диагональная матрица $U_2(t) = \{u_{21}(t), u_{22}(t), \ldots, u_{2n}(t) \} \in \mathbb{R}^{n\times n}$.
Воздействие $u_1$ входит в детерминированную компоненту и влияет на математическое ожидание искомого процесса.
Поскольку $U_2 x d\xi  = (u_{21} x_1 d\xi, u_{22} x_2 d\xi, \ldots, u_{2n} x_n d\xi)$, то
можно считать, что вектор $u_2 = (u_{21}, u_{22}, \ldots, u_{2n})$ регулирует амплитуду случайных помех.
Полагаем, что векторы $u_1$ и $u_2$ принимают значения из априори известных выпуклых компактов $S_{u1}$ и $S_{u2}$,
$u_1(\cdot) \in L_{2}(T;\mathbb{R}^{r})$, $u_2(\cdot) \in L_{2}(T;\mathbb{R}^{n})$.

Под решением системы (\ref{2_1}) понимаем случайный процесс, удовлетворяющий при любом $t$ с вероятностью 1
соответствующему интегральному тождеству, содержащему в правой части стохастический интеграл Ито.
Известно, что сделанные предположения обеспечивают существование единственного решения,
являющегося нормальным марковским процессом с непрерывными реализациями \cite{oksen}.

Рассматриваемая задача формулируется следующим образом. В дискретные, достаточно
частые, моменты времени $\tau _{i} \in T$, $\tau _{i} =i\delta$, $\delta =\vartheta /l$, $i \in [0:l],$
поступает информация о некотором количестве $N$ реализаций случайного
процесса $x(\tau _i)$.
Считаем, что $l=l(N)$ и существуют оценки $m_{i}^{N}$
математического ожидания процесса $m(t)=Mx(t)$ и $D_{i}^{N}$ ковариационной матрицы
$D(t)=M(x(t)-m(t))(x(t)-m(t))'$ (штрих означает транспонирование) такие,
что выполняется соотношение
\begin{equation}\label{2_2}
P\Bigl(\,\max \limits_{i\in [1:l(N)]} \left\{ \left\| m_{i}^{N} - m(\tau _{i})\right\|,
\left\| D_{i}^{N} - D(\tau _{i}) \right\| \right\} \leq h(N)\Bigr) = 1-g(N),
\end{equation}
причем
$h(N),\,g(N) \rightarrow 0$ при $N\rightarrow \infty$.
Символ  $\| \cdot \|$ обозначает евклидову норму в соответствующем пространстве.
Статистические процедуры построения стандартных оценок $m_{i}^{N}$ и $D_{i}^{N}$ \cite{king}
допускают модификации, обеспечивающие выполнение (\ref{2_2}).

Задача состоит в построении алгоритма динамической реконструкции  неизвестных возмущений
$u_1(\cdot)$ и $u_2(\cdot)$, определяющих случайный процесс $x(\cdot)$, по дискретной информации о его реализациях,
причем вероятность сколь угодно малого отклонения приближений от искомых входов
в метрике соответственно пространств $L_{2}(T;\mathbb{R}^{r})$ и $L_{2}(T;\mathbb{R}^{n})$  должна быть близка к~1
при достаточно большом $N$ и специальным образом согласованном с $N$
шаге временной дискретизации $\delta=\delta(N)=\vartheta /l(N)$.

Квазилинейность системы (\ref{2_1}) позволяет, в соответствии с методом моментов \cite{chern},
заменить задачу для СДУ задачей для системы ОДУ, которой удовлетворяют математическое ожидание и
ковариационная матрица искомого процесса.
После чего возможно  организовать процедуру динамической реконструкции
возмущений в детерминированном и стохастическом членах правой части в ситуации,
когда доступны одновременные измерения
достаточно большого количества траекторий (например, движения однотипных частиц).

\setcounter{equation}{0}
\section{Редукция задачи}
Следуя \cite{rozen2007, rozen2016, rozen2018}, заменим задачу восстановления для квазилинейного СДУ
задачей для системы ОДУ. Нам понадобится

\textbf{Лемма.}
\textit{ Стандартные оценки $m_{i}^{N}$ математического ожидания $m(\tau _i)$ и $D_{i}^{N}$ ковариационной матрицы $D(\tau _i)$,
построенные по $N$\ $(N>1)$ реализациям $x^{1}(\tau_{i}),x^{2}(\tau_{i}),\ldots,x^{N}(\tau_{i})$
случайных величин $x(\tau_{i})$ по следующим правилам {\mbox{\em \cite{king}:}}
\begin{equation}\label{estim1}
m_{i}^{N} = \frac{1}{N} \sum\limits_{r=1}^{N}x^{r}(\tau_{i}),
\end{equation}
\begin{equation}\label{estim2}
D_{i}^{N} = \frac{1}{N-1} \sum\limits_{r=1}^{N}(x^{r}(\tau_{i})-m_{i}^{N})(x^{r}(\tau_{i})-m_{i}^{N})',
\end{equation}
обеспечивают выполнение свойства $(\ref{2_2})$.}

Отметим, что оценку (\ref{estim2}) с целью оптимизации вычислений (в частности, для их распараллеливания)
удобно переписать в виде:
\begin{equation}\label{estim3}
D_{i}^{N} = \frac{1}{N-1} \sum\limits_{r=1}^{N}x^{r}(\tau_{i})(x^{r}(\tau_{i}))' -
\frac{N}{N-1}m_{i}^{N}(m_{i}^{N})'.
\end{equation}

Доказательство Леммы приведено в \cite{rozen2018}, где получены следующие
параметры оценки (\ref{2_2}):
$$
\delta(N) = C_{1}/N^{\min \{\alpha, \alpha(1/2+3\epsilon) \} }, \quad
l(N) = \vartheta / \delta(N), \quad
$$
\begin{equation}\label{param}
h(N) = C_{2}/N^{1/2 - \epsilon }, \quad
g(N) = C_{3}/N^{\min \{1-\alpha, (1-\alpha)(1/2+3\epsilon) \} }.
\end{equation}
Здесь $0<\epsilon < 1/2$, $0<\alpha <1$ --- параметры, позволяющие регулировать порядки малости в аппроксимационной и
вероятностной частях оценки (\ref{2_2}). Очевидно, $\delta(N) \rightarrow 0$, а также $h(N),\,g(N) \rightarrow 0$, как требуется
в соотношении (\ref{2_2}), при любых допустимых значениях $\epsilon$ и $\alpha$.
Через $C_i$ обозначаем вспомогательные константы, которые зависят от структуры системы и от заданных ограничений на возмущения
(но не от самих функций, подлежащих реконструкции) и могут быть выписаны явно.

Кратко опишем процедуру получения системы ОДУ, которой удовлетворяют
математическое ожидание и ковариационная матрица исходного процесса
(подробное изложение приведено в \cite{rozen2016, rozen2018}).
Введем обозначения:
$m_{0}=Mx_{0}$, $D_{0}=M(x_{0}-m_{0})(x_{0}-m_{0})'$.
В силу квазилинейности исходной системы и равенства нулю математического
ожидания интеграла Ито величина $m(t)$ не зависит от $u_2(t)$;
ее динамика задается уравнением
$$
\dot{m}(t)=A(t)m(t)+B(t)u_1(t)+f(t),\,  m \in \mathbb{R}^{n}, \, m(0)=m_{0}.
$$
Ковариационная матрица $D(t)$ явно не зависит от $u_1(t)$; в соответствии со схемой
метода моментов \cite{chern}, получаем следующее матричное уравнение:
$$
\dot{D}(t)= A(t)D(t) + D(t)A'(t) + U_2(t)(D(t) + m(t)m'(t))U'_2(t), \, D \in \mathbb{R}^{n\times n}, \, D(0)=D_{0}.
$$
Его удобно переписать в виде более традиционного для рассматриваемых задач
векторного уравнения, размерность которого, с учетом симметричности матрицы $D(t)$, равна ${n_{d}=(n^2+n)/2}$.
Вводится вектор $d(t)=\{d_{s} (t)\},\,s\in [1:n_{d}]$, состоящий из последовательно записанных и пронумерованных
элементов матрицы $D(t)$, взятых построчно, начиная с элемента, расположенного на  главной диагонали;
его координаты находятся по элементам матрицы
$D(t)=\{d_{ij} (t) \},\,i, j\in [1:n]$:
$$
d_{s} (t)=d_{ij} (t),\quad i\leq j, \quad s=(n-i/2)(i-1)+j.
$$
Затем матричное уравнение записывается в виде
$$
\dot{d}(t)=\bar{A}(t)d(t)+\bar{B}(d(t),m(t))u_3(t),\, d(t_0)=d_{0}.
$$
Здесь матрица ${\bar{A}(t):\,T\rightarrow R^{n_{d}\times n_{d}}}$ и
диагональная матрица ${\bar{B}(d(t),m(t)):\,T\rightarrow R^{n_{d}\times n_{d}}}$
могут быть выписаны явно \cite{rozen2018},
вектор $u_3(t)=\{u_{3s} (t)\}$, $s\in [1:n_{d}]$ получается из $u_2(t)$ по формулам
$u_{3s} = u_{2i} u_{2j}$, $i\leq j$, $s=(n-i/2)(i-1)+j$,
он принимает значения из некоторого выпуклого компакта $S_{u3} \in \mathbb{R}^{n_d}$
и принадлежит пространству $L_{2}(T;\mathbb{R}^{n_d})$.
Начальное состояние $d_{0}$ получается из $D_{0}$, а измерения $d_{i}^{N}$ --- из $D_{i}^{N}$.

Таким образом, имеем систему уравнений размерности $n+n_d$:
\begin{equation}\label{3_1}
\dot{m}(t)=A(t)m(t)+B(t)u_1(t)+f(t), \quad m(0)=m_{0},
\end{equation}
\begin{equation}\label{3_2}
\dot{d}(t)=\bar{A}(t)d(t)+\bar{B}(d(t),m(t))u_3(t), \quad d(0)=d_{0}.
\end{equation}
Переформулируем исходную задачу восстановления. Исходя из (\ref{2_2}), считаем, что
по ходу развития процесса в дискретные моменты времени
$\tau _{i} \in T$ поступает неточная информация $m_{i}^{N}, \, d_{i}^{N}$ о фазовом состоянии
системы (\ref{3_1}), (\ref{3_2}) такая, что выполняется соотношение
\begin{equation}\label{3_3}
P\Bigl(\,\max \limits_{i\in [1:l(N)]} \left\{ \left\| m_{i}^{N} - m(\tau _{i})\right\|,
\left\| d_{i}^{N} - d(\tau _{i}) \right\| \right\} \leq h(N)\Bigr) = 1-g(N),
\end{equation}
где $h(N),\,g(N) \rightarrow 0$ при $N\rightarrow \infty$.
Следует указать алгоритм динамической реконструкции неизвестных возмущений
$u_1(\cdot)$ и $u_3(\cdot)$ по информации (\ref{3_3}),
причем вероятность сколь угодно малого отклонения приближений от искомых входов
в метрике соответственно пространств $L_{2}(T;\mathbb{R}^{r})$ и $L_{2}(T;\mathbb{R}^{n_d})$
должна быть близка к 1 при достаточно большом $N$ и специальным образом согласованном с $N$
шаге выполнения измерений $\delta=\delta(N)=\vartheta /l(N)$.
Отметим, что при дополнительных предположениях восстанавливается и функция $u_2(\cdot)$
с помощью координат вектора $u_3(\cdot)$, равных $u^2_{2i}$, $i\in [1:n]$.

Система (\ref{3_1}), (\ref{3_2}) нелинейна по фазовым переменным, но линейна по управлению.
Сформулированная для нее  задача соответствует рассмотренной в \cite{kryazh,oskr2011}.
В настоящей работе модифицируется конечношаговый разрешающий алгоритм,
предложенный ранее для ОДУ, и показывается, что имеет место
конструктивное согласование его параметров с количеством доступных измерению реализаций
исходного случайного процесса.

\setcounter{equation}{0}
\section{Алгоритм реконструкции возмущений}
Опишем поэтапно искомый алгоритм для системы~(\ref{3_1}), (\ref{3_2}).
В начальный момент $\tau _{0} =0$ фиксируется значение $N$, определяются величины ${l^{N} =l(N)}$,
${h^{N} =h(N)}$ и $g^{N} =g(N)$ (см. (\ref{param})) и строится равномерное разбиение промежутка $T$ с
шагом $\delta^{N} =\vartheta /l^N$:
$\tau_{i} \in T, \, \tau_{i}=i\delta^{N}, \, i\in [0:l^N]$.
Вводится дискретная модель, динамика и начальное состояние которой определяются соотношениями
\begin{equation}\label{4_1}
w_m(\tau_{i+1}) = w_m(\tau_{i}) + (A(\tau_{i})m_{i}^{N} + B(\tau_{i})v_{1i}^{N} + f(\tau_{i}))\delta^{N}, \, w_m(0)=m_{0},
\end{equation}
\begin{equation}\label{4_2}
w_d(\tau_{i+1}) = w_d(\tau_{i}) + (\bar{A}(\tau_{i})d_{i}^{N} + \bar{B}(d_{i}^{N},m_{i}^{N})v_{3i}^{N})\delta^{N}, \, w_d(0)=d_{0}.
\end{equation}
Здесь $i\in [0:l^N-1]$, $w_m$, $w_d$ --- модельные переменные, аппроксимирующие фазовое состояние системы ~(\ref{3_1}), (\ref{3_2}),
$v_{1i}^N$, $v_{3i}^N$ --- управляющие воздействия, вычисляемые в момент $\tau_{i}$ по правилам, которые конкретизируются ниже.

Работа алгоритма разбивается на $l^N$ однотипных шагов. На $i$-м шаге,
который выполняется на интервале $(\tau_{i}, \, \tau_{i+1}]$,
исходными данными для вычислений служат оценки $m_{i}^{N}$, $d_{i}^{N}$
и сформированное к этому моменту состояние модели $w_m(\tau_{i})$, $w_d(\tau_{i})$.
Модельные управления определяются следующим образом:
$v_{1i}^{N}$ и $v_{3i}^{N}$ представляют собой  единственные решения экстремальных задач
\begin{equation}\label{4_3}
v_{1i}^{N} = \arg \min \left\{ 2 \langle w_{m}(\tau_{i})-m_{i}^{N}, \, B(\tau_{i})v_1  \rangle + \alpha_{m}^{N}
\left\| v_1\right\|^{2}:\, v_1 \in S_{u1} \right\},
\end{equation}
\begin{equation}\label{4_4}
v_{3i}^{N} = \arg \min \left\{ 2 \langle w_{d}(\tau_{i})-d_{i}^{N}, \, \bar{B}(d_{i}^{N}, m_{i}^{N})v_3 \rangle + \alpha_{d}^{N}
\left\| v_3\right\|^{2}:\, v_3 \in S_{u3} \right\},
\end{equation}
где $\langle \cdot,\,\cdot \rangle$ --- скалярное произведение в соответствующем евклидовом пространстве,
$\alpha_{m}^{N}$, $\alpha_{d}^{N}$ --- параметры регуляризации.
После вычисления управлений (\ref{4_3}) и (\ref{4_4}) пересчитывается согласно (\ref{4_1}) и (\ref{4_2}) состояние модели $w_m(\tau_{i+1})$, $w_d(\tau_{i+1})$.
Процесс заканчивается в конечный момент времени $\vartheta$.

Пусть $U_{*}=U_{*}(m(\cdot), d(\cdot))$ --- множество всех возмущений
$u(\cdot) = (u_1(\cdot), u_3(\cdot)) \in L_2(T;\mathbb{R}^{r+n_d})$, порождающих пару $(m(\cdot),d(\cdot))$.
Минимальный по норме пространства $L_2(T;\mathbb{R}^{r+n_d})$ элемент множества $U_{*}$ обозначим через $u_{*}(\cdot)$,
его существование и единственность обусловлены выпуклостью $U_{*}$ и строгой выпуклостью нормы в
$L_2(T;\mathbb{R}^{r+n_d})$, см. ~\cite{kryazh, oskr2011}.
Функция $u^{N}(t) = (v_{1}^{N}(t), v_{3}^{N}(t))$, $v_{1}^{N}(t)=v_{1i}^{N}$, $v_{3}^{N}(t)=v_{3i}^{N}$, $t\in [\tau _{i}, \tau _{i+1})$, $i\in [0:l^N-1]$
определяется из соотношений (\ref{4_3}), (\ref{4_4}). Приведем основной результат \cite{rozen2018}
о сходимости алгоритма.

\textbf{Теорема.}
\textit{Пусть выполняются условия согласования параметров
\begin{equation}\label{4_5}
h^{N}, \, \delta^{N}, \, \alpha_{m}^{N}, \, \alpha_{d}^{N} \rightarrow 0, \,
\frac{\delta^{N}+h^{N}}{\alpha_{m}^{N}}, \, \frac{\delta^{N}+h^{N}}{\alpha_{d}^{N}} \rightarrow 0 \quad \mbox{при} \,\, N\rightarrow \infty.
\end{equation}
Тогда последовательность $\{u^{N} (\cdot )\}$ компактна в
$L_2(T;\mathbb{R}^{r+n_d})$ и имеет место сходимость
\begin{equation}\label{4_6}
P\left( \| u^{N}(\cdot ) - u_{*}(\cdot) \|_{L_2(T;\mathbb{R}^{r+n_d})} \, \rightarrow \,0\right) \rightarrow 1 \quad
\mbox{при} \quad N\rightarrow \infty.
\end{equation}
В предположении об ограниченности вариации реальных возмущений справедлива следующая оценка точности алгоритма относительно
количества реализаций процесса, доступных измерению:
\begin{equation}\label{4_7}
P\left(\| u^{N}(\cdot ) - u_{*}(\cdot) \|_{L_2(T;\mathbb{R}^{r+n_d})} \, \leq D_{1}\left( \frac{1}{N} \right) ^{2/13} \right)
= 1 - D_{2}\left( \frac{1}{N} \right) ^{2/13},
\end{equation}
где $D_1$ и $D_2$--- некоторые положительные константы.}

Как показано в \cite{rozen2018}, для получения оценки (\ref{4_7}) следует положить
\begin{equation}\label{4_8}
h^{N} = 1/N^{6/13}, \quad \delta^{N} = K_{1}/N^{6/13}, \quad \alpha_{m}^{N} = K_{2}/N^{4/13}, \quad \alpha_{d}^{N} = K_{3}/N^{4/13},
\end{equation}
где $K_1$, $K_2$ и $K_3$--- некоторые положительные константы.

Заметим, что константы в (\ref{4_7}), (\ref{4_8}) зависят от параметров исходной системы, от априорных структурных ограничений,
в частности, на вариацию реальных возмущений, но не зависят от самих восстанавливаемых функций.
Выбор показателей степени величины $1/N$ в аппроксимационной и вероятностной частях оценки типа (\ref{4_7}) допускает определенный произвол
в диапазоне изменения $\epsilon$ и $\alpha$ из соотношений (\ref{param}).
В теореме критерием является равенство показателей, его обеспечивают соотношения (\ref{4_8}).
Очевидно, возможны другие варианты.
В случае, когда вектор $u_2 = (u_{21}, u_{22}, \ldots, u_{2n})$ единственен, а его координаты неотрицательны
(ограничение представляется естественным для характеристики помехи с нулевым средним, ниже считаем эти условия выполненными),
алгоритм $(\ref{4_1})$--$(\ref{4_5})$
восстанавливает и саму функцию $u_2(\cdot)$ в смысле метрики пространства $L_2(T;\mathbb{R}^{n})$
с оценкой точности вида (\ref{4_7}).

\setcounter{equation}{0}
\section{Настройка параметров алгоритма}
Количество $N$ доступных измерению траекторий исходного СДУ (\ref{2_1}) определяет как погрешность $h^{N}$ входных оценок (\ref{2_2}), (\ref{3_3}),
так и параметры алгоритма $\delta^{N}$, $\alpha_{m}^{N}$ и $\alpha_{d}^{N}$, обеспечивающие сходимость (\ref{4_6}), (\ref{4_7}).
Соотношения не прописаны явно, как это необходимо для расчетов, а задаются в виде порядка сходимости при $1/N \rightarrow 0$ с точностью до констант,
для которых оценки сверху могут быть выписаны явно, но сами значения нуждаются в уточнении.
Отметим подтвержденную экспериментально чувствительность алгоритма к ограничениям на функции,
формирующие исходную систему (\ref{2_1}), на множества допустимых значений возмущений
$S_{u1}$ и $S_{u2}$, на вариацию возмущений и др.
Необходимо получить соотношения (\ref{4_8}), обеспечивающие успешную работу алгоритма
по восстановлению в конкретной системе вида (\ref{2_1}) любых возмущений,
удовлетворяющих априорным ограничениям (в этом суть калибровки алгоритма).
На данном этапе полагаем, что <<оптимальные>> зависимости параметров определяются эмпирически через решение следующей экстремальной задачи:
\begin{equation}\label{5_1}
(K^{*}_{1},\,K^{*}_{2},\,K^{*}_{3})
 = \arg \min \left\{ \sum_{N \in I^{N}} \beta^{N}\| u^{N}_{K_{1},K_{2},K_{3}}(\cdot ) - u_*(\cdot) \|_{L_2(T)}:\, (K_{1},\,K_{2},\,K_{3}) \in S_{K} \right\}.
\end{equation}
Здесь $L_2(T)=L_2(T;\mathbb{R}^{r+n_d})$, $u^{N}_{K_{1},K_{2},K_{3}}(\cdot )$ --- выход алгоритма, зависящий от констант из (\ref{4_8});
$I^{N}$ --- множество рассматриваемых значений $N$, $\beta^{N}$ --- весовые коэффициенты, \linebreak $\sum_{N \in I^{N}} \beta^{N} =1$;
$S_{K}$ --- некоторое множество допустимых значений
вектора $(K_{1},\,K_{2},\,K_{3})$, как правило, параллелепипед в положительном октанте.
Фиксируется некоторая реальная функция $u_*(\cdot )$, вводится равномерная сетка по
$K_{1}$, $K_{2}$, $K_{3}$ и полным перебором по ней решается задача (\ref{5_1}),
что требует многочисленных прогонок вычислительного модуля
с моделированием большого числа независимых траекторий СДУ (\ref{2_1}) для формирования оценок (\ref{2_2}) для различных $N$.
Эксперименты подтверждают применимость решения (\ref{5_1}) для восстановления других реальных возмущений при сохранении заданных ограничений.
Для получения динамики уравнения (\ref{2_1}) используется метод Эйлера с заменой винеровского процесса последовательностью случайных импульсов
(аппроксимационная схема детально описана в \cite{mil}); его среднеквадратический порядок точности для квазилинейного СДУ
равен $O(\delta_s)$, где $\delta_s$ --- шаг моделирования. Этот шаг должен быть  мал относительно шага $\delta^{N}$,
с которым поступает  информация о траекториях и, соответственно, работает процедура восстановления неизвестных функций $u_1(t)$ и $u_2(t)$.
Таким образом, в процедуре подбора параметров алгоритма для конкретной системы задействованы две временные сетки,
причем моделирование с более мелким шагом имеет место только на стадии настройки алгоритма,
при решении реальной задачи входная информация должна поступать из <<внешнего мира>>.
Отметим, что вероятностный характер оценки сходимости (\ref{4_7}) предполагает усреднение выхода алгоритма
при фиксированном наборе параметров по серии запусков.

Вычислительные эксперименты показали, что описанная процедура, подразумевая множество прогонов на сетке по нескольким параметрам алгоритма,
требует значительных ресурсов при расчете на последовательной машине, однако допускает естественное (и эффективное) распараллеливание.
Оно базируется на стандартной схеме <<мастер-рабочий>> \cite{gerg} с единым загрузочным MPI-модулем, почти
равномерным распределением вычислительной нагрузки по всем ядрам и
отсутствием обменов между рабочими модулями.
Моделирование большого числа независимых траекторий СДУ с мелким временным шагом является единственной времяемкой частью процедуры,
поэтому оно распределяется поровну между всеми ядрами.
В момент поступления входной информации об оценках (\ref{2_2}), (\ref{estim1}), (\ref{estim3})
все модули проводят расчет одинаковых по размеру порций $x^{r}(\tau_{i})$ и усредняют
как сами величины, так и их квадраты, в соответствии с формулами (\ref{estim1}), (\ref{estim3}).
Затем рабочие отправляют свой вклад в оценки мастеру, который добавляет свою порцию и, усредняя полученные суммы по количеству задействованных ядер
(легитимность операции  обеспечивается одинаковостью порций),
окончательно формирует (\ref{estim1}) и (\ref{estim3}).
Все вычисления в соответствии с алгоритмом реконструкции  $(\ref{4_1})$--$(\ref{4_4})$
выполняются мастером. Блок-схема основного фрагмента параллельной процедуры приведена на рис.~1.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=6cm]{chart.eps}\\
\caption{Схема основного фрагмента параллельной процедуры.
Операции, выполняемые только мастером, помечены символом <<M>>, только рабочими --- символом <<W>>
(остальные операции выполняются на всех ядрах)}
\end{center}
\end{figure}

\setcounter{equation}{0}
\section{Модельный пример}
Рассмотрим квазилинейную систему, описывающую случайный процесс,
близкий к известному средне-возвратному процессу Орнштейна---Уленбека \cite{oksen}:
$$
dx_1(t) = (-2x_1(t) + 0.1x_2(t) + u_{11}(t) + u_{12}(t) + 1)\,dt + u_{21}(t)x_1(t)\,d\xi(t),
$$
$$
dx_2(t) = (0.1x_1(t) - x_2(t) + 2u_{12}(t) + 2)\,dt + u_{22}(t)x_2(t)\,d\xi(t),
$$
\begin{equation}\label{6_1}
t\in T=[0,1], \, x_1(0)=1, x_2(0)=1.
\end{equation}
Здесь $u(t)=(u_1(t),\,u_2(t))$, $u_1(t) = (u_{11}(t),\,u_{12}(t))$ и $u_2(t) = (u_{21}(t),\,u_{22}(t))$ --- неизвестные
возмущения, подлежащие восстановлению, $u_{1i}(t) \in [-2,2]$, $u_{2i}(t) \in [0,2]$, $i=1,2$;
$\xi(t)$ --- стандартный скалярный винеровский процесс.
Уравнения типа (\ref{6_1}) используются, в частности, в некоторых простейших моделях,
описывающих динамику численности относительно устойчивых популяций биологических особей \cite{oksen}.
В таком случае величина $x(t) = (x_1(t),\,x_2(t))$ представляет текущую численность (в условных единицах)
двух взаимодействующих видов; структура детерминированной части уравнения
определяет некоторые численности (например, средние величины за прошедший интервал),
стремление вернуться к которым <<подсознательно>> присутствует у популяций.
Такому возврату могут препятствовать соседи и  воздействие внешних факторов через возмущение
$u_1(t)$, влияющее на математическое ожидание процесса, и возмущение $u_2(t)$,
характеризующее амплитуду случайных шумов.
Измерению (в дискретные моменты времени) доступны реализации $x(t)$, соответствующие различным колониям.
Запишем для (\ref{6_1}) систему ОДУ, соответствующую (\ref{3_1}), (\ref{3_2}):
$$
\dot m_1(t) = -2m_1(t) + 0.1m_2(t) + u_{11}(t) + u_{12}(t) + 1,
$$
$$
\dot m_2(t) = 0.1m_1(t) - m_2(t) + 2u_{12}(t) + 2,
$$
$$
\dot d_1(t) = -4d_1(t) + 0.2d_2(t) + (d_1(t) + m_1^2(t))u_{31}(t),
$$
$$
\dot d_2(t) = 0.1d_1(t) - 3d_2(t) + 0.1d_3(t) + (d_2(t) + m_1(t)m_2(t))u_{32}(t),
$$
$$
\dot d_3(t) = 0.2d_2(t) - 2d_3(t) + (d_3(t) + m_2^2(t))u_{33}(t),
$$
\begin{equation}\label{6_2}
x_1(0)=1, x_2(0)=1, d_1(0)=0, d_2(0)=0, d_3(0)=0.
\end{equation}
Здесь $u_{31}=u_{21}^2$, $u_{32}=u_{21}u_{22}$, $u_{33}=u_{22}^2$.

Целью вычислительных экспериментов являлась адаптация алгоритма
к конкретной динамической системе (\ref{6_1})--(\ref{6_2}), т.е. получение универсального набора
зависимостей его параметров $\delta^{N}$, $\alpha_{m}^{N}$ и $\alpha_{d}^{N}$ от $N$ через $K_{1}$, $K_{2}$, $K_{3}$
(см.~(\ref{4_8})) для восстановления любых возмущений, удовлетворяющих заданным ограничениям.
В первой серии для решения (\ref{5_1}) были выбраны функции
$$
u_{1}(t) = (\sin{10t}, 1), \, u_{2}(t)= (\sqrt{t}, 10\sqrt{t^3/3 - t^2/2 + 3t/16})
$$
и параметры $K_i \in [0.1, 5]$, $i=1,2,3$, шаг сетки по $K_i$ равный $0.05$, $I^{N}=\{10^{3},10^{6}\}$, $\beta^{N}=\{0.5,0.5\}$.
Были получены оптимальные значения констант (в смысле (\ref{5_1})) в соотношениях (\ref{4_8})
для системы (\ref{6_1}) и принятых ограничений на возмущения: $K^{*}_{1}=0.6$, $K^{*}_{2}=0.35$, $K^{*}_{3}=3.5$.
Результаты восстановления функций $u_1(t)$ и $u_2(t)$, полученные для найденных $K_i$ и различных $N$,
в принципе согласуются с основным утверждением статьи, подтверждая сходимость вида (\ref{4_6}): чем больше $N$, тем меньше погрешность реконструкции,
см.~рис.~2, 3.
\begin{figure}[h!]
\begin{center}
\begin{tabular}{p{5.0cm}p{5.0cm}}
\includegraphics[width=4.7cm]{bad_u11.eps} & \includegraphics[width=4.7cm]{bad_u12.eps}\\
\includegraphics[width=4.7cm]{bad_u21.eps} & \includegraphics[width=4.7cm]{bad_u22.eps}
\end{tabular}
\caption{Восстановление возмущений $u_{11}$, $u_{12}$ (верхний ряд) и $u_{21}$, $u_{22}$ (нижний ряд):
реальные функции изображены пунктирной линией, результаты их восстановления~--- сплошной;
горизонтальная ось соответствует времени $t$, вертикальная~--- значениям возмущений. Параметры:
$N=10^{3}$, $\delta ^{N} =0.025$, $\delta_s = \delta ^{N}/10^{2}$, $\alpha_{m}^{N} =0.04$, $\alpha_{d}^{N} =0.4$,
погрешность $\| u^{N}(\cdot ) - u_*(\cdot) \|_{L_2(T)}= 0.868$}
\end{center}
\end{figure}
%\bigskip
Затем полученные соотношения параметров были протестированы на других возмущениях,
удовлетворяющих заданным геометрическим ограничениям.
На рис.~4 приведен пример восстановления $u_{2}(t) = (\sqrt{2}(1-2|t-0.5|), 1+ 0.5\sin{30t})$
со значениями параметров алгоритма из варианта, представленного на рис.~3.
Точность реконструкции сравнима с результатами, полученными
для тестового (в смысле (\ref{5_1})) возмущения.
Отметим, что вторая координата возмущения $u_{2}$ восстанавливается хуже первой, что
объясняется большей амплитудой шума во втором уравнении системы (\ref{6_1}).
\begin{figure}[h!]
\begin{center}
\begin{tabular}{p{5.0cm}p{5.0cm}}
\includegraphics[width=4.7cm]{good_u11.eps} & \includegraphics[width=4.7cm]{good_u12.eps}\\
\includegraphics[width=4.7cm]{good_u21.eps} & \includegraphics[width=4.7cm]{good_u22.eps}
\end{tabular}
\caption{Восстановление возмущений $u_{11}$, $u_{12}$ (верхний ряд) и $u_{21}$, $u_{22}$ (нижний ряд):
реальные функции изображены пунктирной линией, результаты их восстановления~--- сплошной;
горизонтальная ось соответствует времени $t$, вертикальная~--- значениям возмущений. Параметры:
$N=10^{6}$, $\delta ^{N} =0.001$, $\delta_s = \delta ^{N}/10^{2}$, $\alpha_{m}^{N} =0.005$, $\alpha_{d}^{N} =0.05$,
погрешность $\| u^{N}(\cdot ) - u_*(\cdot) \|_{L_2(T)}= 0.118$}
\end{center}
\end{figure}
%\vspace{-1cm}

Для иллюстрации показателей, характеризующих качество распараллеливания,
был выбран вариант из рис.~3 с симуляцией 128000 траекторий уравнения (\ref{6_1}).
Моделирование проводилось в Институте математики и механики им. Н.Н.Красовского
УрО РАН на гибридном суперкомпьютере <<Уран>>, который состоит из 196 вычислительных узлов,
оснащенных процессорами Intel Xeon и графическими ускорителями NVIDIA Tesla,
и имеет пиковую производительность более 215 Тфлопс.
Отметим, что в общей сложности пользователям доступно 1940 вычислительных ядер CPU, 314 платы GPU и 4~Тб
оперативной памяти. В наших экспериментах использовался компактный раздел кластера, называемый
$Apollo$ (16 узлов по два 18-и ядерных процессора Intel Xeon CPU E5-2697 с частотой 2.30~ГГц,
оперативная память 256~Гб, локальный жесткий диск 1~Тб); его узлы объединены
высокоскоростной сетью Infiniband нового поколения со скоростью передачи данных 100 Гбит/с.

\newpage
\begin{figure}[h!]
\begin{center}
\begin{tabular}{p{5.0cm}p{5.0cm}}
\includegraphics[width=4.7cm]{alt_u21.eps} & \includegraphics[width=4.7cm]{alt_u22.eps}
\end{tabular}
\caption{Восстановление альтернативных возмущений $u_{21}$ и $u_{22}$:
реальные функции изображены пунктирной линией, результаты их восстановления~--- сплошной;
горизонтальная ось соответствует времени $t$, вертикальная~--- значениям возмущений. Параметры:
$N=10^{6}$, $\delta ^{N} =0.001$, $\delta_s = \delta ^{N}/10^{2}$, $\alpha_{m}^{N} =0.005$, $\alpha_{d}^{N} =0.05$,
погрешность $\| u^{N}(\cdot ) - u_*(\cdot) \|_{L_2(T)}= 0.169$}
\end{center}
\end{figure}
\vspace{-0.5cm}

Результаты тестирования представлены на рис.~5, где приведены графики зависимостей
ускорения и эффективности параллельного алгоритма от количества ядер.
Напомним, что  ускорение $S_{p}$ и эффективность $E_{p}$ вычисляются,
соответственно, по формулам: $S_{p} =T_{1}/T_{p}$, $E_{p}=S_{p}/p$,
$T_{p}$ --- время выполнения программы (для одного набора параметров) на кластере
для $p$ ядер, $T_{1}$ --- время работы последовательного алгоритма.
Для <<идеального>> параллельного алгоритма ускорение равно
числу ядер; в этом случае имеем единичную эффективность.
\begin{figure}[h]
\begin{center}
\begin{tabular}{p{5.0cm}p{5.0cm}}
\includegraphics[width=4.7cm]{acce.eps} & \includegraphics[width=4.7cm]{effi.eps}
\end{tabular}
\caption{Ускорение $S_{p}$ (слева) и эффективность $E_{p}$ (справа) в зависимости от количества ядер:
реальные характеристики изображены сплошной линией, идеальные величины~--- пунктирной;
горизонтальная ось соответствует количеству ядер $p$, вертикальная~--- значениям $S_{p}$ и $E_{p}$}
\end{center}
\end{figure}
%\vspace{-0.5cm}

Из рис.~5 следует, что эффективность распараллеливания достаточно высока,
причем с ростом числа задействованных процессоров она не падает ниже
разумного уровня (например, при $p=128$ имеем $S_{p}=103.5$ и $E_{p}=0.81$).
Это важно, поскольку для настройки алгоритма на конкретную систему и его тестирования
(в вероятностном смысле) требуется достаточно большое количество запусков программы.

\section*{Заключение}
В работе рассмотрена обратная задача для
квазилинейной стохастической системы с диффузией, зависящей от фазового состояния.
Постановка предполагает динамическую реконструкцию двух неизвестных неслучайных возмущений,
входящих в детерминированную и стохастическую части системы,
на основе точных дискретных измерений некоторого количества реализаций
искомого случайного процесса.
Задача для СДУ заменяется обратной задачей для системы двух ОДУ, которым удовлетворяют
математическое ожидание и ковариационная матрица исходного процесса.
Для конструктивного алгоритма восстановления
получена оценка точности относительно количества доступных измерению реализаций.
Ключевым результатом статьи следует считать разработку процедуры адаптации алгоритма
к конкретной динамической системе и ограничениям на ее структурные элементы;
это обусловило необходимость прогонки вычислительного модуля на сетке параметров и моделирования
большого числа траекторий СДУ для получения статистических оценок моментов.
С целью экономии вычислительных ресурсов  введено естественное и достаточно эффективное
распараллеливание основных времяемких процедур.

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

\begin{biblio}

\bibitem{kryazh}
Кряжимский~А.В., Осипов~Ю.С. О моделировании управления в динамической
системе~// Изв.\ АН СССР. Техн. кибернетика. 1983. №~2. С.~51--60.

\bibitem{oskr2011}
Осипов~Ю.С., Кряжимский~А.В., Максимов~В.И. Методы динамического восстановления входов управляемых систем.
Екатеринбург: Изд-во УрО РАН, 2011. 291~с.

\bibitem{maxim2000}
Максимов~В.И. Задачи динамического восстановления входов бесконечномерных
систем. Екатеринбург: Изд-во УрО РАН, 2000. 305~с.

\bibitem{kras}
Красовский~Н.Н., Субботин~А.И. Позиционные дифференциальные игры. М.:~Наука, 1984. 456~с.

\bibitem{tikh}
Тихонов~А.Н., Арсенин~В.Я. Методы решения некорректных задач. М.:~Наука, 1978. 142~с.

\bibitem{maxbliz}
Близорукова~М.С., Максимов~В.И.
Об одном алгоритме реконструкции траектории и управления в системе с запаздыванием~//
Тр. Ин-та математики и механики  УрО РАН.
Екатеринбург: ИММ УрО РАН, 2012. Т.~18, №~1. C.~109--122.

\bibitem{maxim2007}
Кадиев~А.М., Максимов~В.И. О реконструкции управлений в параболическом уравнении~// Дифференц. уравнения.
2007. Т.~43, №~11. С.~1545--1552.

\bibitem{oskr1989}
Кряжимский~А.В., Осипов~Ю.С.
Об устойчивом позиционном восстановлении управления по измерениям части координат~//
Некоторые задачи управления и устойчивости. Свердловск: Изд-во УрО АН СССР, 1989. С.~33--47.

\bibitem{rozen2007}
Розенберг~В.Л. Задача динамического восстановления неизвестной функции в линейном
стохастическом дифференциальном уравнении~// Автом. и телемех. 2007. №~11. С.~76--87.

\bibitem{rozen2016}
Розенберг~В.Л. Восстановление амплитуды случайной помехи в
линейном стохастическом уравнении по измерениям части координат~//
Журнал вычисл. математики и мат. физики.
2016. Т.~56, №~3. С.~377--386.

\bibitem{rozen2018}
Rozenberg~V.L. Dynamical Input Reconstruction Problem for a Quasi-linear
Stochastic System~// Proceedings of the 17th IFAC Workshop on Control Applications of
Optimization, CAO 2018, October~15--19, 2018, Yekaterinburg, Russia. P.~727--732.

\bibitem{oksen}
Оксендаль~Б. Стохастические дифференциальные уравнения. Введение
в теорию и приложения. М.:~Мир, 2003. 408~с.

\bibitem{king}
Королюк~В.С., Портенко~Н.И., Скороход~А.В., Турбин~А.Ф.
Справочник по теории вероятностей и математической статистике. М.:~Наука, 1985. 640~с.

\bibitem{chern}
Черноусько~Ф.Л., Колмановский~В.Б. Оптимальное управление при случайных возмущениях. М.:~Наука, 1978. 272~с.

\bibitem{mil}
Мильштейн~Г.Н. Численное интегрирование стохастических дифференциальных уравнений.
Свердловск:~Изд-во УрГУ, 1988. 224~с.

\bibitem{gerg}
Гергель~В.П. Теория и практика параллельных вычислений. М.: БИНОМ, 2007. 424~с.

\end{biblio}

% Сведения об авторах (Ф.И.О.(полностью), степень, звание, кафедра (отдел), вуз, e-mail)
Мельникова Лидия Антоновна, математик 1-й кат., Институт математики и механики им. Н.Н.~Красовского УрО РАН (Екатеринбург, Российская Федерация)

Розенберг Валерий Львович, к.ф.-м.н., с.н.с., Институт математики и механики им. Н.Н.~Красовского УрО РАН (Екатеринбург, Российская Федерация)


%\newpage

{\phantom{eee} \hrule \vskip 4 pt}
\classify{} % Код MSC не нужен
%\classify{MSC 00000} % Коды MSC см. http://www.ams.org/msc/

\title{ALGORITHM OF DYNAMICAL INPUT RECONSTRUCTION FOR A STOCHASTIC DIFFERENTIAL EQUATION: TUNING OF PARAMETERS AND NUMERICAL EXPERIMENTS}

\author{\copyright~2019 L.A.~Melnikova, V.L.~Rozenberg}

\address{Krasovskii Institute of Mathematics and Mechanics,
               Ural Branch of the Russian Academy of Sciences (S.~Kovalevskoi str.~16, Yekaterinburg, 620990 Russia )}

\email{E-mail: \href{mailto:meln@imm.uran.ru}{meln@imm.uran.ru},
\href{mailto:rozen@imm.uran.ru}{rozen@imm.uran.ru}}

\received{Received: 16.04.2019}

\maketitle{}

\begin{abstract}%
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 algorithm's parameters in order to have best approximation results
for different disturbances satisfying a priori constraints in a specific dynamical system is proposed.
Desired dependencies of algorithm's parameters on the number of measured realizations are determined
empirically via solving a specific extremal problem, where the deviation of 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.

\keywordsen{
stochastic differential equation, dynamical reconstruction,
controlled model, parameter tuning, parallelization of calculations.}
\end{abstract}

\begin{flushleft}
	\textbf{FOR CITATION}
\end{flushleft}
\justifying

\begin{citationplace}
Melnikova~L.A., Rozenberg~V.L.
Algorithm of Dynamical Input Reconstruction for a Stochastic Differential Equation: Tuning of Parameters and Numerical Experiments.
{\it Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering.} 2019. vol.~X, no.~Y. pp.~Z1--Z2. (in Russian) DOI: 10.14529/cmseXXXXXX.
\end{citationplace}

\begin{biblio_lat}
%\itemsep -2pt

\bibitem{kryazhen}
Kryazhimskii~A.V., Osipov~Yu.S.
Modelling of a Control in a Dynamic System.
Engrg. Cybernetics. 1984. vol.~21, no.~2, pp.~38--47.

\bibitem{oskr2011en}
Osipov~Yu.S., Kryazhimskii~A.V., Maksimov~ V.I.
\textit{Metody dinamicheskogo vosstanovleniya vkhodov upravlyaemykh sistem}
[Dynamic Recovery Methods for Inputs of Control Systems].
Yekaterinburg, Publishing of the UB of RAS, 2011. 291~p.

\bibitem{maxim2000en}
Maksimov~V.I. Dynamical Inverse Problems of Distributed Systems.
VSP, Utrecht--Boston, 2002. 269~p.

\bibitem{krasen}
Krasovskii~N.N., Subbotin~A.I. Game-Theoretical Control Problems.
Springer, New York, 1988. 517~p.

\bibitem{tikhen}
Tikhonov~A.N., Arsenin~V.Ya. Solutions of Ill-Posed Problems.
Wiley, New York, 1981. 258~p.

\bibitem{maxblizen}
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, suppl.~1, pp.~S66--S79.

\bibitem{maxim2007en}
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.

\bibitem{oskr1989en}
Kryazhimskii~A.V., Osipov~Yu.S.
On a Stable Positional Recovery of Control from Measurements of a Part of Coordinates.
\textit{Nekotorye zadachi upravleniya i ustoichivosti}
[Some Problems of Control and Stability].
Sverdlovsk, Publishing of the UB of Academy of Sciences of the USSR, 1989. P.~33--47.

\bibitem{rozen2007en}
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.

\bibitem{rozen2016en}
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.

\bibitem{rozen2018en}
Rozenberg~V.L. Dynamical Input Reconstruction Problem for a Quasi-linear
Stochastic System. Proceedings of the 17th IFAC Workshop on Control Applications of
Optimization (CAO 2018, Yekaterinburg, Russia, October~15--19, 2018). P.~727--732.

\bibitem{oksenen}
\O ksendal~B. Stochastic Differential Equations: an Introduction with Applications.
Springer, Berlin, 1985. 408~p.

\bibitem{kingen}
Korolyuk~V.S., Portenko~N.I., Skorokhod~A.V., Turbin~A.F.
\textit{Spravochnik po teorii veroyatnostei i matematicheskoi statistike}
[Handbook on Probability Theory and Mathematical Statistics].
Moscow, Nauka, 1985. 640~p.

\bibitem{chernen}
Chernous'ko~F.L., Kolmanovskii~V.B. \textit{Optimal'noe upravlenie pri sluchainykh vozmushcheniyakh}
[Optimal Control under Random Perturbation].
Moscow, Nauka, 1978. 272~p.

\bibitem{milen}
Milshtein~G.N. \textit{Chislennoe integrirovanie stokhasticheskikh differentsial'nykh uravnenii}
[Numerical Integration of Stochastic Differential Equations].
Sverdlovsk, Publishing of the Ural State University, 1988. 224~p.

\bibitem{gergen}
Gergel'~V.P. \textit{Teoriya i praktika parallel'nykh vychislenii}
[Theory and Practice of Parallel Computing].
Moscow, Binom, 2007. 424~p.

\end{biblio_lat}

\end{document}
