简介这份程序资源用于电力系统或综合能源系统中源荷场景生成适合做毕业设计、学术仿真及算法研究的读者。程序基于拉丁超立方抽样方法每个时刻抽取200个服从正态分布的样本均值取原始数据方差由0到1随机值乘以原始数据得到再结合概率距离快速削减算法将场景削减至5个并通过各场景概率与对应场景相乘求和量化不确定性出力。压缩包共3个文件包含可直接运行的m源程序、用于测试或比对的xlsx数据文件以及结果示意png图片包体积仅720KB轻量易部署。资源已有152人学习适合正处于建模攻坚或代码实现阶段的电力专业学生与科研人员。通过该程序可快速掌握拉丁超立方抽样与场景削减的完整流程减少从算法公式到代码落地的摸索时间。1. 拉丁超立方抽样在源荷场景生成里到底解决了什么问题做电力系统不确定性分析的人大概率都有过这样的经历用蒙特卡洛随机抽样生成风光出力和负荷场景抽两三千次曲线还是毛刺感十足场景削减之后概率分布又跟原始数据对不上。这个资源里的核心思路是把蒙特卡洛换成拉丁超立方抽样Latin Hypercube Sampling, LHS每个时刻只抽 200 个样本配合概率距离快速削减法砍到 5 个典型场景最后按「场景概率 × 对应场景出力」加权求和得到一条既能反映不确定性、又足够平滑的源荷出力曲线。适合正在做毕业设计、需要生成风光荷联合场景的电气工程学生也适合做鲁棒优化、随机规划但不想在场景生成上耗费太多算力的研究者。这个方法的关键在于LHS 不是随机撒点而是把每个变量的分布区间等概率分层然后在每一层里强制抽取样本。所以 200 个样本的覆盖效率往往比蒙特卡洛抽 2000 个还整齐。后续的削减也不是粗暴聚类而是用概率距离来衡量场景之间的相似度逐步合并距离最近的场景直到剩下 5 个。这套流程兼顾了分布拟合精度和计算开销几分钟就能跑完。2. LHS 分层抽样的数学原理与 yuanhe.m 实现2.1 分层采样的核心思想拉丁超立方抽样的本质是「等概率分层 层内随机」。假设某个时刻的负荷随机变量 X 服从均值为 μ、标准差为 σ 的正态分布要把它的取值范围分成 N 个互不重叠的区间每个区间的概率都是 1/N。然后对每个区间独立抽取一个样本点最后得到 N 个样本。这样一来无论 N 多大样本都保证覆盖了整个分布区间不会出现蒙特卡洛那种大量样本挤在均值附近、尾部稀疏的情况。这个资源里设定的是每个时刻抽取 200 个样本均值取原始数据方差则用一个 0 到 1 之间的随机数乘以原始数据。这种方差设置方式很有讲究直接把方差设成固定比例会让所有时刻的波动幅度都一个样生成出来的场景在时序上显得很假而用随机比例因子相当于对每个时刻的波动程度做了一次随机扰动场景库的整体多样性会明显提高。代价是方差的物理意义变得模糊但用在源荷不确定性建模里这种「半经验」做法其实很常见。代码里涉及的 MATLAB 实现大致遵循以下逻辑对每个时刻 t先调用 LHS 函数生成 200 个服从 N(0,1) 的标准化样本再按 Y μ σ × Z 的线性变换映射到实际出力区间。这里 σ rand × μrand 是 0 到 1 均匀随机数所以每个时刻的变异系数是一个随机值。2.2 核心抽样代码的逻辑拆解假设原始数据 shuju.xlsx 里第一列是时刻编号第二列到第四列分别是光伏、风电、负荷的原始出力那么对应的采样核心片段可以写成下面这样% 读取原始数据 data xlsread(shuju.xlsx); T size(data, 1); % 总时刻数 N 200; % 每个时刻的抽样次数 K 5; % 削减后的场景数 % 为每个时刻预分配样本矩阵 samples_pv zeros(N, T); samples_wind zeros(N, T); samples_load zeros(N, T); for t 1:T mu_pv data(t, 2); mu_wind data(t, 3); mu_load data(t, 4); % 方差取0到1随机数乘以原始数据 sigma_pv rand * mu_pv; sigma_wind rand * mu_wind; sigma_load rand * mu_load; % 生成均匀分层样本并映射到标准正态 u_pv ((1:N) - rand(N, 1)) / N; % LHS核心分层 层内随机偏移 u_wind ((1:N) - rand(N, 1)) / N; u_load ((1:N) - rand(N, 1)) / N; z_pv norminv(u_pv, 0, 1); z_wind norminv(u_wind, 0, 1); z_load norminv(u_load, 0, 1); % 线性变换到实际出力 samples_pv(:, t) mu_pv sigma_pv .* z_pv; samples_wind(:, t) mu_wind sigma_wind .* z_wind; samples_load(:, t) mu_load sigma_load .* z_load; end这里的 LHS 核心在于u ((1:N) - rand(N, 1)) / N这一行。每个秩 i 对应区间[(i-1)/N, i/N]rand是 0 到 1 的随机数所以抽样点在区间内部随机偏移而不是固定在区间中心。这样既保证了覆盖均匀性又保留了随机性。如果用((1:N) - 0.5) / N取区间中点就是确定性的分层采样会损失随机性场景多样性会差很多。norminv 函数的作用是把均匀分布的分位数转换成标准正态分布的分位数这样样本就服从 N(0,1) 了。之后再通过mu sigma * z完成仿射变换。注意这里是逐时刻独立抽样所以 t 循环里每轮生成的 200 个场景之间没有时间相关性属于完全独立的横向场景集。如果你需要场景在时序上有连续性就要引入 Cholesky 分解或 Copula 来处理相关性那是进阶做法。2.3 为什么方差要设计成随机比例因子固定方差的场景生成通常在削减后会出现一个现象概率最大的那个场景几乎就是原始数据的平滑版而小概率场景又偏离太远看起来像异常值。这是因为固定 σ 下大部分样本集中在均值附近削减算法很容易把大量相近样本合并成一个高概率场景导致多样性不足。随机比例因子相当于给每个时刻的分布宽度做了一次扰动有的时刻 σ 是 0.8 倍均值有的时刻是 0.2 倍均值。这样一来即便两个场景在某一时刻的均值接近它们在另一个时刻的波动幅度也可能差得很远削减算法能保留更多有区分度的场景。实际参数调整上如果你发现削减后的 5 个场景曲线太挤可以把 rand 改成 0.3 0.7 * rand强制最小方差比例不低于 0.3如果发现场景太散、出现负出力可以加一行 max(0, x) 截断。3. 基于概率距离的场景削减算法与参数选址3.1 为什么不用 K-means 而是概率距离削减很多课程设计里做场景削减直接上 K-means 聚类把 200 个场景聚成 5 类取每类的均值作为典型场景。这种做法的问题在于K-means 是基于欧氏距离的几何聚类它不感知场景的发生概率。如果某个区域样本点稀疏但物理意义重要K-means 可能会直接忽略它。而概率距离削减法是从场景集合的概率分布出发用 Kantorovich 距离简称 KD衡量两个场景之间的「搬运成本」每次合并距离最近的一对场景并把被合并场景的概率累加到保留场景上。这个资源里描述的是「基于概率距离快速削减算法」它的迭代逻辑可以用下面的伪代码描述% 输入sample_set为N×T矩阵每行是一个场景prob为N×1初始概率等权 % 输出削减后的场景及对应概率 while size(sample_set, 1) K % 计算所有场景两两之间的概率距离 dist_matrix zeros(M, M); % M为当前场景数 for i 1:M for j i1:M % 概率距离 欧氏距离 × 两个场景概率之积 dist_matrix(i,j) prob(i) * prob(j) * norm(sample_set(i,:) - sample_set(j,:)); dist_matrix(j,i) dist_matrix(i,j); end end % 找到距离最小的场景对 (i, j) [min_val, idx] min(dist_matrix(:)); [i_rm, j_keep] ind2sub(size(dist_matrix), idx); % 删除场景i_rm把概率累加到j_keep上 prob(j_keep) prob(j_keep) prob(i_rm); sample_set(i_rm, :) []; prob(i_rm) []; end这里的核心是prob(i) * prob(j) * 欧氏距离它同时考虑了场景之间的几何差异和概率权重。两个距离相近的场景如果它们各自的概率都很高合并的代价就大算法会优先保留它们而不轻易合并两个概率都极小的「边缘场景」哪怕距离稍远也可能被合并掉。这比 K-means 只按几何位置聚类要合理得多因为最终得到的概率分布能更好地保留原始场景集的统计特征。3.2 削减过程中的概率再分配策略上面的伪代码里删除场景 i 后直接把概率累加给 j这种策略叫「最近邻概率累加法」是最常用的一种。它的直观意义是被删掉的场景用离它最近的保留场景来代替因此它的概率应该转移给这个最近邻。实际工程里还有另一种做法是「按距离加权分配」即把被删场景的概率按距离比例拆给多个相邻场景这样做出来的概率分布更平滑但会破坏场景的稀疏性削减后概率值普遍偏小不利于后续计算期望值。我个人的建议是如果后续要做的是两阶段随机规划或机会约束规划用最近邻累加法就好如果只是生成一个代表性的出力曲线用来做确定性替代可以试试加权分配看哪条曲线更接近原始数据的均值曲线。不过资源里的场景概率和场景出力相乘求和本质上是求期望值所以最近邻累加法完全够用不需要额外复杂化。3.3 削减数量 K 的经验选择K 取 5 是个比较保守的工程选择。从数学角度看任何离散分布都至少需要 K 个支撑点才能保留 K 个矩特征。5 个场景可以保留前 4 阶矩的大致趋势包括均值、方差、偏度和峰度的粗粒度信息。对于光伏出力的场景生成来说5 个场景已经能刻画「晴天上午爬坡」「午后云层遮挡」「雨天低出力」这几类典型形态。如果 K 取 3削减后的场景曲线会在个别时刻出现突变因为可用的支撑点太少每个场景必须覆盖更大的概率质量曲线形状会变得比较生硬。K 取 8 或 10场景的多样性更好但后续随机优化模型里的 0-1 变量和连续变量数量会成倍增加求解速度下降。所以 K5 是精度和求解效率的折中。如果你跑出来的削减结果里某个场景的概率超过 0.5说明原始场景集合里存在一个绝对主导形态这时候可以适当调大 K 来看更细的分布。3.4 削减结果如何转成不确定性出力曲线资源和摘要里最关键的输出是「每个场景的概率与每个对应场景相乘求和得到不确定性出力」。这一步在代码里的实现非常直接% 削减后的场景矩阵 reduced_scenes: K×T概率向量 scene_prob: K×1 uncertainty_output scene_prob * reduced_scenes; % 1×T 的期望曲线 % 对比原始均值曲线 original_mean mean(samples, 1); % N×T样本集按列取均值 % 计算误差指标 mae mean(abs(uncertainty_output - original_mean), 2); % 平均绝对误差这行矩阵乘法的含义是第 k 个场景的出力曲线乘以它出现的概率然后对 K 个场景加权求和得到一个时刻 t 上的期望出力。这种不确定性出力的表达方式在随机规划里非常常用——它不直接给出一个确定的预测值而是给出一个考虑到各种可能场景后的加权期望配合场景概率一起输入到优化模型里能够天然处理不确定性对决策的影响。值得留意的是scene_prob * reduced_scenes得到的是期望值曲线它跟原始数据均值曲线之间的差距衡量了场景削减的信息损失。如果 MAE 太大说明削减后的 5 个场景不足以代表原始 200 个场景的分布特征这时候应该调大 K 或者检查方差设置。4. 从原始数据到完整场景的全流程搭建与参数矩阵4.1 数据文件 shuju.xlsx 的结构约定资源里附带的 shuju.xlsx 是输入数据源。从场景生成的角度它至少要包含三列光伏出力序列、风电出力序列、负荷序列。每一行对应一个时间段通常 24 小时就是 24 行96 点就是 96 行。如果文件里还有时间戳列读取时要跳过。我一般会在代码开头加一段数据检查逻辑% 检查数据维度 data xlsread(shuju.xlsx); if size(data, 2) 3 error(数据至少需要3列光伏、风电、负荷); end % 检测并剔除缺失值NaN nan_idx any(isnan(data), 2); if sum(nan_idx) 0 warning(发现%d行缺失数据已剔除, sum(nan_idx)); data(nan_idx, :) []; end % 归一化开关如果数据量纲差异大建议做归一化 % data(:, 2:4) data(:, 2:4) ./ max(data(:, 2:4), [], 1);这里检查了列数、缺失值、量纲三个维度。光伏出力的量纲是 MW 或 kW负荷可能是同一个量纲但不同数据源的量纲可能不一致。如果直接混在一起做场景生成负荷的数值尺度会主导距离计算导致削减结果几乎只由负荷变化决定光伏和风电的形态信息被淹没。碰到这种数据建议先做归一化到 [0,1] 区间削减完成后再反归一化回去。4.2 运行时序维度对场景生成的影响资源描述里说「每个时刻用拉丁超立方抽样函数抽取 200 样本」这句话意味着每个时刻是独立抽样的。独立抽样的最大问题是削减前的 200 个场景在时序上完全不相关时刻 t 处的最大出力样本和时刻 t1 处的最大出力样本大概率不是来自同一条「曲线」。削减后的场景在经济调度上往往表现为单时刻的出力波动而不是连续的趋势性波动。如果你希望场景呈现日出而作、日落而息的光伏特性就需要在抽样时加入时序相关性。常见做法是生成 200×T 个独立的标准正态样本后乘以一个 T×T 的相关系数矩阵的 Cholesky 因子使得同一场景在不同时刻之间的相关系数等于预设值时滞系数。但在毕设场景下独立抽样的简化处理是可接受的因为最终计算期望出力时时序相关性的影响会在平均中被部分抵消。4.3 削减算法运行过程中的矩阵维度细节MATLAB 里删除矩阵行然后继续循环的做法在 M 从 200 递减到 5 的过程中每次都要重新计算距离矩阵复杂度是 O(M²T) 乘以迭代轮数总轮数 195 次。200 个场景时距离矩阵是 200×200也就是 4 万个元素计算还是很快的。但如果把初始抽样数提到 2000距离矩阵变成 400 万元素每次迭代还要重新算一次速度就会明显变慢。工程上优化这个循环有几种做法% 进阶用向量化计算减少循环开销 % 初始化时一次性算出所有场景两两欧氏距离 all_dist zeros(N, N); for i 1:N diff sample_set - sample_set(i, :); all_dist(i, :) sqrt(sum(diff.^2, 2)); end % 后续迭代只需要查表 局部更新但这种做法的问题是样本集删行后序号发生变化查表索引要同步维护代码复杂度上升。对于毕设或项目演示直接删行重算就行200 个场景的规模完全不需要性能优化反而更不容易写错。4.4 削减后场景概率的归一化与校验按照上述合并逻辑每次累加概率后总概率恒等于 1不会出现概率和不为 1 的情况。但如果你修改过代码、加了条件判断分支最好在削减结束后做一个显式校验scene_prob scene_prob / sum(scene_prob); % 强制归一化 % 校验概率和非负性 assert(abs(sum(scene_prob) - 1) 1e-10, 概率和不为1); assert(all(scene_prob 0), 存在负概率); % 打印削减后的场景概率分布 disp(削减后场景概率); disp(scene_prob);概率归一化在后续做期望计算时不是必须的因为本来就归过一但加了 assert 断言能提前发现代码逻辑错误。我在实际调试时经常遇到的问题是误把削减代码里的场景概率初始化为零向量导致距离矩阵全是零削减循环直接停摆。这时候 assert 就会立即捕获。5. 场景削减结果的进阶验证、参数调优与 Python 复现对比5.1 削减质量的两个量化指标削减完成后不能只盯着曲线形状看「像不像」要有两个量化指标。第一个是削减前后的概率分布差异用 Wasserstein 距离来度量本质上就是概率距离削减算法里的 KD 距离。第二个是期望出力曲线与原始样本均值曲线的最大偏差和平均偏差。% 计算削减前后的统计量对比 orig_mean mean(samples_orig, 1); % 原始200场景的均值 orig_std std(samples_orig, 0, 1); % 原始200场景的标准差 red_mean scene_prob * reduced_scenes; % 削减后加权均值 red_std sqrt(sum(scene_prob .* sum((reduced_scenes - red_mean).^2, 2), 1)); % 均值绝对误差 mean_abs_err mean(abs(orig_mean - red_mean)); % 标准差偏差 std_rel_err mean(abs(orig_std - red_std) ./ (orig_std 1e-6)); fprintf(均值MAE %.4f\n, mean_abs_err); fprintf(标准差相对误差 %.2f%%\n, std_rel_err * 100);注意这里标准差的计算用的是场景概率加权而不是等权。削减后的均值曲线大概率跟原始均值曲线很接近但标准差往往偏小因为削减过程天然会丢弃掉远离均值的小概率高波动场景。如果你发现标准差相对误差超过 20%就说明 K5 的粒度不够需要提高到 7 或 8。5.2 参数敏感性抽样数 N 和方差随机因子对结果的影响N 从 200 变到 500削减后的 5 个场景形状不会有太大变化但概率值会微调。这是因为 LHS 在 N 增大时对分布尾部的覆盖更充分原本某些被忽略的极端场景可能获得更精确的概率估计。N 取太小比如 50削减结果会随随机种子明显抖动同一份数据跑两次会有不同的场景概率。方差随机因子的范围影响更大。把rand换成0.2 0.8 * rand场景的波动范围会被压缩削减后的 5 条曲线会更紧凑概率集中度更高。如果你希望场景里包含更多极端情况比如光伏骤降、负荷尖峰可以把上限提高到 1.5方差变成 1.5 倍均值尾部场景的场景会被削减算法保留下来但代价是期望曲线的平滑度变差。这在不同应用场景下没有绝对优劣关键是理解这个旋钮的物理含义——它控制的是场景集的离散程度。5.3 Python 复现的对照实现如果你后续要把这套场景生成算法整合到 Python 的强化学习或深度学习框架里比如做深度强化学习中环境随机性建模MATLAB 版本迁移到 Python 很直接。SciPy 里有现成的 qmc 模块提供拉丁超立方抽样函数import numpy as np from scipy.stats import norm, qmc # 生成 N 个 LHS 样本维度 d 3对应光伏、风电、负荷 sampler qmc.LatinHypercube(d3, seed42) sample sampler.random(n200) # 形状 (200, 3)取值 (0,1) z norm.ppf(sample, loc0, scale1) # 映射到标准正态 # 对每个时刻做均值方差映射 data np.loadtxt(shuju.csv, delimiter,) T data.shape[0] N 200 scenes np.zeros((N, T * 3)) # 展平后的场景矩阵 for t in range(T): mu data[t, :3] sigma np.random.rand(3) * mu # 方差随机比例 scenes[:, t*3:(t1)*3] mu sigma * z这里qmc.LatinHypercube生成的样本已经是分层均匀的不需要手动实现(i-rand)/N的逻辑。norm.ppf等同于 MATLAB 的norminv。削减部分用 while 循环重写时要注意 MATLAB 的norm(A-B)对应 Python 的np.linalg.norm(A-B)而场景矩阵用列表加 pop 操作要比 numpy 删行高效得多尤其原因是每轮只删一行Python list 的 pop 是 O(1) 操作而 numpy 的np.delete会复制整个矩阵200 个场景时没问题但规模上来后性能差距很大。5.4 两个字节能直接影响结果的工程细节第一个细节是随机种子。LHS 的样本生成依赖rand的返回值而 MATLAB 的 rand 状态每次启动都不同。如果不设置rng(2024)这类固定种子两次运行得到的场景曲线和概率会有很大差异这在毕设答辩现场可能造成「前后数据对不上」的尴尬。建议在脚本第一行加上固定种子。第二个细节是负值截断位置。光伏出力理论上不能为负但正态分布抽样一定会产生负值。如果你在抽样后立刻截断会让分布左尾的质量堆积到零上改变分布形态如果削减后再截断削减算法又会在负值区域计算距离。我的做法是在削减完成后对最终 5 个场景做max(0, scene)截断然后重新归一化场景概率。虽然处理不够严格但工程上可接受并且能避免负出力导致优化模型无解的问题。5.5 一个没有写在注释里的坑yuanhe.m 这类文件中常见的一个坑是矩阵维度隐式扩展导致的结果错位。比如抽样时samples_pv(:, t) mu_pv sigma_pv .* z_pv左边是 200×1右边 mu_pv 是标量sigma_pv 是标量z_pv 是 200×1没问题。但如果有人把代码改成从 Excel 的某一行向量里读取均值 mu 为 1×3 向量再和 z 做点乘就会触发 MATLAB 的隐式扩展生成 200×3 的矩阵然后赋值给 200×1 的列向量时报错。建议多用size(mu)和size(z)打印中间维度不要凭感觉写点乘。本文还有配套的精品资源点击获取