云台 PID 控制理论分析¶

Intro¶

本文以 Robomaster 中常见的英雄机器人云台 Pitch 轴控制为例，重点分析常规电机角度控制时的 PID 设计。

在分析时，本文忽略了执行器能力有限，PID输出应当限幅的问题；不考虑微分项、积分项在工程实践中的近似处理；不考虑阻力和外界扰动（假设前馈项能完美抵消）；不考虑工程实践时反馈数据的噪声部分。这些问题或许会留到之后探讨。

本文需要一点常微分方程、 Laplace 变换和复变函数的知识，如果读者没有接触过，可以先阅读一下附录，或许有所帮助。

1 建模¶

我们取云台水平位置为云台角度 \(\theta\) 的零位，向上转为正，向下转为负。假定云台重心位于 \(\theta_0\) 的位置，重力为 \(G\)，与转动轴心连线长度为 \(l\)。假定阻力矩恒定为 \(f\)，正负号与转动方向相反。用 \(J\) 表示云台绕 Pitch 轴转动的转动惯量，用 \(T\) 表示云台 Pitch 轴电机的输入力矩。从而得到

我们采用 PID 控制器 + 前馈的策略计算输入力矩 \(T\)，从而控制云台角度 \(\theta\)。【假设我们的前馈是完美前馈，能完全抵消掉重力矩和阻力矩】（注：这是本文的一个重要假设，所有分析都建立在该假设成立的情况下）。那么，我们有：

分析我们常见的云台 Pitch 控制方案，控制板读取遥控器/键鼠指令后，在原有 Pitch 期望角度信号上加上一个值，因此 Pitch 期望角度信号可以看成是一个阶跃变化的信号。我们只分析一次阶跃变化时的请况，设这个阶跃变化量 \(\Delta\theta_{ref} = r\).

施加新的期望信号的同时，云台的 Pitch 往往同时具有一个不为 0 的期望角度/实际角度/实际角速度/实际角加速度初值。我们记这些初值为 \(\theta_{ref}(0) = A, \theta(0) = a, \dot\theta(0) = a_{d}, \ddot\theta(0) = a_{dd}\) .

2 PID 控制器¶

PID 控制器在复频域上可以表示为

其中，若积分时间 \(T_i \to \infty\)，则相当于没有积分项；若微分时间 \(T_d = 0\)，则相当于没有微分项。我们默认 PID 里的参数都大于 0.

3 单环 PID 稳定性分析¶

如果我们用一个单环的 PID 控制器来控制 Pitch，控制策略可以表示为

对等式两边做拉普拉斯变换，代入相关符号和 PID 控制器表达式，有

把相关符号代入，得到

可以解出 \(\Theta\) 的表达式：

代入 PID 控制器的表达式，得到

这里的 \(\Theta_{ref} = \mathcal{L}(\theta_{ref}) = \mathcal{L}(\theta_{ref}(0) + \Delta\theta_{ref}) = \mathcal{L}(A+r) = (A+r)/s\)，\(\mathcal{L}(\cdot)\) 表示 Laplace 变换 . 因此有

3.1 暂态过程的稳定性¶

(3.2) 式对角相乘后，Laplace 反变换可以得到

其中，\(D^n\theta = \frac{\mathrm{d}^2\theta(t)}{\mathrm{d}t^2}\)，\(C\) 是一个常数项。

根据常微分方程的知识，\(\theta\) 的解可以由两部分相加得到：\(\theta = \theta_{ss} + \theta_{t}\) . 其中，\(\theta_{ss}\) 是 (3.3) 微分方程的一个特解，又叫稳定解。\(\theta_t\) 是 (3.3) 微分方程的一个通解，也叫暂态解，也就是 (3.4) 方程的解：

通常情况下，稳态解描述了 \(t \to \infty\) 时 \(\theta\) 的样子，而暂态解主要描述了从开始变化到最终稳定之间的过渡过程中 \(\theta\) 的样子。我们先来看看暂态解。

根据常微分方程的知识，这个暂态解的形式是 \(\theta(t) = \sum C_ie^{\lambda_it}\) ，\(\lambda_i\) 是复数，如果这个和式解中存在一个 \(\lambda_i = a+bj\) （虚部不为零）的部分，那么必然同时存在一个 \(\lambda_j = a-bj\) 的共轭部分。

