基于状态空间模型的预测控制

1、状态空间模型

线性离散时间系统的状态空间模型如下：
$\begin{aligned} &x(k+1)=Ax(k)+B_uu(k)+B_dd(k)\\[1ex] &y_c(k)=C_cx(k)\tag{1} \end{aligned}$

其中 $x(k)\in\cal{R}^{n_x}$ 是状态变量； $u(k)\in\cal{R}^{n_u}$ 是控制输入变量； $y_c(k)\in\cal{R}^{n_c}$ 是被控输出变量； $d(k)\in\cal{R}^{n_d}$ 是可以测量的外部干扰变量。

上面的离散时间模型与连续时间模型
$\begin{aligned} &\dot{x}=A_cx(t)+B_{cu}u(t)+B_{cd}d(t)\\[1ex] &y_c(t)=C_cx(t)\tag{2} \end{aligned}$
之间有如下的转化关系：
$\begin{aligned} A&=e^{A_cT_s}\\[1ex] B_u&=\int_0^{T_s}e^{A_c\tau}\mathrm{d}\tau\cdot B_{cu}\\[1ex] B_d&=\int_0^{T_s}e^{A_c\tau}\mathrm{d}\tau\cdot B_{cd} \end{aligned}\tag{3}$
其中 $T_s$ 是系统的采样时间。

2、预测方程

将模型（1）改写成增量模型
$\begin{aligned} \Delta x(k+1)&=A\Delta x(k) + B_u\Delta u(k) + B_d\Delta d(k)\\[1ex] y_c(k)&=C_c\Delta x(k)+y_c(k-1)\tag{4} \end{aligned}$

其中
$\begin{aligned} \Delta x(k)=x(k)-x(k-1)\\[1ex] \Delta u(k)=u(k)-u(k-1)\\[1ex] \Delta d(k)=d(k)-d(k-1)\\[1ex] \end{aligned}$

首先以最新测量值为初始条件，基于模型（4）预测系统未来的动态。为此，设定预测时域为 $p$ ，控制时域为 $m$ 且 $m\leq p$ 。为了推导系统的预测方程，我们假设：

控制时域之外，控制量不变，即 $\Delta u(k+i)=0,\quad i=m,m+1,\cdots,p-1$ ；
该假设是因为控制时域有可能小于预测时域，而预测系统未来动态需要在整个预测时域的控制输入。
可测干扰在 $k$ 时刻之后不变，即 $\Delta d(k+i)=0,\quad i=1,2,\cdots,p-1$ ；
该假设是因为在当前时刻（ $k$ 时刻），我们还不知道干扰的未来取值。

在当前时刻 $k$ ，测量值为 $x (k)$ ，可以计算 $\Delta x(k)=x(k)-x(k-1)$ ，这个将作为预测系统未来动态的起点。由（4）可以预测 $k + 1$ 到 $k + 3$ 时刻的状态（实际上是状态增量，简单起见称为状态）如下：
$\begin{aligned} \Delta x(k+1|k)&=A\Delta x(k) + B_u\Delta u(k) + B_d\Delta d(k)\\[1ex] \Delta x(k+2|k)&=A\Delta x(k+1|k) + B_u\Delta u(k+1) + B_d\Delta d(k+1)\\[1ex] &=A^2\Delta x(k) + AB_u\Delta u(k) + B_u\Delta u(k+1) + AB_d\Delta d(k)\\[1ex] \Delta x(k+3|k)&=A\Delta x(k+2|k) + B_u\Delta u(k+2) + B_d\Delta d(k+2)\\[1ex] &=A^3\Delta x(k) + A^2B_u\Delta u(k) + AB_u\Delta u(k+1) + B_u\Delta u(k+2) + A^2B_d\Delta d(k)\tag{5} \end{aligned}$

式中， $k + 1∣ k$ 表示 $k$ 时刻对 $k + 1$ 时刻的预测；符号 $∣$ 后面的 $k$ 表示当前时刻为 $k$ 。进而，可以预测 $k + m$ 至 $k + p$ 时刻的状态
$\begin{aligned} \Delta x(k+m|k)&=A\Delta x(k+m-1|k) + B_u\Delta u(k+m-1) + B_d\Delta d(k+m-1)\\[1ex] &=A^m\Delta x(k) + A^{m-1}B_u\Delta u(k) + A^{m-2}B_u\Delta u(k+1) +\cdots + B_u\Delta u(k+m-1) + A^{m-1}B_d\Delta d(k)\\[1ex] &\vdots\\[1ex] \Delta x(k+p|k)&=A\Delta x(k+p-1|k) + B_u\Delta u(k+p-1) + B_d\Delta d(k+p-1)\\[1ex] &=A^p\Delta x(k) + A^{p-1}B_u\Delta u(k) + A^{p-2}B_u\Delta u(k+1) +\cdots + A^{p-m}B_u\Delta u(k+m-1) + A^{p-1}B_d\Delta d(k)\\[1ex] \end{aligned}\tag{6}$

