频域响应分析

最新推荐文章于 2025-03-29 19:00:25 发布

纯良FOC

最新推荐文章于 2025-03-29 19:00:25 发布

阅读量5k

点赞数 4

分类专栏： MATLAB 控制理论MATLAB教程文章标签： matlab 开发语言

本文链接：https://blog.youkuaiyun.com/m0_46732561/article/details/126533666

版权

MATLAB 同时被 2 个专栏收录

13 篇文章

订阅专栏

控制理论MATLAB教程

6 篇文章

订阅专栏

用MATLAB获取増益裕度、相角裕度、相角穿越频率和增益穿越频率。

用MATLAB求取谐振峰值、谐振频率和带宽的方法。

用MATLAB绘制伯徳图

bode命令可以计算连续线性时不变(定常)系统频域响应的幅值和相角。

在MATLAB中输入指令bode(不带左端参数)，屏幕上就会画出伯德图，图中的幅值单位为分贝(dB)。最常用的bode命令形式有

bode(num,den)

bode(num,den,w)

bode(A,B,C,D)

bode(A,B,C,D,w)

bode(A, B,C,DZ iu,w)

bode(sys)

在使用左端参数时，例如采用

[mag,phase,w] = bode(num,den,w)

时，bode只将系统频域响应的数值赋予矩阵变量mag.phase和w,而不在屏幕上绘制伯德图。矩阵mag和phase分别包含系统频率响应在指定频率点上算得的幅值和相角。相角的单位为度(。)。幅值可以通过如下语句转换为分贝值：

magdB = 20*log10(mag)

带左端参数的bode命令的常用形式有

[mag,phase,w] = bode(num,den)

[mag,phase,w] = bode(num,den,w)

[mag,phase,w] = bode(A,BCD)

[mag,phase, w] = bode(A,B,C,D,w)

[mag,phase,w] = bode(A, B,C,D, iu,w)

[mag,phase, w] = bode(sys)

需指定频率范围，可使用命令logspae(d1, d2)或logspace(d1,de, n)。命令logspace (d1, d2)可以在十倍频程10^d1和10^d2内生成一个包含50个点的、呈对数均匀分布的向量(这50个频率点包含起点和终点;在起点和终点之间有48个点)。为了在频率0.100 rad/s之间产生 50个点，可以输入如下命令：

w = Iogspace(-1,2)

在绘制伯德图时，如果需要采用用户指定的频率点，那么bode命令必须像命令bode

(num, den, w)和[mag, phase, w] =bode(A, B, C, D, w)那样包含频率向量 w。

>> num = [25];
>> den = [1 4 25];
>> bode(num,den)

绘制伯德图为：

对于这个系统：

如果需要在0.01 ~1000rad/s频率范围内绘制伯德图，可以输入命令

w = logspace(—2,3,100)

该命令在0.01 -1000 rad/s范围内产生100个对数均匀分布的频率点(注意：向量w指定了计算频率响应的频率，其单位为rad/s)。

如果使用命令

bode(num,den,w)

绘制伯德图，则频率范围将和用户指定的相同，但频率响应的幅值范围和相角范围将会由

MATLAB自动设定。

>> num = [9 1.8 9];
>> den = [1 1.2 9 0];
>>  w = logspace(-2,3,100);
>>  bode(num,den,w)
>> grid
>> title('Bode Diagram of C(s) = 9(sA2 + 0.2s + 1 )/[s(s^2 + 1.2s + 9)]')

对于系统

$\dot{x}=Ax+Bu$

$y=Cx+Du$

该系统的伯德图可通过输入命令

bode（A,B,C,D）

或本小节开始所列的其他命令获得。

命令bode(A, B, C, D)可以绘制与每个系统输入相对应的一系列伯德图，其频率范围是自动设定的(当系统响应变化很快时，会采用较多的频率点来绘制伯德图)。

命令bode(A, B, C, D, iu)可以绘制从输入iu到该系统的所有输出(为y1,y2,…，)的伯德图，其中iu为系统的第i个输入。伯德图的频率范围是自动设定的(标量iu是系统输入的下标，它指定采用哪个输入来绘制伯德图)。假设控制信号向量u有3个输入，即

$\begin{bmatrix} u1\\ u2\\u3 \end{bmatrix}$

那么iu就必须设定为1, 2或3。

如果该系统仅有一个输入"，那么命令 bode(A,B,C,D) 或

bode(A,B,C,D,1)

都可以用来绘制伯德图。

可以分别求出四个传递函数，也可以直接采用命令：

>> A=[0 1;-25 -4];
>> B=[1 1;0 1];
>> C=[1 0;0 1];
>> D=[0 0;0 0];
>> bode(A,B,C,D)
>> grid

仿真图为:

用MATLAB绘制奈奎斯特图

