О некоторых особенностях поиска оптимального управления на основе принципа максимума для задачи некомпланарного межорбитального перехода Текст научной статьи по специальности «Математика»

УДК 629.783

БОТ 10.18698/2308-6033-2016-03-1475

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

© Е.В. Кирилюк1, М.Н. Степанов1'2

1МГТУ им. Н.Э. Баумана, Москва, 105055, Россия 2ФГБУ «4 ЦНИИ МО РФ», Юбилейный, 141091, Россия

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

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

Введение. В работе рассмотрена задача поиска оптимального управления в случае пространственного выведения орбитального блока (ОБ) с круговой низкой опорной орбиты (НОО) произвольного наклонения на некомпланарную ей орбиту с помощью реактивного двигателя большой ограниченной тяги. В качестве критерия оптимальности принята максимизация выводимой массы. Задача

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

Данная задача актуальна для стадии проектно-баллистических расчетов, связанных с определением максимально возможной массы полезного груза (ПГ), доставляемого средством выведения на целевые орбиты некоторого класса, или же с предварительным определением теоретической возможности отправки ПГ заданной массы на определенную целевую орбиту с помощью существующих средств выведения. Применение принципа максимума позволяет получить оптимальные схему выведения и ориентацию вектора тяги, которые могут быть адаптированы к особенностям функционирования системы управления средствами выведения. Схема выведения и соответствующая ей структура управления могут быть уточнены с учетом ряда факторов, не входящих в модель движения, используемую при постановке задачи оптимального управления, таких как особенности режима работы двигательной установки (выход на режим, импульс последействия [1]), нецентральность поля тяготения Земли [2] и пр. Как показывает практика, изменения в схеме выведения при этом несущественны.

В данной работе представлены результаты численного решения задачи оптимального пространственного межорбитального перехода между круговой НОО и произвольной целевой орбитой для диапазона углов некомпланарности от 10° до 30°. Продемонстрировано влияние эксцентриситета целевой орбиты на схему выведения (структуру активных участков работы двигательной установки), а также рассмотрены особенности решения задачи поиска оптимального управления с точки зрения принципа максимума для переходов на целевые эллиптические орбиты, симметричные относительно плоскости НОО.

Одним из распространенных подходов к определению неизвестных параметров краевой задачи принципа максимума является метод продолжения решения по параметру, сходимость которого зависит от удачного выбора начального приближения [3]. В работе предложен подход к пересчету компонент начального вектора сопряженных переменных, основанный на свойствах первого векторного интеграла задачи оптимального управления, существующего в рамках модели центрального поля тяготения Земли (ЦПТЗ). Интеграл позволяет перейти от имеющегося «исходного» решения к решению задачи выведения на целевую орбиту, симметричную «исходной» целевой относительно НОО. Данный подход в комбинации с методом продолжения решения по параметру в случае значительных отличий характеристик начальной и целевой орбит от орбит для некоторого

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

Постановка задачи. Рассмотрим задачу оптимального пространственного выведения ОБ с НОО на некомпланарную ей целевую орбиту в ЦПТЗ. При этом выведение будем анализировать как переход с орбиты на орбиту: точки старта и окончания выведения заранее не фиксируют, а получают в результате решения задачи.

Пространственное движение центра масс ОБ с двигателем большой тяги в ЦПТЗ описываем в абсолютной геоцентрической экваториальной системе координат (АГЭСК) системой дифференциальных уравнений

т> P о № г

Vx =— cosBcosy —3 = fv ,

m

r

МУ _

. p

Vy = — sin Bcosy= fv ,

mm r y

v =— sin у —3 = fv mm r

m = —,

X = V

Л y x->

y = Vy, z = Vz,

(1)

где х, у, 2 — координаты центра масс ОБ; Ух, Уу, У2 — проекции вектора скорости на соответствующие оси АГЭСК; у, , у — проекции вектора ускорения на соответствующие оси АГЭСК; ц = 3,986005 1014 м3/с2 — гравитационный параметр Земли; 0, у — углы ориентации вектора тяги (рис. 1);

Г = л[х

Ix2 + y2 + z2 — модуль радиус-вектора.

Для удобства и унификации расчетов без привязки к конкретному типу ОБ производят нормирование массы ОБ ()) и массового секундного расхода топлива (в) относительно начальной массы ОБ:

m (t) =

m(t) m(to)

P =

P

m(to)

Рис. 1. Ориентация вектора тяги в АГЭСК

