企业官网建设流程全解析

简介：一套基于YALMIP的两阶段鲁棒优化MATLAB实现，同时涵盖CCG列与约束生成和Benders分解，面向需要求解不确定条件下决策问题的学术研究者与工程师。压缩包共5个文件，其中2个.m脚本为核心算法代码，3个PDF为方法对比与综述文献，整体大小仅1.44MB，便于快速下载与阅读，适合作为学习与二次开发的起点。目前已有2287人学习下载，代码复现自原论文，支持扩展改编，用户可结合YALMIP环境直接运行或根据自身模型调整约束与场景生成。CCG通过逐步添加切割平面逼近最坏情况下的最优解，Benders分解则将不确定性处理与主决策问题分离，两种方法互补性强，适合对比学习。资源内置两类算法的对比案例及参考论文，帮助读者理解迭代切割逻辑与实现细节，是一份兼顾理论参考与工程落地的实用学习资料。

1. 两阶段鲁棒优化不是黑匣子：C&CG 与 Benders 到底在算什么

接手过库存或产能规划的人大概都撞过这样的墙：需求历史数据只有上下界，没有像样的概率分布，随机优化根本无从建模。两阶段鲁棒优化就是为这种场景生的——第一阶段先做“现在必须定”的决策，第二阶段等不确定参数揭晓后再做“还能调整”的决策，目标是最坏情况下的总成本。C&CG（Column-and-Constraint Generation，列与约束生成）和 Benders 分解是求解这类问题的两条主流路线，两套代码其实共用同一个主问题-子问题骨架，只是切“割”的方式不同。这篇文章把模型、最小可复现代码、扩展改编路径和最容易翻车的细节一次讲透，适合准备复现原论文、或者想把手头确定性模型升级成鲁棒版本的人。

2. 先立模型：主问题—子问题结构与两条不同的割生成路径

2.1 两阶段鲁棒优化的标准写法与决策结构

两阶段鲁棒优化问题的标准形式可以写成：

[ \min_{x \in X} \left( c^T x + \max_{u \in U} \min_{y \in \mathcal{F}(x,u)} b^T y \right) ]

其中 (x) 是第一阶段决策变量（here-and-now），比如产能投资、仓库选址、采购合同量；(u) 是不确定参数，限定在不确定集 (U) 内；(y) 是第二阶段决策变量（wait-and-see），在看到 (u) 的实际取值后才能定，它的可行域 (\mathcal{F}(x,u)) 与第一阶段的 (x) 和第二阶段的 (u) 都相关。

之所以说“鲁棒”，核心在于不确定集而不是概率分布。最常见的两种不确定集是：

• 盒式不确定集：(U = {u \mid 0 \le u_j \le 1})，每个参数独立波动。
• 预算不确定集：(U = {u \mid 0 \le u_j \le 1,\ \sum_j u_j \le \Gamma})，用预算参数 (\Gamma) 控制“最多同时有多少个参数能取到坏值”，避免所有参数同时拉满导致结果过于保守。

两阶段鲁棒优化要同时处理两个困难：外层对 (u) 求最大（找最坏场景），内层对 (y) 求最小（在该场景下做最优调整）。直接求解这个 min-max-min 结构非常困难，所以工程上几乎都用“主问题-子问题”迭代：先猜一组场景，解主问题得到一个候选 (x)，再把 (x) 固定住，解子问题找出当前真正最坏的场景，把这个场景“切成约束”加回主问题，循环直到收敛。

2.2 C&CG 主问题—子问题迭代与有限收敛

C&CG 的英文全称是 Column-and-Constraint Generation，名字本身就说明了它的核心操作：每一轮迭代往主问题里加一“列”新变量和一组新约束。假设已经迭代了 (K) 轮，积累了 (K) 个最坏场景 (u_1^, \dots, u_K^)，主问题写成：

[ \min_{x, \eta, y_1, \dots, y_K} c^T x + \eta ]

[ \text{s.t.} \quad Ax \ge d,\quad x \in X ]

