MATLAB进行控制新版系统频域分析.doc

一、基于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(′Nyquist Plot of G(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(′Bode Diagram of G(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(′Bode Diagram of G(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('Real Axis')
ylabel('Imag Axis')
% 第二个图为时域响应图
[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.1429 1.1578
1.7321 1.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 =
　　　　　　　　NaN 100.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=[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)
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
echo off
%演示结束
运行结果：
开环系统幅值裕量及相位裕量：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=[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, pm=38.692, wg=25.847, wc=1.249
kg=50.4272

例5.8 已知一带延迟因子系统开环传函为：，
试求其有理传输函数频率响应，同时在同一坐标中绘制以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);
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--');grid on
title('bode plot');
xlabel('Frequency（rad／sec）');
ylabel('Gain db');
subplot(212);grid
semilogx(w,phase1,w,phase2,'r--');grid on
xlabel('Frequency（rad／sec）');
ylabel('Phase deg');
[gm1,pm1,wcp1,wcg1]=margin(s1);
[gm2,pm2,wcp2,wcg2]=margin(s2);

例5.9 已知知两个单位负反馈系统开环传输函数分别为：
，
试用Bode图法判定闭环系统稳定性。
程序以下：
num1=[0 0 0 2.7];
den1=[1 5 4 0];
sys1=tf(num1,den1);
figure(1);hold on
[gm,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.）