Размещено на http://www.allbest.ru/

63

Размещено на http://www.allbest.ru/

Содержание

• Введение
• 1. Проблема удержания космического аппарата на ограниченной орбите вокруг точки L₂ системы Солнце-Земля
• 2. Математическая модель
• 3. Стратегия удержания КА на гало-орбите вокруг точки L₂ системы Солнце-Земля
• 3.1 Описание стратегии удержания КА
• 3.2 Реализация стратегии удержания КА
• 3.2.1 Алгоритм подбора начальной скорости и величины корректирующего импульса
• 3.2.2 Моделирование технических ограничений
• 3.2.3 Сценарий, моделирующий движение КА на гало-орбите с периодическим применением корректирующих импульсов
• 4. Расчет направления неустойчивости
• 4.1 Методика расчета направления устойчивости
• 4.2 Результаты расчета направлений устойчивости и неустойчивости
• 4.3 Интерполяция направления неустойчивости
• 4.4 Зависимость направления неустойчивости от координаты Z
• 5. Применение разработанных инструментов к моделированию движения КА на гало-орбите
• 5.1 Исследование зависимости энергетики поддержания гало-орбиты от места и направления исполнения импульса
• 5.2 Исследование влияния неточности определения параметров КА на геометрию гало-орбиты
• 5.2.1 Исследование случая неточного определения скорости КА
• 5.2.2 Исследование случая неточного определения координат КА
• 5.3 Имитационное моделирование движения КА на гало-орбите с учетом направления неустойчивости
• Заключение
• Список использованных источников


Введение

Точками либрации в ограниченной задаче трех тел, описывающей движение тела малой массы в гравитационном поле, создаваемом двумя массивными телами, вращающимися вокруг общего барицентра, называются точки равновесия тела с малой массой (космического аппарата) относительно более массивных тел (Солнца и Земли). Существуют 5 точек либрации. 3 из них называются коллинеарными и расположены на прямой, соединяющей массивные тела. 2 точки расположены на орбите одного из массивных тел (Земли). Точка L₂ расположена на прямой, соединяющей Солнце и Землю, за Землей.

Рис.1. Точки либрации.

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

Гало-орбиты впервые были предложены Робертом Фаркуаром для программы Аполлон в его диссертации в 1968 г. [1]. Он предложил поместить КА на гало-орбиту вокруг точки L₂ системы Земля-Луна для связи с аппаратом, находящимся на обратной стороне Луны. Хотя данная концепция не была реализована, впоследствии были запущены многие миссии, включающие в себя гало-орбиты. Наиболее известными из них являются:

· ISEE-3, ACE - запущены на гало-орбиты вокруг точки L₁ системы Солнце-Земля;

· Herschel, Plank - запущены на гало-орбиты вокруг точки L₂ системы Солнце-Земля.

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

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

Комплекс корректирующих импульсов с учетом места их исполнения, периодичности и т.д. называется стратегией удержания КА на гало-орбите (station-keeping). Так же под стратегией понимается непосредственно методика расчета импульсов. Разработка стратегии удержания является одной из важнейших задач при планировании миссии.

Данная работа посвящена разработке и анализу стратегии удержания КА на гало-орбите вокруг точки L₂ системы Солнце-Земля.

Целью данной работы является разработка и реализация стратегии удержания КА и моделирование движения КА на гало-орбите.

Задачами данной работы являются:

· Создание программного обеспечения (ПО) для расчета движения КА на гало-орбите;

· Моделирование движения КА на гало-орбитах с различными амплитудами с использованием разработанного ПО;

· Получение зависимости между энергетикой поддержания орбиты и направлением исполнения импульса для различных амплитуд орбит и анализ полученных зависимостей;

· Расчет направления неустойчивости и применение полученных результатов при имитационном моделировании движения КА на гало-орбите.

Результаты данной работы, в частности, разработанное ПО могут быть применены для моделирования как гало-орбит, так и других ограниченных орбит вокруг точки L₂ системы Солнце-Земля.



1. Современные методики удержания космического аппарата на ограниченной орбите вокруг точки L₂ системы Солнце-Земля

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

Первоначальным этапом планирования миссии к точке либрации является расчет т.н. номинальной орбиты КА. Как правило, для описания движения КА по номинальной орбите используются упрощенные модели. Так, в статьях [1]-[11] приводится ограниченная круговая задача трех тел. В статьях [1]-[3] приводятся уравнения движения КА в инерциальной системе координат; в статьях [1]-[7] также приводятся уравнения движения КА в системе координат, связанной с соответствующей точкой либрации. На начальном этапе исследования наибольший интерес вызывает линеаризованная система уравнений движения; такие уравнения приведены в статьях [1], [8]-[10]. Более подробно ограниченная задача трех тел и уравнения движения КА на гало-орбите будут рассмотрены в соответствующем разделе.

