企业官网建设流程全解析

做供应链、能源、调度这些方向的读者，这两年应该没少见过“两阶段鲁棒优化”这个词。它解决的是这样一类问题：你现在要花钱做决策，但未来的需求、价格、可用容量全都不确定，好在事情发生后还有第二阶段的调整机会。真正上手之后你会发现，两阶段鲁棒模型本身不难写，难的是把那个嵌套的 max-min 结构变成计算机能处理的形式，而这恰恰是分布鲁棒、KKT函数这些工具派上用场的地方。这篇文章用一个小而完整的算例，把两阶段鲁棒优化从建模到代码实现完整过一遍，重点放在内层问题的 KKT 改写、分布鲁棒怎么扩展，以及能照着改的经典代码。适合正在入门鲁棒优化、或者已经有基础但被子问题求解卡住的同学。

1. 两阶段鲁棒优化到底在优化什么

1.1 模型结构与两个决策阶段

先看两阶段鲁棒优化最通用的形式：

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

这个式子看着绕，拆开就三块。

第一阶段的 x 是“现在就要拍板的决策”，也叫 here-and-now 决策。比如建多大容量的仓库、签多少长期运力合同、在哪些点位设产能。这个决策会在看到不确定性之前定下来，所以要给成本 c^T x。

第二阶段的 y 是“看到运气之后还能再调的决策”，叫 wait-and-see 决策。比如某个城市突然需求暴增，可以临时租车、加班生产、调货救急。这个决策针对的是已经实现的不确定性 u，所以可以依据 u 的具体值来调整，成本是 d^T y。

中间的 max_{u∈U} 是鲁棒优化的灵魂。它不是拍一个最可能的场景，而是把所有场景放在集合 U 里，然后找最坏的那个。也就是说，第一阶段决策要让“最坏情况下第二阶段的补救成本”也尽量小。

用一句话概括整个模型：现在怎么定，才能让将来无论遇到什么情况，总成本都可控。

1.2 为什么内层是一个 max-min 而不是一个 min

初学者最容易懵的地方是：为什么 u 和 y 不是同一个决策者？

想象一个场景。你是物流公司老板，第一阶段决定租几个常规仓库，x 是常规仓容。第二阶段是双十一当天，某个区域的需求 u 突然爆掉了。这时候你的调度系统会怎么做？它会在已知 u 的具体值之后，通过临时调车、外部仓补货等手段，找到当前情况下成本最小的补救方案，这就是内层 min_y。但因为你在第一阶段做决策时还不知道 u 是 1 还是 3，你只能假设老天爷会挑一个最刁钻的 u 来坑你，这就是外层 max_u。max 和 min 不是两个决策者在博弈，而是同一个问题里“先不确定性、后适应性优化”的自然顺序。

所以 max-min 不是故弄玄虚，它是两阶段决策机制的真实产物。第一阶段模型看到的不是某个具体 u，而是整个 U 的“恶意选择”。

从算法角度，这个嵌套结构才是真正的麻烦。普通单层优化只要交给求解器就行，但 max-min 是一个双层结构。如果不把它转化掉，Gurobi、CPLEX 这些主流求解器都没法直接吃进去。

2. 内层问题怎么变成可计算的形态

2.1 从一个具体算例开始

后面所有代码都围绕下面这个算例展开，我把参数控制得非常小，方便手算和验证。

第一阶段只做一个决策 x，表示自有资源容量，单位成本为 1：

min_{x ≥ 0} x + max_{u ∈ U} min_{y ≥ 0} (2 y1 + 3 y2 + 100 y3)

第二阶段有三个决策变量：

• y1：使用自有资源，边际成本 2，上限是第一阶段决策 x；
• y2：使用外部稳定资源，边际成本 3，上限取决于外部可用量 u2；
• y3：紧急高价资源，边际成本 100，没有上限，用于兜底保证问题恒可行。

第二阶段约束是：

y1 + y2 + y3 ≥ u1 y1 ≤ x y2 ≤ u2 y ≥ 0

不确定性集合取一个简单的盒子：

U = { u = (u1, u2) | 1 ≤ u1 ≤ 3, 2 ≤ u2 ≤ 4 }

