考虑时序相关性的蒙特卡洛场景生成与削减：理论、实操与避坑指南

1. 场景生成的研究动机与整体设计思路

先说一个我自己的体会：在电力系统里做随机优化、概率潮流、可靠性评估这些方向，绕不开一个核心问题——你手里拿到的风、光、负荷数据，本质上是一条条实测的时间序列，但优化模型需要的是“未来可能发生的一批情况”，而不是“某一天具体发生了什么”。这批“可能发生的情况”在学术圈和工程圈都叫场景。场景生成和场景削减，说白了就是把连续的概率分布变成有限个有代表性的离散样本，让后续的优化、评估、决策能在可接受的计算代价内跑起来。

1.1 为什么要用蒙特卡洛方法生成场景

蒙特卡洛方法在这个领域里之所以是主流，核心原因是它的逻辑足够朴素、足够通用。给定风电场历史出力数据，你可以先拟合出风速或者出力的概率分布，然后从这个分布里反复抽样，抽出一大批可能的出力曲线。每一条曲线就是一个场景，场景数量可以做到几百上千，这比单纯拿历史数据“回放”要灵活得多——因为你能生成出没出现过但概率上完全可能出现的组合。

但纯蒙特卡洛有个很明显的问题：它默认每次抽样是独立的。这在某些场景下没问题，比如不考虑时间先后顺序的静态不确定性分析。可一旦你的研究对象是“未来24小时的风电出力曲线”，独立抽样就坏了事——你抽出来的场景之间完全没有时序关联，可能出现前两天风速低、第三天突然跳到满发的离谱曲线，这种场景虽然单点概率分布是对的，但整条曲线在时序上站不住脚。所以就有了“考虑时序相关性”的MC场景生成，这也是这个标题里最核心的增量点。

1.2 场景生成与削减的完整技术链路

我习惯把整套流程拆成四步，每一步都有明确的输入输出：

• 第一步：数据清洗与分布拟合。拿到历史出力数据，做异常值剔除、缺失值补齐，然后拟合每个时段或者某个时间粒度的概率分布。如果是风速，Weibull分布通常是首选；如果是光伏出力，Beta分布用得比较多。
• 第二步：相关性建模。这一步是“考虑时序相关性”的落点。常见做法有Cholesky分解后对独立样本做相关性变换、Copula函数建模、或者用马尔可夫链直接建模时段间的转移关系。
• 第三步：蒙特卡洛抽样生成初始场景集。大规模抽样，通常是1000到5000个场景，覆盖完整的时序长度。
• 第四步：场景削减。从几千个场景里挑出几十个有代表性的，同时保证削减前后的概率分布特性、统计特征尽量接近。

这四步是一条完整链路，缺一个环节在实际工程里都会出问题。很多初学者只关注第三步的抽样技巧，忽略了第一步分布拟合的质量和第四步削减算法对最终结果的影响——事实上这两个环节往往决定了整个方案的成败。

2. 考虑时序相关性的蒙特卡洛场景生成实现细节

时序相关性之所以要专门拎出来讲，是因为“独立采样、逐点拼接”生成的伪时序场景，在优化模型里会产生系统性偏差。举个最简单的例子：储能调度模型里，如果生成的风电场景是剧烈波动的锯齿状曲线，储能就会被迫频繁充放电，算出来的调度策略和成本评估会明显偏离真实情况。

2.1 基于历史序的时序相关性刻画方式

处理时序相关性，我见过三种工程上真正可用的路线，它们适用的场景和实现代价差别很大。

第一种是Cholesky分解加排序样本法。它的思路是先估计出各时段之间的相关系数矩阵，然后对独立抽样得到的样本矩阵做线性变换，使得变换后的样本满足预设的相关系数结构。具体操作不复杂：假设你要生成T个时段的风电出力场景，先根据历史数据算出T乘T的相关系数矩阵R，做Cholesky分解R = L L^T，然后对独立标准正态样本Z做变换X = L Z，X就带上了R所描述的相关结构。这个方法实现简单、计算量小，但有个前提假设——各时段的边际分布必须能转换到正态空间，实际中一般用分位数变换来处理。

