УДК 629.78.076.6 раСчет приближенно-оптимальных ... · 2020. 12....

94 КОСМИЧЕСКАЯ ТЕХНИКА И ТЕХНОЛОГИИ 4(27)/2019

Салмин В.В., Петрухина К.В., Кветкин А.А

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

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

Самарский национальный исследовательский университет имени академика С.П. Королёва (Самарский университет)

Московское шоссе, 34, г. Самара, Российская Федерация, 443086, e-mail: [email protected]

УДК 629.78.076.6

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

коСмичеСкого аппарата С двигателями малой тяги

С выСокоэллиптичеСкой на геоСтационарнуЮ орбиту

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

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: [email protected]

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

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

95 4(27)/2019 КОСМИЧЕСКАЯ ТЕХНИКА И ТЕХНОЛОГИИ

рАСЧЕТ прИбЛИжЕННО-ОпТИМАЛьНыХ пЕрЕЛЕТОв

Салмин в.в. петрухина к.в. кветкин а.а.

введение

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

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

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

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

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

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.

САЛМИН Вадим Викторович — доктор технических наук, профессор, директор научно-исследовательского института космического машиностроения Самарского университета, e-mail: [email protected] Vadim Viktorovich — Doctor of Science (Engineering), Professor, Director of the Research Institute of space engineering at Samara University, e-mail: [email protected]

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

КВЕТКИН Александр Александрович — аспирант Самарского университета, e-mail: [email protected] Aleksandr Aleksandrovich — Post-graduate at Samara University, e-mail: [email protected]



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

В настоящее время в США проводят-ся запуски геостационарных КА средне-го класса, использующих ЭРДУ для вы-полнения всех динамических операций, включая довыведение КА с геопере-ходной орбиты на ГСО на платформе 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

a > H

ГСО (H

a — высота апогея

орбиты; HГСО

— высота ГСО).

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



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

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

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

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

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

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

dt

dA

(1 – e)2

2p p= [esinϑS + (1 + ecosϑ)T],

µ

dt

de p= sinϑS + T

,

µ 1 + ecosϑ

ecos2ϑ + 2cosϑ + e

dt

di p= W,

µ 1 + ecosϑ

cosu

dt

dω p= – cosϑS + T – W

,

µ 1 + ecosϑ

sinϑ(2 + ecosϑ)

e

1

1 + ecosϑ

esinuctgi

dt

dΩ p= W,

µ sini(1 + ecosϑ)

sinu

dt

du µp= (1 + ecosϑ)2 –

ctgi

sinu

W ,

p2 (1 + ecosϑ)µ

p2

dt

dϑ µ=

,

pp

(1 + ecosϑ)

где p = A(1 – e2) — фокальный параметр; A — большая полуось; е — эксцентри-ситет; ϑ = u – ω — истинная аномалия; u — аргумент широты; ω — угловое расстояние перицентра от узла; Ω — дол-гота восходящего узла; i — наклонение орбиты; t — время; S, T, W — проекции возмущающего ускорения на направ-ление радиус-вектора, на перпенди-кулярное к нему в плоскости орбиты

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

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



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

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

Uсж

= – (3sin2isin2u – 1),3r3

ε

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

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

Sg = (3sin2isin2u – 1);

r4

ε

Тg = – sin2isin2u;

r4

ε

Wg = – sin22isinu.

r4

ε

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

При разложении вектора аэродина-мического ускорения КА по осям орби-тальной системы координат TSW можно не учитывать радиальную (S

f) и боковую

(Wf) составляющие аэродинамического

ускорения в силу их малости. Тогда вы-ражение для трансверсальной составляю-щей T

f примет вид:

T

f = σρ

µ

2p (1 + ecosϑ)2,

Sf = 0, W

f = 0,

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

CXS

mσ = ,

где СХ — коэффициент лобового сопро-

тивления; S — эффективная площадь миделя КА; m — масса КА.

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

S = Sr + S

g;

T = Tr + T

