西安電子科技大學

系統仿真上機作業

指導老師：屈勝利

計算機輔助系統分析:

G0(s)

__________K(l+s)

其中，為單位階躍,G,\,為非線性器件Go（s）

5(105+1)(0.625.?+1)(0.0255+1)

要求：

1.當GN=1，K=40時，用MATLAB畫出開環Bode圖，求出。,、%。由其

估計出/,、小8%

BlFigure1

FileEditViewInsertToolsDesktopWindowHelp

白昌U）41-、踮⑨要□園■口

1

00

50

pm

)50

800

p50

n9-0

u35

u■810

625

70

e■1

s-

0

9

1

P-

)-

1

-2

-

2

10-3

Frequencyrad/s(rad/s)

clear

num=[O004040];

denl=conv([l0],[101]);

den2=conv([0.6251],[0.0251]);

den=conv(denl,den2)

bode(num,den)

gridon;

xlabeU'Frequencyrad/s'/fontsize',10)

title(fBodeDiagramG(s)=40(1+s)/s(10s+1)(0.625s+1)(0.025s+1)r)

[mag,phase,w]=bode(num,den)；

(Gm,Pm,Wcg,Wcp]=margin(mag,phase,w)

所得參數為：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

計算乙、4、6%程序

sys=tf(num,den)

sys=feedback(sys,1);

[yzt]=step(sys);

ytr=find(y>=l);

rise_time=t(ytr(1))

[ymax,tp]=max(y);

peak_time=t(tp);

max_overshoot=ymax-l

s=length(t);

whiley(s)>0.98&y(s)<1.02

s=s-l;

end

settling_time=t(s)

結果：rise_time=0.7772

max_overshoot=0.8168

settling_time=17.0985

可知：rr=0.7772s4=17.0985s8%=81.68%

2.當K=40時，用MATLAB畫出根軌跡圖，并求出K=40時的閉環極點；由

其估計出仆、3%

繪制系統的根軌跡:

num=[0004040];

denl=conv([10]z[101]);

den2=conv([0.6251],[0.0251]);

den=conv(deni,den2);

sys=tf(num,den);

[p,z]=pzmap(sys)

rlocus(sys)

根軌跡圖形：

