企业官网建设流程全解析

居民用电行为分析这个方向，我最早用的是最朴素的Kmeans加肘部法则，跑出来结果也还行，但过几天换一批数据，或者稍微调一下输入特征，聚类就开始“漂移”。后来我在Matlab里把粒子群算法和Kmeans聚类组合到一起，用粒子群去搜初始聚类中心，再交给Kmeans做局部精调，才把这类问题真正压下去。这篇就完整整理一下我做“基于粒子群算法优化Kmeans聚类的居民用电行为分析”这个项目时的思路、公式、代码和踩过的坑，重点偏向Matlab代码实现，希望对正在做负荷聚类、用户画像或者智能用电分析的朋友有帮助。

1. 居民用电行为分析为什么要用粒子群算法优化Kmeans聚类

1.1 居民用电行为分析想解决什么问题

居民用电行为分析，本质上是从海量智能电表数据里提取出用户的用电模式，再做分群和画像。常见的业务场景包括：需求响应目标用户筛选、分时电价套餐推荐、窃电嫌疑排查、配电网台区负荷预测等。很多研究会把每户的日负荷曲线作为样本，特征可能是96点功率序列、24小时电量，也可能是峰谷电量、最大负荷时刻、周末与工作日差异等统计量。

这个问题的难点在于，真实居民用电行为非常杂。同样都是晚高峰型，有人19点做饭用电，有人21点才开空调；同样都是低谷型，有人是上班族白天不在家，有人是老人全天待机。曲线的相似性不是一眼能看穿的。聚类算法就是要把这种“模糊的行为共性”量化出来，把用户分成若干个内部相似、外部差异明显的群体。

1.2 传统Kmeans聚类为什么不够用

Kmeans是大家用得最多的聚类方法，逻辑也最简单：随机选K个中心，然后反复迭代，直到类内距离平方和最小。但这个算法有两个天生短板，在居民负荷数据上尤其致命。

第一，初始中心敏感。Kmeans的迭代只是局部搜索，最后收敛到哪个局部最优，很大程度上取决于一开始那K个中心怎么挑。负荷曲线维度高、样本量大，目标函数存在大量局部极值，所以同一份数据跑二十遍，经常出现三五套不同的聚类结果。

第二，K值不好定。肘部法则给出的“肘点”经常不明显，尤其是负荷曲线之间过渡平滑时，你很难说K=5比K=4好多少。再加上数据里的异常值、缺失值，Kmeans的中心会被拉偏，聚类稳定性更差。

1.3 PSO与Kmeans组合的定位

粒子群算法（PSO）是一种全局随机搜索算法，它不依赖梯度信息，适合处理Kmeans初始中心这种“多点组合的最优化问题”。把PSO和Kmeans组合起来，定位很清楚：PSO负责在大范围搜索里找到一组较好的聚类中心，Kmeans负责把这些中心进一步迭代精调。相当于先派侦察兵摸清地形，再由主力部队精确占领，比闭着眼睛随机扎营靠谱得多。

这套组合在Matlab里实现成本也不高，PSO本身代码量不大，Kmeans可以直接调用内置函数，两者之间只需要打通“粒子位置→聚类中心”的接口就行。下面我把原理和实现一步步拆开说。

2. 粒子群算法优化Kmeans的核心原理与关键参数

2.1 PSO算法的基本思想与公式

粒子群算法的灵感来自鸟群觅食。每个粒子代表一个候选解，它不停的移动，移动速度受两个记忆影响：一个是自己历史最优位置pBest，一个是整个群体的历史最优位置gBest。通俗说，每个粒子既坚持自己的经验，又参考同伴发现的好位置，两种方向加权组合，就构成了下一轮速度。

速度更新公式是：

[ v_{i+1} = w \cdot v_i + c_1 \cdot r_1 \cdot (pBest - x_i) + c_2 \cdot r_2 \cdot (gBest - x_i) ]

位置更新公式是：

[ x_{i+1} = x_i + v_{i+1} ]

其中w是惯性权重，c1、c2是学习因子，r1、r2是0到1之间的随机数。惯性权重w越大，粒子越容易保持原有速度，全局探索能力强；w越小，粒子越容易被个体和群体最优吸引，局部开发能力强。我习惯让w从0.9线性递减到0.4，前期探索，后期收敛，这样比固定w更稳。

2.2 适应度函数怎么设计

在PSO-Kmeans里，粒子的“好坏”必须量化成一个适应度值。最自然的目标函数是所有样本到各自所属聚类中心的距离平方和，也就是SSE，也叫簇内误差平方和。

[ SSE = \sum_{i=1}^{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维的日负荷曲线，K=5，那么每个粒子就是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。如果你设成大于1，Matlab会自动以多次随机初始化来找更优结果，那你给的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精调。两组统一K=5，距离度量都是平方欧氏距离。这里要说明一下，K=5不是拍脑袋，是我先用轮廓系数扫了K从2到10之后定的。两组用完全相同的K和特征，才有可比性。

4.2 评价指标对比

单一指标容易骗人，我每次报告都会给三个指标：SSE、轮廓系数（Silhouette Coefficient）和Calinski-Harabasz指数。

轮廓系数综合了内聚度和分离度，范围从-1到1，越大越好。CH指数是簇间散布与簇内散布的比值，也是越大越好。下面是我项目里的一组代表性结果：

可以看到，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收敛慢最常见的两个原因，一个是粒子维度太高，一个是初始粒子都集中在可行域边缘。居民负荷特征如果做了很多维统计量，再加上K=8，粒子维度可能超过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:K+1); 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和轮廓系数，或者把聚类结果接入一个简单的用户标签系统，这些都是在当前这套代码基础上能低成本展开的方向。