Метод аппроксимации экстремальных решений в задачах оптимального управления с фазовыми ограничениями

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

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

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

Альтернативным подходом является использование базисных функций с локальным носителем, таких как сплайны и B-сплайны [2]. Сплайновая параметризация гарантирует заданную гладкость траекторий, а также порождает разреженные матрицы якобианов, что снижает вычислительную нагрузку. Тем не менее реализация строгих фазовых ограничений на всем непрерывном интервале при использовании сплайнов математически затруднительна и часто требует введения сложных полубесконечных условий либо искусственной релаксации ограничений.

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

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

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

• 1.    Методы штрафных функций, регуляризации и аугментированного лагранжиана. Для сглаживания нелинейностей в ограничениях состояния широкое развитие получили усовершенствованные барьерные методы и алгоритмы регуляризации Моро — Иосиды. При решении сильно дискретизированных задач успешно применяются методы аугментированного лагранжиана, а также проекционные методы на основе алгоритмов расщепления (например, Дугласа — Рекфорда). Данные подходы позволяют эффективно разделять сложную оптимизационную задачу на последовательность более простых подзадач, гарантируя сходимость для широкого класса выпуклых и невыпуклых ограничений [4].

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

• 3.    Машинное обучение и безопасное обучение с подкреплением. Отдельным направлением стало внедрение нейросетевых архитектур . Глубокие нейронные сети и алгоритмы обучения с подкреплением применяются для аппроксимации сопряженных переменных и синтеза субоптимальных законов управления. Формирование барьерных функций управления совместно с обучением с подкреплением позволяет гарантировать соблюдение фазовых ограничений в режиме реального времени [6]. В ряде нелинейных задач машинное обучение сегодня выступает серьезной альтернативой традиционному оптимальному управлению [7].

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

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

• 1    Задача и метод

Рассматривается класс задач оптимального управления с фазовыми ограничениями, приводимых к следующему каноническому виду:

• x(t) = f (x(t),u(t),t),x(10) = x°,u(t)g U c Rm, t gT = [10;t} ],(1)

J0 (u ) = ^0 (x(ti )) + JFo (x (t),u (t),t)dt ^ inf,(2)

T

J1 (u ) = Pi (x (t 1 )) =0, в котором x(t)g Rn — кусочно-дифференцируемая векторная функция, u (t) — m -мерный вектор управляющих функций; U — замкнутое выпуклое множество. Интервал T фиксирован. В качестве доступных управляющих функций рассматривается множество V кусочнонепрерывных на T функций со значениями в множестве U : V = {u g PC (T): u (t )g U, t g T}. Функции p0 (x), p1 (x) непрерывно дифференцируемы на Rn, функции F0 (x,u, t), f (x,u, t) и их частные производные по x, u непрерывны по совокупности аргументов на множестве Rn х U х T. Функция f (x, u, t) удовлетворяет условию Липшица по x в Rn х U х T с константой L > 0 :

||^f ( ^x , ^u , t ) ^{- f} ( y , ^u , t )|| < L||x ^- y ||.

К виду (1)-(3) стандартными способами штрафования за нарушение ограничений могут быть сведены многие задачи оптимального управления с фазовыми и терминальными ограничениями.

Приближенное решение исходной задачи (1)-(3) будем искать в классе кусочно-постоянных управлений. Для этого введем разбиение фиксированного временного интервала T = [ 1 ₀; t ₁ ] на N подынтервалов

Tk — [tk-1;tk],k — 1,2,—,N, таких что tN — t1. Аппроксимирующая задача рассматривается в следующем виде:

X ( t ) = f ( x ( t ) , u_k , t ) , x ( 1 0 ) = x °, u_k g U c R m , t ^G T k = [ t k - 1 ; t k ], ^{k —} ^A-K N ,

N tk inf , u={Ui,—,Un }gQ

^ф о ( ^u ) = P 0 ( ^x ( t N ) ) + E f F 0 ( ^x ( t ) , u k , t ) ^dt ^ ^{k = 1} t_t - i

