Журнал вычислительной математики и математической физики, 2022, T. 62, № 11, стр. 1840-1850

1. ВВЕДЕНИЕ

Рассматриваемая задача реконструкции вкладывается в проблематику обратных задач динамики управляемых систем, которые, как правило, являются некорректными и требуют применения регуляризирующих процедур. Для ее решения используется ставший классическим подход, предложенный в работах Ю.С. Осипова и его коллег (см. [1]–[6] и библиографию в них) и получивший название метода динамического обращения. Он основан на сочетании принципов теории позиционного управления (см. [7]) и идей теории некорректных задач (см. [8]). Суть его в сведении задачи реконструкции к задаче управления по принципу обратной связи (часто на основе принципа экстремального прицеливания Н.Н. Красовского, см. [7]) вспомогательной динамической системой (моделью). Адаптация модельного управления к результатам текущих наблюдений обеспечивает аппроксимацию (в подходящем смысле) неизвестного входа. Метод динамического обращения был многократно реализован для различных детерминированных систем, описываемых обыкновенными дифференциальными уравнениями (ОДУ), дифференциально-функциональными уравнениями, уравнениями и вариационными неравенствами с распределенными параметрами и др. Были созданы устойчивые алгоритмы, работающие для некоторых классов частично наблюдаемых систем (в случае конечномерной системы роль входного сигнала могли играть измерения части координат фазового вектора, а в случае бесконечномерной – значения решения на некоторых подмножествах области определения). Такие задачи были сформулированы и решены, например, в [4], [5]. Касательно стохастических объектов отметим, что впервые задача позиционного моделирования неизвестного стохастического управляющего воздействия в системе, описываемой ОДУ, была рассмотрена в [9]. Впоследствии поле применения методов динамической регуляризации было расширено для систем с неопределенными входами в [10], где рассмотрена задача об устойчивой аппроксимации неизвестного входа управляемой системы по результатам неточных наблюдений ее фазовых состояний при условии, что помехи в канале наблюдения подчинены некоторому вероятностному распределению с малым математическим ожиданием. Впервые в рамках теории динамического обращения использован вероятностный критерий качества и построен динамический алгоритм реконструкции нормального входа, доставляющий сколь угодно высокую точность среднеквадратической аппроксимации с вероятностью, сколь угодно близкой к единице.

Настоящая работа фактически продолжает исследования задач реконструкции неизвестных параметров линейных и квазилинейных стохастических дифференциальных уравнений (СДУ), начатые в [11] с рассмотрения задачи динамического восстановления возмущения, входящего в интеграл Ито и характеризующего амплитуду случайной помехи, на основе измерения реализаций всего фазового вектора линейного СДУ. Идея решения состояла в переходе к задаче для линейной системы ОДУ, которой удовлетворяют математическое ожидание и ковариационная матрица исходного процесса. В [12] исследовалась достаточно общая постановка обратной задачи для линейного СДУ, предполагающая реконструкцию возмущений, входящих и в детерминированный, и в стохастический члены уравнения, на основе дискретной по времени информации о некотором количестве реализаций части координат случайного процесса. Обоснована применимость алгоритма восстановления неизвестных параметров, разработанного ранее для частично наблюдаемой системы ОДУ (см. [4], [5]), и предложена соответствующая модификация. Анализ задачи реконструкции неизвестных входов квазилинейного СДУ при измерении реализаций всего фазового вектора был проведен в [13]. Отличительной особенностью задачи для квазилинейного СДУ является нелинейность системы ОДУ, описывающей динамику математического ожидания и ковариационной матрицы процесса. Настоящая статья призвана дополнить исследование задач динамической реконструкции для квазилинейных СДУ случаем неполной информации, когда измеряется лишь часть координат фазового состояния. Отметим, что обратная задача для системы гибридного типа, состоящей из квазилинейного СДУ и ОДУ, при условии измерения реализаций случайного процесса рассмотрена в [14].

2. ПОСТАНОВКА ЗАДАЧИ

Рассматривается квазилинейное (в соответствии с терминологией [15] и других работ этих авторов) СДУ с диффузией, зависящей от фазового состояния:

Здесь $x = ({{x}_{1}}, \ldots ,{{x}_{n}}) \in {{\mathbb{R}}^{n}}$, $\xi \in \mathbb{R}$; ${{x}_{0}}$ – известный детерминированный или случайный (нормально распределенный, с независимыми координатами) вектор начальных условий; $\omega \in \Omega $, $(\Omega ,F,P)$ – вероятностное пространство (см., например, [16]); $\xi (t,\omega )$ – стандартный скалярный винеровский процесс (т.е. выходящий из нуля процесс с нулевым математическим ожиданием и дисперсией, равной $t$); $f(t) = \{ {{f}_{i}}(t)\} $, $A(t) = \{ {{a}_{{ij}}}(t)\} $ и $B(t) = \{ {{b}_{{ij}}}(t)\} $ – непрерывные матричные функции размерности $n \times 1$, $n \times n$, и $n \times r$ соответственно. На уравнение действуют два внешних возмущения: вектор ${{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}}xd\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}}({\kern 1pt} \cdot {\kern 1pt} ) \in {{L}_{2}}(T;{{\mathbb{R}}^{r}})$ и ${{u}_{2}}({\kern 1pt} \cdot {\kern 1pt} ) \in {{L}_{2}}(T;{{\mathbb{R}}^{n}})$ принимают значения из заданных выпуклых компактов ${{S}_{{{{u}_{1}}}}}$ и ${{S}_{{{{u}_{2}}}}}$.

