经典预测控制算法：动态矩阵控制(DMC)上篇——理论推导

前言

  在自动控制理论领域，基于对象的传递函数，使用劳斯判据、根轨迹、伯德图、奈奎斯特曲线等工具的古典自动控制原理，和基于对象的状态空间方程，使用李雅普诺夫定律、状态观测器等工具的现代控制理论（这里指的是狭义上的现代控制理论，有些比较广义的说法中，现代控制理论还把模糊控制、神经网络控制、预测控制等都包含进去了，个人倾向于狭义的说法）无疑是学术上的两座高峰，它们都有着坚固的数学基础，经过严格的数学公式推导与证明，从理论上保证了所设计的自动控制系统稳定（这里的稳定类似于数学上函数收敛的概念），应用十分广泛。但也有一些专家学者另辟蹊径，提出一些十分有趣的控制理论，这些控制理论的基础架构完全独立于前面两者，自成一派。可惜的是这些控制理论由于某些缺陷（文章的最后我们会稍稍讨论一下缺陷）导致应用范围十分受限，个人认为除非学术上需要，这些控制理论不必做过多研究，实际中应用也很少（工业慢对象上基本是PID一招鲜吃遍天），但是这些算法蕴含的思想还是值得一讲的。
  上个世纪，基于状态空间方程的现代控制理论在航空航天等高精度，需要快速控制的领域成就卓越，解决了很多古典控制理论中的难题，但工业被控对象有非线性，耦合性强，参数还会随时间变化等特点，无法建立精确模型，导致现代控制理论应用效果很差，于是人们开始探索一些更适用于工业被控对象的智能控制算法，预测控制就是在这样的背景下诞生的。预测控制算法是一类算法的统称，比较经典的预测控制算法有模型预测控制(MPC)，动态矩阵控制(Dynamic Matrix Control, DMC)，广义预测控制(GPC)，它们的区别是预测模型有所不同，但这些算法的思想都是一致的：预测模型——滚动优化——反馈校正，通过三个步骤的配合，预测控制在一些行业取得了良好的效果。
  预测控制为什么有效果呢？打个不恰当的比方，两个人洗澡，一个人调节水温用手试，发现太烫了往右拧，又太冷了再往左拧，就这样来回来回好几次，水温终于调好了。另一个人就比较厉害了，他将阀门和水温的关系建立了一个预测模型，提前进行预测，现在水太冷，他把阀门按照预期的水温直接就拧到了对应位置，实际水温可能会比预期的稍高或稍低，他会先修正预测模型，然后再微调阀门位置。第二个人的水温调节过程要比第一个人优雅许多，但是他需要一个更加强大的大脑来构建预测模型。
  动态矩阵控制在经典的预测控制算法中相对好理解一些，下面我会分几篇文章来介绍。预测控制算法的共性就是矩阵运算很多，初次看推导会觉得很绕，但实际应用不需要把公式全部用代码实现，只有几个涉及到结论的公式需要编程，用着用着再回来看看推导，就会有一种豁然开朗的感觉。下面我会尽量讲的通俗一些，但还是需要一些控制理论和线性代数相关的基础，可以先看结论，想理解透彻再去翻翻相关内容。

算法推导

  首先要介绍一些基本概念，我们要控制的东西，专业名词叫被控对象（也称作系统/被控系统，后面我也会混用）。控制理论中的很多概念都和数学有关联，数学中的一个简单函数$y=f(x)$，自变量$x$变化会导致因变量$y$变化，到了控制理论中，函数的$x$和$y$有了具体含义，分别是被控对象的输入和被控对象的输出，函数有多个自变量或多个因变量，对应到控制理论中，这个系统就是多输入/多输出的（输入也叫控制量，输出也叫被控量）。$y=f(x)$如果是完全连续的，称作线性的，如果存在不连续的部分，称作非线性的，与之相似，被控对象同样也分为线性对象和非线性对象。
  早期的控制理论研究的被控对象都是线性的，虽然后来针对非线性对象有了一些研究成果，但只能局限于一些特殊情况或者得到一些特殊解。工程上很多被控对象都是非线性的，但经常采用的方法是不针对非线性本身进行研究，而将非线性对象分解为多个线性对象（非线性对象在局部/一段范围内一般也有很好的线性特征），这样我们就可以使用各种针对线性对象的工具了。

