做新能源随机规划、或者正在被导师安排去搞微电网/配电网优化的人，十有八九绕不开一个环节：风光场景生成与削减。我在帮几个同学调代码的过程中发现，大家卡住的点几乎一模一样——理论公式能看懂，但一落到MATLAB里就不知道怎么写主循环，写出来又慢得离谱，削减完的场景分布也不对劲。今天直接把我用的“基于概率距离快速削减法”的思路、代码和踩过的坑串起来讲一遍。

这套方法说白了就解决一个问题：蒙特卡洛抽样出来的几千上万个风光出力场景，扔进随机机组组合、储能容量配置这类优化模型里，计算机根本算不动。场景削减就是把这些场景压缩成十几个、几十个，但压缩完之后，削减前后概率分布的“距离”要尽量小，也就是不能把代表性的天气模式丢掉。文章适用的人很明确：正在写小论文、做课题的硕士博士，或者实际工程里要做随机规划、时序仿真的工程师。

先说结论：概率距离快速削减法，本质是一种基于Kantorovich距离（也叫Wasserstein距离）度量的贪心后向削减算法。原理不复杂，但实现细节决定成败。

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

1.1 从“一条曲线”到“上千条曲线”

风电场、光伏电站的出力，核心特征就是随机波动。今天中午太阳好，光伏曲线是一道漂亮的钟形；明天多云，中午瞬间掉一半。做确定性分析时，我们拿一条预测曲线就进模型了；但做随机优化时，需要把这种不确定性显式建模——不能只靠“预测值加减几个百分点的误差”对付过去。

最常规的做法是蒙特卡洛抽样。先对风速、辐照度这类随机变量做概率建模，比如风速常用Weibull分布，光照用Beta分布，然后抽样生成N条完整的T时段出力曲线，每条曲线就叫一个“场景”，概率初始化为1/N。N取多少？很多论文里喜欢取1000到5000，图倒是挺好看，但后面接优化模型就麻烦大了：一个24时段、1000场景的随机机组组合，要混合整数规划求解，模型规模直接爆炸，普通电脑可能要跑十几个小时甚至直接内存溢出。

所以“场景生成”是基础，“场景削减”才是让整个链条能跑起来的关键。我见过很多人把重点全放在生成端，堆了很多精细的时序模型，到了削减环节草草用个K-means聚类交差。结果后面优化出来的方案，要么过于乐观，要么过于保守，根源就是削减环节的概率信息丢了太多。

1.2 削减的本质：逼近概率分布，不是“挑几张代表图”

刚接触场景削减的人，最容易陷入一个误区：以为削减就是“在N个场景里挑长得最像的几个留下来”。如果只做K-means聚类，聚类中心确实能代表每簇的“平均形态”，但这里有个致命问题——K-means只最小化几何距离的平方和，它不保证削减后的离散概率分布在“概率距离”意义下逼近原分布。

打个比方：你要给一群人的身高做摘要。K-means的做法是把人按身高分几组，然后报告每组的平均身高；但如果你关心的其实是“人群里有没有一米九以上的人”，K-means聚类中心很可能把这个信息抹平了。场景削减的目标不一样：它希望你用少数几个场景，尽可能还原整个概率空间的分布特征，尤其不能丢掉那些概率不高但影响很大的极端场景。这就是为什么不直接聚类，而要引入概率距离度量。

概率距离快速削减法，走的是另一条路：以概率测度之间的距离最小化为削减准则，每一步削减都明确知道“我删掉这个场景、把概率并给哪个场景，会带来多大的分布损失”。这个方法理论上干净，实现起来也比一般人想的要简单。

2. 快速削减法的数学原理与算法设计

2.1 Kantorovich距离：两个场景集合到底“像不像”

先绕不开一个名词：概率距离。这里通常用的是Kantorovich距离，也叫Wasserstein-1距离。离散概率分布场景集合之间的Kantorovich距离，可以写成一个很直观的式子：

$$ D_K(P,Q) = \sum_{i=1}^{N} p_i \cdot \min_{j=1}^{M} | s_i - t_j | $$

