☰
蒙特卡洛场景生成中的时序相关性建模与场景削减实战解析
2026/10/14 4:32:04 网站建设 项目流程

讲个实际经历。我之前在模拟项目X里做含新能源的随机优化调度,第一步就是把风光出力的不确定性变成一组可计算的场景。当时直接用蒙特卡洛(MC)抽样,独立地对每个时段采样,生成上千条日出力曲线,结果第一眼就觉得不对劲:相邻两个小时的出力像抽签一样高低乱跳,上午还在满发,下午就掉到接近零,完全没有真实出力那种"惯性"。那段研究最后基本是从头返工。问题不在于MC本身,而在于单纯按边缘分布抽样完全丢了时序相关性。后来我把"时序相关性"显式建模进MC场景生成流程,再用场景削减算法压缩场景数量,整个研究才真正能落地。这篇就把整个思路、数学工具、完整流程和踩到的坑梳理一遍。

1. 场景生成里最容易被忽视的"时间箭头"问题

1.1 先搞清楚MC场景生成到底是干什么的

随机优化里最常见的需求就是:不确定性参数(风电出力、光伏功率、负荷曲线)不是一个数,而是一条随时间变化的曲线。一天24个时段,就有24个随机变量,每个场景就是一条24维的时间序列向量,一条"出力轨迹"。蒙特卡洛场景生成做的事情,本质上就是从历史数据中估计出概率结构,然后大量采样,凑出一个能代表未来可能性的场景集合。

听起来直接,但这里有个很容易忽略的前提:采样来的场景必须和真实过程具有相同的时间结构。真实风电场出力曲线上,相邻时段之间有明显的依赖关系。凌晨的风速往往持续稳定,不会一个小时内从满发掉到停机;白天的光伏出力更像一条平滑的单峰曲线,不可能上午10点达到峰值后11点直接变零。

如果只考虑每个时段的边缘分布,用MC对各时段独立抽样,生成出来的场景集合虽然每个时段的统计均值、方差都对得上,但时序上完全失真。这样的场景集喂给优化模型,最直接的后果是调度算法会去"钻空子",例如利用场景中相邻时段出力突变去调整储能策略,计算出的备用量、爬坡需求要么过于乐观要么过于保守,根本没法用于实际决策。

1.2 独立采样为什么必然失真:用一个简化例子说明

假设某个风电场的归一化出力服从均值0.35、方差0.08的正态分布,且相邻时段的一阶自相关系数是0.92。如果忽略时序相关性,两个相邻时段的联合分布被当成独立变量处理,那么真实世界里P_t和P_{t-1}同高同低的概率远大于独立情况下计算的概率。

直观对比一下:

情况P_t = 0.8时,P_{t+1}>0.6的条件概率
真实过程(ACF=0.92)约85%
独立采样(ACF=0)约38%

可以看到,一旦忽略时序相关性,条件概率会直接从85%跌到38%。也就是说,独立采样生成的场景中,有六成以上的"高出力持续到下一时段"的实际情况会被漏掉,取而代之的是大量离谱的高-低切换。这个差异会穿过场景削减环节,最终进入优化模型的约束条件,影响风电渗透率较高的系统中的旋转备用容量、机组组合结果和储能充放电策略。

所以,"时序相关性"不是锦上添花,而是MC场景生成里绕不过去的一环。要么在采样前注入相关结构,要么在建模时就采用天然具备时序依赖的方法。

2. 时序相关性的数学抓手:协方差结构、自回归与转移概率

2.1 从ACF到滞后协方差矩阵:量化"高点之后还是高点"

要建模时序相关,先得有量化的工具。最常用的是自相关函数(ACF),定义为滞后k的相关系数:

ACF(k) = E[(P_t - μ)(P_{t-k} - μ)] / σ²

其中μ是均值,σ²是方差。ACF(k)落在[-1,1],风电出力数据的典型特征是:lag-1很大(0.8以上),随后随滞后阶数缓慢衰减。光伏出力更特殊,由于昼夜周期性,ACF还带有明显的24小时节律,但这不影响下面的建模思路。