这个模型很典型。u1 是需求，需求越大越难；u2 是外部可用量，外部可用量越小越难。最坏场景直觉上就是让 u1 尽量大、u2 尽量小，也就是 u = (3, 2)。

为什么最坏场景恰好是顶点？后面用强对偶一推就清楚了。

2.2 强对偶改写：化简 max-min 的标准动作

只要内层 min_y 是线性规划，就能用强对偶把 min 消掉。

内层 LP 写成标准形式：

min 2 y1 + 3 y2 + 100 y3 s.t. y1 + y2 + y3 ≥ u1 y1 ≤ x y2 ≤ u2 y ≥ 0

对偶变量分别记作 λ1、λ2、λ3，对应前三个不等式约束，且 λ ≥ 0。对偶问题写成：

max u1 λ1 - x λ2 - u2 λ3 s.t. λ1 - λ2 ≤ 2 λ1 - λ3 ≤ 3 λ1 ≤ 100 λ ≥ 0

这里有个符号细节容易错：原问题是 y1 ≤ x、y2 ≤ u2，写成 ≥ 形式的时候要取相反数，所以对偶目标里是 -x λ2 和 -u2 λ3。我看到过不少人在这一步把符号写反，导致最坏场景判断完全错误。

现在整个两阶段问题的嵌套部分就变成了：

max_{u ∈ U, λ ∈ Λ} (u1 λ1 - x λ2 - u2 λ3)

其中 Λ 是上面那三个线性约束加上非负约束定义的集合。

虽然目标里出现了 u1 λ1 和 u2 λ3 这种双线性乘积，但这个结构反而把最坏场景的性质暴露得很清晰：λ1 ≥ 0 且 λ3 ≥ 0，所以对固定的乘子 λ，目标关于 u1 单调递增、关于 u2 单调递减。u1 越大就越坏，u2 越小就越坏。因此最坏场景一定落在盒子的顶点上，具体来说就是 u1 = 3 且 u2 = 2。

这个性质是强对偶赠送给我们的额外红利。当不确定性集合是盒式或者一般多面体，并且内层 LP 对偶可行域不依赖 u 时，最坏场景往往可以在顶点里找。实际操作中，很多团队偷懒直接枚举盒子四个顶点，效果也很好。

换一种更书面的说法：这个转化把两层问题变成了一个在 U 和 Λ 的笛卡尔积上求最大化的问题，虽然双线性项让整体非凸，但 U 的顶点枚举降低了难度。

2.3 KKT 条件重构：YALMIP kkt 函数的使用

强对偶好用，但它有个前提：内层必须是 LP，并且你愿意手工推导对偶。如果内层带非线性约束、或者模型参数一变就要重新推导，强对偶就很麻烦。这时候第二条路是 KKT 条件重构。

KKT 条件的本质是刻画“什么情况下 y 是内层问题的最优解”。把内层最优性用一组方程和不等式写出来之后，max-min 问题就变成一个带平衡约束的优化问题，学术上叫 MPEC（Mathematical Program with Equilibrium Constraints）。

以刚才的算例为例，内层问题约束写成 g(y) ≤ b 形式：

-y1 - y2 - y3 + u1 ≤ 0 y1 - x ≤ 0 y2 - u2 ≤ 0 -y1 ≤ 0 -y2 ≤ 0 -y3 ≤ 0

引入乘子 λ1 到 λ6，都要求非负。KKT 条件包含三块。

第一块是平稳性条件，对目标函数和约束求梯度：

2 - λ1 + λ2 - λ4 = 0 3 - λ1 + λ3 - λ5 = 0 100 - λ1 - λ6 = 0

第二块是原始可行性和对偶可行性，也就是 y 满足所有不等式、λ ≥ 0。

第三块是互补松弛条件，每组乘子和它对应的约束不能同时“松”：

λ1 (u1 - y1 - y2 - y3) = 0 λ2 (y1 - x) = 0 λ3 (y2 - u2) = 0 λ4 y1 = 0 λ5 y2 = 0 λ6 y3 = 0

这六个乘子里，λ4、λ5、λ6是从非负约束来的，看起来多余，但写 KKT 时千万不能漏，否则平稳性条件不完整。我自己调试时吃过这个亏，少了一个乘子导致对偶目标对不上。

