做风光出力随机性研究的同学，大概率都遇到过这个问题：靠蒙特卡洛抽了几千个场景，精度看着不错，但一套随机规划模型跑下来，要么内存爆炸，要么求解时间长得像年终总结。场景削减就是解决这个问题的标准动作，而基于概率距离的快速削减法，是我们在工程里最常用到的折中方案。下面围绕风光场景生成、场景削减、快速削减法，以及整套流程在MATLAB里的落地写法展开，先说清楚为什么要做，再给可以直接抄走的代码。

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

1.1 为什么几千个场景不能直接扔给优化模型

风电和光伏出力的随机性很直观：风速时大时小，云层来来走走，一条24小时的出力曲线，今天和明天可能完全两样。学术界和工程界最常见的做法是用历史数据拟合出力的概率分布，或者用时间序列模型生成大量可能出现的出力曲线，这些曲线在随机规划里叫“场景”。场景越多，对随机性的描述越精细，但代价也很现实：每增加一个场景，优化模型就要多复制一组决策变量和约束条件。

举一个储能容量配置的例子，如果只考虑24个时段，规模很小；但当你把场景数加到1000个，允许每个场景独立决策充放电，模型里的变量和约束会成倍增长。即便用商业求解器，也容易出现内存占用过高、求解时间无法接受的情况。更麻烦的是，很多随机规划模型对场景数极其敏感，5000个场景和10000个场景的求解时间差距不是线性增长，而是接近指数级膨胀。所以实际工作中从来不会把蒙特卡洛抽出来的原始场景直接塞给优化模型，而是先做“场景削减”。

1.2 场景生成与场景削减是两件事

很多时候大家把生成和削减混在一起，其实这两件事分工不同。场景生成要做的是“把不确定空间铺满”：用概率分布、时间序列模型、Copula或者机器学习方法，生成一大批能覆盖可能出力状态的初始场景。场景削减要做的是“从铺满的状态里挑代表”：从原始场景集合中选出一个规模较小的子集，同时给每个保留场景重新分配概率权重，让这个小子集在概率意义上尽可能接近原始大集合。

一个比较生活化的类比是选照片。你拍了一千张旅行照片，不能把所有照片都塞进相册，只能在里面挑十张代表整趟旅行，并且不能简单等权重地放，而是要根据每组照片的相似程度和出现频率分配版面。场景削减也是这个逻辑：相似的出力曲线合并成一条，并保护好它的出现概率，最终用少量“典型场景”还原原始分布的统计特征。

从实现路径上看，场景生成的质量决定了削减的上限。如果生成阶段本身丢失了时序相关性或忽略了风光间的耦合关系，后面削得再漂亮，也只是在一个有偏的样本集合里做文章。所以成熟的做法通常是把生成和削减串成一条流水线：先生成足够多的候选场景，再做削减，最后评估削减前后统计特性是否一致。

1.3 概率距离快速削减法的定位

场景削减的方法不少，最朴素的思路是拿聚类算法硬切，比如K-means。但K-means本质上是找几何中心，它优化的是“点到中心的总距离”，并不天然保证概率质量守恒。实际中很容易出现某个簇里场景数量很多、但簇中心代表不了簇内极端情况的尴尬。随机规划领域更常用的，是一类基于概率距离的场景削减方法。

概率距离的思路，是把带权重的场景看成一个离散概率分布，两个场景集合之间的接近程度，用“把一个分布的概率质量搬运到另一个分布所需的最小代价”来衡量。搬得越省力，说明这两个场景集合越像。基于这个度量，就有了两条经典削减路线：一条是从全集开始，逐步删除“删除后损失最小”的场景，叫后向消除；另一条是从空集开始，逐步选入“加入后收益最大”的场景，叫快速前向选择。两种方法都带有贪心性质，不追求全局最优，但计算量可控、效果稳定，这也是它们在工程里被反复使用的原因。

2. 概率距离快速削减法的原理拆解

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

