最近在做居民用电行为分析时被手头一批智能电表小时级负荷数据折腾得不轻。事情本身不复杂——把几万个用户的用电曲线分类画像但直接用Matlab自带的kmeans跑结果实在没法用同样一份数据十次随机初始化跑出十种结果SSE波动幅度接近百分之十几。这个稳定性拿到科研报告里能挑好看的写但放在实际业务里根本说不过去。后来把粒子群算法PSO和Kmeans绑在一起让PSO先全局搜索出一组好的初始聚类中心再把这个中心作为Kmeans的起点做局部收敛。实验效果不错但把整个流程走通的过程中踩了不少坑。这篇文章就把我的完整思路、Matlab代码实现和排查教训整理出来给正在做居民用电聚类、PSO优化算法或者Kmeans相关课题的朋友一个可复用的参考。1. 为什么居民用电行为分析要动用到聚类算法智能电表普及之后电网公司手里掌握的用户负荷数据量级已经不是人工能处理的了。一个中型城市一天的采集数据就是几十万条曲线每条曲线包含24个甚至96个时点。面对这个规模的数据想要回答“哪些用户适合执行峰谷分时电价”“哪些用户具备可调负荷参与需求响应”这类问题第一步都是做用户分类。分类的方式有两种一种是根据用户台账手动打标签比如小区类型、户型面积、是否安装分布式光伏但这种标签往往滞后而且不完整另一种就是从负荷数据本身出发用无监督聚类把行为模式相似的用户聚到一组这就是用电行为分析的起点。1.1 用电行为聚类能解决什么实际问题居民用电行为聚类不是学术上的自嗨业务价值可以落到几个非常具体的方向。第一个是分时电价套餐设计。不同用户对电价的敏感度差异很大上班族白天空置家庭负荷低晚上回家负荷陡增退休老人家庭白天持续低负荷运行部分用户夏季空调用电占比极高。这些特征通过聚类分组之后就可以针对性地设计尖峰费率、低谷优惠、季节性套餐而不是对所有用户一刀切。第二个是需求响应潜力评估。现在的削峰填谷任务越来越重电网需要知道哪些用户群体在有激励时能主动压低负荷。聚类之后如果某个簇的负荷曲线有明显的用电集中时段、且这些时段恰好与电网峰时段重叠那么这个群体就是需求响应资源池中的优质对象。第三个是异常用电检测和客户细分服务。聚类本质上刻画的是“大多数用户”的用电规律那些偏离规律的样本自然成为重点关注对象。服务方面聚类画像可以作为精准营销的基础比如识别出有储能设备、有电动车充电习惯的用户群体针对性推荐增值服务。可以说聚类是连接原始负荷数据和上层业务决策的中间层。没有这个中间层机器学习、数据挖掘的手段再多也无法落到电力系统的实际管理动作上。1.2 为什么选Kmeans而不是其他聚类算法用电数据分析领域常见的聚类算法不少DBSCAN、层次聚类、谱聚类都有人用但我最终把Kmeans作为基础方案原因有三个。第一是数据规模。Kmeans的时间复杂度接近线性对几十万条样本依然可以在可接受时间内跑完。DBSCAN在实现较好的情况下效率也不错但在簇密度不均匀的数据上参数很难调而居民负荷曲线恰好是密度分布极不均匀——白天用电峰值和夜间低谷对应的样本分布差异非常大。第二是可解释性。Kmeans的聚类中心就是“平均用电曲线”每个簇可以用一条有物理意义的负荷曲线来描述。层次聚类也能做到但树状图在大样本下可读性极差谱聚类的嵌入空间不好向业务方解释。第三是工程实现成熟度。Matlab的kmeans函数集成了多种距离度量、并行计算和重复启动机制可以直接复用。但Matlab自带的Kmeans有一个老问题——初始中心随机选陷入局部最优是家常便饭。这为后面引入PSO优化提供了充分的动机。2. PSO如何给Kmeans找到更好的初始质心Kmeans聚类的核心目标是找到一组聚类中心使所有样本到所属中心距离的平方和最小。这个目标函数非凸自带Kmeans用的是随机初始化和多次重复启动本质上是在“碰运气”。粒子群算法在这里的角色就是替代随机初始化通过群体智能搜索在解空间中找到一个更接近全局最优的起始点。2.1 Kmeans的痛点初始质心对结果的影响远比你想象中大Kmeans迭代过程可以分成交替的两步第一步根据当前聚类中心把每个样本分配到最近的中心第二步根据分配结果重新计算每个簇的均值作为新中心。问题在于这个迭代过程是“局部爬山”。如果初始质心选得不好比如两个初始中心落在同一个数据密集区域那么最终聚类结果就会在局部最优处停下来无法跳到更合理的划分。真实数据里的表现就是随便跑几次kmeans每次的SSE都不一样轮廓系数也有明显起伏。我遇到的最夸张的一次同一份数据K4SSE最小值与最大值相差22%。这种情况下你没法判断当前这个结果是不是可用的。解决初始质心问题有几个常规手段。Kmeans的思路是让初始质心尽量分散这个策略对凸簇或者分布比较规则的簇有效但如果簇形状重叠、存在噪声Kmeans的初始中心仍然可能落在不合适的位置多次随机重启结合挑选最小SSE的做法本质上是用计算时间去换概率但重复次数再多也还是随机搜索没有目标导向。粒子群算法不一样它是带着“方向”去找解的——每个候选解都会计算适应度通过种群协作逐步逼近更优区域这正是Kmeans两点方案中最需要的稀缺能力。2.2 粒子群算法的主要思想与更新过程粒子群算法的灵感来自鸟群觅食。每只鸟就是一个粒子代表解空间中的一个候选解鸟群的飞行过程就是不断搜索最优解的过程。每个粒子具备两个属性位置和速度。如果把当前解空间记作D维空间第i个粒子的位置记为X_i速度记为V_i。算法为每个粒子保存两个记忆粒子自己的历史最优位置pbest以及整个种群目前发现的最优位置gbest。每次迭代时粒子的速度和位置按照下面的方式更新速度更新V_i^{t1} wV_i^t c1r1*(pbest_i - X_i^t) c2r2(gbest - X_i^t)位置更新X_i^{t1} X_i^t V_i^{t1}其中w是惯性权重决定当前速度的保持程度c1是认知学习因子c2是社会学习因子r1和r2是[0,1]之间的随机数。三个分量的物理意义很清楚惯性项让粒子沿当前趋势继续飞行认知项把粒子拉向自己曾经找到过的好位置社会项把整个种群引向全局最优位置。这种机制在连续域搜索问题中表现稳定参数少、实现简单、不需要梯度信息非常适合用来优化Kmeans的初始质心。2.3 粒子编码方式与适应度函数设计用PSO优化Kmeans初始质心最关键的环节是把Kmeans的解映射成PSO粒子。我采用的编码方式是如果聚成K个簇每个簇中心是d维特征向量那么一个粒子的位置就是一个长度为 K*d 的一维向量前d个元素是第1个簇中心的坐标接下来d个元素是第2个簇中心的坐标依次类推。适应度函数的设计直接影响优化效果。我用的适应度是样本到最近聚类中心的欧氏距离平方和即Kmeans的SSE目标函数。这个选择的理由很直接——PSO搜索到的适应度越小就意味着这组初始质心本身已经处在一个比较好的聚类划分附近。为什么不直接用轮廓系数作为适应度因为轮廓系数要计算所有样本两两之间的距离关系复杂度接近O(N²)几万条样本时每次适应度评估都慢得难以接受而SSE的计算是线性的在PSO上百次迭代中每次都算轮廓系数实验基本就等死了。还有一个常见做法值得提一下PSO每评估一次粒子时让这个粒子先跑几步Kmeans再计算适应度相当于是局部搜索和全局搜索的混合。这个方案效果好但时间复杂度翻了几倍我在实际项目中用的是朴素的两阶段法——PSO只负责找初始质心找到之后交给Kmeans做最终收敛。这样工程上简单效果已经足够替代随机初始化了。3. Matlab实现从数据清洗到结果评估整个实验我用的Matlab版本是R2024b代码实现层面没有用任何工具箱之外的特殊功能换到R2020之后的老版本也能跑。实现过程分为数据清洗、特征工程、PSO核心代码、参数配置和结果评估几个环节。3.1 数据准备与特征工程我手上的原始数据是一个居民小区的智能电表负荷记录每一行是一个用户每一小时一个采集点每天形成一条24维的负荷向量。为了研究用电行为需要扩展到一段时间的平均表现而不是单日偶然值。采集到的原始数据不能直接进聚类清洗规则是把零值时段占比超过80%的用户剔除掉这些表计可能已经拆表或者长期空置把峰值负荷超过电表额定容量的样本去掉基本是接线错误或者转供户连续多日缺数的用户也删除。清洗之后我取每户近三个月的有效数据计算特征。我最终用于聚类的特征一共5个维度工作日平均负荷曲线的峰时段用电占比工作日平均负荷曲线的谷时段用电占比休息日平均负荷水平归一化值负载率即平均负荷与最大负荷的比值峰谷差率反映日内用电波动幅度为什么不用完整24维曲线直接聚类从数据维度上考虑居民用电曲线24个点之间有强相关性直接聚类容易受到个别时点噪声的干扰也增加了PSO解空间的维度。比如K4时如果直接用24维数据粒子长度就是96用5个统计特征时粒子长度只有20搜索空间小了一个数量级收敛难度完全不同。数据处理上有一件必须做的事对所有特征做标准化。我用的是zscore函数也就是减去均值除以标准差。这个步骤不做的话量纲大的特征会在距离计算中占据绝对主导地位比如负荷值本身可能是几十个单位而负载率是0到1的小数聚类等价于只用负荷大小说话。% 读取清洗后的数据每一行是一个用户每一列是一个特征 data readmatrix(residential_features.csv); data_z zscore(data);3.2 PSO优化Kmeans的核心代码实现粒子群优化部分的完整代码并不复杂我拆成几个函数来写方便调试和复用。下面这段是PSO主流程框架。%% 参数设置 K 4; % 聚类数 N 30; % 粒子群规模 T 100; % 迭代次数 D size(data_z, 2); % 特征维度 w_max 0.9; % 惯性权重上限 w_min 0.4; % 惯性权重下限 c1 2; % 个体学习因子 c2 2; % 社会学习因子 %% 解空间边界聚类中心应落在数据取值范围内 lb repmat(min(data_z), 1, K); ub repmat(max(data_z), 1, K); %% 种群初始化 X lb rand(N, K * D) .* (ub - lb); % 粒子位置 V zeros(N, K * D); % 粒子速度 pbest X; pbest_fit inf(N, 1); gbest_fit inf; gbest X(1, :); %% 迭代 for iter 1:T w w_max - (w_max - w_min) * iter / T; for i 1:N centers reshape(X(i, :), K, D); [~, fit] computeSSE(data_z, centers); if fit pbest_fit(i) pbest_fit(i) fit; pbest(i, :) X(i, :); end if fit gbest_fit gbest_fit fit; gbest X(i, :); end end for i 1:N V(i, :) w * V(i, :) c1 * rand() * (pbest(i, :) - X(i, :))... c2 * rand() * (gbest - X(i, :)); X(i, :) X(i, :) V(i, :); % 边界还原粒子越界时拉回边界 X(i, :) max(min(X(i, :), ub), lb); end fprintf(iter %d, best fitness %.4f\n, iter, gbest_fit); end %% 用PSO找到的最优质心作为Kmeans的初始中心 init_centers reshape(gbest, K, D); [idx, C] kmeans(data_z, K, Start, init_centers, MaxIter, 500);适应度计算函数computeSSE的实现是关键之一。这里注意不要用双重循环逐个样本算距离要通过矩阵运算一次算出所有样本到某中心的距离function [labels, sse] computeSSE(data, centers) n size(data, 1); nk size(centers, 1); distMat zeros(n, nk); for k 1:nk diff data - centers(k, :); distMat(:, k) sum(diff.^2, 2); end [minDist, labels] min(distMat, [], 2); sse sum(minDist); end这个函数里遍历的是聚类中心数量k而不是样本数量n。k相对于n来说可以忽略不计所以整体效率是可以接受的。如果你的样本量特别大还可以在PSO内部改用pdist2来计算距离矩阵速度会更快一点不过会多占用一些内存。我实测下来这个computeSSE用10000行×5列的数据一次评估只需要几毫秒。PSO运行100代、30个粒子总耗时大概在几十秒级别完全在可接受范围。3.3 参数配置与K值确定PSO本身的参数不是玄学但有几个需要根据实际数据调整的地方。种群规模N。N太小容易漏掉优质解空间N太大计算量线性上涨但收益边际递减。我在这个项目里试过20、30、5030粒子和50粒子的最终SSE相差不到1%运行时间却差了接近一倍。建议先从N30开始如果发现gbest适应度曲线后期还在大幅下降说明粒子多样性不足或迭代次数不够再往上加。惯性权重w的线性递减策略是PSO的常见设定。w大时粒子飞行速度快、全局探索能力强w小时粒子趋于精细局部搜索。从0.9线性降到0.4正好覆盖前期探索、后期收敛的过程。如果你想进一步优化可以改成非线性递减或者自适应惯性权重但对当前Kmeans初始化这个场景线性递减已经够用。K值的确定是整个流程里最影响结果的一步。我用的是SSE肘部法配合轮廓系数交叉验证。具体做法分别运行K2到8的PSO-Kmeans把每次的最终SSE值画成折线图。通常随着K增加SSE单调下降但下降速度会在某个K值之后明显放缓这个点就是“肘部”。同时使用silhouette函数计算轮廓系数轮廓系数越高说明簇内紧凑度越好、簇间分离度越高。肘部法会倾向做参考轮廓系数在做裁决。我实验中最优K在4到5之间考虑业务落地的可解释性最终选择K4。3.4 结果评估与业务映射聚类完成后评价不能只看算法指标。我先把四种聚类模式的SSE和轮廓系数做了对比方案SSE轮廓系数单次运行时间普通Kmeans随机初始化最优一次356.820.3120.8sKmeans自带默认348.610.3311.2sPSO-Kmeans本文方案341.150.36826.5s这个结果符合预期PSO-Kmeans比普通Kmeans的SSE降低了约4.4%轮廓系数从0.312提升到0.368聚类质量有明显改善。代价是运行时间从秒级到几十秒但对于离线分析场景这个时间完全可控。算法指标之外真正让我认可这个方案的是业务侧的可解释性。K4时四个簇的典型负荷曲线特征非常清晰第一簇用户工作日白天负荷低、晚间负荷集中典型的上班族作息第二簇用户全天负荷水平低但稳定多为老人留守家庭空调使用率低第三簇用户峰谷差极大晚高峰负荷显著高于其他时段可能有电动车充电或者电热水器集中使用第四簇用户全天负荷水平高、峰谷差小疑似有分布式光伏或者全天候生产性用电。这种聚类结果可以直接落到营销策略上。比如对第三簇用户可以推荐错峰充电电价对第四簇用户不适合做需求响应邀约因为本身负荷刚性太强。4. 实操中踩过的坑与排查实录项目从零到跑通前前后后折腾了两周。下面这些坑不是从文档里看来的都是实际调试中遇到过的值得认真记一笔。4.1 K值怎么定才靠谱最初我图省事只看了SSE肘部图就定了K5结果聚类出来的第二个簇和第三个簇中心距离极近两个簇的特征几乎一致。这就是只看单一指标的下场。后来我把轮廓系数也加了进来发现K5时轮廓系数反而比K4低不少说明5类划分出现了簇重叠。同时我还计算了戴维斯-布尔丁指数作为参考DBI越低代表类内分散度小、类间分散度大。三项指标一起看K4的结论就非常稳定了。这里给一个建议K值不要只看算法指标。算法指标再好如果聚出来的4个簇在业务上无法描述、无法对应到客户群体类型这个聚类就是失败的。我在实验中把K3到K6的结果都打印成负荷曲线图拿给业务同事看他们最认可的是K4的结果理由是“每类用户都能对上一个可以描述的行为特征”。4.2 PSO搜索效果差、收敛慢的排查有一次实验PSO迭代了200代但适应度曲线从第60代开始就平了而且最终SSE比普通Kmeans还高。这个现象第一反应以为是PSO代码写错了后来排查完发现是数据标准化没有做负荷绝对值特征在距离计算中权重太大导致搜索空间里的大部分区域都是无效区域。标准化之后问题直接就消失了适应度曲线下降正常了。如果做了标准化还是效果差下一步就检查种群初始化和边界设置。我在最初版本里把粒子边界设成了[0,1]没有根据实际数据范围调整结果粒子普遍在边界附近越界被拉回搜索效率极低。后来改成lb repmat(min(data_z), 1, K)先把粒子约束在合理范围收敛速度明显改善。还有一个典型问题是PSO陷入局部最优。判断特征是多次运行结果里gbest_fit每次都一样但数值偏大。这说明粒子群没有很好地覆盖解空间。处理办法有三个方向一个是加大种群规模到50左右另一个是调整惯性权重下限从0.4降到0.3保持后期一定的探索能力还有一个是在速度更新中给速度设置最大限幅避免粒子飞得太猛直接越过优质解区域。我实测中调整权重下限最有效。4.3 聚类结果与真实数据画像对不上有一版聚类结果出现了一个让我困惑的簇这个簇包含的用户数量非常大而且簇心曲线形状跟总体均值几乎一样。排查之后发现问题出在特征工程我只用了工作日特征和负载率休息日的用电行为信息没有进入模型导致很多行为差异只在休息日展现的用户被划到了一起。加上休息日均值特征和峰谷差率特征之后这个现象消失了。这说明聚类效果的上限其实是由特征工程决定的算法只能在给定特征下尽量优化。特征没有区分度再好的优化算法也白搭。另一个和业务对不上的原因是直接用PCA降维后聚类再回头解释业务。PCA降维之后再聚类簇在原始特征上的业务可解释性往往变差。如果需要降维建议降维后仍然回到原始特征上统计每个簇的业务指标而不是只看降维空间的坐标。4.4 Matlab代码层面的几个注意点第一随机数种子。PSO和Kmeans都涉及随机初始化为了让实验可复现务必在脚本开头加rng(42)或者任意固定的种子。很多科研实验把随机性误当作改进效果最后被人复核时对不上这是最尴尬的。第二避免样本级的双重循环。有人喜欢用两重for循环逐样本计算距离数据量小的时候看不出问题几万条样本的时候那个速度真的会怀疑人生。使用向量化矩阵运算代码短、速度快、可读性也更好。第三大数据量的加速技巧。如果你的样本超过几十万行PSO的适应度计算会开始吃力。此时可以先随机抽样一部分代表性样本做PSO寻找初始中心再用全部数据做最终Kmeans。抽样时注意分层抽样保证覆盖全天的负荷特征不要只按序号机械抽。第四聚类完成后要把簇标签回写到原始数据表。我一开始只保存了最终的聚类中心后来要分析每个簇的客户特征时还得重新跑一遍分配白白浪费了一次完整的Kmeans迭代。正确做法是在聚类结束后直接用[idx, C] kmeans(...)把idx存下来匹配到用户ID上后面画图、统计、出报告都不用再重跑。回到最初的问题——居民用电行为聚类值不值得上PSO这种优化算法如果只是写课程报告自带Kmeans多跑几次选最优勉强够用。但如果是做真实业务的研究结果稳定性和可复现性远比那几秒运行时间更重要。PSO-Kmeans用适度的计算代价换取了聚类质量的大幅提升让整个分析链条从“碰运气”变成了“可预期”。在做负荷画像、需求响应潜力和分时电价策略支撑这类场景时这套流程是值得复制到生产环境里的。我个人的体会是算法组合不必追求复杂但每个环节为什么要这样设计必须自己能讲明白——这比代码本身更能决定项目的成败。