如何入门自动控制理论

一、一个典型的系统响应有什么特点

人，对于外界信息的响应，有两个明显特征，一个是“选择性记忆”，另一个是“永久性拖延”。什么是“选择性记忆”——人们只选择记忆自己感兴趣的、认为有价值的信息，而自动筛选掉那些无关紧要的信息，这实际是给大脑减负，减少能量消耗。什么是“永久性拖延呢”？——这个就很好解释了，就一个字：“懒”。

那对于一个更简单普遍一些的系统（线性时不变系统，简称LTI），它的响应是什么呢？——和人类似，它也有“选择性记忆”和“永久性拖延”，只不过改了个名字：“选择性记忆”叫做“幅值有衰减”；“永久性拖延”叫做“相位有延迟”。

二、幅频响应和相频响应

举个简单的例子——弹簧振子来说明一下。

这是一个单输入（质量块激励力）单输出（质量块位移）系统，质量块可以看成一个刚体（质点），两端分别通过一个弹簧和一个阻尼器连接到壁面上。弹簧提供与位移成正比的力，阻尼器提供与速度成正比的力，质量块在有加速下还会产生惯性力，这是一个典型的二阶系统。数学描述如下：

$m\cdot\frac{d^2r}{dt^2}+b\cdot\frac{dr}{dt}+k\cdot r=kF(t)$

惯性力+阻尼力+弹性力=外界激励

现在假设用一个交变的正弦力 $F(t)=F_0sin(\omega t)$ 去激励这个弹簧振子，会出现什么响应呢呢？——俗话说：“龙生龙，凤生凤，老鼠的儿子会打洞”，试验和理论分析结果都表明：当给LTI系统一个正弦激励时，其响应也是一个正弦，而且频率不变。具体见下图：

图中蓝线表示一个正弦的输入，红线表示系统的响应，稳态时也呈现一个正弦曲线，MATLAB代码如下：

clc;clear;
wc=1;eta=0.7;
w=0.02;
s=tf('s');
H=wc^2/(s^2+2*wc*eta*s+wc^2);
h=figure('color','w');
for k=1:15
    [u,t] = gensig('sin',2*pi/(w*10*k),10*pi/(w*10),0.001/(w*10*k));
    y=lsim(H,u,t);
    plot(t,u,t,y,'linewidth',1);
    title('Made by J Pan')
    xlim([0 10*pi/(w*10)]);xlabel('Time(s)')
    ylim([-1.5 1.5]);ylabel('Amplititude(mm)')
    grid on;
    legend('Input','Output');
    drawnow;
end

对于一个正弦函数，有三个变量：幅值、频率和相位。当频率不变时，能变的只剩下幅值和相位了。其中幅值响应随频率的变化情况就叫做“幅频响应”，相位响应随频率的变化情况就叫做“相频响应”。那怎么进行“幅频响应”和“相频响应”的分析呢？

还是前面的弹簧振子模型，两边都进行拉普拉斯变换： $(ms^2+bs+k)R(s)=kF(s)$

则系统的传递函数为： $Y(s)=\frac{R(s)}{F(s)}=\frac{k}{ms^2+bs+k}=\frac{\omega_n^2}{s^2+2\xi\omega_ns+\omega_n^2}$

其中 $\omega_n=\sqrt{\frac{k}{m}}$ ， $\xi=\frac{b}{2\sqrt{km}}$ ，建议大家都自己动手推导一下，理解了二阶系统的特性，一半的自动控制问题都能解决了。关于什么是拉普拉斯变换，可参照文章：

J Pan：从另一个角度看拉普拉斯变换

知道了传递函数，怎么知道幅频响应和相频响应呢？其实也很简单，令 $s=j\omega$ ，则

$\left| Y(j\omega) \right|=\frac{1}{\sqrt{(1-\frac{\omega}{\omega_n})^2+4\xi^2(\frac{\omega}{\omega_n})^2}}$ ， $\angle Y(j\omega)=-artan\frac{2\xi\frac{\omega}{\omega_n}}{1-(\frac{\omega}{\omega_n})^2}$

把这两个函数画出来就长这样：

自动控制领域里面有一个专门的名词来描述上面这张图，那就是——“伯德图”。伯德图分两部分，上面一部分是幅频响应，横坐标为频率，纵坐标为 $20log(A_{out}/A_{in})$ ，单位为分贝，这本质还是响应幅值比上输入幅值，因为这个比值变化量级比较大，所以一般采取对数的表达方式。伯德图的下面一部分是相频响应，横坐标是频率，纵坐标是相位延迟。