В [11] приводится классификация стратегий удержания КА на гало-орбите вокруг точки L₁ системы Солнце-Земля. Стратегии коррекции разделяются на два класса: loose control strategy (техника свободного контроля), в которой варьируется только одна компонента импульса с целью приближения орбиты КА к номинальной и tight control technique (техника строгого контроля), подразумевающая варьирование двух и более компонент корректирующего импульса.

В качестве примера tight control technique в [11] приводится методика расчета импульсов коррекции для КА ISEE-3. На первом этапе обрабатывались данные о траектории КА за время, прошедшее с момента выполнения последнего маневра. Затем орбита аппарата численно интегрировалась на один оборот вперед (примерно 178 дней) в реальной модели сил и в 8 равномерно распределенных точках вычислялось расстояние между номинальной орбитой и полученной. На следующем этапе сумма полученных расстояний минимизировалась за счет варьирования компонент вектора импульса. Эта процедура производилась для различных дат, после чего для исполнения импульса выбиралась дата, которой соответствовали минимальные затраты топлива (хотя иногда, по техническим причинам, выбирался день с неоптимальной величиной импульса).

Всего в ходе миссии было выполнено 15 импульсов коррекции. Импульсы совершались в среднем один раз в 82 дня, а суммарный импульс составил 30 м/с. Согласно [11], суммарный импульс мог бы быть меньше, но т.к. ISEE-3 был первым КА, выведенным на гало-орбиту, ученые, управлявшие его движением, стремились сделать минимальное количество маневров с минимальным риском.

В качестве примера loose control strategy в [11] описана методика, с помощью которой рассчитывались корректирующие импульсы для КА SOHO. Использованная методика называется orbital energy balancing. Ее суть заключается в том, что если энергия аппарата слишком велика (т.е. был исполнен слишком большой импульс), он отклоняется от номинальной орбиты в одну сторону, а если слишком мала - в другую. В данном случае уравнения движения КА интегрировались до оси Солнце-Земля. В процессе управления КА в течение миссии было добавлено условие, чтобы аппарат двигался перпендикулярно данной оси при ее пересечении. Импульсы исполнялись только вдоль прямой Солнце-Земля (в отрицательном или положительном направлении). Единственным варьируемым параметром была величина импульса коррекции (с учетом знака). В результате описанная аппаратом траектория была не периодической гало-орбитой, а квазигало-орбитой. С момента выхода КА на орбиту в марте 1996 было исполнено 9 корректирующих импульсов, не считая экстренные маневры, выполненные в период с сентября 1998 по март 1999 (в тот момент аппарат сошел с гало-орбиты и экстренные маневры выполнялись для его возвращения). Суммарный импульс, без учета экстренных маневров, составил чуть более 5 м/с. Также стратегия удержания КА SOHO была описана в статье [12].

В [13] приведен аналитический подход к расчету величин корректирующих импульсов. В статье приведены 2 основные методики: target point approach (приближение целевой точкой) и Floquet mode approach (приближение модами Флоке). Несмотря на то что статья посвящена орбитам вокруг точки L₂ системы Земля-Луна, данные методики применимы и для орбит вокруг точек либрации системы Солнце-Земля.

Методика target point approach, впервые описанная в [15], использует штрафную функцию, которая зависит от затрат топлива и от отклонения параметров КА от номинальных. Задача метода состоит в минимизации штрафной функции. Для этой методики вводится ряд ограничений. Во-первых, между импульсами должен быть перерыв не менее некоторого заданного времени. Во-вторых, отклонение от номинальной траектории также должно быть не меньше заданного. В-третьих, при моделировании процесса удержания КА, величины отклонений от номинальной траектории сравниваются между собой и если отклонение уменьшается, величина импульса не рассчитывается. Для успешного расчета величины импульса требуется выполнение всех трех условий.

Вторая методика, Floquet mode approach, впервые предложенная в [16], описана в [13] в случае ограниченной круговой задачи трех тел. Для уравнений движения КА выводится матрица монодромии, собственным числам и векторам которой соответствуют различные характеристики. Наиболее важной характеристикой, с точки зрения управления КА, является направление неустойчивости. Оно совпадает с первым собственным вектором матрицы монодромии. Импульсы, совершенные в направлении неустойчивости, являются наиболее эффективными. Соответственно, расчет величин импульса коррекции в дальнейшем ведется в предположении, что импульсы совершаются вдоль направления неустойчивости. Как и в предыдущем случае, для данной методики вводятся ограничения на минимальный интервал между импульсами и минимальное отклонение от номинальной орбиты.

