Автоматика и телемеханика, № 8, 2019

Управление в технических системах

Ю.Н. ЧЕЛНОКОВ, д-р физ.-мат. наук (ChelnokovYuN@gmail.com)

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

ОПТИМАЛЬНЫЙ ПОВОРОТ ПЛОСКОСТИ ОРБИТЫ

КОСМИЧЕСКОГО АППАРАТА ПЕРЕМЕННОЙ МАССЫ

В ЦЕНТРАЛЬНОМ ГРАВИТАЦИОННОМ ПОЛЕ

ПОСРЕДСТВОМ ОРТОГОНАЛЬНОЙ ТЯГИ¹

С использованием кватернионов и принципа максимума решена в нели-

нейной постановке задача об оптимальном переводе орбиты космического

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

движением аппарата производится с помощью ограниченной по модулю

реактивной тяги, ортогональной к плоскости оскулирующей орбиты КА.

Учитывается изменение массы аппарата за счет расхода рабочего тела

на процесс управления. Функционал, определяющий качество процесса

управления, представляет собой линейную свертку с весовыми множите-

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

на процесс управления.

Излагается теория решения задачи. Приводятся результаты расчетов

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

нированном функционале качества процесса управления одновременно

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

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

содержащие до 192 пассивных и активных этапов. Установлены законо-

мерности оптимального управления поворотом плоскости орбиты КА.

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

орбиты, ограниченная (малая) реактивная тяга, оптимальное управление,

кватернион.

DOI: 10.1134/S0005231019080087

1. Введение

В работе с использованием кватернионного дифференциального уравне-

ния ориентации орбиты космического аппарата (КА) и принципа максимума

Понтрягина изучается задача оптимальной переориентации плоскости орби-

ты КА с переменной массой посредством ограниченной по модулю реактив-

ной тяги, ортогональной плоскости оскулирующей орбиты. Частным случаем

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

значение задача коррекции угловых элементов орбиты КА, когда изменения

¹ Исследование выполнено при финансовой поддержке Российского фонда фундамен-

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

угловых элементов плоскости орбиты в процессе управления имеют малые

значения. Использование управления, ортогонального плоскости оскулирую-

щей орбиты КА, позволяет корректировать элементы орбиты КА, сохраняя

форму и размеры орбиты КА неизменными. Это ценное свойство изучаемого

процесса переориентации орбиты КА является полезным как при решении за-

дачи коррекции угловых элементов орбиты КА, так и других задач механики

космического полета, например при управлении конфигурацией группировки

спутников.

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

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

процесс управления. Функционал, определяющий качество процесса управ-

ления, представляет собой линейную свертку с весовыми множителями двух

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

процесс управления. Для выбранного функционала оптимальное управление

состоит из пассивных этапов, на которых тяга отсутствует, и активных эта-

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

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

ниям ориентации орбиты КА и изучаемой проблеме оптимальной переориен-

тации плоскости орбиты и орбиты КА посредством реактивной тяги, ортого-

нальной плоскости оскулирующей орбиты. Излагается кватернионная нели-

нейная теория решения задачи оптимальной переориентации плоскости орби-

ты КА переменной массы в непрерывной постановке с использованием огра-

ниченной (малой) тяги, ортогональной плоскости орбиты КА, для нефиксиро-

ванного числа включений реактивной тяги. В отличие от большинства других

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

ния центра масс КА, а вектор реактивной тяги двигателя КА. В состав ис-

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

цесс переориентации плоскости орбиты КА, входят кватернионное дифферен-

циальное уравнение ориентации оскулирующей орбиты КА, дифференциаль-

ное уравнение для истинной аномалии, характеризующей положение центра

масс КА на орбите, и дифференциальное уравнение, описывающее изменение

массы КА в процессе переориентации плоскости его орбиты.

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

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

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

рабочего тела на процесс управления. Показано, что с уменьшением допусти-

мого максимального значения величины тяги увеличивается продолжитель-

ность процесса управления и число этапов в процессе оптимального управле-

ния. Получены примеры оптимального управления, содержащие до 192 пас-

сивных и активных этапов. Анализ результатов расчетов позволил опреде-

лить закономерности в оптимальном управлении плоскостью орбиты КА.

Отметим, что оптимальное управление посредством ортогональной тяги

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

орбиты КА (коррекции наклонения и долготы восходящего узла) при сохра-

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

ми элементами начальной и конечной орбиты малы. Однако в статье также

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

для немалых разностей угловых элементов орбиты. Это, с одной стороны,

демонстрирует возможности разработанной программы численного решения

задачи (ее работоспособность для широкого диапазона изменения парамет-

ров, входящих в постановку задачи), а с другой стороны, позволяет, более

полно установить закономерности такого управления. Кроме того, в буду-

щем нельзя исключить необходимости решения задачи для случаев больших

изменений значений наклонения и долготы восходящего узла при сохранении

неизменными формы и размеров орбиты КА.

2. Задачи оптимальной переориентации плоскости орбиты КА

Будем считать, что вектор ускорения w центра масс КА от тяги реактив-

ного двигателя во все время управляемого движения КА направлен ортого-

нально плоскости оскулирующей орбиты, т.е. ортогонально радиус-вектору r

и вектору v скорости центра масс КА (коллинеарно вектору c = r × v мо-

мента скорости центра масс КА). Тогда дифференциальные уравнения дви-

жения центра масс КА в ньютоновском гравитационном поле, описывающие

изменение размеров и формы мгновенной орбиты КА, интегрируются, да-

вая уравнение конического сечения. Поэтому управляемое движение центра

масс КА в этом случае описывается дифференциальными уравнениями, ха-

рактеризующими изменение мгновенной ориентации орбиты КА или исполь-

зуемой (например, орбитальной) вращающейся системы координат, в которой

записываются исходные уравнения движения центра масс КА, и дифферен-

