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

Министерство образования  и науки РФ

_________

 

Федеральное агентство по образованию

 

Московский авиационный институт

(государственный технический университет)

_______________________________________________________________________

 

 

 

Факультет №3

Кафедра 301

Группа 03 – 401

 

 

 

 

 

Курсовая работа

 

«Исследование динамики и синтез цифровых алгоритмов управления

боковым движением  летательного аппарата».

 

 

 

 

 

 

 

 

 

 

 

 

Москва

2008г.

 

Содержание:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Задание.

Заданы:

1. уравнения движения объекта управления – самолета в боковом движении.
2. структурная схема системы управления
3. уравнения исполнительных устройств
4. технические требования к системе

 

 

Структурная схема системы  управления боковым движением.

 

 

Уравнения ЛА:

 

Требуется:

Разработать алгоритмы  анализа и синтеза. Провести с  их помощью исследования динамики и  выбор параметров САУ.

 

Перечень задач:

1.  Исследование динамики  объекта управления и выбор периода квантования T₀:

а) провести анализ динамики многосвязного объекта управления:

- построить переходные  процессы объекта управления  при ступенчатом отклонении элеронов  и при ступенчатом отклонении  руля направления;

- определить корни  характеристического уравнения объекта;

- построить частотную  характеристику объекта управления  по координате w_x;

- сделать выводы о  динамических свойствах системы;

б) Выбрать параметр T0 – период дискретизации для реализации цифровых алгоритмов управления:

 

2. Произвести синтез передаточных чисел в канале руля направления исходя из требований к затуханию колебаний самолета по рысканию.

 

а) получить Z-передаточные функции самолета от δ_н к ω_у;от δ_н к β;

б) для выбранного периода  дискретизации Т₀ вычислить z-передаточную функцию соединения ключ-фиксатор-объект;

в) составить характеристическое уравнение дискретной системы в  канале руля направления;

г) построить область  устойчивости в плоскости коэффициентов  К11 и К12 с использованием метода стандартных z-передаточных функций;

д) выбрать оптимальные  значения К11 и К12;

е) построить переходные процессы по координатам ω_у и β.

 

3. Записать уравнение  состояния дискретной системы  с учетом выбранной обратной  связи в канале руля направлений. 

 

4. Произвести синтез  ЦАУ в канале элеронов частотными методами.

а) построить псевдочастотные  характеристики;

б) выбрать параметры  цифрового регулятора;

в) построить переходные процессы по w_х.

 

5. Выбрать коэффициент  электрической связи Кх от  ручки летчика из условия обеспечения  требуемой эффективности управления самолета.

 

6. Записать уравнение  дискретного алгоритма оптимальной фильтрации для вычисления неизменяемой координаты объекта. Составить программу моделирования и исследовать динамику оптимального фильтра.

 

 Требования к системе:

 

1. Не менее чем 2-х кратные запасы устойчивости по всем коэффициентам обратных связей.
2. Затухание короткопериодических колебаний по ω_у и β не менее чем в 10 раз за период.
3. Собственная частота колебаний не менее 3 рад/с.
4. Время переходного процесса при управлении ω_х не более 1с при монотонном характере переходного процесса (перерегулирование не более 5%).

 

 

1. Исследование динамики  объекта управления.

 

а)Найдем матрицу передаточных функций системы.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

б)Построим переходные процессы в системе при ступенчатом отклонении рулей направления и элеронов

 

Переходный процесс по ω_y.

—— - переходный процесс  по ω_y при отклонении рулей направления.

– – – - переходный процесс по ω_y при отклонении элеронов.

 

Переходный процесс по β.

—— - переходный процесс  по β при отклонении рулей направления.

– – – - переходный процесс по β при отклонении элеронов.

Переходный процесс по ω_x.

—— - переходный процесс  по ω_x при отклонении рулей направления.

– – – - переходный процесс по ω_x при отклонении элеронов.

 

