CC BY
810
Наука к Образование
МГТУ им. Н.Э. Баумана
Сетевое научное издание
Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. № 01. С. 53-67.
ISSN 1994-0408
Б01: 10.7463/0115.0755072
Представлена в редакцию: Исправлена:
© МГТУ им. Н.Э. Баумана
УДК 629.764
Оптимизация траектории выведения ракеты-носителя с малым космическим аппаратом на солнечно-синхронную орбиту на основе
псевдоспектрального метода
1 *
Ван Л.1,
10.01.2015 27.01.2015
УапаИг&е_Ьтаи'ЗтаЦ.ш 1МГТУ им. Н.Э. Баумана, Москва, Россия
Рассматривается задача оптимизации траектории выведения ракеты-носителя воздушного базирования с малым космическим аппаратом на солнечно-синхронную орбиту. На основе программы трехэтапного полета с пассивным промежуточным участком, используется математическая модель движения ракеты-носителя с учетом внешних возмущений, вызванных несферичностью Земли, атмосферным сопротивлением и ветром. Для решения задачи проводится процедура дискретизации, базирующаяся на псевдоспектральном методе, позволяющая преобразовать исходную задачу в задачу нелинейного программирования с динамическими ограничениями и критерием максимума конечной массы на целевой орбите. Применение предлагаемой процедуры решения иллюстрируется примером расчета оптимального управления и соответствующей траектории выведения для двухступенчатой жидкостной ракеты-носителя, выводящей малый космический аппарат на солнечно-синхронную орбиту высотой 512 км. В численном результате показывается эффективность предложенного метода, а также анализируется параметры полученной трехэтапной траектории с промежуточным пассивным участком полета.
Ключевые слова: ракета-носитель воздушного базирования, малый космический аппарат, солнечно-синхронная орбита, псевдоспектральный метод, оптимизация траектории выведения
Введение
Исследования экономичных транспортных космических систем для обеспечения запуска малых космических аппаратов (МКА) на орбиты спутников Земли проводятся во многих странах. Одним из вариантов малозатратных технических решений является использование авиационно-космических комплексов (АКК), включающих в себя самолет-носитель и ракету-носитель воздушного базирования. АКК отличается от традиционной системы наземного старта важнейшим преимуществом - возможностью выбора географического места старта, что позволяет исключить длительное время ожидания достижения требуемых условий пуска [1]. Первым из числа реализованных проектов
ракет-носителей воздушного базирования является ракета-носитель (РН) «Пегас», ориентированный на выведение МКА массой до 450 кг [2].
Применение ракет-носителей воздушного базирования для выведения МКА с целью дистанционного зондирования Земли требует высокой точности выхода на заданные солнечно-синхронные орбиты (ССО), поэтому становится задача построения программной оптимальной траектории выведения и реализующего ее оптимального управления. Проблематика построения программного движения РН воздушного базирования была рассмотрена в работах О.М. Алифанова и др. [3], а также И.А. Пышного и В.Е. Чепиги [1]. В указанных работах задача рассматривалась в рамках постановки плоского движения на основе решения двухточечной краевой задачи. В диссертации Д. В. Мазгалина [4] рассмотрены возможные способы решения задачи определения оптимального управления для построения программной траектории выведения РН типа «Союз-2» на эллиптические орбиты. Следует отметить, что в данной диссертации предложены разные законы управления для каждой ступени ракеты-носителя. Задача многоэтапного оптимального управления, в которой осуществляется дискретизация интервала времени полета на основе псевдоспектрального метода и применяется концепция многоэтапного оптимального управления, рассмотрена в работе [5]. Разработанный метод использован для оптимизации траекторий выведения ракеты-носителя «Delta III» на геостационарные орбиты. В работе [5] при оптимизации траекторий выведения ракет-носителей в системах наземного старта использовался ряд упрощающих предположений о модели полета. В частности, в указанных работах не учитывались возмущающие воздействия, обусловленные влиянием атмосферы и несферичностью поля тяготения Земли.
В настоящей работе для ракет-носителей воздушного базирования решена задача оптимального управления полетом (по критерию максимума полезной массы) с учетом атмосферного и гравитационного возмущений. Существенным отличием алгоритма управления полетом является наличие оптимизируемого пассивного участка перелета, а также солнечно-синхронная целевая орбита ракеты-носителя.
1. Постановка задачи непрерывного оптимального управления
1.1. Уравнение движения центра масс ракеты-носителя в инерциальной стартовой
системе координат
При описании движения центра масс ракеты-носителя будем использовать инерциальную стартовую систему координат Oxyz. Начало O в центре общего земного эллипсоида (ОЗЭ). Ось Oy от начала O вертикально вверх направлена на стартовое место в момент времени запуска, а плоскость Oxz параллельна плоскости стартового горизонта, причем Ox указывает направление пуска.
В задаче выведения РН воздушного базирования на заданную ССО используем математическую модель движения цента масс РН на всей траектории полета. Соответствующая система дифференциальных уравнений, описывающая движение РН в
инерциальной стартовой системе координат от момента старта to до момента tf, может
быть записана в следующем виде:
dr dv , N ттг „7 dm .
- = v = g(r)+wR■u+w,, -= -m , (1.1)
где, r и v - радиус-вектор и вектор скорости центра масс РН в СК Oxyz, m - масса РН,
m- кусочно-постоянная функция, характеризующая массовый секундный расход топлива, g - вектор гравитационного ускорения с учетом возмущения, вызванного несферичностью Земли и вычисляемого в соответствии с EGM2008 [6].
u = (cos &cos у sin &cos y, -sin у) - единичный вектор направления вектора тяги в СК Oxyz, О - угол тангажа, у - угол рыскания, Wr - кажущееся ускорение, обусловленное реактивной тягой ДУ РН, W, - вектор сопротивления, создаваемого аэродинамическими силами. При вычислении W, принята атмосферная модель в соответствии с ГОСТ 440181 [7] и данные о величине и направлении ветра в соответствии с ОСТ [8].
Величина кажущегося ускорения, обусловленного реактивной тягой ДУ РН, вычисляется по соотношению:
Wr (í) = (P - Sapa (H))/m(t), (1.2)
где, P - реактивная тяга ДУ в вакууме, Sa - суммарная площадь сечения сопел ДУ, pa (H) - атмосферное давление, H - геодезическая высота полета в текущий момент времени t.
Сопротивление, создаваемое аэродинамическими силами, имеет вид:
W, = -1 ■ Р ■ Cd ■ Smid ■ \Vrel \' Vrel / m(t) , (13)
где, р - атмосферная плотность, Cd - коэффициент аэродинамического сопротивления, Smid - площадь миделя сечения, Vrel - скорость РН относительно атмосферы.
Для системы дифференциальных уравнений (1.1) задаются начальные условия:
r(t0) = r0, v(to) = vq, m(tQ) = mo, (1.4)
Вектор направления реактивной тяги u рассматривается как управляющее воздействие. Поскольку вектор u является единичным, для него справедливо условие:
\u\ = 1. (1.5)
1.2 Программа полета ракеты-носителя
Оптимальная траектория выведения ракеты-носителя включает в себя три участка: начальный активный участок полета, промежуточный пассивный участок полета и конечный активный участок полета [1, 9]. Упрощенно оптимальная траектория выведения на низкую орбиту может быть представлена следующим образом: запуск ракеты-носителя на промежуточную орбиту, апогей которой находится на высоте целевой низкой орбиты;
пассивное движение по промежуточной орбите; выход на целевую орбиту в апогее промежуточной орбиты.
В схеме запуска МКА с использованием ракеты-носителя авиационного базирования «Пегас» реализована аналогичная траектория выведения [2].
В соответствии с участками полета, перепишем систему дифференциальных уравнений (1.1) для каждого из участков:
дг (Г) ш* (Г) . . тт. ш ш
= - , (16)
аЬ
где, - = 1,...,3 - номер участка полета.
Между последовательными участками полета должны удовлетворяться непрерывные условия, т.е. радиус-вектор, вектор скорости, масса ракеты изменяются по следующим соотношениям:
Г+1 ) = г,-), V,-+1 ) = V, ),
т,+1(1,) = т, (1,) -Ащ , (1.7)
где, знаки «,» и «г +1» соответствуют участкам полета до и после момента времени , Ат, - разность массы РН между последовательными участками полета. Если между последовательными участками полета происходит отделение ступени от РН, величина Ат, равна конструктивной массе отделяемой ступени; в противном случае Аmi равна нулю.
1.3 Целевая солнечно-синхронная орбита
Солнечно-синхронная орбита, на которую выводится РН, задается следующими параметрами (параметры оскулирующей орбиты): наклонение орбиты / , эксцентриситет е, большая полуось а .
Постановка задачи: Для системы (1.6) с заданными начальными условиями (1.4) в момент старта, найти управления и , максимизирующие значение функционала
3 = т(^). (1.8)
При этом должны быть выполнены ограничения (1.5), (1.7) и
фЩ),v(tf))-Т <5,, феЩ),v(tf))-е <5е,
Фа (г(tf), v(tf)) - а\<5а, (1.9)
где 5,, 5е, 5а - допустимые отклонения от номинальных элементов орбиты, ф,, фе, фа -соотношения между элементами орбиты и векторами положения и скорости.
2. Представление сформулированной задачи в форме задачи нелинейного программирования
При решении рассматриваемой оптимизационной задачи будем применять псевдоспектральный метод Ье§епёге-Оаш8-Яаёаи (ЬОЯ) [10], который обладает быстрой сходимостью при нахождении оптимального управления. Применяя в рамках этого подхода процедуру дискретизации непрерывной задачи посредством аппроксимации переменных состояния и управляющих переменных интерполяционными полиномами, сводят ее к задаче нелинейного программирования, которая далее решается известными численными методами нелинейного программирования.
2.1 Узловые точки, используемые в псевдоспектральном методе ЬОЯ
Для использования метода ЬОЯ [10] необходимо перейти к дискретной модели полета и определить точки коллокации ЬОЯ ту е(—1,1], ] = 1,К ,N, являющие корнями
разности полиномов Лежандра степени N и N—1: т)- РN—1( т).
На основе псевдоспектрального метода ЬОЯ, с использованием узловых точек, включающих в себя точки коллокации ЬОЯ и дополнительную точку то =—1, можно
представить приближенную оценку переменных состояния интерполяционными функциями Лагранжа:
N
X(т) - Е XуЬу (т), (2.1)
]=0
где, Ьу - базисные полиномы Лагранжа:
N т — т,
Ьу(т) = п -К] = 1,..., N. (2.2)
к=0,кФ] т] —тк
2.2 Преобразование интервала времени на диапазон узловых точек
Поскольку узловые точки ту е [—1,1], ] = 0, К, N, поэтому нам приходится
преобразовать реальный интервал времени г е —1,^](/ = 1,...,3) каждого участка полета на т е [—1,1] путем следующего соотношения:
1 = (Ц —!, — 1)т + (Ц +!, —1) (23)
Подставив выражение (2.3) в уравнения (1.6), получим приближенное дискретное представление системы дифференциальных уравнений на у -ой дискретной точке каждого участка полета:
^(т]) ь — ь—1
й т 2 "
(хj)_ ц - Ц-1
Ш х
{g(гг-) + ЖК • и, + }
(2.4)
в которых знаки «1», ..., « N » соответствуют LGR-точкам х^ е (-1,1], j = 1,..., N.
2.3 Матрица приближенного дифференцирования
Дифференцируя выражение (2.1) относительно аргумента х, получим приближенные производные переменных состояния на LGR-точках:
Обозначаем через Б = {Оу* } матрицу приближенного дифференцирования размерностью N х ^ +1), элементы которой равны производной функции Ь* (хj) относительно аргумента х j, т.е.
)
Используя матрицу Б приближенного дифференцирования, перепишем (2.4) в следующем виде:
• г,(х*) = ^^• V,(х;), к=0 2
Е % • V, (х*) = ^^ • (хj)) + (хj) • и (хj) + ^А/ (хj)}, к=0 2
Следует отметить, что приближенные уравнения (2.6) описаны на N LGR-точках. В этих уравнениях применяются переменные управления тоже только на N LGR-точках, а дискретные переменные состояния используются на N +1 дискретных точках.
2.4 Непрерывное условие между последовательными участками
Поскольку траектория выведения состоит из нескольких участков, между ними может происходить прерывистое изменение переменных состояния. Если мы разделим всю траекторию на участки, на которых непрерывно изменяются переменные состояния, а между участками добавим ограничение для переменных состояния последовательных участков, то сформулируем задачу многоэтапного непрерывного оптимального
2
управления. Такая идея, именуемая концепцией многоэтапного оптимального управления, позволяет относительно просто решать сложную задачу оптимизации.
С учетом условий изменения переменных состояния (1.7), установим ограничение между последовательными участками:
г/+1(т0) = г/ (т N ), У/+1(т0) = У/(т N ),
Щ/+1(т0) = Щ (тN ) — Ат/ . (2 7)
2.5 Остальные ограничения
Преобразуем неравенства (1.9), содержащие абсолютную величину следующим образом:
—5/ <ф/ (г3(т N X У3(тЖ )) — ! <5/,
—5е <фе (Г3(т N ), У3(т N )) — е <5е, (2.8)
—5а < фа (г3(^ ), У3(т N )) — а < 5а.
Отметим, что в программе полета конечные моменты времени пассивного участка ^ и конечного активного участка полета ¿3 являются свободными параметрами, которые удовлетворяют условию:
¿3 > ¿2. (2.9)
В итоге проблему построения оптимальной многоэтапной траектории полета можно сформулировать как задачу нелинейного программирования: найти переменные состояния
г/(ту), V/(ту), щ(ту), ] = 0,...,N, управления и/(т*), / = 1,...,N, / = 1,...,3, а также
моменты времени ¿2 и ¿3 , минимизирующие значение функционала
7 = —т3(т N ) (2.10)
при динамических ограничениях (2.6) - (2.9).
В настоящее время существуют несколько алгоритмов нелинейного программирования, реализованных в виде программных пакетов 1РОРТ [11], БКОРТ [12] и др. В настоящей работе с помощью пакета 1РОРТ задача нелинейного программирования эффективно решается.
3. Численные примеры и анализ
Методика построения оптимальной траектории полета РН воздушного базирования, основанная на использовании псевдоспектрального метода, была протестирована на основе решения задачи выведения жидкостной двухступенчатой ракеты-носителя с МКА.
3.1 Характеристики двухступенчатой ракеты-носителя с ЖРД
В работе [13] приводятся результаты проектирования двухступенчатый РН с ЖРД, основные характеристики которой приведены в табл. 1. и рис. 1. Суммарная масса РН составит 18643,0 кг., секундный массовый расход ЖРД т = 105 кг/с, удельный импульс ЖРД =310 с.
Рис. 1. Геометрические характеристики ракеты-носителя Таблица 1. Массовые характеристики ракеты
«Сухая» масса Масса топлива, кг Длина, м. Диаметр сечения, м.
Первая ступень 1560,0 14040,0 9,25 1,8
Вторая ступень 254,0 2286,0 3,75 1,12
3.2. Результаты и анализ
Допустим, что ракета-носитель запускается из-под авиационного носителя над космодромом Цзюцюань. Географические координаты запуска РН составят
100°17'30"Е. Начальная скорость РН относительно поверхности Земли в направлении запуска равна 270 м/с, высота запуска - 12 км. Номинальные параметры солнечно-синхронной орбиты приведены в табл. 2.
Таблица 2. Параметры целевой орбиты
а , км е 1
Номинальные 6890,00 0,0 97,448°
В численном эксперименте количество LGR-точек на каждом участке полета равно 22, что позволяет, используя известный алгоритм нелинейного программирования, за короткое время (несколько секунд) вычислить оптимальную программу управления и соответствующую траекторию полета.
Зависимость геодезической высоты от времени приведена на рис. 2. Длительность выведения равна 501,07 с, конечная высота орбиты равна 511,13 км (из табл. 2 высота целевой орбиты равна 511,76 км).
Зависимости углов тангажа и рыскания от времени показаны на рис. 3. В результате получено, что интервал времени первого активного участка полета составит 0,0 с ~ 133,7 с, интервал времени конечного активного участка полета - 481,7 с ~ 501,07 с, длительность промежуточного пассивного участка полета равна 348,0 с. На рис. 4
приведена зависимость массы РН от времени, из которой заметим, что конечная масса на целевой орбите после отделения второй ступени составляет 806,8 - 254,0 = 552,8 кг.
Н, км 600 -500 400 300 200 100
0 -С-С-С-С-С-с
0 100 200 300 400 500 600 г, с
Рис. 2. Изменение высоты полета
80 60 40 20 0
0 20 40 60 80 100 120 140 с
а) первой активный участок полета
-24 -24.5 -25 -25.5 -26
-26.5 [— 480
_[_
485
_[_
490
_[_
495
_и_
500
505 г, с
-0 8 е-[-1-1-1-1
480 485 490 495 500 505 г, с
б) конечный активный участок полета
Рис. 3. Изменение углов тангажа и рыскания по времени сплошная линия - угол рыскания, штриховая линия
- угол тангажа
4
х 10
ш, кг. 2 г
Рис. 4. Изменение массы РН
Графики положения и скорости РН, вычисляемые в экваториальной инерциальной СК, от времени показаны на рис. 5 и 6, соответственно. В начальный момент времени (г = 133,7 с) промежуточного пассивного участка полета положение и скорость РН составят г = [4745099,46 -40158,46 4430400,55] м,у=[-1580,3 -711,8 4043,2] м/с. Из этих параметров орбиты следует, что геоцентрическое расстояние до точки перигея пассивной траекторий равно 1201,7 км. Интеграл энергии на суборбитальной траектории будет относительно мал. Поэтому полученная пассивная траектория может считаться более энергоэкономной. А в конечный момент времени (г = 501,07 с) конечного активного участка полета положение и скорость РН составят
г = [3932499,33 - 293176,99 5649154,67] м, у=[ - 6186,8 -1263,9 4241,4 ] м/с. Соответствующие элементы полученной орбиты а = 6890,023 км, е = 0,00010,
I = 97,438°.
3
2
1
0
100
200
300
400
500
600
г, с
0
Рис. 5. Изменение положения г = [X у z ] движения РН сплошная линия - X , штрихпунктирная линия - у, штриховая линия - 2 .
V, м / с 6000 Г-
Рис. 6. Изменение скорости V = ^Ух Уу У2 J движения РН сплошная линия - Ух , штрихпунктирная линия -
V , штриховая линия - vг.
Из графиков, приведенных на рис. 6, следует, что в силу малости V координата у
медленно меняется относительно других составляющих, и поэтому РН движется приблизительно в одной плоскости. На последнем активном участке полета скорость быстро изменяется в силу относительной малости оставшейся массы РН.
Заключение
Приведенные результаты позволяют сделать следующие выводы по работе. В работе предложена дискретная математическая модель движения ракеты-носителя воздушного базирования с учетом несферичности Земли, атмосферного сопротивления и ветра. Для увеличения конечной массы на целевой орбите, использована трехэтапная программа выведения РН, включающая промежуточный пассивный участок полета. На основе псевдоспектрального метода, разработан алгоритм численного решения задачи для оптимизации траектории трехэтапного выведения РН воздушного базирования с пассивным участком полета на ССО. Показано, что предложенный алгоритм позволяет эффективно найти оптимальную траекторию для выведения жидкостного двухступенчатого ракеты-носителя с МКА. Следует отметить, что в полученной траектории выведения промежуточный участок полета является суборбитальной траекторией с малым интегралом энергии, перигей которой находится под поверхностью Земли.
Благодарность
Выражаю глубокую благодарность доктору технических наук, профессору Лысенко Л.Н. за полезные советы, способствовавшие улучшению настоящей статьи.
Список литературы
1. Пышный И.А., Чепига В.Е. Запуск малых искусственных спутников Земли с использованием самолетов-носителей. М.: Машиностроение, 2005. 168 с.
2. Pegasus user's guide // Orbital: company website. Режим доступа: http://www. orbital. com (дата обращения 16.10.2014).
3. Алифанов О.М., Андреев А.Н., Гущин В.Н., Золотов А.А., Матвеев Ю.А., Перелыгин В.П., Хохулин В.С. Баллистические ракеты и ракеты-носители. М.: Дрофа, 2004. 512 с.
4. Мазгалин Д.В. Вопросы построения программной траектории выведения ракеты-носителя с космическим аппаратом: дис. ... канд. техн. наук. Екатеринбург, 2012. 136 с.
5. Benson D. A gauss pseudospectral transcription for optimal control: Ph.D. thesis. Cambridge, 2005. 243 p.
6. Vallado D.A. Fundamentals of astrodynamics and applications. Microcosm Press, New York: Springer, 2007. 1055 p.
7. ГОСТ 4401-81. Атмосфера стандартная. Параметры. Введ. 1982-07-01. М.: Изд-во стандартов, 1981. 180 c.
8. ОСТ 92-9704-95. Ракеты и ракеты-носители. Модели определения горизонтальной скорости ветра и термодинамических параметров атмосферы в диапазоне высот 0 -120 км в районе космодрома Плесецк. Введ. 95-01-01. М., 1995. 94 с.
9. Исаев В.К. Принцип максимума Л.С. Понтрягина и оптимальное программирование тяги ракет // Автоматика и телемеханика. 1961. Т. 22, № 8. С. 986-1001.
10. Garg D., Hager W.W., Rao A.V. Pseudospectral methods for solving infinite-horizon optimal control problems // Automatica. 2011. Vol. 47, no. 4. P. 829-837. DOI: 10.1016/j.automatica.2011.01.085
11. Biegler L.T., Zavala V.M. Large-scale nonlinear programming using IPOPT: An integrating framework for enterprise-wide dynamic optimization // Computers & Chemical Engineering. 2009. Vol. 33, no. 3. P. 575-582. DOI: 10.1016/j.compchemeng.2008.08.006
12. Gill P.E., Murray W., Saunders M.A. SNOPT: An SQP algorithm for large-scale constrained optimization // SIAM Review. 2005. Vol. 47, no. 1. P. 99-131. DOI: 10.1137/S0036144504446096
13. Ingrid M. D. Modeling dispersions in initial conditions for air-launched rockets and their effect on vehicle performance: Master's thesis. Cambridge, 2013. 159 p.
Science and Education of the Bauman MSTU, 2015, no. 01, pp. 53-67.
DOI: 10.7463/0115.0755072
Received: Revised:
10.01.2015 27.01.2015
Science^Education
of the Bauman MSTU
ISS N 1994-0408 © Bauman Moscow State Technical Unversity
Ascent Trajectory Optimization for Air-Launched Launch Vehicle with Small Sun-Synchronous Orbit Satellite Based on Pseudo-spectral Method
L. Wang1' "^^analit5e_bm&tuigmailju
1Bauman Moscow State Technical University, Moscow, Russia
Keywords: air-launched launch vehicle, small satellite, sun-synchronous orbit, pseudospectral
method, ascent trajectory optimization
Economical space transportation systems to launch small satellites into Earth's orbits are researched in many countries. Using aerospace systems, included aircraft and air-launched launch vehicle, is one of the low cost technical solutions. The airborne launch vehicle application to launch a small satellite with the purpose of remote sensing requires high precision exit on specified sun-synchronous orbit. So a problem is stated to construct an optimal ascent trajectory and optimal control.
In this paper, the mathematical motion model of the air-launched launch vehicle with the external disturbances caused by the Earth's non-sphericity, drag and wind is put forward based on the three-stage flight program with passive intermediate section. A discrete process based on pseudo-spectral method is used to solve the problem, which allows converting the initial problem into a nonlinear programming problem with dynamic constraints and aims for the criteria of maximization of the final mass released onto the target orbit.
Application of the proposed solution procedure is illustrated by calculating the optimal control and the corresponding trajectory for two-stage liquid launch vehicle, which places the small spacecraft on the orbit of sun-synchronous at the height of 512 km. The numerical simulation results have demonstrated the effectiveness of the proposed algorithm and allow us to analyze three-stage trajectory parameters with intermediate passive flight phase. It can be noted that in the resulting ascent trajectory, the intermediate passive flight part is a suborbital trajectory with low energy integral, perigee of which is under the surface of the Earth.
References
1. Pyshnyy I.A., Chepiga V.E. Zapusk malykh iskusstvennykh sputnikov Zemli s ispol'zovaniem samoletov-nositeley [Launch of small artificial satellites using aircraft carrier]. Moscow, Mashinostroenie Publ., 2005. 168 p. (in Russian).
2. Pegasus user's guide. Orbital: company website. Available at: http://www.orbital.com , accessed 16.10.2014.
3. Alifanov O.M., Andreev A.N., Gushchin V.N., Zolotov A.A., Matveev Iu.A., Perelygin V.P., Khokhulin V.S. Ballisticheskie rakety i rakety-nositeli [Ballistic missiles and launch vehicles]. Moscow, Drofa Publ., 2004. 512 p. (in Russian).
4. Mazgalin D.V. Voprosy postroeniya programmnoy traektorii vyvedeniya rakety-nositelya s kosmicheskim apparatom. Kand. dis. [Program trajectory construction problems for launch vehicle with spacecraft. Cand. dis.]. Ekaterinburg, 2012. 136 p. (in Russian).
5. Benson D. A gauss pseudospectral transcription for optimal control: Ph.D. thesis. Cambridge, 2005. 243 p.
6. Vallado D.A. Fundamentals of astrodynamics and applications. Microcosm Press, New York, Springer, 2007. 1055 p.
7. GOST 4401-81. Atmosfera standartnaya. Parametry [State Standard 4401-81. Standard atmosphere. Parameters]. Moscow, Standards Publishing House, 1981. 180 p. (in Russian).
8. OST 92-9704-95. Rakety i rakety-nositeli. Modeli opredeleniya gorizontal'noy skorosti vetra i termodinamicheskikh parametrov atmosfery v diapazone vysot 0-120 km v rayone kosmodroma Plesetsk [Industry Standard 92-9704-95. Missiles and launch vehicles. Models for determining the horizontal wind speed and the thermodynamic parameters of the atmosphere at the altitude of 0-120 km in the launch site Plesetsk]. Moscow, 1995. 94 p. (in Russian).
9. Isaev V.K. Pontryagin's maximum principle and optimum programming of rocket thrust. Avtomatika i telemekhanika, 1961, vol. 22, no. 8, pp. 986-1001. (in Russian).
10. Garg D., Hager W.W., Rao A.V. Pseudospectral methods for solving infinite-horizon optimal control problems. Automatica, 2011, vol. 47, no. 4, pp. 829-837. DOI: 10.1016/j.automatica.2011.01.085
11. Biegler L.T., Zavala V.M. Large-scale nonlinear programming using IPOPT: An integrating framework for enterprise-wide dynamic optimization. Computers & Chemical Engineering, 2009, vol. 33, no. 3, pp. 575-582. DOI: 10.1016/j.compchemeng.2008.08.006
12. Gill P.E., Murray W., Saunders M.A. SNOPT: An SQP algorithm for large-scale constrained optimization. SIAMReview. 2005. Vol. 47, no. 1. P. 99-131. DOI:10.1137/S0036144504446096
13. Ingrid M.D. Modeling dispersions in initial conditions for air-launched rockets and their effect on vehicle performance: Master's thesis. Cambridge, 2013. 159 p.
