CC BY
57
УДК 629.7
ОПТИМИЗАЦИЯ УПРАВЛЕНИЯ ТЯГОЙ ДВИГАТЕЛЕЙ РАКЕТЫ-НОСИТЕЛЯ «СОЮЗ 2-1в» НА ОСНОВЕ ПРИНЦИПА МИНИМУМА УПРАВЛЕНИЯ
© В. Ф. Петрищев
АО «Ракетно-космический центр «Прогресс», г. Самара
Статья поступила в редакцию 23.10.2015
Разработанный автором принцип минимума управления, обеспечивающий минимизацию расхода ресурсов на управление, рассмотрен применительно к задаче синтеза дискретной оптимальной замкнутой системы управления модулем вектора тяги двигателей ракеты-носителя «Союз 2-1в» как материальной точки. Задача решается при следующих упрощениях: предположении о сферичности фигуры Земли, использовании аналитических зависимостей изменения плотности атмосферы и скорости звука с высотой, изменения коэффициента лобового сопротивления по траектории. Для моделирования движения центра масс использованы уравнения в конечных разностях. Предполагается возможность форсирования и глубокого дросселирования тяги двигателей обеих ступеней при условии сохранения их удельного импульса. Моделируется старт с космодрома «Плесецк». Движение рассматривается плоским в полярной системе координат. Приведены результаты моделирования движения ракеты-носителя. Показано, что использование предложенного принципа минимума управления повышает высоту орбиты выведения того же полезного груза. Ключевые слова: качество, ковариация, орбита, система, состояние, ступень, тяга, траектория, функция Ляпунова, управление.
Существующие методы синтеза оптимальных систем управления базируются на фундаментальных результатах применения вариационного исчисления, принципа максимума Понтрягина и динамического программирования Беллмана. Во всех этих методах в результате синтеза определяется некоторая вспомогательная система, сопряжённая заданной, которая обеспечивает оптимальное программное движение заданной системы. Для определения параметров сопряжённой системы необходимо, как правило, решать двухточечную краевую задачу.
С целью исключения необходимости решения двухточечной краевой задачи автором в [1] разработан принцип минимума управления. Согласно этому принципу в качестве критерия оптимальности выбран минимум следа ковариационной матрицы управления на каждом шаге. При этом структура вспомогательной системы и структура оптимального регулятора задаются из эвристических соображений.
Принцип минимума управления оказался весьма конструктивным, поскольку заданный критерий качества системы является в то же время положительно определённой функцией Ляпунова, доставляющей системе условие асимптотической устойчивости в целом.
В работе [2] рассмотрено применение принципа минимума управления к задаче подъёма геофизической ракеты на максимальную высоту на примере геофизической ракеты В-2А, разрабо-
Петрищев Владимир Фёдорович, доктор технических наук, ведущий научный сотрудник. E-mail: mail@samspace.ru
танной в 1958 году под руководством С. П. Королёва. Как известно [3], ракета была предназначена для исследования высотных слоёв атмосферы и двигалась по вертикальной траектории.
В работе [2] показано, что применение принципа минимума управления величиной тяги жидкостного ракетного двигателя позволяет увеличить максимально достижимую высоту подъёма ракеты В-2А (212 км) на 10 %.
1. ИСХОДНЫЕ ПОЛОЖЕНИЯ, ИСПОЛЬЗУЕМЫЕ ПРИ ФОРМИРОВАНИИ УПРАВЛЕНИЯ ТЯГОЙ ДВИГАТЕЛЕЙ
Рассматривается плоское движение ракеты-носителя «Союз 2-1в» как материальной точки в неподвижной полярной системе координат, связанной с точкой старта на момент старта и с началом в центре масс Земли, под действием вектора тяги двигателей первой и второй ступеней в центральном поле тяготения с ускорением:
g =
(Ro + h )2
где р =3,986-105 км3/с2 - гравитационный потенциал Земли, =6371 км - средний радиус Земли, Ь - высота полёта над поверхностью Земли.
В качестве возмущающей силы рассматривается лишь одна сила лобового аэродинамического сопротивления^ которая зависит от параметров движения и вычисляется по формуле:
2
К = р(ь)-^-8тех, где р(Ь) - плотность атмосферы, зависящая от
высоты Ь, V - модуль вектора скорости РН, Бт =5,865 м2 - площадь миделевого сечения РН, Сх - безразмерный коэффициент лобового аэродинамического сопротивления.
Модель изменения плотности атмосферы с высотой принята в виде экспоненты
р(Ь) = Ро • е-а\ где р0 =1,225 кг/ м3 , а=1,2-10-4 1/м .
Принятый график изменения коэффициента Сх в зависимости от числа М (Маха) представлен на рис. 1. До числа М=1 Сх линейно возрастает, дальнейшее его изменение описывается законом
1.5-М
Сх=15 -
1.5 + е
1.5-М
Дифференциальные уравнения второго порядка управляемого движения РН как материальной точки в соответствии с [4] записываются
в виде:
р = -р(И)• Бш• Сх • V-р/2М -- 2р • г/г + ит /(г• М); & = -р(Ь)• Бш^Сх • V• Г/2М +
(1)
2
+ р • г - g + иг /М ,
где иг, и т - проекции вектора тяги двигателя на направление радиус-вектора центра масс РН и ортогональное ему в направлении полёта соответственно, М - текущее значение массы РН.
Вводя четырёхмерный вектор состояния системы управления тягой двигателей РН для плоского случая движения
х = [Х1, х2
Х4]Т =
(2)
= [р, р, г, Г]т,
приведём систему (1) двух дифференциальных уравнений второго порядка к системе четырёх дифференциальных уравнений первого порядка и запишем их в матричной форме:
0
1
0 - Е
0
0
0
0 0 1
0
0 - g/x3 - р2
х +
0 0 1/(гМ) 0 0 0
0 1/М
(3)
и т
и
Рис. 1. График зависимости коэффициента Сх от числа М
Число М есть отношение текущей скорости движения РН к скорости звука Узв(Ь) на данной высоте
М = v/v зв (И).
Зависимость величины скорости звука от высоты описана аналитически алгебраическим полиномом четвёртой степени
V (Ь)=а, .Ь4+Ь, .Ь3+с, .Ь2+^ .Ь+е ,
зв^ ' 1 1 1 1 1 >
где коэффициенты полинома получены путём статистической обработки параметров, а именно скорости звука, в стандартной атмосфере Земли по ГОСТ 4401-81:
а1=0,122521798.10-16 ;
Ь1=-0,31247614.10-11;
^=0,234850942.10-6;
^=-0,0058183409 ;
е^340,4959676.
2. УРАВНЕНИЯ ДВИЖЕНИЯ РН
Текущие координаты центра масс обозначим через р и г , а их производные по времени - р и г .
где
^ = р(И)• Бш"Сх • V/2М + 2г/г;
Г2 = р(Ь)• Бш^Сх • V/2М-р2 • г/г.
Интегрируя это матричное уравнение по известной формуле Коши на малом временном интервале А1 = 0,1с от 1 до 1+1, получим уравнение движения РН в конечных разностях в матричной форме:
4+1
1 А1
0
0 1 - ^ 0
0 0
1
0
А1
0 0 -gАt/x31 -Г2А1
(4)
Аt2 / 2(гМ)
0
А^(гМ) - F1Аt2/2(rM) 0
Аt2 / 2М 2
0 Аt/M-Е2Аt2 /2М
В сжатой форме это уравнение записывается в виде:
х
3
0
0
+
0
х +
X
+
г
X1+1 = А; • X; + В; • и; . (5)
К этому уравнению следует добавить соотношение для изменения массы РН на одном шаге М •
расчетов
М;+1 = М; -т •М ;
где т - секундный расход массы топлива, определяемый отношением текущей тяги двигателя
И; =(2 + Иг2 )1
/2
Ь = (
к удельному импульсу.
3. РАСЧЁТ МОДЕЛЬНОЙ ТРАЕКТОРИИ РН
Для проведения расчета модельной траектории РН необходимо задать программу изменения во времени угла наклона касательной к траектории по отношению к местной вертикали. Штатная программа учитывает большое количество различного рода ограничений, которые в настоящей работе не представляется возможным учесть. (Такие, например, как допустимые районы падения отработавшей первой ступени и головного обтекателя). Поэтому положим, что закон изменения составляющей Иг, направленной по радиус-вектору центра масс РН, описывается с использованием функции гиперболического тангенса от текущего шага расчета 1 (и, следовательно, времени):
И г =( - ^ )• и;
и т=(и2 - иг2 )1/2
(6)
% = (((- 2,88 +1 / 63844) + 1)/ 2;
где и - текущее значение модуля вектора тяги двигателя. Тем самым функция задает программу изменения угла наклона касательной к траектории. При этом учитывается, что для двигателя первой ступени тяга возрастает с высотой за счет уменьшения атмосферного давления. График изменения функции 1- Гц по времени представлен на рис. 2.
Функция Г ; получена подбором входящих в нее параметров, исходя из условия приближенного равенства нулю радиальной скорости движения второй ступени РН в момент окончания работы ее двигателя: г = 11 = 0 (РН выводится в точку перигея). Функция гиперболического тангенса весьма успешно обеспечивает выполнение граничных условий, как в точке старта, так и в точке выведения.
С использованием приведенных выше данных и функции Г ; проведен расчет траектории РН « Союз 2-1в», стартующей с космодрома «Плесецк». В результате расчета РН выводит головной блок массой 2100 кг на эллиптическую орбиту с наклонением 82,4о и с высотами перигея 270 км и апогея 595 км. Графики изменения: орбитальной скорости V т = р • г, высоты полета Ь = г - Я 0 ,
Рис. 2. График изменения функции 1-^ по времени
радиальной скорости 1 = 11 на витке выведения, а также модуля вектора тяги и и его составляющих и т и Иг, скоростного напора Ч = р(Ь)• V2 / 2 и изменения массы РН М на активном участке траектории приведены на рис. 3 - рис. 10. При этом отсчет высоты ведется относительно принятого среднего радиуса Земли.
При первом запуске 28.12.2013 г. РН «Союз 2-1в» обеспечила вывод в тех же условиях той же полезной нагрузки на орбиту с высотами перигея 271,4 км и апогея 621,9 км. При этом отсчет высот ведется относительно общеземного эллипсоида. Некоторое отличие в параметрах модельной и фактической орбит выведения вызвано принятыми допущениями.
УТ, м/с | аооо
о 1000 2000 зош мло яюо " с
Рис. 3. График изменения трансверсальной скорости на витке выведения
Рис. 4. График изменения высоты на витке выведения
Рис. 7. График изменения трансверсальной тяги на активном участке
Рис. 5. График изменения радиальной скорости на витке выведения
Рис. 8. График изменения радиальной тяги на активном участке
Рис. 6. График изменения модуля вектора тяги Рис. 9. График изменения скоростного напора на активном участке на активном участке
М, кг
160000 140000 120000' 100000 80000 60000 ¡,0000 20000
4 ■ 10, с
11 1000 2000 3000 4000 5000 6000
Рис. 10. График изменения массы РН на активном участке
4. ЗАКОН ОПТИМАЛЬНОГО УПРАВЛЕНИЯ ТЯГОЙ ДВИГАТЕЛЕЙ РН
В соответствии с принципом минимума управления [1] для рассматриваемой системы (5) формируется вспомогательная система
% 1+1 = С•% 1 + D1 • Иь (7)
где % - вектор состояния системы той же размерности, что и вектор xi, и с теми же компонентами. Матрицы 01 и Б1 имеют размерности матриц Ai и Bi соответственно. Вспомогательная система управляется тем же вектором И, что и система (5). Кроме того, начальное состояние вспомогательной системы совпадает с начальным состоянием заданной системы:
% 0 = х0 .
От матрицы С1 требуется, чтобы невозмущённое движение вспомогательной системы (в случае И = 0) было асимптотически устойчиво в целом. Что касается матрицы Б1, то полагаем, что она удовлетворяет равенству: Б1 = -Б1 .
Закон оптимального управления с обратной связью задаётся в линейной форме:
и 1+1 = и + Р1+1 -(х изм ,1+1 -%1+1), (8)
то есть, управление на текущем шаге определяется в виде алгебраической суммы управления на предыдущем шаге и взвешенной разности векторов состояния заданной и вспомогательной систем на текущем шаге. Здесь Рм - весовая матрица, оптимальным образом взвешивающая разность между векторами состояний заданной и вспомогательной систем. Вектор состояния заданной системы х изм 1+1 определяется в результате измерений соответствующими датчиками на каждом шаге расчётов:
х и
,1+1 = х1+1 + 81+1
(9)
ностей измерений компонент вектора Х1+1 типа дискретного белого шума с заданной ковариационной матрицей погрешностей измерений, диагональными элементами которой являются дисперсии погрешностей измерений.
При некоторых дополнительных предположениях о ковариационных свойствах векторов, входящих в соотношения (7) - (9), о которых будет сказано ниже, оптимальное выражение для весовой матрицы Рм, найденное из условия минимума следа ковариационной матрицы управления
Ки,1+1 = м( • ит+1), ао)
где М - знак математического ожидания, т - знак транспонирования вектора, И;+1 - вектор управления, имеет вид:
Р+1 =-Ки,1 (Б1 - В1 )тх
-1 (11)
+ С1К %,1СТ +
где х1+1 - истинное значение вектора состояния, +1 - вектор случайных аддитивных погреш-
А1Кх,1АТ
х
.+ ( - В1 )Ки,1 (Б1 - )т + К6
Значения входящих в это выражение ковариационных матриц состояния рассматриваемой и вспомогательной систем и ковариационной матрицы погрешностей измерений
Кх,1 = м( • хТ);К%,1 = м(%1 -%т);Кв = м( -8т)
считаются известными. Они могут быть заданы постоянными. Начальное значение ковариационной матрицы управления (10) также известно.
Показатель степени (-1) в (11) означает операцию обращения суммы матриц.
Уточнённое на текущем 1+1 шаге расчёта значение ковариационной матрицы управления с учётом полученного значения оптимальной весовой матрицы (11) на этом шаге имеет вид:
Ки,1+1 = [+1(1 - В1) + Е]-Ки,1, (12)
где Е - единичная матрица соответствующей размерности.
Начальное значение вектора управления и0 в (8) определяется из условий решаемой задачи. В настоящей задаче вектор управления должен быть достаточным для отрыва РН от стартового стола.
5. РАСЧЁТ ОПТИМАЛЬНОЙ ТРАЕКТОРИИ РН
В данной задаче оптимизируется величина модуля тяги двигателя РН. Направление вектора тяги при этом принимается равным принятому при расчёте модельной траектории. Для решения этой задачи вектор управления в (4) с учётом (6) и функции гиперболического тангенса : =((- 3,0 +1/400)+1)/2,
преобразуем к виду:
(13)
где Ц - модуль вектора тяги двигателя на шаге 1. Подставляя это выражение в правую часть (4), получим
"И х" (2 • ^2,1 - ^2,1) • и
.Иг _ 1 . 1 - ^,1 _
xi+1 -
1 At
0
0 1 - F1At 0
00
0
1
•X; +
/2
0 0 -gAt/x31 -F2At (At2 / 2(rM)) -(2 • f2,i - f22; )1 (At /(rM) - F1 At2 /2(rM)) • (2 • f2,; - f22 (At2 / 2M) • (1 - f2,;) (At /M - F2At2 /2M) • (1 - f2,;)
• U;
Здесь матрица В1 превращена в четырехмерный вектор. Управляющим параметром системы служит модуль вектора тяги. Его оптимизация на каждом шаге управления - цель настоящей работы.
Для проведения расчета оптимальной траектории РН в соответствии с материалами предыдущего раздела должны быть заданы значения ковариационных матриц Кх ,К ^ , К Е и матрицы С .
Значения элементов матрицы С1 вспомогательной системы (7) выберем с использованием функций гиперболического тангенса и 141, аналогичных функции (6), примененной при построении моделной траектории:
^ = ((- 2,0 +1/500) + 1)/2;
' / V ч ч (14)
Г4; = (((- 2,2 +1 / 400) +1)/ 2 ;
C; -
At
-f3,i /101,5 1 - f3,;/102
At
0
0 -f4i /103'5 1 -f4i /102
. (15)
Значения начальных векторов состояния систем в исходной системе координат, учитывающие вращение связанной с Землей точки старта, принимаются (с теми же размерностями) следующими:
х0 = £,0 = [0; 9,528-10-6 ; 6371000; 0]Т .
Перенося начало системы координат в точку выведения, исходные начальные векторы состояния основной и вспомогательной систем в новой системе координат, записанные в тех же обозначениях, но с волнистой чертой сверху, будут иметь значения (с упомянутыми размерностями), равные разности начальных и конечных значений:
~0 = ~0 = [-0,194 ; -0,0012144; -280000, 0]Т .
Из простых геометрических построений можно получить выражения для расчета вектора состояния основной системы, записанные в новой системе координат:
- rk - 2 • r • rk • cos „2 2
(pk -P))/2 ;
p - arccos((r 2 + rk2 - r2 )/(2 • r • rk )/;
a - arccos((r2 + r2 - r] //(2 • r • r)/; A — n-a;
Асимптотическая устойчивость невозмущенного движения вспомогательной системы с использованием данной матрицы С. будет подтверждена ниже численным путем в процессе проведения расчета.
Следуя концепции построения замкнутой системы с обратной связью, начало координат принятой полярной системы необходимо поместить в конечную точку выведения РН. В этом случае векторы текущего состояния систем должны определяться относительно нового начала координат.
Желая вывести РН на большую высоту, чем она выводится в штатном варианте, поставим задачу выведения РН на орбиту в точку перигея с высотой перигея 280 км и высотой апогея 650 км. В исходной полярной системе координат конечный вектор состояния в этом случае должен быть равен
хк = [х1,к; х2,к; хз,к; х4,к] =
= [0,194 ; 0,0012239 ; 6651000 ; 0] Т .
Здесь размерности элементов: рад, рад/с, м и м/с соответственно.
~ - -(r • cos A + p • r • sin A); p--(p • r • cos A-r • sin A)/p + p ] .
Здесь Д - меньший угол между векторами г и
r в текущей точке положения РН. У компонентов вектора состояния в исходной системе координат для простоты опущен индекс «измеренный».
Аналогичные соотношения имеют место и для вспомогательной системы.
Таким образом, с использованием введённых обозначений закон управления модулем вектора тяги двигателей (8) можно записать в окончательном виде:
Ui+1 - Ui + Pi+1 (ХИзМ,i+1 - pi+1/.
Поскольку размерность вектора управления равна скалярной единице, единичная матрица Е в равенстве (12) вырождается в скалярную единицу.
Выбраны следующие значения ковариационных матриц:
Kx -
X - K £ -
10 -0 0 0
10 -12 1/с 2 0 0
0 10-2 м 2/ с 2
Ки0 = 100 Н2 .
Ковариационная матрица погрешностей измерений принята нулевой.
Результаты расчетов оптимальной траектории РН приведены на рис. 11 - рис. 20.
0
0
+
1
0
0
0
0
0
0
1
16
0
0
0
0
0
10-4 м 2
0
УТХ, м/с 8000
0 ' 500 1000 1500 2000 2500 3000 " 1 С
Рис. 11. График изменения трансверсальной скорости на активном участке
иТ, Н'
1.2е+07
1000 1500 2000 2500 3000
Рис. 15. График изменения трансверсальной тяги на активном участке
Нх, км
1000 1500 2000 2500
г - ю, с
Рис.12. График изменения высоты на активном участке
и-г, н
5а+0б
500 1000 1500 2000 2500 3000
£ ■ 10, С
Рис. 16. График изменения радиальной тяги на активном участке
Кх, м/с 1600 1А-00 1200 1000 800 600 Л-ОО 200
1000 1500 2000 2500
£10, с
Рис. 13. График изменения радиальной скорости на активном участке
и, н
1.2е+07
500 1000 1500 2000 2500 3000
г ■ ю, с
Рис. 14. График изменения модуля вектора тяги на активном участке
Я, Н/м2 100000
0 500 1000 1500 2000 2500 ЗООО'
£ • 10, с
Рис. 17. График изменения скоростного напора на активном участке
М, кг 160000
0 500 1000 1500 2000 2500 3000
£ ■ 10, с
Рис. 18. График изменения массы РН на активном участке
Рис. 19. График изменения ковариационной функции управления на активном участке
Как и в случае модельной орбиты, выведение осуществляется в точку перигея (рис. 13).
По результатам расчетов оптимальная орбита выведения имеет следующие параметры: высоту перигея 200 км (рис. 12) и высоту апогея 3980 км (это следует из рис. 11). Для получения этих параметров двигатели РН включаются дважды. Начальная тяга двигателей увеличена практически в два с половиной раза по сравнению со штатной величиной тяги, поэтому даже четырехкратное увеличение максимальной величины скоростного напора (рис. 17) не препятствует увеличению высоты апогея.
Функция Ляпунова (рис. 20)
и; = ( -4;)Т -(С; )
на траектории движения вспомогательной системы является монотонно убывающей функцией. Отсюда следует, что необходимое условие теоремы Ляпунова - отрицательная определенность приращения этой функции на каждом шаге управления - выполнено. Следовательно, матрица С обеспечивает асимптотическую устойчивость в целом невозмущенного движения вспомогательной системы.
ЗАКЛЮЧЕНИЕ
Приближенно решена задача оптимального выведения заданной полезной нагрузки ракетой-носителем «Союз 2-1в» на орбиту с космодрома «Плесецк». Рассмотрено плоское движение РН как материальной точки. Программа изменения угла наклона касательной к траектории задавалась подобранной функцией гиперболического тангенса угла в функции времени. Правильность выбранного математического инструментария подтверждена близостью параметров модельной и фактической орбит выведения:
- модельная орбита 270х595 км (над сферической Землей радиуса 6371 км);
- фактическая орбита после запуска 28.12.2013г.
Рис. 20. График изменения функции Ляпунова на активном участке
271,4x621,9 км (над общеземным эллипсоидом).
На основе принципа минимума управления осуществлён синтез оптимальной замкнутой системы управления модулем тяги двигателей РН « Союз 2-1в». В замкнутой системе управление осуществлялось относительно точки выведения. Параметры оптимальной орбиты выведения 200x3980 км.
В результате проведенных расчётов подтверждена эффективность принципа минимума управления.
Вместе с тем, настоящая работа не может рассматриваться в качестве предложения по повышению эффективности ракеты-носителя «Союз 2-1в». Она демонстрирует лишь эффективность разработанного автором принципа минимума управления. В настоящее время ещё не созданы двигатели, допускающие форсирование и дросселирование тяги в широких пределах при условии сохранения величины удельного импульса. Кроме того, форсирование тяги двигателей неразрывно связано с необходимостью увеличения прочности конструкции РН. А это, в свою очередь, отразится на её массовых характеристиках, так как при сохранении общей массы РН неизбежно приведёт к увеличению массы сухой РН и уменьшению массы рабочего тела.
СПИСОК ЛИТЕРАТУРЫ
1. Петрищев В.Ф. Принцип минимума управления в задаче синтеза линейных дискретных систем. // Материалы научно-технической конференции «Актуальные проблемы ракетно-космической техники» (III Козловские чтения). 2013. С. 513 - 521.
2. Петрищев В.Ф. Принцип минимума управления в задаче синтеза дискретных систем / 7-я Российская мультиконференция по проблемам управления // Материалы конференции «Управление в морских и аэрокосмических системах» (УМАС-2014). 7-9 сентября 2014 г. Санкт-Петербург. С 64-73.
3. Кирилин А.Н., Ахметов Р.Н., Тюлевин С.В. и др. Самарские ступени семёрки // Самара, Издательский
дом «Агни».2011. 253 с. 4. Аппазов Р. Ф., Сытин О. Г. Методы проектирования
траекторий носителей и спутников Земли // М.: «Наука», 1987. 440 с.
OPTIMIZATION OF SOYUZ 2-1V LAUNCH-VEHICLE ENGINE THRUST CONTROL BASED ON THE PRINCIPLE OF MINIMUM CONTROL
© 2015 V.F. Petrishchev
JSC Space-Rocket Centre PROGRESS, Samara
The author developed a principle of minimum control, which is aimed to minimize energy consumption necessary for control, and applied it for synthesis of a discrete optimal closed-loop system controlling engine thrust vector magnitude of a Soyuz 2-1v launch-vehicle as a particle. The problem is solved with the following simplifications: assumption of the Earth's figure sphericity, use of analytic dependences of atmospheric density and sound velocity variation with altitude, and drag coefficient variation with trajectory. To model centre-of-mass motion, the author used finite-difference equations. It is assumed that it is possible to force and throttle the engine thrust of both stages given conservation of their specific thrust. A launch from the Plesetsk cosmodrome is simulated. The motion is considered to be two-dimensional in polar coordinate system. The article includes the results of the launch-vehicle motion simulation. It also shows that application of this principle of minimum control results in increase in injection orbit altitude for the same payload.
Key words: quality, covariance, orbit, system, state, stage, thrust, trajectory, Lyapunov function, control.
Vladimir Petrishchev, Doctor of Technics, Leading Research Fellow. E-mail: mail@samspace.ru
