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

1、西安電子科技大學 系統(tǒng)仿真上機作業(yè) 指導老師：屈勝利 一 計算機輔助系統(tǒng)分析： 系統(tǒng)如下圖所示其中為單位階躍，為非線性器件 要求：1 當，時，用MATLAB畫出開環(huán)Bode圖，求出、。由其估計出、clear num=0 0 0 40 40; den1=conv(1 0,10 1); den2=conv(0.625 1,0.025 1); den=conv(den1,den2) bode(num,den)grid on;xlabel(Frequency rad/s,fontsize,10)title(Bode Diagram G(s)=40(1+s)/s(10s+1)(0.625s+1)(0.0

2、25s+1) mag,phase,w=bode(num,den)；Gm,Pm,Wcg,Wcp=margin(mag,phase,w)所得參數(shù)為：den = 0.1563,6.5156 ,10.6500 ,1.0000 ,0 增益裕量Gm = 4.3168 相位裕量Pm = 10.0158 穿越頻率Wcg = 5.1598 增益為0的頻率Wcp =2.3975 所以可知：=5.1598 =10.0158計算、程序 sys=tf(num,den) sys=feedback(sys,1); y,t=step(sys); ytr=find(y=1); rise_time=t(ytr(1) ymax,t

3、p=max(y); peak_time=t(tp); max_overshoot=ymax-1 s=length(t); while y(s)0.98&y(s)=1); rise_time=t(ytr(1) ymax,tp=max(y); peak_time=t(tp); max_overshoot=ymax-1 s=length(t); while y(s)0.98&y(s)1.02 s=s-1; end settling_time=t(s)結(jié)果：rise_time =0.7772 max_overshoot = 0.7438 settling_time = 16.9691 可知：=0.77

4、72s =16.9691s =74.38%3，：仿真之，并由仿真結(jié)果求出、答：建立系統(tǒng)模型用自適應(yīng)變步長方法（ode45）仿真可得：觀察上圖可知：超調(diào)量約為80%、上升時間Tr約為0.5S；調(diào)整時間Ts約為18S。定步長RK-2法H=0.08,仿真結(jié)果發(fā)散h=0.05，仿真收斂H=0.02，收斂H=0.001，收斂，但仿真精度明顯降低綜上可知：仿真的步長必須選擇合理才能更準確的得出仿真結(jié)果，步長太大，截斷誤差大，甚至導致仿真失敗，而步長太小，使舍入誤差逐漸積累，也導致誤差增大。所以步長的合理選擇和仿真精度有很大的關(guān)系。4令圖1的分別為：分別仿真，并計算、當GN為飽和特性時，仿真結(jié)果如下： tr

5、=1.37s ts=15.7s、 %=58%。當GN為死區(qū)特性時，仿真結(jié)果如下： tr=0.49s，ts=34.2s、 %=100%。，在之后，反饋點之前加上，仿真之，并計算、仿真圖形： tr=0.43s，ts=17.3s、 %=72%。 對 3 和4 中、的r t 、s t 、 %比較，并解釋差異的原因。比較自適應(yīng)變步長法和定步長RK-2法的仿真結(jié)果得出，RK-2的ts較長，但兩種方法的 tr、 %相差較小。原因：自適應(yīng)變步長法中步長hk的大小與y的變化率有關(guān)，當系統(tǒng)趨于穩(wěn)定時，y變化較慢，步長hk變大，而定步長時hk為一個較小值，故定步長RK-2法的ts較大。而在仿真前期,y變化較劇烈，兩

6、種方法的步長hk都比較小，故其 tr、 %相差不大。若非線性環(huán)節(jié)加在G0(s)之前，死區(qū)特性僅作用于誤差e(s)上，若非線性環(huán)節(jié)加在G0(s)之后，則死區(qū)特性作用于e(s)G0(s)上，在該系統(tǒng)中，G0（S）增益較大，所以后者對原系統(tǒng)的動態(tài)性能影響較大，從而 tr、ts、 %與原系統(tǒng)相差較多；二、病態(tài)系統(tǒng)（stiff）仿真（simulink）:單位階躍:1 用自適應(yīng)變步長法（RK45）仿真之解：當1=102，2=10-2時。仿真結(jié)果如下：2.用定步長四階龍格庫塔法仿真，并試著搜索收斂的步長h 的范圍；若找不到h，將1 增大，2減小，用定步長四階龍格庫塔法仿真，尋找h。解：當1=102，2=10

7、-2時。收斂的步長h的范圍是：h=0.027。當h=0.027時的仿真結(jié)果如下：3.用病態(tài)仿真算法仿真之以上三問，均打印出仿真曲線，計算暫態(tài)響應(yīng)，并比較討論之。答：Y（s）=G(s)R(s)=1s(10s+1)(0.1s+1)=1s-10099s+0.1+199s+10 y(t)=-10099e-0.1t+199e-10t+1。當t=10s時，可算出y(t)=0.39,用自適應(yīng)變步長法（RK45）仿真得到的y(t)=0.164,用定步長四階龍格庫塔法仿真得到的y(t)50，系統(tǒng)嚴重病態(tài)，而RK4-5和RK4若取低步長則仿真速度極慢，若不如此則誤差較大。三計算機輔助控制器設(shè)計：+-ABB要求：開

8、環(huán),且，1 開關(guān)處于A時，系統(tǒng)性能滿足上述要求否？仿真波形：Bode圖為：由仿真波形可以知道，開關(guān)接A時，不能滿足指標。2.開關(guān)處于B時，計算機輔助設(shè)計，使系統(tǒng)性能滿足上述要求（1）基于頻率法的串聯(lián)超前校正主程序：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

9、);hold onstep(Gx_c,b,5);grid onfigure(2)bode(G,r)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) t(2)ans = 1

