Известия РАН. Теория и системы управления, 2021, № 5, стр. 52-73

Введение. Рассматриваемая в работе задача может быть отнесена к традиционной задаче оптимального управления дифференциальными стохастическими системами. Задача решается классическими методами, равно как и основанными на них приближенными методами поиска оптимальных управлений [1–3]. Ключевые результаты, составляющие классическую теорию управления, получены более 40 лет назад. Наибольшую полноту и результативность имеет задача управления линейно-гауссовскими стохастическими системами по квадратичному критерию (LQG) [4]. Эта задача дает исходный базис, на котором основано множество успешно исследованных и решенных задач стохастического управления и оптимизации. Множество результатов в этой области как в части моделей, так и в части математического аппарата, имевших место в последующие годы, существенно обогатили теорию управления. Но и до настоящего времени линейные модели и квадратичный критерий, несмотря на всю справедливую критику в отношении их адекватности и гибкости, сохраняют исследовательский интерес и находят современные области приложения.

Не претендуя на сколько-нибудь полное обоснование последнего тезиса, приведем несколько примеров, показавшихся наиболее интересными. Так, в [5] решается линейно-квадратичная задача в игровой постановке с запаздыванием. В близкой по модели работе [6] задача управления ставится в терминах ${{H}_{\infty }}$-робастности. Точнее называть эту тематику H₂/H_∞-управлением, и работ по этой теме очень много. Аккуратности ради следует уточнить, что под линейными в этой теории понимаются модели с мультипликативными по состоянию возмущениями. Совсем другой класс моделей, особо популярных в последние годы, составляют скачкообразные процессы. Например, линейные уравнения в сочетании с пуассоновскими скачками в [7] используются в моделях, описывающих различные показатели функционирования сетевых протоколов передачи данных транспортного уровня. Надо заметить, что телекоммуникации представляют в последние годы популярный прикладной материал для исследований, работ по этой проблематике множество, математические техники привлекаются самые разные и самые современные, но и линейным моделям место находится. В этой сравнительно новой прикладной области имеются приложения и для методов стохастического оптимального управления [8].

Близкой к области телекоммуникаций являются задачи управления ресурсами программных систем, в качестве которых могут выступать сайты, базы данных и банки знаний, а также традиционные вычислительные ресурсы, такие, как узлы распределенных систем, компоненты центров обработки данных и т.д. Такие ресурсы управляются программами – менеджерами задач, оптимизаторами запросов в системах управления базами данных, системами балансировки нагрузки. Моделями для постановки задач оптимизации функционирования программных систем являются системы массового обслуживания или, при значительном росте числа состояний, их диффузионные аппроксимации [9]. Соответствующая практическая интерпретация в качестве состояния оптимизируемой системы может использовать число пользователей, запросов, буферов памяти, процессорного времени и т.п., описывая его диффузионным процессом. Объем ресурсов при этом – это линейный выход, учитывающий и состояние, и активное действие программной системы, такое, как объем выделенных или освобожденных ресурсов. Отметим, что такая прикладная интерпретация изучаемой задачи ранее была представлена для случая дискретного времени [10, 11].

Исследуемая в статье модель следует традиционному представлению управляемой системы в форме вход-выход. Входная переменная предполагается неуправляемой, оптимизируется выход. Такого рода задачи известны в числе многих исследуемых постановок задач оптимального управления под терминами управление (регулирование) выходом или управление по выходу. Одно из первых применений этой терминологии обнаруживается в [12], где специальным образом записанная классическая LQG-задача используется для математической формализации инженерной задачи “доведения выхода до нуля”.

В качестве примера более современных постановок в терминологии управления выходом можно привести работу [13] (также указан ряд релевантных ссылок аналогичной тематики). Неопределенность там не имеет стохастической природы, а описывается простыми геометрическими ограничениями, что указывает на робастный характер применяемого метода инвариантных эллипсоидов. Изучаемая система при этом является линейной и автономной. Аналогичную модель исследуют авторы [14] с помощью аппарата линейных матричных неравенств. В указанной и других работах оптимизируется квадратичный критерий и используются свойства линейных регуляторов. Исследований, посвященных нелинейным моделям систем с управляемым выходом, обнаруживается немного. Так, довольно популярны постановки в терминах нелинейного адаптивного управления, предлагающие различные варианты задачи стабилизации [15].

В данной работе, как и во многих других, эксплуатируются привлекательные свойства линейных моделей и квадратичного критерия. Обобщение затронуло модель динамики системы. Оптимизируемая динамическая система описывается двумя системами уравнений. Состояние задается системой нелинейных стохастических дифференциальных уравнений Ито, не содержащих управляемой переменной. Возмущения здесь описываются стандартным векторным винеровским процессом, накладываются простые условия существования и единственности решения. Поскольку состояние не управляется, то возможно его интерпретировать еще и как сложное внешнее возмущение. Вторая переменная – управляемый выход – задается системой линейных стохастических дифференциальных уравнений. Цель оптимизации выхода формируется квадратичным критерием общего вида. Формальная постановка задачи дана в разд. 2. Для решения задачи применяется метод динамического программирования, решается уравнение Беллмана [3], используемые общие положения приведены в разд. 1. Соответственно в результате получаются аналитические выражения и для оптимального управления, и для значения функционала качества. Технически традиционный, стандартный подход к задаче обременен, пожалуй, единственной проблемой – поиском верного представления структуры функции Беллмана. Справиться с этой проблемой в большей степени удается за счет результата, полученного при решении дискретного по времени аналога рассматриваемой постановки [11]. Конечные соотношения для оптимального решения, как и во всех подобных задачах, включая классическую LQG, содержат решения определенных дифференциальных уравнений (обыкновенных и в частных производных). Вывод этих уравнений представлен в разд. 3. Решение дает функция Беллмана вида $V(t,y,z) = {{z}^{T}}\alpha (t)z + {{z}^{T}}\beta (t,y) + \gamma (t,y)$, в котором переменная y соответствует вектору состояния, переменная z – выходу, а коэффициенты α(t), $\beta \left( {t,y} \right),$ $\gamma \left( {t,y} \right)$ представлены решениями полученных дифференциальных уравнений (обыкновенного Риккати и системы линейных уравнений в частных производных параболического типа). При определении функции Беллмана получено и выражение для оптимального управления в виде линейной функции z и коэффициента $\beta \left( {t,y} \right).$ Соответственно естественная численная реализация оптимального решения – это приближенное решение методом конечных разностей соответствующих дифференциальных уравнений. При этом применение сеточных методов для параболических уравнений осложняется тем обстоятельством, что имеются только начальные условия для уравнений, определяющих $\beta (t,y),$ $\gamma (t,y),$ оснований для указания граничных условий, нужных для метода конечных разностей, нет. Из-за этого граничные условия в модельном эксперименте, которому посвящен раздел 6, выбираются довольно волюнтаристски, обсуждается экспериментальный анализ их влияния на решение. И хотя удается получить удовлетворительный результат, поиск альтернативного метода для приближенного вычисления коэффициентов функции Беллмана сохраняет актуальность.

