Оптимизация траекторий с малой тягой в переменных Кустаанхеймо–Штифеля

Abstract

В работе рассматривается регуляризация уравнений движения космического аппарата преобразованием Кустаанхеймо — Штифеля для координат и Сундмана для времени в задаче поиска оптимальной траектории межпланетного перелета с двигателем малой тяги. Из принципа максимума Понтрягина находится оптимальное управление вектором тяги при условии ограниченной мощности двигателя. Задача перелета Земля — Марс решается в регулярных переменных. Проводится сравнение найденных траекторий с траекториями, полученными методом продолжения по параметру, а также исследуется чувствительность решений краевой задачи принципа максимума в декартовых и регулярных переменных.

Full Text

ВВЕДЕНИЕ

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

К прямым методам относятся, например, метод псевдоимпульсов [1] и метод пристрелки, комбинированный с гомотопией [2, 3]. Использование прямых методов в задаче перелета Земля — Луна с малой тягой и баллистическим захватом рассматривается в работе [4]. Одним из популярных непрямых методов является принцип максимума Понтрягина [5] в сочетании с методами гомотопии. Такой подход представлен в работах В. Г. Петухова [6–9] и Д. Переза-Палау [10]. Метод гомотопического продолжения погружает краевую задачу в однопараметрическое семейство и непрерывно трансформирует начальное приближение в решение краевой задачи. Продвинутые методы гомотопии и двойной гомотопии рассматриваются в работах Б. Пана [11, 12], Ф. Джианга [13] и Ч. Джанга [14].

Для увеличения устойчивости непрямых методов особое внимание стоит уделить выбору фазового вектора. Возможные способы выбора переменных состояния рассматриваются в работах Э. Тахери [15–18]. Сравнительное исследование различных вариантов выполняется в работах Дж. Джанкинса [19] и С. Жефруа [20]. Зачастую замена вектора состояния сопровождается репараметризацией траектории, то есть заменой физического времени фиктивным. Первым, кто предложил такой подход был К. Ф. Сундман [21]. В случае репараметризации физическое время может стать компонентой вектора состояния или может быть заменено на временную компоненту. Временные компоненты изучаются в работах П. Накози [22] и Е. Брумберга [23].

Уравнения орбитального движения в декартовых переменных не являются регулярными: правая часть уравнений содержит сингулярность в начале координат. Приведение системы к регулярному виду называется регуляризацией. Один из методов регуляризации — преобразование Кустаанхеймо — Штифеля (KS) [24]. KS-преобразование осуществляет переход от трехмерного пространства к четырехмерному и представляет собой обобщение преобразования Леви — Чивиты [25]. Такое преобразование обсуждается в статье [26]. Стоит учитывать, что условие на физическое существование траекторий в пространстве KS-переменных выражается через билинейное соотношение. Использование KS-переменных вместе с преобразованием Сундмана для задачи оптимального перелета рассматривается в публикации [27], где метод продолжения по параметру применяется следующим образом: сначала траектория в KS-переменных рассчитывается без соблюдения билинейного соотношения, после чего методом продолжения приводится к решению с его соблюдением.

Свойство KS-переменных регуляризировать уравнения движения оказывается полезно в различных задачах небесной механики и астродинамики, где наблюдается высокая чувствительность траекторий. Скажем, в настоящее время KS-переменные находят применение в области предсказания и коррекции движения КА [28] и в моделировании методом Монте-Карло траекторий небесных тел в задачах планетарной защиты [29]. Эти переменные также используются при исследовании топологической устойчивости в задаче N тел [30].

Целью настоящего исследования является изучение потенциальной выгоды от применения KS-преобразования совместно с преобразованием Сундмана к задаче оптимального межпланетного перелета с малой тягой. Этот подход будет сравниваться с гомотопическим методом продолжения по параметру в декартовых переменных. Особое внимание уделяется анализу обусловленности матрицы чувствительности. Билинейное соотношение, выражающее физичность траектории в KS-пространстве, будет учтено в качестве смешанного ограничения в функции Гамильтона — Понтрягина. Терминальные многообразия в KS-переменных выписываются явным образом.

В работе формулируется задача оптимального перелета с ограниченным по мощности идеально регулируемым двигателем малой тяги. После этого вводятся KS-преобразование и преобразование Сундмана. Далее описывается применение принципа максимума Понтрягина в KS-переменных для сведения задачи к двухточечной краевой, которая затем сводится к задаче оптимизации и решается с помощью метода внутренней точки [31]. В качестве модельной задачи рассматривается перелет Земля — Марс. Сравниваются оптимальные решения этой задачи в KS-переменных и в декартовых переменных.

УРАВНЕНИЯ ДВИЖЕНИЯ

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