в)Найдем корни характеристического  многочлена системы:

Построим частотную  характеристику объекта:

A=[-0.6 -5.71 -0.04;1 -0.22 0.065;-0.7 -24 -2.5];

B=[-1.07 0.12;-0.005 0;-2 -6.5];

C=eye(3);

D=zeros(3,2);

WS=ss(A,B,C,D);

WS1=tf(WS(3,2))

margin(WS1)

 

 

Сделаем выводы по полученным переходным процессам:

1. Объект малодемпфированный, с достаточно низкой частотой  собственных колебаний (ω^*≈2.546 рад/с).

2. Переходные процессы по ω_x T_пп≈4÷5с, что не удовлетворяет требованиям, наложенным на систему.

Это требует демпфирования  системы, а также корректировки  системы, для улучшения времени  переходного процесса, что достигается  изменением корней характеристического  уравнения за счет коэффициентов  К11, К12, К23.

Для расчёта идеальных ОС зададим желаемый характеристический многочлен

       

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

Мнимая часть комплексного корня.

Im(s) = ω = 4.5

Вещественная часть  комплексного корня.

Время переходного процесса Т_пп=1с. Тогда действительная часть комплексного корня будет равна:

Re(s) ≤ – 2.303.

s_1,2=  – 2.303±4.5i

Возьмем действительный корень равный Re(s).

s₃= – 2.303

Найдем характеристический полином системы.

 

 

Определим коэффициенты К₁₁, К₁₂, К₂₃.

 

 

 

 

 

 

 

 

 

В результате получаем следующие  значения коэффициентов:

K₁₁=3.1502

K₁₂=14.13

K₂₃=0.016563

Построим переходные процессы с полученными коэффициентами.

 

 

Переходный процесс по ω_y.

—— - переходный процесс по ω_y при отклонении рулей направления.

– – – - переходный процесс по ω_y при отклонении элеронов.

 

Переходный процесс по β.

—— - переходный процесс по β при отклонении рулей направления.

– – – - переходный процесс по β при отклонении элеронов.

Переходный процесс по ω_x.

—— - переходный процесс  по ω_x при отклонении рулей направления.

– – – - переходный процесс по ω_x при отклонении элеронов.

 

Определение периода  квантования T₀

 

Для выбора периода дискретности Т₀ построим частотные характеристики объекта, замкнутого по ω_У и β, но разомкнутого по ω_Х. Также необходимо учесть динамику привода элеронов.

Матрица коэффициентов  обратных связей:

K=[3.1502 14.13 0;0 0 0.016563];

Матрица коэффициентов  обратных связей, разомкнутой по ω_Х системы.

K1=[3.1502 14.13 0;0 0 0];

A1=A+B*K1;

WS2=ss(A1,B,C,D);

W2=tf(WS2(3,2))

Передаточная функция от δэ к ω_X

   -6.5 s^2 - 28.54 s - 147.2

---------------------------------

s^3 + 6.761 s^2 + 35.43 s + 56.01

Wpre=tf([1],[0.11 1])

Передаточная функция привода

   1

----------

0.11 s + 1

Wwxeu=W2*Wpre

Передаточная функция от uэ к ω_X

          -6.5 s^2 - 28.54 s - 147.2

--------------------------------------------------

0.11 s^4 + 1.744 s^3 + 10.66 s^2 + 41.59 s + 56.01

margin(Wwxeu)

 

 

 

 

 

 

 

 

 

 

 

Частотная характеристика системы, замкнутой по каналу направления  и разомкнутой по каналу элеронов.

 

Выберем параметр T₀ из условия обеспечения запаса по фазе замкнутой системы не менее 60⁰:

На частоте ω_с = 5.08 рад/с запас по фазе составляет 88.5⁰

Фиксатор ухудшит нашу систему на ∆φ= ω_сT₀/2 рад

 

 Выберем Т₀=0.2 сек.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.Синтез передаточных  чисел в канале руля направления

 