g + T

f;

W = Wr + W

g ,

где Sr, T

r, W

r — составляющие реактивного

ускорения КА.

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

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

dt

dω

dt

dΩ

dt

du, ,

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

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

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

Введем две правые системы коорди-нат: орбитальную (Onrb) и связанную с КА (OXYZ) (рис. 2). Вектор тяги

–p направлен

вдоль оси OX.

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

Tr = δacosλcosψ;

Sr = δasinλcosψ; (1)

W

r = δasinψ.

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



Здесь a — модуль полного реактивного ускорения (a = a

0 /(1 – a

0t/c)); δ — функ-

ция включения/выключения двигателей (δ = 0, 1); λ — угол между проекцией вектора тяги на мгновенную плоскость орби-ты и трансверсалью

–T (λ ∈ [–180°; 180°]);

ψ — угол отклонения тяги от мгновенной плоскости орбиты (ψ ∈ [–90°; 90°]) (рис. 2). Очевидно, что при постоянной работе двигателя (δ ≡ 1) время перелета равно моторному времени Т

м.

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

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

dt

di, содер-

жащей cosu.Запишем граничные условия: t = t

0 t = t

к

A(t0) = A

0

→ A(t

к) = A

к

e(t0) = e

0 e(t

к) = e

к

i(t0) = i

0 i(t

к) = i

к.

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

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

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

Ψ = (ΨA, Ψ

e, Ψ

i)T.

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

A(1 – e 2)H =

[A

SS + A

TT + A

WW],

µ

где

AS =

(1 – e)

2A(1 + e)esinϑΨ

A + sinϑΨ

e,

AT =

(1 – e)

2A(1 + e) (1 + ecosϑ)Ψ

A +

+ 1 + ecosϑ

ecos2ϑ + 2cosϑ + e Ψ

e,

AW

= 1 + ecosϑ

cosuΨ

i.

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

H = A(1 – e2)

µ×

× a[AScosλcosψ + A

Tsinλcosψ + A

Wsinψ].

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

∂λ

∂H =

A(1 – e2)

µ×

× a[–ASsinλcosψ + A

Tcosλcosψ] = 0;

∂ψ

∂H = A(1 – e2)

µ×

× a[–AScosλsinψ – A

Tsinλsinψ + A

Wcosψ] = 0.

Откуда получаем стационарные точки (2):

AT

AT

2 + AS

2s i n λ = ±

,

AS

AT

2 + AS

2cosλ = ±

.

AW

AT

2 + AS

2 + AW

2sinψ = ±

,

AT

2 + AS

2 + AW

2cosψ = ±

.

AT

2 + AS

2

(2)

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

x11

x12

x21

x22

∆ =

,

где x11

= ∂λ2

∂2H; x

12 = x

21 ∂λ∂ψ

∂2H; x

22 =

∂ψ2

∂2H.



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

(x

11x

22 – x

12x

21) > 0; x

11 < 0. (3)

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

AT

AT2 + A

S2

s i n λ = ,

AS

AT2 + A

S2

c o s λ = .

AW

AT2 + A

S2 + A

W2

sinψ = ,

AT2 + A

S2 + A

W2

cosψ = .

AT2 + A

S2

(4)

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

.Ψ

A = –

∂A

∂H = f

A(A, e, ϑ, u, λ, ψ, a, Ψ

A, Ψ

e, Ψ

i);

.

Ψe = – ∂e

∂H = f

e(A, e, ϑ, u, λ, ψ, a, Ψ

A, Ψ

e, Ψ

i);

.

Ψi = –

∂i

∂H = 0, Ψ

i = const.

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

сопряженных переменных ΨA0

, Ψe0

, Ψi0.

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

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

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

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

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

I = DxкTaDx

к → min, (5)

где Dxк = [DA, De, ∆i]T.

Здесь:

A( t) – Aк

A0

∆A = ;

De = e(t) – eк; i = i(t) – i

к;

α11

0 0

0 α2 2

0