预测模型

  最基本的动态矩阵控制针对的是完全线性的对象，下面我们针对单输入单输出的线性对象进行推导。线性对象有一个特点，满足叠加性，也就是说给对象$y=f(u)$两个不同输入$u_1$和$u_2$，得到的输出$y_1$和$y_2$相加，与给对象输入$u_1+u_2$得到的输出是相等的，这个与数学中的函数特性是一致的：$$\begin{array}{c}{y}_1+y_2=f(u_1)+f(u_2)=f(u_1+u_2)\end{array}$$  式(1)非常简单，却是预测控制得以对未来进行预测的基础。控制系统稳定之后（输入不再变化），再给系统一个新的阶跃输入$\Delta u$，间隔相同采样周期$dt$采集$N+1$个输出（$N$称作建模时域），我们会得到一组输出的采样序列$y_0,y_1,\dots ,y_N$，如图所示：

  类似于深度学习需要先进行数据预处理一样，我们这里也对序列进行平移和标准化（标准化就是所有数据都除以$\Delta u$，也就是无论阶跃输入是多少，都统一成$\Delta u=1$）：
$$\begin{array}{c}\text{a}=\left[\begin{array}{c}{a}_1\\ a_2\\ \dots \\ a_N\end{array}\right]=\left[\begin{array}{c}\frac{y_1-y_0}{\Delta u}\\ \frac{y_2-y_0}{\Delta u}\\ \dots \\ \frac{y_N-y_0}{\Delta u}\end{array}\right]\end{array}$$  上式得到的$\text{a}$称作阶跃响应模型，$\text{a}$用一列数表示了给系统施加单位输入时输出的数值，这就是DMC的预测模型。DMC的预测模型很简单对吧，这个“模型”简单到连公式都没有，哈哈。利用$\text{a}$，就可以对接下来$P$个时刻的系统输出进行预测，$P$称作预测时域，显然$P$一定是小于等于$N$的，因为预测模型只知道到$N$为止输出的情况。补充一下，有些推导中是预测到$N$个未来时刻输出，个人认为这样会导致后面的公式$P$和$N$看的有点乱，我倾向于预测全部用$P$表示，毕竟它叫做“预测时域”。大家知道一下，我后面的$P$全部用$N$代替也是可以的，因为$P$必须小于等于$N$这个规则没有违背。
  有了预测模型，下面开始推导算法，假设现在所处的时刻是$k$时刻，且$k$时刻及之后输入$u$都不再变化（注意是不再变化，即输入的增量$\Delta u$为0，输入$u$本身可能是不为0的）时预测未来的输出为${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})$，${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})$是一组数不是一个数：
$$\begin{array}{c}{\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})=\left[\begin{array}{c}{\tilde{y}}_0(k+1|k)\\ {\tilde{y}}_0(k+2|k)\\ \dots \\ {\tilde{y}}_0(k+P|k)]\end{array}\right]\end{array}$$  后续都用加粗表示这个变量是一个矩阵（包括行/列向量），式(3)里面的变量都只是一种表示方式，一种指代的代号，后面出现类似形式的变量含义也是类似的。后面公式中所有带波浪号~变量表示是一个估计（预测）的值，${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})$中的下标$P$表示对未来进行$P$步预测，下标0表示未来没有施加控制增量，$k$表示现在是$k$时刻，${\tilde{y}}_0(k+1|k)$中的下标0表示未来没有施加控制增量，$k+1|k$表示在$k$时刻预测$k+1$时刻，$k+2|k$表示在$k$时刻预测$k+2$时刻，依此类推。
  你可能会奇怪还没有开始进行预测，怎么就冒出来一堆要预测的值呢？很好的问题！${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})$里面每个数是多少，我们确实还不知道，预测模型里只提供了输入变化时预测未来的输出值，并没有提供输入不变时对未来的预测。先别急，我们知道式(3)中的变量都还没有值，先当作已经知道好了，这是数学中常有的一种思路，虽然某个值未知，但我们可以装作它已知，然后通过一系列推导后，最后得到它的值。
  假设式(3)已知情况下，在$k$时刻施加一个控制增量$\Delta u(k)$，根据叠加性原理，显然可以预测出施加控制增量后的未来$P$个输出${\tilde{\mathbf{y}}}_{\mathbf{P}1}(\mathbf{k})$为：