Поскольку обсуждаемые коэффициенты описываются линейными уравнениями в частных производных второго порядка параболического типа, то перспективным представлялось опираться в исследовании этой задачи на известные факты, ассоциирующие такие уравнения с уравнениями стохастическими. Видимо самая известная ассоциация такого рода – это уравнение А.Н. Колмогорова [16], которое и было использовано в качестве инструмента для реализации альтернативного метода приближенного вычисления коэффициентов $\beta \left( {t,y} \right),$ $\gamma \left( {t,y} \right).$ Точнее применялось утверждение, которое определяет связь между решением задачи Коши для уравнения А.Н. Колмогорова, с одной стороны, и теоретико-вероятностным представлением ее терминального условия. Использованное в работе интегральное представление этого результата известно как формула Фейнмана–Каца [17]. Такой метод для численного решения параболических уравнений и уравнений даже более общей структуры успешно реализован и исследован [18], так что примененный в статье подход можно считать апробированным. Обсуждению этого результата посвящен раздел 4. Соответственно стохастическая интерпретация $\beta \left( {t,y} \right)$ дает возможность определения этого коэффициента как статистической оценки математического ожидания терминальной переменной вспомогательной динамической системы, формируемой из интегрального представления Фейнмана–Каца. Таким образом, альтернативный метод аппроксимации оптимального решения рассматриваемой задачи сводится к использованию компьютерного моделирования и метода Монте-Карло. Результаты применения этого метода и сравнение с методом конечных разностей приведены в разд. 6, представляющем численный эксперимент.

Практически важный частный случай рассматриваемой задачи, случай линейного сноса, но нелинейной диффузии, обычно характеризуемый как модель с мультипликативными возмущениями, рассмотрен в разд. 5. Оптимальное решение в этом случае не требует интегрирования параболического уравнения, управление линейно по состоянию и выходу, а численная реализация сводится к решению только обыкновенных дифференциальных уравнений. Это дает возможность, во-первых, получить подходящий эталонный пример для сравнительного анализа двух других методов аппроксимации оптимального решения – сеточного и имитационного. Во-вторых, позволяет в полном объеме рассмотреть модельный пример, связанный с выбранной целевой прикладной областью телекоммуникаций. Именно в разд. 6 обсуждаются результаты численного эксперимента, выполненного предложенными методами для модели Кокса–Ингерсона–Росса, описывающей показатель времени сетевого отклика RTT (round-trip time) основного сетевого протокола Интернета TCP (transmission control protocol).

1. Общий подход к решению. Предваряя основной результат, кратко опишем положения, привлекаемые далее к решению задачи, следуя обозначениям и терминологии [3]. Будем рассматривать задачу оптимального управления состоянием $x\left( t \right) \in {{\mathbb{R}}^{n}}$ стохастической динамической системы по полной информации, применяя метод динамического программирования. Допустимыми будем считать управления с обратной связью, т.е. функции вида $u = u\left( {t,x} \right) \in {{\mathbb{R}}^{k}},$ $x \in {{\mathbb{R}}^{n}}$, множество допустимых управлений обозначается $U_{0}^{T} = \{ u\left( {t,x} \right),{\text{\;}}0 \leqslant t \leqslant T,~~x \in {{\mathbb{R}}^{n}}\} $, такие, что реализациями управления $u(t) = u(t,x(t))$ обеспечиваются выполненными условия существования x(t) для $u \in U_{0}^{T}$. Заметим, что использование обозначений $u = u\left( {t,x} \right)$ для управления с обратной связью и $u(t) = u(t,x(t))$ для его реализации, соответствующей траектории x(t), далее будет везде понятно из контекста и не приведет к неоднозначности.

Минимизируемый целевой функционал имеет вид

где $\mathbb{E}\left\{ \cdot \right\}$ – оператор математического ожидания, x(t) – случайный процесс, описываемый системой стохастических дифференциальных уравнений Ито: где $W\left( t \right) \in {{\mathbb{R}}^{m}}$ – стандартный векторный винеровский процесс, $X \in {{\mathbb{R}}^{n}}$ – случайный вектор с конечным вторым моментом, $m\left( {t,x,u} \right),~\lambda \left( {t,x} \right)$ – векторная и матричная функции соответствующих размерностей.

Относительно функций сноса $m\left( {t,x,u} \right)$ и диффузии ${{\lambda }}\left( {t,x} \right)$ в (1.2) будем предполагать выполненными какие-либо условия существования и единственности решения для заданного допустимого управления. Так, учитывая, что допустимыми являются управления с обратной связью, будем считать, что $m\left( {t,x,u\left( {t,x} \right)} \right)$ и $\lambda \left( {t,x} \right)$ удовлетворяют условию линейного роста по x и локальному условию Липшица по x равномерно по t (т.е. условиям Ито). Также относительно функций из целевого функционала $L\left( {t,x,u} \right)$ и l(x) предполагаются выполненными типовые условия существования, например ограничения полиномиального роста по x и u.

Для поиска оптимального управления, минимизирующего $J(U_{0}^{T})$, рассматривается функция Беллмана

где $\mathcal{F}_{t}^{x}$ – σ-алгебра, порожденная $x\left( \tau \right),{\text{\;}}0 \leqslant \tau \leqslant t$; $\mathbb{E}\{ \cdot |\mathcal{F}\} $ – оператор условного математического ожидания относительно $\mathcal{F}$; $V\left( {t,x} \right)$ предполагается непрерывно дифференцируемой по t и дважды непрерывно дифференцируемой по x. В качестве достаточного условия оптимальности рассматривается уравнение динамического программирования (используется представление из [3], записанное в векторной форме): где T – операция транспонирования, tr – след матрицы, ${{\partial }^{2}}V{\text{/}}\partial {{x}^{2}}$ – матрица Гессе V, $\partial V{\text{/}}\partial x$ – градиент V, $m = m\left( {t,x,u} \right)$, $\lambda = \lambda \left( {t,x} \right),$ $V = V\left( {t,x} \right)$. Далее также будут опускаться аргументы функций, если это возможно сделать однозначно.

Традиционно в рамках применения метода динамического программирования будем предполагать, что функции $L,l,m,\lambda $ обеспечивают существование хотя бы одного решения уравнения (1.4), а следовательно, и оптимального управления

и доставляющего минимум целевому функционалу (1.1). Задача оптимизации далее решается путем указания конкретных выражений для $L,l,m,\lambda $.

2. Постановка задачи управления выходом. Рассмотрим стохастическую дифференциальную систему, состояние которой представляет диффузионный процесс $y\left( t \right) \in {{\mathbb{R}}^{{{{n}_{y}}}}}$, описываемый системой нелинейных стохастических дифференциальных уравнений Ито

где ${v}(t) \in {{\mathbb{R}}^{{{{n}_{v}}}}}$ – стандартный векторный винеровский процесс, $Y \in {{\mathbb{R}}^{{{{n}_{y}}}}}$ – случайная величина с конечным вторым моментом, векторная функция $A\left( {t,y} \right):{{\mathbb{R}}^{{{{n}_{y}}}}} \to {{\mathbb{R}}^{{{{n}_{y}}}}}$ и матричная функция Σ(t, y) : ${{\mathbb{R}}^{{{{n}_{y}}}}} \to {{\mathbb{R}}^{{{{n}_{y}} \times {{n}_{{v}}}}}}$ удовлетворяют условиям Ито: обеспечивающим существование единственного потраекторного решения уравнения [3] (здесь и далее |⋅| обозначает евклидову норму вектора или матрицы).