Приведенная в [14] методика была разработана для миссии ARTEMIS. В данной статье также рассчитываются направления устойчивости и неустойчивости с использованием матрицы монодромии. Для удержания КА на орбите используется т.н. optimal continuation strategy (стратегия оптимального продолжения). Она разработана для удержания КА на гало-орбите вокруг точки L₂ системы Земля-Луна на протяжении 1-2 оборотов. Это достигается благодаря подбору такой величины импульса, при которой выполняются заданные оператором условия (например, определенная координата и скорость по одной из осей) при достижении аппаратом плоскости орбиты Луны.

При разработке стратегии удержания КА на гало-орбите важно знать направление неустойчивости, т.к. импульсы в этом направлении являются наиболее эффективными. В литературе предлагается методика расчета направления неустойчивости, связанная с матрицей монодромии уравнений движения и модами Флоке. Данный метод был описан в [17]-[19].



2. Стратегия удержания КА на гало-орбите вокруг точки L₂ системы Солнце-Земля

2.1 Математическая модель

Для описания движения КА по ограниченной орбите введем вращающуюся систему координат, связанную с точкой L₂. Центр системы координат расположен в точке L₂, ось X совпадает с осью Солнце-Земля и направлена от Солнца к Земле, ось Z направлена в северный полюс эклиптики, ось Y дополняет систему координат до правой тройки. Иллюстрация системы координат приведена на рис. 2. Далее в работе всегда подразумевается данная система координат, если не оговорено иное.

Рис. 2. Вращающаяся система координат, связанная с точкой L2.

Движение КА в системе n тел в инерциальной системе координат описывается следующими уравнениями:

(1)

где n - количество притягивающих центров, G - гравитационная постоянная, R - радиус-вектор КА, m_i - масса i-го тела, R_i - радиус-вектор i-го тела.

Во вращающейся системе координат, при переходе к ограниченной задаче трех тел () уравнения (1) могут быть приведены к виду [1]:

(2)

где c₂ - параметр, зависящий от масс тел, a_x, a_y, a_z - возмущающие ускорения, которые являются функциями координат КА и эксцентриситета орбиты КА.

Уравнения (2) допускают линеаризацию в окрестности точки L₂. Тогда они принимают следующий вид [2]:

(3)

Система (3) имеет следующее решение:

(4)



где , , , , , - фаза колебаний в плоскости XY, - фаза колебаний по оси Z. Коэффициенты , , , и фазы , зависят от начальных условий.

Решения и представляют собой линейную комбинацию трех компонент: ограниченной (и ), экспоненциально возрастающей и убывающей к нулю . Будем предполагать, что и решение системы (1), записанной во вращающейся системе координат, связанной с точкой L₂, также представимо в виде линейной комбинации ограниченной, возрастающей и убывающей компонент:

где , - возрастающие компоненты, , - убывающие компоненты, , , - ограниченные компоненты.

Коэффициенты и могут быть как положительными, так и отрицательными. Поэтому термины «возрастающая» и «убывающая» компоненты подразумевают возрастание и убывание по модулю.

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

Ограниченные орбиты, периоды обращения КА по которым по осям Y и Z совпадают, называются гало-орбитами. На рис. 3-5 представлены проекции движения КА на гало-орбите на плоскости XY, YZ, XZ.



Рис. 3. Проекция движения КА на гало-орбите на плоскость XY.

Рис. 4. Проекция движения КА на гало-орбите на плоскость YZ.



Рис. 5. Проекция движения КА на гало-орбите на плоскость XZ.

Из рис. 3-5 можно видеть, что гало-орбиты симметричны относительно плоскости XZ и несимметричны относительно плоскости эклиптики и плоскости YZ. В зависимости от того, в положительном или в отрицательном направлении оси Z амплитуда гало-орбиты является большей, гало-орбиты называют северными и южными. Их изображение в проекции на плоскость XZ представлено на рис. 6. Гало-орбита, изображенная на рис. 3-5, является южной. Амплитуды гало-орбиты определяются ее начальными координатами [22].



Рис. 6. Проекция движения КА по северной и южной гало-орбитам на плоскость XZ.

2.2 Описание стратегии удержания КА

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

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

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

