САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
КАФЕДРА КОМПЬЮТЕРНЫХ ТЕХНОЛОГИЙ И СИСТЕМ
Колотов Михаил Евгеньевич
Выпускная квалификационная работа бакалавра
Синтез робастных регуляторов для управления
квадрокоптером
Направление 010400
Прикладная математика и информатика
Научный руководитель,
кандидат физ.-мат. наук,
доцент
Сотникова М. В.
Санкт-Петербург
2016
Содержание
Введение ........................................................................................................ 3
Глава 1. Постановка задачи ......................................................................... 5
1.1. Математическая модель динамики квадрокоптера ......................9
1.2. Формализованная постановка задачи .......................................... 15
Глава 2. Синтез регулятора ........................................................................ 17
2.1. Преобразование линейной системы ............................................17
2.2. Синтез закона управления ............................................................ 19
2.3. Синтез асимптотического наблюдателя ......................................23
Глава 3. Синтез робастных регуляторов для управления движением
квадрокоптера ............................................................................................. 25
3.1. Частотный анализ SISO-систем....................................................25
3.2. Анализ робастной устойчивости регуляторов ............................29
3.3. Имитационное моделирование динамики квадрокоптера .........33
Выводы ........................................................................................................ 37
Список литературы .....................................................................................38
2
Введение
Квадрокоптером называется летающий аппарат крестовидной формы с
моторами, находящимися на концах ребер аппарата (рис 1). Моторы
фиксированы. Контроллер (бортовой компьютер), аккумуляторы, датчики и
прочие периферийные устройства располагаются, по возможности, в центре
рамы. Квадрокоптер имеет 6 степеней свободы. Управление происходит либо
при помощи пульта дистанционного управления, где человек-оператор задает
ориентацию квадрокоптера и его высоту, либо
указанием маршрута,
которому квадрокоптер должен следовать самостоятельно. Также возможны
варианты самостоятельного функционирования, например, по видеоданным с
бортовой камеры.
Рис. 1. Квадрокоптер.
Квадрокперы используются человеком для перемещения любых
нетяжелых предметов (в среднем до трех килограмм) а также для фото- и
видеосъемки. Использование квадрокоптеров для съемки дает возможность
снимать с ракурсов, недоступных ранее или же доступных, но лишь после
специальной подготовки сцены. Использование же беспилотного
летательного аппарата (далее – БПЛА) в качестве курьеров позволяет
задействовать меньшее количество людей и снижает нагрузку на городскую
дорожную сеть.
3
Для того, чтобы летательный аппарат двигался согласно нашим
ожиданиям, необходимо спроектировать систему управления. Основное звено
любой системы управления это регулятор - устройство, которое следит за
состоянием объекта управления как системы и вырабатывает для неё
управляющие сигналы. Регуляторы следят за изменением некоторых
параметров объекта управления и реагируют на их изменение с помощью
некоторых алгоритмов управления. Естественно, что математическая модель
описывает поведение объекта не идеально точно. Требуется,
Данная работа посвещяна синтезу законов управления, способных
управлять перемещением квадрокоптера. На сложность регулятора
накладывается условие выполнимости в условиях реального времени на
бортовом компьютере (который имеет ограниченные вычислительные
ресурсы). Дополнительно требуется, чтобы регулятор, синтезированный для
одной модели квадрокоптера, функционировал и с другими моделями,
отличающимися от изначального по значениям параметров, например, массе,
размерам, мощностями моторов и т.п.. Для достижения поставленной цели
будет использоваться LQR-синтез с последующим приминением частотного
анализа робастности.
4
Глава 1. Постановка задачи
Построим математическую модель квадрокоптера в виде системы
обыкновенных дифференциальных уравнений, согласно [6]:
X́=f (X , X́ ,U )
Вектор управления U =(u1 ,u2 ,u3 , u4 )T
будет представлять собой усилие
каждого из четырех моторов. Вектор X
должен отражать все степени
свободы квадрокоптера, а именно, положение центра масс и ориентацию тела
в пространстве.
Как цель поставим задачу синтеза системы управления, которая
стабилизирует положение равновесия на желаемой высоте с желаемыми
углами ориентации квадрокоптера. При этом устойчивость должна
сохраняться при вариации параметров модели, которые тяжело точно
измерить в реальных условиях.
Рис 2. Строение квадрокоптера.
Структура квадрокоптера, угловые скорости винтов, подъемная сила и
крутящий момент моторов представлены на рис. 2. Положение центра масс
квадрокоптера в неподвижной связанной с землей системе координат
5
обозначим ξ . Для ориентации квадрокоптера в той же системе координат
введем углы Эйлера η . θ
крена φ
– угол вращения вокруг оси y
(тангаж). Угол
определяет вращение вокруг оси x . Угол рысканья, который
обозначим как ψ, определяет вращение вокруг оси z . Вектор q
–
конкатенация векторов ξ и η
[] []
x
ξ= y
z
,
φ
η= θ
ψ
[]
ξ
, q= η
За центр связанной системы координат возьмем центр масс
квадрокоптера. В этой системе линейные скорости обозначим как V B , а
угловые скорости обозначим за ν
[ ] []
υx, B
V B= υ y , B
υz , B
,
p
ν= q
r
Матрицу перехода из связанной системы в абсолютную обозначим за R
[
C ψ C θ C ψ S θ S φ−S ψ C θ C ψ S θ C φ +S ψ S φ
R= Sψ C θ Sψ Sθ S φ +C ψ C θ Sψ Sθ C φ +C ψ S φ
−Sθ
C θ Sφ
CθCφ
]
где Sψ =sin ( ψ ) ,C ψ =cos ( ψ ) . Так как матрица поворота R
−1
R =R
T
ортогональная, то
.
Матрицу преобразования угловых скоростей из абсолютной системы
координат в связную с квадрокоптером обозначим за
[
[][
1 S φ T θ Cφ T θ
0 Cφ
−S φ
W η=
Sψ
Cψ
0
Cθ
Cθ
ν=W η ή
]
][ ]
1 S φ T θ Cφ T θ
p
0 Cφ
−S φ φ́
q =
θ́
Sψ
Cψ
r
ψ́
0
Cθ
Cθ
где T θ=tan ( θ ) . Матрица W η обратима при θ ≠
6
(2 k−1)
φ,k ∈ Z .
2
Предположив, что квадрокоптер имеет симметричную относительно
осей x и y , структуру, получим диагональную матрицу тензора инерции,
состоящую лишь из главных центральных моментов инерции. Также будем
считать, что моменты инерции известны нам неточно.
[
]
I xx 0 0
I = 0 I yy 0 .
0 0 I zz
Угловая скорость винта i -го мотора, обозначенная за ωi ,
п орож дает п одъ ем н ую силу, направленную по о си вращения
соответствующего мотора f i .
Совместно с угловым ускорением, угловая
скорость также порождает крутящие моменты вокруг оси вращения винта
τM .
i
2
2
f i=k ωi , τ M =b ωi + I M ώ i
i
где k
– подъемный коэффициент и b
– коэффициент крутящего
момента моторов, а I M - момент инерции винта. Так как винт очень легок,
эффект ώi
обычно опускается. Будем считать, что величины k
и b
извыестны нам неточно. Объединяя подъемные силы четырех моторов,
получим тягу T , направленную вдоль оси аппликат в связанной с телом
системе координат ( z B ). Общий крутящий момент τ B
крутящих моментов, направленных по углам Эйлера
состоит из
связанной системы
координат τ φ , τ θ и τ ψ .
[]
4
4
0
2
T =∑ f i=k ∑ ωi ,T B= 0 ,
i=1
i=1
T
[][
2
2
]
lk (ω4−ω2 )
τφ
2
2
τ B = τ θ = lk (ω3−ω 1 ) ,
4
τψ
∑ τM
где l
i=1
i
– расстояние от центра масс до винта. Следовательно, увеличение
крена достигается увеличением скорости вращения 4-го мотора
(относительно 2-го), увеличение тангажа достигается увеличением скорости
7
вращения 3-го мотора (относительно 1-го), а изменение угла рысканья
достигается изменением мощностей винтов, вращающихся сонаправленно
(относительно скорости двух других, вращающихся в противоположную
сторону).
8
1.1 Математическая модель динамики квадрокоптера
Рассматривая квадрокоптер как твердое тело, можем описать его
динамику при помощи уравнений Ньютона-Эйлера. В связной системе
координат, сила, необходимая для ускорения массы m V́ B
и центробежная
сила ν ×(m V B ) равны гравитации R T G и общей тяге моторов T B
T
m V́ B +ν × ( m V B ) =R G+T B
В абсолютной системе координат центробежная скорость обнуляется, и,
следовательно, на ускорение квадрокоптера воздействует только
гравитационная сила, величина и направление тяги.
m ξ́=G+ R T B
[] [] [
x́
0 T C ψ Sθ C φ +Sψ Sφ
ý =−g 0 + S ψ Sθ C φ−C ψ Sφ
m
ź
1
Cθ Cφ
]
(1)
В связной системе координат, угловые ускорения инерции I ν́ ,
центростремительные силы ν ×(Iν)
и гироскопические силы δ
равны
внешнему моменту τ
I ν́+ ν ×(Iν)
+ δ =τ ,
([][ ] [] [] )
[][
][ ][ ]
ν́=I
−1
I xx p
p
p
0
− q × I yy q − I r q × 0 ω δ + τ
r
r
1
I zz r
,
( I yy −I zz ) qr /I xx
q/ I xx
τ φ /I xx
ṕ
=
−
I
ω
+
q́
( I zz− I xx ) pr /I yy r − p /I yy δ τ θ /I yy
ŕ
0
τ ψ / I zz
( I xx −I yy ) qp/ I zz
,
Где ω δ =ω1 −ω2 +ω3 −ω4 . Угловые ускорения в абсолютной системе координат
это производные по времени угловых скоростей из связной системы,
преобразованных матрицей трансформации W −1
.
η
[
]
0
φ́ C φ T θ + θ́ S φ /C 2θ
−φ́ S φ C θ + θ́ C φ /C 2θ
d
d −1
−1
−1
−1
ή= ( W η ν ) = W η ν+W η ν́= 0
ν +W η ν́
−φ́ S φ
−φ́ C φ
dt
dt
0 φ́ C φ /C θ + φ́ Sφ T θ /C θ −φ́ S φ /C θ + θ́ C φ T θ /C θ
9
Уравнения Лагранжа второго рода
Лагранжиан L суть сумма вращательной E rot
и переносной E trans
энергий минус потенциальная энергия E pot
L ( q , q́ )= Etrans +E rot −E pot =
m T 1 T
ξ́ ξ́ + ν Iν−mgz
2
2
Уравнения Лагранжа второго рода имеют вид
[] ( )
d ∂L ∂L
f
=
−
τ dt ∂ q́ ∂ q
Линейные и угловые компоненты не влияют друг на друга,
следовательно, они могут быть рассмотрены раздельно. Линейная внешняя
сила суть общая тяга моторов. Линейное уравнение Лагранжа второго рода:
[]
0
f =R T B =m ξ́ +mg 0
1
Что эквивалентно уравнениям (1). Матрица Якоби преобразования J ( η ) от
ν
к ή :
T
J ( η )=W η I W η ,
−I xx S θ
I
I xx
0
¿
I
2
2
0
I yy C φ + I zz S φ
(¿ ¿ yy−I zz )C φ Sφ C θ
(¿ ¿ yy−I zz )C φ Sφ C θ
¿−I xx Sθ
¿
2
2 2
2 2
¿ [ I xx Sθ + I yy C θ S φ + I zz C φ C θ ¿ ]
Следовательно, вращательная энергия E rot
может быть выражена в
абсолютной системе координат как
1 T
1 T
E rot = ν Iν= ή J ή
2
2
Внешняя угловая сила – суть моменты моторов. Уравнения Лагранжа второго
рода для углов это
ή
(¿¿ T J ή)=J ή+C(η , ή) ή
d
1 ∂
τ=τ B=J ή+ ( J ) ή−
¿
dt
2 ∂η
10
(2)
Где матрица C (η , ή) имеет вид
[
C 11 C 12 C 13
C ( η , ή )= C 21 C 22 C 23
C 31 C 32 C 33
]
C 11 =0
C 12= ( I yy − I zz ) ( θ́ C φ Sφ +ψ́ S2φ C θ )+ ( I zz −I yy ) ψ́ C θ C 2φ −I xx ψ́ C θ
2
C 13 =( I zz −I yy ) ψ́ C φ C θ Sφ
C 21= ( I zz −I yy )
C 22= ( I zz −I yy ) ( θ́ C φ Sφ +ψ́ S❑φ C θ )+ ( I yy −I zz ) ψ́ C θ C 2φ +I xx ψ́ C θ
2
C 23 =−I xx ψ́ C θ S θ + I yy ψ́ Sφ C θ Sθ
2
C 31= ( I yy −I zz ) ψ́ C φ C θ Sφ −I xx C θ θ́
C 32= ( I zz −I yy ) ( θ́ C φ Sφ Sθ + φ́ S 2φ C θ )+ ( I yy −I zz ) φ́ C θ C 2φ +I xx ψ́ C θ S θ −I yy ψ́ S2φ C θ Sθ −I zz ψ́ C θ C 2φ Sθ
2
2
2
C 33 =( I yy −I zz ) φ́ C φ Sφ C θ −I yy θ́ Sφ C θ Sθ −I zz θ́ C θ S θ C φ + I xx θ́ C θ Sθ
Выражение (2) приводит к дифференциальным уравнениям для угловых
скоростей
−1
ή=J (τ B−C ( η , ή ) ή)
Аэродинамический эффект
Данная модель получена с
упрощением сложных динамических
взаимодействий. Для того, чтобы добиться более реалистичного поведения
квадрокоптера, учтем силу сопротивления воздуха.
[] [] [
] [
A
x́
0 T C ψ Sθ C φ +Sψ Sφ
1 x
0
ý =−g 0 + S ψ Sθ C φ−C ψ Sφ −
m
m
ź
1
Cθ Cφ
0
Гд е A x , A y
и Az
0
Ay
0
][ ]
0 x́
0 ý
A z ź
(3)
коэффициенты сопротивления воздуха в направлении
соответствующих осей абсолютной системы координат. Будем считать их
сложноизмеримыми.
11
Линеаризация
С
точки зрения анализа устойчивости и синтеза управления нас
б о л ь ш е и н т е р е с у е т л и н е й н о е п р и бл и ж е н и е д а н н о й с и с т е м ы
дифференциальных уравнений. Для того, чтобы при линеаризации
управление полностью не ушло из системы, возьмем за него квадрат угловой
скорости винта - потеря знака нас не беспокоит, поскольку винты могут
вращаться только в одну сторону, а их направление вращения уже учтено в
уравнениях.
[]
u1
u
2
u= 2 ,ui=ωi
u3
u4
За положение равновесия возьмем
[ ] [ ][ ] [ ][ ] [ ]
[ ] [ ][ ] [ ][ ] [ ]
x́
x0
0 x́
0 x
ý = 0 , ý = 0 , y = y 0 ,
ź
0 ź
0 z
z0
φ́
0 φ́
0 φ
0
=
,
=
,
=
0 θ́
0 θ
0 ,
θ́
0 ψ́
0 ψ
ψ0
ψ́
[]
u0
u
u= 0 .
u0
u0
Гд е u0
– квадрат угловой скорости одного мотора такой, что тяга,
создавая четырьмя такими значениями будет полностью компенсировать силу
тяжести. Иначе говоря, управление, при котором дрон неподвижно висит в
воздухе. Значение находится из уравнения динамики по
условии, что ź=0, ź ≈ 0,θ=φ=0 :
ź=−g +
k (u0 +u0 +u0+u0 )
Az
C θ C φ − ź
m
m
относительно переменной u0 .
12
z
(3), при
u0=
mg
4k
Для линеаризации системы разложим правую часть системы в отклонениях
от положения равновесия в ряд Тейлора как функцию нескольких
переменных в окрестности положения равновесия, приняв за управление
квадрат угловой скорости винтов, и отбросим нелинейные слагаемые. В
результате получим следующую математическую модель:
{
Ax
x́
m
A
ý=−gφ− y ý
m
x́=g θ−
A
k
( u1+u2+u3 +u 4 )− z ź
m
m
kl
φ́= (u 4−u2 )
Ix
kl
θ́= (u3 −u1 )
Iy
b
ψ́= (u1 −u2 +u3−u4 )
Iz
ź=
13
(4)
1.2 Формализованная постановка задачи
Как цель поставим задачу синтеза системы управления с законами вида
U =K (s) y ,
где передаточная матрица регулятора K (s)
не изменяется со временем.
Данная система должна стабилизировать положение равновесия
[]
x des
y des
z
X des= des
φdes
θ des
ψ des
.
При этом устойчивость должна сохраняться при отклонении значений
параметров k ,b , I , A
от номинальных на 30%. Если данное требование
невыполнимо, то достичь максимально-возможного диапазона для данных
параметров.
Количество русскоязычных материалов по данной тематике весьма
скудно. И в абсолютном большинстве случаев работы выполнены в виде
небольших статей. В [1] рассмотрена возможность синтеза системы
управления на базе адаптивных ПИД-регуляторов. Цель работы оптимизация переходного процесса и, в результате, повышение устойчивости
беспилотного летательного аппарата.
В статье [2] анализируют вопросы математического моделирования
динамики квадрокоптера, по результатам анализа строится нейросетевой
регулятор. Плюсом данного регулятора является его возможность адекватно
реагировать на аэродинамические эффекты, которые весьма трудно полно
описать математически. Однако, качество данного типа регуляторов сильно
зависит от обучающих данных. Примененный в [3] метод основывается на
линеаризации уравнений обратной связью с применением полиномов
Баттерворта. Полученная система позволяет перемещать аппарат в заданную
точку и поворачивать его на заданный угол вокруг вертикальной оси. Во всех
14
перечисленных статьях отсутствует анализ робастности синтезированных
регуляторов.
В англоязычной литературе множество статей с различными методами
синтеза. По методам, близким к тем, которые используются в этой работе,
есть статья [4]. В ней строится LQR-регулятор, однако, с допущениями, что
вращательные моменты и тяга моров могут управляться независимо. Но это
не соответствует действительности. В статье [5] авторы затрагивают вопрос
робастности и выставляют точно такие же границы неопределенности
параметров в 30%, однако используется нелинейное управление типа sliding и
backstepping.
15
Глава 2. Синтез регулятора
2.1. Преобразование линейной системы
Полученная линейная система (4) имеет несколько входов и несколько
выходов, поэтому её анализ будет достаточно громоздким. Этого можно
избежать, если попытаться разбить её на несколько систем с одним входом и
одним выходом. Для этого уравнения должны быть независимыми. Этого
можно достичь, если допустить следующие ограничения на управление:
величина ui будет складываться из двух составляющих:
{
u1=uz +u Δθ
u2 =u z +u Δφ
u3 =u z −u Δθ
u 4 =u z −u Δφ
(5)
Иначе говоря, первая составляющая отвечает за сумму всех управлений и
будет влиять на динамику z , тогда как вторая составляющая будет влиять
на динамику соответствующего мотору винта.
Тогда система (4) будет выглядеть:
{
Ax
x́
m
A
ý=−gφ− y ý
m
kl
φ́=−2 uΔφ
I xx
¿
kl
θ́=−2
u
I yy Δθ
b
ψ́= 0
I zz
Az
k
ź=4 u z − ź
m
m
x́=gθ−
(6)
Используя подобное допущение мы лишаемся возможности управлять углом
ψ́ , однако, это можно исправить, если для
ui
ввести аналогичную
специальную замену и синтезировать отдельный регулятор. Потребуем,
16
ч т о б ы θ́= φ́=0
следовательно,
{
u4 −u2 =0
. Удовлетворяющая этим
u3 −u1 =0
требованиям замена:
{
u1=u z +u Δψ
u2=uz −u Δψ
u3=u z +u Δψ
u 4=u z −u Δψ
Тогда система (4) будет выглядеть:
{
kl
0
I xx
kl
θ́=−2 0
I yy
b
ψ́=4 u Δψ
I zz
Az
k
ź=4 uz − ź
m
m
φ́=−2
(7)
Рассмотрим систему (6). Очевидно, что она состоит из независимых
уравнений. Тогда разобьем её на системы в пространстве состояний:
{
ź 1= z 2
− Az
k
ź 2 =
z 2 +4 u z
m
m
y=z 1
(8)
θ́1=θ2
kl
θ́2 =−2
u
I yy Δθ
y=θ 1
(9)
{
17
{
x́ 1 =x 2
A
x́ 2 =g θ 1− x x1
m
θ́1 =θ2
kl
θ́2 =−2
u
I yy Δθ
y=x 1
(11)
Система для φ аналогична системе (9) и отличается лишь осью момента
инерции, который, в силу структуры квадрокоптера, одинаковый. Так же,
система для y аналогична системе (11) и отличается лишь тем, что вместо
угла θ , в ней используется угол φ . Поэтому в дальнейшем будем
рассматривать лишь системы (9) и (11). Аналогичными рассуждениями из
системы (7) получим уравнения динамики ψ :
{
ψ́ 1 =ψ 2
b
ψ́ 2 =4 uΔψ
I zz
y=ψ 1
(12)
2.2 Синтез законов управления
Для построения законов управления для систем (8)-(11) воспользуемся
LQR-синтезом. LQR-синтез подразумевает собой модальный синтез с
условием, что корни системы располагаются так, чтобы минимизировать
функционал (12), отвечающий за оптимальное энергопотребление при
желаемом быстродействии:
Пусть имеется SS-модель
{
´ Bu
x=Ax+
y=Cx
с линейной обратной связью
u=−Kx
Синтез управления заключается в решении задачи оптимизации (12)
относительно SS-системы.
18
x
(¿ ¿ T Qx +u Ru)dt → min
T
(12)
∞
J ( u )=∫ ¿
0
Матрицы Q и R
задаются, обычно, диагональными, таким образом, чтобы
регулятор удовлетворял желаемой динамике. Чем больше значения
коэффициентов Q относительно коэффициентов R, тем интенсивнее будет
управляющий сигнал.
Матрица K , в таком случае, имеет вид K =−R−1 BT P ,
где P находится
из матричного уравнения Риккати
T
1
T
A P+ PA−PB R B P+Q=0
Важно отметить, что для возможности синтезировать LQR-регулятор
необходимо измерять весь вектор состояния, а SISO-система по своей сути не
обеспечивает полноты набора измеряемых велечин. Для удовлетворения
этого условия построим и объединим с регулятором асимптотический
наблюдатель. Система наблюдателя
{
ź= ( A−GC ) z+Bu+Gy
y obs= z
Входным сигналом является вектор измерений
(13)
y , вектор состояния это
оценка вектора состояния системы, для которой строится наблюдатель.
Коэффициенты вектора G выбираются таким образом, чтобы обеспечивать
устойчивость системы. Закон управления будет иметь вид
u=−Kz
(14)
Замкнув (13) обратной связью (14), получим систему:
{
ź= ( A−GC−BK ) z+Gy
u=−Kz
(15)
Выпишем матрицы A , B и C для систем (8)-(11) соответственно:
[ ] []
0
A=
1
0
− A z , B= 4 k , C=[ 1 0 ]
0
m
m
19
(16)
[ ]
0
0 1
A=
, B= −2 kl ,C=[ 1 0 ]
0 0
I yy
[ ]
(17)
[ ][]
0
1
−A x
A= 0
m
0
0
0
0
0 0
0
0
g 0 , B=
0 , C=[ 1 0 0 0 ]
−2 kl
0 1
I yy
0 0
(18)
[]
0
0 1
A=
, B= 4 b ,C=[ 1 0 ]
0 0
I zz
[ ]
(19)
Решим, теперь, для каждой системы задачу (12) при помощи пакета
прикладных программ Matlab. Взяв за номинальные значения параметров
БПЛА
значения:
−6
−3
−8
m=2 ,k=3∗10 , I xx =I yy =5∗10 , I zz =0.01,l=0.4, A x = A y = A z =0.25, b=7∗10
K z=
K φ= K θ =
[
[
.
]
3
6
∗10
0.97938
]
−3
6
∗10
−0.1118
[ ]
−0.9487
−0.7564
6
K x =K y =
∗10
−3.2707
−0.1167
[
Kψ=
]
3
6
∗10
0.45847
Таким образом, получим три регулятора. Первый для управления
высотой, креном и тангажом, второй для управления координатами центра
масс, третий – для управления рысканьем.
[]
z−z des
ź− ź des
uz
φ−φ des
=¿
uΔφ =−K
φ́−φ́ des
u Δθ
θ−θdes
θ́−θ́des
[]
20
[]
z− z des
ź− ź des
3 0.97938 0
0
0
0
φ−φdes
6
¿−10 0
0
−3 −1.1180 0
0
φ́− φ́des
0
0
0
0
−3 −1.1180
θ−θ des
θ́−θ́ des
[
]
[]
x− xdes
x́− x́des
y− y des
ý− ý des
uz
φ−φdes
=¿
uΔφ =−K
φ́− φ́des
u Δθ
θ−θ des
θ́−θ́ des
z− z des
ź−ź des
[]
[
0
0
0
0
0
0
0
0 K z1 K z2
¿−10 6 0
0 K y 1 K y 2 K y3 K y 4 0
0
0
0
K x1 K x2 0
0
0
0 K x3 K x 4 0
0
[ ]
[ ]
z−z des
z− z des
uz
ź− ź des
ź− ź des
0
=−K
=−106 3 0.97938 0
0
0
3 0.45847 ψ−ψ des
ψ−ψ des
u Δψ
ψ́−ψ́ des
ψ́− ψ́ des
[ ]
[
]
21
]
[]
x− xdes
x́− x́des
y− y des
ý− ý des
φ−φdes
φ́− φ́des
θ−θ des
θ́−θ́ des
z− z des
ź−ź des
2.3 Синтез асимптотического наблюдателя для
квадрокоптера
Перейдем теперь к пострению асимптотичских наблюдателей для
систем (8)-(11). Характеристические полиномы систем-наблюдатлей:
2
Δφ́ , θ́ , ψ́ ( s )=s +2 g1 s+ g1
s1 = √ g1 −g1 −g1 , s 2 =−√ g1 −g1 −g1
2 g m+ A z
A g +g m
Δ ź ( s ) =s2 + 1
s+ z 1 1
m
m
2
2
−( A z +2 g1 m− √ A z +4 g1 m −4 g1 m )
s 1=
,
2m
2
2 2
2
−( A z +2 g1 m+√ A z +4 g1 m −4 g1 m )
s2=
2m
g
Δx́ , ´y ( s )=s 4 + ( g1 +0.125 ) s3 + 1 + g2 s2 +9.81 g 3 s+9.81 g 4
8
2
(
2
2
2
)
Желая, чтобы вещественная часть корней характеричстических полиномов
была отрицательной, найдем соответствующие значения gi .
g z =gφ ,θ , ψ =4
[] [ ]
g1
−0.1017
g
3 −0.3918
x , y : 2 =10
−7.8893
g3
−6.7187
g4
По результатам моделирования наблюдателя в пакете прикладных программ
Matlab можно сделать вывод, что полученные наблюдатели полностью
удовлетворяют желаемой точности. На рис. 3 представлена разность
наблюдений и действительного значения векторая состояния системы (8).
Порядок погрешности 10−18
и ее вид говорит о том, что это больше
вычислительная ошибка, чем теоретическая.
22
t
Рис. 3. Разность действительного и наблюдаемого векторов состояния
23
Глава 3. Синтез робастных регуляторов для
управления движением квадрокоптера
3.1 Частотный анализ SISO-систем
Проведем частотный анализ систем согласно [7]. Построим АЧХ для
систем (8)-(11) с номинальными значениями параметров и при изменении
параметров на 30%.
|P n( j ω)|
.
ω
Рис. 4. АЧХ системы (8) при различных параметрах A z , k
24
|P n( j ω)|
ω
Рис. 5. АЧХ системы (11) при различных параметрах I zz ,b .
|P n( j ω)|
ω
Рис. 6. АЧХ системы (9) при различных параметрах k , I xx , yy .
25
|P n( j ω)|
ω
Рис. 7. АЧХ системы (10) при различных параметрах k , I xx , A x .
Далее построим АЧХ относительного возмущения данных систем.
|Δ( j ω)|
ω
Рис. 8. АЧХ относительного возмущения системы (9).
26
|Δ( j ω)|
ω
Р
ис. 9. АЧХ относительного возмущения системы (8)
|Δ( j ω)|
ω
Рис. 10. АЧХ относительного возмущения системы (11)
27
| Δ( j ω)|
ω
Рис. 11. АЧХ относительного возмущения системы (10)
28
3.2 Анализ робастной устойчивости регуляторов
Для анализа робастности полученого регулятора воспользуемся
частотным подходом, изложенным в [7]. Найдем передаточные функции
систем управления вида (15), для рассматриваемых SISO-систем (8)-(11).
−12 s−1.5
s +10 s+42.01
12 s−4997
6
W φ ,θ ( s )=10 2
s −466.1 s−3320
−12 s−175
W ψ ( s )=10 6 2
s −31.42 s−193.7
3
2
−2.7 s −8.4 s −12.6 s−6.2
10
W x , y ( s )=10 4
s −45.5 s 3−4525 s 2−255900 s−5395000
W z ( s )=10
6
2
Согласно [7], достаточным условием робастности регулятора является
выполнение неравенства
|[ P ( jω )−P
( jω ) ] P n ( jω )|<
−1
n
1
,∀ ω ∈ ¿
|T ( jω)|
Г д е T ( s ) =W ( s ) [1+P n ( s ) W (s)]−1 P n ( s ) ,
функция, а P ( s )−¿
Pn ( s )
– номинальная передаточная
передаточная функция объекта при измененных
параметрах. Правая часть неравенства трактуется как частотная граница
устойчивости (на изображениях красная кривая). Левая часть неравенства –
относительное возмущение передаточной функции объекта (на изображениях
синяя кривая).
29
|P n( j ω)|
ω
Рис. 12. Граница робастности регулятора по крену и тангажу
|P n( j ω)|
ω
Рис. 13. Граница робастности регулятора высоте
30
|P n( j ω)|
ω
Рис. 14. Граница робастности регулятора по рысканью
|P n( j ω)|
ω
Рис. 15. Граница робастности регулятора по x , y
Как видно по изображениям, синтезированные регуляторы обладают
достаточным запасом робастной устойчивости.
31
3.3 Имитационное моделирование динамики
квадрокоптера
При помощи пакета прикладных программ Matlab (и Simulink)
проведем моделирование нелинейной модели квадрокоптера, полученной в
главе 1, с использованием синтезированных регуляторов. Сравним
результаты при изменении сложноизмеримых параметров модели на 30%.
Желаемыми значениями будем считать x=2 , z=5 ,ψ=60° ,θ=φ=15° .
Рис. 16. Simulink-модель нелинейной модели квадрокоптера.
Рис. 17. Simulink-модель нелинейной модели квадрокоптера,
замкнутая управлением
33
z
t
Рис. 18. Результаты моделирования нелинейных систем. Высота
БПЛА при номинальных параметрах и отклонении значений параметров
на 30%.
φ
ψ
34
Рис. 19-20.
21. Результаты
Результаты
моделирования
моделирования
нелинейных
нелинейных
систем.
систем.
Координаты
Крен(тангаж) xи рысканье
и y БПЛА
БПЛА
при
при
номинальных
номинальных
параметрах
параметрах
и и
t
отклонении
отклонениизначений
значенийпараметров
параметровнана30%.
30%.
t
x,y
t
Рис. 21. Результаты моделирования нелинейных систем. Координаты x , y
БПЛА при номинальных параметрах и отклонении значений параметров на 30%.
Код алгоритма системы управления представлен ниже на лис. 1.
function u = fcn(inp)
xdes=(inp(1));
ydes=(inp(2));
zdes=(inp(3));
fdes=(inp(4));
tdes=(inp(5));
pdes=(inp(6));
X=inp(7:18);
kp1=-3;
kp2=-0.45847;
kz1=-3;
kz2=-0.97938;
kx1=0.9487;
kx2=0.7564;
kx3=3.2707;
kx4=0.1167;
X1=[X(1) X(4) X(2) X(5) X(7) X(10) X(8) X(11) X(3) X(6)]';
X1des=[xdes 0 ydes 0 fdes 0 tdes 0 zdes 0]';
X2=[X(3) X(6) X(9) X(12)]';
X2des=[zdes 0 pdes 0]';
K1=10e+6*[0 0 0 0 0 0 0 0 kz1 kz2;...
0 0 kx1 kx2 kx3 kx4 0 0 0 0;...
kx1 kx2 0 0 0 0 kx3 kx4 0 0];
K2=10e+6*[kz1 kz2 0 0; 0 0 kp1 kp2];
35
U=zeros(4,1);
U1=zeros(3,1);
u2=zeros(2,1);
D=sum(abs(X1(1:8)-X1des(1:8)));
if D>0.01
U1=-K1*(X1des-X1);
U(1)=U1(1)+U1(3);
U(2)=U1(1)+U1(2);
U(3)=U1(1)-U1(3);
U(4)=U1(1)-U1(2);
else
U2=-K2*(X2des-X2);
U(1)=U2(1)+U2(2);
U(2)=U2(1)-U2(2);
U(3)=U2(1)+U2(2);
U(4)=U2(1)-U2(2);
end
u = U;
Лис. 1. Программный код, реализующий алгоритм управления
написанный на зыке Matlab.
36
Выводы
1. Математическая модель квадрокоптера была преобразована для
упрощения работы (линеаризована и разбита на SISO-системы).
2. С помощью технологии LQG-синтеза, была получена система
управления, позволяющая занимать квадрокоптеру желаемую
позицию и ориентацию.
3. Метод частотного анализа робастности [7] гарантирует, что данная
система управления имеет запас устойчивости по параметрам
в 30%.
4. Имитационное моделирование подкрепило вышеизложенные
k ,b , I , A
результаты применительно к нелинейной модели.
Если бы в ходе анализа [7] выяснилось, что данный регулятор не
удовлетворяет поставленным треованиям, то далее, посредством вариации
коэффициентов регулятора (и наблюдателей), нужно было бы добиться
максимально возможного запаса робастной устойчивости.
Таким образом, поставленная задача была полностью решена, однако,
хотелось бы также синтезировать робастную систему управления для
следования вдоль траектории. Для достижения этой цели планируется
изучить подходы к построению управления backstepping и sliding с точки
зрения робастности.
37
Список литературы
1. Юшкин Д. А., Евдокимов С. А. Разработка адаптивного нечеткого ПИДрегулятора системы автоматического управления и стабилизации
мультироторного БПЛА типа квадрокоптер // Актуальные проблемы
современной техники, науки и образования, 2015. Т. 2, № 1. С. 194-198.
2. Павловский В.Е., Яцун С.Ф., Емельянова О.В., Савицкий А.В.
Моделирование и исследование процессов управления квадрокоптером //
Робототехника и техническая кибернетика, 2014. № 4(5). С. 49-57.
3. Гурьянов А. Е. Моделирование управления квадрокоптером // Инженерный
вестник, 2014 №8, С. 4.
4. Cowling, I.D., Whidborne, J.F. and Cooke, A.K. Optimal trajectory planning and
LQR // Proc. UKACC Int. Conf. Control 2006 (ICC2006), September 2006,
Glasgow, UK.
5. Guilherme V. R., Manuel G. O., Francisco R. R., Robust Nonlinear Control for
Path Tracking of a Quad-Rotor Helicopter, 2015 V. 17, I. 1, P. 142-156.
6. Luukkonen T., Modelling and control of quadcopter, 2011, P.2-6.
7. Веремей Е. И., Линейные системы с обратной связью: Учебное пособие.
2013, С. 372-384.
38
Отзывы:
Авторизуйтесь, чтобы оставить отзыв