我们在做电力系统随机优化或者微电网规划的时候，经常要面对一个很实在的问题：风光出力怎么描述？直接拿一整年时序数据丢进去，计算量吃不消；只取典型日，又怕丢掉极端情况。我自己最早是被一个“要生成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分布描述，参数可以通过历史数据的均值和方差去估计。
实际项目中我没法拿到每个风电场的实测风速分布参数，所以代码里可以做一个简化：允许用户直接输入风速的威布尔参数，或者用历史数据来拟合。如果完全没有历史数据，也可以用经验值，比如c=8.5m/s、k=2.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_{t+1} = μ + ρ*(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，规模为M（M < 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内存，N=2000、T=24时距离矩阵约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的最小距离”可以迭代，所以单独处理。
如果你需要处理N=2000，可以改成循环里逐行算距离：dist_to_J = sqrt(sum((scenarios(remaining(i), :) - scenarios(J, :)).^2, 2))，意思一样，但省内存。实际项目中我用过的数据N=3000，M=20，跑完耗时不到两分钟，还是可接受的。
5. 参数设置与结果评估
5.1 关键参数推荐与调参逻辑
场景生成和削减的参数有很多，但最关键的几个必须理解清楚：
| 参数 | 推荐范围 | 选择依据 |
|---|---|---|
| 原始场景数N | 500~2000 | N太小，尾部丢失；N太大，距离矩阵内存爆炸 |
| 削减后场景数M | 5~30 | 取决于后续数学模型复杂度，尽量取10~20 |
| 时段数T | 24（小时级） | 如果做日前调度，24点足够 |
| AR相关系数rho | 0.7~0.95 | 时序平滑度，太接近1会过平滑 |
| 威布尔形状k | 1.5~3 | 越大风速波动越小，按地区季风特征调 |
| 功率曲线参数 | 风机手册 | 简化时v_in=3, v_rated=12, v_out=25 |
调参有个笨办法：先固定N=1000，M=10，然后反复改变时间相关系数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 一个简化的运行示例
我拿一个简化算例演示下效果：原始场景N=500，T=6小时的小规模（为了方便看曲线），削减到M=3。威布尔参数c=8.5、k=2.2，rho=0.85。
运行结果一般是这样的：削减前的500条曲线密密麻麻，削减后3条曲线呈现“高、中、低”三种出力特征，概率大约是0.34、0.41、0.25。3条曲线分别对应原始场景里的大风段、中强风段和低风速段。这符合我们对风电不确定性的直觉：少数几个典型场景就能代表大体分布。
值得注意的是，M=3时均值误差可能在8%左右，M=10时通常能压到1%以内，M=20时进一步降低，但改善幅度递减。这就是典型的拐点，你需要结合自己模型的求解时间来选M。
6. 常见问题与避坑指南
6.1 削减后场景概率分布失真怎么办
如果削减后的概率分布和原始分布差了太多，首先检查距离矩阵是否用了标准化。原始风速和光伏出力尺度不同，比如风电是0~1，光伏是0~0.6，如果不归一化，距离矩阵会被数值大的变量主导。建议先把场景数据归一化到[0,1]，算完再映射回去，或者用加权距离。
第二件要检查的是削减场景数M。当M太小时，削掉一些概率不高的极端场景很正常，但如果连均值都跑偏，那就说明代表场景没选好。可以手动把原始场景中出力最大和最小的场景强制放进保留集合，再跑一次前向选择，这样可以避免尾部丢失。
6.2 程序运行过慢，怎么优化
概率距离快速削减的大头是距离矩阵计算。N=2000，T=24时，矩阵约4000万元素，MATLAB里双精度矩阵要320MB，还能忍；但N=5000就会直接内存爆掉。优化方向：
- 不要一次性生成全量距离矩阵，用“按需计算”的方式在循环里直接算到保留场景的距离，省内存但慢。
- 用并行计算，MATLAB的
parfor在距离计算上改善明显，但要注意不能和parfor嵌套。 - 先做个简单聚类，把明显相同的曲线合并，再跑削减。比如先把风速分成低、中、高三类，各自内部再做前向选择。
我试过把N=3000、T=24、M=30的场景削减从两分钟降到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，多试几组随机种子，选一组在目标函数值上最稳定的结果。
这套方法后续还可以继续扩展：把场景削减和鲁棒优化结合，或者用机器学习的方法生成更“真实”的场景。但只要蒙特卡洛加概率距离这个框架在，风光不确定性刻画这个活儿，你上手就不会慌。