从单变量ACF延伸到多时段,就得到滞后协方差矩阵Σ。对24时段曲线来说,Σ是一个24×24的对称矩阵,第i行第j列是时段i和时段j的协方差。这个矩阵完整刻画了所有时段两两之间的线性关系,是后面注入相关性的核心对象。

构造Σ有两种方式:

  • 直接从历史数据算样本协方差矩阵;
  • 或者先拟合ARMA等时间序列模型,从模型参数推导出理论协方差结构。

实际项目里我建议第一种,因为直接、透明,还能保留数据固有特征。

2.2 用Cholesky分解把相关性注入随机样本

有了Σ,下一步就是在MC采样时把相关性带进样本。标准做法是对Σ做Cholesky分解,得到下三角矩阵L,满足:

Σ = L Lᵀ

然后算法流程就非常简洁:

  1. 生成一个24维独立标准正态随机向量Z,Z~N(0, I);
  2. 计算X = μ + L Z;
  3. 将X从标准正态分位空间逐步变换回原始出力的概率空间。

检查一下:Cov(X) = E[(LZ)(LZ)ᵀ] = L Lᵀ = Σ,相关结构精确还原。这就是为什么Cholesky分解是这类任务的首选:把复杂的相关采样问题,变成一次矩阵分解加一次矩阵乘法。

核心代码用Python写也就十来行:

import numpy as np # 假设 mu 为 24 维均值向量,sigma 为 24x24 协方差矩阵,n_scen 为场景数 mu = power_mean sigma = power_cov # 1. Cholesky 分解,使 sigma = L @ L.T L = np.linalg.cholesky(sigma) # 2. 生成独立标准正态样本 Z = np.random.normal(size=(24, n_scen)) # 3. 注入相关性后的样本 X_corr = L @ Z + mu.reshape(-1, 1) # 4. 把标准正态值映射到 0~1,再用历史分布的分位数逆变换回实际功率 samples_uniform = norm.cdf(X_corr.T) # 形状变为 (n_scen, 24) scenarios = power_quantile_func(samples_uniform)

有一个前提必须说明:上面的操作假设变量服从联合正态或能以Nataf变换近似。真实风电出力往往偏态、含大量零出力时段和中高段厚尾,直接硬套正态会低估极端场景。稳妥的做法是先用正态分布建模秩相关,然后通过分位数逆变换映射回真实边缘分布,形成"正态Copula + 经验边缘分布"的组合。这比单纯假设联合正态稳健得多,也是我在模拟项目X里最终采用的方式。

2.3 马尔可夫链视角:一个与蒙特卡洛互补的替代建模

标题里"MC"我默认指蒙特卡洛(Monte Carlo),但对时序相关性研究来说,马尔可夫链(Markov Chain)其实是另一个强力且互补的建模角度,很多研究里也把两者结合用。

马尔可夫链的思路是把出力水平离散化为若干状态。比如按标幺值把出力分成5档:0~0.2、0.2~0.4、0.4~0.6、0.6~0.8、0.8~1.0。然后从历史数据统计状态转移概率矩阵P,P[i][j]表示当前时段处于状态i时下一时段转移到状态j的概率。生成场景时,从初始状态出发,按转移矩阵逐时段抽样状态,就能得到一整条状态序列。

这种方法的优势在于:

  • 不需要预设联合正态分布,能表达非对称、非线性的时序依赖;
  • 对"出力卡在0附近持续好几小时"这类真实风况,马尔可夫模型能自然描述;
  • 转移矩阵含义直观,后续可以通过调整矩阵来反映季节、天气类型等。

但缺点也明显:状态数一旦增多(比如分成20档),需要估计400个转移概率,历史数据不够时矩阵会布满零值或噪声;另外离散化会丢失状态内的细微差异。所以实际项目里,如果风速分布复杂且设备资源够,可以考虑马尔可夫链蒙特卡洛(MCMC)或隐马尔可夫模型做进一步扩展,入门阶段先把协方差+Cholesky和马尔可夫状态转移这两条路摸熟,就已经能覆盖绝大多数需求了。

3. 场景生成与削减的完整实验流水线

