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

УДК 629.7 Б01: 10.18287/2541-7533-2018-17-1-7-19

ОПТИМИЗАЦИЯ ТРАНСАТМОСФЕРНОГО ДВИЖЕНИЯ ЛЕТАТЕЛЬНОГО АППАРАТА НА ОСНОВЕ МЕТОДА ПРИНЦИПА МАКСИМУМА ПОНТРЯГИНА

© 2018

В. Л. Балакин доктор технических наук, профессор, советник при ректорате;

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

имени академика С.П. Королёва;

balakin@ssau.ru

С. А. Ишков доктор технических наук, профессор, профессор кафедры космического машиностроения;

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

имени академика С.П. Королёва;

ishkov@ssau.ru

А. А. Храмов кандидат технических наук, старший преподаватель кафедры космического машиностроения;

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

имени академика С.П. Королёва;

khramov@ssau.ru

Решается задача максимизации конечной скорости летательного аппарата при пассивном трансатмосферном движении. Заданы конечные значения высоты и угла наклона траектории. Управлением является коэффициент подъёмной силы (угол атаки). Ранее решение данной задачи оптимизации было проведено на основе метода последовательной линеаризации. Получено, что на больших высотах, в зависимости от выбранного начального приближения, программы управления несколько отличаются. Поэтому целью данной работы является определение оптимальной программы управления на основе «строгого» решения задачи оптимизации принципом максимума Понтрягина. Решение задачи оптимизации проведено на примере трансатмосферного набора высоты субгиперзвукового летательного аппарата MPV -первой ступени аэрокосмической системы RASCAL (США). Коэффициент подъёмной силы (угол атаки) на большей части траектории увеличивается для обеспечения заданных конечных значений высоты и угла наклона траектории, а затем уменьшается для обеспечения наибольшей конечной скорости. Подтверждением правильности полученных решений задачи оптимизации методом принципа максимума является нулевое значение Гамильтониана на оптимальной траектории. Полученные результаты показывают близость решений рассматриваемой задачи оптимизации методом принципа максимума и методом последовательной линеаризации.

Летательный аппарат; трансатмосферный набор высоты; оптимальное управление; коэффициент подъёмной силы; угол атаки; принцип максимума Понтрягина.

Цитирование: Балакин В.Л., Ишков С.А., Храмов А.А. Оптимизация трансатмосферного движения летательного аппарата на основе метода принципа максимума Понтрягина // Вестник Самарского университета. Аэрокосмическая техника, технологии и машиностроение. 2018. Т. 17, № 1. С. 7-19. DOI: 10.18287/2541-7533-2018-17-1-7-19

Введение

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

В работе [4] решение данной задачи оптимизации проведено на основе метода последовательной линеаризации [5]. В полученных программах управления коэффициент подъёмной силы (угол атаки) уменьшается к концу участка трансатмосферного движения. Было установлено, что на больших высотах, в зависимости от выбранного началь-

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

Поэтому целью данной работы является определение оптимальной программы управления летательного аппарата (ЛА) на основе «строгого» решения задачи оптимизации принципом максимума Понтрягина [6].

Постановка задачи

Модель движения. Считается, что ЛА движется над поверхностью Земли, имеющей форму сферы; поле притяжения принимается центральным. Движение аппарата относительно Земли происходит под действием силы тяготения и полной аэродинамической силы. Силами, обусловленными неинерциальностью системы отсчёта, пренебрегаем.

Система дифференциальных уравнений в траекторной системе координат имеет вид [7]:

Ж- Т/ 2 • о

= -ах р- -gsmв,

ш

-=Оу рУ+{ совв, (1)

Ж у ^ ЯГ)

ШЯ Г ■ в

— = V Бта. Ж

Здесь V - земная скорость ЛА; в - угол наклона траектории; Я - величина радиус-вектора центра масс ЛА; р - плотность атмосферы на высоте И = Я - Яср, задаваемая