(

-S

?

s

p

u

o

o

a

?o

x

v

A

J

e

u

a

BqQ

E-

-

1o。

系統閉環零極點：

程序：sys=feedback(sys,1);

[p,z]=pzmap(sys)

極點：p=-40.1616

-0.2274+2.4146i

-0.2274-2.4146i

-1.0837

零點：z=-1

求小ts,3%：根據零極點的特點，取主導極點p=?0.2274+2.4146i和p=-0.2274-2.4146i,

可以將此系統化為二階系統個G(s)=sys=

5.882

sA2+0.4548s+5.882,有以下程序可求出4、8%

程序：

num=[005.882];

den=[l0.45480];

sys=tf(num,den)

sys=feedback(sys,1);

[y,t]=step(sys);

ytr=find(y>=l);

rise_time=t(ytr(1))

[ymax,tp]=max(y);

peak_time=t(tp);

max_overshoot=ymax-l

s=length(t);

whiley(s)>0.98&y(s)<1.02

s=sT;

end

settling_time=t(s)

結果：rise_time=0.7772

max__overshoot=0.7438

settling_time=16.9691

可知：tr=0.7772s16.969Is5%=74.38%

3.GN=1，K=40：仿真之，并由仿真結果求出,、4、3%

答：建立系統模型

①用自適應變步長方法(ode45)仿真可得:

2

1.8

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

15202530

觀察上圖可知：

超調量約為80%、上升時間Tr約為0.5S；調整時間Ts約為18S。

②定步長RK-2法

H=0.08,仿真結果發散

6

5

4

3

2

1

0

H=0.02,收斂

綜上可知：

仿真的步長必須選擇合理才能更準確的得出仿真結果，步長太大，截斷誤差

大，甚至導致仿真失敗，而步長太小，使舍入誤差逐漸積累，也導致誤差增大。

所以步長的合理選擇和仿真精度有很大的關系。

4.令圖1的K=40

①G”分別為：

當％為飽和特性時，仿真結果如下:

1.6

tr=l.37s=15.7s＞。%=58%o

當與為死區特性時，仿真結果如下:

*二0.49s-ts=34.2s>o%=100%o

②GN=1，在G°(s)之后，反饋點之前加上

，仿真之，并計算乙、ts>3%

仿真圖形：

J=0.43s,t$=173s、o%=72%O

③對3和4中①、②的“、.”、。％比較，并解釋差異的原因。

比較自適應變步長法和定步長RK-2法的仿真結果得出，RK-2的q較長，但兩

種方法的%、。%相差較小。

原因：

自適應變步長法中步長腦的大小與y的變化率有關，當系統趨于穩定時，y

變化較慢，步長&變大，而定步長時屏為一個較小值，故定步長RK-2法的務較大。

而在仿真前期,y變化較劇烈，兩種方法的步長甌都比較小，故其J、。外相差

不大。

若非線性環節加在GO(s)之前，死區特性僅作用于誤差e(s)上，若非線性環

節加在G0(s)之后，則死區特性作用于e(s)G。(s)上，在該系統中，q(S)增益

較大，所以后者對原系統的動態性能影響較大，從而今、J、。%與原系統相差

較多；

二、病態系統(stiff)仿真(simulink)

r----------y

-------?G(s)--------?

r:單位階躍：

G(s)=-a--o--F7----+-----D------a-o-.--s---+-----D-------

1.用自適應變步長法(RK45)仿真之

解:當H時。仿真結果如下:

2.用定步長四階龍格庫塔法仿真，并試著搜索收斂的步長h的范圍；若找不到h,

將H增大，憶減小，用定步長四階龍格庫塔法仿真，尋找h。

解:當勺=1。2.&=10-2時。

收斂的步長h的范圍是：h<=0.027o當h=0.027時的仿真結果如下:

2

0.5

-0.5

4860488049004920496049805000

Timeoffset:0

3.用病態仿真算法仿真之

以上三問，均打印出仿真曲線，計算暫態響應，并比較討論之。

答：

Y(s)=G(s)R(s)=-

s(ios-nj(ai?-n5aa-tia

y(t)=------s3+一‘3+lo

9999

當t=10s時，可算出y(t)=0.39,用自適應變步長法(RK45)仿真得到的y(t)=0.164,

用定步長四階龍格庫塔法仿真得到的y(t)<0.2,用病態仿真算法仿真得到的

y(t)=0.370o

顯然，結果都比較接近，但病態仿真算法仿真得到的結果最準確，定步長四階龍

格庫塔法最不準確。

原因：

被仿真系統的遼=10000>>50,系統嚴重病態，而RK4-5和RK4若取低步長則仿真

速度極慢，若不如此則誤差較大。

三.計算機輔助控制器設計:

要求：開環％245°,小?4.2,且，，W0.4s,ts<l.5s,8%<25%

1.開關處于A時，系統性能滿足上述要求否？

仿真波形:

1.8

I：:：

0.2-?.....................................................................................:............................................-1

I:::::

Q__I_________|____________|____________|____________|____________|____________

051015202530

Timeoffset:0

Bode圖為:

由仿真波形可以知道，開關接A時，不能滿足指標。

2.開關處于B時，計算機輔助設計Gc(s)，使系統性能滿足上述要求

(1)基于頻率法的串聯超前校正

主程序：

function[]=question3_2(k)

num=k*[1];

den=conv([1,0],[1,1]);

G=tf(num,den);

kc=l;

yPm=45+10;

Gc=plsj(G,kczyPm)%超前校正環節

Gy_c=feedback(G,1)%校正前系統閉環傳遞函數

Gx_c=feedback(G*kc*Gc,1)%校正后系統閉環傳遞函數

figure(1)

step(Gy_c,'r',5);

holdon

step(Gx_c,'b',5);

gridon

figure(2)

1

bode(Gz*r)

holdon

bode(G*kc*Gcz*b*)

gridon

[Gm,Pm,Wcg,Wcp]=margin(G*kc*Gc);

r=Pm

Wc=Wcp

[y,t]=step(Gx_c);

ytr=find(y>=l);

rise_time=t(ytr(1));

tr=rise_time

[ymaxztp]=max(y);

peak_time=t(tp);

max_overshoot=ymax-l;

max__over=max_overshoot

s=length(t);

whiley(s)>0.98&y(s)<1.02

s=s-l;

end

settling_time=t(s+1);

ts=settling_time

子程序：

functionGc=plsj(G,kc,yPm)

G=tf(G);

[mag,pha,w]=bode(G*kc);

Mag=20*logl0(mag);

[Gm,Pm.WcgzWcp]=margin(G*kc);

1

phi=(yPm-getfield(Pmz*Wcg))*pi/180;

alpha=(1+sin(phi))/(1-sin(phi));

Mn=-10*logl0(alpha);

Wcgn=spline(Mag,w,Mn);

T=l/(Wcgn*sqrt(alpha));

Tz=alpha*T;

Gc=tf([Tz,l],[T,l])；

程序運行結果：

Transferfunction:

0.4536s+1

0.1126s+1

Transferfunction:

10

sA2+s+10

Transferfunction:

4.536s+10

0.1126sA3+1.113sA2+5.536s+10

R=49.7706

Wc=4.4248

tr=0.3982s(上升時間)

max_over=0.2235(超調量)

ts=1.2279s(調整時間)

校正前后階躍響應曲線：

..68

8

np

_=

mds

v5o.

6

0.4-葉/-?”

0.2

0*^------L_

00.5

校正前后bode圖:

BodeDiagram

100

m50

)p

c

no

u

u

罩

a

9

p

)36

9一

s

B

o

d

Frequency(rad/sec)

注釋：藍色線表示校正后，紅色線表示校正前

建立如圖系統驗證指標是否符合要求：

ToWcxkspace

可得：tr=0.2604s,f,=0.3712s,3%=22.9534%

滿足系統指標要求。

四.觀測站（0點）為測得的某航班數據，當飛機到達某位置P時開始對飛機的相

關數據進行記錄，在位置P處時間t記為0,當飛機到觀測站的距離達到最短時

飛機所處的位置記為M點，飛機在t時刻所處的位置與觀測站。點的連線到直

線OM的夾角記為0（t）。附表一給出了觀測站某次記錄數據：A列為記錄時間t,

其間隔為0.001S；B列為飛機飛行過程中t時刻飛機相對于P點的距離S（t）；C

列給出了飛機t時刻角度o（t）的理論值theta_theory（參考輸入），D列給出了

觀測站實際上觀測到的飛機在t時刻時角度9（t）的觀測值theta_observation；E

列給出了當把C列數據theta_theory作為某標準二階伺服系統G⑸的輸入信號時

該標準二階伺服系統的輸出值theta_output。

作業要求：

（1）將附表一中的數據導入matlab工作空間，使各列數據都能作為變量使用。

（2）試根據表格中的數據使用MATLAB完成以下問題

1.根據A、B兩列數據確定飛機飛行時的理論運行軌跡和飛機的飛行速度；

2.根據A、B、C三列數據確定當飛機到達M點時觀測站。到M點的距離。

G6)二?——.

