Журнал вычислительной математики и математической физики, 2021, T. 61, № 2, стр. 189-205

Многометодные алгоритмы для решения сложных задач оптимального управления

А. И. Тятюшкин^*

Институт динамики систем и теории управления им. В.М. Матросова СО РАН
664033 Иркутск, ул. Лермонтова, 134, Россия

Поступила в редакцию 19.12.2019
После доработки 20.08.2020
Принята к публикации 16.09.2020

DOI: 10.31857/S0044466921020137

Аннотация

Рассматриваются задачи оптимального управления с терминальными условиями без ограничений на управление, задачи со свободным правым концом траектории с ограничениями на управление и задачи оптимизации с параметрами при ограничениях на параметры и управление. Для каждого из этих классов задач конструируются многометодные алгоритмы, состоящие из наиболее эффективных для данных классов численных методов оптимального управления. Работа предложенных алгоритмов подтверждается численным решением сложных прикладных задач. Библ. 14. Фиг. 4.

Ключевые слова: многометодные алгоритмы оптимизации, оптимальное управление, задачи с параметрами, градиентные методы, принцип максимума, численные методы, итерации.

ВВЕДЕНИЕ

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

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

${{u}_{{{{i}_{0}}}}} = \mathop {argmax}\limits_{i \in \{ 1,2,3\} } \,(I(u_{i}^{k}) - I(u_{i}^{{k - 1}})).$

Затем это приближение передается всем трем методам для выполнения следующей итерации: $u_{i}^{{k + 1}} = {{u}_{{{{i}_{0}}}}}$, $i = 1,2,3$.

Фиг. 1.

Схема выполнения $(k + 1)$-й итерации многометодным алгоритмом для группы из трех методов ${{M}_{1}}$, ${{M}_{2}}$, ${{M}_{3}}$.

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

1. ЗАДАЧА БЕЗ ОГРАНИЧЕНИЙ НА УПРАВЛЕНИЕ

Рассмотрим сначала задачу оптимального управления (см. [1]) с условиями типа равенств и без ограничений на управление

(1.1)

$x = f(x,u,t),\quad t \in T = [{{t}_{0}},{{t}_{1}}],\quad x(t) \in {{R}^{n}},\quad u(t) \in {{R}^{r}},\quad x({{t}_{0}}) = {{x}^{0}},$

(1.2)

${{I}_{0}}(u) \to min,$

(1.3)

${{I}_{j}}(u) = 0,\quad j = \overline {1,m} ,$

где

${{I}_{j}}(u) = {{\varphi }^{j}}(x({{t}_{1}})) + \int_{{{t}_{0}}}^{{{t}_{1}}} {{{F}^{j}}(x,u,t)dt} ,\quad m \leqslant n.$

Градиенты функционалов (1.2), (1.3)

(1.4)

$\nabla {{I}_{j}}(u) = - H_{u}^{j}(\psi ,x,u,t),$

