Quasi-linearization of a nonlinear model of an unmanned aerial vehicle

Мұқаба

Дәйексөз келтіру

Толық мәтін

Аннотация

The article discusses the problem of representing a nonlinear dynamic system in the state space using the quasi-linearization method. The developed method for constructing a quasi-linear model is based on the principles of structural-parametric synthesis. The first stage of construction consists in the analytical output of the structure of the quasi-linear model from nonlinear differential equations with the identification of variable parameters, on which the dynamic characteristics of the control object depend. At the stage of parametric synthesis, identification of unknown elements of the resulting structure is carried out in the form of functional dependencies on the corresponding variables. For this, the task of approximation in the form of a linear generalized model using algorithms of sparse regression and heuristic decomposition of this model on the necessary functional dependencies is solved. As an example, the quasi-linearization method is applied to the nonlinear model of the isolated movement of the unmanned aerial vehicle in the longitudinal channel relative to the center of the masses. The proposed method can be applied to other models with similar properties to solve problems of controlling nonlinear dynamic systems with modern methods of control theory.

Толық мәтін

Введение

Общепринятым подходом к решению задач управления является линеаризация нелинейной динамической системы относительно некоторых опорных точек, характеризующихся рабочими режимами применения объекта. Полученная таким образом линеаризованная модель является лишь приближением исходной системы и не учитывает все аспекты нелинейной динамики. Применительно к беспилотным летательным аппаратам (БПЛА) такими особенностями могут быть: аэродинамические «ложки», срыв потока на больших углах атаки и т.д. Достижение высокого качества управления во всем диапазоне изменения динамических характеристик предполагает обеспечение наиболее точного представления нелинейной модели объекта в необходимой для синтеза регуляторов форме. Решение рассмотренных задач целесообразно производить с помощью методов современной теории управления, что подразумевает представление системы в пространстве состояний.

Перспективным направлением в представлении нелинейных динамических систем в пространстве состояний является преобразование исходной системы в квазилинейную форму с изменяющимися параметрами (quasi-Linear Parameter Varying – qLPV) специальным образом [1]. В зависимости от априорной информации о динамике исследуемого объекта управления могут быть применимы различные методы идентификации систем, например, аналитический вывод, численные методы [2], квазилинеаризация [3], аппроксимация с использованием мягких вычислений (искусственных нейронных сетей, нечеткой логики и т.д.). Перечисленные методы, за исключением последнего, относятся к классу идентификации в узком смысле, что подразумевает наличие достаточно подробной информации о нелинейной динамике. Аппроксимация с помощью мягких вычислений, напротив, относится к классу идентификации в широком смысле и требует только достаточного объема данных в виде «вход-выход». Точность такой аппроксимации существенно зависит от настроек выбранного инструмента: типа и структуры нейронной сети, базы правил нечеткой модели и т.д. В некоторых случаях система нелинейных дифференциальных уравнений рассматриваемого объекта управления полностью известна, но не соответствует нормальной форме Коши. В результате структурного преобразования к указанной форме появляются неизвестные переменные, которые необходимо идентифицировать по имеющимся исходным данным.

Целью настоящей работы является разработка эффективного метода квазилинеаризации нелинейной модели БПЛА с высокой точностью представления динамической системы по исходным данным о характеристиках объекта управления.

1. Структурный синтез квазилинейной модели

Для решения задач стабилизации относительно центра масс удобно использовать нелинейную модель пространственного движения БПЛА как твердого тела в проекциях на оси скоростной системы координат. В некоторых случаях при необходимости можно не учитывать взаимное влияние продольного и бокового каналов и пренебречь перекрестными связями между ними. Так, например, нелинейная модель изолированного движения БПЛА в продольном канале может иметь следующий вид:

α˙(t)=ωz(t)Ya(t)mV,ω˙z(t)=Mz(t)Jz,ϑ˙(t)=ωz(t),ny(t)=Ya(t)mg(H), (1)

где α(t) – угол атаки [рад], ωz(t) – угловая скорость тангажа [рад/с], ϑ(t) – угол тангажа [рад], ny(t) – нормальная перегрузка [g], Ya(t) – подъемная аэродинамическая сила [Н], Mz(t) – аэродинамический момент тангажа [Н∙м], Jz – момент инерции тангажа [кг∙м2], g(H) – ускорение свободного падения [м/с2], H – высота [м], V – модуль воздушной скорости [м/с], m – масса БПЛА [кг].