[ \eta \ge b^T y_k,\quad k = 1,\dots,K ]

[ W y_k \ge h - T x - M u_k^*,\quad y_k \in Y,\quad k = 1,\dots,K ]

这里 (\eta) 是第二阶段成本的代理变量。每加入一个场景，就新增一组第二阶段变量 (y_k) 和对应约束，这就是“列”的由来。子问题是固定第一阶段决策 (x^*) 后求解：

[ Q(x^) = \max_{u \in U} \min_{y \in \mathcal{F}(x^,u)} b^T y ]

迭代流程如下：

1. 初始化下界 (LB = -\infty)、上界 (UB = +\infty)，场景集合为空。
2. 求解主问题，得到 (x^)、(\eta^)，令 (LB = c^T x^* + \eta^*)。
3. 固定 (x^)，求解子问题，得到最坏场景 (u^) 和第二阶段最优值 (Q(x^*))。
4. 令 (UB = \min(UB, c^T x^* + Q(x^*)))。
5. 若 (UB - LB < \epsilon)，停止；否则把 (u^*) 加入场景集合，回到第 2 步。

Zeng 和 Zhao 在 2013 年发表的 C&CG 原论文里证明了：当第二阶段是线性规划时，这个算法在有限步内收敛到全局最优。为什么？因为不确定集 (U) 是多面体，线性规划的最优解一定可以在 (U) 的极点处取到，而极点数量有限，每次加入的场景不会重复。

2.3 Benders 与 C&CG：割的结构差异和选型理由

Benders 分解比 C&CG 早了大半个世纪，思路完全不同。C&CG 往主问题里加“原变量 + 原约束”，Benders 则是对子问题做对偶，往主问题里加一条只含 (\eta) 和 (x) 的标量割。

固定 (u) 后，内层 (\min_{y \in \mathcal{F}(x,u)} b^T y) 的 LP 对偶可以写成：

[ \max_{\alpha, \beta} \ \beta^T (h - M u) - \alpha^T x ]

[ \text{s.t.} \quad W^T \beta \le b,\quad T^T \alpha + W^T \beta \ge 0,\quad \alpha, \beta \ge 0 ]

对偶解 ((\alpha^, \beta^)) 就给出了一条 Benders 割：

[ \eta \ge (\beta^)^T (h - M u^) - (\alpha^*)^T x ]

两条路线的差别很具体：

• 主问题规模：C&CG 每轮加一组变量和约束，主问题越来越大；Benders 每轮只加一条不等式，变量数不变。
• 割的紧度：C&CG 的割对应一个完整的原场景，约束更紧，收敛更快；Benders 的割是对偶空间的线性近似，往往需要很多轮才能逼近真实曲面。
• 适用性：第二阶段含整数变量时，C&CG 的子问题用 MILP 求解器仍可直接跑，割结构不变；Benders 的纯 LP 对偶推导会失效，处理起来麻烦得多。

所以现在的两阶段鲁棒优化论文里，C&CG 几乎是默认算法，Benders 更多出现在随机优化和网络设计这类对偶结构特别漂亮的问题中。两套代码的骨架完全一样，差别只在“SP 返回什么、MP 加什么”，这为后续代码复用留了很大空间。

3. 用 Python + CVXPY 跑通 C&CG：最小可复现代码与关键参数

3.1 问题定义与数据：需求不确定下的产能投资算例

为了不绕弯子，用一个经典的产能投资问题做载体：I 个工厂，J 个客户，第一阶段决定每个工厂建多少产能，单位成本 c；第二阶段客户需求不确定，需要安排每个工厂给每个客户供多少货，单位生产成本 b，受工厂产能约束和客户需求约束。需求在基础值 d0 和波动幅度 delta 之间变化，不确定集先用盒式。

这是两阶段鲁棒优化里最常用的入门算例，结构简单但 C&CG 的每个要素都齐全。数据这样设定：3 个工厂、5 个客户，产能单位成本 0.9～1.2，生产成本随工厂和客户索引略有差异，需求基础值统一 100，波动 25～40。

