☰
粒子群优化模糊C均值聚类:居民用电负荷行为分析实战
2026/10/6 16:46:34 网站建设 项目流程

我拿到一批居民智能电表的日负荷数据时,第一反应就是跑一遍模糊C均值聚类(FCM),想看看能不能把用户的用电习惯自动分成几类。结果Matlab里同一个脚本反复运行,聚类中心每次都不一样,有一户白天高负荷的家庭甚至在不同轮次里被分到了完全不同的类别。后来我才明白,FCM聚类处理这类问题,初值敏感和局部最优这两个毛病几乎是绕不过去的,必须用粒子群算法(PSO)这种全局寻优方法去改。这篇内容就把我完整跑通“粒子群算法优化FCM聚类”的过程拆开讲,包括数据怎么清洗、特征怎么选、Matlab代码怎么写、参数怎么调,以及最后的居民用电行为分型怎么解读,适合正在做负荷聚类、需求响应或用电模式识别相关研究的同学参考。

1. 为什么非要动FCM:居民用电聚类里最让人头疼的两件事

1.1 FCM聚类的软划分逻辑:它天生适合用电行为分析

普通K-means聚类是硬划分,一个样本只能属于某一个类,要么是A类要么是B类。但居民用电行为是一个典型的“不那么硬”的对象:一个家庭白天可能因为老人开电视、空调有较高的基荷,傍晚又因为下班做饭出现一个高峰,你说它到底是“白天型”还是“晚高峰型”?单看整体曲线,它可能两者的特征都有。

FCM聚类用一个隶属度矩阵U来解决这个问题。假设样本规模为n,聚类数为c,每个样本对每个类都有一个介于0到1之间的隶属度u_ij,而且满足每行的隶属度加起来等于1。优化目标是最小化这个函数:

J = ∑ᵢ ∑ⱼ u_ij^m · d_ij²

其中d_ij是第i个样本到第j个聚类中心的欧氏距离,m是模糊指数,一般取2。m越大,隶属度分布越“软”,也就是每个样本越容易被多个类别共同解释;m接近1,结果就越接近硬划分。这个特性用在居民用电行为上非常合适,因为一个用户的用电习惯本来就可能是多面性的,用软隶属度去刻画比硬标签更贴近现实。

1.2 初值敏感和局部最优:为什么每次跑FCM结果都不一样

FCM的求解本质上是交替迭代:先固定聚类中心更新隶属度,再固定隶属度更新聚类中心,反复循环直到目标函数变化小于阈值。问题是,这个目标函数J并不是一个凸函数,它有很多个局部极小点。交替迭代法本质上是一种坐标下降思路,最后收敛到哪个极小点,完全取决于迭代起点——也就是初始聚类中心选在哪里。

可以拿爬山来打比方:FCM相当于把一个登山者随机丢到山里的某个位置,让他只能往最近的高处爬,他最终能站到的山顶取决于出发点。如果初始中心本身选在了一个糟糕的位置,他可能爬到一个小土坡就以为到顶了,永远到不了真正的山峰。反映到实际聚类结果上,就是你用随机初始中心多跑几遍,会发现目标函数J值有大有小,聚类中心位置也会有明显偏移。我实测了一组500户样本的数据,普通FCM循环跑20次,J值波动最大能到8%左右,那段时间我对聚类结果心里是真没底。

1.3 PSO的思路:把“猜初始中心”变成“搜索最优中心”

粒子群算法的核心逻辑其实不复杂:维护一个粒子种群,每个粒子代表一组聚类中心候选解,粒子在解空间里飞来飞去,每次飞行都参考自身历史最佳位置和整个种群的历史最佳位置来更新速度,逐步逼近全局最优解。

把两者结合起来的做法就是:用PSO去搜索更好的聚类中心,再把这些中心交给FCM做精细迭代。PSO负责“看到全局”,FCM负责“局部打磨”,两者配合比单独用任何一个都稳得多。后面我会详细给出这套混合流程的Matlab实现,包括粒子编码方式、适应度函数构造和迭代细节。

2. 负荷数据清洗与特征设计:决定聚类结果的往往不是算法而是特征

2.1 原始负荷数据的清洗:缺失值和异常尖峰怎么处理

智能电表采集的原始负荷数据远没有理论上那么干净。我拿到的那批数据是15分钟一个采样点,一天96点,其中约有3%的采样点存在缺失或为0值。如果直接把0值当作真实负荷参与聚类,个别用户会被塑造成一个“夜间完全断电”的假形象,这会严重干扰聚类中心。

