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