输入整形(Input Shaping)——一种振动抑制方法

最近刚开始了一个机器人的建模与控制方面的课题。其中要求机器人末端到达目标点后对其残留振动做一个抑制。按师兄的建议，去小学了一手工业中常用的输入整形(Input Shaping)的控制方法。

我们在抑制振动之前，首先要明确振动到底是什么。如果对控制原理有一定了解的话，我们知道，对于一个需要精确定位的系统，比如机械臂、塔吊、码头集装箱等，一般由于模型中存在着振荡模态，所以当其到达期望的位置时，并不能立刻停下来，而是会进行一段时间逐渐衰减的摆动，也就是我们所说的残留振动。
抑制振动的方法有很多，有通过增大零部件刚度、阻尼等通过改变系统模型而达到减小振动的被动控制的方法；有通过引入传感器形成闭环的反馈控制的方法；当然也有开环控制方法，我们这里重点探究的输入整形的方法便是开环控制的一种。

而提高零部件的刚度、阻尼势必会增大机构的重量，也会提高成本，而引入传感器不仅增加了系统的复杂性，还使控制难度也增大了。

什么是输入整形？

输入整形从字面意思就能够大致的理解了，就是把原本你期望的输入信号进行一步整形变换后再丢给后面的控制器。一个简单的控制系统是这样的：

而加入输入整形器后系统变成了这样：

显然，相当于只是在控制系统中加入了一个前馈环节。

输入整形器里面其实就是一系列的脉冲，它的作用就是把输入的信号与这一系列脉冲信号进行卷积后输出出来，成为整形信号再发给后面的系统。

输入整形的基本思想

试想有一辆小车放在桌上，车上用绳子挂了一个小球。我现在想让小车往前走1m，假如我直接给小车发的命令是走1m，那么根据常识我们很容易知道小车到1m停下来之后，小球并不会停下来，而是会来回摆动。

那么如果现在我换一种做法，我分两段走，先让小车走0.7m，看到小球摆到车前面时，再发命令让它走剩下的0.3m，这样一来，当小车停在1m的位置时，小球很有可能也就定住了。

其实从这个小例子可以看出输入整形的基本思想，就是一个输入信号不管怎么样都会在系统中产生一个振动，那么我就把原本的输入信号分几次给（最少两次），只要时间间隔合适，他们产生的振动通过线性叠加后，便可以实现相互抵消。比如下图，在0时刻处给一个A1脉冲信号，它产生的振动是红色线所示的衰减振荡，再在0.5s的时刻给一个A2的信号，它也会产生一个蓝色线所示的衰减震荡。然而当把它们俩产生的振动相加后，发现完美的抵消了。

输入整形器的实现

根据前面的介绍，输入整形器的关键就是计算出每一个脉冲信号的幅值以及它们各自的时滞。输入整形器的数学表示为
C(s)=∑Aie−tisC(s)=\sum A_ie^{-t_is} C(s)=∑Aie−tis
输入整形器的计算方法有好几种，这里只介绍一种可以计算任意个脉冲的整形器的计算方法——脉冲响应法。

以二阶系统为例，下面是一个典型的二阶系统的传递函数
G(s)=ωn2s2+2ξωns+ωn2(1)G(s)=\frac{\omega_n^2}{s^2+2\xi \omega_ns+\omega_n^2} \tag 1 G(s)=s2+2ξωns+ωn2ωn2(1)
它的单位脉冲响应为
g(t)=ωn1−ξ2e−ξωntsin⁡ωdt(2)g(t)=\frac{\omega_n}{\sqrt{1-\xi^2}}e^{-\xi\omega_nt}\sin\omega_dt \tag 2 g(t)=1−ξ2ωne−ξωntsinωdt(2)
其中ωn\omega_nωn为系统的自然频率，ωd\omega_dωd为系统的阻尼振荡频率，二者关系为ωd=ωn1−ξ2\omega_d=\omega_n\sqrt{1-\xi^2}ωd=ωn1−ξ2。

