做电力系统不确定性分析的人几乎都绕不开“场景削减”这道坎。风电和负荷的历史数据动辄一年8760个小时直接塞进随机优化模型里求解规模大得吓人内存和计算时间根本扛不住。于是大家普遍的做法就是从海量历史样本里挑出少数几个有代表性的“场景”让模型既能保留不确定性信息又能在可接受的时间内算完。我最近在MATLAB里用DBSCAN密度聚类做了一版风电-负荷场景削减方法实测下来比传统K-means顺手很多尤其是在处理异常样本和带状分布数据时优势非常明显。这篇文章就把我的实现思路、参数调优过程和踩过的坑完整记录下来给同样在做场景削减的朋友一个可直接参考的样板。1. 为什么是DBSCAN从场景削减的本质需求说起1.1 场景削减到底在削什么先理清概念。所谓“场景”本质上是风电出力或者负荷功率在时间维度上的一组采样序列比如一天24个时刻点或者更细粒度到96个点每组序列就是一个场景。历史数据里可能包含几千甚至上万个这样的序列它们之间高度相关、彼此相似直接全量使用既浪费资源还会让随机规划的复杂度爆炸。场景削减的目标就是从原始场景集中选出一批数量少得多、但概率特性接近原始分布的“代表场景”让后续的优化计算例如机组组合、经济调度不必再背负全部场景的负担。聚类是完成这件事的主流手段之一而聚类的核心问题就变成了怎么把相似场景归到同一类再为每一类挑一个代表点并计算这个代表点出现的概率。同样做聚类选的算法不同结果天差地别。1.2 聚类算法选型DBSCAN的直觉与优势最早做场景削减时我第一反应是用K-means因为经典、简单MATLAB自带kmeans函数一行代码就能出结果。但用久了会发现一个很尴尬的问题K-means需要你事先指定K值也就是聚类个数而在一条实际的风电数据面前谁也没法准确说出“到底应该分几类”。你拍脑袋定一个K等于8跑出来的轮廓系数也许还行但换一天的数据K等于6可能更合理。K-means对初始聚类中心极其敏感跑十次可能得到十种结果。DBSCAN最打动我的一点是它不需要指定聚类数量。聚类个数完全由数据密度分布自动决定它天生能识别任意形状的簇不像K-means那样假设簇是凸的、球状的。风电和负荷场景在特征空间里的分布往往不是规整的圆团可能是长条形、弧形甚至带空洞这种时候K-means强行按圆形边界划分会把原本属于同一类别的场景硬生生切开。除了形状DBSCAN还有一记杀手锏——它自带“离群点识别”能力。真实风电数据里总有那么几条异常场景比如某一天风机故障导致出力骤降、负荷数据里混入通信异常产生的跳变这些点在特征空间里往往孤悬在外围。K-means会把这些噪声也当作正常样本参与聚类逼着聚类中心往它们身上靠结果就是代表场景严重失真。DBSCAN则会把低密度区域的点直接标记为噪声场景削减时就自然把这些“脏样本”剔除出去剩下的场景质量高出一大截。1.3 DBSCAN的三个核心概念DBSCAN的全称是Density-Based Spatial Clustering of Applications with Noise基于密度的带噪声空间聚类。它的逻辑一句话就能概括在数据空间里如果一个点周围足够稠密那它就属于某个簇如果周围空荡荡那它就是噪声。这里展开三个关键概念Eps邻域以某个点为中心、Eps为半径的圆形区域。Eps越大能纳入邻域的点越多聚类越“松弛”。MinPts一个点要想成为核心点它的Eps邻域内至少要包含MinPts个点包含自身。MinPts决定了“密度”的下限。核心点、边界点、噪声点邻域内点数超过MinPts的是核心点不到MinPts但落在某个核心点邻域内的是边界点既不是核心点也不在任何一个核心点邻域内的就是噪声点。DBSCAN从一个随机点出发层层扩张核心点的邻域直到再也找不到新的核心点时一个簇就封闭了。这个过程不依赖任何全局参数如聚类个数也不要求簇的形状是凸的这是它和其他聚类算法最本质的区别。1.4 和K-means、层次聚类、GMM的对比做选型时我也认真比较过层次聚类和GMM高斯混合模型这里直接给结论算法是否需要指定K簇形状限制噪声处理稳定性适用性K-means是凸球状无噪声会干扰聚类中心依赖初始化结果波动数据简单时可快速上手层次聚类是较灵活一般不剔除确定性好小样本场景GMM是椭圆状软划分无显式噪声依赖初始化分布近似高斯时效果好DBSCAN否任意形状自动识别噪声参数稳定时结果固定场景数据含离群点和复杂形状时层次聚类虽然不假设凸形但它的计算复杂度通常在O(n²)量级样本上到几千之后矩阵合并开销极大。GMM对初始值敏感而且数据分布严重偏斜时的高斯假设很容易失真。对比下来DBSCAN在风电-负荷这种“大样本、含噪声、形状不确定、不想拍脑袋定K”的场景里是最省心、最稳健的路径。2. 参数选择DBSCAN的关键在哪2.1 Eps和MinPts怎么配才算合理DBSCAN参数只有两个——Eps和MinPts但这两个参数调不好聚类结果就是灾难。Eps设得太大所有点都被并进一个簇削减等于没削设得太小簇被拆得七零八落噪声点暴增代表场景失去统计意义。MinPts设得太低任何一小撮点都能聚成一个簇结果里会出现大量细碎类设得过高真正的小簇又会被误杀成噪声。我先说经验规则再讲严格方法。MinPts一般取大于等于数据维度加1也就是MinPts不小于dim1。对于风电-负荷场景单个场景序列如果按24小时采样维度是24MinPts至少取25实际我常用30到50效果比较稳。为什么有dim1这个下限因为密度定义本身需要足够多的点来支撑维度越高维度空间中的点越稀疏需要的邻居数就越多否则“密度”无从谈起。2.2 k距离图法确定EpsEps的确定是DBSCAN使用里最大的一道坎几乎没人能一次估准。我习惯用k距离图法来定初值简单说就是先算每个点到它的第k个最近邻的距离k取MinPts减1然后把这些距离从小到大排序并画成曲线。曲线会在某个位置出现明显的“肘部”拐点这个拐点对应的距离值就是Eps的理想参考值。为什么这个方法是有效的因为DBSCAN的核心假设是簇内点密度高簇间点密度低。所以每个点到第k近邻的距离在簇内的会普遍偏小而在簇边界的点会急剧变大。把距离升序排列后曲线会先平缓后陡峭拐点恰好划分了“簇内距离”和“簇间距离”的过渡带。我以某风电场8760小时数据为例分24维场景后算k30的距离图拐点大约在3.5附近我就把Eps初设定为3.5然后再根据聚类轮廓系数微调到3.8效果不错。2.3 数据预处理对参数的影响不可忽略场景数据如果不做标准化Eps和MinPts几乎没有任何参考价值。风功率的数值范围一般是0到额定功率负荷在几百到几千兆瓦之间两者量纲差好几倍如果直接拼在一起算欧氏距离负荷维度会彻底压制风电维度聚类结果被负荷牵着鼻子走。我通常先把所有特征做Z-score标准化也就是减去均值除以标准差让每个维度都落在相近的尺度上或者做Min-Max归一化到[0,1]区间两种方式我都试过Z-score在存在极端离群点时更稳因为Min-Max会被极端值拉伸正常数据的分布区间。标准化完成之后再画k距离图Eps值就有了跨天、跨季节数据的可比性。否则你上个月调的Eps换到下一个数据集直接就失效了。2.4 距离度量的选择DBSCAN默认用欧氏距离但这里有个值得讨论的细节。风电和负荷场景都是时间序列欧氏距离度量的是“逐时刻差值平方和的根”它把每个时刻点同等对待不考虑相邻时刻的相关性。这个特点既是优点也是局限。如果场景之间的差异主要体现在整体幅值上欧氏距离就够用如果更关注波形形状的相似性可以先对序列做DTW动态时间规整再算距离但DTW的计算开销远高于欧氏距离样本量一大基本跑不动。我的做法是分阶段考量做基础场景削减时用欧氏距离即可因为风电和负荷场景的主要差异本来就在幅值和时间分布上欧氏距离的结果已经足够清晰如果后续要做更精细的时序形态分析可以考虑用PCA降维后再做DBSCAN既保留主要信息又降低维度灾难风险。MATLAB的pca函数很成熟降维到10到15维时聚类速度会明显提升k距离图的拐点也更尖锐。3. MATLAB代码实现与实操步骤3.1 数据准备与场景矩阵构建假设我们有一年8760小时的风电功率和负荷数据采样间隔1小时。要做日场景削减可以把8760小时的数据按每24小时一段切分得到365个场景矩阵。其中每一行是一个场景包含24小时的风电功率和24小时的负荷功率拼在一起就是48个特征维度的矩阵365行乘48列。下面是MATLAB数据准备的核心代码片段% 读取原始风电和负荷时间序列单位: MW wind_power load(wind_data.mat).wind; % 8760 x 1 load_power load(load_data.mat).load; % 8760 x 1 % 按日切分每天 24 点 days 365; n_pts 24; wind_scenes reshape(wind_power(1:days*n_pts), n_pts, days); load_scenes reshape(load_power(1:days*n_pts), n_pts, days); % 按列拼接: 每行包含 24 维风电 24 维负荷 scenes [wind_scenes, load_scenes]; % 365 x 48 % Z-score 标准化标准化参数需保存削减后要反标准化还原场景 [scenes_std, mu, sigma] zscore(scenes);这段代码里有个容易忽略的细节reshape之后一定要加转置符号因为MATLAB的reshape按列填充直接把行向量重排成矩阵时如果不转置场景顺序和时刻顺序全反了。我一开始就栽在这里画出来的场景曲线全是乱的排查了半天才发现是行列方向的问题。3.2 自定义DBSCAN函数实现MATLAB自带的dbscan函数从R2019a开始提供我平时基本用它但为了讲清原理并适配一些复杂需求我写了一份自定义实现核心逻辑如下function [idx, core_flags] dbscan_custom(X, eps, minpts) n size(X, 1); idx zeros(n, 1); % 0 表示未分类-1 表示噪声 core_flags false(n, 1); cluster_id 0; % 预计算距离矩阵数据量大时可改为 pdist2 分块计算 D pdist2(X, X); for i 1:n if idx(i) ~ 0 continue; end % 找邻域点 neighbors find(D(i, :) eps); if numel(neighbors) minpts idx(i) -1; % 先标记为噪声 continue; end cluster_id cluster_id 1; idx(i) cluster_id; queue neighbors; % BFS 扩展核心点邻域 head 1; while head numel(queue) p queue(head); head head 1; if idx(p) -1 idx(p) cluster_id; % 边界点归入当前簇 end if idx(p) ~ 0 continue; end idx(p) cluster_id; p_neighbors find(D(p, :) eps); if numel(p_neighbors) minpts core_flags(p) true; queue [queue, p_neighbors]; %#okAGROW end end end % 未访问到且邻域不足的点为噪声 idx(idx 0) -1; end这段代码虽然是教学版但逻辑是完整的每个核心点出发通过BFS不断吸收邻域内的边界点最终形成一个连通密度簇。实际工程建议直接用MATLAB内置的dbscan函数它底层做了索引加速大样本下比我的循环实现快至少一个数量级。内置函数调用方式很简单[idx, core_flags] dbscan(scenes_std, eps, minpts);3.3 k距离图绘制的实现细节在跑聚类之前先用k距离图来定Eps初值。这里有一个工程细节k值取值应该是MinPts减1因为k距离本身是指“第k近邻居的距离”而DBSCAN的邻域判断标准是包含自身在内共MinPts个点所以第MinPts-1近邻居的距离就是Eps的临界参考。代码如下k minpts - 1; [n, ~] size(scenes_std); % 用 pdist2 算两两距离样本多时用并行计算加速 D pdist2(scenes_std, scenes_std); D_sort sort(D, 2); % 取第 k 近的距离 k_dist D_sort(:, k 1); % 第1列是自身距离0所以偏移一位 k_dist_sorted sort(k_dist, 1); % 画出k距离曲线找拐点 figure; plot(k_dist_sorted, LineWidth, 1.2); xlabel(样本序号按距离排序); ylabel(sprintf(第%d近邻距离, k)); title(k距离曲线——用于确定Eps初值); grid on;需要说明的是D_sort(:, k 1)里的k1是因为第1列是自身距离0第2列才是最近邻距离第k1列就是第k近邻的距离刚接触的人很容易在这里差一行。拐点位置我用肉眼看再结合聚类评价指标微调。有经验之后你会发现拐点区域的Eps值即使偏差10%到20%聚类结果大体还能接受真正要警惕的是Eps取到拐点以下那会导致簇被严重拆分。3.4 聚类结束后生成代表场景得到聚类标签之后最关键的一步就是为每个簇挑选代表场景并计算对应的概率权重。我的做法是最朴素的“簇内中心点法”计算簇内所有场景的平均曲线然后选离这个平均线最近的原始场景作为代表。% 聚类结束idx 为每个场景的簇标签-1 表示噪声 clusters unique(idx(idx 0)); num_clusters numel(clusters); % 存储代表场景和对应概率 rep_scenes zeros(num_clusters, size(scenes, 2)); scene_probs zeros(num_clusters, 1); for c 1:num_clusters members find(idx clusters(c)); % 簇内平均场景标准化空间 mean_scene_std mean(scenes_std(members, :), 1); % 找到离平均场景最近的原始场景 dist_to_mean pdist2(mean_scene_std, scenes_std(members, :)); [~, min_loc] min(dist_to_mean); rep_scene_std scenes_std(members(min_loc), :); % 反标准化还原到原始物理量纲 rep_scenes(c, :) rep_scene_std .* sigma mu; % 概率权重簇内样本数 / 总样本数 scene_probs(c) numel(members) / n; end % 剔除噪声场景重新归一化概率 scene_probs scene_probs / sum(scene_probs);这里有个概率归一化的细节值得讲清楚DBSCAN会标记噪声点如果噪声点的占比是5%那么剩下95%的样本才被聚成有效簇有效簇的概率总和只有0.95。直接拿每个簇的样本数除以总样本数得到的概率加起来不为1后续代入随机规划模型会出问题必须对有效簇概率做一次归一化处理也就是除以总有效簇概率之和。3.5 可视化验证聚类效果削减完一定要可视化肉眼判断永远是最直接的验证手段。我一般画三张图第一张是聚类结果的散点图用前两个主成分做坐标轴不同簇用不同颜色标出来噪声点单独用黑色第二张是原始场景集和代表场景的对比图把所有原始场景画成浅色细线代表场景画成粗黑线叠加进去第三张是削减前后概率分布对比比如画出总出力分布的经验累计分布函数看两者是否吻合。% 用PCA降维到二维空间把聚类结果画出来 [coeff, score] pca(scenes_std); figure; gscatter(score(:,1), score(:,2), idx, bgrcmyk, ., 8); hold on; noise_idx (idx -1); plot(score(noise_idx,1), score(noise_idx,2), kx, MarkerSize, 6); legend(聚类簇, 噪声点); xlabel(主成分1); ylabel(主成分2); title(DBSCAN聚类结果可视化); grid on;gscatter是MATLAB里很好用的分组散点函数颜色会自动按组分配但是要注意它默认的颜色数量有限当聚类数超过7个时颜色会循环我吃过一次亏后来干脆改用scatter加自定义颜色映射这样在簇数多的时候依然能一眼分辨。4. 参数调优与典型问题排查4.1 聚类数量忽多忽少怎么办我在实际应用中最常见的现象是同一套数据Eps稍微改个0.1聚类数就从8跳到15完全没有稳定性。排查下来主要原因是标准化后的数据分布不均匀某些区域的样本密度恰好处于阈值临界点一点点扰动就会让核心点判定翻转。解决思路有两个。第一个思路是做Eps的“扫描式”调参比如在k距离图拐点附近以0.1为步长遍历Eps每次跑完DBSCAN后计算一个评价指标我推荐轮廓系数MATLAB有现成函数能够直接计算。选轮廓系数最高、且噪声占比不超过10%的Eps这样既保证聚类质量又不会因为追求纯密度结构而丢掉太多样本。第二个思路是提升样本量如果原始场景只有几十个Eps的微小变化会直接导致核心点判定不稳定建议把多个年份的数据拼在一起或者用bootstrap重采样扩充样本量让密度估计更平滑。4.2 噪声点过多场景被削得太狠如果聚类结果里噪声点占比超过20%那场景削减后的概率分布会和原始数据明显偏离优化结果也会失真。噪声多的根源通常是Eps太小或者是MinPts取太高导致很多正常样本缺邻居被误判成离群臂。我遇到过一种特殊情况某天的风电场景确实出了极端气象事件比如台风切出时段整条曲线都和其他天差异巨大这种样本被标成噪声其实是合理的场景削减时本来就该把它剔除因为它们对典型运行方式没有代表性。但剔除过多就需要调整参数。我的经验是把Eps适当往拐点右侧偏移10%到20%让边界更宽松同时把MinPts降回dim1的下限。如果调整后噪声占比还在15%以上那就要怀疑是不是标准化方式出了问题可以改试Min-Max归一化有时Z-score对稀疏尾部数据的扩张效果太强会把正常样本推到密度分布的外围。4.3 聚类结果里出现大量“孤立型单点簇”有时候聚类数会突然变得非常大里面有一堆只有两三个样本的小簇这些簇的代表场景几乎和原始样本没区别削减后的场景数量不但没减少反而增加了。这里的核心原因是MinPts取得太低了。MinPts是密度下限取得太小任何一小撮互相接近的点都能形成合法簇。我建议MinPts至少取到20以上而不是机械地遵循dim1的下限。对于48维的风电-负荷联合场景矩阵我最终把MinPts固定为35这个值既保证簇有足够统计量又不会把小簇误杀成噪声。另一个小技巧是聚类完成后过滤掉样本数小于5个的簇把它们的样本并入最近邻簇或者直接标为噪声从工程角度这相当于是“聚类后清洗”可以有效减少碎片簇。4.4 大数据量下计算效率低下DBSCAN内置函数的计算复杂度虽然比传统实现低但用pdist2预计算完整距离矩阵时n个样本就是n×n的存储压力。365个场景还好但如果你把采样粒度拉高到96个时刻点或者场景个数上到几千完整距离矩阵的内存占用就非常可观。比如5000个场景距离矩阵就是5000×5000用double类型存储需要200MB内存这还只是单精度如果再跑多次调参内存和实时都扛不住。我的解决办法是分块计算。先按k距离图估计Eps然后用内置dbscan函数它不会生成完整距离矩阵而是内部按块处理。如果自定义实现可以在BFS过程中用KD树做邻域搜索MATLAB的knnsearch支持KD树加速最坏复杂度能降到对数级别实测5000样本时速度能提升几十倍。如果你的数据实在太大还有一个降维思路先对场景做PCA降维到10维以下再去跑DBSCAN计算量和存储量都能大幅降低聚类的代价是需要验证降维后的信息保留率是否足够。5. 场景削减结果的校验与在后续优化里的应用5.1 削减前后分布一致性校验场景削减做的到底好不好不能光看聚类图还得做定量校验。我一般用概率分布距离指标来对比原始场景集和削减后场景集最常用的是Wasserstein距离在MATLAB里可以用ecdf估计两个经验分布再去计算一阶Wasserstein距离。削减后的Wasserstein距离越小说明代表场景集和原始集越接近。% 计算原始和削减后的总出力经验分布 raw_profile sum(scenes, 2); % 原始365个场景的日总电量 rep_profile sum(rep_scenes .* scene_probs, 1); % 这个计算方法其实是期望值要对比分布需要加权累积 % 正确做法画出两张经验CDF曲线对比 [f_raw, x_raw] ecdf(raw_profile); % rep_profile_series 由 rep_scenes 和 scene_probs 构造出完整的代表场景集合 rep_profile_series repmat(rep_scenes(:,1), 1, num_clusters) .* rep_scenes(:,1); % 示例非严谨写法 figure; plot(x_raw, f_raw, b-, LineWidth, 1.5); hold on; % 生成加权经验分布 [f_rep, x_rep] ecdf(rep_scenes * scene_probs); plot(x_rep, f_rep, r--, LineWidth, 1.5); legend(原始分布, 削减后分布); xlabel(日总出力MW·h); ylabel(经验累计概率); grid on;这段代码里我特意留了一个注释提醒自己rep_scenes * scene_probs这个操作算的是“代表的加权平均”得到的是一个单条期望曲线不是分布。想对比分布正确做法是把每个代表场景按它对应的概率权重复制成多条样本或者直接用加权经验CDF函数。这是我调试时最容易错的地方值得提醒大家注意。5.2 削减结果在随机规划里的应用削减出来的代表场景和概率最直接的应用场景就是随机机组组合和随机经济调度。传统确定性调度只用一组预测曲线随机优化则把代表场景代入目标函数目标变为所有场景下成本概率加权之和。DBSCAN削减相对K-means有一个特别自然的好处它已经返回了噪声点我们在构建随机优化场景树时可以直接把极端噪声场景作为“极端场景”单独加入约束比如用来做备用容量校核而不需要像K-means那样手动从聚类结果里挑极端场景。这种“常规场景极端场景”的组合建模方式比单纯依赖某一种聚类方法做全量削减更贴合工程实际。实际代入时把rep_scenes的每一行作为一个不确定 scenarioscene_probs作为这个scenario的权重写进Yalmip或Matpower的优化框架里就行。这里有个注意事项削减后的场景数量虽然少了但如果DBSCAN聚类出6个簇每个簇还保留了自己的形态差异那么模型中约束的数量会变为原来的6倍计算量的下降幅度取决于聚类数减了多少。这里我建议对比一下削减前后的求解时间如果聚类数太多导致削减没意义可以考虑调大Eps、降低簇数。5.3 DBSCAN与其他削减方法打组合拳最后分享一个我后期摸索出来的策略——DBSCAN并不一定要单独完成全部削减工作。可以先用DBSCAN做一次粗筛把明显离群的极端场景剔掉同时得到一个大致的簇结构然后在每一簇内部再用K-means或ARNOLD方法做精细削减每个簇内削减出2到3个场景。这样既保留了DBSCAN处理噪声和大形状复杂簇的优势又发挥K-means在局部区域形状规整时聚类高效的优点。具体操作时我先跑DBSCAN得到簇标签然后对每个簇内的样本分别调用kmeans函数设置K等于每个簇内样本数对总场景数的比例乘以目标场景总数最后把各簇削减出的子场景合并再重新归一化概率。实测下来这种组合方法在365个日场景的风电-负荷联合数据上能把场景数量从365削减到8到10个同时Wasserstein距离比纯DBSCAN削减低10%左右计算时间却只比纯DBSCAN多了不到半秒。6. 参数调优速查表与我的最终建议说了这么多我把DBSCAN在场景削减中的参数调优经验整理成一张速查表方便各位直接对照参考参数或环节推荐设置关键说明数据标准化Z-score优先标准化后Eps才具备跨数据集可比性最小邻域数MinPtsdim1~3×(dim1)48维场景建议取35~50太小会出碎片簇Eps初始值k距离图拐点处k取MinPts-1拐点在曲线陡升段Eps调整方向轮廓系数噪声占比双指标噪声占比20%时Eps往大调距离度量欧氏距离为主高维可先PCA降维再聚类碎片簇处理样本数5并入近邻簇防止无效场景过多概率归一化有效簇概率重新归一噪声剔除后概率总和不为1必须归一再代入模型结果校验Wasserstein距离经验CDF分布一致性比散点图美观更重要大规模样本内置dbscan函数避免pdist2全矩阵内存压力大最后再说一点个人体会。DBSCAN在场景削减这个任务里最大的价值不在“聚得好不好看”而在于它把“什么是正常场景、什么是异常场景”的判断标准变成了显式参数。K-means让你纠结选几个聚类中心DBSCAN让你考虑的是“一个场景周围有多少邻居才算正常”后者和物理直觉更贴近也更容易解释给项目里的其他成员听。做可再生能源不确定性建模这么多年我越来越觉得一个好的场景削减方法不只是数学工具更像是在帮你理解数据本身的脾气。踩过几次参数坑之后现在我用DBSCAN跑削减第一件事永远是画k距离图哪怕只是瞄一眼曲线形状心里就有底了。