控制教程 —— 介绍篇：5.频域方法控制器设计(粗略版)

采用频域响应方法来设计控制器可能不如之前研究的其他方法直观。但是，它具有某些优势，尤其是在实际情况下，例如根据实测数据来对系统进行建模的时候。在本教程中，我们将看到如何使用系统的开环频率响应来预测其闭环时间响应。
本教程中使用的主要MATLAB命令为：

增益和相位裕度

在介绍篇：系统分析博文中已经了解，系统的频率响应是指在正弦输入下，评估系统输出对于正弦输入的缩放和平移的表现。正弦输出的缩放和平移随频率变化的方式提供了有关系统时间响应的有用信息，特别地，将系统的频率响应用于确定系统的“鲁棒性”。例如，系统离不稳定有多近？在这里，我们使用增益裕度和相位裕度这两个量来表示系统在不稳定之前具有的裕度。
考虑以下通用的反馈系统：

其中 $K$ 是可变增益，而 $G(s)$ 是被控对象。增益裕度定义为使闭环系统不稳定所需的开环增益的变化。增益裕度更高的系统在闭环不稳定之前，可以承受较大的系统参数变化。
相位裕度定义为使闭环系统不稳定所需的开环相移变化。相位裕度还可以测量系统对时间延迟的容忍度。如果环路中的时间延迟大于 $PM/{\omega }_{gc}$ (其中 ${\omega }_{gc}$ 是 rad/sec为单位，$PM$ 是增益为0时对应的相位欲度，转换成弧度)，则闭环系统将变得不稳定。时间延迟 ${\tau }_d$ 可以看作是框图前馈路径中的一个额外模块，它会给系统增加相位滞后，但对增益没有影响，也就是说，时间延迟可以表示为大小为1且相位为 $-\omega {\tau }_d$ 的模块。
现在，我们将不关心这些延迟来源于哪里，而仅考虑伯德图上的增益和相位裕度。
闭环系统的相位裕度是开环系统的相位在开环系统幅度为0dB的频率下达到 $-180°$ 所需的额外相位滞后量(增益穿越频率 ${\omega }_{gc}$)。同样，增益裕度是开环系统的相位等于 $-180°$ 所需的额外增益裕量(通常以dB为单位)(相位穿越频率 ${\omega }_{pc}$)。

关于相位裕量的一个优势时，在更改环路增益时，无需重新绘制伯德图就可以找到新的相位裕量。如果您还记得，增加增益只会使幅度向上(或向下)移动。这等效与更改幅度图上的y轴。找到相位裕度仅是找到新的增益交叉频率并读取相位裕度的问题。例如，假设您使用以下命令来生成如下所示的伯德图：

s = tf('s');
G = 50/(s^3 + 9*s^2 + 30*s +40);

bode(G)
grid on
title('Bode Plot with No Gain')

您可以看到相位裕度约为100度，现在，假设您通过输入命令bode(100*G)增加了100的增益。您可以得到以下：

bode(100*G)
grid on
title('Bode Plot with Gain = 100')

如图所示，相位图和之前完全一样，并且幅度图向上移动了40dB(增益为100)。现在相位裕度约为-60度。如果幅度图的y轴向下移动40dB，则可以得到相同的结果，这样操作的话，查看第一张伯德图，找到曲线与-40dB相交的位置，然后读出相位裕度，它应约为-60度，与第二个伯德图相同。
我们可以让MATLAB使用margin(G)命令计算并显示增益和相位裕度，该命令返回增益和相位裕度，增益和相位穿越频率，以及这些量在伯德图上的位置，参考以下示例：

margin(100*G)

带宽频率

带宽频率定义为闭环响应的幅度相比于直流幅度(频率接近0的幅度)下降3dB时对应的频率。通常，我们在对频率响应进行系统设计时，我们是希望根据开环响应来预测其闭环响应，这里我们将使用二阶系统来等效闭环。有些教科书上描述带宽等于开环响应幅度在-6到-7.5dB之间的频率，相位在-135°到-225°之间。
为了说明带宽频率的重要性，这里将显示输出如何随不同的输入频率而变化。可以看到，频率小于 ${\omega }_{bw}$（带宽频率）的正弦输入会被系统良好地跟踪。频率大于 ${\omega }_{bw}$ 的正弦输入被衰减(幅度)达到0.707或更大(并且也发生了相移)。
假设我们有以下表示系统的闭环传递函数：

