Побудова математичної моделі та числове її дослідження

Предмет:

Математика

Тип роботи:

Курс лекцій

К-сть сторінок:

Мова:

Українська

Оцінка:

Створимо файл yos_kr41.m:

A=[1 1;

lb(1) lb(2)];

H=[h0;h0p];

C=linsolve(A,H);

C=double(C)

Одержимо константи:

C =1.0e-005 *

-0.61028435734888

0.61028435734888

Для знаходження імпульсної перехідної функції запишемо рівняння:

.Підставимо h(t):

. (3.13)

б) Знайдемо імпульсну перехідну функцію системи (3.8):

(3.14)

Початкові умови , змінені за рахунок дії одиничного імпульсного вхідного сигналу , знайдемо застосовуючи відомі залежності:

Для знаходження розв*язку рівняння (3.14), знайдемо корені його характеристичного рівняння. Для цього створимо скрипт-файл (yos_kr5.m):

kr3;

A3=1; A2=-a11-a22-a33;

A1=a11*a33+a22*a33+a11*a22;

A0=-a11*a22*a33;B2=0;

B1=a32*b2; B0=a12*a31*b2-a32*a11*b2;

Т0=0;Т0p=B1;Т0p2=B0-A2*B1;

lb=roots([A3 A2 A1 A0])

Результати виконаня:

lb = -0.03200000000000

-0.01775609460462

-0.01561419577680

Розв*язок системи має вигляд:

, (3.15)

де ,C3- константи, які знайдемо з початкових умов.

(3.16)

Для знаходження констант створимо файл yos_kr51.m:

A=[1 1 1;

lb(1) lb(2) lb(3);

lb(1)^2 lb(2)^2 lb(3)^2];

H=[Т0;Т0p;Т0p2];

C=linsolve(A,H);

C=double(C);

Результати:

C = 0.00012704856739

0.00088148763791

-0.00100853620530

Для знаходження імпульсної перехідної функції запишемо рівняння:

.Підставимо h(t):

(3.17)

3.3 Порівняння графіків перехідних функцій, отриманих за аналітичним розв*язком та числовим методом.

а) по рівню h

Створимо файл yos_kr6.m:

t=[0:3:500];

yos_kr3;

yos_kr4;yos_kr41;

h=C(1)*exp(lb(1)*t)+C(2)*exp(lb(2)*t);

sys=ss(a,b,c,d);

[y,tt,x]=impulse(sys,500),

plot(t,h,*.k*,tt,x(:,1));grid;

title(*impulsna perexidna f*);

ylabel(*h-h0,m*);

xlabel(*time,sec*);

pause;

r=C(1)/lb(1)*(exp(lb(1)*t)-1)+C(2)/lb(2)*(exp(lb(2)*t)-1);

sys=ss(a,b,c,d);

[y,tt,x]=step(sys,500),

plot(t,r,*.k*,tt,x(:,1));

grid;

title(*perehidna f*);

ylabel(*h-h0,m*);

xlabel(*time,sec*);

В результаті отримуємо:

Рисунок 3.1. Графіки порівняння імпульсної перехідної функції,одержаної:

«-»- за аналітичною залежністю; «·» - із застосуванням функції IMPULSE.

Рисунок 3.2. Графіки порівняння перехідної функції,одержаної: «-»- за аналітичною залежністю; «·» - із застосуванням функції STEP

б)по температурі Т:

Створимо файл yos_kr61.m:

t=[0:8:500];

yos_kr3;

yos_kr5;yos_kr51;

h=C(1)*exp(lb(1)*t)+C(2)*exp(lb(2)*t)+C(3)*exp(lb(3)*t);

sys=ss(a,b,c,d);

[y,tt,x]=impulse(sys,500),

plot(t,h,tt,x(:,3),*.k*);grid;

title(*impulsna perexidna f*);

ylabel(*T-T0,K*);

xlabel(*time,sec*);pause;

r=C(1)/lb(1)*(exp(lb(1)*t)-1)+C(2)/lb(2)*(exp(lb(2)*t)-1)+C(3)/lb(3)*(exp(lb(3)*t)-1);

sys=ss(a,b,c,d);

[y,tt,x]=step(sys,500),

plot(t,r,tt,x(:,3),*.k*);grid;

title(*perehidna f*);

ylabel(*T-T0,K*);

xlabel(*time,sec*);

В результаті отримуємо:

Рисунок.3.3. Графіки порівняння імпульсної перехідної функції,одержаної:

«-»- за аналітичною залежністю; «·» - із застосуванням функції IMPULSE.

