Об оптимальном одноимпульсном маневре уклонения

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

«Московский физико-технический институт (пациопальпый исследовательский университет)», 2026

Среди наиболее примечательных работ по этой теме можно выделить [1]. В ней Патера и Питерсон предложили оптимизационный алгоритм, который работает с мгновенными маневрами уклонения и состоит из двух отдельных этапов. Сначала по градиенту поля вероятности столкновения определяется направление импульса маневра. Затем по уже известному направлению импульса его величина определяется итерационно с помощью метода Ньютона для линеаризованной функции вероятности столкновения с учетом ограничения на ее значения. Разделение таким образом этих этапов позволяет свести оптимизацию к одномерной задаче. Этот алгоритм получил развитие в работе [2]. В ней вероятность столкновения представляется как функция управляющего воздействия, зависящая от произвольного количества параметров — в результате ограничение получается полиномиальным. Кроме того, используется квадратичная целевая функция. Такой подход приводит к решению полиномиальной задачи оптимизации. Этот метод применим как к мгновенным импульсам, так и к малой тяге. Он также учитывает изменение времени наибольшего сближения (Time of Closest Approach — ТСА) КО ввиду совершаемого маневра. В работе [3] этот метод был дополнительно улучшен для учета нескольких ограничений, а также работы в случае множественных опасных сближений.

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

Существует множество методов расчета вероятности столкновения в момент опасного сближения. Большинство методов основываются на приближенном вычислении двумерного интеграла плотности вероятности столкновения по области в плоскости, перпендикулярной вектору относительной скорости объектов. В таких методах часто пренебрегают геометрическими особенностями КО — вместо этого рассматривают сферы с характерными размерами этих объектов. Метод Хуторовского [4] предлагает аналитическую формулу расчета вероятности столкновения. В методе Чана [5] формула вероятности столкновения представляется в виде суммы сходящегося бесконечного ряда за счет аппроксимации двумерного нормального распределения плотности вероятности одномерным распределением Райса. Методы Чена - Ляна - Бая - Ли для круговых и эллиптических орбит [6] основаны на методе Чана, но ограничиваются только двумя первыми членами ряда. Метод Патеры [7] с помощью специальных преобразований переходит от расчета двумерного интеграла вероятности столкновения к интегралу по контуру. В методе Фостера [8] за счет перехода к полярной системе координат (СК) расчет вероятности сводится к вычислению повторного интеграла. Наконец, метод Альфано [9] использует ряд из функций ошибок для аппроксимации интеграла плотности вероятности столкновения КО.

Расчет вероятности столкновения двух КО предполагает определение их ТСА. Делать это можно методом полного перебора расстояний между объектами в зависимости от времени — такой подход является вычислительно затратным из-за необходимости интегрировать орбиты КО, а его точность зависит от шага интегрирования. Другим вариантом является метод ANCAS, описанный в [10, 11]. Этот алгоритм аппроксимирует относительную скорость кубическим полиномом для определения ТСА, а затем аппроксимирует расстояние полиномом пятой степени для оценки дистанции в этот момент. ANCAS привносит ошибки в вычисления ТСА, связанные с приближенным расчетом расстояний. В настоящей работе используется модифицированная версия алгоритма, SBO-ANCAS [12]. В ней применяется метод суррогатной оптимизации (surrogate-based optimization — SBO) [13]. Подобный подход является компромиссом между скоростью вычислений и их точностью.

2.    Постановка задачи

Рассматривается система из двух КО. Считается, что первым объектом можно управлять и именно скороств первого объекта будет изменяться под действием импульсного маневра. Вводится функция вероятности столкновения КО в момент их наибольшего сближения p(Av, t), где Av — изменение скорости под воздействием импульса в момент времени t.

Задачей поиска оптимального одноимпульсного маневра уклонения является нахождение такого импульса скорости Avopt и момента его приложения t* на заданном временном отрезке [tstart, tend], который обеспечит во время наибольшего сближения tcA итоговую вероятность столкновения p(Avopt,t*) <  р* при минимальных затратах характеристической скорости; здесь tstart < tend < tcA, р* — заданное целевое значение вероятности столкновения. Таким образом, минимизация происходит на множестве Q = R³ х [tstart,tend].

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

