考虑需求响应的电热综合能源系统两阶段优化调度及Matlab实现

搞电-热综合能源系统优化调度的人，一开始都很容易陷入一个误区：把模型堆得特别复杂，电、热、气、冷全都塞进去，约束写了几百行，结果求解器跑半个小时不出结果，最后只好把论文里写的“日前-日内两阶段”从PPT里删掉，因为代码根本没搭出两阶段框架。

我刚开始做这个方向的时候也一样，走了不少弯路。直到后来想明白一个问题：所谓“多时间尺度优化调度”，本质上不是把模型做大，而是把决策分层——哪些变量需要在昨天就定下来，哪些变量可以留到今天见招拆招。这个分层逻辑理顺了，Matlab代码的架构自然就清晰了。这篇就基于“考虑需求响应的电-热综合能源系统两阶段日前日内多时间尺度优化调度策略”这个题，把我从建模到代码实现的完整思路、以及调试过程中踩过的坑一次性说清楚，给正在做类似方向的同学一个可以直接参考的框架。

1. 先搞清楚优化对象：这个系统里到底有哪些设备在参与调度

1.1 电-热耦合系统的典型物理架构

先说为什么要“电-热综合”。单独做电力系统调度，或者单独做热力系统调度，模型都能简单很多，但现实中的设备恰恰是跨能流耦合的。最常见的耦合设备就是热电联产机组（CHP），它烧天然气（或煤）产生的高温蒸汽先用来发电，发电后的余热进入热网供热。这种“以热定电”或“以电定热”的运行特性，让电和热两套系统互为约束——如果只看电网，你没法决定CHP的电出力；如果只看热网，你也没法决定CHP的供热量，因为两者共用同一个蒸汽量。

典型的电-热综合能源系统一般包含这么几类设备：

设备类型	作用	调度中的关键特性
CHP机组	联产电和热，系统的核心	电出力与热出力相互耦合，需满足可行运行域约束；有最小开停机时间和爬坡约束
电锅炉	将富余电力转化为热能	效率通常在0.95~0.98之间，启停灵活，秒级响应，是弃风消纳的利器
储热罐	存储热网中的热能	有储/放热功率上限、储热容量上限，充放热时存在热损失
常规纯凝机组	只发电不供热	爬坡约束、出力上下限、最小开停机时间
风电/光伏	不可控新能源	出力具有随机性和波动性，日前预测误差较大
电负荷/热负荷	用户侧需求	部分电负荷具备可调节能力，即需求响应资源

这个结构里最值得关注的是CHP。国内很多文献喜欢把CHP简化成“热电比恒定”的定比模型，就是P_CHP = k * H_CHP，这种模型写起来简单，但实际运行中CHP的热电比是可以在一定范围内变化的，它的运行可行域是下图这样一个多边形区域（AB-BC-CD-DE-EF围起来的区域），包含背压工况和抽汽工况。用可行域约束替代定比约束，更贴近实际，求解难度也并没有增加太多，只是需要在代码里多写几个线性不等式。

1.2 需求响应在这个模型里的角色

需求响应（Demand Response, DR）在这个系统里不是一个“可选项”，而是参与优化调度的虚拟资源。

在电-热系统里，需求响应有两种参与方式。第一种是价格型需求响应，用户根据电价信号的引导调整用电行为——比如电价高的时段减少用电，电价低的时段增加用电。这种DR在模型里表现为价格弹性矩阵，电负荷的调整量是电价变化量的函数，通常写成：

ΔL_t = E * L_0_t * (Δπ_t / π_0_t)

其中E是自弹性系数和交叉弹性系数组成的矩阵。自弹性系数一般为负（电价升高用户减少用电），交叉弹性系数一般为正（其他时段电价升高，用户会把负荷转移到当前时段）。这个矩阵的估计可以从用户调查数据或者历史负荷数据回归得到，实际算例中常用 -0.2 作为自弹性系数、0.03 作为交叉弹性系数的经验值。

第二种是激励型需求响应，直接与用户签订合同，约定在系统需要时可以削减或转移一部分负荷，系统给予用户补偿。在模型里表现为可中断负荷或者可平移负荷约束：

0 ≤ ΔL_DR_t ≤ ΔL_DR_max_t

即每个时段的DR削减量有上限。目标函数里相应要加入DR补偿成本项。

