如何理解MPC模型预测控制理论

MPC的基本原理与特点

基本原理

其实MPC的思想很直白也很简单：在每个采样时刻，用测量到的信息和一个模型，求解一个优化问题，并将优化出来的结果的第一个量作用到系统中，之后如此反复。总的来说分为三步：

• 预测系统未来动态
• (数值)求解优化问题
• 将优化解的第一个元素作用于系统

同时这里有几个点需要思考：

• 模型的形式是什么
• 模型为什么会有预测功能,如何体现的
• 优化问题如何建立以及怎么求解

这些问题在推导过程中会一一体现。

特点

1.基于模型

在MPC中，需要一个描述对象动态行为的模型，作用是用来根据当前的测量值预测未来的动态。但是需要注意的是，我们更注重模型的作用，即预测，而不是模型的具体形式。至于这个模型是如何构建的并不重要，只要它具有根据当前状态和未来控制输入预测未来输出的能力就可以。比如可以使用状态空间模型，传递函数模型，模糊模型等等。

2.滚动优化

简单来说就是如下的表格：

时刻	预测时间区间
k	$[k+1，k+p]系统动态$
k+1	$[k+2，k+p+1]系统动态$
k+2	$[k+3，k+p+2]系统动态$

可以看出，这种优化是在线进行的，虽然不是全局最优，但是这种滚动的处理，可以将模型失配、时变、干扰等不确定考虑在里面，进行实时的弥补。

3.前馈-反馈结构

• 模型中的预测作用考虑到了来自未来参考输入和可测干扰的前馈补偿
• 作用在系统中的控制量考虑到了当前测量的状态值的反馈补偿

从无约束模型预测控制出发理解推导过程

这里我们使用控制里面常见的状态空间模型，并且在不考虑约束的情况下进行一个详细的公式推导。主要分为三个步骤：使用模型预测系统未来动态，求解优化问题，解的第一个元素作用到系统。

1.使用状态空间模型

这里直接给出离散情况下系统的状态空间模型：
$$\begin{array}{rl}x(k+1)& =Ax(k)+B(k)u(k)\\ y_c(k)& =Cx(k)\end{array}$$
这里为了引入积分以减少或消除静态误差，将模型改写成增量形式，设：
$$\begin{gathered} \Delta x(k)=x(k)-x(k-1) \\ \Delta u(k)=u(k)-u(k-1) \end{gathered}$$
这样状态空间模型就可以改写成：
$$\begin{array}{rl}\Delta x(k+1)& =A\Delta x(k)+B\Delta u(k)\\ y(k)& =C\Delta x(k)+y(k-1)\end{array}$$

2.模型的预测功能

