第7章 微分方程模型仿真

1、第第7章章 微分方程模型仿真微分方程模型仿真常微分方程的數(shù)值求解常微分方程的數(shù)值求解微分方程模型的建立及仿真微分方程模型的建立及仿真7.1 微分方程的求解微分方程的求解 在現(xiàn)在數(shù)學(xué)研究和工程實踐中，很多數(shù)學(xué)模型都是在現(xiàn)在數(shù)學(xué)研究和工程實踐中，很多數(shù)學(xué)模型都是用微分方程確定的，很多基本方程本身就是一個微分方用微分方程確定的，很多基本方程本身就是一個微分方程，因此求微分方程非常重要，但是大部分的微分方程程，因此求微分方程非常重要，但是大部分的微分方程目前難以求得其解析解，因此人們只有利用計算機強大目前難以求得其解析解，因此人們只有利用計算機強大的計算功能來求其數(shù)值解。的計算功能來求其數(shù)值解。MAT

2、LAB主要使用龍格主要使用龍格-庫庫塔法求解微分方程。塔法求解微分方程。 在控制系統(tǒng)仿真中，常用的求微分方程數(shù)值解的函在控制系統(tǒng)仿真中，常用的求微分方程數(shù)值解的函數(shù)是數(shù)是ode23和和ode45。1. ode23 在在MATLAB中，函數(shù)中，函數(shù)ode23采用采用2-3階龍格階龍格-庫塔法庫塔法求解微分方程。求解微分方程。t,y=ode23(odefun,tspan,y0)t,y=ode23(odefun,tspan,y0,options)odefun：定義微分方程的形式：定義微分方程的形式y(tǒng)=f(t,y)tspan=t0,tfinal：表示微分方程的積分限從：表示微分方程的積分限從t0（始值

3、）（始值）到到tfinal（終值），該積分限也可以是一些離散的點。（終值），該積分限也可以是一些離散的點。y0:初始狀態(tài)列向量初始狀態(tài)列向量options：積分參數(shù)，包括積分參數(shù)，包括RelTol（相對誤差）和（相對誤差）和AbsTol（絕對誤差），可省略。（絕對誤差），可省略。例例：使用：使用ode23函數(shù)求解常微分方程函數(shù)求解常微分方程y=-y+x2+4x+1，x=1 4， x=1時，時，y=1。解解：首先創(chuàng)建函數(shù)：首先創(chuàng)建函數(shù)fun1.mfunction f=fun1(x,y)f=-y+x2+4*x+1;在命令窗口中輸入在命令窗口中輸入 x,y=ode23(fun1,1,4,1); dy

4、=-y+x.2+4*x+1; plot(x,y,x,dy); legend(y,dy)11.522.533.540510152025ydy2. ode45 在在MATLAB中，函數(shù)中，函數(shù)ode45采用普通采用普通4-5階龍格階龍格-庫庫塔法求解微分方程。其使用方法與塔法求解微分方程。其使用方法與ode23函數(shù)的使用方函數(shù)的使用方法基本相同。法基本相同。 ode45函數(shù)是大部分場合的首選算法，函數(shù)是大部分場合的首選算法，ode23函數(shù)主函數(shù)主要適用于精度較低的場合。要適用于精度較低的場合。例例：解經(jīng)典非線性方程，范得波（：解經(jīng)典非線性方程，范得波（Van der Pol）微分方）微分方程（程（

