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

0. Введение. Задача оптимального управления нелинейным динамическим объектом с сосредоточенными параметрами при случайных возмущениях решается, как правило, в предположении марковости и даже диффузионности (отсутствия скачков) его вектора состояния. Для его описания используется стохастическое дифференциальное уравнение (СДУ) [1], возмущаемое гауссовским белым шумом. Это позволяет свести данную стохастическую задачу к оптимизации другого динамического объекта с распределенными параметрами, заданного уравнением Фоккера–Планка–Колмогорова (ФПК) – детерминированным дифференциальным уравнением (ДУ) в частных производных (ЧП) второго порядка относительно плотности вероятности случайного состояния объекта. В случае же присутствия и скачков вектора состояния объекта необходимо учитывать порождающий их дробовой пуассоновский белый шум, что приводит к использованию более общего интегродифференциального уравнения Колмогорова–Феллера.

Обычно ищется оптимальный в среднем закон ограниченного управления при традиционном терминальном критерии качества, учитывающем потери на управление за весь отрезок времени. Тогда задачу синтеза управления необходимо решать заранее, до начала работы системы, а вид этого закона и сложность его получения существенно зависят как от полноты и точности измерений вектора состояния объекта (вида обратной связи), так и от учета или неучета предыстории изменения выхода объекта [1–5].

0.1. Статические (безынерционные) законы стохастического управления. Они возникают в следующих трех частных случаях.

1. Программное управление появляется при отсутствии измерений вектора состояния (обратной связи нет), когда управление есть функция одного только времени. Она определяется с помощью вероятностного обобщения принципа максимума Понтрягина. При этом необходимо найти частный максимум только по управлению безусловного математического ожидания стохастического аналога функции Гамильтона (параметрическая задача нелинейного программирования) и решить на всем промежутке времени управления двухточечную краевую задачу (ДТКЗ) для системы из двух ДУ в ЧП. Они определяют плотность вероятности случайного состояния управляемого объекта (уравнение ФПК или прямое уравнение Колмогорова с начальным условием) и соответствующую ей сопряженную функцию (аналог обратного уравнения Колмогорова с конечным условием). В частной линейно-квадратично-гауссовской (ЛКГ) задаче нахождение неограниченного программного управления тоже сводится к решению ДТКЗ, но уже для системы обыкновенных ДУ (ОДУ), как линейных, так и квадратичных типа Риккати.

2. Позиционное управление появляется при точном измерении всего вектора состояния объекта (полная обратная связь), когда управление есть функция еще и известного текущего состояния объекта. Она определяется путем решения в обратном времени, от момента окончания процесса управления (конечного горизонта) до его начала, вероятностного обобщения уравнения Беллмана относительно функции оставшихся средних потерь. Последнее тоже есть ДУ в ЧП, но осложненное еще и операцией нахождения частного минимума по управлению стохастического аналога функции Гамильтона. В частной ЛКГ-задаче и при отсутствии ограничений на управление это уравнение удалось решить аналитически. В результате оказалось, что и в этой стохастической задаче оптимальным выступает известный линейный регулятор Летова–Калмана, наилучший в детерминированной задаче линейно-квадратичного неограниченного управления. Коэффициент усиления этого регулятора определяется решением в обратном времени ОДУ Риккати квадратичного типа.

3. Частично позиционное управление появляется при точном измерении только части вектора состояния объекта (неполная обратная связь). В этом случае, отказываясь от процедуры восстановления всего вектора состояния, управление определяется как функция времени и только этой точно измеренной части вектора состояния объекта, т.е. является промежуточным вариантом случаев 1 и 2. Его нахождение сводится к частной максимизации по управлению условного математического ожидания стохастического аналога функции Гамильтона, учитывающего значение измеренной части вектора состояния, и решению ДТКЗ для системы уравнений из двух ДУ в ЧП в виде уравнения ФПК и сопряженного к нему. В частной ЛКГ-задаче поиск неограниченного частично позиционного управления также сводится к решению ДТКЗ, но для системы из линейных и квадратичных ОДУ, существенно более громоздкой, чем в случае 1.

Как видно, процесс нахождения даже этих простых теминально-оптимальных законов управления довольно сложен и непременно должен быть выполнен заранее. Реализация же всех этих алгоритмов безынерционного управнения в реальном масштабе времени осуществляется тривиально. В случае 1 она сводится к простому воспроизведению полученной зависимости управления от времени, а случаях 2 и 3 – к вычислению значения найденного закона управления по поступающему с точного измерителя полному текущему вектору состоянию объекта или его части, что тоже затруднений не вызывает.

0.2. Динамические (инерционные) законы стохастического управления. Существенно более проблемными являются задачи учета всей предыстории изменения выхода объекта, применимые как при неполной обратной связи, так и при стохастической обратной связи, когда измерение вектора состояния объекта осуществляется со случайными погрешностями [2, 5]. В этих случаях учитывают, что вся информация о полном векторе состояния объекта содержится в его апостериорной (условной) плотности вероятности, определяемой по всей предыстории измерений путем решения в прямом времени задачи Коши для стохастического интегродифференциального уравнения (ИДУ) Стратоновича–Кушнера (Дункана–Мортенсена–Закаи) в ЧП. Последнее представляет собой уравнение состояния абсолютно оптимального нелинейного фильтра Стратоновича. При этом различают следующие подходы к нахождению оптимального закона управления.

1. Точное решение задачи доставляет регулятор Стратоновича–Мортенсена, определяющий оптимальное управление как наилучший функционал от текущего значения апостериорной плотности вероятности. В принципе он может быть найден в результате решения в обратном времени весьма сложного функционального уравнения Беллмана–Мортенсена в вариационных производных второго порядка. Однако в процессе реализации такого управления необходимо получать апостериорную плотность (моделировать нелинейный фильтр) по уравнению Стратоновича–Кушнера и вычислять по ней значения найденного оптимального функционала в реальном масштабе времени. Выполнить это в реальном времени, как правило, невозможно.

