Известия РАН. Теория и системы управления, 2022, № 5, стр. 131-151

Введение. В настоящее время разработка современных систем самолетовождения (ССВ) является одним из важнейших направлений в создании бортового оборудования самолетов гражданской авиации. К числу основных задач ССВ относятся:

1) формирование горизонтального маршрута полета, проходящего через заданные точки на поверхности Земли от аэродрома вылета до аэродрома прилета;

2) формирование высотно-скоростного профиля полета воздушного судна (ВС) в соответствии с заданными авиакомпанией (АК) требованиями по эффективности использования его в данном конкретном полете;

3) формирование в бортовую систему управления (БСУ) управляющих сигналов и команд, обеспечивающих реализацию рассчитанных по п.п. 1, 2 горизонтального и высотно-скоростного профилей полета;

4) прогноз требуемых значений отдельных параметров полета в заданных точках маршрута (например, так называемая задача 4D навигации) с расчетом времени прибытия в конечный пункт маршрута и затрат топлива на весь полет;

5) формирование необходимых “подсказок” и сигналов экипажу, помогающих выполнить предписанные ему действия.

Наибольшие сложности в научном плане связаны с решением второй и четвертой задач. Именно они определяют основное научное содержание задачи вертикальной навигации и позволяют отнести ее к классу задач управления и прогноза, а не к навигации в обычном понимании этого термина. Решение этих задач должно проводиться на борту ВС непрерывно в процессе полета по сформированному в п. 1 маршруту при заданных аэродинамических характеристиках конкретного ВС и характеристиках его двигателей, а также с учетом заданного АК критерия эффективности и доступной информации о состоянии атмосферы. При этом в соответствии с требованиями международных стандартов ошибка в прогнозе времени прибытия должна быть меньше 30 с, а в отдельных случаях – 5 с.

Критерием эффективности для самолетов гражданской авиации принимается определяемый АК экономический критерий в виде стоимости конкретного полета. Он является линейной комбинацией стоимости затраченного топлива и времени полета, связанных так называемым индексом стоимости.

Вопросам оптимизации режимов полета ВС и предсказания траекторий уделялось значительное внимание в работах сотрудников различных организаций, например, [1–6]. Их общей чертой является то, что в них не в полной мере учитываются наиболее важные для задачи вертикальной навигации особенности полета ВС. В частности:

используются некоторые абстрактные характеристики ВС и двигателей;

не всегда корректно учитываются изменение массы ВС, скорости ветра на высотах полета, отличие параметров атмосферы от стандартных значений;

роль АК и системы управления воздушным движением (СУВД) при формировании критериев эффективности;

ограничения, накладываемые на режимы полета ВС руководствами по летной эксплуатации (РЛЭ).

Данная статья является продолжением работ [7, 8], где представлены результаты исследования режимов взлет и набор высоты. Здесь применительно к крейсерскому полету предлагается методика и приводятся результаты решения задачи оптимизации полета ВС с учетом существенных, но обычно не учитываемых факторов.

1. Режимы движения ВС (сценарии) на этапе крейсерского полета. Режимы движения ВС на этапе крейсерского полета задаются в РЛЭ каждого конкретного ВС. Кроме того, они определяются АК и СУВД. Чтобы в наибольшей степени охватить возможное множество возникающих при этом режимов в качестве типовых, на этапе крейсерского полета будем рассматривать следующие:

1) полет с оптимальным (либо с заданным М_зад) значением числа Маха М_opt = const на заданной высоте H_зад;

2) полет на оптимальной (либо заданной H_зад) высоте H_opt = const с заданным числом М_зад;

3) полет с оптимальным значением числа М_opt = const на оптимальной высоте H_opt = const.

4) полет с переходом с текущей высоты H на оптимальную высоту H_opt (либо переход с текущего значения числа М на М_opt);

5) полет с переходом на очередной оптимальный эшелон;

6) оптимальный по выбранному критерию полет с обеспечением заданного времени прибытия T_зад в очередную точку маршрута.

Здесь заданные значения высоты и числа Маха формируются экипажем ВС в соответствии с руководящими документами и указаниями АК и СУВД.

Оптимальные значения этих параметров должны рассчитываться на борту ВС (в процессе предполетной подготовки и в полете) по заданным АК критериям оптимальности в зависимости от:

текущих параметров полета;

доступной информации о состоянии атмосферы;

имеющихся в бортовой базе данных по актуальным аэродинамическим характеристикам ВС и его двигателей.

2. Критерии оптимизации. В настоящее время в качестве основного критерия при формировании вертикального профиля полета ВС принимаются прямые эксплуатационные расходы (C), которые складываются из стоимости затраченного топлива и времени полета. АК заинтересована в минимизации C и поэтому требует, чтобы ССВ обеспечила в каждом конкретном полете реализацию траектории с минимизацией прямых эксплуатационных расходов. При этом в стоимости времени учитывается множество факторов, к которым относятся: стоимость погрузоразгрузочных работ, стоимость аэронавигационного сопровождения полета ВС, аренда и амортизация ВС, различные регламентные работы, зарплата и страховка экипажа и др. Очевидно, что критерий C для этапа крейсерского полета можно записать в виде

где q_т – расход топлива в единицу времени, C_т, C_t – стоимость единицы топлива и времени соответственно, t_k – длительность участка крейсерского полета.

Учитывая, что стоимость топлива C_т для конкретного маршрута полета является величиной постоянной, критерий стоимости полета можно заменить критерием относительной стоимости I_C , т.е. отношением стоимости полета C к величине стоимости единицы топлива C_т:

где $CI = {{C}_{t}}{\text{/}}{{C}_{{\text{т}}}}$ и называется индексом стоимости (CostIndex).

Заметим, что критерий I_C имеет размерность расхода топлива (кг, т, …), поэтому его можно рассматривать как приведенный расход топлива.

Наряду с критерием относительной стоимости или приведенного топлива (2.1) будем считать значение фактического расхода топлива I_т и полное время полета I_t = t_k, как имеющее в ряде случаев важное самостоятельное значение.

