粒子群算法优化Kmeans聚类：居民用电行为分析Matlab实战

去年我在做居民用电行为分析时，用Kmeans聚类用户负荷曲线，最头疼的就是每次跑出来的结果都不一样。同样的数据，换一次初始中心就得到一批完全不同的用户分群，跟业务部门对需求响应方案的时候解释成本特别高。后来我用粒子群算法去优化Kmeans的聚类中心选择，整个流程在Matlab里跑通，不仅聚类结果稳定多了，还顺手把用户画像做成了能让非技术背景的同事看懂的东西。这篇文章就围绕这个思路，把从数据预处理、特征构造、PSO-Kmeans融合设计到Matlab代码实现的关键环节完整复盘一遍。

无论你是做电力数据分析、负荷预测，还是在其他行业碰到“聚类不稳定”的问题，这套方案都有直接的参考价值。我会把粒子群算法原理、参数设置、编码方式、适应度函数设计、踩坑记录都写清楚，尽量做到可以直接照着复现。

1. 居民用电行为聚类为什么“不好聚”：Kmeans的先天局限

1.1 用电行为数据的“高维、高噪、强混合”特性

智能电表采集的居民负荷数据，最常见的是15分钟一个点，一天96个点，一个月就是2880维。如果你拿这样的原始数据直接丢给Kmeans，理论上可以跑，但实际效果往往一塌糊涂。

原因有两个层面。第一，维度太高，距离度量失效。在高维空间里，欧氏距离会变得“扁平”，所有样本彼此之间的差距趋同，聚类算法很难区分出真正的结构。第二，居民用电行为本身就不是干净的簇状分布。一个家庭可能是上班族晚上回来用电，也可能是老人全天在家慢悠悠用电，还有可能装了电动车充电桩，深夜才启动大功率充电。这些模式彼此叠加、混合，边界非常模糊。

所以做居民用电聚类，第一步必须做特征工程，把2880维的时序曲线压缩成几十个有业务含义的特征维度，而不是直接对着原始曲线聚类。这一步做得不好，后面用什么算法都白搭。

1.2 初始中心敏感：一个被低估的稳定性问题

Kmeans的本质是坐标下降法，目标函数是非凸的，容易收敛到局部最优。最典型的症状就是：同一个数据集，你跑10次Kmeans，可能得到3种甚至更多种不同的分簇方案。

很多人觉得这个“没问题，多跑几次选SSE最小的呗”。但实际工作里问题很大。第一，如果聚类结果随机波动，意味着你无法稳定地给每个用户打标签，今天跑出来的“夜猫子型”用户，明天可能变成“全天均衡型”，后续的营销策略、负荷预测模型全都跟着飘。第二，单纯比较SSE选最优，并不能保证选到的是业务上有意义的分簇。Kmeans天然倾向于把大簇切碎，有时候SSE小了，用户画像反而乱了。

k-means++初始化算法能缓解一部分问题，但只是“缓解”，不能根治。居民负荷数据特征维度之间往往还存在相关性，比如峰段占比高的人谷段占比就低，这类共线性特征会让Kmeans的局部收敛问题更严重。

1.3 为什么不用层次聚类、DBSCAN或高斯混合模型

先说层次聚类。它对距离矩阵的依赖很强，计算复杂度高，几十万个居民用户根本算不动。就算只取5000户做分析，凝聚层次聚类的计算量也不小。而且层次聚类一旦合并错误就回不了头，对居民用电这种噪声大户不够友好。

DBSCAN适合发现不规则形状的簇，但它依赖密度阈值参数，居民负荷特征空间的密度差异非常大，一部分用户在特征空间里非常密集，另一部分则零散分布，很难找到全局适用的邻域参数。高斯混合模型可以看作软聚类版Kmeans，对重叠簇的处理更好，但它假设数据服从高斯混合分布，而且要估计协方差矩阵，特征维度稍微高一点，参数数量就爆炸。

所以Kmeans结合粒子群优化，是在工程可行性和聚类效果之间比较平衡的选择。Kmeans本身速度快、理解门槛低，粒子群负责解决它对初始值敏感的问题，两边各干各擅长的活。

