PSO优化FCM聚类的居民用电行为分析Matlab实现

居民用电行为分析，这是近两年电力数据挖掘里非常常见的一个题目，也是很多刚入门数据分析方向的研究生和电力从业者绕不开的实战项目。我前阵子做的一个项目就是把粒子群算法优化FCM聚类这套组合用在居民的用电负荷曲线上，通过聚类把用户分成几类，再进一步做用电画像。这项目做完之后我回头整理了一遍代码，发现里面有不少坑和细节值得写出来，所以有了这篇文章。

这套技术路径说白了就一句话：FCM（模糊C均值聚类）负责把形状相似的负荷曲线归成一类，PSO（粒子群优化算法）负责帮FCM找到更好的初始聚类中心，避免FCM自己瞎猜初值然后在局部最优里出不来。结果就是聚类的稳定性明显提升，用户分组也更有业务可解释性。不管你是准备写期末作业、做毕业论文，还是接到类似的“用户分群”需求，这份思路和Matlab代码框架都能直接拿来改。

适合谁来读？首先是正在做FCM相关课题、但被“聚类结果每次跑出来都不一样”困扰的人；其次是想了解PSO怎么和传统聚类算法结合的人；最后就是纯粹想把居民用电数据变成业务结论的数据分析者。下面我把项目从思路到代码到调参踩坑完整过一遍。

1. 项目概述：为什么是“PSO+FCM”而不是直接上FCM

1.1 居民用电行为分析到底分析什么

很多人一听“居民用电行为分析”就想到了负荷预测，其实它们之间还有点距离。负荷预测是预测未来一时间段的用电曲线，而行为分析更偏“把人分群”：谁习惯晚上8点以后用电、谁白天耗电多、谁整天用得很平稳、谁专挑谷时电价时段开蓄热设备。这些信息拿到手之后，可以支撑阶梯电价套餐设计、需求响应潜力评估、台区增容规划等工作。

从数据层面看，分析对象通常是一个个用户在某段时间内的有功功率曲线，常见的是按15分钟或1小时一个采样点。拿1小时采样来说，一天就是24个点，这个用户就对应一个1×24的向量。把所有用户堆起来，就是一个n×24的矩阵。用电行为分析的本质，就是在这个矩阵上做无监督聚类，然后再对每一类中心曲线做解读。

1.2 FCM为什么适合负荷曲线聚类，又卡在哪

FCM在负荷曲线聚类里受欢迎，是因为它比KMeans更贴近真实世界的模糊性。一个用户的用电习惯不可能百分之百属于某个类型，FCM给每个样本算的是一组隶属度，比如“这个用户有70%像晚高峰型，20%像均衡型，10%像白天型”。这种软划分在连续变化的负荷数据上非常自然。

但FCM有个明显软肋：它的求解过程通常从一个随机的聚类中心出发，通过迭代更新隶属度和中心，目标函数是不断下降的。可这个函数本身是非凸的，初始中心选得不好，很容易收敛到局部最优解。表现出来就是你同同样的数据跑好几遍FCM，得到的聚类中心明显不一样，某几类用户群可能被合并掉，或者出现一个异常大的“垃圾桶类”。

1.3 PSO进来之后补上了什么短板

PSO做的事情很直接：它不在随机中心上硬迭代，而是把聚类中心当成一个“粒子”去种群化搜索。每一次迭代都有一群候选聚类中心方案在解空间里飞，所以最终选出来的那组初始中心基本能绕开局部最优，再丢给FCM做精细迭代。

用PSO的优点还有两个。第一是不需要对FCM函数求导，适合这类非光滑的优化问题；第二是参数可解释性强，惯性权重、学习因子都有明确意义，调起参来心里有数。缺点当然也有：增加了不小的计算开销，需要反复评估适应度函数。但如果样本量在几千到一两万这个量级，Matlab跑起来完全能接受。

2. 算法原理拆解：FCM的数学细节和PSO怎么插进来

2.1 模糊C均值的数学模型

FCM的目标函数长这样：

J_m = Σ_{i=1}^{n} Σ_{j=1}^{c} u_{ij}^m · ||x_i - v_j||²

其中x_i是第i个用户的负荷向量（24维），v_j是第j类的聚类中心（也是24维），u_{ij}是第i个用户属于第j类的隶属度，m是模糊指数（大于1，通常取2），c是聚类数。

