MPC模型预测控制理论分析详解

MPC的基本原理与特点
- 基本原理
- 特点
- - 1.基于模型
  - 2.滚动优化
  - 3.前馈-反馈结构
从无约束模型预测控制出发理解推导过程
- 1.使用状态空间模型
- 2.模型的预测功能
- 3.优化问题的建立
- 4.优化问题求解
加入约束之后的计算方法
- 对于约束的解释
- 目标函数抛弃常数项，转化成zTHz−gTzz^THz-g^TzzTHz−gTz
- 控制增量约束转化成Cz≤bCz≤bCz≤b
- 控制量约束转化成Cz≤bCz≤bCz≤b
- 转化为QP问题

MPC的基本原理与特点

基本原理

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

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

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

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

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

特点

1.基于模型

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

2.滚动优化

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

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

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

3.前馈-反馈结构

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

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

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

1.使用状态空间模型

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

2.模型的预测功能

这是MPC求解的第一步：使用模型预测。
我们可以观察状态方程，其实它表达的是：给定一个当前状态和输入，就能计算出下一时刻的状态，那么这里的预测就是计算出的下一时刻状态。而有了下一时刻的状态和输入，就可以计算下下时刻的状态，以此类推，想预测之后多久的都可以。当前状态是测量出来的，而各个时间的输入是需要优化的变量，根据这个思路进行公式推导(这就是使用模型进行预测的体现)。这里设置预测时域为ppp，控制时域为mmm。并且m≤pm≤pm≤p。
首先根据当前时刻kkk拿到的数据和优化变量控制对下一时刻进行预测：
Δx(k+1∣k)=AΔx(k)+BΔu(k)Δx(k+1|k) = AΔx(k)+BΔu(k) Δx(k+1∣k)=AΔx(k)+BΔu(k)
现在进行一步预测之后，就拥有了当前状态和下一时刻状态，同理，用k+1k+1k+1时刻状态也可以对k+2k+2k+2进行预测。
Δx(k+2∣k)=AΔx(k+1∣k)+BΔu(k+1)Δx(k+2|k) = AΔx(k+1|k)+BΔu(k+1) Δx(k+2∣k)=AΔx(k+1∣k)+BΔu(k+1)
进一步，使用kkk时刻测量的量表示k+2k+2k+2时刻的预测值为：
Δx(k+2∣k)=AΔx(k+1∣k)+BΔu(k+1)=A(AΔx(k)+BΔu(k))+BΔu(k+1)=A2Δx(k)+ABΔu(k)+BΔu(k+1)\begin{aligned} Δx(k+2|k) &= AΔx(k+1|k)+BΔu(k+1)\\[2mm] &= A( AΔx(k)+BΔu(k))+BΔu(k+1)\\[2mm] &=A^2Δx(k)+ABΔu(k)+BΔu(k+1) \end{aligned} Δx(k+2∣k)=AΔx(k+1∣k)+BΔu(k+1)=A(AΔx(k)+BΔu(k))+BΔu(k+1)=A2Δx(k)+ABΔu(k)+BΔu(k+1)
以此类推，如果我们希望预测控制步长结束之前最后的第mmm步，那么有：
Δx(k+m∣k)=AmΔx(k)+Am−1BΔu(k)+Am−2BΔu(k+1)+...+BΔu(k+m−1)Δx(k+m|k) = A^mΔx(k)+A^{m-1}BΔu(k)+A^{m-2}BΔu(k+1)+...+BΔu(k+m-1) Δx(k+m∣k)=AmΔx(k)+Am−1BΔu(k)+Am−2BΔu(k+1)+...+BΔu(k+m−1)
当控制步长结束之后，我们假设控制量不变，也就是Δu=0Δu=0Δu=0,那么有：
Δx(k+p∣k)=ApΔx(k)+Ap−1BΔu(k)+Ap−2BΔu(k+1)+...+Ap−mBΔu(k+m−1)Δx(k+p|k) = A^pΔx(k)+A^{p-1}BΔu(k)+A^{p-2}BΔu(k+1)+...+A^{p-m}BΔu(k+m-1) Δx(k+p∣k)=ApΔx(k)+Ap−1BΔu(k)+Ap−2BΔu(k+1)+...+Ap−mBΔu(k+m−1)
同样的思路，对于输出方程在控制步长没有结束之前有：
y(k+m∣k)=∑i=1mCAiΔx(k)+∑i=1mCAi−1BΔu(k)+∑i=1m−1CAi−1BΔu(k+1)+...+CBΔu(k+m−1)+y(k)y(k+m|k)=\sum_{i=1}^mCA^iΔx(k)+\sum_{i=1}^mCA^{i-1}BΔu(k)+\sum_{i=1}^{m-1}CA^{i-1}BΔu(k+1)+...+CBΔu(k+m-1)+y(k)\\ y(k+m∣k)=i=1∑mCAiΔx(k)+i=1∑mCAi−1BΔu(k)+i=1∑m−1CAi−1BΔu(k+1)+...+CBΔu(k+m−1)+y(k)
在控制步长结束之后：
y(k+p∣k)=∑i=1pCAiΔx(k)+∑i=1pCAi−1BΔu(k)+∑i=1p−1CAi−1BΔu(k+1)+...+∑i=1p−m+1CAi−1BΔu(k+m−1)+y(k)y(k+p|k)=\sum_{i=1}^pCA^iΔx(k)+\sum_{i=1}^pCA^{i-1}BΔu(k)+\sum_{i=1}^{p-1}CA^{i-1}BΔu(k+1)+...+\sum_{i=1}^{p-m+1}CA^{i-1}BΔu(k+m-1)+y(k) y(k+p∣k)=i=1∑pCAiΔx(k)+i=1∑pCAi−1BΔu(k)+i=1∑p−1CAi−1BΔu(k+1)+...+i=1∑p−m+1CAi−1BΔu(k+m−1)+y(k)
以上就完成了预测。虽然写成这样就可以，但是为了方便计算和之后的推导，我们把它写成矩阵的形式（其实就是一个简单的列排列）。
为了方便理解这里先用p=2p=2p=2举个例子：
y(k+1∣k)=CAΔx(k)+CBΔu(k)+y(k)y(k+2∣k)=(CA2+CA)Δx(k)+(CAB+CB)Δu(k)+CBΔu(k+1)+y(k)\begin{aligned} y(k+1|k) &= CAΔx(k)+CBΔu(k)+y(k)\\ y(k+2|k) &= (CA^2+CA)Δx(k)+(CAB+CB)Δu(k)+CBΔu(k+1)+y(k) \end{aligned} y(k+1∣k)y(k+2∣k)=CAΔx(k)+CBΔu(k)+y(k)=(CA2+CA)Δx(k)+(CAB+CB)Δu(k)+CBΔu(k+1)+y(k)
这两个式子并列之后，可以整理成矩阵的形式：
(y(k+1∣k)y(k+2∣k))=(CA(CA2+CA))Δx(k)+(CB0CABCB)(Δu(k)Δu(k+1))+(Iny×nyIny×ny)y(k)\left( \begin{array}{c} y(k+1|k)\\ y(k+2|k) \end{array} \right)= \left( \begin{array}{c} CA\\ (CA^2+CA) \end{array} \right)Δx(k)+ \left( \begin{array}{ccc} CB &0\\ CAB &CB \end{array} \right) \left( \begin{array}{c} Δu(k)\\ Δu(k+1) \end{array} \right)+ \left( \begin{array}{c} I_{n_y×n_y}\\ I_{n_y×n_y} \end{array} \right)y(k) (y(k+1∣k)y(k+2∣k))=(CA(CA2+CA))Δx(k)+(CBCAB0CB)(Δu(k)Δu(k+1))+(Iny×nyIny×ny)y(k)
这里需要说明一些，Δx,Δu,yΔx,Δu,yΔx,Δu,y一般都是一个向量，所以后面的yyy前面的系数是一个个单位矩阵。在明白矩阵构成方式之后，将其扩展到k=pk=pk=p时刻,可以写成如下形式：
Y(k+1∣k)=SxΔx(k)+SuΔU(k)+Iy(k)Y(k+1|k)=S_xΔx(k)+S_uΔU(k)+Iy(k) Y(k+1∣k)=SxΔx(k)+SuΔU(k)+Iy(k)
其中：
Y(k+1∣k)=(y(k+1∣k)y(k+2∣k)⋮y(k+p∣k))p×1Y(k+1|k)= \left( \begin{array}{c} y(k+1|k)\\ y(k+2|k)\\ \vdots\\ y(k+p|k) \end{array} \right)_{p×1} Y(k+1∣k)=⎝⎜⎜⎜⎛y(k+1∣k)y(k+2∣k)⋮y(k+p∣k)⎠⎟⎟⎟⎞p×1
ΔU(k)=(Δu(k)Δu(k+1)⋮Δu(k+m−1))m×1ΔU(k)= \left( \begin{array}{c} Δu(k)\\ Δu(k+1)\\ \vdots\\ Δu(k+m-1) \end{array} \right)_{m×1} ΔU(k)=⎝⎜⎜⎜⎛Δu(k)Δu(k+1)⋮Δu(k+m−1)⎠⎟⎟⎟⎞m×1
Sx=(CA∑i=12CAi⋮∑i=1pCAi)p×1S_x= \left( \begin{array}{c} CA\\ \sum_{i=1}^2CA^i\\ \vdots\\ \sum_{i=1}^pCA^i \end{array} \right)_{p×1} Sx=⎝⎜⎜⎜⎛CA∑i=12CAi⋮∑i=1pCAi⎠⎟⎟⎟⎞p×1
I=(Iny×nyIny×ny⋮Iny×ny)p×1I= \left( \begin{array}{c} I_{n_y×n_y}\\ I_{n_y×n_y}\\ \vdots\\ I_{n_y×n_y} \end{array} \right)_{p×1} I=⎝⎜⎜⎜⎛Iny×nyIny×ny⋮Iny×ny⎠⎟⎟⎟⎞p×1
Su=(CB00⋯0∑i=12CAi−1BCB0⋯0⋮⋮⋮⋱⋮∑i=1mCAi−1B∑i=1m−1CAi−1B⋯⋯CB∑i=1pCAi−1B∑i=1p−1CAi−1B⋯⋯∑i=1p−m+1CAi−1B)S_u= \left( \begin{array}{ccc} CB &0&0&\cdots&0\\[2mm] \sum_{i=1}^2CA^{i-1}B &CB&0&\cdots&0\\[2mm] \vdots&\vdots&\vdots&\ddots&\vdots&\\[2mm] \sum_{i=1}^mCA^{i-1}B&\sum_{i=1}^{m-1}CA^{i-1}B&\cdots&\cdots&CB\\[2mm] \sum_{i=1}^pCA^{i-1}B&\sum_{i=1}^{p-1}CA^{i-1}B&\cdots&\cdots& \sum_{i=1}^{p-m+1}CA^{i-1}B \end{array} \right) Su=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛CB∑i=12CAi−1B⋮∑i=1mCAi−1B∑i=1pCAi−1B0CB⋮∑i=1m−1CAi−1B∑i=1p−1CAi−1B00⋮⋯⋯⋯⋯⋱⋯⋯00⋮CB∑i=1p−m+1CAi−1B⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
这里可以看出，SuS_uSu是一个下三角，这直接反映了控制量在时间上的因果关系，即之后的输入对之前的输入没有影响。

