简介基于拉丁超立方抽样与样本削减的场景分析方案面向从事可再生能源出力预测、负荷预测及电力系统运行研究的工程师与科研人员。压缩包为rar格式仅含1个MATLAB脚本.m文件体积约1KB聚焦核心算法实现便于直接阅读与复用。已有2593人学习该资源。脚本围绕风光出力及负荷场景的生成与分析展开先通过拉丁超立方抽样对风速、辐照度、温度等气象条件进行多维均匀采样构建代表性强的不确定性输入样本再结合样本削减策略从大规模候选场景中筛选典型子集降低计算开销而不牺牲预测精度。同时覆盖数据预处理、异常值清洗、缺失值填充及标准化等环节并配有清晰注释适合需要快速上手场景缩减方法或将其嵌入预测模型的读者。该方案兼顾多因素复杂性与计算效率对提升风光并网稳定性、辅助电网供需平衡评估具有实践参考价值。1. 拉丁超立方抽样为什么成了风光场景生成的首选从蒙特卡洛跑不动说起做微电网储能容量配置时最折磨人的不是设备选型而是风光出力说不清风不吹、光不照、负荷还天天变。拿历史一年数据直接跑场景太少等于看运气上蒙特卡洛抽一万个场景每次优化都像重新蹲一个算例算力成本直接失控。拉丁超立方抽样LHS用分层采样的方式把样本量压到十分之一样本削减再把上千个候选场景浓缩成十几个代表性典型场景——这就是「基于拉丁超立方抽样和样本削减的风光出力及负荷场景预测分析」这条线的核心价值。这篇笔记按我实际做过的方案讲分布怎么建、LHS 怎么调、削减怎么不丢极端信息、削减后怎么算出能说服人的置信区间给做微电网和电网规划的工程师当一份抄作业清单。2. 拉丁超立方抽样的落地实现从数据拟合到相关性控制2.1 分布怎么选Weibull、Beta 和正态的拟合细节场景生成的第一步不是抽样而是把历史数据变成可采样的分布。风速我一般用双参数 Weibull固定位置参数 loc0只拟合形状参数 c 和尺度参数 scale光伏辐照度先除以交流侧上限归一化到 0~1再用 Beta 分布拟合负荷直接用正态分布虽然尾部不够厚但对大多数规划研究够用了。选这三个分布不是因为精度最高而是因为它们跟物理约束对得上Weibull 能表达大风这种右偏尾部Beta 自带 0~1 边界正态对负荷这种大量样本的均值波动很稳。import numpy as np from scipy import stats # 历史数据: wind_speed(m/s), irradiance(W/m^2), load(kW) 为三个一维数组 c_w, _, s_w stats.weibull_min.fit(wind_speed, floc0) # 风速: 固定位置项 beta_a, beta_b, _, _ stats.beta.fit(irradiance / 1000.0, floc0, fscale1) mu_l, sigma_l stats.norm.fit(load) n_samples 2000 # 初始候选场景数常见取 2000~5000 rng np.random.default_rng(42) def lhs_from_dist(dist_ppf, n_samples, seed): rng np.random.default_rng(seed) u (np.arange(n_samples) rng.random(n_samples)) / n_samples rng.shuffle(u) return dist_ppf(u) s_wind lhs_from_dist(lambda u: stats.weibull_min.ppf(u, c_w, 0, s_w), n_samples, 0) s_solar lhs_from_dist(lambda u: stats.beta.ppf(u, beta_a, beta_b, 0, 1) * 1000.0, n_samples, 1) s_load lhs_from_dist(lambda u: stats.norm.ppf(u, mu_l, sigma_l), n_samples, 2)这段代码的核心在lhs_from_dist把 [0,1] 区间切成 n_samples 个等宽分层每层随机取一个分位点打乱后过逆 CDF 映射回物理量。LHS 保证每个变量的覆盖范围被强制铺满不像纯随机抽样那样可能整段尾部没采到。参数上floc0是风速拟合里最关键的约束不固定位置项可能拟合出负的起始风速后续采样全是残次品Beta 拟合必须先归一化采样后乘回 1000 W/m² 才能得到真实辐照度n_samples我一般取 2000 起步少于 1000 时尾部极端出力很容易在后续削减中彻底消失。2.2 Iman-Conover 重排把独立 LHS 变成有相关性的样本逐维独立 LHS 是很多新手直接踩坑的地方风速、辐照、负荷各抽各的结果散点图上一片混乱风大的日子光也强完全不合物理规律。现实中同一地区的风电和光伏通常是负相关午后辐照峰值时段风速往往偏低相邻风电场之间有空间正相关负荷与温度、季节有弱相关。要解决这个问题常见做法是先用 LHS 生成独立样本再用 Iman-Conover 方法重排各列的秩把目标秩相关注入进去。def iman_conover(X, target_corr): n, d X.shape # 把原始数据的每列变成秩再映射到标准正态空间 ranks np.argsort(np.argsort(X, axis0), axis0) 1 S stats.norm.ppf(ranks / (n 1)) # 当前 S 的秩相关矩阵及其 Cholesky 分解 corr_S np.corrcoef(S, rowvarFalse) L_S np.linalg.cholesky(corr_S 1e-10 * np.eye(d)) # 目标相关矩阵的 Cholesky 分解 L_T np.linalg.cholesky(target_corr 1e-10 * np.eye(d)) # 消除原相关并注入目标相关 Z (S np.linalg.inv(L_S).T L_T.T) # 按 Z 的列秩重排原始 X保留原始分布但改变秩相关 X_out np.zeros_like(X) for j in range(d): order np.argsort(np.argsort(Z[:, j])) X_out[:, j] np.sort(X[:, j])[order] return X_out target_corr np.array([ [1.0, -0.35, 0.10], [-0.35, 1.0, 0.05], [0.10, 0.05, 1.0] ]) X np.column_stack([s_wind, s_solar, s_load]) X_corr iman_conover(X, target_corr)这里必须注意Iman-Conover 控制的是秩相关而不是 Pearson 相关因为风速、辐照这些变量不服从联合正态Pearson 相关会被非线性变换扭曲。target_corr的数值建议直接从历史数据里算秩相关系数得到比如用scipy.stats.spearmanr如果没有完整历史数据按经验填也行风电-光伏负相关取 -0.3 到 -0.4 是合理的。Cholesky 分解要求相关矩阵半正定手工填的矩阵很容易特征值小于 0代码里加1e-10 * np.eye(d)是防病态真要碰上强负相关组合还得先对矩阵做特征值截断修正。2.3 边界截断和初筛先看截断比例再谈削减LHS 采样完不能直接拿去削减得先过一遍物理边界。Weibull 拟合偶尔会采出负风速Beta 采出的辐照度可能超过光伏逆变器的交流侧上限负荷也可能被正态分布尾部带成负值。常见做法是裁剪后统计截断比例这个比例是判断分布拟合质量的重要信号。wind_cap, solar_cap 60.0, 30.0 # MW按实际电站装机填 s_wind_c np.clip(s_wind, 0, None) s_solar_c np.clip(s_solar, 0, solar_cap) s_load_c np.clip(s_load, 0, None) print(wind 截断比例:, (s_wind_c ! s_wind).mean()) print(solar 超限比例:, (s_solar_c ! s_solar).mean())截断比例在 3% 以下可以接受超过 5% 就说明分布假设或拟合参数有问题而不是简单裁一刀了事。比如风速截断高多半是 Weibull 形状参数偏小、分布右尾太重辐照超限多往往是归一化时用的分母不是逆变器交流侧上限而是直流侧峰值导致 Beta 分布位置整体偏移。做完截断我还会画一张风速-辐照散点图跟历史数据对比如果 LHS 样本的相关性形态和实际明显不同说明 target_corr 该调。初筛这一步不做后面所有场景分析都会带着系统性偏差。3. 样本削减的两种做法同步回代与 K-means 的关键参数3.1 同步回代削减的步骤与复杂度保留真实场景的路线场景削减的目标是去掉冗余、保留代表性。文献里最常见的 scenario reduction 是同步回代削减每次找到距离最近的一对场景把概率较小的那个删掉并把它的概率叠加到另一个上循环到只剩目标数量。这个思路来自 Heitsch 和 Römisch 的 Kantorovich 距离近似好处是留下的都是真实采样点极端出力日期有可能被保留坏处是计算复杂度高2000 个场景起步时跑起来很慢。def backward_reduction(scenes, probs, n_target): idx list(range(len(scenes))) p probs.copy() while len(idx) n_target: d np.sqrt(((scenes[idx][:, None, :] - scenes[idx][None, :, :]) ** 2).sum(axis-1)) np.fill_diagonal(d, np.inf) i, j np.unravel_index(np.argmin(d), d.shape) # 概率小的场景被削减概率并给概率大的场景 if p[idx[i]] p[idx[j]]: p[idx[i]] p[idx[j]] dell j else: p[idx[j]] p[idx[i]] dell i del idx[dell] return scenes[idx], p[idx] / p[idx].sum()这个实现的每次循环都重算完整距离矩阵2000 个初始场景要循环到 10 个场景时耗时可能到分钟级不是不能等但没必要。正式项目里我会先用 K-means 粗聚类定下目标场景数再用同步回代做最终削减。参数上最核心的是n_target一般取 5~20太少了极端概率事件全丢太多了削减失去意义。还要强调一个原则距离计算必须用整个时序向量做欧氏距离不能按时刻逐点独立削减否则会破坏风速-辐照的时序相关性。3.2 K-means 削减聚类中心与概率赋值K-means 是更工程化的选择速度快、评估方便、参数直觉。做法是把每个场景当成一条 T 维时序样本聚类中心就是典型场景每类样本占总样本的比例就是该场景的概率。from sklearn.cluster import KMeans def kmeans_reduction(scenes, n_target, seed42): km KMeans(n_clustersn_target, n_init50, random_stateseed) labels km.fit_predict(scenes) centers km.cluster_centers_ count np.bincount(labels, minlengthn_target) probs count / count.sum() return centers, probs, labelsn_init50是我从翻车里换来的习惯默认值 10 在场景集有重叠簇时很容易陷进局部最优导致两次跑出的典型场景差异很大。random_state固定是为了让报告里的结果可复现评审问起「为什么场景长这样」时能拿出确定性的图。K-means 输出的是聚类中心本质是对场景做平滑这会直接抹掉极端峰值。所以它适合负荷场景这种波动相对平滑的序列风光场景想保留大风、强辐照的极端日我会在聚类之后把原始样本里超过 90% 分位的极端场景单独挑出来强制塞进典型集合并重新归一遍概率。3.3 ADE 与轮廓系数削减质量的量化标准削减得好不好不能靠肉眼两个量化指标我每次必算ADE 衡量削减前后概率分布的期望偏差轮廓系数衡量聚类内聚与分离程度。def ade(orig_scenes, orig_probs, red_scenes, red_probs): orig_mean (orig_probs.reshape(-1, 1) * orig_scenes).sum(axis0) red_mean (red_probs.reshape(-1, 1) * red_scenes).sum(axis0) return np.linalg.norm(orig_mean - red_mean) / (np.linalg.norm(orig_mean) 1e-12) from sklearn.metrics import silhouette_score def quick_check(labels, scenes): print(轮廓系数:, silhouette_score(scenes, labels))ADE 的分子是削减前后概率加权期望向量的欧氏距离分母用于归一化结果在 3% 以下我认为可以接受超过 5% 必须增加场景数或换削减方法。但只盯均值会骗自己均值几乎不会被削减影响真正的损失在分布尾部。我一般会额外比较削减前后每小时出力的 10% 和 90% 分位数如果分位数偏差超过 5%说明极端场景已经被吃掉要用第 5 章的方法补救。轮廓系数只在 K-means 时有意义同步回代没有这个指标。注意ADE 小不等于削减成功分位数对比才是兜底检查。均值是抛物线的最低点谁都能碰上尾部才是风险所在。4. 风光出力及负荷场景预测分析从典型场景到运行评估4.1 从典型场景到加权统计量期望、分位数、极端值削减完成后十个典型场景进入分析阶段。这里的「预测分析」不是指预测明天几点几分出多少电而是给出一个完整的不确定性描述最可能的出力曲线、80% 置信带、极端场景的出现概率。第一步是把削减后的概率权重用起来算加权期望和加权分位数。probs red_probs weighted_mean (probs.reshape(-1, 1) * red_scenes).sum(axis0) def weighted_quantile(scenes, probs, q): order np.argsort(scenes, axis0) wq np.empty(scenes.shape[1]) for t in range(scenes.shape[1]): cum_w np.cumsum(probs[order[:, t]]) wq[t] np.interp(q, cum_w / cum_w[-1], scenes[order[:, t], t]) return wq q10 weighted_quantile(red_scenes, probs, 0.1) q90 weighted_quantile(red_scenes, probs, 0.9)weighted_quantile是按概率权重累积后插值比直接用np.quantile可靠得多。砍到 10 个场景后概率并不均等如果忽略权重直接取 10 个场景的 10% 分位等于默认每个场景概率相同小概率场景的影响力被放大。q10 和 q90 构成 80% 置信带q10 适合做保守出清和储能备电的下界q90 适合做线路不越限的上界。若发现置信带宽到没法用原因多半不是削减数量不够而是初始分布拟合太宽回到第 2 章查分布参数。4.2 概率赋值与置信区间为什么不能用简单平均削减后的场景概率赋值有讲究K-means 的概率就是簇内样本数除以总样本数同步回代的概率来自概率合并过程。但如果强行保留极端场景比如手动加入大风日、阴雨日就需要把极端场景的概率和原典型场景概率一起重新归一化否则所有场景概率加起来不等于 1。# 假设 extreme_scenes 是按规则挑出的 2 个极端场景各占原始样本的 0.5% extreme_idx np.array([True, True]) extreme_probs np.array([0.005, 0.005]) all_scenes np.vstack([red_scenes, extreme_scenes]) all_probs np.concatenate([red_probs * (1 - extreme_probs.sum()), extreme_probs]) all_probs / all_probs.sum() # 重新归一化这样处理后极端场景参与置信区间计算但其概率权重被压缩到真实水平不会把 80% 置信带硬生生撑宽。实际算置信区间时我经常用 Bootstrap 从削减后场景集里重复抽样 1000 次统计分位数的稳定性如果分位数的标准差超过 5%说明场景集还不够支撑结论。4.3 下游运行评估场景集怎么喂给优化和仿真削出来的典型场景最终要进优化模型或时序仿真。这里最常见的错误是把每个典型场景单独跑一次、最后取算术平均这忽略了场景概率的差异。正确的做法是为每个场景算指标后做概率加权。metrics np.zeros(n_target) for s in range(n_target): # run_operation 是任意的运行模拟函数返回成本/容量/越限次数等指标 metrics[s] run_operation(red_scenes[s], ...) expected_metric (red_probs * metrics).sum()run_operation可以是储能容量配置的线性规划、机组组合、或者一次简单的潮流计算具体函数各项目不同但加权逻辑是通用的。把 expected_metric 和单场景的 max/min 一起写进报告比只报期望值更有信息量。实际项目里我还会额外输出「最差场景指标」也就是概率加权后表现最差的 10% 场景对应的结果用这个值指导安全裕度设置。5. 场景削减的常见问题与排查五个高发坑的根因与对策5.1 削减后丢了极端日为什么典型场景比原始数据还保守现象削减后场景集里最大风电出力只有原始历史的 70%储能配置结果偏小后续安全校核过不去。原因是 K-means 以欧氏距离为准则离群点对聚类中心影响极小极端大风日被平均成了普通场景同步回代如果概率相差悬殊小概率极端场景会先被合并掉。解决方法是先聚类再从原始样本中挑出超过 90% 分位的极端场景强制保留重新归一化概率。我一般把极端日分位数阈值写到配置里每次削减自动执行。5.2 LHS 采出负风速和超装机的出力现象采样结果导入优化模型报负风速辐照值超过逆变器容量模型直接不可行。原因是 Weibull 拟合时没固定 loc或者 Beta 归一化分子分母搞混另一种是采样后没做截断就进了削减。解决分两步拟合时风速固定floc0辐照归一化用交流侧上限采样后统一 clip 并打印截断比例。截断比例超过 5% 不是裁一刀的问题是分布拟合错了回头检查参数。5.3 相关性控制后散点图还是不对现象做了 Iman-Conover 重排后风速-辐照散点图仍然看不出负相关或者报了 LinAlgError。原因是目标相关矩阵填的是 Pearson 相关而数据强非高斯Pearson 相关在非线性变换下根本不守恒另外手工填的矩阵特征值可能为负Cholesky 分解直接失败。解决方法是改传spearmanr算出的秩相关矩阵分解前对特征值做修正给矩阵对角线加一个小的正则项代码里那行1e-10 * np.eye(d)就是干这个的。重排后重新计算秩相关来验证而不是看 Pearson 相关。5.4 置信区间宽到没有参考价值现象80% 置信带从 0 一直铺到满发调度拿到手根本没法用。原因是场景数太少又直接用极值当上下界或者极端场景的概率权重在削减时不正常放大。解决是用加权分位数而不是 max/min想要更稳就引入核密度估计或 Bootstrap 平滑。多个项目做下来置信带过宽基本都是分布拟合阶段方差参数没调好削减只是背锅。5.5 被问「为什么选 10 个场景」时答不上来现象评审或领导指着报告问10 个典型场景的依据是什么现场答不出可量化的依据。原因是只按经验取数没做收敛性验证。解决是提前跑一组 ADE-场景数曲线分别削减到 5、10、15、20、30 个场景算对应 ADE选 ADE 下降明显变缓的点作为最终场景数。把这张曲线图直接放进报告附录比任何口头解释都有说服力。6. 削减场景做回测验证的三个实用技巧验证削减结果靠不靠谱我的习惯是做三层回测每个都不复杂但能堵住大部分质疑。第一层是历史留出法。取一年历史数据用前 11 个月拟合分布并生成场景拿最后一个月真实出力作为检验样本统计实际值落在削减场景 80% 置信带里的比例。这个比例应该接近 80%如果只有 50%说明置信带整体偏窄分布参数的方差估计小了。from scipy.stats import ks_2samp, wasserstein_distance # hist_flatten 和 red_flatten 分别是历史观测和削减场景展平的一维数组 ks_stat, p_value ks_2samp(hist_flatten, red_flatten) wd wasserstein_distance(hist_flatten, red_flatten) print(fKS p{p_value:.3f}, Wasserstein{wd:.3f})第二层是分布对比。KS 检验看削减场景与历史的累计分布是否一致Wasserstein 距离量化两个分布之间的搬运成本。p 值小于 0.05 说明削减后场景集和原始数据分布差异显著这时要回头调削减数或分布参数。第三层是决策结果收敛性验证用 5、10、20 个场景分别跑同一个储能配置问题看最优容量变化是否在可接受区间内。如果 10 个和 20 个场景的结果差得很大说明 10 个场景承载不了这个问题的不确定性。我自己翻车最多的就是第二层检验栽在 K-means 聚类中心上聚类中心是平滑点一对比分布就发现尾部全没了。后来养成习惯削减完先画典型场景曲线图肉眼扫一遍再往下走。场景削减不是玄学每一步都能用回测兜底愿你少走这几个坑希望帮到你。本文还有配套的精品资源点击获取