Прикладная математика и механика, 2023, T. 87, № 4, стр. 631-641

1. Введение. Поиск оптимального решения в рамках математической модели, учитывающей множество ограничений, которые связаны с техническими особенностями определенного типа устройств – сложная задача. Можно выделить два основных подхода к ее решению. Первый основан на принципе максимума Л.С. Понтрягина [1], а второй – на методе динамического программирования Р. Беллмана [2]. Первый подход, как правило, оказывается более сложным с точки зрения применяемой математической теории, но требует меньше вычислительных ресурсов, поскольку предполагает решение краевой задачи применительно к системе обыкновенных дифференциальных уравнений, для чего существует много разнообразных численных методов. Второй подход проще с точки зрения математического содержания, но для его использования в большинстве случаев необходимы значительные вычислительные мощности, поскольку в своей основе он представляет собой метод перебора всех возможных решений. Однако с ростом числа и сложности ограничений ситуация меняется, потому что для применения принципа максимума в этом случае нужны чрезвычайно громоздкие системы уравнений, а в рамках метода динамического программирования увеличение количества ограничений приводит, вообще говоря, к существенному уменьшению количества тех вариантов решений, среди которых нужно делать выбор. Следовательно, целесообразно применять оба подхода, опираясь на свойства конкретных решений. Далее на примере вычисления оптимальных по расходу топлива траекторий неманевренного сверхзвукового летательного аппарата будет показано, как можно реализовать такую комбинацию.

В настоящее время применяют, в основном, локальную оптимизацию отдельных участков полета самолета [3]. В качестве таковых часто выделяют взлет, крейсерский режим и посадку [4, 5]. Далее решение для всего перелета формируют на основе комбинации из этих участков [6, 7], для чего используют различные численные методы [8, 9]. С увеличением мощности бортовых компьютеров появилась возможность оптимизировать траекторию по расходу топлива во время полета, учитывая текущие воздушные потоки [10, 11]. В данном исследовании воспользуемся результатами работы [12], где была реализована глобальная оптимизация всей траектории целиком, без выделения отдельных участков.

При поиске параметров движения гражданского самолета необходимо учитывать много ограничений [13], часть из которых основана на технических характеристиках конкретного аппарата, а часть – на особенностях управления воздушным движением в целом. Влияние последних хорошо заметно на примере обычной формы траектории [4] для дозвукового воздушного судна. Эта форма существенно отличается от той, которая необходима для экономии топлива (см., например, [14]) и показана на рис. 1, хотя такая экономия исключительно важна с коммерческой точки зрения. Причина в том, что для снижения расхода горючего необходимо, как видно из рис. 1, набрать максимально возможную высоту сразу после взлета и затем постепенно снижаться по мере приближения к аэропорту назначения, а такое снижение повышает вероятность столкновений воздушных судов друг с другом. По этой причине реальные полеты проходят, как правило, на одной из заранее оговоренных высот, а переходы с одной такой высоты на другую редки. Выяснилось [12], что в случае сверхзвукового гражданского лайнера траектории, оптимальные по расходу топлива, весьма похожи на традиционные траектории дозвуковых воздушных судов, а потому их можно использовать на практике. Следовательно, изучение таких решений актуально.

2. Уравнения движения. Численные значения большинства используемых далее постоянных, описывающих математическую модель сверхзвукового самолета, приведены в работе [12]. В данной статье они опущены для краткости изложения.

Ограничимся рассмотрением траекторий движения самолета, полностью лежащих в вертикальной плоскости. Обозначим через $m$ массу летательного аппарата в текущий момент времени. Пусть $x$ – горизонтальная, а $y$ – вертикальная координаты центра масс аппарата в неподвижной системе отсчета, связанной с землей. Модуль вектора скорости этого центра обозначим через $V$, а величину угла между горизонтальной осью и этим вектором – через $\theta $. Кроме того, пусть данный вектор всегда коллинеарен вектору силы тяги, имеющему модуль $P$. Далее введем вектор, равный отношению суммы векторов тяги и полной аэродинамической силы к величине силы тяжести. Его проекцию на направление вектора скорости центра масс самолета, называемую тангенциальной перегрузкой, обозначим через ${{n}_{x}}$, а на ось, ортогональную вектору скорости и направленную к верхней части самолета – через ${{n}_{y}}$. Эту проекцию называют нормальной скоростной перегрузкой. Тогда уравнения движения аппарата можно записать в форме [15]:

В соотношениях (2.1) функция ${{Q}_{t}}$ от тяги $P$, числа Маха $M$ и высоты полета $y$ определяет величину расхода топлива за секунду, а число $g$ равно модулю ускорения свободного падения. Число Маха $M$ определено как $M = V{\text{/}}{{V}_{*}}$, где скорость звука ${{V}_{*}}$ зависит от высоты $y$ известным образом.

