做居民用电行为分析最头疼的往往不是算法本身而是数据背后的规律看不到。你拿到的是一堆负荷曲线怎么告诉别人这小区里哪些用户是上班族、哪些是全天在家、哪些可能偷偷搞生产经营聚类就是个好帮手。但真去用Kmeans的时候那个初始质心的选择问题足够让人抓狂——同一个数据集换一次初始值跑出来就是另一种分法。我这段时间正好在Matlab里把粒子群算法和Kmeans拼在一起拿居民用电负荷数据做行为分析整个过程踩了不少坑也把关键细节理清楚了。这篇就把PSO-Kmeans聚类的思路、代码实现和实战经验完整走一遍适合正在做负荷分析、用户画像或者刚接触聚类优化的朋友参考。1. 为什么居民用电分析绕不开聚类优化1.1 负荷曲线背后的用户行为差异先看数据。居民用户的用电行为最直观的载体就是日负荷曲线。把一天24小时或者96个采样点的功率值串联起来就构成一条代表用户当天用电习惯的曲线。正常来说上班族用户的工作日负荷曲线会出现明显的两峰一谷——早高峰可能是7点到9点晚高峰是18点到22点白天和夜晚则相对平缓而老人家庭、待业在家的用户白天的负荷往往比上班族高出一截曲线形态更平稳家里有电动汽车的用户可能晚上会出现一个持续的充电功率平台。不同行为模式之间的差异在负荷曲线上是有迹可循的。但问题在于一个城市或一个台区往往有几千几万个用户靠人工去看曲线分门别类根本不现实。聚类就是用来做这件事的——它能把相似形态的负荷曲线自动归到同一组让有相同用电习惯的用户自然聚到一起。只要聚类方法靠谱分出来的每一类用户你都能倒推出一套对应的行为描述这对接下来的需求响应、分时电价策略、台区负荷预测都特别有价值。1.2 标准Kmeans的天然短板在负荷聚类这个场景里最常用的是Kmeans算法。原理其实特别直白先在样本空间里挑K个点当初始质心然后把每个样本分配给离它最近的质心分完以后再重新计算每个簇的中心点反复迭代直到结果稳定。整个过程就是分配—更新—再分配—再更新很符合直觉代码也简单。但它有个致命问题——Kmeans对初始质心的选择极度敏感。初始质心选偏了迭代多少次都可能停留在某个局部最优解上。比如有两类用户一类是白天用电一类是晚上用电如果你的初始质心都落在白天那堆数据里晚上那一类很可能被硬生生拆散最终聚类结果从业务角度怎么解释都不合理。而且Kmeans迭代过程中一旦某个簇被分配为空算法还会出现质心失效的异常情况。实际处理负荷数据时样本量大、曲线波动多Kmeans跑出来的结果经常不稳定同一份数据跑十次能有七八种分法。这就是我要引入粒子群算法的直接原因。粒子群优化算法PSO是一种全局寻优方法它的思路是模拟鸟群觅食——每个粒子代表一个候选解靠个体经验个体最优和群体经验全局最优不断调整自己的位置逐步逼近全局最优解。把PSO和Kmeans结合简单说就是用PSO先把Kmeans的初始质心这个老大难问题解决掉让聚类从一个好的起点开始跑。这样既保留了Kmeans计算快的优点又显著降低了落入局部最优的概率。下面我把这套组合方案的原理和代码一步步拆开讲。2. PSO和Kmeans是怎么配合的2.1 先从Kmeans的目标说起要理解PSO-Kmeans的配合逻辑得先把Kmeans在做什么看透。Kmeans本质上是求解一个最小化问题把n个样本分到K个簇中让所有样本到所属簇质心的距离平方和也就是簇内误差平方和SSE最小。公式写出来是[ SSE \sum_{i1}^{K}\sum_{x \in C_i} |x - \mu_i|^2 ]其中(C_i)是第i个簇(\mu_i)是这个簇的质心。聚类结果好不好直接看SSE——SSE越小说明簇内样本越紧凑同类用户之间的相似度越高。但Kmeans采用的是一种贪心式的交替优化先固定质心分配样本再固定分配更新质心。这种方式求解速度快却很依赖初始质心给得怎么样。初始质心离全局最优解太远交替优化就可能收敛到SSE较大的局部最优解。怎么跳出这个坑我在项目里采用的思路是把Kmeans的初始质心当作粒子群优化算法中的决策变量用PSO去全局搜索一组好的质心位置再把搜到的结果作为Kmeans的起点。换句话说Kmeans负责局部精修PSO负责全局寻优两者分工合作。这里有一个技术路线选择问题。有些人会把PSO直接作为聚类工具来用让每个样本以一定概率归属某个簇最后按隶属度划分。这个方案在高维负荷数据上计算量非常大而且解释性不如Kmeans清晰。我更推荐的是把PSO定位成质心初始化优化器后面照常跟Kmeans迭代。实测下来这种方式既快又稳结果还好用业务语言解释。2.2 粒子编码方式与适应度函数设计要把PSO用于优化Kmeans的初始质心第一个要解决的是粒子怎么编码的问题。假设负荷数据经过特征处理后每个样本是一个D维向量最简单的做法就是24维对应24小时负荷值我实验里也用过96维对应96个采样点聚类数设为K。那么一组完整的初始质心就是K个D维向量把它们按顺序拼接成一个长向量这个长向量就是一个粒子的位置。粒子的维度就是(K \times D)。举个例子K4D24粒子维度就是96。在Matlab里我习惯用矩阵来组织蜂群——一个粒子用一个(K \times D)的矩阵表示整个粒子群用一个三维数组存储这样在计算距离时可以避免频繁的reshape操作。适应度函数的设计是整个算法的灵魂。PSO的迭代方向完全靠适应度值牵引。我这里使用的适应度函数就是Kmeans的SSE[ fitness(x) \sum_{j1}^{K}\sum_{x_n \in C_j} |x_n - z_j|^2 ]其中(z_j)是粒子的第j个质心位置。计算时先按每个样本到各质心的欧氏距离做最近邻分配再累加出SSE。适应度值越小代表该粒子对应的质心组合越好。有人会问PSO时代里的好到底是全局最优还是局部最优这正是PSO的价值所在——粒子群中的每个粒子都在自己的位置附近搜索同时向全局最优粒子靠拢这种信息共享机制让种群不容易卡死在单个局部区域。配合惯性权重和学习因子PSO能在搜索前期保持较强的全局探索能力后期逐渐收敛到精细区域。这比随机撒点选初始质心要靠谱得多。2.3 算法流程梳理我实际跑通的PSO-Kmeans完整流程如下读取并预处理负荷数据缺失值处理、归一化。确定聚类数K用轮廓系数或肘部法则辅助判断后面细说。初始化粒子群体。每个粒子的位置为K个随机的样本点速度为全零或小随机数。对每个粒子计算适应度SSE更新个体最优pbest和群体最优gbest。按标准PSO公式更新粒子速度和位置。速度更新公式为 [ v_{t1} w \cdot v_t c_1 r_1 (pbest - x_t) c_2 r_2 (gbest - x_t) ] 位置更新为简单累加。这里w是惯性权重c1和c2是学习因子r1和r2是0到1之间的随机数。检查是否达到最大迭代次数否则返回第4步。将gbest还原为(K \times D)的质心矩阵作为Kmeans的初始质心。执行标准Kmeans迭代分配样本、更新质心直到收敛。输出聚类标签、质心、SSE并做可视化。整个流程里PSO阶段其实相当于在做全局热身Kmeans阶段在做局部冲刺。我用这个方案对比过纯Kmeans在典型居民负荷数据上SSE能降低15%到25%而且多次运行的结果稳定性明显提升。3. Matlab代码实现与参数配置3.1 数据准备与特征构建代码实现上我建议先把特征工程做成独立脚本不要把数据处理和聚类算法混在一起。原始用电数据通常是这样的每15分钟一个采样点一天96个点连续若干天。但如果直接用96维做聚类维度较高PSO粒子搜索空间的体积会指数增长不仅慢而且效果不一定好。我的做法是折中做日负荷曲线特征压缩。常用的压缩方式有几种我列个表对比一下特征方案维度优点缺点24小时均值负荷24直观、计算快能保留峰谷形态丢失了日内变化细节96点原始负荷96信息完整维度高PSO粒子维度过大易过拟合峰谷特征峰时负荷、谷时负荷、峰谷差、日用电量等5-8业务解释性强维度低需要按当地峰谷时段定义有一定主观性统计特征均值、方差、峰度、偏度、最大负荷时间5-8压缩程度高形态信息流失多我在项目里最终选的是24小时均值负荷几个统计特征的组合总维度约28。折中的原因有两个一是24小时曲线能让聚类结果直接画图解释生成工作族居家型这种标签二是维度控制在30以内PSO搜索效率高很多。作为补充我也做了96维的对照实验后面在问题排查部分会讲这个方案踩了什么坑。数据清洗这一步很关键。居民负荷数据里常见的问题是采集终端偶尔掉线导致整天数据是0或者个别时段出现异常尖峰。我的处理规则是连续3小时以上全为0的用户直接剔除非零时段中超过99.5%分位的数值视为异常尖峰用前后时刻的均值替换。这些规则比单纯用是否大于某阈值判断更鲁棒。归一化也要特别注意。如果不做归一化用电量大的用户比如冬夏开空调日电量几十度甚至上百度会在欧氏距离计算中占据绝对主导聚类结果基本就变成了按用电量分等级而不是按行为模式分类。我的做法是按特征列做Z-score标准化也就是每列减去均值再除以标准差这样每个特征对距离的贡献平等。在Matlab里一行代码就能搞定data_norm zscore(data_raw);处理完以后记得保存一份标准化参数后面做新用户分类或者画原尺度曲线时要用。3.2 PSO-Kmeans主程序编写主程序我分了三个函数块粒子初始化、适应度计算、PSO迭代主循环。这种模块化写法方便调试也便于替换不同的适应度函数或数据集。先看粒子初始化% 输入data为标准化后的样本矩阵(nxD)K为聚类数N为种群规模 % 输出particle为(N, K, D)的三维数组 n size(data, 1); D size(data, 2); particle zeros(N, K, D); velocity zeros(N, K, D); for i 1:N idx randperm(n, K); % 随机选K个样本作为初始质心 particle(i, :, :) data(idx, :); velocity(i, :, :) 0.02 * randn(K, D); end初始化方式选择随机取样本点而不是在整个搜索空间随机撒点。原因是负荷数据做完Z-score标准化后虽然有少数离群点但绝大多数样本都集中在可行区域内。从样本中选初始质心相当于一开始就没有偏离合理区域能明显加快收敛。这个细节我建议一定保留。适应度函数我单独写核心逻辑如下function fitness calcFitness(data, particle_i, K) n size(data, 1); distMat zeros(n, K); for j 1:K centroid squeeze(particle_i(j, :)); diff data - centroid; % n x D distMat(:, j) sqrt(sum(diff.^2, 2)); end [~, assign] min(distMat, [], 2); fitness 0; for j 1:K clusterData data(assign j, :); if ~isempty(clusterData) centroid mean(clusterData, 1); fitness fitness sum(sum((clusterData - centroid).^2, 2)); end end end注意这里我在适应度计算中不是用粒子自带质心算SSE而是按分配结果重新计算实际质心再算SSE。为什么不直接用粒子里的质心因为粒子在PSO迭代中可能移动到远离任何样本的位置用空簇质心算距离会产生虚低的SSE误导搜索方向。重新计算簇质心相当于做了局部投影适应度值更真实。这个细节是我调试过程中对比了几种方案后确定的效果确实更稳。主迭代循环采用标准的PSO公式惯性权重w随迭代次数线性递减maxIter 50; N 30; K 4; c1 1.5; c2 1.5; wMax 0.9; wMin 0.4; pbestScore inf(N, 1); pbestParticle particle; gbestScore inf; gbestParticle squeeze(particle(1, :, :)); for t 1:maxIter w wMax - (wMax - wMin) * t / maxIter; for i 1:N fitness calcFitness(data, squeeze(particle(i, :, :)), K); if fitness pbestScore(i) pbestScore(i) fitness; pbestParticle(i, :, :) particle(i, :, :); end if fitness gbestScore gbestScore fitness; gbestParticle squeeze(particle(i, :, :)); end end for i 1:N r1 rand(K, D); r2 rand(K, D); velocity(i, :, :) w * velocity(i, :, :) ... c1 * r1 .* (squeeze(pbestParticle(i, :, :)) - squeeze(particle(i, :, :))) ... c2 * r2 .* (gbestParticle - squeeze(particle(i, :, :))); particle(i, :, :) particle(i, :, :) velocity(i, :, :); end end最后把gbestParticle作为初始质心送给Kmeans[clusterIdx, centroid] kmeans(data, K, Start, gbestParticle, MaxIter, 1000);如果Matlab版本较老不支持Start参数直接传入矩阵可以先调用类的静态方法设置选项再执行聚类或者自己手写10-20轮Kmeans迭代。老版本其实也完全可以用我后面遇到过一次版本兼容问题在常见问题部分会展开说明。3.3 关键参数的选择依据与调试建议PSO-Kmeans涉及到的参数不少我把我实测下来比较合适的配置整理一下。种群规模N我建议取20到40之间。太小了全局搜索能力不足太大了计算量明显上升。负荷曲线的样本数通常在几千到几万之间每次适应度计算都要遍历所有样本做距离计算N取30不算大但加上50次迭代在几千样本量下Matlab要跑几十秒可以接受。如果样本量超过5万建议先把训练集采样到1万规模做粒子搜索再用跑出来的质心初始化全量Kmeans。最大迭代次数maxIter50次通常够了。我在调试时观察过适应度收敛曲线大约在30次以后下降曲线就趋于平缓50次属于留有余量。如果追求速度25到30次也能得到差不多的结果差别在2%以内。但首次实验我建议还是跑到50次先把算法的稳定基线摸清楚。惯性权重w采用0.9到0.4线性递减。前期w大粒子飞得快、探索范围广不容易陷进局部最优后期w小粒子精细琢磨加速收敛。这个区间是粒子群算法的经典经验值实测在聚类问题上效果稳定。学习因子c1和c2取1.5是比较均衡的组合。也有文献推荐c1c22我试过收敛快一些但偶尔会跳过好的质心区域。1.5加上0.9到0.4的惯性权重搭配探索和开发平衡得更舒服。如果你发现结果波动大可以尝试把c1降到1.2、c2提到1.8增强向群体最优靠拢的趋势。聚类数K用轮廓系数辅助判断。轮廓系数综合考虑了簇内紧密度和簇间分离度取值范围-1到1越大代表聚类效果越好。我在项目里对K2到K8分别跑PSO-Kmeans计算每个K下的平均轮廓系数选峰值对应的K。实际业务上K取4或5比较常见这样每一类用户都有足够明确的画像不会分得过细而失去解释力。4. 实验效果分析与聚类结果解读4.1 与标准Kmeans的对比实验我拿来验证的数据是某市一个台区3000户居民用户30天的用电记录按前文方法清洗和特征化后得到3000×28的特征矩阵聚类目标K4PSO种群取30迭代50次。为了控制变量标准Kmeans我用Matlab自带的kmeans函数跑100次随机初始化取SSE最小的一次作为参照这种多次随机取最优本身就是实践中应对Kmeans不稳定的常见手段但计算开销远高于PSO辅助。最终实验数据如下表方案平均SSE最优SSE波动范围SSE单次运行耗时标准Kmeans单次1846.71752.3160.40.8秒标准Kmeans100次取最优1635.21635.2024秒PSO-Kmeans单次1658.11641.533.218秒PSO-Kmeans3次取最优1642.01641.53.254秒几个结论很直观。PSO-Kmeans单次结果明显优于Kmeans单次SSE从1846.7降到1658.1下降了大约10.2%即使对比Kmeans跑100次取最优的1635.2PSO-Kmeans的最优SSE 1641.5也非常接近差了不到0.4%。更关键的是稳定性——PSO-Kmeans三次运行的最优与最差只差33.2几乎都在同一水平线上这说明算法已经不太受随机初始化的影响而标准Kmeans单次运行的波动范围高达160以上这在工程上非常致命。当然PSO-Kmeans也不是免费的午餐18秒的处理时间比标准Kmeans单次0.8秒慢得多。但对离线用户画像分析这种场景18秒完全可接受。4.2 聚类结果如何映射到用电行为聚类跑完只是第一步更重要的工作是把每一类用户的行为模式描述出来。我是这样做的拿到聚类标签后把原始负荷数据未标准化按类分组计算每类用户的平均24小时负荷曲线然后结合日用电量、峰谷比等业务指标做解读。在我的实验里K4时的四类用户画像如下第一类工作日早、晚双峰特别突出白天负荷很低午间有小幅回落周末曲线相对平缓。结合日用电量处于中低水平可以判定为典型的上班族家庭工作日只有早晚在家用电。第二类白天负荷较高曲线全天相对平稳夜晚略降但不会降到很低日用电量处于中上水平。这是全天居家型用户可能是老人、家庭主妇或自由职业者。第三类夜间和凌晨负荷异常偏高白天反而较低日用电量也比较大。结合当地电价政策这类用户很可能是有意将洗衣机、热水器等大功率设备挪到夜间使用甚至可能有电动汽车充电行为。第四类整体负荷水平低曲线平缓无峰长时间维持很小的用电功率。这种通常是空心户或者出租率较高的房屋用电行为不活跃。每类用户对应的策略建议也不一样第一类适合宣传分时电价引导削峰填谷第二类可以推荐节能设备第三类可以作为需求响应的重点对象第四类则需要在台区管理上排查是否有空置房或者表计异常。这些业务层面的延伸才是分析工作真正产生价值的地方。4.3 可视化技巧如何把聚类结果画得让业务方看懂聚类结果可视化我踩过不少坑。最开始我直接用plot画所有用户的原始曲线3000条线叠在一起密密麻麻根本看不出差异。后来改成每个类画一条平均曲线标准差带效果立刻不一样。Matlab里用fill可以画带meanCurve mean(clusterData, 1); stdCurve std(clusterData, 1); t 1:24; fill([t fliplr(t)], [meanCurvestdCurve fliplr(meanCurve-stdCurve)], ... [0.9 0.9 0.9], FaceAlpha, 0.4, EdgeColor, none); hold on; plot(t, meanCurve, LineWidth, 2);标准差带能够直观表达这一类用户内部的波动程度。如果某类的带很窄说明这类用户的负荷形态高度一致聚类可信度高带很宽则说明这一类内部还存在细分可以考虑是否增加K值。另外一个可视化技巧是降维散点图。高维特征矩阵不好直接展示可以用t-SNE或者PCA降到2维再按聚类标签着色。不过我要提醒一句降维后再看聚类是否分得开只能作为辅助参考因为降维过程会扭曲真实距离关系。业务汇报时这东西很好看内部验证时别太当真。5. 常见问题与排查技巧实录5.1 粒子维度爆炸和计算速度慢怎么办我在96维特征上尝试过直接跑PSO-Kmeans粒子维度是(K \times 96)K取4就是384维。粒子群优化在这么高的维度上进行搜索效果非常差——适应度收敛慢、粒子群容易散开、结果还不稳定。因为高维空间里距离度量变得稀疏随机初始化的粒子互相之间差异很小PSO很难通过对比分辨哪个方向更好。解决思路有两个。第一是在特征层面降维比如用24小时均值替代96点数据或者先用PCA把特征压到15到20维再做聚类。第二是改变PSO的搜索策略比如将速度初始化设置为0限制粒子的搜索半径但这样又会牺牲全局搜索能力。我的建议是优先做特征降维因为居民负荷数据本身的冗余度很高96个采样点之间存在很强的时序相关性强行保留全部维度得不偿失。计算速度问题还有另一层来源适应度函数里频繁的矩阵运算。如果循环写的效率低几千样本都够让Matlab卡上几分钟。我把计算距离的代码从for循环改成矩阵广播后原来45秒一次迭代缩到3秒左右。Matlab效率的关键就是不要让循环套循环多用维度广播和矩阵运算如果还想更快可以把calcFitness写成mex函数或者用parfor并行计算粒子群中不同粒子的适应度。5.2 陷入局部最优的判断与处理有一种情况PSO迭代结束后gbest对应的质心组其实还不是理想解Kmeans再迭代也跳不出来。怎么判断我会把PSO-Kmeans的SSE和多次随机初始化的Kmeans最优SSE做对比如果前者显著大于后者基本可以断定PSO阶段早收敛了。处理办法有这么几种。一是检查粒子群初始化如果初始粒子全都挤在样本集中的区域多样性不够PSO很容易早熟。初始化时除了随机采样样本点我还会刻意加几个远离中心的点。二是增大惯性权重或者调节学习因子如果w从0.9降到0.4太快个体经验权重过大可以在实验中把wMax提到1.0wMin提到0.5让粒子飞得更激进一点。三是重启策略如果一个粒子连续N代都没有改进自己的pbest给它重新初始化到随机位置这是个简单但很有效的办法。我再强调一次PSO-Kmeans不是银弹它只能显著降低落入局部最优的概率不能完全消除。所以在项目落地时我通常跑3次PSO-Kmeans取SSE最小的那次。由于单次已经很稳定3次取最优带来的额外收益也有限更多是买个心理保险。5.3 K值怎么选最合理选择K值最常见的是肘部法则画SSE随K变化的折线图找那个拐点。但实际数据里肘部往往不明显SSE下降曲线保持平滑你很难说出3和4哪个是肘。我用轮廓系数配合业务可解释性一起判断。轮廓系数对第i个样本的定义是[ s_i \frac{b_i - a_i}{\max(a_i, b_i)} ]其中(a_i)是样本i与同簇其他样本的平均距离(b_i)是样本i与最近其他簇的平均距离。把全部样本的轮廓系数平均就是总体轮廓系数。我一般要求总体轮廓系数大于等于0.5如果某个K下只有0.3说明簇内不够紧凑或者簇间分得不清楚这个K值基本不可用。但我也要说业务可解释性有时候比数值指标更关键。比如K5时轮廓系数最高但其中有一类用户曲线形态和另一类非常接近业务上完全无法区分和应对那K5就没有实际意义。我的习惯是先选2到3个候选K轮廓系数比较高的然后把这几个K下的聚类结果拿给业务同事看问哪一版最容易讲故事通常答案很明确。5.4 版本兼容和Matlab环境的坑我在实验过程中遇到过一次运行环境导致的怪问题在Matlab R2021b上能正常运行的脚本换到老版本后kmeans的Start参数传矩阵就报错。Matlab每个版本对聚类函数输入参数的校验机制不一样如果公司或实验室的Matlab版本不统一建议不要依赖版本较新的参数特性。我的做法是手写一个20轮的Kmeans精修函数替代内置的kmeans代码不超过30行却能在所有版本上稳定运行。核心逻辑就是循环分配样本—更新质心和我们第一部分讲的Kmeans原理完全一致。另外如果你跟我一样被工程化逼得没有正版授权也可以考虑用GNU Octave代替Matlab写这个流程。Octave对大部分数值计算和矩阵运算的支持都很好PSO-Kmeans这种以矩阵运算为主的代码迁移成本很低。不过Octave的kmeans函数不是内置的需要自己手写用来替代内置函数时正好省了上面的兼容性问题。5.5 数据质量细节这些坑会影响聚类结论最后分享几个和算法无关但直接影响结论的数据细节。第一归一化必须在缺失值处理之后做否则Z-score会把缺失值当成0参与均值计算扭曲特征分布。第二聚类的输入应该是行为特征不应该直接放日期、用户编号、台区编号这些标识性变量。第三如果用户数据的天数不一致有的用户只有15天记录有的有30天建议先按用户求平均再做聚类否则天数少的用户会被当成异常样本。第四季节因素要重视——冬季和夏季的负荷曲线形态差异很大如果你直接拿一整年数据混在一起聚类得到的分群往往是季节分群而非行为分群。我的做法是按季节分别建模型然后在业务层面对比同一用户的季节归属变化这样既能识别行为差异又能捕捉季节性规律变化。这套组合方案跑下来我最大的体会是算法层面没有太多高大上的东西PSO-Kmeans本质上是把一个简单而顽固的问题——初始质心敏感——用群智能算法解决掉了。居民用电行为分析的价值也不在于把轮廓系数从0.55提高到0.6而在于每一类用户分出来以后你能针对性地做点什么。最后再分享一个小技巧给准备落地的朋友就算聚类结果已经稳定也别直接信任数据去抽查10个用户的原始负荷曲线和聚类标签是否匹配。光看平均曲线会骗人单条曲线才暴露真相。这个步骤花不了十分钟却能避免向业务方汇报时被一句我看这明显不是一类用户问得哑口无言。 PSO-Kmeans聚类优化：居民用电负荷曲线用户行为分析实战
做居民用电行为分析最头疼的往往不是算法本身而是数据背后的规律看不到。你拿到的是一堆负荷曲线怎么告诉别人这小区里哪些用户是上班族、哪些是全天在家、哪些可能偷偷搞生产经营聚类就是个好帮手。但真去用Kmeans的时候那个初始质心的选择问题足够让人抓狂——同一个数据集换一次初始值跑出来就是另一种分法。我这段时间正好在Matlab里把粒子群算法和Kmeans拼在一起拿居民用电负荷数据做行为分析整个过程踩了不少坑也把关键细节理清楚了。这篇就把PSO-Kmeans聚类的思路、代码实现和实战经验完整走一遍适合正在做负荷分析、用户画像或者刚接触聚类优化的朋友参考。1. 为什么居民用电分析绕不开聚类优化1.1 负荷曲线背后的用户行为差异先看数据。居民用户的用电行为最直观的载体就是日负荷曲线。把一天24小时或者96个采样点的功率值串联起来就构成一条代表用户当天用电习惯的曲线。正常来说上班族用户的工作日负荷曲线会出现明显的两峰一谷——早高峰可能是7点到9点晚高峰是18点到22点白天和夜晚则相对平缓而老人家庭、待业在家的用户白天的负荷往往比上班族高出一截曲线形态更平稳家里有电动汽车的用户可能晚上会出现一个持续的充电功率平台。不同行为模式之间的差异在负荷曲线上是有迹可循的。但问题在于一个城市或一个台区往往有几千几万个用户靠人工去看曲线分门别类根本不现实。聚类就是用来做这件事的——它能把相似形态的负荷曲线自动归到同一组让有相同用电习惯的用户自然聚到一起。只要聚类方法靠谱分出来的每一类用户你都能倒推出一套对应的行为描述这对接下来的需求响应、分时电价策略、台区负荷预测都特别有价值。1.2 标准Kmeans的天然短板在负荷聚类这个场景里最常用的是Kmeans算法。原理其实特别直白先在样本空间里挑K个点当初始质心然后把每个样本分配给离它最近的质心分完以后再重新计算每个簇的中心点反复迭代直到结果稳定。整个过程就是分配—更新—再分配—再更新很符合直觉代码也简单。但它有个致命问题——Kmeans对初始质心的选择极度敏感。初始质心选偏了迭代多少次都可能停留在某个局部最优解上。比如有两类用户一类是白天用电一类是晚上用电如果你的初始质心都落在白天那堆数据里晚上那一类很可能被硬生生拆散最终聚类结果从业务角度怎么解释都不合理。而且Kmeans迭代过程中一旦某个簇被分配为空算法还会出现质心失效的异常情况。实际处理负荷数据时样本量大、曲线波动多Kmeans跑出来的结果经常不稳定同一份数据跑十次能有七八种分法。这就是我要引入粒子群算法的直接原因。粒子群优化算法PSO是一种全局寻优方法它的思路是模拟鸟群觅食——每个粒子代表一个候选解靠个体经验个体最优和群体经验全局最优不断调整自己的位置逐步逼近全局最优解。把PSO和Kmeans结合简单说就是用PSO先把Kmeans的初始质心这个老大难问题解决掉让聚类从一个好的起点开始跑。这样既保留了Kmeans计算快的优点又显著降低了落入局部最优的概率。下面我把这套组合方案的原理和代码一步步拆开讲。2. PSO和Kmeans是怎么配合的2.1 先从Kmeans的目标说起要理解PSO-Kmeans的配合逻辑得先把Kmeans在做什么看透。Kmeans本质上是求解一个最小化问题把n个样本分到K个簇中让所有样本到所属簇质心的距离平方和也就是簇内误差平方和SSE最小。公式写出来是[ SSE \sum_{i1}^{K}\sum_{x \in C_i} |x - \mu_i|^2 ]其中(C_i)是第i个簇(\mu_i)是这个簇的质心。聚类结果好不好直接看SSE——SSE越小说明簇内样本越紧凑同类用户之间的相似度越高。但Kmeans采用的是一种贪心式的交替优化先固定质心分配样本再固定分配更新质心。这种方式求解速度快却很依赖初始质心给得怎么样。初始质心离全局最优解太远交替优化就可能收敛到SSE较大的局部最优解。怎么跳出这个坑我在项目里采用的思路是把Kmeans的初始质心当作粒子群优化算法中的决策变量用PSO去全局搜索一组好的质心位置再把搜到的结果作为Kmeans的起点。换句话说Kmeans负责局部精修PSO负责全局寻优两者分工合作。这里有一个技术路线选择问题。有些人会把PSO直接作为聚类工具来用让每个样本以一定概率归属某个簇最后按隶属度划分。这个方案在高维负荷数据上计算量非常大而且解释性不如Kmeans清晰。我更推荐的是把PSO定位成质心初始化优化器后面照常跟Kmeans迭代。实测下来这种方式既快又稳结果还好用业务语言解释。2.2 粒子编码方式与适应度函数设计要把PSO用于优化Kmeans的初始质心第一个要解决的是粒子怎么编码的问题。假设负荷数据经过特征处理后每个样本是一个D维向量最简单的做法就是24维对应24小时负荷值我实验里也用过96维对应96个采样点聚类数设为K。那么一组完整的初始质心就是K个D维向量把它们按顺序拼接成一个长向量这个长向量就是一个粒子的位置。粒子的维度就是(K \times D)。举个例子K4D24粒子维度就是96。在Matlab里我习惯用矩阵来组织蜂群——一个粒子用一个(K \times D)的矩阵表示整个粒子群用一个三维数组存储这样在计算距离时可以避免频繁的reshape操作。适应度函数的设计是整个算法的灵魂。PSO的迭代方向完全靠适应度值牵引。我这里使用的适应度函数就是Kmeans的SSE[ fitness(x) \sum_{j1}^{K}\sum_{x_n \in C_j} |x_n - z_j|^2 ]其中(z_j)是粒子的第j个质心位置。计算时先按每个样本到各质心的欧氏距离做最近邻分配再累加出SSE。适应度值越小代表该粒子对应的质心组合越好。有人会问PSO时代里的好到底是全局最优还是局部最优这正是PSO的价值所在——粒子群中的每个粒子都在自己的位置附近搜索同时向全局最优粒子靠拢这种信息共享机制让种群不容易卡死在单个局部区域。配合惯性权重和学习因子PSO能在搜索前期保持较强的全局探索能力后期逐渐收敛到精细区域。这比随机撒点选初始质心要靠谱得多。2.3 算法流程梳理我实际跑通的PSO-Kmeans完整流程如下读取并预处理负荷数据缺失值处理、归一化。确定聚类数K用轮廓系数或肘部法则辅助判断后面细说。初始化粒子群体。每个粒子的位置为K个随机的样本点速度为全零或小随机数。对每个粒子计算适应度SSE更新个体最优pbest和群体最优gbest。按标准PSO公式更新粒子速度和位置。速度更新公式为 [ v_{t1} w \cdot v_t c_1 r_1 (pbest - x_t) c_2 r_2 (gbest - x_t) ] 位置更新为简单累加。这里w是惯性权重c1和c2是学习因子r1和r2是0到1之间的随机数。检查是否达到最大迭代次数否则返回第4步。将gbest还原为(K \times D)的质心矩阵作为Kmeans的初始质心。执行标准Kmeans迭代分配样本、更新质心直到收敛。输出聚类标签、质心、SSE并做可视化。整个流程里PSO阶段其实相当于在做全局热身Kmeans阶段在做局部冲刺。我用这个方案对比过纯Kmeans在典型居民负荷数据上SSE能降低15%到25%而且多次运行的结果稳定性明显提升。3. Matlab代码实现与参数配置3.1 数据准备与特征构建代码实现上我建议先把特征工程做成独立脚本不要把数据处理和聚类算法混在一起。原始用电数据通常是这样的每15分钟一个采样点一天96个点连续若干天。但如果直接用96维做聚类维度较高PSO粒子搜索空间的体积会指数增长不仅慢而且效果不一定好。我的做法是折中做日负荷曲线特征压缩。常用的压缩方式有几种我列个表对比一下特征方案维度优点缺点24小时均值负荷24直观、计算快能保留峰谷形态丢失了日内变化细节96点原始负荷96信息完整维度高PSO粒子维度过大易过拟合峰谷特征峰时负荷、谷时负荷、峰谷差、日用电量等5-8业务解释性强维度低需要按当地峰谷时段定义有一定主观性统计特征均值、方差、峰度、偏度、最大负荷时间5-8压缩程度高形态信息流失多我在项目里最终选的是24小时均值负荷几个统计特征的组合总维度约28。折中的原因有两个一是24小时曲线能让聚类结果直接画图解释生成工作族居家型这种标签二是维度控制在30以内PSO搜索效率高很多。作为补充我也做了96维的对照实验后面在问题排查部分会讲这个方案踩了什么坑。数据清洗这一步很关键。居民负荷数据里常见的问题是采集终端偶尔掉线导致整天数据是0或者个别时段出现异常尖峰。我的处理规则是连续3小时以上全为0的用户直接剔除非零时段中超过99.5%分位的数值视为异常尖峰用前后时刻的均值替换。这些规则比单纯用是否大于某阈值判断更鲁棒。归一化也要特别注意。如果不做归一化用电量大的用户比如冬夏开空调日电量几十度甚至上百度会在欧氏距离计算中占据绝对主导聚类结果基本就变成了按用电量分等级而不是按行为模式分类。我的做法是按特征列做Z-score标准化也就是每列减去均值再除以标准差这样每个特征对距离的贡献平等。在Matlab里一行代码就能搞定data_norm zscore(data_raw);处理完以后记得保存一份标准化参数后面做新用户分类或者画原尺度曲线时要用。3.2 PSO-Kmeans主程序编写主程序我分了三个函数块粒子初始化、适应度计算、PSO迭代主循环。这种模块化写法方便调试也便于替换不同的适应度函数或数据集。先看粒子初始化% 输入data为标准化后的样本矩阵(nxD)K为聚类数N为种群规模 % 输出particle为(N, K, D)的三维数组 n size(data, 1); D size(data, 2); particle zeros(N, K, D); velocity zeros(N, K, D); for i 1:N idx randperm(n, K); % 随机选K个样本作为初始质心 particle(i, :, :) data(idx, :); velocity(i, :, :) 0.02 * randn(K, D); end初始化方式选择随机取样本点而不是在整个搜索空间随机撒点。原因是负荷数据做完Z-score标准化后虽然有少数离群点但绝大多数样本都集中在可行区域内。从样本中选初始质心相当于一开始就没有偏离合理区域能明显加快收敛。这个细节我建议一定保留。适应度函数我单独写核心逻辑如下function fitness calcFitness(data, particle_i, K) n size(data, 1); distMat zeros(n, K); for j 1:K centroid squeeze(particle_i(j, :)); diff data - centroid; % n x D distMat(:, j) sqrt(sum(diff.^2, 2)); end [~, assign] min(distMat, [], 2); fitness 0; for j 1:K clusterData data(assign j, :); if ~isempty(clusterData) centroid mean(clusterData, 1); fitness fitness sum(sum((clusterData - centroid).^2, 2)); end end end注意这里我在适应度计算中不是用粒子自带质心算SSE而是按分配结果重新计算实际质心再算SSE。为什么不直接用粒子里的质心因为粒子在PSO迭代中可能移动到远离任何样本的位置用空簇质心算距离会产生虚低的SSE误导搜索方向。重新计算簇质心相当于做了局部投影适应度值更真实。这个细节是我调试过程中对比了几种方案后确定的效果确实更稳。主迭代循环采用标准的PSO公式惯性权重w随迭代次数线性递减maxIter 50; N 30; K 4; c1 1.5; c2 1.5; wMax 0.9; wMin 0.4; pbestScore inf(N, 1); pbestParticle particle; gbestScore inf; gbestParticle squeeze(particle(1, :, :)); for t 1:maxIter w wMax - (wMax - wMin) * t / maxIter; for i 1:N fitness calcFitness(data, squeeze(particle(i, :, :)), K); if fitness pbestScore(i) pbestScore(i) fitness; pbestParticle(i, :, :) particle(i, :, :); end if fitness gbestScore gbestScore fitness; gbestParticle squeeze(particle(i, :, :)); end end for i 1:N r1 rand(K, D); r2 rand(K, D); velocity(i, :, :) w * velocity(i, :, :) ... c1 * r1 .* (squeeze(pbestParticle(i, :, :)) - squeeze(particle(i, :, :))) ... c2 * r2 .* (gbestParticle - squeeze(particle(i, :, :))); particle(i, :, :) particle(i, :, :) velocity(i, :, :); end end最后把gbestParticle作为初始质心送给Kmeans[clusterIdx, centroid] kmeans(data, K, Start, gbestParticle, MaxIter, 1000);如果Matlab版本较老不支持Start参数直接传入矩阵可以先调用类的静态方法设置选项再执行聚类或者自己手写10-20轮Kmeans迭代。老版本其实也完全可以用我后面遇到过一次版本兼容问题在常见问题部分会展开说明。3.3 关键参数的选择依据与调试建议PSO-Kmeans涉及到的参数不少我把我实测下来比较合适的配置整理一下。种群规模N我建议取20到40之间。太小了全局搜索能力不足太大了计算量明显上升。负荷曲线的样本数通常在几千到几万之间每次适应度计算都要遍历所有样本做距离计算N取30不算大但加上50次迭代在几千样本量下Matlab要跑几十秒可以接受。如果样本量超过5万建议先把训练集采样到1万规模做粒子搜索再用跑出来的质心初始化全量Kmeans。最大迭代次数maxIter50次通常够了。我在调试时观察过适应度收敛曲线大约在30次以后下降曲线就趋于平缓50次属于留有余量。如果追求速度25到30次也能得到差不多的结果差别在2%以内。但首次实验我建议还是跑到50次先把算法的稳定基线摸清楚。惯性权重w采用0.9到0.4线性递减。前期w大粒子飞得快、探索范围广不容易陷进局部最优后期w小粒子精细琢磨加速收敛。这个区间是粒子群算法的经典经验值实测在聚类问题上效果稳定。学习因子c1和c2取1.5是比较均衡的组合。也有文献推荐c1c22我试过收敛快一些但偶尔会跳过好的质心区域。1.5加上0.9到0.4的惯性权重搭配探索和开发平衡得更舒服。如果你发现结果波动大可以尝试把c1降到1.2、c2提到1.8增强向群体最优靠拢的趋势。聚类数K用轮廓系数辅助判断。轮廓系数综合考虑了簇内紧密度和簇间分离度取值范围-1到1越大代表聚类效果越好。我在项目里对K2到K8分别跑PSO-Kmeans计算每个K下的平均轮廓系数选峰值对应的K。实际业务上K取4或5比较常见这样每一类用户都有足够明确的画像不会分得过细而失去解释力。4. 实验效果分析与聚类结果解读4.1 与标准Kmeans的对比实验我拿来验证的数据是某市一个台区3000户居民用户30天的用电记录按前文方法清洗和特征化后得到3000×28的特征矩阵聚类目标K4PSO种群取30迭代50次。为了控制变量标准Kmeans我用Matlab自带的kmeans函数跑100次随机初始化取SSE最小的一次作为参照这种多次随机取最优本身就是实践中应对Kmeans不稳定的常见手段但计算开销远高于PSO辅助。最终实验数据如下表方案平均SSE最优SSE波动范围SSE单次运行耗时标准Kmeans单次1846.71752.3160.40.8秒标准Kmeans100次取最优1635.21635.2024秒PSO-Kmeans单次1658.11641.533.218秒PSO-Kmeans3次取最优1642.01641.53.254秒几个结论很直观。PSO-Kmeans单次结果明显优于Kmeans单次SSE从1846.7降到1658.1下降了大约10.2%即使对比Kmeans跑100次取最优的1635.2PSO-Kmeans的最优SSE 1641.5也非常接近差了不到0.4%。更关键的是稳定性——PSO-Kmeans三次运行的最优与最差只差33.2几乎都在同一水平线上这说明算法已经不太受随机初始化的影响而标准Kmeans单次运行的波动范围高达160以上这在工程上非常致命。当然PSO-Kmeans也不是免费的午餐18秒的处理时间比标准Kmeans单次0.8秒慢得多。但对离线用户画像分析这种场景18秒完全可接受。4.2 聚类结果如何映射到用电行为聚类跑完只是第一步更重要的工作是把每一类用户的行为模式描述出来。我是这样做的拿到聚类标签后把原始负荷数据未标准化按类分组计算每类用户的平均24小时负荷曲线然后结合日用电量、峰谷比等业务指标做解读。在我的实验里K4时的四类用户画像如下第一类工作日早、晚双峰特别突出白天负荷很低午间有小幅回落周末曲线相对平缓。结合日用电量处于中低水平可以判定为典型的上班族家庭工作日只有早晚在家用电。第二类白天负荷较高曲线全天相对平稳夜晚略降但不会降到很低日用电量处于中上水平。这是全天居家型用户可能是老人、家庭主妇或自由职业者。第三类夜间和凌晨负荷异常偏高白天反而较低日用电量也比较大。结合当地电价政策这类用户很可能是有意将洗衣机、热水器等大功率设备挪到夜间使用甚至可能有电动汽车充电行为。第四类整体负荷水平低曲线平缓无峰长时间维持很小的用电功率。这种通常是空心户或者出租率较高的房屋用电行为不活跃。每类用户对应的策略建议也不一样第一类适合宣传分时电价引导削峰填谷第二类可以推荐节能设备第三类可以作为需求响应的重点对象第四类则需要在台区管理上排查是否有空置房或者表计异常。这些业务层面的延伸才是分析工作真正产生价值的地方。4.3 可视化技巧如何把聚类结果画得让业务方看懂聚类结果可视化我踩过不少坑。最开始我直接用plot画所有用户的原始曲线3000条线叠在一起密密麻麻根本看不出差异。后来改成每个类画一条平均曲线标准差带效果立刻不一样。Matlab里用fill可以画带meanCurve mean(clusterData, 1); stdCurve std(clusterData, 1); t 1:24; fill([t fliplr(t)], [meanCurvestdCurve fliplr(meanCurve-stdCurve)], ... [0.9 0.9 0.9], FaceAlpha, 0.4, EdgeColor, none); hold on; plot(t, meanCurve, LineWidth, 2);标准差带能够直观表达这一类用户内部的波动程度。如果某类的带很窄说明这类用户的负荷形态高度一致聚类可信度高带很宽则说明这一类内部还存在细分可以考虑是否增加K值。另外一个可视化技巧是降维散点图。高维特征矩阵不好直接展示可以用t-SNE或者PCA降到2维再按聚类标签着色。不过我要提醒一句降维后再看聚类是否分得开只能作为辅助参考因为降维过程会扭曲真实距离关系。业务汇报时这东西很好看内部验证时别太当真。5. 常见问题与排查技巧实录5.1 粒子维度爆炸和计算速度慢怎么办我在96维特征上尝试过直接跑PSO-Kmeans粒子维度是(K \times 96)K取4就是384维。粒子群优化在这么高的维度上进行搜索效果非常差——适应度收敛慢、粒子群容易散开、结果还不稳定。因为高维空间里距离度量变得稀疏随机初始化的粒子互相之间差异很小PSO很难通过对比分辨哪个方向更好。解决思路有两个。第一是在特征层面降维比如用24小时均值替代96点数据或者先用PCA把特征压到15到20维再做聚类。第二是改变PSO的搜索策略比如将速度初始化设置为0限制粒子的搜索半径但这样又会牺牲全局搜索能力。我的建议是优先做特征降维因为居民负荷数据本身的冗余度很高96个采样点之间存在很强的时序相关性强行保留全部维度得不偿失。计算速度问题还有另一层来源适应度函数里频繁的矩阵运算。如果循环写的效率低几千样本都够让Matlab卡上几分钟。我把计算距离的代码从for循环改成矩阵广播后原来45秒一次迭代缩到3秒左右。Matlab效率的关键就是不要让循环套循环多用维度广播和矩阵运算如果还想更快可以把calcFitness写成mex函数或者用parfor并行计算粒子群中不同粒子的适应度。5.2 陷入局部最优的判断与处理有一种情况PSO迭代结束后gbest对应的质心组其实还不是理想解Kmeans再迭代也跳不出来。怎么判断我会把PSO-Kmeans的SSE和多次随机初始化的Kmeans最优SSE做对比如果前者显著大于后者基本可以断定PSO阶段早收敛了。处理办法有这么几种。一是检查粒子群初始化如果初始粒子全都挤在样本集中的区域多样性不够PSO很容易早熟。初始化时除了随机采样样本点我还会刻意加几个远离中心的点。二是增大惯性权重或者调节学习因子如果w从0.9降到0.4太快个体经验权重过大可以在实验中把wMax提到1.0wMin提到0.5让粒子飞得更激进一点。三是重启策略如果一个粒子连续N代都没有改进自己的pbest给它重新初始化到随机位置这是个简单但很有效的办法。我再强调一次PSO-Kmeans不是银弹它只能显著降低落入局部最优的概率不能完全消除。所以在项目落地时我通常跑3次PSO-Kmeans取SSE最小的那次。由于单次已经很稳定3次取最优带来的额外收益也有限更多是买个心理保险。5.3 K值怎么选最合理选择K值最常见的是肘部法则画SSE随K变化的折线图找那个拐点。但实际数据里肘部往往不明显SSE下降曲线保持平滑你很难说出3和4哪个是肘。我用轮廓系数配合业务可解释性一起判断。轮廓系数对第i个样本的定义是[ s_i \frac{b_i - a_i}{\max(a_i, b_i)} ]其中(a_i)是样本i与同簇其他样本的平均距离(b_i)是样本i与最近其他簇的平均距离。把全部样本的轮廓系数平均就是总体轮廓系数。我一般要求总体轮廓系数大于等于0.5如果某个K下只有0.3说明簇内不够紧凑或者簇间分得不清楚这个K值基本不可用。但我也要说业务可解释性有时候比数值指标更关键。比如K5时轮廓系数最高但其中有一类用户曲线形态和另一类非常接近业务上完全无法区分和应对那K5就没有实际意义。我的习惯是先选2到3个候选K轮廓系数比较高的然后把这几个K下的聚类结果拿给业务同事看问哪一版最容易讲故事通常答案很明确。5.4 版本兼容和Matlab环境的坑我在实验过程中遇到过一次运行环境导致的怪问题在Matlab R2021b上能正常运行的脚本换到老版本后kmeans的Start参数传矩阵就报错。Matlab每个版本对聚类函数输入参数的校验机制不一样如果公司或实验室的Matlab版本不统一建议不要依赖版本较新的参数特性。我的做法是手写一个20轮的Kmeans精修函数替代内置的kmeans代码不超过30行却能在所有版本上稳定运行。核心逻辑就是循环分配样本—更新质心和我们第一部分讲的Kmeans原理完全一致。另外如果你跟我一样被工程化逼得没有正版授权也可以考虑用GNU Octave代替Matlab写这个流程。Octave对大部分数值计算和矩阵运算的支持都很好PSO-Kmeans这种以矩阵运算为主的代码迁移成本很低。不过Octave的kmeans函数不是内置的需要自己手写用来替代内置函数时正好省了上面的兼容性问题。5.5 数据质量细节这些坑会影响聚类结论最后分享几个和算法无关但直接影响结论的数据细节。第一归一化必须在缺失值处理之后做否则Z-score会把缺失值当成0参与均值计算扭曲特征分布。第二聚类的输入应该是行为特征不应该直接放日期、用户编号、台区编号这些标识性变量。第三如果用户数据的天数不一致有的用户只有15天记录有的有30天建议先按用户求平均再做聚类否则天数少的用户会被当成异常样本。第四季节因素要重视——冬季和夏季的负荷曲线形态差异很大如果你直接拿一整年数据混在一起聚类得到的分群往往是季节分群而非行为分群。我的做法是按季节分别建模型然后在业务层面对比同一用户的季节归属变化这样既能识别行为差异又能捕捉季节性规律变化。这套组合方案跑下来我最大的体会是算法层面没有太多高大上的东西PSO-Kmeans本质上是把一个简单而顽固的问题——初始质心敏感——用群智能算法解决掉了。居民用电行为分析的价值也不在于把轮廓系数从0.55提高到0.6而在于每一类用户分出来以后你能针对性地做点什么。最后再分享一个小技巧给准备落地的朋友就算聚类结果已经稳定也别直接信任数据去抽查10个用户的原始负荷曲线和聚类标签是否匹配。光看平均曲线会骗人单条曲线才暴露真相。这个步骤花不了十分钟却能避免向业务方汇报时被一句我看这明显不是一类用户问得哑口无言。