Получение Z – передаточных функций ОУ по ω_У и β:

Составим уравнения  упрощенной системы, исключив из системы  ω_Х и δ_Э:

 

Aup=[-0.6 -5.71; 1 -0.26]

Bup=[-1.07;-0.005]

Cup=eye(2)

Dup=zeros(2,1)

T0=0.3;

Sup=ss(Aup,Bup,Cup,Dup)

 

Получение Z-передаточных функций по и :

Wup=tf(Sup)

SDup=c2d(Sup,T0,'zoh')

[F,D1,H,G]=ssdata(SDup)

Wz=tf(SDup)

step(Wz)

Передаточная функция по

       -0.1936 z + 0.1846

#1:  ---------------------

      z^2 - 1.631 z + 0.842

Передаточная функция по

 

      -0.02078 z - 0.01789

#2:  ---------------------

      z^2 - 1.631 z + 0.842     

 

 Составление характеристического уравнения дискретной системы в канале руля направления:

 

syms K11 K12 z w s

Kup=[K11 K12];

f=F+D1*Kup

f1=vpa(f,4)

f1=[.7853-.1936*K11,   -1.009-.1936*K12;.1766-.2078e-1*K11, .8453-.2078e-1*K12]

mf=det(z*eye(2)-f1)

h=vpa(collect(mf,z),4)

 

h=z^2+(-1.631+.1936*K11+.2078e-1*K12)*z+.1787e-1*K12-.1846*K11+.8420

Построение  области устойчивости дискретной системы  в плоскости параметров К11 и К12:

Для этого необходимо перейти из Z – области в w – область

 

z=(1+w)/(1-w);

zw=subs(h,z)

zw1=vpa(zw,4)

zw2=simplify(zw1)

zw2=-.1e-4*(-21100.-31600.*w-347300.*w^2-900.*K11+37820.*K11*w^2-3865.*K12+291.*K12*w^2+3574.*K12*w-36920.*K11*w)/(-1.+w)^2

zw5=vpa(collect(zw2,w),4)

 

K11=-2:.01:10;

Составим уравнения  границ области устойчивости:

K121=.3473e6/291.-.3782e5/291.*K11;

K122=.3692e5/3574.*K11+.3160e5/3574.;

K123=-900/3865.*K11-.2110e5/3865.;

plot(K11,K121,'r',K11,K122,'b',K11,K123,'g');grid;zoom yon;xlabel('K11');ylabel('K12')

 

Область устойчивости в  плоскости параметров К11 и К12:

 

 

 

 

 

 

 

 

Выбор коэффициентов  обратных связей в канале руля направления  К11 и К12

 

Зададим желаемый характеристический полином и выберем из области устойчивости оптимальные значения К11 и К12

 

Желаемый характеристический полином:

 

roots([1 2*0.707*4.5 4.5*4.5])

z1=exp(T0*(-3.1815 + 3.1825i))

z2=exp(T0*(-3.1815 - 3.1825i))

syms z

yu=expand((z-z1)*(z-z2))

vpa(collect(yu),3)

 

yu=z^2-.851*z+.280

Приравняем коэффициенты при одинаковых степенях:

-1.631+.1936*K11+.2078e-1*K12=-.851

.1787e-1*K12-.1846*K11+.8420=.280

 

Решим уравнения относительно К11 и К12:

[K11,K12]=solve('-1.631+.1936*K11+.2078e-1*K12=-.851','.1787e-1*K12-.1846*K11+.8420=.28')

В результате получим:

K11=3.511

K12=4.823

 

Построение  переходных процессов по ω_У и β: 

Ko=[3.511 4.822]

Ao=Aup+Bup*Ko

S1=ss(Ao,Bup,Cup,Dup)

SD1=c2d(S1,T0)

 

При ступенчатом отклонении руля направления:

step(S1,4)