Этот факт используется следующим образом. Предположим, что в начальный момент времени КА находится на плоскости XZ (Y=0) и движется перпендикулярно ей (V_x=V_z=0 км/с). Положим также значение скорости вдоль оси Y равным некоторому значению V_y0. Построим виртуальные плоскости X=X_min и X=X_max такие, что гало-орбита лежит между ними, и будем интегрировать уравнения движения КА до того, как КА достигнет любой из указанных плоскостей. Если аппарат достиг левой границы (X=X_min), то значение коэффициента при возрастающей компоненте отрицательно, и наоборот. Т.о. определяется функция - конечная координата X КА от начальной скорости V_y0. Эта функция терпит разрыв, т.е. КА не достигает ни левой, ни правой границы, при некотором значении V_y0. Задача заключается в отыскании этого значения. Данная задача решается методом деления отрезка пополам. Работа алгоритма заканчивается при достижении минимальной вычислительной погрешности, обеспечиваемой средой вычислений . Иллюстрация работы алгоритма представлена на рис. 7.



Рис. 7. Подбор начальной скорости КА.

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

Разработанный алгоритм позволяет рассчитывать начальную скорость КА и величины корректирующих импульсов для любых ограниченных орбит в окрестности точки L₂ системы Солнце-Земля. При этом он позволяет рассчитывать импульсы коррекции в любой точке орбиты в любом направлении.



2.3 Реализация стратегии удержания КА

2.3.1 Алгоритм подбора начальной скорости и величины корректирующего импульса

Описанные алгоритмы были реализованы в программе GMAT (General Mission Analysis Tool). Данная программа позволяет интегрировать уравнения движения КА в реалистичной модели сил различными численными методами, в т.ч. методами Рунге-Кутта различных порядков, Принса-Дорманда и др. В данной работе численное интегрирование проводилось методом Рунге-Кутта 8-9 порядка. В GMAT были созданы сценарии, в которых производился подбор начальной скорости КА и величин корректирующих импульсов и интегрировались уравнения движения КА. Метод деления отрезка пополам позволяет достичь любой наперед заданной точности вплоть до машинного ограничения 10^-16.

Например, подберем начальную скорость для КА, имеющего следующие начальные параметры:

· X = -277548 км;

· Y = 0 км;

· Z = 200000 км;

· V_x = V_z = 0 км/с.

Алгоритмом была найдена скорость V_y, равная -0.372794445417389 км/с. Она обеспечивает наиболее продолжительное время пребывания КА на гало-орбите - 4 оборота (около 700 суток). Для более длительного пребывания КА на гало-орбите требуется исполнение математических коррекций.

Блок-схемы алгоритмов подбора величины начальной скорости и величины корректирующего импульса представлены на рис. 8 и 9.



Рис. 8. Алгоритм подбора начальной скорости КА.

Из рис. 9 можно видеть, что алгоритмы подбора начальной скорости КА и величины импульса коррекции в целом совпадают. Единственным отличием алгоритма подбора величины импульса является учет направления, в котором исполняется импульс. Для этого вводятся переменные DirX и DirY, являющиеся направляющими и орта вектора коррекции, положенного в плоскость XY.



Рис. 9. Алгоритм подбора величины импульса коррекции.

2.3.2 Моделирование технических ограничений

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

Для моделирования отклонения параметров от номинальных был использован метод Монте-Карло, реализованный в пакете Matlab. Были созданы две функции, генерирующие случайные отклонения параметров от входных данных в шаре заданного радиуса или на сфере заданного радиуса.

Входными параметрами функции MatLab являются вектор, нуждающийся в рассеянии, и радиус сферы (шара) рассеяния.

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

2.3.3 Сценарий, моделирующий движение КА на гало-орбите с периодическим применением корректирующих импульсов

Для моделирования движения КА на гало-орбите был разработан сценарий в пакете GMAT. Он позволяет моделировать движение КА по ограниченной орбите с периодическим применением корректирующих импульсов (например, в определенной точке орбиты, либо раз в несколько дней и пр.) в заданном направлении. Ниже приведено описание сценария, моделирующего движения КА на гало-орбите, причем место совершения коррекции циклически изменяется, а направление совершения коррекции - постоянно. Блок-схема сценария представлена на рис. 10.

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



Рис. 10. Сценарий, моделирующий движение КА на гало-орбите с периодическим применением корректирующих импульсов.



