Космические исследования, 2019, T. 57, № 5, стр. 361-372

Методика оптимизации траекторий межпланетных перелетов с гравитационными маневрами при использовании двигателей малой тяги

А. А. Орлов *

Научно-исследовательский институт прикладной механики и электродинамики МАИ
г. Москва, Россия

* E-mail: aa.orlov@rambler.ru

Поступила в редакцию 25.03.2019
После доработки 05.04.2019
Принята к публикации 25.04.2019

Полный текст (PDF)

Аннотация

Предлагается методика поиска оптимальных траекторий с гравитационными маневрами (ГМ) для межпланетных перелетов космических аппаратов (КА) с электроракетной двигательной установкой (ЭРДУ). При этом используется непрямой метод оптимизации. Отличительной особенностью данной методики является объединение условий оптимальности в точке гравитационного маневра в рамках одной краевой задачи для двух случаев, когда высота пролетной гиперболы при гравитационном маневре меньше или равна максимальной. Такой подход позволяет существенно сократить объем необходимых вычислений при проведении оптимизации межпланетных траекторий КА, включающих гравитационные маневры. При этом рассматривается сквозная оптимизация траектории с анализом полного набора условий оптимальности в точке гравитационного маневра. Эффективность предлагаемого подхода демонстрируется на примере оптимизации межпланетных траекторий от Земли к Меркурию с гравитационным маневром у Венеры и от Земли к Юпитеру с гравитационным маневром у Земли.

DOI: 10.1134/S0023420619050054

ВВЕДЕНИЕ

Процесс поиска оптимальных траекторий в задачах межпланетных перелетов с гравитационными маневрами заметно более трудоемок по сравнению с поиском решения для задач прямого межпланетного перелета. Данное обстоятельство объясняется заметным увеличением порядка краевой задачи при добавлении в нее гравитационных маневров. В рамках использования непрямых методов оптимизации, к настоящему времени разработан ряд методик, позволяющих решать указанные задачи [1, 2]. Общим для этих методик является то, что добавление каждого гравитационного маневра приводит к добавлению девяти граничных условий к краевой задаче (при фиксированных датах старта и гравитационного маневра), что приводит к существенному росту ее размерности. Таким образом, построение эффективного алгоритма для решения подобных задач по-прежнему является актуальной проблемой. В связи с этим предлагается методика оптимизации межпланетных траекторий КА, применение которой увеличивает число граничных условий краевой задачи на пять единиц при добавлении каждого гравитационного маневра. Как следствие, использование предлагаемого подхода позволяет обеспечивать более высокую вычислительную устойчивость итерационного процесса поиска решения. Относительно не сильное увеличение порядка краевой задачи при добавлении в нее гравитационных маневров обеспечивается применением в математической модели гравитационного маневра матрицы поворота векторов, за счет чего часть условий оптимальности в точке гравитационного маневра выполняется автоматически.

1. УРАВНЕНИЯ ДВИЖЕНИЯ КА С ЭРДУ

Рассматривается гелиоцентрический участок траектории КА. Математическая модель, описывающая движение КА с ЭРДУ, использует следующие основные допущения: величины удельного импульса и тяги ЭРДУ принимается постоянными на всех участках работы ЭРДУ (рассматривается модель нерегулируемого двигателя); на направление вектора тяги ЭРДУ ограничений не накладывается; гравитационные поля планет и Солнца описываются моделью Ньютона; используются допущения метода грависфер нулевой протяженности: при исследовании гелиоцентрической траектории КА пренебрегается протяженностью грависфер других планет. Время начала гелиоцентрической траектории считается равным времени старта КА с промежуточной орбиты искусственного спутника Земли (ИСЗ), положение КА на гелиоцентрической траектории в начальный и конечный моменты времени совпадает с положениями планеты отправления и цели соответственно. При анализе гелиоцентрического участка скорость КА считается равной сумме скорости планеты, с которой он стартует и гиперболического избытка скорости, которую имеет КА при выходе из грависферы этой планеты. Такой подход оправдан для предварительного проектного расчета, поскольку заметно упрощает математическую модель движения КА и дает погрешность порядка десятых долей процента в затратах топлива по сравнению с решением задачи n-тел [3].

Рассматривается задача нулевой стыковки с планетой назначения (задача выравнивания положения и скорости КА с положением и скоростью планеты назначения).

В расчетах применяется прямоугольная инерциальная система координат. Ее роль играет гелиоцентрическая эклиптическая система координат ICRF (International Celestial Reference Frame) для эпохи J2000.0, являющаяся практической реализацией в радиодиапазоне общих требований и принципов построения системы координат ICRS (International Celestial Reference System) [4]. Стандарт ICRS был принят в августе 1997 года на 23-м съезде Международного Астрономического Союза (International Astronomical Union, IAU) и вступил в силу с 1 января 1998 года. Положения и скорости планет определялись с помощью эфемерид JPL (Jet Propulsion Laboratory) DE405 [5].

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

$\begin{gathered} \frac{{d{\mathbf{r}}}}{{dt}} = {\mathbf{V}}, \\ \frac{{d{\mathbf{V}}}}{{dt}} = - \frac{1}{{{{r}^{3}}}} \cdot {\mathbf{r}} + \frac{{\delta \cdot {{P}_{0}}}}{m} \cdot {\mathbf{e}}, \\ \frac{{dm}}{{dt}} = - \frac{{\delta \cdot {{P}_{0}}}}{{{{w}_{0}}}}, \\ \end{gathered} $
где r – величина радиуса вектора КА, V – величина вектора гелиоцентрической скорости КА, e –орт вектора тяги ЭРДУ, δ – функция включения-выключения двигателя $\left( {\delta \in \left[ {0,1} \right]} \right),$ w0 – величина скорости истечения ЭРДУ, P0 – величина тяги ЭРДУ; m – масса КА; t – время.

В качестве критерия оптимизации рассматривается максимизация конечной массы КА, т.е. рассматривается функционал следующего вида: J = m(tk) → max, где tk – время прибытия КА к цели.

Оптимизационная задача сводится к краевой с помощью принципа максимума Понтрягина. Гамильтониан задачи оптимального управления имеет вид:

$H = {\mathbf{p}}_{{\mathbf{r}}}^{T} \cdot {\mathbf{V}} + {\mathbf{p}}_{{\mathbf{V}}}^{T}\left( { - \frac{1}{{{{r}^{3}}}} \cdot {\mathbf{r}} + \frac{{\delta \cdot {{P}_{0}}}}{m} \cdot {\mathbf{e}}} \right) - {{p}_{m}}\frac{{\delta \cdot {{P}_{0}}}}{{{{w}_{0}}}},$
где ${\mathbf{p}}_{{\mathbf{r}}}^{T}$ = [prx, pry, prz], ${\mathbf{p}}_{{\mathbf{V}}}^{T}$ = [pVx, pVy, pVz], pm – сопряженные переменные к радиус-вектору, вектору скорости и массе КА соответственно.

При этом дифференциальные уравнения оптимального движения КА определяются по следующим соотношениям:

$\begin{gathered} \frac{{d{\mathbf{r}}}}{{dt}} = \frac{{\partial H}}{{\partial {{{\mathbf{p}}}_{r}}}},\,\,\,\,\frac{{d{\mathbf{V}}}}{{dt}} = \frac{{\partial H}}{{\partial {{{\mathbf{p}}}_{{\mathbf{V}}}}}},\,\,\,\,\frac{{dm}}{{dt}} = \frac{{\partial H}}{{\partial {{p}_{m}}}}, \\ \frac{{d{{{\mathbf{p}}}_{{\mathbf{r}}}}}}{{dt}} = - \frac{{\partial H}}{{\partial {\mathbf{r}}}},\,\,\,\,\frac{{d{{{\mathbf{p}}}_{{\mathbf{V}}}}}}{{dt}} = - \frac{{\partial H}}{{\partial {\mathbf{V}}}},\,\,\,\,\frac{{d{{p}_{m}}}}{{dt}} = - \frac{{\partial H}}{{\partial m}}. \\ \end{gathered} $

