去年做园区光储容量配置项目时，我第一次只用平均光伏出力曲线去做优化，结果业主一句话把我问住了：“连续阴雨天怎么办？”我盯着那条光滑的平均曲线算出来的储能容量，确实心里没底。后来才踏实补上这一课：新能源优化里，随机性不是噪声，是必须正面处理的约束。解决办法就是先把不确定性“摊开”——用场景生成产出一大批可能的出力轨迹，再用场景削减收拢到少量典型场景，最后拿这套场景集去跑优化。这套流程用Matlab完整走一遍并不复杂，但每一步都藏着不少坑。这篇文章就记录我用Matlab做新能源场景生成与削减的完整过程，从方法选型、代码实现到削减后的场景校验，适合正在做微电网、光储优化、电力系统随机规划或者刚接触随机优化的同行参考。
1. 随机性这个坑：为什么单条曲线做决策不靠谱
1.1 只拿均值曲线做优化，结果会在极端天气翻车
光伏和风电出力的本质是强随机过程。云层遮挡、辐照度波动、风速随机性这些因素叠加在一起，让出力曲线每一时刻都在一个区间内漂移。很多人做容量配置或调度优化时，习惯直接把历史数据取平均得到一条“典型日曲线”，然后用确定性优化去算。问题是均值曲线丢了方差信息，等于默认预测完全准确。一旦真实天气偏离均值——比如连续三天阴雨、大风天出力骤降——之前算好的配置方案就可能直接失负荷。
我做第一个项目时就是拿全年8760小时数据平均成了一条“平均日出力曲线”，储能容量按这条曲线来定。结果业主追问极端场景时，我根本没数据支撑。后来我意识到，方法是错的：随机优化需要的不再是一条曲线，而是一组能覆盖不同天气模式的曲线族。这个曲线族在学术和工程里就叫“场景集”，每个场景是一条完整的时间序列，代表一种可能发生的新能源出力轨迹。
1.2 场景集：用一堆曲线代替一条曲线，用典型场景代替全部时序
场景生成与削减这对组合拳，基本思路其实很朴素。第一步，通过采样或模拟生成足够多的候选场景，比如500个或1000个，目标是让这些场景能覆盖历史数据体现出来的主要统计特征，包括均值、方差、时序相关性，还有尾部的极端情况。第二步，用场景削减算法把这大批场景压缩到10到20个有代表性的典型场景，并给每个场景分配一个概率。优化模型里不用再面对500个约束组，只用处理10个到20个约束组，计算规模直接降一两个数量级。
这里要讲清楚一个概念：场景削减不是随便抽几个点，而是要在“保留概率分布信息”和“控制计算规模”之间做权衡。削减后的场景集，均值、标准差和分位数应该尽量接近原始场景集，极端场景也不能被平均掉。否则你辛辛苦苦做随机优化，最后算出来的方案跟确定性优化没啥本质区别，还多了一堆计算量。在我的项目实践里，这套流程已经成了固定动作：生成、削减、校验、进优化，四步一步都不能省。
2. 场景生成的三条技术路线与Matlab选型
2.1 路线一：参数分布拟合法（以Beta拟合光伏出力为例）
最简单直接的场景生成方式，是假设新能源出力服从某个参数分布，然后从分布里采样。光伏出力归一化之后的值域是[0,1]，天然适合Beta分布，因为Beta分布恰好定义在[0,1]区间内。风电风速通常用Weibull分布拟合，再通过功率曲线转换成出力。
Matlab里做Beta拟合非常方便，核心函数是betafit和betarnd。我拿历史光伏出力数据来演示：
% 读取历史光伏归一化出力，矩阵维度为 nDays x 24 % 这里提取每天中午12点的出力（第12列），用于拟合Beta分布 noon = P_hist(:, 12); noon = noon(noon > 0.01); % 去掉接近0的样本，避免Beta拟合退化 phat = betafit(noon); nScen = 200; noonScen = betarnd(phat(1), phat(2), nScen, 1);这段代码跑起来很快，几毫秒就能生成几百个随机场景。但Beta分布拟合法有个明显的短板：它只能刻画单个时刻的分布，没法天然表达时序相关性。中午12点的出力分布和下午13点的出力分布是有强相关性的，单独对每个时刻采样会得到大量“时间上乱跳”的伪场景，这种场景扔到储能调度模型里，会让充放电策略变得极其激进，因为系统误以为出力可以随意跳变。
2.2 路线二：历史自举加扰动（工程上我最常用）
如果不想被参数分布束缚，我推荐用历史自举法。思路也很直接：从历史数据里随机抽取一整天的出力曲线作为基底，再叠加一个小幅随机扰动。这样做的好处是不假设任何参数分布，历史数据的真实相关结构能保留下来。如果担心单日块的随机性不够，可以改成块自举，每次随机抽一段连续多天的历史块，这样跨天的持久天气过程也能保留。
我封装了一个简单的块自举场景生成函数，你们可以直接拿来用：
function scen = genScenBlockBootstrap(histData, nScen, blockLen, noiseStd) % histData: nDays x T 的历史出力矩阵 % blockLen: 块长度，保留时间相关性 % noiseStd: 扰动噪声标准差 [nDays, T] = size(histData); nBlocks = ceil(T / blockLen); scen = zeros(nScen, T); for k = 1:nScen picks = randi(nDays, nBlocks, 1); raw = zeros(nBlocks * blockLen, 1); for b = 1:nBlocks raw((b-1)*blockLen+1 : b*blockLen) = histData(picks(b), :); end raw = raw(1:T)'; % 对齐到T个时点 noise = noiseStd * randn(1, T); scen(k, :) = raw + noise; scen(k, scen(k, :) < 0) = 0; % 截断到[0,1] scen(k, scen(k, :) > 1) = 1; end end调用时假设P_hist是nDays x 24的历史光伏归一化出力矩阵，我通常会生成500个场景作为初始候选集：
rng(2025); % 固定随机种子，保证结果可复现 nScen = 500; scen0 = genScenBlockBootstrap(P_hist, nScen, 24, 0.03); prob0 = ones(nScen, 1) / nScen; % 初始等概率噪声标准差noiseStd的取值要看数据本身的波动水平。如果光伏出力本身波动很大，0.03可能偏小，生成出来的场景会过于集中在历史轨迹附近；如果数据本身比较平滑，0.05甚至0.08也行。建议先跑几步试算，看生成场景的标准差和历史数据的标准差是否接近。这一步没有固定答案，跟具体数据集关系很大，多调几次自然有手感。
2.3 路线三：马尔可夫链和随机微分方程（适合做研究）
如果要做更精细的时序建模，可以考虑马尔可夫链或随机微分方程（SDE）。马尔可夫链的思路是把出力值离散化成若干个状态，统计状态之间的转移概率矩阵，然后按转移概率一步步模拟。SDE则是用连续的随机过程来描述出力动态，常见的有均值回归过程加噪声。
这两种方法在Matlab里都有办法实现，但调参成本明显上升。马尔可夫链需要确定状态数和转移概率估计方法，SDE需要估计漂移项和扩散项参数，搞不好还要用最大似然估计去拟合。我的看法是：如果你的目标是写论文，这些方法能提供更漂亮的理论框架；如果目标是工程交付，自举法通常已经够了。我做过对比，用马尔可夫链生成场景再接削减，和用自举法生成场景再接削减，最终优化出来的储能配置差别不大，但自举法的代码量只有马尔可夫链的三分之一。
2.4 选型对比与我的建议
| 生成方法 | 核心原理 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|---|
| 参数分布拟合（Beta） | 用Beta分布拟合归一化出力样本 | 实现简单、速度极快 | 难以刻画时序相关性和多峰分布 | 快速估算、单点分布分析 |
| 历史自举加扰动 | 重采样历史块并叠加噪声 | 不假设分布、保留时序结构 | 对历史数据质量依赖高 | 工程落地首选 |
| 马尔可夫链/SDE | 对状态转移概率或随机过程建模 | 能显式建模时序转移 | 参数调优复杂、代码量大 | 学术研究、精细化模拟 |
有一点要提醒：生成场景的规模不是越大越好。500个场景做削减，Matlab几分钟能跑完；5000个场景做削减，距离矩阵就可能把内存挤爆。我见过有人一上来生成几万个场景，结果削减代码直接OOM。合理做法是先生成500到1000个候选场景，削减到10到20个，既保留多样性又控制计算量。如果确实需要超大场景集，建议先用K-means粗削减到500，再用后面的同步回代消除做精细削减，两级削减法在计算效率和场景质量之间平衡得很好。
3. 同步回代消除与K-means削减：原理、代码和选择
3.1 同步回代消除的原理：每次删除代价最小的场景
场景削减做得最多的算法是同步回代消除，英文缩写SBR。它的核心逻辑很直白：每次从现有场景集里找出一个“删掉后整体概率分布损失最小”的场景，把它删除，同时把这个场景的概率转移到离它最近的保留场景上。重复这个过程直到场景数降到目标值。
这里的关键在于“损失”怎么定义。常用定义是概率距离，具体计算方式是：对每个场景，找到它到最近保留场景的欧氏距离，然后乘以该场景的概率作为删除代价。概率大、离其他场景远、代表性强的场景，删除代价高，会留到最后；概率小、离其他场景近、可有可无的场景，会先被删掉。这套机制保证了削减后的场景集在概率分布上最接近原始场景集。
3.2 SBR的Matlab实现：注意距离矩阵内存优化
SBR实现不复杂，但有几个细节容易写错。我直接给出完整函数：
function [scenOut, probOut] = scenReductionSBR(scen, prob, nKeep) % scen: N x T 场景矩阵 % prob: N x 1 场景概率，初始可全为 1/N % nKeep: 目标保留场景数 N = size(scen, 1); idx = (1:N)'; % 预先计算场景两两欧氏距离矩阵 % 场景多时这块最吃内存，建议分块计算或只存下三角 D = pdist2(scen, scen); while length(idx) > nKeep Dsub = D(idx, idx); Dsub(1:length(idx)+1:end) = Inf; % 对角线置为Inf，排除自身 [dmin, jmin] = min(Dsub, [], 2); % 每个场景到最近保留场景的距离 cost = prob(idx) .* dmin; % 删除场景带来的概率距离损失 [~, rm] = min(cost); % 选择损失最小的场景删除 r = idx(rm); j = idx(jmin(rm)); % 被删场景的最近保留场景 prob(j) = prob(j) + prob(r); % 概率转移 idx(rm) = []; end scenOut = scen(idx, :); probOut = prob(idx); end两个容易出问题的地方。第一，对角线置Inf不能省，否则每个场景到自己的距离是0，删除代价永远会选中自己，算法就废了。第二，概率转移那一步，被删场景的概率要加到离它最近的保留场景上，而不是丢弃。如果只删概率不加，削减后场景概率总和会小于1，后续优化模型的目标函数加权就会整体偏小，算出来的结果没有经济意义。用我上面这个函数，概率和会始终保持为1。
另一个工程细节是距离矩阵的内存问题。pdist2会生成N x N的稠密矩阵，N=500时没什么压力，N=5000时就是5000乘5000乘8字节，约200MB，还能接受；N=10000时直接飙到800MB，很多机器就跑不动了。实际项目里我很少直接用上万个场景跑SBR，真要遇到大场景集，先用K-means粗削减到500再进SBR，内存完全不是问题。
3.3 K-means削减的快速实现与致命短板
K-means也能做场景削减，而且核心代码更短：
rng(7); [idxC, C] = kmeans(scen0, nKeep, 'Replicates', 30, 'MaxIter', 1000); probC = accumarray(idxC, 1) / size(scen0, 1); % C: nKeep x T 的簇中心，即削减后的典型场景 % probC: 每个簇的样本占比，即场景概率注意这里有个概念上的关键区别：K-means输出的是簇中心，也就是一条“虚拟场景”，它不一定存在于原始场景集里，甚至可能不满足物理约束。比如光伏出力本来应该在0到1之间，但簇中心在某些时刻可能取到一个中间值，这个值在实际天气条件下可能不合理。SBR输出的则是原始场景的子集，物理真实性天然满足。
还有一点，K-means按欧氏距离聚类，倾向于把形状相近的曲线归为一类，对极端场景不敏感。一个概率只有0.01但出力极低的极端场景，在聚类时很容易被平均进某个大簇里，簇中心被“稀释”掉。而SBR的删除代价考虑了概率和距离，极端场景只要概率不小，就有机会存活。所以我的建议是：K-means适合做粗筛，SBR适合做精细保留。最终进优化模型的场景，尽量用SBR或至少用SBR校验一遍。
3.4 两种算法怎么配合用
| 对比项 | SBR同步回代消除 | K-means聚类削减 |
|---|---|---|
| 输出类型 | 原始场景子集（物理真实） | 簇中心（虚拟场景） |
| 概率更新 | 被删场景概率转移给近邻 | 按簇样本占比 |
| 极端场景保留能力 | 较强，概率和距离共同决定 | 较弱，容易被均值淹没 |
| 计算复杂度 | 需要NxN距离矩阵，大样本内存高 | 依赖聚类迭代，相对快 |
| 适用阶段 | 精细削减，最终确定典型场景 | 大规模场景粗筛 |
我在实际项目里的固定组合是：先生成1000个候选场景，用K-means粗削减到100个，再用SBR精削减到10到20个。这样内存压力小，最终场景质量也高。如果你只跑一次优化，场景数需求不大，直接用SBR从500个削减到10个也完全没问题，代码更少，省得绕一圈。
4. 完整Matlab流程：从历史数据到可投入优化的场景集
4.1 数据准备与归一化
场景生成的前提是有一份干净的历史出力数据。我的做法是先把原始功率数据除以装机容量，归一化到[0,1]区间，然后剔除异常值。数据里偶尔会出现负值，通常来自传感器零漂或数据错误，直接截断到0。超过1的值大概率是数据重复或口径不一致，也截断掉。另外要注意数据的时间分辨率和缺失值。我用过一个站点数据，里面有几个小时的记录是NaN，如果不处理直接做自举，生成出来的场景会带着NaN一路传到优化模型里，求解器直接报错。
load('solar_hist.mat'); % P_hist: nDays x 24，已归一化 P_hist(P_hist < 0) = 0; P_hist(P_hist > 1) = 1; P_hist(isnan(P_hist)) = 0; % 简单处理缺失值，实际情况建议插值4.2 场景生成与削减的完整主流程
数据准备好之后，依次调用前面写的函数。下面是主流程脚本：
%% 1. 参数设置 rng(2025); nScen = 500; % 初始候选场景数 nKeep = 10; % 最终保留场景数 blockLen = 24; % 块自举块长度，1天 noiseStd = 0.03; % 扰动标准差 %% 2. 场景生成 scen0 = genScenBlockBootstrap(P_hist, nScen, blockLen, noiseStd); prob0 = ones(nScen, 1) / nScen; %% 3. 场景削减（SBR） [scenRed, probRed] = scenReductionSBR(scen0, prob0, nKeep); fprintf('削减后场景数: %d, 概率和: %.4f\n', length(probRed), sum(probRed)); %% 4. 可视化 figure; plot(scen0', 'Color', [0.8 0.8 0.8]); hold on; plot(scenRed', 'LineWidth', 1.5); xlabel('时刻/h'); ylabel('归一化出力'); legend('削减前', '削减后'); title('场景生成与削减结果');这段脚本跑完，你会看到一张图：背景是500条灰色细线，前景是10条彩色粗线。从视觉上就能直观判断削减效果——10条典型场景应该覆盖住灰色线条的主要包络范围，既不是全都挤在中间，也没有明显脱离群体的离谱轨迹。我在项目里每次跑完都会先看这张图，比任何数值指标都直观。
4.3 场景集的保存与结构设计
削减后的场景集我习惯存成Matlab的.mat文件，方便后续直接加载到优化脚本里：
save('scenario_result.mat', 'scenRed', 'probRed', 'nKeep', 'scen0', 'prob0');保存时我会把削减前的scen0和prob0也一起存下来，原因是后续校验时需要对比削减前后的统计差异。如果只存削减后结果，校验时还得重新生成一遍初始场景，纯属浪费时间。
4.4 用parfor提升生成速度
自举法生成场景的循环本质上互相独立，可以并行。如果nScen比较大，或者历史数据矩阵很大，把普通for改成parfor能显著提升速度，前提是装了Parallel Computing Toolbox。方法是在函数里把外层循环改成parfor k = 1:nScen。我第一次用这个优化时，从1000个历史曲线生成10000个场景，耗时从十几秒降到了三四秒，在场景数需求大的场景里非常划算。
5. 削减后场景的工程校验与高频踩坑记录
5.1 三个校验指标：均值、分位数、时序形态
削减不是跑完就完事，必须校验削减后的场景集是否保留了原始场景集的关键统计特征。我一直用的校验方式有三个。
第一个是均值曲线对比。计算削减前所有场景的逐时均值，再计算削减后场景的加权均值，看两条曲线是否接近。加权均值要把每个场景的概率乘进去算，不能简单平均。
mean0 = mean(scen0, 1); meanR = sum(scenRed .* probRed, 1); errMean = max(abs(mean0 - meanR)); fprintf('均值最大偏差: %.4f\n', errMean);第二个是分位数对比。这个指标特别重要，它能告诉你尾部场景有没有被削掉。做法是取几个关键分位点——我常用5%、50%、95%——分别对比削减前后的逐时分位数。如果95%分位数接近但5%分位数显著偏高，说明低出力场景被削掉了，储能配置会偏小，真实场景下容易失负荷。
qtl = [0.05 0.5 0.95]; T = size(scen0, 2); q0 = zeros(T, 3); qR = zeros(T, 3); for t = 1:T q0(t, :) = quantile(scen0(:, t), qtl); qR(t, :) = quantile(scenRed(:, t), qtl); end errQ = max(abs(q0 - qR), [], 'all'); fprintf('分位数最大偏差: %.4f\n', errQ);第三个是时序形态检查。把削减后的场景曲线画出来，人工看一眼有没有剧烈跳变。如果某条场景在相邻时刻从0.8直接掉到0.1再拉回0.7，这大概率是噪声参数设太大，或者自举时块拼接处把不同日期的数据硬接在了一起。这种场景对储能调度模型的影响很大，因为系统会误以为出力波动极快，进而改变充放电策略。
5.2 场景数怎么定：不是越多越好
场景数太少，极端场景覆盖不足；场景数太多，优化模型求解时间指数级上升。我的经验是10到20个场景对大多数容量配置和调度问题都够用。如果目标函数变化随场景数增加开始趋于平缓，就说明再增加场景数对优化结果影响不大了，可以把场景数定在那个拐点上。我习惯做一组敏感性测试：分别用5、10、15、20个场景跑优化，看目标函数值变化。如果从10加到15目标函数只动了不到1%，那10个就够；如果还差3%以上，就继续加。
5.3 我踩过的坑和对应处理
| 坑 | 现象 | 处理办法 |
|---|---|---|
| Beta拟合遇上大量0值 | betafit报错或拟合参数退化 | 剔除极小值样本，或加eps做平滑 |
| 大场景集内存爆炸 | N=10000时pdist2直接OOM | 分块算距离，或先用K-means粗削减到500再SBR |
| 随机种子不固定 | 每次生成削减结果不一致 | 统一用rng(2025)固定种子 |
| 削减后概率和不为1 | 目标函数加权整体偏小 | SBR会自动概率转移，手动实现时别漏这步 |
| 时间分辨率对不上 | 15分钟数据被当成1小时处理 | 生成时保持原始时间分辨率，削减后不要重采样 |
| 只看均值不看分位 | 尾部极端场景被磨平 | 加5%和95%分位数校验，必要时强制保留极端场景 |
这里面我要单独强调一下随机种子。Matlab里如果不固定种子，每次运行randi、randn都会得到不同结果，导致生成的场景集每次都不一样，优化结果自然对不上。工程上要保证可复现，第一行代码固定rng(2025)，这个习惯我从一开始就坚持到现在，省了很多沟通成本。
5.4 场景集接入随机优化模型的两种姿势
削减后的场景集最终要进优化模型。我常用Yalmip工具包做建模，核心思路是把目标函数按场景概率加权，每个场景的约束各自独立：
% Yalmip随机优化模型骨架：目标为场景概率加权 x = sdpvar(nVar, 1); % 决策变量 obj = 0; cons = []; for k = 1:nKeep % 对每个场景分别构造约束，使用场景 scenRed(k, :) cons = [cons, constraint_k]; % 例如功率平衡、储能动态约束 obj = obj + probRed(k) * cost_k; end optimize(cons, obj);如果不熟悉Yalmip，也可以用Matlab优化工具箱自带的optimproblem结构，逻辑类似。核心思想只有一个：概率高的场景对目标函数贡献大，概率低的场景同样在约束里起作用，不能被忽略。这样算出来的方案才是在全场景期望意义下最优的，而不是只在平均曲线上最优。
6. 从项目实践里总结出的几条实在经验
最后说说我做完几个项目之后沉淀下来的判断标准。第一次做场景削减，我总觉得保留得越多越安全，结果20个场景的优化问题求解时间比10个场景多了将近一倍，但目标函数只改善不到1%。从那以后我就养成了先做敏感性测试的习惯，不要凭感觉定场景数。
另一个体会是校验环节坚决不能省。我早期图省事，削减完直接扔进优化器，结果算出来的储能配置在5%分位数场景下失负荷，后来加上分位数校验才发现场景尾部被磨平了。现在我的流程固定是“生成、削减、校验、进优化”，如果均值偏差和分位数偏差都控制在一个很小的范围内，我才会把场景集交给下一步。如果校验不通过，先回头调噪声标准差或场景数，而不是直接修改优化模型去硬凑结果。
如果你刚开始接触新能源场景生成与削减，我强烈建议先用自己手头的数据把上面这套流程完整跑通，感受一下不同噪声参数和场景数下削减结果的变化规律。跑通之后再考虑替换成更复杂的Copula或深度学习生成器也不迟，基础路线的工程收益已经能覆盖大部分实际需求。
