РАСЧЕТ ПРИБЛИЖЕННО-ОПТИМАЛЬНЫХ ПЕРЕЛЕТОВ КОСМИЧЕСКОГО АППАРАТА С ДВИГАТЕЛЯМИ МАЛОЙ ТЯГИ С ВЫСОКОЭЛЛИПТИЧЕСКОЙ НА ГЕОСТАЦИОНАРНУЮ ОРБИТУ Текст научной статьи по специальности «Математика»

УДК 629.78.076.6

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

© 2019 г. Салмин в.в., петрухина К.в., Кветкин А.А.

Самарский национальный исследовательский университет имени академика С.П. Королёва (Самарский университет) Московское шоссе, 34, г. Самара, Российская Федерация, 443086, e-mail: ssau@ssau.ru

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

Ключевые слова: приближенно-оптимальный закон управления, теория локальной оптимизации, электроракетный двигатель, высокоэллиптическая орбита, геостационарная орбита, математическая модель управляемого движения, принцип максимума Понтрягина.

DOI 10.33950/spacetech-2308-7625-2019-4-94-108

calculation of suboptimal high-elliptical orbit to geostationary orbit transfers for spacecraft with low thrusters

Salmin v.v., petrukhina K.v., Kvetkin A.A.

Samara National Research University (Samara University) 34 Moskovskoe shosse, Samara, 443086, Russian Federation, e-mail: ssau@ssau.ru

At present the geostationary orbit is where communication satellites are preferably placed. Conventional orbital insertion profiles using chemical propulsion are insufficiently effective and require the use of heavy launch vehicles. Combining electric thrusters with chemical propulsion increases the mass of payload. On the other hand, orbital injection of a spacecraft using electrical propulsion brings up the problem of looking for an optimal control law. The paper discusses the transfer of a spacecraft with low-thrust electrical thruster from a high-elliptical orbit to geostationary orbit. It proposes a suboptimal control law for the thrust vector. It provides examples of simulations of the transfer using the control law for various initial conditions. Parameters of intermediate high-elliptical orbits were selected to minimize the time of transfer to the final orbit, and estimates were made of the effects of the residual

САЛМИН в.в. ПЕТРУХИНА К.В. КВЕТКИН А. А.

САЛМИН Вадим Викторович — доктор технических наук, профессор, директор научно-исследовательского института космического машиностроения Самарского университета, e-mail: sputnik@ssau.ru

SALMIN Vadim Viktorovich — Doctor of Science (Engineering), Professor, Director of the Research Institute of space engineering at Samara University, e-mail: sputnik@ssau.ru

ПЕТРУХИНА Ксения Вячеславовна — кандидат технических наук, доцент Самарского университета, e-mail: kseniya_10.05@mail.ru PETRUKHlNA Ksenia Vyacheslavovna — Candidate of Science (Engineering), Associate Professor at Samara University, e-mail: kseniya_10.05@mail.ru

КВЕТКИН Александр Александрович — аспирант Самарского университета, e-mail: emancipee163@yandex.ru

KVETKIN Aleksandr Aleksandrovich — Post-graduate at Samara University, e-mail: emancipee163@yandex.ru

atmospheric drag during flight in the vicinity of the perigee of the orbit. Considering the low level of error, simplicity and high computational speed the proposed method can be used for trajectory design calculations.

Key words: suboptimal control law, local optimization theory, electric thruster, high-elliptical orbit, geostationary orbit, math model of controlled motion, Pontryagin's maximum principle.

Введение

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

(ЭРД) на этапах довыведения на ГСО КА AEHF-1, «Экспресс-АМ5», SES-9 [1-3].

Среди множества баллистических схем перелетов КА на ГСО с использованием ЭРД можно выделить три основных:

1) с помощью ракеты-носителя (РН) формируется переходная эллиптическая орбита, при этом радиус апогея орбиты может быть равным радиусу ГСО, а также — больше или меньше его (рис. 1). Функцию довыведения КА выполняет собственная электроракетная двигательная установка (ЭРДУ);

2) формируется низкая круговая орбита с помощью РН, затем с помощью химического разгонного блока (РБ) реализуется переходная орбита. В дальнейшем довыведение осуществляется собственным ЭРД КА, аналогично варианту 1;

3) в третьем варианте РН формирует начальную орбиту, затем с помощью химического РБ формируется переходная орбита. На последнем этапе перелета задействуется многоразовый электроракетный транспортный модуль с большой энерговооруженностью [4], который совершает обратный перелет на начальную орбиту.

Рис. 1. Баллистическая схема перелета с высокоэллиптической на геостационарную орбиту: 1 — выдача импульса последней ступенью ракеты-носителя; 2 — траектория движения КА с электроракетной двигательной установкой (ЭРДУ); 3 — начало работы ЭРДУ; 4 — переходный эллипс; 5 — геостационарная орбита

В настоящее время в США проводятся запуски геостационарных КА среднего класса, использующих ЭРДУ для выполнения всех динамических операций, включая довыведение КА с геопереходной орбиты на ГСО на платформе Boeing-702SP [5]. В Европе находится на стадии разработки схожая платформа — Electra/ARTES-33 [6]. Аналогом зарубежных космических транспортных платформ среднего класса являются отечественные унифицированные спутниковые платформы «Экспресс-1000» [7] и «Экспресс-2000» [8]. При выведении данных спутников предполагается использование химического РБ.

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

