НАУКА и ОБРАЗОВАНИЕ

УДК (629.7.018.7.004.94):006.63

МГТУ имени Н.Э. Баумана

marmotto@rambler.ru

vesel_foto@mail.ru

marmotto@rambler.ru

Введение

В настоящее время сохраняется актуальность синтеза алгоритмов автоматического управления самолетами на аварийноопасных режимах полета, например, при пилотировании вблизи эксплуатационных ограничений,  при маневрировании на малых высотах в целях недопущения столкновения с землей и т.д. Одним из путей решения данной задачи является включение в состав системы управления прогнозирующего модуля, обеспечивающего расчет оценок прогноза  параметров движения при заданном  управлении. Учитывая высокую динамичность процессов движения современных самолетов, можно предположить, что весьма полезным может оказаться прогнозирование даже  на малых интервалах времени длительностью от нескольких десятых долей секунды до 2…3 секунд. Известно [1], что пространственное движение самолета достаточно полно описывается системой  нелинейных  дифференциальных уравнений 10…12 порядка. Аэродинамические коэффициенты, входящие в это уравнение, зависят от параметров полета и в общем случае описываются  таблично задаваемыми нелинейными функциями, содержащими от нескольких сотен  до нескольких тысяч элементов. Такую модель несложно реализовать на современном персональном компьютере, однако для реализации в бортовом вычислителе и расчетов в реальном масштабе времени она является достаточно громоздкой. Другое, более серьезное, ограничение состоит в том, что современная аэродинамика не гарантирует высокой точности параметров модели, получаемой до начала летных испытаний. Опыт показывает, что отдельные аэродинамические коэффициенты, существенно влияющие на параметры движения самолета, могут иметь погрешности  30…50 процентов, а иногда и более [2, 3]. Банк аэродинамических характеристик может быть проверен и уточнен по результатам летных экспериментов методами теории идентификации, однако это требует проведения большого объема работ, которые в отечественной практике выполняются не для всех самолетов и не всегда  в полном объеме.

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

1 Выбор математической модели движения самолета

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

,

,

,

,

,

,                                                             (1)

,

,

,

.

Ускорения вдоль связанных осей:

a_x= qS(–c_x+c_p)/m–gsinu  = g(n_x–sinu) ,

a_y = qSc_y /m–g cosu  cosγ = g(n_y–cosu cosγ),                                             (2)

a_z = qSc_z /m+ gcosu  sinγ = g(n_z+cosu  sinγ ),

где  c_x, c_y, c_z – коэффициенты  аэродинамических сил  в связанной системе координат.

Перегрузки вдоль связанных осей:

n_x= qS(–c_x+c_p)/gm = a_x/g+ sinu ;

n_y = qSc_y /gm=a_y /g+cosu cosγ ;                                                                (3)

n_z = qSc_z /gm = a_z /g–cosu sinγ .

В системах уравнений (1)…(3) используются следующие обозначения:

  – углы атаки и скольжения, рад;

ω_x, ω_y, ω_z – угловые скорости относительно связанных осей, рад/с;

u, γ, ψ – углы тангажа, крена, рыскания, рад;

V – истинная воздушная скорость, м/с;

H – высота полета, м;

m_x, m_y, m_z – коэффициенты аэродинамических моментов;

c_x, c_y, c_z – коэффициенты аэродинамических сил;

J_x, J_y, J_z , J_xy– моменты инерции относительно связанных осей, кг´м²;

m – масса самолета, кг;

l, b_A – размах крыла и длина средней аэродинамической хорды, м;

S – эквивалентная площадь крыла, м²;

q – скоростной напор, Па;

r_H – плотность воздуха на высоте полета, кг/м³;

c_p=P/qS – коэффициент тяги двигателей;

Р_пр, Р_лeв – сила тяги правого и левого двигателей, Н.

k_дв – кинетический момент роторов двигателей, кг м².

 у_дв, z_дв – координаты двигателя относительно связанных осей, м;