import numpy as np import cvxpy as cp from itertools import product # 算例规模与数据 I, J = 3, 5 # 工厂数、客户数 c = np.array([1.0, 1.2, 0.9]) # 第一阶段产能单位成本 b = 1.0 + 0.1 * np.arange(I)[:, None] + 0.05 * np.arange(J)[None, :] d0 = np.full(J, 100.0) # 客户基础需求 delta = np.array([30., 40., 25., 35., 30.]) # 需求波动幅度 # 盒式不确定集顶点：u_j 取 0 或 1 vertices = list(product([0.0, 1.0], repeat=J)) # 收敛参数 eps = 1e-4 max_iter = 30

为什么把不确定集顶点枚举出来？因为第二阶段是线性规划，最坏场景一定落在盒式不确定集的某个顶点上。这里 J=5，顶点一共 (2^5=32) 个，枚举成本很低。如果你的问题顶点数爆炸，后面第 4 章会给替代方案。

3.2 子问题 SP：枚举顶点 + LP 求解

子问题要做的事是：给定第一阶段决策 x_val，找出使第二阶段成本最大的那个需求场景。代码里对 32 个顶点逐一求解一个 LP，取最大值。

def solve_sp(x_val): """给定第一阶段决策 x_val，求最坏场景和第二阶段成本""" best_val = -np.inf best_u = None for u in vertices: u_arr = np.array(u) demand = d0 + u_arr * delta # 当前场景下的需求 y = cp.Variable((I, J), nonneg=True) cons = [ cp.sum(y, axis=0) >= demand, # 每个客户需求必须满足 cp.sum(y, axis=1) <= x_val, # 每个工厂产能不能超 ] prob = cp.Problem(cp.Minimize(cp.sum(b * y)), cons) prob.solve() if prob.status == 'optimal' and prob.value > best_val: best_val = prob.value best_u = u_arr return best_val, best_u

核心逻辑就两层：外层把每个顶点转成具体需求向量，内层用 CVXPY 解一个标准 LP。代码里两个约束的方向别搞反——需求约束是“大于等于”，产能约束是“小于等于”，一旦写反，子问题要么无解要么给出虚假的负成本。参数说明：x_val 是一阶段产能决策数组，长度等于工厂数；返回值 best_val 是最坏场景下的第二阶段成本，best_u 是对应的最坏场景，后续要作为“割”加进主问题。

3.3 主问题 MP 与主循环：三个必须盯住的参数

主问题每轮重建一次。注意一个容易踩的坑：(\eta) 在没有场景约束时必须有一个下界，否则主问题是无界的。第二阶段生产成本非负，所以初始化时直接加 (\eta \ge 0)。

def solve_mp(): """求解主问题，返回一阶段决策、eta 值和主问题目标""" x = cp.Variable(I, nonneg=True) eta = cp.Variable() cons = [eta >= 0] # 防止第一轮无界 obj = cp.sum(c * x) + eta for u in added_u: # added_u 是已发现的最坏场景列表 yl = cp.Variable((I, J), nonneg=True) demand_l = d0 + u * delta cons += [ cp.sum(yl, axis=0) >= demand_l, cp.sum(yl, axis=1) <= x, eta >= cp.sum(b * yl), # C&CG 割：eta 要能覆盖该场景成本 ] prob = cp.Problem(cp.Minimize(obj), cons) prob.solve() if prob.status != 'optimal': raise ValueError('master problem failed: ' + prob.status) return x.value, eta.value, prob.value # 主循环 added_u = [] UB, LB = 1e8, -1e8 for k in range(max_iter): x_val, eta_val, mp_obj = solve_mp() LB = mp_obj # 主问题给出下界 sp_val, u_star = solve_sp(x_val) # 固定 x 找最坏场景 UB = min(UB, np.dot(c, x_val) + sp_val) print(f'iter {k}: LB={LB:.4f}, UB={UB:.4f}, gap={(UB-LB)/UB:.6f}') if UB - LB < eps: break added_u.append(u_star) # 把最坏场景加进主问题