其中$P$是原始分布，有N个场景，每个场景$s_i$对应概率$p_i$；$Q$是削减后的分布，有M个场景$t_j$。翻译成人话就是：对于原始场景集合里的每个场景，找到削减后离它最近的场景，把这个距离乘以原始概率，最后全加起来。

这个式子的妙处在于，它把“分布差异”具体量化成了“每个原始概率质量往上最近的保留下来的场景移动，总共搬了多远”。概率质量搬运的比喻很贴切——我们削减场景，本质上就是在重新分配概率质量，分配之后质量所在的位置变了，变化越小，削减后分布越接近原分布。

还要注意一个容易被忽视的点：理想情况下Kantorovich距离应该是对称双向的，但对于“削减”这个特定问题，我们通常只计算原始分布到削减分布的“单向搬运代价”。因为在削减这个方向，我们关心的是“原始概率质量丢失了多少信息”，不需要再反向搬一遍。代码里实现的时候，直接算一句min距离求和就行。

2.2 同步回代消除法：每一步删掉“损失最小”的场景

同步回代消除（Simultaneous Backward Reduction）是快速削减法的核心实现方式。算法逻辑非常清晰：

第一步，计算所有场景两两之间的距离矩阵。对每个保留场景$i$，找到它最近的保留场景$j$，记为最近距离$d_i$。

第二步，对每个保留场景计算“如果删除它，会产生的概率距离损失”，公式是$p_i \times d_i$。为什么这个乘积就是损失？因为删除场景$i$后，它的概率$p_i$要并给离它最近的场景$j$。这个搬运动作，正是把$p_i$大小的概率质量搬运了$d_i$那么远。损失等于概率乘以距离，跟Kantorovich距离的定义完全对应。

第三步，找损失最小的那个场景，删掉它，把它的概率累加到最近的场景上。

第四步，重复循环，直到剩余场景数等于目标K。

这个算法最吸引人的地方是简单且可解释。每一步都在最小化“当下的”分布损失，所以叫贪心。贪心不保证全局最优——真正的最优削减得做大规模整数规划——但实践里效果非常接近最优，而且速度要快几个数量级。我在实际项目中默认就选这个方案，因为后向削减的每一步保持场景原始形态，不会像K-means那样生成一个“非真实存在”的合成场景。对后续做调度的人来说，场景必须是历史上确实可能出现过的出力形态，合成平均曲线没有物理意义。

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

除了后向削减，还有一类方法是快速前向选择（Fast Forward Selection）。思路反过来：从一个空集合出发，每次选一个使“当前已选集合逼近原始分布”的损失增量最小的场景加入，直到凑满K个。

两者的取舍我直接给结论：

实际工程项目里，我的经验是后向削减用得最多，因为初始场景数一般不超过5000，后向消除完全跑得动。前向选择更适合场景规模极大、后向算法内存吃紧的场合。这篇博文的代码实现以后向削减为主，同时也符合标题里“快速削减法”的普遍认知。

3. MATLAB实现：从数据生成到削减完成的全流程

3.1 前置准备：分布拟合与初始场景生成

场景削减的前提是先有“待削减的场景集合”。如果你手里是某风电场的实测历史数据，可以直接用经验分布；但更多时候手头只有风速、光照的统计参数，那就得先做蒙特卡洛抽样。

以风电为例：风速用两参数Weibull分布建模。

% 风速场景生成参数 % 某风电场风速统计拟合结果：形状参数k=2.1，尺度参数c=8.5 m/s k = 2.1; c = 8.5; N = 1000; % 初始场景数 T = 24; % 一天24个时段的出力序列 % 生成 N 个 T 维风速场景 wind_speed = wblrnd(c, k, N, T);

这里有第一个坑：你要是直接把风速场景拿去削减，削减出来的“场景”是风速曲线，不能直接用于电力调度。必须先把风速通过风机功率曲线，转换成出力标幺值序列。

