Chapter7.2：MATLAB在频率法中的应用及频率法稳定性分析

该系列博客主要讲述Matlab软件在自动控制方面的应用，如无自动控制理论基础，请先学习自动控制系列博文，该系列博客不再详细讲解自动控制理论知识。
自动控制理论基础相关链接：https://blog.csdn.net/qq_39032096/category_10287468.html?spm=1001.2014.3001.5482
博客参考书籍：《MATLAB/Simulink与控制系统仿真》。

2.MATLAB在频率法中的应用及频率法稳定性分析

2.1 MATLAB在频率法中的应用

2.1.1 求取和绘制频率响应曲线相关的函数

1. $\mathrm{Nyquist}$曲线绘制函数$\mathrm{nyquist}()$。
   • $\mathrm{nyquist}(\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d})$：绘制系统的一组$\mathrm{Nyquist}$曲线，每条曲线对应于连续状态空间系统$[\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d}]$的输入/输出组合对，其中频率范围由函数自动选取，且在响应快速变化的位置会自动采用更多取样点；
   • $\mathrm{nyquist}(\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d},\mathrm{iu})$：绘制从系统第$\mathrm{iu}$个输入到所有输出的极坐标图；
   • $\mathrm{nyquist}(\mathrm{num},\mathrm{den})$：绘制以连续时间多项式传递函数表示的系统极坐标图；
   • $\mathrm{nyquist}(\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d},\mathrm{iu},\mathrm{w})$或$\mathrm{nyquist}(\mathrm{num},\mathrm{den},\mathrm{w})$：利用指定的角频率矢量绘制系统的极坐标图；
2. $\mathrm{Bode}$图绘制函数$\mathrm{bode}()$。
   • $\mathrm{bode}(\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d})$：绘制系统的一组$\mathrm{Bode}$图，是针对连续状态空间系统$[\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d}]$的每个输入的$\mathrm{Bode}$图，其中频率范围由函数自动选取，且在响应快速变化的位置会自动采用更多采样点；
   • $\mathrm{bode}(\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d},\mathrm{iu})$：绘制从系统第$\mathrm{iu}$个输入到所有输出的$\mathrm{Bode}$图；
   • $\mathrm{bode}(\mathrm{num},\mathrm{den})$：绘制以连续时间多项式传递函数表示的系统$\mathrm{Bode}$图；
   • $\mathrm{bode}(\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d},\mathrm{iu},\mathrm{w})$或$\mathrm{bode}(\mathrm{num},\mathrm{den},\mathrm{w})$：利用指定的角频率矢量绘制系统的$\mathrm{Bode}$图；
3. $\mathrm{Nichols}$曲线绘制函数$\mathrm{nichols}()$。
   • $\mathrm{nichols}(\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d})$：绘制系统的一组$\mathrm{Nichols}$图，是针对连续状态空间系统$[\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d}]$的每个输入的$\mathrm{Nichols}$图，其中频率范围由函数自动选取，且在响应快速变化的位置会自动采用更多采样点；
   • $\mathrm{nichols}(\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d},\mathrm{iu})$：绘制从系统第$\mathrm{iu}$个输入到所属有输出的$\mathrm{Nichols}$图；
   • $\mathrm{nichols}(\mathrm{num},\mathrm{den})$：绘制以连续时间多项式传递函数表示的系统$\mathrm{Nichols}$图；
   • $\mathrm{nichols}(\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d},\mathrm{iu},\mathrm{w})$或$\mathrm{nichols}(\mathrm{num},\mathrm{den},\mathrm{w})$：利用指定的角频率矢量绘制系统的$\mathrm{Nichols}$图；
4. 绘制等$\mathrm{M}$圆和等$\mathrm{N}$圆的函数$\mathrm{ngrid}()$。

2.1.2 实战1

实验要求：已知一阶环节传递函数为：$G(s)=\frac{5}{3s+1}$，绘制该环节的$\mathrm{Nyquist}$图。

解：

% 实例Chapter7.2.1.2
clc;clear;

% 建立传递函数模型
num=5;den=[3,1];
G=tf(num,den);

% 绘制图形
nyquist(G);grid;
axis([-1,8,-2.5,2.5]);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统Nyquist图','FontSize',15);

2.1.3 实战2