3. Модель движения ВС. При решении любой задачи, связанной с исследованием и оптимизацией движения какого-либо объекта, построение модели его движения – одна из важнейших составляющих методики решения поставленной задачи. Здесь должна быть найдена та золотая середина, которая бы охватывала все существенные особенности рассматриваемого явления и, в то же время, позволяла бы довести решение до конца с использованием предполагаемых методов.

К числу основных особенностей, характерных для этапа крейсерского полета, относятся:

доходящая до нескольких часов длительность и большая протяженность этапа;

наличие на высотах полета значительной горизонтальной составляющей скорости ветра, достигающей сотен километров в час;

существенное изменение в процессе полета массы ВС;

движение ВС в поле тяготения Земли.

Анализ имеющихся в научно-технической литературе моделей движения ВС показывает, что они в полной мере не учитывают эти особенности. Это связано с тем, что основное внимание уделяется либо вопросам динамики полета, либо оценкам потенциальных возможностей и разработке методов оптимизации, в которых неизбежно (в угоду предполагаемым методам исследования) проводятся существенные упрощения модели. Длительность и протяженность полета предполагают использование полных моделей движения ВС с учетом формы Земли и ее вращения, имеющихся, например, в публикациях [9–11].

Значительная продолжительность этапа полета и сравнительно небольшие длительности переходных процессов при отработке углов ориентации ВС позволяют в описании модели движения ВС ограничиться только движением центра масс. Так поступают, например, в задачах управления траекториями аэрокосмических аппаратов [10, 11].

В полной модели движения твердого тела переменной массы, каковым является ВС, относительно Земли помимо активных сил и сил гравитации должны участвовать переносное и кориолисово ускорения, ускорение от изменения массы и систематических составляющих скорости ветра. Как показывают расчеты, дополнительные ускорения от этих факторов могут достигать значений (2–4) × 10^–3 g. Понятно, что при формировании законов управления такими величинами, по сравнению с ускорениями от активных сил, можно пренебречь. Однако в задачах прогноза на длительных интервалах времени, что имеет место в задачах вертикальной навигации, не учет их приводит к ошибкам в несколько десятков километров за час полета или несколько минут в оценке времени прибытия в заданную точку маршрута.

По этим причинам модели для решения задач управления (п. 2) и прогноза (п. 4) должны быть различны. В качестве основы для разработки “прогнозных” моделей могут использоваться соотношения (1.16) публикации [10] после коррекции их в направлении более строгого учета изменения массы ВС и скорости ветра.

3.1. Уравнения движения центра масс ВС. При разработке модели движения ВС для решения задач формирования вертикального профиля полета ввиду достаточной малости значений переносного и кориолисового ускорения можно отказаться от их учета. В этом случае движение ВС рассматривается как движение его центра масс в равномерном гравитационном поле относительно плоской поверхности Земли. Соответствующая этой модели система дифференциальных уравнений, описывающая движение центра масс, получится из теоремы об изменении количества движений материальной точки.

Как показано в [7, 8], при принятых допущениях уравнения движения центра масс ВС при формировании профиля полета имеют вид

где ${{\alpha }_{{\text{т}}}} = \alpha - \varphi $; ${{\alpha }_{{\text{т}}}}$ – угол между вектором тяги Т и вектором воздушной скорости V; $\alpha ,\varphi $ – угол атаки и угол установки двигателя относительно связанной оси ВС; $V$ – воздушная скорость; $\theta $ – угол наклона вектора воздушной скорости к горизонту; $H$ – высота полета; $L$ – пройденная дальность; $m$ – масса ВС; $T$ – сила тяги двигателей ВС; $Y,R$ – аэродинамические подъемная сила и сила сопротивления, соответственно; U_в – горизонтальная составляющая скорости ветра; q_c – секундный расход топлива; g – ускорение силы тяжести.

Входящие в систему уравнений (3.1) силы R, Y определяются соотношениями:

где ${{C}_{X}}$, ${{C}_{Y}}$ – коэффициенты соответствующих аэродинамических сил; $S$ – характерная площадь крыла ВС; $\rho $ – плотности воздуха на высоте H.

Аэродинамические коэффициенты С_X, С_Y в общем случае являются функциями большого числа переменных, однако в траекторных задачах принято использовать их балансировочные значения в виде функций

где Fl – параметр, определяющий конфигурацию крыла.

Входящие в уравнения (3.1) сила тяги Т и секундный расход топлива q_c обычно задаются таблицами или графиками вида

где $M$ – число Маха; δ – эквивалент тяги, в качестве которого могут использоваться обороты турбины или компрессора высокого или низкого давлений, положение рычага управления двигателем и др.; $\Delta {{T}^{0}}$ – отклонение температуры наружного воздуха от стандартного для данной барометрической высоты Н_б.

Из анализа соотношений (3.1)–(3.4) видно, что при известных (таблицы, графики, аналитические выражения и т.д.) зависимостях ${{C}_{X}} = {{C}_{X}}(\alpha ,M,Fl)$, ${{C}_{Y}} = {{C}_{Y}}(\alpha ,M,Fl)$, $T = T(\delta ,{{H}_{{\text{б}}}},M,\Delta {{T}^{0}})$, ${{q}_{{\text{c}}}} = {{q}_{{\text{c}}}}(\delta ,{{H}_{{\text{б}}}},M,\Delta {{T}^{0}})$ система уравнений (3.1) однозначно определяет динамику изменения фазовых переменных H, L, V, θ, m при заданных $\alpha (\xi )$, $\delta (\xi )$, где $\xi $ – какой-либо монотонный скалярный параметр, например время t или дальность L.

Функции

являются в данном случае для системы (3.1) функциями управления или управляющими функциями, на которые наложены ограничения:

Соответствующие ограничения могут быть наложены и на фазовые координаты модели.