$$\begin{array}{c}{\tilde{\mathbf{y}}}_{\mathbf{P}1}(\mathbf{k})={\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})+\text{a}\Delta u(k)\end{array}$$  式(4)是一种精简写法，详细列一下的话其实式(4)由$P$行公式构成：
$$\begin{array}{c}\{\begin{array}{l}{\tilde{y}}_1(k+1|k)={\tilde{y}}_0(k+1|k)+a_1\Delta u(k)\\ {\tilde{y}}_1(k+2|k)={\tilde{y}}_0(k+2|k)+a_2\Delta u(k)\\ \dots \\ {\tilde{y}}_1(k+P|k)={\tilde{y}}_0(k+P|k)+a_P\Delta u(k)\end{array}\end{array}$$  继续叠加，我们将一系列的控制增量一共$M$个施加给对象：
$$\begin{array}{c}\Delta {\mathbf{u}}_{\mathbf{M}}(\mathbf{k})=\left[\begin{array}{c}\Delta u(k)\\ \Delta u(k+1)\\ \dots \\ \Delta u(k+M-1)\end{array}\right]\end{array}$$  强调一下，DMC的运算都是一个时刻一个时刻推演下去的，所以步序非常重要！！！一定要清楚预测$k$时刻还是$k+1$时刻输出，控制增量是$k$时刻还是$k+1$时刻施加给对象的。上面的控制增量在$k$时刻施加给对象的是$\Delta u(k)$，在下个时刻也就是$k+1$时刻观测到输出的变化，$k+1$时刻施加给对象的是$\Delta u(k)+\Delta u(k+1)$，在下一时刻也就是$k+2$时刻观测到输出的变化，一直到$k+M-1$时刻把上面$M$个控制增量全部施加给对象，在下一时刻$k+M$时刻观测到输出的变化，在$k+M$时刻一直到$k+P-1$时刻控制增量不再继续叠加了。
  分开列公式，在$k$时刻及之后依次施加$M$个控制增量，基于$k$时刻预测$k+1$时刻到$k+M$时刻的输出为：
$$\begin{array}{c}\{\begin{array}{l}{\tilde{y}}_M(k+1|k)={\tilde{y}}_0(k+1|k)+a_1\Delta u(k)\\ {\tilde{y}}_M(k+2|k)={\tilde{y}}_0(k+2|k)+a_2\Delta u(k)+a_1\Delta u(k+1)\\ \dots \\ {\tilde{y}}_M(k+M|k)={\tilde{y}}_0(k+M|k)+a_M\Delta u(k)+a_{M-1}\Delta u(k-1)+\cdots +a_1\Delta u(k+M-1)\end{array}\end{array}$$  在$k$时刻及之后依次施加$M$个控制增量，基于$k$时刻预测$k+M+1$时刻到$k+P$时刻的输出为：
$$\begin{array}{c}\{\begin{array}{l}{\tilde{y}}_M(k+M+1|k)={\tilde{y}}_0(k+M+1|k)+a_{M+1}\Delta u(k)+a_M\Delta u(k+1)+\cdots +a_2\Delta u(k+M-1)\\ {\tilde{y}}_M(k+M+2|k)={\tilde{y}}_0(k+M+2|k)+a_{M+2}\Delta u(k)+a_{M+1}\Delta u(k+1)+\cdots +a_3\Delta u(k+M-1)\\ \dots \\ {\tilde{y}}_M(k+P|k)={\tilde{y}}_0(k+P|k)+a_P\Delta u(k)+a_{P-1}\Delta u(k+1)+\cdots +a_{P-M+1}\Delta u(k+M-1)\end{array}\end{array}$$  上面的式(7)和式(8)写成矩阵形式：
$$\begin{array}{c}{\tilde{\mathbf{y}}}_{\mathbf{PM}}(\mathbf{k})={\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})+\mathbf{A}\Delta {\mathbf{u}}_{\mathbf{M}}(\mathbf{k})\end{array}$$  其中：
$$\begin{array}{c}{\tilde{\mathbf{y}}}_{\mathbf{PM}}(\mathbf{k})=\left[\begin{array}{c}{\tilde{y}}_M(k+1|k)\\ {\tilde{y}}_M(k+2|k)\\ \dots \\ {\tilde{y}}_M(k+P|k)]\end{array}\right]\end{array}$$
$$\begin{array}{c}\mathbf{A}={\left[\begin{array}{ccc}{a}_1& \dots & 0\\ ⋮& \ddots & ⋮\\ a_M& \dots & a_1\\ ⋮& \ddots & ⋮\\ a_P& \dots & a_{P-M+1}\end{array}\right]}_{P\times M}\end{array}$$  这里的$\mathbf{A}$叫做动态矩阵，这就是动态矩阵控制的由来，可以看到动态矩阵其实就是预测模型取了前面一部分，重新排列组合构成的。
  前面已经假设过${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})$已知，那么在式(9)基础上，未来给系统输入任何已知控制增量，我们都可以预测出系统输出。到此为止，我们有了预测模型，并在预测模型基础上预测了未来有限个时刻内（$P$个时刻）对象的输出。

