MATLAB进行控制系统频域分析.docx
《MATLAB进行控制系统频域分析.docx》由会员分享，可在线阅读，更多相关《MATLAB进行控制系统频域分析.docx（22页珍藏版）》请在冰豆网上搜索。
MATLAB进行控制系统频域分析
MATLAB进行控制系统频域分析
一、基于MATLAB的线性系统的频域分析基本知识
（１）频率特性函数
。
设线性系统传递函数为：
则频率特性函数为：
由下面的MATLAB语句可直接求出G（jw）。
i=sqrt（-1） % 求取-1的平方根
GW=polyval（num，i*w）./polyval（den，i*w）
其中（num，den）为系统的传递函数模型。
而w为频率点构成的向量，点右除（./）运算符表示操作元素点对点的运算。
从数值运算的角度来看，上述算法在系统的极点附近精度不会很理想，甚至出现无穷大值，运算结果是一系列复数返回到变量GW中。
（２）用MATLAB作奈魁斯特图。
控制系统工具箱中提供了一个MATLAB函数nyquist（），该函数可以用来直接求解Nyquist阵列或绘制奈氏图。
当命令中不包含左端返回变量时，nyquist（）函数仅在屏幕上产生奈氏图，命令调用格式为：
nyquist（num,den）
nyquist（num,den,w）
或者
nyquist（G）
nyquist（G,w）
该命令将画出下列开环系统传递函数的奈氏曲线：
如果用户给出频率向量w,则w包含了要分析的以弧度/秒表示的诸频率点。
在这些频率点上，将对系统的频率响应进行计算，若没有指定的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考虑二阶典型环节:
试利用MATLAB画出奈氏图。
利用下面的命令，可以得出系统的奈氏图，如图5-1所示。
>>num=[0,0,1];
den=[1,0.8,1];
nyquist（num,den）
%设置坐标显示范围
v=[-2,2,-2,2];
axis（v）
grid
图5-1二阶环节奈氏图
title（′NyquistPlotofG（s）=1/（s^2+0.8s+1）′）
（3）用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）
公式logspace（d1，d2）是在指定频率范围内按对数距离分成50等分的，即在两个十进制数
和
之间产生一个由50个点组成的分量，向量中的点数50是一个默认值。
例如要在
弧度/秒与
弧度/秒之间的频区画伯德图，则输入命令时，
在此频区自动按对数距离等分成50个频率点，返回到工作空间中，即
w=logspace（-1，2）
要对计算点数进行人工设定，则采用公式logspace（d1，d2，n）。
例如，要在
与
之间产生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的分割。
例5.2 给定单位负反馈系统的开环传递函数为：
试画出伯德图。
利用以下MATLAB程序，可以直接在屏幕上绘出伯德图如图5-2。
>>num=10*[1,1];
den=[1,7,0];
bode（num,den）
grid
title（′BodeDiagramofG（s）=10*（s+1）/[s（s+7）]′）
该程序绘图时的频率范围是自动确定的，从0.01弧度/秒到30弧度/秒，且幅值取分贝值，
轴取对数，图形分成2个子图，均是自动完成的。
图5-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′）%用红线画
grid
title（′BodeDiagramofG（s）=10*（s+1）/[s（s+7）]′）
xlabel（¹Frequency（rad/s）¹）
ylabel（¹Gain（dB）¹）
%第二个图画伯德图相频部分。
subplot（2,1,2）;
semilogx（w,phase,¹-r¹）;
grid
xlabel（¹Frequency（rad/s）¹）
ylabel（′Phase（deg）′）
图5-3用户指定的频率点画出的伯德图
修改程序后画出的伯德图如5-3所示。
（4）.用MATLEB求取稳定裕量
同前面介绍的求时域响应性能指标类似，由MATLAB里bode（）函数绘制的伯德图也可以采用游动鼠标法求取系统的幅值裕量和相位裕量。
例如，我们可以在图20的幅频曲线上按住鼠标左键游动鼠标，找出纵坐标（Magnitude）趋近于零的点，从提示框图中读出其频率约为7.25dB。
然后在相频曲线上用同样的方法找到横坐标（Frequence）最接近7.25dB的点，可读出其相角为-53.9度，由此可得，此系统的相角裕量为126.1度。
幅值裕量的计算方法与此类似。
此外，控制系统工具箱中提供了margin（）函数来求取给定线性系统幅值裕量和相位裕量，该函数可以由下面格式来调用：
[Gm,Pm,Wcg,Wcp]=margin（G）;（67）
可以看出，幅值裕量与相位裕量可以由LTI对象G求出，返回的变量对（Gm,Wcg）为幅值裕量的值与相应的相角穿越频率，而（Pm,Wcp）则为相位裕量的值与相应的幅值穿越频率。
若得出的裕量为无穷大，则其值为Inf，这时相应的频率值为NaN（表示非数值），Inf和NaN均为MATLAB软件保留的常数。
如果已知系统的频率响应数据，我们还可以由下面的格式调用此函数。
[Gm,Pm,Wcg,Wcp]=margin（mag,phase,w）;
其中（mag,phase,w）分别为频率响应的幅值、相位与频率向量。
例5.3已知三阶系统开环传递函数为：
利用下面的MATLAB程序，画出系统的奈氏图，求出相应的幅值裕量和相位裕量，并求出闭环单位阶跃响应曲线。
>>G=tf（3.5,[1,2,3,2]）;
subplot（1,2,1）;
%第一个图为奈氏图
nyquist（G）;
grid
xlabel（'RealAxis'）
ylabel（'ImagAxis'）
%第二个图为时域响应图
[Gm,Pm,Wcg,Wcp]=margin（G）
G_c=feedback（G,1）;
subplot（1,2,2）;
step（G_c）
grid
xlabel（′Time（secs）′）
ylabel（′Amplitude′）
显示结果为：
ans=1.14291.1578
1.73211.6542
图5-4 三阶系统的奈氏图和阶跃响应图
画出的图形如图5-4所示。
由奈氏曲线可以看出，奈氏曲线并不包围（-1，j0）点，故闭环系统是稳定的。
由于幅值裕量虽然大于1，但很接近1，故奈氏曲线与实轴的交点离临界点（-1，j0）很近，且相位裕量也只有7.1578o，所以系统尽管稳定，但其性能不会太好。
观察闭环阶跃响应图，可以看到波形有较强的振荡。
如果系统的相角裕量γ>45o,我们一般称该系统有较好的相角裕量。
例5.4考虑一个新的系统模型，开环传递函数为：
由下面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.4365
Wcg= Wcp=
NaN100.3285
再输入命令
>>G_c=feedback（G,1）;
step（G_c）
grid
xlalel（′Time（sec）′）
ylalel（′Amplitude′）
图5-5 较理想的系统响应
可以看出，该系统有无穷大幅值裕量，且相角裕量高达85.4365o。
所以系统的闭环响应是较理想的，闭环响应图如图5-5.
（5）.时间延迟系统的频域响应
（1）时间延迟系统的传递函数模型
带有延迟环节e-Ts的系统不具有有理函数的标准形式，在MATLAB中，建立这类系统的模型。
要由一个属性设置函数set（）来实现。
该函数的调用格式为：
set（H,′属性名′,′属性值′）
其中H为图形元素的句柄（handle）。
在MATLAB中，当对图形元素作进一步操作时，只需对该句柄进行操作即可。
例如以下调用格式
h=plot（x,y）
G=tf（num,den）
Plot（）函数将返回一个句柄h，tf（）函数返回一个句柄G，要想改变句柄h所对应曲线的颜色，则可以调用下面命令：
Set（h,color,[1,0,0]）;
即对“color”参数进行赋值，将曲线变成红色（由[1,0,0]决定）
同样，要想对G句柄所对应模型的延迟时间’Td’进行修改，则可调用下面命令
Set（G,′Td′,T）
其中T为延迟时间。
由此修改后，模型G就已具有时间延迟特性。
（2）时间延迟系统的频域响应
含有一个延迟环节的系统，其开环频域响应为
可见，该系统的幅频特性不变，只加大了相位滞后。
例5.5考虑系统的开环模型为:
当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）
grid
figure %建立一个新的绘图窗口
step（G）
图5-6时间延迟系统奈氏图
图5-7时间延迟系统的阶跃响应
例5.6控制系统如图5-8所示，试分析其闭环特性。
图5-8
下面是该系统的演示过程，带符号％语句为注释语句。
%清文字屏
%clc
%清图形屏
%clf
%系统1描述
a1=[130];
b1=2;
%系统2描述
a2=[122];
b2=[12];
%两系统串联后的传递函数：
[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）
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]指将系统2置于反馈通道上：
[bb,ab]=cloop（b121,a12）;pause
[bb1,ab1]=feedback（b1,a1,b2,a2,-1）;pause
%求闭环频率响应
figure（6）
bode（bb1,ab1）;pause
%clg;
%闭环过渡过程
figure（7）
step（bb1,ab1）;pause
%clg
figure（8）
impulse（bb1,ab1）;pause
%clg;
%再求闭环零一极点，作为校核：
[zb,pb,kb]=tf2zp（bb,ab）;pause
printsys（zb,pb,kb）;
format
echooff
%演示结束
运行结果：
开环系统幅值裕量及相位裕量：
Gm=3.5435,Pm=55.022,Wg=1.6159,Wp=0.66898
在根轨迹上选择一点：
selected_point=-4.0158+0.0309i
此时闭环系统幅值裕量及相位裕量：
Gm=0.34677,Pm=-29.162,Wg=1.6159,Wp=2.5192
根据选定根轨迹上一点的位置可计算出相应的闭环零点，闭环极点及闭环增益值。
二、附录基于MATLAB的线性系统的频域分析的其它实例
例5.7已知单位负反馈系统前向通道的传递函数为：
试绘制出Bode图并计算系统的频域性能指标。
程序及结果如下：
num=[00281282];den=[151010510];
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,pm=38.692,wg=25.847,wc=1.249
kg=50.4272
例5.8已知一带延迟因子的系统开环传函为：
，
试求其有理传递函数的频率响应，同时在同一坐标中绘制以Pade近似延迟因子式系统的Bode图，并求此时系统的频域性能指标。
程序及结果如下：
n1=[1050];
d1=conv（[13],[13]）;
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）;
dent=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--'）;gridon
title（'bodeplot'）;
xlabel（'Frequency（rad／sec）'）;
ylabel（'Gaindb'）;
subplot（212）;grid
semilogx（w,phase1,w,phase2,'r--'）;gridon
xlabel（'Frequency（rad／sec）'）;
ylabel（'Phasedeg'）;
[gm1,pm1,wcp1,wcg1]=margin（s1）;
[gm2,pm2,wcp2,wcg2]=margin（s2）;
例5.9已知知两个单位负反馈系统开环传递函数分别为：
，
试用Bode图法判断闭环系统的稳定性。
程序如下：
num1=[0002.7];
den1=[1540];
sys1=tf（num1,den1）;
figure
（1）;holdon
[gm,pm,wcp,wcg]=margin（sys1）;
margin（sys1）;
title（'对数频率特性图1'）;
xlabel（'频率rad/sec'）;
ylabel（'GaindB'）;
num2=[0002.7];
den2=[15-40];
sys2=tf（num2,den2）;
figure
（2）;
[gm,pm,wcp,wcg]=margin（sys2）;
margin（sys2）;
title（'对数频率特性图2'）;
xlabel（'频率rad/sec'）;
ylabel（'GaindB'）;
结论：
系统1稳定，系统2不稳定（Warning:
Theclosed-loopsystemisunstable.）
升级会员 