拓十年匠心定制 · 商业建站与技术教学双线并行 咨询热线:400-886-1026 service@lmnt.cn
ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

两阶段鲁棒优化与CCG算法:Python+Gurobi复现与避坑指南

两阶段鲁棒优化与CCG算法:Python+Gurobi复现与避坑指南

简介:这份资源面向运筹优化方向的研究生、科研人员与工程实践者,聚焦两阶段鲁棒优化中CCG列生成与Benders分解两类主流求解框架,帮助读者在不确定环境下构建最坏情况性能可接受的决策模型。压缩包共5个文件,含3篇PDF与2个MATLAB脚本,整体约1.44MB,PDF用于梳理Benders分解与列生成对比、算法综述及两阶段鲁棒问题求解思路,m文件则给出可运行的代码案例,便于对照论文复现与二次改编。资源已有2285人学习下载,说明其在同类资料中具备一定参考价值。读者可从中获得CCG切割平面迭代逼近最优解的编程实现、Benders切割分离主问题与子问题的建模流程,以及基于YALMIP调用求解器的完整代码结构,既能用于学术复现,也可为电力调度、综合能源等含不确定参数的工程优化问题提供可扩展的求解模板与排错参考。

1. 两阶段鲁棒优化与 CCG:为什么它成了调度与配置问题的标配解法

做电力调度、微网容量配置、供应链选址这类问题,只要数据里带「不确定」,你迟早会撞上两阶段鲁棒优化。它的核心承诺很硬:在最坏的不确定实现下,仍然保证一个可行的、代价可控的决策。而 CCG(Column-and-Constraint Generation,列与约束生成)就是求解这类模型最常用的算法骨架,配合 Benders 分解处理第二阶段的对偶子问题,形成一套能落地、能扩展、能复现的代码结构。

这套方法解决的痛点很具体:场景法需要预先枚举大量不确定场景,规模一上来就爆;而 CCG 只在与「最坏场景」交手的过程中逐步加列、加约束,迭代到收敛。适合谁?适合已经会写基本线性规划、懂一点对偶理论,但被「不确定集 + 两阶段结构」卡住的工程师和研究生。下面我按自己复现论文时的顺序,把模型、代码、参数和坑一次讲透。

2. 两阶段鲁棒优化的模型骨架与 CCG 迭代逻辑

2.1 为什么是 min-max-min,而不是直接场景法

两阶段鲁棒优化的标准形式长这样:

min_x c^T x + max_{u∈U} min_{y∈Ω(x,u)} b^T y

第一层min_x是「here-and-now」决策,比如建多少容量、开哪些机组,必须在不确定量揭晓前定死。第二层max_{u∈U}是自然或对手在不确定集 U 里挑最坏场景。第三层min_y是「wait-and-see」决策,比如实际调度出力,场景揭晓后才做。

场景法是把 U 离散成有限个场景,一次性写进一个大 LP。问题是场景数随不确定维度指数增长,而且很多场景根本不会成为最坏场景,白算。CCG 的思路是:我不枚举,我迭代。主问题(MP)只带当前已知的少数最坏场景,解出一个下界;子问题(SP)在给定 x 下找当前最坏场景,如果找到的场景让目标变差,就把它作为一列(新变量 y 和新约束)加进 MP,再解。上下界收敛即停。

这个「加列加约束」就是名字的由来。Benders 分解在这里的角色是:SP 本身是 max-min 结构,用对偶把内层 min 转成 max,整个 SP 变成一个单层 max 问题,可以直接求解或再对偶化成 MILP。

2.2 不确定集怎么选:盒式、多面体还是预算约束

不确定集 U 的建模直接决定鲁棒性和保守度。常见三种:

类型形式特点适用
盒式u ∈ [u_min, u_max]最简单,最保守快速验证
多面体线性约束集合可刻画相关性有历史数据
预算约束Σu_i - ū_i/Δ_i ≤ Γ

预算约束(Bertsimas-Sim)是我最推荐的:Γ 控制同时偏离预测值的维度数,Γ=0 退化成确定性,Γ=维度数退化成盒式。调 Γ 就是在鲁棒性和经济性之间找平衡,这个参数在论文复现里几乎必调。

