这篇论文我前后复现了小半个月，中间踩了不少坑。核心期刊上这类“计及需求响应的区域综合能源系统双层优化调度”的论文，看起来模型都差不多，可真到自己用Matlab把代码搭起来的时候，卡点全在上下层耦合、KKT条件线性化和求解器配置这些环节上。这篇博文我从头到尾梳理一遍我当时是怎么拆解题目、怎么建模型、怎么求解、怎么调试的，代码框架和参数都会贴出来，给正在复现类似论文的读者一个可以直接参考的路径。

先说清楚这个项目解决什么问题。区域综合能源系统（RIES）里面有很多设备——热电联产机组（CHP）、燃气锅炉、电锅炉、储能、光伏等等，它们既供电又供热，互相耦合。传统调度方式把用户负荷当成一个固定值，但引入需求响应后，用户会根据电价或者激励信号主动调整用能，这时候用户就不再是被动接受者，而是跟运营商有利益博弈的参与者。“双层优化”就是把这种博弈关系写进数学模型：上层是运营商做设备出力和购能决策，下层是用户根据上层给出的价格信号做负荷调整，两边互相影响、迭代求解，最终达到既满足用户利益、又让系统总成本尽可能低的均衡点。Matlab里做这个，主流方案是YALMIP建模加CPLEX/Gurobi求解，这也是核心期刊复现中最常用的技术路线。

这篇内容适合三类人：正在做综合能源系统优化方向的研究生，准备复现或参考核心期刊仿真代码的工程师，以及想把双层优化、需求响应这两个概念真正落成可运行代码的学习者。下面按我当时从零开始的推进顺序来写。

1. 先拆解题目：这套双层调度到底要解决什么问题

1.1 “区域综合能源系统”不是把电和气放在一起就完了

第一步，我得弄清楚区域综合能源系统里到底有哪些元件在参与调度。最典型的配置是：外部电网购电、天然气网购气、CHP机组发电供热、燃气锅炉补热、电储能削峰填谷，用户侧有电负荷和热负荷，可能还有分布式光伏。

这些设备不是独立的，核心耦合点在于CHP：它烧天然气，同时发出电和热，电热出力之间有固定的可行域关系。这个关系在建模时通常用多边形不等式描述，而不是一个简单的效率公式，因为CHP在抽汽凝汽式工况下的电热出力范围是有约束的。这一点很多人忽略，结果模型跑出来调度结果很奇怪——CHP电出力很大但热出力很小，实际机组根本做不到。

调度决策变量包括：各时段CHP的电出力、热出力，燃气锅炉的热出力，电储能的充放电功率，向上级电网购电或者售电的功率，以及必要时给的备用容量。约束条件包括电功率平衡、热功率平衡、设备出力上下限、爬坡约束、储能SOC递推约束、与外网交互功率限制。这些是上层模型的主体。

1.2 “双层”究竟在划清谁的权责

很多读者第一次接触“双层优化”会问：为什么不能把所有东西放进一个目标函数里直接求？答案在于决策主体不同，利益目标也不同。

上层是区域综合能源系统的运营商或调度中心，它决定设备出力、购能策略和售能价格，目标是让系统总运行成本最小。下层是终端用户或者说负荷聚合商，他们根据上层给出的价格和激励信号，决定用多少电、哪些负荷可以转移、哪些负荷可以削减，目标是让自己的用能费用和舒适度损失之和最小。

这种结构本质是Stackelberg博弈：上层是领导者，先给出价格信号；下层是跟随者，基于这个信号做最优反应；上层的决策要考虑到下层的反应，不是单方面拍板。把这种关系写进优化模型，就是双层规划。如果简单把用户负荷当成固定输入，那就退化成单层优化，需求响应等于没做。

1.3 需求响应怎么落进模型里

需求响应有两类常见形式，建模方式也不同。

价格型需求响应：上层给出分时电价，用户根据电价高低调整用电时段，把峰时段的负荷往谷时段挪。这类响应可以用价格弹性系数矩阵来描述，或者在下层模型里显式建模可转移负荷。

激励型需求响应：运营商会跟用户签订协议，约定在特定时段用户需要削减一定量的负荷，用户因此获得补偿。这类响应在下层模型里对应一个“可削减负荷”决策变量，目标函数里有一项削减补偿收益，同时约束削减量不能超过协议上限，还要避免过度削减导致用户不满。

我复现的时候用的是激励型和部分可转移负荷的混合模型，因为纯价格弹性系数矩阵在实际求解时容易把负荷响应算得过于理想，算出来的曲线失真。显式建模可转移负荷和可削减负荷，虽然变量更多，但物理意义清楚，调试起来也方便。

2. 模型搭建：上层、下层和耦合变量的Matlab建模

2.1 上层调度模型：什么算成本、什么算收益

上层目标函数我写成四个部分加起来，再减掉售电收益：

• 购电费用：从上级电网购电，按外部分时电价结算
• 购气费用：CHP和燃气锅炉消耗天然气，按气价结算
• 设备运维费用：各时段各设备的出力乘以单位运维成本
• 需求响应补偿费用：对用户削减负荷和转移负荷的补偿
• 售电收益：如果允许用户直接向电网售电，或者在内部售电，可以有收益项

举个例子，购电费用这一项在Matlab里用向量点乘就能算：

C_buy_elec = sum(price_buy .* P_buy * dt);

其中price_buy是外部购电价格向量，P_buy是各时段购电功率决策变量，dt是时段长度。注意单位统一，我全程用kW和元，时间用小时，这样算出来的成本单位就是元，不容易乱。

设备运维成本类似：

C_om = sum(k_chp * P_chp + k_gb * H_gb + k_es * abs(P_es));

储能那块我用abs(P_es)算充放电损耗成本，这里P_es为正表示放电，为负表示充电。如果求解器对绝对值项支持不好，可以用两个非负变量分别表示充和放。

2.2 下层用户模型：效用函数与满意度

下层用户的优化目标，是让“用电费用减去补偿收益、再加上舒适度损失”最小。舒适度损失就是削减负荷给用户带来的不便，通常用削减量的二次函数表示，系数越大说明用户越不愿意被削减。

关键约束有这么几个：

• 各时段削减量不超过协议上限：0 <= P_cut(t) <= P_cut_max(t)
• 可转移负荷的转移量在各时段的总和保持平衡，比如一天内转移出去的负荷总量等于转移进来的负荷总量
• 削减后的实际负荷不能低于物理最低负荷，避免把负荷削到零这种离谱结果

下层目标函数是二次的，但约束全是线性，所以本质是一个凸二次规划（QP）。这个性质很重要，因为只有当下层问题是凸问题时，才能用KKT条件转成上层约束，这也是后面单层化的前提。

2.3 耦合关系与数据接口

上下层通过什么变量耦合？我需要把这条链路理清楚。

上层传给下层的，是各时段的电价信号和需求响应补偿价格。下层收到这些信号后，求解自己的优化问题，得到各时段的负荷削减量、可转移负荷的调整方案，然后把调整后的负荷曲线返回给上层。上层看到新的负荷曲线后，重新调整设备出力和购能计划。

在代码层面，这个耦合关系体现在两个地方：一是在上层目标函数里，需求响应补偿费用的计算依赖下层返回的削减量变量；二是上层功率平衡约束里的电负荷、热负荷，不再是固定参数，而是“基础负荷减去削减量、再加上转移负荷”的表达式。

如果不用KKT单层化，而是用迭代求解，数据接口就是两层之间的通信函数，循环里每次更新价格和负荷。如果用KKT单层化，这些耦合变量就会变成同一层模型里的变量和约束，求解器一次性算出均衡点。

3. 求解策略：KKT单层化还是迭代嵌套

3.1 三类常用求解思路对比

我复现过程中查了不少相关资料，双层优化的求解方法大体可以分成三类：

第一类是嵌套迭代法。上层用启发式算法（比如粒子群、遗传算法）做，下层对每个个体调用求解器精确求解。优点是思路简单、代码容易理解，缺点是计算量很大，上层每迭代一次，下层就要反复求解几十上百次，而且启发式算法不保证收敛到全局最优。

第二类是KKT单层化。把下层问题的最优性条件（KKT条件）作为约束加到上层模型里，这样双层问题就变成一个单层的带互补约束的数学规划问题，如果是线性或二次规划，还能线性化成混合整数规划，直接用CPLEX/Gurobi求解。优点是可以精确求解，缺点是KKT推导和线性化过程容易出错，大M参数没选好会引发数值问题。

第三类是元模型或者解析反应函数法。先推导下层对上层决策变量的解析反应函数，再代入上层。这个方法理论上最漂亮，但对模型结构要求高，实际问题很难写出解析式。

核心期刊论文里最常出现的是第二类，尤其在计及需求响应的背景下，下层用户的QP问题完全可以KKT化。所以我复现时也走这条路。

3.2 KKT条件转化与线性化细节

把下层问题KKT化，要写四组条件：

• 拉格朗日函数对各变量求导等于零，也就是驻点条件
• 原问题可行性条件，就是下层所有约束都要保留
• 对偶可行性条件，就是所有不等式约束对应的拉格朗日乘子非负
• 互补松弛条件，比如削减量上限约束P_cut <= P_cut_max对应的互补条件

互补松弛是非线性的，这是单层化过程中最麻烦的地方。它的形式是lambda * (P_cut_max - P_cut) = 0，也就是说，乘子和松弛量不能同时为正。非线性项不能直接交给MILP求解器，要用大M法和二进制变量做线性化。

具体地，对不等式约束g(x) <= 0引入松弛变量s = -g(x)，则s >= 0。互补条件lambda * s = 0等价于下面这组线性约束：

% 设 z 是二进制变量，M是大M常数 Constraints = [Constraints, s >= 0]; Constraints = [Constraints, lambda >= 0]; Constraints = [Constraints, s <= M * z]; Constraints = [Constraints, lambda <= M * (1 - z)];

当z=0时，lambda必须为0，s可以取0到M之间；当z=1时，s必须为0，lambda可以取0到M之间。这样就把“两者乘积为0”的非线性条件拆成了混合整数线性条件。

这里最要命的是M的取值。M太小，会砍掉真正的可行解；M太大，会让约束矩阵的病态程度加剧，CPLEX求解时会出现数值警告。我的做法是根据物理边界推算，而不是随便取大数。比如负荷削减量上限如果是基础负荷的20%，那松弛量s的上界就是该时段最大负荷的20%，大M在这个量级的基础上再放大一倍就够了。

3.3 关键代码片段与Matlab实现结构

我最终代码的主干流程是这样的：

1. 加载系统参数和负荷数据
2. 用YALMIP定义上层决策变量（设备出力、购能、补偿价格）
3. 定义下层KKT条件涉及的变量（削减量、转移量、对偶乘子）
4. 把所有约束和线性化后的互补约束拼起来
5. 定义上层目标函数（其中包含需求响应补偿费用）
6. 调用CPLEX求解混合整数规划
7. 后处理：画负荷曲线、设备出力图、分析成本构成

关键代码段大概是这个样子：

%% 上层决策变量 P_chp = sdpvar(1, T, 'full'); H_chp = sdpvar(1, T, 'full'); P_gb = sdpvar(1, T, 'full'); P_buy = sdpvar(1, T, 'full'); P_es = sdpvar(1, T, 'full'); % 正为放电，负为充电 %% 下层用户的变量 P_cut = sdpvar(1, T, 'full'); % 各时段削减负荷 P_trans = sdpvar(1, T, 'full'); % 各时段可转移负荷（正值表示转入） lambda1 = sdpvar(1, T, 'full'); % 削减量上限约束的对偶乘子 lambda2 = sdpvar(1, T, 'full'); % 削减量非负约束的对偶乘子 z1 = binvar(1, T, 'full'); % 互补约束线性化引入的二进制变量 z2 = binvar(1, T, 'full'); %% KKT驻点条件示例（简化形式） % 用户目标对 P_cut(t) 求导等于零 for t = 1:T Constraints = [Constraints, ... -price_sell(t) + alpha * P_cut(t) + lambda1(t) - lambda2(t) == 0]; end %% 互补松弛线性化 for t = 1:T Constraints = [Constraints, ... % 对应 P_cut(t) <= P_cut_max(t) P_cut_max(t) - P_cut(t) >= 0, ... P_cut_max(t) - P_cut(t) <= M1 * z1(t), ... lambda1(t) <= M1 * (1 - z1(t)), ... lambda1(t) >= 0]; Constraints = [Constraints, ... % 对应 P_cut(t) >= 0 P_cut(t) >= 0, ... P_cut(t) <= M2 * z2(t), ... lambda2(t) <= M2 * (1 - z2(t)), ... lambda2(t) >= 0]; end