2. Приближенное решение задачи осуществляют с помощью метода достаточных координат [2]. При этом уравнение Стратоновича–Кушнера заменяется системой из конечного числа СДУ для условных семи инвариантов, квазимоментов или других числовых характеристик апостериорной плотности (конечномерное приближение к фильтру Стратоновича). Тогда соответствующее приближение к оптимальному управлению определяется как некоторая функция текущих значений этих достаточных координат, которая находится в результате выполняемого заранее решения в обратном времени соответствующего уравнения Беллмана второго порядка. Однако размерность (порядок, число уравнений состояния) получаемого таким образом конечномерного динамического регулятора весьма велика даже при грубой гауссовой аппроксимации апостериорной плотности и резко возрастает с увеличением точности приближения. Поэтому реализовать в реальном времени такой субоптимальный алгоритм управления часто не удается. Кроме того, выбор типа достаточных координат осуществляется произвольно, а получаемый в результате алгоритм управления не является оптимальным в своем классе конечномерных динамических регуляторов.

3. Линейный инерционный регулятор справедлив только в частной ЛКГ-задаче, когда апостериорная плотность вероятности Гауссова, и при отсутствии ограничений на управление. Тогда по теореме разделения оптимальный регулятор представляет собой последовательное соединение назависимых друг от друга динамического линейного фильтра Калмана и статического линейного регулятора Летова–Калмана, объединенных обратной связью по управлению. При этом фильтр Калмана можно синтезировать как заранее, найдя решение соответствующего прямого уравнения Риккати на всем отрезке времени управления, так и в темпе со временем, решая это непростое уравнение по мере необходимости получения оценки вектора состояния, хотя скорость вычисления управления при этом существенно уменьшается. Но независимый регулятор Летова–Калмана все равно нужно синтезировать заранее путем решения соответствующего уравнения Риккати в обратном времени. После этого реализация регулятора сводится к решению либо только линейного стохастического дифференциального уравнения для оценки Калмана, либо еще и матричного квадратичного дифференциального уравнения Риккати для получения коэффициента усиления фильтра.

4. Инженерный подход к решению нелинейных, неквадратичных или негауссовских задач основан на волевой нелинейной модернизации регулятора Летова–Калмана. Вместо линейного регулятора используется псевдооптимальный в этом случае нелинейный детерминированный позиционный закон управления, не учитывающий случайные возмущения. Он находится заранее путем решения в обратном времени обычного уравнения Беллмана первого порядка. Требуемый же для его реализации вектор состояния объекта заменяется оценкой последнего, которая доставляется каким-либо из известных субоптимальных нелинейных фильтров. Простейшими из них являются учитывающие только два первых момента ковариационные приближения: линеаризованное вроде нелинейного обобщенного (расширенного) фильтра Калмана, более точное гауссовское приближение (фильтр нормальной аппроксимации) и т.п. Однако скорость работы и такого инженерного регулятора падает при увеличении как размерности фильтра, так и сложности его уравнений. Так размерность даже только ковариационных приближений к фильтру квадратична к размерности объекта управления. Далее степень этой зависимости (кубическая и т.п.) резко растет при увеличении точности фильтра за счет учета достаточных координат старших порядков.

0.3. Конечномерные динамические законы стохастического управления. Строгой альтернативой указанному инженерному подходу являются два следующих подхода к созданию быстрых инерционных алгоритмов управления.

1. Оперативно-оптимальный регулятор [6] определяется как последовательное соединение быстрого динамического условно-оптимального фильтра Пугачева малого порядка [7, 8] и статического преобразователя Казакова состояния фильтра (вектора оценки) в управление. Оба этих блока имеют заданную, заранее выбранную и не оптимальную нелинейную структуру, а их параметры находятся из условия мгновенной (локальной) оптимальности. Последнее позволило удачно обойтись без решения ДТКЗ, так как достаточно либо найти априорную плотность вероятности путем решения задачи Коши для уравнения ФПК относительно общего вектора состояния объекта и фильтра, либо выполнить многократное статистическое моделирование соответствующих СДУ объекта и фильтра последовательным методом Монте-Карло. Подчеркнем, что обе эти операции в принципе возможно выполнить в темпе со временем.

2. Регулятор терминально-оптимальной структуры ищется как более общая зависимость управления от вектора состояния регулятора желаемой размерности [9]. При этом нелинейные структурные функции ДУ состояния регулятора, его начальное состояние и сама функция его выхода находились минимизацией среднего значения традиционного терминального критерия оптимальности. Поэтому их получение требовало заранее решить ДТКЗ для системы из двух нелинейных ИДУ в ЧП относительно совместной плотности вероятности состояний объекта и регулятора, а также сопряженной ей функции тех же переменных. Применение и здесь метода достаточных координат снова приводило к ДТКЗ для системы ОДУ относительно бесконечного (но урезаемого) количества числовых характеристик этой плотности и соответствующих им сопряженных переменных, что также сложно для практического применения.

Итак, параметры быстрого оперативно-оптимального регулятора Пугачева–Казакова могут быть синтезированы в темпе со временем и не используют информацию о будущем виде как объекта, так и цели управления, позволяя изменять их со временем, но его неоптимальную структуру нужно подбирать методом проб и ошибок. Свободный же от этого недостатка регулятор терминально-оптимальной структуры требует такую информацию о будущем и его приходится синтезировать заранее. Отметим, что детерминированная версия синтеза весьма перспективного оперативно-оптимального (локально-оптимального) позиционного закона управления, альтернативная случаю 2 из разд. 0.1, успешно применена для наведения летательных аппаратов [10].

В настоящей работе с целью радикального упрощения процедуры предварительного синтеза терминально-оптимальной структуры быстрого конечномерного динамического регулятора предлагается заменить традиционный критерий его оптимальности на переменный. Последний учитывает только текущие потери и не позволяет управлять будущим поведением объекта. Показано, что нахождение структуры такого регулятора сводится к оперативному решению начальной задачи Коши для ИДУ в ЧП относительно условной плотности вероятности, характеризующей состояние объекта при знании текущих величин его выхода и состояния регулятора. Это можно осуществить численно последовательным методом Монте-Карло или приближенно методом достаточных координат. Пример построения такого регулятора в частном случае линейного объекта и квадратично-биквадратного критерия будет рассмотрен во второй части статьи.

Сначала приведем вариант точных измерений части переменных состояния управляемого объекта (неполные измерения, парная марковская модель системы управления). Альтернативный случай неточных измерений вектора состояния объекта управления (стохастические измерения, скрытая марковская модель системы управления) по отношению к нему формально является частным и будет представлен в разд. 10.