Тяга двигательной установки (ДУ), фигурирующая в системе (1), с учетом нормирования, рассчитана следующим образом:

Р = wв,

где ^ — скорость истечения.

Множество допустимых управлений задано условиями

0 е[0,2п] ; у

; в е[0, в тах ]. (2)

п п ~2' 2.

В начальный момент времени t = t0 ОБ находится на круговой НОО высотой h над поверхностью Земли. Модель Земли — сфера радиусом R3 = 6371 км. Плоскость исходной орбиты наклонена к плоскости экватора под углом i; начальную плоскость определяем единичным вектором нормали n (sin i sin Q, - sin i cos Q, cos i). Граничные условия (ГУ) на левом конце траектории имеют вид:

\г (^)|=я+и, V(^)| = V(И) = т(¡0) = 1,

11 + И (3)

(г,У) (¡0 ) = 0, (г, п) (¡0 ) = о, (V, п) (¡0 ) = 0.

Положение ОБ после его выведения на целевую орбиту в нефиксированный заранее момент времени t = Т задано следующими ГУ в кеплеровых элементах на правом конце траектории:

а(Т) = ак, е(Т) = ек, ¡(Т) = 1к, О(Т) = Ок, ш(Т) = ш*. (4)

Под задачей оптимального управления будем понимать поиск программы изменения вектора тяги Р ^), т. е. такого управления

й(г) = |),y(t),|3^)^ , которое в конце участка выведения при

выполнении ГУ (3), (4) и ограничений (2) обеспечивает минимум функционала вида

I = Т - кт(Т) . (5)

Данный функционал представляет собой компромисс между затратами на выведение массы и затратами на сокращение времени полета [4], [5]. Коэффициент к [с] выступает как регулятор компромисса. При стремлении его к предельным значениям получаем:

при к ^ кш;п, кШщ > 0 — предельный переход к задаче быстродействия: конечная масса будет стремиться к абсолютному минимуму, соответствующему наискорейшему выведению;

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

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

Формализация задачи с точки зрения принципа максимума.

Задачей оптимального управления в принципе Понтрягина называют задачу, формализованный вид которой приведен в [6]. Ее решение сводится к решению краевой задачи для системы 2п обыкновенных дифференциальных уравнений — исходной системы и системы уравнений Эйлера. Математическая формализация задачи поиска оптимального управления с точки зрения применения принципа максимума выглядит следующим образом.

Функция Гамильтона — Понтрягина (ФГП) для системы (1) будет иметь вид:

Уравнения Эйлера для сопряженных переменных согласно условию стационарности по фазовым координатам запишем как

Н = РУХ /у, + РУ, V + РУ2/У2 + РхУх + РуУу + РгУг - РтР. (6)

Рх = "ГГ(г2 - 3x2)Ру, - 3х(УРу, + 2РУ;

Г у

Ру = "Г Г (Г2 - 3У2 )РУу - 3У(ХРУх + 2РУ; ) 1

Г 1_ у

pz = -- 322 )РУ; - 3хРУх + УРУУ )],

^ Г/-.2

(7)

1УХ = -Рх, РУу =-Ру > ^ =-Р;,

