居民用电行为分析这几年被越来越多人捡起来做，原因很简单：智能电表普及之后，营销侧手里攒了一批按15分钟甚至5分钟一条的负荷曲线数据，东西是好东西，但大多数数据躺在数据库里根本没被充分利用。我做这个项目的时候，目标也很直接——用粒子群算法优化Kmeans聚类，对居民用户的用电模式做一次真正能落地的分群，而不是停留在画几张线图、跑个报告交差的程度。
这篇文章我会把完整的思路、算法原理、Matlab代码实现和踩过的坑一次性讲透。内容包括：为什么要用粒子群算法去包一层Kmeans、居民用电数据到底怎么处理才不会出脏结果、粒子怎么编码、适应度函数怎么设计、代码怎么写才能跑得动上万用户，以及聚类出来之后怎么解读、怎么接入实际业务。适合正在做负荷聚类、用户画像、需求响应分析的同学参考，也适合毕业论文想拿这套方法做组合优化方向的读者。
1. 为什么非要用粒子群算法优化Kmeans
1.1 居民用电数据的聚类需求
居民用电数据本质上是一堆高维时间序列。如果按每天96个采样点来组织，一个用户365天就是96乘365维的数据，直接做聚类几乎不可行，所以通常要压缩成“典型日负荷曲线”或“画像特征”。但不管怎么压缩，聚类的目标始终围绕几个实际问题：哪些用户是“上班族晚高峰型”、哪些是“全天候高耗能型”、哪些是“夜间谷电型”、哪些是“冬夏季节性尖峰型”。
这些标签直接关系到需求侧响应、峰谷电价套餐设计、台区线损诊断等业务方向。问题是，居民用户在用电行为上的差异远没有工业用户那么明显，很多曲线长相接近，靠肉眼根本分不出边界，这时候聚类算法就成了刚需。
1.2 传统Kmeans在用电数据上的硬伤
Kmeans的原理不复杂：随机选k个初始中心，迭代更新簇分配和中心坐标，直到收敛。但它在用电数据上翻车是常态，我实测下来主要有三个问题。
第一个是初始中心敏感。Kmeans对初始中心的选择非常敏感，跑十次可能出三四种完全不同的分群结果，这在项目里是致命的——你没法跟业务方解释为什么昨天跑出来是5类、今天变成4类。
第二个是容易陷入局部最优。用电数据的特征空间里，簇的形状往往不是标准球形，有些簇内部密度差异极大，Kmeans这种基于距离均值的算法很容易卡在局部极值点，结果就是某些簇特别大、某些簇只剩几个点。
第三个问题在业务侧暴露得更明显：用户规模一大，Kmeans的稳定性就会急剧下降。我处理过几万户级别的数据，反复随机初始化带来的方差非常吓人。所以项目里需要的不是又一个普通Kmeans，而是一个能保证“每次跑出来结果都差不多”的聚类方案。
1.3 粒子群算法能解决什么
粒子群算法是一种群体智能优化算法，它的思路是让一群“粒子”在解空间里飞行，每个粒子记录自己找到的最优位置（pbest），同时共享群体全局最优位置（gbest），通过速度和位置更新不断逼近最优解。它不依赖梯度信息，对目标函数的连续性要求很低，特别适合用来搜索Kmeans的初始中心。
放到这个场景里，PSO的角色就是帮Kmeans找一组“好”的初始聚类中心。每个粒子代表一组候选中心点，适应度函数用类内距离平方和（SSE）来评价，粒子群进化若干代之后，把全局最优的粒子作为Kmeans的初始中心，再做一次标准Kmeans迭代精调。
这样做的好处有两个层面。第一个是全局搜索能力强，粒子群在整个解空间里撒点，不容易像随机初始化那样扎堆到局部区域，因此能显著降低Kmeans对初始值的敏感性。第二个是结果可复现，只要固定随机种子和参数，PSO-Kmeans每次跑出来的分群基本一致，这在业务交付和数据更新场景里非常重要。
2. 数据处理和特征工程：决定聚类质量的前置环节
2.1 原始用电数据怎么组织才合理
我一开始做这个项目时直接把所有用户每天的96点曲线堆到一起聚类，结果非常惨，聚类出来的簇完全没有可解释性。问题在于：同一个用户工作日和周末的用电曲线完全不同，冬天和夏天也不同，混在一起之后，聚类算法根本分不清这些差异是用户间的还是时间上的。
后来我的做法是先做“场景切片”，把数据按工作日/周末、峰季/平季/谷季切成若干子数据集，每个数据集单独聚类。以工作日为例，先取出某个用户在连续一个月内所有工作日的日负荷曲线，逐点取平均，得到一条“典型工作日负荷曲线”，再按96点或者降维后的特征参与聚类。
还有一个细节是采样点对齐。不同地区的电表冻结时间可能不一样，有的是整点冻结，有的是15分钟偏移。如果直接拿数据算曲线，必须先把所有用户的时间戳统一对齐到同一个时钟序列上，否则后面的特征计算全是歪的。
2.2 特征选择：从原始曲线到画像特征
用96维的曲线直接喂给聚类算法不是不行，但计算开销大、噪声也多。我更推荐先算一组业务可解释的画像特征，再基于这些特征做聚类。表里列的是我项目里常用的特征。
| 特征名称 | 计算方式 | 业务含义 |
|---|---|---|
| 日平均负荷 | 日用电量 / 24 | 整体用电水平 |
| 日负荷率 | 平均负荷 / 最大负荷 | 用电平稳程度 |
| 峰谷差率 | (峰负荷 - 谷负荷) / 峰负荷 | 调峰潜力 |
| 峰时用电占比 | 高峰时段电量 / 日总电量 | 对峰电价的敏感度 |
| 夜间用电占比 | 22:00-06:00电量 / 日总电量 | 储能或蓄冷潜力 |
| 最大负荷出现时刻 | 典型日最大负荷对应时刻 | 生活作息的间接反映 |
| 休息日与工作日差 | (休息日电量 - 工作日电量) / 工作日电量 | 通勤型用户判别 |
这套特征的组合效果明显好于直接用原始曲线，原因是它把维度从96降到了个位数，同时每个维度都有清晰的业务含义。后面聚类出来的每个簇，都可以直接通过这些特征的均值来给用户“画像”，而不是看着一堆曲线去猜。
我自己经常会加一个“温度敏感度”特征，方法是计算用户负荷与当地气温的相关系数。居民用户里空调用电占相当大比例，这个特征能很好地分出“空调主导型”用户。
2.3 归一化与异常值处理
特征归一化这件事看着简单，翻车概率却很高。直接用原始数值做聚类时，日平均负荷的数值范围可能是1到10，而峰谷差率可能是0到1，量纲差异会直接导致欧氏距离被大数值特征主宰，小数值特征直接失去作用。
我试过两种方案：min-max归一化和z-score标准化。在用电特征这个场景里，z-score通常表现更稳，尤其是存在极端高负荷用户时，min-max会把正常用户压到很窄的区间里。z-score对每个特征独立计算均值和方差，能保留分布形态。
异常值处理也要小心。居民用电数据里经常出现“电表停走”“采集失败补零”“错接互感器”等情况，这些会产生一些离谱的零值曲线或尖峰曲线。我的经验是先用规则过滤掉日用电量常年为0或单日最大负荷超过变压器容量上限的用户，再对剩余数据做聚类。规则过滤要放在聚类之前，否则异常用户会自己独成一簇，把所有正常用户挤在一起。
3. PSO-Kmeans的核心原理与关键参数设计
3.1 粒子编码：聚类中心怎么放进粒子
粒子编码是整套方法里最需要想清楚的一步。常见的做法是把k个聚类中心的所有特征值拼接成一个一维向量，作为粒子的位置。
假设聚类数k=4，每个用户用d=7个特征描述，那么粒子长度就是4乘7共28维。粒子在前28维空间里飞行，每一维对应某个聚类中心的某个特征坐标。解码的时候，把粒子的位置向量reshape成4乘7的矩阵，就得到了k个中心点。
还有一个方案是按“样本分配索引”编码，但我不建议用。原因有两个：一是维度跟样本量挂钩，几万用户就是几万维，粒子群在这种超高维空间里几乎不可能有效搜索；二是离散索引空间不连续，粒子更新的连续位置没法直接映射到簇分配。中心点编码是实践中最稳的方案。
3.2 适应度函数怎么设计才算对
适应度函数是粒子群算法的“指挥棒”。如果设计得不好，算法收敛到的东西对Kmeans来说毫无意义。
我用的是类内距离平方和的倒数，即：
SSE = Σ_i Σ_x∈Ci distance(x, c_i)^2 Fitness = 1 / (SSE + eps)
其中eps是一个很小的正数，防止除零。SSE越小，代表簇内样本离中心越近、簇越紧凑，适应度越大。粒子群算法会朝着适应度更大的方向进化。
但只用SSE有一个隐蔽问题：k=4和k=5之间，SSE天然不在一个量级，所以适应度函数不能用来比较不同k的结果，只能在k固定时比较不同粒子。若想让算法自动权衡簇内紧凑度和簇间分离度，我会改用Davies-Bouldin指数（DBI）做适应度，DBI越小越好，同样取倒数。
这里有一个关键点容易被新手忽略：适应度计算内部要完整跑一遍簇分配和中心更新，而不是简单地算粒子位置到所有样本的距离。我的做法是，粒子位置先看作一组中心点，计算所有样本到这些中心的距离，把样本划入最近中心所在簇，然后重新计算每个簇的均值作为新中心，再用新中心重新计算SSE。这样粒子群搜索出来的解，本身就是经过一轮Kmeans修正的局部最优结构。
3.3 算法流程与关键参数
完整流程我放在这里，代码直接照这个逻辑写就行。
PSO-Kmeans核心流程：
- 载入特征矩阵X，确定聚类数k、粒子群规模N、最大迭代次数T。
- 初始化N个粒子，每个粒子的位置从样本中随机抽取k个样本并叠加小幅扰动得到。
- 对每个粒子，调用内部Kmeans一次，计算适应度值。
- 更新每个粒子的个体最优pbest和群体最优gbest。
- 按PSO速度-位置更新公式刷新粒子位置，并对位置做边界约束。
- 重复3到5步直到达到迭代次数。
- 输出gbest对应的中心点，作为初始中心执行标准Kmeans，得到最终簇划分。
参数选择上，我踩了几次坑之后定了一套比较稳的默认值：粒子群规模取30到50，迭代次数取50到100，惯性权重w从0.9线性递减到0.4，学习因子c1和c2都取1.5左右，速度上限vmax设为特征向量范围宽度的20%。
惯性权重这里值得多说一句。w太大，粒子会飞得太野，不容易收敛；w太小，粒子会过早集中到局部区域，失去全局搜索能力。线性递减是一种折中。实测下来，w从0.9降到0.4配合50代迭代，对绝大多数数据集都够用。如果是上万条样本，我建议先随机抽样5000条跑PSO选中心，再全量数据上跑Kmeans，时间能省下一大半，效果几乎不变。
4. Matlab代码实现与核心片段解析
4.1 主程序框架
先贴主程序框架。我把整套逻辑封装成函数，输入是特征矩阵、聚类数、粒子群参数，输出是最终簇标签和聚类中心。
function [label, center, best_fitness] = pso_kmeans(X, k, swarm_size, max_iter) % PSO优化K-means聚类 % 输入: % X - n×d 特征矩阵 % k - 聚类数 % swarm_size - 粒子群规模 % max_iter - 最大迭代次数 % 输出: % label - n×1 聚类标签 % center - k×d 聚类中心 % best_fitness - 全局最优适应度 [n, d] = size(X); lb = min(X, [], 1); ub = max(X, [], 1); range_w = ub - lb; % 参数设置 w_max = 0.9; w_min = 0.4; c1 = 1.5; c2 = 1.5; v_max = 0.2 * range_w; % 初始化粒子群 for i = 1 : swarm_size % 从样本中随机选k个点作为初始中心，叠加小扰动 idx = randperm(n, k); pos{i} = X(idx, :)' + 0.01 * randn(d, k); vel{i} = zeros(d, k); pbest_pos{i} = pos{i}; [~, fitness] = decode_fitness(pos{i}, X); pbest_fit(i) = fitness; end [best_fitness, gbest_idx] = max(pbest_fit); gbest_pos = pbest_pos{gbest_idx}; % 主迭代 for t = 1 : max_iter w = w_max - (w_max - w_min) * t / max_iter; for i = 1 : swarm_size r1 = rand(d, k); r2 = rand(d, k); vel{i} = w * vel{i} + c1 * r1 .* (pbest_pos{i} - pos{i}) ... + c2 * r2 .* (gbest_pos - pos{i}); vel{i} = max(min(vel{i}, v_max), -v_max); pos{i} = pos{i} + vel{i}; pos{i} = max(min(pos{i}, ub'), lb'); [~, fitness] = decode_fitness(pos{i}, X); if fitness > pbest_fit(i) pbest_fit(i) = fitness; pbest_pos{i} = pos{i}; end if fitness > best_fitness best_fitness = fitness; gbest_pos = pos{i}; end end end % 用gbest作为初始中心，跑标准Kmeans精调 center = reshape(gbest_pos, d, k)'; [label, center] = kmeans(X, k, 'Start', center, 'MaxIter', 500); end主程序结构上我用了cell数组来存每个粒子的位置和速度，因为每个粒子实际上是一个d行k列的矩阵。这样后续reshape和计算距离都很方便，比把粒子拉成一行向量更直观，也避免了索引换算的错误。
4.2 适应度计算函数的关键细节
decode_fitness函数是所有粒子的核心评价函数，里面包含了“距离计算、簇分配、中心更新、SSE重算”四个步骤。其中有一点非常关键：计算完簇分配后必须重新计算中心，然后用新中心算SSE，这样粒子的适应度才对应了“经过Kmeans一轮改进后的质量”。
function [center, fitness] = decode_fitness(pos, X) % pos: d×k 矩阵，每列是一个聚类中心 % 返回更新后的中心和适应度 k = size(pos, 2); n = size(X, 1); D = pdist2(X, pos'); [~, label] = min(D, [], 2); % 重新计算各簇中心 center = zeros(size(pos)); for j = 1 : k idx_j = (label == j); if sum(idx_j) == 0 % 空簇保护：随机选一个样本作为该簇中心 center(:, j) = X(randi(n), :); else center(:, j) = mean(X(idx_j, :), 1); end end % 用新中心计算SSE D2 = pdist2(X, center'); min_dist = min(D2, [], 2); SSE = sum(min_dist .^ 2); fitness = 1 / (SSE + 1e-9); end这里必须做空簇保护。粒子群进化过程中，粒子位置在解空间里乱飞，某一维坐标可能会飘到样本稀疏的区域，导致某个中心点周围一个样本都没有，形成空簇。如果不处理，mean函数会对空矩阵求均值然后报错，整个循环卡死。我的处理方法是随机抽一个样本当新中心，虽然对搜索方向略有扰动，但保证了算法不间断运行。
适应度计算是整个算法计算量的主要瓶颈。样本数多的时候，pdist2这一步占了绝大部分时间。一个优化技巧是用max(min_dist)对SSE做归一化，避免样本量差异导致适应度数值过小。另一个更实用的优化是利用Matlab的sum(X.^2, 2)展开欧氏距离计算，避免pdist2的额外内存开销，但代码可读性会下降，我一般保留pdist2版本。
4.3 Kmeans局部精调与结果输出
PSO迭代结束后，全局最优粒子给出的中心已经比较接近全局最优解，但严格来说，PSO在连续空间里的搜索精度不如Kmeans在局部区域的快速收敛能力。所以最后一步一定不能省：以gbest中心为初始值，再调用Matlab自带的kmeans函数做精调。
kmeans函数里有一个很容易踩的坑：必须显式指定'Start'参数为PSO找到的中心，否则Matlab仍会做随机初始化，把PSO的努力全部浪费掉。另外，'MaxIter'建议调大一些到300到500，确保最终收敛。
输出的标签矩阵是一个n行1列的向量，可以直接合并到原始用户信息表里做后续分析。我习惯把聚类中心和每个簇的样本量、典型特征均值一并输出，做成一张“簇画像表”，后面给业务方汇报时非常有用。
5. 聚类结果怎么解读，怎么用到业务里去
5.1 K值选择：别把轮廓系数当唯一标准
K值怎么定是每个做聚类的人都会遇到的老大难问题。我常用的方法有三板斧：肘部法则、轮廓系数、业务可解释性。
肘部法则看SSE随k变化的拐点，操作非常简单：跑k=2到k=10的PSO-Kmeans，记录每次的SSE，画出曲线，找“肘部”位置。轮廓系数可以算每个样本的s值然后取平均，范围在-1到1之间，越接近1代表聚类效果越好。但轮廓系数有个毛病：样本量大时计算很慢，而且它对簇形状比较敏感，不一定能选出业务上最好的k。
我实际项目的判断逻辑是：先看肘部图圈定2到3个候选k，再对比这几个k下每个簇的样本量是否均衡。比如k=5时某个簇占了60%的用户，这个分法通常有问题；如果k=4时各簇占比在15%到35%之间，且每个簇的特征画像能讲出清楚的业务故事，那k=4就是最终选择。算法指标只是参考，业务可解释性和簇的均衡性才是交付时真正卡人的地方。
5.2 用户画像怎么理解
聚类标签出来之后，最重要的一步是把每个簇的特征均值拉出来，生成一张画像表。我以某次项目的实际结果为例，当时k=4，得到四个典型用户群。
第一类簇的用户工作日负荷呈“双峰”形态，早峰出现在7点到9点，晚峰出现在18点到21点，夜间用电占比很低，这类用户基本可以定义为“上班族通勤型”。第二类簇的用户全天负荷平稳，日负荷率非常高，峰谷差率小，多为家中有老人或全天有人活动的家庭。第三类簇的用户夜间用电占比极高，最大负荷出现在22点以后，这类用户要么是蓄能设备用户，要么是夜间充电车主。第四类簇的用户负荷曲线冬夏两季差异巨大，温度敏感度特征值明显偏高，是典型的“空调主导型”。
这些画像不是靠猜的，而是从聚类簇中心的特征值组合里直接读出来的。比如第三类的夜间用电占比均值达到0.62，峰谷差率为0.18，这个组合在业务上基本能锁定“夜间用电偏好的用户”。
5.3 从聚类到落地的方向
聚类结果本身不是终点，能接上的业务方向很多，我列几个在项目实践中验证过的路径。
第一个是需求侧响应。夜间用电占比高的用户群，适合作为蓄热式电采暖或电动汽车有序充电的潜力用户；空调主导型用户群，是空调负荷调控的主要对象。第二个是峰谷电价套餐设计。上班族通勤型用户在工作日晚间时段有明确高峰，推送“晚峰折扣”类套餐的接受度通常更高。第三个是台区线损诊断。某台区如果“全天候高耗能型”用户占比异常高，可以排查是否存在窃电或表计异常。还有一个方向是分布式光伏和储能的配置。白天用电占比高的用户，自发自用比例更高，配套储能的经济性完全不一样。
6. 常见问题和排查技巧实录
6.1 聚类结果每次跑都不一样
这是PSO-Kmeans被问得最多的问题。虽然整体稳定性远好于随机初始化的Kmeans，但粒子群初始化本身有随机性，某些特征分布特殊的数据集上仍可能出现差异。
我的处理办法是规定随机种子。Matlab里用rng(2024)把随机数流固定住，这样每次跑的结果完全一致。如果项目方要求更高，就连续跑20次，取SSE最小的一次作为最终结果，并同时记录20次结果中簇标签的一致性比率。一致性比率低于70%时，我会回头检查特征工程是否存在问题，而不是盲目增加粒子数。
6.2 粒子群早熟收敛
某些数据集上，粒子群在十几代之后就全部收敛到了同一个位置，gbest不再变化。这通常意味着w衰减太快，或者粒子初始化范围太窄。
我试过有效的对策有两个。一个是把w的衰减从线性改成非线性，前期保持较高的惯性权重更久，比如用w = w_min + (w_max - w_min) * exp(-t / max_iter)这种指数衰减方式。另一个是对gbest加扰动，每五代把全局最优粒子的位置叠加一个高斯噪声，相当于给算法一次“跳出局部”的机会。实测下来，指数衰减w对绝大多数用电数据集效果更好。
6.3 数据规模大导致计算时间爆炸
几千样本时PSO-Kmeans跑起来毫无压力，但到了几万甚至几十万用户时，适应度计算会让人等到崩溃。每次适应度计算都要做一次全量距离矩阵运算，50个粒子乘100代迭代，这个计算量是灾难级的。
我的经验是分两步走：第一步随机抽样5000条样本跑完整的PSO搜索，得到一个较优的中心点集；第二步把这个中心点集作为全量数据的Kmeans初始中心，直接跑标准Kmeans。这种方法在理论上有依据，因为聚类中心本质上是由数据密度分布决定的，采样量足够大时中心点的估计误差可以忽略。实践中，我用1%样本量跑PSO得到的中心，在全部数据上复跑Kmeans，聚类结果与全量PSO-Kmeans的一致率在95%以上，耗时却从半小时降到了两分钟。
6.4 特征量纲和季节扰动导致结果漂移
如果用户跨了不同季节的用电数据，直接聚类很容易把“季节差异”误判成“用户差异”。一个所有季节都有空调的用户，在夏天会被分到“空调主导型”，到了冬天可能就被分到“平稳型”，导致标签不稳定。
我的做法是先把用户数据按季节切片，分别聚类，再综合判断。比如夏季聚类能分出空调主导型，冬季聚类能分出采暖型，两者交叉后就能识别出“全年高耗能用户”。也不要一次性用12个月数据做同类特征，需要先把工作日和周末分开。
6.5 代码运行时报错的一些小问题
Matlab代码在实际跑起来的时候有几个高频报错点，提前说一下可以省很多排查时间。
第一个是pdist2的输入维度不对。粒子位置是d乘k矩阵，pdist2要求第二个参数是样本乘特征矩阵，所以必须转置为k乘d。我多次在转置上栽跟头，现在每次写到这里都会刻意检查维度。第二个是kmeans的Start参数格式问题，它接受k乘d矩阵，因此reshape(gbest_pos, d, k)之后需要转置。第三个是版本兼容性，老版本Matlab不支持某些写法，建议用R2016b以上版本跑pdist2和randperm的k参数用法。
还有一个小技巧：把整个PSO迭代写成for t = 1 : max_iter后，建议加一个if mod(t, 10) == 0，fprintf('Iter %d, fitness %.4f\n', t, best_fitness); end，实时观察收敛过程。如果30代以内适应度就不动了，说明粒子群提前收敛，需要调整参数；如果一直到100代还在缓慢上升，说明迭代次数可以加大。
7. 后续还能怎么扩展
这个项目做完之后，我个人觉得最有意思的是把聚类结果跟用电量预测、异常检测串起来。聚类的簇标签本身是很好的监督信息，可以作为特征输入到后续的预测模型里。比如先聚类得到用户类型，再对每类用户分别训练负荷预测模型，比全局一个模型准确率高出一截。也可以用簇分配结果做“行为漂移检测”，连续监测用户每个月的簇标签是否发生变化，一旦用户从夜间型跳到全天型，说明家庭用电结构发生了重大变化，可能新增了大功率设备，这个信号对营销推送和台区管理都有用。
代码方面，如果想把性能再往上推一步，可以把粒子群的速度更新改成向量化写法，摆脱for循环。但Matlab的for循环在粒子群这种场景里其实没那么慢，除非粒子数特别大。真正值得优化的点是适应度函数内部的距离计算，可以预先对特征矩阵做归一化矩阵乘法，减少重复计算。这部分改动比较复杂，我一般是等数据规模真的撑不住了才去动它。
我在实际使用中的一个体会是：粒子群算法和Kmeans的组合，不要神化，也不要忽视。它解决的核心问题是Kmeans初始值敏感和局部最优，但它解决不了特征工程没做好的问题。很多人在项目里跑出来说“PSO-Kmeans效果也不好”，多数情况下是原始数据没有切片、特征没选好、或者k值没找对，算法本身反而没多大毛病。把数据处理做扎实了，再上这套优化方法，结果不会差。
