1. MPC 算法简介
MPC（Model Predictive Control，模型预测控制）算法是一种进阶过程控制方法。本质上是一个二次优化的过程。过程内环动态模型、控制量的历史数值、在预测区间上的一个最优值方程 J；最优控制量可由以上各量求出。MPC 可应用于线性和非线性系统。
2. 系统描述
线性离散系统 $$ X(k+1) = A X(k) + B u(k) $$ $$ y(k) = C X(k) $$ 增量形式 $$ \Delta X(k+1) = A \Delta X(k) + B \Delta u(k) $$ $$ y(k) = C \Delta X(k) + y(k-1) $$ 其中 $$ \Delta X(k) = X(k) - X(k-1) $$ $$ \Delta u(k) = u(k) - u(k-1) $$
3. 预测方程
预测时域 p，控制时域 m，当 t > m 时 $\Delta u(k+i) = 0$；定义 $$ Y_p(k+1|k) = \begin{bmatrix} y(k+1|k) \ y(k+2|k) \ \vdots \ y(k+p|k) \end{bmatrix}{p \times 1}, \quad \Delta U(k) = \begin{bmatrix} \Delta u(k) \ \Delta u(k+1) \ \vdots \ \Delta u(k+m) \end{bmatrix}{m \times 1} $$ 则有，预测方程为 $$ Y_p(k+1|k) = S_X \Delta X(k) + I y(k) + S_u \Delta U(k) $$ 其中 $$ S_X = \begin{bmatrix} C A \ \sum_{i=1}^{2} C A^i \ \vdots \ \sum_{i=1}^{p} C A^i \end{bmatrix}{p \times 1}, \quad I = \begin{bmatrix} I{nc \times nc} \ I_{nc \times nc} \ \vdots \ I_{nc \times nc} \end{bmatrix}{p \times 1} $$ $$ S_u = \begin{bmatrix} C B & 0 & 0 \ \sum{i=1}^{2} C A^{i-1} B & C B & 0 \ \vdots & \vdots & \vdots \ \sum_{i=1}^{p} C B & \cdots & \sum_{i=1}^{p+m-1} C A^{i-1} B \end{bmatrix}_{p \times m} $$
4. 预测控制
4.1 定义目标函数
$$ J = \sum_{i=1}^{p} | \Gamma_{yi} (y(k+i|k) - r(k+i)) |^2 + \sum_{i=1}^{m} | \Gamma_{ui} \Delta u(k+i-1) |^2 $$ 其中 $\Gamma_{yi} = \operatorname{diag}(\Gamma_{y1,i}, \dots, \Gamma_{yp,i})$；$\Gamma_{ui} = \operatorname{diag}(\Gamma_{u1,i}, \dots, \Gamma_{um,i})$。 矩阵化： $$ J(X(k), \Delta U(k), m, p) = | \Gamma_y (Y_P(k+1|k) - R(k+1)) |^2 + | \Gamma_u \Delta U(k) |^2 $$ 其中 $$ \Gamma_y = \operatorname{diag}(\Gamma_{y1}, \dots, \Gamma_{yp}), \quad \Gamma_u = \operatorname{diag}(\Gamma_{u1}, \dots, \Gamma_{um}) $$ $$ R(k+1) = \begin{bmatrix} r(k+1) \ r(k+2) \ \vdots \ r(k+p) \end{bmatrix}_{p \times 1} $$ $\Gamma_y$ 表征调节对误差的追踪能力；$\Gamma_u$ 表征调节运行的稳定性；$R(k+1)$ 表征参考路径，采用误差演化模型时取零即可，这是优点。
4.2 求解二次优化
$$ \Delta U^*(k) = (S_u^T \Gamma_y^T \Gamma_y S_u + \Gamma_u^T \Gamma_u)^{-1} S_u^T \Gamma_y^T \Gamma_y E_p(k+1|k) $$ 其中 $E_p(k+1|k) = R(k+1) - S_X \Delta X(k) - I y(k)$。
4.3 取第一个元素作用于真实系统
$$ \Delta u(k) = [I_{nc \times nc}, 0, 0, \dots]{1 \times m} * \Delta U^* $$ 定义预测控制增益系数 $K{mpc} = [I_{nc \times nc}, 0, 0, \dots] (S_u^T \Gamma_y^T \Gamma_y S_u + \Gamma_u^T \Gamma_u)^{-1} S_u^T \Gamma_y^T \Gamma_y$，则 $\Delta u(k) = K_{mpc} * E_p(k+1|k)$。
5. 算法流程
Step 1：初始化，设定预测时域 p 和控制时域 m，令 $u(-1)=0, X(-1)=0$，计算由上述公式计算 $S_X$、$I$、$S_u$ 和 $K_{mpc}$； Step 2：$k \ge 0$ 时，测量 $X(k)$，计算 $y(k)$，计算 $\Delta X(k) = X(k) - X(k-1)$； Step 3：计算 $E_p(k+1|k) = R(k+1) - S_X \Delta X(k) - I y(k)$； Step 4：计算 $\Delta u(k) = K_{mpc} * E_p(k+1|k)$； Step 5：将 $u(k) = u(k-1) + \Delta u(k)$ 作用于系统； Step 6：在 $k+1$ 时刻，测量 $X(k+1)$ 的值，并且令 $k=k+1$，返回 Step 2； 注：如果采用时变加权因子，则 $K_{mpc}$ 的计算或部分计算也将在循环内进行。