约束条件是每个用户对各类隶属度之和为1：Σ_{j=1}^{c} u_{ij}=1，且每个u_{ij}在0到1之间。

求解过程是一个交替迭代：

第一步，固定聚类中心V，更新隶属度U：

u_{ij} = 1 / Σ_{k=1}^{c} (||x_i - v_j|| / ||x_i - v_k||)^{2/(m-1)}

第二步，固定隶属度U，更新中心V：

v_j = Σ_{i=1}^{n} u_{ij}^m · x_i / Σ_{i=1}^{n} u_{ij}^m

这两个公式实现起来都不难，但需要小心的是距离计算里如果出现d=0的项，会出现分母为0或者除0警告。工程实现里一般会做一个极小值保护。我在代码里统一用max(dist, 1e-10)这样的处理，后面会贴在代码区。

2.2 粒子群算法的更新机制

PSO模拟的是鸟群觅食行为。每个粒子代表解空间里的一个候选解，粒子有位置和速度两个属性。每次迭代里，粒子向自己的历史最优位置（pbest）和群体的历史最优位置（gbest）飞行：

v_i(t+1) = w·v_i(t) + c1·r1·(pbest_i - x_i(t)) + c2·r2·(gbest - x_i(t))

x_i(t+1) = x_i(t) + v_i(t+1)

公式里最值得解释的是w（惯性权重）。w大，粒子飞得猛，搜索范围广，适合前期探索；w小，粒子精细收敛，适合后期局部开发。我习惯用线性递减：迭代刚开始w=0.9，结束时w=0.4。c1和c2分别是“向自己学”和“向群体学”的系数，通常都取1.5附近。r1和r2是0到1的均匀随机数，用来保持随机性。

2.3 PSO和FCM怎么耦合：编码是核心

耦合方式在文献里有好几种，我这套采用的是“粒子位置 = 一组完整的聚类中心矩阵”。假设聚类数是c、每个样本是24维，那每个粒子就是一个长度c×24的向量。把这段向量按行拆成c×24的矩阵，就是FCM的一组初始中心V。

具体流程是：先将粒子位置转换成聚类中心，然后用这个中心去跑一轮FCM的核心迭代（或者不迭代直接算目标函数），得到J_m值作为适应度。PSO负责最小化这个J_m。等PSO收敛，拿gbest对应的一组中心作为FCM的初值，再跑标准FCM直到收敛。

这里有一个取舍需要说清楚：我的适应度函数是“用当前粒子中心跑一次FCM目标函数计算”，而不是把整个FCM迭代放进PSO循环里。如果每次评估都让FCM全部收敛完，PSO每一轮都会非常慢。实际应用里，先让PSO在粗粒度上找好中心，再用FCM精修，效果和耗时都更平衡。这也是很多论文里的做法，不是我自己拍脑袋想的。

2.4 目标函数与聚类评价指标的选择

既然要优化，就要接受“FCM的目标函数J_m不一定完全等价于‘聚类效果好’”。J_m只能反映组内紧凑程度，无法反映类间分离程度。所以我在项目中不只是看J_m，还会算两类评价指标来验证分群质量：

基础版项目只看轮廓系数就够用，如果写论文或答辩，Xie-Beni和DBI最好都拉出来做个表，说服力完全不同。PSO自身的适应度我仍然用J_m，因为它和PSO的搜索方向天然匹配，而评价指标只用来做横向验证。

3. Matlab实现要点：从数据预处理到核心代码

3.1 数据预处理与特征工程

聚类算法对数据尺度极其敏感，这一块别偷懒。直接把原始千瓦数喂进去，负荷高的用户就会主导距离计算，聚类结果基本被用电量大小带跑了，而不是被“用电形状”带跑。

我做这个项目时先做三步处理：

第一步，缺失值填补。对负荷曲线里个别空点，用前后两点的均值补上；连续缺失超过6个点（对应6小时）的用户直接剔除。

第二步，异常值过滤。用电曲线出现负值、单点突变超过正常值3倍以上，要么按前后均值修正，要么删除整个样本，视异常比例决定。这里的“3倍”我是按历史数据的3σ窗口来判定的，不是拍脑袋。

