几年前我第一次做新能源出力场景的时候，是被一句话敲醒的：你采一万条曲线丢进模型，整个调度计算根本跑不完。那时候我还在老老实实把蒙特卡洛（MC）采样出来的一万多条时序曲线全部塞给机组组合模型，结果一个晚上加一个白天都没算完。后来才慢慢想明白，场景生成和场景削减从来不是两件可以分开凑合的事——生成的时候就要考虑时序相关性，削减的时候更要保住这些相关性。这篇文章就把我在这个方向上的完整思路、实操流程、踩过的坑，以及沉淀下来的这套方法，从头到尾拆开聊一遍。

适合看这篇内容的人，主要是做新能源出力不确定性建模、随机规划、可靠性评估的工程师和研究生。如果你已经在用蒙特卡洛采样，但发现生成的风电、光伏曲线怎么看都觉得不对劲，或者场景削减后风险指标漂移得厉害，那这篇文章应该能帮你把整条链路理顺。

1. 先想清楚：时序相关的场景生成到底在解决什么问题

1.1 场景在不确定性决策里的角色

先退一步看场景生成的意义。比如你要做明天的火电开机计划，风电出力不确定，那总不能让机组组合模型只靠一个期望出力曲线做决策，否则半夜无风、晚高峰缺电的情况会被完全平均掉。更合理的做法是：准备一组"可能发生的未来"，让优化模型在每一种可能下都能找到一个不至于太差的开机方案。这一组"可能发生的未来"，就是场景。

场景是一组带概率的多维时序向量。在风电、光伏出力场景里，每个场景就是一条24点（或96点）的出力曲线，所有场景合起来，近似的不是某一条确定曲线，而是一个概率分布。蒙特卡洛方法在这里扮演的角色，就是从这个概率分布里反复抽样，用大量样本去逼近真实分布，再把逼近出来的样本交给下游的优化或评估模型。

这个思路本身没有错，问题出在"从哪个概率分布里抽样"。如果只是对每个时刻分别拟合一个边缘分布、然后独立采样，那采出来的每个时刻的边际分布是对的，但时刻与时刻之间的关系完全失控。

1.2 独立采样为什么在工程里翻车

举个具体的例子。光伏出力在晴天的典型形态是中午高、早晚低，10点和11点的出力具有很强的正相关性，不可能出现10点还有800MW、11点直接掉到50MW这种突变。风速出力也是同样道理，相邻小时的风速变化受大气惯性约束，不会像白噪声一样乱跳。

如果你对每个时刻独立采样，这种"物理上不可能"的曲线一定会批量出现。把这些伪曲线当成真实场景输入给优化模型，模型会为了应对这些不存在的剧烈波动去做过度保守的决策，或者在某些极端组合下得出完全错误的评价结果。

时序相关性不是锦上添花，而是物理约束在统计层面的体现。这也是"考虑时序相关性MC场景生成"这个方向存在的根本原因：蒙特卡洛负责提供随机性，但随机性必须受相关结构约束。

1.3 整条链路的分工：从历史数据到可用场景

在实际项目里，我一般把整条流程拆成六个环节：

1. 历史数据清洗和归一化，解决缺失值、限电数据、异常毛刺；
2. 边际分布建模，对每个时刻单独拟合分布；
3. 相关结构分析，量化时刻之间的依赖关系；
4. 蒙特卡洛采样，生成大量原始场景；
5. 场景削减，把几千条场景压到十几条；
6. 双层校验，从分布和时序两个角度确认结果可用。

很多初学者一上来就研究各种高级生成模型，却把数据清洗和相关矩阵估计这种基础环节糊弄过去，后面一定会出问题。这篇内容会按这个链路逐步展开，每一步都能直接落地。

2. 时序相关性建模：从独立样本到"会呼吸"的曲线

2.1 先分清要控制哪种"相关"

做时序相关性建模，第一步要分清"相关"有几种形态。一种是无功功率空间上的相关，比如相邻两个风电场的出力同步性；另一种是时间上的相关，比如同一个风电场当前时刻和下一时刻出力的依赖。本文主要讨论时间维度上的相关，也就是一条曲线内部各时刻之间的相关结构。

