简介基于拉丁超立方抽样与样本削减的风光出力及负荷场景预测脚本面向电力系统研究人员、风光电站工程师及相关专业学生用于在风速、辐照度等多气象条件变化下模拟新能源出力与负荷的不确定性从而评估不同场景对电网稳定性的影响。压缩包内仅含 1 个 .m 脚本文件体积 1KB但代码流程覆盖数据清洗、拉丁超立方采样、样本削减、模型训练及误差评估等关键步骤结构清晰注释详细。该资源已有 2591 人学习浏览说明其方法在相关领域具有一定关注度与参考价值尤其适合作为新能源出力预测相关课题的入门示例。通过研读与运行该脚本可以直观理解如何构建代表性场景集并在保持预测精度的同时降低计算开销也便于后续在此基础上扩展算法或适配自身数据。1. 场景预测分析风光出力和负荷的不确定性不能靠拍脑袋电网调度和新能源规划中风电出力受风速影响、光伏出力受辐照度影响负荷本身也在波动。三重不确定性叠加后用单个预测值做决策往往会偏保守或偏冒险。蒙特卡洛模拟能刻画不确定性但几千上万个场景让优化计算直接卡死。拉丁超立方抽样用分层机制以少量样本覆盖全部取值区间再通过样本削减把几百个场景压缩到几十个代表场景是工程上兼顾精度与计算成本的常用路线。下面按这条链路把原理、参数和落地代码讲透。适合做新能源消纳分析、储能容量配置和电网规划的技术人员。2. 拉丁超立方抽样的分层机制与风光出力概率建模蒙特卡洛的问题在于“一视同仁”在随机数空间里样本可能扎堆在某段区间另一端没有覆盖到。拉丁超立方抽样把每一维变量的取值区间等分成N份保证每个区间恰好被取到一次再通过排列组合生成多维样本。这使得相同样本量下LHS对边缘分布的覆盖程度远好于简单随机抽样。2.1 从蒙特卡洛到LHS样本覆盖率的差距蒙特卡洛收敛速度与样本量的平方根成反比想要提高一位精度通常要增加百倍样本量。在场景分析中每个场景最终都要代入优化模型样本量直接决定求解时间所以纯蒙特卡洛方案在工程上基本不落地。LHS的核心操作如下把[0,1]分成N个等宽区间在每个区间内取一个均匀随机数对每一维独立做同样操作后把各维数值随机打乱并拼成N行d列的样本矩阵。这样每个变量的边缘分布都得到完整覆盖。提示LHS适用于所有能用逆变换采样的分布。如果拿到的是历史数据而非分布参数可以先做经验分布或者核密度估计再在经验分布上做逆变换采样。2.2 风速、辐照度与负荷的分布假设做风光出力场景前先要对源数据建模。常见做法是风速用两参数威布尔分布形状参数k和尺度参数c用历史风速数据做极大似然估计风电功率曲线再按切入风速、额定风速、切出风速做分段映射。光伏出力的辐照度常用Beta分布拟合负荷则更多直接用历史曲线的均值加减扰动或者用正态分布近似。这里有个工程细节直接对功率做分布拟合误差往往很大因为风电功率在切入风速以上才非零很多时刻出力为0分布函数在0点有概率质量堆积。我一般会先对功率序列做样本筛选把零功率和满发状态单独处理只对中间区间的连续出力做分布拟合再与状态概率组合。否则后面的LHS采样会把这些阶跃特征抹平。2.3 用秩相关系数与Cholesky分解生成相关场景风速、辐照度和负荷并不是独立的。比如夜间风速往往比日间大光伏出力与负荷峰值在白天有正相关。直接用独立LHS采样会切断这些耦合关系导致生成的场景在物理上不成立。标准做法是先做独立标准正态抽样再乘上相关系数矩阵的Cholesky分解得到带相关性的正态样本再转换为均匀数并反变换到目标分布。import numpy as np from scipy.linalg import cholesky from scipy.stats import norm n_samples, n_dims 500, 3 # 独立LHS样本0-1均匀空间 def lhs_basic(n, d, rngNone): rng rng or np.random.default_rng(42) out np.empty((n, d)) for j in range(d): u rng.uniform(0, 1, n) ranks rng.permutation(n) # 分层用排序编号让每列覆盖所有区间 out[:, j] (ranks u) / n return out u_lhs lhs_basic(n_samples, n_dims) # 转标准正态乘Cholesky分解引入相关性 z norm.ppf(u_lhs) corr np.array([[1.00, 0.35, 0.20], [0.35, 1.00, 0.45], [0.20, 0.45, 1.00]]) L cholesky(corr, lowerTrue) z_corr z L.T u_corr norm.cdf(z_corr)这段代码的要点lhs_basic里ranks是0到n-1的随机排列加上u让采样点落在小区间内而不是区间端点这样的分层效果更稳定相关系数矩阵corr需保证半正定否则Cholesky分解会报错u_corr已不再是严格LHS分层但保留了大部分分层特性同时具备了目标秩相关。接下来按变量做逆变换得到风、光、负荷场景。风速按威布尔反函数、辐照度按Beta反函数、负荷按正态反函数分别处理from scipy.stats import weibull_min, beta, norm # 假设历史数据估计出的参数 wind_k, wind_c 2.2, 8.5 # 威布尔形状、尺度 beta_a, beta_b 0.85, 1.30 # 辐照度Beta分布形状参数 load_mu, load_sigma 5200, 650 # 负荷的均值与标准差 scen np.column_stack([ weibull_min.ppf(u_corr[:, 0], wind_k, scalewind_c), beta.ppf(u_corr[:, 1], beta_a, beta_b) * 1000, # 折算为W/m² norm.ppf(u_corr[:, 2], load_mu, load_sigma) ])ppf是逆累积分布函数把均匀数映射到目标分布。场景矩阵三列分别是风速、辐照度、负荷后续需要结合风功率曲线把风速换算为出力并按光伏组件面积和效率把辐照度换算为功率。参数不同会导致场景形态差异明显所以这里的分布参数应该用最大似然估计而不是手填。3. 样本削减的数学逻辑与两种落地算法场景生成完以后500个场景对优化模型来说仍然太多。样本削减的目标是找到数量更少的场景使新旧场景集合的概率分布距离尽量小而不是简单挑几个“最像”的场景。工程上用得最多的是同步回代消除SBR和K-means聚类两类方法前者在场景概率质量的操作上更细腻后者实现简单、可解释性好。3.1 场景削减要保留什么概率距离而不是形态相似场景削减的评价准则是削减前后的概率分布距离。常用的是Kantorovich距离直观解释是新场景集Q替代原场景集P需要付出的“运输代价”。当场景服从离散概率分布时它退化为线性规划中的最优传输问题。SBR算法不求解完整线性规划而是用贪心策略逼近最优解实践效果好且计算开销小。如果只看形态相似性比如把所有场景聚类后取簇心很可能丢掉稀有的高影响场景例如极端大风场景。这类场景概率低但后果严重在调度鲁棒性分析中恰恰是要保留的对象。所以削减算法的选择必须和下游用途绑定做风险分析时建议用SBR因为它的概率转移逻辑会保留边缘场景做快速日前计划时用K-means更划算。3.2 SBR同步回代消除的实现与参数SBR的核心逻辑是反复删除“代价最小”的场景。每轮从当前场景集中找一对场景i和j其中j是i的最近邻删除i后i的概率质量并入j代价是即p_i乘以i到j的距离。重复直到场景数降到目标值。def sbr_reduction(scenarios, probs, n_keep): scenes scenarios.copy() p probs.copy() while len(scenes) n_keep: n len(scenes) # 计算场景间欧氏距离矩阵 dist np.sqrt(((scenes[:, None, :] - scenes[None, :, :]) ** 2).sum(-1)) np.fill_diagonal(dist, np.inf) # 排除自身 min_dist dist.min(axis1) # 每个场景到最近邻的距离 cost p * min_dist # 删除代价 i np.argmin(cost) # 找到代价最小的场景 j dist[i].argmin() # i的最近邻 p[j] p[i] # 概率转移给最近邻 scenes np.delete(scenes, i, axis0) p np.delete(p, i, axis0) return scenes, p这段代码的复杂度是O(N²·d)的距离矩阵加O(N³)的迭代500个场景削减到50个大约需要几秒可以接受。n_keep是保留场景数取决于下游随机优化的求解能力p数组用于记录各场景权重削减后p的和恒为1。注意每次删除后要重新计算距离矩阵不能复用旧矩阵因为删除操作会改变场景间的邻近关系。提示SBR默认使用欧氏距离。如果原始场景包含量纲差异巨大的变量比如风速0到20m/s和负荷3000到8000kW先做标准化再计算距离否则负荷维度会主导整个削减过程。3.3 K-means快速削减与SBR的取舍K-means做场景削减的思路更简单将全部场景聚类聚类中心作为代表场景每个簇的样本占比作为该场景的概率。from sklearn.cluster import KMeans def kmeans_reduction(scenarios, n_keep, random_state42): km KMeans(n_clustersn_keep, random_staterandom_state, n_init10) labels km.fit_predict(scenarios) probs np.bincount(labels, minlengthn_keep) / len(labels) return km.cluster_centers_, probs两者对比下来K-means快但聚类中心是簇内均值很可能落在数据的“空洞”里产生物理上不存在的场景组合SBR保留的是原始采样点每个对应一组实际可发生的风光负荷数值。从下游模型的角度看SBR的输出更可信所以我在做可靠性分析时默认选SBRK-means适合在场景数量特别大、需要快速粗筛的场景先聚类到几百个再交给SBR精修。3.4 削减数量怎么定n_keep的经验区间削减数量没有标准答案一般按下游问题的计算瓶颈反推每个场景对应一次确定性优化求解保留数乘上单场景求解时间就是总耗时以此定上限再根据精度要求定下限。保留场景数适用场景5~10日前调度、规划阶段的粗粒度场景10~30随机经济调度、储能容量优化30~80概率风险评估、机组组合80~150需要输出概率分布细节的可靠性分析用SBR把500个场景削减到20个时削减后场景集与原集合的概率距离应小于原集合总方差的5%到10%。这个值可以在削减时打印出来观察如果削减前后期望出力的偏差超过5%说明削得太狠需要提高n_keep。4. 风光出力及负荷场景预测分析的完整实现把前面几章串起来就是一条可落地的流水线历史数据 → 分布拟合 → LHS采样 → 场景削减 → 场景预测分析。本节给出一套最小可运行实现并解释每个环节的输入输出。4.1 从历史数据到场景生成第一步是从风电场SCADA数据、光照监测站数据、电网负荷曲线中分别估计分布参数。以风速为例用scipy的weibull_min.fit拟合历史风速序列得到形状参数k和尺度参数c。光伏分布参数用历史辐照度做Beta拟合负荷数据如果有明显的日周期特性可以按小时分段拟合而不是全部揉成一个正态分布。4.2 场景削减与概率分配沿用前面生成的场景矩阵先标准化再传入SBR或K-means。削减输出的代表场景要归一化概率随后用于预测分析。def build_scene_pipeline(flow_data, irrad_data, load_data, n_samp500, n_keep30): # 依赖第2节的lhs_basic和第3节的sbr_reduction from scipy.stats import weibull_min, beta, norm # 1) 分布拟合 wind_k, _, wind_c weibull_min.fit(flow_data, floc0) beta_a, beta_b, _, _ beta.fit(irrad_data / 1000, floc0, fscale1) load_mu, load_sigma norm.fit(load_data) # 2) LHS采样与相关性处理 u_lhs lhs_basic(n_samp, 3) z norm.ppf(u_lhs) z_corr z L.T # L来自第2节Cholesky分解 u_corr norm.cdf(z_corr) scen np.column_stack([ weibull_min.ppf(u_corr[:, 0], wind_k, scalewind_c), beta.ppf(u_corr[:, 1], beta_a, beta_b) * 1000, norm.ppf(u_corr[:, 2], load_mu, load_sigma) ]) # 3) SBR削减初始概率均分 scenes, probs sbr_reduction(scen, np.ones(n_samp) / n_samp, n_keep) return scenes, probsweibull_min.fit返回元组(shape, loc, scale)固定floc0避免拟合出不合理的位移项。beta.fit的floc0, fscale1表示把辐照度归一化到[0,1]后再拟合因为Beta分布定义域是[0,1]直接用原始数值拟合容易遇到数值问题。norm.fit返回均值和标准差负荷按高斯拟合是最省事的假设条件允许时建议用历史负荷的分时经验分布替代。4.3 用削减后的场景做预测分析实际场景预测分析按小时或调度时段展开。假设scenes是24个时段拼接后的场景矩阵形状为(n_keep, 72)每三列一组分别表示风电出力、光伏出力、负荷。常见有三种输出期望出力、预测区间、最坏场景。# 期望值按场景概率加权 expected scenes.T probs # 90%预测区间 p10 np.percentile(scenes, 10, axis0) p90 np.percentile(scenes, 90, axis0) # 最坏场景按净负荷最大挑选 net_load (scenes[:, 2::3] - scenes[:, 0::3] * 0.4 - scenes[:, 1::3]) worst_idx net_load.sum(axis1).argmax() worst_scene scenes[worst_idx]scenes.T probs是矩阵乘法把每个场景的取值按概率加权求和结果是一个72维向量。percentile给出经验分位数与概率加权分位数在小样本下略有差异但用于工程判断足够。净负荷计算里的0.4是风功率容量系数的占位值实际应从风电场功率曲线解析式得到这里仅演示最坏场景选取逻辑。4.4 预测结果的判断标准拿到期望和区间后要回归原始数据做校验。常见做法是留出最近一周的真实数据把预测区间与真实曲线叠加画图看覆盖率是否接近90%同时检查区间宽度是否过大。区间过宽说明削减保留的场景过于分散区间过窄则可能丢失了不确定性信息。如果只有一天的验证数据可以计算平均绝对误差和连续等级概率分数。后者能同时评价预测区间的锐度与校准度是比较合适的单指标。5. 样本削减质量验证与预测调参技巧最后一章讲三个在线工作中会用到的技巧。这些经验来自多次用场景削减做分析之后的复盘篇幅不长但都是实际问题。5.1 用概率距离与期望偏差双指标验证削减前记录原始场景集P削减后得到Q。SBR迭代中已经逐步记录了削减代价可以把这个代价作为内部指标输出。工程上更直观的是期望偏差比较削减前后风光出力均值若超过5%就需要增加n_keep或检查分布拟合是否畸形。建议把这两个指标作为流水线固定输出每次跑完先看一眼数值再决定要不要调参。5.2 样本数n_samp的交叉验证500是个常见起点但实际所需样本数与维度相关。一个实用做法是把n_samp设为200、500、1000分别在削减后计算期望和区间宽度看结果变化是否超过1%。如果200和1000的结果差异很小说明200就够了如果差异明显说明LHS还没覆盖到所有关键场景组合。注意每次比较要固定随机种子否则变异性会干扰判断。5.3 边界场景的处理LHS是均衡采样对极端事件如台风风速、连续阴雨天的刻画能力有限因为它均匀覆盖而非重点采样。处理方法是在LHS样本中加入一定比例的极端场景比如把风速分布的99%分位数单独抽取若干场景混入场景集再执行削减。这样最终保留的代表场景里既有常规工况也有小概率高风险场景用于规划类分析特别有用。如果削减后发现极端场景一个都没留下往往不是没采到而是SBR把它当成“边缘场景”消掉了解决办法是削减后额外强制保留最坏场景或调整距离权重。本文还有配套的精品资源点击获取