min n{|Av| : -(Av,t) < -*}

(^ v,t ) е п

3.    Поле вероятности столкновения

В предыдущем разделе для удобства постановки задачи вероятность столкновения была введена как функция p(A v ,t). Теперь же интерес представляет и мгновенное значение, зависящее от положения и скорости КО в момент опасного сближения. Вероятность столкновения можно рассматривать как функцию Р = Р ( r , v ), г де r и v — радиус-векторы относительных положения и скорости КО, спрогнозированные к моменту ТС А с учетом совершаемого в момент t < t cA маневра с импульсом A v.

Из всех перечисленных ранее способов расчета вероятности столкновения [4-9] больше всего для решения поставленной задачи подходит метод Хуторовского [4]. Во-первых, метод прост в реализации. Во-вторых, метод представляет собой аналитическую формулу, которую можно дифференцировать. Далее указана расчетная формула вероятности столк новения.

^р ( r . ^v) = к • exp (^-2-(^k„^- A )) ,

здесь k.k rr , k_vv и k_rv — скалярные величины:

к = ^-lv (kvv • det(Ki) • det(K2) • det(K-1 + K2 1)) 0'5 k„ = rT • (Ki + K2)-1 • r, kvv = vT • (Ki + K2)-1 • v, krv = rT • (Ki + K2)-1 • v, где Ki и K2 — матрицы ковариации ошибок определения положений первого и второго КО в момент ТС A, S — площадь поперечного сечения тела суммарных размеров КО.

Как уже было сказано, формула (1) позволяет вычислять мгновенное значение вероятности столкновения, однако для построения алгоритма поиска оптимального маневра удобно ее представление в виде квадратичной формы (КФ) относительно вектора скорости маневра Av:

p(A v ,t)= p( 0 ,t)+    Р I • A v +— AvT^T • P 1     • A v^j ,       ( 2)

,       , 7   dAv kv=o       2        dAvdAvT kv=o      , где p(0,t) — вероятность столкновения без маневра, рассчитанная по формуле (1). Отсюда и далее считается, что го и vo — положение и скорость в момент времени t < tcA. Дальше в настоящем разделе второе и третье слагаемые выражения (2) будут рассмотрены по отдельности.

• 3.1.    Градиент поля вероятности столкновения

Градиент во втором слагаемом из (2) можно представить следующим образом:

др I _ дР     d( r , v )   d( r o , v_o )

дA v l A v = o д( r , v ) д( r o , v o )     д v₀

Первый множитель из (3) является вектором-строкой и получается прямым дифференцированием формулы (1):

дР д(г, v )

(дРЩ

^_ Р^(0, 1) • ⁽ K 1 + K 2 ) ¹ • ⁽cv ^{- r)}, д r

_ р(0, *) • (iv2 — (K1+K2)-1 • (£v- cr+c2v)), здесь с _ k„/kvv.

Второй множитель из (3) можно переписать так:

д( r , v )     д( r , v )   дЭ     дЭ ₀

д(ro, vo)      дЭ    д'ф  д(ro, vo)’ где Эo и Э — элементы орбиты в момент времени to и tcA соответственно. Вопрос о том, как считать производные из последней формулы, разобран в [14-16].

Наконец, третий множитель из формулы (3) равен блочной матрице 6 хб:

д( r o , v o ) д v o

10з М

где 0з — путевая матрица ЗхЗ. I3 — единична я матрица. ЗхЗ.

• 3.2.    Матрица Гессе поля вероятности столкновения

Матрица Гессе в третьем слагаемом формулы (2) представляется следующим образом:

д 2Р     I дAvдAvT Lv=o

д( r , v ) ^Т ^ H ^ д( r , v ) д v o          д v o

дР д / д ( r , v ) \