циальным уравнением для истинной аномалии, характеризующей положение

центра масс КА на орбите. В процессе такого управления оскулирующая ор-

бита КА поворачивается в пространстве, не изменяя своих размеров и формы.

Движение центра масс КА будем рассматривать в инерциальной системе

координат - геоцентрической экваториальной системе координат OX₁X₂X₃

(X) c началом в центре O притяжения Земли. Ось OX₃ этой системы ко-

ординат направлена вдоль оси суточного вращения Земли, оси OX₁ и OX₂

лежат в плоскости экватора Земли, ось OX₁ направлена в точку весенне-

го равноденствия для Земли, ось OX₂ дополняет систему до правой тройки

векторов.

Введем также в рассмотрение систему координат ξ, связанную с плоско-

стью и перицентром орбиты КА. Начало этой системы координат находится

в центре O (или в перицентре орбиты), ось ξ₁ направлена вдоль радиус-век-

тора перицентра орбиты, ось ξ₃ перпендикулярна плоскости орбиты и име-

ет направление постоянного по модулю вектора c момента скорости центра

масс КА относительно центра O, а ось ξ₂ образует правую тройку с осями ξ₁

и ξ₃. Ориентация системы координат ξ в инерциальной системе координат X

характеризует собой ориентацию орбиты КА в инерциальном пространстве и

традиционно задается тремя угловыми оскулирующими элементами орбиты

[1, 2]: долготой восходящего узла Ω_u, наклоном орбиты i и угловым расстоя-

нием перицентра от узла ω_π.

Дифференциальные уравнения, описывающие мгновенную ориентацию

орбиты КА в инерциальной системе координат в угловых элементах орбиты

в рассматриваемом случае ортогональности вектора реактивной тяги плос-

кости оскулирующей орбиты КА, имеют вид [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/r², r = p/(1 + e cos ϑ), c = const,

где ϑ - истинная аномалия (угловая переменная, отсчитываемая в плоско-

сти орбиты от ее перицентра и характеризующая положение КА на орбите),

r = |r| - модуль радиус-вектора центра масс КА, p и e - параметр и эксцен-

триситет орбиты, c = |c| = |r × v| - постоянная площадей (модуль вектора

момента скорости v центра масс КА), w - проекция вектора реактивного

ускорения w на направление вектора момента скорости центра масс КА (ал-

гебраическая величина реактивного ускорения, перпендикулярного плоско-

сти оскулирующей орбиты КА), t - время.

Задача оптимального поворота плоскости орбиты КА в угловых перемен-

ных формулируется следующим образом: требуется построить управление w,

переводящее орбиту, изменение ориентации которой описывается уравнения-

ми (2.1), из заданного начального положения

Ω_u = Ω_u(t₀) = Ω0u, i = i(t₀) = i⁰, ω_π = ω_π(t₀) = ω0π, i⁰ = 0,π

в требуемое конечное положение

Ω_u = Ω_u(t₁) = Ω^∗u, i = i(t₁) = i^∗, i^∗ = 0,π

при свободном значении углового расстояния перицентра от узла ω_π на пра-

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

ционал качества процесса поворота плоскости орбиты КА.

Частные случаи задачи переориентации плоскости орбиты КА рассматри-

вались в [3-7]. В [8] рассматривается задача коррекции угловых элементов

орбиты Ω_u, i и ω_π с помощью реактивного ускорения, ортогонального плос-

кости оскулирующей орбиты КА, называемого в этой работе “бинормальным

реактивным ускорением”. Задача решается с помощью принципа максимума

и усреднения уравнений. Из усредненных уравнений получен ряд аналити-

ческих соотношений для определения затрат характеристической скорости

в частных случаях коррекции одного или двух элементов орбиты (наклона

орбиты, долготы восходящего узла) при условии малости изменения наклона

орбиты и долготы восходящего узла. По словам авторов [8] уравнения зада-

чи оптимизации в полном объеме не приведены и не проанализированы из-за

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

Замечание 1. Модель (2.1), которая описывает ориентацию оскулирую-

щей орбиты КА при ортогональности вектора ускорения к ее плоскости и име-

ет особенность для наклонения i при значениях i = 0 и i = π, является частью

системы дифференциальных уравнений для всех шести кеплеровых оску-

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

установки [9, стр. 97]. При этом, наряду с указанными особыми значениями

наклонения, здесь имеется особое значение эксцентриситета e = 1, когда ор-

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

углового расстояния перицентра от узла ω_π. Решение этой проблемы получе-

но в начале 1970-х гг. с применением двух векторов в пространстве состояний

орбитального движения КА - вектора эксцентриситета e = {e cos ϖ, e sin ϖ},

где {·} ≡ col (·), ϖ = Ω_u + ω_π, и вектора наклонения i = {i cos Ω_u, i sin Ω_u}

[10, 11]. Этот метод, впервые реализованный в системе навигации и управ-

ления движением советского спутника связи Радуга на геостационарной ор-

бите в 1974 г. [12], в настоящее время модифицирован и успешно применяется

на борту российских спутников связи, навигации и геодезии [13]. Проблема

вырожденности классических орбитальных элементов частично решается в

механике космического полета за счет использования так называемых «невы-

рожденных» орбитальных элементов с элементами Лагранжа P₁ = e sin ϖ,

P₂ = ecos ϖ и Q₁ = tg(i/2)sin Ω_u, Q₂ = tg(i/2)cos Ω_u (иногда для них исполь-

зуют термин “equinoctial elements”) и соответствующих уравнений ориентации

орбиты R. Battin [14, 15]. В этих уравнениях сохраняется лишь особое значе-

ние наклонения i = π.

Решение задач оптимальной переориентации плоскости орбиты и орбиты

космического аппарата посредством реактивной тяги, ортогональной плоско-

сти оскулирующей орбиты, с помощью уравнений (2.1) в классических угло-

вых элементах орбиты в строгой нелинейной постановке достаточно сложно

в силу нелинейности этих уравнений, наличия в них особых точек i = 0, π,

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

этому для решения задачи оптимальной переориентации плоскости орбиты

и орбиты КА вместо угловых элементов орбиты целесообразно использовать

[16-24] параметры Эйлера (Родрига-Гамильтона).

Дифференциальные уравнения ориентации орбиты КА в параметрах Эй-

лера имеют вид [16-20]

2dΛ₀/dt = -Ω₁Λ₁ - Ω₂Λ₂,

2dΛ₁/dt = Ω₁Λ₀ - Ω₂Λ₃,

(2.2)

2dΛ₂/dt = Ω₂Λ₀ + Ω₁Λ₃,

2dΛ₃/dt = Ω₂Λ₁ - Ω₁Λ₂;

(2.3)

dϑ/dt = c/r²

r = p/(1 + ecosϑ), c = const,

Ω₁ = (r/c)w cos ϑ, Ω₂ = (r/c)w sin ϑ,

где Λ_j (j = 0, 3) - параметры Эйлера, характеризующие ориентацию ор-

биты КА (системы координат ξ) в инерциальной системе координат X;

Ω₁,Ω₂,Ω₃ = 0 - проекции вектора Ω мгновенной абсолютной угловой ско-

рости орбиты на связанные с ней координатные оси Oξ_i.

Отметим, что уравнения ориентации орбиты в параметрах Эйлера вида

(2.2) использовались для описания орбитального движения и другими авто-

рами [25, 26].

Параметры Эйлера Λ_j связаны с угловыми элементами орбиты соотноше-

ниями [17]

Λ₀ = cos(i/2)cos((Ω_u + ω_π)/2), Λ₁ = sin(i/2)cos((Ω_u - ω_π)/2),

(2.4)

Λ₂ = sin(i/2)sin((Ω_u - ω_π)/2), Λ₃ = cos(i/2)sin((Ω_u + ω_π)/2).

Соотношения (2.4) получены с использованием кватернионных формул

сложения конечных поворотов [17, 27, 28] (на основе умножения кватернио-

нов, соответствующих трем элементарным поворотам на углы Ω_u, i и ω_π в

указанной последовательности).

Угол наклона орбиты i определяется через параметры Эйлера по формуле

√

Λ²¹ + Λ²²

i = 2arctg

Λ²⁰ + Λ²

Для определения угловых элементов Ω_u и ω_π сначала вычисляются вели-

чиныΩ и ω по формулам

Λ₁Λ₃ - Λ₀Λ₂

Ω= arctgΛ0Λ² +Λ1Λ³,

ω = arctg

Λ₀Λ₁ - Λ₂Λ₃

Λ₀Λ₁ + Λ₂Λ₃

а затем вычисляются Ω_u и ω_π по следующему однозначному алгоритму:

а) еслиΩ ≥ 0 и Λ₀Λ₂ + Λ₁Λ₃ ≥ 0, то Ω_u =Ω, а если Λ₀Λ₂ + Λ₁Λ₃ < 0, то

