PSO优化FCM聚类：居民用电负荷分析原理与Matlab实现

直接把这段经历写出来，是因为我觉得很多做电力负荷分析、用户画像的同学，都在用FCM聚类但总被“初值敏感、容易陷局部最优”折磨。做居民用电行为分析，核心是把用户的负荷曲线分成几类：有人白天用电多，有人晚上用电多，有人周末飙升。传统FCM能分，但每次跑结果都不一样，聚类中心初始化一换，类别就飘。于是我把粒子群算法引入进来，用PSO先去搜索一套好的聚类中心，再交给FCM去细化，整套用Matlab实现，效果稳定不少。这篇就围绕这个方案，把原理、代码、调参踩坑全部讲透，给打算用智能优化算法做负荷聚类分析的朋友一个可复现的参考。

1. 为什么选粒子群优化FCM：从痛点说起

1.1 FCM的先天不足：初值敏感与局部最优

模糊C均值聚类（FCM）在负荷曲线分类里用得非常多，因为它不像K-means那样硬性把每个样本归到某一类，而是给每个样本一个隶属度，描述它“多大程度上属于这一类”。对于居民用电这种边界模糊、行为波动大的数据，FCM的软划分确实比硬聚类合理。但用过的人都知道，FCM的迭代过程是从初始聚类中心出发，不断更新隶属度矩阵和聚类中心，直到目标函数收敛。这个目标函数通常是类内加权平方误差和，它是个非凸函数，山谷很多。初始点落在哪个山谷附近，最后就收敛到哪个局部极小值。换个随机种子，结果就变。对于居民负荷曲线这种本身就高维、噪声多的数据，这个问题尤其明显。

我最早跑FCM时，同一份数据分两次，结果里“工作日白天型”和“夜间型”的比例能差出10多个百分点。后来尝试多次随机初始化取最优，虽然能缓解，但本质上还是在同一个非凸函数里碰运气。于是我把目光转向全局优化算法，想找一个能在解空间里做全局搜索的机制，先给FCM一个靠谱的起点。

1.2 PSO与FCM结合的逻辑

粒子群算法（PSO）模拟鸟群觅食，每个粒子代表解空间中的一个候选解，粒子跟着个体历史最优和群体历史最优飞行，速度和位置迭代更新。它不需要梯度信息，对目标函数形式没有要求，特别适合给FCM这种“内部迭代嵌套复杂”的问题做外层优化。

这里有个关键设计：到底用PSO直接做聚类，还是用PSO优化FCM的初始聚类中心？我试过前者，也就是把每个粒子的位置直接编码成全部聚类中心，用聚类评价指标当适应度，迭代到底。这样做的缺点是PSO在连续空间里搜索高维中心点，收敛慢，而且最后一步缺乏模糊隶属度的精细调整，轮廓系数往往不如“PSO粗搜+FCM细调”的方案。所以我最终采取的策略是：先让PSO在解空间里跑几十代，把每个粒子位置解码成一组聚类中心，用FCM的一轮迭代结果（或者直接用FCM目标函数值）作为适应度，让粒子群往全局较优区域飞；迭代结束后取出历史最优位置作为初始聚类中心，再交给标准FCM迭代到收敛。这样既利用了PSO的全局搜索能力，又保留了FCM的局部精炼能力，速度和精度兼顾。

1.3 适用范围与场景分析

这套方案主要是针对中小规模数据集，比如几千户居民的日负荷曲线、96点采样或者24点采样。样本数几千、维度几十到上百，PSO的粒子数和迭代次数控制在几十这个量级，Matlab跑起来一点压力没有。

如果你的数据是海量用户、分钟级采集，比如上百万户每15分钟一个点，那直接用PSO优化的矩阵运算会比较吃力。这时候更好的做法是抽样做聚类原型，再用MapReduce或者Mini-Batch思路扩展。另外，PSO-FCM并不是万能的，它解决的是初始值敏感问题，如果数据本身没有明显的簇结构，任何聚类算法都分不出有意义的结果。所以做这个分析前，我建议先算一下数据的聚类倾向性，比如Hopkins统计量，再决定要不要上这套组合拳。

