基于PSO-FCM的居民负荷曲线聚类方法及Matlab实现

拿到两三百户居民半年逐时负荷曲线的时候，我脑子里冒出的第一句话是：这些曲线真的会说话。白天家里几乎一条平线、傍晚准时拉高的，大概率是朝九晚五的上班族；中午和傍晚各拱起一个峰、峰形还错落有致的，多半家里有人做饭；凌晨一两点还在高位波动的，作息和我们大多数人不太一样。问题是，数据量一旦到几百上千户，靠眼睛看根本不现实。

把“形状相似”的负荷曲线自动归到一起，本质上是个时间序列聚类问题。我当时用的是模糊C均值聚类（FCM），效果尚可，但有个硬伤让人很抓狂：它太吃初始聚类中心了。同一个数据集，随机初始化跑十次，轮廓系数能在 0.35 到 0.58 之间反复横跳，聚类结果每次都有出入。后来我把粒子群算法（PSO）放在 FCM 前面，先让粒子群在搜索空间里找到一组靠谱的初始中心，再交给 FCM 去做局部精修，结果一下就稳了。整套方案基于 Matlab 完整实现，能直接跑通，也能迁移到交通流量时段划分、商超客流画像这类曲线聚类场景。这篇就把我的完整思路、Matlab 代码骨架、参数依据和一路踩过的坑都写出来。

1. 居民用电行为分析：聚类问题到底难在哪

1.1 从负荷曲线到用户画像的常见闭环

居民用电行为分析不是新鲜课题。智能电表普及之后，供电公司手里攒着海量用户数据，但很多都躺在数据库里没有转化成业务价值。做完数据清洗后，最常见的第一个动作就是聚类：把用户按负荷曲线特征分成几类，每一类给出一个“行为画像”。这个画像往下能接很多东西，需求侧响应潜力评估、分时电价套餐设计、台区负荷预测、异常用电检测都可以基于用户分类结果展开。

这个闭环里，聚类是整个链路的第一个枢纽。前面数据准备得再好，聚类结果一团糟，后面的业务分析全站不稳。所以“聚类算法选型”并不是哪把锤子顺手就抡哪把的问题。

我当时的数据是逐时负荷，也就是每户每天 24 个点。把 30 天的曲线逐日拆分，每行是一条“日负荷曲线”，矩阵规模大概是 N 行乘 24 列。对这种数据聚类，传统做法有 K-means、层次聚类、DBSCAN，各有适用面，但做居民行为画像时 FCM 有明显优势。

1.2 为什么硬聚类不够用，FCM 模糊聚类更贴合用电场景

K-means 给每个样本分配一个硬标签，一就是一、二就是二。可居民用电行为天生就不是“非黑即白”的。我用模拟数据验证时造过一类曲线：白天有一个小鼓包，晚上又有一个大高峰。你要硬把它归到“白天型”还是“晚间型”，怎么归都别扭，因为它俩特征都占了。

FCM 的核心区别在于引入了隶属度矩阵 U。每个样本不再是硬性属于某个簇，而是以不同概率属于各个簇，所有隶属度之和为 1。比如某户曲线 0.6 像上班族、0.25 像居家型、0.15 像夜间活跃型，这个软分配结果对业务解释特别友好。后续做需求响应时，可以针对“高比例属于大功率可调型”的用户做重点筛选，而不是一刀切。

硬聚类和模糊聚类还有一个重要区别，就是抗噪性。负荷数据里噪声和异常点非常多，空调启停、临时出差、节假日作息变化，都会把曲线扯得变形。K-means 的硬划分对异常点敏感，一个离群值可能会拽动整个簇中心；FCM 在隶属度加权下，异常点的影响被模糊机制摊开了。

我在实验里做过一个直观对比：同一份数据，K-means 和 FCM 各跑 20 次，K-means 的轮廓系数波动范围是 FCM 的近两倍。当然 FCM 也不是没毛病，下面这个毛病才是我换 PSO 的原因。

1.3 FCM 的初始化敏感和局部最优困境，需要全局优化器介入

