Курсовая работа: Математическая модель в пространстве состояний линейного стационарного объекта управления
grid on
xlabel('t')
ylabel('k(t)')
title('Function k(t)')
%-------------------------------------------------------------------------%
x1=tf([b1 b0],[a4 a3 a2 a1 a0 0]);
ltiview(x1)
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)
%ВСЕ
ПОДСЧИТАНО ВЕРНО!!!
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
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
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
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
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
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
end
end
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
end
t_perevoda = N * h;
end
N_opt
h
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);
Страницы: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17 |
