第7章 微分方程模型仿真

1、第第7章章 微分方程模型仿真微分方程模型仿真常微分方程的數(shù)值求解常微分方程的數(shù)值求解微分方程模型的建立及仿真微分方程模型的建立及仿真7.1 微分方程的求解微分方程的求解 在現(xiàn)在數(shù)學研究和工程實踐中，很多數(shù)學模型都是在現(xiàn)在數(shù)學研究和工程實踐中，很多數(shù)學模型都是用微分方程確定的，很多基本方程本身就是一個微分方用微分方程確定的，很多基本方程本身就是一個微分方程，因此求微分方程非常重要，但是大部分的微分方程程，因此求微分方程非常重要，但是大部分的微分方程目前難以求得其解析解，因此人們只有利用計算機強大目前難以求得其解析解，因此人們只有利用計算機強大的計算功能來求其數(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ù)主要適用于精度較低的場合。要適用于精度較低的場合。例例：解經典非線性方程，范得波（：解經典非線性方程，范得波（Van der Pol）微分方）微分方程（程（

5、=2）。）。0)1 (222xdtdxxdtxd當當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的值，返回上述的導數(shù)值，并將的值，返回上述的導數(shù)值，并將y（y1和和y2）與）與導數(shù)值以導數(shù)值以列向量列向量的形式給出。的形式給出。function fun2=vdpol(t,y)

6、fun2=y(2) 2*(1-y(1)2)*y(2)-y(1)%輸出結果必須是列向量，輸出結果必須是列向量，w=2（3）計算結果如下：計算結果如下： 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ù)學模型的一種主要形式。當采用微分方程模型是數(shù)學模型的一種主要形式。當采用一階微分方程的數(shù)值積分法進行數(shù)值計算時，應把高階一階微分方程的數(shù)值積分法進行數(shù)值計算時，應把高階微分方程變換成微分方程變換成n個一階微分方程形式。對于微分方程個一階微分方程形式。

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

9、ttdVRCVtVtRI)()()()(使用使用ode函數(shù)時，對微分方程進行如下假設函數(shù)時，對微分方程進行如下假設)(tVyCRCyVdtdyi（2）建立導數(shù)函數(shù)建立導數(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) 當電壓源為

10、直流電壓源時，加載在電容上的電壓當電壓源為直流電壓源時，加載在電容上的電壓隨時間呈拋物線增大，穩(wěn)態(tài)值為電源電壓。電容電壓隨時間呈拋物線增大，穩(wěn)態(tài)值為電源電壓。電容電壓在在t=0時取得最小值，最小值為時取得最小值，最小值為0；電容電壓在；電容電壓在t=0.0023s時達到最大值，為時達到最大值，為20V。例例：用簡單的：用簡單的LC諧振電路組成濾波器電路，其電路方諧振電路組成濾波器電路，其電路方程是二階微分方程。觀察該程是二階微分方程。觀察該RLC電路中，電容電壓電路中，電容電壓VC(t)和線路電流和線路電流I(t)的時域變化情況。假設電源為直流的時域變化情況。假設電源為直流電壓源電壓源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ù)時，對微分方程作出如下假設：)2(*) 1 ()(1)2()2(1) 1 ()()2()() 1 (yRytVLdtdyyCdtdytIy

12、tVyiC（2）創(chuàng)建狀態(tài)導數(shù)函數(shù)創(chuàng)建狀態(tài)導數(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) 在電容電壓與時間曲線中，可以得出如下結論：在電容電壓與時間曲線中，可以得出如下結論：當電壓源為直流電壓源時，加載在電容上的電壓隨時當電壓源為直流電壓源時，加載在電容上的電壓隨時間呈震蕩衰減；在電路電流與時間曲線中，可以得出間呈震蕩衰減；在電路電流與時間曲線中，可以得出結論，當電壓源為直流電壓源時，流經電路的電流隨結論，

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

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

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

17、如值，如此反復直到占空比此反復直到占空比d與輸出電壓都達到穩(wěn)定值。與輸出電壓都達到穩(wěn)定值。(1)當當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ù)學模型變換器的數(shù)學模型0000diLdtdiCdtRuuu 000didtLidtCRCuduu C CS S+ +- -D DR R+ +- -u0uiL Li iT TdTdT(2)當當

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=1時，時，A=A1，B=B1，此時為開關開通的狀態(tài)模型；，此時為開關開通的狀態(tài)模型； 當當0d1 ?d1000) R=10; %在采樣過程中突減負載在采樣過程中突減負載 end T(n)=(n-1)*Tc; %對系統(tǒng)進行計時對系統(tǒng)進行計時 vot(n)=vo0;%