系统仿真上机大作业

西安电子科技大学

系统仿真上机作业

指导老师：屈胜利

计算机辅助系统分析:

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分别仿真这三个过程。

（提示：这三个过程可能会受非线性、二阶谐振等环节的影响）

注：交作业时，写明题目、所用的仿真方法、图形及计算结果。

一.简化模型。展开的