正如伯德图一样，奈奎斯特图也常用于描述线性时不变反馈控制系统的频域响应。奈奎斯特图是极坐标图，而伯德图是直角坐标图。

关于奈奎斯特稳定性的分析，在此不再赘述。

MATLAB命令nyquist可以计算线性时不变连续系统的频域响应。如果在调用时不带左端参数，nyquist就在屏幕上绘制奈奎斯特图。

命令

nyquist(num,den)

绘制传递函数。

其他常用nyquist命令还有

nyquist(num,den,w)

nyquist(A,B,C,D)

nyquist(A, B,C, D,w)

nyquist(A, B, C, D, iu, w)

nyquist(sys)

当命令中包含了用户指定的频率向量w时，例如

nyquist(numz，den，w)

MATLAB就计算系统在以rad/s为单位的指定频率上的响应。

在使用带左端参数的命令时，例如

[re,im,w] = nyquist(num,den)

[re,im,w] = nyquist(num,den,w)

[re,im,w] = nyquist(A,B,C,D)

[re,im,w] = nyquist(A, B,C,D,w)

[re,im,w] = nyquist(A,B,C,D,iu,w)

[re,im,w] = nyquist(sys)

MATLAB将系统的频域响应的数值赋予矩阵re,im和w,但不在屏幕上绘制奈奎斯特图。矩阵re和im分别为系统在由向量w给定的频率点上的频域响应的实部和虚部。应当注意的是, re和im的列数与系统输出的数目相同，其每一行对应于w的一个元素。若要绘制奈奎斯特图，则可以使用plot指令。

代码如下：

>> num = [1 ];
>> den = [1 0.8 1];
>> nyquist(num,den)

如果希望依据自行设定的实轴和虚轴范围绘制奈奎斯特图，例如实轴和虚轴的范围均为从-2到+2,就可以输入命令

v= [-2 2 -2 2];

axis(v);

或者将以上两条命令合并为

axis([-2 2 -2 2]);

注意事项。在绘制奈奎斯特图时，如果MATLAB运算涉及“Divided by zein”(被零除), 所绘制的奈奎斯特图就可能是错误的，或者不适合进一步分析的需要。举例来说，如果传递函数为

G(s)= $\frac{1}{s(s+1)}$

那么MATLAB指令

num = [1];

den = [1 1 0];
nyquist(num,den)

就绘制出图5.14所示的一幅不合要求的奈奎斯特图。如果计算机上出现了这类图形，则可以通过指定axis(v)来修正图形。例如，在程序中输入命令

v = [-2 2 -5 5]; axis(v)

就能获得一个正确的奈奎斯特图。

不正确的奈奎斯特图

正确的奈奎斯特图

代码如下：

>> num = [1];
>> den = [1 1 0];
>> nyquist(num,den)
>>  v = [-2 2 -5 5]; axis(v)

如果只需要频率为正的部分所对应的曲线，就需要使用命令

[re,im,w]=nyquist(num,den,w)

将频率区间划分为几个具有不同增量因子的区间。

具体实现代码如下：

num = [20 20 10];
den = [1 11 10 0];
w1 = 0.1:0.1:10;w2 = 10:2:100; w3 = 100:10:500;
w = [w1 w2 w3];
[re,im,w] = nyquist(num,den,w);
plot(re,im)
v = [-3 3 -5 1]; axis(v)
grid
title('Nyquist Plot of G(s) = 20(s^2 + s + 0.5)/[s(s + 1)(s + 10)]')
xlabel('Real Axis')
ylabel('Imaginary Axis')

上述例程中，可以用ww = logspace(—1,2,100)代替

w1 = 0.1:0.1:10;w2 = 10:2:100; w3 = 100:10:500;
w = [w1 w2 w3];

再看一个例子：

附上源程序：

>> num = [20 20 10];
>> den = [1 11 10 0];
>> ww = logspace(-1,2,100);
>> [re,im,ww] = nyquist (num,den,ww);
>> plot(re,im)
>>  v = [-2 3 -5 0]; axis(v)
>>  grid
>> title('Nyquist Diagram')
>>  xlabel('Real Axis')
>> ylabel('Imaginary Axis')
>> hold
>> w = [0.2 0.3 0.5 1 2 6 10 20];
>>  [re,im,w] = nyquist(num,den,w);
>>  plot(re,im,'o')
>>  text(1.1,-4.8,'w = 0.2')
>> text(1.1,-3.1,'w = 0.3')
>> text(1.25,-1.7,'0.5')
>> text(1.37,-0.4,'1')
>> text(1.8,-0.3,'2')
>> text(1.4,-1.1,'6')
>> text(0.77,-0.8,'10')
>> text(0.037,-0.8,'20')
>> [mag,phase,w] = bode(num,den,w);
>> [w mag phase]

运行结果为：

ans =