滚动优化

  所谓自动控制算法，目的是什么呢？是求得控制量的解，某个范围内，或者某些条件下，这个解对于被控对象来说是最优的。
  最优代表了什么呢？最优代表了我把求得的控制量输入给系统，能够使系统最快且最准地从当前状态移动到我期望的状态，比如我想让水温从现在的20℃升到30℃，要足够快，并且温度上升的过程中要精确，如果温度到达期望状态是这样的过程：从20℃升到35℃，又降到25℃，又升到33℃，又降到27℃，最后才到30℃，这样显然不准。
  在理论研究中，最优确实只需要关注怎样让输出又快又准，但是在实际工程中，最优还包括了控制量的稳定，因为实际工程中，控制量要交给执行机构去动作，比如阀门，不能一下开到顶，又一下关到死，这样大开大合，对实际设备损害极大。比如某一种算法，让水流量从0t/h升到100t/h，只需要一秒钟，一步到位，又快又准，但是观察控制量，从0%一秒到100%，下一秒又一下降到1%，然后下一秒50%，反复的开关开关，很伤设备。
  结合上面的讨论，先给出DMC的性能指标函数：
$$\begin{array}{c}\min (J)=\sum _{i=1}^P{q}_i[\omega (k+i)-{\tilde{y}}_M(k+i|k){]}^2+\sum _{j=1}^M{\lambda }_j\Delta u^2(k+j-1)\end{array}$$  要求出最优的控制增量，就要使得式(6)这$M$个控制增量带入式(12)后，在所有的控制增量中能使性能指标函数最小。式(12)出现了很多新的变量，下面慢慢来解释。
  $\omega (k+i)$叫做参考轨迹，表示的是$P$个公式：
