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

hl=legend('x_2(t)');

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

xlabel('t, cek'); ylabel('x_2(t)');

grid on

% Построение графика x2 = x2(x1);

figure(3)

plot(X1, X2, 'm', 'LineWidth', 2);

hl=legend('x_2 = x_2(x_1)');

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

xlabel('x_1(t)'); ylabel('x_2(x_1(t))');

grid on

% Построение графика u(t)

figure(4)

t = (0 : 1 : length(u)-1) * h;

plot(t, u, 'r', 'LineWidth', 2);

hl=legend('u(t)');

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

xlabel('t, cek'); ylabel('u(t)');

grid on

Optimal_L_problem_moments.m

clc

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

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

RETURN = 1;

while RETURN == 1

disp('L - проблема моментов в пространстве вход-выход: 1')

disp('L - проблема моментов в пространстве состояний : 2')

reply = input('Выберете метод решения [1 или 2]: ', 's');

switch reply

 case '1'

 disp('L - проблема моментов в пространстве вход-выход')

% ------------------------L - проблема моментов---------------------------%

% ----------------------в пространстве вход-выход-------------------------%

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

% Передаточная функция

W_obj_s = 1/(a5*s^5 + a4*s^4 + a3*s^3 + a2*s^2 + a1*s + a0);

% Полюса передаточной функции

polyusa_TF = roots([a5 a4 a3 a2 a1 a0]);

% ИПФ

K_t = simplify (vpa (ilaplace(1 / (a5*s^5 + a4*s^4 + a3*s^3 + ...

 a2*s^2 + a1*s + a0)),50))

% K_t = vpa(K_t,6)

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

% Составление матрицы Вронского

for i = 1 : poryadok

Matrix_Vron (i, 1) = diff (exp (polyusa_TF(1) *t), t, i - 1);

Matrix_Vron (i, 2) = diff (exp (polyusa_TF(2) *t), t, i - 1);

Matrix_Vron (i, 3) = diff (exp (real(polyusa_TF(3))*t) * ...

cos(imag(polyusa_TF(3))*t), t, i - 1);

Matrix_Vron (i, 4) = diff (exp (real(polyusa_TF(4))*t) * ...

sin(imag(polyusa_TF(4))*t), t, i - 1);

Matrix_Vron (i, 5) = diff (exp (polyusa_TF(5) *t), t, i - 1);

end

% Матрица Вронского при t = 0;

Matrix_Vron_t_0 = double(subs(Matrix_Vron,t,0));

% Матрица Вронского при t = T;

T = 3;

Matrix_Vron_t_T = double(subs(Matrix_Vron,t,T));

% vpa(Matrix_Vron_t_0,6)

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

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

% Определение неизвестных коэффициентов C

C_ = inv(Matrix_Vron_t_0) * X_0;

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

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

% Нахождение моментных функций

K_Tt_1 = subs (K_t,t, T - t);

K_Tt = diff (K_t);

K_Tt_2 = subs (K_Tt, t, T - t);

K_Ttt = diff (K_Tt);

K_Tt_3 = subs (K_Ttt, t, T - t);

K_Tttt = diff (K_Ttt);

K_Tt_4 = subs (K_Tttt, t, T - t);

K_Ttttt = diff (K_Tttt);

K_Tt_5 = subs (K_Ttttt, t, T - t);

h1_Tt = K_Tt_1

h2_Tt = K_Tt_2

h3_Tt = K_Tt_3

h4_Tt = K_Tt_4

h5_Tt = K_Tt_5

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

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

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

for i = 1 : poryadok

Matrix_a(i) = X_T(i) - C_' * Matrix_Vron_t_T(i,:)';

end

Matrix_a = Matrix_a'

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

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

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

RETURN = 2;

 case '2'

 disp('L - проблема моментов в пространстве состояний')

% ------------------------L - проблема моментов---------------------------%

% ----------------------в пространстве состояний--------------------------%

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

Matr_Ex_T = subs(MatrEx, t, T);

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

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

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

for i = 1 : poryadok

Matrix_a(i) = X_T(i) - Matr_Ex_T(i,:) * X_0;

end

Matrix_a = Matrix_a'

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

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

% Нахождение моментных функций

Matr_Ex_Tt = subs(MatrEx, t, T - t);

h_Tt = vpa(expand(simplify(Matr_Ex_Tt * B)),50);

h1_Tt = h_Tt(1)

h2_Tt = h_Tt(2)

h3_Tt = h_Tt(3)

h4_Tt = h_Tt(4)

h5_Tt = h_Tt(5)

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

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

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

RETURN = 2;

otherwise

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

 RETURN = 1;

end

end

% h1_Tt = vpa(h1_Tt,6)

% h2_Tt = vpa(h2_Tt,6)

% h3_Tt = vpa(h3_Tt,6)

% h4_Tt = vpa(h4_Tt,6)

% 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 + ...

Страницы: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17