系統(tǒng)仿真上機(jī)大作業(yè)

1、系統(tǒng)仿真上機(jī)大作業(yè) 姓名：吳江平 學(xué)號(hào)：04103043班級(jí)：041031老師：屈勝利 一、計(jì)算機(jī)輔助系統(tǒng)分析：yu+-+re：?jiǎn)挝浑A躍：非線性器件1 ，時(shí)，用MATLAB畫出開環(huán)Bode圖，求出、。由其估計(jì)出、M程序：%*%功能：畫輸入為單位階躍，開環(huán)傳遞函數(shù)G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)的系統(tǒng)的Bode圖及求增益裕量Gm、相位裕量Pm、相位穿越-180°時(shí)的交界頻率Wcg和增益為0時(shí)的頻率Wcp%*num=40 40;den=conv(conv(conv(1 0,10 1),0.625 1),0.025 1);w=logspac

2、e(-3,3,100); mag phase w=bode(num,den,w); %返回系統(tǒng)頻率特性的幅值和相角向量和對(duì)應(yīng)的頻率向量magdB=20*log10(mag); %幅值轉(zhuǎn)換成dB值subplot(2,1,1)semilogx(w,magdB,'k-','LineWidth',2) %畫幅頻特性Bode圖grid onxlabel('Frequecy(rad/s)','FontSize',14)ylabel('Gain(dB)','FontSize',14)h1=title('B

3、ode Diagram of G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)')set(h1,'FontWeigh','bold','FontSize',12)subplot(2,1,2)semilogx(w,phase,'k-','LineWidth',2) %畫相頻特性Bode圖grid onxlabel('Frequecy(rad/s)','FontSize',14)ylabel('Phase(deg)','F

4、ontSize',14)Gm Pm Wcg Wcp=margin(mag,phase,w)%*%功能：求輸入為單位階躍，開環(huán)傳遞函數(shù)G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)的系統(tǒng)的上升時(shí)間，峰值時(shí)間，超調(diào)量和調(diào)整時(shí)間%*sys=tf(num,den); %建立系統(tǒng)的開環(huán)傳遞函數(shù)模型sys=feedback(sys,1); %建立系統(tǒng)的閉環(huán)傳遞函數(shù)模型y,t=step(sys); %求出該系統(tǒng)的單位階躍響應(yīng)ytr=find(y>=1);rise_time=t(ytr(1) %計(jì)算上升時(shí)間ymax,tp=max(y);peak_time=t(

5、tp) %計(jì)算峰值時(shí)間max_overshoot=ymax-1 %計(jì)算超調(diào)量s=length(t);while y(s)>0.98&y(s)<1.02 s=s-1;endsetting_time=t(s+1) %計(jì)算調(diào)整時(shí)間figureh2=plot(t,y,'k',t,ones(length(t),1),'k-.') %繪制響應(yīng)曲線set(h2(1),'LineWidth',2)xlabel('Time(sec)','FontSize',14)ylabel('Reponse',

6、'FontSize',14)h2=title('Plot of Unit-step Reponse of G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)');set(h2,'FontWeigh','bold','FontSize',12)結(jié)果：Gm Pm Wcg Wcp=4.3221 10.0223 5.1653 2.3981rise_time peak_timel max_overshoot setting_time=0.7772 1.2953 0.8168 17.2280所以

7、 =Wcp=2.3981，b=Pm=10.0223 =0.7772，=17.2280，=0.8168截圖：2，時(shí)，用MATLAB畫出根軌跡圖，并求出時(shí)的閉環(huán)極點(diǎn)；由其估計(jì)出、M程序：%*%功能：畫輸入為單位階躍，開環(huán)傳遞函數(shù)G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)的系統(tǒng)的根軌跡，返回閉環(huán)零極點(diǎn)%*num=40 40;den=conv(conv(conv(1 0,10 1),0.625 1),0.025 1);sys=tf(num,den); rlocus(sys)sys=feedback(sys,1); p,z=pzmap(sys) %返回零極點(diǎn)結(jié)果：