⁺ д ( r , v ) д v o V д v o     ’

д( r , v ) _ д( r , v )    д( r o , v o )

^ГД<   д v o      д ( r o , v o )     д v o

H — матрица вторых производных вероятности столкновения из (1) по r 11 v. Формулы вычисления элементов этой матрицы (П1) - (ПЗ) приведены в приложении.

Вторым слагаемым в (4) можно пренебречь — это доказывается результатами тестирования в разделе 5. Тогда д2р    l _ д(r, v)Т ^ H ^ д(r, v)

дAvдAv^Т I a v = o     д v o         д v o

• 4.    Численная минимизация

4.1.    Оптимизация момента приложения импульса

Алгоритм нахождения оптимального маневра уклонения можно разделить на две взаимосвязанные части, одна из которых отвечает за оптимизацию момента приложения импульса, а другая — за оптимизацию скорости этого импульса. Первая из них, алгоритм 1, проводит минимизацию затрачиваемой характеристической скорости относительно момента приложения импульса на интервале [tstart, tend]. Минимизация выполняется методом полного перебора на сетке значений. Вторая часть, алгоритм 2, проводит минимизацию характеристической скорости маневра при фиксированном времени его приложения. Минимизация происходит с условием понижения значения вероятности столкновения до заданного порога р*.

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

Рассматривается первая часть алгоритма. На вход подаются следующие параметры: tcA — ТСА в отсутствие маневра, рассчитанное предварительно (в настоящей работе с помощью алгоритма из работы [12]), tstart и tend — начало и конец расчетного временного отрезка, At — шаг разбиения расчетного отрезка для минимизации методом полного перебора, ri,vi и Г2, V2 — положение и скорость в момент ТСА первого и второго КО соответствешю. Ki 11 K2 — матрицы ковариации ошнб(ж определения их положения. S — площадь поперечного сечения тела суммарных размеров КО. Подразумевается, что на вход еще подается ряд дополнительных параметров, нужных для работы алгоритма 2, однако тут они убраны из рассмотрения, чтобы избежать лишней нагруженности, — вместо них на схеме алгоритма изображено «...», а сами они будут описаны в следующем разделе.

Выходными данными алгоритма являются: P_res — значение вероятности столкновения после маневра, Avopt — вектор скорости оптимального маневра, t* — момент приложения импульса маневра, tcA — новое ТСА с учетом маневра.

Схема первой части алгоритма представлена далее (рис. 1).

Алгоритм 1: Оптимизация момента приложения импульса Вход: t_CA, tstart, tend, At, щ, V_b Г₂, V₂, Ki, K₂, S, ...

Выход: P_ies, Av_opt, t*, t^_A

1 p(O,tcA) «- F(n,vi,r2,v2,Ki,K2,S);

2 P^s <- р(ОАса), Av_opt<- {00,00,00};

3    t < tstartj 4    до тех пор пока t < tend выполнять 5    P, Av, ^дЧ-МИНХССХО^саМсаЛгьУьГз^К!,^ 6    если |Av| < Avopt| то PTes e- P, Avopt <— Av, t* <— t; 7    t <- t + At; 8    конец 9    вернуть Pres, Avopt, t*, t^A;

Рис. 1. Алгоритм оптимизации момента приложения импульса

Здесь P — функция расчета мгновенной вероятности столкновения по формуле (1), МИНХС — функция минимизации характеристической скорости маневра, совершаемого в заданный момент времени t; она же является алгоритмом 2.

4.2.    Минимизация характеристической скорости

Рассматривается вторая часть алгоритма. Как уже было сказано выше, минимизация происходит для конкретного момента времени t < tcA. Сначала прогнозируются положение и скорость первого КО в момент t. Изначально импульс маневра уклонения инициализируется нулевым вектором. После этого запускается сам цикл минимизации. В рамках этого цикла сначала по положениям и скоростям КО в моменты времени t и tcA с учетом приложения импульса маневра вычисляются градиент и матрица Гессе поля вероятности из КФ, описанной в разделе 3. Затем эти вектор и матрица подаются на вход функции минимизации квадратичной формы (МИНКФ) наряду с еще одним параметром s, ответственным за размер области поиска решения — о нем, как и о самой функции, будет подробнее сказано далее. В результате работы функции МИНКФ обновляется вектор импульса маневра, после чего он используется для определения нового ТСА в окрестности старого. По вычисленным вектору импульса маневра, а также новому ТСА с помощью прогнозирования пересчитывается мгновенная вероятность столкновения КО. Дополнительно, если найденный на актуальной итерации импульс находится достаточно близко к границе области поиска решения функцией МИНКФ, что регулируется параметром е С 1, область поиска для следующей итерации расширяется в (1 + а) раз, где а > 0. После этого цикл повторяется, пока значения вероятности столкновения не уменьшатся до заданной пороговой величины.

Параметры, подающиеся на вход алгоритма 2, помимо тех, что уже были обозначены в описании алгоритма 1: t — момент времени, для которого проводится минимизация скорости, д ( ОДса ) — вероятность столкновения в отсутствие маневра, вычисленная в алгоритме 1, р* — целевое значение вероятности столкновения, Avmin и Avmax — нижнее и верхнее ограничения по характеристической скорости.

Выходными данными алгоритма являются: Р — значение вероятности столкновения после совершения рассчитанного маневра, Av — вектор минимизированной скорости импульса, tCA — пересчитанное время ТСА с учетом совершаемого маневра.

Алгоритм 2: Минимизация характеристической скорости (МИНХС)

Вход: t, t_CA, р(ОДсл), P*, ri, v_b r₂, v₂, К_ь К₂, S, Av_min, Av_max

Выход: Р, Av, ^_А

1 Р ^ /ДОДсаМса ^СА, Av t- {0, 0, 0}, s ^ Av_min;

2 го, vo ■«— ПРОГНОЗ(г1^1ДсаД);

з до тех пор пока Р > р* и s < Av_max выполнять

5 в

Ь. А <— РАЗЛОЖЕНИЕ(гц, vq + Av. гр vp r₂, v_2t Kp K₂,5);

Av <^МИНКФД,Ь. А);