S-+2sli+<

(3)設計標準二階伺服系統，要求:

1.確定合適的。值，其單位階躍響應的性能指標滿足：

系統的超調量。％介于4.5%~8.0%之間

2.使用MATLAB仿真確定參數3n,使得：當把C列數據theta_theory作為該二階

伺服系統的輸入時，系統的輸出盡HT能的接近E列所給出的theta_output。求出

系統傳遞函數參數。

3.將C列theta_theory作為設計好的標準二階伺服系統的輸入信號，計算輸出相

對于輸入的相扁誤差；作出輸出隨時間變化的曲線，以及相對誤差隨時間變化的

曲幺戈0

（4）在上述設計好的伺服系統中加入合適的干擾信號和適當的非線性環節，使

得當把C列數據theta_theory作為輸入信號時系統的輸出盡可能接近D列數據

（0（t）的觀測值theta_observatuon）。試：

1.確定干擾信號和東線性環節的相關參數，并畫出系統simulink模型。

2.畫出系統的輸出隨時間變化的曲線。

3.畫出輸出相對于輸Atheta_theory的相對誤差隨時間變化的曲線，以及輸出

相對于D列數據（。⑴的觀測值theta_observatuon）的相對誤差曲線。

附表一：

觀察站測得角度

觀察站測得角度標準二階伺服系

