DBSCAN场景削减法：风电-负荷随机规划从5000条场景到15条，求解加速20倍

上周我把一个跑了很久的随机规划模型彻底跑明白了，瓶颈不在求解器，而在于前期那5000条蒙特卡洛风电-负荷场景。原始场景一多，两阶段随机优化的规模直接失控，内存和求解时间双双爆炸。后来我用DBSCAN密度聚类做了一套确定性场景缩减方法，在Matlab里完整实现了一遍，把5000条场景压到15条，优化模型求解时间从四十多分钟降到两分钟，目标函数期望值的偏差控制在2%以内。这篇整理就把我从选型、参数设置到代码实现的完整过程写出来，供正在做风电-负荷不确定性分析、随机规划，或者单纯需要从大量时序场景里提取代表性场景的同行参考。

先说结论：DBSCAN在场景削减这个任务里真正的优势不是“聚类效果好”，而是它天然不关心“到底该保留几个场景”。你只需要给出密度参数，它自动把场景分成若干个密度相似的簇，剩下的噪声场景要么剔除、要么单独处理。这一点对工程场景极其友好，因为“削减到多少个场景”通常不是你先验决定的，而是需要根据下游求解器的承受能力来回调整的。K-means那种必须先定K的玩法，在这个场景下很容易让人磨掉一下午。

整套方法的链条其实不长：蒙特卡洛生成大量风荷场景 → 提取特征并标准化 → DBSCAN聚类 → 在每个簇内挑一个离中心最近的原始场景作为代表 → 按簇内场景数量分配权重 → 用概率距离和优化目标双重验收。我按这个顺序逐个说清楚细节。

1. 为什么风电-负荷场景一定要先削减再进优化模型

1.1 场景是怎么生成的，以及它为什么会长成那个样子

风电出力场景的生成逻辑，通常是从风速的概率分布出发。风速一般用Weibull分布描述，抽样后用风功率曲线映射成出力，再叠加预测误差；负荷场景则是在历史负荷曲线的基础上叠加不同幅度的预测误差。反复抽样N次，就得到一个n×T的矩阵，每一行是一条场景，每列是一个时段（典型的有24点或96点）。初始情况下，每条场景的概率都是1/n。

这里要强调一个认知：场景不是“未来一定会发生的运行状态”，而是对不确定性的离散近似。你没法把所有可能都装进优化模型，只能挑有代表性的。这个“代表性”三个字，就是整个场景削减问题的出发点。

1.2 场景数量与优化问题的“放大效应”

两阶段随机规划里，每新增一条场景，第二阶段就会新增一组连续变量和约束。场景数从10涨到5000，求解开销不是线性增长的，约束矩阵的稀疏结构会变差，迭代次数也会明显增加。我在某个带储能和可调负荷的两阶段模型中实测过：5000条场景时单次求解要四十多分钟，内存占用逼近机器上限；削减到15条场景后两分钟左右跑完，数值稳定性也好了不少。

这还没算调试成本。每跑一次四十多分钟，调参根本没法进行，更别说做灵敏度分析。因此场景削减不是“锦上添花”的可选项，而是随机规划落地的必经步骤。

1.3 削减的本质：概率分布的近似，而不是“挑几条好看曲线”

很多入门者把场景削减理解成“挑几条形状看起来有代表性的曲线”，这是最大的误区。削减后的m条场景及其概率权重(y_k, w_k)，要在某种概率距离下尽可能近似原始n条场景的联合分布。你要保留的是均值、方差、时序相关性、风电与负荷之间的耦合关系、极端事件出现的概率，而不是某条曲线的形态好不好看。

理解这一点，才能理解为什么评估削减效果要用Wasserstein距离而不是肉眼观察，也才能理解为什么后面我们要花不少功夫在特征构造上。

2. 为什么这次我选了DBSCAN而不是K-means

2.1 K-means/K-medoids在场景削减里的三个硬伤