我认为真正有用的需求响应建模，是把它当作“等效储能”。用户削减了一段负荷，相当于系统从负荷侧获得了一个“负功率”，这个负功率和储热罐的充放热在数学上有着相似的约束形态——都有容量上限、都有时段耦合关系。这样理解，写代码的时候就不容易乱。

2. 两阶段调度框架：怎么划分“日前”和“日内”的决策边界

2.1 为什么需要两阶段：预测精度决定调度精度

如果风电、光伏出力预测是完美的，负荷预测是完美的，那就不需要两阶段了，直接做一次日前优化，把未来24小时每个设备每个时刻的出力定下来就行。但现实中的预测误差——尤其是风电功率预测，日前误差可达15%~25%，日内超短期预测能降到5%以内——导致日前制定的计划根本没法直接执行。

这就是两阶段调度（Two-Stage Scheduling）存在的根本原因。日前阶段基于日前预测数据，确定一些一旦定下来就不方便频繁调整的决策变量；日内阶段基于最新的超短期预测数据和实时运行状态，在日前计划的基础上进行修正，追踪真实运行工况。

一句话概括：日前计划负责“保底”，日内优化负责“纠偏”。

2.2 第一阶段（日前）：慢变量与长周期决策

日前调度的主时间尺度一般是1小时，决策窗口是未来24小时，即T=24个时段。在这个阶段要确定的决策变量包括：

• 机组的启停状态（0-1整数变量，启动/停机后由于最小开停机时间约束，短时间内不能频繁切换）
• CHP机组的电出力/热出力计划值
• 电锅炉的启停与计划出力
• 储热罐的充放热计划
• 需求响应的基线削减计划（即提前通知用户“明天下午3点到5点准备削减负荷”）
• 与大电网的购售电功率计划

这些变量有一个共同特征——都是需要提前一天“承诺”的。尤其是机组启停和DR签约量，一个是物理设备层面的承诺（机组开起来不能立刻停），一个是商业合同层面的承诺（答应用户削减负荷不能随意改），都不适合在日内频繁调整。

日前阶段的优化目标通常是系统总运行成本最小化，包括：

• 机组燃料成本（CHP机组的燃料成本是电出力和热出力的函数）
• 启停成本
• 购电成本（从外部电网购电）
• 需求响应补偿成本

2.3 第二阶段（日内）：快变量与滚动修正

日内调度的时间尺度一般为15分钟或5分钟，决策窗口可以是未来4小时（滚动时域优化）或与日前一样覆盖未来24小时但步长更细。

日内阶段的决策变量包括：

• 各机组出力的修正量（在日前计划基础上调整）
• 电锅炉出力的修正量
• 储热罐充放热功率的修正量
• 实际DR削减量的调整（但通常只在日前签约量上下一定范围内调整）

日内阶段的目标是最小化运行成本的修正量，同时要保证跟踪日前计划的偏差尽量小。在代码实现层面，这是通过在目标函数中加入偏差惩罚项实现的：

min Σ(燃料成本修正 + 购电成本修正 + DR补偿修正) + λ * Σ(实际出力 - 日前计划的偏差)²

λ是惩罚系数，取值大小决定了日内优化对日前计划的“忠诚度”。λ设得太大，日内优化基本不敢偏离日前计划，等于没做；λ设得太小，日前计划形同虚设，机组启停状态可能无法满足实际负荷需求。实际调参经验是：λ的量级应该与边际成本系数相当，这样日内优化既有修正能力又不会任意偏离。

2.4 两阶段之间的衔接逻辑

两阶段之间的数据传递是代码架构里最容易被忽略的部分。

我的做法是：日前阶段结束后，把以下结果保存到一个结构体schedule_day_ahead中——机组启停状态u_onoff、CHP电出力P_chp_da、CHP热出力H_chp_da、电锅炉出力P_eb_da、储热罐充放热功率H_st_da（正为放热，负为充热）、DR削减量L_dr_da、购电功率P_buy_da。

日内滚动优化时，每一轮的输入包括：

1. 上一轮日前计划中对应时段的值
2. 实时更新的超短期预测（风电出力、电负荷、热负荷）
3. 当前设备实际运行状态

然后输出修正量ΔP_chp、ΔH_chp、ΔP_eb等，实际执行值 = 日前计划值 + 修正量。

P_chp_actual_t = P_chp_da_t + ΔP_chp_t

这个衔接关系写清楚后，代码的模块边界就非常清晰了。

