MPC概念简介

Medel predictive control(MPC)是一种预测控制，可以依据未来的信息来控制当下的输入，本文主要介绍一种线性MPC（系统和约束都是线性）的实现方法，大致思路是将控制问题进行数学建模，整理成二次规划（quadratic programming，QP）的形式，而后求解。总体MPC的实现包括以下三步：

估计系统的初始状态
通过优化，计算得到未来NpN_pNp个点的输入序列
只输入第一个控制序列，而后将新状态视作初始状态，循环执行

MPC简单公式推导

系统方程推导

假设系统的离散状态方程为：x(k+1)=Amx(k)+Bumy(k)=Cmx(k)x(k+1) = A_{m}x(k)+Bu_{m} \\ y(k) = C_{m}x(k)x(k+1)=Amx(k)+Bumy(k)=Cmx(k)
则可以先将状态方程整理成增广形式（augmented form）,增广形式更加适用于轨迹跟踪，如下所示：[Δxm(k+1)y(k+1)]⏞x(k+1)=[AmomTCmAm1]⏞A[Δxm(k)y(k)]⏞x(k)+[BmCmBm]⏞BΔu(k)y(k)=[om1]⏞C[Δxm(k)y(k)],\begin{array}{l} \overbrace{\left[\begin{array}{c} \Delta x_{m}(k+1) \\ y(k+1) \end{array}\right]}^{x(k+1)}=\overbrace{\left[\begin{array}{cc} A_{m} & o_{m}^{T} \\ C_{m} A_{m} & 1 \end{array}\right]}^{A} \overbrace{\left[\begin{array}{c} \Delta x_{m}(k) \\ y(k) \end{array}\right]}^{x(k)}+\overbrace{\left[\begin{array}{c} B_{m} \\ C_{m} B_{m} \end{array}\right]}^{B} \Delta u(k)\\ y(k)=\overbrace{\left[\begin{array}{ll} o_{m} & 1 \end{array}\right]}^{C}\left[\begin{array}{c} \Delta x_{m}(k) \\ y(k) \end{array}\right], \end{array}[Δxm(k+1)y(k+1)]x(k+1)=[AmCmAmomT1]A[Δxm(k)y(k)]x(k)+[BmCmBm]BΔu(k)y(k)=[om1]C[Δxm(k)y(k)],
随后可以根据初始状态x(ki)x(k_i)x(ki)和系统的状态方程得到在预测范围NpN_pNp的范围内所有的状态变量，如下所示：x(ki+1∣ki)=Ax(ki)+BΔu(ki)x(ki+2∣ki)=Ax(ki+1∣ki)+BΔu(ki+1)=A2x(ki)+ABΔu(ki)+BΔu(ki+1)⋮x(ki+Np∣ki)=ANpx(ki)+ANp−1BΔu(ki)+ANp−2BΔu(ki+1)+…+ANp−NcBΔu(ki+Nc−1).\begin{aligned} x\left(k_{i}+1 \mid k_{i}\right) &=A x\left(k_{i}\right)+B \Delta u\left(k_{i}\right) \\ x\left(k_{i}+2 \mid k_{i}\right) &=A x\left(k_{i}+1 \mid k_{i}\right)+B \Delta u\left(k_{i}+1\right) \\ &=A^{2} x\left(k_{i}\right)+A B \Delta u\left(k_{i}\right)+B \Delta u\left(k_{i}+1\right) \\ & \vdots \\ x\left(k_{i}+N_{p} \mid k_{i}\right) &=A^{N_{p}} x\left(k_{i}\right)+A^{N_{p}-1} B \Delta u\left(k_{i}\right)+A^{N_{p}-2} B \Delta u\left(k_{i}+1\right) \\ &+\ldots+A^{N_{p}-N_{c}} B \Delta u\left(k_{i}+N_{c}-1\right) . \end{aligned}x(ki+1∣ki)x(ki+2∣ki)x(ki+Np∣ki)=Ax(ki)+BΔu(ki)=Ax(ki+1∣ki)+BΔu(ki+1)=A2x(ki)+ABΔu(ki)+BΔu(ki+1)⋮=ANpx(ki)+ANp−1BΔu(ki)+ANp−2BΔu(ki+1)+…+ANp−NcBΔu(ki+Nc−1).
而所有输出变量（观测变量）的表达如下所示：
y(ki+1∣ki)=CAx(ki)+CBΔu(ki)y(ki+2∣ki)=CA2x(ki)+CABΔu(ki)+CBΔu(ki+1)y(ki+3∣ki)=CA3x(ki)+CA2BΔu(ki)+CABΔu(ki+1)+CBΔu(ki+2)⋮y(ki+Np∣ki)=CANpx(ki)+CANp−1BΔu(ki)+CANp−2BΔu(ki+1)+…+CANp−NcBΔu(ki+Nc−1)\begin{aligned} y\left(k_{i}+1 \mid k_{i}\right) &=C A x\left(k_{i}\right)+C B \Delta u\left(k_{i}\right) \\ y\left(k_{i}+2 \mid k_{i}\right) &=C A^{2} x\left(k_{i}\right)+C A B \Delta u\left(k_{i}\right)+C B \Delta u\left(k_{i}+1\right) \\ y\left(k_{i}+3 \mid k_{i}\right) &=C A^{3} x\left(k_{i}\right)+C A^{2} B \Delta u\left(k_{i}\right)+C A B \Delta u\left(k_{i}+1\right) \\ &+C B \Delta u\left(k_{i}+2\right) \\ & \vdots \\ y\left(k_{i}+N_{p} \mid k_{i}\right) &=C A^{N_{p}} x\left(k_{i}\right)+C A^{N_{p}-1} B \Delta u\left(k_{i}\right)+C A^{N_{p}-2} B \Delta u\left(k_{i}+1\right) \\ &+\ldots+C A^{N_{p}-N_{c}} B \Delta u\left(k_{i}+N_{c}-1\right) \end{aligned}y(ki+1∣ki)y(ki+2∣ki)y(ki+3∣ki)y(ki+Np∣ki)=CAx(ki)+CBΔu(ki)=CA2x(ki)+CABΔu(ki)+CBΔu(ki+1)=CA3x(ki)+CA2BΔu(ki)+CABΔu(ki+1)+CBΔu(ki+2)⋮=CANpx(ki)+CANp−1BΔu(ki)+CANp−2BΔu(ki+1)+…+CANp−NcBΔu(ki+Nc−1)
根据上式，将复杂的线性方程写成矩阵形式：
Y=Fx(ki)+ΦΔUY = Fx(k_i)+\Phi\Delta UY=Fx(ki)+ΦΔU
矩阵表达如下所示：F=[CACA2CA3⋮CANp];Φ=[B00…0ABB0…0A2BABB…0⋮ANp−1BANp−2BANp−3B…ANp−NcB]F=\left[\begin{array}{c} C A \\ C A^{2} \\ C A^{3} \\ \vdots \\ C A^{N_{p}} \end{array}\right] ; \Phi=\left[\begin{array}{ccccc} B & 0 & 0 & \ldots & 0 \\ A B & B & 0 & \ldots & 0 \\ A^{2} B & A B & B & \ldots & 0 \\ \vdots & & & & \\ A^{N_{p}-1} B & A^{N_{p}-2} B & A^{N_{p}-3} B & \ldots & A^{N_{p}-N_{c}} B \end{array}\right]F=⎣⎡CACA2CA3⋮CANp⎦⎤;Φ=⎣⎡BABA2B⋮ANp−1B0BABANp−2B00BANp−3B…………000ANp−NcB⎦⎤