В начале сценария задаются начальные координаты КА и количество импульсов коррекции. Т.к. в данном случае коррекции совершаются один раз в оборот, количество коррекций совпадает с числом оборотов КА на орбите. Также указывается направление совершения коррекций (Beta), которое задается как угол между осью X и направлением совершения коррекции. Место исполнения импульса (Alpha) задается как угол между осью X и радиус-вектором аппарата и меняется в цикле от 180є до -170є с шагом -10є. Для хранения исходных параметров математической модели КА, используемой в цикле, параметры модели копируются во вспомогательную модель КА.

В работе сценария вызывается алгоритм подбора величины импульса коррекции. На первом шаге этот алгоритм работает как алгоритм подбора скорости за счет того, что в разложении направления импульса по осям компоненты X и Z полагаются равными 0, а компонента Y = 1. На следующих шагах производится подбор величины импульса коррекции в направлении Beta.

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

При исполнении заданного количества корректирующих импульсов происходит выход из цикла и расчет статистических характеристик (среднего значения и дисперсии импульса и времени полета, максимального и минимального значения импульса).

По достижении переменной Alpha границы цикла (-170°) производится выход из цикла и работа сценария завершается.

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

2.4 Методика расчета направления неустойчивости

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

В данном исследовании была разработана методика расчета направления устойчивости. Т.к. известно, что направления устойчивости и неустойчивости ортогональны друг другу, зная направление устойчивости, можно отыскать и направление неустойчивости. Методика заключается в следующем: в некоторой точке гало-орбиты производится возмущение движения КА, т.е. скорость аппарата изменяется по некоторому направлению. Задача состоит в отыскании такого направления изменения скорости, при котором КА как можно дольше находится на орбите.

Как и в случае отыскания начальной скорости КА, данная задача решается итерационно. При каждом направлении изменения скорости аппарат отклоняется либо в сторону положительного, либо в сторону отрицательного изменения координаты X, т.е. возникает зависимость направления изменения скорости от конечной координаты X аппарата. Эта функция имеет разрыв, при котором КА не отклоняется ни в одну из сторон, т.е. остается на ограниченной орбите. Требуется отыскать этот разрыв.

Данный алгоритм был реализован в пакете GMAT. Блок-схема алгоритма приведена на рис. 1.

Рис. 11. Алгоритм расчета направления устойчивости.

Таким образом, были разработаны следующие алгоритмы:

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

· Алгоритм подбора величин корректирующих импульсов. Данный алгоритм основан на алгоритме подбора начальной скорости; основное отличие - учет направления исполнения импульса.

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

Также были в пакете Matlab были разработаны функции, моделирующие технические ограничения на точность определения параметров КА и исполнение корректирующих импульсов.

На основе разработанных алгоритмов в пакете GMAT был создан сценарий, позволяющий моделировать движение КА на ограниченных орбитах вокруг точки L₂.

1. 

3 Применение разработанных инструментов к моделированию движения КА на гало-орбите

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

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

· Угол б - определяет место исполнения импульса (см. рис. 12).

· Угол в - определяет направление исполнения импульса (угол между проекцией направления маневра на плоскость эклиптики и осью Солнце-Земля). Данный угол задает ось, вдоль которой совершается маневр в положительном или отрицательном направлении.

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

Рис. 12. Параметры б и в.

Задачей данного раздела является оценить зависимость суммарного импульса от места и направления совершения коррекции. Данное исследование было проведено без учета погрешности определения параметров КА и выдачи импульса коррекции. Коррекция совершается один раз в оборот в точке, описываемой углом б, в направлении, описываемом углом в.

Были исследованы 3 гало-орбиты со следующими начальными координатами:

· X = -277548 км, Y = 0 км, Z = 200000 км;

· X = -373454 км, Y = 0 км, Z = 400000 км;

· X = -566256 км, Y = 0 км, Z = 600000 км.

Проекции движения КА, полученные в результате моделирования движения КА на данных орбитах на протяжении 50 оборотов, представлены на рис. 13-15.

Рис. 13. Проекция движения КА на гало-орбитах с различными амплитудами на плоскость XY.

Рис. 14. Проекция движения КА на гало-орбитах с различными амплитудами на плоскость YZ.

Рис. 15. Проекция движения КА на гало-орбитах с различными амплитудами на плоскость ZX.

На рис. 16 представлена зависимость среднего импульса от места исполнения коррекций (угол б). Для каждого угла б была рассчитана траектория движения КА на гало-орбите на протяжении 350 оборотов. Направление импульса в данном случае совпадает с направлением оси X, т.е. направлением от Солнца к Земле (угол в = 0є).

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

Из рис. 16 видно, что в целом характер зависимости совпадает для различных орбит. Возрастание данных зависимостей при б меньших -100 ° и больших 70° объясняется тем, что в этих точках направление исполнения импульса близко к направлению устойчивости. В связи с этим требуются большие затраты на удержание КА на гало-орбите.

