Математическая модель в пространстве состояний линейного стационарного объекта управления

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

AKOR_stabilizaciya_na_polybeskon_interval.m

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);

% ------------------------------------------------------------------------%

% Решение уравнения Риккати методом обратного интегрирования

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

AKOR_stabilizaciya_na_konech_interval.m

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

Sravnenie_stabilizacii.m

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

AKOR_stabilizaciya_pri_vozmusheniyah.m

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

AKOR_slegenie_na_konech_interval_I_podxod.m

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 = 1;

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

% Получение 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+10;

% Q(2,2) = Q(2,2)*1e+8;

% Q(3,3) = Q(3,3)*1e+6;

% Q(4,4) = Q(4,4)*1e+2;

% Q(5,5) = Q(5,5)*1e+2;

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);

% ------------------------------------------------------------------------%

% Задающее воздействие

A_o = [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];

X_o_0 = [12; 10; 14; 8; 16];

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

% Расширенный вектор состояния и расширенные матрицы A,B,Q

%X_rassh = [X_0; X_o];

NULL_M1 = zeros(size(A));

A_rassh = [A NULL_M1;

NULL_M1 A_o];

NULL_M2 = zeros(length(A(:,1)), 1);

B_rassh = [B; NULL_M2];

Q_rassh = [Q -Q;

-Q Q];

X_rassh_0 = [X_0; X_o_0]

% ------------------------------------------------------------------------%

P_nach = zeros(2*poryadok, 2*poryadok);%+ones(poryadok, poryadok);

% ------------------------------------------------------------------------%

% Решение уравнения Риккати методом обратного интегрирования

P = Solve_Riccati_Method_Revers_Integr(A_rassh,B_rassh,Q_rassh,R,Time,2*poryadok, P_nach)

% ------------------------------------------------------------------------%

% Нахождение переменных коэффициентов регулятора

load Solve_Riccati_Method_Revers_Integr_for_slegenie Time_R P N_str

% ------------------------------------------------------------------------%

% % Формирование матриц P11 и P12

PP = P;

for k = 1 : N_str

P = reshape(PP(k, :), 2*poryadok, 2*poryadok);

for i = 1 : poryadok

for j = 1 : poryadok

P11(i,j,k) = P(i,j);

end

end

for i = 1 : poryadok

for j = (poryadok+1) : (2*poryadok)

P12(i,j-poryadok,k) = P(i,j);

end

end

end

P11(:,:,k)

P12(:,:,k)

% ------------------------------------------------------------------------%

for k = 1 : N_str

K_o(k, :) = -inv(R) * B' * P11(:,:,k);

K_pr(k, :) = -inv(R) * B' * P12(:,:,k);

end

% Формирование вектора коэффициентов регулятора

% в прямом порядке

size(K_o)

size(K_pr)

i = 1;

len_K = length(K_o(:,1))

for j = len_K : -1 : 1

K_o_p(i,:) = K_o(j,:)

K_pr_p(i,:) = K_pr(j,:);

i = i + 1;

end

% ------------------------------------------------------------------------%

% Построение графика переменных коэффициентов регулятора обратной связи

% в прямом времени

figure(2)

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(3)

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;

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

% Нахождение отслеживаемого сигнала

X_o(:,1) = X_o_0;

h = 0.01;

for k = 1 : len_K

X_o(:, k+1) = X_o(:, k) + h * A_o * X_o(:, k);

end

X_o(:, k+1) = [];

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

for k = 1 : len_K

A_(:,:,k) = A + B * K_o_p(k,:);