3.Уравнения  состояния дискретной системы

с учётом обратных связей в канале руля направления.

Составим модель замкнутой системы с учётом динамики приводов:

 

M=[0 0 0 -1/0.11 0; 0 0 0 0 -1/0.11]

Apr=[A B; M]

Bpr=-1*M'

Cpr=eye(5)

Dpr=zeros(5,2)

Spr=ss(Apr,Bpr,Cpr,Dpr)

SDpr=c2d(Spr,0.3)

F=SDpr.a

D1=SDpr.b

H=SDpr.c

G=SDpr.d

Kpr=[3.511 4.823 0 0 0;0 0 0 0 0]

AprZ=F+D1*Kpr

SDprzam=ss(AprZ,D1,H,G,T0)

F1=SDprzam.a

D2=SDprzam.b

H1=SDprzam.c

G1=SDprzam.d

 

F1 =

0.4121   -1.5038   -0.0114   -0.0848    0.0139

    0.1414    0.7829    0.0089   -0.0131   -0.0025

   -1.0484   -4.2267    0.5865   -0.1059   -0.4346

    2.9411    4.0401         0    0.1623         0

         0         0         0         0    0.1623

D2 =

   -0.1068    0.0148

   -0.0096   -0.0019

   -0.1633   -0.5832

    0.8377         0

         0    0.8377

H1 =

     1     0     0     0     0

     0     1     0     0     0

     0     0     1     0     0

     0     0     0     1     0

     0     0     0     0     1

G1 =

     0     0

     0     0

     0     0

     0     0

     0     0

 

 

4.Синтез цифрового алгоритма управления в

канале элеронов частотным методом.

 

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

Найдем z-передаточную функцию ω_Х по uэ:

 

Wprwxue=tf(SDprzam(3,2))

-0.5832 z^4 + 0.5148 z^3 - 0.2439 z^2 - 0.09572 z + 0.08793

------------------------------------------------------------

z^5 - 2.106 z^4 + 2.169 z^3 - 1.251 z^2 + 0.3463 z - 0.03119

 

Pch=-0.5832*z^4 + 0.5148*z^3 - 0.2439*z^2 - 0.09572*z + 0.08793

Pzn=z^5 - 2.106*z^4 + 2.169*z^3 - 1.251*z^2 + 0.3463*z - 0.03119

syms w

z=(1+0.1*w)/(1-0.1*w)

Pch1=vpa(subs(Pch),4)

Pzn1=vpa(subs(Pzn),4)

 

Pch1=-.5832*(1.+.1000*w)^4/(1.-.1000*w)^4+.5148*(1.+.1000*w)^3/(1.-.1000*w)^3-.2439*(1.+.1000*w)^2/(1.-.1000*w)^2-.9572e-1*(1.+.1000*w)/(1.-.1000*w)+.8793e-1

Pzn1=(1.+.1000*w)^5/(1.-.1000*w)^5-2.106*(1.+.1000*w)^4/(1.-.1000*w)^4+2.169*(1.+.1000*w)^3/(1.-.1000*w)^3-1.251*(1.+.1000*w)^2/(1.-.1000*w)^2+.3463*(1.+.1000*w)/(1.-.1000*w)-.3119e-1

 

Для получения передаточной функции разложим многочлены числителя  и знаменателя по степеням:

a1=expand(-0.5832*(1+1/10*w)^4*(1-1/10*w))

a2=expand(0.5148*(1+1/10*w)^3*(1-1/10*w)^2)

a3=expand(-0.2439*(1+1/10*w)^2*(1-1/10*w)^3)

a4=expand(-0.09572*(1+1/10*w)*(1-1/10*w)^4)

a5=expand(0.08793*(1-1/10*w)^5)

b1=expand((1+1/10*w)^5)

b2=expand(-2.106*(1+1/10*w)^4*(1-1/10*w))

b3=expand(2.162*(1+1/10*w)^3*(1-1/10*w)^2)