Будем считать, что y(t) описывает состояние некоторой динамической системы. Под выходом этой системы понимается процесс $z\left( t \right) \in {{\mathbb{R}}^{{{{n}_{z}}}}}$, линейно связанный с состоянием:

где $w(t) \in {{\mathbb{R}}^{{{{n}_{w}}}}}$ – не зависящий от ${v}(t),Y,Z$ стандартный векторный винеровский процесс; $Z \in {{\mathbb{R}}^{{{{n}_{z}}}}}$ – случайная величина с конечным вторым моментом, не зависящая от ${v}(t),Y;$ $u(t) \in {{\mathbb{R}}^{{{{n}_{u}}}}}$ – управление. Задача формулируется в предположении наличия полной информации о состоянии y(t) и выходе z(t) (соответствующую σ-алгебру обозначим $\mathcal{F}_{t}^{{y,z}}$) и управление u(t) $\mathcal{F}_{t}^{{y,z}}$ – измеримо. Класс $U_{0}^{T}$ допустимых управлений составляют управления с полной обратной связью, т.е. функции вида $u = u\left( {t,y,z} \right) \in {{\mathbb{R}}^{{{{n}_{u}}}}},$ $y \in {{\mathbb{R}}^{{{{n}_{y}}}}}$ $z \in {{\mathbb{R}}^{{{{n}_{z}}}}}$, в предположении, что соответствующая реализация $u(t) = u(t,y(t),z(t))$ обеспечивает выполнение условий существования $y\left( t \right),z\left( t \right)$ для $u \in U_{0}^{T}$.

Качество управления $U_{0}^{T}$ определяется целевым функционалом следующего вида:

где $P\left( t \right) \in {{\mathbb{R}}^{{{{n}_{J}} \times {{n}_{y}}}}}$, $Q\left( t \right) \in {{\mathbb{R}}^{{{{n}_{J}} \times {{n}_{z}}}}}$, $R\left( t \right) \in {{\mathbb{R}}^{{{{n}_{J}} \times {{n}_{u}}}}}$, $S\left( t \right) \in {{\mathbb{R}}^{{{{n}_{J}} \times {{n}_{J}}}}}$, $S\left( t \right) \geqslant 0,$ $S\left( t \right) = {{S}^{{\text{T}}}}\left( t \right),$ $0 \leqslant t \leqslant T,$ – заданные ограниченные матричные функции, весовая функция ${\text{||}}x{\text{||}}_{S}^{2} = {{x}^{{\text{T}}}}Sx$ для симметричной неотрицательно определенной матрицы S, единичной матрице S = E соответствует евклидова норма ${\text{||}}x{\text{||}}_{E}^{2} = {{\left| x \right|}^{2}}.$ Для исключения возможности отсутствия штрафа для отдельных компонентов вектора управления u(t), которая делает целевой функционал (2.3) физически некорректным, предполагается выполненным обычное условие невырожденности, в данных обозначениях принимающее вид $R{{\left( t \right)}^{{\text{T}}}}S\left( t \right)R\left( t \right) > 0.$

Матричные функции $a\left( t \right) \in {{\mathbb{R}}^{{{{n}_{z}} \times {{n}_{y}}}}},$ $b\left( t \right) \in {{\mathbb{R}}^{{{{n}_{z}} \times {{n}_{z}}}}},$ $c\left( t \right) \in {{\mathbb{R}}^{{{{n}_{z}} \times {{n}_{u}}}}},$ $\sigma \left( t \right) \in {{\mathbb{R}}^{{{{n}_{z}} \times {{n}_{w}}}}}$ ограничены: |a(t)| + $\left| {b(t)} \right| + \left| {c(t)} \right| + \left| {\sigma (t)} \right| \leqslant C$ для всех $0 \leqslant t \leqslant T$. Таким образом, обеспечивается существование решения уравнения (2.2) для любого допустимого управления u. Более того, будем предполагать, что все используемые функции времени $a\left( t \right),b\left( t \right),c\left( t \right),\sigma \left( t \right),P\left( t \right),Q\left( t \right),R\left( t \right),S\left( t \right)$ кусочно-непрерывны, чтобы обеспечить выполнение типовых условий существования решений дифференциальных уравнений, получаемых далее.

Задачу составляет поиск $(U{\text{*}})_{0}^{T} = \{ u{\text{*}}\left( {t,y,z} \right),{\text{\;}}0 \leqslant t \leqslant T,y \in {{\mathbb{R}}^{{{{n}_{y}}}}},z \in {{\mathbb{R}}^{{{{n}_{z}}}}}\} $ – допустимого управления с обратной связью, доставляющего минимум квадратичному функционалу $J(U_{0}^{T}).$

Функционал (2.3) формализует типичную задачу управления выходом. В [12], как уже упоминалось, в качестве примера практического приложения такого типа целевой функции рассматривается задача “доведения выхода до нуля”. Предлагаемый более общий вид целевого функционала отражает практическое содержание упомянутой выше задачи распределения ресурсов. Так, частный вид функционала (2.3)

позволяет ставить задачи слежения выходной переменной за состоянием, учитывая расходы на управляющее воздействие и/или значение выходной переменной, функционал позволяет ставить задачи слежения управлением за выходной переменной, учитывая расходы на управляющее воздействие.

Поставленная задача очевидным образом формулируется в терминах выражений (1.1)–(1.3), для этого нужно обозначить

для записи уравнения состояния в виде (1.2) и для записи целевого функционала в виде (1.1). В записи блочной матрицы выше через ${{0}_{{n \times m}}}$ обозначена нулевая матрица размерности n × m.

Функция Беллмана (1.3) принимает вид $V = V\left( {t,y,z} \right)$, а уравнение (1.4) с обозначениями $\Sigma = \Sigma (t,y),$ $A = A\left( {t,y} \right)$, $a = a\left( t \right),$ $b = b\left( t \right),$ $c = c\left( t \right),$ $u = u\left( t \right),$ $\sigma = \sigma \left( t \right),$ $P = P\left( t \right),$ $Q = Q\left( t \right),$ $R = R\left( t \right),$ $S = S(t),$ записывается как

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

3. Динамическое программирование и оптимальное управление. В рассматриваемой постановке линейность выхода и квадратичность критерия дают те же преимущества, что и в классической LQG-задаче [4], а именно позволяют сразу определить вид оптимального управления и фактические условия его существования. Действительно, сохраняя в (2.4) под знаком $\mathop {{\text{min}}}\limits_u $ только члены, зависящие от u, получаем

Далее из (3.1) с учетом сделанного предположения ${{R}^{{\text{T}}}}SR > 0$, обозначения “A^–1” для обратной матрицы A и выражения для минимума квадратичной формы получаем равенство для оптимального управления в предположении существования решения (3.1):

Это управление доставляет минимум последнему слагаемому в (3.1):