φ _дв – угол установки двигателей, рад.

Значения коэффициентовподъемной силы и сопротивления, определенные по результатам аэродинамических продувок, обычно задаются в полусвязанной системе координат OX_eY_eZ_e .Для перевода данных коэффициентов из одной системы в другую используются соотношения

c_хe=c_х cosα + c_у sin α ,

c_уe=c_у cosα – c_х sin α ,                                                                                                         (4)

c_х=c_хe cosα – c_уe sin α ,

c_у=c_уe cosα + c_хe sin α ,

где  c_xe, c_ye – коэффициенты составляющих аэродинамической силы в полусвязанной системе координат.

В случае необходимости учета  положения самолета в земной системе координат модель (1)  дополняются следующими уравнениями:

=V [cosα cosβ cosu cosψ+sinα cosβ (sinu cosγ cosψ–

–sinγ sinψ)+sinb (sinu sinγ cosψ+cosγ sinψ)];                                 (5)

= –V [cosα cosβ cosu sinψ+sinα cosβ (sinu cosg sinψ+

+sinγ cosψ)+sinβ (sinu sinγ sinψ–cosγ cosψ)].

Линейная модель аэродинамических коэффициентов продольного и бокового движения имеет вид

,

.                                  (6)

,

,                         (7)

,

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

Рассмотренную систему уравнений можно упростить, рассматривая, например, отдельно продольное движение и линейную модель аэродинамических сил и моментов (6). Далее, выбирая в качестве опорной траектории, например, прямолинейный горизонтальный полет с постоянной скоростью, можно получить уравнении продольного короткопериодического движения в приращениях относительно параметров соответствующего опорного движения:

,

.           (8)

Переход от полных уравнений (1) к упрощенной системе (8) и принимаемые при этом допущения  детально рассматриваются в [4]  и других курсах динамики полета.

         Введем обозначения:           

         Тогда уравнения (8) принимают вид

,

.                                             (9)

2 Преобразование модели к эквивалентному дискретному виду

Уравнения (9) с точки зрения теории динамических систем [5] представляют собой линейную стационарную динамическую систему. В векторно-матричной форме их можно записать следующим образом:

                                                                                    (10)

где                .

Для вычислений в реальном масштабе времени на цифровом процессоре целесообразно рассмотреть вместо дифференциальных уравнений (10), заданных в непрерывном времени, эквивалентные дискретные уравнения. Процедура дискретизации линейной стационарной системы хорошо известна [5]. Эквивалентная дискретная рекуррентная формула для стационарной системы (10) принимает вид

                                                                      (11)

Дискретные матрицы зависят в этом случае только от одного аргумента, а именно от разности  При этом интервал дискретизации  предполагается постоянным

Тогда для нахождения эквивалентных дискретных матриц стационарной системы необходимо решить следующие уравнения:

                                                                                                 (12)

Эти уравнения решаются на интервале    с начальными условиями

                                                                                                                (13)

Эквивалентные дискретные матрицы вычисляются только один раз для интервала  Для остальных интервалов  дискретные матрицы будут такими же, ибо свойства стационарной системы не зависят от времени.

Итак, непрерывной системе (10) соответствует следующая эквивалентная дискретная системa, которую  и предлагается  использовать  для прогноза:

,

 ,                                 (14)

где  - элементы эквивалентой дискретной матрицы ;
  - элементы эквивалентой дискретной матрицы , которая в данном случае редуцируется в вектор, поскольку используется только один входной сигнал – отклонения стабилизатора .

Отметим, что в рассматриваемой задаче не требуется решать матричные уравнения (12), поскольку значения параметров в уравнениях (14) предлагается вычислять в реальном масштабе времени по полетным данным на скользящем интервале по известным формулам множественной линейной регрессии [6].

3 Результаты апробации по данным летного эксперимента

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

