1、约束优化问题描述

考虑如下线性离散时间系统的状态空间增量模型：
$$\begin{array}{rl}\Delta x(k+1)& =A\Delta x(k)+B_u\Delta u(k)+B_d\Delta d(k)\\ y_c(k)& =C_c\Delta x(k)+y_c(k-1)\\ y_b(k)& =C_b\Delta x(k)+y_b(k-1)\end{array}\text{\hspace{1em}(1)}$$

其中
$$\begin{array}{rl}\Delta x(k)& =x(k)-x(k-1)\\ \Delta u(k)& =u(k)-u(k-1)\\ \Delta d(k)& =d(k)-d(k-1)\end{array}\text{\hspace{1em}(2)}$$

模型中 $\Delta x(k)$ 是状态增量；$\Delta u(k)$ 是控制输入增量；$\Delta d(k)$ 是可以测量的外部干扰增量；$y_c(k)$ 是被控输出量；$y_b(k)$ 是约束输出量；$A,B_u,B_d,C_c,C_b$ 是相应维数的系统矩阵。

我们考虑的控制目标是使得被控输出 $y_c$ 跟踪给定的参考输入 $r$。同时，系统的控制量，控制增量和输出量满足下面的控制约束及输出约束：
$$\begin{array}{rl}{u}_{\min }(k)& \le u(k)\le u_{\max }(k),\forall k\ge 0,\\ \Delta u_{\min }(k)& \le \Delta u(k)\le \Delta u_{\max }(k),\forall k\ge 0,\\ y_{\min }(k)& \le y_b(k)\le y_{\max }(k),\forall k\ge 0\end{array}\text{\hspace{1em}(3)}$$

简单起见，我们假设系统的全部状态是可以测量的。如果系统的状态不是全部可以测量的，则需要设计状态观测器，以估计状态作为模型的初始状态，预测系统未来动态。

在 $k$ 时刻有状态测量（估计）值 $x(k)$，根据预测控制的基本原理，约束 MPC 的优化问题描述为

$$\underset{\Delta U(k)}{\min }J(x(k),\Delta U(k))\text{\hspace{1em}(4)}$$
满足系统动力学（$i=0,1,\cdots ,p$）
$$\begin{array}{rl}\Delta x(k+i+1|k)& =A\Delta x(k+i|k)+B_u\Delta u(k+i)+B_d\Delta d(k+i)\\ \Delta x(k|k)& =\Delta x(k)\\ y_c(k+i|k)& =C_c\Delta x(k+i|k)+y_c(k+i-1|k),i\ge 1\\ y_c(k|k)& =y_c(k)\\ y_b(k+i|k)& =C_b\Delta x(k+i|k)+y_b(k+i-1|k),i\ge 1\\ y_b(k|k)& =y_b(k)\end{array}\text{\hspace{1em}(5)}$$

及时域约束
$$\begin{array}{rl}& u_{\min }(k+i)\le u(k+i)\le u_{\max }(k+i),i=0,1,\cdots ,m-1,\\ & \Delta u_{\min }(k+i)\le \Delta u(k+i)\le \Delta u_{\max }(k+i),i=0,1,\cdots ,m-1,\\ & y_{\min }(k+i)\le y_b(k+i)\le y_{\max }(k+i),i=0,1,\cdots ,p\end{array}\text{\hspace{1em}(6)}$$

其中
$$J(x(k),\Delta U(k))=\parallel {\Gamma }_y(Y_{p,c}(k+1|k)-R(k+1)){\parallel }^2+\parallel {\Gamma }_u\Delta U(k){\parallel }^2\text{\hspace{1em}(7)}$$

