УДК 681.51
Оценка угловой скорости линии визирования в процессе сближения космических аппаратов по результатам измерения дальности и скорости продольного движения
© Н.Е. Зубов1,2, Е.А. Микрин1,2, А.С. Олейник1
1 ОАО «Ракетно-космическая корпорация "Энергия" имени С.П. Королёва», г. Королев Московской области, 141070, Россия 2 МГТУ им. Н.Э. Баумана, Москва, 105005, Россия
Рассматривается задача оценки угловой скорости линии визирования в процессе сближения космических аппаратов по результатам измерения дальности и скорости продольного движения. Получено ее аналитическое решение. Приведены численные решения данной задачи.
Ключевые слова: сближение КА, дальность, угловая скорость линии визирования, скорость продольного движения, метод точного размещения полюсов.
Введение. В пилотируемой космонавтике сближение космических аппаратов является штатной операцией, при этом она должна обладать повышенной надежностью реализации. В качестве измерительных средств параметров взаимного сближения в практике отечественных полетов используют радиотехнические системы. Несмотря на высокую надежность этих систем, возможность их отказов существует. В данной статье рассматривается один из возможных таких отказов, связанный с отсутствием измерений угловой скорости линии визирования. В этом случае для продолжения процесса сближения необходимо иметь бортовой алгоритм оценки угловой скорости линии визирования. Решению этой проблемы и посвящена данная работа.
1. Уравнения относительного движения в орбитально-луче-вой системе координат. Рассмотрим орбитально-лучевую систему координат (рис. 1). Ось x направлена вдоль вектора дальности от цели до перехватчика (вдоль линии визирования), ось y лежит в плоскости орбиты цели, а ось z дополняет систему координат до правой тройки.
Векторы угловой скорости вращения орбитально-лучевой системы координат относительно инерциального пространства и относительной дальности имеют следующие компоненты:
Q =
ß ß ß
t
(1.1)
D = [D, 0, 0]T, (1.2)
Плоско сл орбиты
Рис. 1. Орбитально-лучевая система координат
Выполняя последовательно дифференцирование, получим:
Б = [Б, 02Б, -ПуБ]т,
(1.3)
Б = [Ъ - (П2у +П2 ) 5 (0 + 0Х0у )Ъ + 202Ъ, (-0 + 0)Ъ - 2ПуЪ]т.
Векторное уравнение взаимного сближения движения записывается так:
Здесь g — вектор относительного гравитационного ускорения; a — вектор управляющего ускорения.
Подставляя (1.3) в уравнение (1.4) и принимая во внимание тот факт, что на небольших расстояниях между перехватчиком и целью можно пренебречь относительным гравитационным ускорением, получим систему уравнений относительного движения в орбитально-лучевой системе координат:
2. Линеаризация системы уравнений относительного движения. Уравнение невязок. Перепишем систему (1.5) следующим образом (без учета управляющих воздействий):
Б = g + а.
(1.4)
Б - (П2у + 02 )Б = ах, (О г + 0 х О у ) Б + 20 г Б = а у,
(-0у + 0х02)Б - 20уБ = а
(1.5)
(2.1)
Х3 —
Х4 —
2 Хо Хо _
-— + 0 X Х4,
2Х4Х2 -0
(2.1)
где х = [ х1 х2 х3 х4 ]т = [ В В 0.у 02 ]т.
Считаем, что на борту активного КА имеется бортовая вычислительная машина с тактом вычисления At. Для линеаризации системы (2.1) в окрестности значения вектора состояния на момент начала такта работы бортовой машины воспользуемся разложением в ряд Тейлора [1]:
— Х2,
Х2 = ([Хз2 + Х42]Х° + (Хз2 + Х42° (Х! - х°) + (2ХХ1 ^0 (Х - Х°) +
+ (2Х4Х1 ^° (х4 - Х4° X
Х3 -
2 Х3 Х2
V Х1
Л
+ ®х Х4
+
^2х3Х2 >
° V Х1 У
(Х1 - Х° ) -
2 Х3
V Х1 У
(х2 Х2 )
^ 2 Х2 ^
V Х1 У
(х3 - Х3° ) + ®х (х4 -Х40).
Х4 —
2 Х4 Х2
V Х1
Л
-®х Х3
+
'2 Х4 Х2 ^
° V Х1 У
(х1 - х° ) -
' 2 х >
V Х1 У
(х2 Х2)
-®х (Х3 - Х3° ) -
^ 2 Х2 ^ V Х1 У
(х4 Х4 ^
Запишем линеаризованную систему уравнений относительного движения в виде
где
х = Ах + К
о»
(2.2)
А =
0
2 2 Х3 + Х4
2 Х3 Х2
2 Х4 Х2
1 0 0
0 2 Х3 Х1 2 Х4 Х1
2 Х3 2 Х2 О X
х1 Х1
2 х4 -О X 2 Х2
х1 Х1
(2.3)
2
^ = [0 -2х10{(хз0 )2 + (х40 )2} 0 0]т. (2.4)
В дискретном случае система (2.2) будет иметь следующий вид:
х(п +1) = АД (п)х(п) + F(п), (2.5)
где индекс п — номер такта работы дискретной системы, к = 0,2 с — период квантования (длительность такта вычислительной бортовой ЭВМ),
0
АД (п) =
1 н
н(х2 + Х4) 1
2 х3 х2 н 2 н 2 Х
х1 х1
2 х4 х2 н 2 н 2 х1
Х12 х1
2 Хо
н
3 н 1
х1
4 н -он 1 -
0
2 Нх 4 Х1
о н
2 Х2
н
(2.6)
^ (п) = [0 -2 х1к{(х3 )2 + (х4)2} 0 0]т
(2.7)
3. Алгоритм точного размещения полюсов при решении задач наблюдения и идентификации. Далее будем рассматривать динамический объект, описываемый уравнениями вида
х(п +1) = А(п)х(п) + Г(п) , у(п) = Сх(п).
(3.1)
Наблюдатель полного ранга для системы (3.1) определяется уравнением
х(п +1) = (А(п) - ЬС)х(п) + Ьу(п) + п). Объединив (3.1) и (3.2), получим уравнение невязок х(п +1) = (Ль (п) - ЬС)х(п).
(3.2)
(3.3)
В соответствии со сформулированной задачей известными параметрами измерения являются хь х2 и неизвестными — х3, х4. Следовательно, матрица С имеет вид
г 1 0 0 0^ .0 0 0 0^'
с =
Для получения матрицы Ль для системы (2.5) необходимо разложить х(п + 1) и ^(п) в ряд Тейлора в окрестности точки X(п):
х2(п +1) = х2(п +1) + (2ЛХ3х1 + 4кх3х1 )х3 + (2кх4х1 + 4кх4х1 )Х4, (3.4)
п
С 2 х х х ^
х3 (п +1) = Х3 (п +1) + 1 —- к —- к + [1 —- к] \ Х3 + ОхкХ4 — Х3 + ОхкХ4,
V Х
С 2 х Х Х ^
х4 (п +1) = Х4 (п +1) - Ох кХ3 +1 —- к —2 к + [1 —- к] \ Х4 — -Ох кХ3 + Х4,
V Х1 Х1 Х1 /
Р — Р 4кх! Х3 Х3 4кх^ Х4.
Окончательно матрица для системы невязок примет вид
(3.4)
А =
1 к 0 0
к( Х32 + х4) 1 2 кХХз Х1 2 кХ4 Х1
2 Хз Х2 7 —^^ к 2 п Х1 2 ХХ3 х1 к 1 " X к
2 Х4 Х2 к 2 к Х1 2 Х4 Х1 к X к 1
(3.5)
Выбором матрицы коэффициентов Ь при известных матрицах Ль и С всегда можно обеспечить любое заданное размещение на комплексной плоскости корней характеристического полинома
ьщап - (А - ьс)).
Для решения задачи оценки угловой скорости линии визирования воспользуемся следующим алгоритмом синтеза наблюдателя состояния полного ранга [2].
1. Задать матрицы
2. А = Ат, В0 = ст. 2. Вычислить
N = сеП
V т )
(3.6)
3. Задать матрицы Ф = Ф0, Ф1, ..Ф^ такие, что
N+1
и eig(Фг_1) — желаемый спектр наблюдателя состояния.
I=1
4. Вычислить ортогональный аннулятор Вк_х, а затем матрицы
А - В1 А В1 т
Ак - Вк-1 Ак-1 Вк-1,
Вк - Вк-1 Ак-1Вк-1,
(3.7)
к -1, N.
5. Последовательно вычислить матрицы
Ц =Ф N - К АИ,
Вк = Вк - Цк+1 Вк , (з 8)
Ц =ФкВ" - В" Ак, '
к = N -1,0.
Здесь В0, ..., ВN — псевдообратные матрицы Мура — Пенроуза.
В случае если некоторые матрицы Вк (3.7) не являются матрицами
полного ранга, воспользоваться алгоритмом напрямую нельзя. Тогда необходимо воспользоваться скелетным разложением матрицы Вк:
Вк = Вк Т. (3.9)
Если приведенный алгоритм «перезапустить» для пары матриц
(Ак,к-Ь Вк,к-\)
Ак, к-1 = Ак , Вк, к-1 = Вк ,
то получим новый подуровень декомпозиции
(3.10)
А - Я1 А Я1Т
Ак, к - Як, к-1 АкЯк, к-1,
Як, к - Як, к-1 Ак Як, к-1, где ВКк * 0.
Известно, что регулятор Рк, обеспечивающий заданное размещение полюсов управляемой пары матриц (Ак, Вк), т. е.
Ак + ВкРк) с С81аЬ,
дает в преобразованном виде Тк+ Рк такое же размещение полюсов исходной пары матриц (Ак, Вк):
е^ Ак + Вк Тк+ Рк) = е^ Ак + Вк Рк) с С81аЬ. (3.11)
4. Применение алгоритма точного размещения полюсов при решении задачи оценки угловой скорости линии визирования в процессе сближения космических аппаратов. Применим изложенный выше подход к решению задачи оценки угловой скорости линии визирования при сближении КА.
В соответствии с (3.10) имеем:
"1 а21 а31 а41 "1 0"
Ат = к 1 а32 а42 , ст = 0 1
0 а23 1 + а33 - а34 0 0
0 а24 а34 1 + а33 0 0_
(4.1)
В данном случае размерность подпространств состояний, описывающих систему, — п = 4, наблюдаемых переменных — т = 2. Число уровней декомпозиции для каждого канала из выражения (3.6)
N = сеП
п
— I-1 = 2 -1 = 1
V т )
— два (нулевой и первый).
Согласно введенной многоуровневой декомпозиции, нулевой уровень для системы (3.3) с матрицами (4.1) имеет вид:
А = Ат, В0 = ст.
С учетом (4.1) имеем:
В = (Ст )х =
0 0 10 0 0 0 1
Во = (Ст )+ =
10 0 0 0 10 0
(4.2)
(4.3)
Первый уровень выглядит следующим образом:
А = В1АВо1т = (С1 г А (Сí)- =
тч 1 АТ//-Т\±т
1 + а
33
-а
34
а-
34
1 + а.
33
(4.4)
В, = В0 А В0 = (Ст )х Атст =
0 а 0 а
23
24
Анализ матрицы Б\ показывает наличие у нее нарушения полноты ранга по столбцам и необходимости выполнения «скелетного разложения». В соответствии с выражением (3.9) определим
В = Ъ{Т, Ъ =
23
24
, т = [ 0 1].
(4.5)
Тогда:
Ь = 1] к =
23
а
23
а.
24
а2з + а24
а23 + а24
(4.6)
В соответствии с (3.10) для пары матриц (Ьг, А1) получим новый подуровень декомпозиции:
«1,1 = ь1о Аь11 = Ь11 АА1Т =
—24 (1 + а33 ) + а
V а23
34
а23а24
а
24 а34 +1 + а33
+
V а23
23
у
(4.7)
а23 + а24
ь1,1 = bi, оAibi , о = bi Aibi = -—(1 + а33) + а
2 , 2 - а24
Л f
а23 + а24
V а23
34
а23 +
/
— а34 +1 + а33
V а23
а24-
Зададим матрицы Ф0, Ф1 в следующем виде:
фо =
/01 0 0 /о2
, Ф1 =
/11 0 ■ /12.
(4.8)
Выполняя вычисления по формулам (3.8) с учетом (4.2) - (4.8), получим
lt _ fn - а1,1 _ ( - f12 + 1 + а33 )а23 L1,1 ~ '
b
1,1
a34 (a23 + a24 )
V = bi+ - AVi1 »
Z = fnbb - КAx.
Принимая во внимание (3.11) и (4.5), получим:
Ат = [о 1]+ £,
B0 = B0 + - V >
LT =Ф 0 B0~ - B0- ( AT )0 =
/п /j2 Zj3 /14
/21 /22 /23 /24
(4.9)
где:
'11 = foi -1,
'12 = —a21'
'13 = - a31'
'14 = —a41,
'21 = -h,
122 _ ^02 + ^11 + ^12 2а33 3'
-а24 (/о2 + /и + ./12 - /02/11 - /й^ - /1/2 + ^^^ - 1
123 _ '
а34 (а23 + а24 )
+
+ (а24 а34 а24а33 + 2а23а33а34 )(/02 + /11 + /12 3) + а34 (а23 + а24 )
34 а24а33 )(^./02 + 2/11 + 2/12 - "/>2/1 /02/2 - - 3)
+
+
а34 (а23 + а24 )
+
3 3 2 2 2 2
^24^33 + ^23 ^34 ^32^34^23 ^32^34 ^24 3^23^34^33 3^24 ^33^34
а34 (а23 + а24 )
а23 (/02 + /11 + /12 - /^И - ./2/2 - /1/2 + /02/ 1/2 - 1
114 _
а34 (а23 + а24 )
+
+ (а23а33 - а23а34 + 2а24а33а34 )(/02 + /1 + /2 - 3) + а34 (а23 + а24 )
(а23а33 + а24а34)(^./02 + 2/11 + 2/12 - /2/1 - /02/2 - /1/2 - 3)
+
+
а34 (а23 + а24 )
+
3 3 2 2 2 2
^24^34 ^23^33 ^34^42 ^23 ^34 ^42 ^24 + 3^23^33^34 3^24 ^34 ^33
а34 (а23 + а24 )
Соответственно матричное произведение ЬС, входящее в (3.3), для каналов управления будет равно:
/11 /21 0 0
^12 /22 0 0
/13 /23 0 0
/14 /24 0 0
¿С =
Уравнения (3.2) можно записать в следующем виде: Х1 (п + 1) = х1 (п) + (п) + /1181 (п + 1) + /2182 (п + 1), Х2 (п +1) = Н( Х32 + Х4) х (п) + х2 (п) + /12а1 (п +1) + 122а2 (п +1),
(4.10)
Х3 (п +1) =
г 2 .с Х ^ 1 - Ь
V Л1 У Х4 (п +1) = -Юж Х3 (п) +
Х3 (п) + ШхХ4(п) + /13а1 (п +1) + /23а2(п +1), (4.11) Х4 (п) + /14а1 (п +1) + /24а2 (п +1),
г 2 .с Х ^
1 - ±ХЛ А
V
Ч У
где
82 (п) = х1 (п) - х (п), а2 (п) = х2 (п) - Х2 (п) . (4.12)
Рассмотрим несколько численных примеров оценки угловой скорости линии визирования по результатам измерения дальности и скорости продольного движения для различных наборов начальных значений вектора состояния активного КА относительно пассивного. Величина &2у + &2 представлена на конец пятой вычислительной итерации. Результаты сведем в таблицу.
№ В, м В, м/с О , У 1/с О г, 1/с Й У' 1/с 6 г, 1/с +"2; 10-5,1/с2 6 2+6 2, 10-51/с2
1. 200 -1,0 0,006 0,008 0,005 0,005 10,040 9,974
2. 150 -0,7 -0,003 0,007 0,005 0,005 5,821 5,794
3. 100 -0,5 0,003 -0,004 0,005 0,005 2,510 2,524
4. 50 -0,2 -0,009 -0,003 0,005 0,005 9,028 8,977
Изменение величины отклонения квадрата угловой скорости линии визирования за пять вычислительных итераций представлено на рис. 2.
хю-5
Рис. 2. Отклонение квадрата угловой скорости линии визирования
Следовательно, за пять тактов работы бортового алгоритма обеспечивается сходимость итерационного процесса оценки вектора угловой скорости линии визирования. Учитывая, что решение (4.9) представляет собой аналитические выражения, вычислительные затраты незначительны, и он вполне может быть реализован в реальном масштабе времени.
Заключение. В статье рассмотрена задача восстановления вектора состояния при сближении космических аппаратов при наличии отказа в канале измерения угловой скорости линии визирования. Раз-
работан бортовой алгоритм оценки и получено аналитическое решение для этого алгоритма.
ЛИТЕРАТУРА
[1] Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. Определения, теоремы, формулы. Москва, Наука, 1973.
[2] Зубов Н.Е., Микрин Е.А., Мисриханов М.Ш., Рябченко В.Н. Модификация метода точного размещения полюсов и его применение в задачах управления движением КА. Изв. РАН. ТиСУ, 2013, № 2, с .148-162.
Статья поступила в редакцию 28.06.2013
Ссылку на эту статью просим оформлять следующим образом: Зубов Н.Е., Микрин Е.А., Олейник А.С. Оценка угловой скорости линии визирования в процессе сближения космических аппаратов по результатам измерения дальности и скорости продольного движения. Инженерный журнал: наука и инновации, 2013, вып. 10. URL: http://engjournal.ru/catalog/it/nav/1079.html
Зубов Николай Евгеньевич — д-р техн. наук, заместитель руководителя по науке научно-технического центра ОАО «РКК "Энергия" имени С.П. Королёва», профессор кафедры «Системы автоматического управления» МГТУ им. Н.Э. Баумана. Автор более 70 работ в области проблем управления космических аппаратов. e-mail: [email protected]
Микрин Евгений Анатольевич — д-р техн. наук, академик РАН, первый заместитель генерального конструктора ОАО «РКК "Энергия" имени С.П. Королёва», заведующий кафедрой «Системы автоматического управления» МГТУ им. Н.Э. Баумана. Автор более 100 работ в области проблем управления космических аппаратов.
Олейник Алексей Сергеевич — аспирант ОАО «РКК "Энергия" имени С.П. Королёва». Автор 2 работ в области проблем управления космических аппаратов.