锅炉主蒸汽温度 — 预测控制（MPC/DMC）仿真方案
一、被控对象建模：为什么要这样建？
锅炉主蒸汽温度（过热器出口汽温）通过喷水减温来调节，其热工动态特性公认可用 "导前区 + 惰性区" 串联模型来描述：
减温水流量 Wj ──→ [导前区 G₁(s)] ──→ 减温器出口汽温 θ₂ ──→ [惰性区 G₂(s)] ──→ 过热器出口主汽温 θ₁
↑ 导前信号（副环反馈） ↑ 主信号（主环反馈）
| 分区 | 物理含义 | 动态特征 | 典型传递函数形式 |
|---|---|---|---|
| 导前区 G₁(s) | 减温器喷水 → 减温器出口温度 | 惯性小、响应快 | ( |
| 惰性区 G₂(s) | 减温器出口 → 过热器末端出口主汽温 | 大惯性 + 大纯迟延 | ( |
工程上常用的几种等价表达
形式① — 导前区+惰性区分开（便于串级控制设计）：
形式② — 合并为单回路等效对象（便于直接MPC设计）：
取自文献中对600MW级机组不同负荷工况的辨识结果，例如满负荷附近可用
或 等。
二、控制策略选择：为什么用预测控制？
主汽温的痛点是 大惯性 + 大纯迟延——常规PID只看"当前偏差"，相当于闭着眼追一个40~170秒后才反应的系统，必然超调大、调节时间长。
预测控制（MPC/DMC）的核心优势：
- 多步预测——把未来的延迟"看"进去，提前规划控制动作
- 滚动优化——每一步重新求解最优，容忍模型误差
- 天然处理约束——减温水流量有上下限，汽温有安全限值
推荐架构为 串级结构：
┌──────────────┐
T_sp ──►(+)──►│ 外环: MPC/DMC │──► u₂_sp ──►┌──────────┐
│ (主汽温回路) │ │ 内环: PID │──► W_j（减温水流量）
←── T₁ ←──│ │←── θ₂ ←───│ (导前回路)│
└──────────────┘ └──────────┘
- 内回路（PID）：快速稳定导前汽温 θ₂，压制二次扰动
- 外回路（MPC/DMC）：克服主通道的大延迟和大惯性，做设定值跟踪
三、完整 MATLAB 仿真代码
两套实现，你可以按需选用：
方案 A — 用 MATLAB MPC Toolbox
%% ============================================================
%% 锅炉主蒸汽温度 — MPC 预测控制仿真（基于MPC Toolbox）
%% ============================================================
clear; clc; close all;
%% ── 1. 被控对象模型 ──────────────────────────────────────────
% 等效主通道: G(s) = 1.2 * e^{-40s} / ((15s+1)*(40s+1)^2)
% 先建不带延迟的连续模型，再用 TransportDelay 属性挂延迟
s = tf('s');
G_cont = 1.2 / ((15*s + 1)*(40*s + 1)^2);
G_cont.IODelay = 40; % 纯迟延 40s
Ts = 2; % 采样时间 2s（工程上主汽温采样常用1~5s）
Gd = c2d(G_cont, Ts, 'zoh'); % 离散化供参考
fprintf('离散化模型（零极点形式）:\n'); Gd
%% ── 2. 构建 MPC 控制器 ───────────────────────────────────────
% 预测时域 Np = 20~40（需 > 延迟步数 = 40/2 = 20）
% 控制时域 Nu = 2~5
mpcobj = mpc(G_cont, Ts);
mpcobj.PredictionHorizon = 30; % Np
mpcobj.ControlHorizon = 3; % Nu
% 输出权重（主汽温跟踪重要性）
mpcobj.Weights.OutputVariables = 1.0;
% 控制增量权重（让控制动作平滑）
mpcobj.Weights.ManipulatedVariablesRate = 0.05;
% 控制量权重（可选）
mpcobj.Weights.ManipulatedVariables = 0.01;
% ── 约束 ────────────────────────────────────────────────────
% 减温水流量 0~100% (归一化) 或按实际物理量
mpcobj.Manipulated_variables.Min = -0.3; % 允许"反向降温"的增量限
mpcobj.Manipulated_variables.Max = 0.3;
% 汽温安全硬限（例：额定540℃，允许±10℃）
% mpcobj.Output_variables.Min = 530;
% mpcobj.Output_variables.Max = 550;
%% ── 3. 闭环仿真（用 Simulink-free 的 sim 方式：mpcmove） ───
Tsim = 800; % 总仿真 800s
Nstep = Tsim / Ts;
% 参考轨迹：阶跃型设定值
Tsp = 540; % 额定主汽温 ℃
ref = Tsp * ones(Nstep,1);
% （可选）中段变设定值看跟踪
% ref(201:end) = 545;
% 状态容器
xMPC = mpcstate(mpcobj); % MPC 状态
xPlant = zero(G_cont); % 连续对象状态（自己推进）
yHist = zeros(Nstep,1);
uHist = zeros(Nstep,1);
yCurr = 0; % 初态汽温 = 0（偏离量视角）或 540
for k = 1:Nstep
% —— 测量当前主汽温 ——
yMeas = yCurr + 0.2*randn; % 加少量测量噪声
% —— MPC 求解一步 ——
[u, Info] = mpcmove(mpcobj, xMPC, yMeas, ref(k));
% —— 施加到plant，推进一步 ——
% 用 lsimplus 方式手工积分（也可 c2d+filter 替代）
[~, ~, xPlant] = lsim(ss(G_cont), u*ones(1,1), [0 Ts], xPlant, 'foh');
% 读输出
C = tf(1.2/((15*s+1)*(40*s+1)^2));
% 更可靠的做法是直接用 dlsim on Gd:
if k==1
yCurr = 0; % 初始偏离归零（即实际540，建模用偏差量）
end
% 简便写法：直接用 feedback+lsim 闭环（见下方简化版）
uHist(k) = u;
end
%% ── 4. 用 Simulink-free 的高层方式重做（更干净的写法）────────
% 用 closed-loop sim 通过 d2c + step 展现——其实最简单是走 Simulink，
% 但如果你没有 Simulink，下面给一套"手写DMC"完全透明版。
disp(' MPC Toolbox 框架搭建完毕！');
disp(' 如需完整闭环曲线，请运行下方的「方案B：手写DMC」或搭Simulink框图。');
方案 B — 手写 DMC（动态矩阵控制） 完整闭环仿真
DMC 是 MPC 的一种，直接基于阶跃响应，非常适合大迟延热工对象，且不依赖任何工具箱。
%% ============================================================
%% 锅炉主蒸汽温度 — DMC（动态矩阵控制）手写完整仿真
%% ============================================================
clear; clc; close all;
%% ── 0. 参数 ─────────────────────────────────────────────────
Ts = 2; % 采样周期 [s]
Tsim = 1200; % 总仿真时长 [s]
N = Tsim/Ts; % 总步数
Tsp = 540; % 主汽温设定值 [℃]
theta_delay = 40; % 纯迟延 [s] ← 工程上可换成 63~170
N_delay = round(theta_delay/Ts); % 延迟步数
%% ── 1. 对象阶跃响应（解析法：对 G(s)=K/[(T1s+1)(T2s+1)^2] * e^{-θs}）
K = 1.2;
T1 = 15; T2 = 40;
% 留足够长的阶跃响应序列
N_a = N + N_delay + 200;
t_a = (0:N_a)*Ts;
% 三阶惯性阶跃响应解析式（无延迟部分）
% G_noDelay(s) = K / ((T1s+1)(T2s+1)^2)
% 阶跃响应: y(t) = K * [ 1 - A1*exp(-t/T1) - A2*exp(-t/T2) - A3*t.*exp(-t/T2) ]
% 部分分式展开求系数：
% 1/((T1s+1)(T2s+1)^2) = c1/(T1s+1) + c2/(T2s+1) + c3/(T2s+1)^2
c1 = 1/( (T2-T1)^2 );
c2 = -T2^2/( T1*(T2-T1)^2 );
c3 = -1/T1;
a_step = zeros(size(t_a));
for i = 1:length(t_a)
t = t_a(i);
if t <= 0
a_step(i) = 0;
else
term1 = c1 * exp(-t/T1);
term2 = c2 * exp(-t/T2);
term3 = c3*t .* exp(-t/T2);
a_step(i) = K * (1 - term1 - term2 - term3);
end
end
% 加纯迟延 → 真正的阶跃响应 a(k)
a = [zeros(1,N_delay), a_step(1:end-N_delay)];
a = a(1:N+60); % 截断到够用的长度
a = a(:); % 列向量
% 画图看看阶跃响应
figure('Name','对象阶跃响应','Color','w');
subplot(2,2,1); plot((0:length(a)-1)*Ts, a,'b','LineWidth',1.8);
grid on; xlabel('t [s]'); ylabel('幅值');
title('主汽温对象阶跃响应 a(k)');
%% ── 2. 构造动态矩阵 A ───────────────────────────────────────
Np = 35; % 预测时域（须 >> N_delay，建议 25~50）
Nu = 3; % 控制时域（2~5即可）
% 动态矩阵 A ∈ R^{Np×Nu}
A = zeros(Np, Nu);
for i = 1:Np
for j = 1:min(i, Nu)
if (i-j+1) <= length(a)
A(i,j) = a(i-j+1);
end
end
end
%% ── 3. DMC 参数 ─────────────────────────────────────────────
Q = eye(Np); % 输出跟踪权重（可改成 diag 调各步权重）
R = 0.08 * eye(Nu); % 控制增量权重（↑更平滑，↓更激进）
umin = -0.4; umax = 0.4; % 控制量(减温水变化率)约束
y_min = -Inf; y_max = Inf; % 输出约束（可自行加：±5℃）
%% ── 4. 仿真循环 ─────────────────────────────────────────────
% 状态
u_hist = zeros(N,1); % 实际控制量（减温水阀门/流量 归一化）
du_hist= zeros(N,1); % 控制增量 Δu
y_hist = zeros(N,1); % 测量输出（主汽温 偏差量，即 y = T - Tsp_0）
y_model= zeros(N,1); % 模型预测输出
d_hist = zeros(N,1); % 可加外扰（蒸汽流量扰动等）
% 初始条件：对象内部用移位寄存器存历史Δu
buffer_len = length(a);
deltaU_buf = zeros(buffer_len, 1); % 足够长的历史Δu缓冲区
% 误差校正（反馈校正系数）
h = ones(Np,1); % 最简单的：全1向量（可换成 [1,0.8,0.6,...]）
ref_profile = zeros(N,1); % 设定值轮廓（偏差量视角：ref=0 表示跟踪540不变）
ref_profile(:) = 0; % 即目标=540℃
% 如果想看设定值阶跃：
% ref_profile(1:round(50/Ts)) = 0; % 已在540，或从中段阶跃
% ---- 主循环 ----
for k = 1:N
%% 4.1 测量（含噪声）
ym = y_hist(k);
%% 4.2 预测：无未来控制时的"自由响应"
y0 = zeros(Np,1);
for i = 1:Np
% y0(i) = Σ_{j=1}^{Nbuf} a(i+j) * Δu_buf(j) 即已有历史Δu推出来的预测
acc = 0;
for j = 1:buffer_len-i
acc = acc + a(i+j) * deltaU_buf(j);
end
y0(i) = acc;
end
%% 4.3 误差校正
e = ym - y0(1); % 模型-实际 偏差
y_corr = y0 + h * e; % 校正后的预测基准
%% 4.4 求解 ΔU* = argmin ||Q(y_corr+AΔU - ref(k:k+Np-1))||² + ΔU'RΔU
% 无约束解析解：
H = A'*Q*A + R;
g = A'*Q*(ref_profile(k:min(k+Np-1,N)) - y_corr(1:length(ref_profile(k:min(k+Np-1,N)))));
% 对齐维度（ref可能尾部不够长）
ref_vec = ref_profile(k:min(k+Np-1,N));
g = A(1:length(ref_vec),:)' * Q(1:length(ref_vec),1:length(ref_vec)) ...
* (ref_vec - y_corr(1:length(ref_vec)));
du_opt = pinv(H) * g; % Δu*
% 简化写法（更鲁棒）：
du_opt = (A(1:length(ref_vec),:)'*Q(1:length(ref_vec),1:length(ref_vec))*A(1:length(ref_vec),:) + R) \ ...
(A(1:length(ref_vec),:)'*Q(1:length(ref_vec),1:length(ref_vec))*(ref_vec - y_corr(1:length(ref_vec))));
%% 4.5 只取第一个控制增量（滚动）
du_k = du_opt(1);
% 限幅
if du_k > umax, du_k = umax; end
if du_k < umin, du_k = umin; end
%% 4.6 更新控制量
u_k = u_hist(max(k-1,1)) + du_k;
u_hist(k) = u_k;
du_hist(k)= du_k;
%% 4.7 推对象状态（卷积更新）
% shift buffer
deltaU_buf = [du_k; deltaU_buf(1:end-1)];
% 计算当前模型输出
y_pred_k = 0;
for j = 1:buffer_len
y_pred_k = y_pred_k + a(j) * deltaU_buf(j);
end
y_model(k) = y_pred_k;
%% 4.8 推进到"真实"plant（这里用模型+扰动=真值，所以 y_true≈y_model+noise）
noise = 0.15*randn; % 测量噪声
d_k = 0; % 外扰（可在某段注入，见下方）
if k>=150 && k<=200
d_k = -0.8; % 模拟一次蒸汽流量突增扰动
end
y_true = y_pred_k + d_k + noise;
% 准备下一拍（这里 y_hist(k+1) 的真实值需要plant数据——
% 因为我们拿模型当plant，所以：）
if k < N
y_hist(k+1) = y_true;
end
end
%% ── 5. 绘图 ─────────────────────────────────────────────────
t_axis = (0:N-1)*Ts;
figure('Name','DMC主汽温控制仿真','Color','w','Position',[200 100 1100 700]);
% (1) 主汽温跟踪
subplot(3,1,1);
plot(t_axis, y_hist+Tsp, 'b','LineWidth',2); hold on;
plot(t_axis, Tsp*ones(size(t_axis)), 'r--','LineWidth',1.2);
plot(t_axis, (Tsp+5)*ones(size(t_axis)), 'r:','LineWidth',0.8);
plot(t_axis, (Tsp-5)*ones(size(t_axis)), 'r:','LineWidth',0.8);
grid on;
ylabel('主汽温 T_1 [℃]');
title(sprintf('DMC控制 — 主蒸汽温度跟踪 (θ_d=%.0fs, Np=%d, Nu=%d)',theta_delay,Np,Nu));
legend(sprintf('T_1(t) [基准%.0f℃]',Tsp),'Sp=540℃','±5℃安全限','Location','best');
% (2) 控制量：减温水相对变化
subplot(3,1,2);
plot(t_axis, u_hist, 'k','LineWidth',1.5); grid on;
ylabel('u (减温水 Δ%)');
title('控制量（减温水调节动作）');
% (3) 控制增量
subplot(3,1,3);
stairs(t_axis, du_hist, 'm','LineWidth',1.2); grid on;
xlabel('t [s]'); ylabel('Δu');
title('控制增量 Δu(k)');
sgtitle('锅炉主蒸汽温度 — DMC预测控制仿真结果','FontSize',12,'FontWeight','bold');
参考代码 对锅炉主蒸汽温度控制进行预测控制仿真 www.youwenfan.com/contentcnv/81325.html
四、关键参数怎么调
| 参数 | 含义 | 调法 |
|---|---|---|
| Np（预测时域） | "往前看多远" | 必须 > 延迟步数 Nd = τ/Ts（即 Np ≥ 1.5~2×Nd），太小则等于没预测，太大则计算慢 |
| Nu（控制时域） | 优化几个Δu | 2~5 够了，更大→更柔但更迟钝 |
| R（控制权重） | 惩罚控制动作 | R↑ → 阀门动作更稳但响应变慢；R↓ → 激进跟踪 |
| 约束 umin/umax | 减温水阀门物理限 | 按执行器量程设（例 ±30%/step） |
| 误差校正向量 h | 反馈修正模型误差 | 最简单 h=ones，更精细可用指数衰减 h=[1,λ,λ²,...] (λ≈0.7~0.9) |
典型组合（τ≈40s, Ts=2s ⇒ Nd≈20）：Np=30~40, Nu=2~4, R=0.05~0.2
