★ FEATURED ARTICLE 最近在做一个综合能源系统的规划项目内容本身不算复杂——给一个小型园区设计电、热、气三种能源的供能方案设备无非是CHP机组、电锅炉、燃气锅炉、光伏和储能。真正让我头疼的是求解效率。第一次把整个问题建成一个混合整数规划模型直接丢给求解器等了整整一天内存快要吃光最优解还是遥遥无期。后来换成了广义Benders分解Generalized Benders DecompositionGBD重新组织求解流程同样的问题几十次迭代收敛耗时从那一天的绝望变成了一顿饭的功夫。这篇博文就围绕这个项目展开把我从建模到代码实现、再到排查求解坑的完整过程梳理出来。内容包括综合能源系统规划问题为什么难解GBD到底在做什么数学模型怎么拆成主问题和子问题Matlab代码的核心框架长什么样以及一组能跑通的小型算例结果。适合正在做能源系统优化、学过优化算法但没真正动手写过分解算法的同学也适合那些遇到大规模混合整数规划求解困难、想换思路的工程师。1. 解这道题之前先拆清楚它到底在求什么1.1 综合能源系统规划问题的天然结构综合能源系统优化规划说白了就是回答两个问题建什么设备、建多大容量以及建成之后每天怎么运行。前者是投资决策后者是运行调度决策。这两个问题不是独立的——你选的设备容量决定了运行阶段能调度多少资源而运行成本反过来又影响投资方案是否划算。典型的小型园区系统里供能侧可能有CHP机组热电联产、燃气锅炉、电锅炉、光伏板、电储能、热储能有些还会配置电制冷机和吸收式制冷机。需求侧则有电负荷、热负荷、冷负荷某些场景还有气负荷。能量流的方向大致是天然气进入CHP和燃气锅炉CHP同时出电和出热光伏出电电锅炉把电变成热储能负责跨时段的能量平移。对接真实工程时规划周期一般取一年但运行层不能按一整年单时段算因为电负荷、热负荷、光伏出力都有明显的日内波动和季节性差异。常见做法是把一年压缩成几个典型日每个典型日按小时或15分钟一个时段来建模。这样运行层就变成一个多场景、多时段的优化问题。1.2 为什么直接求解会卡死你可以把这个规划问题写成一个混合整数线性规划MILP或者混合整数非线性规划MINLP然后把所有变量一股脑扔给求解器。在项目规模很小时这没问题比如一个单节点、单典型日、设备只有两三种的简化模型求解器秒出结果都正常。但问题一旦扩大组合爆炸就来了。第一投资层有整数变量。每台设备的选型、台数、容量档位都是整数决策比如某个设备有5个可选容量档位一个系统里有8类设备组合数就是5的8次方。虽然有求解器的分支定界法在剪枝但分支树在最坏情况下仍然可能指数增长。第二运行层把问题维度撑大了。假设你有4个典型日每个典型日24个时段每时段有设备出力、储能充放、购能、弃能等几十个连续变量再加上几百条平衡约束运行子问题的规模就已经不小。而投资层每试探一个方案就要把所有典型日的运行调度全部算一遍。第三也是最关键的整数变量和连续变量纠缠在一起。分支定界法在搜索过程中反复把整数变量分为0和1两个分支每走一个分支就要解一个大线性规划。如果这个线性规划本身有几万个变量和约束分支节点一多总计算量就是天文数字。当时我在做这个项目的模拟算例时8类设备、4个典型日、每个典型日24个时段直接求解的MILP大概有12万个变量其中整数变量不到100个但求解器跑了几个小时连一个可行解都没找到。原因就是松弛后的线性规划太大分支节点上每次求LP都很慢。那时候我就意识到必须把问题拆开处理——该分层的分层该迭代的迭代。2. 解法选型的思路广义Benders分解为什么合适2.1 GBD的直观逻辑Benders分解的思想其实不复杂既然整数变量是麻烦的根源那就把整数变量和连续变量分开看待用迭代的方式让它们交替决策。想象一个工程项目负责人要决定买哪些设备整数决策买完之后运行部门根据设备配置去安排生产计划连续决策。负责人不可能一开始就知道哪个设备方案最好他只能先猜一猜然后让运行部门按照这个设备方案给出运行成本同时还要让运行部门告诉他如果你改变设备容量成本变化的方向和速度是怎样的。运行部门给出的这个反馈信息就是Benders割。负责人看了反馈修改方案再问运行部门如此循环直到方案成本不再改善。这个思路在数学上可以严格证明收敛到全局最优解前提是问题满足一定条件——例如连续子问题是凸的整数变量只以线性形式出现在耦合约束中。综合能源系统规划模型正好满足这些条件所以GBD是一个很自然的候选解法。2.2 与其他分解算法的对比可能有人会问分解算法那么多为什么选GBD而不是拉格朗日松弛或者ADMM交替方向乘子法这里有一个很实际的工程经验总结我放在下表里特性GBD拉格朗日松弛ADMM适用场景整数变量连续变量的两阶段结构约束复杂、含整数变量的问题大尺度凸问题、目标可分离结构迭代机制子问题给主问题传割平面乘子更新逼近对偶最优原始-对偶交替更新对整数变量处理保留在主问题中精确求解需启发式恢复可行性需额外处理离散变量收敛速度问题规模合适时收敛快乘子震荡收敛慢线性收敛高精度下较慢实现复杂度中割生成逻辑需细心低中低适用性判别子问题连续凸即可需要可松弛结构需要目标可分离ADMM在分布式优化里非常火但它在处理整数变量时需要额外的映射操作对规划这种强组合问题不如GBD直接。拉格朗日松弛收敛慢而且对偶间隙处理麻烦对工程场景不够稳定。GBD的收敛性在MILP结构下有明确的理论保障实现起来又有清晰的迭代框架所以最终选了它。下面再细说数学模型的具体形式。3. 数学模型与分解逻辑的细节3.1 两阶段规划模型的标准形式先写一个抽象但足够通用的两阶段模型后面所有实现都围绕这个形式展开。第一阶段是投资决策[ \min_{x, \eta} \quad c^T x \eta ][ \text{s.t.} \quad A x \le b, \quad x \in {0,1}^{n_x} ]其中 (x) 是设备投资决策的0-1变量(c^T x) 是投资等年值成本(\eta) 是一个辅助变量表示对第二阶段运行成本的估计在后面的迭代中会被Benders割逐步约束。第二阶段是运行调度给定一组投资变量 (x \hat{x})[ Q(\hat{x}) \min_{y} \quad d^T y ][ \text{s.t.} \quad W y \le r - T \hat{x}, \quad y \ge 0 ]这里 (y) 是每个典型日中每个时段的连续运行变量包括各设备出力、储能充放功率、购电量等。约束 (Wy \le r - T\hat{x}) 把所有投资变量和运行变量耦合起来实质上是在表达只有装了设备设备才能出力以及设备容量限制出力上限。(Q(\hat{x})) 就是在给定投资方案下的最优运行成本。整个问题的最终目标就是让投资成本和运行成本之和最小。写成这种形式的好处是数学上可以直接套GBD的标准框架不用对每个具体模型重新推导。实际项目里的设备模型、能流平衡、爬坡约束全部装进上式的 (A)、(W)、(r)、(T) 矩阵里就可以。3.2 主问题与子问题的划分策略GBD的核心是把上面的两阶段模型拆成两个部分。主问题只包含整数投资变量 (x) 和辅助变量 (\eta)再加上一组不断累积的Benders割约束。它的形式是[ \min_{x, \eta} \quad c^T x \eta ][ \text{s.t.} \quad A x \le b, \quad x \in {0,1} ][ \eta \ge \alpha_k \beta_k^T (x - x_k), \quad k1,2,...,K ][ 0 \ge \gamma_l \delta_l^T (x - x_l), \quad l1,2,...,L ]第一类割叫最优性割当子问题可行且求得最优解时生成作用是收紧对运行成本 (\eta) 的估计。第二类割叫可行性割当子问题不可行时生成作用是砍掉那些导致运行不可行的投资方案。子问题则是固定投资变量 (x \hat{x}) 后求解连续运行问题也就是上面提到的 (Q(\hat{x}))。它的目标是得到运行成本精确值和拉格朗日乘子乘子 (\lambda) 正是计算Benders割的关键材料。3.3 两类Benders割的数学来源很多人写代码时卡在割的推导上这里我展开说一下。先看最优性割。固定投资方案为 (\hat{x})子问题的最优单纯形表中对耦合约束 (Wy \le r - T\hat{x}) 对应的对偶变量拉格朗日乘子记为 (\lambda_k)。根据线性规划对偶理论子问题的最优值可以写成一个关于 (\hat{x}) 的线性函数[ Q(\hat{x}) \lambda_k^T (r - T\hat{x}) ]注意这个等式只在 (\hat{x}) 附近严格成立在其他位置它是一个下界。所以对任意的 (x)有[ Q(x) \ge \lambda_k^T (r - T x) ]把右边展开[ Q(x) \ge \lambda_k^T (r - T\hat{x}) \lambda_k^T (T\hat{x} - T x) ]第一项等于 (Q(\hat{x}))第二项是线性修正项。这样我们就得到了最优性割[ \eta \ge Q(\hat{x}) (T^T \lambda_k)^T (x - \hat{x}) ]记 (\alpha_k Q(\hat{x}))(\beta_k T^T \lambda_k)就对应前面主问题里的第一条割约束。再看可行性割。如果某个投资方案导致子问题无解我们需要生成一条割把这个不可行的投资方案从主问题可行域中切掉。做法是修改子问题引入松弛变量求解一个最小化约束违反量的辅助问题[ \min_{y, s} \quad \mathbf{1}^T s ][ \text{s.t.} \quad W y s \le r - T x, \quad s \ge 0, \quad y \ge 0 ]如果最优值大于0说明原问题在当前投资方案下不可行。同样取这个辅助问题中耦合约束的对偶乘子 (\mu_l)可以得到可行性割[ \mu_l^T (r - T x) 0 ]等价地写成主问题约束[ 0 \ge \gamma_l \delta_l^T (x - x_l), \quad \gamma_l \mu_l^T (r - T x_l), \quad \delta_l T^T \mu_l ]这里 (x_l) 是导致不可行的那个投资方案。这条割加入主问题后主问题会避开不符合运行要求的投资方向。4. Matlab代码实现完整流程与关键点4.1 代码结构总览我的代码框架分为六个模块每个模块一个脚本或函数模块作用main主入口初始化参数启动迭代循环build_problem构建所有矩阵、参数、时段数据生成典型日场景solve_subproblem求解运行子问题返回运行成本与乘子generate_cut根据子问题结果生成最优性割或可行性割solve_masterproblem求解主问题返回新的投资方案与下界check_convergence计算上下界间隙判断是否收敛主循环的逻辑可以浓缩成下面的伪代码% 初始化 x_current initial_investment(); ub inf; lb -inf; cuts_opt []; cuts_feas []; for iter 1:max_iter % 求解子问题 [status, q_value, lambda] solve_subproblem(x_current); if status feasible ub min(ub, investment_cost(x_current) q_value); cuts_opt [cuts_opt, generate_optimality_cut(x_current, q_value, lambda)]; else cuts_feas [cuts_feas, generate_feasibility_cut(x_current, mu)]; end % 求解主问题 [x_new, lb] solve_masterproblem(cuts_opt, cuts_feas); % 收敛检查 gap abs(ub - lb) / abs(ub); if gap tol break; end x_current x_new; end4.2 子问题的关键实现细节子问题是整个算法里计算量最大的部分。在Matlab里我使用YALMIP建模求解器用的是Gurobi。子问题的输入是当前的投资变量取值输出是运行成本和乘子。这里有个很重要的经验子问题必须是线性规划或二次规划不能含整数变量。如果子问题还有整数变量拉格朗日乘子就不存在Benders割就没法直接生成。所幸综合能源系统的运行调度在固定设备容量后通常都是连续的储能充放电连续、机组出力连续唯一可能出现整数变量的是机组启停但规划尺度的模型通常把启停变量简化为容量上限约束或者用更粗糙的固定运行策略处理。子问题求解完成后乘子提取是这个环节最核心的操作。以YALMIP为例% 假设model是已经构建好的子问题模型 % 固定投资变量 x_hat assign(investment_vars, x_hat); % 求解 optimize(model_constraints, objective, sdpsettings(solver,gurobi,verbose,0)); % 提取耦合约束的对偶乘子 lambda dual(coupling_constraints);注意这里dual()返回的乘子正负号取决于约束怎么写。我第一次实现时就是吃了正负号的亏得出的割方向完全反了主问题越迭代下界越离谱。建议大家在实现后用一个小算例验证乘子符号固定一个已知最优投资方案检查 (Q(\hat{x}) \lambda^T (r - T\hat{x})) 是否成立如果差一个负号说明乘子或者矩阵里符号出了问题。4.3 主问题求解与割集累积主问题本质上是一个小规模的MILP变量是投资方案 (x) 加上一个辅助变量 (\eta)约束除了投资自身的约束外全部是累积的Benders割。为什么要累积因为每轮割都是对目标函数的一个局部近似单条割不够必须随着迭代逐渐逼近真实运行成本函数。实际项目中最优性割可能有几十条可行性割相对少一些。主问题规模小求解速度非常快实测每次主问题求解时间在几秒以内完全不是瓶颈。主问题的YALMIP建模大致长这样% x是投资变量binvareta是连续辅助变量 x binvar(n_device, 1); eta sdpvar(1, 1); Constraints []; % 投资变量自身约束例如至少选择一台、总预算限制 Constraints [Constraints, A_ineq * x b_ineq]; % 累积最优性割 for k 1:size(cuts_opt, 2) alpha_k cuts_opt{k}.alpha; beta_k cuts_opt{k}.beta; x_k cuts_opt{k}.x; Constraints [Constraints, eta alpha_k beta_k. * (x - x_k)]; end % 累积可行性割 for l 1:size(cuts_feas, 2) gamma_l cuts_feas{l}.gamma; delta_l cuts_feas{l}.delta; x_l cuts_feas{l}.x; Constraints [Constraints, 0 gamma_l delta_l. * (x - x_l)]; end Objective c. * x eta; optimize(Constraints, Objective, sdpsettings(solver,gurobi,verbose,0));这里提一个容易踩的坑主问题求解得到的目标值就是当前的下界LB但不能直接把主问题最优解对应的 (x) 拿来当最终方案因为主问题里的 (\eta) 只是运行成本的估计只有在收敛时才是真实运行成本。最终可行方案必须从子问题得到。4.4 收敛判据和终止条件收敛判据我采用的是经典的上下界间隙[ \text{gap} \frac{UB - LB}{|UB|} ]当gap小于设定阈值比如 (10^{-3}) 或 (10^{-4}) 时终止迭代并以当前上界对应的方案作为最终解。上下界的计算要特别小心。上界UB必须来自一个真实可行的方案即投资方案 (x) 能让子问题有解且目标函数值是投资成本加精确运行成本。下界LB来自主问题目标值它是原问题最优目标值的下界。还有一个工程细节在迭代初期主问题给出的 (\eta) 可能是个极大的负数或正数导致下界波动很大。所以在收敛判断时我加了两个保护条件至少迭代5次后才开始算gap避免前期下界未成形时误判。如果gap在连续10次迭代中下降幅度小于设定值且UB已经在一段时间内没有更新就接受当前解提前退出防止迭代陷入死循环。提前退出的判断标准用代码表达是if iter min_iter abs(ub - lb) / abs(ub) tolerance break; end if ub ub_old abs(ub - lb) / abs(ub) looser_tolerance break; end5. 算例演示小型园区规划5.1 算例设置我在测试阶段使用了一个虚构的小型园区数据规模刻意做得不大方便验证算法正确性。园区配置了4种可选设备CHP机组、燃气锅炉、电锅炉、电储能。可选的容量档位分别为三档、三档、两档、两档所以投资组合的总组合数为 3×3×2×236 种。这个规模足够小可以直接枚举验证GBD是否找到了全局最优。运行层取了2个典型日分别代表冬季和夏季每个典型日24个时段。负荷曲线是人工生成的冬季热负荷高、夏季冷负荷高电负荷两个季节都相对平稳。光伏出力数据取了典型晴天曲线光伏建设容量为固定参数不参与投资决策纯作为出力资源计入子问题。设备参数如下简化版设备容量选项单位投资成本效率运维成本系数CHP机组200/400/600 kW180万元/台·档电效率0.35热效率0.450.02元/kWh燃气锅炉300/500/800 kW80万元/台·档0.90.01元/kWh电锅炉200/500 kW60万元/台·档0.950.01元/kWh电储能100/200 kWh50万元/台·档充放效率0.920.005元/kWh购电电价取了分时电价峰段1.1元/kWh平段0.7元/kWh谷段0.35元/kWh。天然气价格固定为2.6元/立方米按热值折算。5.2 迭代收敛过程用GBD算法跑这个算例初始投资方案随机设为全部取最小容量。让我记录一下前几轮的数据迭代轮次主问题下界LB/万元子问题上界UB/万元gap15821286不合理初期跳过2966122421.1%3106311789.8%5110511423.2%8112111462.2%12113411410.6%15113811400.2%可以看到迭代到第15轮时gap低于0.5%我设置的是0.5%容差所以算法在第15轮终止。作为验证我用枚举法把所有36种组合逐一算过去最优方案为CHP取400kW、燃气锅炉取500kW、电锅炉取200kW、电储能取200kWh总等年值成本约1138.7万元。GBD给出的终解与之完全一致。5.3 结果里的信息有哪些值得解读最优投资方案反映了一个很典型的规划结论CHP机组是系统中的核心供能设备因为它同时供电和供热能显著减少外购电和锅炉用气整体能量利用效率高。燃气锅炉作为调峰备用容量不宜过大否则投资沉没成本偏高。电锅炉虽然投资便宜但在谷电条件下才划算所以容量取最小档。电储能200kWh发挥了削峰填谷的作用在谷时充电、峰时放电抵消了一部分高电价时段的外购电成本。运行层的出力曲线还说明了一个问题CHP机组在冬季典型日几乎是满发热运行热负荷高热电比刚好匹配在夏季典型日则主要跟随电负荷运行热出力作为富余热量由热储能或散热处理。如果实际项目中出现大量弃热说明CHP装机偏大模型应该自行优化掉但在本算例中由于热负荷曲线温和没有出现明显弃热。6. 实战调试中的坑和排查记录6.1 收敛太慢怎么办如果GBD迭代几十次甚至上百次都不收敛多数情况是割的质量太差而不是算法本身的问题。我遇到过的典型案例是子问题中存在大量强相关的约束导致拉格朗日乘子具有多重最优解某些极端乘子生成的割会非常弱。针对这个问题有两个经验可参考尽量保证子问题的求解精度。Gurobi默认的LP求解精度是1e-6如果为了提高速度调低了精度乘子误差会直接影响割的有效性。对乘子做归一化处理。我在代码里对不同量纲的乘子乘以一个比例系数使割的系数保持在相近数量级有效降低了主问题的数值困难。还有一个注意点是如果投资变量是连续容量而非0-1选择GBD的收敛性会变差因为连续变量对应的割是逐点近似逼近非线性函数时容易出现锯齿现象。对于连续容量问题可以用分段线性化处理或者改用外逼近Outer Approximation方法更合适。6.2 子问题不可行怎么定位问题在哪子问题不可行是GBD实现中最容易出现的情况尤其是模型数据没核对好时。我的排查顺序是看不可行松弛变量的数值。在辅助问题里哪些约束对应的松弛变量大于0就说明是哪些约束导致不可行。这一步能快速指向问题根源比如某个时段负荷太高而设备容量不够或者储能充放能力与削峰需求不匹配。检查耦合约束的符号。(Wy \le r - T\hat{x}) 的形式下所有投资变量应该只在右边缩小约束域也就是设备的容量越大、右边越宽松。如果你发现某个容量增大会使右边变小说明系数矩阵 (T) 的符号写反了。检查基础负荷数据的数量级。有些时候不是模型错而是负荷数据比设备容量大了一个数量级导致任何方案都不可行。遇到这种情况先把负荷降到设备容量的50%~80%范围跑通流程后再换真实数据。6.3 数值问题乘子正负号和矩阵稀疏性关于乘子符号我前面已经提过这是最容易踩的坑没有之一。在YALMIP里如果约束写成A*x b那么dual返回的非负乘子对应的是约束如果写成A*x b乘子就是非正。我建议在代码注释里明确标注约束形式和乘子方向避免调试时反复怀疑人生。另一个常见问题是矩阵稀疏性导致的求解奇慢。子问题里每个典型日的约束是分块对角结构不同典型日之间没有耦合。在Matlab里我把所有典型日合并成一个大矩阵求解最开始用的是稀疏矩阵存储但YALMIP在构建模型时如果没注意稀疏性也会生成大量无用的零元素。优化方式是把每个场景的子问题拆开单独求解然后并行化乘子分别提取最终汇总生成割。这样既提高了求解速度又保持了模型的物理结构清晰。6.4 一些更隐蔽的bug主问题和子问题用同一套变量名导致覆盖在Matlab的YALMIP中如果子问题和主问题都定义名为x的变量第二次建模时如果不加区分约束会混在一起。我的做法是把变量封装成结构体比如mp.x、sp.y建模时严格区分。典型日权重没乘进目标函数多典型日建模时总运行成本必须乘以每个典型日所代表的天数权重否则成本被低估投资方案会偏向大容量设备因为运行成本看着太便宜。储能初始电量取值错误储能跨日循环约束写成“初始电量等于调度周期末电量”时如果初始电量设为一个随机值会扭曲调度行为。正确做法是把初始电量也作为优化变量或者在边界条件里体现“初始电量结束电量”。从实用主义的角度来看GBD这类方法最大的价值不是取代求解器而是给那些直接求解吃力的中型规划问题提供一条可控的、工程上可接受的求解路径。它把一个大而难的MILP拆成了一串连续LP和一个小型MILP让每一次迭代都有清楚的物理解释调试起来也更有方向感。我个人在实际项目里用过几次之后最大的体会是算法本身并不难写难的是建模时的结构设计和割验证的细心程度。如果你也想在自己项目里试试这个方法建议先抄一个极小规模的算例枚举所有组合验证算法的正确性确认没有正负号这类低级问题之后再逐步扩大算例规模。这套流程走通一次后面遇到更大规模的问题时你心里会踏实很多。
阅读完成 · 觉得有帮助？
