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 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); 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Популярное: Как распознать напряжение: Говоря о мышечном напряжении, мы в первую очередь имеем в виду мускулы, прикрепленные к костям ... Почему двоичная система счисления так распространена?: Каждая цифра должна быть как-то представлена на физическом носителе... Почему стероиды повышают давление?: Основных причин три... ©2015-2024 megaobuchalka.ru Все материалы представленные на сайте исключительно с целью ознакомления читателями и не преследуют коммерческих целей или нарушение авторских прав. (178)
|
Почему 1285321 студент выбрали МегаОбучалку... Система поиска информации Мобильная версия сайта Удобная навигация Нет шокирующей рекламы |