% 简化风机功率曲线转换 % 切入风速3m/s，额定风速12m/s，切出风速25m/s，额定出力1.0 v_in = 3; v_r = 12; v_out = 25; P_max = 1.0; power = zeros(N, T); for t = 1:T v = wind_speed(:, t); p = zeros(N, 1); p(v < v_in) = 0; idx_run = v >= v_in & v < v_r; p(idx_run) = P_max .* (v(idx_run).^3 - v_in^3) / (v_r^3 - v_in^3); p(v >= v_r & v < v_out) = P_max; p(v >= v_out) = 0; power(:, t) = p; end

光伏的处理同理，用Beta分布抽样辐照度，再乘上光转电效率。这里要注意：风速或辐照度场景和出力场景，本身可以各自生成一套削减。但如果你要研究风光联合出力的随机性，就需要处理两者之间的相关性，这放在后面的扩展方向里细聊。当前先把独立生成这条链路跑通。

场景生成完毕后，所有场景的概率初始化为：

p = ones(N, 1) / N;

3.2 削减主循环：MATLAB核心代码

这是整篇博文的主角。我直接给出后向削减的核心函数，这个版本是我在大量调试之后确定的稳定形态：距离矩阵预计算一次，用逻辑掩码标记活跃场景，避免反复构造临时数组。

function [S_red, P_red] = reduce_scenarios(S, P, K) % S: N x D 场景矩阵，每行一个场景 % P: N x 1 概率向量 % K: 目标场景数量 % 返回削减后的场景矩阵与概率向量 N = size(S, 1); % 一次预计算所有场景两两欧氏距离 % 对角线置为Inf，避免把自己当作最近邻 D = pdist2(S, S); D(1:(N+1):end) = Inf; % 活跃标记：true表示场景仍保留 active = true(N, 1); while sum(active) > K % 获取当前活跃场景的索引 idx_active = find(active); n_active = length(idx_active); % 对每个活跃场景，找它距离最近的活跃场景 min_dist = Inf(N, 1); nearest_j = zeros(N, 1); for i = idx_active' % 提取D中第i行、活跃列的元素 d_row = D(i, active); % 找最小距离和对应位置 [md, j_loc] = min(d_row); min_dist(i) = md; % 将局部索引转换成全局场景索引 nearest_j(i) = idx_active(j_loc); end % 计算“删除损失”：概率乘以最近距离 loss = P .* min_dist; % 非活跃场景不参与比较 loss(~active) = Inf; % 选出损失最小的场景删除 [~, i_del] = min(loss); j_keep = nearest_j(i_del); % 概率合并：把被删场景的概率累加到保留场景上 P(j_keep) = P(j_keep) + P(i_del); % 标记删除 active(i_del) = false; end % 提取结果 S_red = S(active, :); P_red = P(active); % 浮点累计可能造成轻微偏差，归一化保证概率和为1 P_red = P_red / sum(P_red); end

这个函数有几个细节必须说明白。

第一个，距离矩阵为什么要在循环外算一次而不是每轮重新算？因为场景坐标在削减过程中始终不变，距离矩阵是固定不变的；变的只是“哪些场景还活着”。每轮重新算距离矩阵，白白浪费O(N²D)的时间，N=3000时一轮就要算很久。而我们只需要在每轮遍历活跃场景时，屏蔽掉已删除的列就行。

第二个，为什么用D(i, active)这个索引方式？MATLAB里向量索引传逻辑向量，返回的是true位置对应的元素。这个写法既简洁又不会越界。有人喜欢用idx_active索引子矩阵，每轮重构，也可以，但频繁构建D(:, idx_active)会产生很多临时拷贝，内存开销大，场景多的时候性能有明显差距。

第三个，删除准则用P .* min_dist，这里体现的是“概率距离损失”而不是单纯的距离。初学者最容易漏掉这个概率加权——如果你只用min_dist选最小值，等于默认所有场景概率相等，但削减几轮后概率已经不均等了，权重必须参与决策。

3.3 削减效果的评估与可视化

