机械振动

1. 振动的阻抗表示

机械输入阻抗为输入力复幅值与驱动点速度之比。

1.1 无阻尼系统

对于单惯量无阻尼系统：

$$M\ddot{x}(t)+Kx(t)=f(t)$$

令：$x(t)=Ae^{j\omega t}$，$f(t) = F e^{j\omega t}$

那么：$\dot{x}(t)=j\omega Ae^{j\omega t}$，$\ddot{x}(t)=-\omega^2 Ae^{j\omega t}$

则输入阻抗为：

$$Z_i=\frac{f(t)}{\dot{x}(t)}=\frac{F}{j\omega A}$$

令：$\omega_n=\sqrt{\frac{K}{M}}$

那么：

$$Z_i=\frac{f(t)}{\dot{x}(t)}=\frac{-jK(1-{\frac{\omega}{\omega_n}}^2 )}{\omega}$$

联立得到：

$$\frac{AK}{F}=\frac{1}{1-{\frac{\omega}{\omega_n}}^2}$$

其中$\frac{AK}{F}$为无量纲（归一化）位移，$\frac{\omega}{\omega_n}$为无量纲（归一化）输入频率。

1.2 有阻尼系统

对于单惯量有阻尼系统：

$$M\ddot{x}(t)+C\dot{x}(t)+Kx(t)=f(t)$$

令：$x(t)=Ae^{j\omega t}$，$f(t) = F e^{j\omega t}$

那么：$\dot{x}(t)=j\omega Ae^{j\omega t}$，$\ddot{x}(t)=-\omega^2 Ae^{j\omega t}$l

则输入阻抗为：

$$Z_i=\frac{f(t)}{\dot{x}(t)}=\frac{F}{j\omega A}$$

令$\zeta=\frac{C}{2M\omega_n}$

$$Z_i=\frac{f(t)}{\dot{x}(t)}=\frac{-jK(1+j2\zeta\frac{\omega}{\omega_n}-{\frac{\omega}{\omega_n}}^2 )}{\omega}$$

联立得到：

$$\frac{AK}{F}=\frac{1}{1-{\frac{\omega}{\omega_n}}^2+j2\zeta\frac{\omega}{\omega_n}}$$

其中$\frac{AK}{F}$为无量纲（归一化）位移，$\frac{\omega}{\omega_n}$为无量纲（归一化）输入频率。

有阻尼系统位移是复数，幅值表示最大振幅，相角表示与输入力的相位差。

2. 闭环控制

$$M\ddot{x}(t)+C\dot{x}(t)+Kx(t)=Du(t)+Ef(t)$$

其中D和E分别是控制力输入矩阵和外部扰动力矩阵。

M，C，K是nXn维矩阵，D是nXm维矩阵，E是nXr维矩阵。

认为控制输入$u(t)$是位移$x(t)$，速度$\dot{x}(t)$，外力$f(t)$的线性组合：

$$u(t)=K_1x(t)+C_1\dot{x}(t)+E_1f(t)$$

其中$K_1,C_1,E_1$可以是时变的。

$$M\ddot{x}(t)+(C-DC_1)\dot{x}(t)+(K-DK_1)x(t)=(E+DE_1)f(t)$$

系统闭环控制就是通过输入改变系统的阻尼和刚度。

3. 状态空间方程

令$z(t)=\begin{bmatrix} x(t)\\ \dot{x}(t) \end{bmatrix}$：

$$\dot{z}(t)=A{z}(t)+Bu(t)+Hf(t)$$

其中：$A=\begin{bmatrix} 0 & I\\ -M^{-1}K & -M^{-1}C \end{bmatrix}$，$B=\begin{bmatrix} 0 \\ M^{-1}D \end{bmatrix}$，$H=\begin{bmatrix} 0 \\ M^{-1}E \end{bmatrix}$

这就是系统状态空间方程的形式。

4. 最优控制

取惩罚函数为：

$$J = \int_{0}^{t_f} [z^T(t)Qz(t)+u^T(t)Ru(t)]dt$$

其中Q为半正定矩阵，R为正定矩阵。

拉格朗日方程：

$$L=\int_{0}^{t_f}[z^T(t)Qz(t)+u^T(t)Ru(t)+λ^T(t)(Az(t)+Bu(t)+Hf(t)−\dot{z}(t))]dt$$

令$\mathcal{H}=z^TQz+u^TRu+λ^T(Az+Bu+Hf)$

$$L=\int_{0}^{t_f}(\mathcal{H}−λ^T\dot{z})dt$$

对J取一阶变分：$z(t)→z(t)+δz(t)，u(t)→u(t)+δu(t)$，$λ(t)，f(t)$当做常值，$δz(0)=0$

$$δL=\int_{0}^{t_f}((\frac{∂\mathcal{H}}{∂z})^Tδz+(\frac{∂\mathcal{H}}{∂u})^Tδu−λ^Tδ\dot{z})dt$$

将 $λ^Tδ\dot{z}$ 项进行分部积分：