用最直白的话说，FCM 的迭代过程就是一个“坐标下降”式的局部搜索。它做的事情是反复交替更新隶属度矩阵 U 和聚类中心 V，让目标函数 J_m 逐步下降。数学上可以证明这个交替更新会让 J_m 单调下降，但“单调下降”不等于“降到全局最小”，它收敛到什么水平，相当大程度取决于初始 V 落在哪个“盆地里”。

我实测过一个很典型的现象：对同一个 180×24 的模拟数据集，设置聚类数 C=4，模糊指数 m=2，随机初始化跑 10 次标准 FCM，10 次落点各不相同。有些收敛到漂亮的簇结构，轮廓系数 0.55 以上；有些直接出现两个聚类中心挤在一起，相当于 4 类实际只分出了 2~3 类；还有一次迭代结束后某个簇只剩一个样本，等于分了个空壳类出来。

如果只是做一两次分析，多跑几次随机初始化挑最优结果也能凑合。可当数据量上来、或者要批量处理多个台区时，“碰运气”式初始化在工程上不可接受。所以我才把目光投向了粒子群算法，用一个无梯度的全局优化器去解决初始化问题。

2. 粒子群算法优化FCM的混合方案：设计细节与参数依据

2.1 三条技术路线对比：我为什么选“PSO粗搜+FCM精修”

把 PSO 和 FCM 结合，网上能搜到多种形态。我把它们归纳为三条路线，各有各的道理：

我最后选了路线A，也即“PSO粗搜+FCM精修”。核心原因很务实：FCM 的局部精化能力本来就非常强，它的交替迭代在局部收敛速度和精度上远超 PSO 的游走式逼近。真正缺的是好起点，而不是更好的局部搜索器。让 PSO 做自己擅长的事——全局大范围粗搜，让 FCM 做自己擅长的事——从好起点快速收敛，分工明确，计算量也可控。

这个判断怎么来的？我实测过路线B，粒子维度是 C×D=4×24=96，要让 PSO 在这个连续空间里靠速度位移慢慢逼近局部极值，光是迭代就要 500 代往上。而路线A里 PSO 只跑 80 代给一个差不多的起点，FCM 再迭代 30~50 次就收敛了，总计算量差了一个数量级。

2.2 粒子怎么编码成聚类中心，适应度函数怎么写

选好路线后，第一个要解决的问题是粒子编码。聚类中心矩阵 V 的形状是 C×D，C 是聚类数，D 是负荷曲线维度。粒子是 PSO 里的一维位置向量，所以要把 V 拉平。

% 假设 C = 4, D = 24，粒子维度就是 96 % 粒子位置 pos(1:96) 按行优先顺序 reshape 成 4×24 V_particle = reshape(pos, C, D);

这里编码顺序无所谓，只要所有粒子的编码和解码方式一致就行。我习惯按“簇1的全部维度、簇2的全部维度……”来排，代码里注释写清楚，调试时不至于绕晕。

适应度函数用的是 FCM 的目标函数 J_m：

[ J_m = \sum_{i=1}^{C}\sum_{j=1}^{N} u_{ij}^m \cdot |x_j - v_i|^2 ]

Matlab 实现里，对粒子位置解码得到 V，先算隶属度矩阵 U，再代进目标函数。注意一个经典坑：算隶属度时，某个样本可能到某个中心的距离恰好为 0，此时分母会出现除零，NaN 直接传染整个矩阵。我处理方式很简单，在距离分母上加一个 1e-10 的极小量：

