CC BY
89
УДК 629.191
МЕТОД ОПЕРАТИВНОГО РЕШЕНИЯ
НЕЛИНЕЙНЫХ КРАЕВЫХ ЗАДАЧ ДИНАМИКИ ПОЛЕТА
НА ОСНОВЕ ЭНЕРГЕТИЧЕСКОГО ПРИНЦИПА
И. О. Петров,
канд. техн. наук, доцент
Санкт-Петербургский государственный университет аэрокосмического приборостроения
Рассмотрен новый метод оперативного решения краевых задач динамики полета на основе энергетической теории маневрирования. Предложенный метод не требует решения сложных и неустойчивых вариационных задач и может использоваться для решения краевых задач в темпе реализации маневра.
Ключевые слова — маневры летательного аппарата, оперативное оптимальное управление, энергетическая теория, краевые задачи, граничные условия.
Введение
Поиск оптимальных траекторий полета любых летательных аппаратов (ЛА) в атмосфере на базе традиционных методов оптимизации с математической точки зрения представляет довольно сложную и трудоемкую задачу, обусловленную необходимостью решать вариационные краевые задачи, близкие, как правило, к некорректным из-за отсутствия решения даже при использовании современных ЭВМ и методов вычислительной техники [1-9]. Основные сложности вызваны в первую очередь тем, что для сопряженной системы дифференциальных уравнений не удается получить все необходимые интегралы, и возникают большие трудности с поиском начальных условий интегрирования данной системы. Указанные начальные условия (или их часть) обычно определяют методом простого подбора. Это ставит под большое сомнение оптимальность управления, полученного с помощью принципа максимума Л. С. Понтрягина [10, 11], так как оно напрямую зависит от сопряженных переменных, содержащих в себе волюнтаризм подбора начальных условий интегрирования сопряженной системы дифференциальных уравнений. Кроме этого, «оптимальные» траектории зависят от способа (аналитического или численного) вычисления каждой производной в системе сопряженных переменных. Поэтому в общем случае можно говорить
о некоторых «попадающих» в заданную точку
траекториях, полученных с помощью принципа максимума Л. С. Понтрягина, но не об оптимальных траекториях [1].
В связи с этим в данной статье рассмотрен новый подход к решению краевых задач динамики полета на основе энергетической теории маневрирования без решения сложных и неустойчивых краевых вариационных задач.
Постановка задачи выведения воздушно-космического аппарата
Предметом исследования в статье является управляемое движение центра масс ЛА на участке выведения. В качестве объекта исследования выбран воздушно-космический аппарат (ВКА) авиационного базирования, являющийся второй ступенью двухступенчатой многоразовой воздушно-космической системы (МВКС) 2-го типа на основе авиационно-космического комплекса. Данный выбор объясняется тем, что МВКС 2-го типа являются переходным этапом от существующих МВКС к перспективным одноступенчатым МВКС. Создание таких МВКС реально уже в ближайшие годы, поэтому для проверки работоспособности созданных моделей, методов и алгоритмов целесообразно ориентироваться именно на этот класс МВКС, как наиболее исследованный.
Общая схема выведения полезной нагрузки (ПН) на заданную орбиту с помощью МВКС 2-го типа включает в себя следующие этапы:
1) старт самолета-носителя с ВКА с обычной взлетно-посадочной полосы и вход в заданную плоскость орбиты;
2) набор высоты и скорости самолетом-носителем при полете в плоскости орбиты с целью обеспечить требуемые начальные условия старта для ВКА в момент разделения;
3) разделение ВКА и самолета-носителя;
4) возвращение самолета-носителя к месту старта и посадка на взлетно-посадочную полосу;
5) выведение ВКА вместе с ПН на заданную орбиту.
В данной статье рассматривается выведение ВКА от момента отделения его от самолета-носителя до момента выхода его на заданную орбиту.
Из анализа работ, посвященных теории оптимизации маневров различных ЛА [1-12], следует вывод, что при выведении ВКА в идеальном случае надо все топливо сжигать взрывом с целью обеспечить максимальную тягу двигателей и сокращение потерь характеристической скорости. Но из-за технических ограничений, накладываемых на ВКА, сделать это невозможно. Поэтому в течение всего участка выведения с непрерывной тягой тяга двигателей должна быть максимально возможной или изменяться по какому-либо жестко заданному закону, учитывающему конкретные технические ограничения, накладываемые на ВКА.
С учетом этого получается эквивалентность следующих задач выведения:
— минимизации потерь характеристической скорости;
— минимизации времени выведения (задачи на быстродействие);
— минимизации расхода топлива;
— максимизации массы выводимой ПН.
Это видно из следующих несложных рассуждений. Абсолютные конечная Ка.к и начальная Уа0 скорости, а также характеристики двигательной установки ВКА однозначно определяют минимальное время работы двигателей для достижения требуемой скорости Ка.к при отсутствии потерь характеристической скорости:
г =
то
- _-(Ц.е-Уа0)/Щь
(1)
где то — начальная масса ВКА; тсек — массовый секундный расход двигателей ВКА; Жэф — эффективная скорость истечения двигателей ВКА.
Но поскольку потери характеристической скорости неизбежны, то за время ВКА не достигнет требуемого конечного значения скорости Ка к. За это время потери характеристической скорости
д£ї =/ Д ^аЕЇ ^ г-
Го
(2)
Под ДЦ,понимается сумма потерь абсолютной скорости, расходуемая на преодоление сопротивления атмосферы, гравитационного поля Земли, управление ВКА, а также учитывающая влияние атмосферы на величину тяги. С учетом сказанного достигнутая ВКА за время фактическая абсолютная скорость
Чб(г ) = Чё - ДЦ> „ ■
(3)
Для того чтобы доразогнать ВКА до скорости Ка к, требуется дополнительное время Д^доп, определяемое величиной суммарных потерь скорости
Д^а = V
аЕї аЕї
/ Д
(4)
Второе слагаемое в (4) определяет потери на астке доразгона. Итак, вре делить следующим образом:
участке доразгона. Итак, время Д^доп можно опре-
ДгаТї =
т(г )
те,.
(5)
Общее время полета ВКА до момента набора заданной скорости К, к
(6)
Из выражения (5) следует, что при уменьшении суммарных потерь характеристической скорости Д^аЕ | уменьшается дополнительное время полета Д£доп. Следовательно, задача минимизации потерь характеристической скорости на траектории выведения равносильна задаче минимизации дополнительного времени полета Д£доп-А так как время является константой при заданных величинах Ка0, Уа.к, Ж^, тсек, т0, то задача минимизации потерь характеристической скорости равносильна задаче на быстродействие.
При заданном законе изменения тяги двигателей ВКА минимизация времени выведения #выв приводит к минимизации расхода топлива тт, а это, в свою очередь, к максимизации выводимой массы ПН тпн при условии постоянства начальной массы т0.
Математическую эквивалентность функционалов (функций целевого эффекта операции) задач можно записать таким образом:
ттДVa .. ^ т1п тттй тахту [. (7)
иШ Е| иШ иШ иШ
хеХ хеХ хеХ хеХ
С учетом эквивалентности функционалов (7) целесообразно сформулировать и решить задачу на быстродействие, так как она проще реализуется математически.
В качестве управления при выведении ВКА принимается вектор и = [у, К, иЙ ] , где у — ско-
ростной угол крена; к — полный угол атаки; ид — угол установки двигателей на ВКА.
Фазовые координаты ВКА в момент отделения его от самолета-носителя обозначаются через x(to) = xo, параметры орбиты выведения обозначаются через x(tg) = xe. Точка xo представляет собой начальные условия движения, а точка xR — граничные условия выведения.
Формулировка задачи выведения
Определить оптимальное управление в виде
u = [у, к, Ua ] , переводящее ВКА из точки xo при t = to в точку xR при t = tR, удовлетворяющее дифференциальным связям x = f(x, u) выбранной математической модели, ограничениям
O(t) Є [Omin, Omax] (8)
и доставляющее минимум функционалу
J = !aua. (9)
В выражении (8) O(t) — вектор конструктивных и других параметров, относящихся к самому ВКА, реализуемому управлению u(t) и траектории выведения x(t). Векторы Omin и Omax представляют собой соответственно минимальные и максимальные пределы, в которых должны изменяться физические величины вектора O(t). Другими словами, векторы Omin и Omax — это технические ограничения на ВКА и управление, а также ограничения на траекторию выведения.
Конкретную математическую модель движения центра масс ВКА целесообразно записать во вращающейся системе отсчета, орты которой представляются следующим образом:
xO = Г0; zO = jliVL; yO = zO xx0. (lo)
|rxVa|
С точки зрения простоты и физического толкования реального процесса движения ВКА в атмосфере на участках выведения и спуска удобнее использовать модель в следующем виде [1]:
* 2 1 / \
Vr = ax + ozr; г = Vr; соz = —(y — 2rnzVr);
Г = ox — ю^дфcosr; ф = rozsinr;
І = raz^°sr — Qc, (11)
cos ф
о az
где юх = -^-.
raz
Представленная модель имеет правые части, наиболее близкие к линейным без проведения линеаризации, что позволяет увеличить шаг интегрирования без потери заданной точности расчетов.
Проекции абсолютного ускорения ax, ay, az на оси вращающейся системы отсчета определяют-
ся с учетом вращения атмосферы совместно с Землей:
ax =—(—G + Fx sin 9Г + Fy cos у cos 9r);
f'ff' v У v !
1
ay = — (( cos 9r — Fy^ cos у sin 9r )■ cos(r|r — r) +
+ ~ Fyv sin у sin(rr—r);
m Jv
1
az =—cos 9r — Fy^ cos у sin 9r )■ sin(rr — r) —
—~ Fyv sin у cos(rr—r);
mv
Fx =—Q + Pcos up,
v
Fyv = Y + Psin up, ug = к + ua. (12)
В выражениях (12) силы, входящие в правую часть, вычисляются следующим образом:
2 2 G = ^, Q = CxP^ Sm, Y = CyP^ Sm, (13)
Г 2 2
где р — плотность воздуха; Ув — модуль воздушной скорости; Sm — площадь миделя; cx — коэффициент силы лобового сопротивления; cy — коэффициент подъемной силы.
Следует заметить, что аэродинамические силы будут ориентированы в скоростной (относительной) системе координат, что обусловлено вращением атмосферы вместе с Землей.
Аэродинамические коэффициенты целесообразно представить в виде полиномов второй степени:
cx = kQ к2 + kQ2 к + kQ3;
Cy = kY к2 + k.Y2 к + kY3, (14)
где коэффициенты kq , kQ2, k^, k^i, ky^, ky^ являются функциями числа Маха M.
Полиномы (14) позволяют сгладить табличные данные и уменьшить погрешности измерений, допущенные при проведении продувок модели ВКА в аэродинамических трубах. Кроме этого, полиномы (14) позволяют получать частные производные от аэродинамических сил по полному углу атаки в аналитическом виде при определении оптимального управления.
Полиномы (14) являются универсальной формой представления зависимостей cx и cy от полного угла атаки к. Например, в случае нелинейной зависимости cx = cx (к) коэффициент kQ ^ О, а коэффициенты kQ2, kq могут принимать любые значения; в случае линейной зависимости kQ = О, kQ2 ^ О, kq — произвольный; в случае отсутствия зависимости cx = ^(к) коэффициенты kQ = О, kQ2 = 0 а коэффициент kQ3 = cXo (M), где Cx0 (M) —
коэффициент силы лобового сопротивления при нулевом угле атаки.
В качестве начальных условий движения ВКА приняты фазовые координаты в момент отделения ВКА от самолета-носителя в виде
хо = ['о 140, 0го, Пг0, Фо, ^0 ]Т (15)
Граничные условия на правом конце траектории можно записать в виде вектора
Т
xe = lre, Va.e, 0ё|
(16)
Для того чтобы получить начальные условия для интегрирования системы дифференциальных уравнений во вращающейся системе отсчета, необходимо с использованием координат вектора x0 в виде (15) рассчитать фазовые переменные Vr(, ю, , -По:
Vr = V, sin(
ro э0
г0'
ю, = z0
^ cos0ro cosrr0 + r0^c cosф0
2
+ (( cos0ro sin Гг0 )
Го = 2л — arctg
r0
— V, cos0r sinrr d} r0 {r0
Vs0 cos0ro cosrr0 + ro^c cosфо
(17)
(18)
(19)
где операция вычисления арктангенса соответствует круговому арктангенсу и позволяет вычислять значение угла в пределах от 0 до 360°.
В результате преобразований начальные условия (15) можно переписать в виде
х0 =[ \ , r0, ю2о, ^ Ф0, ^0] . (20)
Метод определения структуры оперативного оптимального управления движением центра масс ЛА
Поскольку начальные (20) и граничные (16) условия определяют потребную величину полной механической энергии ДЕа [1, 4], которую необходимо набрать на участке выведения, то для минимизации времени выведения следует максимизировать скорость набора полной механической энергии в абсолютном движении Еа.
На основании проведенных исследований на экстремум по управлению скоростей изменения полной механической энергии в относительном (Ег) и абсолютном (Еа) движениях [4] получено оперативное оптимальное управление для относительного и абсолютного движения, максимизирующее мощность поверхностных сил. В результате анализа полученных управлений струк-
туру оперативного оптимального адаптивного автономного управления в общем виде можно записать следующим образом:
tg у =
Py cos(rr — r)— P0 sin(rr — г)
P°. sin 0r sin(rr — r) + 0r cos(rr — r)— р°я
Py sin 0r cos
Ua = arctg
cos0
N2
N1
— к;
к=—
kQ2 N1
kY N2
2(kq N1
kY N2
(21)
(22)
(23)
где
N1 =— P°3 sin0r — Py cos0r cos(rr — r) — — Pycos 0r sin(rr—r);
N2 = P°3 cos у cos 0r — Pyg x
(24)
x[cos у sin 9r cos(rr — r) — sin у sin(rr — r)] —
— p0 [cos у sin 9r sin(rr — r) + sin у cos(rr — r); (25)
P°a, P0, P°a — проекции орта вектора импульсов p0 единичной массы на оси вращающейся системы отсчета (ВСО), относительно которой максимизируется мощность поверхностных сил.
Вектор p0 колинеарен либо вектору воздушной скорости, т. е. p0 = V°, либо вектору абсолютной скорости, т. е. p0 = V°. В первом случае управление (21)-(25) обеспечивает максимизацию мощности поверхностных сил в относительном движении, во втором случае — в абсолютном движении.
Структура управления в виде (21)-(25) при условии p0 = V° полностью идентична оптимальному управлению, полученному из решения аналогичной вариационной задачи на основании принципа максимума Л. С. Понтрягина. Характерной особенностью данной структуры является отсутствие необходимости интегрировать сопряженную систему уравнений, так как все величины, входящие в состав управления (21)-(25), известны в каждой точке траектории.
В результате исследований выяснено, что при решении задач динамики полета брать за основу относительное движение нецелесообразно, так как при этом получается неопределенность управления. Продемонстрируем это. Если в управлении (21)-(25) вместо проекций вектора p0 на оси вращающейся системы отсчета подставить проекции орта вектора воздушной скорости V°, т. е.
„О ■ п . „О
Py =sin0r; Py, = cos0rcos(rr—r);
Py = cos 0r sin(rr — r),
(26)
№ 1, 2013
йи^ 37
то после преобразований
о
їду=о; 19ир=°;
кОг
к = — , _ 2 ; иа =—к,
2к
(27)
0\
так как N1 = -1, N2 = 0.
Таким образом, вектор тяги двигателей ВКА должен совпадать с направлением вектора воздушной скорости, угол К обеспечивает минимальный ущерб от аэродинамической силы, т. е. ее максимальную мощность в относительном движении, а скоростной угол крена у остается неопределенным и может принимать любые значения, при этом мощность поверхностных сил не будет изменяться.
Однако, задавая различные значения скоростного угла крена, в абсолютном пространстве можно получать самые различные траектории. Задание любой конкретной программы движения у = у(Ь) в этом случае равносильно заданию начальных значений сопряженных переменных при использовании принципа максимума Л. С. Понтрягина. Отсюда ясно, почему процесс решения вариационных задач с использованием принципа максимума является неустойчивым и неоднозначным.
Для устранения данной неоднозначности необходимо максимизировать мощность поверхностных сил в абсолютном движении, где управление полностью определено. Для этого случая проекции вектора импульсов определяются следующим образом:
р°а = ^п0=уа; Руа = ^0=Чт-; р?а = (28)
V- о
ГЮ-
о
V,
Следует отметить, что разработанное управление (21)-(25), максимизирующее мощность поверхностных сил, не определяет вид траектории, по которой будет двигаться ВКА. Вид траектории будет зависеть от начальных условий старта, характеристик двигательной установки и аэродинамического облика ВКА, ограничений, накладываемых на траекторию и управление, и т. д. Поэтому если требуемые граничные условия не лежат на траектории максимального набора энергии, то уточнение управления производится путем решения специфической нелинейной краевой задачи. Специфика заключается в том, что подбираются параметры управления, определяющие пространственную ориентацию вектора импульсов p0 при выведении ВКА и обеспечивающие удовлетворение заданных граничных условий на правом конце траектории. Количество параметров и их физический смысл зависят от конкретной задачи выведения.
Метод определения параметров управления
Прежде чем приступить к описанию метода решения специфических нелинейных краевых задач выведения, необходимо определить его место среди разнообразных методов выбора оптимального программного управления. Принципиально различают два этапа разработки управления движением ЛА [1, 7]:
1) баллистическое проектирование маневров ЛА;
2) практическую реализацию маневров ЛА.
В данной статье рассматривается этап баллистического проектирования маневров ЛА на участке выведения. На этом этапе применительно к задачам оптимального управления динамическими системами отмечают две основные группы методов — универсальные и широкоспециализированные.
В рамках первой группы выделяют прямые и непрямые методы. Подробная классификация и анализ существующих традиционных методов решения задач оптимального управления приведены в работе [1].
Предлагаемый метод относится к группе непрямых методов, основанных на необходимых условиях оптимальности, и включает в себя три этапа определения управления:
1) определение вида (структуры) оптимального управления;
2) определение количества и состава параметров управления;
3) решение собственно краевой задачи.
Для решения поставленной задачи выведения ВКА используется полученная структура управления (21)-(25) с учетом (28). Структура оперативного оптимального адаптивного автономного управления для абсолютного движения записывается в виде
0
-РУа ^п(|г-г)
руа 0гСО^Гг - Л)- Рха cos0Г
иа = aгctg
N2
N1
к=-
ко2 N1 + ку2 N2 2( N1 + ку N2
(29)
(30)
(31)
где
N1 = — р<оа sin0Г — руа cos0Гcos(rlГ — л); (32)
N2 = р^а cosуcos0Г — р°уа х x[cosуsin0Гcos(rГ — г) — sinуsin(rГ — г); (33)
.о „о
Рха , Р
Уа
проекции орта вектора импульсов p
0
единичной массы на оси ВСО, определяемые по формулам (28).
Управление (29)-(33) является оперативным, так как вычисляется по конечным простым формулам только через текущие фазовые координаты, которые известны в каждый момент времени. Это управление по своей структуре идентично структуре оптимального управления, полученного из решения аналогичной вариационной задачи. Поэтому при использовании управления (29)-(33) отпадает необходимость решать сложные вариационные задачи, не требуются поиск начальных значений и интегрирование сопряженной системы дифференциальных уравнений.
Для удовлетворения граничных условий на правом конце траектории выведения необходимо управлять ориентацией вектора импульсов p0 в процессе маневра в целях выполнения двух требований:
— отсутствия падения высоты ВКА в начале активного участка;
— наведения ВКА в конце активного участка.
Для этого необходимо определить количество
и состав параметров управления, определяющих поведение вектора p0 в течение активного участка.
Из анализа управления (29)-(33) видно, что оно обеспечивает ориентацию вектора тяги двигателей по вектору абсолютной скорости, так как p0 = V0. В этом случае при рассмотрении движения во вращающейся системе отсчета если сумма вертикальной составляющей вектора тяги двигателей и вектора переносной силы инерции будет по модулю меньше, чем модуль вектора силы притяжения, то ВКА начнет через некоторое время терять высоту даже при достаточно больших углах 00 наклона вектора начальной абсолютной скорости Ка0 к линии местного горизонта. Так как в момент старта ВКА с дозвукового самолета-носителя модуль переносной силы инерции мал по сравнению с модулем силы притяжения, а максимальное значение 00 невелико и ограничено возможностями самолета-носителя, то компенсацию силы притяжения необходимо осуществлять за счет отклонения вектора тяги двигателей ВКА от оптимального направления. Следовательно, надо изменить проекции вектора p0 на оси ВСО, которые задают ориентацию вектора тяги двигателей. Управляя ориентацией вектора p0 в конце активного участка, можно осуществлять наведение ВКА в целях удовлетворения граничных условий. Таким образом, из приведенных рассуждений наглядно видно противоречие между требованием оптимальности управления с точки зрения максимизации мощности поверхностных сил в абсолютном движении и требованием удовлетворения граничных условий на правом конце траектории. При выведении ВКА необходимо найти такое управление ориентацией вектора p0, чтобы обеспечить удовлетворение
граничных условий при минимальном отклонении от режима полета с максимальным набором энергии.
Предлагается следующий алгоритм управления ориентацией вектора p0, удовлетворяющий предъявленным требованиям:
pXa = sin9p + во; pya = ^1 —(sin9p+во) ;
рО = О, (34)
' za
где 9р — угол наклона вектора p0 к линии местного горизонта;
Isin9 'Эе Sр <9 sin 9 р = , (35)
Р |sin»ра 'Эе Sра >9
2
Vr
где 0 = arctg
местного горизонта; S Pa = arcsin
V0
G — Fe
угол
Ра
тангажа тяги двигателей, обеспечивающий ком-
2
пенсацию силы притяжения G, Ре = тгю2 — модуль вектора переносной силы инерции Fe;
[о ї5е Г < Г,
єо =] ' , (36)
[єдаааТ г > гї
где є?аааТ е[—1 — sin0p, 1— :5т.р
Как видно из формул (34)-(36), управление ориентацией вектора p0 осуществляется с помощью двух параметров, которые называются время полета в режиме максимального набора энергии с учетом ограничений на траекторию — £п и отход от режима максимального набора энергии — Є0.
Параметр £п очень сильно влияет на величину радиуса гк в конце активного участка, а параметр Є0 — на угол наклона 0к вектора Vа к линии местного горизонта. Заданное значение модуля абсолютной скорости У,ё в конце активного
а- сдаааГ
участка физически обеспечивается выключением двигателей ВКА, когда требуемая скорость набрана, а при моделировании на ЭВМ — окончанием процесса интегрирования.
Таким образом, вектор уточняемых параметров управления можно записать в виде
С = [, єо ]Т, (37)
а вектор заданных граничных условий
іТ
ёдаааГ ' ёдаааГ ' а-ёдаааГ
xraaaf = I re , 0e , Vae I . (38)
* ^caaaf ^caaaf а,сс?5е|г 1
Теперь, варьируя значениями координат вектора ^ можно обеспечить заданные значения координат вектора xзадан, одна из которых удовлетворяется автоматически, а именно модуль Уаё.....
№ 1, 2013
39
Заключение
Предложенный метод оперативного решения специфических нелинейных краевых задач динамики полета ЛА позволяет использовать его для расчета управления движением центра масс непосредственно на борту в темпе полета, так как управление рассчитывается по конечным формулам на основании известных в каждой точке тра-
Литература
1. Насонов В. П. Нетрадиционный подход к решению традиционных задач динамики полета ракет-носителей / ВИКИ им. А. Ф. Можайского. — СПб., 1992. — 64 с.
2. Насонов В. П. Выбор программ движения ракет / МО СССР. — М., 1976. — 89 с.
3. Петров И. О. Проблема определения оперативных оптимальных автономных алгоритмов управления ЛА и пути ее решения // Аэрокосмическое приборостроение. 2012. № 3. С. 14-20.
4. Петров И. О. Энергетическая теория синергетического маневрирования // Аэрокосмическое приборостроение. 2012. № 4. С. 10-20.
5. Петров И. О., Насонов В. П. Принцип применения энергетической теории к определению оптимальных программ движения ЛА относительно центра масс // Аэрокосмическое приборостроение. 2012. № 7. С. 3-11.
6. Петров И. О. Математические модели движения летательных аппаратов во вращающейся атмосфере Земли // Аэрокосмическое приборостроение. 2012. № 9. С. 10-20.
ектории фазовых координат. На основании этого могут быть созданы различные системы терминального наведения в условиях действия возмущающих факторов.
Кроме этого, разработанный метод создает единую методологическую базу для исследования маневров ЛА в атмосфере на участках выведения и спуска без изменения полученной структуры управления и без решения вариационных задач.
7. Аверкиев Н. Ф., Волков В. Ф., Петров И. О. Баллистическое проектирование РН: учеб. пособие / ВИКУ им. А. Ф. Можайского. — СПб., 1999. — 72 с.
8. Шкадов Л. М., Буханова Р. С., Илларионов В. Ф., Плохих В. П. Механика оптимального пространственного движения летательных аппаратов в атмосфере. — М.: Машиностроение, 1972. — 240 с.
9. Школьный Е. П., Майборода Л. А. Атмосфера и управление движением летательных аппаратов. — Л.: Гидрометеоиздат, 1973. — 308 с.
10. Беллман Р., Калаба Р. Динамическое программирование и современная теория управления. — М.: Наука, 1969. — 119 с.
11. Понтрягин Л. С., Болтянский В. Г., Гамкрелидзе Р. В., Мищенко Е. Ф. Математическая теория оптимальных процессов. — М.: Наука, 1983. — 392 с.
12. Баринов К. Н., Насонов В. П. Краевые задачи динамики полета космических аппаратов / МО СССР. — М., 1970. — 211 с.