Ω_u =Ω + π;

б) еслиΩ < 0 и Λ₀Λ₁ - Λ₂Λ₃ ≤ 0, то Ω_u =Ω + π, а если Λ₀Λ₁ - Λ₂Λ₃ > 0,

то Ω_u =Ω + 2π;

в) если ω ≥ 0 и Λ₁Λ₃ - Λ₀Λ₂ ≥ 0, то ω_π = ω, а если Λ₁Λ₃ - Λ₀Λ₂ < 0, то

ω_π = ω + π;

г) если ω < 0 и Λ₀Λ₁ + Λ₂Λ₃ ≤ 0, то ω_π = ω + π, а если Λ₀Λ₁ + Λ₂Λ₃ > 0,

то ω_π = ω + 2π.

Эти формулы получены из соотношений (2.4) с учетом диапазонов изме-

нения углов i, Ω_u и ω_π: 0 ≤ i ≤ π, 0 ≤ Ω_u ≤ 2π, 0 ≤ ω_π ≤ 2π. Такие формулы

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

времени процесса управления через параметры Λ_j.

Используя вектор параметров Эйлера Λ = {Λ₀, Λ_v}, Λ_v = {Λ_i}, i = 1, 2, 3,

уравнения (2.2) можно представить в векторном виде

dΛ₀/dt = -(Ω · Λ_v)/2, dΛ_v/dt = (Λ₀Ω + Ω × Λ_v)/2,

Ω ≡ {Ω_i} = (r/c)w{cos ϑ,sinϑ,0}.

Уравнения (2.2) в кватернионной записи принимают компактный вид

[16-20] (здесь и далее применяются общепринятые обозначения операций с

кватернионами [17, 18, 27, 28])

(2.5)

2dΛ/dt = Λ ◦ Ω_ξ, Ω_ξ = Ω₁i₁ + Ω₂i₂ = (r/c)w(cos ϑ i₁ + sin ϑ i₂

где Λ = Λ₀ + Λ₁i₁ + Λ₂i₂ + Λ₃i₃ - кватернион ориентации орбиты КА (ква-

тернионный оскулирующий (медленно изменяющийся) элемент орбиты КА);

Ω_ξ - отображение вектора Ω на базис ξ (вектор Ω мгновенной абсолютной

угловой скорости орбиты направлен вдоль радиус-вектора r центра масс КА

и определяется формулой: Ω = (w/c)r); i₁, i₂, i₃ - векторные мнимые единицы

Гамильтона, ◦ - символ кватернионного умножения.

Отметим, что уравнения (2.2), (2.3) или (2.5), (2.3) - система пяти нели-

нейных стационарных дифференциальных уравнений первого порядка отно-

сительно параметров Эйлера Λ_j и истинной аномалии ϑ. Эти уравнения, в

отличие от четырех нелинейных дифференциальных уравнений (2.1) ориен-

тации орбиты в угловых элементах орбиты Ω_u, i, ω_π, не имеют особых точек

i = 0, π, к тому же при переходе в них от времени t к новой независимой пе-

ременной ϑ в соответствии с дифференциальным соотношением dϑ = (c/r²)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)) обладает свойством са-

