做分布鲁棒优化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_{i1}^{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_{i1}^{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);这里dualize1是 YALMIP 的一个隐藏性能开关它会把原问题交给求解器时先做一次内部对偶化很多求解器在对偶形式下收敛更快。别小看这个设置我遇到过同一个问题不开这个开关要跑 40 秒打开之后 3 秒出解。5. 数值实验一个简化的两阶段报童库存 Demo5.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第一阶段采购量 x2.302.512.77目标函数值样本内4.124.384.79out-of-sample 平均成本5.314.925.02最差 5% 场景平均成本8.817.206.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 距离本身并不难理解难的是把对偶转化、线性决策、数值实现这三块拼在一起的时候很容易因为某一个环节的小误差导致整体结果失真。如果你也打算在项目里尝试分布鲁棒优化建议从最简单的两阶段线性模型开始先把对偶代码跑顺再加复杂度。这样才能确保每一步调整都能定位到问题来源而不是在一个不可解释的黑盒里反复试参数。