实验要求：已知二阶环节传递函数为：$G(s)=\frac{{\omega }_n^2}{s^2+2\zeta {\omega }_{n}s+{\omega }_n^2}$，其中：${\omega }_n=0.7$，分别绘制$\zeta =0.1、0.4、1.0、1.6、2.0$时的$\mathrm{Bode}$图。

解：

% 实例Chapter7.2.1.3
clc;clear;

% 参数定义
w=[0,logspace(-2,2,200)];
wn=0.7;zeta=[0.1,0.4,1.0,1.6,2.0];

% 建立传递函数模型
for j=1:5
    sys=tf([wn*wn],[1,2*zeta(j)*wn,wn*wn]);
    bode(sys);hold on;
end

set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
legend('ζ=0.1','ζ=0.4','ζ=1.0','ζ=1.6','ζ=2.0');
title('控制系统不同ζ下的Bode图','FontSize',15);

2.1.4 实战3

实验要求：已知一高阶系统的传递函数为：$G(s)=\frac{0.0001s^3+0.028s^2+1.06356s+9.6}{0.0006s^3+0.0286s^2+0.06356s+6}$，绘制系统的$\mathrm{Nichols}$图。

解：

% 实例Chapter7.2.1.4
clc;clear;

% 建立模型
num=[0.0001,0.0281,1.06356,9.6];
den=[0.0006,0.0286,0.06356,6];
G=tf(num,den);

% 绘制等M圆和等N圆，Nichols图
ngrid('new');
nichols(G);

axis([-360,0,-20,40]);
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('等M圆、等N圆及Nichols图','FontSize',15);

2.2 频率法的稳定性分析

2.2.1 对数频率稳定判据

复平面${\Gamma }_{GH}$曲线一般由两部分组成：开环幅相曲线和开环系统存在积分环节和等幅振荡环节时所补作的半径为无穷大的虚圆弧；

半对数坐标下的${\Gamma }_{GH}$曲线确定穿越次数$N$或$N_{+}$和$N_{-}$；$N$的确定取决于$A(\omega )>1$时${\Gamma }_{GH}$穿越负实轴的次数；

