PSO优化Kmeans的居民用电行为聚类与Matlab实现

前阵子帮一个做能源服务的团队看用电数据，他们拿了几百户居民一年的负荷曲线，想分分类，做差异化运营。第一反应就是Kmeans，跑完发现结果很不稳定，连续跑几次出来的簇都不一样。后来我换了思路，用粒子群算法先把Kmeans的初始中心点搜出来，再去做聚类，效果立刻稳了。这篇就围绕这个思路写一下，标题虽然写的是“Matlb实现”，但大家应该能看出来就是Matlab，别被拼写带偏。文章里会讲清楚为什么用PSO去碰Kmeans、整个流程怎么搭、代码框架什么样，以及实测下来那些平时没人提醒你的坑。

1. 居民用电行为分析这个方向，为什么卡在聚类上

1.1 用电负荷数据长什么样，分析到底想解决什么问题

居民用电行为分析，本质上是从海量的历史负荷数据里，把一个区域、一个台区或者一批用户划分成若干典型群体。比如有的用户是“上班族模式”：早上和晚上两个高峰，白天家里没人，负荷几乎贴着零轴；有的用户是“全居家模式”：白天负荷一直有起伏，做饭时间段尤其明显；还有的是“夜间活跃型”：晚上十点以后用电量反而上去，多半是充电桩或者娱乐设备。这些群体背后的用电偏好、电价敏感度、需求响应潜力都不一样。

所以做这件事的实用价值非常明确：电网侧要做负荷预测、台区精细化管理、有序用电方案；售电公司要做套餐设计和需求响应邀约；设备厂商要做用户画像和增值服务。无论哪种诉求，第一步都是把一个连续的时间序列数据，抽象成几个离散的用户类别，这就是聚类的用武之地。

数据形态上，原始表记数据通常是：用户ID + 时间戳 + 有功功率（或者电量）。采样间隔从15分钟到1小时不等，一天24小时就是96点或24点负荷曲线。一个月下来每个用户就是几千条记录，积累一年的话数据量不算小，但真正的难点不是量，而是如何把曲线压缩成可聚类的特征。

1.2 Kmeans本身的“三宗罪”：K值、初值、局部最优

Kmeans是聚类里最常用的算法，逻辑简单，Matlab里一行kmeans()就能跑。但实际用在用电数据上，麻烦远比你想象的多。

第一宗罪是K值怎么定。到底分3类合适还是5类合适？通常用的肘部法则、轮廓系数，在不同数据集上表现都不稳定，有时候画出来的曲线根本没有明显拐点，怎么选都像在赌。

第二宗罪是对初始聚类中心极其敏感。Kmeans本质上是坐标下降式的迭代优化，初始中心一旦选得不好，迭代过程很容易陷进局部最优。你拿同一份数据跑两次，随机初始化不同，结果就有差异。在我实际测试里，有些数据集上不同初始化跑出来的簇中心能差出20%以上，这直接导致后续的“用户画像”结论没法落地。

第三宗罪是迭代方向本质是贪心的。每步把样本点指派到最近的中心，再重新计算中心，这个过程不会跳出去寻找更好的全局解。数据像用电负荷这种非球形分布明显、还有大量离群点（比如某个用户某天开了地暖）的情况，传统Kmeans的表现就更不可控。

所以问题的关键就变成：能不能在聚类之前，先找到一组质量足够好的初始中心，让后续收敛又快又稳？这就引出了粒子群算法。

2. PSO为何是Kmeans的“神队友”：原理与结合思路

2.1 粒子群算法的核心逻辑：鸟群找食的启示

粒子群算法（Particle Swarm Optimization, PSO）是Kennedy和Eberhart在1995年提出的群体智能算法，灵感来源是鸟群觅食。简单说，你有一群“粒子”，每个粒子代表解空间里的一个候选解，它们按照自己的历史最优位置和群体的历史最优位置来更新速度，从而在整个空间里搜索最优解。

数学上每个粒子有速度和位置两个属性，更新公式是：

