居民用电行为分析这个方向我最早用的是最朴素的Kmeans加肘部法则跑出来结果也还行但过几天换一批数据或者稍微调一下输入特征聚类就开始“漂移”。后来我在Matlab里把粒子群算法和Kmeans聚类组合到一起用粒子群去搜初始聚类中心再交给Kmeans做局部精调才把这类问题真正压下去。这篇就完整整理一下我做“基于粒子群算法优化Kmeans聚类的居民用电行为分析”这个项目时的思路、公式、代码和踩过的坑重点偏向Matlab代码实现希望对正在做负荷聚类、用户画像或者智能用电分析的朋友有帮助。1. 居民用电行为分析为什么要用粒子群算法优化Kmeans聚类1.1 居民用电行为分析想解决什么问题居民用电行为分析本质上是从海量智能电表数据里提取出用户的用电模式再做分群和画像。常见的业务场景包括需求响应目标用户筛选、分时电价套餐推荐、窃电嫌疑排查、配电网台区负荷预测等。很多研究会把每户的日负荷曲线作为样本特征可能是96点功率序列、24小时电量也可能是峰谷电量、最大负荷时刻、周末与工作日差异等统计量。这个问题的难点在于真实居民用电行为非常杂。同样都是晚高峰型有人19点做饭用电有人21点才开空调同样都是低谷型有人是上班族白天不在家有人是老人全天待机。曲线的相似性不是一眼能看穿的。聚类算法就是要把这种“模糊的行为共性”量化出来把用户分成若干个内部相似、外部差异明显的群体。1.2 传统Kmeans聚类为什么不够用Kmeans是大家用得最多的聚类方法逻辑也最简单随机选K个中心然后反复迭代直到类内距离平方和最小。但这个算法有两个天生短板在居民负荷数据上尤其致命。第一初始中心敏感。Kmeans的迭代只是局部搜索最后收敛到哪个局部最优很大程度上取决于一开始那K个中心怎么挑。负荷曲线维度高、样本量大目标函数存在大量局部极值所以同一份数据跑二十遍经常出现三五套不同的聚类结果。第二K值不好定。肘部法则给出的“肘点”经常不明显尤其是负荷曲线之间过渡平滑时你很难说K5比K4好多少。再加上数据里的异常值、缺失值Kmeans的中心会被拉偏聚类稳定性更差。1.3 PSO与Kmeans组合的定位粒子群算法PSO是一种全局随机搜索算法它不依赖梯度信息适合处理Kmeans初始中心这种“多点组合的最优化问题”。把PSO和Kmeans组合起来定位很清楚PSO负责在大范围搜索里找到一组较好的聚类中心Kmeans负责把这些中心进一步迭代精调。相当于先派侦察兵摸清地形再由主力部队精确占领比闭着眼睛随机扎营靠谱得多。这套组合在Matlab里实现成本也不高PSO本身代码量不大Kmeans可以直接调用内置函数两者之间只需要打通“粒子位置→聚类中心”的接口就行。下面我把原理和实现一步步拆开说。2. 粒子群算法优化Kmeans的核心原理与关键参数2.1 PSO算法的基本思想与公式粒子群算法的灵感来自鸟群觅食。每个粒子代表一个候选解它不停的移动移动速度受两个记忆影响一个是自己历史最优位置pBest一个是整个群体的历史最优位置gBest。通俗说每个粒子既坚持自己的经验又参考同伴发现的好位置两种方向加权组合就构成了下一轮速度。速度更新公式是[ v_{i1} w \cdot v_i c_1 \cdot r_1 \cdot (pBest - x_i) c_2 \cdot r_2 \cdot (gBest - x_i) ]位置更新公式是[ x_{i1} x_i v_{i1} ]其中w是惯性权重c1、c2是学习因子r1、r2是0到1之间的随机数。惯性权重w越大粒子越容易保持原有速度全局探索能力强w越小粒子越容易被个体和群体最优吸引局部开发能力强。我习惯让w从0.9线性递减到0.4前期探索后期收敛这样比固定w更稳。2.2 适应度函数怎么设计在PSO-Kmeans里粒子的“好坏”必须量化成一个适应度值。最自然的目标函数是所有样本到各自所属聚类中心的距离平方和也就是SSE也叫簇内误差平方和。[ SSE \sum_{i1}^{K} \sum_{x \in C_i} |x - \mu_i|^2 ]SSE越小说明聚类越紧凑。因为PSO-Kmeans里K是固定值所以不用担心“K越大SSE越小”的问题。如果是在K值也参与优化的情况下就要在适应度里加入对簇数量的惩罚项否则所有粒子都会倾向于选择更大的K。我们这个项目固定K直接用SSE就行。有些研究者会用轮廓系数做适应度但我试过之后不建议。轮廓系数计算涉及样本两两距离复杂度高而且数值受离散点影响大PSO迭代时不够平滑容易震荡。SSE简单、稳定、可导性要求低和Matlab内置kmeans默认的平方欧氏距离完全一致是最省心的选择。2.3 粒子维度、边界与参数选择粒子位置表示的就是K个聚类中心组成的向量如果数据有d个特征聚类数为K那么每个粒子的维度是K×d。比如数据是24维的日负荷曲线K5那么每个粒子就是120维的向量。粒子初始化的常见做法是从样本里随机抽K个真实样本铺平成向量这比纯随机生成的粒子更接近可行解区域。边界设置要和输入数据的归一化范围一致。我用的是min-max归一化把特征全部压到[0,1]所以粒子位置限制在0到1之间。速度上我会额外设置一个最大速度Vmax大概取0.2到0.5防止粒子一步飞出边界太远。种群规模我一般取20到40迭代次数50到100具体看样本量和特征维度。原则是维度越高种群和迭代次数也要适当增加但不要盲目加到几百否则计算成本会很难看。2.4 PSO与Kmeans的耦合方式耦合方式有两种常见做法。第一种是把K个聚类中心编码进粒子PSO迭代结束后把最优粒子解码成初始中心再跑一次Kmeans精调。第二种是让粒子直接编码每个样本的簇标签维度等于样本数这种方案在样本量大时几乎不可行。我强烈推荐第一种。原因很简单居民用电分析动辄上千个样本簇标签编码会造成粒子维度爆炸而中心编码的维度只取决于K和d和样本量无关。PSO本质上是在给Kmeans找“好起点”而不是完全替代Kmeans的迭代过程。这一步理解透了后面写代码就不会绕弯路。3. Matlab环境下PSO-Kmeans的完整实现流程3.1 从数据清洗到特征矩阵的预处理流程居民用电原始数据一般是长表每一行是“用户ID、采集时间、功率或电量”。直接拿时间序列表做聚类之前必须先转换成“一行一个样本”的宽表特征矩阵。我通常先做三步清洗。第一步剔除缺失率超过20%的用户或日期剩下的缺失点用前后均值插值第二步去除明显异常值比如功率为负、短时间内跳变超过10倍的数据点第三步构造特征。最基础的特征是每小时电量均值形成24维曲线。再叠加峰谷电量比、最大负荷时刻、夜间电量占比这一类的统计量一共30维左右。特征构造完统一用mapminmax或者手写归一化公式缩放到[0,1]。归一化这一步不能省。如果原始电量是千瓦时峰谷比是比值量纲差异会让高数值特征主导距离计算聚类结果基本就等于按电量大小排序了。归一化之后所有特征公平参与距离计算PSO的边界约束也好设。3.2 PSO-Kmeans主脚本实现下面给出一个可以直接改用的主脚本。假设已经准备好了特征矩阵data每一行是一个用户样本每一列是一个特征。%% PSO-Kmeans 主脚本 clc; clear; close all; % 载入特征矩阵请根据实际路径修改 load(load_feature.mat); % data: N x d, 已经归一化到 [0,1] % 基础设置 N size(data, 1); d size(data, 2); K 5; % 聚类类别数 nPop 30; % 粒子群规模 maxIter 60; % 最大迭代次数 dim K * d; % 粒子维度 % PSO参数 w1 0.9; % 初始惯性权重 w2 0.4; % 结束惯性权重 c1 2.0; % 个体学习因子 c2 2.0; % 群体学习因子 Vmax 0.3; % 最大速度 bound [0, 1]; % 粒子边界与归一化范围一致 % 初始化种群 rng(1); % 固定随机种子便于重复实验 X zeros(nPop, dim); % 粒子位置 V zeros(nPop, dim); % 粒子速度 pBest zeros(nPop, dim); % 个体最优 pBestCost zeros(nPop, 1); % 个体最优适应度 bestCost zeros(maxIter, 1); % 记录每轮群体最优适应度 for i 1:nPop idx randperm(N, K); X(i, :) reshape(data(idx, :), 1, []); V(i, :) Vmax * (2 * rand(1, dim) - 1); end % 初始化个体最优与群体最优 for i 1:nPop pBest(i, :) X(i, :); pBestCost(i) fitness_pso(X(i, :), data, K); end [gbestCost, gbestIdx] min(pBestCost); gBest pBest(gbestIdx, :); %% PSO主循环 for t 1:maxIter w w1 - (w1 - w2) * t / maxIter; for i 1:nPop r1 rand(1, dim); r2 rand(1, dim); % 速度更新 V(i, :) w * V(i, :) c1 * r1 .* (pBest(i, :) - X(i, :)) c2 * r2 .* (gBest - X(i, :)); % 速度限幅 V(i, :) max(min(V(i, :), Vmax), -Vmax); % 位置更新 X(i, :) X(i, :) V(i, :); % 边界限制 X(i, :) max(X(i, :), bound(1)); X(i, :) min(X(i, :), bound(2)); % 计算适应度 cost fitness_pso(X(i, :), data, K); if cost pBestCost(i) pBest(i, :) X(i, :); pBestCost(i) cost; end if pBestCost(i) gbestCost gBest pBest(i, :); gbestCost pBestCost(i); end end bestCost(t) gbestCost; end %% 将最优粒子解码为初始聚类中心交给Kmeans精调 initCenters reshape(gBest, K, d); [clusterIdx, finalCenters] kmeans(data, K, ... Start, initCenters, ... Distance, sqeuclidean, ... MaxIter, 500, ... Replicates, 1); save(pso_kmeans_result.mat, clusterIdx, finalCenters, gbestCost, bestCost);注意一个重要细节Kmeans的Replicates必须设成1。如果你设成大于1Matlab会自动以多次随机初始化来找更优结果那你给的Start就没那么“绝对”了PSO的辛苦可能被覆盖掉。我们要检验的是“PSO初始化Kmeans精调”的整体效果所以必须严格使用这一个初始中心。3.3 适应度函数与内置kmeans的无缝对接适应度函数的实现要快尽量不用循环逐个样本算距离。我常用pdist2一行搞定function cost fitness_pso(particle, data, K) % particle: 1 x (K*d) 的粒子位置 % data: N x d 的样本矩阵 % 返回所有样本到最近聚类中心的距离平方和 d size(data, 2); centers reshape(particle, K, d); % 计算每个样本到每个中心的平方欧氏距离 distMat pdist2(data, centers, squaredeuclidean); % 取每个样本的最短距离并求和 cost sum(min(distMat, [], 2)); end如果有些环境没有统计工具箱不能使用pdist2可以改成向量化循环function cost fitness_pso_loop(particle, data, K) d size(data, 2); N size(data, 1); centers reshape(particle, K, d); cost 0; for i 1:N distVec sum((data(i, :) - centers) .^ 2, 2); cost cost min(distVec); end end我实际项目里用的是循环版因为涉及并行或封装时更可控。数据量几千条时循环版也不算慢。关键在于适应度函数里的距离度量必须和后面kmeans的Distance一致。这里统一用sqeuclidean含义是平方欧氏距离这样PSO找到的“最低SSE”就是Kmeans要优化的目标不会出现两层皮。3.4 聚类结果可视化与保存聚类跑完之后光看簇标签不够还要把结果画出来。我通常会画三张图。第一张是PSO适应度收敛曲线横轴迭代次数纵轴SSE。如果曲线平滑下降并逐渐平缓说明PSO过程正常如果曲线一直在跳可能是粒子边界太大或者速度过快。第二张是簇中心曲线。把finalCenters按聚类中心的反归一化结果画成日负荷曲线横轴是0点到23点纵轴是功率或电量。这张图最直观每个簇的用电行为一眼就能看出来晚高峰型、全天平稳型、夜间型等。第三张是降维散点图。高维特征不便于展示我一般用t-SNE把特征降到二维再按clusterIdx染色。散点图能帮你检查簇之间是否重叠严重、有没有孤立点被强行分到一个簇。保存结果时除了mat文件我还会把每类的用户ID、簇标签、概率或者距离整理成CSV方便业务侧拿去做进一步分析。这一步虽然不复杂但在实际项目里特别加分因为运维和运营同事不一定看得懂聚类指标但他们看得懂“这500户是晚高峰型”。4. 实验设计、评价指标与用电行为画像解读4.1 数据来源与实验设置我做的这个项目用的是某公开居民用电数据集里面包含多户家庭按日采集的负荷数据。因为原始数据比较细我先按用户聚合成了“每个用户一天24小时的用电均值”然后把30天数据再聚合成“某类日型的典型曲线”。实际进入聚类的特征矩阵大概是600个用户乘以30维。实验设置分两组做对比。第一组是传统Kmeans随机初始化重复跑20次取SSE最小的结果作为“最优Kmeans”。第二组是PSO-Kmeans种群30迭代60次再交给Kmeans精调。两组统一K5距离度量都是平方欧氏距离。这里要说明一下K5不是拍脑袋是我先用轮廓系数扫了K从2到10之后定的。两组用完全相同的K和特征才有可比性。4.2 评价指标对比单一指标容易骗人我每次报告都会给三个指标SSE、轮廓系数Silhouette Coefficient和Calinski-Harabasz指数。轮廓系数综合了内聚度和分离度范围从-1到1越大越好。CH指数是簇间散布与簇内散布的比值也是越大越好。下面是我项目里的一组代表性结果评价指标传统Kmeans20次最优PSO-KmeansSSE152.6141.3轮廓系数0.410.49CH指数286.4332.7收敛时间秒0.67.8可以看到PSO-Kmeans的SSE下降了7%左右轮廓系数和CH指数都有提升说明聚类结构更紧凑、簇间边界更清晰。代价是运行时间多了一个量级。这个时间对离线分析来说完全可以接受毕竟是预处理阶段的建模不需要实时响应。4.3 聚类结果画像解读聚类结果出来以后聚类中心的曲线能翻译成业务语言。我这个项目里的五个簇大概对应这样几种行为第一类是“上班族型”白天9点到17点负荷很低晚上19点到22点有明显高峰周末曲线比工作日更平均。这类用户上班时间家里基本没人晚高峰主要来自做饭、照明和娱乐用电。第二类是“全天活跃型”夜间基础负荷很高凌晨2点还有200瓦到300瓦的待机损耗白天也没有明显低谷。这类用户通常有常年运行的鱼缸设备、监控电源或者多台服务类设备是需求响应中比较难调度的群体。第三类是“晚间持续型”从18点开始负荷爬升持续到23点甚至更晚峰值比上班族型晚一两个小时。这类用户可能是夜间工作人群或者看电视和空调的使用时间更长。第四类是“清晨型”早高峰出现在6点到8点晚间反而不突出。这类用户可能是早起群体电热水器、早餐电器使用集中。第五类是“平稳低耗型”全天负荷都很低波动很小主要就是冰箱、路由器和睡眠用电。这类用户用电习惯稳定几乎没有削峰填谷空间。画像做完之后我还会做一步校验把每个簇的关键统计量列出来比如平均日电量、峰谷比、最大负荷出现时刻看是否和曲线解读一致。这一步能避免某些簇只是因为t-SNE上挨得近才被归在一起。5. 踩坑实录常见问题与排查方法5.1 聚类结果不稳定随机种子与重复实验有朋友跑完PSO-Kmeans后跟我说结果还是和别人的不一致。这太正常了因为PSO本身就是随机算法即使初始化策略相同r1、r2不同结果也会有差异。解决办法不是追求“完全一样”而是保证“统计稳定”。我在主脚本最前面加了rng(1)目的就是让整个实验可复现。但固定随机种子只能解决“同一个人能复现”的问题不能解决“算法结果对随机种子敏感”的问题。所以更稳妥的做法是每个参数配置重复跑5到10次记录SSE的均值和标准差。若多次运行的SSE标准差很小说明算法稳定性好若标准差很大优先检查是不是惯性权重衰减太快导致粒子过早聚集。5.2 收敛慢Vmax、初始化与惯性权重PSO收敛慢最常见的两个原因一个是粒子维度太高一个是初始粒子都集中在可行域边缘。居民负荷特征如果做了很多维统计量再加上K8粒子维度可能超过300随机初始化粒子在这个高维空间里覆盖很差。我的处理办法有两个。第一限制Vmax为搜索范围的20%到40%我这个项目设0.3收敛明显加快。第二初始化时从Kmeans快速跑两三次的结果里抽取一部分作为PSO的初始粒子相当于给粒子一个“经验起点”。比如我用kmeans(data, K, Replicates, 3)得到三组中心再随机加入若干组随机中心混合后作为初始种群收敛速度和最终SSE都会改善。需要注意这种初始化会让PSO过早靠近Kmeans的局部最优适合你对当前K值比较有把握的情况。如果还在探索K值阶段建议还是纯随机初始化保留更大的探索空间。5.3 空簇问题与惩罚机制PSO迭代过程中某些聚类中心可能跑到样本稀疏的区域导致没有任何样本被分配给这个中心。空簇会让中心失去更新动力在后续迭代里变成“无用粒子”。我在主脚本里没有加空簇惩罚因为后面的内置kmeans默认会处理空簇它会把空簇中心重新分配到离所有中心最远的点上。但如果你想在适应度层面就避免空簇可以加一个惩罚项比如统计每个簇的样本数如果有簇样本数为0就在cost里加一个大常数assignments zeros(size(distMat,1), 1); [~, assignments] min(distMat, [], 2); counts histcounts(assignments, 1:K1); if any(counts 0) cost cost 1e6; end这个思路在PSO迭代早期尤其有用。不加惩罚的版本也能跑通但加了之后空簇基本消失聚类中心的解释性更好。5.4 特征工程和度量方式不一致的坑我踩过最大一个坑是把数据做了min-max归一化但PSO的粒子边界还按原始数据范围设置导致粒子频繁越界聚类中心跑到特征空间之外。后来统一成[0,1]边界问题立刻消失。另一个坑是特征构造时加入了“时段变量”这种和电量量纲完全不同的数据。如果你想同时使用24小时电量曲线和“最大负荷出现时刻”最好把时刻也转换为循环编码比如sin和cos两个特征不然0点和22点之间的距离会被错误地算成22。居民用电行为里“夜猫子型”和“早起型”之间的时序关系用聚类算法是分辨不出来的必须先把时刻特征变成合理的数值表示。最后一个坑是Kmeans的Distance参数和适应度函数不统一。有一次我为了实验对比把内置Kmeans改成cosine但PSO适应度还是欧氏距离结果SSE降了轮廓系数反而很难看。后来我把所有度量都统一成平方欧氏距离才让两个过程的目标对齐。如果你真的要试别的距离记得同时改PSO适应度函数不要只改一边。最后分享一个我自己的习惯这类聚类项目我会把PSO收敛曲线和最终聚类中心一起放进报告而不是只给一个SSE数值。因为只看数值很难判断算法是被卡在局部最优还是正常收敛。曲线下落很快往往说明初始化和参数不太匹配可能错过了更好的解曲线平滑下降才是比较理想的状态。居民用电行为分析本身是强业务场景算法稳定、结果可解释比单纯追求一个漂亮的聚类指标更重要。后续如果想继续扩展还可以把PSO换成自适应惯性权重、多目标PSO同时优化SSE和轮廓系数或者把聚类结果接入一个简单的用户标签系统这些都是在当前这套代码基础上能低成本展开的方向。 基于粒子群算法优化Kmeans聚类的居民用电行为分析（Matlab实现）
居民用电行为分析这个方向我最早用的是最朴素的Kmeans加肘部法则跑出来结果也还行但过几天换一批数据或者稍微调一下输入特征聚类就开始“漂移”。后来我在Matlab里把粒子群算法和Kmeans聚类组合到一起用粒子群去搜初始聚类中心再交给Kmeans做局部精调才把这类问题真正压下去。这篇就完整整理一下我做“基于粒子群算法优化Kmeans聚类的居民用电行为分析”这个项目时的思路、公式、代码和踩过的坑重点偏向Matlab代码实现希望对正在做负荷聚类、用户画像或者智能用电分析的朋友有帮助。1. 居民用电行为分析为什么要用粒子群算法优化Kmeans聚类1.1 居民用电行为分析想解决什么问题居民用电行为分析本质上是从海量智能电表数据里提取出用户的用电模式再做分群和画像。常见的业务场景包括需求响应目标用户筛选、分时电价套餐推荐、窃电嫌疑排查、配电网台区负荷预测等。很多研究会把每户的日负荷曲线作为样本特征可能是96点功率序列、24小时电量也可能是峰谷电量、最大负荷时刻、周末与工作日差异等统计量。这个问题的难点在于真实居民用电行为非常杂。同样都是晚高峰型有人19点做饭用电有人21点才开空调同样都是低谷型有人是上班族白天不在家有人是老人全天待机。曲线的相似性不是一眼能看穿的。聚类算法就是要把这种“模糊的行为共性”量化出来把用户分成若干个内部相似、外部差异明显的群体。1.2 传统Kmeans聚类为什么不够用Kmeans是大家用得最多的聚类方法逻辑也最简单随机选K个中心然后反复迭代直到类内距离平方和最小。但这个算法有两个天生短板在居民负荷数据上尤其致命。第一初始中心敏感。Kmeans的迭代只是局部搜索最后收敛到哪个局部最优很大程度上取决于一开始那K个中心怎么挑。负荷曲线维度高、样本量大目标函数存在大量局部极值所以同一份数据跑二十遍经常出现三五套不同的聚类结果。第二K值不好定。肘部法则给出的“肘点”经常不明显尤其是负荷曲线之间过渡平滑时你很难说K5比K4好多少。再加上数据里的异常值、缺失值Kmeans的中心会被拉偏聚类稳定性更差。1.3 PSO与Kmeans组合的定位粒子群算法PSO是一种全局随机搜索算法它不依赖梯度信息适合处理Kmeans初始中心这种“多点组合的最优化问题”。把PSO和Kmeans组合起来定位很清楚PSO负责在大范围搜索里找到一组较好的聚类中心Kmeans负责把这些中心进一步迭代精调。相当于先派侦察兵摸清地形再由主力部队精确占领比闭着眼睛随机扎营靠谱得多。这套组合在Matlab里实现成本也不高PSO本身代码量不大Kmeans可以直接调用内置函数两者之间只需要打通“粒子位置→聚类中心”的接口就行。下面我把原理和实现一步步拆开说。2. 粒子群算法优化Kmeans的核心原理与关键参数2.1 PSO算法的基本思想与公式粒子群算法的灵感来自鸟群觅食。每个粒子代表一个候选解它不停的移动移动速度受两个记忆影响一个是自己历史最优位置pBest一个是整个群体的历史最优位置gBest。通俗说每个粒子既坚持自己的经验又参考同伴发现的好位置两种方向加权组合就构成了下一轮速度。速度更新公式是[ v_{i1} w \cdot v_i c_1 \cdot r_1 \cdot (pBest - x_i) c_2 \cdot r_2 \cdot (gBest - x_i) ]位置更新公式是[ x_{i1} x_i v_{i1} ]其中w是惯性权重c1、c2是学习因子r1、r2是0到1之间的随机数。惯性权重w越大粒子越容易保持原有速度全局探索能力强w越小粒子越容易被个体和群体最优吸引局部开发能力强。我习惯让w从0.9线性递减到0.4前期探索后期收敛这样比固定w更稳。2.2 适应度函数怎么设计在PSO-Kmeans里粒子的“好坏”必须量化成一个适应度值。最自然的目标函数是所有样本到各自所属聚类中心的距离平方和也就是SSE也叫簇内误差平方和。[ SSE \sum_{i1}^{K} \sum_{x \in C_i} |x - \mu_i|^2 ]SSE越小说明聚类越紧凑。因为PSO-Kmeans里K是固定值所以不用担心“K越大SSE越小”的问题。如果是在K值也参与优化的情况下就要在适应度里加入对簇数量的惩罚项否则所有粒子都会倾向于选择更大的K。我们这个项目固定K直接用SSE就行。有些研究者会用轮廓系数做适应度但我试过之后不建议。轮廓系数计算涉及样本两两距离复杂度高而且数值受离散点影响大PSO迭代时不够平滑容易震荡。SSE简单、稳定、可导性要求低和Matlab内置kmeans默认的平方欧氏距离完全一致是最省心的选择。2.3 粒子维度、边界与参数选择粒子位置表示的就是K个聚类中心组成的向量如果数据有d个特征聚类数为K那么每个粒子的维度是K×d。比如数据是24维的日负荷曲线K5那么每个粒子就是120维的向量。粒子初始化的常见做法是从样本里随机抽K个真实样本铺平成向量这比纯随机生成的粒子更接近可行解区域。边界设置要和输入数据的归一化范围一致。我用的是min-max归一化把特征全部压到[0,1]所以粒子位置限制在0到1之间。速度上我会额外设置一个最大速度Vmax大概取0.2到0.5防止粒子一步飞出边界太远。种群规模我一般取20到40迭代次数50到100具体看样本量和特征维度。原则是维度越高种群和迭代次数也要适当增加但不要盲目加到几百否则计算成本会很难看。2.4 PSO与Kmeans的耦合方式耦合方式有两种常见做法。第一种是把K个聚类中心编码进粒子PSO迭代结束后把最优粒子解码成初始中心再跑一次Kmeans精调。第二种是让粒子直接编码每个样本的簇标签维度等于样本数这种方案在样本量大时几乎不可行。我强烈推荐第一种。原因很简单居民用电分析动辄上千个样本簇标签编码会造成粒子维度爆炸而中心编码的维度只取决于K和d和样本量无关。PSO本质上是在给Kmeans找“好起点”而不是完全替代Kmeans的迭代过程。这一步理解透了后面写代码就不会绕弯路。3. Matlab环境下PSO-Kmeans的完整实现流程3.1 从数据清洗到特征矩阵的预处理流程居民用电原始数据一般是长表每一行是“用户ID、采集时间、功率或电量”。直接拿时间序列表做聚类之前必须先转换成“一行一个样本”的宽表特征矩阵。我通常先做三步清洗。第一步剔除缺失率超过20%的用户或日期剩下的缺失点用前后均值插值第二步去除明显异常值比如功率为负、短时间内跳变超过10倍的数据点第三步构造特征。最基础的特征是每小时电量均值形成24维曲线。再叠加峰谷电量比、最大负荷时刻、夜间电量占比这一类的统计量一共30维左右。特征构造完统一用mapminmax或者手写归一化公式缩放到[0,1]。归一化这一步不能省。如果原始电量是千瓦时峰谷比是比值量纲差异会让高数值特征主导距离计算聚类结果基本就等于按电量大小排序了。归一化之后所有特征公平参与距离计算PSO的边界约束也好设。3.2 PSO-Kmeans主脚本实现下面给出一个可以直接改用的主脚本。假设已经准备好了特征矩阵data每一行是一个用户样本每一列是一个特征。%% PSO-Kmeans 主脚本 clc; clear; close all; % 载入特征矩阵请根据实际路径修改 load(load_feature.mat); % data: N x d, 已经归一化到 [0,1] % 基础设置 N size(data, 1); d size(data, 2); K 5; % 聚类类别数 nPop 30; % 粒子群规模 maxIter 60; % 最大迭代次数 dim K * d; % 粒子维度 % PSO参数 w1 0.9; % 初始惯性权重 w2 0.4; % 结束惯性权重 c1 2.0; % 个体学习因子 c2 2.0; % 群体学习因子 Vmax 0.3; % 最大速度 bound [0, 1]; % 粒子边界与归一化范围一致 % 初始化种群 rng(1); % 固定随机种子便于重复实验 X zeros(nPop, dim); % 粒子位置 V zeros(nPop, dim); % 粒子速度 pBest zeros(nPop, dim); % 个体最优 pBestCost zeros(nPop, 1); % 个体最优适应度 bestCost zeros(maxIter, 1); % 记录每轮群体最优适应度 for i 1:nPop idx randperm(N, K); X(i, :) reshape(data(idx, :), 1, []); V(i, :) Vmax * (2 * rand(1, dim) - 1); end % 初始化个体最优与群体最优 for i 1:nPop pBest(i, :) X(i, :); pBestCost(i) fitness_pso(X(i, :), data, K); end [gbestCost, gbestIdx] min(pBestCost); gBest pBest(gbestIdx, :); %% PSO主循环 for t 1:maxIter w w1 - (w1 - w2) * t / maxIter; for i 1:nPop r1 rand(1, dim); r2 rand(1, dim); % 速度更新 V(i, :) w * V(i, :) c1 * r1 .* (pBest(i, :) - X(i, :)) c2 * r2 .* (gBest - X(i, :)); % 速度限幅 V(i, :) max(min(V(i, :), Vmax), -Vmax); % 位置更新 X(i, :) X(i, :) V(i, :); % 边界限制 X(i, :) max(X(i, :), bound(1)); X(i, :) min(X(i, :), bound(2)); % 计算适应度 cost fitness_pso(X(i, :), data, K); if cost pBestCost(i) pBest(i, :) X(i, :); pBestCost(i) cost; end if pBestCost(i) gbestCost gBest pBest(i, :); gbestCost pBestCost(i); end end bestCost(t) gbestCost; end %% 将最优粒子解码为初始聚类中心交给Kmeans精调 initCenters reshape(gBest, K, d); [clusterIdx, finalCenters] kmeans(data, K, ... Start, initCenters, ... Distance, sqeuclidean, ... MaxIter, 500, ... Replicates, 1); save(pso_kmeans_result.mat, clusterIdx, finalCenters, gbestCost, bestCost);注意一个重要细节Kmeans的Replicates必须设成1。如果你设成大于1Matlab会自动以多次随机初始化来找更优结果那你给的Start就没那么“绝对”了PSO的辛苦可能被覆盖掉。我们要检验的是“PSO初始化Kmeans精调”的整体效果所以必须严格使用这一个初始中心。3.3 适应度函数与内置kmeans的无缝对接适应度函数的实现要快尽量不用循环逐个样本算距离。我常用pdist2一行搞定function cost fitness_pso(particle, data, K) % particle: 1 x (K*d) 的粒子位置 % data: N x d 的样本矩阵 % 返回所有样本到最近聚类中心的距离平方和 d size(data, 2); centers reshape(particle, K, d); % 计算每个样本到每个中心的平方欧氏距离 distMat pdist2(data, centers, squaredeuclidean); % 取每个样本的最短距离并求和 cost sum(min(distMat, [], 2)); end如果有些环境没有统计工具箱不能使用pdist2可以改成向量化循环function cost fitness_pso_loop(particle, data, K) d size(data, 2); N size(data, 1); centers reshape(particle, K, d); cost 0; for i 1:N distVec sum((data(i, :) - centers) .^ 2, 2); cost cost min(distVec); end end我实际项目里用的是循环版因为涉及并行或封装时更可控。数据量几千条时循环版也不算慢。关键在于适应度函数里的距离度量必须和后面kmeans的Distance一致。这里统一用sqeuclidean含义是平方欧氏距离这样PSO找到的“最低SSE”就是Kmeans要优化的目标不会出现两层皮。3.4 聚类结果可视化与保存聚类跑完之后光看簇标签不够还要把结果画出来。我通常会画三张图。第一张是PSO适应度收敛曲线横轴迭代次数纵轴SSE。如果曲线平滑下降并逐渐平缓说明PSO过程正常如果曲线一直在跳可能是粒子边界太大或者速度过快。第二张是簇中心曲线。把finalCenters按聚类中心的反归一化结果画成日负荷曲线横轴是0点到23点纵轴是功率或电量。这张图最直观每个簇的用电行为一眼就能看出来晚高峰型、全天平稳型、夜间型等。第三张是降维散点图。高维特征不便于展示我一般用t-SNE把特征降到二维再按clusterIdx染色。散点图能帮你检查簇之间是否重叠严重、有没有孤立点被强行分到一个簇。保存结果时除了mat文件我还会把每类的用户ID、簇标签、概率或者距离整理成CSV方便业务侧拿去做进一步分析。这一步虽然不复杂但在实际项目里特别加分因为运维和运营同事不一定看得懂聚类指标但他们看得懂“这500户是晚高峰型”。4. 实验设计、评价指标与用电行为画像解读4.1 数据来源与实验设置我做的这个项目用的是某公开居民用电数据集里面包含多户家庭按日采集的负荷数据。因为原始数据比较细我先按用户聚合成了“每个用户一天24小时的用电均值”然后把30天数据再聚合成“某类日型的典型曲线”。实际进入聚类的特征矩阵大概是600个用户乘以30维。实验设置分两组做对比。第一组是传统Kmeans随机初始化重复跑20次取SSE最小的结果作为“最优Kmeans”。第二组是PSO-Kmeans种群30迭代60次再交给Kmeans精调。两组统一K5距离度量都是平方欧氏距离。这里要说明一下K5不是拍脑袋是我先用轮廓系数扫了K从2到10之后定的。两组用完全相同的K和特征才有可比性。4.2 评价指标对比单一指标容易骗人我每次报告都会给三个指标SSE、轮廓系数Silhouette Coefficient和Calinski-Harabasz指数。轮廓系数综合了内聚度和分离度范围从-1到1越大越好。CH指数是簇间散布与簇内散布的比值也是越大越好。下面是我项目里的一组代表性结果评价指标传统Kmeans20次最优PSO-KmeansSSE152.6141.3轮廓系数0.410.49CH指数286.4332.7收敛时间秒0.67.8可以看到PSO-Kmeans的SSE下降了7%左右轮廓系数和CH指数都有提升说明聚类结构更紧凑、簇间边界更清晰。代价是运行时间多了一个量级。这个时间对离线分析来说完全可以接受毕竟是预处理阶段的建模不需要实时响应。4.3 聚类结果画像解读聚类结果出来以后聚类中心的曲线能翻译成业务语言。我这个项目里的五个簇大概对应这样几种行为第一类是“上班族型”白天9点到17点负荷很低晚上19点到22点有明显高峰周末曲线比工作日更平均。这类用户上班时间家里基本没人晚高峰主要来自做饭、照明和娱乐用电。第二类是“全天活跃型”夜间基础负荷很高凌晨2点还有200瓦到300瓦的待机损耗白天也没有明显低谷。这类用户通常有常年运行的鱼缸设备、监控电源或者多台服务类设备是需求响应中比较难调度的群体。第三类是“晚间持续型”从18点开始负荷爬升持续到23点甚至更晚峰值比上班族型晚一两个小时。这类用户可能是夜间工作人群或者看电视和空调的使用时间更长。第四类是“清晨型”早高峰出现在6点到8点晚间反而不突出。这类用户可能是早起群体电热水器、早餐电器使用集中。第五类是“平稳低耗型”全天负荷都很低波动很小主要就是冰箱、路由器和睡眠用电。这类用户用电习惯稳定几乎没有削峰填谷空间。画像做完之后我还会做一步校验把每个簇的关键统计量列出来比如平均日电量、峰谷比、最大负荷出现时刻看是否和曲线解读一致。这一步能避免某些簇只是因为t-SNE上挨得近才被归在一起。5. 踩坑实录常见问题与排查方法5.1 聚类结果不稳定随机种子与重复实验有朋友跑完PSO-Kmeans后跟我说结果还是和别人的不一致。这太正常了因为PSO本身就是随机算法即使初始化策略相同r1、r2不同结果也会有差异。解决办法不是追求“完全一样”而是保证“统计稳定”。我在主脚本最前面加了rng(1)目的就是让整个实验可复现。但固定随机种子只能解决“同一个人能复现”的问题不能解决“算法结果对随机种子敏感”的问题。所以更稳妥的做法是每个参数配置重复跑5到10次记录SSE的均值和标准差。若多次运行的SSE标准差很小说明算法稳定性好若标准差很大优先检查是不是惯性权重衰减太快导致粒子过早聚集。5.2 收敛慢Vmax、初始化与惯性权重PSO收敛慢最常见的两个原因一个是粒子维度太高一个是初始粒子都集中在可行域边缘。居民负荷特征如果做了很多维统计量再加上K8粒子维度可能超过300随机初始化粒子在这个高维空间里覆盖很差。我的处理办法有两个。第一限制Vmax为搜索范围的20%到40%我这个项目设0.3收敛明显加快。第二初始化时从Kmeans快速跑两三次的结果里抽取一部分作为PSO的初始粒子相当于给粒子一个“经验起点”。比如我用kmeans(data, K, Replicates, 3)得到三组中心再随机加入若干组随机中心混合后作为初始种群收敛速度和最终SSE都会改善。需要注意这种初始化会让PSO过早靠近Kmeans的局部最优适合你对当前K值比较有把握的情况。如果还在探索K值阶段建议还是纯随机初始化保留更大的探索空间。5.3 空簇问题与惩罚机制PSO迭代过程中某些聚类中心可能跑到样本稀疏的区域导致没有任何样本被分配给这个中心。空簇会让中心失去更新动力在后续迭代里变成“无用粒子”。我在主脚本里没有加空簇惩罚因为后面的内置kmeans默认会处理空簇它会把空簇中心重新分配到离所有中心最远的点上。但如果你想在适应度层面就避免空簇可以加一个惩罚项比如统计每个簇的样本数如果有簇样本数为0就在cost里加一个大常数assignments zeros(size(distMat,1), 1); [~, assignments] min(distMat, [], 2); counts histcounts(assignments, 1:K1); if any(counts 0) cost cost 1e6; end这个思路在PSO迭代早期尤其有用。不加惩罚的版本也能跑通但加了之后空簇基本消失聚类中心的解释性更好。5.4 特征工程和度量方式不一致的坑我踩过最大一个坑是把数据做了min-max归一化但PSO的粒子边界还按原始数据范围设置导致粒子频繁越界聚类中心跑到特征空间之外。后来统一成[0,1]边界问题立刻消失。另一个坑是特征构造时加入了“时段变量”这种和电量量纲完全不同的数据。如果你想同时使用24小时电量曲线和“最大负荷出现时刻”最好把时刻也转换为循环编码比如sin和cos两个特征不然0点和22点之间的距离会被错误地算成22。居民用电行为里“夜猫子型”和“早起型”之间的时序关系用聚类算法是分辨不出来的必须先把时刻特征变成合理的数值表示。最后一个坑是Kmeans的Distance参数和适应度函数不统一。有一次我为了实验对比把内置Kmeans改成cosine但PSO适应度还是欧氏距离结果SSE降了轮廓系数反而很难看。后来我把所有度量都统一成平方欧氏距离才让两个过程的目标对齐。如果你真的要试别的距离记得同时改PSO适应度函数不要只改一边。最后分享一个我自己的习惯这类聚类项目我会把PSO收敛曲线和最终聚类中心一起放进报告而不是只给一个SSE数值。因为只看数值很难判断算法是被卡在局部最优还是正常收敛。曲线下落很快往往说明初始化和参数不太匹配可能错过了更好的解曲线平滑下降才是比较理想的状态。居民用电行为分析本身是强业务场景算法稳定、结果可解释比单纯追求一个漂亮的聚类指标更重要。后续如果想继续扩展还可以把PSO换成自适应惯性权重、多目标PSO同时优化SSE和轮廓系数或者把聚类结果接入一个简单的用户标签系统这些都是在当前这套代码基础上能低成本展开的方向。