这是MPC求解的第一步：使用模型预测。
我们可以观察状态方程，其实它表达的是：给定一个当前状态和输入，就能计算出下一时刻的状态，那么这里的预测就是计算出的下一时刻状态。而有了下一时刻的状态和输入，就可以计算下下时刻的状态，以此类推，想预测之后多久的都可以。当前状态是测量出来的，而各个时间的输入是需要优化的变量，根据这个思路进行公式推导(这就是使用模型进行预测的体现)。这里设置预测时域为$p$，控制时域为$m$。并且$m\le p$。
首先根据当前时刻$k$拿到的数据和优化变量控制对下一时刻进行预测：
$$\Delta x(k+1|k)=A\Delta x(k)+B\Delta u(k)$$
现在进行一步预测之后，就拥有了当前状态和下一时刻状态，同理，用$k+1$时刻状态也可以对$k+2$进行预测。
$$\Delta x(k+2|k)=A\Delta x(k+1|k)+B\Delta u(k+1)$$
进一步，使用$k$时刻测量的量表示$k+2$时刻的预测值为：
$$\begin{array}{rl}\Delta x(k+2|k)& =A\Delta x(k+1|k)+B\Delta u(k+1)\\ & =A(A\Delta x(k)+B\Delta u(k))+B\Delta u(k+1)\\ & =A^2\Delta x(k)+AB\Delta u(k)+B\Delta u(k+1)\end{array}$$
以此类推，如果我们希望预测控制步长结束之前最后的第$m$步，那么有：
$$\Delta x(k+m|k)=A^m\Delta x(k)+A^{m-1}B\Delta u(k)+A^{m-2}B\Delta u(k+1)+...+B\Delta u(k+m-1)$$
当控制步长结束之后，我们假设控制量不变，也就是$\Delta u=0$,那么有：
$$\Delta x(k+p|k)=A^p\Delta x(k)+A^{p-1}B\Delta u(k)+A^{p-2}B\Delta u(k+1)+...+A^{p-m}B\Delta u(k+m-1)$$
同样的思路，对于输出方程在控制步长没有结束之前有：
$$y(k+m|k)=\sum _{i=1}^{m}C{A}^i\Delta x(k)+\sum _{i=1}^{m}C{A}^{i-1}B\Delta u(k)+\sum _{i=1}^{m-1}C{A}^{i-1}B\Delta u(k+1)+...+CB\Delta u(k+m-1)+y(k)$$
在控制步长结束之后：
$$y(k+p|k)=\sum _{i=1}^{p}C{A}^i\Delta x(k)+\sum _{i=1}^{p}C{A}^{i-1}B\Delta u(k)+\sum _{i=1}^{p-1}C{A}^{i-1}B\Delta u(k+1)+...+\sum _{i=1}^{p-m+1}C{A}^{i-1}B\Delta u(k+m-1)+y(k)$$
以上就完成了预测。虽然写成这样就可以，但是为了方便计算和之后的推导，我们把它写成矩阵的形式（其实就是一个简单的列排列）。
为了方便理解这里先用$p=2$举个例子：
$$\begin{array}{rl}y(k+1|k)& =CA\Delta x(k)+CB\Delta u(k)+y(k)\\ y(k+2|k)& =(C{A}^2+CA)\Delta x(k)+(CAB+CB)\Delta u(k)+CB\Delta u(k+1)+y(k)\end{array}$$
这两个式子并列之后，可以整理成矩阵的形式：
$$\left(\begin{array}{c}y(k+1|k)\\ y(k+2|k)\end{array}\right)=\left(\begin{array}{c}CA\\ (C{A}^2+CA)\end{array}\right)\Delta x(k)+\left(\begin{array}{cc}CB& 0\\ CAB& CB\end{array}\right)\left(\begin{array}{c}\Delta u(k)\\ \Delta u(k+1)\end{array}\right)+\left(\begin{array}{c}{I}_{n_y\times n_y}\\ I_{n_y\times n_y}\end{array}\right)y(k)$$
这里需要说明一些，$\Delta x,\Delta u,y$一般都是一个向量，所以后面的$y$前面的系数是一个个单位矩阵。在明白矩阵构成方式之后，将其扩展到$k=p$时刻,可以写成如下形式：
$$Y(k+1|k)=S_x\Delta x(k)+S_u\Delta U(k)+Iy(k)$$
其中：
$$Y(k+1|k)={\left(\begin{array}{c}y(k+1|k)\\ y(k+2|k)\\ ⋮\\ y(k+p|k)\end{array}\right)}_{p\times 1}$$
$$\Delta U(k)={\left(\begin{array}{c}\Delta u(k)\\ \Delta u(k+1)\\ ⋮\\ \Delta u(k+m-1)\end{array}\right)}_{m\times 1}$$
$$S_x={\left(\begin{array}{c}CA\\ {\sum }_{i=1}^{2}C{A}^i\\ ⋮\\ {\sum }_{i=1}^{p}C{A}^i\end{array}\right)}_{p\times 1}$$
$$I={\left(\begin{array}{c}{I}_{n_y\times n_y}\\ I_{n_y\times n_y}\\ ⋮\\ I_{n_y\times n_y}\end{array}\right)}_{p\times 1}$$
$$S_u=\left(\begin{array}{cccccc}CB& 0& 0& \cdots & 0& \\ {\sum }_{i=1}^{2}C{A}^{i-1}B& CB& 0& \cdots & 0& \\ ⋮& ⋮& ⋮& \ddots & ⋮& \\ {\sum }_{i=1}^{m}C{A}^{i-1}B& {\sum }_{i=1}^{m-1}C{A}^{i-1}B& \cdots & \cdots & CB& \\ {\sum }_{i=1}^{p}C{A}^{i-1}B& {\sum }_{i=1}^{p-1}C{A}^{i-1}B& \cdots & \cdots & {\sum }_{i=1}^{p-m+1}C{A}^{i-1}B& \end{array}\right)$$
这里可以看出，$S_u$是一个下三角，这直接反映了控制量在时间上的因果关系，即之后的输入对之前的输入没有影响。

3.优化问题的建立

