简介：这是一份面向机器学习、优化算法研究者与电力系统从业者的分布鲁棒优化实践代码包，聚焦以Wasserstein距离构造模糊集来处理风光出力等不确定性。代码完整演示了Wasserstein球半径确定、样本集标准化、复杂分布鲁棒约束的数学转化与求解思路，能将理论公式落地为可执行的Python实现，适合需要快速上手该方法的工程师与科研人员。资源共4个文件，包体仅6KB，包括核心Python脚本（wasserstein_dro.py）、txt说明文档、inscode交互运行页面及html可视化文件，结构精炼便于对照理解。目前已有113人学习，可用于电力需求预测、储能管理、新能源并网等场景，也可迁移至金融风控、供应链管理等相关分布鲁棒优化问题。

1. 把“拍脑袋的置信区间”换成带保护带的决策：Wasserstein距离与分布鲁棒优化解决什么

在做电力系统调度和机器学习模型稳健性这两类工作时，我反复遇到同一个尴尬：历史样本足够多，但“明天”不等于“昨天”。随机优化把历史分布当成真实分布，一旦出现极端风电、突发负荷就翻车；纯鲁棒优化又过于保守，成本高到没法落地。分布鲁棒优化（DRO）用Wasserstein距离在经验分布周围画一个“球”——未来真实分布只要落在球内，优化结果就要顶得住其中最坏的情况。它适合做机组组合、新能源消纳、以及对抗样本训练，核心是用优化算法把分布不确定性显式嵌进目标函数，而不只是加一个扰动项。本文以电力系统场景为主线，把Wasserstein球的构造、对偶求解与参数踩坑讲透。

2. Wasserstein球：把“分布不确定”变成可度量的数学对象

2.1 为什么是Wasserstein而不是KL散度

Wasserstein距离的定义来自最优传输。对两个分布μ和ν，给定地面距离d(ξ1,ξ2)，p阶Wasserstein距离写作

W_p(μ,ν) = ( inf_{π∈Π(μ,ν)} ∫_{Ξ×Ξ} d(ξ1,ξ2)^p dπ(ξ1,ξ2) )^(1/p)

其中Π是μ和ν的所有耦合。p=1且d取范数时，它等于把μ的质量搬到ν所需要的最小运输成本。工程上可以这样理解：一个场景的分布偏移主要是“平移”和“形变”造成的，而不是概率质量凭空增减。

KL散度的定义D_KL(μ‖ν)=∫log(dμ/dν)dμ要求μ关于ν绝对连续。用历史数据的经验分布作为名义分布时，μ是N个离散点的均匀分布，ν是连续分布，这时KL会得出无穷大，无法作为距离。而Wasserstein距离对离散-连续、离散-离散都给出有限值，不需要支撑重叠。

φ散度和总变差同样存在这个问题：它们更多度量“在同一支撑上概率密度形状差多少”，而Wasserstein度量“质量移动了多少”。对电力系统这种场景偏移频繁的问题（风电爬坡时刻前后移、负荷峰谷平移），Wasserstein的语义更贴切。还有一个计算上的优势：当名义分布是经验分布时，Wasserstein球的对偶有显式结构，这是KL和φ散度没有的。所以项目代码里几乎所有模糊集都以Wasserstein距离为度量，这不是偶然。

另一个容易被忽略的点是W1和W2的选择。W1用线性成本，对离群场景的惩罚相对温和；W2对大偏移给更大惩罚，对偶后的约束从线性变成二次锥。论文里常用W1，因为把问题压成LP后求解器成熟、调试方便。我在复现时会先用W1把整个流程跑通，再用W2做敏感性对比，两种距离下的ε语义完全不同，不能拿同一个数值直接套。

2.2 模糊集半径ε：量纲、退化和样本量缩放

模糊集的定义是

B_ε(μ̂_N) = { Q ∈ P(Ξ) | W(Q, μ̂_N) ≤ ε }

其中μ̂_N是经验分布，(1/N)Σ δ_{ξ_i}。ε是半径，决定我们愿意容忍分布偏移多少。语义上，ε与样本空间的距离单位一致，不是无纲量的百分比。如果样本是风电功率MW，ε=0.5意味着允许候选分布相对经验分布有0.5MW的“运输成本”偏移。如果代码在标准化之前就把ε设定成某个“看起来合理”的小数，然后标准化之后又忘记调整，几乎必然踩坑。

