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

1. Методы классической и современной теории автоматического управления: Учебник в 5 – и т. Т.4: Теория оптимизации систем автоматического управления / Под ред. Н.Д. Егупова. – М.: Изд-во МГТУ им. Н.Э. Баумана, 2004. – 748 с.

2. Краснощёченко В.И.: Методическое пособие: «Методы теории оптимального управления».

Приложение.

PlotTimeFrHaract.m

clc

clear all

close all

b1 = 9;

b0 = 5;

a4 = 0.1153;

a3 = 1.78;

a2 = 3.92;

a1 = 14.42;

a0 = 8.583;

% syms s w

% W_s_chislit = b1 * s + b0;

% W_s_znamen = s * (a4 * s^4 + a3 * s^3 + a2 * s^2 + a1 * s + a0);

% W_s_obj = W_s_chislit/W_s_znamen;

%A_w = collect(simplify(abs(subs(W_s_obj, s, i*w))))

%----------------------Построение АЧХ-------------------------------------%

figure('Name', '[0,10]');

w = 0 : 0.01 : 10;

A_w = sqrt((b0^2 + b1^2.*w.^2)./((-a1*w.^2+a3*w.^4).^2+(a0*w-a2*w.^3+a4*w.^5).^2));

plot(w,A_w,'k', 'LineWidth', 2);

grid on

xlabel('w')

ylabel('A(w)')

title('Function ACHX(w)')

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

r_ch = roots([b1 b0])

r_zn = roots([a4 a3 a2 a1 a0 0])

%----------------------Построение ФЧХ-------------------------------------%

figure('Name', '[0,100]');

w = 0 : 0.01 : 100;

fi_w = (atan(w/0.5556)-atan(w/0)-atan(w/13.5832)-atan((w-2.7677)/0.5850)...

-atan((w+2.7677)/0.5850) - atan(w/(0.6848)))*180/pi;

plot(w,fi_w, 'k', 'LineWidth', 2);

grid on

xlabel('w')

ylabel('fi(w)')

title('Function FCHX(w)')

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

%----------------------Построение АФЧХ------------------------------------%

figure('Name', '[0,100]');

w = 0 : 0.01 : 100;

A_w = sqrt((b0^2 + b1^2.*w.^2)./((-a1*w.^2+a3*w.^4).^2+(a0*w-a2*w.^3+a4*w.^5).^2));

fi_w = (atan(w/0.5556)-atan(w/0)-atan(w/13.5832)-atan((w-2.7677)/0.5850)...

-atan((w+2.7677)/0.5850) - atan(w/(0.6848)));

polar(fi_w,A_w, 'k');

grid on

xlabel('Re(W(jw))')

ylabel('Im(W(jw))')

title('Function AFCHX(fi_w,A_w)')

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

%----------------------Построение ЛАЧХ------------------------------------%

figure('Name', '[0,100]');

w = -100 : 0.01 : 100;

LA_w = 20*log(sqrt((b0^2 + b1^2.*w.^2)./((-a1*w.^2+a3*w.^4).^2+(a0*w-a2*w.^3+a4*w.^5).^2)));

plot(w,LA_w,'k', 'LineWidth', 2);

grid on

xlabel('w')

ylabel('L(w)')

title('Function L(w)')

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

%----------------------Построение ФАЧХ------------------------------------%

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

%----------------------Построение h(t)------------------------------------%

figure('Name', '[0,50]');

t = 0 : 0.01 : 50;

h_t = 0.0024 * exp(-13.5832.*t) - 0.2175 * exp(-0.6848.*t)...

+ 0.1452 * exp(-0.5850.*t).* cos(2.7677.*t)...

- 0.2217 * exp(-0.5850.*t).* sin(2.7677.*t)...

+ 0.5825 .* t + 0.0699;

plot(t,h_t, 'k', 'LineWidth', 2);

grid on

xlabel('t')

ylabel('h(t)')

title('Function h(t)')

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

%----------------------Построение k(t)------------------------------------%

figure('Name', '[0,50]');

