概率距离快速削减法实现风光场景生成与削减：MATLAB代码详解

做电力系统不确定性分析的朋友大概都绕不过“场景”这两个字。风电、光伏出力一会儿高一会儿低光伏到了晚上直接归零你要是拿一条确定曲线去做调度、做规划结果基本没法用。所以大家习惯先生成一大堆可能的风光出力场景再用场景削减方法挑出少量有代表性的场景——这一步做好后面计算量能降一个量级精度还不会损失太多。我最近用MATLAB写了一套基于概率距离快速削减法的风光场景生成与削减程序实测效果不错思路也不算复杂。这篇就把算法原理、完整代码、参数调优和踩过的坑一次说清楚适合正在做随机规划、鲁棒优化或者风光并网仿真的同学直接拿去参考。1. 场景生成与削减到底在解决什么问题1.1 为什么需要一堆“可能发生”的场景风电和光伏出力属于典型的强随机过程风速受地形和天气影响光照受云层和昼夜影响单一的“预测曲线”只是期望值远远刻画不了真实波动。在电力系统的日前调度、机组组合、储能容量配置这类决策里我们往往需要把不确定性显式建模进去于是就有了“场景法”。场景法的大致思路是先用蒙特卡洛采样、拉丁超立方采样或者历史数据扰动生成大量可能的出力时间序列每个序列带一个概率权重构成一个离散概率分布。理论上场景数越多越接近真实分布但计算量也会跟着爆炸。比如一个含1000个场景、每个场景24个时段的两阶段随机优化模型求解时间可能是几十倍甚至上百倍的增长。更麻烦的是许多商业求解器在场景数超过一定规模后内存和耗时都让人无法接受。这时候就需要场景削减。场景削减不是简单随机抽几条曲线而是要在“减少场景数量”和“保持原有概率分布特性”之间找平衡。专业一点说就是找一个规模更小的离散概率分布让它与原始场景集合的某种概率距离最小。我们这里用的概率距离快速削减法就是这一类方法里的主流方案。1.2 概率距离快速削减法为什么值得用削减场景有多种路径最常见的是三类K均值聚类、基于启发式筛选的历史场景抽样以及基于概率距离的场景削减。K均值聚类直观但聚类出来的“中心曲线”是平均后的结果可能抹掉极端情况而且聚类结果依赖初始质心跑一遍和跑两遍可能不一样。随机抽样更简单但样本量小时代表性很差容易出现概率失真。概率距离快速削减法在数学上更有底气。它把场景削减看作两个概率分布之间的近似问题用Kantorovich距离也叫Wasserstein-1距离衡量削减前后分布差异通过迭代合并或删除场景使概率距离增量最小。这样做有几个直接好处保留场景来自原始场景集而不是合成出的平均曲线因此场景的物理形态不会失真被删除场景的概率会累加到距离它最近的保留场景上保留了概率质量每一步只做局部最优决策计算复杂度可控适合工程实现在随机优化里削减后的场景集能比较好地保持最优解的稳定性。“快速削减”中的“快速”体现在实现上只需要预先算好场景间距离矩阵每轮迭代更新有限行和列不需要重新做全局优化也不依赖外部工具箱。这也是我选它写MATLAB实现的原因纯代码就能跑不给读者增加额外负担。2. 概率距离快速削减法的数学原理2.1 场景概率与Kantorovich距离先定义最基本的符号。假设原始场景集合有 N 个场景第 i 个场景是长度为 T 的出力序列 ξ_i对应的概率为 p_i且 Σ p_i 1。削减后保留了 J 个场景J N对应的场景和概率分别记为 ψ_j 和 q_j。要衡量两个离散分布之间的差异常用Kantorovich距离公式写出来是这样的D(P, Q) inf{ Σ_{i,j} π_{ij} c(ξ_i, ψ_j) }其中 π 是两个分布之间的联合分布耦合c(ξ_i, ψ_j) 表示地面距离一般取欧氏距离 ||ξ_i - ψ_j||。可以把这个公式理解为“把一个分布的概率质量搬到另一个分布上平均搬运成本最小是多少”。如果削减后的分布与原始分布很接近那么搬运成本就小距离就小。生活里的类比你有两筐不同重量的苹果要从一筐挪到另一筐每个苹果搬一米的成本相同最终搬完的总重量乘以距离就是搬运成本。Kantorovich距离就是找一种搬法让总成本最低。场景削减就是想找到一小筐苹果使得搬运成本最低。对离散场景来说这个最优耦合可以简化。削减后每个保留场景 ψ_j 会“吸收”一部分原始场景概率最理想的情况是把每个被删除场景 ξ_i 的概率加到距离它最近的保留场景上这时总搬运成本就是Σ_{i 被删除} p_i · min_j ||ξ_i - ψ_j||这也是概率距离削减法的目标函数。我们不需要真的求最优耦合只需要按这个原则去删场景。2.2 后向削减的标准步骤与“快速”实现我实现的是后向削减Backward Reduction它从完整场景集合出发每一步删掉一个场景直到剩余场景数为 J。每一步选择删除哪个场景的标准是让上述总搬运成本的增量最小。具体步骤如下计算所有场景两两之间的欧氏距离矩阵 DD(i,j) ||ξ_i - ξ_j||对角线设成无穷大。标记所有场景为“活动”状态。循环直到活动场景数为 J对每个活动场景 i计算它到其他活动场景的最近距离 d_i min_{j ≠ i} D(i,j)计算概率加权距离 w_i p_i · d_i选择 w_i 最小的场景 i* 删除并把它的概率累加到距离它最近的场景 j* 上把场景 i* 从活动集合中移除在距离矩阵中把对应行列置为无穷大。循环结束后活动场景就是保留场景其概率就是累加后的概率。为什么删除 w_i 最小的场景因为 w_i 表示“如果删除场景 i需要把它搬运到最近保留场景的概率加权距离”。选择最小的删除对整体分布扰动最小。这就是“贪心”思想每一步都找当前最优解。虽然不保证全局最优但在实际工程中效果非常好而且计算复杂度可控。“快速”体现在哪里呢主要有两点。第一距离矩阵只在初始时算一次后续只更新被删除场景的行列不需要每次重新计算两两距离。第二每轮删除时只需要对活动场景做一次向量取最小操作避免了嵌套双重循环里的重复搜索。MATLAB 里用 min 和逻辑索引完全可以向量化N1000、J10 时运行时间基本在秒级。2.3 削减效果的评价指标削减完之后不能光说“看着差不多”要用指标验证。我常用三个指标概率保留比例削减后所有保留场景概率之和应该仍为1检查累加后是否接近1避免浮点误差累积。Kantorovich距离计算原始分布与削减后分布之间的距离距离越小说明分布保真度越高。如果没有现成函数可以用删除场景的加权距离来衡量这个值在削减过程中是逐步增大的J 越小距离越大。期望曲线对比原始场景集和缩减场景集分别求每个时刻的期望值也就是概率加权平均两条期望曲线应该足够贴近。这个指标很直观也最容易跟非专业人士解释。另外还可以看方差、极端场景是否存在。比如削减后是否保留了出力为0的场景、是否保留了出力极低的“枯风期”场景这些对决策影响很大不容忽视。3. MATLAB核心实现与代码拆解3.1 初始场景生成以风速为例场景削减需要先有一批初始场景。这里提供一个简单但好用的风速场景生成器以一条基准出力曲线为期望在每个时段叠加一个时间相关的随机扰动。扰动用 AR(1) 模型生成这样相邻时段之间的波动更自然不会出现每时刻独立随机导致的高频毛刺。function scen generate_wind_scenarios(base, nScen, sigma) % base: 1 x nTime 基准出力曲线标幺值范围0~1 % nScen: 生成的场景数量 % sigma: 扰动的标准差比例 nTime length(base); scen zeros(nScen, nTime); % AR(1) 系数控制时间相关性 phi 0.6; for s 1:nScen e zeros(1, nTime); e(1) randn * sigma; for t 2:nTime e(t) phi * e(t-1) randn * sigma * sqrt(1 - phi^2); end scen(s, :) base e; % 限幅到[0,1] scen(s, :) max(min(scen(s, :), 1), 0); end end这里 AR(1) 的扰动序列稳态方差等于 sigma^2因为用了 sqrt(1 - phi^2) 做innovative噪声缩放。你可能会问为什么不用单纯 randn因为风电相邻时段出力不可能完全独立夜里风速也不会瞬间从0跳到满发时间相关性能让场景更贴近真实物理过程。光伏场景的生成逻辑类似只不过基准曲线要带“白天高、夜里为0”的凸峰形状扰动只作用在日照时段夜里直接置零。这里不展开细写方法是一样的。3.2 场景削减函数fast_backward_reduction核心削减函数是整篇博文的重点我把它单独拆出来讲。输入是场景矩阵、概率向量、目标保留场景数输出是削减后的场景矩阵和概率向量。代码使用向量化实现避免了三重循环性能要好很多。function [scen_red, prob_red] fast_backward_reduction(scen, prob, J) % 概率距离快速后向场景削减 % 输入: % scen - N x T 矩阵每行一个场景T个时段 % prob - N x 1 向量初始场景概率要求和为1 % J - 目标保留场景数 % 输出: % scen_red - J x T 保留场景 % prob_red - J x 1 削减后概率 % % 思路: 每轮删除一个场景该场景被删除后 % 其概率加权最小距离最小即对Kantorovich距离增量最小。 N size(scen, 1); if N J scen_red scen; prob_red prob; return; end % 计算场景间欧氏距离矩阵 D pdist2(scen, scen); D(1:N1:end) inf; % 对角线置为inf防止选到自己 deleteFlag false(N, 1); % 为了提高效率在循环里只对活动场景计算 while sum(~deleteFlag) J active find(~deleteFlag); nActive length(active); % 提取活动场景对应的距离子矩阵 Dsub D(active, active); % 对角线置inf diagIdx sub2ind([nActive, nActive], 1:nActive, 1:nActive); Dsub(diagIdx) inf; % 每个活动场景到其他活动场景的最近距离 [minDist, minIdxLocal] min(Dsub, [], 2); % 概率加权距离增量 weighted prob(active) .* minDist; % 选择加权距离最小的场景作为删除对象 [~, delLocal] min(weighted); del active(delLocal); keepLocal minIdxLocal(delLocal); keep active(keepLocal); % 概率累加到最近的保留场景上 prob(keep) prob(keep) prob(del); % 标记删除并更新距离矩阵 deleteFlag(del) true; D(del, :) inf; D(:, del) inf; end scen_red scen(~deleteFlag, :); prob_red prob(~deleteFlag, :); % 重新归一化避免浮点误差导致概率和不为1 prob_red prob_red / sum(prob_red); end几个关键点你需要注意。第一概率累加是加到“最近场景”上而不是平均分配。这符合Kantorovich距离的最优搬运逻辑。第二每轮删除后D 矩阵中删除场景的行列被置为 inf下一轮提取 Dsub 时已经删除的场景不会参与距离计算。第三minIdxLocal 返回的是在 Dsub 子矩阵中的列索引要映射回原始场景编号 active(minIdxLocal)。你可以直接复制这个函数到 MATLAB 脚本里运行。pdist2 在 Statistics and Machine Learning Toolbox 里如果不想依赖工具箱可以用sqrt((scen*scen - 2*scen*scen scen*scen))自己算或者用D sqrt(sum(abs(scen).^2, 2))的展开式。但一般大家手里都有统计工具箱直接用 pdist2 最省事。3.3 主程序串联与可视化验证把生成和削减串起来的主程序如下%% 参数设置 rng(42); T 24; % 24个时段 N 500; % 初始场景数 J 10; % 保留场景数 sigma 0.15; % 扰动标准差 %% 生成基准风速曲线简易日变化波形 t 0:T-1; base 0.3 0.45 * sin(2*pi*t/24 - pi/2); base max(min(base, 1), 0); %% 生成500个风速场景 scen generate_wind_scenarios(base, N, sigma); prob ones(N, 1) / N; %% 场景削减 [scen_red, prob_red] fast_backward_reduction(scen, prob, J); %% 可视化对比 figure(1); plot(scen, Color, [0.75 0.75 0.75]); hold on; plot(scen_red, LineWidth, 2); xlabel(时段/h); ylabel(风速标幺值); title(削减前后的风电场出力场景); legend(原始场景, 削减后场景, Location, best); grid on; %% 期望曲线对比 exp_orig sum(scen .* prob, 1); exp_red sum(scen_red .* prob_red, 1); figure(2); plot(exp_orig, k-, LineWidth, 1.5); hold on; plot(exp_red, r--, LineWidth, 1.5); xlabel(时段/h); ylabel(期望出力); legend(原始期望, 削减后期望); title(剔除前后期望曲线对比); grid on;运行之后你会看到 gray 背景的原始场景逐渐收敛到几条加粗的代表性曲线。期望曲线对比通常非常接近这就说明削减后的场景集很好地保留了整体统计特征。我建议你在此基础上把削减后的概率打印出来看一眼disp([scen_red, prob_red]);削减后的概率一般不会像初始时那样均匀而是有的场景权重高、有的权重低。权重高的场景往往代表“典型状态”权重低的场景代表“边缘状态”。这个概率分布信息后续进入随机优化模型时非常关键。4. 参数选择、常见报错与我的实测心得4.1 关键参数怎么调最核心的参数有三个初始场景数 N、保留场景数 J、生成扰动强度 sigma。初始场景数 N 决定了样本对真实分布的覆盖程度。N 太小比如只有100个削减出来的场景容易漏掉极端情况N 太大比如5000个距离矩阵一次要算12.5 million个元素MATLAB里会明显卡顿。我的经验是对24时段的场景500到2000是一个比较合适的区间。如果你的后续优化模型吃得消1500个初始场景削减成15个效果通常好于300个初始场景削减成15个。保留场景数 J 取决于你后续问题的复杂度和对精度的要求。J 太少分布失真严重J 太多随机优化求解负担大。我一般从 J 5 开始试画期望曲线对比如果 Kanto 距离不再明显下降就停。实际操作中很多论文和工程案例会把初始场景削减到1020个日度调度用10个左右足够。扰动强度 sigma 要结合基准曲线的量级来看。风速标幺值的 sigma 取 0.10.2 比较常见。如果把某个时段的基准出力是0.5sigma0.15那么95%的扰动范围大约在 ±0.3 内看起来是合理的。sigma 太小场景都挤在期望线附近削减出来的场景太“乐观”sigma 太大限幅操作会把大量场景顶到0和1产生不真实的“满发/零出力”堆积。需要根据历史数据的实际波动标定不能拍脑袋。4.2 典型报错与排查对照表我在把代码给其他同学用的时候遇到最多的几个报错和问题整理成一张表给你参考。现象可能原因解决方法pdist2 未定义缺少 Statistics and Machine Learning Toolbox改用自定义距离函数或用sqrt(sum((scen - scen) .^ 2, 3))的三维广播方式注意内存削减后概率之和不为1浮点误差累积在函数末尾加prob_red prob_red / sum(prob_red);削减后保留场景数比设定的多循环条件写错或 while 写成大于等于检查sum(~deleteFlag) J减到 J 就停场景间距全为 inf报 min 输出空值距离矩阵中删除行未正确置inf导致活动集合里 Dsub 全inf确认每次删除后设置D(del, :) inf; D(:, del) inf;削减结果全是同一条曲线初始场景扰动太弱所有曲线太相似增大 sigma或者检查生成函数是否限幅把场景都拉平了代码运行极慢初始场景数过大或用了双层循环找最近对用向量化版本避免在循环里嵌套双重 for期望曲线在零附近波动但异常高光伏场景生成了负值或夜间非零值生成后统一把日落时段置零并做限幅这里最容易被忽视的是场景削減完没有重新归一化概率。因为累加操作会引入浮点误差特别是删了上百个场景后概率和可能是0.9999或1.0001带入优化模型就会触发概率约束不闭合的报错。养成好习惯输出前一定要归一化。4.3 几个容易踩的坑及优化建议第一个坑是“删除场景数太猛导致概率失衡”。比如把500个场景削减到5个个别保留场景概率可能高达0.7这时候某个“关键场景”概率偏差很大会影响优化结果。我的处理办法是如果保留场景概率出现极端集中可以选择稍微增加 J或者改用前向选择法。前向选择法不是删除场景而是从零开始往保留集里逐个加场景能更好地避免一开始把所有概率堆到少数场景上。第二个坑是距离函数的选择。欧氏距离简单直观适合风速、光伏出力这种数值型时序数据。但如果你的场景代表“电价曲线”可能更关心形状而非幅值这时候可以用动态时间规划距离或1减去相关系数。换距离只需要改 Dsub 那一步即可其他逻辑不用动。第三个坑与随机种子有关。场景生成用了 randn如果不固定 rng每次运行结果不一样。工程上写文档、交报告时一定要固定随机种子否则结果不可复现。我习惯在脚本开头写rng(2024)这种固定值方便别人验证。第四个优化建议是在大场景量级下用分阶段削减。比如要削减到10个可以先从2000削减到100再从100削减到10。虽然理论上不符合一次削减的最优性但实践中能缓解距离矩阵内存压力而且最终效果差别不大。如果初始N5000这一步几乎是必须的。第五个建议是输出削减过程的“删除映射表”也就是每个删除场景最终被合并到了哪个保留场景。这个信息在做场景到场景的灵敏度分析时非常有用。代码里deleteFlag只记了状态没有记录最终归属。你可以维护一个 N 长的数组map (1:N)每次删除 del 时执行map(del) map(keep)最后所有原始场景都能映射到保留场景编号上这样后续分析能追踪概率来源。最后再分享一个小技巧在削减完成后用prob_red作为权重做后续随机优化而不是把保留场景当等概率场景。很多初学者削减完就把概率忽略了等概率加权会严重误估期望成本。保留场景的概率本身就是削减算法最重要的输出一定要传进优化模型里。这个细节恰恰是概率距离快速削减法区别于普通聚类的核心价值所在。