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

УДК 519.6

БОТ: 10.18698/2309-3684-2021-3-7487

Моделирование и оптимизация управлением спутника малой массы при перелете с орбиты Земли на орбиту Марса под солнечным парусом

© Т.Ю. Мозжорина, Д.А. Рахманкулов

МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

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

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

Идея использования давления солнечного света для космических перелетов возникла еще в 19-ом веке. Попытка разработки реального проекта была осуществлена в России в 20-х годах двадцатого века. Принадлежит Фридриху Цандеру, одному из основоположников ракетостроения. Давление света очень мало и уменьшается пропорционально квадрату расстояния от Солнца. Солнечный парус, как движитель космического аппарата даже небольшого по массе, должен иметь большую площадь поверхности и быть при этом изготовлен из сверхлегких материалов. Разработка проектов солнечного паруса во многих странах велась с начала 90-х годов прошлого века. Но только

в этом столетии впервые были осуществлены и выведены на орбиту первые спутники с солнечными парусами на борту. Первым успешным проектом был японский IKAROS [1], продемонстрировавшим управляемое маневрирование с помощью солнечного паруса. Запущен в 2010 году. В качестве отражающей поверхности он использует квадратную полиимидную пленку (каптон, производства DuPont), состоящую из четырех трапециевидных фрагментов. Толщина паруса составляет 7,5 микрон, но в нее дополнительно вшиты тонкопленочные солнечные батареи, предназначенные для генерации электричества. Парус представляет собой квадрат со стороной 14 метров. Процесс раскрытия занял семь дней, после чего IKAROS отправился к Венере.

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

Первым наиболее успешным американским проектом стал аппарат LightSail-2 [1]. 23 июля 2019 года успешно раскрыл солнечный парус, а 31 июля ему удалось достигнуть роста максимальной высоты орбиты со скоростью около километра в сутки. Его парус из майлара имеет площадь 32 квадратных метра. Таким образом, реальное воплощение идеи межпланетных перелетов без затрат топлива начала осуществляться только в последние годы. Этим объясняется всплеск интереса к моделированию возможных перелетов под солнечным парусом и с двигателями малых тяг в пределах солнечной системы [2, 3].

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

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

Моделированию движения спутников с солнечным парусом в пределах солнечной системы посвящены и работы последних лет [8-12].

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

1

J = |йг = Т ^ шт.

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

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

• считаем орбиты круговыми и лежащими в одной плоскости.

Уравнения движения в гелиоцентрической системе координат:

^) ОС83 в,

Я )

йи йг V2 Г =У-— а я а ^ Я)

dV т =---а Я Г яо

йг 1 Я

йЯ йг = и,

йр _ V

йг = Я'

Здесь и — радиальная составляющая скорости спутника в гелиоцентрической системе координат (рис. 1), м/с; V — тангенциальная составляющая скорости, м/с; Я — расстояние от Солнца до спутника, м; р — полярный угол спутника, рад; Т — время перелета, с; А0 — ускорение, придаваемое спутнику силой притяжения к Солнцу на орбите Земли, м/с2; а — ускорение, придаваемое спутнику силой давления света Солнца на орбите Земли, м/с2; в — угол установки солнечного паруса, управление, рад, измеряется от тангенциальной составляющей скорости космического аппарата. По часовой стрелке в направлении от тангенциальной скорости идут положительные углы, против часовой стрелки — отрицательные. На рис.1 изображена сила Т , действующая на космический аппарат от солнечного паруса. При таком действии силы парус будет установлен в положение отрицательных углов в. Ограничения на угол установки солнечного паруса к к "2'2 _