ф1 ( u ) = Pi (x ( tN )) = 0.

Допустимое кусочно-постоянное управление определяется N -мерным набором m -мерных векторов u — { u ₁, — , u_N } , u k g U , k — 1,2, — , N . Обозначим О множество допустимых наборов векторов управлений.

Обозначим через x ( t , v ) , t g T = [ 1 ₀; t_N ] решение системы (4) при управлении v — { v ₁, — , v_N } gO . Значения x ( t , v ) , t g T_k , k — 1,2, — , N определяются путем последовательного интегрирования системы (4) на интервалах T_k при u_k — v_k , t g T_k , k — 1,2, — , N .

Задача (4)-(6) может рассматриваться как конечномерная задача математического программирования.

Рассмотрим вспомогательную задачу без ограничений на основе функционала Лагранжа:

■x ( t ) = f ( x ( t ) , U k , t ) , x ( t о ) = x °, U k G U c R^m , t ^{G T} k = [ t k - 1 ; t k ] , k = ^’K ^N ,

L ( u ) — Я ₀ Ф ₀ ( u ) + * Ф ₁ ( u )

^    inf   , Я = ( * 0 , * ) g R ², * ^ 0.

u —{ u 1 , — , un } eQ

Введем функцию Понтрягина с сопряженной переменной у g R n и стандартную сопряженную систему в задаче Лагранжа (7)-(8):

HЯ (у, x, w, t) — (у, f (x, w, t)) - Я, F0 (x, w, t), w g U c Rm, у (t) — -H (у (t), x(t), uk, t), t G Tk, k —1,2,—, N, У (tN)— -P* (x(tN)), (РЯ (x)— AM (x)+ Я191 (x).

Для допустимого управления v — { v ₁, — , v_N } gQ обозначим через у ^Я ( t , v ) , t g T — [ 1 ₀; t_N ] решение стандартной сопряженной системы (9), полученное последовательным интегрированием на интервалах разбиения T_k , k — 1,2, — , N при u k — v k , t g T_k , k — 1,2, — , N и x ( t ) — x ( t , v ), t g T .

В соответствии с известной формулой приращения целевого функционала [8] для задач оптимального управления без ограничений в классе кусочно-непрерывных управлений имеет место формула приращения функции Лагранжа в задаче (7), (8):

А

L ( ^v ) ^- L ( ^u ) ^—-EL\ H u ( у ( t , ^u ) , ^x ( t , ^u ) , ^uk , t ) , ^{A u} k] dt + o |EI^{a u} k |||^,(10) k — 1 ‘^k у k — 1 )

где H обозначает производную функции Понтрягина по переменной управления, A u k = v _k - u _k .

Обозначим через P Y оператор проектирования на множество Y с R m в евклидовой норме:

^Py ( z ) = ^агё ^min (| 1- У — z ||) , z e R m • y e Y

На основе формулы приращения (10) можно получить необходимое условие оптимальности в задаче (4)-(6) для управления u еО в следующем виде при некотором А ^ 0 c параметром a >  0:

^u k = P U ( ^u k + ^a L H u ( ^ ^A ( t , ^u ) , x ( t , ^u ) , ^u k , t ) ^dt ) , ^k = ^1,2, . ^, N ,

\ ^J T                                    )                        (11)

Вырожденный случай (А0 = 0, А ^ 0) необходимого условия оптималь- ности в конкретных задачах оптимального управления, как правило, исследуется аналитически. В невырожденном случае (А0 = 1, А е R) полученная система уравнений (11) относительно пары неизвестных (u,А) е О хR решается численными методами.

В данной работе для невырожденного случая рассматривается подход к поиску экстремальных управлений, основывающийся на представлении необходимого условия оптимальности (11) в форме специальной задачи о неподвижной точке оператора управления в конечномерном пространстве допустимых наборов векторов u = {u1,.,uN}. Такое представление дает возможность применить и модифицировать известную теорию и методы неподвижных точек для конструирования итерационного алгоритма поиска экстремального управления как решения рассматриваемой задачи о неподвижной точке.

