Автоматика и телемеханика, № 2, 2014
© 2014 г. Ю.Я. ОСТОВ, канд. техн. наук А.П. ИВАНОВ, канд. физ.-мат. наук (Санкт-Петербургский государственный университет)
МЕТОД ОПТИМИЗАЦИИ В ЗАДАЧЕ ДИНАМИКИ ПОЛЕТА
Рассматривается и численно решается задача субоптимального управления летательным аппаратом в атмосфере. Уникальность работы состоит в том, что решение сложной нелинейной краевой задачи получено с помощью конечного числа арифметических операций. Приведенные результаты численного моделирования подтверждают эффективность рассмотренной методики.
1. Введение
Оптимальное управление, найденное как результат решения краевой задачи на основе принципа максимума в его классической формулировке, является программным управлением, и при наличии всякого рода возмущений оказывается неэффективным, т.е. не обеспечивает оптимум заданного критерия качества. Поэтому, если говорить о прикладном аспекте решения задачи, достаточно построить управление, при котором значение оптимизируемого функционала отличается от его оптимального значения не более, чем на заданную величину. Такое управление будем называть субоптимальным управлением.
Новый вариационный метод применен к решению следующей задачи: оптимизировать траекторию продольного движения центра масс (ЦМ) летательного аппарата (ЛА), совершающего полет из начальной точки атмосферного пространства в заданную конечную точку на поверхности Земли. Критерием оптимальности траектории (управления) является максимум кинетической энергии ЛА в конечной точке траектории.
2. Уравнения движения летательного аппарата
Движение ЦМ ЛА описывается уравнениями
IV X — д вш в.
- =--
и т
<1в У д сов в
mV V
0(0) = во,
(!) сШ
^ = Я(0)=Яо,
где V(Ь) - модуль скорости, 9{Ь) - угол наклона траектории к горизонту, Н(Ь) - высота над поверхностью Земли, £(Ь) - безразмерная "взвешенная" длина траектории, Ь - текущий момент времени, Ь £ [0,Т], т - масса ЛА, д -модуль ускорения свободного падения, X = (Сх0+CXia2)pV2Б/2 - лобовое сопротивление ЛА, У = C'aapV2Б/2 - подъемная сила ЛА, р = р£ ехр(—(ЗН) -плотность атмосферы на высоте Н, р = р/р0, а - угол атаки (управление).
В данной модели Схо, Схг, С£, Б, т, р0, в, д - заданные константы. Множество А допустимых значений угла атаки а является открытым.
Для этой модели ставится задача: максимизировать кинетическую энергию
(2) Ла) = <р(У{Т)) =
при заданных значениях высоты и "взвешенной" длины траектории полета в конечный момент времени Т, т.е.
(3) Н (Т) — Нт = 0, £ (Т) — £т = 0,
где Нт и £т - заданные константы. Угол 9(Т) и конечный момент времени Т не фиксированы, что не влияет на общность подхода к решению задачи.
На первом этапе решения поставленной задачи исходная система (1) преобразуется к виду [1]
Щ = ~ (+ СуЪ, = Ую,
с2
'У
d,£ "у'1 \ ' ■ ' 2<1
(4) ^ = суЦ - сх0 + --. 9 , Vh(£0) = Vh0,
P\]V? + VÏ
dp Vh
# y/vF+vi?'
p(Cû) = Po,
где
P = "4 = exP h = 9 = ê>
po P
_ _ CxoPqS (CayfplS _ (C»)p*0Sa
0x0 ~ 2mp ' ~ 4т[ЗСхг ' Cy ~ 2mp '
а Vi и Vh - проекции скорости V ЦМ ЛА соответственно на горизонтальную и вертикальную оси инерциальной системы.
Соответственно функционал (2) и ограничения (3) перепишутся в виде:
J(u) = 0{(V2 + V2)(£t )) = V 2(£t ), р(£г) - pT = 0, Ç G [£o,£t], Ço = 0, £t > 0.
Как следует из приведенных выше соотношений, ЛА обладает параболической полярой, т.е. коэффициент лобового сопротивления связан с коэффициентом аэродинамической подъемной силы соотношением сх = сХо + су2/21, где Сх о и (I - положительные константы. Из последнего соотношения следует У = §л/ж ~ а1 гДе приняты обозначения у = су, х = сх, а = схо, Ь = у/2/(I.
Обозначим С(х) = \у/х — а. Пусть ж д. - коэффициент лобового сопротивления. Касательная к кривой С(х) в точке х = хи представляется в виде:
_ dG ^ дх
(х - хк) + G(xk) = + 2ЩЛ.
x=xk by/xk — a b
Обозначим точку пересечения этой прямой с осью Ox через z, т.е.
0 = (z -Xk) 2 у/Хк - а Ьу/хк — a, b
Откуда следует
(5) xk + z = 2a, a = const.
С учетом последних двух соотношений уравнение касательной принимает вид
(x — z)
У = К /--
by/a — z
Точку z будем называть сопряженной точкой по отношению к xk. Таким образом, точке параболы (cx = x,cy = y) сопоставлено уравнение прямой (в пространственном случае - уравнение гиперплоскости, которое соответствует преобразованию Лежандра применительно к параболоиду).
3. Оценочная система и метод оптимизации
Теперь система уравнений (4) перепишется в виде
^ = -ъц - lvh, v;m = vl0,
by/a — z
(6)
dV± = i£^jlVi _ -CxVh _ 9 yh{(o) = у
bsj^z P^Vf + V* dp Vh
Jv,2 + K2
P(£ o) = po,
Система уравнений (6) при переходе к новой независимой переменной р примет вид
йр Ун Ъу/а — г
2 , т/2
(7) дУн _ (сх-гЩ^У? +
(1р ЬУну/а-х V ' п рУн
№ + ^
йр Ун
Гамильтониан системы уравнений (7) и сопряженная система уравнений теперь запишутся так:
(8) я = + +
ар ар ар
¿Фо дН
-—— = -Ф0 Фо - Ф1 - Фг Фг,
а р дУг
н
А
аФ2 дн
(9) = -ЁЕ. = -Ф0Ф3 - Ф:Ф4 - Ф2Ф5,
ар дУн
2 дН
о,
где
а р д£
рхУ рхV2 (Рх - г)Уг
срп = - - - —
Ун У/г^ _ (сх — х)У (ст-г)^2 схУ1 1 ~ ЬУну/а^х ЬУУнф1^~х V ' ф = _А_ ф = . (ёх ~ *)Ун
2 УУн 3 УУ£ ЪУу/а=х'
(ст - х)У? схУн д У,2
ЬУУ?л/сГ^~г У рУ'н 5 УУ'н
Так как V = \^У2 + > 0, то условие стационарности гамильтониана (8) относительно управления Рх можно представить в виде
(10) * = = + 1 ')+»,(_«= + !)=„. У<9ст Ьл/а -г/ \ ЪУн\[а — ^ )
Следует отметить, что соотношение (10) есть необходимое и достаточное условие оптимальности Ру в исходной задаче согласно принципу максимума, т.е. задача оптимизации управления Рх в системе (7) регулярна (не вырождена). Систему (7) будем называть оценочной снизу для исходной системы (1).
Условие стационарности гамильтониана (8) относительно г(р) при ограничении (а — г) > 0 можно записать так:
1 дН^'а-гУ/г = ( Ф0 -Ф1
V dz
2 b(a - zf2 = ^Ф0 - (2а - с, - г) = 0.
Первый сомножитель в правой части этого выражения отличен от нуля (иначе с учетом соотношения (10) Фо = Ф1 = 0), второй сомножитель - модификация соотношения (5) - равен нулю.
В силу нестационарности системы (7) гамильтониан (8) на экстремали является функцией р, т.е. H = H(р). Из соотношений (6), (7) следует, что
Н(р) = —C^jvf + V^/Vh, где С = const. > 0. Поэтому с учетом соотношения (10) выражение (8) можно преобразовать к интегралу следующего вида
(11) + . =0-
Дифференцируя интегралы F1(X(p),z(p)) = 0 и F2(X(р), z(p),p) = 0 по р, где X = (Vi,Vh, Ф0, Ф1) в силу системы (7), (9) и исключив переменную dz/dp, получим интеграл F3 = 0 в виде:
F3 = an й>22 — ai2 &21 = 0,
где
dF1 dX dF2
0,11
dX dp dz
dF1 dF2 dX dF2
dz dX dp dp
В развернутом виде интеграл F3 = 0 можно представить так:
5Ф1 , . уV?
(12> F' = WV + { bv; +фщ)с = 0'
Переменные Фо и Ф1 выражаются из уравнений (10), (11):
Фо = Ф1 =
zV3 p + gb s/a - z Vi + Vhg zV3 p + gb sja, - z Vi + Vhg
Подстановка последнего соотношения в интеграл = 0 с учетом соотношения (5) приводит к уравнению 3-й степени относительно сопряженной переменной г:
(13) Р = Рг у/а^г + Ро = 0,
причем можно положить параметр Ф2 = С + С, а С = 1. В этом случае Ро и Р\ примут следующий вид:
(14)
, „cos2 в a cos2 в sin QC\ „cos2 в a
Po(z) = 3 —-2- г + _у . +-- 2 . ,
sin2 в pV2 sin в a sin2 в
„ , . a V2 cosв a b cos3 в cosв b cosвCl
где 0 = arct.g ( ), F = ,/К2 + V,2, a = cx0-
sin2 вЬд sin^V^ sinвb p
r2 , V2
^ <9Я dX
Третий интеграл может быть вычислен по формуле F¿ = —--— = 0.
dX dp
В этом случае уравнение относительно переменной z принимает следующий вид:
(15) Р = А + Ро = 0,
где полиномы Po и А определяются выражениями:
„ 3 cos2 в д cos2 в C1 sin в 2a cos2 в
р0 =___—z------1--1____
sin2 в pV2 sin в p sin2 в
2cos0 (62Ci sin0 + —^
„ . a / Cibcos в
Pl = - , • 2л-~ - +
b sin2 в p
2 cos 9 2 a b С i cos 9 bg cos3 9
b sin 9 sin 9 V 2p sin2 9
Подстановка z = a—Ir2 приводит соответственно к полиномам третьей степени
P (Ir ) = Pl(Ir )Ir + Po (Ir ), P(Ir ) = A (Ir )Ir + Po (Ir ),
имеющим общий положительный корень I°, который находится за конечное число арифметических операций. Синтезированное по упрощенной модели управление а находится по формуле
<ЛРЛ мто 2тв
(16) а =
при текущих значениях фазовых переменных 9, р, V. Знак выражения (16) зависит от граничных условий (3).
Уникальность работы состоит в том, что решение сложной нелинейной краевой задачи (при подобранных константах) получено с помощью конечного числа арифметических операций.
4. Результаты численного моделирования
Результаты численного моделирования, приведенные ниже, подтверждают эффективность рассмотренной методики. Система уравнений (1) интегрировалась методом Рунге-Кутты (4-го порядка точности на шаге) с шагом hu = 0,01с при следующих значениях параметров ЛА и атмосферы: Cx0 = = 0,1931; Cxi = 5,88; = 0,8548046; m = 422; S = 0,159; g = 9,81.
Вариант 1: p*0 = 2,047; в = 1,5682 • 10"4.
Начальная точка здесь £(0) = 0,00; H(0) = 24054,5; V(0) = 1088,31; 9(0) = = -0,54113.
Терминальные ограничения (3), (4) следующие: H(T) — 0 = 0; £(T) — - 1,42500 = 0.
Результаты счета для этого варианта приведены в табл. 1.
В этом варианте управление a(t) < 0, t G [0; 37,88].
Второй вариант отличается от первого начальной точкой траектории: £(0) = 0,0; H(0) = 30000,0; V(0) = 1200,0; 9(0) = —0,4, а также ограничением (4) на конечную точку траектории: £(T) — 1,6500 = 0. Остальные параметры соответствуют варианту 1. Результаты счёта варианта 2 приведены в табл. 2. В этом варианте a(t) > 0, t G [0; 50,93].
В табл. 3 приведены результаты решения задачи для первого варианта (табл. 1), полученные с использованием уравнения третьей степени (табл. 3, п. 1), и для для этого же варианта с использованием линейного уравнения (табл. 3, п. 2).
Таблица 1
Параметр Упрощенная модель Полная модель
Ъ 13,375710 —
Сг 0,422616910 —
Т 37,88 37,88
t(T) 1,42500 1,42500
Н(Т) 0,00000 0,00000
V(T) 700,84704 700,8499
в(Т) -0,844170 -0,844283
Таблица 2
Параметр Упрощенная модель Полная модель
Ъ 4,87190 —
Сг 0,7479072 —
Т 50,93 51,01
t(T) 1,6500 1,6500
Н(Т) -0,00001 0,00000
V(T) 695,8655 695,93738
в(Т) -0,72382 -0,72557847
№ п/п V(T) в(Т)
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.