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

УДК 519.6

DOI: 10.18698/2309-3684-2022-2-88101

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

© Т.Ю. Мозжорина, А.А. Закуражная МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

В данной работе рассматривается оптимизация перелета космического аппарата с Земной орбиты на орбиту Венеры с помощью ионных двигателей. Первый полет к планете состоялся в 1961 году советской автоматической межпланетной станцией «Венера-1», которая прошла в 100 000 километрах от Венеры. Кроме этого, в 1962 году был совершен полет американской станцией «Маринер-2». Самый последний корабль, запущенный к планете, был «Венера Экспресс» Европейского космического агентства в 2005 году, который долетел до Венеры за 153 дня. При решении текущей задачи были приняты следующие допущения: рассматривается межорбитальный перелет без учета притяжения планет, а орбиты планет считаются круговыми и лежащими в одной плоскости. В качестве управления был выбран угол между касательной скоростью космического аппарата и направлением тяги. Оптимизация управления проводилась с использованием принципа максимума Понтрягина. Полученная краевая задача для системы обыкновенных дифференциальных уравнений решалась численным методом — методом пристрелки. Для решения систем нелинейных алгебраических уравнений использовался метод Ньютона. Программа расчета была написана с использованием языка программирования С++. В результате работы удалось минимизировать время перелета между орбитами, таким образом была показана работоспособность метода пристрелки для решения задач оптимизации

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

Введение. Изучение Венеры всегда занимало важное место в программе планетных исследований Солнечной системы. Венера находится гораздо ближе к Солнцу, чем Марс, поэтому ей достается больше солнечной энергии, помимо этого на границе атмосферы этой планеты есть ионизированный слой, который создает достаточную защиту от солнечной радиации. Исследование этой планеты начались в 1961 году с полета советской автоматической межпланетной станции «Венера-1», которая пролетела в 100 тысячах километрах от планеты. Помимо серии «Венера» к планете направлялись такие аппараты как «Маринер» (Америка), «Магеллан» (Америка), «Венера-экспресс» (Европа), «Акацуки» (Япония). В 1970 году космический аппарат (КА) «Венера-7» впервые совершил на планету мягкую посадку, а в 1975 «Венера-9» передала на Землю панорамные изображения поверхности

Венеры [1]. Для получения изображения поверхности планеты в месте посадки КА была установлена панорамная камера в герметичном приборном отсеке, в котором были учтены суровые венерианские условия: давление до 100 атм. и температура до +500°С. Запечатлён-ные развалы камней могут быть результатом смещений в коре, что подтверждает факт тектонической активности на планете [2]. Несмотря на то, что сегодня ученые в основном занимаются изучением Марса, венерианская одиссея еще не завершена. В настоящий момент в разработке находится российская автоматическая межпланетная станция (АМС) «Венера-Д» («Венера-Долгоживущая»), которая является усовершенствованным аналогом советских аппаратов серии «Венера». Запуск аппарата запланирован на 2029 год для изучения структуры минералогического и химического состава поверхности Венеры, а также процессов взаимодействия атмосферы и других свойств планеты [3].

Ионный двигатель является одним из самых используемых для управления ориентацией и положением на орбите космических аппаратов. Такой двигатель может создать очень большой удельный импульс, но очень маленькую тягу, обусловленную массой выбрасываемых частиц. Его главным достоинством является время его работы, так как он может работать очень долго благодаря низкому потреблению газа, и все время своей работы будет уходить на разгон аппарата. Поэтому ионные двигатели используются только в космосе, где нет сопротивления воздуха и на небольших объектах вроде спутников. Принцип работы данного двигателя основан на создании реактивной тяги на базе ионизированного газа, разогнанного до высоких скоростей в электростатическом поле, и разгона космического корабля согласно третьему закону Ньютона [4].

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

Математическая постановка задачи. Рассматривается задача перелета космического аппарата с Земной орбиты на орбиту Венеры. При решении текущей задачи принимаем следующие допущения:

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

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