1. 穿越点确定

   设$\omega ={\omega }_c$时，
   $$\{\begin{array}{ll}& A({\omega }_c)=|G(\mathrm{j}{\omega }_c)H(\mathrm{j}{\omega }_c)|=1\\ & \\ & L({\omega }_c)=20\lg A({\omega }_c)=0\end{array}$$
   其中：${\omega }_c$称为截止频率；

   对于复平面的负实轴和开环对数相频特性，当取频率为穿越频率${\omega }_x$时，
   $$\phi ({\omega }_x)=(2k+1)\pi ，k=0，\pm 1，\pm 2，\dots$$
   设半对数坐标下${\Gamma }_{GH}$的对数幅频特性曲线和对数相频特性曲线分别为：${\Gamma }_L$和${\Gamma }_{\phi }$，由于${\Gamma }_L$等于$L(\omega )$曲线，则${\Gamma }_{GH}$在$A(\omega )>1$时，穿越负实轴的点等于${\Gamma }_{GH}$在半对数坐标下，对数幅频特性$L(\omega )>0$时对数相频特性曲线${\Gamma }_{\phi }$与$(2k+1)\pi (k=0,\pm 1,\pm 2,\dots )$，平行线的交点；

2. ${\Gamma }_{\phi }$确定

   • 开环系统无虚轴上极点时，${\Gamma }_{\phi }$等于$\phi (\omega )$曲线；
   • 开环系统存在积分环节$\frac{1}{s^{\nu }}(\nu >0)$时，复数平面的${\Gamma }_{GH}$曲线，需从$\omega =0_{+}$的开环幅相特性曲线的对应点$G(\mathrm{j}{0}_{+})H(\mathrm{j}{0}_{+})$起，逆时针补作$\nu \times 90°$半径为无穷大的虚圆弧；需从对数相频特性曲线$\omega$较小且$L(\omega )>0$的点处向上补作$\nu \times 90°$的虚直线，$\phi (\omega )$曲线和补作的虚直线构成${\Gamma }_{\phi }$；
   • 开环系统存在等幅振荡环节$\frac{1}{(s^2+{\omega }_n^2{)}^{{\nu }_1}}({\nu }_1>0)$时，复数平面的${\Gamma }_{GH}$曲线，需从$\omega ={\omega }_{n^{-}}$的开环幅相曲线的对应点$G(\mathrm{j}{\omega }_{n^{-}})H(\mathrm{j}{\omega }_{n^{-}})$起，顺时针补作${\nu }_1\times 180°$半径为无穷大的虚圆弧至$\omega ={\omega }_{n^{+}}$的对应点$G(\mathrm{j}{\omega }_{n^{+}})H(\mathrm{j}{\omega }_{n^{+}})$处；相应地，需从对数相频特性曲线$\phi ({\omega }_{n^{-}})$起，向上补作${\nu }_1\times 180°$的虚直线至$\phi ({\omega }_{n^{+}})$处，$\phi (\omega )$曲线和补作的虚直线构成${\Gamma }_{\phi }$；

3. 穿越次数的计算

   1. 正穿越一次：${\Gamma }_{GH}$由上向下穿越$(-1,\mathrm{j}0)$点左侧的负实轴一次，等价于在$L(\omega )>0$时，${\Gamma }_{\phi }$由下向上穿越$(2k+1)\pi$线一次；
   2. 负穿越一次：${\Gamma }_{GH}$由下向上穿越$(-1,\mathrm{j}0)$点左侧的负实轴一次，等价于在$L(\omega )>0$时，${\Gamma }_{\phi }$由上向下穿越$(2k+1)\pi$线一次；
   3. 正穿越半次：${\Gamma }_{GH}$由上向下止于或由上向下起于$(-1,\mathrm{j}0)$点左侧的负实轴，等价于在$L(\omega )>0$时，${\Gamma }_{\phi }$由下向上止于或由下向上起于$(2k+1)\pi$线；
   4. 负穿越半次：${\Gamma }_{GH}$由下向上止于或由下向上起于$(-1,\mathrm{j}0)$点左侧的负实轴，等价于在$L(\omega )>0$时，${\Gamma }_{\phi }$由上向下止于或由上向下起于$(2k+1)\pi$线；
   5. 补作的虚直线所产生的穿越皆为负穿越；

4. 对数频率稳定判据

   设$P$为开环系统正实部的极点数，反馈控制系统稳定的充要条件是$\phi ({\omega }_c)\ne (2k+1)\pi (k=0,1,2,\dots )$和$L(\omega )>0$时，${\Gamma }_{\phi }$曲线穿越$(2k+1)\pi$线的次数：$N=N_{+}-N_{-}$满足：$Z=P-2N=0$；

5. 奈奎斯特稳定判据

   奈氏判据：反馈控制系统稳定的充要条件是半闭合曲线${\Gamma }_{GH}$不穿过$(-1,\mathrm{j}0)$点，且逆时针包围临界点$(-1,\mathrm{j}0)$点的圈数$R$等于开环传递函数的正实部极点数$P$；

2.2.2 稳定裕度

1. 相角裕度。

   设${\omega }_c$为系统的截止频率，则：
   $$A({\omega }_c)=|G(\mathrm{j}{\omega }_c)H(\mathrm{j}{\omega }_c)|=1$$
   定义相角裕度：
   $$\gamma =180°+\angle [G(\mathrm{j}{\omega }_c)H(\mathrm{j}{\omega }_c)]$$
   相角裕度$\gamma$含义：对于闭环稳定系统，如果系统开环相频特性再滞后$\gamma$度，则系统将处于临界稳定状态；

2. 幅值裕度。

   设${\omega }_x$为系统的穿越频率，则系统在${\omega }_x$处的相角：
   $$\phi ({\omega }_x)=\angle [G(\mathrm{j}{\omega }_x)H(\mathrm{j}{\omega }_x)]=(2k+1)\pi ；k=0,\pm 1,\dots$$
   定义幅值裕度为：
   $$h=\frac{1}{|G(\mathrm{j}{\omega }_x)H(\mathrm{j}{\omega }_x)|}$$
   幅值裕度$h$的含义：对于闭环稳定系统，如果系统开环幅频特性再增大$h$倍，则系统将处于临界稳定状态；

3. 关于相角裕度和幅值裕度的说明。

   • 控制系统的相角裕度和幅值裕度是系统的极坐标图对$(-1+\mathrm{j}0)$点靠近程度的度量，这两个裕量可以用来作为设计准则；
   • 只用幅值裕度或相角裕度，都不足以说明系统的相对稳定性；为了确定系统的相对稳定性，必须同时给出这两个量；
   • 对于最小相位系统，只有当相角裕度和幅值裕度都是正值时，系统才是稳定的，负的裕度表示系统不稳定；
   • 适当的相角裕度和幅值裕度可以防止系统中元件老化造成的影响，且指明了频率值；为了得到满意的性能，相角裕度应当为$30°～60°$，幅值裕度应当大于$6\mathrm{dB}$；
   • 要求相角裕度在$30°$和$60°$之间，即在伯德图中，对数幅值曲线在截止频率处的斜率应大于$-40dB/dec$；在大多数实际情况中，为了保证系统稳定，要求截止频率处的斜率为$-20dB/dec$；如果截止频率上的斜率为$-40dB/dec$，则系统可能是稳定的，可能是不稳定的；如果在截止频率处的斜率为$-60dB/dec$，或更陡，则系统多半是不稳定的；

2.3 MATLAB在稳定性分析中的应用

2.3.1 计算系统稳定裕度的函数margin()

• $\mathrm{margin}(\mathrm{mag},\mathrm{phase},\mathrm{w})$：由$\mathrm{Bode}$指令得到的幅值$\mathrm{mag}$、相角$\mathrm{phase}$及角频率$\mathrm{w}$矢量绘制带有裕度及相应频率显示的$\mathrm{Bode}$图；
• $\mathrm{margin}(\mathrm{num},\mathrm{den})$：计算连续系统传递函数表示的幅值裕度和相角裕度，并绘制相应的$\mathrm{Bode}$图；
• $\mathrm{margin}(\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d})$：计算连续状态空间系统表示的幅值裕度和相角裕度，并绘制相应的$\mathrm{Bode}$图；
• $\mathrm{gm},\mathrm{pm},\mathrm{wcg},\mathrm{wcp}=\mathrm{margin}(\mathrm{mag},\mathrm{phase},\mathrm{w})$：由幅值$\mathrm{mag}$、相角$\mathrm{phase}$及角频率$\mathrm{w}$矢量，计算系统幅值裕度、相角裕度及相应的相角交界频率$\mathrm{wcg}$、截止频率$\mathrm{wcp}$，不直接绘制$\mathrm{Bode}$图；
• $\mathrm{s}=\mathrm{allmargin}(\mathrm{sys})$：计算幅值裕度、相角裕度及对应的频率；