第二种是Copula方法。本质上是在说“相关性结构和边际分布可以分开建模”。你可以用历史数据分别估计出每个时段的边际分布，再用一个Copula函数去描述时段之间的相关结构，最后通过Copula采样生成具有特定相关性的样本。工程上最常用的是Gaussian Copula和t-Copula，前者尾部和正态类似，后者能捕捉更极端的联动情况。这个方法理论上更严谨，但实现复杂度也上了一个台阶。

第三种是马尔可夫链方法。把时序出力划分成若干状态，通过历史数据统计状态转移概率矩阵，然后按转移概率一步步生成未来的状态序列。这个方法对强波动、强非平稳的数据效果出奇地好，特别是在风速这种有明显“持续惯性”的序列上。缺点也很明显——状态划分的粒度直接影响场景质量，状态太多转移矩阵稀疏，状态太少又丢失了细节。

2.2 分布拟合的实操要点

无论选哪种相关性建模方法，第一步的分布拟合都绕不开。这里我踩过几次坑之后总结出的经验是：

• 不要只做单一分布的拟合优度检验。比如风速，Weibull分布是首选，但不同风电场、不同高度的风速分布形态差异很大，有时候Gamma分布、对数正态分布反而拟合得更好。实际操作中我会同时拟合三四类候选分布，用AIC、BIC或者KS检验去挑最优的，而不是拍脑袋直接上Weibull。
• 注意零出力时段。风电里零出力比例很高，光伏夜间出力就是严格的零。这类“零膨胀”特性会让连续分布在零点的拟合效果非常差。我的做法是拆成两个部分建模：先建模“是否为零出力”的二值过程，再对非零出力部分建模连续分布，最后组合成完整的混合分布。
• 时段粒度这个参数容易被忽略。做24小时场景，你可以对每个小时分别拟合一个分布，也可以把全天作为一个整体去拟合。前者更精细但需要的数据量更大，后者简单但对日出、日落这类强时序变化表达不够。折中方案是分时段拟合，每3到4个小时共用一个分布参数，能显著降低参数估计的方差。

2.3 基于蒙特卡洛时序抽样的工程实现

工程实现层面，我给出一个可以直接参考的标准流程，以某风电场96点出力场景生成为例：

1. 取历史出力数据，按每天96个点（15分钟一个点）整理成一条条完整曲线。
2. 对每个时段（这里就是96个时段）分别做分布拟合，得到各时段的边际分布函数。
3. 计算96个时段间的相关系数矩阵。这里注意，直接用原始出力数据算相关系数会被异常值干扰，建议先做一次分位数变换或者用Spearman秩相关系数。
4. 生成独立标准正态随机矩阵Z，维度是96乘N，N是抽样数量（比如2000），用Cholesky分解后的相关矩阵做线性变换，得到带相关性的标准正态样本。
5. 用分位数变换把标准正态样本映射到各时段的实际分布上，得到最终的场景矩阵。

这套流程在Python里实现，核心代码量不超过50行，用numpy和scipy就够。实测生成2000个96点场景的耗时在秒级，效率上完全没有瓶颈。但要注意的是，Cholesky分解要求相关系数矩阵必须是正定矩阵，实际数据算出来的矩阵往往只是半正定，甚至可能出现微小的负特征值，这时候需要做特征值修正或直接用Nearest Positive Definite Matrix近似，否则分解直接报错。

3. 场景削减的核心算法与选型分析

场景削减这一步，直接决定后续优化模型的计算规模。你不可能把2000个场景全部塞进一个混合整数规划模型里跑，场景数量一旦超过100个，求解时间就会指数级增长。削减的目标是用最少的场景数量，保留原始场景集的主要概率信息。