我的做法分三步:先做缺失值填充,对不超过连续三个缺失点的位置用相邻点的线性插值;对长时间段的连续缺失数据,直接用该用户同类型日(比如工作日对工作日)的中位数曲线补齐;再处理异常尖峰,也就是那些超过该用户整体负荷曲线95%分位数加上三倍四分位距的点,这些基本是采集误差或设备抖动,直接替换成局部中位数。清洗之后还需要看一眼整体曲线形状,确保没有离谱突变。

2.2 特征从哪来:96维曲线不能直接进聚类

很多新手拿到96点日负荷曲线就直接丢进FCM,这是个大坑。原因有两个:一是96维数据在欧氏距离计算下,维度灾难会稀释真实结构,导致聚类边界非常模糊;二是原始负荷曲线里包含大量冗余信息——相邻时段相关性极高,真正驱动行为分型的其实是峰谷时段、负荷水平和稳定性这几类指标。

我从原始曲线里提取了8个特征,每个都有明确的用电行为含义:

特征计算方式行为含义
日均电量全天96点之和总用电规模
峰段电量占比18:00-22:00电量/全天电量晚高峰依赖程度
谷段电量占比23:00-6:00电量/全天电量夜间用电活跃程度
最大负荷出现时刻全天最大负荷对应的时刻用电重心位置
负荷率全天平均负荷/全天最大负荷用电平稳程度
峰谷差率(峰段平均-谷段平均)/峰段平均日内波动幅度
夜间平均负荷0:00-5:00平均功率基础性夜间负载
工作日/休息日差异两类日负荷曲线的相关系数差行为规律稳定性

这8个特征把一天的负荷曲线压缩成了行为画像,既降低了维度,又让每个特征都具备业务解释能力。我自己实际跑下来的经验是:特征设计合理的情况下,PSO-FCM的轮廓系数能比直接在96维曲线上聚类高出30%左右,这个提升非常可观。

2.3 归一化:一个小步骤,影响巨大

FCM的距离计算用的是欧氏距离,这意味着特征的量纲直接决定聚类结果的偏向。举个例子,日均电量动辄几十千瓦时,而最大负荷出现时刻是0到24的小时数,两者数值差距太大,如果不做归一化,聚类基本就只看日均电量这一个特征了,其他特征等于白设计。

我推荐z-score标准化:对每个特征减去均值再除以标准差。为什么不用min-max归一化到0到1?因为min-max受异常值影响很大,比如某个用户某天出现一个异常尖峰,会把整个特征的极值拉偏,导致其他用户被压缩到很小的区间里。z-score对异常值更稳健。标准化之后,要不要保留原始曲线信息?可以保留一天的96维归一化曲线作为备选特征源,但聚类主体特征用那8个统计指标就够了。

3. PSO-FCM完整流程与Matlab核心代码实现

3.1 算法整体结构:外层PSO搜中心,内层FCM做打磨

把PSO和FCM组合起来,可以用“粗搜+精修”来理解。外层PSO负责在解空间里搜索聚类中心的位置,每次粒子更新位置后,以内层的FCM迭代作为局部优化器,跑几步让聚类中心在当前候选解附近更优,然后以FCM的目标函数值作为这个粒子的适应度。适应度越小,说明这组聚类中心越好。

%% 主程序参数设置 load('featureData.mat'); % 预处理后的特征矩阵,行表示用户 X = zscore(featureData); % z-score标准化 [n, d] = size(X); % n为样本数, d为特征数 c = 4; % 聚类数 m = 2; % FCM模糊指数 Dim = c * d; % 一个粒子代表 c 个聚类中心,共 c*d 维 nPop = 30; % 种群规模 MaxIter = 100; % PSO最大迭代次数 c1 = 1.5; % 个体学习因子 c2 = 1.5; % 社会学习因子 wMax = 0.9; wMin = 0.4; % 惯性权重线性递减范围 lb = repmat(min(X), 1, c); % 粒子位置下界 ub = repmat(max(X), 1, c); % 粒子位置上界

这里有一个关键设计:粒子的维度是c × d。比如聚类数c=4、特征维度d=8,那么每一个粒子是一个32维的向量,每8个连续维度投射成一行聚类中心,reshape成4×8的矩阵就是这一组候选中心。lb和ub来自特征矩阵的最小值和最大值,保证搜索出来的聚类中心不会跑到数据范围以外。

3.2 粒子编码与适应度函数:核心就一段代码

粒子编码的重点在于reshape的用法。Matlab的reshape按列优先填充,也就是说一个c×d矩阵在向量化时,是先把第一列的所有行拼起来,再拼第二列。所以从粒子向量还原中心矩阵的时候,必须用 reshape(pos, c, d) 而不是直接reshape成d×c,两个结果在数学上是转置关系,很多第一次写这个代码的人都会在这一步翻车,后面我会在踩坑部分专门说。