Пусть $S$ – площадь крыла самолета, а $\rho = \rho (y)$ – зависимость плотности атмосферы от высоты. Тогда

причем в первой формуле из соотношений (2.2) учтено, что вектор скорости всегда направлен против силы лобового сопротивления. Коэффициент лобового сопротивления ${{C}_{x}}$ в соотношении (2.2) зависит от числа Маха $M$, а также коэффициента подъемной силы ${{C}_{y}}$ следующим образом:

В соотношении (2.3) использованы известные постоянные величины ${{C}_{{y0}}}$, ${{k}_{y}}$, ${{M}_{0}}$, ${{k}_{M}}$, ${{c}_{{ij}}}$ и ${{b}_{{ij}}}$.

Функция

задает скорость расхода топлива, причем величины ${{P}_{q}}$, ${{M}_{q}}$, ${{y}_{q}}$, ${{\zeta }_{P}}$, ${{\zeta }_{M}}$, ${{\zeta }_{y}}$, а также ${{\chi }_{{kij}}}$ – известные константы.

Необходимо найти функции ${{C}_{y}}(t)$ и $P(t)$, доставляющие минимум функционалу

где величина $T$ равна заданному времени движения от известной начальной точки до известной конечной, причем начальная и конечная скорости центра масс самолета тоже заданы. Заметим, что значение функционала (2.5) равно массе израсходованного топлива за все время полета.

Заменим в (2.1) независимую переменную $t$ на дальность $x$. Это допустимо, поскольку горизонтальная проекция $Vcos\theta $ скорости центра масс неманевренного самолета никогда не равна нулю. Имеем

Тогда граничные условия можно записать в форме

где ${{x}_{0}}$ – начальная горизонтальная координата центра масс самолета, а ${{x}_{T}}$ – конечная. В дальнейшем будем полагать ${{x}_{0}} = 0$, так как движение летательного аппарата не зависит от положения начала системы координат на оси абсцисс. Функционал (2.5), основываясь на соотношениях (2.6), перепишем в виде

Искомые управления ${{C}_{y}}(x)$ и $P(x)$, доставляющие минимум функционалу (2.8), были получены в статье [12] с помощью классического варианта метода динамического программирования Беллмана. Подробное описание этого варианта применительно к рассматриваемой задаче можно найти в работе [16]. Также были найдены соответствующие оптимальные траектории, удовлетворяющие дифференциальным уравнениям (2.6) с начальными условиями (2.7). В процессе вычислений были учтены сложные и разнообразные фазовые и иные ограничения, которые будут описаны ниже.

3. Ограничения. Разрешенные значения коэффициента подъемной силы ${{C}_{y}}$ принадлежат отрезку $0 \leqslant {{C}_{y}} \leqslant {{C}_{{ymax}}}(M)$, где

В формулах (3.1) величины ${{C}_{{y1}}}$, ${{M}_{c}}$, ${{M}_{{ci}}}$ известны.

Допустимые значения величин $y$ и $V$ образуют замкнутую область, причем высота полета $y$ не должна превышать 14 000 м и не должна быть меньше 100 м, а модуль вектора скорости центра масс самолета $V$ ограничен снизу функцией ${{V}_{{min}}} = {{V}_{{min}}}(y)$, а сверху – функцией ${{V}_{{max}}} = {{V}_{{max}}}(y)$:

где ${{y}_{0}}$, ${{H}_{i}}$ и ${{h}_{j}}$ – известные постоянные величины.

Величина $P$ силы тяги двигателей ограничена снизу и сверху функциями ${{P}_{{min}}}$ и ${{P}_{{max}}}$ соответственно. Они определены следующим образом:

где ${{\lambda }_{M}}$, ${{M}_{p}}$, ${{\lambda }_{y}}$, ${{y}_{p}}$, $\xi _{{0i}}^{*}$, $\xi _{{1i}}^{*}$, ${{\xi }_{{ij}}}$, $\eta _{{ji}}^{*}$ и ${{\eta }_{{ij}}}$ – известные константы.

Значение угла наклона траектории $\theta $ должно удовлетворять ограничению $ - 45^\circ \leqslant \theta \leqslant 45^\circ $.

Величина нормальной скоростной перегрузки должна быть неотрицательной и не превышать значения $n_{y}^{{max}} = 4$.

Заметим, что соотношения (2.3), (2.4), (3.1)–(3.3) представляют собой численные аппроксимации экспериментальных данных.