Решение уравнения (2.1) определяется как случайный процесс, удовлетворяющий при любом $t$ с вероятностью 1 соответствующему интегральному тождеству, содержащему в правой части стохастический интеграл Ито. Как известно, при сделанных предположениях существует единственное решение, являющееся нормальным марковским процессом с непрерывными реализациями (см. [17]). Отметим, что уравнения типа (2.1) описывают простейшие линеаризованные модели, например, изменения численности многовидовой биологической популяции в стохастической среде или динамики цен на товарных рынках при влиянии случайных факторов.

Изучаемая задача состоит в следующем. В дискретные, достаточно частые, моменты времени ${{\tau }_{i}} \in T$, ${{\tau }_{i}} = i\delta $, $\delta = \vartheta {\text{/}}l$, $i \in [0:(l - 1)],$ поступает информация о некотором количестве $N$ реализаций случайного процесса $x({{\tau }_{i}})$, причем измерению доступны только $q$ ($q \leqslant n$) первых координат, т.е. вектор $({{x}_{1}}, \ldots ,{{x}_{q}})$. Полагаем, что $l = l(N)$ и существуют оценки $m_{{qi}}^{N}$ $q$-подвектора ${{m}_{q}}(t) = \{ {{m}_{j}}(t)\} ,$ ${\kern 1pt} j \in [1:q]$, вектора математического ожидания процесса $m(t) = Mx(t)$ и $D_{{qi}}^{N}$ ($q \times q$)-подматрицы ${{D}_{q}}(t) = \{ {{d}_{{jp}}}(t)\} ,$ $j,p \in [1:q]$, ковариационной матрицы $D(t) = M(x(t) - m(t))(x(t) - $ $ - \;m(t)){\kern 1pt} '$ (штрих означает транспонирование) такие, что выполняется соотношение

причем $h(N),{\kern 1pt} g(N) \to 0$ при $N \to \infty $. Стандартные статистические процедуры (см. [18]) построения оценок $m_{{qi}}^{N}$ и $D_{{qi}}^{N}$ допускают модификации, обеспечивающие выполнение (2.2).

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

3. ПЕРЕХОД К ЗАДАЧЕ ДЛЯ СИСТЕМЫ ОДУ

Применяя ставший классическим метод моментов из [19], по аналогии с [12], [13], сведем задачу восстановления для СДУ к задаче для системы ОДУ. Введем обозначения ${{m}_{0}} = M{{x}_{0}}$, ${{D}_{0}} = M({{x}_{0}} - {{m}_{0}})({{x}_{0}} - {{m}_{0}}){\kern 1pt} '$. Отметим, что матрица ${{D}_{0}}$ является диагональной (ввиду независимости координат вектора ${{x}_{0}}$). В силу квазилинейности исходного уравнения и равенства нулю математического ожидания интеграла Ито величина $m(t)$ зависит только от ${{u}_{1}}(t)$; ее динамика описывается уравнением

Напомним, что измерению доступны первые $q$ координат исходного $n$-мерного вектора $x$, что обеспечивает получение оценки (2.2) для первых $q$ координат вектора $m$; количество неизмеряемых координат равно $n - q$. Используя стандартную схему теории динамического обращения (см. [4], [5]), запишем уравнение (3.1) в виде системы с разделением измеряемой и неизмеряемой компонент. Введем следующие матрицы:

Тогда $q$-мерный вектор $C_{1}^{m}m$ – измеряемая часть $m$, а $(n - q)$-мерный вектор $D_{1}^{m}m$ – его неизмеряемая часть. Обозначим ${{y}_{m}} = C_{1}^{m}m$, ${{z}_{m}} = D_{1}^{m}m$. Умножая уравнение (3.1) поочередно на матрицы $C_{1}^{m}$ и $D_{1}^{m}$, учитывая равенство $m = C_{2}^{m}{{y}_{m}} + D_{2}^{m}{{z}_{m}}$, приходим к системе с разделением измеряемой и неизмеряемой компонент: где $A_{{1y}}^{m} = C_{1}^{m}AC_{2}^{m}$, $A_{{1z}}^{m} = C_{1}^{m}AD_{2}^{m}$, $B_{1}^{m} = C_{1}^{m}B$, ${{y}_{{m0}}} = C_{1}^{m}{{m}_{0}}$, $A_{{2y}}^{m} = D_{1}^{m}AC_{2}^{m}$, $A_{{2z}}^{m} = D_{1}^{m}AD_{2}^{m}$, $B_{2}^{m} = D_{1}^{m}B$, ${{z}_{{m0}}} = D_{1}^{m}{{m}_{0}}$.

