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

西安電子科技大學(xué)

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

指導(dǎo)老師：屈勝利

計(jì)算機(jī)輔助系統(tǒng)分析:

G0(s)

__________K(l+s)

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

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

要求：

1.當(dāng)GN=1，K=40時(shí)，用MATLAB畫出開環(huán)Bode圖，求出。,、%。由其

估計(jì)出/,、小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)

所得參數(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ì)算乙、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)

結(jié)果：rise_time=0.7772

max_overshoot=0.8168

settling_time=17.0985

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

2.當(dāng)K=40時(shí)，用MATLAB畫出根軌跡圖，并求出K=40時(shí)的閉環(huán)極點(diǎn)；由

其估計(jì)出仆、3%

繪制系統(tǒng)的根軌跡:

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。

系統(tǒng)閉環(huán)零極點(diǎn)：

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

[p,z]=pzmap(sys)

極點(diǎn)：p=-40.1616

-0.2274+2.4146i

-0.2274-2.4146i

-1.0837

零點(diǎn)：z=-1

求小ts,3%：根據(jù)零極點(diǎn)的特點(diǎn)，取主導(dǎo)極點(diǎn)p=?0.2274+2.4146i和p=-0.2274-2.4146i,

可以將此系統(tǒng)化為二階系統(tǒng)個(gè)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)

結(jié)果：rise_time=0.7772

max__overshoot=0.7438

settling_time=16.9691

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

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

答：建立系統(tǒng)模型

①用自適應(yīng)變步長(zhǎng)方法(ode45)仿真可得:

2

1.8

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

15202530

觀察上圖可知：

超調(diào)量約為80%、上升時(shí)間Tr約為0.5S；調(diào)整時(shí)間Ts約為18S。

②定步長(zhǎng)RK-2法

H=0.08,仿真結(jié)果發(fā)散

6

5

4

3

2

1

0

H=0.02,收斂

綜上可知：