第一个硬伤，需要预设簇数K。在场景削减问题里，K就是“最终保留多少个代表性场景”。这个数往往由下游优化模型的可承受规模决定，而不是由数据分布直接给出。K-means等于逼你先回答一个你自己都回答不了的问题。第二个硬伤，K-means对初始中心敏感，多次运行可能给出不同划分，结果不可复现。第三个硬伤，风电场景里存在不少“异常出力”样本，比如极端大风日、骤降爬坡日。K-means的质心是均值意义上的中心，一条极端离群场景会把簇中心拉向自己，导致整个簇的代表性被污染。

有人会说K-medoids比K-means稳健，但K-medoids仍然需要预设簇数，而且在大样本下计算开销并不低。场景削减场景里，样本量动辄几千上万，这个方法并不占优。

2.2 DBSCAN的密度逻辑对场景数据很友好

DBSCAN的核心定义需要先理清：核心点是在eps半径内至少有minPts个邻居的点；边界点是落在核心点邻域内、但自身邻居数不足的点；噪声点则是既非核心点、也不在任何核心点邻域内的点。用社交圈类比：核心点是圈子里人脉广、有影响力的角色，边界点是被圈子接纳但存在感不高的成员，噪声点则是社交孤岛。

对风电场景来说，这个逻辑有一个直接的好处：一条场景周围如果聚集了大量相似场景，说明它代表了一类高频出现的运行模式；孤立型场景（比如某次极端天气下的出力样本）自然被划为噪声，不会去污染代表场景。这正是场景削减想要的效果——高频模式要细致刻画，低频极端模式要么单独保留、要么剔除后做概率重分配。

2.3 “确定性场景缩减”到底确定在哪里

标题里的“确定性”三个字，我理解有三层含义。第一层：给定eps和minPts后，DBSCAN对同一数据集的聚类结果是唯一且可复现的，多次运行完全一致，不依赖随机初始化。第二层：簇数由数据密度自动涌现，不需要人为指定，也不需要在多个K值之间反复试。第三层：削减流程是确定性的——每条场景要么落在某个簇里，要么成为噪声，归属清晰，权重可以按簇内成员占比精确计算。

与之对比，K-means属于随机初始化方法，结果会有波动；同步回代削减（backward reduction）的结果虽然是确定的，但计算复杂度高，每条新样本都要全局扫描配对。DBSCAN在“大规模场景粗削减”这个定位上，确实很合适。

2.4 和经典同步回代削减的横向对比

实际项目中，我倾向于把DBSCAN当成“粗削”工具，先把几千条砍到几十条，再用同步回代削减把这几十条做进一步精削。两者是互补，不是替代。如果你只需要一个能快速落地、结果稳定的方案，单用DBSCAN这条链路也完全够用。

3. DBSCAN用在风电场景上，参数和特征处理得当才能出效果

3.1 三个关键变量：eps、minPts、距离度量

DBSCAN的可控变量就三个：邻域半径eps、核心点最少邻居数minPts、距离度量。默认距离度量是欧氏距离，但对场景数据来说，欧氏距离并不是唯一选择。相关性距离（1减Pearson相关系数）和DTW能更好捕捉形状相似性，不过计算代价高。我的经验是：如果特征提取做得够好，欧氏距离在标准化后的特征空间里已经够用；如果场景形态特别复杂，再考虑DTW，但那时候往往需要换一套聚类思路了。

3.2 特征构造：别直接拿96点原始序列硬聚类

直接拿96点原始时序做聚类的问题非常明显。一是维度高，样本在稀疏空间里密度估计不稳定；二是风功率序列相邻时段强自相关，绝大多数维度是冗余信息；三是欧氏距离对幅度敏感，聚类结果容易被“平均出力大小”主导，两条形状几乎一样、但幅值略有差异的场景会被直接切开。

我推荐的做法是：先提取统计特征，再做降维。常用的特征组合包括均值、标准差、峰谷差、最大爬坡率等；也可以把96点按时段聚合成24点或8点，降低维度；如果特征维度仍然偏高，再用PCA降到5到10维。联合场景则要把风电子向量和负荷子向量都纳入特征，不能只看风或只看荷。

