УДК 531.3:681.5.01
МЕТОД ПОСЛЕДОВАТЕЛЬНОГО ЗАМЫКАНИЯ
МОД ДВИЖЕНИЯ ДЛЯ МНОГОМЕРНЫХ, МНОГОСВЯЗНЫХ
ДИНАМИЧЕСКИХ СИСТЕМ
С.Н. Тимаков, К.А. Богданов, С.Е. Нефедов
ОАО "Ракетно-космическая корпорация "Энергия" им. С.П. Королева",
Королев, Московская обл., Российская Федерация
e-mail: [email protected]; [email protected]; [email protected]
В настоящее время одной из наиболее быстро развивающихся областей прикладной теории управления является разработка аналитических и численных методов модального управления корнями многомерных динамических систем посредством замыкания многосвязными обратными связями. Интерес к проблеме модального размещения полюсов замыкаемой динамической системы и вычислению компонент матриц обратной связи, а также матриц коэффициентов обратной связи не ослабевает уже на протяжении нескольких десятилетий. Предложен метод модального управления корнями характеристического полинома многосвязной системы, основанный на принципе последовательного замыкания. Приведено детальное описание алгоритма замыкания и вычисления как матриц обратной связи, так и матриц весовых коэффициентов для задач построения управления при неполном составе измерений. Работоспособность алгоритма продемонстрирована на задаче синтеза алгоритмов управления движением космического аппарата с двойным вращением (вращающийся солнечный парус с компенсирующим гироскопом) и задаче поиска и поддержания положения равновесия МКС.
Ключевые слова: проблема размещения полюсов, многомерная многосвязная динамическая система, синтез алгоритмов управления движением космических аппаратов.
SEQUENTIAL CLOSING METHOD FOR MODES OF MOVEMENT FOR MULTIDIMENSIONAL, MULTILINKED DYNAMICAL SYSTEMS
S.N. Timakov, K.A. Bogdanov, S.E. Nefedov
OAO "Korolev Rocket and Space Corporation "Energiya",
Korolev, Moscow region, Russian Federation
e-mail: [email protected]; [email protected]; [email protected]
Currently one of the fastest growing areas of applied control theory is the development of analytical and numerical methods of modal control by roots of multidimensional dynamical systems by means close loop of multilinked feedbacks. Interest to the problem of modal placement poles of close loop dynamical system and to the calculation of the components of feedback matrix as well as the coefficients of feedback matrix hasn't been decreasing already for several decades. In this paper a method of modal control by roots of the characteristic polynomial of multilinked systems, based on the principle of sequential closing is proposed. Detailed description of the closing algorithm and the calculating both feedback matrices and matrices of weight coefficients for problems of control definition with incomplete composition of measurements is carried out. Efficiency of the algorithm has been demonstrated using the problem of synthesis of algorithms control of spacecraft motion with double rotation (rolling solar sail with compensating gyro) and the other problem ofseeking and maintaining the equilibrium attitude position ofthe International Space Station.
Keywords: poles placement problem, multidimensional and multilinked dynamical
system, synthesis of algorithms control of spacecraft motion.
Большой интерес в мире в последнее время приобрели методы модального синтеза многомерных многосвязных динамических систем. Практическая значимость этих методов для разработки алгоритмов управления сложными системами общеизвестна. Далеко не полный список методов модального синтеза включает следующие методы: Ackermann (Аккерманн) [1], Mayne - Murdoch (Мейн - Мардох) [2], Maki-Van de Vegte (Маки-Ван де Вейт) [3], Barnett (Барнетт) [4], Gourinshankar-Ramar (Гоуришанкар - Ремер) [5], Moore (Мур) [6], Klein-Moore (Клейн-Мур) [7], Porter-D'Azzo (Портер-Д'Аццо) [8], Wonham (Уонем) [9], Munro (Мунро, 1979), Flamm (Флемм) [10], Varga (Варга) [11], Fahmy-O'Reilly (Феми и О'Рейли) [12], Kautsky-Nichols-Van Dooren (Каутский - Никольс - Ван Доорен) [13].
Перечисленные методы используют матричные преобразования, которые применяются к матрицам большого размера. В связи с этим с увеличением размерности динамической системы резко возрастает погрешность вычислений компонент матриц обратной связи и, как следствие, не достигается желаемое расположение корней характеристического полинома замыкаемой системы. Предлагаемый метод замыкания использует принцип последовательного замыкания системы, что позволяет применять процедуры матричных преобразований к выделяемым матричным блокам размером не более 2x2. При таком подходе в силу ортогональности модального разложения ошибка численных преобразований матриц почти не накапливается с увеличением размера. В результате матрицы обратной связи и матрицы весовых коэффициентов достаточно точно переводят корни характеристического полинома в желаемое место в левой полуплоскости s-плоскости для континуальных систем и во внутреннюю область единичной окружности z-плоскости для дискретных систем.
Описание метода последовательных замыканий мод движения. Описываемый метод позволяет по заданному эталонному расположению корней замкнутой системы осуществлять поиск как коэффициентов матрицы обратной связи (т.е. построение регулятора), так и поиск матрицы весовых коэффициентов (построение наблюдателя). Главная особенность метода — способность работать с системами multi-input multi-output, т.е. с системами, имеющими несколько входных (каналов управления) и выходных (измерителей) параметров. Суть данного метода заключается в последовательном замыкании мод движения динамической системы, т.е. сначала осуществляется поиск обратной связи, замыкающей первую пару корней (собственных чисел), затем поиск обратной связи, замыкающей вторую пару, и так далее, пока все моды не будут замкнуты. Наглядно это представлено на рис. 1.
Рис. 1. Графическая иллюстрация метода последовательных замыканий:
0 — эталонное расположение корней замкнутой системы; + — расположение корней матрицы переходных состояний на г-м шаге
Подробно опишем метод последовательного замыкания на примере построения регулятора для линейной динамической системы. Рассмотрим систему, описывающуюся следующим уравнением:
x = Ax + Bu, (1)
где x — n-мерный вектор состояния динамической системы (1); u — m-мерный вектор управления; A — матрица системы размером n x n; B — матрица управления размером n x m.
Пусть данная система является полностью управляемой по критерию Калмана, т.е. ранг матрицы управляемости будет равен размерности вектора состояния [14]:
rank ([ B AB A2B ... An-1 B ]) = n.
Зададим эталонное расположение корней характеристического полинома замкнутой системы
(Ai ,А2 ,...,А; } . (2)
Цель — поиск такой обратной связи u = —Dx, которая бы обеспечивала характеристическому полиному | A — BD — AI| =0 замкнутой системы эталонное расположение корней (2).
Подробнее рассмотрим замыкание первой моды движения системы (1).
Этап 1. Находим преобразование подобия Т, отображающее матрицу А в блочно-диагональную матрицу собственных чисел
Л = Т-1АТ,
где
Л=
а1 Ь1
0 0
0 0
-Ь1 ai... 0 0 ... 0 0
Л=
a1 Ь1 1 0 . . . 0 0..
-Ь1 a1 0 1 . . . 0 0..
00 a1 Ь1 . . . 0 0..
00 -Ь1 a1 . . . 0 0..
00 0 0 . . аг 1 .. ..
00 0 0 . . . 0 аг..
Здесь а и Ь — действительная и комплексная части г-й пары собственных чисел; и — пара действительных собственных чисел (а - вещественный аналог диагональной матрицы). В случае, если матрица А не является нормальной и имеет корни кратности два и более, матрица Л (3) берется в виде вещественного аналога жордановой формы (б).
Зная эталонное расположение корней (2) составляем блочно-диагональную матрицу Л*, аналогичную матрице Л.
Л* =
ai Ь1 -Ь1 ai
0 0
0 0
0 0
0 0
0 0
-Ь*
0 0
0 0
0 0
-1-1
0 0
0 0
Ъ'1
0 0
0 0
j2
(3)
Этап 2. Выделяя из матрицы Т— две первые строки, получаем матрицу Т1. Из матриц (3) и (4) выделяем клетки размером 2x2, содержащие соответственно первую передвигаемую и первую желаемую (эталонную) пару собственных чисел, обозначим их как Л1 и ЛЦ.
Этап 3. Вычисляем матрицу В1 = Т1В, где В — матрица управления, затем находим ее псевдообратную матрицу Мура-Пенроуза:
В-1 = (В1 ВТ)-1В[. (4)
a
Ь
0
Находим матрицу
AAi = Ai - Л1. (5)
Этап 4. Перемножая (5), (6) и Ть получаем
D1 = B-1 х АЛ1 х Т1. (6)
Далее находим новую матрицу A1 = A — BD1.
Затем проделываем все упомянутые ранее действия (этапы 1-4) только уже с парой (A1, B), выделяя на этапе 2 уже вторую пару строк из матрицы Т-1 и вторую пару передвигаемых и желаемых пар собственных чисел. Таким образом, находим матрицу D2 и матрицу переходных состояний A2 = A1 — BD2. На следующем шаге находим матрицы D3 и A3 на i-м шаге матриц Dj и Aj. Данный процесс заканчивается на шаге n/2, когда будет "передвинута" в желаемое место последняя пара собственных чисел (если n — четное), или на шаге n/2 + 1, когда будет "передвинуто" последнее непарное собственное число (если n — нечетное).
Суммируя матрицы обратной связи, найденные на каждом шаге, получаем искомую матрицу обратной связи
n/2
D = Е d
i=1
Задача поиска матриц весовых коэффициентов для адаптивного наблюдателя сводится к предыдущей задаче путем составления сопряженной системы. Пусть пара матриц (A, C), где A — матрица системы, C — матрица измерений, удовлетворяет критерию полной наблюдаемости Калмана
rank ([ CT ATCT (AT)2C ... (AT)n-1C ])T = n.
Используя принцип "дальности" построения наблюдателя и регулятора [14], составим сопряженную систему
x' = AT x' + CTu',
где AT и CT — "новые" матрицы переходных состояний и управления. В соответствии с описанным выше алгоритмом находим матрицу W' = WT, где W — искомая матрица весовых коэффициентов.
Построение адаптивного наблюдателя и регулятора для космической платформы с вращающимся солнечным парусом. Описание конструкции космической платформы с солнечным парусом. Космическая платформа с солнечным парусом представляет собой собственно вращающийся солнечный парус в виде пленочного диска с центральной жесткой вставкой, приборный отсек и компенсирующий силовой гироскоп в подвесе Гука (внутреннем кардановом подвесе).
Внутренний карданов подвес обладает управляемыми и контролируемыми углами поворота, что позволяет отклонять ось вращения ротора гироскопа от оси вращения центральной жесткой вставки паруса для создания управляющего гироскопического момента. Центральная вставка выполнена в виде вантовой конструкции и служит для передачи момента импульса солнечного паруса приборному отсеку [15].
Введем систему координат OXYZ так, чтобы оси координат были связаны с осями чувствительности датчиковой аппаратуры (рис.2). Ось ОХ направим в сторону вращения ротора силового гироскопа, ось OY будет лежать в плоскости вращения центральной вставки, а ось OZ будет дополнять систему координат до правой тройки [15].
Уравнения движения. Для описания динамики объекта управления разобьем его на два тела: первое — солнечный парус; второе — приборный отсек с силовым гироскопом в подвесе Гука (см. рис. 2).
Рассмотрим первое тело. Ранее было доказано, что 99,9 % массы мембранного диска солнечного паруса совершают колебания на первых двух гироскопически связанных модах движения, поэтому динамическое поведение солнечного паруса с большой степени точности можно представить динамикой одного гироскопа в упругом подвесе [15]. Исходя из этого, кинетический момент первого тела будет выглядеть следующим образом:
Ъ1 = Мт 31 (П + Ми),
где M =
1 Ду
ßz 1 -Дх Ду Дх 1
матрица малого поворота векто-
ра угловой скорости; Ji =
A 0 0 0 C 0 0 0 C
приведенный момент
инерции паруса (A = 2C); ш =
Ф x
Ф У
Ф z
— угловая скорость КА;
П =
руса [
-П + ßx ß У ß z
6].
относительная угловая скорость вращения па-
Кинетический момент второго тела равен
h2 = J2w + BH,
Jx 0 0
где J2 = 0 Jy 0 — момент инерции КА (примем для упроще
0 0 Jz
ния дальнейших выкладок, что Jy = Jz = J ); B =
1 -ßz ßy ßz 1 -ßx ßy ßx 1
матрица малого поворота вектора угловой скорости ротора гироско-0
па; H =
0 H
кинетический момент ротора силового гироскопа
в связанной с ним системе координат (для дальнейших вычислений примем, что H = AQ согласно закону сохранения кинетического момента) [16].
Применяя теорему об изменении кинетического момента ко всему объекту управления и отдельно к парусу и пренебрегая моментами сил солнечного давления, получаем следующую систему уравнений, описывающую динамику движения КА:
h1 + h2 = const;
h 1 + ш х h1 = —k2Cp,,
где коэффициент жесткости центральной жесткой вставки k считается известным и вычисляется по формуле
(7)
k2
(3 + ß)2a2 о2
П2 и 0,01 c
2
2(1 + ^)Я2
^ = 0,4 — коэффициент Пуассона; а и Я — внутренний и внешний радиусы пленочного диска солнечного паруса [16].
Расписывая систему (1) покомпонентно, проведя линеаризацию с точностью до второго порядка малости, а также полагая, что вокруг оси ОХ система управления достаточно точно удерживает аппарат,
получаем уравнения движения КА вокруг осей OY и OZ [15]:
ПД, + Ду + (1 + С )( + 2П& = 0;
—Пру + Д + (1 + )( - 2П/Зу = 0; (8)
2ПДг + ру г + (у — 2П(г + к2/у = 0; -2П/лУ + г + (¿г + 2ПшУ + к2/ = 0.
Описание объекта в пространстве состояний. Опишем объект управления в пространстве состояний. Для этого введем восьмикомпо-нентный вектор состояния хр = ( шу /у р,у / рг ) ,
содержащий углы и угловые скорости вращения аппарата, а также углы и угловые скорости отклонения плоскости мембранного
диска от положения равновесия. В качестве вектора управления • • т
ир = ( ¡Зу Дг ) примем угловые скорости отклонения ротора гироскопа от осей OY и OZ.
Таким образом, система (2) в пространстве состояний будет описываться следующим уравнением:
хр = Архр + Врир,
где
Ар
' 0 1 0 0 0 0 0 0
0 0 0 —2C Ü/J Ck2/J 0 0 0
0 0 0 1 0 0 0 0
0 2C Ü/J 0 0 0 0 Ck2/J 0
0 0 0 0 0 1 0 0
0 0 0 2Ü(C + J )/J —(C + J )k2/J 0 0 2Ü
0 0 0 0 0 0 0 1
0 —2Ü(C + J )/J 0 0 0 —2Ü —(C + J)k2/J 0
матрица системы;
B„ =
0 0 0
2C П/2 0 0 0
0
-2C П/2 0 0 0
2C П/2 0
-2CQ/2 0
матрица управления.
Поскольку не все параметры вектора состояния доступны прямому измерению, возникает необходимость введения адаптивного наблюдателя [17].
Хр = АрХр + Брир - Wp(yp - СрХр); (9)
Ур = СрХр, (10)
вектор измерений;
где Хр — вектор оцениваемых параметров; ур Ср — матрица измерений;
Wp =
т
Cp —
Wll Wl2 Wlз Wl4 Wl5 Wl6 Wl7 Wl8 W21 W22 W23 W24 W25 W26 W27 W28 Wзl Wз2 Wзз Wз4 Wз5 Wз6 Wз7 Wз8
w41 w42 w43 w44 w45 w46 w47 w48
— матрица весовых коэффициентов, обеспечивающая асимптотическую сходимость вектора оцениваемых параметров Хр к вектору состояния Хр.
Объект управления имеет четыре датчика для измерения компонент углового положения и угловой скорости КА относительно осей OY и OZ. Исходя из этого, матрица измерений имеет следующий вид:
1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0
0 0 1 0 0 0 0 0
0 0 0 1 0 0 0 0
Введем обратную связь по оценке вектора состояния:
ир ОрХр
dll dl2 dlз dl4 dl5 dl6 dп dl8 d21 d22 d23 d24 d25 d26 d2r d28
ной связи, обеспечивающая асимптотическую устойчивость параметров объекта управления.
Введем вектор невязок Хр = Хр—Хр как разность вектора состояния и вектора оценок.
Объединяя (3)—(6), выражая все через Хр и ХХр, получаем следующую систему уравнений, описывающих объект управления с учетом введения обратной связи и адаптивного наблюдателя:
Хр = (Ар БрОр) Хр + БрОрХр; Хр = (Ар Ср^^р) Хр.
где Dp —
(11)
— матрица обрат-
(12)
Объединим восьмикомпонентный вектор невязок и восьмикомпо-нентный вектор состояния в один новый 16-компонентный вектор со-
стояния xe —
Xp
. Тогда систему (6) можно записать в виде следу-
ющего однородного матричного дифференциального уравнения:
Xe —
Ap BpDp BpDp
0
A _ w C
JT\.p »» p^p
(13)
Для обеспечения асимптотической устойчивости системы (7) необходимо, чтобы все собственные числа матрицы
Ap BpDp BpDp
0
Ap ^Vp Cp
лежали в левой полуплоскости (система описана в континуальном пространстве). Характеристический полином системы (14) выглядит следующим образом:
Is - (Ap - BpDp)| • |Is - (Ap - CpWp)
(14)
где I — единичная матрица размером 8x8.
Из характеристического полинома (8) следует, что расположение корней адаптивного наблюдателя не зависит от расположения корней замкнутой системы, в состав которой он включен. Поэтому можно отдельно выполнить построение как адаптивного наблюдателя (поиск матрицы весовых коэффициентов Wp), так и регулятора (поиск числовых значений матрицы обратной связи Юр) [14].
Поиск матрицы весовых коэффициентов, матрицы обратной связи. Результаты математического моделирования. В качестве эталонного полинома для построения наблюдателя и регулятора возьмем полиномы Баттерворта 8-го порядка [18] с частотой среза шс, равной 1,5 и 0,5рад/с соответственно (рис.3).
88 + 5,1шс/ + 13,1шс286 + 21,8^4
+ 25,7^4^4 + 21,8шс58з + 13,1^6«2 + 5,1^7^ + ^8.
Для поиска матрицы обратной связи Юр и матрицы весовых коэффициентов Wp используется метод последовательных замыканий мод движения. В результате получаем следующие числовые значения коэффициентов матриц Юр и Wp:
Dp
0,0096 0,9140 0,3071 -0,9344 0,0228 -0,0838 0,3134 -0,9597 -0,0219 0,32910,1050 0,4615 -0,0234 0,3410 0,1097 -0,5432
Рис. 3. Расположение корней эталонного полинома Баттерворта на комплексной плоскости для построения адаптивного наблюдателя (а) и регулятора (б)
Wp =
0,7508 0,9949 -0,0664 -0,0021
-0,6899 0,5012 -0,1119 -68,1355
0,1131 0,0019 0,7876 0,9940
0,0485 68,1753 -0,6704 0,5233
-16,4562 -1,3974 -15,2103 0,5080
17,9184 -0,9926 27,1065 -67,7065 16,9256 -0,4940 -15,8974 -1,3921
26,7090 -67,7344 18,2892 -1,0343
Подставляя полученные числовые значения матриц в математическую модель объекта управления, проведя математическое моделирование режима активного демпфирования колебаний пленочного диска солнечного паруса, получаем следующие зависимости основных параметров вектора состояния (рис. 4-6).
Выведение и удержание МКС в положении динамического равновесия. Описание системы координат. Уравнения движения. В качестве системы координат используется система LVLH (Local Vertical -Local Horizontal). В данной системе началом координат является центр масс объекта управления, ось OX направлена вдоль вектора орбитальной скорости, ось OZ против радиуса-вектора из центра Земли в центр масс объекта управления, ось OY дополняет систему координат до правой тройки (рис. 7). В качестве связанной системы координат будем использовать систему ISS, оси которой совпадают с основными осями инерции объекта.
Динамика объекта управления, в данном случае МКС, описывается следующими уравнениями.
Рис. 4. Компоненты углового положения и угловой скорости КА вокруг осей ОУ и ОЕ
Рис. 5. Угловые отклонение и скорость плоскости вращения солнечного паруса относительно осей ОУ и ОЕ
Рис. 6. Скорость прецессии ротора силового гироскопа относительно осей ОУ и ОЕ
1. Динамическим уравнением углового движения КА при воздействии гравитационного момента:
п r Г
JÜabs + Üabs X Jüabs + Ii + Üabs X h = 3—- X J-,
r 3 r
(15)
где Л — тензор инерции объекта управления; Ь — вектор суммарно-
^ г тг
го кинетического момента маховиков; 3—- хЛ- — гравитационный
момент;--единичный вектор местной вертикали [17].
r
Рис. 7. Система координат LVLH
2. Кинематическими соотношениями в кватернионном виде:
♦ е
214. = ^ О (и^)е - (и0). ◦ N1®, (16)
где N1 — кватернион разворота из базиса 1 в базис е (1 — орбитальный базис, е — связанная система координат, оси которой расположены вдоль главных осей инерции КА), (иаЬз)е = [ 0 шу иг ]Т — вектор абсолютной угловой скорости в проекциях на связанный базис, (^0). = [000 ]Т — вектор угловой скорости приборного базиса в проекциях на орбитальную систему координат [19].
Проэцируя (1) и (2) на связанный базис, применяя теорему об изменении кинетического момента отдельно к корпусу МКС и отдельно к маховикам, получаем систему из шести уравнений, описывающих динамику объекта управления [20]:
З*7 = (Л - Зу)7 + ^о(З* - Зу + Л- их; ЛV3 = (З* - Зу- ^о(Зх - Зу + Зг)7 - иг;
Зу<33 = 3^0З - З*- иу;
к * = К г ^0 + и*; Кг = -К* ^0 + иг;
ку иу 5
где 7, V, ^ — соответственно углы разворота по каналам крена, рысканья и тангажа: З*, Зу, Зг — моменты инерции МКС относительно
главных осей инерции; — абсолютная угловая скорость МКС; Нх, Ну, Нг — компоненты вектора суммарного кинетического момента маховиков.
Описание объекта управления в пространстве состояний. Поскольку каналы крена и рысканья имеют гироскопическую связь, будем отдельно рассматривать динамику объекта управления по каналам крена и рысканья и отдельно по каналу тангажа.
Вектор состояния и вектор управления для каналов крена и рысканья и отдельно для канала тангажа соответственно будут иметь следующий вид:
х-гФ = 77 ь,х / нх(],г ф ф Нг / Нг¡1 ¡1 ¡2 ¡2 ¡з ¡з
Х —
и7ф = [их иг ]Т ;
ф ф Ну ! Ну (И я1 (¡1 ¡2 ¡¡2 ¡3 ¡3
и^ = [иу];
где [ их иу иг ]Т — момент реакции в подшипниках маховиков, через который осуществляется управляющее воздействие на корпус КА; ¡1, ¡2, ¡3 — компоненты, отвечающие за аэродинамический момент, в данном случае это три главных тона разложения функции плотности атмосферы по витку.
Таким образом, в пространстве состояний объект управления будет описываться следующими уравнениями:
Х-'УФ хгФ + и~{Ф;
где матрицы системы и матрицы управления для системы крен-рысканье и для системы тангажа выглядят соответственно следующим образом:
A
Y Ф =
" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Gx 0 0 0 0 Hx 0 0 0 0 0 0 0 0 -1/Jx 1 0
0 0 0 0 0 0 Wo 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
0 Hz 0 0 Gz 0 0 0 0 0 0 0 0 0 0 -1/Jz 1
0 0 -Wo 0 0 0 0 0 0 0 0 0 0 0 ; В7ф— 0
0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 1 0 0 0 -W2 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0 0 0 -4W2 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
0 0 0 0 1 0 0 0 0 0 0 0 -9w2 0 0 0
A —
0 10 0 0 0
My 0 0 0 0 0
0 000 0 0
0 0 10 0 0
0 000 0 1
1 0 0 0 —ug 0
0 000 0 0
0000 0000 0000 0000 0000 0000
01
00 00
1 0 0 0 0 0 —4ul 0 0 000 0 0 0 0 0 1 1 0 0 0 0 0 0 0 —9ul 0
0
— 1/jy 1 0 0 0 0 0 0 0
Здесь для упрощения записей введены следующие обозначения:
= 4^0 (Зг - Зу) /З*, Н* =^0 (З* - Зу + Зг) /З*, =^0(З* - Зу)/Зг,
Н = ^0 (З* - Зу + Зг) /Зг, Му = 3^0(Зг - З*)/Зу.
Поиск матриц обратной связи. Для придания асимптотической устойчивости объекту управления введем обратную связь:
и7ф Б7фХ7ф;
для систем крен-рысканье и тангажа соответственно. В результате уравнения, описывающие КА в пространстве состояний, приобретут следующий вид:
Х7 ф (А7 Ф Б7ф Б7ф )х7 ф;
^^^•
В качестве эталонного характеристического полинома для матрицы (А - ББ) возьмем полином Баттерворта с частотой среза 2^0 [18]. Очевидно, что для системы крен-рысканье полином Баттерворта будет иметь 14 порядок, для системы тангажа — 10 порядок (рис. 8).
Рис. 8. Расположение корней эталонного полинома Баттерворта для системы крен-рысканье (а) и для системы тангажа (б)
Используя метод последовательного замыкания, получаем следующие числовые значения для матриц обратной связи
т
-9,3782e + 03 5,1506e + 06 -4,6465e - 02 -1,6915e - 05 -6,4861e - 03 -8,8120e - 01 -2,1849e - 17 -1,5962e + 00 4,5627e - 03 -7,1139e + 00
D
ii>=
-9,2642e + 03 5,7158e + 03
2,2672e + 06 -5,8619e + 06
2,7193e - 02 -4,1871e - 02
4,6806e - 06 3,1351e - 05
2,9367e + 03 -2,8407e + 03
-1,1700e + 07 1,0268e + 07
-6,6631e - 02 6,4046e - 02
-3,0921e - 05 1,5060e - 05
3,5352e - 04 8,9431e - 03
1,8997e + 01 -8,8208e + 00
-3,4739e - 04 -1,9733e - 04
7,9223e - 02 -1,2335e - 01
-6,6505e - 03 3,2664e - 03
-5,2919e - 01 -1,7878e - 00
D,
• 106 кг-м2;
Результаты математического моделирования. Для математического моделирования поиска равновесной ориентации МКС в режиме разгрузки кинетического момента инерционных исполнительных органов (ИИО) с учетом аэродинамического момента понадобится знать значения следующих параметров:
— тензор инерции МКС
" 129,978287 3,294386 5,279790 J = 3,294386 83,966791 0,990955
5,279790 0,990955 193,699689 _
— модуль угловой скорости МКС = 0,0667 °/c;
— модуль скорости движения V = 8 • 103 м/с.
Координаты центра масс МКС в системе координат ISS — CMx = = -4,30 м, CMy = -0,54м, CMz = 3,92 м; коэффициент диффузии c = 0,35; средняя плотность атмосферы р = 6 • 10-13 кг/м3.
С учетом всех упомянутых параметров, а также найденных матриц обратных связей, определяющих закон управления КА, получаем следующие зависимости углового положения МКС, ее угловой скорости и кинетического момента (рис. 9-11).
Заключение. Описан модифицированный метод модального управления корнями характеристического полинома многомерной многосвязной динамической системы. Точность метода проиллюстрирована двумя примерами: на задаче синтеза закона управления для КА с двойным вращением, представляющего собой вращающийся солнечный парус с компенсирующим гироскопом в подвесе Гука, и на алгоритме
T
Рис.9. Угловое положение МКС по каналам крена (1), рысканья (2) и тангажа (3)
Витки
Рис. 10. Угловая скорость МКС по каналам крена (1), рысканья (2) и тангажа (3)
Витки
Рис. 11. Кинетический момент по каналам крена (1), рысканья (2) и тангажа (3)
управления МКС в режиме Momentum Management. Необходимо отметить, что процедура place, основанная на широко известном методе [13] и реализованная в среде MATLAB, со второй тестовой задачей не справляется из-за накопленных ошибок в процессе преобразований матриц большой размерности.
Работа выполнена при финансовой поддержке РФФИ (12-08-00254-а).
ЛИТЕРАТУРА
1. Ackermann J. Der entwurf linearer regelungssysteme im zustandsraum // Regelungstechnik und Prozessdatenverarbeitung. 1992. Vol. 7. P. 297-300.
2. Mayne D.Q., Murdoch P. Modal control of linear time invariant systems // Int. J. Control. 1970. Vol. 11. No. 2. P. 223-227.
3. Maki M.C., Van de Vegte / Optimization of multiple-input systems with assigned poles // I.E.E.E. Trans. autom. Control. 1974. Vol. 19. No. 2. P. 130-133.
4. Barnett S. Introduction to mathematical control theory (Oxford University Press). 1975. 264 p.
5. Gourishankar V., Ramar K. Pole assignment with minimum eigenvalue sensitivity to plant parameter variations // Int. J. Control. 1976. Vol. 23. No. 4. P. 493-504.
6. Moore B.C. On the flexibility offered by state feedback in multivariable systems beyond closed loop eigenvalue assignment // I.E.E.E. Trans. autom. Control. 1978. Vol. 21. No. 5. P. 689-692.
7. Klein G., Moore B.C. Eigenvalue-generalized eigenvector assignment with state feedback//I.E.E.E. Trans. autom. Control. 1977. Vol. 22. No. 1. P. 140-141.
8. Porter B., D'Azzo J.J. Closed-loop eigenstructure assignment by state feedback in multivariable linear systems // Int. J. Control. 1978. Vol. 27. No. 6. P. 487-492.
9. Wonham WMLinear Multivariable Control: A Geometric Approach, 2nd ed. N.Y.: Springer, 1979. 322 p.
10. Flamm D.S. A new proof of Rosenbrock's theorem on pole assignment // I.E.E.E. Trans. autom. Control. 1980. Vol. 25. No. 6. P. 1128-1133.
11. Varga A. ASchur method for pole assignment//I.E.E.E. Trans. autom. Control. 1981. Vol. 26. No. 2. P. 517-519.
12. Fahmy M.M., O'Reilly J. On eigenstructure assignment in linear multivariable systems // I.E.E.E. Trans. autom. Control. 1982. Vol. 27. No. 690.
13. Kautsky J., Nichols N.K., Van Dooren P. Robust pole assignment in linear state feedback, Internat. J. Control, 1985. Vol. 41. No. 5. P. 1129-1155.
14. Квакернаак Х., СиванР. Линейные оптимальные системы управления. М.: Мир, 1977. 650 p.
15. Легостаев В.П., Субботин А.В., Тимаков С.Н., Черемных Е.А. Собственные колебания вращающейся мембраны с центральной жесткой вставкой (применение функций Хойна) // Прикладная математика и механика. 2011. Т. 75. Вып. 2. С. 224-238.
16. Легостаев В.П., Субботин А.В., Тимаков С.Н., Зыков А.В. Исследование динамики управляемого углового движения космического аппарата с вращающимся солнечным парусом // Труды МФТИ. 2013. Т. 5, № 2. С. 106-119.
17. Белецкий В.В. Движение искусственного спутника относительно центра масс. М.: Наука, 1965. 416 c.
18. Кузовков Н.Т. Модальное управление и наблюдающие устройства. М.: Машиностроение, 1976. 184 c.
19. Бранец В.Н., Шмыглевский И.П. Применение кватернионов в задачах ориентации твердого тела. М.: Наука, 1973. 320 c.
20. Тимаков Н.С. Исследование управляемого углового движения космического аппарата на высокоэллиптической орбите // Навигация и управление движением. Матер. IX Конф. молодых ученых. СПб. 2007. С. 330-336.
REFERENCES
[1] Ackermann J. Der entwurf linearer regelungssysteme im zustandsraum. Regelungstechnik und Prozessdatenverarbeitung, 1972, vol. 7, pp. 297-300.
[2] Mayne D.Q., Murdoch P. Modal control of linear time invariant systems. Int. J. Control, 1970, vol. 11, no. 2, pp. 223-227.
[3] Maki M.C., Van de Vegte J. Optimization of multiple-input systems with assigned poles. I.E.E.E. Trans. autom. Control, 1974, vol. 19, no. 2, pp. 130-133.
[4] Barnett S. Introduction to Mathematical Control Theory. Oxford University Press, 1975, 264 p.
[5] Gourishankar V., Ramar K. Pole assignment with minimum eigenvalue sensitivity to plant parameter variations. Int. J. Control, 1976, vol. 23, no. 4, pp. 493-504.
[6] Moore B.C. On the flexibility offered by state feedback in multivariable systems beyond closed loop eigenvalue assignment. I.E.E.E. Trans. autom. Control, 1976, vol. 21, no. 5, pp. 689-692.
[7] Klein G., Moore B.C. Eigenvalue-generalized eigenvector assignment with state feedback. I.E.E.E. Trans. autom. Control, 1977, vol. 22, no. 1, pp. 140-141.
[8] Porter B., D'Azzo J.J. Closed-loop eigenstructure assignment by state feedback in multivariable linear systems. Int. J. Control, 1978, vol. 27, no. 6, pp. 487-492.
[9] Wonham W.M. Linear Multivariable Control: A Geometric Approach. 2nd ed. New York, Springer, 1979. 322 p.
[10] Flamm D.S. A new proof of Rosenbrock's theorem on pole assignment // I.E.E.E. Trans. autom. Control, 1980, vol. 25, p. 1128-1133.
[11] Varga A. A Schur method for pole assignment. I.E.E.E. Trans. autom. Control, 1981, vol. 26, no. 2, pp. 517-519.
[12] Fahmy M.M., O'Reilly J. On eigenstructure assignment in linear multivariable systems. I.E.E.E. Trans. Autom. Control, 1982, vol. 27. no. 3, pp. 690-693.
[13] Kautsky J., Nichols N. K., Van Dooren P.. Robust pole assignment in linear state feedback. Int. J. Control, 1985, vol. 41, no. 5, pp. 1129-1155.
[14] Kvakernaak Kh., Sivan R. Linear optimal control systems. Wiley-Interscience, New-York, 1972. (Russ. ed.: Kvakernaak Kh., Sivan R. Lineynye optimal'nye sistemy upravleniya. Moscow, Mir Publ., 1977. 650 p.)
[15] Legostaev V.P., Subbotin A.V., Timakov S.N., Cheremnykh E.A. Natural oscillations of a rotating membrane with a central rigid insert (use of Heun function) Prikladnaja matematika i mehanika [J. Appl. Math. Mech.], 2001, vol. 75, iss. 2, pp. 224-238 (in Russ.).
[16] Legostaev V.P., Subbotin A.V., Timakov S.N., Zykov A.V. Study of the dynamics of controlled angular motion of spacecraft with rotating solar sail. Trudy Moskovskogo fiziko-tekhn. inst. MIPT (SU) [Proc. of the Moscow Inst. Ph. Techn. (State University)], 2013, no. 2, pp. 106-119 (in Russ.).
[17] Beletskiy V.V. Dvizhenie iskusstvennogo sputnika otnositel'no tsentra mass. [Motion of satellite vehicle with respect to the center of mass]. Moscow, Nauka Publ., 1965. 416 p.
[18] Kuzovkov N.T. Modal'noe upravlenie i nablyudayushchie ustroystva. [Modal control and monitoring devices] Moscow, Mashinostroenie Publ., 1976, 184 p.
[19] Branets V.N., Shmyglevskiy I.P. Primenenie kvaternionov v zadachakh orientatsii tverdogo tela. [Applying of quaternions in problems of orientation of solid]. Moscow, Nauka Publ., 1973. 320 p.
[20] Timakov N.S. Study of controlled angular motion of spacecraft on highly elliptical orbit. Tr. IXKonf. molodykh uchenykh "Navigatsiya i upravlenie dvizheniem". [Proc. IX Conf. of young scientists "Navigation and motion control"], St. Petersburg, 2007, pp. 330-336 (in Russ).
Статья поступила в редакцию 21.02.2014
Сергей Николаевич Тимаков — канд. техн. наук, доцент МФТИ, ведущий научный сотрудник ОАО "Ракетно-космическая корпорация "Энергия" им. С.П. Королева". Автор более 25 научных работ в области систем управления.
ОАО "Ракетно-космическая корпорация "Энергия" им. С.П. Королева", Российская Федерация, 141070, Московская обл., Королев, ул.Ленина, д.4а.
S.N. Timakov — Cand. Sci. (Eng.), assoc. professor of the Moscow Institute for Physics and Technology, leading researcher of the OAO "Korolev Rocket and Space Corporation "Energiya". Author of more than 25 publications in the field of control systems. OAO "Korolev Rocket and Space Corporation "Energiya", ul. Lenina 4a, Korolev, Moscow region, 141070 Russian Federation.
Кирилл Андреевич Богданов — инженер ОАО "Ракетно-космическая корпорация "Энергия" им. С.П. Королева". Специализируется в области систем управления. ОАО "Ракетно-космическая корпорация "Энергия" им. С.П. Королева", Российская Федерация, 141070, Московская обл., Королев, ул.Ленина, д.4а.
K.A. Bogdanov — engineer of the OAO "Korolev Rocket and Space Corporation "Energiya". Specializes in the field of control systems.
OAO "Korolev Rocket and Space Corporation "Energiya", ul. Lenina 4a, Korolev, Moscow region, 141070 Russian Federation.
Сергей Евгеньевич Нефедов — инженер ОАО "Ракетно-космическая корпорация "Энергия" им. С.П. Королева". Специализируется в области систем управления. ОАО "Ракетно-космическая корпорация "Энергия" им. С.П. Королева", Российская Федерация, 141070, Московская обл., Королев, ул.Ленина, д.4а.
S.E. Nefedov — engineer of the OAO "Korolev Rocket and Space Corporation "Energiya". Specializes in the field of control systems.
OAO "Korolev Rocket and Space Corporation "Energiya", ul. Lenina 4a, Korolev, Moscow region, 141070 Russian Federation.