CC BY
17УДК 531.01:629.191
БОТ 10.18698/2308-6033-2017-9-1667
Исследование оптимального трехимпульсного перехода на высокую орбиту искусственного спутника Луны
© Е.С. Гордиенко
АО «НПО Лавочкина», г. Химки, Московская обл., 141400, Россия МГТУ им. Н.Э. Баумана, Москва, 105005, Россия
Исследовано оптимальное выведение космического аппарата от Земли на высокую круговую полярную орбиту искусственного спутника Луны радиусом 6 тыс. км при использовании трехимпульсной биэллиптической схемы, так называемого штерн-фельдовского перехода. Задача выполнена с учетом возмущений от нецентральности поля Луны, гравитационных полей Земли и Солнца, а также конечности тяги двигателя. Оптимальная траектория определена варьированием как параметров импульсов, так и точек их приложения. Анализ проведен в два этапа. На I этапе рассмотрен случай импульсной тяги, при этом ориентация вектора тяги сначала, в «плоском» случае, задается в оскулирующей плоскости углами тангажа или атаки, затем, в «пространственном», — добавляется угол рыскания. Показано, что с увеличением максимального расстояния оптимальные моменты приложения импульсов смещаются из апсидальных точек орбит, а второй импульс не сона-правлен с вектором текущей скорости. В «плоском» варианте ориентации импульса получены характеристики, близкие к апсидальному случаю импульсной тяги, в «пространственном» — конечная масса космического аппарата несколько увеличивается. На II этапе рассмотрен случай конечной тяги, при этом получены результаты, близкие к импульсному случаю.
Ключевые слова: космический аппарат, оптимизация траекторий, искусственный спутник Луны, трехимпульсная биэллиптическая схема, штернфельдовский переход, квазиньютоновский метод Пауэлла.
Введение. Оптимальное выведение космического аппарата (КА) от Земли на высокую орбиту искусственного спутника Луны (ИСЛ) — важная научная задача. В обычном случае перелета КА на низкую орбиту ИСЛ оптимальным признано одноимпульсное торможение. Полученные ранее результаты анализа трехимпульсного перехода при сообщении импульсов в апсидальных точках [1, 2] (далее по тексту — апсидальное, или апсидальное импульсное решение) показали преимущество такого перехода над одноимпульсным. В качестве опорного был выбран вариант трехимпульсного перехода с характе-
ристиками:
Расстояние до центра Луны в точке Р1-г1, км ..................................... 1838,57
Расстояние до центра Луны в точке Р2-г2, тыс. км ............................. 45
Расстояние до центра Луны в точке Р3-г3, тыс. км ............................. 6
Импульс характеристической скорости, сообщаемый
в точке Р1- АК1имп, м/с ........................................................................... 187,7974
Импульс характеристической скорости, сообщаемый в точке Р2- АК2имп, м/с ........................................................................... 71,4826
Импульс характеристической скорости, сообщаемый в точке р3- АК3имп, м/с ........................................................................... 290,986
Суммарная характеристическая скорость перехода АУу, м/с........... 550,2086
Начальное наклонение гиперболы подлета 7У/0, град........................ 55,435
Наклонение конечной орбиты ИСЛ у, град ......................................... 90
Масса космического аппарата тКА, кг.................................................. 1690,438
Время подлета КА к Луне = (это время пролета периселения гиперболы Т0 для высоты Ък = 100 км): 29.09.2016 03:31:48 (ИТС).
Цель работы заключается в разработке методики сквозной оптимизации траектории трехимпульсного перехода КА. Сначала на I этапе оптимизация проводится для случая импульсной тяги. Рассматриваются два варианта ориентации импульса: сначала ориентация, в «плоском случае», задается в оскулирующей плоскости углами тангажа у или атаки а, затем, в «пространственном» случае, добавляется угол выхода из плоскости, угол рыскания у. Показано, что в реальном поле с учетом возмущений от Земли, Луны и Солнца моменты оптимального приложения импульсов смещаются из апсидальных точек орбит, при этом второй импульс не сонаправлен с вектором текущей скорости. В «плоском» случае ориентации импульса получаем характеристики, близкие к апсидальному случаю импульсной тяги, в «пространственном» — конечная масса КА несколько увеличивается. На II этапе задача оптимизации трехимпульсного перелета решается для случая конечной тяги. Как и в случае импульсной тяги, результаты показывают, что во второй точке в конце разгона вектор тяги также не сонаправлен с вектором текущей скорости. Конечные массы КА в «плоском» и «пространственном» случаях получаются похожими.
Постановка задачи. Исследуем задачу оптимизации перехода (с точки зрения максимизации конечной массы КА) с селеноцентрической гиперболы подлета на высокую полярную круговую орбиту ИСЛ с заданным значением большой полуоси а^у) = ау = Ям + Ну (6 тыс. км). Начальная масса КА т0 составляет 2039,736 кг.
Методы, используемые при решении задачи оптимизации каждого из маневров. Распишем общий подход к оптимизации каждого из маневров. Как показали Д.Е. Охоцимский и Т.М. Энеев [3], для модели плоскопараллельного поля притяжения тангенс оптимального угла наклона является дробно-линейной, или линейной, функцией времени, поэтому в качестве аппроксимации закона управления рассмотрим двухпараметрическую задачу оптимальной, в рамках линейного закона, ориентации тяги по углам тангажа у(¿) или атаки а(^):
у(0 = Ю + а(0 = а0 + а'Г, (1)
в случае идеального импульсного управления у' = 0 или а' = 0.
В данном случае оптимальные значения двух параметров в рамках законов управления (1), координат двухмерного вектора х (х1, х2),
определяются в два этапа. При этом х1 соответствует г = 1 - х1 = у (или а), а х2 — г = 2 - х2 = у' (или а').
Ориентация вектора тяги задается в оскулирующей плоскости (гп) углами тангажа у или атаки а и в направлении, ортогональном плоскости (гп), — углом рыскания у (рис. 1). Здесь (гп) — плоскость, образованная ортами г° и п°, получаемыми при переходе от текущего вектора состояния в эпоху 12000.0 к вектору в орбитальной системе координат гпЬ. Угол тангажа у отсчитывается в положительном направлении от радиуса-вектора г до вектора тяги Р и определяется как у = л/2 - 9 - а, где 9 — угол наклона вектора скорости к трансвер-сали п°, у — угол наклона вектора тяги к трансверсали п°.
р *
Рис. 1. Схема определения углов ориентации тяги (стрелки показывают положительные направления углов)
Оптимальные значения двух параметров в рамках законов управления (1), координат двухмерного вектора х (х1, х2), определяются в два этапа: сначала используется метод «градиентного спуска», затем значения параметров, полученные после его применения, уточняются методом Давидсона — Флетчера — Пауэлла [4, 5]. В настоящей работе применен вариант Пауэлла; далее по тексту такую разновидность метода будем называть методом Пауэлла (квазиньютоновским методом или квазиньютоновским методом Пауэлла). Параметр х1 соответствует индексу г = 1, тогда х1 = у (или а), а х2 — индексу г = 2, т. е. х2 = у' (или а').
Сначала, в методе «градиентного спуска», минимизируется функционал Г (например, длительность работы двигателя на активном участке ф:
x(/+D = ХИ) - Xû) gradF(x(/)), (2)
где № — шаг спуска, определяемый так, чтобы F(x(/+1))<F(x(/)) или чтобы F(X(/>) была минимизирована в данной точке x(/).
Производные градиента рассчитаны численным способом через двусторонние разности.
Выполненный численный анализ показал, что поскольку изолинии в задачах разгона и торможения довольно вытянуты, а функция F(x) имеет «овражный» характер, требуется весьма много итераций для нахождения минимума. Сходимость в данном случае плохая. Задача сходится за несколько часов, а число итераций может достигать десятков тысяч. Норма градиента от итерации к итерации очень медленно уменьшается, постепенно стремясь к нулю, задаваемому как £, где £ = 10-6.
Далее используется более эффективный в рассматриваемом случае оптимизации маневров метод Пауэлла. В рамках такой разновидности квазиньютоновского метода фазовая точка сдвигается на к-м шаге
хк = хк-1 + Хкрк, (3)
где начальное положение задается, а шаг Хк определяется, как прежде, в результате минимизации функционала (длительности) в данной точке; направление сдвига
рк = -АкgradF(хк-1) = Akwk, к g N, (4)
положительно-определенная матрица Ак итерационно «подправляется» с учетом изменения градиента.
Применительно к методу Пауэлла [5]
. Aieк (AXк )T . г ...
Ак+i = Ак--к-, Ао = E, (5)
к+1 к ( Awk, Aieк )
где E — единичная матрица;
Axк = Aхк + Ак Awk, Aхк = хк - хк-1 ; Awk = gradF(хк-1 ) - gradF(хк ). (6)
Следовательно, в данном случае учитывается информация об изменении градиента, хотя и без вычисления вторых производных. За счет этого сходимость решений улучшается. Задача решается за 100-150 итераций, а время расчета составляет несколько минут.
Использование описанных выше методов оптимизации показало, что для задачи поиска оптимального трехимпульсного перехода годограф градиента gradF(x(/)) от итерации к итерации постепенно уменьшается, стремясь к нулю на плоскости координат ( dF (x)/ ôx,
ЙГ (х)/ йх^), и представляет собой кривую, близкую к прямой. Это позволяет использовать линейную аппроксимацию для параметров закона управления (у, у') или (а, а'):
х(/+2) _ хи) _ (х(/+1) -х(/^^(х>)| (7)
г г ^гасГ (х |§гас1Г (х))|'
где / _ 1, 3, 5...2п +1, п е N.
Использование в расчетах «подшагиваний» (7), улучшает сходимость задачи оптимизации градиентным методом или квазиньютоновским методом Пауэлла. Задача решается за 10-15 итераций, а время поиска решения сокращается до 2-3 мин.
Модель поля и уравнения движения КА при учете возмущений. Расчет траектории трехимпульсного перехода проводится в рамках задачи четырех тел (КА, Луна с учетом ее нецентральности в разложении в ряд 8^8, Земля и Солнце) и определяется численным интегрированием системы дифференциальных уравнений движения точки в невращающейся геоэкваториальной селеноцентрической системе прямоугольных координат 0ХУ2. При этом используется среднее равноденствие и средний геоэкватор стандартной эпохи 12000.0. Дифференциальные уравнения, описывающие движение КА, имеют вид
2— ( — — ~ Л г цм_ г - г г _ Р
__Ш. г +£ ц. г +ам +-, (8)
&2 г3 г_1 ^ | г - г |3 г3 ) т
где г — селеноцентрический радиус-вектор КА; цм — гравитационный параметр Луны; цг и гг — гравитационные параметры и радиусы-векторы возмущающих небесных тел (Земля, Солнце), где г = 1 соответствует возмущениям от притяжения Земли, а г = 2 — от притяжения Солнца; ам _ {амХ, аш, аш } — возмущающее ускорение, вызванное нецентральностью поля тяготения Луны в разложении в ряд 8x8; Р — тяга двигателя, т — масса КА.
Векторы состояния Луны и Солнца определяются из табличных эфемерид ББ-405 [6]. Использовались следующие значения констант: це = 398 600,4481 км3/с2; цм = 4902,79914 км3/с2; ц? = = 1,32712439935-Ю11 км3/с2. эфемеридная поправка ДГ(ТБВ-иТС) составляет 68,184 с.
Методика оптимизации трехимпульсного перехода на орбиту ИСЛ для «плоского» случая ориентации тяги. В реальном полете используем конечную тягу. Тогда задача оптимизации траектории трехимпульсного перехода, в общем случае, имеет 12 параметров: по четыре параметра на каждом из активных участков: два параметра
линейного закона управления вектором тяги [х1?, х2д], время включения двигателя tkTq и длительности выдачи импульсов Д^, где q — номер активного участка (импульса).
На I этапе рассмотрим идеальный импульсный случай (управление № 1), будем искать оптимальный трехимпульсный переход. Ориентацию тяги зададим прямой, т. е. у(0 = у и а(^ = а, а вместо длительностей выдачи импульсов используем их значения ДУЧ. Тогда задача оптимизации упрощается: на каждом участке будет три параметра, а их общее число сократится до девяти.
Из анализа апсидального решения [1, 2] было выявлено, что трехимпульсный переход существенно лучше по энергетическим характеристикам, чем одноимпульсный. Выигрыш апсидального трехим-пульсного перехода по сравнению с одноимпульсным (при учете гравитационных возмущений и конечности тяги) составил по суммарной характеристической скорости ~85.. ,~120 м/с и по конечной массе КА ~49 кг. ,.~70 кг. Учет возмущений приводит к новым характеристикам трехимпульсного перехода по сравнению с анализом в кеплеров-ском поле Луны. Так, например, обнаружено, что часто существует оптимальное расстояние в промежуточной точке Р2: для конечного наклонения в 90° оно равно ~45 тыс. км. Геометрически траектория получается пространственной, начальное наклонение из-за влияния возмущений сильно отличается от конечного. Разница составляет ~35°. Также определены значения импульсов для построения траек-
„ „ „Л Т^имп л тт-имп л тт-имп
тории трехимпульсного перехода: А У! , АУ2 , АУ3 . Оптимизация управления вектором тяги при реализации траектории биэл-липтического перехода заключается в декомпозиции задачи на составляющие, т. е. в отдельной оптимизации каждого из участков.
Оптимизации первого маневра. Рассмотрим задачу оптимизации первого, тормозного, импульса (точка Р1) перехода с гиперболы подлета к Луне Т0 с перицентром п1, равным г1 = 1838,57 км, на высокоэллиптическую орбиту Т1 с апоцентром а1, равным г2 = 45 тыс. км.
1. В качестве начального приближения значения первого импульса используем полученную ранее его оценку АУ1имп для апси-
дального решения.
2. Тягу направляем против нормали п к радиусу-вектору г: у01 = = 270° (а01 = 180°).
3. Задавая время приложения первого импульса Ьт 0 и параметр управления вектором тяги у01, решим однопараметрическую краевую задачу, в которой варьируется значение первого импульса ДУ1, а контролируется радиус апоселения га1 в точке Р2 на орбите Г1. Краевая задача решается с точностью е(га1) = 10-4 м. Время приложения первого импульса выбирается из диапазона ^т0 е - 2Д^0, + 2Д^0}, где —
время пролета КА периселения при движении по гиперболе Т0. Для заданных t1T 0, Yi вычисляется значение функционала F(t1T0, yi), определяемого характеристической скоростью первого импульса AF1.
4. Используя градиентный (или квазиньютоновский метод Пау-элла), итерационным образом определим параметры первого маневра: (t1T0)min, (y1)min, обеспечивающие минимум функционала F(t1T0, y1), что эквивалентно минимизации характеристической скорости первого импульса — (AV1)min.
Оптимизация второго и третьего маневров. Зависимость параметров третьего маневра и конечной массы КА mf от второго — довольно велика, поэтому следует рассматривать их совместную оптимизацию. В ходе такого анализа выяснилось, что оптимальная по конечной массе mf точка приложения второго импульса (точка Р2) сильно «уходит» из апоселения а1 высокоэллиптических орбит Т1 и Т2.
Для оптимизации второго маневра рассмотрим задачу оптимизации второго, разгонного, импульса (точка P2) перехода с орбиты Т1 с апоцентром а1, равным 45 тыс. км, на орбиту Т2 с перицентром п2, равным 6 тыс. км.
1. Зададим время приложения второго импульса t2T0 из диапазона t2T0 е [ta1 - ô, ta1 + S], где ta1 — время пролета апоселения, 5 = = 80 тыс. с.
2. В качестве начального приближения для величины второго
Л т^имп
импульса используем полученную ранее его оценку для апси-
дального решения.
3. Вектор тяги направим по нормали n к радиусу-вектору r, т. е. Y02 = 90° (ао2 = 0°).
4. Задавая время приложения второго импульса t2T0 и параметр управления y2, решим однопараметрическую краевую задачу, в которой варьируется значение второго импульса AV2, а контролируется радиус периселения гп2 в точке P3 на орбите Т2. Краевая задача решается с точностью e(rn2) = 10-4 м. Находим значение функционала F(Y2), задаваемого характеристической скоростью второго импульса AV2.
5. Применяя градиентный (или квазиньютоновский метод Пау-элла), итерационно решаем задачу оптимизации второго маневра, минимизируя функционал F(y2). В результате определяются параметры (t2T0)min, (y2)min, обеспечивающие минимум характеристической скорости второго импульса (AV2)min.
Для оптимизации третьего маневра рассмотрим задачу оптимизации третьего, тормозного, импульса (точка P3) перехода КА с орбиты Т2 на высокую круговую орбиту ИСЛ Tf с большой полуосью af, равной 6 тыс. км.
1. В качестве начального приближения для третьего импульса используем полученную ранее его оценку АК$имп для апсидального решения.
2. Тягу направляем против нормали n к радиусу-вектору r: у03 = = 270° (аоз = 180°).
3. Задавая параметр управления тяги у03, решим краевую задачу перехода с высокоэллиптической орбиты T2 на конечную круговую орбиту Tf. Данный переход можно выполнить по-разному. Однако учитывая то, что вблизи конечной круговой орбиты может возникнуть ситуация недифференцируемости функции (например, функций эксцентриситета e, радиусов пери- r„ и апоселения ra), для решения задачи перехода необходимо использовать непрерывную монотонно возрастающую вблизи перицентра и легко дифференцируемую функцию радиальной скорости Vr.
3.1. Третий импульс AV3 зададим условием достижения заданной величины большой полуоси af конечной орбиты. Итерационно сведем к нулю значение радиальной скорости Vr (с точностью £(Vr) = 0,1 мм/с) после сообщения импульса. В качестве варьируемого параметра выберем время начала работы двигателя t3T0 до пролета перицентра. Сходимость получается хорошей. Эта «внутренняя» задача решается за одну-две итерации.
3.2. Далее определим апсидальное расстояние ra3 полученной орбиты Tf и его отклонение от конечного радиуса af : Ara = ra3 - af. Если Ara больше допустимого, изменим расстояние в периселении высокоэллиптической орбиты Т2, т. е. rn2(i+1) = rn2i + Ara, возвратимся к этапу оптимизации второго импульса, а затем и третьего импульса, и так будем повторять до тех пор, пока в конце этапа оптимизации третьего импульса не выполнится условие: ra3 = af — в пределах заданной точности £(Arf) = 2 м. Сходимость получается хорошая. Эта «внешняя» задача решается за одну-две итерации. Находим значение функционала F(y3, у3), задаваемого характеристической скоростью третьего импульса AV3.
3.3. В результате для заданного времени сообщения второго импульса t2T0 находим параметры: t3T0, (y3)min, (y3)min, обеспечивающие минимум характеристической скорости третьего импульса (AV3)min.
Полученная конечная масса КА mf является максимальной для заданного времени t2T0, т. е. локальным максимумом. Для нахождения «глобального» — необходимо проварьировать время приложения второго импульса t2T 0. Опишем алгоритм его поиска.
Поиск оптимального времени приложения второго импульса — tir. Задавая нулевое приближение t2Ti = tn - 8i с (где 5г- = 0,5At2) и используя приведенный выше алгоритм совместной оптимизации второго
и третьего импульсов, определим конечную массу КА т-р. В смещен-
йт-г
ной точке Р2 вычислим двухстороннюю производную А _ ——:
йХ
йт- _ т- (Х2Тг +5Х )-т/ (Х2Тг _5Х )
dt 2St
(9)
йтг
При вычислении производной - следует учитывать, что если
йХ
массы справа т-(Х2Т0 + 5Х) и слева т-(Х2Т0 - 5Х) меньше начальной т/(Х2Т0), то величина варьирования уменьшается в 2 раза, т. е. 5Х(Ш) =
( йтг Л
= 0,55tk, и определяется заново:
'f
dt Jk
. Ее уменьшение продолжает-
ся до тех пор, пока выполняется условие > 1 с. Если < 1 с, то поиск оптимального Х2Т заканчивается.
Если одна из масс (справа или слева) больше начальной, а также если 5^ > 1 с, то определяется величина, обратная производной А:
A-1 =
(dmf ^ 1 dt
dt
A 1 =-, затем ее приращение 5(г+1) = A 1. Далее
^ j dmf
вычисляется время приложения второго импульса на следующей итерации t2T(i + 1) = At2Ti + + 1).
После этого для заданного времени t2T(ii + 1) находим массу КА mfi + 1). Если mf(i+1) > m^, то определяется приращение для следующей итерации 5(г- + 2), а далее определяется время приложения второго импульса t2T(i+2) = t2T(i+1) + 5(г+2) и конечная масса КА mfv+2). Иначе, если mf(i+1) < mfi, то приращение на текущей итерации уменьшается в 2 раза: 5(i+1)* = 0,5, пересчитывается момент приложения второго импульса t2T(i+1) = t2T(i+1) + 5(i+1) и алгоритм поиска массы повторяется. Итерации продолжаются, пока выполняется условие 5(i+1) > em, где £m = 1 c. В итоге определяем набор параметров, обеспечивающих максимум конечной массы КА (mf)max - [(t2T)min]max, [(AV2)min]max, [(Y2) min max [(t3T) min max [(AV3) min max [(Y3) min max
После решения задачи оптимизации и определения всех необходимых для трехимпульсного перехода параметров: (t1T)min, (AtT1)min,
(Yl)min, [(t2T)min] max [(AtT2)min] max [(Y2)min] max [(t3T)min] max [(AtT3)min] max
[(y3)min[max — следует вычислить значение наклонения орбиты i3, про-
a -m
верить его рассогласование Aif с заданным значением конечного наклонения if. Если выполняется условие Af > £i, где £i = 0,001°, то
a -m -(m+1)
подлетное наклонение i0 изменяется на значение Aif: i0 =
= ¡т + А/, затем алгоритм оптимизации повторяется. Так происходит до тех пор, пока не выполнится условие: А/™ < £/. Задача сходится за пять-шесть глобальных итераций.
Основные результаты оптимизации трехимпульсного перехода на орбиту ИСЛ приведены в табл. 1.
Таблица 1
Результаты оптимизации трехимпульсного перехода на орбиту ИСЛ квазиньютоновским методом Пауэлла при управлении по углу у
Параметр Импульс
первый второй третий
Момент приложения импульса Тимп (итс) 29.09.2016 3:32:32 01.10.2016 9:06:15 03.10.2016 4:48:15
Время до пролета апсидальных точек* ^ с -6,208 29216,805 -0,168
Конечная масса КА ту, кг 1912,9069 1867,2673 1690,6738
Импульс скорости ДУХ, м/с /5 = 52,984256° 188,0468 70,7356 291,0173
Угол тангажа в момент сообщения
импульса у, град Д/ = 0,0009838° 270,1368 91,0197 270,0017
Угол атаки в момент сообщения импульса а, град 180,13553 19,52027 180,00906
Расстояние до центра Луны гу, км /5 = 52,9852398° 1838,572 44121,591 6000,001
* Для импульсов: первого — я1, второго — а1, третьего — я2.
Из табл. 1 видно, что оптимальные моменты приложения импульсов смещаются из апсидальных точек орбит. Для первого и третьего импульсов сдвиг по времени относительно перицентров мал по сравнению с периодом орбиты, для второго он составил ~8 ч. При этом точки приложения первого и третьего импульсов смещаются в область до пролета периселения, а третьего — в область после пролета апоселения. В данном случае при движении КА по траектории трехимпульсного перехода селенографическое наклонение под действием возмущений (см. формулу (8)) растет. Поэтому смещение точки приложения второго импульса Р2 по орбите Т1 в область после пролета апоселения ведет к увеличению наклонения /3 и уменьшению подлетного /0 на ~2,4°. Углы тангажа у и атаки а после приложения первого и третьего импульсов составили у = 270° и а = 180°, а после приложения второго — у = 91° и а = 19,8°.
Следовательно, для оптимальной траектории трехимпульсного перехода точка Р2 сильно «уходит» из апоселения, а тяга в ней не направлена по текущей скорости.
Значение подлетного наклонения составило i0 = 52,985°, время подлета КА к Луне t0: 29.09.2016 3:31:54 (UTC). Характеристические затраты на переход составили 549,799 м/с и конечная масса КА равна 1690,6738 кг.
Описанная выше методика была применена к задаче исследования оптимального трехимпульсного перехода на орбиту ИСЛ в случае импульсной тяги двигателя. При этом расстояние в удаленной точке ra изменялось в диапазоне от 6 до 55 тыс. км (табл. 2).
Методика оптимизации трехимпульсного перехода на орбиту ИСЛ для «пространственного» варианта ориентации тяги. Рассмотрим «пространственный» вариант ориентации тяги в случае импульсной тяги. Для этого при задании ориентации тяги для второго импульса добавим угол выхода из плоскости, угол рыскания у (управление № 2). Число оптимизируемых параметров увеличивается с 9 до 10.
Тогда методика оптимизации трехимпульсного перехода в общих чертах похожа на описанную ранее, за исключением того, что при оптимизации второго импульса вместо однопараметрической решается двухпараметрическая краевая задача. В ней варьируются значения второго импульса AV2 и угла рыскания у, а контролируются значения радиуса периселения r„2 и наклонения i3 в точке P3 орбиты Т2. Краевая задача решается по радиусу с точностью e(rn2) = 10-4 м и по наклонению i3 е (i3) = 10-4 град (см. табл. 2).
Анализ данных табл. 2 показывает, что применение методики оптимизации для «плоского» случая при варьировании ra дает результаты, похожие на апсидальное решение [1, 2]. Масса КА на конечной орбите mf1 получается немного больше массы при апсидальном решении m™n : в зависимости от ra выигрыш составляет от ~1 до 270 г.
Видно, что с увеличением ra возрастает число итераций Ыит, которые нужно провести для решения задачи: при увеличении расстояния ra до 55 тыс. км их число возрастает до 10.
Сравнительный анализ величин mf1 и mf2 показывает, что учет угла рыскания у позволяет увеличить конечную массу КА на ôm е [18; 846] г, сократить число итераций Ыит по конечному наклонению до одной, а также существенно сократить поиск оптимального решения. Угол рыскания у, в зависимости от расстояния в удаленной точке ra, изменяется в пределах от 2 до 4°.
Результаты оптимизации трехимпульсного перехода при варьировании г„ для «плоского» и «пространственного» случаев управления вектором тяги
Параметр Управление вектором тяги Расстояние от центра Луны до удаленной точки г„, тыс. км
6 10 15 20 25 30 35 40 45 50 55
Конечная масса тд, кг Дня «плоского» случая 1597,5411 1632,8938 1654,5681 1666,4437 1674,4162 1680,4145 1685,5289 1689,6399 1690,6738 1686,3804 1679,3002
Подлетное наклонение 'оь град 89,99900 89,72976 89,24187 89,18913 87,76817 84,67780 78,45617 67,34822 52,98524 40,23127 28,45150
Конечная масса та, кг Для «про-стран-ственно-го» случая 1597,7653 1633,0172 1655,0237 1667,2719 1674,4162 1680,4138 1685,5268 1689,7449 1691,1115 1686,9684 1680,1328
Подлетное наклонение г02, град 89,99698 89,95672 89,76762 89,18917 87,76929 84,68870 78,57969 68,19951 55,4350 44,7083 38,5277
Конечная масса имп т ^ , кг Для ап-сидаль-ного решения 1597,5401 1632,7569 1654,3002 1666,4264 1674,3988 1680,3947 1685,4989 1689,6099 1690,4405 1686,3785 1679,2990
Подлетное наклонение г'о, град 89,99698 89,95672 89,76762 89,18917 87,76929 84,68870 78,57969 68,19951 55,4350 44,7083 38,5277
Количество итераций Л^ИТ — 1 1 1-2 1-2 2-3 2-3 3^1 4-5 5-6 6-8 8-10
Методика оптимизации «трехимпульсного» перехода на орбиту ИСЛ для «плоского» варианта ориентации тяги. На II этапе рассмотрим методику оптимизации трехимпульсного перехода при учете конечности тяги двигателя. Такой переход назовем «трехим-пульсным»1. Для расчетов воспользуемся двигательной установкой со следующими характеристиками: значение тяги считается постоянным и равным Р = 4118,793 (~4119) Н (420 кгс), а удельный импульс составляет Руд = 298,7 с. С учетом линейного закона управления (1) число оптимизируемых параметров равно 12.
Алгоритм оптимизации «трехимпульсного» перехода при варьировании точек приложения (управление № 3) похож на приведенный выше для идеального трехимпульсного перехода. Рассмотрим его отличия:
1) для решения краевых задач на этапах оптимизации маневров вместо значений импульсов AVq варьируются длительности торможения Atq на q-м участке, первое приближение для которых определяется по Циолковскому [7];
2) на этапе оптимизации третьего маневра выход на заданное значение большой полуоси конечной орбиты af определяется итерационно, а длительность торможения считается путем вычисления производной da/dt:
da _ 2a2 2 Р
VT-cos
dt MM T mT,
К), (8)
Aa-
At3(i+1) = At3i+AtС, AtQ = ~da , (9)
dt
aT - af
= f 2 f--, (10)
2a2 2 P
4— cos (атг)
V
Mm mT
где аг — угол атаки в конце торможения, г е N.
Итерации заканчиваются, когда Аа, < еа, еа = 10-6 м. Задача также сходится за пять-шесть глобальных итераций. Результаты оптимизации маневров градиентным методом (управление № 3а) и квазиньютоновским методом Пауэлла представлены в табл. 3.
По характеристикам параметры оптимизации (см. табл. 3) при управлении по углу тангажа у для конечной тяги похожи на случай импульсной тяги двигателя. Гравитационные потери на всех активных участках малы и не превышают 5 см/с, т. е. 5Кгр < 0,05 м/с.
1 Выделим кавычками для различия.
Результаты оптимизации маневров градиентным методом (управление № За) и методом Пауэлла (управление № 36) по углу у (/V,,, = 5)
Управление № За Управление № 36
Параметр 1 2 3 1 2 3
29.09.16 03:32:32 01.10.16 09:14:32 03.10.16 04:53:17 29.09.16 03:32:32 01.10.16 09:06:19 03.10.16 04:51:37
Время до пролета апсидальной точки /,,-/■. с -51,782 29278,564 -85,98915 -51,785 29218,064 -83,79376
Длительность маневра А1т, с 90,2121 32,4739 125,5942 90,2117 32,4727 125,5854
Масса КА после пролета активного участка ту, кг 1912,8894 1867,2281 1690,6381 1912,8898 1867,2302 1690,6454
Характеристические затраты на маневр А1м/с 188,0737 70,7700 291,0177 188,0730 70,7676 291,0083
Параметры закона управления ориентацией вектора тяги у0, град 271,9028 91,0226 269,9944 271,9429 91,0378 269,9806
Уо,°/с -0,036787 0,000206 0,007587 -0,036787 0,000193 0,007468
Расстояние до центра Луны в конце активного участка г,. км 1839,374 44091,685 5999,999 1839,374 44119,648 5999,996
Угловая дальность маневра Ли,. град 6,625 0,007 1,262 6,625 0,006 1,261
Примечание. 1, 2, 3 — номер активного участка (а. у.).
Подлетное наклонение /0 уменьшается с 55,435° до 52,9751°, а «выигрыш» по конечной массе КА по сравнению с апсидальным решением (в точке Р2) составляет 0,233 кг.
Результаты применения квазиньютоновского метода Пауэлла (управление № 3а) (см. табл. 3) незначительно отличаются от данных, полученных с помощью градиентного метода (управление № 3б).
На рис. 2, а, рис. 3, а и рис. 4, а приведены угловые зависимости для каждого из маневров после оптимизации траектории квазиньютоновским методом Пауэлла при задании ориентации тяги по углу тангажа у (управление № 3а). Отметим, что оптимизация «трехим-пульсного» перехода с помощью градиентного метода при управлении по углу тангажа у (управление № 3б) дает похожие угловые зависимости.
у, а, град у, а, град
а б
Рис. 2. Основные углы ориентации тяги вдоль первого маневра:
а — при управлении по углу у; б — при управлении по углу а
у, а, град
30 20 10 0 -10 -20 30 -40
у,,,а, град
30 г, с
10
15
б
г —-- а
" у- 90°
_ Ф
— | 1 1111
20 25 30 г, с
Рис. 3. Основные углы ориентации тяги вдоль второго маневра:
а — при управлении по углу у; б — при управлении по углу а
у, а, град
1 0
у, а, град
у - 270°
СТ 0
1 1 1 1 1
0 20 40 60 80 100 а
г, с
- у - 270°
а -к
| ф-к | 1 |
0 20 40 60 80 100 г, с б
Рис. 4. Основные углы ориентации тяги вдоль третьего маневра:
а — при управлении по углу у; б — при управлении по углу а
На рис. 2-4 представлены зависимости углов тангажа у (красный цвет), атаки а (зеленый цвет), угла наклона вектора тяги ф (сиреневый) к трансверсали п° и угла наклона вектора скорости 9 (синий) к трансверсали п° (см. рис. 1). Некоторые зависимости смещены на определенный угол для большей наглядности:
• на рис. 2 у - 90°, 9 + 180°;
• на рис. 3 у - 90°;
• на рис. 4 у - 270°, а - 180°, ф -180°.
Была также проведена оптимизация «трехимпульсного» перехода с помощью квазиньютоновского метода Пауэлла при задании ориентации тяги по углу атаки а, когда а01 = 180°, а'01 = 0 °/с, а02 = 0°, а 02 = = 0 °/с, а03 = 180°, а'03=0 °/с (управление № 3в). Анализ показал схожие результаты. Угловые зависимости для всех активных участков представлены на рис. 3, 5, 7.
Сравнение рис. 2, а, рис. 3, а и рис. 4, а с рис. 2, б, рис. 3, б и рис. 4, б показало, что управление тягой по углу атаки а позволяет получить угловые зависимости, похожие на управление вектором тяги с помощью угла тангажа у. Конечная масса КА при управлении по углу а — т/а = 1690,6445 кг примерно равна таковой при управлении по углу у — да/у = 1690,6454 кг.
Оптимальные углы ориентации тяги вследствие действия возмущений от гравитационных полей Земли и Солнца показывают характеристики, отличающиеся от апсидального решения. В конце первого маневра вектор тяги направлен против скорости (а1 = 180° ). Задание в качестве начального приближения для линейного закона ориентации тяги в смещенной точке Р2 простейшего управления а02 = 0°, а 02 = 0 град/с при дальнейшей оптимизации приводит к значениям а2 = 31,173° а'2 = -0,3835 град/с, а угол атаки в конце второго маневра составляет а2/ = 18,72°. Таким образом, в точке Р2 сообщается существенная положительная радиальная составляющая скорости АУГ > 0,
кото dV„
dt
рая уменьшает производную радиальнои скорости по времени Это ведет к более медленному уменьшению радиальной ско-
рости Vr при движении КА по орбите Т2, изменению эксцентриситета и формы орбиты. В результате время разгона КА в смещенной точке Р2 получается меньше, чем при сообщении второго импульса в апси-дальной точке ai.
В начале третьего активного участка траектории тяга направлена точно против скорости (a3 = 180°), что хорошо согласуется с условием Лоудена о коллинеарности векторов тяги и скорости в начальный момент торможения PffV, а также с результатами, полученными ранее, при анализе задачи торможения в случае прямого перехода на орбиту ИСЛ [4, 7-10].
На рис. 5 представлена схема трехимпульсного перехода при сообщении второго импульса в смещенной точке Р2.
Рис. 5. Схема трехимпульсного перехода при сообщении второго импульса в смещенной точке Р2 (размер и масштаб не соблюдены)
Рз
Параметры конечной орбиты ИСЛ после оптимизации траектории «трехимпульсного» перехода с помощью квазиньютоновского метода Пауэлла при управлении по углу а приведены ниже:
Время подлета КА к Луне /0 ....................................................... 29.09.16 03:31:43
Начальное наклонение гиперболы Т0 /0, град ........................... 52,989
Радиус апоселения г— км ........................................................... 5999,999
Радиус периселения г—, км ......................................................... 6000,001
Долгота восходящего узла орбиты ИСЛ О, град..................... 195,904
Наклонение конечной орбиты ИСЛ ф град .............................. 89,999
Время перехода на орбиту ИСЛ ТФ............................................. 03.10.16 04:58:48
Характеристические затраты на трехимпульсный
переход Ух, м/с ............................................................................. 549,846
Конечная масса КА т,, кг ........................................................... 1690,6454
Методика оптимизации «трехимпульсного» перехода на орбиту ИСЛ для «пространственного» варианта ориентации вектора тяги. Методика оптимизации «пространственного» «трехимпульсного» перехода (управление № 4) похожа на описанную ранее (управление № 2). Оптимизируемых параметров 13. Краевые задачи
на первом и третьем активных участках решаются так же, как и при управлении № 3, а на втором - как и при управлении № 2, но вместо значения импульса ДУ2 варьируется длительность его сообщения Д^т, а ориентация тяги задается линейным законом управления (1). Результаты оптимизации в случае управления № 4 приведены в табл. 4.
Таблица 4
Результаты оптимизации маневров квазиньютоновским методом Пауэлла для «пространственного» варианта ориентации вектора тяги
Параметр 1 2 3
29.09.16 03:32:32 01.10.2016 06:52:41 03.10.2016 04:50:02
Время до пролета апсидальной точки t0T, с -51,785 21617,872 -60,88122
Длительность активного участка AtT, с 90,101 31,417 125,743
Масса КА после пролета активного участка mf, кг 1913,0456 1868,8704 1692,0631
Характеристические затраты на активный участок AVX, м/с 187,8345 68,4342 291,1253
Параметры закона управления ориентацией вектора тяги у0, град у' о, град/с 271,9172 -0,036796 91,8201 0,017993 269,9999 -5,34E-05
Угол рыскания у, град — 2,857 —
Расстояние до центра Луны в конце активного участка f км 1839,37 44515,39 0,006
Угловая дальность маневра Auf, град — 6000,00 1,264
Учет угла рыскания в момент сообщения второго импульса для оптимизации трехимпульсного перехода на высокую орбиту ИСЛ приводит к увеличению конечной массы КА на ~1 кг по сравнению с «плоским» случаем ориентации тяги и на ~1,6 кг по сравнению с ап-сидальным решением.
Энергетические характеристики перехода и масса КА для всех видов управления представлены в табл. 5. Их значения близки друг к другу: скорость отличается в пределах 2 м/с, а конечная масса — в пределах 1,6 кг.
Выводы. Разработана методика сквозной оптимизации трехимпульсного перехода КА. Проведено сравнение разных вариантов управления вектором тяги. При проведении анализа трехимпульсного перехода для случая импульсной тяги двигателя показано, что в реальном поле при учете возмущений от Земли, Луны и Солнца оптимальные моменты приложения импульсов смещаются из апсидаль-ных точек орбит, при этом второй импульс не направлен по вектору текущей скорости. Для точек приложения первого и третьего импульсов сдвиг по времени относительно перицентра линии апсид мал по сравнению с периодом орбиты, а для второго импульса он довольно большой и составляет ~8 ч.
Массоэнергетические характеристики задачи оптимизации трехимпульсного перехода для всех вариантов управления
Параметр 1 2 3 4 5 6
Управление № 1 Управление № 2 Управление № За Управление № 36 Управление № Зв Управление № 4
Характеристические затраты на активный участок \1м/с 549,799 549,041 549,861 549,848 549,846 547,394
Конечная масса КА ту, кг 1690,6738 1691,1115 1690,6381 1690,6454 1690,6456 1692,0631
«Выигрыш» трехимпульсного перехода по сравнению с ап-сидальным решением Ату кг 0,2358 0,6735 0,2001 0,2074 0,2076 1,6251
Расстояние в удаленной точке (а 1) Г2, тыс. км 45
Подлетное наклонение /0, град 52,985 55,435 52,995 52,982 52,989 55,435
Наклонение конечной орбиты ИСЛ //. град Примечание. 1-6 — номер вар 90 ианта.
Задача оптимизации траектории трехимпульсного перелета решена для случая конечной тяги. Проведено сравнение градиентного метода и квазиньютоновского метода Пауэлла для случая конечной тяги, показавших похожие результаты.
Использование градиентного метода и квазиньютоновского метода Пауэлла совместно с линейной аппроксимацией в случае конечной тяги при незначительных гравитационных потерях (ôFrp = = 0,05 м/с) обеспечивает хорошую сходимость задачи оптимизации. В качестве начальных приближений законов управления на всех участках можно задать простейшие законы управления: у0 = 270°, у0 = = 0°/с, а0 = 180°, а0 = 0°/с — для первого и третьего; у0 = 90°, у0 = = 0°/с, а0 = 0°, а0 = 0°/с — для второго. При выполнении оптимизации каждого из активных участков вне зависимости от метода оптимизации задача сходится за 10-15 итераций, время расчета составляет около 2-3 мин.
Действие возмущений от гравитационных полей Земли и Солнца приводит к тому, что траектории трехимпульсного и «трехимпульсного» перехода КА сильно эволюционируют. Особенно сильно они изменяются вблизи апсидальной точки а1. Поиск оптимальных закона управления и точки приложения второго импульса (точки Р2) смещает ее в область пролета апоцентра а1.
Оптимальный закон управления по углу атаки а на втором активном участке a(t) = 31,173 — 0,3835 град/c свидетельствует об изменении формы орбиты в смещенной точке Р2. Таким образом, используя действие возмущений от гравитационных полей Земли и Солнца, можно сократить значение второго импульса по сравнению с апси-дальным решением.
Реальный трехимпульсный переход по орбитам вокруг Луны для «плоского» случая оптимизации при учете конечности тяги можно выполнить практически с теми же затратами, что и для «плоского» варианта с импульсной тягой. Гравитационные потери на всех активных участках малы и не превышают 5 см/с, т. е. 5Кгр < 0,05 м/с. Выигрыш «плоского» варианта ориентации тяги по сравнению с апси-дальным решением по конечной массе КА составил ~0,233 кг. При этом, как и в случае импульсной тяги, результаты показали, что во второй точке, в конце разгона, вектор тяги также не направлен по вектору текущей скорости. Это отличает данный случай оптимизации от классического, описанного Д.Ф. Лоуденом [9].
Параметры оптимизации при «плоском» управлении по углам у и а для конечной тяги по характеристикам похожи между собой и лишь немного отличаются от апсидального решения.
В начале третьего активного участка траектории тяга направлена точно против скорости (а3 = 180°), что хорошо согласуется с усло-
вием Лоудена о коллинеарности векторов тяги и скорости в начальный момент торможения P\\V, а также с результатами, полученными ранее в ходе анализа задачи торможения при прямом переходе на орбиту ИСЛ [4, 7-10].
Варьирование расстояния в удаленной точке а1 - ra показало, что применение методики оптимизации дает результаты, похожие на полученные для апсидального решения [1, 2]. Масса КА mf на конечной орбите получается немного больше. В зависимости от расстояния в удаленной точке а1 «выигрыш» составляет от 1 до 270 г. Замечено, что чем больше расстояние в удаленной точке ra, тем больше нужно сделать глобальных итераций; при увеличении ra до 55 тыс. км количество итераций Ыит возрастает до 10.
Введение угла выхода из плоскости, угла рыскания у для задания ориентации вектора тяги на втором активном участке при учете конечности тяги увеличивает конечную массу КА на ~1 кг по сравнению с «плоским» случаем ориентации тяги и на ~1,6 кг в сопоставлении с апсидальным решением.
Значения энергетических характеристик перехода и конечной массы КА для всех видов управления близки одно к другому: скорость отличается не более чем на 2 м/с, а конечная масса — не более чем на 1,6 кг. Геометрически рациональная траектория при учете возмущений, углов и точек приложения импульсов сильно отличается от траектории, полученной для апсидального импульсного решения [1, 2].
ЛИТЕРАТУРА
[1] Гордиенко Е.С., Ивашкин В.В. Анализ оптимального трехимпульсного перехода на орбиту искусственного спутника Луны. Инженерный журнал: наука и инновации, 2016, вып. 3. DOI 10.18698/2308-6033-2016-3-1472
[2] Гордиенко Е.С., Ивашкин В.В. Использование трехимпульсного перехода для выведения космического аппарата на высокие орбиты искусственного спутника. Космические исследования, 2017. (В печати).
[3] Охоцимский Д.Е., Энеев Т.М. Некоторые вариационные задачи, связанные с запуском искусственного спутника Земли. Успехи физических наук, 1957, т. 63, № 1а, с. 5-32.
[4] Гордиенко Е.С., Ивашкин В.В., Лю В. Анализ оптимальных маневров разгона и торможения космического аппарата при его полете к Луне. Космонавтика и ракетостроение, 2015, вып. 1, с. 37-47.
[5] Аттетков А.В., Галкин С.В., Зарубин В.С. Методы оптимизации. В.С. Зарубин, А.П. Крищенко, ред. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2003, 440 с.
[6] Standish E.M. JPL Planetary and Lunar Ephemerides. Interoffice memorandum: JPL IOM 312, F, 1998, August 26, рр. 98-048.
URL: ftp://ssd.jpl.nasa.gov/pub/eph/planets/ioms/de405.iom.pdf (дата обращения 27.02.2017).
[7] Циолковский К.Э. Труды по ракетной технике. Москва, Изд-во Оборонгиз, 1947, 368 с.
[8] Цандер Ф.А. Проблемы полета при помощи реактивных аппаратов. Москва, Госавиац. и автотрактизд., 1932, 75 с.
[9] Лоуден Д.Ф. Оптимальные траектории для космической навигации. Москва, Мир, 1966, 152 с.
[10] Hohmann W. Die Erreichbarkeit der Himmelskörper. München, Verlag Oldenbourg, 1925. ISBN 3-486-23106-5.
Статья поступила в редакцию 20.02.2017
Статья подготовлена по материалам доклада, представленного на XLI Академических чтениях по космонавтике, посвященных памяти академика С.П. Королёва и других выдающихся отечественных ученых — пионеров освоения космического пространства. Москва, МГТУ им. Н.Э. Баумана,
24-27января 2017 г.
Ссылку на статью просим оформлять следующим образом: Гордиенко Е.С. Исследования оптимального трехимпульсного перехода на высокую орбиту искусственного спутника Луны. Инженерный журнал: наука и инновации. 2017, вып. 9. http://dx.doi.org/10.18698/2308-6033-2017-9-1667
Гордиенко Евгений Сергеевич — инженер НПО им. С. А. Лавочкина; аспирант кафедры «Динамика и управление полетом ракет и космических аппаратов» МГТУ им. Н.Э. Баумана. Автор 10 научных работ в области моделирования, построения и оптимизации лунных траекторий. е-mail: gordienko.evgenyy@gmail.com
Investigating the optimal three-impulse transfer into the high orbit of the artificial lunar satellite
© E.S. Gordienko
Lavochkin Science and Production Association, Khimki town, Moscow region, 141402, Russia Bauman Moscow State Technical University, Moscow, 105005, Russia
The article examines the problem of optimum spacecraft ascent from the Earth into the high circular polar orbit of the artificial lunar satellite with the radius of 6000 km using the three-impulse bi-elliptic scheme of the so called "Shternfeld" transfer. We have carried out this analysis by taking into account the disturbances from the lunar field noncen-trality, the Earth and the Sun gravitational fields as well as the engine thrust finiteness. The optimum trajectory is defined by varying both impulse control parameters and their application points. The analysis consists of two stages. At the first stage we consider two options of impulse orientation in the ideal impulse occurrence: at first the impulse orientation is set in the osculating plane by "pitch" у and "attack" a angles; then the "yaw" angle у (angle of departure from the plane) is added. It is shown that with the increase of the maximum distance rathe optimal points of the impulses application are shifted from the apsidal points of the orbits. At the same time the second, intermediate impulse is not directed along the current velocity vector. In the first, "flat" option of the impulse orientation we get characteristics which are similar to the ideal impulse apsidal case. In the second, "spatial" variant the spacecraft remaining mass slightly increases. At the second stage we take into account the engine thrust finiteness. The results obtained are very close to the impulse case.
Keywords: spacecraft, trajectory optimization, artificial lunar satellite, three-impulse bi-elliptic scheme, "Sternfeld" transfer, Powell's quasi-Newton method
REFERENCES
[1] Gordienko E.S., Ivashkin V.V. Inzhenernyy zhurnal: nauka i innovatsii — Engineering Journal: Science and Innovation, 2016, iss. 3.
DOI: 10.18698/2308-6033-2016-3-1472
[2] Gordienko E.S., Ivashkin V.V. Kosmicheskiye issledovaniya — Cosmic Research, 2017 (is to be published).
[3] Okhotsimskiy D.E., Eneyev T.M. Uspekhi fizicheskikh nauk — Advances in Physical Sciences, 1957, vol. 63, no. 1a, pp. 5-32.
[4] Gordienko E.S., Ivashkin V.V., Lyu V. Kosmonavtika i raketostroenye — Cosmonautics and Rocket Engineering, 2015, no. 1, pp. 37-47.
[5] Attetkov A.V., Galkin S.V., Zarubin V.S. Metody optimizatsii [Optimization procedure]. V.S. Zarubin, A.P. Krischenko, eds. Moscow, BMSTU Publ., 2003, 440 p.
[6] Standish E.M. JPL Planetary and Lunar Ephemerides. Interoffice memorandum, 1998. Available at: ftp://ssd.jpl.nasa.gov/pub/eph/planets/ioms/de405.iom.pdf (accessed February 27, 2017).
[7] Tsiolkovskiy K.E. Trudy po raketnoy tekhnike [Works on rocket technology]. Moscow, Oborongiz Publ., 1947, 368 p.
[8] Tsander F.A. Problemy poleta pri pomoshchi reaktivnykh apparatov [Problems of flight by jet propulsion: interplanetary flights]. Moscow, Gos. Aviats i avto-trakt Publ., 1932, 75 p.
[9] Lawden D.F. Optimal Trajectories for Space Navigation. London, Butterworths, 1963. [In Russ.: Louden D.F. Optimalnye traektorii dlya kosmicheskoy navi-gatsii. Moscow, Mir Publ., 1966, 152 p.]
[10] Hohmann W. Die Erreichbarkeit der Himmelskörper. München, Verlag Oldenbourg, 1925.
Gordienko E.S. (b. 1990) graduated from Bauman Moscow State Technical University in 2013. Engineer of Lavochkin Science and Production Association. Post-graduate student of the Department of Dynamics and Flight Control of Rockets and Spacecraft, Bauman Moscow State Technical University. Author of 10 scientific papers in the field of lunar trajectory modelling and plotting. e-mail: gordienko.evgenyy@gmail.com
