Коэффициенты и могут быть как положительными, так и отрицательными. Поэтому термины «возрастающая» и «убывающая» компоненты подразумевают возрастание и убывание по модулю.
Для длительного удержания аппарата в окрестности точки 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, и затем передавать в цикл подбора величины корректирующего импульса эти параметры).
Рис. 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 Алгоритм расчета направления устойчивости
Таким образом, были разработаны следующие алгоритмы:
· Алгоритм подбора начальной скорости КА. Данный алгоритм позволяет подобрать такую начальную скорость аппарата, при которой он будет находиться на гало-орбите наиболее продолжительное время.
· Алгоритм подбора величин корректирующих импульсов. Данный алгоритм основан на алгоритме подбора начальной скорости; основное отличие - учет направления исполнения импульса.