【预测控制3】动态矩阵控制（DMC）及Matlab仿真

---

相关文章:
【预测控制1】模型预测控制概论
【预测控制2】预测控制基本原理
【预测控制4】基于状态空间模型的预测控制

---

前言

动态矩阵控制（Dynamic Matrix Control，DMC）起源于20世纪70年代，是一种基于数学模型的先进控制算法，在化工、电力、冶金、制药等领域得到了广泛的应用。传统控制算法主要是基于PID控制器，但对于复杂的动态过程，PID控制器并不能提供良好的控制效果。相比之下，DMC控制算法的控制性能更好，可以适应更为复杂的过程控制。

特别说明本文中介绍的是SISO系统且无约束条件的DMC算法

---

一、基于阶跃响应的动态矩阵控制

 DMC是一种基于阶跃响应的预测控制算法，因而适用于渐进稳定的线性系统。对于弱非线性对象，可以先对其线性化，再采用DMC，而对于不稳定对象，可以先用PID控制使其稳定，然后再使用DMC算法。
 简单的说，DMC就是通过阶跃响应模型预测系统未来输出，然后求解优化问题得到控制增量，将控制增量的第一个元素作用于系统。
 DMC算法包括预测模型、滚动优化和反馈校正三个部分组成，接下来分别介绍三个部分。

1.预测模型

 预测模型，就是通过系统的阶跃响应模型和线性系统的叠加性原理来预测对象在未来时刻的输出值。
  在DMC算法中，是通过系统的单位阶跃响应模型来预测对象的未来输出。 首先需要得到预测对象的单位阶跃响应的采样值$a_i=a(iT),i=1,2,...$,其中$T$为采样时间。对于渐进稳定的系统，阶跃响应会在某一时刻$t_N=NT$趋于稳定平稳,在$TN$时刻后的采样值$a_N，a_{N+1},....$都趋于稳态值，所以可以认为$a_N=a_{N+1}=a_{N+2}=....$。这样，对象的动态信息就可以近似地用有限集合 { a 1 , a 2 , . . . . , a N − 1 , a N } \{a_1,a_2,....,a_{N-1},a_N\} {a1​,a2​,....,aN−1​,aN​}加以描述，这个集合的参数构成了DMC的模型参数，向量$a=[a_1,a_2,....,a_{N-1},a_N{]}^T$称为模型向量,也是系统的阶跃响应模型，$N$则称为建模时域。
 举个例子，假如100斤的小明想要把体重增长到120斤，他想制定一个计划，那么他肯定需要结合他之前每天吃多少饭涨了多少斤来计划以后吃多少。假如小明之前连续15天吃3碗饭，每一天都记录他的体重。在第10天的时候，小明的体重首次到达了105斤，第15天，小明的体重依然是105斤，我们可以认为小明吃3碗饭的稳态就是105斤，但这远远不符合他希望达到的120斤，小明就可以根据前10天体重变化为依据，来预测他今后吃多少饭在几天过后能达到120斤。而这前10天的体重就相当于模型向量$a$,10天为建模时域$N$

 下面我们就可以根据系统的模型向量$a$和线性系统的比例和叠加性质来预测对象在未来时刻的输出值。
 首先，在k时刻，假定控制作用保持不变时对未来$N$个时刻的输出具有初始预测值$\widetilde{y_0}(k+1|k),i=1,...,N$，其中$k+i|k$表示在k时刻对$k+i$时刻的预测。
 当$k$时刻输入$u$的增量为$\Delta u(k)$,在其作用下未来时刻输出预测值为$${\tilde{y}}_1(k+i|k)=\widetilde{y_0}(k+1|k)+a_i\Delta u(k).$$当$k+1$时刻输入$u$再变为$\Delta u(k+1)$，则未来输入预测值为$$\begin{array}{rl}{\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)+a_1\Delta u(k+1)\\ ⋮& \\ {\tilde{y}}_1(k+N|k)& ={\tilde{y}}_0(k+N|k)+a_N\Delta u(k)+a_{N-1}\Delta u(k+1)\end{array}$$
