做风光模拟的时候很多人最开始的想法就是多抽几个场景几千个不够就抽一万个保证把不确定性都覆盖到。但真到了做随机优化、做生产模拟的时候一万个场景扔进去求解器直接罢工好不容易跑出个结果等得黄花菜都凉了。这时候就绕不开一个操作场景削减。说白了就是拿少数几个有代表性的场景去近似替代原来那成千上万个场景同时让概率分布特征尽量不丢。这篇文章我就围绕基于正态分布的风光模拟和场景削减把整个技术链路怎么搭、算法怎么选、代码怎么写、坑在哪里一次讲透适合正在做微电网调度、电力系统随机规划或者风光出力不确定性分析的朋友参考。1. 风光模拟与场景削减先把要解决什么问题讲清楚1.1 为什么需要场景不确定性建模到底在建什么风光出力的本质是一个随机过程。风速受气压、地形、昼夜温差影响光照受云层、季节、大气透明度影响这些因素叠加起来让光伏和风电在任意时刻的出力都带有强烈的不确定性。我们在做系统规划或者调度策略时不能只按“晴天满发”这种理想情况算否则实际运行一遇到多云天气功率预测偏差就能让你安排的机组组合全部失效。场景其实就是把这种不确定性离散化。把连续的概率分布变成一组离散的样本每个样本代表一种可能的风光出力时序。比如风电出力的预测值是30MW实际可能是25、28、35服从某种分布我们就按这个分布抽很多个可能的出力曲线每条曲线都是一个场景。理论上抽得越多覆盖的分布空间越完整但优化问题的规模也随之爆炸。1.2 为什么一定需要削减计算代价和质量之间的平衡我见过有人直接用5000个场景跑一个两阶段的随机优化结果两三天没跑完最后只能砍需求。场景规模对计算复杂度的影响是灾难性的因为随机规划问题里每个场景都对应一组决策变量和约束条件。假设你原来确定性模型的变量是1000个用500个场景做随机规划变量数直接膨胀到50万量级普通内存的机器根本扛不住。场景削减就是在这时候上场。它的数学本质是从原始的大规模场景集里选出一个子集再给每个保留场景重新分配概率使得削减前后场景集的概率分布差异最小。好的削减方法能用几十个场景保留几千个场景的主要统计特征均值、方差、相关性都能对得上。1.3 整个技术路线其实就三步我做风光模拟加场景削减的完整路线是固定的先做数据清洗和分布拟合再做蒙特卡洛抽样生成海量场景最后用削减算法挑代表场景。三者缺一不可。分布拟合不准后面抽出来的场景全是垃圾抽样数量不够削减完的代表性也打折扣削减算法不合适得到的场景集在优化问题里会和真实场景相差很大。所以整篇文章的推进思路就是按这个顺序来每一步我会把背后逻辑讲清楚再给出可以直接复制的代码和参数建议。2. 正态分布风光模拟的数学底座2.1 为什么是正态分布不是硬套是有依据的先泼一盆冷水风光出力的原始分布并不是严格的正态分布。风速一般更接近威布尔分布光照辐照度在晴天时接近Beta分布。那为什么标题还叫“基于正态分布的风光模拟”因为在大多数实际工程场景里我们并不直接对出力值本身建模而是对预测误差建模而预测误差在中心极限定理的意义下往往可以合理近似为正态分布。比如风电功率预测系统给出明天某个时刻出力是50MW实际的出力等于预测值加上一个随机误差项。这个误差受很多独立随机因素影响按中心极限定理大量独立随机变量之和趋向正态分布。所以我们建模时写成P_actual P_forecast epsilon, epsilon ~ N(0, sigma^2)注意功率本身的正负是有边界的0到装机容量之间所以对误差做正态抽样之后必须做截断处理把超出边界的情况截掉或者重新采样。这个细节很多初学者会忽略导致模拟出的场景出现负功率一看就是假的。从数据拟合的角度正态分布的好处是参数少、含义直观就均值和标准差两个量工程上非常容易估计。你只需要统计历史预测误差的均值和方差就可以搭起整个模拟框架。2.2 参数估计不能光看整体要按时段拆有些人拿到一整年的历史数据哗啦啦全丢进去算出均值和方差就开始做模拟。这样做出来的东西基本上是废的。因为风光预测误差有明显的时段特征白天光伏出力大误差绝对值和方差都大凌晨光伏不出力误差方差趋近于零。如果混在一起拟合相当于用一个大方差去覆盖所有时段白天的场景过于保守夜间的场景又过于激进。正确做法是按时间段拆开统计。以1小时间隔为例一天24小时就分24个时段或者按季节再加一层拆。对每个时段分别估计误差的均值 mu_t 和标准差 sigma_t。实际操作中可以写成函数对DataFrame按小时分组计算然后做平滑处理避免某些时段样本量太少导致方差畸高或畸低。我一般是用指数加权移动平均做平滑让相近时段的标准差不会突变。这样拟合出来的参数更能反映实际的物理过程光靠白天数据多就把标准差拉大是不合理的。2.3 相关性是个隐秘的坑单独对风、对光分别做正态抽样会忽略风光之间的相关性。比如同一片区域内风电偏大往往对应阴天或大风天气光伏输出可能下降反过来强日照天气光伏大发时风速往往偏低。如果模拟时光是独立抽样就可能出现“又大风又大太阳”的场景这在物理上很罕见但在数学上完全可能发生。处理办法有两个思路。一个是用Copula函数对风光的联合分布建模另一个更常用的是在抽样时引入相关系数矩阵通过Cholesky分解生成相关正态随机数。Cholesky分解的思路是把相关的多维正态分布拆成下三角矩阵乘法对一个独立的标准正态随机向量做线性变换就得到期望的相关结构。这段逻辑在代码里也就几行的事但效果非常显著能直接把模拟场景的物理合理性拉上一个台阶。3. 场景削减的核心算法从K-means到同步回代消除3.1 场景相似性的度量削减的第一步其实是“算距离”不管是哪种场景削减算法都绕不开一个问题怎么判断两个场景像不像。最常用的就是欧氏距离加L2范数。假设每个场景是T个时段的一条曲线场景i和场景j之间的距离就可以定义成所有时段功率差值的平方和再开根号。D(i, j) sqrt( sum_{t1}^{T} (P_i(t) - P_j(t))^2 )这个距离定义简洁、计算高效对时序曲线的整体形状差异很敏感。不过它也有缺陷如果某个时段的误差特别大平方项会把距离拉得很大掩盖其他时段的相似性。所以有的场景集在做距离计算之前会先对数据进行标准化把每个时段的出力除以该时段的容量或者历史标准差相当于给不同时段加了权重避免那些绝对数值大的时段统治距离计算。如果场景携带概率权重比如每个场景本身有不同的发生概率p_i那么在算法迭代中还要用到加权距离。加权距离的计算不是简单乘一下而是要考虑场景代价的累积这个在同步回代消除里体现得尤为明显。3.2 方法一K-means聚类场景削减K-means可能是上手最快、工程里最喜欢用的场景削减方法。核心思路很朴素把N个原始场景通过聚类压缩成K个簇每个簇的中心就是保留场景簇内所有场景的概率之和就是保留场景的新概率。具体步骤分四步从原始场景中随机选择K个场景作为初始聚类中心对每个场景计算到K个中心的距离把它归到距离最近的中心所在的簇对每个簇内所有场景取平均得到新的中心重复第2和第3步直到聚类中心不再变化或者达到最大迭代次数。这里有一个常规文档里不会细说的重点聚类完成后中心曲线是一个平均曲线它并不一定等于原始场景中的某一条。但在场景削减的实际使用场景里我们通常希望保留场景是可回溯的、和原始数据同构的曲线。所以更稳妥的做法是保留每个簇中距离簇中心最近的那个原始场景作为代表而不是直接使用平均中心。K值的选取也是一个权衡。K太小场景多样性不足极端场景全被抹掉K太大削减又没达到降低计算复杂度的目的。我一般先用轮廓系数圈定一个大致范围再结合随机优化问题的具体规模反推K的上限。比如决策变量总数不能超过求解器能处理的量级在这个上限以内尽量选大的K。3.3 方法二同步回代消除后向消除同步回代消除是文献里出现频率很高的另一种削减算法它的思路和K-means完全不同不是聚类求中心而是从原始场景集里不断“删”场景一直删到剩下目标数量为止。每一步都删除一个对整体概率分布影响最小的场景并把被删除场景的概率加到离它最近的场景上去。算法步骤我整理了一下计算原始场景集中所有场景两两之间的距离形成距离矩阵对每个场景i找到距离它最近的场景j这个最近距离记为d_i在所有场景中找d_i最小的那个场景它的意义是删掉它对整体分布扰动最小把这个场景的概率累加到它最近的那个场景上然后从场景集中移除它更新距离矩阵重复步骤2到4直到剩余场景数等于目标K。单看文字有点干但代码实现起来并不复杂。它的一个优势是每一步都在做全局最优的局部决策得到的保留场景一定是原始场景里的真实曲线物理上可解释性非常强。缺点是在场景数量大的时候计算量偏高因为每删一个场景都要更新一次距离矩阵时间复杂度接近O(N^3)量级。我在实际应用里一般把两种方法配合使用先用K-means粗筛从几千个场景快速压到几百个再用同步回代消除精修从几百个压到最终需要的几十个。两步走既能利用K-means的速度优势又能利用回代消除的精度优势。3.4 两种方法的对比与选择K-means和同步回代消除没有绝对的好坏关键看应用场景。K-means的优势是速度快、实现简单、处理超大规模场景集很轻松不足是结果受初始中心影响有随机性而且平均中心曲线会抹平极端情况。同步回代消除的结果更稳定、更能保留原始场景形态但计算复杂度高场景数和维度一旦上去耗时明显增加。如果你的原始场景集只有几百个直接用同步回代消除完全没问题。如果你蒙特卡洛抽了一万个场景建议别直接跑回代先K-means粗降到300到500个再回代精修到目标数量。这个组合拳在很多论文里都有影子工程上确实好用。4. 完整实操Python实现风光模拟与场景削减4.1 准备模拟数据先造一组实验数据演示完整流程。假设有一个48时段的风电出力预测值序列真实出力等于预测出力加正态误差标准差设置为预测值的10%左右同时做0到装机容量的截断。这里的真实出力是我们用来验证的基准实际应用中没有真实出力只有预测值和误差分布。import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.stats import norm from sklearn.cluster import KMeans np.random.seed(42) # 生成风电预测序列48时段容量200MW T 48 capacity 200 t np.arange(T) # 模拟一个峰谷明显的预测曲线 forecast 120 60 * np.sin(t / T * 2 * np.pi - np.pi/2) forecast np.clip(forecast, 20, 180) # 误差标准差为预测值的10% sigma 0.10 * forecast # 生成5000个随机场景 n_scenarios 5000 errors np.random.normal(loc0, scale1, size(n_scenarios, T)) * sigma scenarios forecast errors scenarios np.clip(scenarios, 0, capacity) # 每个场景等概率 probs np.ones(n_scenarios) / n_scenarios注意代码里的一个细节我先用标准正态随机数乘以逐时段的sigma而不是一次性np.random.normal传入整个sigma矩阵。两种写法结果等价但前者更直观方便你以后扩展时对每个时段用不同的分布类型。4.2 从K-means到精修场景集从5000压到50先做K-means快速粗削减到300个场景再用同步回代消除把300个场景精修到50个。K-means部分可以直接调sklearn重点在于聚类完成后按簇内最近原始场景提取保留场景。# K-means粗削减 n_clusters 300 kmeans KMeans(n_clustersn_clusters, random_state42, n_init20) labels kmeans.fit_predict(scenarios) # 每个簇内找离中心最近的原始场景作为保留场景 reduced_scenarios [] reduced_probs [] for c in range(n_clusters): idx np.where(labels c)[0] center kmeans.cluster_centers_[c] dists np.linalg.norm(scenarios[idx] - center, axis1) rep_idx idx[np.argmin(dists)] reduced_scenarios.append(scenarios[rep_idx]) reduced_probs.append(np.sum(probs[idx]))同步回代消除的实现里最核心的就是距离矩阵的维护和删除操作。网上能搜到的版本里很多没有注意“被删除场景概率要加到最近场景”这个细节导致削减后概率总和不等于1后面做优化就全乱了。我提供一个验证过概率守恒的版本def synchronous_backward_reduction(scenarios, probs, target_count): scenarios scenarios.copy() probs probs.copy() n scenarios.shape[0] while n target_count: # 两两场景距离 dist_matrix np.linalg.norm(scenarios[:, None, :] - scenarios[None, :, :], axis2) np.fill_diagonal(dist_matrix, np.inf) # 每个场景找最近邻居的距离 min_dists np.min(dist_matrix, axis1) min_idx np.argmin(min_dists) # 找到被删除场景最近的场景 nearest np.argmin(dist_matrix[min_idx]) # 概率转移 probs[nearest] probs[min_idx] # 删除场景 keep_mask np.ones(n, dtypebool) keep_mask[min_idx] False scenarios scenarios[keep_mask] probs probs[keep_mask] n - 1 return scenarios, probs # 先用K-means的结果做一次概率归一化 reduced_probs np.array(reduced_probs) reduced_probs reduced_probs / reduced_probs.sum() final_scenarios, final_probs synchronous_backward_reduction( np.array(reduced_scenarios), reduced_probs, target_count50 )这段代码有一个性能隐患每次循环都重新计算完整的距离矩阵所以在场景是300个、目标50个的时候要循环250次每次都要算300x300的距离矩阵体感会明显卡顿。如果你要处理上千个场景建议后面升级成预先算距离矩阵然后同步更新的版本或者直接用一些开源库的实现。4.3 画图对比削减前后的场景包络线削减完了一定要肉眼验证。我习惯画出削减前后的“场景包络带”也就是每个时段的5%分位数到95%分位数的范围。如果削减后的包络带和原始场景集的包络带基本吻合说明削减过程保留了主要的波动范围如果明显变窄说明你的保留场景丢掉了极端情况。def plot_scenario_band(ax, scenarios, probs, color, label): # 加权分位数 low [] high [] p_sorted np.sort(probs)[::-1] # 简化处理直接用等权分位数也够用加权更精确 for t_col in range(T): sorted_idx np.argsort(scenarios[:, t_col]) sorted_vals scenarios[sorted_idx, t_col] sorted_probs probs[sorted_idx] cum_probs np.cumsum(sorted_probs) low.append(sorted_vals[cum_probs 0.05][0]) high.append(sorted_vals[cum_probs 0.95][-1] if np.sum(cum_probs 0.95) 0 else sorted_vals[-1]) ax.plot(np.arange(T), low, colorcolor, linestyle--, labelf{label} 5%分位) ax.plot(np.arange(T), high, colorcolor, linestyle-., labelf{label} 95%分位) fig, ax plt.subplots(figsize(12, 6)) plot_scenario_band(ax, scenarios, probs, gray, 原始场景) plot_scenario_band(ax, final_scenarios, final_probs, red, 削减场景) ax.plot(np.arange(T), forecast, colorblack, linewidth2, label预测值) ax.legend() plt.savefig(scenario_reduction.png, dpi150, bbox_inchestight)加权分位数的计算比等权分位数要绕一点但这是必须的否则那些概率极小的极端场景和中概率场景被同等看待包络带会被极大值拉得很宽。4.4 定量评估削减前后统计特征对比除了肉眼看包络带还要算量化指标。我常用的有三类每个时段的均值误差、标准差误差、以及削减前后概率分布之间的Wasserstein距离。均值误差反映系统性的偏差标准差误差反映波动幅度的保留情况Wasserstein距离则从分布层面整体度量削减质量。# 对比均值与标准差 orig_mean np.average(scenarios, axis0, weightsprobs) red_mean np.average(final_scenarios, axis0, weightsfinal_probs) orig_std np.sqrt(np.average((scenarios - orig_mean) ** 2, axis0, weightsprobs)) red_std np.sqrt(np.average((final_scenarios - red_mean) ** 2, axis0, weightsfinal_probs)) mean_mae np.mean(np.abs(orig_mean - red_mean)) std_mae np.mean(np.abs(orig_std - red_std)) print(f削减前后均值绝对误差: {mean_mae:.4f} MW) print(f削减前后标准差绝对误差: {std_mae:.4f} MW)这个评估过程不能省。场景削减不是做完了就完事削减后的场景最终要作为随机优化模型的输入如果统计特征偏差太大优化的结果就会失去参考价值。至少要做均值误差和标准差误差的对比误差控制在单个时段额定容量的1%以内是比较理想的状态少于3%也可以接受超过5%就说明场景数压得太狠了。5. 常见问题与排查技巧实录5.1 模拟出来的场景里出现了负功率怎么办这是最常见的低级错误。用正态分布直接对功率误差抽样无条件截断就没法保证出力落在0到装机容量之间。解决办法是在抽样后加一步clip操作。但直接clip有个隐患大量负值被截到0会让出力等于0的场景概率堆积失真。更严谨的处理是拒绝采样凡是抽到范围外的样本直接丢弃重新抽一次。在样本量特别大的时候拒绝采样会额外增加计算时间一个折中方案是先构造截断正态分布用scipy的truncnorm直接生成合法的样本。from scipy.stats import truncnorm def truncnorm_sample(mean, std, low0, high200, size1): if std 0: return np.full(size, mean) a (low - mean) / std b (high - mean) / std return truncnorm.rvs(a, b, locmean, scalestd, sizesize)这样抽出来的样本天然限定在出力范围内不会因为大量截断导致分布形状改变。5.2 K-means结果不稳定每次都长得不一样K-means的初值敏感性是出了名的每次运行聚类中心都可能落到不同的局部最优解。解决办法有三个层次第一初始化参数从random改成k-means这是sklearn的默认值覆盖了绝大多数需求第二设置n_init参数让算法多次跑然后选最优的一次我习惯至少20次第三固定random_state保证结果可复现这对做研究和发论文来说是一个容易被忽略但很重要的细节。另外如果你在乎的是最终样本集而不仅是聚类标签可以对原始场景做一次标准化再聚类降低极端值对聚类结果的牵引。5.3 同步回代消除的概率归一化问题如果直接用每个场景概率为1/N的等概率假设回代过程中概率转移不会出问题最后概率总和保持为1。但如果原始场景集来自不同置信区间的组合各场景概率本身就不同这时候删除场景的概率累加就一定要做。我见过很多开源代码实现里只删场景不加概率最后得到的场景概率明显失真优化结果天然就偏了。判断概率处理是否正确的方法很简单削减前后分别算概率加权总和的期望值如果对得上说明概率守恒做得对。5.4 到底削减到多少个场景才合适这不是一个纯数学问题而是要看下游任务的需求。如果是做两阶段随机规划场景数K直接决定第二阶段的变量数你得好求解器能跑动。一般几十个场景在学术论文里最常见50到100个相对稳妥。如果只是做概率潮流或者做时序生产模拟需要的场景通常更少20到30个就够覆盖主要的运行场景。给一个经验法则先用原始场景集算一些关键统计量均值、方差、极端分位数然后从K10开始逐渐加大场景数每换一个K都去比对统计量误差。误差曲线在某一个K值后会进入平台期再多增加场景数也不能显著提高精度这个K就是最合适的选择。5.5 不要只用正态分布包打天下虽然这篇文章的标题是“基于正态分布”但在实操中还是要注意适用边界。正态分布对预测误差的建模效果不错但遇到极端天气事件、台风、寒潮这类气象过程误差分布会出现明显的厚尾现象仅仅用正态分布会低估极端场景的出现概率。这种情况下可以考虑混合分布比如把正态分布和广义帕累托分布组合起来或者用经验分布直接做核密度估计。场景削减的无缝之处就在于不管用什么分布生成场景削减算法本身都能直接复用。风光的时序相关性也是值得扩展的点。光伏和风电不仅在同一时刻相关还存在时间维度上的自相关性比如今天风速高明天风速也通常偏高。简单场景抽样生成的序列容易变成白噪音处理办法是在抽样环节加入时间相关性比如用AR模型先对误差做时序拟合再在这个基础上加噪声。说实话场景削减这套东西做了几次以后最大的感触是算法本身并不难难点都在数据特征和损失评估上。你花一小时搞清楚数据应该按几个时段拆、要不要做标准化、概率怎么守恒后面调算法就是水到渠成的事情。反过来如果你数据没准备好就急着跑步聚类、跑回代出来的结果花再多时间调参也很难满意。我现在的固定工作流是先做数据时段拆分和相关性分析再用截断正态分布生成几千个场景K-means粗减到几百个同步回代精修到目标数量最后每次都要画包络线和统计误差对比图做验收。这套流程跑顺了整套场景库从0到1大概半天就能搞定。下次再做类似的风光模拟背景从光伏换成风电、从单站点换到多站点只需要改第一段数据的逻辑后面的削减流程完全不用动。这大概是这类工作里最有价值的复利效应了。
风光场景模拟与削减实战：正态分布、K-means与同步回代消除
做风光模拟的时候很多人最开始的想法就是多抽几个场景几千个不够就抽一万个保证把不确定性都覆盖到。但真到了做随机优化、做生产模拟的时候一万个场景扔进去求解器直接罢工好不容易跑出个结果等得黄花菜都凉了。这时候就绕不开一个操作场景削减。说白了就是拿少数几个有代表性的场景去近似替代原来那成千上万个场景同时让概率分布特征尽量不丢。这篇文章我就围绕基于正态分布的风光模拟和场景削减把整个技术链路怎么搭、算法怎么选、代码怎么写、坑在哪里一次讲透适合正在做微电网调度、电力系统随机规划或者风光出力不确定性分析的朋友参考。1. 风光模拟与场景削减先把要解决什么问题讲清楚1.1 为什么需要场景不确定性建模到底在建什么风光出力的本质是一个随机过程。风速受气压、地形、昼夜温差影响光照受云层、季节、大气透明度影响这些因素叠加起来让光伏和风电在任意时刻的出力都带有强烈的不确定性。我们在做系统规划或者调度策略时不能只按“晴天满发”这种理想情况算否则实际运行一遇到多云天气功率预测偏差就能让你安排的机组组合全部失效。场景其实就是把这种不确定性离散化。把连续的概率分布变成一组离散的样本每个样本代表一种可能的风光出力时序。比如风电出力的预测值是30MW实际可能是25、28、35服从某种分布我们就按这个分布抽很多个可能的出力曲线每条曲线都是一个场景。理论上抽得越多覆盖的分布空间越完整但优化问题的规模也随之爆炸。1.2 为什么一定需要削减计算代价和质量之间的平衡我见过有人直接用5000个场景跑一个两阶段的随机优化结果两三天没跑完最后只能砍需求。场景规模对计算复杂度的影响是灾难性的因为随机规划问题里每个场景都对应一组决策变量和约束条件。假设你原来确定性模型的变量是1000个用500个场景做随机规划变量数直接膨胀到50万量级普通内存的机器根本扛不住。场景削减就是在这时候上场。它的数学本质是从原始的大规模场景集里选出一个子集再给每个保留场景重新分配概率使得削减前后场景集的概率分布差异最小。好的削减方法能用几十个场景保留几千个场景的主要统计特征均值、方差、相关性都能对得上。1.3 整个技术路线其实就三步我做风光模拟加场景削减的完整路线是固定的先做数据清洗和分布拟合再做蒙特卡洛抽样生成海量场景最后用削减算法挑代表场景。三者缺一不可。分布拟合不准后面抽出来的场景全是垃圾抽样数量不够削减完的代表性也打折扣削减算法不合适得到的场景集在优化问题里会和真实场景相差很大。所以整篇文章的推进思路就是按这个顺序来每一步我会把背后逻辑讲清楚再给出可以直接复制的代码和参数建议。2. 正态分布风光模拟的数学底座2.1 为什么是正态分布不是硬套是有依据的先泼一盆冷水风光出力的原始分布并不是严格的正态分布。风速一般更接近威布尔分布光照辐照度在晴天时接近Beta分布。那为什么标题还叫“基于正态分布的风光模拟”因为在大多数实际工程场景里我们并不直接对出力值本身建模而是对预测误差建模而预测误差在中心极限定理的意义下往往可以合理近似为正态分布。比如风电功率预测系统给出明天某个时刻出力是50MW实际的出力等于预测值加上一个随机误差项。这个误差受很多独立随机因素影响按中心极限定理大量独立随机变量之和趋向正态分布。所以我们建模时写成P_actual P_forecast epsilon, epsilon ~ N(0, sigma^2)注意功率本身的正负是有边界的0到装机容量之间所以对误差做正态抽样之后必须做截断处理把超出边界的情况截掉或者重新采样。这个细节很多初学者会忽略导致模拟出的场景出现负功率一看就是假的。从数据拟合的角度正态分布的好处是参数少、含义直观就均值和标准差两个量工程上非常容易估计。你只需要统计历史预测误差的均值和方差就可以搭起整个模拟框架。2.2 参数估计不能光看整体要按时段拆有些人拿到一整年的历史数据哗啦啦全丢进去算出均值和方差就开始做模拟。这样做出来的东西基本上是废的。因为风光预测误差有明显的时段特征白天光伏出力大误差绝对值和方差都大凌晨光伏不出力误差方差趋近于零。如果混在一起拟合相当于用一个大方差去覆盖所有时段白天的场景过于保守夜间的场景又过于激进。正确做法是按时间段拆开统计。以1小时间隔为例一天24小时就分24个时段或者按季节再加一层拆。对每个时段分别估计误差的均值 mu_t 和标准差 sigma_t。实际操作中可以写成函数对DataFrame按小时分组计算然后做平滑处理避免某些时段样本量太少导致方差畸高或畸低。我一般是用指数加权移动平均做平滑让相近时段的标准差不会突变。这样拟合出来的参数更能反映实际的物理过程光靠白天数据多就把标准差拉大是不合理的。2.3 相关性是个隐秘的坑单独对风、对光分别做正态抽样会忽略风光之间的相关性。比如同一片区域内风电偏大往往对应阴天或大风天气光伏输出可能下降反过来强日照天气光伏大发时风速往往偏低。如果模拟时光是独立抽样就可能出现“又大风又大太阳”的场景这在物理上很罕见但在数学上完全可能发生。处理办法有两个思路。一个是用Copula函数对风光的联合分布建模另一个更常用的是在抽样时引入相关系数矩阵通过Cholesky分解生成相关正态随机数。Cholesky分解的思路是把相关的多维正态分布拆成下三角矩阵乘法对一个独立的标准正态随机向量做线性变换就得到期望的相关结构。这段逻辑在代码里也就几行的事但效果非常显著能直接把模拟场景的物理合理性拉上一个台阶。3. 场景削减的核心算法从K-means到同步回代消除3.1 场景相似性的度量削减的第一步其实是“算距离”不管是哪种场景削减算法都绕不开一个问题怎么判断两个场景像不像。最常用的就是欧氏距离加L2范数。假设每个场景是T个时段的一条曲线场景i和场景j之间的距离就可以定义成所有时段功率差值的平方和再开根号。D(i, j) sqrt( sum_{t1}^{T} (P_i(t) - P_j(t))^2 )这个距离定义简洁、计算高效对时序曲线的整体形状差异很敏感。不过它也有缺陷如果某个时段的误差特别大平方项会把距离拉得很大掩盖其他时段的相似性。所以有的场景集在做距离计算之前会先对数据进行标准化把每个时段的出力除以该时段的容量或者历史标准差相当于给不同时段加了权重避免那些绝对数值大的时段统治距离计算。如果场景携带概率权重比如每个场景本身有不同的发生概率p_i那么在算法迭代中还要用到加权距离。加权距离的计算不是简单乘一下而是要考虑场景代价的累积这个在同步回代消除里体现得尤为明显。3.2 方法一K-means聚类场景削减K-means可能是上手最快、工程里最喜欢用的场景削减方法。核心思路很朴素把N个原始场景通过聚类压缩成K个簇每个簇的中心就是保留场景簇内所有场景的概率之和就是保留场景的新概率。具体步骤分四步从原始场景中随机选择K个场景作为初始聚类中心对每个场景计算到K个中心的距离把它归到距离最近的中心所在的簇对每个簇内所有场景取平均得到新的中心重复第2和第3步直到聚类中心不再变化或者达到最大迭代次数。这里有一个常规文档里不会细说的重点聚类完成后中心曲线是一个平均曲线它并不一定等于原始场景中的某一条。但在场景削减的实际使用场景里我们通常希望保留场景是可回溯的、和原始数据同构的曲线。所以更稳妥的做法是保留每个簇中距离簇中心最近的那个原始场景作为代表而不是直接使用平均中心。K值的选取也是一个权衡。K太小场景多样性不足极端场景全被抹掉K太大削减又没达到降低计算复杂度的目的。我一般先用轮廓系数圈定一个大致范围再结合随机优化问题的具体规模反推K的上限。比如决策变量总数不能超过求解器能处理的量级在这个上限以内尽量选大的K。3.3 方法二同步回代消除后向消除同步回代消除是文献里出现频率很高的另一种削减算法它的思路和K-means完全不同不是聚类求中心而是从原始场景集里不断“删”场景一直删到剩下目标数量为止。每一步都删除一个对整体概率分布影响最小的场景并把被删除场景的概率加到离它最近的场景上去。算法步骤我整理了一下计算原始场景集中所有场景两两之间的距离形成距离矩阵对每个场景i找到距离它最近的场景j这个最近距离记为d_i在所有场景中找d_i最小的那个场景它的意义是删掉它对整体分布扰动最小把这个场景的概率累加到它最近的那个场景上然后从场景集中移除它更新距离矩阵重复步骤2到4直到剩余场景数等于目标K。单看文字有点干但代码实现起来并不复杂。它的一个优势是每一步都在做全局最优的局部决策得到的保留场景一定是原始场景里的真实曲线物理上可解释性非常强。缺点是在场景数量大的时候计算量偏高因为每删一个场景都要更新一次距离矩阵时间复杂度接近O(N^3)量级。我在实际应用里一般把两种方法配合使用先用K-means粗筛从几千个场景快速压到几百个再用同步回代消除精修从几百个压到最终需要的几十个。两步走既能利用K-means的速度优势又能利用回代消除的精度优势。3.4 两种方法的对比与选择K-means和同步回代消除没有绝对的好坏关键看应用场景。K-means的优势是速度快、实现简单、处理超大规模场景集很轻松不足是结果受初始中心影响有随机性而且平均中心曲线会抹平极端情况。同步回代消除的结果更稳定、更能保留原始场景形态但计算复杂度高场景数和维度一旦上去耗时明显增加。如果你的原始场景集只有几百个直接用同步回代消除完全没问题。如果你蒙特卡洛抽了一万个场景建议别直接跑回代先K-means粗降到300到500个再回代精修到目标数量。这个组合拳在很多论文里都有影子工程上确实好用。4. 完整实操Python实现风光模拟与场景削减4.1 准备模拟数据先造一组实验数据演示完整流程。假设有一个48时段的风电出力预测值序列真实出力等于预测出力加正态误差标准差设置为预测值的10%左右同时做0到装机容量的截断。这里的真实出力是我们用来验证的基准实际应用中没有真实出力只有预测值和误差分布。import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.stats import norm from sklearn.cluster import KMeans np.random.seed(42) # 生成风电预测序列48时段容量200MW T 48 capacity 200 t np.arange(T) # 模拟一个峰谷明显的预测曲线 forecast 120 60 * np.sin(t / T * 2 * np.pi - np.pi/2) forecast np.clip(forecast, 20, 180) # 误差标准差为预测值的10% sigma 0.10 * forecast # 生成5000个随机场景 n_scenarios 5000 errors np.random.normal(loc0, scale1, size(n_scenarios, T)) * sigma scenarios forecast errors scenarios np.clip(scenarios, 0, capacity) # 每个场景等概率 probs np.ones(n_scenarios) / n_scenarios注意代码里的一个细节我先用标准正态随机数乘以逐时段的sigma而不是一次性np.random.normal传入整个sigma矩阵。两种写法结果等价但前者更直观方便你以后扩展时对每个时段用不同的分布类型。4.2 从K-means到精修场景集从5000压到50先做K-means快速粗削减到300个场景再用同步回代消除把300个场景精修到50个。K-means部分可以直接调sklearn重点在于聚类完成后按簇内最近原始场景提取保留场景。# K-means粗削减 n_clusters 300 kmeans KMeans(n_clustersn_clusters, random_state42, n_init20) labels kmeans.fit_predict(scenarios) # 每个簇内找离中心最近的原始场景作为保留场景 reduced_scenarios [] reduced_probs [] for c in range(n_clusters): idx np.where(labels c)[0] center kmeans.cluster_centers_[c] dists np.linalg.norm(scenarios[idx] - center, axis1) rep_idx idx[np.argmin(dists)] reduced_scenarios.append(scenarios[rep_idx]) reduced_probs.append(np.sum(probs[idx]))同步回代消除的实现里最核心的就是距离矩阵的维护和删除操作。网上能搜到的版本里很多没有注意“被删除场景概率要加到最近场景”这个细节导致削减后概率总和不等于1后面做优化就全乱了。我提供一个验证过概率守恒的版本def synchronous_backward_reduction(scenarios, probs, target_count): scenarios scenarios.copy() probs probs.copy() n scenarios.shape[0] while n target_count: # 两两场景距离 dist_matrix np.linalg.norm(scenarios[:, None, :] - scenarios[None, :, :], axis2) np.fill_diagonal(dist_matrix, np.inf) # 每个场景找最近邻居的距离 min_dists np.min(dist_matrix, axis1) min_idx np.argmin(min_dists) # 找到被删除场景最近的场景 nearest np.argmin(dist_matrix[min_idx]) # 概率转移 probs[nearest] probs[min_idx] # 删除场景 keep_mask np.ones(n, dtypebool) keep_mask[min_idx] False scenarios scenarios[keep_mask] probs probs[keep_mask] n - 1 return scenarios, probs # 先用K-means的结果做一次概率归一化 reduced_probs np.array(reduced_probs) reduced_probs reduced_probs / reduced_probs.sum() final_scenarios, final_probs synchronous_backward_reduction( np.array(reduced_scenarios), reduced_probs, target_count50 )这段代码有一个性能隐患每次循环都重新计算完整的距离矩阵所以在场景是300个、目标50个的时候要循环250次每次都要算300x300的距离矩阵体感会明显卡顿。如果你要处理上千个场景建议后面升级成预先算距离矩阵然后同步更新的版本或者直接用一些开源库的实现。4.3 画图对比削减前后的场景包络线削减完了一定要肉眼验证。我习惯画出削减前后的“场景包络带”也就是每个时段的5%分位数到95%分位数的范围。如果削减后的包络带和原始场景集的包络带基本吻合说明削减过程保留了主要的波动范围如果明显变窄说明你的保留场景丢掉了极端情况。def plot_scenario_band(ax, scenarios, probs, color, label): # 加权分位数 low [] high [] p_sorted np.sort(probs)[::-1] # 简化处理直接用等权分位数也够用加权更精确 for t_col in range(T): sorted_idx np.argsort(scenarios[:, t_col]) sorted_vals scenarios[sorted_idx, t_col] sorted_probs probs[sorted_idx] cum_probs np.cumsum(sorted_probs) low.append(sorted_vals[cum_probs 0.05][0]) high.append(sorted_vals[cum_probs 0.95][-1] if np.sum(cum_probs 0.95) 0 else sorted_vals[-1]) ax.plot(np.arange(T), low, colorcolor, linestyle--, labelf{label} 5%分位) ax.plot(np.arange(T), high, colorcolor, linestyle-., labelf{label} 95%分位) fig, ax plt.subplots(figsize(12, 6)) plot_scenario_band(ax, scenarios, probs, gray, 原始场景) plot_scenario_band(ax, final_scenarios, final_probs, red, 削减场景) ax.plot(np.arange(T), forecast, colorblack, linewidth2, label预测值) ax.legend() plt.savefig(scenario_reduction.png, dpi150, bbox_inchestight)加权分位数的计算比等权分位数要绕一点但这是必须的否则那些概率极小的极端场景和中概率场景被同等看待包络带会被极大值拉得很宽。4.4 定量评估削减前后统计特征对比除了肉眼看包络带还要算量化指标。我常用的有三类每个时段的均值误差、标准差误差、以及削减前后概率分布之间的Wasserstein距离。均值误差反映系统性的偏差标准差误差反映波动幅度的保留情况Wasserstein距离则从分布层面整体度量削减质量。# 对比均值与标准差 orig_mean np.average(scenarios, axis0, weightsprobs) red_mean np.average(final_scenarios, axis0, weightsfinal_probs) orig_std np.sqrt(np.average((scenarios - orig_mean) ** 2, axis0, weightsprobs)) red_std np.sqrt(np.average((final_scenarios - red_mean) ** 2, axis0, weightsfinal_probs)) mean_mae np.mean(np.abs(orig_mean - red_mean)) std_mae np.mean(np.abs(orig_std - red_std)) print(f削减前后均值绝对误差: {mean_mae:.4f} MW) print(f削减前后标准差绝对误差: {std_mae:.4f} MW)这个评估过程不能省。场景削减不是做完了就完事削减后的场景最终要作为随机优化模型的输入如果统计特征偏差太大优化的结果就会失去参考价值。至少要做均值误差和标准差误差的对比误差控制在单个时段额定容量的1%以内是比较理想的状态少于3%也可以接受超过5%就说明场景数压得太狠了。5. 常见问题与排查技巧实录5.1 模拟出来的场景里出现了负功率怎么办这是最常见的低级错误。用正态分布直接对功率误差抽样无条件截断就没法保证出力落在0到装机容量之间。解决办法是在抽样后加一步clip操作。但直接clip有个隐患大量负值被截到0会让出力等于0的场景概率堆积失真。更严谨的处理是拒绝采样凡是抽到范围外的样本直接丢弃重新抽一次。在样本量特别大的时候拒绝采样会额外增加计算时间一个折中方案是先构造截断正态分布用scipy的truncnorm直接生成合法的样本。from scipy.stats import truncnorm def truncnorm_sample(mean, std, low0, high200, size1): if std 0: return np.full(size, mean) a (low - mean) / std b (high - mean) / std return truncnorm.rvs(a, b, locmean, scalestd, sizesize)这样抽出来的样本天然限定在出力范围内不会因为大量截断导致分布形状改变。5.2 K-means结果不稳定每次都长得不一样K-means的初值敏感性是出了名的每次运行聚类中心都可能落到不同的局部最优解。解决办法有三个层次第一初始化参数从random改成k-means这是sklearn的默认值覆盖了绝大多数需求第二设置n_init参数让算法多次跑然后选最优的一次我习惯至少20次第三固定random_state保证结果可复现这对做研究和发论文来说是一个容易被忽略但很重要的细节。另外如果你在乎的是最终样本集而不仅是聚类标签可以对原始场景做一次标准化再聚类降低极端值对聚类结果的牵引。5.3 同步回代消除的概率归一化问题如果直接用每个场景概率为1/N的等概率假设回代过程中概率转移不会出问题最后概率总和保持为1。但如果原始场景集来自不同置信区间的组合各场景概率本身就不同这时候删除场景的概率累加就一定要做。我见过很多开源代码实现里只删场景不加概率最后得到的场景概率明显失真优化结果天然就偏了。判断概率处理是否正确的方法很简单削减前后分别算概率加权总和的期望值如果对得上说明概率守恒做得对。5.4 到底削减到多少个场景才合适这不是一个纯数学问题而是要看下游任务的需求。如果是做两阶段随机规划场景数K直接决定第二阶段的变量数你得好求解器能跑动。一般几十个场景在学术论文里最常见50到100个相对稳妥。如果只是做概率潮流或者做时序生产模拟需要的场景通常更少20到30个就够覆盖主要的运行场景。给一个经验法则先用原始场景集算一些关键统计量均值、方差、极端分位数然后从K10开始逐渐加大场景数每换一个K都去比对统计量误差。误差曲线在某一个K值后会进入平台期再多增加场景数也不能显著提高精度这个K就是最合适的选择。5.5 不要只用正态分布包打天下虽然这篇文章的标题是“基于正态分布”但在实操中还是要注意适用边界。正态分布对预测误差的建模效果不错但遇到极端天气事件、台风、寒潮这类气象过程误差分布会出现明显的厚尾现象仅仅用正态分布会低估极端场景的出现概率。这种情况下可以考虑混合分布比如把正态分布和广义帕累托分布组合起来或者用经验分布直接做核密度估计。场景削减的无缝之处就在于不管用什么分布生成场景削减算法本身都能直接复用。风光的时序相关性也是值得扩展的点。光伏和风电不仅在同一时刻相关还存在时间维度上的自相关性比如今天风速高明天风速也通常偏高。简单场景抽样生成的序列容易变成白噪音处理办法是在抽样环节加入时间相关性比如用AR模型先对误差做时序拟合再在这个基础上加噪声。说实话场景削减这套东西做了几次以后最大的感触是算法本身并不难难点都在数据特征和损失评估上。你花一小时搞清楚数据应该按几个时段拆、要不要做标准化、概率怎么守恒后面调算法就是水到渠成的事情。反过来如果你数据没准备好就急着跑步聚类、跑回代出来的结果花再多时间调参也很难满意。我现在的固定工作流是先做数据时段拆分和相关性分析再用截断正态分布生成几千个场景K-means粗减到几百个同步回代精修到目标数量最后每次都要画包络线和统计误差对比图做验收。这套流程跑顺了整套场景库从0到1大概半天就能搞定。下次再做类似的风光模拟背景从光伏换成风电、从单站点换到多站点只需要改第一段数据的逻辑后面的削减流程完全不用动。这大概是这类工作里最有价值的复利效应了。