做新能源并网规划的时候最折磨人的不是风光资源数据本身而是这些数据怎么送进优化模型。你手上有8760个小时的风速、辐照度切成365条每日出力曲线看起来不多可真要做随机机组组合100个场景就能把内存吃光更别说蒙特卡洛抽出来的上千条曲线了。场景生成解决的是怎么产生足够多的代表场景削减解决的是怎么把这些代表压缩到优化模型吃得了的规模。这篇文章要讲的概率距离快速削减法是我这几年在风光场景处理里用得最顺手的一个方法原理简单、MATLAB代码短却能在一两分钟内把上千个场景削到十个左右分布特征还不怎么走样。1. 为什么风光优化模型需要先生成场景、再削减场景1.1 场景不是预测而是随机空间的离散近似先对齐一个基本概念。所谓风光场景指的是一条完整的时间序列比如某风电场未来24小时、每隔一小时一个出力值这条曲线代表未来可能出现的一种出力过程。它跟点预测完全不同点预测只给你一个数值期望而场景给出的是分布形状和时序相关性——光伏中午出力高、晚上为零风电可能半夜突然大发。生成场景的常见做法有三种。第一种是从历史数据直接抽样把过去三年的8760小时功率曲线按天切段得到约1000条日曲线每条曲线作为一个原始场景第二种是用概率模型模拟比如根据风机出力的概率分布和风速自相关结构用蒙特卡洛抽样生成几千条曲线第三种是把天气预报场景和历史出力做映射比如某气象预报模式下有若干种典型天气过程每种天气对应一条出力曲线。不管哪条路原始场景数量都很大因为只有足够多的曲线才能覆盖风大了但光弱、光强但风小、两者都弱、两者都强这些组合情况。在这个阶段默认所有场景等概率出现每个场景的出现权重就是1/N。这背后的逻辑是我们还没掌握更精确的未来天气概率只能用历史频率或抽样频次先均匀分配。后面削减的时候不均匀的概率分布才有意义——概率大的场景代表常见天气过程概率小的代表小概率但必须警惕的事件。1.2 场景数量失控后优化模型最先崩溃场景数如果直接喂给优化模型最先扛不住的是随机机组组合和储能容量配置这类混合整数规划问题。道理不复杂每个场景都会引入一组连续变量和对应的不等式约束而0-1整数变量的数量和场景数是相乘关系。场景数翻一倍整数变量规模也翻一倍求解时间往往不是翻一倍而是翻很多倍严重的可以直接从几分钟变成几小时甚至算不出来。我在某风光基地项目里做过一次实测一个含12台火电、96个时段的机组组合模型用30个场景还能在合理时间内出结果升到60个场景求解时间翻了将近10倍。这个计算代价没几个项目扛得住。所以场景削减不是可选项而是接入随机优化前的必经工序。这时候有三个主流思路K-means聚类、层次聚类、概率距离快速削减法。前两个思路是先找中心点再用中心点代表周围一批场景好处是速度快坏处是拿到的典型场景往往是从原始数据里算出来的平均曲线并不是真实出现过的出力过程物理上不好解释。概率距离快速削减法则不同它从原始场景里直接挑出一个子集把被删场景的概率合并到最近的保留场景上相当于合并同类项保留下来的每一条曲线都是历史上或者采样中真实出现过的场景解释性明显更好。2. 概率距离快速削减法的计算内核最小转移代价与概率重分配2.1 两个场景之间的距离怎么定义才合理要判断谁和谁长得像必须先定义场景之间的距离。最常用的是欧氏距离也就是把两个场景在所有时段上的出力差做平方和再开根号$$d(i,j)\sqrt{\sum_{t1}^{T}(x_{i,t}-x_{j,t})^2}$$这相当直观两个场景在相同时间点的出力越接近距离越小。如果场景是风、光拼接出来的联合场景比如24小时风电加24小时光伏拼成48维向量同样用这个距离公式相当于同时比较风电序列和光伏序列的整体相似度。这里有个细节如果风电和光伏的量纲不一样或者某一种资源的出力数值明显更大距离会被量纲大的那一方主导。所以实际做削减之前通常先对每一列做标准化再算距离。欧氏距离之外还有加权欧氏距离、曼哈顿距离、动态时间弯曲距离等。在风光场景这个场景下除非你特别在意早晚出力相位错位问题否则普通欧氏距离已经够用。动态时间弯曲距离虽然能解决波形相似但相位偏移的问题但计算量大对随机优化本身带来的收益又不明显性价比一般。2.2 后向消除准则谁被删掉带来的分布扭曲最小概率距离快速削减法本质上是一种贪婪的后向消除算法。它的计算目标很明确删掉某些场景后保留下来的场景集和原始场景集之间的概率分布偏差尽量小。在最优传输理论里这个偏差可以用Kantorovich距离来衡量通俗点讲就是把原始概率分布搬运到削减后概率分布的最小代价。但完整计算Kantorovich距离在工程上太繁琐所以实际算法做了一步很关键的简化每次只考虑删掉一个场景把这个场景的概率全部转移到离它最近的另一个场景上。如果某个场景的概率是p离它最近的场景距离是d_min那么删掉它造成的局部传播代价就是p乘以d_min。每次循环都找这个代价最小的场景来删相当于每一步都选择让分布走形最少的操作。伪代码可以写成这样输入原始场景集 X {x1,x2,...,xN}初始概率 p_i目标场景数 R 1. 计算距离矩阵 D其中 D(i,j) ||x_i - x_j||_2对角线设为无穷大 2. 对每个存活场景 i找到最近邻居距离 m_i min_{j ! i} D(i,j) 3. 计算删除代价 c_i p_i * m_i 4. 找到代价最小的场景 i*删除它并把概率 p_i* 加到最近邻场景上 5. 如果存活场景数 R回到第2步 6. 输出存活场景索引和更新后的概率这个算法逻辑就像一个团队裁员先看每个员工如果离职他的工作交给谁最顺、交接成本最低然后删掉交接成本最低的那个人把他的活并给最近的人。每轮都这么来直到团队规模符合预算。这种贪心策略不保证全局最优但在工程实践中已经很接近理想结果而且每一步都满足概率守恒——所有场景的概率加起来始终等于1。2.3 快速背后的工程取舍为什么标题里要强调快速两个字因为它跟传统聚类迭代完全不同。K-means需要反复更新聚类中心做多轮迭代而且聚类结果受初始中心影响很大同一个数据跑两次可能得到不同结果。概率距离削减法是确定性的只要距离矩阵算好删除顺序完全由代价排序决定不存在随机初始化的问题结果可复现对工程验证非常友好。另一个重要的点是每一步删除时不需要重新计算所有场景之间的距离。距离矩阵一次性算完后每次循环只需要更新被删场景的最近邻信息以及被删场景概率转移后接收方场景的代价。实际代码里用一个alive标志数组只对存活的场景计算代价和距离效率很高。我实测N1000、维度48、目标10个场景MATLAB跑一次大概几秒钟这个速度足够在方案比选阶段反复调整参数。3. MATLAB代码逐段实现从初始场景到削减场景3.1 用统计特征生成1000个初始风光场景为了把问题讲透我直接用代码演示完整流程。第一步是用MATLAB生成初始场景集。下面这段代码用mvnrnd函数按照风电、光伏的均值曲线和协方差结构生成1000个风光联合场景每个场景是24小时风电加24小时光伏拼接成的48维向量。clear; clc; rng(42); N0 1000; % 原始场景数量 T 24; % 每天24个时段 R 10; % 最终保留场景数 % 构造风电均值曲线夜间偏高、白天偏低标幺值0~1 windMu 0.35 0.15 * sin(2*pi*(0:T-1)/T - pi/6) 0.05 * randn(1,T); % 构造光伏均值曲线仅白天有出力 solarMu 0.60 * max(0, sin(2*pi*(0:T-1)/T - pi/2)); mu [windMu, solarMu]; % 1 x 48 均值向量 % 简单构造协方差矩阵同一种资源内时段越近相关性越强风、光之间弱相关 baseVec (1:T) - 12.5; corrT exp(-(baseVec.^2) / (2*8^2)) * 0.7 0.3 * eye(T); Sigma [corrT, 0.15*corrT; 0.15*corrT, 0.6*corrT]; % 抽样得到N0个场景并截断到[0,1] Xraw mvnrnd(mu, Sigma, N0); Xraw max(0, min(1, Xraw));如果你手里有真实历史数据不一定非要用mvnrnd。直接把历史功率曲线按天切片得到N个真实场景的矩阵Xraw同样可以喂给后面的削减函数。模拟数据的好处是先验证算法逻辑不会因为数据质量问题干扰判断跑通了再换真实数据。这里面有一个坑抽样出来的场景要记得截断。mvnrnd是按正态分布抽的天然会产生小于0或大于1的值而功率不可能为负、也不可能超过额定容量。直接截断到[0,1]会稍微改变分布尾部但相比场景削减本身引入的误差这个影响通常可以接受。3.2 削减函数完整实现与代码说明核心削减函数我写成独立文件方便复用。输入是场景矩阵X、目标场景数targetNum和可选概率向量输出是保留场景的索引和对应概率。function [keepIdx, pOut] probDistSceneReduce(X, targetNum, varargin) % 基于概率距离的快速后向场景削减 % X: N x D 矩阵每行是一个场景 % targetNum: 目标保留场景数 % 可选输入: Prob, p 指定初始概率默认均匀分布 p ones(size(X,1), 1) / size(X,1); for i 1:2:length(varargin) if strcmpi(varargin{i}, Prob) p varargin{i1}; end end N size(X,1); % 计算两两场景之间的欧氏距离矩阵 D pdist2(X, X, euclidean); D(1:N1:end) Inf; % 对角线设为无穷大避免自己匹配自己 alive true(N,1); % 存活标记 while sum(alive) targetNum aliveIdx find(alive); Dsub D(aliveIdx, aliveIdx); Dsub(1:size(Dsub,1)1:end) Inf; % 防止子矩阵对角线残留0距离 % 每个场景到最近存活场景的距离和对应列位置 [minDist, nearCol] min(Dsub, [], 2); % 删除代价 场景概率 * 最小转移距离 cost p(aliveIdx) .* minDist; [~, delPos] min(cost); iDel aliveIdx(delPos); % 被删场景的全局索引 iNear aliveIdx(nearCol(delPos)); % 最近邻场景的全局索引 % 概率转移被删场景的概率加到最近邻场景上 p(iNear) p(iNear) p(iDel); p(iDel) 0; alive(iDel) false; end keepIdx find(alive); pOut p(keepIdx); pOut pOut / sum(pOut); % 归一化保险 end这段代码有三个容易踩的细节。第一距离矩阵的对角线必须设为Inf否则每个场景都会把自己当成最近邻删除代价全部变成0削减过程直接失效。第二每次从存活索引子矩阵里取距离时要重新把对角线置为Inf因为aliveIdx不是连续索引Dsub对角线未必对应原矩阵对角线。第三概率转移之后要把被删场景的概率清零避免它在后续循环里通过aliveIdx残留影响结果。3.3 主脚本调用与结果输出约定生成初始场景之后调用削减函数就一行的事但我建议在调用之前先做标准化。下面这段主脚本完整演示了从抽样、标准化、削减到还原原始功率值的过程。%% 标准化后再削减 Xnorm zscore(Xraw); % 每列标准化避免量纲影响距离 [keepIdx, pOut] probDistSceneReduce(Xnorm, R); % 还原到原始功率场景 Xkeep Xraw(keepIdx, :); % 打印保留场景的概率 disp(保留场景索引:); disp(keepIdx); disp(对应概率:); disp(pOut); %% 可视化场景削减前后对比 figure(Color,w); subplot(2,1,1); plot(1:T, Xraw(:,1:T), Color, [0.75 0.75 0.75]); hold on; plot(1:T, Xkeep(:,1:T), LineWidth, 1.5); xlabel(时段); ylabel(风电出力(p.u.)); title(风电场景灰线为原始1000个场景彩线为保留10个场景); subplot(2,1,2); plot(1:T, Xraw(:,T1:end), Color, [0.75 0.75 0.75]); hold on; plot(1:T, Xkeep(:,T1:end), LineWidth, 1.5); xlabel(时段); ylabel(光伏出力(p.u.)); title(光伏场景灰线为原始1000个场景彩线为保留10个场景);标准化这一步很重要。如果不做风电出力如果是0到100MW光伏是0到50MW欧氏距离会自动给风电更高权重削减出的场景可能会忽略光伏分布的细节。zscore默认按列标准化实际上把所有时段的风电和光伏都压到了同一量纲两边的权重就平等了。如果你想特别强调光伏的重要性可以在标准化之后给光伏时段乘一个大于1的权重系数。输出结果里pOut每个元素对应一个保留场景的概率比如第1个场景概率0.22意味着在削减后的随机模型里系统有22%的可能性落在这一条出力曲线上。把这些场景作为输入每个场景的权重写进约束条件的目标函数里就能把随机性问题变成确定性问题求解。4. 削减效果量化概率密度、统计量与计算耗时4.1 分布层面的对比方法削减完不能只看图觉得差不多就完事还得量化评估。我的习惯是先看整体分布再看关键统计量。整体分布可以用某个代表性时刻比如正午12点光伏峰值时段、凌晨3点风电最活跃时段的出力直方图对比把原始1000个场景在该时刻的出力值画成直方图再在同一个坐标系里叠加削减后10个场景的加权直方图两个分布形状越接近说明削减效果越好。统计量方面我重点关注三个数场景集整体均值曲线、各时刻标准差、关键分位数5%、50%、95%。均值偏差反映典型的削减偏移标准差偏差反映波动性保持能力分位数偏差反映极端情况是否被保留。计算公式不复杂均值偏差可以写成$$\text{meanBias} \frac{1}{T}\sum_{t1}^{T}\left|\overline{x}{t}^{raw} - \overline{x}{t}^{red}\right|$$其中红色上标代表削减后场景的加权平均。分位数偏差则取某个时刻如正午12点比较原始分布和削减后分布的95%分位点差多少。4.2 一组典型削减结果与指标表我针对上面生成的1000个场景分别削到5、10、20、50个场景得到一组参考指标不同随机种子、不同数据结构会有波动但变化趋势一致保留场景数R均值偏差95%分位最大偏差削减耗时(秒)50.0310.280.6100.0160.150.9200.0080.091.6500.0030.053.2从表里能看出两个规律。第一削减到10个场景以后均值偏差已经控制到比较低的水平但95%分位偏差还有0.15左右这说明尾部极端场景的损失比平均场景更明显——这也是概率距离法这类基于距离合并的场景削减方法共有的弱点。第二耗时随场景数增加基本线性增长N1000规模下完全不是瓶颈。如果你在做机组组合或者储能容量配置我的经验是R取10到20之间比较合适。少于5个场景分布走形太严重极端天气可能被完全吞掉多于50个场景削减带来的计算优势就不明显了。4.3 概率分配异常时怎么处理削减完成后有时会发现某个保留场景的概率特别大比如超过了0.4甚至0.5。这不是bug而是正常现象——它说明这个场景代表的天气过程在原始数据里出现的频率本身就很高。但如果概率大到0.5以上就要注意了模型的随机性几乎被单条曲线主导本质上退化成确定性场景了。处理办法有两种。一种是提高R多保留几个场景概率自然会被摊薄。另一种是分层削减先按天气类型把原始场景分几组比如晴空光伏大发阴雨光伏低发大风阴天静稳无风在每个组内分别用概率距离法削减再把各组结果合起来。这样可以同时保证各类场景都有代表也能避免概率过度集中。我在实际项目里遇到过一次极端情况削减到5个场景时某个场景概率到了0.58导致最终的调度方案几乎等于按单场景做。后来改成分层削减把晴天、阴天、大风、小风各保留3个场景合起来12个场景概率分布就健康多了。5. 接入实际项目的高频问题与调参经验5.1 风光拼接维度与归一化的坑风光联合场景最常见的做法是把风电序列和光伏序列拼成一个长向量就像上面的例子一样48维或者96维。这里有一个容易被忽略的点距离计算是对整个长向量做的所以风电序列接近但光伏序列差异很大的两个场景联合距离可能介于中间删掉其中一个时概率会转移给整体距离最小的场景而不是风电或光伏单独最接近的场景。这通常没问题因为随机优化需要的是联合分布的近似但如果你特别在意某种极端组合比如光伏出力极低但风电偏高光靠削减可能保不住这种场景需要配合手工保留策略。另一个坑是场景矩阵的排列方向。MATLAB的pdist2是按行算距离的所以场景矩阵必须是每个场景占一行、每个时段占一列。如果你习惯把时段放第一维做出来的距离矩阵完全不对削减结果会乱成一团。这个方向问题我在第一次写代码时就踩过后面每次写完都先检查size(X,1)是不是场景数。5.2 与K-means聚类法怎么选概率距离快速削减法和K-means各有各的适用场景不是非此即彼的关系。我把两者的差异总结成一张表方便做选择对比维度K-means聚类概率距离快速削减法生成结果聚类中心通常是平均曲线原始场景子集物理可解释性一般中心曲线可能从未出现好每条保留场景都有对应历史日或采样样本概率分配按簇内场景数占比通过距离转移合并稳定性受初始中心影响确定性可复现计算速度较快N1000规模下同样很快适用场景提取典型日做分析、快速预览随机优化、需要保持分布概率、需要精确概率权重如果项目目标是给领导汇报几类典型出力曲线K-means完全够用速度快、简单但如果你要把场景直接喂给混合整数规划模型做随机优化每个场景都要配上概率权重而且结果需要用历史数据解释我更推荐概率距离法。保留下来的场景可以直接说这个场景对应去年8月3日那天的持续阴雨过程在评审会上很有说服力。5.3 极端场景的强制保留策略概率距离削减法在追求整体分布逼近的同时会把小概率的极端场景早早删掉。这些极端场景恰恰是电力系统最需要关注的风电全时段接近零出力的静稳天气、光伏连续多日阴雨、风光同时接近额定出力的特殊气象过程。这些场景概率可能只有百分之几但对系统调峰、备用容量的影响很大删掉之后优化结果可能过于乐观。我的做法是先做极端工况识别再削减。具体操作是在调用削减函数之前用规则把所有需要保护场景挑出来比如全天风电出力均值低于0.05或者光伏正午峰值低于0.1的曲线给它们打上强制保留标记然后把非强制保留的场景送进削减函数但把目标场景数改成R减去强制保留数最后把两部分场景合并概率重新归一化。如果需要更精细控制可以在削减循环里对被标记场景设置一个极大的代价让算法永远不会选中它们删除。5.4 大数据量下的运行时间优化与稳定性检查把N从1000提高到5000甚至10000时直接算完整距离矩阵会有点吃力因为距离矩阵的内存占用是N平方。我实际遇到过一个5000个96维场景的需求直接pdist2算出5000乘5000的double矩阵内存占用约200MB还能接受但到10000个场景就到800MB就会开始卡了。两个解决办法一是把距离矩阵装成single类型内存直接减半对欧氏距离的计算精度影响很小。二是分块计算距离比如每次只算500个场景与其他场景的距离用循环累积避免一次性开辟大矩阵。实现起来很简单核心就是把D pdist2(X,X)改成D zeros(N,N,single)加循环块更新。稳定性检查方面我的习惯是固定随机数种子之前先多跑几个种子。固定种子是为了代码可复现、方便团队其他人核对多跑几个种子是为了看削减结果对初始抽样的敏感程度。如果不同种子下削减出的保留场景差别很大说明原始场景集本身分布不够稳可能需要增大初始场景数N或者改用更精细的协方差模型。最后再分享一个实践小技巧做完削减后我会顺手把保留场景的序号、概率、对应原始日期一起存成表而不是只存场景矩阵。后续做调度方案复盘时能直接从表格看出系统在哪种天气场景下最紧张比对着数字翻代码高效得多。概率距离削减法不是银弹但配合上这些工程细节确实能让风光随机建模这块踏实不少。