Также был рассмотрен более общий случай задачи нескольких тел - при расчетах учитывалось не только гравитационное поле Солнца и Земли, но и остальных планет Солнечной системы и Луны. На рис. 17 представлены полученные зависимости.

Рис. 17. Среднее значение импульса в зависимости от места его исполнения.

Поскольку в реальности положение и скорость аппарата невозможно определить с бесконечной точностью, имеет смысл проанализировать, каким образом малое изменение скорости КА влияет на зависимость математического ожидания значения импульса коррекции от места исполнения импульса. Однако, так как исследование данной зависимости требует, чтобы аппарат удерживался на гало-орбите как минимум 1 оборот после исполнения коррекции, использование неопределенностей, возникающих на практике, для данного исследования невозможно. Это связано с тем, что из-за значительных отклонений параметров от номинальных не всегда удается рассчитать величину корректирующего импульса так, чтобы КА оставался на гало-орбите на протяжении 1 оборота. Чтобы качественно оценить влияние неопределенности знания скорости КА на математическое ожидание значения импульса поддержания, было использовано отклонение 10^-9 км/с. Полученная зависимость представлена на рис. 18.

Рис. 18. Среднее значение импульса в зависимости от места его исполнения для случая неточного определения скорости.

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

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

На рис. 19-27 представлены графики зависимости средней величины импульса от направления исполнения для гало-орбиты с начальной координатой по оси Z равной 200000 км.

Рис. 19. Зависимость среднего значения импульса от угла в, 150є?б?180є.

Рис. 20. Зависимость среднего значения импульса от угла в, 110є?б?140є.



Рис. 21. Зависимость среднего значения импульса от угла в, 70є?б?100є.

Рис. 22. Зависимость среднего значения импульса от угла в, 30є?б?70є.



Рис. 23. Зависимость среднего значения импульса от угла в, -10є?б?20є.

Рис. 24. Зависимость среднего значения импульса от угла в, -50є?б?-20є.

космический аппарат орбита



Рис. 25. Зависимость среднего значения импульса от угла в, -90є?б?-60є.

Рис. 26. Зависимость среднего значения импульса от угла в, -130є?б?-100є.



Рис. 27. Зависимость среднего значения импульса от угла в, -170є?б?-140є.

Аналогичные зависимости были рассчитаны для гало-орбит с начальными координатами Z = 400000 км и Z = 600000 км. Их характер совпадает с характером зависимостей, полученных для гало-орбиты с начальной координатой Z = 200000 км.

Т.к. угол в задает ось, вдоль которой совершается импульс (как в положительном, так и в отрицательном направлении), данные графики являются р-периодическими. Это позволяет рассчитывать зависимость средней величины импульса от направления исполнения импульса только для одного отрезка по углу в (в данном случае - для отрезка, содержащего 0є).

Графики имеют ярко выраженные минимумы. Данные минимумы соответствуют направлениям неустойчивости орбиты в точке, соответствующей углу б.

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

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

На рис. 28-30 представлены графики зависимости направления, соответствующего наименьшему среднему импульсу, от места совершения маневра для различных орбит (начальные координаты по Z 200000 км, 400000 км, 600000 км).

Рис. 28. Зависимость угла в от угла б, соответствующего наименьшему среднему значению импульса и направлению неустойчивости, z = 200000 км.



Рис. 29. Зависимость угла в от угла б, соответствующего наименьшему среднему значению импульса и направлению неустойчивости, z = 400000 км.

Рис. 30. Зависимость угла в от угла б, соответствующего наименьшему среднему значению импульса и направлению неустойчивости, z = 600000 км.



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

3.2 Исследование влияния неточности определения параметров КА на геометрию гало-орбиты

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

3.2.1 Исследование случая неточного определения скорости КА

Пусть погрешность определения начальной скорости аппарата составляет 1 см/с. Были рассчитаны 10000 орбит, для которых в начальный момент времени моделировалось отклонение скорости КА на 1 см/с с помощью функции, описанной в п. 3.2.2.

Координаты КА через 1 оборот (175 суток) образуют эллипсоид с центром в точке, лежащей на номинальной траектории. Также были рассмотрены точки через 1/8, 1/4, 3/8 и т.д. оборота (23,125 суток, 46,25 суток и т.д.). На рис. 31 представлен эллипсоид, образованный координатами КА через ј оборота на гало-орбите (46,25 суток).



Рис. 31. Распределение координат КА через 1/4 оборота.