1. Постановка задачи при неполных точных измерениях. 1.1. Уравнения объекта управления. Пусть функционирующий на отрезке времени $t \in \left[ {0,T} \right]$ динамический объект с частично измеряемым вектором состояния $({{X}_{t}},{{Y}_{t}})$, где ${{X}_{t}} \equiv X(t) \in {{\mathbb{R}}^{n}}$ – n-мерный неизмеряемый вектор, ${{Y}_{t}} \in {{\mathbb{R}}^{m}}$ – m-мерный измеряемый вектор (вектор выхода), управляемый l-мерным вектором ${{U}_{t}} \in \Omega $ из в общем случае ограниченной области $\Omega \subset {{\mathbb{R}}^{l}}$ и возмущаемый k-мерным стандартным гауссовским белым шумом ${{V}_{t}} \in {{\mathbb{R}}^{k}}$, описывается системой СДУ:

Здесь шум определяется как производная по времени ${{V}_{t}} = \dot {W}$ центрированного и нормированного винеровского процесса ${{W}_{t}}$, который не зависит от случайных начальных условий ${{X}_{0}},\;{{Y}_{0}}$. Их совместная плотность вероятности ${{p}_{0}}(x,y)$ известна частично: с точностью до условной плотности вероятности ${{\rho }_{0}}(x\,{\text{|}}\,y)$ неизмеряемой компоненты X₀, тогда как маргинальная (частная) плотность вероятности ${{q}_{0}}(y)$ начального измерения Y₀ может быть произвольной. Пусть также эти начальные условия имеют конечные вторые моменты где ${\text{M}}$ – оператор математического ожидания.

Отметим, что каждое из уравнений (1.1) может быть записано в эквивалентной дифференциальной или интегральной форме, например

Здесь первый из интегралов от случайной функции времени понимается как среднеквадратический, а последний – как стохастический интеграл Ито.

Предположим, что определяющие нелинейную структуру объекта (1.1) вектор-функции $a(t,x,y,u) \in {{\mathbb{R}}^{n}}$, $c(t,x,y,u) \in {{\mathbb{R}}^{m}}$ и матричные функции $B(t,x,y,u) \in {{\mathbb{R}}^{{n \times q}}}$, $D(t,x,y,u) \in {{\mathbb{R}}^{{m \times q}}}$ удовлетворяют известным достаточным условиям существования и единственности сильного (с точностью до вероятностной эквивалентности реализаций) непрерывного по времени решения СДУ Ито [11, с. 470]. Это условия липшицевой непрерывности данных функций по переменным x, y и ограниченной скорости их роста при $\left| x \right| + \left| y \right| \to \infty $. В результате гарантируется конечность вторых моментов случайных векторов ${{X}_{t}},\;{{Y}_{t}}$ и в любой другой момент времени

Если, кроме того, эти функции еще и непрерывны по времени t, то из (1.1) следует [11, с. 478], что случайный вектор ${{(X_{t}^{{\text{T}}}\,\,Y_{t}^{{\text{T}}})}^{{\text{T}}}}$ будет диффузионным (непрерывным, без скачков) марковским процессом с вектор-функцией сноса ${{({{a}^{{\text{T}}}}\,\,{{c}^{{\text{T}}}})}^{{\text{T}}}}$ и матрицей диффузии

Следовательно, его совместная плотность вероятности $p(t,x,y)$ существует и является обобщенным решением (в смысле удовлетворения интегральному тождеству, см. разд. 2) соответствующего уравнения ФПК с начальным условием $p(0,x,y) = {{p}_{0}}(x,y)$. Здесь n × n-матрица Q, m × m-матрица R и n × m-матрица S имеют смысл условных интенсивностей (собственных и взаимной соответственно) двух в общем случае зависимых белых шумов системы (1.1), а именно порождающего $B{{V}_{t}}$ и измерительного $D{{V}_{t}}$. Будем считать, что матрица интенсивности шума измерения ${{Y}_{t}}$ положительно определена $R > 0$, так что при любых $t,\;x,\;y,\;u$ существует обратная матрица R^–1.

Отметим, что вариант возмущения двух уравнений (1.1) для векторов ${{X}_{t}}$, ${{Y}_{t}}$ разными независимыми белыми шумами $V_{t}^{{(1)}}$, $V_{t}^{{(2)}}$, когда

является частным случаем системы (1.1), если

1.2. Предлагаемый регулятор. Наиболее эффективно терминально-оптимальное управление состоянием стохастического объекта (1.1) по его выходу ${{Y}_{t}}$ выполняют различные приближения к динамическому регулятору Стратоновича–Мортенсена. Они учитывают информацию обо всех предыдущих измерениях, обеспечивая функциональную зависимость команды управления ${{U}_{t}}$ от них [2, 5]:

Однако их, во-первых, необходимо синтезировать заранее, используя информацию о предполагаемых в будущем на всем отрезке времени управления свойствах объекта в виде его математической модели. Во-вторых, они, как правило, плохо реализуемы в реальном масштабе времени из-за большой размерности их вектора состояния (см. разд. 0.2). Поэтому наложим на эту функциональную зависимость ограничение конечномерности и будем искать ее из условия оперативной оптимальности, предполагающего знание только текущих и прошлых свойств объекта, но не будущих.

Функциональное преобразование всех измерений $Y_{0}^{t}$ в управление ${{U}_{t}}$ определим [9] в виде дифференциального уравнения для запоминающего информацию о прошлых измерениях p‑мерного (p = 1, 2, …) вектора состояния ${{Z}_{t}} \in {{\mathbb{R}}^{p}}$ регулятора

которое обеспечивает совместную марковость ${{Z}_{t}}$ вместе с $({{X}_{t}},{{Y}_{t}})$, и формулы выхода регулятора

Здесь неизвестными являются четыре структурные функции регулятора: начального значения $h(y) \in {{\mathbb{R}}^{p}}$, смещения$f(t,y,z) \in {{\mathbb{R}}^{p}}$, усиления $G(t,y,z) \in {{\mathbb{R}}^{{p \times m}}}$ и выхода $u(t,y,z) \in {{\mathbb{R}}^{l}}$. При этом размерность p вектора состояния ${{Z}_{t}}$ (порядок регулятора) можно выбирать, например, из условия компромисса между эффективностью регулятора, измеряемой получаемым значением критерия оптимальности, и трудоемкостью его практической реализации в темпе со временем.