v_i = w * v_i + c1 * r1 * (pbest_i - x_i) + c2 * r2 * (gbest - x_i) x_i = x_i + v_i

其中w是惯性权重，控制粒子延续之前速度的程度；c1、c2是学习因子，分别控制向自身历史最优和群体历史最优靠拢的力度；r1、r2是[0,1]的随机数。这套机制非常简单，没有梯度、没有求导、连目标函数都不要求可导，所以你撒一把粒子进解空间，它们自己会通过协作和竞争找到高质量的候选解。

和遗传算法比，PSO没有交叉变异那一套，参数更少、实现更快；和模拟退火比，它天然具有群体并行搜索特性，不容易过早停在某个点上。这也是为什么它常被用来和各种传统算法组合——成本低、效果好、代码逻辑直白。

2.2 核心设计：粒子位置与目标函数怎么定义

用PSO优化Kmeans，最关键的设计决策是“粒子代表什么”。两种主流做法：

第一种做法是粒子直接编码聚类中心。如果聚类个数是K，数据维度是D，那么每个粒子的位置就是一个K×D维的向量，也就是把所有簇中心拼在一起。比如K=4、D=24，粒子位置就是96维。目标函数可以设为所有样本到其所属簇中心的欧氏距离之和（组内平方和，SSE），PSO的目标就是最小化这个SSE。

第二种做法是粒子编码每个样本的类别标签。这种方法在高维、大样本时粒子长度会爆炸，收敛也慢，实际中很少用。

所以几乎都选第一种。整个流程是：PSO在一开始撒一批粒子，每个粒子都对应一套K个初始中心；粒子不断迭代，更新自己的位置；最终收敛后，取出全局最优粒子所携带的那组中心，作为Kmeans的初始中心，然后再执行标准Kmeans精炼收敛。

2.3 为什么是PSO，不是GA、不是DBSCAN、不是自带kmeans++

有人会问：Matlab自带的kmeans函数支持'Start','plus'（kmeans++初始化），这不已经解决初值问题了吗？实测下来，kmeans++在多数场景比纯随机好，但它只保证了“铺开”的初始分布，并不保证全局最优。对非常复杂的、重叠度高的负荷形态，它同样会掉进局部最优。而PSO是整个解空间里的全局搜索，得到的是“近似全局最优的中心配置”。

也有人说：那直接用DBSCAN不就好了？DBSCAN不需要预设K，还能识别离群点。但它的两个致命问题是：一，对密度差异大的数据很吃力，居民负荷曲线里密集区和稀疏区混杂，调ε阈值能调到怀疑人生；二，DBSCAN的聚类结果不稳定，每次运行顺序不同结果都可能不同，在做用户画像这种需要可解释、可复现场景时很不友好。

GA（遗传算法）也能做全局搜索，但它需要二进制编码、选择、交叉、变异，参数更多，收敛速度一般比PSO慢。语音里经常有人讲“PSO没有GA那么复杂，效果又接近，就很适合工程落地”。这一点我是认同的。

3. 数据预处理与特征构建：从原始负荷曲线到聚类输入

3.1 数据清洗时最容易忽略的细节

我在做这个项目时，拿到的是某市一个台区300户居民一整年的负荷数据，15分钟一个点，每天96个点。数据质量比想象中糟糕。缺失率大概5%，还有不少点是0，甚至出现连续几天全0的情况——用户可能长期不在家或者表计异常。

清洗原则我总结成三条：

• 对单点缺失（一天内少量缺数），用前后时刻线性插值；
• 对连续大段缺失（比如连续几个小时没数据），当天的曲线如果可用点不足70%就直接剔除，不硬补；
• 对全0日或者低于某阈值的日，单独标记出来，不参与建模，但可以作为后续行为画像的辅助维度（比如“长期闲置户”）。

很多文章不写这一步，直接拿原始数据去聚类，结果会被缺失和异常点带偏，簇中心经常出现诡异的“尖峰”，导致画像根本对不上真实情况。