Рисунок 3.4 Графіки порівняння перехідної функції,одержаної: «-»- за аналітичною залежністю; «·» - із застосуванням функції STEP.

3.4 Аналітичні вирази для визначення реакції системи на вхідний синусоїдальний сигнал ( ) за допомогою інтеграла згортки

Запишемо інтеграл згортки:

а) по рівню h.

Для порівняння графіків за допомогою аналітичних виразів створимо файл yos_kr7.m:

yos_kr3;

c=[1 0 0]; d=[0];

sys = ss(a, b, c, d);

w=0.2;

t=[0:1:500];

U=sin(w*t);

[y,t] = lsim(sys, U,t);

yos_kr4;yos_kr41;

s=C(1)*(w*exp(lb(1)*t)-(w*cos(w*t)+lb(1)*sin(w*t)))/(w^2+lb(1)^2)+C(2)*(w*exp(lb(2)*t)-(w*cos(w*t)+lb(2)*sin(w*t)))/(w^2+lb(2)^2);

figure(1)

plot(t,s,t,y,*.k*);grid;

ylabel(*h,m*);

xlabel(*time,sec*);

yos_kr3;

w=0.2;

t=[0:1:500];

yos_kr5;yos_kr51;

s=C(1)*(w*exp(lb(1)*t)-(w*cos(w*t)+lb(1)*sin(w*t)))/(w^2+lb(1)^2)+C(2)*(w*exp(lb(2)*t)-(w*cos(w*t)+lb(2)*sin(w*t)))/(w^2+lb(2)^2)+C(3)*(w*exp(lb(3)*t)-(w*cos(w*t)+lb(3)*sin(w*t)))/(w^2+lb(3)^2);

c=[0 0 1]; d=[0];

sys = ss(a, b, c, d);

y = lsim(sys, sin(w*t),t);

figure(2)

plot(t,s,t,y,*.k*);grid;

ylabel(*T,K*);

xlabel(*time,sec*);

Отримуємо графіки:

Рисунок 3.5. Графіки порівняння перехідних процесів по рівню h при дії вхідного сигналу u=sinωt, одержані: «-»- за аналітичною залежністю; «·» - із застосуванням функції LSIM.

Рисунок 3.6.Графіки порівняння перехідних процесів по температурі Т при дії вхідного сигналу u=sinωt, одержані: «-»- за аналітичною залежністю; «·» - із застосуванням функції LSIM.

4. Частотні методи аналізу системи

Амплітудочастотна характеристика (А(w)) - це залежність відношення амплітуди вихідного періодичного сигналу до вхідної амплітуди цього сигналу як функції частоти в усталеному режимі роботи.

Фазочастотна характеристика ((w)) - це залежність зміни фази вихідного сигналу по відношенню до фази вхідного сигналу від частоти в усталеному режимі роботи.

а) по рівню h

В усталеному режимі коливань складова прямує до нуля. Отже, при ,

Звідси визначаємо АЧХ і ФЧХ:

Запишимо рівняння (3.5) замінивши ,отримаємо:

(4.1)

З (4.1) виділимо дійсну та уявну частини:

Для порівняння графіків АЧФ і ФЧХ по рівню h створимо файл yos_kr8.m:

yos_kr3;

w=[0:0.01:0.5];

yos_kr4;yos_kr41;

p=j*w;

Hh=B0./(A2.*p.^2+A1.*p+A0);

Aw=sqrt(real(Hh).^2+imag(Hh).^2);

[Avuh,fivuh]=bode([B0],[A2 A1 A0],w);

plot(w,Aw,w,Avuh,*.k*);grid;

ylabel(*Аvuh/Аvh*);

xlabel(*w,rad/s*);

pause;

Fiw=atan(imag(Hh)./real(Hh));

Fiw=Fiw*180/pi;

plot(w,Fiw,w,fivuh,*.k*);grid;

ylabel(*fivuh/fivh*);

xlabel(*w,rad/s*);

Одержуємо такі графіки:

Рисунок 4.3 Порівняння графіків амплітудо-частотної характеристики системи : «-» - за аналітичною залежністю; «·» - із застосуванням функції MatLab BODE

Рисунок 4.4 Порівняння графіків фазо-частотної характеристики системи: «-» - за аналітичною залежністю; «·» - із застосуванням функції MatLab BODE.

б) по температурі Т:

В усталеному режимі коливань складова прямує до нуля. Отже, при ,

Звідси визначаємо АЧХ і ФЧХ:

Запишимо рівняння (3.5) замінивши ,отримаємо:

(4.2)

З (4.2) виділимо дійсну та уявну частини:

Для порівняння графіків АЧФ і ФЧХ по рівню T створимо файл yos_kr81.m:

yos_kr3;

w=[0:0.007:0.5];

yos_kr5;yos_kr51;

p=j*w;

Hh=B1.*p+B0./(A3.*p.^3+A2.*p.^2+A1.*p+A0);

Aw=sqrt(real(Hh).^2+imag(Hh).^2);

[Avuh,fivuh]=bode([B1 B0],[A3 A2 A1 A0],w);

plot(w,Aw,w,Avuh,*.k*);grid;

ylabel(*Аvuh/Аvh*);

xlabel(*w,rad/s*);

pause;

Fiw=atan(imag(Hh)./real(Hh));

Fiw=Fiw*180/pi;

plot(w,Fiw,w,fivuh,*.k*);grid;ylabel(*fivuh/fivh*);xlabel(*w,rad/s*);

Одержуємо такі графіки:

Рисунок. 4.5 Порівняння графіків амплітудо-частотної характеристики системи: «-» - за аналітичною залежністю; «·» - із застосуванням функції MatLab BODE

Рисунок 4.6 Порівняння графіків фазо-частотної характеристики системи :

«-» - за аналітичною залежністю; «·» - із застосуванням функції MatLab BODE.

5. Дослідження моделі в середовищі SimuLink

Дослідимо отриману математичну модель в середовищі SimuLink. Побудуємо, використовуючи блоки SimuLink, структурну модель двоємнісного обєкту.

Для запуску моделі в SimuLink попередньо потрібно ініціалізувати початкові значення параметрів стану та констант з допомогою програми yos_kr9.m

P1=2000;P2=9500;

T1=293; T2=350;

L1=100; r1=0.05;

kv=0.005;

l2=0.4; l3=0.7;

d=1; nu=0.00001;

dz=0.9; ro=1e3; g=9.81;

S=pi*d^2/4;

kl=pi*r1^4/8/L1/nu;

A=L1*kl/pi/r1^2;

h0=0.617917411421; Q10=0.004908738521; T0=317.530396154693;

Рисунок 5.1 Модель системи у вікні SIMULINK

В даній моделі встановимо для блоків наступні параметри:

MATLAB Fcn: kv*l2*sqrt((P2-ro*g*u)/ro)-kv*l3*sqrt(g*u);

MATLAB Fcn1: kl*u/ro;

MATLAB Fcn2: (u(1)*(T1-u(3))+kv*l2*sqrt((P2-ro*g*u(2))/ro)*(T2-u(3)))/u(2);

Gain: 1/S;

Gain1: 1/A;

Gain1: 1/S;

Встановимо початкові умови для інтеграторів:

Integrator: h0;

Integrator1: Q10;

Integrator2: T0;

Параметри блоку STEP:

Step time: 0;

Initial value: 2000;

Final value: 2400;

Параметри симуляції:

Вкладка Solver:

Start time: 0;

Stop time: 500;

Save to workspace:

Time: tout;

States: xout;

Порівняння одержаних перехідних процесів в середовищах Matlab та Simulink.

Для накладання графіків перехідних процесів, отриманих в середовищі SimuLink і за допомогою розв*язування системи рівнянь функцією ode45 потрібно запустити на виконання файл yos_kr10.m :

x0=[0.617917411421 0.004908738521 317.530396154693];

T=[0 500];

tol=odeset(*RelTol*,3e-14);

[t,y]=ode45(*yos_kr2*,T,x0,tol);

plot(t,y(:,1),tout,xout(:,1),*.k*);grid;

ylabel(*h,m*);

xlabel(*t,sec*);pause;

plot(t,y(:,2),tout,xout(:,2),*.k*);grid;

ylabel(*Q1,mkub/sec*);

xlabel(*t,sec*);pause;

plot(t,y(:,3),tout,xout(:,3),*.k*);grid;

ylabel(*T,K*);

xlabel(*t,sec*);

Результатом виконання програми є наступні графіки:

Рисунок 5.1 Графіки порівняння перехідних процесів для h, отримані:

«-» - шляхом розв*язку системи диференціальних рівнянь за допомогою функції ode45; «·» - шляхом побудови схеми моделі у вікні SIMULINK.

Рисунок 5.2 Графіки порівняння перехідних процесів для Q1, отримані: «-» - шляхом розв*язку системи диференціальних рівнянь за допомогою функції ode45; «·» - шляхом побудови схеми моделі у вікні SIMULINK

Рисунок 5.3 Графіки порівняння перехідних процесів для T, отримані: «-» - шляхом розв*язку системи диференціальних рівнянь за допомогою функції ode45; «·» - шляхом побудови схеми моделі у вікні SIMULINK.