3. Matlab代码实现的核心架构：从数据输入到结果输出

3.1 代码模块划分与主要函数

Matlab实现这套调度策略，我推荐按下面的模块划分组织代码文件，各模块职责单一，方便调试和复用：

main_optimization.m % 主程序，依次调用各模块 load_data.m % 读取负荷、风电、电价等基础数据 build_system_params.m % 定义系统设备参数 build_day_ahead_model.m % 构建日前优化模型（YALMIP变量、约束、目标） solve_day_ahead.m % 求解日前问题 build_intraday_model.m % 构建日内滚动优化模型 solve_intraday.m % 求解日内问题 plot_results.m % 结果可视化

我个人习惯把所有设备参数集中放在build_system_params.m里，用一个结构体params统一管理。这样后面做敏感性分析（比如改储热罐容量、改DR弹性系数）时，只需要修改这个文件中的参数值，不需要在模型文件里到处找。

3.2 用电热综合能源系统专用建模框架还是YALMIP

Matlab里做优化调度，主流的建模方式有三种：

• 直接用Matlab优化工具箱：linprog（线性规划）、intlinprog（混合整数线性规划）、fmincon（非线性规划）。优点是零依赖、不需要额外安装；缺点是把庞大的调度模型写成矩阵形式很痛苦，尤其是热网节点网络约束多的情况下，矩阵索引容易写错。
• YALMIP + 外部求解器：YALMIP是一个Matlab建模工具箱，可以用接近数学表达式的形式描述变量和约束，然后调用CPLEX、Gurobi、MOSEK等商业求解器。代码可读性远高于手写矩阵，而且求解效率也不错。这是我最推荐的方式。
• MATPOWER等专用工具：主要针对电力系统潮流计算，虽然也能做优化，但热网部分需要自己扩展，不如YALMIP灵活。

我见过一些同学坚持用intlinprog手写矩阵，最后因为约束索引错了一位，导致结果出现“负热产出”这种根本不可能出现的数值而查了半天。我的建议是：只要模型的约束超过20行，就直接上YALMIP。

YALMIP的安装很简单，把下载的文件夹解压后加入Matlab路径即可。求解器方面，CPLEX和Gurobi是主流选择，如果实验室没有商业求解器授权，可以用SCIP或者CBC代替，Matlab自带的intlinprog也可以作为兜底选项（但求解大规模MILP会慢一些）。

3.3 日前优化模型的核心代码骨架

下面给出日前优化模型的核心代码骨架，基于YALMIP语法。这里展示的是经过我实际调试通过的简化版本，完整版主要是增加更多设备约束和更精细的目标函数项。

% 定义决策变量 u_chp = binvar(1, T); % CHP启停状态 P_chp = sdpvar(1, T); % CHP电出力 H_chp = sdpvar(1, T); % CHP热出力 P_eb = sdpvar(1, T); % 电锅炉输入电功率 H_st = sdpvar(1, T); % 储热罐放热功率（正放热，负充热） SOC_st = sdpvar(1, T); % 储热罐储热状态 L_dr = sdpvar(1, T); % DR削减量 P_buy = sdpvar(1, T); % 从电网购电功率 Constraints = []; % 电功率平衡约束（不考虑网络潮流，简化为单节点） Constraints = [Constraints, P_chp + P_wind - P_eb + P_buy == P_load - L_dr]; % 热功率平衡约束 Constraints = [Constraints, H_chp + H_eb + H_st == H_load]; % 电锅炉产热（效率eta_eb） H_eb = params.eta_eb * P_eb; Constraints = [Constraints, 0 <= P_eb <= params.P_eb_max]; % CHP可行运行域约束（以两段式线性近似为例） % 简化形式：H_chp >= alpha_1 * P_chp + beta_1 % H_chp <= alpha_2 * P_chp + beta_2 Constraints = [Constraints, H_chp >= params.alpha_1 * P_chp + params.beta_1 * u_chp]; Constraints = [Constraints, H_chp <= params.alpha_2 * P_chp + params.beta_2 * u_chp]; % 电出力上下限 Constraints = [Constraints, params.P_chp_min * u_chp <= P_chp <= params.P_chp_max * u_chp]; % 储热罐动态约束 Constraints = [Constraints, SOC_st(:, t+1) == SOC_st(:, t) * (1 - params.eta_loss) - H_st(:, t)]; Constraints = [Constraints, -params.H_st_max <= H_st <= params.H_st_max]; Constraints = [Constraints, params.SOC_min <= SOC_st <= params.SOC_max]; Constraints = [Constraints, SOC_st(:, T+1) == SOC_st(:, 1)]; % 周期约束 % DR削减量约束 Constraints = [Constraints, 0 <= L_dr <= params.L_dr_max]; % 目标函数：总运行成本最小化 Objective = 0; for t = 1:T % 燃料成本（CHP，二次函数近似为分段线性） Objective = Objective + (params.a * P_chp(t)^2 + params.b * P_chp(t) + params.c) * u_chp(t); % 购电成本 Objective = Objective + params.price_buy(t) * P_buy(t); % DR补偿成本 Objective = Objective + params.price_dr * L_dr(t); end % 求解 ops = sdpsettings('solver', 'cplex', 'verbose', 2); optimize(Constraints, Objective, ops); % 提取结果 P_chp_da = value(P_chp); H_chp_da = value(H_chp); % ...