Отметим, что оптимальное управление получается из двух слагаемых: одно слагаемое линейное по состоянию и выходу (как в классической линейно-квадратичной задаче), второе слагаемое нелинейное, определяется производной $\partial V{\text{/}}\partial z$.

Таким образом, функция Беллмана описывается следующим дифференциальным уравнением:

Преобразуя последнее слагаемое в (3.3), перепишем его в виде

Рассматривая полученное уравнение, предположим, что его решение может быть представлено в как

т.е. решением является функция, квадратичная по переменной z, и, таким образом, поиск оптимального решения сводится к поиску уравнений относительно симметричной матричной функции $\alpha = \alpha (t) \in {{\mathbb{R}}^{{{{n}_{z}} \times {{n}_{z}}}}}$, α = α^T, векторной функции $\beta = \beta \left( {t,y} \right) \in {{\mathbb{R}}^{{{{n}_{z}}}}}$ и скалярной функции γ = = $\gamma (t,y) \in {{\mathbb{R}}^{1}}.$ Отметим сразу, что явный вид функции γ для реализации оптимального управления не требуется, однако далее будет предложен вариант вычисления и этой функции. Источником для предложения (3.5) является уже упоминавшаяся аналогичная задача для случая дискретного времени [10, 11]. В той задаче выражение для функции Беллмана получается формально без дополнительных усилий. При этом форма (3.5) обнаруживается как свойство оптимального решения. В рассматриваемом случае непрерывного времени (3.5) является предположением, правильность которого подтверждается далее результирующими уравнениями для $\alpha ,\beta ,\gamma .$ Кроме того, данное предположение представляется вытекающим из линейной структуры задачи в отношении переменной z, в частности тем фактом, что такой вид функции Беллмана обеспечивает линейность оптимального управления (3.2) по z.

Граничное условие при сделанном предположении (3.5) принимает вид

т.е.

При этом поскольку (3.5) дает $\partial V{\text{/}}\partial z = 2\alpha z + \beta $, то само оптимальное управление, определенное выражением (3.2),

Дальнейшие преобразования удобнее выполнить, представив V в виде

обозначая $z = {{({{z}^{{\left( 1 \right)}}}, \ldots ,{\text{\;}}{{z}^{{\left( {{{n}_{z}}} \right)}}})}^{{\text{T}}}},$ $\beta = {{({{\beta }^{{\left( 1 \right)}}}, \ldots ,{{\beta }^{{\left( {{{n}_{z}}} \right)}}})}^{{\text{T}}}}$, ${{\beta }^{{\left( i \right)}}} = {{\beta }^{{\left( i \right)}}}\left( {t,y} \right)$, чтобы не сталкиваться при дифференцировании с необходимостью использования тензорных обозначений для производных векторной функции. Теперь следует подставить V, представленную в виде (3.5), в уравнение (3.4):

Раскрывая скобки и формируя слагаемые при ${{z}^{T}}z,~z$ (каждой компоненте ${{z}^{{\left( i \right)}}}$) и ${{z}^{0}} = 1$, имеем

Окончательно требуемая группировка слагаемых записывается с использованием обозначения ${{\left[ A \right]}^{{\left( {ji} \right)}}}$ для элемента  j-й строки i-го столбца матрицы A:

Поскольку полученное равенство выполняется для всех $z \in {{\mathbb{R}}^{{{{n}_{z}}}}}$, $y \in {{\mathbb{R}}^{{{{n}_{y}}}}}$, то каждое из трех слагаемых должно быть равно 0, что дает следующие уравнения:

− матричное дифференциальное уравнение для функции $\alpha = \alpha (t)$ (с учетом предположения $\alpha = {{\alpha }^{{\text{T}}}}$):

− систему дифференциальных уравнений в частных производных для функций β⁽ⁱ⁾ = = ${{\beta }^{{(i)}}}(t,y),~i = \overline {1,{{n}_{z}}} $:

− дифференциальное уравнение в частных производных для функции $\gamma = \gamma \left( {t,y} \right)$:

Вводя дополнительно обозначения

для $\beta \left( {t,y} \right)$ и $\gamma \left( {t,y} \right)$ вместо (3.9) и (3.10) окончательно получаем

Уравнение (3.8) – это матричное уравнение Риккати, точнее его частный случай для квадратной симметричной матрицы α. Сделанных выше предположений достаточно для существования единственного неотрицательного решения для всех $0 \leqslant t \leqslant T$. Этот факт требует дополнительного комментария. Рассмотрим исходную задачу при дополнительных условиях: $a\left( t \right) = {{0}_{{{{n}_{z}} \times {{n}_{v}}}}},$ $P\left( t \right) = {{0}_{{{{n}_{J}} \times {{n}_{y}}}}},$ $S\left( t \right) = E$ для всех $0 \leqslant t \leqslant T$. Это дает

и сводит рассматриваемую постановку к классической LQG-задаче. Для такого упрощения целевого функционала требуется наличие блочной структуры матриц Q, R из (2.3), чтобы исключить в упрощенном $J(U_{0}^{T})$ слагаемое с множителями $z$ и $u$. В этой упрощенной форме матрицы $Q\left( t \right) \in {{\mathbb{R}}^{{{{n}_{z}} \times {{n}_{z}}}}}$ и $R\left( t \right) \in {{\mathbb{R}}^{{{{n}_{u}} \times {{n}_{u}}}}}$ представляют соответствующие блоки одноименных матриц из (2.3).

Теперь уравнение (3.8) для α(t) принимает хорошо известный вид

В таком случае, как известно [19], существует единственное оптимальное управление – линейное с обратной связью по выходу z(t), с коэффициентом усиления, описываемым уравнением Риккати (3.13). Именно этот результат дают уравнения (3.8)–(3.10) и описываемая ими функция Беллмана (3.5), так как из $a\left( t \right) = 0,P\left( t \right) = 0$ немедленно следует, что ${{\beta }^{{\left( i \right)}}}\left( {t,y} \right) = 0,{\text{\;}}i = \overline {1,{{n}_{z}}} ,$ откуда, в свою очередь, с учетом независимости решения от y следует, что γ не зависит от y и задается уравнением $d{{\gamma /}}dt + {\text{tr}}\{ {{\sigma }^{{\text{T}}}}\alpha \sigma \} = 0,{{\gamma }}\left( T \right) = 0$. Оптимальное управление при этом $u{\text{*}} = - {{({{R}^{{\text{T}}}}R)}^{{ - 1}}}{{c}^{{\text{T}}}}\alpha z,$ т.е. все полностью совпадает с известным классическим решением.