2. 粒子群算法原理与FCM聚类的核心细节

2.1 粒子群算法：一群鸟是怎么找最优解的

PSO的原理说穿了不复杂。每个粒子有两个属性：位置和速度。位置代表一个候选解，速度代表下一步移动的方向和幅度。每次迭代，粒子比较自己历史最好的位置（pbest）和整个群体历史最好的位置（gbest），然后按照下面公式更新速度：

v(i+1) = w * v(i) + c1 * r1 * (pbest - x(i)) + c2 * r2 * (gbest - x(i))

x(i+1) = x(i) + v(i+1)

这里w是惯性权重，控制粒子的“惯性”，w大全局搜索强，w小局部开发强。c1、c2是学习因子，分别代表向个人经验和群体经验学习的程度。r1、r2是0到1之间的随机数，保证搜索的随机性。

我用个生活化的类比，粒子就像一群在山上找最低点的游客，每个人记得自己走过的最低位置，还通过喊话知道团队里谁找到了最低位置。下一步往哪走，一部分取决于自己惯性，一部分朝着自己最低点方向，一部分朝着团队最低点方向。这里的“最低点”就是目标函数最小化，对FCM来说就是类内距离加权和最小。

实际实现中有几个细节必须注意。一是速度钳制，如果不限制v的最大值，粒子可能飞出去，解直接爆掉。二是边界处理，位置超出解空间范围时要拉回来，常见做法有边界吸收和边界反弹。三是惯性权重衰减，我习惯从0.9线性降到0.4，前期大力探索，后期精细搜索。这样一套标准PSO下来，稳定性比固定权重好很多。

2.2 FCM聚类：模糊归属怎么算

FCM的核心是把每个样本到各个聚类中心的距离通过隶属度加权。目标函数是：

J = sum(i=1 to N) sum(j=1 to C) u(i,j)^m * ||x(i) - c(j)||^2

其中u(i,j)是样本i对第j类的隶属度，m是模糊指数，一般取2。隶属度满足sum_{j=1}^{C} u(i,j)=1的约束。求解时通过拉格朗日乘子法，得到隶属度更新公式和聚类中心更新公式：

u(i,j) = 1 / sum_{k=1}^{C} (||x(i)-c(j)|| / ||x(i)-c(k)||)^(2/(m-1))

c(j) = sum_{i=1}^{N} u(i,j)^m * x(i) / sum_{i=1}^{N} u(i,j)^m

理解这两个公式的重点在于：隶属度不是凭空产生的，它是基于样本到各个中心距离的相对比值算出来的。如果样本距第j类中心很近，距其他类中心很远，u(i,j)就接近1；如果距离差不多，隶属度就被摊薄。这个“软”特性对居民用电行为特别合适，比如一个用户工作日白天用电多、周末晚上用电多，硬聚类要么归这类要么归那类，FCM则会给出“70%属于上班族型，30%属于周末宅家型”这种更真实的结论。

2.3 PSO-FCM的目标函数与编码方式

把PSO和FCM结合，首先要定义粒子的编码。假设要聚成C类，负荷曲线归一化后是D维，那么每个粒子的位置就是长度为C*D的一维向量，依次存放第1类的D维中心、第2类的D维中心……第C类的D维中心。解码时reshape成C行D列的矩阵就是一组聚类中心。

适应度函数怎么定？我试过两种。第一种是直接算FCM的目标函数值J，用PSO最小化它。第二种是先用这组中心做一轮FCM迭代，然后把迭代后的J作为适应度。第二种比第一种稳定，因为FCM一轮迭代会把中心拉向一个局部更优的位置，适应度地形更平滑，PSO不容易被噪声带偏。但代价是每评估一次粒子都要跑一轮FCM迭代，计算量多一个数量级。对于几千样本、几十维的负荷数据，这个代价完全可接受。我最终的方案就是第二种：每个粒子先映射成聚类中心，然后执行一次FCM核心更新（即更新隶属度矩阵、再更新中心），返回更新后的目标函数值。

