Курсовая работа: Математическая модель в пространстве состояний линейного стационарного объекта управления
P_nach = zeros(poryadok, poryadok);%+ones(poryadok, poryadok);
%
------------------------------------------------------------------------%
%
Решение уравнения Риккати методом обратного интегрирования
P2 = Solve_Riccati_Method_Revers_Integr(A,B,Q,R,Time,poryadok, P_nach)
%
------------------------------------------------------------------------%
%
Сравнение расхождения методов
Delta_P
= abs(P1-P2)
%
Построение графика коэффициентов регулятора
load
Solve_Riccati_Method_Revers_Integr Time_R P N_str
PP = P;
for i = 1 : N_str
P = reshape(PP(i, :), poryadok, poryadok);
K(i, :) = -inv(R)*B'*P;
end
figure(2)
plot(Time_R,K(:,1),'-',Time_R,K(:,2),'-',Time_R,K(:,3),'-',Time_R,K(:,4),'-',Time_R,K(:,5),'-',
'LineWidth', 2);
xlabel('t')
tit1
= title('Коэффициенты обратной связи в прямом времени');
set(tit1,'FontName','Courier');
hl=legend('k_1_о_с','k_2_о_с','k_3_о_с','k_4_о_с','k_5_о_с',0);
set(hl,'FontName','Courier');
grid on;
%
------------------------------------------------------------------------%
%
Решение уравнения Риккати с помощью встроенной функции
% P = vpa(care(A,B,Q,R), 10)
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Нахождение коэффициентов регулятора
disp('Коэффициенты
регулятора:')
K1 = -inv(R) * B' * P1
K2 = -inv(R) * B' * P2
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
A1_ =
A + B * K1;
A2_ =
A + B * K2;
%
Вычисление матричной экспоненты
syms s t
MatrEx1 = simplify (vpa(ilaplace(inv(s*eye(5) - A1_)), 50));
MatrEx2 = simplify (vpa(ilaplace(inv(s*eye(5) - A2_)), 50));
%
Нахождение координат состояния
X1 =
vpa(simplify(MatrEx1 * X_0), 50);
X2 = vpa(simplify(MatrEx2 * X_0), 50);
%
Нахождение управления
u1 =
vpa(simplify(K1 * X1),50)
u2 = vpa(simplify(K2 * X2),50)
%
------------------------------------------------------------------------%
%
Построение u(t) и X(t)
T_sravneniya = 0.2;
figure(3);
tt = 0 : 0.01 : T_sravneniya;
uu1 = subs(u1,t,tt);
uu2 = subs(u2,t,tt);
plot(tt, uu1, tt, uu2, 'LineWidth', 2)
title ('u(t)');
xlabel('t')
hl=legend('u(t) - управление',0);
set(hl,'FontName','Courier');
grid on
ezplot(X1(1), [0 Time], 4)
hold on
title ('x_1(t)');
xlabel('t')
grid on
ezplot(X1(2), [0 Time], 5)
title ('x_2(t)');
xlabel('t')
grid on
ezplot(X1(3), [0 Time], 6)
title ('x_3(t)');
xlabel('t')
grid on
ezplot(X1(4), [0 Time], 7)
title ('x_4(t)');
xlabel('t')
grid on
ezplot(X1(5), [0 Time], 8)
title ('x_5(t)');
xlabel('t')
grid on
tt = 0 : 0.01 : T_sravneniya;
X21 = subs(X1(1), t, tt);
X22= subs(X1(2), t, tt);
X23= subs(X1(3), t, tt);
X24= subs(X1(4), t, tt);
X25= subs(X1(5), t, tt);
save Sravnenie_stabilizacii_1 X21 X22 X23 X24 X25 uu1
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];
Time
= 0.2;
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Получение max значений из файла
load Sostoyaniya X_max U_max
%
------------------------------------------------------------------------%
%
Нахождение элементов матриц Q и R
% r(1) = 100;
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)*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);
% P_prib = eye(poryadok, poryadok);
% P_prib(1,1) = 100;
% P_prib(2,2) = 10;
% % P_prib(3,3) = 1000;
% % P_prib(4,4) = 10;
% % P_prib(5,5) = 1;
% ------------------------------------------------------------------------%
P_nach = zeros(poryadok, poryadok);% + P_prib;
%
------------------------------------------------------------------------%
%
Решение уравнения Риккати методом обратного интегрирования
P = Solve_Riccati_Method_Revers_Integr(A,B,Q,R,Time,poryadok, P_nach)
%
------------------------------------------------------------------------%
% Нахождение переменных коэффициентов регулятора
load Solve_Riccati_Method_Revers_Integr Time_R P N_str
PP = P;
for i = 1 : N_str
P = reshape(PP(i, :), poryadok, poryadok);
K(i, :) = -inv(R)*B'*P;
end
%
------------------------------------------------------------------------%
%
Формирование вектора коэффициентов регулятора
% и
решения уравнения Риккати в прямом порядке
load Solve_Riccati_Method_Revers_Integr P
size(K)
i = 1;
len_K = length(K(:,1))
for j = len_K : -1 : 1
K_pr(i,:) = K(j,:);
i = i
+ 1;
end
%
------------------------------------------------------------------------%
%
Построение графика переменных коэффициентов регулятора в прямом времени
figure(2)
plot(Time_R,K(:,1),'-',Time_R,K(:,2),'-',Time_R,K(:,3),'-',...
Time_R,K(:,4),'-',Time_R,K(:,5),'-', 'LineWidth', 2);
grid on;
title('K(t)')
xlabel('t')
legend('k_1','k_2','k_3','k_4','k_5');
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
for k = 1 : len_K
A_(:,:,k) = A + B * K(k,:);
end
size(A_);
%
------------------------------------------------------------------------%
%
Нахождение фазовых координат
X(:,1) = X_0;
h = 0.01;
time_X(1) = 0;
for k = 1 : len_K
X(:, k+1) = X(:, k) + h * A_(:,:,k) * X(:, k);
time_X(k+1) = time_X(k) + h;
end
X(:, k+1) = [];
time_X(k+1) = [];
%
------------------------------------------------------------------------%
%
Нахождение управления
for k
= 1 : len_K
u(k) = K_pr(k,:) * X(:,k);
end
%
------------------------------------------------------------------------%
%
Построение u(t) и X(t)
figure(3);
plot(time_X, u, 'r-', 'LineWidth', 2)
title ('u(t)');
xlabel('t')
grid on
figure(4);
plot(time_X, X(1,:), 'LineWidth', 2)
hold on
title ('x_1(t)');
xlabel('t')
grid on
figure(5);
plot(time_X, X(2,:), 'LineWidth', 2)
title ('x_2(t)');
xlabel('t')
grid on
figure(6);
plot(time_X, X(3,:), 'LineWidth', 2)
title ('x_3(t)');
xlabel('t')
grid on
figure(7);
plot(time_X, X(4,:), 'LineWidth', 2)
title ('x_4(t)');
xlabel('t')
grid on
figure(8);
plot(time_X, X(5,:), 'LineWidth', 2)
title ('x_5(t)');
xlabel('t')
grid on
save Sravnenie_stabilizacii_2 time_X X u
close all
load Sravnenie_stabilizacii_1 X21 X22 X23 X24 X25 uu1
load Sravnenie_stabilizacii_2 time_X X u
figure(31);
plot(time_X, u, time_X, uu1, 'LineWidth', 2)
title ('u(t)');
xlabel('t')
hl=legend('u(t)
- управление с перемен. коеф.','u(t) - управление с пост. коеф.');
set(hl,'FontName','Courier');
grid on
figure(41);
plot(time_X, X(1,:), time_X, X21, 'LineWidth', 2)
hold on
title ('x_1(t)');
xlabel('t')
hl=legend('x_1(t)
- с перемен. коеф.','x_1(t) - с пост. коеф.');
set(hl,'FontName','Courier');
grid on
figure(51);
plot(time_X, X(2,:), time_X, X22,'LineWidth', 2)
title ('x_2(t)');
xlabel('t')
hl=legend('x_2(t)
- с перемен. коеф.','x_2(t) - с пост. коеф.');
set(hl,'FontName','Courier');
grid on
figure(61);
plot(time_X, X(3,:), time_X, X23,'LineWidth', 2)
title ('x_3(t)');
xlabel('t')
hl=legend('x_3(t)
- с перемен. коеф.','x_3(t) - с пост. коеф.');
set(hl,'FontName','Courier');
grid on
figure(71);
plot(time_X, X(4,:), time_X, X24,'LineWidth', 2)
title ('x_4(t)');
xlabel('t')
hl=legend('x_4(t)
- с перемен. коеф.','x_4(t) - с пост. коеф.');
set(hl,'FontName','Courier');
grid on
figure(81);
plot(time_X, X(5,:), time_X, X25,'LineWidth', 2)
title ('x_5(t)');
xlabel('t')
hl=legend('x_5(t)
- с перемен. коеф.','x_5(t) - с пост. коеф.');
set(hl,'FontName','Courier');
grid on
clc
clear all
close all
warning off
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];
Time
= 1;
h =
0.01;
%
------------------------------------------------------------------------%
tic
%
------------------------------------------------------------------------%
%
Получение max значений из файла
load Sostoyaniya X_max U_max
%
------------------------------------------------------------------------%
%
Нахождение элементов матриц Q и R
r(1) = 100;
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)*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);
% P_0 = ones(poryadok, poryadok);
% P_0(1,1) = P_0(1,1)*1e12;
% P_0(2,2) = P_0(2,2)*1e8;
% P_0(3,3) = P_0(3,3)*1e7;
% P_0(4,4) = P_0(4,4)*1e0;
% P_0(5,5) = P_0(5,5)*1e2;
%
------------------------------------------------------------------------%
P_nach = zeros(poryadok, poryadok);%+P_0;
%
------------------------------------------------------------------------%
%
Решение уравнения Риккати методом обратного интегрирования
P = Solve_Riccati_Method_Revers_Integr(A,B,Q,R,Time,poryadok, P_nach);
load Solve_Riccati_Method_Revers_Integr_for_slegenie Time_R P N_str
PP = P;
for k = 1 : N_str
P1 = reshape(PP(k, :), poryadok, poryadok);
for i = 1 : poryadok
for j = 1 : poryadok
P2(i,j,k) = P1(i,j);
end
end
end
size_P = size(P2);
%
------------------------------------------------------------------------%
tic
%
------------------------------------------------------------------------%
%
Получение дискретных значений задающего воздействия в обратном времени
% для
нахождения вспомогательной функции q(t)
Vozmyshyayushee_Vozdeistvie_Discrete_Revers(h, 0, Time);
%
------------------------------------------------------------------------%
load Vozmyshyayushee_Vozdeistvie_Discrete_Revers w_discrete_rev
% ------------------------------------------------------------------------%
size(w_discrete_rev);
% Начальное значение q(t)
q = zeros(poryadok,1);
%
Интегрирование q(t) в обратном времени
for k = 1 : N_str
q(:, k+1) = q(:, k) - h * ((P2(:,:,k)*B*inv(R)*B'-A') * q(:, k) -
P2(:,:,k)*w_discrete_rev(:,k));
end
q(:, k+1) = [];
size_q = size(q);
%
------------------------------------------------------------------------%
%
Нахождение переменных коэффициентов регулятора
for k
= 1 : N_str
K_o(k, :) = -inv(R) * B' * P2(:,:,k);
K_pr(k, :) = -inv(R) * B';
end
%
Формирование вектора коэффициентов регулятора, значений задающего
%
воздействия, значений вспомогательной функции в прямом порядке
size(K_o);
size(K_pr);
K_pr_p = K_pr;
i = 1;
len_K = length(K_o(:,1));
for j = len_K : -1 : 1
K_o_p(i,:) = K_o(j,:);
w_discrete(:,i) = w_discrete_rev(:,j);
q_pr(:, i) = q(:, j);
i = i + 1;
end
%
------------------------------------------------------------------------%
%
Построение графика переменных коэффициентов регулятора обратной связи
% в
прямом времени
toc
figure(3)
plot(Time_R,K_o(:,1),'-',Time_R,K_o(:,2),'-',Time_R,K_o(:,3),'-',...
Time_R,K_o(:,4),'-',Time_R,K_o(:,5),'-', 'LineWidth', 2);
xlabel('t')
tit1
= title('Коэффициенты обратной связи в прямом времени');
set(tit1,'FontName','Courier');
hl=legend('k_1_о_с','k_2_о_с','k_3_о_с','k_4_о_с','k_5_о_с',0);
set(hl,'FontName','Courier');
grid on;
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Построение графика переменных коэффициентов регулятора прямой связи
% в прямом
времени
figure(4)
plot(Time_R,K_pr(:,1),'-',Time_R,K_pr(:,2),'-',Time_R,K_pr(:,3),'-',...
Time_R,K_pr(:,4),'-',Time_R,K_pr(:,5),'-', 'LineWidth', 2);
xlabel('t')
tit1
= title('Коэффициенты прямой связи в прямом времени');
set(tit1,'FontName','Courier');
hl=legend('k_1_п_с','k_2_п_с','k_3_п_с','k_4_п_с','k_5_п_с',0);
set(hl,'FontName','Courier');
grid on;
%
------------------------------------------------------------------------%
tic
% ------------------------------------------------------------------------%
for k = 1 : len_K
A_(:,:,k) = A + B * K_o_p(k,:);
end
size_A_ = size(A_);
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Нахождение фазовых координат
X(:,1)
= X_0;
time_X(1) = 0;
for k = 1 : len_K
X(:, k+1) = X(:, k) + h * (A_(:,:,k) * X(:, k) + B * K_pr_p(k,:) *
q_pr(:,k) + w_discrete(:,k));
time_X(k+1) = time_X(k) + h;
end
X(:, k+1) = [];
time_X(k+1) = [];
size_X
= size(X);
%
------------------------------------------------------------------------%
%
------------------------------------------------------------------------%
%
Нахождение управления
for k = 1 : len_K
u(k) = K_o_p(k,:) * X(:,k) + K_pr_p(k,:) * q_pr(:,k);
end
size_u = size(u);
%
------------------------------------------------------------------------%
toc
%
Построение u(t) и X(t)
figure(5);
plot(time_X, u, 'r-', 'LineWidth', 2)
title ('u(t)');
xlabel('t')
hl=legend('u(t) - управление',0);
set(hl,'FontName','Courier');
grid on
figure(6);
plot(time_X, X(1,:),'r-', time_X, w_discrete(1,:), 'LineWidth', 2)
hold on
title ('x_1(t)');
xlabel('t');
hl=legend('X(t) - реальный сигнал','w(t) - возмущающее воздействие',0);
set(hl,'FontName','Courier');
grid on
figure(7);
plot(time_X, X(2,:),'r-', time_X, w_discrete(2,:), 'LineWidth', 2)
title ('x_2(t)');
xlabel('t');
hl=legend('X(t)
- реальный сигнал','w(t) - возмущающее воздействие',0);
set(hl,'FontName','Courier');
grid on
figure(8);
plot(time_X, X(3,:),'r-', time_X, w_discrete(3,:), 'LineWidth', 2)
title ('x_3(t)');
xlabel('t');
hl=legend('X(t)
- реальный сигнал','w(t) - возмущающее воздействие',0);
set(hl,'FontName','Courier');
grid on
figure(9);
plot(time_X, X(4,:),'r-', time_X, w_discrete(4,:), 'LineWidth', 2)
title ('x_4(t)');
xlabel('t');
hl=legend('X(t)
- реальный сигнал','w(t) - возмущающее воздействие',0);
set(hl,'FontName','Courier');
grid on
figure(10);
plot(time_X, X(5,:),'r-', time_X, w_discrete(5,:), 'LineWidth', 2)
title ('x_5(t)');
xlabel('t');
hl=legend('X(t)
- реальный сигнал','w(t) - возмущающее воздействие',0);
set(hl,'FontName','Courier');
grid on
figure(11);
plot(time_X, q(1,:), time_X, q(2,:), time_X, q(3,:), time_X, q(4,:),
time_X, q(5,:), 'LineWidth', 2)
title ('q(t)- vector-function');
xlabel('t');
hl=legend('q_1(t)', 'q_2(t)', 'q_3(t)', 'q_4(t)', 'q_5(t)');
set(hl,'FontName','Courier');
grid on
clc
clear all
close all
poryadok = 5;
%
------------------------------------------------------------------------%
Страницы: 1, 2, 3, 4, 5, 6, 7, 8
|