При заданном a >  0 для решения системы (11) рассматривается итерационный процесс с индексом s >  0 с заданным начальным приближением управления u ⁰ еО при s = 0 :

u k ^{+ 1} = P u ( u k + a f H A [ u ( t , u^s ) , x ( t , u^s ) , u k , t ) dt ) , k = 1,2, . , N ,

Tk

^ 1 ( x ( t N , u^{s + 1} ) ) = 0.

На каждой итерации предлагаемого процесса решается неявно заданное уравнение   ^ (x(tN,us+1

относительно скалярного

множителя Лагранжа A e R • При этом начальное приближение u ⁰ еО при s = 0 можно выбирать недопустимым.

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

• 2    Модель биореактора

Рассмотрим известную непрерывную модельную задачу с тремя двухсторонними фазовыми ограничениями [10–15] (оптимизация подпитываемого биореактора для производства пенициллина путем анаэробной ферментации глюкозы):

u (t) X, (t)

X j ( t ) = hx (t)-- ——^, x (0) = 1,5,

1    ( ) 11 ( ) 500x4 (t) 1( )

^X 2 ( t ) = ^h 2 ^X 1 ( t ) ^{- 0,01 X} 2 ( t ) ^- uTt^rTX , x 2 ( 0 ) = 0, 500 X 4 ( t )

X3 (t) =

h ₁ x ₁ ( t ) h ₂ x ₁ ( t ) 0,029 x ₁ ( t ) x ₃ ( t ) u ( t ) -------1--

0,47

1,2     0,0001 + x 3 ( t )

^x 4 ( t )

1 - x -^) 1 , x (0) = 0, 500 J ’ 3           ’

/ x u ( t )       / X

X4(t ) = 500’X4(0) = 7, hi (Xi (t), xз (t )) =

0,11 x 3 ( t )         , , xx            0,0055 x 3 ( t )

, h ( X, (t)) =

0,006 x ₁ ( t ) + x 3 ( t ) ^Mi'" 0,0001 + x 3 ( t ) ( 1 + 10 x 3 ( t ) )

0 <  x ₁ ( t ) < 40,0 <  x 3 ( t ) < 25,0 <  x 4 ( t ) < 10, t e T = [ 0;132 ] ,

J ( u ) = - x ₂ ( 132 ) x 4 (132) ^ inf, ^{v 7}                                   ueV

V = {u e PC(T): u (t) e U, t e T}, U = {w e R : 0 < w < 50}.

В работе [15] задача решалась c использованием параллельного алгоритма, основанного на декомпозиции области управления, в классе кусочно-постоянных управлений с шестью нефиксированными моментами переключений.

Начальное приближение управления в [15] выбиралось равным u0 =(12,5; 12,5; 37,5; 12,5; 12,5; 12,5; 12,5 )T с моментами переключения управления т0 = (5; 10; 15; 20; 25; 50)T.

В [15] получено расчетное управление:

u * =(14,54; 11,24; 37,05; 13,89; 13,77; 12,75; 9,41)T с моментами переключения т* =(4,07;11,12; 15,33; 19,04; 24,84; 54,9)T.

При начальном приближении u ⁰ автором данной статьи было рассчитано значение функционала J ( u ⁰) = - 100,23 с максимальным нарушением фазовых ограничений M⁰ = 2,25. При расчетном управлении u * было получено значение J ( u *) = - 88,85 с максимальным нарушением M* = 0,061.

Полученные в [15] результаты иллюстрируются на рисунке 1.

— полученное управление u*(t), - - ■ начальное управление м° (?) •

— траекториях* (г), --- траекториях*(г), ........ траектория х|(у), — траекториях^?).

Рис. 1. Графики управления и фазовых траекторий, полученные в работе [15]

Проведем поиск экстремального управления рассматриваемой задачи в классе кусочно-постоянных управлений предлагаемым методом неподвижных точек.

Введением дополнительной фазовой переменной

Х5 (t) = Q1 (Х1 (t)) + Q3 (X3 (t)) + Q4 (X4 (t)) , X5 (0) = 0, t е T с использованием кубических штрафных функций