以此类推，连续M个输入$\Delta u(k),...,\Delta u(k+M-1)$作用下，未来各时刻的输出值为$${\tilde{y}}_M(k+i|k)={\tilde{y}}_0(k+i|k)+\sum _{j=1}^{min(M,i)}{a}_{i-j+1}\Delta u(k+j-1),i=1,...N\text{\hspace{1em}(1-1)}$$
其中，$y$的下标表示控制量变化的次数，即$y_M$表示控制量$\Delta u$连续变化了$M$次。根据上面的式子，只要我们知道对象的初始预测值${\tilde{y}}_0(k+i|k)$，就可以根据系统的模型向量$a$和控制增量来预测对象的未来输出。式(1-1)就被称为预测模型。

2.滚动优化

 当我们建立好预测模型后，现在的问题就是我们要怎么计算$\Delta u(k),...,\Delta u(k+M-1)$的最优值，来让系统的输出达到期望值，这就是滚动优化的核心问题。也就是说在每一时刻$k$，要确定从该时刻起的$M$个控制增量$\Delta u(k),...,\Delta u(k+M-1)$，使我们的系统的输出预测值${\tilde{y}}_M(k+i|k)$在未来P时刻尽可能接近给定的期望值$\omega (k+i),i=1,...,P$。这里的$M、P$分别被称为控制时域和优化时域（$M⩽P⩽N$）
 继续小明的增肥计划，小明根据自己的模型向量$a$，计划5天分布增加饭量$\Delta u(1),....,\Delta u(5)$,然后预测接下7天的体重，让其尽可能的接近120斤。这里的5天相当于控制时域，7天相当于优化时域。那么如何求出$\Delta u(1),....,\Delta u(5)$的最优值呢？
根据最优控制，我们可以写出在$k$时刻的性能指标$$minJ(k)=\sum _{i=1}^P{q}_i[\omega (k+i)-{\tilde{y}}_M(k+i|k){]}^2$$除了要求让输出达到期望值以外，我们还要求不要让$\Delta u(k)$剧烈的变化，因此我们再在性能指标后加一项来软约束$\Delta u(k)$，最终性能指标可取$$minJ(k)=\sum _{i=1}^P{q}_i[\omega (k+i)-{\tilde{y}}_M(k+i|k){]}^2+\sum _{j=1}^M{r}_j\Delta u^2(k+j-1)\text{\hspace{1em}(2-1)}$$其中，$q_i、r_j$是权系数，分别表示对跟踪误差及控制量变化的抑制程度，$q_i$值越大,表明我们期望对应的控制输出越接近期望值；$r_j$的值越大，表明我们期望对应的控制动作变化越小。
 上述问题就是动态模型(1-1)约束下，以$\Delta u_M(k)=[\Delta u(k),...,\Delta u(k+M-1){]}^T$为优化变量，使性能指标(2-1)最小的优化问题。
 现在我们回过头来将式(1-1)写为向量的形式，其中要强调的是，优化时域为P，表示我们要使预测输出${\tilde{y}}_M(k+1|k),...,{\tilde{y}}_M(k+P|k)$尽可能的接近我们的期望输出；控制时域为M，表示我们有连续M个输入$\Delta u(k),...,\Delta u(k+M-1)$作用，因此可以将式(1-1)展开写为$$\begin{array}{rl}{\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)\\ {\tilde{y}}_M(k+3|k)& ={\tilde{y}}_0(k+3|k)+a_3\Delta u(k)+a_2\Delta u(k+1)+a_1\Delta u(k+2)\\ & ......\\ {\tilde{y}}_M(k+P|k)& ={\tilde{y}}_0(k+P|k)+a_P\Delta u(k)+a_{P-1}\Delta u(k+1)+...+a_{P-M+1}\Delta u(k+M-1)\end{array}$$现在定义$$\begin{gathered} {\tilde{Y}}_{PM}(k)=\left[\begin{array}{c}{\tilde{y}}_M(k+1|k)\\ {\tilde{y}}_M(k+2|k)\\ ...\\ {\tilde{y}}_M(k+P|k)\end{array}\right],{\tilde{Y}}_{P0}(k)=\left[\begin{array}{c}{\tilde{y}}_0(k+1|k)\\ {\tilde{y}}_0(k+2|k)\\ ...\\ {\tilde{y}}_0(k+P|k)\end{array}\right], \\ \Delta U_M(k)=\left[\begin{array}{c}\Delta u(k)\\ \Delta u(k+1)\\ ...\\ \Delta u(k+M-1)\end{array}\right],A=\left[\begin{array}{cccc}{a}_1& 0& ...& 0\\ a_2& a_1& ...& 0\\ ⋮& ⋮& \ddots & ⋮\\ a_M& a_{M-1}& ...& a_1\\ ⋮& ⋮& & ⋮\\ a_P& a_{P-1}& ...& a_{P-M+1}\end{array}\right]. \end{gathered}$$其中 $A$被称为动态矩阵。