好消息是，手工写这一堆东西很容易出错，YALMIP 里专门提供了 kkt 函数。它的基本用法是：

[K, details] = kkt(inner_cons, y, inner_obj);

• inner_cons 是内层问题的约束；
• y 是内层决策变量；
• inner_obj 是内层目标函数；
• 输出 K 是一个约束对象，直接包含了平稳性、可行性、互补松弛条件。

于是我可以用下面这段代码把刚才算例的内层最优性全部写出来，不用手工敲任何一个乘子：

u = sdpvar(2, 1); y = sdpvar(3, 1);

inner_cons = [y >= 0, ... y(1) + y(2) + y(3) >= u(1), ... y(1) <= xval, ... y(2) <= u(2)];

inner_obj = 2y(1) + 3y(2) + 100*y(3);

[K, ~] = kkt(inner_cons, y, inner_obj);

但要特别提醒：kkt 函数生成的约束里包含互补条件，也就是乘子和约束的乘积等于 0。这种乘积约束会让问题变成非凸 MPEC，Gurobi、CPLEX 默认不能直接处理。实际工程里要么用专用求解器比如 Baron、KNITRO、bmibnb，要么手动把互补条件用 big-M 和二值变量线性化。所以 kkt 函数更适合内层是二次规划、非线性规划，或者你在做算法验证的时候用。

3. 从最坏场景到最坏分布：分布鲁棒扩展

3.1 鲁棒优化和分布鲁棒的边界

两阶段鲁棒优化处理的是“在 U 里取最坏场景”，它完全不看概率。这种做法的优点是稳健，缺点是太保守：如果某个坏场景出现的概率只有十万分之一，鲁棒优化也会按它来设计。

分布鲁棒优化（DRO）补上了这个缺口。它不假设你手头那个经验分布就是真实分布，而是把真实分布限制在一个模糊集合 D 里，然后优化最坏分布下的期望目标。翻译成两阶段模型就是：

min_{x ∈ X} c^T x + sup_{P ∈ D} E_P[ Q(x, u) ]

其中 Q(x, u) 就是内层 min_y 的最优值。

如果 D 是只包含一个点，DRO 退化成随机规划；如果 D 是所有以 U 为支撑集的分布集合，DRO 就退化回经典鲁棒优化。DRO 本质上是随机规划和鲁棒优化的中间地带。

我用一个生活类比解释：鲁棒优化是“不管天气预报怎么说，都按最大暴雨准备”；随机优化是“天气预报只有一种，就按这个预报算”；分布鲁棒是“手头有十年历史天气数据，但你知道气候模式可能变了，所以把可能的天气分布圈一个范围，在这个范围里找最坏期望损失”。

3.2 Wasserstein 模糊集：把不确定的分布拉回有限维

模糊集合 D 有很多种构造方式，最常见的是矩约束集合、phi-散度球和 Wasserstein 球。

矩约束集合限制分布的均值和协方差只能在一定范围内变化，理论漂亮，但实际求解时经常把模型搞得很复杂。phi-散度球是定义在离散分布上的，处理离散场景很方便，但没法自然刻画连续分布。

工程里这几年最火的是 Wasserstein 模糊集。给定一组历史样本 u_1, u_2, ..., u_N，经验分布是 \hat{P}_N，Wasserstein 模糊集定义为：

D = { P | W(P, \hat{P}_N) ≤ ε }

其中 ε 是模糊半径，W 是 Wasserstein 距离。它表达的意思是：真实分布可以偏离经验分布，但偏离的“搬运成本”不能超过 ε。

为什么 Wasserstein 模糊集好用？因为它有一组漂亮的对偶性质。对于 Lipschitz 连续的 Q(x, u)，可以写成：

sup_{P ∈ D} E_P[ Q(x, u) ] = inf_{λ ≥ 0} { λ ε + (1/N) Σ_{i=1}^N sup_{u ∈ U} [ Q(x, u) - λ d(u, u_i) ] }

这个公式等于把“最坏分布”的搜索转化为：对每个历史样本，在它附近做一个带惩罚的局部最大化解；λ 是一个全局惩罚系数，起到了把概率质量从样本点搬走的“运费单价”的作用。

3.3 分布鲁棒两阶段问题为什么还能用前面的套路