这里有个容易搞错的点：PSO搜索到的gbest位置是“一次FCM迭代后的中心”，不是原始粒子位置。所以在PSO结束后，不能直接拿gbest当最终中心，应该把它重新解码，再执行标准FCM直到迭代次数耗尽或目标函数变化小于阈值。这样才能发挥FCM的收敛精炼能力。

3. Matlab代码实现全流程拆解

3.1 数据准备与预处理

我用的居民用电数据是某地区智能电表采集的日负荷曲线，每天96个采样点，也就是每15分钟一个点。原始数据从CSV读入，每一行是一户一天的数据，行首是用户ID，后面96列是功率值。直接拿原始功率做聚类会有问题，因为不同用户的用电量基数差别很大，有人一天10度电，有人50度，曲线形状相似但幅值不同，聚出来的结果容易按“用电量大小”而不是“用电行为模式”分群。

所以必须做归一化。我选了最大最小归一化，把每户的负荷曲线缩放到0到1之间：

normalized = (x - min(x)) / (max(x) - min(x))

这样处理后，曲线只保留形状特征，幅值信息被抹掉。如果你也想看用电量级别，可以另外把日均用电量作为一个附加特征拼进去。但要注意附加特征的量纲影响，拼之前最好也归一化。我做的这个版本主要关注行为模式，所以只用归一化后的96维负荷形状。

另外，数据清洗也不能省。智能电表数据经常有缺失值、零值连续段、异常跳变。我的处理规矩是：如果某天缺失率超20%，直接删掉；连续超过4个采样点为零且整日用电量非常低，多半是户主外出或者空置房，这类数据如果不删，聚类时会单独成一类“零用电型”，反而干扰正常居民行为分型。保留它其实也可以，但分析时会发现这类用户没有太多需求响应价值，建议第一版先剔除。

3.2 主函数框架与参数设置

Matlab代码整体框架分四段：加载数据、PSO优化、FCM精炼、结果分析。我先把参数设置列出来，这套参数经过多轮试验，针对96维负荷曲线、聚类数C=5、样本数2000左右的情况表现稳定。

clc; clear; close all; % 数据加载与预处理 load trainData.mat; % data: N x 96, 已经归一化 [N, D] = size(data); % 聚类参数 C = 5; % 聚类数，需要根据实际业务确定 m = 2; % 模糊指数，FCM标准取2 maxIter = 100; % FCM最大迭代次数 tol = 1e-5; % FCM收敛阈值 % PSO参数 numParticle = 40; % 粒子数 maxGen = 80; % 迭代代数 wMax = 0.9; % 最大惯性权重 wMin = 0.4; % 最小惯性权重 c1 = 1.5; % 个体学习因子 c2 = 1.5; % 群体学习因子 vMax = 0.1; % 粒子最大速度（按位置范围比例定）

这里重点说一下聚类数C的确定。多数人是拍脑袋定的，我建议跑一遍不同C值的轮廓系数，选轮廓系数大且业务可解释的C。比如C=3、4、5都算一下，C=5时轮廓系数不是最高，但5类行为能对应上“上班族”“晚间活动型”“全天空调型”“周末型”“平稳型”这几种典型画像，业务价值大，就选5。算法指标要结合业务解释，不能唯指标论。

粒子数40、迭代80代，对于96维、5个中心共480维的解空间来说不算大，但实际效果已经够好。如果你发现收敛不稳定，可以加到60粒子、120代，Matlab运行时间从几十秒变成几分钟，还能接受。

3.3 PSO-FCM核心代码段解读

PSO的主体是个双重循环，外层迭代，内层逐个粒子评估。最核心的是适应度函数，我单独写成子函数，方便把FCM一轮迭代封装进去。