5、=2）。）。0)1 (222xdtdxxdtxd當(dāng)當(dāng)t=0時，時，x=1，dx/dt=0。0)1 (222xdtdxxdtxd解解：（：（1）將高階微分方程式等價變換成一階微分方將高階微分方程式等價變換成一階微分方程組。程組。令令y1=x且且y2=dx/dt dy1/dt=y2 dy2/dt=w(1-y12)y2-y1（2）編寫編寫M文件表示該微分方程，該文件給定時間及文件表示該微分方程，該文件給定時間及y1和和y2的值，返回上述的導(dǎo)數(shù)值，并將的值，返回上述的導(dǎo)數(shù)值，并將y（y1和和y2）與）與導(dǎo)數(shù)值以導(dǎo)數(shù)值以列向量列向量的形式給出。的形式給出。function fun2=vdpol(t,y)

6、fun2=y(2) 2*(1-y(1)2)*y(2)-y(1)%輸出結(jié)果必須是列向量，輸出結(jié)果必須是列向量，w=2（3）計算結(jié)果如下：計算結(jié)果如下： t,y=ode45(vdpol,0 30,1;0); y1=y(:,1); y2=y(:,2); plot(t,y1,:b,t,y2,-r) legend(位移位移,速度速度)3. 定積分的數(shù)值解法定積分的數(shù)值解法 MATLAB軟件使用軟件使用quad函數(shù)進行定積分的數(shù)值解函數(shù)進行定積分的數(shù)值解法。使用格式為：法。使用格式為：q = quad(fun,a,b)fun：被積分函數(shù)：被積分函數(shù)a、b：積分上下限：積分上下限例例：計算下列定積分：計算下

7、列定積分xxxd521203function y = myfun(x) y = 1./(x.3-2*x-5); Q = quad(myfun,0,2)Q = -0.4605 F = (x)1./(x.3-2*x-5); Q = quad(F,0,2)Q = -0.46057.2 微分方程模型微分方程模型7.2.1 方法描述方法描述 微分方程模型是數(shù)學(xué)模型的一種主要形式。當(dāng)采用微分方程模型是數(shù)學(xué)模型的一種主要形式。當(dāng)采用一階微分方程的數(shù)值積分法進行數(shù)值計算時，應(yīng)把高階一階微分方程的數(shù)值積分法進行數(shù)值計算時，應(yīng)把高階微分方程變換成微分方程變換成n個一階微分方程形式。對于微分方程個一階微分方程形式。

8、對于微分方程而言，除了少數(shù)可以得到解析解外，大多數(shù)只能采用數(shù)而言，除了少數(shù)可以得到解析解外，大多數(shù)只能采用數(shù)值解法。值解法。 在在MATLAB中，使用中，使用ode函數(shù)求解微分方程模型。函數(shù)求解微分方程模型。7.2.2 簡單電路模型仿真簡單電路模型仿真例例：在：在RC電路中，電阻電路中，電阻R=5，理想電壓源為，理想電壓源為Vi=20V，電容電容C=70F。分析電容元件的電壓特性。分析電容元件的電壓特性。（1）分析：電容電壓和電流的關(guān)系分析：電容電壓和電流的關(guān)系dttdVCtIdttICtVCC)()()(1)(根據(jù)基爾霍夫定律，可得出微分方程根據(jù)基爾霍夫定律，可得出微分方程iCCiCVtVd

9、ttdVRCVtVtRI)()()()(使用使用ode函數(shù)時，對微分方程進行如下假設(shè)函數(shù)時，對微分方程進行如下假設(shè))(tVyCRCyVdtdyi（2）建立導(dǎo)數(shù)函數(shù)建立導(dǎo)數(shù)函數(shù)function dy=cap(t,y)Vi=20;R=5;C=70e-6;dy=(Vi-y)/(R*C);（3）使用使用ode函數(shù)進行仿真，仿真時間函數(shù)進行仿真，仿真時間00.006s，Vc初始值為初始值為0V。 t,y=ode45(cap,0,0.006,0); plot(t,y) axis(0 0.006 0 25) title(Vc-Time) xlabel(Time/sec) ylabel(Vc/V) 當(dāng)電壓源為

10、直流電壓源時，加載在電容上的電壓當(dāng)電壓源為直流電壓源時，加載在電容上的電壓隨時間呈拋物線增大，穩(wěn)態(tài)值為電源電壓。電容電壓隨時間呈拋物線增大，穩(wěn)態(tài)值為電源電壓。電容電壓在在t=0時取得最小值，最小值為時取得最小值，最小值為0；電容電壓在；電容電壓在t=0.0023s時達(dá)到最大值，為時達(dá)到最大值，為20V。例例：用簡單的：用簡單的LC諧振電路組成濾波器電路，其電路方諧振電路組成濾波器電路，其電路方程是二階微分方程。觀察該程是二階微分方程。觀察該RLC電路中，電容電壓電路中，電容電壓VC(t)和線路電流和線路電流I(t)的時域變化情況。假設(shè)電源為直流的時域變化情況。假設(shè)電源為直流電壓源電壓源Vi=2

11、0V，電阻，電阻R=5，電容，電容C=70F，電感，電感L=70mH。（1）分析：根據(jù)電路分析，可以得出微分方程分析：根據(jù)電路分析，可以得出微分方程2222LCiCCLCiCCCiRI(t) V (t) V (t)V (t)dV (t)d V (t)dI(t)I(t)CV (t)LLCdtdtdtdI(t)LRI(t) V (t)V (t)dtd V (t)dV (t)LCRCV (t)V (t)dtdt在利用在利用ode函數(shù)時，對微分方程作出如下假設(shè)：函數(shù)時，對微分方程作出如下假設(shè)：)2(*) 1 ()(1)2()2(1) 1 ()()2()() 1 (yRytVLdtdyyCdtdytIy

12、tVyiC（2）創(chuàng)建狀態(tài)導(dǎo)數(shù)函數(shù)創(chuàng)建狀態(tài)導(dǎo)數(shù)函數(shù)function dy=RLC(t,y)Vi=20;R=5;C=70e-06;L=70e-03;dy=zeros(2,1);dy(1)=y(2)/C;dy(2)=(Vi-y(1)-R*y(2)/L;（3）使用使用ode函數(shù)進行仿真，仿真時間函數(shù)進行仿真，仿真時間00.12s，Vc初始值為初始值為0V，I初始值為初始值為0A。 t,y=ode45(RLC,0 0.12,0;0);figure(1) subplot(2,1,1); plot(t,y(:,1); title(Vc-Time) xlabel(Time/sec) ylabel(Vc/V)

13、subplot(2,1,2); plot(t,y(:,2); title(I-Time) xlabel(Time/sec) ylabel(I/A) figure(2) plot(y(:,2),y(:,1) title(Vc-I) xlabel(I/A) ylabel(Vc/V) 在電容電壓與時間曲線中，可以得出如下結(jié)論：在電容電壓與時間曲線中，可以得出如下結(jié)論：當(dāng)電壓源為直流電壓源時，加載在電容上的電壓隨時當(dāng)電壓源為直流電壓源時，加載在電容上的電壓隨時間呈震蕩衰減；在電路電流與時間曲線中，可以得出間呈震蕩衰減；在電路電流與時間曲線中，可以得出結(jié)論，當(dāng)電壓源為直流電壓源時，流經(jīng)電路的電流隨結(jié)論，

14、當(dāng)電壓源為直流電壓源時，流經(jīng)電路的電流隨時間呈震蕩衰減；在電容電壓與電路電流的曲線圖中時間呈震蕩衰減；在電容電壓與電路電流的曲線圖中可以得出如下結(jié)論：當(dāng)處于穩(wěn)態(tài)時，電容電壓為可以得出如下結(jié)論：當(dāng)處于穩(wěn)態(tài)時，電容電壓為20V，電路電流為電路電流為0A。7.2.3 直流直流/直流變換器系統(tǒng)的仿真直流變換器系統(tǒng)的仿真 直流直流/直流變換器（又稱為直流變換器（又稱為DC/DC變換器）變換器）廣泛應(yīng)用于開關(guān)穩(wěn)壓電源以及直流電動機的控廣泛應(yīng)用于開關(guān)穩(wěn)壓電源以及直流電動機的控制。它主要利用開關(guān)器件進行通斷控制，將直制。它主要利用開關(guān)器件進行通斷控制，將直流電源斷續(xù)地加在負(fù)載上，通過控制開關(guān)開通流電源斷續(xù)地加

15、在負(fù)載上，通過控制開關(guān)開通與關(guān)斷時間的變化來改變負(fù)載電壓平均值，實與關(guān)斷時間的變化來改變負(fù)載電壓平均值，實現(xiàn)輸出電壓的可調(diào)?，F(xiàn)輸出電壓的可調(diào)。常用的主要有以下幾種類型：常用的主要有以下幾種類型：（1）降壓型降壓型（Buck）變換器）變換器（2）升壓型升壓型（Boost）變換器）變換器（3）升升-降壓型降壓型（ Buck- Boost）變換器）變換器通過控制變換器的輸出可得到用戶所期望通過控制變換器的輸出可得到用戶所期望的電壓值，實現(xiàn)直流電壓的有效調(diào)節(jié)，比如用的電壓值，實現(xiàn)直流電壓的有效調(diào)節(jié)，比如用直流直流/直流變換器控制電機，可有效實現(xiàn)直流電直流變換器控制電機，可有效實現(xiàn)直流電機的調(diào)壓調(diào)速。機

16、的調(diào)壓調(diào)速。閉環(huán)閉環(huán)DC/DC變換器的結(jié)構(gòu)框圖如下：變換器的結(jié)構(gòu)框圖如下：基本原理：基本原理：給定弱電信號給定弱電信號Uo*,閉環(huán)調(diào)節(jié)。就可對應(yīng)一閉環(huán)調(diào)節(jié)。就可對應(yīng)一個固定的占空比個固定的占空比d值，也對應(yīng)一個輸出電壓值，也對應(yīng)一個輸出電壓Uo值，首值，首先將給有初始值的輸出電壓先將給有初始值的輸出電壓Uo采樣，得到反饋信號采樣，得到反饋信號與給定弱電信號進行比較得到偏差，經(jīng)過與給定弱電信號進行比較得到偏差，經(jīng)過PI調(diào)節(jié)器調(diào)節(jié)器產(chǎn)生占空比值。調(diào)用產(chǎn)生占空比值。調(diào)用DC/DC變換器微分?jǐn)?shù)學(xué)模型，變換器微分?jǐn)?shù)學(xué)模型，并且使用并且使用龍格龍格-庫塔算法庫塔算法，可得到輸出電壓，可得到輸出電壓Uo值，

17、如值，如此反復(fù)直到占空比此反復(fù)直到占空比d與輸出電壓都達(dá)到穩(wěn)定值。與輸出電壓都達(dá)到穩(wěn)定值。(1)當(dāng)當(dāng)S=on：時時間長度為間長度為d*T:狀態(tài)方程以兩個狀態(tài)元件狀態(tài)方程以兩個狀態(tài)元件 i和和uo作狀態(tài)變量作狀態(tài)變量000idiLdtdiCdtRuuuu000ididtLLdidtCRCuuuuCS+-DR+-u0uiLi00110110iiiLLCRCuuu 1. .建立建立Buck型型DC/DC變換器的數(shù)學(xué)模型變換器的數(shù)學(xué)模型0000diLdtdiCdtRuuu 000didtLidtCRCuduu C CS S+ +- -D DR R+ +- -u0uiL Li iT TdTdT(2)當(dāng)當(dāng)

18、S=off：時間長度為時間長度為: （1-d)*T:00100110iiiLCRCuuu （3）狀態(tài)空間平均）狀態(tài)空間平均(SSA)： 00110110iiiLLCRCuuu UBXAXUBXAX221121AddAA 21BddBB BUAXX1dd S=OFFS=ON1dd d為為PWM占空比，得占空比，得DC/DC的統(tǒng)一的統(tǒng)一SSA模型為：模型為：uuuiLdiRCCLi0111000有上述有上述SSA模型得：模型得：RCCiLdLiuuuui000 當(dāng)當(dāng)d=1時，時，A=A1，B=B1，此時為開關(guān)開通的狀態(tài)模型；，此時為開關(guān)開通的狀態(tài)模型； 當(dāng)當(dāng)0d1 ?d1000) R=10; %在采樣過程中突減負(fù)載在采樣過程中突減負(fù)載 end T(n)=(n-1)*Tc; %對系統(tǒng)進行計時對系統(tǒng)進行計時 vot(n)=vo0;%