注意看上面公式里的 sup_{u ∈ U} [ Q(x, u) - λ d(u, u_i) ]。Q(x, u) 本身又是一个内层 LP 的最优值。所以整体上，Wasserstein DRO 把一个分布层面的不确定性，转化成了 N 个“局部两阶段鲁棒子问题”的加权组合。这正好能把前面的 KKT、强对偶、C&CG 整套技术搬过来。

如果你选的 d(u, u_i) 是范数，比如 l_1 或 l_∞，那么 [ Q(x, u) - λ d(u, u_i) ] 这个 sup 依然可以写成有限维线性约束；如果 Q(x, u) 用 KKT 重构，那整个问题就是一组带互补约束的优化问题。这就是为什么我建议先搞懂两阶段鲁棒，再去看 DRO 论文会轻松很多，因为 DRO 的核心算法并没有完全脱离鲁棒优化的工具箱。

不过说实话，分布鲁棒两阶段模型的实际求解目前仍然不轻松。N 个局部子问题会让约束数量成倍增加，C&CG 主问题规模会膨胀，通常需要配合场景缩减、并行计算和 Benders 分解来加速。入门阶段先把两阶段鲁棒跑通，再往 DRO 走，路径会顺很多。

4. 经典代码：C&CG 框架逐行拆解

4.1 C&CG 主循环：LB、UB 怎么更新

C&CG（Column-and-Constraint Generation）是求解两阶段鲁棒优化最主流的分解算法之一。它的核心思想是：不确定性集合 U 里的场景太多，不可能一下子全放进模型，那就通过子问题每次找出当前最坏场景，把最坏场景“切成一个约束”加回主问题。

主问题（MP）在第 k 次迭代时维护一个已发现的场景集合 S = {u^1, ..., u^k}：

min x + θ s.t. 对每个 u^j ∈ S 都存在对应的 y^j 满足内层约束，且 θ ≥ 2 y1^j + 3 y2^j + 100 y3^j x ≥ 0, θ ≥ 0

求解主问题会得到当前场景集合下的最优解 x 和目标值，这个目标值是一个下界 LB，因为场景集合没有铺满整个 U，问题比原问题更宽松。

子问题（SP）固定 x，求解：

Q(x) = max_{u ∈ U} min_{y} (2 y1 + 3 y2 + 100 y3)

得到最坏场景 u* 和 Q(x)。用 x + Q(x) 更新上界 UB，因为这是原问题的一个可行解。

如果 UB - LB 小于阈值，就认为收敛；否则把 u* 加入场景集合 S，重新求解主问题。

C&CG 和 Benders 分解长得很像，区别在于：Benders 加的是对偶割，C&CG 加的是第二阶段变量和原始约束。C&CG 在处理两阶段鲁棒时收敛更快，尤其当不确定性集合是多面体时，往往迭代几次就能收敛。

4.2 主问题代码：YALMIP + Gurobi 实现

我用的环境是 MATLAB + YALMIP + Gurobi，这是做鲁棒优化原型验证最顺手的组合。Python 用户可以把 YALMIP 换成 Pyomo，逻辑完全一样。

主问题的 YALMIP 代码：

x = sdpvar(1, 1); theta = sdpvar(1, 1);

cons = [x >= 0, theta >= 0]; obj = x + theta;

for k = 1:size(S, 2) y = sdpvar(3, 1); u1 = S(1, k); u2 = S(2, k); cons = [cons, ... y >= 0, ... y(1) + y(2) + y(3) >= u1, ... y(1) <= x, ... y(2) <= u2, ... theta >= 2y(1) + 3y(2) + 100*y(3)]; end

optimize(cons, obj, sdpsettings('solver', 'gurobi', 'verbose', 0));

xk = value(x); LB = max(LB, value(obj));

这段代码有三个细节值得说。

第一，theta 必须有一个下界 0。这个算例第二阶段成本本来就是非负，加上这个约束可以防止主问题在没加任何场景时出现无界。很多初学者省略这一步，结果迭代初期 LB 直接是 -inf，收敛判据直接失效。

第二，每个场景都有自己的 y 变量。这是 C&CG 的标准做法，y 只在对应场景下有意义，不能所有场景共用一个 y，否则会人为地把第二阶段决策限制成场景无关，得到一个过于保守的错误解。