Ковариационная матрица $D(t)$ не зависит от ${{u}_{1}}(t)$ явно. Для описания ее динамики используется схема вывода уравнения метода моментов, примененная в [12], [13]. В результате для $D \in {{\mathbb{R}}^{{n \times n}}}$ имеем следующее уравнение:

Матричное уравнение (3.4) переписываем в виде более традиционного для рассматриваемых задач векторного уравнения, размерность которого, с учетом симметричности матрицы $D(t)$, определяется как ${{n}_{{yz}}} = ({{n}^{2}} + n){\text{/}}2$. Вводится вектор $d(t) = \{ {{d}_{s}}(t)\} ,$ $s \in [1:{{n}_{{yz}}}]$, состоящий из последовательно записанных и пронумерованных элементов матрицы $D(t)$, взятых построчно, начиная с элементов, расположенных на главной диагонали; его координаты находятся по элементам матрицы $D(t) = \{ {{d}_{{ij}}}(t)\} ,$ $i,j \in [1:n]$:

Отметим, что соотношения (3.5) между индексами $s$ и $i,j$ взаимно однозначны ввиду способа нумерации элементов части матрицы $D(t)$. Процедура, аналогичная подробно описанной в [20], позволяет преобразовать уравнение (3.4) к виду

где матрица $\bar {A}(t):T \to {{\mathbb{R}}^{{{{n}_{{yz}}} \times {{n}_{{yz}}}}}}$ выписывается явно (см. [20]) $\bar {m}(t) = F(m(t))$, $\bar {u}(t) = F({{u}_{2}}(t))$, $F:{{\mathbb{R}}^{n}} \to {{\mathbb{R}}^{{{{n}_{{yz}}}}}}$, $\forall a \in {{\mathbb{R}}^{n}}{\kern 1pt} F(a) = \bar {a} = \{ {{\bar {a}}_{s}}\} $, $s \in [1:{{n}_{{yz}}}]$, ${{\bar {a}}_{s}} = {{a}_{i}}{{a}_{j}}$, $i \leqslant j$, $s = (n - i{\text{/}}2)(i - 1) + j$ (фактически отображение $F$ удлиняет $n$-мерный вектор до размерности вектора $d(t) \in {{\mathbb{R}}^{{{{n}_{{yz}}}}}}$), диагональная матрица $\bar {B}(d(t),\bar {m}(t)):T \to {{\mathbb{R}}^{{{{n}_{{yz}}} \times {{n}_{{yz}}}}}}$ находится из формул Из существования, единственности и вида решения уравнения (2.1) следуют существование и единственность решения системы (3.1), (3.6). Очевидно, что вектор $\bar {u}({\kern 1pt} \cdot {\kern 1pt} ) \in {{L}_{2}}(T;{{\mathbb{R}}^{{{{n}_{{yz}}}}}})$ принимает значения из некоторого выпуклого компакта ${{S}_{{\bar {u}}}} \in {{\mathbb{R}}^{{{{n}_{{yz}}}}}}$, который естественным образом строится по компакту ${{S}_{{{{u}_{2}}}}} \in {{\mathbb{R}}^{n}}$ с учетом того, что координаты вектора $\bar {u}({\kern 1pt} \cdot {\kern 1pt} )$ являются попарными произведениями координат вектора ${{u}_{2}}({\kern 1pt} \cdot {\kern 1pt} )$. Полагаем, что начальное состояние ${{d}_{0}}$ и измерения $d_{i}^{N}$ получаются из ${{D}_{0}}$ и $D_{i}^{N}$ в соответствии с формулой (3.5). Таким образом, формально по информации о векторе $d(t)$ будем восстанавливать вектор $\bar {u}(t)$. Затем, ограничиваясь рассмотрением тех его координат, которые равны $u_{{2i}}^{2}$, $i \in [1:n]$, мы покажем, что реальная величина ${{u}_{2}}(t)$ может быть восстановлена при дополнительных, достаточно естественных, предположениях.

Измерения первых $q$ координат вектора $x$ обеспечивают оценку типа (2.2) для $({{q}^{2}} + q){\text{/}}2$ координат вектора $d$, номера которых находятся из (3.5) и, в общем случае, не идут подряд. Обозначим ${{n}_{y}} = ({{q}^{2}} + q){\text{/}}2$, ${{n}_{z}} = {{n}_{{yz}}} - {{n}_{y}}$. Определим ${{I}_{y}}$ как ${{n}_{y}}$-мерный упорядоченный по возрастанию массив индексов, соответствующих измеряемым координатам вектора $d$:

а также, аналогичным образом, ${{I}_{z}}$ как ${{n}_{z}}$-мерный упорядоченный по возрастанию массив индексов, соответствующих неизмеряемым координатам вектора $d$.

Введем следующие матрицы:

Тогда ${{n}_{y}}$-мерный вектор $C_{1}^{d}d$ – измеряемая часть $d$, а ${{n}_{z}}$-мерный вектор $D_{1}^{d}d$ – его неизмеряемая часть. Обозначим ${{y}_{d}} = C_{1}^{d}d$, ${{z}_{d}} = D_{1}^{d}d$. Умножая уравнение (3.6) поочередно на матрицы $C_{1}^{d}$ и $D_{1}^{d}$, учитывая равенство $d = C_{2}^{d}{{y}_{d}} + D_{2}^{d}{{z}_{d}}$, приходим к системе с разделением измеряемой и неизмеряемой компонент: где $A_{{1y}}^{d} = C_{1}^{d}\bar {A}C_{2}^{d}$, $A_{{1z}}^{d} = C_{1}^{d}\bar {A}D_{2}^{d}$, ${{y}_{{d0}}} = C_{1}^{d}{{d}_{0}}$, $A_{{2y}}^{d} = D_{1}^{d}\bar {A}C_{2}^{d}$, $A_{{2z}}^{d} = D_{1}^{d}\bar {A}D_{2}^{d}$, ${{z}_{{d0}}} = D_{1}^{d}{{d}_{0}}$, $B_{1}^{d} = C_{1}^{d}\bar {B}(d,\bar {m})$, $B_{2}^{d} = D_{1}^{d}\bar {B}(d,\bar {m})$, $\bar {m} = F(m) = F(C_{2}^{m}{{y}_{m}} + D_{2}^{m}{{z}_{m}})$.

Отметим, что ввиду специфики матрицы $\bar {B}$, $B_{1}^{d}$ зависит только от ${{y}_{d}}(t)$, ${{y}_{m}}(t)$, а $B_{2}^{d}$ – только от ${{z}_{d}}(t)$, ${{z}_{m}}(t)$.

Теперь для систем (3.2), (3.3) и (3.8), (3.9) можно переформулировать исходную задачу восстановления. По ходу развития процесса в дискретные моменты времени ${{\tau }_{i}} \in T$, ${{\tau }_{i}} = i\delta $, $\delta = \vartheta {\text{/}}l(N)$, $i \in [0:(l(N) - 1)]$ поступает информация, позволяющая оценить части фазового состояния указанных систем, соответственно, векторы ${{y}_{m}}({{\tau }_{i}})$ и ${{y}_{d}}({{\tau }_{i}})$. Полагаем, что выполняется следующее соотношение, соответствующее (2.2):

где оценочные векторы $\xi _{{mi}}^{N} \in {{\mathbb{R}}^{q}}$ и $\xi _{{di}}^{N} \in {{\mathbb{R}}^{{{{n}_{y}}}}}$ естественным образом получаются из оценок $m_{{qi}}^{N}$ и $D_{{qi}}^{N}$, а $h(N) \to 0$ и $g(N) \to 0$ при $N \to \infty $.

Требуется указать алгоритм динамического восстановления неизвестных возмущений ${{u}_{1}}(t)$ и $\bar {u}(t)$ по информации (3.10), причем вероятность сколь угодно малого отклонения приближений от искомых входов в метрике соответственно пространств ${{L}_{2}}(T;{{\mathbb{R}}^{r}})$ и ${{L}_{2}}(T;{{\mathbb{R}}^{{{{n}_{{yz}}}}}})$ должна быть близка к 1 при достаточно большом $N$ и специальным образом согласованном с $N$ шаге временной дискретизации $\delta = \delta (N) = \vartheta {\text{/}}l(N)$.

В такой формулировке задача соответствует задаче, рассмотренной в [4] для случая измерения части координат ОДУ. В настоящей работе показано, что алгоритм, предложенный в [4], применим к решению задачи, полученной для систем (3.2), (3.3) и (3.8), (3.9), поскольку допускает конструктивное согласование своих параметров с количеством доступных измерению реализаций исходного случайного процесса, при этом легко проверяемые достаточные условия разрешимости задачи фактически формулируются в терминах исходного уравнения (2.1).

4. АЛГОРИТМ ВОССТАНОВЛЕНИЯ НЕИЗВЕСТНЫХ ВОЗМУЩЕНИЙ

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

Условие 1. Размерность неизвестной вектор-функции ${{u}_{1}}({\kern 1pt} \cdot {\kern 1pt} )$ не превосходит размерности компоненты ${{y}_{m}}({\kern 1pt} \cdot {\kern 1pt} )$ ($r \leqslant q$), и при всех $t \in T$ матрица $B_{1}^{m}(t)$ размерности $q \times r$ имеет ранг, равный $r$, т.е. является матрицей полного ранга.