ε=0时，模糊集退化为单点集，DRO退化为样本平均逼近（SAA）；ε→∞时，模糊集包含所有与μ̂_N距离有限的分布，DRO退化为最坏情况鲁棒优化。有意义的ε总是在两者之间。工程上，ε应该随样本量的增加而缩小：经验分布对真实分布的Wasserstein收敛速率大致是N^{-1/d}量级，d是不确定性的维度。N越大、支撑点越密集，球半径可以越小。我通常用两类方法标定ε：一是用历史数据做block bootstrap，估计W(Q,μ̂_N)的经验分位数；二是把ε当作超参数，在验证集上做网格搜索。两类方法结果差别很大，但至少能给出一个可解释的范围。

还有一个关于标准化的细节需要单独提出来：数据标准化不是免费的。先对特征做z-score再画Wasserstein球，等价于把各向同性的球变成椭球。标准化本身没问题，但必须在统一的样本集上计算均值和方差，否则会泄漏未来信息。我习惯在写完数据处理流水线后，顺手打印一组标准化前后的W1距离对比，确认半径语义没有漂移。

2.3 先用一个最小实验确认半径量级

在开始求解项目代码里的大规模模型之前，我会用一个一维合成样本验证ε和分布偏移量级的关系。下面的代码把样本随机分半，用一维经验分布之间的W1距离近似估计“分布自然偏移”的尺度：

import numpy as np np.random.seed(7) N = 1000 # 风电出力历史样本，单位 MW xi = np.random.gamma(2.0, 15.0, N) xi = np.clip(xi, 0, 80) def w1_approx_1d(a, b): # 一维经验分布之间的 W1 近似：排序后逐点配对 # 严格 W1 应对 CDF 积分，这里用于量级估算 n = min(len(a), len(b)) return np.abs(np.sort(a)[:n] - np.sort(b)[:n]).mean() # 随机分半，重复 100 次，得到经验分布偏移量的分布 offsets = [] for _ in range(100): idx = np.random.permutation(N) train = xi[idx[:N // 2]] valid = xi[idx[N // 2:]] offsets.append(w1_approx_1d(train, valid)) offsets = np.array(offsets) print("offset 中位数: {:.2f} MW".format(np.median(offsets))) print("offset 90% 分位: {:.2f} MW".format(np.quantile(offsets, 0.9)))

这段代码的核心是先在数据上估算“经验分布与经验分布之间”的偏移量级。如果分半样本之间的W1中位数是2MW，那么ε取1MW可能太紧，取10MW又可能过度保守。真实复现时要调用真正的W1求解器或传输距离库，但建模初期用这种快速近似能省掉大量排错时间。注意这里的w1_approx_1d用排序配对方式近似，只适用于一维且长度相近的样本，多维场景不能用这个函数，要回到运输问题的线性规划。

3. 从min-max到有限维线性规划：强对偶与对偶变量的落地

3.1 写出原问题，再把它拉回有限维

DRO原问题写作

min_{x∈X} sup_{Q∈B_ε(μ̂_N)} E_Q[f(x,ξ)]

表面上是无穷维问题，因为Q是完备概率空间里的分布。当名义分布是经验分布μ̂_N=(1/N)Σδ_{ξ_i}时，可以借Wasserstein球的强对偶把内层sup转换掉。标准结论是：若地面距离d是范数，f(x,·)对每个x是上半连续且下有界，支撑集Ξ是紧集，那么内层sup等价于

min_{λ≥0} [ λε + (1/N) Σ_{i=1}^N sup_{ξ∈Ξ} ( f(x,ξ) - λ d(ξ, ξ_i) ) ]

其中λ是球约束的解耦乘子。这个形式把“Q到μ̂_N距离不超过ε”转成了两层结构：外层是λε罚项，内层是对每个历史样本ξ_i计算它在支撑集上的最坏偏差。λ是优化变量而不是超参数，求解过程中会被自动拉到合适位置。ε固定时，λ越大，相当于对分布偏移的容忍越小，模型越保守。

