Optimum trajectories for an Earth–asteroid–Earth mission with a high thrust flight

Full Text

Экспедиции к астероидам и кометам очень актуальны для космонавтики. Ещё К.Э. Циолковский отмечал это [1]. Уже реализованы два проекта полётов к малым телам Солнечной системы с возвратом к Земле образцов их вещества: “Stardust” (США), “Hayabusa” (Япония) [2]. Сейчас реализуются ещё два таких проекта — “Hayabusa‑2” (Япония), “OSIRIS-REx” (США) [3]. У нас в стране НПО им. С.А. Лавочкина разрабатывает проект полёта к опасному астероиду Апофис [4]. В рамках проблемы построения и анализа траекторий экспедиции к астероиду Апофис с возвратом к Земле нами исследованы две задачи. Это, во‑первых, задача построения энергетически оптимальных межпланетных траекторий полёта от Земли к астероиду и траекторий возврата от астероида к Земле и, во‑вторых, задача анализа и поиска стабильных орбит спутников астероида. Результаты исследования второй задачи приведены в [5, 6]. В настоящей работе исследуется первая задача. Отметим, что обычно при оптимизации межпланетных траекторий рассматривается полёт только до планеты-цели без возврата к Земле. Требование возврата существенно усложняет задачу. Сначала рассматривается задача оптимизации траекторий полёта к произвольному астероиду с возвратом к Земле, затем результаты анализа будут применены для экспедиции к астероиду Апофис.

Рассмотрена следующая схема экспедиции Земля–астероид–Земля. Ракета-носитель (РН) выводит космический аппарат (КА) с разгонным блоком (РБ) большой тяги на опорную орбиту искусственного спутника Земли (ИСЗ). После пассивного движения по ней в некоторый оптимальный момент t₀ РБ сообщает КА импульс скорости ∆V₁, производится разгон КА и КА переводится на орбиту полёта к астероиду. Затем разгонный блок отделяется от КА, и дальнейшие манёвры осуществляются с помощью второй двигательной установки большой тяги (ДУ2). В момент t_1сд аппарат выходит из сферы действия Земли. Далее, в момент t₂ КА подлетает к астероиду. С помощью ДУ2 сообщается импульс скорости ∆V₂, осуществляется торможение КА, и КА переходит на орбиту искусственного спутника астероида. В окрестности астероида КА пребывает некоторое время ∆t₂₃, это “время ожидания” [7]. В течение этого времени возможны посадка на поверхность астероида, взятие образцов его грунта и другие исследования. Затем в момент t₃ основному КА сообщается импульс скорости ∆V₃, КА разгоняется и переходит на траекторию возвращения к Земле. В момент t_4сд КА подлетает к сфере действия Земли. От КА отделяется спускаемый аппарат, и в момент t_f происходит его гиперболический вход в атмосферу Земли, затем торможение, посадка. Для уменьшения энергетических затрат в качестве основной принята схема полёта, когда тормозной импульс скорости двигателем не сообщается при подлете к Земле, как предлагал К.Э. Циолковский [1] и как сделано в проектах “Stardust”, “Hayabusa”. В этом случае энергетические затраты на экспедицию в номинале определяются тремя величинами импульсов скорости ∆V₁, ∆V₂ , ∆V₃.

Обычно при оптимизации траекторий рассматривают минимизацию характеристической скорости Vхар, равной сумме величин импульсов скорости,

$V_{\text{õàð}}\left(t_{\text{1}},t_{\text{2}},t_{\text{3}},t_{\text{4}}\right)=\sum _{i=1}^n\Delta V_i\to \text{\hspace{0.17em}\hspace{0.17em}min,}$  (1)

или максимизацию конечной массы:

$-m_f\left(t_{\text{1}},t_{\text{2}},t_{\text{3}},t_{\text{4}}\right)\to \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}min.}$ (2)

Нами при построении оптимальных межпланетных траекторий перелёта в основном варианте анализа максимизируется полезная масса экспедиции mp:

$-m_p\left(t_{\text{1}},t_{\text{2}},t_{\text{3}},t_{\text{4}}\right)\to \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}min.}$ (3)

При этом полагаем, что полезная масса получается вычитанием из конечной массы mf массы ДУ2, зависящей от массы топлива, т. е. импульсов DV1, DV2, DV3. Этот функционал лучше отражает требование энергетической эффективности траектории.

Задача построения оптимальных межпланетных траекторий решается в два этапа. На первом этапе гелиоцентрические траектории перелёта КА Земля–астероид и астероид–Земля определяются в модели точечных сфер действия Земли и астероида, поэтому орбиты этих перелётов строятся в центральном поле притяжения Солнца. Схема решения задачи будет следующей. При задании граничных времён экспедиции t1 (отлёт с орбиты Земли), t2 (подлёт к орбите астероида), t3 (отлёт с орбиты астероида), t4 (подлёт к орбите Земли) гелиоцентрические орбиты перелёта между небесными телами определяются путём двукратного решения задачи Эйлера—Ламберта (с учётом возможности совершения пассивного витка хотя бы по одной орбите). Это позволяет найти скорости “на бесконечности” V¥1, V¥2, V¥3, V¥4 в граничные времена ti и требуемые импульсы скорости для перелёта DV1, DV2, DV3. По этим скоростям определяются конечная mf и полезная mp массы КА. Полезная масса КА определяется с учётом отделяемых масс РБ и ДУ2, при этом скорости истечения газов с1, с2 из двигательных установок РБ и ДУ2, вообще говоря, различны. Необходимо выбором времён t1, t2, t3, t4 (при заданных областях для этих времён) найти оптимальные траектории с максимальной полезной массой.

При этом рассмотрено несколько задач оптимизации.

1. Основная задача оптимизации сформулирована следующим образом: при заданной общей продолжительности экспедиции ∆t∑ = t4 − t1 и заданном времени пребывания КА у астероида ∆t23 = t3 − t2 оптимизируются время старта t1 и время перелёта от Земли до астероида ∆t12 = t2 − t1, чтобы выполнялось соотношение (3). Это близко к постановке, данной в работе [8].
2. При заданном времени ∆t23 и при ограничении на общую продолжительность экспедиции ∆t∑ (например, ∆t∑£ 2 года) оптимизируются времена t1, ∆t12 и ∆t∑.
3. При заданном времени ∆t∑ оптимизируются времена t1, ∆t12 и ∆t23.
4. Оптимизируются все времена ∆t∑, t1, ∆t12 и ∆t23.
5. Полная четырёхпараметрическая оптимизация времён ∆t∑, t1, ∆t12 и ∆t23 с учётом ограничения на скорость входа КА в атмосферу Земли при возврате от астероида Vвх: Vвх£ Vmax.

Чтобы обеспечить нахождение глобального оптимума, на первом этапе анализа для поиска оптимальных траекторий использовано несколько методов: метод И.М. Соболя [9], генетический алгоритм [10] и квазиньютоновский BFGS (Broyden–Fletcher–Goldfarb–Shanno) метод [11]. Метод И.М. Соболя с точками ЛП_t-последовательностей, которые очень равномерно распространены в пространстве, позволяет найти области, где лежат локальные оптимумы и глобальный оптимум. Запуск генетического алгоритма в этих областях имеет целью найти глобальный оптимум с точностью до суток. Метод BFGS позволяет затем, если надо, быстро (за 2–8 итераций) уточнить оптимум.

После построения на первом этапе оптимальных гелиоцентрических траекторий перелёта КА эти траектории проверяем на выполнение необходимых условий оптимальности в классе многоимпульсных перелётов с помощью сопряжённых функций. Получены выражения для базис-вектора Лоудена p (вектора , сопряжённого к скорос ти КА V) в граничные времена t1, t2, t3, t4 не только для обычных функционалов Vхар и mf, но и для полезной массы КА mp [12]. При минимизации характеристической скорости Vхар (1):