t = 0 : 0.01 : 50;

k_t = - 0.0329 * exp(-13.5832.*t) + 0.1489 * exp(-0.6848.*t)...

- 0.6986 * exp(-0.5850.*t).* cos(2.7677.*t)...

- 0.2721 * exp(-0.5850.*t).* sin(2.7677.*t)...

+ 0.5826;

plot(t,k_t, 'k', 'LineWidth', 2);

grid on

xlabel('t')

ylabel('k(t)')

title('Function k(t)')

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

x1=tf([b1 b0],[a4 a3 a2 a1 a0 0]);

ltiview(x1)

ProstranstvoSostoyanii.m

clc

clear all

%format rational

b1 = 9;

b0 = 5;

a5 = 0.1153;

a4 = 1.78;

a3 = 3.92;

a2 = 14.42;

a1 = 8.583;

a0 = 0;

%1. Матрица Фробениуса

A=[0 1 0 0 0;

0 0 1 0 0;

0 0 0 1 0;

0 0 0 0 1;

0 -a1/a5 -a2/a5 -a3/a5 -a4/a5]

B=[0; 0; 0; 0; 1/a5]

C=[b0 b1 0 0 0]

%Проверка

syms s

W_s = collect(simplify(C*(s.*eye(5)-A)^(-1)*B),s)

pretty(W_s)

%2. Параллельная декомпозиция

b1 = b1/a5;

b0 = b0/a5;

s1 = 0;

s2 = -6615/487;

s3 = -1022/1747 + 4016/1451*i;

s4 = -1022/1747 - 4016/1451*i;

s5 = -415/606;

alfa = real(s3);

beta = imag(s3);

syms s A B C D E

W_s_etal = collect(((b1*s+b0)/((s-s1)*(s-s2)*((s+alfa)^2+beta^2)*(s-s5))),s)

%pretty(W_s_etal)

Slag_1 = simplify(collect(A*(s-s2)*((s+alfa)^2+beta^2)*(s-s5),s));

Slag_2 = simplify(collect(B*(s-s1)*((s+alfa)^2+beta^2)*(s-s5),s));

Slag_3 = simplify(collect(C*(s-s1)*((s+alfa)^2+beta^2)*(s-s2),s));

Slag_4 = simplify(collect((D*s+E)*(s-s1)*(s-s2)*(s-s5),s));

Chislit_W_s =collect(Slag_1 + Slag_2 + Slag_3 + Slag_4,s);

%Решение системы линейных уравнений

MS =double( [1 1 1 1 0;

6753029497/515578134 -513659/1058682 10560977/850789 4210795/295122 1;

77456808434995506239663107/126764366837761533378822144 1874500571398143988939141/260296441145300889894912 -3300780600401725219142291/418364246989311991349248 915075/98374 4210795/295122;

26189071674868424275768861465/253528733675523066757644288 2853037197681682345182805/520592882290601779789824 45476725452203201718998205/418364246989311991349248 0 915075/98374;

6290947020888109571128085025/84509577891841022252548096 0 0 0 0])

PCH = [0; 0; 0; b1; b0];

Koeff = MS^(-1)*PCH

%Проверка

MS*[Koeff(1);Koeff(2);Koeff(3);Koeff(4);Koeff(5)];

Slag_1 = simplify(collect(Koeff(1)*(s-s2)*((s+alfa)^2+beta^2)*(s-s5),s));

Slag_2 = simplify(collect(Koeff(2)*(s-s1)*((s+alfa)^2+beta^2)*(s-s5),s));

Slag_3 = simplify(collect(Koeff(3)*(s-s1)*((s+alfa)^2+beta^2)*(s-s2),s));

Slag_4 = simplify(collect((Koeff(4)*s+Koeff(5))*(s-s1)*(s-s2)*(s-s5),s));

Chislit_W_s =collect((Slag_1 + Slag_2 + Slag_3 + Slag_4),s);

Znamena_W_s = collect((s-s1)*(s-s2)*((s+alfa)^2+beta^2)*(s-s5),s);