10、.000000000000000e-03 s(2)ans = 0.140000000000000 theta_theory(2)ans = -1.520837232808838theta_observation(2)ans = -7.075789556117257e-05 theta_output(2)ans =-7.573799703386787e-05(2) 試根據(jù)表格中的數(shù)據(jù)使用MATLAB完成以下問題l 根據(jù)A、B兩列數(shù)據(jù)確定飛機飛行時的理論運行軌跡和飛機的飛行速度；l 根據(jù)A、B、C三列數(shù)據(jù)確定當飛機到達M點時觀測站O到M點的距離。編寫腳本文件plane_rate.mplot(t,s

11、);grid on;%作飛機位移關(guān)于時間t的圖象，可看出飛機是勻速運動i=2:143001;rate(i)=s(i)./t(i);%求運行軌跡中各點的速率，可看出運動速率恒為140m/si=1:143001;abs_theta_theory(i) = abs(theta_theory(i);%取觀察站測得角度的絕對值min_theta_theory k = min(abs_theta_theory);%當飛機到達M點時，觀察站測得角度應(yīng)為最小角，故求得此最小角度時%對應(yīng)的時間t(k)distance_OM = s(k)*cot(abs_theta_theory(1)%當時間為t(k)時，對應(yīng)的

12、位移為s(k)，tan(theta)=PM/OM%其中PM=s(k),theta=theta_theory的最初的值，由此可求得OM的距離運行結(jié)果： planeratedistance_OM = 5.000030000000002e+02故飛機飛行速率為140m/s,OM間距離為500m。（3）設(shè)計標準二階伺服系統(tǒng)G(s) 解：1確定合適的值，其單位階躍響應(yīng)的性能指標滿足:系統(tǒng)的超調(diào)量%介于4.5%8.0%之間因為%=e-1-2*100%，使%介于4.5%8.0%之間則可取=0.68，此時%= %5.4280251339488，滿足要求。 2.使用MATLAB仿真確定參數(shù)n,使得：當把C列數(shù)據(jù)

13、theta_theory作為該二階伺服系統(tǒng)的輸入時，系統(tǒng)的輸出盡可能的接近E列所給出的theta_output。求出系統(tǒng)傳遞函數(shù)參數(shù)。用simulink設(shè)計標準二階伺服系統(tǒng)，其中w的值待定設(shè)計的框圖如下：再編寫腳本文件，如下：for n=1:20;%取不同的w值時，求輸出的theta_output sim(erjiesifuxitong);%運行simulink中的框圖文件， %獲得輸出的的p_theta_output值，為一個長度為143001的列向量 for j=1:143001; d_value(j)=abs(p_theta_output(j)-theta_output(j); %求每組

14、輸出的p_theta_output與理論值的差 end abs_dif(n)=sum(d_value);%取誤差和，存到向量ads_dif中，再依次比較向量中 %的值，當取最小值時說明p_theta_output與理論值最接近，此時 %的w值即為標準w值end當w取120時，得到的abs_dif各值為： abs_difabs_dif = 1.0e+03 * Columns 1 through 4 6.019400510911621 2.739540882427722 1.644088863073810 1.100263457041429 Columns 5 through 8 0.773241

15、725337181 0.549005514111091 0.384051198965172 0.256977609508468 Columns 9 through 12 0.156862800155772 0.079109612117110 0.036344028042700 0.064234872892925 Columns 13 through 16 0.109981064663135 0.149590657196720 0.183731144417703 0.213275684364825 Columns 17 through 20 0.238998060855472 0.2615401

16、30309759 0.281422897899948 0.299065677805618由數(shù)據(jù)可看出，當w=11時，取得最小的差值，再取w=1012之間的數(shù)，求得其中最小差值。得到的simulink框圖為：3.將C列theta_theory作為設(shè)計好的標準二階伺服系統(tǒng)的輸入信號，計算輸出相對于輸入的相對誤差；作出輸出隨時間變化的曲線，以及相對誤差隨時間變化的曲線。作輸出變化曲線和相對誤差曲線求相對誤差的腳本文件為：sim(erjiesifuxitong);%運行二階伺服系統(tǒng)for i=1:143001;abs_dif(i)=abs(p_theta_output(i)-theta_theory(

17、i);%求輸出p_theta_output的絕對誤差relative_error(i)=abs_dif(i)/abs(theta_theory(i);%求輸出的相對誤差End作輸出變化曲線為：從圖中可見系統(tǒng)輸出的曲線與給出的理論曲線幾乎重合，可見選擇的系統(tǒng)參數(shù)比較準確。當作相對誤差曲線時，發(fā)現(xiàn)為如下圖所示：接近于沖激函數(shù)。因為時間軸選用的過寬，無法看出細節(jié)，故用axis(0,2,0,2)指令截取相應(yīng)的曲線，就能看清楚了：a) 時間軸為02時plot(t,relative_error);axis(0,2,0,1)b) 時間軸為270時，圖形幾乎無變化，即相對誤差近似為0.c) 時間軸為7073時

18、，plot(t,relative_error);axis(70,73,0,40)由圖可看出相對誤差變化圖在7073之間有一個猛增的過程d) 時間軸為73143時，plot(t,relative_error);axis(73,143,0,1)分析圖可看出，在70s73s之間，相對誤差有一個顯著的沖激式上升過程，分析原因：相對誤差的求取公式為Poutput-thetatheorythetatheory*100%而theta_theory的值在70s73s時近似取到了零，而theta_theory在分母上，故此時相對誤差值會變得極大。當theta_theory取值遠離零值時，由于絕對誤差接近于0，故

19、相對誤差也接近于0。而在0s2S區(qū)間可看到一個明顯的振蕩過程，分析原因：因為系統(tǒng)為二階系統(tǒng)，輸入的起始階段響應(yīng)會有超調(diào)量，但當經(jīng)過一段時間系統(tǒng)穩(wěn)定后，就不會再發(fā)生振蕩。4.在上述設(shè)計好的伺服系統(tǒng)中加入合適的干擾信號和適當?shù)姆蔷€性環(huán)節(jié)，使得當把C列數(shù)據(jù)theta_theory作為輸入信號時系統(tǒng)的輸出盡可能接近D列數(shù)據(jù)(t)的觀測值theta_observation)。試：1) 確定干擾信號和非線性環(huán)節(jié)的相關(guān)參數(shù)，并畫出系統(tǒng)simulink模型。l 作出theta_theory的時間曲線和theta_observation的時間曲線如圖，藍色為theta_theory的時間曲線，紅色為theta_