• КА SES-8 [9] и Thaicom-6 [10] на РН Falcon 9 v1.1, перевод которых с переходной орбиты осуществлялся с помощью двигателей большой тяги;

• КА SES-14 [11] и AlYah-3 [12] на РН Ariane-5, довыведение которых осуществлялось с орбиты при помощи ЭРДУ.

Высота апогея переходной орбиты в первом случае составляла ~80-90 тыс. км, во втором — совпадала с высотой ГСО.

26 декабря 2014 г. с космодрома Байконур был выполнен пуск РН «Протон» с РБ «Бриз-М» и КА «Экспресс-АМ5». Спутник отделился от РБ на орбите со следующими параметрами: наклонение 0,21°; высота перигея 33 694,66 км; апогея — 37 782,33 км. Дальнейшее движение в точку стояния на ГСО спутник совершил за счет собственных ЭРД коррекции СПД-100 за 73 дня. Изначально разработчики планировали запустить «Экспресс-АМ5» ракетой «Протон» по «южной» трассе, обеспечивающей выведение орбитального блока на опорную орбиту наклонением 48°, но в 2009 г. использование данной трассы было запрещено Казахстаном после запуска DirectTV-12. «Экспресс-АМ5» пришлось перепроектировать на менее выгодную трассу с наклонением опорной орбиты 51,5°. Масса аппарата получилась 3 400 кг, а «Протон» обеспечивал выведение на ГСО лишь 3 250 кг. Это заставило разработчиков КА пойти на довыведение «Экспресса» двигателями коррекции [13].

Схему выведения с использованием круговой промежуточной орбиты можно условно назвать традиционной, поскольку рациональная схема выведения и закон управления движением КА с ЭРД в данном случае хорошо исследованы. Использование промежуточных эллиптических орбит в настоящее время исследовано меньше. Наличие дополнительного баллистического параметра (величина эксцентриситета переходной орбиты) дает расширенные возможности по оптимизации схемы выведения КА на рабочую орбиту. Исследования показывают, что использование эллиптической промежуточной орбиты для многих проектов выведения КА на ГСО эффективнее, чем использование круговых промежуточных орбит. Промежуточная эллиптическая орбита может быть сформирована третьей ступенью РН тяжелого класса («Протон», «Ангара-А5», Delta-IV, Falcon-9, Arian-5).

Выбору параметров эллиптической орбиты базирования многоразовых буксиров с ЯЭРДУ посвящена работа [14].

Оптимальное значение эксцентриситета промежуточной орбиты может быть достаточно велико, поэтому можно говорить о схеме выведения на рабочую орбиту с использованием высокоэллиптической промежуточной орбиты, в т. ч. и при Hл > HГСО (^ — высота апогея орбиты; HГСО — высота ГСО).

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

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

Эта система уравнений имеет особенности при е = 0 и г = 0. На практике наиболее распространенным приемом является переход к равноденственным элементам. В данном случае при моделировании движения КА с ЭРДУ обычная система дифференциальных уравнений в оскулирующих элементах является предпочтительной — так за счет понятных и простых уравнений движения с традиционными орбитальными элементами

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

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

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

¿А ¿в ¿1

¿ю

¿а

¿и ¿О

2р

(1 - в)2