适应度函数做的事情是:给定一组聚类中心和样本数据,用FCM的目标函数公式计算J值。

function fit = calcFitness(X, centers, m) n = size(X, 1); c = size(centers, 1); dist = zeros(n, c); for j = 1:c diffMat = X - centers(j, :); dist(:, j) = sum(diffMat .^ 2, 2); % 第j个中心的欧氏距离平方 end dist = sqrt(dist); % 转成欧氏距离 invDist = dist .^ (-2 / (m - 1)); % 求隶属度公式的中间量 invDist(dist == 0) = 1e6; % 距离为0时给一个足够大的值防止除零 U = invDist ./ sum(invDist, 2); % 隶属度矩阵 fit = sum(sum((U .^ m) .* (dist .^ 2))); % FCM目标函数J end

要特别留意invDist里距离为0的情况。正常情况下一个样本距离某个聚类中心恰好为0的概率很小,但PSO在边界搜索时偶尔会出现某些粒子位置和数据点重合,这时候如果不做处理,U矩阵里会出现NaN,整个适应度就崩了。我用一个大的哨兵值替换掉0,既防止除零又不会过度扭曲隶属度。

3.3 FCM局部迭代:为什么只跑几步而不是跑到收敛

在PSO内部,对每个粒子都要调用一次FCM迭代。我的做法是控制内部迭代次数在5到10步。理由是:PSO本身的价值在于全局搜索,如果每个粒子都把FCM跑到完全收敛,会消耗大量时间,而且粒子一旦被局部信息主导,种群的多样性会迅速下降,PSO就退化成多起点FCM了。跑5步的意义相当于给当前候选中心做一轮局部“微调”,让适应度评估更准确。

function [U, centers, J] = fcmLocal(X, centers, m, maxIter) for iter = 1:maxIter dist = pdist2(X, centers); % 计算所有样本到中心的欧氏距离 invDist = dist .^ (-2 / (m - 1)); invDist(dist == 0) = 1e6; U = invDist ./ sum(invDist, 2); % 更新隶属度矩阵 centers = (U .^ m)' * X ./ sum(U .^ m, 1)'; % 更新聚类中心 J(iter) = sum(sum((U .^ m) .* (dist .^ 2))); if iter > 1 && abs(J(iter) - J(iter - 1)) < 1e-6 break; end end J = J(end); end

每次迭代由两步组成:固定聚类中心更新隶属度,固定隶属度更新聚类中心。聚类中心更新公式的写法是(U.^m)' * X除以U.^m按列求和,这里用到矩阵乘法的线性代数性质,直接把所有样本的加权平均中心一次性算出来,比写循环快得多。

3.4 PSO速度位置更新与边界反弹策略

PSO主循环里,惯性权重w采用线性递减策略,从0.9逐渐降到0.4。这是最常用的做法:前期w较大,粒子飞得快,有利于全局探索;后期w较小,粒子飞得慢,有利于在最优解附近精细搜索。

%% PSO主循环 vel = zeros(nPop, Dim); pos = lb + rand(nPop, Dim) .* (ub - lb); pbest = pos; pbestFit = inf(nPop, 1); for iter = 1:MaxIter w = wMax - (wMax - wMin) * iter / MaxIter; for i = 1:nPop centers = reshape(pos(i, :), c, d); [~, centers, ~] = fcmLocal(X, centers, m, 5); % 局部迭代5步 fit = calcFitness(X, centers, m); if fit < pbestFit(i) pbestFit(i) = fit; pbest(i, :) = reshape(centers, 1, []); end if fit < gbestFit gbestFit = fit; gbest = pbest(i, :); end end vel = w * vel + c1 * rand(nPop, Dim) .* (pbest - pos) ... + c2 * rand(nPop, Dim) .* (gbest - pos); pos = pos + vel; pos = max(pos, lb); pos = min(pos, ub); end

边界处理我采用截断法,也就是计算完新位置后,直接把超出边界的维度拉回边界值。还有一种做法是速度反向反弹,但实际测试下来截断法最简单稳定,反弹法在某些维度上会让粒子反复震荡,收敛速度反而变慢。

4. 参数设置与收敛性调试:跑不通和结果乱跳的真正原因

4.1 PSO参数怎么定:不是越大越好

粒子群算法的参数选择有很强的经验性,我的建议是从一组保守参数起步:种群规模nPop=30,最大迭代次数MaxIter=100,c1=c2=1.5,惯性权重从0.9线性降到0.4。这个组合在绝大多数中等规模数据集上都能得到稳定结果。

