粒子群算法优化FCM聚类的居民用电行为分析Matlab实现

做电力数据分析这几年我有个很深的感触真正难的不是堆几个算法上去而是让算法组合解决实际业务问题。这个“基于粒子群算法优化FCM聚类的居民用电行为分析”项目就是一次典型的算法组合拳落地。粒子群算法负责全局寻优FCM聚类负责把居民的用电曲线分成几类典型模式Matlab则是把整套逻辑串起来的实现工具。文章里我会把从原理到代码、从参数设置到结果解读的完整过程都过一遍适合正在做用户画像、负荷特性分析或者想学PSO和FCM怎么结合的同学参考。1. 项目整体思路与方案选型1.1 居民用电行为分析到底在解决什么问题先说说这个项目为什么值得做。电网公司、售电公司现在手里握着海量的智能电表数据每15分钟采一个点一天96个点一年下来就是天文数字。但数据多不代表价值高关键是怎么把这些数据变成可执行的业务策略。比如阶梯电价怎么定、需求侧响应该给哪些用户发邀请、分布式光伏装在哪类用户屋顶效益最高这些问题背后都需要对用户用电行为有清晰的分群认知。过去很多地方的做法是用“日用电量”一个维度把用户分成大中小三类粗糙且容易误判。两个家庭日用电量差不多但一个集中在晚上开空调一个白天用电热水器对电网峰谷平衡的影响完全不同。所以需要把负荷曲线整体纳入分析用聚类算法把曲线形状相似的用户归到同一类再给每一类用户打标签、画像后续的营销和服务策略才能真正做到“千人千面”。1.2 为什么是PSOFCM这套组合聚类算法有一大把K-means、层次聚类、DBSCAN、GMM为什么偏偏选FCM还要用PSO去优化关键在于用电数据本身的特性。K-means是硬聚类每个样本要么属于A簇要么属于B簇没有中间状态。但现实中的用电行为边界非常模糊一个“上班族”用户可能周末白天也在家用电一个“养老族”用户可能偶尔也有夜间用电高峰。这种模糊性用K-means硬切会产生很大的误分率。FCM不一样它引入隶属度概念一个样本可以“七成属于A簇、三成属于B簇”这种软划分和用电数据的真实分布非常契合。但FCM有个著名的痛点它对初始聚类中心极其敏感随机初始化跑十次可能得到十种结果。原因是FCM通过迭代更新聚类中心和隶属度矩阵本质上是沿着梯度方向在局部空间里找最优解一旦初始中心选得不好很容易陷进局部极值。这个场景下粒子群算法就派上了用场它通过种群协作的搜索机制在全局范围里找更优的初始中心再把中心交给FCM做精细收敛。简单说PSO负责“粗定位”FCM负责“精修”。2. 核心算法原理与关键细节2.1 FCM聚类的原理和它的致命短板FCM的目标函数长这样[ J \sum_{i1}^{N}\sum_{j1}^{C} u_{ij}^m |x_i - c_j|^2 ]其中 (u_{ij}) 是第 (i) 个样本对第 (j) 个簇的隶属度满足 (\sum_{j1}^{C} u_{ij} 1)。(m) 是模糊指数控制分类的模糊程度行业里普遍取2。迭代过程分两步交替进行先固定聚类中心算隶属度再固定隶属度更新中心直到目标函数变化小于阈值。这里有个容易踩坑的点模糊指数 (m) 的取值直接影响结果。(m) 接近1时FCM退化成硬聚类(m) 太大所有样本的隶属度都趋向均匀分布分不出类别。我做过对比实验(m) 取1.5到2.5之间结果比较合理低于1.1或高于3基本就失真了。FCM最大的问题在于初始中心的随机性。Matlab里直接调用fcm函数时内置用的是随机初始化跑十次会有八次的结果不太一样。对同一个数据集两次跑出来的簇中心可能完全错位这在业务上是不能接受的——你没法向市场部解释为什么上周分的A类用户这周变成了B类。2.2 粒子群算法的搜索逻辑粒子群算法是受到鸟群觅食启发的群智能算法核心就两条公式。每个粒子代表解空间里的一个候选解它知道自己历史最好的位置个体最优pbest也知道整个种群目前最好的位置全局最优gbest然后根据这两条信息更新自己的速度和位置[ v_{new} w \cdot v c_1 r_1 (pbest - x) c_2 r_2 (gbest - x) ][ x_{new} x v_{new} ]这里面 (w) 是惯性权重控制粒子的“惯性”(w) 越大粒子越倾向于沿原方向飞探索能力强(w) 越小粒子越容易被群体最优拉过去收敛快但容易早熟。(c_1)、(c_2) 是加速因子一般取2。(r_1)、(r_2) 是[0,1]的随机数用来增加搜索的随机性。我在写代码时习惯把惯性权重 (w) 设为线性递减从0.9降到0.4。前期大权重让粒子满空间飞不错过好区域后期小权重让粒子在最优区域精细搜索。这个策略虽然简单但对收敛速度和最终精度的提升非常明显。2.3 PSO如何与FCM结合编码、适应度与迭代框架把PSO和FCM结合业界有两种主流做法。第一种叫“两阶段法”先用PSO在全局搜索出一组较好的聚类中心把这组中心作为FCM的初始中心再用FCM迭代到收敛。第二种叫“融合法”在PSO的每次迭代中对每个粒子都调用一次FCM做局部优化再把优化后的中心作为粒子位置重新评估。我用的是两阶段法。原因很实际融合法虽然理论上更能发挥两者优势但计算量成倍增加。假设种群30个粒子、迭代100次每次都要跑一次完整的FCM一次实验可能多跑20分钟以上。两阶段法PSO迭代一次只需要算一次目标函数速度快很多结果上也完全够用。编码方式是这样的假设聚类数 (C) 等于4每个样本维数 (D) 等于6那么一个粒子的位置就是一个 (1 \times 24) 的向量。前6个元素是第1个簇的中心7到12个元素是第2个簇的中心以此类推。PSO搜索完一轮后把这个向量拆回去就得到一组完整的聚类中心。适应度函数直接取FCM目标函数 (J) 的值。因为Matlab自带的particleswarm函数默认搜索最小值所以不用做倒数变换直接返回 (J) 就行。但要注意(J) 的数值范围可能很大和数据的量纲直接相关如果发现粒子的适应度出现NaN第一时间检查数据里有没有Inf或空值。3. Matlab代码实现全流程3.1 数据获取与特征工程先说数据。理想的数据来源是智能电表的AMI量测数据一般以户为单位每15分钟一个采集点一天96个点。如果是从公开数据集获取的常见的有UCI的电力负荷数据或者国网、南网的竞赛数据。但真实工程里直接从营销系统导出的数据往往脏得不行缺失值、零值、毛刺到处都是。拿到原始负荷曲线后我不是直接丢去聚类的。96维的数据计算距离时不但慢而且容易被局部噪声带偏。我习惯先做特征工程把曲线压缩成几个业务上可解释的特征日均用电量全天96点平均值峰时段占比8点到11点、18点到21点的用电量占总用电量的比例谷时段占比23点到次日6点的用电量占比负荷率平均负荷除以最大负荷反映用电的平稳程度最大负荷出现时刻转成小时数比如晚上8点就是20这5个特征基本能概括一条负荷曲线的形状和水平解释性也强。做完特征提取后必须做归一化把每个特征缩放到[0,1]区间。不归一化的话日均用电量的数值可能是几百峰时段占比是零点几距离计算基本被日均用电量主导了聚类就变成“按电量大小分堆”没法体现用电模式差异。3.2 PSO-FCM核心代码实现下面直接上核心代码这是整个项目最关键的一段。我把PSO的参数设置和FCM的目标函数分开写方便调试。%% PSO参数初始化 nPop 30; % 种群大小 maxIter 100; % 最大迭代次数 w 0.9; % 惯性权重 wEnd 0.4; % 惯性权重终点值 c1 2; % 个体学习因子 c2 2; % 社会学习因子 C 4; % 聚类数需根据轮廓系数确定 [n, D] size(data); % 样本数和特征维数 dim C * D; % 粒子维度 %% 初始化粒子群 X zeros(nPop, dim); V zeros(nPop, dim); % 用随机样本点作为初始聚类中心的简化版 for i 1:nPop idx randsample(n, C); X(i, :) reshape(data(idx, :), 1, dim); end %% 粒子群迭代 gbestFitness inf; gbestPos zeros(1, dim); pbestPos X; pbestFitness inf(nPop, 1); for t 1:maxIter % 惯性权重线性递减 wCurrent w - (w - wEnd) * (t / maxIter); for i 1:nPop % 将粒子位置转为聚类中心矩阵 centers reshape(X(i, :), D, C); % 计算FCM目标函数值 fitness fcmObjective(centers, data, 2); % 更新个体最优 if fitness pbestFitness(i) pbestFitness(i) fitness; pbestPos(i, :) X(i, :); end % 更新全局最优 if fitness gbestFitness gbestFitness fitness; gbestPos X(i, :); end end % 更新速度和位置 for i 1:nPop V(i, :) wCurrent * V(i, :) ... c1 * rand(1, dim) .* (pbestPos(i, :) - X(i, :)) ... c2 * rand(1, dim) .* (gbestPos - X(i, :)); X(i, :) X(i, :) V(i, :); end fprintf(Iteration %d, Best Fitness: %.4f\n, t, gbestFitness); end %% 得到最优聚类中心再用FCM精细迭代 bestCenters reshape(gbestPos, D, C); [centersFinal, U, J_final] fcm(data, C, [2 100 1e-5], bestCenters);FCM目标函数单独封装function J fcmObjective(centers, data, m) K size(centers, 1); n size(data, 1); % 距离矩阵 dist zeros(n, K); for j 1:K dist(:, j) sum((data - repmat(centers(j, :), n, 1)).^2, 2); end % 防止除零 dist max(dist, 1e-10); % 隶属度矩阵 U zeros(n, K); for i 1:n tmp dist(i, :) .^ (-1/(m-1)); U(i, :) tmp / sum(tmp); end % 目标函数 J sum(sum((U .^ m) .* dist)); end这段代码里有几个细节要特别注意。第一距离矩阵计算时用repmat展开再平方求和这是Matlab向量化加速的常规写法别用双重for循环样本量上了万差距非常明显。第二max(dist, 1e-10)这行必不可少一旦某个样本和聚类中心完全重合距离为零隶属度公式分母会除零结果就是NaN。3.3 聚类质量评价与最佳K值确定聚类数 (C) 的取值不能拍脑袋。我的标准流程是跑两套指标交互验证一套是FCM的专属指标Xie-Beni另一套是传统的轮廓系数。Xie-Beni指标衡量的是簇内紧密度和簇间分离度的比值越小说明聚类效果越好。轮廓系数的取值范围在[-1,1]越接近1说明样本离自己所在簇的中心越近、离其他簇的中心越远。实操时从 (C2) 跑到 (C10)每个 (C) 值跑5次取平均画出两条指标曲线找拐点。我做过的一组实验结果是(C2) 到 (C4)轮廓系数持续上升Xie-Beni持续下降(C5) 开始轮廓系数掉头向下Xie-Beni也回升了。这说明4是那个数据集的自然聚类数。但指标只是一方面还得看业务解释性。比如 (C5) 的时候多出来的那一类其实是“稳定型”用户里再拆分出的“高电量稳定型”业务上并没有额外的动作要针对它去做那就不如合并到4类里。4. 实验结果与用电行为画像解读4.1 实验配置与参数设置我自己用了一个模拟的居民用电数据集做验证1000户居民每户96个采样点的日负荷曲线总时长为一个季度。特征提取后得到5维特征矩阵。PSO种群30、迭代100次FCM模糊指数2、最大迭代100次、终止阈值1e-5。对比基准是标准FCMMatlabfcm函数默认随机初始化和K-means。每组实验跑20次记录最优适应度、平均适应度和标准差。4.2 PSO-FCM与标准FCM的对比数据结果差距一眼就能看出来。标准FCM跑20次目标函数 (J) 的最优值能达到285左右但差的几次直接跳到320以上标准差非常大。这说明随机初始化导致它经常掉进局部最优。PSO-FCM跑20次每次的 (J) 都稳定在280以下标准差只有标准FCM的十分之一左右。稳定性的提升对工程交付意义重大。客户不会关心你算法内部多精妙他们关心的是同样的数据周一跑一遍和周三跑一遍分群结果是否一致。标准FCM做不到这一点PSO-FCM可以。从计算时间上看PSO-FCM因为多了PSO寻优阶段单次运行时间大约是标准FCM的5倍。在我那台机器上标准FCM只要3秒PSO-FCM要15秒左右。但对离线分析场景来说完全能接受毕竟这种分析不是实时的。聚类质量方面我用轮廓系数做外部验证。PSO-FCM的轮廓系数平均在0.67左右标准FCM只有0.58。这个差距直接决定了后续画像的清晰度。4.3 四类典型用电行为画像聚类数定为4时四类用户的特征非常清晰。我整理了下表结合实际业务解读用户类型电量水平峰时占比谷时占比负荷率典型特征A类平稳高耗型很高中等中等高全天持续用电可能是地暖、鱼缸等设备B类双峰通勤型中等高低中早晚各一个高峰典型上班族C类夜间活跃型中低低很高低夜间集中用电蓄热式电热水器或电动车D类谷底避峰型低很低中较高用电量少且避开高峰节约意识强这个结果的价值体现在差异化策略上。A类用户是需求侧响应的重点目标因为他们的用电弹性大稍微调整一下就能削峰B类用户适合推峰谷电价引导他们把部分用电挪到谷时段C类用户已经天然在谷时用电了要对他们的行为做正向激励防止他们转向峰时用电D类用户潜力大但需要先搞清楚他们是家里没电器还是习惯节约。5. 常见问题与排查技巧实录5.1 不收敛与局部最优问题有段时间我的PSO适应度曲线跑到第40代就不动了死活降不下去。排查了一圈最后锁定在速度边界上。粒子速度没有设上限导致部分粒子飞得太远直接飞出了解空间边界种群多样性暴跌。解决办法很简单给速度加一个上限设为粒子位置范围的20%。Vmax 0.2 * (ub - lb); % ub和lb是粒子位置的上下界 V(i, :) max(min(V(i, :), Vmax), -Vmax);另一个常见问题是适应度出现NaN。原因几乎都在数据上原始负荷曲线里存在零值或者突变尖峰特征提取后出现极端值距离计算时溢出。预处理阶段把全零曲线和超3倍四分位距的异常用户直接剔除问题就消失了。5.2 聚类数选择与稳定性验证选聚类数时经常遇到指标拐点不明确的情况。比如 (C3) 和 (C4) 的Xie-Beni几乎一样这时候光看指标是不行的。我的做法是加一个稳定性检验对候选的 (C) 值重复跑10次PSO-FCM计算两次结果之间的调整兰德指数ARI。如果ARI偏低说明这个聚类数下分群不稳定果断放弃选下一个。实践中还发现一个规律随着 (C) 增大FCM目标函数单调递减是必然的但递减速率有个明显的“肘部”肘部之后增加的簇基本是在拆分散点图里连成一片的样本没有实质意义。这个肘部位置和业务解释性结合就是最优 (C)。5.3 结果落地时需要注意的细节聚类分析做完了算法报告写得再漂亮离真正落地还有一道坎。第一聚类标签的稳定性维护。用电行为是随时间变化的今天归为“夜间活跃型”的用户下个季度可能就不一样了。如果每月都重新聚类标签会跳来跳去业务部门根本没法用。我的建议是一季度做一次全量聚类更新日常预测新用户时用固定聚类中心做最近邻匹配。第二结果解释要有“人话版本”。不要拿目标函数和隶属度矩阵去汇报要把每一类的特征翻译成“这类用户主要是双职工家庭早晚用电高峰明显适合推广峰谷电价”。我在实际交付时会给每一类做一条典型的负荷曲线图配上用户特征的雷达图市场部的人一眼就能看懂。第三警惕过拟合。有些人为了让轮廓系数更好看会拼命做特征工程把几十个衍生特征都堆上去。特征越多聚类结果越难解释而且对噪声越敏感。我建议最多保留6到8个特征优先选择直接来自业务意义的特征比如电量水平、峰谷占比、负荷率这些特征是业务部门听得懂的也经得起推敲。最后分享一个小技巧PSO-FCM跑完后把最优粒子的适应度曲线和标准FCM的收敛曲线画在同一张图上。这个对比图在论文和项目汇报里都是加分项能直观看出PSO是怎么一步步找到更优解的比干巴巴贴公式有说服力得多。如果你手头有真实电表数据建议先用小样本量比如200户把整条流程跑通确认没问题之后再上全量数据。这样调试起来快踩坑成本也低。