削减完不能直接拿去用，先评估一下削得好不好。核心指标就是Kantorovich距离下降曲线。我习惯这样算削减前后的距离：

function KD = kantorovich_dist(S1, P1, S2, P2) % 计算原始分布 P1 到削减分布 P2 的单向Kantorovich距离 % S1: 原始场景矩阵 N1 x D % S2: 削减场景矩阵 N2 x D D = pdist2(S1, S2); minD = min(D, [], 2); % 每个原始场景到削减集合的最近距离 KD = sum(P1 .* minD); % 概率加权求和 end

运行完削减函数后，把这几个数打出来：

k0 = kantorovich_dist(S, p, S, p); % 原始分布与自身的距离，正常为0 k1 = kantorovich_dist(S, p, S_red, P_red); % 削减后分布距离 fprintf('原始场景数: %d -> 削减后场景数: %d\n', N, length(P_red)); fprintf('削减前后Kantorovich距离: %.4f\n', k1);

这个距离不是绝对意义的误差，最好做横向对比：同一组原始场景，分别削减到K=5、10、20、50，看距离怎么随K下降。这本身就是一张很能说明问题的论文图。

可视化方面，我会画三个东西：

第一，削减前后全天均值曲线对比。把削减后场景按概率加权求期望，得到一条期望出力曲线，和原始1000个场景的期望曲线放在一起，肉眼能看出是否吻合。

第二，分位数带对比。画原始场景的5%和95%分位数带，再把削减后的场景散点叠上去，能直观看到削减场景是否覆盖了不确定性范围。

第三，概率柱状图。削减后的每个场景概率条形图，检查是不是有某个场景概率异常高（比如超过0.5），如果有，说明削减得太狠或者初始场景本身有很强的支配模式，需要留意。

4. 关键参数选择与完整算例

4.1 初始场景数N和目标场景数K怎么定

这两个参数是整个流程里最影响结果的上层选择。

初始场景数N，我一般按“维度和分布复杂度”来定。24时段的单风电场出力，1000个初始场景基本够用；如果你要刻画的地区地形复杂、风速多峰分布，可以上到3000到5000。再多就没什么收益了，因为概率距离削减的效果在N超过一定值后会收敛，徒增计算时间。

目标场景数K怎么取？这里没有万能公式，但有一个很实用的工程经验：K的取值取决于你后面要解的优化模型能承受多少整型变量。随机机组组合里，每个场景会乘上机组启停变量，K翻一倍，整数变量规模就翻一倍。所以先看求解器能扛多少，再反推K。比如我的经验是，商用求解器在普通工作站上，带10000个左右整数变量就是舒适上限，如果机组数有10台、时段24个，那K取到10已经接近规模极限了。

还有一个辅助判断方法：画K从2到30的Kantorovich距离曲线。曲线一般在某个K值之后变平，再增大K距离下降很慢，说明边际收益低了。选那个“肘部”附近的K值，既保留分布信息又不至于撑爆模型。

4.2 一个完整的风电场场景削减算例

拿一个我调试过的模拟风电场景举例（数据脱敏，仅演示方法流程）。某风电场的风速Weibull参数k=2.1、c=8.5，切入风速3m/s、额定风速12m/s、切出风速25m/s，额定容量50MW。生成1000个24时段出力场景，目标削减到10个。

核心调用就三行：

S = power; % 1000 x 24 场景矩阵 p = ones(1000,1) / 1000; [S_red, P_red] = reduce_scenarios(S, p, 10);

跑完的结果具有代表性：削减后的10个场景里，大约6个是“全天小风”或者“后半夜大风”这类常规模态，加起来概率占了0.72；剩下4个是“中午突发大风”“傍晚风速急升”这类带有明显时序特性的过渡模态，概率总和0.28。这正是我们希望看到的效果——削减没有把所有概率压在一个“平均形态”上，反而保留了分布尾部的关键模态。

期望出力曲线的平均绝对误差大概在0.03到0.05之间（标幺值），也就是说削减后10个场景的加权期望，和原始1000个场景的期望曲线偏差不超过5%。这个误差水平对大多数随机规划任务都是可接受的。