2.3 CCG 主问题与子问题的数学拆解

主问题 MP(第 k 次迭代):

min_{x, y^1..y^k, η} c^T x + η s.t. Ax ≤ d η ≥ b^T y^l ∀l=1..k Ex + Fy^l ≤ h + Gu^l ∀l=1..k x ∈ X, y^l ∈ Y

注意 η 是第二阶段最坏代价的估计,约束η ≥ b^T y^l保证它不低于任何已知场景的代价。MP 给的是下界 LB。

子问题 SP(给定 x*):

max_{u∈U} min_y b^T y s.t. Fy ≤ h + Gu - Ex*

内层 min 用对偶转成 max,与外层 max 合并(或保持双层用 KKT)。SP 给的是上界 UB。当 UB - LB ≤ ε 收敛。

提示:SP 若保持 max-min 双层,可以用对偶把内层变成 max,两个 max 合并成单层;也可以对偶后线性化 big-M。前者快但要求内层强对偶成立,后者稳但引入大 M 数值问题。

3. 用 Python + Gurobi 复现 CCG 主问题与子问题

3.1 环境与依赖:最小可跑清单

我一般用 Python 3.9+ 配 Gurobi,学术 license 够用。如果不想装 Gurobi,可以换 PuLP + CBC,但子问题对偶后带 big-M,CBC 收敛会慢很多,调试阶段建议还是 Gurobi。

pip install gurobipy numpy pandas matplotlib

代码结构我习惯分四个文件:model.py放参数和集合,mp.py主问题,sp.py子问题,ccg.py主循环。这样扩展改编时改哪块很清楚。

3.2 主问题建模:变量、约束与 η 的处理

import gurobipy as gp from gurobipy import GRB def build_mp(data, scenarios): # scenarios: 已累积的最坏场景列表,每个是 u 向量 m = gp.Model("MP") # 第一阶段变量:容量配置 x = m.addVars(data['I'], lb=0, name="x") # 第二阶段变量:每个已知场景一组 y = {} for k, u in enumerate(scenarios): y[k] = m.addVars(data['J'], lb=0, name="y_%d" % k) eta = m.addVar(lb=-GRB.INFINITY, name="eta") # 目标:第一阶段成本 + 最坏第二阶段代价估计 m.setObjective( gp.quicksum(data['c'][i] * x[i] for i in data['I']) + eta, GRB.MINIMIZE) # 第一阶段约束 for i in data['I']: m.addConstr(x[i] <= data['x_max'][i]) # 每个场景的第二阶段约束 + eta 下界 for k, u in enumerate(scenarios): for j in data['J']: m.addConstr( gp.quicksum(data['F'][j][i] * x[i] for i in data['I']) + y[k][j] >= data['h'][j] + data['G'][j] * u[j]) m.addConstr(eta >= gp.quicksum(data['b'][j] * y[k][j] for j in data['J'])) return m, x, y, eta

逻辑说明:eta是第二阶段最坏代价的代理变量,每加一个场景就加一条eta >= b^T y^k,逼着 eta 不低于所有已知场景的代价。x是共享的第一阶段变量,y[k]是每个场景独立的第二阶段变量——这正是 CCG「加列」的体现:新场景进来就新增一组 y 变量和对应约束。

参数说明:data['c']第一阶段单位成本,data['b']第二阶段单位成本,data['F']是 x 对第二阶段约束的耦合矩阵,data['G']是不确定量对约束的系数。这些矩阵的维度要和 I、J 对齐,复现论文时最容易在这里维度对不上。

3.3 子问题对偶化:把 max-min 变成可解形式

子问题给定 x*,求最坏 u。内层 min 的对偶:

def solve_sp(data, x_star): # 对偶变量 lambda 对应第二阶段约束 m = gp.Model("SP") lam = m.addVars(data['J'], lb=0, name="lam") u = m.addVars(data['J'], lb=data['u_min'], ub=data['u_max'], name="u") # 对偶目标:max lambda^T (h + G u - F x*) obj = gp.quicksum(lam[j] * (data['h'][j] - gp.quicksum(data['F'][j][i] * x_star[i] for i in data['I'])) for j in data['J']) \ + gp.quicksum(lam[j] * data['G'][j] * u[j] for j in data['J']) m.setObjective(obj, GRB.MAXIMIZE) # 对偶约束:lambda^T F <= b for j in data['J']: m.addConstr(lam[j] <= data['b'][j]) # 预算约束不确定集 m.addConstr(gp.quicksum(u[j] for j in data['J']) <= data['Gamma']) m.optimize() u_worst = [u[j].X for j in data['J']] sp_val = m.ObjVal return u_worst, sp_val

逻辑说明:内层min b^T y s.t. Fy ≤ h + Gu - Ex*的对偶是max λ^T(h + Gu - Ex*) s.t. λ^T F ≤ b, λ ≥ 0。因为内层是 min 且约束是 ≤,对偶变量 λ 非负。外层 max 和这个 max 同向,直接合并成单层 max,不用 big-M,数值上干净很多。

参数说明:Gamma是预算约束参数,控制鲁棒保守度。u_min、u_max是不确定量边界。注意G[j]如果是不确定量的系数,u 的符号要和它匹配,否则最坏场景会取反。

3.4 CCG 主循环:上下界收敛与场景累积

def ccg(data, eps=1e-4, max_iter=50): scenarios = [data['u_nominal']] # 初始用标称场景 LB, UB = -GRB.INFINITY, GRB.INFINITY for it in range(max_iter): # 解主问题 mp, x, y, eta = build_mp(data, scenarios) mp.optimize() x_star = [x[i].X for i in data['I']] LB = mp.ObjVal # 解子问题找最坏场景 u_worst, sp_val = solve_sp(data, x_star) # 真实上界 = 第一阶段成本 + 最坏第二阶段代价 first_stage = sum(data['c'][i] * x_star[i] for i in data['I']) UB = first_stage + sp_val print("iter %d: LB=%.4f UB=%.4f gap=%.4f" % (it, LB, UB, UB - LB)) if UB - LB <= eps: break scenarios.append(u_worst) return x_star, LB, UB

逻辑说明:每轮先解 MP 拿 LB 和 x*,再用 x* 解 SP 找最坏场景和 UB。如果 gap 没到 eps,就把最坏场景加进 scenarios,下一轮 MP 会多一组 y 变量和约束。这就是「列与约束生成」的完整闭环。

参数说明:eps收敛容差,工程上 1e-4 到 1e-3 都行,太小会卡在数值噪声上。max_iter是保险丝,正常 10 到 20 轮收敛,超过 50 轮基本是模型或对偶写错了。初始场景用标称值u_nominal,也可以用一个极端场景加速。

4. 复现论文时的避坑与排查清单

4.1 上下界不收敛,gap 来回震荡

现象:LB 和 UB 交替上升下降,gap 不单调缩小,跑几十轮还在 0.1 以上。

原因:最常见的是子问题对偶符号写反,导致 SP 求出的不是真正的最坏场景,UB 算出来偏小甚至低于 LB。其次是 MP 里 eta 的下界约束漏了某个场景,或者 y 变量没按场景独立。

解决:先固定 x 手动验证 SP——把 x* 代进去,用暴力枚举 u 的顶点算真实最坏值,和对偶结果对比。如果对不上,检查对偶约束方向:内层 min 配 ≤ 约束,对偶是 max 配 λ ≥ 0,λ^T F ≤ b。再检查 UB 计算是不是first_stage + sp_val,别把 eta 当 UB。

4.2 预算约束 Γ 设太大导致无解或极保守

现象:Γ 调到接近维度数时,MP 可行但目标巨大,或者直接 infeasible。

原因:Γ 太大等价于盒式不确定集,最坏场景让所有不确定量同时取极端,第二阶段约束被顶穿。如果模型本身没有足够的第二阶段调节能力,就会 infeasible。

