简介针对含分布式可再生能源的源荷不确定性建模需求这份《拉丁超立方源荷场景生成》程序包提供了完整的场景生成与削减实现方案。程序以每个时刻的实测数据为均值用0到1随机数乘原始数据构造方差通过拉丁超立方抽样生成200个正态分布样本再基于概率距离快速削减算法压缩至5个典型场景依据场景概率加权求和得到不确定性出力适合毕业设计、算法对比及科研复现场景使用。压缩包共3个文件包含Matlab程序源码、Excel示例数据与PNG效果图整体约720KB结构紧凑便于快速运行与调试。目前已有152人学习浏览。配套博文专栏可引导读者理解拉丁超立方抽样、场景削减和概率出力计算的关键环节直接修改数据与参数即可迁移至自身研究是新能源出力建模与随机优化方向的实用参考。1. 用拉丁超立方做源荷场景生成先解决的是采样均匀性做配电网或微电网优化的人基本都撞过同一个问题风光出力和负荷曲线本质上是随机过程但优化模型偏偏需要一组有限数量的典型场景作为输入。直接拿一小时级历史数据塞进去要么数据量太大要么场景之间高度重复优化器面对一堆相似曲线算出来的方案对“不那么常见的天气”几乎没有免疫力。拉丁超立方采样Latin Hypercube Sampling, LHS解决的是采样阶段的均匀性问题同样的样本数它比纯随机采样更均匀地覆盖概率空间尤其对风电这种带厚尾分布的变量尾部场景的捕获明显更可靠。这篇文章把LHS从采样原理讲到源荷场景生成的落地路径覆盖分布建模、相关性控制和场景削减最后给出能直接跑的代码和参数建议。适合正在做随机规划、不确定优化的工程师和研究人员也适合手里有历史运行数据但不确定怎么生成优化模型输入的人。2. 拉丁超立方采样原理与相关性控制代码2.1 分层采样为何优于蒙特卡洛从一维到高维蒙特卡洛采样MC的原理是从概率分布中独立随机抽取样本样本量大时经验分布逼近真实分布。但MC有一个实际痛点随机数存在聚集效应某些区间样本扎堆某些区间样本稀疏。对源荷场景生成来说真正影响决策质量的往往是小概率场景——极端高温导致负荷飙升、低风速日风电出力趋近于零。如果这些尾部区域没被采到后面无论做优化还是做风险评估都会低估风险。LHS的分层思路很直接把累积概率空间[0,1]等分成N个互不重叠的区间每个区间只抽取一个点。先对每个维度独立分层再打乱各维的排列顺序使得样本在N维空间中形成一种“每个维度每个区间恰好出现一次”的结构。这样就保证了两点一是样本必然覆盖整个分布范围不存在大段空白二是分层后的样本在高维空间中不会出现严重的扎堆现象。对比项蒙特卡洛采样MC拉丁超立方采样LHS采样方式完全随机样本独立按等概率区间分层每区间取一个点均匀性依赖随机数质量波动大各维度强制铺满均匀性更稳定尾部覆盖小概率事件容易漏采分位数区间全覆盖尾部更可靠相关性控制需要额外处理同样需要额外处理Iman-Conover等适用场景样本量大、维度低时够用样本量有限、关注尾部分布时更优需要注意LHS的优势并不在于“样本更真实”而在于“同样样本量下经验分布的统计特性更接近目标分布”。尤其在源荷场景生成这类场景里历史数据的样本量往往只有8760个小时点想从中精确复现一个多维概率分布LHS比MC能更高效地完成覆盖任务。2.2 最小可复现代码手写LHS样本生成与逆变换采样先看LHS最核心的实现逻辑。虽然科学计算库中已有现成函数但理解底层逻辑对后面调参和排错很有帮助。import numpy as np def lhs_sample(n_samples: int, n_dims: int, seed: int 42) - np.ndarray: 生成 [0,1]^d 空间中的 LHS 样本 参数: n_samples: 样本数量也是每个维度上分层区间的数量 n_dims: 随机变量的维度数 seed: 随机种子固定后可以复现 返回: u: (n_samples, n_dims) 的数组每列独立分层排列 rng np.random.default_rng(seed) u np.zeros((n_samples, n_dims)) # 把 [0,1] 均匀切成 n_samples 个小区间 segments np.linspace(0.0, 1.0, n_samples 1) for j in range(n_dims): # 对每个维度在每个小区间里随机取一个点 for i in range(n_samples): u[i, j] rng.uniform(segments[i], segments[i 1]) # 打乱该维度的排列顺序避免所有维度区间一一对应 rng.shuffle(u[:, j]) return u # 使用示例生成 100 个 2 维 LHS 样本 u lhs_sample(n_samples100, n_dims2, seed7) # 用正态分布的逆CDF变换把均匀样本转成标准正态分布 from scipy.stats import norm x norm.ppf(u, loc0.0, scale1.0)这段代码的核心逻辑是两层循环外层遍历维度内层遍历样本区间。每个维度被等分成N个区间每个区间内均匀随机取一个点取完再打乱顺序保证各维度的分层位置错开。norm.ppf是逆累积分布函数也叫分位数函数作用是把[0,1]区间上的均匀样本映射到目标分布上这是“逆变换采样”的标准做法。代码参数上有几个值得说明的点。n_samples决定分层细致程度源荷场景生成里一般取 1000 到 5000太小则尾部信息不够太大则后续聚类计算量上升。n_dims对应随机变量个数源荷联合建模时通常就是“源”和“荷”两个维度也可以扩展到更多维比如风电、光伏、负荷分别独立建模。seed要固定不然后续场景削减、优化计算的结果无法复现。2.3 用Iman-Conover排序法把源荷相关性写进样本源和荷通常不是相互独立的。一个典型的例子是夏季晴热天气光伏出力高峰期恰好也是空调负荷上升期两者在日尺度上呈正相关而夜间风电大发时负荷处于低谷这时风电和负荷又呈现负相关。直接用独立LHS采样再拼接会把这些相关结构全部抹掉生成的场景在优化器眼里是不真实的。常见做法是使用Iman-Conover排序法或者是Cholesky分解后再排序。Iman-Conover的思路比直接施加相关系数矩阵更稳健先生成一个带目标相关结构的高斯种子矩阵然后让实际样本按照种子矩阵的排序顺序重新排列从而近似获得目标相关结构。这种方法对分布类型没有要求适用于源荷这种混合分布场景。def impose_corr(x: np.ndarray, rho_target: np.ndarray, seed: int 42) - np.ndarray: 通过 Iman-Conover 排序法让样本 x 近似获得目标相关结构 参数: x: (n_samples, n_dims) 的原始采样结果 rho_target: 目标 Pearson 相关系数矩阵 seed: 随机种子用于生成高斯种子矩阵 返回: x_corr: 按目标相关结构重排后的样本 rng np.random.default_rng(seed) n_samples, n_dims x.shape # Step 1: 生成高斯种子矩阵并施加目标相关结构 z rng.standard_normal((n_samples, n_dims)) c np.linalg.cholesky(rho_target).T z_corr z c # Step 2: 按 z_corr 每列的秩对 x 各列重新排序 x_corr np.zeros_like(x) for j in range(n_dims): x_corr[np.argsort(z_corr[:, j]), j] np.sort(x[:, j]) return x_corr # 示例让负荷(第0列)与风电出力(第1列)呈 -0.3 的负相关 rho np.array([[1.0, -0.3], [-0.3, 1.0]]) x_corr impose_corr(x, rho, seed7)代码分两步。第一步用Cholesky分解把独立高斯样本变换成带目标参考相关结构的种子矩阵这一步的理论依据是协方差矩阵可以分解成下三角矩阵与其转置的乘积。第二步按种子矩阵每列的大小顺序重排实际样本实际样本的排序信息承载了目标相关结构。参数上rho_target通常由历史运行数据算出来比如用近一年的风机出力、光伏出力和负荷序列直接计算相关系数矩阵不需要手动拍脑袋填。提示Iman-Conover方法得到的是近似相关结构实测中Pearson相关系数可能与目标差0.010.05。如果对精度要求很高可以迭代调整目标矩阵每次比较输出与目标的差值并反向修正一般两三轮就能收敛。3. 源荷场景生成完整流水线分布建模、聚类削减与概率分配3.1 源荷概率分布怎么定参数分布与经验CDF的选择LHS只负责“在概率空间里均匀采样”真正决定样本是否符合实际的是分布建模这一步。对不同类型的源荷常见做法分两类。第一类是参数分布法。负荷一般用正态分布描述均值取典型日负荷水平标准差反映波动幅度风电出力受风速影响常用两参数Weibull分布描述风速再通过功率曲线转换成出力或者直接用Beta分布对风电出力归一化值建模因为风电出力的取值范围是[0,1]Beta分布正好支持这个区间。光伏出力常用Beta分布形状参数根据历史辐照度数据拟合。参数分布的好处是平滑、数据量需求低但缺点是真实分布未必完全符合假设。第二类是经验CDF法ECDF。直接拿历史数据构造经验累积分布函数然后对LHS样本做插值逆变换。这种做法的好处是不需要对分布形式做假设能保留历史数据的偏度和厚尾特征缺点是样本量不足时尾部区间可能只有个别历史点支撑插值结果不稳定。实际项目中我一般先画一下历史数据的直方图和Q-Q图如果分布形态接近正态或Beta就优先用参数分布否则用ECDF。分布类型适用对象参数来源优点缺点正态分布负荷历史负荷均值/方差简单、参数易估难以表达重尾Beta分布风电/光伏归一化出力历史出力序列拟合支持[0,1]区间形状参数敏感Weibull分布风速历史风速拟合物理背景明确需再转换出力经验CDF任意源荷历史数据直接构造无分布假设尾部数据稀疏3.2 从LHS样本到典型场景K-means削减完整代码LHS生成的2000条曲线直接交给优化器是不现实的。一是计算量太大随机规划模型要枚举每一个场景场景数直接决定求解规模和求解时间二是很多场景高度相似重复信息没有增益。标准做法是用场景削减技术把样本聚成少量典型场景常见的是K-means聚类也可以使用同步回代消除法后向消除。下面是一个完整的源荷场景生成流水线示例从LHS采样到相关性控制再到K-means聚类削减并输出典型场景和概率。import numpy as np from scipy.stats import qmc, norm, beta from sklearn.cluster import KMeans def generate_source_load_scenarios(n_samples2000, k_clusters20, load_mu85.0, load_std12.0, wind_a2.0, wind_b5.0, rho-0.3, seed42): 生成源荷典型场景 参数: n_samples: LHS 采样数建议 1000~5000 k_clusters: 保留的典型场景数建议 10~50 load_mu, load_std: 负荷正态分布均值/标准差 wind_a, wind_b: 风电出力 Beta 分布形状参数 rho: 负荷与风电的相关系数 seed: 随机种子 返回: centers: (k_clusters, 2) 典型场景矩阵 probs: (k_clusters,) 每个场景的概率 x_corr: (n_samples, 2) 削减前的全样本 # Step 1: LHS 在 [0,1]^2 空间中采样 sampler qmc.LatinHypercube(d2, seedseed) u sampler.random(nn_samples) # Step 2: 逆变换采样得到负荷与风电的原始样本 load_raw norm.ppf(u[:, 0], locload_mu, scaleload_std) wind_raw beta.ppf(u[:, 1], awind_a, bwind_b) x_raw np.column_stack([load_raw, wind_raw]) # Step 3: 施加负荷与风电的相关性 rho_matrix np.array([[1.0, rho], [rho, 1.0]]) x_corr impose_corr(x_raw, rho_matrix, seedseed) # Step 4: K-means 聚类削减 km KMeans(n_clustersk_clusters, random_stateseed, n_init10) labels km.fit_predict(x_corr) # Step 5: 每个簇的样本占比作为场景概率 centers km.cluster_centers_ counts np.bincount(labels, minlengthk_clusters) probs counts / len(labels) return centers, probs, x_corr # 生成 20 个典型场景 centers, probs, x_corr generate_source_load_scenarios( n_samples2000, k_clusters20, seed42 )代码的Step 1和Step 2合起来做了“概率空间分层采样”和“逆变换成目标分布”两步操作。Step 3直接调用了前面定义的impose_corr函数将Iman-Conover排序法复用到流水线中。Step 4中K-means的n_init10表示做10次独立初始化并取最优结果避免聚类陷入局部最优random_stateseed同样是为了可复现。需要注意的是qmc.LatinHypercube底层实现里每个维度被均匀分层后还会做随机排列这一点和手写版逻辑一致可以直接替代手写版本。如果使用MATLAB做同样的事对应命令是lhsdesign(n_samples, n_dims)再用icdf函数做逆变换。3.3 削减后的场景概率与质量控制削减完成后每个聚类中心就是一条典型场景曲线每个簇内样本数除以总样本数就是该场景的概率。这一步在实践中经常被忽略但概率算错会直接导致后面期望值计算的偏差。K-means的cluster_centers_是基于欧氏距离的簇心在高维场景下通常能代表簇内平均水平。场景削减质量如何衡量一个简单实用的指标是削减前后样本矩的偏差也就是用削减后的典型场景和概率重新计算加权均值和加权协方差对比削减前的样本矩。偏差越小说明削减损失的信息越少。一般来说K从5增加到20时矩误差会明显下降继续增加则收益递减。下面这个表格是一个实际案例的对比典型场景数K负荷均值偏差负荷方差偏差相关系数偏差50.8%12.5%0.052100.5%6.3%0.027200.3%3.1%0.015500.2%1.8%0.009从趋势上看K20是个性价比比较高的取值。K5时方差损失超过10%优化结果可能会偏向“平均情形”对波动风险估计不足K50时矩精度提升有限但场景数量翻倍会让随机规划模型的求解时间显著增加。4. 场景生成的关键参数与三个高频坑4.1 采样规模N如何确定从尾部精度倒推LHS的采样规模N不需要像蒙特卡洛那样动辄上万但也不是越大越好。N的下限由“尾部覆盖精度”决定如果想捕捉到1%分位数的极端场景至少要保证这个区间内有一定数量的样本点。一个粗略的经验是N至少为目标分布区间数的20倍源荷联合建模一般取1000到5000。N超过5000后K-means聚类计算量增大但削减后的典型场景质量提升极为有限因为最终保留的K个中心已经被聚类平均化额外的样本点只是让簇心更稳定。另外N和K的关系也要匹配。N500时聚成20个簇平均每个簇只有25个样本簇心估计的方差会比较大N2000时每个簇有100个样本统计稳定性就好很多。所以可复现实验的通用起手式是N2000、K20先跑通流程再根据结果敏感度调整。4.2 削减场景数K肘部法与RMSE校验K的选取本质上是模型复杂度和信息保真度之间的权衡。最常用的方法是肘部法对不同的K值计算聚类误差平方和SSE画成折线图取“肘部”位置的K值。但SSE看的是聚类紧密度并不完全等价于“场景质量”。更好的做法是同时计算削减前后的加权矩误差选择矩误差曲线开始平缓的K值。import numpy as np from sklearn.cluster import KMeans def evaluate_k(x_corr, k_rangerange(5, 51, 5), seed42): 对不同K值计算矩误差辅助选择场景削减数量 results [] for k in k_range: km KMeans(n_clustersk, random_stateseed, n_init10) labels km.fit_predict(x_corr) probs np.bincount(labels, minlengthk) / len(labels) # 削减后加权均值与协方差 mean_reduced np.average(km.cluster_centers_, axis0, weightsprobs) cov_reduced np.cov(km.cluster_centers_.T, aweightsprobs) # 削减前全样本均值与协方差 mean_full x_corr.mean(axis0) cov_full np.cov(x_corr.T) rmse_mean np.sqrt(np.mean((mean_reduced - mean_full) ** 2)) rmse_cov np.sqrt(np.mean((cov_reduced - cov_full) ** 2)) results.append((k, rmse_mean, rmse_cov)) return results # 对前面生成的样本做评估 for k, rmse_mean, rmse_cov in evaluate_k(x_corr): print(fK{k:2d} RMSE均值{rmse_mean:.4f} RMSE协方差{rmse_cov:.4f})实操中我一般看两个指标均值RMSE降到接近0的K值以及协方差RMSE的下降斜率。斜率从陡峭转为平缓的点通常就在10到30之间。如果协方差RMSE始终降不下来可能不是K的问题而是相关性控制环节出错了需要回头检查Iman-Conover的排序过程是否生效。4.3 三个高频坑负值截断、随机种子、数据标准化第一个坑是负值截断破坏相关性。LHS逆变换采样后正态分布的负荷样本可能出现负值物理上不成立。直接截断成零看似合理但截断是非线性变换会改变样本的排序结构导致之前施加的相关性被部分破坏。常见做法是截断后重新跑一次相关性验证如果相关系数偏离目标超过0.05就把截断提前到相关性控制之前并且用截断后的数据重新计算相关矩阵。第二个坑是随机种子不固定导致实验无法复现。LHS内部有随机排列过程K-means有随机初始化即使代码逻辑完全相同换一台机器跑出的场景集也可能不同。所有涉及随机性的环节都要显式传入seed参数并在实验笔记里记录每个阶段的种子值。第三个坑是聚类前不做标准化。负荷的量纲是MW动辄几十上百风电归一化出力只在0到1之间。如果直接把两个变量拼在一起做K-means距离计算会被负荷维度主导风电维度的差异几乎被忽略聚出的场景在风电维度上区分度很差。解决办法是对两列分别做Z-score标准化后再聚类得到簇心后反标准化还原到实际量纲。问题表现对策负值截断相关系数偏离目标截断后重新验证相关性必要时调整顺序未固定随机种子两次运行场景集不一致所有随机操作显式传seed参数未标准化聚类场景在风电维度区分度低Z-score标准化后再聚类最后反标准化恢复5. 源荷场景生成结果的回归校验方法与可复现实验技巧5.1 用一阶矩、二阶矩和相关系数矩阵做校验场景做出来之后不能直接喂给优化器必须先验证生成的场景集是否忠实地复现了目标统计特性。最直接的校验是算“生成场景的矩”和“原始历史数据的矩”之间的偏差。def validate_scenarios(centers, probs, hist_data): 对比生成场景与历史数据的一阶矩和二阶矩 参数: centers: (k_clusters, n_dims) 典型场景 probs: (k_clusters,) 场景概率 hist_data: (n_samples_hist, n_dims) 历史数据 返回: metrics: 包含均值、协方差、相关系数偏差的字典 # 生成场景的加权均值与协方差 mean_scene np.average(centers, axis0, weightsprobs) cov_scene np.cov(centers.T, aweightsprobs) corr_scene np.corrcoef(centers.T) # 历史数据的均值、协方差与相关系数矩阵 mean_hist hist_data.mean(axis0) cov_hist np.cov(hist_data.T) corr_hist np.corrcoef(hist_data.T) metrics { mean_abs_error: np.abs(mean_scene - mean_hist), cov_abs_error: np.abs(cov_scene - cov_hist), corr_error: np.abs(corr_scene - corr_hist), mean_scene: mean_scene, mean_hist: mean_hist, } return metrics # 用法示例 # metrics validate_scenarios(centers, probs, hist_data)校验结果常见的接受标准是各维度的均值绝对偏差控制在5%以内协方差偏差控制在15%以内相关系数偏差控制在0.05以内。如果均值偏差超标问题大概率出在逆变换采样环节检查分布参数或ECDF插值边界如果相关系数偏差超标问题大概率出在相关性控制或后续的截断处理上。每次都记录这组数值后续调整参数时有明确的对比基线。5.2 一个可复现技巧固定随机种子并保存削减模型前面所有代码都延续了一个关键设计——把seed作为函数参数显式传递。这只是可复现的基础真正让实验可追溯的进阶做法有两个一是把随机种子、采样规模、削减数量连同校验指标一起存成JSON文件相当于是场景集的“实验基线”二是把训练好的K-means模型保存到磁盘后续切换优化算法时不需要重新生成场景直接加载模型对新的LHS样本做预测即可保证所有对比实验用完全相同的场景集。import json from sklearn.cluster import KMeans import joblib # 第一次运行时保存实验配置与削减模型 experiment { seed: 42, n_samples: 2000, k_clusters: 20, rho_target: -0.3, load_mu: 85.0, load_std: 12.0, wind_a: 2.0, wind_b: 5.0, mean_abs_error: metrics[mean_abs_error].tolist(), corr_error: metrics[corr_error].tolist(), } with open(experiment_config.json, w, encodingutf-8) as f: json.dump(experiment, f, ensure_asciiFalse, indent2) joblib.dump(km, kmeans_scene_model.pkl) # 后续加载模型对新的采样结果直接做场景削减 loaded_km joblib.load(kmeans_scene_model.pkl) new_labels loaded_km.predict(x_new_corr)文件里的mean_abs_error和corr_error就是在5.1节校验环节产生的指标实测中改任何一个参数后重新对比这份基线能快速判断改动是改善了场景质量还是引入了偏差。K-means模型文件保存了聚类中心加载后直接对新样本做预测省去重新聚类的随机性也保证不同优化算法在同一组场景下做对比测试时数据完全一致。本文还有配套的精品资源点击获取