G = 1/(s^2 + 0.5*s + 1);
bode(G)

由于这是闭环传递函数，因此我们的带宽频率将是增益对应于-3dB的频率，从图中可以得到它约为1.4rad/s，我们还可以从图中得知，对于0.3rad的输入频率，输出正弦波的幅值应约为1，相位偏移几度(在输入之后)。对于3rad/s的输入频率，输出幅度应约为-20dB(或输入的1/10)，相位应接近-180(几乎完全反相)。我们可以使用lsim命令来模拟系统对正弦输入的响应。
首先，考虑一个频率低于 ${\omega }_{bw}$ 的正弦输入，这里还要好注意，我们关心的是稳态响应，因此，我们调整横坐标来让稳态响应更加直观(忽略瞬时相应)。

G = 1/(s^2 + 0.5*s + 1);
w = 0.3;
t = 0:0.1:100;
u = sin(w*t);
[y,t] = lsim(G,u,t);
plot(t,y,t,u)
axis([50 100 -2 2])

注意，输出(蓝色)很好地跟踪了输入(红色)。它可能比输入延迟约几度，但是如果将输入的频率设置为高于系统的带宽频率，则会得到非常失真的响应(相对于输入)：

G = 1/(s^2 + 0.5*s + 1);
w = 3;
t = 0:0.1:100;
u = sin(w*t);
[y,t] = lsim(G,u,t);
plot(t,y,t,u)
axis([90 100 -1 1])

同样，正如预期的那样，幅度大约是输入幅度的1/10，并且几乎和输入相位完全相反(相差180°)。

奈奎斯特图

奈奎斯特图使我们能够通过观察开环表现来预测闭环系统的稳定性和性能。无论开环稳定性如何，都可以将Nyquist准则用于设计目的(而伯德图设计方法假设系统在开环中是稳定的)。因此，当伯德图显示令人困惑的信息时，我们使用Nyquist准则来确定系统闭环稳定性。
注意：对于虚轴上具有开环极点或者零点的系统，MATLAB nyquist命令不能提供良好的显示，因此，建议您使用nyquist1.m，该文件可以正确处理虚轴上的极点和零点。
那奎斯特图基本上是 $G(j\omega )$ 的图，其中 $G(s)$ 是开环传递函数， $\omega$ 是频率，包围整个右半平面。在绘制奈奎斯特图时，要同时考虑正频率(从0到无穷大)和负频率(从负无穷大到0)，我们将以红色表示正频率，以绿色表示负频率。绘制奈奎斯特图时使用的频率通常如下所示(如果您可以想象该图延伸到无穷远)：

但是，如果我们在虚轴上有开环极点或者零点，则不会再这些点处定义 $G(s)$ ，并且在绘制轮廓时必须在他们周围画圈。轮廓示例如下所示：

请注意，轮廓围绕虚轴上的极点画圈，如前所述，MATLAB nyquist命令没有考虑虚轴上的极点或零点，因此会产生错误的绘图。要更正该问题，请下载并使用nyquist1.m，如果在虚轴上有一个极点，则必须使用nyquist1。如果虚轴上没有极点或零点，或者我们有零极点对消，则可以使用nyquist命令或者nyquist1命令

在研究反馈控制时，我们对$G(s)$并不那么感兴趣，而是更关注闭环传递函数。
如果$1+G(s)$过原点，则$G(s)$将过点-1，因此，奈奎斯特图在实轴上-1点附近的表现非常重要。但是，标准奈奎斯特图上的轴可能很难看清这一点，若要更正该问题，可以使用lnyquist.m函数。
使用MATLAB查看简单的Nyquist图，我们将定义一下传递：

s = tf('s');
G = 0.5/(s - 0.5);
nyquist(G)
axis([-1 0 -1 1])

现在我们来看一下以下传递函数的奈奎斯特图：

请注意，该传递函数在原点处有一个极点，我们将看到使用nyquist和nyquist1命令之间的区别。