第三步，归一化。这一步决定聚的是“用量”还是“形状”。如果关注用电行为形态，就对每个用户除以其当日平均负荷，让每个用户曲线的总水平都在同一量级。这样聚类出来的是“峰谷形态”不同的人群，而不是“大户”和“小户”的区别。如果业务上需要兼顾用量，也可以选择均值-方差归一化（zscore），让X每一列服从均值0方差1。

处理完的X长这样：每行是用户在一天24小时的相对负荷，值围绕1上下波动。后面聚类的24个特征就是这24个点。你还可以再手工抽取特征，比如峰时电量占比、谷时电量占比、峰谷差、负荷率、早高峰时刻、晚高峰时刻，用这些压缩特征跑聚类，解释性更强。我的主程序写的是24维原始曲线版，特征抽取版作为对照实验。

3.2 PSO参数编码与主程序框架

下面给出的是我在项目里抽取并简化后的主程序，去掉了业务相关注释，保留了核心实现。我把它命名为pso_fcm_demo.m。

%% PSO优化FCM聚类 主程序 clear; clc; close all; %% 加载数据: X为n×24矩阵，n个用户，24个时刻负荷(已归一化) load('user_load_normalized.mat'); %% 参数设置 num_particles = 30; % 粒子数 max_iter = 100; % 最大迭代次数 c = 4; % 聚类数 m = 2; % FCM模糊指数 c1 = 1.5; % 个体学习因子 c2 = 1.5; % 群体学习因子 w_start = 0.9; % 初始惯性权重 w_end = 0.4; % 结束惯性权重 [n_users, d] = size(X); dim = c * d; lb = repmat(min(X), c, 1); ub = repmat(max(X), c, 1); lb = lb(:)'; ub = ub(:)'; %% 粒子群初始化 pos = rand(num_particles, dim) .* (ub - lb) + lb; vel = zeros(num_particles, dim); pbest_pos = pos; pbest_obj = zeros(num_particles, 1); for i = 1:num_particles pbest_obj(i) = fitness_fcm(pos(i,:), X, c, m); end [gbest_obj, gidx] = min(pbest_obj); gbest_pos = pbest_pos(gidx, :); %% PSO主循环 for iter = 1:max_iter w = w_start - (w_start - w_end) * iter / max_iter; for i = 1:num_particles vel(i,:) = w * vel(i,:) ... + c1 * rand(1,dim) .* (pbest_pos(i,:) - pos(i,:)) ... + c2 * rand(1,dim) .* (gbest_pos - pos(i,:)); pos(i,:) = pos(i,:) + vel(i,:); pos(i,:) = max(pos(i,:), lb); pos(i,:) = min(pos(i,:), ub); obj = fitness_fcm(pos(i,:), X, c, m); if obj < pbest_obj(i) pbest_obj(i) = obj; pbest_pos(i,:) = pos(i,:); if obj < gbest_obj gbest_obj = obj; gbest_pos = pos(i,:); end end end fprintf('iter = %d, gbest = %.4f\n', iter, gbest_obj); end %% 用PSO找到的最优中心初始化FCM，再做精调 V_init = reshape(gbest_pos, c, d); [U_final, V_final, Jm_final] = my_fcm(X, V_init, c, m); save('pso_fcm_result.mat', 'U_final', 'V_final', 'Jm_final');

这个主程序里有几个点我需要特意说明。粒子数30不算多，如果用户量两三万、特征24维，粒子数50到60会更稳，但时间也会成倍增长。迭代次数100次是经验值，多数数据集50轮左右就稳了。如果你发现gbest一直下降，说明还没收敛，可以加大迭代；如果早早就平了，说明100次有点浪费。

有人会问为什么不直接用全局最优中心就结束，还要再跑一次my_fcm。原因是PSO把每个粒子当“聚类中心候选解”来评估，但候选解只是中心，还没把隶属度迭代到最终状态。PSO搜索到的中心大概率在最优解附近，再把它们当FCM的初始化，用FCM的标准迭代去精修，这一步能再把目标函数往下压一点，让隶属度分布更干净。

3.3 适应度函数细节

适应度函数看起来就是FCM的目标函数，但实现有细节。我的做法是：给定一组聚类中心V，先算距离矩阵，由距离算隶属度U，再由U和距离算J_m。完整的fitness_fcm.m如下：