我实测过一个例子：直接用96点原始出力聚类，轮廓系数不到0.35，簇边界模糊；换成统计特征加PCA降到8维后，轮廓系数超过0.6，簇结构清晰很多。轮廓系数虽然不是最终目标，但至少能直观反映簇的可分离性。

3.3 eps的工程化选取：k-距离图

确定eps的经典方法是k-距离图。操作步骤：

1. 先给minPts一个候选值，一般取5到10，或者略大于特征维数；
2. 对每个样本，计算它到第(minPts-1)近邻的距离；
3. 把所有距离降序（或升序）排列，画成折线图；
4. 找到曲线由平缓转为陡峭的“拐点”，拐点对应的距离就是eps的候选值。

Matlab里画图代码很简单：

k = minPts - 1; [~, distK] = knnsearch(X, X, 'K', k + 1); % knnsearch返回的K列中包含样本自身，自身距离为0，所以要取最后一列 distK = distK(:, end); distK = sort(distK, 'ascend'); plot(distK); xlabel('样本序号'); ylabel('第k近邻距离');

图里平坦的部分说明绝大多数场景密度相似，拐点之后距离突然变大，那些点就是离群场景。拐点处的距离作为eps，既能保留核心密度，又能把离群点甩出去。这个图是我每次换数据集必跑的第一道工序。

3.4 标准化：这步不做，聚类结果基本没法看

特征矩阵里各列量纲差异很大：出力均值是几十MW，标准差可能只有几MW，峰谷差又是几十MW。如果不做标准化，聚类几乎等于按数值最大的那一列排序，其他特征全部失效。这一步没有太多技巧：

X = zscore(F);

一句话解释：标准化让每个维度在距离计算中拥有同等话语权。需要注意，如果特征里有常量列（标准差为0），zscore会生成NaN，这种情况要提前剔除该列，或者给常量特征加一个极小扰动。

4. Matlab代码实现：从原始场景集到代表性场景与概率权重

4.1 主流程：六步走

整个代码链路分成六步：导入场景矩阵P_wind和P_load → 构造特征矩阵F → zscore标准化 → 调用dbscan聚类 → 从每个簇里挑出离特征中心最近的原始场景 → 按簇内场景数分配概率权重并输出。

4.2 特征提取与标准化示例

假设已经有两条矩阵：P_wind是n×Tw的风电出力场景，P_load是n×Tl的负荷场景。下面这段代码构造了一个8维特征矩阵，前6列分别刻画风电和负荷的平均水平与波动程度，第7列刻画爬坡极端度，第8列刻画风荷耦合关系。

n = size(P_wind, 1); F = zeros(n, 8); F(:,1) = mean(P_wind, 2); % 风功率均值 F(:,2) = std(P_wind, 0, 2); % 风功率波动 F(:,3) = range(P_wind, 2); % 风功率峰谷差 F(:,4) = mean(P_load, 2); % 负荷均值 F(:,5) = std(P_load, 0, 2); % 负荷波动 F(:,6) = range(P_load, 2); % 负荷峰谷差 % 风功率归一化爬坡指标：最大相邻时段变化 ramp = max(abs(diff(P_wind, 1, 2)), [], 2); F(:,7) = ramp ./ (F(:,1) + 1e-6); % 加小常数防止除零 % 风荷相关系数 for i = 1:n c = corrcoef(P_wind(i,:), P_load(i,:)); F(i,8) = c(1,2); end F = zscore(F);

第8列的循环在样本量极大时会比较慢，可以换成矩阵化写法，但语义不如循环清晰。如果特征维度继续增加，建议在这步之后接一个PCA，保留累计解释方差90%以上即可。8到10维对DBSCAN来说已经是很舒服的输入维度。

4.3 DBSCAN聚类与代表场景选取

Matlab从R2019a版本开始，在统计和机器学习工具箱里提供了dbscan函数。调用方式很直接：

