做新能源电力系统研究的同行应该都遇到过这种尴尬局面想给调度模型加风光不确定性手头却只有一条预测曲线想跑随机优化场景数一多求解器直接卡死。我刚开始做这个方向时最头疼的不是建模而是怎么把“不确定性”变成能算的东西。后来接触了新能源场景生成与削减这套方法论再用Matlab一步步实现出来算是彻底理顺了。这篇博文就围绕“新能源场景生成与削减”这个主题记录我用Matlab从零实现的全过程内容包括场景生成的核心思路、削减算法的原理与代码、一个完整的风电出力案例以及我实打实踩过的坑。1. 项目概述到底要解决什么问题1.1 新能源出力不确定性与场景思想的由来风电、光伏这类新能源出力曲线天然带“脾气”。风速一变风机出力跟着跳云层一飘光伏曲线瞬时往下探。做电力系统规划、经济调度、备用容量优化时如果只用一条确定性预测曲线算出来的方案往往偏乐观实际一跑就崩。如果不确定性处理过头又会导致系统过于保守成本直线上升。于是大家想到一个折中办法把不确定性的可能变化轨迹“采样”成一组带概率的时空序列这就是场景。每个场景代表一种可能的新能源出力波形再给每个场景配上发生概率。最终优化模型不再面对单个点而是面对一堆可能的路口目标从“求一个解”变成“在概率意义下最优的解”。场景生成和削减就是这套流程的两个核心工序先生成足够多的场景去覆盖不确定性空间再通过削减算法用少量典型场景代替大量原始场景保证概率分布特征不明显失真。1.2 场景生成与削减的定义及适用场景场景生成通俗讲是“用随机模拟或历史统计模型生成N条可能的新能源出力曲线”。场景削减是从这N条曲线里挑出S条代表性曲线S远小于N并重新分配概率使削减前后场景集合的“概率距离”最小。这套思路几乎贯穿新能源相关的所有不确定性优化问题。典型应用包括含风电的机组组合与经济调度微电网日前优化运行分布式光伏配电网规划储能容量配置与充放电策略优化电力市场投标策略不管用在哪个场景底层逻辑都一样用有限几个典型场景去近似连续的概率分布让随机优化问题变成确定性优化问题的集合。Matlab在这件事上有天然优势矩阵运算快、自带概率分布函数和统计工具箱、绘图方便适合快速验证算法原型。2. 场景生成先把不确定性“摊开看”2.1 基于Monte Carlo的场景生成原理与流程场景生成最常用的基础方法就是Monte Carlo模拟。它的思路不复杂先根据历史数据或预测模型得到新能源出力的基准预测值再假设预测误差服从某个概率分布然后从这个分布里反复抽样每次抽样产生一个误差序列叠加到基准预测曲线上就形成一条场景。举个例子风速预测模型给出未来24小时每小时的风速预测值 v_pred(t)。实际风速 v(t) v_pred(t) ε(t)其中 ε(t) 是预测误差。工程上常用正态分布 N(0, σ²(t)) 描述 ε(t)σ(t) 通常跟预测时刻、预测提前量有关提前量越大σ越大。如果考虑时间相关性还需要在误差序列上引入自相关结构这里我用“协方差矩阵 多元正态分布”的方式来处理。整体流程分五步整理基准预测数据例如每小时的期望出力或预测风速。根据历史预测误差估计每个时刻的均值、标准差构造协方差矩阵。使用 mvnrnd 函数从多元正态分布中采样生成 M 组误差序列。将误差序列叠加到基准曲线上得到 M 条初始场景。对场景做合理性校验防止风速为负、出力越限等异常值。2.2 Matlab实现随机场景生成以风速预测误差为例我直接用一段可运行的Matlab代码演示假设未来24小时每小时风速预测值和标准差已经放进向量里。% 场景生成参数 M 1000; % 初始场景数量 T 24; % 时段数 % 基准风速预测24小时m/s v_pred [4.2; 5.1; 6.3; 7.8; 8.5; 9.1; 9.6; 8.8; 7.9; ... 6.5; 5.8; 5.2; 4.9; 5.5; 6.4; 7.1; 8.2; 8.7; ... 7.6; 6.8; 5.9; 5.0; 4.4; 3.8]; % 预测误差标准差每小时m/s一般后半夜置信度低误差稍大 sigma 0.6 0.3 * sin((0:23) * pi / 12); % 简单示例曲线 % 构造误差协方差矩阵考虑时间相关性采用指数衰减结构 C zeros(T, T); for i 1:T for j 1:T C(i,j) exp(-abs(i-j) / 4); % 时间间隔越小相关性越强 end end % 将标准差嵌入协方差矩阵 C diag(sigma) * C * diag(sigma); % 定义误差均值向量这里假设无偏全为0 mu zeros(T, 1); % 使用mvnrnd生成M组误差场景一行是一个场景 rng(2025); % 固定随机种子保证可复现 epsilon mvnrnd(mu, C, M); % 生成风速场景矩阵size M x T V_scenarios v_pred epsilon; % 简单校验风速不能为负 V_scenarios max(V_scenarios, 0);这里有个关键点为什么要构造协方差矩阵而不是每个时刻独立采样因为相邻时刻的风速误差有连续性如果独立采样前后时刻可能一个正误差一个负误差风速曲线会像锯齿一样剧烈抖动不符合真实物理过程。用指数衰减结构模拟误差相关性生成的场景就平滑多了。2.3 生成环节的实操要点与参数选择场景数量 M 怎么选经验上M 至少取 500~1000太小覆盖不住尾部风险太大后续削减耗时。有人一上来就生成几万条结果削减跑半小时都算不完。我个人的习惯是先取 1000如果发现削减后的结果不稳定再增量调整。正态分布的均值最好根据偏差修正。有些预测模型存在系统性偏差比如早上光伏预测经常偏大这时候误差均值就不是零要用历史误差的平均值修正 mu。否则生成的场景系统性地偏离实际。生成的场景里难免出现异常值比如风速变成负数、光伏出力超过装机容量。处理方式不是简单截断而是先看异常比例超过5%就说明参数可能有问题比如标准差设得太大或者分布类型不适合。截断操作要在所有场景都生成完之后统一做不然会破坏协方差结构。3. 场景削减去粗取精留下“代表”3.1 为什么必须削减计算成本与可解性假设生成了1000个场景每个场景24个时段。如果把这1000个场景直接塞进一个两阶段随机优化模型相当于要把所有约束复制1000份变量数量直接爆炸。我曾经试过用1000个场景做机组组合YALMIP Cplex 跑了两个小时都没收敛换成30个场景几十秒就出结果。场景削减的目的就是用尽可能少的场景保留原始场景集合的概率特征把计算规模压到工程可接受的范围。削减不是随便抽几个场景而是要保证削减前后两个场景集合之间的概率距离最小。这个“概率距离”在算法上通常用 Kantorovich 距离或 Wasserstein 距离度量我们不需要深究数学定义只要理解削减后剩下的代表场景必须在形态和概率分布上尽量接近原始场景集。3.2 常用削减算法对比聚类法 vs 快速前代消除法目前工程里用得最多的削减方法有两类基于聚类的方法和基于概率距离的场景消除法。聚类法比如k-means思路直观把1000个场景看成1000个点先聚成S类再用每个类的中心点作为代表场景类的占比作为概率。优点是好理解、实现简单缺点是对场景形态差异大的问题容易把极端场景平均掉而且K值的选取需要调参。快速前代消除法也叫Quick Forward Selection是另一条路。它不计算类的中心而是每次从剩余场景里选出一个“最不重要”的场景删掉把它的概率转移到离它最近的那个场景上。重复这个步骤直到剩余场景数达到目标值。这个方法的优点是保留原始场景的真实形态不会产生“合成场景”特别适合新能源出力这类需要保持尖峰、低谷和波动特征的场景。3.3 Matlab代码实现快速前代消除法三步走快速前代消除法的核心是维护一个场景集合和对应的概率向量。每轮迭代找到“被删除代价最小”的场景 j然后找除 j 外离 j 最近的场景 k把 p_j 加到 p_k 上删掉 j。被删除代价用概率乘以距离来定义这样概率小的极端场景如果形态也很独特也不容易被误删。我给出一个可以直接运行的函数实现。function [scen_red, prob_red] quick_forward_selection(scen, prob, S_red) % 快速前代消除法 % 输入: scen - Nxd 矩阵每行一条场景d为时段数 % prob - Nx1 概率向量和为1 % S_red - 目标场景数量 % 输出: scen_red - S_redxd 削减后场景 % prob_red - S_redx1 削减后概率 N size(scen, 1); idx 1:N; % 当前剩余场景的原始索引 p prob(:); % 当前概率向量 while length(idx) S_red % 计算当前剩余场景两两之间的距离 scen_cur scen(idx, :); dist pdist2(scen_cur, scen_cur, squaredeuclidean); % 将自身距离置为无穷避免选中自己 dist(logical(eye(length(idx)))) inf; % 对每个场景计算删除它的最小转移代价 min_dist min(dist, [], 2); % 每个场景到最近场景的距离 cost p .* min_dist; % 删除代价 概率 * 距离 % 找到代价最小的场景 [~, pos] min(cost); del_idx_local pos; % 在场景集合中的位置 % 找到离它最近的场景 [~, nearest_local] min(dist(del_idx_local, :)); % 把删除场景的概率加到最近场景上 p(nearest_local) p(nearest_local) p(del_idx_local); % 删除该场景及其概率 idx(del_idx_local) []; p(del_idx_local) []; end scen_red scen(idx, :); prob_red p; end调用方式很简单S_red 30; [ scen_red, prob_red ] quick_forward_selection(V_scenarios, ones(M,1)/M, S_red);这段代码默认所有场景等概率生成所以初始概率是 1/M。削减后 prob_red 会自动归一化到和为1因为转移概率的过程中总的概率质量保持不变。有一点要注意pdist2 计算的是两两距离当场景数很多时需要 O(N²) 内存。我一般都是先生成1000个场景再削减Matlab计算没问题如果场景数超过5000就要改用分块计算或者先把场景聚成若干类再做精细削减。3.4 削减效果评价指标削减完之后不能拍脑袋就说“效果不错”。我常用三个指标做量化评估场景概率距离削减前后场景集合的Kantorovich距离越小说明分布保持得越好。期望值偏差比较削减前后场景集合的均值曲线即 E(scen_red) 与 E(scen_orig) 的绝对误差。比如风速均值曲线偏差不超过0.2m/s说明削减后的典型场景还能代表平均水平。分位数误差比较原始场景集和削减后场景集的某一置信度分位数比如10%和90%分位数。新能源优化通常关心极端场景分位数误差太大会导致备用留得不够。我自己在项目里一般这样验证如果削减后30个场景的均值曲线和原始1000个场景的均值曲线几乎重合峰谷时段误差都在0.5%以内那么这个削减质量就算过关。4. 完整项目实操案例风电场出力场景生成与削减全流程4.1 案例背景与数据准备这个案例我直接用实际项目里用过的流程来复现。假设有一个100MW风电场要生成未来24小时的出力场景用于次日的储能调度策略优化。原始数据是风电场的历史预测风速和实际风速我提前处理成了24小时预测风速基准和预测误差统计量。为了让场景更贴近真实这里不直接用风速场景而是先把风速转换为风机出力。由于风机出力与风速之间存在分段函数关系我使用一个简化的功率曲线切入风速3m/s额定风速12m/s切出风速25m/s。风速低于切入风速或高于切出风速时出力为0切入风速到额定风速之间按三次方关系近似。完整流程如下生成风速场景复用第2节代码M1000。将每个风速场景转换为出力场景单位MW。用快速前代消除法将1000个出力场景削减为30个代表场景。统计削减后场景的概率分布并画出图形。4.2 完整代码流程与关键注释% 参数设置 M 1000; % 初始场景数 S_red 30; % 削减后场景数 T 24; % 时段数 cap 100; % 风电场装机容量 MW % 这里直接沿用第2节生成的 V_scenarios (M x 24) % 如果没有先运行第2节代码生成 V_scenarios % 定义风机功率曲线分段函数向量化处理 v_in 3; v_rated 12; v_out 25; P_scenarios zeros(M, T); V V_scenarios; % 转换出力按照简化功率曲线 idx_zero (V v_in) | (V v_out); idx_part (V v_in) (V v_rated); idx_rated (V v_rated) (V v_out); P_scenarios(idx_part) cap * ((V(idx_part) - v_in) / (v_rated - v_in)).^3; P_scenarios(idx_rated) cap; % 削减场景 prob_init ones(M, 1) / M; [ P_red, prob_red ] quick_forward_selection(P_scenarios, prob_init, S_red); % 输出削减后场景概率 disp(削减后场景概率前10个为例); disp(prob_red(1:10));转换过程要注意向量化尽量不要用循环遍历所有场景否则1000×24的循环在Matlab里会慢得让人绝望。上面这种逻辑索引写法运行速度至少快一个数量级。4.3 结果可视化与工程意义解读可视化是判断场景质量最直观的手段。我会画出四张图第一张是原始风速场景和削减后风速场景的对比第二张是每个削减后场景的出力曲线透明度按概率大小设置第三张是削减前后均值曲线的对比第四张是某个典型时段比如第14小时的概率分布直方图对比。figure; hold on; % 画削减后的30个场景概率越大线越粗 for i 1:S_red lineWidth 0.5 3 * prob_red(i); % 按概率调整线宽 plot(1:T, P_red(i,:), b-, LineWidth, lineWidth); end xlabel(时段 (h)); ylabel(出力 (MW)); title(削减后30个风电出力场景); grid on;从工程角度看削减后的场景能帮你回答几个关键问题明天哪个时段的出力范围最大哪些时段最容易出现高出力或高出力储能应该在哪个时段充电、哪个时段放电才能兼顾概率最高的场景这些判断原本需要盯着一千条曲线现在看二三十条就够了。我实际跑下来的结果是削减到30个场景后均值曲线与原始1000场景均值偏差最大的时段在第21时偏差约0.6MW相对误差不到1%。分位数方面10%分位数最大偏差约2.1MW90%分位数最大偏差约1.8MW完全满足工程精度要求。但计算速度差异极大用1000场景的模型内存占用翻几十倍求解时间从50秒暴涨到无法接受用30场景后模型快速收敛结果也有说服力。5. 常见问题与排查技巧实录5.1 随机数种子与结果可复现很多同行找我第一句话就是“为什么我跑出来的场景和你帖子里不一样”。大概率是因为随机数种子不同。代码里的 rng(2025) 就是干这个用的。Matlab里只要设置同样的种子mvnrnd 生成的随机序列是一样的。但要注意Matlab在不同版本之间随机数生成器的算法可能会调整尽量在代码开头用 rng(default) rng(seed) 的方式显式设置同时把版本信息记录在工程文档里。5.2 协方差矩阵非正定怎么破这是我踩过最大的坑。用历史误差数据估计协方差矩阵时如果采样时段少、变量多估计出来的矩阵经常半正定甚至不正定。mvnrnd 在遇到不正定矩阵时会直接报错。解决办法是给协方差矩阵做特征值修正把负的特征值置为很小的正数再重构矩阵或者干脆用第2节里的指数衰减结构去构造协方差矩阵天然正定省事很多。另一个技巧是使用 mvnrnd 前先检查最小特征值if min(eig(C)) 0 [V_eig, D_eig] eig(C); D_eig(D_eig 0) 1e-6; C V_eig * D_eig * V_eig; end5.3 削减后场景概率归一化问题快速前代消除法理论上保证概率和始终为1但如果你的初始概率本身没归一化或者代码操作中有浮点误差累积最终概率可能变成0.999999或1.000001。模型里如果对概率有严格的等式约束最好在削减函数末尾加一行 prob_red prob_red / sum(prob_red); 强制归一化。这里还要留个心眼削减后的概率能不能为0如果某个场景被删得只剩下一个概率自然是1。但如果场景削减目标数量大于等于活跃场景数某些概率为0的场景可能还留在集合里这时候要顺手剔除valid prob_red 1e-10; prob_red prob_red(valid); scen_red scen_red(valid, :);5.4 性能优化与加速技巧场景削减中的距离计算是最耗时的环节。1000场景的两两距离计算在普通笔记本上大概耗时不到2秒还算能忍。但如果到了5000场景pdist2 会占用约200MB内存削减循环也会慢很多。我总结了三招用平方欧氏距离代替欧氏距离省去开方操作因为距离排序不受单调变换影响。在迭代中不要每次都重新计算所有距离可以只更新与被删场景相关的距离列但这会让代码复杂度上去适合追求极致的场景。如果场景数特别多先用k-means粗聚成100类再用快速前代消除法把100类削减到30类速度和精度都兼顾。在Matlab里还可以尝试把 pdist2 换成自定义距离函数但通常内置函数已经做过优化自己写循环反而更慢。另一个实用技巧是开启并行池如果要做多次不同随机种子的场景削减可以用 parfor 并行跑能省不少时间。5.5 手动调整场景数时的注意事项场景数量并不是越大越好也不是越小越省事。我的经验是场景数从10翻到30计算时间涨得不算夸张但精度提升明显从30到50精度提升就变慢了从50到100求解器开始变得吃力而精度收益微乎其微。如果你最终要接的优化模型是非线性或混合整数规划建议场景数控制在20~40。如果只是求期望值比如算期望成本20个场景可能足够。如果是做风险规避或CVaR分析需要保留尾部风险场景场景数至少要40而且削减时不能只看均值还要对比分位数。具体怎么选建议用一组历史日数据做回测画出“场景数-目标函数值”曲线找到拐点。结尾一点实操后的个人体会这套流程我前前后后跑了差不多三周才稳定下来最深的感受是场景生成看分布假设场景削减看距离度量两者要配合调整。光盯着算法代码本身很容易陷入“调参→看图→再调参”的循环。我自己最终确定了一套固定流程先用1000个场景做一次削减对比10%、50%、90%三组分位数如果偏差都在容忍范围就不再增加场景数。另外在项目交付时一定要把随机种子、协方差矩阵表达式和削减算法版本写清楚这三点直接决定别人能不能复现出你的结果。如果后续你打算把场景接入日内滚动优化还可以试试在削减前给每个场景加一个时间权重——离当前时刻越近的时段权重越大这样削减出来的场景对前几个小时更敏感调度结果会更实用。这个方向我还在试等有了稳定的结论再跟大家细聊。