因此，可以将式(1-1)写为向量形式$${\tilde{Y}}_{PM}={\tilde{Y}}_{P0}+A\Delta U_M(k).\text{\hspace{1em}(2-2)}$$同样，性能指标(2-1)也可以写为向量形式$$minJ(k)={\Vert \begin{array}{c}{\omega }_P(k)-{\tilde{Y}}_{PM}(k)\end{array}\Vert }_Q^2+{\Vert \begin{array}{c}\Delta U_M(k)\end{array}\Vert }_R^2\text{\hspace{1em}(2-3)}$$其中，$$\begin{gathered} {\omega }_P(k)=[{\omega }_P(k+1),...,{\omega }_P(k+P){]}^T, \\ Q=diag(q_1,...,q_P), \\ R=diag(r_1,...,r_M). \end{gathered}$$由权系数构成的对角矩阵$Q、R$分别被称为误差权矩阵和控制权矩阵。
将式(2-2)带入式(2-3)可得，$$minJ(k)={\Vert \begin{array}{c}{\omega }_P(k)-{\tilde{Y}}_{P0}-A\Delta U_M(k)\end{array}\Vert }_Q^2+{\Vert \begin{array}{c}\Delta U_M(k)\end{array}\Vert }_R^2$$为了求$\Delta U_M(k)$满足上式的最小值，对上式关于$\Delta U_M(k)$求导并令其等于0可以得到，$$\Delta U_M(k)=(A^{T}QA+R{)}^{-1}{A}^{T}Q[{\omega }_P(k)-{\tilde{Y}}_{P0}(k)]$$
但我们不把它们全部实施，我们只取$\Delta U_M(k)$的首元素来作用于对象，可以表示为，$$\Delta u(k)=c^T\Delta U_M(k)=c^T(A^{T}QA+R{)}^{-1}{A}^{T}Q[{\omega }_P(k)-{\tilde{Y}}_{P0}(k)]$$其中，$c^T=\left[\begin{array}{cccc}1& 0& ...& 0\end{array}\right]$。然后定义控制增益$$d^T=c^T(A^{T}QA+R{)}^{-1}{A}^{T}Q$$最终可以得到$$\Delta u(k)=d^T[{\omega }_P(k)-{\tilde{Y}}_{P0}(k)]\text{\hspace{1em}(2-4)}$$在求出控制增量$\Delta u(k)$后，当前k时刻的实际控制量为$$u(k)=u(k-1)+\Delta u(k)\text{\hspace{1em}(2-5)}$$然后将$u(k)$作用于对象，到下一时刻，又以$k+1$取代$k$提出同样的优化问题求出$\Delta u(k+1)$,得到$u(k+1)$作用于对象。如此滚动的进行，这就是滚动优化的含义。