Модель рулевого привода для гипотетического БПЛА может быть описана апериодическим звеном и иметь вид:

Tδ˙в+δв=δвзад, (2)

где δв – угол отклонения руля высоты [рад], δвзад – управляющий сигнал в продольном канале [рад], T – постоянная времени привода [с].

Нелинейная динамическая система может быть представлена в пространстве состояний с помощью qLPV модели:

x˙(t)=A(θ(t))x(t)+B(θ(t))u(t),y(t)=C(θ(t))x(t)+D(θ(t))u(t), (3)

где x(t)m – вектор состояния, y(t)l – вектор измерений, u(t)k – вектор управления, θ(t)=[θ1(t)θ2(t)θN(t)]ΩN – вектор варьируемых параметров (в данном случае зависящий в т. ч. от состояния), N – количество варьируемых параметров, Ω=[θ1minθ1max]×[θ2minθ2max]××[θNminθNmax] – замкнутый гиперкуб, θnmin,θnmax – диапазоны изменения варьируемых параметров.

Чтобы представить систему (1) в виде (3) сначала определим варьируемые параметры, от которых зависят характеристики нелинейной модели. Далее будем считать все параметры модели зависимыми от времени t за исключением Jz, g, H, V, m, которые в рамках короткопериодического движения предполагаются постоянными величинами.

Согласно [4], аэродинамические силу и момент можно представить в виде:

Ya=cyqS, Mz=mzqSba, (4)

где cy, mz – аэродинамические коэффициенты силы и момента, q – скоростной напор [Н/м2], S, ba – характерные площадь [м2] и размер [м].

Аэродинамические коэффициенты сил и моментов гипотетического БПЛА в общем случае являются функциями многих переменных:

cy=cy(M,α,β,δв,δн,δэ), mz=mz(M,α,β,δв,δн,δэ), (5)

где M – число Маха, δн, δэ – углы отклонения рулей направления и элерона [рад].

Пренебрежем перекрестными связями между каналами и проведем разложение рассмотренных коэффициентов в следующей форме:

cy(M,α,δв)=cyα(M,α,δв)α+cyδв(M,α,δв)δв+cy0(M),mz(M,α,δв)=mzα(M,α,δв)α+mzδв(M,α,δв)δв+mz0(M). (6)

Подставляя выражения (4) в систему (1) с учетом (6), q = ρHV2 / 2 и V = aHM, нелинейную модель изолированного движения в продольном канале представим в виде:

α˙=cyα(M,α,δв)ρHaHMS2mα+ωzcyδв(M,α,δв)ρHaHMS2mδвcy0(M)ρHaHMS2m,ω˙z=mzα(M,α,δв)ρHaH2M2Sba2Jzα+mzδв(M,α,δв)ρHaH2M2Sba2Jzδв++mz0(M)ρHaH2M2Sba2Jz,ϑ˙=ωz,ny=cyα(M,α,δв)ρHaH2M2S2mg(H)α+cyδв(M,α,δв)ρHaH2M2S2mg(H)δв++cy0(M)ρHaH2M2S2mg(H), (7)

где ρH – плотность воздуха как функция от высоты H [кг/м3], aH – скорость звука на высоте H [м/с].

Как видно из анализа системы уравнений (7), характеристики нелинейной модели зависят от множества переменных, которое можно представить в виде вектора варьируемых параметров для продольного канала:

θп=HMαδв.

Далее, предположив, что переходной процесс БПЛА совершается из установившегося движения (т.е. движение с нулевыми значениями начальных сил и моментов), статические составляющие cy0(M), mz0(M) будем считать равными нулю. Теперь можно записать нелинейную модель изолированного движения в продольном канале (7) вместе с уравнением динамики рулевого привода (2) в qLPV форме (3):

αωzδв.=a4(θп)1a5(θп)a2(θп)0a3(θп)00T1αωzδв+00T1δвзад,nyωz=a6(θп)0a7(θп)010αωzδв, (8)

где

a2(H,M,α,δв)=mzα(M,α,δв)ρHaH2M2Sba2Jz, a3(H,M,α,δв)=mzδв(M,α,δв)ρHaH2M2Sba2Jz,a4(H,M,α,δв)=cyα(M,α,δв)ρHaHMS2m, a5(H,M,α,δв)=cyδв(M,α,δв)ρHaHMS2m,a6(H,M,α,δв)=cyα(M,α,δв)ρHaH2M2S2mg(H), a7(H,M,α,δв)=cyδв(M,α,δв)ρHaH2M2S2mg(H).

