|
Курсовая работа: Математическая модель в пространстве состояний линейного стационарного объекта управления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.mclc 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.mclose 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.mclc 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.mclc 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; % ------------------------------------------------------------------------% |
|
|||||||||||||||||||||||||||||
|
Рефераты бесплатно, реферат бесплатно, сочинения, курсовые работы, реферат, доклады, рефераты, рефераты скачать, рефераты на тему, курсовые, дипломы, научные работы и многое другое. |
||
При использовании материалов - ссылка на сайт обязательна. |