Условие 2. Неизвестная вектор-функции ${{u}_{2}}({\kern 1pt} \cdot {\kern 1pt} )$ имеет следующую структуру: $\exists k \leqslant q$ $\forall i \in [(k + 1):n]$ $\exists j \in [1:k]:{{u}_{{2i}}} = {{u}_{{2j}}}$.

Таким образом, если количество неизвестных компонент функции ${{u}_{2}}({\kern 1pt} \cdot {\kern 1pt} )$ не превосходит $q$, тогда количество неизвестных компонент функции $\bar {u}({\kern 1pt} \cdot {\kern 1pt} )$, равное ${{n}_{k}} = ({{k}^{2}} + k){\text{/}}2$, не превосходит размерности компоненты ${{y}_{d}}({\kern 1pt} \cdot {\kern 1pt} )$, т.е. ${{n}_{k}} \leqslant {{n}_{y}}$.

Отметим, что, в отличие от условий, гарантирующих реконструкцию возмущения в частично наблюдаемой системе линейных СДУ (см. [12]), в случае квазилинейного СДУ можно потребовать наличие полного ранга матрицы $B_{1}^{m}$, но не матрицы $B_{1}^{d}$, которая зависит от решений систем (3.2), (3.3) и (3.8), (3.9).

Лемма 1. Пусть все структурные элементы исходного уравнения (2.1) (именно, матрицы $A$, $B$, $f$, вектор ${{x}_{0}}$, множества ${{S}_{{{{u}_{1}}}}}$ и ${{S}_{{{{u}_{2}}}}}$) таковы, что решение $m(t)$ уравнения (3.1) покоординатно строго больше нуля на всем промежутке времени $T$, а решение $d(t)$ уравнения (3.6) покоординатно неотрицательно на всем промежутке времени $T$. Тогда существует $\bar {N}$ такое, что $\forall N \geqslant \bar {N}$ матрица $B_{1}^{d}(\xi _{{di}}^{N},\xi _{{mi}}^{N})$ с вероятностью $1 - g(N)$ является матрицей полного ранга $\forall i \in [1:(l(N) - 1)]$.

Доказательство. Отметим, что условие леммы выполняется, например, в случае, когда все элементы матриц $A$, $B$, $f$ неотрицательны, все координаты начального вектора ${{m}_{0}}$ строго больше нуля (тогда как элементы матрицы ${{D}_{0}}$ неотрицательны вследствие предположения о независимости координат начального вектора ${{x}_{0}}$), и возмущения ${{u}_{1}}$ и ${{u}_{2}}$ принимают неотрицательные значения. Полагаем, что существует величина $\beta > 0$ такая, что $\forall t \in T$, $\forall j \in [1:n]$ ${{m}_{j}}(t) \geqslant \beta $.

Полнота ранга матрицы $B_{1}^{d}(\xi _{{di}}^{N},\xi _{{mi}}^{N})$ размерности ${{n}_{y}} \times {{n}_{{yz}}}$, ввиду ее структуры, описанной выше, определяется тем, что на местах, соответствующих ${{n}_{y}}$ единицам матрицы $C_{1}^{d}$, расположены ненулевые элементы, являющиеся элементами диагональной матрицы $\bar {B}(d({{\tau }_{i}}),\bar {m}({{\tau }_{i}}))$ с индексами, пробегающими множество ${{I}_{y}}$ (см. (3.7)), и при подстановке $\xi _{{di}}^{N}$ и $\xi _{{mi}}^{N}$. Фактически, для доказательства леммы следует проверить выполнение соотношений

Интересующие нас элементы матрицы $\bar {B}(d({{\tau }_{i}}),\bar {m}({{\tau }_{i}}))$ можно разделить на два типа, согласно (3.7): первый представляет собой сумму дисперсии измеряемой координаты и квадрата ее математического ожидания, т.е. второй – сумму ковариации двух измеряемых координат и произведения их математических ожиданий, т.е. Очевидно, что все ${{\bar {b}}_{1}}({{\tau }_{i}}) \geqslant {{\beta }^{2}}$ и ${{\bar {b}}_{2}}({{\tau }_{i}}) \geqslant {{\beta }^{2}}$, $i \in [1:(l(N) - 1)]$. Для доказательства утверждения леммы рассмотрим, к примеру, соотношение, получаемое для элементов второго типа (индекс $s$ соответствует паре $(j,p)$, см. (3.7), (4.1)): Данное неравенство получено из (3.10) и выполняется с вероятностью $1 - g(N)$. Здесь $\bar {K}$ – некоторая константа, которая может быть выписана явно. Отсюда выводим Полагаем, что последнее равенство выполняется $\forall N \geqslant \bar {N}$; $\bar {N}$ существует, поскольку $h(N) \to 0$ при $N \to \infty $. Соотношение (4.1) выполняется, следовательно, лемма доказана.

