做电力大数据分析的同行应该都有体会，居民负荷数据看起来维度高、规律性又强，但实际上隐藏的用电模式很值得深挖。这是我近期用Matlab完整做过一遍的项目：基于粒子群算法（PSO）优化Kmeans聚类，对居民用电行为做分群分析。为什么要用PSO而不是直接跑Kmeans？一句话说清楚：Kmeans对初始聚类中心太敏感，多跑几次可能跳出完全不同的分群结果，而PSO相当于用全局搜索能力先帮Kmeans找到一组更合适的初始中心，再让Kmeans在这个基础上做局部精调。整套流程包括数据预处理、PSO寻优、Kmeans聚类、结果可视化和指标评价，代码全部在Matlab里实现，适合正在做电力客户画像、需求响应、峰谷电价策略的学生和工程人员参考。这篇文章不是把论文复述一遍，而是把我在实际调试中真正踩过的坑、反复试过的参数、最后验证可行的方案都摆出来，希望能帮你少走弯路。
1. 为什么用粒子群算法去优化Kmeans：核心设计思路拆解
1.1 Kmeans聚类的老毛病
Kmeans应该是做数据分群最常用的算法了，原理简单，Matlab里一句kmeans就能跑出结果。但真正拿它做居民用电行为分析时，问题就冒出来了——它本质上是迭代爬山式的局部搜索，聚类中心初始位置选得好不好，直接决定最终聚类结果。
随便打开一个用户曲线数据集跑Kmeans，设K=4，连续跑十次，十次里面的分群结果可能都不一样。原因就在于目标函数是非凸的，Kmeans只能找到距离初始点最近的那个局部极小值。你换个初始中心，它可能就爬到另一个谷底去了。在实际项目中，这就意味着你的分群结果不稳定，今天是这几个用户归一类，明天跑一遍又变了，做需求响应方案时根本没法交代。
这个问题在居民负荷数据上尤其明显。居民用电曲线不像工业用户那样规律，用户之间的相似度差距没那么大，边界模糊，Kmeans就更容易因为初始点选取不当而陷入不理想的分群方案。比如本该区分出“夜间活跃型”和“白天外出型”，结果因为初始中心都落在某个高负荷区域，两类被合并成一类，这就非常尴尬。
1.2 PSO为什么能干这个活
粒子群算法的思路很朴素：想象一群鸟在一片区域里找食物，每只鸟都不知道食物在哪，但大家能共享当前位置的适应度信息。鸟群会同时朝着自己找到过的最好位置和整个群体找到过的最好位置飞，通过反复迭代逼近全局最优。它的本质是群体智能的全局搜索，不依赖梯度信息，也不需要目标函数可导，并且实现起来非常容易——这正是它适合用来做Kmeans初始中心搜索的原因。
把PSO和Kmeans结合，有两种常见路子。一种是完全用PSO替代Kmeans的迭代过程，让PSO直接搜索聚类中心，这种做法的优点是全局搜索能力强，但缺点是粒子维度比较高——几个中心拼接在一起，每个粒子的位置就是K个中心坐标的串联，维度一大，搜索效率反而下降。另一种做法就是把PSO当作前级优化器，先用PSO搜出一组比较优的初始中心，再用这组中心初始化Kmeans，让Kmeans发挥它局部收敛快的优势。我在实际项目中用的是后者，理由很简单：实际调试下来，PSO搜索出的中心已经接近全局较优区域，Kmeans在这个基础上只需要少量迭代就能收敛，整体稳定性和效率都兼顾了。
1.3 整套流程怎么搭
整个方案按“数据→特征→寻优→聚类→评估→可视化”这条线走，每一环都是有讲究的。
第一步是数据预处理。居民负荷数据通常来自智能电表，15分钟一个采样点，一天96点，一个月下来就是2880维，直接用原始数据聚类会让计算量爆炸，而且噪声也会把聚类中心带偏。我会先做数据清洗，去掉掉线和异常采样的记录，再按用户维度聚合出典型日负荷曲线，最后做特征降维。
第二步是特征设计与归一化。从96点负荷曲线中提取日最大负荷、日最小负荷、日平均负荷、峰谷差、峰时用电占比、谷时用电占比等特征，再配合PCA降维保留主要形态信息。所有特征必须做z-score标准化，这一步不做的话聚类结果基本没法看，量纲差异会主导距离计算。
第三步是PSO寻优。定义好粒子编码方式后，设定种群规模、迭代次数、惯性权重、学习因子，以簇内总距离作为适应度函数，让PSO搜索出一组聚类中心。
第四步是Kmeans精调。用PSO给出的中心初始化Kmeans，迭代到收敛，得到最终分群结果。
第五步是评估与可视化。用轮廓系数、DBI指数这些指标验证聚类效果，再把聚类结果映射回原始负荷曲线，画出各类别的典型日负荷曲线，看分群是否具备清晰的物理含义。
这套流程的价值在于：PSO负责“找到好起点”，Kmeans负责“快速收敛到最优点”，两者分工明确，配合默契。在实际应用中，同一套数据跑PSO-Kmeans和直接跑Kmeans做对比，PSO-Kmeans每次都会在同一个比较优的解上收敛，稳定性提升非常明显。
2. 数据预处理与特征构造：聚类质量的前半场
2.1 原始数据的清洗策略
居民负荷数据从采集端到分析端，中间会有各种脏数据。最常见的几类问题：通信中断导致的整段缺失、计量异常导致的瞬时跳变、数据乱序和重复记录。这些如果不处理，后续聚类会直接歪掉。我在项目里最常用的一套清洗流程是这样的。
先做去重和乱序修正。智能电表数据一般带时间戳，按用户ID和时间排序，遇到完全重复的记录保留一条。然后是缺失值处理——如果某一天缺失的采样点少于总采样点的10%，用前后两天的同时刻均值补上；如果超过30%，说明当天数据质量太差，直接把这一天的曲线剔除，不参与后续聚合。
最后是异常值检测。负荷数据的异常主要看两种情况：数值为负（说明接线错误或计量故障）、数值突变（比如10分钟时间内负荷从0.2kW跳变到8kW再跳回）。我用滑动窗口的方法，计算每个采样点与前后各4个采样点的均值差异，超出3倍标准差的点标记为异常，整体替换为该用户同时间的周平均负荷。这样处理后，数据形态就比较干净了。
2.2 特征选择：到底该用哪些特征
原始负荷曲线维度太高，直接聚类不仅计算慢，而且噪声信息会干扰相似度计算。但特征提取得太狠又会丢失曲线形态信息。我试过几种方案，最终稳定在下面这组特征组合上。
基础统计特征：日最大负荷、日最小负荷、日平均负荷、日负荷率（平均负荷除以最大负荷）、峰谷差。这5个特征反映的是用户用电的整体水平和波动幅度，能区分出“大负荷用户”和“小负荷用户”、“平稳用电”和“剧烈波动”。
时段特征：峰时段电量占比（比如早高峰8点到11点、晚高峰18点到22点）、谷时段电量占比（0点到6点）、夜间负荷占比（22点到次日6点）。这些特征直接关联用户的生活作息，是区分“上班族”“老人”和“夜间活跃用户”的关键。
曲线形态特征：早晚高峰的比值、最大负荷出现时刻、最小负荷出现时刻。这几个特征帮助捕捉曲线的总体形状特点，比如早高峰型和晚高峰型在形状上完全不同。
再用PCA处理一下，把特征维度压到10到20维。实际经验是：PCA不需要保留太多主成分，累计方差贡献率到90%就可以了，后面的成分大多是噪声。这样做既保留了足够信息，又把聚类的计算量降了一大截。
2.3 归一化和典型日曲线计算
归一化是聚类前最重要的一步，没有之一。如果直接拿原始负荷值做距离计算，用户A的日最大负荷是10kW，用户B是2kW，那用电量大的用户天然就会在欧氏距离上占主导，聚类结果基本等于按用电量大小分了层，用曲线形状分群的目的就完全达不到了。
我用的方法是z-score标准化，也就是每个特征减去均值除以标准差，让所有特征处在同一量纲下。这里强调一下：标准化时要用的是全样本的统计量，不能只对某一个用户单独归一化。如果对每个用户单独归一化，这个用户曲线的绝对负荷水平就被抹掉了——“整天大负荷”和“整天小负荷”会被归一化成同一种模式，非常误导。
典型日负荷曲线的计算方法是：对每个用户，把他一个月的日负荷曲线按日期取平均，得到一条96点的平均曲线，代表这个用户在这个月的典型用电模式。之后再去做聚类，就是把每个用户和一条96维的典型曲线对应起来。这样处理的好处是削弱了单日波动的影响，反映的是用户稳定的生活规律，而不是某一天的特殊情况。
3. PSO优化的数学建模与Matlab代码实现
3.1 目标函数设计
PSO优化Kmeans的核心是适应度函数的设计。Kmeans的目标是让同一簇内的样本到聚类中心的总距离最小，那我们就直接用这个目标作为粒子群的适应度。
假设数据集有n个样本，每个样本是m维向量，我们要聚成K类。一个粒子就代表一组聚类中心，它的位置编码成一个K*m维的行向量。解码的时候，把这个行向量前m维取出来作为第1个中心，接下来m维作为第2个中心，以此类推。
适应度函数定义为：
function cost = psoCost(x, data, K, m) % 解码粒子为K个聚类中心 centers = reshape(x, K, m); n = size(data, 1); totalDist = 0; for i = 1:n % 计算样本到所有中心的欧氏距离，取最小值作为该样本到所属中心的距离 dists = sum((centers - data(i, :)).^2, 2); totalDist = totalDist + min(dists); end cost = totalDist; % 总距离越小，粒子适应度越高 end为什么用“总距离”而不是“总距离的均值”？因为聚类时不同K值下样本数和簇大小不一样，均值有时会掩盖某些簇距离很大的问题，总和更直观。实际优化时，PSO的目标就是最小化这个总距离，理论上它等价于Kmeans的收敛目标——最好的结果是，PSO找到的中心让总距离达到一个比随机初始化更小的值，那Kmeans在这个基础上继续迭代就可以收敛到更好的局部最优点。
3.2 PSO参数配置详解
粒子群算法虽然实现简单，但参数调起来还是有门道的。我调试后的一套参数配置如下。
种群规模：20到40之间。粒子数太少，全局搜索能力不足，容易陷入局部最优；粒子数太多，计算量增大。我的数据量是1000个样本、12维特征，K=4，粒子维度是48维，取30个粒子比较合适。
迭代次数：50到100次。对于这个场景，80次足够。迭代曲线会显示总距离在前20次快速下降，后面几十次基本小幅波动，这时候继续增加迭代次数意义不大。
惯性权重w：我从0.9线性递减到0.4。这是粒子群算法里非常经典的做法，前期的较大权重让粒子有较强的探索能力，保持群体多样性；后期权重减小，让粒子收敛到局部精细搜索。刚开始用固定w=0.6跑，结果不理想，后来改成线性递减，收敛效果明显改善。
学习因子c1和c2：都取2.0。c1表示粒子个体经验的权重，c2表示群体共享信息的权重。两者相等可以让粒子在“跟随自己的历史最好位置”和“跟随群体最好位置”之间取得平衡。
速度限制：限制粒子位置在搜索空间范围内，速度最大不超过搜索空间宽度的20%。不限制的话，粒子容易飞出发散。
这些参数在Matlab里的初始化代码：
nParticle = 30; maxIter = 80; wMax = 0.9; wMin = 0.4; c1 = 2.0; c2 = 2.0; dim = K * m; % 粒子维度 % 搜索边界：取数据各维度的最小值和最大值 lb = min(X); ub = max(X); lbAll = repmat(lb, 1, K); ubAll = repmat(ub, 1, K); speedMax = (ubAll - lbAll) * 0.2; % 初始化粒子位置和速度 positions = lbAll + (ubAll - lbAll) .* rand(nParticle, dim); velocities = -speedMax + 2 * speedMax .* rand(nParticle, dim); personalBest = positions; personalBestCost = zeros(nParticle, 1); for i = 1:nParticle personalBestCost(i) = psoCost(positions(i, :), X, K, m); end [globalBestCost, globalBestIdx] = min(personalBestCost); globalBest = personalBest(globalBestIdx, :);3.3 PSO主循环与Kmeans融合
粒子群的主循环在Matlab里写起来非常简洁。每次迭代需要完成的任务是：计算每个粒子的适应度、更新个体历史最优和群体全局最优、更新速度和位置。速度更新公式是经典的那一套：
velocity = w * velocity ... + c1 * rand * (personalBest - position) ... + c2 * rand * (globalBest - position); position = position + velocity;这里有个容易被忽略的地方：边界处理。粒子更新后位置可能超出搜索空间，最常用的做法是把超界的维度拉回到边界值，同时把对应维度的速度设为0，防止它下次又往边界外冲。如果不做这一层处理，粒子发散是早晚的事。
完整主循环：
for iter = 1:maxIter w = wMax - (wMax - wMin) * iter / maxIter; for i = 1:nParticle currentCost = psoCost(positions(i, :), X, K, m); if currentCost < personalBestCost(i) personalBestCost(i) = currentCost; personalBest(i, :) = positions(i, :); end end [currentGlobalCost, currentGlobalIdx] = min(personalBestCost); if currentGlobalCost < globalBestCost globalBestCost = currentGlobalCost; globalBest = personalBest(currentGlobalIdx, :); end for i = 1:nParticle r1 = rand(1, dim); r2 = rand(1, dim); velocities(i, :) = w * velocities(i, :) ... + c1 * r1 .* (personalBest(i, :) - positions(i, :)) ... + c2 * r2 .* (globalBest - positions(i, :)); velocities(i, :) = max(min(velocities(i, :), speedMax), -speedMax); positions(i, :) = positions(i, :) + velocities(i, :); % 边界回拉 positions(i, :) = max(min(positions(i, :), ubAll), lbAll); end endPSO跑完之后，globalBest里就是最优的一组聚类中心。把这组中心解析出来，传给Matlab自带的kmeans函数做初始化。这里用了个小技巧：kmeans函数的'Start'参数可以直接传入初始中心矩阵，不用自己写Kmeans迭代逻辑。
initCenters = reshape(globalBest, K, m); [clusterIdx, clusterCenters, sumD] = kmeans(X, K, 'Start', initCenters, 'MaxIter', 200);这里有个细节值得注意：Kmeans的'MaxIter'设了200次，实际迭代一般十几步就收敛了，因为这个初始中心已经比较优了。我对比过，直接随机初始化时Kmeans经常要迭代50次以上，而且容易收敛到不同结果。
3.4 聚类结果可视化与评估
聚类完成后，不做可视化等于白做。我需要确认两个问题：聚出来的每一类是不是真的有清晰的用电特征？不同类之间的区分度够不够大？
第一张图是典型日负荷曲线图。把每个簇里的所有用户曲线取平均，画96点曲线，叠加在一个图里。如果分得好，你会看到几条曲线形态分得清清楚楚——比如一类是早晚双峰（上班族），一类是平缓少波动的（老人或常驻在家），一类是夜间高耸的（夜间活跃型），还有一类是全天高负荷的（大用电用户）。
第二张图是主成分散点图。用PCA把高维数据降到二维，按聚类结果着色，观察簇在散点图上是否分得开。虽然PCA降维会损失部分信息，但能直观看出分群的紧凑性和分离性。
第三张图是适应度收敛曲线。把PSO迭代过程中globalBestCost的变化画出来，检查是否收敛稳定。如果曲线来回剧烈震荡，说明参数配置有问题，惯性权重和速度限制需要重新调。
指标评估我用两个核心指标。轮廓系数（Silhouette Coefficient）衡量样本自己簇内的紧密度与最近其他簇的分离度，取值从-1到1，越接近1说明分群效果越好。DBI指数（Davies-Bouldin Index）是簇内散度与簇间距离的比值，越小说明簇间离得越开。
% 轮廓系数 s = silhouette(X, clusterIdx); fprintf('轮廓系数均值: %.4f\n', mean(s)); % DBI指数 DBI = daviesbouldin(X, clusterIdx, clusterCenters); % 需要自写函数看轮廓系数时有一个经验：如果某个簇的轮廓系数为负的样本较多，说明这个簇内部结构比较松散，样本分配到其他簇可能更合适。这时候要回头检查K取值是否合理，或者特征是否区分度不够。
4. 实战中踩过的坑与排查方法
4.1 PSO不收敛或者收敛很慢
这是调试过程中最容易遇到也最让人头疼的问题。症状表现是：迭代曲线下降非常缓慢，到了后期还在大幅震荡，或者从一开始就不动弹。
排查思路按优先级来。第一检查粒子初始化范围——如果粒子位置只在很小的子空间里随机生成，那群体多样性不足，怎么迭代都跳不出局部区域。解决办法是把初始化范围扩大到数据取值范围的整个空间。第二检查速度限制——速度上限设得太小，粒子飞得慢，探索效率低；设得太大，粒子容易越过最优区域，来回振荡。我把speedMax设为搜索空间的20%后效果不错。第三检查惯性权重——固定w=0.6时，后期粒子收敛速度慢，群体聚集困难；改用线性递减后才在限定迭代内稳定收敛。
还有一个容易忽略的问题：数据维度过高时，粒子维度也高，适应度函数中的距离计算越是高维，各粒子之间的差异越模糊，PSO的引导作用就越弱。所以特征降维不只是为了省计算时间，它直接影响PSO的搜索有效性。
4.2 聚类数K怎么定
K的取值没有标准答案，但可以用数据驱动的方式辅助判断。常用的是肘部法则：画K与簇内总距离SEE的曲线，看哪里出现拐点。我实际操作时会结合两类指标一起看——SEE曲线的拐点代表“再增加一个聚类数带来的总距离下降明显变缓”；轮廓系数的峰值代表分群效果最好的点。两者取交集，如果指向同一K值，就选它。
还有一个小技巧是从业务出发判断。设K=4时，分出来的四类——双峰型、平缓型、夜间型、全天高负荷型，正好对应四类典型居民用电行为，业务解释性很强。但K=5时，多出来的一类总是在其中两类之间摇摆，业务含义模糊，那就不如用4。聚类不只是数学问题，最终要能解释得通才有价值。
4.3 聚类结果多次运行不一致
直接跑Kmeans时这个问题最严重，换了PSO初始中心后大幅缓解，但也不是完全没有。原因是PSO算法本身是随机搜索，每次运行可能收敛到不同的全局较优解，只是差异远比直接Kmeans小。
如果对可重复性要求很高，两个办法最有效。第一是固定随机种子：
rng(42);所有随机初始化前加上这一句，保证每次运行结果一致，便于复现和交流。第二是多次运行取最优：把整个PSO+Kmeans流程跑20次，记录每次的总距离sumD，选择总距离最小的那次作为最终结果。这相当于在外层再加了一层保险，实际操作下来效果非常稳。
4.4 数据标准化顺序的错误
我刚开始做的时候犯过一个错误：先对原始负荷数据做PCA降维，再对降维后的数据做标准化。这个顺序是错的。PCA本身对数据尺度非常敏感，如果原始特征里有的量纲大有的量纲小，PCA的主成分就会被大量纲特征主导，丢掉的成分反而可能包含有用的形状信息。
正确的做法是：先对所有特征做z-score标准化，再做PCA，最后用降维后的特征做聚类。标准化保证每个维度对距离计算的贡献一致，PCA只负责压缩维度，不改变数据的相对分布，这个顺序千万不能调换。
5. 从聚类结果到业务价值：居民用电行为分析的实际落地
5.1 典型分群结果的业务解读
做完聚类后，最重要的不是把用户标上1、2、3、4的组号，而是把组号和真实的用电行为对应起来。在我分析的那批数据里，最终分出的四类用户呈现出非常清晰的典型特征。
第一类是“典型上班族”，特征是早晚双峰，白天负荷极低，傍晚急剧拉升。这类用户的需求响应策略重点是晚峰时段的可调节负荷，比如空调、热水器、智能家电，可以在晚峰开展削峰。
第二类是“全天平稳型”，负荷曲线波动小，整体维持中等水平，通常是老人、自由职业者或家庭主妇，白天在家用电较多。这类用户适合参与分时电价响应，把部分用电错峰到低谷时段。
第三类是“夜间活跃型”，夜间用电显著高于白天，可能是夜间工作群体或使用电采暖的用户。他们恰恰是填谷策略的理想目标，低谷时段多用电对电网是好事，政策的重点可以放在鼓励蓄热式电暖器、电动汽车夜间充电。
第四类是“高负荷持续型”，全天负荷都在较高水平，可能是电动车充电频繁或运营小型家庭作坊的用户。这类用户单个影响的电网压力很大，需要重点关注。
这种业务解读才是聚类分析的真正价值——不是给用户分堆，而是从分堆结果里读出不同群体的用电行为规律，再让营销部门、调度部门基于这些规律做针对性方案。
5.2 算法的扩展方向
这套PSO+Kmeans的结构本身是通用的，换数据就能迁移到其他场景。比如工商业用户用电行为分析、台区负荷形态辨识、新能源充电行为分群，甚至非电力领域的用户消费行为聚类，逻辑都一样。
后续可以扩展的地方也很多。一是把Kmeans的欧氏距离换成更贴合负荷曲线的距离度量，比如考虑曲线形态的DTW距离，但代价是计算量大幅增加，PSO的适应度评估会变慢。二是引入时间维度，把用户的季节性用电模式也纳入分析，不做典型日曲线而是按周按月聚类。三是把PSO改成多目标优化，同时优化“簇内距离最小”和“簇间距离最大”两个目标，得到帕累托前沿上的一组可选分群方案。
我自己的体会是，这个项目最大的收获不是某一次聚类的准确率提升了多少，而是让我真正理解了“全局搜索+局部优化”这种混合算法框架的强大之处，也体会到了数据预处理环节对最终结果的决定性影响。无论后面用不用PSO，先把数据洗干净、特征选好，你的聚类分析已经成功了一大半。
