目录
- 前文链接
- 基础DMC算法
- Matlab源码
- 代码解析
- 仿真结果展示
- 参数对性能的影响
- 参数P对性能影响
- Matlab源码
- 仿真结果
- 参数M对性能影响
- Matlab源码
- 仿真结果
- 参数q对性能影响
- Matlab源码
- 仿真结果
- 参数lamda对性能影响
- Matlab源码
- 仿真结果
- 讨论
- 算法改进效果验证
- 阶梯式动态矩阵控制与基础DMC对比
- Matlab源码
- 仿真结果
- 反馈校正环节截断误差改进
- Matlab源码
- 仿真结果
- 对象存在纯迟延时控制器参数优化
- Matlab源码
- 仿真结果
前文链接
经典预测控制算法：动态矩阵控制(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-1S(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:MC(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC(i+1,1)=y; % 存储到Y_DMCU_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循环中的内容）。
上面代码中设置参数的部分我们直接跳过，变量名与之前文章中的符号都是对应的。
被控对象为一个不带纯迟延环节的一阶惯性对象，用传递函数来表示是 K T s + 1 \frac{K}{Ts+1} Ts+1K，matlab中的tf()函数可以将分子和分母表示的多项式（从最高次项一直到常数也就是0次项）转换为传递函数形式，例如下面代码会生成一个 A s + 1 B s 2 + C s + 1 \frac{As+1}{Bs^2+Cs+1} Bs2+Cs+1As+1形式的传递函数：
num=[A 1];
den=[B C 1];
sys=tf(num,den);
step()函数可以生成DMC需要的阶跃响应模型，我们在参数中提供之前生成的传递函数，以及起始结束时间和步长，函数返回的结果就是单位阶跃输入下的阶跃响应模型 a \bm{a} a。
得到 a \bm{a} a之后就可以计算动态矩阵 A \bm{A} A了， A \bm{A} A再结合已经确定的参数矩阵，就可以计算出控制律中前半部分。这一点之前文章中也提到过，控制律不需要每一步完全重新计算，只有后半部分需要每一步重新计算，这里固定的前半部分我用变量D表示。
接下来终于到了仿真中关键的循环部分，Y_DMC和U_DMC用来存储每一步的输出和控制量，Y_P0为 y ~ P 0 \bm{\widetilde{y}_{P0}} y P0，Y_P1为 y ~ P 1 \bm{\widetilde{y}_{P1}} y P1，Y_COR为 y ~ COR \bm{\widetilde{y}_\text{COR}} y COR。上篇中提到过，为了让程序第一个时刻可以走下去，第一个时刻先将这些预测值，控制量等全部初始化为0。有人会较真，预测值必须初始化为0吗？澄清一下，预测值初始化为1甚至1000000都可以，由于这里仿真的时候对象输出是从0开始变化的，所以初始化为0保证了第一个时刻预测值和实际值的误差 e e e不会太大，导致不必要的控制过程波动，因为第一步的 e 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 k k时刻，先计算 e e e，然后计算修正后的预测值 y ~ COR \bm{\widetilde{y}_\text{COR}} y COR，移位后得到对 k + 1 k+1 k+1时刻未施加控制增量的预测值 y ~ P 0 \bm{\widetilde{y}_{P0}} y P0，求解控制律得到du,叠加在之前的控制量u之上，得到 k k k时刻的控制量施加给对象，同时计算施加du之后的预测值，再次强调一下，标蓝的文字都是在 k k k时刻的计算过程，得到的是 k k k时刻的控制量，但是 k k k时刻控制量作用到对象之后， k + 1 k+1 k+1时刻才能收到对象反馈值，时序非常关键，哪个时刻发生了哪些事一定不要弄混！y_P1=y_P0+a(1:P,1)*du;这一行倒是无所谓，放在 k k k时刻计算出控制律后，或者放在求 e e e之前都可以，只不过放在计算控制律之后表示求y_P1发生在 k k k时刻，放在求 e e e之前表示求y_P1发生在 k + 1 k+1 k+1时刻。放在求 e e 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 k+1 k+1时刻的对象反馈，u是 k k 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-1S(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:MC(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC(i+1,1)=y; % 存储到Y_DMCU_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-1S(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:MC(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC2(i+1,1)=y; % 存储到Y_DMCU_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-1S(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:MC(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC3(i+1,1)=y; % 存储到Y_DMCU_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-1S(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:MC(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC(i+1,1)=y; % 存储到Y_DMCU_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:MC(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC2(i+1,1)=y; % 存储到Y_DMCU_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:MC(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC3(i+1,1)=y; % 存储到Y_DMCU_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-1S(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:MC(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC(i+1,1)=y; % 存储到Y_DMCU_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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC2(i+1,1)=y; % 存储到Y_DMCU_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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC3(i+1,1)=y; % 存储到Y_DMCU_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-1S(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:MC(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC(i+1,1)=y; % 存储到Y_DMCU_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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC2(i+1,1)=y; % 存储到Y_DMCU_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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC3(i+1,1)=y; % 存储到Y_DMCU_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增大，首先会使得动态矩阵 A \bm{A} A的行数增加， 导致控制律前半部分的值增大， 但是参数P增大，又会使得控制律后半部分设定值SP与预测值y_P0的差构成向量末尾出现更多的0（可以想象到，P足够大时，设定值与预测值到后面必然会相等，导致继续增大P会使得SP-y_P0末尾出现越来越多的0），所以从控制律公式分析，不好确定P对控制效果的影响。于是我们从仿真结果来看，P的增加导致控制量给的更稳了，对于被控对象来说，这就意味着到达设定值的过程变慢了。
参数M增大，会使得动态矩阵 A \bm{A} A的列数增加， 导致控制律前半部分的值增大， 而M对控制律后半部分没有影响，那么会使得根据控制律得到的控制增量增大，也就是说，从控制律公式来看，M增大应该会使得控制量变化更明显，相应地被控对象到达设定值时间会加快。从仿真结果来看，和我们的推测是一致的，M增大时，一开始的控制量会给得更猛。
参数q和参数lamda从它们的定义中就知道，q越大，被控对象到达设定值会越快，控制量变化会更剧烈，lamda越大，控制量变化会更受限制，使得控制过程更稳定，被控对象到达设定值变慢，仿真结果也是完全对应的。但是还有一个小细节，如果把我的代码里q改成2，3或者lamda改成2，3，会看到几条几乎重合的曲线，实际上你把曲线放大放大再放大，可能要到小数点后四五位，确实也能看到符合上述结论，然后你就明白我上面的q和lamda为什么要设置这么大/这么小了，这提醒了我们，q和lamda确实可以影响控制效果，但是要能够真正“明显”地看到影响，需要先对动态矩阵 A \bm{A} 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-1S(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:MC(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC(i+1,1)=y; % 存储到Y_DMCU_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:Mlam(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC2(i+1,1)=y; % 存储到Y_DMCU_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-1S(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:MC(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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC(i+1,1)=y; % 存储到Y_DMCU_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; % 计算dxy=y+dx*dt; % 计算输出Y_DMC2(i+1,1)=y; % 存储到Y_DMCU_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-1S(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:MC(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);endYDealyArray(1,1)=temp;% 存储数据Y_DMC(i+1,1)=y; % 存储到Y_DMCU_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-1S(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:MC(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);endYDealyArray(1,1)=temp;% 存储数据Y_DMC2(i+1,1)=y; % 存储到Y_DMCU_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')
仿真结果
