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

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

В этом случае траекторное управление КА будет осуществляться за счет ориентации вектора тяги в выбранной системе координат, и, соответственно, этот вектор через углы ориентации α 1 и α 2 запишется так:

-→ a

a cos α1 cos α2 a sin α1

a cos α1 sin α2

где a — модуль постоянного ускорения двигателя.

Уравнения относительного движения КА для круговой орбиты в орбитальной системе координат будут иметь вид [1]:

x˙ = Ax + Ba,

где x =[x1 ,x2, ..., x6]Т,

A =

00 02 ω 00

00 0100

0 -2ω 3ω2

00   0  0-ω2 0

B =

x 1 = x , x 2 = x ˙, x 3 = y , x 4 = y ˙ ,x 5 = z , x 6 = z ˙— фазовые координаты активного КА,

ω = const — орбитальная угловая скорость пассивного КА.

Аналогично, как и в [1], предположим, что моменты времени начала маневра сближения t0 и его окончания tk заданы, причем tk - t0 фиксировано. Управление движением ведется по методу свободных траекторий [1] и состоит в определении первого программного импульса в момент t0 , предназначенного для обеспечения встречи активного и пассивного КА за заданный промежуток времени tk - t0 при их дальнейшем свободном движении. Второй программный импульс в конце маневра сближения выравнивает скорости КА. Возможны дополнительные программные импульсы в промежуточных точках траектории. Управление носит импульсный характер, что позволяет находить величину корректирующих импульсов из уравнения (2) при a =0и соответственно углы α1, α2 ориентации импульса. Но поскольку коррекции реализуются только с помощью ограниченной по величине тяги двигателя, точное выполнение условий встречи таким образом невозможно. Аналогично [1] предлагается рассматривать расчет поправок на конечную тягу двигателя как задачу поиска оптимального управления с использованием алгоритма с прогнозирующей моделью [2]. Однако в отличие от [1], где алгоритм управления с прогнозирующей моделью представляет собой чисто релейное управление, в данном случае он будет сочетать в себе совместно как релейное, так и непрерывное управление. Особенности такого алгоритма совместного непрерывного и релейного управлений с использованием прогнозирующей модели рассмотрены в [3]. Здесь в отличие от [3] в силу наличия аналитического решения системы (2) получим редакцию алгоритма управления объектом с непрерывными и релейными рулевыми органами в редакции с аналитическим решением. В качестве релейных исполнительных органов будем рассматривать двигатель постоянной тяги, положение которого (в смысле «включено», «выключено») характеризует вектор θp, компоненты ко- торого выражаются следующим образом:

M

9p(t, ti, ..., tM) = -a ^(—1)M1M(t — t^),    (3)

μ =1

где t — текущее время, 1 _M ( t — t _M ) — единичная ступенчатая функция, соответствующая μ -му по порядку изменения компоненты, θ _p , t μ — моменты времени начала или конца импульса, M — число скачкообразных изменений функции (3). Если N — число импульсов в методе свободных траекторий, то M = 2 N .

В качестве компонент управления релейных исполнительных органов рассмотрим моменты времени t μ функции (3). Для них можно записать

t^ = u р.

Роль непрерывных рулевых органов будут выполнять углы ориентации вектора тяги α 1 N и α 2 N , и, соответственно, уравнение управления для них выглядит так:

a = u н.

Здесь a = [ a 1 n ,a 2 n ] T .

Объединим компоненты t μ вектора релейного управления θ _p и непрерывного α в единый вектор управления y = [ t _M ,a ] T и запишем объект управления (2) в виде

x = f(xyt), y = u.            (4)

Здесь функции V _зад , Q _зад выражаются квадратичными формами:

Vзад[ x (tk ) ,tk ] = 0,5 xT (tk ) px (tk ) ,

Q зад[ x ( T ) ,T ] =0,5 xT ( T ) ex ( T ) ,

ρ , β , K ^{- 1} — положительно заданные квадратные матрицы параметров.

Выражение вектора тяги (1) с учетом (3) будет иметь вид