С уравнениями (3.11) и (3.12) ситуация, естественно, обстоит сложнее. Это линейные уравнения в частных производных второго порядка параболического типа, поскольку ${{\Sigma }^{{\text{T}}}}\Sigma \geqslant 0$. Уравнения неоднородные, причем в уравнении (3.12) для $\gamma \left( {t,y} \right)$ участвуют все ${{\beta }^{{\left( i \right)}}}\left( {t,y} \right)$, т.е. сначала должны быть решены уравнения (3.11) и результат подставлен в (3.12). Также в уравнения (3.11) нужно подставить α(t) – решение уравнения (3.8). Далее, в каждом из уравнений (3.11) неоднородность для ${{\beta }^{{\left( i \right)}}}\left( {t,y} \right)$ формируется неизвестными функциями ${{\beta }^{{\left( j \right)}}}\left( {t,y} \right),{\text{\;}}j \ne i$, т.е. параболические уравнения объединены в систему и по отдельности решаться не могут, так что актуален еще и вопрос их приближенного решения. Кроме того, отсутствуют конструктивные условия, гарантирующие существование решений (требовать, чтобы все фигурирующие в уравнениях коэффициенты были представлены аналитическими функциями на всем пространстве значений, вряд ли целесообразно, а в имеющейся системе уравнений, может быть, и не достаточно). Поэтому далее будем предполагать, что уравнения (3.11) и (3.12) имеют на интервале $0 \leqslant t \leqslant T$ хотя бы одно гладкое (непрерывно дифференцируемое по t и дважды непрерывно дифференцируемое по y) решение и именно эти условия можно считать достаточными условиями существования оптимального решения рассматриваемой задачи.

Таким образом, доказана следующая теорема.

Теорема. Пусть для диффузионного процесса (2.1) выполнены условия Ито, для выхода (2.2) ограничены коэффициенты, уравнения (3.8)–(3.10) имеют гладкие решения для $0 \leqslant t \leqslant T$. Тогда минимум функционалу (2.3) доставляет $(U{\text{*}})_{0}^{T}$ – оптимальное управление с полной обратной связью, определенное (3.7).

4. Стохастическое представление коэффициентов функции Беллмана. Для реализации оптимального управления и оценки его качества потребуется решать полученные уравнения (3.8), (3.11), (3.12) приближенно. Решение уравнения (3.8) не представляет трудностей и может быть выполнено любым из множества эффективных численных методов решения обыкновенных дифференциальных уравнений. Несколько сложнее обстоит ситуация с уравнениями (3.11), (3.12). Для уравнений параболического типа также имеются эффективные численные схемы, основанные на методе конечных разностей. Но в рассматриваемом случае трудности в реализации возникают при определении граничных условий, необходимых для этого метода. В рассматриваемой задаче физические основания для определения границы и, тем более граничных условий как для $\beta \left( {t,y} \right),$ так и для $\gamma \left( {t,y} \right)$ отсутствуют. Вариант возможного преодоления этих сложностей обсуждается далее в рамках представления численного эксперимента. Более подходящим и лишенным таких трудностей представляется решение параболического уравнения на основе его ассоциации с уравнением А.Н. Колмогорова, точнее связи между решением задачи Коши с терминальным условием для уравнения А.Н. Колмогорова и теоретико-вероятностным представлением этого терминального условия [16]. Выражением этого результата является формула Фейнмана–Каца [17], непосредственно применимая к уравнению (3.12). С учетом введенных выше обозначений интегральное представление Фейнмана–Каца решения дифференциального уравнения (3.12) имеет вид

В предположении, что α(t) и $\beta \left( {t,y} \right)$ известны (вычислены), (4.1) дает основание для применения компьютерного моделирования для аппроксимации $\gamma \left( {t,y} \right)$. Многократно моделируя решение уравнения (2.1), проходящее через точку (t, y) (именно это обозначает σ-алгебра $\mathcal{F}_{t}^{y}$ в условном математическом ожидании), значение $\gamma \left( {t,y} \right)$ можно получить как статистическую оценку функции терминального условия, согласно (4.1). Подробнее этот метод описан далее. Следует отметить, что как в формулировке для уравнения А.Н. Колмогорова, так и в формуле Фейнмана–Каца утверждается существование процесса y_t, обеспечивающего равенство (4.1), т.е. формально это другой процесс и задан он на другом пространстве, чем процесс, участвующий в поставленной выше задаче управления. При этом описывается этот процесс таким же уравнением (2.1) и имеет те же вероятностные характеристики, что дает основание использовать уже имеющуюся дифференциальную систему, не вводить новые обозначения.

Такой же подход хотелось бы применить для приближенного вычисления $\beta \left( {t,y} \right).$ Более того, в скалярном случае ${{n}_{z}} = 1$ интегральное представление Фейнмана–Каца решения дифференциального уравнения (3.11) записывается без дополнительных усилий и имеет вид

Далее выводится аналог (4.2) для ${{n}_{z}} > 1.$ Для этого введем новую переменную $\bar {\beta } = \bar {\beta }\left( {t,y} \right) \in {{\mathbb{R}}^{{{{n}_{z}}}}}$ и выполним замену $\beta \left( {t,y} \right) = I\left( t \right)\bar {\beta }\left( {t,y} \right),$ где матрица $I = I\left( t \right) \in {{\mathbb{R}}^{{{{n}_{z}} \times {{n}_{z}}}}}$ является решением уравнения $dI(t){\text{/}}dt + {{N}^{{\text{T}}}}(t)I(t) = 0$. Такая замена позволяет исключить из уравнений (3.11) слагаемые с ${{\beta }^{{\left( j \right)}}}$. Действительно, для слагаемых

записанных в векторной форме $\frac{{\partial \beta }}{{\partial t}} + {{N}^{{\text{T}}}}\beta $, замена переменной на $\bar {\beta }$ дает $I{\text{\;}}\partial{ \bar {\beta }}{\text{/}}\partial t + dI{\text{/}}dt~\bar {\beta } + {{N}^{{\text{T}}}}I\bar {\beta }$ = = $I{\text{\;}}\partial{ \bar {\beta }}{\text{/}}\partial t$. Необходимая для такой замены матрица I(t) существует и равна матричной экспоненте причем I(t) невырожденная [20], что позволяет далее использовать обратную матрицу I^–1(t). Остальные слагаемые в (3.11) с дифференцированием по $t$ не связаны и линейны по β, что дает

Следующим действием левая часть полученного равенства разрешается относительно эллиптической части и $\partial{ \bar {\beta }}{\text{/}}\partial t$ умножением на I^–1(t)

Теперь для каждого из уравнений (4.3) может быть записана формула Фейнмана–Каца:

или в векторном виде

Наконец, обратная замена $\bar {\beta }\left( {t,y} \right) = {{I}^{{ - 1}}}\left( t \right)\beta \left( {t,y} \right)$ дает

Приводя выражения, нужные для расчета по формулам (4.1) и (4.4), в дифференциальный вид, окончательно получаем следующую систему:

Здесь для удобства интегрирования введены дополнительные переменные ${{y}^{\beta }}(\tau ) \in {{\mathbb{R}}^{{{{n}_{z}}}}}$, ${{y}^{\gamma }}(\tau )\, \in \,{{\mathbb{R}}^{1}}$ и учтено, что

Полученная дифференциальная система (4.5) дает альтернативный вариант приближенного расчета оптимального управления (3.7) в отношении расчета коэффициентов $\beta \left( {t,y} \right)$ и определения его качества $J((U{\text{*}})_{0}^{T}) = V\left( {0,Y,Z} \right) = \mathbb{E}\{ {{Z}^{{\text{T}}}}\alpha \left( 0 \right)Z + {{Z}^{{\text{T}}}}\beta \left( {0,Y} \right) + \gamma \left( {0,Y} \right)\} $ в отношении расчета коэффициента $\gamma \left( {0,y} \right).$ Именно многократно моделируя решение системы (4.5) для различных начальных условий $\left( {t,y} \right),$ указанные коэффициенты можно оценивать по методу Монте-Карло, заменяя операторы $\mathbb{E}\left\{ \cdot \right\}$ в (4.5) статистическим средним. Детально соответствующий алгоритм приведен в [21]. Заметим, что альтернативным предлагаемый метод квалифицируется по отношению к традиционному способу решения параболических уравнений (3.11), (3.12) методом конечных разностей.

