1.介绍

MPC(Model Predictive Control),即模型预测控制,是一种进阶的过程控制方法。相比于LQR和PID算法只能考虑输入输出变量的各种约束,MPC算法可以考虑状态空间变量的约束。可以应用于线性和非线性的系统中。
参考资料:

  1. MPC教学
  2. 421施公队的MPC教学视频

2.MPC四要素

(1) 模型:采用阶跃响应

(2) 预测

(3) 滚动优化:

二次规划
J ( x ) = x 2 + b x + c J\left(x\right)=x^2+bx+c J(x)=x2+bx+c
极值点
J ′ ( x 0 ) = 0 J^\prime\left(x_0\right)=0 J(x0)=0
(4) 误差补偿

参数定义

P:预测步长
y ( k + 1 ) , y ( k + 2 ) … y ( k + P ) y\left(k+1\right),y\left(k+2\right)\ldots y(k+P) y(k+1),y(k+2)y(k+P)
M:控制步长
  ∆ u ( k ) , ∆ u ( k + 1 ) … ∆ u ( k + M ) \ ∆u(k),∆u(k+1)…∆u(k+M)  u(k),u(k+1)u(k+M)

3.推导过程

(1)模型

假定在T,2T…PT时刻模型对应的输出为
a 1 , a 2 … a p a_1,a_2\ldots a_p a1,a2ap
则根据线性系统的叠加原理:
y ( k ) = a 1 ∗ u ( k − 1 ) + a 2 ∗ u ( k − 2 ) + a M ∗ u ( k − M ) y\left(k\right)=a_1\ast u\left(k-1\right)+a_2\ast u\left(k-2\right)+a_M\ast u(k-M) y(k)=a1u(k1)+a2u(k2)+aMu(kM)

增量式:
∆ y ( k ) = a 1 ∗ ∆ u ( k − 1 ) + a 2 ∗ ∆ u ( k − 2 ) + a M ∗ ∆ u ( k − M ) ∆y\left(k\right)=a_1\ast ∆u\left(k-1\right)+a_2\ast ∆u\left(k-2\right)+a_M\ast ∆u(k-M) y(k)=a1u(k1)+a2u(k2)+aMu(kM)

(2)预测:共预测P个时刻

∆ y ( k + 1 ) = a 1 ∗ ∆ u ( k ) ∆y(k+1)=a_1\ast∆u(k) y(k+1)=a1u(k)

∆ y ( k + 2 ) = a 1 ∗ ∆ u ( k ) + 1 + a 2 ∗ ∆ u ( k ) ∆y(k+2)=a_1*∆u(k)+1+a_2*∆u(k) y(k+2)=a1u(k)+1+a2u(k)

∆ y ( k + P ) = a 1 ∗ ∆ u ( k + P − 1 ) + a 2 ∗ ∆ u ( k + P − 2 ) … + a M ∗ ∆ u ( k + P − M ) ∆y(k+P)=a_1*∆u(k+P-1)+a_2*∆u(k+P-2)…+a_M*∆u(k+P-M) y(k+P)=a1u(k+P1)+a2u(k+P2)+aMu(k+PM)

转为矩阵表示:
∑ i = 1 M a i ∗ ∆ u ( k + P − i ) \sum_{i=1}^{M}{a_i\ast}∆u(k+P-i) i=1Maiu(k+Pi)