这里所有的 \(\lambda\) ，如果【实部全都小于 0】 ，说明 \(t\to \infty\) 的时候暂态解会趋近于 0，那么这个 \(\theta(t)\) 就是【收敛的、稳定的】；如果【存在一个实部大于 0 的】\(\lambda\)，说明 \(t \to \infty\) 的时候暂态解会趋近于 \(\infty\)，那么这个 \(\theta(t)\) 就是【发散的、不稳定的】；如果【存在两个实部为 0 ，虚部不为 0 的】 \(\lambda_i\) 和 \(\lambda_j\) ，那么暂态解里面就可以凑个三角函数出来，\(\theta(t)\) 就会有一个【等幅振荡】的成分.

取 \(\sum C_ie^{\lambda_it}\) 的其中一项代入 (3.4) 方程，可以得到：

显然，所有的 \(\lambda\) 都是 \(J\lambda^3+K_cT_d\lambda^2+K_c\lambda+K_c/T_i = 0\) 的根。这里需要注意一下，如果 \(\lambda\) 存在一个 p 重根的话，暂态解的形式略有区别，这个 p 重根 \(\lambda\) 会对应于 \(C_1e^{\lambda t} + C_2te^{\lambda t} + ... + C_{p}t^{p-1}e^{\lambda t}\).

不难发现，\(J\lambda^3+K_cT_d\lambda^2+K_c\lambda+K_c/T_i = 0\) 这个关于 \(\lambda\) 的方程，和 【\(\Theta / \Theta_{ref}\) 的分母等于 0 这个关于 \(s\) 的方程】 ，本质上是同一个方程，根都是一样的，我们管这个方程叫作【特征方程】，把特征方程的解叫作【特征根】。

3.2 稳定时的值是多少¶

我们一般利用 Laplace 变换的终值定理（参见附录）分析 \(\theta(t)\) 稳定时的值是多少。

注意，要使用 Laplace 变换的终值定理，我们需要保证 \(s\Theta(s)\) 在包含虚轴的右半平面上是解析的，也就是说，\(s\Theta(s)\) 的分母在右半平面上没有根，其实就是特征方程的解的实部全都小于 0.

这正好和 3.1.1 节的内容对上了。其实想想也知道，\(\theta(t)\) 必须得是稳定的，才能算它稳定时的值是多少。

3.3 单环纯比例控制¶

此时 \(T_d = 0, T_i \to \infty\) .

特征方程为

特征根为 \(s = \pm i\sqrt{K_c/J}\) ，所以此时【 \(\theta(t)\) 是等幅振荡的】。

3.4 单环 PD 控制¶

此时 \(T_d \neq 0, T_i \to \infty\) .

特征方程为

由于特征根的实部都小于 0，所以此时【 \(\theta(t)\) 是稳定的】。接下来利用 Laplace 变换的终值定理算一下稳定时候的值：

\(A+r\) 是期望角度的初始值加上期望角度的变化量，也就是最终的期望角度。也就是说，在前馈足够好的情况下，【只用单环 PD 就能使 Pitch 角度无稳态误差地跟随阶跃信号】。

3.5 单环 PI 控制¶

此时 \(T_d = 0, T_i\)有界 .

特征方程为

解这样一个三次方程有点太复杂了。我们分析一下吧。

首先，根肯定不会在虚轴上。因为倘若 \(s = j\omega\)，那么 \(K_c/(JT_i)\) 这一项肯定无法抵消为 0。

其次，我们考虑这个三次方程的实数解。令 \(h(s) = s^3 + K_c/J \cdot s + K_c/(JT_i)\)，求导得到 \(\dot{h}(s) = 3s^2 + K_c/J > 0\)，所以 \(h(s)\) 单调递增。又因为 \(h(0) = K_c/(JT_i) >0, h(-\infty) = -\infty\)，所以这个三次方程只有一个小于 0 的实数解。

接下来，我们分析剩下的两个复数解。不妨把刚才发现的实数解记为 \(s_1 = \zeta < 0\) . 改写方程为