我做过一组扫参对比实验,用同一样本集跑了不同参数组合:

种群规模迭代次数目标函数J(相对值)运行耗时评价
105094.2约6秒欠收敛,J偏大
3010088.6约18秒稳健,推荐起步
5010088.3约30秒J降幅很小,耗时明显
8020088.1约50秒边际收益极低

可以看到,从30增加到80个粒子,J只下降了0.3%左右,但耗时翻了三倍。这说明参数不是越大越好,关键是找到收敛的阈值。对中小规模负荷聚类数据,nPop=30到40完全够用。

4.2 模糊指数m和聚类数c:两个最容易忽略的参数

模糊指数m默认取2.0,但值得多跑几组对比。m小于1.5时结果接近硬聚类,样本容易被强行分到某一个类,掩盖了用电行为的模糊性;m大于3时隶属度过于扁平,所有样本对每个类的隶属度都接近0.25,聚类中心之间的差异被稀释,类别就不清晰了。我建议在m=1.8、2.0、2.2三档之间对比轮廓系数,通常2.0附近表现最好。

聚类数c的选择是一门独立的学问。常用的辅助判断指标有三个:手肘法看目标函数J随c增大的下降拐点;轮廓系数取最大值的c;DBI指数取最小值的c。三者不一定指向同一个c,需要结合业务解释性综合判断。我用轮廓系数筛出来c=4最合适,因为分到4类时每类用户的负荷曲线特征差异明显,容易给出业务解释;分到5类时第5类样本太少且混合度高,轮廓系数反而下降。

4.3 收敛不稳定怎么办:用随机种子和多次运行一致性检验

PSO本身也带随机性,所以同一份数据跑两次结果有细微差别是正常的。但这种差别应该远小于普通FCM。我验证稳定性的做法是:用rng函数固定随机种子,跑一遍,记录gbestFit和最终聚类中心;然后换一个种子,再跑一遍,对比聚类中心的变化程度。

如果两次运行聚类中心差异很大,大概率是算法没收敛到位,或者粒子弹出边界太多导致搜索失效。这时优先检查三件事:迭代次数是否足够、内部FCM迭代步数是否太少、粒子的初始位置范围是否真的覆盖了整个数据空间。我遇到过一种情况是lb和ub设置错误,lb比ub还大,粒子位置初始化出来全是乱的,适应度全是NaN,这个错误让我排查了很久,后来发现是repmat方向写反了。

rng(42); % 固定随机种子,保证实验可复现 %% 算法运行...记录gbestFit_1 = gbestFit rng(2024); %% 换种子重跑...记录gbestFit_2

如果两次gbestFit相对误差在0.5%以内,这个解就可以放心使用。普通FCM要达到同样的一致性,我试过要加几十次随机重启,耗时比PSO-FCM还长,效果还不一定更好。

5. 聚类结果解读:四类典型居民用电行为画像

5.1 从聚类中心反推用户行为特征

算法跑完后,我用全局最优粒子还原出的聚类中心,再结合每类的特征均值和原始96点负荷曲线,给四类用户做了行为画像。这里注意不要只盯着聚类中心的数值,要把特征表还原回业务语义去解读。

我用的500户样本结果归纳如下:

用户类型日均电量峰段占比谷段占比负荷率峰谷差率典型特征描述
第一类 晚高峰集中型中等高,约0.55低偏低,约0.25高傍晚负荷飙升,夜间基本安静
第二类 全天活跃型较高中中高,约0.5低全天负荷平稳,无明显低谷
第三类 夜间主导型中低低高,约0.4偏低较高深夜和凌晨负荷明显
第四类 双峰均衡型高较高中中等较高早午和晚间各有一个高峰

第一类基本是上班族家庭,工作日下午五六点之后负荷迅速上升,锅碗瓢盆、照明电视一起用电,到夜里又安静下来。第二类更可能是家里全天有人的家庭,老人白天看电视、空调常开,负荷曲线像一条平稳的直线。第三类最值得关注,深夜负荷高通常意味着有新能源车在充电,或者用户作息与此前预设的常规时段完全不同。第四类是典型的双职工有娃家庭,早晨准备早餐有一波用电,晚上下班后再来一波,但白天有一段安静时间。

5.2 PSO-FCM和普通FCM的实际对比

我同样用这500户样本分别跑了K-means、普通FCM和PSO-FCM,对比了几个关键指标:

方法目标函数J轮廓系数多次运行中心变异度耗时
K-means104.50.21低(初始固定)约2秒
普通FCM92.30.33高,中心偏移约15%约3秒
PSO-FCM88.60.42低,中心偏移约2%约18秒

