1. 从蒙特卡洛到代表性场景这个课题到底在做什么先说明一下这个课题是我在一家设计院做新能源并网评估时接手的一个子任务核心一句话可以概括为用蒙特卡洛MC方法生成大量含时序相关性的风光出力场景再用削减算法把场景压缩到随机优化模型能消费的数量级。很多人一听“场景生成”就以为是随便抽个分布、生成一堆曲线但实际踩进去才发现真正让结果有价值的是“时序相关性”这四个字。风电和光伏出力不是独立白噪声而是有强自相关特性的连续过程太阳辐射从早到晚有平滑的抛物线趋势风机出力受气压场影响存在小时级的惯性。如果不把这些时间维度的依赖关系建模进去生成的场景会让后续的调度或规划算出一堆“看似合理、实际离谱”的结论。这个研究对象大致包括三类读者一是刚进入不确定性优化领域的研究生二是做新能源消纳分析的电网规划工程师三是需要在随机优化模型里灌入场景数据的算法工程师。我也算在坑里摸爬过的人今天把整个研究拆开讲讲从场景生成逻辑、相关结构建模到削减算法细节每个环节都聊一些不那么容易写进论文里的实操经验。2. 整体设计拆解为什么“时序相关性”是绕不开的关键点2.1 场景生成与削减的基本逻辑场景生成的核心目标是从随机变量的概率特性出发构造一组能代表未来可能性的时序路径。在电力系统随机规划中场景就是一组风电场站若干时段通常按 15 分钟或 1 小时的有功出力曲线。蒙特卡洛在这里的角色是抽样器先用历史数据估计出边缘分布和相关系数再按分布生成大量样本。常规做法是直接对多变量统计模型抽样生成 2000 甚至 10000 条场景但优化模型不可能把这么多场景全部放进去否则求解时间会变得不可接受。所以后半段必须有场景削减从原始集合里挑出少数代表性场景并给每个场景配一个概率权重使得缩减后的场景集合在统计特征、概率分布和时序结构上尽量接近原集合。不少人以为削减只是“选几条最像平均值的曲线”这是误区。场景削减本质是一个优化问题要在候选子集大小固定的约束下最小化缩减前后概率分布的差异。常见的削减算法包括聚类法和回代消除法我在后面的章节会详细展开。2.2 忽略时序相关会带来什么样的偏差时序相关性包含两层含义一是单条场景内部的时间自相关二是不同随机变量之间的互相关比如相邻风电场因为共用风带出力同时升同时降。很多初版方案只保证边缘分布正确用独立抽样生成数据结果画出来像是高频抖动的噪声甚至出现“上一时段风速为 12 m/s下一时段直接掉到 3 m/s”这种物理上几乎不可能的事件。这种场景送入随机机组组合模型后系统会因为虚假的高波动而过度配置旋转备用也会让储能调度结果严重失真。举个实际项目的教训某省电网在做新能源装机消纳评估时我们最初生成的风电场景没有处理时序相关性导致极端爬坡事件数量比真实历史多出将近 40%削峰填谷结论过于保守。后来加入时序相关建模生成的极端爬坡频次才回落到与历史观测基本一致。这个案例让我确定了课题的优先级先相关再削减顺序不能颠倒。3. 场景生成的核心实现与细节3.1 蒙特卡洛抽样前的准备工作场景生成不是一拍脑袋就能跑的第一步是数据清洗和格式化。我通常把历史数据抽成以日为单位的矩阵每一行是一个时间截面每一列是一个风电场或一个光伏电站。需要检查并处理数据中的异常跳变、缺测值部分站点还要做折算修正。然后对每列去除季节趋势和日周期均值例如光伏出力有明显的日出日落模式直接对原始数据拟合分布会把白天和夜间的混合分布搞得一团糟。我会先把出力按时段归一化到 [0,1] 之间的容量系数再分时段处理夜间光伏数据直接剔除不参与分布拟合。接下来是边缘分布拟合。对于风速韦布尔分布是常用选择对于风电功率由于有功率曲线非线性直接用 Beta 分布或者核密度估计更稳光伏出力用半经验分布也能凑合。我在项目中倾向于用非参数核密度估计一是很少需要对分布形态做假设二是后期逆变换抽样时直接基于累积分布函数插值方便又准确。拟合完成后用概率积分变换把每个变量从原始分布映射到标准正态分布这样便于下一步处理相关性。3.2 带时序相关性的 MC 抽样流程标准的多元相关正态抽样流程是给定相关系数矩阵 R计算其 Cholesky 分解 R L L^T生成独立标准正态随机向量 Z再用 X LZ 得到具有相关性的正态随机向量。但这里有个关键细节我们需要的不仅是空间断面相关还是时空联合相关。以 24 小时、5 个场站为例完整的相关矩阵是 120×120每个位置对应某一时刻某一场站。直接估计这个高维矩阵在历史数据不足时会产生严重的误差而且非正定的概率也很高。我采用的方案是分块处理加 AR(1) 平滑。第一步对每一时刻例如小时 t的多变量空间相关矩阵 R_t 做 Cholesky 分解生成该时段的空间相关正态样本第二步对每个场站的时序残差拟合一阶自回归模型形式是eps_t a * eps_{t-1} u_t其中 a 是自回归系数u_t 是独立正态噪声。把这两步结合生成的正态样本既满足同一时段不同场站的相关结构也满足相邻时段同一个场站的自相关结构。最后用逆累积分布函数变换把样本从正态空间映射回原始出力空间。整个流程用 Python 实现时关键代码如下import numpy as np from scipy.linalg import cholesky from scipy.interpolate import interp1d # 输入: KDE累积分布函数列表 cdfs, 逆CDF函数列表 inv_cdfs # 空间相关矩阵 R, 自回归系数 ar_coef n_scenarios 2000 n_periods 24 n_sites 5 L cholesky(R, lowerTrue) rng np.random.default_rng(2024) # 先抽独立正态变量 Z rng.standard_normal((n_scenarios, n_periods, n_sites)) # 对每个站点应用 AR(1) 系数同时加入空间相关的逐时段变换 corr_space np.zeros_like(Z) for t in range(n_periods): corr_space[:, t, :] Z[:, t, :] L.T scenarios np.zeros_like(corr_space) for t in range(n_periods): temp corr_space[:, t, :] if t 0: temp temp ar_coef * corr_space[:, t-1, :] # 标准正态转均匀再通过逆CDF转回原始分布 uniform_samples _standard_normal_cdf(temp) # 自行实现正态CDF for s in range(n_sites): scenarios[:, t, s] inv_cdfs[s](uniform_samples[:, s])注意代码里的ar_coef不是常数而是一个维度为n_sites的向量用来表示每个场站自己的时序惯性。不同场站的自相关强度差异很大光伏站点尤其明显白天时段相关性强夜间接日出时会产生阶跃所以对光伏站点可以分时段设不同系数。3.3 实操中的参数与工具选择场景生成阶段有四个参数最影响效果场景总数建议至少 2000条件允许可以到 5000。太少了尾部事件覆盖不住硬削之后概率权重会异常集中。随机种子固定种子是必须的否则两次生成结果不同后续对比实验会乱套。时间步长粒度1 小时粒度对中长期调度够用但做调频或储能秒级响应就必须用 15 分钟粒度。粒度越细相关矩阵越大计算时间和内存增长是指数级的建议结合模型需求选。分布拟合窗口如果历史存在明显的季节特征比如夏季和冬季风资源差异大不要把所有季节混在一起拟合应分季节或滑动窗口建模。工具方面Python 生态里numpy、scipy、scikit-learn基本够用statsmodels里的 AR 模型可以直接估计自回归系数。MATLAB 用户可以用mvnrnd生成多元正态样本但手动写 Cholesky 分解更可控。工程化部署时我会把整个流程封装成两个文件scenario_generator.py负责抽样scenario_reducer.py负责削减接口里统一传入参数文件这样多人协作时改动不会互相踩脚。4. 场景削减把上千个场景精简到可用的十几个4.1 削减算法的数学思想场景削减的任务是在原始场景集 Ω 中挑出一个子集 Ω_s并分配概率权重 p_i使缩减后的场景集合与原始集合的某种“概率距离”最小。最常用的两种路径是同步回代消除法和K-medoids 聚类法。同步回代消除法属于贪心策略每次从所有剩余场景中找到距离最近的一对场景删掉其中对整体分布贡献最小的一个并把它的概率转移到另一个上不断重复直到剩余场景数达标。具体步骤可以写成对所有场景初始化等权重 p_i 1/N。对任意 i≠j计算场景 i 与 j 的距离 d(C_i, C_j)。距离可以选欧氏距离、加权欧氏距离或动态时间弯曲DTW。找到概率加权后“对整体分布损失最小”的场景对即找到使p_i * d(C_i, C_j)最小的一对。将其中一个场景从集合中删除把概率累加到另一个场景上。重复步骤 24直到场景数等于预设值。聚类法的思路更直接把所有场景看成向量用 K-means 或 K-medoids 聚成 K 类以每个簇的质心作为典型场景簇内场景数量占比作为概率权重。K-means 实现简单、速度快但得到的质心可能不是真实存在的曲线对后续优化模型没有本质影响K-medoids 得到的场景都是原始场景中的真实样本物理含义更清晰代价是计算量更大。4.2 考虑时序相关性的削减策略这里要特别提醒如果把每条场景当作独立向量跑 K-means虽然保留了“同一时刻多个场站出力”的空间结构但很容易破坏场景内部的时间平滑性。举个例子聚类中心可能是“前半段像风电出力、后半段像光伏出力”的拼接曲线完全不符合真实过程。原因在于普通欧氏距离按时间维逐点求差对不同时间步的误差是等权累加的没有惩罚那种“时间上不连贯”的奇怪组合。针对这个问题我在项目中用了两步修正。第一步在距离计算时引入时序形状约束项比如给相邻时间段差分加入权重把距离改成d(C_i, C_j) sum_t w_t * (C_i[t] - C_j[t])^2 lambda * sum_t (D_i[t] - D_j[t])^2其中D_i[t] C_i[t1] - C_i[t]表示该场景在 t 时刻的爬坡率。这样能让削减后的典型场景不仅数值接近变化趋势也接近。第二步聚类过程不再对完整曲线直接 K-means而是先对每个场景做特征提取提取出“峰谷出现时间、平均爬坡率、日内最大值、波动方差”这四个特征再在特征空间里聚类。特征提取后的 K-means 结果往往比直接全向量聚类更稳健也能规避高维距离失效的问题。4.3 削减效果评估不要只看削减后的场景图“好看”就说方案有效。我习惯用三个量化指标对削减结果做体检均值与标准差相对误差对比削减前后各时段出力均值与标准差的偏差超过 5% 就要警惕。自相关函数ACF偏差画出削减后典型场景集合的平均 ACF 曲线与原始场景集合的平均 ACF 曲线对比主要看滞后 14 小时的相关系数是否接近。极端事件覆盖率把历史中实际出现过的极低出力或极高爬坡事件作为测试样本统计削减后场景集合中能覆盖的比例。覆盖率小于 80% 时说明尾部信息被削掉了。以我之前做的项目为例原始 3000 条场景削减到 20 条后均值偏差 1.8%标准差偏差 3.2%滞后 1 小时 ACF 偏差 4.1%极端爬坡覆盖率 89%完全可以支撑随机规划模型。如果某个指标超差优先调整距离权重 lambda 或增加最终场景数不要一上来就换算法。5. 常见问题与排查技巧实录做这类研究时容易踩的坑很多我整理了一份问题速查表基本覆盖了我自己踩过和帮人排查过的问题。问题现象典型原因排查思路与解决办法生成的场景出现负值或超过装机容量逆变换抽样时边缘分布边界外插不合理在逆 CDF 插值两端设置裁剪阈值将负值截断为 0把大于 1 的归一化值截断为 1相关系数矩阵非正定Cholesky 分解报错历史数据含缺失值或估计误差导致矩阵不满足正定条件对相关系数矩阵做特征值修正把所有小于零的特征值替换为一个小正数后重建矩阵ACF 偏低场景时序波动过大直接用独立正态样本逐时段抽样没有加入 AR 过程按 3.2 节方式加入时间滞后项或改用多元 AR(1) 模型整体抽样削减后典型场景“四不像”距离函数只考虑数值差异忽略了时间趋势在距离中加入差分项或改用 DW 系数动态时间弯曲距离削减后极端事件变少聚类时 K 值太大或概率权重分配不合理增加一个用于评估尾部覆盖率的指标必要时对极值场景单独保留运行结果在不同机器上不一致算法依赖全局随机数而不固定种子所有随机流程设置固定种子并在配置文件里记录当前种子值场景数量增加但优化结果没改善削减后的场景集合已经足够表达关键统计特征继续加数量边际收益很小结合迭代曲线确定拐点一般 2040 个场景对大多数随机 MILP 模型已经足够其中非正定矩阵是最容易让新人卡死的问题。我通常用一个函数做最近正定矩阵修正对特征值分解后的对角矩阵把所有小于 1e-6 的特征值替换成 1e-6然后重构矩阵并重新归一化对角线。这个修正后的相关系数矩阵和原始矩阵在大多数业务场景下几乎无差别但 Cholesky 分解会非常稳定。另一个高频问题是逆变换后的分布尾部分叉。因为核密度估计在样本边界外可能外推不准导致逆变换抽样时出现极端大值。我的处理方式是对每条场景做物理合理性检查风电出力不能超过额定功率光伏夜间出力接近 0爬坡率不能超过场站允许的最大爬坡速率。这些检查放在生成后、削减前成本低且能显著减少后续优化的不可行问题。6. 项目落地的一点经验如果只让我总结一条经验那就是时序相关性的价值在随机优化中会被放大但必须在场景削减前把相关结构注入到场景里而不是指望削减阶段再修复。削减算法是一个压缩过程任何被原始场景丢失的信息都很难在压缩后找回来。所以我在项目排期上把大块时间放在数据清洗和相关矩阵估计上削减本身反而是费时最少的环节。还有一个容易被忽略的小技巧记得把削减后的场景概率保存成单独的字典并且在优化模型里用场景概率做权重系数而不是对每个场景简单求平均。很多人在随机规划建模时把场景当作等权重样本去拟合一个确定性模型这是错误用法。削减算法给出每个典型场景对应的概率集合优化目标中每一项都要乘以它对应的概率才能真实反映不确定性影响。最后整条流水线我建议做成“配置项驱动”。把历史数据路径、分布参数、场景数量、削减数量、距离权重等全部放进一个 YAML 文件跑实验时只改配置不改代码。这样后面对比不同削减方案时能保证除测试变量外的一切条件完全一致。这个方法帮我省了不知道多少重复劳动今天一并分享出来希望后来者少走弯路。