Космические исследования, 2021, T. 59, № 4, стр. 327-338

Решение задачи оптимального вывода на орбиту космического аппарата с использованием реактивного ускорения и солнечного паруса в переменных Кустаанхеймо–Штифеля

Я. Г. Сапунков 1, Ю. Н. Челноков 1*

1 Институт проблем точной механики и управления РАН
Саратов, Россия

* E-mail: ChelnokovYuN@gmail.com

Поступила в редакцию 16.04.2019
После доработки 05.08.2020
Принята к публикации 17.09.2020

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

Аннотация

Решена с помощью принципа максимума Понтрягина и переменных Кустаанхеймо–Штифеля пространственная задача оптимального вывода на заданную орбиту космического аппарата, управляемого с помощью солнечного паруса и ограниченного или импульсного реактивного ускорения центра масс КА. Минимизируемый функционал представляет собой линейную комбинацию с весовыми множителями двух критериев: времени и интегральной суммы величин импульсов реактивного ускорения центра масс КА, затраченных на процесс управления. Приведены первые интегралы уравнений задачи и формулы для определения приращений фазовых и сопряженных переменных под действием сообщаемого импульса реактивного ускорения. Получены численные решения задачи для ограниченного или импульсного ускорения при наличии солнечного паруса или при его отсутствии. Дана оценка влияния наличия солнечного паруса на длительность процесса, на суммарный импульс величины реактивного ускорения и на величину минимизируемого функционала. Рассмотрены случаи коррекции орбиты и случаи, когда элементы новой орбиты существенно отличаются от элементов начальной орбиты КА. Установлена оптимальность реактивного ускорения, ортогонального к плоскости оскулирующей орбиты КА, для рассмотренных малых значений отклонений угловых элементов орбиты от их начальных значений, т.е. оптимальность такого ускорения в рассмотренных примерах задачи коррекции угловых элементов орбиты КА.

ВВЕДЕНИЕ

Задачи об оптимальном управлении движением центра КА с использованием переменных Кустаанхеймо–Штифеля [1] и кватернионных переменных рассматривались в работах [2, 3]. Задачи об оптимальном управлении движением КА с использованием солнечного паруса и различных функционалов качества, с использованием прямоугольных или полярных координат, кватернионных переменных, переменных Кустаанхеймо–Штифеля, кватернионных элементов орбиты рассматривались в работах [410]. Решались задачи об оптимальном управлении для встречи двух КА, задачи быстродействия, а также задачи с комбинированными функционалами качества процесса управления. Солнечный парус использует давление солнечного света на отражающую поверхность для приведения в движение КА [5]. Использование солнечного паруса позволяет уменьшить массу КА. За прошлые десять лет большой опыт использования солнечного паруса был получен космическими агентствами США, Европы и Японии [1116].

В работе [7] Я.Г. Сапункова решена задача об оптимальной встрече двух КА, один из которых, неуправляемый, движется только под действием силы притяжения к Солнцу, а другой управляется с помощью реактивного ускорения, сообщаемого центру масс КА двигателем ограниченной или импульсной тяги и солнечным парусом. В нашей статье решается задача об оптимальном выводе КА на заданную орбиту. Движение КА управляется теми же двигателями, что и в работе [7]. Обе задачи оптимального управления решаются с использованием переменных Кустаанхеймо–Штифеля и с одним и тем же функционалом качества процесса управления.

Системы дифференциальных уравнений для фазовых и сопряженных переменных краевых задач об оптимальной встрече двух КА и об оптимальном выводе КА на заданную орбиту, первые интегралы этих систем уравнений, установленные в [1, 4], которые используются как для переноса части краевых условий с правого конца на левый конец траектории, так и для контроля точности вычислений при численном решении краевых задач, для указанных выше задач оптимального управления одни и те же. Эти уравнения и первые интегралы приводятся в нашей статье. Постановки задач, граничные условия, условия трансверсальности и, естественно, численные решения задач и закономерности оптимального управления движением центра масс КА в работе [7] и в нашей работе различны.

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

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

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

1. ПОСТАНОВКА ЗАДАЧИ

KS-переменные u = (u0, u1, u2, u3), s = (s0, s1, s2, s3) связаны с векторами r и v положения центра масс КА и его скорости соотношениями [13]

${\mathbf{r}} = {{P}^{T}}({\mathbf{u}}){\mathbf{u}};\,\,\,\,{\mathbf{v}} = \frac{2}{r}{{P}^{T}}({\mathbf{u}})s;\,\,\,\,r = \left| {\mathbf{r}} \right| = {{u}^{2}},$
где PT(u) – матрица, транспонированная к матрице

$P({\mathbf{u}}) = \left\| {\begin{array}{*{20}{c}} {{{u}_{0}}}&{ - {{u}_{3}}}&{{{u}_{2}}} \\ {{{u}_{1}}}&{{{u}_{2}}}&{{{u}_{3}}} \\ { - {{u}_{2}}}&{{{u}_{1}}}&{{{u}_{0}}} \\ { - {{u}_{3}}}&{ - {{u}_{0}}}&{{{u}_{1}}} \end{array}} \right\|.$

Отметим, что матрица P(u) была введена в статье [3] Я.Г. Сапункова. Матрица PТ(u) (матрица, транспонированная к матрице P(u)) отличается от квадратной KS-матрицы L(u), которая используется в книге [1], тем, что в матрице L(u) была опущена последняя строка. В результате получается матрица PТ(u), в которой 3 строки и 4 столбца. Матрица P(u), у которой 4 строки и 3 столбца, действуя на трехмерный вектор, производит четырехмерный вектор. Введение матрицы P(u) оказалось удобным приемом при решении задач оптимального управления движением КА в KS-переменных. Кроме того, в отличие от [1] в статье [3] и в нашей статье у компонент KS-переменных индекс “4” заменен на индекс “0” по аналогии с кватернионами и кватернионными матрицами, которые использовались в наших работах [2, 8].

Переменная h определяет полную энергию единицы массы КА

(1.1)
$h = \frac{{{{v}^{2}}}}{2} - \frac{{\gamma M}}{r},$
где М – масса притягивающего центра, γ – гравитационная постоянная, v = |v|.

Переменная τ, функциями которой являются переменные ${\mathbf{u}}\,,\,\,{\mathbf{s}}\,,\,\,h$, связана с временем t дифференциальным соотношением

(1.2)
${{dt} \mathord{\left/ {\vphantom {{dt} {d\tau }}} \right. \kern-0em} {d\tau }} = {{u}^{2}}.$

Движение центра масс КА в центральном гравитационном ньютоновском поле под действием ускорения ${\mathbf{\tilde {p}}}$ центра масс КА произвольной природы, в KS-переменных описывается системой дифференциальных уравнений [1]:

(1.3)
$\begin{gathered} \frac{{d{\mathbf{u}}}}{{d\tau }} = {\mathbf{s}},\,\,\,\,\frac{{d{\mathbf{s}}}}{{d\tau }} = \frac{1}{2}h{\mathbf{u}} + \frac{1}{2}{{{\mathbf{u}}}^{2}}P\left( {\mathbf{u}} \right){\mathbf{\tilde {p}}}, \\ \frac{{dh}}{{d\tau }} = 2\left( {{\mathbf{s}},P\left( {\mathbf{u}} \right){\mathbf{\tilde {p}}}} \right). \\ \end{gathered} $

Ускорение ${{{\mathbf{p}}}_{{{\text{sol}}}}}$ центра масс КА под действием солнечного паруса определяется по формуле [4]

(1.4)
${{{\mathbf{p}}}_{{{\text{sol}}}}} = d\frac{{{\text{co}}{{{\text{s}}}^{2}}\theta }}{{{{r}^{2}}}}{\mathbf{n}} = d{{\left( {{{{\mathbf{u}}}^{2}}} \right)}^{{ - 4}}}{{\left( {{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{u}},{\mathbf{n}}} \right)}^{2}}{\mathbf{n}},$
в которой n единичный вектор нормали к плоскости паруса, обращенной от Солнца, θ – угол между векторами r и n, d – коэффициент, характеризующий площадь паруса, отнесенную к единице массы аппарата.