Из условия максимума гамильтониана получаем оптимальный закон управления ортом вектора тяги.

${{{\mathbf{e}}}^{{{\text{opt}}}}} = \frac{{{{{\mathbf{p}}}_{{\mathbf{V}}}}}}{{{{p}_{V}}}},\,\,\,\,{{p}_{V}} = \left| {{{{\mathbf{p}}}_{{\mathbf{V}}}}} \right|.$

Закон включения-выключения двигателя имеет вид:

$\Psi = \frac{{{{p}_{V}}}}{m} - \frac{{{{p}_{m}}}}{w},\,\,\,\,{{\delta }^{{{\text{opt}}}}} = \left\{ \begin{gathered} 1, {\text{ е с л и }}\psi > 0; \hfill \\ 0, {\text{е с л и }} {\text{ }}\psi \leqslant 0, \hfill \\ \end{gathered} \right.$
особые режимы управления не рассматриваются.

Таким образом, система дифференциальных уравнений оптимального движения КА записывается следующим образом:

$\begin{gathered} \frac{{d{\mathbf{r}}}}{{dt}} = {\mathbf{V}},\,\,\,\,\frac{{d{\mathbf{V}}}}{{dt}} = - \frac{{\mathbf{r}}}{{{{r}^{3}}}} + \frac{{\delta \cdot {{P}_{0}}}}{m}\frac{{{{{\mathbf{p}}}_{{\mathbf{V}}}}}}{{{{p}_{V}}}}, \\ \frac{{dm}}{{dt}} = - \frac{{\delta \cdot {{P}_{0}}}}{{{{w}_{0}}}}, \\ \frac{{d{{{\mathbf{p}}}_{{\mathbf{r}}}}}}{{dt}} = \frac{{{{{\mathbf{p}}}_{{\mathbf{V}}}}}}{{{{r}^{3}}}} - \frac{{3 \cdot {\mathbf{r}} \cdot ({{p}_{{Vx}}} \cdot x + {{p}_{{Vy}}} \cdot y + {{p}_{{Vz}}} \cdot z)}}{{{{r}^{5}}}}, \\ \frac{{d{{{\mathbf{p}}}_{{\mathbf{V}}}}}}{{dt}} = - {{{\mathbf{p}}}_{{\mathbf{r}}}},\,\,\,\,\frac{{d{{p}_{m}}}}{{dt}} = \delta \cdot \frac{{{{P}_{0}} \cdot {{p}_{V}}}}{{{{m}^{2}}}}. \\ \end{gathered} $

Конечные условия для задачи прямого перелета имеют вид:

(1.1)
${\mathbf{r}}({{t}_{k}}) = {{{\mathbf{r}}}_{{\mathbf{k}}}},\,\,\,\,{\mathbf{V}}({{t}_{k}}) = {{{\mathbf{V}}}_{k}},\,\,\,\,{{p}_{m}}\left( {{{t}_{k}}} \right) = 1,$
где tk время прибытия КА к цели.

Начальные условия для анализа гелиоцентрической траектории полета с учетом условия оптимальности на направление начального гиперболического избытка скорости имеют вид:

(1.2)
${\mathbf{r}}({{t}_{0}}) = {{{\mathbf{r}}}_{0}},\,\,\,\,{\mathbf{V}}({{t}_{0}}) = {{{\mathbf{V}}}_{0}} + {{V}_{\infty }}({{t}_{0}}) \cdot \frac{{{{{\mathbf{p}}}_{{\mathbf{V}}}}({{t}_{0}})}}{{{{p}_{V}}({{t}_{0}})}},$
где V(t0) – гиперболический избыток скорости КА при отлете от планеты отправления в момент t0, ${{p}_{V}}({{t}_{0}}) = \left| {{{p}_{V}}({{t}_{0}})} \right|.$

Значения сопряженных переменных в начальной точке являются неизвестными параметрами рассматриваемой краевой задачи: ${{{\mathbf{p}}}_{{\mathbf{r}}}}({{t}_{0}}),$ ${{{\mathbf{p}}}_{{\mathbf{V}}}}({{t}_{0}}),$ ${{p}_{m}}({{t}_{0}}).$

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

– условие оптимальности для даты старта t0 при фиксированном общем времени перелета имеет вид:

(1.3)
$\begin{gathered} \left( {\frac{{{{w}_{0}} \cdot {{p}_{V}}({{t}_{k}})}}{{m({{t}_{k}})}} - {{p}_{m}}({{t}_{k}})} \right) \cdot \frac{{\delta ({{t}_{k}}) \cdot {{P}_{0}}}}{{{{w}_{0}}}} - \\ - \,\,\left( {\frac{{{{w}_{0}} \cdot {{p}_{V}}({{t}_{0}})}}{{m({{t}_{0}})}} - {{p}_{m}}({{t}_{0}})} \right) \cdot \frac{{\delta ({{t}_{0}}) \cdot {{P}_{0}}}}{{{{w}_{0}}}} - \\ - \,\,\frac{{{{V}_{\infty }}({{t}_{0}})}}{{{{p}_{V}}({{t}_{0}})}} \cdot {{{\mathbf{p}}}_{{\mathbf{r}}}}{{({{t}_{0}})}^{T}} \cdot {{{\mathbf{p}}}_{{\mathbf{V}}}}({{t}_{0}}) = 0. \\ \end{gathered} $

– условие оптимальности величины начального гиперболического избытка скорости:

(1.4)
${{{\mathbf{p}}}_{{\mathbf{V}}}}({{t}_{0}}) + {{p}_{m}}({{t}_{0}}) \cdot \frac{{\partial {{m}_{0}}}}{{\partial {{V}_{{\infty 0}}}}} = 0.$

Частная производная начальной массы КА m0 (масса после отделения химического разгонного блока (ХРБ), выводящего КА на гелиоцентрическую траекторию) по величине гиперболического избытка скорости при старте от планеты отправления V0 находится следующим образом:

$\begin{gathered} \frac{{\partial {{m}_{0}}}}{{\partial {{V}_{{\infty 0}}}}} = - \frac{{{{m}_{{00}}} \cdot {{V}_{\infty }}({{t}_{0}})}}{{{{w}_{{{\text{Х Р Б 0}}}}} \cdot \sqrt {\frac{{2 \cdot {{\mu }_{{{\text{п л 0}}}}}}}{{{{r}_{0}}}} + V_{\infty }^{2}({{t}_{0}})} }} \times \\ \times \,\,\exp \left( {\frac{{\sqrt {\frac{{2 \cdot {{\mu }_{{{\text{п л }}0}}}}}{{{{r}_{0}}}} + V_{\infty }^{2}({{t}_{0}})} - \sqrt {\frac{{{{\mu }_{{{\text{п л 0}}}}}}}{{{{r}_{0}}}}} }}{{{{w}_{{{\text{Х Р Б 0}}}}}}}} \right), \\ \end{gathered} $
где m00 – масса, которую выводит РН на базовую орбиту планеты отправления с радиусом r0, wХРБ0 – скорость истечения ДУ ХРБ, µпл0 – гравитационный параметр планеты отправления. При проведении расчетов все величины используются в безразмерном виде. Следует отметить, что последнее соотношение справедливо в случае пренебрежения потерями в скорости при работе химического разгонного блока, обеспечивающего старт с круговой орбиты.

Условия (1.3) и (1.4) добавляются к граничным условиям краевой задачи, а к варьируемым переменным, соответственно, добавляются дата старта t0 и величина вектора начального гиперболического избытка скорости V(t0): ${{t}_{0}},$ ${{V}_{\infty }}({{t}_{0}}).$

2. НЕОБХОДИМЫЕ УСЛОВИЯ ОПТИМАЛЬНОСТИ В ТОЧКЕ ГРАВИТАЦИОННОГО МАНЕВРА

В связи с тем, что анализируется траектория КА с ЭРДУ, а время длительности планетоцентрического участка при пролете планеты относительно мало и ЭРДУ КА не успеет дать существенное приращение скорости, будем рассматривать пассивный гравитационный маневр (ГМ).

Исходя из метода грависфер нулевой протяженности при пассивном ГМ вектор гелиоцентрической скорости КА изменяется мгновенно. Это изменение заключается в повороте вектора гиперболического избытка скорости под действием сил гравитации планеты на угол β, максимальная величина которого βmax задается минимальной высотой гиперболы пролета. Разворот вектора гиперболического избытка скорости происходит в некоторой плоскости, которая определяется векторами подлетного и отлетного гиперболических избытков скорости и называется плоскостью ГМ. Схема поворота вектора гиперболического избытка скорости КА при гравитационном маневре показана на рис. 1.

Рис. 1.

Схема гравитационного маневра.

Плоскость ГМ может разворачиваться на произвольный угол вокруг оси х за счет выбора положения точки пересечения траектории КА с картинной плоскостью пролета планеты.

Угол между векторами подлетного ${\mathbf{V}}_{{\infty - }}^{{}}$ и отлетного ${\mathbf{V}}_{{\infty + }}^{{}}$ гиперболических избытков скорости при ГМ определяется следующим соотношением:

$\beta = 2 \cdot \arcsin \left( {{{{\left( {1 + \frac{{{{r}_{p}} \cdot {{{\left| {{{V}_{\infty }}} \right|}}^{2}}}}{{{{\mu }_{{pl}}}}}} \right)}}^{{ - 1}}}} \right),$
где $\left| {{{V}_{\infty }}} \right|$ – модуль вектора гиперболического избытка скорости (ГИС), rp – радиус перицентра пролетной гиперболы, μpl – гравитационный параметр планеты. Поскольку радиус пролета планеты ограничен некоторым минимально допустимым значением rpmin, то угол поворота ГИС не может превышать величину угла βmax:

$\beta \leqslant {{\beta }_{{\max }}} = 2 \cdot \arcsin \left( {{{{\left( {1 + \frac{{{{r}_{p}} \cdot {{{\left| {{{{\mathbf{V}}}_{\infty }}} \right|}}^{2}}}}{{{{\mu }_{{pl}}}}}} \right)}}^{{ - 1}}}} \right).$

Абсолютная величина гиперболического избытка скорости при гравитационном маневре остается неизменной:

(2.1)
$\left| {{{{\mathbf{V}}}_{{\infty - }}}} \right| = \left| {{{{\mathbf{V}}}_{{\infty + }}}} \right|.$

В равенстве 2.1 и далее индекс “–” обозначает значение какой-либо величины до ГМ, индекс “+” после.

При этом в точке ГМ нужно удовлетворить ряду необходимых условий оптимальности, вид которых меняется в зависимости от того равен угол поворота гиперболического избытка скорости максимальному β = βmax или меньше β < βmax. Необходимые условия оптимальности будут рассматриваться в форме, предложенной в работах [6, 7].

В случае, если угол поворота гиперболического избытка скорости меньше максимального β < βmax, необходимые условия оптимальности принимают вид:

– базис-вектор (вектор, сопряженный вектору гелиоцентрической скорости КА) при подлете к планете до ГМ коллинеарен вектору подлетного гиперболического избытка скорости:

(2.2)
${{{\mathbf{p}}}_{{{\mathbf{V}} - }}}||{{{\mathbf{V}}}_{{\infty - }}},$

– базис-вектор при отлете от планеты после ГМ коллинеарен вектору отлетного гиперболического избытка скорости:

(2.3)
${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}||{{{\mathbf{V}}}_{{\infty + }}},$

– подлетные и отлетные модули базис-векторов равны:

(2.4)
$\left| {{{{\mathbf{p}}}_{{{\mathbf{V}} - }}}} \right| = \left| {{{{\mathbf{p}}}_{{{\mathbf{V}} + }}}} \right|.$

В случае, если угол поворота гиперболического избытка скорости равен максимальному β = βmax, необходимые условия оптимальности принимают вид:

– базис-вектор при подлете к планете для ГМ ${{{\mathbf{p}}}_{{{\mathbf{V}} - }}}$ принадлежит плоскости ГМ, образованной векторами подлетного ${{{\mathbf{V}}}_{{\infty - }}}$ и отлетного ${{{\mathbf{V}}}_{{\infty + }}}$ гиперболического избытка скорости:

(2.5)
${{{\mathbf{p}}}_{{{\mathbf{V}} - }}} \in {{\Pi }_{{{{{\mathbf{V}}}_{{\infty - }}}{{{\mathbf{V}}}_{{\infty + }}}}}},$

– базис-вектор при отлете от планеты после ГМ ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}$ также принадлежит плоскости ГМ:

(2.6)
${{{\mathbf{p}}}_{{{\mathbf{V}} + }}} \in {{\Pi }_{{{{{\mathbf{V}}}_{{\infty - }}}{{{\mathbf{V}}}_{{\infty + }}}}}}.$

Для ввода остальных условий оптимальности ГМ введем в рассмотрение по две компоненты каждого из двух рассматриваемых базис-векторов.

Одна из этих компонент есть проекция базис-вектора на направление вектора гиперболического избытка скорости. Назовем ее коллинеарной компонентой $p_{V}^{{||}},$ она может быть как положительной, так и отрицательной. Вторая компонента $p_{V}^{ \bot }$ перпендикулярна вектору гиперболического избытка скорости и лежит в плоскости ГМ и, соответственно, называется перпендикулярной компонентой (см. рис. 1).

Условие оптимальности ГМ имеют вид:

– величина перпендикулярной компоненты базис-вектора при подлете к планете равна величине перпендикулярной компоненты при отлете от планеты, и они всегда направлены в сторону гравитационного центра:

(2.7)
$p_{{V - }}^{ \bot } = p_{{V + }}^{ \bot },$

– связь величин подлетных и отлетных параллельных компонент базис-векторов имеет вид:

(2.8)
$p_{{V + }}^{{||}} = p_{{V - }}^{{||}} + A \cdot p_{{V - }}^{ \bot },$

где $A = 4 \cdot {\text{tg}}\left( {\frac{\beta }{2}} \right) \cdot \left( {1 - \sin \left( {\frac{\beta }{2}} \right)} \right).$

3. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ГРАВИТАЦИОННОГО МАНЕВРА

Для примера рассмотрим схему межпланетного перелета, включающую один ГМ. Условно разобьем траекторию точкой ГМ на два участка: первый участок от планеты старта до планеты, в окрестности которой совершается ГМ, второй участок от планеты ГМ до планеты назначения.