N 2

52—a 52(—1)M1 m(t—

N=1

N 2

^ ^— a ^^{( —} 1) ^{M 1} M ⁽ t ^—

N=1

N 2

^ ^— a ^^{( —} 1) ^{M 1} M ⁽ t ^—

N=1

t^N) cos a 1 n cos a2n

I mn )sin a 1 n

tMN) cos a 1N sin a2n

Следовательно, система (2) на основании (3) и (8)

запишется так:

x; 1 rc 2 re 3 x; 4 re 5 re 6

x 1 x 2 x 3 x 4

x 5 x 6

+

Для минимизируемого функционала, заданного в виде [2]:

+B

I = Vзад

t k

[ x (tk) ,tk ] + I

Qзад(X(T) ,T) dT+

t

tk

+0,5 I [uT(t)K_ 1 u(t) + uTnm(t)K_ 1 uопт(T)]dT, (5) t при наличии аналитического решения первого уравнения системы (4) вида

X (x,y,t,tk)

оптимальное управление u определится выражением

u (t)   uопт(t)  

∂            T∂            T

— K^ idyX(x,y,t0,tk)     d^Vзад(x,y,tk)| +

+

tk

∂ ^T

J [^X(x,y,t,T)

∂ ^T

d^Q зад( x,y,T)   dT}. (6)

t

Соответственно при Q _зад = 0

u (t) = u оп( t) =

∂ ^T

k^ dyx(x,y,t0,tk)

∂ ^T

dXV3ag(x,y,tk ) }.

N 2

^L —a 52 (—1)M1 m(t — tMN)cos a 1 n cos a2n

N=1

N 2

52 — a 52( — 1) ^M 1 m ( t — t MN )sin a 1 n

N=1

N 2

У2 — a 52 ( — 1) ^M 1 m ( t — t MN )cos a 1 n sin a 2 n

N=1

t 1 N t 2 N α 1 N α 2 N

u t 1 N u t 2 N u α 1 N u α 2 N

При N = 2 размерность вектора u , входящего в (9) равна 8, при N = 3 — соответственно 12 и так далее.

Решение первого уравнения системы (9) при [ 1 1 n ,t ₂ n ,<а 1 n ,<а ₂ n ] T = 0 и произвольных начальных условиях x ( t _u ), 6 _p ( t _u ,t _MN ) имеет вид

	x 1 ( t k ) x з ( t k ) x 5 ( t k ) x 2 ( t k ) x 4 ( t k ) x 6 ( t k )	=
1   6( t — sin t )      0       ; 1 (4sin t — 3 t ) ^ (1 — cos t )     0 0 4 — 3cos t      0      — 2 (1 — cos t )      1 sin t       0 0        0          cos t           0              0         1 sin t 0 6 w (1 — cos t )     0       4cos t — 3       2sin t       0 X 0    3 w sin t        0         — 2sin t         cos t        0 0       0        —w sin t        0             0        cos t

	x 1 ( t u )		Ф 1 q 11)	Ф 1 q 12)	0
	x 3 ( t u )	N	Ф 1(31)	Ф 1(32)	0
X	x 5 ( t u ) X 2 ( t u )	+ E N =1	0 Ф 1(21)	0 Ф 1(22)	Ф 1(53) 0
	x 4 ( t u )		Ф 1(41)	Ф 1(42)	0
	X 6 ( t u )		0	0	Ф 63

X

a 1(t — t1 n ) cos аi n cos а2n a 1(t — t1 N) sin а 1 N

a 1(t — t1 N) cos а 1 N sin а2N

	Ф 2 q 11)	Ф 2 q 12)	0
N	Ф 2(31)	Ф 2(32)	0
E N =1	0 Ф 2(21)	0 Ф 2(22)	Ф 2(53) 0	X
	Ф 2(41)	Ф 2(42)	0
	0	0	Ф 2(63)

функционала, которую определим квадратичной функцией вида

1 ⁶                       ₂

