基于概率距离快速削减法的风光场景生成与削减MATLAB实现

做新能源并网研究的朋友，大概都遇到过这种尴尬：手里攒了一整年的风光出力实测数据，或者用蒙特卡洛抽了几千个场景，结果优化模型一跑就卡死，哪怕勉强跑通，结果也没有参考价值。商用软件里自带的多场景工具要么是黑盒，要么削减效果说不清楚。今天聊的这个MATLAB实现，叫“基于概率距离快速削减法的风光场景生成与削减方法”，一句话说就是把上千个随机出力场景，按概率分布损失最小的原则压缩成几十个甚至几个代表性场景，同时保留每个场景的概率权重。它很适合用在随机优化、风险评估、可靠性分析里。不管你是刚开始接触场景法的新人，还是已经在调场景规模的老手，这篇的内容多少都能给你一些参考。

1. 风光场景生成与削减到底在解决什么问题

1.1 为什么场景数量会失控

风电和光伏出力受天气影响，本质上是一个随时间变化的高维随机过程。做随机优化时，常规思路是用大量离散场景逼近真实分布，比如对风速、光照预测误差采样，或者直接历史整年功率序列抽样。单看每一个场景，就是一条确定性的时序曲线；把几百上千条放在一起，就构成一个场景集。

麻烦也在这里：场景集是蒙特卡洛风格，数量一多，优化问题里的变量数和约束数直接乘以场景数，混合整数规划根本解不动。我之前在一个两阶段随机调度模型里试过，1000个场景规模，内存占用让人崩溃，求解器跑了两小时有余。所以必须做削减，用少量场景尽可能还原原始场景集的概率分布特征。这个工作相当于“把实测数据的体积压缩但保留信息密度”，概率距离快速削减法干的就是这件事。

1.2 场景生成到削减的完整链路

在“场景法”建模流程中，第一件事是生成初始场景集，常见手段包括：基于历史数据重采样、用ARMA或者马尔可夫链拟合随机过程、基于Copula描述多变量相关性、或者直接搭一个预测误差模型。第二件事才是削减。削弱的场景集最终作为随机参数输入到优化模型中。

概率距离快速削减法属于后处理环节，它不关心场景是怎么来的，只关心输入场景集S、每个场景的概率p、削减目标数量K，输出是保留哪些场景以及保留概率。这个位置很关键，因为无论上游生成方法是什么，只要你还要用场景做计算，就几乎绕不开这一步。它削得好不好，直接决定了随机优化结果的可信度。

2. 概率距离快速削减法的原理细节

2.1 场景之间的距离怎么定义

场景削减的前提是“距离”。对两个N维向量的场景来说，最常见的距离是欧氏距离。如果把一条时序曲线看作一个点，两条曲线之间的“空间距离”可以直接用向量二范数算：

d(i,j)=∥S_i−S_j∥₂

在MATLAB里就是vecnorm(S_i - S_j, 2)，或者pdist2一把全算出来。欧氏距离实现简单，但有时不合适。比如风力发电场景时序存在“时间偏移”，同样一股风早一小时和晚一小时通过，欧氏距离可能很大，而动态时间规整DTW能给出更合理的相似度。不过DTW计算量大，在几百上千个场景配对时开销不小，所以不是默认选项。

如果场景集来自不同量纲的数据混合，比如风速、光照、负荷放在一起，必须先做标准化，否则量级大的变量会把距离计算结果主导，导致削减结果几乎只看那个变量。做标准化不是必须的，但绝大多数情况下都推荐。单一功率场景时，因为量纲一致才可以直接计算。

2.2 概率距离与Wasserstein距离

这里说的“概率距离”，笼统讲是衡量两个概率分布之间差异的度量。离散场景集分布可以写成每个场景点上的概率堆。用最优运输理论的语言，削减后的场景集与实际分布之间的Kantorovich距离（也叫Wasserstein距离）有一个性质：在离散情况下，它等于按最优方式把概率从原始分布“搬”到削减后分布的最小成本，而搬移成本是场景间的距离函数。快速削减法就是在迭代过程中，最小化这个搬移损失。

道理其实不难理解。生活里你请十几个朋友吃饭，不能每个菜都点，只能点几个代表菜，但必须保证大家能吃饱且口味类型大致覆盖。每代表一个“场景”被点掉，你把原先分给这个菜的份额转移给另一个最像它的菜，转移成本就是口味差异。场景削减也是一样：删掉某些场景后，把它们的概率叠加到“最相似”的保留场景上，让概率质量不丢失，只是转移了一小段距离。

2.3 后向削减与前向选择的取舍