新的预测输出 = 原预测输出 + 增量输入
Y ^ 0 = A ∗ ∆ u + Y 0 [ Y ^ 0 ( k + 1 ) Y ^ 0 ( k + 2 ) ⋮ Y ^ 0 ( k + P ) ] = [ a 1 0 a 2 a 1 ⋯ 0 ⋯ 0 ⋮ ⋮ a P a P − 1 ⋱ ⋮ ⋯ a P − M + 1 ] [ ∆ u ( k ) ∆ u ( k + 1 ) ⋮ ∆ u ( k + M ) ] + [ Y 0 ( k + 1 ) Y 0 ( k + 2 ) ⋮ Y 0 ( k + P ) ] {\hat{Y}}_0=A\ast∆u+Y_0 \\ \left[\begin{matrix}\begin{matrix}{\hat{Y}}_0(k+1)\\{\hat{Y}}_0(k+2)\\\end{matrix}\\\begin{matrix}\vdots\\{\hat{Y}}_0(k+P)\\\end{matrix}\\\end{matrix}\right]=\left[\begin{matrix}\begin{matrix}a_1&0\\a_2&a_1\\\end{matrix}&\begin{matrix}\cdots&0\\\cdots&0\\\end{matrix}\\\begin{matrix}\vdots&\vdots\\a_P&a_{P-1}\\\end{matrix}&\begin{matrix}\ddots&\vdots\\\cdots&a_{P-M+1}\\\end{matrix}\\\end{matrix}\right] \left[\begin{matrix}\begin{matrix}∆u(k)\\∆u(k+1)\\\end{matrix}\\\begin{matrix}\vdots\\∆u(k+M)\\\end{matrix}\\\end{matrix}\right]+\left[\begin{matrix}\begin{matrix}Y_0(k+1)\\Y_0(k+2)\\\end{matrix}\\\begin{matrix}\vdots\\Y_0(k+P)\\\end{matrix}\\\end{matrix}\right] Y^0=Au+Y0Y^0(k+1)Y^0(k+2)Y^0(k+P)=a1a20a1aPaP100aPM+1u(k)u(k+1)u(k+M)+Y0(k+1)Y0(k+2)Y0(k+P)

(3)滚动优化

注:求解出来的∆u为M*1的向量,取第一个,即 ∆ u 0 ∆u_0 u0作为输入的增量。

  1. 参考轨迹
    使用一阶滤波处理:
    ω ( k + i ) = α i y ( k ) + ( 1 − α i ) y r ( k ) \omega\left(k+i\right)=\alpha^iy\left(k\right)+(1-\alpha^i)y_r(k) ω(k+i)=αiy(k)+(1αi)yr(k)

  2. 代价函数(与最优化控制相似)
    形式一

    目标一:离目标越近越好
    J 1 = ∑ i = 1 P [ y ( k + i ) − ω ( k + i ) ] 2 J_1=\sum_{i=1}^{P}\left[y\left(k+i\right)-\omega\left(k+i\right)\right]^2 J1=i=1P[y(k+i)ω(k+i)]2

    目标二:能量消耗越小越好
    J 2 = i = ∑ i = 1 M [ ∆ u ( k + i − 1 ) ] 2 J_2=i=\sum_{i=1}^{M}[∆u(k+i-1)]^2 J2=i=i=1M[u(k+i1)]2

    合体得到总的代价函数
    J = q J 1 + r J 2 J=qJ_1+rJ_2 J=qJ1+rJ2

    q和r都是权重系数,它们的相对大小表示更看重哪个指标,越大说明越看重。

    矩阵形式:
    J = ( r − ω ) T Q ( R − ω ) + ∆ u T R ∆ u J=\left(r-\omega\right)^TQ\left(R-\omega\right)+∆u^TR∆u\\ J=(rω)TQ(Rω)+uTRu
    其中,Q和R大多数为对角阵(一般不考虑协方差矩阵,即不同变量之间相互影响的情况)
    [ q 1 0 . . . 0 0 q 2 . . . 0 ⋮ ⋮ ⋱ ⋮ 0 0 0 q n ] [ r 1 0 . . . 0 0 r 2 . . . 0 ⋮ ⋮ ⋱ ⋮ 0 0 0 r n ] \begin{bmatrix} q_1& 0 & ... & 0\\ 0& q_2 & ...& 0\\ \vdots & \vdots &\ddots &\vdots \\ 0&0 & 0&q_n \end{bmatrix} \begin{bmatrix} r_1& 0 & ... & 0\\ 0& r_2 & ...& 0\\ \vdots & \vdots &\ddots &\vdots \\ 0&0 & 0&r_n \end{bmatrix} q1000q20......000qnr1000r20......000rn
    最优解:
    ∂ J ( Δ u ) ∂ ( Δ u ) = 0 → Δ u = ( A T Q A + R ) − 1   A T ( ω − Y 0 ) \frac{\partial \mathrm{J}(\Delta u)}{\partial(\Delta u)}=0 \rightarrow \Delta u=\left(\mathrm{A}^{\mathrm{T}} \mathrm{QA}+\mathrm{R}\right)^{-1} \mathrm{~A}^{\mathrm{T}}\left(\omega-\mathrm{Y}_{0}\right) (Δu)J(Δu)=0Δu=(ATQA+R)1 AT(ωY0)
    形式二

    • 基于 u k , u k + 1 . . . u k + N − 1 进 行 最 优 化 u_k,u_{k+1}...u_{k+N-1}进行最优化 uk,uk+1...uk+N1
      J = ∑ k N − 1 ( E k T Q E k + u k T R u k ) + E N T F E N ⏟ t e r m i n a l   p r e d i c t   s t a t u s J = \sum_k^{N-1}(E_k^TQE_k + u_k^TRu_k)+\underbrace{E_N^TFE_N} _{terminal \ predict \ status} J=kN1(EkTQEk+ukTRuk)+terminal predict status ENTFEN

