Прогнозирование параметров движения самолета на основе идентификации упрощенной линейной модели Текст научной статьи по специальности «Математика»

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

Эл № ФС 77 - 30569. Государственная регистрация №042110002.5. 155И 1994-0408

Прогнозирование параметров движения самолета на основе идентификации упрощенной линейной модели 77-30569/290540

# 12, декабрь 2011

Корсун О. Н., Веселов Ю. Г., Гулевич С. П.

УДК (629.7.018.7.004.94):006.63

МГТУ имени Н.Э. Баумана marmotto@rambler.ru уеБе! foto@mail.ru marmotto@rambler.ru

Введение

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

30...50 процентов, а иногда и более [2, 3]. Банк аэродинамических характеристик может быть проверен и уточнен по результатам летных экспериментов методами теории идентификации, однако это требует проведения большого объема работ, которые в отечественной практике выполняются не для всех самолетов и не всегда в полном объеме.

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

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

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

da

dt

1

ю г -

cos в

dв az

dt V

V

а„

-ю у sin в

= —соя в - —яіп в - т

V

sm а +

cos а +

a

\

— + ю х sin в

Vх

соя а

у

ay

— яіп в + тх

V х

sm а

dV

dt

ax cos а cos в - av sin а cos в + az sin в,

dю z _Лх - у

dt

Л

ю х ю у + а

8_Ь

Л

А

т„

^дв юу ( Рпр + Рлев ) У дв

Л Л

