拉丁超立方抽样与场景削减：风光荷场景预测实战解析

简介这份资源围绕风光出力与负荷场景预测提供基于拉丁超立方抽样和样本削减的完整MATLAB实现适合能源电气领域的研究生、工程师及相关从业者用于处理多维气象参数下的不确定性建模。压缩包整体仅1KB包含1个m脚本WT_PV_Load_Scenario.m代码流程清晰且配有注释涵盖数据预处理、LHS代表性样本抽取、样本削减以及预测模型训练与误差评估等核心环节。已有2593人学习下载说明该方案在同类场景中具备一定参考热度。通过学习此代码可快速掌握如何结合拉丁超立方抽样与样本削减方法在保证预测精度的同时降低计算复杂度进而更高效地分析风、光出力与负荷的多场景波动规律提升可再生能源并网研究的实用性与可复现性。1. 风光出力与负荷场景预测为什么用拉丁超立方抽样做可再生能源并网分析的人恐怕都有过这种体验想评估某个微网在不同气象条件下的运行状态于是用蒙特卡洛随机抽样生成风速和光照场景抽了几千组数据结果发现场景扎堆——极端高温天和静风日出现的频率和实际气候规律对不上直接导致后续调度策略在极端工况下翻车。这不是算法不够快而是随机抽样的空间覆盖能力有先天短板。拉丁超立方抽样LHS配合样本削减技术能把风光出力和负荷的不确定性空间压缩到一个可计算的规模同时保留原始分布的关键统计特征。这套方法对做电力系统随机优化、储能容量配置、微电网调度策略验证的工程师和研究人员都实用尤其适合处理二维以上、变量相互耦合的气象-负荷场景生成问题。下面把原理、实现和踩过的坑逐一拆开讲。2. 拉丁超立方抽样多维空间均匀采样的实现与代码2.1 为什么蒙特卡洛在高维场景下不够用蒙特卡洛抽样的逻辑是让随机数自然地在概率空间里散布样本量足够大时统计特征会收敛。但实际操作里有个尴尬当气象变量不止一两个而是风速、辐照度、温度、湿度、负荷波动因子等五六个维度时蒙特卡洛要达到同样的覆盖质量需要的样本量呈指数级增长。而且随机数天生具备“聚集性”——某些局部区域样本特别密另一些区域特别稀。这意味着按蒙特卡洛抽一千个场景做电网潮流计算真正覆盖到高辐照低负荷这种极端组合的样本可能只有几个概率分布被严重低估。LHS的思路完全不同它先把每个维度的累计概率分布均匀切成 N 段每段只取一个样本点然后把这些点随机组合成 N 个场景。这样做保证每个维度上样本的边际分布是均匀的不会出现某个区间过密或过空的情况。换句话说同样的样本量LHS对空间的覆盖更彻底。2.2 手写一个 LHS 抽样函数MATLAB 的统计工具箱提供了lhsdesign和lhsnorm但手写一个会让人对分层逻辑有更清楚的认识。核心就三步分段、段内随机取值、打乱排列。function X lhs_manual(nSamples, nDims, lb, ub) % 手写拉丁超立方抽样 % nSamples: 需要的场景数每维分多少段 % nDims: 气象变量的个数 % lb, ub: 每个变量的物理下界和上界维度是 1 x nDims X zeros(nSamples, nDims); for d 1:nDims % 1) 把该维度的取值范围切成 nSamples 个等宽区间 edges linspace(lb(d), ub(d), nSamples 1); % 2) 在每个区间内随机取一个点 for i 1:nSamples X(i, d) edges(i) (edges(i1) - edges(i)) * rand(); end % 3) 打乱顺序避免维度之间形成人为的排序相关性 X(:, d) X(randperm(nSamples), d); end end参数上要注意nSamples既是场景数也是每个维度的分段数一般场景数取 200 到 1000lb和ub的取值直接影响抽样结果风速的上下界通常取该地点历史数据的 2% 和 98% 分位数而不是硬性的 0 和极大值否则会抽出大量不切实际的场景。第 3 步打乱顺序非常关键如果所有维度都按区间顺序排列场景矩阵的每一行会呈现明显的升序规律统计上等同于人为制造了正相关。2.3 内置函数与相关性注入手写版本适合理解逻辑实际项目里我更倾向于用内置函数省时间。lhsdesign的Criterion参数可以让样本的填充程度进一步优化maximin表示最大化样本点之间的最小距离适合稀疏场景correlation则适合需要指定变量相关性的情况。% 固定随机种子保证结果可复现 rng(2024); % 生成 500 个样本、4 个维度的 LHS 设计 X_unit lhsdesign(500, 4, Criterion, maximin, Iterations, 50); % 映射到实际的物理范围lb 和 ub 是 1x4 的向量 X_phys X_unit .* (ub - lb) lb;但这里有个埋得很深的坑lhsdesign生成的是 0 到 1 之间的均匀设计如果直接把各列映射到不同变量的物理范围得到的场景矩阵在变量间是相互独立的。而实际的风速和辐照度往往存在负相关——晴天通常风小阴天风大但不一定。独立抽样会生成大量“大风大太阳”或“静风阴天”的不合理组合。解决办法是先让样本服从标准正态分布再通过 Cholesky 分解注入目标相关矩阵最后映射回均匀分布。% 目标相关系数矩阵风速 vs 辐照度 -0.45 R_target [1.0, -0.45, 0.15; -0.45, 1.0, 0.1; 0.15, 0.1, 1.0]; L chol(R_target, lower); % 把均匀 LHS 样本转为标准正态 Z norminv(X_phys_normalized); % 注入相关性 Z_corr Z * L; % 转回 [0,1] 后再映射物理范围 X_corr normcdf(Z_corr);这个小步骤属于“不做后面全是坑”的典型。样本削减阶段如果用欧氏距离衡量场景相似度相关性没注入的话削减出来的代表场景在统计特征上是错的。后文第 5 章会详细展开。3. 样本削减从上千场景到几十个代表场景的取舍逻辑3.1 为什么削减是必要的风光出力场景生成的最终目的往往是喂给优化模型做计算。一个含 96 个时段的日场景每多加一个场景混合整数规划或者随机动态规划的求解时间可能翻几倍。一千个原始场景算起来非常吃力几十个场景却能跑动。样本削减就是在保留原始概率分布信息的前提下挑出一组最小的代表性场景集合。这里的理论支撑是概率测度之间的距离工程上常用的是 Kantorovich 距离通俗理解是“把一个概率分布搬运成另一个概率分布的最小代价”。削减的目标就是找一个子集让子集分布与原始分布的 Kantorovich 距离最小。3.2 前向选择与后向消除的实现实际操作中削减有两种经典策略前向选择和后向消除。前向选择从空集出发每次从剩余场景里挑一个加入集合使得新集合与原始集合的 Kantorovich 距离下降最多后向消除反过来从全集出发每次删一个最不重要的场景。前向选择适合最终场景数远小于原始数的场景计算量更可控。function [sel_idx, new_prob] scenario_reduction_forward(Scen, p_orig, n_keep) % Scen: nDim x nScen 的原始场景矩阵每列一个场景 % p_orig: 原始概率向量通常均匀为 1/nScen % n_keep: 最终保留的场景数 nScen size(Scen, 2); selected []; remaining 1:nScen; for k 1:n_keep best_dist inf; best_candidate remaining(1); for i remaining trial [selected, i]; % 计算当前集合与原始集合的 Kantorovich 距离近似 d calc_kan_dist(Scen(:, trial), Scen, p_orig); if d best_dist best_dist d; best_candidate i; end end selected [selected, best_candidate]; remaining(remaining best_candidate) []; end % 概率重分配被剔除场景的概率归并到最近的保留场景 new_prob redistribute_probs(Scen, selected, p_orig); sel_idx selected; endcalc_kan_dist的具体实现对每个原始场景找到当前集合中离它最近的保留场景把二者欧氏距离乘以原始概率再相加即为近似 Kantorovich 距离。这个逻辑做两个循环就能跑通。redistribute_probs则是把每个被剔除场景的概率累加到距离它最近的保留场景上。3.3 削减后的概率重标定这一步是初学者最容易忽略的。很多人削减完直接把保留场景的概率当成均匀分布也就是每个场景概率取1/n_keep严格来说这不正确。原因是保留场景的位置分布并不均匀有些场景周围被合并掉的原始场景多它的概率应当更大。正确的做法是按最近邻归并function new_prob redistribute_probs(Scen, sel_idx, p_orig) nScen size(Scen, 2); new_prob zeros(length(sel_idx), 1); for i 1:nScen if ismember(i, sel_idx) % 保留场景自身概率 new_prob(sel_idx i) new_prob(sel_idx i) p_orig(i); else % 被剔除场景归并到最近保留场景 dists sqrt(sum((Scen(:, i) - Scen(:, sel_idx)).^2, 1)); [~, nearest] min(dists); new_prob(nearest) new_prob(nearest) p_orig(i); end end % 归一化确保概率和为 1 new_prob new_prob / sum(new_prob); end一个直观的经验是削减后场景概率的最大值和最小值能相差 3 到 5 倍。如果嫌手动实现太费劲也可以用现成的场景削减代码包但务必确认其概率重标定逻辑是否正确。4. 风光荷场景生成完整流程从原始气象数据到可计算的场景集4.1 数据清洗与归一化的顺序问题场景生成的第一步是拿到历史气象和负荷数据一般是逐小时或逐 15 分钟的数据。清洗时要去掉异常跳变比如风速传感器故障导致的零值、填充因为通信中断产生的缺失段。这里有一个容易被忽略的细节归一化要在 LHS 抽样之前做还是在抽取之后做。我的建议是先做异常值截断和线性插值然后直接 LHS 抽样最后对抽样结果做归一化映射。如果先归一化再抽样边界值可能会被压缩掉导致极端场景丢失。% 读取历史数据格式为 [风速, 辐照度, 温度, 负荷] 四列 raw_data load(site_hist_data.mat); wind raw_data.wind_speed; % m/s solar raw_data.solar_irrad; % W/m^2 temp raw_data.temperature; % °C load_p raw_data.load; % kW % 1) 用 2%-98% 分位数截断异常值 p_low prctile(wind, 2); p_high prctile(wind, 98); wind(wind p_low | wind p_high) NaN; % 2) 线性插值填充缺失 wind fillmissing(wind, linear);fillmissing的linear方法对渐变的气象时间序列够用但如果缺失段超过 6 小时建议改用spline或者直接丢弃该段避免插值出不符合物理规律的数据。4.2 风速到风电出力的功率曲线映射LHS 抽样得到的是风速的原始物理量要变成场景矩阵里可用的风电功率需要一个风机功率曲线函数。常用的是分段函数模型function P_wt wind_to_power(v, P_rated, v_in, v_r, v_out) % 典型风机功率特性曲线 % v_in, v_r, v_out: 切入、额定、切出风速 % P_rated: 额定功率 (kW) P_wt zeros(size(v)); % 风速在切入和额定之间按三次方关系上升 idx (v v_in) (v v_r); P_wt(idx) P_rated .* (v(idx).^3 - v_in^3) ./ (v_r^3 - v_in^3); % 达到额定风速后满发 idx (v v_r) (v v_out); P_wt(idx) P_rated; % 超过切出风速停机保护出力为零 end这段代码判断了三个区间。参数按常见机型取v_in3,v_r12,v_out25m/sP_rated按实际风机铭牌填。需要注意功率曲线的非线性会改变抽样分布即便风速是均匀分层抽样的经过三次方映射后低风速区间的功率差异会被压缩呈现大量的低出力样本。4.3 光伏出力建模与负荷场景生成光伏出力与辐照度、组件温度的耦合关系可以简化为线性模型function P_pv solar_to_power(G, T_cell, P_stc, beta_temp) % G: 入射辐照度 (W/m^2) % T_cell: 光伏板温度 (°C) % P_stc: STC 条件下额定功率 (kW) % beta_temp: 温度系数通常 -0.004 到 -0.005 /°C P_pv P_stc .* (G ./ 1000) .* (1 beta_temp .* (T_cell - 25)); % 辐照度低于阈值时出力归零 P_pv(G 50) 0; end负荷场景则要叠加上负荷的周期性特征。典型做法是先构造一个 24 点的标准负荷形态再乘以服从正态分布的随机扰动因子function P_load_scen load_scenario(base_load, std_shape, mu, sigma, nScen) % base_load: 日最大负荷基准值 (kW) % std_shape: 24x1 的归一化典型负荷曲线 % mu, sigma: 负荷扰动因子的均值与标准差 rng(default); noise 1 mu sigma * randn(nScen, 24); P_load_scen base_load * std_shape .* noise; end典型曲线std_shape应当取多个历史工作日的归一化负荷均值而不是某个单日数据否则场景会带上单日特有的随机波动。4.4 场景矩阵的组装与削减所有模型都建立之后把每个场景组织成一列然后调用第 3 章的削减函数将 500 个原始场景缩减到 30 到 50 个代表性场景。组装时每一行是一个时刻每一列是一个场景行顺序建议是 [风电功率 24 点; 光伏功率 24 点; 负荷 24 点] 或按时刻交错排列取决于后续优化模型的读数据习惯。我习惯按时刻顺序排列方便画 24 小时场景曲线族。nScen 500; nKeep 40; % 生成 24 时段的场景矩阵 Scen_wind zeros(24, nScen); Scen_pv zeros(24, nScen); Scen_load zeros(24, nScen); for s 1:nScen v_scen lhs_manual(1, 4, lb, ub); % 简化示意实际应一次性生成 Scen_wind(:, s) wind_to_power(v_scen(1) temp_series, P_rated, 3, 12, 25); Scen_pv(:, s) solar_to_power(v_scen(2), v_scen(3), P_stc, -0.004); Scen_load(:, s) load_scenario(base_load, std_shape, 0, 0.05, 1); end % 拼接为 72 x nScen 的完整场景矩阵 Scen_all [Scen_wind; Scen_pv; Scen_load]; % 削减 [sel_idx, new_prob] scenario_reduction_forward(Scen_all, ones(1,nScen)/nScen, nKeep);5. 避坑LHS 与场景削减实战中的五个典型问题5.1 抽样场景的风光相关性为零导致电网平衡分析失真现象用 LHS 独立抽样生成风速和光照生成的场景里出现了大量“大风天 强光照”的组合。放在实际电网里这种场景基本不存在算出来的弃光率或备用容量需求严重偏大。原因LHS 本身没有考虑变量间的统计相关性各维度独立抽取后随机组合等于假设风速和光照互不影响。解决在抽样后加一步 Cholesky 分解注入目标相关性或者使用lhsdesign的Criterion设置为correlation选项。目标相关矩阵用历史气象数据逐月统计得出不要拍脑袋填。负相关系数取 -0.4 到 -0.5 通常是合理的具体要看当地气候。5.2 削减后极端场景被全部丢掉现象削减后的 30 个场景都集中在“温和天气”附近系统的最大爬坡需求被严重低估。原因Kantorovich 距离最小化本质上是概率分布的保守拟合极端场景在原始概率中占比低被合并时对距离的贡献小因而优先被删除。解决削减前先做场景分类比如按“极端高温日”“极端低温日”“大风日”“静风日”分组每组单独削减再组合或者给极端场景附加一个惩罚权重让它们的概率在距离计算中放大 1.5 到 2 倍。这个技巧在储能容量配置里尤其重要极端无风夜连着高负荷日的场景往往决定储能装机上限。5.3 随机种子不固定同一组数据每次算出不同的结果现象同一份历史数据、同样的场景数连续跑两次程序得出的削减场景不一样连优化结果都跟着变。原因LHS 抽样和场景削减里都用了随机数rand和randperm默认每次启动都重新播种。解决在程序入口固定种子。rng(2024); % 固定为某一年数字可换任意整数若做不同场景数的对比实验不要只固定一次种子应该在每个实验中显式调用不同种子避免所有实验共享同一套随机序列。5.4 削减时用了全时段欧氏距离时序偏移被严重惩罚现象两条负荷曲线形状几乎一样只是午高峰时间差了半小时计算出来的“距离”却比一条完全平坦的曲线还大削减后保留了后者导致负荷形态失真。原因欧氏距离直接逐点相减时间错位会被当成两倍的数值差异而不是形状差异。解决改用动态时间弯曲DTW距离作为削减度量或者对欧氏距离加一个时间窗口平滑把每个点的值替换为前后 3 点的加权平均再做距离计算。DTW 的计算成本高对 96 时段的场景矩阵尤其明显一个工程化的折中方案是把曲线先做 12 点或 6 点降采样再算欧氏距离。降采样后时段点错位的影响被逐级吸收速度损失可接受。5.5 归一化不小心把负值负荷变成了零现象场景削减完成后检查时发现某些负荷场景出现全零段概率还不低。原因用了mapminmax这类全局归一化函数对含零值的向量操作时把负值强制截断。或者在做数据清洗时用max(0, x)的方式处理了传感器负值读数。解决负荷数据不应该做全局最小最大值归一化而应该按均值-标准差标准化且标准化系数只从训练段统计。对光伏辐照度的零值要单独保留不与负荷列混在一起做标准化。6. 从场景集到区间预测一个把场景转为置信区间的实用技巧场景削减得到的代表性场景除了丢给优化模型求解还能直接用于生成预测区间。这里分享一个我惯用的做法把削减后的场景按概率从低到高排列计算每个时段的累积分布直接读取 10%、50%、90% 分位数作为预测区间。% 假设 Scen_reduced 是削减后的 72 x nKeep 场景矩阵 % new_prob 是对应概率已经是重标定后的 nKeep length(new_prob); % 对每个时段计算加权分位数 quantiles [0.1, 0.5, 0.9]; interval zeros(24, 3); for t 1:24 % 提取该时段所有场景的值按概率排序 values [Scen_reduced(t, :); Scen_reduced(t24, :); Scen_reduced(t48, :)]; weights repmat(new_prob, 3, 1); % 按值排序并计算累计概率 sorted_vals sort(values(:)); sorted_weights weights(:); [sorted_vals, idx] sort(values(:)); sorted_weights new_prob(idx); cum_prob cumsum(sorted_weights) / sum(sorted_weights); % 线性插值读取分位数 for q 1:3 interval(t, q) interp1(cum_prob, sorted_vals, quantiles(q), linear); end end这套方法比直接用原始场景分位数更稳因为场景削减已经把概率密度平滑过了区间不容易出现跳变。还有一个习惯值得分享在做随机优化时通常要对比“确定性预测”和“随机优化”两组结果确定性预测用的是quantile0.5的场景随机优化则用全部保留场景。两个方案的成本差就能直观量化“不确定性带来的额外运行成本”这也是场景分析在电力系统应用里最有说服力的输出。把场景削减用到滚动预测中时建议每小时更新一次气象数据时重新跑一遍 LHS 和削减但随机种子保持固定这样相邻时段的场景不会出现剧烈跳变。我自己经历过一次尴尬为了追求新鲜数据每次重置了随机种子结果两个相邻时段削减出来的场景形态完全不同调度指令跟着来回切换被折腾得不轻。从那以后我每次跑滚动场景分析都会把rng设置和气象数据更新时间戳写在同一行日志里保证逻辑可追溯。这个方法不算复杂但确实省掉了很多半夜排查问题的辛苦希望帮到你。本文还有配套的精品资源点击获取