b4=expand(-1.251*(1+1/10*w)^2*(1-1/10*w)^3)

b5=expand(0.3463*(1+1/10*w)*(1-1/10*w)^4)

b6=expand(-0.03119*(1-1/10*w)^5)

 

Pch2=a1+a2+a3+a4+a5

Pzn2=b1+b2+b3+b4+b5+b6

Pch21=vpa(collect(Pch2,w),3)

Pzn21=vpa(collect(Pzn2,w),3)

 

Pch21=-.320-.114*w-.102e-1*w^2-.142e-2*w^3+.275e-3*w^4+.116e-4*w^5

Pzn21=.120+.121*w+.435e-1*w^2+.839e-2*w^3+.110e-2*w^4+.690e-4*w^5

WW=tf([.116e-4 .275e-3 -.142e-2 -.102e-1 -.114 -.320],[.690e-4 .110e-2 .839e-2 .435e-1 .121 .120])

margin(WW)

 

 

Псевдочастотная характеристика:

 

 

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

 

Исходя из качества переходного  процесса w_Х по u_Э выберем коэффициент К23=0.25

Kpr1=[3.511 4.823 0 0 0;0 0 0.25 0 0];

AprZ1=F+D1*Kpr1;

SDprzam1=ss(AprZ1,D1,H,G,T0);

step(SDprzam1(3,2))

5. Выбор коэффициента электрической связи от

ручки лётчика.

 

Условие обеспечения  требуемой эффективности управления:

20мм отклонения ручки  лётчика должно соответствовать  установившееся значения ω_X

20 град/с

 

Располагаемый статический  коэффициент усиления канала руля элеронов равен

К=-1.54

Тогда Кх=1/К 

Кх=-0.649

 

Система, включающая в  себя выбранный коэффициент Кх:

Wwxue1=tf(SDprzam1(3,2))

Kx=-.649

D3=SDprzam1.b

D4=D3*[1 0;0 Kx]

SYS1=ss(AprZ1,D4,H,G,T0);

step(-20*SYS1(3,2),5)

 

Передаточная функция  ω_Х по uэ:

-0.5832 z^4 + 0.5148 z^3 - 0.2439 z^2 - 0.09572 z + 0.08793

-----------------------------------------------------------------------------

z^5 - 1.96 z^4 + 2.04 z^3 - 1.19 z^2 + 0.3703 z – 0.05317

 

 

 

 

 

 

Переходный процесс ω_X при отклонении ручки лётчика на 20 мм:

 

SYS=d2c(SDprzam1)

step(-20*SYS(3,2)*Kx,5)

 

График, представленный на рисунке, показывает реакцию на скачкообразное отклонение ручки летчика по ω_x. Время переходного процесса (вход в 10% трубку) составляет:

t_ПП=0.494с,

что удовлетворяет требованиям, налагаемым на систему автоматического управления.

 

6.Уравнение дискретного алгоритма оптимальной фильтрации для вычисления неизменяемой координаты объекта. Программа моделирования и исследование динамики оптимального фильтра.

Исходные данные:

Измеряемые координаты:

σ₁ σ ₁=0,5

 σ ₂ σ ₂=1,5

Случайное возмущение в  уравнении ω_Y с

Управляющий сигнал:

 

for i=1:180, UP(i,1)=10*sign(sin(pi/2*i*0.1)), end

i=1:180;

plot(i*0.1,UP(i,1))

 

Уравнения объекта управления:

Уравнения алгоритма  оптимальной фильтрации:

Рекуррентный алгоритм оценивания имеет вид:

 

где F – собственная матрица объекта управления

 

      - оценка вектора состояния;

R – единичная матрица  3х3;

H – матрица измерений,  которая имеет вид:

     

Структурная схема алгоритма  оценивания:

Воспользовавшись приведённым  выше алгоритмом оценивания, получим  переходные процессы координат состояния  и их оценок:

 