3.1 数据准备:清洗、归一化与参数估计要点

场景生成这个环节,数据质量比算法选择更影响结果。第一步是处理历史序列,有几个细节容易被新手忽略:

  • 剔除限电时段的数据。部分时段出力曲线被人为调度截断,不代表资源的真实随机特性,混进去会让协方差矩阵失真;
  • 缺失值的处理。我习惯用前后插值法补短缺口,超过连续6小时以上的长缺口直接丢弃对应日期段,不要硬填;
  • 归一化。按装机容量把功率归一到0~1标幺值,这样不同规模的场站可以横向比较;
  • 颗粒度统一。15分钟数据和1小时数据混用是常见错误,建议统一到优化模型所需的颗粒度再进入参数估计。

数据清理完后,估计均值向量mu和历史协方差矩阵sigma。如果时间序列有明确季节性,最好按季节或月份分桶估计,否则夏季和冬季样本混在一起,协方差结构会被"平均"掉,生成场景的季节特征不清晰。

3.2 生成规模怎么定:为什么先搞2000个而不是直接10万个

场景生成有个朴素的直觉:样本越多概率分布越准确。但落到实操,MC生成本身很快,瓶颈在后面的场景削减和随机优化求解。

假设用Cholesky注入相关,生成20万个24维场景也就几十秒的事,但后续的快速前向选择或同步回代削减,其复杂度在小样本集上尚可接受,一旦样本数上万,距离矩阵计算量会迅速失控。削减算法通常需要计算所有场景两两之间的距离,N个场景形成的距离矩阵是N×N,N=20000时这个矩阵需要存储约3GB浮点数,计算量相当可观。

我的习惯是先设N=1000~2000,跑通全流程,验证削减效果和优化模型收敛性,再视需要扩展到更大规模。除非你明确要做高精度尾部分析,否则对大多数调度研究来说,生成2000个场景、削减到10~20个代表性场景已经完全够用,且可控性强。

3.3 削减算法选择:同步回代削减与快速前向选择的取舍

场景削减的目标是:从N个原始场景中选出一个数量更少的子集,并重新分配概率,使削减后的场景集与原场景集在概率分布意义上尽量接近。

我在模拟项目X里对比过三种主流方法:

同步回代削减(SBR)

逻辑很像"逐步淘汰制"。每次迭代寻找一个对整体概率距离贡献最小的场景,删除它,并把它携带的概率加到距离它最近的保留场景上,直到场景数降到目标值。因为它每一步都在做全局最优的近似删除,保留下来的场景都是原始场景本身,没有人为"平均"过的痕迹,时序特征保留得比较完整。缺点是大N下迭代速度偏慢。

快速前向选择(FSS)

和回代相反,它是逐步"招人"。起始集合为空,每一步从剩余场景中选出一个,使当前已选集合对全集的距离下降最多,直到选满目标场景数。FSS的优点是当目标场景数较少时,每次选择的边际收益比较明显,不容易出现SBR后期"删谁都差不多"的模糊局面。

K-means聚类

把场景聚类后取每类质心作为典型场景,每个类的样本占比作为概率。实现最简单、速度快,但质心是"平均产物",会抹平时序细节,而且质心不一定是物理上可实现的一条真实出力曲线。如果目标场景数比较大,比如30个以上,K-means的误差常常会小于SBR;但目标场景数很小(5~10个)时,我还是优先选SBR或FSS。

简化版的SBR代码思路如下:

def sbr(scenarios, probs, target_num): # scenarios: (N, 24) 原始场景 # probs: (N,) 概率,通常初始为 1/N survivors = list(range(len(scenarios))) p = probs.copy() while len(survivors) > target_num: dist = compute_pairwise_distances(scenarios[survivors]) best_cost = np.inf best_idx = None for i in range(len(survivors)): nearest_substitute = min_dist_for_deleted[i] # 按距离矩阵预先算 cost = p[survivors[i]] * nearest_substitute if cost < best_cost: best_cost = cost best_idx = i deleted = survivors.pop(best_idx) # 找被删场景最近的保留场景 closest = nearest_survivor(deleted, survivors) p[closest] += p[deleted] return survivors, p

