计量、需求分析的交互方式我在MATLAB里实现时算了这样一笔账给定一天内的风速均值函数按小时从夜间低风速爬升到午后高风速再回落。如果硬生生从原始序列里直接取96个时刻做场景低频段占比集中的特征会被切散很难忠实地保留风速变化规律。改用概率距离消减后优先保留的是变化幅度最有代表性的时段——典型的风速回升段会被重点保留夜间平稳段则可能被合并掉一两个相近的小时。这在潮流计算里差别很大前者算出来的馈线电压波动特征更像实际运行曲线后者可能把最关键的爬坡时段磨平了。再比如并网系统的容量配置研究里如果场景削减后把少数极端低辐照、高负荷的组合丢掉最终算出的储能配置可能偏小一旦碰上连阴天就撑不住。这也是为什么同样一批数据有的论文用随机抽样法做场景削减结果储能的充放电策略会明显偏离典型场景下的最优策略就是因为随机抽样对概率密度的保留不如概率距离法精准。2.4 方法定位它不是万能的有一点需要说透概率距离快速削减法擅长的是“从已有的大量场景里挑出最有代表性的一个子集”它并不擅长“从无到有生成全新场景”。所以实际工程里通常是把两类方法串联起来用——先用ARMA、Copula或者马尔可夫链等方法生成成百上千个备选场景再用概率距离快速削减法降维得到几十个典型场景。先扩充再浓缩这符合“先广撒网后聚焦”的思路。如果数据源本身就比较单薄比如只有一两周的量测数据直接套用这套流程生成全年场景那结果多半不够看因为备选场景的总量不够削减出来的东西说服力有限。3. 代码实现拆解从关键函数到调用链路3.1 主流程设计思路我实现这套MATLAB代码时参考了一个比较经典的程序结构主脚本负责数据读取、参数设置和结果可视化核心削减部分封装成独立函数。这样做的好处是复现或调整研究场景时不需要动主逻辑只需要改输入数据和削减次数跑出来的结果立即可视化对比。整体流程基本是读取或生成原始风光场景序列数据格式采用矩阵每一行代表一个时间断面每一列代表一个备选场景设置目标削减数目比如把500个场景消减成10个调用概率距离快速削减函数内部迭代计算场景之间的概率距离、逐步剔除最不重要的场景输出保留场景的编号、对应的发生概率权重以及削减后的场景数据绘制削减前后的场景对比图、削减后的典型场景曲线和概率权重柱状图。提示很多人刚开始容易把行列方向搞反。MATLAB的惯性习惯是“列向量优先”如果原始数据里每一列是一个时间序列需要先转置成“每行一个断面、每列一个场景”的格式再喂给削减函数否则算出来的“场景距离”语义完全不对。3.2 核心削减函数的实现逻辑概率距离快速削减法的核心函数我用标准结构实现如下关键部分保留示意function [keptScenarios, keptProb, reducedData] probDistReduction(scenarios, probs, targetNum) % scenarios: 矩阵, 每一列是一个备选场景 % probs: 每个场景的初始概率 % targetNum: 期望保留的场景数 numScen size(scenarios, 2); origData scenarios; currProbs probs; removedFlag false(1, numScen); while (numScen - sum(removedFlag)) targetNum % 计算当前保留场景之间的概率距离采用2-范数 activeIdx find(~removedFlag); numActive length(activeIdx); distMat zeros(numActive); for i 1:numActive for j 1:numActive if i ~ j diffVec origData(:, activeIdx(i)) - origData(:, activeIdx(j)); distMat(i, j) norm(diffVec, 2); else distMat(i, j) Inf; end end end % 寻找使概率距离损失最小的场景进行剔除 minLoss Inf; idxToRemove -1; for i 1:numActive loss_i 0; for j 1:numActive if i ~ j loss_i loss_i currProbs(activeIdx(j)) * distMat(i, j); end end if loss_i minLoss minLoss loss_i; idxToRemove i; end end % 剔除场景并把被剔除场景的概率累加到与之距离最近的保留场景上 removedIdx activeIdx(idxToRemove); removedFlag(removedIdx) true; [~, nearestIdx] min(distMat(idxToRemove, :)); currProbs(activeIdx(nearestIdx)) currProbs(activeIdx(nearestIdx)) currProbs(removedIdx); end keptScenarios activeIdx(find(~removedFlag(activeIdx))); keptProb currProbs(keptScenarios); reducedData origData(:, keptScenarios); end这段代码里最关键的几个点概率距离的计算用的是两个场景向量之间的欧氏距离本质上是衡量“这两个场景在时序曲线形状上差多大”。如果想进一步考虑不同时刻的重要性差异可以在这个距离公式上叠加一个权重向量比如对负荷尖峰时段加大权重算出来的削减结果会更贴合调度关心的时段。剔除标准的推导对候选场景 i 计算剔除后“所有其他场景到它”的加权距离总和选总和最小的那个剔除。这个总和越小说明被剔除场景和其他场景越相似去掉它对整体概率分布的影响越小。概率再分配剔除场景后它的概率要累加到离它最近的保留场景上。这个动作相当关键它就是“概率距离削减”这个名字里“概率”二字的落点——削减之后不是简单扔掉一些曲线而是把它们的代表性折到幸存的场景里保证总概率质量守恒。注意这个函数里我用双重循环求距离矩阵在场景数量达到千级以上时会有性能压力。实际工程中可以先预分配distMat矩阵再考虑用向量化或并行计算优化。但在教学和中小规模研究中双重循环的写法更直观也更方便读者把算法逻辑“翻译”成论文里的公式。3.3 主脚本与可视化验证主脚本部分我按“生成原始场景 - 削减 - 对比展示”的顺序组织clear; clc; rng(42); % 生成模拟场景: 以典型日风速曲线为基准叠加上噪声 hours 0:24; baseWind 4 3 * sin(pi * (hours - 6) / 12); numScen 500; windData zeros(length(hours), numScen); probs ones(numScen, 1) / numScen; % 初始等概率 for k 1:numScen windData(:, k) baseWind 1.2 * randn(length(hours), 1); end % 调用削减函数, 保留10个典型场景 targetKeep 10; [keptIdx, keptProb, reducedWind] probDistReduction(windData, probs, targetKeep); % 绘制对比图 figure; subplot(1,2,1); plot(hours, windData, Color, [0.8 0.8 0.8]); hold on; plot(hours, reducedWind, LineWidth, 1.5); xlabel(时刻/h); ylabel(风速/(m/s)); title(削减前后场景对比); subplot(1,2,2); bar(keptIdx, keptProb); xlabel(保留场景编号); ylabel(概率权重); title(典型场景概率权重);这里有一个实操时很有用的细节rng(42) 固定随机种子。为什么强调这一点因为场景生成本身带有随机性如果你不固定种子每次运行主脚本得到的原始场景都不一样削减结果自然也不同。在调试阶段这会非常痛苦——你可能改了一行代码以为结果变化是代码引起的实际上只是随机数变了。固定随机种子能确保实验的可复现性这也是很多论文对算法结果做敏感性分析的基础。另外绘制削减前后的场景对比图时原始场景用浅灰色细线绘制削减后的典型场景用深色实线加粗显示两者叠加后能直观看出削减后的场景是否很好地“笼罩”在原始场景的变化范围内。如果削减后的曲线形态明显偏离原始场景包络线说明削减数目设置太少或者备选场景生成质量不足需要回头检查。4. 风景场景生成与削减的完整实操流程4.1 从原始数据到备选场景一个通用流程方法和框架说清楚了还得聊聊具体的实操闭环。我习惯把整个流程分成四大步数据准备、场景生成、场景削减、结果应用。数据准备阶段需要确定研究区域和时间粒度。时间粒度一般取1小时一天就是24个断面如果做日内滚动优化可能用到15分钟粒度一天就是96个断面。粒度越细场景向量的维度越高削减时的计算量也越大但保留的时序细节越丰富。实际上做全年8760小时的场景削减计算压力会非常大所以更常见的做法是按“典型日”来做——把全年数据按季节或月份聚类成若干类典型日再在每一类内部做场景生成和削减。每个典型日只需处理24或48个时间断面的场景计算量小很多。场景生成阶段根据可用数据情况选择方法如果只有历史风速/辐照度的统计分布参数最省事的做法是用蒙特卡洛抽样从边缘分布里抽数再通过相关性矩阵生成风速-辐照度的联合场景。如果有历史序列的时序相关性可以用ARMA模型或一阶马尔可夫链来生成能更好地保留“今天风速高明天风速也高”这种持续性特征。如果有风电场和光伏电站的地理相关性数据用Copula函数把两个边际分布“粘”起来效果更好。这一步输出的是一大批备选场景比如500到2000个数量多多益善因为削减算法本质上是在一个高维空间里做代表点选取备选集越稠密选出来的代表场景越有说服力。场景削减阶段调用我们第三部分写好的核心函数把1000个场景消减成10到20个。削减后的场景集合就构成了后续仿真优化的输入。结果应用阶段最常见的是把削减后的典型场景及其概率权重代入日前调度模型目标函数是“期望运行成本最小”——每个典型场景对应一组出力曲线和概率把各个场景下的成本按概率加权求和得到的就是考虑了不确定性的期望成本。相比只取确定性预测场景单场景的做法这种多场景期望模型明显更贴合真实运行需求。4.2 风光联合场景削减的特殊细节当题目里的“风光场景”是风速和辐照度联合场景时需要注意一个实际操作细节——量纲。风速的量纲是m/s典型值在0到15之间辐照度的量纲是W/m²典型值在0到1000之间。如果把风速和辐照度直接拼成一个向量去算欧氏距离辐照度的数值会完全主导距离计算风速的差异几乎可以忽略。这样削减出来的场景很可能风速特征保留得不好而辐照度特征保留得“过好”达不到两种能源特性均衡保留的目标。解决办法是在拼装场景向量之前做归一化或标准化处理。我常用的做法是 min-max 归一化把风速和辐照度分别压到 [0,1] 区间再拼接等效于给两种变量相同的初始权重。如果研究更关注某一种资源也可以故意把该变量的缩放比例调大实现加权效果。% 风光联合场景拼装前的归一化示例 windNorm (windData - min(windData)) ./ (max(windData) - min(windData)); irrNorm (irrData - min(irrData)) ./ (max(irrData) - min(irrData)); jointData [windNorm; irrNorm];还有一点值得注意风速和辐照度在很多地区存在天然的负相关——白天辐照度高的时候风速可能偏低或偏高需要用Copula函数保留这种相关结构。如果直接在拼装时独立抽样就会把相关性抹掉生成的联合场景在物理上就不自洽。4.3 削减数目的选择技巧目标削减数目不是拍脑袋定的。它直接影响后续优化问题的规模和精度。从优化角度看保留的场景数越多期望成本模型的精度越高但决策变量的规模和求解时间也越大。对典型日前调度模型来说如果保留场景数从5个增加到20个模型里的变量数量会成倍增加混合整数规划求解时间可能从几十秒膨胀到几十分钟。所以实际中需要做一个权衡先用10个场景跑一遍再逐步增加到20、30个观察目标函数值的变化。如果目标函数值基本稳定比如变化幅度小于1%就没必要保留更多场景了。从统计角度看有一种经验法则是保留场景数大致取备选场景数的平方根。比如从400个备选场景里削减保留20个左右比较合适。这个经验法则来自聚类分析中的类数选择直觉不一定严格最优但作为初始值非常好用。之后再结合目标函数值做敏感性分析找到一个“性价比最高”的削减数目。我在实操中常用这样一组对比分别保留5/10/20/30个场景每个配置跑同一组优化问题记录求解时间和目标函数值做成一张二维对比表。目标函数从5到10个场景通常变化明显从10到20个场景变化收窄从20到30个场景几乎不动这时候就锁定10到20个作为最终配置。5. 常见问题与排查技巧5.1 削减结果出现负值或超界我自己上手时遇到的第一个问题就是削减后的场景曲线出现负值或超过历史最大值的点。排查下来根源在场景生成阶段如果直接用蒙特卡洛从正态分布抽样抽到的风速值可能是负的或者光伏辐照度超过1000W/m²。这在物理上不合理。解决办法有两个层次。第一层是生成阶段就限制抽样后把小于零的值截断为0把大于辐照度上限的值截断为上限。这虽然简单粗暴但对后续削减影响不大因为截断后的分布依然保留了主要的概率质量。第二层是削减后处理如果削减出来的典型场景里还有非物理值可以再做一次边界钳制。需要注意的是钳制操作会轻微改变场景间的距离结构但对于工程研究来说影响完全可以接受。5.2 削减后概率权重总和不为1这是概率再分配实现错误时的典型症状。正常情况下初始概率和为1每次剔除场景时把它的概率累加到最近场景上总和应该恒为1。如果发现削减后概率总和偏离1大概率是代码里“最近邻居”的查找逻辑写错了比如找最近邻居时把被剔除场景自己也包含了进去导致概率累加时加到了0上。调试手段很简单在每次剔除场景后打印所有保留场景的概率总和观察是否恒为1。如果某一轮开始偏离说明问题出现在那一轮的最近邻查找或概率更新逻辑上。这类问题用断点调试比空想快得多。5.3 削减结果不够“典型”这个问题的表现是削减后的场景曲线五花八门很难看出共性或者每个场景都像是加了很多噪声的毛刺曲线而不是光滑的典型日曲线。出现这种情况多半是备选场景生成时的噪声太大削减算法被高维噪声带偏了。解决方案是降低噪声方差或者在生成备选场景之前先做一次平滑滤波。实际操作中我建议把噪声强度的选择与后续研究需求绑定如果是研究极端场景对系统的影响噪声可以适当大一些保留尾部风险如果是研究典型日运行方式噪声应该小一些让削减后的场景更平滑、更接近理想典型曲线。5.4 性能瓶颈与仿真卡顿当备选场景数量超过2000、时间断面超过48时双重循环计算距离矩阵的耗时可能达到几分钟量级这时候主脚本调用削减函数会明显卡顿。几个优化思路先按时间维度降采样比如从15分钟粒度改成30分钟或1小时维度降一半计算量减半在MATLAB里预先分配distMat避免循环内动态扩容使用向量化计算距离矩阵比如利用矩阵乘法和广播特性替代双重循环。% 向量化计算欧氏距离矩阵的示例 numActive size(activeData, 2); a2 sum(activeData.^2, 1); distMat2 sqrt(max(bsxfun(plus, a2, a2) - 2 * activeData * activeData, 0));这段代码里的核心技巧是任意两列向量之间的欧氏距离平方等于两列向量自身模长的和减去两倍的点积。使用bsxfun或隐式扩展可以一次性得到整个距离矩阵避免显式循环。我实测下来500×24的数据量向量化版本比双重循环快上几十倍场景规模越大越明显。5.5 削减结果与预期场景差异大最后一种常见情况是削减出来的场景时间分布与直观经验不符。比如冬季典型日的风速场景削减结果居然没有保留明显的冬季大风过程反而凑了一堆过渡季节曲线。这种时候首先要怀疑的是备选场景生成的季节性偏差而不是削减算法本身。检查生成备选场景时使用的分布参数是不是按季节分别统计的如果全年混在一起生成季节性规律就会被平均掉削减自然也找不回冬季特征。解决办法是按季节或月份分别生成备选场景再分别削减最后按各季节的天数占比加权组合成全年典型场景集。6. 一个小技巧与个人体会最后分享一个我在多次实操中总结出来的小技巧不要只用削减后的场景做一次优化就收工一定要追加一个回代校验。回代校验的做法是把削减后的每个典型场景当作确定性输入算出一组调度决策然后把所有备选场景包括被削减掉的那些代入这批决策计算它们各自的运行成本和约束越限情况最后按原始概率加权统计。如果加权越限率接近零、期望成本与基于削减场景计算的期望成本相差不大说明削减后的场景集很好地代表了原始场景集如果差异明显说明削减数目太少或者削减算法丢失了重要尾部风险需要调整。这个方法从理论上对应了“场景削减误差评估”的思路——削减本身只是手段核心目的是让削减后的场景集合在后续优化中产生和原始场景集足够接近的决策质量。多花五分钟做一次回代校验比反复调算法参数要靠谱得多。我在用这套流程处理一个模拟项目X的风光容量配置研究时最开始削减到10个场景回代校验发现极端高风速场景下的出力越限率达到4%后来把保留场景数提高到20个越限率降到1%以内同时求解时间增加了大约一半。这个代价换取的风险控制收益我认为非常值得。概率距离快速削减法作为一个场景降维工具真正吃透之后你会发现它不只是论文里的一个公式而是可以反复使用、稳定输出的一把好用的尺子。希望这篇文章能帮你在自己的研究或项目中把风光场景生成与削减这件事做得又快又好。