上面的优化问题中 ${\Gamma }_y,{\Gamma }_u$ 是加权矩阵，给定为
Γ y = d i a g { Γ y , 1 , Γ y , 2 , ⋯   , Γ y , p } p × p , Γ u = d i a g { Γ u , 1 , Γ u , 2 , ⋯   , Γ u , p } m × m (8) \begin{aligned} &\Gamma_y=\mathrm{diag}\lbrace\Gamma_{y,1},\Gamma_{y,2},\cdots,\Gamma_{y,p}\rbrace_{p\times p},\\[1ex] &\Gamma_u=\mathrm{diag}\lbrace\Gamma_{u,1},\Gamma_{u,2},\cdots,\Gamma_{u,p}\rbrace_{m\times m}\\[1ex] \end{aligned}\tag{8} ​Γy​=diag{Γy,1​,Γy,2​,⋯,Γy,p​}p×p​,Γu​=diag{Γu,1​,Γu,2​,⋯,Γu,p​}m×m​​(8)

$R(k+1)$ 是给定的控制输入参考序列，为
$$R(k+1)={\left[\begin{array}{c}r(k+1)\\ r(k+2)\\ ⋮\\ r(k+p)\end{array}\right]}_{p\times 1}\text{\hspace{1em}(9)}$$

$\Delta U(k)$ 是控制量增量序列，作为约束优化问题的独立变量，定义为
$$\Delta U(k)\stackrel{\mathrm{def}}{==}{\left[\begin{array}{c}\Delta u(k)\\ \Delta u(k+1)\\ ⋮\\ \Delta u(k+m-1)\end{array}\right]}_{m\times 1}\text{\hspace{1em}(10)}$$

$Y_{p,c}(k+1|k)$ 是 $k$ 时刻基于模型（1）预测的 $p$ 步控制输出，定义为
$$Y_{p,c}(k+1|k)\stackrel{\mathrm{def}}{==}{\left[\begin{array}{c}{y}_c(k+1|k)\\ y_c(k+2|k)\\ ⋮\\ y_c(k+p|k)\end{array}\right]}_{p\times 1}\text{\hspace{1em}(11)}$$

具体的，预测的控制输出 $y_c(k+i|k)$ 和约束输出 $y_b(k+i|k)$ 由（5）计算。

进一步，$Y_{p,c}(k+1|k)$ 可以由预测方程
$$Y_p(k+1|k)=S_x\Delta x(k)+\mathcal{I}{y}_c(k)+{\mathcal{S}}_u\Delta U(k)+{\mathcal{S}}_d\Delta d(k)\text{\hspace{1em}(12)}$$

计算，其中
$$\begin{array}{rl}{S}_x& ={\left[\begin{array}{c}{C}_{c}A\\ {\sum }_{i=1}^2{C}_c{A}^i\\ ⋮\\ {\sum }_{i=1}^p{C}_c{A}^i\end{array}\right]}_{p\times 1},\mathcal{I}={\left[\begin{array}{c}{I}_{n_c\times n_c}\\ I_{n_c\times n_c}\\ ⋮\\ I_{n_c\times n_c}\end{array}\right]}_{p\times 1},{\mathcal{S}}_d={\left[\begin{array}{c}{C}_c{B}_d\\ {\sum }_{i=1}^2{C}_c{A}^{i-1}{B}_d\\ ⋮\\ {\sum }_{i=1}^p{C}_c{A}^{i-1}{B}_d\end{array}\right]}_{p\times 1},\\ {\mathcal{S}}_u& ={\left[\begin{array}{ccccc}{C}_c{B}_u& 0& 0& \cdots & 0\\ {\sum }_{i=1}^2{C}_c{A}^{i-1}{B}_u& C_c{B}_u& 0& \cdots & 0\\ ⋮& ⋮& ⋮& & ⋮\\ {\sum }_{i=1}^m{C}_c{A}^{i-1}{B}_u& {\sum }_{i=1}^{m-1}{C}_c{A}^{i-1}{B}_u& \cdots & \cdots & C_c{B}_u\\ ⋮& ⋮& ⋮& & ⋮\\ {\sum }_{i=1}^p{C}_c{A}^{i-1}{B}_u& {\sum }_{i=1}^{p-1}{C}_c{A}^{i-1}{B}_u& \cdots & \cdots & {\sum }_{i=1}^{p-m+1}{C}_c{A}^{i-1}{B}_u\end{array}\right]}_{p\times m}\end{array}\text{\hspace{1em}(13)}$$