2. Параметрический синтез квазилинейной модели

Пусть исходные данные об аэродинамических коэффициентах гипотетического БПЛА представлены в общей форме (5). Тогда, на основании синтезированной структуры qLPV модели (8), необходимо провести идентификацию неизвестных параметров

cyα(M,α,δв), cyδв(M,α,δв), mzα(M,α,δв), mzδв(M,α,δв)

в виде функциональных зависимостей от соответствующих переменных.

Одним из вариантов получения разложения коэффициентов в необходимой форме (6) является аппроксимация исходных зависимостей (5) специальным образом. Для реализации подобного преобразования заметим, что в первом приближении коэффициенты аэродинамических сил и моментов могут быть представлены гладкими всюду непрерывными функциями, которые зависят лишь от малого числа входных параметров. Следовательно, нет необходимости рассматривать все возможные нелинейные комбинации, а достаточно вести поиск в малом наборе зависимостей требуемого вида [5, 6].

В результате задача параметрического синтеза может быть сведена к задаче аппроксимации с помощью линейной обобщенной модели:

f(x)i=1mθi(x)ξi=Θ(x)Ξ, (9)

где f(x)=f1(x1),,fn(xn)n×1 – наборы значений аппроксимируемой зависимости для соответствующих наборов входных параметров x=x1,,xnn×1xi=p1,,pb1×b – входные параметры i-ого набора для fi(xi), Θ(x)=θ1(x),,θm(x)n×m – библиотека функций кандидатов, θi(x)=θi(x1),,θi(xn)n×1 – наборы значений i-ой функции кандидата, Ξ=ξ1,,ξmm×1 – вектор коэффициентов перекрестных связей.

Модель (9) вычисляется путем нахождения вектора коэффициентов Ξ на основе задачи выпуклой разреженной регрессии c регуляризацией по L1 норме:

Ξ=argminΞf(x)Θ(x)Ξ2+εΞ1, (10)

где ε – настраиваемый параметр, характеризующий разреженность.

Для решения задачи (10) может быть использован последовательный метод наименьших квадратов с порогом [7]. Далее представлено его описание.

Алгоритм 1. Последовательный метод наименьших квадратов с порогом

Вход: f(x) – наборы значений аппроксимируемой зависимости, x – наборы входных параметров, Θ(x) – библиотека функций кандидатов, ε – параметр разреженности модели, kmax – максимальное число итераций.

Выход: Ξ – вектор коэффициентов перекрестных связей.

  1. Вычислить начальное приближение Ξ0 = Θ−1(x)f(x);
  2. Цикл k < kmax
  3. k = k +1;
  4. Найти коэффициенты меньше установленного порога Iε = |Ξk| < ε;
  5. Исключить из отбора малые коэффициенты Ξk(Iε) = 0;
  6. Найти оставшиеся коэффициенты Is = ~ Iε;
  7. Сформировать новую k оценку Ξk(Is) = Θ−1(Is, :)f(x);
  8. Конец цикла

Выбор нужного набора функций кандидатов зависит от вида необходимого представления f(x). Тогда для получения формы (6) удобно использовать следующую полиноминальную модель:

L(x)=i=0Pmaxj=0Pmaxl=0Pmaxξap1ip2jpblθa, a=a+1, a=1,,m, (11)

где Pmax – максимальная степень полинома. Причем каждая уникальная комбинация степеней p1ip2jpbl из (11) соответствует определенной функции кандидату θa из (9).

Для получения необходимого разложения по форме (6) ряд (11) необходимо представить в виде:

L(x)=i=1MLi(x)pli, (12)

где M – число разложений, Li(x) – полиноминальные модели, входящие в состав (11), определяемые в результате разложения, pl1,pl2,,plM – параметры, для которых нужно найти разложение.

Модель по форме (11) допускает разнесение параметров по форме (12). Однако в общем случае их количество может быть достаточно большим, что усложняет эвристический способ разложения. Предположим, что наибольший вклад в изменение аэродинамических коэффициентов вносят элементы малых порядков без перекрестных связей между входными параметрами, тогда процесс разложения может быть автоматизирован для исключения ошибок, связанных с «ручным» разбиением.

Если коэффициент ξa при функции кандидате θa учитывает несколько входных параметров и является достаточно большим, его следует дополнительно разложить на несколько составляющих, которые могут быть равномерно разнесены между всеми Li(x)pli, если pli входит в θa.

