做风光出力随机性研究的同学大概率都遇到过这个问题靠蒙特卡洛抽了几千个场景精度看着不错但一套随机规划模型跑下来要么内存爆炸要么求解时间长得像年终总结。场景削减就是解决这个问题的标准动作而基于概率距离的快速削减法是我们在工程里最常用到的折中方案。下面围绕风光场景生成、场景削减、快速削减法以及整套流程在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步。我在N2000、K10的情况下实测过前向选择通常比后向消除快一个数量级。当然如果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:N1: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:N1: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:N1: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²字节。N5000时约200MBN10000时约800MBN20000时超过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(后向消除, 快速前向选择);实际使用中K10到20在日前的调度和容量配置场景里非常常见。如果场景代表的是多年气象数据时段更多、变量更多K可能需要到50甚至100。不管选多少都必须配合后面的随机规划模型跑一遍看削减后的最优解是否稳定。只盯着概率距离小不一定能保证优化结果贴近全场景模型。4.4 风光联合场景量纲不一致这个问题我在实际项目里几乎每次都遇到。风电功率标幺化一般在0到1之间光伏如果直接取辐照度或者原始功率数值可能到几百。两者拼进同一个向量后欧氏距离几乎只被光伏部分控制距离矩阵看起来算得挺热闹实际上风电形状差异完全被淹没了削减结果里风电典型场景会非常粗糙。解决办法很简单拼接场景矩阵之前把风电和光伏各自归一化再拼成联合场景或者给不同的子向量设置权重系数。比较严谨的做法是按装机容量确定权重让风、光两个分量在距离计算中的实际贡献可比。归一化不改变各自时段形状只是把单位拉齐了。4.5 固定随机种子让结果可复现场景生成通常依赖随机数。如果你不固定随机种子同一套代码跑两次生成的场景集合完全不同削出来的典型场景和权重也会有差异。做研究或者写技术报告时这是很尴尬的事情因为别人无法复现你的图表。我习惯在脚本最前面写rng(2024)或者任何一组固定整数。如果需要对不同算法做公平对比还要保证它们都用同一个场景矩阵不能各自重新生成。对于重要结果我会在MATLAB里把场景矩阵、距离矩阵、削减后的场景序号和权重一起存成.mat文件后面画图、写报告都从这份文件读取避免因为误操作导致数据漂移。4.6 一次削减结果失真排查我在某个风光储容量配置项目里遇到过一次很典型的问题削减后的概率距离明明很小但优化结果和全场景模型对不上投资成本比预期低了很多。刚开始我以为是削减条数太少把K从10调到20问题依旧。后来打印每个保留场景的权重发现光伏场景里出现了几条高权重曲线风电场景却被平均化得很严重几乎看不出典型形状。查到最后原因就是风光量纲没有统一。生成的光伏出力数值范围远大于风电距离计算被光伏主导相当于风电部分的“意见”没有参与投票。把场景矩阵重新按各自最大值归一化后权重分布立刻合理了优化结果也回到正常范围。从那之后我每次计算距离矩阵前都会先看一眼各列的数量级这个检查只需几秒钟但能避免后面整条流程白跑。从个人习惯来说我现在做场景削减一定会先画K曲线确认拐点位置再定K同时固定随机种子保存中间结果。这套流程看着多花几分钟但对后续随机规划模型的可靠性帮助非常大。概率距离快速削减法本身不是特别复杂的算法真正考验人的往往是这些藏在细节里的坑。
