常微分方程的求解及系統(tǒng)數(shù)學(xué)模型的轉(zhuǎn)換實(shí)驗(yàn)報(bào)告

《計(jì)算機(jī)仿真》實(shí)驗(yàn)報(bào)告

實(shí)驗(yàn)一常微分方程的求解及系統(tǒng)數(shù)學(xué)模型的轉(zhuǎn)換一．實(shí)驗(yàn)?zāi)康耐ㄟ^(guò)實(shí)驗(yàn)熟悉計(jì)算機(jī)仿真中常用到的Matlab指令的使用方法，掌握常微分方程求解指令和模型表示及轉(zhuǎn)換指令，為進(jìn)一步從事有關(guān)仿真設(shè)計(jì)和研究工作打下基礎(chǔ)。二.實(shí)驗(yàn)設(shè)備個(gè)人計(jì)算機(jī)，Matlab軟件。三.實(shí)驗(yàn)準(zhǔn)備預(yù)習(xí)本實(shí)驗(yàn)有關(guān)內(nèi)容（如教材第2、3、5章中的相應(yīng)指令說(shuō)明和例題），編寫(xiě)本次仿真練習(xí)題的相應(yīng)程序。四.實(shí)驗(yàn)內(nèi)容1.Matlab中常微分方程求解指令的使用題目一：請(qǐng)用MATLAB的ODE45算法分別求解下列二個(gè)方程。要求:1.編寫(xiě)出Matlab仿真程序；2.畫(huà)出方程解的圖形并對(duì)圖形進(jìn)行簡(jiǎn)要分析；3.分析下列二個(gè)方程的關(guān)系。2． 解：程序：新建f1：functionf1=f1(t,x)f1=-x^2;end新建f2：functionf2=f2(t,x)f2=x^2;end運(yùn)行：[t,x]=ode45(@f1,[0,40],[1])figure(1);plot(t,x);[t,x]=ode45(@f2,[0,40],[-1])figure(2);plot(t,x)運(yùn)行截圖：結(jié)論：方程一和方程二的圖形關(guān)于橫軸對(duì)稱題目二：下面方程組用在人口動(dòng)力學(xué)中，可以表達(dá)為單一化的捕食者-被捕食者模式（例如，狐貍和兔子）。其中表示被捕食者，表示捕食者。如果被捕食者有無(wú)限的食物，并且不會(huì)出現(xiàn)捕食者。于是有，則這個(gè)式子是以指數(shù)形式增長(zhǎng)的。大量的被捕食者將會(huì)使捕食者的數(shù)量增長(zhǎng)；同樣，越來(lái)越少的捕食者會(huì)使被捕食者的數(shù)量增長(zhǎng)。而且，人口數(shù)量也會(huì)增長(zhǎng)。請(qǐng)分別調(diào)用ODE45、ODE23算法求解下面方程組。要求編寫(xiě)出Matlab仿真程序、畫(huà)出方程組解的圖形并對(duì)圖形進(jìn)行分析和比較。解：OED45：程序：新建f3：functionf3=f3(t,x)f3=[x(1)-0.1*x(1)*x(2)+0.01*t;-x(2)+0.02*x(1)*x(2)+0.04*t];end運(yùn)行：[t,x]=ode45('f3',[0,20],[30,20]);plot(t,x);title('ode45×÷í?');xlabel('oì???a2?ê3??￡?à????a±?2?ê3??');grid運(yùn)行截圖OED23：程序;新建f4：functionf4=f4(t,x)f4=[x(1)-0.1*x(1)*x(2)+0.01*t;-x(2)+0.02*x(1)*x(2)+0.04*t];end運(yùn)行：[t,x]=ode23('f4',[0,20],[30,20]);plot(t,x);title('ode23×÷í?');xlabel('oì???a2?ê3??￡?à????a±?2?ê3??');grid運(yùn)行截圖：2.Matlab中模型表示及模型轉(zhuǎn)換指令的使用題目三：若給定系統(tǒng)的的傳遞函數(shù)為請(qǐng)用MATLAB編程求解其系統(tǒng)的極零點(diǎn)模型。解：程序：num=[612610];den=[12311];[z,p,k]=tf2zp(num,den);Gzpk=zpk(z,p,k)[zeropolesk]=zpkdata(Gzpk,'v')運(yùn)行結(jié)果：Gzpk=6(s+1.929)(s^2+0.07058s+0.8638)-----------------------------------------------(s^2+0.08663s+0.413)(s^2+1.913s+2.421)Continuous-timezero/pole/gainmodel.zero=-1.9294+0.0000i-0.0353+0.9287i-0.0353-0.9287ipoles=-0.9567+1.2272i-0.9567-1.2272i-0.0433+0.6412i-0.0433-0.6412ik=6題目四：例題2.9，求其伴隨標(biāo)準(zhǔn)型。解：程序?yàn)椋簄um=[41716];den=[171612];sys=tf(num,den)canon(sys,'companion')運(yùn)行結(jié)果為：timusisys=4s^2+17s+16-----------------------s^3+7s^2+16s+12Continuous-timetransferfunction.ans=A=x1x2x3x100-12x210-16x301-7B=u1x11x20x30C=x1x2x3y14-1129D=u1y10Continuous-timestate-spacemodel.題目五：習(xí)題5.8，采樣周期改變?yōu)?.05s，求系統(tǒng)的數(shù)學(xué)仿真模型。解：程序：dnum=[10];dden=[1-0.30.02];sys=tf(dnum,dden,0.1)D=d2d(sys,0.05)運(yùn)行結(jié)果:timuwusys=z------------------z^2-0.3z+0.02Sampletime:0.1secondsDiscrete-timetransferfunction.D=0.6222z-0.09725-----------------------z^2-0.7634z+0.1414Sampletime:0.05secondsDiscrete-timetransferfunction.五．總結(jié)與體會(huì)通過(guò)本次實(shí)驗(yàn)，學(xué)會(huì)了使用matlab新建和運(yùn)行函數(shù)，和使用ODE23和ODE45函數(shù)，以及如何利用matlab畫(huà)圖，學(xué)會(huì)了matlab中模型表示及模型轉(zhuǎn)換的使用；學(xué)會(huì)了如何查看和使用圖上的信息；學(xué)會(huì)了求解系統(tǒng)的極零點(diǎn)模型。學(xué)會(huì)了函數(shù)伴隨標(biāo)準(zhǔn)型的產(chǎn)生以及系統(tǒng)的數(shù)學(xué)仿真模型，為以后的學(xué)習(xí)和實(shí)驗(yàn)打下了堅(jiān)實(shí)的基礎(chǔ)。