Набор измерений

X(:,:,1)=zeros(3,1);                     

for k=2:180; ,X(:,:,k)=F*X(:,:,k-1)+Dk*[0;UP(k-1,1)]; ,Y(:,:,k)=H*X(:,:,k)+randn*Sigma;,end

Накапливание коэффициентов

R=[1 0 0;0 1 0;0 0 1];

P(:,:,1)=E*p0;

for k=1:180;,M(:,:,k+1)=F*P(:,:,k)*F'; ,P(:,:,k+1)=M(:,:,k+1)-M(:,:,k+1)*H'*(H*M(:,:,k+1)*H'+R)^(-1)*H*M(:,:,k+1); ,end

Оценки 

Xk(:,:,1)=zeros(3,1);

for k=2:180;,Xk(:,:,k)=F*Xk(:,:,k-1)+Dk*[0;UP(k-1,1)]+M(:,:,k)*H'*(H*M(:,:,k)*H'+R)^(-1)*(Y(:,:,k)-H*(F*Xk(:,:,k-1)+Dk*[0;UP(k-1,1)])); ,end  
 Координата ω_Y:

j=1:1:180;

WY(j)=X(1,1,j);

WYk(j)=Xk(1,1,j);

WYr(j)=Y(1,1,j);

plot(j*0.1,WYk(j),j*0.1,WY(j),j*0.1,WYr(j)),grid

 

 

 

 

 

Измерение ω_у

 

 

 

 

 

 

 

 

 

Доступные измерения  ω_у

 

Оценка измерений ω_у фильтром Калмана

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Координата b:

Betar(j)=Y(2,1,j);

Beta(j)=X(2,1,j);

Betak(j)=Xk(2,1,j);

plot(j*0.1,Betak(j),j*0.1,Beta(j),j*0.1,Betar(j)),grid

 

 

 

Измерения β.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Доступные измерения  β.

 

 

Оценка  измерения  β фильтром Калмана.

 

 Координата ω_X:

WXr(j)=Y(3,1,j);

WX(j)=X(3,1,j);

WXk(j)=Xk(3,1,j);

plot(j*0.1,WXk(j),j*0.1,WX(j),j*0.1,WXr(j)),grid

 

 

 

Измерения ω_X.

 

 

 

 

 

 

 

 

 

 

 

Доступные измерения ω_X.

 

Оценка измерений ω_X фильтром Калмана.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

График изменения ошибки измерения координаты ω_у ( ):

WYi(j)=WY(j)-WYk(j);

plot(j*0.1,WYi(j)),grid

 

 

 

График изменения ошибки измерения координаты β ( ):

 

Betai(j)=Beta(j)-Betak(j);

plot(j*0.1,Betai(j)),grid

 

 

 

 

 

 

 

 

 

График изменения ошибки измерения координаты ω_х ( ):

 

Wxi(j)=WX(j)-WXk(j);

plot(j*0.1,Wxi(j)),grid

 

 

В результате моделирования  фильтра Калмана можно сделать  следующие выводы:

1) Оценка получилась  несмещенной, с минимальной дисперсией.

2) В начале измерений  фильтр Калмана дает большую ошибку, однако за несколько тактов ошибка начинает стремиться к нулевому значению.

 

 

Выводы.

 

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

1. Собственная частота колебаний в системе w* = 4.5 рад / с
2. Время переходного процесса по w_х   t_ПП=0.494 с
3. Затухание короткопериодических колебаний более чем в 10 раз за период.
4. Более чем двукратные запасы устойчивости на увеличение коэффициентов К₁₁, К₁₂, К₂₃.

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

 

1. Шамриков Б.М. «Теория цифровых систем управления летательным аппаратом», Москва, , 1987.

2. Ануфриев И.Е. Самоучитель по MatLab 5.3/6.x –СПБ.: БХВ-Петербург, 2003. – 736 c.