2、约束优化问题求解

由于约束条件（6）的存在，一般情况下，我们无法得到优化问题的解析解，因此，需要采用数值求解方法。由于目标函数是二次型的，动力学方程和时域约束条件是线性的，所以该问题是一个二次规划（Quadratic Programming，QP）问题。

QP 问题是下面形式的优化问题：
$$\underset{z}{\min }{z}^{T}Hz-g^{T}z,满足Cz\ge b\text{\hspace{1em}(14)}$$

其中 $H$ 是 $Hessian$ 矩阵，

2.1、目标函数转化为 $z^{T}Hz-g^{T}z$

$z=\Delta U(k)$ 是优化问题的独立变量。将预测方程（12）代入目标函数（7），并定义
$$E_p(k+1|k)\stackrel{\mathrm{def}}{==}R(k+1)-{\mathcal{S}}_x\Delta x(k)-\mathcal{I}{y}_c(k)-{\mathcal{S}}_d\Delta d(k)\text{\hspace{1em}(15)}$$

则目标函数变为
$$\begin{array}{rl}J(x(k),\Delta U(k))& =\parallel {\Gamma }_y(Y_{p,c}(k+1|k)-R(k+1)){\parallel }^2+\parallel {\Gamma }_u\Delta U(k){\parallel }^2\\ & =\parallel {\Gamma }_y({\mathcal{S}}_u\Delta U(k)-E_p(k+1|k)){\parallel }^2+\parallel {\Gamma }_u\Delta U(k){\parallel }^2\\ & =\Delta U(k{)}^T{\mathcal{S}}_u^T{\Gamma }_y^T{\Gamma }_y{\mathcal{S}}_u\Delta U(k)+\Delta U(k{)}^T{\Gamma }_u^T{\Gamma }_u\Delta U(k)\\ & -2E_p(k+1|k{)}^T{\Gamma }_y^T{\Gamma }_y{\mathcal{S}}_u\Delta U(k)+E_p(k+1|k{)}^T{\Gamma }_y^T{\Gamma }_y{E}_p(k+1|k)\end{array}\text{\hspace{1em}(16)}$$

因为 $E_p(k+1|k{)}^T{\Gamma }_y^T{\Gamma }_y{E}_p(k+1|k)$ 与独立变量 $\Delta U(k)$ 无关，所以对优化问题而言，上式等价于
$$\boxed{\tilde{J}=\Delta U(k{)}^{T}H\Delta U(k)-G(k+1|k{)}^T\Delta U(k)},\text{\hspace{1em}(17)}$$

其中
$$\begin{array}{rl}H& ={\mathcal{S}}_u^T{\Gamma }_y^T{\Gamma }_y{\mathcal{S}}_u+{\Gamma }_u^T{\Gamma }_u,\\ G(k+1|k)& =2{\mathcal{S}}_u^T{\Gamma }_y^T{\Gamma }_y{E}_p(k+1|k)\end{array}\text{\hspace{1em}(18)}$$

2.2、控制增量约束转化为 $Cz\ge b$

可以将控制量增量约束直接写成需要的形式
$$\left[\begin{array}{c}-T\\ T\end{array}\right]\Delta U(k)\ge \left[\begin{array}{c}-\Delta u_{\max }(k)\\ ⋮\\ -\Delta u_{\max }(k+m-1)\\ \Delta u_{\min }(k)\\ ⋮\\ \Delta u_{\min }(k+m-1)\end{array}\right],\text{\hspace{1em}(19)}$$

