搞风电、负荷预测或者电力系统规划的朋友应该都体会过场景分析的份量。风电出力的随机性、负荷波动的时变性叠加起来就是一大把不确定的东西蒙特卡洛一抽就是几百上千条曲线。模型本身不挑数据但求解器扛不住尤其做随机优化、可靠性评估或者输电网规划的时候场景一多整数变量组合直接爆炸原来几分钟能算完的问题拖到一晚上都解不动。所以场景缩减这件事本质就是在计算代价和不确定性保真度之间找平衡。我这次用的方案是基于DBSCAN密度聚类做风电-负荷确定性场景缩减全程在Matlab里实现。相比K-means这种需要预指定簇数的算法DBSCAN用密度结构自动决定分几类风电-负荷场景里常见的离群曲线也能单独识别出来最后输出的是一组带概率的确定性代表场景可以直接喂给后续优化模型。需要说明一下项目标题没有给出详细的原始代码或者数据说明下面讲到的参数标定方法、代码框架和工程细节都是基于我在电力系统场景分析里最常见的做法补全的。整体思路和DBSCAN的核心逻辑是不变的你可以直接照搬改数据。1. 为什么场景缩减偏偏选DBSCAN问题背景与算法选型分析1.1 风电-负荷场景缩减到底解决什么问题场景缩减解决的核心矛盾是随机规划里“场景数量-求解规模-精度损失”三者之间的冲突。拿风电并网调度来说一个典型的两阶段随机优化模型如果每个时段风电有10种可能出力24小时就有10的24次方种组合这显然没法直接算。工程上通常的做法是先通过预测误差采样或者历史数据抽样生成几百上千个初始场景然后通过场景缩减技术挑选出少量代表性场景每个场景附带一个概率值。后面的优化、评估、决策全部基于这些代表场景展开。这里有个关键点叫“确定性场景缩减”。它的意思是说缩减之后得到的不是一堆模糊的分布区间而是具体的、确定性的风电出力和负荷曲线每条曲线都代表一类典型运行工况。DBSCAN聚类的做法正好契合这个逻辑每一簇代表一类相似的风-荷组合模式取簇内代表性样本作为典型场景概率按簇内样本占比计算输出结果就是一组“确定性曲线概率”的组合。1.2 K-means、层次聚类与DBSCAN为什么偏偏是密度聚类做场景缩减最常见的三个候选算法是K-means、层次聚类和DBSCAN。先看K-means它最大的问题就是你得先知道最终要分成几类。我最早做场景缩减踩的第一个坑就是拍脑袋定簇数定少了把早晚高峰迥异的负荷场景硬揉成一团定多了又出现一堆只有两三个样本的小簇优化模型里根本没意义。虽然可以用肘部法则、轮廓系数来辅助选K但它对初始质心敏感而且对不规则分布的场景数据适应性一般。层次聚类不需要预指定簇数但它的计算复杂度偏高而且一旦合并或分裂操作完成就没法撤销。DBSCAN的思路完全不同它不关心你要分几类它只问一句话哪些场景的密度是连在一起的密度够的聚成一簇密度不够、周围稀稀拉拉的场景直接标记为噪声。这种思路对风电-负荷数据有天然优势——风电出力经常存在明显的聚类特性比如“大风天”、“小风天”、“爬坡天”这些模式之间的边界并不规整完全不是球形分布K-means处理这种非凸形状的簇结构会比较吃力而DBSCAN可以聚成任意形状。1.3 DBSCAN在风电-负荷场景中的适配点与取舍DBSCAN对场景缩减有三个明显的适配性优势。第一不需要预先指定簇数。场景缩减的簇数由数据密度结构自动决定你可以根据输出的簇数和噪声比例来判断参数标定是否合理而不是先拍脑袋定“我就要10个场景”。第二内置噪声识别机制。风电和负荷数据里经常有极端场景比如台风天出力骤降、负荷尖峰异常这些曲线往往不是主流模式但它们可能恰恰是规划里不能忽略的边界情况。DBSCAN会把这类离群曲线识别出来你可以在缩减时选择丢弃也可以单独保留作为极端场景研究。第三参数可解释性相对强。Eps和MinPts两个参数直观一个控制邻域半径一个控制密度阈值配合k-距离图可以比较清晰地标定。当然也必须承认它的取舍。DBSCAN在数据规模较大的时候如果直接算全量距离矩阵O(n²)的复杂度会比较吃力。另外它处理高维数据的表现受距离度量影响较大。针对风电-负荷这种24小时×2维的曲线数据合理做数据规约和维度处理能有效缓解这些问题。整体上对于中小规模几百到几千条曲线的风电-负荷场景缩减DBSCAN是比K-means更稳的选择。2. 数据准备与预处理场景缩减的第一步往往最容易被忽略2.1 初始场景集的生成方式做场景缩减之前得先有一批初始场景。我这次用的是简化模拟数据生成方式方便演示流程设定一个有典型日特性的风电基准曲线夜间出力高、白天出力低一个带早晚高峰的负荷基准曲线然后叠加上随机扰动。实际工程中你不会这么搞常见的做法是从风电、负荷的历史预测误差分布中做蒙特卡洛采样或者直接对历史场景做重采样。rng(42); nScen 1000; % 初始场景数量 nHour 24; % 调度时段数 t (0:nHour-1); % 风电基准曲线模拟夜间高、白天低的典型风电日特性 windBase 0.55 0.25 * cos(2*pi*t/24*1.5 - 0.8); windBase max(0.05, min(0.95, windBase)); % 负荷基准曲线模拟8点和19点左右两个高峰 loadBase 1000 220 * exp(-((t-8).^2)/10) 280 * exp(-((t-19).^2)/12); % 叠加随机扰动生成初始场景集 windScen windBase 0.18 * randn(nScen, nHour); windScen max(0, min(1, windScen)); % 风电出力限制在[0,1]区间 loadScen loadBase 90 * randn(nScen, nHour);这里需要留意的是随机扰动用了正态分布实际上风电预测误差更常见的是带偏态或重尾的分布比如beta分布或者混合高斯分布。如果你手上的数据是从预测系统里拿到的误差统计最好用实际误差分布来采样这会影响后续聚类对极端场景的识别能力。另外初始场景数不必盲目追求大1000和5000对聚类结果的影响远不如参数标定和预处理大。2.2 量纲归一化风电与负荷必须放在同一度量下这是整个流程里最容易被新手跳过、但对结果影响最大的一个步骤。风电出力是0到1之间的标幺值负荷可能是几百到几千兆瓦的绝对数值如果不做归一化两个量纲放在一个48维向量里欧氏距离几乎会被负荷这一维完全主导。最后的聚类结果实际上就是“负荷曲线聚类”风电信息形同虚设。我常用的做法是分别把风电和负荷各自线性归一化到[0,1]区间然后拼成一个48维向量。拼完之后风电通道和负荷通道在距离计算中的权重基本均衡。不要用全局min-max因为风电和负荷分布区间不同全局缩放依然会让数值大的负荷主导距离。% 分别归一化到[0,1] windMin min(windScen(:)); windMax max(windScen(:)); loadMin min(loadScen(:)); loadMax max(loadScen(:)); windNorm (windScen - windMin) / (windMax - windMin); loadNorm (loadScen - loadMin) / (loadMax - loadMin); % 拼接成 nScen x 48 的特征矩阵 dataMat [windNorm, loadNorm];还有一种思路是给风电和负荷各分配一个权重系数比如你更关心风电不确定性就可以把风电通道放大。这个取决于具体业务需求不是越精细越好保证两者数量级一致、不互相淹没是底线。有人说“权重怎么定”我的经验是先跑一版等权归一化看聚类结果中风电和负荷的信息保留度再决定要不要调整权重。2.3 距离度量与高维曲线场景的特殊性DBSCAN的核心依赖是“距离”所以距离怎么算必须想清楚。默认用的欧氏距离直观也高效但存在两个问题。一个是高维诅咒。48维空间里样本之间的距离会趋向均匀化也就是说很多点之间的距离差异不大Eps参数的敏感性会上升。针对风电-负荷场景我的建议不是盲目上马氏距离或者DTW而是先考虑降维。比如对24小时风电曲线提取均值、峰值、波动幅度等特征或者用PCA把48维压到10维以内。这不仅能缓解高维问题还能降低计算量。我这次为了演示核心流程还是保留48维全曲线但你在实际工程中可以根据场景形态选择先做PCA或者特征提取。另一个是曲线错峰的误判问题。如果两天的风电曲线形态完全一样只是峰值时段错开两三个小时欧氏距离会判定它们相距很远但物理上它们可能代表同一类天气过程。解决思路有两个一是用DTW距离它对曲线伸缩有容忍度但计算复杂度高、Matlab里实现和DBSCAN结合不太方便二是先做维度规约把时序曲线压成特征从根上降低错峰影响。对大多数场景缩减场景欧氏距离配合归一化和降维处理已经够用不必为了完美而牺牲效率。3. 核心参数标定Eps与MinPts的工程化整定方法3.1 先理解密度可达核心点、边界点与噪声点DBSCAN里三个基本角色一定要先弄明白核心点邻域半径Eps内至少包含MinPts个样本点它代表某个密度区域的内核。边界点本身邻域内样本数少于MinPts但它落在某个核心点的邻域内归为该簇的边界。噪声点既不是核心点也不在任何核心点的邻域内不对应任何簇。把一堆曲线想象成一片云层上的不同“雨团”核心点就是雨团中心雨滴密集边界点就是云团边缘稀疏一点但仍连着噪声点是孤立的零星水汽怎么都归不进去。理解了这个模型你就知道参数标定的目标是什么——让核心点恰好覆盖“真正的长尾场景”让边界点自然衔接各个典型模式噪声点则帮你识别出那些需要单独关注的稀有工况。3.2 k-距离图选Eps最有工程价值的一个步骤Eps选多少直接决定簇数多少和噪声占比。选小了场景被打得稀碎本来同类天气的出力曲线被拆成五六个小簇选大了不同物理特征的风-荷模式被强行捏合成一团场景缩减的区分度就没了。工程上最常用的标定方法是k-距离图。具体做法对每个样本点找到它到第k近邻的距离排序后画出来。这里的k取MinPts-1。曲线会出现明显的“肘部”拐点附近的距离值就是合适的Eps候选值。逻辑很简单在密度均匀的区域k-距离变化平缓一旦进入稀疏区距离快速拉大这个“从平缓转陡峭”的位置就是密度结构的边界。% 计算距离矩阵 D pdist2(dataMat, dataMat, euclidean); % 取 k MinPts - 1 5先算出每个样本到第6近邻的距离 k 5; sortedD sort(D, 2); kDist sortedD(:, k 1); % 降序排列并绘制 [kDistSorted, ~] sort(kDist, descend); figure; plot(kDistSorted, LineWidth, 1.2); grid on; xlabel(样本索引按距离降序); ylabel(第 k 近邻距离); title(k-距离图用于选择 Eps 参数); % 自动寻找拐点区域基于相邻差分最大位置 dk abs(diff(kDistSorted)); [~, idxElbow] max(dk); epsEst kDistSorted(idxElbow); hold on; yline(epsEst, r--, Eps估计值);k-距离图的解读不要死抠数学拐点要结合物理场景看。我通常会在图上把曲线的“膝盖”区域放大然后在这个范围内试两三个Eps值分别跑一遍聚类看簇数和噪声占比哪个更符合业务预期。自动找拐点的代码只是一个快速起点不是最终答案。3.3 MinPts的经验取值与参数联调MinPts控制密度的“最小规模”一个区域至少要有多少个点聚集在一起才称得上一个簇。取太小比如MinPts1任何孤立点都自成一簇聚类形同虚设取太大比如MinPts50小簇全部被合并或被标成噪声极端场景很容易全丢。经验上MinPts一般取数据维度的2倍左右数据量为1000左右时取4到10是常见区间。我建议从MinPts2×dim开始试对48维数据就是MinPts≈6到8。如果噪声点比例偏高可以适当减小MinPts让更多边缘场景有资格入簇如果簇数过碎则增大MinPts。参数标定的完整流程我建议这么做先固定MinPts2×dim画k-距离图确定Eps候选区间。在候选区间内取3个Eps值偏小、居中、偏大分别跑DBSCAN。记录每组参数下的簇数、噪声比例、最大簇样本数。结合业务需求选择参数组合希望保留更多典型模式选偏小的Eps希望场景高度概括选偏大的Eps。微调MinPts观察噪声比例变化控制在10%以内。4. Matlab代码实现从聚类到代表性场景提取的完整流程4.1 数据准备与DBSCAN聚类主程序Matlab从R2019a开始自带dbscan函数属于Statistics and Machine Learning Toolbox调用方式很简单输入数据矩阵、Eps、MinPts即可。如果没有这个工具箱需要用自定义实现我也提供一个精简的教学版本逻辑清晰可以直接用。先说自带函数的用法minPts 6; epsVal epsEst; % 用上一节k-距离图标定的值 % 自带dbscan函数返回簇号-1表示噪声非负整数表示簇类 if exist(dbscan, file) [idx, coreD] dbscan(dataMat, epsVal, minPts); fprintf(自动检测到工具箱使用自带dbscan\n); else [idx, ~] mydbscan(D, epsVal, minPts); fprintf(未检测到工具箱使用自定义DBSCAN实现\n); end fprintf(聚类簇数%d\n, max(idx)); fprintf(噪声点数量%d占比%.2f%%\n, sum(idx -1), sum(idx -1) / nScen * 100);如果没有统计工具箱可以用下面这个精简版自定义函数。它基于距离矩阵实现先以邻域内样本数判断核心点再用队列扩散完成簇的合并和标准DBSCAN逻辑一致。function [idx, isCore] mydbscan(D, eps, minPts) n size(D, 1); idx zeros(n, 1); % 0未访问-1噪声0簇编号 isCore false(n, 1); cluster 0; for i 1:n if idx(i) ~ 0 continue; end % 找所有在Eps邻域内的样本 neighbors find(D(i, :) eps); if numel(neighbors) minPts idx(i) -1; % 暂时标记为噪声 continue; end % 新建一个簇 cluster cluster 1; idx(i) cluster; isCore(i) true; seed neighbors; % 待扩散队列 pos 1; while pos numel(seed) j seed(pos); pos pos 1; if idx(j) -1 idx(j) cluster; % 边界点属于当前簇 continue; elseif idx(j) ~ 0 continue; % 已在其他簇中跳过 end idx(j) cluster; nb find(D(j, :) eps); if numel(nb) minPts isCore(j) true; seed [seed, nb]; % 核心点继续扩散 end end end end这段代码有两个细节想说明一下。第一个是边界点被标记为-1后在后续BFS扩散中遇到它时会被重新归类到当前簇这正好符合DBSCAN的定义。第二个是整个流程里没有使用“visited”标志而是直接用idx的非零状态作为访问标记逻辑上是自洽的。4.2 代表性场景提取选质心最近样本而不是直接取质心聚类完成后的关键决策是怎么从每个簇里挑出那一条“代表性场景”。做法一直接取簇内样本均值作为代表场景。优点是曲线平滑、代表了簇的“平均形态”缺点是平均会抹掉一些物理上重要的波动细节比如风电爬坡的斜率、负荷峰值的尖锐程度。如果你后面的优化模型对曲线形态敏感直接用均值可能引入偏差。做法二取簇内离质心最近的原始场景作为代表。这样得到的曲线是真实存在的场景保留了原始曲线形态而且和质心非常接近几乎等价于“平均形态”但又不失真。我一般选这个做法尤其在工程报告里更好解释——“这条曲线代表了一类典型场景”。代码实现如下nCluster max(idx); repIdx zeros(nCluster, 1); prob zeros(nCluster, 1); validMask idx 0; validTotal sum(validMask); for c 1:nCluster members find(idx c); prob(c) numel(members) / validTotal; if numel(members) 1 repIdx(c) members; continue; end % 簇内所有样本到质心的距离 center mean(dataMat(members, :), 1); dm sum((dataMat(members, :) - center).^2, 2); [~, imin] min(dm); repIdx(c) members(imin); end % 取出代表场景的原始曲线未归一化 repWind windScen(repIdx, :); repLoad loadScen(repIdx, :);这里有一个容易错的点在归一化空间里算质心和距离代表场景编号拿到后要去原始数据矩阵里取曲线还原千万别在归一化的数据里取完再反变换。因为反变换时一旦极小值和最大值搞错曲线还原就会出偏差。直接用原始矩阵索引取既简单又不会错。4.3 概率分配噪声点到底算不算分母场景缩减最后输出给优化模型的不只是曲线本身还有每条曲线的概率。我计算概率的分母用的是有效样本数剔除噪声点后而不是原始场景总数。原因在于噪声点被识别为离群场景后它们不代表任何一类“典型模式”把它们留在分母里会让每个代表场景的概率被人为稀释反而扭曲了典型场景的真实占比。如果你的噪声点比例很小5%以内分母用哪个影响不大如果噪声比例超过10%你更应该回头检查参数标定而不是纠结分母。如果你希望保留极端场景有两条路一是把噪声点单独输出为“稀有场景”列表不参与聚类概率分配只在后续分析中单独研究二是调小Eps让那些极端场景各自成簇。但各自成簇的坏处是簇内样本数很少概率权重很低对优化模型影响有限还需要你额外判断是否重要。我通常倾向于前者参数标定尽量让噪声比例控制在一个合理范围内同时单独把噪声场景保留下来做敏感性分析。4.4 缩减效果评价与可视化聚类缩减做完了怎么判断效果好坏不画图直接说“簇数挺合理”是不够的。我至少会做两组评价。第一组是可视化对比。把原始场景集中随机抽若干条灰色细线画出来再把代表场景用彩色粗线叠加能直观看到代表曲线是否覆盖了原始场景的分布主体。第二组是数值指标用概率加权后的缩减场景均值、标准差和原始场景的均值、标准差做对比误差越小说明缩减对统计特征的保真度越高。% 原始场景的统计特征 meanOrigWind mean(windScen, 1); stdOrigWind std(windScen, 0, 1); % 缩减后场景的统计特征概率加权 meanRedWind sum(repWind .* prob, 1); stdRedWind sqrt(sum(((repWind - meanRedWind).^2) .* prob, 1)); % 平均绝对误差 errMean mean(abs(meanOrigWind - meanRedWind)); errStd mean(abs(stdOrigWind - stdRedWind)); fprintf(均值曲线MAE%.4f标准差曲线MAE%.4f\n, errMean, errStd); % 绘图对比风电场景 figure(Position, [100 100 1200 420]); subplot(1,2,1); hold on; plot(1:nHour, windScen(1:200,:), Color, [0.75 0.75 0.75]); plot(1:nHour, repWind, LineWidth, 1.8); xlabel(小时); ylabel(风电出力 (p.u.)); title(风电场景缩减对比); legend({原始场景(抽样显示), 代表场景}, Location, best); grid on; subplot(1,2,2); hold on; plot(1:nHour, loadScen(1:200,:), Color, [0.75 0.75 0.75]); plot(1:nHour, repLoad, LineWidth, 1.8); xlabel(小时); ylabel(负荷 (MW)); title(负荷场景缩减对比); grid on;这里要注意均值误差低不代表场景缩减质量一定好因为它只看一阶矩。更严格的做法是计算缩减前后场景集的Wasserstein距离不过计算复杂度高一些工程上一般用均值方差对比作为快速判据。如果你特别关心极端场景保真度可以额外对比原始场景集和缩减场景集中最大出力、最小出力的分位数差异。5. 场景缩减实践中的常见问题与避坑技巧5.1 高频问题排查速查表实际跑起来你可能会遇到下面这些情况。我把最典型的几类整理成一张对照表排查起来很方便。现象可能原因处理办法簇数量过多、过于零碎Eps取值偏小密度连不上在k-距离图拐点右侧取更大的Eps噪声点占比过大超过15%MinPts偏大或Eps偏小减小MinPts或适度增大Eps场景聚类结果几乎全部合并成一类Eps太大密度过度连通往k-距离图拐点左侧取更小的Eps聚类结果只反映负荷特性风电特征消失数据未归一化或风电权重被淹没分别归一化检查风电通道的方差贡献曲线形态差别很大但被分为同一簇欧氏距离对时序错峰太敏感先做PCA降维或提取均值/方差/峰值等特征计算时间过长距离矩阵O(n²)太大缩减初始场景数或调用自带dbscan的KDTree加速不同批次聚类结果不稳定数据扰动大、参数在临界区微调Eps稳定在拐点明显区域或增大MinPts5.2 噪声点比例是一个极其有用的诊断信号我自己使用的一个经验指标噪声比例超过10%先别急着往下做回去调参数。这不是说噪声点有错而是高噪声比例通常说明Eps或者MinPts没有匹配好数据的密度结构。有一次我做风电场景缩减初始场景集里夹了不少风电突降的极端曲线一开始噪声比例高达18%。我以为是正常的结果一检查发现是Eps偏小那些极端曲线和主流曲线之间距离明显但没被连接起来。把Eps调大之后部分极端曲线成了边界点进入邻近簇噪声降到6%左右场景缩减效果明显更合理。但反过来如果噪声比例已经降到3%以下说明几乎所有场景都被硬塞进了某个簇这可能掩盖了真实的离群模式。做规划时极端出力场景恰恰是系统最紧张的时刻全丢了风险很大。所以我的建议是噪声比例控制在3%-10%之间少于3%要警惕过度拟合多于10%要重新审视参数。5.3 场景缩减结果到底该看什么三条判断标准很多人缩减完看到簇数和噪声比例觉得没问题就直接把结果丢给优化模型了。我建议再多加三个检查。第一条看代表场景的物理合理性。风电曲线是否落在合理出力区间负荷曲线是否保留早晚高峰特性。如果某个代表场景的曲线形态扭曲到不符合物理常识大概率是参数或者数据预处理出了问题。第二条看均值方差的保真度。前面已经给了代码这里再强调一点不光要看风电也要看负荷的一二阶矩误差。因为场景缩减输出的是一组联合场景风电和负荷的协方差结构如果被破坏对随机优化的影响很大。第三条看代表场景的累积概率分布。你可以把原始场景集和缩减场景集的某一个关键指标比如日均风电出力、日峰值负荷的累积分布函数画在一起两条曲线越贴近说明缩减保真度越高。这一步对说服团队里不熟悉算法的同学特别有用直观且有力。5.4 关于“要不要保留噪声点”的一个工程建议直接丢弃噪声点是最简单的做法适合噪声比例小的情况。但遇到以下两类情况我会单独保留噪声点第一类是噪声点里含有负荷尖峰或者风电极端低出力场景。这类场景虽然发生概率低但一旦发生对电力系统的安全运行影响极大。丢弃它们在概率上是合理的但在工程决策里不可接受。我会把它们单独列出来作为“极端场景集”提供给规划人员做校核。第二类是噪声点数量虽然少、但曲线形态有清晰物理含义。比如一整簇“风电骤降负荷攀升”的曲线因为距离别的主流簇太远而被打成噪声但它们实际上代表一个明确的运行风险模式。这种情况下我会调整Eps或者直接用HDBSCAN这类带层次结构的密度聚类算法把这簇曲线单独识出来。5. 结尾一点个人实操体会说实话DBSCAN这个算法我最早接触的时候是在处理交通轨迹聚类后面迁移到电力场景缩减才发现它在风-荷联合场景里格外顺手。整个过程做下来最深的体会是场景缩减的瓶颈从来不在算法本身而在参数标定和对业务场景的理解。k-距离图我建议每次必画哪怕你已经很有经验了因为它能迅速帮你判断当前数据的密度结构是不是和预期一致。最后再分享一个小技巧跑完聚类之后把每个簇的样本数打出来看看是不是均匀分布。如果出现一个簇占了60%以上的样本另一个簇只占2个样本那说明Eps选得过大或过小场景缩减的区分度不够。这种时候不要硬调重新画一次k-距离图往往几分钟就能找到合适区间。场景缩减这个东西说白了就是用最小代价保留最核心的不确定性信息DBSCAN只是实现这个目标的一个趁手工具。