进一步，由增量方程（4）可以预测 $k + 1$ 至 $k + p$ 的被控输出
$\begin{aligned} y_c(k+1|k)&=C_c\Delta x(k+1|k)+y_c(k)\\[1ex] &=C_cA\Delta x(k)+C_cB_u\Delta u(k)+C_cB_d\Delta d(k)+y_c(k),\\[1ex] y_c(k+2|k)&=C_c\Delta x(k+2|k)+y_c(k+1|k)\\[1ex] &=(C_cA^2+C_cA)\Delta x(k)+(C_cAB_u+C_cB_u)\Delta u(k)\\[1ex]&\quad+C_cB_u\Delta u(k+1)+(C_cAB_d+C_cB_d)\Delta d(k)+y_c(k),\\[1ex] &\vdots\\[1ex] y_c(k+m|k)&=C_c\Delta x(k+m|k)+y_c(k+m-1|k)\\[1ex] &=\sum_{i=1}^mC_cA^i\Delta x(k)+\sum_{i=1}^mC_cA^{i-1}B_u\Delta u(k)+\sum_{i=1}^{m-1}C_cA^{i-1}B_u\Delta u(k+1)+\cdots\\[1ex]&\quad+C_cB_u\Delta u(k+m-1)+\sum_{i=1}^mC_cA^{i-1}B_d\Delta d(k)+y_c(k),\\[1ex] &\vdots\\[1ex] y_c(k+p|k)&=C_c\Delta x(k+p|k)+y_c(k+p-1|k)\\[1ex] &=\sum_{i=1}^pC_cA^i\Delta x(k)+\sum_{i=1}^pC_cA^{i-1}B_u\Delta u(k)+\sum_{i=1}^{p-1}C_cA^{i-1}B_u\Delta u(k+1)+\cdots\\[1ex]&\quad+\sum_{i=1}^{p-m+1}C_cA^{i-1}B_u\Delta u(k+m-1)+\sum_{i=1}^pC_cA^{i-1}B_d\Delta d(k)+y_c(k),\\[1ex] \end{aligned}\tag{7}$

定义 $p$ 步预测输出向量和 $m$ 步输入向量如下：
$\begin{aligned} Y_p(k+1|k)&\overset{\mathrm{def}}{=\!=}\left[ \begin{matrix} y_c(k+1|k) \\[2ex] y_c(k+2|k) \\[2ex] \vdots\\[2ex] y_c(k+p|k) \end{matrix} \right]_{p\times 1},\\[1ex] \Delta U(k)&\overset{\mathrm{def}}{=\!=}\left[ \begin{matrix} \Delta u(k) \\[2ex] \Delta u(k+1) \\[2ex] \vdots\\[2ex] \Delta u(k+m-1) \end{matrix} \right]_{m\times 1} \end{aligned}\tag{8}$

那么，对系统未来 $p$ 步预测的输出可以由下面的预测方程计算：
$Y_p(k+1|k)=S_x\Delta x(k)+{\cal{I}}y_c(k)+{\cal{S}_d}\Delta d(k)+{\cal{S}_u}\Delta U(k)\tag{9}$

其中
$\begin{aligned} S_x&=\left[ \begin{matrix} C_cA \\[2ex] \sum_{i=1}^2C_cA^i \\[2ex] \vdots\\[2ex] \sum_{i=1}^pC_cA^i \end{matrix} \right]_{p\times 1} ,\quad{\cal{I}}=\left[ \begin{matrix} I_{n_c\times n_c} \\[1ex] I_{n_c\times n_c} \\[1ex] \vdots\\[1ex] I_{n_c\times n_c} \end{matrix} \right]_{p\times 1},\quad{\cal{S}_d}=\left[ \begin{matrix} C_cB_d \\[2ex] \sum_{i=1}^2C_cA^{i-1}B_d \\[2ex] \vdots\\[2ex] \sum_{i=1}^pC_cA^{i-1}B_d \end{matrix} \right]_{p\times 1},\\[4ex] {\cal{S}_u}&=\left[ \begin{matrix} C_cB_u & 0 & 0 & \cdots & 0\\[2ex] \sum_{i=1}^2C_cA^iB_u & C_cB_u & 0 & \cdots & 0 \\[2ex] \vdots & \vdots & \vdots & &\vdots \\[2ex] \sum_{i=1}^mC_cA^{i-1}B_u & \sum_{i=1}^{m-1}C_cA^{i-1}B_u & \cdots & \cdots & C_cB_u \\[2ex] \vdots & \vdots & \vdots & &\vdots \\[2ex] \sum_{i=1}^pC_cA^{i-1}B_u & \sum_{i=1}^{p-1}C_cA^{i-1}B_u & \cdots & \cdots & \sum_{i=1}^{p-m+1}C_cA^{i-1}B_u \end{matrix} \right]_{p\times m} \end{aligned}\tag{10}$