Как видно из соотношений (3.1)–(3.6), данная модель отличается от большинства известных тем, что в ней корректно учитываются изменение в процессе полета массы ВС и влияние горизонтальной составляющей скорости ветра.

3.2. Модель активных сил. Значения аэродинамических коэффициентов, а также высотно-скоростные и дроссельные характеристики двигателей ВС обычно задаются в виде таблиц на сетке их аргументов. Такой способ задания не позволяет получать с достаточной точностью производные от этих величин, что обычно требуется при использовании методов оптимизации и построении вычислительных алгоритмов.

Как показали исследования, эффективным методом преобразования заданных табличных значений в соотношения с требуемыми свойствами непрерывности является метод последовательной полиномиальной аппроксимации, суть которого на примере С_X, С_Y понятна из рассмотрения соотношения

где B(I, J) – матрица постоянных коэффициентов для приведенного диапазона переменных.

Как можно видеть, в этом методе в зависимости от предполагаемого в дальнейшем применения, выбирается последовательность аргументов, по которым проводятся последовательная полиномиальная аппроксимация исходных таблиц и далее – коэффициентов предыдущей аппроксимации.

Аналогичные выражения получаются при аппроксимации значений тяги двигателей Т и секундного расхода топлива q_c по четырем переменным δ, M, H_б, ΔT⁰.

Такая аппроксимация обеспечивает для системы (3.1) непрерывность правых частей и их производных не ниже второго порядка, что достаточно для применения современных методов оптимизации и вычислительной математики.

3.3. Модель атмосферы. Особое место в модели движения ВС имеет модель атмосферы. Это связано с тем, что полет ВС практически всегда проходит в условиях отличных от стандартных (принятых в ГОСТ 4401-81), а задачи вертикальной навигации должны решаться для текущих условий полета. По этой причине необходимо иметь модель атмосферы, учитывающую отличия реальной атмосферы от стандартной, в частности по давлению на уровне моря и температуре воздуха.

Полагая, что атмосфера Земли в целом находится в равновесии, для построения модели воспользуемся основным уравнением статики атмосферы (3.8), которое справедливо для любого газа. Оно выражает условие равновесия между двумя силами, действующими на единицу массы воздуха по вертикали: силой вертикального барического градиента и силой тяжести:

где $\rho ,\;p$ – массовая плотность и давление воздуха; g – ускорение силы тяжести; h – геометрическая высота относительно среднего уровня моря.

Дополнив (3.8) соотношением Менделеева–Клайперона, определяющим связь между температурой Т, давлением p и плотностью $\rho $ воздуха

где ${{R}_{{{\kern 1pt} {\text{c}}}}}$ = 287.05287 – удельная газовая постоянная сухого воздуха, из (3.8), (3.9) получим дифференциальное уравнение для определения давления p как функции высоты h:

Отсюда следует, что зависимость давления сухого воздуха P от высоты h в предположении о постоянстве g определяется только видом зависимости температуры воздуха от высоты, т.е. $p(h) = F(T(h))$.

Задаваясь различными законами изменения температуры Т по высоте h, из (3.9), (3.10) можно получить различные модели атмосферы. Для задач вертикальной навигации наибольший интерес представляют следующие две модели атмосферы:

стандартная атмосфера по ГОСТ 4401-81;

так называемая квазистандартная атмосфера, учитывающая возможные отклонения ΔT температуры и давления ΔP от стандартных значений.

В первом случае функция F имеет вид

где a = 6.5 град/км, ${{T}_{0}}$ = 288.15 К, T₁₁ – значение температуры на h = 11 км.

Во втором случае функция F сохраняет вид (3.11), но могут меняться значения параметров ${{T}_{0}}$ и a , а также значение давления на уровне моря ${{P}_{0}}$.

При таком способе построения функции F решение уравнения (3.10) имеет единый вид для обеих моделей атмосферы:

где ΔP, ΔT – отклонение давления и температуры от стандартных значений.

Выражения (3.12) однозначно определяют значения давления воздуха на высоте h как в случае стандартной, так и квазистандартной моделях атмосферы. Далее, используя уравнение Менделеева–Клайперона, получаем соотношение для нахождения плотности воздуха на высоте h:

Выражения (3.12), (3.13) полностью определяют обе требуемые модели атмосферы и показывают влияние температуры воздуха на ее параметры.

4. Постановка задачи. Принципиальное влияние на постановку задачи оказывает принятый в гражданской и транспортной авиации подход к формированию высотно-скоростного профиля полета ВС. Суть его заключается в том, что для АК и СУВД важно знать (и задавать) текущие значения параметров режима полета ВС, какими могут служить высота, скорость, угол наклона траектории, вертикальная скорость и др. Поэтому целью данной работы прежде всего является исследование влияния этих параметров на значения критериев качества и нахождение их оптимальных значений, а не поиск оптимальных законов управления при перелете ВС из точки А в точку В. Целесообразность такого подхода связана с тем, что в современных ВС гражданской и транспортной авиации собственно управлением (созданием целенаправленных воздействий на органы управления ВС) занимается бортовая система управления, а ССВ задает ей необходимые значения параметров режима полета (уставки), обеспечивающие выполнение целей управления.

Множество указанных выше конкретных сценариев этапа крейсерского полета приводит к множеству постановок задач для каждого из возможных сценариев, которые рассматриваются далее.

5. Вектор параметров модели. Анализируя уравнения модели (3.1)–(3.6) и определение критерия оптимизации (2.1), можно видеть, что в данном случае значение критерия зависит от следующего набора параметров:

параметры режима полета: М_зад, Н_зад;

параметр, задаваемый авиакомпанией: CostIndex (CI);

вес ВС в начале этапа полета (m);

дальность полета L;

составляющая ветра на направление полета (U_в);

отклонение температуры воздуха на барометрической высоте от стандартного значения (ΔT⁰).

Заметим, что здесь, так же как и в (3.1), под U_в понимается не проекция ветра на направление полета, а составляющая ветра, получаемая из рассмотрения навигационного треугольника.