$$\begin{array}{c}\omega (k+i)={\alpha }^{i}y(k)+(1-{\alpha }^i)y_{\text{SP}}(k+i),i=1,2,\cdots ,P\end{array}$$  $y_{\text{SP}}(k+i)$是从$k+1$时刻到$k+P$时刻每个时刻设定值，设定值是提前给定好的。$\alpha$叫做柔化因子，是一个0到1之间人工调节的参数，不同取值柔化强度不同。式(13)就是对设定值做一个“削弱”，传统的控制算法会把设定值直接提供给控制算法去计算，DMC是将柔化轨迹作为计算用的设定值，显然通过柔化轨迹这样慢慢增加设定值，控制过程会平缓一些，但速度上会减慢。在控制中，快速和稳定往往只能相对去平衡，鱼和熊掌不可兼得。
  ${\tilde{y}}_M(k+i|k)$就是式(10)中的预测值，$\Delta u(k+j-1)$就是式(6)这些控制增量，这是要求解的量，$q_i(i=1,2,\cdots ,P)$和${\lambda }_j(j=1,2,\cdots ,M)$是人工给定的系数，这些系数要一个一个给定太繁琐了，我们一般采用同一个系数，也就是$q_1=q_2=\cdots =q_P=q$和${\lambda }_1={\lambda }_2=\cdots ={\lambda }_M=\lambda$。
  式(12)是DMC最核心的公式，分为两部分，前一部分是未来预测值和设定值（参考轨迹）的误差乘以权重的累加和，后一部分是控制增量乘以权重的累加和，前一部分小，说明足够快和准，后一部分小，说明控制增量足够稳，结合起来，所给出的性能指标函数同时考虑了快准稳。
  还需要说明一下，式(12)中最能体现预测控制独特性的就是${\tilde{y}}_M(k+i|k)$这一项，一般的控制算法计算误差，都是采用当前的反馈值和当前的设定值，基于当前及历史上的已知信息求解最优控制量，而预测控制基于预测模型，预测了不同控制增量下未来的输出，从而可以在当前时刻将未来的控制增量全部规划好。两个人开车，一个人只看眼前的车流量选择最快的路线，而另一个预测到了路途上每个路口的车流量，提前规划好了全局最快速的行驶路线，显然后者更优秀。
  式(12)稍微加工一下，写成更美观的式(14)形式，这两个式子含义是一样的，只是式(14)的矩阵形式在数学上更加美观，符号是一一对应的不再解释了：
$$\begin{array}{c}\min (J)=\Vert {\mathbf{\omega }}_{\mathbf{P}}(\mathbf{k})-{\tilde{\mathbf{y}}}_{\mathbf{PM}}(\mathbf{k}){\Vert }_Q^2+\Vert \Delta {\mathbf{u}}_{\mathbf{M}}(\mathbf{k}){\Vert }_{\Lambda }^2\end{array}$$  上式中$\Delta {\mathbf{u}}_{\mathbf{M}}(\mathbf{k})$是要求解的未知量，其他都已知，这其实是一个求一元二次函数极值的问题，如果是普通的一元二次函数，我们初中就学过，求它的导数为0的点即可得到极值。式(14)比较麻烦的是多个未知量构成了矩阵，需要矩阵求导相关知识，这里就不做介绍了，感兴趣可以找找相关文献，方法上其实是十分类似的，跳过一下，直接得到求解的结果：
$$\begin{array}{c}\Delta {\mathbf{u}}_{\mathbf{M}}(\mathbf{k})=({\mathbf{A}}^{\text{T}}\mathbf{QA}+\mathbf{\Lambda }{)}^{-1}{\mathbf{A}}^{\text{T}}\mathbf{Q}[{\mathbf{\omega }}_{\mathbf{P}}(\mathbf{k})-{\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})]\end{array}$$  其中$\mathbf{Q}$和$\mathbf{\Lambda }$是$P\times P$和$M\times M$大小的矩阵，只有对角线为$q$和$\lambda$：
$$\begin{array}{c}\mathbf{Q}={\left[\begin{array}{ccc}q& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & q\end{array}\right]}_{P\times P}\end{array}$$ $$\begin{array}{c}\mathbf{\Lambda }={\left[\begin{array}{ccc}\lambda & \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & \lambda \end{array}\right]}_{M\times M}\end{array}$$  式(15)右边都是已知量，于是在$k$时刻我们就可以得到未来的$M$步控制增量，接下来在现有控制量上依次叠加控制增量即可，所以$k+1$时刻控制量为：
$$\begin{array}{c}u(k+1)=u(k)+\Delta u(k)\end{array}$$  虽然一次得到了未来的$M$步控制增量，但是如果到了$k+1$时刻会重新计算一遍控制增量，也就是说实际上每个时刻，我们预测了未来$P$步的输出，计算出了$M$步的控制增量，但实际上只采用其中第一个控制增量：
$$\begin{array}{c}\Delta u(k)=c^{\text{T}}\Delta {\mathbf{u}}_{\mathbf{M}}(\mathbf{k})={\mathbf{c}}^{\text{T}}({\mathbf{A}}^{\text{T}}\mathbf{QA}+\mathbf{\Lambda }{)}^{-1}{\mathbf{A}}^{\text{T}}\mathbf{Q}[{\mathbf{\omega }}_{\mathbf{P}}(\mathbf{k})-{\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})]\end{array}$$  这里的${\mathbf{c}}^{\text{T}}={\left[\begin{array}{cccc}1& 0& \dots & 0\end{array}\right]}_{1\times M}$。
  有些控制算法寻找最优控制量，是尽可能达到“全局最优”，而从前面的求解可以看出，预测控制求的最优解是限定在$P$个时刻内的，也就是说我们求的是局部最优解，随着时间推移，求解的范围也会跟着移动，但都是在基于“当前”的未来$P$步之内，不断向前，循环求解，这就是预测控制中滚动优化的思想。

