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

1、西安電子科技大學(xué) 系統(tǒng)仿真上機(jī)作業(yè) 指導(dǎo)老師：屈勝利 一 計(jì)算機(jī)輔助系統(tǒng)分析： 系統(tǒng)如下圖所示其中為單位階躍，為非線性器件 要求：1 當(dāng)，時(shí)，用MATLAB畫(huà)出開(kāi)環(huán)Bode圖，求出、。由其估計(jì)出、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計(jì)算、程序 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)變步長(zhǎng)方法（ode45）仿真可得：觀察上圖可知：超調(diào)量約為80%、上升時(shí)間Tr約為0.5S；調(diào)整時(shí)間Ts約為18S。定步長(zhǎng)RK-2法H=0.08,仿真結(jié)果發(fā)散h=0.05，仿真收斂H=0.02，收斂H=0.001，收斂，但仿真精度明顯降低綜上可知：仿真的步長(zhǎng)必須選擇合理才能更準(zhǔn)確的得出仿真結(jié)果，步長(zhǎng)太大，截?cái)嗾`差大，甚至導(dǎo)致仿真失敗，而步長(zhǎng)太小，使舍入誤差逐漸積累，也導(dǎo)致誤差增大。所以步長(zhǎng)的合理選擇和仿真精度有很大的關(guān)系。4令圖1的分別為：分別仿真，并計(jì)算、當(dāng)GN為飽和特性時(shí)，仿真結(jié)果如下： tr

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

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

7、-2時(shí)。收斂的步長(zhǎng)h的范圍是：h=0.027。當(dāng)h=0.027時(shí)的仿真結(jié)果如下：3.用病態(tài)仿真算法仿真之以上三問(wèn)，均打印出仿真曲線，計(jì)算暫態(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。當(dāng)t=10s時(shí)，可算出y(t)=0.39,用自適應(yīng)變步長(zhǎng)法（RK45）仿真得到的y(t)=0.164,用定步長(zhǎng)四階龍格庫(kù)塔法仿真得到的y(t)50，系統(tǒng)嚴(yán)重病態(tài)，而RK4-5和RK4若取低步長(zhǎng)則仿真速度極慢，若不如此則誤差較大。三計(jì)算機(jī)輔助控制器設(shè)計(jì)：+-ABB要求：開(kāi)

8、環(huán),且，1 開(kāi)關(guān)處于A時(shí)，系統(tǒng)性能滿足上述要求否？仿真波形：Bode圖為：由仿真波形可以知道，開(kāi)關(guān)接A時(shí)，不能滿足指標(biāo)。2.開(kāi)關(guān)處于B時(shí)，計(jì)算機(jī)輔助設(shè)計(jì)，使系統(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完成以下問(wèn)題l 根據(jù)A、B兩列數(shù)據(jù)確定飛機(jī)飛行時(shí)的理論運(yùn)行軌跡和飛機(jī)的飛行速度；l 根據(jù)A、B、C三列數(shù)據(jù)確定當(dāng)飛機(jī)到達(dá)M點(diǎn)時(shí)觀測(cè)站O到M點(diǎn)的距離。編寫(xiě)腳本文件plane_rate.mplot(t,s

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

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

13、theta_theory作為該二階伺服系統(tǒng)的輸入時(shí)，系統(tǒng)的輸出盡可能的接近E列所給出的theta_output。求出系統(tǒng)傳遞函數(shù)參數(shù)。用simulink設(shè)計(jì)標(biāo)準(zhǔn)二階伺服系統(tǒng)，其中w的值待定設(shè)計(jì)的框圖如下：再編寫(xiě)腳本文件，如下：for n=1:20;%取不同的w值時(shí)，求輸出的theta_output sim(erjiesifuxitong);%運(yùn)行simulink中的框圖文件， %獲得輸出的的p_theta_output值，為一個(gè)長(zhǎng)度為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中，再依次比較向量中 %的值，當(dāng)取最小值時(shí)說(shuō)明p_theta_output與理論值最接近，此時(shí) %的w值即為標(biāo)準(zhǔn)w值end當(dāng)w取120時(shí)，得到的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ù)可看出，當(dāng)w=11時(shí)，取得最小的差值，再取w=1012之間的數(shù)，求得其中最小差值。得到的simulink框圖為：3.將C列theta_theory作為設(shè)計(jì)好的標(biāo)準(zhǔn)二階伺服系統(tǒng)的輸入信號(hào)，計(jì)算輸出相對(duì)于輸入的相對(duì)誤差；作出輸出隨時(shí)間變化的曲線，以及相對(duì)誤差隨時(shí)間變化的曲線。作輸出變化曲線和相對(duì)誤差曲線求相對(duì)誤差的腳本文件為：sim(erjiesifuxitong);%運(yùn)行二階伺服系統(tǒng)for i=1:143001;abs_dif(i)=abs(p_theta_output(i)-theta_theory(

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

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

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

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

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

22、正弦形狀的干擾信號(hào)。為求出干擾信號(hào)的具體差值，我再對(duì)diff進(jìn)行一次線性擬合，編寫(xiě)函數(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),即對(duì)diff曲線進(jìn)行輸入為時(shí)間t的線性擬合，擬合結(jié)果如下：x =-0.049815057234210 0.100054669976175 -0.006617837055889resn

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

24、框圖plot(t,p_theta_output,r);%作仿真結(jié)果曲線，為紅色hold on;3) 畫(huà)出輸出相對(duì)于輸入theta_theory的相對(duì)誤差隨時(shí)間變化的曲線，以及輸出相對(duì)于D列數(shù)據(jù)(t)的觀測(cè)值theta_observatuon)的相對(duì)誤差曲線。輸出相對(duì)于輸入theta_theory的相對(duì)誤差隨時(shí)間變化的曲線：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 輸出相對(duì)于D列數(shù)據(jù)(t)的觀測(cè)值theta_observatuon)的相對(duì)誤差曲線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)由電機(jī)驅(qū)動(dòng)電路、電機(jī)及幕布三部分組成。通過(guò)電機(jī)驅(qū)動(dòng)電路給電機(jī)輸入電信號(hào)，電機(jī)軸轉(zhuǎn)動(dòng)帶動(dòng)幕布