这个转换的价值在于：原来的变量是无穷维分布，现在只需要在有限个历史样本上加辅助变量。项目代码里最能出问题的一步就在这里——对偶引理成立需要紧支撑，而很多实现默认ξ落在全空间上，导致内层sup在目标函数线性增长时无上界，求解器直接报unbounded。

3.2 三种典型的f(x,ξ)结构及对应的约束形式

线性形式最常用。电力系统调度里，购电成本、弃风惩罚、失负荷惩罚大多随预测误差线性增长，所以f取关于ξ的线性函数。内层的sup_{ξ} { a_x^T ξ - λ‖ξ - ξ_i‖ }在支撑集是凸多面体时可转成顶点上的线性函数加范数锥约束，整体保持LP结构。

分段线性形式出现在阶梯报价、分段补偿、以及ReLU形状的切负荷惩罚里。对偶后每一段都需要一个辅助变量，约束数量约为N×(K+1)，K是分段数。如果K很大，存储和求解成本会明显上升，常见做法是把分段函数改写成线性目标加额外整数变量，或者用外逼近法逐次加约束。

二次形式通常配合2阶Wasserstein距离，在均值-方差类风险度量里很常见。它的闭式解很漂亮，但求解器从LP变成锥规划，大规模场景下比线性形式吃力得多。我一般先按线性把流程跑通，验证数据流和参数标定没问题后，再根据业务需求升级到分段线性或二次结构。

3.3 把一个线性DRO写成求解器模板

下面的代码是线性DRO对偶形式的通用求解模板，用支持集上的网格近似内层sup：

import cvxpy as cp import numpy as np np.random.seed(0) N = 100 # 历史样本数 K = 20 # 支撑集离散点数 xi_hist = np.random.randn(N, 2) # 二维不确定性：风电和负荷偏差 xi_grid = xi_hist[:K] # 支撑集候选点，工程上用聚类代表场景 # 决策变量：机组出力调整量 x = cp.Variable(2) # 对偶乘子：对应 Wasserstein 球约束 lam = cp.Variable(nonneg=True) # 每个历史样本的最坏场景偏差 S = cp.Variable(N, nonneg=True) # 模糊集半径，先用 Parameter 方便后续网格搜索 eps = cp.Parameter(nonneg=True, value=0.5) # 线性费用：f(x, xi) = c^T x + d^T xi c = np.array([1.0, 1.0]) d = np.array([0.8, -0.5]) constraints = [] for i in range(N): for k in range(K): dist = np.linalg.norm(xi_grid[k] - xi_hist[i]) constraints.append(S[i] >= d @ xi_grid[k] - lam * dist) objective = eps * lam + (1 / N) * cp.sum(S) + c @ x prob = cp.Problem(cp.Minimize(objective), constraints) prob.solve() print("最优成本: {:.3f}".format(prob.value)) print("对偶乘子 lambda: {:.3f}".format(lam.value))

逻辑说明：S[i]是第i个历史样本周围最坏场景下的线性费用，约束对每个网格点k都要成立，才能保证sup被覆盖。dist是地面距离，它直接决定分布偏移的惩罚尺度——如果dist算错，整个切割约束的几何意义就错了。eps被定义成Parameter而不是普通常数，这样之后做ε网格扫描时不用重建整个问题，只需循环赋值再solve。

参数说明：N是历史样本数，K是支撑集候选数，工程上可以取全部历史场景或聚类后的代表场景。c和d是费用系数，d对应随机参数在目标里的敏感程度，这是“预测误差越大会带来多大经济损失”的直接体现。这个模板的约束数量是N×K=2000，小型问题可以直接求解；到N=5000、K=50时就到25万条约束，必须换算法，我在第4章末尾展开讲。

4. 电力系统场景落地：把风电出力不确定性嵌进调度决策

4.1 从点预测到“分布鲁棒安全域”

风电出力预测误差不是固定的单峰分布。气温、气压、风速模式变化会让预测误差曲线整体平移或尾部变形，SAA在这里经常给出“过于自信”的机组启停方案。DRO的做法：预测误差的历史样本构成经验分布，以它为中心画Wasserstein球，调度决策在球内最坏分布下仍然要可行。机组组合里的失负荷概率约束（比如99%可靠性）被改写成对一个期望惩罚的约束，而这个期望是对球内的最坏分布取的。

