做综合能源系统规划的时候,我之前一直习惯把“投资决策”和“运行调度”揉成一个混合整数线性规划直接丢给求解器。直到接手某高校课题组的一个园区级电-热-冷联供系统规划项目,模型里同时出现设备选型、容量配置和多典型日运行耦合调度,问题规模一下子顶到了几万个变量,商业求解器直解要跑几十个小时,还不一定收敛到像样的整数解。后来换成广义Benders分解法在Matlab里重写求解框架,收敛速度快了一个数量级。这篇博客想把整个项目的设计思路、算法原理、Matlab实现细节和踩坑记录完整拆出来,给正在做综合能源系统优化规划、或者想用广义Benders分解解决混合整数大规模规划问题的朋友一个可以直接参考的落地模板。
1. 项目整体设计:综合能源系统规划为什么选了广义Benders
1.1 这个项目到底在解决什么问题
综合能源系统优化规划,核心是回答两个问题:第一,园区、建筑群或区域级的冷热电多种能源系统中,风、光、气、储、冷热转换等设备要建哪些、建多大容量;第二,在典型运行场景下,这些设备怎么配合电网购电、燃气购气来满足负荷需求,使得投资加运行的整体经济性最优。设备选型和容量属于“一次性花钱”的投资决策,用整数变量或者离散档位变量描述;而运行调度是给定设备方案后,每一时刻各设备出力、储能充放的连续决策问题。
这两个层次天然耦合:容量定大了,投资浪费;容量定小了,运行期高电价时段还要大量购电,运行成本居高不下。所以规划模型必须把投资成本和运行成本放在同一个目标函数里联合优化,也就是业内常说的“规划-运行一体化”。这类问题的典型数学形态是一个混合整数线性规划(MILP)或混合整数非线性规划(MINLP),其中整数变量数量取决于候选设备数量和容量档位数,连续变量数量由运行时段数、场景数、设备数量和网络节点数共同决定。
1.2 全模型直接求解的三大瓶颈
一开始我的做法是把所有变量和约束一次性写入求解器,从工程直觉上这最简单。但项目推进到中型规模时,发现三个非常实际的瓶颈。
第一是整数搜索空间爆炸。园区级候选设备可能多达十几类,每类又有多个可选容量档位,整数变量组合数呈指数增长。商业求解器的分支定界法在数万个连续变量的背景下搜索整数解,分支树会深得可怕。
第二是连续变量之间的强耦合。储能设备的荷电状态(SOC)在时间维度上递归耦合,CHP机组的热电出力又有可行域耦合,再加上多能平衡,系数矩阵的稀疏性虽然不错,但在求解器内部预求解和割平面生成过程中仍然非常耗时。
第三是调试体验极差。模型稍微加一个约束或调整一个参数,整套求解从头开始,分支树的搜索模式也变了。想判断某个设备方案好不好、某个运行约束是否过紧,完全无从下手。这种“黑箱”式调试在项目沟通过程中效率很低,因为对方要看的不仅是结论,还有中间方案的经济性敏感性。
1.3 广义Benders分解在这类问题中的定位
广义Benders分解(Generalized Benders Decomposition,GBD)的思路,是把一个难解的大问题拆成“两层的谈判”:一个只管投资决策的主问题,一个在给定投资方案下做运行优化的子问题。主问题是一棵相对小很多的整数搜索树,子问题则是一个纯连续变量的凸优化问题,可以用LP或QP求解器快速求到全局最优并拿到对偶乘子。
之所以选择广义Benders而不是传统的Benders分解,是因为传统Benders主要针对子问题为线性规划的情形,而综合能源系统中常见CHP机组的可行域、储能损耗、非线性热功率方程等会让子问题变成凸非线性规划。广义Benders把子问题的要求放宽到“凸”即可,适用范围明显更广。从项目实操角度讲,即便我的子问题最终通过线性化成为了LP,用广义Benders框架去实现也没有损失,反而保留了直接接入非线性凸子问题的能力。
这个分解策略还天然契合规划问题的规模规律:主问题的整数变量只占总量的一小部分,子问题虽然连续变量多,但没有整数分支,可以被高效求解。于是整体计算负担被均匀拆分,每一轮迭代只需要解一个小规模MILP和几个连续的LP/QP,速度对比直解有本质提升。
2. 广义Benders分解的算法机制与适用边界
2.1 从经典Benders到广义Benders只差一个“凸”字
Benders分解最早是从混合整数线性规划中提炼出来的,核心思想是把整数变量固定下来,把原来复杂的MILP变成一个关于连续变量的LP子问题,再通过对偶理论把子问题的最优值表达成关于整数变量的函数,用割平面逐次逼近。广义Benders是它的自然推广,关键在于子问题从“线性”放宽到“凸”。广泛意义上,只要固定复杂变量后,问题关于其余变量是凸优化,就可以利用凸优化的极值条件和对偶理论构造割平面。
这个放宽在实际中的意义非常大。综合能源系统的运行层经常出现二次能耗、效率曲线、水头损失等凸函数,如果强行线性化会带来精度损失或增加大量辅助变量。在广义Benders框架下,子问题可以直接保留凸二次目标甚至部分凸非线性约束(在求解器支持的前提下),不需要为了套用算法而牺牲模型精度。
2.2 主问题、子问题各自管什么
在一个典型的两阶段规划框架中,第一阶段变量是设备是否建设x和建设容量c,这是主问题的专属变量。第二阶段变量包括各典型日每个时段的设备出力、储能充放功率、外购电/气量等,属于子问题变量。
主问题的结构是:
- 目标由投资成本和一个额外标量变量
theta组成,theta表示对运行成本的逐步逼近; - 约束包括投资逻辑约束(先建设后才能有容量、容量上下限、设备互斥、场地总面积限制等);
- 随着迭代不断加入Benders割,每一条割都要求在某一组已考察的投资方案附近,运行成本不能低于该方案对应的线性下界近似。
子问题的结构则是:
- 输入是主问题决策出来的一组
(x, c); - 需要求解的是固定设备方案后的多场景、多时段运行调度问题;
- 输出包括该方案下的最小运行成本,以及所有涉及投资变量的约束对偶乘子;
- 如果子问题不可行,则输出一组可行性修复信息,用于生成可行割。
整体逻辑就像甲方乙方反复沟通:甲方先表态“我想建这些设备”,乙方按照这个方案算一个“最低运行开销”,同时告诉甲方“你这个方案如果想再省一点,应该在哪个方向调整投资”。甲方收到反馈后更新方案,继续迭代,直到双方对总成本的估计达成一致。
2.3 割平面是怎么算出来的
割平面是整个算法的技术核心。假设主问题给出了某轮投资方案(x_0, c_0),子问题求解得到最小运行成本R(x_0, c_0),并且我们能从子问题的求解器中取到关于投资相关等式/不等式约束的对偶乘子,记作lambda。那么经典的最优割形式为:
theta >= R(x_0, c_0) + lambda^T * (投资变量 - 投资变量_0)
这里的关键直觉是:因为子问题是凸的,运行成本作为投资变量的函数是凸函数,凸函数在某个点处的切线(一次近似)永远位于函数下方。所以把这个切线作为约束加回主问题后,theta的取值一定会小于等于真实运行成本,不会把主问题带偏;随着割平面越加越多,主问题中theta的下界估计会越来越接近真实函数。
如果子问题在当前投资方案下无解,比如储能容量太小导致某些时刻无法满足负荷平衡,那就不能直接使用最优割。此时需要求解一个带松弛变量的可行性子问题,最小化不可行程度,同样能拿到一组对偶乘子,生成一个“可行割”:它告诉主问题,现有投资方案已经落入了不可行区域,必须要么提升设备容量,要么调整设备组合。可行割的意义在于把主问题的搜索空间“削掉”已经验证不可行的区域。
2.4 收敛判据与迭代终止条件
广义Benders分解的迭代过程,始终维护两个界:下界来自主问题当前最优目标值(投资成本加theta的下界逼近),上界来自每一轮可行子问题得到的真实总成本(投资加该方案的真实运行成本,通常取历史最小者)。
两者的差值除以当前上界,就是迭代间隙:
gap = (UB - LB) / UB
当gap小于预设阈值,例如0.5%或1%时,就可以认为找到了经济上最优的规划方案。实际项目中我一般再叠加一个最大迭代次数作为保险,比如200轮,防止数值振荡导致gap迟迟不收敛。停机后取当前最优上界对应的投资方案作为推荐方案,因为上界对应的是“真实可运行且已计算过实际成本”的方案,可信度更高。
2.5 什么情况下不要硬上GBD
虽然GBD方法论很漂亮,但并非所有综合能源规划问题都适合。最大前提是:固定投资变量后,子问题必须确实是凸优化。如果运行模型里有严格的非凸约束,比如离散档位运行状态、设备启停的固定成本和最小启停时间等,仅仅把整数变量放进主问题还不够,因为子问题里如果还有0-1变量,那么对偶乘子可能取不到,割平面的下界性质也不再保证。
另外,如果整数变量和连续变量的耦合非常紧密,每轮主问题和子问题几乎没有信息交换的空间,那么GBD的迭代次数可能很大,反而不如直接使用大规模MILP求解器配合强割平面策略。判断方式很简单:先把整数变量固定后看剩余问题是不是“容易求”的凸问题,如果是且整数候选量不超过几千,GBD通常值得尝试;如果不是,就要考虑其他分解框架,比如拉格朗日松弛或者列约束生成。
3. Matlab代码实现:核心模块拆解与关键代码
3.1 项目文件结构规划
在动手写Matlab代码前,我强烈建议先把文件结构规划清楚。一个可复用、可调试的GBD工程,至少应该包含以下模块:
main.m:主程序入口,负责参数加载、初始化、迭代控制和结果汇总;data_ies.m:所有输入数据的定义,包括负荷曲线、设备候选参数、价格参数、典型日权重等;build_master_problem.m:构建主问题模型,返回变量句柄和约束句柄;solve_subproblem.m:输入特定投资方案,构建并求解子问题,返回运行成本、对偶乘子、可行状态;add_benders_cut.m:根据子问题返回信息生成并往主问题中添加新割;plot_results.m:收敛曲线、设备容量柱状图、运行出力图表的绘制。
这样拆分的好处是,后续做参数敏感性分析或场景扩展时,只需要改动数据模块和子问题模块,不需要动主循环结构。
3.2 数据与基础参数设计
数据设计直接决定模型的物理意义和数值稳定性。我先定义候选设备集合,常见设备类型包括:风力发电、光伏发电、燃气轮机、燃气锅炉、电锅炉、电制冷机、吸收式制冷机、蓄电装置和蓄热装置。每一类设备需要提供单位投资成本、设计寿命、维护成本、效率/COP、出力上下限、爬坡限制或储能充放倍率等参数。
典型日设计是另一个容易忽视但影响巨大的环节。综合能源规划一般不逐日模拟全年8760小时,而是抽取多个典型日,每个典型日按权重缩放至全年代表天数。比如冬季典型日、夏季典型日、过渡季典型日,每个典型日的24小时电负荷、热负荷、冷负荷、光照强度、风力预测数据都要准备好。我习惯把典型日权重放在一个向量里,例如dayWeight = [120, 120, 125],代表三类典型日在一年中占多少天。
所有单位和量纲必须统一。我在项目里统一采用千瓦(kW)作为功率单位,千元(kCNY)作为货币单位。这句话看起来基础,但在实际工程模型中数据单位不对称导致求解器收敛异常,是出现频率最高的低级错误之一。
3.3 主问题建模与割平面的动态添加
主问题使用Yalmip建模,求解器可以选用CPLEX或Gurobi这类商业MIP求解器。下面是一段在Matlab中构建主问题核心框架的代码。
% 主问题:投资决策 + Benders 割(theta逼近运行成本) x = binvar(Ndev, 1); % 设备是否建设 cap = sdpvar(Ndev, 1); % 建设容量 theta = sdpvar(1, 1); % 运行成本逼近变量 % 投资成本:固定投资部分 + 容量相关投资部分(折算为等年值) C_invest = x' * fixCost + cap' * varCost; Constraints = []; % 投资逻辑约束:建设后容量必须在可行范围内 for i = 1:Ndev Constraints = [Constraints, cap(i) <= capMax(i) * x(i)]; Constraints = [Constraints, cap(i) >= capMin(i) * x(i)]; end % 逐轮添加的Benders割,Alpha和Beta是历史割平面的系数 [nCuts, ~] = size(Alpha); for k = 1:nCuts Constraints = [Constraints, theta >= Alpha(k) + Beta(k,:) * [x; cap]]; end ops = sdpsettings('solver', 'cplex', 'verbose', 2, 'mip.tolgap', 1e-4); sol = optimize(Constraints, C_invest + theta, ops); x_val = value(x); cap_val = value(cap); LB = value(C_invest + theta);主问题中的theta变量本质上是运行成本的“占位符”。初始时如果一条割都没有,theta没有约束会被直接压到负无穷,所以第一轮迭代我会先对主问题增加一个任意宽松的初始下界约束,比如theta >= -1000000,保证主问题不是无界状态。
mip.tolgap参数的设置需要谨慎。太大会导致主问题本身求解粗糙,产生无效割;太小会让主问题在每轮迭代花费大量时间。实测下来,规划类问题的主问题MIP gap设置到1e-4比较合适,既保证割质量,又不会过度消耗分支定界时间。
3.4 子问题建模与对偶信息提取
子问题的构建是运行优化的核心。输入为上一轮主问题输出的x_val和cap_val,在这个阶段x_val和cap_val已经是确定数值,因此子问题中不再有任何整数变量。下面展示单个典型日的运行子问题核心代码,实际项目里多个典型日可以通过循环或parfor并行求解。
% 子问题:给定投资方案下的单典型日运行优化 T = 24; % 每小时一个时段 p_GT = sdpvar(T, 1); % 燃气轮机出力 p_buy = sdpvar(T, 1); % 从电网购电 e_ch = sdpvar(T, 1); % 蓄电充电功率 e_dis = sdpvar(T, 1); % 蓄电放电功率 SOC = sdpvar(T+1, 1); % 荷电状态,含初始时刻 Constraints = []; % 电功率平衡 Constraints = [Constraints, p_GT + p_buy + p_PV + p_WT + e_dis ... == e_ch + e_load]; % 设备出力上下限 Constraints = [Constraints, 0 <= p_GT <= cap_val(idx_GT)]; Constraints = [Constraints, 0 <= p_buy <= buyMax]; Constraints = [Constraints, e_ch >= 0, e_dis >= 0]; % 储能约束 Constraints = [Constraints, SOC(1) == s0]; Constraints = [Constraints, SOC(2:T+1) == SOC(1:T) + ... e_ch * eta_ch - e_dis / eta_dis]; Constraints = [Constraints, 0 <= SOC <= cap_val(idx_ES) * 0.9]; % 热功率平衡等其他平衡约束 % ... C_run = sum(p_GT .* gasPrice) * gasPriceScale + ... sum(p_buy .* elecPrice) * elecPriceScale + ...; ops = sdpsettings('solver', 'linprog', 'verbose', 0); sol = optimize(Constraints, C_run, ops); if sol.problem == 0 run_cost = value(C_run); duals = dual(Constraints); % Yalmip 直接取对偶乘子 feasible = true; else run_cost = []; duals = []; feasible = false; end子问题必须是纯LP或QP才能用dual()函数直接拿到对偶乘子。如果使用fmincon求解非线性凸子问题,需要从lambda输出参数中手动提取。Yalmip的dual()在子问题是线性规划时非常方便,返回的就是约束对应的对偶乘子向量。
关键技巧是,最终用于生成割平面的乘子,只关心那些与投资变量直接相关的约束对偶值。比如储能容量约束的对偶值,反映的是“容量每增加一点,运行成本能降低多少”;购电上限约束的对偶值,反映的是“电网容量裕度对成本的边际影响”。把这些乘子整理成Beta向量,配合子问题在(x_0, c_0)处的运行成本和投资方案本身,就可以组装出新的割平面。
3.5 循环控制、gap判断与结果保存
GBD的主循环相对固定,但一些细节值得仔细处理。我的标准循环结构如下:
UB = inf; LB = -inf; gap = inf; iter = 0; Alpha = []; Beta = []; history_LB = []; history_UB = []; while gap > tol && iter < maxIter iter = iter + 1; % 第一步:求解主问题,得到新投资方案 [x_val, cap_val, LB_new, Alpha, Beta] = solve_master(Alpha, Beta); % 第二步:求解所有典型日子问题 [run_cost_total, duals_cons, feasible] = solve_all_subproblems(x_val, cap_val); % 第三步:生成并添加新割 if feasible UB = min(UB, C_invest_val + run_cost_total); [Alpha, Beta] = generate_optimal_cut(...); else [Alpha, Beta] = generate_feasibility_cut(...); end gap = abs(UB - LB_new) / abs(UB); history_LB(end+1) = LB_new; history_UB(end+1) = UB; end当UB初值为无穷大时,gap计算时要先判空。很多初学者在代码里漏掉了“第一轮子问题不可行”的边界情况,导致gap产生NaN。我会在主循环前先执行一轮“预检”,挑一组保守的较大容量方案跑一次子问题确保可行,并把算出来的总成本作为UB初始值。这个预检方案虽然投资偏高,但能保证后续割平面从可靠上界开始迭代,收敛速度明显更快。
3.6 并行化与求解器参数调整
在典型日较多或使用集群计算时,子问题的并行化是最直接的提速手段。多个典型日之间的子问题完全独立,没有任何共享变量,因此可以直接用parfor并行求解,只需要在循环内收集每个典型日的运行成本和对偶向量,再汇总成整体运行成本即可。我在实际项目中把4个典型日子问题并行后,单轮迭代时间从大约40秒降到12秒,效果立竿见影。
求解器参数方面,主问题求解器我通常把NodeLimit或TimeLimit作为兜底,避免某轮主问题在极端情况下陷入超长搜索。子问题由于是线性规划,没有太多参数可调,主要注意关闭输出避免控制台刷屏,以及在矩阵规模大时开启求解器的预求解和缩放选项。如果遇到“Current solution is not feasible”之类的报错,优先检查子问题模型中的平衡约束与容量边界是否有矛盾。
4. 从数学模型到规划方案:结果分析与算例复盘
4.1 某园区IES规划算例复盘
为了验证代码的可靠性,我构造了一个园区级冷热电联供系统的简化算例。候选设备包括一台燃气轮机、一台燃气锅炉、一台电制冷机、一组光伏阵列、一套蓄电装置和一套蓄热装置;典型日选取冬季和夏季两个典型日,每个典型日24时段电、热、冷负荷各不相同。
直接求解完整MILP模型时,官方求解器在两个小时里给出的gap仍然卡在5%左右;而改用GBD分解后,总共迭代54轮,gap降到0.8%,累计耗时约20分钟。最终规划方案选择了燃气轮机+蓄电装置+电制冷机的组合,燃气锅炉被淘汰,原因是冬季热负荷虽然高,但燃气轮机余热回收已经覆盖了大部分需求,额外投资燃气锅炉的经济性不佳。这个结论在传统直解模型中得到了一致验证,充分说明GBD的收敛方案与全局最优方案的偏差在可接受范围之内。
4.2 收敛过程怎么看
看收敛曲线是判断算法是否健康的第一手段。我在项目里记录每一轮迭代的LB和UB,绘制成两条逐渐靠近的曲线。健康的收敛过程应该是:LB从很低的初始值快速上升,UB从预检得到的较保守上界缓慢下降,两条曲线在前20轮内迅速靠拢,之后每次下降的幅度逐渐变小。
如果LB长时间纹丝不动,而UB在波动,说明割平面提供的信息对主问题来说不够有区分度。这时候需要检查是不是运行成本函数在投资方案变化的范围内几乎平坦,比如储能设备容量增大后,运行成本只下降了一点点,导致割平面的斜率很小,主问题没有动力去探索更大的容量。解决办法是对设备容量做离散化,或者增加投资成本的惩罚系数敏感性分析。
4.3 规划方案怎么解读
规划方案不能只看总成本,还要看各设备容量的利用率和典型日运行方式是否合理。我习惯把每个典型日的设备出力堆叠成面积图,观察燃气轮机是否在电价高峰时段满发,储能是否在低谷充电、高峰放电,电制冷机和燃气轮机余热制冷是否实现了冷负荷互补。利用率太低的设备意味着投资冗余,即使它在优化解里被选中,也可能是受逻辑互斥约束或整数组合的限制才被保留,需要人工研判是否确实有必要。
Benders分解的副产品之一是方案经济性分析。每一轮割平面的系数,实际上给出了设备容量的边际价值。当最终线性化收敛后,查看与容量约束相关的对偶乘子,可以知道“如果该设备容量增加1kW,总成本能下降多少”。这个信息对业主做投资决策非常有价值,也是广义Benders方法相比直接丢黑箱求解器分割方案的一大优势。
5. 常见问题与排查经验实录
5.1 子问题不可行:可行割的生成与修复
子问题不可行是GBD实现中遇到最多的问题。原因通常是主问题在某轮探索出了一个过于激进的投资方案,例如储能容量极小但负荷尖峰极高,或者没有建设任何发热设备却要满足冬季热负荷。此时不能直接跳过子问题,而是要求解一个“最小不可行度”的辅助模型:
min sum(所有松弛变量的绝对值) s.t. 原约束每一项都加上对应松弛变量之后成立,松弛变量 >= 0
如果这个辅助模型的松弛量>0,说明当前投资方案确实不可行。从辅助模型的对偶乘子中提取割平面系数,得到一条可行割。这条割作为一个“惩罚”,会让主问题在下一次迭代中避开该投资区域。注意,可行割的系数符号与最优割不同,添加时不能混用,否则会破坏主问题约束的一致性。我在代码里用cutType标志区分两类割,分别存入不同的系数矩阵。
5.2 迭代不收敛或收敛极慢
有几类我会首先检查的原因。第一,子问题中有未被消除的离散变量,导致割平面不满足凸下界条件;第二,主问题MIP求解gap设置过大,主问题返回的“最优解”其实质量很差,生成的割偏离真实包络;第三,割平面系数数值太大或太小,造成主问题数值病态。
有一个非常容易踩的坑是把“松弛惩罚项”误留在正常运行子问题里。有些初学者为了规避不可行,直接在子问题约束中加入高额惩罚的松弛变量,这样虽然避免了不可行问题,但运行成本中混入了人工惩罚成本,Benders割虽然仍能保证一定的近似性,但上界不再对应真实的运行成本,最终方案会扭曲。正确做法是:正常运行子问题绝不允许松弛,可行性子问题才允许。
5.3 对偶乘子提取失败
Yalmip的dual(Constraints)返回的类型是数组或者空值。如果返回的是空值,通常是求解器没有在内部保留对偶信息,或者子问题的求解器不支持返回对偶。此时可以改用export命令将模型导出为线性规划标准形式,再用linprog的lambda输出获取对偶值。另外一种更隐蔽的原因是优化问题中出现了严格等于的约束(等式约束也需要对偶),Yalmip对等式约束同样会返回对偶乘子,但如果该等式约束实际上在预求解阶段被消掉了,就会导致乘子缺失。
我的处理方法是子问题构建完成后先求一次对偶验证,对每条关键约束检查dual值是否为NaN或Inf,一旦异常立即中断并输出约束名,避免带着脏乘子进入割平面。
5.4 数值病态导致求解器异常
典型表现在主问题求解时报“Numerical trouble”或“Infeasible or unbounded”,而模型本身逻辑上并没有问题。排查步骤是先检查系数矩阵的最大数量级和最小数量级差距,是否超过1e10。电网购电价格动辄几百元,而储能SOC增益可能只有0.9,两者数量级差异并不大;真正容易出问题的是把投资成本按“万元”和“元”混用,导致单位容量成本和其他系数的量纲差出几个数量级。
解决方案是对所有连续变量和成本系数做归一化。把功率统一换算为标幺值(比如基准容量1000kW),把成本统一到同一货币单位,必要时将目标函数总体除以某个常数缩小到1~100量级。做完归一化之后,主问题和子问题的求解稳定性会显著提升,我遇到过一个案例,归一化后迭代次数从八十多轮降到四十轮左右。
6. 扩展应用与代码复用建议
6.1 多典型日与随机场景扩展
当问题从确定性规划扩展到多场景随机规划时,GBD的框架几乎不需要变化,只需要把子问题按场景拆分。每一个场景的每个典型日都是一个独立的子问题,可以在每一轮迭代中并行求解,并把每条割都加入主问题,形成“多割Benders”结构。当然,为了更快收敛也可以采用单割聚合方式,把多个场景的割平面按场景概率加权合并,但多割方式的近似更紧,中小规模场景下优先推荐。
6.2 从规划到运行一体化的横向迁移
这个代码框架不仅适用于综合能源系统规划,也能迁移到配电网储能选址定容、综合能源服务商的年度购电计划、冷链园区多能源设备置换等场景。关键通用逻辑是:有“投资-运营”双层结构、整数决策量只占少部分、剩余部分为凸优化的问题,都可以尝试把这个Matlab框架直接改造成求解器。横向迁移时只需要重新定义主问题变量和子问题约束,Benders迭代骨架不需要推翻重写。
6.3 代码结构优化的个人建议
如果让我重新组织这套代码,我会在最初设计时就把子问题封装成独立的“运行成本黑箱函数”,只暴露三个接口:输入投资方案、返回运行成本、返回对偶乘子。这样后续从单典型日扩展到多个典型日、从确定性到随机场景、从线性子问题到凸二次子问题,都只需要改子问题内部,主循环和割平面生成逻辑完全不动。另外,割平面的历史数据最好每次迭代结束后写入磁盘或工作区备份,一旦出现异常中断,可以从未完成轮次继续启动,而不必从头开始。
这个内容后续还可以这样扩展:在主问题中加入碳排约束或可再生能源配额制约束,把割平面改为带碳排放影子价格的增广割,就能把经济规划升级为低碳规划;或者在子问题中加入运行备用约束,让规划方案天然考虑系统韧性。GBD的魅力正在于此:算法框架稳定,扩展空间全看你对模型的刻画需求。