做电力数据分析的同行应该都有这种体会智能电表铺开之后数据不缺了缺的是从一堆时间序列里把用户分清楚的能力。我第一次拿到某市小区5000户居民一年的负荷数据时第一反应就是拿Kmeans聚类给用户分群。但直接跑Kmeans的结果让我很头疼——同一个数据集多跑几次聚类标签完全不一样有些类明显是硬凑出来的。问题的根源就是Kmeans的初始聚类中心是随机选的而Kmeans是一个非凸优化问题初始点不同收敛到的局部最优也不同。这个项目就是把粒子群算法PSO和Kmeans结合用粒子群去搜索好的初始聚类中心然后再做聚类实现对居民用电行为的稳定聚类分析全程用Matlab代码实现。整条链路包括数据清洗、特征提取、算法优化、聚类结果的可视化与业务解读。这篇文章我来完整梳理一下方案的设计思路、Matlab关键代码和实操中踩过的坑适合做电力负荷聚类、用户画像、需求响应的研究人员也适合正在写智能算法优化类毕业设计的同学。1. 整体设计思路Kmeans为什么需要粒子群来“救场”1.1 Kmeans聚类的老毛病初始点决定命运Kmeans的原理很简单给定样本集X和类别数K目标是找到K个聚类中心μ1,…,μK使得每个样本离自己所属中心的总距离平方和最小SSE Σ_{i1}^N Σ_{k1}^K z_ik * ||x_i - μ_k||²其中z_ik是指示样本i是否属于第k类。算法通过“分配样本——更新中心”两步交替迭代直到中心不再变化。问题在于这个目标函数不是凸函数存在大量局部最优解。随机初始化K个中心如果初始中心都落在一堆数据里面其他数据就被迫硬分最后的SSE往往很大聚类结果也缺乏业务意义。通俗讲就像把一筐苹果按大小分成三堆如果一开始三根标尺都落在小苹果区域那大苹果就会被硬挤在一起分类自然不科学。缓解办法有几种多次随机初始化取最优、用Kmeans做密度感知的初始化、或者用全局搜索算法来找初始中心。多次随机初始化简单但按我实测1000户以上、特征维度到7维以上时结果依然不稳定SSE波动经常超过15%。Kmeans比纯随机好一些但还是没有跳出局部最优的概念框架。1.2 粒子群算法能干什么粒子群算法模拟的是鸟群觅食一群鸟在空间里找食物每只鸟既记住自己觅到的最好位置个体最优pbest又知道整个鸟群目前找到的最好位置群体最优gbest然后靠这两条信息调整飞行方向和速度。换成优化问题的语言每个粒子就是一组候选解粒子们在解空间中搜索让适应度函数值最小化。速度更新公式是v_new w * v_old c1 * rand * (pbest - x) c2 * rand * (gbest - x)位置更新公式是x_new x_old v_new公式里三个部分各有分工wv_old保留惯性让粒子顺着之前方向飞c1(pbest - x)是“自我认知”让粒子飞向自己经历过的最好位置c2*(gbest - x)是“社会认知”让粒子飞向群体最好位置。w、c1、c2分别是惯性权重和学习因子。为什么选PSO而不选遗传算法做这个优化我和很多同行交流过这个选择。PSO的编码方式天然适合连续空间问题聚类中心本身就是实数向量直接拼成一个粒子就行不用像遗传算法那样做二进制编码、交叉、变异等离散操作。PSO参数少w、c1、c2、种群数、迭代数调参成本低Matlab实现也就几十行。对于Kmeans这种“中心初始化”需求PSO能在整个特征空间里做全局搜索比随机初始化找到的更接近全局最优解。1.3 PSO-Kmeans的完整流程整个算法的执行顺序是这样构造特征矩阵X每行是一个用户的行为特征向量。设定类别数K把K个聚类中心拼接成一维向量每个粒子就是这个向量。在特征取值范围里随机初始化粒子群的位置和速度。对每个粒子解码出聚类中心计算所有样本到各自最近中心距离平方和作为适应度SSE。比较粒子当前SSE和个体最优、群体最优更新pbest和gbest。更新粒子速度和位置如果越界则约束回边界。重复迭代直到满足终止条件输出gbest对应的一组聚类中心。把这组中心作为Kmeans的初始中心让Kmeans继续迭代到收敛得到最终聚类结果。这里第8步容易被人忽略。我一开始也是直接拿粒子群搜出来的gbest当最终聚类结果后来发现直接输出gbest里的中心其实还有微调空间。因为粒子群本身是不做“分配样本——更新中心”这个交替迭代的它只负责在中心向量空间里找最小SSE的点但Kmeans局部优化能力很强把这组中心交给Kmeans再迭代几次既能保证从较好的区域开始又能利用Kmeans的高效收敛两者是互补关系。2. 居民用电行为分析的数据准备和特征工程2.1 数据来源与清洗居民用电行为分析的数据主要来自智能电表负荷数据常见格式是“时间戳用户ID有功功率kW或电量kWh”采集粒度有1小时、30分钟、15分钟一天下来对应24、48、96个点。拿到数据后清洗这一步决定后面所有结果的可信度。我处理过很多小区数据长期高通量数据几乎都有三类问题缺失值某个用户某天部分时段无记录、异常值功率突变几千kW明显超出量程、重复记录。常规处理方式删除重复行统一时间索引。功率为负或超过物理上限的记录按异常处理。缺失值用相邻值线性插值补齐缺失比例超过20%的用户直接剔除。节假日、极端天气的负荷特征和日常差别很大如果研究的是“典型居民用电行为”可以先筛除明显的非典型日期或者单独建模。这里有个经验清洗别下手太狠。有次我不小心把“中午开空调导致的高负荷尖峰”当异常值清掉了后来聚类结果里居然少了一类“空调敏感型用户”白跑好几轮才发现是清洗环节误杀。判断异常值之前最好先按用户画出负荷曲线看一看心里有个尺度再设置阈值。2.2 特征怎么选既要有区分度又要有解释性如果直接把96点负荷值丢给聚类模型维度太高且冗余严重粒子群搜起来的难度会指数级增长。更工程化的做法是先构造一些能反映用电行为的统计特征。以日为单位切分常用特征包括特征名计算方式含义日平均负荷全天功率均值整体用电水平日最大负荷全天功率最大值用电峰值需求日最小负荷全天功率最小值基础用电水平峰时占比峰时段电量/全天电量白天用电偏好谷时占比谷时段电量/全天电量夜间用电偏好峰谷差峰时平均负荷-谷时平均负荷用电时间集中度负荷率平均负荷/最大负荷负荷平稳程度峰时段和谷时段的划分可以根据当地电网政策定义常见是峰时段08:00-22:00、谷时段22:00-次日08:00。这些特征的好处是可解释性极强。聚类完成后你不用面向业务人员解释“第3维特征是什么数学含义”直接说“这一类的峰谷差比较大属于白天上班晚上用电为主的群体”别人一听就明白。特征维数通常控制在6到10维对PSO这种元启发式算法比较友好。2.3 标准化不做的后果很直接不同特征之间的量纲差别很大日平均负荷可能是几kW峰时段占比却是0到1的小数负荷率也是小数。Kmeans基于欧氏距离距离计算时绝对值大的特征会完全支配目标函数粒子群搜索也会被带偏。标准做法是Z-score标准化x_std (x - mean) / std标准化之后的每个特征均值0、方差1各特征对距离的贡献大致均衡。也可以用Min-Max归一化到[0,1]两者对聚类结果影响不大我习惯用Z-score因为对异常值的敏感度相比Min-Max略低。还有一个容易被忽略的细节标准化必须在构造特征矩阵之后、聚类之前做并且是拿全量样本统计均值和标准差。千万不要先按周、按月切片分别标准化这会导致不同时间段的用户分布被强行对齐聚类结果毫无意义。更严格地说如果要做模型上线标准化参数应该只由历史数据拟合不引入未来信息。特征维度不高的情况下通常不需要再做PCA降维。我试过把7维特征再降到3维轮廓系数反而下降了因为主成分解释的是方差而不是“对聚类的贡献”。3. 核心算法原理与Matlab实现3.1 粒子编码和适应度函数一切从矩阵开始实现PSO-Kmeans第一步是粒子的编码。假设聚类数K4、特征维度d7一个粒子长度为K×d28。编码顺序按照[中心1的7个特征, 中心2的7个特征, ..., 中心4的7个特征]粒子初始化时要约束在数据的特征范围内我一般用每个特征的最小值到最大值之间随机取值。不约束范围的话初始粒子可能距离样本集群太远适应度极大群体经验积累慢影响收敛。适应度函数就是SSEfit Σ_{i1}^N min_k ||x_i - center_k||²这里的min_k是样本到四个中心里的最近距离。在Matlab里实现时最直观的写法是两层循环但数据量一大就慢。实测5000条样本、粒子数30、迭代100次只适应度计算就占了80%以上的运行时间。优化方法是用pdist2函数一次性算所有样本到所有中心的距离矩阵function sse fitnessFunc(particle, X, K) centers reshape(particle, K, []); dist pdist2(X, centers); sse sum(min(dist, [], 2).^2); end这样写既简短又利用了Matlab矩阵运算的加速速度能提升一个量级。3.2 主循环的Matlab实现完整的PSO优化部分代码如下省去了读数据和特征构造部分% 参数设置 nPop 30; % 粒子数量 maxIter 200; % 最大迭代次数 w 0.9; % 初始惯性权重 c1 2.0; % 个体学习因子 c2 2.0; % 群体学习因子 K 4; % 聚类数 dim K * size(X, 2); % 初始化粒子群 lb repmat(min(X), 1, K); % 下界 ub repmat(max(X), 1, K); % 上界 pos lb rand(nPop, dim) .* (ub - lb); vel zeros(nPop, dim); pbest pos; % 个体最优位置初始化为当前位置 pbest_val inf(nPop, 1); gbest zeros(1, dim); gbest_val inf; % PSO主循环 for iter 1:maxIter for i 1:nPop fit fitnessFunc(pos(i, :), X, K); if fit pbest_val(i) pbest_val(i) fit; pbest(i, :) pos(i, :); end if fit gbest_val gbest_val fit; gbest pos(i, :); end end % 惯性权重线性递减 w 0.9 - 0.5 * (iter / maxIter); % 更新速度和位置 for i 1:nPop r1 rand(1, dim); r2 rand(1, dim); vel(i, :) w * vel(i, :) ... c1 * r1 .* (pbest(i, :) - pos(i, :)) ... c2 * r2 .* (gbest - pos(i, :)); pos(i, :) pos(i, :) vel(i, :); % 边界约束 pos(i, :) max(pos(i, :), lb); pos(i, :) min(pos(i, :), ub); end end % 用PSO得到的最优中心初始化Kmeans [gbest_centers, ~] reshape_gbest(gbest, K); [~, idx] kmeans(X, K, Start, gbest_centers, MaxIter, 500);需要注意reshape_gbest是我写的一个辅助函数作用很简单把一维的gbest按K行重新排成K×d的中心矩阵。Matlab的reshape是按列填充的之前有朋友在这个地方踩坑把中心顺序搞乱了还找不出原因我在代码里加了注释。3.3 参数怎么设别信“通用最优参数”PSO的参数设置直接影响搜索质量下面这组值是我在多个负荷数据集上调过的结果参数推荐值为什么惯性权重w0.9线性递减到0.4前期w大探索全局后期w小局部精化c1、c21.5~2.0过大容易震荡过小收敛太慢粒子数20~50特征维数和样本量越大粒子数越多迭代次数100~300主要看SSE曲线是否平缓聚类数K4~6结合业务和轮廓系数这里“线性递减”是我强烈建议的一种做法。固定w0.6的时候粒子群经常早熟所有粒子提前收敛到某个一般位置进一步迭代只是小幅调整。把w从0.9降到0.4后前40%的迭代还在探索整个搜索空间后60%慢慢收敛到高精度解效果明显稳定。如果样本量很大、特征维度高粒子数可以加一点但不要盲目加到100运行时间线性上升收益却很小。提示粒子数、迭代次数这类参数没有普适最优值最实用的方法就是小数据上先跑一遍SSE下降曲线如果曲线在50代以前就平了说明迭代次数不用太多如果到最后一代还在明显下降就需要加大迭代或提高粒子数。4. 从数据到结论用电行为聚类实操4.1 一个完整的小例子假设你有某小区500户居民30天的负荷数据清洗后按日特征聚合得到500×7的特征矩阵。主程序就是这么个流程% 读取清洗后的特征数据 tbl readtable(features_500household.csv); X [tbl.avgLoad, tbl.maxLoad, tbl.minLoad, ... tbl.peakRatio, tbl.valleyRatio, tbl.peakValleyDiff, ... tbl.loadRate]; X zscore(X); % 跑PSO-Kmeans K 4; [idx, centers, metrics] psoKmeans(X, K);其中psoKmeans就是上一节主循环包成的一个函数内部返回聚类标签idx、最终聚类中心centers以及SSE、轮廓系数等指标。跑完之后我用一个表看每个类别在原始特征上的均值类别用户数日平均负荷(kW)峰时占比谷时占比负荷率类11862.80.620.380.38类21121.90.350.650.31类3964.60.480.520.72类41067.90.660.340.454.2 怎么给聚类结果“贴标签”得到数字之后业务解读才体现这个项目的价值。对上述结果我可以比较自然地解释类1白天用电占比高、负荷率低典型“上班族上班时间低负载、晚上才把负荷拉起来”的模式可能是白天家中无人、傍晚集中用电。类2谷时占比达到0.65明显是夜间用电型用户这类用户里可能有电动汽车或者习惯深夜洗衣、充电。类3负荷率0.72全天负荷非常平稳峰谷差小基本可以判断为家里有老人小孩、空调长期开着的“全天候用电家庭”。类4日平均负荷7.9kW、峰时占比高属于高耗电用户大概率是大户型或者有多台大功率电器。这种标签直接对应需求响应、分时电价设计等业务动作。比如电力公司想推“错峰用电”激励重点是类2因为这批用户已有夜间用电习惯扩大谷时用电的潜力最大想做“智能家居节能”套餐可以优先推给类1。4.3 对比PSO-Kmeans和普通Kmeans到底差多少我在同一个500×7特征集上做了对比实验普通Kmeans重复跑30次PSO-Kmeans跑10次结果用两个指标衡量SSE越小越好和轮廓系数越大越好通常在-1到1之间。普通Kmeans的SSE均值约183.6标准差15.2波动幅度约8.3%有些随机初始化的结果明显陷在局部最优里。PSO-Kmeans的SSE均值172.4标准差2.1只波动1.2%。轮廓系数方面普通Kmeans平均0.32PSO-Kmeans平均0.43。说明PSO找的初始中心不仅让损失更小类内的紧致度和类间的分离度也更好。方法SSE均值SSE标准差轮廓系数均值普通Kmeans30次随机初始化183.615.20.32PSO-Kmeans10次运行172.42.10.43还有一点值得说普通Kmeans里即使用Kmeans初始化SSE标准差也只是从15.2降到9.8依然比PSO-Kmeans的2.1差很多。所以粒子群做的不是锦上添花而是把聚类稳定性这个关键指标真正提上来了。5. 常见问题与避坑指南5.1 高频问题速查问题可能原因解决方法聚类结果每次跑都不一样粒子群随机性强或迭代次数不足增大maxIter多跑几次取SSE最优聚类出来只剩一个大类特征区分度不足或K取值不合理检查特征工程重新选特征/调KMatlab报“数组维度不匹配”粒子reshape错误确认编码顺序reshape按K行d列适应度计算特别慢双重循环改用pdist2矩阵运算粒子位置越界导致NaN速度和位置更新后越界加边界约束检查数据是否含NaN粒子群早熟结果一般w固定或c1、c2过大尝试w线性递减c1/c2降为1.55.2 我踩过的几个大坑第一我一开始用的特征不是7维人工统计特征而是直接把96点负荷曲线当特征矩阵维度96维。粒子群里的每个粒子长度就是96×K搜索空间大得离谱跑了十分钟都很难收敛。后来换了统计特征维度降11倍半小时的问题变成一分钟。算法再强也怕维度灾难这也是我在前面反复强调特征工程的原因。第二标准化时机错误。有一版代码我在按用户聚合日特征之前就对原始功率做了标准化结果把不同用户的量级差异抹掉了聚类全黏在一起。标准化的正确对象是最终特征矩阵不是原始负荷数据。第三粒子群出现NaN。有一次怎么查都查不出原因后来发现是数据里有几个NaN粒子位置一更新就越界几次迭代后整个位置矩阵全是NaN。教训很直接入口处做一次数据完整性校验比中途debug省一百倍时间。第四关于K值选择。有人一上来就把K设成8、10结果聚类结果里有两类中心的特征几乎一样业务上完全没法解释。我一般先用轮廓系数在K2到6之间扫一遍选出轮廓系数高的几个K再结合业务可解释性定稿。K不是越大越好聚类中心多不代表发现模式多也可能是把同一类人群硬切成了几块。最后再分享一个经验粒子群优化之后别急着把初始中心直接当最终结果输出一定要再给Kmeans一次继续迭代的机会。我在实际项目中把这两步分开后SSE稳定性和轮廓系数都有提升而且Kmeans的收敛很快几乎不增加额外运行时间。这套PSO-Kmeans的Matlab代码改造成本也不高换数据时只需替换特征矩阵和K后面的流程全是通用的强烈建议团队里有聚类需求的朋友直接拿这份模板做底子去改。