${\mathbf{p}}_1={\lambda }_v\left(t_1\right)=\frac{{{\displaystyle \mathit{V}}}_{\infty 1}}{V_{p1}}, {\mathbf{p}}_2={\lambda }_v\left(t_2\right)=-\frac{{{\displaystyle \mathit{V}}}_{\infty 2}}{V_{\infty 2}}, {\mathbf{p}}_3={\lambda }_v\left(t_3\right)=\frac{{{\displaystyle \mathit{V}}}_{\infty 3}}{V_{\infty 3}}, {\mathbf{p}}_4={\lambda }_v\left(t_4\right)=0.$ (4)

Здесь — скорость в перигее орбиты отлёта от Земли; притяжение астероида для простоты не учитывается. Для случая (2)

${\mathit{p}}_1=\frac{c_{2}m\left(t_0\right){\mu }_1{\mu }_2{\mu }_3}{c_1{m}_f}\frac{{{\displaystyle \mathit{V}}}_{\infty 1}}{V_{p1}},$ (5а) ${\mathit{p}}_2=-\frac{{{\displaystyle \mathit{V}}}_{\infty 2}}{V_{\infty 2}},{\mathbf{\text{p}}}_3=\frac{{{\displaystyle \mathit{V}}}_{\infty 3}}{V_{\infty 3}},{\mathbf{\text{p}}}_4=\mathrm{0,}$  (5б)

где m(t0) — начальная масса КА на опорной орбите ИСЗ;   ,  Δm1 — отделяемая после разгона у Земли масса РБ. При максимизации полезной массы КА mp (3)

$\begin{array}{l}{\mathit{p}}_1=\frac{c_{2}m\left(t_0\right){\mu }_1}{c_1{m}_f}\left({\mu }_2{\mu }_3-\frac{a_{T2}}{1+a_{T2}}\right)\frac{{\mathit{V}}_{\infty 1}}{V_{p1}},\\ {\mathit{p}}_2=-\frac{{\mathit{V}}_{\infty 2}}{V_{\infty 2}},{\mathbf{\text{p}}}_3=\frac{{\mathit{V}}_{\infty 3}}{V_{\infty 3}}, {\mathit{p}}_4=\mathrm{0,}\end{array}$ (6)

где aT2 — коэффициент пропорциональности массы топливных баков ДУ2 массе топлива. Зная эти граничные сопряжённые переменные, можно определить текущие сопряжённые переменные по переходной матрице F, удовлетворяющей уравнениям

$\left[\begin{array}{l}\delta {\mathit{r}}_f\\ \delta {\mathit{V}}_f\end{array}\right]=\left[\begin{array}{l}{\Phi }_1{\Phi }_2\\ {\Phi }_3{\Phi }_4\end{array}\right]\left[\begin{array}{l}\delta {\mathit{r}}_0\\ \delta {\mathit{V}}_0\end{array}\right]$,  (7)

где r, V — гелиоцентрические радиус-вектор и вектор скорости КА, индексы 0 и f соответствуют началу и концу каждой дуги перелёта,

${\Phi }_1=\frac{\partial {\mathit{r}}_f}{\partial {\mathit{r}}_0},{\Phi }_2=\frac{\partial {\mathit{r}}_f}{\partial {\mathit{V}}_0},{\Phi }_3=\frac{\partial {\mathit{V}}_f}{\partial {\mathit{r}}_0},{\Phi }_4=\frac{\partial {\mathit{V}}_f}{\partial {\mathit{V}}_0}$.   (7а)