快速削减法主要有两种实现路线：后向削减（backward reduction）和前向选择（forward selection）。很多教材里默认的“快速削减法”指前向选择，但MATLAB代码里两者都有人写。

后向削减的流程是：

1. 初始化所有场景等概率p_i=1/N。
2. 计算场景间距离矩阵D。
3. 当剩余场景数大于K时，循环：
   • 找到距离最近的一对场景(i,j)；
   • 删除其中一个（比如i），把p_i累加到p_j；
   • 更新距离矩阵，只保留未被删除的场景。
4. 剩余场景就是削减结果。

前向选择则是从空集开始，每次从未选中的场景中挑一个加入保留集，加入标准是“新增这个场景后，保留集到原始场景集的概率距离增量最小”。前向选择计算代价更大，但削减精度通常更稳定，尤其当K很小的时候。

用通俗话讲，后向削减是不断“淘汰”人，前向选择是不断“面试”人。淘汰简单直接，面试更挑剔但选出来的人可能更均衡。实现时先用后向削到给定K，如果发现结果明显偏，再换前向选择对比一下。

2.4 为什么这种方法比聚类削减更常用

可能有人会问：场景削减直接k-means聚类不就行了吗？k-means也能把N个场景聚成K类，取聚类中心。问题是，聚类中心往往是原始数据坐标空间里的“内插点”，它可能不是任何一条真实能看到的时间序列，随机优化里你没法直接把它当作一个真实出力场景，还需要再做修正。另外k-means对初始中心敏感，高维时序场景容易陷入局部最优。

概率距离快速削减法保存的是原始场景集合中的“代表性场景”，这些场景都是真实采样序列，不需要额外重构。概率转移计算也有明确的理论解释，在论文里解释起来更能站得住脚。这就是它成为随机优化场景削减标准方法的原因。

3. MATLAB代码实现与核心模块拆解

3.1 主函数数据结构设计

写MATLAB代码，第一步是定义清楚输入输出。我习惯把主函数设计成：

function [redIdx, redProb, loss] = scenarioReduction(S, p, K, distType) % S: N-by-T 矩阵，每行一个场景，共N个场景，每个场景T个时间点 % p: N-by-1 概率向量，缺省时等概率 ones(N,1)/N % K: 目标保留场景数 % distType: 'euclidean' 或 'weighted'，缺省为 'euclidean' % 返回： % redIdx: K-by-1 保留场景在原矩阵中的行索引 % redProb: K-by-1 对应的概率，和为1 % loss: 削减过程造成的概率距离损失

这里有一个关键设计选择：场景矩阵S的行索引用原始索引，而不是重新编号。好处是后续如果要把削减后的场景和原始数据打通，不需要额外记录映射关系。熟手通常会在函数内把工作副本和原始索引分开维护，避免删行时索引错乱。

3.2 距离矩阵的高效计算

场景削减需要反复使用场景间距离。最直接的方式是两层循环：

D = zeros(N, N); for i = 1:N for j = i+1:N D(i,j) = norm(S(i,:) - S(j,:)); D(j,i) = D(i,j); end end

这样写逻辑清晰，但N=1000时循环约50万次，MATLAB跑起来有些慢。我建议用pdist2一步到位：

D = pdist2(S, S, 'euclidean');

对于N=1000、T=24，内存占用约800万浮点数，64MB左右，完全能接受。但如果N到5000甚至10000，D矩阵已经到400MB以上，容易卡死。这时就不要全量一次性生成，改成按行分块计算，或者直接用pdist算出紧凑矩阵再用squareform，但pdist同样会先在内存里生成完整矩阵，本质没区别。折中办法是分块：例如每次只算D(i,:)，并只保存对角以下的上三角部分？实际上削减过程中每次删掉一个场景后需要更新距离矩阵，分块更新比较麻烦。实践上N超过3000时，我会先随机抽样或者用启发式方法把规模降下来，再用后向削减；或者改用“快速前向选择”来降低计算量。

3.3 后向削减核心循环的写法

给出一个可以跑的骨架代码：

function [redIdx, redProb, loss] = scenarioReduction(S, p, K, distType) if nargin < 3, error('至少需要 S, p, K 三个输入'); end if nargin < 4, distType = 'euclidean'; end [N, ~] = size(S); if nargin < 2 || isempty(p) p = ones(N, 1) / N; end % 初始场景索引 idx = (1:N)'; D = pdist2(S, S, distType); loss = 0; while length(idx) > K % 找距离最小的一对场景 dmat = D(idx, idx); % 去除对角线 dmat(logical(eye(length(idx)))) = inf; [minVal, linIdx] = min(dmat(:)); [r, c] = ind2sub(size(dmat), linIdx); i = idx(r); j = idx(c); % 删除i，概率加到j p(j) = p(j) + p(i); loss = loss + p(i) * minVal; p(i) = 0; % 更新有效索引 idx(r) = []; % 重新规范化p（其实换出来的p(i)已经是0，不需要严格归一，但保险起见） p = p / sum(p); end redIdx = idx; redProb = p(idx); redProb = redProb / sum(redProb); end