6. Оптимизация сценариев 1–3. В этих сценариях требуется найти оптимальные по критерию (2.1) значения какого-то одного (сценарии 1, 2) или обоих (сценарий 3) параметров режима полета ВС при заданных значениях других параметрах. Так как по условиям реализации сценариев предполагается, что как оптимизируемые, так и заданные параметры постоянны на каждом отдельном этапе полета, то задача сводится к процедуре поиска минимальных значений некоторой функции

где F – значение оптимизируемого критерия (I_т, I_С); p_упр, p_зад – вектор управляющих и заданных параметров соответственно.

Дополнив (6.1) условиями

приходим к классической задаче поиска минимального значения функции (6.1) при ограничениях на область определения переменных в виде (6.2). Особенность данной задачи состоит в том, что ее решение должно проводиться в два этапа.

На первом из них для заданного вектора ${{p}_{{{\text{зад}}}}}$ необходимо найти законы формирования управляющих функций (3.5), обеспечивающие требуемые свойства вектору ${{p}_{{{\text{упр}}}}}$, и провести интегрирование системы уравнений (3.1) с учетом соотношений (3.2)–(3.6).

В сценариях 1, 2 вектор ${{p}_{{{\text{упр}}}}}$ состоит из одной компоненты: М_зад в сценарии 1 и Н_зад в сценарии 2. В сценарии 3 вектор p_упр имеет две компоненты М_зад и Н_зад. Остальные параметры модели образуют вектор p_зад. Эту задачу будем называть задачей управления. Решение ее позволяет получить значение функции F, соответствующее выбранной паре векторов ${{p}_{{{\text{упр}}}}}$ и ${{p}_{{{\text{зад}}}}}$.

На втором этапе проводится собственно оптимизация функции F по вектору ${{p}_{{{\text{упр}}}}}$, его будем называть этапом оптимизации.

Каждый из этапов имеет различные постановки задач и методы решения, поэтому рассмотрим их отдельно.

6.1. Этап управления. Как можно понять из описания сценариев 1–3, для них задача управлений формулируется одинаково. Для системы (3.1) при заданном векторе ${{p}_{{{\text{зад}}}}}$ требуется найти такие законы формирования управляющих функций $\alpha = \alpha (\xi )$, $\delta = \delta (\xi )$, при которых после окончания переходного процесса выполняются условия

Интегрирование системы (3.1) здесь ведется из начального вектора фазовых координат до заданного значения конечной дальности L. Заметим, что задача интегрирования системы (3.1) формально поставлена корректно в том смысле, что имеются две управляющие функции (3.5), обеспечивающие требуемое поведение двух параметров режима полета (6.3).

6.1.1. Возможные подходы к получению уравнений управляемой системы. При формировании законов управления (3.5) системы (3.1) можно использовать различные подходы, к числу которых относятся методы построения пропорционально-интегрально-дифференцирующих (ПИД) регуляторов, методы концепции обратных задач динамики и методы теории оптимального управления.

В частности, в первой группе этих методов можно положить

где p – оператор дифференцирования; ${{\tau }_{\alpha }},{{\tau }_{\delta }}$ – постоянные времени, зависящие от динамических свойств ВС и его двигателей; K₁, K₂, K₃, K₄, K₅, K₆ – коэффициенты, определяемые стандартными методами.

При использовании ПИД-регуляторов система уравнений модели (3.1) дополняется двумя уравнениями (6.4), динамика которых определяется значениями постоянных времени ${{\tau }_{\alpha }},{{\tau }_{\delta }}$, оказывающих существенное влияние на допустимый шаг интегрирования расширенной системы уравнений.

Как показали результаты моделирования, для каждого режима полета можно подобрать параметры законов управления (6.4), обеспечивающие в рассматриваемых диапазонах изменения высоты и скорости полета ВС приемлемое качество переходных процессов. Основным недостатком метода ПИД-регуляторов оказалось значительное ограничение допустимого шага интегрирования (единицы секунд) расширенной системы по сравнению с исходной, что приводит к существенному увеличению времени расчета траектории полета.

В методе концепции обратных задач динамики необходимо задаться в системе (3.1) желаемым законом изменения фазовых переменных $V$ и $\theta $, например

где ${{\tau }_{V}},\;{{\tau }_{\theta }}$ – постоянные времени, определяющие динамику изменения переменных V и $\theta $. Значения ${{V}_{{{\text{зад}}}}}$ и ${{\theta }_{{{\text{зад}}}}}$ находим из условия (6.3). Откуда следует, что ${{V}_{{{\text{зад}}}}} = {{M}_{{{\text{зад}}}}}{{а}_{{{\text{зв}}}}}$, ${{\theta }_{{{\text{зад}}}}} = 0$, где ${{a}_{{{\text{зв}}}}}$ – скорость звука на заданной высоте.

С учетом этого система (3.1) принимает вид

Необходимые для интегрирования этой системы значения отклонения рукоятки управления двигателем $\delta = \delta (\xi )$ (а, следовательно, значения расхода топлива q_c) находим из решения системы уравнений:

Сравнивая (3.1) и (6.6), можно видеть, что использование концепции обратных задач динамики обеспечивает существенное упрощение интегрируемых уравнений. Кроме того, так как $({{\tau }_{V}},{{\tau }_{\theta }}) > ({{\tau }_{\alpha }},{{\tau }_{\delta }})$, полученные уравнения, по сравнению с ПИД-регулятором, допускают применение большего шага интегрирования. Однако за эти положительные свойства приходится платить необходимостью решения уравнений (6.7) на каждом шаге интегрирования системы (6.6).

При решении  поставленных задач методом теории оптимального управления Понтрягина Л.С. [12] управляемая система описывается уравнениями (3.1) без какого-либо расширения, а условия (6.3) выступают как ограничения на фазовый вектор. Известные сложности решения таких задач даже для одного варианта условий (6.3) делают бесперспективным применение этого метода на втором этапе решения задачи при поиске оптимальных значений M_зад и H_зад.