После расчета первого участка траектории нам становятся известны компоненты вектора гиперболического избытка скорости подлета к промежуточной планете ${{{\mathbf{V}}}_{{\infty - }}}.$ Зная величину и направление ${\mathbf{V}}_{{\infty - }}^{{}}$ можно однозначно определить величину и направление ${\mathbf{V}}_{{\infty + }}^{{}}$ путем ввода двух параметров ГМ: угла поворота гиперболического избытка скорости КА при гравитационном маневре β и угла поворота плоскости пролетной гиперболы γ в произвольной, например, эклиптической системе координат. При этом $\beta \in [0,{{\beta }_{{\max }}}],$ $\gamma \in [0,{\text{ 2}}\pi ].$

На рис. 2 [8] показана плоскость ГМ, которой принадлежат векторы ${\mathbf{V}}_{{\infty - }}^{{}}$ и ${{{\mathbf{V}}}_{{\infty + }}}.$ Линия пересечения картинной плоскости с плоскостью эклиптики XY обозначена l, линия пересечения картинной плоскости с плоскостью ГМ обозначена m. Анализируемый угол γ есть угол между l и m, отсчитываемый от l к m против часовой стрелки, если смотреть с конца вектора ${{{\mathbf{V}}}_{{\infty - }}}$ (он перпендикулярен картинной плоскости, в которой расположены l и m).

Рис. 2.

Графическое изображение углов β и γ.

Компоненты вектора ${\mathbf{V}}_{{\infty + }}^{{}}$ определяются из следующего соотношения [9]:

(3.1)
$\begin{gathered} {{{\mathbf{V}}}_{{\infty + }}} = \left( {\begin{array}{*{20}{c}} {{{V}_{{\infty - x}}}}&{\frac{{ - {{V}_{{\infty - y}}} \cdot \left| {{{V}_{{\infty - }}}} \right|}}{{\sqrt {V_{{\infty - x}}^{2} + V_{{\infty - y}}^{2}} }}}&{\frac{{ - {{V}_{{\infty - x}}} \cdot {{V}_{{\infty - z}}}}}{{\sqrt {V_{{\infty - x}}^{2} + V_{{\infty - y}}^{2}} }}} \\ {{{V}_{{\infty - y}}}}&{\frac{{{{V}_{{\infty - x}}} \cdot \left| {{{V}_{{\infty - }}}} \right|}}{{\sqrt {V_{{\infty - x}}^{2} + V_{{\infty - y}}^{2}} }}}&{\frac{{ - {{V}_{{\infty - y}}} \cdot {{V}_{{\infty - z}}}}}{{\sqrt {V_{{\infty - x}}^{2} + V_{{\infty - y}}^{2}} }}} \\ {{{V}_{{\infty - z}}}}&0&{\sqrt {V_{{\infty - x}}^{2} + V_{{\infty - y}}^{2}} } \end{array}} \right) \times \\ \times \,\,\left( {\begin{array}{*{20}{c}} {\cos (\beta )} \\ {\sin (\beta ) \cdot \cos (\gamma )} \\ {\sin (\beta ) \cdot \sin (\gamma )} \end{array}} \right), \\ \end{gathered} $
где V∞–x, V∞–y, V∞–z – компоненты вектора гиперболического избытка скорости при подлете к планете для ГМ.

Аналогичным образом, зная компоненты подлетного базис-вектора ${{{\mathbf{p}}}_{{{\mathbf{V}} - }}}$ и величины углов β и γ, можно однозначно определить компоненты вектора ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}.$

В случае, если β < βmax компоненты вектора ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}$ определяются из следующего соотношения:

(3.2)
$\begin{gathered} {{{\mathbf{p}}}_{{{\mathbf{V}} + }}} = \left( {\begin{array}{*{20}{c}} {{{p}_{{V - x}}}}&{\frac{{ - {{p}_{{V - y}}} \cdot \left| {{{{\mathbf{p}}}_{{{\mathbf{V}} - }}}} \right|}}{{\sqrt {p_{{V - x}}^{2}\, + \,p_{{V - y}}^{2}} }}}&{\frac{{ - {{p}_{{V - x}}} \cdot {{p}_{{V - z}}}}}{{\sqrt {p_{{V - x}}^{2}\, + \,p_{{V - y}}^{2}} }}} \\ {{{p}_{{V - y}}}}&{\frac{{{{p}_{{V - x}}} \cdot \left| {{{{\mathbf{p}}}_{{{\mathbf{V}} - }}}} \right|}}{{\sqrt {p_{{V - x}}^{2}\, + \,p_{{V - y}}^{2}} }}}&{\frac{{ - {{p}_{{V - y}}} \cdot {{p}_{{V - z}}}}}{{\sqrt {p_{{V - x}}^{2}\, + \,p_{{V - y}}^{2}} }}} \\ {{{p}_{{V - z}}}}&0&{\sqrt {p_{{V - x}}^{2}\, + \,p_{{V - y}}^{2}} } \end{array}} \right) \times \\ \times \,\,\left( {\begin{array}{*{20}{c}} {\cos (\beta )} \\ {\sin (\beta )\cos (\gamma )} \\ {\sin (\beta )\sin (\gamma )} \end{array}} \right), \\ \end{gathered} $
где pV–x, pV–y, pVz – компоненты базис-вектора при подлете к планете для ГМ.

В рассматриваемом случае нужно выполнить необходимые условия оптимальности (2.2), (2.3). Условия (2.1) и (2.4) согласно (3.1) и (3.2) выполняются автоматически.

Математическая модель ГМ однозначно связывает условия коллинеарности (2.2) и (2.3) через матрицы поворота векторов (3.1) и (3.2). Из этого следует, что достаточно выполнить одно из двух условий (2.2) или (2.3) и второе условие выполнится автоматически.

Таким образом, если использовать для определения ${{{\mathbf{V}}}_{{\infty + }}}$ и ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}$ соотношения (3.1) и (3.2), то в точке ГМ для случая β < βmax граничные условия можно рассматривать в виде:

(3.3)
$\begin{gathered} {\mathbf{r}}\left( {{{t}_{i}}} \right) - {{{\mathbf{r}}}_{{{\text{п л }}i}}}\left( {{{t}_{i}}} \right) = 0, \\ p_{{V + y}}^{i} \cdot {\mathbf{V}}_{{\infty + x}}^{i} - p_{{V + x}}^{i} \cdot {\mathbf{V}}_{{\infty + y}}^{i} = 0, \\ p_{{V + z}}^{i} \cdot {\mathbf{V}}_{{\infty + x}}^{i} - p_{{V + x}}^{i} \cdot {\mathbf{V}}_{{\infty + z}}^{i} = 0, \\ \end{gathered} $
где $r\left( {{{t}_{i}}} \right)$ – радиус-вектор КА, ti – момент времени совершения ГМ, $i = 1,...,n,$ где n – количество ГМ, $p_{{V + x}}^{i},p_{{V + y}}^{i},p_{{V + z}}^{i}$ – компоненты базис-вектора ${\mathbf{p}}_{{{\mathbf{V}} + }}^{i},$ ${\mathbf{V}}_{{\infty + x}}^{i},{\mathbf{V}}_{{\infty + y}}^{i},{\mathbf{V}}_{{\infty + z}}^{i}$ – компоненты вектора ${\mathbf{V}}_{{\infty + }}^{i}.$ Всего получается 5 граничных условий.

Неизвестными параметрами в точке ГМ являются:

(3.4)
${{{\mathbf{p}}}_{{\mathbf{r}}}}({{t}_{i}}),{{\beta }_{i}},{{\gamma }_{i}}.$

Всего получается 5 неизвестных параметров.

При включении в схему перелета ГМ к граничным условиям (1.1) и (1.2) необходимо добавить условия (3.3) и (3.4). Всего задача межпланетного перелета с одним ГМ будет иметь 11 граничных условий и 11 неизвестных, каждый дополнительный ГМ увеличивает порядок краевой задачи на 5 (если не считать неизвестной дату ГМ).

