前文链接

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

基础DMC算法

Matlab源码

clear;clc;                              % 清除缓存变量，清除命令行信息
dt=1;                                   % 采样周期
%################################ 基本DMC #####################################%
%% 被控对象参数
K=1;                                    % 开环增益
T=215;                                  % 惯性时间
%% DMC参数
N=2000;                                 % 建模时域
P=500;                                  % 预测时域
M=1;                                    % 控制时域
q=1;                                    % Q矩阵（预测输出与期望值误差权矩阵）中的一个权值
lamda=1;                                % λ(LAMDA)矩阵（控制增量权矩阵）中的一个权值
h=1;                                    % H矩阵（预测与实际输出误差校正系数矩阵）中的一个系数
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
H=h*ones(P,1);                          % 生成反馈校正H矩阵
S=zeros(P,P);                           % 生成移位矩阵S
for i=1:P-1
    S(i,i+1)=1;                         % 次对角线元素为1
end
S(P,P-1)=0; S(P,P)=1;                   % 右下角元素为1
%% 阶跃响应模型
num=K;                                  % 传递函数的分子部分
den=[T 1];                              % 传递函数的分母部分
sys=tf(num,den);                        % 构成传递函数
a=step(sys,0:dt:N);                     % 参数依次为要施加阶跃输入的传递函数，起始时间：步长(周期)：截止时间
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 仿真参数设置
SimuStepCount=2999;                     % 总仿真步数
sp=1;                                   % 设定值
SP=sp*ones(P,1);                        % 设定值向量
%% 开始仿真
Y_DMC=zeros(SimuStepCount+1,1);         % 保存每一步对象输出
U_DMC=zeros(SimuStepCount+1,1);         % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);             % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC(i+1,1)=y;                     % 存储到Y_DMC
    U_DMC(i,1)=u;                       % 存储到U_DMC
end
U_DMC(SimuStepCount+1,1)=U_DMC(SimuStepCount,1);
%% 绘图
n=size(Y_DMC,1);
subplot(2,1,1);
t=1:1:n;
plot(t,Y_DMC(1:n),'r');
title('y')
subplot(2,1,2);
plot(t,U_DMC(1:n),'b--');
title('u')

代码解析

  DMC算法要设置各种参数，设置参数之后，还要计算相关的参数矩阵，导致上面的代码看上去很长，其实关键的控制量计算几行就可以实现（for循环中的内容）。
  上面代码中设置参数的部分我们直接跳过，变量名与之前文章中的符号都是对应的。
  被控对象为一个不带纯迟延环节的一阶惯性对象，用传递函数来表示是 $\frac{K}{Ts+1}$，matlab中的tf()函数可以将分子和分母表示的多项式（从最高次项一直到常数也就是0次项）转换为传递函数形式，例如下面代码会生成一个$\frac{As+1}{B{s}^2+Cs+1}$形式的传递函数：

num=[A 1];
den=[B C 1];
sys=tf(num,den);

  step()函数可以生成DMC需要的阶跃响应模型，我们在参数中提供之前生成的传递函数，以及起始结束时间和步长，函数返回的结果就是单位阶跃输入下的阶跃响应模型$\mathbf{a}$。
  得到$\mathbf{a}$之后就可以计算动态矩阵$\mathbf{A}$了，$\mathbf{A}$再结合已经确定的参数矩阵，就可以计算出控制律中前半部分。这一点之前文章中也提到过，控制律不需要每一步完全重新计算，只有后半部分需要每一步重新计算，这里固定的前半部分我用变量D表示。
  接下来终于到了仿真中关键的循环部分，Y_DMC和U_DMC用来存储每一步的输出和控制量，Y_P0为${\tilde{\mathbf{y}}}_{\mathbf{P}0}$，Y_P1为${\tilde{\mathbf{y}}}_{\mathbf{P}1}$，Y_COR为${\tilde{\mathbf{y}}}_{\text{COR}}$。上篇中提到过，为了让程序第一个时刻可以走下去，第一个时刻先将这些预测值，控制量等全部初始化为0。有人会较真，预测值必须初始化为0吗？澄清一下，预测值初始化为1甚至1000000都可以，由于这里仿真的时候对象输出是从0开始变化的，所以初始化为0保证了第一个时刻预测值和实际值的误差$e$不会太大，导致不必要的控制过程波动，因为第一步的$e$很大，会导致最初的几个时刻反复大幅度地修正预测值。实际中，例如实际的锅炉温度在500℃左右，那么几个预测值初始化在500℃显然更好。
  在进入for循环前，理论上应该先求一次y_P1：