Рт = """2 Г РУС08 6cOSy + РУ sin0cosy + РУ■ siny ]

Х У 2 -I

(8)

(9)

Н1 (и ) = РУх — cos0cosy + РУ — sin0cosy + Ру — siny - Рт(3,

р

р

р

т

2 т

или

~ р

р

Н(и) = — Ру - ртв = —

т т

Ру

тр„

у

р

= — р; т

ру = ру^ соб 0 cos у + ру^ sin 0 cos у + ру7 sin у;

тРт

Р = Ру

ж,

(10)

(11)

Составляющая функции Н1 (и) - р, определяющая моменты включения и выключения ДУ, является функцией переключения [7]. Для достижения минимума функционала (5) согласно необходимым условиям оптимальности уравнения должны иметь следующий вид:

Р тах, Р > 0

0, р < 0, (12)

[0, Р тах ], Р = 0;

Ру Ру Г~2 2 2"

бш0собу = —-, siny = —-; Ру =. Ру + ру + ру .(13)

Ру Ру \ > у *

сов 0 СОБ у :

Ру, Ру

Системы (1), (7) краевой задачи имеют скалярный первый интеграл Н (^) = К0, являющийся следствием их автономности. Он не зависит от вида рассматриваемого гравитационного поля.

Из условия (9) следует, что Н(^) = А,0. Система также имеет векторный первый интеграл:

К = г х ^Г + V х ¡у, (14)

где рг (Рх, Ру, Р* ), Ру (Рух , Руу , Ру2) — векторы, сопряженные радиус-

вектору и вектору скорости соответственно; К(Кь К2, К3) — вектор констант.

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

'X - хь

Рис. 2. Взаимное расположение ВСК и АГЭСК

Векторный интеграл (14) в проекции на оси данной системы координат запишется в виде следующих соотношений:

Уь (г) Рь (г) - ^^ (г) Руь (г)+Ууь (г) Руь (г) - V(г) Рууь (г) = кх1, (15) ¿ь (г) Рь (г) - хь (г) ргЬ (г)+^ (г) рКь (г) - ^ (г) (г) = Куь, (16) хь (г)Руь (г) - Уь (г)рхь (г) + ухь (г)рууь (г) - Гуь (г)рКь (г) = КгЬ. (17)

Выражение (17) отражает условие трансверсальности по угловой дальности межорбитального перехода [4]. Так как на угловую дальность в рассматриваемой задаче ограничение не накладывается, то константу Кь можно конкретизировать: Кь = 0. Применение выбранной вспомогательной системы координат имеет дополнительные преимущества при определении соотношений между начальными векторами сопряженных переменных для симметричных межорбитальных переходов.

Переход между АГЭСК и ВСК осуществляется с помощью двух поворотов: по долготе восходящего узла и наклонению. Матрица перехода имеет следующий вид:

¿агэск—вск

СОБ & sinQ 0

- СОБI Бт & cos I cos& sin I Бт I Бт & - sin I cos& cos /

(18)

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

с11

агэск = Рхдх + РудУ + + РудУх + Руу д К + Ру,ду + Ртдт;

С1 вск = Рхь дхь + Руь дУь + Р?ь + Руь духь + Рууь дууь + РУь дуь + Ртдт.

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

дУ

сопряженных переменных в ВСК. Учитывая, что —- = 0, V/,

дЛ

дэ

- = {х, У, и —— = 0, V/, - = {х, у, запишем

дуЛ

уг д/

р-, р/ —;

V. дУ,

Ру, РУ- ={x, У, *>•

/ дУЯ

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

ргь = лагэск-вск рг, рУ, = лагэск—вск рУ •

(19)

Одно из ГУ на правом конце необходимо использовать для определения полного времени полета Т. При численном решении краевой задачи удобно использовать условие (8). Таким образом, можно понизить размерности краевой задачи.

При X0 Ф 0 условия (8), (9) можно заменить условием Рт (Т) = Н (Т)к = Н (г0 )к = К0к. Тогда с учетом произвольности величины Х0 подбор ее значения можно заменить проверкой условия Н (^ ) > 0 •

Для понижения размерности краевой задачи используем свойство однородности множителей Лагранжа и примем, что в начальный момент времени функция ру (¿0) = 1 (10), тогда из условий (13) следует

Рух ^0) = СОБ0(Г0)собУ(^0 );

Руу (10) = БШ 0(?0) СОБ у(^0); (20)

Ру* (?0) = «ПУ(^0).

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

Значение рх ) при заданных значениях х(ц), у(^), Ух (^), Уу (/0), ру (/0), рУх (/0), ру (^) определяется из проекции первого

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

векторного интеграла (17). В итоге приходим к краевой задаче размерности q = 5, для решения которой необходимо определить

значения 9(го), у(го), Ру (го), рг (го), Рт (го), обеспечивающие удовлетворение условий (4).

Предположим, что вектор состояния ОБ в начальный момент времени г = го полностью определен. Однако точка старта ОБ (первого включения ДУ) на НОО в данной постановке неизвестна. В качестве точки начала решения примем одну из следующих точек (А) заданной опорной орбиты:

точку А1( и = 27о°), для которой справедливы соотношения

х(г о ) = (Яз + Л^ш&соб I, у(г о ) = = -(Яз + k)cos&cos I, I(г о) = -(Яз + к^т I,

ух (г о) уу (г о) = ую = 0; (21)

]/Яз + к )]Яз + к

точку А2( и = 9о°), для которой справедливы соотношения

х(г о ) = -(Яз +k)sin&cos I, у (г о ) = = (Яз + к)cos О cos I, ю(г о) = (Яз + к^т I,

ух (г о) = уу (г о) = -ТТ^Т sin &, К! = 0. (22)

\Яз + к \Яз + к

Данное допущение должно обеспечить начало решения задачи с заведомо пассивного участка полета; выбор одной из указанных точек зависит от требований к ориентации вектора Лапласа целевой орбиты. Момент включения ДУ автоматически определяется из условия смены знака функции переключения р (11).

Численное решение краевой задачи принципа максимума проводили с применением метода Ньютона [8].

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

к(го ) = 200 км, &(го ) = 0°, ¡(го ) = 51,73° и семейством целевых орбит, задаваемых следующими ГУ:

Ип(Т) = [2000; 26371] км, Ъа (Т) = 26371 км, 0(Т) = 0°,

ю(Т) = 270°, /(Т) = 63,43°.

Иными словами, рассматриваем семейство оптимальных переходов для диапазона значений эксцентриситета целевой орбиты е(Т) = [0; 0,593] при неизменных круговой НОО и точке начала решения на ней — в данном случае точки Л2.

Результаты приведены для значения коэффициента компромисса

к = 10 . Численные исследования показали, что данное значение достаточно велико, чтобы задачу минимизации функционала (5) рассматривали в предельном случае, соответствующем максимизации массы выводимого ПГ без ограничения на полное время выведения. Дальнейшее увеличение коэффициента компромисса приводит к пренебрежимо малому увеличению конечной массы ОБ

(примерно 10-5...10-6 % от начальной массы ОБ).

Изменение полного полетного времени Т и продолжительности пассивного участка траектории (ПУТ) ЛТПут в зависимости от значения Ип целевой орбиты проиллюстрировано соответствующими графиками (рис. 3). Аналогичные зависимости для моментов первого включения ДУ на активном участке траектории (АУТ) Таут1 и времени начала ПУТ ТПУТ приведены на рис. 4; для длительностей включений ЛТАУТ1, Л7АУТ2 — на рис. 5.

Графики изменения аргументов широты соответствующих точек включения и выключения ДУ в зависимости от высоты перигея представлены на рис. 6, 7.

Рис. 3. Зависимости полной продолжительности выведения и длительности ПУТ от значения высоты перигея целевой орбиты при фиксированной высоте апогея:

--Т; - ДТПУТ

О некоторых особенностях поиска оптимального управления . 7\УТЬ ^ПУТ!С

2 000 1 500 1000 500 0

0,0£+00 5,0£+06 1.0Я+07 1,5£+07 2,0£+07 2,5£+07 ИП,м

Рис. 4. Зависимости времени первого включения ДУ и времени начала ПУТ от значения высоты перигея целевой орбиты при фиксированной высоте апогея:

— ТАУТ1;----ТАУТ2

?АУТЪ ^УТ2>С

500 400 300 200 100 0

0,0£+00 5,0£+06 1,0£+07 1.5Д+07 2,0£+07 2ДЕ+07 Ап, м

Рис. 5. Зависимости длительностей АУТ от значения высоты перигея целевой орбиты при фиксированной высоте апогея:

"луть "пут, град

150

100

50 0

0,0£+00 5,0£Ч-06 1,0^+07 1,5£Ч-07 2,0^+07 2,5£+07 йп, м -50Ь

-100 -150 -200

Рис. 6. Зависимости аргументов широты точек первого включения и выключения ДУ от значения высоты перигея целевой орбиты при фиксированной высоте апогея:

--МАУТЬ----иПУТ

«АУТ1. «ПУТ. ГРад

5г О

ОД 1+00 5,0£+06 1,0£+07 1,5£±(Л-ЪМгЪТЪЗШ07 /гп,м -5

у

у

/

/

/

/

/

/

/

/

-10 -15 -20 -25

Рис. 7. Зависимости аргументов широты точек начала и конца второго АУТ от значения высоты перигея целевой орбиты при фиксированной высоте апогея:

---- иАУТ2;--ик

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

Результаты расчетов при различных углах некомпланарности орбит. В качестве примера решения задачи оптимального межорбитального перехода для целевых орбит с различным наклонением рассмотрим переходы с круговой НОО (параметры ) = 200 км,

) = 0°, ¡^о ) = 50°) на два семейства высокоэллиптических целевых орбит одинаковой геометрии ЛП(Т) = 1000 км, Иа(Т) = 39 500 км, 0(Т) = 0°, ю(Т) = 270°:

I орбиты, для которых ¡(Т) > ¡(^ ),

¡(Т) = 60°, 65°, 70°, 75°, 80°;

II орбиты, для которых ¡(Т) < ¡(?0 ),

¡(Т) = 40°, 35°, 30°, 25°, 20°.

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

Для указанных целевых орбит I семейства было произведено численное решение краевой задачи принципа максимума с использованием метода продолжения по параметру; для целевых орбит II семейства применен метод продолжения по параметру и аналитический расчет неизвестных параметров краевой задачи с помощью пересчета начального вектора сопряженных переменных, полученного при численном решении задачи для соответствующих симметричных по наклонению орбит I семейства.

Предположим, что решение поставленной задачи оптимального пространственного межорбитального перехода для некоторого «исходного» положения целевой орбиты уже существует. Тогда необходимо оперативно определить решение задачи для положения целевой орбиты, симметричного «исходной» целевой относительно рассматриваемой плоскости НОО. Очевидно, что оптимальная схема перехода (структура АУТ ДУ) в рамках модели центрального поля тяготения не претерпит изменения от подобного «отражения» положения целевой орбиты относительно НОО.

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

Рг, {к ) = лагэск-вск Рг Рг, (/о ) = лагэск—вскру

где матрица ЛАГЭСК_ВСК имеет вид (18).

Пересчет начальных сопряженных векторов скорости производится исходя из соображения, что начальный сопряженный вектор скорости, согласно (20), есть единичный орт вектора тяги, т. е. определяет его оптимальную ориентацию. Так как симметричные по наклонению переходы отличаются в момент начала решения знаком при угле между ортом вектора тяги и его проекцией на плоскость НОО, то между компонентами векторов, сопряженных векторам скорости, для симметричных по наклонению переходов существует следующая связь:

РУхь (/о) = РУхЬ (/о);

РУуЬ (/о) = РУ, (/о); (23)

РУ, (/о) = -РУ, (/о).

Для определения начального сопряженного радиус-вектора р' (г0)

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

Кхт = -Кхт, 7 _ (24)

КуТ - - Ку1 •

Учитывая, что точка начала решения осталась прежней, из соотношений (23) и (24) следует, что для компонент векторов рг (¿0), р' ) выполняются соотношения

РхЬ (^0 ) = Рхт (^0 );

Рут (к ) = Рут (к ); (25)

РгИ (^0 ) = -Рт (^0 ).

Таким образом, для получения искомых векторов р' (г0), рV (^),

согласно (23), (25), следует заменить знаки при компонентах , рут

на противоположные и произвести обратный перевод сопряженных векторов из ВСК в АГЭСК:

ру (^0 ) = лагэск-вск • рУ1 оо); р' (^0 ) = лагэск-вск • ргь (^0).

Переменная, сопряженная массе, рт (г0) не изменяется в силу равенства (в рамках модели ЦГПЗ) энергетических затрат, необходимых для осуществления симметричных относительно НОО переходов. Таким образом, начальный сопряженный вектор для искомого оптимального межорбитального перехода может быть достаточно легко определен из начального сопряженного вектора для «исходного», симметричного по наклонению перехода.

Знание компонент вектора р' (г0) позволяет определить начальные углы ориентации 90, у0 вектора тяги в АГЭСК, фигурирующие в качестве искомых параметров краевой задачи, с помощью соотношений (20).

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

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

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

Таблица 1

Результаты решения краевой задачи принципа максимума для межорбитальных переходов между круговой НОО и высокоэллиптическими орбитами различного наклонения, I(Т) > I(/0)

1(Т), град рг (¿0 ) ру (А>) 9(?0 ), град У(*0 X град

Л с106 р 0106 Р.. с106 ру, с10' руу е10' ру.. 0106

60 249,02 -2656,21 3236,908 977 727,54 -137 358,58 -158 686,73 -7,997 -9,131

65 192,70 -1746,00 1473,312 986 742,67 -105 280,66 -123 510,70 -6,090 -7,095

70 202,52 -1565,40 1124,945 985 355,62 -110 252,18 -130 072,17 -6,384 -7,474

75 222,39 -1494,40 984,543 982 320,19 -120 848,52 -142 977,91 -7,013 -8,220

80 246,12 -1458,79 911,146 978 304,03 -133 604,89 -158 338,10 -7,777 -9,110

Таблица 2

Результаты решения краевой задачи принципа максимума для межорбитальных переходов между круговой НОО и высокоэллиптическими орбитами различного наклонения, I (Т) > I(¿0)

г(Т ), град рг (¿0 ) ру (70 ) 9(?0 ), град У(*0 ), град

р 0106 -т х 0 р 0106 г х0 р.. 0Ю6 рух 0106 руу 0106 руг 0106

40 249,02 3648,98 -2053,77 977 727,54 -132 423,85 -162 827,46 -7,713 -9,371

35 192,70 1754,12 -1463,63 986 742,67 -103 352,50 -125 128,62 -5,979 -7,188

30 202,52 1379,68 -1346,28 985 355,62 -108 951,00 -131 163,99 -6,310 -7,537

25 222,39 1229,09 -1300,73 982 320,19 -119 820,63 -143 840,41 -6,954 -8,270

20 246,12 1150,62 -1278,41 978 304,03 -132 732,34 -159 070,25 -7,726 -9,153

На рис. 8 и 9 для переходов I и II семейств представлены графики изменения сопряженных переменных в АГЭСК — ру (7), р. (7) ; на рис. 10 и 11 — графики изменения сопряженных переменных в ВСК, основная плоскость которой лежит в плоскости НОО — рУ1 (7), рг (7).

На рис. 12 и 13 — графики изменения ру (7), ру^ (7); на рис. 14 и 15 —

графики изменения ру(7), руь (7). Зависимости для переходов на

целевые орбиты, симметричные по наклонению, выделены одним цветом. При этом на рис. 8, 9, 11-13, 15 непрерывной линии соответствует случай перехода на орбиту I семейства, а пунктирной — на орбиту II семейства.

-2 000 -4 000 -6 000

500 1 000 1 500 2 000 2 500 3 000 3 500;: 4 000" 4 500 I, с

^ xv

v /

V ^

Рис. 8. Графики изменения сопряженной переменной ру для различных

значений наклонения целевой орбиты:

— ¡к - 80°;--¡к - 75°;--¡к - 70°;--¡к - 65°;--¡к - 60°;

--— ¡к - 40°;---¡к - 35°;---¡к - 30°;---¡к - 25°;---¡к - 20°

Рис. 9. Графики изменения сопряженной переменной р2 для различных значений наклонения целевой орбиты:

— ¡к - 80°;--¡к - 70°;--¡к - 75°;--¡к - 65°;--¡к - 60°;

---¡к - 40°;---¡к - 35°;---¡к - 30°;---¡к - 25°; — — ¡к - 20°

Рис. 10. Графики изменения сопряженной переменной руЬ для

различных значений наклонения целевой орбиты:

— ¡к = 60°, 40°;--¡к = 65°, 35°;--¡к = 70°, 30°;--¡к = 75°, 25°;

— ¡к = 80°, 20°

Аь-Юб

Рис. 11. Графики изменения сопряженной переменной ргЬ для различных значений наклонения целевой орбиты:

— 1к = 80°;--¡к = 7°;--¡к = 70°;--¡к = 65°;--¡к = 60°;

---¡к = 40°;---¡к = 35°;---¡к = 30°;---¡к = 25°;---¡к = 20°

Руу'

10°

4-10

2-10° -

-2-10° -

-4-10" -

500 1 ООО 1 500 2 ООО 2 500 3 000 3 500 4 000 4 500 г, с

----

x у

-6-10° ь

Рис. 12. Графики изменения сопряженной переменной руу для различных значений наклонения целевой орбиты:

— ¡к = 80°;--¡к = 75°;--¡к = 70°;--¡к = 65°;--¡к = 60°;

— ¡к = 40°;---¡к = 35°; -- — ¡к = 30°;---¡к = 25°;---¡к = 20°

Рис. 13. Графики изменения сопряженной переменной ру. для

различных значений наклонения целевой орбиты:

— ¡к = 80°;--¡к = 75°;--¡к = 70°;--¡к = 65°;--¡к = 60°;

---¡к = 40°;---¡к = 35°;---¡к = 30°;---¡к = 25°;---¡к = 20°

Ргл-™6

Рис. 14. Графики изменения сопряженной переменной ру 1 для

различных значений наклонения целевой орбиты:

— Ж = 60°, 40°;--Ж = 65°, 35°;--Ж = 70°, 30°;--& = 75°, 25°;

— Ж = 80°, 20°

М' юб

Рис. 15. Графики изменения сопряженной переменной руь для

различных значений наклонения целевой орбиты:

— 1к = 80°;--1к = 75°;--Ж = 70°;--Ж = 65°;--Ж = 60°;

---1к = 40°;---1к = 35°;---Ж = 30°;---Ж = 25°;--— ¿к = 20°

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

В процессе анализа результатов расчета оптимальных некомпланарных переходов на высокоэллиптические орбиты различного наклонения был предложен подход к пересчету начальных значений сопряженных переменных для переходов на целевые орбиты, симметричные относительно плоскости НОО. При симметричных межорбитальных переходах и использовании модели ЦПТЗ оптимальная схема перехода остается прежней, а оптимальная ориентация вектора тяги меняет знак между ортом вектора тяги и его проекцией на плоскость НОО. Значит, результат решения задачи с помощью алгебраического пересчета сопряженных переменных в том случае, когда необходимо найти решение для симметричного перехода, не является важным. Однако представленные соотношения имеют практическое значение при условии необходимости оперативного решения задачи определения оптимального пространственного перехода между двумя некомпланарными орбитами, а также в случае, когда в располагаемой базе решений отсутствует сколько-нибудь близкий переход, который можно использовать в качестве «исходного» для применения метода продолжения по параметру. Описанный выше прием позволяет быстро изменить угол некомпланарности между НОО и целевой орбитой и получить некое промежуточное решение из имеющегося «исходного». Относительно данного промежуточного решения удобно производить доведение до заданных ГУ на левом и правом концах траектории (описывающих НОО и целевую орбиту соответственно) с применением метода продолжения по отличающимся параметрам.

ЛИТЕРАТУРА

[1] Сердюк В.К. Проектирование средств выведения космических аппаратов. Москва, Машиностроение, 2009.

[2] Зеленцов В.В., Казаковцев В.П. Основы баллистического проектирования искусственных спутников Земли. Москва, Издательство МГТУ им. Н.Э. Баумана, 2012.

[3] Ивашкин В.В., Крылов И.В. Комплексный метод оптимизации космических траекторий с малой тягой и его применение к задаче перелета от Земли к астероиду Апофис. Препринты ИПМ им. М.В. Келдыша, 2011, № 56, с. 1-32.

[4] Григорьев К.Г., Федына А.В. Оптимальное пространственное выведение космического аппарата на геостационарную орбиту с орбиты искусственного спутника Земли. Техническая кибернетика, 1993, № 3, с. 116-126.

[5] Григорьев И.С., Данилина И.А. Оптимизация межорбитальных пространственных траекторий перелета ступенчатых космических аппаратов. Автоматика и телемеханика, 2007, № 8, с. 86-105.

[6] Григорьев К.Г., Григорьев И.С., Заплетин М.П. Практикум по численным методам в задачах оптимального управления. Дополнение 1. Москва, Издательство Центра прикладных исследований при механико-математическом факультете МГУ, 2007.

[7] Лысенко Л.Н. Наведение и навигация баллистических ракет. Москва, Издательство МГТУ им. Н.Э. Баумана, 2007.

[8] Волков Е.А. Численные методы. 5-е изд., Санкт-Петербург, Лань, 2008.

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

Ссылку на эту статью просим оформлять следующим образом:

Кирилюк Е.В., Степанов М.Н. О некоторых особенностях поиска оптимального управления на основе принципа максимума для задачи некомпланарного межорбитального перехода. Инженерный журнал: наука и инновации, 2016, вып. 3. URL: http://engjournal.ru/catalog/arse/adb/1475.html DOI 10.18698/2308-6033-2016-03-1475

Статья подготовлена по материалам доклада, представленного на XL Академических чтениях по космонавтике, посвященных памяти академика С.П. Королёва и других выдающихся отечественных ученых — пионеров освоения космического пространства, Москва, МГТУ им. Н.Э. Баумана, 26-29 января 2016 г.

Кирилюк Елена Владимировна родилась в 1993 г., студентка VI курса кафедры «Динамика и управление полетом ракет и космических аппаратов» МГТУ им. Н.Э. Баумана. Автор пяти работ в области динамики движения летательных аппаратов. e-mail: elena_vic_l@mail.ru

Степанов Михаил Николаевич родился в 1963 г., окончил ВИКИ им. А.Ф. Можайского в 1985 г. Канд. техн. наук, доцент кафедры «Динамика и управление полетом ракет и космических аппаратов» МГТУ им. Н.Э. Баумана. Автор более 30 научных работ в области баллистико-навигационного обеспечения и динамики движения летательных аппаратов. e-mail: step_m_n@mail.ru

Some aspects of solving the optimal control problem on the basis of the maximum principle for non-coplanar interorbital transfer

© E.V. Kiriluk1, M.N. Stepanov1, 2

:Bauman Moscow State Technical University, Moscow, 105005, Russia 24th Central Research Institute of the Ministry of Defence of the Russian Federation, Yubileynyy town, 141091, Russia

The article considers the approach to solve the problem of determining the optimal thrust vector control of propulsion system according to the criterion of maximization of the payload mass in the orbital unit launch into an arbitrary target orbit, which is non-coplanar to the original one. The problem is solved using Pontryagin's maximum principle, which reduces the problem of optimal control search to solve a boundary value problem. The successful solution of the problem depends on the quality of the initial approximation of its unknown parameters. In solving the problem the parameter continuation method was applied. The technique of operational conversion of the conjugate variable initial values using the symmetry properties of the target orbit inclination with respect to the plane of the initial reference orbit was also applied. _As a part of the work numerical modeling of the orbital block motion in the central gravitational field of the Earth with the optimal control, obtained on the basis of the maximum principle, was carried out. The dependence of the optimal launch scenario on the form (eccentricity) of the target orbit was analyzed. The behavior of conjugate variables, determining the vector of optimal control, which delivers a maximum output payload mass with set limits on the duration of transfer, was analyzed for a wide range of target orbit's inclination. The analysis revealed the symmetry properties of the behavior of conjugate variables, the use of which improves the efficiency of finding the optimal solution. The obtained results can be applied in the field of ballistic design calculations related to the development of launch vehicles (boosters, space tugs) and spacecraft for various purposes. They can also be applied for assessing the capacity of existing launch vehicles and determining the first approximation in the development of transfer schemes.

Keywords: Pontryagin's maximum principle, non-coplanar interorbital transfer, optimal control, launch.

REFERENCES

[1] Serduk V.K. Proektirovanie sredstv vyvedeniya kosmicheskikh apparatov [Designing Spacecraft Launch Vehicles]. Moscow, Mashinostroenie Publ., 2009, 504 p.

[2] Zelentsov V.V., Kazakovtsev V.P. Osnovy ballisticheskogo proektirovaniya iskusstvennykh sputnikov Zemli [Principles of Ballistic Design of Artificial Earth Satellites]. Moscow, BMSTU Publ., 2012, 176 p.

[3] Ivashkin V.V., Krylov I.V. Kompleksnyy metod optimizatsii kosmicheskikh traektoriy s maloy tyagoy i ego primenenie k zadache pereleta ot Zemli k asteroidu Apofis [Complex Method of Optimizing Space Trajectories with Low Thrust and Its Application to the Trip From The Earth to the Asteroid Apophis]. Preprints of Keldysh Institute of Applied Mathematics, 2011, no. 56, pp. 1-32.

[4] Grigoryev K.G., Fedyna A.V. Tekhnicheskaya kibernetika — Engineering Cybernetics, 1993, no. 3, pp. 116-126.

[5] Grigoryev I.S., Danilina I.A. Avtomatika i telemekhanika — Automation and Remote Control, 2007, no. 8, pp. 86-105.

[6] Grigoryev K.G., Grigoryev I.S., Zapletin M.P. Praktikum po chislennym metodam v zadachakh optimalnogo upravleniya. Dopolnenie 1 [Workshop on numerical methods in optimal control problems. Update 1]. Moscow, Center for Applied Research at the Mechanics and Mathematics Faculty of Moscow State University Publ., 2007, 184 p.

[7] Lysenko L.N. Navedenie i navigatsiya ballisticheskikh raket [Ballistic Missile Guidance and Navigation]. Moscow, BMSTU Publ., 2007, 672 p.

[8] Volkov E.A. Chislennye metody. 5-e izd. [Numerical methods. 5th ed.]. St. Petersburg, Lan Publ., 2008, 256 p.

Kiriliuk E.V. (b. 1993) the 6th year student, Department of Dynamics and Flight Control of Rockets and Spacecraft, Bauman Moscow State Technical University. Author of five articles in the field of flight dynamics. e-mail: elena_vic_l@mail.ru

Stepanov M.N. (b. 1963) graduated from Military Engineering Institute named after A.F. Mozhaisky in 1985. Cand. Sci. (Eng.), Associate Professor, Department of Dynamics and Flight Control of Rockets and Spacecraft, Bauman Moscow State Technical University. Author of more than 30 scientific articles in a field of ballistic navigation support and flight dynamic. e-mail: step_m_n@mail.ru