Анализ рассмотренных подходов к построению модели управляемой системы показывает, что способ, основанный на использовании метода обратных задач динамики, более предпочтителен по сравнению с двумя другими. Это связано с тем, что здесь в рассматриваемой системе (6.6) можно исключить первые два уравнения и решать задачу для установившихся значений V и $\theta $. Возможность такого подхода базируется на результатах работы [13], где показано, что учет переходных процессов в рассматриваемых задачах практически не влияет на критерии качества управления.

С помощью данного предположения все составляющие фазового вектора, кроме m, становятся известными, а в системе (6.6) остаются только два последние уравнения, которые должны интегрироваться с учетом выполнения двух соотношений:

где T, q_c – функции $\delta $; α_т, R, Y – функции α.

Таким образом, задача управления как определение управляющих функций системы (3.1), обеспечивающих выполнение условий (6.3), свелась к решению уравнений (6.8) относительно $\delta = \delta ( \cdot )$ и $\alpha = \alpha ( \cdot )$.

6.1.2. Способы решения системы уравнений (6.8). Как известно, общих методов решения систем подобного типа не существует и успех решения зависит от наличия в них тех или иных особенностей. В данном случае такой особенностью является линейность системы (6.8) относительно переменных T, q_c. Это позволяет получить следующие решения $T{\kern 1pt} {\text{*}},q{\kern 1pt} *$:

где

В соответствии с принятой ранее формой аппроксимации активных сил запишем

Отсюда следует, что левые части соотношений (6.9) являются полиномами третьей степени относительно δ, а правые части – функциями α. Поэтому (6.9) представим как

где

Решение системы алгебраических уравнений (6.10) проведем с использованием теории полиномов [14, 15]. Первым шагом этого метода является построение так называемого результанта полиномов $R(f,g)$:

Вычеркиванием из (6.11) первых и последних строк и столбцов получаем первый субрезультант ${{R}_{1}}(f,g)$ полиномов $f(x),\;g(x)$:

Вычеркиванием из (6.11) первой и последней строки, а также первого и предпоследнего столбцов находим $\det {{M}_{1}}$:

Далее пользуемся теоремой о существовании решения системы алгебраических уравнений, в соответствии с которой для того, чтобы полиномы $f(x),\;g(x)$ имели только один общий корень x* необходимо и достаточно, чтобы $R(f,g) = 0,$ а ${{R}_{{{\kern 1pt} 1}}}(f,g) \ne 0$. При этом

Наличие данной теоремы обуславливает целесообразность использования для решения системы (6.10) рассмотренного способа, хотя возможны и другие решения.

Конкретное применение предлагаемого метода сводится к получению выражений для $R(f,g)$, ${{R}_{1}}(f,g)$ и $\det {{M}_{1}}$ как функций коэффициентов ${{a}_{i}},{{b}_{i}},i = 0,3$. Проведя необходимые преобразования, находим

где

Первые два выражения (6.13) позволяют вычислить $x{\kern 1pt} *$ по соотношению (6.12) и, таким образом, найти значение отклонения рукоятки управления двигателем $\delta {\kern 1pt} *$, обеспечивающее (совместно c углом атаки $\alpha {\kern 1pt} *$) выдерживание заданных значений высоты и скорости. Значение угла атаки $\alpha {\kern 1pt} *$ находим из условия $R(f,g) = 0$. Если имеется несколько решений этого уравнения, то из них выбирается то, которое удовлетворяет заданным ограничениям (3.6).

Как показали результаты моделирований, предлагаемый метод дает возможность существенно увеличить допустимый шаг интегрирования модели движения. Так, если в основу модели принять соотношения (6.6), (6.7), то допустимый шаг интегрирования не превышает 2 с, а предлагаемый метод позволяет увеличить шаг интегрирования до 500 с и более. Одной из особенностей его является необходимость решения уравнения $R(f,g) = 0$ на каждом шаге интегрирования, что предполагает использование методов численного решения подобного рода уравнений. Так как эти методы основаны на выполнении тех или иных итерационных процедур, то данный метод будем называть итерационным.

Помимо рассмотренного, существует и другой метод решения уравнений (6.8), в котором практически отсутствует необходимость численного решения уравнений типа $R(f,g) = 0$.

Суть его заключается в переходе от конечных уравнений (6.8) к их дифференциальным аналогам, поэтому он называется дифференциальным. Для этого заметим, что в (6.8) участвует параметр m (масса ВС), следовательно решения $\alpha {\kern 1pt} *$ и $\delta {\kern 1pt} *$ являются функциями m. Так как уравнения (6.8) должны выполняться при любых значениях m, то производная от левых частей этих уравнений по m равна нулю, т.е.

где ${{F}_{{i\alpha }}},\;{{F}_{{i\delta }}},\;{{F}_{{im}}},$ $i = 1,2$, – частные производные ${{F}_{1}},\;{{F}_{2}}$ по $\alpha ,\;\delta ,\;m$ соответственно. Из (6.14) находим где

Проведя необходимые преобразования, получим

где ${{T}_{\delta }},\;{{q}_{\delta }}$ – производные от силы тяги и расхода топлива по $\delta $; ${{R}_{\alpha }},\;{{Y}_{\alpha }}$ – производные от аэродинамических сил по углу атаки α.

Уравнения (6.15) являются дифференциальными аналогами уравнений (6.8), они определяют искомые управляющие функции $\delta = \delta ( \cdot )$ и $\alpha = \alpha ( \cdot )$ как решения системы уравнений:

Эти уравнения описывают движение центра масс ВС при полете его на заданной высоте $H = {{H}_{{{\text{зад}}}}}$ = const со скоростью $M = {{M}_{{{\text{зад}}}}} = {\text{const}}$. В таком режиме полета при заданной дальности неизвестными являются время и расход топлива, которые и определяются в процессе решения системы уравнений (6.16).

