CC BY
22
УДК 531.3:681.5
Н.Е. Зубов1,2, Е.А. Микрин1,2, С.С. Негодяев2,3, И.Н. Лаврентьев1
1 Ракетно-космическая корпорация «Энергия» им. С.П. Королёва 2 Московский физико-технический институт (государственный университет)
3 Центральный научно-исследовательский институт химии и механики им. Д.И. Менделеева
Синтез управления сближением космического аппарата с полярной схемой двигательной установки по методу свободных траекторий на основе алгоритма оптимального управления с прогнозирующей моделью
Рассмотрена задача сближения космических аппаратов, находящихся на круговых орбитах, с использованием алгоритма оптимального управления с прогнозирующей моделью при полярной схеме построения двигательной установки. Получены простые аналитические формулы для вычисления управлений.
Ключевые слова: сближение космических аппаратов, оптимальное управление, метод свободных траекторий, прогнозирующая модель.
Будем рассматривать сближение маневрирующего космического аппарата (КА) с пассивным аппаратом, находящимся на круговой орбите. Пусть на активном КА постоянной массы установлен по оси х связанной системы координат один двигатель постоянной тяги для управления движением центра масс.
В этом случае траекторное управление КА будет осуществляться за счет ориентации вектора тяги в выбранной системе координат, и, соответственно, этот вектор через углы ориентации а\ и а2 запишется так:
a cos ai cos a a sin ai a cos ai sin a2
(i)
где а — модуль постоянного ускорения двигателя.
Уравнения относительного движения КА для круговой орбиты в орбитальной системе координат будут иметь вид [1]:
X = Ax + Ba, где X = [xi, X2, ..., Хб]Т,
A=
(2)
0 1 0 0 0 0
0 0 0 2ш 0 0
0 0 0 1 0 0
0 -2ш 3ш2 0 0 0
0 0 0 0 0 1
0 0 0 0 -ш2 0
0 0 0
1 0 0
О О 0
B=
0 1 0
0 0 0
1 0 0 1
x хз = У, Х4 = У,Х5 = z,
фазовые координаты активного КА,
ш = const — орбитальная угловая скорость пассивного КА.
Аналогично, как ив [1], предположим, что моменты времени начала маневра сближения ¿о и его окончания ¿ь заданы, причем ¿ь — ¿о фиксировано. Управление движением ведется по методу свободных траекторий [1] и состоит в определении первого программного импульса в момент ¿о, предназначенного для обеспечения встречи активного и пассивного КА за заданный промежуток времени ¿ь — ¿о при их дальнейшем свободном движении. Второй программный импульс в конце маневра сближения выравнивает скорости КА. Возможны дополнительные программные импульсы в промежуточных точках траектории. Управление носит импульсный характер, что позволяет находить величину корректирующих импульсов из уравнения (2) при а = 0 и соответственно углы «1, а2 ориентации импульса. Но поскольку коррекции реализуются только с помощью ограниченной по величине тяги двигателя, точное выполнение условий встречи таким образом невозможно. Аналогично [1] предлагается рассматривать расчет поправок на конечную тягу двигателя как задачу поиска оптимального управления с использованием алгоритма с прогнозирующей моделью [2]. Однако в отличие от [1], где алгоритм управления с прогнозирующей моделью представляет собой чисто релейное управление, в данном случае он будет сочетать в себе совместно как релейное, так и непрерывное управление. Особенности такого алгоритма совместного непрерывного и релейного управлений с использованием прогнозирующей модели рассмотрены в [3]. Здесь в отличие от [3] в силу наличия аналитического решения системы (2) получим редакцию алгоритма управления объектом с непрерывными и релейными рулевыми органами в редакции с аналитическим решением. В качестве релейных исполнительных органов будем рассматривать двигатель постоянной тяги, положение которого (в смысле «включено», «выключено») характеризует вектор 9Р, компоненты ко-
торого выражаются следующим образом:
м
др{Ь, tl, ..., Ьм) = —а ( — — (3)
м=1
где Ь — текущее время, 1М(Ь — ¿м) — единичная ступенчатая функция, соответствующая ¡л-му по порядку изменения компоненты, вр, — момен-
ты времени начала или конца импульса, М — число скачкообразных изменений функции (3). Если N — число импульсов в методе свободных траекторий, то М = 2N.
В качестве компонент управления релейных исполнительных органов рассмотрим моменты времени функции (3). Для них можно записать = ир.
Роль непрерывных рулевых органов будут выполнять углы ориентации вектора тяги а^ и «2№, и, соответственно, уравнение управления для них выглядит так:
а = ин.
Здесь а = [а^,а2м]Т.
Объединим компоненты вектора релейного управления вр и непрерывного а в единый вектор управления у = [Ь^,а]т и запишем объект управ-
ления (2) в виде
X = f (x,y,t), y — u.
(4)
Для минимизируемого функционала, заданного в виде [2]:
í k
1 — ^зад^^к ),tfc] +
^зад(х(г ),Т )dr+
í k
+0,5
[uT(т)K хи(т)+и„пт(т)K хиопт(т)]dT, (5)
при наличии аналитического решения первого уравнения системы (4) вида
х (х,у,г,гк)
оптимальное управление и определится выражением
и(Ь) = иопт(Ь) =
-K {
ífc ' д
д T д T
-^—XyX^y^Qytk ) ^зад \x j У ->t к ) +
+
Q^X(x,y,t,T)
д
~0^Q зад(х ,У ,т)
dT}. (6)
Соответственно при Qзад — 0
u(t) - uon(t)
г я
= -К{
д T д
X\Х->У-¡toytk ) ~о^Уаа Áx>y¿k)
}.
Здесь функции Узад, Qзад выражаются квадратичными формами:
Узад[х(Ьк ),Ьк] 0,5х (Ьк )рх(Ьк ),
Qзад[x(r ),т] = 0,5хТ (т )вх(т ),
р, в, К_1 — положительно заданные квадратные матрицы параметров.
Выражение вектора тяги (1) с учетом (3) будет иметь вид
У~] -a^2(-1)М1М(t - t^N) cosaiN cosa2N
N=1 м=1
N2
— a ( —1)M1^(t — ^N) sill aiN
N=1 m=1
N2
53 -^53 (-1)M1M (t - t^N )cosa1N sin a2N
N=1 m=1
(8)'
Следовательно, система (2) на основании (3) и (8) запишется так:
+
X1 X1
X2 X2
X3 — A X3
X4 X4
X5 X5
. Xe Xe
+B
N N N =1 N
-a^2 (-1)M1^(t - t^N) cos a1N cos a2N
¡1=1
2
'Sy^j ~a^^j (-1)MlM(t - ^N ) sin a1N
N =1 N
E -
N=1 m=1 N2
N =1
a 2_^(-1)M1M(t - t^N) cos a1N sin a2N ^=1
t1N ut1N
t2N ut2N
a1N ua1N
a2N ua2N
(9)
При N = 2 размерность вектора и, входящего в (9) равна 8, при N = 3 — соответственно 12 и так далее.
Решение первого уравнения системы (9) при [Ь1Ы,Ь2М,а 1м,а2м]Т = 0 и произвольных начальных условиях х(Ьи), вр(Ьи,Ьрн) имеет вид
х1(Ьк ) х3 (Ьк ) х5 (Ьк )
х2(Ьк)
х4 (Ьк ) х6 (Ьк )
0
1 6(т - sinт)
0 4- 3cos т
0
0 6и(1 - cos т) 0 3и sin т
0 0
0
0
cost
0
0
-и sin Т
-К1-cost) 0
4 cos т - 3 - 2sinr 0
'1 - cos T) 0
- sinT W 0
0 - sin T W
2sinT 0
cos T 0
0 COST
I
I
х1 (Ьи ) х3 (Ьи) х5 (Ьи) х2 (Ьи) х4 (Ьи) х6 (Ьи)
N
ф1911) Ф1д12) 0
Ф1(з1) Ф1(з2) 0
0 0 Ф1(5з)
Ф1(21) Ф1(22) 0
Ф1(41) Ф1(42) 0
0 0 Ф6з
+ Е
N =1
а1(Ь — tlN) соя аlN соя а2N а1(Ь — tlN) бш аlN а1(Ь — tlN) соя аlN £>ш а2N
N
—
N=1
Ф2д11) Ф2д12) 0
Ф2(з1) Ф2(з2) 0
0 0 Ф2(5з)
Ф2(21) Ф2(22) 0
Ф2(41) Ф2(42) 0
0 0 Ф2(6з)
а1(Ь — t2N) соя аlN соя а2N а1(Ь — t2N) зш аlN
а
1(Ь — t2N) соя аlN 8Ш а2N
(10)
где т = ш(Ьк — Ьи), т1N = ш(Ьк — t1N^
T2N = ш(Ьк — t2N),
4 3
Ф1,2911) = ^2 “ СОЭ 7"127У - £Т1,2м)>
2
Ф1,2д12) = — (тl,2N -втгх^дг),
функционала, которую определим квадратичной функцией вида
16
^зад [х(^к ^ ' Рг \ ''/ (I /, ,) хг/с]^? (12)
г=1
где Рг — весовые коэффициенты, хгк — заданные конечные значения фазовых координат КА. Функцию Qзaд положим равной нулю.
В соответствии с (7) при функции Узад, определенном выражением (12), с использованием записи аналитического решения (11) выражение для оптимального управления будет иметь вид
опт и
(Ьи
дХТ (Ьк ,tu,tlN,t2N,аlN,а2N )
-А--------------------------------X
ду
Р1[х1(Ьк) — х1к ]
Р2 [х2 (Ьк ) — х2к]
Рз[хз(Ьк) — хзк ]
Р4[х4(Ьк ) х4к]
Р5[х5(Ьк ) х5к]
Р6[х6(Ьк ) х6к]
После вычисления частных производных в развернутом виде получим
N
Е
N =1
Ut1N
Ut2N
Ua1N
Ua2N
N
—К
N=1
ф1,2(21) = (4 в1п 7”! 2ТУ - ЗГ1 2Лг),
Ш
2
Ф 1,2(22) = —(1 ~ СОвТх 2ы),
4 ' Ш
2
ф1,2(31) = — (вШТ! 2ЛГ — Т"1,2Дг), Ш2
Ф1,2(32) = —(1 - СОЭП^Дг),
Ф 1,2(41) — — (СОЭГ! 2ДГ — 1), Ш
Ф
1,2(42) — — вт Г1,2ДГ, Ш
Ф1,2(53) = ^2 (1 - СОЭП^Дг),
ф1,2(63) — — втГ12ЛГ.
Ш
В компактном виде аналитическое решение (10) запишем так:
X [х(Ьи),Ьк ,Ьи,Ьш ,t2N ,аш ,а2N ] =
N
Ход[х(Ьи),Ьк,Ьи]+) ^ Хч(Ьк,tu,t1N,t2N,а1N,а2N)-N=1
(11)
Зададим функционал (5) в соответствии с требованием точного выполнения условий встречи за заданное время сближения. За выход в заданную точку встречи «отвечает» терминальная часть
ргш
pt2N р а^ P’a2N
ргш
р t2N
ра1N
P’a2N
Pзt1N
р t2N ра1N ра2N
р t2N ра1N ра2N
P5t1N
р t2N ра1N ра2N
ZPt1N -р6 т?t2N р6
^pа1N
р6
тpа2N
р6
Р1[х1(Ьк) — х1к ] Р2 [х2 (Ьк ) — х2к] Рз[хз(Ьк) — хзк ]
Р4[х4(Ьк ) х4к]
Р5[х5(Ьк ) х5к]
Р6[х6(Ьк ) х6к]
(13)
где К
квадратная диагональная матрица,
К = dІag(fct1N ,kt2N ,каШ ,kа2N ),
77It1,2N а /о Л \
р1 = — СОвацу С08а2лг(^*1,2лг “ 4 81ПГ*1?2Лг) +
р
t1,2N
р
t1,2N
2а
н----81Паш(С08Ге1 2М - 1),
Ш
= а соя аlN соя а2N (3 — 4 соя TtlJ2N) — —2а siп аlN вш TtlJ2N,
2а
= --СОв СКц\Г ^2ЛГ (1 — СОв Та 2аг ^ —
Ш
----БШТа ,2ЛГ?
Ш
77lt1,2N г>
р = 2а соб аlN соб а2N бш TtlJ2N —
—а siп аlN соя TtlJ2N,
т-It1,2N а
Г К =-сов СКЦ\Г 81П (У-ЧЫ ^ШТа 2ЛГ,
5Ш
т-It1,2N
р = —а соб аlN бт а2N соб Ttl?2N,
Plа1N = —аsiп аlN соя а2N^1+
X
X
X
X
X
X
X
X
1
а
2 2
+0 COS «1ДГ [— {TtíN — sin TtlN) — — (Tt2N — sin Tt2N)],
TpalN
f2
—a sin c^in cos a2N D2+
22
+ocosaijv[— (1 - cosTtlN) - —(1 - cos tun)],
ш ш
F3 = —a sin aiN cos &2ND3+
+ 0 COS «1ДГ [—9 (1 — COS TtlN)-9 (1 — COS Tt2N)],
ш2 ш2
Fa1N = —a sin a1N cos a2N D4+
А ■ 1 •
+ocosaiAr(—smrtiAr-sinrt2jv) ,
шш
F5a1N = a sin a1N sin a2N D5,
F¡?1N = —a sin a1N sin a2N D6,
F a2N
= —a cos a1N sin a2N D1, F2a2N = —a cos a1N sin a2ND2, F3a2N = —a cos a1N sin a2N D3, F4a2N = —a cos a1N sin a2N D4, F5a2N = a cos a1N cos a2ND5,
У.'.
6
F6a2N = —a cos a1N cos a2N D6,
4 3
Di = —(1 — —Tt 1дг — cosrtiAr) — ш2 8
--о (1 - ~Tt2N ~ COS Tt2N),
ш2 8
D2 = — (4sin7£1дг - Згпдг)-
ш
--(4sinri2Ar - 3tí2at),
ш
2
D3 = —{sm tun ~ Tíín) —
ш2
^7 (sin Tt2N ~ Tt2N), ш2
22
Di =----(1 - COS TtíN) + -(1 - COS Tt2N)?
шш
D5 = Ду(1 - COS TfíN) - Ду(1 - COS Tf2N)]? /, 12 /, 12
1
1
Dß = — sin TtíN-sin Tf íN•
шш
Следует отметить, что определение управления иопт является циклическим процессом с формированием управляющего сигнала в моменты времени tu(to). Тогда реализация метода свободных траекторий с использованием алгоритма релейного и непрерывного управления с аналитическим решением сводится к последовательному выполнению в каждом цикле формирования управления следующих операций.
1. По результату измерения вектора состояния КА формируются начальные условия x(tu). Одновременно, если это первый цикл, задается количество импульсов, с помощью которых предполагается реализовать метод свободных траекторий. Тогда в предположении импульсного характера
тяги вычисляются величина и направление первого и последнего программных импульсов по методике двухимпульсного маневра [1] и формируются начальные условия для вектора Тр, компонентами которого являются Ьм, и углов ориентации а. Отметим, что для промежуточных импульсов, равномерно распределенных на всем интервале заданного времени сближения, начальные значения моментов времени (Ьз, ..., Ьм—2) выбираются из условия: Ьз = Ь4, ..., Ьм—з = Ьм—2 и т.д., то есть для всех промежуточных импульсов их начальная длительность берется равной нулю так же как и углы ориентации двигателя.
Для всех последующих циклов (кроме первого) начальные условия для векторов Тр, а определяются из решения уравнения (3), (4), которые имеют вид
TM = TP"-1) +Д*цир'-1),
V V— 1
aV = aV 1
дЫт-1).
(14)
2. На основе начальных условий векторов Тр, а,х(Ьи) с использованием решения (10) прогнозируется состояние КА на момент времени окончания сближения. По результатам прогноза на основании выражений (13) вычисляются управления. Эти управления иV используются как постоянные для КА на протяжении очередного ^-цикла формирования управления вплоть до получения но-
- ^+1)
вых значений иопт , и предназначены для вычисления векторов Тр и а в соответствии с (14).
3. При достижении очередного момента переключения Ь^х(у,г) в одном из трех каналов управления производится действительное, а не моделируемое включение (выключение) двигателя, и тем самым уменьшается число слагаемых циклограммы (2) и исключается из вычислений соответствующая компонента управления иоптм.
В качестве оценки работоспособности алгоритма рассмотрим задачу управления сближением КА по двухимпульсной схеме при условии, что двигательная установка имеет следующие параметры: а = 0,5 мс—2.
Пусть множество начальных положений КА ограничено следующей областью значений фазовых координат х0 € [2500 м, 10000 м],
Уо,%о € [500 м,4000 м], х0,у0,х0 € [0, — 5 мс—1]. Поставим задачу перевода КА из любой точки этой области в точку с хк = ук = ¿к = хк = Ук = ¿к = 0 за время Ьк — Ьо = 1060 с.
Для обеспечения точности выполнения конечных условий встречи алгоритм с прогнозирующей моделью устанавливает зависимость параметров функционала Рг, Kj от начальных условий сближения. Для этого дополнительно вводится вспомогательный функционал, который по форме соответствует главной части функционала обобщенной работы и, следовательно, в наших условиях имеет вид функции Узад (12). Варьируемыми параметрами вспомогательного функционала явля-
ются коэффициенты основного функционала. Таким образом, в начале решается задача отыскания минимума вспомогательного функционала по параметрам основного функционала, а затем по результатам исследований эмпирически устанавливается математическая зависимость коэффициентов функционала относительно начальных условий сближения. При отыскании минимума вспомогательного функционала могут использоваться различные методы. Так, в [4] применялся метод Дэвида-Флетчера-Пауэлла [5]. Мы использовали простой перебор параметров функционала. С целью сокращения вычислений введем ряд требований.
1. Потребуем одинаковой точности для каждого из трех каналов управления. Тогда для параметров pi можно записать:
Р1 = Р3 = Р5 = Px, Р2 = Р4 = Р6 = Px
На основании условия равного вклада максимальных отклонений [2] примем следующие числовые значения: px = 0,0024 м~2, px = 0,4 м~2 с~2.
2. Потребуем, чтобы коэффициенты усиления функционала K^x,y,z удовлетворяли соотношениям:
ktll = kt21 = ki * ,
k2 * при t < (tk - to)/2
k3 * при t > (tk - to)/2 ,
■ ka11 — ka11 — ka11 — l0 •
kt 12 = kt22
3. Заданные конечные значения фазовых координат Ж1й, Ж3Й, Х5к, входящие в функцию Узад основного функционала, определим ступенчатыми функциями вида
x1k , x3k , x5k —
Д„, если t ^ (tk — to)/2
0, если t > (tk — to)/2
где Дп(п = 1, 3, 5) будем считать варьируемыми коэффициентами.
Таким образом, минимум вспомогательного функционала необходимо искать по шести параметрам: &!*, к2*, кз*, Д1, Дз, Д5.
Дальнейшее упрощение процедуры установления зависимости коэффициентов функционала от начальных условий сближения заключалось в следующем. Для средней точки области начальных условий сближения при условии Д1 = Дз = Д5 = 0 осуществлялся поиск минимума вспомогательного функционала по параметрам к1 *, к2*, кз*. В результате были получены следующие числовые значения: к1* = 0,00147, к2 * = 0,00042, к% = 0,0608, а вычислительная трудоемкость составила 20 итераций. Данные значения коэффициентов усиления были приняты постоянными, и, следовательно, дальнейшие исследования проводились только путем вариации параметров Д1 = Дз = Д5.
С этой целью было выбрано 25 точек, равномерно распределенных по области начальных
условий сближения; для каждой из них определялись значения величин Ді, Д3, Д5, обеспечивающих минимум вспомогательного функционала. При этом максимальная вычислительная трудоемкость определения параметров Ді = Д3 = Д5 для каждой точки не превышала шести итераций.
По результатам проведенных исследований ставилась задача математического описания зависимости величин Д1, Д3, Д5 от начальных условий сближения. Анализ полученных материалов показал следующее. Поскольку для определения начальных значений вектора Тр, необходимых для работы предлагаемого алгоритма, используется методика расчета [1], основанная на предположении импульсного характера действия тяги двигателей, то и параметры функционала целесообразно настраивать не только по вектору начального положения КА, но и по вектору конечного положения, который получается в результате реализации программных импульсов. Такая зависимость была получена эмпирически:
0,069(хіик - ^іх2ик)
при |жіик - 5іХ2ик| ^ 440,
Д1 = < 0,304(хіик 5іх2ик) ,
Є^^Хіик 5іх2ик)103,8
при |хіик - ^іх2ик| > 440,
0,091(х3ик 5іх4ик)
при |х3ик - ^іх4ик| ^ 440,
Д3 = < 0,122(х3ик 5іх4ик) ,
49,7 8^п(х3ик 5іх4ик) при |х3ик - ^іх4ик| > 440,
0,069(х5ик ^2іх6ик)
при |х5ик - ^2х6ик| ^ 440,
0,375(х5ик 52х6ик) ,
1338^п(х5ик 52х6ик)
при|х5 ик - «2хбик| > 440,
0,233(
Д
5
02 і x
1 + x3
02
02)1/2 02
«1
(xf + x32)1/2 < 1500,
при
0,02(x02 + x32 )1/2 при
3,20
(х°2 + хз2)1/2 > 1500,
0,233 |х'01 при |х5| ^ 1500, в2 = <; 0,02|х0|+ 3,20 .
при |х°| > 1500,
x1ик, x3ик, x5ик, x2ик, x4ик, х6ик конечные зна-
чения компонент вектора состояния КА, получаемые при реализации программных импульсов.
Такой подход к настройке коэффициентов функционала, как показали статистические испытания, позволяет:
во-первых, для выбранной области начальных условий получить точностные характеристики по таким показателям, как средние значения компонент вектора конечного состояния КА и средние квадратические отклонения по этим компонентам, которые равны
х1 (Ьк) = —0,0014, х2 (Ьк) = —0,0035, хз(Ьк) = 0,0114, х4 (Ьк )= 0,017, х5(Ьк) = 0, х6(Ьк ) = —0,11,
<7x1 = 0,3, 0x2 =0,15, Охз = 0,45,
Стх4 = 0,16, 7x5 = 0,36, 7x6 = 0,025; во-вторых, выполнять сближение с указанными точностными характеристиками из любой другой области начальных условий, отличной от рассматриваемой в данной работе. При этом начальные условия должны удовлетворять следующему ограничению:
• максимальная длительность работы двигателей не должна превышать 30% времени сближения.
Работа выполнена при финансовой поддержке Правительства Российской Федерации в рамках контракта с Минобрнауки России № 13, G25.31.0028.
Литература
1. Зубов Н.Е. Синтез управления сближением по методу траекторий на основе алгоритма с про-
гнозирующей моделью // Космические исследования. — 1990. — Т. 28, вып. 4. — С. 506.
2. Справочник по теории автоматического управления / под ред. А.А. Красовского. — М.: Наука, 1987.
3. Зубов Н.Е. Управление объектами с релейно-импульсными и непрерывными рулевыми органами на основе алгоритма с прогнозирующей моделью и его приложение в динамике сближения КА // Космические исследования. — 1989. — Т. 27, вып. 2. — С. 206.
4. Буков В.Н., Зубов Н.Е. Релейное управление на основе алгоритма с прогнозирующей моделью // Автоматика и телемеханика. — 1986. — № 6. — С. 36.
5. Химмельблау В. Прикладное нелинейное программирование. — М.: Мир, 1975.
Поступила в редакцию 24.09.2010.
