Вторая методика, Floquet mode approach, впервые предложенная в [16], описана в [13] в случае ограниченной круговой задачи трех тел. Для уравнений движения КА выводится матрица монодромии, собственным числам и векторам которой соответствуют различные характеристики. Наиболее важной характеристикой, с точки зрения управления КА, является направление неустойчивости. Оно совпадает с первым собственным вектором матрицы монодромии. Импульсы, совершенные в направлении неустойчивости, являются наиболее эффективными. Соответственно, расчет величин импульса коррекции в дальнейшем ведется в предположении, что импульсы совершаются вдоль направления неустойчивости. Как и в предыдущем случае, для данной методики вводятся ограничения на минимальный интервал между импульсами и минимальное отклонение от номинальной орбиты.
Приведенная в [14] методика была разработана для миссии ARTEMIS. В данной статье также рассчитываются направления устойчивости и неустойчивости с использованием матрицы монодромии. Для удержания КА на орбите используется т.н. optimal continuation strategy (стратегия оптимального продолжения). Она разработана для удержания КА на гало-орбите вокруг точки L2 системы Земля-Луна на протяжении 1-2 оборотов. Это достигается благодаря подбору такой величины импульса, при которой выполняются заданные оператором условия (например, определенная координата и скорость по одной из осей) при достижении аппаратом плоскости орбиты Луны.
При разработке стратегии удержания КА на гало-орбите важно знать направление неустойчивости, т.к. импульсы в этом направлении являются наиболее эффективными. В литературе предлагается методика расчета направления неустойчивости, связанная с матрицей монодромии уравнений движения и модами Флоке. Данный метод был описан в [17]-[19].
2. Стратегия удержания КА на гало-орбите вокруг точки L2 системы Солнце-Земля
2.1 Математическая модель
Для описания движения КА по ограниченной орбите введем вращающуюся систему координат, связанную с точкой L2. Центр системы координат расположен в точке L2, ось X совпадает с осью Солнце-Земля и направлена от Солнца к Земле, ось Z направлена в северный полюс эклиптики, ось Y дополняет систему координат до правой тройки. Иллюстрация системы координат приведена на рис. 2. Далее в работе всегда подразумевается данная система координат, если не оговорено иное.
Рис. 2. Вращающаяся система координат, связанная с точкой L2.
Движение КА в системе n тел в инерциальной системе координат описывается следующими уравнениями:
(1)
где n - количество притягивающих центров, G - гравитационная постоянная, R - радиус-вектор КА, mi - масса i-го тела, Ri - радиус-вектор i-го тела.
Во вращающейся системе координат, при переходе к ограниченной задаче трех тел () уравнения (1) могут быть приведены к виду [1]:
(2)
где c2 - параметр, зависящий от масс тел, ax, ay, az - возмущающие ускорения, которые являются функциями координат КА и эксцентриситета орбиты КА.
Уравнения (2) допускают линеаризацию в окрестности точки L2. Тогда они принимают следующий вид [2]:
(3)
Система (3) имеет следующее решение:
(4)
где , , , , , - фаза колебаний в плоскости XY, - фаза колебаний по оси Z. Коэффициенты , , , и фазы , зависят от начальных условий.
Решения и представляют собой линейную комбинацию трех компонент: ограниченной (и ), экспоненциально возрастающей и убывающей к нулю . Будем предполагать, что и решение системы (1), записанной во вращающейся системе координат, связанной с точкой L2, также представимо в виде линейной комбинации ограниченной, возрастающей и убывающей компонент:
где , - возрастающие компоненты, , - убывающие компоненты, , , - ограниченные компоненты.
Коэффициенты и могут быть как положительными, так и отрицательными. Поэтому термины «возрастающая» и «убывающая» компоненты подразумевают возрастание и убывание по модулю.
Для длительного удержания аппарата в окрестности точки L2 требуется, чтобы коэффициент возрастающей компоненты равнялся нулю.
Ограниченные орбиты, периоды обращения КА по которым по осям 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) и движется перпендикулярно ей (Vx=Vz=0 км/с). Положим также значение скорости вдоль оси Y равным некоторому значению Vy0. Построим виртуальные плоскости X=Xmin и X=Xmax такие, что гало-орбита лежит между ними, и будем интегрировать уравнения движения КА до того, как КА достигнет любой из указанных плоскостей. Если аппарат достиг левой границы (X=Xmin), то значение коэффициента при возрастающей компоненте отрицательно, и наоборот. Т.о. определяется функция - конечная координата X КА от начальной скорости Vy0. Эта функция терпит разрыв, т.е. КА не достигает ни левой, ни правой границы, при некотором значении Vy0. Задача заключается в отыскании этого значения. Данная задача решается методом деления отрезка пополам. Работа алгоритма заканчивается при достижении минимальной вычислительной погрешности, обеспечиваемой средой вычислений . Иллюстрация работы алгоритма представлена на рис. 7.
Рис. 7. Подбор начальной скорости КА.
Подбор величин корректирующих импульсов подбирается по схожему алгоритму. Отличием является учет направления исполнения импульса: если в задаче нахождения начальных условий меняется только проекция вектора скорости на ось Y, то в задаче нахождения величины корректирующего импульса требуется изменение проекций вектора скоростей на оси X и Y. Подробнее методика расчета описана в следующем разделе.
Разработанный алгоритм позволяет рассчитывать начальную скорость КА и величины корректирующих импульсов для любых ограниченных орбит в окрестности точки L2 системы Солнце-Земля. При этом он позволяет рассчитывать импульсы коррекции в любой точке орбиты в любом направлении.
2.3 Реализация стратегии удержания КА
2.3.1 Алгоритм подбора начальной скорости и величины корректирующего импульса
Описанные алгоритмы были реализованы в программе GMAT (General Mission Analysis Tool). Данная программа позволяет интегрировать уравнения движения КА в реалистичной модели сил различными численными методами, в т.ч. методами Рунге-Кутта различных порядков, Принса-Дорманда и др. В данной работе численное интегрирование проводилось методом Рунге-Кутта 8-9 порядка. В GMAT были созданы сценарии, в которых производился подбор начальной скорости КА и величин корректирующих импульсов и интегрировались уравнения движения КА. Метод деления отрезка пополам позволяет достичь любой наперед заданной точности вплоть до машинного ограничения 10-16.
Например, подберем начальную скорость для КА, имеющего следующие начальные параметры:
· X = -277548 км;
· Y = 0 км;
· Z = 200000 км;
· Vx = Vz = 0 км/с.
Алгоритмом была найдена скорость Vy, равная -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, и затем передавать в цикл подбора величины корректирующего импульса эти параметры).