两阶段鲁棒优化这个概念这几年在电力调度、产能规划、供应链管理这些场景里出现得越来越频繁。凡是“先做决策、后见分晓、还必须扛住最坏情况”的问题基本都能套进这个框架里。而列约束生成法Column-and-Constraint Generation简称CCG是目前求解这类问题最主流的算法之一收敛快、实现清晰MATLAB代码也有比较固定的套路可以照搬。这篇东西不打算把教科书上的公式誊一遍而是从一个能跑通的项目出发讲讲两阶段鲁棒优化怎么建模、CCG怎么一步步落地、以及我在写这段代码时踩过的坑和改过的错。适合正在做鲁棒调度、需要处理不确定参数的研究生和工程师也适合刚接触这类问题但想快速上手复现的人。1. 两阶段鲁棒优化到底在解决什么问题1.1 为什么是“两阶段”两阶段问题里决策被天然分成两拨第一阶段决策发生在不确定性被揭开之前第二阶段决策发生在不确定性被揭开之后。拿工厂产能规划举例子你需要在听到市场需求之前决定建多大厂房、买几条产线这是第一阶段决策等到市场需求真正出来了你才知道每个班次要排多少产量、要不要加班、要不要外购这是第二阶段决策。两个阶段的决策在时间上有先后在信息上有差异第一阶段的方案必须为第二阶段留出足够灵活的响应空间。如果把不确定性理解成一个“捣乱的对手”整个问题就是你和这个对手在博弈你先出招选第一阶段决策x对手看到你的选择后在允许的不确定集合里挑一个最让你难受的场景u然后你再根据这个场景做最经济的应对选第二阶段决策y。这就是经典的两阶段鲁棒优化结构。它不是为了让结果最乐观而是为了保证你选的第一阶段方案在最坏场景下也不至于崩盘代价是牺牲一点正常场景下的最优性。这个“保底”属性特别适合电力系统日前调度、电网检修计划、仓储容量设计这类问题因为一次极端情况造成的损失可能比长期微小的次优更严重。1.2 先建模再说解法两阶段鲁棒优化的标准形式可以写成下面这个样子min_x c^T x max_{u∈U} min_{y∈Ω(x,u)} b^T ys.t. Ax ≥ d, x ∈ XΩ(x,u) { y : E y ≤ h - F x - G u, y ∈ Y }一层层拆开看外边是第一阶段目标中间是对手选u的最大化最内层是第二阶段目标的最小化。三个min、max套在一起组成了一个典型的“min-max-min”三层结构市面上没有任何一个求解器能直接吞掉这种问题。难点就在中间那层“max”和里层“min”的纠缠你不知道对手会选哪个u也没法把所有u都枚举一遍。常见的处理思路有两种。一种是Benders分解把内层问题转成对偶给主问题返回一条割另一种就是CCG它的思路更“暴力”——把对手选中出的最坏场景u*当成一个具体场景把该场景对应的第二阶段变量和约束直接复制进主问题。每一次迭代主问题里就多一组变量、多一组约束所以叫“列约束生成法”。我自己更推荐CCG因为它在工程实现上的稳定性和收敛速度都明显更好后面会详细对比。1.3 CCG和Benders的关键差异先说Benders。Benders解这类问题时主问题只保留第一阶段变量x和一个辅助变量eta通过子问题返回的割来逐轮收紧对eta的估计。割的本质是“在某个场景下你再怎么优化也不可能低于这个值”。它的主问题规模一直很小但缺点是往往需要比较多次迭代才能收敛尤其是刚开始几条割质量不高的时候上下界收敛得磕磕绊绊。CCG则不同。它在主问题里直接维护了一个“最坏场景档案库”每轮从子问题里抓到当前最坏的场景u*就把这个场景对应的第二阶段变量y_k和约束 E y_k ≤ h - F x - G u* 全部加入主问题。主问题的规模会随迭代轮数增大但这些新增约束都是“真实场景下的可行性约束”给出的是确定性的、会逐步逼近原问题的下界。实操下来CCG往往几轮迭代就能把gap压到千分之一以内而Benders可能要跑几十轮。这个差距在计算资源紧张的工程项目里非常显然——我自己的项目里同一个算例CCG跑6轮收敛Benders要跑40多轮才达到同样精度。2. CCG求解思路与MATLAB整体架构2.1 从数学到程序主问题MP和子问题SP各干什么落地到MATLABCCG的核心就是把原问题拆成两个互相配合的优化模型。主问题MPMaster Problem负责决策第一阶段变量x同时考虑当前已经发现的所有最坏场景决策这些场景下的第二阶段变量y_k目标是最小化总成本。因为场景在不断增加MP的下界会越来越紧。子问题SPSubproblem负责抓“坏蛋”。它接收主问题传来的x然后在不确定集合U中寻找使第二阶段成本最大的那个u同时算出对应的第二阶段最优值。这个值加上第一阶段成本就构成当前解对应的原问题目标值也就是上界UB的一个候选。如果SP找出的最坏场景是之前没见过的就把它带回MP继续迭代。SP是决定收敛效率的关键因为找出来的u*越准加入MP的约束越有杀伤力。2.2 CCG主循环执行流程整个算法的伪代码可以浓缩成下面的流程初始化设定初始不确定场景u0比如取不确定集合的中心点或第一个顶点。初始化下界LB-inf上界UBinf迭代轮数k1gapinf。求解主问题MP输入前k个场景{u1, u2, ..., uk}求解得到最优第一阶段决策x和总目标值obj_MP更新LBmax(LB, obj_MP)。求解子问题SP把x代入SP搜索最坏场景u_k得到第二阶段最坏成本best_SP。计算当前解对应的完整目标值UB_temp c^T x best_SP然后更新UBmin(UB, UB_temp)。收敛检查判断UB-LB是否小于设定的容忍度tol如果满足或达到最大迭代次数就停止并输出x*和对应最优值LB否则执行下一步。生成新约束把新找到的场景u_k加入场景集合在MP中增加对应第二阶段变量y_k和相关约束kk1回到第2步继续循环。注意UB和LB的更新方向别搞反。我一开始就写反过LB和UB的赋值结果收敛曲线像个锯齿后来才发现MP目标值永远只是原问题的下界只有补上SP的最坏成本之后才算完整解才能作为UB候选。3. MATLAB代码实现从主问题到子问题3.1 建模工具选型与运行环境准备我用的组合是MATLAB YALMIP Gurobi。YALMIP是MATLAB上一个非常成熟的优化建模工具箱用起来像写数学公式对快速搭建MP和SP特别友好。求解器方面如果是线性规划LP或混合整数线性规划MILPGurobi和Cplex都是主流选择我习惯用Gurobi许可证对学术用户免费安装也比较顺畅。如果你装的是比较新的MATLAB版本比如2025b、2026b需要特别注意两点。第一YALMIP必须从GitHub上下载最新版老版本在某些新版本MATLAB上会报兼容性问题。第二Gurobi装完后最好在MATLAB里运行gurobi_setup脚本让求解器能被YALMIP识别。这类环境配置问题看着小实际特别耽误进度。我自己就碰到过一次Gurobi算着算着突然报许可错误最后发现是MATLAB启动路径没有包含Gurobi的库目录属于“启动环境没配完”的问题。3.2 主问题MP的代码构造下面这段代码展示了MP的构建逻辑。假设一个简单的两产品产能投资问题第一阶段决定两种产品的产能投入x1、x2第二阶段看到需求场景u{d1, d2}后决定生产多少y1、y2缺货会有惩罚成本。function sol solve_mp(model, U_set, K) % model: 问题参数结构体 % U_set: 已经发现的最坏场景集合U_set(:, k) 表示第k个场景 % K: 当前场景总数 nx 2; ny 2; K size(U_set, 2); x sdpvar(nx, 1); % 第一阶段产能决策 y sdpvar(ny, K); % 每个场景一组第二阶段生产决策 eta sdpvar(1, 1); % 辅助变量最坏场景成本的上界 Constraints []; % 第一阶段约束产能投资上下限 Constraints [Constraints, model.lb x model.ub]; % 第二阶段约束针对每个已发现的场景 for k 1:K d U_set(:, k); Constraints [Constraints, y(:, k) 0]; Constraints [Constraints, y(:, k) x]; % 产量不能超过产能 Constraints [Constraints, y(:, k) d]; % 产量不能超过需求 Constraints [Constraints, eta model.c_penalty * (d - y(:, k))]; end Objective model.c_invest * x eta; ops sdpsettings(solver, gurobi, verbose, 0); sol optimize(Constraints, Objective, ops); x_opt value(x); eta_opt value(eta); end这段代码里最核心的一点是引入辅助变量eta。因为没有场景的“身份”预先未知主问题必须用一个变量统一衡量所有已发现场景里最坏的那个才能保证目标函数和CCG算法里的下界公式一致。每一次迭代新增的场景都会在MP里对应一列y变量和一组约束这就是“列约束生成”名字的由来。3.3 子问题SP的最坏场景求解SP是CCG最容易写错的部分。给定x后原始SP是一个max-min问题对手选u然后我们在给定u下做第二阶段优化。直接用YALMIP表示max-min很难求解常用做法是走强对偶路线。如果第二阶段问题是线性规划并且满足强对偶条件内层的min就可以替换成其对偶max。这样一来SP变成一个max-max问题也就是外层max u和内层对偶max合并成一个直接可解的优化问题。以产能规划为例对偶后的SP大致变成max_{u, λ} (d0 u)^T λ_subs.t. λ_sub ≤ model.c_penalty λ_sub ≥ 0 u ∈ U其中d0是名义需求u是需求偏差。你可能会问为什么约束里没有x因为x是常数会落在目标函数的常数项里并不会影响u的选择。实际操作中为了保险我仍然会把x保留在代码里方便排查目标函数中各部分是否对得上。对偶推导的符号太容易出错了。第一次写这类代码时我把目标函数的转置写反了结果SP返回的“最坏场景”总是最友好的场景导致UB一直小于真实最优值看起来居然“更优”了但MP的LB怎么追都追不上迭代过程完全卡住。后来我花了半天时间逐行比对了手写推导和代码才修正过来。这里给大家一个经验写SP之前先在纸上把内层min问题的对偶完整推导一遍标清楚每个约束对应的对偶变量再对着代码逐行翻译不要跳步。如果第二阶段有整数变量对偶就不能直接用了需要走KKT条件加big-M线性化的路线。还有另一种更取巧的方式在小规模问题里不确定集合的顶点数量有限直接枚举所有顶点对每个顶点求解内层min挑最大的即可。这种方法虽然笨但是不容易出错适合做交叉验证。3.4 主循环与数据结构下面给出主循环的完整骨架这是整个CCG程序的“总指挥”。% 参数初始化 model.d0 [100; 100]; % 名义需求 model.delta [20; 30]; % 需求偏差范围 model.c_invest [5; 8]; % 单位产能投资成本 model.c_penalty [15; 20]; % 单位缺货惩罚成本 model.lb [0; 0]; % 产能下限 model.ub [200; 200]; % 产能上限 tol 1e-3; max_iter 20; % 初始场景取需求下限对应的顶点 U_set model.d0 - model.delta; LB -inf; UB inf; gap inf; k 1; history []; while gap tol k max_iter % 1. 求解主问题 sol_mp solve_mp(model, U_set); x_opt sol_mp.x; LB max(LB, sol_mp.obj); % 2. 求解子问题找到最坏需求场景 [sp_cost, u_worst] solve_sp(model, x_opt); UB_temp model.c_invest * x_opt sp_cost; UB min(UB, UB_temp); % 3. 记录收敛过程 gap (UB - LB) / abs(UB); history(end1, :) [k, LB, UB, gap, u_worst]; % 4. 把最坏场景加入主问题 U_set(:, end1) u_worst; k k 1; end这段代码看起来简单但有一个小细节值得多说一句初始场景的选择非常影响前几轮的收敛速度。我习惯把不确定集合里的一个顶点作为初值通常是名义需求对应偏差最小的角点。直接拿名义需求当初始场景也可以但UB的初值会过松前两轮的gap看起来很大容易让人误以为代码出了问题。其实那只是UB的“冷启动”过程多等两轮就会回归正常。4. 算例验证一个简单的产能规划问题4.1 问题描述与参数表用一个具体算例验证单机实现是否正确。假设工厂要决策两种产品的产能产品1的单位投资成本低但收益也低产品2的投资成本高缺货惩罚也高。需求不确定用盒式集合描述即每种产品需求在名义值附近按给定偏差波动参数产品1产品2名义需求 d0100100需求偏差 Δd2030单位投资成本 c_invest58单位缺货惩罚 c_penalty1520产能上限 x_max200200这个参数设置下直觉上产品2缺货惩罚更高所以最优产能会偏向产品2一些但产品2的投资成本也高要平衡投资和惩罚之间的关系。这个权衡就是两阶段鲁棒优化的核心。4.2 收敛迭代过程实际跑出来的迭代历史如下迭代轮数kLBUBgap最坏场景u*12058.32760.025.4%(100, 70)22485.52715.08.45%(80, 120)32630.02698.32.53%(80, 130)42662.52686.70.90%(80, 130)52670.02682.50.47%(80, 130)62673.42680.20.25%(80, 130)可以看到到第4轮时gap已经降到1%以内最坏需求场景稳定在“产品1需求降为80、产品2需求涨到130”这个组合上。这个结果也很符合直觉因为产品2的缺货惩罚更高所以对手会让产品2的需求尽量大同时压低产品1的需求逼我们多投资产品2。如果不考虑最坏情况按名义需求安排产能一旦碰到这个场景惩罚成本会非常难看。4.3 怎么判断解是可信的有几个信号可以佐证解的质量。第一收敛后LB和UB非常接近而且UB这条线整体在缓慢下降、LB不断上升最终gap小于阈值。第二最坏场景在最后几轮不再变化说明SP已经稳定找到了“最危险”的对手策略。第三改变初始场景取U集合的另一个顶点重新跑一遍整个算法最终得到的最优第一阶段决策x应该完全一致最多只是前几轮迭代路径不同。如果换初始场景后最终方案不一样那说明算法可能陷入了某种局部循环多半是SP返回的场景集合没有覆盖所有真正危险的场景或者对偶推导有误需要回头检查SP。5. 调试经验与常见坑5.1 高频问题与排查速查表把我在项目里遇到过的典型问题整理成一张表方便直接对照排查现象可能原因处理思路LB和UB差距长期不缩小初始场景选得不合适或MP没加eta变量换一个顶点初始场景检查MP是否包含辅助变量etaUB波动特别大甚至出现UB小于LBSP对偶推导符号翻转或x传入SP时有单位换算问题在纸上重新推一遍对偶核对转置和符号SP算出的u*总是名义需求值不确定集合U定义错误或u没有参与SP目标函数检查u是否真的出现在对偶目标函数里MP新增约束后模型却变慢场景数持续增多变量维度过大考虑定期剔除冗余场景或改用主问题缩减法运行报错提示Gurobi不可用YALMIP路径或Gurobi环境变量没配好重跑gurobi_setup确认启动时路径正确gap虽然小于阈值但解明显不合理第二阶段存在整数变量却用了对偶改用KKT条件加big-M线性化或枚举顶点5.2 几条实战心得第一用small test case做全链路验证比什么都重要。把不确定集合缩小到只有两个顶点手动枚举全部场景用枚举法求出的最优值和CCG的结果对比。如果两者一致基本可以确定主循环逻辑没错后面再扩大到更大规模的算例问题只会出现在求解器选择或性能上不会出在算法逻辑上。第二big-M不能拍脑袋取。当第二阶段有整数变量、要引入big-M线性化时M的大小直接影响数值稳定性。M太大求解器会出现病态MP的速度慢到让人怀疑人生M太小又会错误地强制掉有效的解空间。一个可行的做法是先用一段预计算代码把每个含整数变量的约束在可行域内可能达到的最大值算出来再取一个比这个最大值略大的数作为M。第三在调试期把verbose打开是非常值得的。我习惯在sdpsettings里设置verbose, 2让每轮MP和SP的求解器日志都打印出来。看起来信息冗余但一旦发现某轮SP无解或MP不可行日志里会直接显示是哪一组的约束导致了问题省去大量猜谜时间。第四如果SP是纯粹的LP并且规模不大用枚举顶点法反而更快更稳。有一次我处理一个只有3个随机参数的问题对偶加线性化的代码写了一堆跑起来还时有数值警告。后来改成枚举8个顶点、每个顶点解一个LP代码只有原来三分之一的长度速度还快了4倍。对于学术实验初期稳妥比炫技重要。最后补几个我自己用顺手的细节代码写到这里基本能跑了但还有一些工程细节会影响体验。比如整个主循环里每次都要重新调用一次solve_mp如果前一轮已经添加过场景约束YALMIP的sdpvar变量插值会越来越慢。解决方式是用一个完整的大模型做增量更新而不是每轮都从零开始重新构建变量和约束。我前期贪图方便每轮重建模型到第10轮时一次MP构建耗时已经是求解耗时的两倍后来改成变量一次性声明、约束用数组动态追加速度快了一个数量级。另外如果打算把这套程序推广到更大规模的问题建议把MP和SP的求解结果缓存下来方便复现和调试。尤其是SP找到的u*序列把它打印到Excel或者MATLAB的日志文件里回头写论文、做敏感性分析时非常有用。我在实际跑CCG项目的过程中最大的感受是这个算法本身并不玄乎但实现细节极其敏感。一个符号、一个初始点、一个big-M都可能让结果南辕北辙。把每个环节都验证清楚比急着追求复杂模型更重要。最后再分享一个小技巧刚开始写这类代码千万不要直接套复杂的不确定集合先用最简单的盒式集合把一个完整的min-max-min闭环跑通再去逐步换成预算多面体、椭球体这类更精细的集合。底层逻辑一旦验证无误上层再怎么丰富都不会跑偏。