3.1 同步回代消除法的实现逻辑

工程上最常用的削减算法是快速前向选择法和同步回代消除法。我个人用得最多的是同步回代消除法，它的逻辑非常直观：每轮迭代都删掉一个对整体场景集合“影响最小”的场景，然后把它的概率累加到离它最近的留存场景上，直到场景数量削减到目标值。

具体步骤我拆开说：

1. 计算所有场景两两之间的距离，这个距离可以是欧氏距离，也可以是考虑时序形态的动态时间规整（DTW）距离。对于96点出力曲线这种场景，我建议用欧氏距离就够，DTW虽然能处理时间轴上的形变，但计算量大很多，而且削减结果未必更好。
2. 对每个场景，找到离它最近的场景，记录这个最近距离。
3. 每轮迭代中，找出具有最小最近距离的场景，把它删除，同时把它的概率加到那个“最近场景”上。
4. 重新计算距离矩阵中受影响的部分，继续下一轮迭代，直到满足目标场景数。

这个算法的复杂度是O(N^2 M)，N是初始场景数，M是削减后的场景数。当N等于2000、M等于20时，计算时间在几秒到几十秒之间，属于完全可接受的范围。

3.2 基于聚类的场景削减方案对比

除了同步回代消除法，聚类类方法在工业界也用得很广。K-means聚类、K-medoids聚类、层次聚类都能用，其中K-medoids比K-means更适合场景削减——因为K-medoids选出来的簇中心是实际场景中的某一条曲线，而不是所有样本的平均曲线。平均曲线问题是K-means在场景削减里的硬伤：它生成的“虚拟场景”可能平滑得失真，丢了真实场景中该有的波动特征。

层次聚类和K-medoids的削减效果接近，但层次聚类不需要预先指定K值，可以通过树状图观察不同K值下的分簇质量，适合做敏感性分析。不过层次聚类的计算复杂度偏高，场景数超过1000时就有点吃力了。

在实际项目里，我的做法是把同步回代消除法和K-medoids都跑一遍，对比削减后场景的统计特征——均值、方差、相关系数矩阵——哪个与原始场景集更接近就采用哪个。大多数情况下同步回代消除法的表现更稳定。

3.3 削减质量的评价指标

很多资料讲到这里就不往下讲了，没有说清楚“怎么判断削减得好不好”。实际上，评价削减质量是有标准动作的，我常用以下三个指标：

• 概率分布距离：削减前后场景集的累计分布函数之间的KS距离或Wasserstein距离，这个指标反映整体分布保真度。
• 统计特征偏差：对比均值曲线、标准差曲线、相关系数矩阵的偏差。均值偏差控制在一两个百分点以内算合格。
• 目标函数值偏差：如果你的场景最终要喂给优化模型，最硬核的检验方式是分别在削减前和削减后求解同一优化问题，对比目标函数值的差异。如果差异在可接受范围内（比如1%以内），说明削减没有引入不可忽略的偏差。

第三个指标其实很多论文都会用，但我发现实践中的初学者很少主动去做这个验证。这很可惜——因为哪怕前两个指标都很好看，最终优化结果也可能因为场景削减导致了可行域变化而出现偏差。养成“最终验证”的习惯，能让你的场景削减方案在工程评审时更有说服力。

3.4 场景数量的选取策略

场景数量选取是个典型的权衡问题。设得越多，精度越高，但计算代价随之上升。我一般按以下经验规则来定初始值：

• 做规划类问题（比如容量规划、网架扩展），场景数量取10到20个就够。
• 做运行调度类问题（比如日前机组组合、经济调度），场景数量取20到50个。
• 如果涉及极端场景的捕捉，比如可靠性评估、韧性分析，建议在常规场景之外额外保留少量低概率高风险场景，这时候总量可以放到50到100个。