3.反馈校正

 当我们在k时刻将$u(k)$作用于对象之后，相当于我们在$k$时刻给对象的输入端加上了一个幅值为$\Delta u(k)$的阶跃信号，根据预测模型(1-1)，我们可以得到在接下来$N$个时刻的输出${\tilde{y}}_{N1}(k+1|k),...,{\tilde{y}}_{N1}(k+N|k)$，即$${\tilde{Y}}_{N1}(k)={\tilde{Y}}_{N0}(k)+a\Delta u(k),\text{\hspace{1em}(3-1)}$$其中，$$\begin{gathered} {\tilde{Y}}_{N1}(k)=\left[\begin{array}{c}{\tilde{y}}_{N1}(k+1|k)\\ {\tilde{y}}_{N1}(k+2|k)\\ ...\\ {\tilde{y}}_{N1}(k+N|k)\end{array}\right],{\tilde{Y}}_{N0}(k)=\left[\begin{array}{c}{\tilde{y}}_{N0}(k+1|k)\\ {\tilde{y}}_{N0}(k+2|k)\\ ...\\ {\tilde{y}}_{N0}(k+N|k)\end{array}\right], \\ a=[a_1,a_2,....,a_{N-1},a_N{]}^T. \end{gathered}$$${\tilde{y}}_{N0}$表示控制作用保持不变时接下来N个时刻的预测输出，${\tilde{y}}_{N1}$则表示在当前时刻施加了一个控制增量$\Delta u(k)$后接下来N个时刻的预测输出。
 但是在实际情况中，由于环境因素等的干扰，系统的实际输出可能会比我们预测的输出偏小或偏大，就像小明他有可能会感冒生病，导致他吃相同的饭量会比以前增长的体重少。这种情况我们必须要考虑，因为我们下一步的优化是建立在这一次预测输出之上，如果我们这一次预测输出出现了偏差，那么我们下一步优化偏差值可能更大。那么我们要如何减小这个偏差呢?在DMC算法中，我们首先将我们预测到的下一时刻的输出${\tilde{y}}_1(k+1|k)$和系统的实际输出$y(k+1)$相比较,得到偏差值$$e(k+1)=y(k+1)-{\tilde{y}}_1(k+1|k)\text{\hspace{1em}(3-2)}$$在上面的式子中，系统的实际输出$y(k+1)$是已知的，我们的预测输出${\tilde{y}}_1(k+1|k)$也是已知的，说明偏差$e(k+1)$也是已知的，DMC就想，既然偏差都已知了，那我们直接把这个偏差$e(k+1)$补偿到我们的接下来的预测输出${\tilde{Y}}_{N1}(k)$不就可以减小偏差了吗，因此可以得到$${\tilde{Y}}_{cor}(k+1)={\tilde{Y}}_{N1}(k)+he(k+1),\text{\hspace{1em}(3-3)}$$其中，$${\tilde{Y}}_{cor}(k+1)=\left[\begin{array}{c}{\tilde{y}}_{cor}(k+1|k+1)\\ {\tilde{y}}_{cor}(k+2|k+1)\\ ...\\ {\tilde{y}}_{cor}(k+N|k+1)\end{array}\right]$$为校正后的输出预测量，有权系数组成的$N$维向量$h=[h_1,...h_N{]}^T$称为校正量。
 以上就是反馈校正的主要思想。

---

 现在我们都知道，预测模型是需要我们提前准备好的，而滚动优化和反馈校正是每一次循环都要进行的，那我们如何将滚动优化和反馈校正两部分形成"闭环"呢？回顾公式，我们需要把$k$时刻经过修正的预测输出${\tilde{Y}}_{cor}(k+1)$作为$k+1$时刻的初始预测值，也就是初值${\tilde{Y}}_{N0}(k+1)$。但我们不能直接令${\tilde{Y}}_{cor}(k+1)={\tilde{Y}}_{N0}(k+1)$。因为对于$k$时刻而言，我们得到的预测输出${\tilde{Y}}_{cor}(k+1)=[{\tilde{y}}_{cor}(k+1|k+1),...,{\tilde{y}}_{cor}(k+N|k+1)]$，而在$k+1$时刻的初始预测值为${\tilde{Y}}_{N0}(k+1)=[{\tilde{y}}_{N0}(k+2|k+1),...,{\tilde{y}}_{N0}(k+N+1|k+1)]$,所以我们令$$\begin{gathered} {\tilde{y}}_{N0}(k+2|k+1)={\tilde{y}}_{cor}(k+2|k+1) \\ {\tilde{y}}_{N0}(k+3|k+1)={\tilde{y}}_{cor}(k+3|k+1) \\ ... \\ {\tilde{y}}_{N0}(k+N|k+1)={\tilde{y}}_{cor}(k+N|k+1) \end{gathered}$$由于${\tilde{Y}}_{cor}(k+1)$中没有${\tilde{y}}_{cor}(k+N+1|k+1)$这一项，那么我们用${\tilde{y}}_{cor}(k+N|k+1)$取近似。最终可以表示为$${\tilde{Y}}_{N0}(k+1)=S{\tilde{Y}}_{cor}(k+1)$$其中，$$S=\left[\begin{array}{ccccc}0& 1& 0& ...& 0\\ 0& 0& 1& ...& 0\\ 0& 0& 0& \ddots & 0\\ ⋮& ⋮& ⋮& 0& 1\\ 0& ...& & 0& 1\end{array}\right]$$有了${\tilde{Y}}_{N0}(k+1)$,就可以取其前P项得到${\tilde{Y}}_{P0}(k+1)$，然后可以根据式子(2-4)计算出$\Delta u(k+1)$，整个控制过程就是滚动优化结合反馈校正反复在线进行的。