在实际工作中，很多做多元不确定性的人会把空间相关和时间相关混在一起处理。严格一点的做法是：先用空间相关生成多场站的联合场景，再在每个场站的时序维度上施加时间相关。也就是说，你可以把"场景矩阵"看成多个风电场乘以多个时刻的二维结构，空间和时间两套相关矩阵可以分开估计、分层采样。先处理时间相关，是搭建整条链路的基础。

2.2 高斯Copula + Cholesky分解：性价比最高的方案

要把时序相关引入蒙特卡洛采样，我最常用的方法是"高斯Copula + Cholesky分解"。这个组合的好处很多：能灵活对接任意边际分布，计算量可控，而且实现起来只需要基础线性代数和逆变换采样。

思路分三步。第一步，从历史数据里提取一个"正态空间"下的相关矩阵；第二步，生成一组带这个相关结构的标准正态样本；第三步，通过等概率变换，把正态样本映射到各时刻的真实边际分布上。

第二步的核心是Cholesky分解。假如目标相关矩阵是R，Cholesky分解把它拆成R=LL^T；然后生成一个N行T列的独立标准正态矩阵Z，再算Y=ZL^T，得到的Y就是一个N行T列、列与列之间相关矩阵为R的样本。整个过程非常快，几千个场景也就一瞬间。

代码如下：

import numpy as np # R: 目标相关矩阵 (T x T) L = np.linalg.cholesky(R) # Z: 独立标准正态样本, shape (N_samples, T) Z = np.random.randn(N_samples, T) # 列之间带有目标相关结构的正态样本 Y = Z @ L.T

第三步是把Y的每一列，通过正态累积分布函数的逆变换，映射到对应时刻的目标边际分布。假设第t时刻的边际分布是F_t，那么该时刻的最终样本就是F_t的逆函数作用于Φ(Y_t)。这里Φ是标准正态的CDF。这样做能保证每个时刻的边际分布完全等于F_t，同时时刻之间的秩相关结构基本保留下来。

2.3 相关矩阵的估计：秩相关替换Pearson相关

相关矩阵这一步是很多人翻车的重灾区。如果直接拿原始风电出力数据算Pearson相关矩阵，因为出力分布通常是非正态、重尾的，算出来的相关系数很容易被极值拉偏。

更好的做法是先做"正态分数变换"：把每一时刻的历史出力值按大小排序，得到排名，再把排名映射成标准正态分位数，然后再计算Pearson相关矩阵。这个流程的实际意义是，把所有边缘分布都"规范"到同一个正态尺度下，互相之间才有可比性。

具体代码：

from scipy.stats import norm, rankdata # X_hist: (N_days, T) 历史日出力曲线 n = X_hist.shape[0] rank_hist = np.apply_along_axis(rankdata, 0, X_hist) z_hist = norm.ppf(rank_hist / (n + 1)) # 正态空间的相关矩阵 R = np.corrcoef(z_hist, rowvar=False)

除了转换为正态分数再算Pearson，直接用Spearman秩相关矩阵也是一个不错的选项。秩相关的好处是对异常值不敏感，而且在高斯Copula的框架下，Spearman相关和正态空间Pearson相关之间存在简单转换关系。实际经验是：如果数据量只有几十天到一百天，Spearman更稳；如果数据量大到几百天以上，两者差异会缩小。

2.4 边际分布选择：参数分布还是核密度估计

生成场景时，各时刻的边际分布选择直接决定最终曲线的形态。风电出力速度分布常用Weibull或截断正态，光伏出力常用Beta分布，但这些参数分布都有个问题：真实数据往往是双峰、偏斜、带平台段的，单一参数分布很难拟合到位。

我在项目中更倾向于用核密度估计（KDE），每个时刻单独估计。KDE不需要过多假设分布形态，数据长什么样就拟合什么样。需要注意的是KDE的两个坑：一是带宽选择不合理会过光滑或过毛糙，一般用Silverman规则作为起点；二是边界问题，比如光伏出力的边界是0到装机容量，直接用普通KDE会把概率泄漏到负值或超过上限的区域，需要用反射边界或对估计结果做截断处理。

