模型嵌入控制（2）：被控对象建模

被控对象以可控动态为核心，进一步引入了扰动动态，其数值实现称为嵌入模型（Embedded Model），是模型嵌入控制设计的基础。本文针对线性时不变（LTI：Linear Time-Invariant）系统建模，给出可控动态和扰动动态离散时间的状态空间表述。

可控动态建模

可控动态是对被控对象物理特性的描述，我们之前在系统的描述中讨论过线性时不变系统的三种描述方法，在此简要复述为：

• 微分方程：被控对象的特性通常可以根据一定的物理原理建立微分方程，这是描述系统特性最基本的数学手段；
• 传递函数（或传递函数矩阵）：系统微分方程的拉普拉斯变换，直接地给出了输入输出关系，是频域设计和分析的依据；
• 状态空间方程：微分方程的等价描述，由于进一步定义了系统状态而能够对系统进行更加精细的描述。离散形式的状态空间还有给出了数值系统的具体实现。

我们假设读者已经具备连续时间系统的建模能力，设被控对象在连续时间下的状态空间微分方程为：

$$ \left\{\begin{aligned} \dot{\bm{x}} &= A \bm{x} + B \bm{u} \\ \bm{y} &= C \bm{x} + D \bm{u} \end{aligned}\right. $$

其中 $\bm{x}$ 为系统的状态向量，$\bm{y}$ 为输出，$A$、$B$、$C$、$D$ 分别为状态矩阵、输入矩阵、输出矩阵和直馈矩阵。通常情况下 $D=0$，在后面的讨论中将省略。

可控动态是对被控对象在离散时间下的建模，换言之，需要将上述连续时间的状态空间方程依据采样时间 $T$ 离散化。在一阶近似下，状态微分可以由微分代替，为：

$$ \dot{\bm{x}} \approx \frac{1}{T} \bigl( \bm{x}_c(iT+T) - \bm{x}_c(iT) \bigr) $$

记 $\bm{x}_c$ 为可控动态的状态变量。为了表述的方便，后续离散时间状态变量的索引中将省略采样时间 $T$，改写为 $\bm{x}_c(iT) \rightarrow \bm{x}_c(i)$。如此做，可控动态的状态空间方程为：

$$ \left\{\begin{aligned} & \bm{x}_c(i+1) = A_c \bm{x}_c(i) + B_c \Bigl( \bm{u}(i) + \bm{h}(\bm{x}_c(i)) \Bigr) + \bm{d}(i) \\ & \bm{y}(i) = C_c \bm{x}_c(i) + C_d \bm{x}_d(i) \end{aligned}\right. $$

其中，$A_c = I+AT$，$B_c = BT$，$C_c = C$。为了更加准确地对被控对象进行描述，可控动态补充了 $\bm{h}(\cdot)$ 以隐式地描述状态之间可能的未知或非线性耦合和输入端扰动 $\bm{d}$。在一定的情况下（通常是多个传感器同时测量时），状态预测器能够对被控对象输出端的扰动（如传感器零偏）进行预测，因此不失一般性地可以在输出端增加 $C_d \bm{x}_d$ 以表述该扰动的影响，后面将通过控制律设计对这部分扰动进行补偿。

在连续时间系统的离散化中，如果一阶近似带来较大误差，严格的建模应当从微分方程通解的积分中获得。对于线性时不变系统，状态方程的通解为：

$$ \bm{x}(t) = \mathrm{e}^{A(t-t_0)} \bm{x}(t_0) + \int_{t_0}^t \mathrm{e}^{A(t-\tau)} B \bm{u}(\tau) \,\mathrm{d} \tau $$

取 $t_0 = iT$，$t = (i+1)T$，并设在 $iT \le t < (i+1)T$ 时间内指令 $\bm{u}(t)$ 保持不变，因此离散后的矩阵为：

$$ A_c = \mathrm{e}^{AT} ,\quad B_c = \int_0^T \mathrm{e}^{A\tau} B \, \mathrm{d} \tau $$

扰动动态建模

由于扰动具有一定的随机性，扰动动态通常依赖于扰动的统计特性。将扰动看作随机噪声 $\bm{w}$ 经过扰动动态的输出，于是从频域上可以将扰动动态看作噪声滤波器，其幅频响应与扰动的功率谱密度一致。例如通过某种手段测得外部扰动的功率谱密度呈现 $1/f$ 的特性，意味着扰动动态可以建模为一阶积分器。

此外，如果可预测部分的扰动具有较强的时域特性，例如零偏或者震荡等，可以将扰动在局部时间内看作某个微分方程的解，因此总体扰动可以看作该微分方程的通解。例如对于未知的常值扰动，可以看作微分方程 $\dot{d} = 0$ 的某一个解，根据微分方程可以将扰动动态建模为一阶积分器。

不失一般性，扰动动态的状态空间方程可以写为：

$$ \left\{\begin{aligned} &\bm{x}_d(i+1) = A_d \bm{x}_d(i) + G_d \bm{w}_d(i) \\ &\bm{d}(i) = H_c \bm{x}_d(i) + G_c \bm{w}_c(i) \end{aligned}\right. $$

式中 $\bm{x}_d$ 为扰动动态的状态向量。为了对扰动进行充分的预测，扰动动态的极点（$A_d$ 的特征值）应当分布在单位圆上，这存在两种基本形式：

• 极点为 $1$：串联积分形式，频域上可以解释为低频扰动的预测，时域上根据积分阶数可以解释为分段常值（一阶积分）、分段线性（二阶积分）扰动等；
• 共轭极点 $\mathrm{e}^{j 2\pi f_0 T}, \,(f_0 < f_s/2)$：震荡模型，对应呈现周期特性的扰动。

需要说明的是，随机噪声 $\bm{w}$ 分为扰动动态的驱动噪声 $\bm{w}_d$ 和直接作用于可控动态的噪声 $\bm{w}_c$。其中，$\bm{w}_d$ 驱动扰动动态，基于模型的输出部分 $H_c \bm{x}_d$ 对应可预测部分的扰动，而 $\bm{w}_c$ 完全随机而无法预测。

为了进一步解释可预测扰动和不可预测扰动，我们可以考察一个随机游走的例子：考虑有一个正态分布随机数生成器，每生成一次随机数，一个质点就会沿着数轴移动相应的距离。该随机过程可以由 $x(i+1) = x(i) + w(i)$ 描述。假设现在这个质点的坐标值为 $7$，通过数学期望可以得到下一次移动的距离平均为 $0$，当需要对位置进行反馈时，这个值不会提供任何有用的信息，因而说“移动的距离不可预测”；然而，由于位移的均值为 $0$，可以预测下一时刻质点最可能的位置坐标将保持在 $7$，即“下一时刻的位置是可以预测的”。

假设驱动噪声 $\bm{w} = [\bm{w}_c,\,\bm{w}_d]$ 的各个分量相互独立，对于扰动动态，其状态 $\bm{x}_d$ 不失一般性地都会受到 $\bm{w}_d$ 的影响。设计输入扰动 $\bm{d}$ 时需要考虑扰动的注入点，换句话说，被控对象的哪些状态会受到扰动的影响？因此 $G_c$ 的设计应当根据可控动态的物理条件决定。

示例：二阶系统建模

我们以单自由度弹簧-质量系统为例展示可控动态和扰动动态的建模。设弹簧的特征角频率为 $\omega_0$，在控制加速度 $a_u$ 和扰动加速度 $a_d$ 的作用下，质点相对平衡位置的位移 $z$ 的微分方程为：

$$ \ddot{z} + \omega_0^2 z = a_u + a_d $$

选取连续时间的状态变量为 $\bm{x} = [z, \dot{z}]$，连续时间的状态空间方程可以写为：

$$ \left\{\begin{aligned} \dot{\bm{x}} &= \begin{bmatrix} 0 & 1 \\ -\omega_0^2 & 0 \end{bmatrix} \bm{x} + \begin{bmatrix} 0 \\ 1 \end{bmatrix} \left( a_u + a_d \right) \\ y &= \begin{bmatrix} 1 & 0 \end{bmatrix} \bm{x} \end{aligned}\right. ,\quad \bm{x}(0) = \bm{x}_0 $$

如果使用差分近似微分，可控动态的各个系数矩阵可以写为：

$$ A_c = I+AT = \begin{bmatrix} 1 & T \\ -\omega_0^2 T & 1 \end{bmatrix} , \quad B_c = BT = \begin{bmatrix} 0 \\ T \end{bmatrix} $$

如果使用积分做准确计算，可控动态的系数矩阵为：

$$ \begin{gathered} A_c = \mathrm{e}^{AT} = \begin{bmatrix} \cos(\omega_0T) & \sin(\omega_0T)/\omega_0 \\ -\omega_0 \sin(\omega_0T) & \cos(\omega_0T) \end{bmatrix} \\ B_c = \int_0^T \mathrm{e}^{A\tau} B \, \mathrm{d} \tau = \begin{bmatrix} 2 \sin^2(\omega_0T/2) /\omega_0^2 \\ \sin(\omega_0T)/\omega_0 \end{bmatrix} \end{gathered} $$

当 $\omega_0T \ll 1$ 且忽略 $T^2$ 项后，两种方法得到的可控动态一致。

实际中，我们希望将数字积分器改为累加器以节约时间步长 $T$ 引入的乘法计算。因此，取可控动态的状态向量为 $\bm{x}_c = [z, \dot{z}T]$，可控动态变为：

$$ \left\{\begin{aligned} & \bm{x}_c(i+1) = \begin{bmatrix} 1 & 1 \\ -\omega_0^2 T^2 & 1 \end{bmatrix} \bm{x}_c(i) + \begin{bmatrix} 0 \\ 1 \end{bmatrix} u + \bm{d}(i) \\ & y(i) = \begin{bmatrix} 1 & 0 \end{bmatrix} \bm{x}_c(i) \end{aligned}\right. ,\quad \bm{x}_c(0) = \bm{x}_{c0} $$

其中 $u = T^2 a_u$ 为数字形式的指令加速度。

假设质点受到的外部扰动时缓变的信号，可以认为具有分段线性特性，因而扰动动态可以取二阶积分器（实际上是累加器），其状态空间方程为：

$$ \left\{\begin{aligned} & \bm{x}_d(i+1) = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix} \bm{x}_d(i) + \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \bm{w}_d(i) \\ & \bm{d}(i) = \begin{bmatrix} 0 & 0 \\ 1 & 0 \end{bmatrix} \bm{x}_d(i) + \begin{bmatrix} 0 \\ 1 \end{bmatrix} w_c(i) \end{aligned}\right. $$

其中，$\bm{w}_d$ 的系数矩阵为单位矩阵，意味着驱动噪声 $\bm{w}_d$ 的两个分量可以分别对扰动状态 $\bm{x}_d$ 的两个分量产生影响；$\bm{d}$ 的第一个分量为零，是因为外部扰动时加速度形式，不会直接对质点的速度产生影响；噪声 $w_c$ 表示不可预测的随机扰动加速度。