一个比较聪明的做法不是直接定死场景数量，而是做削减数量扫描：从10个到100个，每间隔10个做一次削减，记录Wasserstein距离和优化目标函数值的变化曲线。你会发现存在一个明显的拐点，拐点之后再加场景数量，精度提升就非常有限了。用这个拐点去定场景数量，既科学又高效。

4. 一次完整的场景生成与削减实操全流程

前面讲了不少原理和方法，这一节我完整走一遍实操流程，用一套具体数据说话。为了让复现方便，我用公开的某风电场出力数据来做示例，处理工具是Python加numpy、pandas、scipy、scikit-learn。

4.1 数据准备与分布拟合

我取了一年的历史出力数据，时间分辨率是15分钟，也就是每天96个点。原始数据有缺失和异常值，先做两步清洗：

• 剔除出力大于装机容量1.05倍的点（测量误差导致）。
• 对缺失值做线性插值，如果连续缺失超过4个点，直接用前一天同时段的数据填充。

清洗完成后，按每天96个点的格式重排成365条曲线。然后对96个时段分别做分布拟合。考虑到风电出力的零膨胀特性，我用“零出力概率+非零出力分布”的混合模型。非零出力部分我试了Weibull、Gamma、对数正态三种分布，用AIC做选择，结果大部分时段的非零出力更适合Gamma分布。这一步的结果直接决定后续抽样的准确性，所以值得花时间细调。

4.2 时序相关性与蒙特卡洛生成

清洗和分布拟合完成以后，进入核心的时序场景生成环节。我先生成2000个场景作为初始场景集，具体步骤如下：

1. 计算96个时段之间的Spearman秩相关系数矩阵R。
2. 对R做特征值修正，确保它是正定矩阵，然后做Cholesky分解。
3. 生成一个96乘2000的标准正态随机矩阵Z。
4. 做线性变换X = L Z，让X带上相关性结构。
5. 对X每一行做标准正态CDF变换，得到均匀分布样本，再通过各时段边际分布的逆CDF变换，映射到实际的出力数值空间。
6. 对零出力概率做一次伯努利抽样，把一部分样本的出力置零。

跑完后，我检查了生成场景的相关系数矩阵和历史数据的相关系数矩阵，两者的Spearman相关系数热图高度一致，均值偏差在1%以内。这说明这一套“边际分布+秩相关+分位数变换”的组合拳是有效的。

4.3 场景削减与结果对比

接下来做场景削减。我分别用同步回代消除法和K-medoids聚类跑了一遍，削减目标设定为20个场景。削减后的场景质量对比如下：

可以看出K-medoids在计算速度上有微弱优势，但削减精度明显不如同步回代消除法。这和我之前提到的判断一致——同步回代消除法在保真度上更可靠。实际操作中，如果场景数在1000到3000之间，我会优先选同步回代消除法；如果场景数超过5000，会先用K-medoids粗削减到1000，再用同步回代消除法做精细削减，这样兼顾了计算效率和削减精度。

最后，我把削减后的20个场景喂进一个日前机组组合模型，对比用原始2000个场景求解得到的目标函数值，偏差控制在0.6%，说明削减后的场景集可以放心用于后续的随机优化计算。

5. 实操中绕不开的坑与排查思路

不管论文里写的流程多么顺滑，实际工程里总是会遇到各种幺蛾子。我把这几年在这个方向踩过的坑、还有各路同行交流时提到的高频问题整理一下，直接给你一份避坑指南。

5.1 相关系数矩阵非正定的处理

这是Cholesky分解最经典的问题。历史数据的相关系数矩阵理论上应该半正定，但实际计算中由于数值误差、数据缺失、样本量不足，经常得到非正定矩阵，Cholesky分解直接崩溃。

处理方式有两种。第一种是特征值裁剪：对矩阵做特征值分解，把负特征值置为0或者一个很小的正数，再用特征向量重建矩阵。这个方法的缺点是改变矩阵幅度，导致相关结构被扭曲。第二种是寻找最近的正定矩阵，工程上常用的是Higham算法，scipy里直接有scipy.linalg.nearby_psd可以用。我实测下来第二种效果好很多，保留下来的相关性信息更完整。