G = (s + 2)/(s^2);
nyquist(G)

nyquist1(G)

lnyquist(G)

注意，标准奈奎斯特图不正确，特别地，对于半径接近无穷大($s$接近0)的情况，它无法捕获奈奎斯特图的行为。nyquist1图是正确的，这使我们可以评估-1点的周围并应用Nyquist准则，但是很难看到在-1点附近发生了什么，lnyquist函数采用修改后的比例，该比例在-1点附近保持绘图的正确性，并保持正确的形状。

伯德图的闭环性能

为了从开环频率响应预测闭环性能，我们需要做一些假设：

• 如果我们要使用伯德图进行设计，则系统必须在开环下稳定。
• 对于典型的二阶系统，如果相位裕度在0到60度之间，则闭环阻尼比大约等于相位裕度除以100，如果相位裕度大于60度，我们需要谨慎使用该结论。
• 对于典型的二阶系统，阻尼比，带宽频率和上升时间之间的关系可由附录部分的公式。
• 您可以使用一个非常粗略的估计是带宽大约等于自然频率。
  让我们使用这些概念为以下系统设计控制器：
  
  其中 $C(s)$ 是控制器， $P(s)$是：
  
  设计必须符合以下规范：
• 零稳态误差
• 最大过冲必须小于40%
• 上升时间必须小于2s

解决此问题的方法有两种：一种是画图，另一种是计算，在MATLAB中，画图是比较好的，所有我们采用这个。首先让我们看一下伯德图。

P = 10/(1.25*s + 1);
bode(P)
margin(P)

可以从该伯德图中直接读取几个特性，首先，我们可以看到带宽频率约为10rad/s。由于带宽频率与自然频率大致相同(对于此类一阶系统)，因此上升时间约为$1.8/BW=1.8/10=0.18$秒。这是一个粗略的估计，我们可以说上升时间约为0.2s。
该系统的相位裕度约为95度，阻尼比$=PM/100$的关系仅在$PM<60$时成立(谨慎)。由于系统是一阶的，因此不会有过冲现象。
最后一个需要注意的是稳态误差。稳态误差也可以直接从伯德图上读取。从低频渐近线与 $\omega =1rad/s$ 线的交点可以找到误差常数($K_p$，$K_v$，$K_a$)。只需将低频线扩展到 $\omega =1$ 即可。此时的大小是误差常数，由于该系统的伯德图低频(频率=0)处是水平线，因此我们知道该系统的单位反馈为零型系统。因此，交点很容找到，增益为20dB(幅度10)。这意味着误差函数的常数为10。稳态误差为$1/(1+K_p)=1/(1+10)=0.091$。
如果我们的系统是1型而不是0型，则可以通过类似的方法找到稳态误差常数

让我们通过查看阶跃响应图来检查我们的预测

sys_cl = feedback(P,1);
step(sys_cl)
title('Closed-Loop Step Response, No Controller')

如图所示，我们的预测非常好。该系统的上升时间约为0.2s，没有超调，并且稳态误差约为9%。现在我们需要选择一个控制器，使我们能够满足设计标准。我们选择PI控制器是因为对于阶跃输入它将获得零稳态误差。同样，PI控制器有一个零点，我们可以调整它，这为我们提供了额外的设计灵活性，回想一下，PI控制器具有如下形式：

我们需要找到的第一件事就是阻尼比，对应于40%的超调。将这个值代入与过冲和阻尼比相关的方程式中，控制教程 —— 介绍篇：2.系统分析。
$$\zeta =\frac{-\ln (Mp)}{\sqrt{{\pi }^2+{\ln }^2(Mp)}}$$
我们发现与该过冲对应的阻尼比约为0.28。因此，我们的相位裕度应该至少为30度。如果我们希望上升时间小于1.75秒，则带宽频率必须大于或等于12。
既然我们知道了所需的相位裕度和带宽频率，就可以开始设计了。记住，我们在看开环伯德图，因此，我们的带宽频率将对应大约-7dB增益的频率
让我们看看PI控制器的积分器如何影响我们的响应。
这里仅添加一个积分项