上文中提到的，整形器里面含有一系列的不同幅值和时滞的脉冲。假如一个脉冲的幅值和时滞分别为AiA_iAi和tit_iti，那么它引起的响应为
yi(t)=Aiωn1−ξ2e−ξωn(t−ti)sin⁡ωd(t−ti)(3)y_i(t)=\frac{A_i\omega_n}{\sqrt{1-\xi^2}}e^{-\xi\omega_n(t-t_i)}\sin\omega_d(t-t_i) \tag 3 yi(t)=1−ξ2Aiωne−ξωn(t−ti)sinωd(t−ti)(3)
系统总的响应为所有脉冲引起的响应之和：
y(t)=∑yi(t)=∑Aiωn1−ξ2e−ξωn(t−ti)sin⁡ωd(t−ti)(4)y(t)=\sum y_i(t)=\sum \frac{A_i\omega_n}{\sqrt{1-\xi^2}}e^{-\xi\omega_n(t-t_i)}\sin\omega_d(t-t_i) \tag 4 y(t)=∑yi(t)=∑1−ξ2Aiωne−ξωn(t−ti)sinωd(t−ti)(4)
用三角函数差角公式展开
y(t)=∑Aiωn1−ξ2e−ξωn(t−ti)(sin⁡ωdtcos⁡ωdti−cos⁡ωdtsin⁡ωdti)(5)y(t)=\sum \frac{A_i\omega_n}{\sqrt{1-\xi^2}}e^{-\xi\omega_n(t-t_i)}(\sin\omega_dt\cos \omega_dt_i-\cos \omega_dt\sin \omega_dt_i) \tag 5 y(t)=∑1−ξ2Aiωne−ξωn(t−ti)(sinωdtcosωdti−cosωdtsinωdti)(5)
稍作变形
y(t)=sin⁡ωdt⋅∑Aiωn1−ξ2e−ξωn(t−ti)cos⁡ωdti−cos⁡ωdt⋅∑Aiωn1−ξ2sin⁡ωdti)=ωn1−ξ2e−ξωnt⋅(sin⁡ωdt⋅∑Aieξωnticos⁡ωdti−cos⁡ωdt⋅∑Aieξωntisin⁡ωdti)(6)y(t)=\sin\omega_dt \cdot \sum \frac{A_i\omega_n}{\sqrt{1-\xi^2}}e^{-\xi\omega_n(t-t_i)}\cos \omega_dt_i-\cos\omega_dt\cdot\sum\frac{A_i\omega_n}{\sqrt{1-\xi^2}}\sin \omega_dt_i)\\ =\frac{\omega_n}{\sqrt{1-\xi^2}}e^{-\xi\omega_nt}\cdot (\sin\omega_dt \cdot \sum A_ie^{\xi\omega_nt_i}\cos \omega_dt_i-\cos\omega_dt \cdot \sum A_ie^{\xi\omega_nt_i}\sin \omega_dt_i) \tag 6 y(t)=sinωdt⋅∑1−ξ2Aiωne−ξωn(t−ti)cosωdti−cosωdt⋅∑1−ξ2Aiωnsinωdti)=1−ξ2ωne−ξωnt⋅(sinωdt⋅∑Aieξωnticosωdti−cosωdt⋅∑Aieξωntisinωdti)(6)
使用三角函数辅助角公式
y(t)=ωn1−ξ2e−ξωnt(∑Aieξωnticos⁡ωdti)2+(∑Aieξωntisin⁡ωdti)2⋅sin⁡(ωdt−ϕ)(7)y(t)=\frac{\omega_n}{\sqrt{1-\xi^2}}e^{-\xi\omega_nt}\sqrt{(\sum A_ie^{\xi\omega_nt_i}\cos \omega_dt_i)^2+(\sum A_ie^{\xi\omega_nt_i}\sin \omega_dt_i)^2}\cdot \sin(\omega_dt-\phi) \tag 7 y(t)=1−ξ2ωne−ξωnt(∑Aieξωnticosωdti)2+(∑Aieξωntisinωdti)2⋅sin(ωdt−ϕ)(7)
其中ϕ=arctan⁡∑Aieξωntisin⁡ωdti∑Aieξωnticos⁡ωdti\phi=\arctan \frac{\sum A_ie^{\xi\omega_nt_i}\sin \omega_dt_i}{\sum A_ie^{\xi\omega_nt_i}\cos \omega_dt_i}ϕ=arctan∑Aieξωnticosωdti∑Aieξωntisinωdti

很多文献里面都定义了一个叫做残留振荡百分比的东西：
V(ξ,ωn)=e−ξωntn(∑Aieξωnticos⁡ωdti)2+(∑Aieξωntisin⁡ωdti)2(8)V(\xi,\omega_n)=e^{-\xi\omega_nt_n}\sqrt{(\sum A_ie^{\xi\omega_nt_i}\cos \omega_dt_i)^2+(\sum A_ie^{\xi\omega_nt_i}\sin \omega_dt_i)^2} \tag 8 V(ξ,ωn)=e−ξωntn(∑Aieξωnticosωdti)2+(∑Aieξωntisinωdti)2(8)
（其实这里我并没有完全理解这个VVV的定义方式，我大概能猜到应该是把(7)式与(3)式作比，从而体现总响应与最后一次脉冲产生的振动的比较，然而里面的AnA_nAn去哪了就是我没相通的地方。如果有懂哥能够帮忙解释一下就再好不过了。不过这个倒不影响大局，人家这么定义的就这么拿来用就行了。）