5.2 生成场景的边际分布失真

做完相关性变换后，有时你会发现生成场景的直方图和你拟合的目标分布对不上。这个问题的根源出在步骤顺序上。如果你先对边际分布做抽样，再做线性相关变换，变换过程本身就会破坏边际分布。正确的做法是先在标准正态空间做所有变换，最后一步才做分位数逆变换回到物理空间。这个顺序不能反，反了之后的结果错得毫无补救余地。

5.3 削减后极端场景丢失

同步回代消除法天然倾向于保留“平均形态”的场景，因为它优先删除的是“离其他场景最近”的场景，那些远离群体的极端场景反而更容易活下来。但如果你恰好需要评估极端情况，比如保供场景、大发场景，削减后极端场景的占比可能会明显降低，导致风险评估结果偏乐观。

我的经验是做两类场景集：一组用标准削减得到常规场景集，用于日常调度计算；另一组在削减前人为注入若干条构造极端场景，并给它们设定一个最小保留概率，再执行削减。这样可以得到一个专门用于极端情况分析的场景集。两个场景集分开使用，不要混在一起。

5.4 数据量不足时的替代方案

有些风电场运行历史不到一年，或者数据质量差到没法拟合出可靠的边际分布。这时候蒙特卡洛方法的基础就不牢了。一个可行的替代方案是采用分块Bootstrap直接重采样历史曲线片段，然后拼接生成新场景。比如把每条历史曲线切成长度为4小时的片段，随机抽取片段拼接成96点曲线。这种方式可以保留时序自相关结构，虽然引入的多样性有限，但在数据有限的情况下比盲目拟合分布稳健得多。

6. 从场景生成延伸开去的几点思考

做到后面你会发现，场景生成和削减这套方法论并不只服务于风电、光伏这类新能源不确定性建模，它本质上是一个通用的“概率信息压缩”工具。电网的负荷预测、电力市场价格模拟、极端天气事件建模，甚至是金融领域的资产收益模拟，核心逻辑都是一样的——从大量可能的状态中提取出少量有代表性的状态，保持概率分布特性，降低后续计算负担。

这几年还有一个明显的趋势，就是用深度学习来做场景生成。生成对抗网络（GAN）、变分自编码器（VAE）被越来越多地应用到时序场景生成中，它们的好处是能直接从历史数据里学习复杂的非线性时序相关性，不需要显式建模相关系数矩阵。但这类方法的通病是训练不稳定、可解释性差、生成样本可能违反物理约束。我在实际应用中更倾向于把深度生成模型和传统蒙特卡洛方法结合：用深度学习模型捕捉复杂的条件分布结构，用传统削减方法做后处理，保证最终场景集物理可行、概率保真。

另外说句题外话，很多人容易忽略场景削减和场景生成之间的迭代关系。如果你在做的项目对精度有很高要求，可以试试“生成-削减-重生成”的闭环迭代思路：先削减出一批代表性场景，用它们求解优化问题，然后在优化解附近的区域加密采样生成新场景，重新削减，更新代表场景。这个迭代过程在实际项目中往往两三轮就能明显提升最终优化方案的可靠性。

我个人在实际操作中的体会是，这套方法真正难的地方不在单个算法的实现，而在于怎么把每一步的参数和细节都调整到匹配真实数据的状态。同样的算法，换一个风电场、换一种数据粒度、换一个目标函数，最优点可能完全不一样。所以多留出时间做敏感性分析，多对比几组参数结果，比死磕某一个理论细节更有工程价值。

如果你正在做新能源出力场景生成、随机优化或者可靠性分析相关的工作，希望这篇梳理能帮你少走一些弯路。拿一份真实数据从头到尾跑一遍，把每个环节的中间结果都亲自看一眼，你会很快建立起对这个方法论的整体感觉。