这一步不是简单地把预测区间放宽。传统随机优化的置信区间只覆盖“参数取值的范围”，DRO覆盖的是“分布整体的变化范围”。所以同样的调度成本预算下，DRO给的备用容量更平滑，不会像逐点鲁棒那样同时把每组出力都钉死在最极端场景上。

4.2 两阶段DRO机组组合：变量与约束结构

两阶段结构里，第一阶段变量在观测到风电实际出力之前就要决定，第二阶段变量在每一组风电场景下分别调整。用DRO时，第二阶段约束按场景展开，目标里的期望项被对偶后替换成λε加上对每个历史样本的S_i。用一个单时段简化模型演示核心的组装方式：

import cvxpy as cp import numpy as np np.random.seed(1) N = 200 # 历史场景数 K = 21 # 支撑集网格数 load = 100.0 # 负荷，MW p_max = 120.0 # 机组容量 unit_cost = 25.0 # 发电边际成本 penalty_loss = 500.0 # 切负荷惩罚 # 历史风电预测误差（单位 MW，有正有负） delta_hist = np.random.normal(0, 10, N) # 支撑集网格：从 -40MW 到 40MW delta_grid = np.linspace(-40, 40, K) # 第一阶段：机组出力 p = cp.Variable(nonneg=True) # 第二阶段：每个场景的切负荷 l = cp.Variable(N, nonneg=True) # 对偶层变量 lam_dual = cp.Variable(nonneg=True) S = cp.Variable(N, nonneg=True) eps = 1.5 # 模糊集半径，MW constraints = [] for i in range(N): # 场景功率平衡：机组出力 + 风电实际出力 = 负荷 - 切负荷 constraints.append(p + (50.0 + delta_hist[i]) == load - l[i]) # Wasserstein 切割：对所有支撑点 k，S[i] 要覆盖该点下的最坏损失 for k in range(K): loss_expr = penalty_loss * cp.pos(load - p - 50.0 - delta_grid[k]) constraints.append(S[i] >= loss_expr - lam_dual * abs(delta_grid[k] - delta_hist[i])) objective = unit_cost * p + eps * lam_dual + (1 / N) * cp.sum(S) prob = cp.Problem(cp.Minimize(objective), constraints) prob.solve() print("最优出力 p = {:.1f} MW".format(p.value)) print("lambda = {:.2f}".format(lam_dual.value))

代码里的cp.pos是ReLU函数，切负荷量只有在负荷大于机组出力加风电出力时才为正。S[i]必须对每个网格点k都成立，它代表的是“在第i个样本周围，最坏场景下可能产生的切负荷损失”。eps=1.5的单位是MW，对应第2章里讲过的量纲问题。多时段扩展时，p变成(G,T)矩阵，切割约束按时段分别生成，目标里累加S矩阵的和。

4.3 求解规模的现实反馈

对偶后的约束数量大致是N×K乘以时段数。一个T=24、N=2000、K=30的两阶段机组组合，切割约束超过百万条，直接用通用LP求解器会非常吃力。常见做法是Benders分解：第一阶段整数变量固定后，第二阶段DRO子问题是一个纯LP，每次求解返回一个最优性切割给主问题，主问题再切一块可行域。对于支撑集过大导致的切割爆炸，我一般先把历史场景做k-means聚类，把K压到50以内，聚类中心作为支撑集候选点，这样误差可控，求解规模也下来一个数量级。项目代码里如果直接对原始N×T场景枚举，大概率在小规模Demo上能跑，放大到真实系统就卡死，这也是拿到代码后第一个需要改造的地方。

5. 避坑手册：把DRO从“论文能跑”变成“工程能用”

5.1 ε不是越大越好：调大半径容易把模型调成“破产级保守”

现象：ε取到样本标准差的5倍以上，目标成本迅速逼近最坏情况鲁棒优化的水平，机组组合结果几乎不随场景变化，备用容量高到经济上不可行。

原因：ε描述的是分布偏移的容忍范围，不是安全系数的代名词。ε增大后，模糊集里纳入大量严重偏移的分布，最坏期望成本主要被这些极端尾部主导，调度方案被迫为完全不合理的大偏差兜底。

解决：用第2.3节的分半W1估算先定量纲，再用bootstrap取90%分位数作为ε的参考上限。经验法则上，ε与N的-1/d次方同量级，N增大时ε要同步缩小，不能抱死一个固定值。