我先把符号定下来。假设有N个原始场景，每个场景是一条P维向量，比如24时段的风电出力曲线。场景i和场景j之间的距离 d(x_i, x_j)，可以直接用欧氏距离，也就是L2范数；也可以用绝对值和，也就是L1范数。到底选哪个，取决于你的优化目标更关心什么。如果后续要用来做电力平衡，功率偏差的平方更敏感，L2合适；如果你更关心累计电量的偏差，L1更直接，因为它对单个时段的大偏差不会过分放大。

如果场景是风光联合场景，处理方式也一样。风力功率曲线是24维，光伏功率曲线是24维，拼在一起就是一个48维向量。但这时候一定要小心量纲问题：风电功率可能在0到1之间波动，光伏辐照或者功率可能是几十的数量级，直接算欧氏距离，结果会被数值大的那部分主导。我一般会先把每个分量标准化到[0,1]，再算距离，否则削减出来的典型场景会偏向某一个电源。

有了场景距离之后，概率距离就很好理解了。假设原始场景集合中，场景i的概率权重是 p_i，削减后的保留场景集合是R，那么削减前后两个离散分布之间的概率距离可以写成：

D(S, R) = sum_i p_i * min_{j in R} d(x_i, x_j)

这个式子其实是在说：每个原始场景都把自己全部的概率质量搬到离它最近的保留场景上，搬运成本等于概率乘以距离。把所有成本加起来，就是削减带来的概率距离。这个量越小，说明保留场景越能代表原始场景。实际编程时，这个公式就是整个削减算法的核心目标函数。

2.2 后向消除：最直观的实现思路

后向消除的起点是全集，也就是所有N个场景一开始都在保留集合里。然后每一步只做两件事：先算一算删除每个场景会造成多大的概率距离损失，找到损失最小的那个删掉；再把被删场景的概率权重加到离它最近的保留场景上。重复以上过程，直到只剩K个场景。

为什么删除损失可以只算当前这一步？因为每一步都是贪心选择，只考虑局部最优。假设场景j当前的概率权重是 p_j，它到其他保留场景的最近距离是 d_min，那么一旦删除j，原本属于j的概率质量就要搬运到最近的保留场景上，这一步产生的额外运输成本就是 p_j × d_min。所以每一步只需要找这个值最小的场景删除。

我用一个整理书架的类比来描述：你把一批书看成场景，每本书都有一个内容相似度。后向消除就是反复观察哪本书和书架上的另一本书最接近，然后把这本书撤掉，并把“这本书可能被翻阅的次数”叠加到旁边那本书上。撤掉的一定是信息冗余度最高的，而不是最冷门的。这个方法实现简单，效果稳定，适合场景总数不是特别大的情况。

2.3 快速前向选择：只进不出的工程加速版

快速前向选择正好反过来，它从空集开始，先选一个最能代表全局的场景进入保留集合，然后每一步从未保留场景里挑一个“加入后总概率距离下降最多”的场景进来，直到保留K个。

关键点是怎么判断“加入后下降最多”。我在代码里维护一个向量 dmin，记录每个原始场景到当前保留集合的最近距离。当候选场景c被加进来时，原始场景i到保留集合的距离可能会从 dmin_i 变成 min(dmin_i, d(x_i, x_c))，最大下降量就是 max(dmin_i - d(x_i, x_c), 0)。所有原始场景的下降量按照概率加权求和，就是加入c带来的总收益。找出收益最大的候选场景加入，然后更新dmin，进入下一轮。

之所以叫“快速”，是因为每一步只需要O(N^2)量级的计算，总共跑K步。当K很小、N很大的时候，前向选择明显比后向消除快，因为后向消除要一路删到K，删除步数大约是N-K，在前向选择里只需要选K步。我在N=2000、K=10的情况下实测过，前向选择通常比后向消除快一个数量级。当然，如果K接近N，后向消除的步数少，会更划算一些。