2.3.2 实战1

实验要求：已知一高阶系统的开环传递函数为：$G(s)=\frac{5(0.0167s+1)}{s(0.03s+1)(0.0025s+1)(0.001s+1)}$，计算系统的相角裕度和幅值裕度，绘制系统的$\mathrm{Bode}$图。

解：

% 实例Chapter7.2.3.2
clc;clear;

% 建立模型
num=5*[0.0167,1];
den=conv(conv([1,0],[0.03,1]),conv([0.0025,1],[0.001,1]));
G=tf(num,den);

% 频率范围
w=logspace(0,4,50);

% 绘制Bode图
bode(G,w);grid;
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统Bode图','FontSize',15);

% 计算幅值裕度和相角裕度
[Gm,Pm,Wcg,Wcp]=margin(G)

% 幅值裕度和相角裕度结果
Gm =
  455.2548

Pm =
   85.2751

Wcg =
  602.4232

Wcp =
    4.9620

• 幅值裕度：$K_g=455.2548$，相角裕度：$\gamma =85.2751°$，穿越频率：${\omega }_g=602.4232rad/s$，截止频率：${\omega }_c=4.9620rad/s$。

2.3.3 实战2

实验要求：已知一高阶系统的开环传递函数为：$G(s)=\frac{K(0.0167s+1)}{s(0.03s+1)(0.0025s+1)(0.001s+1)}$，计算当开环增益$K=5、500、800、3000$时，系统稳定裕度的变化。

解：

% 实例Chapter7.2.3.3
clc;clear;

K=[5,500,800,3000];
den=conv(conv([1,0],[0.03,1]),conv([0.0025,1],[0.001,1]));

for j=1:4
    num=K(j)*[0.0167,1];
    G=tf(num,den);
    y(j)=allmargin(G);               % 计算幅值裕度、相角裕度、频率
end

y(1),y(2),y(3),y(4)

% 计算结果：

% K=5相角裕度、幅值裕度、对应频率
ans = 
  包含以下字段的 struct:
     GainMargin: 455.2548
    GMFrequency: 602.4232
    PhaseMargin: 85.2751
    PMFrequency: 4.9620
    DelayMargin: 0.2999
    DMFrequency: 4.9620
         Stable: 1