Выполнение условий 1, 2 и леммы 1 гарантирует единственность функций ${{u}_{1}}({\kern 1pt} \cdot {\kern 1pt} )$ и $\bar {u}({\kern 1pt} \cdot {\kern 1pt} )$, определяющих решения $({{y}_{m}}({\kern 1pt} \cdot {\kern 1pt} ),{{z}_{m}}({\kern 1pt} \cdot {\kern 1pt} ))$ и $({{y}_{d}}({\kern 1pt} \cdot {\kern 1pt} ),{{z}_{d}}({\kern 1pt} \cdot {\kern 1pt} ))$ систем (3.2), (3.3) и (3.8), (3.9) соответственно. Обоснование этого утверждения опирается на свойства псевдообратной матрицы для матрицы полного ранга (см. [4], [21]), т.е. в данной постановке на свойства псевдообратных матриц для $B_{1}^{m}$ и $B_{1}^{d}$, и, за вычетом малозначительных деталей, следует аналогичным рассуждениям из [4], [5].

Алгоритм, приведенный ниже, является приложением вычислительной процедуры из [4] к системам (3.2), (3.3) и (3.8), (3.9). В начальный момент ${{\tau }_{0}} = 0$ фиксируется значение $N$, определяются величины ${{l}^{N}} = l(N)$, ${{h}^{N}} = h(N)$ и ${{g}^{N}} = g(N)$ (см. (3.10)) и строится равномерное разбиение промежутка $T$ с шагом ${{\delta }^{N}} = \vartheta {\text{/}}{{l}^{N}}$: ${{\tau }_{i}} \in T,\;{{\tau }_{i}} = i{{\delta }^{N}},\;i \in [0:{{l}^{N}}]$. Вводится управляемая система-модель, фактически содержащая два блока, относящихся к системам (3.2), (3.3) и (3.8), (3.9). Фазовый вектор модели обозначим через $w(t)$; он состоит из двух троек: (i) $q$-мерного вектора ${{w}_{{my}}}(t)$, $(n - q)$-мерного вектора ${{w}_{{mz}}}(t)$, $q$-мерного вектора ${{w}_{{mu}}}(t)$ и (ii) ${{n}_{y}}$-мерного вектора ${{w}_{{dy}}}(t)$, ${{n}_{z}}$‑мерного вектора ${{w}_{{dz}}}(t)$, ${{n}_{y}}$-мерного вектора ${{w}_{{dv}}}(t)$. Динамика модели и ее начальное состояние определяются соотношениями

Здесь $B_{1}^{{m + }}$, $B_{1}^{{d + }}$ – псевдообратные матрицы, $\bar {u}_{{1i}}^{N}$, $\bar {v}_{i}^{N}$, $\hat {u}_{{1i}}^{N}$, $\hat {v}_{i}^{N}$ – управляющие воздействия соответствующих размерностей, вычисляемые позиционно в момент ${{\tau }_{i}}$ по правилам, которые конкретизируются ниже.

Динамика (4.2) выбирается из следующих соображений. Движение вспомогательных компонент ${{w}_{{my}}}(t)$ и ${{w}_{{dy}}}(t)$ при подходящем выборе модельных управлений

обеспечивает аппроксимацию координат ${{y}_{m}}(t)$ и ${{y}_{d}}(t)$. Тогда, пользуясь оценкой (3.10) и формальным выражением возмущений ${{u}_{1}}(t)$ и $\bar {u}(t)$ из уравнений (3.2) и (3.8) с заменой ${{\dot {y}}_{m}}(t)$ и ${{\dot {y}}_{d}}(t)$ на $\bar {u}_{{1i}}^{N}$ и $\bar {v}_{i}^{N}$ соответственно, ожидаем близость ${{w}_{{mz}}}(t)$ к ${{z}_{m}}(t)$ и ${{w}_{{dz}}}(t)$ к ${{z}_{d}}(t)$, что, в свою очередь, делает возможным отслеживание компонентами ${{w}_{{mu}}}(t)$ и ${{w}_{{dv}}}(t)$ координат ${{y}_{m}}(t)$ и ${{y}_{d}}(t)$ посредством выбора модельных управлений приближающих в нужном смысле функции ${{u}_{1}}(t)$ и $\bar {u}(t)$. Очевидно, нетрудно записать дискретный аналог модели (4.2).

Работа алгоритма разбивается на ${{l}^{N}}$ однотипных шагов. На $i$-м шаге, который выполняется на интервале $({{\tau }_{i}},{\kern 1pt} {{\tau }_{{i + 1}}}]$, исходными данными для вычислений служат оценки $\xi _{{mi}}^{N}$, $\xi _{{di}}^{N}$ и сформированное к этому моменту состояние модели $w({{\tau }_{i}})$. Предполагая покоординатную ограниченность правых частей уравнений (3.2) и (3.8) константой $\bar {K}$ (ее существование очевидно), находим $s$-ю координату $\bar {u}_{{1is}}^{N}$ вектора $\bar {u}_{{1i}}^{N}$ и $s$-ю координату $\bar {v}_{{is}}^{N}$ вектора $\bar {v}_{i}^{N}$ из соотношений

