Решение систем дифференциальных уравнений методом преобразования Лапласа
Выводы
Построим графики движений и . Из графиков следует, что исследуемый объект является устойчивым. Из графика второго окна видно, при отсутствии воздействия значения компонент стремятся к нулю.
Однако экстремумы показывают, что значения функций f(t) - угловое отклонение от заданного курса и g(t) - угол крена во время воздействия отклоняются от нулевых почти на 20% , что может приводить к значительным флуктуациям курсовой стабильности.
Список литературы
1. В.А. Крамарь, В.А .Карапетьян Методические указания к выполнению курсовой работы по дисциплине «Специальные разделы математики» Севастополь. СевГУ. 2015.-30с.
2. Андык В.С. Теория автоматического управления. Учебное пособие практическим занятиям. Томск. Издательство ТПУ 2004.-108с.
3. Бесекерский В.А., Попов Е.П. Теория систем автоматического регулирования. М. Наука, 1972.-762с.
4. Дорф Р. Бишоп Р. Современные системы управления |Перевод с английского Б.И. Копылова| М. Лаборатория базовых знаний.
5. Линейно-квадратичное гауссовское управление
6. Преобразовамние Лапламса
Приложение
Выполнение задания в пакете MatLab
clear;
%Исходные данные модели ЛА
A = [0 2.6 0; 0 0 0.01; 0 0 0.05];
B = [0 0 0.05];
C = [1 0 0];
D = 0;
x0 = [40;0;-1];
syms t s;
%Оригиналы и изображения входных воздействий
Us_t = heaviside(t-30)-heaviside(t-10);
Ug_t = cos(0.25*t);
Us_s = laplace(Us_t); vpa(Us_s,6);
Ug_s = laplace(Ug_t); vpa(Ug_s,6);
%Определение резольвенты матрицы А
As = s*eye(3)-A; vpa(As,6);
Res = As^-1; vpa(Res,6);
As_pol = vpa(factor(det(As)),2);
%Определение изображений движений
Xc = Res*x0;
Xus_s = Res*B*Us_s; vpa(Xus_s,6);
Xug_s = Res*B*Ug_s; vpa(Xug_s,6);
%Определение оригиналов движений
Xc_t = ilaplace(Xc);
Xus_t = ilaplace(Xus_s); vpa(simplify(Xus_t),6);
Xug_t = ilaplace(Xug_s); vpa(simplify(Xug_t),6);
Xcs_t = Xc_t+Xus_t;
Xcg_t = Xc_t+Xug_t;
%Построение графиков движений нестабилизированного самолета
tk = 100;
time = 0:0.001:tk;
subplot(1,3,1);
ezplot(Xcg_t(1),time); grid on
title(x_1(t))
subplot(1,3,2);
ezplot(Xcg_t(2),time); grid on
title(x_2(t))
subplot(1,3,3);
ezplot(Xcg_t(3),time); grid on
title(x_3(t))
%Построение модели стабилизированного самолета
syms Ky Kf Kg
Koc = [Ky Kf Kg];
Aoc = A-B*Koc; vpa(Aoc,4);
pol = collect(det(s*eye(3)-Aoc),s); vpa(pol,8);
koef = poly([-0.02 -0.02 -0.05]); vpa(koef,8);
K = acker(A,B,[-0.02 -0.02 -0.05]);
Aoc = A-B*K;
eigs(Aoc);
Aoc_s=s*eye(3)-Aoc; vpa(Aoc_s,6);
Res_oc=Aoc_s^-1; vpa(Res_oc,6);
%Определение изображений движений
Xc_oc = Res_os*x0;
Xus_s_oc = Res_oc*B*Us_s; vpa(Xus_s_oc,6);
Xug_s_oc = Res_oc*B*Ug_s; vpa(Xug_s_oc,6);
%Определение оригиналов движений
Xc_t_oc = ilaplace(Xc_oc);
Xus_t_oc = ilaplace(Xus_s_oc); vpa(Xus_t_oc,6);
Xug_t_oc = ilaplace(Xug_s_oc); vpa(Xug_t_oc,6);
Xcs_t_oc = Xc_t_oc+Xus_t_oc; vpa(simplify(Xcs_t_oc),6);
Xcg_t_oc = Xc_t_oc+Xug_t_oc; vpa(simplify(Xcg_t_oc),6);
%Построение графиков движений стабилизированного ЛА
figure
subplot(1,3,1);
ezplot(Xcs_t_oc(1),time); grid on
title(x_1(t))
subplot(1,3,2);
ezplot(Xcs_t_oc(2),time); grid on
title(x_2(t))
subplot(1,3,3);
ezplot(Xcs_t_oc(3),time); grid on
A =
0 2.6000 0
0 0 0.0100
0 0 0.0500
B =
0
0
0.0500
C =
1 0 0
D = 0
x0 =
40
0
-1
Us_t = heaviside(t - 30) - heaviside(t - 10)
Ug_t = cos(t/4)
ans = exp(-30.0*s)/s - (1.0*exp(-10.0*s))/s
ans = s/(s^2 + 0.0625)
ans =
[ s, -2.6, 0]
[ 0, s, -0.01]
[ 0, 0, s - 0.05]
ans =
[ 1/s, 2.6/s^2, 0.52/(s^2*(20.0*s - 1.0))]
[ 0, 1/s, 0.2/(s*(20.0*s - 1.0))]
[ 0, 0, 20.0/(20.0*s - 1.0)]
As_pol = [ 0.05, s, s, 20.0*s - 1.0]
Xc =
40/s - 13/(25*s^2*(20*s - 1))
-1/(5*s*(20*s - 1))
-20/(20*s - 1)
ans =
-(0.026*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(s^2*(20.0*s - 1.0))
-(0.01*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(s*(20.0*s - 1.0))
-(1.0*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(20.0*s - 1.0)
ans =
0.026/(s*(20.0*s - 1.0)*(s^2 + 0.0625))
0.01/((20.0*s - 1.0)*(s^2 + 0.0625))
s/((20.0*s - 1.0)*(s^2 + 0.0625))
Xc_t =
(13*t)/25 - (52*exp(t/20))/5 + 252/5
1/5 - exp(t/20)/5
-exp(t/20)
ans =
0.026*heaviside(1.0*t - 10.0)*(10.0*t - 400.0*exp(0.05*t - 0.5) + 0.5*t^2 + 250.0) + 0.026*heaviside(1.0*t - 30.0)*(10.0*t + 400.0*exp(0.05*t - 1.5) - 0.5*t^2 - 250.0)
0.01*heaviside(1.0*t - 30.0)*(20.0*exp(0.05*t - 1.5) - 1.0*t + 10.0) + 0.01*heaviside(1.0*t - 10.0)*(t - 20.0*exp(0.05*t - 0.5) + 10.0)
heaviside(1.0*t - 30.0)*(exp(0.05*t - 1.5) - 1.0) - 1.0*heaviside(1.0*t - 10.0)*(exp(0.05*t - 0.5) - 1.0)
ans =
0.016*cos(0.25*t) - 0.08*sin(0.25*t) + 0.4*exp(0.05*t) - 0.416
0.00769231*exp(0.05*t) - 0.00153846*sin(0.25*t) - 0.00769231*cos(0.25*t)
0.192308*sin(0.25*t) - 0.0384615*cos(0.25*t) + 0.0384615*exp(0.05*t)
Xcs_t =
(13*t)/25 - (52*exp(t/20))/5 + (13*heaviside(t - 10)*(20*t - 400*exp(t/20 - 1/2) + (t - 10)^2/2 + 200))/500 - (13*heaviside(t - 30)*(20*t - 400*exp(t/20 - 3/2) + (t - 30)^2/2 - 200))/500 + 252/5
(heaviside(t - 10)*(t - 20*exp(t/20 - 1/2) + 10))/100 - exp(t/20)/5 + (heaviside(t - 30)*(20*exp(t/20 - 3/2) - t + 10))/100 + 1/5
heaviside(t - 30)*(exp(t/20 - 3/2) - 1) - heaviside(t - 10)*(exp(t/20 - 1/2) - 1) - exp(t/20)
Xcg_t =
(13*t)/25 + (2*cos(t/4))/125 - 10*exp(t/20) - (2*sin(t/4))/25 + 6248/125
1/5 - (5*exp(t/20))/26 - sin(t/4)/650 - cos(t/4)/130
(5*sin(t/4))/26 - (25*exp(t/20))/26 - cos(t/4)/26
Koc = [ Ky, Kf, Kg]
ans =
[ 0, 2.6, 0]
[ 0, 0, 0.01]
[ -0.05*Ky, -0.05*Kf, 0.05 - 0.05*Kg]
ans = s^3 + (0.05*Kg - 0.05)*s^2 + 0.0005*Kf*s + 0.0013*Ky
ans = [ 1.0, 0.09, 0.0024, 0.00002]
K = 0.0154 4.8000 2.8000
Aoc =
0 2.6000 0
0 0 0.0100
-0.0008 -0.2400 -0.0900
ans =
-0.0200 + 0.0000i
-0.0200 - 0.0000i
-0.0500 + 0.0000i
ans =
[ s, -2.6, 0]
[ 0, s, -0.01]
[ 0.000769231, 0.24, s + 0.09]
ans =
[ (20.0*(2500.0*s^2 + 225.0*s + 6.0))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), (1300.0*(100.0*s + 9.0))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), 1300.0/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)]
[ -0.384615/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), (500.0*s*(100.0*s + 9.0))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), (500.0*s)/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)]
[ -(38.4615*s)/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), -(100.0*(120.0*s + 1.0))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), (50000.0*s^2)/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)]
Xc_oc =
(800*(2500*s^2 + 225*s + 6))/(50000*s^3 + 4500*s^2 + 120*s + 1) - 1300/(50000*s^3 + 4500*s^2 + 120*s + 1)
- (500*s)/(50000*s^3 + 4500*s^2 + 120*s + 1) - 200/(13*(50000*s^3 + 4500*s^2 + 120*s + 1))
- (50000*s^2)/(50000*s^3 + 4500*s^2 + 120*s + 1) - (20000*s)/(13*(50000*s^3 + 4500*s^2 + 120*s + 1))
ans =
-(65.0*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)
-(25.0*s*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)
-(2500.0*s^2*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)
ans =
(65.0*s)/((s^2 + 0.0625)*(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0))
(25.0*s^2)/((s^2 + 0.0625)*(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0))
(2500.0*s^3)/((s^2 + 0.0625)*(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0))
Xc_t_oc =
(460*exp(-t/50))/9 - (100*exp(-t/20))/9 + (7*t*exp(-t/50))/15
(25*exp(-t/20))/117 - (25*exp(-t/50))/117 - (7*t*exp(-t/50))/1950
(8*exp(-t/50))/117 - (125*exp(-t/20))/117 + (7*t*exp(-t/50))/975
ans =
65.0*heaviside(1.0*t - 10.0)*(0.444444*exp(0.5 - 0.05*t) + 0.555556*exp(0.2 - 0.02*t) + 0.0333333*exp(0.2 - 0.02*t)*(t - 10.0) - 1.0) - 65.0*heaviside(1.0*t - 30.0)*(0.444444*exp(1.5 - 0.05*t) + 0.555556*exp(0.6 - 0.02*t) + 0.0333333*exp(0.6 - 0.02*t)*(t - 30.0) - 1.0)
25.0*heaviside(1.0*t - 30.0)*(0.0222222*exp(1.5 - 0.05*t) - 0.0222222*exp(0.6 - 0.02*t) + 0.000666667*exp(0.6 - 0.02*t)*(t - 30.0)) - 25.0*heaviside(1.0*t - 10.0)*(0.0222222*exp(0.5 - 0.05*t) - 0.0222222*exp(0.2 - 0.02*t) + 0.000666667*exp(0.2 - 0.02*t)*(t - 10.0))
2500.0*heaviside(1.0*t - 10.0)*(0.00111111*exp(0.5 - 0.05*t) - 0.00111111*exp(0.2 - 0.02*t) + 0.0000133333*exp(0.2 - 0.02*t)*(t - 10.0)) - 2500.0*heaviside(1.0*t - 30.0)*(0.00111111*exp(1.5 - 0.05*t) - 0.00111111*exp(0.6 - 0.02*t) + 0.0000133333*exp(0.6 - 0.02*t)*(t - 30.0))
ans =
1.13944*exp(-0.02*t) - 0.0759527*sin(0.25*t) - 1.11111*exp(-0.05*t) - 0.0283338*cos(0.25*t) - 0.0137785*t*exp(-0.02*t)
0.0027244*sin(0.25*t) - 0.00730314*cos(0.25*t) + 0.0213675*exp(-0.05*t) - 0.0140644*exp(-0.02*t) + 0.000105988*t*exp(-0.02*t)
0.06811*cos(0.25*t) + 0.182579*sin(0.25*t) - 0.106838*exp(-0.05*t) + 0.0387276*exp(-0.02*t) - 0.000211977*t*exp(-0.02*t)
ans =
51.1111*exp(-0.02*t) - 11.1111*exp(-0.05*t) - 65.0*heaviside(1.0*t - 30.0)*(0.444444*exp(1.5 - 0.05*t) - 0.444444*exp(0.6 - 0.02*t) + 0.0333333*t*exp(0.6 - 0.02*t) - 1.0) + 65.0*heaviside(1.0*t - 10.0)*(0.444444*exp(0.5 - 0.05*t) + 0.222222*exp(0.2 - 0.02*t) + 0.0333333*t*exp(0.2 - 0.02*t) - 1.0) + 0.466667*t*exp(-0.02*t)
0.213675*exp(-0.05*t) - 0.213675*exp(-0.02*t) - 25.0*heaviside(1.0*t - 10.0)*(0.0222222*exp(0.5 - 0.05*t) - 0.0222222*exp(0.2 - 0.02*t) + 0.000666667*exp(0.2 - 0.02*t)*(t - 10.0)) + 25.0*heaviside(1.0*t - 30.0)*(0.0222222*exp(1.5 - 0.05*t) - 0.0222222*exp(0.6 - 0.02*t) + 0.000666667*exp(0.6 - 0.02*t)*(t - 30.0)) - 0.00358974*t*exp(-0.02*t)
0.0683761*exp(-0.02*t) - 1.06838*exp(-0.05*t) + 2500.0*heaviside(1.0*t - 10.0)*(0.00111111*exp(0.5 - 0.05*t) - 0.00111111*exp(0.2 - 0.02*t) + 0.0000133333*exp(0.2 - 0.02*t)*(t - 10.0)) + 0.00717949*t*exp(-0.02*t) - 2500.0*heaviside(1.0*t - 30.0)*(0.00111111*exp(1.5 - 0.05*t) - 0.00111111*exp(0.6 - 0.02*t) + 0.0000133333*exp(0.6 - 0.02*t)*(t - 30.0))
ans =
52.2506*exp(-0.02*t) - 0.0759527*sin(0.25*t) - 12.2222*exp(-0.05*t) - 0.0283338*cos(0.25*t) + 0.452888*t*exp(-0.02*t)
0.0027244*sin(0.25*t) - 0.00730314*cos(0.25*t) + 0.235043*exp(-0.05*t) - 0.22774*exp(-0.02*t) - 0.00348376*t*exp(-0.02*t)
0.06811*cos(0.25*t) + 0.182579*sin(0.25*t) - 1.17521*exp(-0.05*t) + 0.107104*exp(-0.02*t) + 0.00696751*t*exp(-0.02*t)
>>