Последние два слагаемые в (3.3) отвечают за выполнение необходимого условия оптимальности (2.3) – коллинеарности отлетных вектора гиперболического избытка скорости и базис-вектора. Многочисленные практические расчеты показали, что это обеспечивает гораздо лучшую сходимость поиска решения краевой задачи по сравнению с выбором условия (2.2) – коллинеарности подлетных вектора гиперболического избытка скорости и базис-вектора. В самом деле, добавляя в краевую задачу условия коллинеарности (2.2) (это два скалярных условия) необходимо добавить две варьируемые переменные, углы β и γ, но при этом условие (2.2) не зависит от этих углов и, следовательно, мы уменьшаем число степеней свободы для выполнения данного ограничения, что отрицательно сказывается на сходимости итерационного процесса.

Теперь рассмотрим случай, когда β = βmax. В точке ГМ должно быть выполнено условие (2.1) и необходимые условия оптимальности (2.5)–(2.8). Так же как и для случая β < βmax, компоненты вектора ${{{\mathbf{V}}}_{{\infty + }}}$ определяются из соотношения (3.1) по известным компонентам подлетного вектора гиперболического избытка скорости ${{{\mathbf{V}}}_{{\infty - }}}$, полученного после расчета первого участка траектории до ГМ, таким образом, мы выполняем условие (2.1). Теперь нужно найти компоненты вектора ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}$ удовлетворяющие необходимым условиям оптимальности. Вектор ${{{\mathbf{p}}}_{{{\mathbf{V}} - }}}$ считается известным после расчета участка траектории до ГМ.

Компоненты подлетного и отлетного базиса векторов можно найти из следующих соотношений:

(3.5)
$p_{{V \pm }}^{ \bot } = {\mathbf{p}}_{{{\mathbf{V}} \pm }}^{T} \cdot {\mathbf{e}}_{ \pm }^{ \bot },\,\,\,\,p_{{V \pm }}^{{||}} = {\mathbf{p}}_{{{\mathbf{V}} \pm }}^{T} \cdot {\mathbf{e}}_{ \pm }^{{||}},$
где ${\mathbf{e}}_{ \pm }^{{||}} = \frac{{{\mathbf{V}}_{\infty }^{ \pm }}}{{\left| {{\mathbf{V}}_{\infty }^{ \pm }} \right|}},$ ${{{\mathbf{e}}}^{b}} = \frac{{{\mathbf{V}}_{\infty }^{ - } \times {\mathbf{V}}_{\infty }^{ + }}}{{\left| {{\mathbf{V}}_{\infty }^{ - } \times {\mathbf{V}}_{\infty }^{ + }} \right|}},$ ${\mathbf{e}}_{ \pm }^{ \bot } = \frac{{{\mathbf{e}}_{{}}^{b} \times {\mathbf{V}}_{\infty }^{ \pm }}}{{\left| {{\mathbf{e}}_{{}}^{b} \times {\mathbf{V}}_{\infty }^{ \pm }} \right|}}.$

По соотношениям (3.5) найдем перпендикулярную $p_{{V - }}^{ \bot }$ и коллинеарную $p_{{V - }}^{{||}}$ компоненты подлетного базис-вектора. Используя (2.7) и (2.8) найдем перпендикулярную $p_{{V + }}^{ \bot }$ и коллинеарную $p_{{V + }}^{{||}}$ компоненты отлетного базис-вектора ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}$ и вычислим модули базис-векторов ${{{\mathbf{p}}}_{{{\mathbf{V}} - }}}$ и ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}$:

$\left| {{{{\mathbf{p}}}_{{{\mathbf{V}} \pm }}}} \right| = \sqrt {p_{{V \pm }}^{{||}} + p_{{V \pm }}^{ \bot }} .$
Таким образом, мы определили абсолютные величины векторов ${{{\mathbf{p}}}_{{{\mathbf{V}} - }}}$ и ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}.$ Как известно, любой вектор – это величина, которая характеризуется своим абсолютным значением и направлением. В соответствии с этим, нам осталось определить направление базис-вектора ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}$ и его компоненты станут нам известны.

Введем вспомогательные углы α1 и α2, характеризующие угловое положение подлетного и отлетного базис-векторов относительно подлетного и отлетного гиперболического избытка скорости соответственно (см. рис. 1).

Найдем углы α1 и α2 по следующим соотношениям:

$\alpha 1 = {\text{arctg}}\left( {\frac{{p_{{V - }}^{ \bot }}}{{{{p}_{{V - }}}}}} \right),$
$\alpha 2 = {\text{arctg}}\left( {\frac{{p_{{V + }}^{ \bot }}}{{{{p}_{{V + }}}}}} \right).$
Найдем угол α, характеризующий разницу в угловом положении подлетного и отлетного базис-векторов относительно соответствующих им гиперболического избытка скорости: $\alpha = \alpha 2 - \alpha 1.$

Повернем вектор ${{{\mathbf{p}}}_{{{\mathbf{V}} - }}}$ на углы β + α и γ, тем самым определив направление вектора ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}$ и умножим полученную величину на отношение модуля отлетного базис-вектора к модулю подлетного базис-вектора:

(3.6)
$\begin{gathered} {{{\mathbf{p}}}_{{{\mathbf{V}} + }}} = \left( {\begin{array}{*{20}{c}} {{{p}_{{V - x}}}}&{\frac{{ - {{p}_{{V - y}}} \cdot \left| {{{{\mathbf{p}}}_{{V - }}}} \right|}}{{\sqrt {p_{{V - x}}^{2} + p_{{V - y}}^{2}} }}}&{\frac{{ - {{p}_{{V - x}}} \cdot {{p}_{{V - z}}}}}{{\sqrt {p_{{V - x}}^{2} + p_{{V - y}}^{2}} }}} \\ {{{p}_{{V - y}}}}&{\frac{{{{p}_{{V - x}}} \cdot \left| {{{{\mathbf{p}}}_{{V - }}}} \right|}}{{\sqrt {p_{{V - x}}^{2} + p_{{V - y}}^{2}} }}}&{\frac{{ - {{p}_{{V - y}}} \cdot {{p}_{{V - z}}}}}{{\sqrt {p_{{V - x}}^{2} + p_{{V - y}}^{2}} }}} \\ {{{p}_{{V - z}}}}&0&{\sqrt {p_{{V - x}}^{2} + p_{{V - y}}^{2}} } \end{array}} \right) \times \\ \times \,\,\left( {\begin{array}{*{20}{c}} {\cos (\beta + \alpha )} \\ {\sin (\beta + \alpha ) \cdot \cos (\gamma )} \\ {\sin (\beta + \alpha ) \cdot \sin (\gamma )} \end{array}} \right) \cdot \frac{{\left| {{{{\mathbf{p}}}_{{{\mathbf{V}} + }}}} \right|}}{{\left| {{{{\mathbf{p}}}_{{{\mathbf{V}} - }}}} \right|}}. \\ \end{gathered} $