Q 1 ( X 1 ) = max(0, X 1 - 40)³ + max(0, - x 1 )³,

Q ₃ ( x ₃) = max(0, x ₃ - 25)³ + max(0, - x ₃ )³,

Q ₄( X ₄) = max(0, x ₄ - 10)³ + max(0, - x ₄)³

задача сводится к эквивалентной задаче с одним терминальным ограни-чени ем-рав енств ом:

J 0 ( и ) = - x ₂ ( 132 ) X ₄ ( 132 ) ^ inf, J 1 ( u ) = х ₅ ( 132 ) = 0.

Рассмотрим дискретно-непрерывную задачу на основе регулярного функционала Лагранжа

L (и) = -X2 (132)X4 (132) + 21 X5 (132) ^ inf, 2, е R в классе кусочно-постоянных управлений с разбиением интервала T = [0;132] на N подынтервалов Tk =[tk-1;tk],к = 1,2,...,N . Управление и = {u1, и 2,., uN} на каждом подынтервале принимает постоянное значение ик е U с R, к = 1,2,.,N :

ф 0 ( и ) = - x 2 ( 132 ) x ₄ ( 132 ) ^ inf,

Ф1 (и) = x5 (132) = 0,

Q = { и = { и 1 ,и ₂, . , u_n } , и_к е U , к = 1,2, . , N } , U = { w е R : 0 <  w < 50 } .

Функция Понтрягина имеет вид:

H ² ( v , x , w,t ) = E V i f ( ^x ( t ) , w , t ) + V 5 E Q ( x ( ¹ ) ) .

i = 1                                           j е { 1,3,4 }

Производная функции Понтрягина по переменной управления равна:

—

II

—¹— ( v x , + V ^Xi ) + ~ I 1 — x ³ 1 + ^ 4-500 x/          ’ x ₄ ( 500 J 500

Сопряженная система принимает вид:

• V 1 (t) = —V1 (t)f- — V2 (t )f — V3 (t )f — V5Q (x (t)) , V1 (132) = 0, dx1          dx1d

^V 2 ( t ) = ^— V 2 ( t ) ^f , V 2 ( ¹³² ) = ^X 4 ( ¹³² ) , d x ₂

• V 3 (t) = —V1 (t)f — V2 (t)f — V3 (t)f-— V5Q33 (x3 (t)), V3 (132) = 0, dx3         dx3

• V 4 (t) = —V1 (t)f — V2 (t)f — V3 (t)f- — V5Q4’ (x4 (t)), V4 (132) = x2 (132), dx4         dx4

V 5 ( t ) = 0, V 5 ( 132 ) = — 2 .

Производные в правой части системы вычисляются аналитически.

Задача о неподвижной точке принципа максимума с a >  0 имеет вид:

^U k = ^ 0;50] ( ^U k ⁺ a S T H u ( V ² ( t , ^U ) , ^x ( t , ^U ) , ^U k , t ) ^dt ) , ^k = ^1,2, ^ ^{, N} , (12) x 5 ( 132, U ) = 0.

Итерационный процесс метода неподвижных точек с индексом s >  0 для численного решения задачи (12) принимает вид:

usk+1 = P0-5Jusk + af H^V2(t,us),x(t,us),usk,t)dt), k = 1,2,...,N, k  [0;50] kTuk gs (2) = Ф1 (U+1 ) = x5 (132,U+1 ) = 0.

Численное решение задач Коши для фазовой и сопряженной систем осуществлялось с помощью подпрограммы DIVPRK из библиотеки IMSL языка программирования Fortran, реализующей метод Рунге — Кутты — Вернера [16]. Значения расчетных фазовых, сопряженных и управляющих переменных запоминались в узлах заданной равномерной сетки T_h с шагом дискретизации h = 10 ^{— 2} на интервале [ 0;132 ] .

Для поиска множителя Лагранжа 2 1 рассматривается вспомогательная задача g^s ( 2 ) ^ inf, численное решение которой относительно парамет