真实性能和收敛性取决于距离度量。风电场景我常用欧氏距离或加了时段权重的加权欧氏距离,也有人用动态时间规整更好地对齐时序形态,但计算量更大,权衡后我多数时候用加权欧氏距离。

3.4 削减后的场景概率怎么重新分配

场景削减完成后的概率重分配是很容易出错的一环,尤其在FSS里要特别小心。

  • SBR的做法:删除场景后,把被删概率累加到最近的保留场景上,最后所有保留场景的概率和为1,不需要额外归一化;
  • K-means做法:直接用每个类中样本数占总样本数的比例作为该场景概率;
  • FSS做法:严格来说,每个保留场景的最终概率应该是它所能代表的子集内原始场景概率之和。实现时通常先对每个原始场景找到它在保留集中的"最近代表",然后汇总概率。

实际操作中我建议算完概率后立即加上断言检查:sum(new_probs) - 1.0 < 1e-9。别觉得多余,聚类实现里概率分配的逻辑一旦写错,后面整个优化问题的目标函数都会偏。

4. 削减效果评估:不止看分布接近,还要看时序特征没丢

4.1 Kantorovich距离和Wasserstein距离到底说明了什么

评估削减质量,一个基本指标是Kantorovich距离,也是Wasserstein距离的一种形式。直观理解就是:要把原分布变成削减后分布,所需的最小"搬运成本"。如果两个分布本身很像,搬运量很小;如果削减后的场景分布和原分布差距很大,搬运成本就高。

公式化一点,单场景之间的距离度量d(ξ₁, ξ₂)先定义好,Kantorovich距离就是所有场景对之间运输量的最优分配成本。实际使用不一定要手动实现它,很多场景削减库直接内置了这个计算。

但必须提醒:Kantorovich距离只看全局分布接近程度,并不能单独衡量时序结构保持度。两个场景集可能在Wasserstein距离上很接近,却一个保留时序惯性、另一个剧烈波动。所以削减去留质量必须同时看两个维度:分布距离 + 时序统计量误差。

4.2 一个验证实验:1000个场景削减到10个,时序特征还剩多少

我在模拟项目X里做过一组对比实验,用某个虚拟风电场全年8760小时的数据,归一化后均值约0.35、方差约0.08,ACF(1)=0.92。先用正态Copula方案生成1000个场景,再用三种方法削减到10个场景,结果见下表。

削减方法Wasserstein距离平均ACF(1)误差95%分位数误差计算耗时
SBR0.0210.0314.8%420秒
FSS0.0180.0453.6%360秒
K-means0.0260.1129.2%25秒

看Wasserstein距离,FSS表现最好;但看ACF误差,K-means比其他两种方法高出一大截,几乎抹掉了时序惯性。原因是K-means把同类场景平均成质心,时间上的陡升陡降、持续爬坡这类形态全被"平均"成了光滑曲线。如果把削减后的场景直接喂给储能优化,K-means方案给出的充放电策略会在真实场景下明显高估可调节空间。

这个实验给我的经验很直接:

  • 场景数较少时,SBR/FSS的时序结构保持能力显著优于K-means;
  • 场景数较多时,K-means的分布逼近能力可以反超,且速度快得多,适合大规模预筛选;
  • 最终的评估报告里,ACF误差和95%分位数误差必须与Wasserstein距离并列展示,否则决策者只看分布接近程度很容易被误导。

4.3 除了ACF,还要盯哪些时序指标

ACF是最基础的,但不够。实际调度对时序的敏感点比ACF更隐蔽,还需要关注:

  • 相邻时段最大变化量:就是ΔP = |P_t - P_{t-1}|的分布。风电爬坡事件往往对应着ΔP超过某个阈值(如额定容量的10%),削减后如果这个分布变了,备用容量和爬坡约束的评估结果就不可靠;
  • 持续低出力时段长度:夜里连续4小时出力接近0的场景,在削减后如果都被替换成"夜间平稳出力0.2",储能系统的充电策略会被误导;
  • 多场景间的相关结构:如果项目包含多风电场,还需检查各场站之间的空间相关性是否被削减破坏,削减是逐点作用于所有场站联合场景的,不应该分别对每个场站单独削减再拼接,那会让联合相关结构崩塌。