Вектор ускорения центра масс КА от тяги двигателя обозначается через p и на него наложено ограничение

(1.5)
$\left| {\mathbf{p}} \right| \leqslant {{p}_{{\max }}}.$

В начальный момент времени КА находится на эллиптической орбите, которая определяется классическими элементами (а – большая полуось, е – эксцентриситет, I – угол наклона орбиты, ${{\Omega }_{u}}$ – долгота восходящего узла, ${{\omega }_{\pi }}$ – угловое расстояние от восходящего узла до перицентра):

(1.6)
$a = {{a}_{n}},\,\,\,\,e = {{e}_{n}},\,\,\,\,I = {{I}_{n}},\,\,\,\,{{\Omega }_{u}} = {{\Omega }_{{un}}},\,\,\,\,{{\omega }_{\pi }} = {{\omega }_{{\pi n}}}.$

Истинная аномалия КА на орбите φ в начальный момент определяется полярным углом $\varphi = {{\varphi }_{n}}$.

Требуется перевести КА на новую эллиптическую орбиту с классическими элементами

(1.7)
$a = {{a}_{k}},\,\,\,\,e = {{e}_{k}},\,\,\,\,I = {{I}_{k}},\,\,\,\,{{\Omega }_{u}} = {{\Omega }_{{uk}}},\,\,\,\,{{\omega }_{\pi }} = {{\omega }_{{\pi k}}}.$

На основе работы [1] по классическим элементам начальной орбиты (1.6) вычисляются кватернионные элементы начальной орбиты ${{{\mathbf{C}}}_{n}},{{{\mathbf{D}}}_{n}}$ со следующими компонентами:

(1.8)
$\begin{gathered} {{C}_{{0n}}} = \sqrt {{{a}_{n}}(1 - {{e}_{n}})} \sin \frac{{{{I}_{n}}}}{2}\cos \frac{{{{\Omega }_{{un}}} - {{\omega }_{{\pi n}}}}}{2}, \\ {{C}_{{1n}}} = \sqrt {{{a}_{n}}(1 - {{e}_{n}})} \sin \frac{{{{I}_{n}}}}{2}\sin \frac{{{{\Omega }_{{un}}} - {{\omega }_{{\pi n}}}}}{2}, \\ {{C}_{{2n}}} = \sqrt {{{a}_{n}}(1 - {{e}_{n}})} \cos \frac{{{{I}_{n}}}}{2}\sin \frac{{{{\Omega }_{{un}}} + {{\omega }_{{\pi n}}}}}{2}, \\ {{C}_{{3n}}} = - \sqrt {{{a}_{n}}(1 - {{e}_{n}})} \cos \frac{{{{I}_{n}}}}{2}\cos \frac{{{{\Omega }_{{un}}} + {{\omega }_{{\pi n}}}}}{2}, \\ {{D}_{{0n}}} = \sqrt {{{a}_{n}}(1 + {{e}_{n}})} \sin \frac{{{{I}_{n}}}}{2}\sin \frac{{{{\Omega }_{{un}}} - {{\omega }_{{\pi n}}}}}{2}, \\ {{D}_{{1n}}} = - \sqrt {{{a}_{n}}(1 + {{e}_{n}})} \sin \frac{{{{I}_{n}}}}{2}\cos \frac{{{{\Omega }_{{un}}} - {{\omega }_{{\pi n}}}}}{2}, \\ {{D}_{{2n}}} = \sqrt {{{a}_{n}}(1 + {{e}_{n}})} \cos \frac{{{{I}_{n}}}}{2}\cos \frac{{{{\Omega }_{{un}}} + {{\omega }_{{\pi n}}}}}{2}, \\ {{D}_{{3n}}} = \sqrt {{{a}_{n}}(1 + {{e}_{n}})} \cos \frac{{{{I}_{n}}}}{2}\sin \frac{{{{\Omega }_{{un}}} + {{\omega }_{{\pi n}}}}}{2}. \\ \end{gathered} $

Кватернионные элементы конечной орбиты ${{{\mathbf{C}}}_{k}},{{{\mathbf{D}}}_{k}}$ вычисляются по формулам аналогичным (1.8), но с использованием соотношений (1.7).

Начальное состояние управляемого КА в KS-переменных определяется по формулам

(1.9)
$\begin{gathered} {{{\mathbf{u}}}_{n}} = {{{\mathbf{C}}}_{n}}{\text{cos}}\frac{{{{E}_{n}}}}{2} + {{{\mathbf{D}}}_{n}}{\text{sin}}\frac{{{{E}_{n}}}}{2}, \\ {{{\mathbf{s}}}_{n}} = k\left( {{{{\mathbf{D}}}_{n}}{\text{cos}}\frac{{{{E}_{n}}}}{2} - {{{\mathbf{C}}}_{n}}{\text{sin}}\frac{{{{E}_{n}}}}{2}} \right), \\ {{h}_{n}} = - Q_{n}^{{ - 1}},\,\,\,\,k = {{(2{{Q}_{n}})}^{{ - 1/2}}},\,\,\,\,{{Q}_{n}} = {\mathbf{C}}_{n}^{2} + {\mathbf{D}}_{n}^{2}, \\ \end{gathered} $
где ${{E}_{n}}$ – эксцентрическая аномалия КА в начальный момент времени, связанная с истинной аномалией ${{\varphi }_{n}}$ соотношением ${\text{tg}}\frac{{{{\varphi }_{n}}}}{2} = \sqrt {\frac{{1 + {{e}_{n}}}}{{1 - {{e}_{n}}}}} {\kern 1pt} {\text{tg}}\frac{{{{E}_{n}}}}{2}.$

Задачу об оптимальном выводе КА на заданную орбиту в KS-переменных удобно решать как задачу о мягкой встрече управляемого КА с неуправляемым КА, который движется по заданной конечной орбите. Но в начальный момент движения управляемого КА эксцентрическая аномалия неуправляемого КА на его орбите не задается и подлежит определению при решении задачи оптимального управления. Закон движения неуправляемого КА в KS-переменных по заданной конечной орбите определяется по формулам

$\begin{gathered} {{{\mathbf{u}}}_{a}} = {{{\mathbf{C}}}_{k}}{\text{cos}}({{k}_{a}}{{\tau }_{a}}) + {{{\mathbf{D}}}_{k}}{\text{sin}}({{k}_{a}}{{\tau }_{a}}), \\ {{{\mathbf{s}}}_{a}} = {{k}_{a}}({{{\mathbf{D}}}_{k}}{\text{cos}}({{k}_{a}}{{\tau }_{a}}) - {{{\mathbf{C}}}_{k}}{\text{sin}}({{k}_{a}}{{\tau }_{a}})), \\ {{h}_{a}} = - Q_{k}^{{ - 1}},\,\,\,\,{{k}_{a}} = {{(2{{Q}_{k}})}^{{ - 1/2}}}, \\ {{Q}_{k}} = {\mathbf{C}}_{k}^{2} + {\mathbf{D}}_{k}^{2},\,\,\,\,{{\tau }_{a}} \geqslant {{\tau }_{{an}}}. \\ \end{gathered} $

Переменная ${{\tau }_{a}}$, функциями которой являются KS-переменные неуправляемого КА ${{{\mathbf{u}}}_{a}},{{{\mathbf{s}}}_{a}}$, связана со временем t соотношением:

(1.10)
${{dt} \mathord{\left/ {\vphantom {{dt} {d{{\tau }_{a}}}}} \right. \kern-0em} {d{{\tau }_{a}}}} = u_{a}^{2}.$

Начальное значение ${{\tau }_{a}} = {{\tau }_{{an}}}$ заранее не задается и подлежит определению.