2.4 场景权重为什么必须重新分配

有些初学者做削减时只会“挑场景”，把不选中的场景直接扔掉，保留的场景等权重处理。这在随机规划里是个大坑。删除一个场景不等于这个场景代表的未来不会发生，它代表的可能性必须由另一个相近场景承接，否则削减后的总概率不为1，后续优化目标里的期望成本就会系统性偏低。

举个最简单的例子：两个几乎相同的场景A和B，各自出现概率0.5。削减为一个场景时，如果只留下A并把权重设成0.5，那么期望出力只算了一半；正确做法是把B的概率也转移给A，让A的权重变成1.0。后向消除里，这个转移动作发生在删除场景时的“就近合并”；快速前向选择里，所有原始场景最后都会归属到离它最近的保留场景，再把概率权重累加起来。所以说，场景削减不只是选点，还包含一个“概率质量再分配”的过程。

3. MATLAB实现：从原始场景到削减场景

3.1 算例设置与场景生成

先搭一个最小算例，方便你自己跑通。假设我们有24个时段，初始生成500个风光联合场景，目标削减到10个。这个规模在普通笔记本上运行非常快，后向消除和前向选择都只要几秒到十几秒。

我用逆变换法按威布尔分布生成风速场景，然后用典型风电功率曲线把风速转成功率。光伏部分用贝塔分布模拟归一化辐照，叠加一个白天的形状系数。这段代码只是为了验证削减算法，不是精细化建模；你把场景矩阵换成自己项目里生成的风光数据，完全不影响后面削减部分的使用。

rng(2024); N = 500; H = 24; K = 10; % 1. 风速场景：威布尔分布，逆变换采样 kShape = 2.0; AScale = 8.0; v = AScale * (-log(rand(N, H))).^(1/kShape); % 2. 风电功率曲线 vci = 3; % 切入风速 vr = 12; % 额定风速 vco = 25; % 切出风速 Pr = 1; % 额定功率标幺化 windPower = zeros(N, H); mask1 = v >= vci & v < vr; mask2 = v >= vr & v < vco; windPower(mask1) = (v(mask1) - vci) ./ (vr - vci); windPower(mask2) = Pr; % 3. 光伏出力场景：贝塔分布 + 白天形状系数 % 如果没装统计工具箱，可以换成任意辐照采样模型，这里重点是削减 betaA = 2.0; betaB = 1.8; solarRand = betarnd(betaA, betaB, [N, H]); hours = 1:H; dayFactor = max(0, sin((hours - 6) / 24 * 2 * pi)); solarPower = solarRand .* dayFactor; % 4. 拼接成联合场景，并做量纲归一化 scen = [windPower, solarPower]; scen = scen ./ max(scen(:));

如果要做更贴近工程场景的生成，建议在风速里加入时间序列相关性，比如AR(1)模型，或者在风光之间引入Copula来描述耦合关系。对于纯场景削减算法来说，输入矩阵只要满足“每一行是一个场景、每一列是一个时段变量”就够了。

3.2 距离矩阵计算

距离矩阵是削减算法的核心输入。先写一个通用的欧氏距离矩阵函数，输入场景矩阵S，输出N×N的距离矩阵，对角线设为Inf，避免削减时把场景匹配给它自己。

function D = computeDistMatrix(S) N = size(S, 1); D = zeros(N, N); for i = 1:N diff = S - S(i, :); D(i, :) = sqrt(sum(diff.^2, 2)); end D(1:N+1:end) = inf; end

如果你更关心累计电量偏差，可以把里面的sqrt(sum(diff.^2, 2))换成sum(abs(diff), 2)，这就是L1距离。还有一个常见需求是计算距离前先标准化场景矩阵，尤其对风光联合场景。如果不在外面做标准化，建议在函数入口增加一个开关，但要注意调用方后续也要用同一套标准化，否则权重分配和评估会乱。