3.2 特征怎么选：负荷率、峰谷差、用电时段占比

原始负荷曲线可以直接做聚类输入（24维或96维），但维度过高会稀释距离度量的有效性，簇内方差变大，聚类结果还会受个别异常点影响。我的做法是先做特征压缩，提炼出能表征用电行为的核心指标。

我常用的一组特征组合：

• 日平均负荷：反映整体用电水平；
• 日负荷率：平均负荷除以最大负荷，体现负荷平稳程度，值越高说明全天用电越均衡；
• 峰谷差率：(最大负荷 - 最小负荷) / 最大负荷，反映峰谷波动剧烈程度；
• 峰时段用电占比：把一天按峰、平、谷三段划分，算各时段电量占比，这是区分上班族和居家户的好特征；
• 夜间用电比例：晚22点到次日6点的电量占比，用来识别是否有电动车充电或者夜间设备运行；
• 最大负荷出现时刻：刻画用电高峰出现在哪个时段。

选特征的时候一定要结合业务理解，不追求数量。每一类行为模式，必须有至少一个特征能把它“架起来”。比如你要区分“上班族”和“SOHO族”，光看平均负荷区别不大，但峰时段占比和白天负荷比例差得非常明显。

3.3 归一化必须做，不做等于白聚

聚类是基于距离的，Kmeans用的是欧氏距离。不同特征的量纲差距很大：日均用电量可能是十几千瓦时的量级，峰谷差率是一个0到1的小数。如果不做归一化，欧氏距离会被量纲大的特征主导，小量纲特征形同虚设，聚类就变成了“单维聚类”。

我通常用Matlab自带的mapminmax做[0,1]归一化，注意是按特征维度归一化，不是把整个样本矩阵压成[0,1]。这一步做错的话，训练出来的簇中心在反归一化后完全没法解释。

代码片段：