因此，两个复根可以表示为

可以发现，这两个复根的实部都是正数。也就是说，此时【\(\theta(t)\) 是发散的】。

并且，【\(h(0) = K_c/(JT_i)\) 越接近 0】, \(\zeta\) 就越接近 0，发散部分的幅度就越小，【\(\theta(t)\) 就越接近于等幅振荡】。

3.6 单环 PID 控制¶

定性理解，相当于 PI + PD，参数选择恰当就能稳定，参数选择不当就会发散。

在稳定的情况下，计算其稳定时的值：

没有稳态误差。

想分析清楚 \(\theta(t)\) 在什么条件下稳定，什么条件下发散太过复杂，超出我能力范围了。

放个 MATLAB 求解程序在这里。

% 定义参数
J = 1;  % Pitch 轴转动惯量
A = 1;  % 初始期望角度
a = 1;  % 初始角度
r = 1;  % 期望角度变化量
ad = 1; % 初始角速度
Kc = 1; % PID 比例项
Td = 4; % PID 微分时间
Ti = 1; % PID 积分时间

% 拉普拉斯变换表达式
syms s t
F = ( J*a*s^3 + ((A+r)*Kc*Td+J*ad)*s^2 + (A+r)*Kc*s+ (A+r)*Kc/Ti ) / ...
  (s*(J*s^3 + Kc*Td*s^2 + Kc * s + Kc/Ti));

% 计算反变换
f_t = ilaplace(F, s, t);

% 绘制f(t)图像
fplot(f_t, [0, 40])  
xlabel('Time (s)')
ylabel('f(t)')
title('f(t) vs Time')
grid on

3.7 小结¶

综上所述，【单环 PID 控制的稳定性相当差，很容易发散】。这也是为什么我们要做双环 PID 控制。

4 双环 PID 分析¶

对于云台 Pitch 角度的双环 PID 控制（又叫串级控制）策略，一般是角度环+速度环，即一个 PID 根据角度误差计算出速度期望，另一个 PID 接收速度期望和实际速度，计算出输出力矩。（顺带一提，内环参数为什么选择“速度”？参见 附录 -2）

我们的双环 PID 控制策略可以表示为

对等式两边做拉氏变换，可以得到

把相关符号代入，得到：

可以解出 \(\Theta\) 的表达式：

4.1 速度环选什么？¶

速度环的控制策略可以表示为

其中带 \(d\) 下标的代表是角速度，\(V\) 代表角速度期望值阶跃变化量。

用第 3 节的方法很容易分析出，【不管怎么选择 PID 参数，速度环都不会发散，最终稳定时都没有稳态误差】。

在 MATLAB 的帮助下，我们可以直接解出 \(\dot\theta\) 的形式：

% 拉普拉斯变换表达式
syms s t J ad Kc Td Ti V
F = ( ((V+ad)*Kc*Td + J*ad)*s^2 + (V+ad)*Kc*s + (V + ad)*Kc/Ti) / ...
  (s*((J+Kc*Td)*s^2 + Kc*s + Kc / Ti));
% 计算反变换
f_t = ilaplace(F, s, t);

稍微整理一下：

先关注微分项参数 \(T_d\) ，它主要和 \(\frac{J V}{J + K_c T_d} \exp\left(-\frac{K_c}{2(J+K_c T_d)}t\right)\) 有关。 随着 \(T_d\) 增大，暂态解的系数 \(JV / (J+K_cT_d)\) 减小，指数中的系数 \(-K_c/\left(2(J+K_c T_d)\right)\) 绝对值减小，这说明【微分项的引入能够抑制过渡过程的振荡幅度，但是会延长过渡过程的持续时间】。 这不好，因为我们希望内环参数变化灵敏，调节时间短。所以我们速度环不用微分项。

再关注积分项参数 \(1/T_i\) ，它主要和 \(\left( \cosh\left( \omega t \right) - \frac{ \sinh\left( \omega t \right)}{\sqrt{1-4/T_i \cdot \left( T_d+J/K_c\right)}} \right)\) 有关。 随着 \(1/T_i\) 增大，\(\omega\) 从 1 减小到 0 再变为一个虚数。先考虑 \(\omega\) 是个正数的情况：