如果你用的是linprog，还需要自己把二次项处理掉，但既然用了YALMIP，二次目标函数是可以直接交给CPLEX处理的（CPLEX支持MIQP），省很多事。

3.4 日内滚动优化的关键代码逻辑

日内部分我采用的是滚动时域优化（Rolling Horizon）。所谓滚动，就是每过15分钟重新求解一次未来4小时的优化问题，并且只执行第一个时段的结果。Matlab里用循环模拟：

%% 日内滚动优化循环 for k = 1:T_intraday_total % 每个滚动时刻 % k时刻的优化窗口为[k, k+H_intraday-1] t_start = k; t_end = min(k + H_intraday - 1, T_intraday_total); % 构建日内优化问题 [Constraints, Objective] = build_intraday_model(t_start, t_end, ...); % 求解 optimize(Constraints, Objective, ops); % 只执行第一个时段结果 P_chp_execute(k) = value(P_chp_intraday(1)); % 更新下一轮初始状态（如储热罐SOC） SOC_current = value(SOC_intraday(1)); % 读取新一轮超短期预测数据 [P_wind_hat, P_load_hat, H_load_hat] = get_forecast(k+1, k+H_intraday); end

这里有一个非常重要的实现细节：日内优化的决策时间段索引必须与日前保持一致。也就是说，如果日前计划是1小时一个时段，日内优化是15分钟一个时段，那么在日内修正日前计划时，需要把日前计划插值或者对应到日内的时间粒度上。我通常的做法是：日前按24时段（1小时），日内按96时段（15分钟），日前计划在每个15分钟时段内取相同值，即P_chp_da_15min(k) = P_chp_da(ceil(k/4))。这样对应关系清楚，不容易出错。

4. 建模与求解中的关键细节：可行域、爬坡约束与非线性处理

4.1 CHP可行运行域的正确描述方式

很多初学者的代码里，CHP的约束就是P_min ≤ P ≤ P_max和H_min ≤ H ≤ H_max两条边界独立约束，这就完全忽略了热电耦合关系。实际上，CHP的运行范围在(P, H)平面上是一个凸多边形，由若干个顶点定义。描述方式是用一组线性不等式：

A_region * [P; H] ≤ b_region

以常见的抽汽式CHP为例，可行运行域的顶点通常有4~6个。在Matlab代码里，最稳妥的做法是直接把顶点坐标写在参数文件里，然后用polytope或者简单的手写线性不等式组描述。

% CHP可行运行域顶点坐标（示例） vertices = [0, 0; 100, 0; 100, 50; 60, 80; 0, 80]; % 求半空间表示 A*x <= b [A_region, b_region] = vert2lcon(vertices); % 需要鲁棒控制工具箱，或者自己写凸包算法

如果你没有vert2lcon函数，可以自己写一个基于凸包算法的函数，或者简单点，直接把顶点代入不等式推导出A、b矩阵。这些不等式中的每一行都表示可行域的一个边界条件，写入约束时注意乘以机组的0-1启停变量，确保停机时可行域收缩为原点：

A_region * [P_chp; H_chp] ≤ b_region * u_chp

4.2 爬坡约束的时间粒度统一

爬坡约束是另一处容易踩坑的地方。CHP机组的爬坡约束写作：

-P_ramp_down ≤ P_chp(t) - P_chp(t-1) ≤ P_ramp_up

