企业官网建设流程全解析

搞两阶段鲁棒优化的人，大概率都绕不开CCG这三个字母。在我接触到的大多数实际决策模型——从生产排产到运输网络——只要不确定性出现在约束右侧，目标又带着“先决策、再应对最坏情况”的味道，最后基本都会落到同一个形式：min_x c^T x + max_{u∈U} min_y d^T y。CCG（Column-and-Constraint Generation，列约束生成）就是专门啃这类问题的算法。近期我重新用MATLAB和Python各实现了一遍两阶段鲁棒问题的CCG求解器，收获不少，也踩了一些文档里不会写的坑。这篇文章我不打算讲太多纯定理推导，而是把它怎么拆、怎么用两套语言写出来、以及我实际调代码时积累的经验一次讲清楚。适合刚把两阶段鲁棒模型建好、正卡在求解这一步的人。

1. 问题重塑：两阶段鲁棒优化到底在求什么

1.1 数学形式与信息结构

先说建模。标准两阶段鲁棒优化可以写成：

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

其中：

• x是第一阶段决策，也就是现在就要拍板的方案，对应英文里的 here-and-now decision；
• u是不确定参数，只要求落在一个集合U里，不假设概率分布；
• y是第二阶段决策，发生在不确定参数被观测到之后，对应 wait-and-see decision；
• F(x,u)是第二阶段可行域，通常由一组含x和u的线性不等式定义。

我习惯拿一个生产库存的例子来理解：第一阶段决定建几条产线、建在哪里，这是马上要花的钱；第二阶段面对某一周的实际需求u，再决定每个工厂生产多少、怎么调配。鲁棒优化的意思是，不要赌需求一定等于预测值，而是假设需求可能落在某个最坏区间内，我们求的是“在最坏需求实现下总成本最小的第一阶段方案”。

很多人第一次接触这个模型会问：为什么不直接 min-max 求一个方案？因为x必须在你看到u之前定下来，而y是可以根据u调整的。如果y也能等到看到u再定，那整个问题就不是两阶段，而是一个带情景的随机规划。这个“看到 u 之后再优化 y”的结构，正是两阶段鲁棒和单阶段鲁棒最大的差异。

1.2 为什么不能直接交给求解器

如果你的U只有几十个离散场景，那确实可以直接枚举，把每个场景对应的一组y变量和耦合约束全部写进主模型，变成一个大MIP，丢给求解器就行。实际情况是u常常是连续区间、带预算约束的不确定集合，比如：

u_j = d_j^0 + d_j^hat * z_j, z_j∈[0,1], Σz_j ≤ Γ

这里的z_j可以是连续的，场景数是无穷的，根本没法枚举。就算把z限定成0-1变量，组合数量也会爆炸式增长，例如50个变量里选20个，枚举量是天文数字。

所以需要一种分解算法：主问题只保留一小部分“已经发现的最坏场景”，子问题负责在给定第一阶段决策下找新的最坏场景，一边迭代一边逼近真实最优。CCG就是这样一套框架，也是我这两年实际项目中用得最多的求解方案。

2. CCG算法核心：主问题加列，子问题找场景

2.1 主问题：固定已知场景下的下界近似

假设现在已经有了K个不确定场景u_1, u_2, ..., u_K，主问题可以写成：

min c^T x + η

约束条件包括：

• x∈X
• 对每个场景k：η ≥ d^T y_k
• 对每个场景k：D y_k ≥ b - A x - B u_k
• 对每个场景k：y_k ≥ 0

这个主问题是一个标准优化问题，可以直接用MIP/LP求解器求解。为什么说它是下界？因为我们只考虑了U中已经发现的有限个场景，相当于把不确定集合缩小了，求出来的目标值一定不会大于真实最优值，所以记为LB。

代码上的关键点是：每轮迭代要往主问题里加入新场景u_{K+1}，同时新增一组第二阶段变量y_{K+1}，以及对应的η ≥ d^T y_{K+1}和耦合约束。这个动态“加列”的动作，就是CCG名字里 Column 的由来。相比Benders分解只加一条割约束，CCG加的是场景变量和对应约束，主问题会更庞大，但换来的收敛性通常更好。

2.2 子问题：寻找最坏场景

给定第一阶段决策x_0，子问题是：

Q(x_0) = max_{u∈U} min_{y∈F(x_0,u)} d^T y