4. Оптимальные траектории. Возьмем десять характерных оптимальных траекторий, отличающихся друг от друга длительностью полета и соответствующих граничным условиям (2.7): ${{V}_{0}} = 140$ м/с, ${{\theta }_{0}} = 0$, ${{m}_{0}} = 6 \times {{10}^{4}}$ кг, ${{x}_{T}}{{ = 10}^{6}}$ м, ${{V}_{T}} = 140$ м/с, ${{\theta }_{T}} = 0$. Вычисления были проведены для площади крыла, равной $S = 110.16$ м².

На рис. 2 изображены зависимости количества времени, необходимого для удаления от начальной точки на определенное расстояние, от этого расстояния. Видно, что длительность полета находится в диапазоне от 48 до 58 мин. Чем выше находится кривая, тем больше экономия топлива на описываемой ею траектории. Заметим, что на рис. 2 показаны только находящиеся вблизи конечной точки участки графиков, что дает возможность лучше увидеть различия между полными временами полета для полученных траекторий. Сами оптимальные траектории представлены на рис. 3. Чем выше находится график на рис. 3, тем меньше время движения самолета для соответствующего решения.

Рассмотрим указанные траектории с точки зрения приведенных выше ограничений. Выберем на каждой из них множество точек, удаленных друг от друга на 5000 м. Для каждой из выбранных точек вычислим значение функции ${{C}_{y}} = {{C}_{y}}(M)$ и отобразим его на рис. 4, где показано множество, состоящее из всех этих значений для всех десяти кривых. Наиболее близкий к данному множеству участок кривой (3.1), ограничивающей максимально возможное значение коэффициента подъемной силы ${{C}_{{ymax}}}$, показан на рис. 4, кривая 1. Видно, что все точки лежат значительно ниже данной кривой и, следовательно, все полученные решения удовлетворяют условию (3.1).

Теперь для той же самой совокупности удаленных друг от друга на 5000 м точек всех десяти траекторий возьмем значения модуля скорости $V$ и высоты $y$ и поместим соответствующие точки на координатную плоскость с абсциссой $V$ и ординатой $y$. Две кривые на рис. 5 показывают внешние границы допустимого множества согласно формулам (3.2), в то время как высота полета $y$ должна удовлетворять неравенству 100 м $ \leqslant y \leqslant $ 14 000 м.

Из рис. 5 следует, что ограничения на модуль скорости и высоту нигде не достигаются, хотя несколько точек, соответствующих дозвуковому полету на низких высотах, находятся вблизи ограничения.

Ограничения (3.3) на величину силы тяги двигателей достигаются, часть каждого из оптимальных решений лежит на них. В качестве примеров на рис. 6–8 показаны эти ограничения и найденные оптимальные зависимости силы тяги от дальности. Данные рисунки соответствуют самой верхней, пятой и самой нижней кривым на рис. 3. На рис. 6–8 обозначено 1 – найденное решение, 2 – график ${{P}_{{max}}} = {{P}_{{max}}}(x)$, 3 – график ${{P}_{{min}}} = {{P}_{{min}}}(x)$. Видно, что решение выходит как на верхнее, так и на нижнее ограничение, причем на отрезках, которые нельзя считать малыми по сравнению с общей дальностью полета. Заметим, что локальные максимумы и минимумы на графиках, обозначенных 1, вызваны вычислительными погрешностями. Что касается ограничений на величины угла наклона траектории $\theta $ и нормальной скоростной перегрузки ${{n}_{y}}$, то для рассматриваемых оптимальных траекторий справедливы неравенства $ - 11^\circ < \theta < 6^\circ $ и $0.8 < {{n}_{y}} < 1.3$. Следовательно, указанные ограничения удовлетворены, причем выход на них отсутствует и полученные значения далеки от предельно допустимых.

Заключение. Из всех многочисленных и сложных ограничений, наложенных на оптимальные решения, существенными оказались лишь условия на высоту полета и силу тяги двигателей. Полученный результат является новым и ранее в научной литературе не упоминался. Из рис. 3 следует, что наибольшее влияние на вид оптимальной траектории оказало ограничение максимальной высоты полета. Что касается условий, наложенных на силу тяги двигателей, то можно использовать следующий стандартный способ (см., например, [17]). Введем скаляр $u$ согласно следующим формулам:

Это преобразование позволяет заменить управление $P$ с ограничениями (3.3), зависящими от фазовых переменных, на управление $u$, ограничения на которое не зависят от фазовых переменных.

Таким образом, при вычислениях согласно принципу максимума Л.С. Понтрягина достаточно принимать во внимание лишь одно фазовое ограничение – на максимальную высоту полета. Затем необходимо убедиться в том, что все остальные условия выполнены. И только в случае, когда хотя бы одно из них нарушено, целесообразно применять методы, требующие больших вычислительных ресурсов.

Работа частично поддержана РФФИ (грант 21-51-12004) и средствами государственного бюджета по теме государственного задания (госрегистрация 123021700055-6).