我们的目标当然是消除残留振动，所以就要让这个VVV等于0。显然要达到这样一个目的，只能是根号等于0。也就是
{∑Aieξωnticos⁡ωdti=0∑Aieξωntisin⁡ωdti=0(9)\begin{cases} \sum A_ie^{\xi\omega_nt_i}\cos \omega_dt_i=0 \tag 9\\ \sum A_ie^{\xi\omega_nt_i}\sin \omega_dt_i=0 \end{cases} {∑Aieξωnticosωdti=0∑Aieξωntisinωdti=0(9)

ZV整形器

ZV整形器是输入整形最简单的实现方法了。它只有两个脉冲，且为了计算方便，让第一个脉冲从0时刻就开始作用。且为了使系统达到目标点，两次脉冲的幅值之和应该为1。

即约束条件如下：
t1=0A1+A2=1t_1=0\\ A_1+A_2=1 t1=0A1+A2=1
将条件带入(9)式
{A1+A2eξωnt2cos⁡ωdt2=0A2eξωnt2sin⁡ωdt2=0\begin{cases} A_1+A_2e^{\xi\omega_nt_2}\cos \omega_dt_2=0\\ A_2e^{\xi\omega_nt_2}\sin \omega_dt_2=0 \end{cases} {A1+A2eξωnt2cosωdt2=0A2eξωnt2sinωdt2=0
由第二行可知ωdt2=π\omega_dt_2=\piωdt2=π，即t2=πωdt_2=\frac{\pi}{\omega_d}t2=ωdπ

再将ωdt2=π\omega_dt_2=\piωdt2=π带入第一行得A1−A2eξωnt2=0A_1-A_2e^{\xi\omega_nt_2}=0A1−A2eξωnt2=0，再结合A1+A2=1A_1+A_2=1A1+A2=1可解出两个脉冲的幅值
A1=eξωnt21+eξωnt2,A2=11+eξωnt2A_1=\frac{e^{\xi\omega_nt_2}}{1+e^{\xi\omega_nt_2}},A_2=\frac{1}{1+e^{\xi\omega_nt_2}} A1=1+eξωnt2eξωnt2,A2=1+eξωnt21
这样一来ZV整形器也就设计完成了。

仿真的效果如下：

未整形：

模型准确时振动抑制效果：

模型有10%的误差时振动抑制效果：

可以看出，ZV整形器在模型准确时对振动的抑制效果非常好，且收敛的也非常快；然而一旦模型不准确，则还是会出现很多个周期的残留振动。

ZVD整形器

当模型参数不准确时，ZV整形器的零点不能有效的对消系统的零点，从而不能有效的消除残留振动。为了提高整形器的鲁棒性，对式V(ξ,ωn)V(\xi,\omega_n)V(ξ,ωn)的微分方程进行约束，为系统配置多重零点，这种整形器称为ZVD整形器，也就新引入了一个约束：
∂∂ωV(ξ,ω)=0\frac{\partial}{\partial \omega}V(\xi,\omega)=0 ∂ω∂V(ξ,ω)=0
那么要实现这个约束条件，就要增加一个脉冲，即3个脉冲实现输入整形。计算过程不在赘述了，直接给出计算结果：
A1=11+2e−ξπ1−ξ2+e−2ξπ1−ξ2A2=2e−ξπ1−ξ21+2e−ξπ1−ξ2+e−2ξπ1−ξ2A3=e−2ξπ1−ξ21+2e−ξπ1−ξ2+e−2ξπ1−ξ2t1=0t2=πω1−ξ2t3=2πω1−ξ2A_1=\frac{1}{1+2e^{\frac{-\xi\pi}{\sqrt{1-\xi^2}}}+e^{\frac{-2\xi\pi}{\sqrt{1-\xi^2}}}}\\ A_2=\frac{2e^{\frac{-\xi\pi}{\sqrt{1-\xi^2}}}}{1+2e^{\frac{-\xi\pi}{\sqrt{1-\xi^2}}}+e^{\frac{-2\xi\pi}{\sqrt{1-\xi^2}}}}\\ A_3=\frac{e^{\frac{-2\xi\pi}{\sqrt{1-\xi^2}}}}{1+2e^{\frac{-\xi\pi}{\sqrt{1-\xi^2}}}+e^{\frac{-2\xi\pi}{\sqrt{1-\xi^2}}}}\\ t_1=0\\ t_2=\frac{\pi}{\omega\sqrt{1-\xi^2}}\\ t_3=\frac{2\pi}{\omega\sqrt{1-\xi^2}} A1=1+2e1−ξ2−ξπ+e1−ξ2−2ξπ1A2=1+2e1−ξ2−ξπ+e1−ξ2−2ξπ2e1−ξ2−ξπA3=1+2e1−ξ2−ξπ+e1−ξ2−2ξπe1−ξ2−2ξπt1=0t2=ω1−ξ2πt3=ω1−ξ22π
仿真效果如下：

系统模型准确时振动抑制效果

模型有10%的误差时振动抑制效果：

和ZV整形器比较而言，当模型准确时，ZV整形器的快速性显然要优于ZVD整形器，二者的整形效果当然都是非常好的；而一旦当我们拿到的系统模型不准确时，ZVD的振动抑制效果就优于ZV了

其他整形器

在文献里面还看到了一些其他的整形器，如极不灵敏整形器、最优整形器等，就有待进一步的探究了

参考文献
[1]陈俊恒. 输入整形减振算法的研究与实现[D].哈尔滨工业大学,2010.
[2]秦望舒. 基于输入整形技术的运动控制算法的研究[D].哈尔滨工业大学,2010.