这是一个 max-min 双层结构。内层关于y是线性规划，可以直接对偶。假设第二阶段问题满足强对偶条件，内层对偶可以写成：

max_{π} π^T(b - A x_0 - B u)

其中π是对偶变量，满足π^T D ≤ d^T, π ≥ 0。再对u做外层最大化，就变成：

max_{π,u} π^T(b - A x_0 - B u)

麻烦的地方是目标里出现π^T B u这样的双线性项，不再是线性规划。解决办法通常有三种：

1. 如果U的极点很少，直接枚举极值点，每个极值点解一个LP，教学和验证都很方便；
2. 如果u是0-1变量且带预算约束，用大M法把双线性项线性化，得到一个MILP，交给Gurobi这样的求解器；
3. 如果问题规模太大，可以考虑启发式生成场景，但会牺牲最优性保证。

CCG的适用前提是：第二阶段问题是连续LP，且不确定集合是多面体。最坏场景一定可以在U的某个极点取到。如果第二阶段有整数变量，内层对偶不再成立，处理起来会复杂很多，一般要改用KKT条件或者专门的分解策略。

2.3 迭代流程与收敛依据

整个CCG循环我一般写成如下伪逻辑：

1. 初始化场景集，通常先放入名义场景u_0；
2. 求解主问题，得到最优解x*，目标值记为LB；
3. 把x*带入子问题，求解得到最坏场景u*和子问题最优值Q(x*)；
4. 计算上界UB = min{UB, c^T x* + Q(x*)}；
5. 若(UB - LB) / abs(LB) ≤ ε，停止；
6. 否则把u*加入主问题场景集，回到第2步。

这里有个容易理解错的地方：LB是每次重新求解主问题得到的目标值，会随场景增加而上升；UB是通过当前x*构造的可行解目标值，是真实目标的一个上界。两个值不断逼近，gap逐渐缩小。如果某次迭代子问题返回的场景已经在主问题里出现过，说明当前x*已经能够应付所有已知最坏场景，继续迭代也不会有新的信息，可以直接收敛。

理论上，当U是有限极点集时，CCG在有限步内收敛；即使U是连续多面体，只要最坏场景在极点取到，算法也有收敛保证。这个性质比Benders分解更让人省心。

3. MATLAB实现要点与关键代码骨架

3.1 用YALMIP组织动态主问题

MATLAB环境下我一般用YALMIP工具箱写模型，求解器接Gurobi或者CPLEX。YALMIP的好处是变量声明和约束拼装非常自然，特别适合研究阶段快速改模型。

主问题动态加列的关键，是用cell数组保存每一轮新增的第二阶段变量y{k}。不要试图提前把所有变量都定义好，那样场景数不定时很难管理。我的习惯是每轮迭代前先把当前已知场景整理出来，然后重新构建整个主问题模型。

一段典型的MATLAB主循环骨架如下：

