模型预测控制（MPC）简介及matlab实现
基本概念MPC算法 整体流程数学建模线性模型 MPC与PID的区别Matlab实现MPC函数
基本概念
全称：Model-based Predictive Control(MPC)—模型预测控制
本质：MPC利用一个已有的模型、系统当前的状态和未来的控制量，来预测系统未来的输出，然后与我们期望的系统输出做比较，得到代价函数，通过优化的方法，优化出未来控制量，使得代价函数最小。优化出来的控制量即算法的输出。
核心思想：以优化方法求解最优控制器，其中优化方法大多时候采用二次规划（Quadratic Programming）
控制输出：MPC控制器优化得到的控制输出也是系统在未来有限时间步的控制序列。 当然，由于理论构建的模型与系统真实模型都有误差，所以，实际上更远未来的控制输出对系统控制的价值很低，故MPC仅执行输出序列中的第一个控制输出。
模型分类：
机理模型基于数据的模型，例如用神经网络训练一个modle，使用基于数据的模型的MPC可以结合model based RL使用预测：根据已有模型，已经未来有限时间步的控制序列，预测出未来的系统状态
控制：优化的控制量，是一个未来有限时间步的控制序列
预测区间（Predictive Horizon）：
对于一般的离散化系统，在k时刻，我们可以测量出系统的当前状态y(k),再通过优化计算得到u(k),u(k+1),u(k+2),u(k+j),根据模型与控制量，得到系统未来状态的估计值y(k),y(k+1),y(k+2),y(k+j)。
其中系统未来状态的估计值y(k),y(k+1),y(k+2),y(k+j)这部分就称为预测区间的系统预测值，指的是一次优化后预测未来输出的时间步的个数
控制区间（Control Horizon）
将控制估计的部分称为控制区间，通过优化计算得到u(k),u(k+1),u(k+2),u(k+j)这部分就称为控制区间的控制量，
在得到最优输入之后，我们只施加当前时刻的输入u(k)，即控制区间的第一位控制输入
区间参数配置：
过小的控制区间，可能无法做到较好的控制，而较大的控制区间，比如与预测区间相等，则会导致只有前一部分的控制范围才会有较好的效果，而后一部分的控制范围则收效甚微，而且将带来大量的计算开销。
约束：对于约束，一般分为：
Hard约束：物理性质的约束，不可以违背，例如方向盘的转向，刹车的深度Soft约束：软件约束，可以违反，例如最大的速度MPC 优点：
善于处理多输入多输出系统（MIMO）可以处理约束，如安全性约束，上下阈值是一种向前考虑未来时间步的有限时域优化方法（一定的预测能力）控制框图：
MPC算法 整体流程
模型预测控制在k时刻共需三步
第一步：获取系统的当前状态；
第二步：基于u(k),u(k+1),u(k+2),u(k+j)进行最优化处理，代价函数为
其中EN 表示误差的终值，也是衡量优劣的一种标准。
第三步：只取u ( k )作为控制输入施加在系统上。
在下一时刻重复以上三步，在下一步进行预测时使用的就是下一步的状态值，我们将这样的方案称为滚动优化控制（Receding Horizon Control）。
数学建模
线性模型
当模型是线性的时候，MPC的设计求解一般使用二次规划方法。
设线性模型为以下形式：
假定未来m步的控制输入已知u(k),u(k+1),u(k+2),u(k+m),根据以上模型与输入，我们可以计算未来m步的状态：
将上面m步写成矩阵向量形式：
其中，
上式B中的下三角形式，直接反映了系统在时间上的因果关系，即k + 1时刻的输入对k 时刻的输出没有影响，k+2时刻的输入对k和k+1时刻没有影响。
假定参考轨迹为
则MPC的一个简单的目标代价函数如下：
uTRu 这一项是为了让控制输入不会太大，因此代价函数中添加了一项对控制量的约束。
将状态方程带入代价函数，变量仅剩u,以上最优化问题可用二次规划方法求解，得到满足目标代价函数的最优控制序列
MPC与PID的区别
PID控制器不具有“前瞻性”：参与计算的各个量，有当前的 ，上个控制周期的 ，以及之前所有的 累计和，没有未来的 。PID属于无模型控制。PID仅仅通过目标和当前状态的差距，以及三个控制参数，就输出控制量Matlab实现MPC函数
function [M,C,Q_bar,R_bar,G,E,H,U_k] = MPC_Zero_Ref(A,B,N,x_k,Q,R,F) %%%%%%%%%%%建立一个以0为参考目标的MPC求解函数%%%%%%%%%%%其中，状态矩阵A，输入矩阵B系统维度N，初始条件x_k,权重矩阵Q,R及终端误差矩阵F为输入%%%%%%%%%%%输出中U_k为所求控制器，其余为简化过程中引入的中间变量 n=size(A,1); %A是n×n矩阵，求np=size(B,2); %B是n×p矩阵，求pM=[eye(n);zeros(N*n,n)];%初始化M矩阵，M矩阵是(N+1)n × n的， %它上面是n × n个“I”，这一步先把下半部分写成0C=zeros((N+1)*n,N*p);%初始化C矩阵，这一步令它有(N+1)n × NP个0%定义M和Ctmp=eye(n);%定义一个n × n的I矩阵for i=1:N%循环，i从1到N rows = i*n+(1:n);%定义当前行数，从i×n开始，共n行 C(rows,:)=[tmp*B,C(rows-n,1:end-p)];%将C矩阵填满 tmp=A*tmp;%每一次将tmp左乘一次A M(rows,:)=tmp;%将M矩阵写满end%定义Q_barS_q=size(Q,1);%找到Q的维度S_r=size(R,1);%找到R的维度Q_bar=zeros((N+1)*S_q,(N+1)*S_q);%初始化Q_bar为全0矩阵for i=0:N Q_bar(i*S_q+1:(i+1)*S_q,i*S_q+1:(i+1)*S_q)=Q;%将Q写到Q_bar的对角线上endQ_bar(N*S_q+1:(N+1)*S_q,N*S_q+1:(N+1)*S_q)=F;%将F放在最后一个位置 %定义R_barR_bar=zeros(N*S_r,N*S_r);%初始化R_bar为全0矩阵for i=0:N-1 R_bar(i*S_r+1:(i+1)*S_r,i*S_r+1:(i+1)*S_r)=R;end %求解G=M'*Q_bar*M;%GE=C'*Q_bar*M;%EH=C'*Q_bar*C+R_bar;%H%最优化f=(x_k'*E')';%定义f矩阵U_k=quadprog(H,f);%用二次规划求解最优化U_kend