Отметим, что, умножив уравнение (1.3) на dt, получим дифференциальную форму записи уравнения состояния искомого регулятора:

Замечание 1. Случай отказа от накопления информации об измерениях $Y_{0}^{t}$ можно рассматривать как частный вариант предлагаемого регулятора порядка p = 0. Действительно, тогда уравнение (1.3) исчезает, а управление (1.4) принимает вид частично позиционного:

А если еще и весь вектор состояния объекта измеряется точно, т.е. ${{Y}_{t}} = {{X}_{t}}$, то имеем позиционную версию этого регулятора:

1.3. Оперативная оптимальность. В отличие от [9], где использовался традиционный постоянный во времени терминальный критерий, теперь неизвестные структурные функции регулятора $f( \cdot )$, $G( \cdot )$, $h( \cdot )$, $u( \cdot )$ найдем из другого условия оптимальности. В качестве показателя эффективности управления выберем лишь накопленные к моменту времени $t$ потери от процессов, происходящих как в объекте (1.1), так и в регуляторе (1.3), (1.4). Эти потери на каждой реализации определим в виде нестационарной версии случайного функционала Больца

и будем искать структуру регулятора из требования минимума его среднего значения:

Здесь интегрант $\varphi (t,x,y,z,u) \geqslant 0$ – непрерывная функция своих аргументов, определяющая прошлые (интегрируемые) потери, тогда как терминант $\psi (t,x,y,z) \geqslant 0$ – дифференцируемая функция $\psi (t,x,y,z) \in {{\mathbb{C}}^{{1,2}}}([0,T],{{\mathbb{R}}^{{n + m + k}}})$, которая характеризует только текущие (терминальные) потери. Будем считать обе эти функции выпуклыми по всем их аргументам, кроме времени t. При этом интегрант $\varphi ( \cdot )$ непременно должен зависеть хотя бы только от переменной управления

что исключает из рассмотрения аналог задачи оптимального быстродействия, когда φ = 1. Также обязательной является и зависимость терминанта $\psi ( \cdot )$ от переменной состояния регулятора z, которая учитывает расходы на его работу. Этим достигается энергетическая экономичность предлагаемого регулятора. В частности, первая производная терминанта по z не может быть нулевой, а матрица его вторых производных по z должна быть положительно определена:

Отметим, что в (1.6) интеграл по времени от случайной функции $\Phi (\tau ) = \varphi (\tau ,{{X}_{\tau }},{{Y}_{\tau }},{{Z}_{\tau }},{{U}_{\tau }})$ тоже понимается как среднеквадратический, а операции усреднения и интегрирования по времени перестановочны:

как следствие из теоремы Фубини о сведении двойного интеграла к повторному.

Из-за известной невозможности обеспечить минимум нестационарного критерия (1.6) в каждый момент времени t [12, 13] рассмотрим основанный на идее припасовывания следующий метод оптимизации структуры предлагаемого регулятора. Он требует минимизации начального значения критерия выбором функции $h(y)$ начального состояния регулятора:

и минимизации приращения критерия I_t за сколь угодно малый промежуток времени $[t,t + \Delta t]$:

Последнее приводит [14–16] к минимизации скорости изменения этого функционала в каждый момент времени:

Это условие называют локальным критерием качества управления. Отметим, что из неотрицательности функционала (1.6) следует ограниченность снизу этого критерия в каждый момент времени:

Замечание 2. Другим частным случаем рассматриваемой задачи является синтез среднеквадратической оптимальной структуры конечномерного фильтра малого порядка [13]. Действительно, пусть исходный объект (1.1) не управляется:

и требуется по накапливаемым измерениям его выхода ${{Y}_{t}}$ оперативно получать оценку $\hat {X}_{t}^{'}$ наиболее интересной $n{\kern 1pt} '$-мерной части $X_{t}^{'}$ неизмеряемого n-мерного вектора ${{X}_{t}}$:

Тогда в качестве такой оценки можно взять решение дифференциального уравнения (1.3)

и выбрать критерий его оптимальности в виде среднего квадрата ошибки оценивания где ${{L}_{t}} > 0$ – положительно определенная $n{\kern 1pt} '\; \times n{\kern 1pt} '$-матрица весовых коэффициентов. Этот критерий является частным случаем функционала (1.6) и требуется найти только три структурные функции $f( \cdot )$, $G( \cdot )$, $h( \cdot )$ уравнения (1.3), причем его произвольный порядок фиксируется как $p = n{\kern 1pt} '$. Дополнительным для фильтра является требование несмещенности оценки M[Z_t] = = ${\rm M}[X_{t}^{'}]$.

2. Сведение стохастической задачи к детерминированной. Подставим формулу выхода регулятора (1.4) в уравнения объекта (1.1) и объединим полученную систему с уравнением состояния регулятора (1.3). Тогда общий вектор состояний объекта и регулятора ${{{\rm K}}_{t}} = {{(X_{t}^{{\text{T}}}\,Y_{t}^{{\text{T}}}\,Z_{t}^{{\text{T}}})}^{{\text{T}}}}$ также удовлетворяет СДУ Ито:

где вектор-функция сноса $\omega ( \cdot )$ и матричная функция диффузии $\Sigma ( \cdot )$ имеют вид

Здесь и далее верхними индексами для краткости отмечены зависимости соответствующих сложных функций от структурных функций регулятора $f( \cdot )$, $G( \cdot )$, $u( \cdot )$, например

Из (2.1) следует, что при условии дифференцируемости по переменным состояния x, y, z функции сноса $\omega ( \cdot )$ один раз, а функции диффузии $\Sigma ( \cdot )$ дважды мгновенная (одноточечная) плотность вероятности $r(t,x,y,z)$ процесса ${{{\rm K}}_{t}}$ также дважды дифференцируема по ним и один раз по времени t и удовлетворяет уравнению ФПК:

Здесь $K_{{xyz}}^{{ufG}}$ – прямой производящий оператор этого процесса

в котором $\nabla = {{(\nabla _{x}^{{\text{T}}},\nabla _{y}^{{\text{T}}},\nabla _{z}^{{\text{T}}})}^{{\text{T}}}}$ – оператор градиента по векторной переменной $\kappa = (x,y,z)$, тогда как tr – оператор следа матрицы. Начальным для уравнения (2.2) является условие где $\delta ( \cdot )$ – функция Дирака, а естественные граничные условия на бесконечности по переменным x, y, z для плотности $r( \cdot )$ и вектора потока вероятности $\pi = \omega r - 0.5\,{{[{{\nabla }^{{\text{T}}}}(\Sigma {{\Sigma }^{{\text{T}}}}r)]}^{{\text{T}}}}$ будут нулевыми.