гр vi ^— ПРОГНОЗ До, v_n + Av, t, tc_A);

^CA ^— СБЛИЖЕНИЕД1, vi, r2, v2, tcA — ^C tcA + •JOi гф V? ^ ПРОГНОЗ (r0, Vo + Av, t,t^A);

гф v^ ^ ИНТЕРПОЛЯЦИЯ^);

P ^ПДф^,гф^КрКрД);

если |Av | ^ (1 — e)s to s t— (1 + a)s;

12 конец

13 вернуть P, Av, tg_A;

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

На рис. 2 представлена схема второй части алгоритма. Функция ПРОГНОЗ(r, v,tpt2) служит для прогнозирования положения и скорости первого КО спустя время t2 — ti в зависимости от его нынешних положения и скорости r(ti) и v(ti). Функция ИНТЕРПОЛЯЦИЯ служит для получения положения и скорости второго КО для заданного момента времени. РАЗЛОЖЕНИЕ — функция расчета градиента и матрицы Гессе поля вероятности столкновения по формулам (3) и (5). Функция МИНКФ для вектора b и матрицы A решает задачу минимизации Trust Region Subproblem (TRS) с помощью метода Double-start FOCM, описанного в работе [17]:

min q( x ), где

| x |d

q( x ) = x^T Ax — 2 b ^T x.                               (6)

Параметр s нужен для приведения КФ (2) к виду (6) и поиска таким образом решения

TRS на единичном шаре. Соответствующие формулы приведения:

*                *            1

где

A = у A , b = -2 • b , x = ^_ • △ v ,

A = d( r , v ) ^T ^ h ^ d( r , v )     b = dP ^ d( r , v )

d v ₀           d v ₀ ,        d ( r , v )    d v ₀