for j = 1:N dist_j = zeros(1, C); for i = 1:C dist_j(i) = norm(X(j,:) - V(i,:)) + 1e-10; end U(j,:) = 1 ./ sum((dist_j' ./ dist_j).^(2/(m-1)), 1); end

2.3 PSO关键参数设置与物理含义

PSO 的标准更新公式大家都熟：

• 速度更新：( v_i = w v_i + c_1 r_1 (pbest_i - x_i) + c_2 r_2 (gbest - x_i) )
• 位置更新：( x_i = x_i + v_i )

关键是参数怎么给才稳定。我在这个项目里用的参数组合经过了几轮调整，最终确定的如下：

为什么用线性递减而不是固定 w？我对比过：固定 w=0.7 时，PSO 前期搜索步伐太小，容易早早陷入一个局部区域；线性递减则让粒子前期大步跨、后期小步挪，正好匹配优化过程“先扫全图、再抠细节”的需求。速度上限这件事同样重要，没有钳位的话，高维空间里粒子一步就能从搜索空间一头蹦到另一头，整个种群很快散开找不到北。

粒子位置初始化也讲究。我是用均匀随机数在边界范围内初始化，但这在 C×D 高维空间里会有个天然问题：全是随机起点，PSO 前几十代都在“探索茫茫大海”。后来我加了启动策略，把 K-means 快速聚类结果当成种子混入初始种群，效果提升非常明显。至于边界约束，我用的是“吸收模式”，粒子飞出去就拉回边界并清零速度，简单可靠。

2.4 FCM精化阶段的参数约定

PSO 给出全局最优位置后，reshape 成初始 V 交给 FCM。FCM 本身的参数需要稳定：模糊指数 m 我固定在 2，这是 FCM 最经典的取值；迭代上限 100 次；终止容差 1e-5。从多次实验来看，从 PSO 给出的起点出发，FCM 基本 30~50 次内就能收敛，设置 100 次是留足余量。

熔断条件写的是 J_m 相邻两次迭代的绝对差小于 1e-5，或者隶属度矩阵最大变化量小于阈值，二者满足其一就结束。有一点要注意：FCM 收敛后要检查是否有“空簇”——某个聚类中心对应的隶属度总和趋近于零。出现空簇通常意味着聚类数 C 给大了，或者初始中心质量仍然不行，这个在下文排错表格里详细展开。

3. Matlab完整可运行实现：数据预处理、代码框架、评价与可视化

3.1 输入数据长什么样，归一化和特征选择怎么选

先明确输入格式：X 是一个 N×D 的矩阵，N 是样本数，D 是曲线维度。我这里用 24 点日负荷曲线，因此 D=24；如果你的数据来自 15 分钟采集间隔，D=96，做法完全不用变。

数据归一化是第一个容易翻车的地方。我试过三种方式，结果差异很大：

• 最大最小归一化到 [0,1]：效果最稳，但会抹掉“用电量绝对值”的信息，只能看形态；
• z-score 标准化：保留了量级差异，高电量用户天然容易被单独聚成一簇；
• 除以当日峰值：突出峰谷形状，适合判断“峰谷时点”，但夜间高负载用户可能被错误分到“平峰型”。

我在做行为画像时，更关心的是形态而不是绝对电量，所以选了第一种。如果想要兼顾形状和电量，建议把“日用电总量”单独做一列特征加进去，而不是让聚类算法在混合尺度里自己摸索。

特征选择上，如果曲线有 96 个点，直接聚类也能跑，但冗余维度会增加计算负担且更容易踩局部最优。可以先用主成分分析压到 10~15 个主成分，或者手动构造统计特征，比如峰时占比、谷时占比、峰谷差、负荷率、早晚高峰位置。我建议的目标是：特征里既有“形态信息”又有“时点信息”，但总维度控制在 15 以内，这样 PSO 的粒子维度和收敛速度都理想。

3.2 主程序架构与核心函数拆解

我的整体代码结构分四层：数据准备、PSO搜索、FCM精化、评价可视化。下面是主程序骨架。

%% 主入口 clc; clear; close all; rng(42); % 固定种子，保证实验可复现 load('load_data.mat', 'X'); % X: N×24 日负荷曲线矩阵 X = normalize(X, 'range'); % 最大最小归一化到 [0,1] % 超参数 C = 4; m = 2; % 聚类数、模糊指数 Npop = 40; MaxIt = 80; % PSO种群与迭代 w_start = 0.9; w_end = 0.4; c1 = 2; c2 = 2; VmaxScale = 0.1; % PSO阶段 [gbest_pos, gbest_val] = pso_search(X, C, m, Npop, MaxIt, ... w_start, w_end, c1, c2, VmaxScale); % FCM精化阶段 V_init = reshape(gbest_pos, C, size(X,2)); [V_final, U_final, Jm_list] = fcm_refine(X, V_init, C, m, 100, 1e-5); % 评价与可视化 sil = silhouette(X, idx_from_U(U_final)); % 轮廓系数 plot_cluster_centers(V_final); % 簇中心曲线 plot_membership_heatmap(U_final); % 隶属度热力图

pso_search 函数内部就是一个标准 PSO 循环，适应度计算调用 fcm_objective。fcm_refine 是标准 FCM 交替迭代，但初始 V 由 PSO 提供。

function Jm = fcm_objective(X, V, m) C = size(V, 1); N = size(X, 1); U = update_membership(X, V, m); % 由当前中心算隶属度 Jm = 0; for j = 1:N for i = 1:C Jm = Jm + U(j,i)^m * norm(X(j,:) - V(i,:))^2; end end end

这里我单独把 update_membership 抽出来是为了两处复用：PSO 适应度计算用一次，FCM 精化迭代里也用一次，避免重复造轮子。这个面向重构的小习惯，在调试时省了不少时间。

3.3 聚类效果评价：轮廓系数、DBI、XB指数

聚类做完了，怎么判断效果好？这不能靠眼睛看散点图自嗨。我固定用三个指标交叉评估：

• 轮廓系数：衡量同类样本的紧密性和异类样本的分散性，范围 [-1,1]，越大越好。Matlab 自带 silhouette 函数，可以直接喂距离矩阵。
• 戴维斯-布尔丁指数（DBI）：类内散布小、类间距离大时 DBI 低，越低越好，适合对比不同 C 取值。
• 谢贝尼指数（XB）：模糊聚类专用指标，同时考虑紧致度和分离度，越小越好。这个指标对 FCM 的软分配结构更敏感，做模糊聚类对比时建议带上。

选聚类数 C 的时候，我的做法是遍历 C=2~6，每个 C 值跑完整套 PSO-FCM，画出三个指标的折线图，用“肘部法则”挑拐点。常见的错误是只看轮廓系数、忽略业务解释性。我曾经把 C=6 的指标调得挺好看，结果某两个簇中心的曲线形状几乎一样，只是电量大小不同，这种“统计上可分、业务上无用”的分法没有实际价值。

3.4 聚类结果可视化与业务解读

评价指标之外一定要可视化。最少做三张图：簇中心曲线叠加图、隶属度热力图、降维散点图。

簇中心曲线叠加图最直观，每条线是一类用户的典型负荷曲线，横轴是 24 个时段。看图说话就能给用户打标签。隶属度热力图是用 imagesc 画 U 矩阵，横轴是样本、纵轴是类别，颜色深浅代表隶属度，能直接看出哪些用户处于“模糊地带”——这类用户在后续业务里恰恰最值得单独关注。PCA 降维散点图则是辅助判断类与类之间的分离度，虽然 PCA 会损失部分信息，但能快速暴露“两类重叠严重”的问题。

4. 实验记录：模拟数据验证、调参踩坑与行为标签翻译

4.1 先造一套可控的模拟负荷数据，验证算法有效性

直接拿真实项目数据调试算法是一件很痛苦的事，因为数据里混杂了季节效应、节假日、缺失值和采集噪声，出了问题很难判断是算法问题还是数据问题。我习惯先造一套仿真数据，把“真实标签”掌握在手里，先让算法跑通，再上真实数据。

模拟数据生成逻辑很简单：手动定义三类曲线的峰谷模式，再叠加随机噪声。

% 三类基线的形状 t = 1:24; base1 = 0.2 + 0.8 * exp(-((t-19).^2)/8); % 晚间单峰型 base2 = 0.3 + 0.5 * exp(-((t-12).^2)/9) ... + 0.5 * exp(-((t-19).^2)/9); % 午晚双峰型 base3 = 0.2 + 0.7 * exp(-((t-23).^2)/14) ... + 0.3 * exp(-((t-4).^2)/10); % 夜间活跃型 % 每类生成 60 条，加噪声 X = [repmat(base1, 60, 1); repmat(base2, 60, 1); repmat(base3, 60, 1)]; X = X + 0.15 * randn(size(X)); % 15% 噪声 X = max(X, 0.05);

三类曲线本身的峰谷位置有明显区分度，但 0.15 的随机噪声已经足以让硬聚类出现不少错分样本。跑完 PSO-FCM 后，我对比了聚类标签与真实标签的匹配度，准确率在 95% 以上。这套代码的价值在于：它让你在完全受控条件下验证“算法逻辑正确”，而不是一上来就和真实数据里的脏问题纠缠不清。

4.2 三组对比实验看PSO带来了什么

为了搞清楚 PSO 到底贡献了什么，我做了三组对比：标准 FCM 随机初始化（重复10次取最优）、K-means 初始化后接 FCM、PSO-FCM。结果有意思的地方在稳定性。

PSO-FCM 的均值不是绝对最高的，但胜在下限高。最差一次也有 0.50，而随机初始化的最差能掉到 0.26，这是完全不可接受的分群结果。从工程实践看，算法稳定比算法惊艳更重要——你要的是每次跑都能复现，而不是偶尔爆发一次。

K-means 初始化方案耗时最低，性能也还可以，但它需要预设 K-means 自己的 K 值和随机状态，仍然不能彻底解决“初值运气论”。PSO-FCM 多花的 5 秒左右计算成本，换来的是稳定的输出质量，对于离线分析场景完全值得。

4.3 从聚类中心到行为标签：上班族、居家型、夜间活跃型的识别逻辑

聚类得到的簇中心是一条平均负荷曲线，要把数值曲线“翻译”成业务标签，我总结了一套规则化的解读思路：

这一步是让算法结果落地的关键。聚类算法只负责告诉你“数据上有几群”，业务人员才负责回答“这群人是谁”。我做了一个小工具函数，输入簇中心矩阵，输出每个簇的峰谷时点、峰值系数、峰谷差等统计量，自动生成描述文字。这个方法也可以挪到别的领域，比如商超客流的时段画像，逻辑完全一致。

5. 常见问题排查与实操避坑速查

5.1 高频报错与处理对策速查表

把我在项目中遇到的典型问题列个表，每个都是我实测踩过的，不是从文档里抄来的。

5.2 三个容易忽略但非常影响结果的细节

第一，粒子初始化时混入 K-means 种子。纯随机初始化在高维空间里简直是灾难，PSO 要花大量代数才能找到稍微像样的区域。我在初始化种群时，把 5 个粒子的位置设置成 K-means 跑出来的簇中心，其余粒子再随机扰动生成。效果立竿见影，最优适应度在头 10 代就明显下降。

第二，距离度量的取舍不要盲目追求高级。我试过用动态时间弯曲（DTW）替换欧氏距离，识别精度提升确实有，但计算复杂度从线性变成平方级，180 个样本跑一次就要等半天，放大到几千户根本不可行。如果你的曲线已经做了时段对齐和归一化，欧氏距离在绝大多数情况下够用，不要为了精度好看牺牲效率。

第三，模糊指数 m 不要默认一个值跑到底。m=2 是经典取值，但我发现对噪声水平不同的数据，m 在1.5到2.5之间的表现差异挺大。m 越小越接近硬聚类，容易过拟合噪声；m 越大分界面越模糊，类间差距被抹平。建议在 1.5、1.8、2.0、2.5 各跑一遍，用 XB 指数选最优，而不是拍脑袋定成 2。

我个人做完这个项目的体会是：PSO 和 FCM 的混合，真正的价值不在“精度比单一 FCM 高多少”，而在于它彻底解决了随机初始化带来的不确定性问题。从业务视角来说，一个每次都能稳定复现、并且能解释成用户画像的结果，远比一次偶然的高分更有工程意义。这套框架本质上还可以平移到任何“曲线聚类+初值敏感”的场景，比如交通流量的日时段划分、零售门店的客流时段画像、设备振动曲线的工况识别，核心思路和代码结构都不用大改。如果你也要在 Matlab 里做类似的事情，建议先从 24 点的小数据规模把 PSO 参数调顺，再把数据规模放大——小规模上暴露的问题，往往就是大规模跑飞的前兆。