function Jm = fitness_fcm(particle, X, c, m) % particle: 1×(c*d)，是聚类中心的向量化 [n, d] = size(X); V = reshape(particle, c, d); % 距离矩阵: n×c，dist(i,j)是第i个样本到第j个中心的欧氏距离平方 dist2 = zeros(n, c); for j = 1:c diff = X - V(j,:); dist2(:,j) = sum(diff.^2, 2); end dist2(dist2 < 1e-10) = 1e-10; % 防除0 % 根据FCM公式计算隶属度 inv_dist = dist2 .^ (-1/(m-1)); U = inv_dist ./ sum(inv_dist, 2); % FCM目标函数 Jm = sum(sum((U.^m) .* dist2)); end

这段代码里最值得注意的就是dist2 < 1e-10的保护。真实数据有时候两个中心会重叠，或者某个用户恰好像某一类中心，这时候距离为0，0^(-1)直接变成Inf，后面的归一化就会出现NaN。加了这行保护后，计算稳定很多。

当然这个版本是纯公式实现，好处是任何工具箱都能跑。如果你装了统计和机器学习工具箱，可以用pdist2提速，尤其在大样本下能快不少。不过对教学和演示，循环版本更容易理解公式的含义。

配套的my_fcm.m也是手写的标准FCM迭代：

function [U, V, Jm] = my_fcm(X, V_init, c, m) V = V_init; for t = 1:200 % 计算距离矩阵 n = size(X,1); dist2 = zeros(n, c); for j = 1:c diff = X - V(j,:); dist2(:,j) = sum(diff.^2, 2); end dist2(dist2 < 1e-10) = 1e-10; inv_dist = dist2 .^ (-1/(m-1)); U = inv_dist ./ sum(inv_dist, 2); V_new = (U.^m)' * X ./ sum(U.^m, 1)'; if norm(V_new - V, 'fro') < 1e-5 V = V_new; break; end V = V_new; end dist2 = zeros(size(X,1), c); for j = 1:c diff = X - V(j,:); dist2(:,j) = sum(diff.^2, 2); end dist2(dist2 < 1e-10) = 1e-10; Jm = sum(sum((U.^m) .* dist2)); end

有人可能直接用Matlab模糊逻辑工具箱里的fcm()函数。那个函数自带初始化，但没办法方便地把PSO找到的中心作为初值灌进去。手写这个版本的目的是完全掌控迭代过程，方便复现论文结果。我建议你自己跑通这个手写版一次，能极大加深对FCM公式的理解。

3.4 结果可视化

聚类结果光看几个矩阵不够直观。我通常会画三张图：

第一张图，聚类中心曲线。横轴是0点到23点，纵轴是归一化负荷，把c条中心曲线画在一个坐标系里，一眼就能看出哪类用户是晚高峰型、哪类是均衡型。

第二张图，轮廓系数图。每个用户按所属类别排序，画轮廓系数竖条图，能直观看到哪类内部一致性差、哪些样本被分得不情愿。

第三张图，PSO收敛曲线。记录每一轮的gbest，画成下降曲线，用来判断参数是否合理。这张图也是写报告时很讨喜的素材。

代码大致这样：

