沈阳理工大学MATLAB论文
>> syms t s c r G R C omega; >> r=sym('Heaviside(t)'); >> R=laplace(r);
>> c=1+1.5*exp(-3*t)+1.2*exp(-5*t); >> C=laplace(c); >> C=factor(C); >> G=C/R;
>> G=subs(G,s,j*omega) G =
1/10*(-37*omega^2+191*i*omega+150)/(i*omega+3)/(i*omega+5) 即系统的频率响应为 G(j?)?1s?2??ns??22n1
10(?37?2?191j??150)(j??3)(j??5)2.试计算振荡环节G(s)?的频率特性。
其实现的程序代码如下：
>> syms zet omega omegan real;
>> G=omegan^2/((j*omega)^2+2*zet*omegan*(j*omega)+omegan^2); >> Aabs=abs(G) Aabs =
omegan^2/(omega^4-2*omega^2*omegan^2+omegan^4+4*zet^2*omegan^2*omega^2)^(1/2)
>> [Re]=simple(real(G)) Re =
(-omega^2*omegan^2+omegan^4)/(omega^4-2*omega^2*omegan^2+omegan^4+4*zet^2*omegan^2*omega^2) >> [Im]=simple(imag(G)) Im =
-2*zet*omegan^3*omega/(omega^4-2*omega^2*omegan^2+omegan^4+4*zet^2*omegan^2*omega^2)
6
沈阳理工大学MATLAB论文
>> Gang=simple(atan(Im/Re)) Gang =
-atan(2*zet*omegan*omega/(-omega^2+omegan^2)) 即振荡环节的实频特性与虚频特性分别为
2?n?a?jb G(j?)?22???2??n(j?)??n222?n(?n??2)?n*2???na?2b?, 2222222(?n??)?(2???n)(?n??)?(2???n)还有振荡环节的幅频特性与相频特性分别为
A(?)?G(j?)?2?n2(?n??2)2?(2???n)2?
????1???????n????21??????2???2?n??2
????2?(?)??G(j?)??arctan2???n??arctan22?n??2???n????2??1?????n
3.在某反馈系统中，G(s)?5 ， H(s)?1?Kns(Kn?0) 。试确定闭环系
s(s?8)统稳定时反馈参数Kn的临界值。
解：（1）求闭环系统传递函数对应的实频与虚频特性。 其实现的程序代码如下： >> syms s g h u v;
>> syms kn omega omegac real; >> s=j*omega; >> G=10/(s*(s-10)); >> H=1+kn*s; >> GH=G*H;
>> U=factor(real(GH)) U =
7
沈阳理工大学MATLAB论文
-10*(1+10*kn)/(omega+10*i)/(omega-10*i) >> v=factor(imag(GH)) v =
-10*(-10+omega^2*kn)/(omega+10*i)/(omega-10*i)/omega
（2）当闭环系统处于临界稳定时，开环系统的频率响应G(j?)H(j?)，即Nyquist曲线将通过[G(j?)H(j?)]平面上的点（-1，j0）,此时???c。那么有
?V(?c)?0 ?U(?)??1C?其实现的程序代码如下： >> syms kn omega omegac real;
>>[kn,omegac]=solve('10*(1+10*kn)/(-omegac+10*j)=-1','10*(-10*kn*omegac^2)/(o megac+10*j)/omegac/(omegac+10*j)=0',kn,omegac) kn =
-1/10-1/10*j omegac = 10+10*j 0
即闭环系统稳定时反馈参数Kn=1，此时?c?10rad/s。
三．频率响应分析的MATLAB实现
前面已经介绍了频率响应分析的几个概念，下面将对频率响应分析的MATLAB实现展开介绍。 1.Bode图的绘制
MATLAB提供了一条直接求解和绘图系统的Bode图的函数bode()和一条直接求解系统的幅值稳定裕度和相位稳定裕度的函数margin()，其基本调用格式为 bode(sys) bode(sys,w)
[mag,phase,w]= bode(sys)
bode()函数用来计算并绘制系统的Bode图，可用于SISO或MIMO的连续时间系
8
沈阳理工大学MATLAB论文
统。当函数命令为等式左边输出变量的格式，函数在当前窗口中直接绘制出系统的Bode图。
LTI对象sys可以是由函数tf()、zpk()、ss()中任何一个函数建立的系统模型。w用来定义绘制Bode图时的频域范围或者频率点。如果定义频率范围，w必须为[wmin,,wmax]格式；如果定义频率点，则w必须是需要频率点构成的向量。第三条语句只计算系统Bode图的输出数据，而不绘制曲线。mag为系统Bode图的振幅值，phase为Bode图的相位值。
margin(sys)
[Gm,Pm,Wcg,Wcp]= margin（sys） [Gm,Pm,Wcg,Wcp]= margin（mag,phase,w）
margin函数可以从频率响应数据中计算出幅值稳定裕度、相位稳定裕度及其对应的角频率。当不带输出变量引用函数时，margin()函数可在当前图形窗口中绘制出带有稳定裕Bode图。margin（mag,phase,w）函数可以在当前窗口中绘制出带有系统幅值裕度与相位裕度的Bode图。其中，mag,phase及w分别为由bode或dbode求出的幅值裕度、相位裕度及其对应的角频率。 4.已知一个高阶系统的开环传递函数为
G(s)?K(0.00167s?1)
s(0.03s?1)(0.0025s?1)(0.001s?1)试计算当开环增益K=5, 500, 800, 3000时，系统稳定裕量的变化。 其实现的程序代码如下： >> k=[5,500,800,3000]; >> for j=1:4
num=k(j)*[0.0167 1];
den=conv(conv([1 0],[0.03 1]),conv([0.0025,1],[0.001,1])); G=tf(num,den); y(j)=allmargin(G); end >> y(1) ans =
9
沈阳理工大学MATLAB论文
GMFrequency: 602.4232 GainMargin: 455.2548 PMFrequency: 4.9620 PhaseMargin: 85.2751 DMFrequency: 4.9620 DelayMargin: 0.2999 Stable: 1 >> y(2) ans =
GMFrequency: 602.4232 GainMargin: 4.5525 PMFrequency: 237.7216 PhaseMargin: 39.7483 DMFrequency: 237.7216 DelayMargin: 0.0029 Stable: 1 >> y(3) ans =
GMFrequency: 602.4232 GainMargin: 2.8453 PMFrequency: 329.9063 PhaseMargin: 27.7092 DMFrequency: 329.9063 DelayMargin: 0.0015 Stable: 1 >> y(4) ans =
GMFrequency: 602.4232 GainMargin: 0.7588
10