做风电和负荷不确定性建模的同学手上大概率都攒着这么一批东西几百上千个场景。说白了每一个场景就是一段风电出力曲线加一段负荷曲线代表未来可能发生的一种“日子”。场景越多随机规划模型越精确可求解也越难受——一个几十节点的机组组合模型叠上一千个场景基本就没法玩了。所以咱们得在精度和算力之间找个平衡这就是场景削减要干的活。这几年做风电-负荷场景削减主流方案无非是K-means、层次聚类、同步回代消除再就是用DBSCAN这类密度聚类。DBSCAN不是新东西1996年就被提出来了但在电力系统里用它做场景削减近几年才慢慢多起来。我最初拿到“基于DBSCAN密度聚类风电-负荷确定性场景缩减”这个需求时第一反应是一个密度聚类算法凭什么适合做场景削减后来真跑数据才发现它和风电场景的契合度确实比传统K-means高不少。这篇文章把整体思路、Matlab实现、参数调法和踩坑记录都写出来用我实际跑过的流程做蓝本希望能帮你少走几步弯路。全文不是算法教材更像是项目笔记。你如果正在做随机优化、机组组合、配电网规划需要把几百上千个场景压到十几个还能保住代表性这篇文章应该能对上路。1. 为什么风电-负荷场景非要“削减”不可——先把这个需求看清1.1 场景多的本质不确定性要用“样本”来逼近风电出力和负荷都有很强的随机性。风速受气象影响负荷受温度、作息、重大事件影响二者叠加起来未来一天的情况有无数种可能。实际工程里我们不可能把这些可能全部塞进优化模型而是先通过采样生成一批典型场景用这批场景去近似真实分布。最常见的做法是拉丁超立方采样或者基于历史数据的Bootstrap重采样。假设风速用Weibull分布描述负荷用正态分布叠加时间相关系数描述采样1000次就能得到1000个“24小时风电出力曲线”和1000个“24小时负荷曲线”。每个曲线就是48维向量风电24维负荷24维或者按15分钟粒度更多维这1000个向量构成了场景集。场景数越多模型对不确定性的刻画越精细但计算复杂度往往成倍增长。两阶段随机规划里决策变量和约束数量随场景数线性膨胀到了求解机组组合、优化调度这类整数规划问题时1000个场景基本是灾难级的求解时间甚至直接内存溢出。所以场景削减不是“要不要做”的问题而是“怎么在信息损失最小的前提下做”的问题。1.2 削减到底削掉什么场景削减的本质是从N个原始场景中选出M个代表场景M远小于N并给每个代表场景分配一个概率权重。削减后的场景集要满足两个基本要求代表性保留的场景能反映原始场景集的整体分布形态包括均值、方差、相关性、极端情况。概率合理性场景削减不是等权平均而是每个代表场景要承担一个“概率质量”让它近似代替原始场景在该区域的密度。DBSCAN在这个需求下有一个天然适配点它按密度划分场景区域落在同一高密度区域的场景合为一个簇簇中心作为代表场景簇内样本数量占总体样本的比例可以直接作为概率权重。同时那些偏离所有簇的孤立场景会被标记为噪声这部分场景在物理上往往对应极端气象或特殊运行工况——它们虽然概率小但可能是调度最紧张的时刻单独拎出来反而更有价值。1.3 “确定性场景”到底指什么标题里的“确定性场景”容易让人犯迷糊。我的理解是场景削减的输入是概率性的海量样本输出则是有限的几条“确定的、不含概率分布细节的典型曲线”每条曲线配一个权重。后续优化模型把每条代表场景当作确定性的边界条件来求解不需要再叠一层随机波动的描述。换句话说不确定性被“离散化”到了一组确定场景里这组场景就是随机规划的输入模板。搞清楚这一点很重要因为它决定了你在代码层面的输出结构削减结果应该是“代表场景矩阵权重向量”而不是一堆聚类标签。很多人跑完聚类就结束了忘了概率权重的处理这就是后面要专门说的坑。2. 方案选型为什么我选了DBSCAN而不是K-means2.1 K-means的几个硬伤K-means是场景聚类最常用的方法但它有几个地方在风电-负荷场景里很别扭。第一个问题是K值必须预设。你到底需要20个场景还是30个场景业务上往往没有先验。K选小了代表性不足K选大了削减效果不明显。虽然可以画肘部图来找K但风电-负荷联合场景的簇结构通常不清晰肘部图经常给出模棱两可的结果。第二个问题是K-means对簇形状的假设太强。它假设每个簇是凸的、球状的样本按欧氏距离分配到最近的质心。而风电场景在二维降维投影后经常出现带状、弧形乃至嵌套形态的簇。硬要用K-means的情况就是在一个细长簇的两端各放一个质心把一个真实的连续状态空间硬切成两半人为地制造出不存在的“模式”。第三个问题是K-means对噪声点完全不设防。初始场景里总会有少数“怪场景”比如凌晨风速骤升、负荷突变叠加外部故障。此类离群点参与均值计算会把质心拉扯偏移导致代表场景“不干不净”既不像正常场景又带走了部分概率质量。2.2 DBSCAN为什么适合这个场景DBSCAN的核心思想很朴素给定半径Eps如果一个点周围半径内的邻居数量超过MinPts它就是一个核心点核心点之间通过密度可达关系连通成簇孤立在密度区域之外的点就是噪声。放到风电-负荷场景里这个逻辑非常直观场景大量集中在“风平浪静”的常规区域密度高自然形成几个大簇。少数场景位于区域边缘密度低要么并入附近簇要么成为噪声。簇的数量和形状由数据自己决定不需要预先指定K。它还有一个隐含好处簇内样本的分布可以是非球形的。风电出力受风速上限限制负荷受峰谷规律限制联合场景在特征空间里天然不是球形分布DBSCAN能比K-means更老实地贴合真实形态。2.3 和其他场景削减方法的横向对比传统场景削减里还有一类很有代表性的方法叫同步回代消除也常叫基于概率距离的快速前向选择。它的思路是反复合并距离最近的两个场景把一个场景的概率加到另一个上面直到场景数降到目标值。这类方法能严格保持概率距离理论性质漂亮但计算复杂度偏高而且每次合并都只能局部操作对全局密度结构不敏感。层次聚类也可以用但它需要指定树状图的截断高度来得到簇数调参一样麻烦而且大规模场景下距离矩阵的计算和树构建都比较慢。我用一个表格总结一下几种方法在我实际使用中的体感方法是否预设簇数对噪声处理计算复杂度场景削减适配度K-means需要不设防噪声拉偏质心线性很快一般适合簇结构简单的场景层次聚类需要截断高度没有显式噪声机制较高样本量大时吃力一般小样本可用同步回代消除需要会把小概率场景合并到邻近场景较高迭代合并较好理论严谨但工程实现麻烦DBSCAN不需要显式标记噪声较低配合索引更快好能自动找簇数适合大量场景单从效果讲同步回代和DBSCAN都是可用的方案。但我个人更倾向于DBSCAN原因很实际代码量小、参数少、不用迭代、一次聚类就能把簇和权重都算出来。尤其是初始场景数量上了5000、10000之后同步回代那种O(N²)的迭代让人等得心烦DBSCAN配合数据索引能快一个数量级。3. Matlab代码实现从场景生成到聚类削减的完整链路3.1 场景数据怎么组织在写聚类代码之前先把数据结构定清楚。假设我们要削减的是“24小时风功率24小时负荷”的联合场景时间粒度取1小时那么每个样本是一个48维的向量。把所有场景堆起来得到一个 N×48 的矩阵XN是初始场景数。这个排列顺序很关键。如果你把风功率放在前24列、负荷放在后24列那么之后做特征归一化、距离计算都要基于这个固定顺序。我习惯在代码开头写清楚每一列的含义注释否则隔两周再看代码自己都容易懵。矩阵X的每一行代表一个场景行与行之间互相独立每一列代表某个时段的风功率或负荷值。聚类算法处理的单位是“行”它计算的是两行之间的相似度也就是两条“24小时曲线”之间的接近程度。这个前提想清楚了后面看代码才不会绕晕。3.2 归一化处理不做就等于白做这是整个流程里最容易翻车的地方。风电功率一般用标幺值范围在0到1之间负荷可能是几百兆瓦范围在几十到上千不等。如果不做归一化直接把两个量拼成48维向量欧氏距离会被负荷维度牢牢控制风功率那24维的变化在距离计算里几乎不起作用。结果就是聚类结果只反映了负荷模式的差异风电模式再怎么变都不会影响分簇。我常用的归一化方式有两种具体选哪种要看你对物理量纲的重视程度全局z-score对X全体做均值0、方差1的标准化。简单但会把风电和负荷混在一套统计量里。分段z-score对风电部分24列和负荷部分24列分别做标准化。这样两类变量对距离的贡献比例更均衡我项目中默认用这个。代码如下% X: N x 48 矩阵前24列为风电后24列为负荷 windPart X(:, 1:24); loadPart X(:, 25:48); windPart (windPart - mean(windPart(:))) / std(windPart(:)); loadPart (loadPart - mean(loadPart(:))) / std(loadPart(:)); X_norm [windPart, loadPart];如果你希望风功率的标幺值物理意义不丢失也可以不做z-score而是把风功率除以额定容量、负荷除以基准负荷让两者都变成0-1区间的标幺值。这种做法在工程汇报里更好解释但要注意两种变量量纲不一致时距离计算里隐含的权重比例仍是人为设定的。还有一个细节如果做的是单纯的时序场景而且你特别关心曲线形状而不仅是幅值可以考虑对每条曲线先做幅值归一化减均值除标准差再做聚类。这样能把“波形相似但幅值不同”的场景聚到一类。但风电场景的幅值本身就有物理意义幅值大小直接关系到发电量所以我不建议对风电部分做这个操作只有负荷曲线的形状聚类才值得这么干。3.3 自己写DBSCAN还是用内置dbscan函数MathWorks的Statistics and Machine Learning Toolbox提供了一套不错的dbscan函数核心调用就一行idx dbscan(X_norm, epsilon, minpts);idx是N×1的标签向量簇编号从1开始噪声点标为-1。函数内部用knnsearch做邻居搜索对普通规模的数据速度还行代码也简短。但我还是建议你在项目里自己实现一遍DBSCAN。原因有两个第一内置函数对NaN非常敏感输入数据只要有一个NaN整个结果就可能变成全噪声或者直接报错排查起来很费劲。而自己实现的时候你可以控制每一环的数据清洗。第二内置函数返回的信息有限拿不到“哪些点是核心点”也没法方便地自定义距离度量。自己实现一遍后续想扩展DTW距离、想加权重、想看某个簇的核心点范围都更自由。自己实现一个基本版DBSCAN并不复杂核心逻辑就三步function idx myDBSCAN(X, eps, minpts) N size(X, 1); idx zeros(N, 1); % 初始化全为未分类0表示未访问 clusterId 0; % 预计算距离矩阵样本量大时可用分块方式替代 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 % 新建簇 clusterId clusterId 1; idx(i) clusterId; seed neighbors; k 1; while k numel(seed) j seed(k); if idx(j) -1 idx(j) clusterId; % 边界点归入当前簇 elseif idx(j) 0 idx(j) clusterId; % 如果这个点也是核心点则扩展种子集 neighborJ find(D(j, :) eps); if numel(neighborJ) minpts seed union(seed, neighborJ); end end k k 1; end end end这段代码是教学版的清晰实现。实际项目中我不会直接拿find(D(i,:) eps)硬算因为当N超过5000时pdist2生成的就是一个5000×5000的double矩阵占200MB内存到10000样本就是800MB。我常用分块策略每次只加载一个分块的距离行或者直接调用knnsearch只保留每个点的前K个近邻信息然后用邻接表形式做密度扩展。不过在功能验证阶段上面的简单版足够你跑通流程、理解算法行为。3.4 从聚类结果到“确定性削减场景”聚类标签只是中间产物真正要输出的是代表场景矩阵和概率权重。操作方法如下对每个簇包括噪声点单独作为一个极端场景簇看你的需求计算簇内所有样本的均值向量作为代表场景权重等于簇内样本数除以总样本数。% idx为聚类标签噪声点为-1 uniqueClusters unique(idx(idx 0)); % 只统计非噪声簇 numClusters numel(uniqueClusters); weights zeros(numClusters, 1); reprScenes zeros(numClusters, size(X_norm, 2)); for c 1:numClusters members (idx uniqueClusters(c)); weights(c) sum(members) / sum(idx 0); reprScenes(c, :) mean(X(members, :), 1); % 注意原始空间里求均值 end这里要提醒一句**最后求代表场景要用原始未归一化的X而不是X_norm。**因为归一化只是为了让聚类算法合理划分簇输出给优化模型的场景必须回到真实的物理量纲否则风功率0.3和负荷500MW对调度代码来说没有任何意义。另外背景场景的权重归一化问题我习惯把所有簇的权重除以非噪声样本占比让去除噪声后的权重和为1。如果你需要把噪声场景也纳入后续随机优化可以把噪声点单独归成一个“极端场景簇”用均值做代表、样本数定权重这样不会丢失小概率但高影响的工况。4. 调参与省钱经验Eps和MinPts到底怎么取4.1 MinPts和Eps的关系DBSCAN两个核心参数一个是邻域半径Eps一个是最小邻居数MinPts。它的物理含义可以理解为一个点周围至少要有多大的“人群密度”它才算是一个有代表性的场景这个密度区域的范围又有多大。MinPts的经验值一般是维度数的2倍也就是2×d。但风电负荷场景维度动辄24、48、96维如果严格按2×d取MinPts会导致邻居数要求过高几乎所有点都变成噪声。实际项目中我会把MinPts当作一个独立的“密度阈值”而不是维度函数来调常用的起始值是5到15。比如初始场景2000个MinPts10表示一个场景周围得至少找到10个邻居它才被认作是区域的中心代表否则就只是个边缘样本。Eps则完全取决于数据的尺度。先做归一化Eps的合理区间通常在0.5到2.0之间。但具体值还是要看数据本身的密集程度不能拍脑袋定。4.2 KNN距离图选Eps选Eps最经典的做法是画KNN距离排序图。思路是把每个点到它第MinPts个最近邻居的距离算出来按从小到大排序画成曲线曲线从平缓到陡峭的拐点处所对应的纵坐标就是Eps的参考值。% MinPts 10 k 10; [idxSort, distK] knnsearch(X_norm, X_norm, K, k); % 取每个点到第k个邻居的距离 kDist sort(distK(:, k)); figure; plot(1:length(kDist), kDist); xlabel(Points sorted by distance); ylabel(K-th nearest neighbor distance);曲线拐点出现的原因很直观处于同一高密度区域的点彼此离得近第k个邻居的距离小离群点或者稀疏区域边界上的点第k个邻居的距离大。排序之后密集区域的距离值连成一个低平的平台稀疏区域的距离值突然拉高这个突变点就对应着合适的密度分界。不过要泼一盆冷水维度高了以后这个拐点往往被抹平了。48维空间里绝大多数点对的距离都挤在一个窄区间里KNN曲线可能从头到尾都没有明显拐点。这时候我会换个思路扫描Eps值计算每次聚类的轮廓系数取轮廓系数最高的那个Eps。这个方法虽然在维度高时不完美但至少能给出一个可复现的参考。4.3 轮廓系数扫描在Matlab里轮廓系数用silhouette函数计算它衡量每个样本落在自己簇内的紧凑度和与其他簇的分离度范围从-1到1。均值越高聚类质量越好。epsRange 0.2:0.1:2.0; bestScore -inf; bestEps epsRange(1); for epsVal epsRange idxTmp dbscan(X_norm, epsVal, minpts); % 需要至少2个簇才有轮廓系数的意义 if numel(unique(idxTmp(idxTmp 0))) 2 continue; end score mean(silhouette(X_norm, idxTmp)); if score bestScore bestScore score; bestEps epsVal; end end这个扫描的逻辑很简单但计算量不小。每个Eps都要跑一次silhouette而轮廓系数本身要算两两距离两千个样本×二十个候选Eps跑起来要起码十几分钟。我的经验是先用KNN距离图缩小Eps搜索范围再在候选区间内做扫描能省掉大半时间。4.4 结合业务场景定参数调参不能只盯着数据指标还得回到物理场景。我实际调参时有个习惯先看目标削减数量。如果你希望最终得到15到20个代表场景那么可以先扫一遍聚类结果看哪个Eps组合下簇的数量恰好在这个范围。如果某个Eps下簇数太少比如只有5个说明密度阈值太松如果出现几十个簇多半是太严了。还有个经验是分两阶段调先用轮廓系数粗选Eps再人工抽查几个簇的代表曲线看簇内场景的物理一致性。如果同一个簇里出现了“风功率早高峰”和“风功率午高峰”这类明显不同模式的曲线说明Eps设大了把本应分开的模式糊在了一起。这种物理检查比任何数学指标都可靠。5. 案例实测1000个场景削减到18个误差多少5.1 测试数据构造为了验证方法我构造了一组模拟数据用Weibull分布生成风速套用典型风机功率曲线转成风功率叠加时序平滑负荷用“基础负荷上下班峰谷随机波动”生成。一共1000个联合场景每个场景48维24小时风功率24小时负荷。归一化后调用dbscan(X_norm, 0.8, 10)第一轮聚出21个簇噪声点数96个。略作调整后把Eps微调到0.75得到18个簇噪声点降到43个。这个过程自然体现了DBSCAN和K-means的一个差异——你不需要提前定K稍微动一下Eps/MinPts簇数会自动适应数据。5.2 削减结果评估削减效果我通常从两个维度看一个是分布保留程度另一个是优化模型的输入偏差。分布保留上最直接的是看削减后场景集在风功率均值和负荷均值上的偏差指标削减前1000场景削减后18场景偏差风电均值标幺0.4230.4171.4%负荷均值标幺0.6180.6220.6%风电标准差0.1870.1718.6%负荷标准差0.1420.1364.2%这个结果在我跑过的项目里算中等偏好。均值保住了方差的损失也在可接受范围内——毕竟要用18条曲线替代1000条曲线的分布方差被低估是必然的。削减比例越高方差损失越明显这是物理规律不是算法缺陷。再看代表性场景的形态。从18个代表场景里拿出3个看一个是“大风大负荷”场景概率权重0.12一个是“小风小负荷”场景权重0.09还有一个是“大风小负荷但早高峰突出”场景权重0.06。这几类场景的划分有明显物理含义后续放进机组组合模型里可以直接对应不同的调度策略。5.3 和K-means对比同样的1000个场景用K-means强制聚成18类对比结果很有意思。K-means在均值偏差上和DBSCAN差不多风电均值偏差1.1%甚至略好一点但它在“簇形态合理性”上差不少——有几个簇的中心曲线明显是两种不同模式被硬平均出来的“中间产物”曲线形状怪两头矮中间平看起来就不像真实场景。计算耗时上dbscan调用大概3秒K-means迭代5到10次大概2秒差距不大。但如果把同步回代消除拉进来1000个场景迭代到18个场景花了将近40秒。所以大规模场景下我的排序是DBSCAN优先K-means备选同步回代只在小样本且预算充足时用。6. 实际运行踩坑记录与排查手册6.1 内存爆炸距离矩阵太大这是最经典的坑。自己实现DBSCAN时最容易写D pdist2(X, X)然后N5000时内存直接飙升。有一次我在一台16GB内存的机器上跑10000样本Matlab当场卡死任务管理器里内存一路飙红。解决思路有几个一是用pdist配合squareform只保存压缩后的距离矩阵二是用knnsearch只查每个点的前K个近邻然后用稀疏矩阵构造邻接关系三是分块计算距离每次只处理一片。实际工程中我优先推荐直接调用内置dbscan它内部用knnsearch做了索引对大数据集友好得多。6.2 聚类结果全是噪声或者全是一个簇这种现象通常不是算法问题而是参数问题。全是噪声说明Eps太小或MinPts太大密度条件太苛刻全是一个簇说明Eps太大密度条件太宽松把整个空间都连成一片了。排查方法很简单先画KNN距离图看看实际的邻居距离分布在哪里。如果大部分点的第10个邻居距离在0.6左右你Eps取0.2那自然全是噪声取4那自然全是一簇。从KNN图的平台区域开始尝试通常一两次就能落在合理区间。还有一个隐蔽原因没有做归一化。风电和负荷量纲差异极大时欧氏距离几乎由负荷维度主导KNN距离图也失真Eps的绝对大小变得不可解释。这种时候即使怎么调Eps聚类结果都看起来不对劲。先归一化再看KNN图顺序不能颠倒。6.3 场景概率权重没有归一化这个坑最隐蔽因为它不算代码报错而是结果悄悄出错。有人聚类完之后直接拿weights 簇内样本数 / 全部样本数当概率但忽略了噪声点。如果噪声点占比5%那么所有非噪声簇的权重之和就是0.95这会导致后续随机规划模型里所有场景概率加起来不为1。正确做法有两种要么把噪声点单独设为一个簇并纳入权重归一化适合极端场景保真要么直接只用非噪声样本重新计算权重适合作为主场景集。我一般选择后者把噪声场景剔除让权重和为1再单独建一个“极端场景档案”供调度人员参考两不耽误。6.4 高维距离失效问题风-荷联合场景维度高48维起步欧氏距离在高维空间里会趋向“均匀化”密度差异不再明显。这是我实际观察到的现象KNN距离图的拐点越来越难找轮廓系数普遍偏低聚类结果对参数极其敏感。应对措施有三个层级最粗暴的是不做处理接受K-means式结果折中做法是先PCA降维到10维左右再做DBSCAN精细做法是改用余弦距离或动态时间弯曲DTW距离但后者计算开销大不适合大规模场景。我的建议是优先尝试PCA降维。风电和负荷时序本身高度相关48维的真实有效维度可能只有5到8维。降维后DBSCAN的密度结构会清晰很多聚类也更加稳定。不过降维后的代表场景要映射回原空间再输出操作上别忘记反变换。6.5 从方法论角度提醒一句跑得多了以后我有个体会场景削减的核心是“先想清楚削减后的场景拿来干什么再选算法和参数”。如果你只把DBSCAN当黑盒调包参数调半天也不知道合不合理如果你先确定下游模型对场景数量、极端场景覆盖度、概率精度的需求再回头定簇数规模和噪声处理策略整个流程就会顺畅很多。我最初在几个项目里折腾了好几天没出结果后来就是把“物理检查”放到调参里去问题反而自动消解了。这个习惯建议你也试试看。
