经典预测控制算法：动态矩阵控制(DMC)中篇——算法改进

前言

  在上篇中介绍到，预测控制有预测模型——滚动优化——反馈校正三个部分，对应的每个部分都可以进行改进，于是基本DMC之上又衍生出了一些改进版本的DMC算法。另外还有一些从其他思路出发进行的算法改进，这里也简单做个介绍。

基于神经网络的预测控制

  上篇中提到，根据预测模型的不同，预测控制分为了不同算法，也就是说，将阶跃响应模型换成其他模型后，确实是进行了算法改进，但应该称之为新的预测控制算法，而不是改进的DMC了。由于是新的算法，这里就不再详细展开了，简单介绍一下。
  DMC的预测模型很简单，按照数学建模中自变量因变量的说法，自变量固定为1，因变量为$N$个预测输出，模型明显太脆弱了，实际的对象稍有误差就会导致模型不准。如果模型能够将历史上的输出和控制量都作为自变量进行未来输出预测，想必会准确很多。广义预测控制的预测模型就蕴含了这个思想，但是它在所使用的每个历史输出和控制量前都设置了一个系数，在求解控制律前，还需要先求解丢番图方程确定这些系数，如果模型使用的历史输出和控制量很多，求解会变得复杂。
  随着神经网络的发展，将神经网络模型和预测控制结合，构成神经网络预测控制，是一个比较热门的方向，神经网络本身有相当多的研究就是围绕数据预测展开的，直接把神经网络搬过来作为预测模型，添加滚动优化和反馈校正环节，自然就可以形成一种新的预测控制。建立神经网络预测模型的方法很明确，使用某时刻之前控制量及系统输出作为神经网络输入，未来系统输出作为神经网络输出，训练模型即可，甚至有一些相关的变量也可以作为神经网络输入。
  神经网络预测控制和广义预测控制最大的变化是前者求出的是数值解，后者可以求出解析解。因为广义预测控制可以使用和DMC类似的矩阵求导使得导数为0得到控制增量，而神经网络预测控制的预测模型是无法直观地表示为一个数学公式的，这就导致无法通过矩阵求导的方式来得到最优控制增量。为了求出最优控制增量，神经网络预测控制还需要加入寻优算法，先规定控制增量的最大最小范围，在该范围内通过穷举，或者黄金分割法等寻优方法来找到最优控制增量，也就是说，要通过代入不同控制增量到性能指标函数，反复比较得到最优控制增量。
  可能有人不太清楚解析解和数值解，我再换个说法表达一下，解析解就是如同DMC上篇中式(15)一样，按照式子经过加减乘除等数学运算就可以得到结果，而数值解的计算过程往往需要if else等逻辑判断才能实现，不能简单地用很少量的数学公式表达。

阶梯式动态矩阵控制

  在上篇中，可以看到DMC的控制律存在矩阵求逆环节，虽然在个人计算机上对几百乘几百的矩阵求逆很快，但工业现场的控制器配置较低（别忘记DMC刚刚提出时还在上个世纪，如果是当时，连普通计算机进行矩阵求逆也很耗时间），有的控制器甚至只是一个单片机，如果矩阵求逆太慢，导致没有及时算出控制量，控制出现滞后，会对整个系统的稳定调节造成恶劣影响。
  PID只需要几步加减乘除做代数运算，需要的算力很小，而这些智能控制算法需要更高的算力，这也是智能控制算法在实际中应用受限的一大原因。
  有些智能控制算法整体运算量大，比如涉及到优化算法，需要进行大规模的计算且无法拆分，减小运算量很难实现。而DMC的运算需求就集中在矩阵求逆这一个环节上，如果可以不进行矩阵求逆运算，DMC运算量就可以大幅减少。有学者通过引入阶梯式因子$\beta$，得到了阶梯式动态矩阵控制(Stair-like Dynamic Matrix Control, SDMC)算法，下面我们介绍介绍这种改进算法。回顾一下基本的DMC算法，它的控制律为(1)式：
