DBSCAN场景削减MATLAB实现：风电-负荷随机优化高效聚类方法-平芜编程栈

做新能源电力系统随机优化的人，大概都经历过这种痛苦：一上不确定性，风机出力、负荷曲线就变成一堆场景，几百上千条，调度模型转头就跑不动了。场景削减要干的活，就是在这堆场景里挑出几个最有代表性的把概率分布特性尽量留下来，同时把计算规模压下去。我最近在MATLAB里把基于DBSCAN密度聚类的风电-负荷场景削减方法完整跑了一遍，从场景生成、参数标定到聚类削减和概率分配，踩了不少坑，也把整套流程理顺了。这篇文章分享给正在做随机规划、新能源并网或鲁棒调度的研究生和工程师，手里有场景不知道怎么削，或者觉得K-means削得不够准的，都可以参考。

1. 风电-负荷场景削减：为什么我选了DBSCAN

做这个方向之前，我自己也用过不少场景削减手段，真正上手之后才发现，方法本身不复杂，难的是搞清楚每种方法的适用边界。下面把问题本质和方法选型的过程捋清楚。

1.1 场景削减到底在解决什么问题

风电出力和负荷都是随机变量，数学上可以用概率分布描述，但在优化模型中，处理连续的随机分布很不方便，所以工程上普遍采用场景法：先通过历史数据或统计模型抽样生成大量离散场景，每个场景是一条24小时时序曲线，并附带发生概率。把期望值、风险指标写成这些场景下的组合，模型就变成一个可计算的确定性问题。

场景越多，模型对不确定性的刻画越精细，但计算量也成倍上涨。以两阶段随机规划为例，场景数从100涨到1000，求解时间往往不是涨10倍，而是涨几十倍。场景削减就是在“精度”和“计算量”之间找平衡——用少量典型场景逼近原始场景集的概率分布。

经典做法有同步回代消除法、快速前向选择法，以及聚类方法。同步回代消除每次删除对概率分布影响最小的场景，并把被删场景的概率累加到距离最近的场景上，精度不错，但处理上千个场景时计算代价不低；快速前向选择反过来逐步挑选新场景，速度快但对初始选择比较敏感。聚类法把场景当样本聚类，每簇提取代表场景，实现简单、可解释性强，也容易嵌入随机规划流程，这就是DBSCAN路线的基本出发点。

1.2 K-means的三个毛病，DBSCAN正好绕开

很多人第一反应是用K-means做场景削减，因为代码短、工具箱现成。但K-means有三个让电力场景很难受的毛病。一是K值要预先给定，实际工程里你很难拍脑袋说风电场景该分成8类还是12类。二是K-means偏向凸形簇，风电场景高发、低发、爬坡、反调峰，形态并不规整，硬切成凸簇会丢掉过渡形态的场景。三是对离群点敏感，个别极端天气场景会把簇中心拉偏，让代表场景失真。

DBSCAN靠密度连通来定义簇，不需要预设簇数，能识别任意形状的簇，还能把稀疏区域的场景自动标记为噪声。这对场景削减其实是一个隐藏优势——极端场景往往落在这堆场景的边缘位置，与其强行把它们并进某个簇里产生一条不伦不类的代表场景，不如单独挑出来，作为小概率极端场景单独考虑，这才是保留风险特性的正确姿势。

当然DBSCAN也有脾气：高维时序场景直接丢进去，距离会变得很稀疏，参数也不好调。所以后面的完整链路是：先降维、再聚类、再回到原始空间取代表场景。这套流程我实测下来，稳定性比K-means高不少。

2. DBSCAN核心原理与参数确定：调好eps和minPts

DBSCAN的原理写进论文就三四句话，但真正调参时坑不少。我尽量用大白话把核心概念和参数选择讲透，并给出能直接在MATLAB里跑的K-距离曲线脚本。

2.1 核心点、边界点、噪声点：一个生活化理解

DBSCAN把每个场景点分成三类。核心点：以它为中心、半径eps的邻域内至少有minPts个点（含自身）；边界点：落在某个核心点的邻域内，但自己邻域内不够minPts个点；噪声点：既不满足核心点条件，也不挨着任何核心点。

聚类过程是这样的：随便找一个未访问的核心点出发，把它邻域内的所有点拉进来，再对这些点的邻域继续扩张，所有密度相连的点就形成一个簇。边界点被并进它所属核心点的簇里，但不再往外扩张。噪声点不参与任何簇。

这里有个生活化的类比：想象把人按“聚在一起的程度”分组。核心点是小区里人缘好、邻居多的住户；边界点是住在核心点隔壁、自己周围却没什么邻居的人；噪声点就是方圆几里只有他一家。DBSCAN把所有挨着核心点的住户连成一片社区，“孤零零”的房子就排除在外。这个直觉比公式好记，也直接决定了参数调整的方向。

2.2 eps和minPts到底怎么取值

两个参数直接决定聚类结果。eps太小，大部分点都变成噪声；eps太大，所有点融成一团，聚类等于没做。

minPts的常见经验取法有两个：一是取维度数加1（dim+1），这是区分真实簇与随机噪声的下限；二是取2倍维度（2dim），适用于较大数据集。对场景削减这种样本量上千、维度已经降过维的情况，我一般先用2dim，再根据噪声比例微调。如果数据维度高、样本又稀疏，minPts要取小一些，比如6到10，否则大部分点都会被判成噪声。

eps不能拍脑袋，要用K-距离曲线来定：对每个点，找它到第minPts个最近邻居的距离，把所有距离升序排列画曲线。曲线从平缓到陡峭的拐点处的纵坐标，就是比较合理的eps。拐点之前的区域里，绝大多数点周围密度都差不多；一旦越过拐点，说明进入了稀疏区，距离迅速拉大，用那里的距离做邻域半径，就能把稀疏点排除成噪声。

2.3 用K-距离曲线确定eps的MATLAB脚本

假设已经完成了降维，feature是N×d的降维特征矩阵，画K-距离曲线的代码如下：

distMat = pdist2(feature, feature); minPts = 2 * size(feature, 2); % 经验取值，可后续微调 sortedDist = sort(distMat, 2); kDist = sortedDist(:, minPts); % 每个点到第minPts近邻的距离 [kDistSorted, ~] = sort(kDist); plot(kDistSorted, 'LineWidth', 1.5); grid on; xlabel('样本序号'); ylabel(['第', num2str(minPts), '近邻距离']);

注意，distMat对角线是0，sort之后第1列永远是0，所以这里其实取到了包括自身在内的前minPts个点，相当于minPts-1个真实邻居。工程上影响不大，如果想要严谨，可以用sortedDist(:, minPts+1)。我两种都试过，曲线形状几乎一致，拐点处差异不明显，不用过于纠结。

画完图找拐点：曲线先是一段接近水平的段，然后在某个位置突然翘头。翘头起点对应的y值就是eps的候选值。实在拿不准，就在这个值附近上下浮动20%各跑一遍聚类，观察簇数和噪声占比的变化，取一个簇数稳定、噪声占比不超过20%到30%的值。

3. MATLAB实现流程：从场景生成到概率分配

这一章是全文的核心，我把完整流程按顺序拆开，每一步都给出可复现的代码和需要注意的细节。代码我尽量保持精简，去掉和核心逻辑无关的修饰，直接能跑。

3.1 风电-负荷联合场景生成

场景削减的前提是先有一批原始场景。真实项目里可以用风电场历史出力数据加上预测误差抽样生成，也可以像我这次一样，先用风速模型模拟一组出力曲线来验证方法。

风速通常用两参数Weibull分布建模，然后经过风电机组功率曲线转换成出力。用wblrnd抽样得到N×T的风速矩阵，功率曲线按三段式简化：低于切入风速、高于切出风速出力为零，切入到额定风速之间线性爬升，额定风速以上满发。这种简化曲线对聚类算法验证完全够用，实际工程中替换成厂商功率曲线查表即可。

负荷场景更简单，以一条典型日负荷曲线为基线，叠加正态分布随机波动，用来模拟预测误差。具体代码如下：

rng(2024); N = 1000; T = 24; % 风速与风电出力场景 k = 2.2; c = 8.0; v = wblrnd(k, c, N, T); vci = 3; vco = 25; vr = 12; Pr = 1.5; Pw = zeros(N, T); idx = (v >= vci) & (v < vr); Pw(idx) = Pr .* (v(idx) - vci) / (vr - vci); idx = (v >= vr) & (v <= vco); Pw(idx) = Pr; % 负荷场景 baseLoad = 30 + 8*sin((0:23)/24*2*pi) + 4*sin((0:23)/24*4*pi); Pl = repmat(baseLoad, N, 1) + randn(N, T) * 2; % 合并为联合场景，每条场景是48维 scenes = [Pw, Pl];

代码里baseLoad刻意构造了“早晚双峰、夜里有低谷”的典型日负荷形状。randn*2表示负荷预测误差标准差约2MW。把风电和负荷拼成联合场景，是为了聚类时同时考虑二者的联合分布，保留风荷之间的耦合关系，这在风-荷相关性强、需要做联合随机优化的场景下格外重要。代价是维度变成48，所以下一步必须降维。

3.2 标准化与PCA降维：别让量纲毁了聚类

48维时间序列直接丢进DBSCAN有两个问题：一是风电出力幅值（1.5MW级别）和负荷幅值（30MW级别）差异太大，欧氏距离完全被负荷分量主导，风电场景的差异被淹没；二是高维空间距离稀疏，密度聚类效果退化。

所以先标准化，再PCA降到低维。标准化我用z-score，逐列减去均值除以标准差，把量纲拉齐。PCA只保留累计方差贡献率超过85%的前若干主成分。这一步看起来只是“数据预处理”，其实是整个流程里最关键的一步。降维做不好，后面的DBSCAN基本就是给白噪声聚类。

scenesStd = (scenes - mean(scenes)) ./ std(scenes); [coeff, score, ~, ~, explained] = pca(scenesStd); dim = find(cumsum(explained) >= 85, 1); feature = score(:, 1:dim); fprintf('保留主成分个数: %d, 方差贡献率: %.2f%%\n', ... dim, sum(explained(1:dim)));

实际跑下来，48维联合场景通常前4到6个主成分就能贡献90%以上的方差，降维效果明显。如果对可解释性有要求，也可以提取几个手工特征替代PCA，比如风电日均值、峰值、峰谷差、负荷峰值时刻等，结果类似，但PCA更省事、更稳定。

3.3 DBSCAN聚类与代表场景提取

聚类直接调用MATLAB的dbscan函数，这个函数从R2019a版本开始在Statistics and Machine Learning Toolbox里提供。输入降维特征矩阵、eps和minPts，返回每个样本的簇编号，-1表示噪声：

minPts = 2 * dim; eps = 0.35; % 由K-距离曲线确定，这里是示例值 idxCluster = dbscan(feature, eps, minPts); numClusters = max(idxCluster); noiseCount = sum(idxCluster == -1); fprintf('簇数: %d, 噪声点: %d\n', numClusters, noiseCount);

聚类完成后，每个簇里的场景是一组在特征空间中位置相近的原始场景。接下来要做的是从每个簇里挑出一条代表场景，并计算它继承的概率。

代表场景有两种选择思路。一种是取簇内所有场景各时段的均值，好处是曲线平滑，坏处是这条“平均曲线”很可能不是任何一个可能发生的场景，放到调度问题里会出现“期望出力”这种现实中不发生的情况。另一种是从簇内挑离簇中心最近的真实场景，保证代表场景是可实现的，我强烈推荐这种。代码如下：

K = numClusters; rngScenes = zeros(K, 2*T); prob = zeros(K, 1); for i = 1:K members = find(idxCluster == i); memberScenes = scenes(members, :); center = mean(memberScenes, 1); distToCenter = sum((memberScenes - center).^2, 2); [~, idxMin] = min(distToCenter); rngScenes(i, :) = memberScenes(idxMin, :); prob(i) = length(members) / N; end

注意这里是在原始场景空间里计算到簇中心的距离，而不是在PCA特征空间。虽然聚类是在降维空间完成的，但代表场景要在原始空间里生成才有物理意义。PCA是线性投影，原始空间里接近簇中心的位置和特征空间通常也接近，结果差别不大。实在不放心，也可以把代表场景索引映射到特征空间做，结果几乎等价。

3.4 概率分配与噪声点处理

概率分配本身很简单：每个簇的样本数除以总样本数。由于噪声点没有进入任何簇，这部分概率直接丢弃。这在场景削减里是合理的，因为我们把稀疏极端场景当作小概率异常剔除，不参与后续随机优化。

如果希望把这些极端场景也纳入分析，就不要丢掉概率，而是把噪声点单独归为一个特殊簇，或者把距离各个簇中心最近的噪声点各保留一条。我个人的习惯是：先看噪声点占比，控制在20%以内，这时丢掉它们对概率分布影响很小；如果占比超过30%，说明eps给太小了，回去调参数，而不是硬着头皮丢数据。

4. 算例实测与结果对比：和K-means差在哪

光讲原理和代码不行，得有实测结果支撑。这里分享一个标准算例的跑数过程，以及和K-means对比时看到的实际差异。

4.1 算例设置与聚类结果

我用上面的参数做了标准算例：1000条联合场景，每条48维，PCA保留前4维，minPts=8，根据K-距离曲线拐点取eps=0.35。跑出来的典型结果是5到15个簇，噪声点占比10%到25%，具体数值会随随机种子和负荷曲线形状变化。

聚类结果可以用散点图直观检查。把降维后的前两维特征画出来，不同颜色标记不同簇，一眼就能看出DBSCAN把数据分成了几坨形状不规则的区域，噪声点散落在稀疏地带。这种可视化检查比看任何数值指标都重要，建议每次跑完都看一眼。

4.2 削减效果怎么评估才靠谱

评估削减质量，最直接的办法是看削减前后场景集合的统计特性是否接近。对风电出力，先比均值曲线：

meanBefore = mean(scenes(:, 1:T), 1); meanAfter = prob' * rngScenes(:, 1:T); figure; plot(meanBefore, 'LineWidth', 2); hold on; plot(meanAfter, '--', 'LineWidth', 2); legend('削减前均值', '削减后均值');

如果两条曲线几乎重合，说明期望水平保留得好。接着看10%和90%分位数边界，检查波动范围的覆盖情况。数学上更严谨的指标是Kantorovich距离，它度量两个概率分布之间的最优传输代价，是场景削减论文里最常用的标准。不过很多工程同学一看到这个名词就头疼，我的建议是：写论文、做对比实验时用Kantorovich距离；工程验证阶段，均值曲线、分位数覆盖和下游优化目标值这三件套就足够说明问题了。

4.3 和K-means对比时的三个实际差异

我把同批数据用K-means跑了一遍，K值分别试了5、8、12。最大的差别有三个。

第一，K-means把一些低风速向高风速过渡的“爬坡场景”强行平均成一条平缓曲线，代表场景失真；DBSCAN靠密度聚类，爬坡场景如果数量足够会单独成簇，保留住了陡峭的爬坡特征，这对后续模拟机组爬坡约束非常重要。

第二，K-means对离群点敏感，簇中心被极端场景拖拽，导致代表场景整体偏高；DBSCAN把极端场景标成噪声，主流簇的中心干净得多，代表场景更贴近大多数场景的真实形态。

第三，K-means的K需要反复试，而DBSCAN只需要调eps和minPts，调参逻辑更贴近数据的自然分布。当然K-means计算速度快，样本量特别巨大的时候仍有价值；DBSCAN的距离矩阵计算复杂度是O(N^2)，N超过5000时会有点吃力，这个后面聊优化。

5. 常见问题与避坑指南：这些坑我替你踩过了

最后这章是最想写给同行看的。DBSCAN做场景削减，代码跑通不难，跑出可靠的结果才是真功夫。以下问题我在实际项目里都遇到过，每一个都对应着具体调整方法。

5.1 噪声点太多，代表场景不够用

这是DBSCAN场景削减里最常翻车的点。我一开始把eps定得很小，1000条场景里500条全是噪声，最后只剩三五个簇，削减后的场景集根本没法用。排查思路很明确：先画K-距离曲线，看看拐点是不是被噪声拖得很早；再检查降维后是不是仍有大量孤立点。噪声占比超过30%，优先调大eps，而不是去调minPts。minPts变大只会让核心点条件更严格，噪声只会更多。

5.2 高维时序场景的稀疏问题

不降维直接跑DBSCAN基本必翻车。48维空间里，任意两个样本的欧氏距离都差不多大，密度概念失去意义，聚类结果随机性极强。我见过有人拿24维风电场景硬跑，画出来的簇毫无规律，换了降维之后马上正常。解决办法就是第3章说的PCA或特征压缩，把维度降到4到6维，DBSCAN的性能立刻恢复。

5.3 没有Statistics Toolbox怎么办

dbscan函数需要Statistics and Machine Learning Toolbox，旧版本或精简版MATLAB可能没有。解决途径有三条：一是升级到带工具箱的新版MATLAB，新版dbscan函数做得挺完善，文档里也支持多种距离度量；二是自己动手实现经典DBSCAN，算法本身不复杂，双重循环加邻域搜索，几十行内能写完，网上开源实现也很多；三是改用K-means++或谱聚类做过渡，但那样又回到预设簇数的问题。如果你只有基础MATLAB环境，我建议直接找一份DBSCAN公开实现改一改距离计算就行。

5.4 场景间距离尺度不统一

这个问题说起来谁都懂，但实操中很多人都栽过：风电和负荷单位不同、幅值不同，直接拼在一起聚类，距离贡献完全被负荷主导。处理方式是先标准化，或者给两个分量分别设权重。如果不想动PCA，也可以在聚类前对风电分量和负荷分量分别归一化后再拼接。这个坑比较隐蔽，因为聚类结果“看起来”能跑出来，但代表场景里的风电特征明显失真。