去年我在做居民用电行为分析时用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 OptimizationPSO模仿鸟群觅食行为。每只鸟就是搜索空间里的一个粒子代表一个候选解。它在空间里飞的时候会参考两个经验自己历史上找到过的最好位置以及整个鸟群目前找到的最好位置。数学上每个粒子的速度和位置更新公式是v_i(t1) w * v_i(t) c1 * r1 * (pbest_i - x_i(t)) c2 * r2 * (gbest - x_i(t))x_i(t1) x_i(t) v_i(t1)其中 w 是惯性权重控制粒子沿原方向飞行的程度c1、c2 是加速度常数分别控制粒子向个体最优和全局最优靠拢的程度r1、r2 是[0,1]之间的均匀随机数保证搜索的随机性。从直觉上理解w 大粒子探索新区域的能力强c1 大粒子倾向于回顾自己走过的路c2 大粒子容易被大家集中到当前最好的位置。三者必须平衡如果 c2 过大而 w 过小粒子群会快速收敛到某个局部区域早期就早熟如果 w 过大而 c1、c2 过小粒子就在空间里乱飞收敛很慢。2.2 参数怎么设一份可以直接用的经验值针对Kmeans聚类中心优化这个具体场景粒子群的参数设置不用太复杂。标准的建议是参数经验取值说明粒子数30~50聚类中心数量小的时候30个够用特征维度高可以加到50惯性权重 w0.9 线性递减到 0.4前期多探索后期多收敛加速度常数 c1. c21.49445经典收缩因子配置允许在参照系内使用速度上限 Vmax每维搜索范围的10%~20%防止粒子飞得太远导致适应度计算失去意义迭代次数150~300居民负荷特征空间比较简单通常200代以内收敛这些参数不是拍脑袋定的。w 从0.9递减到0.4是Shi和Eberhart在1998年提出的经典做法后来几乎所有PSO变体都沿用这个思路。c1c21.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太大分出来的簇太碎也没意义。K5通常能兼顾解释性和区分度。4. PSO-Kmeans融合模型设计与Matlab代码实现4.1 粒子编码方式把聚类中心压成一维向量假设数据有 N 个用户每个用户有 D 个标准化特征我前面建议10到15个总共要聚成 K 类。每个粒子必须代表“一组完整的聚类中心”。编码方式很直接把K个中心按顺序拼接成一个一维向量向量长度等于 K 乘以 D。比如D10K5粒子的位置向量长度就是50。第1到10维是第一个簇的中心第11到20维是第二个簇的中心以此类推。在Matlab的适应度函数里用reshape函数把粒子向量还原成 K 行 D 列的矩阵就能算距离了。这种编码方式的优点是完全贴合Kmeans的计算逻辑不需要额外设计解码规则。粒子群的搜索空间就是聚类中心的连续实数空间每一维的边界就用数据在对应特征上的最小值和最大值圈起来粒子不至于飞到不可能出现的特征值范围外面去。4.2 适应度函数紧凑性加空簇惩罚粒子群优化Kmeans目标是让每个样本到它最近簇中心的距离平方和最小。这个目标就是Kmeans聚类算法自己在迭代时最小化的SSE误差平方和。我会在SSE基础上加一个空簇惩罚项。粒子群搜索过程中会频繁出现几个粒子跑到同一个区域形成重复中心导致某些簇没有样本。如果不对空簇做惩罚粒子群会发现“抛弃几个簇把中心集中在稠密区域”能让SSE下降从而收敛到一个所有用户全被塞进一个簇的极端结果。加了空簇惩罚后每次算完最近中心归属统计实际有样本的簇个数如果少于K就把SSE加上一个很大的惩罚值比如1e6。这样粒子群会自动避开空簇区域。我在代码里也是这么实现的。4.3 用Matlab自带的particleswarm还是自己写PSOMatlab的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 endpdist2的squaredeuclidean选项直接算欧氏距离的平方这样 minDist 求和就是我们要的SSE。如果直接用particleswarmnvars 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维特征、K5的样本上做过快速对比。单纯Kmeans随机初始化跑20次SSE的平均值波动明显PSO-Kmeans跑20次最终SSE基本稳定簇中心也基本一致。轮廓系数方面PSO-Kmeans的均值比纯Kmeans最优值略有提升关键是多次运行的标准差大幅下降这意味着结果可解释、可复现在业务上非常有用。评价指标纯Kmeans多次运行PSO-Kmeans多次运行SSE均值偏高更低SSE标准差明显很小轮廓系数均值中等更高轮廓系数标准差较大很小这些数值具体是多少取决于数据分布但趋势是稳定的。PSO-Kmeans在稳定性上的优势对我来说比SSE绝对值下降更有价值。5. 聚类结果怎么解读从簇中心到用户行为画像5.1 簇中心特征向量还原成负荷曲线聚类跑完之后不要直接对着标准化特征向量做业务判断。标准化后的特征值没有量纲你看到某个簇的第一个特征值是0.63根本不知道反映到实际用电量上是多少千瓦时。正确做法是保存特征构造时的均值、标准差或者最大最小值把簇中心还原到原始特征空间再画负荷曲线。比如簇中心在第2维“晚高峰占比”上的标准化值是0.8还原后对应晚高峰用电占全天用电的45%业务人员一看就能联想到“晚上回家开空调开热水器的上班族”。5.2 典型的居民用电行为画像K5的情况下我常看到这样几类典型画像簇中心曲线的形态差异非常清楚。第一类白天平稳型。日负荷曲线全天波动小负荷率很高特征是白天大部分时间有人在家用电大概率是退休老人或者居家办公人群。他们对电价不敏感但空调负荷占比较高夏季是台区尖峰的主要推手之一。第二类早晚双峰型。特征上出现明显的早高峰和晚高峰白天负荷低典型上班族。这类用户对分时电价有一定响应能力晚高峰的空调、热水器负荷是可削减的潜力来源。第三类深夜活跃型。特征是夜间电量占比极高白天几乎平线大概率是拥有电动车充电桩的用户或者夜间从事生产活动的小作坊。深夜用电对电网来说是填谷资源是在制定低谷电价激励时最值得争取的一类用户。第四类随机大功率型。日用电量波动极大个别天出现很高的尖峰负荷曲线形态每一天都不一样。这类用户可能拥有电采暖设备、即热式热水器等对网格变压器冲击大是台区改造和增容需求的关注对象。第五类空置低量型。用电量极低但又不是常年为零可能是老人偶尔居住或者民宿短租。这类用户不需要营销干预但要及时识别避免在台区线损计算里造成干扰。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参数随便设都能出来合理结果特征做得糙算法再高级也只是把噪声聚出个看似合理的形状。我写这篇文章的初衷就是把这个过程完整复述一遍希望能帮到正在做居民用电聚类或者类似项目的朋友。如果你在自己的数据上试过这套流程欢迎交流你们那边的聚类结果和踩坑经验这种问题在实际项目里确实是常碰常新。