Для ликвидации отмеченных жестких условий гладкости отметим, что при этих граничных условиях имеет место справедливая для любой функции $\eta \left( {t,\kappa } \right)$, дважды непрерывно дифференцируемой по переменной $\kappa \in {{\mathbb{R}}^{{n + m + p}}}$

формула интегрирования по частям [2]

Здесь $K_{{xyz}}^{{*ufG}}$ – обратный производящий оператор процесса K_t, сопряженный оператору $K_{{xyz}}^{{ufG}}$:

который не требует дифференцируемости функций сноса $\omega ( \cdot )$ и диффузии $\Sigma ( \cdot )$ процесса K_t.

Запишем этот оператор подробнее, учитывая вид этих функций. Используя введенные выше обозначения $B{{B}^{{\text{T}}}} = Q,$ $D{{D}^{{\text{T}}}} = R,$ $B{{D}^{{\text{T}}}} = S$ и группируя следы от матриц одинаковых порядков n × n, m × m, p × p, получим

Здесь по-прежнему одинарными или двойными нижними индексами x, y, z обозначены столбцы первых и матрицы вторых частных производных скалярной функции $\eta \left( {t,x,y,z} \right)$ соответственно. Например ${{\eta }_{x}} = {{\nabla }_{x}}\eta $ есть n-мерный столбец, тогда как матрица ${{\eta }_{{xy}}} = {{\nabla }_{x}}\nabla _{y}^{{\text{T}}}\eta $ имеет порядок (размер) n × m.

Отметим, что в обратном операторе $K_{{xyz}}^{{*ufG}}$ переменные u,   f, G являются обычными параметрами, так как никакие действия в отношении соответствующих им функций $u( \cdot )$, $f( \cdot )$, $G( \cdot )$ не производятся. Причем функция $u( \cdot )$ просто входит в аргументы сложных функций сноса ${{a}^{u}}( \cdot ),\,\;{{c}^{u}}( \cdot )$ и диффузии ${{Q}^{u}}\,,\;{{S}^{u}},\;{{R}^{u}}$ системы (1.1).

Используя формулу (2.4), вместо уравнения (2.2) получим известное интегральное тождество, определяющее его обобщенное решение [17]:

где угловыми скобками для краткости обозначен интеграл усреднения функции $\eta ( \cdot )$ по всем ее аргументам с весом в виде совместной плотности вероятности $r( \cdot )$:

Здесь и далее, если не оговорено иначе, интегралы берутся по всему соответствующему евклидову пространству, например

Подчеркнем, что тождество (2.6) не накладывает на функции сноса и диффузии случайного процесса K_t, а также и на саму плотность вероятности отмеченные выше довольно жесткие условия гладкости по переменным состояния x, y, z.

Подставляя формулу выхода регулятора (1.4) также и в функционал (1.6), операцию усреднения М в нем запишем с помощью двух интегралов (2.7):

Начальное значение I₀ этого переменного по времени функционала, согласно (2.3), зависит только от функции $h(y)$ начального состояния регулятора:

Поэтому его минимизация (1.7) при оговоренном условии выпуклости терминанта $\psi ( \cdot )$ и частично заданном распределении начального состояния (1.2) дает параметрическое алгебраическое уравнение для определения наилучшей функции $h(y)$ начального состояния регулятора:

В результате остается решить задачу локально-оптимального в смысле (1.8) управления детерминированным распределенным объектом, определяемым уравнением в частных производных (2.2) и его начальным условием (2.3), путем выбора наилучших функций управляющих полей $u( \cdot )$, $f( \cdot )$, $G( \cdot )$, которые доставляют минимум производной нестационарного функционала (2.8) в любой момент времени $t \in \left[ {0,T} \right]$.

3. Явный вид производной критерия. Для решения задачи минимизации (1.8) продифференцируем нестационарный функционал (2.8) по времени t. Так как

то, опуская далее для краткости аргументы $t,\;x,\;y,\;z$ функций ${{\varphi }^{u}}( \cdot )$, $\psi ( \cdot )$, $r( \cdot )$, получим или, используя во втором слагаемом тождество (2.6),

Искомые функции входят в это выражение только в виде их сечений при фиксированном t:

Выделяя здесь слагаемые, которые зависят или не зависят от этих функций:

приходим к необходимости минимизировать далее, согласно (1.8), такой функционал

При этом, согласно (1.9), и он должен быть ограничен снизу:

Учитывая вид (2.5) сопряженного оператора $K_{{xyz}}^{{*ufG}}$, выделим в (3.1) группы слагаемых, зависящих от разных структурных функций предлагаемого регулятора:

Это такие совместные средние

В последнем выражении вспомогательная m × p-матрица H^u имеет вид

а при его получении были использованы известные свойства следа матрицы

4. Упрощение производной критерия. Так как искомые функции $u(t,y,z)$, $f(t,y,z)$, $G(t,y,z)$ от переменной $x$ не зависят, то функционал (3.3) можно упростить, заменив совместную плотность вероятности $r( \cdot )$ случайного вектора ${{{\rm K}}_{t}} = {{(X_{t}^{{\text{T}}},Y_{t}^{{\text{T}}},Z_{t}^{{\text{T}}})}^{{\text{T}}}}$ произведением маргинальной (частной) плотности $s( \cdot )$ случайных величин ${{Y}_{t}},\;{{Z}_{t}}$ на соответствующую условную плотность $\rho ( \cdot )$:

Тогда интеграл совместного среднего (2.7) можно представить в виде повторного, состоящего из внутреннего интеграла усреднения по условной плотности $\rho ( \cdot )$ и внешнего интеграла усреднения по маргинальной плотности $s( \cdot )$:

Поэтому минимизируемый функционал (3.3) со слагаемыми (3.4) принимает вид маргинального среднего:

от подынтегральной функции условного среднего: слагаемые которой, опуская тройку их очевидных первых аргументов t, y, z, имеют, согласно (3.4), такой вид:

Здесь уже зависимости функций ${{\xi }_{1}}(u),\,\,{{\xi }_{3}}(u,G)$ от переменной $u$ определяются исходными зависимостями от нее функций сноса и диффузии управляемой системы (1.1), например $a(t,x,y,u)$, $Q(t,x,y,u)$ и т.п., а также интегранта $\varphi (t,x,y,z,u)$ критерия (1.6). При этом функция ${{\xi }_{2}}(f)$ линейна по f, тогда как функция ${{\xi }_{3}}(u,G)$ квадратична по G.

Маргинальное представление минимизируемого функционала (4.2) существенно упрощает задачу его минимизации (3.1). Действительно, из неотрицательности маргинальной плотности $s(t,y,z) \geqslant 0$ и свойства монотонности операции интегрирования легко получить следующее утверждение.

Лемма 1 [13]. Для отыскания минимума функционала маргинального среднего (4.2) достаточно найти минимум его усредняемой функции $\xi ( \cdot )$ по части ее аргументов $v = (u,f,G)$ при любых значениях других ее аргументов t, y, z:

В результате задача синтеза регулятора сведена к задаче на частный минимум по переменным u, f, G функции условного среднего (4.3):

При этом, согласно условию (3.2), эта функция тоже должна быть ограничена снизу:

5. Нахождение функции смещения регулятора. Учтем условие ограниченности снизу (4.6) минимизируемой функции (4.3). Так как ее слагаемое ${{\xi }_{2}}(f)$ оказалось неограниченным из-за его линейности по своему аргументу  f, то для этого достаточно, согласно (4.4), обнулить соответствующий коэффициент пропорциональности:

В результате сама производная критерия оказалась независимой (инвариантной) от функции $f( \cdot )$ смещения регулятора

т.е. ее выбор на значение основного критерия (1.6) никак не влияет. Это позволяет задать ее, для простоты, равной нулю: либо искать из условия инвариантности (5.1) или какого-то другого дополнительного условия.

Для использования равенства (5.1) выясним влияние искомой функции смещения $f( \cdot )$ на условную плотности $\rho ( \cdot )$. Нетрудно показать аналогично случаю неуправляемой системы, что из-за произвольности в плотности вероятности (1.2) начального состояния объекта управления (1.1) закона распределения ${{q}_{0}}(y)$ начального измерения ${{Y}_{0}}$ имеют место следующие утверждения.

Лемма 2 [13]. Если условная плотность вероятности $\rho (t,x\,{\text{|}}\,y,z)$ из (4.1) дважды непрерывно дифференцируема по двум своим аргументам-условиям y, z, то она определяется начальным значением $\rho (0,x\,{\text{|}}\,y,z) = {{\rho }_{0}}(x\,{\text{|}}\,y)$ и аналогичным (2.6) интегродифференциальным тождеством

в котором оператор $K_{{xyz}}^{{*ufG}}$ задается выражением (2.5), а новый оператор $L_{{yz}}^{{*ufG}}$ имеет вид причем в коэффициенты последнего входят условные средние $\bar {c}$, $\bar {R}$. Здесь и далее чертой сверху обозначаются функции условного среднего

Следствие 1 [13]. При более жестких условиях гладкости условная плотность вероятности $\rho (t,x\,{\text{|}}\,y,z)$ удовлетворяет следующему интегродифференциальному уравнению в частных производных:

Используя эти соотношения, получим следующее утверждение.

Теорема 1. Если p × p-матрица ${{\nabla }_{z}}{{\left( {\overline {{{\psi }_{z}}} } \right)}^{{\rm T}}} - \overline {{{\psi }_{{zz}}}} $ из апостериорных средних производных терминанта $\psi (t,x,y,z)$ критерия (1.6) невырождена, то функция смещения $f(t,y,z)$ уравнения состояния (1.3) регулятора, удовлетворяющая условию независимости (5.1), может быть определена по формуле

Доказательство теоремы 1. Дифференцируя равенство (5.1) по времени и используя (5.3), получим уравнение

Выделим здесь слагаемые с функцией смещения $f( \cdot )$. Из выражений для операторов (2.5), (5.4) следует, что она входит в них линейно:

В результате имеем линейное уравнение

или, учитывая независимость искомой функции $f(t,y,z)$ от переменной интегрирования x,

Это уравнение при невырожденности p × p-матрицы его коэффициентов ${{\nabla }_{z}}{{(\overline {{{\psi }_{z}}} )}^{{\rm T}}} - \overline {{{\psi }_{{zz}}}} $ действительно имеет решение (5.7). Теорема 1 доказана.

Таким образом условие (5.1), обеспечивающее ограниченность снизу локального критерия, позволило выразить функцию смещения $f( \cdot )$ регулятора через его функции усиления $G( \cdot )$ и выхода $u( \cdot )$, а также через условную плотность вероятности $\rho ( \cdot )$.

6. Минимизация производной критерия. Теперь при выполнении условий (5.1) или (5.2) минимизируемая в (4.5) функция условного среднего (4.3) зависит лишь от двух искомых параметров:

причем ее второе слагаемое квадратично по G.

Для ее минимизации по двум этим параметрам используем известный принцип сечений [18], согласно которому данную операцию будем проводить последовательно:

Сначала для любых допустимых значений переменной u найдем условную (частную) минималь $\xi $ по G:

а затем минимизируем по переменной $u$ функцию частного минимума:

Тогда очевидным является следующее утверждение.

Теорема 2. Если функция (6.3) является выпуклой по переменной u, то оптимальная функция выхода регулятора определяется как единственное решение задачи параметрического нелинейного программирования:

а ее подстановка в условно-оптимальную функцию (6.2) позволяет найти и оптимальную функцию усиления регулятора:

Рассмотрим эти операции подробнее.

7. Нахождение функции усиления регулятора. Из (6.1) следует, что операция (6.2) сводится к минимизации по матричному параметру G квадратичной функции ${{\xi }_{3}}(u,G)$ из (4.4). Аналогично [19] получим следующий результат.

Теорема 3. Если матрицы ${{\psi }_{{zz}}}( \cdot )$, $R( \cdot )$ положительно определены, то условно-оптимальная функция усиления (6.2) конечномерного регулятора находится из линейного алгебраического уравнения

где $ * * $ – знак умножения (свертки) по двум индексам четырехмерной $p \times p \times m \times m$-матрицы $\bar {\Lambda }$ на плоскую p × m-матрицу $\tilde {G}$, которое возвращает p × m-матрицу, тогда как матрицы коэффициентов и свободных членов имеют вид