V зад ^{[ X (} t k ) ,t k ^]    2 ^ ' P i ^[ x i ⁽ t k ) x ik ^] ,      ⁽¹²⁾

i =1

где ρ i — весовые коэффициенты, x ik — заданные конечные значения фазовых координат КА. Функцию Q _зад положим равной нулю.

В соответствии с (7) при функции V _зад , определенном выражением (12), с использованием записи аналитического решения (11) выражение для оптимального управления будет иметь вид

X

a 1(t — t2N) cos а 1 N cos а2N a 1(t — t2N) sin а 1 N

a 1(t — t2N) cos а 1 N sin а2N

где т = ш ( t k — t u ), т i n   = ш ( t k — t 1 n ),

u опт ( t u ) = — K 9 X T ( t k ,t u		,t 1 N ,t 2 N ∂y	,а 1 n ,а 2 n ) x
	Р 1 [ X 1 ( t k ) Р 2 [ X 2 ( t k ) Р 3 [ X 3 ( t k )	— X 1 k ] — X 2 k ] — X 3 k ]
X	Р 4 [ X 4 ( t k ) Р 5 [ X 5 ( t k ) Р 6 [ X 6 ( t k )	— X 4 k ] X 5 k ] X 6 k ]

После вычисления частных производных в развернутом виде получим

т2 N = Ш (tk — t2 N ),

4 _                  3 Q ,

Ф1,2q 11) = Щ2 (1 — cos т1,2N — ^ т1,2N ) ,

^Ф 1 , 2 q 12) = — ^{( т} 1 , 2 N ^{— sin т} 1 , 2 N ) ,

Ф1,2(21) = w (4 sin т1,2N — 3т1,2N) ,

Ф1,2(22) = Ш (1 — cos т1,2N) ,

Ф1,2(31) = — (sin т1,2N — т1,2N) ,

Ф1,2(32) = Ш2 (1 — cos т1,2N),

Ф1,2(41) = Ш (cos т1,2N — 1),

^Ф 1 , 2(42) = Ш Sin ^т 1 , 2 N ,

Ф1,2(53) = — (1 — cos т1,2N),

N E N=1 ut1N ut 2 N uα1N uα2N N = — E N=1 KX X " Ft1N Ft2N F а 1N F α2N F2t1N Ft2N Fα1N α2N F2 Ft1N   Ft1N F3t2N   Ft2N Fα1N  Fα1N α2N    α2N F3     F4 t1N F5 t2N F5 α1N F5 α2N F5 t1N F6 t2N F6 α1N F6 α2N F6 X X Р 1[X 1( tk) — X1 k ] Р 2[ X 2( tk ) — X 2 k ] Р 3[ X 3( tk) — X 3 k ] Р4 [X4 (tk) — X4k ] Р5 [X5 (tk) — X5k ] Р6 [X6 (tk) — X6k ] , (13) где K — квадратная диагональная матрица,

K = diag( kt 1 n ,kt 2 n ,ka 1 n ,ka 2 n ),

Ft1,2N = _ cos а 1 n cos а2n (3тt 1,2n — 4 sin тt 1,2n )+ ω

^Ф 1 , 2(63) = Ш Sin ^т 1 , 2 N .

В компактном виде аналитическое решение (10) запишем так:

+ w sin а 1N (cos Tt 1,2N — 1), t1,2N

F 2    = a cos а 1 n cos а 2 n (3 — 4cos тt 1,2 n ) —

X [ X (tu ) ,tk ,tu,t 1N ,t 2 N ,а 1N ,а 2 N ] =

N

Xод[ X ( tu ) ,tk ,tu] + ^ ' X ч ( tk ^ut 1N ,t 2 N,а 1N,а 2 N ).

N =1

Зададим функционал (5) в соответствии с требо-

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

точку встречи «отвечает» терминальная часть

—2a sin а 1 N sin тt 1,2N,

1 1 , 2 N    ^{2 a}

F3 ,   =—cos а 1 n cos а2n (1 — cos тt 1,2n ) —

ω

— a sin тt 1,2n , ω t1,2N