此时， \(\left( \cosh\left( \omega t \right) - \frac{ \sinh\left( \omega t \right)}{\sqrt{1-4/T_i \cdot \left( T_d+J/K_c\right)}} \right)\) 可以写成 \((1-k)(e^{-kx} - e^{kx})/(2k)\) 的形式，这里的 \(k = \sqrt{1-4/T_i \cdot \left( T_d+J/K_c\right)} < 1\)，\(x = t / [2(J/K_c + T_d)]\) ，是个随着 \(k\) 减小而递增的值。

也就是说，当 \(T_i > 4(T_d+J/K_c)\) 时，\(1/T_i\) 越大，过渡过程的振荡幅度越大。

再考虑 \(\omega\) 是个虚数，即 \(\omega = \pm j W\) 的情况，由于 \(\left( \cosh\left( \omega t \right) - \frac{ \sinh\left( \omega t \right)}{\sqrt{1-4/T_i \cdot \left( T_d+J/K_c\right)}} \right)\) 关于 \(W\) 是个偶函数，所以我们只需要考虑 \(W\) 为正数的情况，此时

显然，这是一个等幅振荡的形式。随着 \(1/T_i\) 增大，\(W\) 也增大，\(\sqrt{1+W^2}/W\) 趋近于 1，\(\cos(Wt+\phi)\) 的角频率增大，相位角减小。

也就是说，当 \(T_i\) 从 \(>4(T_d+J/K_c)\) 转变到 \(<4(T_d+J/K_c)\) 后，过渡过程中出现了等幅振荡的部分，其稳定性完全依赖于 \(\exp\left(-\frac{K_c}{2(J+K_c T_d)}t\right)\)，调节时间延长了。随着 \(1/T_i\) 继续增大，等幅振荡部分的幅度减小，振荡频率增大。

综上所述，【积分项的引入会导致过渡过程振荡加剧、延长】。这也不好。所以我们速度环也不应该使用积分项。

留个速度环的 MATLAB 计算和画图程序在这里。

% 定义参数
J = 1;  % Pitch 轴转动惯量
ad = 1; % 初始角速度
Kc = 1; % PID 比例项
Td = 0.1; % PID 微分时间
Ti = 0.01; % PID 积分时间
V = 5; % 期望角速度变化量

% 拉普拉斯变换表达式
syms s t
F = ( ((V+ad)*Kc*Td + J*ad)*s^2 + (V+ad)*Kc*s + (V + ad)*Kc/Ti) / ...
  (s*((J+Kc*Td)*s^2 + Kc*s + Kc / Ti));

% 计算反变换
f_t = ilaplace(F, s, t);

% 绘制f(t)图像
fplot(f_t, [0, 100])  
xlabel('Time (s)')
ylabel('f(t)')
title('f(t) vs Time')
grid on

4.2 双环 PID 性能分析¶

根据 4.1 节的结论，速度环的 PID，我们只应该使用比例环节，即 \(\text{PID}_{vel} = K_p\) . 从而，我们改写 (4.1) 式，得到：

代入 PID 控制器的表达式，得到：

这里的 \(\Theta_{ref} = \mathcal{L}(\theta_{ref}) = \mathcal{L}(A+r) = (A+r)/s\) . 因此有

2.1 角度环纯比例¶

此时，\(T_d = 0, T_i \to \infty\) .

特征方程为

显然，两个特征根的实部都小于 0，所以此时【 \(\theta(t)\) 是稳定的】。接下来利用 Laplace 变换的终值定理算一下稳定时候的值：

可以发现，Pitch 实际角度最终与期望角度保持一致。也就是说，【只要前馈足够理想，纯比例就能实现 Pitch 角度无稳态误差】。

用 MATLAB 对 \(\Theta(s)\) 做 Laplace 反变换，可以得到

首先，观察暂态解指数项的系数 \(K_p/(2J)\)，可以发现此时暂态过程的持续时间与外环比例项 \(K_c\) 无关，仅与内环比例项 \(K_p\) 有关。【内环比例项 \(K_p\) 越大，暂态过程持续时间越短，暂态过程振荡幅度越小。下同】。