20、observation的時間曲線。編寫M文件，使每選定一種非線性環(huán)節(jié)或添加一種干擾都能得到對應(yīng)的曲線，再將對應(yīng)曲線與給定的觀察曲線作比較，即能找到對應(yīng)的參數(shù)。Thetaobservation.m%plot(t,theta_theory);hold on;%作理論曲線plot(t,theta_observation,r);hold on;%作觀察所得曲線sim(erjiesifuxitong);%運行仿真框圖plot(t,p_theta_output);%作仿真結(jié)果曲線hold on;首先在前向通道中加一個死區(qū)非線性環(huán)節(jié)，得到如下圖所示波形：此時輸出曲線與觀察曲線相差較大，故不正確。再將死區(qū)非線

21、性放到回路外：得到曲線：輸出曲線在中間停滯區(qū)與觀察曲線近似，為得到死區(qū)寬度，將題設(shè)所給的理論曲線與觀察曲線打印出來如圖：從上圖中可看出，死區(qū)寬度為1，故將死區(qū)的起始與截止參數(shù)改為如下所示：此時能得到與觀察曲線非線性彎曲部分重合的曲線：此即找到非線性環(huán)節(jié)。l 為找到干擾環(huán)節(jié)，我首先將觀察曲線與輸出曲線作對比，如下圖：為求出合適的干擾信號使藍色曲線接近紅色曲線，先將兩曲線作差(i=1:143001;diff(i)=p_theta_output(i)-theta_observation(i);plot(t,diff);)得到如下曲線：由圖可看出，差值diff的時間曲線近似為一個正弦波形，故應(yīng)輸入一個