y=0;                                    % 对象当前输出
y_P1=y_P0+a(1:P,1)*du; % ############## 理论上应该在这个位置求一次初始的y_P1
for i=1:1:SimuStepCount

  由于第一个时刻前du为0，所以y_P1和y_P0是相等的，这里就省略了这一行代码，严格按照算法的流程，应该是有这一行的。
  进入for循环，现在是$k$时刻，先计算$e$，然后计算修正后的预测值${\tilde{\mathbf{y}}}_{\text{COR}}$，移位后得到对$k+1$时刻未施加控制增量的预测值${\tilde{\mathbf{y}}}_{\mathbf{P}0}$，求解控制律得到du,叠加在之前的控制量u之上，得到$k$时刻的控制量施加给对象，同时计算施加du之后的预测值，再次强调一下，标蓝的文字都是在$k$时刻的计算过程，得到的是$k$时刻的控制量，但是$k$时刻控制量作用到对象之后，$k+1$时刻才能收到对象反馈值，时序非常关键，哪个时刻发生了哪些事一定不要弄混！y_P1=y_P0+a(1:P,1)*du;这一行倒是无所谓，放在$k$时刻计算出控制律后，或者放在求$e$之前都可以，只不过放在计算控制律之后表示求y_P1发生在$k$时刻，放在求$e$之前表示求y_P1发生在$k+1$时刻。放在求$e$之前的话这段代码就变成这样：

for i=1:1:SimuStepCount
    % DMC计算部分
    y_P1=y_P0+a(1:P,1)*du;
    e=Y_DMC(i,1)-y_P1(1,1);
    y_COR=y_P1+H*e;
    y_P0=S*y_COR;
    du=D*(SP-y_P0);
    u=du+u;

  这里体现了算法和编程的区别，算法的推导一定要求严谨，不能缺少步骤，步骤的顺序更不能随意调换，而实际的编程一些步骤缺少或者步骤间调换顺序并不会影响计算。
  如果是实际的被控对象，计算出每一步的控制量u就可以，但仿真需要我们自己实现被控对象从输入到输出的过程。

dx=K/T*u+(-1/T)*y;                  % 计算dx
y=y+dx*dt;                          % 计算输出

  传递函数表示的是连续的被控对象，计算机的控制都是离散的。（类似于数学中积分的概念，人工计算可以计算连续函数的积分，计算机计算只能将函数切割成一小段一小段累加来得到积分值。）在控制过程中我们要将被控对象进行离散化，这里使用的是离散相似法中的分环节离散的方式，关于被控对象离散化可以参考相关资料，一阶对象只需要修改上面代码的K和T即可。
  Y_DMC(i+1,1)=y;和U_DMC(i,1)=u;这里刻意明确了一下索引，y是$k+1$时刻的对象反馈，u是$k$时刻的控制量。
  循环结束，我们得到了整个控制过程的对象输出和控制量，绘图就可以观察结果。

仿真结果展示

参数对性能的影响

参数P对性能影响

Matlab源码

clear;clc;                              % 清除缓存变量，清除命令行信息
dt=1;                                   % 采样周期
%################################ 基本DMC P=500 #####################################%
%% 被控对象参数
K=1;                                    % 开环增益
T=215;                                  % 惯性时间
%% DMC参数
N=2000;                                 % 建模时域
P=500;                                  % 预测时域
M=1;                                    % 控制时域
q=1;                                    % Q矩阵（预测输出与期望值误差权矩阵）中的一个权值
lamda=1;                                % λ(LAMDA)矩阵（控制增量权矩阵）中的一个权值
h=1;                                    % H矩阵（预测与实际输出误差校正系数矩阵）中的一个系数
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
H=h*ones(P,1);                          % 生成反馈校正H矩阵
S=zeros(P,P);                           % 生成移位矩阵S
for i=1:P-1
    S(i,i+1)=1;                         % 次对角线元素为1
end
S(P,P-1)=0; S(P,P)=1;                   % 右下角元素为1
%% 阶跃响应模型
num=K;                                  % 传递函数的分子部分
den=[T 1];                              % 传递函数的分母部分
sys=tf(num,den);                        % 构成传递函数
a=step(sys,0:dt:N);                     % 参数依次为要施加阶跃输入的传递函数，起始时间：步长(周期)：截止时间
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 仿真参数设置
SimuStepCount=2999;                     % 总仿真步数
sp=1;                                   % 设定值
SP=sp*ones(P,1);                        % 设定值向量
%% 开始仿真
Y_DMC=zeros(SimuStepCount+1,1);         % 保存每一步对象输出
U_DMC=zeros(SimuStepCount+1,1);         % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC(i+1,1)=y;                     % 存储到Y_DMC
    U_DMC(i,1)=u;                       % 存储到U_DMC
end
U_DMC(SimuStepCount+1,1)=U_DMC(SimuStepCount,1);
%% ################################ 基本DMC P=400 #####################################%
P=400;                                  % 预测时域
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
H=h*ones(P,1);                          % 生成反馈校正H矩阵
S=zeros(P,P);                           % 生成移位矩阵S
for i=1:P-1
    S(i,i+1)=1;                         % 次对角线元素为1