仔细观察上面的伯德图可以发现，当频率小于 $1Hz$ 时，幅值响应增益基本为 $1$ ，也就是幅值基本和输入一致，相位落后约 $(0\sim90^o)$ 。当频率大于 $1Hz$ 时，幅值响应开始迅速衰减，当频率增加至 $10Hz$ 时，幅值响应为 $-40dB$ ，也就是输入幅值的1%，相位落后接近180°。

可见，对于一个一般的线性时不变系统（LTI），系统具有低通特性，我们把幅频响应增益为0.707，也就是 $20log0.707=-3dB$ 时对应的频率称之为截止频率，低于这个频率的系统能通过，高于这个频率的，系统会有较大幅度的过滤，输出很小。对于一个二阶系统，当阻尼比 $\xi=1/\sqrt{2}$ 时， $\omega_n$ 就代表了系统的截止频率。

举个例子：输入有两个正弦函数组成： $u=sin(0.5\omega_nt)+sin(10\omega_nt)$ ，一个分量是截止频率的1/2，另一个分量是截止频率的10倍。

由上图可以看出， $0.5\omega_c$ 频率的分量能较好的通过，而 $10\omega_ct$ 频率分量则基本被过滤掉了。这是两个频率分量的情况，那假如更复杂一些的输入呢？比如常见的阶跃信号：