运行时间上，我用的是普通笔记本，1000个场景、24维、削减到10，主循环大约跑4到8秒。这个速度完全够用，配合参数扫描做敏感性分析也不心疼。

4.3 概率初始化与浮点归一化细节

场景概率初始化为均等1/N，这个做法物理意义清楚：蒙特卡洛抽样中每个样本本来就有相同的被抽中概率。削减过程中的概率合并用累加，比如被删场景概率是0.001，合并后保留场景概率从0.002变成0.003，只要最后整体归一化，不影响相对大小。

但我必须提醒一个容易踩的坑：循环几十轮之后，由于浮点数的有限精度，所有概率加总可能不再是严格的1，偏差在小数点后第15位左右，本来无所谓。但如果你后面要用概率作为约束或目标函数的权重，归一化就是必须的。所以主函数最后那行P_red = P_red / sum(P_red);不是可有可无，是保命操作。

还有个小细节：概率合并时如果两个场景本身物理上非常接近（比如距离小于1e-6），合并后会出现两个几乎重复的保留场景，白白浪费一个场景名额。这种情况在初始场景里冗余特别多时会出现，解决办法是在削减前先做一次简单的“去重”预处理，把距离小于阈值的场景先合并掉。

5. 常见问题与调试实录

5.1 削减代码跑得很慢怎么办

这篇文章的标题叫“快速削减法”，但如果你直接照搬网上某些教程的版本，跑起来可能一点都不快速。最常见的原因是没有预计算距离矩阵，而是每轮循环里重新算一遍所有场景的距离。这等于把本来该做一次的事重复做了N-K次，时间复杂度从O(N²D + N·K)变成O(K·N²D)，N=1000时还凑合，N=3000时直接让人崩溃。

优化思路有三步：

第一步，把距离矩阵计算挪到循环外，预计算一次D = pdist2(S, S)。pdist2是多维矩阵两两距离计算的最快内置函数，远好过自己写双层for循环套欧氏距离公式。

第二步，在主循环内部，不要再对全量N做运算，尽量通过逻辑索引只触碰活跃场景。比如D(i, active)天然只取活跃列，运算量随活跃场景数递减。

第三步，如果N实在太大，可以对最内层遍历使用parfor并行。但要小心：并行循环里对min_dist和nearest_j的赋值，最好是每个worker只写自己负责的索引位，避免数据竞争。我的经验是一般N不超过5000时没必要上并行，上了反而有开销。

5.2 削减后出现NaN或者概率变成负数

出现NaN，百分之九十是距离矩阵对角线的问题。如果某一行距离全是无穷大，min的结果自然是Inf，Inf乘以0得到NaN。什么样的行会全是无穷大？当场景矩阵里有重复行，且对角线被置为Inf后，重复场景之间的真实距离也是0，但如果你的代码里连非对角也误置为Inf了，就会出现整个活跃集合里找不到近邻的情况。

排查方法很直接：削减前打印min(D(i, :))，看有没有行全部是Inf；削减过程中在while循环里加断言assert(all(P >= 0))。我调试的时候一般把N调小为50，打印每一步的删除索引、保留索引和合并后的概率，一行行盯，问题马上就暴露了。

概率变负数的情况很少见，但如果你的概率累加逻辑写成了P(j_keep) = P(j_keep) - P(i_del)，那不出三轮就会出现负概率。记住这里永远是加号：被删场景的概率转移给邻居，不删除任何概率质量。

5.3 削减完场景“挤成一团”，极端场景丢了

这是另一种常见抱怨：削减后剩下的场景长得都很像，分布特征消失。根本原因通常是K值取得太小，或者初始场景本身就缺乏多样性。

K=5的时候，算法必然优先保概率质量集中的区域，尾部极端场景的单个概率太小，删掉它们带来的损失也小，所以会被优先淘汰。这其实是正确的数学行为——概率距离准则在K过小时本来就会牺牲尾部。解决办法不是改代码，而是增加K，或者改“强制保留极端场景”的带约束削减。