мосопряженности: оно с точностью до обозначения кватернионной перемен-

ной совпадает с кватернионным сопряженным уравнением, что позволяет по-

низить размерность краевой задачи оптимизации на 4 единицы с использо-

ванием новой кватернионной переменной, являющейся мультипликативной

композицией кватернионных фазовой и сопряженной переменных (в виде их

произведения). Таким свойством классические дифференциальные уравне-

ния (2.1) ориентации орбиты в угловых элементах орбиты Ω_u, i, ω_π и уравне-

ния Battin не обладают, причем соответствующие им сопряженные уравнения

гораздо сложнее фазовых.

3. Кватернионная постановка задачи оптимальной

переориентации плоскости орбиты КА

Изменение ориентации орбиты КА переменной массы с помощью реактив-

ной тяги, ортогональной плоскости оскулирующей орбиты КА, описывается

системой дифференциальных уравнений

2dΛ/dt = Λ ◦ Ω_ξ, Ω_ξ = (r/c)(u^∗/m^∗)(cos ϑ i₁ + sin ϑ i₂),

(3.1)

dϑ

dm^∗

= -|u^∗|β^∗, r =

r²

1 + ecosϑ

где m^∗ - масса аппарата, u^∗ - проекция вектора реактивной тяги u^∗ на на-

правление вектора момента скорости центра масс КА (алгебраическая ве-

личина реактивной тяги, перпендикулярной плоскости оскулирующей орби-

ты КА), β^∗ - коэффициент пропорциональности, равный обратной величине

скорости истечения рабочего тела из сопла реактивного двигателя [31].

На величину тяги u^∗ наложено ограничение

(3.2)

|u^∗| ≤ u^∗m.

Компоненты Λ_j (j = 0, 3) кватерниона ориентации орбиты Λ выражаются

через угловые элементы орбиты Ω_u, i и ω_π по формулам (2.4).

Перейдем к безразмерным переменным ρ, τ, m, u, u_m, β по формулам

c²m∗0

(3.3)

r = pρ, t =

τ, ,m∗ = m∗0m, u∗ =

u, u^∗m =

u_m, β^∗ =

β,

p³

где m∗0 - начальная масса аппарата.

Уравнения движения КА (3.1) в безразмерных переменных примут вид

dΛ

Λ ◦ (i1 cosϑ + i2 sinϑ),

dτ

2(1 + ecos ϑ)m

(3.4)

dϑ

= (1 + e cos ϑ)² ,dm

= -β |u| .

dτ

На управляющий параметр - безразмерную тягу u наложено ограничение

(3.5)

|u| ≤ u_m.

В начальный момент времени состояние управляемой системы (3.4) опре-

деляется соотношениями

(3.6)

ϑ=ϑ_n, Λ=Λ_n

m = 1,

где компоненты кватерниона ориентации орбиты Λ_n выражаются однознач-

но через угловые элементы начальной орбиты КА i_n , Ω_un, ω_πn по форму-

лам (2.4).

Требуется перевести КА на новую орбиту, плоскость которой определяется

угловыми элементами i_k, Ω_uk.

Таким образом, при оптимизации поворота плоскости орбиты краевое на-

чальное условие Λ_n определяется однозначно согласно (2.4), но при назначе-

нии правого краевого условия Λ_k значение углового расстоянием перицентра

от узла ω_π является свободным.

В фазовом пространстве ϑ × Λ × m многообразие, на которое требуется

перевести управляемую систему (3.4), определяется соотношениями

i_k

Λ₀ sinΩ_uk sin

- Λ2 cos

- Λ3 cosΩuk sin

= 0,

(3.7)

i_k

Λ₁ sin Ω_uk cos

- Λ3 sin

- Λ2 cosΩuk cos

= 0,

получающимися после подстановки угловых элементов орбиты i_k, Ω_uk в (2.4)

и последующего исключения из них незаданного ω_πk. Для этого необходимо

из соотношений, определяющих Λ₂ и Λ₃, выразить sin(ω_πk/2) и cos(ω_πk/2) и

подставить их в соотношения для Λ₀ и Λ₁.

Качество процесса перевода КА на новую орбиту будем определять мини-

мальным значением функционала

τ_k

∫

(3.8)

J = (α₁ + α2 |u|)dτ, α₁ ≥ 0, α₂

≥ 0,

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

жителями α₁ и α₂ двух критериев: времени и импульса величины тяги, за-

траченных на процесс управления. Изменяя величины весовых множителей

в функционале (3.8), можно усиливать влияние того или иного критерия на

процесс управления.

Таким образом, требуется найти оптимальное управление безразмерной

реактивной тягой u = u (τ), удовлетворяющее ограничению (3.5), которое пе-

реводит управляемую систему (3.4) из начального состояния (3.6) на много-

образие (3.7) и сообщает функционалу (3.8) минимальное значение.

Промежуток времени процесса управления заранее не задается. Ниже рас-

сматривается случай, когда α₂ > 0, чтобы учесть влияние импульса величины

тяги на качество процесса управления.

4. Решение задачи с помощью принципа максимума Понтрягина

Функция Гамильтона для поставленной задачи оптимального управления

имеет вид

H = -(α₁ + α2 |u|) + χ(1 + ecosϑ)² +

(4.1)

(M, Λ ◦ (i₁ cos ϑ + i₂ sin ϑ)) - ηβ |u| .

2(1 + e cos ϑ)m

Здесь и далее для кватернионов a = a₀ + a_v и b = b₀ + b_v использует-

