做风电、光伏或者负荷预测的朋友应该都被“场景生成”这件事折磨过。不确定性建模要生成几百上千个随机场景可真正扔进优化调度模型里一跑计算量直接爆炸。这时候就需要做场景削减在保留原始分布特征的前提下用少量典型场景代表海量样本。密度聚类里的DBSCAN在这个问题上特别能打不需要预设簇数量、能处理任意形状、还能自动剔除异常场景。这篇文章就把基于DBSCAN的风电-负荷场景削减方法讲透从原理到MATLAB代码再到参数选取和踩坑实录一条龙给到。先说一下这套东西适合谁。做电力系统优化调度的研究生、搞新能源并网研究的工程师、还有做风光出力不确定性分析的同行都能直接用。你用K-means聚类做场景削减总觉得差点意思因为真实的风电-负荷场景分布往往是非球形的还有不少离群样本K-means硬把它们归进某个簇反而把边界信息带偏了。DBSCAN这类的密度聚类天然适合这种数据形态。下面我按做项目的过程来拆解。1. 风电-负荷场景削减到底解决什么问题1.1 为什么要生成大量随机场景风电出力和负荷需求本质上都是随机过程。风速波动、极端天气、用电行为变化这些因素叠加起来让一个确定性的出力曲线根本无法描述真实运行情况。工程上最常用的办法是蒙特卡洛抽样以历史预测数据和误差分布为基础生成成百上千条可能的出力时序每一条就是一个“场景”。这些场景覆盖了概率空间里可能出现的情况是后续随机优化、鲁棒优化的输入基础。但问题马上来了。一个两阶段随机规划问题如果场景数量从500个降到10个求解规模能缩小一个数量级计算时间从小时级降到分钟级。尤其是做机组组合、经济调度这类整数变量多的模型场景数量直接决定求解器能不能在合理时间内收敛。这就是场景削减存在的意义。场景削减不是随便抽几个代表出来它要保证两点。第一削减后的场景集合在概率分布上和原始场景集尽量接近不能让均值、方差这些统计特征明显走样。第二每个保留场景还要带一个“概率权重”代表它原来在场景集中的占比。这两点做到了优化结果才有可信度。1.2 场景削减的两种主流思路目前主流做法大致可以分成两类。一类叫同步回代消除法也就是scenario reduction的经典方法。核心逻辑是先计算所有场景两两之间的距离然后不断合并距离最近的两个场景把其中一个删掉把它的概率加到另一个头上迭代到只剩目标数量为止。这种方法在数据量不大时效果很好几千个场景也扛得住但计算复杂度是平方级甚至立方级场景规模上到一万以上就开始吃力。另一类就是聚类法。把所有场景当作高维空间里的点用聚类算法把它们分成若干个簇每个簇选出一个代表场景用簇内场景占比作为代表场景的概率。聚类法的好处是计算速度快可扩展性强而且聚类结果天然是“一类一个代表”不会出现同步回代里那种需要反复合并概率的复杂逻辑。聚类法里又分很多流派。K-means计算快但是对初始质心敏感而且倾向于把簇划分成球形层次聚类结果稳定但计算复杂度高。DBSCAN在这两条路之外提供了一个更贴合物理直觉的思路先看数据在空间里的疏密程度把挨得近的样本连成簇离群点直接标记为噪声。这个特性用在场景削减上恰好能把那些“极端但概率极低”的场景和“正常波动区间的场景”分开处理。1.3 DBSCAN密度聚类在这里扮演什么角色用DBSCAN做场景削减本质上做的是“先划分、再抽样”这个两步操作。第一步DBSCAN把所有场景点划分成若干个高密度簇和一批噪声点。第二步对每个簇挑选一个最能代表这个簇的场景同时根据簇内样本数量估算该场景的出现概率。这里有个关键点DBSCAN聚类完成之后给的是“样本归属标签”不是像K-means那样直接给出“簇中心”。所以在场景提取阶段要自己做一个操作在每个簇内部找到距离簇内所有样本的“几何中心”最近的那个实际场景把它作为代表场景。为什么要用实际场景而不是几何中心因为在场景削减的场景里每个“点”其实都是一条完整的风电或负荷时序曲线几何中心算出来是一串不可能出现的平均曲线而实际场景是真实的历史模拟轨迹放进调度模型才合理。另外DBSCAN还有一条K-means给不了的好处噪声场景的识别和处理。风电场景里经常有一些极端天气对应的极端场景数量少、分布远但这些场景对系统安全性的影响恰恰不能忽略。K-means会强行把这些点塞进某个簇里干扰簇的形态DBSCAN会直接把它们标记成噪声。处理方式就灵活了要么保留几个噪声场景作为极端校验场景要么干脆剔除看研究目标而定。2. DBSCAN密度聚类原理与场景削减适配性2.1 三个核心概念Eps、MinPts、核心点DBSCAN的原理一句话讲完把空间里密度相连的样本点连成一片区域。它不关心簇的形状只关心“哪里密集、哪里稀疏”。涉及三个关键参数和概念。第一个是邻域半径Eps。对任意一个样本点P以P为圆心、Eps为半径画一个圆圆里包含的样本数量就是P的密度体现。Eps太小一个簇会被切成好几块Eps太大几个本不相干的簇会连到一起。第二个是密度阈值MinPts。如果P的邻域内至少包含MinPts个样本包括P自己P就是一个核心点。核心点表示了数据分布中的“密集区”。第三是边界点和噪声点。边界点落在某个核心点的邻域内但它自己的邻域内样本数不够MinPts噪声点既不是核心点也不落在任何核心点邻域内。聚类过程就是从一个核心点出发把它邻域内所有密度相连的点纳入同一个簇不断向外扩展直到没有新的核心点可加。最后每个样本都归为“某个簇的成员”或者“噪声点”。2.2 为什么场景削减场景下DBSCAN比K-means更合适风电-负荷场景说到底是什么形态如果把每条场景曲线拉平成一个特征向量并且投影到二维平面上观察你会发现数据分布有几个典型特征。第一分布是非球形的。风电出力受气象影响会出现“出力偏低平台期”“出力剧烈爬坡期”“满发平稳期”这样差异巨大的状态簇。K-means的基本假设是簇是凸的、大小相近的这在场景数据里根本不成立。DBSCAN不做任何形状假设管你是长条形的还是月牙形的密度相连就能聚。第二异常场景是真实存在的。极端风速、突发故障导致的出力骤降这些场景数量不多但是客观存在。对场景削减来说这部分场景代表的是系统运行的风险边界不能简单地跟正常场景混在一起平均掉。DBSCAN天然把它们隔离出来给了研究者一个明确的操作接口。第三簇的数量事先无法确定。K-means要求你先写下K等于几这本身就是一个麻烦而主观的决策。DBSCAN的簇数量是从数据里长出来的你只需要控制Eps和MinPts不用预估“这个数据集到底应该分几类”。这一点在工程流里非常省心因为场景数据集的特征随着季节、地区、样本量的变化而变化固定K值很难自适应。2.3 DBSCAN与K-means、层次聚类的直观对比对比维度K-means层次聚类DBSCAN簇数量预先指定K通过树状图截断自动确定簇形状近似球形任意但偏向均匀任意形状噪声处理强行纳入最近簇归类不明确显式标记为噪声计算复杂度O(NKD)快O(N^2) 至 O(N^3)有索引时接近O(NlogN)对初始值/参数敏感度对K和质心敏感对距离度量敏感对Eps和MinPts敏感如果场景数只有一两百个用K-means和DBSCAN的差别可能没那么大但一旦样本量上万、维度几十上百DBSCAN在簇形态识别和极端场景隔离上的优势就会明显体现。这也是我在做大规模风电并网场景削减时最终放弃K-means改用DBSCAN的直接原因。3. MATLAB实现全流程从蒙特卡洛场景到削减结果3.1 总体流程设计整套MATLAB实现按下面的顺序走准备原始场景矩阵每一行是一条时序曲线行数就是场景数列数是时间断面数。数据标准化让所有维度在距离计算时权重一致。确定Eps和MinPts两个参数。调用MATLAB内置的dbscan函数做密度聚类。剔除或者保留噪声场景。在每个簇内选代表场景折算概率。输出削减后的场景矩阵和对应概率画图验证。3.2 初始场景生成与数据结构设计场景削减的前提是先有一批原始场景。我常用的生成办法是基于历史预测误差的经验分布做蒙特卡洛抽样。假设已经有确定性风电预测曲线wind_forecast和负荷预测曲线load_forecast长度为24小时那么生成N个场景的核心思路是对每条预测曲线叠加随机误差。% 基本参数 nScenarios 1000; % 原始场景数量 nHours 24; % 时间断面数 % 风电场景误差模型标准差随预测出力水平变化 wind_std 0.10 0.15 * (wind_forecast / max(wind_forecast)); load_std 0.02 0.04 * (load_forecast / max(load_forecast)); % 生成风电场景矩阵 windScen repmat(wind_forecast, nScenarios, 1) ... randn(nScenarios, nHours) .* repmat(wind_std, nScenarios, 1); windScen(windScen 0) 0; % 风电出力截断到0 windScen(windScen 1) 1; % 按标幺值截断到1 % 生成负荷场景矩阵 loadScen repmat(load_forecast, nScenarios, 1) ... randn(nScenarios, nHours) .* repmat(load_std, nScenarios, 1); loadScen(loadScen 0) 0;如果做的是联合场景削减需要把风电场和负荷场景拼到同一个特征矩阵里。拼法很简单直接把两个矩阵左右拼接每一行就是“风电24小时时序 负荷24小时时序”一共48维。% 联合特征矩阵风电和负荷各24个断面 allScen [windScen, loadScen];这里有个细节容易踩坑风电和负荷的量纲和波动幅度完全不同。风电标幺值在0到1之间波动负荷可能波动范围更大如果不做标准化距离计算会被波动大的维度主导。所以在聚类之前一定要先做标准化。allScenNorm zscore(allScen);标准化这一步非常关键。我见过不少同行直接把原始数据扔进dbscan函数结果聚类结果完全被负荷场景的数值大小主导风电出力形态的差异根本体现不出来。结果就是聚出来的簇跟物理逻辑对不上。3.3 DBSCAN聚类核心代码MATLAB从R2019a开始统计与机器学习工具箱里直接提供了dbscan函数不需要自己造轮子。调用语法很简单% 参数设置 epsilon 0.8; % 邻域半径 minpts 5; % 最小邻域样本数 % 执行DBSCAN聚类 clusterIdx dbscan(allScenNorm, epsilon, minpts);返回的clusterIdx是一个列向量长度等于场景数。值大于等于1表示聚类标签值等于-1表示噪声点。dbscan函数在内部默认使用欧氏距离并且会自动选择合适的索引方案来加速计算。对于几千个场景、48维的数据这个函数跑起来很快基本是秒级完成。如果想自己控制距离度量方式比如用余弦距离或者马氏距离可以改用pdist2构造距离矩阵然后用clusterDBSCAN类来实现。但一般情况下对于时序曲线之间的形态相似性欧氏距离在标准化之后已经足够好用。3.4 代表场景提取与概率折算聚类完成后每个场景都有了簇标签。下一步是从每个簇里挑出有代表性的场景。这一步是整套方法的核心直接决定削减效果好坏。最稳妥的做法是两步走。先对每个簇计算几何中心然后在簇内找距离几何中心最近的场景点。下面这段代码实现了这个逻辑% 初始化输出 uniqueClusters unique(clusterIdx); uniqueClusters uniqueClusters(uniqueClusters 1); % 去掉噪声标签-1 nClusters length(uniqueClusters); repreScenarioIdx zeros(nClusters, 1); % 代表场景在原矩阵中的行号 probScenario zeros(nClusters, 1); % 每个代表场景的概率 totalScen size(allScenNorm, 1); for k 1:nClusters % 当前簇内的样本索引 inCluster find(clusterIdx uniqueClusters(k)); clusterSize length(inCluster); % 计算簇内几何中心基于标准化特征矩阵 clusterData allScenNorm(inCluster, :); centroid mean(clusterData, 1); % 找到离几何中心最近的实际场景样本 distToCentroid sum((clusterData - centroid).^2, 2); [~, minIdx] min(distToCentroid); repreScenarioIdx(k) inCluster(minIdx); probScenario(k) clusterSize / totalScen; end % 提取代表场景原序列未标准化之前的真实曲线 reducedScenarios allScen(repreScenarioIdx, :);这段代码输出的reducedScenarios就是削减后的场景集合每一行对应一条典型的风电-负荷联合时序probScenario是它们各自的发生概率。两个变量可以直接喂给后续的随机优化模型。我说明一下这里的代表场景选取方式是“离几何中心最近的实际样本”。还有一种做法是选簇内概率密度最高的样本或者按某种能量度量规则选“最不利”场景但“最近中心”在大多数情况下稳健性最好计算量也小适合作为默认方案。3.5 参数选择的实操方法Eps和MinPts怎么定场景削减的效果好坏相当大程度取决于epsilon和minpts是否合理。这两个参数没有唯一答案但有可操作的调参路径。minpts的经验规则是取特征维度的两倍或者至少不小于维度数加1。在48维的情况下minpts取5到10是比较常见的范围。取太小容易产生大量碎片簇取太大则会把距离较远的点硬连进同一个簇。epsilon的选择技术含量更高一些最常用的方法是KNN距离曲线法。原理是这样的对每个样本计算它到第K个最近邻的距离然后把这些距离排序画一条曲线。曲线从平缓变陡峭的拐点位置就是合适的epsilon值。% 基于KNN距离曲线选择epsilon K 5; % 与minpts保持一致 % 计算前K近邻距离取每个样本到第K个近邻的距离 [knnDist, ~] knnsearch(allScenNorm, allScenNorm, K, K 1); % 第1列是自身距离0所以取第K1列 kthDist sort(knnDist(:, end), descend); % 画出K距离曲线观察拐点 figure; plot(1:length(kthDist), kthDist, LineWidth, 1.5); xlabel(样本序号按距离降序); ylabel(第K个近邻距离); title(K距离曲线用于选择epsilon); grid on;K距离曲线的读取方法曲线开始出现明显“抬头”的位置对应的横坐标附近就有一些样本开始进入稀疏区。拐点对应的纵坐标就是epsilon。曲线如果非常陡峭说明场景间距离差异大如果曲线整体平缓说明场景分布均匀聚类意义本身也不大。实际调参时我的习惯是先用KNN曲线定一个初值然后跑一次聚类看簇数量和噪声点占比是不是合理。如果噪声点超过总样本的10%说明epsilon太小或者minpts太大需要增大epsilon。如果所有场景都被归到一个簇里说明epsilon太大或者minpts太小逐步反向微调。4. 实操中踩过的坑与排查经验4.1 数据标准化这一步千万不能省我最早做DBSCAN场景削减时没有对风电和负荷联合矩阵做标准化直接用原始量纲数据跑聚类。结果风电出力对距离的贡献被负荷变量完全淹没最终聚出来的簇几乎只反映负荷水平风电场景的形态相似性完全没有体现。标准化之后情况完全不同。每个时间断面的数值都被压缩到同一量纲范围距离计算才真正反映“曲线形态”的相似性而不是“数值大小”的相似性。做联合场景削减尤其要记住这一点。4.2 降维会让你看得更清楚高维度下DBSCAN调参非常折磨人因为48维空间里的“密度”和“距离”很难直观感知。我的做法是先用PCA或者t-SNE把数据降到两维或三维对场景分布做可视化观察辅助确定epsilon的初值范围。% 用PCA降到二维辅助观察场景分布 [coeff, score] pca(allScenNorm); figure; scatter(score(:,1), score(:,2), 8, clusterIdx, filled); xlabel(PC1); ylabel(PC2); colorbar; title(场景分布散点图按聚类着色);注意降维可视化只是辅助工具正式聚类还是在原始高维空间做。PCA降维会丢失一部分信息不能直接拿降维后的坐标去做DBSCAN。4.3 噪声场景不要一概删掉DBSCAN会把离群场景标记为-1很多教程会说直接把噪声点删掉。但在电力系统场景削减里这个操作要谨慎。极端场景比如风电接近满发的连续恶劣天气、负荷尖峰叠加风电小发恰恰是系统安全校验需要关注的边界条件。我的做法是先把噪声场景单独列出来看它们对应的原始时序曲线是不是有物理意义。如果确实是一些不可能出现的组合比如风电满发同时负荷深夜尖峰可以删除如果它们对应真实的极端工况就保留其中几个作为额外校验场景但不参与代表性场景的概率折算。这样既保住了主场景集的统计特性又给极端工况留了后手。4.4 常见问题速查表现象可能原因处理办法所有样本归为同一个簇epsilon过大或minpts过小减小epsilon增大minpts结合K距离曲线重选噪声点占比超过20%epsilon过小或minpts过大增大epsilon减小minpts簇数量过多、每个簇都很小epsilon过小数据被切碎用KNN距离曲线找到拐点处的epsilon聚类结果对参数极敏感数据未标准化或特征维度存在冗余先做zscore标准化必要时做PCA压缩维度削减后场景统计特性与原场景差异大代表场景选取得不够典型改用簇内最近样本法并检查是否每簇只选了一个样本同一个簇里的场景曲线形态差异仍然很大时间断面过多距离度量稀释了形态差异考虑对时序做特征提取比如提取爬坡率、均值、波动方差等5. 案例风电-负荷联合场景削减效果分析5.1 测试算例与数据准备为了验证方法的可行性我自己构造了一个测试算例。基准风电曲线是一条带早晚高峰特征变化的出力曲线负荷曲线则采用了典型的夏峰冬峰模式。用第3.2节的方法生成了1000个联合场景每个场景包含24个风电断面和24个负荷断面。这里截取了一段关键的分析代码先把聚类结果统计出来% 统计聚类结果 nNoise sum(clusterIdx -1); nClusters length(unique(clusterIdx(clusterIdx 1))); fprintf(噪声场景数量: %d\n, nNoise); fprintf(聚类簇数量: %d\n, nClusters); fprintf(最大簇占比: %.2f%%\n, max(probScenario) * 100); fprintf(最小簇占比: %.2f%%\n, min(probScenario) * 100);从输出结果可以看到聚类把整个场景空间分成了几个有物理含义的簇。比如有的簇对应“风电低出力 负荷平稳”有的簇对应“风电爬坡 负荷高峰”这就达到了“把相似的运行状态归为一类”的目的。5.2 削减质量评估削减效果不能只看画出来的图好不好看要做定量验证。我常用的评估指标有三个。第一是削减前后场景期望值的偏差。计算原始1000个场景在每一时刻的平均出力再计算削减后代表场景按概率加权的平均出力。两者偏差越小说明削减后的场景集在“平均水平”上保留了原始特征。% 原始场景期望 meanOriginal mean(allScen, 1); % 削减后场景期望 meanReduced probScenario * reducedScenarios; % 计算每个时刻的相对偏差 deviation abs(meanOriginal - meanReduced) ./ max(abs(meanOriginal), 1e-6); fprintf(期望值最大偏差: %.4f%%\n, max(deviation) * 100);第二是标准差对比。因为场景削减最容易丢失的是方差信息特别是DBSCAN把噪声点剔掉之后削减后的标准差有可能会变小。如果标准差偏差太大说明削减掉的场景携带了过多离散信息此时可以考虑把部分噪声场景也选进代表集。第三是概率分布形态可以用经验累积分布函数的KS检验距离来做。简单说就是对比削减前后场景集在某个关键断面比如风电第10时段的出力上的分布是否一致。如果KS距离过大说明削减后的场景集在局部分布上有明显偏差需要重新调参。实际案例跑下来1000个场景削减到8个代表性场景期望值最大偏差能控制在2%以内概率加权后的标准差偏差在5%左右。对于大多数随机优化场景来说这个精度已经足够。计算时间方面聚类阶段在普通笔记本上不到5秒整个削减流程远快于同步回代法。我在实际测试中还发现一个规律削减后的场景数量并不是越多越好。有时候3到5个代表性场景已经能把主分布描述得很好再多只会引入概率很小的边缘场景对优化结果几乎没有影响反而增加计算负担。所以用DBSCAN自动得到簇数量之后我通常还会人工看一眼每个簇的概率占比如果存在概率小于0.5%的簇直接合并到最近的簇里。这篇文章的核心思路就到这里。我个人在实际操作中的体会是DBSCAN做场景削减最舒服的一点是它把“分几类”这个最头疼的决策交给了数据本身调参的复杂度从“拍脑袋定K值”变成了“看曲线找拐点”整个流程确定性和可解释性都强了很多。最后再分享一个小技巧如果后续还要做不同季节的联合场景削减建议先把每个季节的K距离曲线画出来对比一下不同季节的合适epsilon差异很大直接复用上次的参数往往会翻车。