反馈校正

  到这里为止，我们得到了DMC的控制律式(19)，那么建立预测模型，设定好参数，DMC算法就可以运行起来了吗？不，别忘记我们一开始提出的假设${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})$还是个未知量呢！前面为了得到控制律先假设它存在，现在该知道它到底是多少了，这就是DMC最后一个核心环节，反馈校正的工作。
  我们在$k$时刻将$\Delta u(k)$施加给对象，同时计算出未来的预测输出${\tilde{\mathbf{y}}}_{\mathbf{P}1}(\mathbf{k})$，注意因为只施加了一个控制增量，后面的控制增量还都没有施加也就是都为0，所以预测输出应该是${\tilde{\mathbf{y}}}_{\mathbf{P}1}(\mathbf{k})$，而不是${\tilde{\mathbf{y}}}_{\mathbf{PM}}(\mathbf{k})$。时间来到$k+1$时刻，现在可以观测到对象的实际输出$y(k+1)$，取出${\tilde{\mathbf{y}}}_{\mathbf{P}1}(\mathbf{k})$中的第一行${\tilde{y}}_1(k+1|k)$，求一下误差$e(k+1)$：
$$\begin{array}{c}e(k+1)=y(k+1)-{\tilde{y}}_1(k+1|k)\end{array}$$  这个误差表明$k+1$时刻预测模型得到的预测输出与实际输出的误差，这个误差不为0表示我们建立的预测模型不是完全准确，需要校正预测模型，校正的方法很简单，就是在${\tilde{\mathbf{y}}}_{\mathbf{P}1}(\mathbf{k})$基础上加上误差，为了灵活调整，这里有一个参数$\mathbf{h}$用来调节校正强度，和前面的$q$，$\lambda$一样，本身$\mathbf{h}$是$P$行1列的可调整参数，但是太繁琐了，同样令$h_1=h_2=\cdots =h_P=h$，一般没有特殊需要$h=1$即可，先分开列一下公式：
$$\begin{array}{c}\{\begin{array}{l}{\tilde{y}}_{COR}(k+1|k+1)={\tilde{y}}_1(k+1|k)+he(k+1)\\ {\tilde{y}}_{COR}(k+2|k+1)={\tilde{y}}_1(k+2|k)+he(k+1)\\ \dots \\ {\tilde{y}}_{COR}(k+P|k+1)={\tilde{y}}_1(k+P|k)+he(k+1)\end{array}\end{array}$$  然后同样用矩阵形式表达：
$$\begin{array}{c}{\tilde{\mathbf{y}}}_{\text{COR}}(\mathbf{k}+1)={\tilde{\mathbf{y}}}_{\mathbf{P}1}(\mathbf{k})+\text{h}e(k+1)\end{array}$$  ${\tilde{\mathbf{y}}}_{\text{COR}}(\mathbf{k}+1)$就是$k+1$时刻校正后的未来预测输出：
$$\begin{array}{c}{\tilde{\mathbf{y}}}_{\text{COR}}(\mathbf{k}+1)=\left[\begin{array}{c}{\tilde{y}}_{COR}(k+1|k+1)\\ {\tilde{y}}_{COR}(k+2|k+1)\\ \dots \\ {\tilde{y}}_{COR}(k+P|k+1)]\end{array}\right]\end{array}$$  上面的过程求实际输出和预测输出的误差，然后用来修正预测模型，这就是反馈校正的思想。
  得到$k+1$时刻校正后的未来预测输出，我们发现，此时未来的控制增量就是都为0的！也就是说，式(21)除了第一行，从第二行到最后一行，表示的就是在$k+1$时刻，未来控制增量全为0时对$k+2$到$k+P$时刻的输出预测！