（觀測值

時間t(/s)位移s(/m)（理論值統輸出

theta_observat

theta_theory）（theta_output）

ion）

00-1.52083793100

0.0010.14-1.520837233-7.07579E-05-7.5738E-05

0.0020.28-1.520836535-0.000281898-0.000301739

0.0030.42-1.520835836-0.000631723-0.000676187

nnnanRA-1A9HP2R12Q-nnm11MdG-nnm1Q797R

一、完成內容

(1)將附表一中的數據導入matlab工作空間，使各列數據都能作為變量使

用。

編寫腳本readxlsx.m,讀取各列數據的值：

%將附表?中的數據導入matlab工作空間，使各列數據都能作為變量使用

[typ,desc,fmt]=xlsfinfoC5finaltest.xlsx?)

%讀取表格文件的信息

t=xlsread(，finaltest,xlsx',1,'A2:A143002');

%讀取excel文件中第一個工作表中A列的全部數據，用數組t()保存數據

s=xlsreacK,finaltest,xlsx',1,'B2:B143002,);

%讀取excel文件中第一個工作表中B列的全部數據，用數組s()保存數據

theta_theory=xlsread(，finaltest,xlsx',1,'C2:C143002,);

%讀取excel文件中第一個工作表中C列的全部數據，用數組theta_thoory()保存數據

theta_observation=xlsread(，finaltest,xlsx',1,'D2:D143002,);

%讀取excel文件中第一個工作表中D列的全部數據，用數組thetaobsorvation()保存

數據

theta_output=xlsreadCfinaltest,xlsx',1,'E2:E143002,);

%讀取excel文件中第一個工作表中E列的全部數據，用數組theta_output()保存數據

formatlong;

%設定顯示格式為long,因為表格中數據有9位有效值

運行結果：依次可取得各數據的值

?t(2)

ans=

1.000000000000000e-03

?s(2)

ans=

0.140000000000000

?theta_theory(2)

ans=

-1.520837232808838

?theta_observation(2)

ans=

-7.075789556117257e-05

?theta_output(2)

ans=

-7.573799703386787e-05

(2)試根據表格中的數據使用MATLAB完成以下問題

?根據A、B兩列數據確定飛機飛行時的理論運行軌跡和飛機的飛行速度;

?根據A、B、C三列數據確定當飛機到達M點時觀測站0到M點的距離。

編寫腳本文件plane_rate.m

plot(t,s);gridon;

%作飛機位移關于時間t的圖象，可看出飛機是勻速運動

i=2:143001;

rate(i)=s(i)./t(i);

為求運行軌跡中各點的速率，可看出運動速率恒為140m/s

i=l:143001;

abs_theta_theory(i)=abs(theta_theory(i));

%取觀察站測得角度的絕對值

[min_theta_theoryk]=min(abs_theta_theory);

%當飛機到達M點時，觀察站測得角度應為最小角，故求得此最小角度時

為對應的時間t(k)

distanceOM=s(k)*cot(absthetatheory(1))

%當時間為t(k)時,對應的位移為s(k),tan(theta)=PM/OM

%其中PM=s(k),theta=thetatheory的最初的值，由此可求得0M的距離

運行結果：

?planerate

distance_OM=

5.000030000000002e+02

故飛機飛行速率為140m/sQM間距離為500m。

2

(3)設計標準二階伺服系統G(s)=-——%-----------

S2+2gs+說

解：1.確定合適的C值，其單位階躍響應的性能指標滿足：

系統的超調量。％介于4.5%~8.0%之間

—

l_z2

因為O%=ed1*100%,使O%介于4.5%~8.0%之間則可取。=0.68,此時

。％=%5.4280251339488,滿足要求。

2.使用MATLAB仿真確定參數3n,使得：當把C列數據theta_theory作為

該二階伺服系統的輸入時，系統的輸出盡可能的接近E列所給出的

theta_outputo求出系統傳遞函數參數。

