基于Matlab控制系统频率特性分析法

本文主要介绍了基于Matlab控制系统的频率特性分析方法、频域稳定性判据以及开环频域性能分析，并获得频率响应曲线等。通过本章的学习，可以利用MATLAB对各种复杂控制系统进行频率分析，以此获得系统稳定性及其它性能指标。

一、频率特性基本概念

如果将控制系统中的各个变量看成是一些信号，而这些信号又是由许多不同频率的正弦信号合成的，则各个变量的运动就是系统对各个不同频率信号响应的总和。系统对正弦输入的稳态响应称频率响应。利用这种思想研究控制系统稳定性和动态特性的方法即为频率响应法。频率响应法的优点为：

⑴物理意义明确；

⑵可利用试验方法求出系统的数学模型，易于研究机理复杂或不明的系统，也适用于某些非线性系统；

⑶采用作图方法，非常直观。

1. 频率特性函数的定义

对于稳定的线性系统或者环节，在正弦输入的作用下，其输出的稳态分量是与输入信号相同频率的正弦函数。输出稳态分量与输入正弦信号的复数比，称为该系统或环节的频率特性函数，简称为频率特性，记作G（jω）=Y(jω)/R(jω)

对于不稳定系统，上述定义可以作如下推广。

在正弦输入信号的作用下，系统输出响应中与输入信号同频率的正弦函数分量和输入正弦信号的复数比，称为该系统或环节的频率特性函数。

当输入信号和输出信号为非周期函数时，则有如下定义。

系统或者环节的频率特性函数，是其输出信号的傅里叶变换像函数与输入信号的傅里叶变换像函数之比。

2. 频率特性函数的表示方法

系统的频率特性函数可以由微分方程的傅里叶变换求得，也可以由传递函数求得。这三种形式都是系统数学模型的输入输出模式。

当传递函数G(s)的复数自变量s沿复平面的虚轴变化时，就得到频率特性函数

G(jω)=G(s)|s=jω。所以频率特性是传递函数的特殊形式。

代数式：G(jω)=R(w)+jI(ω)

R(w)和I(w)称为频率特性函数G(jw)的实频特性和虚频特性。

指数式：G(jω)=A(w)eΦ(ω)

式中A(ω)=| G(jω)|是频率特性函数G(jw)的模，称为幅频特性函数。

Φ(w)=arg G(jω)是频率特性函数G(jω)的幅角，称为相频特性函数。

3．频率响应曲线

系统的频率响应可以用复数形式表示为G(jω)，常用的频率响应表示方法是图形表示法。根据系统频率响应幅值、相位和频率之间的不同显示形式，有伯德（Bode）图、奈魁斯特（Nyquist）图和尼柯尔斯（Nichols）图。

3.1 伯德图

伯德（Bode）图又称对数频率特性图，由对数幅频特性图和相频特性图组成。伯德图的横坐标为角频率ω，按常对数lgω分度。对数幅频特性的纵坐标是对数幅值。

L(ω)=20lg A(ω)，单位为分贝（dB），线性分度。

对数相频特性的纵坐标为φ(ω)，单位为度，线性分度。

一般情况下，控制系统开环对数频率特性图的绘制步骤如下：

⑴将开环频率特性按典型环节分解，并写成时间常数形式；

⑵求出各转角频率（交接频率），将其从小到大排列为ω1，ω2，ω3，……，并标注在ω轴上；

⑶绘制低频渐近线（ω1左边的部分），这是一条斜率为-20dB/decade（为系统开环频率特性所含1/jw因子的个数）的直线，它或者它的延长线应通过点（1,20lgK）；

⑷各转角频率间的渐近线都是直线，但自最小的转角频率ω1起，渐近线斜率发生变化，斜率变化取决于各转角频率对应的典型环节的频率特性函数。

例1 绘制一阶惯性环节G(s)=1/(4s+1)的伯德图。

程序代码如下：

num=1;

den=[4 1];

G=tf(num,den);

bode(G,'r')

3.2奈魁斯特图

奈魁斯特图又称为极坐标图或者幅相频率特性图。频率特性函数G(jω)的奈魁斯特图是角频率ω由0变化到∞时，频率特性函数在复平面上的图像。它以ω为参变量，以复平面上的向量表示G(jω)的一种方法。G(jω)曲线的每一点都表示与特定ω值相应的向量端点，向量的幅值为|G(jω)|，相角为argG(jω)；向量在实轴和虚轴上的投影分别为实频特性R(ω)和虚频特性I(ω)。

一般情况下，系统开环频率特性函数奈魁斯特图的绘制步骤如下：

⑴将系统的开环频率特性函数G0(jω)写成G(jω)=A(w)eΦ(ω)；

⑵确定奈魁斯特图的起点（ω＝0+）和（ω→+∞）。起点与系统所包含的积分环节个数（）有关，终点的A(ω)与系统开环传递函数分母和分子多项式阶次的差有关；

