模型预测控制（MPC）解析（一）：模型
摘要
本文探讨了基于预测控制（MPC）的方法及其在工业自动化中的应用。MPC通过建立系统的动态预测模型，在有限的时间窗口内优化未来的输入轨迹以实现目标输出。文章详细介绍了MPC的工作流程：从建立系统的数学表达式开始，通过状态空间模型描述系统行为，并将其转换为增广状态空间形式以便于后续计算；随后讨论了如何利用优化算法求解最优控制序列以最小化预测误差与控制量之间的平衡关系。此外，文章还展示了如何使用Matlab语言实现上述计算过程，并通过具体示例验证了算法的有效性与正确性。最终结果表明，在不同参数设置下（如调整权重矩阵r_w），系统能够快速收敛至目标值并保持稳定运行。
一、MPC简介
1.1 预测控制的日常应用
该模型的主要目的是确定未来时间段内各时间段内各变量的变化趋势和最优调节策略以实现系统的稳定运行与预期目标的有效达成这一过程通常会在预设的时间范围内展开并且依赖于初始时刻的状态信息作为决策依据这有助于实现对复杂动态系统的有效管理与实时响应
工作将在上午9时启动，在此期间团队将致力于开发一个液态容器模型预测控制系统的设计方案，并最终实现这一控制系统的功能。我们制定了详细的8小时工作安排表，在执行过程中，我们仅依据下一个小时的计划进行操作。在完成第一个小时的工作后，我们将重新审视并制定新的8小时工作计划。不断重复上述步骤直至项目完成。
在9点整时安排工作计划时
综合考虑了各项因素后，在接下来的8个小时中将任务转化为可调节变量的函数模型，并详细计算出完成任务所需的时间分配方案。在制定计划时需考虑到自身能力的限制，在有限的时间内寻找最优的工作策略。通过系统分析得出从上午9点到下午5点的工作安排后，在第一小时内实施该时间段的具体工作内容。
在上午十时整，在过去的一个小时内我们需要对所完成的工作情况进行一次全面评估。这些数据将被用来规划我们下一阶段的工作计划。然而，在实际操作中可能出现各种因素干扰可能导致这一小时内无法达成既定目标。到了十点钟整，在此之前我们已经完成了当前一小时的工作安排，并根据实际表现对下一阶段的工作进行了重新规划。最终目标可能会保持不变状态或者根据实际情况进行调整制定计划的时间段固定在8小时不会改变按照上述步骤重复执行下一个阶段的计划过程仅执行下一个时间段内的工作内容到了下午一点整我们会再次评估当前的表现并据此制定新的八小时工作计划每个时间段内的实施过程都将遵循相同的流程直到最终目标得以实现
制定计划时需要重点关注的三个要素包括预判可能发生的事件（即建立理论框架）、监测当前状态与进度（即应用动态反馈机制）、以及执行计划中的各项任务（即实现调控机制）。规划工作主要包括以下内容：
- 工作时间被固定为8小时；
- 在进行规划之前需要了解当前的工作状态；
- 基于约束条件分析和应用的方法，在8小时的工作周期内制定最优计划，并结合最新的可用数据，在动态的时间窗口内进行实时调整和优化。
所述内容中的规划行为与MPC的核心思想一致。在这一实例中，提出了若干关键概念，在后续讨论中会频繁引用这些概念。
- 移动时间窗口：从任意时刻
到
的时间窗口。窗口的长度
改写说明
，
是时间窗口的开始点，每小时进行递增，规划从上午9点开始的，所以开始时
。
- 预测时域：规定了希望未来能被预测到多久。这个参数与时间窗口长度
滚动滑动窗口控制：其最优轨迹完全由移动的时间窗口所限定，在规划阶段仅采用该时间段内的初始控制信号而不考虑后续时段的影响。
实时状态：在规划过程中，实时获取必要的信息来预测未来的动态变化。
这一信息被用来构建相应的控制系统模型。
这是一个包含相关因素的向量，在测量或预测中通常被使用。
- 模型：在预测控制领域中构建系统的动态模型至关重要。
- 为了实现最佳决策目标
其核心作用是量化期望结果与实际输出之间的差异程度
其数学表达式称为代价函数J
能够使得优化窗口内的代价函数值最小化的控制轨迹即为最优规划轨迹。
1.2 预测控制中的模型
该文本已经按照要求进行了同义改写以降低重复率：对原文进行了词汇替换、句式变换以及语序调整等处理方式，并且保持了原有段落结构和数学公式的准确性。改写后的文本更加简洁流畅，在不改变原意的基础上提升了可读性。
二、从连续时间模型到增广模型
2.1 连续模型到离散模型
连续时间状态空间模型可以表示为：
对应的离散时间状态空间模型可以表示为：
在MATLAB环境中，我们可以利用c2dm函数来设定相应的采样时间间隔Δt，并实现对连续系统建模到离散系统的转化。
已知连续时间的状态空间模型如下所示；请建立相应的离散时间状态空间模型，并采用均匀采样周期为1的时间间隔。
MATLAB程序如下：
Ac = [0 1 0; 3 0 1; 0 1 0 ];
Bc= [1; 1; 3];
Cc=[0 1 0];
Dc=zeros(1,1);
Delta_t=1;
[Ad,Bd,Cd,Dd]=c2dm(Ac,Bc,Cc,Dc,Delta_t);得到的离散模型矩阵如下：
2.2 离散模型到增广模型
增广模型如下：
其中，
增广模型的矩阵分别为：
例： 离散时间状态空间模型的系统矩阵如下，求对应的增广模型矩阵(A,B,C)。
首先要确定
的列数，由增广模型可知，
的列数与
的行数相同，那么
， 将
代入增广模型矩阵方程中，就可以得到增广模型矩阵：
用MATLAB实现上述过程如下：
Ad = [1 1;0 1];
Bd = [0.5;1];
Cd = [1 0];
[m1,n1]=size(Cd);
[n1,n_in]=size(Bd);
A_e=eye(n1+m1,n1+m1);
A_e(1:n1,1:n1)=Ad;
A_e(n1+1:n1+m1,1:n1)=Cd*Ad;
B_e=zeros(n1+m1,n_in);
B_e(1:n1,:)=Bd;
B_e(n1+1:n1+m1,:)=Cd*Bd;
C_e=zeros(m1,n1+m1);
C_e(:,n1+1:n1+m1)=eye(m1,m1);得到的结果如下：
与直接代入计算的结果完全相同，验证了MATLAB程序的正确性。
三、状态预测与优化
3.1 状态和输出变量的预测
假设当前时刻为
，
为通过检测得到的系统状态变量，这里假设状态变量可以检测到。那么
个未来的控制变量序列表示为：
（
不是表示当前的时刻吗，为什么
它代表了后续时刻的控制量？由于当前时刻的控制量在经过优化后才能确定。
也是未来的控制量）
个预测的状态变量表示为：
代入增广模型中，得到
个预测的状态变量为：
那么，
个预测的输出变量为：
可以看出，所有预测的变量都只与当前观测到的系统状态
以及未来的控制变量序列有关：
将未来的输出变量序列和未来的控制变量序列写成：
在单输入，单输出系统中，Y的维数是
，
的维数是
。将上式合起来写得到：
其中，
3.2 优化
在
时刻的期望输入值为
，需通过MPC实现系统的控制目标；希望使系统的输出尽可能贴近系统的输入。假设在优化窗口期间, 系统的输入保持恒定。问题就转化为寻求能够在优化窗口中确定出一组最优控制变量序列; 以使系统的输出与输入之间的误差达到最小化 。
将输入
改写成优化窗口内的输入序列：
序列长度为
，定义代价函数为：
代价函数的第一项表征了输入信号与预期输出之间的误差；而第二项则旨在约束控制变量。
，避免控制量太大。其中
是对角矩阵：
其中，
反映了闭环性能的调节变量。在闭环系统中进行设置时。
，则说明不用考虑控制变量
有多大，只考虑控制误差越小越好。如果
取值比较大，则表示需要仔细考虑
的大小。
将代价函数写成：
要求上述代价函数的最小值，最简单粗暴的方法就是求一阶导：
当一阶导为零时，可以求得控制变量的最优解：
其中，假设了矩阵
存在。将优化的
写成与
和
有关的形式为：
其中，
长度为
。至此，就将优化的控制序列求解出来了。
例： 一阶系统状态空间方程如下：
其中，
，求(1)增广状态空间模型；(2)计算矩阵
；(3)假设在时刻
，
，当
分别为0和10的时候，求解最优
，并比较结果。
解： 增广状态空间方程为：
矩阵
如下：
矩阵
中的系数计算方式如下：
其中，
矩阵
中的系数计算方式如下：
将
代入得到：
当
在该时刻时，在这种情况下（即控制系统只关注减少输入输出之间的偏差），并未考虑控制量的大小时，则会计算得到的结果
个控制量序列为：
可以看出控制量变化比较大。
当
时，即限制了控制量变化大小的情况下，计算得到的
个控制量序列为：
与之前一种情况相比，控制量平稳很多。
将这两种情况下的输出量和状态量曲线画出来如下图所示。可以看到
时状态变化量减少到0时，预测输出可以到达期望输入。而在
时，预测输出缓慢向期望输入靠近，在预测窗口内误差比较大，收敛速度比
慢，但是状态变化曲线更加平稳。
图1
和
时状态变化及输出量曲线
由对比结果可知，在调节过程中若希望实现更加平滑的波动，则必须经过较长的时间段。当调整控制域至特定位置时
增加到
，在
的时候得到的最优控制量序列为：
观察得知，在前四个数据中各不相同，并且均有所下降。根据不同情境，应采取适当的策略以应对不同的情况。
来进行MPC控制。
3.3 MATLAB实现
用MATLAB函数实现上述计算各矩阵的过程，最关键的是求矩阵
。其中
属于托普利兹矩阵类别的矩阵具有特定结构特征：其每一列（除了第一列）都可以通过将前一列向右移动一位来生成。构建一个数学模型用于描述离散时间的状态变化过程：该模型将被用作函数的输入参数。
以及预测时域和控制时域长度
。MATLAB程序如下：
function [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;为了验证程序的正确性，在MATLAB的work space输入：
Ap = 0.8;
Bp = 0.1;
Cp = 1;
Nc = 4;
Np = 10;
[Phi_Phi, Phi_F, Phi_R, A_e, B_e, C_e] = mpcgain(Ap, Bp, Cp, Nc, Np);运行后得到：
可以看出计算结果与上例中的相同，表明程序正确。
四**、**参考文献
[1] L. Wang. Model Predictive Control System: Design and Implementation, Implemented in MATLAB Software.