Функция СБЛИЖЕНИЕ ( r i, v i, r 2, v 2 ,t i ,t 2 ) рассчитывает по заданным положениям и скоростям обоих КО ТС А в пределах временного отрезка [ti,t2] — конкретно для алгоритма 2 ТС А пересчитывается в окрестности tcA, поэтому этот отрезок задавался как [tcA - At, tcA + At ], В рамках реализации алгоритма для настоящей работы полагалосв, что е = 0.01, _а = 0.05.

5.    Исследование алгоритма

В целях проверки работоспособности алгоритма были выбраны два показательных случая: искусственный случай столкновения двух КО на заданных орбитах и реальное задокументированное столкновение спутников «Космос-2251» и «Iridium 33» [18]. Каждый из случаев будет рассмотрен подробно в настоящем разделе. Дополнительно будут представлены результаты работы алгоритма для большого количества случайно сгенерированных орбит.

5.1.    Искусственный случай столкновения

В качестве первого случая рассматривается столкновение КО, один из которых находится на круговой орбите, а другой — на околокруговой орбите. Орбиты расположены так, чтобы иметь одну точку пересечения, в которой объекты встретятся в момент столкновения. В таблице 1 представлены параметры орбит КО. Отсюда и далее конкретно для этого примера все величины, связанные с расстоянием, временем и скоростью, представлены обезразмеренными на a i — большую полуось орбиты первого КО, T i — период его обращения на этой орбите, а также vi — соответствующую скорость.

Таблица 1

Параметры орбит

Орбита	a	е	i, рад	ш, рад	Д рад	Т
1	1.0	0.0	0.0	5.655351	0.0	1.0
2	1.064849	0.060899	1.570796	3.141593	3.141593	1.098834

Здесь а, е, i, ш и Q- стандартные элементы кеплеровой орбиты: большая полуось, эксцентриситет, наклонение, аргумент перицентра и долгота восходящего узла соответственно; T — период обращения КО на орбите.

Таблица 2

Дополнительные параметры

КО	v, рад	d, 10-7	aN, 10-6	ат, 10-6	aR, 10-6
1	3.879763	1.47	4.65	4.65	4.65
2	5.648888	1.47	4.65	4.65	4.65

В таблице 2 представлен ряд дополнительных параметров: v — истинные аномалии объектов в начальный момент времени t = 0.0, которые были заданы так, что КО оказались разведены на сутки от момента достижения точки пересечения их орбит, d — габариты объектов, а также aN, ат а ад — среднеквадратичные отклонения ошибок определения их положений в момент предполагаемого столкновения.

ТС А без маневра составило tcA = 15.482440 х Т1, а расстояние во время наибольшего сближения, ожидаемо, было равно rmin = 0.0. Значение вероятности столкновения в отсутствии маневра составило 2.5 • 10^-4, что больше заданного целевого значения вероятности столкновения р* = 1.0 • 10^-6.

После выполнения маневра уклонения, определенного предлагаемым алгоритмом, вероятность столкновения уменьшилась до 9.3 • 10-7 за счет совершения маневра с характеристической скоростью |Avopt | = 1.1 • 10-7 х vi в момент времени I * = 0.012544 х Т 1. Минимальное расстояние между КО при этом увеличилось до rm in = 2.1 • 10^-5 х а1, а новое ТС А составило t CCA = 15.482437 х Т1.

Была проведена проверка оптимальности полученного решения методом полного перебора возможных маневров в момент времени t* — отдельно по направлению и величине импульса. Сначала рассматривалось множество всевозможно направленных векторов скорости импульса фиксированной величины |Avopt | — для них рассчитывались соответствующие вероятности столкновения в момент ТС А. Результирующий набор значений вероятности столкновения был отображен в виде тепловой карты зависимости вероятности столкновения в момент ТСА от направления импульса маневра фиксированной величины в сферической системе координат, построенной на орбитальной СК первого КО.

Рис. 3. Тепловая карта вероятности столкновения для искусственного случая