这个函数的时间复杂度是O(N^2×P)，N是场景数，P是向量维度。N在几千以内问题不大；N到一两万，内存和时间都会非常吃紧。后文会专门讲怎么处理大场景。

3.3 后向消除法主函数实现

后向消除的输入是距离矩阵D和保留场景数K。我给出一个便于理解的循环版本，再给一个矩阵化加速思路。循环版每一步都会扫描当前活跃场景，算法逻辑清晰，适合用来对照原理。

function [redIdx, redProb] = backwardReduction(D, K) N = size(D, 1); active = true(N, 1); p = ones(N, 1) / N; while sum(active) > K losses = inf(N, 1); act = find(active); % 计算删除每个活跃场景造成的概率距离损失 for j = act' drow = D(j, :); drow(~active) = inf; drow(j) = inf; losses(j) = p(j) * min(drow); end % 删除损失最小的场景 [~, del] = min(losses); active(del) = false; % 把被删场景的概率转移给最近的活跃场景 drow = D(del, :); drow(~active) = inf; [~, nearest] = min(drow); p(nearest) = p(nearest) + p(del); p(del) = 0; end redIdx = find(active); redProb = p(redIdx); redProb = redProb / sum(redProb); end

每一步里，losses(j) = p(j) * min(drow)就是删除j带来的概率距离增量的贪心估计。删除后必须立刻把概率转移给最近邻，否则后续步骤里的p(j)会漏掉已经合并过的概率质量。

如果想加速，可以减少内层循环，把“找损失最小”那一步变成矩阵运算：每次从D(active, active)中取子矩阵，对角线置Inf，对每一行取最小值，再乘以对应权重。实现起来更紧凑，但需要小心active索引位置。我一般先用循环版调通逻辑，再视场景规模决定是否优化。

3.4 快速前向选择法主函数实现

快速前向选择在工程里更适合大批量场景，核心是维护dmin这个“每个原始场景到保留集合的最近距离”。这里有一个非常容易踩的坑：距离矩阵对角线一开始通常是Inf，但保留场景本身到自己的距离应当视为0。我在函数里先复制了一个Dloc，并把对角线清零，这样初始选择、增益计算、最后概率分配都能保持一致。

function [redIdx, redProb] = fastForwardSelection(D, K) N = size(D, 1); p = ones(N, 1) / N; % 自身距离清零，供前向选择和权重分配使用 Dloc = D; Dloc(1:N+1:end) = 0; R = false(N, 1); C = true(N, 1); % 初始选一个总距离最小的场景 [~, j0] = min(sum(Dloc, 1)); R(j0) = true; C(j0) = false; dmin = Dloc(:, j0); % 逐步加入场景 for k = 2:K gain = -inf(N, 1); for c = 1:N if C(c) improved = max(dmin - Dloc(:, c), 0); gain(c) = p' * improved; end end [~, add] = max(gain); R(add) = true; C(add) = false; dmin = min(dmin, Dloc(:, add)); end redIdx = find(R); % 按最近保留场景累加概率权重 assign = zeros(N, 1); for i = 1:N [~, pos] = min(Dloc(i, redIdx)); assign(i) = redIdx(pos); end redProb = zeros(K, 1); for r = 1:K redProb(r) = sum(p(assign == redIdx(r))); end redProb = redProb / sum(redProb); end

这里的improved = max(dmin - Dloc(:, c), 0)是整个函数的关键。它表达的是：加入候选场景c后，原始场景i到保留集合的距离会下降多少。如果下降量是正数，说明这个候选场景确实补上了当前代表性不足的区域；如果下降量为0，说明它对那些场景几乎没有贡献。把所有下降量乘上概率权重再求和，就是c的“增量收益”。

一个容易出错的细节是，候选场景c本身在加入后，它到保留集合的距离应从原来的某个值变成0，所以必须保证Dloc对角线是0，否则max(dmin - Dloc(:,c), 0)会在c这一行被Inf截断，导致c永远拿不到应有的增益，选出来的场景会整体偏向初始集合。

