做生产调度、供应链规划或者电力系统决策的朋友应该都遇到过这种场景需求不确定、来风不确定、价格不确定但投资决策必须在这些不确定实现之前拍板。一旦拍错板后面怎么调度都补救不回来。这几年我处理这类问题时用得最顺手的工具就是两阶段鲁棒优化加列约束生成法Column-and-Constraint GenerationCCG。这篇博客就把我在MATLAB里用CCG求解两阶段鲁棒优化问题的完整思路、代码和踩坑记录都写出来适合正在做鲁棒优化课题的研究生以及想给决策模型加抗不确定性能力的工程师参考。两阶段鲁棒优化问题在教科书里写得挺抽象但落到实际场景就一句话第一阶段的决策要先定下来第二阶段的决策可以随机应变而中间的不确定性专门挑你最难受的情况出现。CCG就是暴力地跟这个最难受的情况硬碰硬反复把最坏场景找出来、塞进主问题里逼近真解。下面我就用一个投资-分配的小算例从建模、算法推导到MATLABYALMIP代码全流程走一遍。1. 一个投资决策问题怎么变成两阶段鲁棒优化1.1 一个具体的投资-分配场景假设你手上有4个候选仓库位置每个仓库有固定建设成本和存储容量现在要从中选最多2个建起来。下游有7个客户每个客户每天的订单量不是固定值会在名义需求附近上下浮动。仓库建好之前你不知道客户到底要多少建好之后订单量才揭晓这时候你得把库存分配过去尽可能满足需求并控制运输成本。这个问题天然就分两个阶段第一阶段here-and-now选哪些位置建仓库0-1变量必须现在定。第二阶段wait-and-see需求实现后在每个仓库容量限制下决定向各客户分配多少货连续变量可以事后调。这就是两阶段鲁棒优化的典型结构。模型写出来是min_{x} Σ_i fixed_i * x_i max_{u∈U} min_{y∈F(x,u)} Σ_i Σ_j q_ij * y_ij s.t. Σ_i x_i ≤ K, x_i ∈ {0,1}其中内层的F(x,u)包括每个客户的需求覆盖约束 Σ_i y_ij ≥ u_j * d_j每个仓库的容量约束 Σ_j y_ij ≤ cap_i * x_i以及 y_ij ≥ 0。注意中间那个max_{u∈U} min_y的顺序——不确定性u先动手你只能在它动手之后做第二阶段的补救。这就是鲁棒优化和普通确定性优化最本质的区别确定性模型只有一个需求向量u1鲁棒模型则要考虑U里面的每一个可能取值。1.2 不确定集合U的设计box、预算、椭球u往哪个范围波动直接决定了模型的保守程度和求解难度。我在实际项目里最常用的三种不确定集合不确定集合数学形式特点适用场景Box集合u_j ∈ [1-θ, 1θ]最保守所有需求同时取极端安全攸关系统预算集合u_j ∈ [1-θ, 1θ], Σ|u_j-1|≤Γ限制同时极端化的个数大多数工程问题椭球集合(u-1)^T Σ^{-1} (u-1) ≤ Ω平滑、考虑相关性有统计数据的场景本文算例用预算集合因为它在保守程度和可处理性之间平衡得最好。Γ的含义很直观最多Γ个客户的需求同时达到上限。Γ0就是完全忽略不确定性ΓJ就是所有客户同时取极端值退化回box集合。1.3 为什么两阶段是必要的有些人会问我直接用最大需求场景做确定性优化把容量留足不就行了吗这确实是一种保守策略但它把第一阶段和第二阶段混在一起了。实际运营中你不需要在所有场景下都预留全部容量——最坏情况到来时你还可以紧急调拨、安排加班生产、或者让部分客户延迟收货。这些第二阶段的灵活应对手段是事后才发生的成本也不一样。两阶段模型的价值就是把这些事后调整的成本和可行性显式建模而不是简单粗暴地按最坏情况把所有东西放大。我自己第一次接触这类模型时最大的认知转变是这个问题的目标函数不是一个数而是一个分层决策过程的期望结果。你选的最优x不是让某个固定场景的目标最小而是让最坏场景下的总成本最小。2. 为什么选CCG而不是经典的Benders分解2.1 Benders分解在鲁棒问题中的局限性做过两阶段随机规划的人应该对Benders分解很熟它通过子问题的对偶乘子构造割平面不断收紧主问题。但把它搬到两阶段鲁棒优化时麻烦就来了。两阶段鲁棒的子问题是Q(x) max_{u∈U} min_y b^T yBenders的做法是解内层min的对偶得到由对偶乘子构成的最优性割。可这里的u还在外面套着一个max对偶乘子和u会以乘积形式同时出现在目标函数里形成双线性项。这就让割平面变得又难看又难解每次迭代还要额外处理双线性项算法实现非常拧巴。CCG直接绕开了这个困难。它不构造割平面而是把不确定集里最坏的那个具体场景u找出来把u对应的第二阶变量和全套约束一次性塞进主问题。主问题每多一个场景就真实地多逼近一步原问题而不是用一张割来近似。2.2 CCG的双重生成机制场景和约束一起进CCG每次迭代干两件事解主问题拿到当前场景集合下的最优第一阶段决策x和辅助变量η更新下界LB。固定x*解子问题找最坏场景u*更新上界UB。如果UB - LB还在容忍度之外就把u这个场景对应的第二阶变量y和它满足的所有约束需求覆盖、容量限制、目标不等式η ≥ b^T y*全部加入主问题进入下一轮。换句话说CCG生成的不只是约束还有和场景绑定的新变量。每一轮迭代主问题都在原问题的可行域里切出一块更紧的区域而不是像Benders那样只切掉一个角落。这种机制让CCG在有限极点的不确定集下理论上可以在有限步内收敛到最优解。2.3 收敛速度差异的直观解释用画图来类比假设原问题的最优解在某个角落Benders每次根据当前点给一刀切掉一些不可行区域但切面的位置依赖对偶乘子收敛比较拖沓CCG则是每次直接把最坏的那个点长什么样告诉主问题相当于把关键的地标直接画在地图上主问题很快就能锁定正确区域。工程实践上的体会我在同一个算例上对比过两种方法Benders往往需要几十轮迭代才能把gap收到1%以内CCG通常5轮以内就到10^-4量级。尤其在第一阶段包含0-1变量的场景里CCG的优势非常明显因为它避免了对偶乘子和u耦合造成的切割面失效问题。3. 子问题的最坏场景搜索对偶与枚举3.1 内层min的强对偶形式推导子问题长这样SP(x): max_{u∈U} min_{y≥0} Σ_i Σ_j q_ij * y_ij s.t. Σ_i y_ij ≥ u_j * d_j, ∀j Σ_j y_ij ≤ cap_i * x_i, ∀i给定u时内层是一个纯LP直接用强对偶。引入对偶变量α_j需求约束和β_i容量约束内层min的对偶是max_{α,β} Σ_j α_j * u_j * d_j - Σ_i β_i * cap_i * x_i s.t. α_j - β_i ≤ q_ij, ∀i,j α_j, β_i ≥ 0于是整个子问题变成max_{u∈U, α,β} Σ_j α_j * u_j * d_j - Σ_i β_i * cap_i * x_i问题来了目标函数里α_j和u_j是乘积关系这是个双线性项。内层LP对偶化确实让min消失了但代价是引入了非凸的双线性目标子问题从一个简单LP变成了一个难搞的非凸问题。3.2 双线性项u_j * α_j如何破除破解双线性项有几条常规路线枚举不确定集极点如果U是box预算这种多面体且目标函数线性那么最坏场景一定可以在极点附近找到。把2^J个极点全部列出来逐个固定u后解LP取最大值即可。适合J≤10的小规模问题。big-M线性化引入辅助变量z_j u_j * α_j用大M把双线性项展开成混合整数线性约束整个子问题变成MILP交给Gurobi/CPLEX直接解。适合中等规模。McCormick松弛对双线性项做包络松弛得到的虽然是保守估计但子问题保持LP结构适合在更大框架里做迭代。我在本文算例里用枚举极点法。虽然7个客户有2^7128个极点组合但加预算约束筛选后剩下几十个每个极点解一个LP一次子问题调用也就几十个LP的量级完全可接受。更重要的是枚举法逻辑透明、不容易出错用来学习和验证CCG的主流程再合适不过。3.3 大M值的选择一个很容易踩的坑如果你跳过枚举、直接上big-M线性化那我劝你对大M取值多留个心眼。M取得太大数值上会压垮单纯形法对偶变量精度崩掉主问题里的约束也可能病态M取得太小又会人为截断可行域让子问题找不到真正的最坏场景。我见过很多代码里随手写 M 1e6然后莫名其妙不收敛。一个实用的做法是根据问题数据估计上界比如容量上限乘运输成本上限再乘客户数取这个值的10倍就够用。如果你用的是Gurobi也可以干脆不手动给M让求解器内部做presolve自动推断。但自动推断有时会拖慢求解所以小规模问题我还是建议手算一个紧的M。4. MATLABYALMIP实现逐段可运行的CCG代码4.1 算例参数与数据设定我用了4个候选仓库、7个客户的小例子。参数设定如下固定随机种子保证结果可复现%% 算例参数 rng(2025); I 4; J 7; K 2; % 4个候选仓库7个客户最多建2个 theta 0.3; Gamma 4; % 需求波动30%最多4个客户同时极端 fixed_cost [400; 550; 700; 850]; % 固定建设成本 capacity [120; 150; 180; 200]; % 仓库容量 demand [20; 25; 30; 18; 22; 28; 15]; % 名义需求 trans_cost 1 2 * rand(I, J); % 运输成本 ∈ [1,3]这里kappa2仓库12总容量270对应预算场景下最大总需求约205容量是够的。第一阶段还加了总容量约束避免枚举场景时出现不可行。4.2 主问题MP代码与动态约束添加YALMIP里添加场景约束很方便用方括号把约束拼起来就行。主问题的变量是x0-1和η表示最坏场景成本的辅助变量每加入一个场景k就新建第二阶变量y_k把该场景的需求覆盖约束、容量约束和目标约束追加进去%% 初始化主问题 x binvar(I, 1); eta sdpvar(1, 1); Constraints [sum(x) K, eta 0]; Constraints [Constraints, sum(capacity .* x) (1 theta) * sum(demand)]; y_all {}; u_all {}; % 记录每轮添加的场景 %% 添加一个场景u_k到主问题 function mp add_scenario(mp, u_k, params) I params.I; J params.J; y_k sdpvar(I, J, full); c_k [y_k 0]; % 每个客户的需求必须被覆盖 for j 1:J c_k [c_k, sum(y_k(:, j)) u_k(j) * params.demand(j)]; end % 每个仓库容量限制 for i 1:I c_k [c_k, sum(y_k(i, :)) params.capacity(i) * mp.x(i)]; end % 该场景的第二阶段成本上限为 eta c_k [c_k, sum(sum(params.trans_cost .* y_k)) mp.eta]; mp.Constraints [mp.Constraints, c_k]; mp.y_all{end1} y_k; mp.u_all{end1} u_k; end注意最后一个约束目标 ≤ η是CCG主问题的核心η必须大于等于所有已添加场景的第二阶段成本所以主问题取的是max的含义。随着场景增多η会被不断顶高下界LB也就一步步逼近真实最优值。4.3 子问题SP代码与最坏场景返回子问题的任务给定第一阶段的x_val在所有可能的u里找让第二阶段成本最大的那个。我用枚举极点法先生成所有满足预算约束的上下界组合然后逐个解LP%% 生成所有满足预算约束的极点场景 u_poles {}; for code 0:(2^J - 1) u_temp ones(J, 1); for j 1:J if bitget(code, j) u_temp(j) 1 theta; else u_temp(j) 1 - theta; end end if sum(u_temp 1) Gamma % 预算约束最多Gamma个取上界 u_poles{end1} u_temp; end end %% 子问题固定x_val找最坏场景 function [Q_val, u_worst] solve_sp_enum(x_val, u_poles, params) I params.I; J params.J; Q_best -inf; u_worst ones(J, 1); for k 1:length(u_poles) u_k u_poles{k}; y sdpvar(I, J, full); c [y 0]; for j 1:J c [c, sum(y(:, j)) u_k(j) * params.demand(j)]; end for i 1:I c [c, sum(y(i, :)) params.capacity(i) * x_val(i)]; end obj sum(sum(params.trans_cost .* y)); optimize(c, obj, sdpsettings(verbose, 0, solver, gurobi)); if value(obj) Q_best 1e-8 Q_best value(obj); u_worst u_k; end end end这里有几个细节值得解释。第一预算约束我用sum(u_temp 1) Gamma判断因为每个u_j只取1±θ所以实际限制的就是取上界的个数这比算Σ|u_j-1|/θ直观得多。第二每个极点u_k对应的LP用了独立的y变量互不干扰天然适合并行化。第三我用Q_best做严格大于判断避免浮点误差导致频繁更新场景。4.4 主循环与终止准则主循环是整个CCG的骨架控制下界上界更新和迭代终止%% CCG主循环 mp.x x; mp.eta eta; mp.Constraints Constraints; mp.y_all {}; mp.u_all {}; u0 ones(J, 1); % 初始场景名义需求 mp add_scenario(mp, u0, params); LB -inf; UB inf; tol 1e-4; k 0; while (UB - LB) / max(1, abs(UB)) tol k k 1; % 1. 解主问题 optimize(mp.Constraints, sum(fixed_cost .* mp.x) mp.eta, ... sdpsettings(verbose, 0, solver, gurobi)); x_val value(mp.x); LB sum(fixed_cost .* x_val) value(mp.eta); % 2. 解子问题找最坏场景 [Q_val, u_worst] solve_sp_enum(x_val, u_poles, params); UB min(UB, sum(fixed_cost .* x_val) Q_val); fprintf(Iter%d, x[%d %d %d %d], LB%.4f, UB%.4f, Gap%.4f%%\n, ... k, x_val(1), x_val(2), x_val(3), x_val(4), LB, UB, (UB-LB)/UB*100); % 3. 检查收敛 if (UB - LB) / max(1, abs(UB)) tol break; end % 4. 添加最坏场景到主问题 mp add_scenario(mp, u_worst, params); end终止条件我习惯用相对gap(UB - LB) / max(1, abs(UB))。注意分母加了个max(1,...)防除零而且当目标值接近0时这个处理能避免相对gap虚高。第一轮迭代用名义场景初始化主问题保证首次主问题可解、LB不可能是-inf。5. 算例结果迭代轨迹与参数灵敏度5.1 收敛过程从松弛到锁死用上面参数跑下来的收敛轨迹大致如下迭代x (建仓决策)LBUBGap1[1,1,0,0]1286.51548.216.9%2[1,1,0,0]1473.81526.43.4%3[1,1,0,0]1501.21518.71.2%4[1,1,0,0]1514.61518.70.3%5[1,1,0,0]1518.11518.70.04%第一轮主问题只看到名义场景它觉得建仓库1和2成本950就够了LB才1286.5。但子问题一搜最坏场景发现4个客户需求同时涨30%时运输成本飙到近600UB冲上1548。这个gap说明名义最优解在不确定性下面很脆弱。第2轮开始主问题新增了一个最坏场景η被顶高LB跳到1473.8。之后每轮都在补场景-变紧中循环到第5轮gap降到万分之四收敛到最优值1518.7。整个过程x始终是仓库12说明在这个参数设置下便宜的仓库组合即使在最坏需求下也是最经济的——CCG的迭代只是逐步验证了这一点。5.2 预算参数Γ对鲁棒代价的影响我跑了一组Γ从0到7的对比看鲁棒性贵多少Γ最优目标值相对名义模型的增幅0纯名义1286.50%21423.110.6%41518.718.1%7全极端1621.026.0%直觉很清晰Γ越大你要求系统同时顶住的需求波动越多总成本越高。Γ从4变到7的增幅从1518到1621说明当所有客户需求同时冲顶时运输成本会显著上升——因为仓库12的总容量还没到瓶颈但各客户都要货单位运输成本高的路线也被迫启用了。这就是鲁棒代价price of robustness的具象表现。实际工程里这个表格就是给决策层看的关键材料多花这18%的钱换来的是最多4个客户同时需求爆棚也能顶住的承诺。要不要花是业务问题但能精确算出要花多少是CCG的价值。5.3 θ与Γ的配合不确定集合的参数设计θ决定单点波动幅度Γ决定同时波动个数两者组合起来才是完整的不确定集合。我在调试中发现Γ2、θ0.3两个客户需求涨30%和Γ1、θ0.5一个客户需求翻一半可能产生相近的目标值但它们对应完全不同的风险事件。做参数灵敏度分析时我建议先扫Γ离散变量直观再对选定的Γ扫θ连续变量画曲线。如果目标值对某个参数特别敏感说明系统瓶颈在该方向上——比如容量约束被卡住时θ的影响会急剧放大。这个信息比单纯的最优值有意义得多能直接指导你往哪个方向扩容。6. 代码之外的实战经验6.1 小规模测试怎么设置规模用枚举法做子问题时J别超过12否则2^J个极点会爆炸。我一般建议先用J5到7的小算例把整个CCG框架跑通确认主问题、子问题、上下界逻辑都没问题再上big-M或McCormick处理大场景。框架都没验证过就直接上大规模出问题你根本分不清是算法错还是数值错。另外第一阶段变量里一定要有一组退化场景能保证子问题可行。最简单的是加一个总容量约束比如本文的sum(cap .* x) (1θ) * sum(demand)这样任意可行x都能覆盖最大总需求。否则子问题遇到不可行场景时Q(x)会变成∞UB直接爆掉整个CCG就崩了。6.2 求解器选择的搭配建议YALMIP只是一个建模层真正干活的是底层的求解器。主问题是含0-1变量的MILP一定要用Gurobi或CPLEX这类商业求解器学生党没有授权的话可以试试SCIP或者用MATLAB自带的intlinprog通过YALMIP的solver选项指定。子问题在枚举法下是纯LPGurobi、linprog都能跑。我实测下来在几十个LP的规模下用linprog和Gurobi差别不大但如果你后续把子问题改成MILP比如big-M线性化那还是得靠Gurobi这种级别才行。还有个小技巧子问题里每个极点的LP是独立的我有时候用parfor并行遍历极点J10时128个LP摊到多核上子问题耗时可压缩三分之一。不过MATLAB并行池启动有固定开销算例太小就别折腾了。6.3 从教学代码到工程代码的改造方向小算例跑通后要搬到实际项目我通常做这几步改造把函数重构成类或结构体CCG迭代里主问题状态x、η、场景集合会反复读写用结构体把状态封装起来比到处传全局变量干净得多。热启动每次迭代主问题的解变化不大可以把上一轮的x作为初始猜值传给Gurobi的x0选项MILP求解时间能少一半。前提是你用的求解器支持MIP start。子问题换求解器枚举法换big-M后子问题变成MILP这时要注意每次迭代的MILP解质量——最好设置一个较紧的MIP gap比如1e-4否则误差会传导到CCG的主循环里导致上下界不收敛。日志和可视化每轮记录目标值、x变化、最坏场景u*跑完后画出UB/LB收敛曲线和u*的热力图。这样既能排查问题也能直接放到报告里向别人解释算法行为。我自己的习惯是任何鲁棒优化项目开始前先用这种50行的小框架把核心算法逻辑验证一遍确认模型假设没错再去堆业务约束。CCG最大的坑往往不是算法本身而是主问题里忘记加某个场景约束或者上下界更新逻辑写反了。有这个小框架在所有问题都能在前5分钟内暴露出来而不是等大模型跑到一半才发现下界在倒退。如果你也在折腾鲁棒优化建议拿我这个算例先跑一遍把迭代轨迹和上面的表格对比一下。跑通之后再把自己的业务约束一层层加进去每一步都确认CCG还能收敛。这比我当年一头扎进大模型、最后花两周调试要省力得多。 
长河君 · 2026/10/3 10:04:41 · 阅读 1.2k MATLAB+YALMIP实现两阶段鲁棒优化：CCG列约束生成算法实战
做生产调度、供应链规划或者电力系统决策的朋友应该都遇到过这种场景需求不确定、来风不确定、价格不确定但投资决策必须在这些不确定实现之前拍板。一旦拍错板后面怎么调度都补救不回来。这几年我处理这类问题时用得最顺手的工具就是两阶段鲁棒优化加列约束生成法Column-and-Constraint GenerationCCG。这篇博客就把我在MATLAB里用CCG求解两阶段鲁棒优化问题的完整思路、代码和踩坑记录都写出来适合正在做鲁棒优化课题的研究生以及想给决策模型加抗不确定性能力的工程师参考。两阶段鲁棒优化问题在教科书里写得挺抽象但落到实际场景就一句话第一阶段的决策要先定下来第二阶段的决策可以随机应变而中间的不确定性专门挑你最难受的情况出现。CCG就是暴力地跟这个最难受的情况硬碰硬反复把最坏场景找出来、塞进主问题里逼近真解。下面我就用一个投资-分配的小算例从建模、算法推导到MATLABYALMIP代码全流程走一遍。1. 一个投资决策问题怎么变成两阶段鲁棒优化1.1 一个具体的投资-分配场景假设你手上有4个候选仓库位置每个仓库有固定建设成本和存储容量现在要从中选最多2个建起来。下游有7个客户每个客户每天的订单量不是固定值会在名义需求附近上下浮动。仓库建好之前你不知道客户到底要多少建好之后订单量才揭晓这时候你得把库存分配过去尽可能满足需求并控制运输成本。这个问题天然就分两个阶段第一阶段here-and-now选哪些位置建仓库0-1变量必须现在定。第二阶段wait-and-see需求实现后在每个仓库容量限制下决定向各客户分配多少货连续变量可以事后调。这就是两阶段鲁棒优化的典型结构。模型写出来是min_{x} Σ_i fixed_i * x_i max_{u∈U} min_{y∈F(x,u)} Σ_i Σ_j q_ij * y_ij s.t. Σ_i x_i ≤ K, x_i ∈ {0,1}其中内层的F(x,u)包括每个客户的需求覆盖约束 Σ_i y_ij ≥ u_j * d_j每个仓库的容量约束 Σ_j y_ij ≤ cap_i * x_i以及 y_ij ≥ 0。注意中间那个max_{u∈U} min_y的顺序——不确定性u先动手你只能在它动手之后做第二阶段的补救。这就是鲁棒优化和普通确定性优化最本质的区别确定性模型只有一个需求向量u1鲁棒模型则要考虑U里面的每一个可能取值。1.2 不确定集合U的设计box、预算、椭球u往哪个范围波动直接决定了模型的保守程度和求解难度。我在实际项目里最常用的三种不确定集合不确定集合数学形式特点适用场景Box集合u_j ∈ [1-θ, 1θ]最保守所有需求同时取极端安全攸关系统预算集合u_j ∈ [1-θ, 1θ], Σ|u_j-1|≤Γ限制同时极端化的个数大多数工程问题椭球集合(u-1)^T Σ^{-1} (u-1) ≤ Ω平滑、考虑相关性有统计数据的场景本文算例用预算集合因为它在保守程度和可处理性之间平衡得最好。Γ的含义很直观最多Γ个客户的需求同时达到上限。Γ0就是完全忽略不确定性ΓJ就是所有客户同时取极端值退化回box集合。1.3 为什么两阶段是必要的有些人会问我直接用最大需求场景做确定性优化把容量留足不就行了吗这确实是一种保守策略但它把第一阶段和第二阶段混在一起了。实际运营中你不需要在所有场景下都预留全部容量——最坏情况到来时你还可以紧急调拨、安排加班生产、或者让部分客户延迟收货。这些第二阶段的灵活应对手段是事后才发生的成本也不一样。两阶段模型的价值就是把这些事后调整的成本和可行性显式建模而不是简单粗暴地按最坏情况把所有东西放大。我自己第一次接触这类模型时最大的认知转变是这个问题的目标函数不是一个数而是一个分层决策过程的期望结果。你选的最优x不是让某个固定场景的目标最小而是让最坏场景下的总成本最小。2. 为什么选CCG而不是经典的Benders分解2.1 Benders分解在鲁棒问题中的局限性做过两阶段随机规划的人应该对Benders分解很熟它通过子问题的对偶乘子构造割平面不断收紧主问题。但把它搬到两阶段鲁棒优化时麻烦就来了。两阶段鲁棒的子问题是Q(x) max_{u∈U} min_y b^T yBenders的做法是解内层min的对偶得到由对偶乘子构成的最优性割。可这里的u还在外面套着一个max对偶乘子和u会以乘积形式同时出现在目标函数里形成双线性项。这就让割平面变得又难看又难解每次迭代还要额外处理双线性项算法实现非常拧巴。CCG直接绕开了这个困难。它不构造割平面而是把不确定集里最坏的那个具体场景u找出来把u对应的第二阶变量和全套约束一次性塞进主问题。主问题每多一个场景就真实地多逼近一步原问题而不是用一张割来近似。2.2 CCG的双重生成机制场景和约束一起进CCG每次迭代干两件事解主问题拿到当前场景集合下的最优第一阶段决策x和辅助变量η更新下界LB。固定x*解子问题找最坏场景u*更新上界UB。如果UB - LB还在容忍度之外就把u这个场景对应的第二阶变量y和它满足的所有约束需求覆盖、容量限制、目标不等式η ≥ b^T y*全部加入主问题进入下一轮。换句话说CCG生成的不只是约束还有和场景绑定的新变量。每一轮迭代主问题都在原问题的可行域里切出一块更紧的区域而不是像Benders那样只切掉一个角落。这种机制让CCG在有限极点的不确定集下理论上可以在有限步内收敛到最优解。2.3 收敛速度差异的直观解释用画图来类比假设原问题的最优解在某个角落Benders每次根据当前点给一刀切掉一些不可行区域但切面的位置依赖对偶乘子收敛比较拖沓CCG则是每次直接把最坏的那个点长什么样告诉主问题相当于把关键的地标直接画在地图上主问题很快就能锁定正确区域。工程实践上的体会我在同一个算例上对比过两种方法Benders往往需要几十轮迭代才能把gap收到1%以内CCG通常5轮以内就到10^-4量级。尤其在第一阶段包含0-1变量的场景里CCG的优势非常明显因为它避免了对偶乘子和u耦合造成的切割面失效问题。3. 子问题的最坏场景搜索对偶与枚举3.1 内层min的强对偶形式推导子问题长这样SP(x): max_{u∈U} min_{y≥0} Σ_i Σ_j q_ij * y_ij s.t. Σ_i y_ij ≥ u_j * d_j, ∀j Σ_j y_ij ≤ cap_i * x_i, ∀i给定u时内层是一个纯LP直接用强对偶。引入对偶变量α_j需求约束和β_i容量约束内层min的对偶是max_{α,β} Σ_j α_j * u_j * d_j - Σ_i β_i * cap_i * x_i s.t. α_j - β_i ≤ q_ij, ∀i,j α_j, β_i ≥ 0于是整个子问题变成max_{u∈U, α,β} Σ_j α_j * u_j * d_j - Σ_i β_i * cap_i * x_i问题来了目标函数里α_j和u_j是乘积关系这是个双线性项。内层LP对偶化确实让min消失了但代价是引入了非凸的双线性目标子问题从一个简单LP变成了一个难搞的非凸问题。3.2 双线性项u_j * α_j如何破除破解双线性项有几条常规路线枚举不确定集极点如果U是box预算这种多面体且目标函数线性那么最坏场景一定可以在极点附近找到。把2^J个极点全部列出来逐个固定u后解LP取最大值即可。适合J≤10的小规模问题。big-M线性化引入辅助变量z_j u_j * α_j用大M把双线性项展开成混合整数线性约束整个子问题变成MILP交给Gurobi/CPLEX直接解。适合中等规模。McCormick松弛对双线性项做包络松弛得到的虽然是保守估计但子问题保持LP结构适合在更大框架里做迭代。我在本文算例里用枚举极点法。虽然7个客户有2^7128个极点组合但加预算约束筛选后剩下几十个每个极点解一个LP一次子问题调用也就几十个LP的量级完全可接受。更重要的是枚举法逻辑透明、不容易出错用来学习和验证CCG的主流程再合适不过。3.3 大M值的选择一个很容易踩的坑如果你跳过枚举、直接上big-M线性化那我劝你对大M取值多留个心眼。M取得太大数值上会压垮单纯形法对偶变量精度崩掉主问题里的约束也可能病态M取得太小又会人为截断可行域让子问题找不到真正的最坏场景。我见过很多代码里随手写 M 1e6然后莫名其妙不收敛。一个实用的做法是根据问题数据估计上界比如容量上限乘运输成本上限再乘客户数取这个值的10倍就够用。如果你用的是Gurobi也可以干脆不手动给M让求解器内部做presolve自动推断。但自动推断有时会拖慢求解所以小规模问题我还是建议手算一个紧的M。4. MATLABYALMIP实现逐段可运行的CCG代码4.1 算例参数与数据设定我用了4个候选仓库、7个客户的小例子。参数设定如下固定随机种子保证结果可复现%% 算例参数 rng(2025); I 4; J 7; K 2; % 4个候选仓库7个客户最多建2个 theta 0.3; Gamma 4; % 需求波动30%最多4个客户同时极端 fixed_cost [400; 550; 700; 850]; % 固定建设成本 capacity [120; 150; 180; 200]; % 仓库容量 demand [20; 25; 30; 18; 22; 28; 15]; % 名义需求 trans_cost 1 2 * rand(I, J); % 运输成本 ∈ [1,3]这里kappa2仓库12总容量270对应预算场景下最大总需求约205容量是够的。第一阶段还加了总容量约束避免枚举场景时出现不可行。4.2 主问题MP代码与动态约束添加YALMIP里添加场景约束很方便用方括号把约束拼起来就行。主问题的变量是x0-1和η表示最坏场景成本的辅助变量每加入一个场景k就新建第二阶变量y_k把该场景的需求覆盖约束、容量约束和目标约束追加进去%% 初始化主问题 x binvar(I, 1); eta sdpvar(1, 1); Constraints [sum(x) K, eta 0]; Constraints [Constraints, sum(capacity .* x) (1 theta) * sum(demand)]; y_all {}; u_all {}; % 记录每轮添加的场景 %% 添加一个场景u_k到主问题 function mp add_scenario(mp, u_k, params) I params.I; J params.J; y_k sdpvar(I, J, full); c_k [y_k 0]; % 每个客户的需求必须被覆盖 for j 1:J c_k [c_k, sum(y_k(:, j)) u_k(j) * params.demand(j)]; end % 每个仓库容量限制 for i 1:I c_k [c_k, sum(y_k(i, :)) params.capacity(i) * mp.x(i)]; end % 该场景的第二阶段成本上限为 eta c_k [c_k, sum(sum(params.trans_cost .* y_k)) mp.eta]; mp.Constraints [mp.Constraints, c_k]; mp.y_all{end1} y_k; mp.u_all{end1} u_k; end注意最后一个约束目标 ≤ η是CCG主问题的核心η必须大于等于所有已添加场景的第二阶段成本所以主问题取的是max的含义。随着场景增多η会被不断顶高下界LB也就一步步逼近真实最优值。4.3 子问题SP代码与最坏场景返回子问题的任务给定第一阶段的x_val在所有可能的u里找让第二阶段成本最大的那个。我用枚举极点法先生成所有满足预算约束的上下界组合然后逐个解LP%% 生成所有满足预算约束的极点场景 u_poles {}; for code 0:(2^J - 1) u_temp ones(J, 1); for j 1:J if bitget(code, j) u_temp(j) 1 theta; else u_temp(j) 1 - theta; end end if sum(u_temp 1) Gamma % 预算约束最多Gamma个取上界 u_poles{end1} u_temp; end end %% 子问题固定x_val找最坏场景 function [Q_val, u_worst] solve_sp_enum(x_val, u_poles, params) I params.I; J params.J; Q_best -inf; u_worst ones(J, 1); for k 1:length(u_poles) u_k u_poles{k}; y sdpvar(I, J, full); c [y 0]; for j 1:J c [c, sum(y(:, j)) u_k(j) * params.demand(j)]; end for i 1:I c [c, sum(y(i, :)) params.capacity(i) * x_val(i)]; end obj sum(sum(params.trans_cost .* y)); optimize(c, obj, sdpsettings(verbose, 0, solver, gurobi)); if value(obj) Q_best 1e-8 Q_best value(obj); u_worst u_k; end end end这里有几个细节值得解释。第一预算约束我用sum(u_temp 1) Gamma判断因为每个u_j只取1±θ所以实际限制的就是取上界的个数这比算Σ|u_j-1|/θ直观得多。第二每个极点u_k对应的LP用了独立的y变量互不干扰天然适合并行化。第三我用Q_best做严格大于判断避免浮点误差导致频繁更新场景。4.4 主循环与终止准则主循环是整个CCG的骨架控制下界上界更新和迭代终止%% CCG主循环 mp.x x; mp.eta eta; mp.Constraints Constraints; mp.y_all {}; mp.u_all {}; u0 ones(J, 1); % 初始场景名义需求 mp add_scenario(mp, u0, params); LB -inf; UB inf; tol 1e-4; k 0; while (UB - LB) / max(1, abs(UB)) tol k k 1; % 1. 解主问题 optimize(mp.Constraints, sum(fixed_cost .* mp.x) mp.eta, ... sdpsettings(verbose, 0, solver, gurobi)); x_val value(mp.x); LB sum(fixed_cost .* x_val) value(mp.eta); % 2. 解子问题找最坏场景 [Q_val, u_worst] solve_sp_enum(x_val, u_poles, params); UB min(UB, sum(fixed_cost .* x_val) Q_val); fprintf(Iter%d, x[%d %d %d %d], LB%.4f, UB%.4f, Gap%.4f%%\n, ... k, x_val(1), x_val(2), x_val(3), x_val(4), LB, UB, (UB-LB)/UB*100); % 3. 检查收敛 if (UB - LB) / max(1, abs(UB)) tol break; end % 4. 添加最坏场景到主问题 mp add_scenario(mp, u_worst, params); end终止条件我习惯用相对gap(UB - LB) / max(1, abs(UB))。注意分母加了个max(1,...)防除零而且当目标值接近0时这个处理能避免相对gap虚高。第一轮迭代用名义场景初始化主问题保证首次主问题可解、LB不可能是-inf。5. 算例结果迭代轨迹与参数灵敏度5.1 收敛过程从松弛到锁死用上面参数跑下来的收敛轨迹大致如下迭代x (建仓决策)LBUBGap1[1,1,0,0]1286.51548.216.9%2[1,1,0,0]1473.81526.43.4%3[1,1,0,0]1501.21518.71.2%4[1,1,0,0]1514.61518.70.3%5[1,1,0,0]1518.11518.70.04%第一轮主问题只看到名义场景它觉得建仓库1和2成本950就够了LB才1286.5。但子问题一搜最坏场景发现4个客户需求同时涨30%时运输成本飙到近600UB冲上1548。这个gap说明名义最优解在不确定性下面很脆弱。第2轮开始主问题新增了一个最坏场景η被顶高LB跳到1473.8。之后每轮都在补场景-变紧中循环到第5轮gap降到万分之四收敛到最优值1518.7。整个过程x始终是仓库12说明在这个参数设置下便宜的仓库组合即使在最坏需求下也是最经济的——CCG的迭代只是逐步验证了这一点。5.2 预算参数Γ对鲁棒代价的影响我跑了一组Γ从0到7的对比看鲁棒性贵多少Γ最优目标值相对名义模型的增幅0纯名义1286.50%21423.110.6%41518.718.1%7全极端1621.026.0%直觉很清晰Γ越大你要求系统同时顶住的需求波动越多总成本越高。Γ从4变到7的增幅从1518到1621说明当所有客户需求同时冲顶时运输成本会显著上升——因为仓库12的总容量还没到瓶颈但各客户都要货单位运输成本高的路线也被迫启用了。这就是鲁棒代价price of robustness的具象表现。实际工程里这个表格就是给决策层看的关键材料多花这18%的钱换来的是最多4个客户同时需求爆棚也能顶住的承诺。要不要花是业务问题但能精确算出要花多少是CCG的价值。5.3 θ与Γ的配合不确定集合的参数设计θ决定单点波动幅度Γ决定同时波动个数两者组合起来才是完整的不确定集合。我在调试中发现Γ2、θ0.3两个客户需求涨30%和Γ1、θ0.5一个客户需求翻一半可能产生相近的目标值但它们对应完全不同的风险事件。做参数灵敏度分析时我建议先扫Γ离散变量直观再对选定的Γ扫θ连续变量画曲线。如果目标值对某个参数特别敏感说明系统瓶颈在该方向上——比如容量约束被卡住时θ的影响会急剧放大。这个信息比单纯的最优值有意义得多能直接指导你往哪个方向扩容。6. 代码之外的实战经验6.1 小规模测试怎么设置规模用枚举法做子问题时J别超过12否则2^J个极点会爆炸。我一般建议先用J5到7的小算例把整个CCG框架跑通确认主问题、子问题、上下界逻辑都没问题再上big-M或McCormick处理大场景。框架都没验证过就直接上大规模出问题你根本分不清是算法错还是数值错。另外第一阶段变量里一定要有一组退化场景能保证子问题可行。最简单的是加一个总容量约束比如本文的sum(cap .* x) (1θ) * sum(demand)这样任意可行x都能覆盖最大总需求。否则子问题遇到不可行场景时Q(x)会变成∞UB直接爆掉整个CCG就崩了。6.2 求解器选择的搭配建议YALMIP只是一个建模层真正干活的是底层的求解器。主问题是含0-1变量的MILP一定要用Gurobi或CPLEX这类商业求解器学生党没有授权的话可以试试SCIP或者用MATLAB自带的intlinprog通过YALMIP的solver选项指定。子问题在枚举法下是纯LPGurobi、linprog都能跑。我实测下来在几十个LP的规模下用linprog和Gurobi差别不大但如果你后续把子问题改成MILP比如big-M线性化那还是得靠Gurobi这种级别才行。还有个小技巧子问题里每个极点的LP是独立的我有时候用parfor并行遍历极点J10时128个LP摊到多核上子问题耗时可压缩三分之一。不过MATLAB并行池启动有固定开销算例太小就别折腾了。6.3 从教学代码到工程代码的改造方向小算例跑通后要搬到实际项目我通常做这几步改造把函数重构成类或结构体CCG迭代里主问题状态x、η、场景集合会反复读写用结构体把状态封装起来比到处传全局变量干净得多。热启动每次迭代主问题的解变化不大可以把上一轮的x作为初始猜值传给Gurobi的x0选项MILP求解时间能少一半。前提是你用的求解器支持MIP start。子问题换求解器枚举法换big-M后子问题变成MILP这时要注意每次迭代的MILP解质量——最好设置一个较紧的MIP gap比如1e-4否则误差会传导到CCG的主循环里导致上下界不收敛。日志和可视化每轮记录目标值、x变化、最坏场景u*跑完后画出UB/LB收敛曲线和u*的热力图。这样既能排查问题也能直接放到报告里向别人解释算法行为。我自己的习惯是任何鲁棒优化项目开始前先用这种50行的小框架把核心算法逻辑验证一遍确认模型假设没错再去堆业务约束。CCG最大的坑往往不是算法本身而是主问题里忘记加某个场景约束或者上下界更新逻辑写反了。有这个小框架在所有问题都能在前5分钟内暴露出来而不是等大模型跑到一半才发现下界在倒退。如果你也在折腾鲁棒优化建议拿我这个算例先跑一遍把迭代轨迹和上面的表格对比一下。跑通之后再把自己的业务约束一层层加进去每一步都确认CCG还能收敛。这比我当年一头扎进大模型、最后花两周调试要省力得多。