ю Jz Лх ю ю . т . kдвюz . (РпР Рлев ^z дв

------=------юхю7 + а— ту +------------------------1--------------------

Л Лу Х Z Лу У Лу Лу

dю х Лу - Л 81

—- = —---------------ю у ю z + а — тх,

dt Л У Z Л Х

du

dt

= ту яіп у + т2 соя у.

dy

dt

тх - tg и (т у соя у - тг яіп у),

dy 1

dt соя и

(ту соя у - т2 яіп у).

н

dt

= V [ соя а соя в яіп и - яіп а соя в соя и соя у - яіп в соя и яіп у ].

(1)

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

ах = а8(-сх+ср)/т^ яіпи = ^(пх-яіпи) , ау = д8су /т^ соя и сояу = g(nv-coяu сояу),

az = /т+ gcoяu яіпу = g(nZ+coяu яіпу ),

где сх, су, с:, - коэффициенты аэродинамических сил в связанной системе

координат.

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

пх = д8(-сх+ср)/§т = ахяіпи;

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

ах, ау, а2 - угловые скорости относительно связанных осей, рад/с; и, у, ц - углы тангажа, крена, рыскания, рад;

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

Н - высота полета, м;

тх, ту, ш2 - коэффициенты аэродинамических моментов; сх, су, с2 - коэффициенты аэродинамических сил;

Зх, Jy, Jz , Зху- моменты инерции относительно связанных осей, кгхм ; т - масса самолета, кг;

I, ЪА - размах крыла и длина средней аэродинамической хорды, м;

S - эквивалентная площадь крыла, м2; д - скоростной напор, Па;

рН - плотность воздуха на высоте полета, кг/м3;

Ср=Р/д$> - коэффициент тяги двигателей;

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

кдв - кинетический момент роторов двигателей, кг м2.

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

Ф дв - угол установки двигателей, рад.

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

пу = д8су/^=ау^+сояи сояу ; п = /^ = az ^-сояи яіпу .

(3)

cxe=cx cosa + cy sin a .

cye=cy cosa - cx sin a , (4)

cx=cxe cosa - cye sin a

cy=cye cosa + cxe sin a

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

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

dX

----=V [cosa cosв cos и cos^+sina cosв (sinu cosy cos^-

dt

-siny sin^+sin/? (sinu siny cos^+cosy sin^)]; (5)

dZg

----= -V [cosa cosв cosu sin^+sina cosв (sinu cos^ sin^ +

dt

+siny cos^+s^ (sinu siny sin^-cosy cos^)].

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

Суе = Су 0 + <> + СУВ Фв ,

а 8& У Ъ А а Ъ А

т2 = тг0 + т2 а + тгв Фв + тУ '-АУ + тг ■ . (6)

V V ш

сz = свв + 4- 8н,

IV 2V

в п а 1 ау 1 д ? д

ту = ту в + тух — ах + туу — ау + ту-дн + туэдэ, у уг у 2V х у 2V у у н у э

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

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

Аа(,) - С,Ч^

тУ тУ

Ла'(і) = Лог (і ) - еау—— Ла(і) - с4’—— Лрв (і),

Л02 (і) = < Ла(і) + тр Лрв (і) + .Ьа Лю2 (і). (8)

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

Введем обозначения: У а = сау; УР = сР;

у у т^ у уmV

л,а а Ч^ЪА р ц8ЪА т0’7-ц8ЪА Ьа

ма = та -—А; МР = тР -—А;М0 =—А .—;

^ г ^ г ^ г *

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

Ла'(і) = Ла>г (і) - У аЛа(і) - У РЛрв (і),

Л0г (і ) = М а Ла(і) + М р Лрв (і ) + М 00 Лог (і ). (9)

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

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

У( і) =Ау( і )+Ви( і),

(10)

где у(і )=

Ла(і) Л0 (і )

"- Уа 1 " "- У р

, А= , в=

Маг М°г 1 М 1

и (і )=ЛРв (і )

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

Дискретные матрицы зависят в этом случае только от одного аргумента, а именно от разности ^+1 - ^ = Т. При этом интервал дискретизации Т предполагается постоянным

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

Эквивалентные дискретные матрицы вычисляются только один раз для интервала [0, Т]. Для остальных интервалов [їк, ік + ] дискретные матрицы будут такими же, ибо

свойства стационарной системы не зависят от времени.

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

где /п, /12, /21, /22 - элементы эквивалентой дискретной матрицы Ф(Т); g1, g2 - элементы эквивалентой дискретной матрицы Г(Т), которая в данном случае редуцируется в вектор, поскольку используется только один входной сигнал - отклонения стабилизатора рв (?к ) .

У(ік+і) = Ф(Т)У(ік ) + Г(Т)и(ік ).

(11)

ік+і - ік = Т, к = 0,1,2,...N - 1.

(12)

Эти уравнения решаются на интервале і є [0, Т] с начальными условиями

Ф(0) = Е, Г(0) = 0.

(13)

Ла(ік+1) = /цЛа(ік ) + /12Л0г (ік ) + §1ЛРв (ік ) >

Л0г (ік+1 ) = /21Ла(ік ) + /22Л0г (ік ) + ^2ЛРв (ік ) ,

(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]:

*к + у

У(1к+у ) = Ф(*к+у , *к ) У(*к ) + |Ф0к+у В(Т> (т)Т (15)

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

В первом случае (известный входной сигнал) при прогнозировании использовались оба слагаемых правой части формулы (15). В рассматриваемом случае, полагая приращения входного сигнала нулю, мы обнуляем и второе слагаемое в (15), поскольку используемая модель (10) сформулирована в приращениях. Это означает, что прогноз

основывается только на значении вектора состояния у(їк) в момент начала

экстраполяции и на свойствах системы, учитываемых в матрице Ф(?к+у., їк).

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

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

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

пилотирования, которое, как правило, заключается в том, чтобы самолет в определенном смысле «следовал» за рукояткой управления.

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

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

8

6

4

2 0 -2 -4 -6 -8 -10

0 20 40 60 80 100 120 140 160 180 200

Рисунок 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 электрон. опт. диск (СБ-КОМ).

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

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 г.

electronic scientific and technical periodical

SCIENCE and EDUCATION

___________El. .Vs KS 77 - 30569. -V«042l 100025. ISSN 1994-0408_

Aircraft movement prediction using a linear model identification 77-30569/290540

# 12, December 2011

Korsun O.N., Veselov Yu., G., Gulevitch S., P.

Bauman Moscow State Technical University marmotto@rambler.ru vesel foto@mail.ru marmotto@rambler.ru The authors consider a simplified linear model of aircraft pitch movement for real-time

prediction to support control design. The model is derived out of the full nonlinear system of

aircraft space movement. The model parameters estimates are achieved through the least square

regression estimator using moving interval. Some examples of flight test data processing are also

presented.

Publications with keywords: system identification, mathematical modeling, flight tests, aircraft control design, control system, parameter identification, real-time prediction Publications with words: system identification, mathematical modeling, flight tests, aircraft control design, control system, parameter identification, real-time prediction

Reference

1. In: G.S. Biushgens (Ed.), Aerodynamics, stability and controllability of supersonic aircraft, Moscow, Nauka. Fizmalit, 1998, 811 p.

2. Korsun O.N., Mekhatronika, avtomatizatsiia, upravlenie. Supplement 06 (2008) 2-7.

3. Korsun O.N., Poplavskii B.K., in: Proc. of the International scientific-technical conference of the New frontiers of the aviation science ASTEC’07, Moscow, TsAGI, 2007.

4. Lebedev A.A., Chernobrovkin L.S., Dynamics of flight of unmanned aerial vehicles, Moscow, Mashinostroenie, 1982, 635 p.

5. Korsun O.N., Methods of mathematical modeling of technical systems: electronic textbook, Moscow, MGTU im. N.E. Baumana - BMSTU, 2010 <http://db.inforeg.ru/deposit/Catalog/mat.asp?id= 284494>.

6. Korsun O.N., Methods of parametric identification of technical systems: electronic textbook, Moscow, MGTU im. N.E. Baumana - BMSTU, 2011 <http://db.inforeg.ru/deposit/Catalog/mat.asp?id= 286062>.