W_s = collect(simplify(Koeff(1)/(s-s1)+Koeff(2)/(s-s2)+(Koeff(4)*s+Koeff(5))/((s+alfa)^2+beta^2)+Koeff(3)/(s-s5)),s)

pretty(W_s)

%Расчет матриц состояния

A = [s1 0 0 0 0;

0 s2 0 0 0 ;

0 0 0 1 0;

0 0 -(alfa^2+beta^2) -2*alfa 0;

0 0 0 0 s5]

B = [Koeff(1); Koeff(2); 0; 1; Koeff(3)]

C = [1 1 Koeff(5) Koeff(4) 1]

%Проверка

W_s = collect(simplify(C*(s.*eye(5)-A)^(-1)*B),s)

pretty(W_s)

%ВСЕ ПОДСЧИТАНО ВЕРНО!!!

SimplexMetod2.m

function SimplexMetod2

clc

clear all

close all

format short

% Матрицы системы

A = [0 2;

-3 0];

B = [0; 2];

% Координаты начальной и конечной точки

X_0 = [14; 0];

X_N = [0; 0];

% Ограничение на управление

u_m = -3;

u_p = 5;

eps = 1e-10;% погрешность сравнения с нулем

N = 195;% число шагов

%h = t1/N;% шаг дискретизации

h = 0.0162;

alfa = 1;

a = 0;

b = 0;

%t1 = 796/245;% время перехода в конечное состояние

n = size(A);

n = n(1);% порядок системы

% Нахождение матричного экспоненциала

syms s t

MatrEx = ilaplace((s*eye(n)-A)^(-1));

MatrEx_B = MatrEx*B;

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

F = subs(MatrEx, t, h);

G = double(int(MatrEx_B, t, 0, h));

ФОРМИРОВАНИЕ ЗАДАЧИ БЫСТРОДЕЙСТВИЯ КАК ЗАДАЧИ

ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ

for index = 1 : 1e+10

% Вычисление правой части

PravChast = X_N - F^N * X_0;

% Вычисление произведения F на G

FG = zeros(n, N);% формирование матрицы для хранения данных

for j = 1 : n

for i = 0 : N - 1

fg = F^(N-i-1) * G;

if PravChast(j) < 0

fg = -fg;

end

FG(j, i+1) = fg(j);

end

% Построение z-строки

z_stroka = zeros(1, 4*N+n+2);% формирование матрицы для хранения данных

% Первый элемент z-строки

z_stroka(1) = 1;

% Суммирование правых частей

for j = 1 : n

z_stroka(4*N+n+2) = z_stroka(4*N+n+2) + abs(PravChast(j));

end

% Формирование элементов z-строки между 1-м и последним элементами

%при 2N небазисных переменных, т.е. при управлениях

for i = 2 : 2 : 2 * N

for j = 1 : n

z_stroka(i) = z_stroka(i) + FG(j, i/2);

end

for j = 1 : n

z_stroka(i+1) = z_stroka(i+1) - FG(j, i/2);

end

% Формирование симплекс-таблицы

CT = zeros(n+2*N+1, 4*N+n+2);

% Построение симплекс-таблицы начиная с z-строки

CT(1,:) = z_stroka(1,:);

% Формирование R-строк в симплекс-таблице

for j = 2 : n + 1

% Формирование правой части в R-строках

CT(j, 4*N+n+2) = abs(PravChast(j-1));

% Формирование элементов R-строк между 1-м и последним элементами

%при 2N небазисных переменных, т.е. при управлениях

for i = 2 : 2 : 2 * N

CT(j, i) = FG(j-1, i/2);

CT(j, i+1) = -FG(j-1, i/2);

end

% Формирование S-строк в симплекс-таблице

l = 2;

for j = n + 2 : 2 : n + 2 * N + 1

% Формирование правой части в S-строках

CT(j, 4*N+n+2) = u_p;

CT(j+1, 4*N+n+2) = abs(u_m);

% Формирование элементов S-строк между 1-м и последним элементами