% 初始化粒子位置（在[0,1]区间均匀随机，负荷归一化后中心也落在[0,1]） dim = C * D; lb = zeros(1, dim); ub = ones(1, dim); x = rand(numParticle, dim) .* (ub - lb) + lb; v = -vMax + 2 * vMax * rand(numParticle, dim); pbest = x; pbestFitness = inf(numParticle, 1); [gbest, gbestFitness] = deal(zeros(1, dim), inf); for gen = 1:maxGen w = wMax - (wMax - wMin) * gen / maxGen; for i = 1:numParticle center = reshape(x(i, :), C, D); fitness = fcmOneStepFitness(data, center, C, m); if fitness < pbestFitness(i) pbestFitness(i) = fitness; pbest(i, :) = x(i, :); end if fitness < gbestFitness gbestFitness = fitness; gbest = x(i, :); end end % 更新速度与位置 for i = 1:numParticle v(i, :) = w * v(i, :) ... + c1 * rand * (pbest(i, :) - x(i, :)) ... + c2 * rand * (gbest - x(i, :)); v(i, :) = max(min(v(i, :), vMax), -vMax); x(i, :) = x(i, :) + v(i, :); x(i, :) = max(min(x(i, :), ub), lb); % 边界吸收 end fprintf('Generation %d, gbestFitness = %.4f\n', gen, gbestFitness); end

那个fcmOneStepFitness函数是精髓：

