做新能源并网的人大概都有过这种体验早上看气象预报风速曲线很漂亮功率预测模型也给出了一个看着很稳的出力曲线。结果下午风向一转风电场实际出力直接打了七折。调度电话紧跟着就来了——你的预测偏差怎么这么大备用容量不够了。说白了风光出力在数值天气预报解算完成的那一刻并不是一个确定的值而是一团“可能的状态”。这跟薛定谔的猫很像盒子里的猫在揭晓之前既是活着也是死的某个风电场在 t 时刻的出力在实测数据回传之前也同时叠加着“高出力”与“低出力”两种可能。蒙特卡洛场景生成就是把这团“叠加态”展开成一组离散的、可计算的可能世界再拿这组场景去跑优化、评估风险、配置备用。这篇文章我会完整复盘一个我经手的风光联合场景生成项目为什么要用场景而不是单点预测、蒙特卡洛抽样背后的数学直觉、以及一套可以直接跑起来的 Python 实现从误差建模、相关性处理到场景削减最后聊聊我在工程落地时踩过的几个坑。适合新能源场站、电力市场、储能配置和电网调度方向的朋友参考。1. 为什么风光功率预测需要一张“场景集合”1.1 “薛定谔的猫”与功率预测的叠加态在风光发电的功率预测任务里我们通常只拿得到一个输出一条未来 96 个点15 分钟分辨率的确定性出力曲线。它来自数值天气预报NWP数据经过物理模型或统计模型折算成功率。问题在于这条曲线本质上只是“给定当前输入信息下的条件期望”——也就是说它是在用平均值去代表一堆分散的可能性。真实世界里有太多东西让预测偏离期望NWP 对风速的预报天然有误差地形对气流的扰动很难建模光伏组件温度和辐照度的换算存在偏差风机可用率受检修和故障影响甚至电网限电都会让出力曲线变形。这些因素叠加起来使得实际出力围绕预测值呈现一个分布而不仅仅是一个点。这就引出场景生成的核心动机既然无法把“猫”的生死状态精确预言那就把“活着”和“死了”这两种结果以及它们出现的概率都摆到台面上。对功率预测而言就是给出一组未来出力曲线每条曲线附带一个发生权重共同构成对未来不确定性的完整描述。1.2 确定性预测的短板一次调度复盘我印象很深的一次复盘来自某个风电场某天的晚高峰时段。当日 NWP 显示 20 点风速会起来预测出力 100MW调度按这个值预留了火电备用。结果实际风速比预报低了近 2m/s20 点实发只有 40MW缺口 60MW。好在当时联络线备用够宽裕火电紧急加出才没酿成事故。事后我重新用场景方法做了一次回测。基于当日的预测误差分布采样生成 500 条场景后发现 20 点出力低于 50MW 的场景占比大约是 15%。这个信息如果提前摆在调度面前意味着至少有 15% 的概率需要额外加备 50MW 以上的火电容量——准备预案的力度和方式都会完全不同。单点预测丢掉的东西恰恰是风险管理者最关心的分布尾部。尾部描述的是“坏情况有多坏、出现概率有多大”而场景集合恰好把这些信息完整保留了下来。储能容量配置、电力现货市场申报、随机机组组合、备用容量评估这些场景下游能用到的地方几乎都与“如何处理不确定性”有关。2. 蒙特卡洛场景生成的核心思路2.1 蒙特卡洛到底在“蒙”什么很多人第一次听到“蒙特卡洛”四个字会以为是在搞随机乱猜。实际上它做的事情非常有章法当解析方法求不出概率分布时就通过大量随机抽样用样本统计量逼近真实分布。大数定律保证了当抽样次数足够多时样本均值会收敛到期望样本分位数会逼近真实分位数。举个生活化的例子想知道一个池塘里鲤鱼和草鱼的比例又不想把水抽干就随机撒一网捞上来一千条数一下比例。只要撒网位置足够随机、捞得足够多这个比例就非常接近真实比例。蒙特卡洛场景生成做的就是类似的事——把未来的“不确定性池塘”反复抽样得到一组可能出力轨迹用它们的分布代表真实的出力分布。在场景生成任务里蒙特卡洛三要素缺一不可一是高质量随机数二是准确的误差分布模型三是变量间的相关性结构。随机数只是工具误差模型才是灵魂相关性结构决定场景是否“像真实世界”。2.2 误差分布建模场景生成的地基误差分布建模是整个场景生成流程里最不该省功夫的一步。常见的建模方法分为加性误差和乘性误差两类加性误差实际出力 预测出力 误差项。适用于出力较大、误差尺度相对稳定的时段。乘性误差/相对误差实际出力 预测出力 ×1 相对误差项。更适用于低出力时段尤其是光伏夜间出力接近零时。如果对光伏的午前时段使用加性正态误差很容易抽出负出力这是最典型的低级错误。我建议先转换为相对误差序列即 e_t (P_actual,t - P_forecast,t) / P_capacity把误差按装机容量归一化再做分布拟合效果会稳健很多。风功率误差通常接近正态分布但尾部偏厚光伏功率误差受云层影响往往呈现双峰形态——一个峰对应晴空、一个峰对应云遮挡。面对这类数据用单一正态分布拟合显然不够可以考虑混合高斯分布、偏正态分布或者干脆用核密度估计这类非参数方法。原则上没有“万能分布”一切以数据说话务必做分布检验。2.3 从独立抽样到“像那么回事”时空相关性如果对每个时刻独立抽样误差生成的场景会像一把高频抖动的毛刷子——前一秒 80MW下一秒 20MW再下一秒 70MW。现实中的风电场出力曲线是平滑变化的风速不会在一分钟内翻天覆地。所以必须对时间相关性建模。一个非常实用的方案是用一阶自回归模型 AR(1)ε_t φ × ε_{t-1} δ_t。其中 φ 是相邻时刻误差的相关系数可以从历史误差序列里用最小二乘估算出来。φ 接近 1 时场景曲线很平滑接近 0 时恢复成白噪声。另一个问题是风与光的空间相关性。同一片区域里风电场和光伏电站往往受同一个天气系统影响风速大时云通常也少风、光出力可能存在某种此消彼长的关系而相邻的风电场之间出力又有很强的正相关。要刻画这种“联动”可以用 Cholesky 分解如果已知变量间的相关矩阵 C先对 C 做乔列斯基分解得到 C L × L^T。然后先生成一组独立标准正态样本 Z再用 Z 乘以 L 的转置得到的样本就带上了矩阵 C 描述的相关性。这个过程可以理解为给一群独立行动的合唱队员下发统一指挥指令相关矩阵就是那个指挥棒让所有声部按约定的节奏协同。2.4 场景削减从一千个世界到十个世界蒙特卡洛生成的初始场景数量通常是几百甚至上千个。下游如果要做随机优化1000 个场景 × 96 个时段 × 几十台机组的组合计算规模会迅速爆炸。因此需要用场景削减算法从大场景集中挑出少量有代表性的场景并给每个场景分配权重使削减后的场景集合与原集合在概率分布意义上尽量接近。方法原理概要适用场景优点缺点K-means 聚类按距离把场景聚成 K 类取簇中心均值型优化简单、快速容易丢失尾部极端场景快速后向削减FFS基于 Wasserstein 距离逐步删除“最不孤单”的场景并将其权重转移给最近邻风险约束评估保留概率分布尾部实现稍复杂、计算耗时分层抽样按分位数分层采样需要极少场景数速度快对相关结构表达能力弱我的经验是下游关注的是期望收益或平均成本选 K-means 够用下游有风险约束、需要评估极端场景老老实实用 FFS。场景削减不是把所有场景“平均”成几条中庸曲线而是把概率分布的骨架保留下来。3. 实战一个用 Python 实现的场景生成全流程3.1 数据准备与误差序列构造先列出我用的技术栈pandas 做数据处理numpy 做随机抽样和矩阵运算scipy.stats 做分布拟合matplotlib 做可视化。数据只需两张表每张表的每一行对应一个时段实际出力表风电场 01 / 风电场 02 / 光伏站 A 的实际功率。预测出力表与上述场站对应的日前预测功率。数据清洗这一步非常关键。限电时段要剔除因为限电期间实际出力被调度指令“钳住”在上限附近不反映真实的气象驱动波动异常零值要排查传感器故障或通讯中断会产生成片的 0 值混入误差序列会把分布严重拉偏。当时我处理的风光联合项目里有 3 个风电场和 2 个光伏电站时间分辨率 15 分钟历史数据一整年。场景生成的时间范围是未来 24 小时也就是 96 个点。import pandas as pd import numpy as np # 假设 capacity 是各场站装机容量(单位MW) capacity { wind_01: 150, wind_02: 100, wind_03: 120, pv_01: 80, pv_02: 60 } # 读取数据 actual pd.read_csv(actual_power.csv, index_col0, parse_datesTrue) forecast pd.read_csv(forecast_power.csv, index_col0, parse_datesTrue) # 剔除限电时段实际出力连续多点为装机上限附近且预测也偏高时很可疑 # 此处简化为剔除实际出力达到装机容量98%以上的长时段 limit_mask (actual / pd.Series(capacity)).rolling(3h).mean() 0.98 actual_clean actual[~limit_mask] forecast_clean forecast[~limit_mask] # 构造相对误差序列误差/装机容量 err (actual_clean - forecast_clean) / pd.Series(capacity) err err.replace([np.inf, -np.inf], np.nan).dropna(howall) print(err.describe())相对误差做归一化后不同装机规模的场站放在同一尺度下比较才有意义。后面抽样时再把误差乘回对应装机容量还原成功率偏差值。3.2 误差分布拟合与参数标定拿到清洗后的误差序列下一步是选分布、拟合参数。我一般先画直方图和 Q-Q 图看一眼峰度、偏态再试几个候选分布用 Kolmogorov-Smirnov 检验或 AIC/BIC 做筛选。from scipy import stats # 以某风电场为例 err_wind err[wind_01].dropna() # 正态拟合 mu, sigma stats.norm.fit(err_wind) print(Normal fit:, mu, sigma) # t分布拟合 df, loc, scale stats.t.fit(err_wind) print(t fit:, df, loc, scale) # 偏态正态拟合适合光伏 skew_mu, skew_sigma, skew_alpha stats.skewnorm.fit(err_wind) # 简单检验 ks_stat, ks_p stats.kstest(err_wind, norm, args(mu, sigma)) print(KS p-value:, ks_p)实际项目中风电场的误差序列我最终选了 t 分布因为尾部更厚能覆盖风速骤变的极端情形光伏电站则用了混合高斯模型——一个分量代表晴空状态另一个代表云遮挡状态。混合模型没有现成的 scipy 封装我用 EM 算法手写了一个简版但如果你不想折腾直接用偏态正态拟合也够用。拟合完成后务必画一下拟合曲线和真实直方图的叠加图。数值检验指标再好眼睛这一关过不了后面场景出来的形态也会有惊喜。3.3 生成带相关性的初始场景误差模型标定好之后开始蒙特卡洛抽样。生成场景的核心分四步先按各场站的误差分布抽样再施加时间相关性平滑然后用 Cholesky 分解施加空间相关性最后还原成功率曲线并做边界裁剪。np.random.seed(42) n_scenes 1000 n_steps 96 n_sites len(capacity) # 步骤1: 从各误差分布抽取独立性标准化样本 # 这里为了演示,以正态或t分布为例,实际项目按场站用不同分布 base_samples np.zeros((n_scenes, n_steps, n_sites)) site_params { wind_01: (t, df, loc, scale), wind_02: (norm, mu, sigma), # 其余同理此处省略 } for s_idx, site in enumerate(capacity.keys()): dist_name site_params[site][0] if dist_name t: base_samples[:, :, s_idx] stats.t.rvs( site_params[site][1], locsite_params[site][2], scalesite_params[site][3], size(n_scenes, n_steps) ) else: base_samples[:, :, s_idx] stats.norm.rvs( locsite_params[site][1], scalesite_params[site][2], size(n_scenes, n_steps) ) # 步骤2: AR(1)时间平滑, 施加在时间维度上 phi 0.85 # 从历史误差序列估计得到 smooth_samples np.zeros_like(base_samples) smooth_samples[:, 0, :] base_samples[:, 0, :] for t in range(1, n_steps): smooth_samples[:, t, :] phi * smooth_samples[:, t-1, :] base_samples[:, t, :]注意 AR(1) 递推会让误差序列的标准差发生变化严格来说要除以 sqrt(1 - φ²) 做归一化修正。实际操作中如果 φ 是从历史残差序列估出来的那个残差序列本身就已经包含了自相关所以直接用历史残差的分布参数去抽样反而无需额外修正——这也是我后来踩过坑才想明白的细节。空间相关性的处理放在平滑之后。用各场站误差序列的历史数据估计相关矩阵再 Cholesky 分解然后对同一时刻所有场站的样本施加相关结构。# 估计相关矩阵用历史误差序列 corr_matrix err.corr().values # 步骤3: Cholesky分解施加空间相关性 L np.linalg.cholesky(corr_matrix) # 把平滑后的样本在站点维度上施加相关性 for t in range(n_steps): # 抽取当前时段的样本矩阵 (n_scenes, n_sites) X_t smooth_samples[:, t, :] # 左乘 L.T: 让列与列之间带上相关结构 # 注意: 如果 X ~ N(0, I), 则 X L.T ~ N(0, L L.T) N(0, corr) smooth_samples[:, t, :] X_t L.T # 步骤4: 还原成功率和边界裁剪 scenes {} for s_idx, site in enumerate(capacity.keys()): site_scenes forecast[site].values.reshape(1, -1) smooth_samples[:, :, s_idx] * capacity[site] site_scenes np.clip(site_scenes, 0, capacity[site]) scenes[site] site_scenes # 合并成 (n_scenes, n_steps, n_sites) 数组 scenes_array np.stack([scenes[site] for site in capacity.keys()], axis-1) print(scenes_array.shape)关于整体负荷出力——如果你要的是“风光联合出力”的单条场景直接把所有场站的场景逐时点相加即可。但要注意总出力不一定能简单相加后裁剪到 0 到总装机之间因为光伏和风电的构成不同极端叠加时仍可能出现超过总装机的情况这时候裁剪逻辑要放到求和之后再做一次。3.4 场景削减基于后向削减算法实现1000 个初始场景直接扔给下游优化器会被骂回来我通常削减到 20 个左右。下面给出快速后向削减FFS的 Python 实现核心逻辑是场景之间的距离用欧氏距离或 Wasserstein 距离衡量迭代删除“离最近邻距离 × 权重”最小的场景把它的权重转移给最近邻。def scenario_reduction(scenes_2d, n_keep): scenes_2d: 形状为 (n_scenes, n_features)每个场景展平成一维向量 n_keep: 削减后保留的场景数 返回: 保留场景的索引数组和对应权重 n scenes_2d.shape[0] weights np.ones(n) / n idx list(range(n)) while len(idx) n_keep: min_cost np.inf min_i -1 min_j -1 for i in range(len(idx)): best_d np.inf best_j -1 for j in range(len(idx)): if i j: continue d np.linalg.norm(scenes_2d[idx[i]] - scenes_2d[idx[j]]) if d best_d: best_d d best_j j cost weights[idx[i]] * best_d if cost min_cost: min_cost cost min_i i min_j best_j # 把被删场景的权重加到最近邻上 weights[idx[min_j]] weights[idx[min_i]] # 从集合中删除该场景 idx.pop(min_i) # 归一化权重 final_weights weights[idx] final_weights final_weights / final_weights.sum() return np.array(idx), final_weights # 使用示例: 把 (n_scenes, n_steps, n_sites) 展平成二维 n_scenes scenes_array.shape[0] scenes_2d scenes_array.reshape(n_scenes, -1) keep_idx, final_weights scenario_reduction(scenes_2d, n_keep20) print(保留场景索引:, keep_idx[:5], 权重示例:, final_weights[:5])这个实现是 O(n²) 复杂度1000 个场景要循环 980 轮每轮都算全量两两距离实测大概几分钟可以接受。如果场景规模到 5000 以上记得用空间索引加速或者降维后再削减。关于场景之间的距离度量直接用 96×5 维欧氏距离够用但如果你想更强调分布形状可以先把场景按分位数归一化再算距离。历史回测中我发现直接用原始功率值计算时出力大的场景天然距离远容易被优先删除导致低估出力水平。后来改成先对每条场景做标准化减自身均值、除自身标准差再算距离效果好了很多。3.5 结果可视化与导出削减完的场景要画图给业务方看。一张灰线为全部初始场景、彩线为削减后代表场景的图比任何解释都有说服力。另外我习惯再加两张图一张是风光联合出水的分位带图把 5%、25%、50%、75%、95% 分位数画成色带另一张是削减前后概率分布的对比图证明削减后的场景集“长得很像”原始集。import matplotlib.pyplot as plt # 以某时段的风光联合出力为例 total_scenes scenes_array.sum(axis-1) # 总出力场景 keep_total total_scenes[keep_idx] # 置信带 low_5 np.percentile(total_scenes, 5, axis0) low_25 np.percentile(total_scenes, 25, axis0) high_75 np.percentile(total_scenes, 75, axis0) high_95 np.percentile(total_scenes, 95, axis0) plt.figure(figsize(12, 5)) t np.arange(n_steps) * 15 / 60 # 小时 plt.plot(t, keep_total.T, alpha0.8, linewidth1.2, colorC0) plt.fill_between(t, low_5, high_95, alpha0.15, colorC0) plt.fill_between(t, low_25, high_75, alpha0.25, colorC0) plt.xlabel(Hour) plt.ylabel(Total Power (MW)) plt.title(Generated Scenarios after Reduction) plt.grid(alpha0.3) plt.savefig(scenarios.png, dpi150)导出环节我建议用 CSV 分两列场景 ID、权重。场景曲线本身存成一个二维表每行一条场景每列一个时段。务必不要丢掉权重很多下游工程师默认场景是等权的实际上削减后的场景权重差异很大等权处理会让期望值计算完全错误。4. 场景质量校验与工程落地经验4.1 我踩过的三个坑第一个坑是负出力。最初我图省事直接对全时段用加性正态误差建模光伏夜间场景里出现了一大堆负值裁剪到 0 后又让凌晨时段的场景分布整体左偏。后来改成相对误差建模并且拟合前先按功率区间分段处理负出力问题才彻底根治。第二个坑是削减过度造成尾部丢失。某次为了下游计算速度把场景从 1000 个直接削到 5 个。结果恰巧把最坏的一组极端场景全部剔除导致备用容量评估偏于乐观。复盘后我把削减数量定在 1520 个并专门检查削减前后 5% 分位数和 95% 分位数是否显著漂移。如果漂移超过 5%就说明削多了或者选的距离度量不对。第三个坑是时间相关性没处理好。第一版场景生成每个时点的误差完全独立抽样出来的风光联合出力曲线一天里面上下跳动几十次看着就像噪声。下游机组组合优化直接不收敛我一开始还以为是求解器的问题后来画出场景曲线才发现是输入数据质量太差。引入 AR(1) 平滑后问题立刻消失。4.2 场景质量怎么验收场景生成不是画完图就完事了要建立量化验收标准。我常用的三个指标经验覆盖率把历史实际出力序列拿出来统计它落在场景集合区间内的比例。比如构建 90% 置信区间理论上应有约 90% 的实际点落入区间。覆盖率明显低于标称值时说明场景集合太“窄”低估了不确定性。平均带宽/锐度区间不能为了覆盖率无限放宽。平均带宽越小表示场景越“锐利”信息量越大。经验上同样的覆盖率下带宽越小越好。PIT 检验把实际值代入场景集合的经验 CDF得到 PIT 值。如果场景概率分布校准良好PIT 序列应近似服从均匀分布Kolmogorov-Smirnov 检验可以直接验证。除此之外我还有一个笨办法拿生成场景回测历史上的随机优化决策对比实际执行结果。中间也许会暴露一些概率分布之外的问题比如相关性结构建模不准确导致联合出力的极端组合被低估。4.3 工程层面的几点建议随机种子一定要可配置。场景生成是随机算法如果每次跑出来的结果都不同下游工程师没法复现你的结果也不会信任你的系统。我在项目里用统一的 seed 管理模块每个场站、每次抽样都记录种子保证回测可复现。场景数量和削减数量不必贪多。1000 个初始场景通常是足够的再多边际收益很低削减到 20 个左右已经能覆盖绝大多数优化需求。如果下游只需要期望值配合 K-means 聚类削减10 个场景都很稳。最后是数据闭环。场景生成做完后要持续跟踪实际出力定期重新标定误差分布参数。气象模式换版本、场站加装新风机、光伏组件衰减都会让误差分布漂移。我习惯每个月重新拟合一次分布参数每季度做一次全量回测校验确保整个流程始终贴合最新运行数据。项目上线后的几点体会这个项目做完我最大的感受是场景生成最难的部分其实不在算法而在让下游团队真正信任“集合”这个东西。单点预测错了大家都会说“预测不准”可场景集合哪怕再准第一次拿到的人总会问一句“你到底准备用哪条曲线”。我的经验是先把分位带和代表场景推给调度侧让他们把思维从“一个值”切换到“一个区间”。用顺手之后再做完整的随机优化、储能配置、市场申报。这套框架后续可以扩展的方向很多把气象集合预报成员直接作为种子场景把场景生成嵌入储能容量配置的随机优化甚至配合交易策略做鲁棒性回测。但这些都是后话先把场景生成跑起来、用起来比什么都重要。