figure; plot(0:23, V_final', 'LineWidth', 1.2); xlabel('时刻(h)'); ylabel('归一化负荷'); legend(arrayfun(@(k) sprintf('类型%d', k), 1:c, 'UniformOutput', false)); grid on;

4. 实验设计与聚类效果分析

4.1 模拟数据构造与场景设计

我自己手头并没有可以直接公开的用户负荷数据，项目里用的是一套由业务规律生成的仿真数据，模拟了四类典型用电用户：

生成方式就是在典型负荷曲线基础上加上高斯噪声，再加一点点随机峰时偏移，比如晚高峰型用户的峰值可能在19点到22点之间随机浮动。这样生成的模拟数据带有真实场景里的噪声和类内差异，用来验证算法能不能稳健分离四类，效果很直观。

模拟数据生成代码的核心其实就是对每个类型定义一条“基准曲线”，然后叠加噪声：

base_a = [0.35,0.33,0.32,0.31,0.33,0.40,0.55,0.65,0.55,0.42,0.40,0.45,0.50,0.48,0.42,0.38,0.40,0.60,0.85,0.95,0.80,0.60,0.45,0.38]; X1 = repmat(base_a, 80, 1) + 0.06*randn(80,24);

注意加噪声后用max(0.05, X1)把下限夹一下，防止出现负数。这种模拟数据很适合在代码阶段快速验证聚类效果，等验证完再换成真实台区数据就行。

4.2 聚类数怎么定才不会拍脑袋

聚类数c是最影响结果的决定性参数之一。我的做法不是凭感觉选4，而是先跑一遍不同c下的对比：对c从2到8，各跑10次PSO-FCM，取平均轮廓系数和DBI，然后画折线图。

选择标准通常有两个维度：

• 轮廓系数：越大越好，但要注意类数增加后轮廓系数通常会单调上升，不能只看峰值；
• DBI：越小越好，DBI的肘部位置往往对应合适的c；
• 业务约束：聚类结果最终要能解读成不同营销策略和电价套餐，类别太多后面没法落地。

我这个仿真数据的肘部出现在c=4或c=5附近，结合DBI数值，最终定4类。实际项目中，如果业主方明确说“我们套餐就三种”，也可以直接c=3，然后在解释上下工夫。

4.3 FCM与PSO-FCM的效果对比

同样的数据，我分别跑了标准FCM（随机初始化100次取最好）和PSO-FCM，对比两个维度的表现。

一是收敛性。标准FCM从不同随机初值出发，J_m最终在46.8到52.3之间波动，说明多次运行会得到质量不同的划分；PSO-FCM的gbest稳定落在45.9左右，最后再经FCM精调得到45.7，明显更优也基本固定。

二是类别可分性。标准FCM容易出现某两类中心靠得很近，甚至D类（谷电型）被整体并入A类的现象，导致业务上没法解释“为什么夜里用电的人没了”。PSO-FCM则能稳定还原四类，各类轮廓系数都在0.35以上，其中均衡型最高，因为这类用户靠类内一致性就很好区分。

从运行时间看，标准FCM100次随机初始化大约用时8秒，PSO-FCM整体用时约150秒（粒子数30，迭代100次）。时间多了近20倍，但换来的是稳定性和可解释性。如果跑十几次结果都不稳定，后面所有业务分析都是空中楼阁，我个人认为这个时间成本值得。

4.4 聚类结果如何翻译成业务语言

聚类完成之后，真正有价值的是给每一类用户打标签、写画像。我按PSO-FCM得到的四类中心曲线，在项目报告里是这样表述的：

• A类“晚高峰型”：中心曲线在19点到22点出现明显高峰，负荷率低，峰谷差大。适合参与晚峰需求响应，或推送“夜间谷电优惠套餐”。
• B类“白天型”：高峰集中在9点到17点，多半是白天在家人群，错峰充电桩对这个群体意义大。
• C类“均衡型”：全天曲线平缓，峰谷差小，用电调节弹性小，适合稳定的阶梯电价。
• D类“谷电型”：夜间1点到5点用电显著高于白天，说明可能有蓄热或家电定时启动，是谷时移动电量的优质潜力用户。

每一类画像里还可以附上统计指标，比如类内用户的平均日用电量、峰谷差、负荷率、用户占比。这些数字放到PPT里，业务方一眼就懂了。

5. 常见问题与排查记录

5.1 适应度不收敛或来回震荡

如果你发现PSO迭代曲线不降反升，或者上下剧烈震荡，先别改算法，优先检查三个地方。

第一个是边界越界处理。粒子位置超出[lb, ub]后，如果直接不处理，聚类中心就可能跑到数据范围之外，距离矩阵变大，J_m反而上升。上面主程序里的max/min夹取操作必须保留。

第二个是惯性权重设置。w固定为0.9不下降，粒子飞得太随意，收敛会很慢；w一开始就0.2，粒子容易陷入局部。线性递减是我最常用的策略，简单有效。

第三个是随机种子。Matlab不设固定种子时，每次运行结果都不同，少量实验可能让你误以为算法震荡。排查阶段先rng(42)固定随机种子，等调好参数再放开对比多次结果。

5.2 每次运行聚类结果不一样

这个问题的根源有两个：FCM自身初值敏感，以及PSO是随机搜索算法。只要不是固定种子，第二次运行gbest粒子大概率是另一个位置，聚类结果可能完全分得不一样。

我的应对策略是“多次运行+中心排序”。大数据量下把整个PSO-FCM跑10次，每次记录最终聚类中心和目标函数，保留目标函数最好的一次。同时聚类中心顺序在每次运行中可能互换，比如上次叫“类型1”的这次变成了“类型2”，单纯看类型编号没有意义。处理方法是对最终聚类中心按24点曲线的某类特征排序，比如按峰值时刻排序，峰越晚的排前面，这样跨次对比就稳定了。

5.3 代码报错与工具箱依赖问题

Matlab跑这段代码时最容易报错的点有这么几个：

其中reshape维度问题是最容易绊倒新手的。一个粒子长度是c×d，reshape成c行d列时，顺序必须和初始化时repmat的排列方式一致。我建议统一用“把c个聚类中心按行堆叠成一个向量”，初始化的lb和ub也按同样方式生成，这样从粒子到矩阵的转化就不会乱。

5.4 聚类结果无法解释

跑出来的聚类中心曲线形状差不多，或者有一类样本数特别少，这类问题通常不是算法bug，而是数据或者预处理导致。

最常见的原因是没做归一化，聚类被“用电量大的用户”主导，结果类型全是“大户1、大户2、小户”，没有峰谷形态差异。解决办法就是按日平均功率归一化，重新跑。

还有一种情况是数据里确实混入了异常用户，比如整日负荷为0或负荷恒定不变。这类样本放进去会把某个类拉成一个“平直线类”，业务价值很低。我在预处理时会把日方差小于阈值的样本单独剔除或单独归为特殊类。判断标准可以是我用的std(X,0,2) < 0.05，这类用户不是没上报数据就是有长期空置。

5.5 Matlab版本兼容问题

老版本Matlab对arrayfun匿名函数、pdist2这些支持程度不一。如果你用的是比较老的版本，建议把arrayfun调用改成最朴素的for循环，代码性能差别不大，但兼容性会好很多。我发布的版本尽量不用高版本独有的函数，保证大多数Matlab都能跑。

6. 我自己踩过的坑和扩展思路

6.1 参数选择心得

做这类项目，最怕一上来就想着“调一个最优参数”。我的实际体会是：先固定一批合理默认值，把流程跑通，再逐个调参。

先调的永远是聚类数c，因为它对业务影响最大；其次调模糊指数m。m取2是主流选择，m越大聚类边界越模糊，小样本下不建议超过2.5；最后调PSO参数。PSO参数其实很皮实，w在0.4到0.9线性衰减、c1=c2=1.5，大多数场景都够用，刻意微调收益不大。

粒子数倒是值得根据数据规模调整。数据量小（几百个样本），粒子数20就够；样本几千到几万，粒子数建议50到80；如果超过十万，单机Matlab跑PSO会比较吃力，建议先对数据抽样或降维。我试过在四万用户24维数据上跑粒子数100迭代200次，大概是十几分钟量级，还能接受，但再大就得考虑并行或换工具了。

6.2 后续可以往哪些方向扩展

这个项目做完，扩展空间其实很大，每一个方向都足够单独开一个新项目。

可以把24维负荷曲线换成更高频的数据（15分钟96点）甚至多维特征，比如温度、节假日类型、用电类型标签，做联合聚类；可以把PSO换成其他元启发式算法，比如灰狼优化、鲸鱼算法，和PSO做对比实验，这类题材在论文里非常吃香；也可以在聚类结果基础上继续做负荷预测，每一类用户单独训练一个预测模型，通常比全局模型精度高不少。

如果往工程化方向走，还可以把这段Matlab代码封装成函数，用Matlab Compiler做成独立程序，或者翻译成Python版本部署到数据挖掘平台上。我在实际项目里就把这套逻辑翻译成了Python的numpy版本，核心部分也就一百多行，维护起来反而更轻松。

最后分享一个我每次做聚类项目都会保留的小操作：把最终聚类结果里每一类的典型用户曲线和该类中心曲线放在同一张图里，画成半透明细线加粗中心线的对比图。这张图能非常直观地展示“类内差异大不大”“有没有需要拆分的次级模式”，比任何指标都更能说服业务方。很多时候算法上看着没问题，业务上一看图就发现某一类内部其实还能再分，这种细节才是项目落地的关键。你可以拿到自己的数据上试试，一定会对聚类结果有新的认识。