这段代码有个细节：删除索引时，如果直接idx(r)=[]，那么原场景中r位置之后的所有索引在idx里会整体前移，所以不要使用原始矩阵S的行号做后续操作。更稳妥的方式是用逻辑掩码，或者先找到最要删除的那个场景编号i，再从idx中删除idx==i的元素。否则容易产生bug。

还要注意距离矩阵D永远是原始场景之间的距离，不随删除而更新。循环里每次从剩余索引中取出子矩阵找最近对，复杂度大约是O((N-K)N²)。K比较小时，重复计算子矩阵和min查找并不经济。更高效的实现是在每轮更新距离矩阵自身，但当N不大时，子矩阵min足够用。

3.4 前向选择快速实现思路

如果决定用前向选择，代码略有不同。常见的快速前向选择算法需要维护原始场景集到保留场景集的“最优传输”关系，复杂度高。简化版本是从空集开始，循环K次，每次遍历所有未选场景，计算把该场景加入后削减损失的最小值，复杂度O(K N²)。N=1000时通常也能接受。写法类似：

selected = []; unselected = (1:N)'; for k = 1:K bestLoss = inf; for candIdx = unselected' % 将cand加入selected，计算概率距离损失 [~, lossTmp] = calcSelectedLoss(S, p, [selected; candIdx]); if lossTmp < bestLoss bestLoss = lossTmp; bestCand = candIdx; end end selected = [selected; bestCand]; unselected(unselected == bestCand) = []; end

其中calcSelectedLoss是计算当前保留集与原始集之间的Wasserstein距离，这里需要额外实现。为了简化，很多代码直接用“最近代表”原则估算损失：

for i = 1:N distToSel = min(D(i, selected)); totalLoss = totalLoss + p(i) * distToSel; end

这种做法是后向削减的单步损失在保留集上的推广，虽然不是严格的Kantorovich距离，但作为快速估算足够用。

3.5 可视化模块怎么写

削减完不能光给个索引和概率，还要能画图。我的习惯是输出三张图：

第一张，原始场景集和削减后场景集叠加。用原始场景的浅色线做背景，削减后的代表性场景用粗线突出。第二张，削减后场景概率条形图，直接用bar(1:K, redProb)。第三张，削减前后出力均值、标准差和累积概率分布的对比曲线。

可视化模块虽然简单，但能帮你直观判断削减是否把边缘分布和时序波动特征保留住。我经常跑完之后先看均值曲线能不能重合，再看概率分布尾部是否有明显缺失，这两点通过了，优化结果才敢拿出去用。

4. 参数选择、结果评估与常见问题

4.1 目标场景数K怎么选

K是最关键的参数。它不是一个需要过拟合的数字，而是取决于下游问题的可计算性。两阶段随机优化里，场景数每翻一倍，对偶变量的规模都会增加，所以K要控制在可解范围内。

一个比较实用的做法是画“削减损失-场景数”曲线。横轴是K，纵轴是削减后的概率距离损失。把K从2一直试到30，计算每个K下的loss，你会看到曲线先快速下降，到某个值后变平缓。选择“肘部”附近的K，比如从10到15，既保留了主要概率信息，又不至于让优化模型爆炸。

我之前在某个风电场出力随机调度项目里，N=500条初始场景，看到K=8时loss已经降到0.1以内，K=15后loss基本不下降，最后选K=12，调度模型两分钟能跑完。如果没有这条曲线，我可能随手选K=5，结果就会偏保守。

4.2 概率距离度量的选择与影响

很多人直接默认用欧氏距离，但不同场合有不同考量。

实践里我遇到过用欧氏距离削完以后，保留场景都是“平典型”的出力曲线，没有波动剧烈的场景。后来把距离改成“时段加权”：白天光伏时段权重高，夜晚风电时段权重高，削减出来的场景更符合调度关心的关键时段。这个技巧在论文里可以写成“时段自适应距离”。

4.3 常见问题与排查速查表

我在跑这个代码的过程中踩过几个坑，整理成了一张速查表：