22、正弦形狀的干擾信號。為求出干擾信號的具體差值，我再對diff進行一次線性擬合，編寫函數(shù)文件如下：function fun_sin = xianxingnihe(x,t)for i=1:143001 fun_sin(i)=x(1)*sin(x(2)*t(i)+x(3)*pi);end再輸入指令行x(0)=0.1 0.1 0.1x,resnorm=lsqcurvefit(xianxingnihe,x0,t,diff),即對diff曲線進行輸入為時間t的線性擬合，擬合結(jié)果如下：x =-0.049815057234210 0.100054669976175 -0.006617837055889resn

23、orm = 1.820876780757486再作出擬合所得曲線與真實的diff曲線作比較：plot(t,diff);hold on;plot(t,xianxingnihe(x,t),r);得到圖為：由圖可知，擬合結(jié)果近似接近于原曲線，因為原曲線幅值較小，故將此擬合結(jié)果作為干擾信號輸入。輸入系統(tǒng)框圖如圖：其中正弦波參數(shù)采用擬合曲線得到的參數(shù)，如下輸出此框圖的結(jié)果為：由圖可看出，兩條曲線近似重合，故此系統(tǒng)即為標準的輸入系統(tǒng)。2) 畫出系統(tǒng)的輸出隨時間變化的曲線。plot(t,theta_observation);hold on;%作理論觀察曲線sim(erjiesifuxitong);%運行仿真

24、框圖plot(t,p_theta_output,r);%作仿真結(jié)果曲線，為紅色hold on;3) 畫出輸出相對于輸入theta_theory的相對誤差隨時間變化的曲線，以及輸出相對于D列數(shù)據(jù)(t)的觀測值theta_observatuon)的相對誤差曲線。輸出相對于輸入theta_theory的相對誤差隨時間變化的曲線：i=1:143001;diff_a(i)=abs(p_theta_output(i)-theta_theory(i)./abs(theta_theory(i);plot(t,diff_a);hold on;axis(0 143 0 1);l 輸出相對于D列數(shù)據(jù)(t)的觀測值theta_observatuon)的相對誤差曲線i=1:143001;diff_a(i)=abs(p_theta_output(i)-theta_observation(i)./abs(theta_observation(i);plot(t,diff_a);hold on;axis(0 143 0 0.1);五. 幕布系統(tǒng)由電機驅(qū)動電路、電機及幕布三部分組成。通過電機驅(qū)動電路給電機輸入電信號，電機軸轉(zhuǎn)動帶動幕布