控制系统建模与仿真-基于MATLABSimulink的分析与实现（第2版）课件 第6讲MATLAB在频率及根轨迹分析中的应用

第6讲MATLAB在频率及根轨迹分析中的应用主讲内容一、控制系统的频域分析二、Matlab绘制Bode图三、绘制Nyquist图四、绘制Nichols图五、控制系统频域设计六、频域法校正设计七、根轨迹八、根轨迹校正设s=j,系统的频域响应可由下式直接求出：一、控制系统的频域分析频域分析法是应用频率特性研究控制系统的一种经典方法。采用这种方法可直观地表达出系统的频率特性，分析方法比较简单，物理概念比较明确。1、频域的表示方法2、计算系统的频率特性的方法（1）计算实频响应和虚频响应real(G(jω));%取频率特性的实部imag(G(jω));%取频率特性的虚部（2）计算幅频响应和相频响应20*log(abs(G(jω));%计算幅频值angle(G(jω));%计算相频值3、Matlab频域分析常用函数bode

画出系统伯特图(幅频/相频特性曲线)margin画曲线获得系统幅值和相位裕度

nyquist画出系统乃奎斯特曲线nichols画出频率图二、Matlab绘制Bode图[mag,phase,w]=bode(num,den,w)bode(num,den)bode(num,den,w)[mag,phase,w]=bode(sys)bode(sys1,sys2,...,w)1、bode图语法格式其中：mag表示幅值向量，注意要得到幅频用分贝表示dB,MAGDB=20*log10(mag)。phase表示相角向量，w角频率，它用对数表示,即：w=logspace(a,b,n)表示10a与10b之间的n点num=[0.55];d1=[0.510];d2=[1/25000.0121];den=conv(d1,d2);Bode(num,den)已知传递函数：例12、求增益和相角裕度语法格式[Gm,Pm,Wcg,Wcp]=margin(num,den)[Gm,Pm,Wcg,Wcp]=margin(sys)[Gm,Pm,Wcg,Wcp]=margin(mag,phase,w)margin(num,den)margin(sys)其中：Wcg=相角-180对应频率，Wcp=增益0dB对应频率(穿越频率）；

Gm=增益裕量(dB)，Pm=相角裕量(deg)num=[0.55];d1=[0.510];d2=[1/25000.0121];den=conv(d1,d2);margin(num,den)[Gm,Pm,Wcg,Wcp]=margin(num,den)从图上看出：已知传递函数：Gm（增益裕量）=28.6(dB)，Pm(相角裕量)=48.6(deg)Wcp=2.94(rad/s)例2三、

绘制Nyquist图[re,im,w]=nyquist(num,den,w)%绘制角频率w下的nyquist图[re,im]=nyquist(sys,w)nyquist(sys)nyquist(sys,w)nyquist(sys1,sys2,...,w)；其中：re表示实部向量，im表示虚部向量。2、绘制Nyquist图语法格式奈奎斯特图是用图解法表现系统频率特性的方法，将频率响应通过其幅频特性及相频特性表示在极坐标中的图形，横坐标表示实部，纵坐标为虚部的图，称为幅相图，或奈奎斯特(Nyquist)图。1、什么是Nyquist图num=[0.5];den=[1210.5];[re,im]=nyquist(num,den);subplot(1,2,1);plot(re,im),gridaxis([-1,1,-1.5,0.2])subplot(1,2,2);nyquist(num,den)axis([-1,1,-1.5,0.2])

%设置坐标轴的显示尺度已知传递函数：绘制系统的Nyquist图不包围(-1，j0)点，系统稳定例3num=[0.55];d1=[0.51];d2=[10.6/51];den=conv(d1,d2);nyquist(num,den)包围(-1，j0)点，系统不稳定已知传递函数：绘制系统的Nyquist图例4四、绘制Nichols图[mag,phase,w]=nichols(num,den,w)[mag,phase]=nichols(num,den)nichols(sys)nichols(sys,w)nichols(sys1,sys2,...,w)[mag,phase,w]=nichols(sys)单位负反馈系统的幅相频率特性语法格式：其中：mag表示幅值向量，phase表示相角向量，w表示自定义角频率num=[0.55];d1=[0.51];d2=[10.6/51];den=conv(d1,d2);nichols(num,den)已知传递函数：绘制系统的Nichols图程序代码例5自动控制理论研究中，分析最透彻的是二阶系统，一般高阶系统都用标准二阶系统等价进行研究，因为，高阶系统确定其动态性能指标是比较复杂的，工程上常采用闭环主导极点的概念对高阶系统进行近似分析，或从频域分析中获得频率特性参数，利用公式将其对应为二阶标准系统时域指标，达到降阶的目的。从而得出高阶系统的动态性能指标。五、控制系统频域设计

已知雕刻机系统开环传递函数：使用闭环频域响应法选择控制器增益K的值，完成控制指标。例6要求：（1）根据闭环系统的谐振峰值和谐振频率值确定时域指标，使系统超调量小于15%、在稳态误差2%的情况下，稳态时间小于15秒；（2）将高阶系统等价为标准二阶系统，输出等价后的二阶系统传递函数和动态特性参数；（3）画出控制前Bode图和控制后等价的二阶系统阶跃响应曲线。程序命令：forK=0.2:0.1:20G=tf(K,[1320]);G1=feedback(G,1);[mag,pha,w]=bode(G1);magn=mag(1,:);phase=pha(1,:);[M,i]=max(magn);Mr=20*log10(M);Pr=phase(1,i);Wr=w(i,1);temp=sqrt(1-1/M^2);zeta=sqrt((1-temp)/2);Wn=Wr/sqrt(1-2*zeta^2);ts=4.4/(zeta*Wn);Mp=exp((-pi*zeta)/sqrt(1-zeta^2));ifMp<0.2&ts<15break;endendifMp<0.2&ts<15G2=tf(Wn^2,[1,2*zeta*Wn,Wn^2]);disp(['K=',num2str(K),'系统满足了设计指标；'])disp(['放大系数K=',num2str(K)]);disp(['谐振峰值=',num2str(M)])disp(['谐振频率=',num2str(Wr)]);disp(['超调量=',num2str(Mp*100),'%'])disp(['稳态时间=',num2str(ts),'，

稳态误差2%'])elsedisp(['K=',num2str(K),'不满足设计指标，仅修改K值不能达到给定指标'])endsubplot(1,2,1);margin(G1);subplot(1,2,2);step(G2);结果：0.2694G2=

----------------------------

s^2+0.7105s+0.2694Continuous-timetransferfunction.K=0.7系统满足了设计指标；放大系数K=0.7谐振峰值=1.002谐振频率=0.13036超调量=5.2377%稳态时间=12.3858，

稳态误差2%原闭环系统的Bode图和控制后等价的二阶系统阶跃响应曲线如图所示。六、频域法校正设计1、超前校正设计（1）校正原理：超前校正的两个转折频率应分设在未校正系统的剪切频率的两侧。相频特性具有正相移，幅频特性具有正斜率。校正后，低频段不变，剪切频率比原系统大，使得快速性提高。（2）设计方法1）essK画出校正前原系统的Bode2）计算出稳定裕度及剪切频率PMk

c3）由m=PMd-PMk+(5°~10o),求值.4）L()=–10lgdB.确定校正后的系统的剪切频率m.5）根据m计算校正器的零极点的转折频率.6.画出校正后的系统Bode图,校验系统性能指标7.使系统构成闭环,验证闭环系统的响应已知单位负反馈系统被控对象的开环传函,设计校正器使之满足:（2）斜坡信号作用下,稳态误差小于0.001;（3）相角裕度在43°~48°之间（4）画出校正前后的阶跃响应曲线。（1）例7k0=1000;n1=1;d1=conv(conv([10],[0.11]),[0.0011]);G=tf(k0*n1,d1);[Gm,Pm,Wcg,Wcp]=margin(G);figure(1);margin(G);figure(3);step(feedback(G,1))fm=45-Pm+8;a=(1-sin(fm*pi/180))/(1+sin(fm*pi/180));[mag,pha,w]=bode(G);Lg=-10*log10(1/a);wmax=w(find(20*log10(mag(:))<=Lg));wmax1=min(wmax);wmin=w(find(20*log10(mag(:))>=Lg));wmin1=max(wmin);wm=(wmax1+wmin1)/2;T=1/(wm*sqrt(a));T1=a*T;Gc=tf([T,1],[T1,1])G1=Gc*G[Gm1,Pm1,Wcg1,Wcp1]=margin(G1);ifPm1>=45disp(['设计后相位裕量是：',num2str(Pm1),'相位裕量满足了设计要求'])elsedisp(['设计后相位裕量是：',num2str(Pm1),'相位裕量不满足设计要求'])endfigure(2);margin(G1);figure(4);

step(feedback(G1,1))校正前Bode图Gc=

0.01727s+1-------------0.00194s+1

Continuous-timetransferfunction.G1=

17.27s+1000---------------------------------------------1.94e-07s^4+0.0002959s^3+0.1029s^2+s

Continuous-timetransferfunction.设计后相位裕量是：46.6259相位裕量满足了设计要求校正后Bode图：校正前后的阶跃响应曲线（1）校正原理：

滞后校正的两个转折频率应比未校正系统的剪切频率小很多。相频特性具有负相移，幅频特性具有负斜率。校正后，低频段不变，剪切频率比原系统小，说明快速性变差，意味着牺牲系统的快速性换取稳定性。（2）设计方法2、滞后校正设计1）essK画出校正前原系统的Bode2）由(

gc)=-180°+PMd+5°,确定校正后的系统的剪切频率

gc3）L(

gc)=20lg,get=10(L/20).4）根据gc计算校正器的零极点的转折频率。（5）画出校正后的系统Bode图,校验系统性能指标（6）使系统构成闭环,验证闭环系统的响应已知单位负反馈系统被控对象的传函,设计校正器使之满足:（1）斜坡信号作用下,Kv30s-1;（2）相角裕度>45°;（3）穿越频率

c

2.3s-1例8K=30;num=1;den=conv(conv([10],[0.11]),[0.21]);Gp=tf(K*num,den);[mag,phase,w]=bode(Gp);fm=-180+45+10;figure(1);margin(Gp);figure(2);step(feedback(Gp,1))Wc1=spline(phase(:),w,fm);magdb=20*log10(mag(:));Lg=spline(w,magdb,Wc1);B=10^(-Lg/20);w1=0.1*Wc1;T=1/(B*w1);nc=[B*T,1];dc=[T,1];Gc=tf(nc,dc);printsys(nc,dc,'s');G=Gp*Gc;figure(3);margin(G)num1=G.num{1};den1=G.den{1};printsys(num1,den1,'s')[Gm,Pm1,Wcg,Wcp]=margin(G);Gm1=20*log10(Gm);ifGm1>=10&Pm1>=40&Wc1>=2.3disp(['设计后相位裕量：',num2str(Pm1),'，幅值裕量：',num2str(Gm1),'，满足了设计要求'])elsedisp(['设计后相位裕量是：',num2str(Pm1),'幅值裕量：',num2str(Gm1),'相位裕量或幅值裕量或穿越频率不满足设计要求'])end;figure(4);step(feedback(G,1))程序代码校正前Bode图num/den=

4.067s+1-------------43.2408s+1

num/den=

122.0107s+30

----------------------------------------0.86482s^4+12.9922s^3+43.5408s^2+s结果：设计后相位裕量：44.7078，幅值裕量：13.9098，满足了设计要求。校正前后的阶跃响应曲线（1）校正原理：实现滞后校正与超前校正的综合。滞后校正把剪切频率左移。从而减小了系统在剪切频率处的相位滞后，超前校正的作用是新剪切频率提供一个相位超前量，用以增大一个相位超前量，用以增大系统相位稳定裕度，使其满足动态性能要求3、滞后-超前校正设计（2）设计方法1）根据稳态性能要求，确定开环增益K；2）绘制未校正系统的对数幅频特性，求出未校正系统的截止频率ωc、相位裕量Pm和幅值裕量Gm；3）使中频段斜率为-20dB/dec，确定ωb，通常，在未校正系统对数幅频特性上，选择斜率从-20dB/dec变为-40dB/dec的转折频率作为校正网络超前部分的转折频率为ωb，这种选法可以降低已校正系统的阶次，且可保证中频区斜率为-20dB/dec，并占据较宽的频带。4）根据响应速度要求和校正后相位裕量Pm1的要求，选择校正后系统的开环截止频率Wc1。使校正网络中的1/Tb及α/Tb位于

Wc1的两侧。由公式先计算Kr，ts,（ts稳态时间一般系统指标值给出）

再计算校正后系统的开环截止频率，方法是：5）根据计算Tb和Wc1的值，计算衰减因子1/α：

其中：20lgα表示网络的最大幅值衰减量，L(Wc1)表示在角频率Wc1处的幅值；20lg(TbWc1)表示超前部分在Wc1处的幅值；6）根据相角裕量的Pm1要求，估算校正网络滞后部分的转折频率Wa。7）最后校验已校正系统的各项能指标。已知单位负反馈系统被控对象的传函,设计校正器使之满足:例9要求：Kv=20(1/s)，相位裕量Pm1=50(度)，调节时间ts≤4s，设计串联滞后-超前校正装置，使系统满足性能指标。

（1）（2）绘制原系统Bode图求出Wc1和αG0=tf(20,[0.1251.12510]);[gm,pm,wcg,wcp]=margin(G0)Pm1=50;ts=4;r=Pm1*pi/180;Kr=2+1.5*(1/sin(r)-1)+2.5*(1/sin(r)-1)^2Wc1=Kr*pi/ts

Wb=1;Tb=1/Wb;[mag,phase,w]=bode(G0);magdb=20*log10(mag);Lg=spline(w,magdb,Wc1)a1=(Lg+20*log10(Tb*Wc1))/20;a=10^(a1)结果：gm=0.4500pm=-13.9063wcg=2.8284wcp=4.1541Kr=2.6913Wc1=2.1137Lg=11.8483a=8.2693figure(1);step(feedback(G0,1)figure(2)bode(G0);程序代码为校正的阶跃响应及Bode图为：计算结果Wc1=2.1137,a=8.269，为计算方便，取Wc1=2.2,α=8。取Tb=1；（3）代入到校正传递函数：（4）令s=jWc1代入到总传递函数中：（5）令整个系统相位裕量为55o（取50+∆=55），则系统相角为：上式为：令：C=2*atan(0.275)-35*pi/180=-0.074计算的：ωa=0.2（6）代入到校正传递函数：（7）代入到总传递函数验证结果：G=tf([100,20],[0.625,10.156,40.25,1,0])[Gm,Pm,wg,wp]=margin(G)Margin(G)幅值裕量Gm=15.8DB；相位裕量Pm1=53.1o截止频率：Wcg=2.3rad/s交接频率：Wcp=7.8rad/s（8）验证各项性能指标：

由式5-15计算调节时间ts：Pm1=53.1>50;将Pm1=53.1，Wc1=2.3代入式5-15：

Bode图ts=(2+1.5*(1/sin(Pm1*pi/180)-1)+2.5*(sin(Pm1*pi/180)-1)^2)*pi/Wc1得到ts=3.5<4满足系统的指标要求。结果：38七、根轨迹1、根轨迹定义：

随着传递函数里某个参数从0变到无穷大，系统特征方程的特征根(也就是闭环极点)变化的轨迹。这里的某个参数通常是开环增益K，这时候是常规根轨迹，如果变化的参数不是K而是其他参数，那么叫做广义根轨迹，或者参数根轨迹。根轨迹的起点是开环传递函数的极点，终点是开环函数的零点，一个极点对一个零点，一个起点对应一个终点。如果极点个数大于零点个数，那么有部分根轨迹的终点终点变为无穷远处，这时称为无限零点，开环极点比零点多几个，就有几个无限零点。2、

绘制根轨迹的基本规则（1）根轨迹的方向、起点和终点

根据根轨迹定义，根轨迹起始于K*=0，终止于K

*→∞。由幅值条件得当K*=0时，s→-pj(j=1，2，…，n)为系统的开环极点；当K*→∞时，s→-zi(i=1，2，…，m)为系统的开环零点。结论：根轨迹起始于开环极点，终止于开环零点。如果开环零点数目m小于开环极点数目n，则有n-m条根轨迹终止于无穷远处。（2）根轨迹的分支数每个闭环特征根的变化轨迹都是整个根轨迹的一个分支，因此根轨迹的分支数与闭环特征方程的根的数目相同。

结论：根轨迹的分支数等于特征方程的阶次，也即开环零点数m和开环极点数n中的较大者。（3）根轨迹的连续性和对称性

K*的无限小增量与s平面上的长度|s+pj|及|s+zi|的无限小增量相对应，即复变量s在n条根轨迹上均有一个无限小的位移。当K*从零到无穷大连续变化时，根轨迹在s平面上一定是连续的。

特征根可以是实数根或复数根，而复数根又必然是成对出现的共轭复数，所以这些根必然对称于实轴。结论：根轨迹是连续的，且以实轴为对称的曲线。

（4）实轴上根轨迹的分布

若点s0右边零、极点个数之和为奇数，则s0点所在线段是根轨迹一部分；若点s0右边零、极点个数之和为偶数，则s0点所在的线段不是根轨迹。结论：实轴上属于根轨迹的部分，其右边开环零、极点的个数之和为奇数。（5）根轨迹的渐近线结论：如果系统的有限开环零点数m少于其开环极点数n，则当根迹增益K*→∞时，趋向无穷远处根轨迹的渐近线共有n-m条。这些渐近线与实轴上的交点坐标为与实轴正方向的夹角为设系统开环传递数为：

根据开环传递函数，得知开环极点数n=3，开环零点数m=0，首先将开环极点0，-1和-2标注在s平面上。由规则1）根轨迹有三条分支，分别起始于0，-1和-2，且这三条根轨迹都将趋向无穷远处；

2）实轴上根轨迹分布在０～-１以及-2～-∞之间；

3）根轨迹的渐近线共有n-m=3条，与实轴的交点和夹角计算公式如下：根轨迹在实轴上的分布、渐近线

（6）根轨迹的分离、会合点分离点：根轨迹分支在实轴上某点相遇又向复平面运动会合点：根轨迹分支从复平面运动到实轴上某点结论：根轨迹分离点或会合点的坐标，可通过求解方程得到（7）根轨迹与虚轴的交点与虚轴的交点可利用下面两种方法之一：方法一：用s=jω

代入特征方程求解方法二：根据系统临界稳定的条件，利用劳斯判据法求解结论：根轨迹与虚轴的交点坐标及临界根迹增益，可以通过用s=jω代入系统闭环特征方程求取，也可用劳斯判据列表方法确定。或44

在极点0和极点-1之间的根轨迹上一定有分离点存在，令d[G(s)H(s)]/ds＝0，整理后求得s1=-0.42（在根轨迹上，是分离点），s2=-1.58（不在根轨迹上，舍），所对应K*由幅值条件确定：

用劳斯判据法与虚轴交点：设系统开环传递数为:

(8)根轨迹的出射角和入射角出射角：根轨迹离开开环复数极点处的切线方向与正实轴之间的夹角入射角：根轨迹进入开环复数零点处的切线方向与正实轴之间的夹角结论：开环复数极、零点的出射角与入射角由下面公式计算

出射角：入射角：rlocus(G) %绘制根轨迹rlocus(G1,G2,…)%绘制多个系统的根轨迹

[r,k]=rlocus(G)

%得出闭环极点和对应的Kr=rlocus(G,k)%根据K得出对应的闭环极点[r,k]=rlocus(num,den)%返回闭环根矩阵r和对应的开环增益向量k。poles=rlocus(num,den,k1)

%增益k1对应的根轨迹极点。程序代码系统的开环传递函数为:：

num=1;den=[conv([1,4],conv([1-2+4i],[1-2-4i])),0]

G=tf(num,den)

rlocus(G) [r,k]=rlocus(G);%闭环极点和增益

[k,p]=rlocfind(G)%定位点的增益和极点sgrid(0.707,10)

%绘制ζ=0.707线和ωn=10线程序代码例1047系统的开环传递函数：

绘制系统的根轨迹如图。

num=[10.125];den=conv([1,0],conv([1,0],conv([1,5],conv([120],[150]))))[z,p,k]=tf2zp(num,den);G=zpk(z,p,k)rlocus(G)[r,k]=rlocus(G);

[k,p]=rlocfind(G)

sgrid(0.707,10)

程序代码例11八、根轨迹校正如果性能指标以单位阶跃响应的峰值时间、调整时间、超调量、阻尼系统、稳态误差等时域特征量给出时，一般采用根轨迹法校正。根轨迹法校正的基本思路为借助根轨迹曲线进行校正。由根轨迹的理论可知，如果系统的期望主导极点不在系统的根轨迹上。可添加开环零点或极点校正环节，适当选择零、极点的位置，就能够使得系统的根轨迹经过期望主导极点，且使在主导极点处满足给定要求，此时相当于相位超前校正。1、根轨迹校正的作用2、根轨迹超前校正（1）依据要求的系统性能指标，求出主导极点的期望位置。观察期望的主导极点是否位于校正前的系统根轨迹上。（2）确定校正网络零点方法是直接在期望的闭环极点位置下方（或在头两个实极点的左侧）增加一个相位超前网络的实零点。（3）确定校正网络极点方法是使期望的主导极点位于校正后的根轨迹上。利用校正网络极点的相角，使得系统在期望主导极点上满足根轨迹的相角条件。使用几何作图法来确定零极点位置的方法为步骤：过主导极点A|(si)与原点作直线OA，过主导极点A(si)作水平线PA；平分两线夹角<PAO作直线AB交负实轴于B点，由直线AB两边各分φ/2作射线交负实轴于C、D点

，左边交点C为极点-P，右边交点D为为零点-Z。如图所示。估计在期望的闭环主导极点处总的系统开环增益。计算稳态误差系数。如果稳态误差系数不满足要求，重复上述步骤。4）利用根轨迹设计相位超前网络时，超前网络的传递函数可表示为：其中：|z|<|p|，设计超前网络时，首先应根据系统期望的性能指标确定系统闭环主导极点的理想位置，然后通过选择校正网络的零、极点来改变根轨迹的形状，使得理想的闭环主导极点位于校正后的根轨迹上。某典型二阶系统的开环传递函数为：要求性能指标：

σ%≤25%，ts≤2s,试用根轨迹法确定超前校正装置。并用阶跃响应曲线比较校正前后系统的动态响应指标。步骤：（1）绘制原系统的根轨迹num=1;den=[1,2,0];rlocus(num,den);title('超前校正前根轨迹')axis([-3,0.5,-5,5]);校正前根轨迹如图所示。

例12（2）根据超调量、稳态时间指标，计算闭环主导极点：由自动控制理论中二阶系统阻尼比与超调量的关系表：

ζ00.10.150.20.250.30.40.50.707

σ%100%72.962.152.744.437.325.316.34.32根据公式和表，计算ωnclc;num=1;den=[1,2,0];rlocus(num,den);sgrid([0.5],[4])title(‘ζ=0.5，ωn=4‘);axis([-3,0.5,-5,5]);（3）画出原系统的根轨迹程序代码（4）设计校正网络，计算希望极点（5）根据期望的闭环极点，计算校正相角φ（6）在闭环主导极点处s1点为A，连接与原点直线OA，过A作横轴水平线AP，作角<PAO平分线AB交负实轴于B点，在AB两侧分别作角ϕ/2(15o)射线交负实轴，左边交点为极点P=-5.4，右边交点为零点Z=为-2.9，如图所示。（7）从图上得到校正网络参数校正后总传递函数：（8）根据校正后总传递函数，绘制根轨迹clc;clear;num=1;den=[1,2,0];sys1=tf(num,den);

num=[1,2.9];den=conv([1,2,0],[1,5.4]);sys2=tf(num,den);rlocus(sys1,sys2);

sgrid([0.5],[4])legend(‘超前校正前’,’超前校正前后’)title(‘超前校正前后的根轨迹')axis([-3,0.5,-5,5]);（9）根据模角条件确定K

4K=sqrt(4+3.46^2)*3.46*sqrt(3.4^2+3.46^2)/sqrt(3.46^2+0.9^2)=18.8K=4.7将s=-2+3.46j代入得：校正网络传递函数：（10）校正前后的时域分析：校正前闭环传函：校正后闭环传函：（11）绘制校正前后的单位阶跃响应曲线clc;num=4;den=[1,2,4];sys1=tf(num,den);num2=18.8*[1,2.9];den2=conv([1,2,0],[1,5.4]);den2=den2+[0,0,num2];sys2=tf(num2,den2);step(sys1,sys2)

校正前校正后结论：虽然校正后超调量略有提高，但小于25%，但上升时间大提幅度高。稳态时间达到了给定的2s要求。3、根轨迹滞后校正

滞后校正引入一对靠近原点的开环负实数偶极子，使根轨迹形状基本不改变，