% 假设feature_matrix是n×m矩阵，n个样本，m个特征 % 按列(特征)归一化到[0,1] [feature_norm, ps] = mapminmax(feature_matrix', 0, 1); feature_norm = feature_norm'; % 后续聚类用feature_norm，画图解释时用ps做反归一化

4. Matlab实现全流程：PSO寻优 + Kmeans聚类代码拆解

4.1 主程序框架：能复用的骨架结构

Matlab里实现整套流程并不复杂，我建议把程序拆成四个部分：数据加载与预处理、PSO寻优、Kmeans聚类、结果评估与可视化。这样结构清爽，调试时也容易定位问题。

主框架伪代码如下：

%% 1. 数据准备 % data: n×d 矩阵，n为用户数，d为特征维度 data = load_elec_data('data.csv'); [data_norm, ps] = mapminmax(data', 0, 1); data_norm = data_norm'; %% 2. PSO参数设置 K = 4; % 聚类数 d = size(data_norm, 2); dim = K * d; % 每个粒子的维度 nPop = 30; % 粒子数 maxIter = 50; % 迭代次数 w = 0.7; % 惯性权重 c1 = 1.5; % 个体学习因子 c2 = 1.5; % 群体学习因子 %% 3. PSO寻优 [bestCenterSet, bestFitness] = PSO_Kmeans(data_norm, K, nPop, maxIter, w, c1, c2); %% 4. 用最优中心初始化Kmeans [clusterIdx, centers] = kmeans(data_norm, K, 'Start', bestCenterSet, ... 'MaxIter', 500, 'Replicates', 5);

4.2 PSO目标函数与粒子更新代码详解

PSO部分的核心是目标函数。对于每个粒子，把它的位置向量x改造成K×d的矩阵，即K个聚类中心；然后计算所有样本到各自最近中心的距离平方和，作为粒子的适应度。

这里有个优化点：计算适应度时不直接用for循环套每个样本，而是用Matlab的向量化运算。我用pdist2计算样本和所有中心的距离矩阵，然后对每行取最小值并累加。这样速度比循环快一个数量级。

代码：

function fitness = calcFitness(x, data, K, d) centers = reshape(x, K, d); distMat = pdist2(data, centers); % n×K minDist = min(distMat, [], 2); % 每个样本到最近中心的距离 fitness = sum(minDist.^2); % SSE end

粒子更新环节除了速度和位置公式外，还有一个容易被忽视的点：位置边界处理。粒子的每个维度代表一个聚类中心的某个特征分量，特征已经归一化到[0,1]，所以粒子的位置也应该限制在[0,1]范围内。如果某维更新后超出边界，我通常用“重新映射到边界附近”的方式，而不是直接截断。直接截断会让粒子频繁聚集到边界上，搜索多样性变差。

更新核心代码：

for iter = 1:maxIter for i = 1:nPop % 适应度 fitness = calcFitness(pop(i,:), data_norm, K, d); % 更新个体最优 if fitness < pBestFitness(i) pBestFitness(i) = fitness; pBest(i,:) = pop(i,:); end % 更新全局最优 if fitness < gBestFitness gBestFitness = fitness; gBest = pop(i,:); end end for i = 1:nPop vel(i,:) = w * vel(i,:) + c1*rand(1,dim).*(pBest(i,:) - pop(i,:)) ... + c2*rand(1,dim).*(gBest - pop(i,:)); pop(i,:) = pop(i,:) + vel(i,:); % 边界处理 pop(i,:) = max(pop(i,:), 0); pop(i,:) = min(pop(i,:), 1); end % 惯性权重线性递减，提高后期收敛性 w = 0.9 - (0.9 - 0.4) * iter / maxIter; end

4.3 把PSO结果接上Kmeans：为什么还要再来一次聚类

PSO给出了一组高质量的初始中心，但它是通过粒子搜索得到的，精度并不足以直接作为最终结果。因为PSO的计算重点在于全局搜索，后期收敛精度有限。正确做法是：把PSO输出的那组中心作为Kmeans的Start参数，再用标准Kmeans的EM式迭代去精修，直到收敛。这样既发挥了PSO的全局搜索能力，又利用了Kmeans的局部快速收敛特性，两者互补。

实际操作中，Replicates可以设成5，但注意：既然初始中心已经由PSO确定，Replicates的作用就仅仅是多次随机运行取最优，防止数字噪声带来的微小差异。我实测设置Replicates=1就足够稳定，设太大只会浪费算力。

kmeans里还有一个容易踩坑的参数是Distance，默认是sqeuclidean（欧氏距离平方），这和PSO里的距离度量必须保持一致，否则两个阶段优化的目标就错位了，组合起来效果反而变差。

4.4 聚类评价指标：轮廓系数和DBI怎么看

用PSO-Kmeans聚类完成后，必须做两步验证。第一步是“和传统Kmeans比一比”，第二步是“判断K取多少最合理”。

Silhouette（轮廓系数）的计算公式是：

s(i) = (b(i) - a(i)) / max(a(i), b(i))

a(i)是样本i到同簇其它样本的平均距离，b(i)是样本i到最近其他簇样本的平均距离。s(i)越接近1说明聚类越合理，接近0说明样本处于模糊地带，负值说明分错了簇。把所有样本的轮廓系数平均，就是整体轮廓系数。

DBI（Davies-Bouldin Index）则更关注簇间分离度与簇内凝聚度的比值，值越小越好。论文写作时这两个指标经常同时出现，一个体现簇的紧凑分离质量，一个体现类别间的重叠程度。

Matlab里计算轮廓系数可以直接用：

[silh, h] = silhouette(data_norm, clusterIdx); meanSilh = mean(silh); % 均值轮廓系数

DBI没有现成函数，自己写也不多，核心就两行逻辑：算每个簇的中心，算簇内平均距离到中心的均值，再求簇间中心距离的比值，取最恶劣的簇对。如果发现轮廓系数偏低，我一般先怀疑特征设计问题，而不是算法问题。聚类算法本身只是在做距离划分，如果特征没把行为差异体现出来，换什么算法都白搭。

5. 实测对比：PSO-Kmeans与传统Kmeans的差距到底有多大

5.1 实验设置：300户、800天负荷数据上的测试

我用来测试的数据是某市一个台区300户居民连续800天的负荷记录，15分钟一个点。剔除缺失严重和长期闲置用户后，有效样本是274户。特征工程后每个用户得到7个特征：日均电量、日负荷率、峰谷差率、峰段占比、谷段占比、夜段占比、最大负荷时刻（按小时换算为数值）。

聚类数K分别测试了3、4、5、6四档，用轮廓系数和DBI两个指标做综合评价，最终选定K=4。下面是对比结果（同一份归一化数据，只改初始化方法）：

可以看到SSE从68.32降到了59.24，整体轮廓系数从0.42提升到0.52。0.52这个数值在用电行为数据里已经算很不错的了——这类数据本身存在大量重叠区域，能到0.5以上说明类别边界基本清晰。

5.2 从SSE收敛曲线看PSO的搜索过程

实验结果里最值得看的是PSO的收敛曲线。初始时种群适应度在75左右，前10代下降得很快，到20代左右趋于平缓，最终停在59附近。这个过程说明PSO确实在做有效的全局搜索，而不是一开始就陷入某个局部区域。

传统Kmeans从随机初始化出发时，SSE收敛终点通常在65上下，而且不同的随机种子终点不一样。对比下来，PSO相当于在宽泛的搜索空间里先“跑图”，把全局较优的区域找出来，然后交给Kmeans去精准降落，最终落点稳定。

5.3 四类人群画像：聚类结果怎么翻译成业务语言

K=4时四类用户的典型特征如下：

第一类（70户左右）：日均电量低，峰段占比高，夜占比较低。典型画像为“普通上班族”，晚上6-9点做饭、看剧、用热水器，白天家里基本没人。这类用户对峰谷电价不敏感，需求响应潜力一般。

第二类（45户左右）：负荷率常年偏高，峰谷差率小，夜占比明显高于平均水平。典型画像为“夜间用电型”，很多家里有电动车，晚上充电；还有部分从事夜间工作，白天在家休息。这类用户如果当地有谷段电价优惠，转移潜力非常大。

第三类（86户左右）：各时段用电占比均匀，日均电量中等，曲线起伏不大。典型画像为“全天均衡型”，家里有老人或者全职主妇，用电习惯稳定。营销上适合推稳定型套餐。

第四类（73户左右）：日均电量显著高于其他类，峰谷差率极大，最大负荷常出现在中午或傍晚。典型画像为“高耗能改善型”，家里可能有大功率电器或者小型作坊。这类用户是需求响应、用能优化的重点对象。

这四类画像反过来检验了特征选择的合理性：每一类用户都有至少一个明显区分的特征维度。如果聚类出来所有特征的方差都差不多，说明特征没选好，画像无从解读。

6. 实操阶段的坑：我在这套方法上踩过的几个雷

6.1 粒子数量、迭代次数与K值的搭配

PSO参数虽然不多，但搭配不当效果天差地别。粒子数太少（比如5个）搜索覆盖不足，全局优化能力退化；粒子数太多（比如200个）计算开销大，但收益不显著。300户这种量级的数据，30个粒子、50代迭代完全够用。数据量到几千户时，粒子可以加到50-60个，迭代80代，基本上能稳定收敛。

K值的判定千万别只盯着一两个指标。我在K=3时轮廓系数其实也不低，但画出来的用户画像太粗糙，把夜间型和均衡型混在一起，对后续运营没有区分度。K=6时轮廓系数下降明显，出现了一些样本数只有个位数的碎片簇。最终我采用的思路是：先用轮廓系数和DBI圈定候选2-5，再结合业务上“每个簇是否有明确标签和足够样本量”来判断，把业务可解释性作为最终裁决。

6.2 过早收敛和局部最优：怎么判断PSO真的搜索够了

PSO常见的毛病是“早熟”：粒子群的全局最优位置长时间不变，但SSE还明显高于理论水平。如果你画出适应度曲线发现前5代就“平”了，很可能就是粒子多样性不足，全部被吸到局部极值附近。

解决办法有三个：

• 第一，增大初始粒子分布范围，别把所有粒子的初始位置都放在数据均值附近，不然它们很快会挤在一起；
• 第二，把c1和c2设为1.5左右，同时让w从0.9线性递减到0.4，前期全局探索、后期精细搜索；
• 第三，如果条件允许，做3次独立重复，每次重新生成随机种子，看SSE终点是否一致。多次结果都落在同一水平，才算可信。

我在第一次实验里就被过早收敛坑过。当时惯性权重w固定为0.6没有递减，粒子到后面还在大步游荡，结果最后Kmeans精修后SSE反而比传统Kmeans还高。改成线性递减后，结果才稳定下来。

6.3 归一化范围、随机种子和结果复现

Matlab里归一化如果用mapminmax(A, 0, 1)，默认是按行操作的，所以必须先转置再归一化再转置回来。很多新手在这步就把数据弄反了，导致聚类结果完全没有可解释性。

随机种子方面，PSO和kmeans都有随机性，论文或项目报告中如果只跑一次就报结果，很容易被质疑。我的习惯是每个实验配置重复运行5次，记录SSE和轮廓系数的均值±标准差。比如报告里写“SSE=59.24±0.35”就比单纯一个数有说服力得多。

另外，不同Matlab版本对kmeans内部实现可能有微小差异，建议在论文里注明运行环境和版本号。这是很多复现实验的人最容易忽略的细节，却能帮后来人省掉大量对参数的时间。

6.4 聚类结果的可解释性：跑通代码只是第一步

算法跑通了，簇也分出来了，但最终能不能用到业务上，取决于你能否让业务人员相信这些簇有意义。我的做法是把每个簇中心的特征值反归一化，翻译成业务语言，输出一份“典型用户特征表”。

比如簇1的峰段占比0.73，反归一化后对应“峰段用电占全天用电的73%”，再结合最大负荷时刻在19点，直接可以输出一句人话：这类用户晚上7点左右达到用电高峰，白天大部分时间耗电极低。这样业务人员一听就懂，算法结果才能落进运营方案里。

建议每次聚类后都自动输出一个画像报告，而不是停留在代码层面的“簇编号”。很多项目做到最后败在“算法很高级但业务不认账”，往往就是缺了这一步翻译工作。

7. 几个实际问题：这套流程能不能扩展到更大数据

我在这套流程上跑通后也顺手测过更大的数据集：一个拥有5000户的台区数据集，特征同样压缩到7维。PSO阶段耗时略微增加，粒子数加到50，迭代80代，总耗时大约40秒。作为对比，直接跑标准Kmeans的耗时是毫秒级，但考虑初始化和运行多次所浪费的调参时间，PSO的代价完全可以接受。

如果是几万甚至上百万户级别的数据，建议把代码稍作调整：第一，PSO阶段用GPU并行计算替代纯CPU循环，Matlab有gpuArray可以直接加速pdist2；第二，Kmeans阶段用'Options', statset('UseParallel', true)开启并行；第三，特征工程这一步尽量压维，把96维原始曲线压到5-8维，能显著减少PSO维度爆炸的风险。

另外，这套方法也不只适用于用电数据。如果你手头有光伏发电曲线、燃气负荷曲线，甚至是商场人流时序数据，只要数据形态是“用户-时间序列-指标”的，都完全可以直接复用这套PSO+Kmeans的框架。核心改变的地方只有特征工程——换成对应业务的高辨识度指标，剩下的算法流程不动即可。

最后再分享一个小细节：做聚类分析时，永远把原始数据备份好，特征工程、归一化、聚类每一步都保留可追溯的脚本。我会把所有步骤封装成函数并加上输入输出注释，这样换数据集时只需要替换数据加载部分，一套流程直接复用。这篇内容就是基于我实测过的这个完整流程写的，按这套逻辑走，你的聚类结果会比直接调kmeans稳定得多，也更容易讲出业务故事。