CC BY
62
УЧЕНЫЕ ЗАПИСКИ ЦАГИ Т ом XXII 199 1
мз
УДК 629.735.33.016
МЕТОД .ВТОРОГО ПОРЯДКА ОПТИМИЗАЦИИ УПРАВЛЕНИЯ НЕЛИНЕЙНЫХ ДИНАМИЧЕСКИХ СИСТЕМ И ЕГО ПРИМЕНЕНИЕ ДЛЯ РАСЧЕТА ОПИМАЛЬНЫХ ТРАЕКТОРИЙ САМОЛЕТА
О. Е. Ефимов
Рассматривается методика численной оптимизации нелинейных динамических систем, о снованная на дискретном варианте прямого метода второго порядка. Улучшающая вариация управления находится из решения краевой задачи для линеаризированных прямой и сопряженной систем. Методика позволяет получить слабый локальный оптимум с высокой точностью, а также достаточно быстро построить семейство соседних экстремальных решений. Методика используется для решения задач оптимизации траекторий самолета, в которых учитываются ограничения на нормальную перегрузку и тягу, а также на некоторые траекторные параметры. Рассматриваются задачи вывода самолета в заданные кон ечные условия по высоте и скорости за минимальное время и с минимальным расходом топлива.
Известные градиентные методы первого порядка [1—4] применительно к задачам оптимизации управления нелинейных динамических систем имеют плохую сходимость в окрестности оптимального решения. Вследствие этого обычно трудно судить о том, насколько получаемое-этими методами решение близко к оптимальному. Эта трудность может быть преодолена с помощью методов второго порядка, применение-которых вблизи оптимума обеспечивает хорошую сходимость. К преимуществам этих методов можно отнести также возможность достаточно, просто контролировать не только необходимые, но и достаточные условия оптимальности, большую точность получаемых решений. Однако для использования этих методов, как правило, требуется выпуклое [1 ]■ начальное приближение. Это условие, а также значительные объемы вычислений и/или памяти ЭВМ, необходимые для расчета и хранения матриц вторых производных, сдерживают широкое применение этих методов для задач оптимального управления, в том числе для задач -оптимизации траекторий самолета.
В настоящей статье рассматривается модификация предложенного? в [1] метода второго порядка, которая позволяет учитывать ограничения на управляющие и фазовые переменные для задач с фиксирован-
ным и нефиксированным временем окончания процесса. Учет ограничений проводится с использованием множителей Лагранжа. При решении краевой задачи для линеаризованных прямой и сопряженной систем используется алгоритм обратной прогонки.
В качестве иллюстрации эффективности метода были рассмотрены задачи об оптимальном по времени и расходу топлива полете самолета в заданные конечные условия по высоте и скорости. Расчеты проводились для самолета, характеристики которого заданы полиномиальными и дробно-рациональными функциями, обеспечивающими непрерывную дифференцируемость правых частей уравнений движения.
I. Постановка задачи. Рассмотрим динамическую систему, описываемую дифференциальным уравнением
X =/ (х, в), х(О) = х0, О < Ь< Т. (1)
Здесь Х — п-мерный фазовый вектор, и— т-мерный вектор управления, Т—время окончания процесса, вектор х0 задан.
На движение системы могут быть наложены ограничения вида
<р (х, и) > 0,
где у — р-мерный вектор.
Задача состоит в том, чтобы найти вектор-функцию и (/) и время Т, максимизирующие функционал
т
У =Г (х (Т)) +\ Ь(х, и) dt
о
при условии
ф(х(Т)) = О,
г де Ф — 9-мерный вектор, причем 9 -< п — 1, если I == О и 9 <; п, если £фО.
При численном решении поставленной задачи на ЭВМ непрерывная система (1) заменяется конечно-разностными уравнениями вида
Х*+1 = £*(**, "*, I
Д^+1 = д**, (2)
к = 0, ... , Ы — 1, )
где вектор Х0 — задан; верхний и нижний индексы к соответствуют к-му узлу. Для получения условия оптимальности на шаг интегрирования Д* введено второе уравнение (2), из которого следует, что Д£* = Д*, к = О, ... , Ы, а последовательность Д**, к = = О, ... — 1, дополнена членом Д^. Векторыи* и х* должны
удовлетворять ограничениям вида
?*(■**. **а)>0. (3)
В дальнейшем будем считать, что все компоненты вектора ограничений (3) зависят от управления явно. Если какой-либо компонент <р» явно от управления не зависит, то, подставляя в него выражение для
X к из (2), получим <?[ (£*-' (Хк-1, и*-1, Д**-!)) > О, т. е. в дискретной
постановке любое ограничение на движение системы можно представить в виде
**(■**. Щу Д^)>0. (4)
Таким образом, требуется найти последовательность и0, •.. , «N-1 и значение М, удовлетворяющие условиям (4) и
Ф (х„) = 0 (5)
и максимизирующие критерий качества
N—1
У = Е(*„)+ Е ^(хк, «к-
(6)
к =0
Здесь сумма представляет собой дискретный аналог интегрального члена в исходном функционале.
2. Условия оптимальности. Используя множители Лагранжа, запишем первую вариацию критерия качества с учетом условий (2), (4), (5): ,
'О/ = Рх Чу + Е (^5** + ^ 8ик + ЬЪШ") + 'У' Фх 8х„ + к = 0
N—1 _ N-1
+ Е Лк+1 х 0хк + gu 8ик + gДt оЫк — 8хк+{) + 2 ^*+1 (8А!а — 8Д/к+0 + к=0 ' к=0
+ 2 (I)* (?* 8Хк + 9а 8»к + ?д< 8Мк - 8ф*).
Здесь Л, 11-, (1), 'У — множители Лагранжа, соответствующие условиям (2), (4), (5), верхний штрих обозначает операцию транспонирования.
Введем обозначение
Нк = Ьк(хк, ик, Мк) + Хк+^к(хк, ик, М//) + 11-к+1 М// +
+ (I)* ?к (Хк, ик, Дt//).
Полагая
получаем
Лк= нХ, 1
11-// = Н«, к = 0, ... , N -1,
^ = (Рх + 'У' Фх)" ,
^N = 0,
N-1 N-1
ЬJ— х0 8лг0 + 11-0 Ш0 4- Е н« 8и* —■ Е 0>* 8?*'-
(7)
й = 0
к=0
Гамильтониан Н" и соотношения (7) представляют собой дискретные аналоги функции Гамильтона и уравнений, определяющих сопряженные переменные, соответственно для непрерывного случая. Если последовательность «о, ..., «N-1 и № максимизируют критерий качества (6), то для всех допустимых 6иь, к=0, ... , N — 1, и ЬМ должно быть б!-<О. Можно показать [1], что при 8хо = О:
н: =о, ®А>0, ?к = 0,
ык = °» ?к > О, к = 0, ... , N— 1,
1-'-о = 0.
(8)
Следует отметить, что для задачи с фиксированным временем окончания процесса все соотношения остаются справедливыми за исключением последнего условия в (8), так как вариация шага интегрирования 6Ato = O и J.to свободно. Таким образом, для получения оптимальных последовательности и0, •.. , Un-i и At необходимо решить нелинейную-.1 двухточечную краевую задачу, описываемую соотношениями (2), (5), (7), (8). .
3. Вариация решения. Предположим, что имеется номинальное решение, которое удовлетворяет условиям (2), (5), (7), (8). Рассмотрим теперь малые отклонения от этого решения, возникающие вследствие малых возмущений в начальном состоянии 8хо, в конечных условиях 8ф, в ограничениях на движение системы 8^k, в условии оптимальности управления 8H* и в переменной Лагранжа, соответствующей шагу интегрирования, 6J.t0. Эти возмущения приведут к появлению вариаций &Xk, OAk, 6J.tk, БЮк, 6At, б'У и Oв//, удовлетворяющих линеаризированным в окрестности номинального решения уравнениям (2), (5), (7), (8):
OXft+i = g*0xk + gka Ьвк + gl№tk,
Mtft+i = Mtk,
8?k = ?x OXk + Ф* OBk + фд i Mtk,
OAk — H.*x SXk + gX OAhi + H*u OBk + HxmMtk + фх
°t*k — 0fAk+i + Hitx^Xk + H:ta OUk + Hit f.t OAtk + g:t OAhi + фд< 8Wki
(9)
3И* = Н** 8х* + И*„8и* + н£д, ОА ** + $ 8Ан1 + 9и Ощ, к = О, ... , N-1,
Оф ^^х^
8А* = (Рхх + 'У' Ф**У 8х* + (Фх)' 5'1,
8И„ = °
Варьирование условий оптимальности (2), (5), (7), (8), следующих из рассмотрения первой вариации критерия качества, соответствует рассмотрению второй вариации критерия качества (6).
Представим уравнения (9) в более компактном виде:
8^+1 = А* 8у* + а* 8«*, (10)
8гс* = А*' оя*+1 + В\8у* + Ь* 8«* + а*' 8шъ (11)
оНИ = Н*а 8«* + Ь* 8у* + а* 8я*+1 + ^ 8ю*, (12)
= ?и 8«к + а* 8уь (13)
8яа = ^ 8ум, (14}
где
gi gbt Qnq
Qm 1 Qlq
Qqn Qg 1 Igg
\
gu \
Q1 m ,
Qqm J
m Н*д, Q
B>!=={ Нш QXl
Qqn Qq 1 Qgg
dk = (9* фд< Q
pq!
tqm
.... Рл 1 ^ \
^ = ( Q.„ 0 Q1? , (15)
Qg 1 Qqq /
/mm — единичная матрица размеров тХт, Qmi—матрица размеров т Х/, все элементы которой — нули.
Соотношения (10)—(14) описывают линейную краевую задачу относительно Syk, 81:ъ 8ик, В начальной точке заданы компоненты 8х0 вектора 8у0 и 8[10, Оф вектора 81:0, в конечной — граничное условие (14). Для решения линейной краевой задачи это условие переносится в начальную точку с использованием метода
обратной прогонки [1]. Предположим, что оно от шага к шагу
трансформируется согласно выражению
81:k = Sk 8yft + hft- (16)
Подставим уравнение (10) в (16):
81:k+i = Sk+i Ак8^к + Sk+i ak 8«ft + Ajt+i.
Затем это соотношение подставим в (11) и (12). Получим
01:k = Zyy O).'k + Zyu 8«k + dk' 8«k + Ak' hk+i, (17)
SH = ZUu 8«k + ZUy 8yk + СРИ 8wk + a,k hk+i, (18)
где
Zyy = Bk ~j- Ak Sk+i Ak,
ZyU = bk + Ak' Sk+lak,
zkau = HL + ak' Sk+1 a* . ^
/ 9k »*''
Предполагая, что матрица ( ““ I невырождена, разрешим
\¥* 0 /
уравнение (13) и (18) относительно Ьак, Ъык:
(Ьи* u _ [VК -1zL V1 (ak'hk+1 - 5Н^ (19>
[*>) U 0) [d> I у* 0 ) [ -V У V '
С помощью соотношения (19) исключим Ьак и шк из (17):
8яь = J zk
к — \*W
\Ги о
1 / Л'
Ak' hk+,-{zkyudk')(Zua ) Г Нк+1 ~Ш“). (20>
. ?и
Отсюда следует, что представление (16) справедливо при
............................. } (21)
А// = А А*+1
По этим рекуррентным соотношениям, используя условия (15) и й№ = 0, можно получить 5с и Ас" Условие (16) при к = О можно разрешить относительно неизвестных компонентов 81.0, 8Д^, 8.., векторов 8уо и 8яс, а затем получить решение линейной задачи (10)-(14), по рекуррентным соотношениям (10) с использованием (19).
4. Итерационная процедура. Предположим, что имеется последовательность «*, к = О, ... , N -—1, величина 1: t и удовлетворяющая (2) последовательность х6, к =1, .. , которые близки к оптимальному решению в том смысле, что начальное условие и условия (5), (8) выполняются с небольшими погрешностями Дх0, Дф, Ду*, 1:Н:, Д110- Тогда, полагая
8х0 — — Дх0, 8ф = — Дф, 8у* = — Ду&,
8И* = - ДН* , Д11о = -Д!,-
и пользуясь результатами п. 3, можно получить улучшающую вариацию решения. Если полученные в соответствии с (10) вариации 8х к k= 1, ... ,#, сложить с исходной последовательностью X*, k= 1, ... ,#, то вследствие конечности погрешностей уравнения (2) точно удовлетворяться не будут. Поэтому вариации 8х6, которые подставляются в ('19), предлагается вычислять не по первой формуле (9), а как разность
к = 1, .,. ,^,
где 8хо = 8хо, а 8«к— находится по формуле (19). Для получения приближения в том случае, если погрешности Дхс, Дф, Ду*, ДН* , Д110 не являются малыми, необходимо ограничить вариации 8хо, .оф, 8у*, 8И*, 8!'-о.
Например, вариация б^хо может быть задана следующим образом:
' 8!'-о шах при 1:110 < — 8!'-о шах ,
8!'-о := 1:1А0 при | 1:11о1-< 811Ошах,
, 8!'-Ошах при > 8!'-Ош1х.
Аналогично можно определить компоненты векторов 8хо. оф, 8^к , ЗН* . Предельные значения вариаций выбираются так, чтобы они не приводили к значительным отклонениям вариаций решения от линейных представлений. Подобную' процедуру можно повторять до тех пор, пока условия (5), (8) не будут выполнены с заданной точностью.
Одно из ограничений применения описанного метода состоит в том, что начальное приближение должно удовлетворять одному из достаточных условий оптимальности — условию выпуклости. В частности, при
отсутствии ограничений (4), оно означает отрицательную определенность матриц г£и, к= 1, ... ,М — 1. Для получения выпуклого начального приближения можно использовать градиентные методы первого порядка [1], например, метод проекции градиента [4], либо воспользоваться одним из вариантов метода обратной прогонки, предложенным Д. Н. Джекобсоном [5] —дифференциальным динамическим программированием (ДДП). В работе был применен подход, связанный с построением такой короткой траектории, вообще, говоря, с другими начальными и конечными условиями, на которой выполняется условие выпуклости. Искомое решение можно получить затем варьированием начальных и конечных условий, т. е. методом соседних экстремалей.
Решение задачи (2), (5), (7), (8) может привести к значению Ы, при котором точность аппроксимации конечно-разностными уравнениями (2) исходных уравнений (1) оказывается неудовлетворительной. В этом случае задачу (2), (5), (7), (8) необходимо решать при новом значении количества шагов N соответствующем шагу А/, при котором обеспечивается приемлемая точность аппроксимации. Такое решение можно получить на основе прежнего с использованием (19).
5. Задачи об оптимальном по времени и расходу топлива полете самолета в заданные конечные условия по высоте и скорости. Рассмотрим применение описанной методики для задач оптимизации траекторий самолета, движение которого описывается традиционными уравнениями движения его центра масс в вертикальной плоскости, в которых фазовыми координатами являются величина скорости, угол наклона траектории к местному горизонту, дальность, высота и вес самолета. В качестве управляющих переменных выбраны коэффициент аэродинамической подъемной силы, который ограничивается своим максимальным допустимым значением и значением, соответствующим максимальной эксплуатационной перегрузке, а также тяга двигателя, ограниченная тягой малого газа снизу и тягой, соответствующей режиму полного форсажа, сверху. Учитывается ограничение по минимальной высоте полета.
Были получены аналитические выражения для первых и вторых производных от правых частей уравнений движения по фазовым и управляющим переменным, что существенно повысило быстродействие и точность расчетов. В качестве метода интегрирования был выбран метод Эйлера, т. е. в формуле (2)
ё*(хк, ик, А^) = хк + Мк/(хк, и*).
Расчеты проводились для самолета с большой тяговооруженностью, в качестве начальных были выбраны малая высота и малая скорость. Конечная высота была зафиксирована, а конечная скорость V/ варьировалась.
Текущие величины М и Н на рис. 1, 3 отнесены соответственно к Мтах и Нтах, представляющим собой максимальные значения числа М и высоты Н установившегося горизонтального полета. Стрелками на всех рисунках показано направление увеличения конечного значения скорости V/. На рис. 1 представлено семейство оптимальных траекторий для задачи, в которой функционалом являлось время полета самолета. Оно обладает закономерностями, сходными с полученными в (4], где решалась аналогичная задача с использованием метода проекции градиента, и соответствующими структуре решения вырожденной задачи оптимизации [2, 3]:
Рис. 2
— на всех траекториях реализуется максимальная форсажная тяга, а полная механическая энергия монотонно возрастает;
— оптимальные траектории формируются из трех основных участ-
ков: участка входа в окрестность некоторой базовой кривой, которая является аналогом особого решения вырожденной задачи оптимизации [2, 3], движения вблизи этой кривой и схода с нее с приведением самолета в соответствующие условия по высоте и скорости. Отметим, что зависимости коэффициента аэродинамической подъемной силы от времени (рис. 2) являются гладкими и носят регулярный характер при изменении краевых условий, что является результатом использования принятых метода оптимизации и гладких аппроксимаций характеристик атмосферы и .самолета. ,
На рис. 3 представлено семейство оптимальных траекторий для задачи, в которой функционалом являлся расход топлива. Это семейство также характеризуется монотонным возрастанием полной механичес-
кой энергии и наличием базовой кривой. Для полученных оптимальных траекторий зависимости коэффициента дросселирования (отношения тяги к ее величине, соответствующей режиму полного форсажа) от времени (рис. 4) практически совпадают между собой, а зависимости коэффициента аэродинамической подъемной силы от времени (рис. 5) различаются только на конечных участках, соответствующих переходу от базовой кривой в конечные условия. Следует отметить, что, начиная с некоторой конечной скорости, реализуются отрицательные значения конечных углов наклона траектории, т. е. на конечном участке полет происходит со снижением. Максимальная высота полета достигается в промежуточной точке, а ее значение и длина участка снижения увеличиваются с возрастанием конечной скорости. Как и в предыдущей задаче, на всех полученных траекториях зависимости управления от времени (см. рис. 4, 5) являются гладкими.
ЛИТЕРАТУРА
1. Б райс о н А., Х о-Ю-Ш и. Прикладная теория оптимального управления. — М.: Мир, 1972.
2. И л ь и н В. А. Оптимальные режимы движения летательных аппаратов. (Вариационные задачи механики полета), обз ор БНИ ЦАГИ, N 40, 1961.
3. Исследование оптимальных режимов движения ракет. Сб. переводов иностранных статей. — М.: Оборонгиз, 1959.
4. Про к о п е ц Н. Б., С а в е л ь е в А. В. Приближенный синтез управления самолетом для набора высоты с заданной конечной скоростью за минимальное время. —- Ученые записки ЦАГИ, 1986, т. 17, № 4_
5. J а с о b s о n D. Н. New second order and first order algorithms for
determining optimal control: а differential dynamic programming ар-
proach. — J. Optimization Theory and Applkatoin (Dec., 1968).
Рукопись поступила ЗО/ІІІ 1989 г.