Как показывают результаты моделирования, оба предлагаемых способа достаточно близки по затратам времени моделирования на движения центра масс ВС, причем первый из них дает завышенное (на несколько килограмм) значение расхода топлива, а второй – такого же порядка заниженное значение.

6.2. Этап оптимизации. Использование итерационного или дифференциального методов решения соотношений (6.8) обеспечивает построение модели объекта управления и позволяет провести моделирование режима полета ВС при любых заданных значениях высоты и скорости, а следовательно, и рассчитать значения критериев стоимости ${{I}_{C}}$, расхода топлива I_т и времени полета ${{I}_{t}}$ на дальность L. При заданном значении CI эти критерии, как следует из (2.1), связаны соотношением ${{I}_{C}} = {{I}_{{\text{т}}}} + {{I}_{t}}CI$. Поэтому фактически для расчета всех критериев достаточно найти только значение I_т, что можно выполнить, например, путем интегрирования системы уравнений (6.16) до заданной дальности L_зад. Заметим, что начальные условия $\alpha (0)$ и $\delta (0)$ для системы (6.16) всегда получаются с помощью итерационного метода решения соотношений (6.8).

Как можно видеть, в сценариях 1, 2 задача оптимизации сводится к поиску оптимального значения управляющего параметра M или H при заданных значениях всех остальных параметров сценария. В сценарии 3 оптимизацию критериев необходимо провести по двум параметрам. Таким образом, задача оптимизации параметров в сценариях 1–3 свелась к классической задаче нелинейного программирования, для которой имеется множество эффективных методов решения. Выбор конкретного из них определяется целью проводимого исследования и практическим использованием результатов. В данном случае работа направлена на создание методической основы для разработки алгоритмического обеспечения бортового вычислителя ССВ. Это предъявляет ряд требований к вычислительному методу, к числу которых относятся: предсказуемость времени выполнения расчетов, глобальность полученного решения и простота контроля процесса вычислений. С другой стороны, требования к методу диктуются свойствами оптимизируемых функций, такими, как непрерывность, дифференцируемость и количество локальных оптимумов.

Проведенное исследование показало, что в рассматриваемом диапазоне высот и скоростей полета среднемагистральных ВС типа SSJ-100 типичный вид зависимости ${{I}_{{\text{т}}}} = {{I}_{{\text{т}}}}(M{\text{/}}m)$ для значений массы m = 33–44 т имеет вид, представленный на рис. 1.

Показанный на рис. 1 характер зависимости ${{I}_{{\text{т}}}} = {{I}_{{\text{т}}}}(M{\text{/}}m)$ сохраняется для других значений параметров сценария, что позволяет сделать вывод о непрерывности, выпуклости и унимодальности критерия расхода топлива ${{I}_{{\text{т}}}}$, а следовательно, и критерия стоимости ${{I}_{C}}$, во всей области определения параметров. Кроме того, на рис. 1 просматривается одна специфическая особенность зависимости ${{I}_{{\text{т}}}} = {{I}_{{\text{т}}}}(M{\text{/}}m)$ от числа М. Она проявляется в том, что в диапазонах чисел М = = 0.778–0.799 ее значения изменяются практически по линейному закону. Это обстоятельство оказывает существенное влияние на выбор метода поиска оптимальных значений М_opt и H_opt и характер их зависимостей от параметров модели.

Учитывая сформулированные выше требования к методу оптимизации и показанный на рис. 1 характер зависимости критериев ${{I}_{{\text{т}}}}$ и ${{I}_{C}}$ от числа М в статье в качестве рабочего используется сочетание прямого метода поиска минимума методом перебора с первой итерацией метода квадратичной параболы. Главными достоинствами такого способа является возможность получения глобального минимума и постоянное количество расчетов оптимизируемой функции (что существенно при организации вычислительного процесса в бортовых вычислителях). Единственным параметром такого способа поиска минимума выступает количество интервалов N, на которое разбивается область оптимизируемых переменных (6.4).

Будем считать, что ошибка в определении М_opt не должна превосходить нескольких тысячных долей единицы, а ошибка в определении H_opt допускается в пределах сотни метров. Учитывая, что ошибка в методе перебора задается соотношением $\varepsilon \leqslant (b - a){\text{/}}N$, где a, b – границы области определения, а ошибка первой итерации метода квадратичной параболы пропорциональна ${{\varepsilon }^{3}}$, в данном случае с достаточной степенью надежности можно принять N = 10.

Таким образом, все необходимые модели и методы оптимизации подготовлены для практического решения поставленных задач. Остался последний вопрос: нужно ли этим вообще заниматься и что с этого можно получить? Для ответа на этот вопрос рассмотрим показанные на рис. 2, 3 результаты расчетов значений критериев стоимости ${{I}_{С}}$ как функции числа Маха и высоты для различных значений массы m.

Результаты обработки представленных на рис. 1–3 графиков показаны в табл. 1, 2, где первая относится к рис. 2, а вторая – к рис. 3.

В этих таблицах для каждого значения массы ВС приведены минимальные (Min), максимальные (Max) и среднеарифметические (Cр.ар.) значения критериев (I_C, ${{I}_{{\text{т}}}}$, I_t), а также величины разброса значений в абсолютном (Дельта) и процентном (%) измерениях. Как можно видеть из табл. 1, экономия топлива за счет оптимального выбора числа М может достигать 15–20%, сокращение времени – порядка 9% и выигрыш по стоимости составляет 4.5–5% (последние цифры соответствуют значению CI = 0.5). За счет оптимального выбора высоты экономия топлива, как видно из табл. 2, может составлять 10–30%, времени – порядка 3% и стоимости – порядка 3–10%. Эти результаты говорят о практической значимости рассматриваемых далее задач оптимизации сценариев полета ВС.

