控制系統(tǒng)仿真復(fù)習(xí)題及答案

#Simulation—IAMCS)中國(guó)系統(tǒng)仿真學(xué)會(huì)(ChineseAssociationforSystemSimulation—CASS)2MATLAB語(yǔ)言及編程2.1求解下列線性方程，并進(jìn)行解得驗(yàn)證：-721-2-721-2-「4-(1)9153-2x二7,(2)-2-2115-113213057651由A*X=B得：X=A\B「2496_34136x=3614435140156076511087281093791042345解：>>a=[721-2;9153-2;-2-2115;13213]a=TOC\o"1-5"\h\z721-29153-2-2-211513213>>b=[47-10]'b=-10>>x=a\bx=0.49790.14450.0629-0.0813(2)解：>>a=[5765171087268109357910412345]a=5765171087268109357910412345>>b=[2496136361441401560]b=249634136144351401560>>x=a\bx=1.00004.00001.00004.00001.00004.00001.00004.00001.00004.0000進(jìn)行下列計(jì)算，給出不使用for和while等循環(huán)語(yǔ)句的計(jì)算方法。⑴k仝2ii=0解:根據(jù)等比數(shù)列求和方法，在利用matlab中的m文件，編寫程序求解。M文件為n=64;q=2;k=(l-qM)/(l-q);disp('k的值為');disp(k);保存文件q1.m在matlab命令框中輸入>>q1k的值為1.8447e+019(2)求出y=x*sin(x)在0<x<100條件下的每個(gè)峰值解：畫出圖形>>x=0:0.01:100;>>y=x.*sin(x);>>plot(x,y);>>gridon>>title('y=x*sin(x)')>>xlabel('x')>>ylabel('y')

方法1。從圖形中不難看出峰值點(diǎn)取決于函數(shù)sin(x),即在sin(x)為峰值時(shí)，y就得到峰值。所以求取函數(shù)的峰值轉(zhuǎn)化為求取正弦函數(shù)波峰問(wèn)題。而sin(x)8060-40200y-20-40'I-60-80i102030405060708090xI'uTT在x=L+2k“(k為整數(shù)8060-40200y-20-40'I-60-80i102030405060708090xI'uTT在x=L+2k“(k為整數(shù))，所以求取y在上述x時(shí)刻的數(shù)值就是峰值。2y=x*sin(x)在matlab命令行里鍵入100100-1000>>x=pi/2:pi*2:100;>>y=x.*sin(x)%注意是。*不是*%得到結(jié)果y=1.57087.859814.148120.435026.719833.001939.280445.554951.824558.088764.346770.597876.841483.076989.30395.5204方法2.a=size(y)a=11001b=([y(2:1000)]>[y(1:999)])&([y(2:1000)]>[y(3:1001)]);at=find(b==1);disp(y(at))就可以找到最大值點(diǎn)繪制下面的圖形。(1)sin(1/t),-1<t<1(2)1cos3(7t)-1<t<1解:>>t=-1:0.01:1;>>y=sin(1./t);%注意是./不是/%Warning:Dividebyzero.>>plot(t,y)>>gridon>>xlabel('t')>>ylabel('y')>>title('y=sin(1/t)')解：>>t=-1:0.01:1;>>y=l-(cos(7.*t))43;%注意是.*與4%>>plot(t,y)>>gridon>>xlabel('t')>>ylabel('y')>>title('y=1-cos(7t)A3')

已知元件的實(shí)驗(yàn)數(shù)據(jù)如下，擬合這一數(shù)據(jù)，并嘗試給出其特性方程。X0.01001.01002.01003.01004.0100Y2.54377.88849.624211.607111.9727X5.01006.01007.01008.01009.0100y13.218914.267914.613415.404515.0805解：采用最小二乘曲線擬合>>x=0.01:1:9.01;>>y=[2.54377.88849.624211.607111.972713.218914.267914.613415.404515.0805];>>p=polyfit(x,y,3);%選定曲線的階數(shù)為3階，階數(shù)<5,否則曲線不光滑，有數(shù)據(jù)振蕩%>>xi=0:0.01:9.01;>>yi=polyval(p,xi);>>plot(x,y,xi,yi)>>gridon3動(dòng)態(tài)系統(tǒng)數(shù)學(xué)模型及轉(zhuǎn)換線性時(shí)不變動(dòng)態(tài)系統(tǒng)LTI的數(shù)學(xué)模型的微分方程，狀態(tài)方程，傳遞函數(shù)，零極點(diǎn)增益和部分分式五種形式，各有什么特點(diǎn)？答：微分方程是直接描述系統(tǒng)輸入和輸出量之間的制約關(guān)系，是連續(xù)控制系統(tǒng)其他數(shù)學(xué)模型表達(dá)式的基礎(chǔ)。狀態(tài)方程能夠反映系統(tǒng)內(nèi)部各狀態(tài)之間的相互關(guān)系，既適用于單輸入單輸出（SISO）系統(tǒng)，既適用于多輸入多輸出（MIMO）系統(tǒng)。傳遞函數(shù)是零極點(diǎn)形式和部分分式形式的基礎(chǔ)。零極點(diǎn)增益形式可用于分析系統(tǒng)的穩(wěn)定性和快速性。利用部分分式形式可直接分析系統(tǒng)的動(dòng)態(tài)過(guò)程。數(shù)學(xué)模型各種形式之間為什么要互相轉(zhuǎn)換？答：不同的控制系統(tǒng)的分析和設(shè)計(jì)方法，只適用于特定的數(shù)學(xué)模型形式?？刂葡到y(tǒng)建模的基本方法有哪些？他們的區(qū)別和特點(diǎn)是什么？答：控制系統(tǒng)的建模方法大體有三種：機(jī)理模型法，統(tǒng)計(jì)模型法和混合模型法。機(jī)理模型法就是對(duì)已知結(jié)構(gòu)，參數(shù)的物理系統(tǒng)運(yùn)用相應(yīng)的物理定律或定理，經(jīng)過(guò)合理的分析簡(jiǎn)化建立起來(lái)的各物理量間的關(guān)系。該方法需要對(duì)系統(tǒng)的內(nèi)部結(jié)構(gòu)和特性完全的了解，精度高。統(tǒng)計(jì)模型法是采用歸納的方法，根據(jù)系統(tǒng)實(shí)測(cè)的數(shù)據(jù)，運(yùn)用統(tǒng)計(jì)規(guī)律和系統(tǒng)辨識(shí)等理論建立的系統(tǒng)模型。該方法建立的數(shù)學(xué)模型受數(shù)據(jù)量不充分，數(shù)據(jù)精度不一致，數(shù)據(jù)處理方法的不完善，很難在精度上達(dá)到更高的要求?；旌戏ㄊ巧鲜鰞煞N方法的結(jié)合。

控制系統(tǒng)計(jì)算機(jī)仿真中的實(shí)現(xiàn)問(wèn)題是什么含意？答：實(shí)現(xiàn)問(wèn)題就是根據(jù)建立的數(shù)學(xué)模型和精度，采用某種數(shù)值計(jì)算方法將模型方程轉(zhuǎn)換為適合在計(jì)算機(jī)上運(yùn)行的公式和方程，通過(guò)計(jì)算來(lái)使之正確的反映系統(tǒng)各變量動(dòng)態(tài)性能，得到可靠的仿真結(jié)果。用matlab語(yǔ)言求下列系統(tǒng)的狀態(tài)方程、傳遞函數(shù)、零極點(diǎn)增益、和部分分式形式的模型參數(shù)，并分別寫出其相應(yīng)的數(shù)學(xué)模型表達(dá)式：(1)G（s）=s(1)G（s）=s3+7s2+24s+24s4+10s3+35s2+50s+24(2)「2.25-5-1.25-0.5--4「2.25-4.25-1.25-0.252X+0.25-0.5-1.25-121.25-1.75-0.25-0.750uy=[0202]X（1）解：（1）狀態(tài)方程模型參數(shù):編寫matlab程序如下>>num=[172424];>>den=[110355024];>>[ABCD]=tf2ss(num,den)「-10-35-50-24「丁「-10-35-50-24「丁A=1000,B=0010000010_0_得到結(jié)果：,C=h72424〕,D=[0]「-10-35-50-24-丁1000x+0010000010__0_所以模型為:X.=u,y=h72424]x(2)零極點(diǎn)增益:編寫程序>>num=[172424];>>den=[110355024];>>[ZPK]=tf2zp(num,den)得到結(jié)果Z=-2.7306+2.8531,-2.7306-2.8531i,-1.5388P=-4,-3,-2,-1K=1部分分式形式:編寫程序>>num=[172424];>>den=[110355024];>>[RPH]=residue(num,den)得到結(jié)果R=4.0000,-6.0000,2.0000,1.0000TOC\o"1-5"\h\zP=-4.0000,-3.0000,-2.0000,-1.0000H=[]G(s)4—621G(s)=+++—s+4s+3s+2s+1(2)解：(1)傳遞函數(shù)模型參數(shù)：編寫程序>>A=[2.25-5-1.25-0.5-4.25-1.25-0.250.25-0.5-1.25-1-1.75-0.25-0.75]；>>B=[4220]';>>C=[0202];>>D=[0];>>[num>>den]=ss2tf(A,B,C,D)得到結(jié)果num=04.000014.000022.000015.0000den=1.00004.00006.25005.25002.25004s3+14s2+22s+153)S4+4S3+6.25S2+5.25s+2.25(2)零極點(diǎn)增益模型參數(shù):編寫程序>>A=[2.25-5-1.25-0.5-4.25-1.25-0.250.25-0.5-1.25-1-1.75-0.25-0.75]；>>B=[4220]';>>C=[0202];>>D=[0];>>[Z,P,K]=ss2zp(A,B,C,D)得到結(jié)果Z=-1.0000+1.2247i-1.0000-1.2247i-1.5000P=-0.5000+0.8660i-0.5000-0.8660i-1.5000-1.5000K=4.0000表達(dá)式G(s)二4表達(dá)式G(s)二(s+0.5-0.866i)(s+0.5+0.866i)(s+1.5)(3)部分分式形式的模型參數(shù)：編寫程序>>A=[2.25-5-1.25-0.5-4.25-1.25-0.250.25-0.5-1.25-1-1.75-0.25-0.75]；>>B=[4220]';>>C=[0202];>>D=[0];>>[numden]=ss2tf(A,B,C,D)>>[R,P,H]=residue(num,den)得到結(jié)果R=4.0000-0.00000.0000-2.3094i0.0000+2.3094iP=-1.5000-1.5000-0.5000+0.8660i-0.5000-0.8660iH=[]42.3094i2.3094iG(s)=—+—s+1.5s+0.5—0.866is+0.5+0.866i3.6單位反饋系統(tǒng)的開環(huán)傳遞函數(shù)已知如下5s+100s(s+4.6)(s2+3.4s+16.35)

用matlab語(yǔ)句、函數(shù)求取系統(tǒng)閉環(huán)零極點(diǎn)，并求取系統(tǒng)閉環(huán)狀態(tài)方程的可控標(biāo)準(zhǔn)型實(shí)現(xiàn)。解：已知開環(huán)傳遞函數(shù)，求得閉環(huán)傳遞函數(shù)為5s+100s(s+4.6)(s2+3.4s+16.35)+5s+100在matlab命令行里鍵入>>a=[10];>>b=[14.6];>>c=[13.416.35];>>d=conv(a,b)；>>e=conv(d,c)e=1.00008.000031.99000>>f=[0005100];>>g=e+fg=1.00008.000031.9900100.0000%以上是計(jì)算閉環(huán)傳遞函數(shù)的特征多項(xiàng)式%>>p=roots(g)%計(jì)算特征多項(xiàng)式的根，就是閉環(huán)傳遞函數(shù)的極點(diǎn)%p=-0.9987+3.0091i-0.9987-3.0091i-3.0013+0.9697i-3.0013-0.9697i>>m=[5100];

>>z=roots(m)綜上：當(dāng)閉環(huán)傳函形如z=-20%計(jì)算零點(diǎn)%綜上：當(dāng)閉環(huán)傳函形如bs—i+…+b—s+b時(shí)，可控標(biāo)準(zhǔn)型為:sn+asn-i+...+as+a1n-1n廠010...0廠010...0、001...0001—a—a1/nA二；Bbn-1所以可控標(biāo)準(zhǔn)型是x1x.2x.3所以可控標(biāo)準(zhǔn)型是x1x.2x.3x.400一一1001:X1()X2+[X3(9—A—8_—100-80.2131-.400uY二Y二[-100500]+x33.7已知系統(tǒng)的閉環(huán)傳遞函數(shù)?(s)=6s3+26s2+6s+20,試分析該系統(tǒng)的穩(wěn)s4+3s3+4s2+2s+2定性。解：由穩(wěn)定性判據(jù)：當(dāng)閉環(huán)傳遞函數(shù)所有極點(diǎn)都位于虛軸左半平面時(shí)，該系統(tǒng)穩(wěn)定。傳遞函數(shù)的特征方程為:s4+3s3+4s2+2s+2=0,解此方程，得到特征根，即閉環(huán)極點(diǎn)。在matlab命令行里鍵入>>p=[13422];>>r=roots(p)%求多項(xiàng)式等于零的根%得到-1.4734+1.0256i-1.4734-1.0256i-0.0266+0.7873i-0.0266-0.7873i閉環(huán)極點(diǎn)的實(shí)部都小于零，即都位于虛軸左半平面所以系統(tǒng)穩(wěn)定。4動(dòng)態(tài)系統(tǒng)仿真數(shù)值算法動(dòng)態(tài)系統(tǒng)仿真中常用的數(shù)值算法有哪幾類，分別是什么？答：主要有求解線性和非線性微分方程的數(shù)值積分法和計(jì)算線性時(shí)不變動(dòng)態(tài)系統(tǒng)的離散相似法。其中，數(shù)值積分法主要有：歐拉(Euler)法、梯形法、龍格一庫(kù)塔(Runge-Kutta)法和阿達(dá)姆斯(Adams)法；離散相似法主要有：置換法和相似變換法。數(shù)值積分法的選用應(yīng)遵循哪幾條原則？答：數(shù)值積分法應(yīng)該遵循的原則是在滿足系統(tǒng)精度的前提下，提高數(shù)值運(yùn)算的速度和并保證計(jì)算結(jié)果的穩(wěn)定。4.3用圖示法說(shuō)明歐拉(Euler)法的計(jì)算原理答：?jiǎn)栴}設(shè)定設(shè)系統(tǒng)方程為=f(t，y(t))，y(°)=y0