如果历史数据特别少，比如只有30天，KDE会拟合出很多假峰。这种情况下参数分布反而更稳。我的经验是：样本少于100天，优先用参数分布；样本大于100天，优先用KDE。

2.5 为什么不直接用ARMA或马尔可夫链

很多人会问，时序相关直接用ARMA模型或者马尔可夫链不就行了？ARMA模型对边际分布的正态假设很强，虽然可以扩展成ARMA-GARCH，但要同时保证非正态边际分布和指定相关结构，做起来非常绕。马尔可夫链的思路是把出力离散成有限状态，再估计转移概率矩阵，但离散化本身会损失连续出力的精度，尤其是出力接近0和接近满发这两个区间，状态划分稍粗一点就会扭曲尾部分布。

Copula方法之所以更适合这类场景，是因为它把"边际分布"和"相关结构"分开建模：边际分布随便选，参数分布或KDE都行；相关结构统一塞进一个Copula函数里。灵活性高很多，蒙特卡洛采样也顺理成章。

3. 场景削减：两千条曲线如何留下十张王牌

3.1 场景数的代价：不止是算力

假设你生成3000个场景，下游是一个含二进制变量的机组组合模型，场景数每多10倍，模型规模就膨胀好几倍，求解时间可能从分钟级跳到小时级甚至不收敛。场景削减的核心目标，就是用少量带权重的场景去逼近原始场景集的概率分布，让下游决策模型算得动，同时决策质量不明显下降。

这里要注意，削减后的每个场景不再等概率，而是带一个权重。最后你手里的是"10条曲线+10个概率"，而不再是"10条等概率曲线"。这个权重的准确性直接关系到优化结果的期望成本计算。

3.2 同步回代削减：被删场景的概率要"交接"

最常用的场景削减算法是同步回代削减（Simultaneous Backward Reduction，SBR）。直觉很简单：每一轮从场景集中删掉一条"对整体分布影响最小"的场景，然后把被删场景的概率加到离它最近的保留场景上，重复直到剩下指定数量。

删掉一条场景的代价怎么算？一般用概率乘以它到最近场景的距离来近似。距离通常取欧氏距离，概率乘以距离表示"把它并到邻近场景后带来的分布扭曲程度"。每一轮选代价最小的场景删除。

算法逻辑如下：

from scipy.spatial.distance import cdist def sbr(S, probs, K): # S: (N, T) 场景矩阵, probs: (N,), K: 目标场景数 while len(S) > K: D = cdist(S, S) np.fill_diagonal(D, np.inf) # 排除自己 costs = [] for i in range(len(S)): j = int(np.argmin(D[i])) costs.append(probs[i] * D[i, j]) idx = int(np.argmin(costs)) # 要删除的场景 j = int(np.argmin(D[idx])) # 最近的保留场景 probs[j] += probs[idx] S = np.delete(S, idx, axis=0) probs = np.delete(probs, idx) return S, probs

SBR的理论背景是Wasserstein距离的启发式求解。每一步删除都是在尽量减小原有场景集与剩余场景集之间的Wasserstein距离。这个算法效果稳定，也是我默认的削减去。

3.3 快速前向选择与K-means的取舍

除了SBR，还有两类常见方案。快速前向选择（FBS）是反向思路：从空集开始，每一步从候选集中选入一个场景，使得已选场景集对原始分布的近似程度最高，直到达到K个。FBS对极端场景的保留能力强于SBR，但实现更复杂，每一步都要评估选入候选者对分布的影响，计算量更大。

K-means聚类是最简单粗暴的：把所有场景当成高维向量聚类，用簇中心代表一个场景，簇大小占比作为概率。优点是速度快、实现容易，缺点也很明显——簇中心是平均形态，会把一些尾部特征抹平。如果场景分布里有明显的"极端大风日"和"平稳小风日"，聚类中心很可能介于两者之间，两头都不像。

我的选型建议是：常规场景削减用SBR；业务特别关注极端场景和尾部风险时，用FBS或对极端场景进行分组强制保留；数据规模非常大且对精度要求不高时，再考虑K-means。

3.4 距离度量和削减数量：两个容易被忽略的参数

