控制系统仿真

1、-an-1-an-2-a1-a0011001000011 00nxxuaaabyx0)0(,)(,()(yytytfty1)(,()t ()t (11nnttnnnndttytfyyyy1( , ( ) ( , ( )( )nntnnnnntf t y t dthfff ty ty t，11( , ( )nntnntyyf t y t dty(t) t0yn+1y(tn+1)yn n n+1f(tn,yn)f(t,y(t)t0fn+1fn n n+1hf(tn,yn)1)0(02yyy 1( )1y tt1111( ,) 1( ( ,)(,)2nnnnnnnnnnyyf ty hyyf tyf

2、 tyh预估：校正：ffn+1fntn+1tnht00.511.522.530.20.30.40.50.60.70.80.91ty(t) 精 确 解欧 拉 法梯 形 法多 步 梯 形 法% 欧拉法、梯形法求解微分方程比较欧拉法、梯形法求解微分方程比较close all;clear all;h = 0.4; t = 0:h:3; tn,tm = size(t); y = zeros(size(t); yo = y; yt = y; ytn = y; for i = 1:tm y(i) = 1/(1+t(i);end %欧拉法欧拉法yo(1) = 1;for i = 2:tm yo(i) = yo

3、(i-1)+ h*(-yo(i-1)2 );end %梯形法梯形法yt(1)=1;for i = 2:tm yy = yt(i-1)+h*(- yt(i-1)2 ); yt(i) = yt(i-1)+ h*(-yt(i-1)2+(-yy2) )/2;end %多步梯形法多步梯形法ytn(1)=1;yy = 0;for i = 2:tm yy = ytn(i-1)+h*(- ytn(i-1)2 ); while abs(ytn(i)-yy) 0.00001 yy = ytn(i-1)+ h*(-ytn(i-1)2+(-ytn(i)2) )/2; ytn(i) = ytn(i-1)+ h*(-yt

4、n(i-1)2+(-yy2) )/2; endend plot(t,y,t,yo,t,yt,t,ytn)Legend(精确解精确解,欧拉法欧拉法,梯形法梯形法,多步梯形法多步梯形法) xlabel(t); ylabel(y(t)yyx_Euler01 .m，并且，并且y在在a,b内连续可微。内连续可微。0( )( , ( ),( )y tf t y tatby ay 1212( ,)( ,)f t yf t yL yy（1）11),(kkkkdtytfyy011NNattttbiy1iiihtt23111.,0,1,2,.2!3!kkkkkyyyhyhyhk（2）（3）（3）23111.,0,

5、1,2,.2!3!kkkkkyyyhyhyhk00)()(,(ytytytfdtdy（4）1kkkyyyh （5）001000021111( ,)( ,).( ,)t ty ykkkkyyyhyf ty hyyf t y hyyf ty h（6）23211.2!3!()nnRyhyhRO h1()pO h()pO h111( ( ,)(,)2kkkkkkhyyf tyf ty1111( ,) 1( ( ,)(,)2kkkkkkkkkkyyf ty hyyf tyf tyh预估：校正：(8)(9)(0)1( ,)kkkkyyf ty h(1)(0)111( )(1)1111( ( ,)(,)2.

6、1( ( ,)(,)2kkkkkkiikkkkkkyyf tyf tyhyyf tyf tyh直到( )(1)11iikkyy为止212121( ,)( ,)2kkkkkkkkkkkyyyhyhffyhf tyf tyhty1212111122( ,)(,)()kkkkkkKf tyKf tbh yb K hyyh a Ka K(12)(11)11212(,)(,)(,)kkkkkkkkkkffyya hf tya h f tyb hb hf tyty2121( ,)kkkkffKf tyb hb K hty122 12211212aaa ba b(13)(14)(15)211( ,)( ,)

7、2kkkkkkkkffyyhf tyf tyhty(16)(17)(18)1212121bbaa112121()2( ,)(,)kkkkkkhyyKKKf tyKf th yK h3411.3!4!nRyhyh3()nRO h或一次梯形法一次梯形法111(,)211(,)22223 (,)433kkkkkkkkkkkkkkkkkhyyff thyhfyyhf thyhfhyyff thyhftktk+1yk+1ykhf(t,y(t)f(k+1,yk+1)f(k,yk)k k+h/2 k+1t0hf(k+h/2, y(k+h/2)f(t,y(t)f(k+1,yk+1)f(k,yk)k h k+1

8、t0hf(k, y(k)+f(k+1, y(k+1)/20)1 (2xxxmx xyxy21,1221221)1 (,yyymyyy0)0(, 1)0(21yy051015202530-4-3-2-101234 Time , Second Y(1) and Y(2) Van der Pol Solution for mu=2 yyx_euler03.mglobal MU % define MU as a global variable in the Command WorkspaceMU=2; % set global parameter to desired valuet,y=ode23(y

9、yxvdpol,0 ;30,1 ; 0); % to=0 , tf=30 , yo=1 ; 0y1=y(:,1); % first column is y(1) versus time points in ty2=y(:,2); % second column is y(2)plot(t , y1 , t , y2 , - )xlabel( Time , Second ) , ylabel( Y(1) and Y(2) )title(Van der Pol Solution for mu=2 )function yprime=yyxvdpol(t , y)global MU % choose

10、0MU10 in Command workspaceyprime=y(2); MU*(1-y(1)2)*y(2)-y(1); 11113223(,)(3),(,)43322(,)33kkkkkkkkkf tyhkhhyykkkf tyhkhkf ty112112342343(,)(,)22(22),6(,)22(,)kkkkkkkkkkkf tyhkhkf tyhyykkkkhkhkf tykf thyhk1111(,)1,2,.,rkkiiirikikijjjyyK hKf ta h yb Kir（19）1( ,)kkkkyyf ty h11122122211( ,)(,)kkkkkkyyK

11、 hK hKf tyKf tc h ya K（20）1123121311(4)6( ,)1(,)22(,2)kkkkkkmmhyyKKKKf tyhKf tyK hKf th yK hK h（21）112341213243(22)6( ,)1(,)221(,)22(,)kkkkkkkkkkhyyKKKKKf tyhKf tyhKhKf tyhKKf th yhK12121( ,)(,)22kkkkkkyyhKKf tyhhKf tyK（22）112341213124123(33)8( ,)(,)332(,(3)33(,)kkkkkkkkkkhyyKKKKKf t yhhKf tyKhKf th

12、 yKKKf th yhKhKhK11kkkEyy（23）（24）611kkiiiyyhc K11(,)iikikijijKf ta h yhb K6,.,2, 1i（25）11),(kkkkdtytfyy11( )kktkkrtyyP t dt21012(.)rkkkkkrkyyh b fbfbfbf （26）（27）12112322112323kkkkkkkkkkkkkkkkffffffffffffffff 一阶向后差分：二阶向后差分：三阶向后差分：1011(.)kkkkrkryyhfff（28）（27）11),(kkkkdtytfyy10111(.)rkkkkrkyyh c fcfcf（

13、29）111011(.)kkkkrkryyh r fr frf （30）(2)2( ),(31)rrkrEB yh)()2(ry(0)101(1)( )( )( )101111(.),(32)(.),(33)rkkkkrkiiirikkkkrkyyh b fbfbfyyh c fc fcf显示预估：隐式预估：( )( )111( )1( )111( ,),(,),kkkiikkkjijijkkkff tyff tyfff （34）（35）（36）(0)1011(1)( )111011(.)(.)kkkkrkriikkkkkryyhfffyyh rfr frf ，（37），（38）(0)1123

14、(0)( )1112(5559379)24(9195)24kkkkkkikkkkkkhyyffffhyyffff（39）（40）10rth 或ch51，（41），（42）maxkkEy，（43）minkminmaxkmaxk1max1()Kkkkhh10.110,kkkhhh（44）（45）maxkmaxk)( tydtdy0j（46）（47）)(1mmmthfyymmmhyyy1)()1 ()(zYhzzY，（50）11hz11)(hja，（53） ，（55）2222221111)2()(hhahahh1h1ah2Th 2，（54），（57）aT1111hz（6060）即22211)(hha

15、（61）式（）表明，在 平面上，欧拉后差法稳定区域是一个圆外区域，该圆的圆心为)0 ,1(h，半径为h1，如图 5 所示。 （61）),.,(.),.,(),.,(2121222111nnnnnxxxtfxxxxtfxxxxtfx(68),(XtFX (69)TnxxxX,.,21),.,(.),.,(),.,(),(21212211nnnnxxxtfxxxtfxxxtfXtF112121()2( ,)(,)kkkkkkhXXKKKF tXKF th XK h112341213243(22)2( ,)1(,)22(,)22(,)kkkkkkkkkkhXXKKKKKF tXhKF tXhKhhK

16、F tXKKF th XhKDUCXYBUAXX（70）TnxxxX,.,21112341213243(22)6( )()()22()()22()()kkkkkkkkkkhXXKKKKKAXBU thhKA XKBU thhKA XKBU tKA XhKBU th( )kkkYCXDU t例：例：yyx_siso_step1.ms = tf(s)G1 = 5/(2*s2+s+1)y,t = yyx_siso_step2(G1,25,0.2,5);), 0, (,)0, (,optionsytspanFsolverYTytspanfsolverYT123121323123(0)0(0)10.51

17、(0)1yy yyyy yyyy yy options = odeset(RelTol,1e-4,AbsTol,1e-4 1e-4 1e-5);T,Y = ode45(rigid,0 12,0 1 1,options);figure;grid onplot(T,Y(:,1),-,T,Y(:,2),-.,T,Y(:,3),.)figure;grid onplot3(Y(:,1),Y(:,2),Y(:,3),r)function dy = rigid(t,y)dy = zeros(3,1); % a column vectordy(1) = y(2) * y(3);dy(2) = -y(1) *

18、y(3);dy(3) = -0.51 * y(1) * y(2);024681012-1.5-1-0.500.511.5-101-2020.70.80.911.1)()()()()() 1(kDUkCXkYkGUkFXkX)()(112ln11zGsGzzTzTsorsTezsTTzssTz1112 . 01)()()(2sssUsYsG)2 . 01 ()2 . 02()()()(222TTzTzTzUzYzG)()()2 . 01 () 1()2 . 02()2(22kuTkyTTkyTkyyyx_ZH01.mTkTyTkydttdykTt)() 1()(T1T减小减小S平面和稳定域平面和

19、稳定域平面和稳定域平面和稳定域Tz1121ln1zszTTz112zzTs2121TsTsz)() 1(2)2()()4 . 04() 1()82()2()4 . 04(222222kuTkuTkuTkyTTkyTkyTT)()()()(2211zGsGzGsG)()()()(2121zGzGsGsGsesGsT1)(0TktkTkTutuk) 1(),()(g(t)01T211)1 ()( TseTsTsGsTTktkTkTtTTkukTukTutuk) 1(),() 1()()()(零阶保持器对信号的响应零阶保持器对信号的响应一阶保持器对信号的响应一阶保持器对信号的响应ssesGsGsHs

20、T11)()()(0sTseTsTsGsGsHsT11)1 ()()()(21 1)(zTzH) 1(1-32)(zzzTzH)()() 1(kTukyky) 1(32)() 1(kukuTkyky）)()1 ()(aTaTezaeAzH(1)( )(1) (aTaTAy key keu ka）1010)(AsABsBsG)()1 () 1()(101010aTaTezAAeBAzABzH01/ AAa )(1()() 1(kCukBukAyky）aTeA10AAD DABB/10DABeBABACaT/ )(101010541)(2ssssG0.6703 + z 1.629 - z0.085

21、57 - z 0.07736H(z) 20.6703 + z 1.629 - z0.03954 - z 0.01093 - z 0.04226H(z) 22Time (sec.)AmplitudeStep Response00.511.522.53-0.25-0.2-0.15-0.1-0.0500.050.10.15From: U(1)To: Y(1)Time (sec.)AmplitudeStep Response00.511.522.53-0.25-0.2-0.15-0.1-0.0500.050.10.15From: U(1)To: Y(1)yyx_C2D01.mclose all;cle

22、ar all;G = tf(1 -1,1 4 5)H = c2d(G,0.1)step(G,H);figureH = c2d(G,0.1,foh)step(G,H);)()()()()(tCxtytButAxtx kkkkkHxyGuFxx1采样周期采样周期T) 1( ,),()(TkkTtkTutuTTAATCHBdeGeF0)(,) 1 ()()() 1()1()1(1dBuekTxeTkxxTkkTTkAATk作业()()()()()10001()()A TA TTTTA TATA TA TATATdeAeddeIeA ededAeIdGAeI B 例：例：yyx_siso_step2.

23、ms = tf(s)G1 = 5/(2*s2+s+1)y,t = yyx_siso_step1(G1,25,0.2);)()()1()()()()1(,)1()2()3()3()()0()1()2()()0()()()0()(0)()1(0)1()1(0)(0)(kTudBekTxeTkxkTuuTkkTedBuexeTkxdBuexekTxdBuexetxTTAATATTkTkATkAkTkTAAkTttAAT采用零阶保持器上，在式。得F = 0.54881 0 0 -0.18373 0.65023 -0.39463 -0.052061 0.39463 0.89981G = 0.22559

24、0.21666 0.057712C = 0 0 1.2649D = 0Time (sec.)AmplitudeStep Response02468101200.10.20.30.40.50.60.7From: U(1)To: Y(1)Time (sec.)AmplitudeStep Response02468101200.10.20.30.40.50.60.7From: U(1)To: Y(1)yyx_C2D02.mclose all;clear all; Ts = 0.3;Z = -1;P=-2,-0.5+1.5i,-0.5-1.5j;K=2;G = zpk(Z,P,K)A,B,C,D =

25、ZP2SS(Z,P,K) sys = ss(A,B,C,0)sysd = c2d(sys,Ts) step(sys,sysd); figuresysd = c2d(sys,Ts,foh) step(sys,sysd);%impulse(sys,-,sysd,-);BuAXXXAXBu()0( )( )( )0A ttAtX te XeBudduBeXekTXkTAkTAkT)()0()()(0duBeXeTkXTkATkTkA)()()(01)1()1(0)1(duBekTXeTkXTkATkkTAT)(1)1()1()()()0()()()()1A tTATX kTe X kTeBdu kTFX kTGu kTTkTuTkukTuu)()()()(1BuAXXTkTuTk