где (функция Понтрягина из [1]) ${{H}^{j}}(\psi ,x,u,t) = \psi _{j}^{'}(t)f(x,u,t) - {{F}^{j}}(x,u,t);$ ${{\psi }_{j}}(t)$ – решение сопряженной системы

(1.5)

${{\psi }_{j}} = - f_{x}^{'}(x,u,t){{\psi }_{j}} + F_{x}^{j}(x,u,t),\quad \psi ({{t}_{1}}) = - \psi _{x}^{j}(x({{t}_{1}})).$

Рассмотрим численный метод решения задачи (1.1)–(1.3), основанный на применении первой и второй вариаций. На первой фазе этого метода, когда итерационный процесс реализуется с сильным нарушением терминальных условий, минимизируется штрафной функционал

(1.6)

$J(u) = \varphi (x({{t}_{1}})) + \int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,{{F}^{0}}(x,u,t)dt,$

где

$\varphi (x({{t}_{1}})) = {{\varphi }^{0}}(x({{t}_{1}})) + \sum\limits_{j = 1}^m \,{{K}_{j}}{{[{{\varphi }^{j}}(x({{t}_{1}})) + {{x}_{{n + j}}}({{t}_{1}})]}^{2}},$

${{K}_{j}} \geqslant 0,$ ${{x}_{{n + j}}}({{t}_{1}}),$ $j = \overline {1,m} $ – решения дополнительных к системе (1.1) уравнений

${{\dot {x}}_{{n + j}}} = {{F}^{j}}(x,u,t),\quad {{x}_{{n + j}}}({{t}_{0}}) = 0,\quad j = \overline {1,m} .$

После прекращения сходимости первой фазы выполняется вторая фаза метода, на итерациях которой минимизируется исходный функционал (1.2), а вариация $\delta u(t)$ строится уже с учетом линеаризованных краевых условий.

Пусть теперь $H(\psi ,x,u,t) = \psi {\kern 1pt} {\text{'}}f(x,u,t) - {{F}^{0}}(x,u,t),$

(1.7)

$\dot {\psi } = - {{H}_{x}}(\psi ,x,u,t),\quad \psi ({{t}_{1}}) = - {{\varphi }_{x}}(x({{t}_{1}}))$

и для заданного управления ${{u}^{k}}(t)$ найдены решения ${{x}^{k}}$ и ${{\psi }^{k}}$ систем (1.1) и (1.7). Тогда проблему построения подходящей вариации $\delta {{u}^{k}}(t)$ сформулируем в виде следующей линейно-квадратичной задачи:

(1.8)

$I(\delta u) = \frac{1}{2}\delta x{\kern 1pt} '({{t}_{1}}){{\varphi }_{{xx}}}({{x}^{k}}({{t}_{1}}))\delta x({{t}_{1}}) - \int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,H_{u}^{'}\delta udt - \frac{1}{2}\int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,[\delta u{\kern 1pt} '{\kern 1pt} {{H}_{{uu}}}\delta u + 2\delta u{\kern 1pt} '{\kern 1pt} {{H}_{{ux}}}\delta x + \delta x{\kern 1pt} '{\kern 1pt} {{H}_{{xx}}}\delta x]{\kern 1pt} dt \to min,$

(1.9)

$\delta \dot {x} = {{f}_{x}}({{x}^{k}},{{u}^{k}},t)\delta x + {{f}_{u}}({{x}^{k}},{{u}^{k}},t)\delta u,\quad \delta x({{t}_{1}}) = 0.$

Здесь ${{H}_{u}}$, ${{H}_{{uu}}}$, ${{H}_{{ux}}}$, ${{H}_{{xx}}}$ – $r \times 1$-, $r \times r$-, $r \times n$-, $n \times n$-матрицы частных производных функции $H$, вычисленные на управлении ${{u}^{k}}(t)$ и траекториях ${{x}^{k}}(t)$, ${{\psi }^{k}}(t)$.

Для этой задачи построим гамильтониан

$\mathcal{H}(\delta \psi ,\delta x,\delta u,t) = \delta \psi {\kern 1pt} '{\kern 1pt} {{f}_{x}}\delta x + \delta \psi {\kern 1pt} '{\kern 1pt} {{f}_{u}}\delta u + H_{u}^{'}\delta u + \frac{1}{2}(\delta u{\kern 1pt} '{\kern 1pt} {{H}_{{uu}}}\delta u + 2\delta u{\kern 1pt} '{\kern 1pt} {{H}_{{ux}}}\delta x + \delta x{\kern 1pt} '{\kern 1pt} {{H}_{{xx}}}\delta x)$

и сопряженную систему

(1.10)

$\delta \dot {\psi } = - f_{x}^{'}\delta \psi - H_{{ux}}^{'}\delta u - {{H}_{{xx}}}\delta x,\quad \delta \psi ({{t}_{1}}) = - {{\varphi }_{{xx}}}\delta x({{t}_{1}}).$

Из условия ${{\mathcal{H}}_{{\delta u}}} = 0$ найдем решение вариационной задачи (1.8), (1.9):

(1.11)

$\delta u = - H_{{uu}}^{{ - 1}}(f_{u}^{'}\delta \psi + {{H}_{{ux}}}\delta x + {{H}_{u}}).$

Подставив последнюю формулу в уравнения (1.9), (1.10), получим следующую линейную двухточечную краевую задачу:

(1.12)

$\delta \dot {x} = C\delta x - {{f}_{u}}H_{{uu}}^{{ - 1}}f_{u}^{'}\delta \psi - {{f}_{u}}H_{{uu}}^{{ - 1}}{{H}_{u}},$

(1.13)

$\delta \dot {\psi } = - ({{H}_{{xx}}} - H_{{ux}}^{'}H_{{uu}}^{{ - 1}}{{H}_{{ux}}})\delta x - C{\kern 1pt} '{\kern 1pt} \delta \psi + {{H}_{{ux}}}H_{{uu}}^{{ - 1}}{{H}_{u}},$

где

(1.14)

$\begin{gathered} C = {{f}_{x}} - {{f}_{u}}H_{{uu}}^{{ - 1}}{{H}_{{ux}}}, \\ \delta x({{t}_{0}}) = 0,\quad \delta \psi ({{t}_{1}}) = - {{\varphi }_{{xx}}}\delta x({{t}_{1}}). \\ \end{gathered} $

Стандартный способ решения задачи состоит в применении формулы Коши, устанавливающей связь между краевыми условиями с помощью переходной матрицы $\Phi ({{t}_{0}},{{t}_{1}})$ размерности $2n \times 2n$:

$\left( {\begin{array}{*{20}{c}} {\delta x({{t}_{1}})} \\ {\delta \psi ({{t}_{1}})} \end{array}} \right) = \Phi ({{t}_{0}},{{t}_{1}})\left( {\begin{array}{*{20}{c}} {\delta x({{t}_{0}})} \\ {\delta \psi ({{t}_{0}})} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {\rho ({{t}_{1}})} \\ {\eta ({{t}_{1}})} \end{array}} \right),$

где $(\rho (t),\;\eta (t))$ – решение задачи Коши для системы (1.12), (1.13) при $\delta x({{t}_{1}}) = \delta {{x}_{0}} = 0$, $\delta \psi ({{t}_{0}}) = \delta {{\psi }_{0}} = 0$.

Матрицу $\Phi ({{t}_{0}},{{t}_{1}})$ разобьем на четыре равных блока и, учитывая равенство $\delta x({{t}_{0}}) = 0$, последнее уравнение перепишем в следующем виде:

(1.15)

$\begin{array}{*{20}{c}} {\delta x({{t}_{1}}) = {{\Phi }_{{12}}}\delta {{\psi }_{0}} + \rho ({{t}_{1}}),} \\ {\delta \psi ({{t}_{1}}) = {{\Phi }_{{22}}}\delta {{\psi }_{0}} + \eta ({{t}_{1}}).} \end{array}$

Подставляя краевое условие (1.14) в систему (1.15), получаем уравнение для определения начальных условий $\delta {{\psi }_{0}}$:

(1.16)

$({{\Phi }_{{22}}} + {{\varphi }_{{xx}}}{{\Phi }_{{12}}})\delta {{\psi }_{0}} = - {{\varphi }_{{xx}}}\rho ({{t}_{1}}) - \eta ({{t}_{1}}).$

Блоки ${{\Phi }_{{12}}}$ и ${{\Phi }_{{22}}}$ можно вычислить, проинтегрировав матричное уравнение

$\dot {\Phi } = A(t)\Phi ,\quad \Phi ({{t}_{0}},{{t}_{0}}) = E,$

где $A(t)$ есть $2n \times 2n$-матрица коэффициентов системы (1.12), (1.13).

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

Полагая $\delta x({{t}_{0}}) = 0,$ $\delta \psi ({{t}_{0}}) = {{e}^{i}},$ $i = \overline {1,n} $ (${{e}^{i}}$ – векторы стандартного ортонормированного базиса), $n$ раз интегрируем $2n$-мерную систему (1.12), (1.13). Каждый из полученных в результате интегрирования векторов возьмем в качестве $(n + i)$-го столбца матрицы $\Phi $; в результате получим ее блоки ${{\Phi }_{{12}}}$ и ${{\Phi }_{{22}}}$.

Решив систему линейных алгебраических уравнений (1.16), определим вектор $\delta {{\psi }_{0}}$. Затем, интегрируя систему (1.12), (1.13) в прямом времени, на решении ($\delta x(t),\;\delta \psi (t)$) по формуле (1.11) найдем искомую вариацию $\delta u(t)$. Формула (1.11) применима только при достаточно хорошей обусловленности матрицы ${{H}_{{uu}}}$. Чтобы сохранить вычислительную устойчивость метода для общего случая, вместо ${{H}_{{uu}}}$ на практике применяется матрица ${{H}_{{uu}}} + W$, где $W$ – положительно-определенная матрица. Найденная вариация используется для построения нового приближения ${{u}^{k}} + {{\alpha }_{k}}\delta {{u}^{k}}$ при ${{\alpha }_{k}} = \mathop {argmin}\limits_{\alpha \geqslant 0} \,J({{u}^{k}} + \alpha \delta {{u}^{k}})$. Итерационный процесс первой фазы метода прекращается при выполнении неравенства $J({{u}^{k}}) - J({{u}^{{k + 1}}}) \leqslant \varepsilon ,$ $\varepsilon \geqslant 0$. Поскольку при этом минимизировался штрафной функционал (1.6), требуемая точность выполнения краевых условий может оказаться не достигнутой, тогда выполняется вторая фаза метода, которая учитывает линеаризованные краевые условия (1.3) при решении задачи (1.12), (1.13).

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

(1.17)

${{\varphi }^{i}}(x({{t}_{1}})) = 0,\quad i = \overline {1,m} ,\quad m \leqslant n,$

или $\varphi (x({{t}_{1}})) = 0$, где $\varphi = ({{\varphi }^{1}},{{\varphi }^{2}}, \ldots ,{{\varphi }^{m}}){\kern 1pt} '$, а функционал

${{I}_{0}}(u) = {{\varphi }^{0}}(x({{t}_{1}})).$

Линеаризуем краевые условия (1.17) в окрестности точки ${{x}^{k}}({{t}_{1}})$:

(1.18)

$\varphi ({{x}^{k}}({{t}_{1}})) + {{\varphi }_{x}}({{x}^{k}}({{t}_{1}}))\delta x({{t}_{1}}) = 0.$

Часть компонент $\delta x_{i}^{f},i = \overline {1,m} $, вектора $\delta x({{t}_{1}})$ связана соотношениями (1.18), а оставшиеся $n - m$ свободных компонент должны удовлетворять терминальным условиям

(1.19)

$\delta {{\psi }^{c}}({{t}_{1}}) = - \varphi _{{xx}}^{0}({{x}^{k}}({{t}_{1}}))\delta {{x}^{0}}({{t}_{1}}).$

Разобьем введенные выше матрицы ${{\Phi }_{{12}}}$ и ${{\Phi }_{{22}}}$ на блоки в соответствии с составляющими векторов $\delta x({{t}_{1}}) = (\delta {{x}^{f}},\delta {{x}^{c}}),$ $\delta \psi ({{t}_{1}}) = (\delta {{\psi }^{f}},\delta {{\psi }^{c}})$ и перепишем уравнения (1.15) в следующем виде:

$\Phi _{{12}}^{{11}}\delta \psi _{0}^{f} + \Phi _{{12}}^{{12}}\delta \psi _{0}^{c} + {{\rho }^{f}} = \delta {{x}^{f}},$

$\Phi _{{12}}^{{21}}\delta \psi _{0}^{f} + \Phi _{{12}}^{{22}}\delta \psi _{0}^{c} + {{\rho }^{c}} = \delta x_{1}^{c},$

$\Phi _{{12}}^{{11}}\delta \psi _{0}^{f} + \Phi _{{12}}^{{12}}\delta \psi _{0}^{c} + {{\eta }^{f}} = \delta \psi _{1}^{f},$

$\Phi _{{12}}^{{21}}\delta \psi _{0}^{f} + \Phi _{{12}}^{{22}}\delta \psi _{0}^{c} + {{\eta }^{c}} = \delta \psi _{1}^{c}.$

Подставив в эти уравнения краевые условия (1.18), (1.19), получим систему линейных алгебраических уравнений относительно вектора начальных условий $\delta {{\psi }_{0}} = (\delta \psi _{0}^{f},\delta \psi _{0}^{c}){\kern 1pt} :$

(1.20)

$\begin{gathered} (\Phi _{{22}}^{{21}} + \varphi _{{xx}}^{0}\Phi _{{12}}^{{21}})\delta \psi _{0}^{f} + (\Phi _{{22}}^{{22}} + \varphi _{{xx}}^{0}\Phi _{{12}}^{{22}})\delta \psi _{0}^{c} = - \varphi _{{xx}}^{0}{{\rho }^{0}} - {{\eta }^{c}}, \\ ({{\varphi }_{{{{x}^{f}}}}}\Phi _{{12}}^{{11}} + {{\varphi }_{{{{x}^{c}}}}}\Phi _{{12}}^{{21}})\delta \psi _{0}^{f} + ({{\varphi }_{{{{x}^{f}}}}}\Phi _{{12}}^{{12}} + {{\varphi }_{{{{x}^{c}}}}}\Phi _{{12}}^{{22}})\delta \psi _{0}^{c} = - {{\varphi }_{x}}\rho - \varphi , \\ \end{gathered} $

где

$\varphi _{{xx}}^{0} = \varphi _{{xx}}^{0}({{x}^{k}}({{t}_{1}})),\quad {{\varphi }_{x}} = ({{\varphi }_{{{{x}^{f}}}}},{{\varphi }_{{{{x}^{c}}}}}) = {{\varphi }_{x}}({{x}^{k}}({{t}_{1}})),\quad \varphi = \varphi ({{x}^{k}}({{t}_{1}})).$

Следовательно, вариация $\delta u(t)$, построенная на решении $(\delta x(t),\delta \psi (t))$ системы (1.12), (1.13) с начальными условиями $(0,\delta {{\psi }_{0}})$, будет выполнять в линейном приближении краевые условия (1.17) и одновременно минимизировать квадратичное приближение функционала (1.2).

Алгоритм 1. Приведем вычислительную схему, включающую обе фазы описанного метода. На итерациях первой фазы выполняются следующие операции.

1. При заданном управлении ${{u}^{k}}(t)$ интегрируется уравнение (1.1), в узлах интегрирования запоминается траектория ${{x}^{k}}(t)$.

2. В обратном времени интегрируется уравнение (1.7), в узлах интегрирования запоминается решение ${{\psi }^{k}}(t)$.

3. С коэффициентами, вычисленными на управлении ${{u}^{k}}(t)$ и решениях ${{x}^{k}}(t)$, ${{\psi }^{k}}(t)$, $n + 1$ раз интегрируется система (1.12), (1.13) при начальных условиях

$\delta x({{t}_{0}}) = 0,\quad \delta \psi ({{t}_{0}}) = {{e}^{i}},\quad i = \overline {1,n} ;\quad \delta x({{t}_{0}}) = 0,\quad \delta \psi ({{t}_{0}}) = 0;$

$(n + j)$-е столбцы матрицы $\Phi $ принимаются равными $(\delta x{\kern 1pt} '({{t}_{1}}),\delta \psi {\kern 1pt} '({{t}_{1}})){\kern 1pt} ',$ $j = \overline {1,n} $.

4. С помощью блоков ${{\Phi }_{{12}}}$ и ${{\Phi }_{{22}}}$ матрицы $\Phi $ формируется система (1.16) и находится ее решение $\delta {{\psi }^{0}}$.

5. Интегрируется система (1.12), (1.13) при $\delta x({{t}_{0}}) = 0,$ $\delta \psi ({{t}_{0}}) = \delta {{\psi }^{0}}$, в каждом узле интегрирования по формуле (1.11) вычисляется и запоминается вариация $\delta {{u}^{k}}(t)$.

6. С помощью процедуры одномерного поиска за $\nu $ задач Коши (1.1) находится параметр ${{\alpha }_{k}} = \mathop {argmin}\limits_{\alpha \geqslant 0} \,J({{u}^{k}} + \alpha \delta {{u}^{k}})$.

7. Строится новое приближение ${{u}^{{k + 1}}}(t) = {{u}^{k}}(t) + {{\alpha }_{k}}\delta {{u}^{k}}(t),$ $t \in T$, и если $J({{u}^{k}}) - J({{u}^{{k + 1}}}) > \varepsilon $, то повторяется итерация при $k = k + 1$ с п. 1. В противном случае выполняется вторая фаза алгоритма.

Основное отличие второй фазы алгоритма состоит в том, что вместо системы (1.16) для нахождения вектора $\delta {{\psi }^{0}}$ применяется система (1.20). На итерациях второй фазы выполняются те же операции 1–7, кроме 4, где вместо (1.16) формируется система (1.20) и операции 2, где вместо (1.7) интегрируется уравнение (1.5) при $j = 0$.

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

Графически убывание функционала $I(u)$ на итерациях многометодного алгоритма будет изображаться ломаной линией, составленной из графиков отдельных методов. На фиг. 2 показана работа многометодного алгоритма в случае использования двух методов ${{M}_{1}}$ и ${{M}_{2}}$. Приведенные здесь графики показывают убывание функционала на итерациях этих методов. Затем из них составлен график убывания функционала на итерациях многометодного алгоритма – кривая $ABC$, участок $BC$ которой получен параллельным переносом кривой $EL$. Согласно фиг. 2, нулевое значение функционала методом ${{M}_{1}}$ достигается за ${{k}_{1}}$ итераций, а методом ${{M}_{2}}$ – за ${{k}_{2}}$ итераций. Многометодный алгоритм до $\bar {k}$-й итерации работает по методу ${{M}_{1}}$ (кривая $AB$), а далее – по методу ${{M}_{2}}$ (кривая $BC$), так как, начиная с $\bar {k}$, скорость убывания функционала по методу ${{M}_{2}}$ будет выше. В результате нулевое значение функционала многометодным алгоритмом достигается за $k{\kern 1pt} *$ итераций, т.е. за существенно меньшее число итераций, чем каждым из методов ${{M}_{1}}$ и ${{M}_{2}}$.

Фиг. 2.

Убывание функционала на итерациях методов ${{M}_{1}}$, ${{M}_{2}}$ и l многометодного алгоритма.

2. ЗАДАЧИ С ОГРАНИЧЕНИЯМИ НА УПРАВЛЕНИЕ

Для решения наиболее важного с прикладной точки зрения класса задач с ограничениями на управление

(2.1)

$u(t) \in U,\quad t \in T,$

где $U$ – компактное множество из ${{R}^{r}}$, построим многометодные алгоритмы, основанные на принципе максимума (см. [3], [4]) и методах градиентного типа (см. [5]–[8]).

Будем предполагать, что с помощью введения штрафного функционала уже имеем задачу со свободным правым концом и что при некотором допустимом управлении ${{u}^{k}}(t) \in U,$ $t \in T$, найдено решение уравнения (1.1) ${{x}^{k}}(t)$. Решая уравнение (1.5) при $u = {{u}^{k}}(t),$ $x = {{x}^{k}}(t),$ $j = 0$, найдем ${{\psi }^{k}} = {{\psi }_{0}}(t)$ и вычислим

(2.2)

${{\bar {u}}^{k}}(t) = \mathop {argmax}\limits_{u \in U} \,H({{\psi }^{k}},{{x}^{k}},u,t),\quad t \in T.$

Построим скалярную функцию

(2.3)

${{w}_{k}}(u(t),t) = H({{\psi }^{k}},{{x}^{k}},u,t) - H({{\psi }^{k}},{{x}^{k}},{{u}^{k}},t),$

которая при $u = \bar {u}(t)$, очевидно, удовлетворяет неравенству

(2.4)

${{w}_{k}}({{\bar {u}}^{k}}(t),t) \geqslant 0,\quad t \in T.$

Пусть ${{\tau }_{k}} \in T$ – точка максимума функции ${{w}_{k}}({{\bar {u}}^{k}}(t),t)$:

(2.5)

${{w}_{k}}({{\bar {u}}^{k}}({{\tau }_{k}}),{{\tau }_{k}}) = \mathop {max}\limits_{t \in T} {{w}_{k}}({{\bar {u}}^{k}}(t),t).$

Тогда необходимое условие первого порядка (принцип максимума (см. [1], [2])) формулируется так: если ${{u}^{k}}(t)$ – оптимальное управление в задаче (1.2), (1.1), (2.1), то

(2.6)

${{w}_{k}}({{\bar {u}}^{k}}({{\tau }_{k}}),{{\tau }_{k}}) = 0.$

Предположим, что для заданного ${{u}^{k}}(t) \in U$ и найденных ${{x}^{k}}(t),\;{{\psi }^{k}}(t),\;{{\bar {u}}^{k}}(t)$ условие (2.6) не выполняется:

${{w}_{k}}({{\bar {u}}^{k}}({{\tau }_{k}}),{{\tau }_{k}}) > 0.$

Тогда можно найти новое управление, на котором значение функционала (1.2) будет меньше, чем ${{I}_{0}}({{u}^{k}})$.

Построим отрезок $T_{\varepsilon }^{k} \subseteq T$ по следующему правилу:

(2.7)

$T_{\varepsilon }^{k} = [{{\tau }_{k}} - \varepsilon ({{\tau }_{k}} - t_{0}^{k}),{{\tau }_{k}} + \varepsilon (t_{1}^{k} - {{\tau }_{k}})],\quad \varepsilon \in [0,1],$

где $t_{0}^{k}$ и $t_{1}^{k}$ – ближайшие слева и справа точки разрыва функции ${{w}_{k}}({{\bar {u}}^{k}}(t),t)$. Этот отрезок с мерой $\operatorname{mes} T_{\varepsilon }^{k} = {{\alpha }_{k}}(\varepsilon ) = {{\alpha }_{k}}\varepsilon ,$ $0 \leqslant {{a}_{k}} \leqslant {{t}_{1}} - {{t}_{0}},$ $\varepsilon \in [0,1]$, обладает следующими свойствами:

1) ${{\alpha }_{k}}(\varepsilon ) \to 0$, когда $\varepsilon \to 0$;

2) при $\varepsilon \to 0$ он стягивается в точку ${{\tau }_{k}}$;

3) при всех $\varepsilon \in [0,1]$ функция ${{w}_{k}}({{\bar {u}}^{k}}(t),t)$ непрерывна на ${{T}_{\varepsilon }}$.

Далее находим параметр

(2.8)

${{\varepsilon }_{k}} = \mathop {argmin}\limits_{\varepsilon \in [0,1]} \,{{I}_{0}}(u_{\varepsilon }^{k}),$

где

(2.9)

$u_{\varepsilon }^{k} = \left\{ {\begin{array}{*{20}{l}} {{{{\bar {u}}}^{k}}(t),\quad t \in {{T}_{\varepsilon }},} \\ {{{u}^{k}}(t),\quad t \in T{\backslash }{{T}_{\varepsilon }},} \end{array}} \right.$

и определяем новое приближение

(2.10)

${{u}^{{k + 1}}}(t) = u_{{{{\varepsilon }_{k}}}}^{k}(t),\quad t \in T,\quad k = 0,1, \ldots .$

Алгоритм 2. Построим вычислительную схему описанного метода, сходимость которого доказана в [4].

1. Задаем граничное управление ${{u}^{0}}(t),t \in T$; полагаем $k = 1$.

2. В прямом времени интегрируем систему (1.1) при $u = {{u}^{k}}(t)$, запоминая в узлах интегрирования ${{x}^{k}}(t)$.

3. В обратном времени интегрируем сопряженную систему (1.5). При этом в каждом узле интегрирования находим и запоминаем управление (2.2), вычисляем значение функции ${{w}_{k}}({{\bar {u}}^{k}}(t),t)$, определяем точку максимума ${{\tau }_{k}}$.

4. Если ${{w}_{k}}({{\bar {u}}^{k}}({{\tau }_{k}}),{{\tau }_{k}}) \leqslant \varepsilon $, то процесс прекращаем. В противном случае выполним операцию 5.

5. С помощью процедуры одномерного поиска решаем задачу (2.8), (2.9), (2.7). Поскольку в точках $t \in T$, удовлетворяющих неравенству $t < {{\tau }_{k}} - \varepsilon ({{\tau }_{k}} - t_{0}^{k})$, управление (2.9) будет равным ${{u}^{k}}(t)$, интегрирование системы (1.1) следует начинать с ближайшего к ${{\tau }_{k}} - \varepsilon ({{\tau }_{k}} - t_{0}^{k})$ левого узла ${{t}_{k}}$, в котором $x({{t}_{k}}) = {{x}^{k}}({{t}_{k}})$ было найдено в п. 2. Учитывая структуру управления (2.9), в качестве начальной точки можно брать также “наибольший” узел ${{t}_{k}} \in \{ {{t}_{k}} \in T_{\varepsilon }^{k}:{{u}^{k}}({{t}_{k}}) = {{\bar {u}}^{k}}({{t}_{k}})\} $. Так как при поиске ${{\varepsilon }_{k}}$ система (1.1) интегрируется несколько раз, то такой выбор начальной точки может существенно сократить время решения задачи.

6. Если при найденном ${{\varepsilon }_{k}}$, точности $\delta $ вычисления ${{I}_{0}}({{u}^{k}})$, шаге интегрирования $h$ выполняются неравенства ${{I}_{0}}(u_{{{{\varepsilon }_{k}}}}^{k}) \geqslant {{I}_{0}}({{u}^{k}}) - \delta $ и $t_{1}^{s} - t_{0}^{s} > h$, то сокращаем отрезок $T_{\varepsilon }^{k}$, полагая $\varepsilon = {{2}^{{ - s}}}$ ($s$ – число сокращений $T_{\varepsilon }^{k}$), и вновь выполняем операцию 5.

Если ${{I}_{0}}(u_{{{{\varepsilon }_{k}}}}^{k}) < {{I}_{0}}({{u}^{k}}) - \delta $, то, приняв ${{u}^{{k + 1}}} = u_{{\varepsilon k}}^{k},$ $k = k + 1$, повторяем цикл с п. 2. В противном случае итерационный процесс прекращается.

Поскольку итерации алгоритма идут в классе кусочно-постоянных управлений, на полученном управлении условие оптимальности (2.8) может не выполняться с заданной точностью (хотя величина $\int_T {{{w}_{k}}({{u}^{k}}(t),t)dt} $ будет достаточно мала). В отличие от градиентных методов данный алгоритм применим и к таким задачам оптимального управления, в которых вектор-функция $f(x,u,t)$ не дифференцируема по $u$, а множество $U$ не является выпуклым или даже связным.

Как уже отмечалось в [4], [8], [9], алгоритмы принципа максимума нередко приводят к “прилипанию” управления к границе, что является причиной ухудшения сходимости. Этот эффект обусловлен тем, что в некоторых системах (например, линейных по управлению) решение задачи (1.2) достигается на границе и, следовательно, приближения по управлению вида (2.9), (2.10) также имеют граничные значения. При таком приближении сходимость алгоритма на последних итерациях обеспечивается очень малым шагом дискретизации и влечет большие затраты процессорного времени. Восстановить или улучшить сходимость итерационного процесса в этой ситуации часто удается более простым способом – построением выпуклой комбинации управлений ${{u}^{k}}(t)$ и ${{\bar {u}}^{k}}(t)$:

(2.11)

${{u}^{{k + 1}}}(t) = {{u}^{k}}(t) + {{\alpha }_{k}}[{{\bar {u}}^{k}}(t) - {{u}^{k}}(t)],\quad {{\alpha }_{k}} \in [0,1],$

на отрезке $T_{\varepsilon }^{k}$.

Причиной прекращения сходимости может служить и то обстоятельство, что при сокращении отрезка $T_{\varepsilon }^{k}$ малое число узлов интегрирования, оставшихся в нем, не обеспечивает построение такой вариации управления, на которой приращение функционала достигает значения, большего (по модулю) величины погрешностей интегрирования.

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

(2.12)

${{T}_{\varepsilon }} = \{ t \in T:{{w}_{k}}({{\bar {u}}^{k}}(t),t) \geqslant \varepsilon {{w}_{k}}({{\bar {u}}^{k}}({{t}_{k}}),{{\tau }_{k}})\} ,\quad \varepsilon \in [0,1].$

Множество $T$ включает все точки $t \in T$ нарушения принципа максимума и состоит из нескольких непересекающихся отрезков, если функция ${{w}_{k}}({{\bar {u}}^{k}}(t),t)$ многоэкстремальна. Варьируя $\varepsilon $, можно найти такое значение ${{\varepsilon }_{k}}$, при котором управление (2.9) доставит наименьшее значение функционалу. Применение множества (2.12) вместо отрезка (2.7) во многих задачах улучшает сходимость алгоритма 2. В случае, когда максимум ${{w}_{k}}({{\bar {u}}^{k}}(t),t)$ достигается на множестве ${{T}_{{{{\varepsilon }_{k}}}}}$ с положительной мерой, возможно прекращение сходимости алгоритма. Тогда итерации алгоритма 2 продолжаются с выбором $T_{\varepsilon }^{k}$ по формуле (2.7).

Алгоритм 3. Таким образом, итоговая вычислительная схема для построения нового приближения будет иметь следующий вид.

1. Интегрируется уравнение (1.1) при $u = {{u}^{k}}(t)$, в узлах интегрирования запоминается траектория ${{x}^{k}}(t)$.

2. В обратном времени интегрируется сопряженная система (1.5), в каждом узле интегрирования вычисляются и запоминаются управление ${{\bar {u}}^{k}}(t)$ и скалярная функция

${{w}_{k}}(t) = {{w}_{k}}({{\bar {u}}^{k}}(t),t).$

3. Находится точка ${{\tau }_{k}} = \mathop {argmax}\limits_{t \in T} \,{{w}_{k}}(t)$. Если ${{w}_{k}}({{\tau }_{k}}) \leqslant \varepsilon $, то итерационный процесс прекращается.

4. Решается задача одномерного поиска $I(u_{\varepsilon }^{k}) \to min$. При этом для каждого $\varepsilon \in [0,1]$, используемого методом одномерного поиска, интегрирование уравнения (1.1) начинается с узла ${{t}_{k}}$, в котором впервые выполняется неравенство ${{w}_{k}}(t) > \varepsilon {{w}_{k}}({{\tau }_{k}})$, а управление выбирается по формуле

$u_{\varepsilon }^{k}(t) = \left\{ {\begin{array}{*{20}{l}} {{{{\bar {u}}}^{k}}(t),\quad {\text{если}}\quad {{w}_{k}}(t) \geqslant \varepsilon {{w}_{k}}({{\tau }_{k}}),} \\ {{{u}^{k}}(t),\quad {\text{если}}\quad {{w}_{k}}(t) < \varepsilon {{w}_{k}}({{\tau }_{k}}),\quad t \in T.} \end{array}} \right.$

5. Если ${{I}_{0}}(u_{{{{\varepsilon }_{k}}}}^{k}) \geqslant {{I}_{0}}({{u}^{k}}) - \delta $ и $\operatorname{mes} {{T}_{{{{\varepsilon }_{k}}}}} > h$, выполняется операция 5 алгоритма 2, где $T_{\varepsilon }^{k}$ строится по формуле (2.7). При этом в качестве первого приближения берется ${{T}_{{{{\varepsilon }_{k}}}}}$.

6. Если итерация алгоритма 2 тоже не улучшает управление ${{u}^{k}}(t)$, то на полученном в п. 4 отрезке $T_{{{{\varepsilon }_{k}}}}^{k}$ выполняется итерация метода условного градиента (2.11), где для определения ${{\alpha }_{k}}$ снова применяется одномерный поиск.

7. Если в п. 4–6 будет найдено значение ${{\varepsilon }_{k}}$ или ${{\alpha }_{k}}$, при которых новое приближение для управления обеспечивает меньшее значение функционалу, то выполняется новая итерация (при $k: = k + 1$) с п. 1. В противном случае итерационный процесс прекращается.

Заметим, что данный алгоритм, как и описанный выше, ориентирован на решение задач с ограничениями на управление, но со свободным правым концом.

3. ОБЩАЯ ЗАДАЧА ОПТИМАЛЬНОГО УПРАВЛЕНИЯ С ПАРАМЕТРАМИ

Рассмотрим теперь более общую задачу оптимального управления – с фазовыми ограничениями, и когда правая часть системы зависит не только от управлений, но и от параметров. Начальные условия системы также могут зависеть от параметров, и их выбором обычно обеспечивается, например, оптимальный “старт” процесса.

Для решения этой сложной задачи применим сначала редукцию к конечномерной задаче, а затем построим многометодный алгоритм поиска оптимального управления.

Пусть задан управляемый процесс, зависящий от параметров,

(3.1)

$\begin{gathered} \dot {x} = f(x,u,w,t),\quad x(t) \in {{E}^{n}},\quad u(t) \in {{E}^{r}},\quad t \in T = [{{t}_{0}},{{t}_{1}}], \\ x({{t}_{0}}) = \Theta (v),\quad w \in {{R}^{p}},\quad v \in {{R}^{n}}, \\ \end{gathered} $

с терминальными условиями

(3.2)

${{I}_{i}}(u) = {{h}_{i}}(x({{t}_{1}})) = 0,\quad i = \overline {1,m} ,$

и фазовыми ограничениями

(3.3)

${{J}_{i}}(u,v) = {{g}_{i}}(x(t),t) = 0,\quad t \in T,\quad i = \overline {1,s} .$

Управление и параметры стеснены следующими ограничениями:

(3.4)

${{c}_{i}}(u,t) = 0,\quad t \in T,\quad i = \overline {1,l} ,$

(3.5)

${{u}^{{\text{н}}}}(t) \leqslant u(t) \leqslant {{u}^{{\text{в}}}}(t),\quad t \in T,$

(3.6)

${{v}^{{\text{н}}}} \leqslant v \leqslant {{v}^{{\text{в}}}},\quad {{w}^{{\text{н}}}} \leqslant w \leqslant {{w}^{{\text{в}}}},$

где ${{c}_{i}}(u,t),\;i = \overline {1,l} $ – непрерывно дифференцируемые по $u$ и кусочно-непрерывные по $t$ функции; $\Theta (v)$ – непрерывно дифференцируемая вектор-функция. Относительно функций, определяющих условия (3.1)–(3.3), справедливы предположения, оговоренные ранее, к которым добавляется также их непрерывная дифференцируемость по параметрам.

Требуется среди управлений и параметров, удовлетворяющих ограничениям (3.4)–(3.6), найти такие, которые обеспечивают выполнение условий (3.3) для управляемого процесса (3.1) и приводят его в точку фазового пространства, где с заданной точностью будут выполнены условия (3.2), а функционал

(3.7)

${{I}_{0}}(u) = \varphi (x({{t}_{1}}))$

достигнет наименьшего значения.

3.1. Редукция к конечномерной задаче

Для построения конечномерной задачи на заданном интервале $T$ вводится сетка дискретизации с узлами ${{t}_{0}},{{t}^{1}}, \ldots ,{{t}^{N}}$ такими, что

(3.8)

${{t}_{0}} = {{t}^{0}} < {{t}^{1}} < \ldots < {{t}^{N}} = {{t}_{1}}.$

Эта сетка может быть и неравномерной.

Управляющие функции ${{u}^{i}}(t),\;i = \overline {1,r} $, ищутся только в узлах (3.8), а для получения промежуточных значений ${{u}^{i}}(t),\;i = \overline {1,r} $, используется либо кусочно-постоянная аппроксимация

${{u}^{i}}(t) = {{u}^{i}}({{t}^{j}}) = u_{j}^{i},\quad t \in [{{t}^{j}},{{t}^{{j + 1}}}],$

либо кусочно-линейная

(3.9)

${{u}^{i}}(t) = [({{t}^{{j + 1}}} - t){{u}_{j}} + (t - {{t}^{j}}){{u}_{{j + 1}}}]({{t}^{{j + 1}}} - {{t}^{j}}),\quad t \in [{{t}^{j}},{{t}^{{j + 1}}}].$

Тогда конечномерная задача, аппроксимирующая задачу (3.1)–(3.7), будет иметь следующий вид:

$\dot {x} = f(x,u,w,t),\quad t \in T = [{{t}_{0}},{{t}_{1}}],\quad x({{t}_{0}}) = \Theta (v),$

${{h}_{i}}(x({{t}^{N}})) = 0,\quad i = \overline {1,m} ,$

${{g}_{i}}(x({{t}^{j}}),{{t}^{j}}) = 0,\quad i = \overline {1,s} ,\quad j = \overline {0,N} ,$

(3.10)

${{c}_{i}}({{u}_{j}},{{t}^{j}}) = 0,\quad i = \overline {1,l} ,\quad j = \overline {0,N} ,$

${{v}^{{\text{н}}}} \leqslant v \leqslant {{v}^{{\text{в}}}},\quad {{w}^{{\text{н}}}} \leqslant w \leqslant {{w}^{{\text{в}}}},$

$\varphi (x({{t}^{N}})) \to min,\quad u_{j}^{{\text{н}}} \leqslant {{u}_{j}} \leqslant u_{j}^{{\text{в}}},\quad j = \overline {0,N} ,$

где

$u_{j}^{{\text{н}}} = {{u}^{{\text{н}}}}({{t}^{j}}),\quad u_{j}^{{\text{в}}} = {{u}^{{\text{в}}}}({{t}^{j}}),\quad j = \overline {0,N} .$

Заметим, что в аппроксимирующей задаче (3.10) управляемый процесс (3.1) остается непрерывным, а в процессе счета он с требуемой точностью (достаточно высокой) моделируется численным методом интегрирования.

3.2. Численное решение конечномерной задачи

Градиенты функционалов ${{I}_{j}}(u),\;j = \overline {0,m} $, с помощью функций ${{H}^{j}}({{\psi }_{j}},x,u,t) = \psi _{j}^{'}(t)f(x,u,t)$ и сопряженной системы

${{\dot {\psi }}_{j}} = - {{f}_{x}}(x,u,t){\kern 1pt} '{{\psi }_{j}}(t),\quad {{\psi }_{j}}({{t}_{1}}) = - \varphi _{x}^{j}(x({{t}_{1}})$

традиционно определяются по формулам

$\nabla {{I}_{j}}(u) = - H_{u}^{j}({{\psi }_{j}},x,u,t),\quad j = \overline {0,m} .$

Для каждого $t \in T$ можно аналогично вычислить градиенты ${{J}_{j}}(u,t),\;j = \overline {1,s} $:

$\nabla {{I}_{j}}(u,t) = - \bar {H}_{u}^{j}({{\Phi }_{j}},x,u,t,\tau ),\quad {{t}_{0}} \leqslant \tau \leqslant t \leqslant {{t}_{1}},$

где $\bar {H}_{{}}^{j}({{\Phi }_{j}},x,u,t,\tau ) = \Phi _{j}^{'}(t,\tau )f(x,u,\tau ),{{\Phi }_{j}}(t,\tau ),$ $j = \overline {1,s} $ – решения сопряженной системы

$\frac{{\partial {{\Phi }_{j}}(t,\tau )}}{{\partial \tau }} = - \frac{{\partial f(x,u,\tau )}}{{\partial x}}{{\Phi }_{j}}(t,\tau ),\quad \tau \in [{{t}_{0}},t],$

с краевыми условиями

${{\Phi }_{j}}(t,t) = - \frac{{\partial {{g}^{j}}(x(t))}}{{\partial x}},\quad j = \overline {1,s} .$

Линеаризуем ограничения в аппроксимирующей задаче. Матрица-якобиан линеаризованных составляется из градиентов $\nabla {{I}_{i}},i = \overline {1,m} $, и $\nabla {{J}_{j}}(t),\;j = \overline {1,s} ,$ $t \in T$. Так как правые части и начальные условия системы (3.1) зависят еще и от параметров, то необходимо иметь также градиенты функционалов ${{I}_{i}},i = \overline {1,m} $, и ${{J}_{j}}(t),\;j = \overline {1,s} ,$ $t \in T$, по этим параметрам (см. [3], [8]):

(3.11)

$\begin{gathered} {{\nabla }_{v}}{{I}_{i}}({{u}^{k}},{{w}^{k}},{{v}^{k}}) = - {{\psi }_{i}}({{t}_{0}}){\kern 1pt} '{{\Theta }_{v}}({{v}^{k}}),\quad i = \overline {1,m} , \\ {{\nabla }_{w}}{{I}_{i}}({{u}^{k}},{{w}^{k}},{{v}^{k}}) = - \int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,{{\psi }_{i}}(t){\kern 1pt} '{\kern 1pt} {{f}_{w}}({{x}^{k}},{{u}^{k}},{{w}^{k}},t)dt, \\ \end{gathered} $

(3.12)

${{\nabla }_{w}}{{J}_{i}}({{u}^{k}},{{w}^{k}},{{v}^{k}},{{t}^{j}}) = - \int\limits_{{{t}_{0}}}^{{{t}^{j}}} \,{{\Phi }_{i}}(t){\kern 1pt} '{\kern 1pt} {{f}_{w}}({{x}^{k}},{{u}^{k}},{{w}^{k}},t)dt,$

(3.13)

${{\nabla }_{v}}{{J}_{i}}({{u}^{k}},{{w}^{k}},{{v}^{k}},{{t}^{j}}) = - {{\Phi }_{i}}({{t}_{0}}){\kern 1pt} '\Theta ({{v}^{k}}),\quad i = \overline {1,s} ,\quad j = \overline {1,N} .$

Пусть теперь на $k$-й итерации внешнего метода на сетке (3.8) найдено ${{u}^{k}}({{t}^{j}})$ и соответствующее ему ${{x}^{k}}({{t}^{j}}),\;j = \overline {1,N} $. Для расчета градиентов по управлению ${{\nabla }_{u}}{{I}_{i}},i = \overline {1,m} $, система (3.10) $m$ раз интегрируется от ${{t}_{1}}$ до ${{t}_{0}}$ с разными начальными условиями. Попутно вычисляются градиенты (3.11) с использованием квадратурных формул для расчета интегралов. Далее ищутся градиенты функционалов ${{J}_{i}}({{t}^{j}}),\;j = \overline {1,N} ,$ $i = \overline {1,s} $. Для этого нужно $s$ раз решить задачу Коши для каждого узла сетки (3.8), т.е. проинтегрировать систему $s \cdot N$ раз в среднем на половине отрезка $T$.

На полученных решениях вычисляются компоненты градиентов ${{\nabla }_{u}}{{I}_{i}},i = \overline {1,m} $, и ${{\nabla }_{u}}{{J}_{i}}({{t}^{j}}),\;i = \overline {1,s} ,$ $j = \overline {1,N} $; с учетом аппроксимации управления их значения равны

$\int\limits_{{{t}^{j}}}^{{{t}^{{j + 1}}}} \,{{\psi }^{k}}(t){\kern 1pt} '{\kern 1pt} {{f}_{u}}({{x}^{k}},u_{j}^{k},{{w}^{k}},t)dt$

в случае кусочно-постоянной аппроксимации и

(3.14)

$\frac{1}{{{{t}^{{j + 1}}} - {{t}^{j}}}}\left[ {\int\limits_{{{t}^{{j - 1}}}}^{{{t}^{j}}} \,{{\psi }^{k}}(t){\kern 1pt} '{\kern 1pt} {{f}_{u}}({{x}^{k}},{{{\bar {u}}}^{k}}(t),{{w}^{k}},t)(t - {{t}^{{j - 1}}})dt} \right. + \left. {\int\limits_{{{t}^{j}}}^{{{t}^{{j + 1}}}} \,{{\psi }^{k}}(t){\kern 1pt} '{\kern 1pt} {{f}_{u}}({{x}^{k}},{{{\bar {u}}}^{k}}(t),{{w}^{k}},t)({{t}^{{j + 1}}} - t)dt} \right]$

в случае кусочно-линейной аппроксимации (3.9). При этом ${{\bar {u}}^{k}}(t)$ вычисляется по формуле (3.9) при ${{u}_{j}} = u_{j}^{k},$ ${{u}_{{j + 1}}} = u_{{j + 1}}^{k}$.

Из полученных значений градиентов по управлению$\nabla {{I}_{i}},i = \overline {1,m} $, и $\nabla {{J}_{j}}(t)$, $j = \overline {1,s} ,$ $t \in T$, и вычисленных по формулам (3.11)–(3.13) градиентов по параметрам составляется матрица коэффициентов линеаризованных ограничений. Она дополняется также блоком элементов $\partial {{c}_{i}}{\text{/}}\partial {{u}_{j}}$, соответствующим ограничениям на управление, и приобретает вид матрицы специальной блочной структуры, которую обозначим через $A$.

3.3. Алгоритм метода приведенного градиента из [6]

Вводя векторные обозначения для равенств (3.2)–(3.4), построим модифицированную функцию Лагранжа (см. [7]) для задачи (3.1)–(3.7):

(3.15)

$\begin{gathered} L = \varphi (x({{t}_{1}})) - {{\lambda }^{{k{\kern 1pt} '}}}[h(x({{t}_{1}})) - {{{\bar {h}}}^{L}}] + \frac{\rho }{2}[h(x({{t}_{1}})) - {{{\bar {h}}}^{L}}]{\kern 1pt} '[h(x({{t}_{1}})) - {{{\bar {h}}}^{L}}] - \int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,{{\mu }^{{k{\kern 1pt} '}}}(t)[g(x(t),t) - {{{\bar {g}}}^{L}}]{\kern 1pt} dt + \\ + \;\frac{\rho }{2}\int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,[g(x(t),t) - {{{\bar {g}}}^{L}}]{\kern 1pt} '[g(x(t),t) - {{{\bar {g}}}^{L}}]{\kern 1pt} dt - \int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,{{\gamma }^{k}}(t)[c(u,t) - {{{\bar {c}}}^{L}}]{\kern 1pt} dt + \frac{\rho }{2}\int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,[c(u,t) - {{{\bar {c}}}^{L}}]{\kern 1pt} '[c(u,t) - {{{\bar {c}}}^{L}}]{\kern 1pt} dt, \\ \end{gathered} $

где

${{\bar {h}}^{L}} = h({{x}^{k}}({{t}_{1}})) + {{h}_{x}}({{x}^{k}}({{t}_{1}}))\delta x({{t}_{1}}),\quad {{\bar {g}}^{L}} = g({{x}^{k}}(t),t) + {{g}_{x}}({{x}^{k}}(t),t)\delta x(t),$

${{\bar {c}}^{L}} = c({{u}^{k}}(t),t) + {{c}_{u}}({{u}^{k}}(t),t)\delta u(t),\quad \delta u = u - {{u}^{k}},\quad \delta x = x - {{x}^{k}}.$

Далее линеаризуем ограничения (3.2), (3.3) на $k$-м приближении:

(3.16)

${{I}^{k}} + \sum\limits_{j = 0}^N \,{{\nabla }_{u}}{{I}^{k}}({{t}^{j}}){\kern 1pt} '({{u}_{j}} - u_{j}^{k}) + {{\nabla }_{w}}{{I}^{{k{\kern 1pt} '}}}(w - {{w}^{k}}) + {{\nabla }_{v}}{{I}^{{k{\kern 1pt} '}}}(v - {{v}^{k}}) = 0,$

(3.17)

$J_{j}^{k} + \sum\limits_{i = 0}^j \,[{{\nabla }_{u}}{{J}^{k}}({{t}^{j}}){\kern 1pt} '({{u}_{i}} - u_{i}^{k}) + {{\nabla }_{w}}{{J}^{k}}({{t}^{j}}){\kern 1pt} '(w - {{w}^{k}}) + {{\nabla }_{v}}{{J}^{k}}({{t}^{j}}){\kern 1pt} '(v - {{v}^{k}})] = 0,\quad j = \overline {0,N} .$

Здесь $I = ({{I}_{1}},{{I}_{2}}, \ldots ,{{I}_{m}})$, $J = ({{J}_{1}},{{J}_{2}}, \ldots ,{{J}_{s}})$. Следовательно, имеем $m$ ограничений (3.16) и $(N + 1)s$ ограничений (3.17), которые представляют собой явную форму (через $u$, $w$, $v$) линеаризованных $({{h}^{L}},g_{j}^{L})$ ограничений (3.2), (3.3), причем вместо равенств (3.3), заданных для каждого момента $t \in T$, имеем $N$ равенств, определенных в узлах сетки (3.8).

Линеаризуем также условия (3.4):

(3.18)

$c({{u}^{k}},{{t}^{j}}) + {{\nabla }_{u}}c({{u}^{k}},{{t}^{j}}){\kern 1pt} '({{u}_{j}} - u_{j}^{k}) = 0,\quad j = \overline {0,N} ,$

где $c = ({{c}_{1}},{{c}_{2}}, \ldots ,{{c}_{l}})$. Прямые ограничения на управление и параметры оставим без изменений:

(3.19)

$u_{j}^{{\text{н}}} \leqslant {{u}_{j}} \leqslant u_{j}^{{\text{в}}},\quad j = \overline {1,N} ,$

(3.20)

$v_{j}^{{\text{н}}} \leqslant {{v}_{j}} \leqslant v_{j}^{{\text{в}}},\quad j = \overline {1,n} ,\quad w_{i}^{{\text{н}}} \leqslant {{w}_{i}} \leqslant w_{i}^{{\text{в}}},\quad i = \overline {1,p} .$

Конечномерная аппроксимация функционала (3.15), в котором переменные $x(t)$ определены через систему (3.1) по заданному $u(t),\;t \in T$, имеет следующий вид:

(3.21)

$\begin{gathered} L = \varphi ({{x}^{N}}) - {{\lambda }^{{k{\kern 1pt} '}}}[h(x({{t}^{N}})) - {{{\bar {h}}}^{L}}] + \frac{\rho }{2}[h(x({{t}^{N}})) - {{{\bar {h}}}^{L}}]{\kern 1pt} '[h(x({{t}^{N}})) - {{{\bar {h}}}^{L}}] - \sum\limits_{j = 0}^N \,\mu _{j}^{{k{\kern 1pt} '}}[g(x({{t}^{j}}),{{t}^{j}}) - {{{\bar {g}}}^{L}}] + \\ + \;\frac{\rho }{2}\sum\limits_{j = 0}^N \,[g(x({{t}^{j}}),{{t}^{j}}) - {{{\bar {g}}}^{L}}]{\kern 1pt} '[g(x({{t}^{j}}),{{t}^{j}}) - {{{\bar {g}}}^{L}}] - \sum\limits_{j = 0}^N \,\gamma _{j}^{{k'}}[c({{u}_{j}},{{t}^{j}}) - {{{\bar {c}}}^{L}}] + \\ + \;\frac{\rho }{2}\sum\limits_{j = 0}^N \,[c({{u}_{j}},{{t}^{j}}) - {{{\bar {c}}}^{L}}]{\kern 1pt} '[c({{u}_{j}},{{t}^{j}}) - {{{\bar {c}}}^{L}}]{\kern 1pt} {\kern 1pt} , \\ \end{gathered} $

где

${{\bar {h}}^{L}} = h({{x}^{k}}({{t}^{N}})) + {{h}_{x}}({{x}^{k}}({{t}^{N}}))(x({{t}^{N}}) - {{x}^{k}}({{t}^{N}}))$,

${{\bar {g}}^{L}} = g({{x}^{k}}({{t}^{j}}),{{t}^{j}}) + {{g}_{x}}({{x}^{k}}({{t}^{j}}),{{t}^{j}})(x({{t}^{j}}) - {{x}^{k}}({{t}^{j}})),$

${{\bar {c}}^{L}} = c({{u}^{k}}({{t}^{j}}),{{t}^{j}}) + {{c}_{u}}({{u}^{k}}({{t}^{j}}),{{t}^{j}})({{u}_{j}} - u_{j}^{k}),\quad j = \overline {0,N} .$

Для минимизации функционала (3.21), который теперь по сути является функцией многих переменных, при линейных ограничениях (3.16)–(3.20) применяется метод приведенного градиента (см. [6]). Заметим, что функционал (3.21) предполагает использование исходной системы (3.1) для расчета траектории $\{ x({{t}^{1}}),x({{t}^{2}}), \ldots ,x({{t}^{N}})\} $ по заданным параметрам $v$, $w$ и управлению $u({{t}^{0}}),u({{t}^{1}}),$ $ \ldots ,u({{t}^{N}})$, т.е. полная модель вспомогательной задачи описывается соотношениями (3.1), (3.16)–(3.21).

Обозначив через $A[m + (l + s)(N + 1)] \times [r(N + 1) + p + n]$ матрицу коэффициентов линейных равенств (3.16)–(3.18), через $b$ – вектор их свободных членов размерности $m + (l + s)(N + 1)$ и через $z$ – вектор искомых переменных (${{u}_{j}},j = \overline {0,N} $; $v$; $w$) размерности $r(N + 1) + p + n$ соответственно, поставленную задачу запишем в виде

(3.22)

$\begin{gathered} L(z) \to min, \\ Az = b, \\ {{z}^{{\text{н}}}} \leqslant z \leqslant {{z}^{{\text{в}}}}. \\ \end{gathered} $

Для решения задачи (3.22) применяется метод приведенного градиента (см. [6]), который отличается от известного в линейном программировании симплекс-метода тем, что в силу нелинейности целевой функции его последовательные приближения необязательно будут находиться в вершинах многогранника линейных ограничений, а могут быть и его внутренними точками.

3.4. Алгоритм метода спроектированного лагранжиана (см. [6]–[9])

Рассмотрим теперь полный алгоритм решения исходной задачи (3.1)–(3.6).

1. С заданным управлением $u_{j}^{k},\;j = \overline {0,N} $, интегрируется система (3.1), и в узлах сетки (3.8) запоминаются точки фазовой траектории $x_{j}^{k},\;j = \overline {0,N} $. Здесь $k$ – номер итерации (первый раз $k = 0$).

На полученном решении линеаризуются ограничения задачи (3.10) и строится вспомогательная задача (3.16)–(3.21).

2. Методом приведенного градиента решается вспомогательная задача минимизации модифицированной функции Лагранжа (3.21) при линейных ограничениях (3.16)–(3.20).

В результате будут найдены новые приближения для управления $u_{j}^{{k + 1}},\;j = \overline {0,N} $, параметров ${{w}^{{k + 1}}}$ и ${{v}^{{k + 1}}}$, а также для двойственных переменных ${{\lambda }^{{k + 1}}}$ и $\mu _{j}^{{k + 1}},\;j = \overline {0,N} $.

3. Проверяется критерий окончания итерационного процесса как по прямым, так и по двойственным переменным:

$\left| {{{I}_{i}}({{u}^{{k + 1}}},{{w}^{{k + 1}}},{{v}^{{k + 1}}})} \right|{\text{/}}(1 + {{\alpha }^{{k + 1}}}) \leqslant \varepsilon ,\quad i = \overline {1,m} ,$

$\left| {{{J}_{i}}({{u}^{{k + 1}}},{{w}^{{k + 1}}},{{v}^{{k + 1}}})} \right|{\text{/}}(1 + {{\alpha }^{{k + 1}}}) \leqslant \varepsilon ,\quad i = \overline {1,s} ,$

где

${{\alpha }^{{k + 1}}} = max\{ {\kern 1pt} {\text{|}}{\kern 1pt} {\text{|}}u_{j}^{{k + 1}}{\text{|}}{\kern 1pt} {\text{|}},j = \overline {0,N} ;\left| {{{w}_{i}}} \right|,i = \overline {1,p} ;\left| {{{v}_{l}}} \right|,l = \overline {1,n} \} $;

$\left| {\lambda _{j}^{k} - \lambda _{j}^{{k + 1}}} \right|{\text{/}}(1 + {{\Theta }^{{k + 1}}}) \leqslant \varepsilon ,\quad j = \overline {1,m} ;$

$\left| {\mu _{{ij}}^{k} - \mu _{{ij}}^{{k + 1}}} \right|{\text{/}}(1 + {{\Theta }^{{k + 1}}}) \leqslant \varepsilon ,\quad i = \overline {1,s} ,\quad j = \overline {0,N} ;$

${{\Theta }^{{k + 1}}} = max\{ {\text{|}}\lambda _{j}^{{k + 1}}{\text{|}},j = \overline {1,m} ;{\text{|}}\mu _{{ij}}^{{k + 1}}{\text{|}},i = \overline {1,s} ,j = \overline {0,N} \} .$

При нарушении хотя бы одного из этих условий выполняется новая $(k + 1)$-я итерация с п. 1. Если же эти неравенства выполняются для заданного $\varepsilon > 0$, то итерационный процесс прекращается, а найденные $u_{j}^{{k + 1}},\;j = \overline {0,N} ,\;{{w}^{{k + 1}}}$ и ${{v}^{{k + 1}}}$ выдаются в качестве приближенного решения задачи оптимального управления.

4. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ. ПРИКЛАДНЫЕ ЗАДАЧИ С РЕШЕНИЯМИ

Приведем примеры прикладных задач, относящихся к рассмотренным в разд. 1–3 классам задач оптимизации, решения которых найдены с помощью многометодных алгоритмов, изложенных в этих разделах.

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

Задача сформулирована М.М. Свининым и приведена в [10]. Рассматривается мобильный сферический робот, перемещающийся по плоскости. Конструкция робота представляет собой оболочку с размещенными в ней тремя роторами-двигателями. Динамика робота, редуцированная к контактным координатам, описывается следующей системой обыкновенных дифференциальных уравнений:

$\dot {x} = G(x){{J}^{{ - 1}}}(x){{J}_{r}}\sum\limits_{k = 1}^n \,{{n}_{k}}(x){{u}_{k}},$

где векторы состояний и управлений определены как

$x \triangleq {{[{{u}_{a}},{{v}_{a}},{{u}_{o}},{{v}_{o}},\psi ]}^{{\text{T}}}},\quad u \triangleq {{[{{\dot {\varphi }}_{1}},{{\dot {\varphi }}_{2}},{{\dot {\varphi }}_{3}}]}^{{\text{T}}}},$

а ${{\varphi }_{i}},\;i = \overline {1,3} $, обозначают углы поворота двигателей.

Положение точки контакта на плоскости задается координатами ${{u}_{a}},\;{{v}_{a}}$, а ее координаты на сфере задаются углами ${{u}_{o}},\;{{v}_{o}}$.

Матричные и векторные величины определяются следующим образом:

$G = \left[ {\begin{array}{*{20}{c}} 0&{ - R}&0 \\ R&0&0 \\ {sin\psi {\text{/}}\cos {{v}_{o}}}&{\cos \psi {\text{/}}\cos {{v}_{o}}}&0 \\ {cos\psi }&{ - sin\psi }&0 \\ {sin\psi \operatorname{tg} {{v}_{o}}}&{cos\psi \operatorname{tg} {{v}_{o}}}&1 \end{array}} \right],$

а векторы ${{n}_{1}},\;{{n}_{2}},\;{{n}_{3}}$ – столбцы матрицы

$R = \left[ {\begin{array}{*{20}{c}} {cos{{u}_{o}}cos\psi + sin{{u}_{o}}sin{{v}_{o}}sin\psi }&{cos{{{v}}_{o}}sin\psi }&{ - sin{{u}_{o}}cos\psi + cos{{u}_{o}}sin{{v}_{o}}sin\psi } \\ { - cos{{u}_{o}}sin\psi + sin{{u}_{o}}sin{{v}_{o}}cos\psi }&{cos{{v}_{o}}cos\psi }&{sin{{u}_{o}}sin\psi + cos{{u}_{o}}sin{{v}_{o}}cos\psi } \\ {sin{{u}_{o}}cos{{v}_{o}}}&{ - sin{{v}_{o}}}&{cos{{u}_{o}}cos{{v}_{o}}} \end{array}} \right].$

Матрица инерции системы определяется как

$J = \left[ {\begin{array}{*{20}{c}} {(2{\text{/}}3{{m}_{o}} + M){{R}^{2}}}&0&0 \\ 0&{(2{\text{/}}3{{m}_{o}} + M){{R}^{2}}}&0 \\ 0&0&{2{\text{/}}3{{m}_{o}}{{R}^{2}}} \end{array}} \right] + (2{{J}_{p}} + {{J}_{r}})E,$

где $M$ – общая масса робота, а ${{m}_{o}}$ и ${{m}_{r}}$ обозначают, соответственно, массу сферической оболочки и массу отдельного ротора.

Задача оптимального управления сферическим роботом состоит в переводе из точки $x(0)$ в точку $x(T)$ при условии минимизации энергии управления $J = \int_0^T {{{u}^{T}}udt} $.

Пусть, например, $x(0) = [0,0,0,0,0]$, а $x(10) = \left[ {0,2,0,3,0,0,\frac{\pi }{6}} \right]$, тогда полученные решения представлены на фиг. 3 и являются физически реализуемыми. Точность выполнения ограничений порядка $\mathop {10}\nolimits^{ - 5} ,$ значение функционала $J(10) = 5769001$.

Фиг. 3.

Найденные управления и траектории для мобильного робота.

4.2. Задача об оптимальном управлении манипулятором робота

Динамика движения манипулятора промышленного робота описывается системой дифференциальных уравнений

${{\dot {x}}_{1}} = {{x}_{2}},$

${{\dot {x}}_{2}} = \frac{{[{{M}_{1}}(x,u) - {{F}_{1}}(x)]{{a}_{{22}}} - [{{M}_{2}}(x,u) - {{F}_{2}}(x)]{{a}_{{12}}}(x)}}{{{{a}_{{11}}}{{a}_{{22}}} - {{a}_{{12}}}(x){{a}_{{21}}}(x)}},$

${{\dot {x}}_{3}} = {{x}_{4}},$

${{\dot {x}}_{4}} = \frac{{[{{M}_{2}}(x,u) - {{F}_{2}}(x)]{{a}_{{11}}} - [{{M}_{1}}(x,u) - {{F}_{1}}(x)]{{a}_{{21}}}(x)}}{{{{a}_{{11}}}{{a}_{{22}}} - {{a}_{{12}}}(x){{a}_{{21}}}(x)}},$

где

${{M}_{1}}(x,u) = - {{c}_{1}}({{x}_{1}} - {{u}_{1}}),\quad {{M}_{2}}(x,u) = - {{c}_{2}}({{x}_{3}} - {{x}_{1}} - {{u}_{2}}),$

${{F}_{1}}(x) = - {{m}_{2}}{{l}_{1}}{{R}_{2}}sin({{x}_{3}} - {{x}_{1}})x_{2}^{2},\quad {{F}_{2}}(x) = {{m}_{2}}{{l}_{2}}{{R}_{2}}sin({{x}_{3}} - {{x}_{1}})x_{4}^{2},$

${{a}_{{11}}} = {{m}_{1}}\rho _{1}^{2} + {{m}_{2}}l_{1}^{2},\quad {{a}_{{12}}} = {{a}_{{21}}} = {{m}_{2}}{{R}_{1}}{{l}_{1}}cos({{x}_{3}} - {{x}_{1}}),\quad {{a}_{{22}}} = {{m}_{2}}\rho _{2}^{2}.$

Для рассматриваемой модели робота ${{m}_{1}} = 7.62$, ${{m}_{2}} = 8.73$, ${{R}_{1}} = 0.239$, ${{R}_{2}} = 0.251$, ${{\rho }_{1}} = 0.968$, ${{\rho }_{2}} = 0.973$, ${{l}_{1}} = 0.5$, ${{l}_{2}} = 0.67$, ${{c}_{1}} = {{c}_{2}} = 10$. На траекторию движения накладываются ограничения ${\text{|}}{{M}_{i}}(x,u){\kern 1pt} {\text{|}} \leqslant 10,$ $i = 1,2$, $\pi {\text{/}}6 \leqslant {{x}_{1}}(t) \leqslant 5{\text{/}}6\pi $, $\pi {\text{/}}3 \leqslant {{x}_{1}}(t) - {{x}_{3}}(t) \leqslant 5{\text{/}}6\pi ,$ $t \in [0,T]$.

Необходимо найти управление, переводящее систему из точки $x{\kern 1pt} '(0) = (\pi {\text{/}}6,0, - \pi {\text{/}}6,0)$ в точку $x{\kern 1pt} '(T) = (5{\text{/}}6\pi ,0,\pi {\text{/}}3,0)$ за минимальное время $T$.

Начиная с приближения $T = 4$, ${{u}_{1}}(t) = 0$, ${{u}_{2}}(t) = 0,$ $t \in [0,T]$, было получено решение, на котором невязки по ограничениям не превысили ${{10}^{{ - 3}}}$, а значение функционала составило 2.88. Вид оптимального управления и соответствующих ему фазовых координат приведен на фиг. 4.

Фиг. 4.

Графики управлений и фазовых координат для задачи 4.2.

4.3. Задача оптимизации электроэнергетической системы (ЭЭС)

Математическая модель ЭЭС разработана в Институте систем энергетики им. Л.А. Мелентьева СО РАН и представляет собой совокупность подсистем, описывающих генерирование и потребление электроэнергии, которые объединены в единую систему уравнениями электрической сети, и много лет успешно применяется для расчета различных режимов работы проектируемых ЭЭС. Рассмотрим небольшую модель ЭЭС, состоящую из $m$ синхронных генераторов и $m$ паровых турбин. Каждая синхронная машина описывается дифференциальными уравнениями Парка–Горева (без учета переходных процессов в обмотке статора), которые после приведения их к стандартной форме принимают следующий вид:

$\mathop {\dot {x}}\nolimits_i = {{x}_{{m + i}}},\quad i = \overline {1,m} ,$

$\begin{gathered} \mathop {\dot {x}}\nolimits_{m + i} = {{\omega }_{0}}{\text{/}}({{T}_{{ji}}}{{P}_{{{{n}_{i}}}}})[{{x}_{{5m + i}}} - U_{i}^{1}{\text{/}}X_{{{{d}_{i}}}}^{1}({{x}_{{2m + i}}}sin{{x}_{i}} + {{x}_{{3m + i}}}cos{{x}_{i}}) + \\ \, + U_{i}^{2}{\text{/}}X_{{{{d}_{i}}}}^{1}({{x}_{{2m + i}}}cos{{x}_{i}} - {{x}_{{3m + i}}}sin{{x}_{i}}) - {{D}_{i}}{{x}_{{m + i}}}],\quad i = \overline {1,m} , \\ \end{gathered} $

${{\dot {x}}_{{2m + i}}} = 1{\text{/}}{{T}_{{{{d}_{0}}i}}}[{{x}_{{4m + i}}} - {{x}_{{2m + i}}} + ({{X}_{{{{d}_{i}}}}} - X_{{{{d}_{i}}}}^{1}){\text{/}}X_{{{{d}_{i}}}}^{1}(U_{i}^{2}sin{{x}_{i}} - {{x}_{{2m + i}}} + U_{i}^{1}cos{{x}_{i}})],\quad i = \overline {1,m} ,$

$\begin{gathered} {{{\dot {x}}}_{{4m + i}}} = 1{\text{/}}T{{R}_{i}}[{{E}_{{{{R}_{0}}i}}} - {{x}_{{4m + i}}} + {{K}_{{Ui}}}({{U}_{{0i}}} - \sqrt {{{{(U_{i}^{1})}}^{2}} + {{{(U_{i}^{2})}}^{2}}} + {{K}_{{{{I}_{i}}}}}{\text{/}}X_{{{{d}_{i}}}}^{1}( - {{I}_{{0i}}}X_{{{{d}_{i}}}}^{1} + ({{({{x}_{{2m + i}}})}^{2}} + \\ \, + {{({{x}_{{3m + i}}})}^{2}} + {{(U_{i}^{1})}^{2}} + {{(U_{i}^{2})}^{2}} - 2cos{{x}_{i}}(U_{i}^{2}{{x}_{{3m + i}}} + U_{i}^{1}{{x}_{{2m + i}}}) + 2sin{{x}_{i}}(U_{i}^{1}{{x}_{{3m + i}}} - U_{i}^{2}{{x}_{{2m + i}}}){{)}^{{1/2}}}) + \\ \, + {{K}_{{{{f}_{i}}}}}{\text{/}}2\pi {{x}_{{m + i}}} + K_{{{{f}_{i}}}}^{1}{\text{/}}2\pi {{x}_{{2m + i}}}],\quad i = \overline {1,m} . \\ \end{gathered} $

Уравнения динамики паровых турбин следующие:

$\mathop {\dot {x}}\nolimits_{5m + i} = 1{\text{/}}{{T}_{{si}}}[ - {{P}_{{ni}}}{\text{/}}{{\omega }_{0}}{{\sigma }_{{1i}}}{{x}_{{m + i}}} + {{x}_{{6m + i}}} - {{x}_{{5m + i}}}],\quad i = \overline {1,m} ,$

$\mathop {\dot {x}}\nolimits_{6m + i} = 1{\text{/}}{{T}_{{p2i}}}[ - {{P}_{{ni}}}{\text{/}}{{\omega }_{0}}{{\sigma }_{{2i}}}{{x}_{{m + i}}} + {{u}_{i}}],\quad i = \overline {1,m} .$

Переменными состояния ${{x}_{{jm + i}}},\;j = \overline {0,6} $, являются угол ротора генератора, скольжение, составляющие переходной ЭДС машины в продольной и поперечной осях и напряжение обмотки возбуждения соответственно для каждого $i = \overline {1,m} $. Управление ${{U}_{i}}$ изменяет установку регулятора скорости, чтобы обеспечить устойчивый динамический переход в заданное послеаварийное состояние при аварийных сбросах нагрузки. В правые части уравнений входят также технические параметры генераторов и турбин, смысл которых здесь приводить не будем. Число генераторов и турбин было задано равным пяти. Таким образом, при $m = 5$ число дифференциальных уравнений равно 35, т.е. $n = 35$.

Модель электрической сети составляют алгебраические уравнения на узловые напряжения, в которые входят и переменные состояния. Эти уравнения обычно задаются в комплексных переменных, а при численном решении выполняется переход к действительным переменным. Так, например, для эксперимента было задано $N = 14$ уравнений в комплексных переменных, а после перехода была получена система из 28 алгебраических уравнений

$C_{i}^{1}U_{i}^{1} - C_{i}^{2}U_{i}^{2} + \sum\limits_{k = 1,k \ne i}^N \,( - U_{k}^{1}Y_{{ik}}^{1} + U_{k}^{2}Y_{{ik}}^{2}) = 1{\text{/}}X_{{{{d}_{i}}}}^{1}({{x}_{{2m + i}}}sin{{x}_{i}} + {{x}_{{3m + i}}}cos{{x}_{i}}),\quad i = \overline {1,N} ,$

$C_{i}^{1}U_{i}^{2} + C_{i}^{2}U_{i}^{1} + \sum\limits_{k = 1,k \ne i}^N \,( - U_{k}^{1}Y_{{ik}}^{2} - U_{k}^{2}Y_{{ik}}^{1}) = 1{\text{/}}X_{{{{d}_{i}}}}^{1}({{x}_{{2m + i}}}cos{{x}_{i}} + {{x}_{{3m + i}}}sin{{x}_{4}}),\quad i = \overline {1,N} ,$

где $C_{i}^{1} = {{P}_{{ni}}}{\text{/}}{{({{U}_{i}})}^{2}} + A_{i}^{1},$ $C_{i}^{2} = A_{i}^{2} - {{Q}_{{ni}}}{\text{/}}{{({{U}_{i}})}^{2}},$ ${{({{U}_{i}})}^{2}} = {{(U_{i}^{1})}^{2}} + {{(U_{i}^{2})}^{2}}$.

Кроме того, заданы ограничения на фазовые координаты и управление

${\text{|}}{\kern 1pt} {{x}_{i}}(t) - {{x}_{j}}(t){\kern 1pt} {\text{|}} \leqslant {{\delta }_{{max}}},\quad i,j = \overline {1,m} ,$

${{x}_{{mi{{n}_{i}}}}} \leqslant {{x}_{{4m + i}}}(t) \leqslant {{x}_{{ma{{x}_{i}}}}},\quad i = \overline {1,2m} ,$

${{U}_{{mi{{n}_{i}}}}} \leqslant {{U}_{i}}(t) \leqslant {{U}_{{ma{{x}_{i}}}}},\quad i = \overline {1,m} .$

Целевым функционалом является функция конечного состояния системы (при $t = 10$ с), измеряющая отклонение некоторых фазовых координат от заданных величин (например, мощностей).

При численном решении поставленной задачи методом спроектированного лагранжиана было сделано 11 внешних итераций, каждая из которых содержала около 20 внутренних итераций метода приведенного градиента. Заданные равенства были выполнены с точностью ${{10}^{{ - 6}}}$, при этом ни одна из переменных не вышла на заданные ограничения. Полученное оптимальное управление обеспечивает вывод ЭЭС в требуемый режим работы за 10 с после резкого сброса нагрузки.

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

ЗАКЛЮЧЕНИЕ

Практика показывает, что последовательность приближений многометодного алгоритма в сложных задачах управления, как правило, состоит из приближений нескольких (3–5) численных методов, которые выбирались по заданному критерию автоматически в процессе оптимизации. Проведенные вычислительные эксперименты подтверждают эффективность предложенной технологии при решении прикладных задач оптимального управления. Установлено, что применение многометодной технологии нередко является единственным способом получения численного решения в сложной задаче оптимального управления, так как сходимость каждого из методов в отдельности прекращалась до получения оптимального решения. Современные информационные технологии и многопроцессорная вычислительная техника допускают достаточно эффективную реализацию многометодных алгоритмов. Программное обеспечение, разработанное на основе данного подхода и реализующее многометодную технологию расчета оптимального управления и оптимальных параметров (см. [9]–[12]), успешно применяется для решения сложных прикладных задач оптимального управления из различных областей науки и техники (см. [10]–[14]). Применение эффективной технологии расчета управления особенно актуально в управляемых системах реального времени, например, в системах управления летательными аппаратами, обладающими высокой маневренностью. Например, при проектировании СУ-57 (мирового лидера по маневренности) для решения серии задач оптимального маневрирования использовалось данное программное обеспечение (см. [11]). Известно также, что наличие программы выбора оптимального начального приближения в бортовом компьютере беспилотного космического аппарата Буран позволило ему выбрать наименее зависимую от бокового ветра стартовую точку для начала посадки на аэродром и успешно завершить свой беспрецедентный полет, реализовав программу оптимальной посадки (глиссаду) с высокой точностью.

Список литературы

Понтрягин Л.С., Болтянский В.Н., Гамкрелидзе Р.В., Мищенко Е.Ф. Математическая теория оптимальных процессов. М.: Наука, 1969.
Беллман Р., Калаба Р. Квазилинеаризация и нелинейные краевые задачи. М.: Мир, 1968.
Габасов Р., Кириллова Ф.М. Качественная теория оптимальных процессов. М.: Наука, 1971.
Васильев О.В., Тятюшкин А.И. Об одном методе решения задач оптимального управления, основанном на принципе максимума // Ж. вычисл. матем. и матем. физ. 1981. Т. 21. № 6. С. 1376–1384.
Габасов Р., Кириллова Ф.М., Тятюшкин А.И. Конструктивные методы оптимизации. Ч. 1: Линейные задачи. Минск: Университетское, 1984.
Гилл Ф., Мюррей У., Райт М. Практическая оптимизация. М.: Мир, 1985.
Евтушенко Ю.Г. Методы решения экстремальных задач и их применение в системах оптимизации. М.: Наука, 1982.
Федоренко Р.П. Приближенное решение задач оптимального управления. М.: Наука, 1978.
Тятюшкин А.И. ППП КОНУС для оптимизации непрерывных управляемых систем // Пакеты прикладных программ: Опыт использования. М.: Наука, 1989. С. 63–83.
Горнов А.Ю., Тятюшкин А.И., Финкельштейн Е.А. Численные методы для решения терминальных задач оптимального управления // Ж. вычисл. матем. и матем. физ. 2016. Т. 56. № 2. С. 224–237.
Тятюшкин А.И., Федунов Б.Е. Возможности защиты от атакующей ракеты задней полусферы самолета вертикальным маневром // Изв. РАН, ТиСУ. 2006. № 1. С. 111–125.
Тятюшкин А.И. Многометодная технология оптимизации управляемых систем. Новосибирск: Наука, 2006.
Тятюшкин А.И. Численные методы решения задач оптимального управления с параметрами // Ж. вычисл. матем. и матем. физ. 2017. Т. 57. № 10. С. 1615–1630.
Тятюшкин А.И. Многометодная оптимизация управления в сложных прикладных задачах // Ж. вычисл. матем. и матем. физ. 2019. Т. 59. № 2. С. 235–246.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал