1. 双层优化调度问题的拆解与建模思路很多刚接触综合能源系统方向的同学看到核心期刊论文里那个双层优化模型第一反应通常是这代码该怎么写Yalmip里怎么表达下层问题的解作为上层问题的约束我一开始复现这篇计及需求响应的区域综合能源系统双层优化调度策略时也一样后来发现只要把模型拆明白、把求解思路理顺Matlab代码实现并不玄乎关键在于搞清楚双层优化到底在做什么。1.1 双层优化到底在优化什么先说概念。所谓双层优化本质上是两个决策主体之间存在的上下级博弈关系。在区域综合能源系统里上层通常是系统运营商手里握着各类设备的调度权目标是让整个园区的运行成本最低——买气、买电、设备启停、储能充放都是它说了算。下层呢是用户侧用户会根据运营商发布的电价、热价调整自己的用电用热行为也就是需求响应目标是让自己的用能成本最低。这里有个关键点上下层不是各干各的。下层的行为会改变负荷曲线进而影响上层的设备出力安排而上层的价格策略又会直接影响下层用户的响应幅度。这种你决策影响我、我决策又反过来影响你的结构用普通单层优化根本描述不了必须老老实实建双层模型。那为什么核心期刊爱用双层因为单层优化隐含的假设是所有人听调度中心指挥——这在严格的电网环境下成立但综合能源系统里用户是有选择权的电价定高了用户就少用电纯靠单一目标去优化结果在现实中往往跑不通。双层优化把价格信号→用户响应→系统再调整这条反馈链条显式建模了更贴近实际。数学上双层优化可以写成下面这种形式min F(x, y)s.t. G(x, y) ≤ 0y ∈ argmin { f(x, y) : g(x, y) ≤ 0 }外层是上层问题内层是下层问题。y是下层决策变量但它的取值必须来自下层问题的最优解集。这个y必须是最优解的约束就是双层优化让无数人头疼的地方——你不能随便给y赋个可行值就完事。1.2 区域综合能源系统的设备建模与能量流搞清楚双层结构后接下来就是把区域综合能源系统这个物理实体映射成约束方程。常见的区域综合能源系统包含电、热、气三种能源形式核心设备包括热电联产机组CHP、燃气锅炉、电锅炉、电储能、热储能可能还有光伏和风电。我做复现时用的是最经典的能源集线器思路把所有设备当成一个转换枢纽输入是电网购电和天然气输出是电负荷和热负荷。每个设备的建模都有固定的套路这里直接给出我用的数学表达后面代码段会一一对应CHP机组是最重要的设备它的核心特征是吃气、出电、出热。用简化模型表示P_CHP(t) η_e * F_CHP(t)H_CHP(t) η_h * F_CHP(t)其中F_CHP(t)是t时段消耗的天然气功率η_e是发电效率η_h是制热效率。这里要注意CHP的电热出力之间存在耦合关系不是两个独立变量所以约束里通常加一项电热比范围限制我还见过一种做法是直接固定热电比H_CHP(t) r_CHP * P_CHP(t)r_CHP就是热电比这个参数在原文里一般会给复现时直接用就行。燃气锅炉就简单多了纯粹是气→热的单向转换H_GB(t) η_GB * F_GB(t)电锅炉则是电→热H_EB(t) η_EB * P_EB(t)储能设备是另一个必须处理的模块。电储能和热储能都遵循同一个动态约束——能量状态随时间递推SOC(t1) SOC(t) (P_ch(t) * η_ch - P_dis(t) / η_dis) * ΔtP_ch和P_dis分别是充能和放能功率η_ch和η_dis是充放效率。储能建模有三个坑一个是SOC上下限别漏另一个是同一时刻不能既充又放——这个需要引入0-1变量来约束还有一个是始末SOC要衔接很多复现结果对不上原文就是末时段SOC没处理好。能量平衡约束是整个模型的骨架。电功率平衡和热功率平衡都要写成等式约束P_buy(t) P_CHP(t) P_dis(t) P_PV(t) P_load(t) P_EB(t) P_ch(t) P_sell(t)H_CHP(t) H_GB(t) H_EB(t) H_dis(t) H_load(t) H_ch(t)P_buy和P_sell还有一个互斥约束——同一时刻要么从电网买电要么向电网卖电这个和储能充放互斥的处理方式一样用0-1变量加上大M法搞定。搞明白了设备和约束模型的上层框架就搭起来了剩下的核心就是下层用户的需求响应建模了。2. 需求响应机制与负荷侧的数学建模需求响应是这篇论文的关键词之一也是双层优化里绕不开的一环。在复现时我发现需求响应建模如果做得粗糙整个模型可能直接变成单层——那就失去了复现的意义。需求响应的核心思想很简单用户不是死负荷电价涨了我就少用点或者把部分负荷挪到便宜时段去用。问题是怎么把这个行为量化成可优化的数学表达式。2.1 价格型需求响应弹性矩阵怎么算价格型需求响应学术点说叫Price-Based Demand ResponsePBDR用的最多的工具是需求价格弹性矩阵。基本逻辑是电价的相对变化会引起电负荷的相对变化两者之比就是弹性系数。e_tt是自弹性表示t时段电价变化对t时段负荷的影响一般为负值——电价涨用电降e_ts(t≠s)是互弹性表示s时段电价变化对t时段负荷的影响一般为正值——比如谷段电价低了用户把峰段的负荷挪到谷段来。响应后的负荷计算公式是P_load(t) P_load0(t) × [1 Σ e_ts × (ρ(s) - ρ0(s)) / ρ0(s)]其中P_load0是原始负荷ρ0是原始电价ρ是优化后的电价。这个公式实现起来不复杂但有两个细节容易踩坑。第一是弹性系数的取值。这个一般是论文里给定的或者参考行业经验值自弹性我见过用-0.2的也见过-0.3的互弹性一般在0.01到0.05之间。弹性系数取得太大响应后的负荷可能变成负值取得太小需求响应的效果又看不出来。我复现时先按原文参数试如果原文没给就用自弹性-0.2、互弹性0.03做基准再根据结果微调。第二是价格变量的问题。注意下层优化中价格ρ是上层的决策变量但下层的负荷P_load又依赖于ρ。这意味着上下层之间形成一个隐式的耦合关系。处理思路是这样在下层模型里用户的优化变量是自己各时段的负荷调整量ΔP(t)目标是用能成本最小化约束是调整后的负荷满足弹性关系。这个处理在Matlab里实现时要把弹性关系当作约束方程代入而不是作为公式直接算好代入。2.2 可削减负荷与激励型响应的接入方式价格型响应之外还有一类常见的需求响应叫激励型需求响应Incentive-Based Demand ResponseIBDR。这类响应的特征是用户与运营商签订合同允许运营商在特定时段削减一部分负荷运营商给用户补偿。在我的复现模型里这部分以可削减负荷的形式接入。可削减负荷的数学模型比价格型简单但约束条件更琐碎。核心约束有三个第一削减量有上限不能把用户负荷削到零0 ≤ ΔP_DR(t) ≤ P_DR_max(t)第二用户侧削减后不能倒送电即削减量不能超过原始负荷ΔP_DR(t) ≤ P_load0(t) - P_min第三削减行为要付出补偿成本计入上层目标函数C_DR Σ comp_price × ΔP_DR(t)补偿价格一般取一个高于电价的值比如1.2倍目录电价太低用户不签合同太高运营商成本太高——这个参数的敏感性分析我后面调试时会重点说。这里有个容易犯的错可削减负荷虽然声明在下层但它本质上更像上层的一个调度资源。很多复现代码把可削减负荷当成上层直接控制的变量这其实也没问题——只要目标函数里体现了补偿成本就行。关键是不能把价格型响应的负荷变化量和激励型削减量重复计算。我见过一种做法是电负荷先经过价格型响应修正再叠加激励型削减两者相乘等于把用户响应算了两次成本翻倍结果自然和原文对不上。需要注意需求响应的建模方式决定了双层的耦合紧密度。如果只做价格型响应那么上下层的耦合关系在电价变量上如果加激励型响应则在削减量变量上。二者的约束分组、变量归属在写代码之前就要规划清楚不然约束矩阵堆在一起求解器分分钟报错。3. Matlab代码框架从数学模型到可运行实现模型搭完终于到了写代码这一步。我对这些核心期刊复现代码的最大感受就是数学模型是一回事能跑出结果是另一回事中间隔着一个求解器脾气的距离。这一节把完整的技术路线和关键代码段梳理一遍。3.1 环境配置与工具箱选型先说我用的环境Matlab R2023b Yalmip工具箱 CPLEX求解器。这三者是目前复现能源系统优化调度最主流的搭配。Yalmip是建模层面的语法糖专门把优化问题翻译成求解器能懂的格式CPLEX负责真正求解混合整数线性规划。版本这个问题确实折腾人。Yalmip版本太老可能不支持最新的Matlab语法CPLEX版本太老可能和64位Windows系统不兼容而且CPLEX从12.10开始对Matlab版本的适配也有变化。我的建议是Matlab用2020以后的版本Yalmip直接用GitHub上的最新版CPLEX用12.9或者22.1版本。装好之后在Matlab里运行yalmiptest看到输出一堆solver installed就是成功了。有些同学电脑上装的是Gurobi那也没问题Yalmip对Gurobi的支持同样很好。只要注意一点同一个模型不要频繁切换求解器因为不同求解器对MIP的默认参数有差异跑出来的结果会有细微差别复现时容易造成我怎么和原文对不上的困惑。选定一个求解器后面所有的调试都在这个求解器上进行。3.2 双层模型转单层KKT条件与互补松弛线性化双层优化不能直接丢给求解器处理这是所有复现者面对的最大拦路虎。Yalmip再厉害也没法原生处理下层问题最优解这种约束。常见的处理办法有三种KKT条件转换法、强对偶转换法、以及启发式迭代法。这篇论文我用的是KKT条件转换法——因为底层用户的优化问题是一个线性规划它的KKT条件既是必要的也是充分的可以无损替换。KKT条件替换的思路是这样的先把下层问题的拉格朗日函数写出来然后列出三条条件——平稳性条件拉格朗日函数对各变量求导为零、原始可行性条件、以及对偶可行性条件。这三组约束再配合互补松弛条件就完整等价于下层优化问题。具体来说下层的负荷调整问题可以写作min f(ΔP)s.t. h(ΔP) ≤ 0引入对偶变量后KKT条件包含平稳性∂f/∂ΔP Σ λ × ∂h/∂ΔP 0原始可行h(ΔP) ≤ 0对偶可行λ ≥ 0互补松弛λ × h(ΔP) 0麻烦的是互补松弛条件——它是一个双线性乘积项λ和h(ΔP)至少有一个必须为0。这是个非线性约束不能直接丢给MILP求解器。标准做法是用大M法线性化λ ≤ M × zh(ΔP) ≤ M × (1 - z)z ∈ {0, 1}大M法的M值选择有讲究。M取得太大求解器数值容易出问题MIP求解时间明显变长M取得太小可能会错误地把可行解剪掉。我调试的经验是M取目标函数中各变量典型量级的10到100倍比如补偿价格是0.8元/kWh、削减量是500kW那么M取1000左右就够用了。可以先跑一版看看结果如果出现infeasible再调整M。3.3 目标函数与约束条件的Yalmip实现理论讲完直接上一个精简版的Yalmip代码框架这是我从完整的复现代码里抽出来的骨干结构读者可以直接参照这个骨架去填充自己的设备参数。%% 决策变量定义 % 上层变量 P_CHP sdpvar(24, 1, full); % CHP电出力 H_CHP sdpvar(24, 1, full); % CHP热出力 F_CHP sdpvar(24, 1, full); % CHP耗气量 P_GB sdpvar(24, 1, full); % 燃气锅炉耗气量 H_GB sdpvar(24, 1, full); % 燃气锅炉热出力 P_EB sdpvar(24, 1, full); % 电锅炉耗电量 H_EB sdpvar(24, 1, full); % 电锅炉热出力 P_buy sdpvar(24, 1, full); % 电网购电 P_sell sdpvar(24, 1, full); % 电网售电 u_ch binvar(24, 1); % 储能充电状态 u_dis binvar(24, 1); % 储能放电状态 % 下层用户的负荷调整变量经过KKT转换后并入本层 dP sdpvar(24, 1, full); % 需求响应削减量 lam sdpvar(24, 1, full); % 下层对偶变量 z binvar(24, 1); % 互补松弛线性化辅助变量 %% 目标函数 % 上层购气成本 购电成本 - 售电收益 DR补偿 运维成本 C_gas gas_price * sum(F_CHP F_GB) * dT; C_grid sum(price_buy .* P_buy - price_sell .* P_sell) * dT; C_dr comp_price * sum(dP) * dT; C_om sum(k_chp * P_CHP k_gb * P_GB k_eb * P_EB) * dT; Objective C_gas C_grid C_dr C_om; %% 约束条件 Constraints []; % 设备耦合约束 Constraints [Constraints, P_CHP eta_e * F_CHP]; Constraints [Constraints, H_CHP eta_h * F_CHP]; Constraints [Constraints, H_GB eta_gb * F_GB]; Constraints [Constraints, H_EB eta_eb * P_EB]; % CHP/GB/EB出力上下限 Constraints [Constraints, 0 P_CHP P_CHP_max, 0 H_GB H_GB_max, 0 P_EB P_EB_max]; % 联络线功率约束 Constraints [Constraints, 0 P_buy P_buy_max, 0 P_sell P_sell_max]; % 储能约束动态充放互斥 Constraints [Constraints, SOC(2:24) SOC(1:23) (P_ch * eta_ch - P_dis / eta_dis) * dT]; Constraints [Constraints, P_ch M * u_ch, P_dis M * u_dis, u_ch u_dis 1]; % 电/热功率平衡 Constraints [Constraints, P_buy P_CHP P_dis P_pv P_load0 - dP P_EB P_ch P_sell]; Constraints [Constraints, H_CHP H_GB H_EB H_dis H_load H_ch]; %% 下层KKT条件 % 平稳性条件 Constraints [Constraints, 1 lam - comp_price 0]; % 原始可行 对偶可行 Constraints [Constraints, dP - P_load0 0, -dP 0]; Constraints [Constraints, lam 0]; % 互补松弛线性化大M法 Constraints [Constraints, lam M * z, P_load0 - dP M * (1 - z)]; %% 求解 ops sdpsettings(solver, cplex, verbose, 2, mip.tolerances.mipgap, 1e-4); optimize(Constraints, Objective, ops);这段代码有几个重点要解释。第一个是互补松弛条件的处理这里的例子比较简村实际论文中下层问题会有多个约束每个约束都要配一个互补松弛条件和一个0-1变量这部分是代码量的主要来源。第二个是SOC的递推写法用向量的方式一步完成24个时段的递推比写for循环干净得多也不容易出错。第三个是储能充放互斥用u_ch u_dis 1来杜绝同时充放。很多人第一次写完这个代码跑出来结果发现目标函数是负的或者某个设备出力一直在边界上跳来跳去多半是某个约束忘写了。我调试的顺序是先检查所有0-1变量相关的约束再检查KKT条件里互补松弛的成对约束最后再查功率平衡——这三类问题占了复现调试的八成。4. 复现用的数据准备与参数调试代码写完下一步就是喂数据。很多人以为数据不重要随便编一组就往上跑结果怎么调都不对其实大部分时候问题出在参数没有自洽。这一节聊复现过程中数据到底要怎么准备哪些参数是关键中的关键。4.1 典型日数据与设备参数表区域综合能源系统最常见的分析场景是典型日调度一般取冬季典型日和夏季典型日各一组。我复现用的是一组24小时的电负荷、热负荷、光伏出力数据时间分辨率1小时。数据来源可以是原文附录找不到就用相似的典型负荷曲线代替但要注意量级匹配——比如系统总电负荷峰值如果是10MW那CHP容量大概率在5到8MW附近不能差太远。我用的关键设备参数整理成一张表大家可以对照着检查自己的参数自洽性参数名称数值单位说明CHP发电效率 ηe0.35-气转电效率CHP制热效率 ηh0.45-气转热效率燃气锅炉效率 ηgb0.9-气转热效率电锅炉效率 ηeb0.95-电转热效率电储能容量1000kWh额定容量电储能充放效率0.9 / 0.9-充电/放电效率储能SOC范围0.1 - 0.9-防止过充过放电网购电上限2000kW联络线容量天然气价格2.5元/m³折算成功率计价约0.25元/kWh分时电价峰/平/谷1.1 / 0.6 / 0.3元/kWh上层制定的初始价格自弹性系数-0.2-价格型响应关键参数互弹性系数0.03-价格型响应关键参数DR补偿价格0.8元/kWh激励型响应补偿这个表里的天然气价格要特别说一下天然气一般是按m³计量的但模型里耗气量F_CHP和F_GB用的是功率单位kW所以要把天然气价格换算成元/kWh——1 m³天然气完全燃烧大约释放9.5到10 kWh热量按2.5元/m³算相当于0.25元/kWh左右。这个换算关系如果不做目标函数里气成本会比电成本低一个数量级整个调度结果会严重失真。4.2 求解器参数与收敛性调试数据喂进去之后求解器参数是影响复现结果能否对上原文的最后一环。CPLEX求解MIP问题时我用的核心参数有三个第一个是MIP gap容忍度。Yalmip里对应的设置是mip.tolerances.mipgap我一般设1e-4表示允许最优解和当前最优整数解的相对差距在0.01%以内。设太大比如1e-2求解速度快但结果精度不够成本对不上原文设太小比如1e-6求解时间可能从几分钟膨胀到几小时得不偿失。第二个是大M值。前文说过互补松弛线性化用的大M值过大会引起数值问题。这里再补充一个检查方法求解完成后把z变量的值和λ、h的值拉出来看看如果发现λ和h同时非零逼近M量级说明M设置不合理需要缩小。第三个是初始可行解。双层问题转单层后变量多、约束多CPLEX有时要在分支定界树里搜索很久才找到第一个可行解。这时候可以先用松弛模型把0-1变量全放开为连续变量跑一遍把结果作为MIP的初始解传入能显著缩短求解时间。在Yalmip里可以用assign函数给变量赋初值再用optimize求解。% 先解松弛问题获取初值 ops0 sdpsettings(solver, cplex, verbose, 0); relax relaxdouble(Constraints, Objective); optimize(Constraints, Objective, ops0); val_dP value(dP); % 把松弛解赋给原始MIP变量 assign(dP, val_dP); % 正式求解 optimize(Constraints, Objective, ops);这个松弛→初始化→正式求解的流程非常实用特别是当模型规模大、约束多而且变量离散程度高的情况下效果立竿见影。5. 复现过程中的常见问题与排查技巧实录复现的核心期刊代码几乎不可能一次跑通。我在这个项目里前前后后调了两周把典型的坑都踩了一遍。这一节直接给出排查思路和解决办法希望后来的人少走弯路。5.1 常见问题速查表先把我在复现中遇到的几类高发问题整理成表格每个问题后面附排查方向问题现象可能原因排查方向求解器报infeasible约束条件冲突常见于功率平衡与设备上下限矛盾逐个检查约束先注释掉功率平衡跑一次再看哪类约束导致无解目标函数为负售电价格高于购电价格或补偿价格设置不合理检查分时电价参数是否出现倒挂需求响应的削减量恒为0补偿价格设置过低削减不经济提高补偿价格或检查互补松弛条件的M值是否太小储能SOC持续维持在边界SOC递推约束有误或始末SOC约束缺失拉出SOC曲线检查是否满足末时段衔接条件求解时间过长0-1变量过多或大M值过大增加初始解缩小大M设置MIP gap容忍度结果对不上原文原始负荷数据量纲或弹性系数取值不同仔细核对原文参数表特别是量纲换算这个表浓缩了我踩坑的经验遇到问题先对照一下基本能快速锁定位置。下面挑三个最有代表性的坑展开讲讲。5.2 三个最值得展开的坑第一个坑互补松弛条件M值导致的伪最优解。我当时把M设成1e6想着反正比所有变量都大就行结果求解器返回了最优解成本也比原文低不少但拉出负荷曲线一看dP居然超过了负荷总量——这已经是明显的物理不可行了。原因就是M太大互补松弛条件在数值上失去了约束力。后来把M改到1000结果完全变正常。所以M值不是越大越好够用就行。第二个坑价格型需求响应的负荷循环引用。我在写响应后负荷时一开始直接用dP作为决策变量加进电功率平衡但是弹性矩阵中的ρ又是另一个上层变量两个变量之间存在非线性乘积——结果是模型变成了非线性的求解器直接报错。后来我把弹性关系拆成约束把响应后的负荷写成原始负荷减去一个决策变量ΔP才彻底绕开这个非线性项。这里想提醒大家凡是遇到两个决策变量相乘第一反应应该是能不能线性化或者重新建模不能硬塞给求解器。第三个坑也是大家最容易忽略的始末SOC的衔接。储能调度问题如果只给了SOC初始状态0.2没约束末状态优化器会把SOC在最后一个时段全部放空——因为放空意味着多出一部分能量来满足负荷成本降低了。但实际运行中储能不可能每天结束时电量为零。解决办法是在约束里加一个等式SOC(24) SOC(1)让储能每天回到初始状态。这个约束一加结果就平稳了成本和原文对上了。这些排查经验换来一个教训复现核心期刊代码第一步不是看代码而是把论文里的每个公式都写成自己的符号系统逐个核对约束条件是否在代码里有对应实现。代码与公式的一一对应关系一旦建立调试就是查漏补缺的事情不会出现从头到尾找不到问题在哪的窘境。最后补一句个人感受这个双层优化模型跑通之后我把价格型响应和激励型响应的参数各做了几组敏感性分析发现需求响应带来的成本降幅通常在5%-8%之间而且弹性系数超过一定阈值后收益会逐渐饱和。这个结论如果写进论文的算例分析部分是很好的加分项——既说明了需求响应的价值又不至于夸大它的作用。复现论文的最高境界不是原样跑通而是在跑通的基础上有自己的发现这才是复现两个字真正的价值所在。