Московський державний технічний університет ім. Н.Е. Баумана
Калузький філія
Кафедра "САУ і Електротехніки"
ЕІУ3-КФ
Розрахунково-пояснювальна записка до курсової роботи
на тему:
"МАТЕМАТИЧНЕ МОДЕЛЮВАННЯ В ЗАДАЧАХ РОЗРАХУНКУ ТА ПРОЕКТУВАННЯ СИСТЕМ АВТОМАТИЧНОГО УПРАВЛІННЯ"
по курсу:
Системи аналітичних обчислень
Калуга
Зміст
1 Постановка завдання
2 Аналіз стійкості
3 Рішення диференціального рівняння інтерполяційним методом Адамса
4 Синтез
Висновок
Список літератури
Додаток 1 (Лістинг скрипта для знаходження коренів полінома)
Додаток 2 (Лістинг скрипта для вирішення диференціального рівняння
чисельним методом)
Додаток 3 (Лістинг скриптів для знаходження коефіцієнтів регулятора)
1 Постановка завдання
Потрібно:
1. Виконати аналіз стійкості роботи нескорректированной системи управління.
2. Виконати аналіз функціонування системи
3. Синтезувати регулятор для системи управління.
4. Виконати аналіз роботи скоригованої системи управління.
Структурна схема системи приведена на рис. 1.
Рис. 1. Структурна схема контуру стабілізації кута тангажу
Параметри системи мають таке значення:
Вимоги до системи:
2 Аналіз стійкості
Виконаємо аналіз нескорректированной системи з використанням критеріїв Михайлова і Гурвіца.
Знайдемо передатну функцію всієї системи
Складемо матрицю Гурвіца
a 0 = 1; a 1 = 7.4; a 2 = 19; a 3 = 10;
За критерієм Гурвіца для того, щоб система була стійка необхідно і достатньо, щоб всі визначники на головній діагоналі були більше нуля Знайдемо всі мінори на головній діагоналі:
Система стійка.
Критерій Михайлова:
З умови
Одержуємо, що система стійка.
Побудуємо годограф розімкнутої системи і знайдемо запас стійкості.
На рис. 2 наведено графік АФЧХ розімкнутої системи і одинична коло.
Рис. 2.Годограф АФЧХ розімкнутої системи
За рис. 2 легко визначити запас стійкості замкнутої системи.
Знаходження коренів характеристичного рівняння методом градієнтів.
Знайдемо коріння передавальної функції за допомогою методу градієнтів.
Робоча формула використовуваного методу має такий вигляд
де
і
вектори невідомих на кроці k +1 і k.
- Транспонована матриця Якобі, обчислена на кроці k.
Невязка на кроці k
Кроковий множник
Знаходимо полюса для передавальної функції, що має вигляд
Текст програми приведений у додатку 1.
Результат наведено на рис.3
Рис. 3. Приклад знаходження полюсів ПФ W (s)
Аналітичні вирази для перехідної та імпульсної перехідної функцій, АЧХ, ФЧХ, АФЧХ
Знайдемо імпульсну перехідну функцію.
Графік k (t) наведено на рис. 4.
Рис. 4. Графік імпульсної перехідної функції.
Знайдемо перехідну функцію.
Графік h (t) наведено на рис. 5.
Рис. 5. Графік перехідної функції.
Знайдемо амплітудно-частотну характеристику.
Графік АЧХ наведено на рис. 6.
Рис. 6. Графік АЧХ
Знайдемо ФЧХ:
Графік ФЧХ наведено на рис. 7.
Рис. 7. Графік ФЧХ
Знайдемо АФЧХ.
Графік ФЧХ наведено на рис. 8.
Рис. 8. Графік АФЧХ
Висновок: Система є стійкою, перерегулювання дорівнює 0, час керування приблизно дорівнює 5с.
3 Рішення диференціального рівняння інтерполяційним методом Адамса
Так як ДУ заданої системи має третій порядок, то його необхідно звести до системи рівнянь, кожне з яких повинна мати перший порядок, тобто має місце нормальна форма Коші:
Запишемо нормальну форму Коші в наступному вигляді:
Наведемо рівняння до нормального формі Коші:
Інтерполяційний метод Адамса 3:
, Точність
Для того, щоб використовувати цей неявний метод, потрібно знати
Отримаємо методом Ейлера:
точність
Для отримання точності на першому кроці, візьмемо
Текст програми знаходиться у додатку 2.
Результати роботи програми при h рівних 0.5, 0.2, 0.01 наведені на рис. 9.
Рис. 9. Відгуки на одиничне поетапне вплив
4 Синтез
Введемо в пряму ланцюг ПІД регулятор, а у зворотний ПД.
Вид скоригованої системи наведено на рис. 10.
Рис.10. Структурна схема скоректованої системи
Знайдемо передавальну функцію системи
Передавальна функція розімкнутої ланцюга має вигляд:
Передавальна функція розімкнутої ланцюга має вигляд:
Для вирішення завдання синтезу необхідно знайти параметра регулятора, при яких реальний вихідний сигнал, що є реакцією на одиничне поетапне вплив, буде близький до заданого еталонного сигналу.
В якості еталонного вихідного сигналу використовуємо наступний сигнал:
,
Коефіцієнт знаходимо за такою формулою:
Знайдемо параметри регулятора методом квадратичної апроксимації.
Робоча формула методу має вигляд:
Де,
градієнт функції.
матриця Гессе функції
знаходимо за допомогою методу Золотого перетину.
Текст програми знаходиться в додатку 3.
Результат роботи програми приведений на рис. 11.
Рис. 11. Приклад отримання коефіцієнтів регулятора.
Перехідна функція скоригованої системи зображена на рис. 12.
Час управління скоригованої системи виходячи з графіка приблизно дорівнює 2.4с.
Рис. 12. Порівняння еталонною і реальною перехідних функцій
Висновок
У цій роботі був синтезований регулятор САУ, знайдені його параметри чисельним методом. Також було вирішено диференціальне рівняння неявним чисельним методом.
Список літератури
1. Методи класичної та сучасної теорії автоматичного керування: Підручник у 5-ти т.; 2-е вид., Перераб. і доп. Т.3: Синтез регуляторів систем автоматичного керування / За редакцією К.А. Пупкова і Н.Д. Егупова. - М.: Видавництво МГТУ ім. Н.Е. Баумана, 2004. - 616с.; Мул.
2. Н.Д. Єгупов, Ю.П. Корнюшин, Ю.І. Мишляев. Навчальний посібник з виконання курсового проектування з дисципліни «Системи аналітичних обчислень» для студентів спеціальності 160403 «Системи управління літальними апаратами»
Додаток 1
Метод градієнтів
function M _ Gradientov
clc
% Вирішимо рівняння s ^ 3 +7,4 * s ^ 2 +19 * s +10 = 0
e = 10 ^ -4;
s = 0;
A1 = 1;
A2 = 7.4;
A3 = 19;
A4 = 10;
r0 = 1;
i = 0;% кількість ітерацій
while abs (r0)> e
i = i +1;
s0 = s;
r0 = A1 * s ^ 3 + A2 * s ^ 2 + A3 * s + A4;% невязка
Ar = (A1 + A2 + A3) * r0;
AAr = (A1 ^ 2 + A2 ^ 2 + A3 ^ 3) * r0;
m = r0 * AAr / AAr ^ 2;
s = s0-m * Ar;
end
S 1 = s;% Знайшли вешественний корінь
Тепер вирішуємо ураненіе: A1 * s ^ 2 + (A2 + A1 * S1) * s + (A3 + A2 * S1 + A1 * s ^ 2) = 0
% Коріння комплексні
D = (A2 + A1 * S1) ^ 2-4 * A1 * (A3 + A2 * S1 + A1 * s ^ 2);
S2 = (- (A2 + A1 * S1) + sqrt (D)) / 2 * A1;
S3 = (- (A2 + A1 * S1)-sqrt (D)) / 2 * A1;
disp (S1);
disp (S2);
disp (S3);
disp ('Кількість ітерацій'); disp (i);
Додаток 2
Інтерполяційний метод Адамса
function Int_Adams_3
clc
% Час перехідного процесу
T = 10;
%-----------------%
% Матриця А (X '= AX + BY)
A = [0 1 0;
0 0 1;
-10 -19 -7.4];
% Матриця B
B = [0 5 10] ';
Y = [0 0 1] ';
k = 1;
% Початкові умови
X (1,1:3) = [0 0 0];
I = [1 0 0; 0 1 0, 0 0 1];
while (k <= 3)
% Крок
if (k == 1) h = 0.1; end;
if (k == 2) h = 1; end;
if (k == 3) h = 0.01; end;
%---------------------------%
n = 1;
F (1,1:3) = (A * (X (1,1:3)) '+ B. * Y)';
X (n +1,1:3) = (X (n, 1:3) '+ h/10 * (F (n, 1:3 ))')';% Метод Ейлера
n = n +1;
while (n <= T / h)
F (n, 1:3) = (A * (X (n, 1:3)) '+ B. * Y)';
X (n +1,1:3 )=((( I-5 * h/12 * A) ^ -1) * (X (n, 1:3) '+ h/12 * (5 * B. * Y +8 * (F (n, 1:3 ))'-( F (n-1, 1:3 ))')))';
n = n +1;
end
t = 0: h: 10;
% K = t / h +1;
i = 1;
while (i <= n)
if (k == 1) t1 = t; x1 (i) = X (i, 1); Xa1 = 1-0.9202 * exp (-0.6983 * t) -0.4636 * exp (-3.3508 * t) .* cos ( 1.7584 * t +4.382) +0.2433 * exp (-3.3508 * t) .* sin (1.7584 * t +4.382); end;
if (k == 2) t2 = t; x2 (i) = X (i, 1); Xa2 = 1-0.9202 * exp (-0.6983 * t) -0.4636 * exp (-3.3508 * t) .* cos ( 1.7584 * t +4.382) +0.2433 * exp (-3.3508 * t) .* sin (1.7584 * t +4.382); end;
if (k == 3) t3 = t; x3 (i) = X (i, 1); Xa3 = 1-0.9202 * exp (-0.6983 * t) -0.4636 * exp (-3.3508 * t) .* cos ( 1.7584 * t +4.382) +0.2433 * exp (-3.3508 * t) .* sin (1.7584 * t +4.382);
end;
i = i +1;
end
k = k +1;
end
t = 0:0.01:10;
Xa = 1-0.9202 * exp (-0.6983 * t) -0.4636 * exp (-3.3508 * t) .* cos (1.7584 * t +4.382) +0.2433 * exp (-3.3508 * t) .* sin (1.7584 * t +4.382);
plot (t, Xa, t1, x1, t1, (Xa1-x1), t2, x2, t2, (Xa2-x2), t3, x3, t3, (Xa3-x3)), grid on
Пріложеніе3
Оптимізація методом квадратичної апроксимації
function minK
% Задамо точність і крок
eps = 0.1;
h = 0.1;
% Визначимо матрицю K = [Kp, Kd, Ki, Kp 2, Kd 2] ';
T = 4;
K0 = [26 6 50 1 0.2] ';
% Знайдемо J0
J0 = Xr5 (26, 6, 50, 1, 0.2, T);
%------------------------
% Шукаємо матрицю G
a = Xr 5 (K 0 (1), K 0 (2), K 0 (3), K 0 (4), K 0 (5), T);
g11 = (Xr5 (K0 (1) +2 * h, K0 (2), K0 (3), K0 (4), K0 (5), T) -2 * Xr5 (K0 (1) + h, K0 ( 2), K0 (3), K0 (4), K0 (5), T) + a) / h ^ 2;
g12 = (Xr5 (K0 (1) + h, K0 (2) + h, K0 (3), K0 (4), K0 (5), T)-Xr5 (K0 (1) + h, K0 (2) , K0 (3), K0 (4), K0 (5), T)-Xr5 (K0 (1), K0 (2) + h, K0 (3), K0 (4), K0 (5), T) + a) / h ^ 2;
g21 = g12;
g13 = (Xr5 (K0 (1) + h, K0 (2), K0 (3) + h, K0 (4), K0 (5), T)-Xr5 (K0 (1) + h, K0 (2) , K0 (3), K0 (4), K0 (5), T)-Xr5 (K0 (1), K0 (2), K0 (3) + h, K0 (4), K0 (5), T) + a) / h ^ 2;
g31 = g13;
g14 = (Xr5 (K0 (1) + h, K0 (2), K0 (3), K0 (4) + h, K0 (5), T)-Xr5 (K0 (1) + h, K0 (2) , K0 (3), K0 (4), K0 (5), T)-Xr5 (K0 (1), K0 (2), K0 (3), K0 (4) + h, K0 (5), T) + a) / h ^ 2;
g41 = g14;
g15 = (Xr5 (K0 (1) + h, K0 (2), K0 (3), K0 (4), K0 (5) + h, T)-Xr5 (K0 (1) + h, K0 (2) , K0 (3), K0 (4), K0 (5), T)-Xr5 (K0 (1), K0 (2), K0 (3), K0 (4), K0 (5) + h, T) + a) / h ^ 2;
g51 = g15;
g22 = (Xr5 (K0 (1), K0 (2) +2 * h, K0 (3), K0 (4), K0 (5), T) -2 * Xr5 (K0 (1), K0 (2) + h, K0 (3), K0 (4), K0 (5), T) + a) / h ^ 2;
g23 = (Xr5 (K0 (1), K0 (2) + h, K0 (3) + h, K0 (4), K0 (5), T)-Xr5 (K0 (1), K0 (2) + h , K0 (3), K0 (4), K0 (5), T)-Xr5 (K0 (1), K0 (2), K0 (3) + h, K0 (4), K0 (5), T) + a) / h ^ 2;
g32 = g23;
g24 = (Xr5 (K0 (1), K0 (2) + h, K0 (3), K0 (4) + h, K0 (5), T)-Xr5 (K0 (1), K0 (2) + h , K0 (3), K0 (4), K0 (5), T)-Xr5 (K0 (1), K0 (2), K0 (3), K0 (4) + h, K0 (5), T) + a) / h ^ 2;
g42 = g24;
g25 = (Xr5 (K0 (1), K0 (2) + h, K0 (3), K0 (4), K0 (5) + h, T)-Xr5 (K0 (1), K0 (2) + h , K0 (3), K0 (4), K0 (5), T)-Xr5 (K0 (1), K0 (2), K0 (3), K0 (4), K0 (5) + h, T) + a) / h ^ 2;
g52 = g25;
g33 = (Xr5 (K0 (1), K0 (2), K0 (3) +2 * h, K0 (4), K0 (5), T) -2 * Xr5 (K0 (1), K0 (2) , K0 (3) + h, K0 (4), K0 (5), T) + a) / h ^ 2;
g34 = (Xr5 (K0 (1), K0 (2), K0 (3) + h, K0 (4) + h, K0 (5), T)-Xr5 (K0 (1), K0 (2), K0 (3) + h, K0 (4), K0 (5), T)-Xr5 (K0 (1), K0 (2), K0 (3), K0 (4) + h, K0 (5), T) + a) / h ^ 2;
g43 = g34;
g35 = (Xr5 (K0 (1), K0 (2), K0 (3) + h, K0 (4), K0 (5) + h, T)-Xr5 (K0 (1), K0 (2), K0 (3) + h, K0 (4), K0 (5), T)-Xr5 (K0 (1), K0 (2), K0 (3), K0 (4), K0 (5) + h, T) + a) / h ^ 2;
g53 = g35;
g44 = (Xr5 (K0 (1), K0 (2), K0 (3), K0 (4) +2 * h, K0 (5), T) -2 * Xr5 (K0 (1), K0 (2) , K0 (3), K0 (4) + h, K0 (5), T) + a) / h ^ 2;
g45 = (Xr5 (K0 (1), K0 (2), K0 (3), K0 (4) + h, K0 (5) + h, T)-Xr5 (K0 (1), K0 (2), K0 (3), K0 (4) + h, K0 (5), T)-Xr5 (K0 (1), K0 (2), K0 (3), K0 (4), K0 (5) + h, T) + a) / h ^ 2;
g54 = g45;
g55 = (Xr5 (K0 (1), K0 (2), K0 (3), K0 (4), K0 (5) +2 * h, T) -2 * Xr5 (K0 (1), K0 (2) , K0 (3), K0 (4), K0 (5) + h, T) + a) / h ^ 2;
G = [g11, g12, g13, g14, g15; g21, g22, g23, g24, g25; g31, g32, g33, g34, g35; g41, g42, g43, g44, g45; g51, g52, g53, g54 , g55;];
% G1 = G. ^ -1;
G1 = inv (G);
% Побудуємо градієнт
gr1 = (Xr5 (K0 (1) + h, K0 (2), K0 (3), K0 (4), K0 (5), T)-a) / h;
gr2 = (Xr5 (K0 (1), K0 (2) + h, K0 (3), K0 (4), K0 (5), T)-a) / h;
gr3 = (Xr5 (K0 (1), K0 (2), K0 (3) + h, K0 (4), K0 (5), T)-a) / h;
gr4 = (Xr5 (K0 (1), K0 (2), K0 (3), K0 (4) + h, K0 (5), T)-a) / h;
gr5 = (Xr5 (K0 (1), K0 (2), K0 (3), K0 (4), K0 (5) + h, T)-a) / h;
grad = [gr1 gr2 gr3 gr4 gr5] ';
L = lambdamin (K0, G1, grad);
K = K0 + L * G1 * grad;
G10 = G1;
grad0 = grad;
J = Xr5 (K (1), K (2), K (3), K (4), K (5), T);
% Квадратична апроксимація: X (i +1) = X (i)-L (i) G ^ -1 (i) GRAD (x (i))
while (J0> J)
J0 = J;
% Шукаємо матрицю G
a = Xr5 (K (1), K (2), K (3), K (4), K (5), T);
g11 = (Xr5 (K (1) +2 * h, K (2), K (3), K (4), K (5), T) -2 * Xr5 (K (1) + h, K ( 2), K (3), K (4), K (5), T) + a) / h ^ 2;
g12 = (Xr5 (K (1) + h, K (2) + h, K (3), K (4), K (5), T)-Xr5 (K (1) + h, K (2) , K (3), K (4), K (5), T)-Xr5 (K (1), K (2) + h, K (3), K (4), K (5), T) + a) / h ^ 2;
g21 = g12;
g13 = (Xr5 (K (1) + h, K (2), K (3) + h, K (4), K (5), T)-Xr5 (K (1) + h, K (2) , K (3), K (4), K (5), T)-Xr5 (K (1), K (2), K (3) + h, K (4), K (5), T) + a) / h ^ 2;
g31 = g13;
g14 = (Xr5 (K (1) + h, K (2), K (3), K (4) + h, K (5), T)-Xr5 (K (1) + h, K (2) , K (3), K (4), K (5), T)-Xr5 (K (1), K (2), K (3), K (4) + h, K (5), T) + a) / h ^ 2;
g41 = g14;
g15 = (Xr5 (K (1) + h, K (2), K (3), K (4), K (5) + h, T)-Xr5 (K (1) + h, K (2) , K (3), K (4), K (5), T)-Xr5 (K (1), K (2), K (3), K (4), K (5) + h, T) + a) / h ^ 2;
g51 = g15;
g22 = (Xr5 (K (1), K (2) +2 * h, K (3), K (4), K (5), T) -2 * Xr5 (K (1), K (2) + h, K (3), K (4), K (5), T) + a) / h ^ 2;
g23 = (Xr5 (K (1), K (2) + h, K (3) + h, K (4), K (5), T)-Xr5 (K (1), K (2) + h , K (3), K (4), K (5), T)-Xr5 (K (1), K (2), K (3) + h, K (4), K (5), T) + a) / h ^ 2;
g32 = g23;
g24 = (Xr5 (K (1), K (2) + h, K (3), K (4) + h, K (5), T)-Xr5 (K (1), K (2) + h , K (3), K (4), K (5), T)-Xr5 (K (1), K (2), K (3), K (4) + h, K (5), T) + a) / h ^ 2;
g42 = g24;
g25 = (Xr5 (K (1), K (2) + h, K (3), K (4), K (5) + h, T)-Xr5 (K (1), K (2) + h , K (3), K (4), K (5), T)-Xr5 (K (1), K (2), K (3), K (4), K (5) + h, T) + a) / h ^ 2;
g52 = g25;
g33 = (Xr5 (K (1), K (2), K (3) +2 * h, K (4), K (5), T) -2 * Xr5 (K (1), K (2) , K (3) + h, K (4), K (5), T) + a) / h ^ 2;
g34 = (Xr5 (K (1), K (2), K (3) + h, K (4) + h, K (5), T)-Xr5 (K (1), K (2), K (3) + h, K (4), K (5), T)-Xr5 (K (1), K (2), K (3), K (4) + h, K (5), T) + a) / h ^ 2;
g43 = g34;
g35 = (Xr5 (K (1), K (2), K (3) + h, K (4), K (5) + h, T)-Xr5 (K (1), K (2), K (3) + h, K (4), K (5), T)-Xr5 (K (1), K (2), K (3), K (4), K (5) + h, T) + a) / h ^ 2;
g53 = g35;
g44 = (Xr5 (K (1), K (2), K (3), K (4) +2 * h, K (5), T) -2 * Xr5 (K (1), K (2) , K (3), K (4) + h, K (5), T) + a) / h ^ 2;
g45 = (Xr5 (K (1), K (2), K (3), K (4) + h, K (5) + h, T)-Xr5 (K (1), K (2), K (3), K (4) + h, K (5), T)-Xr5 (K (1), K (2), K (3), K (4), K (5) + h, T) + a) / h ^ 2;
g54 = g45;
g55 = (Xr5 (K (1), K (2), K (3), K (4), K (5) +2 * h, T) -2 * Xr5 (K (1), K (2) , K (3), K (4), K (5) + h, T) + a) / h ^ 2;
G = [g11, g12, g13, g14, g15; g21, g22, g23, g24, g25; g31, g32, g33, g34, g35; g41, g42, g43, g44, g45; g51, g52, g53, g54 , g55;];
% G1 = G. ^ -1;
G1 = inv (G);
% Побудова градієнта
gr1 = (Xr5 (K (1) + h, K (2), K (3), K (4), K (5), T)-a) / h;
gr2 = (Xr5 (K (1), K (2) + h, K (3), K (4), K (5), T)-a) / h;
gr3 = (Xr5 (K (1), K (2), K (3) + h, K (4), K (5), T)-a) / h;
gr4 = (Xr5 (K (1), K (2), K (3), K (4) + h, K (5), T)-a) / h;
gr5 = (Xr5 (K (1), K (2), K (3), K (4), K (5) + h, T)-a) / h;
grad = [gr1 gr2 gr3 gr4 gr5] ';
if (Xr5 (K (1), K (2), K (3), K (4), K (5), T)> Xr5 (K0 (1), K0 (2), K0 (3), K0 (4), K0 (5), T))
L = lambdamin (K0, G10, grad0);
end
K0 = K;
G10 = G1;
grad0 = grad;
K = K0 + L * G1 * grad;
J = Xr5 (K (1), K (2), K (3), K (4), K (5), T);
end
disp (K0);
disp (J0);
Метод Золотого Перетини
function L = lambdamin (K, G1, grad)
Xzs = (-1 + sqrt (5)) / 2;% золотий перетин
a = 0;
b = 1;
while (abs (ba)> 0.01)
x1 = a + (ba) * Xzs;
x2 = b-(ba) * Xzs;
F1 = K-x1 * G1 * grad;
F2 = K-x2 * G1 * grad;
if ((x1 <x2) & & (Xr5 (F1 (1), F1 (2), F1 (3), F1 (4), F1 (5), 2) <= Xr5 (F2 (1), F2 (2 ), F2 (3), F2 (4), F2 (5), 2)))
b = x2;
end
if ((x1> x2) & & (Xr5 (F1 (1), F1 (2), F1 (3), F1 (4), F1 (5), 2) <= Xr5 (F2 (1), F2 (2 ), F2 (3), F2 (4), F2 (5), 2)))
a = x2;
end
if ((x1 <x2) & & (Xr5 (F1 (1), F1 (2), F1 (3), F1 (4), F1 (5), 2)> Xr5 (F2 (1), F2 (2) , F2 (3), F2 (4), F2 (5), 2)))
a = x1;
end
if ((x1> x2) & & (Xr5 (F1 (1), F1 (2), F1 (3), F1 (4), F1 (5), 2)> Xr5 (F2 (1), F2 (2) , F2 (3), F2 (4), F2 (5), 2)))
b = x1;
end
end
L = abs ((x2-x1) / 2);
Побудова еталонного і реального вихідного сигналу, пошук значення функціоналу.
function J = Xr5 (Kp, Kd, Ki, Kp2, Kd2, t)
% Коефіцієнти ДУ
a4 = 1;
a3 = (7.4 +5 * Kd * Kp2 +5 * Kp * Kd2 +10 * Kd * Kd2) / (1 +5 * Kd * Kd2);
a2 = (14 +5 * Kp * Kp2 +10 * Kp2 * Kd +10 * Kd2 * Kp +5 * Ki * Kd2) / (1 +5 * Kd * Kd2);
a1 = (10 * Kp * Kp2 +10 * Ki * Kd2 +5 * Kp2 * Ki) / (1 +5 * Kd * Kd2);
a0 = (10 * Ki * Kp2) / (1 +5 * Kd * Kd2);
b3 = 5 * Kd / (1 +5 * Kd * Kd2);
b2 = (5 * Kp +10 * Kd) / (1 +5 * Kd * Kd2);
b1 = (10 * Kp +5 * Ki) / (1 +5 * Kd * Kd2);
b0 = 10 * Ki / (1 +5 * Kd * Kd2);
% Крок
h = 0.01;
% Початкові умови
X (1,1:4) = [0 0 0 0];
% Матриця А (X '= AX + BY)
A = [0 1 0 0;
0 0 1 0;
0 0 0 1;
-A0-a1-a2-a3];
% Матриця B
B = [b3 b2 b1 b0] ';
Y = [0 0 0 1] ';
n = 1;
k = 1;
while (k <= 10)
F (n, 1:4) = (A * (X (n, 1:4)) '+ B. * Y)';
X (n +1,1:4) = (X (n, 1:4) '+ h/10 * (F (n, 1:4 ))')';% Метод Ейлера
n = n +1;
k = k +1;
end
X (2,1:4) = X (n, 1:4);
n = 2;
I = [1 0 0 0; 0 1 0 0, 0 0 1 0, 0 0 0 1];
while (n <= t / h)
F (n, 1:4) = (A * (X (n, 1:4)) '+ B. * Y)';
X (n +1,1:4 )=((( I-5 * h/12 * A) ^ -1) * (X (n, 1:4) '+ h/12 * (5 * B. * Y +8 * (F (n, 1:4 ))'-( F (n-1, 1:4 ))')))';
n = n +1;
end
i = 1;
while (i <= n)
Xr (i) = X (i, 1);
i = i +1;
end
% Знайдемо J0____________
i = 0;
while (i <= t / h)
Xe (i +1) = 1-exp ((log (0.05) / 2) * i * h);
i = i +1;
end
% Знайшли еталон
J = 0;
i = 1;
while (i <= t / h)
J = J + (Xr (i)-Xe (i)) ^ 2;
i = i +1;
end
% T = 0:0.01: t;
% Plot (t, Xr, t, Xe), grid on