3.优化问题的建立

一般来说，目标函数的选取没有标准形式，这是一种艺术而不是数学。考虑一种最简单的形式，我们希望被控输出接近参考输入：
J=∑i=1pQi(y(k+i∣k)−r(k+i))2J=\sum_{i=1}^pQ_i(y(k+i|k)-r(k+i))^2 J=i=1∑pQi(y(k+i∣k)−r(k+i))2
这里QiQ_iQi是权重，表示的是在时刻iii对各个分量误差的加权。通常来说，还希望控制的变化不要太大：
J=∑i=1pQi(y(k+i∣k)−r(k+i))2+∑i=1mHi(Δu(k+i−1))2J=\sum_{i=1}^pQ_i(y(k+i|k)-r(k+i))^2+\sum_{i=1}^mH_i(Δu(k+i-1))^2 J=i=1∑pQi(y(k+i∣k)−r(k+i))2+i=1∑mHi(Δu(k+i−1))2
设参考序列：
R(k+1)=(r(k+1)r(k+2)⋮r(k+p))p×1R(k+1)=\left( \begin{array}{c} r(k+1)\\ r(k+2)\\ \vdots\\ r(k+p) \end{array} \right)_{p×1} R(k+1)=⎝⎜⎜⎜⎛r(k+1)r(k+2)⋮r(k+p)⎠⎟⎟⎟⎞p×1
这样我们就可以将目标函数写成矩阵的形式：
J=Q(Y(k+1∣k)−R(k+1))2+H(ΔU(k))2J=Q(Y(k+1|k)-R(k+1))^2+H(ΔU(k))^2 J=Q(Y(k+1∣k)−R(k+1))2+H(ΔU(k))2
这里可以看到目标函数都是平方项，如果我们定义：
ρ=(Q(Y(k+1∣k)−R(k+1))HΔU(k))ρ=\left( \begin{array}{c} \sqrt{Q}(Y(k+1|k)-R(k+1))\\ \sqrt{H}ΔU(k)\\ \end{array} \right) ρ=(Q(Y(k+1∣k)−R(k+1))HΔU(k))
那么目标函数就可以简化成：
J=ρTρJ=ρ^Tρ J=ρTρ
这样我们就可以将开环优化问题描述成如下形式：
min⁡ΔU(k)J(x(k),ΔU(k),m,p)s.t.Y(k+1∣k)=SxΔx(k)+SuΔU(k)+Iy(k)\min_{ΔU(k)}J(x(k),ΔU(k),m,p)\\[2mm] \begin{aligned} s.t.　Y(k+1|k)=S_xΔx(k)+S_uΔU(k)+Iy(k) \end{aligned} ΔU(k)minJ(x(k),ΔU(k),m,p)s.t.　Y(k+1∣k)=SxΔx(k)+SuΔU(k)+Iy(k)