экспоненциальным законом р = р0 ехр {-/И) ; р0 = 1,225 кг/м3 - плотность атмосферы у поверхности Земли; / = 0,141 -10-3м 1 - логарифмический градиент плотности; Яср = 6371,22км - средний радиус Земли; g = /л/Я2 - гравитационное ускорение; / = 398600,44 км3/с2 - гравитационная постоянная Земли.

Баллистические коэффициенты ах, ау определяются по соотношениям:

<>■.=ст, =ст£. (2)

2т 2т

где сха, суа - соответственно коэффициенты силы лобового сопротивления и аэродинамической подъёмной силы; £ - характерная площадь ЛА; т - масса ЛА.

Число Маха рассчитывается как отношение скорости аппарата к скорости звука а на данной высоте:

V

М = -. (3)

а

Скорость звука связана с температурой атмосферы Та соотношением:

а = 20,046796^, (4)

где температура атмосферы Та определяется согласно ГОСТ 4401-81 [8] как кусочно-линейная функция геопотенциальной высоты Игп, равная:

И = -^И. (5)

гп Я + И v;

Граничные условия движения. Для интегрирования дифференциальных уравнений движения (1) при решении задачи оптимизации в начальный момент времени t0

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

^ = 0, V = ¥0, в = в0, Я = Я = И0 + Яср . (6)

В конце участка движения аппарата в трансатмосфере заданы конечная высота и конечный угол наклона траектории:

t = Ч , И = Ик = Як - Яср, в = вк. (7)

Управление. ЛА на рассматриваемом участке движения управляется за счёт изменения угла атаки а и связанного с ним коэффициента подъёмной силы суа.

На угол атаки а и соответствующий ему коэффициент подъёмной силы суа наложены ограничения:

а ■ <а <а ,

mm max >

c ■ < c < c , (8)

ya mm ya ya max > \ /

гДе cya min и cya max определяются соответственно минимальным amin и максимальным amax углами атаки и зависят от числа M .

Задачу оптимизации движения сформулируем следующим образом: для системы дифференциальных уравнений (1) при заданных начальных (6) и конечных (7) условиях определить оптимальную программу изменения коэффициента подъёмной силы cya (t) при заданных ограничениях (8), обеспечивающую максимум конечной скорости ЛА.

Метод решения

Для определения оптимального управления воспользуемся принципом максимума Понтрягина [6].

Гамильтониан для системы (1) записывается в виде

dV d9 dR dt dt dt

или

H = ~¥v | ^ PV2 + ¥o {cmS PV + (V ~ W ] C0S °] + ¥rV sln 9. (9)

2m ' R2

V

Согласно принципу максимума оптимальное управление для принятого критерия должно обеспечивать минимум функции Гамильтона (9) в каждой точке траектории [6].

Примем параболическую зависимость коэффициента лобового сопротивления от коэффициента подъёмной силы:

с„ = с„0 (M) + A (M) c2ya,

(10)

где коэффициент лобового сопротивления cxa0 (М) и коэффициент отвала поляры A (M) есть функции числа Маха.

Тогда необходимое условие оптимальности управления запишется в виде

дH S

дcya 2ш