P = 10/(1.25*s + 1);
C = (s + 1)/s;
bode(C*P, logspace(0,2))

我们的相位裕度和频率太小，我们将增加一个零点，并将零点放置在-1处，再来看看发生了什么：

P = 10/(1.25*s + 1);
C = 5*((s + 1)/s);
bode(C*P, logspace(0,2))

可以看到新增零点为-1时，我们的相位裕度大于60度(甚至比预期的过冲更少)，我们的带宽频率约为11rad/s，这时将尝试提高贷款频率，而又不会产生更多的相位裕量，让我们尝试将增益增加5，看看会发生什么。

P = 10/(1.25*s + 1);
C = 5*((s + 1)/s);
bode(C*P, logspace(0,2))

看起来变好了，让我们看看阶跃响应：

如图所示，闭环响应比我们期望的要好，但是这种方法需要我们尝试使用极点和/或零点来调整系统的增益和相位。

奈奎斯特图的闭环稳定性

考虑负反馈系统：

请记住，根据柯西幅角定理， $G(s)H(s)$ 的图包围 $-1+j0$ 的次数(N) 等于 $1+G(s)H(s)$ 的零点个数(Z)，减去 $1+G(s)H(s)$ 的极点个数§，（N = P - Z）。

• $1+G(s)H(s)$的零点是闭环传递函数的极点
• $1+G(s)H(s)$的极点是开环传递函数的极点

根据以上定义：

• P = $G(s)H(s)$ 的开环极点。
• N = 奈奎斯特图绕 -1 的次数。
• -1的顺时针包围圈算作正包围圈。
• -1的逆时针包围圈算作负包围圈。
• Z = 闭环系统的右半平面(正，实)极点。

与这三个量相关的奈奎斯特未定性判据为：

注意：这只是奈奎斯特准则的一种约束。另一个约束指出，正数N表示逆时针或逆时针包围圈 -1 ，这种情况下方程式变成Z = P -N。在这个教程中，我们将顺时针包围使用正。
知道开环的右半平面(不稳定)极点数§，以及由奈奎斯特图(N)得出的 -1 的包围圈数，我们可以确定系统的闭环稳定性。如果 Z = P + N为非0正数，则闭环系统不稳定。
我们也可以使用Nyquist图来找到使闭环单位反馈系统稳定的增益范围，我们将测试如下系统：


该系统具有增益 $K$，可以改变该增益来修改闭环系统的响应。但是，由于必须确保闭环稳定，我们将看到只能在特定范围内修改该增益。这就是我们要寻找的增益范围。
我们需要做的第一件事是在开环传递函数中找到正实极点数。

roots([1 -8 15])

开环传递函数的两个极点均为正。因此，我们需要Nyquist图的两个逆时针(N = -2) 包围，以便具有稳定的闭环系统。如果包围数量少于两个或包围数量不是逆时针，我们的系统将不稳定。
让我们看一下奈奎斯特图，增益为1：

有两个逆时针 -1 。因此，系统在增益为 1 时使稳定的。现在，我们看到将增益增加到20时系统的表现：

nyquist(20*G)

从图中可以看出，圆圈是增大的，说明无论增加多少增益，闭环系统都是稳定的。但是，如果我们减小增益，该图将缩小，系统可能变得不稳定。

nyquist(0.5*G)

系统现在不稳定，通过反复试验，我们发现该系统在增益小于0.80时将变得不稳定，我们可以通过放大奈奎斯特图以及查看闭环阶跃响应(增益为0.79,0.80,0.81)来验证结果。
增益裕度
我们已经将增益裕度定义为分贝(dB)，表示开环增益的变化，在-180度相移下要求开环增益等于0dB。现在，我们将找出来源。首先，假设如果没有 -1 的奈奎斯特包围，我们的系统将是稳定的，例如：

可以发现右半平面没有开环极点，因此，如果没有 -1 的奈奎斯特环绕，则右半平面中没有闭环极点，现在，在闭环系统变得不稳定之前，我们可以改变增益多少，让我们看下面的图：

