MPC轨迹跟踪控制器推导及Simulink验证

MPC轨迹跟踪控制器推导及Simulink验证

MPC的特点

  1. 多变量控制

一步都需要解决一个优化问题,特别是当系统变量和预测范围增大时,计算量会显著增加。
7. 需要精确模型:虽然MPC可以设计来应对模型不确定性,但其性能很大程度上依赖于模型的准确性。
8. 实现复杂:MPC的实现相对复杂,需要较高的专业知识和经验。
9. 与LQR相比:

  • 计算较复杂:LQR控制器的求解相对简单,只需要解决一个线性二次调节器问题。
  • 考虑约束:传统的LQR不考虑输入和状态约束。
  • 灵活性:LQR通常优化一个固定的二次型目标,不如MPC灵活。

MPC轨迹跟踪控制器推导

一 系统离散化

对连续线性状态方程:
{ x ˙ ( t ) = A x ( t ) + B u ( t ) y ( t ) = C x ( t ) + D u ( t ) \left\{ \begin{aligned} \dot{x}(t) &= Ax(t) + Bu(t) \\ y(t) &= Cx(t) + Du(t) \end{aligned} \right. {x˙(t)y(t)=Ax(t)+Bu(t)=Cx(t)+Du(t)
其离散化状态方程为:
{ x [ k + 1 ] = A d x [ k ] + B d u [ k ] y [ k ] = C d x [ k ] + D d u [ k ] \left\{ \begin{aligned} & x[k+1] = A_dx[k] + B_du[k] \\ &y [k] = C_dx[k] + D_du[k] \end{aligned} \right. {x[k+1]=Adx[k]+Bdu[k]y[k]=Cdx[k]+Ddu[k]
其中,
{ A d = e A T B d = A − 1 ( e A T − I ) B C d = C D d = D \left\{ \begin{aligned} A_d &= e^{AT}\\ B_d &= A^{-1}(e^{AT}-I)B \\ C_d &= C\\ D_d &= D \end{aligned} \right. AdBdCdDd=eAT=A1(eATI)B=C=D

二 预测区间状态和变量推导

X = [ x ( k ∣ k ) x ( k + 1 ∣ k ) x ( k + 2 ∣ k ) ⋮ x ( k + N ∣ k ) ] ( N + 1 ) n × 1 X = \begin{bmatrix} x(k|k) \\ x(k+1|k) \\ x(k+2|k) \\ \vdots \\ x(k+N|k) \end{bmatrix}_{(N+1)n\times1} X= x(kk)x(k+1∣k)x(k+2∣k)x(k+Nk) (N+1)n×1

U = [ u ( k ∣ k ) u ( k + 1 ∣ k ) u ( k + 2 ∣ k ) ⋮ u ( k + N − 1 ∣ k ) ] N p × 1 U = \begin{bmatrix} u(k|k) \\ u(k+1|k) \\ u(k+2|k) \\ \vdots \\ u(k+N-1|k) \end{bmatrix} _{Np\times1} U= u(kk)u(k+1∣k)u(k+2∣k)u(k+N1∣k) Np×1
根据离散化的状态空间方程
{ x ( k ∣ k ) = x k x ( k + 1 ∣ k ) = A d x ( k ∣ k ) + B d u ( k ∣ k ) = A d x k + B d u ( k ∣ k ) x ( k + 2 ∣ k ) = A d x ( k + 1 ∣ k ) + B d u ( k + 1 ∣ k ) = A d 2 x k + A d B d u ( k ∣ k ) + u ( k + 1 ∣ k ) ⋮ x ( k + N ∣ k ) = A d N x k + A d N − 1 B d u ( k ∣ k ) + ⋯ + u ( k + N − 1 ∣ k ) \left\{ \begin{aligned} & x(k|k) = x_k\\ & x(k+1|k) = A_dx(k|k)+B_du(k|k) = A_dx_k+B_du(k|k)\\ & x(k+2|k) = A_dx(k+1|k)+B_du(k+1|k) = A_d^2x_k+A_dB_du(k|k)+u(k+1|k)\\ \vdots\\ & x(k+N|k) = A_d^Nx_k+A_d^{N-1}B_du(k|k)+\dots +u(k+N-1|k) \end{aligned} \right. x(kk)=xkx(k+1∣k)=Adx(kk)+Bdu(kk)=Adxk+Bdu(kk)x(k+2∣k)=Adx(k+1∣k)+Bdu(k+1∣k)=Ad2xk+AdBdu(kk)+u(k+1∣k)x(k+Nk)=AdNxk+AdN1Bdu(kk)++u(k+N1∣k)
化简得:
X k = M ( N + 1 ) n × n x k + C ( N + 1 ) n × N p U k X_k = M_{(N+1)n\times n} x_k+C_{(N+1)n\times Np}U_k Xk=M(N+1)n×nxk+C(N+1)n×NpUk
其中:
M ( N + 1 ) n × n = [ I A d A d 2 ⋮ A d N ] C ( N + 1 ) n × N p = [ O B d A d B d B d ⋮ ⋮ ⋱ A d N − 1 B d A d N − 2 B d … B d ] M_{(N+1)n\times n} = \begin{bmatrix} I\\ A_d \\ A_d^2 \\ \vdots \\ A_d^N \end{bmatrix} \\ C_{(N+1)n\times Np}= \begin{bmatrix} O\\ B_d \\ A_dB_d & B_d \\ \vdots &\vdots &\ddots \\ A_d^{N-1}B_d & A_d^{N-2}B_d & \dots & B_d \end{bmatrix} M(N+1)n×n= IAdAd2AdN C(N+1)n×Np= OBdAdBdAdN1BdBdAdN2BdBd

三 代价函数推导

定义代价函数为状态误差加权和、控制输入加权和以及终端误差加权三者之和
J = ∑ k N − 1 ( E k T Q E k + u T k R u k ) + E N T F E N J = \sum_{k}^{N-1}(E_{k}^{T}QE_k + u_{T}^{k}Ru_{k})+E_{N}^{T}FE_N J=kN1(EkTQEk+uTkRuk)+ENTFEN
其中, E k = y − r E_{k}=y-r Ek=yr, y y y为输出变量,r为参考向量,在倒立摆中,取 y = x y = x y=x,故 E k = x − r e f E_k = x-ref Ek=xref,化简代价函数为:
J = [ E ( k ∣ k ) E ( k + 1 ∣ k ) E ( k + 2 ∣ k ) ⋮ E ( k + N ∣ k ) ] T Q ‾ [ E ( k ∣ k ) E ( k + 1 ∣ k ) E ( k + 2 ∣ k ) ⋮ E ( k + N ∣ k ) ] + [ u ( k ∣ k ) u ( k + 1 ∣ k ) u ( k + 2 ∣ k ) ⋮ u ( k + N − 1 ∣ k ) ] T R ‾ [ u ( k ∣ k ) u ( k + 1 ∣ k ) u ( k + 2 ∣ k ) ⋮ u ( k + N − 1 ∣ k ) ] J = \begin{bmatrix} E(k|k)\\ E(k+1|k)\\ E(k+2|k)\\ \vdots\\ E(k+N|k) \end{bmatrix}^{T}\overline{Q} \begin{bmatrix} E(k|k)\\ E(k+1|k)\\ E(k+2|k)\\ \vdots\\ E(k+N|k) \end{bmatrix}+ \begin{bmatrix} u(k|k)\\ u(k+1|k) \\ u(k+2|k)\\ \vdots\\ u(k+N-1|k) \end{bmatrix}^{T}\overline{R} \begin{bmatrix} u(k|k)\\ u(k+1|k) \\ u(k+2|k)\\ \vdots\\ u(k+N-1|k) \end{bmatrix} J= E(kk)E(k+1∣k)E(k+2∣k)E(k+Nk) TQ E(kk)E(k+1∣k)E(k+2∣k)E(k+Nk) + u(kk)u(k+1∣k)u(k+2∣k)u(k+N1∣k) TR u(kk)u(k+1∣k)u(k+2∣k)u(k+N1∣k)
其中:

{ E ( k + i ∣ k ) = x ( k + i ∣ k ) − r e f , i = 0 , 1 , 2 , … , N Q ‾ = d i a g ( Q , Q , Q , … , F ) R ‾ = d i a g ( R , R , R , … , R ) \left\{ \begin{aligned} & E(k+i|k) = x(k+i|k)-ref, i= 0,1,2,\dots,N\\ & \overline{Q} = diag(Q,Q,Q,\dots,F)\\ & \overline{R} = diag(R,R,R,\dots,R) \end{aligned} \right. E(k+ik)=x(k+ik)ref,i=0,1,2,,NQ=diag(Q,Q,Q,,F)R=diag(R,R,R,,R)
进一步展开得:
J = X k T Q ‾ X k − X T Q ‾ R e f − R e f T Q ‾ X + R e f T Q ‾ R e f + U k T R ‾ U k J = X_k^T\overline{Q}X_k - X^T\overline{Q}R_{ef} - R_{ef}^T\overline{Q}X + R_{ef}^T\overline{Q}R_{ef} + U_k^T\overline{R} U_k J=XkTQXkXTQRefRefTQX+RefTQRef+UkTRUk
易知:
{ R e f = [ r e f ; r e f ; …   ; r e f ] X T Q ‾ R e f = R e f T Q ‾ X X k = M x k + C U k \left\{ \begin{aligned} & R_{ef} =[ref;ref;\dots;ref] \\ & X^T\overline{Q}R_{ef} = R_{ef}^T\overline{Q}X \\ & X_k = Mx_k+CU_k \\ \end{aligned} \right. Ref=[ref;ref;;ref]XTQRef=RefTQXXk=Mxk+CUk
继续化简:
J = ( M x k + C U k ) T Q ‾ ( M x k + C U k ) − 2 ∗ ( M x k + C U k ) T Q ‾ R e f + R e f T Q ‾ R e f + U k T R ‾ U k = x k T M T Q ‾ M x k + 2 x k T M T Q ‾ C U k + U k T C T Q ‾ C U k − 2 x k T M T Q ‾ R e f − 2 U k T C T Q ‾ R e f + R e f T Q ‾ R e f + U k T R ‾ U k = x k T M T Q ‾ M x k + 2 x k T M T Q ‾ C U k − 2 x k T M T Q ‾ R e f − 2 U k T C T Q ‾ R e f + R e f T Q ‾ R e f + U k T ( C T Q ‾ C + R ‾ ) U k = x k T G x k + 2 x k T L U k − 2 x k T P R e f − 2 U k T S R e f + R e f T Q ‾ R e f + U k T T U k = x k T G x k − 2 R e f T P T x k + R e f T Q ‾ R e f + 2 ( x k T L − R e f T S T ) U k + U k T T U k \begin{aligned} J & = (Mx_k+CU_k)^T\overline{Q}(Mx_k+CU_k) - 2*(Mx_k+CU_k)^T\overline{Q}R_{ef} + R_{ef}^T\overline{Q}R_{ef} + U_k^T\overline{R} U_k \\[2ex] & = x_k^TM^T\overline{Q}Mx_k + 2x_k^TM^T\overline{Q}CU_k + U_k^TC^T\overline{Q}CU_k - 2x_k^TM^T\overline{Q}R_{ef} \\ & \hspace{1.5em} -2U_k^TC^T\overline{Q}R_{ef} + R_{ef}^T\overline{Q}R_{ef} + U_k^T\overline{R} U_k \\[2ex] & = x_k^TM^T\overline{Q}Mx_k + 2x_k^TM^T\overline{Q}CU_k - 2x_k^TM^T\overline{Q}R_{ef} - 2U_k^TC^T\overline{Q}R_{ef} + \\ & \hspace{1.5em} R_{ef}^T\overline{Q}R_{ef} + U_k^T(C^T\overline{Q}C+\overline{R})U_k \\[2ex] & = x_k^TGx_k + 2x_k^TLU_k - 2x_k^TPR_{ef} - 2U_k^TSR_{ef} + R_{ef}^T\overline{Q}R_{ef} + U_k^TTU_k\\[2ex] % & = x_k^TGx_k - 2R_{ef}^TLx_k + 2x_k^THU_k - 2R_{ef}^TPU_k + R_{ef}^T\overline{Q}R_{ef} + U_k^TSU_k\\[2ex] & = x_k^TGx_k - 2R_{ef}^TP^Tx_k + R_{ef}^T\overline{Q}R_{ef} + 2(x_k^TL - R_{ef}^TS^T)U_k + U_k^TTU_k\\ \end{aligned} J=(Mxk+CUk)TQ(Mxk+CUk)2(Mxk+CUk)TQRef+RefTQRef+UkTRUk=xkTMTQMxk+2xkTMTQCUk+UkTCTQCUk2xkTMTQRef2UkTCTQRef+RefTQRef+UkTRUk=xkTMTQMxk+2xkTMTQCUk2xkTMTQRef2UkTCTQRef+RefTQRef+UkT(CTQC+R)Uk=xkTGxk+2xkTLUk2xkTPRef2UkTSRef+RefTQRef+UkTTUk=xkTGxk2RefTPTxk+RefTQRef+2(xkTLRefTST)Uk+UkTTUk
其中:

{ G = M T Q ‾ M L = M T Q ‾ C P = M T Q ‾ S = C T Q ‾ T = C T Q ‾ C + R ‾ \left\{ \begin{aligned} & G = M^T\overline{Q}M\\ & L = M^T\overline{Q}C \\ & P = M^T\overline{Q}\\ & S = C^T\overline{Q}\\ & T = C^T\overline{Q}C+\overline{R}\\ \end{aligned} \right. G=MTQML=MTQCP=MTQS=CTQT=CTQC+R

四 优化求解

通过3. 代价函数推导,采用二次优化方法使得代价 J J J最小,进而进行控制,定义优化问题如下:
min ⁡ U k J = U k T T U k + 2 ( x k T L − R e f T S T ) U k subject to U m i n ≤ U k ≤ U m a x \begin{aligned} \underset{U_k}\min \quad & J = U_k^TTU_k+2(x_k^TL-R_{ef}^TS^T)U_k \\ \text{subject to}\quad & U_{min}\leq U_k \leq U_{max}\\ \end{aligned} Ukminsubject toJ=UkTTUk+2(xkTLRefTST)UkUminUkUmax
在MATLAB中可使用函数quadprog(quadprog函数介绍)求解:

根据MATLAB中quadprog要求,定义:
{ H = 2 T = 2 ( C T Q ‾ C + R ‾ ) f = 2 ( x k T L − R e f T S T ) T A = [ ] b = [ ] A e q = [ ] b e q = [ ] l b = U m i n u b = U m a x \left\{ \begin{aligned} H & = 2T = 2(C^T\overline{Q}C+\overline{R})\\ f & = 2(x_k^TL-R_{ef}^TS^T)^T\\ A & = [\quad]\\ b & = [\quad]\\ A_{eq} & = [\quad]\\ b_{eq} & = [\quad]\\ lb & = U_{min}\\ ub & = U_{max}\\ \end{aligned} \right. HfAbAeqbeqlbub=2T=2(CTQC+R)=2(xkTLRefTST)T=[]=[]=[]=[]=Umin=Umax
调用 quadprog(H,f,A,b,Aeq,beq,lb,ub) 即可求解代价函数 J J J取最小值时,最优控制输入 u k u_k uk

基于MPC的倒立摆控制系统

采用这篇文章建立的倒立摆模型进行仿真测试。
simulink中搭建控制仿真模型如下:
simulink模型

输入计算出的反馈增益 K K K仿真结果如下:
仿真结果

相关资料

MATLAB代码和simulink模型点击这里

NMPC轨迹跟踪控制器推导及Simulink验证点击这里

欢迎大家与我交流!!!

Reference:

  1. 推导倒立摆模型并使用LQR进行位置跟踪控制:https://blog.csdn.net/weixin_55249340/article/details/140421091?spm=1001.2014.3001.5501
  2. MATLAB quadprog函数:https://ww2.mathworks.cn/help/optim/ug/quadprog_zh_CN.html?requestedDomain=cn
  • 11
    点赞
  • 30
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

有头发的垃圾猿

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值