6.2.1. Исследование и оптимизация сценария 1. В этом сценарии требуется найти оптимальное значение числа Маха при заданных значениях всех других перечисленных в разд. 5 параметрах. Для этого интервал возможных значений числа М разбивается на 10 равных частей и для каждой точки ${{M}_{i}}$ проводится интегрирование системы уравнений (6.16). При этом все остальные значения параметров модели остаются постоянными. По полученным в процессе интегрирования на заданную дальность L значениям времени и расходу топлива формируются массивы значений ${{I}_{t}}(i) = {{I}_{t}}({{M}_{i}}{\text{/}}p)$ и ${{I}_{т}}(i) = {{I}_{т}}({{M}_{i}}{\text{/}}p)$, $i = \overline {1,11} $, где p – множество других (кроме ${{М}_{i}}$) параметров модели. Далее по заданному значению $CI$ рассчитывается массив ${{I}_{C}}(i) = {{I}_{C}}({{M}_{i}}{\text{/}}p)$. В соответствии с предлагаемым выше способом оптимизации на множестве 11 значений ${{I}_{C}}(i)$ находится минимальное, и затем методом квадратичной параболы с учетом ограничений рассчитывается М_opt. Очевидно, что при этом получаем

Первые четыре параметра будем называть параметрами сценария 1, а два других – возмущающими параметрами.

Как можно понять, среди всего множества параметров особое место занимает дальность полета L, так как она определяет длительность интегрирования системы уравнений (6.16) и получающиеся значения критериев I_C, ${{I}_{{\text{т}}}}$. Расчеты показали, что при $L < 1000{\kern 1pt} - {\kern 1pt} 1200$ км значения М_opt c точностью порядка 0.001 практически не зависят от L. Поэтому в множестве параметров сценария 1 остались три параметра H, m, CI, влияние которых на значения М_opt необходимо исследовать. Разработанный алгоритм поиска М_opt позволил решить задачу исследования влияния на значение М_opt всей совокупности его параметров.

Результаты расчетов зависимости ${{M}_{{opt}}} = {{M}_{{opt}}}(m,CI{\text{/}}H)$ для H = 9; 10; 11 км показаны на рис. 4, а зависимость ${{M}_{{opt}}} = {{M}_{{opt}}}(H,CI{\text{/}}m)$ для m = 44; 38; 32 т – на рис. 5.

Представленные на рис. 4, 5 графики зависимостей М_opt от параметров сценария m, H,CI характеризуются наличием в них “изломов”. Согласно исследованиям, эти изломы являются следствием дискретного отображения соответствующих зависимостей, в которых имеются разрывы первого рода. Эти разрывы связаны с существованием показанных на рис. 1 (в диапазоне изменения числа М = 0.778–0.799) площадок с практически “постоянным” углом наклона. Заметим, что разрывы в зависимостях М_opt  от параметров сценария приводят к соответствующим разрывам в значениях времени полета ${{I}_{t}}$ и затрат топлива ${{I}_{{\text{т}}}}$, в то время как критерий стоимости ${{I}_{С}}$ сохраняет непрерывность.

Анализ приведенных результатов позволяет сделать следующие выводы:

рабочий диапазон значений CI, в котором обеспечивается выполнение заданного ограничения на число М < 0.82, ограничен значением CI < 2 кг/с;

зависимости М_opt от m и Н в диапазоне значений 0 < CI < 2 имеют разрывы первого рода, в которых происходит скачкообразное возрастание значения М_opt при увеличении m, Н или CI;

значения М_opt существенно зависят от заданного значения CI;

при CI < 1 значения М_opt возрастают с увеличением массы m и высоты Н полета ВС;

при CI > тенденция возрастания М_opt снижается и даже переходит в противоположную.

При этом экономия топлива за счет оптимального выбора числа М может достигать 15–20%, а сокращение времени – порядка 9%.

Исследование влияния на значения М_opt возмущающих факторов в виде ветра и отклонения температуры воздуха от стандартного значения показало следующее:

наличие ветра не меняет характер поведения зависимостей, а влияет только на значения М_opt;

встречный ветер по сравнению с попутным оказывает большее влияние на М_opt;

степень влияния ветра на М_opt зависит от значений CI;

при встречном ветре М_opt увеличивается, а при попутном ветре уменьшается с увеличением модуля скорости ветра;

при встречном ветре порядка 100 км/ч увеличение М_opt (в зависимости от m, H, CI) может достигать значения 0.02, а при попутном уменьшение М_opt может достигать значения –0.015;

изменение температуры воздуха в пределах ±20° практически не влияет на значение М_opt.

6.2.2. Результаты исследования и оптимизации сценариев 2–5. В сценарии 2 задача заключается в определении оптимального значения высоты Н_opt при заданных всех остальных значениях параметров полета. Параметрами сценария здесь являются M, m, CI, L, а внешними (возмущающими) служат U_в, ΔT, поэтому в общем случае можем записать

Проверка параметров сценария по степени влияния их на значения критериев оптимальности показала, что на дальностях полета 1000–1200 км с точностью $ \pm $50 м влиянием дальности можно пренебречь. Другие параметры сценария существенно влияют на критерии оптимальности. При этом типичный вид зависимости ${{H}_{{opt}}} = H(M,m{\text{/}}CI)$ для различных значений возмущающих параметров имеет вид, представленный на рис. 6.

Рассмотренные на рис. 6 зависимости получены при значении CI = 2 и стандартных значениях состояния атмосферы. Влияние CI на вид зависимости ${{H}_{{opt}}} = H(M,m{\text{/}}CI)$ прежде всего проявляется в величине скачка в точках разрыва, который примерно пропорционален значению CI и отсутствует при CI = 0. Ветер оказывает существенное влияние на вид зависимостей H_opt = = $H(M,m{\text{/}}CI)$, которое отражается в величине скачка и смещении области его расположения. При этом имеет место такая тенденция, что с увеличением (алгебраическим) скорости ветра величина скачка уменьшается и он смещается в область больших значений веса ВС.

Завершая рассмотрение сценария 2, заметим, что за счет оптимального выбора высоты экономия топлива может составлять 10–30%, а времени – порядка 3%.