$$\begin{array}{c}{\tilde{y}}_0(k+1+i|k+1)={\tilde{y}}_{\text{COR}}(k+1+i|k+1),i=1,2,\cdots ,P-1\end{array}$$  但是上面的${\tilde{y}}_0(k+1+i|k+1)$就等于${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k}+1)$吗?不是的，因为${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})$表示$k$时刻对$k+1$到$k+P$时刻的预测，那么${\tilde{y}}_0(k+1+i|k+1)$显然还得再加一行对$k+P+1$时刻的预测才能等价于${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k}+1)$，但是式(21)只计算到了$k+P$时刻，并没有得到$k+P+1$时刻预测输出。为了保证预测能持续进行下去，DMC简单粗暴地采用了近似，认为${\tilde{y}}_0(k+P+1|k+1)={\tilde{y}}_0(k+P|k+1)$，于是得到了完整的${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k}+1)$：
$$\begin{array}{c}{\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k}+1)=\mathbf{S}\cdot {\tilde{\mathbf{y}}}_{\text{COR}}(\mathbf{k}+1)\end{array}$$  其中$\mathbf{S}$叫做移位矩阵：
$$\begin{array}{c}\mathbf{S}={\left[\begin{array}{cccccc}0& 1& 0& \cdots & \cdots & 0\\ 0& 0& \ddots & \ddots & & ⋮\\ ⋮& ⋮& & \ddots & \ddots & 1\\ 0& 0& 0& \cdots & \cdots & 1\end{array}\right]}_{P\times P}\end{array}$$  这样$k+1$时刻就可以由基础预测值继续求解控制律，完成下一时刻的控制器输出和误差校正。其实到这里细心的你会发现，我们虽然得到了${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k}+1)$，但是${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})$还是没有嘛！确实，因为因果关系是先有$k$时刻的预测值，才能有$k+1$时刻的预测值，一个人没有迈出第一步，怎么迈出第二步第三步呢?为了在因果上能实现，我们回到起点，也就是0时刻，或者叫DMC计算的第一步，这时让${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})$全部为0，当然这个预测是和对象完全不符的，但是因为DMC有反馈校正环节，通过反馈校正，我们的${\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})$会逐渐逼近对象实际输出。自此，DMC就完整实现了闭环，可以编程来解决实际问题了。
  实际的编程并不需要把上面推导的公式全部实现，只需要能够求出控制增量即可，也就是式(19)右边全部要先求出来，其中${\mathbf{c}}^{\text{T}}({\mathbf{A}}^{\text{T}}\mathbf{QA}+\mathbf{\Lambda }{)}^{-1}{\mathbf{A}}^{\text{T}}\mathbf{Q}$在得到预测模型并且确定参数后就可以固定，不需要每次运算循环都计算一遍（除非参数变动或者模型有变化），只有${\mathbf{\omega }}_{\mathbf{P}}(\mathbf{k})-{\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})$这部分要每次循环都进行求解。

相关讨论

  通过上面的推导，我们感受了预测控制算法的独特思想，预测模型——滚动优化——反馈校正，DMC作为一种经典的预测控制算法，得到了一些实际应用，但是在理论上很难站稳脚跟，因为在它之前的控制算法，都是在数学稳定的基础上进行算法设计，DMC并没有在理论上证明稳定（这也与DMC算法本身基于离散的而不是连续的过程有关，古典控制理论很多工具都依赖于连续对象，离散对象可用的工具很少），而是先提出了算法本身。DMC算法到现在已经提出五六十年了，理论上的稳定性证明依然停留在一些特殊情况，在通用情况下的稳定性并没有得到证明，也就是说，虽然我们可以使用DMC算法做一些控制，但是前面的理论推导算是“想当然”的结果，可能预测着预测着模型与实际偏差越来越大，这都是说不准的。随着人工智能的发展，更先进的控制算法不断提出，DMC本身也一直没有得到继续发展——它的数学基础太薄弱了，但是算法蕴含的思想依然让人眼前一亮，希望大家看过之后能有一些启发。