做综合能源系统IES规划的人应该都有同感模型写起来容易算起来要命。设备选型是0-1整数变量容量和运行功率是连续变量电、气、热网络耦合在一起再加几十个典型日的生产模拟直接丢给求解器要么内存爆掉要么几天不出结果。我最近在做一个园区级综合能源系统规划项目最后用广义Benders分解法GBD把问题拆成“投资决策主问题”和“运行模拟子问题”两层配合Matlab实现把一个混合整数非线性规划MINLP变成若干个线性规划LP的迭代求解整个过程收敛很稳定代码量也不算吓人。这篇文章就把建模、算法、Matlab实现细节和调试经验完整记录下来给正在做IES规划、或者被非凸大模型劝退的同学一个可直接上手的参考。广义Benders分解法不是新东西1962年Benders就提出了经典版本后来广义框架又补上了非线性凸问题这一环但我在项目里见过太多论文把重点放在定理证明上一到实际落地就卡在割平面怎么加、对偶乘子怎么提、子问题不可行该怎么办这些细节。所以我这篇不是复述教材而是从工程视角把每一步讲透尤其适合已经会建模、但不想把时间耗在求解器玄学上的读者。1. 综合能源系统优化规划问题到底在优化什么1.1 规划问题的物理对象与决策变量综合能源系统规划说白了就是回答三件事设备买不买、买多大、怎么运行。常见候选设备包括燃气轮机、燃气锅炉、电锅炉、热泵、光伏、风电、电储能、热储能甚至冰蓄冷。它们通过电母线、热母线、气母线耦合向负荷侧供应电、热、冷。规划问题的决策变量可以分三层。第一层是0-1选型变量表示第i类候选设备是否被建设第二层是连续容量变量表示设备安装功率或容量第三层是运行变量包括每个典型日、每个小时的设备出力、储能充放电功率、从电网购电、气网购气等。如果做的是多阶段扩展规划还要加一个“何时建设”的时间维度模型规模迅速膨胀。所以写之前先想清楚你的模型是年度一次性规划还是多阶段规划典型日选几个每个典型日按24小时还是按96刻录这些直接决定后续Benders分解的粒度。我常用的做法是候选设备不超过10类时先按3个典型日模拟每一类典型日代表春夏季或冬季某类天气的负荷曲线这样既能抓住多能互补特性又不会让子问题LP规模失控。1.2 目标函数与约束的数学表达目标函数一般是全生命周期总成本最小包括投资年值、年运行维护成本、燃料成本、外购电成本有时还考虑排放成本和可再生优先消纳。写成公式就是min C_inv C_om C_fuel C_gridC_inv是设备投资的年等值成本通常用“容量×单位容量投资×折现系数”计算折现系数和寿命强相关。C_om可以是可变运维成本按出力线性计费。C_fuel是气耗成本燃气轮机和燃气锅炉的耗气量与出力和效率挂钩。C_grid是向电网购电的费用按分时电价计算。约束条件有几类必写。第一类能量平衡约束电母线、热母线、气母线在每个时段的供应等于负荷加损耗。第二类设备运行范围约束出力不超过容量爬坡速率有上限启停状态和出力耦合。第三类储能约束荷电状态SOC动态递推、充放电功率限制、周期始末SOC相等。第四类投资与容量耦合约束比如某类设备不选时其容量必须为0表达式是 cap ≤ M×xM是大数。难点在于设备效率不是恒值。燃气轮机满负荷和部分负荷效率不同光伏出力受辐照度影响热泵COP随环境温度变化。这些非线性特性如果直接写进MILP求解器会很难受。我的建议是先用分段线性化把它们变成一系列线性约束因为Benders分解的收敛性证明对凸性有要求尽量保持模型线性化才能保证割平面有效。分段线性化的误差控制在3%以内对规划结果影响很小但求解稳定性提升一个量级。1.3 为什么这类问题需要分解算法直接单层求解整一个大规模MILP不是不行但效率很玄学。以我的小型算例为例5类候选设备、3个典型日×24小时二进制变量加连续变量大约2万个约束3万个用Cplex直接求解分支定界树可能铺到几百万个节点很多情况下几个小时内也拿不到好gap。为什么因为整数选型变量和运行连续变量高度耦合。设备选定了运行优化其实是个相对简单的LP设备没选定分支定界就要反复猜测整数组合。Benders分解的核心思想正是把整数变量“投影”出去让主问题只负责猜测选型然后用运行子问题的对偶信息生成切割平面逐步排除掉坏的整数组合。这样真正大规模求解的是纯连续LP子问题Cplex/Gurobi求解一个2万变量LP通常只要几百毫秒到几秒比硬啃MILP划算得多。当然分解不是免费的。它需要反复迭代每次都要重新求解主问题MILP。如果整数变量特别多比如超过100个主问题本身可能成为瓶颈。所以GBD适合“外层整数决策相对少、内层连续问题规模大”的结构这也正是IES规划最常见的形态。2. 广义Benders分解法把大问题拆成能解的小问题2.1 从经典Benders到广义Benders经典Benders分解是1962年提出的最初针对混合整数线性规划。它要求子问题是线性规划通过子问题对偶乘子构造线性割。广义Benders分解在1972年扩展了适用范围利用拉格朗日对偶和凸分析把Benders割推广到了一类凸非线性规划问题也就是“固定复杂变量后剩余子问题是凸优化问题”的情形。IES规划里的模型如果是线性化的那么经典Benders和广义Benders在实际操作上没有本质区别。但广义Benders的框架更灵活比如热管网水力方程、设备效率曲线做凸松弛后依然可以用同样的对偶割思想。因此很多文献直接写“广义Benders”即使最终代码里求解的还是LP。从我的经验看理解清楚“投影”比记住公式重要。原问题有两类变量复杂变量x选型0-1变量和简单变量y运行变量。如果我先拍脑袋给定x剩下关于y的子问题就是一个普通LP很容易求解。反过来看原问题关于x的最优值函数是某个凸函数的下包络Benders割就是用子问题对偶信息去逼近这个包络。每一次迭代都在给这个包络加一个线性支持面直到逼近足够紧。2.2 主问题、子问题与割平面的数学导出为了后面讲代码先把符号固定。原问题写成min c^T x d^T ys.t. A x B y ≥ bx ∈ X, y ≥ 0其中x是0-1整数变量X是选型可行域y是连续运行变量。给定一组x x^k运行子问题是SP(x^k): min d^T ys.t. B y ≥ b - A x^ky ≥ 0这个子问题的对偶最优解为λ^k。如果子问题可行且有界根据对偶理论原问题关于x的真实运行成本下界可以由下面的最优性Benders割逼近η ≥ (λ^k)^T (b - A x)如果子问题不可行说明x^k这个选型方案无法满足约束需要生成可行性割。此时可以求一个极射线ν^k使得B^T ν^k ≤ 0且ν^k ≥ 0满足ν^k^T (b - A x^k) 0于是可行性割为0 ≥ (ν^k)^T (b - A x)大白话版本可行割告诉主问题“这个选型组合行不通离它远一点”最优割告诉主问题“即使你选了某个方案运行成本也不会低于这条线性函数”。主问题则变成min c^T x ηs.t. 对每个最优割 kη ≥ (λ^k)^T (b - A x)对每个可行割 l0 ≥ (ν^l)^T (b - A x)x ∈ X, η ∈ R主问题的目标就是投资成本加上一个对运行成本的当前下界估计。每迭代一次往主问题里加一条线性割割越多η越接近真实运行成本主问题的解越接近全局最优。2.3 算法流程与收敛判断整个GBD流程可以写成下面几步初始化选一个可行的初始x^0或者用启发式构造置上界UB∞下界LB-∞割集为空迭代次数k1。求解主问题得到当前最优选型x^k、容量参数以及η^k。将主问题目标值c^T x^k η^k赋给LB。固定x^k求解运行子问题SP(x^k)。子问题不可行时生成可行性割并加入主问题子问题可行时得到运行成本d^T y^k将c^T x^k d^T y^k作为当前上界候选更新UB同时生成最优性割并加入主问题。判断UB - LB是否小于给定容差tol。如果是停止并输出当前最优方案否则kk1回到第2步。每一步的割平面都让主问题的松弛边界更贴近原问题。理论上如果模型是凸的算法有限次或很快接近收敛。实际工程碰到非凸模型比如效率曲线没有凸化GBD可能收敛到局部最优或震荡这时候就要靠分段线性化和初始化技巧补救。3. Matlab实现整体框架与核心代码解读3.1 建模工具选型Yalmip还是纯Matlab我强烈建议用Yalmip来搭模型后端接Gurobi或Cplex。Yalmip在Matlab里写约束非常接近数学表达而且对偶乘子提取特别方便Benders割直接可以拿到数组型乘子。如果只用Matlab自带的linprog和intlinprog也可以实现但提取对偶乘子的流程要手工管理约束顺序代码会又臭又长。不少刚入门的同学被Matlab下载安装和环境配置折腾够呛这里多说一句Matlab配Gurobi其实只需要装好Gurobi在Matlab里执行一波gurobi_setup然后Yalmip的sdpsettings里指定solver为gurobi就行。Cplex新版本需要额外的MATLAB支持包Gurobi兼容性更省心。我做这个算例用的是YalmipGurobiWindows下没有任何问题。主问题涉及MILP子问题是纯LP所以求解器至少要有连续LP和MILP能力。Gurobi和Cplex都满足。如果没有商业License用学术License或退而求其次用intlinproglinprog也能跑通只是千万规模以上的算例可能吃力。3.2 数据组织先用struct把设备参数装好写Benders代码最忌讳所有参数散落在一堆全局变量里。建议先把算例数据封装成几个struct示例% 候选设备参数 dev(1).name GT; % 燃气轮机 dev(1).inv 4200; % 单位容量投资元/kW dev(1).fix 120; % 固定运维元/kW/年 dev(1).eff 0.35; % 额定发电效率 dev(1).life 20; % 寿命 dev(1).capacity_min 0; dev(1).capacity_max 2000; dev(2).name GB; % 燃气锅炉 dev(2).inv 900; dev(2).fix 30; dev(2).eff 0.92; dev(2).life 20; dev(2).capacity_min 0; dev(2).capacity_max 5000;负荷和能源价格也按典型日存好这里不再展开。关键点在于每个典型日都要有独立的运行变量和约束主问题只有跨典型日共享的选型变量和容量变量。运行时子问题可以直接对三个典型日并行求解因为典型日之间解耦Benders割按典型日分别生成后再一起返回主问题。3.3 主问题与子问题的核心代码结构下面是核心框架我按适用性做了精简。先写主问题x binvar(n_dev,1); % 0-1选型 cap sdpvar(n_dev,1); % 连续容量 eta sdpvar(1,1); % 运行成本估计 mp_obj sum(devCost.*cap) eta; mp_con [cap_min cap cap_max]; for i 1:n_dev mp_con [mp_con, cap(i) dev(i).capacity_max * x(i)]; end mp_con [mp_con, cuts, sum(mp_cost) budget]; % 再加预算约束 mp_opt sdpsettings(solver,gurobi,verbose,0);其中cuts是之前迭代生成的所有Benders割用Yalmip约束数组保存。每轮迭代在切割旧约束的同时把新割加进去。子问题部分以第k次迭代为例。为了可读性先把x_k和cap_k转成double值固定进子问题约束x_val value(x); cap_val value(cap); % 典型日循环生成子问题 sp_cost sdpvar(n_dev, T); % 设备出力 sp_charge sdpvar(n_dev, T); % 充放电功率 sp_con []; for t 1:T sp_con [sp_con, sum(sp_cost(:,t)) grid(:,t) ... elec_load(t) heat_load(t)]; end % 其他约束设备出力上下限、储能递推、母线平衡 sp_obj sum(sum(opCost .* sp_cost)) sum(gridCost .* grid); sp_opt sdpsettings(solver,gurobi,verbose,0); solvesp optimize(sp_con, sp_obj, sp_opt);如果solvesp.problem 0表示子问题可行这时提取对偶乘子lambda dual(sp_con(1)); % 第1组约束的对偶乘子 % 构造最优性Benders割 newCut (eta value(sp_obj) - lambda*(A*[x;cap] - A*[x_val;cap_val])); cuts [cuts, newCut]; UB min(UB, invCost(x_val, cap_val) value(sp_obj));这里的dual提取需要小心。Yalmip中dual是对整个约束数组返回一个数组如果约束被拆成多条最好用索引对应。最省事的方法是把目标与约束整体写成标准型然后用Yalmip的dual函数逐一取出。如果子问题不可行则需要生成可行性割。在Matlab中可以直接调用Gurobi返回的Farkas乘子或使用recover加dual。一个稳妥的代替方案是给子问题加个无穷大惩罚松弛变量先算出“缺额”再用缺额对应的对偶信息构造割。我用这个方法避开了Farkas乘子提取的麻烦效果不错把子问题写成最小化松弛惩罚的形式若最优解中松弛量大于容忍阈值说明当前选型不可行生成的割就是该松弛子问题的“最优割”同样能引导主问题避开不可行区域。3.4 迭代主循环与收敛监控迭代循环是最容易写出bug的部分。我的建议是把收敛判据、割数量、上下界历史都记录下来方便事后查问题。核心循环如下UB inf; LB -inf; gap inf; iter 0; cuts []; while gap 1e-3 iter max_iter iter iter 1; optimize(mp_con, mp_obj, mp_opt); LB value(mp_obj); x_val value(x); cap_val value(cap); [subCost, feas, newCut] solveSubproblem(x_val, cap_val); if feas UB min(UB, sum(invCost(x_val, cap_val)) subCost); cuts [cuts, newCut]; else cuts [cuts, newFeasCut]; end gap (UB - LB) / max(1e-6, UB); record(iter) struct(iter,iter,LB,LB,UB,UB,gap,gap); end写这段的时候有两点经验。第一初始UB非常重要。可以先取一个很保守的“全选”方案把所有候选设备都选上求一次运行子问题得到一个真实可行运行成本UB就有了。如果初始UB给的是inf前几次迭代gap巨大容差判断会很难看。第二主问题求解时给一个时间限制比如sdpsettings里设置gurobi的timelimit为120秒配合MIP Gap设置。否则每次迭代主问题都要花大量时间证明最优性整体迭代次数虽然不多总时间却不划算。3.5 结果输出与可行性检查程序跑完别忘了做“后验”。至少输出四样东西选型结果、容量结果、典型日运行曲线、收敛曲线。选型结果可以用表格打印容量结果要检查是否超过最大容量、是否与选型变量逻辑一致。运行曲线要检查能量平衡残差Yalmip求解器在数值容差内可能允许微小不平衡规划方案里要强调这是可接受的。我习惯把每次主问题解出的x_mp也保存下来因为同一组割下可能出现多个整数解都满足割平面。最终上界对应的x_val才是真正通过子问题验证的方案输出时以上界对应方案为准不要用主问题最后一次解直接拿去出报告。4. 算例测试一个小型园区IES的规划结果4.1 算例基础数据为了演示我构造一个简化但完整的园区算例。候选设备5类燃气轮机、燃气锅炉、电锅炉、光伏、电储能。负荷为冬季典型日的电负荷与热负荷24小时分辨率。电价峰平谷峰值1.2元/kWh18:00-22:00、平段0.75元/kWh、谷段0.35元/kWh天然气价格2.8元/m3热值为9.7 kWh/m3。设备参数如表1。设备单位投资(元/kW)效率/COP运行维护(元/kWh)寿命(年)容量上限(kW)燃气轮机GT4200发电0.35余热回收0.450.0320800燃气锅炉GB9000.920.005201500电锅炉EB15000.980.01201000光伏PV3500按辐照曲线0.005251200电储能BESS1800元/kWh充放效率0.920.0210500kWh典型日负荷峰值电负荷780kW热负荷950kW光伏辐照最高在13点。规划目标是满足全年典型日运行约束下年总成本最小取折现率6%。4.2 求解过程与收敛曲线把数据塞进上面的Matlab框架子问题用Gurobi求解主问题也是Gurobi。初始UB用“全选”方案得到大约249.6万元/年LB从-∞开始第一次主问题解出来大概是196.8万元gap约21%。之后每次迭代添加1条最优割在迭代第8次时gap降到2.1%第14次降到0.4%第17次gap小于0.1%程序终止。总耗时约7.2秒其中主问题耗时约4.1秒子问题耗时约2.3秒剩余是I/O和Yalmip生成时间。收敛曲线呈现典型的“LB阶梯上升、UB平稳下降”的形态下界被割一步步抬高上界在刚开始就落到一个比较准的位置。这说明初始化全选方案虽然不是最优但距离真实方案不远所以UB没有剧烈波动。如果初始方案给得太差UB会从很大值缓慢下降gap收敛曲线会拖得很长。4.3 规划结果与单层求解对比优化结果燃气轮机选型为1台容量650kW燃气锅炉选型为1台容量520kW电锅炉不选光伏选择容量480kW电储能选择90kWh容量。项目的年化总投资约82.6万元年运行成本约158.4万元总年成本241万元。光伏在白天压低购电成本燃气轮机在晚高峰大量出力电储能利用谷电充电、峰时放电总体趋势符合预期。作为对照我把同样模型用单层MILP直接交给Gurobi求解设置MIP Gap0.1%。由于模型规模不大单层也能解出来耗时约58秒内存峰值是GBD的3倍以上。两者得到的总成本相差不到1.5%基本验证了GBD解的可靠性。算例规模扩大后单层法很可能提前内存爆炸GBD的分解优势会更明显。5. 常见问题与调试经验5.1 子问题不可行最容易翻车的环节固定x后子问题不可行最常见原因是容量给得太小某个时段的负荷无法满足。比如主问题给出了不选燃气轮机的方案偏偏冬季热负荷峰值很高热平衡约束无论如何都凑不齐。此时如果不生成可行性割主问题会反复出现同一个不可行方案白跑很多轮。我的处理方式是用“松弛变量法”。给子问题每个平衡约束加一个松弛变量目标改为最小化惩罚值。如果最优解里松弛总量大于1e-4就判定不可行并把当前松弛子问题的对偶乘子拿来构造可行性割。这个做法比直接提Farkas射线更容易在Matlab里稳定实现。具体来说把子问题写成min d^T y M*(s1 s2)s.t. B y s1 - s2 ≥ b - A xs1, s2 ≥ 0, y ≥ 0M取一个足够大的数比如1e6但不要太大以免数值爆炸。这样即使原问题不可行松弛子问题也有最优解对偶乘子照样可用可行割的推导方式和最优割统一。5.2 割平面数值不稳定把系数控制住我遇到过割平面系数达到1e7量级主问题每次都在找各种各样的极端整数解迭代十几轮还不收敛。主要原因是对偶乘子太大叠加到割上之后数值尺度失衡。解决思路有几个第一把量纲统一比如费用单位从“元/年”换成“万元/年”割系数立刻小两个数量级第二给割约束加一个上限例如η ≤ 1e6避免主问题在无界处乱跳第三对充放电、SOC这类变量做归一化。另外一个容易被忽略的点可行性割和最优割要分开存放不要全都塞进同一个cut数组。因为最优割里有η变量可行性割没有。如果不小心把一个不含η的可行性割写成η≥0的形式模型会莫名其妙不能收敛。5.3 收敛慢时的几个有效“提效”操作如果迭代超过50轮还在磨可以按优先级检查四件事初始化UB是不是太差先跑一次全选方案或贪心方案给UB一个合理初值。主问题每次解出的x是否在震荡如果是就在主问题里加一个no-good割排除掉已经验证过的较差整数组合。子问题是不是被求解器提前终止了子问题必须解到最优否则产生的对偶乘子不对割质量会很差。给子问题设置最优性容差为1e-9。是不是没有充分利用典型日并行多个典型日的子问题可以并行求解切割一次回传大幅缩短单轮耗时。5.4 从确定性规划到不确定性规划GBD天然适配两阶段随机规划IES规划正好常用。把典型日换成多场景加上概率权重主问题仍是选型决策子问题变成每个场景的运行优化每个场景都可以生成一组Benders割。这也是为什么很多文献用GBD做综合能源系统随机规划而不是用单纯MILP。如果要做分布鲁棒优化可行性割还会承担更多职责但整体框架不用推翻。这里要提醒一句如果不是线性或凸模型GBD不能保证全局最优。不要为了省事把非凸汽轮机效率曲线直接扔给GBD。先分段线性化或者做凸松弛否则你得到的所谓“最优解”可能连局部最优都算不上。结尾小记我在实际跑这个算例时最大的体会是GBD在IES规划里的价值不只在省内存而在“可解释性”主问题告诉你应该买什么设备子问题告诉你这样跑要花多少运行成本两者交替非常直观。调试过程中真正起作用的不是高深的理论而是一次又一次老老实实检查初始化、对偶乘子方向和割的数值尺度。最后再分享一个小技巧把每次迭代记录下来的LB、UB和gap画在同一个图上如果发现LB长期不涨先别怀疑算法去看主问题的割是否有重复或数值异常如果UB来回跳动就去给主问题加no-good割。这个小动作能帮你省掉至少一半的调试时间。