0 0 α33

α = [αĳ] =

, Σα

ĳ = 1,

где aA = a

11, a

e = a

22, a

i = a

33 — весовые

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

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

dt

dI в каждый

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



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

I =

T

0 dt

dI

dt + I

0 → min.

Пусть I0 > 0. Тогда потребуем выполне-

ния условий:

sign(dt

dI) = const;

dt

dI

→ max.

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

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

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

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

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

.Э = f(Э, au), (6)

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

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

.Э = f(Э, 0) + aBu + 0(a2), (7)

где B =

∂f(Э, y)

∂y y = 0

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

производных.Отбросим в уравнении (7) малые вто-

рого порядка 0(a2) и рассмотрим урав-нение

.Э = f(Э, 0). Пусть его полный

интеграл имеет вид:

Э = g(t, C). (8)

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

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

∂C

∂gC = a

–1

Bu..

(9)

Рассмотрим функционал:

I(Э(T)) = I(g(T, C(T))) = I(C(T)). (10)

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

C = C0 + aC

1 + 0(a2),

причем C1 удовлетворяет уравнению (11):

C1 =

∂C

∂g –1

Bu, (11)

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

I*(C(T)) = I(C

0) + a

dC

dI*

C = C0

С1 + 0(a2).

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

H = ΨT

∂C

∂g –1

Bu = –a dC

dI*,

∂C

∂g –1

Bu .

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

dC

dI*,

∂C

∂g –1

Bu . (12)

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

(С(T))=

dC

dI*, С = a

dC

dI*,

∂C

∂g –1

Bu . (13)

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



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

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

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

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

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

dt

dI = 2a

1DA

A0

1

dt

dA +

+ 2a2De

dt

de + 2a

3Di

dt

di → max.

(14)

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

dt

dA,

dt

de,

dt

di:

AS* = a

1DA

A0(1 – e)

2A(1 + e) esinϑ + a

2Desinϑ;

AT* = a

1DA

A0(1 – e)

2A(1 + e) (1 + ecosϑ) +

+ a2De

1 + ecosϑ

ecos2ϑ + 2cosϑ + e;

AW* = a

3Di

1 + ecosϑ

cosu

.

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

dt

dI = 2

A(1 – e2)

µ [A

S*S + A

T*T + A

W*W].

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

dλ

∂

dt

dI = 0;

dψ

∂

dt

dI = 0.

Результатом поиска максимума

dt

dI =

dt

dI(λ(t), ψ(t)) по двум переменным

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

AT*

( AT*) 2 + ( A

S*) 2

s i n λ = ,

AS*

( AT*) 2 + ( A

S*) 2

c o s λ = .

(15)

AW*

(AT*)2 + (A

S*)2 + (A

W*)2

sinψ = ,

cosψ = .

(AT*)2 + (A

S*)2

(AT*)2 + (A

S*)2 + (A

W*)2

(16)

Полученный закон управления ~ψ(t), ~λ(t) имеет достаточно простую структуру и позволяет провести расчет динамиче-ского маневра без процедуры решения краевой задачи.

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

A,

ae, a

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

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

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

A = a

e = a

i.

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

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



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

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

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

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

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

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

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

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

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

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

a = [40 000; 80 000] для i = 28°

и ra = [70 000; 110 000] для i = 51,6°. Ис-

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

0 = const).

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

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

Способ решения

Параметры перелета

r0 = 20 000 км

rк = 23 350 км

Di = 19,022°

a = 0,00498 м/с2

r0 = 50 000 км

rк = 58 375 км

Di = 19,022°

a = 0,00080 м/с2

r0 = 80 000 км

rк = 93 400 км

Di = 19,022°

a = 0,00031 м/с2

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

Приближенный 5,2416 20,648 41,760

Точный 5,1580 20,389 41,264

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

Примечание. r0 — радиус начальной орбиты; r

к — ра-

диус конечной орбиты; Di — изменение наклонения;

а — модуль полного реактивного ускорения.