5. Частный случай линейного сноса. Прежде чем перейти к вопросам практической реализации и обсуждению модельного примера, рассмотрим один частный случай, когда оптимальное управление (3.7) может быть найдено аналитически, т.е. без численного решения параболического уравнения (3.11). Очевидно, что частным случаем рассматриваемой задачи является классическая линейно-гауссовская постановка, о чем уже упоминалось выше. В таком случае функция Беллмана (3.5) является квадратичной, что означает линейность коэффициентов $\beta (t,y)$ и квадратичность коэффициента $\gamma \left( {t,y} \right)$ по y, т.е. решения (3.11), (3.12) могут быть найдены в виде $\beta (t,y) = {{\beta }^{a}}(t)y + {{\beta }^{b}}(t)$ и $\gamma (t,y) = {{y}^{{\text{T}}}}{{\gamma }^{a}}(t)y + {{y}^{{\text{T}}}}{{\gamma }^{b}}(t) + {{\gamma }^{c}}(t)$ соответственно, где матрица ${{\beta }^{a}}\left( t \right) \in {{\mathbb{R}}^{{{{n}_{z}} \times {{n}_{y}}}}}$, вектор ${{\beta }^{b}}\left( t \right) \in {{\mathbb{R}}^{{{{n}_{z}}}}}$, матрица ${{\gamma }^{a}}\left( t \right) \in {{\mathbb{R}}^{{{{n}_{y}} \times {{n}_{y}}}}}$, вектор ${{\gamma }^{b}}\left( t \right) \in {{\mathbb{R}}^{{{{n}_{y}}}}}$ и функция ${{\gamma }^{c}}\left( t \right) \in {{\mathbb{R}}^{1}}$ не зависят от y. Такой случай с точки зрения новизны интереса не представляет.

Однако рассматриваемая постановка включает другой более интересный частный вариант. Из-за независимости состояния y(t), заданного уравнением (2.1), от управления u(t), которое входит линейно только в уравнение (2.2) выхода z(t), для линейности коэффициентов $\beta \left( {t,y} \right)$ достаточно линейного сноса $A\left( {t,y} \right)$ в (2.1), а линейность диффузии $\Sigma \left( {t,y} \right)$ не требуется. В отличие от случая LQG функция Беллмана при этом не будет квадратичной, решение уравнения (3.12) в виде $\gamma \left( {t,y} \right) = {{y}^{{\text{T}}}}{{\gamma }^{a}}\left( t \right)y + {{y}^{{\text{T}}}}{{\gamma }^{b}}\left( t \right) + {{\gamma }^{c}}\left( t \right)$ не представляется именно из-за нелинейности $\Sigma \left( {t,y} \right),$ но простое линейное представление оптимального управления (3.7) в этом практически важном частном случае получается.

Итак, в предположении, что $A\left( {t,y} \right) = {{A}^{a}}\left( t \right)y + {{A}^{b}}\left( t \right)$, где матрица ${{A}^{a}}\left( t \right) \in {{\mathbb{R}}^{{{{n}_{y}} \times {{n}_{y}}}}}$ и вектор ${{A}^{b}}\left( t \right) \in {{\mathbb{R}}^{{{{n}_{y}}}}}$ не зависят от y, найдем решение (3.11) в виде $\beta \left( {t,y} \right) = {{\beta }^{a}}\left( t \right)y + {{\beta }^{b}}\left( t \right)$. При этом предположении

Кроме того, учтем, что для слагаемых с $\partial {{\beta }^{{\left( i \right)}}}{\text{/}}\partial y$ в (3.11) имеет место

так что а уравнения (3.11), записанные в векторном виде, дают откуда получаются уравнения для ${{\beta }^{a}}\left( t \right)$ и ${{\beta }^{b}}(t)$: т.е. система обыкновенных линейных дифференциальных уравнений. Оптимальное управление (3.7) принимает вид

6. Вопросы численной реализации и модельный пример. В этом разделе обсуждаются некоторые аспекты практической реализации полученного оптимального управления. С этой целью рассматривается прикладной пример из области телекоммуникаций – популярного современного источника приложений для различных теоретических задач. Далее использована простая модель для показателя RTT сетевого протокола TCP [22]. Рассматривается модельный пример, поэтому параметры подобраны для лучшей выразительности иллюстративного материала. Переменная состояния RTT описывается следующей диффузией:

Здесь $\mathbb{N}\left( {M,D} \right)$ – нормальное распределение со средним M и дисперсией D.

Надо отметить, что такое уравнение (с точностью до варьируемых параметров) изначально известно как модель Кокса–Ингерсона–Росса (Cox–Ingersoll–Ross model) [23] и описывает эволюцию процентных ставок. Конечно, для описания реальных данных RTT модель типа (6.1) представляет не более чем начальное приближение, но оно полезно уже потому, что дает основу для дальнейших обобщений, уточнений модели (см., например, [7]). Свойства процесса y(t) из (6.1) хорошо изучены, в частности имеется свойство неотрицательности выборочных значений $y\left( t \right) > 0$, наличие эргодичности, известно предельное распределение и переходная вероятность. Отметим, что начальные условия в (6.1) выбраны так, что M, D значительно отличаются от моментных характеристик предельного распределения, т.е. управление ведется в рамках переходного процесса, а не в стационарном режиме.

Выход для состояния (6.1) задается уравнением

целевой функционал

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

Моделировалось $N = 1000$ траекторий $y(t),z{\text{*}}(t),~u{\text{*}}(t)$ для T = 5. Здесь через z*(t) обозначен выход (6.2), рассчитанный для оптимального управления с обратной связью $u{\text{*}} = u{\text{*}}(t,y,z)$. Кроме того, моделировались траектории ${{z}^{{(prog)}}}(t),{{u}^{{(prog)}}}(t)$ для программного управления ${{u}^{{(prog)}}}(t) = \mathbb{E}\{ u{\text{*}}(t)\} $ и ${{z}^{{(0)}}}(t),{{u}^{{(0)}}}(t)$ для неуправляемого выхода (6.2), т.е. для $u\left( t \right) = 0$. В каждом случае путем осреднения по пучку траекторий оценивались величины $J(U_{0}^{t})$. Соответственно в результате можно увидеть не только конечные значения целевой функции для разных вариантов управления $J((U{\text{*}})_{0}^{T})$, $J(({{U}^{{\left( {prog} \right)}}})_{0}^{T})$, $~J(({{U}^{{\left( 0 \right)}}})_{0}^{T})$, но и их формирование в динамике.