$f(t)=\left\{ \begin{array}{rcl} 1 & & {t\geq0}\\ 0 & & t<0\\ \end{array} \right.$

阶跃信号的傅里叶变换长什么样呢？这个函数的傅里叶变换需要一些小技巧，不过这不是关键，我们直接给出结论：

$F(\omega)=\pi\delta(\omega)+\frac{1}{j\omega}$

画出来就是：

阶跃函数在所有频率都有分量，而且随着频率的增加，其幅值越来越小，也就是低频下的分量贡献更多。假如现在输入是阶跃函数，那输出会是什么样？

输出整体来说，基本呈“阶跃”样貌，但是在细节上又有不同。现在，假如改变系统的截止频率：

MATLAB代码如下：

clc;clear;
wc=0.5;eta=0.7;
s=tf('s');
h=figure('color','w');
for k=1:50
    t=0:0.001:10;
    u=ones(length(t),1);
    H=(wc*k)^2/(s^2+2*(wc*k)*eta*s+(wc*k)^2);
    y=lsim(H,u,t);
    plot(t,u,t,y,'linewidth',1);
    title('Made by J Pan')
    xlim([0 10]);xlabel('Time(s)')
    ylim([0 1.5]);ylabel('Amplititude(mm)')
    grid on;
    string=['Output (\omega_c=',num2str(wc*k),'rad/s)'];
    legend('Input (step)',string);
    drawnow;
end

可以看出，随着截止频率的增加，越来越多频率的分量能够通过系统，那输出也就更接近输入信号了。

三、传递函数的零、极点是怎么个回事

还是前面的弹簧振子模型，当输入 $kF(t)$ 为单位阶跃函数时，系统的输出为：

$y(t)=1-\frac{1}{\sqrt{1-\xi^2}}e^{-\xi\omega_nt}sin(\omega_dt+\beta)$ ，其中 $\omega_d=\omega_n\sqrt{1-\xi^2}$ ， $\beta=arcos\xi$ 。

可见，系统的响应分两部分，一部分恒为1，也就是输入，我们称之为稳态响应；另一部分为幅值不断衰减的正弦函数，我们称之为瞬态响应。瞬态响应不影响最终稳态结果，因为它最终会衰减为零，但是它影响瞬态过程，接下来我们研究瞬态响应的影响因素有哪些。

考虑一下传递函数的分母： $s^2+2\xi\omega_ns+\omega_n^2=0$ ，这是一个一元二次方程，解起来很容易，当 $0<\xi<1$ 时，它的解为：

$s_{1,2}=-\xi\omega_n\pm j\omega_n\sqrt{1-\xi^2}=-\xi\omega_n\pm j\omega_d$

在复平面，实轴分量 $\sigma=-\xi\omega_n$ ，虚轴分量为 $\omega_d=\omega_n\sqrt{1-\xi^2}$ ，瞬态响应部分是由系统传递函数分母的根在复平面的分布决定的，因此，一般称传递函数的分母为特征方程。我们知道瞬态响应是系统模态的叠加，换句话说，特征方程的根就代表了系统的“模态”。关于什么是模态及模态分析，参照文章：

J Pan：模态分析——一个找“基”的过程

可以很容易看出：特征方程的根的幅值 $\left| s_1 \right|=\left| s_2 \right|=\omega_n$ ，根与实轴的夹角 $\beta$ 与阻尼比有关。当阻尼比大的时候，实轴分量 $\sigma$ 大，虚轴分量 $\omega_d$ 小，系统的截止频率低，系统响应慢。当阻尼比小的时候，实轴 $\sigma$ 分量小，虚轴分量 $\omega_d$ 大，瞬态衰减慢，系统稳定性差。综合来说，当阻尼比在 $0.65\sim0.8$ 的时候，系统的快速性和稳定性能达到一个较好的平衡，一般推荐阻尼比 $\xi=0.707$ ，也就是 $\beta$ 取45°。

MATLAB代码如下：

clc;clear;
wc=1;eta=0.05;
s=tf('s');
h=figure('color','w');
for k=1:20
    t=0:0.001:10;
    u=ones(length(t),1);
    H=(wc)^2/(s^2+2*(wc)*eta*k*s+(wc)^2);
    y=lsim(H,u,t);
    plot(t,u,t,y,'linewidth',1);
    title('Made by J Pan')
    xlim([0 10]);xlabel('Time(s)')
    ylim([0 2]);ylabel('Amplititude(mm)')
    grid on;
    string=['Output (\xi=',num2str(eta*k)];
    legend('Input (step)',string);
    drawnow;
end

传递函数的极点代表了系统的模态，也就是说，当激励的频率与极点代表的频率接近或相同时，系统将产生共振，系统的响应达到最大值。

那零点代表什么呢？——与极点相反，代表系统能够屏蔽的的模态，也就是反共振点，换句话说，当激励的频率与零点代表的频率相同时，系统将没有任何响应或响应非常小。

综合起来来说：系统瞬态响应取决于极点代表的模态；稳态响应取决于激励频率离零点和极点的距离，离极点越近，稳态响应越大，离零点越近，稳态响应越小。

举个例子：假设一个传递函数是

$G(s)=\frac{s^2+2\xi\omega_zs+\omega_z^2}{s^2+2\xi\omega_ps+\omega_p^2}$ ，其中 $\xi=0.0005,\omega_z=1,\omega_p=5$ ，即 $G(s)=\frac{s^2+0.001s+1}{s^2+0.005s+25}$ ，零点对应的频率为 $1rad/s$ ，极点对应的频率为 $5rad/s$ 。

可见，从幅频响应看，零点对应的稳态响应很小，而极点对应的稳态响应很大。

clc;clear;figure('color','w');
a = [1 0.001 1];b = [1 0.005 25];
w = logspace(-0.1,1,10000);
[h,w]=freqs(a,b,w);
loglog(w,abs(h));xlim([0.8 10]);
xlabel 'Frequency (rad/s)', ylabel 'Magnitude';
grid on

四、什么是闭环控制系统

由前面的分析，我们知道，只要系统的截止频率足够高，就会有足够多频率的分量通过，输出就能很好的跟随输入。然而，不幸的是，有很多系统，其截止频率没有足够高，那怎么办呢？举个例子：

假设一个系统的传递函数是： $G(s)=\frac{1}{s+1}$

这个系统的阶跃响应为：

响应到稳定值的63.2%时所需要的时间约为1s，很慢，原因就是截止频率太低了，其伯德图为

可见，当输入大于 $1rad/s$ 时，系统幅值响应急剧衰减，那我们假如输入是 $10rad/s$ ，幅值衰减到 $-20dB$ ，也就是输入的10%，输出很小！假如我们想输入维持在 $10rad/s$ ，输出又不怎么衰减，怎么办呢？——负反馈系统！

这时候，万能的PID就该上场了！关于PID，有很多不错的帖子从定性方面进行了描述，如：

如何通俗地解释 PID 参数整定？确定有穷自动机：PID控制算法原理（抛弃公式，从本质上真正理解PID控制）

下面我们从定量的方面再来讨论一下。前面分析了，只要截止频率足够高，输出就能很好的跟随输入，接下来我们就可以看看比例环节P是怎么改变截止频率的。

假定 $G(s)=\frac{1}{\tau s+1}$ ，很容易得出： $(R(s)-C(s))\frac{k_p}{\tau s+1}=C(s)$

则整个系统的传递函数变为了： $\frac{C(s)}{R(s)}=\frac{k_p}{\tau s+k_p+1}$ ，与之前的开环传递函数相比： $\frac{C(s)}{R(s)}=\frac{1}{\tau s+1}$ ，截止频率由 $\frac{1}{\tau}$ 变为了 $\frac{k_p+1}{\tau}$ 。假如是令 $\tau=1$ , $k_p=100$ ，则闭环系统的伯德图为：

可见，其截止频率增加到了 $100rad/s$ ，闭环系统的阶跃响应为：

响应到稳定值63.2%的时间缩减到了约0.01s，响应的快速性得到明显的提升。

前面介绍了比例环节 $k_p$ 的影响，那如果加入 $k_i/s$ 积分环节会有什么变化？

$(R(s)-C(s))\frac{k_p+k_i/s}{\tau s+1}=C(s)$

前项通道传递函数为：

$\frac{k_p+k_i/s}{\tau s+1}=\frac{k_ps+k_i}{\tau s^2+s}$

系统的闭环传递函数变为：

$\frac{C(s)}{R(s)}=\frac{k_ps/\tau+k_i/\tau}{s^2+(k_p+1)s/\tau +k_i/\tau}$

这就是一阶系统经过PI控制器校正之后的闭环传递函数，计算这个传递函数单位截止频率还是有一点小复杂的，不适合工程计算。怎么办呢？

在说解决方法之前，先补充一点知识点：前面我们说，只要截止频率足够高，输出就能很好的跟随输入，这种说法不是很严谨，因为文章开头就说了，各个频率分量通过系统时，不仅仅是幅值有衰减，相位也是有延迟了，如果相位延迟过多，信号重新叠加后也有可能和输入差别比较大。因此严谨的说法是：当相位延迟在一定范围内，截止频率越高，输出和输入就越接近。怎么定量说这个事呢？——相位裕度。

前面提到，加入PI控制器后，系统的闭环传递函数变得很复杂了，既有极点，又有零点，不太容易分析，因此，自动控制领域的前辈呢想了一个办法，那就是先设计开环前项通道的截止频率和相位裕度，再来分析闭环系统的响应。

当校正后的系统（控制器与G(s)组成的开环前向通道）的幅值响应开始衰减时（输出幅值小于输入，也就是0dB），信号的相位延迟与-180°的距离就是相位裕度。可见，相位裕度越大越好，理论上如果相位裕度达到180°，那就没有任何相位延迟，当然实际上是达不到的，一般来说相位裕度在40°以上，当然根据不同的应用场合会不一样。总结一下：

系统的幅值响应可以用截止频率来约束，相位响应用相位裕度来约束。

有了这两个约束，我们就可以设计很多控制系统了。

和之前一样，还是假设 $\tau=1$ ，截止频率改为 $75rad/s$ ，相位裕度设为65°，我们可以偷点懒——借助于MATLAB来计算。

打开PI controller模块，选择PI并点击tuning，输入截止频率和相位裕度即可计算出所需要的比例增益和积分增益：

上图左侧为闭环传递函数伯德图，右侧为开环前向通道传递函数伯德图，可见，闭环传递函数截止频率更高，由 $75rad/s$ 增加到了 $100rad/s$ 。比例增益和积分增益分别为 $k_p=67.5$ ， $k_i=2445$ ，将 $k_p$ ， $k_i$ 代入闭环传递函数可得：

$\frac{C(s)}{R(s)}=\frac{67.5s+2445}{s^2+68.5s+2445}$

其伯德图对比为：

可见，PI控制器在截止频率附件分量幅值有所放大，而P控制器没有，这会导致PI控制器有一些超调；在截止频率之前，PI控制器的相位延迟更小，因此，PI控制器的稳态误差更小一些，具体见其阶跃响应对比：

有的童鞋说，我不想用MATLAB，能不能手算PI参数呢？也是可以的，但是需要那么一丁点的经验。前面说了，系统的闭环传递函数是：

$\frac{C(s)}{R(s)}=\frac{k_ps/\tau+k_i/\tau}{s^2+(k_p+1)s/\tau +k_i/\tau}=\frac{k_ps/\tau+k_i/\tau}{s^2+2\xi\omega_ns +\omega_n^2}$

第三部分我们说到，传递函数特征方程的根的分布决定了系统的瞬态响应，这里我们先考察特征方程 $s^2+2\xi\omega_ns+\omega_n^2$ ，令 $\omega_n=50rad/s$ ， $\xi=0.7$ ，则 $k_p=2\xi\omega_n\tau-1$ ， $k_i=\tau\omega_n^2$ ，代入数据很容易得到： $k_p=69$ , $k_i=2500$ 。计算结果和matlab给出结果是比较接近的，这里的小经验就是截止频率 $\omega_n$ 的选择，因为有零点的存在，闭环系统的截止频率会更大一些（此处选择50rad/s计算，实际闭环截止频率是100rad/s左右），不过截止频率的选择本来范围就很大，在工程上这点误差都是可以接受的。