我在实际项目里用过一种简单有效的变体：先在原始场景里选出若干“极端场景”（比如最大出力超过额定90%的场景、连续低出力超过6小时的场景），把这些场景直接从候选删除池里豁免，只对所有非豁免场景做概率距离削减，最后把豁免场景加回结果集。这种做法在储能容量配置问题里很管用，因为储能容量往往由极端场景决定，你把尾部削没了，容量配置就容易偏小。

5.4 独立抽样忽略了风光的耦合关系

如果做的是风光联合场景，独立抽样会产生一个问题：自然界的风电和光伏出力在同一时段有相关性——比如多云天气光伏低、同时风速往往较大。独立抽样会让两个随机变量各抽各的，把耦合关系抹掉。

处理方案有两条路。工程上最简单的是用历史场景数据做“再采样”：从历史同期数据里直接抽取时段曲线，天然带相关性。学术上更讲究的可以用Copula函数建模相关性，先生成相关的风速和光照随机数，再转换出力场景。但要注意：无论哪种方式，生成的场景矩阵到了削减阶段，逻辑完全不变——削减算法不关心场景是怎么生成的，它只认场景矩阵和概率向量。所以如果你只是想跑通削减流程，先独立抽样没问题；如果要做严谨研究，把相关性处理放进生成端即可。

6. 几个值得继续尝试的扩展方向

6.1 多阶段场景树与动态削减

这篇博文讲的是一次性削减，适合单阶段随机规划。但实际调度往往要处理多阶段决策，比如日前计划加日内滚动修正。这时需要的是场景树，而不是一组平面场景。

场景树的生成方式可以递归使用削减算法：先在根节点生成大量场景，削减出一组子场景；再把时间轴切半，对每个子场景的未来时段继续生成和削减。代码复用度很高，把本博文的reduce_scenarios封装成一个子函数，对每个节点的场景矩阵分别调用就能搭出树结构。唯一要留心的是父节点的概率要沿着树的分叉逐级相乘，保证叶节点概率总和还是1。

6.2 和其他削减方法做横向对比

有人会问：既然K-means那么快，为什么不用K-means？我做对比实验时的结论是：K-means的优点是极快，缺点是聚类中心不是原始场景，且对初始聚类中心敏感；k-medoids算法保留原始场景作为中心，概率分配也算合理，但计算量比K-means大；谱聚类能处理非凸分布，但参数难调。

概率距离快速削减法在这些方法里处于很实用的平衡点：严格以概率距离为准则、保留原始场景形态、参数少且直接可解释。如果你的场景分布呈明显的多峰形态，并且你关心尾部风险，这个方法基本是首选。只有当场景数量上到十万级别，计算资源吃紧时，我才会考虑先K-means粗选一遍，再用概率距离削减做精化。

6.3 从传统削减到生成式方法的过渡

最近这两年，基于扩散模型和生成对抗网络的风光场景生成成了热门方向。这些方法能生成比蒙特卡洛抽样更加逼真的时空相关场景，但它们的输出依然是“大量样本”，最终还是要靠场景削减收敛到可计算的规模。所以本文这套削减方法在生成式方法的流程里反而变得更重要了——不管你前面用什么花哨生成器，后面接的还是这个概率距离削减主函数。

我个人在实际操作中的体会是：不要被“生成”环节漂漂亮亮的论文图迷惑，真正决定随机优化质量上限的，往往是削减环节有没有把概率信息保留住。如果你也在做相关课题，建议先把这篇文章里的核心函数跑熟，把坎托罗维奇距离的评估逻辑刻在脑子里，再去折腾更复杂的生成模型，能省掉非常多的返工时间。

最后再分享一个小技巧：削减后的场景集合，如果不是马上用，最好连同概率一起存成mat文件存档。因为后面做敏感性分析需要重复调用优化模型，每次重新生成和削减场景既浪费时间，又引入了额外的随机性，前后结果的差异会干扰对算法本身的判断。存好了场景，后续对比实验就干净多了。