$$\int_{0}^{t_f}λ^Tδ\dot{z}dt={λ(t)^Tδz(t)}\mid^{t_f}_{0}−\int_{0}^{t_f}\dot{λ}^Tδz dt$$ $$δL=-{λ(t)^Tδz(t)}\mid^{t_f}_{0}+\int_{0}^{t_f}[(\frac{∂\mathcal{H}}{∂z}+\dot{λ}^T)^Tδz+(\frac{∂\mathcal{H}}{∂u}) ^Tδu]dt$$

设$δz(0)=0$，要使得$δL=0$:

$$\dot{λ}=-\frac{∂\mathcal{H}}{∂z}，λ(t_f)=0$$ $$\frac{∂\mathcal{H}}{∂u}=0$$

得到最终结果：

$$\dot{λ}=-A^Tλ-2Qz，λ(t_f)=0$$ $$u=-\frac{1}{2}R^{-1}B^Tλ$$

4.1 闭环控制（LQR）

令$λ(t)=P(t)z(t)$：

$$[\dot{P}(t) + P(t )A -\frac{1}{2} P(t)BR^{-1}B^TP(t) +A^TP(t) + 2Q] z(t) + P(t )Hf(t) = 0，P(t_f)=0$$

当$f(t) = 0$时：

$$\dot{P}(t) + P(t )A -\frac{1}{2} P(t)BR^{-1}B^TP(t) +A^TP(t) + 2Q，P(t_f)=0$$ $$u(t)=G(t)z(t)=-\frac{1}{2}R^{-1}B^TP(t)z(t)$$ $$G(t)=-\frac{1}{2}R^{-1}B^TP(t)$$

即为控制增益。

若$f(t) \ne 0$，此时并不是最优解，需要结合前馈。在结构控制中$P(t)$一般在过程中不变，只在快结束时快速下降。因此可以简化为常数矩阵：

$$PA+A^TP−PBR^{−1}B^TP+2Q=0$$ $$G=-\frac{1}{2}R^{-1}B^TP$$

原状态空间方程变为：

$$\dot{z}(t)=(A+BG){z}(t)+Hf(t)$$

由开环的矩阵$A$变为闭环矩阵$A+BG$。

4.1.1 状态观测

这种控制方式需要了解全部的状态$z(t)$，实际很少能测得全部的状态:

$$y(t)=Cz(t)+\eta$$

其中$\eta$是测量噪声。需要进行状态估测得到$\hat{z}(t)$。当信噪比足够高时使用Luenberger观测器，否则使用Kalman滤波器：

$$\hat{\dot{z}}(t)=A\hat{z}(t)+Bu(t)+L(C\hat{z}(t)-y(t))$$

其中$L$如果使用Luenberger观测器设计就是一个常数矩阵，如果使用Kalman滤波器是一个时变矩阵，需要通过解另一个Riccati方程实时计算，计算量比较大。

4.2 闭开环控制和开环控制

如果能够测得外部力，那么可以更好地控制。

令$λ(t)=P(t)z(t)+S(t)f(t)$：$[\dot{P}(t) + P(t )A -\frac{1}{2} P(t)BR^{-1}B^TP(t) +A^TP(t) + 2Q] z(t) + \dot{S}(t)f(t)+S(t)\dot{f}(t)-(\frac{1}{2} P(t)BR^{-1}B^T-A^T)S(t)f(t)+P(t )Hf(t) = 0，P(t_f)=0， S(t_f)=0$
$P(t)$部分依然按照前述处理，剩余部分：

$$[\dot{S}(t)-(\frac{1}{2} P(t)BR^{-1}B^T-A^T)S(t)+P(t )H]f(t) +S(t)\dot{f}(t)= 0，S(t_f)=0$$

由于$S(t)$依赖于$f(t)，\dot{f}(t)$无法预先得知，因此无法计算最优前馈。

5. 极点配置

$$\dot{z}(t)=(A+BG){z}(t)+Hf(t)$$

设$\eta$是特征根：

$$\left | \eta I-A-BG \right |=0$$

可以写作：

$$(\left | \eta I-A \right |)(\left | I-\Psi(\eta)G \right |)=0$$ $$\Psi(\eta)=(\eta I-A)^{-1}B$$

因为特征值发生了变化，所以：

$$\left | \eta I-A \right | \ne0$$

因此令$\Delta(\eta)$：

$$\left |\Delta(\eta)\right |=\left | I_{2n}-\Psi(\eta)G \right |=0$$

由行列式等式：

$$\left |\Delta(\eta)\right |=\left | I_{m}-G\Psi(\eta) \right |=0$$

这样就将2n维转化为m维。

对于第i个特征值$\eta_i$，需要让$\Delta(\eta_i)$的行列式为0，只需要让$\Delta(\eta_i)$的第j列为0即可。

选择第j列，$e_j$是一个单位列向量，只有第j个值为1：

$$e_j = G\psi_j(\eta_i)$$

将选择的i个向量拼起来：

$$E=G \Psi$$

即可求解G：

$$G=E\Psi^{-1}$$

因此，为了获得相同的特征根G的选择不唯一。