距离度量是削减算法里最容易想当然的地方。对24点等间隔时序曲线，直接用欧氏距离就是合理选择。有人喜欢用动态时间规整（DTW）来衡量曲线形状差异，但DTW允许时间轴拉伸，会把"9点最大出力"和"10点最大出力"视作相似——在风电出力场景里这种时间错位没有物理意义，反而会让削减结果失真。所以除非场景是不同长度的事件序列，否则不建议在等间隔时序上用DTW。

削减数量怎么定？经验上从2000到5000个原始场景削减到10到20个是比较常见的范围。但具体取多少，不能拍脑袋。我一般会画一条敏感性曲线：分别削减到5、10、15、20、30个场景，计算削减后场景集的均值、方差、分位数与原始场景集的偏差，看偏差在哪个数量级开始平缓。如果10个场景已经让均值误差小于1%、分位数误差小于2%，那就是够用的。

4. 怎么证明生成和削减的结果真的能用

4.1 分布层面对比：看形状、看分位点

场景生成和削减做完了，不能只看曲线好看，要用指标证明"削减后的少数场景仍然代表原来那个分布"。

分布层面最直接的对比是均值、标准差、0.05/0.5/0.95分位数。我曾经做过一个某区域电网的风电场景项目，初始生成3000个场景，削减到10个后，各时刻均值最大偏差约1.2%，标准差最大偏差约3.6%，0.05分位数偏差约5%，整体可以接受。但如果0.05分位数偏差达到10%以上，我通常会怀疑削减数量太少或者距离度量选错了。

更严谨一点的指标可以用Wasserstein距离和KS检验。Wasserstein距离在分布对比里比KL散度更有意义，因为它同时考虑概率质量的空间位置，对场景削减这种"删掉某些场景"的操作更敏感。KS检验适合单时刻的分布对比，逐时刻做一遍，能找出哪些时段分布被削减扭曲得最严重。

4.2 时序相关层面对比：ACF和相关矩阵

分布检验只能证明单时刻的形态没跑偏，还要证明时序相关性没丢。两个常用工具：自相关函数（ACF）和相关矩阵。

ACF可以看滞后1小时、2小时、24小时的相关系数。把历史数据的ACF作为基准线，生成大量场景后画出ACF的分布包络，看历史ACF是否落在包络内部。如果生成模型根本不含时序相关，ACF会在滞后1阶就掉到接近0，一眼就能看出来。

相关矩阵对比更精确定量：先算历史出力在正态空间下的相关矩阵R_hist，再对生成场景也做同样变换、求相关矩阵R_syn。比较两者对应元素的最大绝对误差。如果最大误差小于0.05到0.1，说明时序相关保真度不错。这个指标我每次必做。

削减后的场景也要做同样的检查。因为SBR本质上是"删曲线"，它完全有可能会削减掉一些相关性信息。如果一个削减后场景集的相关矩阵和历史相关矩阵差异突然变大，那就说明削减数量太少，或者距离矩阵没有很好体现时序相关结构。

4.3 端到端验证：把场景塞进决策模型里看误差

统计指标都合格之后，最有说服力的验证是把场景塞进一个真实的下游决策模型里。以机组组合为例，做法是：

1. 用削减后的场景集求解，得到一阶段决策变量（比如机组开停机计划）；
2. 把这个计划固定下来，代入原始的3000个全场景，计算期望运行成本；
3. 同时用全场景直接求解（如果能算得动的话），得到理论最优成本；
4. 对比两步结果的偏差。

如果削减后得到的期望成本比全场景直接求得的成本偏差在1%到2%以内，说明场景削减对决策质量几乎没有伤害。这种做法在论文里常见，在工程项目里也同样适用。我这里强调一下：验证环节不能省，因为有时候统计指标都好看，但恰恰丢了某个关键约束触发的场景，决策结果就会截然不同。

4.4 结果可视化与汇报经验

做研究和做工程有一点是相通的：图比表格有说服力。我常用的三种可视化方式：

相关矩阵热力图，历史样本和生成样本各一张，并排放在一起，色彩差异一目了然。这是向同事证明"时序相关性被保住了"最直观的图。