二、Matlab仿真示例

1.算法结构

根据上面整个DMC算法的流程，我们可以得到它的算法结构，如下图所示。图中粗箭头表示向量流，细箭头表示数值流。

2.计算流程

流程图中的算式依次对应于式(3-2)、式(3-3)式(2-4)、式(2-5)、式(3-1)

3.仿真示例

例： 考虑如下SISO系统，其传递函数为$$G(s)=\frac{3}{4s+1}$$

程序代码：

clear all;
close all;

P = 6;M = 2;N = 25;%优化时域、控制时域和建模时域
ts = 1;%采样周期
k_step = 50;
times = zeros(1,k_step);

den = [4,1];
num = [3];
order = length(den);
sys = tf(num,den);%建立传递函数
sysd = c2d(sys,ts);
[numd,dend] = tfdata(sysd,'v');
[a,~] = step(num,den,ts:ts:N*ts);%获得模型向量a

S = diag(ones(1,N-1),1);S(N,N)=1;%移位矩阵S
yr = 3;%期望输出

dum = zeros(M,1);%优化变量
du = 0;%控制增量
yn0 = zeros(N,1);%前一时刻预测输出
yn1 = zeros(N,1);%当前预测输出
yp0 = zeros(P,1);%优化输出
us = zeros(1,k_step);%实际控制量
ys = zeros(1,k_step);%实际输出
yk_history = zeros(order-1,1);
uk_history = zeros(order-1,1);

R = eye(M);Q = eye(P);%控制权矩阵和误差权矩阵
h = 0.5*ones(N,1);
h(1) = 1;%校正系数矩阵

A = zeros(P,M);
for i=1:M       
    A(i:P,i) = a(1:P-i+1);
end  %建立动态矩阵A
c = zeros(M,1);c(1) = 1;
d = c'*inv(A'*Q*A+R)*A'*Q;d = d';%控制增益

%第一步预测
times(1) = ts;
ys(1) = 0;
ek = ys(1)-yn1(1);
yn0 = S*(yn1+h*ek);%预测值校正
yp0 = yn0(1:P);
du = d'*(yr-yp0);%计算这一时刻控制增量
us(1) = 0+du;%施加控制增量，得到实际控制量
uk = us(1);
uk_history(1) = uk;
yn1 = yn0+a(1:N)*du;%得到施加控制增量后的预测输出

%滚动优化
for k=2:k_step
    times(k) = k*ts;
    ys(k) = -dend(2:end)*yk_history + numd(2:end)*uk_history;%得到系统的实际输出
    yk_history(2:end) = yk_history(1:end-1);
    yk_history(1) = ys(k);
    ek = ys(k)-yn1(1);%计算偏差,公式(3-2)
    yn0 = S*(yn1+h*ek);%预测值校正,公式(3-3)
    yp0 = yn0(1:P);
    du = d'*(yr-yp0);%计算这一时刻控制增量,公式(2-4)
    us(k) = us(k-1)+du;%施加控制增量，得到实际控制量,公式(2-5)
    uk_history(2:end) = uk_history(1:end-1);
    uk_history(1) = us(k);
    yn1 = yn0+a(1:N)*du;%得到施加控制增量后的预测输出,公式(3-1)
end

figure(1)
plot(times,ys');
title('输出曲线')
figure(2)
stairs(times,us');
title('控制量变化曲线')

运行结果：

---

总结

DMC算法（Dynamic Matrix Control）是一种经典的模型预测控制算法，适用于线性、纯时滞、开环渐进稳定的非最小相位系统，其优缺点如下：
优点：

1. 较好的控制性能：DMC算法能够很好地解决传统PID控制算法难以处理的复杂控制问题，能够提供更好的控制性能。
2. 可以处理多变量控制问题：DMC算法可以对多个变量同时进行控制，能够处理多变量控制问题。
3. 具有鲁棒性：DMC算法对于模型不确定性、测量噪声等影响具有很强的鲁棒性。

缺点：

1. 计算量大：DMC算法计算比较复杂，需要大量的计算资源，处理速度较慢。
2. 模型依赖性强：DMC算法的控制性能依赖于模型的精度，如果模型不准确，控制效果会受到影响。
3. 参数调节困难：DMC算法需要调节多个参数，参数调节比较困难，需要一定的经验和技巧。