广义Benders分解求解综合能源系统优化规划：建模、算法与Matlab实现

这几年我前前后后接触过几轮综合能源系统规划项目最深的体会是真正难的不是把风、光、气、电、热放在同一个模型里而是模型写完之后求解器不给力。设备选型、容量配置这类决策天生带了0-1整数变量运行调度又需要连续变量两套变量耦合在一起形成一个大规模混合整数规划。此前某课题组做一个园区算例30个候选设备、8760小时时序运行约束直接用通用求解器跑分支定界树膨胀得让人绝望。后来改用广义Benders分解法把原问题拆成投资主问题和运行子问题计算时间一下子压到了可接受范围。这篇文章围绕“基于广义Benders分解法的综合能源系统优化规划”这个主题完整复盘一遍这类问题的处理流程从综合能源系统优化规划的数学建模到广义Benders分解的迭代逻辑再到Matlab代码实现和调试中真正踩过的坑。适合刚接触能源系统优化、且准备自己写分解算法而不是只会调包求解器的朋友。文章里给的代码是我在实际项目中反复改过之后的框架可以直接拿去做教学算例也能作为工程项目的起点。1. 先看清问题综合能源系统规划模型到底在优化什么1.1 投资变量与运行变量两种时间尺度的一次性耦合综合能源系统优化规划核心任务是在一堆候选设备里选出需要建设的电、热、冷、气设备确定它们的容量同时安排它们在典型日里的运行方式。因此模型里天然存在两类变量一类是设备是否投建、投建容量等级的0-1整数变量另一类是每个时刻设备出力、购能功率、储能充放电功率这类连续变量。两者通过设备容量上下限等约束绑定在一起。一个典型的最小化年总费用模型可以写成[ \min \quad F \sum_{i \in \mathcal{I}} \alpha_i c_i y_i \sum_{s \in \mathcal{S}} w_s \sum_{t \in \mathcal{T}} \left( c_{t}^{e} P_{t,s}^{buy} c_{t}^{g} G_{t,s}^{buy} \right) ]其中 (\mathcal{I}) 是候选设备集合(y_i) 表示第 (i) 个设备是否投建(\alpha_i) 是年化投资折算系数(c_i) 是设备单位容量投资成本(\mathcal{S}) 是典型日集合(w_s) 是场景权重(\mathcal{T}) 是每个典型日内的时间段集合(P_{t,s}^{buy}) 和 (G_{t,s}^{buy}) 分别是从电网购电和购买天然气的功率。这里要注意如果单独看投资部分这是一个典型的“选不选、选多大”的问题单独看运行部分又是一个带时序约束的调度问题。真正的麻烦在于投资决策要先于运行调度做出但运行调度又反过来决定运行成本和设备容量是否满足需求。这种“先决策、后评估”的结构是后面用分解法的直接原因。1.2 约束条件能量平衡、设备特性、储能状态一个都不能少模型约束大致分成四类。首先是电、热、冷、气等能量平衡约束比如电力平衡可以写成[ \sum_{i \in \mathcal{I}{elec}} P{i,t}^{out} P_{t}^{buy} P_{t}^{PV} L_{t}^{elec} P_{t}^{sell} \sum_{i \in \mathcal{I}{elec}} P{i,t}^{self} ]意思是各类电源出力、外购电力、可再生能源出力加总要满足电负荷以及自身消耗和向外售电。热、冷、气同理每个母线上都得有相应的平衡方程。第二类是设备物理特性约束最常见的是出力上下限[ y_i P_i^{\min} \le P_{i,t} \le y_i P_i^{\max} ]当 (y_i0) 时设备不出力当 (y_i1) 时出力被限制在最小和最大技术出力之间。如果同时考虑爬坡还要加上相邻时段出力差值的约束。第三类是储能设备的状态约束例如[ E_{i,t1} E_{i,t} \eta_i^{ch} P_{i,t}^{ch} - \frac{P_{i,t}^{dis}}{\eta_i^{dis}} ]储能末状态往往要求等于初状态保证一个运行周期内能量自洽。第四类是网络约束。如果是冷热电联供园区还会包含热网管道的容量、电压或潮流、天然气管道流量上限等。这类约束会让模型更贴近实际但也会显著增加求解难度。做规划时可以在初版先忽略网络或者用线性化网络等分解框架跑通后再加入非线性项。1.3 难在哪里组合爆炸和时序耦合一起出现有了上面的模型结构直接丢给求解器会怎样对于 (I) 个候选设备如果每个只有“投”和“不投”两种状态理论上就有 (2^I) 种组合。通用MILP求解器采用分支定界法在搜索树每个节点上解一个线性规划。当设备数量超过20、运行时段达到几十甚至上百时节点数量可能指数级膨胀并且每个节点上的LP还带着整段时间序列的耦合约束。我见过很多人第一反应是“用Gurobi或CPLEX直接算不就行了吗”但实际上商业求解器往往只适合中小规模的算例。在多典型日、多能源耦合、含有网络约束的工程算例里求解器经常会长时间卡在分支定界树的中间层输出gap迟迟不降。这时候自然想到能不能不盲目搜索所有投资组合而是通过迭代用少量关键组合逼近最优解广义Benders分解法解决的就是这个问题。2. 广义Benders分解的核心思路与迭代链条2.1 经典Benders和广义Benders的差别经典Benders分解最初是为求解大型混合整数线性规划设计的。它把所有变量分成“复杂变量”和“简单变量”两类复杂变量通常是整数或0-1变量固定它们之后剩余的子问题变成一个相对容易的线性规划。然后利用线性规划对偶理论把子问题的最优值信息“投影”回复杂变量空间形成一组割平面逐步逼近原问题最优解。广义Benders分解是经典版本的推广。经典Benders要求子问题是线性规划而广义Benders只要求子问题是凸优化问题甚至可以是凸二次规划、二阶锥规划等。在综合能源系统规划中如果子问题里的非线性约束经过凸松弛或分段线性化广义Benders分解就可以直接使用。如果子问题是非凸的那只能把Benders当成启发式算法不能保证收敛到全局最优。这也是为什么很多教程都强调用GBD前先确认子问题的凸性。综合能源系统里最常见的非凸来源是天然气管道流量方程、热网节点温度和压力耦合方程这些约束通常可以写成凸松弛形式。工程上比较稳妥的做法是先用线性化版本跑通框架再逐步加入更精确的非线性模型。2.2 主问题、子问题与割平面三个角色的协作机制把原问题抽象地写成[ \min_{y} ; C^{cap}(y) G(y) ]其中 (y) 是投资变量(C^{cap}(y)) 是投资成本(G(y)) 是在给定投资方案 (y) 下运行子问题的最优运行成本[ G(y) \min_{x} ; C^{op}(x,y) \quad s.t. \quad g(x,y) \le 0 ]固定 (y) 后子问题是一个只包含连续变量 (x) 的优化问题。如果子问题有解假设目标函数是凸的约束也是凸的那么可以通过拉格朗日乘子构造一条Benders割[ \theta \ge G(y_k) \lambda_k^T (y - y_k) ]这条割的含义是运行成本的下界可以表示成一个线性函数斜率由对偶乘子决定。如果子问题在当前 (y_k) 下不可行就需求解一个最小化不可行量的可行性子问题得到另一组乘子生成可行性割把不可行的投资方案排除掉。主问题则变成[ \min_{y, \theta} ; C^{cap}(y) \theta ]约束是所有已经生成的割平面。主问题的目标函数值给出了原问题的一个下界因为割平面是对 (G(y)) 的外部逼近而每次子问题求出的“投资成本运行成本”则给出原问题的一个上界。上下界不断收缩直到它们之间的相对差满足容忍度。迭代流程可以用一句话概括子问题负责评估给定投资方案的真实成本和对偶信息主问题负责根据这些信息找到更有希望的投资方案。两者轮流求解互相传递信息直到收敛。2.3 一个生活类比食堂备菜计划把Benders分解理解成一个做菜计划问题。主问题是食堂经理他决定未来一周的“菜单模板”投资方案比如哪些窗口开放、哪些食材按什么标准采购。子问题是每天的后厨固定菜单模板后后厨需要安排每个时段的备菜、烹饪和售卖节奏运行调度并计算出人力、能耗成本。如果后厨发现今天某种食材采购量太少导致不能出菜就会反馈“这个菜单模板行不通”要求经理调整。如果菜能做出来后厨还会告诉经理“如果把某类食材的采购上限提高一些每天的运营成本还能再降多少”这个反馈就是Benders割。经理根据这些反馈重新调整菜单模板再交给后厨试跑。经过若干轮菜单模板和每日成本都逐渐稳定在最优解附近。这个比喻虽然简单但“评估—反馈—调整”三个动作和GBD完全对应。3. Matlab实现拆到代码级别的可复现方案3.1 工具选型YALMIP加外部求解器还是自带优化工具箱Matlab本身提供的intlinprog和linprog可以用于小规模教学算例。但在GBD框架里主问题要反复加入割约束并重新求解子问题也要不断提取对偶信息用自带工具箱需要手工维护稀疏矩阵非常容易出错。实际项目中我更推荐用 YALMIP 建模后端接 Gurobi 或 CPLEX。YALMIP 的好处在于约束和目标的写法非常接近数学表达式可以在循环中动态添加约束并且可以直接用dual()拿到约束对应的拉格朗日乘子。主问题里的0-1变量用binvar定义连续变量用sdpvar定义割约束本质上就是一些线性不等式动态添加很直观。后端求解器方面Gurobi 和 CPLEX 的稳定性比开源求解器好得多尤其在数值尺度复杂的能源系统模型里开源的glpk很难保证对偶信息的准确性。下面我给出一个相对通用的代码骨架它包含主循环、子问题求解、主问题求解三个部分。我的习惯是先把整个流程写成一个函数输入是封装好的IES数据体输出是最优投资方案和成本上界。这样后期换算例、换数据都很方便。3.2 主循环上下界迭代的骨架代码function [best_y, best_obj] run_gbd(IESdata) tol 1e-4; max_iter 200; UB 1e8; LB -1e8; y_k ones(IESdata.nDevice, 1); % 初始方案先全投再看是否可行 cuts struct(type,{}, y,{}, dual,{}, val,{}); gap 1.0; for k 1:max_iter % 1. 子问题固定 y_k求解运行优化 [sp_obj, is_feas, lambda] solve_subproblem(IESdata, y_k); if is_feas ub_now compute_cap_cost(IESdata, y_k) sp_obj; if ub_now UB UB ub_now; best_y y_k; end cuts(end1) struct(type,opt, y,y_k, dual,lambda, val,sp_obj); else % 子问题不可行生成可行性割 mu solve_feasibility_subproblem(IESdata, y_k); cuts(end1) struct(type,feas, y,y_k, dual,mu, val,0); end % 2. 主问题求解下界和新的投资方案 [y_new, LB] solve_master_problem(IESdata, cuts); % 3. 更新迭代数据 y_k y_new; gap abs(UB - LB) / max(1, abs(UB)); fprintf(iter%3d, UB%.4f, LB%.4f, gap%.6f\n, k, UB, LB, gap); if gap tol break; end end end这里有几个容易错的地方。第一子问题的目标函数不要包含与固定变量无关的常数项最好把运行成本单独返回投资成本在主问题或上界计算里再加。第二如果初始方案全投仍不可行说明原问题本身可能无解要先检查约束和数据是否一致。第三UB不是简单取sp_obj而是“固定投资方案下的投资成本运行成本”。我在第一次写代码时就是在这里把上界算错了导致迭代一直不收敛。3.3 子问题求解与对偶信息提取子问题通常写成YALMIP模型function [obj_val, is_feas, lambda] solve_subproblem(IESdata, y_current) x sdpvar(IESdata.nContVar, 1); constraints []; % 设备出力上下限受 y_current 控制 for i 1:IESdata.nDevice constraints [constraints, ... y_current(i) * IESdata.MinPower(i) x(i), ... x(i) y_current(i) * IESdata.MaxPower(i)]; end % 能量平衡等约束 constraints [constraints, IESdata.Aeq * x IESdata.beq]; options sdpsettings(solver, gurobi, verbose, 0); objective IESdata.Cost * x; sol optimize(constraints, objective, options); if sol.problem ~ 0 is_feas false; lambda []; obj_val []; return; end is_feas true; obj_val value(objective); % 提取第一个等式约束对应的对偶乘子 lambda dual(constraints(1)); % 实际使用时要根据约束顺序逐个提取而不是只取第一个 end需要注意的是YALMIP 的dual()返回的乘子符号与约束形式有关。如果约束写成A*x b乘子非负如果写成A*x b乘子非正。Benders割里的系数符号直接受这个影响。我习惯在写子问题时统一把不等式整理成A*x b的形式这样对偶乘子就是非负向量构造割的时候不容易搞混。另外如果你的子问题里面用了逻辑约束、if-else、implies这类YALMIP高级表达提取对偶信息就会很麻烦因为它们会引入辅助变量并改变对偶结构。做GBD时子问题尽量保持单纯的线性或凸二次模型所有整数变量都放进主问题。3.4 主问题中割约束的组装方式主问题每次迭代都会新增一条割所以代码里要用一个结构体数组cuts保存历史割。组装割约束时最优性割写成function [y_new, LB] solve_master_problem(IESdata, cuts) y binvar(IESdata.nDevice, 1); theta sdpvar(1, 1); constraints []; % 设备容量逻辑约束 constraints [constraints, IESdata.CapCon * y IESdata.CapRhs]; % 遍历所有割 for i 1:length(cuts) yk cuts(i).y; dual_k cuts(i).dual; val_k cuts(i).val; if strcmp(cuts(i).type, opt) % 最优性割: theta val_k dual_k * (y - yk) constraints [constraints, theta val_k dual_k * (y - yk)]; else % 可行性割: 0 dual_k * (y - yk) 或类似形式 constraints [constraints, 0 dual_k * (y - yk)]; end end objective compute_cap_cost(IESdata, y) theta; options sdpsettings(solver, gurobi, verbose, 0); optimize(constraints, objective, options); y_new value(y); LB value(objective); end可行性割的写法比最优性割更灵活。子问题不可行时我们通常先求解一个带非负松弛变量的最小化不可行问题然后取松弛变量对应的对偶乘子。这个乘子的物理意思是“当前投资方案下哪些约束稍微放松一点就能让运行可行”。把它作用到投资变量上得到的线性约束会让主问题避开不可行区域。如果直接套用最优性割的形式反而可能生成错误的割这一点务必小心。4. 调试实录不可行子问题、收敛振荡和数值尺度4.1 子问题不可行第一个拦路虎做GBD的人第一周基本都在跟“子问题不可行”搏斗。最常见的原因是初始投资方案没有给足容量导致某些时段能量平衡无解。比如只有一个电锅炉和一个燃气轮机电负荷高峰同时还要满足热负荷但燃气轮机的爬坡速率不够就会出现无法同时满足电平衡和热平衡的情况。处理不可行问题的标准做法是构造可行性子问题。给所有硬约束加上非负松弛变量最小化松弛量的和[ \min ; \sum_{l} s_l^ s_l^- \quad s.t. \quad A x s^ - s^- \ge b(y_k), ; s^\pm \ge 0 ]这个问题的目标就是在当前投资方案下看看至少要“破坏”多少需求才可行。求解后取松弛变量对应的对偶乘子就可以生成可行性割。实际操作中我不建议给所有约束都加松弛最好只给最核心的能量平衡约束加否则可行性割会变得太平淡主问题很难有效调整方案。还有一个更容易忽略的情况储能约束的初末电量不等导致不可行。设计模型时如果设置了一个24小时周期且要求末电量等于初电量那么初始方案中储能容量不足或者负荷模式偏斜时子问题就会无解。此时可以把储能初末电量写成一个软约束允许微小偏差并计入运行成本这样至少能保证子问题始终有解GBD迭代也能正常启动。4.2 上下界不收敛问题不一定出在算法可能在模型调试GBD时最令人头大的现象是每个子问题都能求解主问题也能求解但上下界之间就是有一块怎么都缩不回去的gap。这种时候先不要怀疑算法大概率是模型本身有问题。第一个常见原因是主问题中的辅助变量 theta 缺少下界。若前期割比较少theta可能被主问题取成很大的负值导致下界低得离谱。解决办法是给 theta 加一个物理下界比如theta 0甚至用原问题LP松弛的最优值作为更紧的下界。初始下界越紧迭代次数越少。第二个常见原因是割平面太弱。Benders割虽然是有效割但不一定是Pareto最优割。如果子问题的最优解不唯一随便取一个对偶乘子得到的割可能非常弱导致边界更新缓慢。工程上有一个实用技巧在解完最优子问题后再解一个“正则化”的目标比如在原目标基础上加上一个很小的二次项让解更稳定。很多商业求解器默认的对偶信息实际上不稳定用正则化可以让割更有效。第三个原因是数值尺度过差。投资成本的量纲可能是“万元”运行成本是“元”一个设备功率用“MW”另一个用“kW”这种数量级的混搭会让割平面的系数差好几个数量级。体现在收敛曲线上就是上下界的相对gap死活降不到1e-4以下。解决方式是统一单位费用统一到万元功率统一到kW时间统一到小时最后再检查目标函数的绝对数量级。4.3 常见现象速查表现象可能原因处理建议子问题第一次求解就无解初始投资方案容量不足或数据单位不统一先用启发式生成一套可行方案再启动GBD无休止迭代但gap不下降割平面太弱或theta缺少下界加theta初始下界尝试Pareto最优割主问题求解非常慢累积割数量过多模型膨胀定期删除冗余割或使用多切割Benders对偶乘子符号与预期不一致约束写成形式导致乘子符号反向统一把不等式整理成并验证符号UB和LB出现振荡Benders割系数差数量级或求解器数值不稳定统一单位、放宽求解器容忍度、使用双精度计算主问题得到的y方案始终不变割没有真正限制投资变量检查可行性割生成方式确认割作用在正确变量上这张表是我在实际项目里积累出来的不是教科书里的标准答案。很多问题并不是孤立的而是多个原因叠加在一起。遇到不收敛时我习惯先看一条迭代过程中的“UB/LB曲线”如果UB单调下降而LB几乎不动那多半是割太弱如果UB和LB都在波动那就要回头查对偶符号和数值尺度。5. 从基础框架到工程化扩展5.1 多典型日和多场景GBD天生适合并行规划模型里很少只用一个典型日通常是冬季、夏季、过渡季各选若干个典型日再按天数加权。固定投资方案后每个典型日的运行子问题相互独立可以并行求解。比如有10个典型日那每轮迭代就能同时产生10条子问题结果如果采用多切割Benders每个典型日对应一个独立的theta变量主问题目标变成“投资成本各典型日运行成本加权之和”。我曾经在一个风、光、气、热多能互补算例里用了24个典型日。直接构建单层MILP求解内存占用高得惊人用GBD后子问题按场景并行每轮迭代的计算量反而比完整模型的单次松弛求解还小。需要注意并行时要小心随机数种子和求解器线程配置保证同一场景在不同的迭代轮次里结果可复现否则收敛条件会乱掉。多切割Benders的好处是每个典型日的运行成本互不干扰割平面数量变多但每一条割都更紧。缺点也很明显主问题变量和约束数量会随典型日数量线性增加。如果典型日特别多比如几百个场景聚合为一条期望割反而更稳定。具体选哪种取决于主问题的求解时间占比。5.2 与随机规划和鲁棒优化的结合点综合能源系统规划里负荷和新能源出力天然带不确定性。常用处理方式是把不确定性建模成多场景构成两阶段随机规划第一阶段是投资决策第二阶段是每个场景下的运行决策。这个结构与GBD的分工几乎完全一致——投资变量作为第一阶段变量场景运行变量作为第二阶段变量。做随机扩展时在第2节的主循环里需要增加的只是一个“场景维”子问题从单场景变成多场景。如果采用机会约束或分布鲁棒约束则还需要在割平面之外额外加入与风险度量相关的约束通常和Benders割结合在一起。我曾经尝试过将条件风险价值约束写入主问题效果还可以但主问题会从一个纯MILP变成一个带二阶锥约束的MIP求解速度明显下降。还有一个方向是鲁棒优化。两阶段鲁棒优化常用列约束生成算法而不是Benders分解。但如果把不确定集表示成有限离散场景列约束生成和GBD在某种意义上是等价的。理解了GBD之后再去读列约束生成的代码几乎不会觉得费力。5.3 非线性能源网络如何处理综合能源系统里最麻烦的非线性来源是热网管道的流量方程和天然气网的节点压力方程。天然气管道流量与两端压力平方差之间存在非线性关系通常在工程上会被改写成二阶锥约束或分段线性近似。如果采用二阶锥约束子问题仍然是凸的广义Benders分解在理论上仍然适用。不过现实没有这么顺利。热网节点温度与流量之间有时会包含 bilinear 项这种约束在一般情况下是非凸的。某课题组在初期直接把非凸约束放进子问题再用GBD求解结果迭代数次后上下界明显振荡收敛到了一个非常奇怪的局部解。后来他们把热网水力方程线性化温度动态过程保留在凸约束范围内模型才稳定下来。这个教训说明了GBD是一个“凸优化友好”的框架非凸模型不能指望它神奇地找到全局最优。如果非要保留非精度较高的非凸子问题我建议把GBD当作一种“热启动工具”。先用线性化模型跑出一个种子解再用这个种子解去初始化原非凸模型的求解器这样通常能得到比单纯用求解器从零开始更快的计算速度。最后分享两点个人经验一是调试GBD时一定先把子问题的对偶信息打印出来和有限差分或用扰动法验证一下很多时候程序“跑得通”只是表象割的符号错了照样收敛不了二是养成画UB/LB曲线的习惯曲线几乎能告诉你所有调试方向。我第一次写GBD代码时就是靠这两条经验把一个原本收敛到0.05误差都费劲的算例调到了稳定收敛在1e-5以内。这套框架到现在我已经用在好几个园区类的规划模型上依然很顺手。