Краевые условия: V, и =0, Я, (=0 — значения составляющих орбитальной скорости Земли, радиус Земной орбиты и полярный угол в выбранной системе координат в начальный момент времени,

ве

V, и = о

Я

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

Марса, радиус орбиты Марса в момент окончания перелета.

0

2

2

<

Орбита Марса

Орбита Земли

Солнце

Рис. 1. Схема перелета между орбитами Функция Понтрягина, соответственно, будет равна:

Н = Wu

+ Wy

R - * i RJ Ч RJ *

f - — -aiRY sin*cos2 * R l R J

y

Л

+wru + wç--1.

Здесь Wu, Wy, WR, Wv — сопряженные переменные.

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

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

Функция Понтрягина не линейна по управлению, поэтому можем воспользоваться уравнением Эйлера по управлению:

ÔH

дв

■ = О

^ -3 • \уиа

Яо

v я у

• СОБ2 в • БШ в —

-Wvа

Я0) •( СОБ3 в — 2 СОБ в • бш2 в) = 0.

Отсюда имеем два решения относительно ного угла установки паруса):

(тангенс оптималь-

* 3Уи ±>/

*§в =-1Т~2-.

Из соображений практического характера в начале перелета угол установки солнечного паруса должен находиться в четвертой четверти (отсчет угла ведется по часовой стрелке, см. рис.1 и систему уравнений движения), то есть быть отрицательным. Окончательно получим:

в = аг^ Сопряженная система:

й^и _ WvV

=-л—2-;

3¥и — V 9Уи 2 + 8Wv V

Ч2

йг

dWv

Я

чи ¥р

--1----.

йг Я Я Я

йг ¥и

Су2

Я2

2АЯ? 2аЯ2 СОБ3 в

Я3

Я

(

Wv

йг

Ш 2аЯ2 • СОБ2 в • МП в

Я2

= 0.

Я3

+ ¥(

V

р Я2

По условию транверсальности имеем ч (Т) = 0. Учитывая

постоянство этой сопряженной переменной по времени, можем исключить последнее уравнение из системы и исключить все слагаемые, содержащие у. Окончательно Пи-система (система ДУ, содержащая уравнения движения и сопряженную систему с

подставленным в нее оптимальным управлением) будет содержать 7 уравнений и примет вид:

'du.=VI _ 4 Г r. Y+a\ R. YCos3 <r,

dt R 01R J \R J

dV UV Г R Л2 . ..

— =---a\ — I sin0 cos20 ,

dt R I R J

dR

= U,

dt V

dt ~ R' d^u _ VvV

dt

R

■¥r ,

d^ _ 2^uV + WvU

dt

dt WU

R

2

R

V2 2AR 2aR2 cos3 0 R2 R3 R3

Г

~Wv

UV 2aR cos2 0 sin0 R2 R3

Л

Алгоритм численной реализации математической модели.

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

<

Ньютоновском направлении использовалась не стандартная евклидова норма по невязкам, а так называемая локальная норма, введенная Федоренко Р.П. [5], рассчитываемая по формуле:

N =

^ %

1 Г , „Л2

=1

I

з =1

д%

кдхз у

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

Был введен дополнительный параметр пристрелки — Т, время перелета. Система ДУ была переписана с учетом введения нового аргумента для более корректного выхода из метода численного интегрирования: г1 = г / Т, г1 е [0,1].

Параметры пристрелки в рассматриваемой задаче: щ, щ, щ, Т. Невязки:

% = и(Т) ^ 0, % = (V(Т) — Vf) / Vf ^ 0,

% =( Я(Т)—я ) / Я ^ 0, % = Н(Т) ^ 0.

Программа расчета была написана на языке программирования С++ с описанием переменных двойной степени точности. Точность по невязкам была выбрана £ = 10—6. Точность интегрирования была проверена путем обратного интегрирования системы ДУ от конца к началу для сошедшегося варианта. Совпадение составило по всем переменным по всем шагам интегрирования (всего выбрано 1000 шагов) не менее 6 значащих цифр.

Результаты расчета. В качестве начального приближения параметров пристрелки были выбраны следующие значения: щ = 695, щ =1383, щ= —10—5, Т = 400. Значения подобраны из соображения равенства нулю функции Понтрягина в начальный момент, анализа системы ДУ и реальным данным по перелету спутников с орбиты Земли на орбиту Марса. На рис. 2-5 представлены графики изменения фазовых переменных в процессе перелета с одной орбиты на другую. На рис. 6-8 — графики изменения сопряженных переменных. На рис. 9 — оптимальное значение угла установки солнечного паруса. Параметры пристрелки, найденные в процессе решения, равны: щ (0) = 715,83, щ =1319,82, щ = 0,0002485. Время перелета составило

403,5 суток. Решение найдено за 16 итераций метода Ньютона. Полученные результаты близки к графическим данным, приведенным в [4]. На рисунке 10 указана траектория спутника в Солнечной системе. Следует отметить, что орбиту Марса спутник пересекает в двух местах. В первый раз процесс не остановлен в силу несоответствия краевым условиям составляющих скорости космического аппарата. В таблице 1 приведены значения фазовых переменных в конце процесса и в точке пересечения орбиты Марса в первый раз в оптимальном процессе.

Функция Понтрягина в течении всего перелета сохраняет нулевое значение с точностью до 10-6, что соответствует теории: в автономных системах со свободным временем на правом конце функция Понтря-гина равна нулю в течении всего процесса.

и-103,

7 6 5 4 3 2 1

-1

0

100

400

500

Рис.

200 300

г, сутки

2. Радиальная составляющая скорости полета в оптимальном процессе

V-103, м/с

32

30 28 26

24

22

20

100

200 г, сутки

300

400

500

Рис. 3. Тангенциальная составляющая скорости полета в оптимальном процессе

0

0

Я-101-, м 2,4

2,3

2,1

1,9

1,7

1,5

р, рад 4,5

4,0

3,5

3,0

2,5

2,0

1,5

1,0

0,5

100

300

400

200 г, сутки

Рис. 4. Радиус в гелиоцентрической системе координат в оптимальном процессе

500

200 300

г, сутки

500

Рис. 5. Изменение полярного угла в гелиоцентрической системе координат в оптимальном процессе

1,0 0,5 0,0 -0,5 -1,0 -1,5 -2,0

0

100

400

500

200 300

г, сутки

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

2

2

1

1

0

0

Рис. 7. Изменение сопряженной переменной, соответствующей тангенциальной скорости полета, в оптимальном процессе

¥я -ю-3

Рис. 8. Изменение сопряженной переменной, соответствующей радиусу, в оптимальном процессе

в, град 0

-80 -- —- —- —

0 100 200 300 400 500

t, сутки

Рис. 9. Управление (угол установки солнечного паруса) в оптимальном процессе

Рис. 10. Схема перелета с орбиты Земли на орбиту Марса:

• — Солнце;--орбиты планет; — — траектория полета;

х — точки пересечения орбиты Марса

Таблица 1

Сравнение фазовых переменных в первой точке пересечения орбиты Марса со значениями в конце процесса

Переменные U, м/с V, м/с R, м t, сутки р, рад

Первая точка пересечения орбиты Марса 1096 22087 227,9-109 248,6 2,972

Конец перелета 0 24100 227,9-109 403,5 4,312

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

ЛИТЕРАТУРА

[ 1 ] Андреев А.А. Регата под солнечным парусом. STIмул. Журнал об инновациях в России [Электронный ресурс], 2019. URL: https://stimul.online/articles/sti-ence-and-technology/regata-pod-solnechnym-parusom/ (дата обращения: 14.06.2021)

[2] Григорьев И.С., Заплетин М.П., Самохин А.С., Самохина М.А. Оценка возможного выигрыша по массе при использовании двигателей малой тяги в экспедиции к Марсу. Сборник докладов XI Всероссийского съезда по фундаментальным проблемам теоретической и прикладной механики, Казань, 2015, с. 1063-1065.

[3] Ишков С.А., Старинова О.Л. Оптимизация и моделирование движения космического аппарата с солнечным парусом. Известия Самарского научного центра Российской академии наук, 2005, т. 7, № 1, с. 99-106.

[4] Лейтман Дж. Методы оптимизации с приложениями к механике космического полета. Москва, Наука, 1965, 538 с.

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

[6] Мозжорина Т.Ю. Численное решение задач оптимального управления с переключением методом пристрелки. Математическое моделирование и численные методы, 2017, № 2 (14), с. 94-106.

[7] Мозжорина Т.Ю., Осипов В.В. Численное решение задачи о мягком приземлении методом пристрелки. Инновационное развитие, 2018, № 8 (25), c. 11-15.

[8] Федоренко А.Н. Математическое моделирование переориентации орбитального космического аппарата со сферическим солнечным парусом. Дисс. канд. техн. наук. Москва, 2014, 126 с.

[9] Чернякина И.В. Программы локально-оптимального управления и траектории гелиоцентрических перелетов космического аппарата с солнечным парусом с учетом возмущений. Дисс. канд. техн. наук. Самара, 2020, 138 с.

[10] Старинова О.Л., Горбунова И.В. Оптимизация гелиоцентрического движения космического аппарата с солнечным парусом. Сборник трудов XVII Всероссийского семинара по управлению движением и навигации летательных аппаратов, Самара, 2015, с. 168-171.

[11] Starinova O.L., Kurochkin D.V., Materova I.L. Optimal control choice of non-keplerian orbits with low-thrust propulsion. AIP Conference Proceedings, 2012, no. 1493 (1), pp. 964-971.

[12] Сомов Е.И., Бутырин С.А., Старинова О.Л. Регулируемый солнечный парус в системе управления движением мини-спутника землеобзора: модели и оценки эффективности. Сборник трудов 7-ой Российской мультиконферен-ции по проблемам управления «Управление в морских и аэрокосмических системах (УМАС-2014)», Санкт-Петербург, 2014, с. 545-556.

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

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

Мозжорина Т.Ю., Рахманкулов Д.А. Моделирование и оптимизация управлением спутника малой массы при перелете с орбиты Земли на орбиту Марса под солнечным парусом. Математическое моделирование и численные методы, 2021, № 3, с. 74-87.

Мозжорина Татьяна Юрьевна — канд. техн. наук, доцент кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. e-mail: [email protected]

Рахманкулов Дмитрий Александрович — студент кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. e-mail: [email protected]

Modeling and optimization of low-mass satellite control when flying from Earth orbit to Mars orbit under a solar sail

© T.Yu. Mozzhorina, D.A. Rakhmankulov Bauman Moscow State Technical University, Moscow, 105005, Russia

T.W. Moswopuna, ff.A. PaxManKyrne

In this paper, the optimization of the transfer of a low-mass satellite from the Earth's orbit to the Mars orbit under a solar sail is considered. Optimization of the control of the pitch angle of the solar sail is carried out using the Pontryagin maximum principle while minimizing the flight time. In contrast to previous works on this topic, the solution of the boundary value problem, to the solution of which the maximum principle is reduced, was obtained by the false position method. The calculation program is written in the C++ programming language. Despite the computational difficulties arising when using the false position method, it was possible to achieve good convergence of the Newton method underlying the algorithm. The analysis of the accuracy of the results obtained is carried out and the possibility of using the false position method in solving such problems is shown. A comparison is made with the data of previously published works. Despite some assumptions used in the development of the calculation algorithm, the work has its value in terms of assessing the possibility of using the false position method, which gives the most accurate numerical optimization results.

Keywords: solar sail, the false position method, boundary value problems of ordinary differential equations, optimal control, Pontryagin's maximum principle, flight between the orbits of the Earth and Mars

REFERENCES

[1] Andreev A.A. Regata pod solnechnym parusom [Regatta under a solar sail]. STImul. Zhurnal ob innova-ciyah v Rossii [Stymul. Journal of Innovations in Russia] [Electronic resource], 2019. URL: https://stimul.online/articles/science-and-technology/regata-pod-solnechnym-parusom / (accessed: 06.14.2021)

[2] Grigoriev I.S., Zapletin M.P., Samokhin A.S., Samokhina M.A. Ocenka vozmozhnogo vyigrysha po masse pri ispol'zovanii dvigatelej maloj tyagi v ek-spedicii k Marsu [Estimation of the possible gain in mass when using low-thrust engines in an expedition to Mars]. SbornikdokladovXI Vserossijskogo s"ezdapo fundamental'nym problemam teoreticheskoj i prikladnoj mekhaniki [Collection of reports of the XI All-Russian Congress on Fundamental problems of Theoretical and Applied Mechanics], Kazan, 2015, pp. 1063-1065.

[3] Ishkov S.A., Starinova O.L. Optimization and modelling of movement with the solar sail. Izvestia of Samara Scientific Center of the Russian Academy of Sciences, 2005, vol. 7, no. 1, pp. 99-106.

[4] Lietmann G. Optimization techniques: with applications to aerospace systems. Academic Press, 1962, 453 p.

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

[6] Mozzhorina T.Yu. Numerical solution to problems of optimal control with switching by means of the shooting method. Mathematical modeling and Computational Methods, 2017, no. 2 (14), pp. 94-106.

[7] Mozzhorina T.Yu., Osipov V.V. Numerical solution of the problem of a soft landing by false position method. Innovative development, 2018, no. 8 (25), pp. 11-15.

[8] Fedorenko A.N. Matematicheskoe modelirovaniepereorientacii orbital'nogo kos-micheskogo apparata so sfericheskim solnechnym parusom [Mathematical modeling of the reorientation of an orbiting spacecraft with a spherical solar sail]. Diss. Cand. (Eng.) Sc. Moscow, 2014, 126 p.

[9] Chernyakina I.V. Programmy lokal'no-optimal'nogo upravleniya i traektorii geli-ocentricheskih pereletov kosmicheskogo apparata s solnechnym parusom s uchetom vozmushchenij [Programs of locally optimal control and trajectories of heliocentric flights of a spacecraft with a solar sail, taking into account disturbances]. Cand. (Eng.) Sc. Samara, 2020, 138 p.

[10] Starinova O.L., Gorbunova I.V. Optimizaciya geliocentricheskogo dvizheniya kosmicheskogo apparata s solnechnym parusom [Optimization of the heliocentric motion of a spacecraft with a solar sail]. Sbornik trudov XVII Vserossijskogo seminara po upravleniyu dvizheniem i navigacii letatel'nyh apparatov [Proceedings of the XVII All-Russian Seminar on Motion Control and Navigation of aircraft], Samara, 2015, pp. 168-171.

[11] Starinova O.L., Kurochkin D.V., Materova I.L. Optimal control choice of non-keplerian orbits with low-thrust propulsion. AIP Conference Proceedings, 2012, no. 1493 (1), pp. 964-971.

[12] Somov E.I., Butyrin S.A., Starinova O.L. Reguliruemyj solnechnyj parus v sisteme upravleniya dvizheniem minisputnika zemleobzora: modeli i ocenki effek-tivnosti [Adjustable solar sail in the motion control system of a mini-satellite Earth survey: models and efficiency estimates]. Sbornik trudov 7-oj Rossijskoj mul'tik-onferencii po problemam upravleniya « Upravlenie v morskih i aerokosmicheskih sistemah (UMAS-2014)» [Proceedings of the 7th Russian Multi-conference on Management Problems "Management in Marine and Aerospace Systems (UMAS-2014)"], St. Petersburg, 2014, pp. 545-556.

Mozzhorina T.Yu., Cand. Sc. (Eng.), Assoc. Professor of Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University. e-mail: [email protected]

Rakhmankulov D.A., Student of Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University. e-mail: [email protected]