Описание численных процедур начнем с моделирования. Траектории y(t), z(t) моделировались приближенно, для чего уравнения (6.1), (6.2) дискретизировались по методу Эйлера с шагом ${{\delta }_{t}} = 0.001$ (тем же, что далее использовался при численной аппроксимации $\beta \left( {t,y} \right)$). Соответствующие смоделированные значения $y\left( {t - {{\delta }_{t}}} \right),z{\text{*}}\left( {t - {{\delta }_{t}}} \right)$ применялись в (3.7) для вычисления реализаций управления $u{\text{*}}\left( t \right) = u{\text{*}}(t,y\left( t \right),z{\text{*}}\left( t \right))$.

Далее заметим, что уравнение Риккати (3.8) для α(t) с учетом постоянных коэффициентов в выбранной модели выхода (6.2) и целевом функционале (6.3) решается точно. Именно

Наконец, оптимальное управление u* рассчитывалось тремя способами. Во-первых, поскольку снос в модели (6.1) линейный, то в рассматриваемом примере можно вычислить решение, не интегрируя общее уравнение (3.11) для $\beta \left( {t,y} \right)$, а решая уравнения (5.1) и вычисляя далее управление, согласно (5.2). Система (5.1) линейных уравнений для ${{\beta }^{a}}\left( t \right)$ и ${{\beta }^{b}}\left( t \right)$ решалась неявным методом Эйлера с тем же шагом ${{\delta }_{t}} = 0.001$. Соответствующее приближенное оптимальное управление обозначается далее ${{u}^{{\left( {ref} \right)}}}$, решение уравнение (3.11) – ${{\beta }^{{\left( {ref} \right)}}}\left( {t,y} \right)$. Эти аппроксимации в данном примере носят референсный (эталонный) характер.

Второй способ аппроксимации $\beta \left( {t,y} \right)$ – использование традиционного метода конечных разностей для решения параболического уравнения (3.11). Реализованные варианты и особенности применения этого метода ниже обсуждаются подробнее, а соответствующее приближенное решение обозначается ${{u}^{{\left( {net} \right)}}}$, решение уравнения (3.11) – ${{\beta }^{{\left( {net} \right)}}}\left( {t,y} \right).$

Третий способ – моделирование решений системы (4.5) и замена математического ожидания статистической оценкой (метод Монте-Карло). Соответствующие приближенные решения – ${{u}^{{\left( {MK} \right)}}}$ и ${{\beta }^{{\left( {MK} \right)}}}\left( t \right).$

Оказалось, что на графических иллюстрациях далее визуализировать разницу между реализациями управлений ${{u}^{{\left( {ref} \right)}}}\left( t \right),~{{u}^{{\left( {net} \right)}}}\left( t \right),~{{u}^{{\left( {MK} \right)}}}\left( t \right)$ не удается, поэтому в этих иллюстрациях использовано общее обозначение u*(t).

Вопросы применения метода конечных разностей для уравнения параболического типа давно и хорошо изучены [24]. Первое, что надо сделать для расчета ${{\beta }^{{\left( {net} \right)}}}\left( {t,y} \right),$ – это ограничить область значения аргумента y для формирования численной схемы расчета. С этой целью использовались N смоделированных траекторий (6.1): выборочные значения использовались для оценивания функций математического ожидания и дисперсионной, с помощью которой определялась 3σ-трубка. Рисунок 1 иллюстрирует характерные траектории (6.1), среднее значение процесса и вычисленные границы (закрашенная серым область). В расчетах граница значений $y$ задавалась добавлением еще 10% к границе 3σ-трубки, чтобы уменьшить влияние выбора граничных условий, а также учитывалась положительность y(t), поэтому далее область интегрирования $y \in [\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{y} ,\bar {y}]$ = [0, 40].

Определение границы, нужное для инициирования любой схемы метода конечных разностей, выявляет проблему отсутствия граничных условий для решаемого уравнения (3.11). Этих условий нет, и нет физических оснований для их обоснованного выбора. Из-за этого граница $\bar {y}$ была сдвинута дополнительно, а сами граничные условия приходилось выбирались волюнтаристски в надежде на нечувствительность решения в отношении этих граничных условий.

Выбор граничного условия и анализ его влияния на решение выполнялся следующим образом. Были взяты два традиционных варианта граничных условий – в задаче Дирихле: ${{\left. {\beta (t,y)} \right|}_{{y = \bar {y}}}}$ = 0, т.е. условие поглощения, и в задаче Неймана: ${{\left. {\partial \beta \left( {t,y} \right){\text{/}}\partial y} \right|}_{{y = \bar {y}}}} = 0$, т.е. условие отражения. Причем наличие в данном примере эталонного решения ${{\beta }^{{\left( {ref} \right)}}}\left( {t,y} \right)$ позволяет видеть, что оба варианта далеки от точного. Далее уравнение (3.11) решалось методом конечных разностей с использованием явной и неявной численных схем сначала для одного граничного условия, затем для второго и полученные решения сравнивались. Устойчивость в таком случае означает совпадение (близость) решений в области, выделенной по результатам моделирования, т.е. для $y \in \left[ {0,30} \right].$

Для  численного  интегрирования  в  выбранной области $\left[ {0,T} \right] \cup [\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{y} ,\bar {y}]$ формировалась сетка $\left\{ {{{t}_{k}},{{y}_{i}}} \right\}$ с шагами δ_t  (для  интегрирования по переменной t)  и  δ_y (для интегрирования по переменной y). С использованием обозначения $\beta _{i}^{k}$ для значения $\beta \left( {t,y} \right)$ в узле $\left( {{{t}_{k}},{{y}_{i}}} \right)$ производные аппроксимировались следующим образом:

где $~K = k$ для явной схемы и $K = k + 1$ – для неявной.

Шаги интегрирования δ_t, δ_y при анализе влияния граничных условий и устойчивости схем аппроксимации выбирались разными, например явная схема обсчитывалась для ${{\delta }_{t}} = 0.0002$ и ${{\delta }_{y}} = 0.1$, в том числе с учетом того, что для нее самой вопрос устойчивости ограничивает вариации шагов. Главное, что в итоге подтвердилось отсутствие влияния выбора граничных условий на результаты численного интегрирования уравнения (3.11) для $\beta \left( {t,y} \right),$ по крайней мере, для выбранной модели, а также для комбинации шагов δ_t, δ_y, обеспечивающей устойчивость схем аппроксимации. Итоговый расчет выполнен для неявной схемы с шагами ${{\delta }_{t}} = 0.001$ и ${{\delta }_{y}} = 0.01,$ и именно это решение обозначено ${{\beta }^{{\left( {net} \right)}}}\left( {t,y} \right)$, т.е. сетка содержала ${{N}_{T}} = 5000$ точек разбиения отрезка [0, T] и ${{N}_{y}} = 4000$ точек разбиения отрезка $[\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{y} ,\bar {y}]$.

Наконец, последний третий метод расчета на основании стохастического представления коэффициентов функции Беллмана дал приближенное решение ${{\beta }^{{\left( {MK} \right)}}}\left( {t,y} \right)$. Использовалась та же сетка $\left\{ {{{t}_{k}},{{y}_{i}}} \right\}$ с ${{\delta }_{t}} = 0.001$ и ${{\delta }_{y}} = 0.01$, что и в итоговом расчете методом конечных разностей, и выборка из $N \cdot {{N}_{T}} \cdot {{N}_{y}} = 2 \times {{10}^{{10}}}$ фрагментов траекторий y(t), аппроксимирующих (4.5) согласно явному методу Эйлера.