Вторая пара модельных управлений определяется следующим образом: $\hat {u}_{i}^{N}$ и $\hat {v}_{i}^{N}$ суть единственные решения экстремальных задач где $\langle \cdot , \cdot \rangle $ – скалярное произведение в соответствующем евклидовом пространстве, ${{\alpha }^{N}} = \alpha ({{h}^{N}})$ – параметр регуляризации. После вычисления управлений по формулам (4.3) и (4.4) состояние модели $w({{\tau }_{{i + 1}}})$ пересчитывается согласно дискретному аналогу (4.2). Процесс заканчивается в конечный момент времени $\vartheta $.

Теорема 1. Пусть выполняются условия согласования параметров

Тогда для модельных управлений $\hat {u}_{1}^{N}({\kern 1pt} \cdot {\kern 1pt} )$ и ${{\hat {v}}^{N}}({\kern 1pt} \cdot {\kern 1pt} )$, формируемых согласно (4.4), имеет место сходимость при $N \to \infty $: При дополнительных предположениях об ограниченности вариации исходных возмущений ${{u}_{1}}({\kern 1pt} \cdot {\kern 1pt} )$ и ${{u}_{2}}({\kern 1pt} \cdot {\kern 1pt} )$ справедлива следующая оценка точности алгоритма относительно количества реализаций процесса, доступных измерению: где ${{C}_{1}}$ и ${{C}_{2}}$ – некоторые константы, не зависящие от $N$, ${{u}_{1}}( \cdot )$ и $\bar {u}( \cdot )$.

Доказательство теоремы, за вычетом некоторых технических деталей, повторяет доказательство соответствующего утверждения для алгоритма реконструкции двух возмущений в случае измерения части координат линейного СДУ (см. [12]). В частности, там показано, что стандартные оценки $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}})$, $i \in [1:{{l}^{N}}],$ по следующим правилам (см. [18]):

с учетом измерения $q$ первых координат, обеспечивают выполнение свойства (2.2) (и, следовательно, (3.10)).

Теорема 2. Если искомый вектор ${{u}_{2}}( \cdot )$ является единственным, и его координаты неотрицательны, т.e. ${{S}_{{{{u}_{2}}}}}$ таково, что ${{u}_{{2i}}}(t) \geqslant 0$ $\forall t \in T$ $\forall i \in [1:n]$, то алгоритм (4.2)–(4.5) восстанавливает ${{u}_{2}}({\kern 1pt} \cdot {\kern 1pt} )$ в смысле ${{L}_{2}}(T;{{\mathbb{R}}^{n}})$-метрики.

Доказательство. Отметим, что ограничение ${{u}_{{2i}}}(t) \geqslant 0$ представляется естественным для шума с нулевым средним. Рассмотрим $n$-мерный вектор $\bar {\bar {u}}({\kern 1pt} \cdot {\kern 1pt} )$, состоящий из координат вектора $\bar {u}({\kern 1pt} \cdot {\kern 1pt} )$, которые равны $u_{{2i}}^{2}$, $i \in [1:n]$, и вектор ${{\bar {\bar {v}}}^{N}}({\kern 1pt} \cdot {\kern 1pt} )$, состоящий из соответствующих координат вектора ${{\hat {v}}^{N}}({\kern 1pt} \cdot {\kern 1pt} )$. Из структуры множества ${{S}_{{\bar {u}}}}$ и формул (4.4) для нахождения ${{\hat {v}}^{N}}({\kern 1pt} \cdot {\kern 1pt} )$ следует, что $\bar {\bar {v}}_{i}^{N}(t) \geqslant 0$ $\forall t \in T$ $\forall i \in [1:n]$. Поэтому возможно положить $v_{{2i}}^{N}(t) = \sqrt {\bar {\bar {v}}_{i}^{N}(t)} $. Тогда из (4.6) получаем

Используя неравенства

выводим Таким образом, имеем сходимость, соответствующую (4.6) и доказывающую аппроксимацию в ${{L}_{2}}(T;{{\mathbb{R}}^{n}})$-метрике возмущения ${{u}_{2}}({\kern 1pt} \cdot {\kern 1pt} )$ вектором $v_{2}^{N}({\kern 1pt} \cdot {\kern 1pt} )$. Порядок полученной при дополнительных предположениях оценки точности алгоритма относительно величины $1{\text{/}}N$ (см. (4.7)) сохраняется.

5. ЧИСЛЕННЫЙ ПРИМЕР