一般来说，目标函数的选取没有标准形式，这是一种艺术而不是数学。考虑一种最简单的形式，我们希望被控输出接近参考输入：
$$J=\sum _{i=1}^p{Q}_i(y(k+i|k)-r(k+i){)}^2$$
这里$Q_i$是权重，表示的是在时刻$i$对各个分量误差的加权。通常来说，还希望控制的变化不要太大：
$$J=\sum _{i=1}^p{Q}_i(y(k+i|k)-r(k+i){)}^2+\sum _{i=1}^m{H}_i(\Delta u(k+i-1){)}^2$$
设参考序列：
$$R(k+1)={\left(\begin{array}{c}r(k+1)\\ r(k+2)\\ ⋮\\ r(k+p)\end{array}\right)}_{p\times 1}$$
这样我们就可以将目标函数写成矩阵的形式：
$$J=Q(Y(k+1|k)-R(k+1){)}^2+H(\Delta U(k){)}^2$$
这里可以看到目标函数都是平方项，如果我们定义：
$$\rho =\left(\begin{array}{c}\sqrt{Q}(Y(k+1|k)-R(k+1))\\ \sqrt{H}\Delta U(k)\end{array}\right)$$
那么目标函数就可以简化成：
$$J={\rho }^T\rho$$
这样我们就可以将开环优化问题描述成如下形式：
$$\begin{gathered} \underset{\Delta U(k)}{\min }J(x(k),\Delta U(k),m,p) \\ \begin{array}{r}s.t.Y(k+1|k)=S_x\Delta x(k)+S_u\Delta U(k)+Iy(k)\end{array} \end{gathered}$$

4.优化问题求解

这是MPC求解的第二和第三步：求解优化问题，将第一个解作用到系统。
这里举一个例子，首先看一个很简单的优化问题：
$$\begin{array}{rl}& min{x}_1^2+x_2^2\\ & s.t.x_1+x_2=1\end{array}$$
如果是你你怎么求解，当然可以拉格朗日，但是我选择让$x_1=1-x_2$，然后带到目标函数中，进行求解。这个思路有了之后，其实MPC的推导，也是一个道理，只不过是将模型放进了目标函数里面。
我们现在将$Y$带到$\rho$的第一项中：
$$\begin{array}{rl}\rho & =\left(\begin{array}{c}\sqrt{Q}(S_x\Delta x(k)+S_u\Delta U(k)+Iy(k)-R(k+1))\\ \sqrt{H}\Delta U(k)\end{array}\right)\\ & =\left(\begin{array}{c}\sqrt{Q}{S}_u\\ \sqrt{H}\end{array}\right)\Delta U(k)-\left(\begin{array}{c}\sqrt{Q}(R(k+1)-S_x\Delta x(k)-Iy(k))\\ 0\end{array}\right)\\ & =\left(\begin{array}{c}\sqrt{Q}{S}_u\\ \sqrt{H}\end{array}\right)\Delta U(k)-\left(\begin{array}{c}\sqrt{Q}E(k+1|k)\\ 0\end{array}\right)\\ & =Az-b\end{array}$$
对${\rho }^T\rho =(Az-b{)}^T(Az-b)$求导等于零：
$$\frac{\mathrm{d}{\rho }^T\rho }{\mathrm{d}z}=2A^T(Az-b)=0$$
得到最优解为：
$$z^{\ast }=(A^{T}A{)}^{-1}{A}^{T}b$$
这样我们就可以求解出最优解了，这里需要强调一下，对于
$$E(k+1|k)=R(k+1)-S_x\Delta x(k)-Iy(k)$$
这里面都是计算或者观测出来的量。如果将结果$z^{\ast }$写成：
$$z^{\ast }=\Delta U^{\ast }=K(R(k+1)-S_x\Delta x(k)-Iy(k))$$
可以看出，这里面的前馈反映在未来的参考输入$R(k+1)$,反馈则反映在当前时刻获取的$\Delta x(k)和y(k)$上。

加入约束之后的计算方法

对于约束的解释

可能有人有疑问，刚才上面的模型构建中其实是有约束的(模型约束)。这里说的约束一般指的是对自变量范围的约束，而不是模型的约束。那又会有人有疑问，目标函数里面的加权项考虑了对控制动作变化不能太大的要求，这其实是一个软约束。MPC中的约束，是指在优化问题中包含时域硬约束：
$$\begin{array}{rl}{u}_{min}\le & u(k+i)\le u_{max},i=0,1,...,m-1\\ \Delta u_{min}\le & \Delta u(k+i)\le \Delta u_{max},i=0,1,...,m-1\end{array}$$

目标函数抛弃常数项，转化成$z^{T}Hz-g^{T}z$

