做新能源规划、微电网调度或电力系统随机优化的同学，大概率都绕不开这两个动作：先生成一批新能源出力场景，再把它们削减到可计算的数量。简单说，“场景生成”是根据历史风速、光照数据模拟未来可能出现的一组组出力曲线，“场景削减”则是从这成千上百条曲线里挑出少数几条，还能代表原来的不确定性分布。Matlab在这两个环节里都很顺手——有统计工具箱、聚类函数、绘图又方便，所以我把自己实现新能源场景生成与削减的一整套流程和踩过的坑整理出来，希望能当成直接抄的作业，也帮新手省点弯路。

这篇文章不是入门科普，更像是一份带代码、带参数的实操记录。适合正在做风电/光伏出力不确定性建模、机组组合、微电网规划，或者刚接触场景法但被一堆概率公式劝退的同学。我会先讲清楚“为什么不能只用一条典型曲线”，再给出场景生成的具体建模思路，然后重点对比K-means聚类和同步回代削减两种打法，最后附上一段完整的Matlab实现和问题排查经验。

1. 为什么需要“先生成，再削减”？

1.1 一条确定性曲线解决不了随机优化

很多人一开始会犯一个直觉性错误：既然调度需要负荷预测、风光预测，那把每个时段的风速、光照取个平均值，不就能算了吗？

实际上，随机优化看重的不是“平均情况”，而是“极端组合”。比如某天风速整体偏低，但正好赶上午高峰负荷，光伏又碰上多云，三者叠加就是系统最紧张的时刻。这种组合在平均值曲线里根本看不到，但真实运行中出现的概率并不等于零。新能源的波动性、间歇性和相关性，决定了我们不能用一个“典型日”去代表一整年的可能性。

场景法要做的，就是把随机变量离散化为一个场景集合，每个场景是一条完整的时序曲线，附带一个出现概率。这一组场景在概率上逼近历史数据的分布，又保留了时间上的前后关联，可以直接放进优化模型里，让决策者考虑“如果明天是这种天气，我的机组怎么安排最稳”。

1.2 削减不是偷懒，而是让计算变得可行

既然场景越多越精确，为什么不生成几千条直接用？因为优化模型的计算量会随场景数量直线上升。假设你要做24时段机组组合，100条场景意味着决策变量在非场景模型基础上乘100，求解时间和内存占用根本不是一回事。到了微电网双层规划、鲁棒优化这些场景，500条甚至1000条场景直接会把求解器拖到难以容忍。

所以场景削减的意义在于：在损失极小信息量的前提下，把场景规模降下来。理想情况下，削减后的几十条场景，其整体概率分布、期望值、方差甚至分位数，都尽量贴近原始场景集。换句话说，削减不是“取个平均”，而是在场景集中挑出一组具有代表性的子集，并重新分配它们的概率权重。

理解这一点后，就知道削减算法的好坏不能只看“像不像”，还要看“分布失真大不大”。所以后文所有比较，我都会围绕概率分布保持度和实际使用效果去评价。

2. 场景生成：用概率模型造出一批“未来曲线”

2.1 新能源时序的统计特征

在动手写代码前，必须先想清楚要模拟的是什么。风电出力和光伏出力统计特征差异很大：

风电出力是非负的，有上限（额定功率），往往呈明显的右偏分布，低出力时段居多，而且具有较强的时序自相关性——今天风大，明天大概率也不会突然完全没风；光伏出力则相对规律，白天从零升到峰值再回落，夜晚恒为零，主要随机性来自云层遮挡，表现为“晴天基准曲线”上的高频波动。

这些特征意味着不能直接拿白噪声生成曲线，也不能简单套用纯高斯分布。较好的做法是抓住两个要素：边际分布（每个时段出力的概率形状）和时间相关性（上一时刻出力对下一时刻的影响）。

2.2 用AR(1)模型生成风电出力序列

最简单也最容易解释的方法，是一阶自回归模型AR(1)：

x_t = μ + φ·(x_{t-1} - μ) + ε_t

其中μ是历史平均出力，φ是自回归系数，ε_t是均值为零、标准差为σ的随机扰动项。

用Matlab实现时，核心就是按时间步长递推：

% 场景生成参数 T = 24; % 一天24个时段 nScen = 500; % 生成500条场景 mu = 0.35; % 历史平均出力（标幺值） phi = 0.85; % 自回归系数，体现时序连续性 sigma = 0.12; % 扰动标准差 rng(42); % 固定随机数种子，保证结果可复现 scen = zeros(nScen, T); for s = 1:nScen scen(s, 1) = mu + sigma * randn; for t = 2:T scen(s, t) = mu + phi * (scen(s, t-1) - mu) + sigma * randn; end end % 出力限制为[0,1]标幺值 scen(scen < 0) = 0; scen(scen > 1) = 1;

这里我把φ取0.85，意味着相邻时段出力相关系数很高，模拟出的曲线不会像噪声一样乱跳。如果φ取0.1，曲线会非常毛糙，明显不符合真实风速的惯性特征。σ则控制波动幅度，需要根据实际数据残差来估计。

不过这个版本有明显缺陷：强制截断会让分布失真，大量值被压到0或1。实际项目中我一般不直接用这种简单模型，而是把它当成测试代码，验证流程能否走通。

改进一点的做法是先产生高斯场景，再用正态累积分布函数转换到均匀分布，最后通过Beta分布逆变换得到指定边际分布：

u = normcdf(scen); % 高斯序列转为[0,1]均匀分布 aBeta = 2.5; % Beta分布形状参数，需要拟合 bBeta = 1.8; scenBeta = betainv(u, aBeta, bBeta);

这样可以保证边际分布符合Beta形状，又保留时间相关性，缺点是计算速度稍慢。风电数据充裕时，更推荐直接用核密度估计拟合历史出力的分布，然后做同样的概率映射。

2.3 光伏出力场景的简化生成方式

光伏出力的生成没法照搬风电模型，因为它的确定性部分太强了。比较实用的做法是分成“确定基底”和“随机波动”两层。

先用天文公式计算一年中每天每个时段的晴空理论出力，得到一个基准曲线；再通过历史数据分析天气折扣系数，比如晴天为0.9、多云为0.5、阴天为0.2。生成场景时，对每个场景随机抽取一个天气类型，然后用Beta分布或马尔可夫链模拟云层带来的逐时段波动。

% 假设clearSky是24维基准出力曲线，weatherScale是晴天折扣系数 T = 24; nScen = 300; rng(7); p = 0.6; % 晴空概率 cloud = rand(nScen, T) < p; % 简化：云层遮挡逻辑 scenPV = repmat(clearSky, nScen, 1) .* (0.9 - 0.3 * cloud) ... + 0.03 * randn(nScen, T); scenPV(scenPV < 0) = 0; scenPV(scenPV > 1) = 1;

实际做风光互补场景时，我习惯把风电场景和光伏场景分别生成，再按同一随机种子下的天气状态关联起来，避免风电大、光伏也大这类离谱组合。这个话题可以单独写一篇，这里点到为止。

3. 场景削减：从几百条到几十条的关键一步

场景生成只是第一步。真正决定优化模型规模和精度的，是场景削减这一步。削减方法很多，常用的两类是聚类法和同步回代法。

3.1 K-means聚类削减：直观、快速，适合初步尝试

K-means的思路很直接：把每条场景看成T维空间里的一个点，用聚类算法聚成K簇，每簇中心作为代表场景。场景概率取簇内成员数量除以总数。

在Matlab里，统计和机器学习工具箱自带kmeans函数：

% scen维度为 [nScen, T] K = 10; [idx, C] = kmeans(scen, K, 'Distance', 'sqeuclidean', 'Replicates', 10); % 计算每个簇的样本数量比例 clusterProb = histcounts(idx, K)' / nScen; % 以聚类中心作为削减后的代表场景 cutScene = C; cutProb = clusterProb;

注意kmeans默认输入是样本×特征，所以我们传场景矩阵scen，每个样本是24维向量。Replicates设成10是为了避免陷入局部最优。距离用欧氏距离平方，对大维度场景聚类比较自然。

这个方法的优点是快、简单，缺点也很明显：聚类中心是簇内平均值，会把尖峰削平，还会模糊时序形态。比如风电场景中那种“先小风后大风突变”的形态，平均后可能变成一条温和上升曲线，丢失真实运行中的爬坡信息。

3.2 同步回代削减：保留曲线形态，分布更准

同步回代法（Backward Reduction）是电力系统文献里最常用的场景削减方法，核心思路是迭代删除场景，每次删除一条让整体概率分布变化最小的曲线，并把被删场景的概率累加到离它最近的保留场景上。

算法流程是：

1. 初始给每条场景赋相等概率1/N；
2. 计算所有场景之间的概率距离D(i,j)；
3. 对每个待删场景i，找到离它最近的保留场景j，计算删除代价p_i·D(i,j)；
4. 找出代价最小的场景i删除，把概率p_i加到场景j上；
5. 重复2-4，直到保留场景数达到目标K。

我用Matlab写过一版简化实现：

function [redScen, redProb] = backwardReduction(scen, K) [N, ~] = size(scen); p = ones(N,1) / N; keep = true(N,1); D = squareform(pdist(scen)); % 欧氏距离矩阵 while sum(keep) > K cand = find(keep); bestCost = inf; bestRemove = -1; bestAdd = -1; for a = 1:length(cand) i = cand(a); for b = 1:length(cand) if b == a, continue; end j = cand(b); cost = p(i) * D(i,j); if cost < bestCost bestCost = cost; bestRemove = i; bestAdd = j; end end end p(bestAdd) = p(bestAdd) + p(bestRemove); p(bestRemove) = 0; keep(bestRemove) = false; end redScen = scen(keep, :); redProb = p(keep) / sum(p(keep)); % 归一化，防止浮点误差 end

这段代码的复杂度是O(N²·(N-K))，N比较大时会很慢。我一般把N控制在500以内，目标K控制在15左右，实际运行还能接受。如果生成2000条又想削减到5条，建议先用K-means粗聚到100条，再用同步回代精减，效率会高很多。

3.3 两种削减方式的对比选型

我的选型经验是：如果最终目的是做规划或者需要展示几个“典型日”给非技术同事看，K-means更友好；如果要用于机组组合、经济调度这类对不确定性敏感的场景，强烈推荐同步回代，或者用K-means预聚、同步回代精减的两步方案。

4. 完整实现：把生成、削减、评价串起来

讲完原理，我们串起一个完整的Matlab小工程。以一座100MW风电场为背景，用24时段风电出力数据做场景生成和削减。

4.1 参数设置与整体流程

% 主流程 clear; clc; close all; % 1. 场景生成 T = 24; nScen = 500; mu = 0.35; phi = 0.85; sigma = 0.12; rng(2024); scen = generateAR1Scenarios(nScen, T, mu, phi, sigma); % 2. 场景削减（两种方法都用一下，方便对比） K = 10; [scenKmeans, probKmeans] = kmeansReduction(scen, K); [scenBack, probBack] = backwardReduction(scen, K); % 3. 计算期望出力曲线并对比 avgOriginal = mean(scen, 1); avgKmeans = sum(probKmeans .* scenKmeans, 1); avgBack = sum(probBack .* scenBack, 1); % 4. 可视化 figure('Color','w'); subplot(2,1,1); plot(scen', 'Color', [0.7 0.7 0.7]); hold on; plot(avgOriginal, 'k-', 'LineWidth', 2); title('原始500场景与平均出力'); subplot(2,1,2); plot(scenBack', 'r-', 'LineWidth', 1.2); hold on; plot(avgBack, 'b-', 'LineWidth', 2); title('同步回代削减10场景与加权平均出力');

配合函数文件：

function scen = generateAR1Scenarios(nScen, T, mu, phi, sigma) scen = zeros(nScen, T); for s = 1:nScen scen(s,1) = mu + sigma * randn; for t = 2:T scen(s,t) = mu + phi * (scen(s,t-1) - mu) + sigma * randn; end end scen(scen < 0) = 0; scen(scen > 1) = 1; end function [center, prob] = kmeansReduction(scen, K) [idx, center] = kmeans(scen, K, 'Replicates', 10); prob = accumarray(idx, 1) / length(idx); end

4.2 削减质量评价指标

光看图不够，建议用数值指标评价削减好坏。常用的两个指标：期望曲线平均绝对误差和累计分布差异。

期望曲线平均绝对误差MAE：

maeKmeans = mean(abs(avgKmeans - avgOriginal)); maeBack = mean(abs(avgBack - avgOriginal)); fprintf('K-means MAE: %.4f\n', maeKmeans); fprintf('Backward MAE: %.4f\n', maeBack);

累计分布差异可以用一维KS统计量近似。做法是把每个时段的值全拉成一维向量，比较原始场景和削减后场景的经验累计分布：

[~, ksStat] = kstest2(scen(:), scenBack(:)); fprintf('KS统计量: %.4f\n', ksStat);

我自己的实测数据里，同步回代削减到10条时，MAE通常在0.01~0.03之间，KS统计量在0.05左右，完全能满足后续随机优化的精度需求。K-means的MAE可能差不多，但KS统计量经常更大，原因就是尾部场景被平均抹平了。

4.3 削减数量的经验选择

削减到多少条才合适？这不是拍脑袋决定的，要结合下游模型类型。

• 5条：用于概念验证、教学演示，速度极快，但可能丢极端场景；
• 10~15条：用于机组组合、经济调度，兼顾速度与精度，是大多数论文的常用值；
• 30~50条：用于需要精细刻画风险和价值损失的场合，比如储能容量配置；
• 100条以上：基本没必要，除非你跑分布式鲁棒优化，否则求解压力太大。

我的习惯是先用10条跑通模型，再逐步增加到15、20条，看目标函数值变化是否小于某个阈值（比如0.5%）。如果变化很小，说明当前削减数量够了；如果还在明显变化，就继续增加。

5. 实操中踩过的坑与规避方法

5.1 削减后概率必须重新归一化

同步回代迭代过程中，浮点累加可能让最终概率之和出现细微偏差。有些代码写完后概率和是0.9999或者1.0001，直接带入优化模型，约束条件会出现问题。所以我在函数最后总会加一行：

redProb = redProb / sum(redProb);

5.2 K-means中心不是真实场景，要有心理预期

聚类中心是平均值，出来的曲线平滑得很“干净”，峰值被压低，爬坡变缓。如果下游模型对爬坡约束敏感，直接使用均值中心会低估系统风险。解决方法是“最近样本法”：在K-means计算完簇后，从每个簇中挑一个离中心最近的真实场景作为代表。这样既保留了真实施工形态，又只需要在kmeans基础上加几行判断，强烈推荐。

5.3 随机数种子决定结果，发布代码要固定rng

我在写Matlab脚本时，凡是涉及随机数的部分都会显式设置rng。否则不同机器上跑出来的结果差异很大，复现论文数据时会被审稿人问死。设置rng(固定值)之后，生成的场景集和削减结果完全可复现，调试时也能定位问题。

5.4 同步回代法在小规模场景下也要注意内存

pdist函数会生成N×N矩阵，N=500时是25万个距离值，内存不到2MB，不痛不痒。但N=5000时矩阵有2500万个元素，约200MB，计算也开始吃力。这种情况下建议先降维：用PCA把24维降到5~8维，再跑削减，或者直接分块计算距离，不要一次全算。

5.5 工具箱缺失的替代方案

kmeans属于Statistics and Machine Learning Toolbox，如果没有这个工具箱，可以用自己写的简单聚类函数代替（Lloyd迭代10行左右）。同步回代法只用到了pdist和squareform，pdist在基础Matlab里也有，严格说不需要额外工具箱。所以完整流程里最依赖工具箱的还是可视化部分，plot、subplot这些谁都有。

我遇到过一种情况：换了新版本Matlab后，kmeans函数的默认选项变了，提示‘Distance’参数不再接受某个值。遇到这类报错，优先用help kmeans查看版本说明，或者干脆自己写简化版聚类，避免环境依赖。

6. 一点个人体会

做新能源场景生成与削减这件事，技术上不难，难的是理解每个参数背后的物理含义。AR模型的φ到底取多少，Beta分布的形状参数怎么拟合，削减后要不要保留爬坡信息——这些都不是数学公式能直接告诉你的，需要对着历史数据反复摸索。

我在实际项目中，会把削减后的场景和概率存成一个mat文件，给后续的优化模型当输入。存储格式也很简单：一个矩阵scenarios，一行一条曲线；一个向量probabilities，对应每条曲线的概率。这两个变量就是整个不确定性分析的核心资产。

一个小技巧：在跑完整套流程前，先花十分钟把原始场景画出来，观察时序的波动范围、相关性和极端值。很多时候，削减算法看起来效果不好，其实不是算法的问题，而是生成环节的分布就模拟歪了。先保证“生成得合理”，再讨论“削减得漂亮”，顺序不能反。