这里爬坡速率单位是MW/h，如果你的日前模型时间步长是1小时，那么不等式右边直接用爬坡速率即可。但如果日内模型步长是15分钟，就需要写：

-P_ramp_down / 4 ≤ P_intraday(k) - P_intraday(k-1) ≤ P_ramp_up / 4

很多同学在这里忘记除以4，导致日内模型爬坡约束过紧，求解结果出现机组无法跟踪预测负荷的情况。同样的逻辑适用于储热罐、电锅炉的功率变化约束。

4.3 非线性项的处理：二次成本函数与双线性约束

CHP燃料成本一般是二次函数，如C = a * P² + b * P + c * H + d。在MILP里直接处理二次函数有两条路：一是用CPLEX/Gurobi直接求解MIQP（如果你的模型规模不大，这是最简单的）；二是做分段线性近似（Piecewise Linear Approximation），把二次曲线切成若干段，每段用线性函数近似，加上SOS2约束（Special Ordered Set Type 2）或用二进制变量选择线性段。

我的建议是：先用MIQP跑通流程，如果求解速度不满足要求，再换成分段线性近似。YALMIP对这两种方式都支持，切换成本不高。用sdpvar建二次目标后直接求解，如果换成PWL，只需要把目标函数里的二次项替换成sdpvar+binvar的组合。

还有一个常见的双线性约束是：电锅炉的耗电量和产热量之间的关系，如果效率恒定，是线性关系H = η * P，没问题。但如果考虑效率随负荷率变化，就变成了双线性约束H = η(P) * P，这时候要么拟合线性关系，要么分段线性化。实际算例中，电锅炉效率随负荷变化不大，用恒定效率足够。

4.4 储热罐的动态建模与周期约束

储热罐的模型是典型的离散时间状态转移方程：

SOC(t+1) = (1 - η_loss) * SOC(t) - H_st(t) * Δt

这里H_st(t)是放热功率，正为放热，负为充热。值得注意的是，如果同时考虑充热效率和放热效率，约束形式会变为分段形式：

SOC(t+1) = (1 - η_loss) * SOC(t) - H_st_discharge(t) / η_dis + H_st_charge(t) * η_ch

这就引入了一对互补变量H_st_discharge和H_st_charge，加上两者不同时为正的约束。YALMIP里可以用implies写逻辑约束，或者通过引入二进制变量来处理。如果只是做简化研究，用单效率形式SOC(t+1) = (1-η_loss)*SOC(t) - η_st * H_st(t)也可行。

周期约束SOC(T+1) = SOC(1)也是一个容易忽略但很重要的约束，它保证调度周期结束时储热罐回到初始状态，使方案具备可持续性。如果不加这个约束，优化器倾向于在周期结束时把储热罐的热量全部放完，获取最大收益。

5. 实际调试过程中踩过的坑与解决手段

5.1 模型不可行：不是约束写错了，就是平衡方程错了

做调度优化，最常见的异常就是求解器返回infeasible problem。如果你是第一次跑通代码，出现不可行不要太惊讶，这是常态。

我的排查顺序是这样的：

1. 检查功率平衡约束：把各自的出力上下限求和，看是否覆盖负荷变化范围。比如夜间热负荷高、风电出力大而电负荷低，如果电锅炉最大耗电功率不足以消纳全部风电，电功率平衡约束就可能无法满足。这时候需要放宽约束（允许弃风）或者提高电锅炉容量。

2. 检查储热罐SOC范围：如果SOC初始值设定了，而周期约束要求SOC(T+1)=SOC(1)，那么在储热罐容量很小的情况下可能找不到可行解。解决办法是引入松弛变量，允许周期约束有小偏差，并在目标函数中加惩罚项。

3. 检查0-1变量与连续变量的耦合：比如启停状态为0时出力必须为0，这个逻辑如果漏写了某个变量，求解器可能给出“机组停机但仍有出力”的荒谬结果，甚至导致平衡约束崩掉。

YALMIP有一个好用的小技巧，在求解前用check命令检查每一个约束的残差：

optimize(Constraints, Objective, ops); check(Constraints);

它会把每个约束的边界值打出来，残差不为0的约束就是问题所在。

5.2 求解时间过长：整数变量太多是罪魁祸首

日前24时段的MILP模型，如果有2台CHP机组，每台机组每时段1个0-1变量，那就是48个二进制变量。这个规模对CPLEX来说是小意思，几秒钟就能求出来。