%при 2N небазисных переменных, т.е. при управлениях

CT(j, l : l+1) = [1 -1];

CT(j+1, l : l+1) = [-1 1];

l = l + 2;

end

% Формирование базиса в симплекс-таблице, т.е коэффициентов, стоящих при

%базисных переменных от 2N небазисных переменных до правой части (до 4*N+n+1)

CT(2 : n+2*N+1, 2*N+2 : 4*N+n+1) = eye(n+2*N, n+2*N);

РЕШЕНИЕ ЗАДАЧИ БЫСТРОДЕЙСТВИЯ

% Цикл смены базисных переменных

nn = size(find(CT(1,2:2*N+1) >= eps));

while nn > 0

[znach, N_stolb] = max(CT(1, 2 : 2*N+1));

N_stolb = N_stolb + 1; % т.к. при небазисн. перемен.

PravChast = CT(:, 4*N+n+2);

for j = 2 : n + 2 * N + 1

if CT(j, N_stolb) > 0

PravChast(j) = PravChast(j) / CT(j, N_stolb);

else

PravChast(j) = inf;

end

[znach, N_str] = min(PravChast(2 : n+2*N+1));

N_str = N_str + 1;

% Формирование матрицы перехода B

B = eye(n+2*N+1, n+2*N+1);

B(:, N_str) = CT(:, N_stolb);

% Обращение матрицы B

RE = B(N_str, N_str);

for j = 1 : n + 2 * N + 1

if j == N_str

B(j, N_str) = 1 / RE;

else

B(j, N_str) = -B(j, N_str) / RE;

end

%B = inv(B);

% Получение новой симплекс таблицы

CT = B * CT;

nn = size(find(CT(1,2:2*N+1) >= eps));

end

u = zeros(1,N);

% Формирование управления

for j = 2 : n + 2 * N + 1

for i = 2 : 2 * N + 1

if CT(j, i) >= eps

if mod(i, 2) < eps

u(i/2) = CT(j, 4*N+n+2);

else

u((i-1)/2) = -CT(j, 4*N+n+2);

end

% Формирование x1 и x2

X = zeros(n, N);

X(:, 1) = F * X_0 + G * u(1);

for i = 2 : N

X(:, i) = F * X(:, i-1) + G * u(i);

end

% Объединение с начальными условиями

X1 = [X_0(1) X(1, :)];

X2 = [X_0(2) X(2, :)];

% проверка на окончание выбора количества шагов

XX = [X_0 X];

% Вычисление нормы вектора состояния

normaXX = norm(XX(:,N))

% Вычисление значения переменной R

R = abs(X_N - F^N * X_0) - FG * u';

R = R';

z = sum(R);

% Погрешность приближения к точному решению

pogresh = 0.3;

if (normaXX < pogresh)

N_opt = N;

break;

else

if (z > h)

if a == 1

alfa = ceil(alfa/2);

end

N = N + alfa;

a = 0;

b = 1;

else

if b == 1

alfa = ceil(alfa/2);

end

N = N - alfa;

a = 1;

b = 0;

end

t_perevoda = N * h;

end

N_opt

t_perevoda

ОФОРМЛЕНИЕ ПОЛУЧЕННЫХ РЕЗУЛЬТАТОВ

В ГРАФИЧЕСКОМ ВИДЕ

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

figure(1)

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

plot(t, X1, 'b', 'LineWidth', 2);

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

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

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

grid on

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

figure(2)

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

plot(t, X2, 'b', 'LineWidth', 2);

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

% h1_Tt = vpa(h1_Tt,6)

% h2_Tt = vpa(h2_Tt,6)

% h3_Tt = vpa(h3_Tt,6)

% h4_Tt = vpa(h4_Tt,6)

Страницы: 1, 2, 3, 4, 5, 6, 7

МЕНЮ

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

PlotTimeFrHaract.m

ProstranstvoSostoyanii.m

SimplexMetod2.m

Optimal_L_problem_moments.m

ИНТЕРЕСНОЕ