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

-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

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