考虑 \(K_c\) 从 \(<K_p/(4J)\) 的情况下逐渐增大，\(\omega\) 将从一个正数逐渐减小至 0 ，再转变为一个虚数，虚部逐渐增大。

这里我不知道怎么分析，但是结果是，【外环比例项 \(K_c\) 越大，暂态过程振荡越剧烈】。

2.2 角度环 PD¶

此时，\(T_d \neq 0, T_i \to \infty\)。

特征方程为

显然，两个特征根的实部都小于 0，\(\theta(t)\) 稳定，并且可以用 Laplace 变换的终值定理算出稳定值为 \(A+r\) .

可以发现，微分项的引入不改变控制器的稳态特性。我们分析一下暂态特性。

这个部分是最开始写的。刚开始不知道 MATLAB 能直接形式化做 Laplace 反变换，自己去分类讨论算了两天，结果还有错的哈哈

自己做 Laplace 反变换还是太吃操作了，还好我有简单有强势的 MATLAB ：

% 拉普拉斯变换表达式
syms s t J a Kp Kc Td Ti A r ad
F = ( J*a*s^2 + (Kp*Kc*Td*(A+r) + Kp*a + J*ad)*s + Kp*Kc*(A+r) ) / ...
  ( s * (J*s^2 + Kp*(Kc*Td+1)*s + Kp*Kc ) );

% 计算反变换
f_t = ilaplace(F, s, t);

随着 \(T_d\) 增大，\(\omega\) 增大，暂态解系数 \(1/(J^2\omega)\) 减小，暂态解指数项里的系数绝对值 \((K_p+K_cK_pT_d)/(2J)\) 增大，说明【微分项的引入能有效抑制暂态过程的振荡，但会延长暂态过程的持续时间】。

2.3 角度环 PI¶

此时，\(T_d = 0, T_i\)有界。

特征方程

令 \(m(s) = Js^3+K_ps^2+K_pK_cs+K_pK_c/T_i\)，求个导得到 \(\dot{m}(s) = 3Js^2+2K_ps+K_pK_c\)，分别考虑有两个实根的情况和没有两个实根的情况。

有两个实根时，这两个根分别为 \(s = \frac{-K_p \pm K_p\sqrt{1-3K_c/K_p}}{3}\)，都为负数，所以 \(m(s)\) 在先增再减再增，又因为 \(m(-\infty) = -\infty, m(0) > 0\)，所以 \(m(s)\) 的三个根都是负实数，\(\Theta\) 是稳定的。

\(\dot{m}(s)\) 没有两个实根时，\(m(s)\) 单调递增，因为 \(m(0) > 0\)，所以 \(m(s)\) 只有一个负实根 \(\delta\)，从而可以将 \(m(s)\) 写成 \(Js^2 + (K_p+J\delta)s+K_p(K_c+\delta)+J\delta^2\) ，显然另外两个复根的实部也是负数，因此 \(\Theta\) 是稳定的。

稳态值是 \((A+r)\)，无误差。

暂态过程分析不出来了。大概测试结果是【积分项的引入会导致暂态过程振动幅度增大，暂态过程时间延长】。

在不考虑扰动的情况下，积分的引入带来的似乎只有负面作用。

附录¶

0 前置知识¶

0.1 描述线性系统的基本方式（微分方程）¶

对于一个线性系统 \(G\)，输入 \(r(t)\)，输出 \(y(t)\)，我们通常用微分方程的形式来描述这一系统：

其中，\(a_{i}, b_{j}\) 是常系数，\(D^{k}y\) 代表 \(\frac{\mathrm{d}^ky(t)}{\mathrm{d}t^k}\)。

如果是现实世界中存在的系统，考虑到因果性，必然有 \(n > m\)。

这样一个系统，输入信号后，输出信号可以分成两个部分，一个是暂态输出，指的是 \(t\to \infty\) 后会变成 0 的部分；另一个是稳态输出，指的是 \(t\to \infty\) 后仍然存在的部分。

0.2 拉普拉斯变换¶

