几年前我第一次做新能源出力场景的时候是被一句话敲醒的你采一万条曲线丢进模型整个调度计算根本跑不完。那时候我还在老老实实把蒙特卡洛MC采样出来的一万多条时序曲线全部塞给机组组合模型结果一个晚上加一个白天都没算完。后来才慢慢想明白场景生成和场景削减从来不是两件可以分开凑合的事——生成的时候就要考虑时序相关性削减的时候更要保住这些相关性。这篇文章就把我在这个方向上的完整思路、实操流程、踩过的坑以及沉淀下来的这套方法从头到尾拆开聊一遍。适合看这篇内容的人主要是做新能源出力不确定性建模、随机规划、可靠性评估的工程师和研究生。如果你已经在用蒙特卡洛采样但发现生成的风电、光伏曲线怎么看都觉得不对劲或者场景削减后风险指标漂移得厉害那这篇文章应该能帮你把整条链路理顺。1. 先想清楚时序相关的场景生成到底在解决什么问题1.1 场景在不确定性决策里的角色先退一步看场景生成的意义。比如你要做明天的火电开机计划风电出力不确定那总不能让机组组合模型只靠一个期望出力曲线做决策否则半夜无风、晚高峰缺电的情况会被完全平均掉。更合理的做法是准备一组可能发生的未来让优化模型在每一种可能下都能找到一个不至于太差的开机方案。这一组可能发生的未来就是场景。场景是一组带概率的多维时序向量。在风电、光伏出力场景里每个场景就是一条24点或96点的出力曲线所有场景合起来近似的不是某一条确定曲线而是一个概率分布。蒙特卡洛方法在这里扮演的角色就是从这个概率分布里反复抽样用大量样本去逼近真实分布再把逼近出来的样本交给下游的优化或评估模型。这个思路本身没有错问题出在从哪个概率分布里抽样。如果只是对每个时刻分别拟合一个边缘分布、然后独立采样那采出来的每个时刻的边际分布是对的但时刻与时刻之间的关系完全失控。1.2 独立采样为什么在工程里翻车举个具体的例子。光伏出力在晴天的典型形态是中午高、早晚低10点和11点的出力具有很强的正相关性不可能出现10点还有800MW、11点直接掉到50MW这种突变。风速出力也是同样道理相邻小时的风速变化受大气惯性约束不会像白噪声一样乱跳。如果你对每个时刻独立采样这种物理上不可能的曲线一定会批量出现。把这些伪曲线当成真实场景输入给优化模型模型会为了应对这些不存在的剧烈波动去做过度保守的决策或者在某些极端组合下得出完全错误的评价结果。时序相关性不是锦上添花而是物理约束在统计层面的体现。这也是考虑时序相关性MC场景生成这个方向存在的根本原因蒙特卡洛负责提供随机性但随机性必须受相关结构约束。1.3 整条链路的分工从历史数据到可用场景在实际项目里我一般把整条流程拆成六个环节历史数据清洗和归一化解决缺失值、限电数据、异常毛刺边际分布建模对每个时刻单独拟合分布相关结构分析量化时刻之间的依赖关系蒙特卡洛采样生成大量原始场景场景削减把几千条场景压到十几条双层校验从分布和时序两个角度确认结果可用。很多初学者一上来就研究各种高级生成模型却把数据清洗和相关矩阵估计这种基础环节糊弄过去后面一定会出问题。这篇内容会按这个链路逐步展开每一步都能直接落地。2. 时序相关性建模从独立样本到会呼吸的曲线2.1 先分清要控制哪种相关做时序相关性建模第一步要分清相关有几种形态。一种是无功功率空间上的相关比如相邻两个风电场的出力同步性另一种是时间上的相关比如同一个风电场当前时刻和下一时刻出力的依赖。本文主要讨论时间维度上的相关也就是一条曲线内部各时刻之间的相关结构。在实际工作中很多做多元不确定性的人会把空间相关和时间相关混在一起处理。严格一点的做法是先用空间相关生成多场站的联合场景再在每个场站的时序维度上施加时间相关。也就是说你可以把场景矩阵看成多个风电场乘以多个时刻的二维结构空间和时间两套相关矩阵可以分开估计、分层采样。先处理时间相关是搭建整条链路的基础。2.2 高斯Copula Cholesky分解性价比最高的方案要把时序相关引入蒙特卡洛采样我最常用的方法是高斯Copula Cholesky分解。这个组合的好处很多能灵活对接任意边际分布计算量可控而且实现起来只需要基础线性代数和逆变换采样。思路分三步。第一步从历史数据里提取一个正态空间下的相关矩阵第二步生成一组带这个相关结构的标准正态样本第三步通过等概率变换把正态样本映射到各时刻的真实边际分布上。第二步的核心是Cholesky分解。假如目标相关矩阵是RCholesky分解把它拆成RLL^T然后生成一个N行T列的独立标准正态矩阵Z再算YZL^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, rowvarFalse)除了转换为正态分数再算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 ReductionSBR。直觉很简单每一轮从场景集中删掉一条对整体分布影响最小的场景然后把被删场景的概率加到离它最近的保留场景上重复直到剩下指定数量。删掉一条场景的代价怎么算一般用概率乘以它到最近场景的距离来近似。距离通常取欧氏距离概率乘以距离表示把它并到邻近场景后带来的分布扭曲程度。每一轮选代价最小的场景删除。算法逻辑如下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, axis0) probs np.delete(probs, idx) return S, probsSBR的理论背景是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 端到端验证把场景塞进决策模型里看误差统计指标都合格之后最有说服力的验证是把场景塞进一个真实的下游决策模型里。以机组组合为例做法是用削减后的场景集求解得到一阶段决策变量比如机组开停机计划把这个计划固定下来代入原始的3000个全场景计算期望运行成本同时用全场景直接求解如果能算得动的话得到理论最优成本对比两步结果的偏差。如果削减后得到的期望成本比全场景直接求得的成本偏差在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里的距离矩阵计算因为在不断删除场景的过程中距离矩阵需要反复更新。几个工程落地小建议用float32存储场景矩阵3000乘24的矩阵也能从72MB降到36MB对几十个风电场的场景集来说意义明显SBR里用numpy的cdist函数计算距离不要用纯Python循环场景削减是一次性的离线计算结果可以缓存下来反复用于不同下游模型不必每次重跑如果场景数到上万级别距离矩阵内存会吃不消这时考虑KDTree加速找最近邻或者分块处理距离计算。最后再分享一点个人体会。整套流程做下来我最大的感受是场景生成和削减从来不是两个步骤而是一整套可控的近似。生成阶段在近似真实概率分布削减阶段在近似生成出来的大量样本每一步都在牺牲一些细节关键是得知道自己在牺牲什么、牺牲得能不能接受。如果让我再来一次我会先把生成3000个、削减到10个这个端到端流程在小规模数据上完全跑通把所有检验指标建立起来再考虑t-Copula、DTW、深度生成模型这些更进阶的道具。基础链路稳了后面的一切才有可能谈。