这里alpha是用户舒适度损失函数的二次项系数。price_sell是上层给用户的售电价格，在双层模型里它可以是固定参数，也可以是上层决策变量。如果是决策变量，上层目标函数里会多出一个二次项，模型从MILP变成MIQP，CPLEX也能处理，但速度会慢一些，建议先跑通固定价格的版本再加码。

YALMIP定义变量的时候要注意'full'参数，默认的方形矩阵定义不适合一维决策变量，不写这个参数很容易在后面拼接矩阵约束时出现维度错误。

求解设置方面：

ops = sdpsettings('solver', 'cplex', 'verbose', 2, ... 'cplex.mip.tolerances.integrality', 1e-5, ... 'cplex.mip.tolerances.mipgap', 1e-3); sol = optimize(Constraints, Objective, ops);

跑完之后务必检查sol.problem是否为0，我用一行断言或者打印来确认：

if sol.problem ~= 0 error('求解失败: %s', yalmiperror(sol.problem)); end

这一步看似多余，实际很有用。我遇到过好几次求解器返回了结果但其实是不可行条件下的近似解，如果不检查problem字段，后面画图和分析全在拿错误结果做文章。

4. 参数设置与场景设计

4.1 电价、气价和需求响应补偿参数

这类项目的仿真结果很大程度取决于价格数据，我建议不要随便编一组数就上。我采用的典型场景参数是：分时电价峰平谷三个时段，峰时段电价大约是谷时段的2到3倍，这样才看得出需求响应削峰的效果。购气价格按热值折算后，要让CHP在热电联产工况下有经济性优势，否则调度结果会极端——全部用外购电和燃气锅炉，CHP变成摆设。

需求响应补偿价格的设计更讲究。补偿价格定得低，用户不愿意削减负荷，需求响应等于没起作用；定得高，运营商成本反而增加，甚至比直接购电还贵。我的经验是让补偿价格处于用户舒适度损失成本的1.2到1.8倍之间。这样下层用户会积极响应，但又不至于响应过度。

4.2 设备参数与负荷数据的组织方式

设备参数我推荐用一个结构体集中管理，而不是在脚本里到处写数字。一个典型的数据结构是这样：

param.CHP.P_min = 30; % kW param.CHP.P_max = 200; % kW param.CHP.H_min = 20; % kW param.CHP.H_max = 180; % kW param.CHP.eta_elec = 0.35; param.CHP.eta_heat = 0.45; param.GB.eta = 0.9; param.ES.P_max = 50; % 储能最大充放电功率 param.ES.eta = 0.95; param.ES.E_max = 200; % 储能容量

负荷数据我是用Matlab里的表对象读入，可以是Excel文件或者CSV，24个时段的电负荷、热负荷、可转移负荷比例和可削减负荷比例。基础负荷曲线设计时可以故意在晚高峰设一个突出尖峰，这样才能在结果对比图里明显看到需求响应把它削下来的效果。

4.3 对照实验与灵敏度分析的常用做法

复现论文不能只跑一个场景，核心期刊的论文一般要有对照实验。我这里至少跑了四个场景：

一是无需求响应单层优化，相当于把用户负荷当固定值，作为基准线。二是有激励型需求响应的双层优化，就是本文模型。三是不同补偿价格下的需求响应结果，观察负荷削减量和总成本的趋势。四是让储能容量或光伏容量变化，做灵敏度分析。

每个场景跑完之后，我会记录三个指标：系统总运行成本、峰时段最大负荷削减率、用户总费用变化。这些指标整理成表，再画负荷曲线对比图和设备出力堆叠图。核心期刊论文里的图基本都是这个套路——上面是电负荷平抑曲线，下面是各设备出力堆叠图，加个成本对比表。这套后处理代码建议一开始就写好，不要等模型通了再补，不然中间调参时完全靠肉眼比较曲线，效率太低。

5. 常见问题与排错实录

5.1 YALMIP/CPLEX环境问题

我在环境配置上浪费了至少一天。首先声明，Matlab版本直接影响YALMIP的兼容性。我一开始用的是Matlab R2021b，装的YALMIP版本比较老，对binvar的维度处理有问题，后来换到R2023b并升级了YALMIP才解决。CPLEX我用的是12.10学术版，需要注意CPLEX的官方安装包里路径不能有中文，否则Matlab调用动态链接库时会报Unable to load cplexmex之类的错误。