仿真的步長(zhǎng)必須選擇合理才能更準(zhǔn)確的得出仿真結(jié)果，步長(zhǎng)太大，截?cái)嗾`差

大，甚至導(dǎo)致仿真失敗，而步長(zhǎng)太小，使舍入誤差逐漸積累，也導(dǎo)致誤差增大。

所以步長(zhǎng)的合理選擇和仿真精度有很大的關(guān)系。

4.令圖1的K=40

①G”分別為：

當(dāng)％為飽和特性時(shí)，仿真結(jié)果如下:

1.6

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

當(dāng)與為死區(qū)特性時(shí)，仿真結(jié)果如下:

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

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

，仿真之，并計(jì)算乙、ts>3%

仿真圖形：

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

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

比較自適應(yīng)變步長(zhǎng)法和定步長(zhǎng)RK-2法的仿真結(jié)果得出，RK-2的q較長(zhǎng)，但兩

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

原因：

自適應(yīng)變步長(zhǎng)法中步長(zhǎng)腦的大小與y的變化率有關(guān)，當(dāng)系統(tǒng)趨于穩(wěn)定時(shí)，y

變化較慢，步長(zhǎng)&變大，而定步長(zhǎng)時(shí)屏為一個(gè)較小值，故定步長(zhǎng)RK-2法的務(wù)較大。

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

不大。

若非線性環(huán)節(jié)加在GO(s)之前，死區(qū)特性僅作用于誤差e(s)上，若非線性環(huán)

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

較大，所以后者對(duì)原系統(tǒng)的動(dòng)態(tài)性能影響較大，從而今、J、。%與原系統(tǒng)相差

較多；

二、病態(tài)系統(tǒng)(stiff)仿真(simulink)

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

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

r:單位階躍：

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

1.用自適應(yīng)變步長(zhǎng)法(RK45)仿真之

解:當(dāng)H時(shí)。仿真結(jié)果如下:

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

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

解:當(dāng)勺=1。2.&=10-2時(shí)。

收斂的步長(zhǎng)h的范圍是：h<=0.027o當(dāng)h=0.027時(shí)的仿真結(jié)果如下:

2

0.5

-0.5

4860488049004920496049805000

Timeoffset:0

3.用病態(tài)仿真算法仿真之

以上三問，均打印出仿真曲線，計(jì)算暫態(tài)響應(yīng)，并比較討論之。

答：

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

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

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

9999

當(dāng)t=10s時(shí)，可算出y(t)=0.39,用自適應(yīng)變步長(zhǎng)法(RK45)仿真得到的y(t)=0.164,

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

y(t)=0.370o

顯然，結(jié)果都比較接近，但病態(tài)仿真算法仿真得到的結(jié)果最準(zhǔn)確，定步長(zhǎng)四階龍

格庫塔法最不準(zhǔn)確。

原因：

被仿真系統(tǒng)的遼=10000>>50,系統(tǒng)嚴(yán)重病態(tài)，而RK4-5和RK4若取低步長(zhǎng)則仿真

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

三.計(jì)算機(jī)輔助控制器設(shè)計(jì):

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

1.開關(guān)處于A時(shí)，系統(tǒng)性能滿足上述要求否？

仿真波形:

1.8

I：:：

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

I:::::

Q__I_________|____________|____________|____________|____________|____________

051015202530

Timeoffset:0

Bode圖為:

由仿真波形可以知道，開關(guān)接A時(shí)，不能滿足指標(biāo)。

2.開關(guān)處于B時(shí)，計(jì)算機(jī)輔助設(shè)計(jì)Gc(s)，使系統(tǒng)性能滿足上述要求

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

主程序：

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)%超前校正環(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);

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])；

程序運(yùn)行結(jié)果：

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(上升時(shí)間)

max_over=0.2235(超調(diào)量)

ts=1.2279s(調(diào)整時(shí)間)

校正前后階躍響應(yīng)曲線：

..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)

注釋：藍(lán)色線表示校正后，紅色線表示校正前

建立如圖系統(tǒng)驗(yàn)證指標(biāo)是否符合要求：

ToWcxkspace

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

滿足系統(tǒng)指標(biāo)要求。

四.觀測(cè)站（0點(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è)站。點(diǎn)的連線到直

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

其間隔為0.001S；B列為飛機(jī)飛行過程中t時(shí)刻飛機(jī)相對(duì)于P點(diǎn)的距離S（t）；C

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

觀測(cè)站實(shí)際上觀測(cè)到的飛機(jī)在t時(shí)刻時(shí)角度9（t）的觀測(cè)值theta_observation；E

列給出了當(dāng)把C列數(shù)據(jù)theta_theory作為某標(biāo)準(zhǔn)二階伺服系統(tǒng)G⑸的輸入信號(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.根據(jù)A、B、C三列數(shù)據(jù)確定當(dāng)飛機(jī)到達(dá)M點(diǎn)時(shí)觀測(cè)站。到M點(diǎn)的距離。

G6)二?——.

S-+2sli+<

(3)設(shè)計(jì)標(biāo)準(zhǔn)二階伺服系統(tǒng)，要求:

1.確定合適的。值，其單位階躍響應(yīng)的性能指標(biāo)滿足：

系統(tǒng)的超調(diào)量。％介于4.5%~8.0%之間

2.使用MATLAB仿真確定參數(shù)3n,使得：當(dāng)把C列數(shù)據(jù)theta_theory作為該二階

伺服系統(tǒng)的輸入時(shí)，系統(tǒng)的輸出盡HT能的接近E列所給出的theta_output。求出

系統(tǒng)傳遞函數(shù)參數(shù)。

3.將C列theta_theory作為設(shè)計(jì)好的標(biāo)準(zhǔn)二階伺服系統(tǒng)的輸入信號(hào)，計(jì)算輸出相

對(duì)于輸入的相扁誤差；作出輸出隨時(shí)間變化的曲線，以及相對(duì)誤差隨時(shí)間變化的

曲幺戈0

（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ù)

（0（t）的觀測(cè)值theta_observatuon）。試：

1.確定干擾信號(hào)和東線性環(huán)節(jié)的相關(guān)參數(shù)，并畫出系統(tǒng)simulink模型。

2.畫出系統(tǒng)的輸出隨時(shí)間變化的曲線。

3.畫出輸出相對(duì)于輸Atheta_theory的相對(duì)誤差隨時(shí)間變化的曲線，以及輸出

相對(duì)于D列數(shù)據(jù)（。⑴的觀測(cè)值theta_observatuon）的相對(duì)誤差曲線。

附表一：

觀察站測(cè)得角度

觀察站測(cè)得角度標(biāo)準(zhǔn)二階伺服系

（觀測(cè)值

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

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

一、完成內(nèi)容

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

用。

編寫腳本readxlsx.m,讀取各列數(shù)據(jù)的值：

%將附表?中的數(shù)據(jù)導(dǎo)入matlab工作空間，使各列數(shù)據(jù)都能作為變量使用

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

%讀取表格文件的信息

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

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

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

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

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

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

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

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

數(shù)據(jù)

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

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

formatlong;

%設(shè)定顯示格式為long,因?yàn)楸砀裰袛?shù)據(jù)有9位有效值

運(yùn)行結(jié)果：依次可取得各數(shù)據(jù)的值

?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)試根據(jù)表格中的數(shù)據(jù)使用MATLAB完成以下問題

?根據(jù)A、B兩列數(shù)據(jù)確定飛機(jī)飛行時(shí)的理論運(yùn)行軌跡和飛機(jī)的飛行速度;

?根據(jù)A、B、C三列數(shù)據(jù)確定當(dāng)飛機(jī)到達(dá)M點(diǎn)時(shí)觀測(cè)站0到M點(diǎn)的距離。

編寫腳本文件plane_rate.m

plot(t,s);gridon;

%作飛機(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/s

i=l:143001;

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

%取觀察站測(cè)得角度的絕對(duì)值

[min_theta_theoryk]=min(abs_theta_theory);

%當(dāng)飛機(jī)到達(dá)M點(diǎn)時(shí)，觀察站測(cè)得角度應(yīng)為最小角，故求得此最小角度時(shí)

為對(duì)應(yīng)的時(shí)間t(k)

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

%當(dāng)時(shí)間為t(k)時(shí),對(duì)應(yīng)的位移為s(k),tan(theta)=PM/OM

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

運(yùn)行結(jié)果：

?planerate

distance_OM=

5.000030000000002e+02

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

2

(3)設(shè)計(jì)標(biāo)準(zhǔn)二階伺服系統(tǒng)G(s)=-——%-----------

S2+2gs+說

解：1.確定合適的C值，其單位階躍響應(yīng)的性能指標(biāo)滿足：

系統(tǒng)的超調(diào)量。％介于4.5%~8.0%之間

—

l_z2

因?yàn)镺%=ed1*100%,使O%介于4.5%~8.0%之間則可取。=0.68,此時(shí)

。％=%5.4280251339488,滿足要求。

2.使用MATLAB仿真確定參數(shù)3n,使得：當(dāng)把C列數(shù)據(jù)theta_theory作為

該二階伺服系統(tǒng)的輸入時(shí)，系統(tǒng)的輸出盡可能的接近E列所給出的

theta_outputo求出系統(tǒng)傳遞函數(shù)參數(shù)。

用simulink設(shè)計(jì)標(biāo)準(zhǔn)二階伺服系統(tǒng)，其中w的值待定設(shè)計(jì)的框圖如下：

再編寫腳本文件，如下：

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

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

%獲得輸出的的P_theta_oulput值，為一個(gè)長(zhǎng)度為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中,再依次比較向量中

%的值，當(dāng)取最小值時(shí)說明p_theta_output與理論值最接近，此時(shí)

研枷值即為標(biāo)準(zhǔn)w值

end

當(dāng)w取1?20時(shí)，得到的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

由數(shù)據(jù)可看出，當(dāng)W=ll時(shí)，取得最小的差值,再取w=10~12之間的數(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)

fori=l:143001;

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

為求輸出pthetaoutput的絕對(duì)誤差

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

%求輸出的相對(duì)誤差

End

作輸出變化曲線為：

從圖中可見系統(tǒng)輸出的曲線與給出的理論曲線兒乎重合，可見選擇的系統(tǒng)參數(shù)比

較準(zhǔn)確。

當(dāng)作相對(duì)誤差曲線時(shí)，發(fā)現(xiàn)為如下圖所示：

接近于沖激函數(shù)。因?yàn)闀r(shí)間軸選用的過寬，無法看出細(xì)節(jié)，故用axis([0,2,0,2])指令截取相應(yīng)

的曲線，就能看清楚了：

a)時(shí)間軸為0?2時(shí)

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

b)時(shí)間軸為2?70時(shí),

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

時(shí)間軸為時(shí),

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

40

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

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

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

對(duì)誤差的求取公式為產(chǎn)警造鬻呻卜1。。路而theta_theory的值在70s?73s時(shí)近似取到

了零，而theta_theory在分母匕故此時(shí)相對(duì)誤差值會(huì)變得極大。當(dāng)theta_theory取值遠(yuǎn)離零

值時(shí)，由于絕對(duì)誤差接近于0,故相對(duì)誤差也接近于0。而在0s?2s區(qū)間可看到一個(gè)明顯的

振蕩過程，分析原因：因?yàn)橄到y(tǒng)為二階系統(tǒng)，輸入的起始階段響應(yīng)會(huì)有超調(diào)量，但當(dāng)經(jīng)過一

段時(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ù)(0⑴的觀測(cè)值

theta_observation)o試：

1)確定干擾信號(hào)和非線性環(huán)節(jié)的相關(guān)參數(shù)，并畫出系統(tǒng)simulink模型。

■作出theta_theory的時(shí)間曲線和theta_observation的時(shí)間曲線如圖,藍(lán)色為

theta_theory的時(shí)間曲線，紅色為theta_observation的時(shí)間曲線。

編寫M文件，使每選定一種非線性環(huán)節(jié)或添加一種干擾都能得到對(duì)應(yīng)的曲線，再將對(duì)應(yīng)曲

線與給定的觀察曲線作比較，即能找到對(duì)應(yīng)的參數(shù)。

Thetaobservation.m

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

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

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

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

holdon;

首先在前向通道中加一個(gè)死區(qū)非線性環(huán)節(jié)，得到如下圖所示波形：

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

再將死區(qū)非線性放到回路外：

111*2

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

I36Tl

From

DeadZoneIntegratorTransferFenToWorkspace

Workspace

得至I」曲線:

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

線打印出來如圖：

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

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

?為找到干擾環(huán)節(jié)，我首先將觀察曲線與輸出曲線作對(duì)比，如下圖:

為求出合適的干擾信號(hào)使藍(lán)色曲線接近紅色曲線，先將兩曲線作差

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

由圖可看出，差值diff的時(shí)間曲線近似為一個(gè)正弦波形，故應(yīng)輸入一個(gè)正弦形狀的干擾

信號(hào)。

為求出干擾信號(hào)的具體差值，我再對(duì)diff進(jìn)行一次線性擬合，編寫函數(shù)文件如下：

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),即對(duì)diff曲線進(jìn)行

輸入為時(shí)間t的線性擬合，擬合結(jié)果如下：

x=

-0.0498150572342100.100054669976175-0.006617837055889

resnorm=

1.820876780757486

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

得到圖為:

由圖可知，擬合結(jié)果近似接近于原曲線，因?yàn)樵€幅值較小，故將此擬合結(jié)果作為干擾信

號(hào)輸入。輸入系統(tǒng)框圖如圖：

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

輸出此框圖的結(jié)果為:

由圖可看出，兩條曲線近似重合，故此系統(tǒng)即為標(biāo)準(zhǔn)的輸入系統(tǒng)。

2)畫出系統(tǒng)的輸出隨時(shí)間變化的曲線。

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

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

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

holdon;

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

對(duì)于D列數(shù)據(jù)(。(t)的觀測(cè)值theta_observatuon)的相對(duì)誤差曲線。

輸出相對(duì)于輸入theta_theory的相對(duì)誤差隨時(shí)間變化的曲線：

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]);

?輸出相對(duì)于D列數(shù)據(jù)(0(t)的觀測(cè)值theta_observatuon)的相對(duì)誤差曲線

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]);

五.幕布系統(tǒng)由電機(jī)驅(qū)動(dòng)電路、電機(jī)及幕布三部分組成。通過電機(jī)驅(qū)動(dòng)電路給

電機(jī)輸入電信號(hào)，電機(jī)軸轉(zhuǎn)動(dòng)帶動(dòng)幕布上升或下降。通常狀況下，長(zhǎng)按“升”按

鈕會(huì)使幕布順利升起；但有時(shí)長(zhǎng)按“升”按鈕會(huì)在底端出現(xiàn)“卡死”（振蕩）現(xiàn)

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

試用Matlab分別仿真這三個(gè)過程。

（提示：這三個(gè)過程可能會(huì)受非線性、二階諧振等環(huán)節(jié)的影響）

注：交作業(yè)時(shí)，寫明題目、所用的仿真方法、圖形及計(jì)算結(jié)果。

一.簡(jiǎn)化模型。展開的