我们在做电力系统随机优化或者微电网规划的时候经常要面对一个很实在的问题风光出力怎么描述直接拿一整年时序数据丢进去计算量吃不消只取典型日又怕丢掉极端情况。我自己最早是被一个“要生成500个风电场出力场景”的需求逼着去折腾蒙特卡洛和场景削减的搞完才发现这套流程其实可以通用到风电、光伏甚至负荷预测不确定性建模上。这个项目的核心就是两件事先用蒙特卡洛“凭空造”出大量符合统计规律的风光出力场景再用基于概率距离的快速削减算法把几百上千个场景砍到几十个同时保证削减后的场景集合还能代表原始分布。如果你正被“场景怎么生成”“场景怎么削”卡住这篇内容值得收藏。1. 为什么需要生成与削减场景从随机性问题说起1.1 单条时序曲线不够用我们必须面对“不确定性”电网里风光接入以后调度、规划、可靠性分析都要考虑“明天风电到底发多少”。传统做法是拿历史数据均值或者典型日曲线当输入但实际出力波动很大极端天气、云层移动都会让实际值和预测值差出一大截。于是大家开始用“场景法”用一组带概率的确定性曲线去近似随机过程。每种可能的出力情况叫一个场景所有场景带上概率就是场景集。但这里有个麻烦如果直接用蒙特卡洛或者随机抽样生成场景要想把分布描述准确往往要生成上千个场景。把这些场景全部放进随机优化模型比如经济调度、储能优化配置每个场景对应一组约束和变量求解器直接就跑不动了。所以必须做“削减”用少量场景保留原始概率分布的关键信息。1.2 蒙特卡洛负责“生成”概率距离负责“削减”项目标题里有两个关键词蒙特卡洛和概率距离快速削减算法。前者负责生成足够多样本覆盖概率空间后者负责在样本集合中挑出一部分代表同时重新计算每个代表的概率。这两步连起来就是一套完整的“场景不确定性刻画”工具链。很多人会问为什么不直接用场景聚类K-means聚类当然也能削减但传统聚类往往只考虑场景之间的“空间距离”忽略了概率权重而且维数高的时候收敛很慢。概率距离削减方法则把概率分布的距离作为衡量标准削出来的场景集合更贴近原始分布。其中常见的是Kantorovich距离还有Wasserstein距离的变体。这套方法在SCENRED工具GAMS里那个场景削减工具里用得很多本质上就是后向削减或者快速前向选择。2. 场景生成基本原理与建模2.1 风速和光照的统计特征怎么描述做风光场景生成第一步是给随机变量选概率分布。风电出力通常先转换成风速风速的经典分布是双参数的Weibull分布概率密度函数是f(v; k, c) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)其中k是形状参数c是尺度参数。k一般在1.5到3之间c和平均风速有关。光伏出力则和光照辐照度强相关辐照度通常用Beta分布描述参数可以通过历史数据的均值和方差去估计。实际项目中我没法拿到每个风电场的实测风速分布参数所以代码里可以做一个简化允许用户直接输入风速的威布尔参数或者用历史数据来拟合。如果完全没有历史数据也可以用经验值比如c8.5m/s、k2.0。2.2 从风速到风电出力再到时序场景风速样本抽出来后要转换成风电出力。这里需要风机功率曲线工程上常用分段线性化或者用指数逼近的公式P_w(v) 0, v v_in 或 v v_out P_rated * (v - v_in) / (v_rated - v_in), v_in ≤ v v_rated P_rated, v_rated ≤ v ≤ v_out类似地光伏出力可以用辐照度和温度修正但简化场景生成时可以只考虑辐照度到出力的线性关系最多加一个光照转换效率。重要的不是模型多精细而是保持时序相关性。场景不是单点抽样而是一条时间曲线。直接每个时刻独立抽样生成的场景会非常“毛糙”相邻时刻出现剧烈跳变这不符合实际风光的平滑变化特征。所以代码里一般要引入时序相关性比如用自回归模型AR(1)生成风速序列或者先抽一个基础风速趋势再加噪声。简单版本可以假设风速在相邻时段内满足v_{t1} μ ρ*(v_t - μ) ε_t其中ρ是时间相关系数通常0.7到0.95ε是服从正态分布的随机扰动。我之前做的一个工程里用AR(1)模型生成风速序列再卷积一个风机功率曲线出来的风电场景比单纯独立抽样合理多了。2.3 蒙特卡洛抽样的实现逻辑蒙特卡洛的核心就是“大量随机抽样”。我们可以通过MATLAB内置的rand和wblrnd威布尔分布随机数函数等函数来生成基础随机变量。以风电场景为例假设我们要生成N个场景每个场景包含T个时段给每个场景i初始化一个基准扰动用AR模型生成风速序列v(i, t)。将风速通过功率曲线转为风电出力P_w(i, t)。如果是风光联合场景再独立生成光伏出力序列P_pv(i, t)。把风电和光伏出力组合成总出力序列或者直接生成多维场景向量。这步生成的场景数不能太少太少会缺失尾部风险也不能太多否则后续削减压力大。一般原始场景数在500到2000之间比较合适具体看你的模型计算能力。3. 概率距离快速削减算法详解3.1 场景削减的本质是“概率测度的近似问题”我们原本有一个经验概率测度每个场景概率都是1/N。削减的目标是找到一个新的场景子集J规模为MM N并重新分配概率让新旧概率测度之间的某种“距离”最小。这里最常用的是Kantorovich距离它在两个离散分布的累积分布函数之间定义计算上又可以通过场景间欧氏距离的线性指派问题来实现。直觉理解的话可以把场景看成平面上的很多点每个点都有权重。现在要删掉一些点并把删掉的点的权重累加到离它最近保留点上让整体形状分布变化最小。这有点像在地图上去掉小城市把人口并入附近的中心城市尽量不改变全国人口分布。3.2 快速前向削减一个可用的具体步骤“概率距离快速削减算法”这个词与GAMS/SCENRED里经典的“fast forward selection”一脉相承。基本思路是从原始场景集合里逐个挑选“最重要”的场景进入保留集合直到数量达到M。步骤如下计算任意两个场景之间的欧氏距离矩阵D(i, j)这里距离可以按整个时间序列T维的欧氏距离也可以加权比如对高峰时段权重大一点。初始化保留集合J为空候选集合I包含所有N个场景。第一步选择一个使其他所有场景到它的距离之和最小的场景放入J。对剩余场景计算它们到J中所有场景的最小距离每次选择能使“所有候选场景到J的距离总和最小”的场景加入J。重复直到J的大小达到M。概率重分配对每个被删除的场景k找到J中距离它最近的保留场景j把概率加到j上最终保留场景概率pk 1/N sum_{k deleted assigned to j} (1/N)。这个算法的复杂度大约是O(M*N^2)当原始场景2000个、保留50个时也就是百万级距离计算MATLAB跑起来还行。如果要更快可以先用区块距离矩阵预计算。3.3 后向削减怎么选什么时候用前向还有一种相反思路后向削减。先从N个场景开始每轮删除一个对总概率距离影响最小的场景并把概率转移给最近邻居直到剩下M个。后向削减在小规模场景N 200时效果直观但每轮都要重新算N大了会很慢。快速前向削减更适合大规模场景生成之后的离线处理。我在实际使用中更推荐“先快速前向选择再用局部概率再分配”的组合因为前向选择由简到繁不容易出现后向删除时早删了不该删的场景导致不可逆的问题。4. MATLAB代码结构与核心实现4.1 主程序整体框架项目代码一般由几个模块组成参数设置、场景生成、场景削减、结果绘图。为了让后续扩展方便我习惯写成函数调用形式%% 主脚本 clear; clc; close all; % 1. 参数设置 N 1000; % 原始场景数 T 24; % 时段数小时 M 10; % 削减后场景数 wind_params [8.5, 2.0]; % Weibull(尺度c, 形状k) pv_params [2.0, 2.0]; % Beta分布alpha beta % 2. 生成原始场景 scenarios generate_wind_pv_scenarios(N, T, wind_params, pv_params); % 3. 削减场景 [reduced_scenarios, reduced_probs] fast_forward_reduction(scenarios, M); % 4. 绘图对比 plot_scenarios(scenarios, reduced_scenarios, reduced_probs);主程序里最需要调的是N、M和T。N多了削减耗时上升N少了尾部信息不够。M具体取多少要看后续优化模型的负担通常取5到30个就够了。4.2 场景生成函数实现这里以风电场景为例写一个简化版的风电场景生成函数function wind_scenarios generate_wind_scenarios(N, T, c, k, rho, P_rated) wind_scenarios zeros(N, T); v_in 3; v_out 25; v_rated 12; % 典型风机参数 for i 1:N v zeros(1, T); mu c * gamma(1 1/k); % 威布尔均值 v(1) wblrnd(c, k); for t 2:T eps 0.2 * wblrnd(c, k); % 简化噪声 v(t) mu rho * (v(t-1) - mu) eps; end % v - P P zeros(1, T); P(v v_in v v_rated) P_rated * (v(v v_in v v_rated) - v_in) / (v_rated - v_in); P(v v_rated v v_out) P_rated; wind_scenarios(i, :) P; end end函数里用gamma函数计算威布尔均值这步是让AR模型有一个稳定基准值。噪声项我故意加了0.2倍尺度否则序列太平滑会显得假。如果你要更精细可以改成用标准正态随机数乘上一个波动量。光伏场景类似可以用betarnd生成辐照度基础值再按时序平滑。我建议不要把风光分开跑而是写成一个联合函数这样后面可以加相关系数矩阵。4.3 场景削减函数实现下面是快速前向削减的核心函数距离用欧氏距离矩阵。注意MATLAB内存N2000、T24时距离矩阵约2万乘2万需要64GB内存不可行所以要用分块或者逐点计算。这里先展示一个针对中小规模N≤500的清晰版本function [reduced, probs] fast_forward_reduction(scenarios, M) [N, T] size(scenarios); % 计算距离矩阵 D zeros(N, N); for i 1:N D(i, :) sqrt(sum((scenarios - scenarios(i, :)).^2, 2)); end selected false(N, 1); J []; % 保留场景索引 % 第一步选择与其他场景总距离最小的场景 total_dist sum(D, 2); [~, idx] min(total_dist); selected(idx) true; J idx; % 逐步挑选 for ii 2:M remaining find(~selected); dist_to_selected D(remaining, J); min_dist_to_J min(dist_to_selected, [], 2); % 选择使“候选场景到J最小距离之和最小”的场景 [~, rel_idx] min(sum(min_dist_to_J)); new_idx remaining(rel_idx); selected(new_idx) true; J [J; new_idx]; end reduced scenarios(J, :); % 概率重分配 probs zeros(M, 1); for i 1:N if selected(i) local find(J i); probs(local) probs(local) 1/N; else dist_to_J D(i, J); [~, local] min(dist_to_J); probs(local) probs(local) 1/N; end end end这段代码里有个细节概率初始每个场景都是1/N保留场景本身自带1/N删除场景把它最近保留场景的概率加上1/N所以最终所有概率和为1。local用来定位保留场景在J里的位置。第一次挑选时没有“候选场景到J的最小距离”可以迭代所以单独处理。如果你需要处理N2000可以改成循环里逐行算距离dist_to_J sqrt(sum((scenarios(remaining(i), :) - scenarios(J, :)).^2, 2))意思一样但省内存。实际项目中我用过的数据N3000M20跑完耗时不到两分钟还是可接受的。5. 参数设置与结果评估5.1 关键参数推荐与调参逻辑场景生成和削减的参数有很多但最关键的几个必须理解清楚参数推荐范围选择依据原始场景数N500~2000N太小尾部丢失N太大距离矩阵内存爆炸削减后场景数M5~30取决于后续数学模型复杂度尽量取10~20时段数T24小时级如果做日前调度24点足够AR相关系数rho0.7~0.95时序平滑度太接近1会过平滑威布尔形状k1.5~3越大风速波动越小按地区季风特征调功率曲线参数风机手册简化时v_in3, v_rated12, v_out25调参有个笨办法先固定N1000M10然后反复改变时间相关系数rho看削减后的场景是否还保留原始曲线的波动规律。如果削减后曲线普遍比原始曲线更平滑说明大概率丢了边缘场景。5.2 削减效果的评价指标场景削减不是削完就完事了要会评价。常用评价指标包括削减前后总出力均值的相对误差mean(sum(reduced_scenarios, 2))与原始场景均值对比。概率分布距离比如用Wasserstein距离计算削减后与原始分布的差异这个可以直接用之前距离矩阵的加权和近似。峰值出力覆盖削减后场景的最大峰值是否接近原始场景的最大值。时序相关性保留度计算削减后各时段间相关系数与原场景对比。我一般会在主脚本里输出这些指标orig_mean mean(mean(scenarios, 2)); red_mean sum(sum(reduced_scenarios .* reduced_probs, 2), 1); fprintf(原始总出力均值%f\n, orig_mean); fprintf(削减后加权总出力均值%f\n, red_mean); fprintf(均值误差%f%%\n, abs(red_mean - orig_mean) / orig_mean * 100);如果均值误差超过5%我基本会怀疑削减算法或者M取太少。5.3 一个简化的运行示例我拿一个简化算例演示下效果原始场景N500T6小时的小规模为了方便看曲线削减到M3。威布尔参数c8.5、k2.2rho0.85。运行结果一般是这样的削减前的500条曲线密密麻麻削减后3条曲线呈现“高、中、低”三种出力特征概率大约是0.34、0.41、0.25。3条曲线分别对应原始场景里的大风段、中强风段和低风速段。这符合我们对风电不确定性的直觉少数几个典型场景就能代表大体分布。值得注意的是M3时均值误差可能在8%左右M10时通常能压到1%以内M20时进一步降低但改善幅度递减。这就是典型的拐点你需要结合自己模型的求解时间来选M。6. 常见问题与避坑指南6.1 削减后场景概率分布失真怎么办如果削减后的概率分布和原始分布差了太多首先检查距离矩阵是否用了标准化。原始风速和光伏出力尺度不同比如风电是0~1光伏是0~0.6如果不归一化距离矩阵会被数值大的变量主导。建议先把场景数据归一化到[0,1]算完再映射回去或者用加权距离。第二件要检查的是削减场景数M。当M太小时削掉一些概率不高的极端场景很正常但如果连均值都跑偏那就说明代表场景没选好。可以手动把原始场景中出力最大和最小的场景强制放进保留集合再跑一次前向选择这样可以避免尾部丢失。6.2 程序运行过慢怎么优化概率距离快速削减的大头是距离矩阵计算。N2000T24时矩阵约4000万元素MATLAB里双精度矩阵要320MB还能忍但N5000就会直接内存爆掉。优化方向不要一次性生成全量距离矩阵用“按需计算”的方式在循环里直接算到保留场景的距离省内存但慢。用并行计算MATLAB的parfor在距离计算上改善明显但要注意不能和parfor嵌套。先做个简单聚类把明显相同的曲线合并再跑削减。比如先把风速分成低、中、高三类各自内部再做前向选择。我试过把N3000、T24、M30的场景削减从两分钟降到30秒方法就是把距离循环改成矩阵分块乘加利用MATLAB的向量化指令快速求欧氏距离。6.3 风光联合场景怎么处理相关性风电和光伏在一个时段内往往有互补性比如夜里没光但风可能大白天风小但光强。如果独立生成风光场景再拼接会高估系统净负荷波动。更科学的做法是在生成随机数时加一个相关系数矩阵corr_matrix [1, 0.3; 0.3, 1]; % 风电-光伏相关性假设 R mvnrnd([0, 0], corr_matrix, T);然后用这些相关的正态随机数去对应风、光的随机变量比如通过逆变换。这块做起来相对复杂需要你对Copula理论有些了解但实际项目里效果显著。一个小技巧如果你不想用Copula也可以先独立生成风、光场景再按时间段做一个排序重配对让风大的时刻尽量配光小的时刻实现相关性的粗略拟合。6.4 MATLAB版本相关小问题不少朋友私信问过MATLAB版本问题或者安装问题这里多提一嘴这套代码用到的都是很基础的函数比如wblrnd、betarnd、gamma、parfor从R2016a到R2023b都能跑。如果你用的是较新版本注意wblrnd可能需要Statistics and Machine Learning Toolbox没装的话可以用逆变换手写u rand(1); v c * (-log(1 - u))^(1/k); % 威布尔逆变换这个方法完全不需要工具箱适合备无患。另外如果你的机器上MATLAB启动卡在“Setup没反应”可能是旧版本和系统兼容问题建议优先安装官方最新版旧脚本用兼容模式打开。最后再分享一点个人体会我最初做场景削减的时候疯狂调模型参数后来发现最有用的调试手段就是把削减前后的场景画在同一张图上用透明度区分一眼就能看出分布是否变形。削减算法说到底是为优化模型服务的不要一味追求削减误差小而忽略了后续模型的求解时间。你要是遇到M取30求解器还顶不住那不如把M压到10多试几组随机种子选一组在目标函数值上最稳定的结果。这套方法后续还可以继续扩展把场景削减和鲁棒优化结合或者用机器学习的方法生成更“真实”的场景。但只要蒙特卡洛加概率距离这个框架在风光不确定性刻画这个活儿你上手就不会慌。