Способ

решения

Исходные данные

ra0 = 42

171

км

rп0 = 6

871

км

i0 = 75°

ra0 = 42

378

км

rп0 = 6

578

км

i0 = 7°

ra0 = 46

500 км

rп0 = 6

642,9

км

i0 = 7°

ra0 = 34

171

км

rп0 = 6

595

км

i0 = 63,17°

rк = 42

165

км r

к = 42

378

км r

к = 42

378 км r

к = 42

160

км

m0 = 1

320 кг

Iуд

= 1 500 с

P = 0,332 Н

m0 = 2

000 кг

Iуд

= 2 000

с

P = 0,350 Н

m0 = 1

500 кг

Iуд

= 1 994,06 с

P = 0,200 Н

m0 = 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

Примечание. ra0 — начальный радиус апогея; r

п0 —

начальный радиус перигея; rк — радиус конечной ор-

биты; m0 — начальная масса КА; I

уд — удельный

импульс; Р — тяга.

Таблица 1

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

Таблица 2

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



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

В соответствии с правой частью дифференциального уравнения

dt

di с воз-

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

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

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

Наклонение орбиты i, °

Высота перигея Н

п, км

Стартовая масса КА

m0, кг

Масса спутника на ГСО

mКА

, кг

Масса рабочего

тела mРТ

, кг

Тяга двигателя

P, мН

Удельный импульс

Iуд

, с

Диапазон перебора r

a, тыс. км

28

200 3 500 2 600 900 360 1 600

40…80

51,6 70…110

i, ° Элемент орбиты

Необходимое значение

Результат

28

e 0 291,22 сут

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

i 0° 291,72 сут

Tк

— 291,72 сут

Масса рабочего тела

— 567,11 кг

51,6

e 0 335,15 сут

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

i,° 0 335,56 сут

Tк

— 336,19 сут

Масса рабочего тела

— 653,56 кг

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

Таблица 3

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

a) б)

Рис. 3. Результаты анализа оптимального радиуса апогея: а — для наклонения i = 28°; б — для наклонения i = 51,6°

Таблица 4

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



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

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

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

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

Данный факт можно объяснить тем, что сила аэродинамического сопро-тивления F

A в нижней точке траекто-

рии (наибольшее сопротивление на пер-вом витке) примерно равна силе тяги P (F

Aπ = 320 мН; P = 360 мН), а длитель-

ность и протяженность участка воздей-ствия крайне малы. За счет этого высота перигея орбиты выходит за пределы дей-ствия атмосферы за 3–4 витка.

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

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

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

28 51,6

aA

0,082 0,073

ae

0,285 0,456

ai

0,633 0,471

Таблица 5

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

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

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

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

0 = 28°



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

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

Зависимость управляющего угла λ для случая перелета с i

0 = 51,6° носит

колебательный характер, угол меняется в пределах –180…180°.

выводы

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

п = 200 км), приведены графики, ото-

бражающие характер изменения орби-тальных параметров и углов ориентации вектора тяги ЭРДУ.

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

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

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

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

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

0 = 51,6°

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

0 = 51,6°



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

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

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

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

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

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

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 г.

reference1. Polyarny P. Spasenie AEHF-1 [Saving AEHF-1]. Novosti kosmonavtik, 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 (accessed 23.04.2019).3. Cherny I. S pyatoj popytki [The fifth attempt]. Novosti kosmonavtiki, 2016, vol. 26, no. 05(400), pp. 27–28.



4. Khamits I.I., 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. 67 p.8. Ecoruspace. Ekspress-2000. Platforma apparata [Ecoruspace. Ekspress-2000. Apparatus platform].

Available at: https://www.ecoruspace.me/экспресс-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 (accessed 23.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. 108 p.

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.

Date post:	11-Aug-2021
Category:	Documents
Upload:	others
View:	2 times
Download:	0 times

УДК 629.78.076.6 раСчет приближенно-оптимальных ... · 2020. 12....

Documents