— [ввтМ + (1 + всов^Т],

—

всов2^ + 2сов& + в втМ +-:-:- Т

1 + всовО

W,

8т&(2 + всовО)

совМ + —:-Т

1 + всовО

ввти^г 1 + всовО

■ W

вти

ц 8тг(1 + всовО)

W,

р2

(1 + всовО)2

р2

(1 + всовО)^

вти W

(1 + всовО) Гц р V Р

где р = А(1 - е2) — фокальный параметр; А — большая полуось; е — эксцентриситет; & = и - ш — истинная аномалия; и — аргумент широты; ш — угловое расстояние перицентра от узла; О — долгота восходящего узла; г — наклонение орбиты; Ь — время; S, Т, W — проекции возмущающего ускорения на направление радиус-вектора, на перпендикулярное к нему в плоскости орбиты

и на перпендикулярное к плоскости орбиты, соответственно; ц = /Ы — произведение гравитационной константы и массы притягивающего центра.

В качестве возмущающих факторов, помимо реактивного ускорения, создаваемого ЭРД, рассматриваются возмущения, вызываемые нецентральностью поля тяготения Земли, а также возмущения, создаваемые верхней атмосферой Земли

(наличие теневых участков на траектории выведения не учитывается).

В качестве основной возьмем вторую зональную гармонику для потенциала сил притяжения, который при этом имеет вид [15]:

U =

Эг3

(3sin2¿sm2M - 1),

где е = 2,634-1010 км5/с — константа, определяющая сжатие Земли; г — текущий радиус КА.

Составляющие возмущающего ускорения, обусловленного потенциалом второй зональной гармоники, определяются следующими соотношениями [15]:

е

5 = — (38ш2шп2и - 1);

г г4

s

Т =--sin2¿sin2M;

g ,4 '

s

W =--sin22¿sinM.

в r

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

При разложении вектора аэродинамического ускорения КА по осям орбитальной системы координат TSW можно не учитывать радиальную (Sf) и боковую (W) составляющие аэродинамического ускорения в силу их малости. Тогда выражение для трансверсальной составляющей Tf примет вид:

ц

T = ар -г- (1 + ecosS)2,

2р

Sf = 0, = 0,

где р — плотность атмосферы; с — баллистический коэффициент КА, который вычисляется в соответствии с выражением:

а =

m

где СХ — коэффициент лобового сопротивления; S — эффективная площадь миделя КА; т — масса КА.

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

S, T, W будут

S = S + S ;

r g'

T = T + T + T;

rgf

Ж = Ж + Ж ,

Г ё '

где Sr, Т, — составляющие реактивного ускорения КА.

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

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

dw ¿О йы

ваются, поэтому уравнения -, -, -

¿Т ¿Т ¿Т

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

Сформулируем задачу об оптимальном изменении основных элементов орбиты. В качестве этих элементов примем большую полуось, наклонение и эксцентриситет орбиты.

В качестве критерия оптимальности примем продолжительность перелета Т.

Введем две правые системы координат: орбитальную (ОптЬ) и связанную с КА (ОХУ2) (рис. 2). Вектор тяги P направлен вдоль оси ОХ.

Рис. 2. К определению углов ориентации вектора тяги КА с ЭРДУ: 1 — мгновенная плоскость орбиты (МПО); 2 — мгновенная плоскость местного горизонта (МПМГ)

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

Т = ЗасозЯсозу;

Sr = 5«sinA,cosy; W = 5«siny.

(1)

s

Здесь а — модуль полного реактивного ускорения (а = а0 /(1 - а0г/с)); 5 — функция включения/выключения двигателей (5 = {0, 1}); Я — угол между проекцией вектора тяги на мгновенную плоскость орбиты и трансверсалью T (Я е [-180°; 180°]); у — угол отклонения тяги от мгновенной плоскости орбиты (у е [-90°; 90°]) (рис. 2). Очевидно, что при постоянной работе двигателя (5 = 1) время перелета равно моторному времени Т .

Управлениями в указанной задаче являются углы ориентации вектора тяги Я, у.

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

¿1

ренциальных уравнений для -, содер-

жащей ео8и.

Запишем граничные условия:

г = А(д = А0 е(д = е0

г(Ь0) = г0

г = г

к

А(Ьк) = Ак

е(К) = ек г(гк) = V

методика, основанная на принципе максимума л.С. понтрягина

Решение вариационной задачи можно проводить в соответствии с принципом максимума Л.С. Понтрягина [16] при условии, что ЭРДУ работает без выключений (5 = 1), поскольку в этом случае обеспечивается минимум общей продолжительности перелета.

Введем сопряженную вектор-функцию ¥, составим гамильтониан и найдем его максимум по управляющим функциям Я и у:

¥ = (¥А, ¥ , ¥)т.

^ А' е' и

Запишем выражение для гамильтониана в следующем виде:

где

н =

А =

А ( 1 - е2 ) Ц

[А/ + ЛтТ + Л„Щ,

2А(1 + е)

-е8т&¥. + 8т&¥ ;

(1 - е) А е

2А(1 + е) АТ =-(1 + еео8&)¥. +

+

(1 - е) есо82& + 2сов& + е 1 + есовО

р>

сочи

Аш =-¥.

W 1 + есочО !

С учетом уравнений (1) выражение для гамильтониана принимает вид:

н = АИ-еИ х

х а[А^08Яе08у + АТ8тЯео8у + AWsinу].

В соответствии с принципом максимума, необходимое условие оптимальности имеет вид:

дН = /А(1 - е2) х

ЭХ

Ц

а[-А^тЯео8у + Атео8Яео8у] = 0;

дН

= / А(1 - е2)

ду V Ц

х a[-ASеo8Я8inу - АТ8тЯ8ту + AWеo8у] = 0. Откуда получаем стационарные точки (2):

А

<

81пЯ = ±

ео8Я = ±

А+А

А

А

(2)

81Пу = ±

ш

ш

ео8у = ±

№+А2 + А

К+А2 +

Проверим стационарные точки на экстремум. Для этого составим матрицу производных второго порядка:

д =

21

где х

11

д2н

ЭХ2

; х12 = х21 =

^12 ^22

д2н д2ы ; Х22

дЯду ду2

Ц

х

X

11

Условием экстремума будет отрицательная определенность матрицы производных второго порядка, т. е. одновременное выполнение неравенств

(Х11Х22 Х12Х21

) > 0

; Х11 < 0. (3)

Таким образом, условиям (3) удовлетворяют стационарные точки (4), в которых

достигается максимум гамильтониана. (

A

sinX =

cqs^ =

{aT+AS

(4)

A

siny

cosy :

W

)A+As+AW

iAT + A

№+A2 + AW

В общем виде сопряженная система уравнений записывается следующим образом:

У, =

У в =

дН дА

дн

де

= fA(A, в, S, u, X, у, а, У a, ^в, ^i);

= f (A, в, S, u, X, у, а, Уа, Ув,

У. =

дн

д1

= 0, У = const.

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

сопряженных переменных УА, Уе, У.

Отыскание начальных значений обычно проводится с помощью модифицированного метода Ньютона [17], при этом необходимо задавать начальное приближение, которое может быть найдено, исходя из физического смысла этих значений или основываясь на результатах решения упрощенных задач. Применяемый метод позволяет найти оптимальное решение, однако это связано с большими вычислительными трудностями и требует высокой квалификации исследователя.

Например, многочисленные результаты решения задач оптимизации межорбитальных перелетов с малой тягой получены В.Г. Петуховым с применением принципа максимума Л.С. Понтряги-на [18]. Отмечается, что на первый план в вычислительных схемах решения задач выходит проблема сходимости и устойчивости алгоритма решения краевой задачи и единственности решения.

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

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

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

I = ДхтаДх ^ min,

к к

(5)

где Ахк = [ДА, Де, Здесь:

т - ак

ДА = ■

А

Дв = в(0 - вк; i = i(t) -

'au 0 0 a = [a.] = \ 0 «22 0 >, 2a.. = 1, 1 '0 0 a 1 ;

~33 ,

весовые

где аА = «Ш ае = a22, ai = а33 коэффициенты (элементы диагональной

матрицы) по большой полуоси, эксцентриситету и наклонению, соответственно.

Локально-оптимальными управлениями [19] в дальнейшем будем называть такие управления t, х), которые минимизируют не функционал динамической задачи I (интегральный), а подынтегральное

¿I

выражение, т. е. производную - в каждый

&

момент времени.

Рассмотрим задачу о минимуме функционала:

т ¿I

I = — ¿г + 10 ^ шт. 0 ¿г 0

Пусть 10 > 0. Тогда потребуем выполнения условий:

dI

sign(—-—) = const; dt

dI

dt

^ max.

Очевидно, что полученное решение при этом является еще и решением исходной задачи о минимуме функционала I.

Решение при этом получается в виде конечных соотношений, не содержащих неопределенные величины (сопряженные переменные в принципе максимума Л.С. Понтрягина [16]).

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

В общем случае синтез локально-оптимальных управлений не гарантирует абсолютного оптимума в исходной постановке задач. Оценка близости локально-оптимального управления к оптимальному приведена в работе Н.Н. Моисеева [19].

Представим систему уравнений в оску-лирующих элементах, описывающую динамику движения КА с двигателем малой тяги, в общем виде:

Э = /(Э, au),

(6)

где а — малый параметр (реактивное ускорение, развиваемое двигателем малой тяги); и — управление, которым в данном случае являются углы Я(г), у(г).

В силу дифференцируемости функции / перепишем выражение (6) в виде:

Э = /(Э, 0) + aBu + 0(a2),

(7)

где B =

д/(Э, y) dy

матрица частных

y - 0

производных.

Отбросим в уравнении (7) малые второго порядка 0(а2) и рассмотрим уравнение Э = /(Э, 0). Пусть его полный интеграл имеет вид:

Э = g(t, C).

(8)

Уравнение (8) может быть рассмотрено в качестве формулы замены переменных. Переходя от переменной Э к С,

уравнение (7) можно с точностью до величины 0(a2) заменить:

/ \-1

C = а

дс

Bu.

(9)

Рассмотрим функционал: 1(Э(Т)) = 1^(Т, С(Т))) = 1(С(Т)). (10)

Задача сводится к поиску минимума выражения (10) при условии (9). Решение этой задачи можно искать в следующем виде:

С = С0 + а С1 + 0(а2), причем С1 удовлетворяет уравнению (11):

Ci =

дс

Bu,

(11)

а функционал имеет вид:

(

I*(C(T)) = I(Co) + a — l dC

С + 0(a2).

с = C

Запишем функцию Гамильтона:

H = ¥7

дс

Bu = -a

dl* dC

v dC j

-1

Bu

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

dl* dC

dC

-1

Bu

У У

(12)

Теперь найдем локально-оптимальное управление для системы (9). Управление будем выбирать из условия минимума производной

dI*(C(T) dt

'dl* —, С

dC

= a

f dI* dC

dC

1 > Bu

(13)

Видно, что выражения (12) и (13) совпадают с точностью до множителя, т. е. их минимизирует одна и та же функция и(г). Поскольку этот результат был получен для уравнений, в которых были отброшены числа порядка 0(а2), то можно сделать вывод: локально-оптимальные управления тем ближе к оптимальным, чем меньше величина параметра а.

Приближенно-оптимальные схемы управления вектором тяги ЭРД, основанные на комбинации управлений, максимизирующих производные элементов орбиты А, е, i, описаны в работе [20].

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

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

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

Заменим этот функционал локальным критерием, обеспечивающим максимальную скорость изменения I:

di 1 DA = 2atAA--+

dt

Л0 dt

de di + 2a„Ae--+ 2aAi-^ max.

(14)

dt

dt

Введем обозначения, учитывая правые

dA de di

части уравнений движения

A* = a1AA

2A(1 + e) Ao(1 - e)

dt dt dt

esinS + anAesin&;

2A(1 + e)

A* = a.DA-(1 + ecos») +

1 A0(1 - e)

ecos2O + 2cosO + e

+ a2De-;

1 + ecosO

cosm

AW = a. Di-.

W 3 1 + ecosO

Тогда функционал записывается в виде:

di I A( 1 - e2)

= 2 1 -(--L [A*S + A*T + A*W\.

dt

Найдем стационарные точки, для чего решим уравнения:

5 Г DI dX \ dt j

д ( diЛ dy

dt

= 0;

= 0.

поиска

максимума

Результатом ¿I ¿I

- = -(Я(0, у(^) по двум переменным

dt dt

являются аналитические выражения (15), (16) для углов ориентации вектора тяги Я и у, где у — угол отклонения тяги от мгновенной плоскости орбиты; Я — угол между проекцией вектора тяги на плоскость орбиты и трансверсалью.

sin^ =

cqs^ =

А *

A*

)¡(A*yT(Ajj2

А*

(15)

)J(A*y + (А**)2

А * Aw

sinw V (A t)2 + A)2 + (A w)2 i (A*)2 + (A*)

cosy :

(16)

Полученный закон управления у(0, Я(^ имеет достаточно простую структуру и позволяет провести расчет динамического маневра без процедуры решения краевой задачи.

Как следует из выражения (14), от значений весовых коэффициентов аА,

а , а. зависит скорость изменения боль-

е i ^

шой полуоси, эксцентриситета и наклонения орбиты.

За счет подбора значений весовых коэффициентов можно добиться одновременности выполнения конечных условий. В первом приближении можно принять а . = а = а..

А е i

Сравнение точного и приближенного методов расчетов на модельных примерах

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

решения задач. На первом этапе апробация метода проводилась на примере перелетов типа «круговая орбита-круговая орбита». В качестве эталона были взяты результаты, приведенные в монографии В.Н. Лебедева [17].

В результате, решение с применением локально-оптимальной программы управления дает несколько отличные значения времени перелета относительно результатов, полученных по методу принципа максимума. Разница составляет 1,2-1,6% (табл. 1).

Таблица 1

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

Способ решения Параметры перелета

гО = 2О ООО км гк = 23 35О км Di = 19,О22° a = О,ОО498 м/с2 гО = 5О ООО км гк = 58 375 км Di = 19,О22° a = О,ООО8О м/с2 гО = 8О ООО км гк = 93 4ОО км Di = 19,О22° a = О,ООО31 м/с2

Время перелета, сут

Приближенный 5,2416 2О,648 41,76О

Точный 5,158О 2О,389 41,264

Отклонение, % 1,62 1,27 1,2О

Примечание. г0 — радиус начальной орбиты; гк — радиус конечной орбиты; Аг — изменение наклонения; а — модуль полного реактивного ускорения.

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

На втором этапе сравнение проводилось на примере перелетов типа «эллиптическая орбита - круговая орбита». Источником для сравнения послужили результаты, приведенные в статье В.Г. Петухова [18]. Конечные значения эксцентриситета и наклонения для всех случаев приняты близкими к нулю, чтобы избежать особенностей при интегрировании уравнений.

Анализ результатов показал, что разница по времени перелета между решениями, полученными с помощью предлагаемого метода, и точными результатами довольно мала и составляет ~0,02-1,00% (табл. 2).

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

Таблица 2

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

Способ решения Исходные данные

r1 = 42 171 км гпО = 6 871 км «о = 75° ra0 = 42 378 км гп0 = 6 578 км i = 7° «0 ' г1 = 46 500 км r° = 6 642,9 км i = 7° «0 ' r° = 34 171 км r° = 6 595 км i0 = 63,17°

гк = 42 165 км гк = 42 378 км гк = 42 378 км гк = 42 160 км

т0 = 1 320 кг I = 1500 с уд P = 0,332 Н т0 = 2 000 кг I =2000с уд P = 0,350 Н т0 = 1 500 кг I = 1 994,06 с уд P = 0,200 Н т0 = 776 кг I = 1 500 с уд P = 0,166 Н

Время перелета, сут

Приближенный 171,7303 139,0683 178,1134 193,3796

Точный 170,117 139,0382 177,3602 191,406

Отклонение, % 0,58 0,02 0,42 1,03

Примечание. r0 — начальный радиус апогея; r0 — начальный радиус перигея; гк — радиус конечной орбиты; т0 — начальная масса КА; I — удельный импульс; Р — тяга.

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

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

Для автоматизации расчетов баллистических параметров перелета КА с ЭРДУ между круговыми и эллиптическими некомпланарными орбитами была разработана программа NEOS.

Для определения оптимальных параметров переходной орбиты с точки зрения минимального времени перелета проведем моделирование с перебором высоты апогея промежуточной орбиты на отрезке r = [40 ООО; 80 ООО] для i = 28° и ra = [7О ООО; 11О ООО] для i = 51,6°. Исходные данные для моделирования приведены в табл. 3. В расчетах примем массу спутника для всех переходных орбит постоянной (тО = const).

Результаты моделирования приведены на рис. 3.

На следующем этапе исследований планируется провести уточненный анализ с учетом зависимости выводимой массы КА от радиуса апогея переходной орбиты.

Таблица 3

Исходные данные для перелета

Наклонение орбиты i, ° Высота перигея Нп, км Стартовая масса КА m0, кг Масса спутника на ГСО ™КА, кг Масса рабочего тела mpт, кг Тяга двигателя P, мН Удельный импульс I , с уд' Диапазон перебора г , тыс. км

28 200 3 500 2 600 900 360 1 600 40...80

51,6 70...110

301

,300

299

гй 298

Г"

О Ч 297

а о

Й

И я 295

о

И 294

293

292

291

С

\

80

Рис

40 50 60 70

Высота апогея, тыс. км а)

3 ■ Результаты анализа оптимального радиуса апогея: а

70 75 80 85 90 95 100 105 Высота апогея, тыс. км б)

для наклонения i = 28°; б — для наклонения i = 51,6°

Результаты перебора высот апогея орбиты для двух наклонений выявили оптимальные значения. В случае с наклонением I = 28° оптимальная высота апогея составляет 60 000 км, для I = 51,6° это значение возрастает до 93 000 км. Столь большое различие этого параметра для разных наклонений связано с тем, что весомую часть энергетики перелета на ГСО составляет изменение наклонения орбиты.

В соответствии с правой частью

¿1

дифференциального уравнения —-— с воз-

т

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

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

Таблица 4

Результаты моделирования

i, ° Элемент орбиты Необходимое значение Результат

e 0 291,22 сут

А, км 42 164 км 291,61 сут

28 i 0° 291,72 сут

Г — 291,72 сут

Масса рабочего тела — 567,11 кг

e 0 335,15 сут

А, км 42 164 км 336,19 сут

51,6 0 335,56 сут

Г — 336,19 сут

Масса рабочего тела — 653,56 кг

Примечание. Обозначения приведены в тексте.

Можно заметить, что в результате перелета был истрачен не весь запас рабочего тела. В итоге для случая с наклонением I = 28° остаток составляет 332,89 кг, а для I = 51,6° он составляет 246,44 кг. Излишки рабочего тела могут быть потрачены на продление срока

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

Также в результате расчета двух перелетов были подобраны весовые коэффициенты по правилу одновременного достижения элементами орбиты своих конечных значений; они приведены в табл. 5.

Таблица 5

Значения весовых коэффициентов

Весовой коэффициент i, °

28 51,6

aA 0,082 0,073

ae 0,285 0,456

a. 0,633 0,471

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

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

Данный факт можно объяснить тем, что сила аэродинамического сопротивления ВА в нижней точке траектории (наибольшее сопротивление на первом витке) примерно равна силе тяги Р (*А = 320 мН; Р = 360 мН), а длительность и протяженность участка воздействия крайне малы. За счет этого высота перигея орбиты выходит за пределы действия атмосферы за 3-4 витка.

В результате моделирования были получены зависимости элементов орбиты и управляющих переменных от времени, представленные на рис. 4-9.

Моделирование показало, что производная функционала не меняет свой знак, а сам функционал монотонно убывает до нуля, минимизируя невязки по большой полуоси, эксцентриситету и наклонению орбиты (рис. 4). Таким образом, обосновано применение метода локальной оптимизации.

0,70 0,65 0,60 0,55 >4 0,50 | 0,45 | 0,40 | 0,35 0,30 0,25 0,20 0,15 0,10 0,05

\ ...................i...................;.............

\..............

\

\

\

\ .................\

Г......................................................

\ для i -51,6'

it !■■

для i ~

0 50 100 150 200 250 300 Время, сут

Рис. 4. Изменение функционала I во времени

о Ей >1

50 45 40 35 30 25 20 15 10 5

О

t.................

\

\............ \

\ ..........\

\

\

\ .дляi= 51,6°

. 28" — -

--•1Я < " —-— ~ "

50 100 150 200 250 300 Время, сут

Рис. 5. Зависимость наклонения орбиты от времени

w

65 60 55 50 45 40 35 30 25 20 15 10 0

Радиус a nor

x.

_____ ......................

_____........._..4p

50

250

100 150 200

Время, сут

Рис. 6. Зависимость изменения большой полуоси, радиуса апогея и перигея от времени для i0 = 28°

90

80 70

60

50

АО

30

20

10

Jim ^С'д,

Попыпая nnTivrtb

--

——-—

0 50 100 150 200 250 300 Время, сут

Рис. 7. Зависимость изменения большой полуоси, радиуса апогея и перигея от времени для = 51,6°

0,8 0,7

0,6

0,5

0,4

0,3

0,2

0,1

= 51. .............

ДЛЯ I

,„ О.ЙЕ............\ .

.......ДЛЯ7

О

50

100 150 200 250 300 Время, сут

Рис. 8. Зависимость эксцентриситета орбиты от времени

На рис. 9 представлен график зависимости амплитуды управляющего угла у от времени для случая перелета с i = 51,6°. В период 0...250 сут полета видно постепенное монотонное уменьшение амплитудных значений угла у, связанное со скоростью уменьшения наклонения. Начиная от 250 и до 330 сут, наклонение орбиты начинает сначала медленно, а с 300 сут — резко убывать, с чем связано возрастание амплитудных значений угла у. Спад величины амплитудных значений в период с 328 до 336 сут коррелируется с замедлением скорости уменьшения эксцентриситета орбиты.

Зависимость амплитуды управляющего угла у от времени для i = 51,6° приведена на рис. 9.

£

20 lil

Л

о _—— ' 1

10

20

-30

.■in

-50 -60 70

80

90 i:

Л 50 100 150 200 250 300

Время, сут

Рис. 9. Линии амплитудных значений управляющего угла для = 51,6°

Зависимость управляющего угла X для случая перелета с ¡0 = 51,6° носит колебательный характер, угол меняется в пределах -180...180°.

Выводы

Показана принципиальная возможность применения метода локальной оптимизации для расчета перелетов с высокоэллиптических орбит на геостационарную. При этом получены наиболее выгодные значения высоты апогея для случаев перелета с различными наклонениями (I = 28°; I = 51,6°), рассчитано время перелета и потребный запас рабочего тела ЭРДУ. Также проведена оценка величины влияния атмосферного сопротивления на движение спутника с учетом низких значений высоты перигея (Нп = 200 км), приведены графики, отображающие характер изменения орбитальных параметров и углов ориентации вектора тяги ЭРДУ.

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

Основные задачи, решаемые с использованием довыведения КА на ГСО с помощью ЭРДУ:

• обеспечение выведения на ГСО КА увеличенной массы с использованием существующих средств выведения

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

• обеспечение выведения на ГСО КА с использованием средств выведения среднего класса;

• снижение стоимости выведения на ГСО;

• обеспечение конкурентоспособности российских геостационарных КА (компенсация северного расположения космодромов);

• обеспечение резерва массы геостационарных КА для повышения их надежности и эффективности.

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

1. Полярный П. Спасение AEHF-1 // Новости космонавтики. 2011. Т. 21. № 12(347). С. 47.

2. Официальный сайт АО «ИСС» имени академика М. Ф. Решетнёва». Режим доступа: https://www.iss-reshetnev.ru/ projects/telecommunication (дата обращения 23.04.2019 г.).

3. Черный И. С пятой попытки // Новости космонавтики. 2016. Т. 26. № 05(400). С. 27-28.

4. Хамиц И.И., Филиппов И.М., Буры-лов Л. С., Тененбаум С.М., Перфильев А.В., Гусак Д.И. Концепция космической транс-портно-энергетической системы на основе солнечного межорбитального электроракетного буксира // Космическая техника и технологии. 2017. № 1(16). С. 32-40.

5. Gunter's Space Page. Hughes/Boeing HS-702. Режим доступа: https://space. skyrocket.de/doc_sat/hs-702.htm (дата обращения 23.04.2019 г.).

6. European Space Agency. Electra. Режим доступа: https://artes.esa.int/news/ electra (дата обращения 23.04.2019 г.).

7. Ермолаев В.И. Спутниковая платформа «Экспресс-1000». Уч. пос. / Под ред. В.А. Бабука, Н.А. Тестоедова. СПб.: Балтийский государственный технический университет, 2015. 67 с.

8. Ecoruspace. Экспресс-2000. Платформа аппарата. Режим доступа: https:// www.ecoruspace.me/Экспресс-2000.html (дата обращения 23.04.2019 г.).

9. Gunter's Space Page. SES-8. Режим доступа: https://space.skyrocket.de/doc_sdat/ ses-8.htm (дата обращения 23.04.2019 г.).

10. Gunter's Space Page. Thaicom-6. Режим доступа: https://space.skyrocket.de/ doc_sdat/ thaicom-6.htm (дата обращения 23.04.19 г.).

11. Gunter's Space Page. SES-14. Режим доступа: https .//space.skyrocket. de/doc_sdat/ses-14.htm (дата обращения 23.04.2019 г.).

12. Ecoruspace. Спутник связи Al Yah 3. Режим доступа: https://www.ecoruspace.me/ Al+Yah+3.html (дата обращения 23.04.2019 г.).

13. Завершено довыведение спутника «Экспресс-АМ5». Режим доступа: http ://www.iss-reshetnev.ru/media/news/ news-110314 (дата обращения 08.07.2019 г.).

14. Архангельский Н.И., Акимов В.Н., Кувшинова Е.Ю., Синицын А.А. Выбор параметров эллиптической орбиты базирования для повышения безопасности применения многоразовых ядерных буксиров // Космическая техника и технологии. 2016. № 2(13). С. 45-54.

15. Иванов Н.М., Лысенко Л.Н. Баллистика и навигация космических аппаратов // М.: Дрофа. 2004. 544 с.

16. Понтрягин Л.С., Болтянский А.Г., Гамкрелидзе Р.В., Мищенко Е.Ф. Математическая теория оптимальных процессов / Под ред. Л.С. Понтрягина. М.: Наука, 1976. 392 с.

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

18. Петухов В.Г. Оптимизация много-витковых перелетов между некомпланарными эллиптическими орбитами // Космические исследования. 2004. Т. 42. № 3. С. 260-279.

19. Моисеев Н.Н. Элементы теории оптимальных систем. М.: Наука, 1975. 528 с.

20. Kluever C. Simple Guidance Scheme for Low-Thrust Orbit Transfers // Journal of Guidance, Control, and Dynamics. 1998. V. 21. № 6. P. 1015-1017.

Статья поступила в редакцию 15.07.2019 г.

Reference

1. Polyarny P. Spasenie AEHF-1 [Saving AEHF-1]. Novosti kosmonavtiki, 2011, vol. 21, no. 12(347),p. 47.

2. AO «ISS» imeni akademika M.F. Reshetnyova» [JSC Academician M.F. Reshetnev Information Satellite Systems]. Available at: https://www.iss-reshetnev.ru/projects/telecommunication (accessed23.04.2019).

3. Cherny I. Spyatojpopytki [The fifth attempt]. Novosti kosmonavtiki, 2016, vol. 26, no. 05(400),pp. 27-28.

4. Khamits 1.1., Filippov I.M., Burylov L.S., Tenenbaum S.M., Perfil'ev A.V., Gusak D.I. Kontseptsiya kosmicheskoi transportno-energeticheskoi sistemy na osnove solnechnogo mezhorbital'nogo elektroraketnogo buksira [A concept of space transportation and power generating system based on a solar electric propulsion orbital transfer vehicle]. Kosmicheskaya tekhnika i tekhnologii, 2017, no. 1(16), pp. 32-40.

5. Gunter's Space Page. Hughes/Boeing HS-702. Available at: https://space.skyrocket.de/doc_sat/hs-702.htm (accessed 23.04.2019).

6. European Space Agency. Electra. Available at: https://artes.esa.int/news/electra (accessed 23.04.2019).

7. Ermolaev V.I. Sputnikovaya platforma «Ekspress-1000» [Spacecraft platform Ekspress-1000]. Ed. by V.A. Babuk, N.A. Testoedov. Saint Petersburg, Baltiiskii gosudarstvennyi tekhnicheskii universitet publ., 2015. 67p.

8. Ecoruspace. Ekspress-2000. Platforma apparata [Ecoruspace. Ekspress-2000. Apparatus platform]. Available at: https://www.ecoruspace.me^Kcnpecc-2000.html (accessed 23.04.2019).

9. Gunter's Space Page. SES-8. Available at: https://space.skyrocket.de/doc_sdat/ses-8.htm (accessed 23.04.2019).

10. Gunter's Space Page. Thaicom-6. Available at: https://space.skyrocket.de/doc_sdat/ thaicom-6.htm (accessed 23.04.2019).

11. Gunter's Space Page. SES-14. Available at: https://space.skyrocket.de/doc_sdat/ses-14.htm (accessed23.04.2019).

12. Ecoruspace. Sputnik svyazi Al Yah 3 [Ecoruspace. Communication satellite Al Yah 3]. Available at: https://www.ecoruspace.me/Al+Yah+3.html (accessed 23.04.2019).

13. Zaversheno dovyvedenie sputnika «Ekspress-AM5» [Satellite Ekspress-AM5 completion completed]. Available at: http://www.iss-reshetnev.ru/media/news/news-110314 (accessed 08.07.2019).

14. Arkhangel'skii N.I., Akimov V.N., Kuvshinova E.Yu, Sinitsyn A.A. Vybor parametrov ellipticheskoi orbity bazirovaniya dlya povysheniya bezopasnosti primeneniya mnogorazovykh yadernykh buksirov [Selecting parameters of elliptical basing orbit to improve safety of nuclear reusable tugs]. Kosmicheskaya tekhnika i tekhnologii, 2016, no. 2(13), pp. 45-54.

15. Ivanov N.M., Lysenko L.N. Ballistika i navigatsiya kosmicheskikh apparatov [Ballistics and spacecraft navigation]. Moscow, Drofa publ., 2004. 544 p.

16. Pontryagin L.S., Boltyanskii A.G., Gamkrelidze R.V., Mishchenko E.F. Matematicheskaya teoriya optimal'nykh protsessov [Mathematical theory of optimal processes]. Ed. by L.S. Pontryagin. Moscow, Nauka publ., 1976. 392 p.

17. Lebedev V.N. Raschet dvizheniya kosmicheskogo apparata s maloi tyagoi [The calculation of the motion of a spacecraft with low thrust]. Moscow, VTs AN USSR publ., 1968. 108p.

18. Petukhov V.G. Optimizatsiya mnogovitkovykh pereletov mezhdu nekomplanarnymi ellipticheskimi orbitami [Multi-revolution transfer optimization between non-coplanar elliptic orbits]. Kosmicheskie issledovaniya, 2004, vol. 42, no. 3, pp. 260-279.

19. Moiseev N.N. Elementy teorii optimalnykh system [Elements of the theory of optimal systems]. Moscow, Nauka publ., 1975. 528 p.

20. Kluever C. Simple guidance scheme for low-thrust orbit transfers. Journal of Guidance, Control, and Dynamics, 1998, vol. 21, no. 6, pp. 1015-1017.