Далее для всех изображений  красная линия показывает исходные данные, полученные в летном эксперименте,  синяя линия – экстраполированные данные. Красными маркерами-окружностями выделен участок исходных данных, используемый для расчета параметров модели. По оси абсцисс отложено число отсчетов.  В эксперименте частота регистрации составляла 16 Гц, следовательно, каждые 16 отсчетов соответствуют интервалу времени 1 с.

Рассмотрим несколько вариантов поставленной задачи.

Известный входной сигнал.  В уравнениях (14) в качестве входного сигнала модели рассматриваются отклонения стабилизатора . Предположим, что  вычислитель системы  управления рассчитал входной сигнал на некоторый интервал времени вперед и требуется определить, в каком состоянии окажется объект при таком управлении. Если состояние окажется неудовлетворительным, управление необходимо скорректировать. Получим оценку точности предлагаемого алгоритма прогнозирования в этом случае. На рис. 1 и 2 показаны результаты экстраполяции на 6,25 с (100 отсчетов) для сигналов угла атаки  и угловой скорости тангажа . Значения среднеквадратического отклонения (с.к.о.) погрешности экстраполяции оценивалось как  с.к.о. разности между вычисленным по модели прогнозом и полученным в эксперименте значением соответствующего сигнала. Графики показывают весьма близкое соответствие по углу атаки. Для угловой скорости имеет место постоянная погрешность величиной  приблизительно 2 градуса. Из рис. 2 следует, что после исключения этой погрешности ( например, путем расширения вектора идентифицируемых параметров) точность прогноза станет сопоставимой с точностью канала угла атаки. 

Рисунок 1 – График изменения угла атаки во времени при экстраполяции на 100 отсчетов (6,25 с), с.к.о. погрешности экстраполяции  0,8517 градусов

Рисунок 2 – График изменения угловой  скорости тангажа во времени при экстраполяции на 100 отсчетов (6,25 с), с.к.о. погрешности экстраполяции 2,1057 градус/с

На рис. 3 и 4 показаны графики зависимости погрешности экстраполяции от  числа отсчетов, что эквивалентно длительности участка экстраполяции. Из графика на рис.  3 следует, что при прогнозировании на интервале длительностью до 3 с  (48 отсчетов) с.к.о. погрешности находится в удовлетворительных пределах 0,4…0,6 градусов. Погрешность канала угловой скорости достигает при этом 2,5 градусов/с за счет упомянутой выше неисключенной постоянной погрешности, что указывает на необходимость коррекции модели.

Рисунок 3 – График зависимости с.к.о. погрешности экстраполяции угла атаки от количества отсчетов экстраполяции

Рисунок 4 – График зависимости с.к.о. погрешности экстраполяции угловой  скорости тангажа от количества отсчетов экстраполяции

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

В теоретическом смысле это означает следующее. Известно, что выходной сигнал системы (10) на произвольном интервале времени определяется интегралом Дюамеля [5]:

,                          (15)

где  входной сигнал. Отметим, что рекуррентная формула (11) является дискретным эквивалентом операции интегрирования  в (15).

В первом случае (известный входной сигнал) при прогнозировании использовались оба слагаемых правой части формулы (15). В рассматриваемом  случае, полагая приращения входного сигнала нулю, мы  обнуляем и второе слагаемое в (15), поскольку используемая модель  (10) сформулирована в приращениях. Это означает, что прогноз основывается только на значении вектора состояния  в момент начала экстраполяции и на свойствах системы, учитываемых в матрице . Следовательно, удовлетворительная точность прогнозирования будет обеспечена, пока движение определяется начальными условиями, то есть первым слагаемым (15). Очевидно, что интервал точного прогнозирования будет больше при менее энергичном маневрировании, то есть при меньших значениях второго слагаемого, которое считается неизвестным. Результат прогнозирования для угла атаки на 50  отсчетов (3,125 с) показан на рис. 5. Высокая точность сохраняется в течение 2 с, то есть до начала второй серии отклонений РУС. После появления нового энергичного управляющего сигнала модель теряет адекватность, поскольку в ней это сигнал не учитывается.