y(t)是一連續(xù)變量t的函數(shù)，現(xiàn)在要以一系列離散時(shí)刻的近似值yl,y2,…,y來(lái)代替，即微分方程的數(shù)值解問(wèn)題。fn+1f(t,y(t))fn+1f(t,y(t))用圖示法說(shuō)明梯形法的計(jì)算原理答：■'預(yù)估:用圖示法說(shuō)明梯形法的計(jì)算原理答：■'預(yù)估:梯形法’校正:y=y+f(t,y)hn+1nnny=y+(f(t,y)+f(t,y))hF+ln2nnn+1n+1ttttttnn+1答：4.5請(qǐng)寫出典型的二階龍格——庫(kù)塔(Runge-Kutta)法公式答：y=y+hf(t+1h,y+1hf)n+1nn2n2ny=y+-r[/+3f(t+-h,y+-對(duì))]n+1n4nn3n3nhy=y+n[f+f(t+h,y+hf)]4.6n+1n2nnnn4.6請(qǐng)寫出至少一種典型的四階龍格——庫(kù)塔(Runge-Kutta)法公式答：hy=答：hy=y+—(k+2k+2k+k),<n+1n61234k=f(t,y)nnhk=f(t+，n2hk=f(t+，3n2k=f(t+h,4nhky+〒)

n2hk

y+22)

n2y+hk)n34.7簡(jiǎn)述數(shù)值積分算法中阿達(dá)姆斯(Adams)法的基本原理答：為了近似計(jì)算積分=yk+「f(t，y)dt,取r+1個(gè)點(diǎn)(tm,ym),(tm-1,ym-1)，?…(tm-r,ym-r)構(gòu)成多項(xiàng)式Pr(t)逼近f(t,y),則可得y=y+1tm+iP(t)dtm+1mtrmy=y+h(bf+bVf+bV2f+...+bVrf)m+1m0m1m2mrm或式中，系數(shù)bj,▽為向后差分。Vf=f-fmmm-1V2f=Vf-Vf=f-2f+fmmm-1mm-1m-2V3f=V2f-V2f=f-3f+f-fmmm-1mm-1m-2m-3在數(shù)值積分算法中，單步法和多步法、顯式法和隱式法、定步長(zhǎng)和變步長(zhǎng)的概念是什么？答：■單步法和多步法：?單步法是指計(jì)算某時(shí)刻數(shù)值yk+1,只需前一時(shí)刻tk有關(guān)信息，它是一種能自啟動(dòng)的算法。?多步法是指計(jì)算某時(shí)刻數(shù)值yk+1需要tk，tk+1時(shí)刻有關(guān)信息，它是一種不能自啟動(dòng)的算法。顯式法和隱式法：顯式是指計(jì)算yk+1時(shí)所需數(shù)據(jù)均已算出。隱式是指計(jì)算yk+1的算式中含有tk+1時(shí)刻的數(shù)據(jù)。因此在使用隱式公式中，需要用另一公式估計(jì)這里未知數(shù)據(jù)，然后用隱式公式進(jìn)行迭代，這叫預(yù)估－校正法。顯然這種方法不能自啟動(dòng)。定步長(zhǎng)和變步長(zhǎng)：?定步長(zhǎng)為積分步長(zhǎng)h在仿真運(yùn)行過(guò)程中始終不變。而積分步長(zhǎng)在仿真運(yùn)行過(guò)程中自動(dòng)修正改變?yōu)樽儾介L(zhǎng)。4.9用歐拉法求下面系統(tǒng)的輸出響應(yīng)y(t)在0<t<上,h=0.1時(shí)的數(shù)值。y'=-y,y(0)=1要求保留4位小數(shù)，并將結(jié)果與真解y(t)=e-t比較。'y=y+h*f(t,y)k豐1kkk解：歐拉法L'=f(t,y)(前向歐拉法,可以自啟動(dòng))其幾何意義：把kky(t)=y00f(t,y)在［t,y］區(qū)間內(nèi)的曲邊面積用矩形面積近似代替。利用matlab提供的mkk文件編程，得到算法公式。如下所示(1)m文件程序?yàn)閔=0.1;disp('函數(shù)的數(shù)值解為')；％顯示中間的文字％disp('y=');%同上%y=1；fort=0:h:1m=y；disp(y)；%顯示y的當(dāng)前值%y=m-m*h；end保存文件q2.m在matalb命令行中鍵入>>q2得到結(jié)果函數(shù)的數(shù)值解為y=10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487(2)另建一個(gè)m文件求解y二e-t在tG［0,1］的數(shù)值(%y二e-1是y'二—y,y(0)二1的真解％)程序?yàn)閔=0.1；disp('函數(shù)的離散時(shí)刻解為');disp('y=')；