實(shí)驗(yàn)二Matlab優(yōu)化工具箱的使用一．實(shí)驗(yàn)?zāi)康耐ㄟ^(guò)上機(jī)操作熟悉Matlab優(yōu)化工具箱的主要功能及其使用方法，掌握優(yōu)化工具箱中常用函數(shù)的功能和語(yǔ)法，并利用其進(jìn)行極值運(yùn)算、求解線性和非線性問(wèn)題等，為進(jìn)一步的仿真設(shè)計(jì)和研究打下基礎(chǔ)。二.實(shí)驗(yàn)設(shè)備個(gè)人計(jì)算機(jī)，Matlab軟件。三.實(shí)驗(yàn)準(zhǔn)備預(yù)習(xí)本實(shí)驗(yàn)有關(guān)內(nèi)容（如教材第6章中的相應(yīng)指令說(shuō)明和例題），編寫(xiě)本次仿真練習(xí)題的相應(yīng)程序。四.實(shí)驗(yàn)內(nèi)容1.應(yīng)用Matlab優(yōu)化工具箱求解優(yōu)化問(wèn)題例題6.6原函數(shù)：新建myfun211functionF=myfun(x)F=[2*x(1)-x(2)-exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];調(diào)用優(yōu)化函數(shù)：x0=[-5;-5];options=optimset('Display','iter');[x,fval]=fsolve(@myfun211,x0,options)結(jié)果為：NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius0347071.22.29e+0411612003.415.75e+031293147.0211.47e+031312854.45213881415239.5271107151867.0412130.8162116.704219.0517242.4278812.2618270.0326580.7595110.2062.59307.03149e-060.1119270.002942.510333.29525e-130.001691326.36e-072.5Equationsolved.fsolvecompletedbecausethevectoroffunctionvaluesisnearzeroasmeasuredbythedefaultvalueofthefunctiontolerance,andtheproblemappearsregularasmeasuredbythegradient.<stoppingcriteriadetails>x=0.56710.5671fval=1.0e-06*-0.4059-0.4059更改參數(shù)后：新建myfun2111:functionF=myfun(x)F=[3*x(1)-1*x(2)-exp(-x(1));-2*x(1)+4*x(2)-exp(-x(2))];調(diào)用優(yōu)化函數(shù)：x0=[-5;-5];options=optimset('Display','iter');[x,fval]=fsolve(@myfun2111,x0,options)運(yùn)行結(jié)果為：yunxing2111NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius0350189.52.4e+041161337016.17e+031293740.911.67e+0313121108.1414751415343.11911451518104.483147162126.215115.417242.6223414.0318270.009749870.5367630.2152.59301.25882e-070.03695340.0007372.510332.06617e-170.0001337829.08e-092.5Equationsolved.fsolvecompletedbecausethevectoroffunctionvaluesisnearzeroasmeasuredbythedefaultvalueofthefunctiontolerance,andtheproblemappearsregularasmeasuredbythegradient.<stoppingcriteriadetails>x=0.35170.3517fval=1.0e-08*-0.3800-0.2495結(jié)論：更改參數(shù)前后，運(yùn)行過(guò)程中使用的方法是一樣的，運(yùn)行結(jié)果結(jié)構(gòu)也類(lèi)似。例題6.8:原函數(shù)：新建fun212:functionf=fun212(x)f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);運(yùn)行函數(shù)：x0=[-11];options=[];[x,options]=fminunc('fun212',x0,options)結(jié)果為：>>yunxing212Localminimumfound.Optimizationcompletedbecausethesizeofthegradientislessthanthedefaultvalueoftheoptimalitytolerance.<stoppingcriteriadetails>x=0.5000-1.0000options=3.6609e-16更改后：新建fun2121：functionf=fun212(x)f=exp(x(1))*(5*x(1)^2+6*x(2)^2+7*x(1)*x(2)+8*x(2)+1);運(yùn)行：x0=[-11];options=[];[x,options]=fminunc('fun2121',x0,options)結(jié)果為：>>yunxing2121Localminimumfound.Optimizationcompletedbecausethesizeofthegradientislessthanthedefaultvalueoftheoptimalitytolerance.<stoppingcriteriadetails>x=1.2671-1.4058options=-10.04832.應(yīng)用Matlab優(yōu)化工具箱求解極值問(wèn)題已知函數(shù)f(x)=10*exp(-x)*cos(x)，求函數(shù)的極值。（1）x∈[2，5]時(shí)，求函數(shù)的最小值，并畫(huà)出函數(shù)的曲線。解：程序?yàn)椋簒=fminbnd('10*exp(-x)*cos(x)',2,5)f='10*exp(-x)*cos(x)';fplot(f,[2,5])y=10*exp(-x)*cos(x)運(yùn)行結(jié)果為：x=2.3562y=-0.6702（2）x∈[3，9]時(shí)，求函數(shù)的最大值，并畫(huà)出函數(shù)的曲線。解：程序?yàn)椋簒=fminbnd('-(10*exp(-x)*cos(x))',3,9)f='10*exp(-x)*cos(x)';fplot(f,[3,9])y=10*exp(-x)*cos(x)結(jié)果為：x=5.4978y=0.0290五．總結(jié)與體會(huì)通過(guò)本實(shí)驗(yàn)熟悉了Matlab優(yōu)化工具箱的主要功能及其使用方法，掌握了優(yōu)化工具箱中常用函數(shù)的功能和語(yǔ)法，能利用其進(jìn)行極值運(yùn)算、求解線性和非線性問(wèn)題，為接下來(lái)的學(xué)習(xí)和實(shí)驗(yàn)進(jìn)一步打下了基礎(chǔ)。