F4    =2a cos а 1 n cos а2n sin тt 1,2n —

—a sin а 1 N cos тt 1,2N, t1,2N       a

F5     =--cos а 1 n sin а2n sin тt 1,2n,

ω t1,2N

F6 ,   = —a cos а 1 N sin а2N cos тt 1,2N,

Fа1N = —a sin а 1 N cos а2ND 1+

+ a cos α 1 N [ ² ( τ t 1 N - sin τ t 1 N ) - ² ( τ t 2 N - sin τ t 2 N )] , ω ²                  ω ²

F ₂ ^{α 1 N} = - a sin α 1 N cos α 2 N D 2 +

+ a cos α 1 N [— (1 - cos τ t 1 N ) - —(1 - cos τ t 1 N )] , ωω

F ₃ ^{α 1 N} = - a sin α 1 N cos α 2 N D 3 +

+ a cos α 1 N [ ¹ ₂ (1 - cos τ t 1 N ) - ¹ ₂ (1 - cos τ t 2 N )] ,

F ₄ ^{α 1 N} = - a sin α 1 N cos α 2 N D 4 +

+ a cos α 1 N ( — sin τ t 1 N

— sin τ t 2 N )] , ω

F ₅ ^{α 1 N} = a sin α 1 N sin α 2 N D 5 , F ₆ ^{α 1 N} = - a sin α 1 N sin α 2 N D 6 , F ₁ ^{α 2 N} = - a cos α 1 N sin α 2 N D 1 , F ₂ ^{α 2 N} = - a cos α 1 N sin α 2 N D 2 , F ₃ ^{α 2 N} = - a cos α 1 N sin α 2 N D 3 , F ₄ ^{α 2 N} = - a cos α 1 N sin α 2 N D 4 , F ₅ ^{α 2 N} = a cos α 1 N cos α 2 N D 5 ,

F ₆ ^{α 2 N} = - a cos α 1 N cos α 2 N D 6 ,

D 1 = ⁴ (1 ω ²

— τt 1 N

- cos τ _t 1 N ) -

- _{ω 2} (1 - ₈ τ _t 2 N - cos τ _t 2 N ) , D 2 = — (4 sin τ t 1 N - 3 τ t 1 N ) -

—(4 sin τ t 2 N - 3 τ t 2 N ) , ω

D 3 = ² (sin τ t 1 N - τ t 1 N ) - ω ²

- 2 (sin τt2N - τt2N), ω2

D 4 = - — (1 - cos τ t 1 N ) + ^^^^^ (1 - cos τ t 2 N ) ,

D 5 = _ω ¹ ₂ (1 - cos τ t 1 N ) - _ω ¹ ₂ (1 - cos τ t 2 N )] ,

D 6 = ^^^^^™ sin τ t 1 N - ^^^^^ sin τ t 1 N .

Следует отметить, что определение управления u _опт является циклическим процессом с формированием управляющего сигнала в моменты времени t _u ( t 0 ). Тогда реализация метода свободных траекторий с использованием алгоритма релейного и непрерывного управления с аналитическим решением сводится к последовательному выполнению в каждом цикле формирования управления следующих операций.

• 1.    По результату измерения вектора состояния КА формируются начальные условия x ( t _u ). Одновременно, если это первый цикл, задается количество импульсов, с помощью которых предполагается реализовать метод свободных траекторий. Тогда в предположении импульсного характера

• 2.    На основе начальных условий векторов T _p , α,x ( t _u ) с использованием решения (10) прогнозируется состояние КА на момент времени окончания сближения. По результатам прогноза на основании выражений (13) вычисляются управления. Эти управления u _o ^ν используются как постоянные для КА на протяжении очередного ν -цикла формирования управления вплоть до получения но- ( ν +1)

• 3.    При достижении очередного момента переключения t _μx ( _y,z ) в одном из трех каналов управления производится действительное, а не моделируемое включение (выключение) двигателя, и тем самым уменьшается число слагаемых циклограммы (2) и исключается из вычислений соответствующая компонента управления u _опт μ .

