北京科技大学控制实验报告4

实验四线性系统的频域分析一.实验目的（1）熟练掌握运用MATLAB命令绘制控制系统伯德图的方法；（2）了解系统伯德图的一般规律及其频域指标的获取方法；（3）熟练掌握运用伯德图分析控制系统稳定性的方法；（4）设计超前校正环节并绘制Bode图；（5）设计滞后校正环节并绘制Bode图。二.实验原理及内容1、频率特性函数)(jG。设线性系统传递函数为：nnnnmmmmasasasabsbsbsbsG1101110)(则频率特性函数为：nnnnmmmmajajajabjbjbjbjwG)()()()()()()(1101110由下面的MATLAB语句可直接求出G(jw)。i=sqrt(-1)%求取-1的平方根GW=polyval(num，i*w)./polyval(den，i*w)其中（num，den）为系统的传递函数模型。而w为频率点构成的向量，点右除（./）运算符表示操作元素点对点的运算。从数值运算的角度来看，上述算法在系统的极点附近精度不会很理想，甚至出现无穷大值，运算结果是一系列复数返回到变量GW中。2、用MATLAB作伯德图控制系统工具箱里提供的bode()函数可以直接求取、绘制给定线性系统的伯德图。当命令不包含左端返回变量时，函数运行后会在屏幕上直接画出伯德图。如果命令表达式的左端含有返回变量，bode()函数计算出的幅值和相角将返回到相应的矩阵中，这时屏幕上不显示频率响应图。命令的调用格式为：[mag,phase,w]=bode(num,den)[mag,phase,w]=bode(num,den,w)或[mag,phase,w]=bode(G)[mag,phase,w]=bode(G,w)矩阵mag、phase包含系统频率响应的幅值和相角，这些幅值和相角是在用户指定的频率点上计算得到的。用户如果不指定频率w,MATLAB会自动产生w向量，并根据w向量上各点计算幅值和相角。这时的相角是以度来表示的，幅值为增益值，在画伯德图时要转换成分贝值，因为分贝是作幅频图时常用单位。可以由以下命令把幅值转变成分贝：magdb=20﹡log10(mag)绘图时的横坐标是以对数分度的。为了指定频率的范围，可采用以下命令格式：logspace(d1，d2)或logspace(d1，d2，n)第一种格式是在指定频率范围内按对数距离分成50等分的，即在两个十进制数1101d和2102d之间产生一个由50个点组成的分量，向量中的点数50是一个默认值。例如要在1.01弧度/秒与1002弧度/秒之间的频区画伯德图，则输入命令时，)(log1101d,)(log2102d在此频区自动按对数距离等分成50个频率点，返回到工作空间中，即w=logspace(-1，2)要对计算点数进行人工设定，则采用第二种格式。例如，要在11与10002之间产生100个对数等分点，可输入以下命令：w=logspace(0，3，100)在画伯德图时，利用以上各式产生的频率向量w，可以很方便地画出希望频率的伯德图。由于伯德图是半对数坐标图且幅频图和相频图要同时在一个绘图窗口中绘制，因此，要用到半对数坐标绘图函数和子图命令。（1）对数坐标绘图函数利用工作空间中的向量x，y绘图，要调用plot函数，若要绘制对数或半对数坐标图，只需要用相应函数名取代plot即可，其余参数应用与plot完全一致。命令公式有：semilogx(x，y，s);只对x轴进行对数变换，y轴仍为线性坐标。semilogy(x，y，s);y轴取对数变换的半对数坐标图。Loglog(x，y，s);全对数坐标图，即x轴和y轴均取对数变换。（2）子图命令MATLAB允许将一个图形窗口分成多个子窗口，分别显示多个图形，这就要用到subplot()函数，其调用格式为：subplot(m，n，k)该函数将把一个图形窗口分割成m×n个子绘图区域，m为行数，n为列数，用户可以通过参数k调用各子绘图区域进行操作，子图区域编号为按行从左至右编号。对一个子图进行的图形设置不会影响到其它子图，而且允许各子图具有不同的坐标系。例如，subplot(4，3，6)则表示将窗口分割成4×3个部分。在第6部分上绘制图形。MATLAB最多允许9×9的分割。【范例6-3】给定单位负反馈系统的开环传递函数为：)7()1(10)(ssssG，试画出伯德图。利用以下MATLAB程序，可以直接在屏幕上绘出伯德图如图6-2。num=10*[1,1];den=[1,7,0];bode(num,den)gridtitle(′BodeDiagramofG(s)=10*(s+1)/[s(s+7)]′)该程序绘图时的频率范围是自动确定的，从0.01弧度/秒到30弧度/秒，且幅值取分贝值，轴取对数，图形分成2个子图，均是自动完成的。如果希望显示的频率范围窄一点，则程序修改为：num=10*[1,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(′BodeDiagramofG(s)=10*(s+1)/[s(s+7)]′)xlabel(¹Frequency(rad/s)¹)ylabel(¹Gain(dB)¹)%第二个图画伯德图相频部分。subplot(2,1,2);semilogx(w,phase,¹-r¹);gridxlabel(¹Frequency(rad/s)¹)ylabel(′Phase(deg)′)图6-2自动产生频率点画出的伯德图修改程序后画出的伯德图如6-3所示。【范例6-4】已知控制系统开环传递函数：210()()210GsHsss，绘制其Bode图并确定谐振峰值Mr和谐振频率r。绘制Bode图（略），插值函数spline(m,w,mr)求谐振频率5、用MATLAB求取稳定裕量同前面介绍的求时域响应性能指标类似，由MATLAB里bode()函数绘制的伯德图也可以采用游动鼠标法求取系统的幅值裕量和相位裕量。此外，控制系统工具箱中提供了margin()函数来求取给定线性系统幅值裕量和相位裕量，该函数可以由下面格式来调用：margin（num,den）;给定开环系统的数学模型，作Bode图，并在图上方标注幅值裕度Gm和对应频率ωg，相位裕度Pm和对应的频率ωc。[Gm,Pm,Wcg,Wcp]=margin(G);可以看出，幅值裕量与相位裕量可以由LTI对象G求出，返回的变量对（Gm,Wcg）为幅值裕量的值与相应的相角穿越频率，而（Pm,Wcp）则为相位裕量的值与相应的幅值穿图6-3用户指定的频率点画出的伯德图越频率。若得出的裕量为无穷大，则其值为Inf，这时相应的频率值为NaN(表示非数值)，Inf和NaN均为MATLAB软件保留的常数。如果已知系统的频率响应数据，我们还可以由下面的格式调用此函数。[Gm,Pm,Wcg,Wcp]=margin(mag,phase,w);其中(mag,phase,w)分别为频率响应的幅值、相位与频率向量。【范例6-6】考虑一个新的系统模型，开环传递函数为：)9)(1()5(100)(22sssssG，由下面MATLAB语句可直接求出系统的幅值裕量和相位裕量:G=tf(100*conv([1,5],[1,5]),conv([1,1],[1,1,9]));[Gm,Pm,Wcg,Wcp]=margin(G)结果显示Gm=Pm=Inf85.4365Wcg=Wcp=NaN100.3285再输入命令G_c=feedback(G,1);step(G_c)gridxlalel(′Time(sec)′)ylalel(′Amplitude′)可以看出，该系统有无穷大幅值裕量，且相角裕量高达85.4365o。所以系统的闭环响应是较理想的，闭环响应图如图6-5。【自我实践6-1】某单位负反馈系统的开环传递函数()(1)(2)kGssss，求(1)当k=4时，计算系统的增益裕度，相位裕度，在Bode图上标注低频段斜率，高频段斜率及低频段、高频段的渐近相位角。(2)如果希望增益裕度为16dB，求出响应的k值，并验证(1)图6-5较理想的系统响应clcnum=[4];den=[1,3,2,0];G=tf(num,den)[Gm,Pm,Wcg,Wcp]=margin(G)bode(num,den)gridTransferfunction:4-----------------s^3+3s^2+2sGm=1.5000Pm=11.4304Wcg=1.4142Wcp=1.1431低频段斜率-20dB/10dec高频段-60dB/10dec低频段渐进相位角-90度，高频段-270度。增益裕度Gm=1.5000dB/10dec相位裕度Pm=11.4304度（2）K=0.951clcnum=[0.951];den=[1,3,2,0];G=tf(num,den)[Gm,Pm,Wcg,Wcp]=margin(G)bode(num,den)gridTransferfunction:0.951-----------------s^3+3s^2+2sGm=6.3091Pm=54.7839Wcg=1.4142Wcp=0.42766、系统对数频率稳定性分析【自我实践6-2】系统开环传递函数()(0.51)(0.11)kGssss，试分析系统的稳定性。经计算，当k=12时系统的增益裕度、相角裕度为零。clcnum=[12];den=[0.05,0.6,1,0];G=tf(num,den)[Gm,Pm,Wcg,Wcp]=margin(G)bode(num,den)gridTransferfunction:12----------------------0.05s^3+0.6s^2+sGm=1Pm=9.5374e-006Wcg=4.4721Wcp=4.4721当k12时clcnum=[1];den=[0.05,0.6,1,0];G=tf(num,den)[Gm,Pm,Wcg,Wcp]=margin(G)bode(num,den)gridTransferfunction:1----------------------0.05s^3+0.6s^2+sGm=12.0000Pm=60.4231Wcg=4.4721Wcp=0.9070系统产生衰减震荡，增益裕度和相角裕度大于0，系统稳定。当K12时clcnum=[50];den=[0.05,0.6,1,0];G=tf(num,den)[Gm,Pm,Wcg,Wcp]=margin(G)bode(num,den)gridTransferfunction:50----------------------0.05s^3+0.6s^2+sGm=0.2400Pm=-27.5774Wcg=4.4721Wcp=8.5945增益裕度相角裕度均小于0，系统不稳定。【自我实践6-3】某单位负反馈系统的开环传递函数31.6()(0.011)(0.11)Gssss，求(1)绘制Bode图，在幅频特性曲线上标出低频段斜率、高频段斜率、开环截止频率和中频段穿越频率；在幅频特性曲线标出：低频段渐近相位角、高频段渐近相位角和-180线的穿越频率。(2)计算系统的相位裕度和幅值裕度h，并确定系统的稳定性（1）低频段斜率-20dB/10dec高频段-60dB/10dec,幅值穿越频率Wcp=16.3053低频段渐进相位角-90度，高频段渐进相位角-270度；-180度线穿越频率为Wcg=31.6228clcnum=[31.6];den=[0.001,0.11,1,0];G=tf(num,den)[Gm