3.5 主流程串联与结果评估

把前面几个函数串起来，就是一个完整的“生成—削减—评估”流程。评估时最直接的两个指标：一个是削减后场景的加权均值曲线与原始均值是否接近；另一个是削减前后概率距离是否足够小。

% 生成场景矩阵后，计算距离并削减 D = computeDistMatrix(scen); [ridxB, pB] = backwardReduction(D, K); [ridxF, pF] = fastForwardSelection(D, K); % 概率距离评估函数 function dRed = evalReduction(D, redIdx, p) N = size(D, 1); Dloc = D; Dloc(1:N+1:end) = 0; dRed = 0; for i = 1:N dRed = dRed + p(i) * min(Dloc(i, redIdx)); end end dRedB = evalReduction(D, ridxB, ones(N,1)/N); dRedF = evalReduction(D, ridxF, ones(N,1)/N); % 加权平均出力曲线对比 meanOrig = mean(scen, 1); meanRedB = pB' * scen(ridxB, :); meanRedF = pF' * scen(ridxF, :);

如果meanRedB和meanRedF与meanOrig的趋势一致，峰值时段偏差不大，说明削减结果可用。我一般还会再对比几个分位数，比如5%和95%出力曲线，因为均值曲线接近不代表极端场景保留得好。分位数对比对后续鲁棒优化更重要。

3.6 削减效果可视化

代码跑完了，一定要画图，图比任何数值指标都直观。我习惯画三样东西：原始场景集合、削减后的典型场景、原始均值与削减后加权均值的对比。

figure('Color', 'w'); subplot(2, 1, 1); plot(scen', 'Color', [0.82 0.82 0.82]); hold on; plot(scen(ridxB, :)', 'LineWidth', 1.2); title('原始场景与后向消除保留场景'); subplot(2, 1, 2); plot(meanOrig, 'k--', 'LineWidth', 2); hold on; plot(meanRedB, 'r-', 'LineWidth', 2); plot(meanRedF, 'b-', 'LineWidth', 2); legend('原始期望', '后向消除期望', '前向选择期望', 'Location', 'best');

也可以单独画一下保留场景的概率权重柱状图，看权重分配是不是合理。如果某个保留场景权重特别高，说明它吸收了很大一片场景的概率质量；如果权重都特别平均，而原始场景分布又明显不均匀，就要怀疑是不是削减算法没有按概率合并。

4. 实操踩坑与常见问题排查

4.1 削减后期望曲线对不上：权重归一化

最常见的问题就是加权期望曲线整体偏低，甚至总出力曲线明显“缺了一块”。原因多半是削减时只删了场景，没有转移概率，或者转移了概率后忘记做归一化。你要记住，削减不是筛选，是合并。每删除一个场景，它的概率权重必须合并到最近的保留场景上。在后向消除的代码里，p(nearest) = p(nearest) + p(del)这一行不能省；在快速前向选择里，最后必须按最近保留场景重新累加原始权重。

我习惯在每一步之后加一个校验：sum(active) 对应的总概率是不是恒等于1。如果发现概率不为1，基本上就是转移逻辑出问题了。另外，有些场景削减方法会故意保留常数权重，那是另一种建模风格，但和概率距离削减的初衷不一致，不能混用。

4.2 距离矩阵内存爆炸与减速

距离矩阵是N×N的double数组，占用的内存大约是8×N²字节。N=5000时约200MB，N=10000时约800MB，N=20000时超过3GB。很多人在生成上万场景后直接跑computeDistMatrix，会发现电脑卡成PPT。

我的建议是，初始场景数控制在1000到3000之间，对大多数风光出力不确定性建模已经够用。如果实在需要更大样本，可以用分块计算。快速前向选择每轮只需要某个候选场景到所有场景的距离列，完全可以按块计算，没必要把整个D常驻内存。后向消除也可以用同样的思路，但实现会麻烦一些。更省事的折中是先对原始场景做一次粗筛，把高度重复的场景合并掉，再进削减流程，样本量能迅速降下来。