2. 粒子群算法：一只鸟和一群鸟的全局搜索博弈

2.1 粒子群更新机制的核心逻辑

粒子群算法（Particle Swarm Optimization，PSO）模仿鸟群觅食行为。每只鸟就是搜索空间里的一个粒子，代表一个候选解。它在空间里飞的时候，会参考两个经验：自己历史上找到过的最好位置，以及整个鸟群目前找到的最好位置。

数学上，每个粒子的速度和位置更新公式是：

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

x_i(t+1) = x_i(t) + v_i(t+1)

其中 w 是惯性权重，控制粒子沿原方向飞行的程度；c1、c2 是加速度常数，分别控制粒子向个体最优和全局最优靠拢的程度；r1、r2 是[0,1]之间的均匀随机数，保证搜索的随机性。

从直觉上理解，w 大，粒子探索新区域的能力强；c1 大，粒子倾向于回顾自己走过的路；c2 大，粒子容易被大家集中到当前最好的位置。三者必须平衡，如果 c2 过大而 w 过小，粒子群会快速收敛到某个局部区域，早期就早熟；如果 w 过大而 c1、c2 过小，粒子就在空间里乱飞，收敛很慢。

2.2 参数怎么设：一份可以直接用的经验值

针对Kmeans聚类中心优化这个具体场景，粒子群的参数设置不用太复杂。标准的建议是：

这些参数不是拍脑袋定的。w 从0.9递减到0.4是Shi和Eberhart在1998年提出的经典做法，后来几乎所有PSO变体都沿用这个思路。c1=c2=1.49445对应Clerc的收缩因子方法，保证算法收敛性比默认的2.0更稳。

2.3 为什么是PSO而不是遗传算法或模拟退火

遗传算法（GA）和模拟退火（SA）也能做全局优化，但在这个场景里PSO有几个实际优势。

第一，PSO没有交叉、变异那些算子，代码量小很多，Matlab里写核心循环不到30行。第二，PSO直接利用适应度函数值本身进行位置调整，不需要求梯度，也不需要做复杂的编码解码。聚类中心本身就是一组连续实数，直接压成一维向量就能当粒子，非常自然。第三，PSO对计算资源的消耗相对可控，每次迭代只需要算一遍所有粒子对应Kmeans划分的SSE，配合向量化代码很快。遗传算法通常需要维护种群、做选择交叉变异，调参维度更多，收敛速度也没优势。

说白了，不是GA不好，而是处理“给Kmeans找一组好初始中心”这个问题，PSO是性价比最高的工具箱。

3. 从原始负荷到特征向量：预处理决定聚类质量上限

3.1 数据清洗：先把“假用户”和“异常用户”摘出去

居民用电原始数据拿回来，不能直接进入聚类流程。我一般按下面几步做清洗。

第一步，删除无效用户。日用电量几乎常年为0的用户，比如长期空置房，直接单独分成“空置户”，不参与聚类。判断标准可以用：某个用户在统计周期内日用电量超过0的天数不足30%，就认为这个用户本质上不活跃。

第二步，筛掉极端大电量用户。居民用户里偶尔混着家庭作坊、违规商业用电，日电量可能是普通用户的几十倍。这些极端值对Kmeans的簇中心影响很大，会把几个簇的中心全部拉偏。我习惯用箱线图或者分位数法把日电量高于99.5%分位的用户单独拎出来，归为“大电量异常户”，留给稽查团队去核查。

第三步，处理缺失值和零值。智能电表偶尔离线，数据采集有缺失。对于短时间缺失（少于连续4个点），用前后线性插值补上；对于长时间缺失，直接判为数据质量不合格，不参与聚类。这里要注意，不要把所有零值都当缺失值。居民夜间用电可能真的是0，这时候硬插值会把夜间负荷抬得虚高，导致峰谷特征失真。

3.2 特征构造：不只取平均，要构造“行为形态”

原始负荷曲线经过清洗后，按用户聚合成特征向量。我常用的特征分四类。

