简介面向电气工程、电力系统方向的毕业设计或论文复现场景这份MATLAB源程序包聚焦两阶段鲁棒优化在含电热联合调度的配电网中的应用用于处理分布式电源与负荷不确定性。压缩包仅321KB共12个文件主体为10个.m源程序另有1张系统示意图和1份docx模型说明结构清晰便于按模块阅读。程序完整覆盖拉丁超立方采样、K-means场景聚类、手肘法聚类数确定、1-范数与∞-范数概率置信区间构建以及两阶段分解与列与约束生成算法求解同时结合盒式分解消除储能与机组组合的时间耦合将第二阶段转化为单时间尺度问题可复现动态重构与无功优化协同的综合调度模型。已有168人浏览学习适合电力系统方向研究生和毕设学生用于算法对比、论文复现及进一步扩展研究。1. 两阶段鲁棒优化的电热范数程序先想清楚它到底在优化什么如果你是从知网摸到《两阶段鲁棒优化-电热-范数程序》这篇论文又顺手拿到了对应的94号源程序第一件事别急着跑代码先把“两阶段鲁棒优化”这个词拆干净。它解决的是这么一类问题电负荷、热负荷和风电出力都是波动的你必须在波动发生之前把机组启停、购电合同这类决策定下来等实际值出来以后再调整火电、热电联产机组、电锅炉和蓄热罐出力。前者叫第一阶段决策后者叫第二阶段重调度两个阶段合起来的目标是让最坏情况下的总成本也可接受。电热系统因为热网、蓄热罐的存在天然给第二阶段留了调整余地。而标题里的“范数”不是修辞它决定你把“最坏情况”圈在多大一个边界里。这篇笔记就按“建模、求解、避坑、验证”的顺序讲清楚这类程序怎么写、怎么调、怎么判断结果没写错。2. 电热两阶段鲁棒优化建模用向量范数约束不确定集拆出主问题与子问题2.1 为什么电热系统特别适合两阶段鲁棒热网慢动态提供了“后悔药”电力平衡要求瞬间完成风电一波动火电和储能就得立刻跟上。热力不一样。热水在管道里流动有延迟建筑物墙体本身在蓄热蓄热罐可以把多余的热存起来热负荷的波动不会像电负荷那样直接冲击机组。这意味着第二阶段做出调整时可行域比纯电力系统大得多鲁棒优化做出来的解不会过于保守。常见做法是把电、热两套平衡方程同时放进第二阶段约束风电出力偏差进入电平衡热负荷偏差进入热平衡蓄热罐的充放热功率作为第二阶段变量参与调节。这样一来热力系统的惯性被显式建模成“灵活性资源”而不是被忽略掉。两阶段鲁棒优化在综合能源系统里流行的另一个原因是参数少。随机优化要生成上百个场景才能覆盖风电的分布特性两阶段鲁棒只需要一个不确定集和几个范数参数就能刻画最坏情况。对工程确认来说这种“少参数、保底结论”的思路更容易向调度人员解释。2.2 用向量范数把不确定集圈起来从∞-范数到1-范数预算约束不确定参数一般写成预测值加偏差的形式w_k w_k^0 δ_k·Δw_k其中 δ_k 在 [-1,1] 之间取值Δw_k 是最大偏差幅度。把所有 δ_k 拼成一个向量 δ 后需要限制它的“总扰动预算”不然每个参数都取到最坏值解就会保守到没法用。典型做法是用向量的范数来定义不确定集W { δ | ‖δ‖∞ ≤ 1, ‖δ‖1 ≤ Γ }这个集合的意思是每个参数单独最多偏差一个单位但所有参数的偏差绝对值之和不能超过 Γ。Γ 的取值从 0 到 KK 是不确定参数的个数。Γ0 退化为确定性模型ΓK 等价于盒式不确定集是最保守的情况。实际工程里 Γ 取 K 的一半左右就能兼顾鲁棒性和经济性。与之对比还有纯盒式集合只用 ∞-范数和椭球集合用 2-范数。盒式集合建模简单但最坏场景所有参数同时发生过于极端椭球集合更接近统计置信域但会引入二次约束线性求解器就不能直接用了。所以“∞-范数 1-范数预算”的组合最流行它本质是一个多面体可以保持线性规划结构。范数的阶数选择直接决定源程序里约束怎么写、求解器选哪个。2.3 两阶段结构拆解第一阶段定合同第二阶段做电热再调度两阶段鲁棒优化的通用形式可以写成min_x { c^T x max_{δ∈W} min_{y∈F(x,δ)} d^T y }外层 x 是第一阶段决策比如日前签订的购电电量、机组的最小技术出力、启停状态。内层对每个 δ 求第二阶段最小成本再对外层取最坏情况。这里的 y 包含火电出力、CHP 电热出力、电锅炉功率、蓄热罐充放热功率。第二阶段可行域 F(x,δ) 通常包含电平衡、热平衡、机组上下限、爬坡约束、蓄热罐容量约束。注意δ 只出现在约束右端项里也就是平衡方程的资源侧。这样处理有个好处子问题对 δ 是线性的枚举候选场景、求最坏场景都只涉及线性规划。第一阶段决策不能等 δ 出来再改这正是“鲁棒”二字的含义。它不像随机优化那样给每个场景一组独立的一阶段变量而是强制所有场景共享同一组 x再让 y 去适配不同场景。两个阶段之间通过约束 A·x B·y b0 D·δ 耦合。2.4 地基代码先写一个确定性电热调度LP验证第二阶段模型能跑不要一上来就写 CCG 主循环。我一般先把第二阶段的确定性 LP 写出来跑通确认电平衡、热平衡、机组耦合约束的系数矩阵没写错。以下是一个极简电热调度模型变量按固定顺序排列火电出力、CHP 电出力、CHP 热出力、电锅炉耗电、电锅炉热出力、蓄热罐充热、蓄热罐放热。import numpy as np from scipy.optimize import linprog # 变量顺序 # 0: p_g 火电出力 # 1: p_chp CHP电出力 # 2: q_chp CHP热出力 # 3: p_eb 电锅炉耗电 # 4: q_eb 电锅炉热出力 # 5: hs_chg 蓄热罐充热功率 # 6: hs_dis 蓄热罐放热功率 P_LOAD 60.0 # 电负荷基准 WIND 20.0 # 风电预测出力 H_LOAD 25.0 # 热负荷基准 def solve_dispatch(wind, heat_load): # 成本系数火电70CHP电80CHP热成本很低 # 电锅炉耗电成本35蓄热罐充放热象征性收一点维护费 c np.array([70, 80, 1, 35, 0.1, 4, 4]) # 不等式约束 A_ub np.array([ [0, 1, 0.6, 0, 0, 0, 0], # 背压式CHPp_chp 0.6*q_chp 100 [0, 0, 0, -0.9, 1, 0, 0], # 电锅炉q_eb 0.9*p_eb写成 0 [0, 0, 0, 0, 0, 6, 6], # 蓄热罐总容量hs_chg hs_dis 6 ]) b_ub np.array([100.0, 0.0, 6.0]) # 等式约束电平衡和热平衡 A_eq np.array([ [1, 1, 0, -1, 0, 0, 0], # p_g p_chp - p_eb P_LOAD - wind [0, 0, 1, 0, 1, -1, 1], # q_chp q_eb hs_dis - hs_chg H_LOAD ]) b_eq np.array([P_LOAD - wind, heat_load]) bounds [ (0, None), # p_g (0, 120), # p_chp (0, 80), # q_chp (0, 30), # p_eb (0, 30), # q_eb (0, 3), # hs_chg (0, 3), # hs_dis ] res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq, boundsbounds, methodhighs) if not res.success: raise RuntimeError(LP infeasible) return res.fun, res.x # 验证一个确定场景 cost, x solve_dispatch(WIND, H_LOAD) print(dispatch cost:, cost) print(p_g %.2f, p_chp %.2f, q_chp %.2f % (x[0], x[1], x[2]))这里的 A_ub 第二行写成-0.9*p_eb q_eb 0是为了用线性不等式逼近等式关系避免重复定义变量。如果后续要把电锅炉效率从 0.9 调成 0.85只需改这个系数。蓄热罐约束里的6是容量上限3是充放功率上限这两个参数直接影响第二阶段的调节空间。跑通这段代码以后才能接着写最坏场景搜索否则一旦 CCG 不收敛你分不清是模型错还是算法错。3. 用 Python 复现两阶段鲁棒 CCG 求解不确定集枚举、最坏场景回填与主循环3.1 为什么选 CCG 而不是 Benders列约束生成在鲁棒调度里更稳两阶段鲁棒优化常用的求解方法是 Benders 分解和列约束生成。Benders 通过对偶变量构造割平面把子问题的信息反馈给主问题列约束生成则直接在主问题里新增一组第二阶段变量和约束每迭代一次就把一个最坏场景“钉”进主问题。我的经验是在电热两阶段鲁棒里优先用 CCG。原因有两个第一第二阶段含电平衡和热平衡只有原变量比只有对偶变量更容易诊断不可行第二CCG 的收敛速度对场景数量不敏感通常十次以内就能把下界抬到上界附近而 Benders 遇到子问题对偶退化会来回振荡。CCG 的流程可以概括为四步。第一步初始化一个名义场景 δ0第二步求解主问题得到第一阶段决策和当前下界第三步把第一阶段决策代入第二阶段在所有候选场景里找最坏场景得到上界第四步如果上下界之差小于阈值停止否则把最坏场景对应的变量和约束加入主问题回到第二步。这里的“场景”不是随机抽样而是不确定集里的候选点。3.2 子问题怎么求最坏场景对偶单层化与顶点枚举的取舍子问题形式是 max_{δ∈W} min_{y∈F(x,δ)} d^T y。直接对 δ 和 y 同时求极值是个 max-min 问题标准做法是拿内层 min 的对偶来替换把子问题变成单层 max。对偶变换本身不复杂内层 LP 对偶后目标变成 max λ^T b(δ)约束是 A^T λ ≤ d其中 b(δ) 是含 δ 的右端项。问题在于 λ 和 δ 会在目标里相乘出现 λ^T Dδ 这种双线性项。如果 δ 是连续多面体变量就需要用大 M 法线性化如果 δ 的候选集是有限顶点问题立刻退化成枚举求最大值。我一般采用的务实方案是把不确定集 W 的顶点枚举出来子问题只在顶点集合里找最坏场景。对“∞-范数 1-范数预算”这种多面体顶点组合有限。比如只有风电偏差和热负荷偏差两个不确定参数时候选顶点就是 (±1,0)、(0,±1)再加一些满足 |δ1||δ2|≤Γ 的组合点。这样做绕开了双线性项代码也更容易排查。3.3 主循环代码候选场景枚举、SP 求解、MP 回填的骨架下面这段代码是上一节 LP 的延伸实现不确定集候选枚举、第二阶段成本求解和最坏场景搜索。它不替代完整主问题但给出了 CCG 里最核心的两个环节找最坏场景和判断要不要继续迭代。import numpy as np from scipy.optimize import linprog from itertools import product # 不确定参数0风电偏差1热负荷偏差 P_LOAD 60.0 WIND_BASE 20.0 WIND_DEV 8.0 H_LOAD 25.0 H_DEV 5.0 GAMMA 1.5 # 1-范数预算 # 候选顶点枚举每个维度取 -1/0/1再按范数预算过滤 def build_candidates(K2, GammaGAMMA): cands [] for vec in product([-1.0, 0.0, 1.0], repeatK): if np.abs(vec).sum() Gamma 1e-9: cands.append(np.array(vec)) return cands # 第二阶段右端项让偏差进入电、热平衡 def make_rhs(delta): wind WIND_BASE delta[0] * WIND_DEV heat H_LOAD delta[1] * H_DEV return np.array([P_LOAD - wind, heat]) # 第二阶段 LP给定偏差场景求最小调度成本 def solve_sp(delta): c np.array([70, 80, 1, 35, 0.1, 4, 4]) A_ub np.array([ [0, 1, 0.6, 0, 0, 0, 0], [0, 0, 0, -0.9, 1, 0, 0], [0, 0, 0, 0, 0, 6, 6], ]) b_ub np.array([100.0, 0.0, 6.0]) A_eq np.array([ [1, 1, 0, -1, 0, 0, 0], [0, 0, 1, 0, 1, -1, 1], ]) b_eq make_rhs(delta) bounds [(0, None), (0, 120), (0, 80), (0, 30), (0, 30), (0, 3), (0, 3)] res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq, boundsbounds, methodhighs) if not res.success: return np.inf, delta return res.fun, delta # 在候选集里找最坏场景 def find_worst_case(candidates): worst_cost -np.inf worst_delta None for delta in candidates: cost, _ solve_sp(delta) if cost worst_cost: worst_cost cost worst_delta delta return worst_cost, worst_delta candidates build_candidates() worst_cost, worst_delta find_worst_case(candidates) print(worst cost:, worst_cost) print(worst delta:, worst_delta)这段代码的核心是make_rhs函数。风电偏差 δ0 出现在电平衡右端热负荷偏差 δ1 出现在热平衡右端二者互不干扰有利于定位问题。build_candidates里用product生成三维笛卡尔积再过滤是因为 K 很小如果 K 超过 10组合数会爆炸那时应该用单纯形法或直接对多面体求顶点不能枚举。运行结果里最坏场景通常不是 (±1,±1)而是预算约束 Γ 限制下的边界组合。如果打印出来的 worst_delta 一直是某个固定顶点说明你的 Γ 或偏差幅度写错了。3.4 CCG 主循环主问题怎么回填上下界怎么收敛完整主问题要加入 η 和每轮最坏场景对应的第二阶段变量副本然后重新求解第一阶段决策。主问题写成线性规划时需要按场景复制 y 变量。以下代码给出 CCG 主循环的骨架并用 PuLP 风格说明约束生成逻辑实际生产中建议用商业求解器的 API 动态加约束。# CCG 主循环示意结构 # 完整 MP 需要用 PuLP/Gurobi 按场景复制第二阶段变量这里简化为循环伪码 def ccg_main(max_iter10, tol1e-3): discovered [np.zeros(2)] # 初始场景取预测值δ0 lower_bound -np.inf upper_bound np.inf for k in range(max_iter): # 1. 在当前 discovered 场景集合上求解主问题。 # 主问题目标min 第一阶段成本 η # 对每个场景 j复制一组第二阶段变量 y_j # 并加入电平衡、热平衡、上限约束得到 # η cost_j(y_j) # 这里省略具体 LP 构造只记录返回结果 mp_cost, x_stage1 solve_mp(discovered) lower_bound mp_cost # 2. 固定 x_stage1在候选场景里求最坏情况。 # 如果第二阶段约束中有 x_stage1 参与需要把它代入 solve_sp worst_cost, worst_delta find_worst_case(candidates) upper_bound min(upper_bound, worst_cost) gap upper_bound - lower_bound print(iter, k, lb, lower_bound, ub, upper_bound, gap, gap) if gap tol: return x_stage1, worst_delta, discovered # 3. 把最坏场景加入 discovered下一次主问题会把该场景的变量列入 discovered.append(worst_delta) raise RuntimeError(CCG not converged in %d iterations % max_iter)主问题的正确回填是整个算法最难的地方。常见的错误是只把最坏场景的成本值加进主问题没有复制对应的第二阶段变量和约束导致主问题算出来的下界偏小。每加入一个场景主问题的变量数和约束数都会增加这是 CCG 和 Benders 最直观的区别。收敛阈值 tol 我一般取 1e-3 或 1e-4 相对值太小会让主问题在最后几次迭代反复增加场景计算时间翻倍。需要说明上面代码里find_worst_case是枚举候选顶点如果候选集很大应改写成带对偶变量的单层 LP。但无论哪种写法核心都是找到让第二阶段成本最大的 δ并把它回填给主问题。4. 落地电热两阶段鲁棒源程序的避坑清单四个高频翻车点4.1 范数预算没写进约束Gamma 形同虚设现象把 Γ 从 0 调到 K算出来的总成本几乎不动鲁棒优化结果和确定性模型一样。原因最常见的是不确定集只写了 ‖δ‖∞ ≤ 1也就是盒式约束忘了加 ‖δ‖1 ≤ Γ。这时每个参数都可以独立取最坏值Γ 根本没有进入模型。另一种可能是make_rhs里没把偏差乘进右端项子问题对 δ 完全不敏感。解决在程序里单独打印build_candidates的结果确认候选顶点数量随 Γ 变化。然后把 Γ0 的成本作为基准逐个增加 Γ画出成本-Γ曲线。正常情况成本应该单调不降。如果曲线走平去查不确定约束矩阵和右端项别急着调算法。4.2 子问题对偶转换翻车报错却指向主问题现象CCG 迭代两三轮后 gap 变成负数或者求解器报“quadratic constraint”错误程序直接中断。原因把 max-min 转成单层 max 时目标里出现 λ^T Dδ 双线性项如果没做线性化就直接丢给 linproglinprog 不支持二次项就会报错。更多时候是强对偶条件不满足比如第二阶段原问题可行域为空子问题返回 inf上界变成负无穷。解决先跑确定性场景 δ0确认第二阶段 LP 可行性。做对偶变换前检查所有变量是否有界、约束是否满足 Slater 条件。如果只是想把程序跑通直接枚举顶点候选集跳过对偶单层化。4.3 热平衡模型写得太“松”鲁棒解反而比确定性解便宜现象理论上鲁棒优化成本应不低于确定性成本但程序输出反过来了第一阶段成本加最坏场景成本比确定性还低。原因热平衡方程里没有接入 δ2热负荷偏差只在目标函数里做了一点惩罚或者蓄热罐的容量约束没有和 δ2 联立。此时第二阶段可以用蓄热罐“免费”吸收热偏差鲁棒约束形同虚设成本自然更低。解决把热负荷偏差写进热平衡右端项并设置蓄热罐充放功率上限和容量上限。检查时用 δ(0,1) 跑一次第二阶段观察蓄热罐出力是否超过了允许范围如果没超说明约束写漏了。4.4 对不上论文里的曲线先是怀疑参数实际是范数阶数不同现象复现结果和论文图线趋势对不上形状一样但成本整体偏低或者最坏场景不同。原因论文可能用的是 2-范数椭球集合你用的是 1-范数预算集合或者论文里 Γ 的定义是相对值你当成了绝对值。94号源程序的注释里如果写了范数阶数先对一下。源程序里如果默认用无预算约束的盒式集合你硬加预算曲线自然不一样。解决先从论文正文里找到不确定集定义的公式确认是 ∞ 1 组合还是纯 2-范数还是用历史数据协方差矩阵生成的椭球。然后把源程序里对应的norm_order参数改成一致再比较确定性模型的成本。差 1% 以内说明实现基本对齐差很多回去查右端项系数。5. 出结果后怎么验证用保守度曲线和场景回放确认两阶段鲁棒优化没写错两阶段鲁棒优化跑通只是第一步真正有价值的是验证结果可信。我拿到一个新模型时首先做成本-Γ曲线。Γ0 是确定性场景成本Γ 从 0 逐渐增大到 K成本曲线应当是一条缓慢上升的单调曲线。如果某一段突然跳变说明那个 Γ 值正好让某个候选顶点进入最坏场景集合这是正常的但如果整体平得像一条直线回到第 4 章第一个坑去查不确定约束。第二个验证手段叫场景回放。取 CCG 收敛后找到的最坏 δ把它对应的风电和热负荷实际值代入确定性模型单独求解一次。这个解的成本应该等于或小于 CCG 给出的上界如果大于说明第二阶段模型在主问题和子问题之间不一致最常见的原因是主问题里复制变量时漏了某个约束。回放时我会把每个决策变量的值打印出来对比主问题里同一个场景的 y 值能精确定位是哪个平衡方程对不上。第三个习惯是保存每一轮迭代的上下界历史。CCG 收敛快不代表每次迭代都有效如果下界连续三轮没有变化而上界在跳说明子问题找到的最坏场景没有改变主问题的最优解问题可能出在目标函数系数上。我会把上下界差异除以确定性成本得到相对 gap控制在 1% 以内再拿去画图否则论文里的曲线很容易被审稿人质疑。这套流程走完你再看“两阶段鲁棒优化-电热-范数程序”这个标题就会发现关键不在代码能跑而在范数集合选得对不对、最坏场景回填准不准、上下界有没有闭合。就我自己的经验来说先跑确定性 LP再枚举顶点找最坏场景最后才上完整 CCG是翻车率最低的路径。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/26 9:11:23 · 阅读 1.2k Python实现电热两阶段鲁棒优化：CCG与范数约束建模
简介面向电气工程、电力系统方向的毕业设计或论文复现场景这份MATLAB源程序包聚焦两阶段鲁棒优化在含电热联合调度的配电网中的应用用于处理分布式电源与负荷不确定性。压缩包仅321KB共12个文件主体为10个.m源程序另有1张系统示意图和1份docx模型说明结构清晰便于按模块阅读。程序完整覆盖拉丁超立方采样、K-means场景聚类、手肘法聚类数确定、1-范数与∞-范数概率置信区间构建以及两阶段分解与列与约束生成算法求解同时结合盒式分解消除储能与机组组合的时间耦合将第二阶段转化为单时间尺度问题可复现动态重构与无功优化协同的综合调度模型。已有168人浏览学习适合电力系统方向研究生和毕设学生用于算法对比、论文复现及进一步扩展研究。1. 两阶段鲁棒优化的电热范数程序先想清楚它到底在优化什么如果你是从知网摸到《两阶段鲁棒优化-电热-范数程序》这篇论文又顺手拿到了对应的94号源程序第一件事别急着跑代码先把“两阶段鲁棒优化”这个词拆干净。它解决的是这么一类问题电负荷、热负荷和风电出力都是波动的你必须在波动发生之前把机组启停、购电合同这类决策定下来等实际值出来以后再调整火电、热电联产机组、电锅炉和蓄热罐出力。前者叫第一阶段决策后者叫第二阶段重调度两个阶段合起来的目标是让最坏情况下的总成本也可接受。电热系统因为热网、蓄热罐的存在天然给第二阶段留了调整余地。而标题里的“范数”不是修辞它决定你把“最坏情况”圈在多大一个边界里。这篇笔记就按“建模、求解、避坑、验证”的顺序讲清楚这类程序怎么写、怎么调、怎么判断结果没写错。2. 电热两阶段鲁棒优化建模用向量范数约束不确定集拆出主问题与子问题2.1 为什么电热系统特别适合两阶段鲁棒热网慢动态提供了“后悔药”电力平衡要求瞬间完成风电一波动火电和储能就得立刻跟上。热力不一样。热水在管道里流动有延迟建筑物墙体本身在蓄热蓄热罐可以把多余的热存起来热负荷的波动不会像电负荷那样直接冲击机组。这意味着第二阶段做出调整时可行域比纯电力系统大得多鲁棒优化做出来的解不会过于保守。常见做法是把电、热两套平衡方程同时放进第二阶段约束风电出力偏差进入电平衡热负荷偏差进入热平衡蓄热罐的充放热功率作为第二阶段变量参与调节。这样一来热力系统的惯性被显式建模成“灵活性资源”而不是被忽略掉。两阶段鲁棒优化在综合能源系统里流行的另一个原因是参数少。随机优化要生成上百个场景才能覆盖风电的分布特性两阶段鲁棒只需要一个不确定集和几个范数参数就能刻画最坏情况。对工程确认来说这种“少参数、保底结论”的思路更容易向调度人员解释。2.2 用向量范数把不确定集圈起来从∞-范数到1-范数预算约束不确定参数一般写成预测值加偏差的形式w_k w_k^0 δ_k·Δw_k其中 δ_k 在 [-1,1] 之间取值Δw_k 是最大偏差幅度。把所有 δ_k 拼成一个向量 δ 后需要限制它的“总扰动预算”不然每个参数都取到最坏值解就会保守到没法用。典型做法是用向量的范数来定义不确定集W { δ | ‖δ‖∞ ≤ 1, ‖δ‖1 ≤ Γ }这个集合的意思是每个参数单独最多偏差一个单位但所有参数的偏差绝对值之和不能超过 Γ。Γ 的取值从 0 到 KK 是不确定参数的个数。Γ0 退化为确定性模型ΓK 等价于盒式不确定集是最保守的情况。实际工程里 Γ 取 K 的一半左右就能兼顾鲁棒性和经济性。与之对比还有纯盒式集合只用 ∞-范数和椭球集合用 2-范数。盒式集合建模简单但最坏场景所有参数同时发生过于极端椭球集合更接近统计置信域但会引入二次约束线性求解器就不能直接用了。所以“∞-范数 1-范数预算”的组合最流行它本质是一个多面体可以保持线性规划结构。范数的阶数选择直接决定源程序里约束怎么写、求解器选哪个。2.3 两阶段结构拆解第一阶段定合同第二阶段做电热再调度两阶段鲁棒优化的通用形式可以写成min_x { c^T x max_{δ∈W} min_{y∈F(x,δ)} d^T y }外层 x 是第一阶段决策比如日前签订的购电电量、机组的最小技术出力、启停状态。内层对每个 δ 求第二阶段最小成本再对外层取最坏情况。这里的 y 包含火电出力、CHP 电热出力、电锅炉功率、蓄热罐充放热功率。第二阶段可行域 F(x,δ) 通常包含电平衡、热平衡、机组上下限、爬坡约束、蓄热罐容量约束。注意δ 只出现在约束右端项里也就是平衡方程的资源侧。这样处理有个好处子问题对 δ 是线性的枚举候选场景、求最坏场景都只涉及线性规划。第一阶段决策不能等 δ 出来再改这正是“鲁棒”二字的含义。它不像随机优化那样给每个场景一组独立的一阶段变量而是强制所有场景共享同一组 x再让 y 去适配不同场景。两个阶段之间通过约束 A·x B·y b0 D·δ 耦合。2.4 地基代码先写一个确定性电热调度LP验证第二阶段模型能跑不要一上来就写 CCG 主循环。我一般先把第二阶段的确定性 LP 写出来跑通确认电平衡、热平衡、机组耦合约束的系数矩阵没写错。以下是一个极简电热调度模型变量按固定顺序排列火电出力、CHP 电出力、CHP 热出力、电锅炉耗电、电锅炉热出力、蓄热罐充热、蓄热罐放热。import numpy as np from scipy.optimize import linprog # 变量顺序 # 0: p_g 火电出力 # 1: p_chp CHP电出力 # 2: q_chp CHP热出力 # 3: p_eb 电锅炉耗电 # 4: q_eb 电锅炉热出力 # 5: hs_chg 蓄热罐充热功率 # 6: hs_dis 蓄热罐放热功率 P_LOAD 60.0 # 电负荷基准 WIND 20.0 # 风电预测出力 H_LOAD 25.0 # 热负荷基准 def solve_dispatch(wind, heat_load): # 成本系数火电70CHP电80CHP热成本很低 # 电锅炉耗电成本35蓄热罐充放热象征性收一点维护费 c np.array([70, 80, 1, 35, 0.1, 4, 4]) # 不等式约束 A_ub np.array([ [0, 1, 0.6, 0, 0, 0, 0], # 背压式CHPp_chp 0.6*q_chp 100 [0, 0, 0, -0.9, 1, 0, 0], # 电锅炉q_eb 0.9*p_eb写成 0 [0, 0, 0, 0, 0, 6, 6], # 蓄热罐总容量hs_chg hs_dis 6 ]) b_ub np.array([100.0, 0.0, 6.0]) # 等式约束电平衡和热平衡 A_eq np.array([ [1, 1, 0, -1, 0, 0, 0], # p_g p_chp - p_eb P_LOAD - wind [0, 0, 1, 0, 1, -1, 1], # q_chp q_eb hs_dis - hs_chg H_LOAD ]) b_eq np.array([P_LOAD - wind, heat_load]) bounds [ (0, None), # p_g (0, 120), # p_chp (0, 80), # q_chp (0, 30), # p_eb (0, 30), # q_eb (0, 3), # hs_chg (0, 3), # hs_dis ] res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq, boundsbounds, methodhighs) if not res.success: raise RuntimeError(LP infeasible) return res.fun, res.x # 验证一个确定场景 cost, x solve_dispatch(WIND, H_LOAD) print(dispatch cost:, cost) print(p_g %.2f, p_chp %.2f, q_chp %.2f % (x[0], x[1], x[2]))这里的 A_ub 第二行写成-0.9*p_eb q_eb 0是为了用线性不等式逼近等式关系避免重复定义变量。如果后续要把电锅炉效率从 0.9 调成 0.85只需改这个系数。蓄热罐约束里的6是容量上限3是充放功率上限这两个参数直接影响第二阶段的调节空间。跑通这段代码以后才能接着写最坏场景搜索否则一旦 CCG 不收敛你分不清是模型错还是算法错。3. 用 Python 复现两阶段鲁棒 CCG 求解不确定集枚举、最坏场景回填与主循环3.1 为什么选 CCG 而不是 Benders列约束生成在鲁棒调度里更稳两阶段鲁棒优化常用的求解方法是 Benders 分解和列约束生成。Benders 通过对偶变量构造割平面把子问题的信息反馈给主问题列约束生成则直接在主问题里新增一组第二阶段变量和约束每迭代一次就把一个最坏场景“钉”进主问题。我的经验是在电热两阶段鲁棒里优先用 CCG。原因有两个第一第二阶段含电平衡和热平衡只有原变量比只有对偶变量更容易诊断不可行第二CCG 的收敛速度对场景数量不敏感通常十次以内就能把下界抬到上界附近而 Benders 遇到子问题对偶退化会来回振荡。CCG 的流程可以概括为四步。第一步初始化一个名义场景 δ0第二步求解主问题得到第一阶段决策和当前下界第三步把第一阶段决策代入第二阶段在所有候选场景里找最坏场景得到上界第四步如果上下界之差小于阈值停止否则把最坏场景对应的变量和约束加入主问题回到第二步。这里的“场景”不是随机抽样而是不确定集里的候选点。3.2 子问题怎么求最坏场景对偶单层化与顶点枚举的取舍子问题形式是 max_{δ∈W} min_{y∈F(x,δ)} d^T y。直接对 δ 和 y 同时求极值是个 max-min 问题标准做法是拿内层 min 的对偶来替换把子问题变成单层 max。对偶变换本身不复杂内层 LP 对偶后目标变成 max λ^T b(δ)约束是 A^T λ ≤ d其中 b(δ) 是含 δ 的右端项。问题在于 λ 和 δ 会在目标里相乘出现 λ^T Dδ 这种双线性项。如果 δ 是连续多面体变量就需要用大 M 法线性化如果 δ 的候选集是有限顶点问题立刻退化成枚举求最大值。我一般采用的务实方案是把不确定集 W 的顶点枚举出来子问题只在顶点集合里找最坏场景。对“∞-范数 1-范数预算”这种多面体顶点组合有限。比如只有风电偏差和热负荷偏差两个不确定参数时候选顶点就是 (±1,0)、(0,±1)再加一些满足 |δ1||δ2|≤Γ 的组合点。这样做绕开了双线性项代码也更容易排查。3.3 主循环代码候选场景枚举、SP 求解、MP 回填的骨架下面这段代码是上一节 LP 的延伸实现不确定集候选枚举、第二阶段成本求解和最坏场景搜索。它不替代完整主问题但给出了 CCG 里最核心的两个环节找最坏场景和判断要不要继续迭代。import numpy as np from scipy.optimize import linprog from itertools import product # 不确定参数0风电偏差1热负荷偏差 P_LOAD 60.0 WIND_BASE 20.0 WIND_DEV 8.0 H_LOAD 25.0 H_DEV 5.0 GAMMA 1.5 # 1-范数预算 # 候选顶点枚举每个维度取 -1/0/1再按范数预算过滤 def build_candidates(K2, GammaGAMMA): cands [] for vec in product([-1.0, 0.0, 1.0], repeatK): if np.abs(vec).sum() Gamma 1e-9: cands.append(np.array(vec)) return cands # 第二阶段右端项让偏差进入电、热平衡 def make_rhs(delta): wind WIND_BASE delta[0] * WIND_DEV heat H_LOAD delta[1] * H_DEV return np.array([P_LOAD - wind, heat]) # 第二阶段 LP给定偏差场景求最小调度成本 def solve_sp(delta): c np.array([70, 80, 1, 35, 0.1, 4, 4]) A_ub np.array([ [0, 1, 0.6, 0, 0, 0, 0], [0, 0, 0, -0.9, 1, 0, 0], [0, 0, 0, 0, 0, 6, 6], ]) b_ub np.array([100.0, 0.0, 6.0]) A_eq np.array([ [1, 1, 0, -1, 0, 0, 0], [0, 0, 1, 0, 1, -1, 1], ]) b_eq make_rhs(delta) bounds [(0, None), (0, 120), (0, 80), (0, 30), (0, 30), (0, 3), (0, 3)] res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq, boundsbounds, methodhighs) if not res.success: return np.inf, delta return res.fun, delta # 在候选集里找最坏场景 def find_worst_case(candidates): worst_cost -np.inf worst_delta None for delta in candidates: cost, _ solve_sp(delta) if cost worst_cost: worst_cost cost worst_delta delta return worst_cost, worst_delta candidates build_candidates() worst_cost, worst_delta find_worst_case(candidates) print(worst cost:, worst_cost) print(worst delta:, worst_delta)这段代码的核心是make_rhs函数。风电偏差 δ0 出现在电平衡右端热负荷偏差 δ1 出现在热平衡右端二者互不干扰有利于定位问题。build_candidates里用product生成三维笛卡尔积再过滤是因为 K 很小如果 K 超过 10组合数会爆炸那时应该用单纯形法或直接对多面体求顶点不能枚举。运行结果里最坏场景通常不是 (±1,±1)而是预算约束 Γ 限制下的边界组合。如果打印出来的 worst_delta 一直是某个固定顶点说明你的 Γ 或偏差幅度写错了。3.4 CCG 主循环主问题怎么回填上下界怎么收敛完整主问题要加入 η 和每轮最坏场景对应的第二阶段变量副本然后重新求解第一阶段决策。主问题写成线性规划时需要按场景复制 y 变量。以下代码给出 CCG 主循环的骨架并用 PuLP 风格说明约束生成逻辑实际生产中建议用商业求解器的 API 动态加约束。# CCG 主循环示意结构 # 完整 MP 需要用 PuLP/Gurobi 按场景复制第二阶段变量这里简化为循环伪码 def ccg_main(max_iter10, tol1e-3): discovered [np.zeros(2)] # 初始场景取预测值δ0 lower_bound -np.inf upper_bound np.inf for k in range(max_iter): # 1. 在当前 discovered 场景集合上求解主问题。 # 主问题目标min 第一阶段成本 η # 对每个场景 j复制一组第二阶段变量 y_j # 并加入电平衡、热平衡、上限约束得到 # η cost_j(y_j) # 这里省略具体 LP 构造只记录返回结果 mp_cost, x_stage1 solve_mp(discovered) lower_bound mp_cost # 2. 固定 x_stage1在候选场景里求最坏情况。 # 如果第二阶段约束中有 x_stage1 参与需要把它代入 solve_sp worst_cost, worst_delta find_worst_case(candidates) upper_bound min(upper_bound, worst_cost) gap upper_bound - lower_bound print(iter, k, lb, lower_bound, ub, upper_bound, gap, gap) if gap tol: return x_stage1, worst_delta, discovered # 3. 把最坏场景加入 discovered下一次主问题会把该场景的变量列入 discovered.append(worst_delta) raise RuntimeError(CCG not converged in %d iterations % max_iter)主问题的正确回填是整个算法最难的地方。常见的错误是只把最坏场景的成本值加进主问题没有复制对应的第二阶段变量和约束导致主问题算出来的下界偏小。每加入一个场景主问题的变量数和约束数都会增加这是 CCG 和 Benders 最直观的区别。收敛阈值 tol 我一般取 1e-3 或 1e-4 相对值太小会让主问题在最后几次迭代反复增加场景计算时间翻倍。需要说明上面代码里find_worst_case是枚举候选顶点如果候选集很大应改写成带对偶变量的单层 LP。但无论哪种写法核心都是找到让第二阶段成本最大的 δ并把它回填给主问题。4. 落地电热两阶段鲁棒源程序的避坑清单四个高频翻车点4.1 范数预算没写进约束Gamma 形同虚设现象把 Γ 从 0 调到 K算出来的总成本几乎不动鲁棒优化结果和确定性模型一样。原因最常见的是不确定集只写了 ‖δ‖∞ ≤ 1也就是盒式约束忘了加 ‖δ‖1 ≤ Γ。这时每个参数都可以独立取最坏值Γ 根本没有进入模型。另一种可能是make_rhs里没把偏差乘进右端项子问题对 δ 完全不敏感。解决在程序里单独打印build_candidates的结果确认候选顶点数量随 Γ 变化。然后把 Γ0 的成本作为基准逐个增加 Γ画出成本-Γ曲线。正常情况成本应该单调不降。如果曲线走平去查不确定约束矩阵和右端项别急着调算法。4.2 子问题对偶转换翻车报错却指向主问题现象CCG 迭代两三轮后 gap 变成负数或者求解器报“quadratic constraint”错误程序直接中断。原因把 max-min 转成单层 max 时目标里出现 λ^T Dδ 双线性项如果没做线性化就直接丢给 linproglinprog 不支持二次项就会报错。更多时候是强对偶条件不满足比如第二阶段原问题可行域为空子问题返回 inf上界变成负无穷。解决先跑确定性场景 δ0确认第二阶段 LP 可行性。做对偶变换前检查所有变量是否有界、约束是否满足 Slater 条件。如果只是想把程序跑通直接枚举顶点候选集跳过对偶单层化。4.3 热平衡模型写得太“松”鲁棒解反而比确定性解便宜现象理论上鲁棒优化成本应不低于确定性成本但程序输出反过来了第一阶段成本加最坏场景成本比确定性还低。原因热平衡方程里没有接入 δ2热负荷偏差只在目标函数里做了一点惩罚或者蓄热罐的容量约束没有和 δ2 联立。此时第二阶段可以用蓄热罐“免费”吸收热偏差鲁棒约束形同虚设成本自然更低。解决把热负荷偏差写进热平衡右端项并设置蓄热罐充放功率上限和容量上限。检查时用 δ(0,1) 跑一次第二阶段观察蓄热罐出力是否超过了允许范围如果没超说明约束写漏了。4.4 对不上论文里的曲线先是怀疑参数实际是范数阶数不同现象复现结果和论文图线趋势对不上形状一样但成本整体偏低或者最坏场景不同。原因论文可能用的是 2-范数椭球集合你用的是 1-范数预算集合或者论文里 Γ 的定义是相对值你当成了绝对值。94号源程序的注释里如果写了范数阶数先对一下。源程序里如果默认用无预算约束的盒式集合你硬加预算曲线自然不一样。解决先从论文正文里找到不确定集定义的公式确认是 ∞ 1 组合还是纯 2-范数还是用历史数据协方差矩阵生成的椭球。然后把源程序里对应的norm_order参数改成一致再比较确定性模型的成本。差 1% 以内说明实现基本对齐差很多回去查右端项系数。5. 出结果后怎么验证用保守度曲线和场景回放确认两阶段鲁棒优化没写错两阶段鲁棒优化跑通只是第一步真正有价值的是验证结果可信。我拿到一个新模型时首先做成本-Γ曲线。Γ0 是确定性场景成本Γ 从 0 逐渐增大到 K成本曲线应当是一条缓慢上升的单调曲线。如果某一段突然跳变说明那个 Γ 值正好让某个候选顶点进入最坏场景集合这是正常的但如果整体平得像一条直线回到第 4 章第一个坑去查不确定约束。第二个验证手段叫场景回放。取 CCG 收敛后找到的最坏 δ把它对应的风电和热负荷实际值代入确定性模型单独求解一次。这个解的成本应该等于或小于 CCG 给出的上界如果大于说明第二阶段模型在主问题和子问题之间不一致最常见的原因是主问题里复制变量时漏了某个约束。回放时我会把每个决策变量的值打印出来对比主问题里同一个场景的 y 值能精确定位是哪个平衡方程对不上。第三个习惯是保存每一轮迭代的上下界历史。CCG 收敛快不代表每次迭代都有效如果下界连续三轮没有变化而上界在跳说明子问题找到的最坏场景没有改变主问题的最优解问题可能出在目标函数系数上。我会把上下界差异除以确定性成本得到相对 gap控制在 1% 以内再拿去画图否则论文里的曲线很容易被审稿人质疑。这套流程走完你再看“两阶段鲁棒优化-电热-范数程序”这个标题就会发现关键不在代码能跑而在范数集合选得对不对、最坏场景回填准不准、上下界有没有闭合。就我自己的经验来说先跑确定性 LP再枚举顶点找最坏场景最后才上完整 CCG是翻车率最低的路径。希望帮到你。本文还有配套的精品资源点击获取