干这行时间长了你会发现但凡涉及新能源出力、负荷预测、电力市场这类问题最后都会撞到同一个坎上不确定性怎么量化。风光出力不是一条确定的曲线而是成千上万种可能的曲线。这个项目标题“考虑时序相关性MC的场景生成与削减研究”说白了就是在解决一件事——把“不确定性”变成“算得动的东西”。蒙特卡洛负责生成足够多的可能性削减负责把这些可能性压缩到优化模型能接受的规模而“时序相关性”是中间最容易翻车的那个环节。这篇文章我从头到尾拆一遍包括为什么时序相关性不能丢、MC场景生成怎么把相关性揉进去、削减为什么不能用普通聚类、以及一套我实测下来还算顺手的代码流程。适合正在做随机规划、鲁棒优化、新能源消纳分析或者储能容量配置的同行参考无论你是刚入门还是已经写了一堆场景生成的脚本里面应该都有能直接用上的东西。1. 场景生成与削减到底在解决什么问题1.1 先搞清楚“场景”是哪来的场景Scenario本质上是一个高维样本常见的是24维24小时出力曲线或者96维15分钟间隔。比如描述风电出力的不确定性不能只说“明天平均出力30%”因为优化模型需要知道的是完整的时间序列——几点到几点出力高、爬坡有多快、低谷维持多久这些都是一个场景里隐含的信息。蒙特卡洛的思路很简单根据历史数据的统计特征随机抽几千条时序曲线每条曲线在“统计意义上”都有可能出现。抽得够多这些曲线就能覆盖住真实出力的分布范围。但注意这里藏着一个关键问题如果抽样时把每个时段当成独立的生成出来的场景就是一团乱码。我见过不少人第一步就栽在这——每个时段单独采样最后画出来的曲线像心电图完全没法用。这就要引出时序相关性了。1.2 时序相关性为什么是生死线风电出力一个最典型的特征就是惯性凌晨的出力高它不会瞬时间突然跌到零因为风速变化是有过程的。光伏也好中午辐照度达到峰值它不会在10分钟内从100%掉到5%。如果你生成场景时忽略了这个特性就会出现“前一个小时满发、后一个小时零出力”这种物理上不存在的序列。从数学角度看时序相关性体现为自相关函数相邻时段、间隔几个时段的数值之间存在统计依赖。在一个正确的场景集里你应该能看到一定的“平滑惯性”——上一时段出力高的下一时段大概率也不会太低。一旦你生成了违背时序规律的场景后面所有计算都会跟着错。比如用这个场景集去算储能系统容量模型会误以为爬坡需求极其频繁把储能配置得过大算备用容量的时候又会高估系统的调峰压力。所以时序相关性不是锦上添花而是底线要求。1.3 削减是“算得动”的前提场景生成动辄几千条但随机优化模型处理不了这么多。两阶段随机规划里场景数量直接决定决策变量的规模和求解时间。1000个场景可能几分钟能解5000个场景可能几十分钟都出不来而实际枚举所有可能性又不可能。模拟某课题组的算例场景数从2000削减到20求解时间从50分钟压到了40秒削减前后的期望成本误差不到1%。削减的本质是找“最能代表全体”的少数几个场景让场景集整体概率分布尽量不畸变。这里面衡量“不畸变”的标准是什么、怎么保证时序相关性在削减后不丢我放到第3节讲。2. 生成带时序关联的MC场景三种路线实测2.1 朴素MC为什么一测就翻车最简单的蒙特卡洛做法先拟合每个时段的边际分布比如风速用Weibull分布、光伏功率用Beta分布然后每个时段独立抽样拼成一条时序曲线。我自己第一次跑这个流程的时候生成的场景大概是这样的曲线看起来“合理”边际分布也对得上绘制分位数包络线90%区间和观测数据极为吻合。结果一算自相关函数就露馅了——间隔1小时的自相关系数勉强有个0.2间隔2小时直接归零。真实数据的自相关系数间隔1小时能到0.9以上间隔6小时还有0.5。这说明独立抽样把原始数据的时序惯性彻底抹掉了。为啥边际分布拟合得再好都不行因为边际分布描述的是“单点”的概率它不管这个点跟前后点怎么串成线。时序场景的关键恰恰在“怎么串”这是时序模型和MC配合要处理的事。2.2 ARIMA MC残差重抽样最稳的入门路线我最常用的方案是对历史序列拟合ARIMAp,d,q模型然后在MC抽样时对模型残差进行重采样逐时间点迭代生成新序列。原理不复杂。ARIMA里的自回归项AR本质就是在用前几个时刻的值预测当前时刻你把残差抽得再随机AR系数也会把“惯性”重新拉回来。假设你拟合出一个AR(1)模型系数是0.85那么下一个时刻的出力总是有0.85倍上一时刻的底子在自相关自然保住了。实操中可以再简化一点不必硬套ARIMA直接估计历史序列的自相关结构用“时序Bootstrap”的思路生成。具体来说先对历史序列做一阶差分或去趋势得到平稳残差序列然后从残差序列里随机抽一个起始点按序向后拼接加上趋势项还原。我实测这个方法生成的光伏场景相邻时段自相关系数在0.92~0.98之间跟原始数据几乎一致。代码示意像这样import numpy as np def generate_scenario(residuals, trend, n_steps, start_idxNone): # residuals: 历史残差序列已平稳 # trend: 确定的趋势分量如季节均值曲线 if start_idx is None: start_idx np.random.randint(0, len(residuals) - n_steps) sample residuals[start_idx:start_idx n_steps] return trend sample这个思路之所以稳是因为它把“相关结构”打到了采样机制里而不是事后补救。2.3 COPULA方法当边际分布“不听话”的时候还有一类情况风电功率的边际分布不是标准的Weibull或者多个风电场之间还有空间相关性。这时候单条ARIMA不够用我会推荐用Nataf变换或者高斯Copula来做。思路是三步走先用历史数据把每个变量每个时段的经验分布函数求出来然后通过等概率变换把原始变量映射到标准正态空间在正态空间里估计协方差矩阵用多元正态分布做MC抽样最后再等概率映射回原始分布空间。高斯Copula本质上就是把“时序相关性”和“边际分布”解耦处理。你先分别把每个时间点的分布形态管好再用一个相关结构把时间维度串起来。这样生成的场景边际分位数和自相关结构都能同时对上而且实现不复杂核心就一个多元正态抽样常用的科学计算库都自带。有一点需要说明Copula只能捕捉“整体相关趋势”如果数据里有明显的厚尾、极端段纯高斯Copula会稍微低估极端事件频率。这时可以换t-Copula自由度设小一点比如5~8尾部的联合概率能捕捉得更紧。我实测过对风电场景来说t-Copula比高斯Copula的极端爬坡场景多生成约12%更贴近历史极端事件密度。3. 场景削减从几千条压到十条信息不丢3.1 削减的数学目标和评价标准场景削减的目标很明确从原始N个场景中选出一个子集让削减前后的概率分布距离最小。这里最关键的工具是Kantorovich距离也叫Wasserstein距离最优传输里的概念。说人话如果我有2000个原始场景选10个出来当代表那“好与不好”的评判标准就是这10个场景替代2000个场景时任意一个原始场景跟最近的代表场景之间的“加权距离”加起来要最小。权重就是这个场景的概率。实际工程里你不用手写最优传输求解直接用快速前向选择QFS的启发式逻辑就够了每轮挑一个场景加入保留集要求加入后所有未保留场景到保留集的距离总和乘以概率下降最多。这样生成的保留集虽然不一定全局最优但业内经验表明已经足够用了而且复杂度低很多。用Kantorovich距离做指标还有一个好处可以定量评估削减损失。我这边习惯在削减完成后同时输出削减前后的Kantorovich距离变化一旦距离增长曲线出现突变就说明削减过头了。3.2 三种主流削减方法对比方法核心思路优点缺点适用场景快速前向选择逐步挑最佳代表保留原始样本特征严格按概率距离逼近每轮都要算距离矩阵N大时慢时序场景、样本量在2000以内K-medoids聚类代表必须是样本点实现简单sklearn生态有现成库初始中心敏感距离度量需自定义一般场景削减层次聚类自底向上合并相近场景结果稳定可输出树状图内存占用随N²增长N过大吃力中等规模样本我强烈建议时序场景优先用快速前向选择退而求其次用K-medoids。K-means不要直接用因为K-means聚类中心的“均值”会把时序形态磨平——几个有不同爬坡时段的场景一平均出来的代表场景就变成一条平坦曲线相关性结构彻底破坏。3.3 时序约束下削减的核心距离度量要改这是我在实操中最深的体会削减算法跑不动、结果失真八成不是算法问题而是距离度量没选对。如果你把每条场景当成普通向量逐点算欧氏距离相当于把时间轴压扁成一个点时段间的先后顺序信息全被丢弃了。举个直观例子A场景是早上爬坡、晚上平稳B场景是早上平稳、晚上爬坡。逐点欧氏距离会认为两者“中等相似”因为平均下来都是全天的中度出力水平。可实际上它们的时序形态完全相反——如果这是光伏出力一个对应晴天一个对应上午阴天下午晴物理含义差了十万八千里。用欧氏距离削减时这两个场景很可能会被分到一类最终取平均生成一个“不存在”的典型场景。解决方式有两种第一种用Dynamic Time WarpingDTW距离替代欧氏距离。DTW允许时间轴伸缩对齐专门匹配形状。我测试过一个真实数据集用欧氏距离削出来的10个场景中有3条爬坡段明显被“磨平”换成DTW后爬坡率分布跟原始场景集基本重合。代价是DTW距离计算比欧氏距离慢不少N2000时一轮距离矩阵大概要算几十秒还在可接受范围。第二种加时间权重的欧氏距离。给相邻时段赋更高权重或者用一阶差分特征把场景相邻时段差值作为附加特征并到距离度量里。这个方法省CPU但效果不如DTW干净。4. 实操记录一次完整的生成—削减流程4.1 数据准备与参数设定我用某风电场一整年8760个小时的历史出力数据来做演示时间分辨率1小时T24作为单条场景的维度。注意不要跨季节混着建模风资源在冬夏两季的概率特征差异很大混在一起会出来一堆“四不像”场景。按四季拆开分别建模、生成、削减最后每季留典型日使用时按天数加权组合。关键参数按常见实践的推荐值来场景生成数N 2000削减后场景数K 10时段数T 24相似性指标DTW削减算法快速前向选择为什么要N2000采样太少了比如500极端场景根本露不出来采得太多比如10000距离矩阵的计算量呈平方增长收益却递减。2000是我测试下来性价比最高的平衡点。4.2 环境依赖与生成实现依赖库就不多说了都是科学计算和机器学习生态里常用的那几个。核心生成逻辑如下import numpy as np from scipy.spatial.distance import cdist from scipy.cluster.hierarchy import linkage, fcluster # 假设 hist_data: shape (365, 24)历史日场景 # 1. 拟合每个时段均值曲线与残差序列 mean_curve hist_data.mean(axis0) residuals hist_data - mean_curve # 2. 时序Bootstrap生成2000个场景 N, T 2000, 24 gen_scenarios [] for _ in range(N): rng_start np.random.randint(0, len(residuals) - T) gen_scenarios.append(mean_curve residuals[rng_start:rng_start T]) gen_scenarios np.array(gen_scenarios) # shape (2000, 24)这一步生成的数据你可以检查一下任意取几条画出来视觉上应该是连续的、有惯性而不是锯齿状。我是习惯再算一下自相关取生成矩阵的列间相关系数矩阵的次对角线均值应该在0.85以上才算合格。4.3 快速前向选择削减实现快速前向选择的逻辑按下面这种写法就可以跑通def forward_selection_reduction(scenarios, probs, K): n len(scenarios) selected [] remaining list(range(n)) # 预计算两两距离用DTW或带权欧氏 dist_mat compute_dtw_matrix(scenarios) # 初始选概率加权后到所有场景距离最小的场景 cost probs dist_mat first np.argmin(cost) selected.append(first) remaining.remove(first) # 每轮让“总Kantorovich距离下降量最大”的场景入列 while len(selected) K: best_gain -np.inf best_candidate None for j in remaining: # 当前未保留场景到保留集的最短距离等于场景j作为新代表时损失下降量 delta 0.0 for i in remaining: if i j: continue old_dist min(dist_mat[i, s] for s in selected) new_dist dist_mat[i, j] if new_dist old_dist: delta probs[i] * (old_dist - new_dist) if delta best_gain: best_gain delta best_candidate j selected.append(best_candidate) remaining.remove(best_candidate) return selected这段代码朴素直白N2000时在普通办公电脑上跑一轮大约2到3分钟完全能接受。想提速的话可以改成向量化写法用cdist直接算所有点到当前保留集的距离矩阵替代内层循环。K的取值也别拍脑袋定。我常用的方法是分别跑K1到K30画出Kantorovich距离随K的下降曲线。你会看到曲线在某个K值后趋于平滑这个拐点就是“性价比最优”的削减规模。拿这个项目的实际数据拐点大概在K8到K12之间所以我最后选了10。4.4 削减效果如何验证削减完不能直接拿去用按下面三件事过一遍心里才有底。第一件看代表场景的时间形态是否有物理意义。我要求10个典型场景里不能出现“在第一小时满发、第二小时归零”这种怪异序列。如果出现优先去检查距离度量和K值。第二件对比削减前后的关键统计量。这里用表格列一下这个案例的实测结果表格数据为模拟该案例得到的结果仅用于流程演示指标原始2000场景削减后10场景偏差期望出力标幺值0.3520.349-0.85%90%分位数出力0.6110.605-0.98%10%分位数出力0.1130.1162.65%相邻时段自相关系数均值0.910.88-3.3%日最大爬坡率90%分位数0.420.40-4.8%期望和分位数都在1%以内说明整体分布保住了自相关系数只掉了3个百分点说明时序惯性也基本没坏。但日最大爬坡率掉了接近5%这说明极端爬坡事件在削减里还是会被稍微低估——如果你的后续研究重点在爬坡/调峰K值建议再加大或者在削减阶段对“爬坡剧烈的场景”提高抽样权重。第三件把削减后的10个典型场景丢回随机优化模型里跑出一组决策再用全部2000个场景做回测看目标函数期望值的偏差。工程上小于3%就算合格。这个“最终用途偏差”才是最有说服力的指标因为场景削减本身不是目的目的是让优化结果不偏。5. 常见问题与排查心得5.1 削减后场景太“像”多样性消失K5时会出现这种现象削减出来的几个典型日长得都差不多形态差异很小。原因有两个一是你用的距离度量对“形状差异”不敏感——比如纯欧氏距离就会这样换成DTW立刻改善二是原始场景集本身就高度同质说明MC生成阶段没有把极端场景抽样充分需要检查生成阶段是否丢失了尾部。我踩过这个坑之后现在的习惯是削减后先目测10条曲线要求这10条里至少有2条有明显高峰、2条有明显低谷、若干条是中间过渡形态。没有这个多样性后面场景驱动出来的决策就会“过于平均”失去随机优化的意义。5.2 蒙特卡洛场景生成时忽略“跨天相关性”我的演示里场景维度是24默认假设昨天24点和今天0点之间不存在强行约束这在大多数日前调度场景中说得通。但你如果做的是多日连续优化比如储能跨日套利、中长期检修计划场景必须按48小时、72小时甚至168小时来生成。否则会人为砍掉跨日持续性特征——比如连续两三天低风速的“枯风期”事件完全丢失这会严重低估储能量需求。解决办法是生成场景时的块大小从24扩展到需要的窗口长度。应急之下也可以把两天的场景拼接前做约束——把第一天的最后6小时作为第二天的初始段强制代入这样能保住跨天衔接。5.3 距离矩阵计算太大跑不动N到5000以上时欧氏距离矩阵也要占掉约100MB内存DTW矩阵直接上GB级别。我现在的习惯是先做一次速筛——用欧氏距离做预聚类把2000个场景粗分成30到50个簇然后在每簇内部各选代表这样既省内存又保留多样性。或者直接降维用PCA把24维压到6到8维再算距离也能提5倍以上速度代价是可解释性稍微差一点。5.4 削减结果对初始随机种子敏感快速前向选择本身是确定的但生成MC场景的随机种子一变削减出的场景组合可能会有细微差别。这个不是bug而是场景生成阶段的随机性传导过来的。我的处理方式固定种子生成一次基础场景库作为项目迭代的“基准版本”。后续调参、写报告、复现结果都用同一份场景库避免每个结果都对不上。如果要做更严谨的概率评估那就多跑几个种子比如10个每个种子独立做一遍生成—削减—优化最后看优化结果的均值和方差。这比单次场景库的“精细结果”有意义得多。5.5 简化版速查表常见问题典型现象首要排查方向场景时序太毛糙自相关系数低于0.7生成阶段是否独立抽样改用残差Bootstrap或ARIMA削减后爬坡被磨平最大爬坡率明显偏低距离度量换成DTW或差分加权欧氏距离代表场景物理失真出力曲线像锯齿或阶跃检查削减的K和聚类方式禁止直接用K-means场景库不够多样形状趋同、极端场景稀少增大N检查采样是否覆盖尾部优化结果对场景集过于敏感换种子结果波动大固定基准场景库或多次采样看统计分布6. 最后的经验之谈做场景生成与研究这类项目时间长了我最强烈的感受是方法论背后的“物理直觉”比算法本身重要得多。场景生成不是靠一个复杂模型堆一堆数字就能解决的。画个图把生成场景叠在观测数据上用眼睛扫一眼能发现一多半问题。如果生成的东西画出来都不像真的那后面算得再怎么花哨结果也指望不上。再分享一个小技巧削减后的场景权重别忘了做归一化。快速前向选择输出的场景总体概率之和可能不是1丢进优化模型之前要把每个代表场景的权重设为“该场景原始概率”加上“它所代表的场景概率总和”。很多从公式到代码的细节都是在这个环节出错的。先算总概率再归一化最后再带入模型顺序一定不能乱。这套“MC生成快速前向选择削减DTW距离度量”的组合我在多个新能源出力、负荷不确定性项目里反复用过整体稳当。你要是刚开始做相关研究照着这个流程搭一套出来再根据手里的数据形态去微调应该能少走不少弯路。