做新能源规划、微电网调度或电力系统随机优化的同学大概率都绕不开这两个动作先生成一批新能源出力场景再把它们削减到可计算的数量。简单说“场景生成”是根据历史风速、光照数据模拟未来可能出现的一组组出力曲线“场景削减”则是从这成千上百条曲线里挑出少数几条还能代表原来的不确定性分布。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/N计算所有场景之间的概率距离D(i,j)对每个待删场景i找到离它最近的保留场景j计算删除代价p_i·D(i,j)找出代价最小的场景i删除把概率p_i加到场景j上重复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更友好如果要用于机组组合、经济调度这类对不确定性敏感的场景强烈推荐同步回代或者用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); end4.2 削减质量评价指标光看图不够建议用数值指标评价削减好坏。常用的两个指标期望曲线平均绝对误差和累计分布差异。期望曲线平均绝对误差MAEmaeKmeans 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矩阵N500时是25万个距离值内存不到2MB不痛不痒。但N5000时矩阵有2500万个元素约200MB计算也开始吃力。这种情况下建议先降维用PCA把24维降到5~8维再跑削减或者直接分块计算距离不要一次全算。5.5 工具箱缺失的替代方案kmeans属于Statistics and Machine Learning Toolbox如果没有这个工具箱可以用自己写的简单聚类函数代替Lloyd迭代10行左右。同步回代法只用到了pdist和squareformpdist在基础Matlab里也有严格说不需要额外工具箱。所以完整流程里最依赖工具箱的还是可视化部分plot、subplot这些谁都有。我遇到过一种情况换了新版本Matlab后kmeans函数的默认选项变了提示‘Distance’参数不再接受某个值。遇到这类报错优先用help kmeans查看版本说明或者干脆自己写简化版聚类避免环境依赖。6. 一点个人体会做新能源场景生成与削减这件事技术上不难难的是理解每个参数背后的物理含义。AR模型的φ到底取多少Beta分布的形状参数怎么拟合削减后要不要保留爬坡信息——这些都不是数学公式能直接告诉你的需要对着历史数据反复摸索。我在实际项目中会把削减后的场景和概率存成一个mat文件给后续的优化模型当输入。存储格式也很简单一个矩阵scenarios一行一条曲线一个向量probabilities对应每条曲线的概率。这两个变量就是整个不确定性分析的核心资产。一个小技巧在跑完整套流程前先花十分钟把原始场景画出来观察时序的波动范围、相关性和极端值。很多时候削减算法看起来效果不好其实不是算法的问题而是生成环节的分布就模拟歪了。先保证“生成得合理”再讨论“削减得漂亮”顺序不能反。