做新能源出力场景的人都知道一句话：场景好不好，不只看分布像不像，还要看曲线“像不像真实发生的事”。我之前在处理某个光伏基地的出力场景生成任务时，第一版直接用蒙特卡洛对每个时段的分布独立抽样，生成了1000条日曲线。分布检验全过——每个时刻的均值、方差、分位数都对得上历史。可把任意一条曲线画出来，立刻露馅：上午十点还是满发，十一点突然掉到10%，下午两点又拉回满发。相邻时段之间完全没有“记忆”，而这种锯齿在真实光伏出力里几乎不会出现。这就是典型的“抛弃了时序相关性”的蒙特卡洛场景生成。
这篇内容就围绕“考虑时序相关性MC的场景生成与削减”这个主题展开，讲讲为什么时序相关性这么重要、怎么把相关性“装进”MC抽样，以及生成几千条场景之后如何用削减算法收缩成十几条有代表性的场景，供随机规划、生产模拟、容量评估使用。整理的都是我在实际项目中反复验证过的做法，适合做电力系统不确定性分析、新能源功率建模、随机规划基础场景库构建的同行参考。
1. 时序相关性为什么是场景生成里最容易被忽略的坑
1.1 独立抽样满足的是“分布”，不是“过程”
蒙特卡洛最直觉的用法，是把一天24个时段当成24个独立的随机变量，对每个时段单独拟合一个边缘分布，然后独立抽一条随机数，再逐时段映射回功率。问题是，这种做法的潜台词是“t时刻取什么值与t-1时刻毫无关系”，这在很多场景中都不成立。真实的风电出力是一条连续受惯性约束的物理过程，光伏出力受云层和地表辐照度的转移规律影响，相邻时段的功率变化率天然受限。
从数据上讲，我统计过某光伏基地历史日曲线的相邻时段Pearson相关系数，白天时段普遍在0.75到0.9之间，而独立抽样生成的场景集，同样指标只有0.2左右。分布形态完全一致，自相关结构却被抹平了。这种“分布正确、过程失真”的问题，恰恰是独立MC抽样最容易踩的坑。
1.2 过程失真之后，下游优化模型会被“骗”
场景最终是要喂给下游随机规划的，最常见的是机组组合、储能调度、容量配置。这些模型都包含爬坡约束或者电量平衡约束，因为物理上火电、储能、可调负荷的爬坡都不是无限的。如果场景里的时序关系被破坏，模型会看到一个“一会满发一会骤降”的世界，于是为应对高频大幅波动去预留更多向上/向下备用，或者反过来低估了连续几个时段低出力带来的电量缺口风险。
我做过一个小实验：同一组历史数据生成两个场景集，一个独立抽样，一个考虑时序相关性，分别代入一个模拟风电的随机调度模型。结果独立抽样那组多了约一倍的备用需求。原因是锯齿曲线把相邻时段的变化率拉大，模型的爬坡约束被极大地激活。可以这么理解：场景集是模型认知世界的“训练样本”，样本过程不对，模型学出来的策略自然不对。
1.3 什么时候可以不用管时序相关性
这不是所有场景生成都需要考虑时序。如果你的目标只是评估每个时段的边际分布，比如做概率潮流时想得到节点电压的概率分布，各时段独立抽样就够了。但只要你关心跨时段行为，比如储能充放电策略、备用容量、连续阴天/无风段的最大缺额，时序相关性就必须进入建模视野。
有没有一个“量化判断标准”？可以算一下相邻时段互信息的占比，或者直接画滞后1阶自相关函数的衰减曲线，对比独立抽样下自相关衰减速度，差异一目了然。以光伏为例，如果滞后1阶自相关系数低于0.4，对调度类应用的失真就比较明显了；如果把目标改成“连续两小时低出力频次”，独立抽样和带相关性抽样可以差出一个量级，这种差距不是靠增加抽样次数能补回来的。
2. 给MC抽样加上“时间记忆”：三种可落地的相关结构建模方案
2.1 路线A：马尔可夫链转移矩阵
最早想到的是马尔可夫链。把每个时段出力离散成若干个状态，比如把额定功率按0~0.1、0.1~0.3、0.3~0.6、0.6~1.0分成4个档，再用历史数据统计相邻时段的状态转移概率矩阵。生成场景时，从初始状态分布出发，按转移矩阵逐时段抽样整条路径。
这个方案的优点是直白，代码量很小，而且天然保留了一阶时序相关性。缺点是状态划分边界敏感，状态数太少会丢失出力精度，状态数太多转移矩阵又变得稀疏。在光伏场景里，云层遮挡带来的多状态转移有时远超一阶，马尔可夫链容易把“中间阴天、后面转晴”的长期依赖模拟成随机游走。
2.2 路线B：ARIMA模型加残差重抽样
第二种路线是直接对历史序列拟合ARIMA模型，核心思路是把时间序列分解成确定性的自回归结构加白噪声残差。具体做法：对每一日的历史出力曲线做季节差分或去趋势处理，估计AR(p)项的系数，然后从历史残差中重抽样，再通过模型递推合成新的完整序列。
这条路线的好处是保留了线性自相关结构，而且残差重抽样比纯正态假设灵活。问题在于ARIMA对非线性相关、异方差特征适配一般。光伏出力的方差在正午大、早晚小，如果不对方差做额外建模，容易在峰谷时段生成不合理的波动。如果手里是风电数据、时序线性特征比较强，ARIMA可以作为一个快速基线。
2.3 路线C：t-copula加多元Student-t抽样（我最推荐）
我最终保留在项目里的方案是t-copula。它的核心思想是“边缘分布与相关结构解耦”：先不管曲线整体形状，只用概率积分变换把每个时段的历史功率转换到[0,1]均匀空间，再用均匀分位数去估计时段间的相关结构；抽样时先生成带相关结构的多元随机数，再逆变换回功率空间。
这样做的好处非常明显：边缘分布可以任意复杂（哪怕逐时段用非参数分布拟合），而相关性由一张相关系数矩阵单独表达，两者互不污染。下面这段是我在项目里的核心抽样代码，刻意注释得很详细：
import numpy as np from scipy.stats import t as student_t # X_hist: 历史日曲线矩阵, shape=(n_days, T)，T为日内时段数 T = 24 # 1) 对每个时段t，用历史数据拟合边缘分布F_t（核密度或参数分布均可） # 2) 概率积分变换：把每列的功率值转为[0,1]均匀分位数 U = np.zeros_like(X_hist) for t in range(T): U[:, t] = fitted_cdf_list[t](X_hist[:, t]) # 3) 估计相关结构：把均匀分位数映射到t分布空间 nu = 5 Y = student_t.ppf(U, df=nu) # 每列是t_nu分布的分位数 Sigma = np.corrcoef(Y, rowvar=False) # 时段间的相关矩阵 # 4) 生成MC场景：先抽独立t分布随机数，再用Cholesky分解叠加相关结构 n_scen = 1000 Z = student_t.rvs(df=nu, size=(n_scen, T)) L = np.linalg.cholesky(Sigma) S = Z @ L.T # 行向量带上相关矩阵 U_scen = student_t.cdf(S, df=nu) # 回到均匀空间 X_scen = np.empty_like(U_scen) for t in range(T): X_scen[:, t] = fitted_ppf_list[t](U_scen[:, t]) # 逆变换回功率几个设定说一下：自由度nu取5而不是更大的数，是为了保留“尾部相关性”，让连续极端天气有合理出现概率。nu太大会退化成正态copula。相关矩阵用的是Spearman秩相关还是Pearson相关？对于光伏功率这类强非高斯变量，我建议用秩相关（对均匀分位数算Pearson相关等价于Spearman秩相关），避免被极端值带偏。
三条路线的取舍我也整理成了表格：
| 方案 | 实现成本 | 时序相关性保留 | 非线性/异方差适配 | 适用场景 |
|---|---|---|---|---|
| 马尔可夫链 | 低 | 一阶 | 一般 | 快速原型 |
| ARIMA+残差重抽样 | 中 | 线性多阶 | 弱 | 长期序列建模 |
| t-copula+多元t抽样 | 中高 | 秩相关全时段 | 强 | 工程主力方案 |
2.4 还有一个“笨办法”：时序块抽样
如果历史数据量非常充足，还可以直接做块抽样：从历史日曲线里随机抽一段连续长度为L的子块，拼接成新曲线。因为子块来自真实数据，它的非线性相关、异方差、尖峰特性全都被原样保留。
缺点是需要数据覆盖足够多的天气类型，否则容易反复抽到同几条模板曲线，多样性不够。在我这个项目里，历史有近3年日曲线，我会把块抽样生成的场景作为t-copula方案的“对照组”，两者分位数带和极值覆盖基本一致时，就认为场景库是可靠的。
3. 场景削减：把千条曲线收缩到十几条，本质是加权取舍
3.1 为什么非削减不可
MC场景生成很容易一把生成几千条曲线，这本身不是问题，问题是下游优化。以两阶段随机机组组合为例，每多一个场景，不代表仅仅多一个决策分支，而是整个目标函数里多一组带权重的约束和变量。一千个场景可以直接让一个原本几分钟能解完的混合整数规划，在数小时内都无法收敛。
所以工程上几乎必然要做一个场景削减步骤：用少量典型曲线带权重地逼近原始场景集合。削减的目标不是“挑长得好看的曲线”，而是最小化削减后的概率分布与原始分布之间的某种距离。工程里最常用的是Kantorovich概率距离，它本质上是带权重的Wasserstein距离在多维欧氏空间中的近似。简单理解：在原始场景集里找到一些“代表点”，把周边场景的概率归并给代表点，让整体概率分布尽量不塌缩。
3.2 快速前向选择算法：从零开始挑典型场景
我项目里最常用的是快速前向选择。它从空集开始，第一个场景选“离所有场景加权距离之和最小”的那个，即最居中的曲线；随后每一步迭代，从剩余场景中选一个“使得整体概率距离下降幅度最大”的场景加入集合。等价地说，每次挑一个能把最没有代表的区域重新覆盖住的场景。
下面是核心骨架代码，为了可读性我没做任何复杂度优化，思路很清楚：
import numpy as np from scipy.spatial.distance import cdist def forward_selection(scenarios, probs, n_select): n = scenarios.shape[0] D = cdist(scenarios, scenarios, metric='euclidean') selected = [] remain = set(range(n)) # 选第一个：加权距离和最小的那个 total_dist = np.sum(probs[:, None] * D, axis=0) first = int(np.argmin(total_dist)) selected.append(first) remain.remove(first) for _ in range(n_select - 1): best_gain = -1.0 best_cand = None for cand in remain: # 计算每个未选场景在加入cand后，距离已选集合的“下降增益” gain = 0.0 for i in remain: if i == cand: continue d_before = min(D[i, k] for k in selected) d_after = min(D[i, k] for k in selected + [cand]) gain += probs[i] * (d_before - d_after) if gain > best_gain: best_gain = gain best_cand = cand selected.append(best_cand) remain.remove(best_cand) # 归并权重：把每个未选场景的概率加到离它最近的已选场景上 weights = np.zeros(n) for i in range(n): if i in selected: weights[i] += probs[i] else: nbr = min(selected, key=lambda k: D[i, k]) weights[nbr] += probs[i] return selected, weights注意遍历里有一个for i in remain嵌套，n=1000时大约是百万量级的min操作，计算机能跑但偏慢，实际工程我会用numpy矩阵化重写，上面版本是为了讲清楚过程。
3.3 同步回代削减：从满集出发反向删场景
与快速前向相对的是同步回代，即从完整的场景集出发，每一轮找出“删除后导致概率距离增加最小”的场景，把它删掉，同时把它被归并的概率加到离它最近的剩余场景上。前向选择的优点是最终场景之间相距较远，“居中代表”明显；回代削减的优点是计算过程中一直保留原始场景的全局分布信息。
我个人的偏好是：如果目标场景数量与原始数量相差一个数量级以上（比如1000削到10），前向选择往往表现更稳；如果只是从100削到50这种“轻削减”，同步回代会保留更多细节。两个方法的时间复杂度都比较高，但都有成熟的高效实现，动态规划版的快速前向在实际工程里很常用。
| 方法 | 思路 | 适用削减比例 | 优点 | 缺点 |
|---|---|---|---|---|
| 快速前向选择 | 从空集挑代表 | 大比例（10~20倍以上） | 代表性强 | 初始迭代慢 |
| 同步回代削减 | 从满集删除冗余 | 小比例（2~5倍） | 保留全局分布细节 | 大比例时计算量大 |
3.4 削减到多少个场景合适：找概率距离曲线的“拐点”
削减数量不是一个拍脑袋参数。我会固定原始MC场景集，然后依次削减到5、10、20、50、100个，分别计算削减后的场景集与原始场景集之间的概率距离。画出来的曲线通常在10到30之间出现明显拐点：30个场景以后，继续增加场景数量，概率距离下降已经微乎其微。
我的经验是：对于日内24时段的风光功率场景，500到1000条原始场景削到10到30条，概率距离可以控制在15%以内；如果削到5条，概率距离会跳到30%以上，这时候除非你对计算速度极度敏感，否则不建议。
提示：概率距离曲线不是单调下降的，因为每个场景被选中后还会带走周边概率，个别迭代步可能出现距离短暂回升。要判断拐点，看整体趋势而不是单点波动。
4. 削减效果验证：不能只盯一张图，要看下游闭环
4.1 三个局部指标：概率距离、统计矩、自相关结构
削得好不好，我习惯从三个层面看。第一是概率距离本身，这已经由削减算法优化过了，但要看绝对值有没有失控。第二是逐时段统计矩，把原始MC场景集和削减场景集分别算均值、标准差、P10和P90分位数，画在同一张图上。如果削减后的分位数带宽明显变窄，说明场景多样性被砍掉了。第三是时间结构，很多人在这个环节失守——削减后虽然每个时刻分布都对，但场景数量太少，时段间的自相关结构会被拉平，这是削减方法用欧氏距离作为相似度度量时很常见的副作用。
我在项目里的实际对照数据大致如下（只列几个代表指标）：
| 指标 | 原始MC场景（1000条） | 削减场景（15条） |
|---|---|---|
| 平均时段均值误差 | 基准 | 0.8%以内 |
| P90分位数偏移 | 基准 | 4.2% |
| 相邻时段平均自相关 | 0.82 | 0.71 |
| 极端低出力场景覆盖 | 100% | 93% |
自相关从0.82降到0.71，绝对值看起来还行，但如果你拿着这15条场景去做储能调度，低谷段的连续时长会被低估。所以我在削减后还会做一次“自相关校验”，如果退化严重，就改用带一阶差分距离加权的相似度矩阵，或者直接额外保几条“过程响应”的极端场景。
4.2 最关键的一步：把场景喂进下游模型做闭环测试
所有分布指标都是间接证据，真正要看的，是场景削减前后的策略质量差异。我在模拟项目中做过一次完整闭环：先把场景集A（1000条）代入随机机组组合模型，求出的最优期望成本作为“全场景基准”；再把削减后的15条场景代入同一个模型，比较期望成本和第一阶段的决策变量差异。
结果是期望成本偏差在1.5%左右，第一阶段机组启停方案的重合度在90%以上。作为对照，如果从原始场景里“随手挑15条”而不是用削减算法，期望成本偏差会跳到8%到10%，启停方案也频频变化。这组实验强烈说明：削减算法不是可有可无的锦上添花，而是保证随机规划结果稳定的必要步骤。
我还会做一个更苛刻的鲁棒性测试：用削减场景训练出来的调度方案，拿到另外独立生成的1000条“未见场景”上回测，看是否频繁出现切负荷或弃风。在我的项目里，约束违反率控制在3%以下，我才敢把场景库交给下游。
4.3 可视化的呈现方式
分享一个比较实用的展示方法：削减后的15条场景不要全画实线，按概率权重给透明度，画成“瀑布图”，一眼能看出哪些是主线、哪些是小概率的尾部场景。再叠加上原始MC场景的P10-P90阴影带，能直观看到削减场景是不是恰好把分位数带的界覆盖住了。
这种方法在写报告、给非算法背景的同事评审时特别管用，比贴一堆数值表有效得多。尤其是场景库要移交到别的团队使用时，一张“权重透明度图加分位数带图”就能把代表性讲清楚。
5. 踩坑记录与参数建议
5.1 坑一：削减算法输出的权重被下游丢掉
这是我最常见到的低级错误，严格说不是算法的坑，而是接口的坑。很多开源工具削减完只导出“场景曲线列表”，权重单独存一个数组；下游优化程序如果只读曲线没读权重，就等于默认15条场景等概率。这会彻底毁掉削减的意义——为极小概率的尾部场景分配了过高权重，导致结果整体偏移。
建议在交付场景库时，把场景和权重打包成结构化文件，并在接口测试里专门做一个“权重和为1”的断言。如果用的还是老式CSV交接，建议加一列权重，并在文档里用小例子说明：削减前100条场景里只有2条极端大风场景，权重加起来不超过0.05；如果不读权重，这两条极端场景就占了两份的“话语权”，结果会明显失真。
5.2 坑二：t-copula的自由度选太大
自由度nu是t-copula的敏感旋钮。nu趋近于无穷时，t分布退化为正态分布，尾部相关性完全丢失。对风光这种经常出现连续极端天气的过程，低估尾部相关性会让场景库漏掉最危险的组合。
我测试过nu=50和nu=5两组场景，后者的P1分位数的极值深度明显更深。建议从nu=3到10之间遍历测试，用历史极端事件的发生频率来反推，哪个值能让场景库的极值频率对得上，就用哪个。不要因为“t分布在线性相关结构上没有解析解”这类理由回避它，工程里更关心的是极值能不能被覆盖，而不是解析优雅。
5.3 坑三：每时段单独拟合分布时过拟合
逐时段拟合边缘分布听起来很合理，实际上很容易被小样本坑到。一年365天、一天24点，白天可用时段可能只有几千个样本，用核密度估计去拟合会出现明显的“尖峰毛刺”，体现在场景里就是某些时段的极值被反复抽中。
我在项目里的处理是：先把历史时段做聚类，比如把日升、日落、平谷段分成3到4类，在同一类时段内共享分布形状，再乘一个逐时段的尺度和偏移量。这样既保留逐时段差异，又避免分布形状被噪声支配。样本量少于1000天时，这种聚类共享分布的做法几乎是必须的。
5.4 坑四：相似度度量只用欧氏距离
经典削减算法几乎都默认欧氏距离，但它对“时间轴上的小错位”不敏感。两条光伏曲线，一条在12:00到达峰值、一条在12:15到达峰值，本质上很接近，但欧氏距离会判成完全不同。遇到这种情况，可以用动态时间规整（DTW）距离替代欧氏距离，或者对一阶差分后的序列做削减。
后者有一个额外好处：直接维护相邻时段变化率的分布，能有效缓解削减后自相关结构退化。我自己常做的一个变通是，在欧氏距离基础上加10%到20%的一阶差分距离权重，既能保持大部分熟悉的行为，又能起到“过程相似性”的引导作用。
5.5 一套我实测靠谱的参数基线
整理一套可以直接当启动配置的参数，方便快速复现：
| 参数 | 建议值 | 备注 |
|---|---|---|
| MC场景数 | 500~2000 | 视计算资源，越多越稳 |
| t-copula自由度nu | 3~8 | 用极值频率反推 |
| 边缘分布拟合 | 非参数核密度 | 数据量不足时用时序聚类共享 |
| 削减方法 | 快速前向选择 | 大比例削减优选 |
| 相似度度量 | 欧氏距离 | 自相关退化时改用差分距离 |
| 削减目标数量 | 10~30 | 看概率距离拐点 |
| 权重验证 | 权重和=1 | 接口集成时强制校验 |
最后分享两个落地经验。第一，我在实际工程版本里不会“先削减、后使用”就完事，而是先生成4000条MC场景，削减出20条典型场景形成基础场景库，再根据项目关注的尾部风险，手工追加5到8条“极端场景”并调整权重。比如连续48小时低出力、正午前后剧烈先降后升、设备降额导致的出力折减。这样既保住了随机规划的计算速度，又不会丢掉长尾风险，比单纯加大场景数量划算得多。
第二，场景库交付后一定要留“版本号”——分布参数、相关性矩阵、削减随机种子都记录下来，否则下游模型每一次跑出不同结果时，你很难判断是场景变了还是调度模型本身的问题。这套方法论在我参与的模拟光伏与储能容量配置项目里实测下来很稳，也算是我做场景生成与削减这几年最值钱的一套经验。