fort=0:h:1y=exp(-t);disp(y);end保存文件q3.m在matalb命令行中鍵入>>q3函數(shù)的離散時(shí)刻解為y=10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679比較歐拉方法求解與真值的差別歐10.9000.8100.7290.6560.5900.5310.4780.4300.3870.348拉0001543547直10.9040.8180.7400.6700.6060.5480.4960.4490.4060.367值8783586369誤0-0.00-0.00—0.01—0.01—0.01—0.01—0.01—0.01-0.01-0.01差48071842607483889292顯然誤差與h2為同階無(wú)窮小，歐拉法具有一階計(jì)算精度，精度較低，但算法簡(jiǎn)單。4.10用二階龍格庫(kù)塔法求解4.9的數(shù)值解，并于歐拉法求得的結(jié)果比較。y=y+h(k+k)解：我們經(jīng)常用到預(yù)報(bào)-校正法的二階龍-格庫(kù)塔法，TOC\o"1-5"\h\zk+1k2解：我們經(jīng)常用到預(yù)報(bào)-校正法的二階龍-格庫(kù)塔法，<k=f(t,y)kkk=f(t+h,y+hk)kk1、f(t,y)=y'此方法可以自啟動(dòng)，具有二階計(jì)算精度，幾何意義：把f（t,y）在[t,y]kk區(qū)間內(nèi)的曲邊面積用上下底為f和f、高為h的梯形面積近似代替。kk1利用matlab提供的m文件編程，得到算法公式。如下所示1）m文件程序?yàn)閔=0.1;disp（'函數(shù)的數(shù)值解為'）;disp('y=');y=1;fort=0:h:1disp(y);k1=-y;k2=-(y+k1*h);y=y+(k1+k2)*h/2;end保存文件q4.m在matlab的命令行中鍵入>>q4顯示結(jié)果為函數(shù)的數(shù)值解為y=10.90500.81900.74120.67080.60710.54940.49720.45000.40720.36852）比較歐拉法與二階龍格-庫(kù)塔法求解.（誤差為絕對(duì)值）直/、10.900.810.740.670.600.540.490.440.400.36值48870803658866936679龍10.900.810.740.670.600.540.490.450.400.36庫(kù)50901208719472007285