Рисунок 5 – График изменения угла атаки во времени при экстраполяции на 50  отсчетов (3,125 с), с.к.о. погрешности экстраполяции  6,1 градусов

Использование отклонений РУС в качестве входного сигнала. Современные самолеты оснащены достаточно сложными системами

автоматического управления, поэтому отклонения стабилизатора не совпадают с отклонениями РУС по тангажу вследствие наличия обратных связей. С другой стороны, алгоритмы автоматизации пилотажного контура синтезируются из соображений удобства пилотирования, которое, как правило, заключается в том, чтобы самолет в определенном смысле «следовал» за рукояткой управления.

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

Результаты по углу атаки для времени прогнозирования 6,25 с представлены на рис. 6. Сравнение с рис. 1 показывает, что уровень погрешностей увеличивается незначительно.

Рисунок 6 – График изменения угла атаки при экстраполяции на 100 отсчетов (6,25 с), с.к.о. погрешности 1,03 градуса (входной сигнал - РУС)

Сравнение обобщенного графика зависимости с.к.о. погрешностей от длительности участка экстраполяции  (рис. 7) с аналогичным графиком на рис. 3 показывает, что отличия в с.к.о. не превышают 0,2 градуса, а для интервала времени до 3 с даже несколько меньше, чем при использовании в качестве входного сигнала отклонений стабилизатора.

Выполнение расчетов для неизвестного входного сигнала при использовании РУС по тангажу также показало результаты, близкие к полученным выше для случая использования отклонений стабилизатора.

Рисунок 7 – График зависимости с.к.о. погрешности экстраполяции угла атаки от количества отсчетов (входной сигнал - РУС)

Заключение

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

Апробация алгоритма по результатам летных испытаний одного из современных самолетов показала, что на участке прогнозирования длительностью до 3 с при известном управляющем сигнале  с.к.о. погрешности прогнозирования угла атаки не превышает 0,4…0,6 градусов.

Аналогичные показатели точности достигаются при использовании в качестве входного сигнала отклонений РУС по тангажу вместо отклонений стабилизатора. Этот результат обусловлен наличием сильной корреляционной связи между отклонениями РУС и стабилизатора, создаваемой современными системами управления.

Работа выполнена при поддержке РФФИ, проект 11-08-00292-а.

Литература

1.   Аэродинамика, устойчивость и управляемость сверхзвуковых  самолетов. / Под ред. Бюшгенса Г.С. – М.: Наука. Физмалит. 1998. 811 с.

2.   Корсун О.Н. Принципы параметрической идентификации  математических моделей самолетов по данным летных испытаний. – М.: Мехатроника, автоматизация, управление. Приложение-2008-06, с. 2-7.

3.   Корсун О.Н., Поплавский Б.К. Структура методологии идентификации  математических моделей самолетов по результатам летных испытаний. - Новые рубежи авиационной науки ASTEC’07. Международная научно-техническая конференция: Сб. научн. тр. М., ЦАГИ, 2007. 1 электрон. опт. диск (CD-ROM).

4. Лебедев А.А., Чернобровкин Л.С.  Динамика полета беспилотных летательных аппаратов. – М.: Машиностроение, 1982. 635 c.

5. Корсун О.Н.  Методы математического моделирования технических  систем: электронное учебное пособие - М., МГТУ им. Н.Э. Баумана. http://db.inforeg.ru/deposit/Catalog/mat.asp?id= 284494, зарегистрировано во ФГУП "Информрегистр" № 0321002457, 2010 г.

6.   Корсун О.Н. Методы параметрической идентификации технических  систем: электронное учебное пособие. М., МГТУ им. Н.Э. Баумана. http://db.inforeg.ru/deposit/Catalog/mat.asp?id= 286062, зарегистрировано во ФГУП "Информрегистр" № 0321100941, 2011 г.