第一类是总量水平特征：日均用电量、日最大用电量、日最小用电量，代表用户的基本盘子。

第二类是时间分布特征：峰段（8点到22点，以常见居民峰谷划分为例）用电占比、谷段（22点到次日8点）占比、早高峰占比、晚高峰占比。这些特征直接反映用户用电的时间偏好，是区分“上班族”和“居家型”的关键。

第三类是形态特征：负荷率（平均负荷除以最大负荷）、负荷变异系数、高峰时段的位置和宽度。形态特征描述用户用电行为的平滑程度和波动特征。

第四类是行为稳定性特征：工作日与休息日平均用电量的比值，或者一周内日用电量的标准差。能区分“作息规律型”和“随机波动型”。

我一般控制在10到15维特征，既能保留信息，又不会让特征空间过于稀疏。特征构造完，做一次零均值和单位方差的标准化。标准化这一步非常重要，如果不做，日均用电量几千瓦时这种量级会直接压过占比类特征，导致聚类结果基本就是按用电量大小排序，而不是按行为模式分群。

3.3 K值怎么选：肘部法则和轮廓系数的组合判断

Kmeans和PSO-Kmeans都需要预先指定K。我通常先跑纯Kmeans，用肘部法则看SSE随K变化的拐点，再用轮廓系数交叉验证。

肘部法则不是自动的，要人工判断哪里“肘部明显”。居民用电数据叠了很多噪声，有时候拐点不明显。这种情况下我会画两条曲线：一条是SSE，一条是轮廓系数平均值。轮廓系数兼顾类内紧密度和类间分离度，取值在-1到1之间，越大越好。一般选择轮廓系数开始下降或者趋于平缓之前的K值。

业务侧的经验也很重要。做过几个项目之后，我发现居民用电聚类K取4到7比较合适。K太小，一个簇里混了好几种行为模式，业务上没法用；K太大，分出来的簇太碎，也没意义。K=5通常能兼顾解释性和区分度。

4. PSO-Kmeans融合模型设计与Matlab代码实现

4.1 粒子编码方式：把聚类中心压成一维向量

假设数据有 N 个用户，每个用户有 D 个标准化特征（我前面建议10到15个），总共要聚成 K 类。每个粒子必须代表“一组完整的聚类中心”。

编码方式很直接：把K个中心按顺序拼接成一个一维向量，向量长度等于 K 乘以 D。

比如D=10，K=5，粒子的位置向量长度就是50。第1到10维是第一个簇的中心，第11到20维是第二个簇的中心，以此类推。在Matlab的适应度函数里，用reshape函数把粒子向量还原成 K 行 D 列的矩阵，就能算距离了。

这种编码方式的优点是完全贴合Kmeans的计算逻辑，不需要额外设计解码规则。粒子群的搜索空间就是聚类中心的连续实数空间，每一维的边界就用数据在对应特征上的最小值和最大值圈起来，粒子不至于飞到不可能出现的特征值范围外面去。

4.2 适应度函数：紧凑性加空簇惩罚

粒子群优化Kmeans，目标是让每个样本到它最近簇中心的距离平方和最小。这个目标就是Kmeans聚类算法自己在迭代时最小化的SSE（误差平方和）。

我会在SSE基础上加一个空簇惩罚项。粒子群搜索过程中，会频繁出现几个粒子跑到同一个区域形成重复中心，导致某些簇没有样本。如果不对空簇做惩罚，粒子群会发现“抛弃几个簇，把中心集中在稠密区域”能让SSE下降，从而收敛到一个所有用户全被塞进一个簇的极端结果。

加了空簇惩罚后，每次算完最近中心归属，统计实际有样本的簇个数，如果少于K，就把SSE加上一个很大的惩罚值，比如1e6。这样粒子群会自动避开空簇区域。我在代码里也是这么实现的。

4.3 用Matlab自带的particleswarm还是自己写PSO

Matlab的Global Optimization Toolbox自带particleswarm函数，可以直接用，也可以自己写一个标准PSO循环。两种方式我都试过，各有利弊。