$$\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]=({\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}$$  SDMC将需要求解的$M$个控制增量表示为式(2)形式，$\beta$人工给定（$\beta >0$），那么控制律求最优的$M$个未知数构成向量，就简化为了求单个未知数δ。
$$\begin{array}{c}\Delta u(k+j-1)={\beta }^{j-1}\delta ,j=1,2,\cdots ,M\end{array}$$  显然$0<\beta <1$时控制增量会逐渐减小，$\beta >1$会像指数函数一样逐渐增大。由于每次运算我们只采用计算出控制增量的第一步，$\beta >1$会导致第一步控制增量一直很小，这必然会使系统稳定的过程拉长，所以讨论的较少，有些文章中干脆限定了$0<\beta <1$。这里大家明白一下，$\beta >1$时算法是成立的，只是从调节过程“快速”的角度来说，$\beta >1$可以想象到会让控制变慢。下面画一下$0<\beta <1$时控制增量示意图，可以看到控制增量像阶梯一样，这便是阶梯式动态矩阵控制的名称由来。

  由预测模型得到施加阶梯式的$M$个控制增量后，未来输出的预测值为：
$$\begin{array}{c}{\tilde{\mathbf{y}}}_{\mathbf{PM}}\left(\mathbf{k}\right)={\tilde{\mathbf{y}}}_{\mathbf{P}0}(\mathbf{k})+\mathbf{A}\Delta {\mathbf{u}}_{\mathbf{M}}(\mathbf{k})={\tilde{\mathbf{y}}}_{\mathbf{P}0}\left(\mathbf{k}\right)+\left[\begin{array}{c}\begin{array}{c}{a}_1\\ a_2+a_1\beta \end{array}\\ \begin{array}{c}\cdots \\ a_P+a_{P-1}\beta +\cdots +a_{P-M+1}{\beta }^{M-1}\end{array}\end{array}\right]\delta \end{array}$$  于是最优目标函数也变为了下面的形式：
$$\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=\Vert {\mathbf{\omega }}_{\mathbf{P}}(\mathbf{k})-{\tilde{\mathbf{y}}}_{\mathbf{PM}}(\mathbf{k}){\Vert }_Q^2+\lambda (1+{\beta }^2+{\beta }^4+...+{\beta }^{2(M-1)}){\delta }^2\end{array}$$  同样进行求导，导数为0就可以得到δ的解：
$$\begin{array}{c}\Delta u(k)=\delta =\frac{{\mathbf{G}}^{\text{T}}\mathbf{Q}({\mathbf{\omega }}_{\mathbf{P}}\left(\mathbf{k}\right)-{\tilde{\mathbf{y}}}_{\mathbf{P}0}\left(\mathbf{k}\right))}{{\mathbf{G}}^{\text{T}}\mathbf{QG}+\lambda (1+{\beta }^2+\dots +{\beta }^{2(M-1)})}\end{array}$$  这里的$\mathbf{G}$是为了公式表达简练进行的简写：
$$\begin{array}{c}\mathbf{G}=\left[\begin{array}{c}{a}_1\\ a_2+a_1\beta \\ \cdots \\ a_P+a_{P-1}\beta +\cdots +a_{P-M+1}{\beta }^{M-1}\end{array}\right]\end{array}$$  其实$\mathbf{G}$在式(3)中已经出现了。
  原始DMC的控制律再搬出来对比一下，见式(7)，可以看到SDMC将矩阵求逆变为了普通的分式，能够大大加快运算速度（矩阵的加减法乘法和求逆完全不是一个数量级，比求逆要快得多）。
$$\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}$$  虽然DMC和SDMC实际使用的时候每次只需要第一步$\Delta u(k)$，但是我们还是把完整的$\Delta {\mathbf{u}}_{\mathbf{M}}(\mathbf{k})$列出来对比一下。
  DMC的$\Delta {\mathbf{u}}_{\mathbf{M}}(\mathbf{k})$如下：