4.优化问题求解

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

加入约束之后的计算方法

对于约束的解释

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

目标函数抛弃常数项，转化成zTHz−gTzz^THz-g^TzzTHz−gTz

这一步的思路很简单，这里先给个简单的例子,比如你想求x2+3x+1x^2+3x+1x2+3x+1的最小化的那个点x∗x^*x∗，完全可以不管后面的111，求x2+3xx^2+3xx2+3x的x∗x^*x∗就可以了。初中就知道−b/2a-b/2a−b/2a对吧。这里一个道理。对于刚才的目标函数省去最后一项，直接可以推导出：
J1=ΔU(k)THΔU(k)−G(k+1∣k)TΔU(k)\begin{aligned} J_1 = ΔU(k)^THΔU(k)-G(k+1|k)^TΔU(k) \end{aligned} J1=ΔU(k)THΔU(k)−G(k+1∣k)TΔU(k)
其中：
H=SuTQSu+HG(k+1∣k)=2SuTQE(k+1∣k)H = S_u^TQS_u+H\\[2mm] G(k+1|k)=2S_u^TQE(k+1|k) H=SuTQSu+HG(k+1∣k)=2SuTQE(k+1∣k)

控制增量约束转化成Cz≤bCz≤bCz≤b

对于控制的增量：
Δumin≤Δu(k+i)≤Δumax,i=0,1,...,m−1Δu_{min}≤Δu(k+i)≤Δu_{max},　i=0,1,...,m-1 Δumin≤Δu(k+i)≤Δumax,　i=0,1,...,m−1
可以将其竖着排列，就会有：
(ΔuminΔumin⋮Δumin)≤ΔU(k)≤(ΔumaxΔumax⋮Δumax)\left( \begin{array}{c} Δu_{min}\\ Δu_{min}\\ \vdots\\ Δu_{min}\\ \end{array} \right)≤ΔU(k)≤ \left( \begin{array}{c} Δu_{max}\\ Δu_{max}\\ \vdots\\ Δu_{max}\\ \end{array} \right) ⎝⎜⎜⎜⎛ΔuminΔumin⋮Δumin⎠⎟⎟⎟⎞≤ΔU(k)≤⎝⎜⎜⎜⎛ΔumaxΔumax⋮Δumax⎠⎟⎟⎟⎞
将左边的小于乘一个负号，这样可以改写成：
ΔU(k)≤(ΔumaxΔumax⋮Δumax)−ΔU(k)≤(−Δumin−Δumin⋮−Δumin)ΔU(k)≤ \left( \begin{array}{c} Δu_{max}\\ Δu_{max}\\ \vdots\\ Δu_{max}\\ \end{array} \right)\\ -ΔU(k)≤ \left( \begin{array}{c} -Δu_{min}\\ -Δu_{min}\\ \vdots\\ -Δu_{min}\\ \end{array} \right) ΔU(k)≤⎝⎜⎜⎜⎛ΔumaxΔumax⋮Δumax⎠⎟⎟⎟⎞−ΔU(k)≤⎝⎜⎜⎜⎛−Δumin−Δumin⋮−Δumin⎠⎟⎟⎟⎞
这样合并成：
(T−T)ΔU(k)≤(ΔumaxΔumax⋮Δumax−Δumin−Δumin⋮−Δumin)\left( \begin{array}{c} T\\ -T\\ \end{array} \right)ΔU(k)≤ \left( \begin{array}{c} Δu_{max}\\ Δu_{max}\\ \vdots\\ Δu_{max}\\ -Δu_{min}\\ -Δu_{min}\\ \vdots\\ -Δu_{min}\\ \end{array} \right) (T−T)ΔU(k)≤⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛ΔumaxΔumax⋮Δumax−Δumin−Δumin⋮−Δumin⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
其中：
T=(Inu×nu0⋯00Inu×nu⋯0⋮⋮⋮00⋯Inu×nu)T=\left( \begin{array}{ccc} I_{n_u×n_u}&0& \cdots&0\\[2mm] 0&I_{n_u×n_u}&\cdots&0\\[2mm] \vdots&\vdots&&\vdots\\[2mm] 0&0&\cdots&I_{n_u×n_u}\\[2mm] \end{array} \right) T=⎝⎜⎜⎜⎜⎜⎜⎜⎛Inu×nu0⋮00Inu×nu⋮0⋯⋯⋯00⋮Inu×nu⎠⎟⎟⎟⎟⎟⎟⎟⎞