这一步的思路很简单，这里先给个简单的例子,比如你想求$x^2+3x+1$的最小化的那个点$x^{\ast }$，完全可以不管后面的$1$，求$x^2+3x$的$x^{\ast }$就可以了。初中就知道$-b/2a$对吧。这里一个道理。对于刚才的目标函数省去最后一项，直接可以推导出：
$$\begin{array}{r}{J}_1=\Delta U(k{)}^{T}H\Delta U(k)-G(k+1|k{)}^T\Delta U(k)\end{array}$$
其中：
$$\begin{gathered} H=S_u^{T}Q{S}_u+H \\ G(k+1|k)=2S_u^{T}QE(k+1|k) \end{gathered}$$

控制增量约束转化成$Cz\le b$

对于控制的增量：
$$\Delta u_{min}\le \Delta u(k+i)\le \Delta u_{max},i=0,1,...,m-1$$
可以将其竖着排列，就会有：
$$\left(\begin{array}{c}\Delta u_{min}\\ \Delta u_{min}\\ ⋮\\ \Delta u_{min}\end{array}\right)\le \Delta U(k)\le \left(\begin{array}{c}\Delta u_{max}\\ \Delta u_{max}\\ ⋮\\ \Delta u_{max}\end{array}\right)$$
将左边的小于乘一个负号，这样可以改写成：
$$\begin{gathered} \Delta U(k)\le \left(\begin{array}{c}\Delta u_{max}\\ \Delta u_{max}\\ ⋮\\ \Delta u_{max}\end{array}\right) \\ -\Delta U(k)\le \left(\begin{array}{c}-\Delta u_{min}\\ -\Delta u_{min}\\ ⋮\\ -\Delta u_{min}\end{array}\right) \end{gathered}$$
这样合并成：
$$\left(\begin{array}{c}T\\ -T\end{array}\right)\Delta U(k)\le \left(\begin{array}{c}\Delta u_{max}\\ \Delta u_{max}\\ ⋮\\ \Delta u_{max}\\ -\Delta u_{min}\\ -\Delta u_{min}\\ ⋮\\ -\Delta u_{min}\end{array}\right)$$
其中：
$$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)$$

控制量约束转化成$Cz\le b$

由于有如下关系：
$$\Delta u(k+i)=u(k+i)-u(k+i-1)$$
因此当$i=0$有$u(k)=u(k-1)+\Delta u(k)$,进而有：
$$u_{min}\le u(k)\le u_{max}⟹\{\begin{array}{l}\Delta u(k)\le u_{max}-u(k-1)\\ -\Delta u(k)\le -(u_{min}-u(k-1))\end{array}$$
对$i=1$有$u(k+1)=u(k-1)+\Delta u(k)+\Delta u(k+1)$，因此有
$$u_{min}\le u(k+1)\le u_{max}⟹\{\begin{array}{l}\Delta u(k)+\Delta u(k+1)\le u_{max}-u(k-1)\\ -(\Delta u(k)+\Delta u(k+1))\le -(u_{min}-u(k-1))\end{array}$$
以此类推，对于任何$i=0,1,2...,m-1$：
$$u_{min}\le u(k+i)\le u_{max}⟹\{\begin{array}{l}{\sum }_{j=0}^i\Delta u(k+j)\le u_{max}-u(k-1)\\ -{\sum }_{j=0}^i\Delta u(k+j)\le -(u_{min}-u(k-1))\end{array}$$
使用同样的方法，可以将其写成如下形式：
$$\left(\begin{array}{c}L\\ -L\end{array}\right)\Delta U(k)\le \left(\begin{array}{c}{u}_{max}-u(k-1)\\ u_{max}-u(k-1)\\ ⋮\\ u_{max}-u(k-1)\\ -(u_{min}-u(k-1))\\ -(u_{min}-u(k-1))\\ ⋮\\ -(u_{min}-u(k-1))\end{array}\right)$$
其中：
$$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)$$

转化为QP问题

其实个人认为MPC最后转化为什么问题求解绝大多数情况取决于目标函数的设计。我们可以看到目标函数设计成了二次型的，动力学方程和时域约束条件是线性的，所以问题就变成了一个二次规划(QP)问题:
$$\begin{gathered} \underset{z}{\min }{z}^{T}Hz-g^{T}z \\ s.t.Cz\le b \end{gathered}$$
一般对于这种问题最常用的两种方法是：积极集法(active-set method)和内点法(interior-point method)。这样我们就把MPC问题变成了QP问题进行求解即可。之后将求出的第一个变量作用上去，并循环这个过程。对应于现实问题比如：

• 速度跟踪，可以参考：【Carsim Simulink自动驾驶仿真】基于MPC的速度控制
• 轨迹跟踪，可以参考：【Carsim Simulink自动驾驶仿真】基于MPC的轨迹跟踪控制