做随机优化的人应该都有过这种体验风电场景刚生成了一大堆两阶段规划还没来得及跑内存先撑不住了。我手头有个含风电接入的机组组合算例初始场景数4000每场景24个时段叠加负荷不确定性后直接求解的耗时从十几分钟飙到接近两小时。后来改用场景削减把4000个场景压到二十几个典型场景计算时间降到原来的十分之一误差却控制在可接受范围内。这篇文章用完整MATLAB代码讲讲我基于DBSCAN密度聚类实现风电-负荷联合场景削减的方法包括聚类参数怎么选、典型场景怎么提取、削减质量怎么评估以及我在实际调试中踩过的一堆坑。1. 项目核心场景削减到底在削什么1.1 新能源随机优化里为什么离不开场景削减风电出力受气象条件影响负荷也有峰谷波动和随机偏差。工程上常用蒙特卡洛抽样生成大量场景来描述这种不确定性每个场景是一条风电功率曲线和一条负荷曲线的组合。场景越多对概率分布的刻画就越精细但随机规划问题的规模也随之急剧膨胀。以两阶段随机机组组合为例每个场景对应一组二阶段决策变量和约束条件。场景数从100增加到2000约束矩阵规模可能翻几十倍单纯靠求解器硬解基本不可行。场景削减就是在这个矛盾里做折中从原始场景集合中挑出一部分代表性场景并重新赋予概率权重使得削减后的离散概率分布在某种距离度量下尽量逼近原始分布同时把场景规模降下来。这个需求在配电网规划、储能容量配置、电力市场出清里都很常见尤其适合需要反复求解随机优化模型的场景。场景削减不是一个新概念但不同方法的削减质量差别很大直接影响优化结果的可靠性。1.2 削减效果怎么量化评估光说“差不多”不行工程上需要可量化的指标来判断场景削减方法是否可用。我一般关注三类指标指标类型具体计算方式判据建议期望曲线偏差削减前后风电/负荷期望曲线的2-范数相对误差控制在5%以内分布尾部偏差削减前后P5/P50/P95分位数的相对误差P95误差控制在10%以内概率一致性削减后场景概率之和是否为1各场景概率非负必须严格满足分位数指标容易被忽略但很重要。风电场次极端出力场景如果被削减掉调度方案对极端工况的适应性就会变差。所以我在评估时不仅看期望值还会专门统计削减前后极端分位数的变化情况。1.3 场景削减在算法链路中的位置场景生成、场景削减、随机优化三者是一整条流水线先由历史数据或预测误差模型生成大量原始场景再通过削减得到少量典型场景及对应概率最后把这些典型场景离散概率输入随机优化模型求解。DBSCAN密度聚类是削减环节的实现工具它解决的是“原始场景该保留哪些、各占多少概率”的问题不涉及后续优化模型的构建这也是这篇代码能独立运行并复用的原因。2. 方案选型为什么用DBSCAN而不是K-means2.1 主流场景削减思路对比学术和工程里最常遇到的场景削减方法大概有三条技术路线基于K-means的场景削减先指定削减后的场景数量K迭代聚类后取簇中心作为典型场景概率由簇内样本数占比决定。实现简单收敛快但K需要提前给定对初始中心敏感且每个场景都被强制归属到某个簇异常场景也会被拉进某个簇参与均值计算。基于层次聚类的后向削减从全场景集合开始每次合并距离最近的两个场景迭代到剩余指定数量为止。优点是无需预设初始中心缺点是计算量大且合并过程不可逆。基于最优概率距离的削减算法理论上在Kantorovich距离意义下做最优削减但实现复杂度高通常也只用于小规模场景。2.2 DBSCAN在场景削减中的差异优势DBSCAN不太一样。它依据样本分布的密集程度自动发现簇核心参数就两个——邻域半径eps和最小样本数MinPts。把场景看作高维空间里的点集DBSCAN能自动识别出“哪些场景经常一起出现”同时把周围稀疏分布的异常场景标记成噪声点。体现在场景削减里就是三个特点不需要预先指定场景削减个数。这在工程里非常实用因为很多时候我真的不知道最后该留多少个场景只知道希望某个数量范围DBSCAN直接告诉我“你这个数据集天然形成了这么多簇”。可以剔除非典型场景。原始场景里经常混入一些极端抽样的离群场景K-means会把它们硬塞进某个簇典型场景就容易被污染。DBSCAN直接把这些离群点标为噪声并在概率重构时剔除掉得到典型场景更纯。能发现任意形状的簇。风电场景不是简单球形分布古德曼分布、多峰场景在空间里的形状不规则K-means对球形簇假设较强DBSCAN没有这个限制。2.3 要清楚DBSCAN的问题再下手DBSCAN不是银弹特别是场景维数较高时。假设每个场景有48个时段特征24时风电24时负荷这个维数下空间距离趋于稀疏密度概念被稀释这时候需要用以下手段缓解先做标准化或PCA降维让每个特征量纲一致、距离集中在有效方向上eps的选择要更谨慎靠k-distance曲线结合实验效果综合判断如果数据量太大pdist2计算距离矩阵是O(N²)级别的需要注意内存和时间开销。这些点我会在第四部分展开讲具体处理方式。3. DBSCAN核心原理与MATLAB调用细节3.1 eps和MinPts的物理意义与选择思路DBSCAN的密度定义很直白对某个点如果它的eps半径邻域内至少有MinPts个点它就算核心点。由核心点的密度可达关系连成簇落在任何簇外且不满足核心条件的点就是噪声。eps决定邻域范围。太小密度高的地方被割裂成很多碎簇甚至大片点被标为噪声太大不同模式的场景被并成一坨削减后的场景太粗糙。MinPts决定核心点门槛。调大MinPts会要求每个簇必须有足够多的样本支撑碎片簇会被吞并调小则容易产生很多小簇概率权重被稀释。设置MinPts可以参考一个粗经验不小于数据特征维数加1。例如48维场景可以设MinPts在10到15之间。实际产业项目里不需要追求理论最优直接在合理区间里跟eps配合着试就行这一点我在第四部分给出一套可操作的标定流程。3.2 场景间距离度量怎么定义才合理DBSCAN聚类前必须先定义场景间“距离”。不同距离度量直接改变聚类结果的结构这里我建议结合场景物理意义做选择欧氏距离最常用计算快适合数值尺度均匀的场景序列。使用前需要标准化否则风电和负荷数量级不同距离会被量纲大的变量主导。加权欧氏距离如果实际业务更看重风电误差或者负荷峰值可以给不同时段、不同变量加权重。DTW动态时间规整适合时间轴不对齐的场景但风电负荷场景都是固定时段序列对齐性良好DTW的收益很小计算量却成倍增加不推荐。夹角余弦距离适合场景形状相似但幅值缩放明显的情况但由于忽略了出力绝对值在容量约束和调度场景中往往不合适。实际项目中我的默认选择是“标准化后的欧氏距离”。这里有个容易踩坑的细节如果直接在原始功率值上算欧氏距离风电和负荷各自波动幅度的量级完全不对等负荷动辄上千兆瓦风电只有几百兆瓦DBSCAN会把距离的差异几乎完全由负荷决定。所以必须先对每个特征做零均值单位方差标准化再进入距离计算。3.3 MATLAB的dbscan函数使用细节MATLAB从R2019a开始提供dbscan函数位于Statistics and Machine Learning Toolbox工具箱里。调用格式是idx dbscan(X, eps, MinPts)X是N行P列的矩阵每行是一个样本即一个场景如果场景同时包含风电和负荷就把两个序列拼接成一行。返回的idx是N行1列向量每个值代表该样本所属的簇编号编号从1开始噪声点统一为-1。额外支持Distance参数可以指定euclidean、squaredeuclidean、mahalanobis甚至函数句柄。也可以传入预先算好的距离矩阵配合Distance,precomputed来用适合在自定义距离度量的场景下复用。需要提醒一点dbscan函数内部虽然对输入做了缓存优化但和所有基于距离密度的方法一样最坏情况计算复杂度是O(N²)场景规模上万时要做好心理准备。4. 完整实践MATLAB实现风电-负荷场景削减全流程这一节是整篇的核心。我把完整流程拆成五个部分从数据生成到削减评估每一段代码都是可以直接拿走的。4.1 第一步生成原始场景数据原始场景的生成方式不唯一可以基于历史出力数据重采样也可以基于预测误差做蒙特卡洛抽样。为了演示完整链路我用一段简化但符合实际统计规律的代码模拟生成2000个风电场景和同数量负荷场景每个场景包含24个时段。风电出力生成用典型日过程叠加AR(1)噪声和随机扰动负荷用双峰日负荷曲线叠加扰动。%% 参数设置 clear; clc; close all; rng(2024); N 2000; % 原始场景数 T 24; % 时段数 Pw_max 300; % 风电场额定容量 MW Pl_max 1200; % 最大负荷 MW %% 1. 生成风电原始场景 t (1:T); % 模拟一个日尺度风电出力趋势夜间大、午间小 wind_profile 0.55 - 0.35*cos(2*pi*t/T) 0.15*sin(2*pi*(t-6)/24); wind_profile wind_profile / max(wind_profile); Pw_raw zeros(N, T); for i 1:N base wind_profile * Pw_max; ar_noise filter(0.8, [1, -0.4], randn(1, T) * 25); % AR(1)扰动模拟时间相关性 random_error randn(1, T) * 15; % 非相关误差 Pw_raw(i, :) max(0, min(Pw_max, base ar_noise random_error)); end %% 2. 生成负荷原始场景 load_profile 0.6 0.25*sin(2*pi*(t-8)/24) 0.15*sin(2*pi*(t-18)/24); load_profile load_profile / max(load_profile); Pl_raw zeros(N, T); for i 1:N base load_profile * Pl_max; Pl_raw(i, :) max(0, min(Pl_max, base randn(1, T) * 45)); end %% 3. 合并为联合场景矩阵 X [Pw_raw, Pl_raw]; % 每行是一个场景前T列风电、后T列负荷AR(1)噪声是很多人在做场景生成时容易漏掉的一环。风电出力在相邻时段有强自相关性如果每个时刻独立加高斯误差生成的场景曲线“毛刺感”很强缺乏物理合理性。加一层AR(1)滤波后场景曲线更贴近真实波动轨迹。4.2 第二步标准化与eps参数确定标准化这一步很关键原因前面说过。这里我不做PCA降维保持48维原始特征因为48维还不算过高但标准化必须做%% 4. 场景标准化 mu mean(X, 1); sigma std(X, 0, 1); sigma(sigma 1e-6) 1; % 防止某些时段零方差导致除零 X_norm (X - mu) ./ sigma; %% 5. 用k-distance曲线辅助选择eps MinPts 10; % 密度阈值先给一个经验初始值 k MinPts - 1; % 计算第k近邻距离 % 计算每个样本到其他样本的第k近邻距离 % 使用pdist2的Smallest参数一次得到全部近邻距离 [~, Dist] pdist2(X_norm, X_norm, euclidean, Smallest, MinPts); % Dist是MinPts行N列矩阵第1行是自己距离0第MinPts行是第MinPts近邻 kdist sort(Dist(MinPts, :), ascend); % 第k近邻距离k MinPts - 1 % 画出k-distance曲线 figure; plot(1:N, kdist, LineWidth, 1.5); grid on; xlabel(场景序号 (按距离排序)); ylabel(第MinPts近邻距离); title(k-distance 曲线辅助选择eps);k-distance曲线的横坐标是样本序号纵坐标是每个样本到第MinPts近邻的距离。曲线从平坦到急剧上升的转折点一般俗称“拐点”或“肘部”对应的距离值可以作为eps的起始参考。实际操作经验是拐点往往不明显尤其高维场景。我会把排序后第10%到30%位置的kdist值作为候选eps区间然后取中间值初跑一次DBSCAN根据聚类结果再做微调。4.3 第三步运行DBSCAN聚类并检查结果拿到eps候选值后直接调用dbscan函数完成聚类同时做一次检查性统计%% 6. DBSCAN聚类 eps_value 2.8; % 根据k-distance曲线拐点附近取值后续可调 [idx, corepts] dbscan(X_norm, eps_value, MinPts); num_clusters max(idx); num_noise sum(idx -1); fprintf(聚类结果: 簇个数%d, 噪声点个数%d (%.2f%%)\n, ... num_clusters, num_noise, num_noise/N*100); % 看一眼每个簇的样本规模避免出现极端小簇 cluster_sizes accumarray(idx(idx 0), 1); fprintf(各簇样本数: %s\n, mat2str(cluster_sizes(:)));这里结果评估要看几个方面。簇个数是否落在预期的削减规模区间内比如我希望后续随机优化保留10到30个场景簇个数在这个区间就比较理想。噪声点占比不应太高超过20%说明eps偏小需要适当增大。极端小簇比如只有一两个样本过多说明密度分裂严重也需要调整参数。注意MATLAB的dbscan函数里MinPts参数可以直接传整数也可以传为一个比例值比如0.05表示样本总数的5%。我一般用固定整数因为固定值在场景削减里语义更清楚一个典型场景至少需要多少个原始场景支撑。4.4 第四步典型场景提取与概率重估聚类完成后每个簇就对应一个典型场景。这里有两种提取典型场景的方式取簇内样本均值作为典型场景曲线平滑但可能把极端出力平滑掉取簇内离重心最近的原始样本作为典型场景保留真实物理模式可追溯性强。我强烈推荐第二种也就是在簇内找距离质心最近的样本作为典型场景。这样生成的典型场景是“真实存在的场景曲线”和上游场景生成环节的数据保持一致性后续调度方案也能追溯到原始数据。概率的计算逻辑是每个簇的样本数除以总样本数扣除噪声点后并按比例归一化。%% 7. 提取典型场景并重新赋概率 num_typical num_clusters; % 典型场景数等于簇个数 rep_scen zeros(num_typical, 2*T); prob zeros(num_typical, 1); valid_count N - num_noise; % 剔除噪声后的有效样本数 for c 1:num_clusters members find(idx c); centroid mean(X(members, :), 1); % 簇内质心原始值空间 d_inner pdist2(centroid, X(members, :), euclidean); [~, pos] min(d_inner, [], 2); % 离质心最近的原始场景 rep_scen(c, :) X(members(pos), :); % 保留原始场景作为典型场景 prob(c) length(members) / valid_count; % 概率重估 end % 归一化确保概率和为1 prob prob / sum(prob); fprintf(削减后典型场景个数: %d\n, num_typical); fprintf(各场景概率: %s\n, mat2str(prob(:), 3));这里有个细节概率归一化时用valid_count而不是N。原因是把噪声场景看作被削减掉的不重要场景它的概率应当分摊到保留场景上而不是让概率和小于1。如果直接用N做分母削减后场景概率之和不为1后面做随机规划等式约束就不守恒。4.5 第五步削减效果可视化评估削减效果不能只靠眼睛看我用定量指标加可视化两条线并行。可视化方面我画四个子图原始场景抽样曲线、削减后的典型场景曲线、风电期望曲线对比、负荷期望曲线对比。%% 8. 削减效果评估 % 削减前后风电/负荷期望曲线 wind_orig_mean mean(Pw_raw, 1); load_orig_mean mean(Pl_raw, 1); % 削减后典型场景按概率加权求期望 wind_red_mean rep_scen(:, 1:T) * prob; load_red_mean rep_scen(:, T1:end) * prob; % 期望曲线相对误差 err_wind norm(wind_red_mean - wind_orig_mean, 2) / norm(wind_orig_mean, 2) * 100; err_load norm(load_red_mean - load_orig_mean, 2) / norm(load_orig_mean, 2) * 100; fprintf(风电期望曲线相对误差: %.2f%%\n, err_wind); fprintf(负荷期望曲线相对误差: %.2f%%\n, err_load); % 分位数误差 q_wind_orig quantile(Pw_raw(:), [0.05 0.5 0.95]); wind_red_all rep_scen(:, 1:T); q_wind_red quantile(wind_red_all(:), [0.05 0.5 0.95]); fprintf(风电P5/P50/P95原始: %s\n, mat2str(q_wind_orig, 3)); fprintf(风电P5/P50/P95削减: %s\n, mat2str(q_wind_red, 3));注意分位数这块我没有用概率加权而是直接统计所有典型场景的值。如果要更严格应该用经验分布函数加权计算分位数。但对于粗略评估直接用典型场景集合统计也够用了。可视化部分我用下面这段代码%% 9. 绘图 figure(Position, [100 100 1000 650]); % 子图1原始场景抽样 subplot(2,2,1); plot(1:T, Pw_raw(1:100, :), Color, [0.7 0.7 0.7]); hold on; plot(1:T, rep_scen(1:min(end,8), 1:T), LineWidth, 1.2); xlabel(时段/h); ylabel(风电功率/MW); title(原始风电场景(灰)与典型场景(彩)); grid on; % 子图2负荷场景对比 subplot(2,2,2); plot(1:T, Pl_raw(1:100, :), Color, [0.7 0.7 0.7]); hold on; plot(1:T, rep_scen(1:min(end,8), T1:end), LineWidth, 1.2); xlabel(时段/h); ylabel(负荷/MW); title(原始负荷场景(灰)与典型场景(彩)); grid on; % 子图3风电期望曲线对比 subplot(2,2,3); plot(1:T, wind_orig_mean, k-, LineWidth, 1.8); hold on; plot(1:T, wind_red_mean, r--, LineWidth, 1.8); xlabel(时段/h); ylabel(期望风电/MW); title(sprintf(风电期望曲线对比 误差%.2f%%, err_wind)); legend(原始, 削减后, Location, best); grid on; % 子图4负荷期望曲线对比 subplot(2,2,4); plot(1:T, load_orig_mean, k-, LineWidth, 1.8); hold on; plot(1:T, load_red_mean, r--, LineWidth, 1.8); xlabel(时段/h); ylabel(期望负荷/MW); title(sprintf(负荷期望曲线对比 误差%.2f%%, err_load)); legend(原始, 削减后, Location, best); grid on;我实际跑下来标准化的48维场景数据在MinPts10、eps取2.8左右时2000个场景通常能削减出20到30个簇噪声占比在5%上下风电和负荷期望曲线误差都能控制在3%以内。这个规模和精度对随机优化来说是相当理想的输入。5. 实战中反复踩坑的经验总结5.1 eps参数定不准场景个数不受控DBSCAN不那么依赖人为指定簇个数但eps和MinPts组合其实间接决定削减规模。最常见的问题是期望保留15个场景结果聚类出来40个簇或者只有3个簇。我处理这个问题的方式是“双参数联动扫描”固定MinPts10把eps从1.0到5.0按0.2步长扫一遍记录每个eps下的簇个数和噪声占比画一条曲线。设备运行几秒钟就能跑完。然后根据“簇个数落在目标区间”的原则反选eps。更保守的做法是加一个“粗聚类转细聚类”的兜底逻辑如果DBSCAN结果簇个数少于预期下限把簇中心提出来再跑一次K-means分到目标数量。这本质上是两阶段聚类一般应对工业场景里的业务硬约束时才会用到。5.2 聚类结果出现大面积噪声点噪声点占比超过20%通常不是“异常场景真的很多”而是eps太小密度可达范围不足很多原本应该聚类在一起的场景被孤立成离群点。我的排查顺序是先看k-distance曲线上第10%位置的值如果明显小于当前eps说明给的eps过小接着画一个二维投影比如用t-SNE或者PCA前两主成分看数据分布确认场景是否真的存在密度分离结构最后尝试把eps放大30%左右再看聚类结果。有些数据集本身没有清晰的密度分离结构这时候DBSCAN反而不如K-means。你可以先用PCA投影快速判断一下数据形态。5.3 削减后典型场景代表性不足即使期望曲线误差很小偶尔也会出现某几个典型场景出现概率极小的现象导致后续随机优化的某些场景权重过低失去意义。这类小概率簇通常对应数据中的小众运行工况不是纯噪声。我建议在概率重估时设置一个概率下限阈值比如1%。低于阈值的簇可以选择强制合并进最近的簇重新计算典型场景和概率或者直接标记为噪声并从样本空间中剔除。具体阈值根据你后级优化对最小场景概率的敏感性确定。5.4 原始数据量大导致聚类太慢DBSCAN的低效根源在于距离矩阵计算。2000个场景时pdist2计算约400万对距离很快上万场景就明显卡顿。我先说一个在实际项目里很好用的提速技巧用DBSCAN的Distance,precomputed模式把距离矩阵一次性算好后传给dbscan。乍一听更慢但实际上如果你反复调参跑同一个数据集预计算距离矩阵可以反复复用节省大量重复计算时间。更大的数据量可以做两步削减先用随机抽样或K-means快速把场景压到初始规模比如2万压到5000再用DBSCAN做精细化密度聚类。两步法在工程上非常实用能兼顾计算速度和削减质量。5.5 关于工具链的补充实现DBSCAN除了MATLAB自带的dbscan函数也可以用Python的sklearn.cluster.DBSCAN逻辑一致。MATLAB方案的好处是在电力系统仿真和优化求解器集成上更方便数据流不用跨语言对接。如果你的团队以Python为主我个人建议直接统一用sklearn版本参数定义更细社区资料也多。6. 一点扩展思路场景削减之后如果发现典型场景数量还是不够用或者希望进一步提高削减质量可以考虑对DBSCAN聚类结果做三层处理第一层用DBSCAN剔除噪声并保留密度核心第二层对核心点做加权K-means获得目标数量的场景第三层用概率最优传输算法把场景概率微调至与原始分布更为匹配。这套组合在需要高精度削减的场景比如跨省送电计划随机评估里很有效。另外风电与负荷的联合场景中如果风电和负荷之间存在相关性建议在标准化之后、聚类之前加入一个旋转操作把联合分布映射到主成分空间。这不仅缓解高维距离稀疏问题还能部分保留变量间相关结构的信息。DBSCAN clustering之后回到原始功率空间提取典型场景即可。我在实际项目中感受到DBSCAN在场景削减里最大的价值不是“比K-means准确多少”而是“不需要提前指定场景个数”和“自动识别异常场景”这两点让削减结果更符合数据的真实结构。和K-means做对比同一批风电负荷数据K-means在K15时某些簇中心明显偏向离群场景DBSCAN聚类后噪声点被剥离开典型场景曲线明显更贴近主流出力区间。最后分享一个用得上的细节无论用什么聚类方法保存典型场景时最好顺便把簇内样本索引也存下来。后面做灵敏度分析、追溯某个典型场景对应的原始工况时这份索引能帮你省大量时间。我一开始没保存后来需要从典型场景反查原始数据时差点把原始场景重新生成一遍。这套代码实测下来在2000个场景、48维特征的数据规模下从聚类到完成评估总耗时不超一分钟完全满足工程迭代需要。建议你拿到代码后先不改动参数完整跑一遍再根据自己数据的实际情况去调eps和MinPts。 基于DBSCAN密度聚类的风电-负荷联合场景削减MATLAB实现
做随机优化的人应该都有过这种体验风电场景刚生成了一大堆两阶段规划还没来得及跑内存先撑不住了。我手头有个含风电接入的机组组合算例初始场景数4000每场景24个时段叠加负荷不确定性后直接求解的耗时从十几分钟飙到接近两小时。后来改用场景削减把4000个场景压到二十几个典型场景计算时间降到原来的十分之一误差却控制在可接受范围内。这篇文章用完整MATLAB代码讲讲我基于DBSCAN密度聚类实现风电-负荷联合场景削减的方法包括聚类参数怎么选、典型场景怎么提取、削减质量怎么评估以及我在实际调试中踩过的一堆坑。1. 项目核心场景削减到底在削什么1.1 新能源随机优化里为什么离不开场景削减风电出力受气象条件影响负荷也有峰谷波动和随机偏差。工程上常用蒙特卡洛抽样生成大量场景来描述这种不确定性每个场景是一条风电功率曲线和一条负荷曲线的组合。场景越多对概率分布的刻画就越精细但随机规划问题的规模也随之急剧膨胀。以两阶段随机机组组合为例每个场景对应一组二阶段决策变量和约束条件。场景数从100增加到2000约束矩阵规模可能翻几十倍单纯靠求解器硬解基本不可行。场景削减就是在这个矛盾里做折中从原始场景集合中挑出一部分代表性场景并重新赋予概率权重使得削减后的离散概率分布在某种距离度量下尽量逼近原始分布同时把场景规模降下来。这个需求在配电网规划、储能容量配置、电力市场出清里都很常见尤其适合需要反复求解随机优化模型的场景。场景削减不是一个新概念但不同方法的削减质量差别很大直接影响优化结果的可靠性。1.2 削减效果怎么量化评估光说“差不多”不行工程上需要可量化的指标来判断场景削减方法是否可用。我一般关注三类指标指标类型具体计算方式判据建议期望曲线偏差削减前后风电/负荷期望曲线的2-范数相对误差控制在5%以内分布尾部偏差削减前后P5/P50/P95分位数的相对误差P95误差控制在10%以内概率一致性削减后场景概率之和是否为1各场景概率非负必须严格满足分位数指标容易被忽略但很重要。风电场次极端出力场景如果被削减掉调度方案对极端工况的适应性就会变差。所以我在评估时不仅看期望值还会专门统计削减前后极端分位数的变化情况。1.3 场景削减在算法链路中的位置场景生成、场景削减、随机优化三者是一整条流水线先由历史数据或预测误差模型生成大量原始场景再通过削减得到少量典型场景及对应概率最后把这些典型场景离散概率输入随机优化模型求解。DBSCAN密度聚类是削减环节的实现工具它解决的是“原始场景该保留哪些、各占多少概率”的问题不涉及后续优化模型的构建这也是这篇代码能独立运行并复用的原因。2. 方案选型为什么用DBSCAN而不是K-means2.1 主流场景削减思路对比学术和工程里最常遇到的场景削减方法大概有三条技术路线基于K-means的场景削减先指定削减后的场景数量K迭代聚类后取簇中心作为典型场景概率由簇内样本数占比决定。实现简单收敛快但K需要提前给定对初始中心敏感且每个场景都被强制归属到某个簇异常场景也会被拉进某个簇参与均值计算。基于层次聚类的后向削减从全场景集合开始每次合并距离最近的两个场景迭代到剩余指定数量为止。优点是无需预设初始中心缺点是计算量大且合并过程不可逆。基于最优概率距离的削减算法理论上在Kantorovich距离意义下做最优削减但实现复杂度高通常也只用于小规模场景。2.2 DBSCAN在场景削减中的差异优势DBSCAN不太一样。它依据样本分布的密集程度自动发现簇核心参数就两个——邻域半径eps和最小样本数MinPts。把场景看作高维空间里的点集DBSCAN能自动识别出“哪些场景经常一起出现”同时把周围稀疏分布的异常场景标记成噪声点。体现在场景削减里就是三个特点不需要预先指定场景削减个数。这在工程里非常实用因为很多时候我真的不知道最后该留多少个场景只知道希望某个数量范围DBSCAN直接告诉我“你这个数据集天然形成了这么多簇”。可以剔除非典型场景。原始场景里经常混入一些极端抽样的离群场景K-means会把它们硬塞进某个簇典型场景就容易被污染。DBSCAN直接把这些离群点标为噪声并在概率重构时剔除掉得到典型场景更纯。能发现任意形状的簇。风电场景不是简单球形分布古德曼分布、多峰场景在空间里的形状不规则K-means对球形簇假设较强DBSCAN没有这个限制。2.3 要清楚DBSCAN的问题再下手DBSCAN不是银弹特别是场景维数较高时。假设每个场景有48个时段特征24时风电24时负荷这个维数下空间距离趋于稀疏密度概念被稀释这时候需要用以下手段缓解先做标准化或PCA降维让每个特征量纲一致、距离集中在有效方向上eps的选择要更谨慎靠k-distance曲线结合实验效果综合判断如果数据量太大pdist2计算距离矩阵是O(N²)级别的需要注意内存和时间开销。这些点我会在第四部分展开讲具体处理方式。3. DBSCAN核心原理与MATLAB调用细节3.1 eps和MinPts的物理意义与选择思路DBSCAN的密度定义很直白对某个点如果它的eps半径邻域内至少有MinPts个点它就算核心点。由核心点的密度可达关系连成簇落在任何簇外且不满足核心条件的点就是噪声。eps决定邻域范围。太小密度高的地方被割裂成很多碎簇甚至大片点被标为噪声太大不同模式的场景被并成一坨削减后的场景太粗糙。MinPts决定核心点门槛。调大MinPts会要求每个簇必须有足够多的样本支撑碎片簇会被吞并调小则容易产生很多小簇概率权重被稀释。设置MinPts可以参考一个粗经验不小于数据特征维数加1。例如48维场景可以设MinPts在10到15之间。实际产业项目里不需要追求理论最优直接在合理区间里跟eps配合着试就行这一点我在第四部分给出一套可操作的标定流程。3.2 场景间距离度量怎么定义才合理DBSCAN聚类前必须先定义场景间“距离”。不同距离度量直接改变聚类结果的结构这里我建议结合场景物理意义做选择欧氏距离最常用计算快适合数值尺度均匀的场景序列。使用前需要标准化否则风电和负荷数量级不同距离会被量纲大的变量主导。加权欧氏距离如果实际业务更看重风电误差或者负荷峰值可以给不同时段、不同变量加权重。DTW动态时间规整适合时间轴不对齐的场景但风电负荷场景都是固定时段序列对齐性良好DTW的收益很小计算量却成倍增加不推荐。夹角余弦距离适合场景形状相似但幅值缩放明显的情况但由于忽略了出力绝对值在容量约束和调度场景中往往不合适。实际项目中我的默认选择是“标准化后的欧氏距离”。这里有个容易踩坑的细节如果直接在原始功率值上算欧氏距离风电和负荷各自波动幅度的量级完全不对等负荷动辄上千兆瓦风电只有几百兆瓦DBSCAN会把距离的差异几乎完全由负荷决定。所以必须先对每个特征做零均值单位方差标准化再进入距离计算。3.3 MATLAB的dbscan函数使用细节MATLAB从R2019a开始提供dbscan函数位于Statistics and Machine Learning Toolbox工具箱里。调用格式是idx dbscan(X, eps, MinPts)X是N行P列的矩阵每行是一个样本即一个场景如果场景同时包含风电和负荷就把两个序列拼接成一行。返回的idx是N行1列向量每个值代表该样本所属的簇编号编号从1开始噪声点统一为-1。额外支持Distance参数可以指定euclidean、squaredeuclidean、mahalanobis甚至函数句柄。也可以传入预先算好的距离矩阵配合Distance,precomputed来用适合在自定义距离度量的场景下复用。需要提醒一点dbscan函数内部虽然对输入做了缓存优化但和所有基于距离密度的方法一样最坏情况计算复杂度是O(N²)场景规模上万时要做好心理准备。4. 完整实践MATLAB实现风电-负荷场景削减全流程这一节是整篇的核心。我把完整流程拆成五个部分从数据生成到削减评估每一段代码都是可以直接拿走的。4.1 第一步生成原始场景数据原始场景的生成方式不唯一可以基于历史出力数据重采样也可以基于预测误差做蒙特卡洛抽样。为了演示完整链路我用一段简化但符合实际统计规律的代码模拟生成2000个风电场景和同数量负荷场景每个场景包含24个时段。风电出力生成用典型日过程叠加AR(1)噪声和随机扰动负荷用双峰日负荷曲线叠加扰动。%% 参数设置 clear; clc; close all; rng(2024); N 2000; % 原始场景数 T 24; % 时段数 Pw_max 300; % 风电场额定容量 MW Pl_max 1200; % 最大负荷 MW %% 1. 生成风电原始场景 t (1:T); % 模拟一个日尺度风电出力趋势夜间大、午间小 wind_profile 0.55 - 0.35*cos(2*pi*t/T) 0.15*sin(2*pi*(t-6)/24); wind_profile wind_profile / max(wind_profile); Pw_raw zeros(N, T); for i 1:N base wind_profile * Pw_max; ar_noise filter(0.8, [1, -0.4], randn(1, T) * 25); % AR(1)扰动模拟时间相关性 random_error randn(1, T) * 15; % 非相关误差 Pw_raw(i, :) max(0, min(Pw_max, base ar_noise random_error)); end %% 2. 生成负荷原始场景 load_profile 0.6 0.25*sin(2*pi*(t-8)/24) 0.15*sin(2*pi*(t-18)/24); load_profile load_profile / max(load_profile); Pl_raw zeros(N, T); for i 1:N base load_profile * Pl_max; Pl_raw(i, :) max(0, min(Pl_max, base randn(1, T) * 45)); end %% 3. 合并为联合场景矩阵 X [Pw_raw, Pl_raw]; % 每行是一个场景前T列风电、后T列负荷AR(1)噪声是很多人在做场景生成时容易漏掉的一环。风电出力在相邻时段有强自相关性如果每个时刻独立加高斯误差生成的场景曲线“毛刺感”很强缺乏物理合理性。加一层AR(1)滤波后场景曲线更贴近真实波动轨迹。4.2 第二步标准化与eps参数确定标准化这一步很关键原因前面说过。这里我不做PCA降维保持48维原始特征因为48维还不算过高但标准化必须做%% 4. 场景标准化 mu mean(X, 1); sigma std(X, 0, 1); sigma(sigma 1e-6) 1; % 防止某些时段零方差导致除零 X_norm (X - mu) ./ sigma; %% 5. 用k-distance曲线辅助选择eps MinPts 10; % 密度阈值先给一个经验初始值 k MinPts - 1; % 计算第k近邻距离 % 计算每个样本到其他样本的第k近邻距离 % 使用pdist2的Smallest参数一次得到全部近邻距离 [~, Dist] pdist2(X_norm, X_norm, euclidean, Smallest, MinPts); % Dist是MinPts行N列矩阵第1行是自己距离0第MinPts行是第MinPts近邻 kdist sort(Dist(MinPts, :), ascend); % 第k近邻距离k MinPts - 1 % 画出k-distance曲线 figure; plot(1:N, kdist, LineWidth, 1.5); grid on; xlabel(场景序号 (按距离排序)); ylabel(第MinPts近邻距离); title(k-distance 曲线辅助选择eps);k-distance曲线的横坐标是样本序号纵坐标是每个样本到第MinPts近邻的距离。曲线从平坦到急剧上升的转折点一般俗称“拐点”或“肘部”对应的距离值可以作为eps的起始参考。实际操作经验是拐点往往不明显尤其高维场景。我会把排序后第10%到30%位置的kdist值作为候选eps区间然后取中间值初跑一次DBSCAN根据聚类结果再做微调。4.3 第三步运行DBSCAN聚类并检查结果拿到eps候选值后直接调用dbscan函数完成聚类同时做一次检查性统计%% 6. DBSCAN聚类 eps_value 2.8; % 根据k-distance曲线拐点附近取值后续可调 [idx, corepts] dbscan(X_norm, eps_value, MinPts); num_clusters max(idx); num_noise sum(idx -1); fprintf(聚类结果: 簇个数%d, 噪声点个数%d (%.2f%%)\n, ... num_clusters, num_noise, num_noise/N*100); % 看一眼每个簇的样本规模避免出现极端小簇 cluster_sizes accumarray(idx(idx 0), 1); fprintf(各簇样本数: %s\n, mat2str(cluster_sizes(:)));这里结果评估要看几个方面。簇个数是否落在预期的削减规模区间内比如我希望后续随机优化保留10到30个场景簇个数在这个区间就比较理想。噪声点占比不应太高超过20%说明eps偏小需要适当增大。极端小簇比如只有一两个样本过多说明密度分裂严重也需要调整参数。注意MATLAB的dbscan函数里MinPts参数可以直接传整数也可以传为一个比例值比如0.05表示样本总数的5%。我一般用固定整数因为固定值在场景削减里语义更清楚一个典型场景至少需要多少个原始场景支撑。4.4 第四步典型场景提取与概率重估聚类完成后每个簇就对应一个典型场景。这里有两种提取典型场景的方式取簇内样本均值作为典型场景曲线平滑但可能把极端出力平滑掉取簇内离重心最近的原始样本作为典型场景保留真实物理模式可追溯性强。我强烈推荐第二种也就是在簇内找距离质心最近的样本作为典型场景。这样生成的典型场景是“真实存在的场景曲线”和上游场景生成环节的数据保持一致性后续调度方案也能追溯到原始数据。概率的计算逻辑是每个簇的样本数除以总样本数扣除噪声点后并按比例归一化。%% 7. 提取典型场景并重新赋概率 num_typical num_clusters; % 典型场景数等于簇个数 rep_scen zeros(num_typical, 2*T); prob zeros(num_typical, 1); valid_count N - num_noise; % 剔除噪声后的有效样本数 for c 1:num_clusters members find(idx c); centroid mean(X(members, :), 1); % 簇内质心原始值空间 d_inner pdist2(centroid, X(members, :), euclidean); [~, pos] min(d_inner, [], 2); % 离质心最近的原始场景 rep_scen(c, :) X(members(pos), :); % 保留原始场景作为典型场景 prob(c) length(members) / valid_count; % 概率重估 end % 归一化确保概率和为1 prob prob / sum(prob); fprintf(削减后典型场景个数: %d\n, num_typical); fprintf(各场景概率: %s\n, mat2str(prob(:), 3));这里有个细节概率归一化时用valid_count而不是N。原因是把噪声场景看作被削减掉的不重要场景它的概率应当分摊到保留场景上而不是让概率和小于1。如果直接用N做分母削减后场景概率之和不为1后面做随机规划等式约束就不守恒。4.5 第五步削减效果可视化评估削减效果不能只靠眼睛看我用定量指标加可视化两条线并行。可视化方面我画四个子图原始场景抽样曲线、削减后的典型场景曲线、风电期望曲线对比、负荷期望曲线对比。%% 8. 削减效果评估 % 削减前后风电/负荷期望曲线 wind_orig_mean mean(Pw_raw, 1); load_orig_mean mean(Pl_raw, 1); % 削减后典型场景按概率加权求期望 wind_red_mean rep_scen(:, 1:T) * prob; load_red_mean rep_scen(:, T1:end) * prob; % 期望曲线相对误差 err_wind norm(wind_red_mean - wind_orig_mean, 2) / norm(wind_orig_mean, 2) * 100; err_load norm(load_red_mean - load_orig_mean, 2) / norm(load_orig_mean, 2) * 100; fprintf(风电期望曲线相对误差: %.2f%%\n, err_wind); fprintf(负荷期望曲线相对误差: %.2f%%\n, err_load); % 分位数误差 q_wind_orig quantile(Pw_raw(:), [0.05 0.5 0.95]); wind_red_all rep_scen(:, 1:T); q_wind_red quantile(wind_red_all(:), [0.05 0.5 0.95]); fprintf(风电P5/P50/P95原始: %s\n, mat2str(q_wind_orig, 3)); fprintf(风电P5/P50/P95削减: %s\n, mat2str(q_wind_red, 3));注意分位数这块我没有用概率加权而是直接统计所有典型场景的值。如果要更严格应该用经验分布函数加权计算分位数。但对于粗略评估直接用典型场景集合统计也够用了。可视化部分我用下面这段代码%% 9. 绘图 figure(Position, [100 100 1000 650]); % 子图1原始场景抽样 subplot(2,2,1); plot(1:T, Pw_raw(1:100, :), Color, [0.7 0.7 0.7]); hold on; plot(1:T, rep_scen(1:min(end,8), 1:T), LineWidth, 1.2); xlabel(时段/h); ylabel(风电功率/MW); title(原始风电场景(灰)与典型场景(彩)); grid on; % 子图2负荷场景对比 subplot(2,2,2); plot(1:T, Pl_raw(1:100, :), Color, [0.7 0.7 0.7]); hold on; plot(1:T, rep_scen(1:min(end,8), T1:end), LineWidth, 1.2); xlabel(时段/h); ylabel(负荷/MW); title(原始负荷场景(灰)与典型场景(彩)); grid on; % 子图3风电期望曲线对比 subplot(2,2,3); plot(1:T, wind_orig_mean, k-, LineWidth, 1.8); hold on; plot(1:T, wind_red_mean, r--, LineWidth, 1.8); xlabel(时段/h); ylabel(期望风电/MW); title(sprintf(风电期望曲线对比 误差%.2f%%, err_wind)); legend(原始, 削减后, Location, best); grid on; % 子图4负荷期望曲线对比 subplot(2,2,4); plot(1:T, load_orig_mean, k-, LineWidth, 1.8); hold on; plot(1:T, load_red_mean, r--, LineWidth, 1.8); xlabel(时段/h); ylabel(期望负荷/MW); title(sprintf(负荷期望曲线对比 误差%.2f%%, err_load)); legend(原始, 削减后, Location, best); grid on;我实际跑下来标准化的48维场景数据在MinPts10、eps取2.8左右时2000个场景通常能削减出20到30个簇噪声占比在5%上下风电和负荷期望曲线误差都能控制在3%以内。这个规模和精度对随机优化来说是相当理想的输入。5. 实战中反复踩坑的经验总结5.1 eps参数定不准场景个数不受控DBSCAN不那么依赖人为指定簇个数但eps和MinPts组合其实间接决定削减规模。最常见的问题是期望保留15个场景结果聚类出来40个簇或者只有3个簇。我处理这个问题的方式是“双参数联动扫描”固定MinPts10把eps从1.0到5.0按0.2步长扫一遍记录每个eps下的簇个数和噪声占比画一条曲线。设备运行几秒钟就能跑完。然后根据“簇个数落在目标区间”的原则反选eps。更保守的做法是加一个“粗聚类转细聚类”的兜底逻辑如果DBSCAN结果簇个数少于预期下限把簇中心提出来再跑一次K-means分到目标数量。这本质上是两阶段聚类一般应对工业场景里的业务硬约束时才会用到。5.2 聚类结果出现大面积噪声点噪声点占比超过20%通常不是“异常场景真的很多”而是eps太小密度可达范围不足很多原本应该聚类在一起的场景被孤立成离群点。我的排查顺序是先看k-distance曲线上第10%位置的值如果明显小于当前eps说明给的eps过小接着画一个二维投影比如用t-SNE或者PCA前两主成分看数据分布确认场景是否真的存在密度分离结构最后尝试把eps放大30%左右再看聚类结果。有些数据集本身没有清晰的密度分离结构这时候DBSCAN反而不如K-means。你可以先用PCA投影快速判断一下数据形态。5.3 削减后典型场景代表性不足即使期望曲线误差很小偶尔也会出现某几个典型场景出现概率极小的现象导致后续随机优化的某些场景权重过低失去意义。这类小概率簇通常对应数据中的小众运行工况不是纯噪声。我建议在概率重估时设置一个概率下限阈值比如1%。低于阈值的簇可以选择强制合并进最近的簇重新计算典型场景和概率或者直接标记为噪声并从样本空间中剔除。具体阈值根据你后级优化对最小场景概率的敏感性确定。5.4 原始数据量大导致聚类太慢DBSCAN的低效根源在于距离矩阵计算。2000个场景时pdist2计算约400万对距离很快上万场景就明显卡顿。我先说一个在实际项目里很好用的提速技巧用DBSCAN的Distance,precomputed模式把距离矩阵一次性算好后传给dbscan。乍一听更慢但实际上如果你反复调参跑同一个数据集预计算距离矩阵可以反复复用节省大量重复计算时间。更大的数据量可以做两步削减先用随机抽样或K-means快速把场景压到初始规模比如2万压到5000再用DBSCAN做精细化密度聚类。两步法在工程上非常实用能兼顾计算速度和削减质量。5.5 关于工具链的补充实现DBSCAN除了MATLAB自带的dbscan函数也可以用Python的sklearn.cluster.DBSCAN逻辑一致。MATLAB方案的好处是在电力系统仿真和优化求解器集成上更方便数据流不用跨语言对接。如果你的团队以Python为主我个人建议直接统一用sklearn版本参数定义更细社区资料也多。6. 一点扩展思路场景削减之后如果发现典型场景数量还是不够用或者希望进一步提高削减质量可以考虑对DBSCAN聚类结果做三层处理第一层用DBSCAN剔除噪声并保留密度核心第二层对核心点做加权K-means获得目标数量的场景第三层用概率最优传输算法把场景概率微调至与原始分布更为匹配。这套组合在需要高精度削减的场景比如跨省送电计划随机评估里很有效。另外风电与负荷的联合场景中如果风电和负荷之间存在相关性建议在标准化之后、聚类之前加入一个旋转操作把联合分布映射到主成分空间。这不仅缓解高维距离稀疏问题还能部分保留变量间相关结构的信息。DBSCAN clustering之后回到原始功率空间提取典型场景即可。我在实际项目中感受到DBSCAN在场景削减里最大的价值不是“比K-means准确多少”而是“不需要提前指定场景个数”和“自动识别异常场景”这两点让削减结果更符合数据的真实结构。和K-means做对比同一批风电负荷数据K-means在K15时某些簇中心明显偏向离群场景DBSCAN聚类后噪声点被剥离开典型场景曲线明显更贴近主流出力区间。最后分享一个用得上的细节无论用什么聚类方法保存典型场景时最好顺便把簇内样本索引也存下来。后面做灵敏度分析、追溯某个典型场景对应的原始工况时这份索引能帮你省大量时间。我一开始没保存后来需要从典型场景反查原始数据时差点把原始场景重新生成一遍。这套代码实测下来在2000个场景、48维特征的数据规模下从聚类到完成评估总耗时不超一分钟完全满足工程迭代需要。建议你拿到代码后先不改动参数完整跑一遍再根据自己数据的实际情况去调eps和MinPts。