企业官网建设流程全解析

简介：基于拉丁超立方抽样与样本削减的场景分析方案，面向从事可再生能源出力预测、负荷预测及电力系统运行研究的工程师与科研人员。压缩包为rar格式，仅含1个MATLAB脚本（.m），文件体积约1KB，聚焦核心算法实现，便于直接阅读与复用。已有2593人学习该资源。脚本围绕风光出力及负荷场景的生成与分析展开：先通过拉丁超立方抽样对风速、辐照度、温度等气象条件进行多维均匀采样，构建代表性强的不确定性输入样本；再结合样本削减策略，从大规模候选场景中筛选典型子集，降低计算开销而不牺牲预测精度。同时覆盖数据预处理、异常值清洗、缺失值填充及标准化等环节，并配有清晰注释，适合需要快速上手场景缩减方法或将其嵌入预测模型的读者。该方案兼顾多因素复杂性与计算效率，对提升风光并网稳定性、辅助电网供需平衡评估具有实践参考价值。

1. 拉丁超立方抽样为什么成了风光场景生成的首选：从蒙特卡洛跑不动说起

做微电网储能容量配置时，最折磨人的不是设备选型，而是风光出力说不清：风不吹、光不照、负荷还天天变。拿历史一年数据直接跑，场景太少等于看运气；上蒙特卡洛抽一万个场景，每次优化都像重新蹲一个算例，算力成本直接失控。拉丁超立方抽样（LHS）用分层采样的方式把样本量压到十分之一，样本削减再把上千个候选场景浓缩成十几个代表性典型场景——这就是「基于拉丁超立方抽样和样本削减的风光出力及负荷场景预测分析」这条线的核心价值。这篇笔记按我实际做过的方案讲：分布怎么建、LHS 怎么调、削减怎么不丢极端信息、削减后怎么算出能说服人的置信区间，给做微电网和电网规划的工程师当一份抄作业清单。

2. 拉丁超立方抽样的落地实现：从数据拟合到相关性控制

2.1 分布怎么选：Weibull、Beta 和正态的拟合细节

场景生成的第一步不是抽样，而是把历史数据变成可采样的分布。风速我一般用双参数 Weibull，固定位置参数 loc=0，只拟合形状参数 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, floc=0) # 风速: 固定位置项 beta_a, beta_b, _, _ = stats.beta.fit(irradiance / 1000.0, floc=0, fscale=1) 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 保证每个变量的覆盖范围被强制铺满，不像纯随机抽样那样可能整段尾部没采到。参数上，floc=0是风速拟合里最关键的约束，不固定位置项可能拟合出负的起始风速，后续采样全是残次品；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, axis=0), axis=0) + 1 S = stats.norm.ppf(ranks / (n + 1)) # 当前 S 的秩相关矩阵及其 Cholesky 分解 corr_S = np.corrcoef(S, rowvar=False) 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, seed=42): km = KMeans(n_clusters=n_target, n_init=50, random_state=seed) labels = km.fit_predict(scenes) centers = km.cluster_centers_ count = np.bincount(labels, minlength=n_target) probs = count / count.sum() return centers, probs, labels

n_init=50是我从翻车里换来的习惯，默认值 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(axis=0) red_mean = (red_probs.reshape(-1, 1) * red_scenes).sum(axis=0) 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(axis=0) def weighted_quantile(scenes, probs, q): order = np.argsort(scenes, axis=0) 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 归一化分子分母搞混；另一种是采样后没做截断就进了削减。解决分两步：拟合时风速固定floc=0，辐照归一化用交流侧上限，采样后统一 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(f"KS p={p_value:.3f}, Wasserstein={wd:.3f}")

第二层是分布对比。KS 检验看削减场景与历史的累计分布是否一致，Wasserstein 距离量化两个分布之间的搬运成本。p 值小于 0.05 说明削减后场景集和原始数据分布差异显著，这时要回头调削减数或分布参数。第三层是决策结果收敛性验证：用 5、10、20 个场景分别跑同一个储能配置问题，看最优容量变化是否在可接受区间内。如果 10 个和 20 个场景的结果差得很大，说明 10 个场景承载不了这个问题的不确定性。

我自己翻车最多的就是第二层检验栽在 K-means 聚类中心上，聚类中心是平滑点，一对比分布就发现尾部全没了。后来养成习惯：削减完先画典型场景曲线图，肉眼扫一遍再往下走。场景削减不是玄学，每一步都能用回测兜底，愿你少走这几个坑，希望帮到你。

本文还有配套的精品资源，点击获取