做风光出力场景的人多半都撞过同一堵墙蒙特卡洛随便一抽就是上千个场景扔进随机优化模型里求解器直接被迭代次数和约束规模拖垮。我在某园区微网的日前调度项目里第一次遇到这个问题时光是处理一台风机和一个光伏阵列的联合场景就跑了一整晚没出结果。后来把场景数从一千砍到几十计算时间从小时级掉到分钟级优化结果却几乎没有肉眼可见的差异——关键就在场景削减这一环。这篇文章就围绕“基于概率距离快速削减法的风光场景生成与削减”展开我会讲清楚场景生成怎么建模、概率距离削减的原理为什么站得住、MATLAB代码怎么一步一步写出来以及削减效果如何量化评估。内容主要面向正在做电力系统随机优化、微电网能量管理、风光出力不确定性分析的研究生和工程师如果你只是想把场景削减模块接进自己的项目里也可以直接抄后面的代码和参数设置。1. 场景削减这件事为什么绕不开1.1 随机优化里的“场景风暴”问题风光出力的随机性是个客观事实不是靠提高预测精度就能完全消灭的。无论是风电的湍流强度、光伏的云层遮挡还是负荷侧的随机波动都需要用一组带概率的离散场景去近似连续的概率分布。蒙特卡洛抽样是最直观的做法对风速、辐照度等随机变量做大量抽样每个样本生成一条完整的时间序列就得到一个场景。问题恰恰出在“大量”上。假设风电出力用24小时分辨率建模每个场景是24维向量如果抽样1000个场景随机优化模型里的约束条件就要对1000个场景逐一展开变量规模直接膨胀三到四个数量级。对于混合整数规划问题这往往是压垮求解器的最后一根稻草——我在用某商业求解器解一个带储能的微网日前调度模型时场景数从200加到1000求解时间从4分钟涨到了三个多小时而最优目标值只变动了不到2%。这时候就需要场景削减。本质上它要做的事情是在“保留原始概率分布特征”和“降低计算规模”之间找一个平衡点。场景削减不是简单地从1000个场景里随机挑50个而是要系统地选择一组最能代表原始分布的典型场景并给每个保留场景重新分配一个合理的概率权重。1.2 削减方法选型的三个常见方向目前主流的场景削减方法大致可以归为三类我在不同项目里都试过简单说说它们的适用边界。第一类是聚类方法比如K-means、层次聚类、谱聚类。这类方法把场景当作向量做距离聚类然后取每个簇的质心作为典型场景。它的优点是速度快、易实现缺点是聚类结果只保证簇内距离最小不保证削减后场景集和原始场景集的概率分布距离最小——换句话说K-means更适合数据压缩不完全适合做概率意义上的场景代表。第二类是启发式搜索方法比如快速前向选择Fast Forward Selection和同步回代消除Simultaneous Backward Reduction。这类方法直接面向概率距离指标做优化每步选择“删除后对分布影响最小”的场景保留场景概率按邻近场景吸收合并。这一类就是在学术文献里被称为“基于概率距离的快速削减法”的核心思路。第三类是优化建模方法把场景削减写成混合整数规划或运输问题求解。理论上能得到全局最优但计算开销对于在线应用场景不太友好。综合下来对于绝大多数风光出力削减任务同步回代消除这类概率距离方法是效率和效果最均衡的选择——实现上不依赖额外工具箱只需矩阵运算和排序而且削减过程完全可复现。2. 概率距离快速削减法的原理拆解2.1 从Kantorovich距离说起“概率距离快速削减法”这个名称里的“概率距离”最关键的理论支撑就是Kantorovich距离它衡量的是两个概率分布之间的“最小搬运代价”。可以这样理解把原始场景集想象成一堆沙子每个场景是一小堆沙子削减就是要把若干小沙堆搬走、合并到剩余的沙堆上。Kantorovich距离就是完成这个搬运过程的最小工作量。严格定义下给定两个离散分布P和QKantorovich距离可以写成[ D_K(P, Q) \inf \sum_{i \in P} \sum_{j \in Q} \eta_{ij} , c(i, j) ]其中(\eta_{ij})是从原始场景(i)搬运到保留场景(j)的概率流量(c(i,j))是两个场景之间的单点代价函数一般取欧氏距离或其加权形式。削减的目标就是要找一个保留场景子集和对应的新概率分布使得这个搬运代价尽可能小。同步回代消除就是对这个目标的贪心近似每次迭代都找一个“当前删除后损失最小”的场景删掉并把它的概率按距离远近分配给剩余场景直到保留场景数量达到预设值。虽然贪心不一定保证全局最优但在实际的风光场景削减问题中效果已经非常接近理论最优而且运行速度比全局优化快几个数量级。2.2 同步回代消除的迭代逻辑同步回代消除的基本流程并不复杂我把它拆成五步第一步计算所有场景两两之间的距离得到距离矩阵(d_{ij}) 第二步对于每个场景(i)找到离它最近的场景(j^)记录距离(d_{i,j^}) 第三步遍历所有场景找到“最小的最近距离”对应的那个场景也就是删除后信息损失最小的场景把它删掉 第四步把被删除场景的概率累加到距离它最近的保留场景上 第五步重复第二到第四步直到场景数量降到目标值。这个逻辑和我第一次接触时的直觉有点反着来。我原以为应该优先删除“最不典型”的场景但算法实际上优先删除的是“最近距离最小”的场景——这其实很合理因为如果一个场景和另一个场景靠得足够近说明它们携带的信息高度冗余删除它不会造成很大的分布偏差。2.3 为什么比单纯聚类更贴合概率语义K-means之类的聚类方法有一个隐含预设所有样本的重要性是均等的。但场景生成里的每个场景是带概率权重的不同的场景可能对应不同的概率水平。比如极端天气场景出现的概率可能只有0.02但它在风险评估中却很重要。如果直接用K-means取质心大概率会把极端场景的独特信息磨平。概率距离削减法则不一样。概率在算法中始终扮演核心角色删除场景时看的是“加权后的距离损失”概率分配时按距离吸收合并。因此削减后得到的场景集其概率分布和原始分布在Kantorovich距离意义下尽可能地接近。这一点在后续评估中会体现得非常明显——保留场景的均值、方差、分位数特性都能保持得很好。3. MATLAB分步实现从时序样本到削减程序3.1 场景生成的建模细节在写削减程序之前先要把场景生成这一环做扎实。我的做法是基于风电出力和光伏出力的典型概率模型来抽样。风电出力场景通常从风速概率分布出发。风速可以假设服从两参数Weibull分布[ f(v) \frac{k}{\lambda} \left( \frac{v}{\lambda} \right)^{k-1} \exp \left[ -\left( \frac{v}{\lambda} \right)^k \right] ]其中(\lambda)是尺度参数(k)是形状参数。得到风速序列后还要用风电功率曲线把风速映射成出力。光伏出力场景我习惯用Beta分布拟合光照强度的随机性再结合光电转换效率模型计算出力。对于24小时的时间序列场景需要在每个时段分别抽样还要考虑时段之间的自相关性——直接用独立抽样会生成大量“锯齿状”的不合理场景。处理自相关性的一个实用方法是用AR(1)模型生成预测误差序列再把误差叠加到确定性预测曲线上。这样生成的风光出力场景在时间维度上有合理的连续性后续削减出来的场景才不会出现“风力出力每小时剧烈跳变”这种脱离物理实际的情况。下面是一段生成风电场景的MATLAB代码示例% 风速Weibull分布参数 lambda_w 8.5; % 尺度参数 k_w 2.2; % 形状参数 numScen 1000; % 场景数 T 24; % 24小时 % 基准风速时序某典型日趋势 AR(1)相关系数 baseWind 7 3 * sin((0:T-1) / T * 2 * pi); phi 0.85; % AR(1)自回归系数 windScen zeros(numScen, T); rng(2024); % 固定随机种子保证可复现 for s 1:numScen v zeros(1, T); for t 1:T if t 1 % 首时段直接按Weibull抽样 v(t) wblinv(rand, lambda_w, k_w); else % 后续时段用AR(1)叠加随机扰动 eta wblinv(rand, lambda_w, k_w) - lambda_w * gamma(1 1/k_w); v(t) baseWind(t) phi * (v(t-1) - baseWind(t-1)) eta; end % 限制风速上下限 v(t) max(v(t), 0); end % 通过功率曲线映射为风电出力简化线性模型 vRated 12; Prated 100; % 额定风速和额定功率 pWind zeros(1, T); pWind(v vRated) Prated / vRated * v(v vRated); pWind(v vRated) Prated; windScen(s, :) pWind; end这里要说明一点风速抽样里引入AR(1)模型的目的不是追求统计学上的完美而是要避免场景在时间维度上“失真”。有没有自相关处理生成的场景削减结果在均值上差别不大但在时序曲线形态和爬坡特性上差别很明显——后者对含储能项目的调度优化结果影响很大。3.2 同步回代消除函数的完整代码下面这段是场景削减的核心实现输入是一组场景矩阵和对应的概率向量输出是削减后的场景子集和新的概率权重function [redScen, redProb] backwardReduction(scen, prob, K) % scen: 原始场景矩阵每行一个场景每列一个时段 % prob: 每个场景的概率向量和为1 % K: 需要保留的场景数量 % redScen: 削减后的场景矩阵 % redProb: 削减后的概率向量 [N, M] size(scen); % N个场景M维 idx 1:N; % 记录原始场景索引 % 计算距离矩阵欧氏距离如需要可改为加权距离 D zeros(N, N); for i 1:N for j i1:N d sqrt(sum((scen(i,:) - scen(j,:)).^2)); D(i,j) d; D(j,i) d; end end while length(idx) K n length(idx); % 对每个场景找到最近场景的距离 minDist zeros(n, 1); for i 1:n dRow D(idx(i), idx); dRow(i) inf; % 排除自身 minDist(i) min(dRow); end % 找到最近距离最小的场景删除 [~, delPos] min(minDist); delIdx idx(delPos); % 找到与它最近且仍被保留的场景 dRow D(delIdx, idx); dRow(delPos) inf; [~, nearestPos] min(dRow); % 被删除场景的概率转移给最近保留场景 prob(idx(nearestPos)) prob(idx(nearestPos)) prob(delIdx); % 删除场景 idx(delPos) []; end redScen scen(idx, :); redProb prob(idx) / sum(prob(idx)); end这段代码不依赖任何工具箱直接用基础矩阵运算就能跑。在我的实际测试中处理1000个24维场景削减到50个场景大约需要几秒钟。不过要注意距离矩阵的计算复杂度是(O(N^2))当原始场景数超过5000时内存占用和计算时间会明显上升需要改用分块计算或并行化。3.3 场景削减程序的调用流程场景削减一般不会单独用而是作为随机优化的前置模块。我通常把完整流程组织成一个脚本% 生成风电和光伏场景假设已经生成好 windScen, pvScen % 组合为风光联合场景 numScen size(windScen, 1); jointScen [windScen, pvScen]; % 横向拼接每个场景维度48 % 等概率场景概率为1/numScen prob ones(numScen, 1) / numScen; % 调用同步回代消除 K 30; % 目标保留场景数 [redScen, redProb] backwardReduction(jointScen, prob, K); % 分离削减后的风电和光伏场景 redWind redScen(:, 1:24); redPv redScen(:, 25:48);这里有个容易踩的细节风光联合场景的变量范围差别可能很大。风电出力可能是0到100MW光伏出力可能是0到80MW如果直接拼接距离矩阵会被量纲更大的那部分变量主导。我的做法是先对每个变量维度做标准化即减去均值再除以标准差削减完成后再还原到原始量纲。4. 削减效果评估怎么样才算削减得好4.1 数值指标概率距离与统计量保持度削减得好不好不能靠眼睛看必须有量化指标。我常用的四个指标如下第一个是削减前后的Kantorovich距离这是方法的原生优化目标。计算方式是% 计算削减前后场景集合的Kantorovich距离 function dk kantorovichDistance(origScen, origProb, redScen, redProb) % 对于每个原始场景找到最近削减场景的距离按概率加权求和 N size(origScen, 1); K size(redScen, 1); dk 0; for i 1:N minD inf; for j 1:K d sqrt(sum((origScen(i,:) - redScen(j,:)).^2)); if d minD minD d; end end dk dk origProb(i) * minD; end end这个值越小说明削减后的场景集在概率分布意义上对原始场景集的代表性越好。我一般会把不同K值下的Kantorovich距离做成曲线看在哪里出现“膝盖点”——即继续增加场景数收益开始显著递减的位置。第二个指标是统计量保持度。对比削减前后场景集的加权均值曲线、方差曲线、以及特定分位数比如5%和95%计算相对误差[ \varepsilon_{\mu} \frac{| \mu_{orig} - \mu_{red} |2}{| \mu{orig} |_2} ]第三个指标是边界场景保留度。在随机优化里极端场景往往决定备用容量配置和储能容量规划。我会专门检查削减后的场景集是否还包含原始场景中出力最低和最高的那些场景。如果极致场景被磨平了这个削减方案就需要调整。第四个指标是计算收益。记录削减前后随机优化模型的求解时间比和最优值偏差。毕竟场景削减的最终目的是让优化跑得动这个直接收益指标最具有说服力。4.2 可视化对比一张图看穿削减质量我拿到一组新的场景数据后第一件事永远是画三张图。第一张图是原始场景和削减后场景的覆盖对比用半透明的颜色画出所有原始场景形成的“管子”再叠加削减后的场景曲线。如果削减后的曲线基本落在原始管子的高密度区域内说明代表性没问题如果出现明显跑到管子外面的曲线就要怀疑是标准化或距离权重出了问题。第二张图是概率分布对比。把某个特定时段比如午间光伏高发时段的原始和削减后的出力分布分别画成直方图观察削减后分布是否保持了原始分布的双峰或偏态特征。第三张图是均值与时序曲线对比把削减前后的加权均值曲线和±1倍标准差带画在一起。这是最直观的验收图我通常在论文和项目报告里直接用这一张。4.3 一个实际评估案例从1000到30用上面这套流程我做过一个典型的评估测算原始场景1000个逐步削减到500、200、100、50、30、20每一档都计算Kantorovich距离和均值相对误差。结果大致如下保留场景数Kantorovich距离均值相对误差求解时间某MIP模型100000%约3小时2000.0721.1%约35分钟500.1512.3%约8分钟300.2083.5%约4分钟200.2855.2%约2分钟从表里可以清楚看到从50个场景减到20个场景Kantorovich距离增幅明显加大均值相对误差也突破了3%——对于大多数日前调度问题3%以内的误差是可以接受的再往下砍就要谨慎评估对优化结果的影响了。5. 实操中绕不开的坑与经验5.1 场景标准化与距离函数选择第一个坑我已经提过不同变量量纲不一致。风电出力是几十到一百的量级光伏出力也是类似量级好像不大但如果联合场景里加入负荷变量负荷可能是几百上千的量级直接算欧氏距离场景削减实际上只考虑负荷差异风光出力信息完全没有参与判定。解决办法有两种。一种是做z-score标准化把所有维度统一到同一尺度另一种是给距离矩阵的每个维度赋予权重比如[ d(i,j) \sqrt{\sum_{t1}^{T} w_t \left( x_{i,t} - x_{j,t} \right)^2} ]权重可以根据对时段重要性的先验设定比如峰时段的权重高于谷时段。我在做含储能的微网项目时削峰填谷时段权重设为1.5其他时段设为1.0这样削减后的场景能更好地保留负荷和出力的峰谷联动特征。5.2 目标保留场景数怎么定目标K值没有万能答案。我见过有人不管什么数据一律减到10个场景这在很多问题上会出严重偏差。比较可靠的做法是画“K值-Kantorovich距离”曲线找到距离下降速度从陡峭转为平缓的拐点。如果这批数据本身比较平稳比如晴天为主的光伏场景K值在20左右就能取得很好的效果如果数据包含台风、寒潮这类强随机事件建议K值取到50以上否则极端场景信息很容易被磨掉。另一个参考维度是下游优化模型的预算每个场景对应一组整数变量场景数翻倍分支定界的复杂度可能不止翻倍。所以K值的选择也是在精度和计算代价之间做的权衡。5.3 概率再平衡与场景同质化问题同步回代消除有一个副作用随着迭代推进保留场景的概率会越来越不均衡。有的场景吸收了周围大量被删场景的概率权重可能涨到0.3以上有的场景则一直保持很小的权重。这本身不一定是坏事因为概率权重的差值反映了场景在原始分布中的真实重要性但如果某个场景的概率被推到特别高就容易让优化结果过度依赖这一个场景。遇到这种情况我的处理方式是给概率分布做一个正则化检查如果最大概率超过0.25就说明削减力度过大需要增加K值或者考虑在原算法基础上加入约束限制单个场景的最大概率上限。同时还要注意场景同质化问题。削减后的场景如果很多高度相似虽然Kantorovich距离不大但在随机优化里高度相关的场景会让约束条件冗余降低优化的可靠性。我会额外检查削减后的场景两两相关系数如果相关系数超过0.95的场景占比过高就需要从原始场景中强制保留一些“异质”场景。5.4 随机种子复现与并行加速做学术研究或方案比选时场景削减结果必须可复现。我在所有场景生成的抽样环节都固定随机种子否则每次跑出来的削减场景都不一样下游优化结果无法横向比较。代码里养成一个习惯在脚本最开始写上rng(2024);如果需要跑参数扫描比如不同K值、不同距离权重可以用parfor并行加速。场景削减程序的核心计算是距离矩阵和最近邻查找这部分天然可以并行。我用Parallel Computing Toolbox做过一次测试4核并行下1000个场景削减到50个时间大约能缩短到原来的三分之一到四分之一。5.5 极端场景的主动保护机制概率距离削减法是“数据驱动”的它只会按照概率加权距离做判断没有物理含义。如果原始场景集中有一个超低风速的极端场景出现概率很低它很可能在迭代早期就被删掉并被邻近场景吸收。但在电力系统的可靠性分析里这种极端场景恰恰是关键场景。针对这个问题我的做法是在削减前加一个“受保护场景集合”的约束把原始场景中出力最低和最高的各2到3个场景标记为强制保留削减算法只对未被保护的场景执行删除操作。这样既保留了概率距离削减的全局特性又确保极端情况不会被轻易抹掉。6. 概率距离削减法在项目里的典型应用延伸6.1 风光联合场景的两阶段削减对于包含风电场和光伏电站的联合系统一个常见做法是分别生成风电场景和光伏场景然后同时削减。但有一种更精细的做法是两阶段削减先对风电场景削减到K1个再对光伏场景削减到K2个然后把二者的组合场景做笛卡尔积得到K1×K2个联合场景最后再做一次全局削减。这样做的好处是能分别控制风电和光伏的精度比如在某个地区光伏的不确定性远大于风电就可以让K2 K1。我第一次用这个方法时是把某区域的风光场景分别减到20个和40个联合后再减到50个效果比一次性联合削减到50个要好不少尤其体现在光伏午间时段的分位数匹配上。6.2 多时段滚动优化里的场景在线更新在日内滚动调度里每4小时或每15分钟就要刷新一次场景。用同步回代消除做在线削减时完全没必要重新从原始蒙特卡洛场景开始跑。可以把上一轮的削减结果作为候选池结合最新的预测误差增量生成一批新场景再做一次局部削减。这样一方面保持了历史场景的连续性另一方面也大大降低了重复计算的开销。我在某微网能量管理平台上做的实测数据是完整重削减每轮需要7秒在线增量削减只需要0.8秒几乎不影响调度周期。6.3 与机会约束规划的衔接场景削减最经典的搭配产物就是机会约束规划。在风光不确定性下的储能容量配置中机会约束通过给定置信水平把概率约束转化为可用场景表示的确定性约束。削减后的场景集就是这个转化的依据每个保留场景满足约束其概率权重之和落在置信区间内整体约束就成立。这里有一点值得注意削减误差会在机会约束中被非线性放大。置信水平取95%时如果场景削减的分布误差有3%实际约束满足概率可能掉到90%以下。所以在对接机会约束规划时我通常会把置信水平往上提1到2个百分点作为裕度同时也用上一节提到的边界场景保护机制保证尾部概率不被低估。7. 一套完整的MATLAB削减脚本参考把前文提到的东西整合成一套完整可跑的脚本供需要直接上手的朋友参考%% 场景生成与削减主脚本 clc; clear; close all; rng(2024); % 固定随机种子 % 参数设置 T 24; % 调度周期小时 numScen 1000; % 原始场景数 K_target 30; % 目标保留场景数 vRated 12; % 额定风速 Prated 100; % 风电额定功率MW lambda_w 8.5; k_w 2.2; %% 1. 生成风电场景 baseWind 7 3 * sin((0:T-1) / T * 2 * pi); phi 0.85; windScen zeros(numScen, T); for s 1:numScen v zeros(1, T); for t 1:T if t 1 v(t) wblinv(rand, lambda_w, k_w); else eta wblinv(rand, lambda_w, k_w) - lambda_w * gamma(1 1/k_w); v(t) baseWind(t) phi * (v(t-1) - baseWind(t-1)) eta; end v(t) max(v(t), 0); end pWind zeros(1, T); pWind(v vRated) Prated / vRated * v(v vRated); pWind(v vRated) Prated; windScen(s, :) pWind; end % 场景标准化 mu_w mean(windScen(:)); std_w std(windScen(:)); windScenStd (windScen - mu_w) / std_w; %% 2. 执行同步回代消除 prob ones(numScen, 1) / numScen; [redScenStd, redProb] backwardReduction(windScenStd, prob, K_target); % 还原量纲 redScen redScenStd * std_w mu_w; %% 3. 评估指标 dk kantorovichDistance(windScen, prob, redScen, redProb); mu_orig mean(windScen); mu_red sum(redScen .* redProb); err_mu norm(mu_orig - mu_red) / norm(mu_orig); fprintf(Kantorovich距离: %.4f\n, dk); fprintf(均值相对误差: %.4f%%\n, err_mu * 100);脚本里调用的两个函数backwardReduction和kantorovichDistance在上一节已经给出。如果你需要把光伏场景也加进来只需额外生成光伏场景矩阵按3.3节的方式拼接即可。提示所有场景生成的随机性来自rand和wblinv函数。固定rng(2024)后每个人跑出来的结果完全一致方便复现和对比。这个算法我前后在不同项目里用了很多次最大的体会是不要迷信复杂方法也不要迷信默认参数。同步回代消除的原理看起来简单但把它用对、用稳需要对数据特征、距离定义和概率权重有实实在在的理解。你可以从这篇文章给出的代码出发先跑通一组自己的数据再根据Kantorovich距离曲线和可视化结果慢慢调整参数——这个过程本身就是对“不确定性”最直观的一课。 
长河君 · 2026/10/13 15:52:45 · 阅读 1.2k 基于概率距离快速削减法的风光场景生成与削减实战
做风光出力场景的人多半都撞过同一堵墙蒙特卡洛随便一抽就是上千个场景扔进随机优化模型里求解器直接被迭代次数和约束规模拖垮。我在某园区微网的日前调度项目里第一次遇到这个问题时光是处理一台风机和一个光伏阵列的联合场景就跑了一整晚没出结果。后来把场景数从一千砍到几十计算时间从小时级掉到分钟级优化结果却几乎没有肉眼可见的差异——关键就在场景削减这一环。这篇文章就围绕“基于概率距离快速削减法的风光场景生成与削减”展开我会讲清楚场景生成怎么建模、概率距离削减的原理为什么站得住、MATLAB代码怎么一步一步写出来以及削减效果如何量化评估。内容主要面向正在做电力系统随机优化、微电网能量管理、风光出力不确定性分析的研究生和工程师如果你只是想把场景削减模块接进自己的项目里也可以直接抄后面的代码和参数设置。1. 场景削减这件事为什么绕不开1.1 随机优化里的“场景风暴”问题风光出力的随机性是个客观事实不是靠提高预测精度就能完全消灭的。无论是风电的湍流强度、光伏的云层遮挡还是负荷侧的随机波动都需要用一组带概率的离散场景去近似连续的概率分布。蒙特卡洛抽样是最直观的做法对风速、辐照度等随机变量做大量抽样每个样本生成一条完整的时间序列就得到一个场景。问题恰恰出在“大量”上。假设风电出力用24小时分辨率建模每个场景是24维向量如果抽样1000个场景随机优化模型里的约束条件就要对1000个场景逐一展开变量规模直接膨胀三到四个数量级。对于混合整数规划问题这往往是压垮求解器的最后一根稻草——我在用某商业求解器解一个带储能的微网日前调度模型时场景数从200加到1000求解时间从4分钟涨到了三个多小时而最优目标值只变动了不到2%。这时候就需要场景削减。本质上它要做的事情是在“保留原始概率分布特征”和“降低计算规模”之间找一个平衡点。场景削减不是简单地从1000个场景里随机挑50个而是要系统地选择一组最能代表原始分布的典型场景并给每个保留场景重新分配一个合理的概率权重。1.2 削减方法选型的三个常见方向目前主流的场景削减方法大致可以归为三类我在不同项目里都试过简单说说它们的适用边界。第一类是聚类方法比如K-means、层次聚类、谱聚类。这类方法把场景当作向量做距离聚类然后取每个簇的质心作为典型场景。它的优点是速度快、易实现缺点是聚类结果只保证簇内距离最小不保证削减后场景集和原始场景集的概率分布距离最小——换句话说K-means更适合数据压缩不完全适合做概率意义上的场景代表。第二类是启发式搜索方法比如快速前向选择Fast Forward Selection和同步回代消除Simultaneous Backward Reduction。这类方法直接面向概率距离指标做优化每步选择“删除后对分布影响最小”的场景保留场景概率按邻近场景吸收合并。这一类就是在学术文献里被称为“基于概率距离的快速削减法”的核心思路。第三类是优化建模方法把场景削减写成混合整数规划或运输问题求解。理论上能得到全局最优但计算开销对于在线应用场景不太友好。综合下来对于绝大多数风光出力削减任务同步回代消除这类概率距离方法是效率和效果最均衡的选择——实现上不依赖额外工具箱只需矩阵运算和排序而且削减过程完全可复现。2. 概率距离快速削减法的原理拆解2.1 从Kantorovich距离说起“概率距离快速削减法”这个名称里的“概率距离”最关键的理论支撑就是Kantorovich距离它衡量的是两个概率分布之间的“最小搬运代价”。可以这样理解把原始场景集想象成一堆沙子每个场景是一小堆沙子削减就是要把若干小沙堆搬走、合并到剩余的沙堆上。Kantorovich距离就是完成这个搬运过程的最小工作量。严格定义下给定两个离散分布P和QKantorovich距离可以写成[ D_K(P, Q) \inf \sum_{i \in P} \sum_{j \in Q} \eta_{ij} , c(i, j) ]其中(\eta_{ij})是从原始场景(i)搬运到保留场景(j)的概率流量(c(i,j))是两个场景之间的单点代价函数一般取欧氏距离或其加权形式。削减的目标就是要找一个保留场景子集和对应的新概率分布使得这个搬运代价尽可能小。同步回代消除就是对这个目标的贪心近似每次迭代都找一个“当前删除后损失最小”的场景删掉并把它的概率按距离远近分配给剩余场景直到保留场景数量达到预设值。虽然贪心不一定保证全局最优但在实际的风光场景削减问题中效果已经非常接近理论最优而且运行速度比全局优化快几个数量级。2.2 同步回代消除的迭代逻辑同步回代消除的基本流程并不复杂我把它拆成五步第一步计算所有场景两两之间的距离得到距离矩阵(d_{ij}) 第二步对于每个场景(i)找到离它最近的场景(j^)记录距离(d_{i,j^}) 第三步遍历所有场景找到“最小的最近距离”对应的那个场景也就是删除后信息损失最小的场景把它删掉 第四步把被删除场景的概率累加到距离它最近的保留场景上 第五步重复第二到第四步直到场景数量降到目标值。这个逻辑和我第一次接触时的直觉有点反着来。我原以为应该优先删除“最不典型”的场景但算法实际上优先删除的是“最近距离最小”的场景——这其实很合理因为如果一个场景和另一个场景靠得足够近说明它们携带的信息高度冗余删除它不会造成很大的分布偏差。2.3 为什么比单纯聚类更贴合概率语义K-means之类的聚类方法有一个隐含预设所有样本的重要性是均等的。但场景生成里的每个场景是带概率权重的不同的场景可能对应不同的概率水平。比如极端天气场景出现的概率可能只有0.02但它在风险评估中却很重要。如果直接用K-means取质心大概率会把极端场景的独特信息磨平。概率距离削减法则不一样。概率在算法中始终扮演核心角色删除场景时看的是“加权后的距离损失”概率分配时按距离吸收合并。因此削减后得到的场景集其概率分布和原始分布在Kantorovich距离意义下尽可能地接近。这一点在后续评估中会体现得非常明显——保留场景的均值、方差、分位数特性都能保持得很好。3. MATLAB分步实现从时序样本到削减程序3.1 场景生成的建模细节在写削减程序之前先要把场景生成这一环做扎实。我的做法是基于风电出力和光伏出力的典型概率模型来抽样。风电出力场景通常从风速概率分布出发。风速可以假设服从两参数Weibull分布[ f(v) \frac{k}{\lambda} \left( \frac{v}{\lambda} \right)^{k-1} \exp \left[ -\left( \frac{v}{\lambda} \right)^k \right] ]其中(\lambda)是尺度参数(k)是形状参数。得到风速序列后还要用风电功率曲线把风速映射成出力。光伏出力场景我习惯用Beta分布拟合光照强度的随机性再结合光电转换效率模型计算出力。对于24小时的时间序列场景需要在每个时段分别抽样还要考虑时段之间的自相关性——直接用独立抽样会生成大量“锯齿状”的不合理场景。处理自相关性的一个实用方法是用AR(1)模型生成预测误差序列再把误差叠加到确定性预测曲线上。这样生成的风光出力场景在时间维度上有合理的连续性后续削减出来的场景才不会出现“风力出力每小时剧烈跳变”这种脱离物理实际的情况。下面是一段生成风电场景的MATLAB代码示例% 风速Weibull分布参数 lambda_w 8.5; % 尺度参数 k_w 2.2; % 形状参数 numScen 1000; % 场景数 T 24; % 24小时 % 基准风速时序某典型日趋势 AR(1)相关系数 baseWind 7 3 * sin((0:T-1) / T * 2 * pi); phi 0.85; % AR(1)自回归系数 windScen zeros(numScen, T); rng(2024); % 固定随机种子保证可复现 for s 1:numScen v zeros(1, T); for t 1:T if t 1 % 首时段直接按Weibull抽样 v(t) wblinv(rand, lambda_w, k_w); else % 后续时段用AR(1)叠加随机扰动 eta wblinv(rand, lambda_w, k_w) - lambda_w * gamma(1 1/k_w); v(t) baseWind(t) phi * (v(t-1) - baseWind(t-1)) eta; end % 限制风速上下限 v(t) max(v(t), 0); end % 通过功率曲线映射为风电出力简化线性模型 vRated 12; Prated 100; % 额定风速和额定功率 pWind zeros(1, T); pWind(v vRated) Prated / vRated * v(v vRated); pWind(v vRated) Prated; windScen(s, :) pWind; end这里要说明一点风速抽样里引入AR(1)模型的目的不是追求统计学上的完美而是要避免场景在时间维度上“失真”。有没有自相关处理生成的场景削减结果在均值上差别不大但在时序曲线形态和爬坡特性上差别很明显——后者对含储能项目的调度优化结果影响很大。3.2 同步回代消除函数的完整代码下面这段是场景削减的核心实现输入是一组场景矩阵和对应的概率向量输出是削减后的场景子集和新的概率权重function [redScen, redProb] backwardReduction(scen, prob, K) % scen: 原始场景矩阵每行一个场景每列一个时段 % prob: 每个场景的概率向量和为1 % K: 需要保留的场景数量 % redScen: 削减后的场景矩阵 % redProb: 削减后的概率向量 [N, M] size(scen); % N个场景M维 idx 1:N; % 记录原始场景索引 % 计算距离矩阵欧氏距离如需要可改为加权距离 D zeros(N, N); for i 1:N for j i1:N d sqrt(sum((scen(i,:) - scen(j,:)).^2)); D(i,j) d; D(j,i) d; end end while length(idx) K n length(idx); % 对每个场景找到最近场景的距离 minDist zeros(n, 1); for i 1:n dRow D(idx(i), idx); dRow(i) inf; % 排除自身 minDist(i) min(dRow); end % 找到最近距离最小的场景删除 [~, delPos] min(minDist); delIdx idx(delPos); % 找到与它最近且仍被保留的场景 dRow D(delIdx, idx); dRow(delPos) inf; [~, nearestPos] min(dRow); % 被删除场景的概率转移给最近保留场景 prob(idx(nearestPos)) prob(idx(nearestPos)) prob(delIdx); % 删除场景 idx(delPos) []; end redScen scen(idx, :); redProb prob(idx) / sum(prob(idx)); end这段代码不依赖任何工具箱直接用基础矩阵运算就能跑。在我的实际测试中处理1000个24维场景削减到50个场景大约需要几秒钟。不过要注意距离矩阵的计算复杂度是(O(N^2))当原始场景数超过5000时内存占用和计算时间会明显上升需要改用分块计算或并行化。3.3 场景削减程序的调用流程场景削减一般不会单独用而是作为随机优化的前置模块。我通常把完整流程组织成一个脚本% 生成风电和光伏场景假设已经生成好 windScen, pvScen % 组合为风光联合场景 numScen size(windScen, 1); jointScen [windScen, pvScen]; % 横向拼接每个场景维度48 % 等概率场景概率为1/numScen prob ones(numScen, 1) / numScen; % 调用同步回代消除 K 30; % 目标保留场景数 [redScen, redProb] backwardReduction(jointScen, prob, K); % 分离削减后的风电和光伏场景 redWind redScen(:, 1:24); redPv redScen(:, 25:48);这里有个容易踩的细节风光联合场景的变量范围差别可能很大。风电出力可能是0到100MW光伏出力可能是0到80MW如果直接拼接距离矩阵会被量纲更大的那部分变量主导。我的做法是先对每个变量维度做标准化即减去均值再除以标准差削减完成后再还原到原始量纲。4. 削减效果评估怎么样才算削减得好4.1 数值指标概率距离与统计量保持度削减得好不好不能靠眼睛看必须有量化指标。我常用的四个指标如下第一个是削减前后的Kantorovich距离这是方法的原生优化目标。计算方式是% 计算削减前后场景集合的Kantorovich距离 function dk kantorovichDistance(origScen, origProb, redScen, redProb) % 对于每个原始场景找到最近削减场景的距离按概率加权求和 N size(origScen, 1); K size(redScen, 1); dk 0; for i 1:N minD inf; for j 1:K d sqrt(sum((origScen(i,:) - redScen(j,:)).^2)); if d minD minD d; end end dk dk origProb(i) * minD; end end这个值越小说明削减后的场景集在概率分布意义上对原始场景集的代表性越好。我一般会把不同K值下的Kantorovich距离做成曲线看在哪里出现“膝盖点”——即继续增加场景数收益开始显著递减的位置。第二个指标是统计量保持度。对比削减前后场景集的加权均值曲线、方差曲线、以及特定分位数比如5%和95%计算相对误差[ \varepsilon_{\mu} \frac{| \mu_{orig} - \mu_{red} |2}{| \mu{orig} |_2} ]第三个指标是边界场景保留度。在随机优化里极端场景往往决定备用容量配置和储能容量规划。我会专门检查削减后的场景集是否还包含原始场景中出力最低和最高的那些场景。如果极致场景被磨平了这个削减方案就需要调整。第四个指标是计算收益。记录削减前后随机优化模型的求解时间比和最优值偏差。毕竟场景削减的最终目的是让优化跑得动这个直接收益指标最具有说服力。4.2 可视化对比一张图看穿削减质量我拿到一组新的场景数据后第一件事永远是画三张图。第一张图是原始场景和削减后场景的覆盖对比用半透明的颜色画出所有原始场景形成的“管子”再叠加削减后的场景曲线。如果削减后的曲线基本落在原始管子的高密度区域内说明代表性没问题如果出现明显跑到管子外面的曲线就要怀疑是标准化或距离权重出了问题。第二张图是概率分布对比。把某个特定时段比如午间光伏高发时段的原始和削减后的出力分布分别画成直方图观察削减后分布是否保持了原始分布的双峰或偏态特征。第三张图是均值与时序曲线对比把削减前后的加权均值曲线和±1倍标准差带画在一起。这是最直观的验收图我通常在论文和项目报告里直接用这一张。4.3 一个实际评估案例从1000到30用上面这套流程我做过一个典型的评估测算原始场景1000个逐步削减到500、200、100、50、30、20每一档都计算Kantorovich距离和均值相对误差。结果大致如下保留场景数Kantorovich距离均值相对误差求解时间某MIP模型100000%约3小时2000.0721.1%约35分钟500.1512.3%约8分钟300.2083.5%约4分钟200.2855.2%约2分钟从表里可以清楚看到从50个场景减到20个场景Kantorovich距离增幅明显加大均值相对误差也突破了3%——对于大多数日前调度问题3%以内的误差是可以接受的再往下砍就要谨慎评估对优化结果的影响了。5. 实操中绕不开的坑与经验5.1 场景标准化与距离函数选择第一个坑我已经提过不同变量量纲不一致。风电出力是几十到一百的量级光伏出力也是类似量级好像不大但如果联合场景里加入负荷变量负荷可能是几百上千的量级直接算欧氏距离场景削减实际上只考虑负荷差异风光出力信息完全没有参与判定。解决办法有两种。一种是做z-score标准化把所有维度统一到同一尺度另一种是给距离矩阵的每个维度赋予权重比如[ d(i,j) \sqrt{\sum_{t1}^{T} w_t \left( x_{i,t} - x_{j,t} \right)^2} ]权重可以根据对时段重要性的先验设定比如峰时段的权重高于谷时段。我在做含储能的微网项目时削峰填谷时段权重设为1.5其他时段设为1.0这样削减后的场景能更好地保留负荷和出力的峰谷联动特征。5.2 目标保留场景数怎么定目标K值没有万能答案。我见过有人不管什么数据一律减到10个场景这在很多问题上会出严重偏差。比较可靠的做法是画“K值-Kantorovich距离”曲线找到距离下降速度从陡峭转为平缓的拐点。如果这批数据本身比较平稳比如晴天为主的光伏场景K值在20左右就能取得很好的效果如果数据包含台风、寒潮这类强随机事件建议K值取到50以上否则极端场景信息很容易被磨掉。另一个参考维度是下游优化模型的预算每个场景对应一组整数变量场景数翻倍分支定界的复杂度可能不止翻倍。所以K值的选择也是在精度和计算代价之间做的权衡。5.3 概率再平衡与场景同质化问题同步回代消除有一个副作用随着迭代推进保留场景的概率会越来越不均衡。有的场景吸收了周围大量被删场景的概率权重可能涨到0.3以上有的场景则一直保持很小的权重。这本身不一定是坏事因为概率权重的差值反映了场景在原始分布中的真实重要性但如果某个场景的概率被推到特别高就容易让优化结果过度依赖这一个场景。遇到这种情况我的处理方式是给概率分布做一个正则化检查如果最大概率超过0.25就说明削减力度过大需要增加K值或者考虑在原算法基础上加入约束限制单个场景的最大概率上限。同时还要注意场景同质化问题。削减后的场景如果很多高度相似虽然Kantorovich距离不大但在随机优化里高度相关的场景会让约束条件冗余降低优化的可靠性。我会额外检查削减后的场景两两相关系数如果相关系数超过0.95的场景占比过高就需要从原始场景中强制保留一些“异质”场景。5.4 随机种子复现与并行加速做学术研究或方案比选时场景削减结果必须可复现。我在所有场景生成的抽样环节都固定随机种子否则每次跑出来的削减场景都不一样下游优化结果无法横向比较。代码里养成一个习惯在脚本最开始写上rng(2024);如果需要跑参数扫描比如不同K值、不同距离权重可以用parfor并行加速。场景削减程序的核心计算是距离矩阵和最近邻查找这部分天然可以并行。我用Parallel Computing Toolbox做过一次测试4核并行下1000个场景削减到50个时间大约能缩短到原来的三分之一到四分之一。5.5 极端场景的主动保护机制概率距离削减法是“数据驱动”的它只会按照概率加权距离做判断没有物理含义。如果原始场景集中有一个超低风速的极端场景出现概率很低它很可能在迭代早期就被删掉并被邻近场景吸收。但在电力系统的可靠性分析里这种极端场景恰恰是关键场景。针对这个问题我的做法是在削减前加一个“受保护场景集合”的约束把原始场景中出力最低和最高的各2到3个场景标记为强制保留削减算法只对未被保护的场景执行删除操作。这样既保留了概率距离削减的全局特性又确保极端情况不会被轻易抹掉。6. 概率距离削减法在项目里的典型应用延伸6.1 风光联合场景的两阶段削减对于包含风电场和光伏电站的联合系统一个常见做法是分别生成风电场景和光伏场景然后同时削减。但有一种更精细的做法是两阶段削减先对风电场景削减到K1个再对光伏场景削减到K2个然后把二者的组合场景做笛卡尔积得到K1×K2个联合场景最后再做一次全局削减。这样做的好处是能分别控制风电和光伏的精度比如在某个地区光伏的不确定性远大于风电就可以让K2 K1。我第一次用这个方法时是把某区域的风光场景分别减到20个和40个联合后再减到50个效果比一次性联合削减到50个要好不少尤其体现在光伏午间时段的分位数匹配上。6.2 多时段滚动优化里的场景在线更新在日内滚动调度里每4小时或每15分钟就要刷新一次场景。用同步回代消除做在线削减时完全没必要重新从原始蒙特卡洛场景开始跑。可以把上一轮的削减结果作为候选池结合最新的预测误差增量生成一批新场景再做一次局部削减。这样一方面保持了历史场景的连续性另一方面也大大降低了重复计算的开销。我在某微网能量管理平台上做的实测数据是完整重削减每轮需要7秒在线增量削减只需要0.8秒几乎不影响调度周期。6.3 与机会约束规划的衔接场景削减最经典的搭配产物就是机会约束规划。在风光不确定性下的储能容量配置中机会约束通过给定置信水平把概率约束转化为可用场景表示的确定性约束。削减后的场景集就是这个转化的依据每个保留场景满足约束其概率权重之和落在置信区间内整体约束就成立。这里有一点值得注意削减误差会在机会约束中被非线性放大。置信水平取95%时如果场景削减的分布误差有3%实际约束满足概率可能掉到90%以下。所以在对接机会约束规划时我通常会把置信水平往上提1到2个百分点作为裕度同时也用上一节提到的边界场景保护机制保证尾部概率不被低估。7. 一套完整的MATLAB削减脚本参考把前文提到的东西整合成一套完整可跑的脚本供需要直接上手的朋友参考%% 场景生成与削减主脚本 clc; clear; close all; rng(2024); % 固定随机种子 % 参数设置 T 24; % 调度周期小时 numScen 1000; % 原始场景数 K_target 30; % 目标保留场景数 vRated 12; % 额定风速 Prated 100; % 风电额定功率MW lambda_w 8.5; k_w 2.2; %% 1. 生成风电场景 baseWind 7 3 * sin((0:T-1) / T * 2 * pi); phi 0.85; windScen zeros(numScen, T); for s 1:numScen v zeros(1, T); for t 1:T if t 1 v(t) wblinv(rand, lambda_w, k_w); else eta wblinv(rand, lambda_w, k_w) - lambda_w * gamma(1 1/k_w); v(t) baseWind(t) phi * (v(t-1) - baseWind(t-1)) eta; end v(t) max(v(t), 0); end pWind zeros(1, T); pWind(v vRated) Prated / vRated * v(v vRated); pWind(v vRated) Prated; windScen(s, :) pWind; end % 场景标准化 mu_w mean(windScen(:)); std_w std(windScen(:)); windScenStd (windScen - mu_w) / std_w; %% 2. 执行同步回代消除 prob ones(numScen, 1) / numScen; [redScenStd, redProb] backwardReduction(windScenStd, prob, K_target); % 还原量纲 redScen redScenStd * std_w mu_w; %% 3. 评估指标 dk kantorovichDistance(windScen, prob, redScen, redProb); mu_orig mean(windScen); mu_red sum(redScen .* redProb); err_mu norm(mu_orig - mu_red) / norm(mu_orig); fprintf(Kantorovich距离: %.4f\n, dk); fprintf(均值相对误差: %.4f%%\n, err_mu * 100);脚本里调用的两个函数backwardReduction和kantorovichDistance在上一节已经给出。如果你需要把光伏场景也加进来只需额外生成光伏场景矩阵按3.3节的方式拼接即可。提示所有场景生成的随机性来自rand和wblinv函数。固定rng(2024)后每个人跑出来的结果完全一致方便复现和对比。这个算法我前后在不同项目里用了很多次最大的体会是不要迷信复杂方法也不要迷信默认参数。同步回代消除的原理看起来简单但把它用对、用稳需要对数据特征、距离定义和概率权重有实实在在的理解。你可以从这篇文章给出的代码出发先跑通一组自己的数据再根据Kantorovich距离曲线和可视化结果慢慢调整参数——这个过程本身就是对“不确定性”最直观的一课。