注意：式中的A,B,C是增广矩阵中的系数，不是原方程的系数

而后进入优化部分，定义待追踪的向量是RsR_sRs，给输入的误差矩阵是Rˉ\bar{R}Rˉ，可以将代价函数写成以下格式（大多数书中对误差也有权重矩阵，但是通过调整Rˉ\bar{R}Rˉ也可以简介调整两项之间的权重差距，因此没有设置追踪误差权重）：
J=(Rs−Y)T(Rs−Y)+ΔUTRˉΔUJ=\left(R_{s}-Y\right)^{T}\left(R_{s}-Y\right)+\Delta U^{T} \bar{R} \Delta UJ=(Rs−Y)T(Rs−Y)+ΔUTRˉΔU，带入前面得到的数学方程，可以将代价函数整理成以下形式：J=(Rs−Fx(ki))T(Rs−Fx(ki))−2ΔUTΦT(Rs−Fx(ki))+ΔUT(ΦTΦ+Rˉ)ΔU.J=\left(R_{s}-F x\left(k_{i}\right)\right)^{T}\left(R_{s}-F x\left(k_{i}\right)\right)-2 \Delta U^{T} \Phi^{T}\left(R_{s}-F x\left(k_{i}\right)\right)+\Delta U^{T}\left(\Phi^{T} \Phi+\bar{R}\right) \Delta U .J=(Rs−Fx(ki))T(Rs−Fx(ki))−2ΔUTΦT(Rs−Fx(ki))+ΔUT(ΦTΦ+Rˉ)ΔU.
至此可以使用二次规划的包，得到优化的ΔU\Delta UΔU

约束推导