第三，主问题的 x 是唯一的全局变量，theta 也是。不要在每个场景循环里重复定义 theta，否则 YALMIP 会给你建立一个全部变量的和 theta1 + theta2 + ... 的目标，结果完全不对。

4.3 子问题代码：顶点枚举版本

子问题理论上可以用强对偶转成单层问题，但在盒子不确定性下有个更简单粗暴的方法：枚举顶点。这个算例 U 的顶点只有四个：u1 ∈ {1, 3}，u2 ∈ {2, 4}。

子问题枚举代码：

function [Q, u_worst] = subproblem_vertices(xval) vertices = [1 2; 1 4; 3 2; 3 4]; Q = -inf; u_worst = []; for i = 1:size(vertices, 1) u = vertices(i, :)'; y = sdpvar(3, 1); cons = [y >= 0, ... y(1) + y(2) + y(3) >= u(1), ... y(1) <= xval, ... y(2) <= u(2)]; optimize(cons, 2y(1) + 3y(2) + 100y(3), ... sdpsettings('solver', 'gurobi', 'verbose', 0)); q = value(2y(1) + 3y(2) + 100y(3)); if q > Q Q = q; u_worst = u; end end end

用这个子问题函数，主循环就是标准的 C&CG：

S = [1; 2]; % 初始场景 LB = -1e9; UB = 1e9;

for iter = 1:20 % 求解主问题，得到 xk、LB ... % 求解子问题 [Q, u_worst] = subproblem_vertices(xk);

UB = min(UB, xk + Q); if UB - LB < 1e-4 break; end % 去重后加入场景集合 if ~ismember(u_worst', S', 'rows') S = [S, u_worst]; end

end

运行结果会看到迭代次数非常少。因为最坏场景恰好是顶点 (3, 2)，初始场景 (1, 2) 被加入一次后，第二次迭代就能收敛。具体数字是：x* 在 [1, 3] 内部时总成本都是 9，代码会收敛到 x 接近某个可行点，UB 和 LB 都稳定在 9 附近。

4.4 子问题代码：KKT 函数版本

如果内层是二次规划或非线性规划，顶点枚举就失效了。这时候可以用 YALMIP 的 kkt 函数重写子问题。

子问题的思路是：给定 x，搜索 u 和内层 KKT 满足条件的最优 y，目标是最大化内层最优值。

function [Q, u_worst] = subproblem_kkt(xval) u = sdpvar(2, 1); y = sdpvar(3, 1);

inner_cons = [y >= 0, ... y(1) + y(2) + y(3) >= u(1), ... y(1) <= xval, ... y(2) <= u(2)]; inner_obj = 2*y(1) + 3*y(2) + 100*y(3); [K, ~] = kkt(inner_cons, y, inner_obj); sp_cons = [K, u >= [1; 2], u <= [3; 4]]; sp_obj = -(2*y(1) + 3*y(2) + 100*y(3)); optimize(sp_cons, sp_obj, sdpsettings('solver', 'bmibnb', 'verbose', 0)); Q = value(2*y(1) + 3*y(2) + 100*y(3)); u_worst = value(u);

end

注意这段代码本质上是一个求解 MPEC 的问题，因为 K 里面带了互补松弛条件。bmibnb 是 YALMIP 自带的全局非线性求解器，适合小规模问题；规模一大就很慢。更工程化的做法是把互补松弛条件用 big-M 线性化，让 Gurobi 去解一个混合整数线性规划。

big-M 线性化的思路长这样：对于 λ2 * (y1 - x) = 0 这种条件，引入二值变量 z，写成：

y1 - x ≤ M z λ2 ≤ M (1 - z)

M 的取值是关键。M 太大数值稳定性差，太小会丢掉可行解。通常先求解一次松弛问题拿到乘子和变量的量级，再据此设置 M。这个环节非常考验经验，我一般会在代码里留一个可配置参数 M，做敏感性实验。

我想再强调一遍：kkt 函数版本在教学和验证算法思路时非常好用，但真拿去求解大规模工程问题前，一定要评估求解器对 MPEC 的支撑能力。很多商业求解器连“这个模型里有互补约束”都识别不出来，容易给出荒谬结果，这时候你想Debug 都无从下手。