Таким образом, по соотношению (3.6) мы определили компоненты базис-вектора ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}.$ При этом необходимые условия оптимальности (2.7), (2.8) выполняются автоматически. Достаточно выполнить одно из двух необходимых условий оптимальности (2.5) или (2.6) и другое условие выполнится автоматически. В самом деле, векторы ${\mathbf{V}}_{{\infty - }}^{{}}$ и ${\mathbf{V}}_{{\infty + }}^{{}}$ образуют плоскость ГМ и, если мы повернем вектор ${\mathbf{V}}_{{\infty - }}^{{}}$ на углы β и γ, используя соотношение (3.1), то вектор ${\mathbf{V}}_{{\infty + }}^{{}}$ попадет в плоскость ГМ по определению. Аналогично, предположим, что базис-вектор ${{{\mathbf{p}}}_{{{\mathbf{V}} - }}}$ принадлежит плоскости ГМ (т.е. выполняется условие (2.5)) и если его повернуть на углы β + α и γ, используя соотношение (3.6), то полученный вектор ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}$ так же будет принадлежать плоскости ГМ, поскольку углы β и α находятся в одной плоскости – плоскости ГМ, а угол γ, использующийся в соотношениях (3.1) и (3.6), разумеется, один и тот же.

Практика показала, что желательно выбирать в качестве граничного условие (2.6), поскольку в нем два вектора из трех зависят от углов β и γ (в условии (2.5) только один вектор), которые являются выбираемыми параметрами, что обеспечивает большую гибкость процесса поиска решения и, следовательно, его лучшую сходимость.

Таким образом, в точке ГМ для случая, когда угол поворота гиперболического избытка скорости КА при гравитационном маневре β = βmax, мы имеем следующие граничные условия:

(3.7)
$\begin{gathered} {\mathbf{r}}\left( {{{t}_{i}}} \right) - {{{\mathbf{r}}}_{{{\text{п л }}i}}}\left( {{{t}_{i}}} \right) = 0, \\ {\mathbf{p}}_{{V + i}}^{T} \cdot \left[ {{{{\mathbf{V}}}_{{\infty - i}}} \times {{{\mathbf{V}}}_{{\infty + i}}}} \right] = 0, \\ {{\beta }_{{\max i}}} - {{\beta }_{i}} = 0, \\ \end{gathered} $
где $i = 1,...,n,$ n – количество ГМ. Всего получается 5 граничных условий.

Неизвестными параметрами в точке ГМ являются параметры (3.4). Всего получается 5 неизвестных параметров.

При включении в схему перелета ГМ к граничным условиям (1.1) и (1.2) необходимо добавить условия (3.7) и (3.4). Всего задача межпланетного перелета с одним ГМ будет иметь 11 граничных условий и 11 неизвестных, каждый дополнительный ГМ увеличивает порядок краевой задачи на 5 (если считать дату ГМ известной).

Представленная методика позволяет проводить сквозную оптимизацию всей траектории с выполнением всех необходимых условий оптимальности в точке ГМ.

В представленных ниже численных примерах так же использовались необходимые условия оптимальности для даты проведения ГМ ti:

(3.8)
${\mathbf{p}}_{{r + }}^{i}{{({{t}_{i}})}^{T}}{\mathbf{V}}_{{\infty + }}^{i} - {\mathbf{p}}_{{r - }}^{i}{{({{t}_{i}})}^{T}}{\mathbf{V}}_{{\infty - }}^{i} = 0.$

Условие (3.8) может быть добавлено к остальным граничным условиям краевой задачи, а дополнительной неизвестной переменной станет, соответственно, дата проведения ГМ.

4. ИСПОЛЬЗОВАНИЕ ОГРАНИЧЕНИЙ СМЕШАННОГО ТИПА В ТОЧКЕ ГРАВИТАЦИОННОГО МАНЕВРА

Мы рассмотрели два отдельных случая в точке ГМ, когда β < βmax и β = βmax. Объединение этих случаев позволит заметно снизить трудоемкость оптимизации траекторий, включающих ГМ. Действительно, при одном гравитационном маневре нам необходимо рассмотреть два случая, когда β < βmax и β = βmax и выбрать из них тот, при котором не нарушаются граничные условия β ≤ βmax и условия оптимальности. При оптимизации перелета с произвольным числом ГМ, нам необходимо рассмотреть серию краевых задач, количество которых определяется из следующего соотношения: $k = {{2}^{n}},$ где n – количество ГМ.

В общем случае угол поворота ГИС β при гравитационном маневре согласно (2.2) должен удовлетворять условию:

(4.1)
${{\beta }_{{\max }}} - \beta \geqslant 0.$

Перепишем условия коллинеарности базис-вектора и вектора ГИС КА для случая $\beta < {{\beta }_{{\max }}}$ (два последних равенства в (3.3)) в следующем виде:

(4.2)
$\begin{gathered} {\mathbf{p}}_{{{\mathbf{V}} + i}}^{T} \cdot \left[ {{{{\mathbf{V}}}_{{\infty - i}}} \times {{{\mathbf{V}}}_{{\infty + i}}}} \right] = 0, \\ p_{{V + }}^{ \bot } = 0. \\ \end{gathered} $

В самом деле, если отлетный базис-вектор ${{{\mathbf{p}}}_{{{\mathbf{V}} + }}}$ и отлетный вектор ГИС ${\mathbf{V}}_{{\infty + }}^{{}}$ лежат в одной плоскости и при этом $p_{{V + }}^{ \bot } = 0$ (см. рис. 1), то эти векторы коллинеарны и, как было показано выше в методике расчета точки ГМ для случая β < βmax, векторы ${{{\mathbf{p}}}_{{{\mathbf{V}} - }}}$ и ${\mathbf{V}}_{{\infty - }}^{{}}$ так же становятся коллинеарны в силу соотношений (3.1) и (3.2).

Для общего случая $\beta \leqslant {{\beta }_{{\max }}}$ мы имеем 4 граничных условия в точке ГМ типа равенства – первое векторное равенство в (3.3) или в (3.7) и первое равенство в (4.2) или второе в (3.7). При этом в общем случае перпендикулярная компонента базис-вектора $p_{{V + }}^{ \bot }$ должна удовлетворять неравенству:

(4.3)
$p_{{V + }}^{ \bot } \geqslant 0.$

В итоге, в общем случае в точке ГМ мы имеем два граничных условия типа неравенства (4.1) и (4.3). Таким образом, мы получили краевую задачу с ограничениями смешанного типа (ограничениями в виде равенств и неравенств).

Воспользуемся приемом ввода в граничные условия и неизвестные параметры краевой задачи дополнительных ослабляющих переменных. Такой подход широко применяется в построении математических моделей различных экономических процессов с целью записи поставленной задачи в канонической форме основной задачи линейного или нелинейного программирования. В данных задачах ослабляющие переменные имеют простой физический смысл неиспользуемого остатка сырья того или иного вида. Так же он применяется в других областях науки, где необходимо решать оптимизационные задачи с ограничениями смешанного типа, например, в задачах планирования (организации работ) или проектирования технического объекта (системы автоматического управления, радиотехнический комплекс и т.д.) [1013].

Применим данных подход к задаче оптимизации межпланетного перелета КА. В соответствии с этим методом преобразуем граничные условия в виде неравенств в граничные условия равенства путем добавления к каждому из них неотрицательной ослабляющей переменной $b_{j}^{2}$ $(j = 1,...,2){\text{:}}$

$\begin{gathered} {{\beta }_{{\max i}}} - {{\beta }_{i}} - b_{{1i}}^{2} = 0, \\ p_{{V + }}^{ \bot } - b_{{2i}}^{2} = 0, \\ \end{gathered} $
где $i = 1,...,n,$ n – количество ГМ.

При этом согласно условиям трансверсальности или b1i = 0 или b2i = 0. В соответствии с этим запишем еще одно граничное условие: ${{b}_{{1i}}} \cdot {{b}_{{2i}}} = 0.$

В итоге в точке ГМ для случая $\beta \leqslant {{\beta }_{{\max }}}$ должны выполняться следующие граничные условия:

$\begin{gathered} {\mathbf{r}}\left( {{{t}_{i}}} \right) - {{r}_{{{\text{п л }}i}}}\left( {{{t}_{i}}} \right) = 0, \\ {\mathbf{p}}_{{_{{V + i}}}}^{T} \cdot \left[ {{{{\mathbf{V}}}_{{\infty - i}}} \times {{{\mathbf{V}}}_{{\infty + i}}}} \right] = 0, \\ {{\beta }_{{\max i}}} - {{\beta }_{i}} - b_{{1i}}^{2} = 0, \\ p_{{V + }}^{ \bot } - b_{{2i}}^{2} = 0, \\ b_{{1i}}^{2} \cdot b_{{2i}}^{2} = 0, \\ \end{gathered} $
где $i = 1,...,n$, n – количество ГМ. Всего получается 7 граничных условий. Неизвестных параметров в точке ГМ тоже 7: ${{{\mathbf{p}}}_{{\mathbf{r}}}}({{t}_{i}}),{{\beta }_{i}},{{\gamma }_{i}},{{b}_{{1i}}},{{b}_{{2i}}}.$

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

Несмотря на то, что в итоге мы получили краевую задачу с ограничениями равенствами, физически эта задача полностью соответствует краевой задаче с ограничениями смешанного типа. Такой подход оказался очень эффективным, важным является то, что он не ухудшает сходимости процесса решения краевой задачи по сравнению с решением серии краевых задач, в которых фиксировался в точках ГМ один из случаев β < βmax или β = βmax.

5. ЧИСЛЕННЫЙ ПРИМЕР МЕЖПЛАНЕТНОЙ ТРАЕКТОРИИ С ГРАВИТАЦИОННЫМ МАНЕВРОМ

Проведем сравнение работы предлагаемого подхода поиска оптимального решения, соответствующего краевой задаче с ограничениями смешанного типа с обычно используемыми методами на примерах оптимизации траектории КА с ЭРДУ с одним ГМ. Вариант № 1 соответствует случаю, когда при ГМ высота пролетной гиперболы равна минимально заданной высоте, при этом β = βmax. Вариант № 2 соответствует случаю, когда высота пролетной гиперболы не ограничена. При этом, на угол поворота ГИС β ограничений не накладывается. Постановку задачи, использующую ограничения смешанного типа, назовем вариантом № 3. Варианты расчетных схем оптимизации траектории КА представлены в табл. 1.

Таблица 1.  

Варианты расчетных схем оптимизации траекторий КА

  Вариант № 1 Вариант № 2 Вариант № 3
Условие для угла β при гравитационном маневре β = βmax β ∈ R $\beta \leqslant {{\beta }_{{\max }}}$
Условие для перпендикулярной составляющей базис-вектора при гравитационном маневре $p_{V}^{ \bot } \in R$ $p_{V}^{ \bot } = 0$ $p_{V}^{ \bot } \geqslant 0$

В качестве транспортной системы рассматривается ракета-носитель (РН) Ангара-А5, ХРБ КВТК и маршевая ЭРДУ КА, использующая ионные двигатели типа RIT-22. Ракета-носитель “Ангара-А5” выводит КА на базовую орбиту вокруг Земли. ХРБ “КВТК” переводит КА с базовой орбиты на отлетную гиперболическую траекторию. Далее ХРБ отделяется от КА и дальнейшее движение КА на гелиоцентрическом участке траектории осуществляется с помощью ЭРДУ. В качестве источника энергии используется ядерная энергоустановка. Тяга, массовый расход и удельный импульс работающего двигателя считались постоянными (не зависимыми от условий полета) на всех участках работы двигателя. Оптимизируются закон включения – выключения двигательной установки, программа управления вектором тяги, величина начального гиперболического избытка скорости, даты старта и гравитационного маневра. Общее время перелета фиксируется. Рассматривается нулевая стыковка КА с планетой назначения.

РН “Ангара-А5” обеспечивает выведение на низкую околоземную орбиту (круговая орбита высотой 200 км) КА массой 24 235 кг.

ХРБ “КВТК” в одном из вариантов имеет следующие характеристики, которые используются в настоящей работе: конечная масса ХРБ, включая массу адаптера разгонного блока с КА 3330 кг; удельная тяга двигателя 470 с.

Маршевая ЭРДУ использует двенадцать двигателей типа “RIT-22”.

Характеристики одного двигателя типа “RIT-22” приняты следующими [14]: тяга 0.2 Н; удельный импульс 4650 с; потребляемая электрическая мощность 5.785 кВт.

Общее количество двигателей “RIT-22” принималось равным 13. Минимальная высота пролетной гиперболы при ГМ принималась равной 400 км.

Рассмотрим перелет от Земли к Меркурию с ГМ у Венеры. Общее время перелета принималось равным 720 суткам. Окно старта рассматривалось в диапазоне с 1.I.2026 года по 31.XII.2026 года. Численные результаты представлены в табл. 2.

Таблица 2.  

Результаты оптимизации траектории для различных вариантов расчетных схем. Перелет Земля–Венера–Меркурий

  Вариант № 1 Вариант № 2 Вариант № 3
Дата старта 29.IV.2026 2.V.2026 29.V.2026
Юлианская дата старта 2 461 160.43 2 461 163.51 2 461 160.43
Начальный ГИС, м/с 1120.23 1229.92 1120.23
Масса КА в момент старта от Земли, кг 8549.52 8519.56 8549.52
Дата гравитационного маневра 17.III.2027 13.III.2027 17.III.2027
Юлианская дата гравитационного маневра 2 461 482.19 2 461 478.37 2 461 482.19
Величина ГИС при ГМ у Венеры, м/с 6233.66 6266.13 6233.66
Угол поворота ГИС при ГМ β, градусы 68.73 89.85 68.73
Максимальный угол поворота ГИС при ГМ βmax, градусы 68.73 68.37 68.73
Дата подлета к Меркурию, сутки 18.IV.2028 21.IV.2028 18.IV.2028
Юлианская дата подлета к Меркурию 2 461 880.43 2 461 883.51 2 461 880.43
Конечная масса КА, кг 6531.95 6597.05 6531.95

Как видно из табл. 2, оптимальным оказался случай, когда угол поворота гиперболического избытка скорости КА при гравитационном маневре равен максимальному. При этом результаты расчетов для Вариантов № 1 и 3 совпали. Таким образом, использование ограничений смешанного типа в точке гравитационного маневра позволяет получить оптимальную траекторию в рамках решения одной краевой задачи. В варианте № 2 конечная масса КА оказалась самой большой, но при этом угол поворота ГИС при ГМ оказался больше максимально допустимого, и полученная траектория не представляет практического интереса.

На рис. 3 приведена проекция на плоскость эклиптики межпланетной траектории КА Земля–Венера–Меркурий.

Рис. 3.

Проекция на плоскость эклиптики гелиоцентрической траектории перелета Земля–Венера–Меркурий.

Черным сплошным цветом показаны активные участки траектории, на которых маршевая ЭРДУ включена, черным пунктиром пассивные. Орбиты планет показаны условно круговыми серым цветом. На траектории располагаются 6 активных и 5 пассивных участков. Гелиоцентрический перелет начинается с активного участка и заканчивается пассивным при подлете к Венере для совершения ГМ, всего на этапе перелета Земля–Венера 2 активных и 2 пассивных участка. Этап перелета Венера–Меркурий начинается с пассивного движения КА при отлете от Венеры и заканчивается активным участком при подлете к Меркурию, всего на данном участке 4 активных и 4 пассивных участка.