function fitness = fcmOneStepFitness(data, center, C, m) % 计算隶属度矩阵 dist = pdist2(data, center); % N x C，欧氏距离 dist(dist < 1e-6) = 1e-6; % 防止除零 % 隶属度更新公式 u = dist.^(-2/(m-1)); u = u ./ sum(u, 2); % 一步中心更新（用FCM规则精炼中心） um = u.^m; newCenter = (um' * data) ./ sum(um, 1)'; % 用更新后的中心重新计算距离与目标函数 distNew = pdist2(data, newCenter); fitness = sum(sum(um .* distNew.^2)); end

这里有个小技巧：对距离加了个下界1e-6，防止样本和聚类中心完全重合时出现0除以0。虽然负荷连续值几乎不会完全重合，但数据里有大量相同值（比如归一化后很多0），保底一下总没错。

PSO跑完拿到gbest后，把它reshape成C*D的中心矩阵，作为FCM的初始中心，然后执行标准FCM迭代到收敛。这一步也可以用Matlab自带的fcm函数，但为了保持代码自主可控，我一般自己写，逻辑一样。

% FCM精炼迭代 center = reshape(gbest, C, D); u = zeros(N, C); for iter = 1:maxIter dist = pdist2(data, center); dist(dist < 1e-6) = 1e-6; u = dist.^(-2/(m-1)); u = u ./ sum(u, 2); um = u.^m; oldCenter = center; center = (um' * data) ./ sum(um, 1)'; if max(abs(center - oldCenter), [], 'all') < tol break; end end [~, label] = max(u, [], 2); % 每个样本划入隶属度最大的类

整个流程跑下来，我对比过纯FCM随机初始化跑20次取最优的结果，PSO-FCM的目标函数值更小，而且不同随机种子跑出来的最终聚类中心几乎一致。这就是全局搜索带来的稳定性收益。

3.4 结果可视化与行为分型

聚类完成后，最好把每类的平均负荷曲线画出来，一眼就能看出行为差异。Matlab的绘图代码很简单：

figure; t = 0:15/60:23.75; % 96点对应时间 colors = lines(C); for j = 1:C idx = find(label == j); avgCurve = mean(data(idx, :), 1); plot(t, avgCurve, 'Color', colors(j, :), 'LineWidth', 1.8); hold on; end xlabel('时间 (小时)'); ylabel('归一化负荷'); legend(arrayfun(@(x) sprintf('Class %d', x), 1:C, 'UniformOutput', false)); grid on;

画出来后，我会进一步统计每类的用户数占比、平均日用电量、峰谷时段特征，形成一张行为画像表。比如第1类平均曲线在8点到18点之间有高峰，负荷值稳定在0.6以上，晚上偏低，可以标注为“工作日白天型”，这类用户大概率是家中有白天用电设备或者居家办公人群。第2类曲线在19点到23点出现高峰，白天平缓，是典型的“晚间高峰型”。第3类曲线几乎全天在0.4到0.6之间窄幅波动，可能是“持续平稳型”，比如家里有长期运行的电器或者老人常住。有了这些标签，后面的需求响应、分时电价策略就有了数据支撑。

4. 居民用电行为分析的应用落地

4.1 聚类结果如何落到用电画像

很多文章讲到聚类出几个簇就停了，但实际业务要的是从簇到画像。我会为每个簇计算三个层面的特征：时间维度上的峰谷位置、形态复杂度、以及用户属性层面的交叉信息。

时间维度比较好算，把每类平均曲线找峰值点对应的时间窗，比如峰值出现在上午9点到11点，还是晚上7点到9点。形态复杂度可以用负荷曲线的标准差或者信息熵来度量，标准差大说明这类用户用电波动剧烈，标准差小说明用电平稳。交叉信息则需要结合营销系统的用户档案，比如用户是普通住宅、合租房、老人家庭等，虽然聚合平台不一定有，但如果能拿到一部分样本标签，可以做混淆验证，看看聚类结果是否对应实际人群。

我实际做过一次这样的分析，把某市3000户居民按96点曲线聚类成5类，再和用户档案中的“家庭人数”做交叉统计，发现“白天型”用户里三人及以上家庭占比明显高于平均值，而“平稳型”用户里一人户占比很高。这说明聚类结果不仅形状可解释，背后还有真实的行为驱动因素。

4.2 典型居民用电行为模式解释

下面是我在多次分析中见到的典型模式，不一定完全覆盖所有地区，但抽样数据里普遍会出现：

• 白天工作型：早8点后负荷上升，下午有所回落，晚高峰不明显。常见于家中有白天需要用电的设备，比如新能源车充电、净水器、鱼缸加热棒，也可能是有居家办公的人。
• 晚间高峰型：白天负荷低，18点以后快速攀升，20点到22点达到峰值。这是大多数上班族的特征，下班回家做饭、开灯、看电视、洗澡。
• 全天平稳型：负荷全天在均值附近小幅波动，没有明显峰谷。常见于老人家庭、有长期运行电器的住户，或者家里有多台智能设备在跑。
• 周末突出型：工作日曲线相对平缓，周末负荷明显高且持续时间长。这类用户平时早出晚归，周末在家时间长，可能是双职工家庭。
• 深夜活跃型：22点以后到凌晨2点负荷仍然高，其他时段一般。这类行为可能对应夜班人群或者习惯晚睡的年轻人。

如果你是把居民负荷用于电网规划，白天工作型和深夜活跃型对光伏消纳、储能配置的意义完全不同。白天工作型的午间负荷可以和光伏发电曲线匹配，提高自消纳率；深夜活跃型则要考虑夜间低谷填谷策略。

4.3 从聚类到需求响应、错峰管理

把用户分好型之后，最重要的应用是差异化需求响应。比如夏季高温天的晚高峰，电网要削峰，首选对象肯定是“晚间高峰型”用户，因为他们的晚高峰负荷集中，稍微调整一下空调温度或者推迟使用洗碗机，就能压出可观的容量。而“白天工作型”用户在晚高峰本来负荷就低，动员他们没有任何意义。

另一个应用是分时电价套餐设计。比如针对“深夜活跃型”用户推出深谷电价套餐，鼓励他们把洗衣、充电等活动进一步挪到凌晨，既帮电网填谷，也帮用户省钱。“全天平稳型”用户对峰谷电价不敏感，更适合固定电价或者阶梯电价。

我在做分析报告时，通常会输出一个“用户行为迁移潜力”的估值：计算每类用户在高峰时段的可转移负荷比例，结合该类用户的总户数，估算出理论削峰容量。这个数字对营销部门制定活动方案很有说服力。当然，实际响应率还要考虑用户意愿和激励力度，但聚类分析把“对谁推、推什么”这个问题解决了。

5. 常见问题与排错实录

5.1 FCM聚类结果不稳定

如果你跑了PSO之后发现结果还是不稳定，不要直接怀疑PSO失效，先检查两个地方。一是看PSO的适应度是不是早熟收敛了，即gbestFitness在迭代到第20代之后就不再变化；这时候加大粒子数或者调高wMax试试。二是看数据是不是存在极端离群点，离群点即使聚类中心不落在它附近，也会通过隶属度计算把某些正常样本拉偏。我在处理时先做了离群点剔除，用箱线图找每户曲线中超过上限的明显异常点，把包含大量异常跳变曲线的样本删掉，聚类稳定度立刻上来了。

5.2 粒子群早熟收敛

PSO早熟是原厂问题，表现为群体多样性消失，所有粒子挤在一起。排查方法很简单：每代输出gbest的位置标准差，如果标准差在迭代中期就跌到接近0，说明粒子全部聚集了。我的处理手段有三个。第一，把惯性权重w改成非线性递减，或者加入随机扰动，比如w = wMin + (wMax - wMin) * exp(-gen/maxGen)。第二，在速度更新里加入一个“反向扰动”因子，当gbest连续多代不更新时，给所有粒子一个随机速度冲击。第三，更简单粗暴，把粒子数翻倍，50粒子变100粒子，虽然时间翻倍，但收敛质量明显提升。

5.3 Matlab版本与工具箱问题

PSO-FCM代码只用到了基础矩阵运算和pdist2函数，pdist2在Statistics and Machine Learning Toolbox里。如果没有这个工具箱，可以用dist或自己写距离矩阵：

dist = sqrt(max(bsxfun(@plus, sum(data.^2,2), sum(center.^2,2)') - 2 * (data * center'), 0));

注意bsxfun在R2016b之后可以用隐式扩展替代，直接写成：

dist = sqrt(max(sum(data.^2,2) + sum(center.^2,2)' - 2 * data * center', 0));

这个写法不依赖任何工具箱。我遇到过用的Matlab版本比较老，没有pdist2，用这个替代版解决。另外，如果你的Matlab有并行工具箱，可以把粒子适应度评估的循环改成parfor，但要注意parfor里嵌套函数要支持。我实测4核并行只加速了1.7倍左右，因为PSO的更新公式是串行依赖gbest的，只有粒子评估能并行，加速比有限，性价比一般。

下面是我整理的排错速查表：

5.4 调参与避坑集锦

最后一个实战心得：不要急着用全部96维负荷曲线聚类。96维直接算距离，计算量大，而且相邻时间点高度相关，本质上是冗余信息。我一般先做降维，最简单的就是降采样成24点（每小时一个点），更专业一点用PCA或者负荷形态特征（峰谷差、平段占比、峰段占比、晚高峰占比、负荷率）来替代原始曲线。用24个点做PSO-FCM，维数从96降到24，收敛速度提升好几倍，聚类形状基本不变。

模糊指数m也不是只能取2。m越小，模糊性越低，聚类结果越接近硬聚类；m越大，各类边界越模糊。如果数据噪声大，m取2.2甚至2.5更稳。但m太大信息量损失也大，我一般用“分类系数”检查m取值，分类系数接近1说明隶属度很分明，接近1/C说明基本没法分类。m取2时，分类系数通常在0.8左右就够用了。

如果你用的是Matlab 2026b或者更早版本，代码在R2016b之后都能直接跑，不需要额外安装工具箱。Matlab下载安装这些事不细说了，遇到license报错基本都是环境配置问题，换个激活方式就好。真正麻烦的还是数据的质量和业务解释，代码只是工具。

这个项目做完后，我又往里面加了一个功能：每类用户单独画箱线图，看同类内部的一致性。如果某类用户的箱线图特别宽，说明这个类内还有细分空间，可以把C调大再跑一次。一步步试下来，最终出来的5类画像相当扎实。如果读者也有类似需求，建议从我的代码框架出发，把数据换成你手头的负荷表，调整一下C和m，大概率能直接落地。