0.2000 4.9176 -78.9571
0.3000 3.2426 -72.2244
0.5000 1.9975 -55.9925
1.0000 1.5733 -24.1455
2.0000 1.7678 -14.4898
6.0000 1.6918 -31.0946
10.0000 1.4072 -45.0285
20.0000 0.8933 -63.4385

用MATLAB绘制尼科尔斯图。

用于绘制尼科尔斯图的MATLAB命令有

nichols(num,den)

nichols(num,den,w)

nichols(A,B,C,D)

nichols(A,B,C,D,w)

nichols(A, B,C, D, iu,w)

nichols(sys)

[mag,phase,w]=nichols(num,den)

[mag,phase,w]=nichols(num,den,w)

[mag,phase,w]=nichols(A,B,C,D)

[mag,phase,w]=nichols(A,B,C,D,w)

[mag,phase,w]=nichols(A, B,C, D, iu,w)

[mag, phase, w]=nichols(sys)

以上任一个命令都能绘制尼科尔斯图。使用方法与前面大同小异。如果所用的尼科尔斯命令没有涉及任何左端参数，就会自动生成尼科尔斯图。如果带有左端参数，这些命令就不绘制尼科尔斯图，而是返回幅值和相角(单位为度)的数值。为了绘制曲线，必须使用plot命令。命令ngrid可用于添加尼科尔斯线图网格线。

相角裕度和增益裕度

下图给出了 G(jw)在3种不同开环增益K下的极坐标图。对于较大的增益K,系统不稳定。当增益降低到一定数值时，G(js)轨迹穿越-1+j0点。这意味着当增益取该数值时, 系统濒临不稳定边缘，表现为持续振荡。对于较小的增益K,系统是稳定的。

一般而言，G(js)的轨迹离-1+j0点越近，系统响应就越振荡。G(jw)轨迹与- 1 +j0点的接近程度可以作为系统稳定裕度的度量(不过，这并不适用于条件稳定系统)。通常用相角裕度和增益裕度来描述这种接近程度。

关于相角裕度和增益裕度的概念不再赘述，可以参考自动控制原理课本。

对于稳定的最小相位系统，增益裕度表明增益要增加多少才会使系统变得不稳定。对于不稳定系统，增益裕度表明增益要下降多少才能使系统变得稳定。

一阶或二阶系统的增益裕度为无穷大，这是因为这些系统的极坐标图不穿越负实轴。

因此,从理论上讲，一阶或二阶系统不可能不稳定（但值得注意的是，所谓的-一阶或二阶系统仅仅是实际系统的近似，因为在推导系统方程时忽略了一些小的时延。如果考虑小的时延，这些一阶或二阶系统就可能变得不稳定）。

条件稳定系统有两个或更多的相角穿越频率，一些具有复杂动态特性的高阶系统可能有两个或更多的增益穿越频率。对于有两个或更多增益穿越频率的稳定系统，相角稳定裕度应在最高增益穿越频率点加以测量。

用MATLAB获取増益裕度、相角裕度、相角穿越频率和增益穿越频率。

MATLAB可以通过

如下命令方便地求取增益裕度、相角裕度、相角穿越频率和增益穿越频率：

[Gm, Pm, wcp, wcg] = margin(sys)

其中Gm是增益裕度，Pm是相角裕度，wcp是相角穿越频率，wcg是增益穿越频率，下例给出了该指令的详细使用方法。

>>  num = [20 20];
>>  den = conv([1 5 0],[1 2 10]);
>> sys = tf(num,den);
>>  w = logspace(-1,2,100);
>> bode(sys,w)
>>  grid;
>> [Gm, Pm,wcp,wcg] = margin(sys);
>>  GmdB = 20*log10(Gm);
>> [GmdB Pm wcp wcg]

运行结果为：

ans =

9.9301 103.6573 4.0132 0.4426

用MATLAB求取谐振峰值、谐振频率和带宽的方法。

谐振峰值是闭环系统频域响应幅值(以分贝为单位)的最大值。谐振频率是产生该最大幅值的频率点。以下MATLAB命令可用于获取谐振峰值和谐振频率：

[mag,phase,w] = bode(num,den,w); 或 [mag,phase,w] = bode(sys,w);

[Mp,k] = max(mag);

resonant_peak = 20*log10(Mp);

resonant_frequency = w(k)

可以输入如下几行语句来求取带宽：

n = 1;

while 20*log10(mag(n)) > = -3; n = n + 1;

end

bandwidth = w(n)

下面给出一个具体的MATLAB程序：

>> num = [1];
>>  den = [0.5 1.5 1 0];
>>  sys1 = tf(num,den);
>> sys = feedback(sys1,1);
>> w = logspace(-1,1);
>> bode(sys,w);
>> grid
>> [imag,phase,w] = bode(sys,w);
>> [Mp,k] = max(mag);
>> [Mp,k] = max(mag);
>>  resonant_peak = 20*log10(Mp)