На рис. 32 представлены зависимости максимального отклонения от номинальной траектории для изменений скорости, равных 1 см/с и 5 см/с. Как видно из графиков, максимальное отклонение аппарата от номинальной траектории зависит экспоненциально от времени полета, при этом, коэффициент, стоящий при t, практически не зависит от величины начального отклонения скорости.

Рис. 32. Зависимость отклонения КА от номинальной траектории от погрешности определения скорости КА.

3.2.2 Исследование случая неточного определения координат КА.

Рис. 33 иллюстрирует эволюцию максимального отклонения от номинальной траектории при изменении начального положения аппарата. На рисунке представлены данные для начальных отклонений = 1, 2, 3, 4 и 5 км. Как и в случае с изменением начальной скорости аппарата, максимальное отклонение от номинальной траектории растет экспоненциально с увеличением времени:

, (6)

при этом начальная ошибка влияет только на коэффициент, стоящий при экспоненте .

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

(7)



Рис. 33. Зависимость отклонения от номинальной траектории от погрешности определения координат КА.

Таблица 1.

Зависимость коэффициента k_r от погрешности определения координат КА.

Погрешность, r (км)	kr	Величина достоверности аппроксимации, R2
1	0.9167	0.9994
2	1.8398	0.9994
3	2.7557	0.9994
4	3.6609	0.9994
5	4.5913	0.9994

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

3.3 Результаты расчета направлений устойчивости и неустойчивости

Расчет направления устойчивости производился для 244 плоских орбит Ляпунова, имеющих следующие начальные координаты:

· X = X₀ км, -1200000? X₀?-10000 км;

· Y = 0 км;

· Z = 0 км.

На каждой орбите направление неустойчивости было рассчитано в 360 точках, которые описывались параметром Alpha - углом между радиус-вектором КА и осью OX в плоскости эклиптики.

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

Рис. 34. Карта направлений устойчивости.



По рассчитанным направлениям устойчивости можно получить направления неустойчивости, т.к. они перпендикулярны между собой. Карта направлений неустойчивости изображена на рис. 35.

Рис. 35. Карта направлений неустойчивости.

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



3.4 Интерполяция направления неустойчивости

Зависимость направления неустойчивости от координат X, Y КА образует поверхность, проекции которой представлены на рис. 36-38.

Рис. 36. Точки, для которых рассчитано направление неустойчивости.

Рис. 37. Зависимость направления неустойчивости от координаты X.



Рисунок 38. Зависимость направления неустойчивости от координаты Y.

Для использования этих данных при моделировании движения КА необходимо произвести их интерполяцию. Для этого была разработана функция Matlab. Входными параметрами данной функции являются текущие координаты X, Y КА и угол между радиус-вектором КА и осью X в плоскости эклиптики (Alpha). Выходной параметр - это направление неустойчивости в точке с координатами X, Y.

Точки, для которых рассчитано направление неустойчивости, образуют при проекции на плоскость XY 360 лучей. По параметру Alpha в функции определяются 2 луча, между которыми лежит исследуемая точка. Затем на этих лучах ищутся 3 ближайшие к ней точки, через которые проводится плоскость, после чего в уравнение данной плоскости подставляются координаты исследуемой точки. Полученное значение является интерполированным значением направления неустойчивости.

3.5 Зависимость направления неустойчивости от координаты Z

Орбиты, для которых были рассчитаны направления неустойчивости в предыдущем разделе, лежат в плоскости эклиптики (плоскости XY). Однако также необходимо исследовать, влияет ли координата Z на направление неустойчивости и если влияет, то насколько велико отклонение рассчитанного направления неустойчивости от аппроксимированного с использованием данных из раздела 4.3.

Для исследования зависимости направления неустойчивости от координаты Z по методике, описанной в 4.2, было рассчитано направление устойчивости для точек 57 гало-орбит с различными амплитудами (начальная координата X КА лежит в пределах от -570000 км до -10000 км, начальные координаты КА однозначно определяют гало-орбиту). На рис. 39-41 представлены проекции полученной поверхности на различные плоскости.

Рис. 39. Зависимость направления неустойчивости от координаты X.

Рис. 40. Зависимость направления неустойчивости от координаты Y.

Рис. 41. Зависимость направления неустойчивости от координаты Z.