Поскольку переходная матрица для (p, p¢) идентична матрице для вариаций (δr, δV), производная базис-вектора $\mathit{p}\mathbf{\text{'}}=\frac{d\mathit{p}}{dt}$ в начальный момент

$\mathit{p}{\mathbf{\text{'}}}_0={\Phi }_2^{-1}\left({\mathit{p}}_f-{\Phi }_1{\mathit{p}}_0\right)$. (8)

Тогда можно определить базис-вектор на всей траектории

$\mathbf{\text{p}}\text{(}t\text{)}={\Phi }_1(t,t_0){\mathit{p}}_0+{\Phi }_2(t,t_0)\mathit{p}{\mathbf{\text{'}}}_0$, (9а)

при этом ${\Phi }_1, {\Phi }_2$,  определяются по (7а) в каждый момент t. Можно также получить базис-вектор p(t) интегрированием системы уравнений [8]:

$\begin{array}{l}\frac{d\mathit{r}}{dt}=\mathit{V},\frac{d\mathit{V}}{dt}=-{\mu }_{\text{s}}\frac{r}{r^3},\\ \frac{d\mathit{p}}{dt}=\mathit{p}\text{'},\frac{d\mathit{p}\text{'}}{dt}=-{\mu }_{\text{s}}\frac{p}{r^3}+\left(\mathit{p}\mathbf{,}\mathit{r}\right)\frac{3{\mu }_{\text{s}}\mathit{r}}{r^5},\end{array}$ (9б)

где µs — гравитационный параметр Солнца. Для оптимальности траекторий в классе многоимпульсных перелётов необходимо выполнение условия

$p\left(t\right)=\left|\mathbf{\text{p}}\left(t\right)\right|\le \text{ 1}$ . (10)

Если это условие, играющее здесь роль “принципа максимума”, нарушается на некотором участке, то траекторию можно улучшить введением дополнительных импульсов или вариацией граничных времён [12].

Также получены производные от функционалов (1)–(3) по граничным временам траектории [12], они могут быть использованы для проверки выполнения условий трансверсальности и для улучшения траектории по функционалу, например, градиентным или квазиньютоновским методами. Для функционала J3 = −mp (3) имеем

$\begin{array}{c}\frac{\partial J_3}{\partial t_1}=-d_3{p^{\text{'}}}_1{V}_{\infty 1},\frac{\partial J_3}{\partial t_2}=d_3{p^{\text{'}}}_2{V}_{\infty 2},\\ \frac{\partial J_3}{\partial t_3}=-d_3{p^{\text{'}}}_3{V}_{\infty 3},\frac{\partial J_3}{\partial t_4}=d_3\left({\mathit{p}}_4, {\mathit{V}}_{\infty 4}\right),\end{array}$ (11)

 $\begin{array}{l}{d}_3=m_f\left(1+a_{T2}\right)/c_2,{\mathbf{\text{p}}}_i=\mathit{p}\left(t_i\right),\\ {p^{\text{'}}}_i=\frac{d\left|{\mathit{p}}_i\right|}{dt}=\frac{d}{dt}\left[{({\mathit{p}}_i\cdot {\mathit{p}}_i)}^{\frac{1}{2}}\right]=\frac{{\mathbf{\text{p'}}}_i\cdot {\mathit{p}}_i}{\left|{\mathbf{\text{p}}}_i\right|},i=\mathrm{1,}\mathrm{2,}3.\end{array}$

Если времена ti (i = 1, 2, 3, 4) не зависят друг от друга и лежат внутри допустимых областей, то для оптимальной траектории должно быть $\frac{\partial J_3}{\partial t_i}=0.$ Если между временами ti есть связи, то переходим к независимым переменным xj.

На втором этапе анализа уточняются характеристики полученных на первом этапе оптимальных межпланетных траекторий. При этом движение КА рассчитывается с учётом возмущений от притяжения небесных тел, сжатия Земли J2, давления солнечного света, координаты небесных тел определяются по эфемеридам JPL, а траектории КА определяются численным интегрированием системы дифференциальных уравнений движения КА и решением краевых задач для выполнения граничных условий. Выполнена также коррекция массово-энергетических характеристик. Для этого учтены гравитационные потери при разгоне КА у Земли из-за конечности тяги; предусмотрены дополнительные импульсы скорости на коррекцию траекторий и т. д. На множестве данных возмущённых траекторий выполняется их оптимизация по граничным временам экспедиции, методом покоординатного спуска [13]. При этом, в частности, получаются времена t0 отлёта от Земли и tf входа в атмосферу Земли при возврате.

 

Рис. 1. Изолинии полезной массы в плоскости времен t1 и ∆t12 для варианта Dt∑ = 690 сут, ∆t23 = 7 сут. в окрестности оптимума.

 

Рис. 2. Зависимости скорости Vхар и масс mf , mp от времени ∆t23.

 

Рис. 3. Изменение модуля базис-вектора p(t) на траектории № 19: (а) — перелёт от Земли до Апофиса; (б) — перелёт от Апофиса до Земли.

 

В соответствии с разработанным методом для экспедиции “Земля–Апофис–Земля” определены и исследованы энергетически оптимальные по максимуму полезной массы КА межпланетные траектории с запуском КА в течение 2019–2022 гг. Исследованы варианты с использованием РН “Союз‑2.1a”, “Союз‑2.1б”, “Зенит” и разгонного блока “Фрегат”. Для двигательной установки ДУ2 удельная тяга 304 с; постоянная составляющая массы m20 = 100 кг;в (6) коэффициент массы топливных баков аT2 = 0,15.

Сначала приведём данные для РН “Союз‑2.1a”. Для задачи оптимизации 1 фиксировано время ожидания Dt23 = t3 − t2 = 7 сут. При этом суммарное время Dt∑ выбиралось из множества TS = [390; 420; 450; 510; 540; 570; 600; 630; 660; 690; 730] сут. Решение этой задачи на первом этапе анализа дало оптимальные траектории для разных времен Dt∑. Среди них максимальная полезная масса КА mp = 272 кг получена для траектории № 16, для которой: Dt∑ = 690 сут, t1 = 24.05.2019, Dt12 = t2 − t1 = 335 сут, t4 = 13.04.2021, Vхар = = 6,618 км/с, mf = 527 кг (см. рис. 1). На второй дуге полёта делается пассивный виток, прилёт к Земле — у восходящего узла орбиты Апофиса.

Решена задача 2 трёхмерной оптимизации при заданном времени ожидания Dt23 Î [7; 30; 60; 90; 120; 130] сут, при Dt∑ ≤ 2 года. Получено, что оптимальное время ожидания КА у Апофиса ∆t23opt достигается на интервале 90–120 сут (см. рис. 2). Так, для траектории с Dt23 = 120 сут: t1 = 06.05.2020, Dt12 = = 297 сут, DtΣ = 716 сут., Vхар = 6,35 км/с, mp = 328 кг.

В задаче 3 трёхмерной оптимизации при задании Dt∑ = 690 сут для оптимального решения (траектория № 19): t1 = 23.05.2019, Dt12 = 336 сут, Dt23opt = 93 сут, Vхар = 6,519 км/с, mf = 544 кг, mp = 293 кг.

 

Таблица 1. Характеристики траекторий № 16а, 19а, 20а

Номер траектории	Dts, сут	t0	t2	t3	tf	Vхар, км/c	mf , кг	mp, кг
16а	692	21.05.2019	24.04.2020	01.05.2020	12.04.2021	6,721	492	226
19а	691	21.05.2019	24.04.2020	22.07.2020	11.04.2021	6,624	509	245
20а	715	06.05.2020	02.03.2021	23.06.2021	21.04.2022	6,447	517	290

 

Таблица 2. Конечная и полезная масса КА для траекторий № 16а, 19а, 20а с использованием разных РН

Траектория	Вариант Фрегата	mf, кг	mp, кг
№ 16a Союз-2.1а Союз-2.1б Зенит	1 2 СБ	492  604 1193	226 301 700
№ 19a Союз-2.1а Союз-2.1б Зенит	1 2 СБ	509  624 1233	245 325 748
№ 20a Союз-2.1а Союз-2.1б Зенит	1 3 СБ	517  611 1287	290 362 877

 

Таблица 3. Изменение полезной массы со скоростью V_вх

Vвх max, км/с	mp, кг
	Союз-2.1а	Союз-2.1б	Зенит
12,8   12,9 13  13,1   13,2    13,26	275,3 280 284   287,5 289 290	343,6 349 354 358 361 362	839,6 852 861   869,8  876 877

 

Решена также задача 4 полной четырёхмерной оптимизации при условиях Dt∑ ≤ 2 года, Dt23 ≥ 7 сут. Для неё получена траектория № 20 с характеристиками: ∆t∑ = 716 сут, t1 = 05.05.2020, Dt12 = 300 сут, Dt23opt = 112 сут, Vхар = 6,343 км/с, mf = 545 кг, mp = 329 кг.

Для полученных траекторий построены сопряженные функции, в частности, модуль базис-вектора. Для траектории № 16 ∆t23 = 7 сут < ∆t23opt, и на некотором начальном отрезке второй дуги, после приложения импульса ΔV3: p3¢ =  p(t) > 1, условие (10) не выполняется, траекторию можно улучшить. Из выражений (11) для частных производных от функционала видно, что если dt3 > 0, то уменьшим функционал, dJ3 < 0. Это подтверждает улучшение характеристик траектории № 19, у которой большее время t3 и время ожидания, Δt23 = 93 сут. На рис. 3 представлено изменение p(t) для этой траектории. Здесь p £ 1, условие оптимальности выполняется. Для траектории № 20 условие p(t) £ 1 также выполняется.

Далее, на втором этапе анализа уточнены полученные на первом этапе характеристики оптимальных траекторий. В окрестности граничных времён полученных траекторий выполнена оптимизация задачи на множестве уточненных траекторий методом покоординатного спуска. Получены оптимальные траектории № 16а, 19а, 20а. Характеристики этих траекторий приведены в табл. 1. Здесь Dts = tf − t0. Уточнение на втором этапе привело к некоторому уменьшению полезной массы.

В табл. 2 даны значения масс mf и mp для траекторий № 16a, 19а, 20а при выведении с РН “Союз‑2.1а”, “Союз‑2.1б”, “Зенит”. При разгоне у Земли для РН “Союз” делается два включения разгонного блока, для РН “Зенит” — три включения.

Для полученных траекторий при возврате к Земле оценена скорость входа КА в атмосферу Земли Vвх: Vвх = 12,74 км/c, 12,32 км/c, 13,26 км/c соответственно для траекторий № 16а, 19а, 20а. Отметим, что в проекте “Stardust” Vвх = 12,9 км/c, в проекте “Hayabusa” Vвх = 12,5 км/c, в проекте “OSIRIS-REx” Vвх = 12,2 км/c, т. е. современные технологии позволяют приземлять КА при скорости Vвх до 12,9 км/с. Учитывая это, решена ещё задача 5 полной четырёхпараметрической оптимизации с ограничением на скорость входа, Vвх. £ 12,8 км/с — для траектории 20а. В таблице 3 приведены результаты анализа зависимости полезной массы от Vвх (12,8 – 13,26 км/с). Отметим, что при этом для РН “Союз‑2.1б” и РН “Зенит” полезная масса остаётся большей, чем сухая масса КА “Stardust” (300 кг).

Таким образом, выполненный анализ показывает, что существует принципиальная возможность осуществить экспедицию к астероиду Апофис в течение 2019–2022 гг. на основе существующих РН типа “Союз”.

Авторы выражают икреннюю признательность сотрудникам НПО им. С.А. Лавочкина В.Г. Полю и А.В. Симонову, а также И.В. Крылову за поддержку и полезные обсуждения работы.

About the authors

V. V. Ivashkin

Institute for Applied Mathematics of the Russian Academy of Sciences; Bauman Moscow State Technical University

Author for correspondence.
Email: ivashkin@keldysh.ru
Russian Federation

Anqi Lang

Bauman Moscow State Technical University; Xi’an Jiaotong University

Email: ivashkin@keldysh.ru
Russian Federation, Moscow; China

References

Supplementary files