约束通常有三种形式，分别是约束控制变化速率，约束控制幅值和约束输出变量（观测变量），可以用以下矩阵统一表达：
[M1M2M3]ΔU≤[N1N2N3]\left[\begin{array}{l} M_{1} \\ M_{2} \\ M_{3} \end{array}\right] \Delta U \leq\left[\begin{array}{l} N_{1} \\ N_{2} \\ N_{3} \end{array}\right]⎣⎡M1M2M3⎦⎤ΔU≤⎣⎡N1N2N3⎦⎤式中：
M1=[−C2C2];N1=[−Umin⁡+C1u(ki−1)Umax⁡−C1u(ki−1)];M2=[−II];N2=[−ΔUmin⁡ΔUmax⁡];M3=[−ΦΦ];N3=[−Ymin⁡+Fx(ki)Ymax⁡−Fx(ki)].\begin{array}{c} M_{1}=\left[\begin{array}{c} -C_{2} \\ C_{2} \end{array}\right] ; N_{1}=\left[\begin{array}{c} -U^{\min }+C_{1} u\left(k_{i}-1\right) \\ U^{\max }-C_{1} u\left(k_{i}-1\right) \end{array}\right] ; M_{2}=\left[\begin{array}{c} -I \\ I \end{array}\right] ; \\ N_{2}=\left[\begin{array}{c} -\Delta U^{\min } \\ \Delta U^{\max } \end{array}\right] ; M_{3}=\left[\begin{array}{c} -\Phi \\ \Phi \end{array}\right] ; N_{3}=\left[\begin{array}{c} -Y^{\min }+F x\left(k_{i}\right) \\ Y^{\max }-F x\left(k_{i}\right) \end{array}\right] . \end{array}M1=[−C2C2];N1=[−Umin+C1u(ki−1)Umax−C1u(ki−1)];M2=[−II];N2=[−ΔUminΔUmax];M3=[−ΦΦ];N3=[−Ymin+Fx(ki)Ymax−Fx(ki)].
M1M_1M1对应约束控制幅值，C1,C2C_1,C_2C1,C2分别对应全为1矩阵和下三角矩阵M2M_2M2对应约束控制幅值M3M_3M3约束输出变量（观测变量）

MPC实例与Matlab代码

main

以下代码从连续系统出发，而后离散化，根据自定义函数计算mpcgain，而后使用自定义二次优化函数优化，本例中对变量无约束。

clc
clear
close alldt = 0.1; % 采样时间Ac = [1 1;0 1];
Bc = [0;1];
Cc = [1 0];
Dc = [0];
%% 离散方程
[Ap,Bp,Cp,Dp]=c2dm(Ac,Bc,Cc,Dc,dt);Np=20;
Nc=4;[Phi,F,Phi_Phi,Phi_F,Phi_R,A_e, B_e,C_e] = mpcgain(Ap,Bp,Cp,Nc,Np);
[n,n_in]=size(B_e);
xm=[0;0];
Xf=zeros(n,1);
N_sim=100;
r=ones(N_sim,1);
u=0; % u(k-1) =0
y=0;
R_bar = 0.1*eye(Nc,Nc);
H = 2*(Phi_Phi + R_bar);
R_ref = ones(Np,1);
A_cons = 0; b = 0;for kk=1:N_sim
%     DeltaU=inv(Phi_Phi+R_bar)*(Phi_R*r(kk)-Phi_F*Xf);f = -2*Phi'*(R_ref-F*Xf);DeltaU=QPhild(H,f,A_cons,b);deltau=DeltaU(1,1);u=u+deltau;u1(kk)=u;y1(kk)=y;xm_old=xm;xm=Ap*xm+Bp*u;y=Cp*xm;Xf=[xm-xm_old;y];
end
% 将系统输出、控制量曲线显示出来：
k=0:(N_sim-1);
figure
subplot(211)
plot(k,y1)
xlabel('Sampling Instant')
legend('Output')
subplot(212)
plot(k,u1)
xlabel('Sampling Instant')
legend('Control')

mpcgain

说明，本函数和优化函数都来自参考文献的代码，我只进行了小修改，本来书中是使用求导的方法进行优化，所以本函数返回的矩阵并不都有用。

