做新能源随机规划、或者正在被导师安排去搞微电网/配电网优化的人十有八九绕不开一个环节风光场景生成与削减。我在帮几个同学调代码的过程中发现大家卡住的点几乎一模一样——理论公式能看懂但一落到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_{i1}^{N} p_i \cdot \min_{j1}^{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个。两者的取舍我直接给结论对比维度后向削减回代消除前向选择逐步添加适用初始场景数中等规模几百到几千大规模上万甚至更多极端场景保留能力较好因为初始场景都在候选集里较弱被选中的场景偏“中庸”实现难度简单主循环容易理解中等需要维护已选集合距离增量速度每轮都要遍历剩余场景O(N³)量级类似复杂度但初始集合增长快前面几轮较慢实际工程项目里我的经验是后向削减用得最多因为初始场景数一般不超过5000后向消除完全跑得动。前向选择更适合场景规模极大、后向算法内存吃紧的场合。这篇博文的代码实现以后向削减为主同时也符合标题里“快速削减法”的普遍认知。3. MATLAB实现从数据生成到削减完成的全流程3.1 前置准备分布拟合与初始场景生成场景削减的前提是先有“待削减的场景集合”。如果你手里是某风电场的实测历史数据可以直接用经验分布但更多时候手头只有风速、光照的统计参数那就得先做蒙特卡洛抽样。以风电为例风速用两参数Weibull分布建模。% 风速场景生成参数 % 某风电场风速统计拟合结果形状参数k2.1尺度参数c8.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:(N1):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)的时间N3000时一轮就要算很久。而我们只需要在每轮遍历活跃场景时屏蔽掉已删除的列就行。第二个为什么用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);这个距离不是绝对意义的误差最好做横向对比同一组原始场景分别削减到K5、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参数k2.1、c8.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)N1000时还凑合N3000时直接让人崩溃。优化思路有三步第一步把距离矩阵计算挪到循环外预计算一次D pdist2(S, S)。pdist2是多维矩阵两两距离计算的最快内置函数远好过自己写双层for循环套欧氏距离公式。第二步在主循环内部不要再对全量N做运算尽量通过逻辑索引只触碰活跃场景。比如D(i, active)天然只取活跃列运算量随活跃场景数递减。第三步如果N实在太大可以对最内层遍历使用parfor并行。但要小心并行循环里对min_dist和nearest_j的赋值最好是每个worker只写自己负责的索引位避免数据竞争。我的经验是一般N不超过5000时没必要上并行上了反而有开销。5.2 削减后出现NaN或者概率变成负数出现NaN百分之九十是距离矩阵对角线的问题。如果某一行距离全是无穷大min的结果自然是InfInf乘以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值取得太小或者初始场景本身就缺乏多样性。K5的时候算法必然优先保概率质量集中的区域尾部极端场景的单个概率太小删掉它们带来的损失也小所以会被优先淘汰。这其实是正确的数学行为——概率距离准则在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文件存档。因为后面做敏感性分析需要重复调用优化模型每次重新生成和削减场景既浪费时间又引入了额外的随机性前后结果的差异会干扰对算法本身的判断。存好了场景后续对比实验就干净多了。