拉氏变换可以将连续的时域信号/函数转化为复频域的信号/函数。

对于一个时域上的函数 \(f(t)\)，定义 拉普拉斯变换（以下简称拉氏变换）：

其中，\(\mathcal{L}(f(t))\) 表示对 \(f(t)\) 做拉氏变换。（注：拉氏变换要求 \(f(t)\) 在 \((0,+\infty)\) 上有定义）

拉氏变换实际上是将时域上的函数转化到复频域，有一些特殊的性质，例如：

时域位移相当于复频域乘以 \(e^{-sT}\)。

证明：

只要 \(f(t)\) 在 \(t < 0\) 的部分均等于 0，那么上式就可以简化为

复频域上乘以 s 相当于时域求微分。

证明：

其中，\(F(s) = \mathcal{L}(f(t))\)。在以上推导过程中，我们假定了 \(f(0_+) = f(0_-)\)。如果 \(f(t)\) 的性质特别好，满足 \(t = 0\) 时，\(f(t)\) 的任意阶导数值均为 0，那么上式可以进一步简化为

拉氏变换可以用来分析微分方程，大大简化我们的分析难度。

比如，对于 1.1 节中的系统，我们输入 \(r(t) = 1\)，求输出 \(y(t)\)？直接解微分方程不太容易，但是用拉氏变换分析就会简单不少：

首先，我们在零初始条件下，将系统的微分方程等式两边做拉氏变换，得到

上面这个比值形式就是系统的传递函数。我们考虑一个简单一点的情况：

因此，我们有

对上式做拉氏反变换，即可得到 \(y(t) = 1 - e^{-t}\)。（注：\(\mathcal{L}(1) = \frac{1}{s}\) 很容易证明；拉氏反变换一般是直接查表）

0.3 Laplace 变换初值定理和终值定理¶

若 \(\mathcal{L}[f(t)] = F(s)\)，\(\mathrm{d}f/\mathrm{d}t\) 是可以 Laplace 变换的，且 \(f(+\infty)\) 存在，\(sF(s)\) 在包含虚轴的右半平面上解析（原点除外），则有终值定理：

若 \(\mathcal{L}[f(t)] = F(s)\)，\(\mathrm{d}f/\mathrm{d}t\) 是可以 Laplace 变换的，且 \(\lim_{s\to \infty}sF(s)\) 存在，则有初值定理

-1 复变函数相关知识¶

复变函数中最重要的一个式子是欧拉公式 \(e^{j\theta}=\cos\theta+j\sin\theta\)，其中 \(\cos\theta\) 是复数的实部，\(j\sin\theta\) 是复数的虚部。根据欧拉公式，我们可以很简单地推导出，\(j = e^{j\frac{\pi}{2}}\)。

事实上，在复变函数的理论中，所有的复数 \(\mathbf{A}\) 都可以由 \(Ae^{j\theta}\) 的形式表示。\(\mathbf{A}\) 可以理解为一个复平面上的矢量，从原点指向 \((A\cos\phi,A\sin\phi)\) 的位置。通常，我们用幅度 \(A\) 和相位角 \(\phi\) 来描述这个矢量。定义一个求幅度的运算：\(|\mathbf{A}| = |Ae^{j\theta}| = A\)。

任意一个复数乘以 \(Ae^{j\theta}\)，表示将该复数对应的复平面矢量的幅度乘以 \(A\)，角度转 \(\theta\)。相反地，除以 \(Ae^{j\theta}\)，就是幅度除以 \(A\)，角度反转 \(\theta\)。

根据上述知识，我们可以得到：

-2 串级控制选取内环参数的原则¶

双环 PID 中，外环参数是固定的，我们想控制 Pitch 的角度，那么外环就必须接收角度的期望值和反馈量；但内环参数（比如”速度环“的”速度“）不是固定的，可以由设计者选择。一般来说，选取内环参数需要满足以下几点要求：

1. 可以被测量；
2. 比外环参数能更快地感知干扰；
3. 直接受执行器输出的影响；
4. 对外环参数有明确且较大的影响；
5. 比外环参数的反应更加灵敏，调节时间更短；
6. 最好包含系统的非线性或者时变部分。