% K=500相角裕度、幅值裕度、对应频率
ans = 
  包含以下字段的 struct:
     GainMargin: 4.5525
    GMFrequency: 602.4232
    PhaseMargin: 39.7483
    PMFrequency: 237.7216
    DelayMargin: 0.0029
    DMFrequency: 237.7216
         Stable: 1

% K=800相角裕度、幅值裕度、对应频率
ans = 
  包含以下字段的 struct:
     GainMargin: 2.8453
    GMFrequency: 602.4232
    PhaseMargin: 27.7092
    PMFrequency: 329.9063
    DelayMargin: 0.0015
    DMFrequency: 329.9063
         Stable: 1

% K=3000相角裕度、幅值裕度、对应频率
ans = 
  包含以下字段的 struct:
     GainMargin: 0.7588
    GMFrequency: 602.4232
    PhaseMargin: -6.7355
    PMFrequency: 690.5172
    DelayMargin: 0.0089
    DMFrequency: 690.5172
         Stable: 0

2.3.4 实战3

实验要求：已知系统的开环传递函数为：$G(s)=\frac{100(s+5)}{(s-2)(s+8)(s+20)}$，绘制该系统的极坐标图，利用$\mathrm{Nyquist}$稳定判据判断闭环系统的稳定性。

解：

% 实例Chapter7.2.3.4
clc;clear;

% 建立模型
k=100;z=[-5];p=[2,-8,-20];
sys=zpk(z,p,k);

nyquist(sys);grid;
axis equal;
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统Nyquist图','FontSize',15);

• 由开环传递函数可知，开环系统有一个$\mathrm{s}$右半平面的极点，即$p=1$，因此开环系统是不稳定的；
• 开环系统的$\mathrm{Nyquist}$图逆时针包围$(-1,\mathrm{j}0)$点$1$次，根据$\mathrm{Nyquist}$稳定判据可知，闭环系统稳定；

2.3.5 实战4

实验要求：已知系统的开环传递函数为：$G(s)=\frac{100K}{s(s+5)(s+10)}$，分别绘制$K=1、7.8、20$时，系统的极坐标图，利用$\mathrm{Nyquist}$稳定判据判断闭环系统的稳定性。

解：

% 实例Chapter7.2.3.5
clc;clear;

% 建立控制系统模型
z=[];p=[0,-5,-10];k=100.*[1,7.8,20];

G=zpk(z,p,k(1));[re1,im1]=nyquist(G);
G=zpk(z,p,k(2));[re2,im2]=nyquist(G);
G=zpk(z,p,k(3));[re3,im3]=nyquist(G);

% 绘制图形
plot(re1(:),im1(:),'r-',re2(:),im2(:),'b-',re3(:),im3(:),'g-');
legend('K=1','K=7.8','K=20');
axis([-5,1,-5,1]);grid;
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('不同K值控制系统Nyquist图','FontSize',15);

• 由控制系统开环传递函数可知，开环系统没有$s$右半平面的极点，即$p=0$，因此，开环系统稳定；
• 由$\mathrm{Nyquist}$图可知，当$K=1$时，开环系统的$\mathrm{Nyquist}$图不包围$(-1,\mathrm{j}0)$点，根据$\mathrm{Nyquist}$稳定判据，系统闭环稳定；当$K=7.8$和$K=20$时，开环系统的$\mathrm{Nyquist}$图包围$(-1,\mathrm{j}0)$点，根据$\mathrm{Nyquist}$稳定判据，系统闭环不稳定；

2.4 综合实例及MATLAB/SIMULINK应用

实验要求：已知晶闸管-直流电机开环系统的结构图如下图所示，用$\mathrm{SIMULINK}$动态结构图进行频域分析并求频域性能指标。

解：

【$STEP1$】：$\mathrm{SIMULINK}$中建立控制系统动态模型，文件名为：$experiment7\_2\_4.mdl$。

【$STEP2$】：求系统线性状态空间模型，并求频域性能指标。

>> [A,B,C,D]=linmod('experiment7_2_4');
>> sys=ss(A,B,C,D);
>> margin(sys);

• 控制系统幅值裕度：$\mathrm{GM}=26.4dB$，频率：$152rad/s$；
• 控制系统相位裕度：$\mathrm{PM}=54deg$，频率：$25.5rad/s$；