На рисунке 3 представлена тепловая карта р(р, в ) в логарифмическом масштабе. Оси соответствуют стандартным сферическим углам для орбитальной СК — направление (0, 0) соответствует направлению вектора скорости КО. Маркером «pmin» обозначено направление импульса, рассчитанного с помощью алгоритма, которое соответствует минимальному значению на карте. Таким образом, оптимальное направление импульса выбрано корректно.

Затем была проведена проверка оптимальности характеристической скорости совершаемого маневра. Для каждого значения характеристической скорости в диапазоне от 0 до 5 х |Avopt | была найдена минимальная по направлению приложения импульса в момент времени t* реализуемая вероятность столкновения — в результате была получена кривая зависимости pmin(|Av|). Соответствующий график представлен далее (рис. 4), на нем также показаны графики подобной зависимости для разных отрезков [tcA — At, tcA], на которых минимизация проводилась для своих значений t*. Нумерация графиков соотносится с различными значениями At в легенде (в периодах Т1). Пунктирной линией очерчено целевое значение вероятности столкновения р* = 1.0 • 10 6, а точки на каждой из кривых — решения, полученные при помощи алгоритма.

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

Ниже, в табл. 3, прилагающейся к рис. 4, представлены подробные данные для каждого из решений: номер кривой, длина расчетного отрезка At, результирующая вероятность столкновения р, время совершения оптимального маневра t , величина его импульса |Av_opt |, а также sv н E_a — относительное по величине и абсолютное по углу отклонения скорости импульса, полученного с помощью алгоритма, от вектора скорости, полученного методом полного перебора для итогового значения вероятности столкновения р.

Таблица 3

Результаты проверки решений для искусственного случая

№	At	р, 10-7	t*	|AvOpt |, 10-8	E v ,10-6	Ec, рад
1	0.242	6.0	0.0	937.9	15.5	0.007
2	0.484	6.2	0.0	290.8	10.9	0.009
3	0.968	7.6	0.045	161.9	8.4	0.006
4	1.935	7.8	0.005	81.8	151.2	0.007
5	3.871	7.2	0.0	41.3	8.2	0.008
б	7.741	7.5	0.009	20.9	6.8	0.010
7	15.482	9.3	0.013	10.5	16.1	0.001
8	30.965	7.9	0.005	5.3	160.3	0.007

По результатам проверки решений (табл. 3) видно, что максимальное отклонение по величине импульса составило порядка 0.01%, а максимальное отклонение по направлению его вектора не превысило 0.01 рад.

5.2.    Реальный случай столкновения

Вторым случаем, на котором испытывался алгоритм, стало реальное столкновение КО «Космос-2251» и «Iridium 33» на низкой околоземной орбите (НОО) в 2009 г. [18]. Необходимые параметры положений и скоростей объектов в момент ТС А, а также соответствующие матрицы ковариации были взяты из [19]. Полагалось, что момент времени t = 0.0 с на сутки опережал столкновение. ТСА в отсутствии маневра, таким образом, составило t cA = 86400.0 с — соответствующее значение вероятности столкновения оказалось равным 3.8 • 10^-4. Дальнейший расчет проводился для целевой вероятности столкновения 1.0 • 10^-6.

Решением алгоритма стал импульс величиной 4.7 мм/с с моментом приложения t* = 10.0 с, который понизил вероятность столкновения до 9.7 • 10^-7.

Для «Космос-2251» и «Iridium 33» была проведена такая же проверка оптимальности полученного решения, как и для первого, искусственного, случая. Соответствующие графики, а также таблица представлены ниже.

Рис. 5. Тепловая карта вероятности столкновения для «Космос-2251» и «Iridium 33»

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

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

Данные из табл. 4, прилагающейся к рис. 6, указывают на то, что решения получены с точностью не хуже, чем в табл. 3. Поэтому можно утверждать, что при помощи алгоритма удалось отыскать маневр уклонения, предотвращающий реальное столкновение.

Рис. 6. График минимальной вероятности столкновения для «Космос-2251» и «Iridium 33»

Таблица 4