自带particleswarm的好处是稳定、有完备的迭代退出条件和并行计算支持，代码量小。缺点是它是一个通用优化器，没有专门为聚类场景优化粒子的初始化策略，需要你自己把适应度函数封装好。另外，particleswarm对函数句柄的调用模式有要求，适应度函数里做聚类计算时可以传入额外参数。

自己写PSO的好处是灵活性高，可以在初始化、速度更新、边界处理里注入针对Kmeans的特殊策略，比如用Kmeans++生成初始粒子的一部分。缺点是需要自己处理收敛判断、速度边界、粒子越界映射这些问题，代码稍长。

我的建议是：如果只是想把流程快速跑通，直接用particleswarm；如果需要在算法层面做改进、发表论文或者做详细实验对比，自己写PSO更顺手。下面两种方式的核心适应度函数代码是一致的。

4.4 核心Matlab代码框架

先给适应度函数写法的参考。假设X是标准化后的特征矩阵，每一行是一个用户，每一列是一个特征：

function sse = psoClusterObjFun(x, X, K) % x: 粒子位置，行向量，长度为 K * D % X: 标准化后的特征矩阵，N * D % K: 聚类簇数 D = size(X, 2); Cent = reshape(x, K, D); % 还原成 K * D 的聚类中心矩阵 Dist = pdist2(Cent, X, 'squaredeuclidean'); % K * N 距离矩阵 [minDist, labels] = min(Dist, [], 1); sse = sum(minDist); % 空簇惩罚。unique(labels)得到实际有样本的簇编号 nFilled = numel(unique(labels)); if nFilled < K sse = sse + 1e6; end end

pdist2的'squaredeuclidean'选项直接算欧氏距离的平方，这样 minDist 求和就是我们要的SSE。

如果直接用particleswarm：

nvars = K * D; lb = repmat(min(X), K, 1); % 每维特征的下界 ub = repmat(max(X), K, 1); % 每维特征的上界 lb = lb(:)'; ub = ub(:)'; options = optimoptions('particleswarm', ... 'SwarmSize', 40, ... 'MaxIterations', 200, ... 'HybridFcn', [], ... 'Display', 'iter'); fun = @(x) psoClusterObjFun(x, X, K); [bestX, bestSSE] = particleswarm(fun, nvars, lb, ub, options); % 用PSO得到的最优中心跑一次Kmeans精调 finalCent = reshape(bestX, K, D); [labels, finalCent] = kmeans(X, K, 'Start', finalCent, 'MaxIter', 1000);

注意kmeans的'Start'参数直接传入聚类中心初值。PSO结果已经处在全局较优区域，再交给标准Kmeans迭代几步收敛，能得到更好的SSE。这一步是我做实验时的标准配置，效果比PSO直接输出好一些。

如果想自己写标准PSO并嵌入Kmeans，核心循环是这样的：

% 初始化 nParticles = 40; dim = K * D; pos = repmat(lb, nParticles, 1) + rand(nParticles, dim) .* repmat(ub - lb, nParticles, 1); vel = zeros(nParticles, dim); pbestPos = pos; pbestVal = inf(nParticles, 1); gbestVal = inf; wMax = 0.9; wMin = 0.4; c1 = 1.49445; c2 = 1.49445; maxIter = 200; for iter = 1:maxIter w = wMax - (wMax - wMin) * iter / maxIter; for i = 1:nParticles val = psoClusterObjFun(pos(i, :), X, K); if val < pbestVal(i) pbestVal(i) = val; pbestPos(i, :) = pos(i, :); end if val < gbestVal gbestVal = val; gbestPos = pos(i, :); end end % 速度更新和位置更新 r1 = rand(nParticles, dim); r2 = rand(nParticles, dim); vel = w * vel + c1 * r1 .* (pbestPos - pos) + c2 * r2 .* (gbestPos - pos); pos = pos + vel; % 边界处理：越界的粒子拉回边界，速度也限制在Vmax范围内 Vmax = 0.1 * (ub - lb); vel = max(min(vel, Vmax), -Vmax); pos = max(min(pos, ub), lb); end