8、p = -40.1616 -0.2274 + 2.4146i -0.2274 - 2.4146i -1.0837 =0.7772，=17.2280，=0.8168截圖： 3、，時(shí)，仿真之，并由仿真結(jié)果求出，自適應(yīng)變步長(zhǎng)法M程序：sim('f01031'); %進(jìn)行模型仿真plot(ScopeData(:,1),ScopeData(:,2) %ScopeData為示波器輸出到工作空間的矩陣，其中第一列為仿真時(shí)間，第二列為仿真波形grid onh=title('自適應(yīng)變步長(zhǎng)ode45法');set(h,'FontWeigh','bold

9、9;,'Color','red,'FontSize',20)結(jié)果： =0.685，=17.14，=0.8173截圖：定步長(zhǎng)二階RK法截圖：當(dāng)步長(zhǎng)大于或等于0.050時(shí)發(fā)散，否則收斂。4令圖1的分別為：分別仿真之，并計(jì)算、（從仿真結(jié)果）結(jié)果： =0.895，=15.14，=0.7142結(jié)果： =0.746，=19.60，=0.8193，在之后，反饋點(diǎn)之前加上，仿真之，并計(jì)算、結(jié)果： =0.721，=14.94，=0.8973對(duì)3；4、的、比較，并解釋差異的原因 自適應(yīng)變步長(zhǎng)比用定步長(zhǎng)法tr,ts,小，因?yàn)樽赃m應(yīng)變步長(zhǎng)可自動(dòng)調(diào)整步長(zhǎng)。受非線性環(huán)節(jié)的影響，tr

10、,ts,都相對(duì)變大。二、病態(tài)系統(tǒng)（stiff）仿真（simulink）:單位階躍1 用自適應(yīng)變步長(zhǎng)法（RK45）仿真之2 用定步長(zhǎng)四階RK法仿真，并試著搜索收斂的步長(zhǎng)h的范圍步長(zhǎng)0.1時(shí)發(fā)散0.01發(fā)散0.0020.00294穩(wěn)定。3 用病態(tài)仿真之以上三問(wèn)，均打印出仿真曲線，計(jì)算暫態(tài)響應(yīng)，并比較討論之1.>,(=0.421，=0.981，=0.941)2.>,(=0.421，=1.012，=0.952)3.>,(=0.426，=0.981，=0.948)三計(jì)算機(jī)輔助控制器設(shè)計(jì)：+-ABB要求：開環(huán),且，1 開關(guān)處于A時(shí)，系統(tǒng)性能滿足上述要求否？M程序：num=10;den=1

11、 1 0;sys=tf(num,den); sys=feedback(sys,1); y,t=step(sys); ytr=find(y>=1);rise_time=t(ytr(1) %計(jì)算上升時(shí)間ymax,tp=max(y);max_overshoot=ymax-1 %計(jì)算超調(diào)量s=length(t);while y(s)>0.98&y(s)<1.02 s=s-1;endsetting_time=t(s+1) %計(jì)算調(diào)整時(shí)間bode(num,den)grid onmag phase w=bode(num,den); %返回系統(tǒng)頻率特性的幅值和相角向量和對(duì)應(yīng)的頻率向量

12、Gm Pm Wcg Wcp=margin(mag,phase,w) %求增益裕量Gm、相位裕量Pm、相位穿越-180°時(shí)的交界頻率Wcg和增益為0時(shí)的頻率Wcp結(jié)果：Gm Pm Wcg Wcp= 4.0007e+003 18.0184 200.0277 3.0832rise_time max_overshoot setting_time= 0.5961 0.6042 7.3516所以 =Wcp=3.0832，b=Pm=18.0184 =0.5962，=7.3516，=0.6042由此可以看出，開關(guān)處于A時(shí)，系統(tǒng)性能顯然不能滿足以下要求：開環(huán),且，截圖：2 開關(guān)處于B時(shí)，計(jì)算機(jī)輔助設(shè)計(jì)

13、，使系統(tǒng)性能滿足上述要求?；陬l率法的串聯(lián)超前校正M程序：主程序：function=question3_2(k)num=k*1;den=conv(1,0,1,1);G=tf(num,den);kc=1;yPm=45+10;Gc=plsj(G,kc,yPm) %超前校正環(huán)節(jié)Gy_c=feedback(G,1) %校正前系統(tǒng)閉環(huán)傳遞函數(shù)Gx_c=feedback(G*kc*Gc,1) %校正后系統(tǒng)閉環(huán)傳遞函數(shù)figure(1)step(Gy_c,'r',5);hold onstep(Gx_c,'b',5);grid onfigure(2)bode(G,'r&

14、#39;)hold on bode(G*kc*Gc,'b')grid onGm,Pm,Wcg,Wcp=margin(G*kc*Gc);r=PmWc=Wcpy,t=step(Gx_c);ytr=find(y>=1);rise_time=t(ytr(1);tr=rise_timeymax,tp=max(y); peak_time=t(tp); max_overshoot=ymax-1; max_over=max_overshoot s=length(t);while y(s)>0.98&y(s)<1.02 s=s-1;endsettling_time=t(

15、s+1);ts=settling_time 子程序：function Gc=plsj(G,kc,yPm)G=tf(G);mag,pha,w=bode(G*kc);Mag=20*log10(mag);Gm,Pm.Wcg,Wcp=margin(G*kc);phi=(yPm-getfield(Pm,'Wcg')*pi/180;alpha=(1+sin(phi)/(1-sin(phi);Mn=-10*log10(alpha);Wcgn=spline(Mag,w,Mn);T=1/(Wcgn*sqrt(alpha);Tz=alpha*T;Gc=tf(Tz,1,T,1);結(jié)果：Transfe

16、r function:0.4536 s + 1-0.1126 s + 1 Transfer function: 10-s2 + s + 10 Transfer function: 4.536 s + 10-0.1126 s3 + 1.113 s2 + 5.536 s + 10R = 49.7706Wc = 4.4248tr = 0.3982s（上升時(shí)間）max_over = 0.2235（超調(diào)量）ts =1.2279s（調(diào)整時(shí)間）截圖：校正前后階躍響應(yīng)曲線：校正前后bode圖：注釋：藍(lán)色線表示校正后，紅色線表示校正前上機(jī)作業(yè)題四仿真在位移x處加階躍干擾的倒立擺非線性模型（已建好，但仿真不成功）

17、，可以修改各種仿真參數(shù)，如算法、誤差限、步長(zhǎng)、時(shí)間等，但是不能修改模型結(jié)構(gòu)，要求仿真成功，說(shuō)明修改過(guò)程，解釋物理意義，并說(shuō)明如何驗(yàn)證結(jié)果不穩(wěn)定是系統(tǒng)本身不穩(wěn)定還是仿真不當(dāng)造成的？備注：以下是倒立擺系統(tǒng)建模過(guò)程。倒立擺系統(tǒng)如圖1所示，它由在導(dǎo)軌上運(yùn)動(dòng)的小車和與小車相連的擺桿組成。圖 1 倒立擺系統(tǒng)我們假設(shè)：M：小車的質(zhì)量 0.5kgm：擺桿的質(zhì)量 0.2kgb：小車與基座之間的摩擦系數(shù) 0.1N/m/sl：小車與擺桿連接點(diǎn)到擺桿重心之間的距離 0.3mI：擺桿的轉(zhuǎn)動(dòng)慣量 0.006F：作用在小車上的力 x：小車的位置座標(biāo)（水平）：擺桿與豎直方向的夾角 將倒立擺系統(tǒng)分為兩個(gè)獨(dú)立的部分：小車和擺桿，分別如圖2和圖3所示。圖 2 小車圖 3 擺桿由小車水平方向的受力情況可以得到下面等式： Error! No bookmark name given.Error! No bookmark name given.由擺桿水平方向的受力情況可以得到下面等式： 將代入中得：Error! No b