end

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,:) * X_o(:,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_o_p(k,:) * X(:,k) + K_pr_p(k,:) * X_o(:,k);

end

% ------------------------------------------------------------------------%

% Построение u(t) и X(t)

figure(4);

plot(time_X, u, 'r-', 'LineWidth', 2)

title ('u(t)');

xlabel('t')

hl=legend('u(t) - управление',0);

set(hl,'FontName','Courier');

grid on

figure(5);

plot(time_X, X(1,:),'r-', time_X, X_o(1,:), 'LineWidth', 2)

hold on

title ('x_1(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

figure(6);

plot(time_X, X(2,:),'r-', time_X, X_o(2,:), 'LineWidth', 2)

title ('x_2(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

figure(7);

plot(time_X, X(3,:),'r-', time_X, X_o(3,:), 'LineWidth', 2)

title ('x_3(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

figure(8);

plot(time_X, X(4,:),'r-', time_X, X_o(4,:), 'LineWidth', 2)

title ('x_4(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

figure(9);

plot(time_X, X(5,:),'r-', time_X, X_o(5,:), 'LineWidth', 2)

title ('x_5(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

AKOR_slegenie_na_konech_interval_II_podxod.m

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 = 45;

h = 0.01;

H = 0.8;

% ------------------------------------------------------------------------%

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_nach = zeros(poryadok, poryadok);%+ones(poryadok, poryadok);

% ------------------------------------------------------------------------%

% Решение уравнения Риккати методом обратного интегрирования

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)

Zadayushee_Vozdeistvie_Discrete_Revers_Modern(h, 0, Time);

% ------------------------------------------------------------------------%

load Zadayushee_Vozdeistvie_Discrete_Revers X_o_discrete_rev

% ------------------------------------------------------------------------%

size(X_o_discrete_rev);

% Нахождение q(t)

for i = 1 : poryadok

qq = -P_nach(:,:,1) * X_o_discrete_rev(i,1);

q(i,1) = qq(i,1);

end

% Интегрирование q(t) в обратном времени

for k = 1 : N_str

q(:, k+1) = q(:, k) - h * ((P2(:,:,k)*B*inv(R)*B'-A') * q(:, k) + Q*X_o_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,:);

X_o_discrete(:,i) = X_o_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));

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, X_o_discrete(1,:), time_X, X_o_discrete(1,:)-0.8,'LineWidth', 2)

hold on

title ('x_1(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон', 'уровень',0);

set(hl,'FontName','Courier');

grid on

figure(7);

plot(time_X, X(2,:),'r-', time_X, X_o_discrete(2,:), 'LineWidth', 2)

title ('x_2(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

figure(8);

plot(time_X, X(3,:),'r-', time_X, X_o_discrete(3,:), 'LineWidth', 2)

title ('x_3(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

figure(9);

plot(time_X, X(4,:),'r-', time_X, X_o_discrete(4,:), 'LineWidth', 2)

title ('x_4(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

figure(10);

plot(time_X, X(5,:),'r-', time_X, X_o_discrete(5,:), 'LineWidth', 2)

title ('x_5(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

AKOR_slegenie_so_skolz_intervalami_Modern.m

function AKOR_slegenie_so_skolz_intervalami_Modern

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 = 45;

Kolvo_intervalov = 3;

h = 0.01;

H = 0.8;

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

% Получение 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+13;

% Q(2,2) = Q(2,2)*1e+10;

% Q(3,3) = Q(3,3)*1e+8;

% Q(4,4) = Q(4,4)*1e+5;

% Q(5,5) = Q(5,5)*1e+2;

R(1,1) = R(1,1);

% ------------------------------------------------------------------------%

% ------------------Скользящие интервалы----------------------------------%

shag = Time/Kolvo_intervalov;

Time1 = shag

Time2 = 2*shag

Time3 = Time

% ------------------------------------------------------------------------%

P_nach = zeros(poryadok, poryadok);%+ones(poryadok, poryadok);

% ------------------------------------------------------------------------%

% Решение уравнения Риккати методом обратного интегрирования

P = Solve_Riccati_Method_Revers_Integr(A,B,Q,R,Time1,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)

% ------------------------------------------------------------------------%

% Нахождение переменных коэффициентов регулятора

for k = 1 : N_str

K_o(k, :) = -inv(R) * B' * P2(:,:,k);

K_pr(k, :) = -inv(R) * B';

end

% ------------------------------------------------------------------------%

tic

% 1 интервал

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Solve_Interval(P_nach, N_str, h, P2, A,B,Q,R, 0, Time1, X_0, poryadok, K_o, K_pr);

load Solve_Interval time_X X u X_o_discrete

time_X1 = time_X;

X1 = X;

u1 = u;

X_o_discrete1 = X_o_discrete;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 2 интервал

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Solve_Interval(P_nach, N_str, h, P2, A,B,Q,R, Time1, Time2, X1(:,N_str), poryadok, K_o, K_pr);

load Solve_Interval time_X X u X_o_discrete

time_X2 = time_X;

X2 = X;

u2 = u;

X_o_discrete2 = X_o_discrete;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 3 интервал

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Solve_Interval(P_nach, N_str, h, P2, A,B,Q,R, Time2, Time3, X2(:,N_str), poryadok, K_o, K_pr);

load Solve_Interval time_X X u X_o_discrete

time_X3 = time_X;

X3 = X;

u3 = u;

X_o_discrete3 = X_o_discrete;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

toc

% ------------------------------------------------------------------------%

% Объединение интервалов

time_X = [time_X1 time_X2 time_X3];

u = [u1 u2 u3];

X = [X1 X2 X3];

X_o_discrete = [X_o_discrete1 X_o_discrete2 X_o_discrete3];

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

% Построение u(t) и X(t)

figure(3);

plot(time_X, u, 'r-','LineWidth', 2);

title ('u(t)');

xlabel('t')

hl=legend('u(t) - управление',0);

set(hl,'FontName','Courier');

grid on

figure(4);

plot(time_X, X(1,:),'r-', time_X, X_o_discrete(1,:), time_X, X_o_discrete(1,:)-0.8,'LineWidth', 2)

hold on

title ('x_1(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон', 'уровень',0);

set(hl,'FontName','Courier');

grid on

figure(5);

plot(time_X, X(2,:),'r-', time_X, X_o_discrete(2,:), 'LineWidth', 2)

title ('x_2(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

figure(6);

plot(time_X, X(3,:),'r-', time_X, X_o_discrete(3,:), 'LineWidth', 2)

title ('x_3(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

figure(7);

plot(time_X, X(4,:),'r-', time_X, X_o_discrete(4,:), 'LineWidth', 2)

title ('x_4(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

figure(8);

plot(time_X, X(5,:),'r-', time_X, X_o_discrete(5,:), 'LineWidth', 2)

title ('x_5(t)');

xlabel('t');

hl=legend('X(t) - слежение','X_o(t) - эталон',0);

set(hl,'FontName','Courier');

grid on

function Solve_Interval(P_nach, N_str, h, P2, A,B,Q,R, T_nach, T_konech, X_0, poryadok, K_o, K_pr)

Zadayushee_Vozdeistvie_Discrete_Revers_Modern(h, T_nach, T_konech);

load Zadayushee_Vozdeistvie_Discrete_Revers X_o_discrete_rev

% ------------------------------------------------------------------------%

% Нахождение q(t)

for i = 1 : poryadok

qq = -P_nach(:,:,1) * X_o_discrete_rev(i,1);

q(i,1) = qq(i,1);

end

% Интегрирование q(t) в обратном времени

for k = 1 : N_str

q(:, k+1) = q(:, k) - h * ((P2(:,:,k)*B*inv(R)*B'-A') * q(:, k) + Q*X_o_discrete_rev(:,k));

end

q(:, k+1) = [];

size_q = size(q)

% ------------------------------------------------------------------------%

% Формирование вектора коэффициентов регулятора, значений задающего

% воздействия, значений вспомогательной функции в прямом порядке

K_pr_p = K_pr;

i = 1;

for j = N_str : -1 : 1

K_o_p(i,:) = K_o(j,:);

X_o_discrete(:,i) = X_o_discrete_rev(:,j);

q_pr(:, i) = q(:, j);

i = i + 1;

end

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

for k = 1 : N_str

A_(:,:,k) = A + B * K_o_p(k,:);

end

size_A_ = size(A_)

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

% Нахождение фазовых координат

X(:,1) = X_0;

time_X(1) = T_nach;

for k = 1 : N_str

X(:, k+1) = X(:, k) + h * (A_(:,:,k) * X(:, k) + B * K_pr_p(k,:) * q_pr(:,k));

time_X(k+1) = time_X(k) + h;

end

X(:, k+1) = [];

time_X(k+1) = [];

size_X = size(X)

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

% Нахождение управления

for k = 1 : N_str

u(k) = K_o_p(k,:) * X(:,k) + K_pr_p(k,:) * q_pr(:,k);

end

size_u = size(u)

save Solve_Interval time_X X u X_o_discrete

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Sintez_nablyud_polnogo_poryadka.m

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 = 10;

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

% Получение 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);

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);

% ------------------------------------------------------------------------%

P_nach = zeros(poryadok, poryadok);%+ones(poryadok, poryadok);

% ------------------------------------------------------------------------%

% Решение уравнения Риккати методом обратного интегрирования

P1 = 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

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;

% ------------------------------------------------------------------------%

% Нахождение коэффициентов регулятора

disp('Коэффициенты регулятора:')

K = -inv(R) * B' * P1

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

A_ = A + B * K;

% ------------------------------------------------------------------------%

% Нахождение фазовых координат

X(:,1) = X_0;

h = 0.01;

time_X(1) = 0;

for k = 1 : N_str

X(:, k+1) = X(:, k) + h * A_ * X(:, k);

time_X(k+1) = time_X(k) + h;

end

X(:, k+1) = [];

time_X(k+1) = [];

% ------------------------------------------------------------------------%

% Нахождение управления

for k = 1 : N_str

u(k) = K * X(:,k);

end

% ------------------------------------------------------------------------%

% Нахождение коэффициентов наблюдателя

M_n = [C' A'*C' (A^2)'*C' (A^3)'*C' (A^4)'*C']

rank_M_n = rank(M_n)

A_r = A_

disp('Спектр матрицы регулятора:')

spektr_A_r = eig(A_r)

koeff = 1;

min_lyamda_A_r = min(real(spektr_A_r))

% lyamda = min_lyamda_A_r * koeff;

lyamda = -5;

disp('Спектр матрицы наблюдателя эталонный:')

lyamda_A_n = [lyamda - koeff * 4; lyamda - koeff * 3; lyamda - koeff * 2;...

lyamda - koeff; lyamda]'

syms k_n1 k_n2 k_n3 k_n4 k_n5 lyam

K_n = [k_n1; k_n2; k_n3; k_n4; k_n5];

Koeff_poly_n_etalon = poly(lyamda_A_n)

disp('Характеристический полином наблюдателя эталонный:')

poly_n_etalon = poly2sym(Koeff_poly_n_etalon, lyam)

disp('Характеристический полином наблюдателя реальный:')

poly_n_real = collect(expand(simplify(det(lyam*eye(poryadok) - (A - K_n*C)))),lyam)

raznost_poly = collect(poly_n_etalon-poly_n_real,lyam)

for i = 1 : poryadok

Koeff_raznost_poly(i) = subs(diff(raznost_poly,poryadok-i,lyam)/factorial(poryadok-i),lyam,0);

end

Koeff_raznost_poly

[Kn1 Kn2 Kn3 Kn4 Kn5]= solve(Koeff_raznost_poly(5), Koeff_raznost_poly(4),...

Koeff_raznost_poly(3), Koeff_raznost_poly(2), Koeff_raznost_poly(1), ...

k_n1, k_n2, k_n3, k_n4, k_n5)

Kn = [Kn1; Kn2; Kn3; Kn4; Kn5];

Kn = vpa(Kn,50)

% Проверка

Proverka = solve(det(lyam*eye(poryadok)-(A-Kn*C)))

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

% Нахождение x и x_оценочного

X_ocen_0 = [0 0 0 0 0]';

A_rash = [A B*K;

Kn*C A-Kn*C+B*K]

X_rash_0 = [X_0;X_ocen_0]

X_rash(:,1) = X_rash_0;

for k = 1 : N_str

X_rash(:,k+1) = X_rash(:,k) + h * A_rash * X_rash(:,k);

end

X_rash(:,k+1) = [];

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

% Разделение x и x_оценочного

for i = 1 : poryadok

X_n(i,:) = X_rash(i,:);

end

for i = poryadok + 1 : 2*poryadok

X_n_ocen(i - poryadok,:) = X_rash(i,:);

end

% ------------------------------------------------------------------------%

% ------------------------------------------------------------------------%

% Нахождение управления

for i = 1 : N_str

u_n(i) = K * X_n_ocen(:,i);

end

% Построение u(t) и X(t)

figure(3);

plot(time_X, u, 'r-', time_X, u_n, 'b-', 'LineWidth', 2)

title ('u(t)');

xlabel('t')

hl=legend('управление без наблюдателя','управление c наблюдателем');

set(hl,'FontName','Courier');

grid on

figure(4);

plot(time_X, X(1,:), time_X, X_n(1,:), time_X, X_n_ocen(1,:),'LineWidth', 2)

hold on

title ('x_1(t)');

xlabel('t')

hl=legend('x_1(t) без наблюдателя','x_1(t) c наблюдателем', 'x_о_ц_е_н_1(t)');

set(hl,'FontName','Courier');

grid on

figure(5);

plot(time_X, X(2,:), time_X, X_n(2,:), time_X, X_n_ocen(2,:),'LineWidth', 2)

title ('x_2(t)');

xlabel('t')

hl=legend('x_2(t) без наблюдателя','x_2(t) c наблюдателем', 'x_о_ц_е_н_2(t)');

set(hl,'FontName','Courier');

grid on

figure(6);

plot(time_X, X(3,:), time_X, X_n(3,:), time_X, X_n_ocen(3,:),'LineWidth', 2)

title ('x_3(t)');

xlabel('t')

hl=legend('x_3(t) без наблюдателя','x_3(t) c наблюдателем', 'x_о_ц_е_н_3(t)');

set(hl,'FontName','Courier');

grid on

figure(7);

plot(time_X, X(4,:), time_X, X_n(4,:), time_X, X_n_ocen(4,:),'LineWidth', 2)

title ('x_4(t)');

xlabel('t')

hl=legend('x_4(t) без наблюдателя','x_4(t) c наблюдателем', 'x_о_ц_е_н_4(t)');

set(hl,'FontName','Courier');

grid on

figure(8);

plot(time_X, X(5,:), time_X, X_n(5,:), time_X, X_n_ocen(5,:),'LineWidth', 2)

title ('x_5(t)');

xlabel('t')

hl=legend('x_5(t) без наблюдателя','x_5(t) c наблюдателем', 'x_о_ц_е_н_5(t)');

set(hl,'FontName','Courier');

grid on

Solve_Riccati_Method_Diag.m

% ------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Метод диагонализации для решения алгебраического уравнения Риккати

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function P = Solve_Riccati_Method_Diag(A,B,Q,R)

% Расширенная матрица системы

Z = [A B*inv(R)*B';

Q -A']

% Нахождение собственных векторов и собственных чисел матрицы Z

[V,D] = eig(Z)

% ------------------------------------------------------------------------%

% Построение матрицы S

% Индексы столбцов собственных значений Re(lyamda) > 0

Ind_Re_plus = find(sum(real(D)) > 0);

% Индексы столбцов собственных значений Re(lyamda) < 0

Ind_Re_minus = find(sum(real(D)) < 0);

% Формирование матрицы D в виде Re(lyamda) > 0 -> Re(lyamda) < 0

D1 = sum(D(:, Ind_Re_plus));

D2 = sum(D(:, Ind_Re_minus));

D = [D1 D2];

% Формирование матрицы S в виде Re(lyamda) > 0 -> Re(lyamda) < 0

S1 = V(:, Ind_Re_plus);

S2 = V(:, Ind_Re_minus);

S = [S1 S2];

% Поиск столбцов с комплексными корнями в матрице D

Ind_Complex_D = find(imag(D) ~= 0);

% Формирование конечной матрицы S

for i = 1 : 2 : length(Ind_Complex_D)

S (:, Ind_Complex_D(i) + 1) = imag(S(:, Ind_Complex_D(i)));

S (:, Ind_Complex_D(i)) = real(S(:, Ind_Complex_D(i)));

end

S = S

% ------------------------------------------------------------------------%

poryadok = length(A(1,:));

S12 = S(1 : poryadok, poryadok+1 : 2*poryadok);

S22 = S(poryadok+1 : 2*poryadok, poryadok+1 : 2*poryadok);

% ------------------------------------------------------------------------%

% Вычисление матрицы P

P = -S22 * inv(S12);

Solve_Riccati_Method_Revers_Integr.m

% ------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Решение уравнения Риккати интегрированием в обратном времени

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function P = Solve_Riccati_Method_Revers_Integr(A,B,Q,R,Time,poryadok, P1)

save For_Riccati A B Q R poryadok

% Решение дифференциального уравнения Риккати

P1 = reshape(P1, poryadok^2, 1);

[Time_R, P] = ode45(@Riccati, [Time : -0.01 : 0], P1);

[N_str, N_stolb] = size(P);

% Построение полученного решения

figure(1)

for i = 1 : poryadok^2

plot(Time_R, P(:,i),'-')

hold on

end

% plot(Time_R,P(:,1),'-',Time_R,P(:,2),'-',Time_R,P(:,3),'-',Time_R,P(:,4),'-',Time_R,P(:,5),'-',Time_R,P(:,6),'-',...

% Time_R,P(:,7),'-',Time_R,P(:,8),'-',Time_R,P(:,9),'-',Time_R,P(:,10),'-',Time_R,P(:,11),'-',Time_R,P(:,12),'-',...

% Time_R,P(:,13),'-',Time_R,P(:,14),'-',Time_R,P(:,15),'-',Time_R,P(:,16),'-',Time_R,P(:,17),'-',Time_R,P(:,18),'-',...

% Time_R,P(:,19),'-',Time_R,P(:,20),'-',Time_R,P(:,21),'-',Time_R,P(:,22),'-',Time_R,P(:,23),'-',Time_R,P(:,24),'-',...

% Time_R,P(:,25),'-', 'lineWidth', 2);

grid on;

tit1 = title('Решения уравнения Риккати');

set(tit1,'FontName','Courier');

xlabel('t');

% legend('p_1','p_2','p_3','p_4','p_5','p_6','p_7','p_8','p_9','p_1_0','p_1_1','p_1_2','p_1_3','p_1_4','p_1_5','p_1_6',...

% 'p_1_7','p_1_8','p_1_9','p_2_0','p_2_1','p_2_2','p_2_3','p_2_4','p_2_5');

save Solve_Riccati_Method_Revers_Integr Time_R P N_str

save Solve_Riccati_Method_Revers_Integr_for_slegenie Time_R P N_str

P = reshape(P(N_str,:), poryadok, poryadok);

function dP = Riccati(Time,P)

load For_Riccati A B Q R poryadok

P = reshape(P, poryadok, poryadok);

% Дифференциальное уравнение Риккати

dP = -P*A - A'*P + P*B*inv(R)*B'*P - Q;

dP = reshape(dP, poryadok^2, 1);

Vozmyshyayushee_Vozdeistvie_Discrete_Revers.m

% Получение дискретных значений возмущающего воздействия в обратном времени

% для нахождения вспомогательной функции q(t)

function Vozmyshyayushee_Vozdeistvie_Discrete_Revers(h, T_nach, T_konech)

% ------------------------------------------------------------------------%

% Возмущающее воздействие

A = 1;

w = 4*pi;

k = 1;

RETURN = 1;

while RETURN == 1

disp('Возмущающее воздействие - const: 1')

disp('Возмущающее воздействие - A*sin(w*t): 2')

reply = input('Выберете возмущающее воздействие [1 или 2]: ', 's');

switch reply

case '1'

disp('Возмущающее воздействие - const')

for t = T_konech: -h : T_nach

w_discrete_rev(:, k) = [A + 0 * t; 0; 0; 0; 0];

k = k + 1;

end

RETURN = 2;

case '2'

disp('Возмущающее воздействие - A*sin(w*t)')

for t = T_konech: -h : T_nach

w_discrete_rev(:, k) = [A * sin(w * t); 0; 0; 0; 0];

k = k + 1;

end

RETURN = 2;

otherwise

disp('Неизвестное воздействие.')

RETURN = 1;

end

end

figure(2)

t = T_konech : -h : T_nach;

plot(t, w_discrete_rev(1,:), 'r-', 'LineWidth', 2);

xlabel('t')

tit1 = title('Возмущающее воздействие');

set(tit1,'FontName','Courier');

hl=legend('Возмущающее воздействие',0);

set(hl,'FontName','Courier');

grid on;

save Vozmyshyayushee_Vozdeistvie_Discrete_Revers w_discrete_rev

% ------------------------------------------------------------------------%

Zadayushee_Vozdeistvie_Discrete_Revers_Modern.m

% Получение дискретных значений задающего воздействия в обратном времени

% для нахождения вспомогательной функции q(t)

function Zadayushee_Vozdeistvie_Discrete_Revers_Modern(h, T_nach, T_konech)

% ------------------------------------------------------------------------%

% Задающее воздействие

alfa = 0.2;

beta = 10;

H = 0.8;

k = 1;

for t = T_konech: -h : T_nach