这个自写版本只保留了标准PSO的核心机制，没有加复杂的拓扑结构，足够处理居民用电聚类问题。实测下来，对几千户级别的数据，迭代200次很快，通常几十秒内完成。

4.5 融合策略的两种变体及实验结果对比

PSO和Kmeans融合并不只有一种方式。我常用的是“PSO提供初值+标准Kmeans精调”，还有一种是在PSO迭代里嵌Kmeans局部搜索：每迭代几次，把当前全局最优粒子的中心矩阵传给kmeans跑一步，再把精调后的结果反写回粒子位置。后者的收敛更稳，但计算代价更大。

我在一个3000户、10维特征、K=5的样本上做过快速对比。单纯Kmeans随机初始化跑20次，SSE的平均值波动明显；PSO-Kmeans跑20次，最终SSE基本稳定，簇中心也基本一致。轮廓系数方面，PSO-Kmeans的均值比纯Kmeans最优值略有提升，关键是多次运行的标准差大幅下降，这意味着结果可解释、可复现，在业务上非常有用。

这些数值具体是多少取决于数据分布，但趋势是稳定的。PSO-Kmeans在稳定性上的优势，对我来说比SSE绝对值下降更有价值。

5. 聚类结果怎么解读：从簇中心到用户行为画像

5.1 簇中心特征向量还原成负荷曲线

聚类跑完之后，不要直接对着标准化特征向量做业务判断。标准化后的特征值没有量纲，你看到某个簇的第一个特征值是0.63，根本不知道反映到实际用电量上是多少千瓦时。

正确做法是保存特征构造时的均值、标准差或者最大最小值，把簇中心还原到原始特征空间，再画负荷曲线。比如簇中心在第2维“晚高峰占比”上的标准化值是0.8，还原后对应晚高峰用电占全天用电的45%，业务人员一看就能联想到“晚上回家开空调开热水器的上班族”。

5.2 典型的居民用电行为画像

K=5的情况下，我常看到这样几类典型画像，簇中心曲线的形态差异非常清楚。

第一类：白天平稳型。日负荷曲线全天波动小，负荷率很高，特征是白天大部分时间有人在家用电，大概率是退休老人或者居家办公人群。他们对电价不敏感，但空调负荷占比较高，夏季是台区尖峰的主要推手之一。

第二类：早晚双峰型。特征上出现明显的早高峰和晚高峰，白天负荷低，典型上班族。这类用户对分时电价有一定响应能力，晚高峰的空调、热水器负荷是可削减的潜力来源。

第三类：深夜活跃型。特征是夜间电量占比极高，白天几乎平线，大概率是拥有电动车充电桩的用户，或者夜间从事生产活动的小作坊。深夜用电对电网来说是填谷资源，是在制定低谷电价激励时最值得争取的一类用户。

第四类：随机大功率型。日用电量波动极大，个别天出现很高的尖峰，负荷曲线形态每一天都不一样。这类用户可能拥有电采暖设备、即热式热水器等，对网格变压器冲击大，是台区改造和增容需求的关注对象。

第五类：空置低量型。用电量极低但又不是常年为零，可能是老人偶尔居住或者民宿短租。这类用户不需要营销干预，但要及时识别，避免在台区线损计算里造成干扰。

5.3 让非技术背景的决策者理解聚类结果

算法工程师容易陷入“聚类结果看起来挺合理”的自我满足里，但业务部门要的是决策依据。我在汇报时习惯做三张图。

第一张是聚类中心的典型负荷曲线对比图，把五类用户的日均曲线画在一个坐标系里，用不同颜色区分。第二张是特征雷达图，选几个核心特征比如峰段占比、谷段占比、负荷率，把五类用户的特征值标准化后画在雷达图上。第三张是用户分布地图或柱状图，展示每类用户的数量和电量占比。

这三张图讲完之后，业务部门能清楚地知道“哪类用户贡献了尖峰负荷”“哪些用户适合参与削峰响应”“哪些台区存在空置户干扰”。聚类算法本身不是目的，聚类结果能推动决策落地才是目的。

6. 实战踩坑记录：粒子群优化Kmeans并不总是一帆风顺

