西電系統(tǒng)仿真大作業(yè)

1、-系統(tǒng)仿真大作業(yè)指導(dǎo)老師:屈勝利班 級(jí):學(xué) 號(hào):姓 名: 系統(tǒng)仿真上機(jī)作業(yè)一、 計(jì)算機(jī)輔助系統(tǒng)分析:系統(tǒng)如下圖所示: 其中 r 是單位階躍, G N 是非線性器件,G 0 (s = 1025. 0(1625. 0(110(s 1+s s s s K ( 1.當(dāng) G N = 1、 K = 40 時(shí),用 MATLAB 畫出開環(huán) Bode 圖,求出 c 、 B 。由其估計(jì)出 t r 、 t s 、 % 。解:matlab 程序如下: 運(yùn)行結(jié)果如下: 由上圖知:den = 0.1563,6.5156 ,10.6500 ,1.0000 ,0增益裕量 Gm = 4.3168 相位裕量 Pm = 10.01

2、58 穿越頻率 Wcg = 5.1598 增益為 0的頻率 Wcp =2.3975所以得到:c =5.1598 B =10.0158仿真波形如下: 求解 t r 、 t s 、 %:由 bode 圖:可將系統(tǒng)的開環(huán)傳遞函數(shù)化成 G(s=(11040+s s , 其相應(yīng)的閉環(huán)傳遞函數(shù)為 G B (s=10s2+s+40,此二階系統(tǒng)可求出相應(yīng)的 =0.025, w n =2。 所以可以估計(jì)出: 18.67smatlab 程序: 運(yùn)行結(jié)果為: 所以得到:t r =0.7772s t s =17.0985s % =81.68%經(jīng)檢驗(yàn):由 bode 圖估算出 tr, % ,ts 的結(jié)果與正確值的差距不大

3、 .2.當(dāng) G N = 1, K =40時(shí), 用 MATLAB 畫出根軌跡圖, 并求出 K =40時(shí)的閉環(huán)極點(diǎn); 由其估計(jì)出 t r 、 t s 、 % 。解:matlab 程序?yàn)? 運(yùn)行結(jié)果: 求系統(tǒng)閉環(huán)零極點(diǎn) matlab 程序: 運(yùn)行結(jié)果: 仿真圖: 由上圖:閉環(huán)系統(tǒng)的所有極點(diǎn)為 -40, -0.2274+ 2.4146j,-0.2274-2.4146j , -1.083,閉環(huán)系統(tǒng)的零點(diǎn)為 -1。根據(jù)零極點(diǎn)對(duì)消 原理,可將 -1.083與 -1對(duì)消掉,取主導(dǎo)極點(diǎn),傳遞函數(shù)為 G(s=4146. 22274. 0(4146. 22274. 0(s40j s j -+ =882. 54548

4、. 02s 5.882+s 由上述 G(s可求得: tr=0.69,ts=15.39s,%=74.38%Matlab 程序如下 : 運(yùn)行結(jié)果: (r t (%=74.38% (s t 與估算結(jié)果差距不大。3.當(dāng) GN = 1, K = 40 時(shí),仿真之,并由仿真結(jié)果求出 tr 、 ts 、 % 。用自適應(yīng)變步長(zhǎng)法用定步長(zhǎng) RK-2 法, 并觀察步長(zhǎng)與仿真結(jié)果的收斂于發(fā)散關(guān)系, 以及仿真精度。 解:G 0(s =1025. 0(1625. 0(110(s 1+s s s s K ( K=40時(shí): =40(.61(10. (s 1256+s s s s (采用自適應(yīng)變步長(zhǎng)法: simulink 仿

5、真模塊連接如下: 仿真波形: 計(jì)算求得: r t =0.69,s t =17.14s,%=81.72%Matlab 程序:ytr=find(simout(:,1>=1;rise_time=tout(ytr(1ymax,tp=max(y;peak_time=tout(tp;max_overshoot=ymax-1s=length(tout;while y(s>0.98&y(s<1.02s=s-1;endsettling_time=tout(s運(yùn)行結(jié)果:rise_time = 0.7630s (rt max_overshoot =0.8169 (%=81.69%settl

6、ing_time =17.1290 s (st 與估算結(jié)果相差不大。用定步長(zhǎng) RK-2 法:h 取 0.01: r t =0.69s, s t =17.14s, % =81.72%h 取 0.02: r t =0.70s, s t =17.14s, % =81.69%h 取 0.03: r t =0.69s, s t =17.13s, % =81.73% h 取 0.04: r t =0.72s, s t =17.16s, %=81.74% h 取 0.05: r t =0.71s, s t =17.15s, %=81.70%h 取 0.052: 綜上可知:使用 RK 2定步長(zhǎng)法時(shí),當(dāng) h 取

7、 0.05以上時(shí),將會(huì)出現(xiàn) 發(fā)散,因此,增大步長(zhǎng)會(huì)使系統(tǒng)穩(wěn)定性下降。在 h 小于等于 0.05時(shí),對(duì)應(yīng)的系統(tǒng)的性能如下所示:4.令圖中的 K = 40 。 G N 分別為:。分別仿真之,并由仿真結(jié)果求出 t r 、 t s 、 % 。解 :1 加飽和非線性特性曲線:simulink仿真模塊連接如下: 仿真波形如下:計(jì)算求得:rt =0.91,st =16.36s,%=72.76%2加死區(qū)非線性特性曲線:simulink 仿真模塊連接如下: 仿真波形: 計(jì)算求得:r t =0.7424,s t =29.35s,%=80.67% G N = 1,在 G 0(s 之后,反饋點(diǎn)之前加 。 仿真之,并由

8、仿真結(jié)果求出 t r 、 t s 、 % 。解:simulink 仿真模塊連接如下:仿真波形: 計(jì)算如下:r t =0.7270s, s t =17.27s, %=89.42%。對(duì) 3 和 4 中、 的 tr 、 ts 、 % 比較, 并解釋差異的原因。 答:3中仿真方法不同,從而精度不同最后導(dǎo)致性能差異。 原因:自適應(yīng)變步長(zhǎng)法中步長(zhǎng) hk 的大小與 y 的變化率有關(guān),當(dāng)系統(tǒng)趨于穩(wěn) 定時(shí), y 變化較慢, hk 變大; hk 為較小值時(shí),定步長(zhǎng) RK-2法的 ts 較大。而在仿真前期, y 變化較劇烈,兩種方法的 hk 都比較小,故 其 tr ,超調(diào)量相差不大。4中由下表知:中加入死區(qū)非線性與

9、飽和非線性相比,加入飽 和非線性的系統(tǒng)性能更好, 的指標(biāo)小于, 中延遲大, 二者不同 說明在系統(tǒng)不同位置加入相同的非線性環(huán)節(jié),對(duì)系統(tǒng)的影響不同。 二、病態(tài)系統(tǒng)(stiff 仿真(simulink 1.用自適應(yīng)變步長(zhǎng)法(RK45仿真之解:選取 1 =1000, 2=0.001(1 simulink 連接圖如下: 仿真波形: 2.用定步長(zhǎng)四階龍格庫(kù)塔法仿真,并試著搜索收斂的步長(zhǎng) h 的范 圍;若找不到 h ,將 1 減小, 2 增大,用定步長(zhǎng)四階龍格庫(kù)塔 法仿真,尋找 h。(2解:simulink 設(shè)置為定步長(zhǎng) RK-4仿真法仿真波形:h=0.279 h=0.28 h=0.3 結(jié)論:由上圖分析可知

10、,當(dāng) h<=0.279時(shí),系統(tǒng)仿真收斂,當(dāng) h>0.279時(shí) , 仿真將發(fā)散。3. 用病態(tài)仿真算法仿真之。 解 : 設(shè)置 simulink 為病態(tài)仿真 分析比較: 當(dāng)設(shè)置 1=1000, 2=0.001時(shí): 系統(tǒng)輸出 1G( (101(0.11s s s =+ 則求出 y(t=0.1101-10099199t t e e -+由上面最后一張圖 :當(dāng) t=6.2時(shí),求出輸出 y(t=0.3985用自適應(yīng)變步長(zhǎng)法 RK45仿真得到的 y(t=0.3994用定步長(zhǎng) RK-4(h=0.279仿真得到的 y(t=0.3981用病態(tài)仿真算法仿真得到的 y(t=0.3987通過比較發(fā)現(xiàn)用病態(tài)仿真

11、算法得到的暫態(tài)響應(yīng)與理論計(jì)算值最接近, 而定步長(zhǎng) RK-4其次,但步長(zhǎng)過大時(shí),系統(tǒng)穩(wěn)定性降低,自適應(yīng)變步 長(zhǎng) 法 精 度 最 低 。 并 且 該 仿 真 系 統(tǒng) 本 身 就 是 一 個(gè) 病 態(tài) 系 統(tǒng) m a x m i n10050s =>,所以用適應(yīng)它的病態(tài)仿真算法最精確。三.計(jì)算機(jī)輔助控制器設(shè)計(jì): 要求:開環(huán) B 45o , c 4.2 ,且 t r 0.4s , t s 1.5s , % 25%。 1. 開關(guān)處于 A 時(shí),系統(tǒng)性能滿足上述要求否?解:(1開關(guān)出于 A 時(shí), simulink 連接圖如下: 仿真結(jié)果: 系統(tǒng) bode 圖:Matlab 程序?yàn)? 程序運(yùn)行結(jié)果: bo

12、de 圖為: 求得:增益裕量 gm=400.07 ,相位裕量:pm =17.9930穿越頻率 wcg=200.0277 ,截止頻率(0dB wcp =3.0837 求 tr,ts, % 的 matlab 源程序: 運(yùn)行結(jié)果: 上升時(shí)間 tr =0.6447s , 峰值時(shí)間 tp =1.0541s%=0.6045s ,調(diào)整時(shí)間 ts =7.2762s所以不符合要求。2. 開關(guān)處于 B 時(shí),計(jì)算機(jī)輔助設(shè)計(jì) G c (s ,使系統(tǒng)性能滿足上述 要求。解 :將開關(guān)處于 B ,參照自動(dòng)控控原理(劉丁著,需采用頻率 法超前矯正,使系統(tǒng)滿足指標(biāo)要求Matlab 程序?yàn)?num=10;den=1 1 0;G=

13、tf(num,den;kc=1.05;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(1step(Gy_c,'r' ,10; %畫校正前系統(tǒng)階躍相應(yīng)曲線 hold onstep(Gx_c,'b-.' ,10 %畫校正后系統(tǒng)階躍相應(yīng)曲線 figure(2bode(G,'r' %畫校正前系統(tǒng)波特圖 hold onbode(G*Gc*kc,'b' %畫校正前系統(tǒng)波

14、特圖 子函數(shù):function Gc=plsj(G,kc,yPmG=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;End運(yùn)行結(jié)果得超前矯正環(huán)節(jié)為 :校正后系統(tǒng) simulink 模塊連接圖為:

15、仿真波形為: 較粗的那條為校正后的波形,其性能指標(biāo):r t =0.2604s, s t =0.3712s, % =22.9534% 滿足系統(tǒng)指標(biāo)要求。四 .觀測(cè)站 (O 點(diǎn) 為測(cè)得的某航班數(shù)據(jù),當(dāng)飛機(jī)到達(dá)某位置 P 時(shí)開始對(duì)飛機(jī)的相關(guān)數(shù)據(jù)進(jìn)行記錄,在位置 P 處時(shí)間 t 記為 0,當(dāng)飛機(jī)到觀測(cè)站的距離達(dá)到最短時(shí)飛機(jī)所處的位置 記為 M 點(diǎn),飛機(jī)在 t 時(shí)刻所處的位置與觀測(cè)站 O 點(diǎn)的連 線到直線 OM 的夾角記為 (t。附表一給出了觀測(cè)站某次 記錄數(shù)據(jù):A 列為記錄時(shí)間 t,其間隔為 0.001s; B 列為 飛機(jī)飛行過程中 t 時(shí)刻飛機(jī)相對(duì)于 P 點(diǎn)的距離 S(t; C 列給出了飛機(jī) t 時(shí)

16、刻角度 (t的理論值 theta_theory (參考輸入, D 列給出了觀測(cè)站實(shí)際上觀測(cè)到的飛機(jī)在 t 時(shí)刻時(shí)角度 (t的觀測(cè)值 theta_observation; E 列給出 了當(dāng)把 C 列數(shù)據(jù) theta_theory 作為某標(biāo)準(zhǔn)二階伺服系統(tǒng) G(s的 輸 入 信 號(hào) 時(shí) 該 標(biāo) 準(zhǔn) 二 階 伺 服 系 統(tǒng) 的 輸 出 值 theta_output。作業(yè)要求:(1 將附表一中的數(shù)據(jù)導(dǎo)入 matlab 工作空間, 使各列數(shù)據(jù)都能作 為變量使用。解:(2 試根據(jù)表格中的數(shù)據(jù)使用 MATLAB 完成以下問題 1. 根據(jù) A、 B 兩列數(shù)據(jù)確定飛機(jī)飛行時(shí)的理論運(yùn)行軌跡和飛機(jī) 的飛行速度;2.

17、根據(jù) A、 B 、 C 三列數(shù)據(jù)確定當(dāng)飛機(jī)到達(dá) M 點(diǎn)時(shí)觀測(cè)站 O 到 M 點(diǎn)的距離。解:1.matlab 程序 ; 飛行曲線: 求解觀測(cè)站 O 點(diǎn)到 M 點(diǎn)的距離。Matlab 代碼如下: 運(yùn)行結(jié)果:先擬合飛機(jī)位移與實(shí)踐的關(guān)系,求解出飛機(jī)的飛行速度為 140m/s, 在根據(jù)三角形 OMP 求解出 OM 變得長(zhǎng)度。 其中 data1(71430 對(duì)應(yīng)飛機(jī)飛行到 M 的時(shí)刻。求解得 OM 的長(zhǎng)度為 500.0030m(3 設(shè)計(jì)標(biāo)準(zhǔn)二階伺服系統(tǒng) G(s=n 2, 要求: S + 2n +n1.確定合適的 值,其單位階躍響應(yīng)的性能指標(biāo)滿足 :系統(tǒng)的超調(diào)量 %介于 4.5%8.0%之間解:matlab

18、 程序:dete=;for yipsin=0:0.1:1dt=100*exp(-yipsin*pi/(sqrt(1-yipsin2; dete=dete,dtendyipsin=0:0.1:1;plot(yipsin,dete,'-b'grid onxlabel('取值 'ylabel('取值 %'仿真波形: 由上圖可確定取 為 0.65時(shí),可滿足 %介于 4.5%8.0%之間。 驗(yàn)證,當(dāng) =0.65=6.8%時(shí), ,滿足要求2. 使用 MATLAB 仿真確定參數(shù) n, 使得:當(dāng)把 C 列數(shù) theta_theory 作為該二階伺服系統(tǒng)的輸入時(shí),

19、系統(tǒng)的輸出盡可能的接近 E 列所給 出的 theta_output。求出系統(tǒng)傳遞函數(shù)參數(shù)。 解:simulink 仿真模塊連接為:仿真波形:當(dāng) wn=1時(shí): Wn=2時(shí) Wn=5由上圖:wn=5時(shí),兩條曲線基本吻合,只是在起始階段存在誤差, 故 wn 大于 5時(shí),只需觀察兩條曲線在起始階段的差異即可。Wn=7 Wn=9 Wn=10 Wn=12由上圖可知,當(dāng) wn 大于 9時(shí),兩條曲線的差異開始增大,故當(dāng) wn 取 9時(shí),系統(tǒng)的輸出更接近 theta_output。=0.65, wn=9傳遞函數(shù)為:G(s=n / S 2 + 2n + n 2=9/(s2+11.7+813. 將 C 列 theta_theory 作為設(shè)計(jì)好的標(biāo)準(zhǔn)二階伺服系統(tǒng)的輸入 信號(hào),計(jì)算輸出相對(duì)于輸入的相對(duì)誤差;作出輸出隨時(shí)間變化的曲 線,以及相對(duì)誤差隨時(shí)間變化的曲線。解:用 matlab 作圖:。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。