两阶段分布鲁棒优化：Wasserstein距离与线性决策Matlab实现-北京尧图网络

先把结论摆在前面这套模型不需要你拥有顶尖的优化理论功底但一定得把“分布不确定”这几个字放在脑子里反复摩擦。两阶段分布鲁棒优化DRO本身不算新概念可一旦引入Wasserstein距离构造模糊集再套上线性决策规则整套模型的对偶转化和Matlab实现就容易把人绕晕。我最近用一个简易版本把整条链路跑通了从建模、Wasserstein对偶转化到线性决策近似再落到Matlab代码整个过程比想象中曲折也比想象中有意思。这篇文章就把我的拆解思路和踩坑过程完整分享出来适合正在复现DRO、被对偶推导劝退、或者想快速跑通一个小算例对比效果的读者。1. 为什么是“两阶段”“Wasserstein距离”1.1 决策问题里天然存在“现在决定”和“看到结果后再调整”很多运营优化问题都可以拆成两阶段第一阶段必须在不确定参数实现之前拍板第二阶段则等到不确定参数逐渐明朗后还能做一些补救性调整。最简单的例子是产能规划第一阶段的x是提前建设的常规产能第二阶段可以针对实际需求临时加急生产、外购或者调度库存。这类问题用数学语言写出来就是一个典型的两阶段随机规划min c’x E_P[ Q(x, ξ) ]其中第一个决策x不受随机参数影响第二个阶段的值函数Q(x, ξ)内部还可以优化。之所以用这种结构是因为现实里我们很少能把所有决策拖到最后一刻才做总有一部分投入需要提前落定。这种“先决策、后观察、再调整”的结构也叫here-and-now与wait-and-see的划分。如果直接把这句话写成随机规划理论上很顺但真正动手求解时难点就全部转移到第二阶段Q(x, ξ)对ξ的依赖上。特别是当ξ的维度稍微高一点、约束复杂一点整个值函数就变成一个高维分段函数切起来很痛。这也是为什么后来大家会考虑给第二阶段的调整策略加一层结构比如线性决策规则让问题从“嵌套优化”变成一个规模更大的单层优化。1.2 只有历史样本时点估计是最容易翻车的做法刚开始做这类问题的人最自然的想法是把E_P[Q(x, ξ)]替换成样本均值也就是随机规划里常说的SAA。假设手头有N个历史样本ξ̂₁…ξ̂_N就把目标写成(1/N)ΣQ(x, ξ̂_i)。代码写起来简单跑起来也快但它有一个致命前提样本分布必须和真实分布足够接近。现实情况往往不是这样。我试过用50个历史样本训练一个订货策略同样的策略放到另一批时间段的真实数据上回测尾端成本比SAA预测的高出一大截。原因也不难理解SAA把所有概率质量都放在观测到的几个样本点上对没出现过的极端情况完全没有抵抗力。你做的决策只是在讨好历史样本而不是在应对“未来可能和过去不完全一样”这件事。于是就有了另一个极端经典鲁棒优化把ξ放进一个确定的集合里比如盒子区间要求所有约束在集合内都满足。这种方法保护性确实强但它完全不区分概率把所有可能点都当成同等重要结果就是解出来的策略往往保守得让人没法用。我需要的是介于两者之间的东西既承认历史样本提供的信息又允许真实分布和经验分布存在一定偏差。这就是分布鲁棒优化登场的时机。1.3 Wasserstein距离本质上是“搬运分布需要花多少钱”分布鲁棒优化里最关键的零件是如何构造一个包含真实分布的模糊集。一种很流行的做法是在经验分布附近画一个Wasserstein球。Wasserstein距离和KL散度这类度量最大的不同在于它描述的是把一个概率分布“搬”成另一个概率分布所需的最小代价。你可以把它想象成有两堆沙Wasserstein距离关心的是每粒沙子平均要搬多远而不是单纯比较两堆沙的形状有多接近。用Wasserstein距离构造模糊集的好处很直观如果真实分布只是相对经验分布做了一个小范围的扰动那么它们的Wasserstein距离就该很小反过来真实分布在某个局部方向出现比较大的偏移就一定会付出相应的距离代价。这个性质让它比“把所有概率都限制在支撑集内”的散度类模糊集更贴合实际也不要求真实分布必须和经验分布有完全相同的支撑。从我复现的角度看Wasserstein距离还有一个隐藏优势它对偶形式更规整可以和第二阶段优化问题做二次对偶最终把原本无穷维的概率分布优化问题转化成有限维凸优化问题。这正是标题里“对偶转化”这四个字的分量所在。做Matlab实现时对偶转化不是一个可选项而是让模型落地成线性规划或二阶锥规划的必经之路。2. 建模与对偶转化把“找最坏分布”变成可求解的优化问题2.1 两阶段DRO问题的完整数学结构带Wasserstein模糊集的两阶段DRO标准形式可以写成min_x c’x sup_{P ∈ B_ε(μ̂_N)} E_P[ Q(x, ξ) ]其中μ̂_N是经验分布B_ε是半径ε的Wasserstein球。第二阶段的值函数Q(x, ξ)本身是一个优化问题Q(x, ξ) min_y d’y s.t. Ay ≥ h - Tx - Mξ, y ≥ 0换句话说第二阶段还有自己的约束和决策变量。放到具体业务里y就是“看到实际需求ξ后能做的补救动作”比如紧急采购、加急运输、临时租用设备。这个结构与普通随机规划最大的区别是最外层的期望不是按照某个固定分布来算而是要在Wasserstein球里寻找那个让期望成本达到最大的分布。这个max对很多初学者来说第一反应是害怕因为它涉及概率分布空间上的优化。理论上这是一个无限维线性规划经典的对偶理论依然成立但需要引入测度论的细节。好消息是对于Wasserstein模糊集已经有成熟的对偶结论内部的最大化可以通过一个带标量对偶变量λ的表达式来替换最终变成对一个有限维问题求最小。2.2 对偶转化到底在做什么如果只记住一句话我会说对偶转化把“寻找最坏概率分布”从显式建模中消掉了代价是多引入一个标量变量λ和一组在每个样本点上的支撑函数上确界。其大致形态如下sup_{P ∈ B_ε(μ̂_N)} E_P[Q(x,ξ)] ≈ inf_{λ≥0} λ ε (1/N) Σ_{i1}^N sup_{ξ∈Ξ} [ Q(x,ξ) - λ d(ξ, ξ̂_i) ]这里的d(ξ, ξ̂_i)是Wasserstein距离背后的运输成本。λ可以理解成“允许分布偏离经验样本的单位价格”。当λ很小时模糊集半径ε几乎不起约束作用模型退化成接近SAA当λ变大惩罚也随之上升模型会更倾向于照顾那些离样本较远的极端点。注意上式里每一项还要再算一个sup_{ξ∈Ξ}这就把困难从“找分布”转移成了“在支撑集Ξ上找一个让方括号内取最大的ξ”。如果Ξ是连续集合这个子问题依然是个半无限规划并不好解如果Ξ是离散的代表点集那这个sup就只是一个有限枚举代码实现立刻简单一个数量级。我的简易Matlab版本正是从这一步切入的先把支撑集离散化用有限个候选场景点代表所有可能出现的ξ再在这个支撑上用线性决策规则近似第二阶段策略从而把整个模型落成一个单个线性规划。2.3 对偶后第二层问题怎么处理上面只是处理了外层“找最坏分布”的sup但Q(x,ξ)内部还有一个min。面对这种“外层分布最大、内层决策最小”的嵌套结构通常有两种处理路线。第一种是保留内层min只在每个离散支撑点上直接求解一个确定性的第二阶段线性规划把Q(x,ξ̂_j)当作一个黑盒函数返回。这种做法在Matlab里最容易实现适合验证对偶公式但由于Q是分段线性的整个外层问题会变成非光滑优化不能用简单的线性规划求解器一把梭。第二种是给内层决策y套一个线性决策规则y(ξ) Y₀ Y₁ξ并把它和第一阶段决策x一起作为待优化变量。这样做的好处是第二阶段不再嵌套求解Q也变成一个关于x和Y的显式函数整条模型就能整体丢给一个大规模线性规划求解器。这也是标题里“线性决策”的技术定位牺牲一点第二阶段的表示能力换取计算上的便利和理论上的可处理性。在标准两阶段DRO论文里这两层通常都要对偶展开推导量很大。但如果你一开始就选线性决策规则就可以绕开内部min的对偶直接将近似的决策函数代进目标函数让所有约束在支撑集上逐点成立。我的体验是对刚入门的人而言这个路线是把DRO从“看不懂的论文”变成“能跑的代码”的最佳捷径。3. 线性决策规则为第二阶段的响应函数加形状约束3.1 为什么不能直接把y当作独立变量理论上第二阶段最优决策y应该在看到实测值ξ之后随ξ变化而自动取到当前最优值。但在分布式鲁棒框架里有个麻烦最坏的那个P到底长什么样我们并没有显式知道。如果你把每个场景下的y都定义为独立变量等于把“分布鲁棒”的信息结构破坏了甚至可能让目标函数里出现可以利用场景概率权重的投机行为。线性决策规则的做法是主动限制第二阶段的决策形式y(ξ) Y₀ Y₁ (ξ - ξ₀)其中ξ₀可以是样本均值或某个参考点Y₀和Y₁变成新的优化变量。这是一个仿射函数不是把所有场景独立处理而是让决策随着ξ平滑变化。这样做的核心价值在于第二阶段变量不再依赖某个特定分布而是一个对所有可能ξ都成立的函数因而能整体放进鲁棒优化框架里考察。用工程上的话说你相当于给第二阶段的控制器指定了一个线性反馈律。第一阶段拍板的是“常规项Y₀和反馈增益Y₁”一旦真正观测到ξy就按照Y₀ Y₁ξ自动执行。这种“先定反馈律再按反馈执行”的思路在鲁棒模型预测控制里非常常见。3.2 代进模型后问题如何降维把y(ξ)代入第二阶段的表达式原本的min_y d’y会变成d’Y₀ d’Y₁ξ这是一个关于ξ的仿射函数。于是原问题的目标里第一阶段成本、第二阶段仿射成本都能很自然地写在一起。对于约束条件比如B y ≥ h - Tx - Mξ只需要把每个代表支撑点ξ_j代进去验证B (Y₀ Y₁ξ_j) ≥ h - Tx - Mξ_j, ∀ξ_j ∈ Ξ这实际上是把“所有可能的ξ都必须可行”这个半无限约束转换成有限个线性约束。只要支撑集Ξ选得足够能代表真实情况线性决策就能在这个近似意义下达到较好的性能。我的经验是这个近似对成本函数相对平滑的问题效果相当好尤其是库存补货、产能分配这类场景第二阶段最优决策本来就很接近需求量的线性函数。3.3 线性决策也不是万能药用久了之后你会发现线性决策最大的猫腻在于它掩盖了第二阶段决策的非线性。如果第二阶段本身带有0-1变量、开机停机逻辑、阶梯价格或者约束中存在max/min这种不可微项那么纯仿射决策会偏保守甚至可能根本没有可行解。一个从实践里来的判断标准是先小规模精确求解第二阶段LP画出不同ξ下的最优y曲线如果y随ξ接近线性再用线性决策规则就非常安全如果曲线明显折线、台阶状那就需要用分段线性决策增加表达力。另外一个容易被忽略的点是线性决策规则与Wasserstein模糊集结合时约束通常是在支撑集上逐点验证的。如果支撑集选得太稀疏那些没被采样到的关键位置可能刚好违反约束最终结果就会翻车。所以我会在支撑集里额外加入方向性极端点让线性函数在这些极端位置也被拦住这和经典鲁棒优化里的“不确定集合顶点枚举”思路异曲同工。4. Matlab代码实现一个可以跑通的简化示例4.1 选一个能说明问题的业务场景为了不让代码淹没在抽象记号里我选了一个很小但结构完整的例子某企业需要在需求不确定的环境下决策常规产能x。第一阶段每建设一单位常规产能的成本是c实际需求ξ发生后如果常规产能不够可以用成本较高的临时产能补足缺口这部分由第二阶段的线性决策y(ξ)决定。这个场景的两阶段结构一清二楚第一阶段决定x第二阶段根据实际需求决定y并保证x y ≥ ξ。作为教学版我把需求ξ的支撑集简化为从历史样本生成的离散代表点Wasserstein距离用一维绝对值距离模糊集半径记为ε。这样对偶转化后整个模型可以用线性规划求解。4.2 核心代码骨架与公式对照下面这段YALMIP代码是模型最核心的部分。我故意不把代码写得像工程仓库那样庞杂希望你能从每一行看到它对应的公式% 基础参数 N 40; % 历史样本数 / 代表点数量 xi randn(1, N) * 2 10; % 生成历史需求样本可换成自己的数据 c 3; % 第一阶段单位产能成本 q 8; % 第二阶段临时产能单位成本 epsRadius 0.3; % Wasserstein球半径后续可扫描 % 计算代表点两两之间的运输距离一维问题直接用绝对值 D abs(xi - xi); % D(j, i) 表示从候选点j搬到样本i的代价 % 决策变量 x sdpvar(1, 1); % 第一阶段常规产能 y0 sdpvar(1, 1); % 线性决策: y(xi) y0 y1 * xi y1 sdpvar(1, 1); lambda sdpvar(1, 1); % Wasserstein距离的对偶变量 t sdpvar(N, 1); % 每个样本点的上确界辅助变量 % 约束集合 Constraints [x 0, lambda 0]; % 线性决策需要在所有支撑点上满足非负和供需约束 for j 1:N y_j y0 y1 * xi(j); Constraints [Constraints, y_j 0]; Constraints [Constraints, x y_j xi(j)]; end % 对偶转化后的核心约束t(i) sup_j [cost_j - lambda * D(j,i)] for i 1:N for j 1:N cost_j q * (y0 y1 * xi(j)); % 支撑点j对应的二阶段成本 Constraints [Constraints, t(i) cost_j - lambda * D(j, i)]; end end % 目标函数一阶段成本对偶化后的Wasserstein惩罚 objective c * x epsRadius * lambda (1 / N) * sum(t); % 求解 optimize(Constraints, objective); % 提取结果 x_opt value(x); y0_opt value(y0); y1_opt value(y1); lambda_opt value(lambda); fprintf(最优常规产能 x %.4f\n, x_opt); fprintf(二阶段线性决策 y(xi) %.4f %.4f * xi\n, y0_opt, y1_opt);这段代码对应的公式是min_{x,Y,λ≥0} c x ε λ (1/N) Σ_i t_i s.t. t_i ≥ q y(ξ_j) - λ d(ξ_j, ξ̂_i), ∀i,j如果你把支撑点直接取成历史样本自身代码里的j循环和i循环会用到同一组xi但这不意味着它们在逻辑上完全重复。i对应的是“经验分布的中心样本”而j对应的是“最坏分布可能把质量转移到的目标位置”。理解这个区别后你后面把代码扩展到二维需求、多维参数时

两阶段分布鲁棒优化：Wasserstein距离与线性决策Matlab实现

相关推荐

opencode 开源 AI 编程助手：从安装配置到实战技巧全解析

AI Skills实战：在腾讯云上让Agent按流程稳定干活

大模型推理服务全流程拆解：从显存规划到性能优化实战

更多资讯

前端html2canvas结合print-js打印页面，并解决多余空白页问题

基于Django的StockMing股票预测系统：从ORM建模到ECharts可视化全解析

Puppeteer 请求拦截完全指南：从基础拦截到 Cooperative Intercept Mode 协作式拦截

从JSON数据清洗到数据库导入：背单词词库建设的实战记录

并发编程（二）—synchronized和volatile

C# 未排序数组中第 k 个最小/最大元素 | 最坏情况下的线性时间

今日更新

能源系统对接碳交易的路径：从碳核算到碳资产运营

智能硬件数据驱动运营：订单日记与预测性维护实践

逆变器下垂控制与虚拟同步机对比：Simulink仿真建模与参数整定

本周精选

超人会飞不算本事：系统稳定依赖清晰规则与边界设计

超人VS蜘蛛侠：拆解超级IP的影响力与传播方法论

基于CNN的调制信号识别：MATLAB实现时频图分类实战

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？