Курсовая работа: Математическая модель в пространстве состояний линейного стационарного объекта управления
ks4*h4_Tt + ks5*h5_Tt))), 50)
% u = vpa(u,6)
u_0 = subs(u,t,0)
u_T = subs(u,t,T)
ezplot(u, [0 T], 1)
hl=legend('u(t)');
set(hl, 'FontName', 'Courier');
title ('u(t)');
xlabel('t')
grid on
%
------------------------------------------------------------------------%
% ------------------------------------------------------------------------%
% Нахождения X
% Вычисление матричной экспоненты
MatrEx = simplify (vpa(ilaplace(inv(s*eye(5) - A)), 50));
syms t tay
X_svob = MatrEx * X_0;
X_vinyg = int ((subs(MatrEx, t, t - tay))*B*(subs (u, t, tay)), tay,
0,t);
X_real = X_svob + X_vinyg;
save Sostoyaniya X_real u
X_real = vpa (expand (simplify(X_real)), 50)
X_real_0 = double(subs (X_real, t, 0))
X_real_T = double(subs (X_real, t, T))
% Погрешность X
delta_X_T = double(vpa(X_T - X_real_T, 50))
delta_X_0 = double(vpa(X_0 - X_real_0, 50))
% Нахождение Y
for i = 1 : poryadok - 1
Y_real(i) = B_(i,:) * X_real;
end
Y_real = vpa (expand(simplify(Y_real')), 50)
Y_real_0 = double(subs (Y_real, t, 0))
Y_real_T = double(subs (Y_real, t, T))
% Погрешность Y
delta_Y_T = double(vpa(Y_T - Y_real_T, 50))
delta_Y_0 = double(vpa(Y_0 - Y_real_0, 50))
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Вычисление max значений для задачи АКОР
h = 0.01;
tic
tt = 0 : h : T;
for i = 1 : poryadok
X_max(i) = max(abs(subs(X_real(i),t,tt)));
end
U_max = max(abs(subs(u,t,tt)));
toc
save Sostoyaniya X_max U_max
%
------------------------------------------------------------------------%
% ------------------------------------------------------------------------%
%
Построение результатов X(t)
ezplot
(X_real(1), [0 T],2)
title ('x_1(t)');
grid on
ezplot (X_real(2), [0 T],3)
title ('x_2(t)');
grid on
ezplot (X_real(3), [0 T],4)
title ('x_3(t)');
grid on
ezplot (X_real(4), [0 T],5)
title ('x_4(t)');
grid on
ezplot (X_real(5), [0 T],6)
title ('x_5(t)');
grid on
%
Построение результатов Y(t)
ezplot (Y_real(1), [0 T],7)
title ('y_1(t)');
grid on
ezplot (Y_real(2), [0 T],8)
title ('y_2(t)');
grid on
ezplot (Y_real(3), [0 T],9)
title ('y_3(t)');
grid on
ezplot (Y_real(4), [0 T],10)
title ('y_4(t)');
grid on
%
------------------------------------------------------------------------%
clc
close all
clear all
format long
%
------------------------------------------------------------------------%
b_0 =
5;
b_1 =
9;
%
Укороченная система данного объекта
a_5 = 0.1153;
a_4 = 1.78;
a_3 = 3.92;
a_2 = 14.42;
a_1 = 8.583;
a_0 = 0;
% ------------------------------------------------------------------------%
% Приведение системы
b0 = b_0/a_5;
b1 = b_1/a_5;
a5 = a_5/a_5;
a4 = a_4/a_5;
a3 = a_3/a_5;
a2 = a_2/a_5;
a1 = a_1/a_5;
a0 = a_0/a_5;
%
------------------------------------------------------------------------%
%
Порядок системы
poryadok
= 5;
%
Начальные и конечные условия относительно вектора Y
Y_0 =
[3 2 1 5]';
Y_T =
[0 -1 0 3]';
%
Конечное время перехода
T =
3;
%
Матрица перехода от Н.У. Y к Н.У. X
B_ =
[b0 b1 0 0 0;
0 b0
b1 0 0;
0 0 b0
b1 0;
0 0 0
b0 b1];
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Начальные условия для упорядоченной системы
X_0 = B_' * inv(B_ * B_') * Y_0
X_T = B_' * inv(B_ * B_') * Y_T
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Представление системы в пространстве состояний
A =
[0 1 0 0 0;
0 0 1
0 0;
0 0 0
1 0
0 0 0
0 1;
-a0
-a1 -a2 -a3 -a4];
B =
[0; 0; 0; 0; 1];
C =
[b0 b1 0 0 0];
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Вычисление матричной экспоненты
syms s t
MatrEx = simplify (vpa(ilaplace(inv(s*eye(5) - A)), 50));
MatrEx_T = vpa(subs(MatrEx, t, T),50);
MatrEx_Tt = vpa(subs(MatrEx, t, T-t),50);
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Вычисление матрицы управляемости
M_c =
[B A*B A^2*B A^3*B A^4*B]
rank_M_c
= rank(M_c); %ранк = 5 - система управляема
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Вычисление грамиана управляемости
W_Tt
= double(vpa(simplify(int(MatrEx_Tt*B*B'*MatrEx_Tt',t,0,T)),50))
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Формирование управления
u =
vpa(expand(simplify(B'*MatrEx_Tt'*inv(W_Tt)*(X_T-MatrEx_T*X_0))),50)
u_0 = subs(u,t,0)
u_T = subs(u,t,T)
u = vpa(u,6)
%
------------------------------------------------------------------------%
ezplot(u, [0 T], 1)
title ('u(t)');
xlabel('t')
grid on
tt = 0 : 0.01 : T;
u2 = -20.605579750692850622177761310569*exp(-40.749492463732569440253455897187+13.583164154577523146751151965729*t)+19.011167813350479567880663060491*exp(-2.0544534472800777280645828326668+.68481781576002590935486094422228*t)+1.3356706538317879679656856470126*exp(-1.7550088311372150108106250409710+.58500294371240500360354168032368*t)*cos(-8.3032397968812277095785721047505+2.7677465989604092365261907015835*t)+7.2830359327562658520685140088852*exp(-1.7550088311372150108106250409710+.58500294371240500360354168032368*t)*sin(-8.3032397968812277095785721047505+2.7677465989604092365261907015835*t)-8.6096491449877801097840179781687;
u1 = subs(u2, t, tt);
u2 = subs(u, t, tt);
figure(2)
plot(tt,u1,'r',tt,u2,'b','LineWidth',2)
hl=legend('u(t)
при решении оптимальной L-проблемы моментов','u(t) с использованием грамиана
управляемости');
set(hl, 'FontName', 'Courier');
xlabel('t, cek'); ylabel('u(t)');
title('u(t)')
grid on
clc
clear all
close all
poryadok = 5;
%
------------------------------------------------------------------------%
b_0 =
5;
b_1 =
9;
%
Укороченная система данного объекта
a_5 = 0.1153;
a_4 = 1.78;
a_3 = 3.92;
a_2 = 14.42;
a_1 = 8.583;
a_0 = 0;
% ------------------------------------------------------------------------%
% Приведение системы
b0 = b_0/a_5;
b1 = b_1/a_5;
a5 = a_5/a_5;
a4 = a_4/a_5;
a3 = a_3/a_5;
a2 = a_2/a_5;
a1 = a_1/a_5;
a0 = a_0/a_5;
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Представление системы в пространстве состояний
A =
[0 1 0 0 0;
0 0 1 0 0;
0 0 0 1 0;
0 0 0 0 1;
-a0 -a1 -a2 -a3 -a4]
B = [0; 0; 0; 0; 1]
C = [b0 b1 0 0 0]
% Начальные условия
X_0 =
[10; 0; 6; 4; 8]
%T =
1;
Time
= 1;
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Получение max значений из файла
load Sostoyaniya X_max U_max
%
------------------------------------------------------------------------%
%
Нахождение элементов матриц Q и R
r(1) = 0.1;
q(1) = 1/poryadok * r(1) * (U_max)^2 / (X_max(1))^2;
for i = 2 : poryadok
q(i) = q(1) * (X_max(1))^2 / (X_max(i))^2;
end
Q = diag(q)
R = diag(r)
% Для
изменения коэффициентов
%
Q(1,1) = Q(1,1);
%
Q(2,2) = Q(2,2);
%
Q(3,3) = Q(3,3);
%
Q(4,4) = Q(4,4);
%
Q(5,5) = Q(5,5);
Q(1,1)
= Q(1,1)*1e+12;
Q(2,2)
= Q(2,2)*1e+8;
Q(3,3)
= Q(3,3)*1e+7;
Q(4,4)
= Q(4,4)*1e+0;
Q(5,5)
= Q(5,5)*1e+2;
R(1,1)
= R(1,1);
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Решение уравнения Риккати методом диагонализации
P1 = Solve_Riccati_Method_Diag(A,B,Q,R)
%
------------------------------------------------------------------------%
P_nach = zeros(poryadok, poryadok);%+ones(poryadok, poryadok);
%
------------------------------------------------------------------------%
Страницы: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17 |