随机场景带状图，把生成的上百条曲线画成透明细线，形成一个置信带，再把历史典型日的曲线叠加进去，看历史曲线是否落在带内。

削减结果图，把削减后的10条场景画在同一个坐标轴上，线宽或透明度按权重区分，旁边标注每个场景的概率。这张图直接交给调度运行人员看，他们能快速分辨哪些是高风险场景、哪些是常规场景。

5. 实战踩坑与调参心得

5.1 相关矩阵不正定：常见但别慌

Cholesky分解要求相关矩阵是正定阵。实际估计出来的R经常不是正定阵，尤其是当样本量少、或者多个时刻高度共线的时候，特征值可能出现负的小值，导致np.linalg.cholesky直接报错。

这不是实现错了，而是估计噪声。解决办法我一般用两种：

一是对角加载，把R替换成(1-ε)R+εI，ε取0.01到0.05，往对角线方向拉一把。二是特征值裁剪，把负特征值置为0再重新组合成协方差矩阵，最后归一化变成相关矩阵。更讲究一点可以用Higham最近相关矩阵算法。实际项目里对角加载最简单，而且对生成结果影响很小。

5.2 逆变换后相关性缩水：用秩相关和迭代校准

高斯Copula里，中间层的正态样本相关矩阵是可控的，但经过最后一步对边际分布的逆变换后，由于非线性映射，实际的Pearson相关矩阵会发生变化。这个现象在边际分布严重非正态时特别明显，比如风电出力大量堆积在0附近的时间段。

如果你一开始就是为了控制Pearson相关，那结果多半会低于目标值。解决办法是：目标相关矩阵改用秩相关（Spearman），因为秩相关在单调变换下保持不变，高斯Copula流程能基本保护秩相关。如果对线性相关的数值有硬性要求，可以做迭代校准：生成一批样本，算出实际秩相关，修改目标矩阵再重新生成，重复几次直到实际相关逼近目标。

5.3 削减丢掉了极端场景：风险视角的修正

SBR削掉的是"对整体分布影响最小"的场景，而影响最小经常意味着风险最温和。如果你下游做的是可靠性评估或者风险评估，那恰恰是最关心尾部风险的场景，比如持续无风、极端高温引发光伏满发但负荷也飙升。结果SBR可能把最关键的1%风险场景削掉了，期望成本计算出来虚低。

对策有三种：一是削减前把历史极端场景单独挑出来，强制加入保留集合，再对剩余场景做削减；二是采用分组削减，比如按季节或按天气聚类，每组独立削减，保证每个风险类别里都保留场景；三是改用FBS这类"从零选入"的算法，它对极端场景的覆盖能力更好。

我通常在项目里把三种办法都跑一遍，再对比尾部指标，选择最稳定的一种。

5.4 计算加速与工程落地要点

场景生成阶段可以并行化，因为蒙特卡洛采样本质上是可并行的计算。几千个场景，用Numpy向量化其实已经很快，根本不需要上多进程。真正耗时的反而是SBR里的距离矩阵计算，因为在不断删除场景的过程中，距离矩阵需要反复更新。

几个工程落地小建议：

1. 用float32存储场景矩阵，3000乘24的矩阵也能从72MB降到36MB，对几十个风电场的场景集来说意义明显；
2. SBR里用numpy的cdist函数计算距离，不要用纯Python循环；
3. 场景削减是一次性的离线计算，结果可以缓存下来反复用于不同下游模型，不必每次重跑；
4. 如果场景数到上万级别，距离矩阵内存会吃不消，这时考虑KDTree加速找最近邻，或者分块处理距离计算。

最后再分享一点个人体会。整套流程做下来，我最大的感受是场景生成和削减从来不是"两个步骤"，而是一整套"可控的近似"。生成阶段在近似真实概率分布，削减阶段在近似生成出来的大量样本，每一步都在牺牲一些细节，关键是得知道自己在牺牲什么、牺牲得能不能接受。如果让我再来一次，我会先把"生成3000个、削减到10个"这个端到端流程在小规模数据上完全跑通，把所有检验指标建立起来，再考虑t-Copula、DTW、深度生成模型这些更进阶的道具。基础链路稳了，后面的一切才有可能谈。