用simulink設計標準二階伺服系統，其中w的值待定設計的框圖如下：

再編寫腳本文件，如下：

forn=l:20;%取不同的w值時,求輸出的thetaoutput

sim('erjiesifuxitong，);%運行simulink中的框圖文件,

%獲得輸出的的P_theta_oulput值，為一個長度為143001的列向量

forj=l:143001;

d_value(j)=abs(p_theta_output(j)-theta_output(j));

%求每組輸出的p_theta_output與理論值的差

end

abs_dif(n)=sum(d_value);%取誤差和,存到向量ads_dif中,再依次比較向量中

%的值，當取最小值時說明p_theta_output與理論值最接近，此時

研枷值即為標準w值

end

當w取1?20時，得到的abs_dif各值為：

?abs_dif

abs_dif=

1.0e+03*

Columns1through4

6.0194005109116212.7395408824277221.6440888630738101.100263457041429

Columns5through8

0.7732417253371810.5490055141110910.3840511989651720.256977609508468

Columns9through12

0.1568628001557720.0791096121171100.0363440280427000.064234872892925

Columns13through16

0.1099810646631350.1495906571967200.1837311444177030.213275684364825

Columns17through20

0.2389980608554720.2615401303097590.2814228978999480.299065677805618

由數據可看出，當W=ll時，取得最小的差值,再取w=10~12之間的數，求得其中最小差

值。

得到的simulink框圖為:

3.將C列theta_theory作為設計好的標準二階伺服系統的輸入信號，計算輸

出相對于輸入的相對誤差；作出輸出隨時間變化的曲線，以及相對誤差隨時間

變化的曲線。

作輸出變化曲線和相對誤差曲線

求相對誤差的腳本文件為：

sim(,erjiesifuxitong,);%運行二階伺服系統

fori=l:143001;

abs_dif(i)=abs(p_theta_output(i)-theta_theory(i));

為求輸出pthetaoutput的絕對誤差

relative_error(i)=abs_dif(i)/abs(theta_theory(i));

%求輸出的相對誤差

End

作輸出變化曲線為：

從圖中可見系統輸出的曲線與給出的理論曲線兒乎重合，可見選擇的系統參數比

較準確。

當作相對誤差曲線時，發現為如下圖所示：

接近于沖激函數。因為時間軸選用的過寬，無法看出細節，故用axis([0,2,0,2])指令截取相應

的曲線，就能看清楚了：

a)時間軸為0?2時

plot(t,relative_error);axis([0,2,0,1])

b)時間軸為2?70時,

圖形幾乎無變化，即相對誤差近似為0.

時間軸為時,

c)70?73plot(t,relative_error);axis([70,73,0,40])

40

由圖可看出相對誤差變化圖在70?73之間有一個猛增的過程

d)時間軸為73?143時，plot(t,relativeerror);axis([73,143,0,1])

分析圖可看出，在70s?73s之間，相對誤差有一個顯著的沖激式上升過程，分析原因：相

對誤差的求取公式為產警造鬻呻卜1。。路而theta_theory的值在70s?73s時近似取到

了零，而theta_theory在分母匕故此時相對誤差值會變得極大。當theta_theory取值遠離零

值時，由于絕對誤差接近于0,故相對誤差也接近于0。而在0s?2s區間可看到一個明顯的

振蕩過程，分析原因：因為系統為二階系統，輸入的起始階段響應會有超調量，但當經過一

段時間系統穩定后，就不會再發生振蕩。

4.在上述設計好的伺服系統中加入合適的干擾信號和適當的非線性環節，使得當把C列

數據theta_theory作為輸入信號時系統的輸出盡可能接近D列數據(0⑴的觀測值

theta_observation)o試：

1)確定干擾信號和非線性環節的相關參數，并畫出系統simulink模型。

■作出theta_theory的時間曲線和theta_observation的時間曲線如圖,藍色為

theta_theory的時間曲線，紅色為theta_observation的時間曲線。

編寫M文件，使每選定一種非線性環節或添加一種干擾都能得到對應的曲線，再將對應曲

線與給定的觀察曲線作比較，即能找到對應的參數。

Thetaobservation.m