ся аналог их скалярного произведения (a, b) = scal((a₀ - a_v) ◦ (b₀ + b_v)) =

= a₀b₀ + a₁b₁ + a₂b₂ + a₃b₃, где a_j и b_j - компоненты кватернионов a и b (a₀,b₀

и a_v,b_v - скалярная и векторная части кватернионов a и b).

Сопряженные переменные χ, M = M₀ + M₁i₁ + M₂i₂ + M₃i₃, η, соответ-

ствующие фазовым переменным ϑ, Λ, m, удовлетворяют системе уравнений

dχ

= 2eχ(1 + e cos ϑ) sin ϑ +

dτ

(4.2)

(M, Λ ◦ (i₁ sin ϑ - i₂(e + cos ϑ))) ,

2(1 + e cos ϑ)²m

(4.3)

M ◦ (i1 cosϑ + i₂

sin ϑ),

dτ

2(1 + e cos ϑ)m

dη

(4.4)

(M, Λ ◦ (i₁ cos ϑ + i₂

sin ϑ)) .

dτ

2(1 + e cos ϑ)m²

Введем новую кватернионную переменную

(4.5)

N=N₀ +N₁i₁ +N₂i₂ +N₃i₃ ≡N₀ +N_v =^Λ

◦ M,

компоненты которой вычисляются по явным соотношениям (^Λ - сопряжен-

ный кватернион:^Λ = Λ₀ - Λ₁i₁ - Λ₂i₂ - Λ₃i₃).

Функция H и уравнения (4.2), (4.4), которым удовлетворяют сопряженные

переменные χ и η, примут вид

H = -(α₁ + α2 |u|) + χ(1 + ecosϑ)² +

(4.6)

(N₁ cos ϑ + N₂ sin ϑ) - ηβ |u| ,

2(1 + e cos ϑ)m

dχ

= 2eχ(1 + e cos ϑ) sin ϑ +

dτ

(4.7)

(N₁ sin ϑ - N₂(e + cos ϑ)) ,

2(1 + e cos ϑ)²m

dη

(4.8)

(N₁ cos ϑ + N₂

sin ϑ) .

dτ

2(1 + e cos ϑ)m²

Из условия максимума для функции H следует, что оптимальное управ-

ление определяется через фазовые и сопряженные переменные по формулам

