Wasserstein距离型分布鲁棒优化在电力调度中的MATLAB复现

简介本资源是爱思唯尔论文《Energy and reserve dispatch with distributionally robust joint chance constraints》的MATLAB复现代码面向电力系统调度、分布鲁棒优化方向的研究生与科研人员帮助理解基于Wasserstein距离的联合机会约束建模与求解流程。压缩包共35个文件约5.04MB以m脚本、mat数据文件为主辅以zbak备份、pdf论文与md说明文档涵盖主程序、结果数据与参考材料便于对照论文逐步调试。资源基于MATLAB与yalmip建模调用gurobi求解器完成优化计算可帮助读者掌握目标函数与约束的构建方式、概率分布集合的处理思路以及鲁棒调度策略的验证方法。目前已有124人浏览学习适合作为分布鲁棒优化入门到进阶的实践参考但仅供学习交流请勿用于商业用途。1. 从一篇《Energy and reserve dispatch》说起Wasserstein 距离到底在解决什么电力系统能量与备用联合调度里风电、光伏的预测误差是绕不开的变量。传统随机规划要先假设一个概率分布再用场景法采样可真实误差分布既不是正态也不是某个已知族样本量小的时候场景法给出的调度方案在样本外经常翻车。分布鲁棒优化DRO的思路是不赌单一分布而是在一个以经验分布为中心的模糊集里找最坏分布下的最优决策。模糊集怎么定义直接决定模型是保守到没法用还是激进到没意义。Wasserstein 距离型模糊集是这几年被反复验证过的一种选择它用经验分布和真实分布之间的搬运代价来度量偏离程度半径 ε 控制鲁棒程度ε 取 0 退化成样本平均近似ε 取大就逼近全局最坏情况。好处是它同时吃到了样本信息和结构信息且对偶之后能把无限维的分布优化转成有限维的线性规划或二阶锥规划MATLAB 里用 YALMIP 或 CVX 就能落地。这篇要复现的《Energy and reserve dispatch》正是这条路线上的典型工作适合做电力调度、综合能源、机组组合方向、手里有 MATLAB 又想把 DRO 真正跑起来的人。下面按模型怎么搭 → 对偶怎么推 → 代码怎么写 → 坑在哪的顺序讲透。2. 把 Wasserstein 模糊集写进调度模型从原始问题到可解形式2.1 两阶段结构与不确定量的建模边界能量与备用联合调度的标准结构是两阶段第一阶段决定机组基荷出力、备用容量和机组启停第二阶段在不确定量实现后决定再调度量。不确定量通常取风电出力的预测误差记作 ξ维度等于风电场数量。经验分布由 N 个历史样本构成每个样本权重 1/N。Wasserstein 距离的定义是给定两个分布 μ、ν一阶 Wasserstein 距离等于所有联合分布中运输代价 E[||ξ - ξ||] 的最小值。基于它构造的模糊集是以经验分布为中心、半径 ε 的球。这里有个容易忽略的点范数选择会改变对偶形式。用 1-范数或 2-范数对偶后分别得到线性约束和二阶锥约束MATLAB 里求解器支持程度不同我一般先用 1-范数跑通再换 2-范数对比保守度。半径 ε 不是拍脑袋定的。常见做法是用样本量 N 和置信水平 β 反推ε 大致正比于 sqrt(1/N)·log(1/(1-β))具体常数取决于支撑集直径。样本少的时候 ε 必须放大否则模糊集覆盖不到真实分布样本外表现会崩。2.2 对偶推导把 min-max 问题变成线性规划原始问题是 min_x max_P E_P[cost(x, ξ)]P 落在 Wasserstein 球内。直接求解不可行因为 P 是无限维的。标准做法是引入对偶变量 λ对应 Wasserstein 约束和一组辅助变量把内层 max 转成 min。推导结果大致是外层决策 x 不变内层多出一组约束形式为 cost(x, ξ_i) γ, ξ_i - ξ_j ≤ λ·ε ...其中 γ 是对偶变量受范数约束。1-范数下 γ 的每个分量落在 [-1, 1] 区间2-范数下 ||γ||₂ ≤ 1。这一步是整个复现的核心推错一个符号结果就会偏保守或偏激进。提示对偶推导建议手推一遍再写代码。我见过不少人直接抄公式结果范数约束写反跑出来的备用容量比确定性模型还小明显不合理。2.3 在 MATLAB 里搭出可求解的模型骨架下面这段是模型骨架用 YALMIP 描述求解器选 Gurobi 或 CPLEX。变量命名和论文里的符号尽量对齐方便对照。% 参数T 时段数G 机组数W 风电场数N 样本数 % Pmax/Pmin 机组出力上下限Rmax 备用上限C 成本系数 % xi 为 W×N 的风电误差样本矩阵eps_w 为 Wasserstein 半径 yalmip(clear); P sdpvar(G, T); % 基荷出力 R sdpvar(G, T); % 备用容量 gamma sdpvar(W, T); % 对偶变量 lambda sdpvar(1, T); % 对偶变量 Constraints [Pmin P Pmax, 0 R Rmax]; for t 1:T % 功率平衡确定性部分 Constraints [Constraints, sum(P(:,t)) load(t)]; % 范数约束1-范数下 gamma 分量有界 Constraints [Constraints, -1 gamma(:,t) 1]; for n 1:N % 对偶后的内层约束cost 用线性成本近似 Constraints [Constraints, ... C * (P(:,t) R(:,t)) gamma(:,t) * xi(:,n) lambda(t)]; end Constraints [Constraints, lambda(t) 0]; end Objective sum(sum(C .* P)) sum(sum(C .* R)) sum(lambda) * eps_w; ops sdpsettings(solver, gurobi, verbose, 1); diagnostics optimize(Constraints, Objective, ops);逻辑说明外层目标包含基荷成本、备用成本和鲁棒项 λ·ε。内层约束对每个样本 n 都要成立这是对偶后把分布优化转成有限约束的直接结果。gamma 的范数约束对应 Wasserstein 距离里的运输代价范数改成 2-范数时换成 norm(gamma(:,t)) 1。参数说明eps_w 是模糊集半径从 0 开始逐步增大观察备用容量和总成本的变化曲线通常会在某个值之后趋于平缓那个拐点就是比较合理的鲁棒水平。N 是样本数样本越多 eps_w 可以越小但约束数量线性增长N 超过 500 后求解时间会明显上升建议先做场景削减。3. 复现时的数据、求解器与结果校验让曲线对得上论文3.1 测试系统与数据准备论文一般用 6 机系统或 IEEE 30 节点系统做算例。复现时如果拿不到原始数据常见做法是用 MATPOWER 里的 case6 或 case30 改造成本和备用参数风电误差样本用历史预测误差或按经验分布生成。样本生成要注意不要用正态分布生成完就当真值那样 Wasserstein 半径的意义会被削弱最好用带偏的分布或真实误差数据。% 用 MATPOWER 读取 6 机系统并提取机组参数 mpc loadcase(case6ww); G size(mpc.gen, 1); Pmax mpc.gen(:, 9); % 最大出力 Pmin mpc.gen(:, 10); % 最小出力 C mpc.gencost(:, 5); % 线性成本系数 % 风电误差样本用 t 分布生成带厚尾 rng(1); xi trnd(3, W, N) * 0.1 * mean(Pmax);逻辑说明case6ww 自带风电接口适合做小规模验证。trnd 生成厚尾样本比正态更接近真实预测误差。参数 0.1*mean(Pmax) 控制误差幅度可按实际风电渗透率调整。3.2 求解器选择与求解时间控制YALMIP 支持 Gurobi、CPLEX、MOSEK 等。1-范数模型是纯线性规划Gurobi 和 CPLEX 都能秒解2-范数模型是二阶锥规划Gurobi 9 以后支持MOSEK 更稳。如果只有 MATLAB 自带 linprog1-范数模型也能跑但大规模算例会慢。求解时间主要受 N 和 T 影响。N100、T24、G6 的规模Gurobi 通常在 10 秒内。N 到 500 以上建议先做场景削减用 k-means 或快速前向选择把样本压到 100 以内再跑 DRO。削减会引入近似误差但比直接跑不动要好。3.3 结果校验三条必须对的曲线复现是否成功看三条曲线。第一总成本随 ε 增大而单调上升上升速度先快后慢。第二备用容量随 ε 增大而增加但不会超过机组备用上限之和。第三样本外测试用另一组未参与训练的误差样本计算实际再调度成本DRO 应该优于样本平均近似ε0和确定性模型。% 样本外测试固定第一阶段决策用新样本算再调度成本 xi_test trnd(3, W, 200) * 0.1 * mean(Pmax); cost_out zeros(200, 1); for k 1:200 % 给定 P、R求解第二阶段再调度 [cost_out(k), ~] solve_recourse(value(P), value(R), xi_test(:,k), C, Pmax, Pmin); end fprintf(DRO 样本外平均成本: %.2f\n, mean(cost_out));逻辑说明solve_recourse 是第二阶段子问题给定第一阶段决策和误差实现求最小再调度成本。对比 ε0 和确定性模型的样本外成本如果 DRO 没有优势说明 ε 选得不对或对偶推导有问题。注意样本外测试的样本必须和训练样本独立。我见过有人用同一组样本既训练又测试结果 DRO 看起来完美实际部署时完全不是那么回事。4. 避坑与排查复现 Wasserstein DRO 调度时最容易翻车的五处现象一模型无界或求解器报 infeasible。原因通常是对偶约束里 lambda 没有下界或者 gamma 的范数约束写成了等式。解决给 lambda 加 0 约束gamma 用不等式约束检查每个时段的功率平衡是否和负荷数据对齐。现象二ε 增大成本反而下降。这是对偶符号写反的典型症状。原因在于内层 max 转 min 时约束方向搞错导致模糊集越大模型越松。解决手推一遍对偶确认约束是 cost gamma, xi lambda 而不是 。现象三样本外表现和样本内差不多DRO 没优势。原因可能是 ε 太小模糊集退化成经验分布或者样本本身没有厚尾Wasserstein 距离区分度不够。解决逐步增大 ε观察样本外成本曲线找到最低点同时检查样本生成方式确保误差分布不是退化的。现象四N 增大后求解时间爆炸。原因是对偶约束数量随 N 线性增长N1000 时约束上万条。解决先做场景削减或者用随机对偶上升法做分解求解。MATLAB 里可以用 parfor 并行构造约束但 YALMIP 的约束拼接本身是瓶颈建议改用 Gurobi 的矩阵接口直接建模。现象五换 2-范数后结果和 1-范数差很多。这是正常的2-范数模糊集更紧保守度更低。但如果差异大到备用容量差一倍说明范数约束的缩放没做好。解决检查 gamma 的范数约束是否和运输代价的范数一致2-范数下应该是 ||gamma||₂ 1不是 ||gamma||₂ sqrt(W)。5. 进阶技巧用对偶变量反推模糊集半径与场景削减的联动跑通基础模型之后真正影响落地效果的是 ε 的选取和场景削减的配合。我一般用对偶变量 lambda 的取值来反推如果最优解处 lambda 全部贴着下界说明 ε 太小模糊集没起作用如果 lambda 很大且多个时段接近说明 ε 偏大模型过于保守。一个实用的做法是二分搜索 ε以样本外成本为目标通常 5 到 8 次迭代就能找到拐点。场景削减和 ε 是联动的。削减后样本数 N 变小理论上 ε 要按 sqrt(1/N) 放大否则模糊集覆盖不足。但放大 ε 又会增加保守度所以削减比例不宜超过 50%。我习惯先保留 100 个样本跑基准再削减到 50 个并相应放大 ε对比样本外成本如果差异在 2% 以内就认为削减可接受。% 二分搜索最优 epsilon eps_lo 0; eps_hi 1.0; for iter 1:8 eps_mid (eps_lo eps_hi) / 2; cost_out evaluate_out_of_sample(eps_mid, xi_train, xi_test); if cost_out best_cost best_cost cost_out; eps_lo eps_mid; else eps_hi eps_mid; end end fprintf(推荐 epsilon: %.4f, 样本外成本: %.2f\n, eps_lo, best_cost);逻辑说明evaluate_out_of_sample 内部先求解 DRO 得到第一阶段决策再用测试样本算再调度成本。二分搜索假设样本外成本关于 ε 是单峰的实践中大多数算例满足这个假设如果不满足改用网格搜索。参数说明eps_hi 的上界取 1.0 是经验值对应误差幅度 10% 左右的算例。如果风电渗透率更高上界要相应放大。迭代次数 8 次对应精度约 0.004足够工程使用。最后说个血泪经验复现这类论文最耗时间的不是写代码而是对齐符号和约束方向。我现在的习惯是先把对偶推导写在纸上每一步标注范数和约束方向再照着写代码能省掉大量调试时间。另外样本外测试一定要做样本内指标再漂亮样本外崩了就是白搭。希望帮到你。本文还有配套的精品资源点击获取