% 需要 Statistics and Machine Learning Toolbox，R2019a及以上版本 [idx, corepts] = dbscan(F, eps, minPts); clusterIds = unique(idx(idx > 0)); % 噪声点在idx中被标为-1 nClusters = length(clusterIds); repIdx = zeros(nClusters, 1); weights = zeros(nClusters, 1); for c = 1:nClusters members = find(idx == clusterIds(c)); center = mean(F(members, :), 1); dists = sum((F(members, :) - center).^2, 2); [~, pos] = min(dists); repIdx(c) = members(pos); % 代表场景在原场景集中的索引 weights(c) = length(members) / n; end % 提取代表性原始场景 repWind = P_wind(repIdx, :); repLoad = P_load(repIdx, :);

这里有一个容易忽略的决策点：为什么选“离簇内特征中心最近的原始场景”，而不是直接拿簇中心当代表场景？因为簇中心是特征空间里的合成点，映射回原始场景空间很可能是一条被“平均”得很平滑的假曲线，现实中未必会发生。场景削减要忠于原始样本，所以应该在簇内选一条真实存在的场景作为代表。这也是我在多个项目里坚持的做法。

权重分配也要注意：上面的代码假设每条初始场景等概率，权重就是簇内成员占全体样本的比例。如果你的原始场景本身带非均匀权重（比如用重要抽样法生成），那weights(c)应该改为该簇内场景权重之和，最后再归一化。

4.4 噪声场景的取舍：一个值得想清楚的问题

DBSCAN会把孤立极端场景标记为噪声。噪声场景直接丢掉的代价，是丢失极端天气或极端负荷事件的发生概率，而这些事件在可靠性分析里恰恰不能忽视。我在实际项目里用过“保留并压权”的策略：噪声场景不参与聚类，但也不直接丢弃，而是合并成一个“边界场景”保留下来，权重取噪声场景数量占全体样本的比例。另一种常见做法是彻底剔除噪声场景，再把剩余各簇权重重新归一化，得到“去掉异常后的典型运行状态”，这对纯经济调度问题一般够用。

两种策略在评估阶段会暴露差异：如果削减后极端事件的概率被明显低估，就选保留策略；如果模型对尾部风险不敏感，直接剔除问题也不大。

4.5 老版本Matlab没有dbscan函数怎么办

R2019a以下版本没有内置dbscan。替代方案有三个：自己基于距离矩阵写一个，核心逻辑也就几十行；先用pdist2计算距离矩阵，再对每个样本找邻域；或者调用第三方实现。场景量在几千条规模时，自己实现完全没压力。唯一要注意的是距离矩阵的规模：5000条场景的距离矩阵是5000×5000，double类型大约占200MB，这个内存花销在普通台式机上还能接受。

5. 削减效果评估：用数据说话，而不是用眼睛看图

5.1 概率距离指标：一阶Wasserstein距离

场景削减的最终目标是分布近似，评估的黄金指标是概率距离。一阶Wasserstein距离俗称“推土机距离”，物理含义是把一个分布“搬运”成另一个分布的最小代价，它天然考虑概率权重。对每个时段t，把原始场景在该时段的取值排成经验分布，削减后场景带权重，然后计算两个分布在[0,1]区间上分位数函数的绝对偏差积分。Matlab里可以用一段很短的函数近似：

function wd = wasserstein_1d(orig, rep) q = linspace(0, 1, 101)'; f1 = quantile(orig, q); f2 = quantile(rep, q); wd = sum(abs(f1 - f2)) / 100; end

这里的orig和rep分别是原始场景集和削减场景集在某一时段的取值列向量。如果削减场景带权重且想精确一点，可以用加权分位数近似，或者把代表场景按权重扩展后再调用quantile。边缘Wasserstein距离虽然只是近似联合分布，但工程上足够暴露大部分问题。更完备的联合分布距离需要计算最优传输，代价高，通常只在论文级对比里使用。

5.2 统计量保真度检查清单

除了概率距离，我还会做一组统计量对比，按下面的清单逐项核对：

方向上有一个很实用的经验：如果某个统计量偏差过大，不用推翻整个方案，回到特征矩阵补一列对应的特征再聚类就行。比如发现峰谷差分布失真，就在特征里加强峰谷差相关维度；发现相关性失真，就在特征里显式加入耦合维度。