$$\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}$$  SDMC的$\Delta {\mathbf{u}}_{\mathbf{M}}(\mathbf{k})$如下：
$$\begin{array}{c}\Delta {\mathbf{u}}_{\mathbf{M}}(\mathbf{k})=\left[\begin{array}{c}\delta \\ \beta \delta \\ \cdots \\ {\beta }^{M-1}\delta \end{array}\right]=\left[\begin{array}{c}\frac{{\mathbf{G}}^{\text{T}}\mathbf{Q}({\mathbf{\omega }}_{\mathbf{P}}\left(\mathbf{k}\right)-{\tilde{\mathbf{y}}}_{\mathbf{P}0}\left(\mathbf{k}\right))}{{\mathbf{G}}^{\text{T}}\mathbf{QG}+\lambda (1+{\beta }^2+\dots +{\beta }^{2(M-1)})}\\ \frac{\beta {\mathbf{G}}^{\text{T}}\mathbf{Q}({\mathbf{\omega }}_{\mathbf{P}}\left(\mathbf{k}\right)-{\tilde{\mathbf{y}}}_{\mathbf{P}0}\left(\mathbf{k}\right))}{{\mathbf{G}}^{\text{T}}\mathbf{QG}+\lambda (1+{\beta }^2+\dots +{\beta }^{2(M-1)})}\\ \cdots \\ \frac{{\beta }^{M-1}{\mathbf{G}}^{\text{T}}\mathbf{Q}({\mathbf{\omega }}_{\mathbf{P}}\left(\mathbf{k}\right)-{\tilde{\mathbf{y}}}_{\mathbf{P}0}\left(\mathbf{k}\right))}{{\mathbf{G}}^{\text{T}}\mathbf{QG}+\lambda (1+{\beta }^2+\dots +{\beta }^{2(M-1)})}\end{array}\right]\end{array}$$  DMC必须一次将$M$个控制增量算出，然后才能取第一行作为$\Delta u(k)$，而SDMC直接计算出了$\Delta u(k)$，并且将矩阵求逆简化为了算术除法，在滚动优化环节尝试进行了改进。