Перейдем к описанию модельного примера. Рассмотрим квазилинейное уравнение, описывающее случайный процесс, который можно трактовать как “испорченный” средневозвратный процесс Орнштейна–Уленбека (см. [17]):

Уравнения типа (5.1) используются, в частности, в некоторых простейших моделях, описывающих динамику численности относительно устойчивых популяций живых организмов (см. [17]). В таком случае величины ${{x}_{1}}(t)$ и ${{x}_{2}}(t)$ представляют текущие численности (в условных единицах) двух взаимодействующих популяций; структура детерминированной части уравнения определяет численности (например, средние за некоторую предысторию), стремление вернуться к которым “подсознательно” присутствует у популяции. Такому возврату могут препятствовать как влияние соседей, так и воздействие внешних факторов через функцию ${{u}_{1}}(t)$ в детерминированной части и функцию ${{u}_{2}}(t)$, входящую в амплитуду случайных колебаний. Именно внешние возмущения ${{u}_{1}}(t)$ и ${{u}_{2}}(t)$ подлежат восстановлению. Измерению в дискретные моменты времени доступны траектории координаты ${{x}_{1}}(t)$, соответствующие различным колониям организмов первого типа.

Запишем “редуцированные” системы вида (3.2), (3.3) и (3.8), (3.9), описывающие, соответственно, динамику математического ожидания $m(t)$ и ковариационной матрицы $D(t) = \left( {\begin{array}{*{20}{c}} {{{d}_{1}}}&{{{d}_{2}}} \\ {{{d}_{2}}}&{{{d}_{3}}} \end{array}} \right)$:

Фактически в дискретные моменты времени в системе (5.2) измеряется координата ${{m}_{1}}(t)$, в системе (5.3) – координата ${{d}_{1}}(t)$ (в соответствии с оценкой (3.10)), восстанавливаются величины ${{u}_{1}}(t)$ и ${{u}_{2}}(t)$. В данном случае $\bar {u}(t) = u_{2}^{2}(t)$, $\bar {u} \in [1,4]$, и выполняются условия реконструкции самой функции ${{u}_{2}}(t)$. Аналог модели (4.2) для систем (5.2), (5.3) выписывается естественным образом. Кроме того, условие 1 выполняется, так как $B_{1}^{m}(t) \equiv 1$, а справедливость леммы 1 для $B_{1}^{d}(\xi _{{di}}^{N},\xi _{{mi}}^{N}) = \xi _{{d1i}}^{N} + {{(\xi _{{m1i}}^{N})}^{2}}$ обеспечивается положительностью решений (5.2), (5.3) и достаточной точностью оценок (3.10). В вычислительном эксперименте были выбраны следующие реализации неизвестных функций:

Для моделирования динамики системы (5.1) использовался метод Эйлера с заменой винеровского процесса последовательностью случайных импульсов (подробное описание аппроксимационной схемы можно найти, например, в [22]); его среднеквадратический порядок точности для квазилинейного СДУ равен $O({{\delta }_{s}})$, где ${{\delta }_{s}}$ – шаг моделирования. Этот шаг, как правило, очень мал относительно шага ${{\delta }^{N}}$, с которым поступает информация о траекториях системы, обрабатываемая согласно (4.8), и, соответственно, работает процедура восстановления неизвестных функций ${{u}_{1}}(t)$ и ${{u}_{2}}(t)$. Для получения оценок (3.10) необходимо отслеживание $N$ независимых траекторий системы (5.1).

Результаты восстановления функций ${{u}_{1}}(t)$ и ${{u}_{2}}(t)$, полученные для различных наборов параметров алгоритма (4.2)–(4.5), приведены на фиг. 1 и фиг. 2, где реальные функции ${{u}_{1}}(t)$ и ${{u}_{2}}(t)$ изображены штриховой линией, а результаты их восстановления – сплошной. Они согласуются с основным утверждением статьи, подтверждая сходимость типа (4.6)–(4.7): чем больше количество $N$ измеряемых траекторий, тем меньше (с вероятностью, близкой к 1) погрешность восстановления.

6. ЗАКЛЮЧЕНИЕ

Таким образом, исследована постановка обратной задачи для квазилинейного СДУ, предполагающая динамическую реконструкцию двух неизвестных неслучайных возмущений, входящих в детерминированный и стохастический члены уравнения. В качестве входной информации используются точные измерения некоторого количества реализаций части координат случайного процесса в дискретные моменты времени. Решение задачи состоит в переходе к обратной задаче для двух систем ОДУ специального вида с последующим использованием модификации алгоритма восстановления, разработанного ранее для частично наблюдаемых систем ОДУ. Доказанная оценка точности алгоритма относительно количества доступных измерению реализаций проиллюстрирована на модельном примере. Ряд вопросов, связанных с рассматриваемой постановкой, остается открытым: например, об улучшении порядка точности оценки (4.7) при дополнительных предположениях и о возможности введения погрешности измерения траекторий исходного СДУ.