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

先解释一个容易让人走错片场的词MC。在游戏圈这个缩写大概率指向Minecraft但在电力系统、能源优化和随机规划这个圈子里MC只有一个含义——Monte Carlo。这篇文章要聊的是蒙特卡洛场景生成与场景削减中的时序相关性问题。标题里的“MC场景生成与削减”说白了就是用蒙特卡洛采样生成大量可能的风光出力时序曲线再用削减算法挑出少数几条有代表性的曲线作为随机优化、概率潮流或储能规划的前置输入。为什么不能直接从历史数据里随便挑几天来用为什么生成的场景必须强调时序相关性削减又是怎么做到“砍掉95%的数据但基本不损失信息”的这些问题我会结合自己的实际项目一个个拆开讲。正在做随机调度、机组组合、源网荷储协同优化、储能容量规划或者刚接手不确定性分析项目的同学应该能从这篇文章里省掉不少弯路。1. 为什么随机规划需要“先铺路再修路”场景生成与削减到底解决什么问题1.1 风光出力本质上是一个随机过程不是一个随机数序列很多新手入坑随机规划时第一反应是把每个时刻的不确定性当成独立的随机变量然后对每个时段采样取值拼成一条“场景”。这种做法看着简单实际上从根上就错了。风电和光伏出力不是一串互不相关的随机数而是有明确时间记忆的随机过程今天的风速大概率影响明天的出力连续几天晴天后光伏午间出力会更高风机出力在几个小时内持续爬坡或者持续跌落才是常态一小时内从满载掉到零这种毛刺在实际物理世界里几乎不会发生。所谓“时序相关性”就是这种时刻与时刻之间的依赖关系。它不能用某时刻单独的均值、方差来描述。你可以把一条风电出力曲线理解成一段连续的故事前后情节是有因果的。如果在建模时忽略这种因果生成出来的场景集合就会充满不可能出现的剧烈跳变而优化模型恰恰会利用这些假跳变来“钻空子”最后算出来的机组组合、储能容量自然就失真了。所以在做场景生成之前必须先正视一个前提我们要生成的不是一堆独立样本而是一组时序样本。这决定了后面所有步骤的选型和实现方式。1.2 为什么不能拿全部历史数据直接算场景树与组合爆炸有人可能会问既然历史数据是最真实的场景为什么不直接把过去几年的8760小时曲线全部丢进优化模型里非要先生成、再削减这不是多此一举吗答案很简单组合爆炸。随机规划里每个场景对应一条完整的时间路径目标函数要对所有场景取期望。如果直接用全部小时级历史数据或者每个时段离散取10个可能值30个时段就会有10的30次方种分支组合。这个规模扔给任何求解器都是灾难。更实际一点即使你只把每一年的历史曲线当作一个场景几十年样本也远远不够覆盖可能的不确定性空间。所以行业里的标准做法是两步走第一步用蒙特卡洛或者其他采样方法生成大量能覆盖不确定性的场景第二步用场景削减算法把成千上万条曲线缩编到几十条甚至十几条。削减后的场景集合会作为场景树的分支进入后续的随机优化模型。这样做既保留了不确定性信息的丰富度又把求解复杂度控制在了可接受范围内。1.3 一个小算例感受削减的价值我自己做过一个某风电场额定容量100MW的随机机组组合项目。用蒙特卡洛生成了10000条48小时出力场景直接丢进MILP模型求解器跑了两小时还没有收敛内存已经逼近极限。后来用同步回代削减把场景压到20条求解时间缩短到15秒左右目标成本与几百条场景下的结果偏差只有不到0.4%。场景数量求解时间目标成本万元与全场景最优解偏差10000超过2小时未收敛——500约23分钟113.2—100约4分钟113.00.18%20约15秒112.80.35%这个算例透露了一个关键信息只要能保证削减后的概率测度逼近原始测度几十条场景的优化结果完全可以逼近上千条场景的结果。但这句话有一个大前提——削减算法必须照顾到时序结构。如果距离度量选错时序相关性被破坏哪怕只削减到500条结果都可能比正确削减到20条差得多。这也是本篇要反复强调的核心。2. 先把时序相关性讲透为什么独立采样生成的场景是“假场景”2.1 时序相关性与空间相关性是两回事做多风电场课题的人对“相关性”这个词往往先想到空间相关性同一时刻不同风电场出力之间的相关结构。比如台风过境时一群风电场会同时满发这就是空间相关。但标题里说的“考虑时序相关性MC”指的是另一件事同一个风电场不同时刻之间的自相关结构。这两种相关性很容易被混淆。有的论文声称“考虑了相关性”仔细一看只是生成了带有交叉相关矩阵的样本但每个电场的时间序列依然是一堆独立毛刺。这对单场随机优化影响还不大一旦进入需要刻画连续爬坡、持续低出力或持续高出力过程的问题问题就暴露了。描述时序相关性的工具通常是自相关函数ACF、偏自相关函数PACF和爬坡率分布。ACF可以告诉你滞后1小时、2小时、24小时的出力之间到底有多大关联爬坡率分布则反映了变化过程的统计特征。把这些指标画出来之后你会很清楚自己的场景集是否丢失了原始数据的时序节奏。2.2 丢掉自相关结构之后优化器会怎么“钻空子”丢掉时序相关性的后果不是“场景看起来不够平滑”这么简单。它在不同优化问题里的破坏方式完全不同。在机组组合问题里独立采样会产生相邻时段出力剧烈跳变的假场景。优化器一看系统调峰需求极大但每条场景的爬坡约束都被这些假跳变推得很紧最后要么超配大量备用容量要么因为平均化效应低估真实爬坡风险两种结果都会导致决策失真。在储能容量规划问题里独立采样会把真实的持续高风速过程打散成随机起伏的片段。原本应该是一段连续48小时的大风过程被拆成了一个个孤立的“高出力点”储能“低充高放”的套利空间被严重误判最后算出来的容量配置既可能偏大也可能偏小全看随机种子给不给面子。在输电通道规划问题里持续多日的区域性强出力过程才是决定通道容量的关键约束场景。独立采样把这种长时序事件稀释成概率极低的“巧合”通道容量规划自然偏向乐观。用个直白的类比独立采样就像把一部电影的每一帧随机打乱单帧画面都是正常的但连起来看完全不是一个故事。风光出力过程也一样时序上下文本身就是信息丢掉上下文等于把事故隐患埋进了优化结果。2.3 描述时序相关性的常用数学工具在工具选择上不同场景有不同做法我用过一个表格整理常见工具的适用边界工具描述对象优点局限ACF/PACF自相关结构直观、易画图只描述线性相关lag相关矩阵多时刻联合相关可直接用于蒙特卡洛采样需要保证矩阵正定ARIMA/AR时序动态方程能生成任意长度样本线性假设较强Copula边际分布相依结构灵活处理非正态边际拟合与采样复杂度高块自举对原始序列重采样保留全部统计特征样本多样性受历史限制我的经验是对于单个风电场、短期出力场景AR(1)过程加上经验边际分布通常就够用如果涉及多个风电场、多个能源品种联合出力或者需要严格刻画极端事件建议转向秩相关矩阵或Copula路线。不必一上来就上复杂的Copula先看数据和问题的需求够用就好。3. Monte Carlo场景生成实操从独立正态采样到带时序相关的样本3.1 最朴素的生成流程以及它错在哪先看一段很多人会写的“基础MC”代码。思路很简单拟合一个边际分布比如风电功率的Weibull分布或经验分布对每个时段独立采样再用逆变换得到功率值。import numpy as np from scipy.stats import norm T 48 # 时段数 N 10000 # 场景数 # 假设历史出力已经拟合了经验分布逆函数 F_inv # 错误示范逐时段独立采样 U np.random.rand(T, N) X_wrong F_inv(U)这段代码只保证了边际分布正确但没有引入任何时序依赖。此时每一条场景曲线都是前后时段完全独立的毛刺ACF在滞后1阶的位置就掉到接近0与实际出力过程的自相关特征完全对不上。问题根源在于U矩阵的各行之间是独立的。要生成带时序相关的样本第一步是让底层随机数带上相关结构再做边际变换。3.2 用Cholesky分解注入相关性的标准做法Cholesky分解是生成相关高斯样本最常用的工具核心公式很简单如果U是独立标准正态向量相关矩阵R可以分解为RL·L^T那么YL·U的协方差矩阵就是R。也就是说只要我们能构造出描述时序相关性的相关矩阵R就能通过线性变换把独立的随机数变成带相关性的随机数。具体操作分几步走根据历史数据的ACF估计滞后相关矩阵R。常用做法是假设指数衰减结构例如R[s,t]ρ^|s-t|ρ是滞后1阶自相关系数对R做Cholesky分解得到下三角矩阵L生成T×N的独立标准正态矩阵U计算YL·U得到带相关性的正态样本对Y的每个元素求标准正态CDF得到(0,1)均匀分布的样本再用F_inv把均匀样本变换到目标边际分布。代码如下from scipy.stats import norm rho 0.9 T 48 N 10000 # 构造指数衰减相关矩阵 R np.array([[rho ** abs(i - j) for j in range(T)] for i in range(T)]) L np.linalg.cholesky(R) # 生成相关正态样本 U np.random.randn(T, N) Y L U # 逆变换到均匀分布再到目标边际分布 V norm.cdf(Y) X F_inv(V)这套流程在很多开源库里都有封装但理解底层逻辑还是很重要因为后面所有排查问题都要回到这五步上。3.3 非正态边际下的秩相关校正Iman-Conover法上面这个方法有一个潜在陷阱风光出力明显不是正态分布所以我们加了逆变换这一步。但非线性变换会改变Pearson相关的大小直接导致最终生成序列的相关矩阵和设定值不吻合。比如你设了ρ0.9经过Weibull或经验分布逆变换后实际Pearson相关可能掉到0.7甚至更低。这不是Cholesky失效了而是因为非线性单调变换只保持Spearman秩相关不变不保持Pearson相关。所以对于非正态边际更稳妥的做法是用Spearman秩相关矩阵来定义时序相关然后用Iman-Conover法做采样。Iman-Conover法的核心思路先生成一个秩相关结构正确的正态样本Z然后按照目标边际分布样本的排序方式替换Z中每个位置的数值。具体来说用Cholesky生成秩相关结构正确的正态样本Z对Z的每一列做排序记录每个位置的秩次从目标边际分布中抽取样本或者从历史数据中重采样也做排序把目标分布排序后的值按Z中对应位置的秩次填回去。这样得到的新矩阵边际分布是目标分布秩相关结构也保持住了。这种方法在处理风电功率、光伏功率这类强偏态分布时非常实用。如果在实际项目中只是直接用Cholesky逆变换而不做秩相关校正削减前的场景就已经失真了后面再怎么削减都是给别人擦屁股。3.4 更省事的替代历史数据块自举如果不想陷在分布拟合和相关矩阵构造这些细节里还有一个几乎不会出错的替代方案块自举。做法是把历史出力序列按固定长度切成块然后随机抽取这些块并拼接成新场景。因为块内就是真实历史数据时序相关性、爬坡特征、日周期变化都会被原封不动地保留下来。块自举的优点是实现简单、不需要建模、对非平稳数据也基本有效缺点是生成场景的多样性受历史样本量限制如果历史数据只有一年翻来覆去就是那365天的片段重组极端场景可能被严重低估。我自己的用法是用参数化方法生成大批量场景做主力用块自举做交叉验证两边对不上时再回头检查模型假设。4. 场景削减10000条曲线砍到20条怎么砍才不心疼4.1 削减的本质是概率测度逼近不是挑几条“代表曲线”很多新手理解场景削减时以为就是选几条看起来不一样的曲线留着、其他删掉。这个理解太浅了。场景削减本质上是一个概率测度逼近问题原始MC场景集可以看成一个离散经验分布每个场景权重相等削减的目标是找到一个小支撑集的经验分布使它与原始分布之间的某种距离最小同时给保留的每个场景分配新的概率权重。这个角度很重要因为它决定了削减算法的每一步都在做什么。同步回代削减里被删场景的权重要加到最近保留场景上就是在做“概率质量转移”快速前向选择里每个代表场景的权重等于它管辖的所有原始场景权重之和也是在重建一个新的离散概率测度。如果只是“挑几条像样的”权重分配这步就会被忽略优化模型里的期望目标函数就会算错。4.2 同步回代缩减最经典的路线同步回代是实际项目里最常被优先尝试的算法思路非常直观每次都把“删掉后对整体概率测度损害最小”的那个场景删掉把它的权重转移到离它最近的场景上一直重复到剩余场景数满足要求。教学版的简化实现长这样def backward_reduction(X, K): 同步回代缩减简版 N len(X) w np.ones(N) / N D pairwise_distances(X) # 需要提前实现距离矩阵 idx list(range(N)) while len(idx) K: min_cost np.inf to_remove None nearest None for i in idx: neighbors [j for j in idx if j ! i] j_i min(neighbors, keylambda j: D[i, j]) cost w[i] * D[i, j_i] if cost min_cost: min_cost cost to_remove i nearest j_i # 被删场景的权重转移到最近场景 w[nearest] w[to_remove] idx.remove(to_remove) return X[idx], w[idx] / w[idx].sum()这个简化版适合教学和中小规模场景真正处理上万条场景时复杂度偏高。实际项目里我会用Heitsch-Römisch等人的快速实现或者直接调用现成库但脑子里保留这个简化流程有助于理解每一步在干什么调参数时不容易出错。4.3 快速前向选择大数据量下的务实选择与同步回代的“删”不同快速前向选择是“选”不断从候选场景中挑一个加入代表集使代表集的覆盖能力提升最大。它的特点是即使场景数量达到几万条性能依然可以接受因此在大规模MC采样后的削减环节非常实用。权重分配逻辑也清楚每个代表场景的最终权重等于所有离它最近的原始场景权重之和。这个“比赛划地盘”的过程本质上就是生成一个离散测度。两种方法在实操中的选择原则我整理了这张表方法方向适合规模代表场景性质典型限制同步回代删除数千以内保留真实场景距离矩阵计算O(N^2)快速前向选择数万也可保留真实场景迭代次数多时耗时上升K-means聚类海量产生平均场景可能不物理4.4 聚类路线K-medoids通常比K-means更对味还有一大类削减思路是聚类。K-means计算效率高但它有个致命问题中心点是簇内样本的均值在时序场景里均值会产生“既不高也不低”的假曲线。风电场景尤其不能接受这个因为平均曲线把爬坡过程平滑掉了接入机组组合模型后原本约束紧张的爬坡时段可能会被抹平结果偏向乐观。K-medoids聚类在这点上明显更适合场景削减它的代表点必须是真实存在的场景不会凭空造出一条物理上不可能出现的曲线。配合DTW距离使用时还能保留时间序列的形态特征。如果项目时间紧也可以在原始序列基础上拼上一阶差分序列构成增强特征向量再用K-medoids聚类既保留部分爬坡信息又控制计算量。5. 削减后时序相关性还在吗距离度量与检验不能省5.1 欧氏距离为什么在时序场景上翻车场景削减算法里距离度量决定了哪些场景会被合并、哪些场景会被删除。很多人默认选欧氏距离但这对时序场景来说是个很容易埋雷的坑。欧氏距离逐点比较对时间轴漂移极其敏感。两条曲线形状几乎一样只是整体平移了一个小时欧氏距离会很大但它们在物理上代表的是非常相似的过程。反过来两条曲线每个时刻数值接近但一条在爬坡、一条在下坡欧氏距离可能很小却被削减算法判断为“相似”而合并。这样削减出的代表场景ACF、爬坡率分布都会和原始集合明显偏离。所以在做时序场景削减时距离度量本身就是“时序相关性”的载体选错了度量后面的所有努力都会被抵消。5.2 用DTW增强时序感知想要让削减“看懂”时序形态有两个方向可以走。第一个方向是用动态时间规整DTW距离替代欧氏距离。DTW允许在时间轴上做非线性对齐衡量的是两条曲线的整体形态相似度对相位漂移、伸缩变形都比较鲁棒。代价是计算量比欧氏距离大很多上万条场景的距离矩阵可能要算到怀疑人生。我的做法是先用增强特征向量做一次预削减缩小候选集再用DTW做最终削减。第二个方向是在特征层面动手脚。把原始曲线和它的一阶差分爬坡率序列拼接起来构成一个特征向量再做距离计算。这样保留了相对丰富的时序动态信息又能继续使用欧氏距离计算效率高很多适合项目第一版快速出结果。5.3 削减结果的三层检验削减完不是看一眼曲线形状就完事了。我自己的习惯是做三层检验缺一不可。第一层是概率分布层对比削减前后出力均值、标准差、分位数和极值确认边际分布没有明显漂移。第二层是时序结构层对比ACF、PACF、爬坡率分布和持续出力时间确认时序相关性保住了。第三层是决策层把削减前后场景集分别放入同一个优化模型比较目标值和方案解的差异这一层最直观也最能说明问题。一个示例检验结果长这样指标原始10000场景削减后20场景出力均值p.u.0.4120.409标准差p.u.0.2180.221滞后1阶自相关0.870.85最大小时爬坡p.u.0.630.61优化目标值万元113.2112.8如果决策层的目标值偏差能控制在1%以内说明削减损失基本可接受如果偏差明显更大不要着急增加场景数先回头检查距离度量是不是选错了往往问题出在这里。6. 实操中踩过的坑非正定矩阵、零概率场景与极端场景流失6.1 Cholesky分解报错协方差矩阵非正定带时序相关性的MC几乎每个人都会撞上同一个报错np.linalg.cholesky抛出LinAlgError: Matrix is not positive definite。我第一次遇到时查了半天后来发现这根本不是罕见问题而是高维lag相关矩阵的常态。排查链路一般是这样的先看R矩阵的特征值如果有负特征值或者接近0的特征值说明矩阵不满秩。常见原因有三个第一时段数比样本数还多相关矩阵估出来秩不足第二多个时刻之间的相关性太强导致行向量线性相关第三用经验ACF塞进相关矩阵时估计噪声太大矩阵已经偏离正定。处理方式通常是特征值裁剪eigvals, eigvecs np.linalg.eigh(R) eigvals[eigvals 1e-8] 1e-8 R_fixed eigvecs np.diag(eigvals) eigvecs.T L np.linalg.cholesky(R_fixed)注意裁剪量不要设太大否则会把相关矩阵拉向单位阵等于削弱了时序相关。另外用Spearman秩相关矩阵代替Pearson相关矩阵通常会更稳定因为秩相关对异常值和分布形态不敏感。6.2 削减后出现0概率场景怎么处理场景削减完成后清理阶段经常发现一部分代表场景的概率权重几乎为0。比如快速前向选择选出了一个孤立场景但它周围没有任何原始场景归属于它最终权重小到可以忽略。这些场景留着对目标函数没有贡献却白白占用模型中的场景索引和决策变量数量。处理办法很直接削减后把所有概率小于1e-8的场景剔除再把剩余概率重新归一化。如果剔除的恰好是极端场景需要额外确认它是真的概率极低还是算法误删导致的假象。前者可以直接删后者要考虑手动保留。6.3 K-means把极端场景平均成“普通场景”用K-means做削减最容易出现的问题是低出力场景和高出力场景被分到同一个簇后中心点落到中间位置变成一条“普通场景”。在充裕性评估里这种平均化会直接掩盖失负荷风险让人误以为系统很安全其实只是极端场景被算法抹掉了。对策有三个按优先级排列换用K-medoids代表点必须是真实场景按出力分位数把场景分层在每个层内分别削减保证极端层有代表场景幸存在特征向量或距离函数里加入极值惩罚项让远离平均状态的场景更难被合并。6.4 样本规模与削减数目的经验法则最后说说规模问题。生成多少条MC场景、削减到多少条这两个数字没有统一标准但有一些经验法则可以省去反复试错。生成场景数量通常是削减后数量的50到200倍。比如最终目标20条先生成1000到4000条目标200条先生成1万到4万条。太少原始概率测度覆盖不足太多削减阶段的计算开销成倍增长。削减后的数量可以通过“场景数-目标值”曲线确定从5、10、20、50、100逐步增加当目标值变化小于0.5%就说明已经逼近收敛不需要再加场景。多风电场联合场景的时序结构比单场复杂得多需要保留的场景数也会成倍增加这一点在做方案设计时要提前留出余量。最后分享一个我自己的操作习惯每次项目开跑前哪怕求解时间再长我也会先用几百个场景做一次全场景优化拿到基准目标值之后再用削减后的场景集做同样的优化。如果两者偏差超过1%我不会急着去调削减数量而是先回头检查距离度量和相关矩阵因为这种偏差通常意味着时序相关性在生成或削减阶段就已经丢了。这个习惯帮我拦下了不少“看着削减效果很好、一进优化模型就翻车”的情况。