Результаты проверки решений для «Космос-2251» и «Iridium 33»

№	At, сут	р, 10-7	t*, с	|AvOpt |, мм/с	S v ,10-5	Sc, рад
1	0.125	9.4	1180.0	31.3	17.5	0.002
2	0.250	6.6	0.0	18.3	3.4	0.005
3	0.500	9.9	630.0	9.3	10.3	0.001
4	1.000	9.7	10.0	4.7	1.8	0.001
5	2.000	7.0	70.0	2.4	2.8	0.004

5.3.    Анализ для большого количества орбит

Дополнительно работоспособность алгоритма была проверена на наборе из 100 000 пар случайно сгенерированных конфигураций КО, для которых изначальная вероятность столкновения была выше целевой р* = 1.0 • 10^-6. Сделано это было следующим образом.

Сначала для каждой пары с помощью алгоритма находился оптимальный маневр A v_opt, время его совершения t* и результирующая вероятность столкновения р < р*. Затем с помощью метода полного перебора в момент времени t* рассчитывались всевозможные маневры с характеристическими скоростями |Av| ^ |Avopt | — среди них отбирались те, которые дают вероятность столкновения, отличающуюся от р не более чем на 5%, и не превышающую р*. Далее, если соответствующая характеристическая скорость оказывалась меньше, чем |Avopt |, для нее считалось относительное отклонение.

Таким образом, для каждого значения характеристической скорости рассчитывалось максимальное по всем направлениям импульса отклонение от оптимального решения — результатом подобного анализа стала гистограмма, представленная ниже (рис. 7). Горизонтальная ось соответствует максимальному отклонению от найденного решения, а вертикальная — процентному соотношению случаев с таким значением отклонения.

Как можно видеть на рис. 7, в основном погрешность определения характеристической скорости составляет не более 1 % — наиболее частое отклонение, почти 25 % от всех случаев, и вовсе находится вблизи нуля. Глобальное максимальное значение ошибки по всем 100 000 парам при этом не превышает 5 %.

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

6.    Заключение

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

Алгоритм был протестирован. Работоспособность алгоритма была продемонстрирована на примере как искусственно заданных орбит, так и реального случая столкновения спутников на НО О. Дополнительно анализ большого количества орбит методом полного перебора по всевозможным импульсам маневров показал высокую эффективность алгоритма в нахождении оптимального маневра уклонения.

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

7.    Научный вклад авторов

• •    В. В. Чугунов. Разработка алгоритма поиска оптимального одноимпульсного маневра уклонения для КО на околоземных орбитах. Численная реализация алгоритма в виде программного кода. Подготовка и написание текста статьи.

• •    В. Н. Чиняев. Обзор существующей литературы по теме маневрирования КО на околоземных орбитах. Написание введения статьи. Проведение тестирования алгоритма.

• •    А. А. Кузнецов. Исследование методов решения задачи TRS. Реализация метода Double-start FOCM.

8.    Приложение

Формулы вычисления элементов матрицы H из (5):

%? = же^ s 0 s^ • q + v ' ^Q - 1 з d r                                  k vv

—— = p( 0 ,t) - Q -( B + C ), drd v

d²P d v ²

dP _    .   .

=    0 D +p( 0 ,t)

dv

1        2• v 0 v^T

| v |^{2 I3-}      | v |⁴

-

Д • E - c ² Q k vv

(П1)

(П2)

(ПЗ)

где, B

= s 0

IV

v^T — —• v^T • Q — cs^T • Q^j , kvv                   J

C = v 0 f r — 2^ v }  • -Q + с 1 з ,   D = 12 v^T — f7“ ‘ ^v + c s } - Q ,

\     kvv J    kvv                 | v |²            kvv         J

^E = ^I 3 +

22cr 0 v^T + 2cv 0 ( r — 2c v ) ^T —

2 • v 0 vT vv

-

r 0 r^T^

• Q ,

Q = ( K i + K 2 ) 1 , s = c v — r , c = k_n/k vv .