% shadows: 已发现的不确定场景，cell数组，每个元素是一列向量 shadows = {}; shadows{1} = nominal_u; % 名义场景 T = 1; LB = -1e6; UB = 1e6; eps = 1e-4; while (UB - LB) / abs(LB) > eps % 1. 重新构建主问题 x = binvar(nx,1); eta = sdpvar(1); cons = []; y = {}; for k = 1:T y{k} = sdpvar(ny,1); % 该场景对应的第二阶段变量 cons = [cons, eta >= d' * y{k}]; % 场景成本的epigraph约束 cons = [cons, D * y{k} >= b - A * x - B * shadows{k}]; % 耦合约束 cons = [cons, y{k} >= 0]; end obj = c' * x + eta; optimize(cons, obj, sdpsettings('verbose',0)); LB = value(obj); x0 = value(x); % 2. 求解子问题，返回最坏场景和对应Q值 [u_star, sp_val] = solve_subproblem(x0, params); % 3. 更新上界 UB = min(UB, c' * x0 + sp_val); % 4. 加入新场景 if norm(u_star - shadows{end}) > 1e-6 T = T + 1; shadows{T} = u_star; else break; end end

这段代码里有几个点需要特别注意。第一，主问题里必须有eta这个标量变量，它代表“当前已知场景中的最大第二阶段成本”，否则目标函数写不成c^T x + η的形式。第二，shadows里每一列要是一个列向量，和B * u的维度保持一致，否则YALMIP经常会报维度不匹配的错误。第三，每次循环都重新构建主问题，模型规模越来越大，这很正常，不要尝试去“增量更新”同一个YALMIP模型，容易把自己绕晕。

3.2 子问题：写一个可以验证的枚举版本

对于教学和中小规模的案例，我强烈建议先写一个枚举极值点的子问题版本，用来验证主问题逻辑。等确认整体流程没问题，再换大规模使用的线性化版本。

枚举版子问题的MATLAB代码可以这样写：

function [u_best, sp_obj] = solve_subproblem(x0, params) U = params.U; % 预先生成的极值点矩阵，每行一个极值点 sp_obj = -inf; u_best = []; for i = 1:size(U,1) ui = U(i,:)'; y = sdpvar(params.ny,1); cons = [y >= 0, params.D * y >= params.b - params.A * x0 - params.B * ui]; optimize(cons, params.d' * y, sdpsettings('verbose',0)); val = value(params.d' * y); if val > sp_obj sp_obj = val; u_best = ui; end end end

为什么先用枚举版？因为它能帮你直接定位问题到底出在主问题还是子问题。如果枚举版都跑不对，那线性化版本更不可能对。我在实际调试中，经常先用U的几个手算极值点测试子问题，确认sp_obj和手算结果一致，再进主循环。

3.3 MATLAB实现中的几个注意点

• YALMIP默认求解器选择有时不是最优，建议显式指定，例如sdpsettings('solver','gurobi')。
• 如果主问题是大规模MIP，尽量开求解器的MIP gap控制，比如设置相对MIP gap为1e-3，不要追求精确到1e-6，否则每轮主问题求解时间会拖垮整个CCG。
• value(x0)取出的是数值解，子问题里用的必须是这样一组具体数值，不能直接把sdpvar变量传给子问题。
• 枚举极值点时，如果U是连续集合，要先手动生成所有极点。区间预算集合的极点可以用0-1组合加预算约束生成，维度不高时完全可行。

4. Python实现要点与Gurobi对接

4.1 用Pyomo动态扩展主问题

Python环境下我常用Pyomo建模，求解器接Gurobi。Pyomo的ConstraintList非常适合CCG这种动态加约束的场景。场景变量存在字典里，第k个场景对应m.y[k]。

主问题构建函数可以这样写：

def build_master(known_u, data): m = ConcreteModel() nx, ny, nc = data['nx'], data['ny'], data['nc'] m.x = Var(range(nx), domain=Binary) m.eta = Var(domain=NonNegativeReals) m.y = {} m.Cuts = ConstraintList() m.Cuts.add(expr=sum(data['c'][i] * m.x[i] for i in range(nx)) + m.eta >= 0) for k, uk in enumerate(known_u): m.y[k] = Var(range(ny), domain=NonNegativeReals) m.Cuts.add(expr=m.eta >= sum(data['d'][j] * m.y[k][j] for j in range(ny))) for i in range(nc): expr = sum(data['D'][i][j] * m.y[k][j] for j in range(ny)) \ >= data['b'][i] \ - sum(data['A'][i][l] * m.x[l] for l in range(nx)) \ - sum(data['B'][i][j] * uk[j] for j in range(len(uk))) m.Cuts.add(expr=expr) m.obj = Objective( expr=sum(data['c'][i] * m.x[i] for i in range(nx)) + m.eta, sense=minimize ) return m

这里每个新场景都会在m.y里增加一个变量，并在m.Cuts里增加若干约束。Pyomo的一个好处是，约束表达式可以直接写Python的sum生成器，代码清晰度比MATLAB高不少。

主循环里，每轮求解完主问题后，把x0传给子问题函数，得到u_star，再检查是否已经存在，不存在就加入known_u列表，继续迭代。Pyomo每次从build_master重建模型，虽然会损失一点点构建时间，但换来了逻辑简单的可维护性。

4.2 子问题用Gurobi建模并对偶线性化

当U采用预算不确定集合且z是0-1变量时，子问题可以线性化成MILP用Gurobi直接求解。核心思路是先把内层min写成对偶形式，再用大M法处理双线性项。

具体来说，内层对偶后子问题形式为：

max π^T(b - A x0) - Σ_j w_j

其中：

• π是对偶变量，满足π^T D ≤ d^T，且π ≥ 0；
• w_j = z_j * (π^T B_j)，z_j是0-1不确定参数；
• 额外约束Σ_j z_j ≤ Γ。

w_j是连续变量和0-1变量的乘积，标准大M线性化需要测试四个约束。Gurobi建模子问题的代码如下：

import gurobipy as gp from gurobipy import GRB def solve_sp(x0, data): with gp.Model("sp") as sp: nc, nx, ny, nu = data['nc'], data['nx'], data['ny'], data['nu'] M = data['BigM'] pi = sp.addVars(nc, lb=0, ub=GRB.INFINITY, name="pi") z = sp.addVars(nu, vtype=GRB.BINARY, name="z") # 对偶可行约束 for j in range(ny): sp.addConstr( sum(pi[i] * data['D'][i][j] for i in range(nc)) <= data['d'][j] ) # 预算不确定集合 sp.addConstr(z.sum() <= data['Gamma']) # w_j = z_j * (pi^T B_j) w = {} for j in range(nu): v = sp.addVar(lb=-M, ub=M, name=f"v_{j}") w[j] = sp.addVar(lb=-M, ub=M, name=f"w_{j}") sp.addConstr(w[j] <= M * z[j]) sp.addConstr(w[j] >= -M * z[j]) sp.addConstr(w[j] <= v + M * (1 - z[j])) sp.addConstr(w[j] >= v - M * (1 - z[j])) sp.addConstr( v == sum(pi[i] * data['B'][i][j] for i in range(nc)) ) # 目标函数 obj = sum( pi[i] * (data['b'][i] - sum(data['A'][i][l] * x0[l] for l in range(nx))) for i in range(nc) ) - sum(w[j] for j in range(nu)) sp.setObjective(obj, sense=GRB.MAXIMIZE) sp.Params.MIPGap = 1e-4 sp.Params.TimeLimit = 300 sp.optimize() z_star = [round(z[j].X) for j in range(nu)] sp_obj = sp.ObjVal return z_star, sp_obj

BigM的取值是这里最大的隐患。我在调试时吃过几次亏：M取小了，线性化约束会错误地限制变量范围；M取大了，Gurobi求解MIP的数值会变差，甚至出现错误的整数解。一个稳妥的取法是先看对偶变量π的可行域边界，再结合B矩阵的元素上界估算|π^T B_j|的上限，然后取这个上限的3到5倍。

4.3 Python实现和MATLAB实现的差异

用两套语言写完同一个CCG后，最大体会是：MATLAB+YALMIP适合快速原型验证，Pyomo+Gurobi更适合做批量实验和工程化。YALMIP的约束拼接非常自由，写抽象模型速度快，但动态扩展变量时对坐标系和变量命名的要求比较高，错了很难定位。Pyomo用ConstraintList和字典管理场景变量，逻辑上更清晰，适合模型规模变大后维护。

另外，MATLAB如果不用YALMIP，直接写对偶模型会比较痛苦；Python生态里除了Pyomo，还可以直接用Gurobi的addVars和addConstrs构建主问题和子问题，代码会紧凑很多。我个人建议是：要写论文算例用MATLAB，要做工程应用和参数扫描用Python。

5. 调参、收敛判据与常见坑

5.1 收敛判据和gap设置

收敛判据不要用绝对gap，除非你的目标值本身只有几十。我一般用相对gap：

(UB - LB) / abs(LB) ≤ ε

默认ε取1e-4，但实际中如果主问题是MIP，求解器自身还要设置一个MIP gap，内外两层gap叠加可能会导致最终结果并不那么“精确”。我的习惯是：外层CCG取1e-4，内层求解器MIP gap取1e-3或1e-4之间，避免内层求解过度浪费时间。

需要特别提醒的是，如果第一阶段变量本身是整数变量，主问题每次都是一个MIP，求解时间会明显变长。实际项目中，我会给主问题设置一个合理的求解时间上限，比如300秒，然后用当前的可行解x*继续跑子问题。虽然这样可能丢掉严格下界，但工程上往往能接受。

5.2 双线性项与大M取值的坑

子问题里π^T B u这个双线性项，是所有实现细节里最容易被忽视的地方。很多初学CCG的人会在子问题里直接尝试用非线性求解器解π^T B u，结果又慢又不稳定。更常见的错误是把0-1变量z当成连续变量，然后强行用大M线性化，导致取到的“最坏场景”根本不是极值点。

大M线性化只适用于二值变量和连续变量的乘积。如果u是连续区间变量，比如z_j∈[0,1]，严格来说需要把连续变量做二进制展开后线性化，或者利用不确定集合极点性质转成枚举。不要在连续u上直接套用0-1的大M公式，那会得到错误结果。

还有预算约束的方向：很多人写Σz_j = Γ，这表示不确定参数必须正好偏离Γ项；但鲁棒优化里更常见的是Σz_j ≤ Γ，允许最坏情况“最多偏离Γ项”。这两个模型解出来不一样，写代码前要先明确建模需求。

5.3 求解器参数与数值稳定性

Gurobi求解MIP时，我通常会设以下参数：

• MIPFocus=1：偏向快速找到好可行解，因为CCG外层更关心UB质量；
• MIPGap=1e-4：控制内层精度；
• TimeLimit=300：防止某个超大MIP拖死整个流程；
• Threads=0：让求解器自动选线程数。

还有一点是关于场景去重的。子问题返回的u_star如果已经存在于主问题场景集中，应该直接停止迭代。但数值上两个场景可能只差1e-9，直接==判断不可靠。我是按norm(u_star - uk, inf) < 1e-6来判断是否重复，阈值根据问题量级调整。

5.4 第二阶段不可行怎么办

标准的CCG一般假设相对完全补偿，也就是对U里的任何u，第二阶段问题都有可行解y。如果实际模型做不到，比如产能不足导致需求无法满足，那子问题求解时会报不可行。这种情况需要引入可行性割：在子问题不可行时，添加一组辅助变量，找到能恢复可行性的最小调整方向，并把对应约束加入主问题。这在Benders分解里很常见，CCG也可以类似处理，但要额外小心。我的建议是先在模型层面加上惩罚变量或松弛变量，避免不可行问题打断迭代，等算法逻辑通了，再考虑要不要换成严格可行性割。

6. 缩微算例与迭代记录

6.1 测试模型与参数

为了验证代码，我常用一个非常简单的两阶段生产分配问题作为测试算例。第一阶段决定两个工厂是否启用，第二阶段面对不确定需求做生产分配。模型参数故意设成只有三个不确定极值点，这样可以手算验证。

不确定集合取名义场景加两个偏差场景：需求点1偏高，或需求点2偏高，预算约束保证两者不能同时取满。这个集合只有三个极点，枚举版子问题很方便。

主问题第一轮只用名义场景，CCG求解后得到第一阶段方案，带入子问题会返回一个最坏场景；第二轮把该场景加入主问题，再求解得到新的方案；第三轮如果子问题返回的场景已经在主问题里，算法结束。

6.2 迭代记录

一次典型的迭代数据如下：

第二轮的gap直接归零，是因为不确定集合一共只有三个极值点，而第一个名义场景在构建主问题时已经加入。第二轮加入一个偏差场景后，哪些场景该选已基本确定。这说明CCG在场景结构简单时收敛极快。如果换成大模型，比如50个需求点、预算Γ=20，极值点数量极大，CCG通常也只需要几十次迭代，而不会遍历所有组合，这正是它的价值所在。

6.3 什么时候该用CCG，什么时候该考虑其他方法

就我个人的实际体会，CCG最适合下面几类情况：

• 不确定集合是多面体，第二阶段是连续LP，最优场景在极值点取到；
• 第一阶段是整数决策，主问题可以接受MIP求解开销；
• 不确定变量数量多，但预算约束让极值场景“有规律”。

如果第二阶段本身含有整数变量，比如运维决策里涉及设备启停，内层对偶不再成立，CCG就不能直接套，需要结合KKT条件、列生成或启发式场景筛选来做。另外，如果主问题每轮都是一个几万变量的大MIP，CCG几十轮迭代会非常慢，这时候可能需要考虑Benders分解或者场景聚类的加速技巧。

我还想分享一个小经验：不管用MATLAB还是Python实现，第一次跑通整个CCG，我都会把热节点全部关闭，打印每一轮的主问题目标值、子问题返回场景编号、UB、LB和gap。等确认输出符合手算预期，再逐步放开求解器参数。这个方法帮我省了至少一周的调试时间。CCG看起来只是一个循环加两个模型，真正写顺手，里面还是有不少细节值得慢慢打磨。