Автоматика и телемеханика, № 8, 2019
Управление в технических системах
© 2019 г. Я.Г. САПУНКОВ, канд. физ.-мат. наук (ChelnokovYuN@gmail.com),
Ю.Н. ЧЕЛНОКОВ, д-р физ.-мат. наук (ChelnokovYuN@gmail.com)
(Институт проблем точной механики и управления РАН, Саратов)
ОПТИМАЛЬНЫЙ ПОВОРОТ ПЛОСКОСТИ ОРБИТЫ
КОСМИЧЕСКОГО АППАРАТА ПЕРЕМЕННОЙ МАССЫ
В ЦЕНТРАЛЬНОМ ГРАВИТАЦИОННОМ ПОЛЕ
ПОСРЕДСТВОМ ОРТОГОНАЛЬНОЙ ТЯГИ1
С использованием кватернионов и принципа максимума решена в нели-
нейной постановке задача об оптимальном переводе орбиты космического
аппарата (КА) с переменной массой на заданную плоскость. Управление
движением аппарата производится с помощью ограниченной по модулю
реактивной тяги, ортогональной к плоскости оскулирующей орбиты КА.
Учитывается изменение массы аппарата за счет расхода рабочего тела
на процесс управления. Функционал, определяющий качество процесса
управления, представляет собой линейную свертку с весовыми множите-
лями двух критериев: времени и суммарного импульса тяги, затраченных
на процесс управления.
Излагается теория решения задачи. Приводятся результаты расчетов
оптимального управления для случаев, когда в минимизируемом комби-
нированном функционале качества процесса управления одновременно
учитываются оба критерия, и для случаев, когда минимизируется лишь
суммарный импульс тяги. Получены примеры оптимального управления,
содержащие до 192 пассивных и активных этапов. Установлены законо-
мерности оптимального управления поворотом плоскости орбиты КА.
Ключевые слова: космический аппарат, ориентация орбиты и плоскости
орбиты, ограниченная (малая) реактивная тяга, оптимальное управление,
кватернион.
DOI: 10.1134/S0005231019080087
1. Введение
В работе с использованием кватернионного дифференциального уравне-
ния ориентации орбиты космического аппарата (КА) и принципа максимума
Понтрягина изучается задача оптимальной переориентации плоскости орби-
ты КА с переменной массой посредством ограниченной по модулю реактив-
ной тяги, ортогональной плоскости оскулирующей орбиты. Частным случаем
этой задачи является хорошо известная и имеющая большое практическое
значение задача коррекции угловых элементов орбиты КА, когда изменения
1 Исследование выполнено при финансовой поддержке Российского фонда фундамен-
тальных исследований в рамках научного проекта № 19-01-00205.
87
угловых элементов плоскости орбиты в процессе управления имеют малые
значения. Использование управления, ортогонального плоскости оскулирую-
щей орбиты КА, позволяет корректировать элементы орбиты КА, сохраняя
форму и размеры орбиты КА неизменными. Это ценное свойство изучаемого
процесса переориентации орбиты КА является полезным как при решении за-
дачи коррекции угловых элементов орбиты КА, так и других задач механики
космического полета, например при управлении конфигурацией группировки
спутников.
При решении задачи оптимального управления плоскостью орбиты КА
учитывается изменение массы аппарата за счет расхода рабочего тела на
процесс управления. Функционал, определяющий качество процесса управ-
ления, представляет собой линейную свертку с весовыми множителями двух
критериев: времени и суммарного импульса величины тяги, затраченных на
процесс управления. Для выбранного функционала оптимальное управление
состоит из пассивных этапов, на которых тяга отсутствует, и активных эта-
пов, на которых тяга принимает максимальное значение.
В статье приводится краткий обзор работ по дифференциальным уравне-
ниям ориентации орбиты КА и изучаемой проблеме оптимальной переориен-
тации плоскости орбиты и орбиты КА посредством реактивной тяги, ортого-
нальной плоскости оскулирующей орбиты. Излагается кватернионная нели-
нейная теория решения задачи оптимальной переориентации плоскости орби-
ты КА переменной массы в непрерывной постановке с использованием огра-
ниченной (малой) тяги, ортогональной плоскости орбиты КА, для нефиксиро-
ванного числа включений реактивной тяги. В отличие от большинства других
работ в качестве управления рассматривается не вектор реактивного ускоре-
ния центра масс КА, а вектор реактивной тяги двигателя КА. В состав ис-
пользуемой для решения задачи математической модели, описывающей про-
цесс переориентации плоскости орбиты КА, входят кватернионное дифферен-
циальное уравнение ориентации оскулирующей орбиты КА, дифференциаль-
ное уравнение для истинной аномалии, характеризующей положение центра
масс КА на орбите, и дифференциальное уравнение, описывающее изменение
массы КА в процессе переориентации плоскости его орбиты.
Приводятся результаты расчетов оптимального управления в случаях, ко-
гда учитываются оба критерия оптимальности и когда минимизируется лишь
импульс величины реактивной тяги, что равносильно минимизации расхода
рабочего тела на процесс управления. Показано, что с уменьшением допусти-
мого максимального значения величины тяги увеличивается продолжитель-
ность процесса управления и число этапов в процессе оптимального управле-
ния. Получены примеры оптимального управления, содержащие до 192 пас-
сивных и активных этапов. Анализ результатов расчетов позволил опреде-
лить закономерности в оптимальном управлении плоскостью орбиты КА.
Отметим, что оптимальное управление посредством ортогональной тяги
удобно использовать в основном для решения задачи коррекции плоскости
орбиты КА (коррекции наклонения и долготы восходящего узла) при сохра-
нении формы и размеров орбиты в случаях, когда разности между угловы-
ми элементами начальной и конечной орбиты малы. Однако в статье также
приведены примеры численного решения задачи оптимального управления
88
для немалых разностей угловых элементов орбиты. Это, с одной стороны,
демонстрирует возможности разработанной программы численного решения
задачи (ее работоспособность для широкого диапазона изменения парамет-
ров, входящих в постановку задачи), а с другой стороны, позволяет, более
полно установить закономерности такого управления. Кроме того, в буду-
щем нельзя исключить необходимости решения задачи для случаев больших
изменений значений наклонения и долготы восходящего узла при сохранении
неизменными формы и размеров орбиты КА.
2. Задачи оптимальной переориентации плоскости орбиты КА
Будем считать, что вектор ускорения w центра масс КА от тяги реактив-
ного двигателя во все время управляемого движения КА направлен ортого-
нально плоскости оскулирующей орбиты, т.е. ортогонально радиус-вектору r
и вектору v скорости центра масс КА (коллинеарно вектору c = r × v мо-
мента скорости центра масс КА). Тогда дифференциальные уравнения дви-
жения центра масс КА в ньютоновском гравитационном поле, описывающие
изменение размеров и формы мгновенной орбиты КА, интегрируются, да-
вая уравнение конического сечения. Поэтому управляемое движение центра
масс КА в этом случае описывается дифференциальными уравнениями, ха-
рактеризующими изменение мгновенной ориентации орбиты КА или исполь-
зуемой (например, орбитальной) вращающейся системы координат, в которой
записываются исходные уравнения движения центра масс КА, и дифферен-
циальным уравнением для истинной аномалии, характеризующей положение
центра масс КА на орбите. В процессе такого управления оскулирующая ор-
бита КА поворачивается в пространстве, не изменяя своих размеров и формы.
Движение центра масс КА будем рассматривать в инерциальной системе
координат - геоцентрической экваториальной системе координат OX1X2X3
(X) c началом в центре O притяжения Земли. Ось OX3 этой системы ко-
ординат направлена вдоль оси суточного вращения Земли, оси OX1 и OX2
лежат в плоскости экватора Земли, ось OX1 направлена в точку весенне-
го равноденствия для Земли, ось OX2 дополняет систему до правой тройки
векторов.
Введем также в рассмотрение систему координат ξ, связанную с плоско-
стью и перицентром орбиты КА. Начало этой системы координат находится
в центре O (или в перицентре орбиты), ось ξ1 направлена вдоль радиус-век-
тора перицентра орбиты, ось ξ3 перпендикулярна плоскости орбиты и име-
ет направление постоянного по модулю вектора c момента скорости центра
масс КА относительно центра O, а ось ξ2 образует правую тройку с осями ξ1
и ξ3. Ориентация системы координат ξ в инерциальной системе координат X
характеризует собой ориентацию орбиты КА в инерциальном пространстве и
традиционно задается тремя угловыми оскулирующими элементами орбиты
[1, 2]: долготой восходящего узла Ωu, наклоном орбиты i и угловым расстоя-
нием перицентра от узла ωπ.
Дифференциальные уравнения, описывающие мгновенную ориентацию
орбиты КА в инерциальной системе координат в угловых элементах орбиты
в рассматриваемом случае ортогональности вектора реактивной тяги плос-
89
кости оскулирующей орбиты КА, имеют вид [1, 2]
dΩu/dt = (r/c)w sin(ωπ + ϑ) cosec i,
di/dt = (r/c)w cos(ωπ + ϑ), dωπ/dt = -(r/c)w sin(ωπ + ϑ) ctg i,
(2.1)
dϑ/dt = c/r2, r = p/(1 + e cos ϑ), c = const,
где ϑ - истинная аномалия (угловая переменная, отсчитываемая в плоско-
сти орбиты от ее перицентра и характеризующая положение КА на орбите),
r = |r| - модуль радиус-вектора центра масс КА, p и e - параметр и эксцен-
триситет орбиты, c = |c| = |r × v| - постоянная площадей (модуль вектора
момента скорости v центра масс КА), w - проекция вектора реактивного
ускорения w на направление вектора момента скорости центра масс КА (ал-
гебраическая величина реактивного ускорения, перпендикулярного плоско-
сти оскулирующей орбиты КА), t - время.
Задача оптимального поворота плоскости орбиты КА в угловых перемен-
ных формулируется следующим образом: требуется построить управление w,
переводящее орбиту, изменение ориентации которой описывается уравнения-
ми (2.1), из заданного начального положения
Ωu = Ωu(t0) = Ω0u, i = i(t0) = i0, ωπ = ωπ(t0) = ω0π, i0 = 0
в требуемое конечное положение
Ωu = Ωu(t1) = Ω∗u, i = i(t1) = i, i = 0
при свободном значении углового расстояния перицентра от узла ωπ на пра-
вом конце поворота. При этом должен минимизироваться выбранный функ-
ционал качества процесса поворота плоскости орбиты КА.
Частные случаи задачи переориентации плоскости орбиты КА рассматри-
вались в [3-7]. В [8] рассматривается задача коррекции угловых элементов
орбиты Ωu, i и ωπ с помощью реактивного ускорения, ортогонального плос-
кости оскулирующей орбиты КА, называемого в этой работе “бинормальным
реактивным ускорением”. Задача решается с помощью принципа максимума
и усреднения уравнений. Из усредненных уравнений получен ряд аналити-
ческих соотношений для определения затрат характеристической скорости
в частных случаях коррекции одного или двух элементов орбиты (наклона
орбиты, долготы восходящего узла) при условии малости изменения наклона
орбиты и долготы восходящего узла. По словам авторов [8] уравнения зада-
чи оптимизации в полном объеме не приведены и не проанализированы из-за
большой громоздкости уравнений для сопряженных переменных.
Замечание 1. Модель (2.1), которая описывает ориентацию оскулирую-
щей орбиты КА при ортогональности вектора ускорения к ее плоскости и име-
ет особенность для наклонения i при значениях i = 0 и i = π, является частью
системы дифференциальных уравнений для всех шести кеплеровых оску-
лирующих элементов орбиты при произвольном векторе тяги двигательной
установки [9, стр. 97]. При этом, наряду с указанными особыми значениями
90
наклонения, здесь имеется особое значение эксцентриситета e = 1, когда ор-
бита становится круговой и исчезает возможность однозначного определения
углового расстояния перицентра от узла ωπ. Решение этой проблемы получе-
но в начале 1970-х гг. с применением двух векторов в пространстве состояний
орбитального движения КА - вектора эксцентриситета e = {e cos ϖ, e sin ϖ},
где {·} ≡ col (·), ϖ = Ωu + ωπ, и вектора наклонения i = {i cos Ωu, i sin Ωu}
[10, 11]. Этот метод, впервые реализованный в системе навигации и управ-
ления движением советского спутника связи Радуга на геостационарной ор-
бите в 1974 г. [12], в настоящее время модифицирован и успешно применяется
на борту российских спутников связи, навигации и геодезии [13]. Проблема
вырожденности классических орбитальных элементов частично решается в
механике космического полета за счет использования так называемых «невы-
рожденных» орбитальных элементов с элементами Лагранжа P1 = e sin ϖ,
P2 = ecos ϖ и Q1 = tg(i/2)sin Ωu, Q2 = tg(i/2)cos Ωu (иногда для них исполь-
зуют термин “equinoctial elements”) и соответствующих уравнений ориентации
орбиты R. Battin [14, 15]. В этих уравнениях сохраняется лишь особое значе-
ние наклонения i = π.
Решение задач оптимальной переориентации плоскости орбиты и орбиты
космического аппарата посредством реактивной тяги, ортогональной плоско-
сти оскулирующей орбиты, с помощью уравнений (2.1) в классических угло-
вых элементах орбиты в строгой нелинейной постановке достаточно сложно
в силу нелинейности этих уравнений, наличия в них особых точек i = 0, π,
а также в силу громоздкости уравнений для сопряженных переменных. По-
этому для решения задачи оптимальной переориентации плоскости орбиты
и орбиты КА вместо угловых элементов орбиты целесообразно использовать
[16-24] параметры Эйлера (Родрига-Гамильтона).
Дифференциальные уравнения ориентации орбиты КА в параметрах Эй-
лера имеют вид [16-20]
2dΛ0/dt = -Ω1Λ1 - Ω2Λ2,
2dΛ1/dt = Ω1Λ0 - Ω2Λ3,
(2.2)
2dΛ2/dt = Ω2Λ0 + Ω1Λ3,
2dΛ3/dt = Ω2Λ1 - Ω1Λ2;
(2.3)
dϑ/dt = c/r2
,
r = p/(1 + ecosϑ), c = const,
Ω1 = (r/c)w cos ϑ, Ω2 = (r/c)w sin ϑ,
где Λj (j = 0, 3) - параметры Эйлера, характеризующие ориентацию ор-
биты КА (системы координат ξ) в инерциальной системе координат X;
Ω1,Ω2,Ω3 = 0 - проекции вектора Ω мгновенной абсолютной угловой ско-
рости орбиты на связанные с ней координатные осиi.
Отметим, что уравнения ориентации орбиты в параметрах Эйлера вида
(2.2) использовались для описания орбитального движения и другими авто-
рами [25, 26].
Параметры Эйлера Λj связаны с угловыми элементами орбиты соотноше-
ниями [17]
Λ0 = cos(i/2)cos((Ωu + ωπ)/2), Λ1 = sin(i/2)cos((Ωu - ωπ)/2),
(2.4)
Λ2 = sin(i/2)sin((Ωu - ωπ)/2), Λ3 = cos(i/2)sin((Ωu + ωπ)/2).
91
Соотношения (2.4) получены с использованием кватернионных формул
сложения конечных поворотов [17, 27, 28] (на основе умножения кватернио-
нов, соответствующих трем элементарным поворотам на углы Ωu, i и ωπ в
указанной последовательности).
Угол наклона орбиты i определяется через параметры Эйлера по формуле
Λ21 + Λ22
i = 2arctg
Λ20 + Λ2
3
Для определения угловых элементов Ωu и ωπ сначала вычисляются вели-
чиныΩ и ω по формулам
Λ1Λ3 - Λ0Λ2
Ω= arctgΛ0Λ21Λ3,
ω = arctg
,
Λ0Λ1 - Λ2Λ3
Λ0Λ1 + Λ2Λ3
а затем вычисляются Ωu и ωπ по следующему однозначному алгоритму:
а) еслиΩ 0 и Λ0Λ2 + Λ1Λ3 0, то Ωu =Ω, а если Λ0Λ2 + Λ1Λ3 < 0, то
Ωu =Ω + π;
б) еслиΩ < 0 и Λ0Λ1 - Λ2Λ3 0, то Ωu =Ω + π, а если Λ0Λ1 - Λ2Λ3 > 0,
то Ωu =Ω + 2π;
в) если ω ≥ 0 и Λ1Λ3 - Λ0Λ2 0, то ωπ = ω, а если Λ1Λ3 - Λ0Λ2 < 0, то
ωπ = ω + π;
г) если ω < 0 и Λ0Λ1 + Λ2Λ3 0, то ωπ = ω + π, а если Λ0Λ1 + Λ2Λ3 > 0,
то ωπ = ω + 2π.
Эти формулы получены из соотношений (2.4) с учетом диапазонов изме-
нения углов i, Ωu и ωπ: 0 ≤ i ≤ π, 0 Ωu 2π, 0 ≤ ωπ 2π. Такие формулы
применяются здесь для вычисления значений этих углов в конечный момент
времени процесса управления через параметры Λj.
Используя вектор параметров Эйлера Λ = {Λ0, Λv}, Λv = {Λi}, i = 1, 2, 3,
уравнения (2.2) можно представить в векторном виде
dΛ0/dt = - · Λv)/2, dΛv/dt = (Λ0Ω + Ω × Λv)/2,
Ω ≡ {Ωi} = (r/c)w{cos ϑ,sinϑ,0}.
Уравнения (2.2) в кватернионной записи принимают компактный вид
[16-20] (здесь и далее применяются общепринятые обозначения операций с
кватернионами [17, 18, 27, 28])
(2.5)
2dΛ/dt = Λ Ωξ, Ωξ = Ω1i1 + Ω2i2 = (r/c)w(cos ϑ i1 + sin ϑ i2
),
где Λ = Λ0 + Λ1i1 + Λ2i2 + Λ3i3 - кватернион ориентации орбиты КА (ква-
тернионный оскулирующий (медленно изменяющийся) элемент орбиты КА);
Ωξ - отображение вектора Ω на базис ξ (вектор Ω мгновенной абсолютной
угловой скорости орбиты направлен вдоль радиус-вектора r центра масс КА
92
и определяется формулой: Ω = (w/c)r); i1, i2, i3 - векторные мнимые единицы
Гамильтона, - символ кватернионного умножения.
Отметим, что уравнения (2.2), (2.3) или (2.5), (2.3) - система пяти нели-
нейных стационарных дифференциальных уравнений первого порядка отно-
сительно параметров Эйлера Λj и истинной аномалии ϑ. Эти уравнения, в
отличие от четырех нелинейных дифференциальных уравнений (2.1) ориен-
тации орбиты в угловых элементах орбиты Ωu, i, ωπ, не имеют особых точек
i = 0, π, к тому же при переходе в них от времени t к новой независимой пе-
ременной ϑ в соответствии с дифференциальным соотношением = (c/r2)dt
получаем (при w = w(ϑ)) систему четырех линейных нестационарных диф-
ференциальных уравнений относительно параметров Эйлера Λj (в то время
как дифференциальные уравнения в угловых элементах орбиты остаются су-
щественно нелинейными).
Указанные обстоятельства делают использование уравнений (2.2), (2.3)
или (2.5), (2.3) для решения задач переориентации орбиты и ее плоскости
с помощью реактивной тяги, ортогональной плоскости оскулирующей ор-
биты КА, более удобным и эффективным в сравнении с использованием
уравнений в угловых оскулирующих элементах (2.1). Такое решение зада-
чи переориентации орбиты КА постоянной массы в инерциальной системе
координат в непрерывной постановке (с использованием в качестве управле-
ния вектора реактивного ускорения центра масс КА от ограниченной (малой)
тяги и принципа максимума) рассмотрено в [21, 22]. В [23, 24] эта задача изу-
чается в кватернионной постановке с помощью импульсной реактивной тяги.
Исследованию задачи оптимальной переориентации орбиты КА постоян-
ной массы в непрерывной постановке (с использованием в качестве управ-
ления вектора реактивного ускорения центра масс КА) и с использовани-
ем кватернионного дифференциального уравнения ориентации орбитальной
системы координат посвящены работы [19, 29, 30]. Использование кватерни-
онного дифференциального уравнения ориентации орбитальной системы ко-
ординат более удобно при аналитическом исследовании задачи оптимальной
переориентации орбиты КА в непрерывной постановке, однако использова-
ние кватернионного дифференциального уравнения ориентации орбиты КА
имеет преимущество при численном решении задачи оптимальной переориен-
тации орбиты КА и ее плоскости, так как кватернион ориентации орбиты КА
является оскулирующим (медленно изменяющимся) элементом орбиты. Ква-
тернион ориентации орбитальной системы координат таким свойством не об-
ладает, так как является быстро меняющейся переменной.
Отметим также, что уравнения Battin [14, 15], также как и используемые
здесь уравнения (2.2), (2.3) (или (2.5), (2.3)) ориентации орбиты КА в пара-
метрах Эйлера, не имеют особой точки (деления на ноль) при i = 0. Однако
уравнения Battin и сопряженные к ним уравнения значительно сложнее ква-
тернионных регулярных фазовых уравнений (2.2), (2.3) (или (2.5), (2.3)) и
ниже приводимых сопряженных уравнений решаемой задачи как с анали-
тической, так и с вычислительной точек зрения. К тому же кватернионное
фазовое уравнение (2.5) (кватернионное уравнение ориентации орбиты КА,
эквивалентное четырем скалярным уравнениям (2.2)) обладает свойством са-
мосопряженности: оно с точностью до обозначения кватернионной перемен-
93
ной совпадает с кватернионным сопряженным уравнением, что позволяет по-
низить размерность краевой задачи оптимизации на 4 единицы с использо-
ванием новой кватернионной переменной, являющейся мультипликативной
композицией кватернионных фазовой и сопряженной переменных (в виде их
произведения). Таким свойством классические дифференциальные уравне-
ния (2.1) ориентации орбиты в угловых элементах орбиты Ωu, i, ωπ и уравне-
ния Battin не обладают, причем соответствующие им сопряженные уравнения
гораздо сложнее фазовых.
3. Кватернионная постановка задачи оптимальной
переориентации плоскости орбиты КА
Изменение ориентации орбиты КА переменной массы с помощью реактив-
ной тяги, ортогональной плоскости оскулирующей орбиты КА, описывается
системой дифференциальных уравнений
2dΛ/dt = Λ Ωξ, Ωξ = (r/c)(u/m)(cos ϑ i1 + sin ϑ i2),
(3.1)
c
dm
p
=
,
= -|u, r =
,
dt
r2
dt
1 + ecosϑ
где m - масса аппарата, u - проекция вектора реактивной тяги u на на-
правление вектора момента скорости центра масс КА (алгебраическая ве-
личина реактивной тяги, перпендикулярной плоскости оскулирующей орби-
ты КА), β - коэффициент пропорциональности, равный обратной величине
скорости истечения рабочего тела из сопла реактивного двигателя [31].
На величину тяги u наложено ограничение
(3.2)
|u| ≤ u∗m.
Компоненты Λj (j = 0, 3) кватерниона ориентации орбиты Λ выражаются
через угловые элементы орбиты Ωu, i и ωπ по формулам (2.4).
Перейдем к безразмерным переменным ρ, τ, m, u, um, β по формулам
2
p
c2m0
c2m0
p
(3.3)
r = pρ, t =
τ, ,m = m0m, u =
u, u∗m =
um, β =
β,
c
p3
p3
c
где m0 - начальная масса аппарата.
Уравнения движения КА (3.1) в безразмерных переменных примут вид
dΛ
u
=
Λ (i1 cosϑ + i2 sinϑ),
2(1 + ecos ϑ)m
(3.4)
= (1 + e cos ϑ)2 ,dm
= -β |u| .
На управляющий параметр - безразмерную тягу u наложено ограничение
(3.5)
|u| ≤ um.
94
В начальный момент времени состояние управляемой системы (3.4) опре-
деляется соотношениями
(3.6)
ϑ=ϑn, Λ=Λn
,
m = 1,
где компоненты кватерниона ориентации орбиты Λn выражаются однознач-
но через угловые элементы начальной орбиты КА in , Ωun, ωπn по форму-
лам (2.4).
Требуется перевести КА на новую орбиту, плоскость которой определяется
угловыми элементами ik, Ωuk.
Таким образом, при оптимизации поворота плоскости орбиты краевое на-
чальное условие Λn определяется однозначно согласно (2.4), но при назначе-
нии правого краевого условия Λk значение углового расстоянием перицентра
от узла ωπ является свободным.
В фазовом пространстве ϑ × Λ × m многообразие, на которое требуется
перевести управляемую систему (3.4), определяется соотношениями
ik
ik
ik
Λ0 sinΩuk sin
- Λ2 cos
- Λ3 cosΩuk sin
= 0,
2
2
2
(3.7)
ik
ik
ik
Λ1 sin Ωuk cos
- Λ3 sin
- Λ2 cosΩuk cos
= 0,
2
2
2
получающимися после подстановки угловых элементов орбиты ik, Ωuk в (2.4)
и последующего исключения из них незаданного ωπk. Для этого необходимо
из соотношений, определяющих Λ2 и Λ3, выразить sin(ωπk/2) и cos(ωπk/2) и
подставить их в соотношения для Λ0 и Λ1.
Качество процесса перевода КА на новую орбиту будем определять мини-
мальным значением функционала
τk
(3.8)
J = (α1 + α2 |u|)dτ, α10, α2
0,
0
который представляет собой линейную свертку с постоянными весовыми мно-
жителями α1 и α2 двух критериев: времени и импульса величины тяги, за-
траченных на процесс управления. Изменяя величины весовых множителей
в функционале (3.8), можно усиливать влияние того или иного критерия на
процесс управления.
Таким образом, требуется найти оптимальное управление безразмерной
реактивной тягой u = u (τ), удовлетворяющее ограничению (3.5), которое пе-
реводит управляемую систему (3.4) из начального состояния (3.6) на много-
образие (3.7) и сообщает функционалу (3.8) минимальное значение.
Промежуток времени процесса управления заранее не задается. Ниже рас-
сматривается случай, когда α2 > 0, чтобы учесть влияние импульса величины
тяги на качество процесса управления.
95
4. Решение задачи с помощью принципа максимума Понтрягина
Функция Гамильтона для поставленной задачи оптимального управления
имеет вид
H = -(α1 + α2 |u|) + χ(1 + ecosϑ)2 +
(4.1)
u
+
(M, Λ (i1 cos ϑ + i2 sin ϑ)) - ηβ |u| .
2(1 + e cos ϑ)m
Здесь и далее для кватернионов a = a0 + av и b = b0 + bv использует-
ся аналог их скалярного произведения (a, b) = scal((a0 - av) (b0 + bv)) =
= a0b0 + a1b1 + a2b2 + a3b3, где aj и bj - компоненты кватернионов a и b (a0,b0
и av,bv - скалярная и векторная части кватернионов a и b).
Сопряженные переменные χ, M = M0 + M1i1 + M2i2 + M3i3, η, соответ-
ствующие фазовым переменным ϑ, Λ, m, удовлетворяют системе уравнений
= 2(1 + e cos ϑ) sin ϑ +
(4.2)
u
+
(M, Λ (i1 sin ϑ - i2(e + cos ϑ))) ,
2(1 + e cos ϑ)2m
dM
u
(4.3)
=
M (i1 cosϑ + i2
sin ϑ),
2(1 + e cos ϑ)m
u
(4.4)
=
(M, Λ (i1 cos ϑ + i2
sin ϑ)) .
2(1 + e cos ϑ)m2
Введем новую кватернионную переменную
(4.5)
N=N0 +N1i1 +N2i2 +N3i3 ≡N0 +Nv =Λ
M,
компоненты которой вычисляются по явным соотношениям (Λ - сопряжен-
ный кватернион:Λ = Λ0 - Λ1i1 - Λ2i2 - Λ3i3).
Функция H и уравнения (4.2), (4.4), которым удовлетворяют сопряженные
переменные χ и η, примут вид
H = -(α1 + α2 |u|) + χ(1 + ecosϑ)2 +
(4.6)
u
+
(N1 cos ϑ + N2 sin ϑ) - ηβ |u| ,
2(1 + e cos ϑ)m
= 2(1 + e cos ϑ) sin ϑ +
(4.7)
u
+
(N1 sin ϑ - N2(e + cos ϑ)) ,
2(1 + e cos ϑ)2m
u
(4.8)
=
(N1 cos ϑ + N2
sin ϑ) .
2(1 + e cos ϑ)m2
96
Из условия максимума для функции H следует, что оптимальное управ-
ление определяется через фазовые и сопряженные переменные по формулам
{ umsign(νu), если εu 0;
(4.9)
u=
u = 0,
если εu < 0,
где νu = N1 cos ϑ + N2 sin ϑ и εu =u| / [2m(1 + e cos ϑ)] - α2 - ηβ.
Особый режим управления в настоящей работе не рассматривается.
Из первого уравнения системы (3.4) и из соотношений (4.9) следует, что
на активном этапе управления орбита КА вращается с безразмерным векто-
ром угловой скорости Ωξd = umsign(νu) (i1 cos ϑ + i2 sin ϑ) /[m (1 + e cos ϑ)], а
на пассивном этапе ориентация орбиты не изменяется.
Так как правый конец траектории в фазовом пространстве ϑ × Λ × m яв-
ляется подвижным, то в конце движения, т.е. в момент τ = τk, должны выпол-
няться условия трансверсальности. На истинную аномалию ϑ и на массу m
в конце движения не налагается никаких условий, поэтому соответствующие
им сопряженные переменные χ и η при τ = τk должны обращаться в нули,
т.е.
(4.10)
χ(τk) = 0, η (τk
) = 0.
Из соотношений (3.7) с помощью метода неопределенных множителей
Лагранжа [32] выводятся условия трансверсальности для сопряженной ква-
тернионной переменной M
ik
M0 = λ1 sinΩuk sin
,
2
ik
M1 = λ2 sinΩuk cos
,
2
(4.11)
ik
ik
M2 =1 cos
- λ2 cosΩuk cos
,
2
2
ik
ik
M3 =1 cos Ωuk sin
- λ2 sin
2
2
Если из первых двух соотношений формул (4.11) выразить неопределен-
ные множители λ1 и λ2 через M0 и M1 и подставить их в третье и четвертое
соотношения этих формул, то в результате получатся условия трансверсаль-
ности, свободные от множителей Лагранжа:
ik
ik
ik
M2 sin Ωuk sin
+ M0 cos
+ M1 cosΩuk sin
= 0,
2
2
2
(4.12)
ik
ik
ik
M3 sin Ωuk cos
+ M1 sin
+ M0 cosΩuk cos
= 0.
2
2
2
Из соотношений (3.7) и (4.12) можно получить, что при τ = τk
(4.13)
(M, Λ) = 0.
97
Для этого из первого соотношения системы (3.7) находится cos(ik/2) и
подставляется в первое соотношение системы (4.12). Аналогично, из второго
соотношения системы (3.7) находится sin(ik/2) и подставляется во второе со-
отношение системы (4.12). Соотношения, полученные из системы (4.12), при-
водятся к общему знаменателю. Далее первое из полученных соотношений
умножается на cos(ik/2), а второе умножается на sin(ik/2), левые и правые
части полученных уравнений складываются. В результате слагаемые, в ко-
торые входят произведения фазовых и сопряженных переменных с разными
индексами, уничтожаются, а произведения с одинаковыми индексами скла-
дываются. В итоге получается соотношение (4.13).
Условием (4.13) можно заменить одно из условий (4.12). Для уравнений
системы (3.4) и (4.3) соотношение (M, Λ) = c = const является первым ин-
тегралом, причем согласно (4.13) значение c = 0. По этой причине соотно-
шение (4.13) имеет место во все время процесса управления. Тогда согласно
(4.5) и (4.13) компонента N0 0.
Для компонент N1, N2, N3 кватерниона N можно получить на основе (3.4),
(4.3), (4.5) систему дифференциальных уравнений
dN1
u sin ϑ
dN2
u cos ϑ
=-
N3,
=
N3,
(1 + e cos ϑ)m
(1 + e cos ϑ)m
(4.14)
dN3
u
=
(N1 sin ϑ - N2 cos ϑ) .
(1 + e cos ϑ)m
Сопряженную кватернионную переменную M можно исключить из рас-
смотрения, если вместо уравнения (4.3) в постановке краевой задачи опти-
мального управления использовать систему уравнений (4.14). Так как ска-
лярная составляющая N0 кватерниона N равна нулю, то можно ввести век-
тор Nv с координатами N1, N2, N3, который представляет собой векторную
часть кватерниона N. Тогда систему дифференциальных уравнений (4.14)
можно записать в виде векторного дифференциального уравнения
(4.15)
dNv/dτ = (u/m)(Nv × (cos ϑ i1 + sin ϑ i2
))/(1 + e cos ϑ).
С помощью (4.5) и с учетом того, что N0 = 0, второе условие трансвер-
сальности из системы (4.12) можно переписать в виде
0N3 + Λ1N2 - Λ3N1) sinΩuk cos(ik/2) -
(4.16)
-1N1 + Λ2N2 + Λ3N3)cos Ωuk cos(ik/2) +
+ (Λ0N1 + Λ2N3 - Λ3N2) sin(ik/2) = 0.
Так как время момента окончания процесса не задается заранее, то в ко-
нечный момент времени функция Гамильтона (4.6) равна нулю:
(4.17)
Hopt(τk
) = 0.
Отметим, что Hopt(τ) = 0 является первым интегралом изучаемой задачи.
98
Таким образом, принцип максимума Понтрягина сводит поставленную за-
дачу оптимального управления к решению краевой задачи для системы диф-
ференциальных уравнений (3.4), (4.7), (4.8), (4.14) (или (4.15)) 11-го порядка
(в скалярном исчислении) с граничными условиями (3.6) в начальный момент
времени и условиями (3.7), (4.10), (4.16), (4.17) в конце процесса управления
(в количестве, равном 12). При этом в каждый момент времени оптималь-
ное управление определяется через фазовые и сопряженные переменные из
соотношений (4.9).
Подчеркнем, что в рассматриваемой задаче оптимальной переориентации
плоскости орбиты КА с переменной массой управляющим параметром явля-
ется тяга реактивного двигателя КА, в то время как в работах [21-24, 29,
30] авторов статьи и их коллег управляющим параметром являлось ускоре-
ние КА от тяги реактивного двигателя, к тому же в этих работах рассмат-
ривалось решение в различных постановках задачи оптимальной переори-
ентации орбиты КА с постоянной массой, а не плоскости орбиты КА с пе-
ременной массой. Формула (4.9) показывает, что на структуру оптимально-
го управления существенно влияет изменение массы КА в процессе управ-
ления и скорость истечения рабочей массы из сопла двигателя, связанная
с параметром β, фигурирующим в законе управления, обратной зависимо-
стью. В дифференциальных фазовых и сопряженных уравнениях, а также
в дифференциальных уравнениях линии переключения управления (4.14) в
рассматриваемой задаче присутствует в явном виде переменная масса КА,
что оказывает качественное влияние на процесс переориентации плоскости
орбиты КА, свойства которого изучаются в следующем разделе.
5. Результаты расчетов и их анализ
Сформулированная краевая задача решалась численно с помощью комби-
нации модифицированного метода Ньютона и метода градиентного спуска.
Для улучшения сходимости итерационного процесса производилось уточне-
ние моментов разрыва управления по специальному алгоритму, описанному
ниже. Оптимальное управление состоит из нескольких этапов: активных и
пассивных. На активных этапах управляющий параметр принимает значе-
ния u = um или u = -um (на активных этапах происходит изменение ориен-
тации орбиты). На пассивных этапах u = 0 (на пассивных этапах ориентация
орбиты не изменяется и совпадает с ориентацией орбиты, соответствующей
ориентации орбиты конца предыдущего этапа).
Численное решение краевой задачи для системы дифференциальных урав-
нений связано с большими трудностями, так как приходится подбирать недо-
стающие начальные условия, а также промежуток времени процесса управ-
ления. На сходимость итерационного процесса при решении краевой задачи
оптимизации влияет выбор начального приближения. Дополнительные труд-
ности со сходимостью процесса численного решения возникают при наличии
разрывов управляющих параметров. В рассматриваемой задаче при малых
значениях максимальной тяги число этапов, а следовательно, и число раз-
рывов управления возрастает. По этой причине необходимо, чтобы алгоритм
решения задачи учитывал наличие разрывов управления.
99
Таблица 1
um
τk
ϑk
mk
J
ωπk
net
0,25
1,720485
130,8325
0,954649
0,781403
64,9135
2
0,125
2,301065
158,9216
0,947368
0,970003
65,2837
2
0,075
4,786863
283,8775
0,940712
1,640374
64,8113
4
0,05
7,277088
86,7048
0,935770
2,301095
64,7793
6
При решении краевой задачи методом Ньютона или методом градиентно-
го спуска необходимо вычислять невязки выполнения краевых условий и их
производных по начальным условиям сопряженных переменных и по проме-
жутку времени процесса. Невязки будут точнее вычисляться при интегриро-
вании системы дифференциальных уравнений методом Рунге - Кутты, если
моменты разрыва управляющего параметра не будут оказываться внутри ша-
га интегрирования. По этой причине на каждом шаге интегрирования надо
проверять, был ли разрыв управления внутри шага. Для этого вычисляется
согласно (4.9) в левой и правой точке шага выражение
εu =u|/[2m(1 + ecos ϑ)] - α2 - ηβ.
Если знаки выражения одинаковые, то разрыв управления внутри шага
интегрирования отсутствовал. Если же знаки разные, то разрыв был. На
шаге, внутри которого был разрыв управления, производится интегрирование
с более мелким шагом. Мелкий шаг, в котором окажется разрыв, делится на
два пропорциональных шага интегрирования, разделенных точкой разрыва
управления.
В примерах, результаты расчетов для которых приводятся в таблицах 1-5,
начальное состояние системы (3.4) определяется параметрами
(5.1)
e = 0,1, ϑn = 30,0, in = 7,0, Ωun = 30,0, ωπn = 50,0, mn
= 1,0.
Ориентация плоскости конечной орбиты определяется углами
(5.2)
ik = 20,0, Ωuk = 15,0.
Примеры различаются максимальными значениями тяги um, значениями
коэффициента β, который определяет отношение характерного значения ско-
рости движения КА по орбите к скорости истечения рабочего тела из сопел
двигателя КА, и значениями весовых множителей α1, α2 в функционале (3.8).
Результаты расчетов представлены в безразмерных переменных, углы в гра-
дусах.
В табл. 1 для значений α1 = 0,25, α2 = 1,5, β = 0,2 и различных значе-
ний um приведены значения длительности процесса управления (вторая ко-
лонка), значения истинной аномалии КА и массы КА в конце процесса (тре-
тья и четвертая колонки соответственно), значения функционала качества
процесса (пятая колонка), значения углового расстояния до перицентра ор-
биты в конце процесса управления (шестая колонка) и число этапов net
в процессе управления (седьмая колонка). Из расчетов, представленных в
100
Таблица 2
β
τk
ϑk
mk
J
ωπk
net
0,2
5,491449
329,0163
0,949182
0,254090
64,8251
4
0,5
5,445041
325,8902
0,879008
0,241984
64,8136
4
1,0
5,377398
321,3655
0,775613
0,224387
64,7971
4
Таблица 3
β
τk
ϑk
mk
J
ωπk
net
0,2
8,586485
154,5172
0,949954
0,250230
64,8388
6
0,5
8,535597
152,0992
0,880727
0,238546
64,8316
6
1,0
8,461673
148,5677
0,778526
0,221474
64,8215
6
табл. 1, видно, что с уменьшением максимального значения тяги, которая
может создаваться двигательной установкой аппарата, увеличивается про-
должительность процесса управления по переводу орбиты КА на заданную
плоскость, незначительно увеличивается расход рабочего тела (топлива) и со-
ответственно уменьшается масса аппарата в конце процесса, увеличиваются
значения функционала качества процесса управления, увеличивается число
этапов процесса управления от двух при um, равном 0,25 или 0,125, до че-
тырех при um = 0,075 и до шести при um = 0,05. Для всех примеров первый
этап является пассивным. Важно отметить, что угловое расстояние до пе-
рицентра орбиты в конце процесса управления слабо зависит от величины
максимальной тяги um. Отсюда следует, что при всех значениях um аппарат
в результате процесса перевода его на заданную плоскость орбиты каждый
раз оказывается вблизи одного и того же положения на орбите.
Проводились расчеты для других значений весовых множителей в функ-
ционале качества (3.8). Отметим, что если число этапов в процессе управле-
ния равно двум, то в этом случае оптимальный процесс управления не за-
висит от значений весовых множителей в функционале (3.8). В этом случае
от весовых множителей зависят лишь сопряженные переменные. Это связано
с тем, что в этом случае для определения двух моментов окончания этапов
(пассивного и активного) имеются два условия: (2.4) (при Ωu = Ωun, i = in
и ωπ = ωπn) и (3.7). Благодаря этому краевую задачу для системы диффе-
ренциальных уравнений для фазовых переменных можно решать без урав-
нений для сопряженных переменных. По этой причине далее рассматрива-
ются решения задач оптимального управления для um = 0,075, um = 0,05 и
um = 0,025, для которых при заданных граничных условиях (5.1) и (5.2) чис-
ло этапов управления больше двух.
Ниже в виде таблиц приводятся результаты решения задачи оптимального
управления для α1 = 0,0 и α2 = 1,0 и нескольких значений коэффициента β.
В этой задаче требуется перевести КА из начального состояния (5.1) на орби-
ту, плоскость, которой определяется условиями (5.2), с минимальным значе-
нием импульса величины тяги двигателя (с минимальным расходом рабочего
тела). В табл. 2 приводятся результаты решения задачи для um = 0,075, а в
табл. 3 - для um = 0,05. В отличие от табл. 1 в табл. 2 и 3 в первой колонке
приводятся значения коэффициента β.
101
Таблица 4
№ u
τ
i
Ω◦u
ω◦π
ϑ
M
1
0,0
0,356075
7,0
30,0
50,0
53,5317
1,0
2
-0,025
2,375833
9,5604
25,4348
54,5403
162,4352
0,949506
3
0,0
3,788566
9,5604
25,4348
54,5403
229,3201
0,949506
4
0,025
5,557657
11,6824
19,9718
59,9271
333,5040
0,905279
5
0,0
6,764364
11,6824
19,9718
59,9271
55,4346
0,905279
6
-0,025
8,716633
14,5198
18,5903
61,2971
160,6597
0,856472
7
0,0
10,202816
14,5198
18,5903
61,2971
231,1084
0,856472
8
0,025
11,903208
16,9037
15,4391
64,3495
331,2583
0,813962
9
0,0
13,173439
16,9037
15,4391
64,3495
57,37776
0,813962
10
-0,025
15,056534
20,0
15,0
64,7932
158,8377
0,766885
Таблица 5
um
τk
J = Δmk
ϑk
nob net
ωπk
0,075
5,377398
0,224387
321,3655
0
4
64,7971
0,05
8,461673
0,221474
148,5677
1
6
64,8215
0,025
15,056534
0,233115
158,8377
2
10
64,7932
0,0125
31,069513
0,235341
332,7422
4
20
64,7170
0,005
75,758654
0,239085
335,9647
11
48
64,6988
0,0025
152,301870
0,239065
335,9453
23
96
64,6912
0,00125
305,388734
0,239055
335,9358
47
192
64,6873
Из табл. 2 и 3 видно, что с увеличением максимальной тяги um и с умень-
шением коэффициента β повышается эффективность двигательной установ-
ки КА, так как при этом увеличивается конечная масса аппарата и, следова-
тельно, уменьшаются затраты рабочего тела на процесс управления. С увели-
чением коэффициента β уменьшается значение функционала качества про-
цесса. Наиболее значительно изменение коэффициента β влияет на поведение
массы аппарата. На время процесса управления, на угловое расстояние до пе-
рицентра в конце процесса изменение коэффициента β влияет незначительно.
Для тех же граничных условий (5.1), (5.2) и α1 = 0,0, α2 = 1,0, β = 1,0 при
условии, что максимальное значение тяги um = 0,025, число этапов в опти-
мальном управлении оказалось равным десяти. В табл. 4 приводятся резуль-
таты расчета для каждого этапа управления. В первой колонке указывается
номер этапа управления, во второй - управление на этапе, в третьей колон-
ке по восьмую - момент времени, угол наклона орбиты, угол восходящего
узла, угловое расстояние до перицентра, истинная аномалия и масса аппа-
рата в конце каждого этапа соответственно. Значение функционала качества
J = 0,233115.
В табл. 5 приводятся для различных значений максимальной тяги двига-
теля КА результаты расчета оптимального управления поворотом плоскости
орбиты КА при тех же начальных условиях и начальной ориентации орби-
ты (5.1) и той же ориентации плоскости конечной орбиты (5.2), что и выше,
при условии, что α1 = 0,0, α2 = 1,0, β = 1,0. В первой колонке табл. 5 указано
102
Таблица 6
um
τk
J = Δmk
ϑk
nob net
ωπk
0,0025
12,246738
0,020574
354,8571
1
8
50,9921
0,001
31,382276
0,020571
354,8337
4
20
50,9918
0,0005
63,275304
0,020570
354,8270
9
40
50,9916
Таблица 7
um
τk
J = Δmk
ϑk
nob net
ωπk
0,0025
11,577105
0,020536
356,9952
1
7
49,0069
0,001
30,707830
0,020500
356,6382
4
19
49,0067
0,0005
62,599522
0,020491
356,5389
9
39
49,0064
0,00025
126,385110
0,020486
356,4930
19
79
49,0064
Таблица 8
β
τk
mk
ϑk
Δmk
J
ωπk
0,2
126,400366
0,995848
357,5861
0,004162
0,020810
49,0064
0,5
126,394885
0,989673
357,1705
0,010327
0,020654
49,0065
1,0
126,385110
0,979514
356,4930
0,020486
0,020486
49,0064
максимальное значение допустимой тяги, для которой производится расчет,
во второй колонке - продолжительность процесса управления, в третьей ко-
лонке - значение функционала качества процесса, который определяет расход
рабочего тела на процесс управления Δmk, в четвертой колонке - значение
истинной аномалии, определяющей положение КА на орбите, в пятой и ше-
стой колонках указано число полных витков nob, выполненных КА в процессе
движения, и число этапов net в оптимальном управлении соответственно, в
седьмой колонке указано значение углового расстояния до перицентра орби-
ты в конце процесса управления.
Из табл. 5 видно, что наименьший расход рабочего тела наблюдается при
um = 0,05, угловое расстояние до перицентра орбиты в конце процесса слабо
зависит от значения максимальной тяги.
Ниже приводятся результаты решения задачи об оптимальном повороте
плоскости орбиты при условии, что плоскости начальной и конечной орбит
близки друг к другу. Различия между соответствующими угловыми элемен-
тами составляют один градус. Расчеты проводились при значениях α1 = 0,0,
α2 = 1,0, β = 1,0, e = 0,1. В начальный момент положение орбиты и поло-
жение КА на орбите определяется условиями (5.1). Ориентация плоскости
конечной орбиты определяется углами
(5.3)
ik = 8,0, Ωuk = 29,0.
В табл. 6 представлены результаты расчетов для трех значений макси-
мальной тяги. На первом этапе пассивное движение, на втором этапе u =
= -um, на третьем этапе пассивное движение, на четвертом этапе u = um и
далее этапы повторяются.
103
Таблица 9
um
τk
J = Δmk
ϑk
nob net
ωπk
0,00025
27,665243
0,001740
151,8025
4
18
49,9008
0,0001
59,646801
0,001788
156,0003
9
38
49,9007
0,00005
120,127999
0,001794
317,0792
18
76
49,9007
0,000025
241,419823
0,001792
317,0605
37
148
49,9007
В табл. 7 приводятся расчеты с другим начальным значением истинной
аномалии, когда начальные условия определяются параметрами
(5.4)
e = 0,1, ϑn = 75,0, in = 7,0, Ωun = 30,0, ωπn = 50,0, mn
= 1,0,
а конечные условия (5.3) сохраняются. В этом случае первый этап активный
u = um, второй этап пассивный, третий этап активный u = -um, четвертый
этап пассивный и далее этапы повторяются. За счет увеличения начальной
аномалии, характеризующей положение КА на начальной орбите, первый
этап в отличие от приведенных ранее примеров оказался активным.
Из табл. 6 и 7 видно, что с уменьшением максимального значения тя-
ги существенно увеличивается продолжительность процесса и незначительно
уменьшается расход рабочего тела (в пятом знаке).
В табл. 8 для различных значений коэффициента β приводятся результа-
ты расчета оптимального поворота плоскости орбиты КА для случая, когда
um = 0,00025, α1 = 0,0, α2 = 1,0, начальные условия определяются соотноше-
ниями (5.4), а ориентация плоскости конечной орбиты - соотношениями (5.3).
В этом случае КА совершает 19 полных витков, оптимальное управление со-
держит 39 этапов, первый этап активный u = um, второй этап пассивный,
третий этап активный u = -um, четвертый этап пассивный и далее этапы
повторяются, на последнем этапе u = -um.
Из табл. 8 видно, что от значения коэффициента β существенно зависит
только расход рабочего тела на процесс управления Δmk, другие показатели
процесса меняются незначительно.
В табл. 9 для нескольких малых значений максимальной тяги um приво-
дятся результаты расчета оптимального поворота плоскости орбиты КА на
малые углы для случая, когда α1 = 0,0, α2 = 1,0, β = 0,5, начальные условия
определяются соотношениями (5.1), а ориентация плоскости конечной орби-
ты параметрами (5.5):
(5.5)
ik = 6,9, Ωuk = 30,1.
Во всех случаях, представленных в табл. 9, первый этап оптимального
управления пассивный, второй этап активный u = um, третий этап пассив-
ный, четвертый этап активный u = -um и далее этапы повторяются. Из
табл. 9 видно, что при малых углах поворота плоскости орбиты и малых зна-
чениях максимальной тяги расход рабочего тела на оптимальный процесс и
конечное значение углового расстояния до перицентра слабо зависят от мак-
симального значения тяги.
104
На основе анализа большого количества расчетов можно сделать выводы о
некоторых закономерностях оптимального управления, которые, в частности,
иллюстрируются табл. 4.
Совокупность следующих друг за другом четырех этапов, которая начи-
нается с пассивного этапа и оканчивается активным этапом, укладывается
в одном витке и длительность таких совокупностей при переходе к следую-
щей аналогичной совокупности с той же последовательностью этапов умень-
шается. Уменьшение длительности таких совокупностей этапов ослабевает
с уменьшением максимальной тяги. В табл. 4 имеются две такие совокуп-
ности, следующие друг за другом: этапы 3-6 и этапы 7-10. В момент на-
чала третьего этапа истинная аномалия местонахождения КА была равна
162,4352, а в момент окончания шестого этапа истинная аномалия стала
равной 160,6597. Следовательно, за четыре этапа этой совокупности ис-
тинная аномалия нахождения КА изменилась на 358,2245, меньше 360, и,
следовательно, все четыре этапа совокупности 3-6 оказались внутри одного
витка. Началу совокупности 3-6 соответствовало τ = 2,375833, а окончанию
совокупности τ = 8,716633. Следовательно, длительность совокупности эта-
пов 3-6 равна 6,340800. Аналогично, на совокупности этапов 7-10 истинная
аномалия нахождения КА изменилась на 358,1784, а длительность совокуп-
ности этапов равна 6,339901, меньше, чем длительность предыдущей сово-
купности.
Совокупность четырех следующих друг за другом этапов, которая начина-
ется с активного этапа и оканчивается пассивным этапом, выходит за рамки
одного витка, и длительность таких совокупностей при переходе к следую-
щей аналогичной совокупности с той же последовательностью этапов увели-
чивается. Увеличение длительности таких совокупностей этапов ослабевает с
уменьшением максимальной тяги. В табл. 4 имеются две такие совокупности,
следующие друг за другом: 2-5 и 6-9. На совокупности этапов 2-5 истинная
аномалия местонахождения КА изменилась на 361,9029, больше 360, и, сле-
довательно, на совокупности этапов 2-5 КА проходит больше чем один виток.
Длительность нахождения КА на этой совокупности этапов равна 6,408289.
На совокупности этапов 6-9 истинная аномалия нахождения КА изменилась
на 361,9330, а длительность нахождения КА на этой совокупности этапов
равна 6,409072 и, следовательно, дольше, чем на предыдущей совокупности
этапов.
Длительности соответствующих активных этапов и разности между ис-
тинными аномалиями КА в конце и начале этапа для следующих друг за
другом совокупностей из четырех этапов уменьшаются при переходе к сле-
дующей совокупности. По мере уменьшения максимальной тяги этот эффект
ослабевает. В табл. 4 в совокупности этапов 2-5 второй этап активный. Дли-
тельность этого этапа равна 2,019758, а разность между истинными анома-
лиями равна 108,9035. В совокупности 6-9 соответствующим второму этапу
является активный шестой этап. Длительность шестого этапа равна 1,952269,
а разность между истинными аномалиями равна 105,2251. Видно, что дли-
тельность этапа и разность между истинными аномалиями для шестого эта-
па меньше, чем для второго этапа. Аналогичный вывод можно сделать из
105
сравнения характеристик соответствующих друг другу четвертого и восьмо-
го этапов.
Длительности соответствующих пассивных этапов и разности между ис-
тинными аномалиями КА в конце и начале этапа для следующих друг за
другом совокупностей из четырех этапов увеличиваются при переходе к сле-
дующей совокупности. По мере уменьшения максимальной тяги этот эф-
фект ослабевает. В табл. 4 в совокупности этапов 2-5 третий этап пассив-
ный. Длительность этого этапа равна 1,412733, а разность между истинными
аномалиями равна 66,8849. В совокупности 6-9 соответствующим третьему
этапу является активный седьмой этап. Длительность седьмого этапа рав-
на 1,486183, а разность между истинными аномалиями равна 70,4487. Вид-
но, что длительность этапа и разность между истинными аномалиями для
седьмого этапа больше, чем для третьего этапа. Аналогичный вывод можно
сделать из сравнения характеристик соответствующих друг другу пятого и
девятого этапов.
Отметим, что длительность первого этапа существенно зависит от началь-
ного значения истинной аномалии и по этой причине длительность первого
этапа не подчиняется перечисленным закономерностям.
6. Заключение
В статье получено с использованием кватернионного дифференциально-
го уравнения ориентации орбиты КА и принципа максимума Понтрягина в
нелинейной постановке новое решение задачи об оптимальном переводе орби-
ты КА с переменной массой на заданную плоскость. Управление движением
аппарата производится с помощью ограниченной по модулю реактивной тяги,
ортогональной к плоскости оскулирующей орбиты КА. Минимизируются вре-
мя и суммарный импульс тяги, затраченные на процесс управления. Расчеты
оптимального управления выполнены для случая минимизации комбиниро-
ванного функционала качества процесса управления и для случая миними-
зации суммарного импульса тяги. Получены примеры численного решения
задачи, содержащие до 192 активных и пассивных этапов управления (боль-
шое количество активных и пассивных этапов управления связано с мало-
стью величины тяги реактивного двигателя, характерной для современных
двигателей КА).
Исследовано влияние максимальной величины тяги, скорости истечения
рабочего тела из сопла двигателя, весовых множителей в функционале ка-
чества процесса управления на процесс перевода орбиты КА на заданную
плоскость: на длительность процесса управления, на число и длительность
этапов управления, а также на суммарный расход топлива. Установлено, что
угловое расстояние до перицентра орбиты в конце процесса управления слабо
зависит от перечисленных выше параметров. Длительность процесса управ-
ления, а также число этапов управления слабо зависят от скорости истечения
рабочего тела, но существенно зависят от максимального значения тяги.
Выявлены закономерности в поведении длительностей совокупностей из
четырех активных и пассивных этапов управления, следующих друг за дру-
гом в процессе оптимального управления ориентацией плоскости орбиты КА.
106
СПИСОК ЛИТЕРАТУРЫ
1.
Дубошин Г.Н. Небесная механика. Основные задачи и методы. M.: Наука, 1968.
2.
Абалакин В.К., Аксенов Е.П., Гребенников Е.А. и др. Справочное руководство
по небесной механике и астродинамике. М.: Наука, 1976.
3.
Копнин Ю.М. К задаче поворота плоскости орбиты спутника // Космич. иссле-
дования. 1965. Т. 3. Вып. 4. C. 22-30.
4.
Лебедев В.Н. Расчет движения космического аппарата с малой тягой. М.: ВЦ
АН СССР, 1968. 108 с.
5.
Борщевский М.З., Иослович М.В. К задаче о повороте плоскости орбиты спут-
ника при помощи реактивной тяги // Космич. исследования. 1969. Т. 7. Вып. 6.
C. 8-15.
6.
Гродзовский Г.Л., Иванов Ю.Н., Токарев В.В. Механика космического полета.
Проблемы оптимизации. М.: Наука, 1975.
7.
Охоцимский Д.Е., Сихарулидзе Ю.Г. Основы механики космического полета:
Уч. пос. М.: Наука, 1990.
8.
Ишков С.А., Романенко В.А. Формирование и коррекция высокоэллиптической
орбиты спутника земли с двигателем малой тяги // Космич. исследования. 1997.
Т. 35. № 3. С. 287-296.
9.
Иванов Н.М., Лысенко Л.Н. Баллистика и навигация космических аппаратов.
М.: Дрофа, 2004.
10.
Чернявский Г.М., Бартенев В.А., Малышев В.А. Управление орбитой стацио-
нарного спутника. М.: Машиностроение, 1984.
11.
Решетнёв М.Ф., Лебедев А.А., Бартенев В.А. и др. Управление и навигация ис-
кусственных спутников Земли на околокруговых орбитах. М.: Машиностроение,
1988.
12.
Bartenev V., Malyshev V., Rayevsky V., et al. Attitude and orbit control systems of
Russian communication, geodetic and navigation spacecraft // Space Technol. 1999.
V. 19. No. 3-4. P. 135-147.
13.
Testoyedov N., Rayevsky V., Somov Ye., et al. Attitude and orbit control systems
of Russian communication, navigation and geodesic satellites: History, Present and
Future // IFAC PapersOnLine. 2017. V. 50. No. 1. P. 6422-6427.
14.
Battin R.H. An Introduction to the Mathematics and Methods of Astrodynamics.
N.Y.: AIAA Press, 1987.
15.
Huntington G.T., Rao A.V. Optimal Reconfiguration of a Tetrahedral Formation via
a Gauss Pseudospectral Method // Proc. AAS/AIAA Astrodynam. Special. Conf.
AS 05-338. 2005. P. 1-22.
16.
Челноков Ю.Н. Применение кватернионов в теории орбитального движения ис-
кусственного спутника. Ч. 2 // Космич. исследования. 1993. T. 31. Вып. 3.
C. 3-15.
17.
Челноков Ю.Н. Кватернионные и бикватернионные модели и методы механики
твердого тела и их приложения. Геометрия и кинематика движения. М.: Физ-
матлит, 2006.
18.
Челноков Ю.Н. Кватернионные модели и методы динамики, навигации и управ-
ления движением. М.: Физматлит, 2011.
19.
Челноков Ю.Н. Оптимальная переориентация орбиты космического аппарата
посредством реактивной тяги, ортогональной плоскости орбиты // Прикладная
математика и механика. 2012. Т. 76. Вып. 6. С. 897-914.
107
20.
Челноков Ю.Н. Кватернионная регуляризация в небесной механике и астродина-
мике и управление траекторным движением. II // Космич. исследования. 2014.
T. 52. № 4. C. 322-336.
21.
Сергеев Д.А., Челноков Ю.Н. Оптимальное управление ориентацией орбиты
космического аппарата // Математика. Механика: Сб. науч. тр. Саратов: Изд-во
Сарат. ун-та, 2001. Вып. 3. С. 185-188.
22.
Панкратов И.А., Сапунков Я.Г., Челноков Ю.Н. Об одной задаче оптимальной
переориентации орбиты космического аппарата // Изв. Сарат. ун-та. Нов. сер.
Сер. Математика. Механика. Информатика. 2012. Т. 12. № 3. С. 87-95.
23.
Сапунков Я.Г., Челноков Ю.Н. Исследование задачи оптимальной переориента-
ции орбиты космического аппарата посредством ограниченной или импульсной
реактивной тяги, ортогональной плоскости орбиты. Ч. 1 // Мехатроника, авто-
матизация, управление. 2016. Т. 17. № 8. С. 567-575.
24.
Сапунков Я.Г., Челноков Ю.Н. Исследование задачи оптимальной переориента-
ции орбиты космического аппарата посредством ограниченной или импульсной
реактивной тяги, ортогональной плоскости орбиты. Ч. 2 // Мехатроника, авто-
матизация, управление. 2016. Т. 17. № 9. С. 633-643.
25.
Deprit A. Ideal frames for perturbed keplerian motions // Celestial Mechan. 1976.
V. 13. No. 2. P. 253-263.
26.
Брумберг В.А. Аналитические алгоритмы небесной механики. М.: Наука, 1980.
27.
Бранец В.Н., Шмыглевский И.П. Применение кватернионов в задачах ориента-
ции твердого тела. М.: Наука, 1973.
28.
Бранец В.Н., Шмыглевский И.П. Введение в теорию бесплатформенных инер-
циальных навигационных систем. М.: Наука, 1992.
29.
Челноков Ю.Н. Оптимальная переориентация орбиты космического аппарата
посредством реактивной тяги, ортогональной плоскости орбиты // Математика.
Механика: Сб. науч. тр. Саратов: Изд-во Сарат. ун-та, 2006. Вып. 8. С. 231-234.
30.
Панкратов И.А., Сапунков Я.Г., Челноков Ю.Н. Решение задачи оптималь-
ной переориентации орбиты космического аппарата с использованием кватер-
нионных уравнений ориентации орбитальной системы координат // Изв. Сарат.
ун-та. Нов. сер. Сер.: Математика. Механика. Информатика. 2013. Т. 13. № 1-1.
С. 84-92.
31.
Ильин В.А., Кузмак Г.Е. Оптимальные перелеты космических аппаратов с дви-
гателями большой тяги. M.: Наука, 1976.
32.
Моисеев Н.Н. Элементы теории оптимальных систем. М.: Наука, 1975.
Статья представлена к публикации членом редколлегии Е.Я. Рубиновичем.
Поступила в редакцию 25.04.2018
После доработки 05.03.2019
Принята к публикации 25.04.2019
108