但如果你的模型规模变大——比如细化到96时段、机组数量增加到5台、再加上储热罐的充放二进制变量、DR的分段线性二进制变量——整数变量数量会轻松突破500个，求解时间可能从几秒飙升到十几分钟甚至几小时。

我的经验是：

• 优先优化模型结构而非依赖求解器：检查是否存在冗余约束，是否存在可以合并的变量。
• 减少不必要的0-1变量：储热罐充放热如果通过H_st的正负自然区分，其实可以不引入二进制变量，只要目标函数对充放热的“成本”或“收益”处理得当，优化器自然会选择合理的充放策略，不需要强制约束“不能同时充放”。但在某些场景下，充放热价格不同，就需要二进制变量。这个要根据具体模型权衡。
• 设置求解器参数：CPLEX的MIP gap容忍度从默认的1e-4放宽到0.01，求解速度可以有数量级的提升。对调度问题来说，1%的精度损失完全可接受。

ops = sdpsettings('solver', 'cplex', 'cplex.mip.tolerances.mipgap', 0.01);

• 使用热启动（warm start）：日内滚动优化时，上一轮的优化结果可以作为下一轮优化的初始解，大幅缩短求解时间。YALMIP里通过assign和optimize(..., 'solver', 'cplex', 'cplex.mip.limits...')的方式不太好做，但CPLEX本身提供了通过x0参数指定初始解的功能，利用滚动优化相邻时段解相似的特点，能显著加速。

5.3 需求响应模型调参的微妙之处

需求响应环节的参数对结果影响非常大，调节的时候要格外小心。

价格弹性矩阵的对称性：实际论文中通常假设交叉弹性是对称的，即E_ij = E_ji。但在代码里，如果不加对称性约束，优化结果会出现很奇怪的负荷转移路径——比如用户把负荷从电价高的时段转移到电价更低的时段，这本身是合理的，但如果弹性矩阵不对称，可能出现循环转移，导致负荷总量不守恒。所以建模时一定要保证：

sum(L_0_t) - sum(ΔL_t) = 总负荷约束

DR削减量的时段分布：如果DR补偿价格设置过高（比如高于机组的边际发电成本），优化器就会倾向于大量削减负荷，导致结果失真。合理的DR补偿价格应该参考系统的边际成本。我调试时一般先不加DR，跑一遍系统边际电价，然后以边际电价的1.2~1.5倍作为DR补偿价格的初始值。

5.4 数值尺度问题：兆瓦和千瓦混用会让求解器崩溃

在综合能源系统里，电功率用MW、热功率用GJ/h，它们之间的量级可能差出几个数量级。如果某些变量的系数在1e0量级，另一些在1e6量级，YALMIP会给出数值警告，CPLEX的求解效率也会大幅下降。

一个非常实用的技巧是：建模时全部采用统一的标幺值或者统一自然单位，避免MW和kW混用。我一般统一用MW和MWh，热功率统一用MWth（兆瓦热），这样电和热在数值上处于同一量级，求解器的数值稳定性好很多。储热罐容量单位也换算成MWh，而不是保留GJ，以免出现1e3量级的跳跃。

如果仍然遇到数值问题，可以在sdpsettings里开启'cplex.emphasis.numerical', 1，让求解器重点处理数值稳定性。

6. 结果分析与可视化：让调度策略“看得见”

6.1 结果分析的三张核心图

代码跑通后，分析结果时我一般会先画三张图：

第一张：电力平衡图。横轴是时间（24h或96时段），纵轴是功率（MW）。堆叠显示CHP电出力、风电出力、购电、DR削减、电锅炉耗电，与电负荷曲线对比。这张图能直观看出调度策略是怎么平衡电力的——哪些时段在买电，哪些时段在消纳风电。

第二张：热力平衡图。显示CHP热出力、电锅炉产热、储热罐充放热与热负荷的匹配情况。重点关注储热罐的充放热动作是否与电负荷的峰谷时段有逻辑关系——典型策略是：夜间风电大发时段，电负荷低谷，用多余电力通过电锅炉产热并储存在储热罐中；日间电负荷高峰时段，储热罐放热补充热负荷，让CHP减少供热、增加发电。

第三张：需求响应削减量分布图。直接从图中能看出来DR资源是否被合理利用，是否过度集中在某几个时段。

6.2 用case对比体现两阶段的价值

两阶段模型相比单阶段模型的好处，不能光靠嘴说，要设计定量对比。我的做法是设计三个场景：