Последние зависят от переменной управления u через входящие в них функции $c( \cdot ),S( \cdot ),R( \cdot )$.

Доказательство теоремы 3. Для нахождения частного минимума функции ${{\xi }_{3}}( \cdot )$ по матричному параметру $G$ воспользуемся необходимым и достаточным условиями экстремума. Применяя известные формулы матричного дифференцирования, имеем

Учитывая симметричность матриц R и ${{\psi }_{{zz}}}$, находим необходимое условие экстремума в виде равенства нулю матрицы частных производных:

а также выражение для четырехмерной матрицы второй производной:

Здесь $ \otimes $ – символ прямого (внешнего) произведения матриц, а положительная полуопределенность второй производной следует из положительной определенности матриц ${{\psi }_{{zz}}}( \cdot ) > 0$, $R( \cdot ) > 0$ и неотрицательности условной плотности $\rho ( \cdot ) \geqslant 0$. Следовательно, за исключением точек с нулевыми значениями условной плотности, на решениях уравнения (7.3) функция ${{\xi }_{3}}( \cdot )$ действительно достигает искомого минимума.

Теперь рассмотрим уравнение (7.3) подробнее. Используя независимость искомой функции $G(t,y,z)$ от переменной интегрирования, вынесем ее из-под знака интеграла. В результате действительно получим линейное уравнение (7.1) с четырехмерной матрицей коэффициентов. Теорема 3 доказана.

8. Нахождение функции выхода регулятора. Конкретизируем функцию частного минимума (6.3).

Лемма 3. Целевая функция задачи программирования (6.4) имеет вид условного математического ожидания стохастического аналога функции Гамильтона:

Доказательство леммы 3. Действительно, вид (8.1) функции $\mu ( \cdot )$ устанавливается, согласно (6.3), простой подстановкой первого и третьего из выражений (4.4) в их сумму (6.1), а также заменой в полученном выражении параметра G на условно-оптимальную функцию усиления $\tilde {G}( \cdot )$. Лемма 3 доказана.

При наличии ограничения на управление $u \in \Omega $ найти оптималь (6.4) аналитически вряд ли возможно. Поэтому найдем градиент функции (8.1), который позволит применить какой-либо из известных численных методов минимизации первого порядка, а в случае отсутствия ограничения на управление при $\Omega = {{\mathbb{R}}^{l}}$ даст уравнение для определения точки минимума. Дифференцируя скалярную функцию (8.1) по векторной переменной u, получим такой результат.

Следствие 2. Если минимизируемая по переменной управления u функция (8.1) дифференцируема по ней, то ее частный градиент по этой переменной ${{\mu }_{u}}$ определяется как условное среднее:

следующей l-мерной вектор-функции:

Здесь нижним индексом u обозначены матрицы Якоби первых частных производных по этой переменной известных функций исходной системы (1.1), а также условно-оптимальной функции усиления регулятора $\tilde {G}( \cdot )$ и вспомогательной функции $H( \cdot )$. При этом ${{a}_{u}}$ – $l \times n$-матрица, ${{c}_{u}}$ – $l \times m$-матрица, ${{Q}_{u}}$ – $n \times n \times l$-матрица, ${{S}_{u}}$ – $n \times m \times l$-матрица, ${{R}_{u}}$ – $m \times m \times l$-матрица, ${{\tilde {G}}_{u}}$и ${{H}_{u}}$ – $p \times m \times l$-матрицы, тогда как $ * $ – знак операции умножения (свертки) трехмерной матрицы по двум ее индексам на согласованную с ней по размерности плоскую матрицу, которая возвращает l-мерный вектор-столбец.

Следствие 3. В случае отсутствия ограничения на управление, когда ${{U}_{t}} \in \Omega = {{\mathbb{R}}^{l}}$, при условии выпуклости интегранта $\varphi ( \cdot )$ и терминанта $\psi ( \cdot )$ функционала (1.6) по переменной u из (8.2) следует необходимое и достаточное условие оптимальности функции выхода регулятора:

9. Алгоритмы синтеза регулятора. Таким образом, функция $h(y)$ начального состояния оптимального в смысле (1.6) регулятора (1.3), (1.4) находится из алгебраического уравнения (2.9), тогда как остальные его структурные функции $f(t,y,z)$, $G(t,y,z)$, $u(t,y,z)$ парами формул (5.2), (5.7), а также (7.1), (7.2) и (6.4), (8.1) соответственно выражены через сечение условной плотности вероятности $\rho (t, \cdot \,{\text{|}}\,y,z)$. Поэтому при неизвестном законе распределения ${{q}_{0}}(y)$ начального измерения Y₀ определение оптимальной структуры регулятора сводится к нахождению плотности $\rho ( \cdot )$ путем решения задачи Коши для нелинейного ИДУ в ЧП (5.6). Предварительно в это уравнение необходимо подставить указанные выражения функций $f( \cdot )$, $G( \cdot )$, $u( \cdot )$ через $\rho ( \cdot )$, что усложняет процедуру его решения.

Если же распределение ${{q}_{0}}(y)$ задано, то из (1.2) известна и совместная плотность вероятности ${{p}_{0}}(x,y)$ всей начальной пары ${{X}_{0}},\;{{Y}_{0}}$. Тогда появляется возможность найти совместную плотность $r(t,x,y,z)$ интегрированием более простого ФПК-уравнения (2.2), но с указанными выше подстановками и представлением (4.1) условной плотности через совместную. Следовательно, становится известным и минимальное значение $I_{t}^{{\min }}$ критерия (1.6), которое определяется по (2.8).

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

Удобная численная альтернатива этому – возможность непосредственного определения структурных функций регулятора последовательным методом Монте-Карло. Необходимо осуществить многократное и пошаговое во времени статистическое моделирование СДУ объекта (1.1) и регулятора (1.3) с помощью известных разностных схем вроде Эйлера–Маруямы. При этом шаги интегрирования $k = 0,1,2, \ldots $ следует чередовать с нахождением по каждому полученному в момент времени ${{t}_{k}}$ ансамблю из довольно большого числа реализаций случайных величин $X({{t}_{k}}),\;Y({{t}_{k}}),\;Z({{t}_{k}})$ выборочных приближений к сечениям при $t = {{t}_{k}}$ требуемых функций условных средних (5.5), что аналогично получению гистограммы сечения частной плотности $s(t,y,z)$. Более подробно эта процедура описана в [13]. Тогда и минимальное значение $I_{t}^{{\min }}$ критерия (1.6) легко найти как среднее арифметическое набора соответствующих реализаций случайного функционала Больца (1.5).