Из рис. 41 видно, что координата Z значительно влияет на направление неустойчивости. Это означает, что использование данных о направлении неустойчивости из п. 4.3 при моделировании движения КА на гало-орбите менее эффективно. Также был произведен анализ отклонения рассчитанного направления неустойчивости от направления неустойчивости, полученного для плоских орбит в разделе 4.3. Наибольшее отклонение при таком сравнении составляет 12°. Данное отклонение является существенным. Для моделирования движения КА на гало-орбите с исполнением корректирующих импульсов вдоль направления неустойчивости требуется производить интерполяцию данных о направлении неустойчивости, рассчитанных для гало-орбит (т.е. с учетом зависимости от координаты Z).

Для реализации данной интерполяции удобно рассматривать зависимость направления неустойчивости от параметра Alpha, а не координат X и Y. На рис. 42 представлены эти зависимости для гало-орбит с различными начальными координатами.



Рис. 42. Зависимость направления неустойчивости от положения КА на орбите.

Интерполяция полученных данных производилась в Matlab. Входными параметрами функции являются начальная координата X КА и его текущее положение на орбите (Alpha). По начальной координате X КА находятся ближайшие орбиты, для которых рассчитано направление неустойчивости. По текущему значению параметра Alpha находятся ближайшие к нему точки на кривых, соответствующих известным орбитам. Затем по полученным 4 точкам проводится линейная интерполяция. Полученное значение является интерполированным направлением неустойчивости для исследуемой точки и подается на выход функции.

Для проверки данной интерполяции были рассчитаны направления неустойчивости для гало-орбит со следующими начальными координатами:

· X = -277549 км, Y = 0 км, Z = 200000 км;

· X = -373454 км, Y = 0 км, Z = 400000 км;

· X = -566256 км, Y = 0 км, Z = 600000 км;

Затем для 360 точек каждой орбиты были интерполированы направления неустойчивости исходя из данных, полученных выше, после чего была вычислена разница между рассчитанным значением направления неустойчивости и интерполированным Наибольшие отклонения для данных орбит составили 0.008°, 0.0097° и 0.007° соответственно.

3.6 Имитационное моделирование движения КА на гало-орбите с учетом направления неустойчивости

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

Как было сказано выше, в реальности существуют технические ограничения на точность определения вектора состояния КА и выдачу корректирующего импульса. Поэтому при планировании миссии важно рассчитать, как изменится импульс после учета технических ограничений. В табл. 2 приведены результаты расчетов для различных гало-орбит с начальными координатами Z=200000 км, 400000 км и 600000 км и различными направлениями выдачи импульса (в=0°, в=28°, в совпадает с направлением неустойчивости, рассчитанном в п.4.2, в совпадает с направлением неустойчивости, рассчитанном в п.4.4). Считалось, что импульс совершается 1 раз в 40 дней, всего было исполнено 100 коррекций.

Таблица 2.

Суммарный импульс для различных направлений выдачи импульса.

	0°	28°	Направление неустойчивости (п.4.2)	Направление неустойчивости (п.4.4)
Z = 200000 км	0,0374 км/с	0,032 км/с	0,0319 км/с	0,0313 км/с
Z = 400000 км	0,0362 км/с	0,0307 км/с	0,0304 км/с	0,0298 км/с
Z = 600000 км	0,0357 км/с	0,0295 км/с	0,0293 км/с	0,028 км/с

Из табл. 2 видно, что наименее эффективно совершать импульсы в направлении оси Солнце-Земля. Использование интерполяции направления неустойчивости, описанной в п. 4.2, несмотря на то что она не учитывает зависимость направления неустойчивости от координаты Z, позволяет достичь заметной экономии суммарного импульса. Использование направления неустойчивости, рассчитанного для гало-орбит (п. 4.4) позволяет достичь существенной экономии топлива (по сравнению с исполнением корректирующих импульсов в направлении Солнце-Земля): 16%, 18% и 22% для орбит с начальной координатой Z = 200000 км, Z = 400000 км, Z = 600000 км соответственно.

Также вызывает интерес исследование гало-орбит с большими амплитудами. В частности, через несколько лет планируется запуск КА «Спектр-М» («Миллиметрон»). Этот аппарат является космическим комплексом для астрофизических исследований в миллиметровом, субмиллиметровом и инфракрасном диапазонах электромагнитного спектра. «Миллиметрон» планируется разместить на квазипериодической орбите вокруг точки L₂ системы Солнце-Земля.

Для успешной реализации научной миссии проекта требуется, чтобы аппарат находился на орбите, обеспечивающей выход КА из плоскости эклиптики более чем на 1 млн км [21]. Также для корректной работы аппарата требуется, чтобы аппарат пребывал в конусе тени Земли не более одного часа [19]. В связи с этими ограничениями разработчиками миссии было решено вывести КА на гало-орбиту.