5.3 端到端验证：把削减场景扔进优化模型

只谈分布误差还不够，最终要看优化模型的目标函数期望值。把原始场景集和削减场景集分别代入同一个两阶段随机规划模型，比较目标函数期望值的相对偏差，以及关键决策变量（比如储能充放电功率、联络线功率）的分布区间。

我在某项目里的验收标准是：目标期望偏差小于2%，关键决策变量的95%置信区间重合度不低于80%。如果这两个条件都满足，这套削减方案就可以直接替换全场景集投入使用。这个方法比任何聚类内部指标都有说服力。

5.4 为什么我不建议只用“轮廓系数”宣传效果

很多聚类文章喜欢汇报轮廓系数、DBI这些内部指标，但在场景削减这个任务里，聚类内部指标再好看，和下游优化结果之间也没有严格对应关系。轮廓系数的用途是参数粗筛——网格搜索eps时快速打分是可以的，但不能作为最终验收信号。项目汇报时如果只拿轮廓系数说事，大概率会被追问“这对优化结果有什么影响”，场面会很难看。

6. 实操中踩过的坑：DBSCAN做场景削减的四点经验

6.1 eps的敏感度远超你想象

稍微改一点eps，簇数就可能从15跳到8，或者从15跳到46。原因是场景数据的密度分布是连续的，没有天然断崖。解决办法是把k-距离图当作初筛工具，缩小范围后做小规模网格搜索：对每个eps候选值调用dbscan，再用概率距离或轮廓系数打分。两轮下来基本能锁定一个稳定区间。我当时就是靠这个办法，把eps从“拍脑袋的0.5”调到了真正合适的0.82。

6.2 距离矩阵的内存与算力

dbscan内部需要计算样本间距离。5000条场景的距离矩阵约2500万个元素，double类型占200MB，如果样本量升到10000，那就是800MB，很多机器直接卡死。对策有三个：一是先降维再算距离，特征矩阵保持8到10维，计算成本会明显下降；二是先抽样预热，用500到1000条场景粗聚类估计密度参数，再全量聚类；三是样本量过万时，考虑换用基于邻域图的变体或者先分块处理。场景削减的数据规模通常还在5000以内，但提前意识到这个问题能省去很多麻烦。

6.3 标准化和特征选择错一步，聚类结果就失真

特征里如果混入量纲差异巨大的列，聚类基本等于按最大的那一列排序。zscore标准化之后还要检查常量列，标准差为0时会生成NaN，需要先剔除。另一个常见坑是聚类的输入特征空间和代表场景空间不一致：你按特征空间聚类，却在原始曲线空间按欧氏距离找“簇内最中心”，两者错位会让代表场景的“中心性”定义失真。我的习惯是统一在特征空间里选代表点，再映射回原始场景矩阵。

6.4 削减后的相关性失真很隐蔽

有的场景集削减之后，单时段边缘分布误差不大，但风荷相关系数明显变化，原因就是特征矩阵里没有等价的联合特征。解决办法是在聚类特征里显式加入耦合维度，比如风荷相关系数、峰谷差比、风电出力占比。我在4.2节那8个特征里专门加入第8列，就是为了避免这个问题。另外，如果噪声场景被直接剔除，削减后场景集的方差和尾部分位数通常会偏低，凡是关注可靠性或极端事件成本的项目，一定要检查尾部分位数或VaR类指标，不要只盯着均值。

我在这个项目里最终定型的组合是：96点场景 → 统计特征加风荷耦合特征 → z-score标准化 → dbscan聚类（eps用k-距离图初选再网格精调，minPts取8） → 簇内最接近中心的真实场景作代表 → 权重按占比计算 → Wasserstein距离和目标期望值双重验收。整套流程换数据集时，基本只需要重跑一遍k-距离图，参数不用大调。

最后再分享一个判断标准：场景削减不是聚类精度比赛，轮廓系数再高也没人关心。下游优化模型跑出来的结果误差，才是衡量这套削减方法好不好的终极裁判。你要是能把削减前后的目标函数期望值误差控制在2%以内，这套方案就可以放心交付了。