CC BY
17УДК 517.977.58
DOI 10.18698/2308-6033-2017-7-1639
Оптимизация экспедиции к Фобосу космического аппарата с комбинированной тягой с возвращением к Земле
© И.С. Григорьев1, М.П. Заплетин1,2, А С. Самохин1,2,3, М.А. Самохина2,3
1МГУ им. М.В. Ломоносова, Москва, 119991, Россия 2РУДН, Москва, 117198, Россия 3НИУ ВШЭ, Москва, 101000, Россия
Исследован актуальный вопрос разработки методики оптимизации межпланетных траекторий космических аппаратов с проведением сквозной оптимизации всей миссии. Рассмотрена задача оптимизации экспедиции пространственного перелета космического аппарата, оснащенного двигателями большой и малой тяги, к Фобосу с возвращением к Земле с учетом эфемерид, приуроченная к реальной миссии, которую планируется осуществить в недалеком будущем. Представлена методика построения экстремалей Понтрягина с учетом притяжения планет на планетоцентрических участках траектории. Задача космодинамики формализована как задача оптимального управления и решена численно методом стрельбы. Построены конкретные траектории, оценен возможный выигрыш при использовании комбинированной тяги, а не только двигателей большой тяги.
Ключевые слова: малая тяга, комбинированная тяга, сквозная оптимизация, оптимизация межпланетных перелетов, оптимизация траектории космического аппарата, Фобос, полет к Фобосу, экспедиция к Фобосу, экспедиция к Марсу, задача космодинамики
Введение. В настоящее время управление космическим аппаратом (КА) посредством только реактивных двигателей большой тяги не позволяет доставить к планете-цели удовлетворительную массу полезного груза. Значительно большую массу полезной нагрузки можно доставить, используя только двигатели малой тяги. Но при этом время ухода КА от Земли и его торможения у планеты-цели увеличится и может составить большую часть времени перелета, став критическим с точки зрения отказоустойчивости всей системы. Достоинством использования двигателей большой тяги является малое время перелета, а малой тяги — большая масса полезной нагрузки. Использование их комбинации объединяет достоинства маневров с использованием каждого из этих двигателей, позволяет избежать указанные недостатки и таким образом реализовать доставку необходимой массы полезной нагрузки за приемлемое время.
Указанная проблема исследована в данной статье на примере актуальной задачи изучения дальнего космоса — забора проб грунта с Фо-
боса, спутника Марса [1, 2]. Рассмотрены сложные задачи оптимизации траекторий, для решения которых необходим синтез методов локальной и многоэкстремальной оптимизации, оптимального управления, космо-динамики, механики космического полета, небесной механики и численных методов.
Доставке проб грунта с Фобоса, различным аспектам этой проблемы и смежных задач марсианских экспедиций посвящено значительное количество работ, в частности [2, 3].
В первой версии проекта «Фобос — Грунт» планировалось оснастить КА комбинированной двигательной установкой, т. е. уход от Земли должен был осуществляться с помощью разгонного блока. Затем аппарат перестраивался и на протяжении перелета к Марсу управлялся электрореактивной двигательной установкой малой тяги, а маневры у Марса и возвращение обеспечивали двигатели большой тяги. При реализации проекта от перспективной схемы комбинирования двигателей большой и малой тяги отказались [4-7], однако проблема оптимизации комбинированных схем межорбитальных перелетов, в том числе построение траекторий таких экспедиций, осталась актуальной.
Проблематика оптимизации траекторий перелетов КА с комбинированным управлением двигателями большой и малой тяги широко известна [8, 9]. Эффективность использования электроракетной двигательной установки (ЭРДУ) в межпланетных экспедициях неоднократно отмечалась в различных источниках [3, 10-12]. В частности, указано, что в проекте «Фобос — Грунт» использование двигательной установки малой тяги в дополнение к двигательной установке большой тяги позволило увеличить массу полезной нагрузки на 150200 кг по сравнению с массой полезной нагрузки при перелете с применением только двигателя большой тяги [3].
В статьях [10-13] рассмотрена проблема оптимизации межпланетных перелетов КА с комбинированным управлением двигателями большой и малой тяги. При этом работа двигательной установки большой тяги либо аппроксимируется импульсным воздействием, либо считается, что КА имеет гиперболический избыток скорости при старте от Земли и при подлете к Марсу [10-12]. И.С. Григорьевым решается задача оптимизации плоского перелета между низкими круговыми орбитами искусственных спутников Земли (ИСЗ) и Марса [10], авторами статьи [12] — между орбитами Земли и Марса. В статьях [10, 12] высказано предположение, что Земля и Марс движутся по компланарным круговым орбитам. В источнике [11] в качестве орбит Земли и Марса взяты средние эллиптические орбиты этих планет, при этом притяжение Марса не учитывается.
В работе [13] рассмотрен полет КА, оснащенного двигательными установками большой и малой тяги. При этом максимизируется по-
лезная масса, доставляемая к астероиду или комете. Принято, что на первом участке траектории, внутри сферы действия Земли, работает только двигательная установка большой тяги, разгоняющая КА до заданного значения полной геоцентрической энергии /2. Кроме того, учтена зависимость характеристической скорости маневра от скорости ухода от Земли, которую нужно набрать на первом участке. После выключения двигательной установки большой тяги начинается второй участок, на котором работает двигательная установка малой тяги, а возмущения от Земли не учитываются. Расчет ведется последовательно: сначала для первого, а затем для второго участка. Масса КА после маневра определяется по формуле Циолковского. Выполняется массовый расчет для двигателя малой тяги на основе траектории перехода к цели при выбранной начальной скорости . Высказана надежда в том, что переход от идеально регулируемой малой тяги к кусочно-постоянной не приведет к существенным потерям массы полезной нагрузки.
Постановка задачи. В настоящей работе рассматривается задача оптимизации экспедиции пространственного перелета КА к Фобосу с возвращением к Земле.
Определен временной период перелета: с 2020 по 2030 гг. В первой части миссии (задача перелета «вперед») КА стартует с круговой орбиты ИСЗ, соответствующей выведению с космодрома «Байконур» при фиксированных угле наклона к экватору ¡в = 51,6° и высоте ИЕ = 200 км, а также при свободных долготе восходящего узла ^0 и положении КА на орбите ф0 .
Прилетев на Фобос, КА с помощью доставленной научной аппаратуры проводит исследования, в том числе забор проб грунта. Продолжительность пребывания КА на Фобосе составит не менее 30 сут. Во второй части миссии (задача перелета «назад») КА полетит в обратном направлении — от Фобоса к Земле.
Общая продолжительность экспедиции ограничена 1500 сут, положение Земли и Марса соответствует эфемеридам БЕ424, орбита Фобоса — МАВ.097 [14]. Гравитационные поля Солнца, Земли и Марса считаются центральными ньютоновскими. Гравитационное поле Земли учитывается только в первой части миссии. Поскольку в конечный момент времени КА тормозится об атмосферу Земли и их скорости не выравниваются. Можно предположить, что КА и Фобос представляют собой не притягивающие материальные точки, их координаты и скорости в конечный момент первой части и в начальный момент второй части миссии совпадают. Угловое положение КА на исходной стартовой орбите, моменты старта и финиша оптимизируются.
Начальная масса КА на круговой орбите ИСЗ М0 = 10 т. Предполагается, что КА последовательно управляется четырьмя разными дви-
гательными установками [2, 15]. Двигательная установка «Фрегат» с тягой Т = 19613 Н, удельной тягой Руд = 331 с и сухой массой МФ = 735 кг осуществляет разгон КА у Земли. Затем на протяжении перелета от Земли к Марсу может работать двигательная установка малой тяги СПД-230 с Т = 0,784532 Н вблизи Земли и Т = 0,392266 Н вблизи Марса, Руд = 2700 с, сухой массой ММТ = 350 кг.
Торможение у Марса обеспечивается двигательной установкой перелетного модуля (ДУПМ) с T = 1529,8374 Н, удельной тягой Руд = 304 с и сухой массой МдУПМ = 180 кг. Разгон у Марса при старте в обратном направлении (к Земле) выполняется двигательной установкой возвращаемого аппарата (ДУВА) с T = 125 Н, удельной тягой Руд = 306 с и сухой массой Мдува = 139 кг. В первой части миссии
каждая очередная ДУ после отработки сбрасывается. На Фобосе также остается перелетный модуль массой Мпм = 600 кг (550 кг — сам перелетный модуль и 50 кг — доставленная научная аппаратура для проведения дальнейших исследований).
Во второй части миссии двигательная установка возвращаемого аппарата сбрасывается в конечный момент времени, так как на протяжении перелета за ее счет можно проводить корректировочные включения тяги. Они есть в схеме миссии, представленной в работе [2], но здесь они не учитываются. Предполагается, что сбросы двигательной установки вместе с баками происходят мгновенно, управление КА осуществляется модулем и направлением вектора реактивной тяги, отлет от Земли и подлет к ней и отлет от Марса — безвитковые, реализующиеся за одно включение двигателей большой тяги. Моменты включения, выключения двигательных установок малой тяги и возвращаемого аппарата оптимизируются, а двигательных установок «Фрегат» и перелетного модуля — заданы продолжительностями их работы — параметрами рассматриваемой задачи.
Минимизируемым функционалом задачи является конечная масса:
M к = M 0
(T Ах ^
д.у д.у + M
С ^^д.у
^ min, (1)
д.у V "Д.у J
где Сд.у = gßPyn — скорость истечения реактивной струи соответствующей двигательной установки; Ах д.у — суммарное время работы соответствующей двигательной установки; gE = 9,81 м/с2 — ускорение свободного падения.
В рамках постановки этой задачи проведем сквозную оптимизацию траектории всей экспедиции по параметрам задачи, оценим возможный выигрыш при использовании комбинированной тяги по сравнению с использованием только двигателей большой тяги.
Формализация. В связи с потерей точности [16], связанной с тем, что при использовании какой-либо одной системы координат (СК) вычисление правых частей дифференциальных уравнений движения КА вблизи Земли или вблизи Марса имеет значительные погрешности вследствие вычитания близких значений координат КА и притягивающего центра. Движение КА на каждом из пяти участков траектории рассматривается в одной из трех различных СК (рис. 1).
Орбита Марса
Рис. 1. Схема перелета от Земли к Марсу
На схеме (см. рис. 1) первой части экспедиции зеленым цветом выделен участок использования геоцентрической СК, синим — гелиоцентрической и красным (при подлете к Марсу) — марсоцентри-ческой. Задача возвращения рассматривается сначала в марсоцентри-ческой, а затем в гелиоцентрической СК. Причем эти СК — невра-щающиеся декартовы, их оси параллельны осям СК 12000. Ниже во всех формулах г = 1.. .5 — номер участка траектории.
Отметим, что на 1-м, 3-м и 4-м участках СК являются неинер-циальными.
Данная задача космодинамики как задача оптимального управления формализуется совокупностью динамических систем, поскольку вид дифференциальных связей изменяется в зависимости от участков траектории вследствие разных СК, и решается на основе соответствующего принципа максимума [17].
На каждом из пяти участков траектории компоненты векторов положения и скорости КА обозначим хг, уг, и иг, , wi соответ-
ственно; компоненты векторов положения Земли и Марса в гелиоцентрической СК — xse , ySE, zse и xsm , ySM, zsm ; компоненты разности этих векторов — xME, yME, zME, а также введем другие °б°значения: xei = xi, yEi = Уь zEI = zi, xSi = xi + xSE , ySi= yi + ySE, zS1 = z + zSE , xMi = xi + xME , yMi= yi + yME , zMi= z + zME ; xEj = xj - xSE , yEj = yj - ySE, zEj = zj - zSE , xSj = xj , ySj = yj , zSj = zj , xMj = xj - xSM , yMj = yj - ySM, zMj = zj - zSM для j = 5; xEk = xk - xME, yEk = = yk - yME, zEk = zk - zME, xSk = xk + xSM, ySk = yk + ySM, zSk = = zk + zSM, xMk = xk, yMk = yk, zMk = zk для k = 3, 4, где j и k —
номера соответствующего участка экспедиции.
Тогда на каждом из пяти участков расстояния от КА до центра
Земли, Солнца и Марса можно записать так: rEi =*JxE, + yEi + zEi,
rSi xSi + ySi + zSi, rMi = *J xMi + yMi + zMi.
Системы дифференциальных уравнений управляемого движения центра масс КА на каждом из пяти участков траектории имеют следующий вид:
x, = ы,, u, = -gxl b^T + — cost,cos^i = уы,,
B rBi m,
< у, = v, , V, = - gyi - Xmb ^t + — sin t, cos л, = Уг, , (2)
B rBi Ш,
. . ^ ZBl P . . P
z, = w,, w, = - gzi B—T +-sin t[, = уш, ш , = ^ —,
B rBi Шi C,
где gxi, gyi, gzi — компоненты вектора ускорения Земли; gx3, gy3, gz3 и gx4, gy4, gz4 — компоненты вектора ускорения Марса в гелиоцентрической СК; gx2 = gy2 = gz2 = gx5 = gy5 = gz5 = 0; Ш (t) = m(t)/mg —
обезразмеренная и M(t) — абсолютная масса КА; P (t) — модуль вектора ускорения за счет реактивной тяги, изменяющийся в пределах от G до Pmax = Гд.у /Mg; и л, — углы, определяющие направление вектора тяги.
Первое слагаемое в выражениях изменения скоростей в системе уравнений (2) обусловлено неинерциальностью рассматриваемых СК, остальные — отражают влияние закона всемирного тяготения: на КА действуют силы притяжения со стороны Солнца, Земли и Марса, равные ^BM (t) / rBi (t), где B — космическое тело (E — Земля, S — Солн-
це, М— Марс), дв — гравитационный параметр тела В = уМВ, —11 2 —2
(у = 6,669 -10 Н • м кг — гравитационная постоянная, МВ — масса тела). Гравитационный параметр Земли равен дЕ = 3,986 013-105 км3/с2, Марса — дМ = 4,297780-104 км3/с2, Солнца — ^ = 1,32712440018 •1011км3/с2. Для задачи перелета «вперед» суммирование в системе уравнений (2) ведется по В е {Е, Б,М}, для задачи перелета «назад» — по В е {Б,М}.
Переход из первой во вторую СК происходит в момент времени 1е пересечения КА сферы влияния Земли, из второй в третью СК — в момент времени 1м и из четвертой в пятую СК — в момент времени ^м 2, в которые КА пересекает сферу влияния Марса. Поскольку на каждом из участков 1, 2, 3 или 4, 5 учитывается гравитационное воздействие одних и тех же притягивающих центров, система уравнений (2) при переходе от одного к другому качественно не изменяется, и на выбор моментов времени 1е , tм влияет только потеря точности, а значит, их можно зафиксировать.
До некоторого момента времени t = te движение КА описывается в СК, связанной с Землей. В этот же момент выполняется и переход к описанию движения КА в гелиоцентрической СК:
Г х2 = х1 + хж, У2 = У1 + Убе, ¿2 = ¿1 + ¿бе, 1 ^2 = Щ, (3)
[ П2 = Щ + Пе , = V + Уе , W2 = + We,
где {ще , Уе ,we} — компоненты вектора скорости Земли в гелиоцентрической СК.
Аналогичные формулы имеют место для связи гелиоцентрической и марсоцентрической СК. До некоторого момента времени t = tм движение КА описывается в СК, связанной с Солнцем, и в тот же момент происходит переход к описанию движения КА в марсо-центрической СК:
| х3 = х2 — хбм, у3 = у2 — убм, ¿3 = ¿2 — 2бм, ..ч
1 т = т2, (4)
[ Щ3 = Щ2 — Щм, У3 = У2 — Ум , Wз = W2 — Wм,
где {щМ,УМ, wм ] — компоненты вектора скорости Марса в
гелиоцентрической СК.
В задаче перелета «назад» движение КА сначала описывается в марсоцентрической СК. В момент времени t = tм2 осуществляется переход к описанию движения КА в гелиоцентрической СК:
I
х5 = х4 + хбм, у5 = у4 + убм, ¿5 = ¿4 + 2бм, Щ5 = Щ4 + щм , У5 = У4 + ум , W5 = W4 + Wм,
т5 = т4. (5)
Средний экваториальный радиус Земли обозначим Яе = 6378,16 км. В начальный момент t = to, до включения двигательной установки «Фрегат», КА находится на круговой орбите ИСЗ:
Представленные соотношения означают, что масса КА задана, расстояние от КА до центра Земли равно радиусу круговой орбиты, вектор скорости КА относительно Земли ортогонален радиусу-вектору КА, величина вектора скорости КА относительно Земли равна круговой скорости, наклон орбиты соответствует выведению КА с космодрома «Байконур».
В конечный момент времени первой части миссии t = ^, после отработки двигательной установки перелетного модуля, координаты и скорости КА должны совпасть с координатами и скоростями Фобоса хр,ур, гр и пр, Ур, Wp соответственно:
х3 = хР, у3 = ур, ¿3 = ¿р, щ3 = щр, у3 = Ур , wз = ^. (7)
В задаче возвращения при t = t2 координаты и скорости КА также должны совпадать с координатами и скоростями Фобоса:
х4 = хр, у4 = ур, ¿4 = ¿р, Щ4 = Щр, У4 = Ур, W4 = Wp. (8)
В конечный момент времени задачи возвращения при t = tз координаты КА должны совпасть с координатами центра Земли. Предполагается, что скорость КА гасится об атмосферу Земли:
Формально обозначим: tф—, ^.т—, Тдупм—, ^.м—, — — моменты
времени до сброса соответствующей двигательной установки или перелетного модуля (Ф — «Фрегат», м.т — малая тяга, в.а — возвращаемый аппарат); ^ +, ^.т+, ^УПМ +, ^.м+, + — моменты времени после
сброса. При этом предполагается, что tф+ = tф_ = tф, ^.т+ = tм.т_ = ,
^ДУПМ+ = ^ДУПМ— = tп.м+ = tп.м_ = tв.а + = ^.а — = t3. Тогда 8 Инженерный журнал: наука и инновации # 7 2017
т1 = 1, х2 + у12 + ¿12 = ( + Ье )2 , хщ + у1У1 + = 0,
(6)
х5 = хе, у5 = уе, ¿5 = ¿е.
(9)
т (/ф+ ) = т (tф_ )(1—Мф /М0), т3 (tм.т+ ) = т3 (/м.т— )(1—Мм.т /М0), 1 т3 (tДУПМ+ ) = т3 (tДУПМ_ )(1—Мдупм /М0), отл ((,■ + ) = ап ($} — ), (10)
т4(См+ ) = m4(tп.м_ )(1 — мп.м /М0), т5^в.а + ) = т^а — )(1—Мв.а /М0),
где условие ап +) = ап —) обозначает 30 условий вида х^ф+) = х^ф)
непрерывности всех фазовых переменных, кроме массы в момент времени сброса двигательной установки или перелетного модуля.
При решении задачи требуется определить неизвестные фазовые переменные xi (•), yi (•), (•), щ (•), Уi (•), wi (•), mi (•), моменты времени
t1, t2, t3, параметры управления р (•), ^ (•), ^ (•), удовлетворяющие системам дифференциальных уравнений (2), краевым условиям (6)-(10) и минимизирующие функционал (1).
Краевая задача. Применим к рассматриваемой задаче оптимального управления принцип максимума. Пусть существует оптимальный в сильном смысле процесс в задаче, описываемой (1)—(10), такой, что Зв> 0: гтЦ) >в Vt е[to, ^] и т51 (t) >е, м(t) >е, Vt е[to, tз], и выполняются условия гладкости: правые части систем дифференциальных уравнений (2) и их частные производные по хг, уг, , Щ, у г, wi, тг, i = 1...5 непрерывны в окрестности оптимальной траектории, функции, определяющие функционал, и ограничения (1), (6)-(10) непрерывно дифференцируемы на каждом отрезке [^, ^+1 ], к = 1...8,
где tk последовательно «пробегает» через токи /ф, tE, tм, /мт, tм2, Тогда согласно принципу максимума Понтрягина для задач управления совокупностью динамических систем [17] существуют множители Лагранжа: константы Xг0, Xгу0, Xу0, ХгБ0,
Ххр, Хур, Х¿р, Хир, Хур, Х wР, Х х1Е, Ху(Е, Х , ХШЕ, ХУЕ, Х wíE, ХхМ, Ху(М, Х¿/М, ХшМ, ХуМ, ХwíM, Ххр2, Хур2, Х¿р2, Хщр2, Хур2, XwР2, Х хЕ, Х уЕ, Х ¿Е, XxtM2, XytM2, XztM2, XutM2, XytM2, Х wíM2,
ХТа11, ХТр, Хт0, ХтЕ , ХтМ, ХтМ 2, 5 констант Xту и 30 констант
Х ап] для условий (10) и функции рх1, Руг, р21, рш , Руг, Pwi, Ртг , не равные одновременно нулю, такие, что на оптимальной траектории выполняется система необходимых условий оптимальности. Выпишем основные задачи конструкции принципа максимума — функции Понтрягина:
= РхгЩ + РуУг + Р^г + РшгФш + Р«Фу + РwiФwi — Ртгр / , и терминант:
I = х(0Мк + ^0 + Ье +1М + 1Л +^2 + 4м2 ++1т + 1а ,
где для задачи перелета «вперед»
1г 0 = Кбо
ху + ухщ -(яЕ + ъЕ ))—соб гБ +^„0 (+ УУ + С '\яе + ъЕ )
++ Кго ^щ2 + у? + (2 - ЯВ^+Ъ^) + ХгО (х2 + У1 + - ( + Ъе)2) +
+Кто(то -1) при г = Го;
Ъе = Кв (х2 - XI - хж ) +Ке (У2 - У1 - Узв ) + Кв ( - ¿1 -¿ЪЕ ) + + КгЕ (и? -щ -ив ) + Кв (У2 -У1 -Ув ) + К„гв ((2 -WE ) + Ктв (Ш2 -Ш1 )
при г = гв;
1гм = Ххгм (х3 - х2 + ХЪМ ) + Хугм (уз — У2 +Уъм ) + Хггм ( ¿3 - ¿2 + ¿ъм ) +
+ Км (из - и2 +им) + Км (Уз-У2 + Ум ) + Км (( - (2 + (м ) + +Хтм (т2 - т1) при г = гм;
11 = КхЕ (х3 - хЕ ) +К (Уз - Уе ) + К (¿3 - ) + КиЕ (и3 - ) +
+ К (Уз - ) + К ((з - () при г = Ч; для задачи перелета «назад»
К2 = КхЕ2 (х4 - хЕ ) + К2 (У4 - Уе ) + К2 (¿4 - ¿Е ) + КиЕ2 (Щ4 - ) + +К 2 (У4 - ) + К 2 (( - ) пРи г = г2; 1м2 = Ххгм2 (х5 — х4 + хъм ) +
+ Км2 (У5 - У4 + Ум ) +Км2 (¿5 - ¿4 + ) +Км 2 ( - и4 + Щм ) + + Км 2 (У5 - У4 + Ум ) + Км 2 ((5 - (4 + (м )+Ктм 2 (( - т4 ) ^ Г = Гм 2; Ьз = КхЕ (х5 - хЕ ) + Ке (У5 - Ув ) + КЕ (¿5 - 2в ) при Г =
для сброса двигательной установки или перелетного модуля
1т = 2 ] Кщ С тк (г]+ ) - тк (г]— ) + м); 4 = Е} Е п Кап] (а„ (] ) - ап (г]— )), V м О У
где индекс ] последовательно «пробегает» набор {Ф, м.т, ДУПМ, п.м, в.а}; индексу к при этом соответствуют числа {1, з, з, 4, 5}; ап для каждого фиксированного к последовательно пробегает хк,
Ук, ¿к, ик, ук, (к .
Применение принципа максимума сводит решение исходной задачи оптимального управления к решению краевой задачи, состоящей из исходных уравнений (2)-(10) и следствий условий оптимальности, а именно: к системе дифференциальных уравнений (2) добавляется сопряженная система (11):
Рм - 2 (ри ( - 3хВг ) - 3 (ру1хВ12В1 + р^1хВ1уВ1
г, Гг,, \ 4 '
в V ГВ1 г
Л
Рyi - 2 ^Т"( ( - 3У2В1) - 3 (РшУв^Вг + РшУВ1хВ1))
в V ГВ1 х
Р* - 2 ^г( ( - 3*Вг) - 3 (в.ув. + ру1*в1хв1
\'В1
• -- • -- • -- • -рш рхг, ру. ру., р т , ротг 2
т
(11)
где для задачи перелета «вперед» суммирование ведется по В = {Е, 5,М}, «назад» — по В = {5, М}, р. ^ •
Из условий трансверсальности и стационарности к краевым условиям (6)-(8), (10) добавляются следующие соотношения: при t -10
рх1 - 2Хг0х1 + Хгу0и1 + хж0у1, ри1 - Хгу0х1 + 2Ху0и1 - хж0у1-.
< ру1 - 2Хг0У1 + Хгу0у1 - хж0иь ру1 - Хгу0У1 + 2Ху0у1 + хж0x1, (12) Р^1 - 2Хг0гу0^ъ Р^1 - Хгу0+ 2Ху0^ъ Рт1 - 0, н1 - 0, при t - tl
Н3 - -Рх3иЕ - Ру3уЕ - - РиЗЯх¥ - Руз^уЕ - Р^3, (13)
где { gxF, ^уЕ, gzF) — компоненты вектора ускорения Фобоса в
марсоцентрической СК; при t -12
Н4 - -Рх4и¥ - Ру4 - Р*4^ - Ри4gxF - ру4gyF - Р^4gzF; (14) при t - tз
ри5 - 0, Ру5 - 0 Рм>5 - 0, Н5 - -рх5иЕ - Ру5УЕ - р*5™Е• (15)
Кроме того, в задачу добавляются условия дополняющей нежесткости и неотрицательности:
\Taill (3 - to -1500) = 0, XTF ( - h - 30) = 0,
(16)
^0 > 0, ^Taill > 0, Xtf > 0,
а также условия непрерывности сопряженных переменных в моменты сброса двигательной установки или перелетного модуля и формулы их пересчета при переходах между разными СК:
Pi+i(tk) = ^TL(tk )Pt (tk), (17)
ох,
где tk — момент склейки участков i и i +1.
В фиксированный момент времени, когда КА «сидит» на Фобосе, перестает учитываться притяжение Земли. При этом функция Понт-рягина терпит разрыв, фазовые и сопряженные переменные остаются непрерывными.
Л mtPmt
Функция переключения управления имеет вид %t = pt--——,
направление тяги определяется соотношениями
cosacos ^ = PuL, cos £t sin ^ = PvL, sin ^ = EwL. (18)
pt pt pt
При xi > 0 тяга двигателя максимальна, при хi < 0 — выключена.
Использование для фазовых переменных формул пересчета (3) в момент времени t = tE, формул (4) в момент времени t = tM , формул (5) при t = tm2 и системы уравнения (10) в моменты сброса двигательной установки или перелетного модуля, а также условий непрерывности при t e{to,¿м.т, t1, t2, ¿з } и склейки (17) при t e{tE,tM, t1, t2, tM2} для сопряженных переменных позволяет свести решение совокупности шеститочечной краевой задачи перелета «вперед» и трехточечной задачи перелета «назад» к решению двухточечных задач.
Вычислительная схема. Полученную нелинейную краевую задачу принципа максимума 70-го порядка (2)—(18) можно решить только численно. В настоящей работе задача решена по вычислительной схеме с использованием метода стрельбы [18, гл. 2] и модифицированного метода Ньютона [19, 20]. Задачи Коши интегрировались методом Рунге — Кутты порядка 8(7) с автоматическим выбором шага [21]. Системы линейных уравнений решались методом Гаусса с выбором главного элемента [22] и повторным пересчетом. Производная вектора-функции невязок по параметрам пристрелки вычислялась численно с помощью центральных разностей [16].
Для решения задачи методом Ньютона следует подобрать близкое к решению локальное начальное приближение. Поскольку исход-
ная задача многоэкстремальна, для поиска области, в которой следует искать локально оптимальное решение, необходимо предварительно провести глобальную оптимизацию. Для этого сначала задача исследовалась на основе решения серии задач Ламберта с учетом одного притягивающего центра на каждом участке перелета [23, 24], а затем решалась в импульсной постановке на основе принципа Лагранжа. В результате проведенного анализа для импульсного случая (рис. 2) получена окрестность экстремалей с лучшим значением функционала, в которой далее уже можно решать задачу с ограниченной тягой [25]. К тому же из рисунка видно, что окна старта к Марсу есть каждые два года, время старта ^ отсчитывается от 1 января 2020 г., а лучший функционал оказался у экспедиции со стартом в 2026 г.
Рис. 2. Окна старта к Марсу в период 2020-2030 гг.: Р(*о) — сумма модулей импульса у Земли и двух импульсов у Марса в импульсной задаче оптимизации экспедиции; /0 — время старта
При исследовании задачи в упрощенной постановке выявлено, что КА на Фобосе выгодно «сидеть» около одного земного года. А значит, в условии (16) обращается в нуль, и исходную задачу можно в окрестности данного предполагаемого локального минимума легко разбить на две независимые: задачу «вперед» и задачу «назад», причем совокупность их решений будет экстремалью в общей задаче.
После получения начального приближения траектории строились следующим образом. Задачу перелета «вперед» просчитывали с двух концов (^ и в середину траектории ^ (рис. 3), затем задачу перелета «назад» считали от t2 к tз. В качестве параметров пристрелки выбирались 12-14 значений: Xг0, ^0, Ф0, г^), Гз(^),
Рх3(к), РУ3(А), Р*3(^), хк, где хк — моменты включения и отключения малой тяги в задаче «вперед», при исследовании задачи на траекториях реализовывалось управление с двумя или тремя активными участками малой тяги, тогда к - 2 или к =4 соответственно. И семь параметров пристрелки: Ахдува, ^ и рх4, ру4, р*4, £4, г4 при
t -12 в задаче «назад». По данным величинам с учетом эфемерид можно найти значения всех фазовых и сопряженных переменных, кроме рт1, в моменты времени to, t2.
Рис. 3. Вычислительная схема. Просчет с двух концов в середину траектории
Невязками выступают 12 условий склейки фазовых и сопряженных переменных в момент времени {', 7 краевых условий на правом конце Ц и, если необходимо, 2 условия хг(х2) - Хг(х3) - 0. Значения рт для интегрирования определялись из условий хг (х1) - Хг (хк) --Х4(АхдУВА) - 0. В качестве нормировки принималось условие непрерывности рт в точке t.
При построении экстремалей Понтрягина время работы двигательной установки большой тяги Лхф и Дхдупм и моменты времени
к, t2 считались заданными, далее проводилась внешняя оптимизация задачи по этим величинам с использованием градиентных методов и продолжения решения по параметру. При заданных значениях времени Лхф , Дхдупм, Ь, ^, t2 и выбранных значениях парамет-
ров пристрелки Тк, Атдува , гз массу следует вычислять аналитически по системам (2), (6), (10) и исключить из неизвестных.
В промежуточных расчетах при продолжении решения по параметру положение Фобоса на орбите вокруг Марса в момент времени г1 не фиксировалось из эфемерид, а рассматривалось как параметр. Для существенного ускорения счета при продолжении решения по параметрам задачи использовалась линейная и параболическая аппроксимация параметров пристрелки.
Результаты расчетов. Предложенная методика позволяет строить экстремали Понтрягина в исходной задаче. Численные расчеты проводились для большого набора значений параметров, в том числе при различных двигательных установках малой тяги, ммт, мо, м(г2), ДтФ, Дтдупм, г0, г1, г2, с помощью метода продолжения
решения по параметру.
Приведем пример конкретной экстремали. Для представленной на рис. 4 проекции траектории экспедиции на плоскость эклиптики суммарное время работы двигательной установки малой тяги составляет 1з0,85 сут. Стартует КА в момент времени Го (2з октября 2026 г. 0з:56:06), при этом х1 = 879,17 004 851 км, У1 =-4809,14 287 з2з км, 21 = - 4401,294 56з 67 км, и1 = 6,91407 517 км/с, у =-1,64 з04 414 км/с, (1 = з, 17 640 222 км/с, рх1 =-7,94 50з0з8, ру1 =-52,07 78з 556, рг1 =-60,02 088 770, ри1 = 0,84з20 9з6, ру1 =-0,05 977 з25, р(1 = 0,53 425195, Кг0 = 8774,8з082 009, =-1,02з 14187,
ф0 = 5,68775з45, £ =-0,07 076941, % = 0,56 з62 257.
Двигательная установка «Фрегат» работает 1088,04 541 645 с, реализуя за это время избыток скорости з, 477 з67 км/с. В момент
времени г1 (1 сентября 2027 г., 16:29:02) хз =-4724,89 72з 949 км,
уз = 5692,78 565 з09 км, ¿з = 55з9,12 7з0912 км, из =-1,56 9з7 954 км/с,
уз =-1,48 650 262 км/с, (з = 0,17 102 з65 км/с, рхз = 8,з8зз6 441,
руз =-8,98 2з7 741, рг3 =-5,17 004 776, риз = 0,78211515,
ру3 = 0,62007 751, р(3 = 0,06164 224, £ з = 0,67 0з48з0,
% = 0,06168135.
Двигательная установка перелетного модуля работает 1502,36195 216 с, реализуя за это время избыток скорости 1,352682 км/с. Невязки в задаче «вперед» вычислялись в момент времени г * = г0 +150 сут.
Рис. 4. Проекция траектории экспедиции на ОХУ плоскость эклиптики:
¡о, ¿ъ Ь — области работы двигательной установки большой тяги, 1 — работа двигательной установки малой тяги; 2 — пассивные участки; 3 — орбита Земли; 4 — орбита Марса
В обратном направлении КА улетает с Фобоса в момент времени t2 (3 сентября 2028 г. 18:49:31). При этом х4 --7438,14156 580 км, у4 --5620,52029 090 км, г4 -1097,27 061370 км, и4 - 0,93517617 км/с, у4 - -1,43 311 628 км/с, - -1,27 677 228 км/с, рх4 --12,38 649 475, ру4 --3,75 652 430, рг4 - 4,29 612 288, ри4 - 0,63 345 872, Ру4 - -0,49 066 283, рк4 - -0,59 831 432 и прилетает к Земле в момент ¿з (9 августа 2029 г. 03:51:04). Двигательная установка в работает 6213,57162 695 с, реализуя избыток скорости 1,848 938 км/с. Моменты времени перехода между различными СК: tЕ - ^ + 3,12124 399 сут, гм - ^ -хм, ¿м2 - t2 +хм, где хм - 2,71974141 сут.
Основные расчетные единицы: а.е./100 -1495 978,70 691 км и Б.Б. - 86 400 с. При использовании других величин приведенные значения сопряженных переменных необходимо пересчитать.
График функции переключения на участке работы двигательной установки малой тяги представлен на рис. 5. На данной траектории имеем один пассивный и два активных участка ее работы. Около Земли двигательная установка малой тяги работает 46,54 311 508 сут, вблизи Марса — 84,305 953 28 сут.
Рис. 5. Функция переключения управления (время считается в днях в системе наблюдения J2000)
График изменения массы представлен на рис. 6.
В таблице для сравнения приведены значения массы КА после посадки на Фобос, перед отлетом с Фобоса и доставленной к Земле в конечный момент времени после сброса двигательной установки возвращаемого аппарата в различных вариантах миссии. Выигрыш по массе при использовании большой и малой тяги на данной экстремали при посадке на Фобос составляет 128 кг, или 12,3 %, при отлете от Фобоса — 29,3 %, при подлете к Земле — 66,58 кг, или 67,7 %, что позволяет сделать вывод о целесообразности использования в миссии комбинации двигательных установок большой и малой тяги.
m{t) 1,0
0,8 0,6 0,4 0,2 0
»Ф
Ti
мт
'1
1 г2+Лт
дува
9800 10000 10200 10400 10600 г
Рис. 6. Изменение массы в зависимости от времени:
1 — уменьшение массы за счет работы двигательных установок большой тяги: «Фрегата» от ^ до ГФ, перелетного модуля от Гм.т до Гь возвращаемого аппарата от Г2 до Г2 + Лтд,
2 —
участки работы двигательных установок малой тяги от ГФ до т и от т2 до Гм.т; 3 — уменьшение массы в моменты времени ГФ, Гм.т, Г2, Н за счет сброса двигательных установок и перелетного модуля; 4 — пассивные участки траектории перелета КА
Масса, доставляемая к Фобосу и обратно к Земле, в различных вариантах экспедиции
Обозначение Масса, кг, при тяге
большой и малой, работавшей 130,85 сут большой большой и малой, работавшей 0 сут
М(Ь) 1162,69 655 986 1035,04 885 519 858,50 068 947
U(h) 562,69 655 986 435,04 885 519 258,50 068 947
М(Гз) 164,86 977 357 98,28 969 192 3,53 554 911
График окна старта представлен на рис. 7, на нем видно, что в окрестности момента времени to есть примерно месяц для старта, когда масса, доставленная на Фобос, будет меньше массы на представленной экстремали не более чем на 5 %.
Рис. 7. Окно старта:
по оси ординат — абсолютное отличие ДГ даты старта от Г0; по оси абсцисс — уменьшение конечной массы на такой экстремали
Характерный вид выигрыша в результате применения комбинированной тяги при полете к Марсу по сравнению с использованием только двигателей большой тяги представлен на рис. 8. При построении этого графика задача была исследована в упрощенной постановке — два продолжительных участка работы двигательной установки большой тяги аппроксимировались импульсными воздействиями, направление которых оптимизировалось, а значение было задано и изменялось пропорциональными шагами в пределах от Ае = 3,65 и Ам = 2,00 км/с, что соответствует решению без использования двигательной установки малой тяги до значений 3,22 и 0,87 км/с, минимально необходимых для безвиткового ухода от Земли и подлета к Марсу. Лучшей траектории без малой тяги соответствует точка А. Видно, что в первый момент времени, когда двигательная установка малой тяги установлена на КА, но не работает, происходит резкий проигрыш по функционалу из-за ее массы, а с уменьшением АЕ и Ам продолжительность работы установки постепенно увеличивается и в итоге получается существенный выигрыш массы.
МЩ, кг
1500
1700
1900
1300
ае,% 99 97 95 93 91 89
Ал/,% 95 85 75 65 55 45
Рис. 8. Изменение конечной массы при увеличении времени работы двигательной уста-
новки малой тяги:
по осям абсцисс — прирост значений импульсов у Земли и у Марса де и ДМ по сравнению с теми же параметрами у КА без двигателей малой тяги; по оси ординат — масса, доставленная КА
на Фобос
Отметим характерное «насыщение» на графике — в определенный момент кривая массы достигает максимального значения, а соотношение времени работы двигательной установки малой тяги и значений импульсов становится локально оптимальным.
Заключение. Разработанная методика вносит существенный вклад в теорию решения сложных задач траекторной оптимизации. Она позволяет проводить сквозную оптимизацию экспедиции КА с построением единой экстремали Понтрягина с учетом эфемерид и притяжения планет, с комбинированным управлением КА. В рассматриваемых задачах требуется синтез методов локальной и многоэкстремальной оптимизации. Представленная методика построения начального приближения дает возможность проведения глобальной оптимизации с помощью последовательной аппроксимации и исследования миссии на основе серии задач Ламберта и принципа Лан-гранжа.
На основе принципа максимума Понтрягина для управления совокупностью динамических систем в данной работе проанализирована актуальная задача оптимизации миссии к Фобосу с возвращением к Земле. Представленные результаты численных расчетов девятиточечной краевой задачи 70-го порядка позволяют рассматривать целесообразность использования комбинированного управления КА с помощью двигательных установок большой и малой тяги.
[1] Энеев Т.М. Актуальные задачи исследования дальнего космоса.
Космические исследования, 2005, т. 43, № 6, с. 403-407.
ЛИТЕРАТУРА
[2] Фобос-Грунт. Проект космической экспедиции. В 2 т. Т. 1, 2. Москва, ФГУП «НПО им. С. А. Лавочкина», 2011, 519 с.
[3] Авдуевский В.С., Аким Э.Л., Кремнев Р.С., Куликов С.Д. Космический проект «Фобос-Грунт»: основные характеристики и стратегия развития. Космонавтика и ракетостроение, 2000, т. 19, с. 8-21.
[4] Аким Э.Л., Заславский Г.С., Морской И.М., Рузский Е.Г. Доставка на Землю реликтового вещества с Фобоса — проект «Фобос-Грунт»: баллистика, навигация и управление полетом. Астрономический вестник, 2010, т. 44, № 1, с. 29-40.
[5] Маров М.Я. Космический проект «Фобос-Грунт»: новый этап российской планетной программы. Астрономический вестник, 2010, т. 44, № 1, с. 3-б.
[6] Галимов Э.М. Научное обоснование проекта доставки грунта с Фобоса. Астрономический вестник, 2010, т. 44, № 1, с. 7-1б.
[7] Симонов А.В., Морской И.М., Степаньянц В.А., Тучин А.Г. Баллистическая схема полета КА «Фобос-Грунт». Вестник ФГУП НПО им. С.А. Лавочкина, 2011, № 3, с. бб-73.
[8] Edelbaum T.N. The Use of High- and Low-Thrust Propulsion in Combination for Space Missions. Journal of the Astronautical Sciences, 19б2, vol. 9, no. 2, pp. 49-б9.
[9] Иванов Ю.Н. Оптимальное сочетание двигательных систем. Изв. АН СССР. Механика и машиностроение, 19б4, № 2, с. 3-14.
[10] Григорьев И.С. Оптимизация траекторий межпланетных перелетов космического аппарата, управляемого посредством двигателей большой и малой тяги. Фундаментальные и прикладные проблемы космонавтики, 2000, № 3, с. 11-18.
[11] Федотов Г.Г. Об использовании возможностей комбинации большой и малой тяги при полетах к Марсу. Космические исследования, 2001, т. 39, № б, с. б13-б21.
[12] Bolle A., Circi C., Corrao G. Optimal Mars transfers for small payload transportation. Celest. Mech. Dyn. Astron. 10б (200б), pp. 183-19б.
DOI: 10.1007/s10569-009-9250-1 (дата обращения 30.12.201б).
[13] Ахметшин Р.З., Белоглазов С.С., Белоусова Н.С., Глазков А.И. Оптимизация перелетов к астероидам и кометам КА с комбинированием двигателей большой и малой тяги. Препринт ИПМ им. М.В. Келдыша, 1985. № 144.
[14] An Observation Geometry System for Space Science Missions. Navigation and Ancillary Information Facility of JPL of NASA. NASA, 2017, 13 Apr. URL: http://naif.jpl.nasa.gov/naif (дата обращения 30.04.2017).
[15] Федеральное государственное унитарное предприятие "Опытное Конструкторское Бюро «Факел»". Характеристики стационарных плазменных двигателей. URL: http://www.fakel-russia.com/ (дата обращения 30.12.201б).
[16] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. Москва, Изд-во БИНОМ. Лаборатория знаний, 2008, б3б с.
[17] Григорьев И.С., Григорьев К.Г. Об условиях принципа максимума в задачах оптимального управления совокупностью динамических систем и их применении к решению задач оптимального управления движением космических аппаратов. Космические исследования, 2003, т. 41, № 3, с. 307-331.
[18] Григорьев И.С. Методическое пособие по численным методам решения краевых задач принципа максимума в задачах оптимального управления. Москва, Изд-во Центра прикладных исследований при механико-математическом факультете МГУ, 2005, 159 с.
[19] Федоренко Р.П. Введение в вычислительную физику. Москва, Изд-во Моск. физ.-техн. ин-та, 1994, 526 с.
[20] Исаев В.К., Сонин В.В. Об одной модификации метода Ньютона численного решения краевых задач. Журнал вычисл. матем. и матем. физ., 1963, т. 3, № 6, с. 1114-1116.
[21] Хайрер Э., Нёрсетт С.П., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Москва, Мир, 1989, 512 с.
[22] Богачёв К.Ю. Практикум на ЭВМ. Методы решения линейных систем и нахождения собственных значений. Москва, Изд-во механико-математического факультета МГУ, 1999, 137 с.
[23] Суханов А.А. Астродинамика. Сер. Механика, управление, информатика. Москва, Ротапринт ИКИ РАН, 2000, 202 с.
[24] Самохин А. С. Оптимизация экспедиции КА к Фобосу при управлении импульсными воздействиями с использованием решения задач Ламберта и учетом притяжения Земли и Марса. Вестник Московского университета. Сер. 1. Математика. Механика, 2014, № 2, с. 62-66.
[25] Григорьев И.С., Григорьев К.Г. Об использовании решений задач оптимизации траекторий КА импульсной постановки при решении задач оптимального управления траекториями КА с реактивным двигателем ограниченной тяги. II. Космические исследования, 2007, т. 45, № 6, с. 1-11.
Статья поступила в редакцию 15.03.2017
Ссылку на эту статью просим оформлять следующим образом: Григорьев И.С., Заплетин М.П., Самохин А.С., Самохина М.А. Оптимизация экспедиции к Фобосу с комбинированной тягой с возвращением к Земле. Инженерный журнал: наука и инновации, 2017, вып. 7. http://dx.doi.org/10.18698/2308-6033-2017-7-1639
Статья подготовлена по материалам доклада, представленного на XLI Академических чтениях по космонавтике, посвященных памяти академика С.П. Королёва и других выдающихся отечественных ученых — пионеров освоения космического пространства.
Москва, МГТУ им. Н.Э. Баумана, 24-27января 2017 г.
Григорьев Илья Сергеевич — канд. физ.-мат. наук, доцент кафедры «Вычислительная математика» механико-математического факультета МГУ им. М.В. Ломоносова. e-mail: iliagri@mech.math.msu.su
Заплетин Максим Петрович — канд. физ.-мат. наук, доцент МГУ им. М.В. Ломоносова, профессор РУДН. e-mail: Zapletin_m@mail.ru
Самохин Александр Сергеевич — ассистент кафедры «Вычислительная математика» механико-математического факультета МГУ им. М.В. Ломоносова, старший преподаватель РУДН, преподаватель НИУ ВШЭ, ведущий специалист департамента баллистического проектирования космических систем в Научно-техническом агентстве «Космоэкспорт». e-mail: SamokhinAlexander@yandex.ru
Самохина Марина Александровна — старший преподаватель РУДН, преподаватель НИУ ВШЭ, ведущий специалист департамента баллистического проектирования космических систем в Научно-техническом агентстве «Космоэкспорт». e-mail: kipt35@gmail.com
Optimization of Phobos mission with hybrid propulsion returning to the Earth
© IS. Grigoryev1, M.P. Zapletin1,2, A.S. Samokhin1'2'3, M.A. Samokhina2,3
1Lomonosov Moscow State University, Moscow, 119991, Russia 2Peoples' Friendship University of Russia, Moscow, 117198, Russia 3Higher School of Economics, Moscow, 101000, Russia
The article considers the problem of optimizing the spacecraft recovery space flight to Phobos. On the one hand, this task is confined to a real mission to Phobos, which the Russian Federation plans to accomplish in the coming years. On the other hand, the development of techniques for optimizing the spacecraft interplanetary trajectories is the question ofpresent interest. Usually, in such tasks, the circumplanetary legs of trajectory are neglected and there is no end-to-end optimization of the entire mission. The paper shows a technique for constructing Pontryagin extremals in similar problems, taking these features into account, using the example of an expedition to Phobos. The positions of the Earth, Mars and Phobos correspond to the ephemerides DE424 and MAR097. Gravitational fields of the Sun and Mars are considered to be central Newtonian. Spacecraft is equipped with high- and low-thrust engines. The problem of cosmodynamics is formalized as an optimal control task and then it is solved numerically by shooting method. The paper shows the particular constructed trajectories, estimates the possible gain due to using hybrid propulsion instead of only high-thrust engines.
Keywords: low-thrust, hybrid propulsion, end-to-end optimization, interplanetary flight optimization, spacecraft trajectory optimization, Phobos, flight to Phobos, Phobos mission, Mars mission, cosmodynamics problem
REFERENCES
[1] Eneev T.M. Kosmicheskie issledovaniya — Cosmic Research, 2005, vol. 43, no. 6, pp. 403-407.
[2] Fobos-Grunt. Proekt kosmicheskoy ekspeditsii. V 2 tomakh. Tom 1, 2. [Phobos-Grunt. The space mission project. Scientific publication in 2 volumes. Vol. 1, 2]. Moscow, FSUE "Lavochkin Association", 2011, 519 p.
[3] Avduevskiy V.S., Akim E.L., Kremnev R.S., Kulikov S.D. Kosmonavtika i raketostroenie — Cosmonautics and Rocket Engineering, 2000, vol. 19, pp. 8-21.
[4] Akim E.L., Zaslavskiy G.S., Morskoy I.M., Ruzskiy E.G. Astronomicheskiy vestnik — Solar System Research, 2010, vol. 44, no. 1, pp. 29-40.
[5] Marov M.Ya. Astronomicheskiy vestnik — Solar System Research, 2010, vol. 44, no. 1, pp. 3-6.
[6] Galimov E.M. Astronomicheskiy vestnik — Solar System Research, 2010, vol. 44, no. 1, pp. 7-16.
[7] Simonov A.V., Morskoy I.M., Stepanyants V.A., Tuchin A.G. Vestnik FGUP NPO im. SA. Lavochkina — Bulletin of Lavochkin Association, 2011, no. 3, pp. 66-73.
[8] Edelbaum T.N. Journal of the Astronautical Sciences, 1962, vol. 9, no. 2, pp. 49-69.
[9] Ivanov Yu.N. Izvestiya AN SSSR. Mekhanika i mashinostroenie — Proceedings of the AS USSR. Mechanics and Mechanical Engineering, 1964, no. 2, pp. 3-14.
[10] Grigoryev I.S. Fundamentalnye i prikladnye problemy kosmonavtiki — Fundamental and Applied Problems of Cosmonautics, 2000, no. 3, pp. 11-18.
[11] Fedotov G.G. Kosmicheskie issledovaniya — Cosmic Research, 2001, vol. 39, no. 6, pp. 613-621.
[12] Bolle A., Circi C., Corrao G. Celestial Mechanics and Dynamical Astronomy, 2006, no. 106, pp. 183-196. DOI: 10.1007/s10569-009-9250-1 (accessed December 30, 2016).
[13] Akhmetshin R.Z., Beloglazov S.S., Belousova N.S., Glazkov A.I. Optimizatsiya pereletov k asteroidam i kometam KA s kombinirovaniem dvigateley bolshoy i maloy tyagi [Optimization of flights to asteroids and comets of spacecraft with the combination of high- and low-thrust engines]. Keldysh Institute Preprint, 1985, no. 144.
[14] An Observation Geometry System for Space Science Missions. Navigation and Ancillary Information Facility of JPL of NASA. NASA, 2017, 13 Apr. URL: http://naif.jpl.nasa.gov/naif (accessed April 04, 2017).
[15] Federal gosudarstvennoe unitarnoe predpriyatie "Opytnoe konstruktorskoe Buro "Fakel". Kharakteristiki statsionarnykh plazmennykh dvigateley [Federal State Unitary Enterprise "Experimental Design Bureau "Fakel". Characteristics of stationary plasma thrusters]. Available at: http://www.fakel-russia.com/ (accessed December 30, 2016).
[16] Bakhvalov N.S., Zhidkov N.P., Kobelkov G.M. Chislennye metody [Numerical methods]. Moscow, BINOM. Knowledge Laboratory Publ., 2008, 636 p.
[17] Grigoryev I.S., Grigoryev K.G. Kosmicheskie issledovaniya — Cosmic Research, 2003, vol. 41, no. 3, pp. 307-331.
[18] Grigoryev I.S. Metodicheskoe posobie po chislennym metodam resheniya kraevykh zadach printsipa maksimuma v zadachakh optimalnogo upravleniya [Methodical manual on numerical methods for solving boundary value problems of the maximum principle in optimal control problems]. Moscow, the Center for Applied Research at the Faculty of Mechanics and Mathematics of Moscow State University Publ., 2005, 159 p.
[19] Fedorenko R.P. Vvedenie v vychislitelnuyu fiziku [Introduction to Computational Physics]. Moscow, MIPT Publ., 1994, 526 p.
[20] Isaev V.K., Sonin V.V. Zhurnal vychislitelnoy matematiki i matematicheskoy fiziki — Journal of Computational Mathematics and Mathematical Physics, 1963, vol. 3, no. 6, pp. 1114-1116.
[21] Hairer E., Nersett S.P., Wanner G. Solving Ordinary Differential Equations. Springer-Verlag Publ., 1987 [In Russ.: Hairer E., Nersett S.P., Wanner G. Reshenie obyknovennykh differentsialnykh uravneniy. Moscow, Mir Publ., 1989, 512 p.].
[22] Bogachev K.Yu. Praktikum na EVM. Metody resheniya lineynykh system i nakhozhdeniya sobstvennykh znacheniy [Workshop on a computer. Methods for solving linear systems and finding eigenvalues]. Moscow, Faculty of Mechanics and Mathematics of Moscow State University Publ., 1999, 137 p.
[23] Sukhanov A. A. Astrodinamika [Astrodynamics]. Seriya "Mekhanika, upravlenie, informatika" [Series "Mechanics, Control, Informatics"]. Moscow, Rotaprint. IKI RAS Publ., 2000, 202 p.
[24] Samokhin A.S. Vestnik Moskovskogo universiteta. Seriya 1. Mekhanika. Matematika — Moscow University Mechanics Bulletin, 2014, no. 2, pp. 62-66.
[25] Grigoryev I.S., Grigoryev K.G. Kosmicheskie issledovaniya — Cosmic Research, 2007, vol. 45, no. 6, pp. 1-11.
Grigoryev I.S. (b. 1969), Cand. Sc. (Phys.-Mat.), Associate Professor, Department of Computational Mathematics, Faculty of Mechanics and Mathematics, Lomonosov Moscow State University. e-mail: iliagri@mech.math.msu.su
Zapletin M.P. (b. 1966) graduated from the Faculty of Mechanics and Mathematics of Lomonosov Moscow State University. Cand. Sc. (Phys.-Mat), Associate Professor, Lomonosov Moscow State University, Professor, Peoples' Friendship University of Russia. e-mail: Zapletin_m@mail.ru
Samokhin A.S. (b. 1989) graduated from the Faculty of Mechanics and Mathematics of Lomonosov Moscow State University. Assistant Lecturer, Department of Computational Mathematics, Faculty of Mechanics and Mathematics, Lomonosov Moscow State University, Senior Lecturer, Peoples' Friendship University of Russia, Lecturer, Higher School of Economics, Leading Specialist, Department of Ballistic Design of Space Systems, Cosmoexport Aerospace Research Agency. e-mail: SamokhinAlexander@yandex.ru
Samokhina M.A. (b. 1987) graduated from the Faculty of Mechanics and Mathematics of Lomonosov Moscow State University, Senior Lecturer, Peoples' Friendship University of Russia, Lecturer, Higher School of Economics, Leading Specialist, Department of Ballistic Design of Space Systems, Cosmoexport Aerospace Research Agency. e-mail: kipt35@gmail.com
