☰
Wasserstein距离型分布鲁棒优化在电力调度中的MATLAB复现
2026/10/26 7:54:54 网站建设 项目流程

简介:本资源是爱思唯尔论文《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 自带 linprog,1-范数模型也能跑,但大规模算例会慢。

求解时间主要受 N 和 T 影响。N=100、T=24、G=6 的规模,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 线性增长,N=1000 时约束上万条。解决:先做场景削减,或者用随机对偶上升法做分解求解。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,足够工程使用。

最后说个血泪经验:复现这类论文,最耗时间的不是写代码,而是对齐符号和约束方向。我现在的习惯是先把对偶推导写在纸上,每一步标注范数和约束方向,再照着写代码,能省掉大量调试时间。另外,样本外测试一定要做,样本内指标再漂亮,样本外崩了就是白搭。希望帮到你。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询