很多做新能源调度、储能配置、微电网规划的朋友一开始接触“场景生成与削减”这个概念时容易把它当成一个单纯的统计工具觉得无非就是抽样、聚类、算距离。但真正上手用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这段代码逻辑上能跑通但如果场景量大——比如N2000——每轮循环都要重新算一次完整的 (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提供的工具箱恰好给了你把各种方案都试一遍的便利条件。最后说一个很多人会忽略的小细节场景生成前先花十分钟看看历史数据的时序图确认有没有明显的季节趋势和日内周期性。如果存在明显的季节性差异直接把全年数据混在一起拟合一个分布出来的场景很容易在冬夏边界出现不伦不类的过渡形状。分季节或者分典型天气类型分别建场景最后按天数加权合并效果会好很多。这个是任何算法代码都替代不了的“先看数据”的功夫。