另外一个容易被忽视的坑是：如果场景是由蒙特卡洛方法生成，初始概率不是等概率，而是和抽样频率有关。削减时要正确传入p向量，不要默认等概率。当时我直接用等概率削，结果某个区域概率被整体低估，优化结果明显偏离实际情况。

4.4 如何验证削减结果是否有效

削减完事不能只看损失数值，还要看下游应用的稳定性。我一般做两个验证：

第一是统计验证。把原始场景集和削减后场景集各自算均值、方差、分位数、时序相关系数，放在一张表里对比。如果均值误差在1%以内，方差误差在5%以内，基本合格。第二是优化验证。用一个固定的决策问题，分别用原始场景集和削减后场景集求解，比较目标函数最优值的偏差。偏差小说明削减没有把决策偏好带偏。做这个验证确实多花时间，但能让你在发表结果时有底气。

有一次我甚至发现，直接用原始1000个场景跑出来的目标函数值和削减后10个场景跑出来的差不到2%，但求解时间从一小时降到一分钟。这种对比图放在报告里非常有说服力。

5. 双风、风水联合场景削减的实际扩展

5.1 联合削减与分别削减的差别

实际系统里大概率同时有风电场和光伏电站，场景对象就不是单条功率曲线，而是两条时序“绑在一起”的联合场景。两种处理办法：第一种，对风、光分别削减，然后做笛卡尔积组合。这种办法麻烦在于如果风电保留8个场景，光伏保留8个场景，组合起来就是64个场景，数量呈平方增长。第二种，直接把风、光时序拼接成一个向量，每行S_i = [风功率序列，光伏功率序列]，再做一次联合削减。这样保留场景直接的联合相关性不会被破坏，而且场景数量可控。

我从实际项目里的体会是，联合削减更适合风、光互补性强的区域。比如白天光伏有出力、晚上风电大，用联合削减可以保证场景集中的“早晚转换”形态仍在，而分别削减可能丢掉风、光之间此消彼长的关系。

5.2 削减结果怎么与随机优化衔接

削减后的场景集和概率可以直接作为随机规划的输入。以两阶段调度模型为例，第二阶段变量按场景分别建，目标函数写成期望形式：

sum(redProb(k) * obj_k(redIdx(k)))

这里redProb和redIdx正是削减函数返回的结果。注意场景数K要匹配优化模型可解范围。如果追求鲁棒性，可以把削减后的场景作为不确定集合的基点，再在周围生成一个凸包做鲁棒优化，能兼顾概率信息和极端场景。

如果还要进一步做CVaR评估，削减后的场景集同样适用：取每个场景的损失值，按概率从大到小排序，再计算尾部平均。这里更考验削减算法是否保留尾部场景，所以前文提到的“不要把尾部削掉”就格外重要。

5.3 一个完整的迷你案例

为了让你直观看到代码怎么用，我写一个最小可复现案例。假设已有风、光两列数据，构建S矩阵后调用削减函数：

% 假设windData: N-by-T, pvData: N-by-T S = [windData, pvData]; % 横向拼接 N = size(S, 1); p = ones(N, 1) / N; % 等概率 K = 10; [redIdx, redProb, ~] = scenarioReduction(S, p, K, 'euclidean'); % 可视化 figure; hold on; for k = 1:length(redIdx) plot(S(redIdx(k), :), 'LineWidth', 1.5); end title('削减后的代表场景（风电+光伏时序拼接）');

这个例子合并了风、光的时序，直接对联合场景削减。实际输出里你可以明显看到保留场景在形态上有差异，不会出现冗余。

6. 最后说一些实在的体会

这套基于概率距离快速削减法的MATLAB代码，我在不同项目中反复用了很多次，最大的感受是：场景削减算法本身不是瓶颈，真正影响结果的是“场景怎么生成”和“距离怎么定义”。生成端如果本身质量就差，削减算法再精巧也救不回来；距离定义如果不贴合业务关注点，削减结果就是“距离上最优但业务上不佳”。所以拿到代码不要急着抄完就跑，先把自己的数据摸底，看看有哪些时段是你真正关心的，再决定是否对距离加权，以及K取多少。

另外，如果你以后会在小场景数量下用前向选择，也别急着优化速度，先保证逻辑对：拿一个N=5,K=2的小例子手算一遍，理解“每次选一个场景降低损失”的过程，比直接跑大矩阵出结果要可靠得多。代码里多写一点注释，把每步的公式和变量含义标清楚，后面回来看不会想“这是我写的吗”。

把削减函数和可视化、loss曲线放到一个统一入口里，遇到新数据直接一键出图，能省掉很多重复工作。这也是我推荐任何科研工具代码走的方向：先能用，再顺手，最后通用。