三个必须盯的参数：

• eps 收敛容差。设成 1e-4 对绝大多数产能问题够了，设太小容易在最后几轮反复震荡，因为数值求解器本身有 1e-8 量级的误差。
• max_iter 迭代上限。正常盒式不确定集几十轮内必然收敛，如果顶到上限还没停，先检查 SP 是不是真的求到了全局最坏场景，而不是盲目加大上限。
• η 的初始下界。这个最容易翻车，删掉eta >= 0，第一轮主问题直接报 unbounded。

3.4 同一框架下切换到 Benders 割

Benders 和 C&CG 的主循环几乎一样，差别只在 SP 返回值。Benders 的 SP 除了要最坏场景值，还要取对偶乘子，用来构造割。

def solve_sp_benders(x_val): """返回最坏场景值、场景以及该场景对应的对偶乘子""" best_val = -np.inf best_u = best_alpha = best_beta = None for u in vertices: u_arr = np.array(u) demand = d0 + u_arr * delta y = cp.Variable((I, J), nonneg=True) cons = [ cp.sum(y, axis=0) >= demand, # 对应乘子 beta >= 0 cp.sum(y, axis=1) <= x_val, # 对应乘子 alpha >= 0 ] prob = cp.Problem(cp.Minimize(cp.sum(b * y)), cons) prob.solve() if prob.status == 'optimal' and prob.value > best_val: best_val = prob.value best_u = u_arr best_alpha = cons[1].dual_value.sum(axis=1) # 每个工厂的产能影子价格 best_beta = cons[0].dual_value # 每个客户的需求影子价格 return best_val, best_u, best_alpha, best_beta

拿到对偶乘子后，主问题里不再加 y 变量，而是加一条标量割：

# Benders 割：eta >= beta^T * demand - alpha^T * x demand_star = d0 + u_star * delta benders_cut = eta >= best_beta @ demand_star - best_alpha @ x

这里 alpha 来自产能约束cp.sum(y, axis=1) <= x_val，beta 来自需求约束cp.sum(y, axis=0) >= demand。符号必须和 LP 对偶推导一致，写成beta^T demand - alpha^T x，方向反了主问题会一路降到负无穷。Benders 每轮只加一条不等式，主问题规模小，但收敛轮数通常比 C&CG 多，在这个算例上大概要多 1.5～2 倍的迭代。

4. 把骨架改成你自己的业务问题：四个可扩展改编切入点

4.1 不确定集改动：从盒式到预算不确定集

盒式不确定集的毛病是过度保守：它假设所有客户需求同时拉满，现实中这几乎不可能。预算不确定集用 (\Gamma) 控制同时变坏的需求个数：(\sum_j u_j \le \Gamma)。代码改动非常小——在枚举顶点后加一个过滤条件。

Gamma = 2 # 最多 2 个客户需求同时取坏值 vertices = [u for u in product([0.0, 1.0], repeat=J) if np.sum(u) <= Gamma]

这个过滤必须在枚举阶段做，不能在 SP 内部做。如果在 SP 内对每个顶点判一下再跳过，逻辑上没错，但会浪费大量时间在无效 LP 求解上。Gamma 的取值玄学在于：太小体现不出鲁棒性，太大又回到盒式的保守。常见做法是取 J 的 40%～60%，然后用敏感性分析检验结果对 Gamma 的响应是否平滑。

如果 J 很大，枚举顶点不现实，需要改写成对偶子问题直接求解。此时 SP 会变成一个含 (u_j \beta_j) 双线性项的非凸问题，常见做法是针对小规模场景枚举、大规模场景用大 M 线性化或启发式迭代。这个切换没有统一模板，但骨架里 SP 函数是唯一需要动的地方，主循环不用改。

4.2 第二阶段从 LP 改成 MILP：C&CG 仍然成立

实际业务里第二阶段常常有整数变量：是否启用某条产线、是否租用某辆卡车、订单是否拆分。此时把第二阶段 y 的一部分变量声明为cp.Variable(shape, integer=True)或boolean=True，SP 用支持 MILP 的求解器（如 GLPK_MI、CBC）求解即可。

C&CG 割的结构完全不受影响，因为割只是“把场景原样加入主问题”，不依赖对偶推导。这也是 C&CG 比 Benders 更适合做扩展的根本原因。Benders 到这里就比较痛苦了：整数变量的 LP 对偶不再有效，需要做 MILP 对偶或 Gomory 割，工程难度立刻上一个台阶。

如果你只是给现有确定性 MILP 加鲁棒性，优先走 C&CG，别在 Benders 上死磕。

4.3 多周期与多产品扩展：变量维度增加

从单周期扩到多周期是生产计划类问题最常见的需求。改动核心是把 y 增加时间维度，并在主问题里加入库存平衡约束。

T = 12 # 12 个月 y = cp.Variable((T, I, J), nonneg=True) inv = cp.Variable((T + 1, J), nonneg=True) # 库存平衡：上月库存 + 本月生产 - 本月需求 = 本月库存 for t in range(T): cons += [inv[t + 1] == inv[t] + cp.sum(y[t], axis=0) - demand_l[t]] cons += [inv[0] == init_inv, inv[T] >= final_inv]

注意不确定参数从向量变成矩阵 (u_{t,j})，最坏场景的枚举复杂度会指数上升。实际操作中，更常见的做法是按周期共享一个预算约束 (\sum_{t,j} u_{t,j} \le \Gamma)，而不是每期独立。多产品只是把 I、J 的维度拆细，骨架不需要动。

4.4 数据接口与代码解耦：从硬编码到通用参数

小算例里硬编码没问题，上业务数据前最好把数据访问和求解逻辑拆开。我一般会做一个RobustModel类，构造时接收 c、b、d0、delta、U_type、Gamma，把 SP 和 MP 作为内部方法。数据来源换成 pandas 读 CSV 或者直接从数据库拉，工厂数、客户数用len()推断，而不是写死。

代码解耦的收益在调参时最明显：你要试不同 Gamma、不同波动幅度、不同收敛精度，如果参数散落在代码各处，每次改动都要通读一遍；集中成配置后，敏感性分析只是循环改参数的事。另外，建议把每次迭代的 LB、UB、gap 记录成表格输出，这是后面验证结果正确性的基础。

5. C&CG 与 Benders 复现避坑：5 条最容易被数值和符号坑死的细节

5.1 收敛与间隙：最容易让 CVXPY 直接报 unbounded 的两个场景

现象一：主问题第一次求解就报Problem unbounded，或者返回的 x 全是 0、eta 是负无穷。

原因：主问题里 (\eta) 没有任何约束时可以自由取负，目标函数趋向负无穷。这个坑在几乎所有 C&CG 新手代码里都会出现，因为模型公式里根本没提 eta 的边界。

解决：初始化时加eta >= 0。如果你的第二阶段成本可以为负（比如有残值回收），就用一个足够小的下界，比如eta >= -1e6，而不是硬编码 0。

现象二：迭代了十几轮，gap 一直在 1e-3 附近震荡，怎么都压不到设定阈值。

原因：UB 更新用了np.dot(c, x_val) + sp_val，但 LB 直接取了主问题目标值。主问题里的 x 每轮都在变，主问题目标并不是单调递增的，LB 可能出现回退；同时 SP 如果没求到真正的全局最坏场景，UB 会偏小，gap 被虚假压缩。

解决：UB 必须用min(UB, candidate)缓存，LB 用max(LB, candidate)缓存，这是两阶段鲁棒优化的标准做法。另外把每轮迭代的 x、场景、对偶乘子都打印出来，人工核对一下最坏场景是否合理——SP 枚举漏了顶点时，gap 会异常小，但解是错的。

5.2 符号与建模：Benders 割方向和不确定集边界问题

现象三：Benders 版本的主问题目标一路下降，很快变成负值，而 SP 返回的成本明明是正的。

原因：Benders 割的符号反了。对偶子问题目标写成 (\beta^T demand - \alpha^T x)，如果主问题里写成 (\beta^T demand + \alpha^T x)，割的方向就错了。这个错误用眼睛很难看出来，因为单看每条割都像模像样。

解决：写一个 x=0 的退化测试。x=0 时产能为 0，SP 只有在需求也为 0 的场景下可行，最优值应为 0；此时 Benders 割必须满足 (\eta \ge 0)。如果割给出负值，说明对偶符号推导有问题。再进一步，用一个小规模确定性算例，把不确定集固定为单点，鲁棒最优应等于确定性最优，凡是 Benders 版本结果对不上，优先怀疑符号。

现象四：加了预算不确定集后，鲁棒最优成本反而比确定性最坏情况还高，甚至不可行。

原因：枚举顶点时忘了过滤 (\sum u_j \le \Gamma)。顶点枚举基于全组合生成，里面包含了超过预算的场景，这些场景在真实不确定集中根本不存在。SP 一旦用越界场景构造割，主问题就被过度约束，结果过于保守，严重时直接不可行。

解决：在生成顶点时过滤，而不是在 SP 循环里判断。代码片段见 4.1 节。如果 J 很大无法枚举，改用 SP 直接建模不确定集约束，并检查 SP 返回的 u_star 是否满足 (\sum u_j \le \Gamma)。

现象五：用大 M 线性化双线性项时，求解器要么报数值警告，要么结果对 M 极其敏感。

原因：M 的尺度没标定。M 太小会切掉真实最优解，M 太大则系数矩阵条件数恶化，单纯形法和内点法都会出现精度问题。这基本是玄学，但有个工程经验：M 取变量实际范围乘积的量级，而不是凭感觉给 1e6。比如 (u_j \in [0,1])、(\beta_j) 的合理上界是产品成本的最大值，M 取这个上界乘以 1.2 就够了。

解决：能不用大 M 就不用。盒式不确定集枚举顶点是首选，预算不确定集对偶离散化处理，只有 J 很大且问题结构复杂时才考虑大 M，并且一定要做 M 的敏感性扫描。

6. 验证复现结果的两个低成本办法：退化算例与完全等价比对

复现别人论文里的鲁棒优化代码，最怕的是“跑通了但结果是错的”。两个不需要真实数据的验证办法，能在半小时内暴露 90% 的隐藏问题。

第一个办法：让不确定集退化。把 delta 全部设成 0，不确定集只剩一个点，两阶段鲁棒优化退化成普通确定性 LP。此时 C&CG 应该在第一轮就收敛——主问题第一次求解得到 x，SP 在最坏场景（唯一场景）下算出的成本应该正好等于主问题的 eta，UB 等于 LB。如果这个测试不通过，说明 SP 和 MP 的场景对应关系有问题，比如需求方向的约束写反了，或者 eta 的割没生效。

第二个办法：利用“所有系数同号”的单调性。当生产成本 b 和波动幅度 delta 全为正时，最坏场景必然是每个客户需求都取上限，即 (u_j=1)。那么鲁棒问题的最优解应该恰好等于把需求直接替换成 (d0 + delta) 的确定性模型的最优解。两个模型独立求解，比较最优值和 x，完全一致才算通过。这个办法对盒式不确定集尤其有效，因为它不依赖迭代过程，只比较最终结果。

验证脚本里再留一个一致性检查：取 C&CG 最后一轮的 x* 和 u*，重新调用一次 SP，返回值应该等于最后一轮的 UB，误差在 1e-6 量级。对不上说明主循环的 UB 更新逻辑有 bug，比如用了旧场景而不是新场景。

这类检查现在是我每次复现的固定动作。早期翻车最狠的一次，就是预算不确定集忘了过滤顶点，结果所有验证都对得上，唯独真实业务数据下结果比预期保守 20%，查了两天才定位到是枚举越界。从那以后，退化算例先跑、单调性比对再做、SP 一致性最后查，三步走完我才敢把代码接到业务数据上。希望帮到你。

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