K-means虽然快,但轮廓系数明显偏低,硬划分对用电行为这种模糊对象确实水土不服。普通FCM轮廓系数比K-means好,但结果太不稳,每次跑都可能给出不同的用户分型。PSO-FCM在J值和轮廓系数上都优于普通FCM,更关键的是多次运行结果稳定,这一点在工程或研究场景里特别重要——你总不希望因为一次随机初始化的运气好坏而得出完全不同的结论。

5.3 软划分在业务里的价值:不止是贴标签

FCM输出的隶属度矩阵在业务上比硬标签更有用。比如第三类“夜间主导型”用户,其中有一批样本对第一类的隶属度也在0.4以上,说明这类用户既有夜间充电特征,同时也保留了一部分晚高峰用电习惯。如果只给他们贴一个“夜间型”标签,后续做需求响应时可能漏掉他们晚高峰的削峰潜力。

在实际应用时,我通常按最大隶属度作为硬标签做统计,同时保留完整隶属度向量做后续的精细化分析。比如可以设定一个“混淆系数”,如果某个样本的最大隶属度只比第二大的高不到0.1,就把它标记为混合型用户,单独建一个待观察名单。这种思路在电力客户分群、套餐推荐和有序用电方案制定时非常有价值。

6. 我用Matlab跑这套代码踩过的坑与效率优化建议

6.1 三个具体的坑:维度、距离矩阵和局部迭代步数

第一个坑就是前面提到的reshape方向问题。粒子向量还原聚类中心矩阵时,Matlab的reshape是列优先填充。假设粒子向量是[1,2,3,4,5,6,7,8]且c=2、d=4,那么reshape(pos,2,4)得到的矩阵第一行是[1,3,5,7],第二行是[2,4,6,8],而不是直觉中的[1,2,3,4]和[5,6,7,8]。这个顺序一旦搞错,聚类中心会乱掉,而且算法不会报错,只有在你回头对比中心数值时才发现问题,非常隐蔽。

第二个坑是pdist2在数据量大的时候内存飙升。如果样本量是十万级、特征维度几十维,pdist2会直接生成一个十万乘聚类数的中间矩阵,加上多次迭代调用,内存很快就见底了。我后来在大规模场景改用自己写的循环欧氏距离,或者分块计算,速度虽然慢一点但内存可控。

第三个坑是内部FCM迭代步数。我开始设置成20步,结果每个粒子都在做深度局部优化,粒子群多样性快速下降,最终结果跟普通FCM的多次随机重启差不多,PSO的全局优势完全被埋没。后来我把内部迭代步数降回5步,全局搜索能力才恢复过来。这个参数需要根据你的数据集规模调整,数据量大的话内部迭代1到3步就够了。

6.2 效率优化:向量化是第一优先级

Matlab的循环效率远低于向量化矩阵运算,尤其是在样本量大的时候。上面代码里聚类中心更新和隶属度更新都已经向量化,全部避开了逐样本循环。如果还想继续提速,可以考虑对每个粒子单独用parfor代替for做并行计算,因为不同粒子的适应度计算互不依赖,是典型的可并行任务。

还有一个提速技巧是给PSO一个好的起点,而不是完全随机初始化。我试过先用K-means或普通FCM快速跑几次,取最优结果作为初始种群中一部分粒子的位置,其余粒子随机生成。这种“半随机初始化”让算法前期的收敛速度快了很多,J值在迭代20轮时就已经接近完全随机初始化60轮的水平。

6.3 这套方法的扩展方向

PSO-FCM做居民用电行为分析,核心价值不在算法本身多新奇,而在于它把“聚类不稳定”这个实际问题解决掉了。后续可以扩展的方向也很多:如果数据粒度为15分钟且有多天连续数据,可以引入负荷曲线的形态特征做时间序列聚类;如果样本量到了几万户甚至百万户,可以把PSO-FCM作为离线分型工具,线上再用轻量模型做快速匹配;还可以把聚类结果作为输入特征,接入用电量预测或需求响应潜力评估模型,让行为分型真正落到业务决策里。

在实际操作中我还有一个体会:用PSO-FCM跑完聚类之后,一定要保留每次运行的目标函数J值和聚类中心,做成一条收敛曲线图。这张图不仅能帮你判断参数选得对不对,写论文或者做汇报的时候也是很有说服力的材料。判断标准很简单,收敛曲线应该在20到30轮之内快速下降,然后趋于平缓,如果到后期还在大幅跳动,就得回去检查边界处理和随机种子是不是有问题了。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询