Далее представлен эвристический алгоритм разложение ряда, облегчающий переход модели из формы (11) в (12). Для удобства работы полиноминальная модель L(x), множество участвующих в разложении аргументов P, библиотека функций кандидатов θa реализованы в алгоритме в виде строковых массивов. В результате, задача параметрического синтеза сводится к решению задачи (10) на основе алгоритмов 1 и 2.

Алгоритм 2. Алгоритм эвристического разложения функциональной зависимости

Вход: L(x) – полиноминальная модель в форме (11), P – множество аргументов p, участвующих в разложении, M – количество параметров, участвующих в разложении.

Выход: Li(x) – полиноминальные модели из (12).

  1. Инициализация Li(x) = [], i = 1, ..., M;
  2. Цикл по всем элементам θa из L(x);
  3. Цикл по всем элементам pli из P;
  4. θa'=θa;
  5. Если θa'pli=
  6. Перейти на новую итерацию цикла по pli;
  7. Конец если
  8. Понизить степень элемента pli в θa';
  9. Сформировать множество R, не включающее в себя pli: R=P\pli;
  10. Если θa'R= или степень pli больше степеней ∀piR входящих в θa'
  11. Li(x)=Li(x)+ξaθa';
  12. Конец если
  13. Если степень pli равна степеням ∀piR входящих в θa'
  14. ξa'=ξa/M;
  15. Цикл по всем элементам plk из P
  16. θa'=θa;
  17. Понизить степень элемента plk в θa';
  18. Lk(x)=Lk(x)+ξa'θa';
  19. Конец цикла
  20. Конец если
  21. Перейти на новую итерацию цикла по θa;
  22. Конец цикла
  23. a = a +1;
  24. Конец цикла 

3. Численное моделирование

Проведем численное моделирование для проверки адекватности разработанного метода. В качестве параметров модели, соответствующих (8), использованы характеристики объекта из работы [8]. Тогда аэродинамические коэффициенты подъемной силы и момента тангажа могут быть представлены следующими функциональными зависимостями:

cy(M,α,δв)=0,1741α+0,004766α20,01185Mα++0,1155δв+0,0008479αδв0,0182Mδв,mz(M,α,δв)=0,0294α0,0006682α20,005962Mα0,04525δв0,0005829αδв0,007453Mδв. (13)

Выбранные диапазоны определения параметров модели (13) представлены в табл. 1.

 

Таблица 1. Диапазоны выходных параметров

Table 1. Output parameter ranges

Параметр

Наименование

Диапазон

Единицы измерения

1

2

3

4

M

Число Маха

[3 6]

α

Угол атаки

[-12 12]

[град]

δв

Угол отклонения руля высоты

[-15 15]

[град]

 

Проведем разложение исходных функций (13) в соответствии с формой (6). Максимальная степень полиноминального члена, участвующего в разложении, составляет Pmax = 4. Применив алгоритмы 1 и 2, получим следующие функциональные зависимости:

cy(M,α,δв)cyα(M,α)α+cyδв(M,α)δв,mz(M,α,δв)mzα(M,α)α+mzδв(M,α)δв, (14)

где

cya(M,α)=0,17410,01185M+0,004766α,cyδв(M,α)=0,11550,0182M+0,0008479α,mzα(M,α)=0,0294+0,005962M0,0006682α,mδв(M,α)=0,04525+0,007453M0,0005829α.

Точность полученной аппроксимации (14) относительно исходных данных оценивается по критериям (15) и для равномерно распределенной случайной выборки ns =106 представлена в табл. 2.

Emax=maxi=1,...,nmaxk0kАk01100%,Emean=1nsi=1nsmaxk0kАk01100%, (15)

где ns – размер случайной выборки, k0, kА – исходный и полученный в результате аппроксимации аэродинамические коэффициенты.

 

Таблица 2. Точность аппроксимации аэродинамических коэффициентов

Table 2. Approximation accuracy of aerodynamic coefficients

Ошибка

cy, [%]

mz, [%]

1

2

3

Максимальная (Emax)

 

 

Средняя (Emean)

 

 

Согласно результатам, полученное разложение (14) по форме (6) обеспечивает высокую точность аппроксимации исходных функций (13). Проведем сравнительный анализ путем численного моделирования исходной нелинейной модели (1) и ее qLPV формы (8). В качестве условий моделирования выбран фиксированный режим с высотой полета H = 23 км и числом Маха M = 5. Соотношения параметров объекта для такого режима представлены в табл. 3.

 