⑶确定奈魁斯特图与坐标轴的交点；

⑷根据以上的分析并且结合开环频率特性的变化趋势绘制奈魁斯特图。

例5-3 绘制一阶惯性环节G（s)=3/(5s+1)的奈魁斯特图。

程序代码如下：

G=tf(3,[5 1]);

nyquist(G);

hold on;

set(G,'inputdelay',5);

nyquist(G);

hold on;

set(G,'inputdelay',10);

nyquist(G);

hold on;

title('Nyquist图');

3.3 尼柯尔斯图

尼科尔斯图又称为对数幅频率特性图，它以开环频率特性函数的对数幅值为纵坐标，以相角值为横坐标，以角频率为参变量绘制的频率特性图。采用直角坐标。纵坐标表示20lg| G(jω)|，单位是dB，线性刻度。横坐标表示∠G(jω)，单位是度，线性分度。在曲线上一般标注角频率ω的值作为参变量。通常是先画出Bode图，再根据Bode图绘制尼科尔斯图。

4频率响应分析

时域分析中的性能指标直观反映控制系统动态相应的特征，属于直接性能指标，而系统频率特性函数的某些特征可以用作间接性能指标。

⑴开环频率特性的性能分析

基于开环频率特性函数的性能分析指标有如下两个：一是相角裕量γ，反映系统的相对稳定性；另一个是截止频率，反映系统的快速性。是A()=1所对应的角频率，或对数幅频特性图上L(ω)穿越0分贝线的斜率，在采用渐近线作图时，两者略有不同。

⑵闭环频率特性的性能分析

基于闭环频率特性函数的常用指标有两个：一是谐振峰值Mr，反映系统的相对稳定性；另一个是频带宽度或者带宽频率，定义为闭环幅频特性幅值M(ω)下降到0.707M(0)时对应的角频率，它反映了系统的快速性。

例2 用直接计算法，确定系统的谐振振幅和谐振频率。

已知一控制系统开环传递函数Go (s)=5.5/(s2+3s+5)，试求此系统的谐振振幅Mr和谐振频率ωr。

程序代码如下：

Go=tf(5.5,[1 3 5]);

[Mr,Pr,Wr]=mwr(Go)

mwr函数程序如下：

function [Mr,Pr,Wr]=mwr(G)

[mag,pha,w]=bode(G);

magn=mag(1,:);

phase=pha(1,:);

[M,i]=max(magn);

Mr=20*log10(M);

Pr=phase(1,i);

Wr=w(i,1);

运行结果：

Wr =    0.6915 

Mr =    0.8714 

Pr =  -24.46

结果中的单位分别是分贝（dB）、rad/s和度。 

例3 利用LTIView工具，获得系统频率响应的谐振振幅和谐振频率。

以例2中的传递函数为例，确定系统的谐振振幅Mr和谐振频率ωr。

在命令窗口界面键入如下命令：

>> Go=tf(5.5,[1 3 5])

>> ltiview

进入LTIView工具箱界面，对此系统进行分析。

Matlab控制仿真的简单例子

sys=tf(1,[3 1 0],'inputdelay',6/100); %10/s(s+2) tf-----建立传递函数模型

Ts=0.25; 

sysd=c2d(sys,Ts,'z'); % c2d ---连续的时间系统转化为离散的时间系统,即离散化.'t'表示双线性变换法

% Find critical valuecrs Kcr,Wcr % 'z'表示带零阶保持器的Z变换法,可缺省

% =============================

figure(1);

bode(sysd);grid on       % 画伯德图

[Gm,Pm,Wc,Wg]=margin(sysd)   % margin---计算增益和相位裕度函数

% Gm--增益裕度(即相位-180度处所对应的幅频特性绝对值的倒数)

% Wg--增益裕度对应的频率(即相位-180度处的频率)

% Wc--剪切频率(即幅频特性增益为0时的频率)

% Pm--相角裕度(即剪切频率处,使系统达到临界稳定状态时所需的附加相移)

Kcr=Gm;

Wcr=Wg;

Tcr=2*pi/Wcr;    % 求临界状态时对应的周期

pause % 暂停 

% Design PID Controller

% =====================

disp('parameters of PID')

Kp=0.6*Kcr;

Ti=0.5*Tcr;

Td=0.125*Tcr;

Ki=Kp*Ts/Ti;

Kd=Kp*Td/Ts;

disp('Model of PID Controller:')

sysc=Kp+tf(Ki*[1 0],[1 -1],Ts)+tf(Kd*[1 -1],[1 0],Ts)

% sysc=Kp+Ki*z/(z-1)+Kd*(z-1)/z 

% Check the Controller

% ====================

sysopen=sysc*sysd;

sysclose2=feedback(sysopen,1); % feedback-----求闭环反馈

figure(2);

step(sysclose2,'b') % step----------作单位阶跃响应函数曲线

grid on

hold on