5.2 W1和W2混用，对偶结构变了还不知道

现象：代码里只改了一个norm参数，结果求解时间暴涨、数值稳定性变差，甚至结果出现明显偏差。

原因：W1的切割约束是线性的，W2的地面距离平方之后，切割里的距离项变成凸二次，问题从LP变成二阶锥规划。很多实现没有区分这两种结构，把W1的模板直接套到W2上，约束语义已经完全不是原问题。

解决：建模文档里显式写明用的是W1还是W2，地面距离的范数阶数单独用常量定义。跑通W1之后再升级W2，不要在同一个代码文件里来回改参数。

5.3 支撑集不紧，内层sup直接无界

现象：求解器返回unbounded或infeasible，检查约束时发现S[i]对应的那一组切割约束永远无法满足。

原因：支撑集Ξ没有加物理边界。风电出力理论上限、负荷上下限、爬坡约束都没有写进支撑集定义，线性目标函数在无界区域上sup趋向无穷，对偶问题自然无解。

解决：给ξ加箱式边界，或者用截断距离替代完整范数。常见做法是ξ∈[ξ_low, ξ_high]，网格点从边界内采样。加了边界之后，对偶问题的紧性条件才成立。

5.4 把“分布鲁棒”实现成“逐点鲁棒”

现象：回测结果比SAA好，但比传统区间鲁棒好不了太多，成本从第2章开始就一路偏高。

原因：实现时把sup_Q E_Q[f]错误地实现了max_i f(x, ξ_i)，也就是把场景最坏值直接挑出来，而不是让分布Q在Wasserstein球内变化。这两种做法的结果差异巨大：逐点鲁棒要求每个场景都可接受，DRO只要求最坏分布下期望可接受。

解决：检查目标函数里是否还有ε和λ的结构。真正的DRO目标在ε=0时退化为SAA，在ε较大时才是接近鲁棒的；如果ε=0时结果就和逐点鲁棒一样，说明切割约束写错了，需要回到第3.3节的模板逐条核对。

5.5 标准化时偷看了未来，导致离线回测失真

现象：离线测试集上表现优异，上线之后连续几天偏差放大，误差统计和离线完全对不上。

原因：标准化用的均值和方差是在全样本上算的，包括验证集和测试集，相当于把未来信息掺进了特征变换。Wasserstein距离对尺度敏感，标准化参数的错误会被放大进球半径语义。

解决：只用训练集计算mean和std，训练集、验证集、测试集分别用同一组参数变换。并且在代码里把标准化参数缓存下来，部署时加载同一份参数，避免在线预测时重新拟合分布统计量。

6. 上线前的一步验证：用合成分布偏移做回测

在把DRO模型放进真实调度流程之前，我会做一个合成验证，专门检验“当真实分布真的偏移时，DRO是否比SAA更稳”。做法是：用一个分布生成训练样本，再用另一个偏移后的分布生成测试样本，看两种方法在测试样本上的表现差距。下面这段代码是评估函数的核心部分：

import numpy as np np.random.seed(3) N_train = 500 N_test = 2000 # 训练分布：中心在 50MW 的偏态分布 train_wind = np.random.beta(2, 5, size=N_train) * 100 # 测试分布：中心偏移 + 右尾加重 test_wind = np.random.beta(2, 4, size=N_test) * 100 + 5 def simulate_out_of_sample(decision, wind_sample, load=100.0, penalty=500.0): shortages = np.maximum(0.0, load - decision - wind_sample) return penalty * shortages.mean()

分别用SAA和DRO求两个决策：SAA决策用训练样本直接优化，DRO决策用第4章的模型在ε=1.5下求解，然后都放到test_wind上评估。典型结果趋势如下：

解读这个表：ε=0的SAA在训练分布上最便宜，测试分布一偏移就崩；ε=1.5时训练成本略升，但测试损失大幅下降；ε=4.0时训练成本明显偏高，测试损失却只降了一点点。这说明存在一个合理ε区间，过低不够稳，过高过度保守。那张风险-成本曲线就是上线与否的决策依据。

从那以后，我每次调ε都强制先画一遍风险-成本曲线，再谈能不能上线。这个习惯帮我避开了至少三次“论文指标漂亮、现场一跑就崩”的尴尬。希望帮到你。

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