5. 实战经验：求解器、数值坑与效率优化

5.1 求解器配置与模型规模取舍

两阶段鲁棒优化像样的模型规模都很容易膨胀。主问题每加一个场景，就会新增一组第二阶段变量。场景加到几百个之后，YALMIP 建模的开销都会肉眼可见地变大。

我的建议是先把子问题做对、做快，再考虑主问题并行。子问题如果只是 LP，顶点枚举或者强对偶 LP 都比 KKT 重构快一个数量级。子问题只有在内层不是 LP 时才上 kkt，且要预期到求解时间会显著上升。

求解器版本也会影响结果。Gurobi 9 之后对 LP 和 MILP 都稳，但 MPEC 类模型还是要靠 Baron、KNITRO 这类非线性规划求解器，或者老老实实做 big-M 线性化再扔回 Gurobi。

5.2 我踩过的几个坑

第一个坑是对偶符号。第一次写这个算例时，我把 y1 ≤ x 直接当作 ≤ 约束放进对偶，结果对偶目标里变成了 +x λ2，最坏 u2 被判断成了 4，整个 UB 比真实解低了整整 2。排查了很久才发现，原式是 ≤ 约束，转成标准 ≥ 形式时必须加负号。

第二个坑是互补松弛遗漏。手工写 KKT 的时候，非负约束 y ≥ 0 的乘子经常被顺手忽略。但对于像 y3 这种高成本紧急变量，恰恰是它的非负约束乘子在平稳性条件里起了决定性作用。少一个乘子，平稳性方程就不闭合，得到的 y 根本不是内层最优解。现在我只用 kkt 函数自动生成，不再手写。

第三个坑是场景去重。C&CG 迭代后期，子问题可能返回同一个最坏场景，如果不加去重，主问题会反复加入完全相同的约束，白白增加求解时间。我的去重判断用 ismember 按行比较，容差取 1e-5，实测很稳。

第四个坑是最坏场景在顶点之外的情况。很多人做完盒式不确定性后，以为所有两阶段鲁棒问题都能靠枚举顶点求解。这是不对的。当内层对偶可行域也依赖 u、或者不确定性集合不是多面体时，最坏场景可能出现在边界内部，顶点枚举会得到低估的 UB。一定要先验证内层对偶的可行域是否确实不依赖 u。

5.3 把算法跑得更快的小技巧

第一阶段可以先跑一个只含初始场景的主问题，得到 x 的初值，再进 C&CG，这样能省掉前几轮无用功。

子问题如果和顶点枚举等价，可以预先算好每个顶点的对偶最优解缓存起来。C&CG 每轮只改变 x 的值，不改变问题和顶点结构，缓存能带来明显加速。

如果对精度要求不高，收敛阈值不用设到 1e-6。事实上两阶段鲁棒模型的数据精度通常到 1e-3 就够了，阈值太小只会让迭代次数无谓地增加。

还有一个细节：YALMIP 里 S = [S, u_worst] 这种动态扩充矩阵，迭代次数多了会反复复制内存，慢但不致命。如果场景数量真的很大，最好预分配一个足够大的矩阵，再用计数器维护有效列数。

6. 常见问题速查表

这些坑我都在实际调试中踩过。尤其是对偶符号问题，看起来很基础，一旦错，整个优化结果全反，而且由于模型本身是嵌套结构，很难一眼看出来。我的习惯是每次改完模型，立刻构造一个手工可解的小算例去验证子问题，不验证清楚不进入下一步。

最后聊点我自己的习惯。两阶段鲁棒优化的代码实现，从模型搭建到 C&CG 跑通，其实半天时间就够；真正费时间的是那些看不见的数值细节和建模约定。我强烈建议新手不要直接拿一个大模型起步，而是先把我上面那个只有三个第二阶段变量的最小算例完整跑通，包括主问题、子问题、收敛判断、场景去重四个模块。等这套最小框架跑通之后，再把你的真实业务约束一层层加进去，每加一类约束，就回去验证一次最坏场景的合理性。分布鲁棒也一样，别一上来就上 Wasserstein 球，先用离散模糊集合的简单版本把 C&CG 的代码框架跑顺，再切换到连续分布集合。这样走，你踩的坑会少很多。