Таблица 3. Соотношение параметров объекта управления на выбранном режиме полета 

Table 3. The ratio of the parameters of the control object in the selected flight mode

Параметр

Значение

QSLJz

3643,4

QSmV

0,6053

Vg

151,1792

T

0,015

 

Управляющий сигнала на рассматриваемом режиме определен следующим законом:

δвзад=0,5730sin(6,52πt)+0,5730cos(0,22πt)+1,4325cos(22πt).

Результат моделирования нелинейной модели (1) с учетом динамики привода (2) и ее qLPV формы (8) представлен на рис. 1. Параметры с индексом «н» соответствуют состоянию и выходу нелинейной модели, а с индексом «q» – qLPV модели. В качестве начальных условий использовалось следующее положение на фазовой плоскости:

α0 = −0,5 град, ωz0 = −5,5 град/с, δв0 = 0 град.

 

Рис. 1. Моделирование динамик нелинейной модели и ее qLPV формы

Fig. 1. Modeling of the dynamic of the nonlinear model and its qLPV form

 

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

Заключение

Разработан метод квазилинеаризации нелинейной модели БПЛА с использованием подхода структурно-параметрического синтеза. Предложенный подход основан на аналитическом выводе структуры qLPV модели с последующей идентификацией неизвестных переменных в виде аппроксимации линейной обобщенной моделью посредством решения задачи разреженной регрессии и эвристического разложения. В ходе разработки метода проведена квазилинеаризация нелинейной модели изолированного движения БПЛА в продольном канале относительно центра масс. Результаты численного моделирования показали минимальные ошибки по соответствующим фазовым координатам исходной модели и ее qLPV формы. Полученные результаты позволяют заключить, что разработанный метод обеспечивает высокую точность представления нелинейной динамической системы в виде qLPV модели.

×

Авторлар туралы

A. Shabashov

TEMP-AVIA Arzamas Scientific and Production Enterprise

Хат алмасуға жауапты Автор.
Email: aa.shabashov@mail.ru
ORCID iD: 0009-0000-8320-5201

инженер 1 категории

Ресей, Arzamas

A. Plotnikov

TEMP-AVIA Arzamas Scientific and Production Enterprise

Email: artyom152rus@yandex.ru
ORCID iD: 0009-0007-6688-5517

инженер

Ресей, Arzamas

Әдебиет тізімі

  1. P. Baranyi, Y. Yam, P. Varlaki. Tensor Product Model Transformation in Polytopic Model-Based Control. – CRC Press, Florida, 2013.
  2. Гайдук, А.Р. Численный метод синтеза квазилинейных моделей нелинейных объектов // Мехатроника, автоматизация, управление. 2021. Т. 22. № 6. С. 283-290.
  3. Гроп, Д. Методы идентификации систем / Д. Грон. – М.: Мир, 1979. – 302 с.
  4. Ефремов, А.В. Динамика полета: учебник для студентов высших учебных заведений / А.В. Ефремов, В.Ф. Захарченко, В.Н. Овчаренко. – М.: Машиностроение, 2011. – 775 с.
  5. K. Kaheman, J.N. Kutz, S.L. Brunton SINDy-PI: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics. Proceedings of the Royal Society A. 2020. Т. 476. №. 2242. С. 20200279.
  6. Lee J.D. et al. Communication-efficient sparse regression //Journal of Machine Learning Research. 2017. Т. 18. № 5. С. 1-30.
  7. S.L. Brunton, J.N. Kutz. Data-driven science and engineering: machine learning, dynamical systems, and control. – Cambridge University Press, 2020.
  8. Xu H. et al. Four-loop feedback control system with integrator design for hypersonic cruise missile. 21st AIAA international space planes and hypersonics technologies conference. 2017. С. 2112.

Қосымша файлдар

Қосымша файлдар
Әрекет
1. JATS XML
2. Fig. 1. Modeling of the dynamic of the nonlinear model and its qLPV form

Жүктеу (50KB)

© Shabashov A.A., Plotnikov A.A., 2024

Creative Commons License
Бұл мақала лицензия бойынша қолжетімді Creative Commons Attribution 4.0 International License.

СМИ зарегистрировано Федеральной службой по надзору в сфере связи, информационных технологий и массовых коммуникаций (Роскомнадзор).
Регистрационный номер и дата принятия решения о регистрации СМИ: серия ПИ № ФС 77 - 56417 от 11 декабря 2013.