誤00.000.000.000.000.000.000.000.000.000.00差02030405060606070606明顯誤差為h3得同階無(wú)窮小，具有二階計(jì)算精度，而歐拉法具有以階計(jì)算精度，二階龍格-庫(kù)塔法比歐拉法計(jì)算精度高。4.11用四階龍格-庫(kù)塔法求解題4.9數(shù)值解，并與前兩題結(jié)果相比較。hy=y+—(k+2k+2k+k)k+1k61234k=f(t,y)kk解：四階龍格-庫(kù)塔法表達(dá)式解：四階龍格-庫(kù)塔法表達(dá)式J，其截?cái)嗾`差為h5k=f(t+，y+k)k2k21——k=f(t+，y+k)k2k22k=f(t+—,y+—k)kk3同階無(wú)窮小，當(dāng)h步距取得較小時(shí)，誤差是很小的.(1)編輯m文件程序h=0.1;disp('四階龍格-庫(kù)塔方法求解函數(shù)數(shù)值解為');disp('y=');y=1;fort=0:h:1disp(y);k1=-y;k2=-(y+k1*h/2);k3=-(y+k2*h/2);k4=-(y+k3*h);y=y+(k1+2*k2+2*k3+k4)*h/6;end保存文件q5.mend在matlab命令行里鍵入>>q5得到結(jié)果四階龍格-庫(kù)塔方法求解函數(shù)數(shù)值解為y=10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679（2）比較這幾種方法：對(duì)于四階龍格-庫(kù)塔方法直10.900.810.740.670.600.540.490.440.400.36值48870803658866936679龍10.900.810.740.670.600.540.490.440.400.36庫(kù)48870803658866936679誤00000000000差顯然四階龍格-庫(kù)塔法求解精度很高，基本接近真值。三種方法比較可以得到精度（四階）〉精度（二階）〉精度（歐拉）4.12已知二階系統(tǒng)狀態(tài)方程為XaaXbx(0)X1=11121+1u;1=10XaaXbx(0)X2212222220寫出取計(jì)算步長(zhǎng)為h時(shí)，該系統(tǒng)狀態(tài)變量X=[x,X]的四階龍格-庫(kù)塔法遞12推關(guān)系式。

hy=y+—(k+2k+2k+k)k+1k61234k=f(t,y)kk解：四階龍格-解：四階龍格-庫(kù)塔法表達(dá)式丿k=f(t+，y+k)k2k21——k=f(t+,y+k)k2k22k=f(t+—,y+—k)kk3所以狀態(tài)變量的遞推公式可以寫作：A=a11A=a11a12,B=,x=x1aax212222可以寫成X二AX+Bu—X=X+-*(k+2k+2k+k)TOC\o"1-5"\h\zk+1k61234k=AX+Bu則遞推形式1k則遞推形式k=A(X+k*—/2)+Buk1k=A(X+k*—/2)+Buk2k=A(X+k*—)+Buk34.13用matlab語(yǔ)言編制單變量系統(tǒng)三階龍格-庫(kù)塔法求解程序，程序入口要求能接收狀態(tài)方程各系數(shù)陣(A,B,C,D),和輸入階躍函數(shù)r(t)=R*l(t);程序出口應(yīng)給出輸出量y(t)的動(dòng)態(tài)響應(yīng)數(shù)值解序列y,y,,y。01n解：m文件為：functiony=hs(A,B,C,D,R,T,h)%T為觀測(cè)時(shí)間,h為計(jì)算步長(zhǎng)，R為輸入信號(hào)幅值％disp('數(shù)值解為')；y=0;r=R;x=[0;0;0;0];N=T/h;fort=1:N;k1=A*x+B*R;

k2=A*(x+h*k1/3)+B*R;k3=A*(x+2*h*k2/3)+B*R;x=x+h*(k1+3*k3)/4;y(t)=C*x+D*R;end在命令行里鍵入A=B=C=D=R=T=h=y=hs(A,B,C,D,R,T,h)得到結(jié)果。4.14用梯形法求解試驗(yàn)方程y'一1y，分析對(duì)計(jì)算步長(zhǎng)h有何限制，說(shuō)明Th對(duì)數(shù)值穩(wěn)定性的影響。y二y+h(k+k)k+1k212解：編寫梯形法程序?yàn)榻猓壕帉懱菪畏ǔ绦驗(yàn)閗=-—yTkk二一-(y—-yh)TkTk得到y(tǒng)二得到y(tǒng)二y(1-h+竺)k+1kT2T2穩(wěn)定系統(tǒng)最終漸近收斂。系統(tǒng)穩(wěn)定則1-h+邑T2t2計(jì)算得0<h<2t。h的選取不能超出上述范圍，否則系統(tǒng)不穩(wěn)定。4.15簡(jiǎn)述用數(shù)值積分法仿真LTI動(dòng)態(tài)系統(tǒng)的過(guò)程。答：■將線性時(shí)不變系統(tǒng)模型（n階微分方程、傳遞函數(shù)等）轉(zhuǎn)化為狀態(tài)空間表達(dá)式；X=AX+BUY=CX+DU■由于狀態(tài)方程實(shí)質(zhì)上是一階微分方程組，則可用數(shù)值積分法，如四階RK法可改寫為hX=X+—(K+2K+2K+K)TOC\o"1-5"\h\zm+1m61234K=AX+BU(t)1mmhh{K=A(X+K)+BU(t+)m21m2K=A(X+-K)+BU(t+-)m22m2K=A(X+hK)+BU(t+h)4m3m■系統(tǒng)輸出：二CXm+DSP4.16簡(jiǎn)述置換法的基本原理和簡(jiǎn)單置換法、雙線性置換法的公式。答：■基本原理：用離散算子z去置換連續(xù)傳函中的s，12z_1G(s)口sT*sT+...ors芍吸〒市>G(z)z?1+sT二s=-Z_l■簡(jiǎn)單置換法s=-?二■雙線性置換法Tz+1

4.17簡(jiǎn)述時(shí)域離散相似法基本原理相應(yīng)公式。答：基本原理：時(shí)域離散相似法是在LTI連續(xù)系統(tǒng)狀態(tài)空間表示中，相鄰兩點(diǎn)計(jì)算區(qū)間內(nèi)引入零階保持器或一階保持器后轉(zhuǎn)換為狀態(tài)空間表示的數(shù)值仿真方法。x(t)=x(t)=Ax(t)+Bu(t)y(t)=Cx(t)采樣周期Tu(t)=u(kT),tg[kT,(k+1)T]x=Fx+Gu<k+1kky=HxkkF=eAT,G=feA(T-gdi-B,H=C0(k+1)Tx=x((k+1)T)=eATx(kT)+JeA((k+i)t-i)Bu(T)dT(1)k+1kT5MATLABSimulink應(yīng)用5.1已知系統(tǒng)結(jié)構(gòu)圖如下圖所示。試建立Simulink仿真模型，并確定K值，使阻尼系數(shù)Z=0.5畫出此時(shí)系統(tǒng)的單位階躍響應(yīng)曲線并求出動(dòng)態(tài)性能指標(biāo)（超調(diào)量、上升時(shí)間、調(diào)節(jié)時(shí)間和峰值時(shí)間）。symssG1G2H1H2phi1phiK;G1=16/(s+0.8);H1=K;phi1=G1/(1+G1*H1);G2=phi1*1/s;H2=1;phi=factor(G2/(1+G2*H2))wn=sqrt(16);zeta=0.5;k=(2*zeta*wn-0.8)/16zeta=0.5;wn=4;beta=acos(zeta);tr=(pi-beta)/(wn*sqrt(l-(zeta)A2))tp=pi/(wn*sqrt(l-(zetaF2))ts=3/(zeta*wn)sigma=exp(-pi*zeta/(1-(zeta)A2)A(1/2))分別使用解微分方程方法、控制工具箱、simulink求解具有如下閉環(huán)傳遞函數(shù)的系統(tǒng)的階躍響應(yīng)。0(s)二-10———__—s4+8s3+36s2+40s+10解：(1)用解微分方程方法：將0(s)轉(zhuǎn)化為狀態(tài)方程，利用matlab語(yǔ)句>>num=[10];>>den=[18364010];>>[ABCD]=tf2ss(num,den)得到結(jié)果：A=-8-36-40-10100001000010B=1

10D=0得到狀態(tài)方程x1x.2x.3x.4_-810D=0得到狀態(tài)方程x1x.2x.3x.4_-8-36-40-10_x1x2T1000+00100x30_0010_x4_0_ux110]x2x3x4致！%編寫m文件求解微分方程組functiondx=wffc(t,x)u=1;%階躍響應(yīng)，輸入為1%dx=[-8*x(1)-36*x(2)-40*x(3)-10*x(4)+u;x(1);x(2);x(3)];保存文件wffc.m%注意：保存文件的名字與函數(shù)名一在命令行鍵入>>[t,x]=ode45('wffc',[0,8],[0;0;0;0]);>>y=10*x(:,4);>>plot(t,y);>>grid

得到結(jié)果為下圖所示：2)控制工具箱:在matlab命令行中鍵入>>num=[10];>>den=[18364010];>>sys=tf(num,den);>>step(sys);>>grid得到階躍響應(yīng)結(jié)果如圖所示：0.30.20.10StepResponse用控制〕[具箱求解02468Time(sec)100.30.20.10StepResponse用控制〕[具箱求解02468Time(sec)101214163)simulink求解:在simulink模型窗口中建立如下模型，鍵入該題的傳遞函數(shù)。start后，觀察scope中的仿真波形如下：5.3某小功率隨動(dòng)系統(tǒng)動(dòng)態(tài)結(jié)構(gòu)如圖所示，已知：T二0.01,T二0.05,K二1,K二300,K二120121,K=0.08.c若系統(tǒng)輸入分別為0(t)二1(t),0二t,e二［1(t)-1(1.5)］，適用simulink分析系統(tǒng)srsrsr的輸出0sc(t)分別如何？解：(1)輸入為1(t):輸出為：2）輸入為t時(shí)：輸出為：輸出為：(3)輸入為[1(t)-1(1?5)]:5.4設(shè)典型閉環(huán)結(jié)構(gòu)控制系統(tǒng)如圖所示，當(dāng)階躍輸入幅值R二20時(shí)，用sp4_1.m求取輸出y(t)的響應(yīng)。解：用sp4_1.m求解過(guò)程如下：在MATLAB語(yǔ)言環(huán)境下，輸入以下命令語(yǔ)句>>a=[0.0160.8643.273.421];>>b=[3025];>>X0=[0000];%系統(tǒng)狀態(tài)向量初值為

>>V=2;>>n=4;>>T0=0;Tf=10;>>h=0.01;R=20;躍輸入幅值R二20>>sp4_1>>plot(t,y)運(yùn)行結(jié)果為：%反饋系數(shù)v二2%>>V=2;>>n=4;>>T0=0;Tf=10;>>h=0.01;R=20;躍輸入幅值R二20>>sp4_1>>plot(t,y)運(yùn)行結(jié)果為：%反饋系數(shù)v二2%仿真步長(zhǎng)h=0.01,階%調(diào)用sp4_1.m函數(shù)18I■1614121086421234567891000附：sp4_l.m函數(shù)為b=b/a(1);a=a/a(1);A=a(2:n+1);A=[