function [Phi,F,Phi_Phi,Phi_F,Phi_R,A_e, B_e,C_e]=mpcgain(Ap,Bp,Cp,Nc,Np)
[m1,n1]=size(Cp);
[n1,n_in]=size(Bp);
A_e=eye(n1+m1,n1+m1);
A_e(1:n1,1:n1)=Ap;
A_e(n1+1:n1+m1,1:n1)=Cp*Ap;
B_e=zeros(n1+m1,n_in);
B_e(1:n1,:)=Bp;
B_e(n1+1:n1+m1,:)=Cp*Bp;
C_e=zeros(m1,n1+m1);
C_e(:,n1+1:n1+m1)=eye(m1,m1);n=n1+m1;
h(1,:)=C_e;
F(1,:)=C_e*A_e;
for kk=2:Np
h(kk,:)=h(kk-1,:)*A_e;
F(kk,:)= F(kk-1,:)*A_e;
end
v=h*B_e;
Phi=zeros(Np,Nc); %declare the dimension of Phi
Phi(:,1)=v; % first column of Phi
for i=2:Nc
Phi(:,i)=[zeros(i-1,1);v(1:Np-i+1,1)]; %Toeplitz matrix
end
BarRs=ones(Np,1);
Phi_Phi= Phi'*Phi;
Phi_F= Phi'*F;
Phi_R=Phi'*BarRs;
end

QPhild

自定义优化函数，也可以使用matlab自带的优化函数：quadprog

function eta=QPhild(H,f,A_cons,b)
% E=H;
% F=f;
% M=A_cons;
% gamma=b;
% eta =x
[n1,m1]=size(A_cons);
eta=-H\f;
kk=0;
for i=1:n1
if (A_cons(i,:)*eta>b(i))kk=kk+1;
else
kk=kk+0;
end
end
if (kk==0) return;
end
P=A_cons*(H\A_cons');
d=(A_cons*(H\f)+b);
[n,m]=size(d);
x_ini=zeros(n,m);
lambda=x_ini;
al=10;
for km=1:38
%find the elements in the solution vector one by one
% km could be larger if the Lagranger multiplier has a slow
% convergence rate.
lambda_p=lambda;
for i=1:n
w= P(i,:)*lambda-P(i,i)*lambda(i,1);
w=w+d(i,1);
la=-w/P(i,i);
lambda(i,1)=max(0,la);
end
al=(lambda-lambda_p)'*(lambda-lambda_p);
if (al<10e-8); break;
end
end
eta=-H\f -H\A_cons'*lambda;
end

输出

参考文献

Wang L. Model predictive control system design and implementation using MATLAB®[M]. Springer Science & Business Media, 2009.

8,MPC的简单matlab实现相关推荐

MPC入门与Matlab实现
本文为B站视频<你还在用PID?MPC模型预测控制,从公式到代码!>的学习笔记,强烈推荐去看这位大佬的视频,链接放在了最后,别忘了给大佬一键三连哈 MPC入门与Matlab实现前言 1. ...
MPC模型预测控制+matlab代码实现+simulink仿真实现
MPC模型预测控制+matlab代码实现 MPC模型预测控制原理步骤 matlab代码实现程序主程序函数参数调试 MPC模型预测控制原理步骤 matlab代码实现程序主程序 cle ...
模型预测控制(MPC)的简单实现 — Matlab
参考文章:原理推导基于离散模型,首先复述其首先原理,然后基于简单的无人机模型进行仿真. 原理问题描述考虑离散系统: xk+1=Axk+Buk{x_{k + 1}} = A{x_k} + B{u_ ...
F-OFDM 系统简单Matlab搭建
1.源论文陈勇. 面向5G的F-OFDM关键技术研究[D].电子科技大学,2020. 本文相关参数均使用了论文中的参数,使用AWGN信道,但为了简单起见调制方式使用BPSK. 2.代码文件(Matl ...
简单Matlab的Gui设计——电子琴
输入guide,拖动控件修改前景色和背景色(黑底白字) 将代码拷贝到对应按钮的回调函数 Fs = 44100; dt = 1.0/Fs; T = 1; N = T/dt; t = linspace( ...
大M法的简单matlab程序
题目如上: 将其化为: min f=-50x1-40x1 -x1-x2<=-50 大M法函数代码如下: function [ A,b,c ] = M ( a,B,C) %使用前提:化为min z ...
一个简单的MATLAB脚本——快速行进算法（FMM））
一个简单的MATLAB脚本--快速行进算法(FMM) 介绍快速行进算法(FMM)的简单MATLAB脚本,不到20行代码实现快速行进算法的运算结果,而且计算速度非常快.给了两个实例模型来说明计算结果. ...
学习随笔#15 MPC控制MATLAB代码详解
MPC控制的详细数学推导可以参照文章:一个MPC详细建模的例子.MATLAB/SIMULINK中自带有MPC相关的工具,但本文给出MPC控制的MATLAB程序. function [M, C, Q ...
Matlab: Adaptive MPC Design, Time-Varying MPC Design, Nonlinear MPC Design
Matlab里的自适应mpc.(线性)时变MPC.非线性MPC有什么联系和区别? 看看帮助文档里Adaptive MPC Design: Adaptive MPC controllers adjust ...

8,MPC的简单matlab实现

线性MPC