6.4 Розв’язок нелінійної математичної моделі об’єкта числовим методом
Процес проектування автоматичних систем керування передбачає створення математичних моделей технологічних об’єктів. Такі моделі у більшості випадків мають вигляд диференціальних рівнянь. Для отримання відомостей про властивості об’єктів необхідно знати розв’язки їх моделей. Отримати аналітичний розв’язок таких моделей можна тільки в окремих випадках. У переважній більшості для розв’язку диференціальних рівнянь. Застосовують числові методи такі як метод Ейлера, метод Гюна, метод Рунге-Кутта та його модифікації. Всі перераховані методи дають можливість розв’язувати математичні моделі, які мають форму рівняння Коші
(6.12)
з початковою умовою
, (6.13)
де - деяка відома у загальному випадку, нелінійна функція двох аргументів.
Метод Ейлера має обмежене застосування із-за великої похибки, яка накопичується в процесі обчислень.
Метод Гюна має точність на порядок вищу, ніж метод Ейлера. Якщо метод Ейлера забезпечує точність розв’язку нелінійного рівняння порядку , то метод Гюна - .
Метод Рунге-Кутта. Найпоширенішим із сімейства методів Рунге-Кутта є метод четвертого порядку (р=4) або просто метод Рунге-Кутта, який породжує таку ітераційну процедуру:
,
,
,
,
.
Метод Рунге-Кутта дає похибку накопичення четвертого порядку - .
У випадку, коли математична модель подана у вигляді системи диференціальних рівнянь
,
,
(6.14)
з початковою умовою
,
де
- відомі функції двох аргументів.
Якщо рівняння (6.14) записати у векторній формі
(6.15)
з початковою умовою
,
(6.16)
то
отримаємо рівняння, яке аналогічне за
своєю структурою рівнянню (6.12). Отже,
для розв’язку рівняння (6.15) з початковою
умовою (6.16) модна застосувати процедуру
Рунге-Кутта, у якій скалярні величин
,
замінені векторними величинами
,
.
У випадку, що розглядається
,
а початкові умови будуть такими:
.
Розв'язок
задачі отримаємо за умови, що вхідними
величинами є зміна витрати
і командного тиску
,
тобто
,
,
де
,
;
,
- задані величини, які визначають величину
збурення. Програма розв'язку задачі,
яка написана у середовищі MatLab,
наведена на рис. 6.4, а
М-файл fun_RK
– на рис. 6.5.
%====================================================
%Метод Рунге-Кутта розв'язку Y'=F(t,y) задачі
%Коші на інтервалі [t0;tk] з початковою
%умовою y(0)=Y0
%====================================================
%Вхід: F(t,y)-вектор-функція, значення якої
% задати у М-файлі
% fun_RM_Sys_5
% to,tk-початкова і кінцева точки
% інтегрування
% Y0-вектор початкових умов
% n-число ітерацій
%Вихід: t-вектор абсцис
% y-матриця ординат
%----------------------------------------------------
%Вхідні дані
%----------------------------------------------------
%Параметри нелінійної моделі alfa, alfa20,alfa3
%a, V0 та коефіцієнти інтерполяційного aLag
%поліному записні у mat-файлі andrew2
%Значення величин - витрати q10, тиску P2,
%площі поперечного перерізу S1, S2, початкові
%умови H10,H20,значення командного тиску U0
%записні у mat-файлі andrew1
load andrew1
load andrew2
%----------------------------------------------------
%Значення змін вхідних величин А1, А2
A1=0.11;
A2=0.13;
%Вхідні дії
u1=(1+A1)*q10;
u2=(1+A2)*U0;
%----------------------------------------------------
%Обчислення коефіцієнта alfa2
N=length(aLag);
alfa2=0;
for i=1:N
alfa2=alfa2+aLag(i)*u2^(N-i);
end
%----------------------------------------------------
n=1000;
t0=0;
tk=1000;
y10=H0;
y20=T0;
%----------------------------------------------------
%Обчислюємо крок h та ініціалізуємо вектори t,y
%----------------------------------------------------
h=(tk-t0)/n;
%----------------------------------------------------
%Ітераційна процедура методу Рунге-Кутта
%----------------------------------------------------
t=t0:h:tk;
Y0=[y10 y20];
%У перший рядок матриці y занесені початкові
%умови
y(1,:)=Y0;
for k=1:n
f1=fun_RK(t(k),y(k,:),alfa20,alfa3,...
u1,u2,V,P2);
f2=fun_RK(t(k)+h/2,y(k,:)+h*f1/2,...
alfa20,alfa3,u1,u2,V,P2);
f3=fun_RK(t(k)+h/2,y(k,:)+h*f2/2,...
alfa20,alfa3,u1,u2,V,P2);
f4=fun_RK(t(k)+h,y(k,:)+h*f3,...
alfa20,alfa3,u1,u2,V,P2);
y(k+1,:)=y(k,:)+(f1+2*f2+2*f3+f4)*h/6;
end
%----------------------------------------------------
%Графік розв'язку задачі
%----------------------------------------------------
plot(t,y(:,1))
figure
plot(t,y(:,2))
save('ch','y')
save('RevoltingAction','A1','A2','n','tk')
Рисунок 6.4 – Програма розв'язку нелінійної математичної моделі методом Рунге-Кутта
function f=fun_RK(t,y,alfa20,alfa3, u1,u2,V,P2)
%====================================================
%Функція задачі Коші Y'=F(t,Y) (Метод Рунге-
%Кутта)
%====================================================
%Вхід:
% ro- густина рідини
% g - земне прискорення
%----------------------------------------------------
ro=1000;
g=9.807;
%----------------------------------------------------
F1=(u1-alfa3*sqrt(ro*g*(y(1)-y(2)))/(ro*V);
F2=(alfa3*sqrt(ro*g*(y(1)-y(2))-Y)+alfa20*...
sqrt(P2-ro*g*y(2)-Y)-...
alfa3*sqrt(ro*g*y(2))/(ro*V);
f=[F1 F2];
Рисунок 6.5 – М-файл до програми рис. 6.4
На
рис. 6.6 і 6.7 показані графіки змін рівнів
рідин у першій та другій ємностях після
змін вхідних величин на значення
і
,
із яких видно, що рівень рідини як у
першій, так і у другій ємностях збільшився
на певне значення і при
набувають
нового постійного значення.
Рисунок 6.6 – Зміна рівня рідини у першій ємності
Рисунок 6.7 – Зміна рівня рідини у другій ємності