Результаты расчетов иллюстрируют рис. 2–4, на которых приведены примеры траекторий для выхода и управлений $z{\text{*}}\left( t \right),~u{\text{*}}\left( t \right)$, ${{z}^{{\left( {prog} \right)}}}\left( t \right),~{{u}^{{\left( {prog} \right)}}}\left( t \right)$ и ${{z}^{{\left( 0 \right)}}}\left( t \right),~{{u}^{{\left( 0 \right)}}}\left( t \right)$.

Рисунок 4 демонстрирует ожидаемый проигрыш неуправляемой системы ${{z}^{{\left( 0 \right)}}}\left( t \right),{{u}^{{\left( 0 \right)}}}\left( t \right)$ и хорошее качество программного управления ${{u}^{{\left( {prog} \right)}}}\left( t \right)$. Последнее обстоятельство, конечно, объясняется исключительно выбранными параметрами рассматриваемого модельного примера.

Остается сравнить результаты трех методов аппроксимации оптимального управления u*. Лучшей иллюстрацией является сравнение аппроксимаций $\beta \left( {t,y} \right)$, т.е. поверхностей ${{\beta }^{{\left( {ref} \right)}}}\left( {t,y} \right),$ ${{\beta }^{{\left( {net} \right)}}}\left( {t,y} \right),$ ${{\beta }^{{\left( {MK} \right)}}}\left( {t,y} \right).$ Выполненные расчеты иллюстрируются на рис. 5–7, представляющих эти поверхности в точках $t = $ 0, 1, 2, 3, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 4.95, 5.0 (точки выбраны там, где изменения поверхности более выражены и иллюстрации более выразительны) и целочисленных точках отрезка $y \in \left[ {0,30} \right].$

Как уже отмечалось, визуально оценить данные на этих рисунках затруднительно, но качественно можно отметить, что, с одной стороны, поверхности похожи, с другой – поверхность ${{\beta }^{{\left( {MK} \right)}}}\left( {t,y} \right)$ отличается некоторым “дрожанием” в сравнении с гладкими ${{\beta }^{{\left( {ref} \right)}}}\left( {t,y} \right)$ и ${{\beta }^{{\left( {net} \right)}}}\left( {t,y} \right),$ что является очевидным следствием аппроксимации значений $\beta \left( {t,y} \right)$ статистическими оценками по моделируемым данным.

Несколько более информативны графические материалы, приведенные на рис. 8–10. Это сечения поверхностей ${{\beta }^{{\left( {ref} \right)}}}\left( {t,y} \right),$ ${{\beta }^{{\left( {net} \right)}}}\left( {t,y} \right)$ и ${{\beta }^{{\left( {MK} \right)}}}\left( {t,y} \right)$ в тех же точках t. Причем в точном решении, представленном ${{\beta }^{{\left( {ref} \right)}}}\left( {t,y} \right),$ это прямые линии. Также прямой $\beta \left( {T,y} \right) = - 2y$ отвечает терминальное условие.

Числовые значения, иллюстрирующие отклонения $\left| {{{\beta }^{{(net)}}}(t,y) - {{\beta }^{{(ref)}}}(t,y)} \right|$ и $\left| {{{\beta }^{{(MK)}}}(t,y) - {{\beta }^{{(ref)}}}(t,y)} \right|$ в некоторых точках t, y, представлены в табл. 1 и 2. В табл. 3 приведены значения ${{\beta }^{{\left( {ref} \right)}}}\left( {t,y} \right)$, чтобы можно было оценить относительные отклонения. По всей совокупности полученных данных можно сказать, что отклонения ${{\beta }^{{\left( {MK} \right)}}}\left( {t,y} \right)$ от эталонной поверхности ${{\beta }^{{\left( {ref} \right)}}}\left( {t,y} \right)$ несколько (2–3%) меньше в большинстве точек, чем отклонения ${{\beta }^{{\left( {net} \right)}}}\left( {t,y} \right)$. Но поверхность, рассчитанная методом конечных разностей, и ее отклонения от эталона носят более регулярный характер. Точнее, в некоторых узлах относительные отклонения $\left| {{{\beta }^{{\left( {MK} \right)}}}\left( {t,y} \right) - {{\beta }^{{\left( {ref} \right)}}}\left( {t,y} \right)} \right|$ составляют до 10% и в целом, как и самой поверхности ${{\beta }^{{\left( {MK} \right)}}}\left( {t,y} \right)$, им свойственно некоторое “дрожание”, характерное для оценок Монте-Карло. Объективную оценку проводимого сравнения, конечно, дает значение целевого функционала (6.3), доставляемое управлениями ${{u}^{{\left( {ref} \right)}}},$ ${{u}^{{\left( {net} \right)}}}$ и ${{u}^{{\left( {MK} \right)}}}$. Эти данные приведены в табл. 4 и подтвердили не только высокую точность расчетов, но и превосходство метода аппроксимации, основанного на стохастическом представлении коэффициентов функции Беллмана. И хотя последний вывод сделан в отношении цели решаемой задачи управления, т.е. качества аппроксимации управления u*, данные в табл. 4 можно интерпретировать и как косвенное свидетельство о более высокой точности аппроксимации $\beta \left( {t,y} \right).$

Заключение. Представленные в работе результаты в полном объеме описывают решение задачи управления линейным выходом векторного диффузионного процесса Ито по квадратичному критерию качества. Оптимальное решение – управление с полной обратной связью – задано замкнутой системой уравнений, приближенное решение которых может быть получено известными численными методами решения обыкновенных дифференциальных уравнений и уравнений в частных производных параболического типа. Свойства оптимального управления описаны стохастическим представлением коэффициентов найденной функции Беллмана на основе формулы Фейнмана–Каца. Это стохастическое представление позволило реализовать альтернативный метод приближенного решения задачи, базирующийся на имитационном компьютерном моделировании оптимизируемой системы. Численный пример показал работоспособность обоих методов аппроксимации и небольшое преимущество метода, основанного на стохастическом представлении коэффициентов функции Беллмана.

Для дальнейшего изучения важны следующие два направления. Во-первых, интересна задача с неполной информацией, т.е. постановка, в которой диффузионный процесс (2.1) предполагается ненаблюдаемым непосредственно, а доступным только косвенным наблюдениям (2.2). В этом направлении решение в общем виде получить не удастся. Причиной являются трудности в решении возникающей в такой постановке вспомогательной задачи оценивания – оптимальной фильтрации состояния y(t) по наблюдениям z(t). Субоптимальная фильтрация позволяет предлагать приближенные решения, практически пригодные, но не очень содержательные для теоретического изучения [25].

Второе перспективное направление развития исследованной задачи – это обобщение модели диффузии на случай смешанных непрерывно-дискретных возмущений. Это важное обобщение сделает постановку более практически содержательной. Так в использованном здесь модельном примере обобщается и становится более соответствующей потребностям практики модель RTT.