做风电和负荷场景缩减，很多人第一反应就是用K-means或者层次聚类，但实际处理过预测场景集的人应该都有体会：K-means得提前给定簇数，而且遇到那些离群的极端场景，会被强行拉进某个簇里，最后缩减出来的“典型场景”既不典型，也不极端，反而把调度决策需要的边界信息给抹掉了。
这篇博文聊一聊我最近在Matlab里实现的一套基于DBSCAN密度聚类的风电-负荷确定性场景缩减方法。整套方法解决的核心问题是：给定一组由风电功率和负荷功率组合而成的确定性场景集合，如何在尽量保留原始概率分布特征的前提下，用少量几个代表场景替代完整场景集，从而减轻后续随机优化、风险评估的计算压力。文章会从问题建模、算法原理、Matlab代码实现、参数调优到效果评估完整走一遍，适合正在做电力系统随机规划、新能源消纳分析或者搞场景生成与缩减的同学参考。
1. 场景缩减到底在解决什么问题
1.1 为什么电网分析需要“场景”这个概念
风电和光伏出力具有天然的不确定性，负荷侧也有很强的日内波动和季节差异。在研究含高比例新能源的电力系统调度策略时，不能只拿一条预测曲线做决策，因为预测误差和波动性本身就会直接影响备用容量配置、储能充放电策略以及机组组合结果。
工程上常见的做法是场景法：通过历史数据统计、蒙特卡洛采样或者概率预测模型，生成大量可能出现的风电-负荷组合场景，每个场景都是一条完整的时间序列，代表一种可能出现的情况。比如某个地区风电出力有96点日曲线，负荷也有96点日曲线，那么一个“组合场景”就是一个192维的向量——前96维是风电，后96维是负荷。
场景数量一多，问题就来了。随机规划、鲁棒优化这些方法在求解时，计算复杂度基本随着场景数量线性甚至超线性增长。如果初始场景集有2000个，直接拿去优化，内存和时间开销都非常可观。场景缩减就是在保证决策精度不显著下降的前提下，把场景集压缩到几十个以内，让优化模型能真正跑得动。
1.2 确定性场景集合与概率性场景的区别
这里需要先厘清概念。所谓“确定性场景”，指的是每个场景本身是一条确定性的时序曲线，场景与场景之间的差异体现了不确定性。与之相对的是模糊场景或区间场景，描述的是一组可能取值范围的包络。
本文处理的对象就是前者：一堆已经生成好的、确定性的风电-负荷组合曲线。我们可以把它们看成一个多维空间中的散点集合，每个点是一个场景，点与点之间的距离反映了场景之间的相似程度。场景缩减本质上就是一个聚类问题——把彼此相似的场景归为一类，从每一类中挑出最有代表性的一个场景，用代表场景和对应概率去近似整个场景分布。
用生活化的方式理解，这件事类似于人口普查中的“抽样调查”。不挨家挨户问所有人，而是把人口按年龄、区域、职业特征聚成几类，再选几个典型样本，用他们的统计数据代表一大群人的情况。只要分类合理、代表选得好，整体统计特征就能保住。
1.3 场景缩减的核心评价标准
评价一个场景缩减方法好不好，一般看三个维度：
第一是分布保真度。缩减前后场景集合的统计特征应该尽量一致，通常用均值、方差、分位数以及概率分布距离来度量。
第二是代表性场景的质量。每个代表场景应当真正反映簇内场景的共同形状，不能是簇内场景的简单平均“糊”出来的一条无意义曲线，更不能选到离群点上。
第三是概率重估的合理性。缩减后每个场景都带一个概率权重，这个权重应该与原始场景集中该类场景出现的频率一致。
传统K-means能解决前两个维度的一部分问题，但第三个维度它天然不擅长——因为离群点会被强制归类，导致某些簇的概率权重虚高。而DBSCAN因为支持噪声点判定，能更真实地刻画场景分布形态，这也是我在本次实现中选择它的根本原因。
2. DBSCAN用于场景聚类的原理与选型依据
2.1 K-means、层次聚类与DBSCAN的直观对比
先看一个实际场景的三维示意。假设我们有500个风电-负荷组合场景，投影到二维特征空间后，大致呈现三个密集区域和十几个零散点。如果用K-means，K值设成3，算法会把零散点也分配到三个簇里，簇的边界会被拉得很奇怪，质心位置也会被零散点拖偏。如果K值设置不当，比如设成10，又会在密集区域内过度切分，产出大量相似度极高的冗余代表场景。
层次聚类不用预设K值，但计算量相对大，而且它的分裂/聚合过程对距离定义比较敏感。更重要的是，层次聚类不提供“噪声点”概念，所有场景最终都会被纳入某个簇。从场景缩减的角度说，有些极端场景恰恰是不该纳入任何典型簇的，它们应该被保留下来作为独立场景，代表极端工况。
DBSCAN的核心思想是基于密度连通性。它把场景空间划分为三类点：核心点、边界点和噪声点。落在高密度区域的场景被聚成一簇，边缘稀疏区域的场景被标记为噪声。这种处理方式更贴近场景缩减的真实需求——我们既要典型场景，也要保留一定数量的极端场景作为边界情况，而不是把所有情况都平均化。
2.2 DBSCAN的数学定义与聚类流程
DBSCAN有两个关键参数：邻域半径ε和最小邻域点数MinPts。对于任一场景点p，定义：
Nε(p) = {q ∈ D | dist(p, q) ≤ ε}
这是p点在半径ε内的邻域点集合。如果|Nε(p)| ≥ MinPts，点p就是核心点；如果某个点不是核心点，但落在某个核心点的ε邻域内，它是边界点；如果一个点既不是核心点，也不属于任何核心点的邻域，它被标记为噪声点。
聚类形成的过程是：任意选一个未被访问的核心点，把它邻域内所有密度相连的点归入同一个簇，不断扩张，直到没有新的点可以加入，再从未访问的点中找下一个核心点开始新的簇。
参数的含义也很好理解。ε决定了场景之间“相似到何种程度才算同一类”，MinPts则控制了一个簇至少要有多密集才能被承认。这两个参数配合起来，相当于在问：半径ε内至少要有MinPts个场景，我们才认为这一带确实存在一种典型模式。
2.3 为什么DBSCAN天然适配风速场景结构
风电出力场景有一个很明显的物理特征：受大气活动影响，出力曲线往往呈现多个典型模态，比如“全天平稳大风”“白天弱风夜晚强风”“爬坡-高峰-陡降”等。不同模态之间的过渡区场景数量比较少，密度低。模式内部因为气象条件的小幅波动，场景彼此非常相似，密度高。
这种“高内聚、低耦合、中间带稀疏”的数据结构，正好是DBSCAN最擅长处理的类型。K-means假设簇是凸形的，对于非凸或不规则形态的模态无能为力；而DBSCAN不关心簇的形状，只要密度足够，就可以聚成任意形态的簇。
此外，场景缩减往往还希望保留一些信息，即哪些场景是不典型的边缘情况。这些边缘场景可能对应极端风速、极端负荷，对调度决策中的备用安排和切负荷风险分析有重要价值。DBSCAN的噪声点标记给了我们一个明确的口子，让这些场景保留下来，而不是被粗糙地平均掉。
2.4 距离度量的选择：欧氏距离与DTW的取舍
DBSCAN需要一个距离定义。对场景时序曲线，最朴素的是逐点欧氏距离。
dist(p, q) = sqrt(Σ(p(t) - q(t))²)
欧氏距离的优点是计算简单、可解释性强，Matlab里pdist2一行就能算完。但对时序场景来说，它有一个问题：两段形状相同但在时间轴上错位的曲线，欧氏距离会很大，被判定为不相似。比如风电出力峰值的时刻差了一小时，实测场景和预测场景可能本质上是同一种天气过程，但欧氏距离会把它们分开。
更精细的做法是用动态时间弯曲（DTW）距离。DTW允许时间轴非线性对齐，能更好地度量形状相似性。但DTW的计算复杂度高，对场景数量大的场景集，构建距离矩阵会非常耗时。我的经验是：如果场景的时间分辨率较粗（如小时级、96点日曲线），且风电出力的峰谷时移不那么明显，欧氏距离完全够用；如果是分钟级采样且需要精细刻画爬坡形状，建议考虑DTW。本次实现以欧氏距离为主，后面会提到如何在此基础上扩展。
3. Matlab代码实现与关键环节解析
3.1 数据准备与预处理
动手写代码前，第一步是把场景数据整理成标准矩阵格式。假设我们有N个场景，每个场景由T_w个风电时刻和T_l个负荷时刻组成，数据矩阵的尺寸就是N×(T_w+T_l)。在Matlab中，用矩阵按行存储场景，每行是一个完整场景，列是各时间断面的功率值。
有一个细节不要忽略：风电出力和负荷功率在数值量级上可能差很多倍。比如某地区风电装机容量500MW，负荷峰值5000MW，如果不做归一化，聚类时欧氏距离会被负荷分量主导，风电场景的差异几乎不影响聚类结果，这样缩减出来的典型场景对风电的刻画就是糊的。
所以必须先做归一化。常用的做法是min-max归一化，把每个维度的数据映射到[0,1]区间：
norm_data = (data - min(data)) / (max(data) - min(data))
这里要注意max和min是按列计算的。归一化之后再做聚类，得到的簇结构才能真正同时反映风电和负荷两个维度的相似性。聚类完成后，在选择代表场景时是基于归一化空间选的，但输出结果时要映射回原始量纲。
3.2 DBSCAN主函数的实现思路
Matlab从2019b版本开始，Statistics and Machine Learning Toolbox内置了dbscan函数，可以直接调用。但如果想灵活控制聚类细节，或者需要输出每个点的核心点/边界点/噪声点身份信息，自己实现一遍也很有价值。我先把内置函数用法说清楚，再给出一个可运行的自主实现。
内置函数用法很简洁：
[idx, corepts] = dbscan(norm_data, eps, minpts);其中idx是每个样本的簇标签，-1表示噪声点；corepts是逻辑数组，1表示核心点，0表示非核心点。这个函数支持'Distance'参数，可以指定距离度量方式，比如'euclidean'、'squaredeuclidean'等。
下面是我实现的精简版DBSCAN，核心逻辑与原始论文一致：
function [clusterIdx, isCore] = dbscanCustom(data, eps, minPts) % data: N x D 矩阵, 每行一个场景 % eps: 邻域半径 % minPts: 最小邻域点数 n = size(data, 1); D = pdist2(data, data); clusterIdx = zeros(n, 1); isCore = false(n, 1); clusterId = 0; for i = 1:n if clusterIdx(i) ~= 0 continue; end neighbors = find(D(i, :) <= eps); if numel(neighbors) < minPts clusterIdx(i) = -1; % 临时标记为噪声 continue; end clusterId = clusterId + 1; clusterIdx(i) = clusterId; isCore(i) = true; seedSet = setdiff(neighbors, i); while ~isempty(seedSet) j = seedSet(1); seedSet(1) = []; if clusterIdx(j) == -1 clusterIdx(j) = clusterId; end if clusterIdx(j) ~= 0 continue; end clusterIdx(j) = clusterId; jNeighbors = find(D(j, :) <= eps); if numel(jNeighbors) < minPts continue; end isCore(j) = true; for k = 1:numel(jNeighbors) if clusterIdx(jNeighbors(k)) == 0 seedSet = [seedSet; jNeighbors(k)]; end end end end end这段代码的逻辑分三层。外层循环扫描所有点，遇到未归类点就开始新的簇。第二层是核心点的判断，如果邻域点数不够，先标记为噪声，后续如果被某个核心点邻域包含，再改标为边界点。第三层是种子扩张循环，用队列方式不断吸收新的密度连通点，直到该簇完全生长完毕。
有一个细节值得注意，这种实现方式将距离矩阵一次性计算并存储。当场景数量到5000以上时，距离矩阵会占用约5000×5000×8字节的内存，大约200MB，在普通电脑上虽可运行但已偏卡。如果场景规模到数万级别，推荐用分块计算距离或在需要时动态计算邻域点，不要一次性放全矩阵。
3.3 核心码段的参数敏感性解析
实现DBSCAN不难，真正让人头疼的是ε和MinPts的选取。MinPts有一个业界常用的经验下限：对于D维数据，MinPts应不小于D+1。如果场景维度是192，那MinPts至少是193？显然不合理——场景分类和低维点聚类不同，数据点虽处于高维空间，但实际流形结构往往很低维。因此实践中一般取MinPts = 2到10之间，我通常取5，然后通过ε调节簇的粒度。
ε的取值决定了聚类的尺度。ε太小，大量点会成为噪声，场景被过度切碎；ε太大，所有场景融成一团，缩减失去意义。有一种经典的参数标定方法是k距离图：
- 计算所有场景到其第k个最近邻点的距离；
- 将距离值降序排列，画成折线图；
- 图中斜率骤变的拐点位置，就是ε的合理候选值。
Matlab实现如下：
k = minPts - 1; [~, distK] = knnsearch(norm_data, norm_data, 'K', k+1); kDist = sort(distK(:, end), 'descend'); plot(1:length(kDist), kDist);画出来之后，观察曲线由陡变缓的位置，对应横坐标约在总点数的某个百分比处，纵坐标即为ε的候选值。实际使用时，我会在这个值附近左右各扩展2~3个候选值，分别聚类，对比结果，取场景数量更符合需求的ε。
3.4 代表场景提取与概率权重分配
聚类完成后，每个簇内有一批相似场景。下一步是从每个簇中选出一个代表场景，并计算该簇对应场景的发生概率。
代表场景的选择有两种常见策略。第一种是簇内平均法：把簇内所有场景逐点取均值，得到一条平均曲线。优点是平滑、噪声少，缺点是可能产生一条物理上不够“真实”的曲线，比如把峰值削平了。第二种是最近中心点法：先计算簇内所有点的平均向量，然后从簇中选择与该平均向量欧氏距离最近的真实场景作为代表。这种方法保证代表场景是一个实际存在的、物理上可出现的场景，在工程调度里更具参考价值。
我倾向于第二种。因为后续要拿代表场景去算机组组合或储能策略，一个现实中真实发生过（或真实预测过）的场景，更容易被工程人员接受，也便于追溯来源。
概率权重的分配则简单很多。统计每个簇内的场景数量占总场景数的比例，作为该代表场景的缩减后概率。注意噪声点怎么办？噪声点没有被归入任何簇，它们不应当被直接丢弃。我通常给每个噪声点单独分配一个很小的概率权重（如总概率的0.5%均分），或者将它们单独归入一个“极端场景集”，在后续分析中单独考察。不处理直接丢弃会导致风险低估。
% 输入: clusterIdx, norm_data, original_data % 输出: repScen(代表场景), repProb(概率权重) clusters = unique(clusterIdx); clusters = clusters(clusters ~= -1); nClust = numel(clusters); repScen = zeros(nClust, size(original_data, 2)); repProb = zeros(nClust, 1); totalPts = size(clusterIdx, 1); for i = 1:nClust members = find(clusterIdx == clusters(i)); centroid = mean(norm_data(members, :), 1); [~, nearestIdx] = min(sum((norm_data(members, :) - centroid).^2, 2)); repScen(i, :) = original_data(members(nearestIdx), :); repProb(i) = numel(members) / totalPts; end这里有个容易踩的坑：在寻找簇内最近点时，必须用归一化空间的距离。如果用原始空间的均值去算最近点，会再次被量纲差异支配。上述代码在归一化空间里计算质心并选点，但输出的代表场景用原始数据，规避了量纲问题。
4. 场景缩减完整流程与效果评估
4.1 从生成场景到缩减落地：一整条可复用的流水线
我按实际项目经验整理了一条完整可执行的操作流程，供参考：
第一步，数据装载。从Excel或CSV读取场景数据，用matlab的readmatrix函数一次性读入。注意列顺序，如果要拼接风电和负荷，务必保证各场景序列长度一致。
第二步，结构组装与归一化。将风电矩阵和负荷矩阵水平拼接成一个总矩阵，然后做按列min-max归一化，记录归一化参数以便后续还原。
第三步，参数标定。设定MinPts（推荐从5开始），画k距离图，用拐点法确定ε候选值。
第四步，聚类与评估。调用dbscan或自定义函数进行聚类，查看簇数、噪声比例、每簇点数分布。正常情况下期望得到5~20个非噪声簇，噪声比例在5%~20%之间。如果噪声比例过高，多半是ε太小；如果簇数太少，可能是ε太大或数据本身确实高度集中。
第五步，代表场景提取。按最近中心点法提取每个簇的代表场景，计算概率权重。
第六步，可视化与对比。把原始场景集、聚类结果、缩减场景集画在同一张图上，直观观察缩减效果。
4.2 缩减效果的物化指标
量化评估缩减效果不能只靠肉眼看图，还需要数值指标支撑。我在实际工作中常用以下三个指标：
第一个是均值绝对误差MAE。计算缩减前后场景集的逐时刻平均功率差，评估整体期望是否被保住。
第二个是标准差误差。对比缩减前后，每个时间断面功率标准差的变化，这个指标能直接看出波动特征是否失真。如果缩减后标准差明显缩小，说明聚类过程把多样性抹掉太多了。
第三个是K-L散度。把场景集按功率区间统计成直方图，再计算缩减前后概率分布的K-L散度。这个指标更严格，适合学术出图。
K-L散度计算公式为：
D_KL(P||Q) = Σ P(x) log(P(x) / Q(x))
实际操作中，P取原始场景集的离散分布，Q取缩减后场景集的分布，值越小说明分布保留越好，通常小于0.1就可以认为缩减质量良好。
4.3 完整的Matlab示例：50个场景缩减到6个典型场景
下面用一个实际案例把整个流程串起来。假设我们有一组包含40个风电场景和40个负荷场景的数据，组合后为50个192维的联合场景，目标缩减到约5~8个代表场景。
代码片段如下：
% 1. 加载数据 data = readmatrix('wind_load_scenarios.csv'); windData = data(:, 1:96); loadData = data(:, 97:192); % 2. 拼接并归一化 rawScen = [windData, loadData]; n = size(rawScen, 1); normScen = (rawScen - min(rawScen)) / (max(rawScen) - min(rawScen)); % 3. 参数标定 minPts = 5; k = minPts - 1; [~, distK] = knnsearch(normScen, normScen, 'K', k+1); kDist = sort(distK(:, end), 'descend'); figure; plot(1:n, kDist, 'LineWidth', 1.5); xlabel('点序号'); ylabel('第k近邻距离'); title('k距离图'); grid on; % 从图中读拐点, 假设 eps = 0.18 eps = 0.18; % 4. 聚类 clusterIdx = dbscan(normScen, eps, minPts); % 5. 提取代表场景与权重 % ...使用之前给出的代码块... % 6. 可视化 figure; plot(rawScen(:, 1:96)', 'Color', [0.7 0.7 0.7]); hold on; plot(repScen(:, 1:96)', 'r-', 'LineWidth', 1.5); xlabel('时刻'); ylabel('风电功率/MW'); title('风电部分: 原始灰色, 代表场景红色');这个案例跑完后，预期得到如下的典型场景集：几个簇对应“平稳大风”“夜间大风白天小风”“白天爬坡晚高峰”“日均低出力”等典型模式，外加可能1~2个噪声场景被单列。整体压缩比10:1左右，K-L散度能控制在0.05~0.1之间，均值和方差误差在个位数百分比内。
4.4 可视化分析：如何判断聚类结果是“真的合理”
聚类结束后，第一件事不是看指标，而是画图。
我在每个聚类迭代后必画三张图。第一张是原始场景曲线叠加图，全部灰度线条，直观感受数据的分散程度。第二张是聚类着色图，按簇标签用不同颜色画曲线，观察簇之间是否重叠严重。第三张是代表场景与簇内成员对比图，选一两个簇放大，看代表场景是否贴合簇内曲线的整体走势。
如果第三张图中代表场景明显偏离簇内曲线的中间带，那说明最近中心点法出了问题，多半是归一化参数不对。如果第二张图中不同簇的曲线大面积交织，说明ε取太大或MinPts太小，需要调参。可视化是调参最直接的反馈，比任何指标都灵敏。
5. 常见问题、避坑指南与工程拓展
5.1 数据量纲差异把聚类结果带偏
这个问题出现的频率极高。风电和负荷量级不同，甚至负荷在夏季和冬季的峰谷值也差异巨大。如果不先归一化，DBSCAN的ε参数会被量纲大的维度带走。
判定方法：聚类完成后，分别计算每个簇在风电维度上的均值和在负荷维度上的均值，看它们是否都各自紧凑。如果负荷维度紧凑而风电维度松散，说明归一化环节没有做好。
解决方式：对每一列分别做归一化；更稳妥的是采用z-score标准化，即减去均值除以标准差，这样可以消除离群极值的影响。
5.2 MinPts取值对结果的影响机理
MinPts调大，聚类标准变严，更少的点会被认作核心点，簇数可能减少。MinPts太大（比如场景总数的一半），会导致几乎所有点都成噪声。MinPts调小，簇数增大、噪声减少，但可能把一些偶然聚集的点也识别成簇。
实操中，我先固定MinPts=5跑一遍，看噪声比例，再尝试MinPts=3和7。多数场景下，结果变化不大就说明参数区间合理；如果变化剧烈，说明数据本身密度不均匀，可能需要在预处理环节先做一次粗筛或降维。
5.3 高维场景距离“被稀释”的问题
当场景维度很高，比如192维，欧氏距离的区分度会下降，因为高维空间中所有点对距离趋于接近。这是维数灾难的一个直接体现。
一个有效的缓解手段是先用主成分分析（PCA）降维到10~20维再做聚类。这样做有两个好处：一是距离度量更稳定，二是计算速度提升。代价是损失少量信息，但对于场景缩减这种应用，主成分保留的前若干维度已经解释了大部分风电-负荷波动特征，实际影响很小。
Matlab里的PCA调用很简单：
[coeff, score, ~, ~, explained] = pca(normScen); reducedScen = score(:, 1:15);然后在reducedScen上做DBSCAN，聚类效果通常比在高维原空间更好。
5.4 噪声场景的处置策略
噪声点是最容易被初学者忽略的部分。注意，这里的“噪声”不是数据坏点，而是那些与所有典型模式都不相似的边缘场景。强行丢弃会导致缩减后的场景集丢失尾部风险信息，在调度中表现为备用容量不足。
我的经验是：将噪声场景保留，赋予均匀分布的小概率权重，或者直接单独列出“低概率但高影响场景集”，专门用于可靠性校核。这两种做法都比丢弃更稳妥。具体操作中，我会先看噪声点数量占总场景数的比例，如果超过25%，优先考虑调整参数，因为过高的噪声比例往往意味着聚类并没有捕捉到数据的真实结构。
5.5 进阶思路：考虑时序相关性的改进方向
以上实现基于逐点独立的距离度量。如果想把风电出力的时序相关性也考虑进来，可以考虑在距离矩阵中引入趋势惩罚项，或者改用DTW距离。另一个方向是在聚类前对风速场景做波动率特征提取，将特征向量加入聚类输入，再配合图示化分析，效果会更精细。
考察DTW方案时，可以先小规模试验，因为其计算复杂度高。若场景数在500以内，计算DTW距离矩阵尚可接受；若超过1000，建议先用PCA或特征提取降低维度，再考虑DTW。实际项目中我也见过有人用分段聚合近似（PAA）先把时序压缩为若干个时间段的均值，再跑DBSCAN，效果也很不错，计算成本大为降低。
5.6 与后续随机优化的衔接要点
场景缩减的最终目的是服务后续优化。缩减后的场景集不只是一个聚类结果，它需要以特定格式输入优化模型。其中最重要的是给每个代表场景配上概率权重，权重之和必须为1。如果保留了噪声场景，记得在权重分配里包含它们，不能只对簇内场景做归一化。
另外，代表场景的数值应当还原到原始量纲并检查合理性，比如风电功率不能为负、负荷不能超过区域峰值上限等。因为在归一化空间选完代表点后，直接映射回原始空间可能引入微小的越界情况，做一次物理约束校核是必要的。
6. 完整代码汇总与使用说明
为了方便直接拿去用，我把整套流程整理成一个脚本模板。用户只需替换数据文件路径，并按自己的场景维度和目标簇数微调几个参数即可。核心结构保持清晰，模块间用注释分隔，命中即跑。
%% 基于DBSCAN的风电-负荷确定性场景缩减 % 输入: wind_load_scenarios.csv, 每行一个场景, 前96列风电, 后96列负荷 % 输出: repScen(代表场景矩阵), repProb(概率权重向量), clusterIdx(聚类标签) clear; clc; close all; % 加载 data = readmatrix('wind_load_scenarios.csv'); windData = data(:, 1:96); loadData = data(:, 97:192); rawScen = [windData, loadData]; % 归一化 minVal = min(rawScen, [], 1); maxVal = max(rawScen, [], 1); normScen = (rawScen - minVal) ./ (maxVal - minVal); % 画k距离图, 手动确定eps minPts = 5; [~, distK] = knnsearch(normScen, normScen, 'K', minPts); kDist = sort(distK(:, end), 'descend'); figure; plot(1:length(kDist), kDist, 'LineWidth', 1.5); grid on; title('k距离图'); eps = 0.18; % 根据拐点调整 % DBSCAN clusterIdx = dbscan(normScen, eps, minPts); % 统计信息 nClust = max(clusterIdx); nNoise = sum(clusterIdx == -1); fprintf('簇数: %d, 噪声点数: %d\n', nClust, nNoise); % 代表场景提取 clusters = unique(clusterIdx); clusters = clusters(clusters ~= -1); nC = numel(clusters); repScen = zeros(nC, size(rawScen, 2)); repProb = zeros(nC, 1); for i = 1:nC members = find(clusterIdx == clusters(i)); centroid = mean(normScen(members, :), 1); [~, nearestIdx] = min(sum((normScen(members, :) - centroid).^2, 2)); repScen(i, :) = rawScen(members(nearestIdx), :); repProb(i) = numel(members) / size(rawScen, 1); end % 可视化 figure; subplot(2,1,1); plot(rawScen(:, 1:96)', 'Color', [0.75 0.75 0.75]); hold on; plot(repScen(:, 1:96)', 'r-', 'LineWidth', 1.5); xlabel('时刻'); ylabel('风电/MW'); title('风电出力场景缩减'); subplot(2,1,2); plot(rawScen(:, 97:192)', 'Color', [0.75 0.75 0.75]); hold on; plot(repScen(:, 97:192)', 'b-', 'LineWidth', 1.5); xlabel('时刻'); ylabel('负荷/MW'); title('负荷场景缩减'); % 输出缩减结果 csvwrite('reduced_scenarios.csv', [repScen, repProb]);这个脚本中，k距离图和聚类结果的可视化是整个调参过程的核心反馈。每调整一次ε，建议把上面最后的可视化重新跑一遍，同时对比噪声比例和簇数的前后变化。多试几轮，你基本能把握住数据集的聚类尺度。
最后再分享一个调试技巧：先不要直接追求“缩减到多少个场景”，而是先用默认参数跑一遍，看看自然聚成几簇、噪声比例多少。多数情况下，数据会自然展现出3~10个明显模态。如果自然簇数和你的目标缩减规模差距不大，说明这套数据适合用DBSCAN做缩减；如果自然簇数远超目标场景数，可以考虑增大ε合并小簇，但要注意合并后代表场景的质量。如果自然簇数太少，则考虑减小ε进一步细分，或者先用PCA等方法提取更多特征再聚类。总之，参数调整的每一步都应该基于可视化反馈去判断，而不是盲目追求某个数字指标。
