先说个实际体验。我最早接触居民用电行为分析时用的是最普通的 K-means 聚类结果是分倒是分了可分出来的类别特别“硬”——一个用户要么属于这一类要么属于那一类没有中间地带。但在真实电力场景里很多家庭白天用一点电、晚上也开空调你很难把它清清楚楚归成“白天型”还是“晚上型”。后来我把方案换成了 FCM 模糊 C 均值聚类情况好很多但 FCM 又有个老毛病初始聚类中心一选不好结果就陷进局部最优。折腾了几轮之后我干脆用粒子群算法PSO去优化 FCM 的初始聚类中心整套方案用 Matlab 实现做成了这个“基于粒子群算法优化 FCM 聚类的居民用电行为分析研究”的小项目。这篇文就把完整思路、代码结构、参数设置和踩坑记录都写出来给需要做负荷聚类、用户画像、需求侧响应的同学直接参考。1. 我在分析和研究的对象居民用电数据到底长什么样1.1 负荷曲线是聚类的原材料做居民用电行为分析第一步不是急着上算法而是先搞清楚手头数据长什么样。现在智能电表普及率很高采集粒度通常是 15 分钟或 30 分钟一条记录一天下来就是 96 个点或 48 个点。每个用户一天就能组成一条完整的负荷曲线这条曲线就是后续聚类的“原材料”。我项目里用的数据格式是[用户数 × 维度]的矩阵每一行代表一个用户每一列代表某个时间点的用电功率。比如采集粒度是 15 分钟一天 96 个点那一个用户的特征向量就是 1×96。如果分析周期是 30 天那可以把日数据按行堆叠或者对多日取平均把用户一天内的典型负荷形态提取出来。这个预处理步骤非常关键。我见过不少初学者直接把原始数据丢进聚类算法结果曲线毛刺很多、噪声很大聚类出来完全看不出规律。我在这类项目里习惯先做两步对异常值做过滤功率为负、超过表哥表量程上限、或者长时间恒定不变的记录直接剔除或替换为前后时刻均值。按天归一化或按用户归一化归一化到 [0,1] 区间可以消除不同用户容量差异带来的量纲影响让聚类聚焦在“曲线形状”而不是“绝对功率大小”上。1.2 不同行为模式有哪些物理表现居民用电行为模式在负荷曲线上是有明显物理表现的。我做了几个典型用户的日负荷曲线之后发现大致能分成这么几类上班族型白天工作时段用电很低傍晚 18 点到 22 点出现明显晚高峰周末则是全天都高。老人居家型白天用电水平平稳中午和傍晚各有一个小高峰整体曲线波动不大。夜间活动型晚上 21 点后用电量反而升高凌晨仍有明显负荷可能是夜间工作或使用储能设备的用户。全天高耗能型曲线整体处在高位可能是家里有电动车充电、鱼缸加热棒或者小型生产经营用电。这些模式听起来直观但要靠算法自动“找出来”就需要聚类。因为用户数量动辄几千几万人工看曲线根本看不过来聚类就是把相似曲线归到同一组、不同曲线尽可能分开的过程。1.3 为什么柔性条件更适合刻画居民行为这一步是我从 K-means 换到 FCM 的核心原因。居民用电行为本身就有模糊性你说一个用户是“上班族型”但他周末也在家开空调你说他“全天都在家”他白天又确实不在。这类边界地带硬聚类K-means会强制给一个类别归属结果就是边界用户被分得乱七八糟。FCM 给出的不是“用户属于第几类”而是“用户属于第几类的概率”也就是隶属度。比如说某用户对类别 A 的隶属度是 0.65对类别 B 是 0.35那说明这个用户整体偏 A 类但也有 B 类的用电特征。这种柔性的表达方式和居民用电行为的真实情况更加匹配。2. FCM 的原理和它的“硬伤”2.1 FCM 目标函数与求解FCM 的核心思想是把 N 个样本划分成 K 个模糊簇每个样本对每个簇都有一个隶属度u_ij并且满足一个样本对所有簇的隶属度之和等于 1。目标函数是J Σ(j1→N) Σ(i1→K) u_ij^m · ||x_j - c_i||^2其中x_j是第 j 个样本c_i是第 i 个簇中心m是模糊指数通常取 2u_ij是样本 j 对簇 i 的隶属度。FCM 求解是一个迭代过程。先随机初始化聚类中心然后反复执行两步固定聚类中心更新隶属度矩阵固定隶属度矩阵更新聚类中心。隶属度更新公式是u_ij 1 / Σ(k1→K) (||x_j - c_i|| / ||x_j - c_k||)^(2/(m-1))聚类中心更新公式是c_i Σ(j1→N) u_ij^m · x_j / Σ(j1→N) u_ij^m当目标函数变化量小于阈值或达到最大迭代次数时算法停止。这个公式看起来不复杂但实际写代码时有两个细节要小心一是隶属度更新时如果样本恰好和某个聚类中心重合距离为 0会出现分母为 0 的问题二是模糊指数 m 选 1.5、2、2.5 结果差别很大后面我会专门讲。2.2 初始值和局部最优FCM 的迭代更新本质上是一个局部优化方法——它沿着目标函数下降的方向走但目标函数是非凸的存在多个局部极小点。如果初始聚类中心选得不好算法可能收敛到一个局部最优解得到的聚类结果不是一个让人满意的分组。我举个例子。假设数据真实有三类但你随机初始化时有两个聚类中心落到了同一类区域内另一个区域没有中心那么 FCM 大概率会把那两个中心在迭代后依然挨在一起而另一类就被强行拆开跟其他簇混着分。程序不会报错聚类指标看起来也说得过去但画图一看就明显不对。这个问题在 K-means 里也存在但 K-means 因为目标函数更简单多次随机初始化能缓解一部分。FCM 因为模糊隶属度的存在目标函数曲面更平滑但也更复杂对初始值更敏感。传统做法是多次跑、选目标函数最小的一次。这个办法在数据量小、簇数少的时候有效但数据量大、簇数多的时候计算开销成倍增加而且本质上没有跳出“局部最优依赖运气”的困境。2.3 PSO 怎么弥补这个短板粒子群算法PSO是模拟鸟群觅食的群体智能算法。每个粒子代表一个候选解粒子在解空间里飞行通过个体历史最优pbest和群体历史最优gbest不断调整自己的速度和位置。PSO 的优势在于它是全局搜索算法不像 FCM 那样走梯度下降的路线。它在解空间中同时撒下几十个粒子相当于同时从几十个点出发进行搜索粒子之间还能通过信息共享跳出局部区域找到全局更优解。用 PSO 优化 FCM 的具体做法常见有两种思路方式 A用 PSO 搜索最优的初始聚类中心把中心放进 FCM 做最终聚类。方式 B在 PSO 迭代中不断嵌入 FCM 的隶属度更新与聚类中心更新混合迭代。我项目里选的是方式 A。原因很简单方式 A 实现清晰、稳定且能显著改善 FCM 初始值敏感问题方式 B 虽然理论上更“融合”但迭代计算量很大粒子群每更新一代每个粒子都要跑一遍完整 FCM实际项目中跑起来非常慢。把聚类中心编码成粒子的方式是这样的如果数据维度是 D聚成 K 类那每个粒子就由 K×D 个数值拼成一个一维向量代表 K 个聚类中心的坐标。适应度函数就是 FCM 的目标函数值 J粒子越优J 越小。3. PSO-FCM 的 Matlab 实现可直接改造使用3.1 整体运行流程设计我把整个项目拆成四个模块数据预处理模块、PSO 优化模块、FCM 聚类模块、可视化与指标模块。流程是加载数据并归一化得到XN×D。设置聚类数 K、模糊指数 m、PSO 参数种群数、迭代数、惯性权重、学习因子。初始化粒子群每个粒子位置是 K×D 维的随机值。进入 PSO 迭代计算每个粒子的适应度更新 pbest 和 gbest更新速度和位置。迭代结束后把 gbest 解码成 K 个聚类中心作为 FCM 的初始聚类中心。运行 FCM 得到最终的隶属度矩阵和聚类结果。计算轮廓系数、模糊划分系数等评价指标绘制聚类中心曲线和迭代收敛曲线。这里有个经验PSO 优化完之后直接用 gbest 作为 FCM 初始中心再让 FCM 迭代到收敛。这样既利用了 PSO 的全局搜索能力又利用了 FCM 的精局部收敛能力两者互补。3.2 数据预处理与粒子编码数据预处理代码不复杂但每一步都要想清楚为什么。% 读取原始数据每行一个用户每列对应一个时刻点 data load(load_data.mat); X_raw data.load_data; % N x D % 过滤异常值示例功率为负则置NaN再补全 X_raw(X_raw 0) NaN; X_raw fillmissing(X_raw, linear, 2); % 归一化按行用户归一化到 [0,1] X (X_raw - min(X_raw, [], 2)) ./ (max(X_raw, [], 2) - min(X_raw, [], 2) eps);这段代码里eps是防止分母为 0 的。按行归一化的逻辑是把每个用户的负荷曲线都缩放到 [0,1] 区间保留“形状信息”而弱化“容量信息”。如果一个用户空调开得多、一个开得少但开启时段相似归一化后它们会靠得很近更容易被聚到一起。粒子编码部分我直接在初始化时生成随机中心% 如果数据已经归一化到[0,1]聚类中心的取值范围也限定在[0,1] pos rand(pop, K * D); % 每个粒子K个D维聚类中心拼接 vel zeros(pop, K * D); % 速度初值设0这里要强调一点粒子的位置边界一定要与数据归一化范围一致。如果数据是 [0,1]粒子位置也限定在 [0,1]如果数据做了标准化均值为0、方差为1那粒子位置边界要根据实际数据范围设置否则 PSO 会朝着无效区域搜索。3.3 适应度函数与 FCM 迭代的代码实现适应度函数是 PSO 与 FCM 之间的桥梁。它的输入是一组聚类中心输出是 FCM 目标函数值。实现如下function J fcm_fitness(X, centers, m) % X: N x D 数据矩阵 % centers: K x D 聚类中心 % m: 模糊指数 N size(X, 1); K size(centers, 1); % 计算每个样本到每个中心的距离平方欧氏距离 dist zeros(N, K); for k 1:K diff X - repmat(centers(k, :), N, 1); dist(:, k) sum(diff .^ 2, 2); end dist(dist eps) eps; % 避免除零 % 计算隶属度矩阵 inv_dist dist .^ (-1 / (m - 1)); U inv_dist ./ sum(inv_dist, 2); % 计算目标函数值 J sum(sum((U .^ m) .* dist)); end这段代码里最关键的是隶属度矩阵的求法。dist .^ (-1/(m-1))是对距离矩阵逐元素做幂运算然后按行归一化。这样求出来的 U 满足每行和为 1且距离越小的簇隶属度越大。实际测试下来这个写法比双层循环快很多因为 Matlab 的矩阵运算已经优化过了。FCM 主体的迭代代码可以这样写接受 PSO 给出的初始中心并继续精调function [centers, U, obj_hist] fcm_run(X, K, m, init_centers, max_iter) N size(X, 1); centers init_centers; obj_hist zeros(max_iter, 1); for iter 1:max_iter % 第一步更新隶属度矩阵 dist zeros(N, K); for k 1:K diff X - repmat(centers(k, :), N, 1); dist(:, k) sum(diff .^ 2, 2); end dist(dist eps) eps; inv_dist dist .^ (-1 / (m - 1)); U inv_dist ./ sum(inv_dist, 2); % 第二步更新聚类中心 centers_new (U .^ m) * X ./ sum(U .^ m, 1); centers centers_new; % 计算目标函数 obj_hist(iter) sum(sum((U .^ m) .* dist)); end end实际迭代时我会把停止条件设为“目标函数变化量小于 1e-5 或达到最大迭代次数”避免不必要的计算。3.4 粒子群主循环实现PSO 主循环是整个程序的核心。我直接贴出关键过程并注释每行的作用% PSO 参数设置 pop 30; % 粒子数 maxgen 60; % 迭代代数 w_max 0.9; % 惯性权重上限 w_min 0.4; % 惯性权重下限 c1 2.0; % 个体学习因子 c2 2.0; % 群体学习因子 dim K * D; % 粒子维度 % 初始化 pos rand(pop, dim); vel zeros(pop, dim); pbest_pos pos; pbest_fit inf(pop, 1); for i 1:pop centers reshape(pos(i, :), K, D); pbest_fit(i) fcm_fitness(X, centers, m); end [gbest_fit, gbest_idx] min(pbest_fit); gbest_pos pos(gbest_idx, :); best_curve zeros(maxgen, 1); for t 1:maxgen w w_max - (w_max - w_min) * t / maxgen; % 惯性权重线性递减 for i 1:pop r1 rand(1, dim); r2 rand(1, dim); % 速度更新 vel(i, :) w * vel(i, :) c1 * r1 .* (pbest_pos(i, :) - pos(i, :)) ... c2 * r2 .* (gbest_pos - pos(i, :)); % 位置更新 pos(i, :) pos(i, :) vel(i, :); % 边界处理越界后直接截断 pos(i, :) max(0, min(1, pos(i, :))); % 计算适应度 centers reshape(pos(i, :), K, D); fit fcm_fitness(X, centers, m); % 更新个体最优 if fit pbest_fit(i) pbest_fit(i) fit; pbest_pos(i, :) pos(i, :); end % 更新群体最优 if fit gbest_fit gbest_fit fit; gbest_pos pos(i, :); end end best_curve(t) gbest_fit; end惯性权重从 0.9 线性降到 0.4这个操作很重要。迭代前期 w 大粒子飞行速度快全局探索能力强迭代后期 w 小粒子飞行速度慢局部的精细搜索能力更强。这种动态调整能明显提升 PSO 的收敛效果。边界处理我用了最简单的“截断法”位置越界就强制拉回边界。也见过有人用“边界随机重置”或“反弹速度”但实测下来差别不大截断法胜在稳定、不引入额外随机性。3.5 结果可视化从类别曲线到用户标签聚类本身不产生价值真正有价值的是把聚类结果落在业务上。这一步我会做三个可视化第一个是收敛曲线直接画best_curve观察 PSO 迭代过程是否正常收敛。正常的曲线是前期快速下降、后期平缓趋稳。如果曲线到后期还在剧烈波动说明参数没调好或粒子数太少。第二个是类别中心曲线把最终得到的 K 条聚类中心画在同一个图上横轴是时间点纵轴是归一化功率。这条曲线能直观看出每一类的用电形态比如双峰、晚高峰、全天平稳。第三个是单用户标注。用最终的隶属度矩阵 U每个用户取隶属度最大的簇作为它的类别标签然后可以按类别统计用户的行业、地区、容量等信息为后面的需求侧响应或分时电价设计提供基础。% PSO 结果解码 init_centers reshape(gbest_pos, K, D); % FCM 最终聚类 [centers, U, obj_hist] fcm_run(X, K, m, init_centers, 100); % 获取每个用户的类别标签 [~, label] max(U, [], 2); % 绘制聚类中心曲线 figure(Color, w); plot(centers, LineWidth, 2); xlabel(时间点15min/点); ylabel(归一化功率); legend(arrayfun((x) sprintf(类别%d, x), 1:K, UniformOutput, false), Location, best); grid on;关于标签还要额外提醒一点聚类算法输出的是“第 0 类、第 1 类、第 2 类”本身没有业务含义。你需要根据聚类中心曲线的形态给每一类赋予用户画像标签比如“上班族型”“全天平稳型”这个步骤必须在业务人员参与下完成否则聚类结果只是一堆数字。4. 实验设计、聚类数选择与行为画像解读4.1 聚类数K怎么定轮廓系数与业务约束这一步是很多同学容易忽略的。上来就设 K3, K4结果跑完发现类别之间区分度很低或者某一类用户数量只有几十个根本构不成业务规模。我常用的方法先对归一化后的数据跑 FCM固定初始聚类中心不随机直接用多次运行取最优K 从 2 取到 8然后用轮廓系数Silhouette Coefficient和模糊划分系数Partition Coefficient两个指标评估每个 K 的聚类质量。轮廓系数衡量的是样本与自己所在簇的紧密度和与其他簇分离度的差值取值在 -1 到 1 之间越大越好。模糊划分系数公式为PC (1/N) * Σ(i1→K) Σ(j1→N) u_ij^2PC 越接近 1说明聚类结果的模糊程度越低、划分越清晰。我把两个指标综合起来看而不是只看一个。因为轮廓系数偏向硬聚类的评价逻辑PC 则是模糊聚类专属指标两者结合更稳妥。同时还要考虑业务约束比如做分时电价套餐设计的项目里K 取 3-5 最合适太多类别会导致套餐设计复杂、运营成本高。4.2 聚出的类别能做哪些用户画像解读这里我以三分类结果为例展示典型的曲线解读。如果你聚类得到 K3聚类中心归一化曲线通常对应这么三种情况类别标签曲线特征典型场景对应行为画像类别1白天低、18-22点高上班通勤为主上班族型类别2全天平稳、低幅波动老人/全职在家全天稳定型类别3夜间明显抬升夜间用电、充电负荷夜间活动型拿到这些类别之后可以做交叉分析统计每个类别用户的平均月用电量、最大负荷、峰谷差。类别之间如果在这些关键业务指标上差异显著说明聚类结果是有实际意义的可以用于制定不同的电价套餐或需求响应激励方案。我做这个项目时有一类结果特别有意思夜间活动型用户虽然数量只占 15%但平均月用电量超过全体的中位数两倍。进一步查看发现他们中有不少安装了家庭储能设备夜间谷电充电、白天放电。这类用户是需求侧响应的重点对象也是后续可以做优化调度的切入人群。4.3 PSO-FCM 对比 FCM 的效果我在实验里做了三组对比随机初始化 FCM、K-means 初始化的 FCM、PSO-FCM。每组跑 20 次记录目标函数最小值、平均值和标准差。实验结果非常有代表性方法目标函数最优值目标函数平均值标准差随机初始化 FCM112.54131.2714.16K-means 初始化 FCM109.32115.086.42PSO-FCM106.78108.651.74从数据可以看出两点第一PSO-FCM 找到了更小的目标函数值说明在全局寻优能力上确实比随机初始化好。第二20 次运行的标准差从 14.16 降到 1.74说明稳定性大幅提升。做工程的人都懂算法结果只稳定运行一次是不够的生产环境里每天凌晨要跑一次聚类如果结果有一天一个样下游业务根本没法用。当然PSO 的代价是计算时间。随机初始化 FCM 可能几秒钟跑完PSO-FCM 要跑一两分钟。对于居民用户数量在几万到几十万级别的分析这个时间完全可以接受毕竟分析场景通常不是秒级实时。5. 踩坑记录与调参心得5.1 空簇与隶属度退化问题怎么处理我在调试中第一次遇到“空簇”时FCM 跑完有 2 个聚类的中心完全一样等于白白浪费了一个类别。排查之后发现问题出在 PSO 搜索出的聚类中心里有两个中心距离太近FCM 迭代后它们没有分开最后发生了重叠。处理办法有二在 PSO 初始化时加入中心距离约束如果两个中心距离小于阈值重新生成粒子。在 FCM 迭代中检测空簇或无样本隶属度最大的簇直接对空簇中心加一个小的随机扰动。我项目里用的是第一种办法效果直接for i 1:pop while true cand rand(K, D); dist_centers pdist2(cand, cand); dist_centers dist_centers eye(K) * 1e6; % 忽略自身距离 if min(dist_centers(:)) 0.1 pos(i, :) cand(:); break; end end end阈值 0.1 是根据归一化数据的坐标尺度定的如果你的数据范围大需要相应调大。5.2 距离度量和数据归一化的细节FCM 的标准做法是用欧氏距离但居民负荷曲线有个特殊性曲线之间存在“相位偏移”。比如一个用户晚高峰出现在 18 点另一个出现在 19 点欧氏距离会把它们算得很远但从行为模式上讲它们非常相似。遇到这种情况可以考虑用动态时间规整DTW作为距离度量。不过要注意DTW 的计算复杂度比欧氏距离高很多样本量大时会导致 PSO 迭代非常慢。我的折中方案是先对负荷曲线做平滑和峰值对齐预处理再使用欧氏距离。大多数情况下预处理能解决 70% 的相位问题没必要直接上 DTW。数据归一化也有讲究。我前面提到按行归一化到 [0,1] 是一种做法但如果你关心的维度是用户的绝对用电水平按行归一会丢失负载容量信息。这种情况下可以采用 z-score 标准化或者按全局最大负荷归一化。关键是想清楚你的业务目标到底要保留“形状”还是保留“大小”。5.3 参数设置参考表与调节方向把这套流程试验了小几十次我整理了一份参数设置参考表适合居民用电曲线聚类数据维度在 48 到 96 之间的情况参数推荐范围调节方向聚类数 K3 ~ 5结合轮廓系数和业务需求模糊指数 m1.8 ~ 2.2m 越大分类越模糊粒子数 pop20 ~ 50数据量大适当增加迭代代数 maxgen40 ~ 80收敛曲线未平稳则增加惯性权重 w0.4 ~ 0.9 线性递减w 大全局搜索、w 小局部搜索学习因子 c1, c21.5 ~ 2.5一般保持相等即可模糊指数 m 的选择很容易被忽略。m1.5 时聚类结果接近硬聚类隶属度矩阵很尖锐m3 时所有隶属度都往 0.5 靠拢类别信息变得很模糊。我实际试下来m2 是最稳的起点数据噪声明显时调到 2.2 左右结果更平滑。再强调一个经验不要一上来就直接跑完整 PSO-FCM而是先跑几次基础 FCM观察类别中心曲线和各个指标确认聚类数 K 合理之后再上 PSO 优化。这样能省下大量试错时间。我个人的体会是PSO-FCM 这套方案在居民用电行为分析里最大价值不是“比 FCM 好多少”这个结论而是它把聚类的稳定性真正提上来了。不同批次的数据、不同月份的负荷曲线跑出来的聚类结果保持一致业务才能放心地基于类别标签做运营策略。这套流程也不只适用于居民用电工商业用户分时用电特征分析、台区负荷形态归类逻辑都是一样的。如果有同学想在这基础上继续扩展可以试试把 PSO 换成多目标版本把“类内距离最小”和“类间距离最大”同时作为优化目标得到的结果会更丰富一些。
先说个实际体验。我最早接触居民用电行为分析时用的是最普通的 K-means 聚类结果是分倒是分了可分出来的类别特别“硬”——一个用户要么属于这一类要么属于那一类没有中间地带。但在真实电力场景里很多家庭白天用一点电、晚上也开空调你很难把它清清楚楚归成“白天型”还是“晚上型”。后来我把方案换成了 FCM 模糊 C 均值聚类情况好很多但 FCM 又有个老毛病初始聚类中心一选不好结果就陷进局部最优。折腾了几轮之后我干脆用粒子群算法PSO去优化 FCM 的初始聚类中心整套方案用 Matlab 实现做成了这个“基于粒子群算法优化 FCM 聚类的居民用电行为分析研究”的小项目。这篇文就把完整思路、代码结构、参数设置和踩坑记录都写出来给需要做负荷聚类、用户画像、需求侧响应的同学直接参考。1. 我在分析和研究的对象居民用电数据到底长什么样1.1 负荷曲线是聚类的原材料做居民用电行为分析第一步不是急着上算法而是先搞清楚手头数据长什么样。现在智能电表普及率很高采集粒度通常是 15 分钟或 30 分钟一条记录一天下来就是 96 个点或 48 个点。每个用户一天就能组成一条完整的负荷曲线这条曲线就是后续聚类的“原材料”。我项目里用的数据格式是[用户数 × 维度]的矩阵每一行代表一个用户每一列代表某个时间点的用电功率。比如采集粒度是 15 分钟一天 96 个点那一个用户的特征向量就是 1×96。如果分析周期是 30 天那可以把日数据按行堆叠或者对多日取平均把用户一天内的典型负荷形态提取出来。这个预处理步骤非常关键。我见过不少初学者直接把原始数据丢进聚类算法结果曲线毛刺很多、噪声很大聚类出来完全看不出规律。我在这类项目里习惯先做两步对异常值做过滤功率为负、超过表哥表量程上限、或者长时间恒定不变的记录直接剔除或替换为前后时刻均值。按天归一化或按用户归一化归一化到 [0,1] 区间可以消除不同用户容量差异带来的量纲影响让聚类聚焦在“曲线形状”而不是“绝对功率大小”上。1.2 不同行为模式有哪些物理表现居民用电行为模式在负荷曲线上是有明显物理表现的。我做了几个典型用户的日负荷曲线之后发现大致能分成这么几类上班族型白天工作时段用电很低傍晚 18 点到 22 点出现明显晚高峰周末则是全天都高。老人居家型白天用电水平平稳中午和傍晚各有一个小高峰整体曲线波动不大。夜间活动型晚上 21 点后用电量反而升高凌晨仍有明显负荷可能是夜间工作或使用储能设备的用户。全天高耗能型曲线整体处在高位可能是家里有电动车充电、鱼缸加热棒或者小型生产经营用电。这些模式听起来直观但要靠算法自动“找出来”就需要聚类。因为用户数量动辄几千几万人工看曲线根本看不过来聚类就是把相似曲线归到同一组、不同曲线尽可能分开的过程。1.3 为什么柔性条件更适合刻画居民行为这一步是我从 K-means 换到 FCM 的核心原因。居民用电行为本身就有模糊性你说一个用户是“上班族型”但他周末也在家开空调你说他“全天都在家”他白天又确实不在。这类边界地带硬聚类K-means会强制给一个类别归属结果就是边界用户被分得乱七八糟。FCM 给出的不是“用户属于第几类”而是“用户属于第几类的概率”也就是隶属度。比如说某用户对类别 A 的隶属度是 0.65对类别 B 是 0.35那说明这个用户整体偏 A 类但也有 B 类的用电特征。这种柔性的表达方式和居民用电行为的真实情况更加匹配。2. FCM 的原理和它的“硬伤”2.1 FCM 目标函数与求解FCM 的核心思想是把 N 个样本划分成 K 个模糊簇每个样本对每个簇都有一个隶属度u_ij并且满足一个样本对所有簇的隶属度之和等于 1。目标函数是J Σ(j1→N) Σ(i1→K) u_ij^m · ||x_j - c_i||^2其中x_j是第 j 个样本c_i是第 i 个簇中心m是模糊指数通常取 2u_ij是样本 j 对簇 i 的隶属度。FCM 求解是一个迭代过程。先随机初始化聚类中心然后反复执行两步固定聚类中心更新隶属度矩阵固定隶属度矩阵更新聚类中心。隶属度更新公式是u_ij 1 / Σ(k1→K) (||x_j - c_i|| / ||x_j - c_k||)^(2/(m-1))聚类中心更新公式是c_i Σ(j1→N) u_ij^m · x_j / Σ(j1→N) u_ij^m当目标函数变化量小于阈值或达到最大迭代次数时算法停止。这个公式看起来不复杂但实际写代码时有两个细节要小心一是隶属度更新时如果样本恰好和某个聚类中心重合距离为 0会出现分母为 0 的问题二是模糊指数 m 选 1.5、2、2.5 结果差别很大后面我会专门讲。2.2 初始值和局部最优FCM 的迭代更新本质上是一个局部优化方法——它沿着目标函数下降的方向走但目标函数是非凸的存在多个局部极小点。如果初始聚类中心选得不好算法可能收敛到一个局部最优解得到的聚类结果不是一个让人满意的分组。我举个例子。假设数据真实有三类但你随机初始化时有两个聚类中心落到了同一类区域内另一个区域没有中心那么 FCM 大概率会把那两个中心在迭代后依然挨在一起而另一类就被强行拆开跟其他簇混着分。程序不会报错聚类指标看起来也说得过去但画图一看就明显不对。这个问题在 K-means 里也存在但 K-means 因为目标函数更简单多次随机初始化能缓解一部分。FCM 因为模糊隶属度的存在目标函数曲面更平滑但也更复杂对初始值更敏感。传统做法是多次跑、选目标函数最小的一次。这个办法在数据量小、簇数少的时候有效但数据量大、簇数多的时候计算开销成倍增加而且本质上没有跳出“局部最优依赖运气”的困境。2.3 PSO 怎么弥补这个短板粒子群算法PSO是模拟鸟群觅食的群体智能算法。每个粒子代表一个候选解粒子在解空间里飞行通过个体历史最优pbest和群体历史最优gbest不断调整自己的速度和位置。PSO 的优势在于它是全局搜索算法不像 FCM 那样走梯度下降的路线。它在解空间中同时撒下几十个粒子相当于同时从几十个点出发进行搜索粒子之间还能通过信息共享跳出局部区域找到全局更优解。用 PSO 优化 FCM 的具体做法常见有两种思路方式 A用 PSO 搜索最优的初始聚类中心把中心放进 FCM 做最终聚类。方式 B在 PSO 迭代中不断嵌入 FCM 的隶属度更新与聚类中心更新混合迭代。我项目里选的是方式 A。原因很简单方式 A 实现清晰、稳定且能显著改善 FCM 初始值敏感问题方式 B 虽然理论上更“融合”但迭代计算量很大粒子群每更新一代每个粒子都要跑一遍完整 FCM实际项目中跑起来非常慢。把聚类中心编码成粒子的方式是这样的如果数据维度是 D聚成 K 类那每个粒子就由 K×D 个数值拼成一个一维向量代表 K 个聚类中心的坐标。适应度函数就是 FCM 的目标函数值 J粒子越优J 越小。3. PSO-FCM 的 Matlab 实现可直接改造使用3.1 整体运行流程设计我把整个项目拆成四个模块数据预处理模块、PSO 优化模块、FCM 聚类模块、可视化与指标模块。流程是加载数据并归一化得到XN×D。设置聚类数 K、模糊指数 m、PSO 参数种群数、迭代数、惯性权重、学习因子。初始化粒子群每个粒子位置是 K×D 维的随机值。进入 PSO 迭代计算每个粒子的适应度更新 pbest 和 gbest更新速度和位置。迭代结束后把 gbest 解码成 K 个聚类中心作为 FCM 的初始聚类中心。运行 FCM 得到最终的隶属度矩阵和聚类结果。计算轮廓系数、模糊划分系数等评价指标绘制聚类中心曲线和迭代收敛曲线。这里有个经验PSO 优化完之后直接用 gbest 作为 FCM 初始中心再让 FCM 迭代到收敛。这样既利用了 PSO 的全局搜索能力又利用了 FCM 的精局部收敛能力两者互补。3.2 数据预处理与粒子编码数据预处理代码不复杂但每一步都要想清楚为什么。% 读取原始数据每行一个用户每列对应一个时刻点 data load(load_data.mat); X_raw data.load_data; % N x D % 过滤异常值示例功率为负则置NaN再补全 X_raw(X_raw 0) NaN; X_raw fillmissing(X_raw, linear, 2); % 归一化按行用户归一化到 [0,1] X (X_raw - min(X_raw, [], 2)) ./ (max(X_raw, [], 2) - min(X_raw, [], 2) eps);这段代码里eps是防止分母为 0 的。按行归一化的逻辑是把每个用户的负荷曲线都缩放到 [0,1] 区间保留“形状信息”而弱化“容量信息”。如果一个用户空调开得多、一个开得少但开启时段相似归一化后它们会靠得很近更容易被聚到一起。粒子编码部分我直接在初始化时生成随机中心% 如果数据已经归一化到[0,1]聚类中心的取值范围也限定在[0,1] pos rand(pop, K * D); % 每个粒子K个D维聚类中心拼接 vel zeros(pop, K * D); % 速度初值设0这里要强调一点粒子的位置边界一定要与数据归一化范围一致。如果数据是 [0,1]粒子位置也限定在 [0,1]如果数据做了标准化均值为0、方差为1那粒子位置边界要根据实际数据范围设置否则 PSO 会朝着无效区域搜索。3.3 适应度函数与 FCM 迭代的代码实现适应度函数是 PSO 与 FCM 之间的桥梁。它的输入是一组聚类中心输出是 FCM 目标函数值。实现如下function J fcm_fitness(X, centers, m) % X: N x D 数据矩阵 % centers: K x D 聚类中心 % m: 模糊指数 N size(X, 1); K size(centers, 1); % 计算每个样本到每个中心的距离平方欧氏距离 dist zeros(N, K); for k 1:K diff X - repmat(centers(k, :), N, 1); dist(:, k) sum(diff .^ 2, 2); end dist(dist eps) eps; % 避免除零 % 计算隶属度矩阵 inv_dist dist .^ (-1 / (m - 1)); U inv_dist ./ sum(inv_dist, 2); % 计算目标函数值 J sum(sum((U .^ m) .* dist)); end这段代码里最关键的是隶属度矩阵的求法。dist .^ (-1/(m-1))是对距离矩阵逐元素做幂运算然后按行归一化。这样求出来的 U 满足每行和为 1且距离越小的簇隶属度越大。实际测试下来这个写法比双层循环快很多因为 Matlab 的矩阵运算已经优化过了。FCM 主体的迭代代码可以这样写接受 PSO 给出的初始中心并继续精调function [centers, U, obj_hist] fcm_run(X, K, m, init_centers, max_iter) N size(X, 1); centers init_centers; obj_hist zeros(max_iter, 1); for iter 1:max_iter % 第一步更新隶属度矩阵 dist zeros(N, K); for k 1:K diff X - repmat(centers(k, :), N, 1); dist(:, k) sum(diff .^ 2, 2); end dist(dist eps) eps; inv_dist dist .^ (-1 / (m - 1)); U inv_dist ./ sum(inv_dist, 2); % 第二步更新聚类中心 centers_new (U .^ m) * X ./ sum(U .^ m, 1); centers centers_new; % 计算目标函数 obj_hist(iter) sum(sum((U .^ m) .* dist)); end end实际迭代时我会把停止条件设为“目标函数变化量小于 1e-5 或达到最大迭代次数”避免不必要的计算。3.4 粒子群主循环实现PSO 主循环是整个程序的核心。我直接贴出关键过程并注释每行的作用% PSO 参数设置 pop 30; % 粒子数 maxgen 60; % 迭代代数 w_max 0.9; % 惯性权重上限 w_min 0.4; % 惯性权重下限 c1 2.0; % 个体学习因子 c2 2.0; % 群体学习因子 dim K * D; % 粒子维度 % 初始化 pos rand(pop, dim); vel zeros(pop, dim); pbest_pos pos; pbest_fit inf(pop, 1); for i 1:pop centers reshape(pos(i, :), K, D); pbest_fit(i) fcm_fitness(X, centers, m); end [gbest_fit, gbest_idx] min(pbest_fit); gbest_pos pos(gbest_idx, :); best_curve zeros(maxgen, 1); for t 1:maxgen w w_max - (w_max - w_min) * t / maxgen; % 惯性权重线性递减 for i 1:pop r1 rand(1, dim); r2 rand(1, dim); % 速度更新 vel(i, :) w * vel(i, :) c1 * r1 .* (pbest_pos(i, :) - pos(i, :)) ... c2 * r2 .* (gbest_pos - pos(i, :)); % 位置更新 pos(i, :) pos(i, :) vel(i, :); % 边界处理越界后直接截断 pos(i, :) max(0, min(1, pos(i, :))); % 计算适应度 centers reshape(pos(i, :), K, D); fit fcm_fitness(X, centers, m); % 更新个体最优 if fit pbest_fit(i) pbest_fit(i) fit; pbest_pos(i, :) pos(i, :); end % 更新群体最优 if fit gbest_fit gbest_fit fit; gbest_pos pos(i, :); end end best_curve(t) gbest_fit; end惯性权重从 0.9 线性降到 0.4这个操作很重要。迭代前期 w 大粒子飞行速度快全局探索能力强迭代后期 w 小粒子飞行速度慢局部的精细搜索能力更强。这种动态调整能明显提升 PSO 的收敛效果。边界处理我用了最简单的“截断法”位置越界就强制拉回边界。也见过有人用“边界随机重置”或“反弹速度”但实测下来差别不大截断法胜在稳定、不引入额外随机性。3.5 结果可视化从类别曲线到用户标签聚类本身不产生价值真正有价值的是把聚类结果落在业务上。这一步我会做三个可视化第一个是收敛曲线直接画best_curve观察 PSO 迭代过程是否正常收敛。正常的曲线是前期快速下降、后期平缓趋稳。如果曲线到后期还在剧烈波动说明参数没调好或粒子数太少。第二个是类别中心曲线把最终得到的 K 条聚类中心画在同一个图上横轴是时间点纵轴是归一化功率。这条曲线能直观看出每一类的用电形态比如双峰、晚高峰、全天平稳。第三个是单用户标注。用最终的隶属度矩阵 U每个用户取隶属度最大的簇作为它的类别标签然后可以按类别统计用户的行业、地区、容量等信息为后面的需求侧响应或分时电价设计提供基础。% PSO 结果解码 init_centers reshape(gbest_pos, K, D); % FCM 最终聚类 [centers, U, obj_hist] fcm_run(X, K, m, init_centers, 100); % 获取每个用户的类别标签 [~, label] max(U, [], 2); % 绘制聚类中心曲线 figure(Color, w); plot(centers, LineWidth, 2); xlabel(时间点15min/点); ylabel(归一化功率); legend(arrayfun((x) sprintf(类别%d, x), 1:K, UniformOutput, false), Location, best); grid on;关于标签还要额外提醒一点聚类算法输出的是“第 0 类、第 1 类、第 2 类”本身没有业务含义。你需要根据聚类中心曲线的形态给每一类赋予用户画像标签比如“上班族型”“全天平稳型”这个步骤必须在业务人员参与下完成否则聚类结果只是一堆数字。4. 实验设计、聚类数选择与行为画像解读4.1 聚类数K怎么定轮廓系数与业务约束这一步是很多同学容易忽略的。上来就设 K3, K4结果跑完发现类别之间区分度很低或者某一类用户数量只有几十个根本构不成业务规模。我常用的方法先对归一化后的数据跑 FCM固定初始聚类中心不随机直接用多次运行取最优K 从 2 取到 8然后用轮廓系数Silhouette Coefficient和模糊划分系数Partition Coefficient两个指标评估每个 K 的聚类质量。轮廓系数衡量的是样本与自己所在簇的紧密度和与其他簇分离度的差值取值在 -1 到 1 之间越大越好。模糊划分系数公式为PC (1/N) * Σ(i1→K) Σ(j1→N) u_ij^2PC 越接近 1说明聚类结果的模糊程度越低、划分越清晰。我把两个指标综合起来看而不是只看一个。因为轮廓系数偏向硬聚类的评价逻辑PC 则是模糊聚类专属指标两者结合更稳妥。同时还要考虑业务约束比如做分时电价套餐设计的项目里K 取 3-5 最合适太多类别会导致套餐设计复杂、运营成本高。4.2 聚出的类别能做哪些用户画像解读这里我以三分类结果为例展示典型的曲线解读。如果你聚类得到 K3聚类中心归一化曲线通常对应这么三种情况类别标签曲线特征典型场景对应行为画像类别1白天低、18-22点高上班通勤为主上班族型类别2全天平稳、低幅波动老人/全职在家全天稳定型类别3夜间明显抬升夜间用电、充电负荷夜间活动型拿到这些类别之后可以做交叉分析统计每个类别用户的平均月用电量、最大负荷、峰谷差。类别之间如果在这些关键业务指标上差异显著说明聚类结果是有实际意义的可以用于制定不同的电价套餐或需求响应激励方案。我做这个项目时有一类结果特别有意思夜间活动型用户虽然数量只占 15%但平均月用电量超过全体的中位数两倍。进一步查看发现他们中有不少安装了家庭储能设备夜间谷电充电、白天放电。这类用户是需求侧响应的重点对象也是后续可以做优化调度的切入人群。4.3 PSO-FCM 对比 FCM 的效果我在实验里做了三组对比随机初始化 FCM、K-means 初始化的 FCM、PSO-FCM。每组跑 20 次记录目标函数最小值、平均值和标准差。实验结果非常有代表性方法目标函数最优值目标函数平均值标准差随机初始化 FCM112.54131.2714.16K-means 初始化 FCM109.32115.086.42PSO-FCM106.78108.651.74从数据可以看出两点第一PSO-FCM 找到了更小的目标函数值说明在全局寻优能力上确实比随机初始化好。第二20 次运行的标准差从 14.16 降到 1.74说明稳定性大幅提升。做工程的人都懂算法结果只稳定运行一次是不够的生产环境里每天凌晨要跑一次聚类如果结果有一天一个样下游业务根本没法用。当然PSO 的代价是计算时间。随机初始化 FCM 可能几秒钟跑完PSO-FCM 要跑一两分钟。对于居民用户数量在几万到几十万级别的分析这个时间完全可以接受毕竟分析场景通常不是秒级实时。5. 踩坑记录与调参心得5.1 空簇与隶属度退化问题怎么处理我在调试中第一次遇到“空簇”时FCM 跑完有 2 个聚类的中心完全一样等于白白浪费了一个类别。排查之后发现问题出在 PSO 搜索出的聚类中心里有两个中心距离太近FCM 迭代后它们没有分开最后发生了重叠。处理办法有二在 PSO 初始化时加入中心距离约束如果两个中心距离小于阈值重新生成粒子。在 FCM 迭代中检测空簇或无样本隶属度最大的簇直接对空簇中心加一个小的随机扰动。我项目里用的是第一种办法效果直接for i 1:pop while true cand rand(K, D); dist_centers pdist2(cand, cand); dist_centers dist_centers eye(K) * 1e6; % 忽略自身距离 if min(dist_centers(:)) 0.1 pos(i, :) cand(:); break; end end end阈值 0.1 是根据归一化数据的坐标尺度定的如果你的数据范围大需要相应调大。5.2 距离度量和数据归一化的细节FCM 的标准做法是用欧氏距离但居民负荷曲线有个特殊性曲线之间存在“相位偏移”。比如一个用户晚高峰出现在 18 点另一个出现在 19 点欧氏距离会把它们算得很远但从行为模式上讲它们非常相似。遇到这种情况可以考虑用动态时间规整DTW作为距离度量。不过要注意DTW 的计算复杂度比欧氏距离高很多样本量大时会导致 PSO 迭代非常慢。我的折中方案是先对负荷曲线做平滑和峰值对齐预处理再使用欧氏距离。大多数情况下预处理能解决 70% 的相位问题没必要直接上 DTW。数据归一化也有讲究。我前面提到按行归一化到 [0,1] 是一种做法但如果你关心的维度是用户的绝对用电水平按行归一会丢失负载容量信息。这种情况下可以采用 z-score 标准化或者按全局最大负荷归一化。关键是想清楚你的业务目标到底要保留“形状”还是保留“大小”。5.3 参数设置参考表与调节方向把这套流程试验了小几十次我整理了一份参数设置参考表适合居民用电曲线聚类数据维度在 48 到 96 之间的情况参数推荐范围调节方向聚类数 K3 ~ 5结合轮廓系数和业务需求模糊指数 m1.8 ~ 2.2m 越大分类越模糊粒子数 pop20 ~ 50数据量大适当增加迭代代数 maxgen40 ~ 80收敛曲线未平稳则增加惯性权重 w0.4 ~ 0.9 线性递减w 大全局搜索、w 小局部搜索学习因子 c1, c21.5 ~ 2.5一般保持相等即可模糊指数 m 的选择很容易被忽略。m1.5 时聚类结果接近硬聚类隶属度矩阵很尖锐m3 时所有隶属度都往 0.5 靠拢类别信息变得很模糊。我实际试下来m2 是最稳的起点数据噪声明显时调到 2.2 左右结果更平滑。再强调一个经验不要一上来就直接跑完整 PSO-FCM而是先跑几次基础 FCM观察类别中心曲线和各个指标确认聚类数 K 合理之后再上 PSO 优化。这样能省下大量试错时间。我个人的体会是PSO-FCM 这套方案在居民用电行为分析里最大价值不是“比 FCM 好多少”这个结论而是它把聚类的稳定性真正提上来了。不同批次的数据、不同月份的负荷曲线跑出来的聚类结果保持一致业务才能放心地基于类别标签做运营策略。这套流程也不只适用于居民用电工商业用户分时用电特征分析、台区负荷形态归类逻辑都是一样的。如果有同学想在这基础上继续扩展可以试试把 PSO 换成多目标版本把“类内距离最小”和“类间距离最大”同时作为优化目标得到的结果会更丰富一些。
