考虑时序相关性的蒙特卡洛场景生成与削减实战指南
2026/9/8 3:24:42 网站建设 项目流程

先解释一个容易让人走错片场的词: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可以分解为R=L·L^T,那么Y=L·U的协方差矩阵就是R。也就是说,只要我们能构造出描述时序相关性的相关矩阵R,就能通过线性变换把独立的随机数变成带相关性的随机数。

具体操作分几步走:

  1. 根据历史数据的ACF估计滞后相关矩阵R。常用做法是假设指数衰减结构,例如R[s,t]=ρ^|s-t|,ρ是滞后1阶自相关系数;
  2. 对R做Cholesky分解,得到下三角矩阵L;
  3. 生成T×N的独立标准正态矩阵U;
  4. 计算Y=L·U,得到带相关性的正态样本;
  5. 对Y的每个元素求标准正态CDF,得到(0,1)均匀分布的样本;
  6. 再用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中每个位置的数值。具体来说:

  1. 用Cholesky生成秩相关结构正确的正态样本Z;
  2. 对Z的每一列做排序,记录每个位置的秩次;
  3. 从目标边际分布中抽取样本(或者从历史数据中重采样),也做排序;
  4. 把目标分布排序后的值,按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, key=lambda 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%,我不会急着去调削减数量,而是先回头检查距离度量和相关矩阵,因为这种偏差通常意味着时序相关性在生成或削减阶段就已经丢了。这个习惯帮我拦下了不少“看着削减效果很好、一进优化模型就翻车”的情况。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询