反馈校正环节截断误差改进

  回顾DMC的反馈校正部分，简单粗暴地采用了近似，认为${\tilde{y}}_0(k+P+1|k+1)={\tilde{y}}_0(k+P|k+1)$，使得算法能够成立，如果预测时域的长度覆盖了阶跃响应模型的整个动态变化过程时，这种近似不会产生误差。
  但是实际上，预测时域很可能不会这么长，比如一个惯性时间几百秒的对象，按一秒的步长采样，预测时域只设为一百多两百多，因为DMC大量矩阵规模都与$P$密切相关，预测时域的长度增加，很多矩阵规模都会增大，计算会明显变慢。
  那么当预测时域不好再继续增大，近似又会产生误差时，怎样弥补这个误差呢？其实就是对${\tilde{y}}_0(k+P+1|k+1)={\tilde{y}}_0(k+P|k+1)$这个式子进行改进，对$k+P+1$时刻的预测不采用$k+P$时刻近似，而是基于$k+P-1$时刻到$k+P$时刻的预测值变化趋势来推测$k+P+1$时刻的值。
  用一幅图来展示改进方案：
  图中$c$为真正的$y_{COR}\left(k+1+P|k+1\right)$，由于在$P$处预测模型截断无法得知其真实值，$d$是直接近似得到的${\tilde{y}}_{COR}\left(k+1+P|k+1\right)$，和$c$的误差为$c-d$。$e$点是$a，b$两点构成直线顺延得到的，$e=\left[\begin{array}{cc}-1& 2\end{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right]$，如果能够在$d，e$两点之间取足够接近$c$的$f$，必然能减小截断误差的累积：
$$\begin{array}{c}f=\left[\begin{array}{cc}-\xi & 1+\xi \end{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right](0<\xi <1)\end{array}$$  换一种表述方式，改进$k+P-1$时刻预测的近似值，其实就是将移位矩阵$\mathbf{S}$改进为了$\bar{\mathbf{S}}$：
$$\begin{array}{c}\bar{\mathbf{S}}=\left[\begin{array}{cccccc}0& 1& 0& \cdots & 0& 0\\ 0& 0& \ddots & \ddots & & ⋮\\ ⋮& ⋮& & \ddots & \ddots & 1\\ 0& 0& 0& \cdots & -\xi & 1+\xi \end{array}\right]\end{array}$$  DMC需要调节的参数太多了，这里又多了一个参数$\xi$，真是让人头大，还好已经有研究总结出了规律，按照经验设置即可：如果0 $<\frac{a_P-a_{P-1}}{a_{P-1}-a_{P-2}}<$ 0.95，则取$\xi =\frac{a_P-a_{P-1}}{a_{P-1}-a_{P-2}}$，0.95 $\le \frac{a_P-a_{P-1}}{a_{P-1}-a_{P-2}}\le$ 1取$\xi =$ 0.95~0.99均可，如果$\frac{a_P-a_{P-1}}{a_{P-1}-a_{P-2}}>$ 1说明在截断处输出仍在快速上升，按照$a_{P-1}$和$a_P$外推得到的趋势与到达新稳态的趋势方向相反，$\xi$并不能改善截断误差，此时$P$如果不能继续增大则只能使用移位矩阵$\mathbf{S}$。
  可能有人会疑惑，反馈校正不是用来消除预测模型误差的吗？那为什么要对近似产生的误差进行改进呢？我们再来明确一下，近似产生的误差和预测模型误差无关，是由于算法的“不合理”人为添加进来的误差。对于反馈校正来说，它是无法校正的随机误差，相当于每次校正时叠加了一个无法预知的随机数，而且由于算法一步步往前推进，这一误差会逐渐向当前时刻移动，同时新的截断误差进入反馈校正环节，预测模型会不断积累截断误差。以上就是对反馈校正环节尝试进行的改进。

其他改进思路

对象存在纯迟延时控制器参数优化

  一个存在$\tau$步迟延的对象，其$k$时刻计算出的控制增量只能影响未来$k+\tau$时刻之后的输出，不会影响到$k$到$k+\tau$输出，换言之$k$时刻仅能对$k$到$k+\tau$这一段时间的输出进行预测，但无法对输出进行控制，因此需要特别针对此类对象进行参数改进。
  假设去掉纯迟延环节后设计的DMC参数分别为$N,P,M,Q,\Lambda$，那么存在$\tau$步迟延时，参数修改为：$\bar{N}=N+\tau ,\bar{P}=P+\tau ,\bar{M}=M,\bar{Q}=\left[\begin{array}{cc}{0}_{\tau \times \tau }& 0_{\tau \times P}\\ 0_{P\times \tau }& Q\end{array}\right],\bar{\Lambda }=\Lambda$。
  如果没有纯迟延时阶跃响应模型为DMC上篇中的$\mathbf{a}$形式，存在纯迟延时阶跃响应模型$\bar{\mathbf{a}}$为下面形式：
$$\begin{array}{c}{\bar{a}}_i=\{\begin{array}{l}0,i\le \tau \\ a_i,\tau \le i\le \bar{N}\end{array}\end{array}$$

加入多模型思想构成多模型预测控制

  DMC算法是一种典型的依赖模型的控制算法，没有模型，那么后面的一切计算都是空谈，并且模型一定要足够准确。从这个方面来说，PID因为不依赖模型，任何场景都可以拿来即用，这也是它能够被广泛使用的原因。
  实际的工业过程中，只是建模准确还不能保证算法效果好，因为我们所建立的模型，只是在特定条件下准确。比如一个锅炉，不同温度下材料性能完全不同，导致了锅炉的整体模型会发生很大变化。
  那么这些基于模型的控制算法要怎么才能更适应实际呢？有学者提出多模型控制的思想，就是在不同条件下建立多个模型，比如400℃建立一个锅炉模型，500℃建立一个锅炉模型，600℃建立一个锅炉模型，凡是正常运行可能处在的状态都要尽可能建立模型。那么你自然而然就会想到，如果精确到每1℃，甚至每0.1℃都建立一个模型，那我的控制算法岂不是可以适用于任何情况？理论上当然如此，但是实际过程，建立一个模型可能就要几十分钟几小时，要建立成百上千个模型，锅炉设备都老化掉了！所以我们采取折中的办法，在比较典型的几个条件下建立模型，在这些条件之间类似于数学中的插值来得到模型，比如有200℃和300℃两个模型，那么250℃的结果用他俩加权来得到。这里的加权也分为两种，可以对模型的参数加权得到250℃的参数模型再计算输出，也可以先用200℃和300℃的模型得到各自输出，对输出加权，得到250℃输出。插值除了基本的线性插值，还有多项式插值，样条插值等等等等，所以加权除了线性进行加权，也可以用多项式等方式计算加权系数。
  多模型的控制除了加权的方式，还有一种更粗暴的选择的方式，也就是判断当前条件和哪个建立好模型的条件接近，直接把近似条件的模型搬过来用，这就需要控制算法本身对模型不匹配有一定自适应性了（控制术语来说的话，叫“鲁棒性”）。

提升DMC的抗扰动能力

  DMC的预测模型，反应的是控制量和对象输出的关系，但实际的控制对象，往往会受到一些干扰因素影响。前面反馈校正改进中我们提到，反馈校正环节只能校正模型的误差，干扰产生的误差也是一个随机数，DMC本身无法消除它。
  减小甚至消除干扰的影响，有不同的思路。
  首先是将扰动包含在预测过程中，我们计算${\tilde{\mathbf{y}}}_{\mathbf{PM}}\left(\mathbf{k}\right)$时，如果能够对干扰如何影响到未来输出也预测到，那么控制律计算中也就可以考虑到干扰的作用进行调整，也就是说控制量可以和干扰配合，共同作用使得对象输出达到期望值（这时由于预测模型及控制律已经包含了对扰动的预测，反馈校正环节就不需要考虑这个干扰了）。我想把一块石头挪动两米，我知道过一会有个人来捣乱挪一米，我可以利用这一点，只需要提前挪一米就可以了。这种思想是很好的，但是要求未来的干扰可以预测，这实际上几乎不可能实现，所以这种思路退而求其次，只要求当前的干扰能够测量到即可，在计算控制律时考虑一下当前的扰动。总结一下，这种思路针对的是可以测量到，但是我们无法提前防止它输入到对象的可测不可控干扰。
  实际生活中，干扰往往非常复杂，很难简单概括它的来源并测量大小，也就是说实际面临的干扰很多是不可预知并且无法测量的不可测不可控干扰，这种干扰我们要尽可能在它影响到对象前提前消除掉。干扰主动抑制的思路，利用对象的输入输出和模型估计干扰，并通过补偿环节抵消，这个方向研究成果很多，扰动观测器，扩张状态观测器，不确定性扰动估计器，等价输入干扰估计器，未知输入观测器等等等等。
  这两种抗扰动的思路一种像是在生病之后对症下药，一种像是在生病前提前预防，防止它影响到身体健康。

结语

  上篇介绍了基本的DMC算法，中篇提到一些算法的改进，接下来在下篇中，我们将通过Matlab编程进行仿真，来观察不同参数对控制性能的影响，并实现部分改进算法，与基本的DMC进行比较。