PSO优化Kmeans聚类：居民用电行为分析的Matlab实践

开头我想先讲一个场景前阵子我拿到一份某小区居民用电负荷数据要做用户行为分类第一反应就是上Kmeans。结果试了几次同样的K值聚类结果每次都不一样甚至有一次把明显是高耗能的用户和一户几乎不用电的老人家庭分到了一组。问题就出在Kmeans的初始中心是随机选的算法一进局部最优就出不来。后来我改用粒子群算法PSO去优化Kmeans的初始中心先把全局搜一遍再把结果作为Kmeans的起点做精细聚类稳了很多轮廓系数也明显提升。这篇文章就完整记录一下这套思路以及对应的Matlab实现包含数据特征处理、PSO参数整定、核心代码、结果解读和踩坑记录想直接复现的同学照着抄即可。这个项目特别适合正在做用电行为分析、用户画像、负荷预测分组的同学参考也适合刚接触群智能算法优化聚类的人拿来做入门练手。核心就一句话用粒子群算法替Kmeans找一个好的起点让聚类结果不再飘忽不定。原理不复杂Matlab里全部代码加起来也就一两百行不需要额外工具箱统计和优化相关的基础函数就够用。1. 项目思路与方案选型1.1 普通Kmeans到底哪里不够用Kmeans是最经典的划分式聚类算法它的目标函数非常明确最小化各样本点到其所属簇中心的距离平方和也就是常说的簇内误差平方和SSE。算法执行时先随机初始化K个中心然后反复执行“分配样本→更新中心”这两个步骤直到中心不再变化。问题恰恰出在“随机初始化”上。Kmeans对初始中心非常敏感初始点选得不好算法很容易收敛到一个局部最优解而不是全局最优解。比如下面这几种情况在实际数据里特别常见两个初始中心同时落在同一个真实簇里另一个真实簇始终没有被抢占最后两个簇被迫分裂或合并初始中心落在离群点上导致一个簇只有一个样本其他簇却拥挤不堪因为迭代本质是贪心的一旦某次分配确定了后续很难翻盘。有人会说用Kmeans可以解决这个问题它确实比纯随机初始化好它会让初始中心之间尽量分散但不能保证全局最优。而且Kmeans只做了一次启发式初始化后面依然是逐步贪婪收敛的过程。换句话说它改善了起点却没有改变Kmeans“贪心收敛”的本质。1.2 为什么选粒子群算法来优化粒子群算法是一种群智能全局优化算法模拟鸟群觅食行为每个粒子代表一个候选解通过个体经验和群体经验不断调整自己的位置。它的核心优势在于不依赖梯度信息、实现简单、参数少非常适合与Kmeans这种“目标函数明确但没有解析解”的问题结合。在这个项目里我把Kmeans的初始中心编码成粒子的位置把Kmeans的SSE作为适应度函数。粒子群在解空间并行搜索寻找使SSE最小的那组初始中心。找到后把这组中心交给Kmeans再次迭代收敛。这就是典型的“全局粗搜索 局部精搜索”组合策略。用这种方式聚类结果的稳定性会大幅提升。因为粒子群本身是多起点并行搜索并且粒子之间有信息共享某个粒子找到的好位置会通过gbest引导其他粒子靠拢这比单纯跑很多次Kmeans选最优更聪明也更有方向性。1.3 整体技术路线整个项目的流程我整理成一条链各位跟着走就行读取居民用户24点负荷数据真实数据或模拟数据从原始负荷曲线中提取聚类特征比如日均用电量、峰段占比、谷段占比、峰谷差率、负荷率对特征矩阵做Z-score标准化消除量纲影响用肘部法则和轮廓系数综合确定K值运行PSO优化算法得到Kmeans的初始中心适应度函数为SSE将PSO得到的中心作为Kmeans的Start参数重新聚类计算轮廓系数、簇内SSE等指标评估聚类效果绘制各类用户的平均负荷曲线解读用电行为画像。这套流程最大的好处是每个环节都是独立的出问题了可以单独排查后面遇到坑也好定位。2. 数据准备特征提取与预处理2.1 居民用电数据长什么样居民用电负荷数据最常见的格式是一张二维表行是用户列是按时间序列排列的负荷值。24点数据就是每天24个整点时刻的负荷96点数据就是每15分钟采一个点一天96个点。实际工程中15分钟采集粒度更常见但为了演示方便我先把数据重采样到24点聚类结论不受影响。需要说明的是我不打算用某个特定小区的真实数据而是构造了4类典型负荷曲线来做算法验证每一类对应一种用电习惯的用户双峰型用户早晚各有一个用电高峰典型上班族家庭早晨吹风机、电热水壶晚上电视、空调、照明集中夜间型用户用电集中在晚8点到11点属于典型的“晚睡型”用户可能是下班晚或者习惯夜间活动平稳型用户全天负荷比较均衡白天也有一定用电量常见于白天有人在家的家庭低负荷型用户全天用电量很低可能是一户独居老人或者极少用电的临时住所。每一种曲线都叠加了随机噪声模拟真实的波动。等流程全部跑通后把读取数据的部分替换成你手上的真实数据文件即可算法部分完全不用改。2.2 聚类的核心特征怎么选直接用24个时刻的负荷值做聚类不是不行但会有两个问题一是维度偏高24维数据在欧氏距离下区分度反而不如低维特征明显二是PSO优化Kmeans时粒子的维度等于K乘以特征维度如果直接拿24维特征做每个粒子的维度就是24乘K搜索空间会爆炸收敛速度大幅下降。所以我建议先提取5个工程意义明确的特征特征名称计算公式含义日均用电量全天24点负荷均值反映整体用电水平峰段占比早晚高峰时段用电量占总用电量的比例反映用电在高峰时段的集中程度谷段占比凌晨低谷时段用电量占总用电量的比例反映夜间基础用电量占比峰谷差率晚高峰均值与低谷均值的差值除以日均负荷反映负荷波动幅度负荷率日均负荷除以日最大负荷反映负荷平稳程度值越接近1越平稳这5个特征从“总量、时段集中度、波动性”三个角度刻画了用户的用电习惯彼此之间的冗余度不高解释性还很强。聚类结束后直接看每个簇在这5个维度上的均值就能写出清楚的行为画像。2.3 Z-score标准化为什么必须做如果不做标准化日均用电量这个特征的数值可能是零点几到几之间峰段占比同样在零到一之间看起来差距不大。但如果将来换成包含电量绝对值的数据量纲差异马上会显现。数值范围大的特征在欧氏距离计算中会占据主导地位聚类结果就变成了“只看某个特征聚类”这是典型需要避开的错误。Z-score标准化的公式很简单z (x - mean) / std也就是每个特征减去均值后除以标准差变换后每个特征的均值为0标准差为1所有特征在距离计算中处于平等地位。Matlab里直接调用zscore函数一行就能搞定。需要注意一点标准化用到的均值和标准差必须在全部样本上计算不能对每个样本单独标准化。样本单独标准化后所有样本的均值都变成0等于把相对差异抹掉了聚类就没有意义。3. PSO-Kmeans核心原理与参数设计3.1 粒子编码方式PSO里的每个粒子代表Kmeans的一组初始聚类中心。假设特征维度是d聚类数是K那么每个粒子是一个长度为K×d的一维向量前d个元素是第1个簇的中心接下来d个元素是第2个簇的中心以此类推。在Matlab里我习惯把一维向量先reshape成d行K列的矩阵再转置成K行d列这样每个行向量就是一个聚类中心。这里有一个特别容易踩的坑reshape是按列填充的如果直接写成reshape(pos, K, d)排列顺序会和预期完全反掉。正确写法是C reshape(pos(i, :), d, K);先按列填充成d行K列保证每一列是一个中心点转置后得到K行d列一行一个中心。3.2 适应度函数设计适应度函数直接决定优化方向这里我选择Kmeans的标准目标函数SSESSE Σ 各簇内所有样本到簇中心的欧氏距离平方之和SSE越小说明聚类越紧凑簇内一致性越好。用SSE而不是轮廓系数做适应度的原因很务实轮廓系数的计算复杂度更高而且它的评价逻辑和Kmeans的迭代目标不完全一致用它当适应度可能导致“按A指标优化、按B指标收敛”的错位。最终评估时再用轮廓系数做外部验证这个组合最合理。适应度函数的具体实现也很简单Matlab自带的pdist2可以一次算出所有样本到所有中心的距离矩阵取每行最小值得到所属簇编号再按簇汇总距离平法和。3.3 核心参数设置PSO算法的四个关键参数是粒子数、迭代次数、学习因子和惯性权重。我常用的参数范围如下参数设置值说明粒子数 N30维度不高时20到30足够多了计算量大但提升有限最大迭代次数 MaxIt100一般50到100就能收敛看收敛曲线调整学习因子 c1, c22.0, 2.0经典取法个体经验和群体经验权重均衡惯性权重 w0.9线性递减到0.4前期重在探索后期重在收敛惯性权重的线性递减是我实测下来最稳定的做法。公式是w w_max - (w_max - w_min) × iter / MaxIt这样前期权重高粒子飞行速度快、探索范围广迭代后期权重降低粒子逐渐围绕最优区域精细搜索。速度限幅也是必须的。如果不限制速度粒子可能直接飞出边界位置更新的步长过大导致搜索失去方向性。一般把最大速度设为每维搜索范围的20%左右。对这个项目来说标准化后数据基本落在[-3, 3]区间内vmax设成0.8到1比较合适。3.4 粒子位置初始化的小技巧粒子位置初始化有一个比均匀随机采样更好的方法对每个粒子从样本集中随机抽取K个样本直接把这K个样本的特征向量作为该粒子的初始位置。这样做的好处是初始中心一定落在实际样本点上不会出现空簇搜索起点更合理。有人会问这和Kmeans随机初始化有什么区别区别在于Kmeans只初始化一组中心然后一路贪心走到底而PSO初始化的是N组中心并且这N组中心之间通过pbest和gbest建立信息共享可以在迭代中互相学习。4. Matlab完整代码实现与关键点说明4.1 数据生成与特征提取部分完整代码我按模块分开贴注意下面的代码可以直接复制保存成一个m文件运行。%% 基于PSO优化Kmeans的居民用电行为分析 clc; clear; close all; rng(42); % 固定随机种子保证结果可复现 % 1. 生成模拟负荷数据24点 numUsers 100; hours 24; loadData generate_load_profile(numUsers, hours); % 2. 特征提取 feat zeros(numUsers, 5); feat(:, 1) mean(loadData, 2); % 日均用电量 feat(:, 2) (sum(loadData(:, 9:11), 2) sum(loadData(:, 18:21), 2)) ./ sum(loadData, 2); % 峰段占比 feat(:, 3) sum(loadData(:, 1:7), 2) ./ sum(loadData, 2); % 谷段占比 feat(:, 4) (mean(loadData(:, 18:21), 2) - mean(loadData(:, 1:7), 2)) ./ mean(loadData, 2); % 峰谷差率 feat(:, 5) mean(loadData, 2) ./ max(loadData, [], 2); % 负荷率 % 3. Z-score标准化 featZ zscore(feat);代码里的时间段说明一下索引1到7大致对应0点到6点属于低谷段9到11对应早高峰18到21对应晚高峰这两个区间合并作为峰段。不同地区的峰谷时段划分可能不一样做真实数据时先读清楚采集点对应的时间标签再改索引不要照搬。模拟数据生成函数放在文件末尾function ld generate_load_profile(N, h) % 生成4类典型居民负荷曲线并叠加随机噪声 t 0:h-1; base zeros(4, h); % 双峰型早8点小峰 晚8点大峰 base(1, :) 0.30 * exp(-((t - 8).^2) / 1.5) 0.55 * exp(-((t - 20).^2) / 3); % 夜间型晚9点为主 base(2, :) 0.35 * exp(-((t - 21).^2) / 4) 0.15 * exp(-((t - 13).^2) / 8); % 平稳型全天较均衡14点左右略高 base(3, :) 0.35 * ones(1, h) 0.20 * exp(-((t - 14).^2) / 6); % 低负荷型几乎无高峰19点少量用电 base(4, :) 0.10 * ones(1, h) 0.15 * exp(-((t - 19).^2) / 5); base base * 1.2; % 统一放大到合理功率范围 typeIdx randi(4, N, 1); ld zeros(N, h); for i 1:N ld(i, :) base(typeIdx(i), :) .* (0.8 0.4 * rand(1, h)); end end4.2 确定K值肘部法则加轮廓系数K值怎么定不能拍脑袋。我的习惯是先跑一个K从2到6的循环计算每个K下的SSE画肘部图再结合轮廓系数一起判断。%% 肘部法则选K Kset 2:6; sseK zeros(size(Kset)); for ii 1:length(Kset) [~, ~, sumd] kmeans(featZ, Kset(ii), Replicates, 10); sseK(ii) sum(sumd); end figure; plot(Kset, sseK, -o); xlabel(聚类数K); ylabel(SSE); grid on;看SSE曲线时找“肘部”位置也就是曲线斜率从陡峭变平缓的拐点。SSE一定随着K增大而下降但下降到一定程度后收益会明显变小这个拐点就是最合适的K。与此同时计算不同K下的平均轮廓系数轮廓系数接近1说明样本离自己簇近、离其他簇远大于0.5就已经是可用的聚类结果。本项目的模拟数据天然分为4类我用K4但不代表真实数据也必然是4类真实用户群的分类数必须靠数据说话。4.3 PSO主循环实现这一部分是整个项目的核心我把粒子初始化、速度更新、边界处理、适应度计算、pbest和gbest更新全部写在一个循环里。%% PSO优化Kmeans初始中心 K 4; d size(featZ, 2); Npop 30; MaxIt 100; wMax 0.9; wMin 0.4; c1 2.0; c2 2.0; dim K * d; lb min(featZ, [], 1); ub max(featZ, [], 1); vmax 0.8 * (ub - lb); % 速度限幅 % 初始化粒子位置每个粒子随机抽取K个样本作为初始中心 pos zeros(Npop, dim); for i 1:Npop idxP randperm(numUsers, K); pos(i, :) reshape(featZ(idxP, :), 1, []); end vel zeros(Npop, dim); % 初始个体最优与全局最优 pbest pos; pbestCost inf(Npop, 1); gbestCost inf; for i 1:Npop C reshape(pos(i, :), d, K); pbestCost(i) sse_cost(featZ, C); if pbestCost(i) gbestCost gbestCost pbestCost(i); gbest pos(i, :); end end % 迭代优化 gbestSeq zeros(MaxIt, 1); for iter 1:MaxIt w wMax - (wMax - wMin) * iter / MaxIt; 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, :)); vel(i, :) max(min(vel(i, :), vmax), -vmax); % 位置更新 边界钳制 pos(i, :) pos(i, :) vel(i, :); pos(i, :) max(min(pos(i, :), repmat(ub, 1, K)), repmat(lb, 1, K)); % 计算适应度并更新pbest和gbest C reshape(pos(i, :), d, K); cost sse_cost(featZ, C); if cost pbestCost(i) pbest(i, :) pos(i, :); pbestCost(i) cost; end if pbestCost(i) gbestCost gbestCost pbestCost(i); gbest pbest(i, :); end end gbestSeq(iter) gbestCost; end % 绘制收敛曲线 figure; plot(gbestSeq, LineWidth, 1.5); xlabel(迭代次数); ylabel(SSE); title(PSO优化过程收敛曲线); grid on;代码里有两个细节提醒一下。边界钳制我直接用了截断处理也就是说粒子某一维超出上下界时直接把它拉回边界。更讲究的写法是“越界反射”让粒子回到界内后速度反向但实测两种方式对这个项目影响不大截断处理更简单也更稳。pbest和gbest的更新逻辑要严格只有cost小于当前pbest时才更新pbest然后用更新后的pbest去比较gbest。很多初学版本会在同一个循环里拿cost直接和gbest比虽然最终结果差别不大但语义上不够严谨。4.4 用PSO结果初始化Kmeans并评估PSO迭代完成后把gbest还原成K行d列的聚类中心矩阵直接作为Kmeans的Start参数。这一步就是“精搜索”环节让Kmeans在PSO给出的良好起点上快速收敛到局部最优而且这个局部最优已经非常接近全局最优。%% 用PSO最优中心初始化Kmeans C0 reshape(gbest, d, K); [lab, Cfinal] kmeans(featZ, K, Start, C0, MaxIter, 500); %% 聚类评估 sil silhouette(featZ, lab); meanSil mean(sil); fprintf(平均轮廓系数: %.4f\n, meanSil); %% 绘制聚类后的平均负荷曲线 figure; colors lines(K); for k 1:K members lab k; plot(mean(loadData(members, :), 1), Color, colors(k, :), LineWidth, 1.5); hold on; end hold off; xlabel(时刻(h)); ylabel(平均负荷(kW)); legend(arrayfun((x) sprintf(类别%d, x), 1:K, UniformOutput, false)); title(各类用户平均日负荷曲线); grid on;跑完这段代码最重要的是看两个东西一是收敛曲线是否平缓下降并趋于稳定二是平均轮廓系数能否到0.6以上。如果收敛曲线锯齿明显说明粒子速度限幅太大或者惯性权重衰减过快。5. 聚类结果分析与用户行为画像5.1 如何解读聚类结果聚类算法只能给你分组标签不会直接告诉你每组是什么人。解读的核心方法是计算每个簇在原始特征和原始负荷曲线上的均值再结合业务常识给每个簇命名。我用模拟数据跑一遍后的典型结果大致如下类别1日均用电量最高峰段占比高负荷曲线呈现明显的早晚两个尖峰。这类用户对应的就是“双峰型”上班族家里电器多早晚集中用电。类别2日均用电量中等夜间负荷曲线从18点后持续走高21点左右达到峰值。这类用户可以叫“夜间活跃型”用电习惯整体后移。类别3日均用电量中等偏高负荷率最高白天晚上差别不大。这类是“全天均衡型”白天有人在家可能是退休老人或者居家办公人群。类别4日均用电量最低全天几乎平稳且数值很低。这类是“低耗型”用电需求弱节能空间小。有了这个画像后面不管是做需求响应分组、电价套餐设计还是异常用电识别都有了业务抓手。比如“夜间活跃型”用户适合参与错峰激励而“双峰型”用户是削峰填谷的重点关注人群。5.2 为什么轮廓系数比SSE更适合做最终评估PSO优化过程中用SSE做适应度是因为它和Kmeans的迭代目标一致。但最终评估聚类效果时我更倾向于看轮廓系数。原因是SSE存在一个天然的问题K越大SSE必然越小它没法衡量“簇间分离度”。轮廓系数同时考虑了簇内紧密度和簇间分离度数值更有参考价值。轮廓系数的计算思路是对每个样本计算它到自己簇内其他样本的平均距离a再计算它到最近的其他簇所有样本的平均距离b轮廓系数等于(b-a)/max(a,b)。结果范围在-1到1之间越接近1表示聚类越合理。实际应用中平均轮廓系数0.5以上就说明聚类结构基本清晰0.7以上就是非常好的结果。我看到很多论文里硬吹自己的聚类有多好结果轮廓系数只有0.3那其实是数据本身就没分群不是算法不行。5.3 这个分析能用到哪里居民用电行为聚类不是停留在画几张曲线图就完事的科研玩具它在实际工程里至少有三种落地场景。第一是电网负荷预测。把用户按用电行为分组后针对每一类用户单独建预测模型比把所有用户混在一起预测精度要高得多。第二是需求响应。不同行为特征的用户对电价的敏感度不同聚类结果可以直接作为制定差异化激励方案的依据。第三是异常用电识别。每一类用户都有自己的正常用电模式偏离模式太远的用户就要重点关注可能是设备故障或者计量异常。6. 常见问题与踩坑记录6.1 粒子群早熟收敛怎么办所谓早熟收敛就是所有粒子在迭代初期就聚集到某个局部最优附近后续迭代基本没有改善收敛曲线前面跌得很快后面完全走平。遇到这种情况我的排查顺序是检查粒子数是否太少如果只有10个建议增加到30到40个检查惯性权重是否衰减太快w从0.9衰减到0.4是比较稳妥的经验区间不要直接从0.6开始检查速度限幅是否过小vmax太小会限制粒子的探索能力检查初始化范围是否合理粒子初始位置必须落在样本分布范围内。还有一个实操技巧如果运行结果一直不理想可以尝试增大c2相对c1的比例比如c11.5、c22.5让群体经验的影响力更大加快向好的区域靠拢。6.2 聚类结果每次运行不一样这是Kmeans类算法最常见的问题解决方案非常简单设置随机种子。在脚本开头加一行rng(42)42可以替换成任意整数就能保证每次运行结果一致。工程上做可复现实验是基本要求不固定随机种子的话你论文里的结果没法被复现审稿人一跑就露馅。6.3 K值怎么确定才靠谱不要只看肘部图也不要只看轮廓系数两个结合着看。我的操作步骤是先用肘部法画出SSE曲线找出拐点候选范围然后在这个范围内分别计算平均轮廓系数选择轮廓系数最高的K。如果两个方法指示的K不一致优先看轮廓系数因为它是更直接的聚类质量量化指标。另外确定K之后建议在业务层面验证一下K值是否合理。比如K4得到的四类用户能不能用业务语言解释清楚如果有一类用户的意义根本无法归纳说明K可能选多了。6.4 数据量大了之后跑得很慢大数据量下PSO的适应度函数计算会成为瓶颈。因为每个粒子的每次迭代都要计算一次全量样本到K个中心的距离。如果样本数有几万甚至几十万跑100次迭代会非常慢。我的处理办法是距离计算改成向量化形式避免循环内逐样本计算文中代码的pdist2就是向量化的先对样本做一次粗聚类或随机抽样用代表样本参与PSO优化得到初始中心后再用全量数据做Kmeans限制PSO的最大迭代次数很多情况下30到50次迭代就已经收敛不需要跑满100次。6.5 真实数据缺测值怎么处理模拟数据不存在缺测真实数据一定会遇到。居民负荷数据最常见的缺测形式是某几个时刻点在数据库里是空值或0值。处理思路是先看缺测比例如果单用户有效数据不足80%直接删除该用户如果只是个别时刻缺测用前后相邻时刻的均值插值填充或者用同类型用户的同期均值填充。千万不要把0值当缺测直接删低谷时段负荷为0是真实存在的正常现象。6.6 常见问题速查表现象可能原因解决方法PSO收敛曲线锯齿波动大速度限幅过大调小vmax比如设为搜索范围的10%到20%聚类结果每次都不同没有固定随机种子脚本开头加rng(任意整数)轮廓系数偏低0.4K值不合适或特征区分度不够调整K值或补充峰谷差、负荷率等特征Kmeans出现空簇初始中心离样本太远用样本点初始化粒子位置别用均匀随机数运行极慢样本量大、粒子数多先抽样优化再全量聚类出现中文乱码或编码问题Matlab文件编码不符将m文件另存为UTF-8编码6.7 环境与版本小坑代码中用到的主要函数包括kmeans、silhouette、pdist2、zscore这些在Matlab基础平台和统计与机器学习工具箱里都有不需要安装额外的第三方包。如果你的Matlab版本比较老个别函数可能存在兼容性问题。比如silhouette函数目前在基础平台里已可用如果搜索不到检查是否安装了Statistics and Machine Learning Toolbox或者直接提示路径错误那就需要重新安装对应工具箱组件。另外遇到许可证或者激活报错比如license manager error -8之类的提示多数情况是旧版本许可服务残留和新版本激活信息冲突先彻底卸载旧的许可后台服务再重新激活这类问题直接查官方支持文档处理即可。从我个人的实操体会来说这个项目最有价值的不是PSO本身而是“用优化算法去解决另一个算法初始化痛点”的组合思路。我在后续做负荷预测和其他聚类项目时也沿用了“全局搜索做初始解经典算法做精细收敛”的模式效果一直很稳。如果你试着把PSO换成遗传算法或者灰狼算法只是适应度函数不变、更新公式换一下就能比较不同群智能算法的差异这个扩展思路也可以拿来继续玩。