• абсолютное значение массы КА не учитывается.

В рамках этих допущений движение КА задается следующей системой уравнений:

= ^ - А (—)2 + Тотн МП(в), Ж Я ^ Я ' (1 - ботн^) ( ^

^ = - ™ + 008(0),

* Я (1 - Оотн<) (1)

йЯ

— = и,

йг

й< _ V

* ~ Я'

где г — время, с; и — скорость движения КА вдоль радиуса, м/с; V — скорость движения КА по касательной к окружности с центром в Солнце, м/с; Я — расстояние до Солнца, м; <р — полярный угол положения КА, рад; Я — начальное положение радиуса Земной орбиты, м; в — угол между касательной к окружности с центром в

Солнце и направлением тяги, рад; Д = — = ОМс, м/с2; ^ — сила

т0 ^

гравитационного притяжения в начальный момент времени, Н; т — масса КА, кг; О — гравитационная постоянная, м3 / с2 • кг; Мс — масса Солнца, кг; = Q / т0 — отношение расхода топлива двигателя к исходной массе, сутки-1; Тоотн = Т / т0 — отношение тяги

двигателя к исходной массе, м/с2.

На рис. 1 представлена схема перелета КА с орбиты Земли на орбиту Венеры с учетом допущения о круговых орбитах. Эксцентриситет планет достаточно близок к круговым. Так эксцентриситет Венеры равен 0,0067, а эксцентриситет Земли — 0,0167. В качестве управления при решении данной задачи был выбран угол в между касательной скоростью и направлением тяги.

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

г

3 = | йг = г ^ тт . (2)

о

Солнце

Рис. 1. Схема перелета между орбитами

Метод решения задачи. Задача ОУ решалась с использованием принципа максимума Понтрягина, который дает возможность среди всех допустимых процессов выделить те, которые могут претендовать на роль оптимальных [8]. Тогда функция Понтрягина для текущей задачи будет иметь вид:

Н = ¥ (_ ^ . (^ +-Тот-

' ( Я ( Я (1 _ От» • <)

Т

+ш (_ — +-т-

( Я (1 _ От. • О

вт(в)) +

(3)

у

К

где ¥и> У, ¥к > '

сопряженные переменные.

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

(

дв " 1 _ а

■ СОБ(в))

т.'

1_ а

-вт(в)) = 0.

Таким образом, получаем оптимальное управление:

tg (в*) =

¥и

'у

из которого выводятся следующие выражения:

'и ч '

Бт(в ) =

4

, СОБ(в)=

(4)

(5)

(6)

Данные знаки у выражений выбраны с учетом максимума функции Понтрягина:

д2 Н

Т

- (-щ 8т(в) - щ 008(в)) < 0 .

дв2 (1 - Qomн • г) Таким образом, получаем, что:

5/£и(8т(в)) = ),

Запишем сопряженную систему:

(7)

(8)

йШи V — — 1

йг Я

и — — 1

йг яш

йШя = Ши

йг

йш< = 0.

, йг

2у Я

¥и

V

Я2

Я2

Я3

Я

ш

Я2

V

< Я2

(9)

Пи-система (система уравнений движений и сопряженная система) будет иметь вид:

Т„

отн

*ч

008(в )

= - а • (^0)2 +

йг Я А (Я; (1 - Qomнг)

йу ш Т

— =--+-о

йг Я (1 - аотнг) йЯ

— = и, йг

й< _ V

йг ~ Я'

йХи V

— = —Ш йг Я

йш и

= я ш

йщ1

8т(в*),

(10)

Оя , 2v

Я

Ши ,

Я2

йг

Я V г, л Я

= ШЯ2 - 2 А •

По условию трансверсальности получаем, что щ (^) = 0. Поэтому

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

Краевые условия для данной задачи:

<

и0 = 0 ; у0 = 29,8; Д0 = 149,6 • 10б; р0 = 0 ; и1 = 0; у1 = 35 ; Д = 108,2• 106.

Из условий трансверсальности: Н= 0, ) = 0.

Алгоритм численного решения задачи. Алгоритм численного решения данной задачи состоит из нескольких этапов: во внутреннем цикле проводится интегрирование системы обыкновенных дифференциальных уравнений методом Рунге-Кутта четвертого порядка [9], при котором образуются так называемые невязки — расхождения решения на правом конце. Далее во внешнем цикле применяется модифицированный метод Ньютона для решения системы нелинейных алгебраических уравнений невязок [10]. Для улучшения сходимости метода Ньютона использовалась локальная норма, введенная Федоренко Р.П. в [11]. Полученная СЛАУ решалась с использованием Ш-разложения [12]. Схема алгоритма решения данной задачи представлена на рис. 2.

Начало

т

Ввод началь данных

iteration< max iteration

ных

Вычисление невязок

Метод Рунге-Кутта

Вычисление производных

< S

Шаг по параметрам пристрелки

Рис. 2. Блок-схема численного алгоритма

Результаты численного моделирования. Для расчета поставленной задачи была написана программа на языке программирования С++. Были использованы следующие исходные данные: • начальные условия (известные значения для Земли):

2

3

„ км ___ км

м0 = 0 — ; ^=29,8 — ; с с

Я =149,6• 10бкм; р0 =0 рад ;

• граничные условия (известные значения для Венеры):

„ км км

и = 0— ; V = 35— ; с с

Я =108,2 • 106 км;

• постоянные:

Т^ =8,3 • 10-4 4;

с

Qomн =1,29• 10 3сутки 1;

у • М м

А =уМ^ = 0,00593-^.

А Я2 , с2

Параметрами пристрелки будут недостающие начальные условия:

Ш(О = Ш(О = -852; Шя(О = -10-5;

г = 150 — время перелета.

Уравнения невязок:

=и(гх) ^ 0, v(г1)- V!

б2 = 0, V

¿ =яго-я ^ 0,

3 я ,

¿4 =н (г) ^ 0.

(11)

Точность решения в методе Ньютона была выбрана равной £ = 10 6. Начальный угол в = -45°. Для численного решения системы ДУ методом Рунге-Кутта было выбрано 1000 шагов. Таким образом, для нахождения решения потребовалось 10 итераций метода Ньютона. По итогу параметры пристрелки получились равными:

Ши (0) = -62 1,04 ; щ (0) = -9 1 1,0б;

О (0) = -0,0022 ; г =139,9 .

На рис. 3-4 изображены графики изменения фазовых переменных, полученных при решении данной задачи.

и, м/с 0 -2000 -4000 -6000

0 25 50 75 100 125 сутки

а

V, м/с 36000

34000

32000

30000

0 25 50 75 100 125 сутки

б

Рис. 3. График изменения скорости космического аппарата: а — радиальная скорость; б — касательная скорость

Л, м 1,5 -1011 1,4 -1011 1,3 -1011 1,2 -1011 1,1-1011

р, градусы 175 150 125 100 75 50 25 0

0 25 50 75 100 125 сутки

а

100 125 /, сутки

0 25 50 75 б

Рис. 4. График изменения: а — радиуса от Солнца; б — полярного угла

На рис. 5-6 изображены графики изменения сопряженных переменных.

¥и 400 200

-200

-400

-600

125 /, сутки

0

-200

-400

-600

-800

50 75

б

100 125 г, сутки

Рис. 5. График изменения: а — сопряженной переменной ц/и; б — сопряженной переменной ^

ц/к -0,000100 -0,000125 -0,000150 -0,000175 -0,000200 -0,000225 -0,000250

0 25 50 75 100 125 /, сутки Рис. 6. График изменения сопряженной переменной

0

а

Для проверки оптимальности решения воспользуемся тем, что функция Понтрягина в конце процесса равна нулю, а также частная производная от функции Понтрягина равна ее полной производной:

ТТ, Л Л дН жн

Н (г1) = 0, -=-.

дг Жг

(12)

Частную производную можно вычислить аналитически по формуле:

дН Т • 0

_ отн х--отн

-(щи вт^ + щ С08(^)).

(13)

дг (1 - О-отн • г)2

Полную производную вычисляем численно, используя левую разностную схему:

ЖН Н - Н

Жг

Аг

(14)

На рис. 7а представлен график изменения функции Понтрягина в оптимальном процессе, на рис. 7 б — совместный график производных функции Понтрягина. Изменение управления в течение всего процесса изображено на рис. 8.

Н 0 -0,02 -0,04 -0,06 -0,08

Производные 0,00150 0,00125 0,00100 0,00075 0,00050 0,00025

0 25 5 0 75 100 125 г, сутки

0

25 50 75 100 125 и сутки б

Рис. 7. График изменения: а — функции Понтрягина; б — производных функции Понтрягина (--полная производная;----частная производная)

а

в, градусы 150 100 50 0

-50 -100

0 25 50 75 100 125 сутки Рис. 8. График изменения управления

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

У-1011, м 1,5 1,0 0,5 0

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

-1,5 -1,0 -0,5 0 0,5 1,0 1,5 х-1011, м

Рис. 9. Траектория перелета космического аппарата с орбиты Земли на орбиту Венеры:

----траектория спутника,--орбита Земли;

--орбита Венеры;--Солнце

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

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

ЛИТЕРАТУРА

[1] Моклецов А. Исследование планеты Венера космическими аппаратами. РИА Новости [Электронный ресурс], 2015. URL: https://ria.ru/20150611/1069063562.html (дата обращения: 01.03.2022)

[2] Fuglesang C., Zetterling C.-M., Wilson C.F. Venus long-life surface package (VL2SP). Proceedings of the International Astronautical Congress, IAC, 2017, vol. 5, pp. 3035-3043.

[3] Российская Газета. Когда состоится полет на Венеру? [Электронный ресурс]. URL: https://rg.ru/2018/10/09/kogda-sostoitsia-pervyj-polet-na-veneru.html (дата обращения: 10.02.2022).

[4] Hi - News. Как работает ионный двигатель и где он применяется [Электронный ресурс]. URL: https://hi-news.ru/eto-interesno/kak-rabotaet-ionnyj-dvigatel-i-gde-on-primenyaetsya.html (дата обращения: 12.02.2022)

[5] Моршнева И. В., Овчинникова С. Н. Численное решение краевых задач для обыкновенных дифференциальных уравнений. Метод стрельбы. Методические указания для студентов 3 и 4 курсов мехмата. Ростов-на-Дону, УПЛ РГУ, 2003, 29 с.

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

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

[8] Ванько В.И., Ермошина О.В., Кувыркин Г.Н. Вариационное исчисление и оптимальное управление: учеб. для вузов. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2006, 488 с.

[9] Adesola O.A., Samson A.A., Ayomide O.A., Adekunle O.A. Dynamic computation of runge-kutta's fourth-order algorithm for first and second order ordinary. International Journal of Computer Science Issues, 2015, vol. 12, iss. 3, pp. 211-218.

[10] Иванов А.П. Метод Ньютона решения нелинейных уравнений и систем уравнений. Санкт-Петербург, Изд-во СПбГУ, 2007, 9 с.

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

[12] Семушин И.В. Вычислительные методы алгебры и оценивания: учебное пособие. Ульяновск, Изд-во УлГТУ, 2011, 366 с.

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

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

Мозжорина Т.Ю., Закуражная Д.А. Моделирование и оптимизация управления полетом космического аппарата с орбиты Земли на орбиту Венеры с помощью ионных двигателей. Математическое моделирование и численные методы, 2022, № 2, с. 88-101.

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

Закуражная Алина Андреевна — студент кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. e-mail: azakurazhnaya@mail.ru

Simulation and optimization of spacecraft flight control from Earth orbit to Venus orbit using ion thrusters

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

In this paper, optimization of the flight of a low-mass satellite from Earth orbit to the orbit of Venus using ion engines is considered. The first flight to the planet took place in 1961 by the Soviet automatic interplanetary station "Venus-1", which passed 100,000 kilometers from Venus. In addition, in 1962, the American station "Mariner-2" was flown. The most recent spacecraft launched to the planet was the European Space Agency's Venus Express in 2005, which flew to Venus in 153 days. When solving the current problem, the following assumptions were made: an interorbital flight is considered without taking into account the attraction of the planets, and the orbits of the planets are considered circular and lying in the same plane. The angle between the tangential velocity of the spacecraft and the thrust direction was chosen as the control. Optimization of satellite control was carried out using the Pontryagin maximum principle. The resulting boundary value problem for a system of ordinary differential equations was solved by a numerical method — the targeting method. Newton's method was used to solve systems of nonlinear algebraic equations. The calculation program was written using the C+ + programming language. As a result of the work, it was possible to minimize the flight time between orbits, thus the operability of the shooting method for solving optimization problems was shown.

Keywords: optimal control, flight between the orbits of Earth and Venus, shooting method, ion thrusters, Pontryagin maximum principle, numerical methods, boundary value problems of ordinary differential equations

REFERENCES

[1] Mokletsov A. Issledovanie planety Venera kosmicheskimi apparatami [Exploration of the planet Venus by spacecraft]. RIA Novosti [RIA News] [Electronic resource], 2019. URL: https://ria.ru/20150611/1069063562.html (accessed: 01.03.2022)

[2] Fuglesang C., Zetterling C.-M., Wilson C.F. Venus long-life surface package (VL2SP). Proceedings of the International Astronautical Congress, IAC, 2017, vol. 5, pp. 3035-3043.

[3] Rossiyskaya Gazeta. Kogda sostoitsa polet na Veneru [When will the flight to Venus take place?] [Electronic resource]. URL: https://rg.ru/2018/10/09/kogda-sostoitsia-pervyj-polet-na-veneru.html (accessed: 10.02.2022)

[4] Hi - News. Kak rabotaet ionnyj dvigatel' i gde on primenyaetsya? [How the ion engine works and where it is applied?] [Electronic resource]. URL: https://hi-news.ru/eto-interesno/kak-rabotaet-ionnyj-dvigatel-i-gde-on-primenyaetsya.html (accessed: 12.02.2022).

[5] Morshneva I.V., Ovchinnikova S.N. Chislennoe reshenie kraevyh zadach dlya obyknovennyh differencial'nyh uravnenij. Metod strel'by. Metodicheskie ukaza-niya dlya studentov 3 i 4 kursov mekhmata [Numerical solution of edge problems for ordinary differential equations. The method of shooting. Methodological guidelines for students of the 3rd and 4th courses of the Faculty of Mechanics]. Rostov-on-Don, UPL RSU Publ., 2003, 29 p.

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

[7] 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, pp. 94-106.

[8] Vanko V.I., Ermoshina O.V., Kuvyrkin G.N. Variacionnoe ischislenie i optimal'noe upravlenie: ucheb. dlya vuzov [Calculus of variations and optimal control: Studies for universities]. Moscow, BMSTU Publ., 2006, 488 p.

[9] Adesola O.A., Samson A.A., Ayomide O.A., Adekunle O.A. Dynamic computation of runge-kutta's fourth-order algorithm for first and second order ordinary. International Journal of Computer Science Issues, 2015, vol. 12, iss. 3, pp. 211-218.

[10] Ivanov A.P. Metod N'yutona resheniya nelinejnyh uravnenij i sistem uravnenij [Newton's method for solving nonlinear equations and systems of equations]. St. Petersburg, SPBU Publ., 2007, 9 p.

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

[12] Semushin I.V. Vychislitel'nye metody algebry i ocenivaniya: uchebnoe posobie [Computational methods of algebra and evaluation: textbook]. Ulyanovsk, UlSTU Publ., 2011, 366 p.

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

Zakurazhnaya A.A., Student of Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University. e-mail: azakurazhnaya@mail.ru