最新用MATLAB实现线性系统的频域分析

实验二用MATLAB实现线性系统的频域分析

[实验目的]

1．掌握MATLAB平台下绘制典型环节及系统开环传递函数的Bode图和Nyquist图（极坐标图）绘制方法；

2．掌握利用Bode图和Nyquist图对系统性能进行分析的理论和方法。

[实验指导]

一、绘制Bode图和Nyquist图

1．Bode图绘制

采用bode()函数，调用格式：

①bode(sys)；bode(num,den)；

系统自动地选择一个合适的频率范围。

②bode(sys，w)；

其中w(即ω)是需要人工给出频率范围，一般由语句w=logspace(a,b,n)给出。logspace(a,b,n)：表示在10a到10b之间的 n个点,得到对数等分的w值。

③bode(sys,{wmin,wmax})；

其中{wmin,wmax}是在命令中直接给定的频率w的区间。

以上这两种格式可直接画出规范化的图形。

④[mag,phase,ω]=bode(sys)或[m,p]=bode(sys)

这种格式只计算Bode图的幅值向量和相位向量，不画出图形。

m为频率特性G(jω )的幅值向量;

p为频率特性G(jω )的幅角向量，单位为角度（°）。

w为频率向量，单位为[弧度]/秒。

在此基础上再画图，可用：

