考虑时序相关性的蒙特卡洛场景生成与削减实践

做新能源出力场景模拟的时候我一直挺烦一种说法只要把每个时刻的分布拟合好然后蒙特卡洛随机抽就行。可真这么做生成的场景不光看着别扭拿去跑随机优化还会得到过于乐观的结果。原因很简单——风电、光伏的出力序列在时间上是有惯性的今天下午三点出力高四点大概率也不低凌晨三点低谷四点也不太可能突然拉满。单时刻独立抽样会把这层时序相关性完全打散生成的风光序列像随机噪声一样上下乱蹦跟真实出力完全不是一个物种。这篇内容就围绕“考虑时序相关性的蒙特卡洛MC场景生成与削减”展开。我会把我实际做过的流程、踩过的坑、调过的参数全部摆在台面上从怎么把时序相关性塞进蒙特卡洛抽样到场景削减时如何避免把相关性切坏再到削减前后怎么量化效果。适合正在做随机优化、电力系统规划、新能源功率预测相关研究或工程的人参考手头有历史出力序列就能直接套这套流程。1. 为什么场景生成必须处理时序相关性1.1 场景的本质是“时间序列”不是“截面数据”很多初学者容易把场景生成理解成“对每个时刻的随机变量分别抽样”然后把这些截面拼起来就完了。这确实能生成一条条曲线但这些曲线的自相关几乎为零相邻时段的值完全随机。用这种场景去做机组组合、储能充放电策略优化算出来的结果往往会偏乐观——因为系统好像拥有了“每时每刻都独立波动”的天然平滑能力。实际中的风光出力不是这样的。拿风速来说大气惯性决定了风速不能在秒级或分钟级内剧烈突变相邻时段的风速高度相关。用专业语言说就是时间序列存在明显的自相关结构。场景生成要还原的不仅是边际分布每个时刻出力的概率分布还要还原时序上的变化路径。否则生成场景的和真实序列在“动态特性”上是完全不同的两码事。我在实际项目中对比过两种场景集一种考虑时序相关性一种不考虑各自生成1000条24小时风光出力序列然后喂给同一个微电网容量配置模型。结果不考虑时序相关性时最优储能容量直接少了大概12%。听着好像“更经济”可真把方案放到历史数据里回测缺电概率马上超标。这就是典型的相关性失真导致优化结果“漂移”。1.2 忽略时序相关性会带来哪些典型问题首先是最直观的——生成序列的高频抖动。真实风电出力曲线通常是一段段较为平滑的爬坡和波动过程但独立抽样得到的序列会频繁出现“前一小时零出力、后一小时满发、再后一小时又几乎为零”这种在物理上几乎不可能的剧烈反转。其次是爬坡速率的失真。风电出力爬坡约束是实际调度中非常关注的安全指标如果生成场景不考虑时序关系爬坡速率的分布就会比真实情况偏大很多极端场景甚至会超出机组最大爬坡限制被优化模型误判为“系统无法应对的恐怖场景”从而被迫配置过量的旋转备用。还有一个隐蔽问题空间相关性和时序相关性是纠缠在一起的。比如同一区域多个风电场它们出力序列之间的互相关性本身就通过气象过程耦合。如果每个场站单独做独立抽样场站间的时空相关性就丢了这对联合优化问题的影响比单站点时序失真更严重。总之不考虑时序相关性的场景生成就像用一堆散落的照片冒充录像一段段单帧都成立连在一起就成了鬼畜视频。2. 基于蒙特卡洛的场景生成从独立抽样到引入时序结构2.1 蒙特卡洛抽样的基础流程蒙特卡洛MC场景生成的基本逻辑很简单用随机抽样的方式从概率模型中产生大量可能的未来序列用这些序列近似真实的不确定性空间。最原始的MC流程是对预测误差或出力时间序列的每个时段用历史数据拟合一个概率分布比如正态分布、Beta分布或核密度估计。在每个时段独立地从对应分布中抽样组成一条完整时间序列。重复大量次数比如几千次得到场景池。这个流程胜在简单、通用但最大的缺陷就是前面说的——忽略了时段之间的关联性。在MC框架内我们可以通过两个方向修复第一在抽样阶段引入相关性结构第二在生成后对场景矩阵做重新洗牌或校正。第二个方向通常更麻烦不如直接在抽样阶段用多元相关抽样。2.2 时序相关性的数学表达协方差矩阵或自相关函数要从数学上刻画时序相关性最简单的是用协方差矩阵。假设共有T个时段我们要生成一条长度为T的出力序列那么就把这T个时段的出力看成T维随机变量它们之间的相关性用一个 T×T 的协方差矩阵 Σ 描述。Σ 的第 i 行第 j 列元素就是第 i 个时段和第 j 个时段出力的协方差。如果序列平稳这个矩阵还是典型的Toeplitz结构对角线附近值大离得越远值越小。另一种表达是自相关函数。对长度为T的时间序列自相关系数 ρ(k) 表示间隔k个时段的两个点之间的相关性。自相关函数和协方差矩阵是等价的描述方式只是后者在多元抽样中操作起来更方便。实际建模时我们可以从历史序列中估计出这个协方差矩阵然后把它作为目标相关结构。2.3 常规做法使用Nataf变换/多元正态分布进行相关抽样在MC框架中引入相关性的最成熟方法是Nataf变换也叫多元正态变换。基本思路是先把每个时段的边缘分布转换为标准正态分布再在正态空间中构建相关模型抽样后逆变换回原始分布。具体来说如果时段 i 的边缘分布函数是 F_i(x)则令 z_i Φ^{-1}(F_i(x))其中 Φ 是标准正态分布的CDF。这样每个时段的值都变成了标准正态变量。此时需要在标准正态空间中给定一个相关矩阵 R_z它和原始空间相关矩阵 R_x 并不完全相等需要经过修正。不过如果我们直接用历史数据估计原始空间的相关矩阵再逐点逆变换通常得到的效果已经足够。实际操作中有一个更直接的近似做法直接对原始出力序列建立多元正态分布模型或者对出力归一化后的残差建立多元正态分布然后用多元正态抽样的方式生成场景。虽然出力通常不是正态分布但如果我们先做分位数变换或直接用核密度估计边缘分布再配合Nataf变换就能保留任意边缘分布同时引入线性相关性。这也是我最后推荐的做法。3. 考虑时序相关性的MC场景生成完整流程3.1 输入数据的处理风电/光伏功率序列的分布拟合整套流程的第一步是准备历史出力数据。别急着直接拿去算协方差先把数据做几个预处理剔除明显的坏数据比如通讯中断导致的零值段、停机检修段对出力做归一化通常除以装机容量变成0到1之间的标幺值如果预测场景生成中要考虑预测值另算预测误差序列而不是直接对出力抽样。做完预处理后分别估计每个时段的边缘分布。以小时级场景为例24个时段就有24个分布。对于天气过程造成的偏态分布可用核密度估计来拟合避免强拉硬套正态分布。核密度估计的好处是不提前假设分布形状数据多的时候拟合得相当好。样本量少的话用Beta分布或者混合正态分布也可以。3.2 构造时序相关结构基于历史场景的经验协方差关键的一步是构造协方差矩阵。把每一天或某一段时期的出力曲线看成一条 T 维样本把所有历史样本拼成一个 N×T 的矩阵然后直接计算列与列之间的协方差得到 T×T 的经验协方差矩阵。这里有一个经验不要直接用原始出力求协方差最好先做分位数变换到正态空间。否则因为原始出力常常在0附近堆积无风、夜间光伏为零线性协方差会低估强相关性。我习惯的做法是对每个时段的边缘分布求分位数映射把历史数据点都变换到标准正态空间在正态空间计算协方差矩阵等价于相关矩阵后续抽样也在正态空间进行最后再通过边缘分布的分位数函数逆变换回原始量纲。这套流程是Nataf变换的标准操作能明显改善低出力时段的相关性保持效果。3.3 抽样与后处理生成候选场景集有了边缘分布和正态空间协方差矩阵就可以做蒙特卡洛抽样了。用Python的numpy很容易实现import numpy as np from scipy.stats import norm # T: 时段数, n: 场景数 rng np.random.default_rng(42) # 生成标准正态空间下多元相关样本 z rng.multivariate_normal(meannp.zeros(T), covR_z, sizen) # 逆变换到原始分布F_i^{-1}(Phi(z_i)) X np.empty_like(z) for i in range(T): # 用核密度估计/拟合好的CDF的逆函数 X[:, i] marginal_inv[i](norm.cdf(z[:, i]))逆变换完成之后得到的就是n个长度为T的原始出力场景。注意这里marginal_inv[i]是每个时段的边缘分布CDF的反函数。如果没有现成逆函数用插值也能做。抽样完成后通常还有一个后处理步骤把极端出力和正常出力混合。因为纯随机抽样虽然从分布上讲是准确的但样本里可能缺少特定类型的拐点场景例如连续爬坡、持续零出力然后暴涨。我一般会额外加一些人工构造的边界场景再一起送往下游削减模块。3.4 数据保障如何验证生成场景的时序相关性生成完场景不要急着用先验证一下相关性和历史数据是否一致。常用的验证指标很直观相邻时段相关系数lag-1自相关把历史序列的lag-1自相关算出来再算生成场景集每条序列的lag-1自相关看均值或中位数是否接近半小时/小时级别的自相关函数曲线画出自相关随滞后步长的衰减曲线对比历史与生成场景的平均曲线方差-协方差矩阵的最大偏差计算生成场景的样本协方差矩阵与目标协方差矩阵的Frobenius范数误差。我用这套验证方法检查过很多次关键点在于当样本量足够大时例如5000条场景生成场景的协方差矩阵会非常接近目标矩阵。如果发现偏差特别大大概率是边缘分布变换或协方差矩阵估计的步骤出了错。此时可以先用正态数据测试整套抽样逻辑再换成真实分布。4. 场景削减在保留时序相关性前提下的场景精简4.1 场景削减为什么不能只看某时刻生成几千条场景后很多优化模型根本扛不住这么大的规模。以机组组合为例场景数量每多一倍求解时间往往指数级增长。场景削减的目的就是在一个可接受精度损失的前提下把场景数量从几千压到几十甚至几个。但削减这件事远比想象中容易“削坏”。最朴素的削减方法是为每个时刻单独去做聚类——比如对24个时刻分别跑K-means然后按各时刻最近的簇中心拼成新场景。这种做法的毛病是拼出来的场景在时序上不连贯可能从第4时刻的簇中心跳到第5时刻的完全另一个簇相关性彻底被切碎。所以必须把一条场景当成一个整体来度量距离。4.2 基于同步回代削减的常规步骤同步回代削减是目前最常用的场景削减方法之一也是一种启发式的场景挑选算法。基本思想是从候选场景集中不断去掉那些“离其他场景最近”的场景每去掉一个就把被去掉场景的概率累加到离它最近的保留场景上。重复执行直到场景数降到目标数量。具体步骤可以定义为计算所有候选场景两两之间的距离对每个场景找到最近邻距离并取出其中最小者最容易被替代的场景删掉该场景将其概率权重加到最近邻场景上更新距离矩阵重复步骤2-4直到场景数量达标。这个算法实现简单、性能稳定在电力系统中应用很广泛。需要注意距离度量要能反映时序相关性常见做法是直接使用欧式距离但欧式距离对每个时段同等加权如果时序变化剧烈可能削掉关键爬坡信息。4.3 时序距离度量动态时间弯曲DTW或加权欧式距离为了在削减过程中保留时序相关性距离度量应当对“时间偏移”有更强的鲁棒性。动态时间弯曲DTW是一个好选择它可以度量四条曲线之间形状上的相似性即便两条曲线存在相位偏移比如同一天风过程提前半小时DTW也能给出较小的距离。不过DTW计算复杂度高场景数量多时速度比较慢。我处理过千级场景DTW可以上但要注意内存和时间。另一种折衷方案是加权欧式距离——对不同时段之间的差乘以一个权重矩阵。权重矩阵可以利用协方差矩阵的逆类似于马氏距离来设计这样高方差的时段权重要低些强相关的时段差异会被适当惩罚。实际操作中我会先用普通欧式距离做初始削减再用DTW做精细挑选两阶段混合能兼顾速度和效果。5. 削减方案的改进与参数选择5.1 场景数量削减幅度的权衡削减到多少场景才合适这个问题没有统一答案完全取决于下游优化模型对场景的敏感度。我自己的经验是用于随机规划、机组组合这类问题2050个场景通常够了如果只关心期望成本1530个场景就足够如果关心极端风险如VaR、CVaR需要保留更多尾部场景建议50个以上。削减幅度过大场景集代表性下降优化结果会出现较大的近似误差。过小则计算负担重失去了削减意义。建议用“削减误差”来监测计算削减前后期望成本或典型指标的偏差如果偏差超过可接受阈值就需要增加保留场景数。5.2 保留极值场景的方法标准同步回代算法天然会把“孤立的极端场景”优先淘汰掉因为极端场景往往离谁都远是最容易被删掉的。但这会导致削减后场景集失去极端风险信息。解决办法是在削减之前先把极值场景挑出来直接强制保留。平时怎么定义极值场景我常用两个指标最大出力极值序列峰值是否接近历史最大值持续极端持续时间比如风速持续超过某阈值的时长。把这些场景挑出来单独放到保留集里其余场景再做同步回代。这个过程要小心概率权重分配被强制保留的场景要分配合理的概率不能把所有概率都压到少数场景上。5.3 削减后相关性保持效果的量化评估削减完成后同样要验证时序相关性有没有被破坏。具体做法是将削减后场景集作为带权重的抽样集合重新计算它们的样本协方差矩阵和自相关系数与削减前的场景集以及历史数据作对比。我遇到过一种情况削减到10个场景后整体协方差矩阵还勉强接近但lag-1自相关系数明显偏低。原因是保留的场景都偏“平滑”波动大的场景被删光了。后来在距离度量里增加了一项“相邻时段差”的惩罚项强制让削后场景保留爬坡波动特性。这提醒我评估指标不能只看单一统计量最少得同时看目标协方差矩阵和自相关函数两条线。6. 实际项目中的应用与坑点6.1 一个虚拟案例某微电网随机优化中的场景生成为了说清楚整套流程怎么串起来我构造了一个模拟案例。某微电网系统包含一台柴油机组、一组锂电池储能和一座风电场需要用典型日场景做容量配置优化。手头有过去两年的逐小时风电出力数据希望生成20个代表性场景。我的处理路径是数据清洗按季节分组。因为我发现冬夏季风况差异大混在一起会低估季节性相关性对每个季节分别做边缘分布拟合和协方差矩阵估计用Nataf变换抽样每个季节生成2000条候选场景用前述两阶段消除法削减到20条保留每个季节各自的典型曲线将削减后的场景作为随机优化输入得到配置方案。最终模拟结果比直接用2000条场景跑优化计算时间从40分钟降到了不到2分钟目标函数误差控制在2%以内。最关键的是把削减后的20条场景放回历史数据中做模拟时序甜度比如极端爬坡事件会不会漏掉比不削减的场景集稍差一些但误差还在可控范围内。这个案例说明场景生成和削减本身不是目的支撑下游优化决策才是。6.2 常用程序实现注意事项代码不能乱写有几个细节我吃了不少苦头协方差矩阵必须做对称正定检查。实际估计出来的矩阵有时候会因为样本不足或舍入误差出现非正定多元正态抽样会直接报错。解决方案是用特征值分解把负特征值截断到很小的正数再重构矩阵。逆变换时边界值要处理。正态空间里如果出现很大的值比如z4逆变换后得到的边缘分布分位数可能是0或1和真实出力的上下限对齐即可不要让它越过物理边界。Nataf变换的“中间层”是标准正态空间不要偷懒直接在原始空间做多元抽样。否则边缘分布被扭曲相关性也保持不当。下面是一个极简实现框架仅示意核心逻辑# 假设已有历史数据 hist_data: N x T from scipy.stats import norm, gaussian_kde # 1. 边缘分布核密度估计 求CDF逆函数 marginal_inv [] for i in range(T): kde gaussian_kde(hist_data[:, i]) # 这里从0.001到0.999逐点构建CDF插值 u np.linspace(0.001, 0.999, 2000) q kde.ppf(u) # 或者用数值积分后插值得到逆函数 marginal_inv.append(interp1d(u, q)) # 2. 变换到正态空间 z_data np.zeros_like(hist_data) for i in range(T): z_data[:, i] norm.ppf(np.clip(hist_data[:, i].map_cdf... , 1e-6, 1-1e-6)) # 用户可自行实现经验CDF R np.corrcoef(z_data.T) # 3. 抽样并逆变换 z_samples rng.multivariate_normal(meannp.zeros(T), covR, sizen) scenarios np.array([[marginal_inv[i](norm.cdf(val)) for i, val in enumerate(z)] for z in z_samples])代码省略了细节但核心思路就是“正态空间抽样逆变换”。6.3 常见坑点与解决建议坑点一只用一个季节的数据估计全年时序相关性。风光的时序相关性在季节间差异很大比如夏季光伏晴空过程与冬季完全不同。解决方法是分季节、分天气类型建立场景模型。坑点二削减时用欧式距离把“时间偏移”当成“幅度差”惩罚了。这会导致削掉延后1小时的爬坡场景。改用DTW或加权的改进距离能显著改善。坑点三概率权重分配不合理。同步回代会把被削减场景的概率叠加到保留场景上但如果直接对强制保留的极值场景也做同样操作可能导致某个权重过大。我一般对极值场景单独限制权重上限。比如每个极值场景权重不超过5%剩余概率再分配给其他场景。坑点四验证环节不足。很多人削减完直接扔进优化模型根本不知道时序相关性丢了。建议把“场景集总体统计量”的验证当成和“生成-削减”流程一样重要的一部分。宁可多花半小时检查也别让优化结果悄悄漂移。另有一个跨领域的小建议如果不是电力系统而是其他时序随机场景比如交通流量、气象要素这套“多元分布相关抽样整体削减”框架同样适用只要调整边缘分布和距离度量的语义即可。最后分享一个个人习惯我会在生成场景时额外存一份“极端场景索引”把每次抽样中峰值最大、持续最小值、连续爬坡最大那几条单独标出来。不是所有项目都需要但只要涉及安全校验这份索引能帮你快速定位最危险的时间序列不需要重新跑一堆抽样。这套东西我从一开始就吃了不少亏现在算是长教训了。