做综合能源系统规划项目的朋友应该都有同感模型越来越大设备从冷热电三联供到储能、光伏、地源热泵越加越多约束条件层层嵌套最后维度冲到几万个变量里面还混着大量0-1整数变量。这时候直接丢给求解器整体硬解不是内存爆掉就是计算时间离谱。我这两年处理园区级综合能源系统优化规划问题时用得最顺手的方案就是广义Benders分解法GBD在MATLAB里用YALMIP建模把投资决策和运行优化拆成两层交替迭代求解程序稳定、收敛可控、改参数也方便。这篇文章就把这个项目的完整思路、数学模型、MATLAB代码实现和我在调试过程中踩过的坑一次性写清楚适合正在做能源系统规划、微电网容量配置、能源站优化设计的硕博生和工程师参考。整套程序的核心思想其实很朴素先把建哪些设备、建多大容量这个投资决策和设备建好之后全年怎么运行最省这个运行决策拆开分别用两个规模小得多的问题去求解再用割平面把两层决策之间的耦合信息传回去反复迭代直到收敛。下面我从问题本质、数学模型、算法原理、代码实现到调优经验一层一层展开。1. 先把这个项目要解决的问题说清楚1.1 综合能源系统规划的本质是什么综合能源系统Integrated Energy System, IES简单地讲就是把电、气、热、冷多种能源在源侧、网侧、荷侧打通通过燃气轮机、燃气锅炉、电制冷机、吸收式制冷机、热泵这类能量转换设备以及电储能、热储能这类缓冲设备实现多能互补与梯级利用。规划问题要回答的核心问题是在给定负荷需求曲线、能源价格和候选设备参数的条件下每类设备到底建不建、建多大容量同时按这个方案建成之后一年下来系统怎么运行经济性最好。这看起来是两个问题实际上是一个问题设备容量选大了投资成本高但运行灵活运行费用能省下来容量选小了投资省了但运行阶段可能被迫高价购电购气甚至出现供能缺口。所以投资决策和运行决策必须放在同一个框架里联合优化这就是规划-运行耦合的由来。1.2 为什么整体建模直接硬解行不通如果当前面说的所有变量全部建模到一起得到的是一个大规模混合整数规划MILP问题考虑非线性效率曲线的话就是混合整数非线性规划MINLP。做一个粗略的规模估算假设一年8760小时、5类候选设备、4种能源载体运行变量的数量级就是8760乘以设备数和能源种类随便一建模就是几万到几十万个变量再加上设备选型的0-1变量和几百个储能状态约束整体模型的复杂度立刻就不一样了。即便用Gurobi或者CPLEX这种商业求解器直接求解这种大规模MILP的耗时可能从几十分钟到几十个小时不等更难受的是很多情况下根本没办法在可接受的时间内找到可行解。我自己测试过一个中等规模园区算例整体建模直接求解内存峰值接近8GB跑了三个小时才收敛到2%的MIP gap基本没法用于参数灵敏度分析。这个问题的结构特点是投资决策变量数量很少撑死几十个但运行决策变量数量巨大成千上万。这种结构天生适合分解——先把少的投资变量定下来再在给定容量下求解运行为主的优化问题运行优化的结果反过来告诉投资层这个方案哪里不好、往哪个方向调能省钱反复迭代逼近全局最优。说白了这就是Benders分解的思想。1.3 广义Benders和经典Benders到底差在哪经典Benders分解最早是解决混合整数线性规划的它要求子问题必须是线性规划这样才能从对偶解中提取割平面。广义Benders分解Generalized Benders Decomposition, GBD是Geoffrion在1972年做的推广——子问题可以是凸非线性规划只要能从子问题的拉格朗日对偶信息中构造出割平面就行。综合能源系统规划里正好用得上这个推广。如果模型里考虑了设备效率随负荷率变化的非线性曲线、储能损耗的非线性项模型就是MINLP此时GBD几乎是唯一既理论严谨又工程可用的分解框架。当然如果像我程序默认版本这样把非线性环节用分段线性近似处理掉模型退化成MILP那用到的其实就是经典Benders的特例但整个代码框架完全复用只是子问题求解器从NLP求解器换成LP/MILP求解器而已。这也是我为什么坚持用GBD而不是其他分解方法比如拉格朗日松弛的原因GBD的收敛性是理论上有保障的而且实现起来就是主问题-子问题交替迭代对工程师非常友好。2. 数学建模怎么把规划问题写成主问题子问题的结构2.1 变量分层与能量集线器建模要在代码里实现分解第一步是建模时就把变量天然分成两层这非常关键。我习惯用能量集线器Energy Hub模型来抽象综合能源系统把系统看成一个多输入多输出的转换节点输入侧是外购电力和天然气中间经过CHP机组、燃气锅炉、电制冷机、热泵等设备转换输出侧是电、热、冷三种负荷储能设备挂在集线器内部作为缓冲。在这个框架下变量分成两块规划变量主问题变量每个候选设备的选型0-1变量 x_i 和连续容量变量 C_i二者通过 C_i ≤ x_i · C_i_max 关联。运行变量子问题变量每个典型时段下各设备的出力 p_{i,t}、储能充放功率和荷电状态 SOC_{s,t}、与外网的交换功率等。一句话总结主问题决定买哪些设备、买多大子问题回答在给定设备容量下全年怎么运行最省钱。建模的时候就要刻意保持这个分层结构不要把投资成本直接和运行成本揉进同一个表达式否则后续分解很别扭。2.2 目标函数投资年化与运行成本怎么算规划问题的目标函数通常取年等值总成本最小包含两部分第一部分是投资年化成本。设备投资不是一次性算进目标就行要按设备寿命和折现率换算成年等值成本用到的系数叫资本回收因子CRFCRF r(1r)^n / ((1r)^n - 1)其中 r 是折现率n 是设备寿命。这个细节很多初学者会忽略直接用总建设成本放进目标函数结果设备寿命不同导致对比失真优化结果明显偏离工程实际。我程序里每个设备都有独立的寿命参数和对应的CRF这是必须的。第二部分是年运行维护成本包括外购电费、天然气费以及设备单位出力的运维费用。如果课题涉及碳约束或者可再生能源配额还会额外加碳排放成本和惩罚项。目标函数总体写出来就是min ∑ (CRF_i · C_i · inv_cost_i) ∑_t ∑_s w_s · (购电成本 购气成本 运维成本)注意这里的运行成本是分场景加权求和后面细说。2.3 约束体系容量、运行、储能三类约束约束条件按层次可以分成四类建模时建议严格分开写方便后面取对偶和排错。第一类设备容量上下限约束每台设备容量要么为0不选型要么在最小可建容量和最大可建容量之间用 C_min · x ≤ C ≤ C_max · x 表达。加最小可建容量很重要否则优化器可能会给你解出一个0.01kW这种毫无工程意义的蚊子腿容量。第二类设备出力与容量的耦合约束任何时段的出力不能超过已配置容量。比如CHP机组电出力 0 ≤ p_chp,t ≤ C_chp热出力通过热电比和电出力耦合。注意这类约束就是后面Benders割平面系数的主要来源所以建模时必须把这些约束的对偶信息单独保留。第三类能量平衡约束这是子问题的核心电平衡外购电 CHP电出力 光伏出力 电池放电 电负荷 电制冷机耗电 电池充电热平衡CHP热出力 燃气锅炉热出力 蓄热罐放热 热负荷 蓄热罐蓄热冷平衡电制冷机制冷 吸收式制冷机制冷 冷负荷气平衡外购天然气 CHP耗气 燃气锅炉耗气第四类储能约束充放功率上下限、SOC状态转移方程、调度周期始末SOC相等约束。储能约束里SOC的初末相等是规划类问题里容易被漏掉的一条不加上会导致储能被白嫖优化结果明显失真。2.4 典型场景日怎么选直接拿8760小时全塞进约束里子问题规模也扛不住。工程上通用的做法是场景缩减用k-means聚类把全年负荷、光伏出力、能源价格曲线聚成若干个典型日比如我常用的组合是4个季节乘以工作日/节假日再乘以典型日类型取12个场景每个场景带一个权重 w_s表示该类型在全年的天数占比。做了场景缩减之后子问题可以进一步拆成每个场景独立的运行优化问题。这个拆分非常重要——它直接让子问题变成多个小规模LP并行求解每一轮迭代几秒钟就能跑完。而且每个场景独立求解后各自生成各自的割平面最后汇总到主问题逻辑上非常干净。3. 广义Benders分解算法迭代逻辑与割平面原理3.1 割平面是怎么来的很多朋友第一次接触Benders分解最迷惑的就是割平面那串公式。我用直白的方式解释一下。把原问题抽象成 min f(x,y)其中 x 是运行变量y 是规划变量。固定 y y* 之后子问题是关于 x 的凸优化问题。求解子问题除了得到最优运行成本和最优运行方案还会得到约束对应的对偶乘子拉格朗日乘子。这些对偶乘子揭示了如果规划方案 y 稍微变动运行成本会怎么变——这不就是梯度信息吗Benders割就是把当前运行成本和这个梯度信息组合起来形成一条关于运行成本的切平面α ≥ 当前运行成本 梯度系数 · (y - y*)主问题里 α 是运行成本的代理变量。把这条切平面加到主问题主问题就能在下一轮决策时预判运行成本的变化方向从而避免给出过差的投资方案。迭代多次这些切平面会逐渐逼近真实的运行成本函数主问题的最优解也就收敛到原问题的最优解。如果子问题不可行说明当前容量方案根本满足不了运行约束那就不能生成最优性割而要解一个最小化约束违背量的可行性子问题生成可行性割把这些问题区域直接切掉。3.2 主问题的具体形式主问题在第 k 次迭代时的形式是min ∑ CRF_i · C_i · inv_cost_i α约束包括设备选型和容量上下限约束以及已经累积的所有Benders割最优性割α ≥ obj_sub^(j) coef^(j) · (C - C^(j))j 1, ..., k可行性割feas_coef^(j) · (C - C^(j)) ≤ -viol_const^(j)j 1, ..., k主问题是小规模MILP决策变量只有选型0-1变量、容量变量和一个标量 αGurobi几乎秒解。主问题求解得到的目标函数值是原问题的下界LB因为 α 是运行成本的松弛估计而投资成本固定容量下真实运行成本给出上界UB。上下界之间的gap就是当前迭代的收敛指标。3.3 子问题与对偶提取的关键点子问题是整个程序里最需要小心的部分。固定容量 C^(k) 后每个场景单独求解运行优化。YALMIP在求解完成后可以用dual(约束句柄)直接提取对偶值非常方便但有几个坑要注意一是必须在建模时保存容量相关约束的句柄。容量参数 C 通过出力不超过容量这类不等式出现在子问题中只有这些约束的对偶值才是割平面系数需要的。能量平衡约束、储能约束的对偶虽然也有意义但它们不直接对容量求导不能混进割系数里。二是如果子问题是MILP含设备启停0-1变量取出的对偶信息可能不是有效的割系数。整数规划的对偶不是LP那种标准对偶直接用会出问题。我的做法是把设备启停这类离散决策也放到主问题的离散层或者固定整数变量后求解LP松弛再取对偶。这样割更干净收敛也快。三是符号问题。不同的人把约束写成p ≤ C还是C - p ≥ 0YALMIP给出的对偶符号可能不同。判断割方向对不对有个简单有效的现场验证方法把某台设备容量微调1%重新求解子问题看运行成本的实际变化量是否和割系数预测的方向一致。如果方向反了代码里对偶取错了要马上修正。3.4 收敛判据和整体流程整个算法迭代框架如下初始化k0LB-∞UB∞割集合为空。求解主问题得到容量方案 C^(k)主问题目标值就是当前的LB。固定 C^(k)逐个场景求解子问题累加全年运行成本加上投资成本得到本轮UB并取历史最小值。若 (UB - LB)/UB ε一般取1%或0.5%停止。根据子问题对偶生成最优性割或可行性割加入割集合kk1回到第2步。这个流程在MATLAB里就是一个while循环逻辑非常简单。但实际调试中割的质量、数值规模、初始方案这些细节会显著影响收敛轮数这些坑我在第5节专门讲。4. MATLAB程序实现与核心代码4.1 程序文件怎么组织整个项目我按功能拆成几个文件保持数据-主问题-子问题-割生成-主循环的关注点分离main_GBD.m主程序负责数据读取、参数初始化和迭代循环。data_ies.m定义系统参数、候选设备参数、负荷曲线、能源价格。build_master.m构建并求解主问题输入割集合输出容量方案和下界。solve_subproblem.m构建并求解单个场景的运行子问题返回目标值和对偶信息。gen_cut.m根据子问题对偶信息生成最优性割或可行性割。post_process.m结果整理与绘图输出容量配置、成本构成和逐时运行曲线。这样组织最大的好处是以后加储能模型、加碳约束、改目标函数只需要改对应的函数迭代框架完全不用动。我在好几个课题里复用这套框架只改数据和约束函数非常省事。4.2 数据准备与量纲处理数据这块最容易踩坑先说两个最重要的经验。第一单位要统一量级要控制在舒服的范围。我的习惯是成本统一用万元功率统一用kW时间用小时。如果你直接用元做单位投资成本可能上亿运行成本几十万量级差了三个数量级求解器数值稳定性会很差。程序里我统一把投资成本换算成万元、运行成本也以万元计主问题和子问题的目标函数都在10^2到10^3量级数值上很舒服。第二调试阶段不要一上来就用完整8760小时数据。先用3到5个典型日跑通整个GBD迭代流程确认收敛趋势正常再切换到完整场景集。我调试期靠这个习惯节省了一半时间。4.3 主问题的YALMIP实现主模型的核心代码如下思路完整可以直接作为模板function [C_opt, x_opt, alpha_opt, LB] build_master(Cuts, params) n_equip params.n_equip; % 主问题变量 x binvar(n_equip, 1); % 选型0-1变量 C sdpvar(n_equip, 1); % 容量变量 alpha sdpvar(1, 1); % 运行成本代理变量 % 容量约束不选则容量为0选型则容量在上下限之间 Constraints [C 0, C params.C_max .* x, C params.C_min .* x]; Constraints [Constraints, alpha 0]; % 运行成本非负改善初始下界 % 投资年化成本 inv_cost_annual sum(params.CRF .* params.unit_inv_cost .* C); % 累加所有Benders割 for j 1:length(Cuts) if strcmp(Cuts(j).type, optimality) Constraints [Constraints, alpha Cuts(j).obj ... Cuts(j).coef * (C - Cuts(j).C)]; else Constraints [Constraints, Cuts(j).coef * (C - Cuts(j).C) ... -Cuts(j).viol]; end end Objective inv_cost_annual alpha; ops sdpsettings(solver, gurobi, verbose, 0, mipgap, 1e-4); sol optimize(Constraints, Objective, ops); C_opt value(C); x_opt value(x); alpha_opt value(alpha); LB value(Objective); end几个细节说明一下。alpha ≥ 0 这条约束很有用第一轮主问题没有割的时候α 会被压到负值导致LB严重失真加了这条约束后初始下界合理很多。容量约束里的 C 综合能源系统规划：广义Benders分解MATLAB/YALMIP求解实践
做综合能源系统规划项目的朋友应该都有同感模型越来越大设备从冷热电三联供到储能、光伏、地源热泵越加越多约束条件层层嵌套最后维度冲到几万个变量里面还混着大量0-1整数变量。这时候直接丢给求解器整体硬解不是内存爆掉就是计算时间离谱。我这两年处理园区级综合能源系统优化规划问题时用得最顺手的方案就是广义Benders分解法GBD在MATLAB里用YALMIP建模把投资决策和运行优化拆成两层交替迭代求解程序稳定、收敛可控、改参数也方便。这篇文章就把这个项目的完整思路、数学模型、MATLAB代码实现和我在调试过程中踩过的坑一次性写清楚适合正在做能源系统规划、微电网容量配置、能源站优化设计的硕博生和工程师参考。整套程序的核心思想其实很朴素先把建哪些设备、建多大容量这个投资决策和设备建好之后全年怎么运行最省这个运行决策拆开分别用两个规模小得多的问题去求解再用割平面把两层决策之间的耦合信息传回去反复迭代直到收敛。下面我从问题本质、数学模型、算法原理、代码实现到调优经验一层一层展开。1. 先把这个项目要解决的问题说清楚1.1 综合能源系统规划的本质是什么综合能源系统Integrated Energy System, IES简单地讲就是把电、气、热、冷多种能源在源侧、网侧、荷侧打通通过燃气轮机、燃气锅炉、电制冷机、吸收式制冷机、热泵这类能量转换设备以及电储能、热储能这类缓冲设备实现多能互补与梯级利用。规划问题要回答的核心问题是在给定负荷需求曲线、能源价格和候选设备参数的条件下每类设备到底建不建、建多大容量同时按这个方案建成之后一年下来系统怎么运行经济性最好。这看起来是两个问题实际上是一个问题设备容量选大了投资成本高但运行灵活运行费用能省下来容量选小了投资省了但运行阶段可能被迫高价购电购气甚至出现供能缺口。所以投资决策和运行决策必须放在同一个框架里联合优化这就是规划-运行耦合的由来。1.2 为什么整体建模直接硬解行不通如果当前面说的所有变量全部建模到一起得到的是一个大规模混合整数规划MILP问题考虑非线性效率曲线的话就是混合整数非线性规划MINLP。做一个粗略的规模估算假设一年8760小时、5类候选设备、4种能源载体运行变量的数量级就是8760乘以设备数和能源种类随便一建模就是几万到几十万个变量再加上设备选型的0-1变量和几百个储能状态约束整体模型的复杂度立刻就不一样了。即便用Gurobi或者CPLEX这种商业求解器直接求解这种大规模MILP的耗时可能从几十分钟到几十个小时不等更难受的是很多情况下根本没办法在可接受的时间内找到可行解。我自己测试过一个中等规模园区算例整体建模直接求解内存峰值接近8GB跑了三个小时才收敛到2%的MIP gap基本没法用于参数灵敏度分析。这个问题的结构特点是投资决策变量数量很少撑死几十个但运行决策变量数量巨大成千上万。这种结构天生适合分解——先把少的投资变量定下来再在给定容量下求解运行为主的优化问题运行优化的结果反过来告诉投资层这个方案哪里不好、往哪个方向调能省钱反复迭代逼近全局最优。说白了这就是Benders分解的思想。1.3 广义Benders和经典Benders到底差在哪经典Benders分解最早是解决混合整数线性规划的它要求子问题必须是线性规划这样才能从对偶解中提取割平面。广义Benders分解Generalized Benders Decomposition, GBD是Geoffrion在1972年做的推广——子问题可以是凸非线性规划只要能从子问题的拉格朗日对偶信息中构造出割平面就行。综合能源系统规划里正好用得上这个推广。如果模型里考虑了设备效率随负荷率变化的非线性曲线、储能损耗的非线性项模型就是MINLP此时GBD几乎是唯一既理论严谨又工程可用的分解框架。当然如果像我程序默认版本这样把非线性环节用分段线性近似处理掉模型退化成MILP那用到的其实就是经典Benders的特例但整个代码框架完全复用只是子问题求解器从NLP求解器换成LP/MILP求解器而已。这也是我为什么坚持用GBD而不是其他分解方法比如拉格朗日松弛的原因GBD的收敛性是理论上有保障的而且实现起来就是主问题-子问题交替迭代对工程师非常友好。2. 数学建模怎么把规划问题写成主问题子问题的结构2.1 变量分层与能量集线器建模要在代码里实现分解第一步是建模时就把变量天然分成两层这非常关键。我习惯用能量集线器Energy Hub模型来抽象综合能源系统把系统看成一个多输入多输出的转换节点输入侧是外购电力和天然气中间经过CHP机组、燃气锅炉、电制冷机、热泵等设备转换输出侧是电、热、冷三种负荷储能设备挂在集线器内部作为缓冲。在这个框架下变量分成两块规划变量主问题变量每个候选设备的选型0-1变量 x_i 和连续容量变量 C_i二者通过 C_i ≤ x_i · C_i_max 关联。运行变量子问题变量每个典型时段下各设备的出力 p_{i,t}、储能充放功率和荷电状态 SOC_{s,t}、与外网的交换功率等。一句话总结主问题决定买哪些设备、买多大子问题回答在给定设备容量下全年怎么运行最省钱。建模的时候就要刻意保持这个分层结构不要把投资成本直接和运行成本揉进同一个表达式否则后续分解很别扭。2.2 目标函数投资年化与运行成本怎么算规划问题的目标函数通常取年等值总成本最小包含两部分第一部分是投资年化成本。设备投资不是一次性算进目标就行要按设备寿命和折现率换算成年等值成本用到的系数叫资本回收因子CRFCRF r(1r)^n / ((1r)^n - 1)其中 r 是折现率n 是设备寿命。这个细节很多初学者会忽略直接用总建设成本放进目标函数结果设备寿命不同导致对比失真优化结果明显偏离工程实际。我程序里每个设备都有独立的寿命参数和对应的CRF这是必须的。第二部分是年运行维护成本包括外购电费、天然气费以及设备单位出力的运维费用。如果课题涉及碳约束或者可再生能源配额还会额外加碳排放成本和惩罚项。目标函数总体写出来就是min ∑ (CRF_i · C_i · inv_cost_i) ∑_t ∑_s w_s · (购电成本 购气成本 运维成本)注意这里的运行成本是分场景加权求和后面细说。2.3 约束体系容量、运行、储能三类约束约束条件按层次可以分成四类建模时建议严格分开写方便后面取对偶和排错。第一类设备容量上下限约束每台设备容量要么为0不选型要么在最小可建容量和最大可建容量之间用 C_min · x ≤ C ≤ C_max · x 表达。加最小可建容量很重要否则优化器可能会给你解出一个0.01kW这种毫无工程意义的蚊子腿容量。第二类设备出力与容量的耦合约束任何时段的出力不能超过已配置容量。比如CHP机组电出力 0 ≤ p_chp,t ≤ C_chp热出力通过热电比和电出力耦合。注意这类约束就是后面Benders割平面系数的主要来源所以建模时必须把这些约束的对偶信息单独保留。第三类能量平衡约束这是子问题的核心电平衡外购电 CHP电出力 光伏出力 电池放电 电负荷 电制冷机耗电 电池充电热平衡CHP热出力 燃气锅炉热出力 蓄热罐放热 热负荷 蓄热罐蓄热冷平衡电制冷机制冷 吸收式制冷机制冷 冷负荷气平衡外购天然气 CHP耗气 燃气锅炉耗气第四类储能约束充放功率上下限、SOC状态转移方程、调度周期始末SOC相等约束。储能约束里SOC的初末相等是规划类问题里容易被漏掉的一条不加上会导致储能被白嫖优化结果明显失真。2.4 典型场景日怎么选直接拿8760小时全塞进约束里子问题规模也扛不住。工程上通用的做法是场景缩减用k-means聚类把全年负荷、光伏出力、能源价格曲线聚成若干个典型日比如我常用的组合是4个季节乘以工作日/节假日再乘以典型日类型取12个场景每个场景带一个权重 w_s表示该类型在全年的天数占比。做了场景缩减之后子问题可以进一步拆成每个场景独立的运行优化问题。这个拆分非常重要——它直接让子问题变成多个小规模LP并行求解每一轮迭代几秒钟就能跑完。而且每个场景独立求解后各自生成各自的割平面最后汇总到主问题逻辑上非常干净。3. 广义Benders分解算法迭代逻辑与割平面原理3.1 割平面是怎么来的很多朋友第一次接触Benders分解最迷惑的就是割平面那串公式。我用直白的方式解释一下。把原问题抽象成 min f(x,y)其中 x 是运行变量y 是规划变量。固定 y y* 之后子问题是关于 x 的凸优化问题。求解子问题除了得到最优运行成本和最优运行方案还会得到约束对应的对偶乘子拉格朗日乘子。这些对偶乘子揭示了如果规划方案 y 稍微变动运行成本会怎么变——这不就是梯度信息吗Benders割就是把当前运行成本和这个梯度信息组合起来形成一条关于运行成本的切平面α ≥ 当前运行成本 梯度系数 · (y - y*)主问题里 α 是运行成本的代理变量。把这条切平面加到主问题主问题就能在下一轮决策时预判运行成本的变化方向从而避免给出过差的投资方案。迭代多次这些切平面会逐渐逼近真实的运行成本函数主问题的最优解也就收敛到原问题的最优解。如果子问题不可行说明当前容量方案根本满足不了运行约束那就不能生成最优性割而要解一个最小化约束违背量的可行性子问题生成可行性割把这些问题区域直接切掉。3.2 主问题的具体形式主问题在第 k 次迭代时的形式是min ∑ CRF_i · C_i · inv_cost_i α约束包括设备选型和容量上下限约束以及已经累积的所有Benders割最优性割α ≥ obj_sub^(j) coef^(j) · (C - C^(j))j 1, ..., k可行性割feas_coef^(j) · (C - C^(j)) ≤ -viol_const^(j)j 1, ..., k主问题是小规模MILP决策变量只有选型0-1变量、容量变量和一个标量 αGurobi几乎秒解。主问题求解得到的目标函数值是原问题的下界LB因为 α 是运行成本的松弛估计而投资成本固定容量下真实运行成本给出上界UB。上下界之间的gap就是当前迭代的收敛指标。3.3 子问题与对偶提取的关键点子问题是整个程序里最需要小心的部分。固定容量 C^(k) 后每个场景单独求解运行优化。YALMIP在求解完成后可以用dual(约束句柄)直接提取对偶值非常方便但有几个坑要注意一是必须在建模时保存容量相关约束的句柄。容量参数 C 通过出力不超过容量这类不等式出现在子问题中只有这些约束的对偶值才是割平面系数需要的。能量平衡约束、储能约束的对偶虽然也有意义但它们不直接对容量求导不能混进割系数里。二是如果子问题是MILP含设备启停0-1变量取出的对偶信息可能不是有效的割系数。整数规划的对偶不是LP那种标准对偶直接用会出问题。我的做法是把设备启停这类离散决策也放到主问题的离散层或者固定整数变量后求解LP松弛再取对偶。这样割更干净收敛也快。三是符号问题。不同的人把约束写成p ≤ C还是C - p ≥ 0YALMIP给出的对偶符号可能不同。判断割方向对不对有个简单有效的现场验证方法把某台设备容量微调1%重新求解子问题看运行成本的实际变化量是否和割系数预测的方向一致。如果方向反了代码里对偶取错了要马上修正。3.4 收敛判据和整体流程整个算法迭代框架如下初始化k0LB-∞UB∞割集合为空。求解主问题得到容量方案 C^(k)主问题目标值就是当前的LB。固定 C^(k)逐个场景求解子问题累加全年运行成本加上投资成本得到本轮UB并取历史最小值。若 (UB - LB)/UB ε一般取1%或0.5%停止。根据子问题对偶生成最优性割或可行性割加入割集合kk1回到第2步。这个流程在MATLAB里就是一个while循环逻辑非常简单。但实际调试中割的质量、数值规模、初始方案这些细节会显著影响收敛轮数这些坑我在第5节专门讲。4. MATLAB程序实现与核心代码4.1 程序文件怎么组织整个项目我按功能拆成几个文件保持数据-主问题-子问题-割生成-主循环的关注点分离main_GBD.m主程序负责数据读取、参数初始化和迭代循环。data_ies.m定义系统参数、候选设备参数、负荷曲线、能源价格。build_master.m构建并求解主问题输入割集合输出容量方案和下界。solve_subproblem.m构建并求解单个场景的运行子问题返回目标值和对偶信息。gen_cut.m根据子问题对偶信息生成最优性割或可行性割。post_process.m结果整理与绘图输出容量配置、成本构成和逐时运行曲线。这样组织最大的好处是以后加储能模型、加碳约束、改目标函数只需要改对应的函数迭代框架完全不用动。我在好几个课题里复用这套框架只改数据和约束函数非常省事。4.2 数据准备与量纲处理数据这块最容易踩坑先说两个最重要的经验。第一单位要统一量级要控制在舒服的范围。我的习惯是成本统一用万元功率统一用kW时间用小时。如果你直接用元做单位投资成本可能上亿运行成本几十万量级差了三个数量级求解器数值稳定性会很差。程序里我统一把投资成本换算成万元、运行成本也以万元计主问题和子问题的目标函数都在10^2到10^3量级数值上很舒服。第二调试阶段不要一上来就用完整8760小时数据。先用3到5个典型日跑通整个GBD迭代流程确认收敛趋势正常再切换到完整场景集。我调试期靠这个习惯节省了一半时间。4.3 主问题的YALMIP实现主模型的核心代码如下思路完整可以直接作为模板function [C_opt, x_opt, alpha_opt, LB] build_master(Cuts, params) n_equip params.n_equip; % 主问题变量 x binvar(n_equip, 1); % 选型0-1变量 C sdpvar(n_equip, 1); % 容量变量 alpha sdpvar(1, 1); % 运行成本代理变量 % 容量约束不选则容量为0选型则容量在上下限之间 Constraints [C 0, C params.C_max .* x, C params.C_min .* x]; Constraints [Constraints, alpha 0]; % 运行成本非负改善初始下界 % 投资年化成本 inv_cost_annual sum(params.CRF .* params.unit_inv_cost .* C); % 累加所有Benders割 for j 1:length(Cuts) if strcmp(Cuts(j).type, optimality) Constraints [Constraints, alpha Cuts(j).obj ... Cuts(j).coef * (C - Cuts(j).C)]; else Constraints [Constraints, Cuts(j).coef * (C - Cuts(j).C) ... -Cuts(j).viol]; end end Objective inv_cost_annual alpha; ops sdpsettings(solver, gurobi, verbose, 0, mipgap, 1e-4); sol optimize(Constraints, Objective, ops); C_opt value(C); x_opt value(x); alpha_opt value(alpha); LB value(Objective); end几个细节说明一下。alpha ≥ 0 这条约束很有用第一轮主问题没有割的时候α 会被压到负值导致LB严重失真加了这条约束后初始下界合理很多。容量约束里的 C