CC BY
68
возрастанию опасного с точки зрения закрутки интервала. Например, сравнивая А = 5000 кг-м2 и А = 10000 кг-м2, можно отметить, что в первом случае наибольший интервал опасного заклинивания составлял 1860 с, а во втором — 1920 с;
- увеличение скорости выпуска троса приводит к увеличению интервала заклинивания, соответствующего разрыву троса, а также к увеличению интервала, опасного с точки зрения закрутки КА. Для систем с рассматриваемыми параметрами выпуск троса со скоростью, большей 4.3 м/с, приводит к намотке троса на КА независимо от момента заклинивания.
Авторы выражают благодарность профессору В. С. Асланову за постановку задачи и постоянное внимание к работе.
Работа выполнена при финансовой поддержке РФФИ (проект 12-01-00317-а). Библиографический список
1. Cartmell M. P., McKenzie D. J. A review of space tether research // Progress in Aerospace Sciences. 2008. Vol. 44. P. 1-21.
2. Smith H. F. The First and Second Flights of the Small Expendable Deployer System(SEDS) // Proceedings of the Fourth Intern. Conf. on Tethers in Space, Smithsonian Inst., Washington, DC. 1995. P. 43-55.
3. Kruijff M, Heide E. J. van der. Qualification and in-flight demonstration of a European tether deployment system on YES2 // Acta Astronautica. 2009. Vol. 64, № 9-10. P. 882-905.
4. Белецкий В. В., Левин Е. М. Динамика космических тросовых систем. М. : Наука, 1990. 329 с.
5. Zimmermann F, Schottle U. M, Messerschmid E. Optimization of the tether-assisted return mission of a guided re-entry capsule // Aerospace Science and Technology. 2005. № 9. P. 713-721.
6. Асланов В. С., Волошенюк О. Л., Кислов А. В., Ящук А. В. Определение времени выживания космиче-
УДК 629
ОБ ОДНОЙ ЗАДАЧЕ ОПТИМАЛЬНОЙ ПЕРЕОРИЕНТАЦИИ ОРБИТЫ КОСМИЧЕСКОГО АППАРАТА
И. А. Панкратов, Я. Г. Сапунков, Ю. Н. Челноков
Саратовский государственный университет,
E-mail: PankratovIA@info.sgu.ru, ChelnokovYuN@info.sgu.ru
С помощью принципа максимума Понтрягина и кватернионных уравнений решается задача оптимальной переориентации орбиты космического аппарата (КА). Управление (вектор реактивной тяги, ортогональной плоскости орбиты) ограничено по модулю. Функционал, определяющий качество процесса управления, равен взвешенной интегральной сумме времени переориентации орбиты КА и модуля (или квадрата) управления. Сформулированы дифференциальные краевые задачи переориентации орбиты КА. Получены законы оптимального управления, построены условия трансверсальности, не содержащие неопределенных множителей Лагранжа. Приведены примеры численного решения задачи.
Ключевые слова: оптимальное управление, космический аппарат, кватернион.
ской тросовой системы // Изв. Самарского науч. центра РАН. 2010. Т. 12, № 4. С. 138-143.
7. Асланов В. С., Ледков А. С., Стратилатов Н. Р. Пространственное движение космической тросовой системы, предназначенной для доставки груза на Землю // Общероссийский науч.-техн. журн. Полет. 2007. № 2. С. 28-33.
8. Асланов В. С., Ледков А. С., Пироженко А. В., Храмов Д. А. Исследование влияния обрыва тросовой системы на возможность доставки груза на Землю // Сборник трудов XIV Всероссийского научно-технического семинара по управлению движением и навигации летательных аппаратов. Самара : СГАУ, 2011. С. 36-39.
9. Асланов В. С. Влияние упругости орбитальной тросовой системы на колебания спутника // Прикладная математика и механика. 2010. Т. 74, № 4. С. 582-593.
10. Белецкий В. В. Движение спутника относительно центра масс в гравитационном поле. М. : Изд-во Моск. гос. ун-та, 1975. 308 с.
About a Problem of Spacecraft's Orbit Optimal Reorientation
I. A. Pankratov, Ya. G. Sapunkov, Yu. N. Chelnokov
The problem of optimal reorientation of the spacecraft's orbit is solved with the help of the Pontryagin maximum principle and quaternion equations. Control (thrust vector, orthogonal to the orbital plane) is limited in magnitude. Functional, which determines a quality of control process is weighted sum of time and module (or square) of control. We have formulated a differential boundary problems of reorientation of spacecraft's orbit. We have obtained optimal control laws, built the transversality conditions, not containing Lagrange multipliers. Examples of numerical solution of the problem are given.
Keywords: optimal control, spacecraft, quaternion.
1. ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ ОРИЕНТАЦИИ ОРБИТЫ КОСМИЧЕСКОГО АППАРАТА
Будем считать, что вектор ускорения и от тяги реактивного двигателя во все время управляемого движения КА направлен ортогонально плоскости его орбиты. Тогда орбита КА в процессе управления движением центра масс КА не меняет своей формы и своих размеров, а поворачивается в пространстве под действием управления как неизменяемая (недеформируемая) фигура.
Движение центра масс КА будем рассматривать в инерциальной системе координат Х — геоцентрической экваториальной системе координат ОХ1Х2Х3(X) с началом в центре О притяжения Земли. Ось ОХ3 этой системы координат направлена вдоль оси суточного вращения Земли, оси ОХ1 и ОХ2 лежат в плоскости экватора Земли, ось ОХ^ направлена в точку весеннего равноденствия для Земли, ось ОХ2 дополняет систему до правой тройки векторов.
Введем также в рассмотрение систему координат £, связанную с плоскостью и перицентром орбиты КА. Начало этой системы координат находится в центре О, ось £1 направлена вдоль радиуса-вектора перицентра орбиты, ось £3 перпендикулярна плоскости орбиты и имеет направление постоянного по модулю вектора с момента скорости центра масс КА, а ось £2 образует правую тройку с осями £1 и £3. Ориентация системы координат £ в инерциальной системе координат X характеризует собой ориентацию орбиты КА в инерциальном пространстве и задается тремя угловыми оскулирую-щими элементами орбиты: долготой восходящего узла Ои, наклоном орбиты I и угловым расстоянием перицентра от узла .
Дифференциальные уравнения ориентации орбиты КА в инерциальной системе координат в угловых элементах орбиты в рассматриваемом случае имеют вид [1,2]
г . , „ ¿I
r • , n т- di r , ,
, = и- sin(wn + w)cosec i, — = и-cos(wn + w), dt c dt c
dwv r dw
-—= -и- sin(Wn + w)ctgI, -7-dt c dt
(1)
r=
p
1 + e cos w'
c = const,
где ^ — истинная аномалия (угловая переменная, отсчитываемая в плоскости орбиты от ее перицентра и характеризующая положение КА на орбите); г — модуль радиуса-вектора г центра масс КА; р и е — параметр и эксцентриситет орбиты, с — постоянная площадей (модуль вектора момента скорости V = ¿г/^ центра масс КА); и — проекция вектора ускорения и на направление вектора момента скорости центра масс КА (алгебраическая величина реактивного ускорения, перпендикулярного плоскости орбиты КА).
Задача переориентации орбиты КА в угловых переменных формулируется следующим образом: требуется построить управление и, переводящее орбиту КА, изменение ориентации которой описывается уравнениями (1), из заданного начального положения
Ou = Ou (ta)=OU, I = I (to)= I°,
W,
= Wn (t°)= , I0 = 0, П,
в требуемое конечное положение Ou = Ou(t* ) = ou,
I = I (t* )= I * , Wn = Wn (t*)= wn, I * = 0, П.
Решение сформулированной задачи достаточно сложно в силу нелинейности дифференциальных уравнений (1) движения центра масс КА в угловых переменных и наличия в этих уравнениях особых точек I = 0, п. Задача решается гораздо проще, если использовать дифференциальные уравнения ориентации орбиты КА в параметрах Эйлера (Родрига-Гамильтона), имеющие вид [3-5]
dA ^
2— = Л о О£, dt
dw
dt r
= и- (cos w ¿1 + sin w ¿2), c = const,
(2) (3)
1 + e cos w'
где Л = Л° + Л1 ¿1 +Л2¿2 + Л3¿з — кватернион ориентации орбиты КА (кватернионный оскулирующий (медленно изменяющийся) элемент орбиты КА); — отображение вектора О на базис £ (вектор О
r
мгновенной абсолютной угловой скорости орбиты направлен вдоль радиуса-вектора г центра масс КА и определяется формулой: О = (п/с)г); г1, ¿2, ¿3 — векторные мнимые единицы Гамильтона, о — символ кватернионного умножения; Лj, = 0,3) — параметры Эйлера, характеризующие ориентацию орбиты КА (системы координат £) в инерциальной системе координат X.
Параметры Эйлера Лj связаны с угловыми элементами орбиты соотношениями
I
Л0 = сое - соэ 0 2
• I •
Л2 = эт - эт 2 2
+ Шп 2
2
. . 1 /"и —
Л1 = эт - соэ -
1 2 2
I . ( Ои +
Л3 = соэ — эт з2
2
(4)
Отметим, что скалярные уравнения ориентации орбиты в параметрах Эйлера использовались для описания орбитального движения и другими авторами [6, 7].
Отметим также, что если (1) — система 4 нелинейных стационарных дифференциальных уравнений первого порядка относительно угловых переменных Ои, I, , у, то (2), (3) — система 5 нелинейных стационарных дифференциальных уравнений первого порядка относительно параметров Эйлера Лj ^ = 0, 3), и истинной аномалии у. Однако система (2), (3) в отличие от (1) не имеет особых точек, к тому же при переходе в ней от времени Ь к новой независимой переменной у в соответствии с дифференциальным соотношением йу = (с/г2) йЬ мы получаем (при и = п(у)) систему 4 линейных нестационарных дифференциальных уравнений относительно параметров Эйлера Лj.
Отметим также, что системы уравнений (1) и (2), (3) могут рассматриваться как нестационарные системы дифференциальных уравнений третьего и четвертого порядков соответственно относительно переменных Ои, I, и Лj = 0, 3), так как последнее уравнение в этих системах для истинной аномалии у интегрируется в квадратурах независимо от других уравнений, в силу чего переменная у может рассматриваться как известная функция времени Ь. При таком рассмотрении (1) является нелинейной системой, а (2), (3) — линейной.
Указанные обстоятельства делают использование уравнений (2), (3) для решения задачи переориентации орбиты более удобным и эффективным в сравнении с использованием уравнений (1). Такое решение задачи переориентации орбиты КА с использованием комбинированного функционала, равного взвешенной сумме времени переориентации и интегрального квадратичного (в отношении управления) функционала качества, рассмотрено в работах [8-10].
2. ПОСТАНОВКА ЗАДАЧИ ПЕРЕОРИЕНТАЦИИ ОРБИТЫ КОСМИЧЕСКОГО АППАРАТА
Требуется определить ограниченное по модулю управление и:
-Птах < П < Птах < ГО, П = ±|и|,
(5)
ортогональное плоскости орбиты КА, переводящее орбиту КА, движение центра масс которого описывается уравнениями (2), (3), из заданного начального состояния
Ь = Ь0 =0, у(0) = у0, Л(0) = Л0
(6)
в конечное состояние
Ь = Ь*, у(Ь*) = у*, Л(Ь* ) = ±Л* и минимизирующее функционал
(7)
= j (а1 + а2п2) йЬ 0
или = J (а1 + а21и|) йЬ, а1, а2 = сопб1 > 0.
0
При а1 = 1, а2 = 0 имеем задачу переориентации орбиты, оптимальную в смысле быстродействия; при а1 =0, а2 = 1 минимум функционала 12 означает минимум характеристической скорости [1].
г
г
Фигурирующая в краевых условиях кватернионная переменная Л характеризует ориентацию орбиты КА, а переменная ^ — положение КА на орбите. Величины с, р, е, , Л0 и Л* заданы (начальное и конечное значения кватерниона Л могут быть найдены через заданные значения угловых элементов орбиты Ои, I, по формулам (4)). Подлежат определению оптимальный закон управления и = и(Ь) и величины Ь*, .
Четыре компоненты Лj кватерниона Л удовлетворяют условию Л2 + А^Ч—ЬЛ2 + Л2 = 1, поэтому краевое кватернионное условие (7), эквивалентное четырем скалярным, заменим на условие
vect |A(í*) о А* =0,
(8)
эквивалентное трем скалярным (в (8) и далее верхняя волна означает сопряженный кватернион). Такая замена повышает эффективность численного решения задачи оптимальной переориентации орбиты КА.
3. ЗАКОНЫ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ
Поставленную задачу будем решать с помощью принципа максимума [11]. Для этого введем дополнительные переменные M = M0 + M1i1 + M2i2 + M3i3 и х, сопряженные по отношению к фазовым переменным А и Функция Гамильтона-Понтрягина имеет вид
c r
H = + х^г + u— (Ni cos ^ + N2 sin , r2 2c
где в соответствии с принципом максимума постоянная -0о < 0, N1, N2 — компоненты кватерниона N = А о M; для функционала J1 a = a1 + a2u2, для функционала J2 a = a1 + a2|u|, в случае быстродействия a = 1.
В дальнейшем рассматривается случай, когда постоянная -0о < 0. В этом случае в силу однородности функции Гамильтона-Понтрягина относительно сопряженных переменных константу -0о можно выбрать произвольным образом, поэтому в дальнейшем нами полагается -0о = -1. Система уравнений для сопряженных переменных примет вид
dM = M о ^, (9)
2-
dt
dx ^ydr r ,AT . , r2 dr ,AT AT . ,
~dt =2Г dt + u2c Sin ^ - N2 cos - u2c2 dt °OS ^ + N2 sin •
(10)
Отметим, что сопряженное уравнение (9) совпадает по своей форме с фазовым уравнением (2), поскольку кватернионное уравнение (2) обладает свойством самосопряженности.
Законы оптимального управления (т. е. законы управления, удовлетворяющие необходимым условиям оптимальности) находятся из условий максимума функции Н по переменной и с учетом наложенного ограничения (5) и имеют вид:
1. В случае a = 1 (быстродействия)
2. В случае a = a + a2u
u — umax 2
sign k, k = N1 cos ^ + N2 sin
f 1 r k
4a2 c
если
1r
1--|k| < Ur
4a2 c
sign k, если
1r 4a2 c
(11)
(12)
|k| > ur
3. В случае a = a1 + a2|u|
umax sign k, 0,
V u G [ umax jumax]i
1r
если ---|k| — 1 > 0,
2a2 c
1r
если ---|k| — 1 < 0,
2a2 c
1r
если--|k| — 1 = 0.
2«2 c 1 1
(13)
Случай особого управления, когда r|k| = 2a2c, в работе не рассматривается.
uo = ^
u
uo = ^
4. УСЛОВИЯ ТРАНСВЕРСАЛЬНОСТИ
Вводя неопределенные множители Лагранжа Ai, A2, A3, получим условия трансверсальности, соответствующие многообразию конечных состояний (8):
при t = t* M - А* о A = 0, A = A1i1 + A2i2 + A3i3, x = 0. (14)
Из (14) получаем следующие условия трансверсальности, не содержащие неопределенных множителей Лагранжа:
при t = t* Л* Mo + Л* Mi +Л2 M2 + Л* M3 =0, x = 0. (15)
Таким образом, задача оптимальной переориентации орбиты КА сведена к краевой задаче с подвижным правым концом траектории, описываемой системой нелинейных дифференциальных уравнений (2), (3), (9), (10), (11) (или (12), или (13)) десятого порядка и восемью краевыми условиями (6), (8), которые необходимо дополнить двумя условиями трансверсальности (15) и равенством
Ho 1= Ho (Л, M, x, uo) 1=0, имеющим место для оптимального управления uo и оптимальной траектории.
5. УРАВНЕНИЯ В БЕЗРАЗМЕРНЫХ ПЕРЕМЕННЫХ
Для численного решения краевой задачи оптимальной переориентации орбиты КА запишем уравнения этой задачи в безразмерных переменных. Фазовые Л^ и сопряженные Mj переменные являются безразмерными. Безразмерные переменные rb, tb и управление ub связаны с размерными переменными r, t и управлением u соотношениями: r = Rrb, u = umaxub, t = Ttb, где R — характерное расстояние (в его качестве принималась величина, близкая к длине большой полуоси орбиты управляемого КА); V, T — характерные скорость и время соответственно, определяемые соотношениями: V = c/R, T = R2/с.
Отметим, что при переходе к безразмерным переменным в уравнениях для фазовых и сопряженных переменных появляется характерный безразмерный параметр N = umaxR3/с2.
Отметим также, что при численном решении задачи удобно перейти к новой независимой переменной — истинной аномалии. При этом порядок краевой задачи понижается на две единицы.
Таким образом, система фазовых уравнений в безразмерных переменных имеет следующий вид:
= 1N И3 ^Л о (cos шi1 + sin шi2), rb =-1-. (16)
аш 2 4 у 1 + e cos ш
Начальные условия интегрирования этой системы
Ш = шо, Л = Л0 (17)
являются заданными.
Для правого конца траектории КА имеем условия
Ш = ш* =?, Л = ±Л* ^ vect [Л о Л*J = 0, (18)
где Л* — заданная кватернионная величина.
Ограничение по управлению в безразмерном виде |ub| < 1.
Система сопряженных уравнений в безразмерных переменных имеет вид
dM = 1N (rb)3 ubM о (cos шi1 + sin шi2). (19)
аш 2 к ;
Таким образом, в безразмерных переменных задача оптимальной переориентации орбиты КА сведена к краевой задаче с подвижным правым концом траектории, описываемой системой нелинейных дифференциальных уравнений (16), (19) восьмого порядка, и семью краевыми условиями (17), (18), которые необходимо дополнить первым из условий трансверсальности (15) и равенством гамильтониана нулю в конце движения. При этом законы оптимального управления аналогичны (11) (или (12), или (13)). (В дальнейшем верхние индексы у безразмерных переменных опускаются.)
6. ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ ПЕРЕОРИЕНТАЦИИ ОРБИТЫ КОСМИЧЕСКОГО АППАРАТА
Задачи оптимального управления решались для функционалов J1 и J2. Алгоритмы численного решения задач реализуют комбинацию метода Рунге-Кутта 4-го порядка точности и двух методов решения краевых задач: модифицированного метода Ньютона и метода градиентного спуска.
Величины, характеризующие форму, размеры орбиты КА, начальное и конечное положения КА на орбите, начальную и конечную ориентации орбиты КА, полагались равными [4] (aor — большая полуось орбиты):
aor = 25500000 м, umax = 0.101907 м/с2, N = 0.35;
для начального положения КА (<0 = 3.940323 рад.) Л0 = 0.679417, Л? = -0.245862, Л0 = -0.539909, Лз = -0.353860; для конечного положения КА вариант 1 (малое отличие в ориентациях орбит КА): Л0 = 0.678275, Л1 = -0.268667, Л* = -0.577802, Л* = -0.366116; вариант 2 (большое отличие в ориентациях орбит КА): Л0 = -0.440542, Л* = -0.522476, Л* = -0.125336, Л* = -0.719189.
Значения выбранных масштабирующих множителей равны: Я = 26000000 м, V = 2751.405874 м/с, Т = 9449.714506 с. Указанные значения этих величин отвечают значениям декартовых координат и проекций вектора скорости центра масс КА, приведенным в [12].
Ориентации начальной и конечной орбит КА характеризуются параметрами Эйлера Л0 и Л*, 3 = 0,3. Если в варианте 1 эти значения близки (отличие ориентаций орбит по долготе восходящего узла, наклону, угловому расстоянию перицентра от узла составляет единицы градусов: ДОи = ^и - Пи = -3.30°, Д! = I0 - I* = -1.51°, Дшп = шП - < = -1.59°), то в варианте 2 они существенно отличаются (отличие ориентаций орбит в угловой мере составляет десятки градусов: ДПи = -32.00°, ДI = -117.57°, Дшп = 39.96°).
На рис. 1-4 приведены законы изменения фазовых, сопряженных переменных, управления и фазовые портреты для случая быстродействия (е = 0), случая минимизации характеристической скорости (е = 0.25) для варианта 2 и случая минимизации функционала (е = 0.5, а1 = 1.0, а2 = 4.2) для вариантов 1 и 2.
10 12 14 t
10 12 14 t
No
-7 -6.5 -6 -5.5 -5 -4.5 -4 М
U 1 0.5 0
-0.5 -1
0
10 12 14 t
Рис. 1. Круговая орбита, вариант 2, быстродействие
б
а
в
г
Aji
А0
—
10 12 t
-11 -10 -9 -8 -7 -6 -5 -4 Ni
M-J10
5 0 -5 -10
мх
мъ
ч М0
____ —
U 1 0.5 0
-0.5 -1
2 4 6
б
2 4 6
г
t*
Рис. 2. Эллиптическая орбита (e = 0.25), вариант 2, / |u| dt ^ min
0
10 12 t
10 12 t
л,.,
0.5 0
-0.5 -1
A,
A-
012345678 9t
N-, 40
-5 0 5 10 15 20 25 30 35 Nx
Mj 200 150 100 50 0 -50 -100 -150 -200
Mn
M,
Mo
M,
012345678 9t
б
012345678 t
Рис. 3. Эллиптическая орбита (e = 0.5), вариант 1, /(1 + 4.2u2) dt ^ min
0
а
в
а
в
г
0123456789 t
Mj Mj__
200
100 ■
0 Му
-100 мъ /Ч»
-200 ^____
0 1
-200 -150 -100 -50 0 50 Nx
U 1 0.5 О
-0.5 -1
4 5
б
9 t
0 1 2 3 4 5 6 7
9 t
Рис. 4. Эллиптическая орбита (e = 0.5), вариант 2, /(1 + 4.2u2) dt ^ min
о
В случае быстродействия управление носит релейный характер, что доставляет большие неудобства при численном решении задачи. Длительность процесса переориентации орбиты КА составила 14.209544 безразмерных единиц, или 134276.134060 с (37.298926 ч). Имеется 5 переключений знака управления на всем интервале движения КА. Из графиков (см. рис. 1) видно, что в момент переключения управления происходит резкая смена законов изменения фазовых и сопряженных переменных.
Длительность процесса переориентации орбиты КА (e = 0.25) в случае минимизации характеристической скорости составила 13.723368 безразмерных единиц, или 129681.91 с (36.022753 ч). Значение минимизируемого функционала при этом оказалось равным 8.081946.
Длительность процесса переориентации орбиты КА (e = 0.5) в случае минимизации функционала Ji составила 8.553301 безразмерных единиц, или 80826.252534 с (22.451737ч) для варианта 1 и 9.317650 безразмерных единиц, или 88049.132367 с (24.458092 ч) для варианта 2. Значение минимизируемого функционала при этом оказалось равным 10.817037 и 42.600058 соответственно.
Отметим, что при смене знака управления с +1 на -1, и наоборот, резко меняется закон изменения переменной N2(N1) на фазовых портретах.
Отметим также и выявленную неединственность численного решения краевых задач оптимальной переориентации орбиты КА, связанную с нелинейностью дифференциальных уравнений задачи. При одних и тех же граничных условиях в постановке краевой задачи оптимального управления (с одним и тем же функционалом) получены различные решения для законов движения, управления и поведения сопряженных переменных. Из них было выбрано то, при котором функционал минимален.
В ходе численного решения задачи были выявлены особенности и закономерности построенных оптимальных траекторий и управлений. Так, установлено, что при увеличении эксцентриситета орбиты от 0 до 0.5 уменьшается количество активных участков движения, их длительности и общее время переориентации. При увеличении эксцентриситета орбиты от 0.5 до 0.7 и далее пропадают участки, на которых управление принимает свое максимальное по модулю значение в случае минимизации затрат времени и энергии. В случае малого различия начальной и конечной ориентаций орбит (единицы градусов) время переориентации, а также диапазоны изменения фазовых и сопряженных переменных меньше, чем в случае, когда различие начальной и конечной ориентаций орбит составляет десятки градусов.
Работа выполнена при финансовой поддержке РФФИ (проект 12-01-00165).
а
в
г
t
Библиографический список
1. Абалакин В. К., Аксенов Е. П., Гребенников Е. А., Демин В. Г., Рябов Ю. А. Справочное руководство по небесной механике и астродинамике. М. : Наука, 1976. 864 с.
2. Дубошин Г. Н. Небесная механика. Основные задачи и методы. М. : Наука, 1968. 799 с.
3. Челноков Ю. Н. Применение кватернионов в теории орбитального движения спутника. I, II // Космические исследования. 1992. Т. 30, вып 6. С. 759-770; 1993. Т. 31, вып. 3. C. 3-15.
4. Челноков Ю. Н. Применение кватернионов в задачах оптимального управления движением центра масс космического аппарата в ньютоновском гравитационном поле. I-III // Космические исследования. 2001. Т. 39, вып 5. С. 502-517; 2003. Т. 41, вып. 1. С. 92-107; 2003. Т. 41, вып. 5. С. 488-505.
5. Челноков Ю. Н. Кватернионные и бикватернионные модели и методы механики твердого тела и их приложения. М. : Физматлит, 2006. 512 с.
6. Deprit A. Ideal frames for perturbed keplerian motions // Celestial Mechanics. 1976. Vol. 13, № 2. P. 253-262.
7. Брумберг В. А. Аналитические алгоритмы небесной механики. М. : Наука, 1980. 208 с.
8. Ненахов С. В., Челноков Ю. Н. Кватернионное решение задачи оптимального управления ориентацией орбиты космического аппарата // Бортовые интегрированные комплексы и современные проблемы управления : сб. тр. междунар. конф. М. : МАИ, 1997. С. 5960.
9. Сергеев Д. А., Челноков Ю. Н. Оптимальное управление ориентацией орбиты космического аппарата // Проблемы точной механики и управления: сб. науч. тр. ИПТМУ РАН. / Сарат. гос. техн. ун-т. Саратов, 2002. С. 64-75.
10. Афанасьева Ю. В., Челноков Ю. Н. Оптимальное управление ориентацией орбиты космического аппарата // Математика. Механика : сб. науч. тр. Саратов : Изд-во Сарат. ун-та, 2005. Вып. 7. С 153-155.
11. Понтрягин Л. С., Болтянский В. Г., Гамкрелид-зе Р. В., Мищенко Е. Ф. Математическая теория оптимальных процессов. М. : Наука, 1983. 393 с.
12. Бордовицына Т. В. Современные численные методы в задачах небесной механики. М. : Наука, 1984. 136 с.