%plot(t,theta_theory);holdon;%作理論曲線

plot(t,theta_observation,'r);holdon;%作觀察所得曲線

sim(，erjiesifuxitong，);%運行仿真框圖

plot(t,p_theta_output)；%作仿真結果曲線

holdon;

首先在前向通道中加一個死區非線性環節，得到如下圖所示波形：

此時輸出曲線與觀察曲線相差較大，故不正確。

再將死區非線性放到回路外：

111*2

:theta_theor\------->>_theta_outpi

I36Tl

From

DeadZoneIntegratorTransferFenToWorkspace

Workspace

得至I」曲線:

輸出曲線在中間停滯區與觀察曲線近似，為得到死區寬度，將題設所給的理論曲線與觀察曲

線打印出來如圖：

從上圖中可看出，死區寬度為±1,故將死區的起始與截止參數改為如下所示:

此時能得到與觀察曲線非線性彎曲部分重合的曲線:

?為找到干擾環節，我首先將觀察曲線與輸出曲線作對比，如下圖:

為求出合適的干擾信號使藍色曲線接近紅色曲線，先將兩曲線作差

(i=l:143001;diff(i)=p_theta_output(i)-theta_observation⑴;plot(t,diff);)得至妝口下曲線:

由圖可看出，差值diff的時間曲線近似為一個正弦波形，故應輸入一個正弦形狀的干擾

信號。

為求出干擾信號的具體差值，我再對diff進行一次線性擬合，編寫函數文件如下：

functionfun_sin=xianxingnihe(x,t)

fori=l:143001

fun_sin(i)=x(l)*sin(x(2)*t(i)+x(3)*pi)；

end

再輸入指令行

?x(0)=[0.10.10.1]

>>[x,resnorm]=lsqcurvefit(@xianxingnihe,xO,t,diff),即對diff曲線進行

輸入為時間t的線性擬合，擬合結果如下：

x=

-0.0498150572342100.100054669976175-0.006617837055889

resnorm=

1.820876780757486

再作出擬合所得曲線與真實的diff曲線作比較：plot(t,diff);holdon;pk)t(t,xianxingnihe(x,t),T);

得到圖為:

由圖可知，擬合結果近似接近于原曲線，因為原曲線幅值較小，故將此擬合結果作為干擾信

號輸入。輸入系統框圖如圖：

其中正弦波參數采用擬合曲線得到的參數，如下

輸出此框圖的結果為:

由圖可看出，兩條曲線近似重合，故此系統即為標準的輸入系統。

2)畫出系統的輸出隨時間變化的曲線。

plot(t,theta_observation);holdon;%作理論觀察曲線

sim(,erjiesifuxitong，);%運行仿真框圖

plot(t,p_theta_output,,r);與作仿真.結果曲線，為紅色

holdon;

3)畫出輸出相對于輸入theta.theory的相對誤差隨時間變化的曲線，以及輸出相

對于D列數據(。(t)的觀測值theta_observatuon)的相對誤差曲線。

輸出相對于輸入theta_theory的相對誤差隨時間變化的曲線：

i=l:143001;

diff_a(i)=abs(p_theta_output(i)-theta_theory(i))./abs(theta_theory(i));

plot(t,diff_a);holdon;

axis([014301]);

?輸出相對于D列數據(0(t)的觀測值theta_observatuon)的相對誤差曲線

i=l:143001;

diff_a(i)=abs(p_theta_output(i)-theta_observation(i))./abs(theta_observation

⑴)；

plot(t,diff_a);holdon;

axis([014300.1]);

五.幕布系統由電機驅動電路、電機及幕布三部分組成。通過電機驅動電路給

電機輸入電信號，電機軸轉動帶動幕布上升或下降。通常狀況下，長按“升”按

鈕會使幕布順利升起；但有時長按“升”按鈕會在底端出現“卡死”（振蕩）現

象；如果間歇性按下幕布開關，幕布就能越過“卡死”點，順利升起。

試用Matlab分別仿真這三個過程。

（提示：這三個過程可能會受非線性、二階諧振等環節的影響）

注：交作業時，寫明題目、所用的仿真方法、圖形及計算結果。

一.簡化模型。展開的