做用户侧数据分析这几年有一个问题几乎每周都会被问到这一万多户居民到底该怎么分群看24小时负荷曲线每一户都不一样但业务上不可能一户一户去设计方案必须归纳成几类典型行为。我最后用的方案是粒子群算法优化Kmeans聚类Matlab一套脚本跑通既解决了标准Kmeans对初始中心太敏感的问题又能把结果拿去做业务解读。这套思路适合刚接触聚类算法、想做用电行为分析的研究生也适合电力营销、需求响应岗位的从业者直接参考。数据清洗、特征构造、PSO和改进Kmeans的融合方式、代码怎么拆、参数怎么调我都会讲到尽量给出一套能直接复现的流程。1. 为什么居民用电行为分析绕不开聚类1.1 业务视角用电行为数据到底能给谁用智能电表普及以后电网侧能拿到的数据粒度已经非常细了。常见的是按15分钟或1小时冻结一次电量一天下来就是96点或24点负荷曲线。数据量一多难题就来了这些数字本身不能直接指导业务必须转化成“哪类用户、什么特征、适合什么策略”的结论。比如同一个台区里A用户晚上7点到10点负荷很高B用户白天9点到17点负荷很高两户的平均日用电量可能差不多但它们在分时电价、需求响应、有序用电场景里的价值完全不同。如果用“平均负荷”这样的单一指标去描述信息的损失太大。聚类就是用来解决这个问题的把N个用户按负荷曲线特征分成K个簇每个簇内部用户在用电模式上尽量相似簇与簇之间尽量不同。聚类结果一旦出来用途非常直接。营销侧可以把“晚高峰型”用户单独挑出来设计峰谷套餐配网侧可以用不同簇的负荷曲线叠加估算台区尖峰负荷需求响应侧则可以针对“可转移负荷占比高”的簇做削峰邀约。所以这不是一个纯算法题而是一个从数据到业务决策的完整链路。1.2 算法视角为什么这类表格型数据普遍用Kmeans居民用电负荷数据经过特征提取后形态上就是一个二维矩阵行是用户列是特征比如24个时刻的功率值、峰谷电量占比等。在这种低维度、稠密矩阵上做无监督分群Kmeans几乎是首选。原因很直接。第一Kmeans的时间复杂度接近O(N·K·D)在几万用户、几十个特征的数据集上跑起来很快第二Matlab里内置了kmeans函数还有silhouette等评价函数工程实现成本很低第三Kmeans输出的簇中心本身就是“平均负荷曲线”可以直接画出来看可解释性很强。但这不代表Kmeans可以直接用。我用标准Kmeans跑同一份居民负荷数据时连续试了五六次每次聚类结果都不一样有时候轮廓系数高一些有时候低不少。问题就出在Kmeans的迭代机制上它对初始聚类中心的选择极其敏感初始中心选得不好算法就收敛到一个局部较优解而不是全局理想解。2. Kmeans的固有缺陷以及PSO能帮什么忙2.1 Kmeans迭代求解的本质与三大软肋先看一下Kmeans在做什么。它本质上在最小化一个目标函数也就是所有样本到它所属簇中心的距离平方和J sum(i1..K) sum(x in C_i) || x - u_i ||^2其中u_i是第i个簇的中心C_i是第i个簇的样本集合。算法流程很简单先随机选K个初始中心然后每个样本就近归类再重新计算每个簇的均值作为新中心重复两步直到中心不再明显变化。这个过程有一个很直观的类比一个人站在山区里想找最低点Kmeans只会沿着当前能看到的下坡方向走走到一个山谷底部就停了但这个山谷未必是整片山区的最低点。初始中心不同相当于把这个人放在不同位置最后掉进的山谷可能完全不同。具体到居民用电数据问题集中在三点一是初始中心敏感。如果随机抽出来的K个用户恰好都集中在某一种用电模式里比如全是低电量用户那么高电量用户所在的簇可能没人去“占领”最后聚类结果会严重失衡。二是容易陷入局部最优。Kmeans每一步都在下降目标函数但每一步都是贪心的局部调整缺少跳出局部谷底的能力。三是K值还需要单独确定。这个不是Kmeans本身能解决的通常靠肘部法则、轮廓系数或业务经验去判断。2.2 粒子群算法为什么适合做Kmeans的前端优化粒子群算法PSO是模拟鸟群觅食行为的群体智能优化算法。它不要求目标函数可导也不要求问题有连续解析解只需要能对每个候选解算出一个适应度值就能在这个解空间里搜索。PSO的核心描述很简单一群粒子在解空间里飞行每个粒子都有一个位置代表候选解和速度代表移动方向和步长。每轮迭代粒子会根据自己的历史最优位置pbest和整个群体的历史最优位置gbest来更新速度再更新位置。速度更新公式是v(t1) w·v(t) c1·r1·(pbest - x(t)) c2·r2·(gbest - x(t)) x(t1) x(t) v(t1)这里w是惯性权重控制粒子保持原来运动趋势的程度c1和c2是学习因子分别控制个体认知和社会认知的权重r1和r2是两个介于0和1之间的随机数。把这个机制对比Kmeans的“贪心下山”PSO的特点是每个粒子不仅仅在局部下降它还会被“拉向”自己曾经到过的好位置和群体发现的好位置。所以在搜素前期粒子比较分散相当于一群人在山区不同位置同时找低洼处搜素后期大家逐渐靠拢到目前发现的最低点附近再在这个区域精细搜索。这里有一个很关键的思路转换要把PSO用在Kmeans上必须把“聚类问题”转换成“优化问题”。聚类中心u1, u2, ..., uK本身就是一组参数我们把它们拼接成一个向量就是PSO里一个粒子的位置。比如特征维度是D簇数是K那么一个粒子的维度就是K×D。对这个粒子解码还原成K个簇中心然后对所有样本计算距离平方和得到的J值就是适应度。J越小说明这组簇中心越好。2.3 两种融合路径质心初始化 vs 质心精调实际做实验时PSO和Kmeans的融合有两条路径可选概念上容易混分开说清楚。第一条路径是把PSO当作“找初始中心”的工具。先用PSO在K×D维空间里搜索一组比较好的簇中心然后把这组中心作为Kmeans的初始质心再让Kmeans继续迭代收敛。这样做的好处是Kmeans保留了它成熟高效的局部聚类能力PSO则负责提供一个更合理的起点。第二条路径是完全用PSO代替Kmeans的迭代。也就是在PSO每轮解出当前最优粒子后直接把所有样本按最近距离原则划分到对应簇中不再调用Kmeans做局部迭代。这种做法理论上更符合“优化”的语义但在实际数据上往往没必要。因为Kmeans本身在给定较好初始中心后局部收敛非常快再让它多走几步并不会破坏PSO的搜索结果反而能把边界样本的分簇质量修得更细一点。我在代码里采用的是第一条路径。整体框架是数据标准化后用PSO搜索最优质心位置再把最优质心喂给Matlab内置kmeans函数作为Start参数最后输出标签和轮廓系数。这样既稳定又符合大多数论文和业务报告里“PSO优化Kmeans”的表述。3. 从原始计量数据到可聚类特征数据准备与预处理3.1 数据清洗和特征怎么构造原始数据不能直接进算法。智能电表采回来的数据常见问题有某天采集失败导致整行缺失、个别时点出现跳变比如负数、超过表计倍率的极大值、用户搬家导致一段时间无数据等等。我的处理习惯是先按用户ID和日期排序对单条缺失时点做线性插值如果某用户一天的完整采集点数不够80%或者连续多日没有数据直接剔除这个用户对超过每日平均负荷3倍以上且持续多个时点的异常值用前后两天的同一时点值替换。清洗之后面临的第二个问题是特征怎么选。最直观的做法是把一天的24小时平均功率作为24个特征。这个方式信息最全聚出来的簇可以直接画平均日负荷曲线业务上非常好解释。缺点是特征维度高PSO粒子维度K×D会变大计算量上升。我建议在24维负荷曲线基础上根据业务目标再补几类关键特征特征名计算方式业务含义日均负荷全天功率均值用户总体用电水平峰时电量占比峰时段电量 / 全天电量用户负荷是否集中在高峰时段谷时电量占比谷时段电量 / 全天电量用户对谷电时段的利用程度峰谷差率(峰值功率 - 谷值功率) / 峰值功率负荷波动程度夜间负荷占比23点到次日6点电量 / 全天电量是否存在夜间刚需负荷负载率平均功率 / 最大功率负荷曲线的平坦程度最大负荷出现时刻一天中功率最大值所在小时用户高峰习惯这些特征本质上都是对24点曲线的压缩能降低维度也方便后期解释。如果原始数据没有峰谷时段定义可以按本地分时电价时段来切或者先用固定规则峰时段取8点到22点谷时段取22点到次日8点。3.2 标准化和K值估计特征构造完成后标准化这一步不能省。Kmeans基于欧氏距离计算相似度如果某个特征的量纲特别大距离就会被它主导。比如日均负荷是几千瓦峰谷差率是0到1之间的小数如果不处理几千瓦的特征在距离计算里会直接把其他特征压死。Matlab里直接调用zscore函数即可把每个特征变成均值为0、标准差为1的序列。标准化后的矩阵记为X每一行一个用户每一列一个特征后续聚类和PSO优化都用X。K值的选择可以用两种方式交叉验证。一种是肘部法则画出K从2到10时目标函数J的变化曲线找拐点另一种是轮廓系数取轮廓系数最大的K值。但在实际业务里我一般不会盲信这两个指标而是把业务约束也放进来。比如电网项目往往希望分3到5类因为太多类业务执行不过来太少类又区分不出差异。多数情况下K4是一个比较合理的起点后期再根据聚类质量和业务反馈微调。4. Matlab代码实现PSO-Kmeans从建模到出图4.1 主脚本框架与参数设置下面这段代码是我整理出来的核心主脚本。为了方便讲解假设已经清洗完数据得到了一个N×24的负荷矩阵loadData每一行对应一个用户的24小时平均功率。% PSO_Kmeans_demo.m % 基于粒子群算法优化Kmeans聚类的居民用电行为分析 % 输入: loadData N×24 原始日负荷曲线 % 输出: idx 聚类标签 % silAvg 平均轮廓系数 rng(1); % 固定随机种子保证可复现 loadData fillmissing(loadData, linear); % 线性插值补缺失 X zscore(loadData); % 标准化消除量纲影响 N size(X, 1); D size(X, 2); K 4; % 聚类数可由肘部法则/业务确定 % PSO参数 nP 30; % 粒子数 maxIter 60; % 迭代代数 wStart 0.9; % 惯性权重初始值 wEnd 0.4; % 惯性权重结束值 c1 1.5; % 个体学习因子 c2 1.5; % 社会学习因子 dim K * D; % 粒子维度: K个簇中心PSO参数不是随便拍的。粒子数一般取20到40特征是24维、簇数取4的时候粒子维度是9630个粒子在96维空间里搜索覆盖能力已经够用。迭代代数取50到80太多收益递减太少容易搜索不充分。w从0.9线性降到0.4是为了实现“前期大范围探索、后期精细收敛”这是PSO调参里比较通用的做法。初始化粒子位置时我不用纯随机数而是从样本里随机抽K个用户把这K个用户的特征向量拼接成一个粒子。这样初始质心一定落在真实数据范围内不会出现一开始就跑到特征空间外部、距离完全失真的情况。pos zeros(nP, dim); vel zeros(nP, dim); pbestPos zeros(nP, dim); pbestVal inf(nP, 1); for p 1:nP idxSample randperm(N, K); % 随机抽K个用户 initCenter X(idxSample, :); % K×D pos(p, :) initCenter(:); % 拼接成1×(K*D) end gbestPos pos(1, :); gbestVal inf;4.2 适应度函数与PSO主循环适应度函数是整个优化过程的核心。给定一个粒子位置先reshape成K×D的质心矩阵然后计算所有样本到最近质心的距离平方和。这个值就是PSO要最小化的目标。function [fit, label] psoKmeansFitness(X, pos, K, D) % 把粒子位置解码成K个聚类中心返回距离平方和与最近中心标签 center reshape(pos, K, D); N size(X, 1); label zeros(N, 1); distSum 0; for i 1:N dist sum((X(i, :) - center).^2, 2); % K×1 [minD, minIdx] min(dist); label(i) minIdx; distSum distSum minD; end fit distSum; end这个函数写法比较直观但注意循环遍历所有样本如果用户量到几十万单次适应度计算会偏慢。单机实验几万用户以内问题不大后续做大规模数据时可以用矩阵化重写后面我会说到。PSO主循环的逻辑是先遍历所有粒子计算适应度更新每个粒子的个体最优和全局最优再按速度公式更新粒子的速度和位置位置更新后要做边界截断防止质心飞到特征空间外面。xMin min(X); % 各特征最小值 xMax max(X); % 各特征最大值 for iter 1:maxIter w wStart - (wStart - wEnd) * iter / maxIter; % 评估当前所有粒子 for p 1:nP [fit, ~] psoKmeansFitness(X, pos(p, :), K, D); if fit pbestVal(p) pbestVal(p) fit; pbestPos(p, :) pos(p, :); end if fit gbestVal gbestVal fit; gbestPos pos(p, :); end end % 更新速度和位置 for p 1:nP r1 rand(1, dim); r2 rand(1, dim); vel(p, :) w * vel(p, :) ... c1 * r1 .* (pbestPos(p, :) - pos(p, :)) ... c2 * r2 .* (gbestPos - pos(p, :)); pos(p, :) pos(p, :) vel(p, :); % 边界截断 pos(p, :) max(pos(p, :), repmat(xMin, 1, K)); pos(p, :) min(pos(p, :), repmat(xMax, 1, K)); end end迭代结束以后gbestPos就包含了PSO找到的最优质心。注意optimum位置可能落在特征空间的边界上如果经常贴边通常说明K值或者特征构造有问题值得回头查数据。4.3 把最优中心交给Kmeans并可视化得到gbestPos之后我习惯再交给kmeans函数跑一次。这样做的原因是PSO搜索到的是“质心放在哪里整体距离更小”而Kmeans在给定这个起点后还能通过局部迭代把质心调整到更精细的位置。bestCenter reshape(gbestPos, K, D); % 把PSO找到的中心作为Kmeans初始质心 idx kmeans(X, K, Start, bestCenter, ... MaxIter, 300, ... Replicates, 1, ... Distance, sqeuclidean); % 轮廓系数 sil silhouette(X, idx); silAvg mean(sil); fprintf(平均轮廓系数: %.4f\n, silAvg);这一步有一个容易踩的坑kmeans函数里的Start参数在传入矩阵时矩阵每一行就是一个初始中心所以传入的必须是K×D矩阵而不是拼接后的向量。因此必须先reshape。结果可视化我一般分两张图。第一张是聚类标签下的PCA降维散点图适合快速看簇的重叠程度第二张是每一类用户的平均日负荷曲线适合业务解读。% 图1: PCA降维后看聚类分布 [coeff, score] pca(X); figure; gscatter(score(:,1), score(:,2), idx); xlabel(PC1); ylabel(PC2); title(PCA降维后的聚类分布); % 图2: 每类的平均日负荷曲线 figure; for k 1:K subplot(2, 2, k); plot(mean(loadData(idx k, :), 1), LineWidth, 1.5); xlabel(小时); ylabel(平均功率/kW); title([第 num2str(k) 类用户平均日负荷曲线]); grid on; endissueloadData还在前面定义过所以在subplot里能用。如果原始数据不是24点而是96点把横轴刻度改成15分钟间隔即可。5. 实验效果怎么看稳定性、轮廓系数与业务解读5.1 评价一个聚类结果不能只看目标函数我试过很多次单独看目标函数J会有一个问题PSO把J降得很低但聚类结果未必是业务上最舒服的。因为J本质上衡量的是“紧凑程度”没有充分考虑到簇之间的分离程度。所以做效果对比时我至少看三个角度。第一个是稳定性。同样一份数据把随机种子换几个跑多次看簇中心和标签是否变化。标准Kmeans如果随机初始化几次结果经常忽大忽小有时甚至某个簇只有几个人。PSO-Kmeans在这点上改善非常明显因为PSO的搜索不是纯随机撒点而是有记忆、有方向的搜索最终质心落点更稳定。第二个是轮廓系数。Matlab里直接算值域在-1到1之间越接近1说明样本离自己簇的中心越近、离邻近簇越远。我个人经验是居民用电行为数据能做到0.2到0.35就算不错毕竟用户行为本身是连续过渡的没有绝对清晰的边界。如果有人报告说轮廓系数0.8大概率是特征或者K值选择有问题比如把完全重复的数据当成了特征。第三个是运行时间。PSO-Kmeans增加了适应度计算这一步时间开销肯定比单次Kmeans大但一般不会失控。拿1万户、24维特征、30个粒子、60代迭代来说跑一趟通常在几十秒到两三分钟取决于机器和循环写法。如果数据到了几十万用户就需要考虑矩阵化加速或者特征降维。我自己的经验是把普通Kmeans单次、Kmeans重复10次取最优、PSO-Kmeans放在同一份数据上对比时差距主要体现在波动性上。普通Kmeans不同随机种子跑出来的轮廓系数可能从0.18到0.26之间乱跳重复10次取最优会稳定在较高值但用户耗时要乘10PSO-Kmeans的轮廓系数我一般能稳定在0.25左右而且耗时通常低于重复10次Kmeans的总时间。5.2 聚类结果的行为画像怎么解读聚类算法只是中间产品真正交付给业务的是“每类用户是谁、有什么特征、该怎么运营”。对居民用户来说K4是比较常见的分法我在实际项目里经常会看到下面几类画像。第一类是晚高峰主导型。平均负荷曲线在18点到22点有明显高峰白天低谷偏低。这类用户大多是上班族傍晚回家后集中用电包括做饭、洗浴、娱乐负荷。峰谷差率通常比较大是分时电价和削峰填谷的重点目标。第二类是全天平稳型。曲线整体偏高且平坦负载率很高夜间负荷也不低。这类用户可能是家里有老人小孩常驻或者有持续运行的设备。它们用电量基数大参与需求响应的空间也大但响应意愿不一定高需要考虑不影响正常生活为前提的互动策略。第三类是白天工作型。工作日白天负荷明显周末和夜间下降。这类用户很多是家庭作坊、小店、或在家办公人群。它们对电价比较敏感更可能响应白天的激励政策。第四类是夜间用电型。曲线在22点到次日6点之间偏高。常见原因包括电动汽车充电、蓄热式电暖器等。这类用户在谷电阶段用能比例高适合谷电套餐也是电网消纳新能源、填谷的主力。分完类以后我一般会把每类用户占总用户数的比例、户均电量、峰谷差率做一张透视表再结合空间位置画到地图上这样业务部门可以直接拿去用。6. 调参、踩坑和可复现性建议6.1 PSO参数调参经验PSO的参数不算多但每个参数对结果都有影响分享一下我的调参手感。粒子数nP特征维度越高粒子数要多一点。24维特征、K4时30个粒子够用如果特征升到40维以上建议加到50。粒子数再多时间开销增长明显但精度收益会边际递减。惯性权重w我习惯用0.9到0.4线性递减。开始阶段w大粒子速度快负责全局搜索后期w小粒子速度慢负责局部精细搜索。固定w也可以但容易在前期收敛太快或后期震荡所以动态w通常更稳。学习因子c1、c2一般取1.5到2.0并且让c1和c2相等。c2如果明显大于c1会让所有粒子过早被gbest吸引种群多样性下降容易早熟c1过大则粒子各自为政收敛变慢。遇到结果发散的情况先把c1、c2都调到1.5再观察。边界处理粒子位置越界后简单粗暴地截断到边界。注意这只是硬约束更平滑的做法是让粒子在越界方向上的速度分量归零但代码复杂度会上去。样本量不大时直接截断问题不大。6.2 Matlab实现中的常见坑第一个坑是维度拼接顺序。initCenter(:)是按列取元素reshape成K×D时也是按列填充。如果中间不小心用了initCenter(:)或者转置质心矩阵的数据布局会完全错乱聚类结果自然不对。调试时可以用一个很小的测试矩阵跑一遍适应度函数打印reshape后的center确认一下。第二个坑是标准化对象搞错。特征工程阶段标准化后的X进聚类但画平均日负荷曲线时要用原始的loadData。如果把标准化的曲线拿去做均值画出来的形状虽然一样但Y轴含义变成了“标准差倍数”业务人员根本看不懂。第三个坑是K值选错。PSO可以把初始质心优化得很好但如果K本身不符合数据结构特点结果依然没有业务价值。比如把K设成2所有用户被硬分成“高用电”和“低用电”虽然看起来轮廓系数还行但行为细节全丢K设成7会出现两个空簇或者两个几乎重叠的簇。建议先在标准化后的数据上跑一遍肘部法则再结合业务确定K。第四个坑是随机种子。PSO、Kmeans、PCA里都有随机成分科研和报告里只要强调复现性就一定要在脚本开头加rng(1)之类的固定种子。不然今天跑出来四类明天跑出来四类但用户编号完全变了后面做案例分析会很被动。第五个坑是样本量变大后的性能问题。前面那个适应度函数为了好懂用了for循环1万用户没问题但如果是10万用户每次适应度计算就要做10万次距离计算60代乘以30个粒子总计算压力就上来了。这时候建议重写适应度函数把距离计算改成矩阵运算计算N×K距离矩阵再对行求最小值和索引速度能快一个数量级。更激进的做法是先PCA降到5到8维再做PSO粒子维度一下子降下来收敛速度会快很多。6.3 从K4开始跑通第一版最后给一个实操建议第一次跑的时候不要急着做复杂的特征工程。直接拿24点负荷曲线标准化K固定为4PSO参数用上面代码里的默认值先看聚类结果和平均日负荷曲线是否合理。如果分出来的某类曲线和业务常识完全对不上先查数据清洗再查K值然后才是调PSO参数。我曾经在真实数据上遇到过一个典型问题有一类用户的平均曲线在凌晨5点突然蹿高后来查数据发现是热水器定时加热这类用户聚在一起其实是“凌晨加热型”并不是数据错误。这就是聚类的价值——它会把你不曾预期但有真实业务逻辑的规律挖出来。这套流程跑顺之后后面要做的扩展无非是三件事一是把适应度函数向量化支撑更大规模的用户量二是尝试用轮廓系数、DBI指标代替距离平方和作为适应度函数看哪种目标更符合业务预期三是把聚类结果与用户档案、电价方案做关联分析让“分群”真正变成“运营动作”。在居民用电行为分析这个方向上算法从来不是终点把聚类结果用起来才是终点。