實(shí)驗(yàn)三利用Matlab和Simulink進(jìn)行系統(tǒng)仿真設(shè)計(jì)一．實(shí)驗(yàn)?zāi)康耐ㄟ^(guò)實(shí)驗(yàn)對(duì)一個(gè)汽車(chē)運(yùn)動(dòng)控制系統(tǒng)進(jìn)行實(shí)際設(shè)計(jì)與仿真，掌握控制系統(tǒng)性能的分析和仿真處理過(guò)程，熟悉用Matlab和Simulink進(jìn)行系統(tǒng)仿真的基本方法。二.實(shí)驗(yàn)設(shè)備個(gè)人計(jì)算機(jī)，Matlab軟件。三.實(shí)驗(yàn)準(zhǔn)備預(yù)習(xí)本實(shí)驗(yàn)相關(guān)說(shuō)明，復(fù)習(xí)PID控制器的原理和作用，明確汽車(chē)運(yùn)動(dòng)控制系統(tǒng)問(wèn)題的描述及其模型表示，編寫(xiě)本次仿真練習(xí)的相應(yīng)程序。四.實(shí)驗(yàn)說(shuō)明本實(shí)驗(yàn)是對(duì)一個(gè)汽車(chē)運(yùn)動(dòng)控制系統(tǒng)進(jìn)行實(shí)際設(shè)計(jì)與仿真，其方法是先對(duì)汽車(chē)運(yùn)動(dòng)控制系統(tǒng)進(jìn)行建摸，然后對(duì)其進(jìn)行PID控制器的設(shè)計(jì)，建立了汽車(chē)運(yùn)動(dòng)控制系統(tǒng)的模型后，可采用Matlab和Simulink對(duì)控制系統(tǒng)進(jìn)行仿真設(shè)計(jì)。注意：設(shè)計(jì)系統(tǒng)的控制器之前要觀察該系統(tǒng)的開(kāi)環(huán)階躍響應(yīng)，采用階躍響應(yīng)函數(shù)step()來(lái)實(shí)現(xiàn)，如果系統(tǒng)不能滿足所要求達(dá)到的設(shè)計(jì)性能指標(biāo)，需要加上合適的控制器。然后再按照仿真結(jié)果進(jìn)行PID控制器參數(shù)的調(diào)整，使控制器能夠滿足系統(tǒng)設(shè)計(jì)所要求達(dá)到的性能指標(biāo)。五.實(shí)驗(yàn)內(nèi)容1.問(wèn)題的描述如下圖所示的汽車(chē)運(yùn)動(dòng)控制系統(tǒng)，設(shè)該系統(tǒng)中汽車(chē)車(chē)輪的轉(zhuǎn)動(dòng)慣量可以忽略不計(jì)，并且假定汽車(chē)受到的摩擦阻力大小與汽車(chē)的運(yùn)動(dòng)速度成正比，摩擦阻力的方向與汽車(chē)運(yùn)動(dòng)的方向相反，這樣，該汽車(chē)運(yùn)動(dòng)控制系統(tǒng)可簡(jiǎn)化為一個(gè)簡(jiǎn)單的質(zhì)量阻尼系統(tǒng)。根據(jù)牛頓運(yùn)動(dòng)定律，質(zhì)量阻尼系統(tǒng)的動(dòng)態(tài)數(shù)學(xué)模型可表示為：系統(tǒng)的參數(shù)設(shè)定為：汽車(chē)質(zhì)量m=1000kg，比例系數(shù)b=50N·s/m，汽車(chē)的驅(qū)動(dòng)力u=500N。根據(jù)控制系統(tǒng)的設(shè)計(jì)要求，當(dāng)汽車(chē)的驅(qū)動(dòng)力為500N時(shí)，汽車(chē)將在5秒內(nèi)達(dá)到10m/s的最大速度。由于該系統(tǒng)為簡(jiǎn)單的運(yùn)動(dòng)控制系統(tǒng)，因此將系統(tǒng)設(shè)計(jì)成10％的最大超調(diào)量和2％的穩(wěn)態(tài)誤差。這樣，該汽車(chē)運(yùn)動(dòng)控制系統(tǒng)的性能指標(biāo)可以設(shè)定為：上升時(shí)間：tr<5s；最大超調(diào)量：σ％<10％；穩(wěn)態(tài)誤差：essp<2％。2.系統(tǒng)的模型表示假定系統(tǒng)的初始條件為零，則該系統(tǒng)的Laplace變換式為：即則該系統(tǒng)的傳遞函數(shù)為：如果用Matlab語(yǔ)言表示該系統(tǒng)的傳遞函數(shù)模型，相應(yīng)的程序代碼如下：num=1;den=[100050];G=tf(num,den)結(jié)果為：G=1-----------1000s+50同時(shí)，系統(tǒng)的數(shù)學(xué)模型也可寫(xiě)成如下的狀態(tài)方程形式：如果用Matlab語(yǔ)言表示該系統(tǒng)狀態(tài)空間模型，相應(yīng)的程序代碼如下：A=-1/20;B=1/1000;C=1;sys=ss(A,B,C,0)sys=A=x1x1-0.05B=u1x10.001C=x1y11D=u1y103.系統(tǒng)的仿真設(shè)計(jì)利用Matlab進(jìn)行仿真設(shè)計(jì)=1\*ROMANI．求系統(tǒng)的開(kāi)環(huán)階躍響應(yīng)在Matlab命令窗口輸入相應(yīng)的程序代碼，得出該系統(tǒng)的模型后，接著輸入下面的指令：step(u*sys)可得到該系統(tǒng)的開(kāi)環(huán)階躍響應(yīng)曲線，如下圖所示：從圖上可看出該系統(tǒng)不能滿足系統(tǒng)設(shè)計(jì)所要求達(dá)到的性能指標(biāo)，需要加上合適的控制器。=2\*ROMANII．PID控制器的設(shè)計(jì)PID控制器的傳遞函數(shù)為：在PID控制中，比例（P）、積分（I）、微分（D）這三種控制所起的作用是不同的（請(qǐng)注意在實(shí)驗(yàn)總結(jié)中進(jìn)行歸納）。下面分別討論其設(shè)計(jì)過(guò)程。（1）比例（P）控制器的設(shè)計(jì)增加比例控制器之后閉環(huán)系統(tǒng)的傳遞函數(shù)為：由于比例控制器可以改變系統(tǒng)的上升時(shí)間，現(xiàn)在假定Kp＝100，觀察一下系統(tǒng)的階躍響應(yīng)。在MATLAB命令窗口輸入指令：num=100;den=[1000150];sys=tf(num,den)step(sys)可得系統(tǒng)階躍響應(yīng)如下：由此仿真結(jié)果，分析系統(tǒng)的穩(wěn)態(tài)值是否滿足設(shè)計(jì)要求，系統(tǒng)的穩(wěn)態(tài)誤差和上升時(shí)間能不能滿足設(shè)計(jì)要求？由上圖知系統(tǒng)的穩(wěn)態(tài)誤差為σ%=(0.667-0.6)/0.6*100%=11.17%上升時(shí)間為t,=14.6s均不滿足設(shè)計(jì)要求若減小汽車(chē)的驅(qū)動(dòng)力為10N，重新進(jìn)行仿真，仿真結(jié)果為：如果所設(shè)計(jì)的比例控制器仍不能滿足系統(tǒng)的穩(wěn)態(tài)誤差和上升時(shí)間的設(shè)計(jì)要求，則可以通過(guò)提高控制器的比例增益系數(shù)來(lái)改善系統(tǒng)的輸出。例如把比例增益系數(shù)Kp從100提高到10000重新計(jì)算該系統(tǒng)的階躍響應(yīng)，結(jié)果為：num=10000;den=[100010050];sys=tf(num,den)step(sys)此時(shí)系統(tǒng)的穩(wěn)態(tài)誤差接近為零，系統(tǒng)上升時(shí)間也降到了0.5s以下。這樣做雖然滿足了系統(tǒng)性能要求，但實(shí)際上該控制過(guò)程在現(xiàn)實(shí)中難以實(shí)現(xiàn)。因此，引入比例積分（PI）控制器來(lái)對(duì)系統(tǒng)進(jìn)行調(diào)節(jié)。（2）比例積分（PI）控制器的設(shè)計(jì)采用比例積分控制的系統(tǒng)閉環(huán)傳遞函數(shù)可表示為：增加積分環(huán)節(jié)的目的是減小系統(tǒng)的穩(wěn)態(tài)誤差，假設(shè)比例系數(shù)Kp=600，積分系數(shù)KI=1，編寫(xiě)相應(yīng)的MATLAB程序代碼如下：num=[6001];den=[10006501];sys=tf(num,den)step(sys)運(yùn)行上述程序后可得系統(tǒng)階躍響應(yīng)曲線為：可以調(diào)節(jié)控制器的比例和積分系數(shù)來(lái)滿足系統(tǒng)的性能要求。例如選擇比例系數(shù)KP＝800，積分系數(shù)KI=40時(shí)，可得系統(tǒng)階躍響應(yīng)曲線為：可見(jiàn)，此時(shí)的控制系統(tǒng)已經(jīng)能夠滿足系統(tǒng)要求達(dá)到的性能指標(biāo)設(shè)計(jì)要求。但此控制器無(wú)微分項(xiàng)，而對(duì)于有些實(shí)際控制系統(tǒng)往往需要設(shè)計(jì)完整的PID控制器，以便同時(shí)滿足系統(tǒng)的動(dòng)態(tài)和穩(wěn)態(tài)性能要求。（3）比例積分微分（PID）控制器的設(shè)計(jì)采用PID控制的系統(tǒng)閉環(huán)傳遞函數(shù)為：假設(shè)該控制器的比例系數(shù)KP=1，積分系數(shù)KI=1，微分系數(shù)KD=1，編寫(xiě)MATLAB程序代碼如下：num=[111];den=[10015