На рис. 4 показаны зависимости функций переключения и включения-выключения ЭРДУ от времени перелета.

Рис. 4.

Функции переключения и включения-выключения двигателя вдоль гелиоцентрической траектории перелета КА при перелете Земля–Венера–Меркурий.

По рис. 4 легко определить последовательность и продолжительность активных и пассивных участков. Режим работы маршевой ЭРДУ полностью соответствует режиму работы нерегулируемого двигателя.

Рассмотрим перелет от Земли к Юпитеру с ГМ у Земли. Общее время перелета 1200 суток. Окно старта рассматривалось в диапазоне с 1.I.2018 года по 31.XII.2018 года. Численные результаты представлены в табл. 3.

Таблица 3.  

Результаты оптимизации траектории для различных вариантов расчетных схем. Перелет Земля–Земля–Юпитер

  Вариант № 1 Вариант № 2 Вариант № 3
Дата старта 22.II.2018 28.I.2018 22.I.2018
Юлианская дата старта 2 458 172.95 2 458 147.49 2 458 172.95
Начальный ГИС, м/с 963.49 1044.17 963.49
Масса КА в момент старта от Земли, кг 8587.55 8568.64 8587.55
Дата гравитационного маневра 25.II.2019 19.II.2019 25.II.2019
Юлианская дата гравитационного маневра 2 458 540.44 2 458 534.25 2 458 540.44
Величина ГИС при ГМ у Земли, м/с 8317.17 9030.94 8317.17
Угол поворота ГИС при ГМ β, градусы 54.71 74.83 54.71
Максимальный угол поворота ГИС при ГМ βmax, градусы 54.71 49.54 54.71
Дата подлета к Юпитеру, сутки 6.VI.2021 12.V.2021 6.VI.2021
Юлианская дата подлета к Юпитеру 2 459 372.95 2 459 347.49 2 459 372.95
Конечная масса КА, кг 6442.99 6526.29 6442.99

Из табл. 3 видно, что оптимальным оказался случай, когда угол поворота гиперболического избытка скорости КА при гравитационном маневре равен максимальному, так же, как и при перелете Земля–Венера–Меркурий. Результаты расчетов для вариантов № 1 и 3 совпали. Для данной схемы перелета использование ограничений смешанного типа в точке гравитационного маневра так же позволило получить оптимальную траекторию в рамках решения одной краевой задачи. В варианте № 2 конечная масса КА оказалась самой большой, при этом угол поворота ГИС при ГМ оказался больше максимально допустимого, и полученная траектория не представляет практического интереса.

На рис. 5 приведена проекция на плоскость эклиптики межпланетной траектории КА Земля–Земля–Юпитер.

Рис. 5.

Проекция на плоскость эклиптики гелиоцентрической траектории перелета Земля–Земля–Юпитер.

Черным сплошным цветом показаны активные участки траектории, на которых маршевая ЭРДУ включена, черным пунктиром пассивные. Орбиты планет показаны условно круговыми серым пунктиром. На траектории располагаются 5 активных и 4 пассивных участка. Гелиоцентрический перелет начинается с активного участка и заканчивается пассивным при подлете к Земле для совершения ГМ, всего на этапе перелета Земля–Земля 3 активных и 3 пассивных участка. Этап перелета Земля–Юпитер начинается с активного движения КА и заканчивается активным участком при подлете к Юпитеру, всего на данном этапе 2 активных участка и 1 пассивны.

На рис. 6 показаны зависимости функций переключения и включения-выключения ЭРДУ от времени перелета.

Рис. 6.

Функция переключения двигателя вдоль гелиоцентрической траектории перелета КА при перелете Земля–Земля–Юпитер.

По рис. 6 легко определить последовательность и продолжительность активных и пассивных участков. Режим работы маршевой ЭРДУ полностью соответствует режиму работы нерегулируемого двигателя.

ВЫВОД

В статье представлена новая методика расчета траекторий с гравитационными маневрами. Отличительной особенностью данной методики является сравнительно несильное увеличение порядка краевой задачи при ее применении. Данный результат достигнут за счет использования математической модели гравитационного маневра, в которой часть условий оптимальности выполняется автоматически и нет необходимости записывать их в граничные условия. Уменьшение порядка краевой задачи позволяет повысить устойчивость итерационного процесса поиска решения и сократить требования к вычислительным ресурсам, необходимым для проведения расчета. Представленная методика позволяет проводить сквозную оптимизацию траекторий с выполнением всех необходимых условий оптимальности в точке гравитационного маневра.

На основе предложенной новой методики расчета ГМ проведено объединение двух случаев, когда β < βmax и β = βmax. В результате объединения получается краевая задача с ограничениями смешанного типа, которую предлагается решать путем ввода в граничные условия дополнительных ослабляющих переменных. Такой подход позволяет получать оптимальные траектории для межпланетных перелетов КА с гравитационными маневрами в рамках решения одной краевой задачи.

Эффективность представленной методики показана на примере поиска оптимальных траекторий перелета от Земли к Меркурию с гравитационным маневром у Венеры и от Земли к Юпитеру с гравитационным маневром у Земли.

Список литературы

  1. Петухов В.Г. Оптимизация траекторий космических аппаратов с электроракетными двигательными установками методом продолжения. Дис. ... докт. тех. наук. Москва, МАИ. 2013. С. 78–93.

  2. Мин Тейн. Оптимизация траекторий космических аппаратов с использованием эволюционной стратегии с адаптацией ковариационной матрицы. Дис. ... докт. тех. наук. Москва, МАИ. 2018. С. 195–205.

  3. Strack W.C. Some numerical comparisons of three-body trajectories with patched two body trajectories for low thrust rockets. NASA Lewis Research Center. 1968. P. 18.

  4. McCarthy D.D. IERS Conventions (2003) // IERS Technical Note. IERS Conventions Centre. № 32. 2003.

  5. Souchay J. The Celestial Reference System I.C.R.S. Principles & Present Realization // IERS Technical Note. 2002. № 29. P. 115.

  6. Casalino L., Colasurdo G., Pastrone D. Optimization Procedure for Preliminary Design f Opposition-Class Mars Missions // J. Guidance, Control, and Dynamics. 1998. V. 21. № 1. P. 134–140.

  7. Casalino L., Colasurdo G., Pastrone D. Optimization of ΔV Earth-Gravity-Assist Trajectories // J. Guidance, Control, and Dynamics. 1998. V. 21. № 1. P. 991–995.

  8. Konstantinov M.S., Thein M. Method of Interplanetary Trajectory Optimization for the Spacecraft with Low Thrust and Swing-bys // Acta Astronautica. 2017. V. 136. P. 297–311.

  9. Константинов М.С., Мин Тейн. Анализ одной схемы полета космического аппарата для исследования Солнца // Труды МАИ. 2013. № 71.

  10. Акулич И.Л. Математическое программирование в примерах и задачах. М.: Высш. школа, 1986.

  11. Андрейчиков Н.П. Сборник задач по экономике, организации и планированию предприятий химической промышленности. М.: Высш. школа, 1980.

  12. Корячко В.П., Курейчик В.М., Норенков И.П. Теоретические основы САПР. М.: Энергоатомиздат, 1987.

  13. Полунин И.Ф. Курс математического программирования. М.: Высшая школа, 2008.

  14. Leiter H.J. et al. Evolution of the AIRBUS DS GmbH Radio Frequency Ion Thruster Family // Joint Conference of 30th International Symposium on Space Technology and Science 34th International Electric Propulsion Conference and 6th Nano-satellite Symposium, Hyogo-Kobe, Japan, July 4–10. 2015. P. 9.

Дополнительные материалы отсутствуют.