企业官网建设流程全解析

简介：这份程序资源用于电力系统或综合能源系统中源荷场景生成，适合做毕业设计、学术仿真及算法研究的读者。程序基于拉丁超立方抽样方法，每个时刻抽取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 = i+1: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 变量和连续变量数量会成倍增加，求解速度下降。所以 K=5 是精度和求解效率的折中。如果你跑出来的削减结果里某个场景的概率超过 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 处的最大出力样本和时刻 t+1 处的最大出力样本大概率不是来自同一条「曲线」。削减后的场景在经济调度上往往表现为单时刻的出力波动，而不是连续的趋势性波动。

如果你希望场景呈现日出而作、日落而息的光伏特性，就需要在抽样时加入时序相关性。常见做法是生成 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%，就说明 K=5 的粒度不够，需要提高到 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(d=3, seed=42) sample = sampler.random(n=200) # 形状 (200, 3)，取值 (0,1) z = norm.ppf(sample, loc=0, scale=1) # 映射到标准正态 # 对每个时刻做均值方差映射 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:(t+1)*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)打印中间维度，不要凭感觉写点乘。

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