Однако построение с достаточной точностью гистограммы распределения системы случайных величин является технически довольно громоздкой процедурой, особенно при их большом количестве. Второй недостаток этого способа состоит в получении структуры регулятора лишь в виде сеточных функций, для хранения которых потребуется большой объем памяти или их дополнительная обработка. Поэтому актуально использование известных способов аппроксимации законов распределения типа гауссовского приближения или применение тейлоровской линеаризации нелинейностей объекта управления (1.1). В частном случае применение гауссовского приближения будет продемонстрировано во второй части данной статьи.

10. Случай неточных измерений состояния объекта. Пусть теперь выход ${{Y}_{t}}$ не влияет на состояние ${{X}_{t}}$. Тогда система уравнений (1.1) распадается на независимое уравнение состояния объекта управления

и зависимое уравнение измерителя, в общем случае управляемого и возмущаемого тем же шумом: причем условие ${{Y}_{0}} = 0$ общности измерений не ограничивает. Формально эти уравнения отличаются от уравнений (1.1) лишь независимостью функций $a( \cdot ),\,\,B( \cdot ),\,\,c( \cdot ),\,\,D( \cdot ),\,\,{{p}_{0}}( \cdot )$ от переменной выхода y. Поэтому и уравнения предлагаемого регулятора (1.3), (1.4) будем искать в соответствующем более простом виде: а его структурные функции найдем из условия минимума частного вида критерия (1.6): функции потерь которого $\varphi ( \cdot ),\psi ( \cdot )$ теперь также не зависят от переменной выхода.

В результате процедуру нахождения структуры регулятора (10.3), оптимальной в смысле критерия (10.4), получим из приведенных выше общих соотношений простым удалением из них переменной y.

Так из (2.9) имеем, что начальное состояние h определяется из алгебраического уравнения

В свою очередь функция смещения $f(t,z)$ либо обнуляется, согласно (5.2):

либо находится по такой модификации формулы (5.7): но уже с аналогично упрощенными операторами (2.5), (5.4): и с условной плотностью $\rho (t,x\,{\text{|}}\,z)$, которая при указанных в разд. 2 условиях гладкости удовлетворяет подобной версии уравнения (5.6):

Также и условно-оптимальную функцию усиления $\tilde {G}(t,z;u)$ найдем, согласно (7.1), (7.2), из линейного матричного уравнения

в котором по-прежнему $\Lambda = {{\psi }_{{zz}}} \otimes R,$ но теперь $H = c\,\psi _{z}^{{\text{T}}} + {{S}^{{\text{T}}}}{{\psi }_{{xz}}}$.

Наконец, функция выхода u(t, z) в общем случае ограниченности управления в соответствии с (6.4), (8.1) определяется как частный минимум

следующей функции трех аргументов:

Градиент последней по u, как следует из (8.2), (8.3), также имеет вид условного среднего ${{\mu }_{u}} = \int {\lambda \rho \,dx} $, но с более простой усредняемой функцией:

Тогда оптимальная функция усиления регулятора, согласно (6.5), получается простой подстановкой его функции выхода (10.6) в результат решения уравнения (10.5):

Отметим, что матрица S взаимной интенсивности зависимых шумов объекта (10.1) и измерителя (10.2) участвует в этих вычислениях только в составе вспомогательной матрицы H. В случае независимости этих шумов, как показано в разд. 1.1, ее следует положить равной нулю.

Заключение. Предложен способ синтеза быстрого, благодаря его конечномерности, дифференциального регулятора состояния стохастического объекта по его выходу, который учитывает всю предысторию измерений, но не требует восстановления состояния объекта. Регулятор может быть реализован в реальном времени на простом вычислителе и обладает наивысшей точностью в своем классе динамических регуляторов с конечной оперативной памятью, определяемой произвольно назначаемым порядком (размерностью) регулятора. Этот порядок можно выбрать из условия компромисса между достигаемой эффективностью управления и доступной скоростью пересчета измерений в управление. Действительно, указанная скорость определяется требуемым количеством вычислений, а последние в основном состоят в интегрировании полученных уравнений состояния регулятора.

При этом синтез структуры оптимального в смысле (1.6) конечномерного регулятора (1.3), (1.4) частично наблюдаемым стохастическим объектом (1.1) сводится к следующим операциям, которые в принципе тоже можно выполнить оперативно, в темпе со временем. Требуется решать прямую задачу Коши для ИДУ в ЧП (5.6) относительно условной плотности вероятности $\rho (t,x\,{\text{|}}\,y,z)$, через текущее сечение которой формулами (5.7) и (7.1) явным образом выражены структурные функции уравнения состояния регулятора (1.3), тогда как функция его выхода (1.4) выражается через него же в результате решения еще и параметрической задачи нелинейного программирования (6.4) для функции (8.1). Эти операции выгодно отличаются от выполняемых только заранее операций синтеза такого же регулятора при терминальном подходе [9], когда необходимо либо решить на всем отрезке времени управления ДТКЗ типа Лагранжа–Понтрягина, либо интегрировать в обратном времени аналог уравнения Беллмана для системы с распределенными параметрами.

Указанное упрощение достигается, естественно, за счет получения в конечный момент времени t = T несколько более плохого значения ${{I}_{T}}$ нестационарного критерия (1.6), чем в результате применения на всем отрезке $t \in [0,T]$ существенно более трудоемкого терминального подхода. Но этот проигрыш может оказаться не столь значительным, особенно если учесть расходы на построение самого регулятора, и может быть уменьшен определенным согласованием интегральных и терминальных слагаемых этих двух критериев.

Во второй части статьи планируется продемонстрировать применение полученных соотношений к системе со стохастической обратной связью, когда объект управления (10.1) и неточный измеритель его состояния (10.2) линейные, а критерий оптимальности (10.4) квадратично-биквадратный. Последний в подобной детерминированной задаче синтеза оперативно-оптимального позиционного регулятора приводит к кубической зависимости управления от состояния [10].