很多做新能源调度、储能配置、微电网规划的朋友，一开始接触“场景生成与削减”这个概念时，容易把它当成一个单纯的统计工具，觉得无非就是抽样、聚类、算距离。但真正上手用Matlab实现一遍后会发现，这套流程实际上是整个不确定性建模链条里最容易被低估的一环——场景生成决定了你输入给优化问题的“信息质量”，场景削减则直接决定了随机优化能不能在规定时间内算完、算得稳。我最近完整走了一遍风电出力场景的生成与削减流程，把蒙特卡洛采样、概率密度拟合、同步回代消除、K-means聚类这些方法在Matlab里逐个落地，踩了不少坑，也理清了每个环节背后的逻辑。这篇文章就是把这次实现之旅的完整过程、关键代码思路和实际经验整理出来，给正准备做这块的朋友一条可以直接参考的路径。
这篇文章适合三类人：一是做电力系统随机优化、鲁棒调度的研究生，需要准备输入场景集；二是做分布式能源、微电网容量规划或者储能配置的工程师，想把不确定性的影响量化到模型里；三是刚接触Matlab概率建模，想用一个完整案例把采样、拟合、聚类串起来的初学者。文中不会只贴代码，我会把每个方案“为什么这么做”也讲清楚，这才是照着写完代码之后真正能带走的东西。
1. 为什么需要“生成”和“削减”这两把刀
1.1 所谓“不确定性问题”到底长什么样
新能源出力最让人头疼的地方，不是它忽多忽少，而是你根本不知道明天下午三点它到底发多少。风电的随机性来自风速的波动和间歇，光伏则受云层遮挡、日出日落节奏影响。在做调度或者规划的时候，如果用确定性的“一条功率曲线”去算，结果大概率偏乐观——因为那条曲线只是众多可能性里的一种，可能还是均值的那种。
想要把不确定性合理地塞进优化模型，主流做法是“随机优化”或者“分布式鲁棒优化”。这类模型的输入不是一个确定值，而是一组可能出现的场景。这组场景就像对“明天风电出力的多种剧情推演”：有些情景风速平稳，有些情景午后有阵风，有些情景几乎无风。每个情景都带一个概率权重，整体构成对风电随机性的一个离散近似。
1.2 生成与削减在整个建模链条中的位置
理解这套流程，建议把三个环节分开看：场景生成、场景削减、场景应用。
场景生成的目的是还原不确定性本身的形态，通常的做法是拟合历史数据的概率分布，再按分布采样出大量场景。这一步讲究“够密够全”，你要能捕捉到尾部风险，不能只围绕均值附近打转。场景削减的目的则相反，是把上千个场景压缩到几十个有代表性的场景，同时尽量保留下原本的概率信息。这一步讲究“少而精”，因为优化模型里每多一个场景，就多一倍的决策变量和约束条件——如果是混合整数规划，场景数直接决定求解器会不会卡死。
生成是加法，削减是减法，加出来的东西越丰富，减完之后留下的精华才越有代表性。这两个环节是配合关系，不是选一个的问题。
2. 场景生成：从历史数据到概率出力曲线的Matlab落地
2.1 两类最常用的出力分布模型
先说风电。风电出力通常先从风速入手，风速的概率密度一般用两参数Weibull分布描述：
[ f(v)=\frac{k}{c}\left(\frac{v}{c}\right)^{k-1}\exp\left(-(v/c)^k\right) ]
其中 (k) 是形状参数，决定分布形态的“胖瘦”；(c) 是尺度参数，表征风速平均水平。在Matlab里拟合历史风速数据时，最省事的是用wblfit，但要对数组维度做检查；也可以用mle自定义概率密度函数做最大似然拟合，更灵活。
% 假设 wind_hist 是历史风速序列（m/s） [k_hat, c_hat] = wblfit(wind_hist);拿到 Weibull 参数之后，用wblrnd(k_hat, c_hat, Nscen, T)就能生成一个 (Nscen) 行、(T) 列的场景矩阵——每一行就是一个风速场景，每一列对应一个时段。
但风速不同，风机出力不是简单的线性关系。常见做法是用分段函数把风速转成功率：
% 风机功率曲线：vin切入风速, vr额定风速, vout切出风速, Pr额定功率 P = zeros(size(v)); P(v < vin | v > vout) = 0; P(v >= vin & v < vr) = Pr .* (v(v >= vin & v < vr) - vin) / (vr - vin); P(v >= vr & v <= vout) = Pr;这段代码逻辑很简单，但有一个隐蔽问题：直接用风速分布套功率转换，会改变最终功率的分布形态。也就是说，就算风速拟合得很准，功率输出场景的分布也不一定符合实际情况。如果手里有历史功率数据，更稳妥的做法是直接对功率数据做经验分布或者 Beta 分布拟合，别绕道风速。
2.2 蒙特卡洛采样与拉丁超立方的实际差异
很多初学者上来就用randn或者wblrnd直接蒙塔卡洛采样，代码跑了半天，回头一看场景集里全是“中不溜”的剧情，极端场景一个都没采到。这是纯随机采样的通病——样本在概率空间的覆盖不均匀。
拉丁超立方采样（Latin Hypercube Sampling, LHS）解决的正是这个问题。它的核心思想是把每个输入变量的取值范围等分成 (N) 个区间，在每个区间里各抽取一个样本，随机打乱后组合起来。这样采出来的 (N) 个点能均匀覆盖整个概率空间，尾部事件被采到的概率大幅提升。
Matlab内置了lhsdesign和lhsnorm，但lhsdesign默认生成的是 [0,1] 均匀拉丁超立方，需要用逆变换法映射到目标分布：
% 生成 [0,1] 的拉丁超立方设计 u = lhsdesign(Nscen, T, 'criterion', 'correlation', 'smooth', 'off'); % 通过 Weibull 逆CDF转换到风速域 wind_scen_lhs = wblinv(u, k_hat, c_hat);用'criterion', 'correlation'可以最小化各列之间的相关性，这个选项很有用，因为后续要做时段间的相关性控制时，这一步能减少不必要的伪相关。
从我的实测对比看，纯蒙特卡洛采1000个场景的覆盖效果，大致相当于LHS采300到500个——也就是说LHS能帮你在同样计算代价下，把分布信息描述得更饱满。如果你要做的是极端天气事件相关的可靠性评估，这个差异很关键。
2.3 相关性建模别漏掉
场景生成最容易翻车的地方，不是分布拟合，而是忽略时段之间的自相关性以及风电场之间的空间相关性。
风速是连续变化的物理量，上午十点的风速和十一点的风速通常不会突然从2 m/s跳到15 m/s。但独立采样出来的场景矩阵，每一列都不管前一列，相邻时段之间的关系完全是“断片”的。这样的场景投到调度模型里，会出现一些现实中根本不可能出现的出力跳变，优化结果自然会偏差。
常用的修正办法是对随机源施加相关性。以多元正态分布为例，先构造时段相关系数矩阵 (R)，再做 Cholesky 分解：
% R 是 T×T 的相关系数矩阵，可用指数衰减模型构造：R(i,j)=exp(-|i-j|/L) L = 4; % 相关长度，越大时段间连续性越强 for i = 1:T for j = 1:T R(i,j) = exp(-abs(i-j)/L); end end C = chol(R, 'lower'); % Z 是 Nscen×T 的标准正态样本 Z_corr = Z * C'; % 再将 Z_corr 映射回 Weibull 分布 wind_scen_corr = wblinv(normcdf(Z_corr), k_hat, c_hat);这里顺序是：先产生独立标准正态样本，乘 Cholesky 因子得到带相关性的正态样本，再通过正态CDF映射到 [0,1]，最后用目标分布的逆CDF取回风速。这套“正态桥”是工程里最常用的相关随机场构造方式，实现起来也简单。
实测下来，相关长度 (L) 取3到5个小时比较合理，太小场景还是碎，太大则会把所有时段绑得太死，场景之间区分度下降。
3. 场景削减：同步回代消除与聚类两条路线的横向对比
生成完一两千个场景之后，下一步就是削减。这里有两个主流路线：基于概率距离的同步回代消除（Fast Backward Reduction，简称SBR），以及基于聚类的场景归并。两个路线我在Matlab里都实现了，直观感受就是：SBR数学上更干净，聚类实现上更灵活。下面分开细说。
3.1 同步回代消除的贪心合并逻辑与实现
同步回代消除的核心思想非常直观：每一次迭代，找出“被删除后整体信息损失最小”的那个场景，把它删掉，同时把它原本的概率权重叠加到离它最近的场景头上，然后更新距离矩阵，继续找下一个，直到场景数降到目标值。
算法步骤：
- 初始化每个场景的权重 (p_i = 1/N)。
- 计算场景两两之间的距离矩阵（通常是欧氏距离，也可以加权）。
- 对于每个场景 (i)，找到它与其它场景的最小距离 (d_{i,min})。
- 找出 (d_{i,min}) 最小的场景（也就是“最不孤独”的场景）作为待删场景，把它删除。
- 把被删场景的权重加到距离它最近的那个场景上。
- 重复3-5，直到剩余场景数等于预设目标。
用Matlab写核心循环：
function [red_scen, red_w, idx_kept] = sbr_reduction(Scen, target_num) N = size(Scen, 1); weights = ones(N, 1) / N; idx_active = (1:N)'; % 记录仍存活的场景索引 while numel(idx_active) > target_num S = Scen(idx_active, :); w = weights(idx_active); D = pdist2(S, S); D(1:size(S,1)+1:end) = inf; % 清零对角线，避免自己算自己 % 每个场景到最近邻居的距离 [min_dist_self, min_idx] = min(D, [], 2); % 找到最近邻距离最小的场景 [~, victim_local] = min(min_dist_self); victim_global = idx_active(victim_local); % 把它概率权重加到最近邻居上 neighbor_local = min_idx(victim_local); neighbor_global = idx_active(neighbor_local); weights(neighbor_global) = weights(neighbor_global) + weights(victim_global); % 删除受害者 idx_active(victim_local) = []; end red_scen = Scen(idx_active, :); red_w = weights(idx_active) / sum(weights(idx_active)); % 重新归一化 idx_kept = idx_active; end这段代码逻辑上能跑通，但如果场景量大——比如N=2000——每轮循环都要重新算一次完整的 (pdist2)，复杂度是 (O(N^2))，2000个场景削减到50个，循环将近2000次，运行时间会让人崩溃。我在实测2000×24的矩阵时，纯for循环跑SBR跑了二十多分钟。
优化空间在两步：一是去掉一个场景后，只需要更新和它相关的部分距离，不必整体重算；二是考虑先做一次粗聚类（比如分成3倍目标数的类），在各簇内部做精确SBR，最后再跨簇合并。后一种方式速度提升非常明显，而且精度损失可以接受。
3.2 K-means聚类的概率化改造
K-means路线更简单：把每个场景当作一个 (T) 维空间里的点，然后直接聚类。每个簇中心就是典型的代表性场景，每个簇里原始样本的概率之和，就是该典型场景的权重。
% Scen: Nscen×T 的原始场景矩阵 % k: 目标削减后的场景数 [idx, C] = kmeans(Scen, k, 'Distance', 'sqeuclidean', 'Replicates', 10); red_scen = C; red_w = zeros(k, 1); for i = 1:k red_w(i) = sum(1/Nscen * (idx == i)); endK-means的关键问题不是怎么调用，而是两点：
第一，初始质心对结果影响很大。Matlab的kmeans默认用k-means++初始化，这个必须保留。Replicates参数建议设5到10次，每次独立初始化，保留误差最小的结果，能有效避免陷入局部最优。
第二，距离度量。如果你场景里有多个变量（风速、光照、负荷同时作为维度），量纲差异会让欧氏距离被数值大的那个变量主导。这时最好对每一列做Z-score标准化之后再聚类。具体做法是z = (Scen - mean(Scen)) ./ std(Scen)，聚类完成后，再把得到的簇中心还原回原始量纲。
3.3 两条路线怎么选
两个方法的差异其实很明显。
同步回代消除每一步都在做“最小概率损失合并”，理论上能最大化保留原概率分布的信息，而且在概率距离意义下有明确的收敛保证。缺点是计算复杂度高，而且规则固定，无法人为控制“每个典型场景代表的物理含义”。
K-means聚类的计算效率高得多，而且灵活性好，可以方便地加入混合变量、带权重的距离、甚至约束“每个场景簇最小规模”。代价是它没有那么强的概率理论基础——最小化的是几何距离，不一定等于最小化分布信息损失。
实际上，很多做配电网规划的团队是混着用的：先用K-means把2000个场景粗聚到200个，再用同步回代消除从200个恢复到30个。这个两级方案在“计算效率+概率保真度”上做到了较好的平衡，我自己实测的效果也很不错。
4. 削减结果怎么量化评价
削减完场景，不能只看图觉得“差不多”，一定要算指标。做研究要被审稿人追问，做工程也要给团队一个“凭什么信这套场景”的证据。
4.1 常用评价指标与实现
最常用的三个指标：
一是相对误差，比较削减前后场景集合的均值曲线和多条分位线偏差。公式是：
% 原始场景均值曲线 mean_ori = mean(Scen, 1); % 削减后场景加权均值 mean_red = sum(red_w .* red_scen, 1); err_rel = norm(mean_ori - mean_red) / norm(mean_ori);二是分位数偏移，分别求原始场景集和削减场景集在5%、50%、95%三个分位上的曲线，逐时段比较偏差。这个指标能看出削减方法有没有丢掉尾部风险。90%分位的极端出力场景在调度里往往决定备用容量的配置，如果削减后95%分位线明显变低，说明极端场景被削没了。
三是概率分布保持度，对场景集合按时间断面切分，看每个断面上削减后场景的概率质量分布和原始分布有没有显著差异。实际操作中可以用每个时段的均值、方差和偏度来对比。
4.2 一个完整算例的指标复盘
我拿某风电场的历史风速数据跑了一组测试：原始场景2000个，分别用SBR和K-means削减到30个，然后对比指标。
SBR的结果是相对误差约0.6%，95%分位数最大偏移约4%；K-means的相对误差约1.2%，95%分位数最大偏移约6%。SBR全面占优，这是符合预期的——它本来就是按概率距离设计的。但K-means的计算时间只有SBR的1/15，如果只是前期做趋势分析，这个精度完全够用。
值得提醒的是，如果场景用于随机优化，还必须做一步“稳定性验证”：把削减前后的场景输入同一个优化模型，比较决策变量的差异。有时候平均指标很漂亮，但个别时段的最优解完全不对，说明关键场景被削减掉了，这种情况要回头检查是不是标准化时把峰值削平了。
5. 参数调优与踩坑记录
5.1 场景数量选择
削减到什么程度才算合适，没有标准答案，但有一条经验可以参考：场景数至少要和决策变量的复杂度匹配。如果一个两阶段随机优化里第二阶段约束很多，20个场景可能连Scheffe场景都撑不不起来；反过来，如果只是做全年8760小时数据的分段聚合分析，50个场景和30个场景的差异并不明显。
我常用的操作是画一条“误差-场景数”曲线，分别算场景数取10、20、30、50、80、100时的相对误差，选在曲线“拐弯”的位置。一般来说，误差随场景数增加快速下降，但到了某一点后边际收益骤减，那个点就是性价比最高的场景数。这种做法比拍脑袋定个数字稳健得多。
5.2 标准化与权重归一化
聚类前做标准化，聚类后还原量纲，这个顺序很容易错。K-means聚出来的簇中心如果是在标准化空间里算的，直接套回原始空间做物理解释会偏差很大。正确做法是先存储标准化时的均值和标准差，还原时：
C_raw = C_std .* std_vals + mean_vals;权重归一化也一样。SBR每次合并权重之后，理论上概率总和始终为1，但浮点误差会让总和轻微漂移，所以最后统一red_w = red_w / sum(red_w)。K-means的权重是基于每个类内样本数占比，天然归一到1，但如果原始场景带权重（比如LHS采样后每个样本权重并不完全相同），那就要用加权计数。
5.3 相关矩阵非正定问题
Cholesky分解要求相关矩阵必须正定，但实际构造的相关系数矩阵经常会因为舍入误差或者经验数据的缺失值变成半正定甚至不正定。报错信息是Matrix must be positive definite，遇到别慌，首选处理办法是做一个特征值修正：把小于阈值的特征值拉到一个极小的正数，再重新组装矩阵。
[V, D] = eig(R); d = diag(D); d(d < 1e-10) = 1e-10; R_fixed = V * diag(d) * V'; C = chol(R_fixed, 'lower');这个方法对轻微的非正定问题很有效。但要是相关矩阵本身因为数据缺失有大片NaN，那就得先做数据插补，不能靠这个补丁兜底。
5.4 多风电场场景的空间相关性
如果项目里不只是一个风电场的出力，而是多个风电场同时建模，那还要考虑空间相关性。同一片风区里的几个风电场，出力往往存在明显的正相关性——大风天大家都多，静风天大家都少。独立采样会低估这种联动，导致优化模型认为“多场出力互补”的空间比实际更大。
处理方法和时段相关性类似，对多风电场构造一个块状相关矩阵，把场间相关系数和时段自相关系数都考虑进去，再用Cholesky分解控制随机源。这里有一个容易忽略的坑：场间相关系数矩阵必须保持正定，多个场两两相关系数填出来经常不正定，需要用上面的特征值修正做预处理。
做完这一整套流程之后，我最大的体会是：场景生成和削减这个活儿，表面上是个统计学问题，实际上是个工程权衡问题。纯随机采样容易丢尾巴，LHS能补这个短板；直接蒙塔卡洛好写，但跑优化算不动；削减太狠丢失风险信息，削减太少模型喘不过气。每一步都像是在“精度”和“计算代价”之间走钢丝，而Matlab提供的工具箱恰好给了你把各种方案都试一遍的便利条件。
最后说一个很多人会忽略的小细节：场景生成前，先花十分钟看看历史数据的时序图，确认有没有明显的季节趋势和日内周期性。如果存在明显的季节性差异，直接把全年数据混在一起拟合一个分布，出来的场景很容易在冬夏边界出现不伦不类的过渡形状。分季节或者分典型天气类型分别建场景，最后按天数加权合并，效果会好很多。这个是任何算法代码都替代不了的“先看数据”的功夫。