pV ((ЩуУА (М )суа + щв) = 0,

откуда:

c

¥в

ya ор1

2 A (М '

Достаточное условие оптимальности (минимума Гамильтониана) имеет вид:

= - * рУ 2¥гА (М )> 0.

дс„„ m

(11)

ya

Поэтому сопряжённый множитель на оптимальной траектории должен быть отрицательным:

щ < 0.

Запишем соотношения для сопряжённых множителей:

(12)

dt

-¥в р 2т V

ё¥в дН

ёх

дН

ёх дК

S +— - {

дН S дс т т2 пТТ

--н2т +2Ус-

S + с'

ya

дУ

+ 1 — + cos в

К У2К

2г,2

К2

Щу -У¥к Iсоэв + щ

- у/К sin в,

V А

2т

А

К3

V К УК2 А V

sin в,

(13)

УК3 К2

др

дИ

■¥в

др ~дйс>' +Р

дс }

ya

Системы уравнений (1) и (13) с управлением (11) определяют семейство оптимальных траекторий. Искомая оптимальная траектория будет определяться граничными условиями задачи. В соответствии с постановкой задачи в начале движения значения фазовых координат зафиксированы (6), в конце движения зафиксированы только высота (радиус-вектор) и угол наклона траектории (7), а скорость (функционал) и время движения являются свободными.

Тогда согласно [6] искомая экстремаль должна в конечный момент времени удовлетворять следующим условиям:

t = tk = var, И = кк, в = вк, уУ =-1, Н = 0.

(14)

Конечное время будем определять моментом достижения требуемого угла наклона траектории в конце полёта: в = вк. Это условие используется для окончания расчётов при математическом моделировании трансатмосферного движения ЛА.

Таким образом, задача оптимизации сводится к трёхпараметрической двухточечной краевой задаче. Требуется определить начальные значения сопряжённых множителей уу0, ув0, уЯ0, обеспечивающие для условий (6) выполнение конечных условий

(14). Для подбора параметров уу 0, ув0, уЯ 0 используются отклонения (невязки) в конечный момент времени параметров уу (tk ) , Н (tk ), И (tk ) от требуемых значений.

В систему уравнений для сопряжённых множителей (13) входят производные аэродинамических коэффициентов по скорости и высоте полёта и производная плотности атмосферы по высоте. Используя (11), для производной коэффициента подъёмной силы по скорости можно записать:

ув ( дЛ дЫ

дсу

2уУ \дЫ дУ

дУ

V + Л

(лу )2

(15)

Здесь с учётом (3) и (4):

дЫ

дУ 20,046796^

Используя (10), для производной коэффициента лобового сопротивления по скорости получим:

дсха дЫ

дУ дУ

дс

ха0

■ + с

дЫ уа

+ 2 Л ^

дЫ уа дЫ

(16)

Аналогично получаем производные аэродинамических коэффициентов по высоте:

дсуа = ув дЛ дЫ / Л2 дИ 2УуУ дЫ дИ '

1

dcxa DM

Dh Dh

Dcx

- + с

DM ya

dA-Ca + 2A dCya ^

DM

DM

J

(18)

где с учётом (3), (4) и (5):

DM Dh

V

DT„ Dh

V

R

2 • 20,046796^ Dh

Dh

40,093592^ V R + h.

DDT

Dhn

Производные температуры атмосферы по геопотенциальной высоте определяются по ГОСТ 4401-81 [8].

Выражения для производных коэффициента подъёмной силы по скорости (15) и по высоте (17) справедливы в случае, если значение коэффициента подъёмной силы, вычисленное по соотношению (11), лежит внутри области допустимого управления (8). При выходе из упомянутой области коэффициент подъёмной силы приравнивается граничному значению, а его производные вычисляются следующим образом:

DCya Dc . ya min DM

DV DM DV

DCya Dc ■ ya min DM

Dh DM Dh J

DCya _ Dc ya max DM ~

DV DM DV

DCya _ Dc ya max DM

Dh DM Dh J

ПРИ Cya < Cyamin(M, «mrn) ,

ПРИ Cya > Cyamax(M, «max )

Для производной плотности атмосферы по высоте для принятой экспоненциальной модели имеем:

= ~РРо ехр (-РИ).

оп

Результаты решения задачи

Решение задачи оптимизации трансатмосферного движения проведено на примере пассивного набора высоты субгиперзвукового ЛА MPV (MIPCC powered vehicle - аппарат с охлаждаемыми двигателями; Mass Injection Pre-CompressorCooling - система дополнительной подачи кислорода и охлаждения двигателя), являющегося первой ступенью АэрКС RASCAL (Responsive Access, Small Cargo, Affordable Launch - доступная система для запуска малых грузов по требованию, США) [1-3].

Исходные данные. Данные для ЛА приняты согласно [1], масса равна 30000 кг,

характерная площадь равна 251 м2.

Заданы следующие начальные: V0 = 1174 м/с, hK = 63160м, вк = 200 условия движения.

h0 = 26800 м и конечные:

Коэффициент лобового сопротивления сха0 (Ы) и коэффициент отвала поляры Л (Ы) рассчитывались по данным [3] (табл. 1) с использованием линейной интерполяции.

Таблица 1. Коэффициенты пассивного лобового сопротивления и отвала поляры

М 2 2,4 3 4

0xa 0 0,05 0,04 0,013 0,013

A 0,36332 0,55147 0,69157 0,85034

Ограничения на управление определяются в зависимости от числа Маха при минимальном amin = 0 и максимальном amax = 200 угле атаки:

min max

S« min (M , amin ) , S« max (M, Onax ) •

Зависимость граничных значений коэффициента подъёмной силы от числа Маха и

до . до

yamin yamax m / г

их производных —-- и —-- рассчитывались по данным [3J (табл. 2) с использо-

дМ дМ

ванием линейной интерполяции.

Таблица 2. Коэффициенты подъёмной силы при минимальном и максимальном угле атаки

М 2 2,4 3 4

0ya min (amin ) 0 0 0 0

0ya max (amax ) 0,87 0,68 0,52 0,42

В силу линейности зависимости eya (a) угол атаки, соответствующий коэффициенту подъёмной силы eya, рассчитывается по соотношению

a = (e _ e ) (amax ~ amin)

\ ya ya min / / \ *

(e _ e ■ )

\ ya max ya min J

Краевые задачи для всех вариантов решались модифицированным методом Ньютона [9]. При определении начального приближения для сопряжённых множителей использовались зависимости между их значениями. Гамильтониан (9), который должен быть равен нулю на оптимальной траектории, определяет одну из зависимостей. Соотношение для управления (11) даёт вторую зависимость в предположении, что коэффициент подъёмной силы задаётся (выбиралось значение, соответствующее максимальному аэродинамическому качеству). Таким образом, неизвестным остаётся один из трёх сопряжённых множителей, при задании которого определяются значения двух других. Поскольку множитель цу на оптимальной траектории согласно (12) должен быть отрицательным и его конечное значение равно цу (tk ) = _1, то целесообразно принять

Ц/у о <

Основной вариант. Этот вариант соответствует начальному углу наклона траектории ЛА согласно [1] - 00 = 38,4°. На рис. 1-3 показаны зависимости от времени соответственно параметров оптимальной траектории и угла атаки, аэродинамических коэффициентов, сопряжённых множителей.

70

V /

—, м / с, 20 60

в, град., Н, км, 50 а, град.

40

30

20

10

10

20

30 40

Время, с

п

в

а

50

60

70

Рис. 1. Зависимости параметров оптимальной траектории и угла атаки от времени при начальном угле наклона траектории в0 = 38,40

0

0

0,35

С с

ха1 уа 03

0,25 0,2 0,15 0,1 0,05 0

СУа^

С ха

10

20

30 40

Время, с

50

60

70

Рис. 2. Зависимости аэродинамических коэффициентов от времени при начальном угле наклона траектории в0 = 38,40

0

¥у,

щЛОТ щк -102

Wv

Wr

We

10

20

30 40

Время, с

50

Рис. 3. Зависимости сопряжённых множителей от времени при начальном угле наклона траектории в0 = 38,40

60

70

0

2

0

Отметим, что полученная с использованием принципа максимума траектория практически совпадает с траекторией, полученной с использованием метода последовательной линеаризации [4]. Оптимальный угол атаки к концу трансатмосферного участка движения несколько уменьшается (рис. 1), что качественно совпадает с результатами работы [4], но на всей траектории он близок к 100.

При использовании метода последовательной линеаризации угол атаки на большей части траектории близок к 100, но в конце траектории существенно зависит от начального приближения а0 (t) = const [4]. Например, угол атаки уменьшается до 00

при таких значениях начального приближения как а0 (t) = 00, а0 (t) = 10°, а0 (t) = 200, а

при а0 (t) = 150 уменьшается до 50.

Дополнительный вариант. Было проведено решение задачи оптимизации для несколько меньшего (00 = 350) начального угла наклона траектории.

Для меньшего значения начального угла наклона траектории должны быть несколько большие значения оптимального угла атаки. Это объясняется необходимостью увеличения подъёмной силы для увеличения угла наклона траектории и выполнения заданного для него конечного условия. Из-за увеличения подъёмной силы будет увеличиваться сила лобового сопротивления и уменьшится конечная скорость ЛА, время полёта будет больше.

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

Сравнение программ управления углом атаки, параметров траектории (рис. 1, 4) и коэффициентов аэродинамических сил (рис. 2, 5), полученных для разных значений начального угла наклона траектории, подтверждает приведённые ранее положения.

V / —, м/с. 20

в, град., Н, км, а, град.

70

60

50

40

30

20

10

10

20

п

в"

а

30

40 Время, с

50

60

70

80

Рис. 4. Зависимости параметров оптимальной траектории и угла атаки от времени при начальном угле наклона траектории в0 = 350

0

0

0,6

С , С 0,5

ха^ уа

0,4

0,3

0,2

0,1

С ^ уа

^ха

10

20

30 40 50

Время, с

60

70

80

Рис. 5. Зависимости аэродинамических коэффициентов от времени при начальном угле наклона траектории в0 = 350

0

0

0

Wy,

WR -102

-4 -6 -8 -10 -12

0 10 20 30 40 50 60 70 80

Время, с

Рис. 6. Зависимости сопряжённых множителей от времени при начальном угле наклона траектории в0 = 350

Подтверждением правильности полученных решений задачи оптимизации методом принципа максимума является и нулевое значение Гамильтониана (9) на оптимальной траектории.

Заключение

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

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

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

Библиографический список

1. Young D.A., Olds J.R. Responsive Access Small Cargo Affordable Launch (RASCAL) Independent Performance Evaluation // A Collection of Technical Papers - 13 th AIAA/CIRA International Space Planes and Hypersonic Systems and Technologies Conference. 2005. V. 1. P. 346-368. D01:10.2514/6.2005-3241

2. Urschel P.H., Cox T.H. Launch Condition Deviations of Reusable Launch Vehicle Simulations in Exo-Atmospheric Zoom Climbs // AIAA Atmospheric Flight Mechanics Conference and Exhibit. 2003. DOI: 10.2514/6.2003-5544

3. Young D. Responsive Access Small Cargo Affordable Launch (RASCAL) Independent Performance Evaluation.

https://www.yumpu.com/en/document/view/11944862/responsive-access-small-cargo-affordable-launch-rascal-.

4. Балакин В.Л., Ишков С.А., Храмов А.А. Оптимизация трансатмосферного движения летательного аппарата на основе метода последовательной линеаризации // Вестник Самарского университета. Аэрокосмическая техника, технологии и машиностроение. 2017. Т. 16, № 3. С. 17-26. DOI: 10.18287/2541-7533-2017-16-3-17-26

5. Федоренко Р.П. Приближённое решение задач оптимального управления. М.: Наука, 1978. 488 с.

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

7. Ярошевский В.А. Вход в атмосферу космических летательных аппаратов. М.: Наука, 1988. 336 с.

8. ГОСТ 4401-81. Атмосфера стандартная. Параметры. М.: Издательство стандартов, 1981. 180 с.

9. Салмин В.В., Ишков С.А., Старинова О.Л. Методы решения вариационных задач механики космического полёта с малой тягой. Самара: Самарский научный центр РАН, 2006. 162 с.

OPTIMIZATION OF FLIGHT VEHICLE TRANS-ATMOSPHERIC MOTION ON THE BASIS OF PONTRYAGIN'S MAXIMUM PRINCIPLE

© 2018

V. L. Balakin Doctor of Science (Engineering), Professor, Adviser to the University Administartion;

Samara National Research University, Samara, Russian Federation; balakin@ssau.ru

S. A. Ishkov Doctor of Science (Engineering), Professor, Professor of the Department of Space Engineering;

Samara National Research University, Samara, Russian Federation; irkt@ssau.ru

A. A. Khramov Candidate of Science (Engineering), Lecturer of the Department of Space Engineering;

Samara National Research University, Samara, Russian Federation; khramov@ssau.ru

The task of optimizing trans-atmospheric motion of a flight vehicle in order to maximize its final velocity with prescribed finite values of the height and flight path angle is considered. The angle of attack acts as control in passive motion of a vehicle. Previously, the sequential linearization method was used to solve this optimization task. It is shown that at great altitudes the control programs are slightly different depending on the chosen initial approximation. Therefore, the aim of this work is to determine the optimum control program on the basis of a "strict" solution of the optimization task using the Pontryagin's maximum principle. Solving the problem of optimizing trans-atmospheric motion of a flight vehicle is illustrated by passive climb of the sub-hypersonic vehicle MPV (the first stage of the aerospace system RASCAL designed in the USA). The coefficient of lift (angle of attack) increases in the greater part of the trajectory to provide the prescribed finite values of height and path inclination and then decreases to provide maximum final velocity. The correctness of the obtained solutions of the optimization task using the maximum principle is confirmed by the zero Hamiltonian value in the optimum trajectory. The results of vehicle motion simulation with optimal control and various initial conditions of motion and the vehicle mass are discussed. The results obtained show that the solutions of the optimization task under consideration using the maximum principle and the sequential linearization principle are in close agreement.

Flight vehicle; trans-atmospheric motion; optimal control; lift coefficient; angle of attack; Pontryagin's maximum principle.

Цитирование: Balakin V.L., Ishkov S.A., Khramov A.A. Optimizing a vehicle trans-atmospheric motion using Pontryagin's maximum principle // Вестник Самарского университета. Аэрокосмическая техника, технологии и машиностроение. 2018. Т. 17, № 1. С. 7-19. DOI: 10.18287/2541-7533-2018-17-1-7-19

References

1. Young D.A., Olds J.R. Responsive Access Small Cargo Affordable Launch (RASCAL) Independent Performance Evaluation. A Collection of Technical Papers - 13th AIAA/CIRA International Space Planes and Hypersonic Systems and Technologies Conference .V. 1. P. 346-368. D01:10.2514/6.2005-3241

2. Urschel P.H., Cox T.H. Launch Condition Deviations of Reusable Launch Vehicle Simulations in Exo-Atmospheric Zoom Climbs. AIAA Atmospheric Flight Mechanics Conference and Exhibit. 2003. DOI: 10.2514/6.2003-5544

3. Young D. Responsive Access Small Cargo Affordable Launch (RASCAL) Independent Performance Evaluation. Available at:

https://www.yumpu.com/en/document/view/11944862/responsive-access-small-cargo-affordable-launch-rascal-.

4. Balakin V.L., Ishkov S.A., Khramov A.A. Optimization of space vehicle trans-atmospheric motion by using the method of sequential linearization. Vestnik of Samara University. Aerospace and Mechanical Engineering. 2017. V. 16, no. 3. P. 17-26. DOI: 10.18287/2541-7533-2017-16-3-17-26

5. Fedorenko R.P. Priblizhennoe reshenie zadach optimal'nogo upravleniya [Approximate solution of optimal control problems]. Moscow: Nauka Publ., 1978. 488 p.

6. Pontryagin L.S., Boltyanskiy V.G., Gamkrelidze R.V., Mishchenko E.F. Matematich-eskaya teoriya optimal'nykh protsessov [Mathematical theory of optimal processes]. Moscow: Nauka Publ., 1976. 392 p.

7. Yaroshevskiy V.A. Vkhod v atmosferu kosmicheskikh letatel'nykh apparatov [Spacecraft atmospheric re-entry]. Moscow: Nauka Publ., 1988. 336 p.

8. GOST 4401-81. Standard atmosphere. Parameters. Moscow: Izdatel'stvo Standartov Publ., 1981. 180 p.

9. Salmin V.V., Ishkov S.A., Starinova O.L. Metody resheniya variatsionnykh zadach mekhaniki kosmicheskogo poleta s maloy tyagoy [Methods of solving variational problems of low-thrust mission mechanics]. Samara: Samarskiy Nauchnyy Tsentr RAN Publ., 2006. 162 p.