在Matlab里验证求解器是否被YALMIP正确识别，就跑一行命令：

yalmiptest

输出的列表里CPLEX和Gurobi那一行必须是found，否则后面求解时YALMIP会悄悄换成内嵌求解器，算得又慢又不对。

5.2 模型不可行与数值病态

我最常遇到的错误是求解器返回infeasible。这时候我不会直接去翻约束，而是用YALMIP的assign和check命令逐条检查约束的残差：

assign(Constraints, sol) % 不能这么直接用，正确做法是 % 手算出每组约束的松弛量，看哪条违反最严重

更实用的办法是先把模型拆成几块分别求解：只求上层不考虑下层KKT约束，或者只求下层不问上层目标，确定哪一块已经不可行。我曾经发现问题是电平衡约束里的负荷表达式写错了，基础负荷减去削减量之后，晚高峰负荷变成负数，导致功率平衡永远无法满足。这种问题单纯看求解器的报错信息根本看不出来，必须自己追踪数据流。

数值病态方面，单位不统一是最常见的。有的变量用的是kW，有的地方我一开始用MW，导致矩阵条件数巨大，CPLEX警告数值问题。后来我统一全部用kW和元，所有约束的量级都控制在三位数以内，求解顺利多了。

5.3 双层迭代不收敛，以及KKT单层化的坑

如果读者选择用粒子群嵌套求解下层，最常见的问题是双层迭代不收敛，每次都震荡。我建议不要盲目加大迭代次数，而是看上下层交互变量的收敛轨迹。如果负荷削减量在两个数值之间来回跳，大概率是上层决策给下层的信号不连续，比如价格信号在某个取值附近导致下层解剧烈变化。处理办法是给迭代过程加一个松弛因子，每次只更新一部分变量：

lambda_user_new = lambda_past + rho * (lambda_candidate - lambda_past);

这个思路参考了ADMM的更新思想，虽然不是严格意义的ADMM，但确实能改善震荡。

走KKT单层化路线也有坑。第一个坑是大M值设置不当，我在前面已经强调过。第二个坑是互补松弛条件的二进制变量数量会很大，一个不等式约束对应一个二进制变量，24时段乘上几个约束，可能产生上百个二进制变量，加上设备启停变量，模型规模会膨胀。我优化的小技巧是：能合并的约束就合并，比如负荷削减量上下界可以合写成一个区间约束的两个互补条件，用同一个二进制变量不同取值去表示，变量数量能省不少。

5.4 结果异常时的排查思路

最后补充一条排查经验。有时候求解器报告的optimal value看似合理，但画出来的负荷曲线或者设备出力曲线明显不合理，比如储能持续满充满放，或者CHP热出力突破了可行域边界。这时候要检查的是后处理代码里的取值逻辑。用value(P_chp)取决策变量数值的时候，如果变量声明时用了矩阵形式，取值顺序可能和预期不一致。我是通过打印每个时段的数值，跟输入的负荷数据逐个对照才发现的。

另一个常见问题是双层目标函数里的补偿费用和实际计算出来的削减量不匹配。原因是上层目标函数里我用的是sum(dr_price .* P_cut)，但下层KKT转化后的P_cut变量在拼接约束时被重复声明了两份，其中一份没有参与目标函数计算。这种重复声明错误YALMIP不会报错，只会在结果里悄悄出错。排查时我用了size检查每个sdpvar变量的维度，发现P_cut被声明成了2×T的矩阵，问题才暴露出来。

我个人实际操作中的体会是，复现这类双层优化论文，最考验人的不是数学推导，而是把推导结果落到Matlab代码里时对变量、约束、数值量级和求解器特性的把握。建议第一次做的时候，先把一个不含需求响应的简化单层模型完整跑通，再逐步加入下层KKT条件和需求响应环节。每加一部分就验证一次，不要一口气把完整模型堆在一起调，否则出错时根本定位不了原因。这个思路我每复现一次论文都用，几乎成了固定流程。后续如果想把模型扩展到多区域或者加入碳交易机制，代码框架其实不用大改，多加几个区域模块和对应的耦合约束就行。