5. 实操落地中的避坑经验与后续扩展

5.1 削减场景数量选择的试错经验:看"肘部曲线"

场景数削减到多少最合适,没有一个普适答案。直接拍脑袋选10个,很可能要么分布误差过大,要么优化时间过长。比较实用的做法是:固定原始场景集,依次削减到5、10、15、20、30个,画出"保留场景数 vs Wasserstein距离"曲线,找"肘部"。

我实际跑过很多次,场景数从5增加到15时,Wasserstein距离通常快速下降;而从15增加到30时,距离改善就很平缓了。这个拐点意味着再增加场景数量,分布精度提升已赶不上计算成本的上升。这比固定场景数科学得多,而且能向外行解释得清楚。

5.2 时序相关性的损失会在优化调度里"传导放大"

这是一个重点。很多人以为场景削减造成的误差最多让目标函数值偏几个百分点,但实际在含储能的机组组合问题里,时序相关性的损失会被放大。

举个我在模拟项目X里遇到的例子:削减后场景数从20降到5,ACF(1)误差从0.03上升到0.12,分布式能源渗透率约45%。同样一个优化问题,目标函数值只变化了大约7%,但最优解中储能的自循环充电次数增加了25%——算法在多出的"伪场景跳变"里发现了套利空间,用短时充放去吸收那些在真实风况下根本不存在的出力波动。

结论是:如果优化结果对场景时序很敏感,一定要把削减前的验证场景集单独保留一份,在最优解出来后做后验回放:把削减场景求出的决策放到原1000个完整场景上评估期望成本,观察是否有异常。这个步骤成本很低,却能提前暴露时序失真问题。

5.3 几个容易踩的统计错误和对应的处理方式

协方差矩阵不正定或接近奇异

当两个时段出力几乎完全线性相关时,Σ的行列式趋近于零,Cholesky分解会直接报错。处理办法是先对Σ做特征值分解,把过小的特征值截断或加上一个很小的正则项εI,通常ε取1e-6即可。

直接对功率采样,导致小于0或大于1的物理不可行值

如果用了高斯假设但没做分位数逆变换,削减后的场景可能出现负出力。解决方式要么改用截断正态,要么像我前面说的,用正态Copula方案,把采样拉到均匀空间再做物理边界映射。

固定随机数种子

MC场景生成有随机性,同一套参数每次运行出来的场景集可能不同。如果实验对比不同削减方法时没有固定种子,很难判断结果差异来自算法还是来自采样波动。我习惯在实验开始就固定np.random.seed(42),并且把每个阶段的种子记录到日志里,方便后续复现。

5.4 后续扩展:时空相关、在线生成与极端场景增强

这套方法本身有很好的扩展性。

  • 多场站时空相关性:把多个风电场、光伏电站的出力拼成一个高维向量,协方差矩阵变成块状结构,Cholesky分解依然适用,只是矩阵规模变大,需要通过空间相关模型或采样降维方法控制计算量;
  • 在线滚动场景生成:把时序相关性参数按滑动窗口在线更新,生成未来24小时场景集,供模型预测控制使用,只要每次估算Σ后重新做Cholesky即可,计算量完全可以接受;
  • 极端天气增强:研究可靠性时常规样本往往不够极端。可以在场景生成阶段对尾部区域做重要性采样或条件采样,人为增加极低出力、极高出力场景的占比,再通过场景削减保留部分极端典型场景。

最后分享一个我个人很依赖的小技巧:场景削减的目的从来不是"把数量变少",而是"让少数场景尽可能多地承载原场景集的决策相关信息"。所以在完成削减后,别急着跑优化,先花十分钟把削减前后的典型日曲线画出来对比一下。如果削减后的10条曲线里,有哪一条看起来明显"假"——比如前后时段跳变幅度远超历史数据范围,那就说明相关结构在生成或削减环节出了问题,及时回头修正参数,比后期分析优化结果再返工要省事得多。

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

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

立即咨询