тяги вычисляются величина и направление первого и последнего программных импульсов по методике двухимпульсного маневра [1] и формируются начальные условия для вектора T _p , компонентами которого являются t μ , и углов ориентации α . Отметим, что для промежуточных импульсов, равномерно распределенных на всем интервале заданного времени сближения, начальные значения моментов времени ( t 3 , ..., t M _- 2 ) выбираются из условия: t 3 = t 4 , ..., t M _- 3 = t M _- 2 и т.д., то есть для всех промежуточных импульсов их начальная длительность берется равной нулю так же как и углы ориентации двигателя.

Для всех последующих циклов (кроме первого) начальные условия для векторов T _p , α определяются из решения уравнения (3), (4), которые имеют вид

T ( _pν ) = T ( _pν - 1) + Δ t _ц u ( _pν - 1) ,

α ^ν = α ^{ν - 1} + Δ t ц u ⁽ _н ^{ν - 1)} .          (14)

вых значений u опт , и предназначены для вычисления векторов T p и α в соответствии с (14).

В качестве оценки работоспособности алгоритма рассмотрим задачу управления сближением КА по двухимпульсной схеме при условии, что двигательная установка имеет следующие параметры: a =0 , 5 мс ^{- 2} .

Пусть множество начальных положений КА ограничено следующей областью значений фазовых координат xo ∈   [2500 м,10000 м], yo,zo ∈ [500 м,4000 м], x˙ o,y˙o,z˙o ∈ [0, - 5 мс-1]. Поставим задачу перевода КА из любой точки этой области в точку с xk = yk = zk =x˙k =y˙k =z˙k =0 за время tk - t0 = 1060 с.

Для обеспечения точности выполнения конечных условий встречи алгоритм с прогнозирующей моделью устанавливает зависимость параметров функционала ρi, Kj от начальных условий сближения. Для этого дополнительно вводится вспомогательный функционал, который по форме соответствует главной части функционала обобщенной работы и, следовательно, в наших условиях имеет вид функции Vзад (12). Варьируемыми параметрами вспомогательного функционала явля- ются коэффициенты основного функционала. Таким образом, в начале решается задача отыскания минимума вспомогательного функционала по параметрам основного функционала, а затем по результатам исследований эмпирически устанавливается математическая зависимость коэффициентов функционала относительно начальных условий сближения. При отыскании минимума вспомогательного функционала могут использоваться различные методы. Так, в [4] применялся метод Дэвида–Флетчера–Пауэлла [5]. Мы использовали простой перебор параметров функционала. С целью сокращения вычислений введем ряд требований.

• 1.    Потребуем одинаковой точности для каждого из трех каналов управления. Тогда для параметров p i можно записать:

• 2.    Потребуем, чтобы коэффициенты усиления функционала K μ_x,y_,z удовлетворяли соотношениям:

p 1 = p з = p 5 = Px,  p 2 = p 4 = p 6 = px

На основании условия равного вклада максимальных отклонений [2] примем следующие числовые значения: p x = 0 , 0024 м ^{_ 2} , p x = 0 , 4 м ^{_ 2} с ^{_ 2} .

kt 11 = kt21 = k 1 * , условий сближения; для каждой из них определялись значения величин Δ1 , Δ3, Δ5, обеспечивающих минимум вспомогательного функционала. При этом максимальная вычислительная трудоемкость определения параметров Δ1 =Δ3 =Δ5 для каждой точки не превышала шести итераций.

По результатам проведенных исследований ставилась задача математического описания зависимости величин Δ 1 , Δ 3 , Δ 5 от начальных условий сближения. Анализ полученных материалов показал следующее. Поскольку для определения начальных значений вектора T p , необходимых для работы предлагаемого алгоритма, используется методика расчета [l], основанная на предположении импульсного характера действия тяги двигателей, то и параметры функционала целесообразно настраивать не только по вектору начального положения КА, но и по вектору конечного положения, который получается в результате реализации программных импульсов. Такая зависимость была получена эмпирически:

Δ 1

