风电出力的间歇性和负荷的随机抖动是电力系统规划调度中绕不开的麻烦。做多场景随机优化时我们往往先通过蒙特卡洛或者拉丁超立方采样生成上千个「风电-负荷确定性场景」可真正把它塞进混合整数规划模型后问题规模瞬间膨胀到求解器直接“摆烂”。这时候就需要场景削减从上千个原始场景里挑出几十个有代表性的场景并给它们配上发生概率让后续优化既保持精度又跑得动。这篇想聊的就是如何用DBSCAN密度聚类实现风电-负荷联合场景缩减并给出可复现的Matlab代码适合正在做新能源并网分析、随机机组组合或储能容量配置的同学参考。场景削减的方法很多比如快速前向选择、k-means聚类、层次聚类而DBSCAN是一个经常被忽视的选项。我试用下来的感受是它能把离群的极端场景识别出来不需要预先指定簇数但参数调起来也确实让人头疼。所以这篇文章不只是丢代码而是把“什么时候用DBSCAN”“怎么调参”“怎么评估削减效果”一起讲透。1. 场景削减为什么是随机优化的“必经之路”1.1 一张场景集背后是爆炸的决策变量假设我们需要做未来一天的机组组合决策每小时一个时段24个时段。如果只考虑风电和负荷的不确定性对每个不确定源采样1000次那就会得到1000个可能出现的“未来”。如果把这1000个场景全部保留机组组合模型中的每个场景都要复制一份变量和约束变量规模直接放大上千倍。对于含整数变量的混合整数规划这基本是灾难——求解器可能几天都收敛不了。但如果我们先把1000个场景缩减到20个代表场景模型规模就缩小了50倍计算时间大大缩短而解的质量损失通常可以控制在可接受范围内。很多人以为场景削减只是为了“省算力”其实它也是为随机优化提供离散概率分布的必然手段。随机规划模型本质上是求一个决策在所有可能出现的不确定场景下都表现良好的方案而这个“所有可能出现”必须用一组有限的带概率场景来逼近。场景削减就是在这组有限场景的质量和数量之间做平衡削得太狠误差大削得太少又回到算力灾难。1.2 所谓“确定性场景缩减”到底在缩什么标题里的“确定性场景”并不是说风光出力变成确定值而是指每一个采样出来的场景内部是一条确定的时间序列。比如一个场景就是“第t时段风电出力等于0.45 p.u.、负荷等于87MW”这样一条完整曲线。场景缩减做的事情是把大量的离散确定场景合并成少量场景并为每个保留场景分配一个概率值。最终得到的是一组新的离散分布每个保留场景都是一条确定曲线概率之和为1。后续随机优化模型只需要对这个离散分布取期望或者做风险度量计算量小得多。这个概念很关键因为很多人会把“场景生成”和“场景削减”混在一起。场景生成是把连续的概率分布离散化场景削减是在已经生成的离散集合里做二次压缩两者解决的问题完全不同。1.3 聚类法削减的通用套路用聚类做场景削减时无论采用哪种聚类算法流程几乎都是四步把每个场景转成一个高维向量比如24小时风电曲线加24小时负荷曲线拼起来就是48维向量定义场景之间的距离最常用的是欧氏距离用聚类算法把相似场景归堆从每个簇里选出代表场景用簇内样本数占总样本数的比例作为该代表场景的概率。DBSCAN只是改变了第三步的做法前两步和第四步是共通的。所以如果你之前用过k-means做削减换成DBSCAN只需要把聚类函数换掉代表场景选取逻辑完全可以复用。我之所以推荐大家先理解这个套路是因为很多实现细节其实和聚类算法关系不大核心难题反而在距离度量、代表场景选择和概率分配上。2. DBSCAN密度聚类凭什么能用于场景缩减2.1 密度聚类的原理用“抱团”代替“画圈”DBSCAN的核心思想可以概括成一句话如果一个点周围足够稠密它就拉着附近的人一起成为一个簇密度不相连的区域就分开孤零零的点则标记为噪声。它只有两个参数邻域半径eps和最小样本数minPts。只要某个点周围eps范围内至少有minPts个点这个点就是核心点核心点互相密度可达的部分连成一个簇不满足核心条件又没有核心点可达的点就是噪声。放到场景削减的语境里“点”就是一个场景“密集”就是大量场景在特征空间中离得很近。风电出力受天气系统影响某个风速层级的场景经常成片出现可能形成不规则的高密度区域而一些极端天气样本数量少、偏离中心自然就成了离群点。DBSCAN不会强行把所有点都塞进簇里这对保留“低概率高影响”的极端场景是有意义的。2.2 相比k-meansDBSCAN的直觉优势k-means做场景削减有三个让使用者不太舒服的地方必须事先指定簇数簇数k与实际场景分布不匹配时误差很大算法假设每个簇近似球形而风电-负荷场景在高维空间里的分布往往是条带、曲线甚至分叉的形式对离群点敏感极端场景会被硬拖进某一个簇扭曲簇中心。DBSCAN则不需要指定簇数能发现任意形状的簇并且会主动把离群点标记为噪声让我们有机会单独处理它们。这些特性让它在处理带有厚尾分布、聚类结构不规整的场景集时比k-means更贴近数据本身的形状。如果你做过电力系统不确定性分析应该能理解极端天气下的风电出力场景往往不是正态分布的尾巴而是完全独立的一小撮点强行用k-means去分会把它们平均成一条平庸的曲线这恰恰是规划人员最不想看到的。2.3 聚类完成之后簇怎么变成缩减场景组DBSCAN跑完会得到若干个簇和一堆噪声点。每个簇需要选出一个代表场景。这里有个容易犯的错误直接用簇内所有场景的平均值作为代表。风电曲线平均之后往往会产生一条光滑、但现实中几乎不可能出现的曲线而且会丢失风电场有功出力时的突变信息。更稳妥的做法是选择簇内“最小平均距离点”也就是该点到簇内其他所有场景距离之和最小的那个真实场景这类点在无监督学习中常被称为簇的medoid。每个簇的代表场景概率等于簇内场景数除以原始场景总数。至于噪声点一般不建议直接丢弃可以单独作为极低概率场景保留也可以根据“离哪个簇最近就并入哪个簇”的原则分配到最近的簇里具体取舍要看后续优化是否重视极端场景。3. Matlab代码落地从造数据到选代表场景3.1 先造一个能复现的风-荷联合场景集在讲解核心代码之前我先给出一段场景生成脚本。真实项目中你手头的场景可能来自历史功率曲线叠加预测误差抽样或者来自气象再分析数据为了演示算法流程这里用基底曲线叠加扰动的方式生成1000个场景每个场景包含24时点的风电标幺值和负荷值。rng(42); nScn 1000; nHour 24; % 基底风电曲线白天风电偏弱早晚偏强 t (1:nHour) - 12; base_wind 0.3 0.5*exp(-(t.^2)/18); % 基底负荷曲线简单日负荷形状晚高峰明显 base_load 60 20*sin(t*2*pi/24) 30*exp(-((t-2).^2)/8); base_load base_load 5; % 基线抬升 % 生成风电场景在基底上叠加噪声 wind_scn repmat(base_wind, nScn, 1) randn(nScn, nHour)*0.08; wind_scn max(0, min(1, wind_scn)); % 生成负荷场景叠加逐小时噪声和全局平移扰动 load_scn repmat(base_load, nScn, 1) randn(nScn, nHour)*5 ... randn(nScn,1)*8; load_scn max(10, load_scn);基底曲线选得不太漂亮也没关系关键是你手头要有两个nScn x nHour矩阵wind_scn里每一行是一条风电日曲线load_scn里每一行是一条负荷日曲线。后面所有处理都基于这两个矩阵。3.2 特征拼接、标准化与距离矩阵有两个常见坑要提前说第一风电是0到1的标幺值负荷是几十到几百MW如果不做标准化负荷分量会完全主导欧氏距离DBSCAN聚出来的簇基本只反映负荷差异。第二标准化应当在拼接之后、计算距离之前完成而不是在生成场景时预先处理成无物理意义的数值。下面的代码把风电和负荷拼成48维特征矩阵做z-score标准化然后计算两两场景的欧氏距离矩阵。X [wind_scn, load_scn]; % 1000 x 48 mu mean(X); sigma std(X); sigma(sigma0) 1; % 防止全零分量导致NaN Xn (X - mu) ./ sigma; D pdist2(Xn, Xn, euclidean); % 1000 x 1000 距离矩阵这里得到的D是对称矩阵对角线全0。如果场景数量很大pdist2会占用大量内存这时候可以先随机抽2000个子样本散点化成簇后再把剩余场景按距离最近归属到已有簇但本文先不展开。3.3 核心缩减函数DBSCAN聚类与代表场景提取Matlab从R2019a开始内置了dbscan函数支持直接输入距离矩阵并指定Distance,precomputed。我用它封装了一个dbscan_reduce函数输入原始场景矩阵、eps和minPts输出代表场景、概率以及每个簇对应的原始场景索引。function [repX, prob, cluster_idx, labels] dbscan_reduce(X, eps, minPts) % X: nScn x nFeature可以是[风电,负荷]拼接后的矩阵 % 标准化 mu mean(X); sigma std(X); sigma(sigma0) 1; Xn (X - mu) ./ sigma; D pdist2(Xn, Xn, euclidean); % DBSCAN需要统计与机器学习工具箱R2019a labels dbscan(D, eps, minPts, Distance, precomputed); cl setdiff(unique(labels), 0); % 0 代表噪声 nClusters numel(cl); repX zeros(nClusters, size(X,2)); prob zeros(nClusters,1); cluster_idx cell(nClusters,1); for i 1:nClusters idx find(labels cl(i)); cluster_idx{i} idx; subD D(idx, idx); subD(1:length(idx)1:end) Inf; % 不考虑自身 avgd mean(subD, 2); [~, j] min(avgd); repX(i,:) X(idx(j),:); prob(i) numel(idx) / size(X,1); end % 概率归一化如果所有点都是噪声概率为空 if ~isempty(prob) prob prob / sum(prob); end end要注意一个细节我在计算簇内平均距离时把对角线置为了Inf。如果不这样做每个点与自己的距离0会把平均距离拉低导致簇内几乎所有场景都有机会成为代表削弱选出“中心点”的效果。这个小改动我第一次调试时没注意结果选出来的代表场景明显偏离簇中心后来才发现是自身距离在作怪。3.4 主脚本一天内跑通整个流程现在就可以把前面的生成函数和缩减函数串起来。假设我们想把1000个场景削减到约30个左右可以先从eps1.0, minPts8开始跑完之后看簇数和噪声比例再调整参数。下面的主脚本会输出缩减前后的曲线图以及每个代表场景的概率。% 假设你已经按照3.1节生成了 wind_scn 和 load_scn X [wind_scn, load_scn]; % 调用削减函数 [repX, prob, cluster_idx, labels] dbscan_reduce(X, 1.0, 8); % 提取代表场景的风电与负荷 rep_wind repX(:, 1:24); rep_load repX(:, 25:48); % 画图风电与负荷的代表场景曲线 figure; subplot(2,1,1); plot(rep_wind, LineWidth, 1); hold on; ylabel(风功率(p.u.)); xlabel(时段); subplot(2,1,2); plot(rep_load, LineWidth, 1); ylabel(负荷(MW)); xlabel(时段); % 输出概率列表 disp(table((1:size(repX,1)), prob, VariableNames, {Scenario,Probability}));这里需要注意dbscan的eps参数作用的距离矩阵是在标准化后计算的所以eps1.0大致对应标准化空间里相差一个标准差的场景会被认为是邻居。如果你的场景数据分布非常稀疏可能需要把eps放到1.5甚至2.0如果非常稠密eps0.6也可能翻车。下一节我会专门讲怎么调这个参数。4. 削减效果怎么看误差指标与参数敏感性4.1 一个简单又直观的损失函数概率加权平均绝对误差削减得好不好不能只看画出来的曲线像不像。我常用的一个指标是把每个原始场景按“它属于哪个簇”归堆用该簇的代表场景代替它然后计算整个场景集上风电和负荷的概率加权平均绝对误差。具体来说对每个原始场景找到它对应的簇中心用它与该簇代表场景之间的逐时段绝对偏差加权平均权重就是该场景的出现概率如果原始场景等概率权重就是1/N。这个指标可以写成function err reduction_error(X, labels, repX, cl) cl setdiff(unique(labels), 0); err 0; N size(X,1); for i 1:numel(cl) idx find(labels cl(i)); for j 1:numel(idx) d abs(X(idx(j),:) - repX(i,:)); err err mean(d) / N; end end end这个函数里我对噪声点的处理偷了个懒如果噪声点存在上面的循环会漏掉它们。实际使用时建议把噪声点合并到最近的簇里或者单独作为极小概率场景否则误差会被低估。我在实际测试时发现把噪声点一律丢弃会让评估误差很好看但后续优化遇到小概率高影响场景时会明显“踩空”这不合算。4.2 一组参数敏感性示意结果DBSCAN削减结果对eps和minPts非常敏感。我拿前面生成的1000个场景做了一组快速测试结果大致如下不同随机种子下数值会有变化仅供参考趋势epsminPts簇数噪声占比概率加权MAE0.680100%无法削减0.88535%0.421.082812%0.231.28415%0.311.04923%0.18观察这些趋势有两点值得留意第一eps过小会导致大量场景被当作噪声簇数骤减甚至归零这时候削减基本无效第二eps增大虽然能减少噪声但会把不同风况的簇并到一起代表场景数量少误差反而变大。minPts控制的是簇的最小规模minPts越小越容易被划成小簇簇数变多误差下降但保留场景过多会削弱削减的意义。所以调参时不能只盯簇数或误差必须把“场景削减后的总数”也放进约束里。4.3 从误差到优化目标误差不是越低下限越高从工程角度看场景削减的误差最终要反映到优化目标上。比如随机机组组合里最关心的是预期发电成本和备用不足概率场景削减造成的误差会转化为目标函数偏差。有一次我用k-means把500个场景削减到10个概率加权MAE很小但一跑备用不足约束发现削减后的期望失负荷指标比原始场景高了近20%。原因就是10个簇太少丢失了尾部极端场景。用DBSCAN时我特意把噪声场景单独保留下来相当于把削减后场景数从10个增加到13个结果失负荷指标被校准了。这说明评估指标不能只用“平均误差”还要看下游优化关心的风险指标。平均误差小只说明“大多数情况下代表场景贴近原始场景”但少数极端场景可能被抹平了。如果你做的是可靠性评估或安全约束经济调度那必须把尾部场景单独拿出来看一遍。5. 调参经验与工程实践中的几个坑5.1 高维时间序列下eps怎么选k距离图DBSCAN最让人头大的就是eps。好在有一种经典方法k距离图。把每个场景到它第minPts近邻居的距离全部算出来升序排列画成一条曲线。曲线从平缓到突然上扬的“拐点”对应的距离就是比较合理的eps候选值。原因是拐点附近的点意味着从稠密区域进入稀疏区域以这个距离为半径正好能区分簇内与簇外。我建议每次调参前先做这一步而不是凭感觉拍一个eps。高维数据下的距离分布经常是“要么完全连成一片要么全是噪声”看到拐点再动手能省很多时间。k 8; % 与minPts保持一致 [idx_dist, ~] knnsearch(Xn, Xn, K, k1); kdist idx_dist(:, end); kdist_sorted sort(kdist); plot(kdist_sorted); % 找“拐点”附近的值作为eps初值注意k1是因为要排除自身距离。拐点越明显说明场景的簇结构越清晰如果整条曲线都接近直线DBSCAN很可能不适合这组数据这时候可以考虑换层次聚类。5.2 标准化顺序和物理意义的一致性前面提到过标准化要在拼接后统一做。这里再多说一句如果只标准化风电部分而负荷部分不标准化或者用各自的均值和方差分开标准化都会改变风电和负荷在整个距离中的相对权重。假设负荷变异系数远小于风电分开标准化会把负荷差异放大聚类结果可能过度关注小时级负荷波动反而忽略了风电不确定性。我见过有同学先在场景生成阶段把风电序列减去均值再拼上负荷结果削减出来全是负荷变化不大的簇就是因为标准化对象搞错了。正确的做法是把风电和负荷作为一个整体样本去缩放这样两个物理量在距离空间中处于同一比较尺度。标准化的另一个好处是让eps有了相对可解释的含义eps1.0大致表示“相差一个标准差内的场景视为邻居”。5.3 欧氏距离 vs DTW要形状还是要计算量风电和负荷曲线都带有明显的时间特性。欧氏距离按同一时刻逐点相减一旦两个场景的峰值时间错开一个小时欧氏距离就被判得很远哪怕它们形状一模一样。如果场景集里存在明显的时移现象可以考虑用动态时间规整DTW距离替代欧氏距离聚类效果会好很多。但DTW有两个代价一是计算复杂度高1000个场景两两算DTW可能要跑到天荒地老二是eps的物理含义变得模糊k距离图不那么好解释了。我的建议是先用欧氏距离和DBSCAN跑通流程如果你发现同类场景看起来“很像”但因为错位被切碎成很多小簇再考虑换DTW。场景削减毕竟是追求整体代表性不是做模式识别大多数情况下欧氏距离够用。5.4 和k-means、层次聚类怎么取舍最后把这几种常见场景削减方法放在一起对比。方法是否需要指定簇数离群点处理高维稳定性主要适用场景k-means需要强行并入簇较好簇结构规整、分布近似球形层次聚类需要看树状图可以单独分支较稳定需要多粒度观察DBSCAN不需要自动识别噪声参数敏感簇形状不规则、噪声明显实操中我倾向于把DBSCAN当作“体检工具”先用k距离图和DBSCAN看看场景集有多少个天然密簇、噪声比例多大再决定用哪种方法做最终削减。如果簇结构清晰且噪声占比低k-means完全够用如果噪声占比超过10%说明场景集中存在大量离群气象事件这时才值得用DBSCAN并把噪声场景单独保留到最终场景组里。5.5 一个容易忽略的细节代表场景要回到原始物理空间聚类过程中我们用了标准化后的特征做距离计算但从簇里选代表场景时一定要选原始空间里的真实场景点而不是标准化空间里的“平均向量”。哪怕你是用medoid方法也要注意repX(i,:)取自原始矩阵X而不是标准化后的Xn。我在第一次封装函数时直接在Xn里选出了中心点然后把这个48维标准化向量还原成风电和负荷曲线结果发现失真非常明显——标准化把原本20倍量纲差异强行压平了还原误差被放大了好几倍。这个坑很好检查输出代表场景的负荷值如果出现负几十兆瓦那就是忘了还原。标准化只服务于距离计算代表场景必须从原始场景里挑这是场景削减的一条底线。我在调试这套DBSCAN场景削减流程时最大的体会是聚类算法永远只是工具理解数据分布、明确下游优化需要什么才是核心。DBSCAN给我们的不是“最佳削减方案”而是一个能自动识别离群场景、不需要预先拍板簇数的观察窗口。如果你打算把场景削减嵌进随机规划代码里我建议保留这个观察窗口先用k距离图看看数据再小步调整eps和minPts最后一定要用下游优化指标做一次端到端验证而不仅仅是看削减误差。这套Matlab代码我后续还在扩展比如加入DTW距离支持以及把离群场景自动并入最近簇的逻辑等稳定下来再继续分享。