解决:Γ 从 0 开始逐步加,观察目标增长曲线,选拐点。工程上 Γ 取维度数的 20% 到 50% 比较常见。另外确认不确定集里 u 的上下界和预算约束不冲突,比如 u_min 全正但预算约束写成 Σu ≤ Γ 就会出问题。

4.3 big-M 线性化后数值不稳定

现象:SP 用 KKT + big-M 时,求解器报 numerical trouble,或者解出来的 u 在边界上跳。

原因:big-M 取太大,LP 松弛的系数矩阵条件数爆炸。或者对偶变量和原变量乘积项线性化时 M 没取紧。

解决:优先用 3.3 的对偶合并法,根本不需要 big-M。如果必须用 KKT,M 取约束实际最大可能值的 1.5 倍左右,别拍脑袋写 1e6。Gurobi 可以开m.setParam('NumericFocus', 3)缓解。

4.4 场景累积后 MP 规模膨胀

现象:迭代到 30 轮以上,MP 变量数和约束数线性增长,求解时间从秒级涨到分钟级。

原因:每个场景一组 y 变量,场景多了 MP 就大。这是 CCG 的固有代价。

解决:加场景筛选——如果新场景和已有场景的 u 向量距离小于阈值,就不加。或者用 bundle 方法,只保留对 eta 约束起作用的活跃场景。工程上迭代 20 轮内收敛的话,规模通常还能接受。

4.5 复现论文时目标函数差一个常数

现象:代码跑出来的最优值和论文表格对不上,差一个固定量。

原因:论文里的目标函数可能省略了常数项,或者第一阶段成本的定义范围不同(比如是否含固定成本)。也可能是单位不一致,论文用万元你用了元。

解决:先对齐目标函数的每一项,把论文的公式逐项抄下来和代码对照。特别注意求和范围——论文写 Σ_{i∈I} 你写成 Σ_{i∈I∪I0} 就会差。单位统一到同一量纲再比。

5. 从能跑到好用:CCG 的扩展改编与验证技巧

复现只是起点,真正值钱的是把这套骨架改成自己的问题。我一般做三件事:换不确定集、加整数变量、做收敛加速。

换不确定集最直接。把 3.3 里的预算约束换成多面体或数据驱动的不确定集,SP 结构不变,只改 u 的可行域约束。如果要用 Wasserstein 球这类数据驱动集,SP 会多一层,需要重新对偶,但 CCG 主循环不动。

加整数变量要小心。如果第二阶段有整数变量,内层 min 不能直接对偶,强对偶不成立。常见做法是把第二阶段整数变量松弛,或者用嵌套 CCG。我一般先松弛跑通,看 gap 能不能接受,不行再上嵌套。

收敛加速我常用两个技巧。一是 warm start:第一轮 MP 除了标称场景,再加两个极端场景(u 全取上界、全取下界),通常能少迭代 3 到 5 轮。二是自适应 eps:前几轮 eps 放宽到 1e-2,接近收敛再收紧到 1e-4,避免早期在数值噪声上浪费轮次。

验证方法我固定用三个:

验证项方法通过标准
SP 正确性暴力枚举 u 顶点对比误差 < 1e-6
鲁棒可行性蒙特卡洛采样 1000 场景全部可行
保守度对比 Γ=0 和 Γ=max目标单调递增

蒙特卡洛那步特别重要。CCG 只保证在最坏场景下可行,但实际中不确定量不会总取最坏。采样验证能看出方案在典型场景下的表现,也能发现不确定集建模是否过松。

最后说个血泪经验:别一上来就调 Γ 和 eps,先把 SP 的对偶验证对。我见过太多人 gap 不收敛,折腾半天参数,结果是子问题符号写反。对偶验证花十分钟,能省两天。另外代码里所有矩阵维度用 assert 卡住,复现论文时维度对不上是最常见的翻车点,早报错早安心。

这套 CCG + Benders 的骨架我从微网配置做到供应链选址,改的只是矩阵和不确定集,主循环几乎没动过。如果你也在做两阶段鲁棒优化,建议先把 3.4 的主循环跑通,再逐步替换成自己的模型,比从头写稳得多。希望帮到你。

本文还有配套的精品资源,点击获取

返回列表