$$\begin{gathered} \dot{r}=v, \\ \dot{v}=-\frac{\mu }{r^3}r+a, \end{gathered}$$ (1)

где r и v — это векторы положения и скорости; a — вектор реактивного ускорения и μ — гравитационный параметр центрального тела. Реактивное ускорение выражается как

$a=\frac{F_T}{m\left(t\right)}=\frac{{\dot{m}}_{s}u}{m\left(t\right)}e$, (2)

где ${\dot{m}}_s$ — массовый расход (масса рабочего тела, затрачиваемого в единицу времени); m –масса космического аппарата в момент времени t; e — вектор направления тяги; и u — скорость истечения рабочего тела. Если мощность двигателя ограничена доступной электрической мощностью N_max с известным КПД η, тогда

$\frac{\dot{m}{u}^2}{2}\le \eta N_{\max }$. (3)

Предполагается, что двигатель идеально регулируемый: можно выбрать u и ${\dot{m}}_s$ произвольным образом, лишь бы выполнялось условие. Очевидно, что при этом величина тяги может принимать любые значения. Выражая из равенства u через a, m, ${\dot{m}}_s$ и подставляя это выражение в формулу, получим

$\frac{a^2\left(t\right)m^2\left(t\right)}{2\dot{m}}\le \eta N_{\max }$. (4)

После интегрирования неравенство принимает вид

$\frac{1}{m_0-m\left(t\right)}-\frac{1}{m_0}\le {\int }_{t_0}^t\frac{a^2\left(\tau \right)d\tau }{2\eta N_{\max }}$. (5)

Минимум потребления рабочего тела достигается при использовании всей доступной электрической мощности, то есть превращении неравенства (5) в равенство. Предполагается, что доступная мощность не меняется со временем, что соответствует, например, случаю использования ядерных бортовых источников энергии постоянной мощности. Минимальное потребление топлива требует минимизации функционала

$J=\frac{1}{2\eta N_{\max }}{\int }_{t_0}^{t_f}{a}^2\left(t\right)dt$. (6)

Преобразование времени Сундмана

Одним из важных методов регуляризации является замена независимой переменной [32]. В этом случае траектория репараметризуется: вместо физического времени t вводится новая переменная s, называемая фиктивным временем. Репараметризация записывается как

$dt=c(r,\dot{r})ds$. (7)

Можно задать масштабирующий множитель c так, чтобы правая часть уравнений движения менялась практически равномерно. Сундман предложил выбирать масштабирующий множитель в виде dt = rds. Если поставить требование для фиктивного времени, чтобы оно менялось на 2π за орбитальный период, масштабирующий множитель немного модифицируется:

$dt=\frac{r}{\sqrt{-2h}}ds$. (8)

Здесь r — расстояние до центрального тела, а h — кеплерова энергия. Можно показать, что в этой нормированной версии преобразования Сундмана роль фиктивного времени s играет эксцентрическая аномалия E. Для отслеживания физического времени необходимо добавить к системе временное уравнение