其中
$$T={\left[\begin{array}{cccc}{I}_{n_u\times n_u}& 0& \cdots & 0\\ 0& I_{n_u\times n_u}& \cdots & 0\\ ⋮& ⋮& & ⋮\\ 0& 0& \cdots & I_{n_u\times n_u}\end{array}\right]}_{m\times m}\text{\hspace{1em}(20)}$$

2.3、控制量约束转化为 $Cz\ge b$

由于 $$\Delta u(k+i)=u(k+i)-u(k+i-1)\text{\hspace{1em}(21)}$$

因此对 $i=0$ 有 $u(k)=u(k-1)+\Delta u(k)$，进而有
$$u_{\min }(k)\le u(k)\le u_{\max }(k)⟹\{\begin{array}{l}-\Delta u(k)\ge u(k-1)-u_{\max }(k)\\ \Delta u(k)\ge u_{\min }(k)-u(k-1)\end{array}\text{\hspace{1em}(22)}$$

对 $i=1$ 有 $u(k+1)=u(k-1)+\Delta u(k)+\Delta u(k+1)$，因此有
$$u_{\min }(k+1)\le u(k+1)\le u_{\max }(k+1)⟹\{\begin{array}{l}-(\Delta u(k)+\Delta u(k+1))\ge u(k-1)-u_{\max }(k+1)\\ \Delta u(k)+\Delta u(k+1)\ge u_{\min }(k+1)-u(k-1)\end{array}\text{\hspace{1em}(23)}$$

对任意 $i=0,1,2,\cdots ,m-1$ 有 $u(k+i)=u(k-1)+{\sum }_{j=0}^i\Delta u(k+j)$，因此有
$$u_{\min }(k+i)\le u(k+i)\le u_{\max }(k+i)⟹\{\begin{array}{l}-{\sum }_{j=0}^i\Delta u(k+j)\ge u(k-1)-u_{\max }(k+i)\\ {\sum }_{j=0}^i\Delta u(k+j)\ge u_{\min }(k+i)-u(k-1)\end{array}\text{\hspace{1em}(24)}$$

综上，将控制量约束写成如下的矩阵形式：
$$\left[\begin{array}{c}-L\\ L\end{array}\right]\Delta U(k)\ge \left[\begin{array}{c}u(k-1)-u_{\max }(k)\\ ⋮\\ u(k-1)-u_{\max }(k+m-1)\\ u_{\min }(k)-u(k-1)\\ ⋮\\ u_{\min }(k+m-1)-u(k-1)\end{array}\right],\text{\hspace{1em}(25)}$$

其中
$$L={\left[\begin{array}{cccc}{I}_{n_u\times n_u}& 0& \cdots & 0\\ I_{n_u\times n_u}& I_{n_u\times n_u}& \cdots & 0\\ ⋮& ⋮& & ⋮\\ I_{n_u\times n_u}& I_{n_u\times n_u}& \cdots & I_{n_u\times n_u}\end{array}\right]}_{m\times m}\text{\hspace{1em}(26)}$$

2.4、输出约束转化为 $Cz\ge b$

对系统未来 $p$ 步约束输出的预测可以由下面的预测方程计算：
$$Y_{p,b}(k+1|k)={\mathcal{S}}_{x,b}\Delta x(k)+{\mathcal{I}}_b{y}_b(k)+{\mathcal{S}}_{u,b}\Delta U(k)+{\mathcal{S}}_{d,b}\Delta d(k)\text{\hspace{1em}(27)}$$

其中，${\mathcal{S}}_{x,b},{\mathcal{I}}_b,{\mathcal{S}}_{d,b},{\mathcal{S}}_{u,b}$ 与式（13）类似，只是将 $C_c$ 更换为 $C_b$。

记
$$Y_{\min }(k+1)={\left[\begin{array}{c}{y}_{\min }(k+1)\\ y_{\min }(k+2)\\ ⋮\\ y_{\min }(k+p)\end{array}\right]}_{p\times 1},Y_{\max }(k+1)={\left[\begin{array}{c}{y}_{\max }(k+1)\\ y_{\max }(k+2)\\ ⋮\\ y_{\max }(k+p)\end{array}\right]}_{p\times 1}\text{\hspace{1em}(28)}$$

可以将输出约束描述为如下的向量形式：
$$Y_{\min }(k+1)\le Y_{p,b}(k+1|k)\le Y_{\max }(k+1)\text{\hspace{1em}(29)}$$

将约束输出预测方程（27）代入上式，则输出约束转换为
$$\begin{array}{r}\left[\begin{array}{c}-{\mathcal{S}}_{u,b}\\ {\mathcal{S}}_{u,b}\end{array}\right]\Delta U(k)\ge \left[\begin{array}{c}({\mathcal{S}}_{x,b}\Delta x(k)+{\mathcal{I}}_b{y}_b(k)\Delta U(k)+{\mathcal{S}}_{d,b}\Delta d(k))-Y_{\max }(k+1)\\ -({\mathcal{S}}_{x,b}\Delta x(k)+{\mathcal{I}}_b{y}_b(k)\Delta U(k)+{\mathcal{S}}_{d,b}\Delta d(k))+Y_{\min }(k+1)\end{array}\right]\end{array}\text{\hspace{1em}(30)}$$

综上所述，约束 MPC 的优化问题可以转换为如下的 QP 问题描述：
$$\begin{array}{rl}& \underset{\Delta U(k)}{\min }\Delta U(k{)}^{T}H\Delta U(k)-G(k+1|k{)}^T\Delta U(k),\\ & 满足C_u\Delta U(k)\ge b(k+1|k)\end{array}$$

其中 $H,G(k+1|k)$ 由（18）给出，
$$\begin{array}{rl}{C}_u& ={\left[\begin{array}{cccccc}-T^T& T^T& -L^T& L^T& -{\mathcal{S}}_{u,b}^T& {\mathcal{S}}_{u,b}^T\end{array}\right]}_{(4m+2p)\times 1}^T,\\ b(k+1|k)& ={\left[\begin{array}{c}-\Delta u_{\max }(k)\\ ⋮\\ -\Delta u_{\max }(k+m-1)\\ \Delta u_{\min }(k)\\ ⋮\\ \Delta u_{\min }(k+m-1)\\ u(k-1)-u_{\max }(k)\\ ⋮\\ u(k-1)-u_{\max }(k+m-1)\\ u_{\min }(k)-u(k-1)\\ ⋮\\ u_{\min }(k+m-1)-u(k-1)\\ ({\mathcal{S}}_{x,b}\Delta x(k)+{\mathcal{I}}_b{y}_b(k)\Delta U(k)+{\mathcal{S}}_{d,b}\Delta d(k))-Y_{\max }(k+1)\\ -({\mathcal{S}}_{x,b}\Delta x(k)+{\mathcal{I}}_b{y}_b(k)\Delta U(k)+{\mathcal{S}}_{d,b}\Delta d(k))+Y_{\min }(k+1)\end{array}\right]}_{(4m+2p)\times 1}\end{array}$$

3、约束 MPC 的闭环解

由（18）知 $H\ge 0$，因此 QP 问题对任何加权矩阵 ${\Gamma }_y\ge 0,{\Gamma }_u\ge 0$ 有解，记为 $\Delta U^{\ast }(k)$。根据预测控制的基本原理，得到的开环控制序列的第一步作用于被控系统。在下一个采样时刻，将用新的测量值刷新约束优化问题，即 QP 问题，并重新求解。因此，约束 MPC 的闭环控制率定义为
$$\Delta u(k)={\left[\begin{array}{cccc}{I}_{n_u\times n_u}& 0& \cdots & 0\end{array}\right]}_{1\times m}\Delta U^{\ast }(k)$$

下面给出了约束 MPC 控制器的实现流程。