1、第 5 章 用 MATLAB 进行控制系统频域分析一、基于 MATLAB 的线性系统的频域分析基本知识（）频率特性函数 。)(jG设线性系统传递函数为：nnmmassabbs110)(则频率特性函数为：nnnn mmmjjjjwG)()()()( 110 由下面的 MATLAB 语句可直接求出 G(jw)。i=sqrt(-1) % 求取-1 的平方根 GW=polyval(num，i*w)./polyval(den，i*w) 其中（num，den）为系统的传递函数模型。而 w 为频率点构成的向量，点右除（./）运算符表示操作元素点对点的运算。从数值运算的角度来看，上述算法在系统的极点附近精度不
2、会很理想，甚至出现无穷大值，运算结果是一系列复数返回到变量 GW 中。（）用 MATLAB 作奈魁斯特图。控制系统工具箱中提供了一个 MATLAB 函数 nyquist( )，该函数可以用来直接求解Nyquist 阵列或绘制奈氏图。当命令中不包含左端返回变量时，nyquist（）函数仅在屏幕上产生奈氏图，命令调用格式为：nyquist(num,den) nyquist(num,den,w) 或者nyquist(G) nyquist(G,w) 该命令将画出下列开环系统传递函数的奈氏曲线：)()(sdenumG如果用户给出频率向量 w,则 w 包含了要分析的以弧度/秒表示的诸频率点。在这些频率点上
3、，将对系统的频率响应进行计算，若没有指定的 w 向量，则该函数自动选择频率向量进行计算。w 包含了用户要分析的以弧度/秒表示的诸频率点,MATLAB 会自动计算这些点的频率响应。当命令中包含了左端的返回变量时，即:re,im,w=nyquist(G) 或re,im,w=nyquist(G,w) 函数运行后不在屏幕上产生图形，而是将计算结果返回到矩阵 re、im 和 w 中。矩阵re 和 im 分别表示频率响应的实部和虚部，它们都是由向量 w 中指定的频率点计算得到的。在运行结果中，w 数列的每一个值分别对应 re、im 数列的每一个值。例 5.1 考虑二阶典型环节:18.0)(2ssG试利用
4、MATLAB 画出奈氏图。利用下面的命令，可以得出系统的奈氏图，如图 5-1 所示。 num=0,0,1;den=1,0.8,1;nyquist(num,den)% 设置坐标显示范围v=-2,2,-2,2;axis(v)gridtitle(Nyquist Plot of G(s)=1/(s2+0.8s+1) )（3）用 MATLAB 作伯德图控制系统工具箱里提供的 bode()函数可以直接求取、绘制给定线性系统的伯德图。当命令不包含左端返回变量时，函数运行后会在屏幕上直接画出伯德图。如果命令表达式的左端含有返回变量，bode()函数计算出的幅值和相角将返回到相应的矩阵中，这时屏幕上不显示频率响
5、应图。命令的调用格式为：mag,phase,w=bode(num,den) mag,phase,w=bode(num,den,w) 或mag,phase,w=bode(G) 图 5-1 二阶环节奈氏图mag,phase,w=bode(G,w) 矩阵 mag、phase 包含系统频率响应的幅值和相角，这些幅值和相角是在用户指定的频率点上计算得到的。用户如果不指定频率 w,MATLAB 会自动产生 w 向量，并根据 w 向量上各点计算幅值和相角。这时的相角是以度来表示的，幅值为增益值，在画伯德图时要转换成分贝值，因为分贝是作幅频图时常用单位。可以由以下命令把幅值转变成分贝：magdb=20log1
6、0(mag) 绘图时的横坐标是以对数分度的。为了指定频率的范围，可采用以下命令格式：logspace(d1，d2) 或 logspace(d1，d2，n) 公式 logspace(d1，d2) 是在指定频率范围内按对数距离分成 50 等分的，即在两个十进制数 和 之间产生一个由 50 个点组成的分量，向量中的点数 5010d20d是一个默认值。例如要在 弧度/秒与 弧度/秒之间的频区画伯德图，则输1.102入命令时， , 在此频区自动按对数距离等分成 50 个频率点，)(log10d)(log02d返回到工作空间中，即w=logspace(-1，2)要对计算点数进行人工设定，则采用公式 log
7、space(d1，d2，n)。例如，要在与 之间产生 100 个对数等分点，可输入以下命令：102w=logspace(0，3，100)在画伯德图时，利用以上各式产生的频率向量 w，可以很方便地画出希望频率的伯德图。由于伯德图是半对数坐标图且幅频图和相频图要同时在一个绘图窗口中绘制，因此，要用到半对数坐标绘图函数和子图命令。1）对数坐标绘图函数利用工作空间中的向量 x，y 绘图，要调用 plot 函数，若要绘制对数或半对数坐标图，只需要用相应函数名取代 plot 即可，其余参数应用与 plot 完全一致。命令公式有：semilogx(x，y，s) 上式表示只对 x 轴进行对数变换，y 轴仍为线
8、性坐标。semilogy(x，y，s) 上式是 y 轴取对数变换的半对数坐标图。Loglog(x，y，s) 上式是全对数坐标图，即 x 轴和 y 轴均取对数变换。2）子图命令MATLAB 允许将一个图形窗口分成多个子窗口，分别显示多个图形，这就要用到subplot()函数，其调用格式为：subplot(m，n，k)该函数将把一个图形窗口分割成 mn 个子绘图区域，m 为行数，n 为列数，用户可以通过参数 k 调用各子绘图区域进行操作，子图区域编号为按行从左至右编号。对一个子图进行的图形设置不会影响到其它子图，而且允许各子图具有不同的坐标系。例如，subplot(4，3，6)则表示将窗口分割成
9、43 个部分。在第 6 部分上绘制图形。 MATLAB 最多允许 99 的分割。例 5.2 给定单位负反馈系统的开环传递函数为：)7(10)sG试画出伯德图。利用以下 MATLAB 程序，可以直接在屏幕上绘出伯德图如图 5-2。 num=10*1,1;den=1,7,0;bode(num,den)gridtitle(Bode Diagram of G(s)=10*(s+1)/s(s+7) )该程序绘图时的频率范围是自动确定的，从 0.01 弧度/秒到 30 弧度/秒，且幅值取分贝值， 轴取对数，图形分成 2 个子图，均是自动完成的。如果希望显示的频率范围窄一点，则程序修改为： num=10*1
10、,1;den=1,7,0;w=logspace(-1,2,50); % 从 0.1 至 100，取 50 个点。mag, phase, w=bode(num, den, w);magdB=20*log10(mag) % 增益值转化为分贝值。% 第一个图画伯德图幅频部分。subplot(2,1,1);semilogx(w,magdB, -r) % 用红线画gridtitle(Bode Diagram of G(s)= 10*(s+1)/s(s+7) )图 5-2 自动产生频率点画出的伯德图xlabel(Frequency(rad/s)ylabel(Gain(dB)% 第二个图画伯德图相频部分。s
11、ubplot(2,1,2); semilogx(w,phase, -r); gridxlabel(Frequency(rad/s)ylabel(Phase(deg) ) 修改程序后画出的伯德图如 5-3 所示。（4）. 用 MATLEB 求取稳定裕量同前面介绍的求时域响应性能指标类似，由 MATLAB 里 bode()函数绘制的伯德图也可以采用游动鼠标法求取系统的幅值裕量和相位裕量。例如，我们可以在图 20 的幅频曲线上按住鼠标左键游动鼠标，找出纵坐标（Magnitude）趋近于零的点，从提示框图中读出其频率约为 7.25dB。然后在相频曲线上用同样的方法找到横坐标（Frequence）最接近
12、 7.25dB的点，可读出其相角为-53.9 度，由此可得，此系统的相角裕量为 126.1 度。幅值裕量的计算方法与此类似。此外，控制系统工具箱中提供了 margin()函数来求取给定线性系统幅值裕量和相位裕量，该函数可以由下面格式来调用：Gm, Pm, Wcg, Wcp=margin(G); (67)可以看出，幅值裕量与相位裕量可以由 LTI 对象 G 求出，返回的变量对（Gm, Wcg）为幅值裕量的值与相应的相角穿越频率，而（Pm, Wcp）则为相位裕量的值与相应的幅值穿越频率。若得出的裕量为无穷大，则其值为 Inf，这时相应的频率值为 NaN(表示非数值)，Inf 和 NaN 均为 MA
13、TLAB 软件保留的常数。图 5-3 用户指定的频率点画出的伯德图如果已知系统的频率响应数据，我们还可以由下面的格式调用此函数。Gm, Pm, Wcg, Wcp=margin(mag, phase, w);其中(mag, phase, w)分别为频率响应的幅值、相位与频率向量。例 5.3 已知三阶系统开环传递函数为： )23(27)23ssG利用下面的 MATLAB 程序，画出系统的奈氏图，求出相应的幅值裕量和相位裕量，并求出闭环单位阶跃响应曲线。 G=tf(3.5,1,2,3,2);subplot(1,2,1);% 第一个图为奈氏图nyquist(G);gridxlabel(Real Axi
14、s)ylabel(Imag Axis)% 第二个图为时域响应图Gm,Pm,Wcg,Wcp=margin(G)G_c=feedback(G,1);subplot(1,2,2);step(G_c)gridxlabel(Time(secs) )ylabel(Amplitude) 显示结果为：ans=1.1429 1.15781.7321 1.6542图 5-4 三阶系统的奈氏图和阶跃响应图画出的图形如图 5-4 所示。由奈氏曲线可以看出，奈氏曲线并不包围（-1，j0）点，故闭环系统是稳定的。由于幅值裕量虽然大于 1，但很接近 1，故奈氏曲线与实轴的交点离临界点（-1，j0）很近，且相位裕量也只有 7
15、.1578o，所以系统尽管稳定，但其性能不会太好。观察闭环阶跃响应图，可以看到波形有较强的振荡。如果系统的相角裕量 45 o,我们一般称该系统有较好的相角裕量。例 5.4 考虑一个新的系统模型，开环传递函数为： )9)(150(2ssG由下面 MATLAB 语句可直接求出系统的幅值裕量和相位裕量: G=tf(100*conv(1,5,1,5), conv(1,1,1,1,9);Gm, Pm, Wcg, Wcp=margin(G)结果显示 Gm = Pm =Inf 85.4365Wcg = Wcp =NaN 100.3285再输入命令 G_c=feedback(G,1); step(G_c)gr
16、idxlalel(Time(sec) )ylalel(Amplitude)可以看出，该系统有无穷大幅值裕量，且相角裕量高达 85.4365o。所以系统的闭环响应是较理想的，闭环响应图如图 5-5.（5）.时间延迟系统的频域响应(1) 时间延迟系统的传递函数模型带有延迟环节 e-Ts的系统不具有有理函数的标准形式，在 MATLAB 中，建立这类系统的图 5-5 较理想的系统响应模型。要由一个属性设置函数 set()来实现。该函数的调用格式为：set(H, 属性名, 属性值) 其中 H 为图形元素的句柄（handle） 。在 MATLAB 中，当对图形元素作进一步操作时，只需对该句柄进行操作即可。
17、例如以下调用格式h=plot(x,y)G=tf(num,den)Plot()函数将返回一个句柄 h，tf()函数返回一个句柄 G，要想改变句柄 h 所对应曲线的颜色，则可以调用下面命令：Set(h,color,1,0,0);即对“color”参数进行赋值，将曲线变成红色（由1,0,0决定）同样，要想对 G 句柄所对应模型的延迟时间T d进行修改，则可调用下面命令Set(G, T d,T)其中 T 为延迟时间。由此修改后，模型 G 就已具有时间延迟特性。(2) 时间延迟系统的频域响应含有一个延迟环节的系统，其开环频域响应为TjjTee)()(可见，该系统的幅频特性不变，只加大了相位滞后。例 5.
18、5 考虑系统的开环模型为:TsesG1)(当 T=1 时，我们可以由下面的 MATLAB 命令绘出系统的奈氏图，如图 5-6 所示,此系统对应的时域响应图为 5-7。 G=tf(1,1, 1); T=1;w=0, logspace(-3, 1, 100), logspace(1,2,200);set(G,Td, T); % 延迟 1 秒。nyquist(G,w)gridfigure % 建立一个新的绘图窗口step(G)图 5-6 时间延迟系统奈氏图 图 5-7 时间延迟系统的阶跃响应例 5.6 控制系统如图 5-8 所示，试分析其闭环特性。图 5-8 下面是该系统的演示过程，带符号语句为注释
19、语句。%清文字屏%clc%清图形屏%clf %系统 1 描述a1=1 3 0;b1=2;%系统 2 描述a2=1 2 2;b2=1 2;%两系统串联后的传递函数：b12,a12=series(b1,a1,b2,a2); %串联后频率响应（波特图）figure(1)bode(b12,a12);pause%clf%求振幅及相位裕度figure(2)m,p,w=bode(b12,a12)%求开环振幅 Gm 和相位裕度 Pm 及其对应频率margin(m,p,w);Gm,Pm,Wg,Wp=margin(m,p,w);pause%画 nichols 曲线，先用 ngrid 函数画坐标网figure(3)
20、ngrid(new);nichols(b12,a12);pause%clf%画根轨迹，先消除画 nichols 曲线时冻结的坐标系，%用 axis 命令恢复自动定标。figure(4)axis;rlocus(b12,a12);%用鼠标器找适当的极点k1,p1= rlocfind(b12,a12);pause%将环路闭合b1=k1*b1;b121=k1*b12;%再看闭环振幅及相位裕度figure(5)m,p,w=bode(b121,a12);margin(m,p,w);Gm,Pm,Wg,Wp=margin(m,p,w);pause%求闭环传递函数， bb，ab指单位反馈情况，% bb1，ab1
21、指将系统 2 置于反馈通道上：bb,ab=cloop(b121,a12);pausebb1,ab1=feedback(b1,a1,b2,a2,-1);pause%求闭环频率响应figure(6)bode(bb1,ab1);pause%clg;%闭环过渡过程figure(7)step(bb1,ab1);pause%clgfigure(8)impulse(bb1,ab1);pause%clg;%再求闭环零一极点，作为校核：zb,pb,kb=tf2zp(bb,ab);pause printsys(zb,pb,kb);formatecho off%演示结束运行结果：开环系统幅值裕量及相位裕量：Gm=3
22、.5435,Pm=55.022,Wg=1.6159,Wp=0.66898 ) ) 0 0 0 1 2 0 ) ) 0 ) 0 1 2 0 t ) ) 0 0 0 0 0 0 B B B B B B B B B B s 0 1 023在根轨迹上选择一点：selected_point = -4.0158 + 0.0309i此时闭环系统幅值裕量及相位裕量：Gm= 0.34677,Pm= -29.162,Wg= 1.6159,Wp= 2.5192 0 0 1 2 0 ) 0 0 0 0 0 1 2İ 5 0 5 e 0 1 2 3 4 5 65 0 05 0 1 2 3 4 5 60 0 0 0 根据
23、选定根轨迹上一点的位置可计算出相应的闭环零点，闭环极点及闭环增益值。二、附录 基于 MATLAB 的线性系统的频域分析的其它实例例 5.7 已知单位负反馈系统前向通道的传递函数为： s510s5s282)(G346234试绘制出 Bode 图并计算系统的频域性能指标。程序及结果如下：num=0 0 2 8 12 8 2;den=1 5 10 10 5 1 0;sys=tf(num,den);mag,phase,w=bode(sys);gm,pm,wcp,wcg=margin(mag,phase,w);margin(mag,phase,w)kg=20*log(gm)结果：gm=332.17, p
24、m=38.692, wg=25.847, wc=1.249kg=50.4272 ) 0 1 例 5.8 已知一带延迟因子的系统开环传函为： ，0.8s-2e)3s(51)G试求其有理传递函数的频率响应，同时在同一坐标中绘制以 Pade 近似延迟因子式系统的 Bode 图，并求此时系统的频域性能指标。程序及结果如下：n1=10 50;d1=conv(1 3,1 3);s1=tf(n1,d1);w=logspace(-1,2);mag1,phase1=bode(n1,d1,w);np,dp=pade(0.8,4);%s1=tf(n1,d1)*tf(np,dp);numt=conv(n1,np);d
25、ent=conv(d1,dp);s2=tf(numt,dent);mag2,phase2=bode(numt,dent,w);subplot(211);semilogx(w,20*log10(mag1),w,20*log10(mag2),r-);grid ontitle(bode plot);xlabel(Frequency（radsec）);ylabel(Gain db);subplot(212);grid semilogx(w,phase1,w,phase2,r-);grid onxlabel(Frequency（radsec）);ylabel(Phase deg);gm1,pm1,wcp
26、1,wcg1=margin(s1);gm2,pm2,wcp2,wcg2=margin(s2);例 5.9 已知知两个单位负反馈系统开环传递函数分别为：，s45s7.2)(G31s45s7.2)(G32试用 Bode 图法判断闭环系统的稳定性。程序如下：num1=0 0 0 2.7;den1=1 5 4 0;sys1=tf(num1,den1);figure(1);hold ongm,pm,wcp,wcg=margin(sys1);margin(sys1);title(对数频率特性图 1);xlabel(频率 rad/sec);ylabel(Gain dB);num2=0 0 0 2.7;den2=1 5 -4 0;sys2=tf(num2,den2);figure(2);gm,pm,wcp,wcg=margin(sys2);margin(sys2);title(对数频率特性图 2);xlabel(频率 rad/sec);ylabel(Gain dB);结论：系统 1 稳定，系统 2 不稳定（Warning: The closed-loop system is unstable.） 0 0 1 2 c 0 0 1 2 5 0