end
S(P,P-1)=0; S(P,P)=1;                   % 右下角元素为1
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 仿真参数设置
SP=sp*ones(P,1);                        % 设定值向量
%% 开始仿真
Y_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步对象输出
U_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC2(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC2(i,1)=u;                      % 存储到U_DMC
end
U_DMC2(SimuStepCount+1,1)=U_DMC2(SimuStepCount,1);
%% ################################ 基本DMC P=300 #####################################%
P=300;                                  % 预测时域
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
H=h*ones(P,1);                          % 生成反馈校正H矩阵
S=zeros(P,P);                           % 生成移位矩阵S
for i=1:P-1
    S(i,i+1)=1;                         % 次对角线元素为1
end
S(P,P-1)=0; S(P,P)=1;                   % 右下角元素为1
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 仿真参数设置
SP=sp*ones(P,1);                        % 设定值向量
%% 开始仿真
Y_DMC3=zeros(SimuStepCount+1,1);        % 保存每一步对象输出
U_DMC3=zeros(SimuStepCount+1,1);        % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC3(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC3(i,1)=u;                      % 存储到U_DMC
end
U_DMC3(SimuStepCount+1,1)=U_DMC3(SimuStepCount,1);
%% 绘图
n=size(Y_DMC,1);
subplot(2,1,1);
t=1:1:n;
plot(t,Y_DMC(1:n),'r');
hold on;
plot(t,Y_DMC2(1:n),'b--');
hold on;
plot(t,Y_DMC3(1:n),'k-.');
legend('P=500','P=400','P=300')
title('y')
subplot(2,1,2);
plot(t,U_DMC(1:n),'r');
hold on;
plot(t,U_DMC2(1:n),'b--');
hold on;
plot(t,U_DMC3(1:n),'k-.');
legend('P=500','P=400','P=300')
title('u')

仿真结果

参数M对性能影响

Matlab源码

clear;clc;                              % 清除缓存变量，清除命令行信息
dt=1;                                   % 采样周期
%################################ 基本DMC M=1 #####################################%
%% 被控对象参数
K=1;                                    % 开环增益
T=215;                                  % 惯性时间
%% DMC参数
N=2000;                                 % 建模时域
P=500;                                  % 预测时域
M=1;                                    % 控制时域
q=1;                                    % Q矩阵（预测输出与期望值误差权矩阵）中的一个权值
lamda=1;                                % λ(LAMDA)矩阵（控制增量权矩阵）中的一个权值
h=1;                                    % H矩阵（预测与实际输出误差校正系数矩阵）中的一个系数
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
H=h*ones(P,1);                          % 生成反馈校正H矩阵
S=zeros(P,P);                           % 生成移位矩阵S
for i=1:P-1
    S(i,i+1)=1;                         % 次对角线元素为1
end
S(P,P-1)=0; S(P,P)=1;                   % 右下角元素为1
%% 阶跃响应模型
num=K;                                  % 传递函数的分子部分
den=[T 1];                              % 传递函数的分母部分
sys=tf(num,den);                        % 构成传递函数
a=step(sys,0:dt:N);                     % 参数依次为要施加阶跃输入的传递函数，起始时间：步长(周期)：截止时间
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 仿真参数设置
SimuStepCount=2999;                     % 总仿真步数
sp=1;                                   % 设定值
SP=sp*ones(P,1);                        % 设定值向量
%% 开始仿真
Y_DMC=zeros(SimuStepCount+1,1);         % 保存每一步对象输出
U_DMC=zeros(SimuStepCount+1,1);         % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC(i+1,1)=y;                     % 存储到Y_DMC
    U_DMC(i,1)=u;                       % 存储到U_DMC
end
U_DMC(SimuStepCount+1,1)=U_DMC(SimuStepCount,1);
%% ################################ 基本DMC M=2 #####################################%
M=2;
%% 根据DMC参数计算DMC参数矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 开始仿真
Y_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步对象输出
U_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC2(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC2(i,1)=u;                      % 存储到U_DMC
end
U_DMC2(SimuStepCount+1,1)=U_DMC2(SimuStepCount,1);
%% ################################ 基本DMC M=5 #####################################%
M=5;
%% 根据DMC参数计算DMC参数矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 开始仿真
Y_DMC3=zeros(SimuStepCount+1,1);        % 保存每一步对象输出
U_DMC3=zeros(SimuStepCount+1,1);        % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC3(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC3(i,1)=u;                      % 存储到U_DMC
end
U_DMC3(SimuStepCount+1,1)=U_DMC3(SimuStepCount,1);
%% 绘图
n=size(Y_DMC,1);
subplot(2,1,1);
t=1:1:n;
plot(t,Y_DMC(1:n),'r');
hold on;
plot(t,Y_DMC2(1:n),'b--');
hold on;
plot(t,Y_DMC3(1:n),'k-.');
legend('M=1','M=2','M=5')
title('y')
subplot(2,1,2);
plot(t,U_DMC(1:n),'r');
hold on;
plot(t,U_DMC2(1:n),'b--');
hold on;
plot(t,U_DMC3(1:n),'k-.');
legend('M=1','M=2','M=5')
title('u')

仿真结果

参数q对性能影响

Matlab源码

clear;clc;                              % 清除缓存变量，清除命令行信息
dt=1;                                   % 采样周期
%################################ 基本DMC q=1 #####################################%
%% 被控对象参数
K=1;                                    % 开环增益
T=215;                                  % 惯性时间
%% DMC参数
N=2000;                                 % 建模时域
P=500;                                  % 预测时域
M=1;                                    % 控制时域
q=1;                                    % Q矩阵（预测输出与期望值误差权矩阵）中的一个权值
lamda=1;                                % λ(LAMDA)矩阵（控制增量权矩阵）中的一个权值
h=1;                                    % H矩阵（预测与实际输出误差校正系数矩阵）中的一个系数
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
H=h*ones(P,1);                          % 生成反馈校正H矩阵
S=zeros(P,P);                           % 生成移位矩阵S
for i=1:P-1
    S(i,i+1)=1;                         % 次对角线元素为1
end
S(P,P-1)=0; S(P,P)=1;                   % 右下角元素为1
%% 阶跃响应模型
num=K;                                  % 传递函数的分子部分
den=[T 1];                              % 传递函数的分母部分
sys=tf(num,den);                        % 构成传递函数
a=step(sys,0:dt:N);                     % 参数依次为要施加阶跃输入的传递函数，起始时间：步长(周期)：截止时间
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 仿真参数设置
SimuStepCount=2999;                     % 总仿真步数
sp=1;                                   % 设定值
SP=sp*ones(P,1);                        % 设定值向量
%% 开始仿真
Y_DMC=zeros(SimuStepCount+1,1);         % 保存每一步对象输出
U_DMC=zeros(SimuStepCount+1,1);         % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC(i+1,1)=y;                     % 存储到Y_DMC
    U_DMC(i,1)=u;                       % 存储到U_DMC
end
U_DMC(SimuStepCount+1,1)=U_DMC(SimuStepCount,1);
%% ################################ 基本DMC q=0.0001 #####################################%
q=0.0001;
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 开始仿真
Y_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步对象输出
U_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC2(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC2(i,1)=u;                      % 存储到U_DMC
end
U_DMC2(SimuStepCount+1,1)=U_DMC2(SimuStepCount,1);
%% ################################ 基本DMC q=0.00005 #####################################%
q=0.00005;
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 开始仿真
Y_DMC3=zeros(SimuStepCount+1,1);        % 保存每一步对象输出
U_DMC3=zeros(SimuStepCount+1,1);        % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC3(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC3(i,1)=u;                      % 存储到U_DMC
end
U_DMC3(SimuStepCount+1,1)=U_DMC3(SimuStepCount,1);
%% 绘图
n=size(Y_DMC,1);
subplot(2,1,1);
t=1:1:n;
plot(t,Y_DMC(1:n),'r');
hold on;
plot(t,Y_DMC2(1:n),'b--');
hold on;
plot(t,Y_DMC3(1:n),'k-.');
legend('q=1','q=0.0001','q=0.00005')
title('y')
subplot(2,1,2);
plot(t,U_DMC(1:n),'r');
hold on;
plot(t,U_DMC2(1:n),'b--');
hold on;
plot(t,U_DMC3(1:n),'k-.');
legend('q=1','q=0.0001','q=0.00005')
title('u')

仿真结果

参数lamda对性能影响

Matlab源码

clear;clc;                              % 清除缓存变量，清除命令行信息
dt=1;                                   % 采样周期
%################################ 基本DMC lamda=1 #####################################%
%% 被控对象参数
K=1;                                    % 开环增益
T=215;                                  % 惯性时间
%% DMC参数
N=2000;                                 % 建模时域
P=500;                                  % 预测时域
M=1;                                    % 控制时域
q=1;                                    % Q矩阵（预测输出与期望值误差权矩阵）中的一个权值
lamda=1;                                % λ(LAMDA)矩阵（控制增量权矩阵）中的一个权值
h=1;                                    % H矩阵（预测与实际输出误差校正系数矩阵）中的一个系数
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
H=h*ones(P,1);                          % 生成反馈校正H矩阵
S=zeros(P,P);                           % 生成移位矩阵S
for i=1:P-1
    S(i,i+1)=1;                         % 次对角线元素为1
end
S(P,P-1)=0; S(P,P)=1;                   % 右下角元素为1
%% 阶跃响应模型
num=K;                                  % 传递函数的分子部分
den=[T 1];                              % 传递函数的分母部分
sys=tf(num,den);                        % 构成传递函数
a=step(sys,0:dt:N);                     % 参数依次为要施加阶跃输入的传递函数，起始时间：步长(周期)：截止时间
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 仿真参数设置
SimuStepCount=2999;                     % 总仿真步数
sp=1;                                   % 设定值
SP=sp*ones(P,1);                        % 设定值向量
%% 开始仿真
Y_DMC=zeros(SimuStepCount+1,1);         % 保存每一步对象输出
U_DMC=zeros(SimuStepCount+1,1);         % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC(i+1,1)=y;                     % 存储到Y_DMC
    U_DMC(i,1)=u;                       % 存储到U_DMC
end
U_DMC(SimuStepCount+1,1)=U_DMC(SimuStepCount,1);
%% ################################ 基本DMC lamda=10000 #####################################%
lamda=10000;
%% 根据DMC参数计算DMC参数矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 开始仿真
Y_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步对象输出
U_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC2(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC2(i,1)=u;                      % 存储到U_DMC
end
U_DMC2(SimuStepCount+1,1)=U_DMC2(SimuStepCount,1);
%% ################################ 基本DMC lamda=50000 #####################################%
lamda=50000;
%% 根据DMC参数计算DMC参数矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 开始仿真
Y_DMC3=zeros(SimuStepCount+1,1);        % 保存每一步对象输出
U_DMC3=zeros(SimuStepCount+1,1);        % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC3(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC3(i,1)=u;                      % 存储到U_DMC
end
U_DMC3(SimuStepCount+1,1)=U_DMC3(SimuStepCount,1);
%% 绘图
n=size(Y_DMC,1);
subplot(2,1,1);
t=1:1:n;
plot(t,Y_DMC(1:n),'r');
hold on;
plot(t,Y_DMC2(1:n),'b--');
hold on;
plot(t,Y_DMC3(1:n),'k-.');
legend('lamda=1','lamda=10000','lamda=50000')
title('y')
subplot(2,1,2);
plot(t,U_DMC(1:n),'r');
hold on;
plot(t,U_DMC2(1:n),'b--');
hold on;
plot(t,U_DMC3(1:n),'k-.');
legend('lamda=1','lamda=10000','lamda=50000')
title('u')

仿真结果

讨论

  对照仿真结果，结合公式推导，我们就可以进一步理解各个参数的作用。
  参数P增大，首先会使得动态矩阵$\mathbf{A}$的行数增加， 导致控制律前半部分的值增大， 但是参数P增大，又会使得控制律后半部分设定值SP与预测值y_P0的差构成向量末尾出现更多的0（可以想象到，P足够大时，设定值与预测值到后面必然会相等，导致继续增大P会使得SP-y_P0末尾出现越来越多的0），所以从控制律公式分析，不好确定P对控制效果的影响。于是我们从仿真结果来看，P的增加导致控制量给的更稳了，对于被控对象来说，这就意味着到达设定值的过程变慢了。
  参数M增大，会使得动态矩阵$\mathbf{A}$的列数增加， 导致控制律前半部分的值增大， 而M对控制律后半部分没有影响，那么会使得根据控制律得到的控制增量增大，也就是说，从控制律公式来看，M增大应该会使得控制量变化更明显，相应地被控对象到达设定值时间会加快。从仿真结果来看，和我们的推测是一致的，M增大时，一开始的控制量会给得更猛。
  参数q和参数lamda从它们的定义中就知道，q越大，被控对象到达设定值会越快，控制量变化会更剧烈，lamda越大，控制量变化会更受限制，使得控制过程更稳定，被控对象到达设定值变慢，仿真结果也是完全对应的。但是还有一个小细节，如果把我的代码里q改成2，3或者lamda改成2，3，会看到几条几乎重合的曲线，实际上你把曲线放大放大再放大，可能要到小数点后四五位，确实也能看到符合上述结论，然后你就明白我上面的q和lamda为什么要设置这么大/这么小了，这提醒了我们，q和lamda确实可以影响控制效果，但是要能够真正“明显”地看到影响，需要先对动态矩阵$\mathbf{A}$以及控制律前半部分的数量级有个估计，进而确定q和lamda在什么数量级。实际上这太麻烦了，实际中q和lamda一般都设成1就懒得调了，如果需要，再自行解决吧O(∩_∩)O。

算法改进效果验证

阶梯式动态矩阵控制与基础DMC对比

  下面的代码如果把M改成1，可以看到DMC和SDMC结果是完全一致的，这里M设置为5方便观察效果，从仿真结果可以看出SDMC将控制量规划为阶梯形式后，削弱了控制量，控制过程更平缓（之前说过，SDMC更重要的作用是取消求逆环节来加快运算速度，用其他编程语言实现SDMC时，可以打印一下运算时间，MATLAB两种算法计算速度都很快，个人认为参考意义不大就没有打印运算时间）。

Matlab源码

clear;clc;                              % 清除缓存变量，清除命令行信息
dt=1;                                   % 采样周期
%################################ 基本DMC #####################################%
%% 被控对象参数
K=1;                                    % 开环增益
T=215;                                  % 惯性时间
%% DMC参数
N=2000;                                 % 建模时域
P=500;                                  % 预测时域
M=5;                                    % 控制时域
q=1;                                    % Q矩阵（预测输出与期望值误差权矩阵）中的一个权值
lamda=1;                                % λ(LAMDA)矩阵（控制增量权矩阵）中的一个权值
h=1;                                    % H矩阵（预测与实际输出误差校正系数矩阵）中的一个系数
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
H=h*ones(P,1);                          % 生成反馈校正H矩阵
S=zeros(P,P);                           % 生成移位矩阵S
for i=1:P-1
    S(i,i+1)=1;                         % 次对角线元素为1
end
S(P,P-1)=0; S(P,P)=1;                   % 右下角元素为1
%% 阶跃响应模型
num=K;                                  % 传递函数的分子部分
den=[T 1];                              % 传递函数的分母部分
sys=tf(num,den);                        % 构成传递函数
a=step(sys,0:dt:N);                     % 参数依次为要施加阶跃输入的传递函数，起始时间：步长(周期)：截止时间
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 仿真参数设置
SimuStepCount=2999;                     % 总仿真步数
sp=1;                                   % 设定值
SP=sp*ones(P,1);                        % 设定值向量
%% 开始仿真
Y_DMC=zeros(SimuStepCount+1,1);         % 保存每一步对象输出
U_DMC=zeros(SimuStepCount+1,1);         % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC(i+1,1)=y;                     % 存储到Y_DMC
    U_DMC(i,1)=u;                       % 存储到U_DMC
end
U_DMC(SimuStepCount+1,1)=U_DMC(SimuStepCount,1);
%% ################################ SDMC #####################################%
beita=0.9;
%% 根据DMC参数计算DMC参数矩阵
beita_all=0;
for i=1:M
    lam(i,1) = beita^(i-1);
    beita_all=beita_all+beita^(2*i-2);
end
G=A*lam;
D=G'*Q/((G'*Q*G)+lamda*beita_all);
%% 开始仿真
Y_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步对象输出
U_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC2(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC2(i,1)=u;                      % 存储到U_DMC
end
U_DMC2(SimuStepCount+1,1)=U_DMC2(SimuStepCount,1);
%% 绘图
n=size(Y_DMC,1);
subplot(2,1,1);
t=1:1:n;
plot(t,Y_DMC(1:n),'r');
hold on;
plot(t,Y_DMC2(1:n),'b--');
legend('DMC','SDMC')
title('y')
subplot(2,1,2);
plot(t,U_DMC(1:n),'r');
hold on;
plot(t,U_DMC2(1:n),'b--');
legend('DMC','SDMC')
title('u')

仿真结果

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

  其实大家从上面一副图中就可以看到，DMC出现了一个明显的“小尖尖”，刚好是在第500步（与预测时域相等）出现的，这就是中篇中提到的截断误差累积导致的。采用改进的移位矩阵，我们可以看到截断误差明显减小了。

Matlab源码

clear;clc;                              % 清除缓存变量，清除命令行信息
dt=1;                                   % 采样周期
%################################ 基本DMC #####################################%
%% 被控对象参数
K=1;                                    % 开环增益
T=215;                                  % 惯性时间
%% DMC参数
N=2000;                                 % 建模时域
P=500;                                  % 预测时域
M=5;                                    % 控制时域
q=1;                                    % Q矩阵（预测输出与期望值误差权矩阵）中的一个权值
lamda=1;                                % λ(LAMDA)矩阵（控制增量权矩阵）中的一个权值
h=1;                                    % H矩阵（预测与实际输出误差校正系数矩阵）中的一个系数
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
H=h*ones(P,1);                          % 生成反馈校正H矩阵
S=zeros(P,P);                           % 生成移位矩阵S
for i=1:P-1
    S(i,i+1)=1;                         % 次对角线元素为1
end
S(P,P-1)=0; S(P,P)=1;                   % 右下角元素为1
%% 阶跃响应模型
num=K;                                  % 传递函数的分子部分
den=[T 1];                              % 传递函数的分母部分
sys=tf(num,den);                        % 构成传递函数
a=step(sys,0:dt:N);                     % 参数依次为要施加阶跃输入的传递函数，起始时间：步长(周期)：截止时间
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 仿真参数设置
SimuStepCount=2999;                     % 总仿真步数
sp=1;                                   % 设定值
SP=sp*ones(P,1);                        % 设定值向量
%% 开始仿真
Y_DMC=zeros(SimuStepCount+1,1);         % 保存每一步对象输出
U_DMC=zeros(SimuStepCount+1,1);         % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC(i+1,1)=y;                     % 存储到Y_DMC
    U_DMC(i,1)=u;                       % 存储到U_DMC
end
U_DMC(SimuStepCount+1,1)=U_DMC(SimuStepCount,1);
%% ################################ 改进移位矩阵 #####################################%
%% 反馈校正环节改进：如果a存在截断误差，移位矩阵需要调整
kesai=(a(P,1)-a(P-1,1))/(a(P-1,1)-a(P-2,1));
if (kesai>0)&&(kesai<1)
    S(P,P-1)=0-kesai;S(P,P)=1+kesai;
end
%% 开始仿真
Y_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步对象输出
U_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    dx=K/T*u+(-1/T)*y;                  % 计算dx
    y=y+dx*dt;                          % 计算输出
    Y_DMC2(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC2(i,1)=u;                      % 存储到U_DMC
end
U_DMC2(SimuStepCount+1,1)=U_DMC2(SimuStepCount,1);
%% 绘图
n=size(Y_DMC,1);
subplot(2,1,1);
t=1:1:n;
plot(t,Y_DMC(1:n),'r');
hold on;
plot(t,Y_DMC2(1:n),'b--');
legend('DMC','改进移位矩阵')
title('y')
subplot(2,1,2);
plot(t,U_DMC(1:n),'r');
hold on;
plot(t,U_DMC2(1:n),'b--');
legend('DMC','改进移位矩阵')
title('u')

仿真结果

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

  要在自己的离散化中加入纯迟延过程，可以参照下面的代码，增加一个“缓冲序列”，在经过tao步后才输出没有纯迟延时当前步数的值。可以看到，优化参数后，提升了控制性能。

Matlab源码

clear;clc;                              % 清除缓存变量，清除命令行信息
dt=1;                                   % 采样周期
%################################ 存在纯迟延时用原始控制器参数 #####################################%
%% 被控对象参数
K=1;                                    % 开环增益
T=215;                                  % 惯性时间
%% DMC参数
N=2000;                                 % 建模时域
P=500;                                  % 预测时域
M=3;                                    % 控制时域
q=1;                                    % Q矩阵（预测输出与期望值误差权矩阵）中的一个权值
lamda=1;                                % λ(LAMDA)矩阵（控制增量权矩阵）中的一个权值
h=1;                                    % H矩阵（预测与实际输出误差校正系数矩阵）中的一个系数
%% 根据DMC参数计算DMC参数矩阵
Q=q*eye(P);                             % 生成Q矩阵
LAMDA=lamda*eye(M);                     % 生成LAMDA矩阵
H=h*ones(P,1);                          % 生成反馈校正H矩阵
S=zeros(P,P);                           % 生成移位矩阵S
for i=1:P-1
    S(i,i+1)=1;                         % 次对角线元素为1
end
S(P,P-1)=0; S(P,P)=1;                   % 右下角元素为1
%% 阶跃响应模型
num=K;                                  % 传递函数的分子部分
den=[T 1];                              % 传递函数的分母部分
tao=150;                                % 迟延时间
sys=tf(num,den,'iodelay',tao); 
a=step(sys,0:dt:N);                     % 参数依次为要施加阶跃输入的传递函数，起始时间：步长(周期)：截止时间

%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 仿真参数设置
SimuStepCount=2999;                     % 总仿真步数
sp=1;                                   % 设定值
SP=sp*ones(P,1);                        % 设定值向量
%% 开始仿真
Y_DMC=zeros(SimuStepCount+1,1);         % 保存每一步对象输出
U_DMC=zeros(SimuStepCount+1,1);         % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
YDealyArray=zeros(floor(tao/dt),1);     % 保存迟延的队列，长度为实际迟延除以采样周期
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    y=YDealyArray(tao,1);               % 实际输出取延迟后的Y,队列的最后一项   
    % 实现纯迟延过程
    dx=K/T*u+(-1/T)*YDealyArray(1,1);   % 按照没有纯迟延计算dx，队列中第一项YDealyArray(1,1)是没有纯迟延对象的前一时刻输出
    temp=YDealyArray(1,1)+dx*dt;        % 按照没有纯迟延计算输出
    for j=floor(tao/dt):-1:2            % 根据迟延逐个移动队列，实现经过tao时间后再输出没有延迟时计算输出    
        YDealyArray(j,1)=YDealyArray(j-1,1);
    end
    YDealyArray(1,1)=temp;
    % 存储数据
    Y_DMC(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC(i,1)=u;                      % 存储到U_DMC
end
U_DMC(SimuStepCount+1,1)=U_DMC(SimuStepCount,1);
Q_0=Q;
P_0=P;
N_0=N;
%% ################################ 存在纯迟延时用优化后控制器参数 #####################################%
tao=150;                                % 迟延时间
if tao>0 
   N=N_0+tao;                                 % 建模时域
   Q=[zeros(tao,tao) zeros(tao,P_0);zeros(P_0,tao) Q_0];
   P=P_0+tao;
end
sys=tf(num,den,'iodelay',tao); 
a=step(sys,0:dt:N);                     % 参数依次为要施加阶跃输入的传递函数，起始时间：步长(周期)：截止时间
%% 控制律求解参数的改进：如果带纯迟延需要改进P和Q
%% 根据DMC参数计算DMC参数矩阵
H=h*ones(P,1);                          % 生成反馈校正H矩阵
S=zeros(P,P);                           % 生成移位矩阵S
for i=1:P-1
    S(i,i+1)=1;                         % 次对角线元素为1
end
S(P,P-1)=0; S(P,P)=1;                   % 右下角元素为1
SP=sp*ones(P,1);                        % 设定值向量
%% 离线计算控制律前一部分的矩阵
A=zeros(P,M);                           % 生成动态矩阵A，A是P行M列
for i=1:M                               % 从1到M列
    for j=i:P                           % 从i到P行
        A(j,i)=a(j-i+1,1);
    end
end
C(1,1)=1;
for i=2:M
    C(1,i)=0;                           % 生成C矩阵
end
D=C*inv(A'*Q*A+LAMDA)*A'*Q;             % 得到控制律前一部分的矩阵
%% 开始仿真
Y_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步对象输出
U_DMC2=zeros(SimuStepCount+1,1);        % 保存每一步控制量
y_P1=zeros(P,1);                        % 当前时刻施加一个控制增量后未来P步预测值
y_P0=zeros(P,1);                        % 当前时刻没有施加控制增量未来P步预测值
y_COR=zeros(P,1);                       % 校正后的y_P1
YDealyArray=zeros(floor(tao/dt),1);     % 保存迟延的队列，长度为实际迟延除以采样周期
du=0;                                   % k时刻计算出的控制增量
u=0;                                    % k时刻计算出的控制量
y=0;                                    % 对象当前输出
for i=1:1:SimuStepCount
    % DMC计算部分
    e=y-y_P1(1,1);                      % 计算实际值与预测值误差
    y_COR=y_P1+H*e;                     % 误差反馈校正
    y_P0=S*y_COR;                       % 校正后预测值作为下一时刻未叠加控制增量的预测值
    du=D*(SP-y_P0);                     % 由控制律计算下一个时刻控制增量
    y_P1=y_P0+a(1:P,1)*du;              % 预测下一时刻施加du后输出
    u=du+u;                             % 这里的u在k时刻算出并施加到对象后，k+1时刻才能得到反馈效果    
    % u作用到对象后计算对象输出(如果是实际对象，程序只需要算出u即可，这里是仿真，对象输出也要计算)
    y=YDealyArray(tao,1);               % 实际输出取延迟后的Y,队列的最后一项   
    % 实现纯迟延过程
    dx=K/T*u+(-1/T)*YDealyArray(1,1);   % 按照没有纯迟延计算dx，队列中第一项YDealyArray(1,1)是没有纯迟延对象的前一时刻输出
    temp=YDealyArray(1,1)+dx*dt;        % 按照没有纯迟延计算输出
    for j=floor(tao/dt):-1:2            % 根据迟延逐个移动队列，实现经过tao时间后再输出没有延迟时计算输出    
        YDealyArray(j,1)=YDealyArray(j-1,1);
    end
    YDealyArray(1,1)=temp;
    % 存储数据
    Y_DMC2(i+1,1)=y;                    % 存储到Y_DMC
    U_DMC2(i,1)=u;                      % 存储到U_DMC
end
U_DMC2(SimuStepCount+1,1)=U_DMC2(SimuStepCount,1);
%% 绘图
n=size(Y_DMC,1);
subplot(2,1,1);
t=1:1:n;
plot(t,Y_DMC(1:n),'r');
hold on;
plot(t,Y_DMC2(1:n),'b--');
legend('存在纯迟延时用原始控制器参数','存在纯迟延时用改进控制器参数')
title('y')
subplot(2,1,2);
plot(t,U_DMC(1:n),'r');
hold on;
plot(t,U_DMC2(1:n),'b--');
legend('存在纯迟延时用原始控制器参数','存在纯迟延时用改进控制器参数')
title('u')

仿真结果