(4)反馈校正,误差补偿

  • k时刻,预测P个输出:
    Y ^ 0 ( k + 1 ) , Y ^ 0 ( k + 2 ) . . . Y ^ 0 ( k + P ) {\hat{Y}}_0\left(k+1\right),{\hat{Y}}_0\left(k+2\right)...{\hat{Y}}_0\left(k+P\right) Y^0(k+1)Y^0(k+2)...Y^0(k+P)

  • k+1时刻,当前输出为:
    y ( k + 1 ) y(k+1) y(k+1)

  • 误差:
    e ( k + 1 ) =   y ( k + 1 ) − Y ^ 0 ( k + 1 ) e\left(k+1\right)=\ y\left(k+1\right)-{\hat{Y}}_0\left(k+1\right) e(k+1)= y(k+1)Y^0(k+1)

h:补偿系数,取值范围0-1,习惯0.5

补偿:
y c o r ( k + 1 ) =   Y ^ 0 ( k + 1 ) + h 1 ∗ e ( k + 1 ) y c o r ( k + 2 ) =   Y ^ 0 ( k + 2 ) + h 2 ∗ e ( k + 1 ) ⋮ y c o r ( k + P ) =   Y ^ 0 ( k + P ) + h P ∗ e ( k + 1 ) {y}_{cor}\left(k+1\right)=\ {\hat{Y}}_0\left(k+1\right)+h_1\ast e\left(k+1\right) \\y_{cor}\left(k+2\right)=\ {\hat{Y}}_0\left(k+2\right)+h_2\ast e\left(k+1\right) \\\vdots \\y_{cor}\left(k+P\right)=\ {\hat{Y}}_0\left(k+P\right)+h_P\ast e\left(k+1\right) ycor(k+1)= Y^0(k+1)+h1e(k+1)ycor(k+2)= Y^0(k+2)+h2e(k+1)ycor(k+P)= Y^0(k+P)+hPe(k+1)
矩阵表示:
Y c o r = Y ^ 0 + H ∗ e ( k + 1 ) Y_{cor}={\hat{Y}}_0+H\ast e(k+1) Ycor=Y^0+He(k+1)

滚动更新预测输出:

K+1时刻:
Y ^ 0 = S ∗ Y c o r {\hat{Y}}_0=S\ast Y_{cor} Y^0=SYcor

S:移位矩阵(PxP)
[ 0 1 0 0 ⋯ 0 ⋯ ⋮ ⋮ ⋮ 0 0 ⋱ 1 ⋯ 1 ] \left[\begin{matrix}\begin{matrix}0&1\\0&0\\\end{matrix}&\begin{matrix}\cdots&0\\\cdots&\vdots\\\end{matrix}\\\begin{matrix}\vdots&\vdots\\0&0\\\end{matrix}&\begin{matrix}\ddots&1\\\cdots&1\\\end{matrix}\\\end{matrix}\right] 001000011

补充:代价函数的推导

在这里插入图片描述
在这里插入图片描述

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