В сценарии 3 задача заключается в определении оптимальных значений высоты Н_opt и числа M_opt при заданных всех остальных значениях параметров полета. Параметрами этого сценария являются $m,CI,L$, а возмущающими служат U_в и ΔT. Отличие этой задачи от двух предыдущих состоит в том, что здесь рассматривается оптимизация критериев I_C, I_т по двум переменным H и M.

Полученные по результатам расчетов (при отсутствии ветра и в стандартной атмосфере) зависимости M_opt = M(m, CI) и H_opt = H(m, CI) без ветра и при стандартной атмосфере показаны на рис. 7.

Представленные на рис. 7 результаты позволяют сделать следующие выводы:

при выполнении условия 0 < CI < 1.2 значение M_opt практически не зависит от массы ВС, а значение Н_opt  линейно уменьшается с увеличением массы и практически не зависит от значения CI;

при CI > 1.2 значения M_opt и Н_opt проявляют зависимость от обоих аргументов (m, CI), которая усиливается с увеличением значения CI;

при CI = 0, что соответствует оптимизации по критерию расхода топлива, оптимальное значение числа Маха изменяется в узком диапазоне 0.775 < M_opt < 0.78, а оптимальное значение высоты с увеличением веса ВС плавно и без скачков уменьшается с 12 500 до 10 700 м.

Исследование влияния на значения M_opt и Н_opt (полученных в отсутствие возмущений) скорости ветра и температуры воздуха показало, что:

в целом характер зависимостей M_opt = M(m, CI) и H_opt = H(m, CI) не меняется;

наиболее существенное влияние на M_opt и Н_opt оказывает встречный ветер, и оно эквивалентно увеличению значения CI, т.е. изменение скорости ветра от нуля до –100 м/с приводит к таким же изменениям, как и увеличение CI от 0 до 2.0 (кроме того, следует отметить, что влияние ветра на значение Н_opt сказывается только в области больших весов ВС);

температура воздуха практически не влияет на оптимальное значение числа Маха, влияние ΔT на оптимальное значение Н_opt ограничено значением $ \approx {\kern 1pt} 500{\kern 1pt} - {\kern 1pt} 700\;{\text{м}}$ и имеет место в области средних и больших весов ВС.

Оптимизация сценариев 4–6, связанных с переходом на оптимальную высоту либо на оптимальное значение числа М в зависимости от текущего изменения веса ВС, проводится в соответствии с результатами, приведенными на рис. 2–4. В частности, при заданном значении шага изменения числа М по зависимостям рис. 4, 5 можно найти те значения массы, при которых целесообразно увеличивать число М, а по рис. 6 можно определить значения m, при которых следует переходить на следующий эшелон полета.

6.2.3. Оптимизация полета с выполнением заданного времени прибытия в очередную точку маршрута. Формально эту задачу можно свести к классической изопериметрической задаче, где в качестве управлений выступают высота H_зад и скорость M_зад, на которые наложены ограничения:

Однако в данном случае более удобно она решается графо-аналитическим методом, суть которого при заданных значениях дальности полета L и времени прибытия T_зад заключается в следующем:

1. Разобьем диапазоны (6.18) возможного изменения переменных H_зад и M_зад на некоторое количество отрезков точками H_i, $i = \overline {1,{{i}_{h}}} $, M_j, $j = \overline {1,{{j}_{m}}} $.

2. Для каждого значения высоты Н_i и всех значений M_j, $j = \overline {1,{{j}_{m}}} $, получим путем интегрирования системы уравнений (6.16) сечения функции ${{I}_{t}} = F({{M}_{{{\text{зад}}}}}{\text{/}}{{H}_{{{\text{зад}}}}})$ плоскостями H_зад = H_i. Эти сечения показаны на рис. 6, а.

3. Если T_зад принадлежит диапазону значений ${{I}_{t}}$, образованному сечениями ${{I}_{t}} = F({{M}_{{{\text{зад}}}}}{\text{/}}{{H}_{{{\text{зад}}}}})$, то проводим линию ${{I}_{t}} = {{T}_{{{\text{зад}}}}}$ (на рис. 8, а T_зад = 2100 с) и находим точки пересечения ее с полученными сечениями.

4. Для этих точек находим соответствующие им значения числа Маха M_j, $j = \overline {1,{{j}_{m}}} $, которые вместе со значениями H_i, $i = \overline {1,{{i}_{h}}} $ образуют линию  f (H_зад, M_зад) = 0, показанную на рис 8, б. Каждая точка этой линии (по построению) определяет пару значений H_зад, М_зад, при полете с которыми обеспечивается требуемое время прибытия в заданный пункт маршрута.

5. Для каждой точки этой линии путем интегрирования системы уравнений (6.16) находим значения критерия (стоимости) I_C = I_C (H) и значение H = H_opt, в котором стоимость минимальна. Эта точка показана на рисунке 8, в.

6. С этим значением H_opt поднимаемся на рис. 8, б, где по функции  f (H_зад, M_зад) = 0 находим значение M_opt.

Таким образом, проведена оптимизация всех рассмотренных сценариев полета ВС и исследовано влияние параметров полета на значения критериев оптимизации.

Заключение. Как показали результаты моделирования, предлагаемый подход к оптимизации этапа крейсерского полета по критериям стоимости и затратам топлива позволяет получить решение задач для множества типовых сценариев полета ВС. Характер зависимостей оптимального числа М от веса ВС хорошо коррелируется с данными, полученными в работе [12]. Результаты расчетов, проведенные для ВС типа SSJ-100, показывают, что за счет оптимизации высоты полета и числа М можно сократить затраты топлива более чем на 10–15% и время полета – на 3–9%. При этом для обеспечения выполнения ограничений на число М диапазон значений CI следует рассматривать в пределах 0–2 кг/с. Учитывая, что в гражданской авиации экономия 1% топлива считается хорошим результатам, полученные оценки показывают практическую значимость предлагаемой методики решения задачи оптимизации и целесообразность реализации ее при разработке бортовых алгоритмов ССВ.