$t^{\text{'}}=\frac{r}{\sqrt{-2h}}$, (9)

где штрих обозначает производную по фиктивному времени. Однако временное уравнение уменьшает эффект регуляризации из-за короткопериодических колебаний масштабирующего множителя. Решением этой проблемы будет использование временного элемента — переменной, которая растет линейно по отношению к фиктивному времени в невозмущенном случае, но в то же время явным образом связана с физическим временем. Временной элемент определяется как

$\tau =t_0+\sqrt{\frac{{\mu }^2}{8h^3}}E$. (10)

Как видно, временной элемент τ линеен по отношению к эксцентрической аномалии E.

Преобразование Кустаанхеймо — Штифеля

Преобразование Кустаанхеймо — Штифеля — это гладкое вложение трехмерного конфигурационного пространства (x, y, z)^T в четырехмерное параметрическое пространство (u₁, u₂, u₃, u₄)^T. Это отображение записывается в матричной форме как

$\left(\begin{array}{c}x\\ y\\ z\\ 0\end{array}\right)=\left(\begin{array}{cccc}{u}_1& -u_2& -u_3& u_4\\ u_2& u_1& -u_4& -u_3\\ u_3& u_4& u_1& u_2\\ u_4& -u_3& u_2& -u_1\end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right)=L\left(u\right)u$, (11)

или же, в краткой форме,

$r=L\left(u\right)u$. (12)

Четвертая компонента радиус-вектора обнуляется при условии

$u_b\cdot u\text{'}=0$, (13)

где u_b — это четвертая строка матрицы L. Для того, чтобы движение в параметрическом пространстве было физичным (соответствующим движению в физическом пространстве), требуется выполнение билинейного соотношения

$u_b\cdot u\text{'}=0$, (14)

в каждый момент времени. Четырехмерное параметрическое пространство расслаивается: каждый слой

$\widehat{u}=\left\{u:L\left(u\right)u=r,r\in {ℝ}^3\right\}$. (15)

представляет собой одномерное многообразие точек KS-пространства, которые отображаются в одну точку в физическом пространстве.

Уравнения движения в KS-переменных

Преобразование Сундмана в KS-переменных принимает вид

$t^{\text{'}}=\frac{u\cdot u}{\sqrt{-2h}}$, (16)

где кеплерова энергия выражается как

$h=\frac{{\nu }^2}{2}-\frac{\mu }{r}=-\frac{\mu }{u\cdot u+4(w\cdot w)}$. (17)

В этом уравнении w — параметрическая скорость. После применения преобразования Сундмана и KS-преобразования уравнения движения могут быть записаны как

$u^{\text{'}\text{'}}+\frac{u}{4}=-\frac{u\cdot u}{4h}{L}^T\left(u\right)a-\frac{h^{\text{'}}}{2h}{u}^{\text{'}}$, (18)

где производная кеплеровой энергии по фиктивному времени имеет вид

$h^{\text{'}}=2w\cdot L^{T}a$. (19)

Для использования KS-переменных расширим вектор реактивного ускорения до четырехмерного $a\in {ℝ}^4$. Позже будет наложено ограничение для его физической осуществимости. Также обратим внимание, что дифференциальные уравнения задачи двух тел принимают форму регулярных уравнений возмущенного осциллятора. В отсутствие возмущений существует точное решение: гармонические колебания с частотой 1/2. Вывод уравнения (18) детально представлен в работе [24].

Для заданной точки в KS-пространстве параметрическая скорость w связана с производной радиус-вектора по фиктивному времени r' = vt' согласно формуле

$w=\frac{L^T\left(u\right)r^{\text{'}}}{2(u\cdot u)}$. (20)

Это выражение позволяет вычислять начальное значение скорости, необходимое для численного интегрирования. Многообразие параметрических скоростей, которые соответствуют некоторой физической скорости, выражается схожим с формулой (15) образом:

$\widehat{w}=\left\{w:\frac{2\sqrt{-2h}L\left(u\right)w}{(u\cdot u)}=v,v\in {ℝ}^3\right\}$. (21)

Наконец, выражение для временного элемента τ в параметрическом пространстве записывается как

$\tau =t_0+2\frac{u\cdot w}{\sqrt{-2h}}$. (22)

Дифференциальное уравнение на τ принимает вид

${\tau }^{\text{'}}=\frac{\left[\mu +4(u\cdot w)h^{\text{'}}+(u\cdot u)u\cdot L^{T}a\right]}{{(-2h)}^{3/2}}$. (23)

Итого для описания движения КА в параметрическом пространстве требуется 9 скалярных дифференциальных уравнений: к уравнению (23) добавляются еще восемь из векторного уравнения второго порядка (18).

ОПТИМИЗАЦИЯ ТРАЕКТОРИИ В РЕГУЛЯРНЫХ ПЕРЕМЕННЫХ

В этом разделе задача оптимального управления в KS-переменных решается с помощью принципа максимума Понтрягина. Условия трансверсальности выражаются в параметрическом пространстве. Получившаяся двухточечная краевая задача преобразуется в задачу нелинейного программирования с граничными условиями и затем решается с помощью метода внутренней точки.

Функционал стоимости, дифференциальные и смешанные ограничения

Система дифференциальных уравнений орбитального движения в KS-переменных записывается как

$$\begin{gathered} u^{\text{'}}=w, \\ {w}^{\text{'}}=-\frac{u}{4}-\frac{u^2}{4h}{L}^{T}a-\frac{2w\cdot L^{T}a}{2h}w, \\ {\tau }^{\text{'}}=\frac{\mu +\left(8\right(u\cdot w)w+(u\cdot u\left)u\right)\cdot L^{T}a}{{(-2h)}^{3/2}}. \end{gathered}$$ (24)

Функционал стоимости также должен быть выражен через KS-переменные и фиктивное время:

$J=\frac{1}{2\eta N_{\max }}{\int }_{s_0}^{s_f}{a}^2\left(t\right)\frac{u^2\left(s\right)}{\sqrt{-2h\left(s\right)}}ds$. (25)

Физическая осуществимость управления требует, чтобы четвертая компонента вектора управления тождественно равнялась нулю:

$g(u,a)=u_b\cdot L^{T}a=0$. (26)

Этот тип ограничений называется смешанным, так как включает в себя и вектор состояния, и вектор управления. Несложно показать, что в силу

${\left(u_b\cdot u^{\text{'}}\right)}^{\text{'}}=u_b\cdot u^{\text{'}\text{'}}+{u^{\text{'}}}_b\cdot u^{\text{'}}=u_b\cdot w^{\text{'}}$ (27)

функция g пропорциональна производной билинейного соотношения:

$u_b\cdot w^{\text{'}}=-\frac{u^2}{4h}g(u,a)$. (28)

В общем случае в условиях возмущенной динамики при численном интегрировании траектории часто проявляется топологическая неустойчивость [30], что может привести к нарушению ее физической осуществимости. Явное добавление в систему смешанного ограничения, в котором вектор-функция g пропорциональна скорости изменения билинейного параметра, представляет собой, по сути, способ контроля локальной точности интегрирования. Такой подход к стабилизации системы с помощью добавления в нее производных от интегральных соотношений был предложен еще Штифелем [24].

Принцип максимума в регулярных переменных

Теперь возможно применить принцип максимума Понтрягина к автономной системе (24) с функционалом стоимости (25) и смешанным ограничением (26). Процедура учета смешанного ограничения описана, к примеру, в работе [33]: оно должно быть добавлено в функцию Гамильтона — Понтрягина с множителем Лагранжа β:

$H=-\varphi (x,a)+p\cdot f(s,x)-\beta g(s,x,a)$. (29)

Здесь φ — это подынтегральная функция в функционале стоимости; p — это вектор сопряженных переменных; и f — функция правых частей системы. Значения β и a могут быть определены из условия стационарности гамильтониана и смешанного ограничения:

$$\begin{gathered} H_a(x,a)=0, \\ g(x,a)=0. \end{gathered}$$ (30)

Здесь введено обозначение $H_a=\partial H/\partial a$. Подстановка подынтегральной функции и функции правых частей системы в гамильтониан дает

$$\begin{gathered} H=-\frac{u^2}{\sqrt{-2h}}\frac{a^2}{2}+p_u\cdot w+p_w\cdot \left(-\frac{u}{4}-\frac{u^2}{4h}{L}^{T}a-\frac{w}{h}w\cdot L^{T}a\right)+ \\ +p_{\tau }{\tau }^{\text{'}}-\beta u_b\cdot L^{T}a. \end{gathered}$$ (31)

Введем вспомогательный вектор

$\Lambda =L\left(-\frac{u^2}{4h}{p}_w-\frac{w}{h}{p}_w\cdot w\right)$, (32)

который упрощает выражение до

$H=-\frac{u^2}{\sqrt{-2h}}\frac{a^2}{2}+\Lambda \cdot a+p_u\cdot w-\frac{p_w\cdot u}{4}+p_{\tau }{\tau }^{\text{'}}-\beta u_b\cdot L^{T}a$. (33)

Также учтем тот факт, что система (24) не зависит от временного элемента τ. Это приводит к тому, что сопряженная переменная p_τ является константой вследствие выполнения $p_{\tau }^{\text{'}}=0$. Забегая немного вперед, отметим, что ограничений на τ в конечный момент времени s_f не накладывается, а значит, из условий трансверсальности следует p_τ(s_f) = 0, что позволяет исключить слагаемое p_ττ' из рассмотрения. Условие стационарности, также известное как уравнение Эйлера–Лагранжа, записывается как

$H_a=-\frac{u^2}{\sqrt{-2h}}a+\Lambda -\beta L{u}_b=0$. (34)

Если уравнение (34) умножить на u_bL^T слева, а потом использовать условие (26) для сокращения первого слагаемого, то можно выразить множитель Лагранжа:

$\beta =\frac{{\Lambda }_4}{u^2}$. (35)

После подстановки в (34) удобно ввести обозначение

$\tilde{\Lambda }=\Lambda -\frac{{\Lambda }_4}{u^2}L{u}_b$. (36)

Тогда для оптимального управления имеем:

$a=\frac{\sqrt{-2h}}{u^2}\tilde{\Lambda }$. (37)

Фактически это выражение влечет выполнение смешанного ограничения (26): четвертая компонента вектора управления всегда будет равна нулю. Подстановка оптимального управления в гамильтониан дает

$H^{*}=\frac{1}{2}\frac{\sqrt{-2h}}{u^2}{\tilde{\Lambda }}^2+p_u\cdot w-\frac{p_w\cdot u}{4}$. (38)

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

$$\begin{gathered} u^{\text{'}}=\frac{dt}{ds}\frac{\partial }{\partial p_u}\left(\frac{ds}{dt}{H}^{*}\right), \\ {w}^{\text{'}}=\frac{dt}{ds}\frac{\partial }{\partial p_w}\left(\frac{ds}{dt}{H}^{*}\right), \\ {p}_u^{\text{'}}=-\frac{dt}{ds}\frac{\partial }{\partial u}\left(\frac{ds}{dt}{H}^{*}\right), \\ {p}_w^{\text{'}}=-\frac{dt}{ds}\frac{\partial }{\partial w}\left(\frac{ds}{dt}{H}^{*}\right). \end{gathered}$$ (39)

Здесь выражение H* ds/dt представляет собой, с точностью до константы, гамильтониан в физическом времени $\tilde{H}$ [34], а множитель dt/ds появляется при переходе от производных по физическому времени к производным по фиктивному времени. Другими словами, в физическом времени уравнения (39) имели бы следующий вид:

$$\begin{gathered} \dot{u}=\frac{\partial \tilde{H}}{\partial p_u}, \\ \dot{w}=\frac{\partial \tilde{H}}{\partial p_w}, \\ {\dot{p}}_u=-\frac{\partial \tilde{H}}{\partial u}, \\ {\dot{p}}_w=-\frac{\partial \tilde{H}}{\partial w}. \end{gathered}$$ (40)

Условия трансверсальности

Чтобы определить двухточечную краевую задачу, нужно записать условия трансверсальности. В общем случае они имеют вид

$p\left(s_{\alpha }\right)={\left(\frac{\partial g_{\alpha }\left(x\right(s_{\alpha }\left)\right)}{\partial x\left(s_{\alpha }\right)}\right)}^T\cdot C$, (41)

где индекс α принимает значения 0 или f; g_α определяет соответствующее терминальное многообразие; x = (u, w)^T — это параметрический вектор состояния; p = (p_u, p_w)^T — сопряженный вектор и C — вектор произвольных постоянных.

Будем рассматривать задачу встречи, когда КА должен достичь требуемых положения r(s_f) и скорости v(s_f) в конечный момент времени. Для этой задачи в KS-переменных терминальное многообразие в точке s_f выражается как $u_f\in \widehat{u}\left(s_f\right)$ и $w_f\in \widehat{w}\left(s_f\right)$. Это дает шесть уравнений для g_f:

$$\begin{gathered} r\left(s_f\right)-\tilde{L}\left(u_f\right)u_f=\mathrm{0,} \\ v\left(s_f\right)-\frac{2\sqrt{-2h_f}\tilde{L}\left(u_f\right)w_f}{u_f\cdot u_f}=0. \end{gathered}$$ (42)

Здесь $\tilde{L}$ — это первые три строки L. Вектор x₀ известен на определенную дату старта t(s₀). Он должен удовлетворять условиям $u_0\in \widehat{u}\left(s_0\right)$ и $w_0\in \widehat{w}\left(s_0\right)$, но в остальном может быть выбран произвольным образом и зафиксирован для дальнейших вычислений. Дополнительные ограничения для u(s_f) не требуются, потому что условие (13) выполняется в каждой точке траектории, если оно выполняется в момент s₀. Значения r(s_f) и v(s_f) определяются как положение и скорость небесного тела в момент времени $t_f=t({\tau }_f,x_f)=t(s_f,x_0,p_0)$, где конкретное значение t_f зависит не только от конечного фиктивного времени, но и от самой траектории КА.

Дифференцирование выражений (42) по x дает матрицу 6 × 8. Следовательно, размер вектора C равен 6. Это значит, что вектор p на концах траектории должен принадлежать шестимерному линейному подпространству восьмимерного пространства сопряженных переменных. Двумерное ортогональное дополнение этого подпространства — это линейная оболочка двух векторов, которые могут быть записаны в виде матрицы

$K\left(s_{\alpha }\right)=\left(\begin{array}{cc}{u}_1{u}_4& u_2{u}_4\\ -u_1{u}_3& -u_2{u}_3\\ u_1{u}_2& u_2^2\\ -u_1^2& -u_1{u}_2\\ u_1{w}_4-u_4{w}_1& u_2{w}_4-u_4{w}_2\\ u_3{w}_1-u_1{w}_3& u_3{w}_2-u_2{w}_3\\ u_1{w}_2-u_2{w}_1& 0\\ 0& u_1{w}_2-u_2{w}_1\end{array}\right)$. (43)

Эти векторы могут быть найдены как фундаментальные решения линейной системы

${\left(\frac{\partial g_{\alpha }\left(x\right(s_{\alpha }\left)\right)}{\partial x\left(s_{\alpha }\right)}\right)}^T\cdot C=0$. (44)

Суммируя вышесказанное, условия трансверсальности в случае задачи встречи могут быть записаны как

$K^T\left(s_{\alpha }\right)p\left(s_{\alpha }\right)=0$. (45)

Между граничными условиями в начальной и конечной точке существует разница. В начальной точке t = t₀, декартовы положение r₀ и скорость v₀ известны. Например, они могут быть равны положению и скорости небесного тела в точке вылета, которые могут быть взяты из эфемерид. Точное значение u₀ неважно, оно может быть выбрано из ${\widehat{u}}_0$ произвольно. После этого параметрическая скорость w₀ вычисляется согласно уравнению (20). Наконец, условия на сопряженные переменные выражаются как.

Граничные условия в конечной точке более сложны: в случае фиксирования фиктивного времени продолжительность полета в физическом времени становится функцией от всей траектории. Тогда, параметрические векторы положения и скорости u_f и w_f каждый принадлежат двумерному тору — семейству одномерных слоев, которое параметризуется длительностью перелета в физическом времени.

Постановка задачи оптимизации

Краевая задача может быть преобразована в задачу оптимизации. Во-первых, необходимо записать усеченное KS-преобразование:

$$\begin{gathered} r\left(s_f\right)=\tilde{L}\left(u_f\right)u_f, \\ v\left(s_f\right)=\frac{2\sqrt{-2h_f}\tilde{L}{w}_f}{u_f\cdot u_f}. \end{gathered}$$ (46)

Вектор физического положения в конечной точке зависит от начального состояния и сопряженных переменных $r_f=r\left(t_f\right(s_f,x_0,p_0\left)\right)$. Аналогичные зависимости будут соблюдаться для физической скорости: $v_f=v\left(t_f\right(s_f,x_0,p_0\left)\right)$. В задаче встречи вектор-функция невязки принимает вид

$\delta =\left(\begin{array}{c}r\left(s_f\right)-\tilde{L}\left(u_f\right)u_f\\ v\left(s_f\right)-\frac{2\sqrt{-2h_f}\tilde{L}\left(u_f\right)w_f}{(u_f\cdot u_f)}\\ K^T\left(s_0\right)p\left(s_0\right)\end{array}\right)$, (47)

где начальные значения r₀ и v₀ уже известны. Значение функции невязки δ — это восьмимерный вектор. Требования δ_i = 0, i = 1...8 можно трактовать как ограничения типа равенства в задаче оптимизации с тождественно постоянным функционалом стоимости. Такой формальный трюк позволяет свести нашу краевую задачу к задаче конечномерной оптимизации и применить эффективные оптимизационные методы. В частности, в данной работе использовался метод внутренней точки [31]. Этот метод сводит исходную задачу к задаче выпуклой оптимизации, сначала линеаризуя ограничения типа равенства, а потом преобразуя их в выпуклую барьерную функцию. В случае тождественно постоянного функционала минимизируется только барьерная функция, что в свою очередь приводит к минимизации невязки (47). В данной задаче метод внутренней точки оказался более эффективен, чем метод последовательного квадратичного программирования [31] или метод Левенберга — Марквардта [31].

Функция невязки зависит от вектора сопряженных переменных в начальной точке и времени интегрирования s_f, то есть δ = δ(s_f, p₀). Также, метод внутренней точки требует наложить ограничения на сопряженный вектор p₀ ∈ [p_a, p_b], то есть задать область поиска решений. Начальное приближение для вектора сопряженных переменных выбирается равным нулю.

Определение параметра s_f является ключевым для успешного решения задачи. Для рассматриваемых траекторий можно записать выражение $n+\Delta E-s_f=\epsilon (x_0,p_0)$, где n — число витков траектории; ΔE — изменение эксцентрической аномалии в результате перелета, а ε(x₀, p₀) — некое малое отклонение, которое зависит от траектории КА. Выбирая параметры n и ΔE, можно контролировать угловую дальность перелета напрямую. При отсутствии априорного знания о функции рекомендуется использовать значение $s_f=n+\Delta E$.

Улучшение сходимости может быть достигнуто при помощи преобразования терминального многообразия в задачу со свободным концом

$s_f\in [s_a,s_b]$, (48)

где $s_a,s_b=s_f^0\pm C$. Константа $C\in ℝ, C\ge 0$ может быть выбрана в широком диапазоне значений, но рекомендуется использовать C ≤ π / 2. Эта модификация позволяет методу сходиться быстрее, что дает начальное приближение для p₀. Имея некоторую дополнительную информацию об оптимальных траекториях, можно сузить диапазон значений фиктивного времени. Например, для задач перелета к внешним планетам Солнечной системы (более удаленным от Солнца, чем Земля) рекомендуется брать $s_a=s_f^0, s_b>s_f^0$.

МОДЕЛЬНАЯ ЗАДАЧА: ПЕРЕЛЁТ ЗЕМЛЯ — МАРС

Этот раздел описывает численные эксперименты с предложенным подходом и его параметрами. Кроме того, проведено сравнение с методом продолжения по параметру [8]. Сравниваемыми характеристиками являются сходимость и численная устойчивость методов.

Параметры численных расчетов

Алгоритм реализован в среде MATLAB, и запускался на ноутбуке со следующими характеристиками: 12 виртуальных ядер с частотой 2.6 ГГц и объемом RAM 32 Гб. Выражения для правых частей уравнений, так же, как и их частные производные, вычислены символьно и затем оптимизированы для уменьшения времени вычислений метода. Для решения оптимизационной задачи использовалась функция fmincon, которая реализует метод внутренней точки [35]. Это метод распараллеливается на 6 процессов, одновременно вычисляя барьерную функцию в нескольких пробных точках. Точность выполнения условий оптимальности составляла 10^–10, точность соблюдения граничных условий 10^–12, максимальное число вызовов функции взято равным 10¹⁰, размер шага был ограничен снизу величиной 10^–10, а максимальное число шагов оптимизации равнялось 250. Дополнительно добавлено условие τ' ≥ 0.

Для нахождения положений и скоростей Земли и Марса в начальный и конечный момент времени использовалась эфемеридная модель DE430 [36]. Начальный момент времени соответствует 00:00 UTC1.I.2022.

Решение оптимизационной задачи состоит из двух этапов. Для получения решения время интегрирования s_f принадлежит определенному диапазону во время первого этапа и фиксируется во время второго этапа. При этом физическое время t_f вычисляется как функция от s_f и траектории КА. Стоит заметить, что все переменные предварительно обезразмериваются, а после нахождения решений их размерность восстанавливается. Кроме того, для вектора p₀ установлены граничные условия вида |p_i,0| '' 1. Также вектор p₀ масштабируется путем домножения на гиперпараметр $C_p={10}^{-4-\sqrt{s_f}}$, который явно зависит от времени интегрирования s_f. Это масштабирование выполняется при подаче p₀ в оптимизатор, чтобы избежать слишком малых значений, однако при интегрировании траекторий КА используется не отмасштабированное значение. Здесь и далее под гиперпараметрами будем иметь в виду некоторые параметры настройки метода, выбираемые исследователем, например, точность интегрирования. Параметр C_p был подобран, исходя из типичных величин вектора сопряженных переменных p₀ для выбранного фиктивного времени.

Рассматривалась задача встречи. На первом этапе минимизировалась невязка (47) на области поиска [s_f, s_f + π / 8] с нулевым начальным приближением p₀. Полученные решения на этом этапе могут являться как локально оптимальными, так и глобально оптимальными. Для определения глобальной оптимальности необходимо выяснить, принадлежит ли данное решение Парето-фронту. На втором этапе в качестве начальных приближений использовались p₀ из первого этапа, а область поиска сужалась до точки [s_f, s_f]. Выбирались минимальные по расходу топлива части семейств, а потом проводилась оптимизация вдоль этих семейств влево и вправо. В обоих случаях невязка домножалась на гиперпараметр C_δ = 10¹⁰. Помимо этого, точность интегрирования контролировалась третьим гиперпараметром: C_tol = 10⁻¹⁰ на первом этапе и C_tol = 10⁻¹² на втором.

После получения оптимального решения находится t_f. Далее применялся метод продолжения с модификацией гравитационного параметра в физических координатах для нахождения решения, которое будем называть эталонным. Затем оба подхода сравнивались с точки зрения затрат рабочего тела и геометрии оптимальных решений. Отличие между ними оценивалось по формуле

$d=\underset{r_0\in {\gamma }_0}{\max }\left(\underset{r_1\in {\gamma }_1}{\min }||r_0-r_1||\right)$, (49)

где γ₀ и γ₁ — две траектории, полученные разными подходами.

В рамках численного эксперимента были выбраны точки s_f = (3/4 + k/16)π, где k = 0...168. Другими словами, рассматривались перелеты с угловой дальностью от трех четвертей витка до шести витков.

Были выбраны следующие параметры КА и двигательной установки, аналогичные параметрам аппарата Европейского космического агентства SMART-1 (запущен к Луне в 2003 г.) и его двигателя PPS-1350 [37]. Начальная масса КА составляла m₀ = 367 кг, а максимальная мощность, доступная двигательной установке, равнялась N = 1350 Вт. Для этого двигателя были получены оценки уровня тяги 73 мН и скорости истечения 16.4 км/с. Судя по этим значениям, эффективный КПД ЭРДУ достигал η = 0.45.

Сравнение двух методов

Части нескольких семейств решений, обеспечивающие для фиксированного времени полета минимум затрат топлива, можно объединить в Парето-фронт в координатах фиктивное время — затраты топлива. Вместо фиктивного времени более наглядно ввести угловую дальность перелета, скажем, в терминах изменения эксцентрической аномалии. На рис. 1 представлено сравнение решений, получаемых двумя подходами, в координатах угловая дальность — затраты топлива. Для сохранения представления о поведении семейств решений на графике оставлены и те их участки, которые не принадлежат Парето-фронту; они показаны штриховыми линиями. Сплошной линией обозначены решения, полученные методом внутренней точки и образующие Парето-фронт, а крестиками — решения, полученные стандартным методом продолжения в декартовых координатах.

Рис. 1. Разные семейства решений и Парето-фронт, восстановленный двумя разными подходами

Среднее значение разницы в затратах топлива составляет 10^–6 кг, а среднее расстояние между траекториями равно 90 км. Такой результат показывает хорошее совпадение между решениями, полученными разными методами. Заметим, что значения угловой дальности перелета (изменения эксцентрической аномалии) оказываются очень близкими к значениям длительности перелета в фиктивном времени. На рис. 2 показана их разность в точках Парето-фронта.

Рис. 2. Разница между фиктивным временем перелета и изменением эксцентрической аномалии для различных Парето-оптимальных решений

На рис. 3 представлено время вычисления решений с разной угловой дальностью двумя рассмотренными методами. Среднее время поиска решения составляет 16 с для метода продолжения и 3 с для метода внутренней точки.

Рис. 3. Время вычисления локально оптимальных решений разными подходами

Сходимость и гиперпараметры

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

Теперь перейдем к смыслу и проблеме выбора значений гиперпараметров. Множитель невязки C_δ позволяет избежать слишком малых ее значений, которые могут достигать машинного нуля, что приводит к ранней остановке оптимизации. На основе ряда численных экспериментов было выбрано значение 10¹⁰.

Точность интегрирования C_tol особенно важна на первом этапе поиска решения и в меньшей степени на втором. При необходимости ее можно уменьшать до 10^–10 на первом этапе, и это повлияет на близость начального приближения второй стадии к оптимальному.

Третий гиперпараметр C_p — это множитель начальных значений сопряженных переменных, к значению которого задача оказывается наиболее чувствительна. Несмотря на отсутствие влияния на получающиеся траектории, этот параметр влияет на процедуру оптимизации. Меньшие его значения будут неявным образом уменьшать шаг оптимизатора, но увеличивать время сходимости. В идеальном случае он должен быть таким, чтобы отмасштабированные значения p₀, подаваемые в качестве аргумента оптимизатора fmincon, были близки к единице.

Устойчивость траекторий

Матрица чувствительности

$A=\frac{\partial \delta }{\partial p_0}$ (50)

может быть использована для оценки устойчивости решений. Это матрица вычисляется с помощью центральной разности при интегрировании траектории. Точность вычисления матрицы чувствительности приблизительно равна 10^–7. Это значение также было определено численно. Чувствительность решений к начальным условиям характеризуется числом обусловленности матрицы чувствительности. Для квадратных матриц это число определяется как

$cond\left(A\right)={||A||}_2\cdot {||A^{-1}||}_2$. (51)

Оно равно соотношению между максимальным и минимальным сингулярными значениями матрицы A. На рис. 4 показаны числа обусловленности для семейства оптимальных траекторий перелета с различной угловой дальностью, полученных двумя подходами. Среднее значение чисел обусловленности для решений, найденных методом продолжения в декартовых координатах, равно 1.3·10⁵, в то время как для решений, полученных методом внутренней точки в KS-переменных, число обусловленности в среднем составляет 4.3·10³, что означает лучшую обусловленность матрицы чувствительности.

Рис. 4. Числа обусловленности матрицы чувствительности для двух разных подходов

ЗАКЛЮЧЕНИЕ

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

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

Метод внутренней точки, используемый для решения задачи оптимизации, к которой сводится краевая задача принципа максимума Понтрягина в KS-переменных, хорошо распараллеливается. Это помогает уменьшить время вычислений на поиск оптимальных траекторий.

В дальнейшем планируется протестировать эффективность разработанного подхода в случае более реалистичных моделей среды и двигателя (то есть, с другими видами функционала стоимости и в присутствии внешних возмущений). Целесообразно провести сравнение KS-переменных с другими наборами фазовых переменных в смысле численных свойств получающихся краевых и оптимизационных задач.

ФИНАНСИРОВАНИЕ РАБОТЫ

Исследование выполнено за счет гранта Российского научного фонда (проект № 19-11-00256).