^{ u_msign(ν_u), если ε_u ≥ 0;

(4.9)

u = 0,

если ε_u < 0,

где ν_u = N₁ cos ϑ + N₂ sin ϑ и ε_u = |ν_u| / [2m(1 + e cos ϑ)] - α₂ - ηβ.

Особый режим управления в настоящей работе не рассматривается.

Из первого уравнения системы (3.4) и из соотношений (4.9) следует, что

на активном этапе управления орбита КА вращается с безразмерным векто-

ром угловой скорости Ω_ξd = u_msign(ν_u) (i₁ cos ϑ + i₂ sin ϑ) /[m (1 + e cos ϑ)], а

на пассивном этапе ориентация орбиты не изменяется.

Так как правый конец траектории в фазовом пространстве ϑ × Λ × m яв-

ляется подвижным, то в конце движения, т.е. в момент τ = τ_k, должны выпол-

няться условия трансверсальности. На истинную аномалию ϑ и на массу m

в конце движения не налагается никаких условий, поэтому соответствующие

им сопряженные переменные χ и η при τ = τ_k должны обращаться в нули,

т.е.

(4.10)

χ(τ_k) = 0, η (τ_k

) = 0.

Из соотношений (3.7) с помощью метода неопределенных множителей

Лагранжа [32] выводятся условия трансверсальности для сопряженной ква-

тернионной переменной M

i_k

M₀ = λ₁ sinΩ_uk sin

i_k

M₁ = λ₂ sinΩ_uk cos

(4.11)

i_k

M₂ = -λ₁ cos

- λ2 cosΩuk cos

i_k

M₃ = -λ₁ cos Ω_uk sin

- λ2 sin

Если из первых двух соотношений формул (4.11) выразить неопределен-

ные множители λ₁ и λ₂ через M₀ и M₁ и подставить их в третье и четвертое

соотношения этих формул, то в результате получатся условия трансверсаль-

ности, свободные от множителей Лагранжа:

i_k

M₂ sin Ω_uk sin

+ M0 cos

+ M1 cosΩuk sin

= 0,

(4.12)

i_k

M₃ sin Ω_uk cos

+ M1 sin

+ M0 cosΩuk cos

= 0.

Из соотношений (3.7) и (4.12) можно получить, что при τ = τ_k

(4.13)

(M, Λ) = 0.

Для этого из первого соотношения системы (3.7) находится cos(i_k/2) и

подставляется в первое соотношение системы (4.12). Аналогично, из второго

соотношения системы (3.7) находится sin(i_k/2) и подставляется во второе со-

отношение системы (4.12). Соотношения, полученные из системы (4.12), при-

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

умножается на cos(i_k/2), а второе умножается на sin(i_k/2), левые и правые

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

торые входят произведения фазовых и сопряженных переменных с разными

индексами, уничтожаются, а произведения с одинаковыми индексами скла-

дываются. В итоге получается соотношение (4.13).

Условием (4.13) можно заменить одно из условий (4.12). Для уравнений

системы (3.4) и (4.3) соотношение (M, Λ) = c = const является первым ин-

тегралом, причем согласно (4.13) значение c = 0. По этой причине соотно-

шение (4.13) имеет место во все время процесса управления. Тогда согласно

(4.5) и (4.13) компонента N₀ ≡ 0.

Для компонент N₁, N₂, N₃ кватерниона N можно получить на основе (3.4),

(4.3), (4.5) систему дифференциальных уравнений

dN₁

u sin ϑ

dN₂

u cos ϑ

N₃,

dτ

(1 + e cos ϑ)m

dτ

(1 + e cos ϑ)m

(4.14)

dN₃

(N₁ sin ϑ - N₂ cos ϑ) .

dτ

(1 + e cos ϑ)m

Сопряженную кватернионную переменную M можно исключить из рас-

смотрения, если вместо уравнения (4.3) в постановке краевой задачи опти-

мального управления использовать систему уравнений (4.14). Так как ска-

лярная составляющая N₀ кватерниона N равна нулю, то можно ввести век-

тор N_v с координатами N₁, N₂, N₃, который представляет собой векторную

часть кватерниона N. Тогда систему дифференциальных уравнений (4.14)

можно записать в виде векторного дифференциального уравнения

(4.15)

dN_v/dτ = (u/m)(N_v × (cos ϑ i₁ + sin ϑ i₂

))/(1 + e cos ϑ).

С помощью (4.5) и с учетом того, что N₀ = 0, второе условие трансвер-

сальности из системы (4.12) можно переписать в виде

(Λ₀N₃ + Λ₁N₂ - Λ₃N₁) sinΩ_uk cos(i_k/2) -

(4.16)

-(Λ₁N₁ + Λ₂N₂ + Λ₃N₃)cos Ω_uk cos(i_k/2) +

+ (Λ₀N₁ + Λ₂N₃ - Λ₃N₂) sin(i_k/2) = 0.

Так как время момента окончания процесса не задается заранее, то в ко-

нечный момент времени функция Гамильтона (4.6) равна нулю:

(4.17)

H_opt(τ_k

) = 0.

Отметим, что H_opt(τ) = 0 является первым интегралом изучаемой задачи.

Таким образом, принцип максимума Понтрягина сводит поставленную за-

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

ференциальных уравнений (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 = u_m или u = -u_m (на активных этапах происходит изменение ориен-

тации орбиты). На пассивных этапах u = 0 (на пассивных этапах ориентация

орбиты не изменяется и совпадает с ориентацией орбиты, соответствующей

ориентации орбиты конца предыдущего этапа).

Численное решение краевой задачи для системы дифференциальных урав-

нений связано с большими трудностями, так как приходится подбирать недо-

стающие начальные условия, а также промежуток времени процесса управ-

ления. На сходимость итерационного процесса при решении краевой задачи

оптимизации влияет выбор начального приближения. Дополнительные труд-

ности со сходимостью процесса численного решения возникают при наличии

разрывов управляющих параметров. В рассматриваемой задаче при малых

значениях максимальной тяги число этапов, а следовательно, и число раз-

рывов управления возрастает. По этой причине необходимо, чтобы алгоритм

решения задачи учитывал наличие разрывов управления.

Таблица 1

u_m

τ_k

ϑ_k

m_k

ω_πk

n_et

0,25

1,720485

130,8325^◦

0,954649

0,781403

64,9135^◦

0,125

2,301065

158,9216^◦

0,947368

0,970003

65,2837^◦

0,075

4,786863

283,8775^◦

0,940712

1,640374

64,8113^◦

0,05

7,277088

86,7048^◦

0,935770

2,301095

64,7793^◦

При решении краевой задачи методом Ньютона или методом градиентно-

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

производных по начальным условиям сопряженных переменных и по проме-

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

вании системы дифференциальных уравнений методом Рунге - Кутты, если

моменты разрыва управляющего параметра не будут оказываться внутри ша-

га интегрирования. По этой причине на каждом шаге интегрирования надо

проверять, был ли разрыв управления внутри шага. Для этого вычисляется

согласно (4.9) в левой и правой точке шага выражение

ε_u = |ν_u|/[2m(1 + ecos ϑ)] - α₂ - ηβ.

Если знаки выражения одинаковые, то разрыв управления внутри шага

интегрирования отсутствовал. Если же знаки разные, то разрыв был. На

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

с более мелким шагом. Мелкий шаг, в котором окажется разрыв, делится на

два пропорциональных шага интегрирования, разделенных точкой разрыва

управления.

В примерах, результаты расчетов для которых приводятся в таблицах 1-5,

начальное состояние системы (3.4) определяется параметрами

(5.1)

e = 0,1, ϑ_n = 30,0^◦, i_n = 7,0^◦, Ω_un = 30,0^◦, ω_πn = 50,0^◦, m_n

= 1,0.

Ориентация плоскости конечной орбиты определяется углами

(5.2)

i_k = 20,0^◦, Ω_uk = 15,0^◦.

Примеры различаются максимальными значениями тяги u_m, значениями

коэффициента β, который определяет отношение характерного значения ско-

рости движения КА по орбите к скорости истечения рабочего тела из сопел

двигателя КА, и значениями весовых множителей α₁, α₂ в функционале (3.8).

Результаты расчетов представлены в безразмерных переменных, углы в гра-

дусах.

В табл. 1 для значений α₁ = 0,25, α₂ = 1,5, β = 0,2 и различных значе-

ний u_m приведены значения длительности процесса управления (вторая ко-

лонка), значения истинной аномалии КА и массы КА в конце процесса (тре-

тья и четвертая колонки соответственно), значения функционала качества

процесса (пятая колонка), значения углового расстояния до перицентра ор-

биты в конце процесса управления (шестая колонка) и число этапов n_et

в процессе управления (седьмая колонка). Из расчетов, представленных в

100

Таблица 2

τ_k

ϑ_k

m_k

ω_πk

n_et

0,2

5,491449

329,0163^◦

0,949182

0,254090

64,8251^◦

0,5

5,445041

325,8902^◦

0,879008

0,241984

64,8136^◦

1,0

5,377398

321,3655^◦

0,775613

0,224387

64,7971^◦

Таблица 3

τ_k

ϑ_k

m_k

ω_πk

n_et

0,2

8,586485

154,5172^◦

0,949954

0,250230

64,8388^◦

0,5

8,535597

152,0992^◦

0,880727

0,238546

64,8316^◦

1,0

8,461673

148,5677^◦

0,778526

0,221474

64,8215^◦

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

может создаваться двигательной установкой аппарата, увеличивается про-

должительность процесса управления по переводу орбиты КА на заданную

плоскость, незначительно увеличивается расход рабочего тела (топлива) и со-

ответственно уменьшается масса аппарата в конце процесса, увеличиваются

значения функционала качества процесса управления, увеличивается число

этапов процесса управления от двух при u_m, равном 0,25 или 0,125, до че-

тырех при u_m = 0,075 и до шести при u_m = 0,05. Для всех примеров первый

этап является пассивным. Важно отметить, что угловое расстояние до пе-

рицентра орбиты в конце процесса управления слабо зависит от величины

максимальной тяги u_m. Отсюда следует, что при всех значениях u_m аппарат

в результате процесса перевода его на заданную плоскость орбиты каждый

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

Проводились расчеты для других значений весовых множителей в функ-

ционале качества (3.8). Отметим, что если число этапов в процессе управле-

ния равно двум, то в этом случае оптимальный процесс управления не за-

висит от значений весовых множителей в функционале (3.8). В этом случае

от весовых множителей зависят лишь сопряженные переменные. Это связано

с тем, что в этом случае для определения двух моментов окончания этапов

(пассивного и активного) имеются два условия: (2.4) (при Ω_u = Ω_un, i = i_n

и ω_π = ω_πn) и (3.7). Благодаря этому краевую задачу для системы диффе-

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

нений для сопряженных переменных. По этой причине далее рассматрива-

ются решения задач оптимального управления для u_m = 0,075, u_m = 0,05 и

u_m = 0,025, для которых при заданных граничных условиях (5.1) и (5.2) чис-

ло этапов управления больше двух.

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

управления для α₁ = 0,0 и α₂ = 1,0 и нескольких значений коэффициента β.

В этой задаче требуется перевести КА из начального состояния (5.1) на орби-

ту, плоскость, которой определяется условиями (5.2), с минимальным значе-

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

тела). В табл. 2 приводятся результаты решения задачи для u_m = 0,075, а в

табл. 3 - для u_m = 0,05. В отличие от табл. 1 в табл. 2 и 3 в первой колонке

приводятся значения коэффициента β.

101

Таблица 4

№ u

i^◦

Ω^◦u

ω^◦π

ϑ^◦

0,0

0,356075

7,0

30,0

50,0

53,5317

1,0

-0,025

2,375833

9,5604

25,4348

54,5403

162,4352

0,949506

0,0

3,788566

9,5604

25,4348

54,5403

229,3201

0,949506

0,025

5,557657

11,6824

19,9718

59,9271

333,5040

0,905279

0,0

6,764364

11,6824

19,9718

59,9271

55,4346

0,905279

-0,025

8,716633

14,5198

18,5903

61,2971

160,6597

0,856472

0,0

10,202816

14,5198

18,5903

61,2971

231,1084

0,856472

0,025

11,903208

16,9037

15,4391

64,3495

331,2583

0,813962

0,0

13,173439

16,9037

15,4391

64,3495

57,37776

0,813962

-0,025

15,056534

20,0

15,0

64,7932

158,8377

0,766885

Таблица 5

u_m

τ_k

J = Δm_k

ϑ_k

n_ob n_et

ω_πk

0,075

5,377398

0,224387

321,3655^◦

64,7971^◦

0,05

8,461673

0,221474

148,5677^◦

64,8215^◦

0,025

15,056534

0,233115

158,8377^◦

64,7932^◦

0,0125

31,069513

0,235341

332,7422^◦

64,7170^◦

0,005

75,758654

0,239085

335,9647^◦

64,6988^◦

0,0025

152,301870

0,239065

335,9453^◦

64,6912^◦

0,00125

305,388734

0,239055

335,9358^◦

192

64,6873^◦

Из табл. 2 и 3 видно, что с увеличением максимальной тяги u_m и с умень-

шением коэффициента β повышается эффективность двигательной установ-

ки КА, так как при этом увеличивается конечная масса аппарата и, следова-

тельно, уменьшаются затраты рабочего тела на процесс управления. С увели-

чением коэффициента β уменьшается значение функционала качества про-

цесса. Наиболее значительно изменение коэффициента β влияет на поведение

массы аппарата. На время процесса управления, на угловое расстояние до пе-

рицентра в конце процесса изменение коэффициента β влияет незначительно.

Для тех же граничных условий (5.1), (5.2) и α₁ = 0,0, α₂ = 1,0, β = 1,0 при

условии, что максимальное значение тяги u_m = 0,025, число этапов в опти-

мальном управлении оказалось равным десяти. В табл. 4 приводятся резуль-

таты расчета для каждого этапа управления. В первой колонке указывается

номер этапа управления, во второй - управление на этапе, в третьей колон-

ке по восьмую - момент времени, угол наклона орбиты, угол восходящего

узла, угловое расстояние до перицентра, истинная аномалия и масса аппа-

рата в конце каждого этапа соответственно. Значение функционала качества

J = 0,233115.

В табл. 5 приводятся для различных значений максимальной тяги двига-

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

орбиты КА при тех же начальных условиях и начальной ориентации орби-

ты (5.1) и той же ориентации плоскости конечной орбиты (5.2), что и выше,

при условии, что α₁ = 0,0, α₂ = 1,0, β = 1,0. В первой колонке табл. 5 указано

102

Таблица 6

u_m

τ_k

J = Δm_k

ϑ_k

n_ob n_et

ω_πk

0,0025

12,246738

0,020574

354,8571^◦

50,9921^◦

0,001

31,382276

0,020571

354,8337^◦

50,9918^◦

0,0005

63,275304

0,020570

354,8270

50,9916^◦

Таблица 7

u_m

τ_k

J = Δm_k

ϑ_k

n_ob n_et

ω_πk

0,0025

11,577105

0,020536

356,9952^◦

49,0069^◦

0,001

30,707830

0,020500

356,6382^◦

49,0067^◦

0,0005

62,599522

0,020491

356,5389^◦

49,0064^◦

0,00025

126,385110

0,020486

356,4930^◦

49,0064^◦

Таблица 8

τ_k

m_k

ϑ_k

Δm_k

ω_π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

49,0064^◦

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

во второй колонке - продолжительность процесса управления, в третьей ко-

лонке - значение функционала качества процесса, который определяет расход

рабочего тела на процесс управления Δm_k, в четвертой колонке - значение

истинной аномалии, определяющей положение КА на орбите, в пятой и ше-

стой колонках указано число полных витков n_ob, выполненных КА в процессе

движения, и число этапов n_et в оптимальном управлении соответственно, в

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

ты в конце процесса управления.

Из табл. 5 видно, что наименьший расход рабочего тела наблюдается при

u_m = 0,05, угловое расстояние до перицентра орбиты в конце процесса слабо

зависит от значения максимальной тяги.

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

плоскости орбиты при условии, что плоскости начальной и конечной орбит

близки друг к другу. Различия между соответствующими угловыми элемен-

тами составляют один градус. Расчеты проводились при значениях α₁ = 0,0,

α₂ = 1,0, β = 1,0, e = 0,1. В начальный момент положение орбиты и поло-

жение КА на орбите определяется условиями (5.1). Ориентация плоскости

конечной орбиты определяется углами

(5.3)

i_k = 8,0^◦, Ω_uk = 29,0^◦.

В табл. 6 представлены результаты расчетов для трех значений макси-

мальной тяги. На первом этапе пассивное движение, на втором этапе u =

= -u_m, на третьем этапе пассивное движение, на четвертом этапе u = u_m и

далее этапы повторяются.

103

Таблица 9

u_m

τ_k

J = Δm_k/β

ϑ_k

n_ob n_et

ω_πk

0,00025

27,665243

0,001740

151,8025^◦

49,9008^◦

0,0001

59,646801

0,001788

156,0003^◦

49,9007^◦

0,00005

120,127999

0,001794

317,0792^◦

49,9007^◦

0,000025

241,419823

0,001792

317,0605^◦

148

49,9007^◦

В табл. 7 приводятся расчеты с другим начальным значением истинной

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

(5.4)

e = 0,1, ϑ_n = 75,0^◦, i_n = 7,0^◦, Ω_un = 30,0^◦, ω_πn = 50,0^◦, m_n

= 1,0,

а конечные условия (5.3) сохраняются. В этом случае первый этап активный

u = u_m, второй этап пассивный, третий этап активный u = -u_m, четвертый

этап пассивный и далее этапы повторяются. За счет увеличения начальной

аномалии, характеризующей положение КА на начальной орбите, первый

этап в отличие от приведенных ранее примеров оказался активным.

Из табл. 6 и 7 видно, что с уменьшением максимального значения тя-

ги существенно увеличивается продолжительность процесса и незначительно

уменьшается расход рабочего тела (в пятом знаке).

В табл. 8 для различных значений коэффициента β приводятся результа-

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

u_m = 0,00025, α₁ = 0,0, α₂ = 1,0, начальные условия определяются соотноше-

ниями (5.4), а ориентация плоскости конечной орбиты - соотношениями (5.3).

В этом случае КА совершает 19 полных витков, оптимальное управление со-

держит 39 этапов, первый этап активный u = u_m, второй этап пассивный,

третий этап активный u = -u_m, четвертый этап пассивный и далее этапы

повторяются, на последнем этапе u = -u_m.

Из табл. 8 видно, что от значения коэффициента β существенно зависит

только расход рабочего тела на процесс управления Δm_k, другие показатели

процесса меняются незначительно.

В табл. 9 для нескольких малых значений максимальной тяги u_m приво-

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

малые углы для случая, когда α₁ = 0,0, α₂ = 1,0, β = 0,5, начальные условия

определяются соотношениями (5.1), а ориентация плоскости конечной орби-

ты параметрами (5.5):

(5.5)

i_k = 6,9^◦, Ω_uk = 30,1^◦.

Во всех случаях, представленных в табл. 9, первый этап оптимального

управления пассивный, второй этап активный u = u_m, третий этап пассив-

ный, четвертый этап активный u = -u_m и далее этапы повторяются. Из

табл. 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

СПИСОК ЛИТЕРАТУРЫ

Дубошин Г.Н. Небесная механика. Основные задачи и методы. M.: Наука, 1968.

Абалакин В.К., Аксенов Е.П., Гребенников Е.А. и др. Справочное руководство

по небесной механике и астродинамике. М.: Наука, 1976.

Копнин Ю.М. К задаче поворота плоскости орбиты спутника // Космич. иссле-

дования. 1965. Т. 3. Вып. 4. C. 22-30.

Лебедев В.Н. Расчет движения космического аппарата с малой тягой. М.: ВЦ

АН СССР, 1968. 108 с.

Борщевский М.З., Иослович М.В. К задаче о повороте плоскости орбиты спут-

ника при помощи реактивной тяги // Космич. исследования. 1969. Т. 7. Вып. 6.

C. 8-15.

Гродзовский Г.Л., Иванов Ю.Н., Токарев В.В. Механика космического полета.

Проблемы оптимизации. М.: Наука, 1975.

Охоцимский Д.Е., Сихарулидзе Ю.Г. Основы механики космического полета:

Уч. пос. М.: Наука, 1990.

Ишков С.А., Романенко В.А. Формирование и коррекция высокоэллиптической

орбиты спутника земли с двигателем малой тяги // Космич. исследования. 1997.

Т. 35. № 3. С. 287-296.

Иванов Н.М., Лысенко Л.Н. Баллистика и навигация космических аппаратов.

М.: Дрофа, 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