6.1 空簇惩罚权重设得太小，粒子群会“躺平”

我最早调试PSO-Kmeans时，空簇惩罚只加了SSE的10%，想着稍微给点压力就够了。结果迭代到后半段，粒子群发现干脆丢掉一个簇，把五个中心都挤在密集区域，SSE反而更小。于是最终结果只有四个簇有样本，另一个簇中心落在数据边缘的奇怪位置。

后来我把惩罚项直接改成固定大常数，空簇直接判死刑。实际效果证明，Kmeans场景下空簇几乎永远不可能是最优解，所以惩罚宁可大不要小。当然，如果你故意想找K个簇以外的特殊模式，那是另一回事。

6.2 归一化之后，簇中心要“翻译”回原始空间

有一段时间我直接拿标准化特征聚类，然后画簇中心曲线，发现图像完全看不懂。标准化后日均用电量为0.8之类的数值，根本没法跟业务人员解释“0.8是什么意思”。

解决方案是保存标准化参数，在可视化前把簇中心反向还原。还有一个细节：特征标准化时用到的均值和标准差，必须来自训练集本身。如果后面要接入新用户做标签预测，需要用同样的均值和标准差做标准化，不能重新计算，否则新用户的特征分布跟训练聚类时的分布不一致，分类结果会系统性偏掉。

6.3 粒子群早熟收敛：问题往往出在速度和惯性权重

PSO-Kmeans收敛得太快，并不总是好事。我遇到过跑20代就停止更新，所有粒子都堆在一个局部区域的情况，最终SSE还不如多跑几次随机初始化的标准Kmeans。

检查下来主要是两个原因。一个是Vmax设太小，粒子飞两步就被钳制住，失去探索能力；另一个是w衰减速度太快，从0.9降到0.4的曲线太陡，后期粒子几乎没有“惯性”，完全被全局最优牵着走，陷入局部最优出不来。调整方式是Vmax设置为每维范围的0.2倍，w衰减改为按迭代次数的线性衰减并延长迭代次数到200代以上。

6.4 特征维度和样本量的性能平衡

Matlab的pdist2在样本量达到几十万、特征维度几十维的时候，内存消耗会明显上升，特别是算K*N的距离矩阵。居民用户量级如果是百万级，不建议一次性全部聚类，按台区或者按网格分批聚类更合理。每批几万户，聚类速度很快，结果在业务上也更贴合实际，因为不同台区的居民用电行为受区域气候和经济水平影响很大，强行全城聚成一套画像反而没有指导意义。

另外，Matlab的parpool并行池可以加速particleswarm的粒子评估。不过要注意，适应度函数里如果用了匿名函数捕获大矩阵X，并行worker传递数据会有额外开销，数据量不大的时候反而不如串行快。实践下来，样本量小于5万，串行更省时间。

6.5 多次运行取最优的“土办法”依然有效

即便用了PSO-Kmeans，由于PSO本身也依赖随机初始粒子，每次运行结果还是会有一点点波动。正规做法是设置随机种子固定复现，工程上我更习惯另一种思路：跑10次，每次随机种子不同，取SSE最低的那次作为最终模型。

这不是耍滑头，而是利用高维非凸问题的特点——多次采样取最优本质上是一种朴素但有效的全局优化补充。跟“PSO+精调Kmeans”配合起来，稳定性已经足够满足业务复现要求。如果连这10次结果的SSE都还有明显差异，那就要回头检查数据预处理和特征构造了，大概率是数据里混进了不该有的异常结构。

写在最后：一点个人实操体会

把粒子群算法和Kmeans结合，技术细节不算复杂，真正的难点在于全程保持对数据的敏感。特征做得好，PSO参数随便设都能出来合理结果；特征做得糙，算法再高级也只是把噪声聚出个看似合理的形状。我写这篇文章的初衷就是把这个过程完整复述一遍，希望能帮到正在做居民用电聚类或者类似项目的朋友。如果你在自己的数据上试过这套流程，欢迎交流你们那边的聚类结果和踩坑经验，这种问题在实际项目里确实是常碰常新。