subplot(211);semilogx(w,20*log10(m) %对数幅频曲线

subplot(212);semilogx(w,p) %对数相频曲线

⑤bode(sys1,sys2,…,sysN) ;

⑥bode((sys1,sys2,…,sysN，w);

这两种格式可在一个图形窗口同时绘多个系统的bode图。

2. Nyquist曲线的绘制

采用nyquist()函数调用格式：

① nyquist(sys) ;

② nyquist(sys,w) ;

其中频率范围w由语句w=w1:Δw:w2确定。

③nyquist(sys1,sys2,…,sysN) ;

④nyquist(sys1,sys2,…,sysN,w);

⑤ [re,im,w]=nyquist(sys) ;

re—频率响应实部im—频率响应虚部

使用命令axis()改变坐标显示范围，例如axis([-1,1.5,-2,2])。

⑥当传递函数串有积分环节时ω=0处会出现幅频特性为无穷大的情况，可用命令axis()，自定义图形显示范围，避开无穷大点。

二、系统分析

1．计算控制系统的稳定裕度

采用margin( )函数可以直接求出系统的幅值裕度和相角裕度。

调用格式为：

①[Gm,Pm,Wcg,Wcp]= margin(num,den) ;

[Gm,Pm,Wcg,Wcp]= margin(A,B,C,D) ;

[Gm,Pm,Wcg,Wcp]= margin(sys) ;

Gm--- 幅值裕度；

Pm---相位裕度；

wcg ---幅值裕度处对应的频率ωc；

wcp ---相位裕度处对应的频率ωg。

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

③ margin(sys)

在当前图形窗口中绘制出系统裕度的Bode图。

2．用幅值裕度和相角裕度判断闭环系统稳定性与相对稳定性

3．用Nyquist图判断闭环系统稳定性

由Nyquist曲线包围（-1，j0）点的情况，根据Nyquist稳定判据判断闭环系统稳定性。

三、举例

例1：振荡环节如下：16

()21016

G s s s =

++，做出该环节的Bode 图和Nyquist 图。

程序：

>>n=[16];d=[1 10 16];sys=tf(n,d);figure(1);bode(sys);figure(2);nyquist(sys)

运行结果：

Bode Diagram

Frequency (rad/sec)

P h a s e (d e g )

M a g n i t u d e (d B )

10

10

10

10

Nyquist Diagram

Real Axis

I m a g i n a r y A x i

s

例2：振荡环节如下：16()216

n G s s s ξω=

++，做出该环节的Bode 图和Nyquist 图。

ξ变化，取[0.05,0.1,0.2,0.5,0.7,1,2]。

1．Bode 图程序：

>> wn=8;znb=[0.05,0.1,0.5,0.7,2];w=logspace(0,2,10000);figure(1);n=[wn^2]; for k=znb d=[1 2*k*wn wn^2];sys=tf(n,d);bode(sys,w);hold on; end 运行结果：

Bode Diagram

Frequency (rad/sec)

P h a s e (d e g )

M a g n i t u d e (d B

)-60-40-200204010

10

10

-180

-135-90-450

1．Nyquist 图程序：

>> wn=8;znb=[0.05,0.1,0.5,0.7,2];w=logspace(0,2,10000);figure(1);n=[wn^2]; for k=znb d=[1 2*k*wn wn^2];sys=tf(n,d); nyqiust(sys,w); hold on; end 运行结果：

Nyquist Diagram

Real Axis

I m a g i n a r y A x i s

例3：系统开环传递函数如下：

①020

()(0.51)

G s s s =+，

②0.231

()0.0551

c s G s s +=+，

③20(0.231)

()(0.0551)(0.51)

s G s s s s +=

++，

做出各自的Bode 图，并求①、③幅值裕度和相角裕度 1．Bode 图程序：

>> n1=20;d1=conv([1,0],[0.5,1]);sys1=tf(n1,d1); figure(2);bode(sys1); n2=[0.23 1];d2=[0.055,1];sys2=tf(n2,d2);hold on; figure(2);bode(sys2);

n=[4.6 20];d=conv([1,0],conv([0.055,1],[0.5,1]));sys=tf(n,d);hold on; figure(2); bode(sys)

运行结果：

Bode Diagram

Frequency (rad/sec)

P h a s e (d e g )

M a g n i t u d e (d B )

10

10

10

10

2．求①②幅值裕度和相角裕度程序（图形与数据） >> n1=20;d1=conv([1,0],[0.5,1]);sys1=tf(n1,d1); figure(1);margin(sys1)

运行结果：

Bode Diagram

Frequency (rad/sec)

P h a s e (d e g )

M a g n i t u d e (d B

)-1001020304050Gm = Inf, P m = 17.964 deg (at 6.1685 rad/sec)

10

10

10

-180

-135

-90

>> n=[4.6 20];d=conv([1,0],conv([0.055,1],[0.5,1]));sys=tf(n,d);hold on; figure(2); margin(sys)

运行结果：

Bode Diagram

Frequency (rad/sec)

P h a s e (d e g )

M a g n i t u d e (d B )

Gm = Inf, P m = 50.472 deg (at 8.9542 rad/sec)

10

10

10

10

2．求①②幅值裕度和相角裕度程序和结果（数据）

>> [Gm,Pm,Wcg,Wcp]= margin(sys1)

Gm = Inf Pm = 17.9642 Wcg = Inf

Wcp = 6.1685

>> [Gm,Pm,Wcg,Wcp]= margin(sys)

Gm = Inf Pm = 50.4719 Wcg = Inf Wcp = 8.9542

例4：系统开环传递函数为：0431)

()(21)

s G s s s +=+（ 做出nyquist 图，按nyquist 稳定判据判

断闭环系统的稳定性。

程序与结果：

>> n=conv([4],[3 1]);d=conv([1 0],[2 1]);sys2=tf(n,d) Transfer function: 12 s + 4 --------- 2 s^2 + s

>> figure(4);nyquist(sys2);v=[-1,6,-60,60];axis(v)

Nyquist Diagram

Real Axis

I m a g i n a r y A x i s

分析判断：p=0，nyquist 曲线没有包围(-1，j0)点，闭环系统是稳定的。 下面通过闭环系统时域阶跃响应来验证闭环系统的稳定性：

>> n=conv([4],[3 1]);d=conv([1 0],[2 1]);G1=tf(n,d); G2=1;G=feedback(G1,G2,-1) Transfer function: 12 s + 4 ---------------- 2 s^2 + 13 s + 4 >> figure(7);step(G)

Time (sec)

A m p l i t u d e

0.10.20.30.40.50.6

0.70.80.91

ω=0-