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

1、系統(tǒng)仿真上機大作業(yè) 姓名：吳江平 學號：04103043班級：041031老師：屈勝利 一、計算機輔助系統(tǒng)分析：yu+-+re：單位階躍：非線性器件1 ，時，用MATLAB畫出開環(huán)Bode圖，求出、。由其估計出、M程序：%*%功能：畫輸入為單位階躍，開環(huán)傳遞函數(shù)G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)的系統(tǒng)的Bode圖及求增益裕量Gm、相位裕量Pm、相位穿越-180°時的交界頻率Wcg和增益為0時的頻率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)頻率特性的幅值和相角向量和對應的頻率向量magdB=20*log10(mag); %幅值轉換成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)的上升時間，峰值時間，超調量和調整時間%*sys=tf(num,den); %建立系統(tǒng)的開環(huán)傳遞函數(shù)模型sys=feedback(sys,1); %建立系統(tǒng)的閉環(huán)傳遞函數(shù)模型y,t=step(sys); %求出該系統(tǒng)的單位階躍響應ytr=find(y>=1);rise_time=t(ytr(1) %計算上升時間ymax,tp=max(y);peak_time=t(

5、tp) %計算峰值時間max_overshoot=ymax-1 %計算超調量s=length(t);while y(s)>0.98&y(s)<1.02 s=s-1;endsetting_time=t(s+1) %計算調整時間figureh2=plot(t,y,'k',t,ones(length(t),1),'k-.') %繪制響應曲線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)結果：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，時，用MATLAB畫出根軌跡圖，并求出時的閉環(huán)極點；由其估計出、M程序：%*%功能：畫輸入為單位階躍，開環(huán)傳遞函數(shù)G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)的系統(tǒng)的根軌跡，返回閉環(huá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) %返回零極點結果：

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

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

10、,ts,都相對變大。二、病態(tài)系統(tǒng)（stiff）仿真（simulink）:單位階躍1 用自適應變步長法（RK45）仿真之2 用定步長四階RK法仿真，并試著搜索收斂的步長h的范圍步長0.1時發(fā)散0.01發(fā)散0.0020.00294穩(wěn)定。3 用病態(tài)仿真之以上三問，均打印出仿真曲線，計算暫態(tài)響應，并比較討論之1.>,(=0.421，=0.981，=0.941)2.>,(=0.421，=1.012，=0.952)3.>,(=0.426，=0.981，=0.948)三計算機輔助控制器設計：+-ABB要求：開環(huán),且，1 開關處于A時，系統(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) %計算上升時間ymax,tp=max(y);max_overshoot=ymax-1 %計算超調量s=length(t);while y(s)>0.98&y(s)<1.02 s=s-1;endsetting_time=t(s+1) %計算調整時間bode(num,den)grid onmag phase w=bode(num,den); %返回系統(tǒng)頻率特性的幅值和相角向量和對應的頻率向量

12、Gm Pm Wcg Wcp=margin(mag,phase,w) %求增益裕量Gm、相位裕量Pm、相位穿越-180°時的交界頻率Wcg和增益為0時的頻率Wcp結果：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由此可以看出，開關處于A時，系統(tǒng)性能顯然不能滿足以下要求：開環(huán),且，截圖：2 開關處于B時，計算機輔助設計

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);結果：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（上升時間）max_over = 0.2235（超調量）ts =1.2279s（調整時間）截圖：校正前后階躍響應曲線：校正前后bode圖：注釋：藍色線表示校正后，紅色線表示校正前上機作業(yè)題四仿真在位移x處加階躍干擾的倒立擺非線性模型（已建好，但仿真不成功）

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