В качестве масштаба длины вводится параметр начальной орбиты КА $p = {{a}_{n}}{\text{(}}1 - e_{n}^{2}{\text{)}}$, тогда связь между размерными и безразмерными переменными, которые обозначаются верхним индексом “*”, будет определяться соотношениями:

(1.11)
$\begin{gathered} {\mathbf{u}} = {{p}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}{\mathbf{u}}*;\,\,\,\,{\mathbf{s}} = {{\left( {\gamma M} \right)}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}{\mathbf{s}}*;\,\,\,\,h = \frac{{\gamma M}}{p}h*; \\ \tau = {{\left( {\frac{p}{{\gamma M}}} \right)}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}\tau *;\,\,\,\,t = p{{\left( {\frac{p}{{\gamma M}}} \right)}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}t*;\,\,\,\,d = \gamma Md*, \\ {\mathbf{r}} = p{\mathbf{r}}*,\,\,\,\,{\mathbf{v}} = {{\left( {\frac{{\gamma M}}{p}} \right)}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}{\mathbf{v}}*,\,\,\,\,{\mathbf{p}} = \frac{{\gamma M}}{{{{p}^{2}}}}{\mathbf{p}}*, \\ {{p}_{{\max }}} = \frac{{\gamma M}}{{{{p}^{2}}}}p_{{\max }}^{*}. \\ \end{gathered} $

Безразмерные величины $d{\text{*}}$ и $p_{{\max }}^{*}$ определяют отношения максимальных величин ускорения от тяги солнечного паруса и ускорения от тяги двигателя к ускорению от силы притяжения аппарата к центру на орбите радиуса p. Далее будут использоваться только безразмерные величины и верхний индекс “*” над ними опускается.

Уравнения движения управляемого КА (1.3), в которых единичный вектор нормали n и вектор ускорения p от тяги двигателя являются управляющими параметрами, с учетом уравнения для τa и уравнений (1.2), (1.4) и (1.10) в безразмерных KS‑переменных можно записать в виде:

(1.12)
$\begin{gathered} \frac{{d{\mathbf{u}}}}{{d\tau }} = {\mathbf{s}},\,\,\,\,\frac{{d{\mathbf{s}}}}{{d\tau }} = \frac{1}{2}h{\mathbf{u}} + \frac{1}{2}d{{({{{\mathbf{u}}}^{2}})}^{{ - 3}}} \times \\ \times \,\,{{({{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{u}},{\mathbf{n}})}^{2}}P\left( {\mathbf{u}} \right){\mathbf{n}} + \frac{1}{2}{{{\mathbf{u}}}^{2}}P\left( {\mathbf{u}} \right){\mathbf{p}}, \\ \frac{{dh}}{{d\tau }} = 2d\left( {{\mathbf{s}},P\left( {\mathbf{u}} \right){\mathbf{n}}} \right){{({{{\mathbf{u}}}^{2}})}^{{ - 4}}}{{({{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{u}},{\mathbf{n}})}^{2}} + \\ + \,\,2\left( {{\mathbf{s}},P\left( {\mathbf{u}} \right){\mathbf{p}}} \right),\,\,\,\,\,\frac{{d{{\tau }_{a}}}}{{d\tau }} = \frac{{{{{\mathbf{u}}}^{2}}}}{{{\mathbf{u}}_{a}^{2}\left( {{{\tau }_{a}}} \right)}}. \\ \end{gathered} $

Последнее уравнение системы (1.12) связывает между собой переменные $\tau $ и ${{\tau }_{a}}$, которые являются переменными фиктивного времени для управляемого и неуправляемого аппаратов соответственно. Ограничение (1.5) на вектор p ускорения центра масс КА от тяги двигателя при переходе к безразмерным переменным сохраняет свой вид.

В безразмерных KS-переменных качество процесса управления определяется функционалом с весовыми множителями ${{\alpha }_{0}} \geqslant 0,{{\alpha }_{1}} \geqslant 0$

(1.13)
$I = \int\limits_0^{{{t}_{k}}} {\left( {{{\alpha }_{0}} + {{\alpha }_{1}}\left| {\mathbf{p}} \right|} \right)} dt = \int\limits_0^{{{\tau }_{k}}} {\left( {{{\alpha }_{0}} + {{\alpha }_{1}}\left| {\mathbf{p}} \right|} \right)} {{{\mathbf{u}}}^{2}}d\tau .$

Функционал (1.13) представляет собой свертку с весовыми множителями двух критериев: времени и суммарного импульса величины ускорения КА от тяги двигателя, затраченных на процесс управления.

Требуется найти оптимальное управление КА с солнечным парусом и реактивным ускорением, удовлетворяющим ограничению (1.5), которое обеспечивает перевод КА из начального состояния до мягкой встречи с неуправляемым аппаратом, т.е. обеспечивает вывод КА на заданную орбиту и сообщает минимальное значение функционалу (1.13).

Для системы (1.12) известны первые интегралы [1]

(1.14)
$\begin{gathered} ~l\left( {{\mathbf{u}},{\mathbf{s}}} \right) = {{u}_{0}}{{s}_{1}} - {{u}_{1}}{{s}_{0}} + {{u}_{2}}{{s}_{3}}--{{u}_{3}}{{s}_{2}} = 0; \\ h{{u}^{2}} - 2{{s}^{2}} + 1 = 0. \\ \end{gathered} $

Начальное состояние управляемой системы (1.12) в пространстве (u, s, h, τa) при τ = 0 согласно (1.1) и (1.9) определяются соотношениями

(1.15)
${\mathbf{u}} = {{{\mathbf{u}}}_{n}},\,\,\,\,{\mathbf{s}} = {{{\mathbf{s}}}_{n}},\,\,\,\,h = {{h}_{n}}.$

В конечный “момент времени” τ = τk, который заранее не задается, управляемая система (1.12) в пространстве (u, s, h, τa) в момент выхода на заданную орбиту, т.е., когда у управляемого и неуправляемого аппаратов совпадают их положения и скорости, должна находиться на многообразии

(1.16)
$\begin{gathered} {{P}^{T}}({\mathbf{u}}({{\tau }_{k}})){\mathbf{u}}({{\tau }_{k}}) = {{P}^{T}}({{{\mathbf{u}}}_{a}}({{\tau }_{a}}({{\tau }_{k}}))){{{\mathbf{u}}}_{a}}({{\tau }_{a}}({{\tau }_{k}})), \\ {{P}^{T}}({\mathbf{u}}({{\tau }_{k}})){\mathbf{s}}({{\tau }_{k}}) = {{P}^{T}}({{{\mathbf{u}}}_{a}}({{\tau }_{a}}({{\tau }_{k}}))){{{\mathbf{s}}}_{a}}({{\tau }_{a}}({{\tau }_{k}})). \\ \end{gathered} $

Таким образом: требуется найти оптимальное управление солнечным парусом и реактивным ускорением от тяги двигателя с учетом ограничения (1.5), которое переводит управляемую систему (1.12) из начального состояния (1.15) на многообразие (1.16) и сообщает минимальное значение функционалу (1.13).

2. РЕШЕНИЕ ЗАДАЧИ С ПОМОЩЬЮ ПРИНЦИПА МАКСИМУМА ПОНТРЯГИНА

Функция Гамильтона–Понтрягина для управляемой системы (1.12) с учетом функционала (1.13) имеет вид

$\begin{gathered} H = - \left( {{{\alpha }_{0}} + {{\alpha }_{1}}\left| {\mathbf{p}} \right|} \right){{{\mathbf{u}}}^{2}} + \left( {{\mathbf{\mu }},{\mathbf{s}}} \right) + \frac{1}{2}h\left( {{\mathbf{\nu }},{\mathbf{u}}} \right) + \\ + \,\,d{{u}^{{ - 8}}}{{({{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{u}},{\mathbf{n}})}^{2}}\left( {{\mathbf{q}},P\left( {\mathbf{u}} \right){\mathbf{n}}} \right) + \\ + \,\,\left( {{\mathbf{q}},P\left( {\mathbf{u}} \right){\mathbf{p}}} \right) + \vartheta \frac{{{{{\mathbf{u}}}^{2}}}}{{{{{\left( {{{{\mathbf{u}}}_{a}}\left( {{{\tau }_{a}}} \right)} \right)}}^{2}}}},\,\,\,\,{\mathbf{q}} = \frac{{{{u}^{2}}}}{2}{\mathbf{\nu }} + 2\eta {\mathbf{s}}{\kern 1pt} . \\ \end{gathered} $

Сопряженные переменные ${\mathbf{\mu }} = {\text{(}}{{\mu }_{0}},{{\mu }_{1}},{{\mu }_{2}},{{\mu }_{3}}{\text{)}},$ ${\mathbf{\nu }} = {\text{(}}{{\nu }_{0}},{{\nu }_{1}},{{\nu }_{2}},{{\nu }_{3}}{\text{)}},\eta ,\vartheta $, соответствующие фазовым переменным u, s, h, τa, удовлетворяют системе уравнений

    (2.1)
$\begin{gathered} \frac{{d{\mathbf{\mu }}}}{{d\tau }} = 2\left( {{{\alpha }_{0}} + {{\alpha }_{1}}\left| {\mathbf{p}} \right|} \right){\mathbf{u}} - \frac{h}{2}{\mathbf{\nu }} - \\ - \,\,d\left\{ { - \frac{8}{{{{u}^{2}}}}} \right.\left( {{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{u}},{\mathbf{n}}} \right)\left( {{\mathbf{q}},P\left( {\mathbf{u}} \right){\mathbf{n}}} \right){\mathbf{u}} + \\ + \,\,4\left( {{\mathbf{q}},P\left( {\mathbf{u}} \right){\mathbf{n}}} \right)P\left( {\mathbf{u}} \right){\mathbf{n}} + \left( {{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{u}},{\mathbf{n}}} \right) \times \\ \left. {\frac{{}}{{}} \times \,\,\left[ {\left( {{\mathbf{\nu }},P\left( {\mathbf{u}} \right){\mathbf{n}}} \right){\mathbf{u}} + P\left( {\mathbf{q}} \right){\mathbf{n}}} \right]} \right\}\left( {{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{u}},{\mathbf{n}}} \right){{u}^{{ - 8}}} - \\ - \,\,\left( {{\mathbf{\nu }},P\left( {\mathbf{u}} \right){\mathbf{p}}} \right){\mathbf{u}} - P\left( {\mathbf{q}} \right){\mathbf{p}} - 2\vartheta \frac{{\mathbf{u}}}{{{{{\left( {{{{\mathbf{u}}}_{a}}\left( {{{\tau }_{a}}} \right)} \right)}}^{2}}}}, \\ \end{gathered} $
$\begin{gathered} \frac{{d{\mathbf{\nu }}}}{{d\tau }} = \, - {\mathbf{\mu }} - 2\eta d{{\left( {{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{u}},{\mathbf{n}}} \right)}^{2}}{{u}^{{ - 8}}}P\left( {\mathbf{u}} \right){\mathbf{n}} - 2\eta P\left( {\mathbf{u}} \right){\mathbf{p}}, \\ \frac{{d\eta }}{{d\tau }} = - \frac{1}{2}\left( {{\mathbf{\nu }},{\mathbf{u}}} \right),\,\,\,\,\frac{{d\vartheta }}{{d\tau }} = 2\vartheta {{u}^{2}}\frac{{\left( {{{{\mathbf{u}}}_{a}},{{{\mathbf{s}}}_{a}}} \right)}}{{u_{a}^{4}}}. \\ \end{gathered} $

Из последнего уравнения системы (2.1) следует, что $\vartheta = Q{{\left( {{{{\mathbf{u}}}_{a}}({{\tau }_{a}})} \right)}^{2}}$, где Q – произвольная константа.

Согласно условию максимума для функции Гамильтона–Понтрягина оптимальное управление для единичного вектора нормали n к плоскости солнечного паруса, определяется по формуле

(2.2)
$\begin{gathered} {\mathbf{n}} = z\frac{{\mathbf{r}}}{r} + {{b}^{{ - 1}}}\left( {z - \frac{2}{{3z}}} \right)\frac{{{{{\mathbf{R}}}_{1}}}}{r},\,\,\,\,{{{\mathbf{R}}}_{1}} = {{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{q}} - b{\mathbf{r}}, \\ b = \frac{{\left( {{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{q}},{\mathbf{r}}} \right)}}{{{{r}^{2}}}},\,\,\,\,z = {{\left[ {\frac{1}{6}\left( {4 - {{a}^{2}} + a\sqrt {8 + {{a}^{2}}} } \right)} \right]}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}, \\ a = b\frac{r}{{\left| {{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{q}}} \right|}}, \\ \end{gathered} $
а для ограниченной тяги оптимальное управление определяется по формуле

(2.3)
$\begin{gathered} {{{\mathbf{p}}}_{{{\text{opt}}}}} = 0,\,\,\,\,{\text{если}}\,\,\,\,\left| {{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{q}}} \right| \leqslant {{\alpha }_{1}}{{u}^{2}}, \\ {{{\mathbf{p}}}_{{{\text{opt}}}}} = \frac{{{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{q}}}}{{\left| {{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{q}}} \right|}}{{p}_{{{\text{max}}}}},\,\,\,\,{\text{если}}\,\,\,\,\left| {{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{q}}} \right| > {{\alpha }_{1}}{{u}^{2}}. \\ \end{gathered} $

Из соотношений (2.3) видно, что управление реактивным ускорением от двигателя малой тяги может состоять из нескольких этапов, на которых двигатель выключен или вырабатывает максимальную тягу.

На левом подвижном конце траектории в пространстве (u, s, h, τa), где значение ${{\tau }_{a}}$ не задается заранее должно выполняться условие трансверсальности

(2.4)
$\vartheta (0) = 0.$

Из последнего уравнения системы (2.1) и условия трансверсальности (2.4) следует, что сопряженная переменная $\vartheta $ во все время движения равна нулю: $\vartheta \equiv 0$.

На правом подвижном конце траектории при τ = τk должны выполняться условия трансверсальности:

(2.5)
$\begin{gathered} l\left( {{\mathbf{\mu }},{\mathbf{u}}} \right) + l\left( {{\mathbf{\nu }},{\mathbf{s}}} \right) = 0;\,\,\,\,l\left( {{\mathbf{\nu }},{\mathbf{u}}} \right) = 0;\,\,\,\,\eta = {\text{ }}0; \\ \left( {{\mathbf{s}},{\mathbf{\mu }}} \right) + h\left( {{\mathbf{u}},{\mathbf{\nu }}} \right) = 0{\text{ }}. \\ \end{gathered} $

Соотношения (2.5) получены из (1.16) методом неопределенных множителей Лагранжа с последующим исключением множителей Лагранжа.

Кроме того, так как τk заранее не задается, то для оптимального процесса при $\tau = {{\tau }_{k}}$ должно выполняться условие для функции Гамильтона–Понтрягина

(2.6)
${{\left. {{{H}_{{{\text{opt}}}}}} \right|}_{{{{\tau }_{k}}}}} = 0.$

Принцип максимума сводит решение задачи оптимального вывода КА на заданную орбиту к решению краевой задачи для системы дифференциальных уравнений по определению фазовых и сопряженных переменных. Краевая задача оптимального управления формулируется для системы дифференциальных уравнений (1.12), (2.1), в которых управляющие параметры n и р согласно условию максимума для функции Гамильтона–Понтрягина определяется по формулам (2.2) и (2.3) соответственно, с начальными условиями (1.15) и (2.4) при τ = 0 и граничными условиями (1.16), (2.5), (2.6) при τ = τk.

3. ПЕРВЫЕ ИНТЕГРАЛЫ

Кроме билинейной формы l(a, b), которая определена первым из соотношений (1.14), введем согласно [4] билинейные формы следующих видов:

l1(a, b) = a0b2a2b0 + a1b3a3b1; l2(a, b) = a0b3a3b0 + a2b1a1b2;

l3(a, b) = a0b1a1b0a2b3 + a3b2.

Система уравнений (1.12), (2.1), в которых управляющие параметры определяется согласно соотношениям (2.2) и (2.3), кроме первых интегралов (1.14), имеет также следующие интегралы [4]:

l(μ, u) + l(ν, s) = 0; li (μ, u) + li (ν, s) = ci = const,

i = 1, 2, 3;

(3.1)
$\frac{\vartheta }{{{{{\left( {{{{\mathbf{u}}}_{a}}({{\tau }_{a}})} \right)}}^{2}}}} = Q = {\text{const}};\,\,\,{{H}_{{{\text{opt}}}}}(\tau ) \equiv 0.$

Из (2.5) и (3.1) следует, что первое условие трансверсальности из условий (2.5) можно заменить условием l(μ, u) + l(ν, s) = 0 при τ = 0, т.е. перенести это условие с правого конца траектории на левый конец.

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

При неограниченном возрастании pmax решение рассмотренной выше задачи будет стремиться к решению задачи об оптимальном управлении движением КА с импульсным реактивным ускорением [4]. Если через Ki обозначить безразмерную величину импульса реактивного ускорения от тяги двигателя, то она будет связана с размерной величиной изменения вектора скорости $\left| {\Delta {{{\mathbf{v}}}_{i}}} \right|$ под действием импульса соотношением $\left| {\Delta {{{\mathbf{v}}}_{i}}} \right| = {{(\gamma M{{R}^{{ - 1}}})}^{{1/2}}}{{K}_{i}}$.

Для импульсной задачи в безразмерных переменных аналогом функционала (1.13) будет функционал $I = {{\alpha }_{0}}{{t}_{k}} + {{\alpha }_{1}}\sum\nolimits_{i = 1}^n {{{K}_{i}}} $, где tk – время процесса управления, n – число импульсов. Величины tk, n заранее не задаются.

Момент фиктивного времени сообщения импульса ускорения, согласно (2.3), определяется из условия $\left| {{{P}^{T}}\left( {\mathbf{u}} \right){\mathbf{q}}} \right| = {{\alpha }_{1}}{{u}^{2}}$, которое можно привести к виду

(4.1)
${{q}^{2}} - \frac{1}{4}{{u}^{2}}{{\left( {l{\text{(}}{\mathbf{u}},{\mathbf{\nu }}{\text{)}}} \right)}^{2}} - \alpha _{1}^{2}{{u}^{2}} = 0.$

Направление вектора ускорения от тяги двигателя в момент сообщения импульса совпадает с направлением вектора PT(u)q.

Из уравнений систем (1.12), (2.1) следует, что при предельном переходе к импульсному ускорению в момент сообщения импульса переменные u, τa, η, ϑ, q, l(ν, u) остаются непрерывными, а переменные s, h, μ, ν, (s, q), (ν, q) испытывают разрыв первого рода или скачок. Обозначим через Δs, Δh, Δμ, Δν изменение переменных s, h, μ, ν при сообщении импульса Ki. Вычисляя интегралы от правых частей уравнений систем (1.12) и (2.1) на промежутке действия реактивного ускорения и выполняя предельный переход к импульсному ускорению, получим

(4.2)
$\begin{gathered} \Delta {\mathbf{s}} = \frac{1}{{2{{\alpha }_{1}}}}\left[ {{\mathbf{q}} - \frac{1}{2}l{\text{(}}{\mathbf{\nu }},{\mathbf{u}}{\text{)}}{{{\mathbf{u}}}^{d}}} \right]{{K}_{i}}; \\ \Delta h = 2\frac{{{\text{(}}{\mathbf{s}},{\mathbf{q}}{\text{)}}}}{{{{\alpha }_{1}}{{u}^{2}}}}{{K}_{i}} + \frac{1}{2}{{K}_{i}}^{2}; \\ \Delta {\mathbf{\mu }}\, = \,\frac{1}{{{{\alpha }_{1}}{{u}^{2}}}}\left\{ {\left[ {\frac{{{{q}^{2}}}}{{{{u}^{2}}}} - {\text{(}}{\mathbf{\nu }},{\mathbf{q}}{\text{)}} + {{\alpha }_{1}}\eta {{K}_{i}}} \right]{\mathbf{u}} - \frac{1}{2}l{\text{(}}{\mathbf{\nu }},{\mathbf{u}}){{{\mathbf{q}}}^{d}}} \right\}{{K}_{i}}; \\ \Delta {\mathbf{\nu }} = - 2\frac{\eta }{{{{\alpha }_{1}}{{u}^{2}}}}\left[ {{\mathbf{q}} - \frac{1}{2}l{\text{(}}{\mathbf{\nu }},{\mathbf{u}}{\text{)}}{{{\mathbf{u}}}^{d}}} \right]{{K}_{i}}. \\ \end{gathered} $

Векторы ud и qd ортогональны соответственно к векторам u и q и выражаются через них по формулам

${{{\mathbf{u}}}^{d}} = \left( {{{u}_{1}}, - {{u}_{0}},{{u}_{3}}, - {{u}_{2}}} \right){\text{;}}\,\,\,\,{{{\mathbf{q}}}^{d}} = \left( {{{q}_{1}}, - {{q}_{0}},{{q}_{3}}, - {{q}_{2}}} \right).$

Значения всех величин в правых частях соотношений (4.2) вычисляются непосредственно перед сообщением импульса. Из условия, что для оптимального процесса Hopt ≡ 0, следует, что во все время движения, исключая моменты сообщения импульса ускорения, выполняется соотношение

(4.3)
$( - {{\alpha }_{1}} + Q{\text{)}}{{{\mathbf{u}}}^{2}} + {\text{(}}{\mathbf{\mu }},{\mathbf{s}}{\text{)}} + \frac{1}{2}h{\text{(}}{\mathbf{\nu }},{\mathbf{u}}{\text{)}} \equiv 0.$

Если момент сообщения импульса не совпадает с начальным или конечным моментом времени, то перед моментом сообщения импульса, кроме соотношения (4.1), должно выполняться соотношение

(4.4)
$\begin{gathered} \frac{{{{q}^{2}}}}{{{{u}^{2}}}}{\text{(}}{\mathbf{u}},{\mathbf{s}}{\text{)}} + \frac{{{{u}^{2}}}}{2}l{\text{(}}{\mathbf{\nu }},{\mathbf{u}}{\text{)}}l{\text{(}}{\mathbf{\nu }},{\mathbf{s}}{\text{)}} + \\ + \,\,\left( {{\mathbf{q}},\left[ {\frac{1}{2}{{u}^{2}}{\mathbf{\mu }} + {\text{(}}{\mathbf{\nu }},{\mathbf{u}}{\text{)}}{\mathbf{s}} - {\text{(}}{\mathbf{u}},{\mathbf{s}}{\text{)}}{\mathbf{\nu }} - \eta h{\mathbf{u}}} \right]} \right) = 0. \\ \end{gathered} $

Между двумя последовательными моментами сообщения импульса, т.е. на участках движения КА только под действием солнечного паруса, фазовые и сопряженные переменные определяются интегрированием систем уравнений (1.12) и (2.1), в которых вектор ускорения от тяги двигателя ${\mathbf{p}} = 0$.

Первые интегралы, приведенные выше для случая ограниченного ускорения, имеют место и для импульсного ускорения.

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

Рассмотрим случай, когда в задаче о выводе КА на орбиту импульсы сообщаются только в начальный и конечный моменты времени. В этой задаче неизвестными являются значения сопряженных переменных μ, ν, η при τ = 0, величины ${{\tau }_{{an}}}$, τk и значения импульсов ${{K}_{n}}$ и Kк. То есть, имеются две векторные и пять скалярных неизвестных величин. В скалярной форме имеем всего 13 неизвестных. Для определения этих величин имеются следующие условия: граничные условия (1.16), условия трансверсальности (2.5), условие (4.1) при τ = 0 и τ = τk, условие (4.3), которое можно выполнить при любом τ ∈ (0,τk). Всего в скалярной форме получается 13 условий.

Если решение задачи оптимального управления содержит импульсы ускорений во внутренних точках, то для определения “момента” сообщения импульса τi и значения величины предыдущего импульса служат условия (4.1) и (4.4).

5. ПРИМЕРЫ ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧИ

Для численного решения краевой задачи, к которой принцип максимума Понтрягина сводит решение поставленной задачи оптимального управления, используется комбинация методов Ньютона и градиентного спуска. Система дифференциальных уравнений для определения фазовых и сопряженных переменных интегрируется методом Рунге–Кутта четвертого порядка точности. Использование KS-переменных приводит к улучшению сходимости итерационного процесса при решении краевой задачи оптимального управления. Это связано с тем, что при малой тяге нелинейные слагаемые в системе дифференциальных уравнений (1.12), (2.1) оказываются малыми величинами. Для улучшения сходимости итерационного процесса производилась специальная процедура уточнения моментов разрыва управляющих параметров, чтобы момент разрыва не находился внутри шага интегрирования. Для этого в окрестности момента разрыва управления шаг интегрирования уменьшался. Затем уменьшенный шаг, внутри которого оказывался разрыв управления, делился на пропорциональные отрезки моментом разрыва.

Ниже приводятся результаты численного решения восьми задач об оптимальном выводе КА на заданную орбиту.

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

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

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

(5.1)
$\begin{gathered} {{a}_{n}} = 1.0,\,\,\,\,{{e}_{n}} = 0.1,\,\,\,\,{{I}_{n}} = 5.0^\circ , \\ {{\Omega }_{{un}}} = 30.0^\circ ,\,\,\,\,{{\omega }_{{\pi n}}} = 50.0^\circ . \\ \end{gathered} $

Истинная аномалия КА на начальной орбите во всех расчетах ${{\varphi }_{n}} = 30.0^\circ $. В начальном условии (5.1) и в последующих соотношениях большая полуось представлена в безразмерном виде, т.е. большая полуось отнесена к значению большой полуоси на начальной орбите. Далее результаты расчетов (координаты положения и проекции вектора скорости и вектора импульса ускорения от тяги) представлены в безразмерных переменных согласно формулам (1.11).

Расчеты выполнялись в декартовой системе координат $O{{x}_{1}}{{x}_{2}}{{x}_{3}}$, начало которой находится в центре Солнца, а плоскость $O{{x}_{1}}{{x}_{2}}$ совпадает, например, с плоскостью орбиты Земли.

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

$\begin{gathered} {{x}_{1}} = - 0.313037,\,\,\,\,{{x}_{2}} = 0.861812,\,\,\,\,{{x}_{3}} = 0.078991, \\ {{v}_{1}} = - 1.037721,\,\,\,\,{{v}_{2}} = - 0.325439,\,\,\,\,{{v}_{3}} = 0.003261. \\ \end{gathered} $

В первых двух задачах максимальное значение реактивного ускорения ${{p}_{{{\text{max}}}}} = 0.1$, безразмерная величина d, входящая в формулу (1.4), определяющую тягу солнечного паруса, полагалась равной $d = 0.1$, в пятой и шестой задачах ${{p}_{{{\text{max}}}}} = 0.01$, $d = 0.005$.

В первых шести задачах КА необходимо перевести на конечную орбиту со следующими классическими элементами:

$\begin{gathered} {{a}_{k}} = 1.0,\,\,\,\,{{e}_{k}} = 0.1,\,\,\,\,{{I}_{k}} = 8.0^\circ , \\ {{\Omega }_{{uk}}} = 32.0^\circ ,\,\,\,\,{{\omega }_{{\pi k}}} = 46.0^\circ , \\ \end{gathered} $
при этом весовые множители в функционале (1.13) принимали значения

(5.2)
${{\alpha }_{0}} = 0.05,\,\,\,\,{{\alpha }_{1}} = 1.0.$

Результаты численного решения 8-ми задач приведены в табл. 1–5. В таблицах приняты следующие обозначения: ОРУ – ограниченное (малое) реактивное ускорение, ИРУ – импульсное реактивное ускорение, СП – солнечный парус, ИВУ – импульс величины ускорения, ВИУ – величина импульса ускорения, ПП – пассивный полет, J = … значение функционала (1.13), ∆ti – длительность i-го этапа в безразмерном времени, x1, x2, x3 и v1, v2, v3 – координаты центра масс КА и проекции его вектора скорости в конце соответствующего этапа.

Таблица 1.
Задача Этапы ti Координаты x1, x2, x3 Проекции вектора скорости
v1, v2, v3
1: ОРУ, СП
ИВУ = 0.025454
J = 0.118252
1, СП, d = 0.1 1.601397 –0.974937, –0.419786, 0.018548 0.292358, –0.882518, –0.095430
2, ОРУ, CП
p ≤ 0.1, d = 0.1
0.254545 –0.875570, –0.632390, –0.010163 0.483441, –0.781827, –0.129186
2: ОРУ
ИВУ = 0.055832
J = 0.151610
1, popt = pmax 0.071438 –0.386099, 0.865774, 0.080466 –1.006714, –0.402968, 0.020566
2, popt = 0 1.357228 –1.286831, –0.286831, 0.029121 0.170042, –0.941068, –0.077432
3, popt = pmax
pmax = 0.1
0.486883 –1.829351, –0.700602, –0.021735 0.543819, –0.734631, –0.128058
Таблица 2.
Задача Этапы ti Проекции impi
вектора ИРУ
Координаты x1, x2, x3 Проекции вектора
скорости v1, v2, v3
3: ИРУ, СП
J = 0.115453
1, ИРУ 0 0.000589, 0.000140, 0.000619 –0.313037, 0.861812,
0.078991
–1.037131, –0.325580,
0.021356
2, CП, d = 0.1 1.692836 0 –0.942932, –0.502999,
0.010275
0.362565, –0.851710,
–0.099735
3, ИРУ 0 0.004904, –0.002468,
–0.029437
–0.942932, –0.502999,
0.010275
0.367470, –0.854177,
–0.129173
4: ИРУ
J = 0.138683
1, ИРУ 0 –0.000727, –0.000399,
0.011423
–0.313037, 0.861812,
0.078991
–1.038445, –0.325838,
0.032159
2, пасс. полет 1.585268 0 –0.970723, –0.430755,
0.020955
0.301596, –0.893179,
–0.080642
3, ИРУ 0 0.000200, 0.006965,
–0.047458
–0.970723, –0.430755,
0.020955
0.301795, –0.886214,
–0.128100
Таблица 3.
Задача Этапы ti Координаты x1, x2, x3 Проекции вектора скорости v1, v2, v3
5: ОРУ, СП
p ≤ 0.01, d = 0.005
ИВУ = 0.056279
J = 0.426228
1, ОРУ, CП 0.159372 x1n = –0.313037,
x2n = 0.861812, x3n = 0.078991
v1n = –1.037721,
v2n = –0.325439, v3n = 0.003261,
2, CП 0.406840
3, ОРУ, CП 3.666200 x1k = –0.987523,
x2k = 0.159018, x3k = 0.092499
v1k = –0.265088,
v2k = –0.960093, v3k = –0.094687
4, CП 3.949749
5, ОРУ, CП 6.587174
6, CП 6.837068
7, ОРУ, CП 7.398773
Таблица 4.
Задача Этапы ti Элементы орбиты a, e, I, Ωu, ωπn в градусах
(an = 1.0, en = 0.1, In = 5.0°, Ωun = 30.0°, ωπn = 50.0°)
a e I u ωπ
6: ОРУ
p = pmax = 0.01
ИВУ = 0.072707
J = 0.489530
1, ОРУ 0.161611 0.999657 0.099791 5.0075 30.9694 48.8825
2, пасс. полет 0.497889 0.999657 0.099791 5.0075 30.9694 48.8825
3, ОРУ 3.638735 0.997838 0.099523 6.3435 34.7778 43.4996
4, пасс. полет 4.023727 0.997838 0.099523 6.3435 34.7778 43.4996
5, ОРУ 6.515021 0.999515 0.100276 7.2898 34.8000 43.0052
6, пасс. полет 6.859527 0.999515 0.100276 7.2898 34.8000 43.0052
7, ОРУ 8.336463 1.0 0.1 8.0 32.0 46.0
Таблица 5.
Задача Этапы ti Проекции impi
вектора ИРУ
Координаты x1, x2, x3 Проекции вектора
скорости v1, v2, v3
7: ИРУ, СП
d = 0.05
J = 0.415972
1, ИРУ 0 –0.100790, 0.027571,
–0.017959
–0.313037, 0.861812,
0.078991
–1.138511, –0.297869,
0.000942
2, CП 3.174236 0 –0.754568, –1.400625,
–0.103665
0.543272, –0.416694,
–0.048537
3, ИРУ 0 0.123220, 0.038680,
–0.078064
–0.754568, –1.400625,
–0.103665
0.666492, –0.378014,
–0.126601
8: ИРУ
J = 0.420310
1, ИРУ 0 –0.104026, 0.029748,
–0.021136
–0.313037, 0.861812,
0.078991
–1.141747, –0.295691,
–0.000399
2, пасс. полет 3.150247 0 –0.773453, –1.389731,
–0.100053
0.537191, –0.426636,
–0.047275
3, ИРУ 0 0.123934, 0.038821,
–0.080041
–0.773453, –1.389731,
–0.100053
0.661125, –0.387815,
–0.127316

Задача 1 (табл. 1). Движение КА управляется ограниченным (малым) реактивным ускорением и солнечным парусом. Оптимальное управление состоит из двух этапов. На первом этапе аппарат управляется только солнечным парусом, на втором этапе – ограниченным реактивным ускорением и солнечным парусом.

Задача 2 (табл. 1). Движение КА управляется только ограниченным реактивным ускорением. Оптимальное управление состоит из трех этапов. На первом этапе КА движется с реактивным ускорением, величина которого максимальна. В конце этого этапа орбита КА определялась следующими значениями классических элементов:

(5.3)
$\begin{gathered} {{a}_{{pas}}} = 1.002171,\,\,\,\,{{e}_{{pas}}} = 0.102049, \\ {{I}_{{pas}}} = 5.06244^\circ ,\,\,\,\,{{\Omega }_{{upas}}} = 34.17732^\circ , \\ {{\omega }_{{\pi pas}}} = 46.01705^\circ . \\ \end{gathered} $

На втором этапе ускорение от тяги двигателя ${\mathbf{p}} = 0$ и аппарат совершает пассивный полет по орбите, которая определяется классическими элементами (5.3). На третьем этапе КА движется с реактивным ускорением, величина которой максимальна.

Суммарная длительность активного полета составила $\Delta {{t}_{1}} + \Delta {{t}_{3}} = 0.558321$, величина импульса реактивного ускорения на первом этапе равна 0.007144, на третьем этапе – 0.048688.

Задача 3 (табл. 2). Движение КА управляется импульсным реактивным ускорением и солнечным парусом. Оптимальное управление состоит из двух импульсов ускорений, первый из которых выполняется в начальный момент времени, второй в конечный момент времени, и одного этапа движения КА под управлением солнечного паруса. Величина первого импульса 0.000866, величина второго – 0.029945.

Задача 4 (табл. 2). Движение КА управляется только импульсным реактивным ускорением центра масс КА. Оптимальное управление состоит из двух импульсов ускорения, первый из которых выполняется в начальный момент времени, второй в конечный момент времени, и одного этапа пассивного движения КА под действием притяжения к Солнцу. В результате действия первого импульса КА оказывается на орбите, которая определяется следующими классическими элементами:

(5.4)
$\begin{gathered} {{a}_{{pas}}} = 1.002396,\,\,\,\,{{e}_{{pas}}} = 0.102269, \\ {{I}_{{pas}}} = 5.13801^\circ ,\,\,\,{{\Omega }_{{upas}}}_{{}} = 36.60555^\circ , \\ {{\omega }_{{\pi pas}}} = 43.57026^\circ . \\ \end{gathered} $

Далее аппарат совершает пассивный полет длительностью 1.585268 по орбите с элементами (5.4). Под действием второго импульса КА оказывается на заданной конечной орбите. Величина первого импульса ускорения 0.011452, величина второго – 0.047967.

В следующих четырех задачах КА необходимо перевести на конечную орбиту с классическими элементами ${{a}_{k}} = 1.52,$ ${{e}_{k}} = 0.05,$ ${{I}_{k}} = 10.0^\circ ,$ ${{\Omega }_{{uk}}} = 40.0^\circ ,$ ${{\omega }_{{\pi k}}} = 60.0^\circ $, с весовыми множителями (5.2) в функционале (1.13) и новыми значениями ${{p}_{{\max }}} = 0.01\,,$ $d = 0.005$.

Задача 5 (табл. 3). КА управляется ограниченным реактивным ускорением и солнечным парусом. Управление состоит из семи этапов. На всех нечетных этапах КА управляется реактивным ускорением и солнечным парусом, на всех четных – только солнечным парусом. В двух последних столбцах табл. 3 указаны координаты и проекции вектора скорости КА в начальный и конечный моменты времени. Суммарная продолжительность действия реактивного ускорения составляет 5.617862.

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

$\begin{gathered} {{x}_{1}} = - 0.796080,\,\,\,\,{{x}_{2}} = - 0.743070,\,\,\,\,{{x}_{3}} = 0.029275, \\ {{v}_{1}} = 0.581179,\,\,\,{{v}_{2}} = - 0.701451, \\ {{v}_{3}} = - 0.126887. \\ \end{gathered} $

Суммарная продолжительность активных этапов управления составляет 7.270687. Проекции вектора ускорения на каждом активном этапе изменялись незначительно. В частности, в начале первого этапа эти проекции при t = 0 имели следующие значения:

$\begin{gathered} {{p}_{1}} = 0.000878,\,\,\,\,{{p}_{2}} = 0.000107, \\ {{p}_{3}} = 0.009961,\,\,\,\,\left| {\mathbf{p}} \right| = 0.01, \\ \end{gathered} $
а в конце седьмого этапа при $t = {{t}_{k}}$ эти проекции имели значения

$\begin{gathered} {{p}_{1}} = - 0.000419,\,\,\,{{p}_{2}} = 0.001024, \\ {{p}_{3}} = - 0.009939,\,\,\,\,\left| {\mathbf{p}} \right| = 0.01. \\ \end{gathered} $

Задача 7 (табл. 5). Движение КА управляется с помощью импульсного реактивного ускорения центра масс КА и солнечного паруса. Оптимальное управление состоит из двух импульсов ускорений, сообщаемых центру масс КА в начальный и конечный моменты времени, и одного этапа движения КА под действием солнечного паруса. Под действием первого импульса изменяется только скорость КА. Далее КА управляется только тягой от солнечного паруса. На этом этапе координаты КА принимают требуемые значения, проекции его скорости изменяются. Под действием второго импульса ускорения изменяется скорость КА, принимая требуемое значение. Величина первого импульса ускорения 0.106352, величина второго – 0.150908.

Задача 8 (табл. 5). Движение КА управляется только импульсным реактивным ускорением центра масс КА. Оптимальное управление состоит из двух импульсов, сообщаемых в начальный и конечный моменты времени, и одного этапа пассивного движения КА под действием притяжения к Солнцу. Под действием первого импульса изменяется только скорость КА, аппарат оказывается на орбите, которая определяется следующими классическими элементами:

(5.5)
$\begin{gathered} {{a}_{{pas}}} = 1.265685,\,\,\,{{e}_{{pas}}} = 0.294470, \\ {{I}_{{pas}}} = 4.9479^\circ ,\,\,\,\,\,\,{{\Omega }_{{upas}}} = 14.2956^\circ , \\ {{\omega }_{{\pi pas}}} = 71.5112^\circ . \\ \end{gathered} $

Далее аппарат совершает пассивный полет по орбите с элементами (5.5). В конце этого промежутка положение и скорость движения КА определяется следующими величинами:

$\begin{gathered} {{x}_{1}} = - 0.773453,\,\,\,\,{{x}_{2}} = - 1.389731, \\ {{x}_{3}} = - 0.100053,\,\,\,\,{{v}_{1}} = 0.537191, \\ {{v}_{2}} = - 0.426636,\,\,\,{{v}_{3}} = - 0.047275. \\ \end{gathered} $

Под действием второго импульса КА оказывается на заданной конечной орбите. Величина первого импульса 0.110241, величина второго – 0.152556.

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

В задачах 1–6 у начальной и конечной орбиты совпадают большие оси и эксцентриситеты соответственно. На малые величины различаются угловые элементы орбит. В этих задачах рассматривается коррекция угловых элементов орбиты КА. Углы отклонения начальной и конечной орбиты от плоскости $O{{x}_{1}}{{x}_{2}}$ малы. Углы отклонения оскулирующих орбит в процессе управления от плоскости $O{{x}_{1}}{{x}_{2}}$ также малы. Это подтверждается соотношениями (5.3) и (5.4), где приводятся значения классических элементов орбит пассивного полета. В численном решении задачи 2 компоненты ограниченного ускорения на первом этапе мало изменяются, в частности, при t = 0

(5.6)
$\begin{gathered} {{p}_{1}} = - 0.011913,\,\,\,\,{{p}_{2}} = - 0.002214, \\ {{p}_{3}} = 0.09963,\,\,\,\,\left| {\mathbf{p}} \right| = 0.1. \\ \end{gathered} $

На втором этапе реактивное ускорение отсутствует, КА совершает пассивный полет. На третьем этапе компоненты ограниченного реактивного ускорения также мало изменяются, в частности, при $t = {{t}_{k}} = 1.915561$

(5.7)
$\begin{gathered} {{p}_{1}} = 0.0019136,\,\,\,\,{{p}_{2}} = 0.014058, \\ {{p}_{3}} = - 0.098988,\,\,\,\,\left| {\mathbf{p}} \right| = 0.1. \\ \end{gathered} $

Из соотношений (5.6) и (5.7) видно, что направление вектора реактивного ускорения почти параллельно оси $O{{x}_{3}}$, так как третья проекция этого вектора по своей величине близка к величине вектора. Отсюда следует, что во время процесса управления направление вектора тяги близко к нормали к плоскости оскулирующей орбиты. Этим свойством обладает вектор реактивного ускорения для задач 1–6. Для задач 3 и 4 это видно из табл. 2, в которой приводятся значения проекций вектора импульса реактивного ускорения. Наличие солнечного паруса несколько ослабляет это свойство.

Из сравнения решений задач 5 и 6, в которых ${{p}_{{{\text{max}}}}} = 0.01,$ $d = 0.005$, видно, что наличие солнечного паруса приводит к уменьшению импульса величины реактивного ускорения, промежутка времени процесса управления, величины функционала качества. Из сравнения решений задачи 6 для ${{p}_{{{\text{max}}}}} = 0.01$ и задачи 2 для ${{p}_{{\max }}} = 0.1$ видно, что уменьшение максимальной величины реактивного ускорения приводит к увеличению импульса величины реактивного ускорения, затраченного на процесс управления, а также к увеличению времени процесса управления и к увеличению числа этапов управления.

ЗАКЛЮЧЕНИЕ

Задача о выводе КА на заданную орбиту с помощью ограниченного или импульсного реактивного ускорения центра масс КА и солнечного паруса сведена в случае использования комбинированного функционала качества процесса управления и KS-переменных к решению краевой задачи оптимального управления, которая решена численно с помощью комбинации модифицированного метода Ньютона и метода градиентного спуска. Приведены формулы, определяющие изменения фазовых и сопряженных переменных под действием импульса реактивного ускорения. Первые интегралы системы дифференциальных уравнений краевой задачи оптимального управления позволяют перенести часть краевых условий с правого конца траектории на левый конец, а также позволяют осуществить контроль точности численного решения задачи. Приведены результаты расчетов и дан их анализ для восьми вариантов задачи оптимального управления. Сделаны выводы о влиянии на характеристики процесса управления (на величины импульсов реактивного ускорения, длительность процесса управления) наличия солнечного паруса, ограниченного и импульсного реактивного ускорения центра масс КА, справедливые для построенных численных решений 8-ми нелинейных пространственных краевых задач оптимального управления, имеющих высокую размерность.

Исследование выполнено при частичной финансовой поддержке РФФИ в рамках научного проекта № 19-01-00205.

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

  1. Штифель Е., Шейфеле Г. Линейная и регулярная небесная механика. М.: Наука, 1975.

  2. Челноков Ю.Н., Сапунков Я.Г. Построение оптимальных управлений и траекторий космического аппарата на основе регулярных кватернионных уравнений задачи двух тел // Космич. исслед. 1996. Т. 34. № 2. С. 150–158. (Cosmic Research. P. 137–145).

  3. Сапунков Я.Г. Применение KS-переменных к задаче оптимального управления космическим аппаратом // Космич. исслед. 1996. Т. 34. № 4. С. 428–433.

  4. Сапунков Я.Г. Решение задач оптимального управления космическим аппаратом с ограниченной и импульсной тягой в KS-переменных // Мехатроника, автоматизация, управление. 2010. № 3. С. 73–78.

  5. Поляхова Е.Н. Космический полет с солнечным парусом: проблемы и перспективы. М.: Наука, 1986.

  6. Сапунков Я.Г. Оптимальное управление движением космического аппарата с комбинированной тягой // Математика. Механика. 2009. Вып. 11. С. 129–132.

  7. Сапунков Я.Г. Оптимальное управление космическим аппаратом с двигателем ограниченной или импульсной тяги и солнечным парусом // Мехатроника, автоматизация, управление. 2014. № 4. С. 55–61.

  8. Сапунков Я.Г., Челноков Ю.Н. Построение оптимальных управлений и траекторий центра масс космического аппарата снабженного солнечным парусом и двигателем малой тяги с использованием кватернионов и переменных Кустаанхеймо–Штифеля // Космич. исслед. 2014. Т. 52. № 6. С. 489–499. (Cosmic Research. P. 450–460).

  9. Жуков А.Н., Лебедев В.Н. Вариационная задача о перелете между гелиоцентрическими круговыми орбитами с помощью солнечного паруса // Космич. исслед. 1964. Т. 2. № 1. С. 46–50.

  10. Khabibullin R.M. Control program for noncoplanar heliocentric flight to Venus of non-perfectly reflecting solar sail spacecraft // Aerospace and Mechanical Engineering. 2019. V. 18. № 4. P. 117–128.

  11. Macdonald M. Advances in Solar Sailing // Materials of the Third International Symposium on Solar Sailing Glasgow. 2013.

  12. Johnson L., Whorton M., Heaton A. et al. NanoSail-D: A solar sail demonstration mission // Acta Astronautica. 2011. V. 68. P. 571–575.

  13. Mori O., Sawada H., Funase R. et al. First Solar Power Sail Demonstration by IKAROS // Transactions of the Japan Society for Aeronautical and Space Sciences, Aerospace Technology. 2010.

  14. Biddy C., Svitek T. LightSail-1 Solar Sail Design and Qualification // Materials of the 41th Aerospace Mechanisms Symposium. 2012. P. 451–463.

  15. Heiligers J., Diedrich B., Derbes B., McInnes C.R. Sunjammer: Preliminary End-to-End Mission Design // Materials of AIAA/AAS Astrodynamics Specialist Conference. 2014.

  16. McInnes C.R. Solar sailing: technology, dynamics and mission applications. Springer Science & Business Media, 2013.

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