0 , 069( x 1ик — s 1 x 2ик )

при | x 1ик — s 1 x 2ик I <  440 ,

<  0,304( x1ик — S 1x 2ик) —

— sign( x 1ик - 1 x 2ик)103,8

при | x 1ик — s 1 x 2ик I >  440 ,

= ,   = / k2 *  при  t < (tk — to) / 2

12- t 22- | k 3 *  при   t> ( t k — t o ) / 2

^k a 11 ^k a 11    ^k a 11    ^k a 11    ^k a 11    ¹⁰

Δ 3

3. Заданные конечные значения фазовых координат x 1 k , x 3 k , x 5 k , входящие в функцию V зад основного функционала, определим ступенчатыми функциями вида

0 , 091( x 3ик — s 1 x 4ик )

при | x 3ик ^- 1 x 4ик I <  440 ,

* ⁰ , 1 ^{22( x} 3ик s 1 ^x 4ик )

—49,7 sign( x 3ик   s 1x 4ик)

при | x 3ик ^- 1 x 4ик I >  440 ,

Δ 5

^x 1 k , ^x 3 k , ^x 5 k

Ап,

0,

если t < (tk — to) / /2 если t> (tk — to) / 2   ,

0 , 069( x 5ик — s 21 x бик )

при | x 5ик — s 2 x бик I <  440 , <   ⁰ , 3 ^{75( x} 5ик s 2 ^x бик )

— 133sign( x 5ик — s 2 x бик )

при I x 5ик ^- 2 x бик I >  440 ,

где А n ( n = 1 , 3 , 5) будем считать варьируемыми коэффициентами.

Таким образом, минимум вспомогательного функционала необходимо искать по шести параметрам: k 1 * , k ₂ * , k 3 * , А 1 , А ₃ , А 5 .

Дальнейшее упрощение процедуры установления зависимости коэффициентов функционала от начальных условий сближения заключалось в следующем. Для средней точки области начальных условий сближения при условии Δ 1 =Δ 3 =Δ 5 = 0 осуществлялся поиск минимума вспомогательного функционала по параметрам k 1 * , k ₂ * , k 3 * . В результате были получены следующие числовые значения: k 1 * = 0 , 00147, k ₂ * = 0 , 00042, k 3 = 0 , 0608, а вычислительная трудоемкость составила 20 итераций. Данные значения коэффициентов усиления были приняты постоянными, и, следовательно, дальнейшие исследования проводились только путем вариации параметров Δ 1 =Δ 3 =Δ 5 .

С этой целью было выбрано 25 точек, равномерно распределенных по области начальных s1 =

0,233( x02 + x32)1//

при (x 12 + x32 )1 /2 < 1500, 0,02( x 12 + x 32)1 / 2 + 3,20

при ( x 12 + x 32 ) 1 / ² >  1500 ,

( 0,233Ix51 при Ix01 < 1500,

s 2 = < 0,02 Ix 51 +3,20 .

[ при Ix01 > 1500, x 1ик, x3ик, x5ик, x2ик, x4ик, xбик — конечные значения компонент вектора состояния КА, получаемые при реализации программных импульсов.

Такой подход к настройке коэффициентов функционала, как показали статистические испытания, позволяет:

во-первых, для выбранной области начальных условий получить точностные характеристики по таким показателям, как средние значения компонент вектора конечного состояния КА и средние квадратические отклонения по этим компонентам, которые равны

x 1( tk ) = — 0,0014, x 2 (tk ) = — 0,0035, x s( tk )=0,0114,

x 4 (tk ) = 0,017, x 5( tk ) = 0, x б (tk ) = — 0,11,

ax i = 0,3, ax 2 = 0,15, ax 3 = 0,45, ax 4 = 0,16, ax 5 = 0,36, ax6 = 0,025;

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

• •    максимальная длительность работы двигателей не должна превышать 30% времени сближения.

Работа выполнена при финансовой поддержке Правительства Российской Федерации в рамках контракта с Минобрнауки России № 13, G25.31.0028.