控制量约束转化成Cz≤bCz≤bCz≤b

由于有如下关系：
Δu(k+i)=u(k+i)−u(k+i−1)Δu(k+i)=u(k+i)-u(k+i-1) Δu(k+i)=u(k+i)−u(k+i−1)
因此当i=0i=0i=0有u(k)=u(k−1)+Δu(k)u(k)=u(k-1)+Δu(k)u(k)=u(k−1)+Δu(k),进而有：
umin≤u(k)≤umax⟹{Δu(k)≤umax−u(k−1)−Δu(k)≤−(umin−u(k−1))u_{min}≤u(k)≤u_{max} \Longrightarrow\begin{cases} Δu(k)≤u_{max}-u(k-1)\\ -Δu(k)≤-(u_{min}-u(k-1)) \end{cases} umin≤u(k)≤umax⟹{Δu(k)≤umax−u(k−1)−Δu(k)≤−(umin−u(k−1))
对i=1i=1i=1有u(k+1)=u(k−1)+Δu(k)+Δu(k+1)u(k+1)=u(k-1)+Δu(k)+Δu(k+1)u(k+1)=u(k−1)+Δu(k)+Δu(k+1)，因此有
umin≤u(k+1)≤umax⟹{Δu(k)+Δu(k+1)≤umax−u(k−1)−(Δu(k)+Δu(k+1))≤−(umin−u(k−1))u_{min}≤u(k+1)≤u_{max} \Longrightarrow\begin{cases} Δu(k)+Δu(k+1)≤u_{max}-u(k-1)\\ -(Δu(k)+Δu(k+1))≤-(u_{min}-u(k-1)) \end{cases}umin≤u(k+1)≤umax⟹{Δu(k)+Δu(k+1)≤umax−u(k−1)−(Δu(k)+Δu(k+1))≤−(umin−u(k−1))
以此类推，对于任何i=0,1,2...,m−1i=0,1,2...,m-1i=0,1,2...,m−1：
umin≤u(k+i)≤umax⟹{∑j=0iΔu(k+j)≤umax−u(k−1)−∑j=0iΔu(k+j)≤−(umin−u(k−1))u_{min}≤u(k+i)≤u_{max} \Longrightarrow\begin{cases} \sum_{j=0}^iΔu(k+j)≤u_{max}-u(k-1)\\[3mm] -\sum_{j=0}^iΔu(k+j)≤-(u_{min}-u(k-1)) \end{cases}umin≤u(k+i)≤umax⟹⎩⎨⎧∑j=0iΔu(k+j)≤umax−u(k−1)−∑j=0iΔu(k+j)≤−(umin−u(k−1))
使用同样的方法，可以将其写成如下形式：
(L−L)ΔU(k)≤(umax−u(k−1)umax−u(k−1)⋮umax−u(k−1)−(umin−u(k−1))−(umin−u(k−1))⋮−(umin−u(k−1)))\left( \begin{array}{c} L\\ -L\\ \end{array} \right)ΔU(k)≤ \left( \begin{array}{c} u_{max}-u(k-1)\\ u_{max}-u(k-1)\\ \vdots\\ u_{max}-u(k-1)\\ -(u_{min}-u(k-1))\\ -(u_{min}-u(k-1))\\ \vdots\\ -(u_{min}-u(k-1))\\ \end{array} \right) (L−L)ΔU(k)≤⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛umax−u(k−1)umax−u(k−1)⋮umax−u(k−1)−(umin−u(k−1))−(umin−u(k−1))⋮−(umin−u(k−1))⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
其中：
L=(Inu×nu0⋯0Inu×nuInu×nu⋯0⋮⋮⋮Inu×nuInu×nu⋯Inu×nu)L=\left( \begin{array}{ccc} I_{n_u×n_u}&0& \cdots&0\\[2mm] I_{n_u×n_u}&I_{n_u×n_u}&\cdots&0\\[2mm] \vdots&\vdots&&\vdots\\[2mm] I_{n_u×n_u}&I_{n_u×n_u}&\cdots&I_{n_u×n_u}\\[2mm] \end{array} \right) L=⎝⎜⎜⎜⎜⎜⎜⎜⎛Inu×nuInu×nu⋮Inu×nu0Inu×nu⋮Inu×nu⋯⋯⋯00⋮Inu×nu⎠⎟⎟⎟⎟⎟⎟⎟⎞

转化为QP问题

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

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

如何理解MPC模型预测控制理论相关推荐

mpc模型预测控制原理详解
mpc模型预测控制原理详解前言 mpc算法步骤 mpc算法推导前言本文是对mpc模型预测控制学习的记录,主要参照了DR_CAN老师的视频进行学习.视频专栏链接:DR_CAN老师mpc视频专栏.在 ...
Matlab/simulink 风储调频，风电调频，模糊控制，mpc模型预测，虚拟惯性控制，下垂控制。
风电采用虚拟惯性控制,储能采用下垂控制. 风电利用模糊控制改变系数,根据风速和频率变化系数.做到了自适应控制. MPC通过状态空间表达式通过fk预测下一时刻fk+1频率,仿真结果表明预测频率和实际频率 ...
Matlab/simulink MPC模型预测风电调频，风储调频，一次调频，频域模型。
Mpc预测风电出力,进而使得风电参与一次调频,风电内部结构详细.风储联合调频.
MPC模型预测控制及在Matlab中实现函数定义
基于b站DR_CAN老师的MPC控制视频[MPC模型预测控制器]4_数学建模推导--Matlab代码详解_哔哩哔哩_bilibili的学习分享如下: 一.研究目的在约束条件(物理限制)下达到最优的系 ...
预测电流FCS-MPC模型预测电流控制
预测电流FCS-MPC模型预测电流控制级联(链式):5个单元下载地址[https://download.csdn.net/download/a_zxswer/20009112](https://d ...
MPC模型预测控制原理和Matlab以及Python代码实现
MPC模型预测控制原理和代码一. 介绍模型预测控制(MPC)原理简要解释一下最优控制最优控制的目标是在一定的约束条件下达到最优的系统表现,那么要让系统达到最优表现,一般是通过定义损失函数J,通过最 ...
控制算法-模型预测MPC
本文记录一下MPC控制算法的学习过程和自己的理解,初步接触控制算法,理解肯定不是很完善,重在记录思考的过程. 背景随着自动驾驶技术以及机器人控制技术的不断发展及逐渐火热,模型预测控制(MPC)算法作 ...
ADRC，自抗扰控制器，扩张状态观测器，ESO，模型预测控制算法MPC
ADRC,自抗扰控制器,扩张状态观测器,ESO,模型预测控制算法MPC,自适应模型预测控制算法,时变模型预测控制算法,H无穷算法,混合灵敏度,鲁棒控制算法,四旋翼,直升机,控制算法设计,仿真模型,算法 ...
模型预测控制器(MPC)系列: 2.求解车辆横向控制中的MPC
求解MPC: 在滚动时间窗内建立并求解QP问题该小节旨在根据上一节得到的离散车辆横向动力学模型,在滚动时间窗内建立并求解QP问题,实现MPC横向控制. 在开始求解之前,我们需要回答以下两个问题: 1 ...

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