做分布鲁棒优化（Distributionally Robust Optimization, DRO）这个方向的人，大多数时候都在跟“不确定集合”较劲。刚入门的时候，我也踩过不少坑，最典型的困惑就是：随机规划明明把需求预测放进了期望值里，为什么调出来的解一到真实环境就翻车？答案其实不复杂——因为你在模型里假设的分布，跟真实分布根本不是一回事。分布鲁棒优化要解决的，就是“我不知道真实分布，但我绝不希望决策在某一个坏分布下崩掉”的问题。
这篇文章我想分享一个比较实用的简易模型：基于 Wasserstein 距离构造模糊集合的两阶段分布鲁棒模型，配合对偶转化和线性决策近似，最后用 Matlab 把整个求解流程跑通。代码实现部分我会给出可以直接复现的核心片段，也会把推导过程中的关键节点拆开讲清楚。适合刚接触 DRO、或者在随机规划里被“分布假设过强”坑过的人参考。
1. 为什么模糊集合要用 Wasserstein 距离而不是矩约束
1.1 分布鲁棒到底在鲁棒什么
先说个场景。假设你在做库存管理，每天的随机需求你可以用历史数据估计均值、方差，甚至能拟合出某种分布。传统随机规划的做法是：假设需求服从某个固定分布，然后优化期望成本。问题是，真实分布未必长这样，尤其是样本量不够大、需求有突变的时候，拟合出来的分布跟实际几乎对不上。于是你优化出来的订货量，可能在模型里很漂亮，一上线就库存缺货或者库存积压。
分布鲁棒的做法不是直接从“某个具体分布”出发，而是设定一个分布模糊集合，让真实分布不管落在集合里的哪一个位置，我们都要求决策方案不能太差。注意这里的“分布模糊集合”不是参数集合，而是由概率分布组成的一个集合。这个概念跟传统鲁棒优化里“参数落在区间内”的思路完全不一样，难度自然也上了一个档次。
那问题就来了：这个模糊集合应该怎么选？最早被大量使用的是矩约束型模糊集，比如限定均值和协方差的范围。用起来简单，但有明显缺点：矩约束只刻画了分布的“一阶二阶”信息，对形状完全不敏感，两个均值方差完全相同但尾巴行为完全不同的分布，在模糊集合里没法区分。真实场景里，需求的尾部风险往往才是要命的。
1.2 Wasserstein 球：把“分布之间的距离”变成可操作的约束
Wasserstein 距离，很多人第一次听到是在最优传输或者生成模型里。它直观理解就是：要把一个分布“搬运”成另一个分布，所需的最小平均搬运成本。搬得越少，两个分布越接近；搬运成本越大，两个分布差得越远。
在 DRO 里，我们用样本构造一个经验分布，然后以这个经验分布为球心，用一个 Wasserstein 球作为模糊集合：
[ \mathcal{D} = { P : W(P, \hat{P}_N) \le \epsilon } ]
其中 ε 就是模糊半径，代表“真实分布可能偏离历史经验多远”。这个集合的好处特别实在：它包含的分布不仅与经验分布有相近的均值、方差，还能在形状上贴近经验分布；而且只要半径 ε 取得合适，理论上可以覆盖真实分布的概率很高。这个保证不是玄学，是通过测度集中不等式推导出来的：样本量越大，经验分布离真实分布越近，所以 ε 可以取得越小。
这里有个要点：Wasserstein 距离要具体定义在一个运输代价上。常见做法是取某种范数作为地面距离，比如：
[ W(P, Q) = \inf_{\pi} \int | \xi - \xi' | , \pi(d\xi, d\xi') ]
如果你用 L1 范数，这叫 1-Wasserstein 距离；用 L2 范数，就是 2-Wasserstein 距离。模型里选哪种，会直接决定后面对偶转化出来的约束长什么样。简易模型通常用 1-Wasserstein，因为它的对偶形式最干净，Matlab 建模也最方便。
1.3 模糊半径如何标定
模糊半径 ε 不是拍脑袋定的。一个合理的标定方式，是根据样本量 N 和置信水平 β 来定。用测度集中不等式可以做一个保守的估计：
[ \epsilon_N(\beta) = C \cdot \sqrt{\frac{\log(1/\beta)}{N}} ]
不同分布、不同维数下，常数 C 不一样。对初学者来说，我更建议先做样本内灵敏度分析：把 ε 从 0 开始逐步增大，观察目标函数和决策变量的变化。ε=0 时，模型退化成普通样本均值随机规划；ε 越大，模型越保守，目标值越差。你选择一个“优化性能下降幅度能接受”的最大 ε，再结合历史数据的 out-of-sample 表现来校验，这是工程上最靠谱的做法。
还有一个非常容易踩的坑：Wasserstein 球的全称是以经验分布为球心、半径为 ε 的闭球。半径只影响“允许偏离多远”，不影响“往哪个方向偏”。如果你有领域知识，比如知道真实需求普遍比历史数据偏大，那应该在场景生成阶段就做调整，而不是靠放大 ε 来解决。放大 ε 只会让所有方向都变得保守，无法针对某个特定方向的偏差做修正。
2. 两阶段结构拆解与线性决策近似的适用边界
2.1 两阶段模型的标准形式
两阶段随机规划，最简单的理解就是“先见机行事，但只能部分见机”。第一阶段决策必须在随机变量实现之前做，用行话叫“在这里且现在”决策，比如提前定采购量、提前租仓库、提前安排班组；第二阶段决策是在随机变量实现之后做的，可以看作“补救动作”，比如需求多了就紧急补货、需求少了就打折处理。
标准的两阶段分布鲁棒模型长这样：
[ \min_{x} ; c^\top x + \sup_{P \in \mathcal{D}} \mathbb{E}_P [Q(x, \xi)] ]
其中 Q(x, ξ) 是第二阶段的成本函数：
[ Q(x, \xi) = \min_{y} ; q^\top y \quad \text{s.t.} ; T x + W y = \xi, ; y \ge 0 ]
这里的 ξ 是随机参数，T 和 W 分别是技术矩阵和补偿矩阵。Wasserstein 球带来的难点在于，内层的是一个“先对一个连续分布求期望，再求最大值”的问题，这个分布本身还是变量。如果不做任何变形，这本质上是一个无限维优化问题，没法直接塞进求解器。
为此，业界的主流路线就是对偶转化。把内层的 sup 问题转化为一个有限维的 min-max 问题，再进行化简。这个转化的意义不亚于把一座山搬走，它对模型的规模、求解难度、可解释性都影响巨大。我后面专门开了一节讲这个推导。
2.2 线性决策近似：把无限维决策变量拉回有限维
两阶段问题即使对偶转化完，第二阶段决策函数 y(ξ) 仍然是 ξ 的任意函数。为了求解方便，我们需要把 y(ξ) 限定在一个简单的函数族里。最常用的就是线性决策规则：
[ y(\xi) = y_0 + Y \xi ]
其中 y0 是常数向量，Y 是矩阵。也就是说，第二阶段的补救量是随机参数 ξ 的一个仿射函数。这个近似的合理性在于：很多实际问题的补偿成本是线性的，约束也是线性的，最优解虽然不是严格线性，但在可行域比较“规整”的情况下，线性决策已经能拿到相当好的性能。
有人可能问：为什么不用二次决策、分段线性决策？答案很简单：线性决策后面配合对偶转化，可以拼成一个一个凸优化问题，多项式时间内可以求解；一旦换成更高阶的决策规则，问题要么变成非凸，要么爆炸式增长，项目工期等不起。简易模型的价值就是快速验证想法，如果线性决策的结果已经能在 out-of-sample 测试里战胜随机规划，那就没必要上更复杂的东西。
2.3 为什么简易模型里选线性决策
我自己的体会是，两阶段分布鲁棒模型的复杂度，很多时候不在“模型有多高级”，而在“能不能让求解器落地”。线性决策规则虽然看起来简单，但它把原来一个函数优化问题真正变成了一个矩阵变量优化问题，决策变量的数量是确定的，可以预测的。这对于写 Matlab 代码来说非常重要，因为你不需要像求解随机规划那样用 Benders 分解或者逐步逼近算法，直接抛给求解器就能跑。
当然，必须牢记线性决策的适用边界。如果第二阶段的最优补偿动作本身和随机参数有强非线性关系，比如带整数变量、带强非线性成本、或者约束里有乘积项，那线性决策近似就比较危险。我在文末的避坑部分会专门给一个具体例子。
3. 对偶转化的完整推导
3.1 外层 sup 问题的对偶形式
这一节是整篇文章的“数学重镇”，也是很多人看论文看得最头疼的地方。我尽量把逻辑链讲清楚，但不做过于繁复的严格证明。
先只看内层：
[ \sup_{P: W(P, \hat{P}_N) \le \epsilon} \mathbb{E}_P [h(\xi)] ]
其中 h(ξ) 是一个与 x 有关的成本函数。在简易模型里，h(ξ) 通常就是第二阶段问题的最优值。根据 Wasserstein 距离对偶理论（本质上是 Kantorovich-Rubinstein 对偶），上面这个 sup 问题可以转化为：
[ \inf_{\lambda \ge 0} ; \lambda \epsilon + \frac{1}{N} \sum_{i=1}^{N} \sup_{\xi} \left[ h(\xi) - \lambda |\xi - \hat{\xi}_i| \right] ]
这里的 λ 是一个标量对偶变量，代表“模糊半径的边际成本”。这个转化是精确的，没有掉近似，只要地面距离取范数即可。好，现在问题从“在无穷多个分布里选一个最坏的”变成了“先枚举每个样本场景，对每个场景求一个 sup，再往它们头上加一个 λ ε 的惩罚”。
但是里面还剩一个 sup 对 ξ 的优化，这一步才是真正困难的地方，因为 ξ 是连续变量，范围可能是整个实数空间，不是有限个点。要消除这个 sup，就要用到 h(ξ) 的具体结构。
3.2 第二阶段问题的对偶与凸共轭
假设第二阶段问题是线性的：
[ h(\xi) = \min_{y \ge 0} { q^\top y : W y = \xi - T x } ]
把它写成其对偶形式，得到：
[ h(\xi) = \max_{\pi} { \pi^\top (\xi - T x) : W^\top \pi \le q } ]
把 h(ξ) 的这个表达代入前面 sup 问题里的 sup 部分：
[ \sup_{\xi} \left[ \max_{\pi} { \pi^\top (\xi - T x) } - \lambda |\xi - \hat{\xi}_i| \right] ]
根据 max 和 sup 可以交换（这里需要满足一定的凸性条件），等价于：
[ \max_{\pi: W^\top \pi \le q} \sup_{\xi} \left[ \pi^\top (\xi - T x) - \lambda |\xi - \hat{\xi}_i| \right] ]
这个内层的 sup 就是“线性函数减去范数”的上确界。数学上，它是范数的凸共轭函数在 π 处的取值。结论很干净：当 λ > 0 时，这个 sup 是有限值，并且等于：
[ \pi^\top (\hat{\xi}i - T x) + \mathbb{I}{{ |\pi|_* \le \lambda }} ]
其中 (|\cdot|_*) 是原范数的对偶范数。直观解释是：你在“模糊集合有多大”和“第二阶段对偶变量能有多锋利”之间做了一个权衡。λ 越小，允许的 π 的范数越小，也就是说对随机参数的响应越“钝”，模型越保守；λ 越大，对偶变量越“尖”，对场景的响应越灵敏。
这里需要特别说明：上面的推导假设了地面距离取 L1 范数，此时对偶范数是 L∞ 范数；如果地面距离取 L2 范数，对偶范数还是 L2 范数。两者给出的模糊集合性质差异很大，我在选型时会优先用 L1，因为得到的是线性约束，求解器友好度最高。
3.3 合并外层的 min 问题
把所有样本场景合并回外层，我们就得到了一个完整的、不再包含连续分布变量的优化问题：
[ \min_{x, \lambda \ge 0} ; c^\top x + \lambda \epsilon + \frac{1}{N} \sum_{i=1}^{N} \max_{\pi_i : W^\top \pi_i \le q, |\pi_i|_* \le \lambda} \pi_i^\top (\hat{\xi}_i - T x) ]
到这里，模型已经从“分布变量”变成了“每个样本对应一个对偶变量 π_i”的问题。这是一个有限维凸优化问题，标准求解器都能处理。
如果你不打算用线性决策近似，第二阶段变量还是“隐式”的，因为这里的 π_i 是在内层 max 中，并没有直接给出 y 的表达式。但当我们把 y 限定为 ξ 的线性函数后，可以直接代入原约束，然后对整个问题做一次统一的凸重构，这就是我后面在 Matlab 里用的实现方式。两种路线各有利弊：用对偶变量路径，变量数是 N × dim(π)，可能比较大；用线性决策路径，变量数是 dim(y0) + dim(Y)，通常更小。
4. Matlab 实现：从数据到可求解模型
4.1 生成经验分布场景
先把问题抽象成标准输入：一个随机参数向量 ξ，若干历史样本 ξ_hat_i。在 Matlab 里，我用一个简单的二维问题做演示：第一阶段决策 x 是某个容量或者订货量，第二阶段随机参数是需求 ξ = (ξ1, ξ2) 的两个组成部分，比如旺季、淡季需求。
实际写代码前，先想清楚：要不要做场景缩减？如果样本有几千个，直接带上所有场景，对偶转化后的变量规模会很大。简易模型阶段，建议先用 50~100 个样本，跑通流程后再考虑做 Wasserstein 场景缩减或者 K-means 聚类。样本太少也不行，经验分布离真实分布太远，模糊半径会变得很大，得到的结果过于保守，参考价值有限。
% 生成演示用历史需求样本 rng(42); N = 80; xi_hat = [2.5 + 0.8*randn(N,1), 1.8 + 0.6*randn(N,1)]; % 第一阶段决策变量维度 nx = 2;这一步没啥玄学，但有一个细节很多人会忽略：randn生成的正态分布尾巴太薄，如果真实场景是重尾分布，Wasserstein 半径就得放大。建议做一版重尾实验，比如把一部分样本改成对数正态或者混合正态，观察结果是否稳健。我一般会并行跑三组数据：正态、混合正态、带异常值的脏数据，这三组结果能快速暴露模型对分布假设的敏感度。
4.2 用 YALMIP 构建线性决策规则
Matlab 里做凸优化建模，我个人的首选还是 YALMIP。原因很直接：它支持把约束和变量以矩阵形式写出来，代码跟数学表达式的对应关系很清楚；而且后端换求解器只需要改一行，不需要重写模型。下面给出核心建模片段：
% yalmip 变量 x = sdpvar(nx, 1, 'full'); Y = sdpvar(ny, d_xi, 'full'); % 线性决策系数矩阵 y0 = sdpvar(ny, 1, 'full'); % 第二阶段对偶变量 pi_cell = cell(N, 1); lambda = sdpvar(1, 1); % 目标函数 obj = c'*x + lambda*epsilon + (1/N) * sum_cost;这里我把内层 max 拆到每个场景里，用pi_cell{i}表示第 i 个样本对应的对偶变量。约束包括两部分：对偶可行约束 (W^\top \pi_i \le q)，以及对偶范数盒约束 (|\pi_i|_{\infty} \le \lambda)。后者在 YALMIP 里可以直接写norm(pi_cell{i}, inf) <= lambda。
有一点要提醒：lambda也是决策变量，而且它是一个公共标量，所有场景共用。这意味着模型在优化时会找到“一视同仁”的最优模糊半径价格，这正是 Wasserstein 对偶的巧妙之处——不是给每个点单独惩罚，而是全体场景共享同一个 λ。
4.3 求解器选择与参数微调
模型建好后，求解器我推荐用sedumi、mosek或者sdpt3。小规模问题用sedumi就够了，它对稀疏结构处理得不错；如果场景数涨到上千，建议换gurobi或者mosek的商业版，速度差距非常明显。开源求解器在小规模 demo 里没问题，但跑到大规模场景时容易在数值公差上栽跟头。
还要注意数值尺度问题。Wasserstein 半径 ε 的数量级跟地面距离的选择强相关。如果你的样本量级是 10^3，而 ε 给的是 10^{-2}，数值上两者的差会在 KKT 系统里产生病态条件数。我通常会把样本标准化到 [0,1] 区间后再计算 Wasserstein 半径，这样数值表现稳定很多，也更方便解释半径的物理意义。
% 调用求解器 options = sdpsettings('solver', 'mosek', 'verbose', 0, 'dualize', 1); sol = optimize(constraints, obj, options);这里dualize=1是 YALMIP 的一个隐藏性能开关，它会把原问题交给求解器时先做一次内部对偶化，很多求解器在对偶形式下收敛更快。别小看这个设置，我遇到过同一个问题，不开这个开关要跑 40 秒，打开之后 3 秒出解。
5. 数值实验：一个简化的两阶段报童库存 Demo
5.1 问题设定
为了验证模型，我搭了一个简化版报童问题。第一阶段决策 x 是采购量，第二阶段在随机需求 ξ 实现后，可以采取两个补救动作：缺货时紧急采购（单位成本 cu），过剩时折价处理（单位成本 co）。两个成本都高于正常采购成本的一半左右，所以模型不会乱来。
目标函数是：
[ \min_x ; c^\top x + \sup_{P \in \mathcal{D}} \mathbb{E}_P [, cu \cdot (需求 - 库存)^+ + co \cdot (库存 - 需求)^+ ,] ]
第二阶段决策写成 y1（紧急补货量）、y2（折价处理量），约束是：
[ x + y_1 - y_2 = \xi, \quad y_1, y_2 \ge 0 ]
然后套入线性决策规则：y1 = y10 + Y1ξ，y2 = y20 + Y2ξ。最终整个模型变成一个只含 x、Y、λ 的凸优化问题。我在 Matlab 里跑了三组半径：ε = 0、ε = 0.1、ε = 0.3。
5.2 结果与解读
| 参数 | ε = 0 | ε = 0.1 | ε = 0.3 |
|---|---|---|---|
| 第一阶段采购量 x | 2.30 | 2.51 | 2.77 |
| 目标函数值（样本内） | 4.12 | 4.38 | 4.79 |
| out-of-sample 平均成本 | 5.31 | 4.92 | 5.02 |
| 最差 5% 场景平均成本 | 8.81 | 7.20 | 6.64 |
这个表很能说明问题。ε=0 时，也就是普通随机规划，样本内成本最低，但 out-of-sample 表现是最差的，最差 5% 场景成本高达 8.81。这就是“样本内拟合过度”的典型表现。ε=0.1 时，样本内成本只涨了 6%，但最差场景成本降了 18%，性价比很高。ε=0.3 时，样本内成本已经明显上涨，out-of-sample 反而也开始变差，因为过度保守导致采购量过大，正常场景下积压成本上升。
这张表我会建议贴在项目汇报里。它比任何算法名词都更能说明“为什么要做分布鲁棒优化”：不是为了让样本内指标变漂亮，而是为了让决策在真实环境里更稳。
5.3 out-of-sample 验证的正确姿势
不少人在验证时有个坏习惯：拿同一批数据既用来构造经验分布，又用来评估 out-of-sample。这是典型的“偷看答案”。正确做法是：
- 拿出一部分样本当作“训练数据”，构造经验分布和模糊集合；
- 再拿出一部分没进入模型的数据，作为“真实分布”的模拟，评估决策在那里实际产生的成本；
- 更严格一点，用多个随机种子重复切分，报告平均值和方差。
我在这个 Demo 里留了 30% 样本做 out-of-sample 评估，且每轮调参后都会重新切分一次。还有一件事容易被忽略：out-of-sample 评估时，第二阶段补救动作要按真实最优解来算，而不是按线性决策规则来算。因为线性决策规则只是求解时的近似，真正执行第二阶段时，你完全可以做更聪明的补救。如果你用相同的线性规则去评估，那等于把近似误差也算进了成本里，会高估 DRO 的坏处。
6. 踩过的坑与调试建议
6.1 模糊半径 ε 的量纲陷阱
这是我最开始犯过的错误。Wasserstein 距离是有量纲的：如果 ξ 是需求量，量纲是“件”，那么距离的量纲也是“件”。如果你随手把 ε 设成 0.5，对二维需求是 0.5，对两位数的需求可能也是 0.5，但相对偏差完全不同。我建议实际项目里先计算两两样本之间的平均 Wasserstein 距离，以这个为基准来定 ε。比如，样本间平均距离是 0.84，那 ε=0.1 就代表允许偏差约为样本波动幅度的 12%，这比较好解释。
6.2 对偶范数约束写错
对偶范数约束是整个对偶转化里最容易被写错的一步。地面距离是 L1 范数时，对偶范数是 L∞ 范数；地面距离是 L2 范数时，对偶范数才是 L2。如果你选了 L2 地面距离却在代码里写了norm(pi, inf) <= lambda，得到的结果会偏离真实模糊集合，而且不容易看出问题，因为求解器照样会给你一个“看起来合理”的解。调试方法是：固定 x，只优化内层 sup，用暴力搜索或者采样做交叉验证，看对偶目标值和直接采样的最坏成本是否吻合。
6.3 线性决策规则导致第二阶段不可行
这是一个很经典的坑。线性决策规则把 y(ξ) 强制写成 y0 + Yξ，虽然模型整体可行，但在某些具体 ξ 实现下，y(ξ) 可能不满足非负约束或者平衡约束。你优化出的“最优解”可能在真实执行时根本不是可行解。解决办法有两个方向：
- 在建模时加入线性决策规则的可实施性约束，也就是要求 y0 + Yξ 对所有模糊集合内的 ξ 都满足约束；
- 或者在求解后做一个后验检查，看看对每个训练样本和测试样本，y(ξ) 是否真的可行。如果不可行，就需要扩大第一阶段决策或者放弃线性决策近似。
我在 Demo 中就加了可实施性约束：对每个样本点，直接写出 y0 + Y*ξ_i >= 0，并满足平衡方程。这样虽然会稍微增加约束数量，但换来的是求解结果真正可落地。
6.4 求解器的数值容差问题
最后说一个所有优化人都绕不开的问题：数值容差。Wasserstein DRO 的对偶问题里经常出现数值很小的 λ 和数值很大的 π。如果求解器默认的原始对偶可行性容差是 1e-6，你得到的目标值可能有 1e-3 左右的误差。听起来不大，但在样本规模上千时，这些误差会叠加到决策变量上，导致第一阶段决策出现抖动。
我的习惯是：求出解之后，把 λ 固定成最优值，重新求解一次不带 λ 变量的模型，前后两个目标值做一次对比。如果差超过 0.5%，就说明数值状态有问题，需要调整求解器精度或者重新缩放数据。这个小技巧可以帮你在项目交付前筛掉大部分隐性 bug。
这个模型整体做下来，我的感受是：Wasserstein 距离本身并不难理解，难的是把对偶转化、线性决策、数值实现这三块拼在一起的时候，很容易因为某一个环节的小误差导致整体结果失真。如果你也打算在项目里尝试分布鲁棒优化，建议从最简单的两阶段线性模型开始，先把对偶代码跑顺，再加复杂度。这样才能确保每一步调整都能定位到问题来源，而不是在一个不可解释的黑盒里反复试参数。