如果增益增加到一定范围内，则该图表示的开环系统将在闭环中变得不稳定。负实轴区域介于 -1/a (定义为开环具有-180度相位的点，即该图与实轴交叉的点)与 -1 之间的负值，代表在闭环不稳定之前可以容忍的增益量。
可以意识到如果增加增益a，该图将达到 -1：
$$G(j\omega )=\frac{-1}{a}$$
或者：
$$aG(j\omega )=-1$$因此，我们说增益裕度是一个单位，但是，我们之前提到，增益裕度通常以分贝为单位，因此，增益裕度为：
$$GM=20{\log }_{10}(a)\mathrm{d}\mathrm{B}$$现在，我们将找到之前看过的稳定开环函数的增益裕度：
$$G(s)=\frac{50}{s^3+9s^2+30s+40}$$并且可以通过输入一下命令来看奈奎斯特图

G = 50/(s^3 + 9*s^2 + 30*s + 20);
nyquist(G)

如前所述，找到增益裕度所需要的就是找到a，为此，我们需要找到恰好有-180度相移的点。着意味着此时的传递函数为实数(没有虚部)。分子已经是实数，所以我们只需要看分母，当 $s=j\omega$时，分母中唯一具有虚部的项是s的奇次幂。因此，要想 $G(jw)$ 都是实数，则必须有：
$$-j{\omega }^3+30j\omega =0$$着意味着 $\omega =0$ (这是Nyquist图中最右边额点) 或 $\omega =\sqrt{30}$，然后，我们可以使用polyval找到 $G(j\omega )$的值：

w = sqrt(30);
polyval(50,j*w)/polyval([1 9 30 40],j*w)

结果是 -0.2174 + 0i，虚部为0，可以看出答案是正确的。我们还可以通过再次查看奈奎斯特图来验证这一点。
我们发现 -180度相移发生在 -0.2174 + 0i。该点以前定义为 - 1/a 。 可以得出：
$$\frac{-1}{a}=-0.2174$$$$\Rightarrow a=4.6$$$$\Rightarrow GM=20{\log }_{10}(4.6)=13.26\mathrm{d}\mathrm{B}$$
现在，我们有了增益裕度。让我们看一下使用a=4.6的增益并放大奈奎斯特图来看准确性如何。

a = 4.6;
nyquist(a*G)

相位裕度
我们已经讨论了相位裕度的重要性。因此，我们只会谈论这个概念的来源，我们将相位裕度定义为单位增益所需的开环相移变化，以使闭环系统不稳定，让我们看一下此概念用图表示，

让我回顾之前的分析，如果奈奎斯特图围绕 -1，则该系统在闭环中将不稳定。但是，我们还必须意识到，如果将图移动 $\theta$ ，它将在负实轴上碰到 -1。因此，使该系统在闭环中临界稳定所需要的角度称为相位裕度(以度为单位)。为了找到测量该角度的点，我们绘制了一个半径为 1 的圆，在奈奎斯特图中找到幅度为1(增益为0dB)的点，并测量该点到-180度的相移的角度。

附录

下面是函数wbw.m，给定响应的阻尼比和上升时间或稳定时间，此函数将返回系统的近似带宽频率。要在MATLAB中使用wbw.m函数，请在命令行中输入wbw。首先将提示您输入阻尼比，然后如果您使用稳定时间，则输入”0“，如果您使用稳定时间，则输入”1“。最后，将要求您提供所需的建立时间和峰值时间，然后将返回带宽频率。

   function[] = wbw()
         DR = input('Damping Ratio?    ');
         n = 0;
         n = input('Settling Time(0) or Peak Time(1)?   ');
         if n == 0
         ts = input('Settling Time?   ');
         temp1 = sqrt((4*DR^4) - (4*DR^2) +2);
         temp2 =  1- (2*DR^2);
         temp3 = 4/(ts*DR);
         Wbw = temp3*sqrt(temp1 + temp2)
         end
         if n ==1
         ts = input('Peak Time?   ');
         temp1 = sqrt((4*DR^4) - (4*DR^2) +2);
         temp2 =  1- (2*DR^2);
         temp3 = ts*sqrt( 1- DR^2);
         temp4 = pi/temp3;
         Wbw = temp4*sqrt(temp1 + temp2)
         end