• 场景A：不考虑需求响应，单阶段日前调度。
• 场景B：考虑需求响应，单阶段日前调度。
• 场景C：考虑需求响应，两阶段日前日内调度。

对比三个场景的总运行成本和风电消纳率。如果场景C的成本低于场景B、弃风率也低于场景B，就说明两阶段框架和需求响应机制确实发挥了作用。这种定量对比是论文审稿人最看重的内容。

6.3 Matlab绘图的几个实用设置

画功率平衡图时，用area函数比用plot堆叠效果好很多：

figure; h = area(t, [P_chp; P_wind; P_buy]', 'LineWidth', 0.5); hold on; plot(t, P_load, 'k-', 'LineWidth', 2); % 负荷曲线单独画在上面 legend('CHP', '风电', '购电', '电负荷'); xlabel('时间 (h)'); ylabel('功率 (MW)');

注意area需要按列输入堆叠数据，且矩阵顺序决定了堆叠的上下次序。如果显示结果发现堆叠图顺序不对，把数据矩阵转置一下就好。

还有一个非常容易被忽视的问题——坐标轴字体大小。Matlab默认的坐标轴字体在论文插图中偏小，一般需要手动设置：

set(gca, 'FontSize', 11, 'FontName', 'Times New Roman');

另外，打印图片时建议用exportgraphics(gca, 'fig1.png', 'Resolution', 300)，导出300dpi的图片，保证期刊印刷清晰。

7. 代码扩展方向：怎么把这个框架改造成你自己的项目

7.1 加入网络约束

我上面这套模型是单节点模型，即不考虑电网潮流和热网管道损耗，只考虑功率平衡。如果你的项目需要加入网络约束，可以在Matlab里引入电力系统的直流潮流模型（DC Power Flow）和热力系统的节点法模型（Heat Network Model）。注意这两种网络约束的引入都会显著增加变量数量和约束数量，求解时间也会变长。建议先用单节点模型跑通整体逻辑，再逐步加网络约束。

电力的直流潮流模型约束形式为：

P_branch = B * θ -P_branch_max ≤ P_branch ≤ P_branch_max

热网部分稍微复杂一些，需要引入节点温度变量、管道流量变量和热力平衡方程，通常表示为非线性约束。在MATLAB里处理热网非线性约束时，建议先做线性化近似，否则求解速度难以接受。

7.2 考虑不确定性：场景法或鲁棒优化

日前阶段的预测误差影响，可以通过引入不确定性集合来处理。最常见的做法是场景法（Scenario-Based Stochastic Programming），即生成多个风电出力场景，让日前优化在期望意义下最优。Matlab里生成风电场景可以用scenario generation的方法，比如基于历史数据的蒙特卡洛抽样或拉丁超立方抽样。

另一种思路是鲁棒优化，构造风电出力的盒式不确定性集合，让方案在所有可能的风电出力下都可调。鲁棒优化的代码实现主要是在目标函数和约束里引入对偶变量和不确定性预算约束，逻辑上会比场景法复杂一些。

7.3 把代码封装成可复用的函数库

如果这个方向准备长期做，建议把模型构建、求解、结果分析都封装成独立函数，通过参数配置来切换不同的算例场景。这样后续做多场景对比、参数敏感性分析时，只需要写一个循环调用这些函数，而不是每次重新改主程序。

一个典型的参数敏感性分析代码长这样：

for cap_st = 100:100:1000 % 遍历储热罐容量 params.SOC_max = cap_st; run_optimization_scenario(); % 调用主优化函数 total_cost(cap_st) = result.cost; wind_curtailment(cap_st) = result.curtailment; end

然后用一句plot就能画出储热罐容量对总成本和弃风率的影响曲线，这就是很好的论文素材。

写作并调试这套电-热综合能源系统两阶段调度模型的过程中，我最大的体会是：“两阶段”不是简单的两次优化，而是决策变量在时间维度上的分层。代码实现上的难点不在于求解器能解多大规模的问题，而在于你是否把“哪些变量应该昨天定”、“哪些变量应该今天调”这个逻辑理清楚。理清了这一点，YALMIP代码的架构就是水到渠成的事。另外，不要迷信复杂的模型——把CHP可行域、储热罐动态、DR约束这几个核心环节做扎实，用一个适中的算例跑通，得到的调度结果已经足够说明问题。