1. 从实际需求聊起为什么要做风光场景生成与削减做新能源电力系统优化的人大概率都撞过同一堵墙风电和光伏的出力曲线波动太大不确定性太强。拿一组实测数据直接塞进优化模型可能算出来的是“运气最好”或“运气最差”的情况结果压根不具备参考价值塞太多组数据进去模型规模又爆炸求解器直接罢工。所以业内默认的做法是**先把大量的可能出力情况“模拟”出来再用数学方法把那些长得差不多的场景合并掉留下一小撮有代表性的“典型日”。**这个过程就是标题里说的风光场景生成与场景削减。这套流程往上接随机规划、鲁棒优化往下接机组组合、经济调度、储能容量配置。你做任何涉及风光不确定性的MATLAB仿真大概率都绕不开这段“先蒙特卡洛、再概率距离削减”的组合拳。我最早接触这个思路是在做微电网容量配置的时候当时直接用一年的实测数据分365个场景往里灌结果模型规模大到连商用求解器都要跑十几个小时。后来换成先模拟再削减的处理方式几十个场景就能逼近原来的精度求解时间骤降到几分钟。这篇文章就把整个流程掰开揉碎讲一遍从数学原理到MATLAB代码实现再到实际运行中我踩过的坑全部写出来。2. 整体设计思路与两大核心模块拆解2.1 为什么是“蒙特卡洛生成 概率距离削减”的组合很多初学者会有个疑问既然场景削减能减少计算量那直接少生成一点不就行了比如一共算50个场景就不削减了直接把这50个丢进模型。这在理论上可行但实际操作中问题很大。蒙特卡洛法生成场景的本质是随机采样它的精度高度依赖于采样规模。你只抽50组随机数很可能把极端出力情况比如大风天气下风电满发、连续阴雨天光伏接近零出力漏掉但如果你抽5000组甚至10000组能覆盖各种工况却又没法直接用于计算。所以标准做法是两段式第一段用蒙特卡洛法大规模生成场景保证“可能性”都被覆盖第二段用概率距离削减法把海量场景压缩到目标数量保留最典型的几种情况并给每个典型场景配上发生概率。这个思路和拍电影的流程有点像——海选时广撒网找几百个候选人试镜后只留下少数几个风格差异最大的演员再根据市场偏好分配戏份权重。海选对应蒙特卡洛生成试镜对应场景削减。2.2 蒙特卡洛法与概率距离削减法的本质区别很多人把这两个概念混在一起其实它们分工完全不同。蒙特卡洛法是生成工具负责产出大量模拟场景。它依赖概率分布模型——风速通常用威布尔分布描述光照强度用贝塔分布描述然后通过随机抽样把这些分布变成一条条具体的出力时间序列。概率距离削减法是筛选工具负责把海量场景压缩成有限个代表场景。它的核心逻辑是两个场景长得越像就越没必要同时保留把多余场景的概率叠加到它最像的那个场景上这样既能减少场景数量又能保持整体概率分布特性不严重失真。打个比方蒙特卡洛是“制造硬币”概率距离削减是“找出最有代表性的几枚硬币把其他硬币的价值折算进去”。2.3 这个方法解决的最大痛点算力与精度的平衡场景削减这件事的价值归根结底一句话——让随机优化模型从“算不动”变成“算得动”同时让结果从“偏乐观”变成“可信”。我实测过一个算例用一组真实的风光数据生成3000个原始场景削减到30个典型场景。优化计算结果和用3000个场景直接算的期望值相比偏差控制在3%以内但求解时间只用了原来的1/10还不到。这样的性价比在实际工程中非常可观。对于研究新能源消纳、储能配置、微电网调度的同学来说这套方法基本属于标配技能。数学基础要求不高概率论基本概念加一点点最优化思想就能理解但要把它写成高效可复用的MATLAB代码还是有不少细节讲究。3. 核心算法详解与MATLAB实现3.1 风光出力概率模型威布尔分布与贝塔分布场景生成的底层是概率分布模型。风电和光伏的随机性来源不同模型也不同。风速模型风速本身通常用两参数威布尔分布来拟合概率密度函数为% 威布尔分布概率密度函数 % f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k) % k: 形状参数决定分布形态一般取2左右 % c: 尺度参数与平均风速相关风速转换成风电出力要看风电功率曲线。工程上常用简化的分段函数切入风速以下出力为0额定风速到切出风速之间出力饱和中间段近似线性或二次关系。光照强度模型光照强度G通常用贝塔分布描述概率密度函数为% 贝塔分布概率密度函数 % f(g) (g^(a-1) * (1-g)^(b-1)) / B(a, b) % a, b: 形状参数由历史数据的均值与标准差估计光伏出力与光照强度近似线性关系再乘上光伏板面积和光电转换效率即可。温度的影响一般模型中会忽略或者做简单修正。3.2 时序相关性不是所有随机数都能直接拼成场景蒙特卡洛生成风光场景时最容易犯的错误是每天的风速和光照随机数之间没有任何相关性生成出来的场景根本不符合实际。真实的风速存在明显的连续性和日周期性——今天风速大明天大概率也不小白天和夜晚的差异、季节差异都很大。如果忽视这些特性生成的场景就会像“白噪声”一样没有实际意义。我在代码里处理这块的思路是**先建立时序骨架再叠加随机扰动。**具体分三步根据历史数据的每小时均值构造一条“典型日曲线”作为基线。用AR(1)模型或在随机抽样时引入时段间协方差让相邻时刻的风速/光照变化平滑一些。把随机抽样的对象从“每个时刻独立抽取”改为“提取整日特征、再按特征生成曲线”。最终我实测下来最实用的做法是保持每小时独立抽样做原始场景是可以接受的因为削减阶段会把时间序列整体距离作为衡量标准最终保留下来的场景自然具备时序连续性。但如果你做的是需要精确评估爬坡率的调度问题那就必须在生成阶段加入时序相关性。3.3 场景削减算法快速前向削减 vs 同步回代削减概率距离削减法是一族方法最常用的是快速前向削减Fast Forward Selection和同步回代削减Simultaneous Backward Reduction。两者本质都是迭代法但思路相反同步回代削减初始场景集合S每次找一对“距离最近的场景”删掉其中一个把它的概率加到另一个身上。不断重复直到场景数达到目标值。快速前向削减从空集出发每次从未选中的场景里挑一个“使现有集合代表性最强”的场景添加进保留集合反复迭代直到场景数达到目标值。两种方法效果差不多但适用范围有区别。同步回代削减的计算量随场景数平方级增长适合场景数不是特别大的情况快速前向削减在迭代过程中计算更密集但往往能保留更多“极端场景”更适合需要覆盖尾部风险的决策问题。我自己的实践中同步回代削减在风电场景上用得多一些实现也直观适合讲解算法原理。下面给出核心代码框架。3.4 同步回代削减算法的MATLAB核心实现这一段是整篇博文的重点。我直接给出一份可以运行的结构化代码并逐步解释。function [reduced_scenes, reduced_probs] scenario_reduction(scenes, probs, target_num) % 输入: % scenes - 原始场景矩阵大小为 [T, S] % 行表示时间点小时列表示场景编号 % probs - 原始场景概率向量大小为 [1, S]一般初始化为等概率 1/S % target_num - 目标保留场景数量 % 输出: % reduced_scenes - 削减后的场景矩阵 [T, target_num] % reduced_probs - 削减后的概率向量 [1, target_num] [T, S] size(scenes); active ones(1, S); % 场景是否仍然保留1表示保留0表示已删除 probs probs(:); % 确保是行向量 % 迭代删除直到剩余场景数等于目标值 while sum(active) target_num min_dist Inf; del_idx 0; % 找到距离最近的一对场景 for j 1:S if active(j) 0, continue, end for i j1:S if active(i) 0, continue, end % 计算两个场景之间的欧氏距离 dist_val sqrt(sum((scenes(:,j) - scenes(:,i)).^2)); % 用概率加权的距离作为删除标准 weighted_dist probs(j) * dist_val; if weighted_dist min_dist min_dist weighted_dist; del_idx i; % 通常删除概率较小或索引靠后的场景 keep_idx j; end end end % 被删除场景的概率叠加到保留场景上 probs(keep_idx) probs(keep_idx) probs(del_idx); active(del_idx) 0; % 删除场景 end % 提取保留的场景与概率 keep_scenes find(active 1); reduced_scenes scenes(:, keep_scenes); reduced_probs probs(keep_scenes); reduced_probs reduced_probs / sum(reduced_probs); % 重新归一化 end这段代码的逻辑只有几层循环但它背后有两个细节值得大写特写**第一删除标准为什么用“加权距离”而不是纯粹的距离**因为不同场景发生的概率不同两个距离很近的场景如果其中一个概率特别大那么合并时“吞并”动作对整体分布的影响也大。用概率加权的距离可以优先删除“低概率且与其他场景相似的”目标让高概率场景尽量保留原样。**第二为什么要把被删场景的概率加到保留场景上**这是场景削减的本质——不是简单地把场景丢掉而是用“一个场景代表多个场景”的降维思想。被删场景原来占的概率必须体现在某个保留场景上否则削减后概率之和不为1后续优化模型就没法用了。3.5 基于概率距离的快速前向削减MATLAB实现快速前向削减的实现稍复杂但保留极端场景的能力更强。核心代码如下function [selected_scenes, selected_probs] fast_forward_selection(scenes, probs, target_num) % 快速前向削减法 % 每次从未选中场景中选取一个使得当前选中集合的“覆盖能力”最大 [T, S] size(scenes); selected false(1, S); % 选中标记 remain true(1, S); % 剩余场景标记 selected_probs []; % 选中场景概率 % 先选第一个场景与所有场景距离之和最小的场景 dist_mat pdist2(scenes., scenes.); % S x S 距离矩阵 total_dist sum(dist_mat, 1); [~, first_idx] min(total_dist); selected(first_idx) true; remain(first_idx) false; while sum(selected) target_num best_candidate 0; best_gain -Inf; for j find(remain) % 只遍历仍未选中的场景 % 计算将场景j加入选中集后的“收益” % 收益 所有未选中场景到选中集的最短距离的平均值越小越好 % 用负值表示便于比较 cand_sel selected; cand_sel(j) true; % 对每个未选中场景计算到选中集的最小距离 min_dists zeros(1, S); for i 1:S if selected(i), continue, end dist_to_set min(dist_mat(i, cand_sel)); min_dists(i) probs(i) * dist_to_set; end gain -sum(min_dists); % 负的总距离开销越大越好 if gain best_gain best_gain gain; best_candidate j; end end selected(best_candidate) true; remain(best_candidate) false; end selected_scenes scenes(:, selected); selected_probs probs(selected); selected_probs selected_probs / sum(selected_probs); end这份代码我刻意写成循环实现方便小白理解每一步在干什么。实际工程里如果处理上万场景强烈建议用向量化或矩阵运算优化否则计算量会非常大。后面我会给出性能优化建议。4. 完整实操流程从原始数据到削减后的场景集合4.1 数据准备与参数初始化为了演示完整流程我假设你要生成一个包含风电场光伏电站的微电网系统时间分辨率为1小时规划周期为24小时一天。目标是从2000个原始场景削减到30个典型场景。首先要定义地区特性参数% 风机参数 Wind.rated_power 1.5e6; % 额定功率 1.5MW Wind.cut_in 3; % 切入风速 m/s Wind.rated_wind 12; % 额定风速 m/s Wind.cut_out 25; % 切出风速 m/s % 光伏参数 PV.efficiency 0.16; % 光电转换效率 PV.area 6000; % 有效面积 m^2 PV.k 0.9; % 综合效率系数 % 威布尔分布参数根据历史数据拟合 Weibull.k 2.1; % 形状参数 Weibull.c 7.5; % 尺度参数 % 贝塔分布参数根据历史数据拟合 Beta.a 2.3; % 形状参数 Beta.b 2.1; % 形状参数 % 蒙特卡洛采样规模 num_samples 2000; % 原始场景数 num_target 30; % 削减目标场景数这里威布尔分布的 k 和 c 参数以及贝塔分布的 a、b 参数在真实工程中应该通过历史气象数据的极大似然估计拟合得到。演示用数据是我拍脑袋给的但代码逻辑完全一致。4.2 蒙特卡洛法生成风光出力场景% 预分配场景矩阵 wind_scenes zeros(24, num_samples); pv_scenes zeros(24, num_samples); % 波动系数用来控制每个时刻的随机波动幅度简化的时序相关性处理 fluctuation 0.3; for s 1:num_samples % 生成风速时间序列每小时一个随机值 wind_speed wblrnd(Weibull.c, Weibull.k, [24, 1]); % 简单时序平滑处理让相邻小时风速变化不会太突兀 for h 2:24 wind_speed(h) wind_speed(h-1) ... fluctuation * (wind_speed(h) - wind_speed(h-1)); end % 风速转风电功率 wind_power zeros(24, 1); for h 1:24 if wind_speed(h) Wind.cut_in || wind_speed(h) Wind.cut_out wind_power(h) 0; elseif wind_speed(h) Wind.rated_wind wind_power(h) Wind.rated_power * ... (wind_speed(h) - Wind.cut_in) / (Wind.rated_wind - Wind.cut_in); else wind_power(h) Wind.rated_power; end end wind_scenes(:, s) wind_power; % 生成光照强度时间序列0-1标幺值 daylight_mask (7 0:23) (0:23) 18; % 7点到18点有光照 solar_irr zeros(24, 1); for h 1:24 if daylight_mask(h) solar_irr(h) betarnd(Beta.a, Beta.b); if solar_irr(h) 1, solar_irr(h) 1; end end end % 光照转光伏出力 pv_power PV.efficiency * PV.area * PV.k * solar_irr * 1000; % W pv_scenes(:, s) pv_power; end % 合并成总场景矩阵 total_scenes wind_scenes pv_scenes; % 初始场景概率等概率 initial_probs ones(1, num_samples) / num_samples;这段代码里有几个可能会被忽略但实际至关重要的细节白天窗口的处理。光照在夜间就是0不能把贝塔分布的随机数往夜里塞。我简单用了一个固定窗口[7, 18]实际项目里要根据纬度和季节动态调整。4.3 调用削减算法并分析结果% 削减到30个典型场景 [reduced_wind, red_wind_probs] scenario_reduction(wind_scenes, initial_probs, num_target); [reduced_pv, red_pv_probs] scenario_reduction(pv_scenes, initial_probs, num_target); % 也可以对总出力场景直接削减 [reduced_total, red_total_probs] scenario_reduction(total_scenes, initial_probs, num_target);削减完之后建议做两件事来验证结果质量。第一件事比较削减前后的期望出力曲线。% 削减前期望2000个场景均值 mean_expect mean(total_scenes, 2); % 削减后加权期望 weighted_mean sum(reduced_total .* red_total_probs, 2); % 对比误差 error_curve abs(weighted_mean - mean_expect); max_error max(error_curve); fprintf(最大期望误差: %.2f kW\n, max_error);如果误差在5%以内说明削减质量不错如果误差偏大可以适当增加保留场景数量或者检查原始场景生成时分布参数是否合理。第二件事可视化对比。把原始场景集合和削减后的场景集合画在同一张图上直观感受削减效果。通常你会看到30条曲线能很好地勾画出2000条曲线的“包络带”。4.4 随机数与可复现性的处理蒙特卡洛法最大的坑就是随机数种子的问题。你没设置好种子每次运行结果都不同论文里的图表没法复现审稿人直接质疑你的结果。我习惯在每次运行前固定随机数种子rng(2024); % 固定随机数种子这样同一台机器上每次都生成一样的场景。如果你在做算法对比实验务必固定种子保证不同方案在同一组场景下比较否则结果差异可能来自随机性而非算法优劣。4.5 直接可用的整合脚本框架最后给一个整合版的 publish-friendly 脚本结构方便你直接套用%% 主脚本风光场景生成与削减 clear; clc; close all; %% Step 1: 定义参数 rng(2024); % ... 参数与4.1节一致 ... %% Step 2: 蒙特卡洛生成原始场景 % ... 与4.2节一致 ... %% Step 3: 场景削减 % ... 调用削减函数 ... %% Step 4: 结果可视化与误差分析 figure; plot(total_scenes, Color, [0.8 0.8 0.8]); hold on; plot(reduced_total, r, LineWidth, 1.5); xlabel(时间 (h)); ylabel(出力 (kW)); title(场景削减效果对比灰原始2000场景红削减后30场景);5. 常见问题与排查技巧实录这部分是我最想写的因为这些坑我基本都趟过。很多问题教科书上根本不提只有真跑代码才会遇到。5.1 削减后概率不归一或出现异常值现象削减后所有概率加起来不是1或者某些概率出现负值。原因分析多半在代码某处把概率覆盖了或者把已经被删除的场景又重新纳入计算。排查方法用disp打印每一步迭代后的sum(probs)检查是否始终为1。如果中间某一步不满足基本可以确定是索引错误导致场景重复或概率赋值逻辑有问题。我的建议始终在函数末尾做一次归一化处理作为兜底。但这只是补救真正的排查还得靠逐步打印。5.2 削减后场景少了一条“大风光天”现象削减后的场景集里找不到风电满发持续一天的极端场景导致储能配置结果偏保守。原因分析这是同步回代削减的一个固有缺陷——它倾向于删除远离其他场景的“极端点”因为标准化距离衡量下极端场景往往和谁都离得远在“找最近对”的过程里一直轮不到被删除但极端场景数量少、概率低可能在迭代初期就被一个高概率场景吞并了。解决方案有两个改用快速前向削减法它天然倾向于保留极端场景。在削减前做场景聚类预处理先把2000个场景按K-means聚成50类每类里挑一个代表场景并赋予该类总概率再做一次削减。这个思路能显著减少极端场景被误删的概率。5.3 双参数分布拟合不准确生成场景失真现象生成的风速曲线整体偏低/偏高和实际历史数据对不上。原因分析威布尔分布的两个参数直接决定生成数据的样子。如果 k 和 c 没拟合好后面一切都白搭。我实测有效的拟合方法用MATLAB自带的fitdist函数对历史风速数据做最大似然估计% 假设 wind_hist 是历史小时风速数据列向量 pd fitdist(wind_hist, wbl); Weibull.k pd.B; % 威布尔形状参数 Weibull.c pd.A; % 威布尔尺度参数光照贝塔分布也可以用类似方式% 假设 solar_hist 是历史光照数据标幺化到0~1 pd fitdist(solar_hist, beta); Beta.a pd.a; Beta.b pd.b;如果fitdist报错或拟合效果差常见原因是历史数据里0值太多夜晚导致贝塔分布无法收敛。解决办法是只提取光照大于0的时刻来拟合成形状参数在生成场景时再叠加上“白天窗口”和“天气概率”效果比硬拟合完整数据好得多。5.4 循环计算太慢2000个场景跑半小时现象scenario_reduction函数跑起来慢得离谱尤其当S2000时双重循环计算距离矩阵每轮迭代都要重新计算复杂度接近O(N^3)完全不能忍。解决方案使用pdist2提前一次性算好所有场景两两间的距离矩阵之后每次迭代只需查表不再重复计算。用向量化逻辑替代内层循环用布尔索引直接找出active的场景对。如果场景数量上万考虑改用K-means预聚类把规模降到500以内再削减。优化后的距离矩阵版本% 提前计算距离矩阵这是最大的优化点 dist_mat pdist2(scenes., scenes.); while sum(active) target_num min_dist Inf; del_idx 0; keep_idx 0; for j find(active) for i 1:S if i j || ~active(i), continue, end weighted_dist probs(j) * dist_mat(i, j); if weighted_dist min_dist min_dist weighted_dist; del_idx i; keep_idx j; end end end ... end这样改动后2000个场景的削减时间通常能从几分钟降到几秒效果立竿见影。5.5 削减数量选多少合适这可能是大家最纠结的另一个问题——到底削减成多少个场景合适我的经验公式看你优化模型的变量规模。如果你的随机规划模型每个场景带一套连续变量和若干0-1变量那场景数每增加10个求解时间可能翻倍。所以要从计算资源和精度需求两头挤。实测一般规律保留场景数精度误差求解时间混合整数线性规划模型108%~10%十几秒303%~5%2~3分钟501%~2%20~30分钟1001%数小时通常取20~30个是性价比较好的区间。如果你需要覆盖极端天气大风光日、无风无光日可额外手动添加几个“工程典型场景”再和削减后的场景合并保证优化结果不因极端场景缺失而失真。5.6 风、光场景是否要单独削减还是合并削减这是一个经常被问起但在文献里很少细讲的细节问题。从数学角度风、光出力可以合成一个总出力向量来削减也可以各自单独削减再配对。场景削减的核心是场景间的距离而距离定义方式直接影响削减结果。我的实测体会如下单独削减风场景适合做风电为主、光伏为辅的系统或者需要对风电和光伏分别建模的场景比如储能配置要分别看两部分波动。合并成总出力再削减适合做微网总体调度追求的是总出力画像的典型性。合并后的削减能保留风光互补特性但可能牺牲对单一电源极端情况的覆盖。最稳妥的建议是两套方案都跑一版比较削减后的期望出力曲线与原始期望的误差选误差更小的那个。如果时间不够优先选单独削减并配对因为后续可以把风光间的相关系数也注入配对过程调整空间更大。6. 我实测下来的一些补充心得最后一段说说我在实际项目里用这套方法的体感。最值得注意的一个点**场景削减不是一个“做完就完了”的步骤它和后面优化模型的目标函数有耦合。**比如我的目标是最小化系统总成本那么削减时用欧氏距离衡量场景相似性本质上假设了“出力曲线差异和成本差异成正比”。但储能系统的成本和非线性充放电特性会让这个假设不完全成立——两个出力曲线看起来像的场景配上不同储能策略后成本差异可能很大。如果目标函数里含有强非线性项比如电池寿命衰减、多阶段决策建议在削弱阶段改用动态时间弯曲距离或考虑成本函数的一阶近似距离不要固定死在欧氏距离上。这一点在很多教程里不会讲但对最终方案的可行性影响不小。还有一个小习惯削减完一定要检查代表性场景的累计概率分布函数和原始场景的差异。直接对比均值误差是达标的但尾部概率分布可能已经严重失真。负荷侧关注可靠性指标的务必检查削减后场景集中是否保留了风电/光伏低于20%额定的连续天场景这类场景直接决定切负荷风险。实际操作中我做了一套可视化脚本把原始场景的累计概率分布和削减后场景的累计概率分布画在一起直观观察两曲线有没有明显偏移。如果尾部对不上我会手动调低极端场景的“被合并优先级”比如在距离计算中给极端出力场景的加权系数翻倍虽然这在纯数学上不太严谨但工程上非常好用。这套代码我后来封装成了一个通用的wind_solar_scenario_analysis工具把分布拟合、数据导入、场景生成、削减、质量分析、可视化全部串在一起。整个过程说不上复杂但每一步的细节都决定了最终结果能不能被审稿人或工程验收认可。如果你也是做新能源不确定性建模的建议把这套流程彻底吃透自己动手跑一遍比只看代码要记忆深刻得多。我踩过的坑你不需要再踩一遍照着上面的排查思路能省下好几个星期的时间。