4.3 保留场景数K怎么选

K不是越大越好，也不是固定拍脑袋。正确做法是画“概率距离随K变化的曲线”。随着K增大，概率距离会单调下降，但下降速度会越来越慢。曲线出现明显拐点的位置，通常就是精度和计算量的平衡点。

Ks = 5:2:50; errsB = zeros(size(Ks)); errsF = zeros(size(Ks)); pUniform = ones(N,1) / N; for t = 1:numel(Ks) k = Ks(t); [~, pB] = backwardReduction(D, k); [~, pF] = fastForwardSelection(D, k); errsB(t) = evalReduction(D, ridxB, pUniform); errsF(t) = evalReduction(D, ridxF, pUniform); end plot(Ks, errsB, 'o-', Ks, errsF, 's-'); legend('后向消除', '快速前向选择');

实际使用中，K=10到20在日前的调度和容量配置场景里非常常见。如果场景代表的是多年气象数据，时段更多、变量更多，K可能需要到50甚至100。不管选多少，都必须配合后面的随机规划模型跑一遍，看削减后的最优解是否稳定。只盯着概率距离小，不一定能保证优化结果贴近全场景模型。

4.4 风光联合场景量纲不一致

这个问题我在实际项目里几乎每次都遇到。风电功率标幺化一般在0到1之间，光伏如果直接取辐照度或者原始功率，数值可能到几百。两者拼进同一个向量后，欧氏距离几乎只被光伏部分控制，距离矩阵看起来算得挺热闹，实际上风电形状差异完全被淹没了，削减结果里风电典型场景会非常粗糙。

解决办法很简单：拼接场景矩阵之前，把风电和光伏各自归一化，再拼成联合场景，或者给不同的子向量设置权重系数。比较严谨的做法是按装机容量确定权重，让风、光两个分量在距离计算中的实际贡献可比。归一化不改变各自时段形状，只是把单位拉齐了。

4.5 固定随机种子让结果可复现

场景生成通常依赖随机数。如果你不固定随机种子，同一套代码跑两次，生成的场景集合完全不同，削出来的典型场景和权重也会有差异。做研究或者写技术报告时，这是很尴尬的事情，因为别人无法复现你的图表。

我习惯在脚本最前面写rng(2024)或者任何一组固定整数。如果需要对不同算法做公平对比，还要保证它们都用同一个场景矩阵，不能各自重新生成。对于重要结果，我会在MATLAB里把场景矩阵、距离矩阵、削减后的场景序号和权重一起存成.mat文件，后面画图、写报告都从这份文件读取，避免因为误操作导致数据漂移。

4.6 一次削减结果失真排查

我在某个风光储容量配置项目里遇到过一次很典型的问题：削减后的概率距离明明很小，但优化结果和全场景模型对不上，投资成本比预期低了很多。刚开始我以为是削减条数太少，把K从10调到20，问题依旧。后来打印每个保留场景的权重，发现光伏场景里出现了几条高权重曲线，风电场景却被平均化得很严重，几乎看不出典型形状。

查到最后，原因就是风光量纲没有统一。生成的光伏出力数值范围远大于风电，距离计算被光伏主导，相当于风电部分的“意见”没有参与投票。把场景矩阵重新按各自最大值归一化后，权重分布立刻合理了，优化结果也回到正常范围。从那之后，我每次计算距离矩阵前都会先看一眼各列的数量级，这个检查只需几秒钟，但能避免后面整条流程白跑。

从个人习惯来说，我现在做场景削减一定会先画K曲线，确认拐点位置，再定K；同时固定随机种子，保存中间结果。这套流程看着多花几分钟，但对后续随机规划模型的可靠性帮助非常大。概率距离快速削减法本身不是特别复杂的算法，真正考验人的，往往是这些藏在细节里的坑。