' '       21е R ра 2 е R осуществлялось с помощью подпрограммы DUMPOL из библиотеки IMSL языка программирования Fortran, реализующей метод деформируемого многогранника [17], с параметром ftol, равным 10"4. Досрочный выход из подпрограммы на s -й итерации для поиска множителя

Лагранжа Л₁ с итерационным индексом I >  0 осуществлялся при выполнении условия g^s ( Л ) <  £ 1 , где £ 1 >  0 — заданное малое число.

Критерий остановки поиска управления:

max uS+1 - uS < £ max uS , k =1,2,...,УI k         k'      2 k =1,2,.,N' kl где £2 > 0 — заданное малое число.

Проводились расчеты для разного числа N подынтервалов разбиения и разных начальных управлений; в представлении результатов расчетов использованы следующие обозначения:

— n_iter — количество проведенных итераций,

— M — максимальное нарушение ограничений.

Ниже представлены результаты расчета для N = 7 неравномерных подынтервалов разбиения; в качестве моментов переключения управления были использованы значения, полученные в [15]:

т = ( 4,07;11,12; 15,33; 19,04; 24,84; 54,9 ) ^т .

Начальное управление:

u

0 k

к = 3, к * 3,

к = 1,2,...,7,

где к = 3 соответствует подынтервалу 7 3 = [ 11,12;15,33 ] .

При начальном управлении получены следующие значения:

Ф 0 = - 98,80, Ф 0 = 0,68, M⁰ = 0,51, x ⁰ ( 132 ) = ( 32,82; 9,40; 0,0032; 10,51 ) T .

В таблице 1 представлены результаты расчета со следующими параметрами: £ ₁ = 10 ^{- 6}, £ 2 = 5 - 10 ^{- 6}, a = 3 - 10 ^{- 3}.

Таблица 1

n iter	Расчетное управление		Ф0 ( u )	ф1 ( u )	M	x ( 132 ,u )
	k	uk				i	xi
1182	1	12,50	-87,61	8 , 91 - 10 - 7	0,016	1	28,33
	2	12,51
	3	37,51				2	8,75
	4	12,51
						3	0,0014
	5	12,51
	6	12,09				4	10,016
	7	9,45

На рисунке 2 представлены графики расчетного управления и соответствующих фазовых траекторий.

— полученное управление и = {u\,U2,..., иу}, - - - начальное управление u^Q.

— траекториях^), --- траектория^(О, ........ траектория %з(г), — траекториях^).

Рис. 2. Графики управления и фазовых траекторий при N = 7

Ниже представлены результаты расчета для N = 66 равномерных подынтервалов разбиения. Интервалы вычислялись по формуле:

T =ГМ 3 1 к =   , N.

к N N

Начальное управление:

u_k ⁰

37,5,

12,5,

T. ^ [ 10;14 ] ,

T ф[10;14 ] ,

к = 1,2,..., N .

При начальном управлении получены следующие значения:

Ф 0 = - 98,42, Ф 0 = 0,63, M ⁰ = 0,50, x ⁰ ( 132 ) = ( 32,81; 9,37; 0,0033; 10,50 ) ^T .

В таблице 2 представлены результаты расчета со следующими пара- метрами: £1 = 10 6, £2 = 5 • 10 6, a = 2,12 *10 2.

Таблица 2

n iter	Расчетное управление		Ф 0 ( u )	ф1 ( u )	M	x ( 132 ,u )
	k	uk				i	xi
13996	1	12,77	-87,83	4 , 58 - 10 - 7	0,014	1	28,11
	2	12,87				2	8,77
						3	0,0018
	66	9,60				4	10,014

На рисунке 3 представлены графики расчетного управления и соответствующих фазовых траекторий.

полученное управление и = {wi, И2,..., Ыбб},             — траектория xi (7), начальное управление и°.                               ’'' траектория *2(0,

........ траектория %з(г), — траектория X4(f).

Рис. 3. Графики управления и фазовых траекторий при N = 66

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

Заключение

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

Указанные свойства метода являются важными для повышения эффективности поиска решения рассматриваемых задач.