
做新能源电力系统随机优化的人大概都经历过这种痛苦一上不确定性风机出力、负荷曲线就变成一堆场景几百上千条调度模型转头就跑不动了。场景削减要干的活就是在这堆场景里挑出几个最有代表性的把概率分布特性尽量留下来同时把计算规模压下去。我最近在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和minPtsDBSCAN的原理写进论文就三四句话但真正调参时坑不少。我尽量用大白话把核心概念和参数选择讲透并给出能直接在MATLAB里跑的K-距离曲线脚本。2.1 核心点、边界点、噪声点一个生活化理解DBSCAN把每个场景点分成三类。核心点以它为中心、半径eps的邻域内至少有minPts个点含自身边界点落在某个核心点的邻域内但自己邻域内不够minPts个点噪声点既不满足核心点条件也不挨着任何核心点。聚类过程是这样的随便找一个未访问的核心点出发把它邻域内的所有点拉进来再对这些点的邻域继续扩张所有密度相连的点就形成一个簇。边界点被并进它所属核心点的簇里但不再往外扩张。噪声点不参与任何簇。这里有个生活化的类比想象把人按“聚在一起的程度”分组。核心点是小区里人缘好、邻居多的住户边界点是住在核心点隔壁、自己周围却没什么邻居的人噪声点就是方圆几里只有他一家。DBSCAN把所有挨着核心点的住户连成一片社区“孤零零”的房子就排除在外。这个直觉比公式好记也直接决定了参数调整的方向。2.2 eps和minPts到底怎么取值两个参数直接决定聚类结果。eps太小大部分点都变成噪声eps太大所有点融成一团聚类等于没做。minPts的常见经验取法有两个一是取维度数加1dim1这是区分真实簇与随机噪声的下限二是取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对角线是0sort之后第1列永远是0所以这里其实取到了包括自身在内的前minPts个点相当于minPts-1个真实邻居。工程上影响不大如果想要严谨可以用sortedDist(:, minPts1)。我两种都试过曲线形状几乎一致拐点处差异不明显不用过于纠结。画完图找拐点曲线先是一段接近水平的段然后在某个位置突然翘头。翘头起点对应的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维minPts8根据K-距离曲线拐点取eps0.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也可以在聚类前对风电分量和负荷分量分别归一化后再拼接。这个坑比较隐蔽因为聚类结果“看起来”能跑出来但代表场景里的风电特征明显失真。5.5 最后分享几个片段式的经验第一次做联合场景削减时我把风、负荷分开各自聚类再合并场景结果削减后的联合场景完全破坏了风荷相关性下游调度模型天天出现违反爬坡约束的结果。后来改成一个联合向量一起聚类问题立刻消失。所以做成联合场景的时候就别偷懒拼在一起降维、聚类、选代表一次搞定。另外MATLAB的dbscan函数默认把所有点分成若干簇加一个噪声簇但如果不留意它会把两个靠得比较近的簇合并成一个。遇到这种情况把eps稍微调小0.01到0.02试试边界就会分开。这个微调手法比重新调一大堆参数省事得多我经常这么干。概率分配我一直用簇内占比代表场景的原始概率为1/N每个簇概率相加总概率加上噪声概率还是1逻辑自洽。个别文献里会再按Kantorovich距离优化权重能把削减精度再提一些不过常规工程场景必要性不大。这套流程我现在已经沉淀成项目里的标准工具了新来的人拿着就能用基本不用再调。