简介这是一份面向机器学习、优化算法研究者与电力系统从业者的分布鲁棒优化实践代码包聚焦以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)其中Π是μ和ν的所有耦合。p1且d取范数时它等于把μ的质量搬到ν所需要的最小运输成本。工程上可以这样理解一个场景的分布偏移主要是“平移”和“形变”造成的而不是概率质量凭空增减。KL散度的定义D_KL(μ‖ν)∫log(dμ/dν)dμ要求μ关于ν绝对连续。用历史数据的经验分布作为名义分布时μ是N个离散点的均匀分布ν是连续分布这时KL会得出无穷大无法作为距离。而Wasserstein距离对离散-连续、离散-离散都给出有限值不需要支撑重叠。φ散度和总变差同样存在这个问题它们更多度量“在同一支撑上概率密度形状差多少”而Wasserstein度量“质量移动了多少”。对电力系统这种场景偏移频繁的问题风电爬坡时刻前后移、负荷峰谷平移Wasserstein的语义更贴切。还有一个计算上的优势当名义分布是经验分布时Wasserstein球的对偶有显式结构这是KL和φ散度没有的。所以项目代码里几乎所有模糊集都以Wasserstein距离为度量这不是偶然。度量支撑不重叠时对场景偏移的惩罚经验分布球的对偶友好度Wasserstein距离有限按距离线性W1强对偶可转有限维LPKL散度无穷大指数级、易爆表名义分布需有密度总变差/φ散度有限但度量整体差异饱和快忽略小偏移约束结构不直观另一个容易被忽略的点是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) Σ_{i1}^N sup_{ξ∈Ξ} ( f(x,ξ) - λ d(ξ, ξ_i) ) ]其中λ是球约束的解耦乘子。这个形式把“Q到μ̂_N距离不超过ε”转成了两层结构外层是λε罚项内层是对每个历史样本ξ_i计算它在支撑集上的最坏偏差。λ是优化变量而不是超参数求解过程中会被自动拉到合适位置。ε固定时λ越大相当于对分布偏移的容忍越小模型越保守。这个转换的价值在于原来的变量是无穷维分布现在只需要在有限个历史样本上加辅助变量。项目代码里最能出问题的一步就在这里——对偶引理成立需要紧支撑而很多实现默认ξ落在全空间上导致内层sup在目标函数线性增长时无上界求解器直接报unbounded。3.2 三种典型的f(x,ξ)结构及对应的约束形式f(x,ξ)结构内层sup的显式结果最终问题类型线性f a_x^T ξ b_xsup落在支撑集顶点配合对偶范数线性规划分段线性f max_k (a_k^T ξ b_k)每段取凸包引入辅助变量线性/MILP凸二次f ξ^T A_x ξ 线性项用二次共轭与2-范数闭式二阶锥规划线性形式最常用。电力系统调度里购电成本、弃风惩罚、失负荷惩罚大多随预测误差线性增长所以f取关于ξ的线性函数。内层的sup_{ξ} { a_x^T ξ - λ‖ξ - ξ_i‖ }在支撑集是凸多面体时可转成顶点上的线性函数加范数锥约束整体保持LP结构。分段线性形式出现在阶梯报价、分段补偿、以及ReLU形状的切负荷惩罚里。对偶后每一段都需要一个辅助变量约束数量约为N×(K1)K是分段数。如果K很大存储和求解成本会明显上升常见做法是把分段函数改写成线性目标加额外整数变量或者用外逼近法逐次加约束。二次形式通常配合2阶Wasserstein距离在均值-方差类风险度量里很常见。它的闭式解很漂亮但求解器从LP变成锥规划大规模场景下比线性形式吃力得多。我一般先按线性把流程跑通验证数据流和参数标定没问题后再根据业务需求升级到分段线性或二次结构。3.3 把一个线性DRO写成求解器模板下面的代码是线性DRO对偶形式的通用求解模板用支持集上的网格近似内层supimport 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(nonnegTrue) # 每个历史样本的最坏场景偏差 S cp.Variable(N, nonnegTrue) # 模糊集半径先用 Parameter 方便后续网格搜索 eps cp.Parameter(nonnegTrue, value0.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×K2000小型问题可以直接求解到N5000、K50时就到25万条约束必须换算法我在第4章末尾展开讲。4. 电力系统场景落地把风电出力不确定性嵌进调度决策4.1 从点预测到“分布鲁棒安全域”风电出力预测误差不是固定的单峰分布。气温、气压、风速模式变化会让预测误差曲线整体平移或尾部变形SAA在这里经常给出“过于自信”的机组启停方案。DRO的做法预测误差的历史样本构成经验分布以它为中心画Wasserstein球调度决策在球内最坏分布下仍然要可行。机组组合里的失负荷概率约束比如99%可靠性被改写成对一个期望惩罚的约束而这个期望是对球内的最坏分布取的。这一步不是简单地把预测区间放宽。传统随机优化的置信区间只覆盖“参数取值的范围”DRO覆盖的是“分布整体的变化范围”。所以同样的调度成本预算下DRO给的备用容量更平滑不会像逐点鲁棒那样同时把每组出力都钉死在最极端场景上。4.2 两阶段DRO机组组合变量与约束结构阶段决策内容变量类型约束对象第一阶段机组启停、备用容量二进制/整数启停逻辑、最小启停时间第二阶段出力调整、弃风量、切负荷量连续功率平衡、爬坡、线路对偶层λ、S_i连续非负每个历史样本的Wasserstein切割两阶段结构里第一阶段变量在观测到风电实际出力之前就要决定第二阶段变量在每一组风电场景下分别调整。用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(nonnegTrue) # 第二阶段每个场景的切负荷 l cp.Variable(N, nonnegTrue) # 对偶层变量 lam_dual cp.Variable(nonnegTrue) S cp.Variable(N, nonnegTrue) eps 1.5 # 模糊集半径MW constraints [] for i in range(N): # 场景功率平衡机组出力 风电实际出力 负荷 - 切负荷 constraints.append(p (50.0 delta_hist[i]) load - l[i]) # Wasserstein 切割对所有支撑点 kS[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个样本周围最坏场景下可能产生的切负荷损失”。eps1.5的单位是MW对应第2章里讲过的量纲问题。多时段扩展时p变成(G,T)矩阵切割约束按时段分别生成目标里累加S矩阵的和。4.3 求解规模的现实反馈对偶后的约束数量大致是N×K乘以时段数。一个T24、N2000、K30的两阶段机组组合切割约束超过百万条直接用通用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, sizeN_train) * 100 # 测试分布中心偏移 右尾加重 test_wind np.random.beta(2, 4, sizeN_test) * 100 5 def simulate_out_of_sample(decision, wind_sample, load100.0, penalty500.0): shortages np.maximum(0.0, load - decision - wind_sample) return penalty * shortages.mean()分别用SAA和DRO求两个决策SAA决策用训练样本直接优化DRO决策用第4章的模型在ε1.5下求解然后都放到test_wind上评估。典型结果趋势如下ε训练损失测试损失最坏损失(测试)0SAA8.241.789.30.59.124.655.21.511.319.841.04.018.921.438.6解读这个表ε0的SAA在训练分布上最便宜测试分布一偏移就崩ε1.5时训练成本略升但测试损失大幅下降ε4.0时训练成本明显偏高测试损失却只降了一点点。这说明存在一个合理ε区间过低不够稳过高过度保守。那张风险-成本曲线就是上线与否的决策依据。从那以后我每次调ε都强制先画一遍风险-成本曲线再谈能不能上线。这个习惯帮我避开了至少三次“论文指标漂亮、现场一跑就崩”的尴尬。希望帮到你。本文还有配套的精品资源点击获取 Wasserstein距离与分布鲁棒优化：原理、对偶与电力系统应用
简介这是一份面向机器学习、优化算法研究者与电力系统从业者的分布鲁棒优化实践代码包聚焦以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)其中Π是μ和ν的所有耦合。p1且d取范数时它等于把μ的质量搬到ν所需要的最小运输成本。工程上可以这样理解一个场景的分布偏移主要是“平移”和“形变”造成的而不是概率质量凭空增减。KL散度的定义D_KL(μ‖ν)∫log(dμ/dν)dμ要求μ关于ν绝对连续。用历史数据的经验分布作为名义分布时μ是N个离散点的均匀分布ν是连续分布这时KL会得出无穷大无法作为距离。而Wasserstein距离对离散-连续、离散-离散都给出有限值不需要支撑重叠。φ散度和总变差同样存在这个问题它们更多度量“在同一支撑上概率密度形状差多少”而Wasserstein度量“质量移动了多少”。对电力系统这种场景偏移频繁的问题风电爬坡时刻前后移、负荷峰谷平移Wasserstein的语义更贴切。还有一个计算上的优势当名义分布是经验分布时Wasserstein球的对偶有显式结构这是KL和φ散度没有的。所以项目代码里几乎所有模糊集都以Wasserstein距离为度量这不是偶然。度量支撑不重叠时对场景偏移的惩罚经验分布球的对偶友好度Wasserstein距离有限按距离线性W1强对偶可转有限维LPKL散度无穷大指数级、易爆表名义分布需有密度总变差/φ散度有限但度量整体差异饱和快忽略小偏移约束结构不直观另一个容易被忽略的点是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) Σ_{i1}^N sup_{ξ∈Ξ} ( f(x,ξ) - λ d(ξ, ξ_i) ) ]其中λ是球约束的解耦乘子。这个形式把“Q到μ̂_N距离不超过ε”转成了两层结构外层是λε罚项内层是对每个历史样本ξ_i计算它在支撑集上的最坏偏差。λ是优化变量而不是超参数求解过程中会被自动拉到合适位置。ε固定时λ越大相当于对分布偏移的容忍越小模型越保守。这个转换的价值在于原来的变量是无穷维分布现在只需要在有限个历史样本上加辅助变量。项目代码里最能出问题的一步就在这里——对偶引理成立需要紧支撑而很多实现默认ξ落在全空间上导致内层sup在目标函数线性增长时无上界求解器直接报unbounded。3.2 三种典型的f(x,ξ)结构及对应的约束形式f(x,ξ)结构内层sup的显式结果最终问题类型线性f a_x^T ξ b_xsup落在支撑集顶点配合对偶范数线性规划分段线性f max_k (a_k^T ξ b_k)每段取凸包引入辅助变量线性/MILP凸二次f ξ^T A_x ξ 线性项用二次共轭与2-范数闭式二阶锥规划线性形式最常用。电力系统调度里购电成本、弃风惩罚、失负荷惩罚大多随预测误差线性增长所以f取关于ξ的线性函数。内层的sup_{ξ} { a_x^T ξ - λ‖ξ - ξ_i‖ }在支撑集是凸多面体时可转成顶点上的线性函数加范数锥约束整体保持LP结构。分段线性形式出现在阶梯报价、分段补偿、以及ReLU形状的切负荷惩罚里。对偶后每一段都需要一个辅助变量约束数量约为N×(K1)K是分段数。如果K很大存储和求解成本会明显上升常见做法是把分段函数改写成线性目标加额外整数变量或者用外逼近法逐次加约束。二次形式通常配合2阶Wasserstein距离在均值-方差类风险度量里很常见。它的闭式解很漂亮但求解器从LP变成锥规划大规模场景下比线性形式吃力得多。我一般先按线性把流程跑通验证数据流和参数标定没问题后再根据业务需求升级到分段线性或二次结构。3.3 把一个线性DRO写成求解器模板下面的代码是线性DRO对偶形式的通用求解模板用支持集上的网格近似内层supimport 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(nonnegTrue) # 每个历史样本的最坏场景偏差 S cp.Variable(N, nonnegTrue) # 模糊集半径先用 Parameter 方便后续网格搜索 eps cp.Parameter(nonnegTrue, value0.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×K2000小型问题可以直接求解到N5000、K50时就到25万条约束必须换算法我在第4章末尾展开讲。4. 电力系统场景落地把风电出力不确定性嵌进调度决策4.1 从点预测到“分布鲁棒安全域”风电出力预测误差不是固定的单峰分布。气温、气压、风速模式变化会让预测误差曲线整体平移或尾部变形SAA在这里经常给出“过于自信”的机组启停方案。DRO的做法预测误差的历史样本构成经验分布以它为中心画Wasserstein球调度决策在球内最坏分布下仍然要可行。机组组合里的失负荷概率约束比如99%可靠性被改写成对一个期望惩罚的约束而这个期望是对球内的最坏分布取的。这一步不是简单地把预测区间放宽。传统随机优化的置信区间只覆盖“参数取值的范围”DRO覆盖的是“分布整体的变化范围”。所以同样的调度成本预算下DRO给的备用容量更平滑不会像逐点鲁棒那样同时把每组出力都钉死在最极端场景上。4.2 两阶段DRO机组组合变量与约束结构阶段决策内容变量类型约束对象第一阶段机组启停、备用容量二进制/整数启停逻辑、最小启停时间第二阶段出力调整、弃风量、切负荷量连续功率平衡、爬坡、线路对偶层λ、S_i连续非负每个历史样本的Wasserstein切割两阶段结构里第一阶段变量在观测到风电实际出力之前就要决定第二阶段变量在每一组风电场景下分别调整。用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(nonnegTrue) # 第二阶段每个场景的切负荷 l cp.Variable(N, nonnegTrue) # 对偶层变量 lam_dual cp.Variable(nonnegTrue) S cp.Variable(N, nonnegTrue) eps 1.5 # 模糊集半径MW constraints [] for i in range(N): # 场景功率平衡机组出力 风电实际出力 负荷 - 切负荷 constraints.append(p (50.0 delta_hist[i]) load - l[i]) # Wasserstein 切割对所有支撑点 kS[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个样本周围最坏场景下可能产生的切负荷损失”。eps1.5的单位是MW对应第2章里讲过的量纲问题。多时段扩展时p变成(G,T)矩阵切割约束按时段分别生成目标里累加S矩阵的和。4.3 求解规模的现实反馈对偶后的约束数量大致是N×K乘以时段数。一个T24、N2000、K30的两阶段机组组合切割约束超过百万条直接用通用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, sizeN_train) * 100 # 测试分布中心偏移 右尾加重 test_wind np.random.beta(2, 4, sizeN_test) * 100 5 def simulate_out_of_sample(decision, wind_sample, load100.0, penalty500.0): shortages np.maximum(0.0, load - decision - wind_sample) return penalty * shortages.mean()分别用SAA和DRO求两个决策SAA决策用训练样本直接优化DRO决策用第4章的模型在ε1.5下求解然后都放到test_wind上评估。典型结果趋势如下ε训练损失测试损失最坏损失(测试)0SAA8.241.789.30.59.124.655.21.511.319.841.04.018.921.438.6解读这个表ε0的SAA在训练分布上最便宜测试分布一偏移就崩ε1.5时训练成本略升但测试损失大幅下降ε4.0时训练成本明显偏高测试损失却只降了一点点。这说明存在一个合理ε区间过低不够稳过高过度保守。那张风险-成本曲线就是上线与否的决策依据。从那以后我每次调ε都强制先画一遍风险-成本曲线再谈能不能上线。这个习惯帮我避开了至少三次“论文指标漂亮、现场一跑就崩”的尴尬。希望帮到你。本文还有配套的精品资源点击获取