做居民用电行为分析的人应该都遇到过同一种尴尬手里握着几千户的负荷曲线却不知道该按什么标准给用户分层。“上班族”“夜猫子”“全天候用电”这些词大家都会说真到了算法层面往往第一反应就是K-Means跑了之后发现边界样本分哪类都别扭。后来我换成了FCM模糊C均值聚类让每个用户以不同隶属度归属于多个用电模式效果立刻自然了很多。但FCM有个老毛病——聚类中心初值选不好目标函数很容易掉进局部最优。这个项目里我把粒子群算法PSO和FCM组合到一起用PSO去全局搜一份靠谱的初始聚类中心再用FCM精细迭代最后在Matlab里完整跑通了居民用电行为分析流程。文章会从原理、公式、代码到调参踩坑一条线讲清楚适合正在做电力数据分析、用户画像或者刚接触PSO、聚类的同学做参考。1. 从用电数据到用户画像这个项目到底在做什么1.1 居民用电行为分析的核心问题居民用电行为分析本质上是把“一堆用户”拆成“几类特征鲜明的人”。电网公司拿到的是海量负荷曲线每15分钟或者每1小时一个采样点一户居民一天就是96点或24点一个月下来数据量大到没法直接看。我们不可能逐个用户去问“你家是夜里用电多还是白天用电多”只能靠数据特征去猜。常用的做法是把时序数据压缩成若干业务特征比如日平均负荷、峰段电量占比、谷段电量占比、负荷率、最大负荷出现时刻、夜间用电占比等等。有了这些特征之后才能讨论“哪些用户适合参与需求响应”“哪些用户可能装了电动汽车”“哪些用户是典型的上班族生活方式”。这里有个很实际的问题用户类型之间并不是一刀切。比如一个白天在家办公的人他的负荷曲线跟普通上班族有重叠又跟全天候住宅用户相似一个晚上开空调的用户可能同时具备“夜间高耗能”和“季节性波动”两种属性。如果用K-Means这种硬聚类每个样本只能被硬塞进一个类别结果往往是类别边缘的用户被分得特别勉强。所以我在这个项目里选了FCM让每个用户对每个类别都有一个0到1之间的隶属度允许“白天办公晚上娱乐”这种交叉状态存在分群结果更贴近真实生活。1.2 FCM聚类为什么比K-Means更贴合用电数据FCM的全称是Fuzzy C-Means模糊C均值聚类。它的核心思想很简单一个样本可以同时属于多个类只是隶属度不同。假设有N个用户、C个用电模式FCM要优化的目标函数是J Σ(i1,N) Σ(j1,C) u(i,j)^m * ||x(i) - c(j)||²其中u(i,j)是第i个样本对第j个聚类中心的隶属度m是模糊指数一般取2x(i)是第i个样本的特征向量c(j)是第j个聚类中心。隶属度满足两个条件每个样本对所有类的隶属度之和等于1每个隶属度都在0到1之间。这个目标和K-Means很像区别在于K-Means的u(i,j)只能是0或1而FCM是连续的。用电数据恰恰需要这种连续感。举个例子一个用户的谷段用电占比是52%既不像纯“夜间型”的80%也不像“白天型”的20%用K-Means可能直接被分到某一类而FCM会给出一个比较平均的隶属度比如“夜间型0.45”“全天均衡型0.4”“白天型0.15”。这种输出对业务判断很有价值后续做需求响应时可以按隶属度加权计算潜力而不是只用一个非黑即白的标签。1.3 粒子群算法解决了FCM的哪个死穴FCM虽然好用但它跟K-Means一样是迭代爬山算法。迭代过程从一组初始聚类中心开始通过反复更新隶属度和聚类中心来降低目标函数J问题在于J是一个非凸函数存在很多局部极小点。初始聚类中心选得不好算法很可能收敛到某个局部最优也就是“差不多能用但不够好”的分群结果。尤其当用户特征维数高、数据分布不均匀时FCM对初值非常敏感同一个数据集跑十次可能出来十种不同的效果。粒子群算法Particle Swarm Optimization, PSO是一种群体智能搜索算法模拟鸟群觅食行为。每个粒子代表解空间里的一个候选解通过追踪个体历史最优和群体历史最优来调整自己的速度与位置从而在全局范围内搜索最优解。把PSO和FCM结合以后PSO负责在大范围搜索一组较好的聚类中心FCM负责在局部做精细迭代。这样既保持了FCM的模糊聚类优势又绕开了“初值乱选、结果随缘”的问题。我在项目里实测过同样一份居民用电数据随机初始化FCM的轮廓系数在0.5到0.7之间波动PSO-FCM连续跑五次结果基本稳定在0.7左右这种稳定性对实际项目很重要。1.4 整体技术路线与模块划分整个项目的处理流程我按下面这条线来组织数据读取与清洗从CSV或Excel读入负荷数据剔除抄表异常值处理缺失数据。特征构造与归一化按用户聚合出峰谷占比、负荷率、最大负荷等特征再统一做标准化。PSO全局搜索把C个聚类中心拼成一个粒子位置向量用FCM目标函数作为适应度函数迭代搜索最优聚类中心组合。FCM精细聚类以PSO搜索到的结果作为FCM的初始中心继续迭代到收敛得到最终隶属度矩阵和聚类标签。效果评估与画像计算轮廓系数、Xie-Beni指数绘制各类用户的特征对比图并给每个类命名。这套路线里PSO不是直接替代FCM而是给FCM找一个好起点。听起来简单但真正在Matlab里实现时编码方式、参数设置、数据处理细节都会影响最终效果。下面我把原理和代码拆开讲。2. 粒子群算法与FCM聚类原理先把底层的数学逻辑铺平2.1 PSO速度-位置更新公式与参数含义粒子群算法的核心公式就两个。对于第i个粒子在第k1次迭代时速度更新v(i,k1) w * v(i,k) c1 * r1 * (pbest(i) - x(i,k)) c2 * r2 * (gbest - x(i,k))位置更新x(i,k1) x(i,k) v(i,k1)其中w是惯性权重表示粒子保持先前速度的程度c1是自我学习因子c2是社会学习因子r1和r2是[0,1]之间的随机数pbest(i)是第i个粒子历史最优位置gbest是整个种群历史最优位置。参数怎么理解我习惯用生活类比w像是你滑冰时的惯性惯性大就冲得远、探索性强惯性小就容易停下来细看c1像是“回想自己以前走过的路”c2像是“跟着队伍里最强的人走”。这两个系数如果太大粒子容易在最优解附近震荡如果太小又容易很早聚集到局部区域。在PSO-FCM项目里我把w设为从0.9线性递减到0.5前段侧重全局搜索后段侧重局部收敛c1和c2都设为2.0。这个配置不是拍脑袋而是实测中对FCM适应度函数比较稳的组合。2.2 FCM的模糊隶属度与聚类中心迭代FCM的迭代更新公式也要吃透。先算隶属度u(i,j) 1 / Σ(k1,C) (||x(i)-c(j)|| / ||x(i)-c(k)||)^(2/(m-1))再算聚类中心c(j) Σ(i1,N) u(i,j)^m * x(i) / Σ(i1,N) u(i,j)^m从公式能看出来隶属度取决于样本到每个聚类中心的相对距离。距离某个中心越近隶属度越大。m越大聚类结果越“模糊”也就是隶属度越接近均匀分布m1时FCM就退化成K-Means。通常的经验值是m2用电行为这种边界模糊的数据m在1.8到2.2之间问题都不大。迭代停止条件有两个一是达到最大迭代次数二是目标函数J的变化量小于某个阈值比如1e-6。我在实际代码里两个条件都写上避免数据量大的时候跑太久。2.3 PSO-FCM融合编码与适应度函数设计把两个算法融合在一起第一个要解决的问题是粒子的位置怎么表示聚类中心我的做法很简单假设有C个类、D维特征每个粒子位置就是一个长度为C*D的一维向量前D个元素是第一个聚类中心再D个是第二个聚类中心以此类推。比如C3、D5的时候粒子长度就是15。每次计算适应度时把这个向量reshape成C×D的矩阵当作FCM的初始聚类中心。适应度函数选择就很有讲究了。我在项目里直接用FCM的目标函数J作为适应度因为PSO的搜索方向是让J最小化。但这里有个细节如果在每次适应度评估时都让FCM完整迭代几百次计算成本非常高。粒子数30、迭代100次就意味着要做几千次FCM迭代数据量稍大一点就跑得很慢。所以我做了一个折中在PSO评估阶段只对每个候选中心做10到20次FCM迭代不需要完全收敛只要能对比出不同粒子位置的优劣等PSO结束拿到全局最优位置后再用完整的FCM迭代去做最终聚类。这个技巧能省下不少时间实测效果和完整评估差别不大。2.4 Matlab环境准备、数据读取与特征归一化Matlab环境方面这个项目不依赖太多工具箱。自写FCM迭代函数的情况下只需要基础Matlab就够了如果要用现成的fcm函数则需要Fuzzy Logic Toolbox。我建议自写原因有两个一是fcm函数不允许直接指定初始聚类中心二是自写代码方便打印中间过程、调整模糊指数。数据读取方面我通常把负荷数据整理成一张宽表每一行是一个用户每一列是一个特征。读入代码很简单data readtable(resident_load.csv); disp(head(data))如果原始数据是长表结构也就是一行一个时间点需要先用unstack或者自定义循环把数据重排成“用户×特征”格式。特征归一化是必须做的一步因为峰谷占比可能是0到1的小数而最大负荷可能是几十千瓦的大数如果不归一化欧氏距离会被大数值特征主导聚类结果基本失真。我用zscore做标准化把每个特征变成均值0、方差1的分布。需要注意的是归一化参数只需要在训练数据上计算后续如果有新用户要复用旧的均值和标准差而不是重新算。3. Matlab代码实现从数据读取到聚类结果可视化3.1 特征工程代码与处理细节特征工程这一步直接决定聚类效果的上限。我的项目里用了五个核心特征日平均负荷、峰段电量占比、谷段电量占比、负荷率、夜间电量占比。特征不是越多越好太多会稀释聚类结构太少又区分不出用户类型。下面是一段特征构造的示例代码假设原始数据表raw中包含用户编号、时刻、功率三列% 假设 raw 是长表: UserID, TimeID, Power users unique(raw.UserID); features zeros(length(users), 5); for i 1:length(users) idx raw.UserID users(i); p raw.Power(idx); t raw.TimeID(idx); features(i,1) mean(p); % 日平均负荷 features(i,2) mean(p(t8 t20)); % 峰段电量占比 features(i,3) mean(p(t8 | t20)); % 谷段电量占比 features(i,4) mean(p) / max(p); % 负荷率 features(i,5) mean(p(t22 | t6)); % 夜间电量占比 end这里注意峰谷时段的划分不同地区、不同季节可能不一样要先用业务规则把时段边界定清楚。如果数据里存在个别用户最大负荷为0负荷率就会出现除零问题所以我一般会在前面加一个判断最大负荷小于某个阈值时该用户标记为异常样本直接剔除。缺失值处理上对于单个特征缺失我用该特征的中位数填充因为中位数对离群值不敏感。比均值填充更稳。3.2 PSO优化FCM核心模块主程序的框架长这样clear; clc; close all; rng(42); % 固定随机种子保证实验可复现 data readtable(resident_load.csv); features [data.AvgLoad data.PeakRatio data.ValleyRatio data.LoadRate data.NightRatio]; X zscore(features); C 3; % 聚类数 m 2; % 模糊指数 popsize 30; % 粒子数 maxiter 100; % PSO最大迭代次数 wmax 0.9; wmin 0.5; c1 2; c2 2; [bestCenter, bestVal] pso_fcm(X, C, m, popsize, maxiter, wmax, wmin, c1, c2); [U, centers, J] fcm_iter(X, bestCenter, m, 200, 1e-6); [~, label] max(U, [], 2);PSO函数体我写成了独立子函数方便复用function [gbest, gbestFit] pso_fcm(X, C, m, popsize, maxiter, wmax, wmin, c1, c2) [N, D] size(X); dim C * D; lb repmat(min(X), 1, C); ub repmat(max(X), 1, C); positions rand(popsize, dim) .* (ub - lb) lb; velocities zeros(popsize, dim); fitness inf(popsize, 1); pbest positions; pbestFit fitness; for i 1:popsize centers reshape(positions(i,:), C, D); [~, ~, fitness(i)] fcm_iter(X, centers, m, 15, 1e-3); end [gbestFit, idx] min(fitness); gbest positions(idx,:); for iter 1:maxiter w wmax - (wmax - wmin) * iter / maxiter; for i 1:popsize r1 rand(1, dim); r2 rand(1, dim); velocities(i,:) w * velocities(i,:) ... c1 * r1 .* (pbest(i,:) - positions(i,:)) ... c2 * r2 .* (gbest - positions(i,:)); positions(i,:) positions(i,:) velocities(i,:); positions(i,:) max(min(positions(i,:), ub), lb); centers reshape(positions(i,:), C, D); [~, ~, fitness(i)] fcm_iter(X, centers, m, 15, 1e-3); if fitness(i) pbestFit(i) pbestFit(i) fitness(i); pbest(i,:) positions(i,:); end end [minFit, idx] min(fitness); if minFit gbestFit gbestFit minFit; gbest positions(idx,:); end end end层位置约束那块我特意加了一句把粒子位置限制在数据特征的最小值和最大值之间防止速度太大导致粒子飞到完全没意义的区域。边界处理是PSO工程实现里很容易被忽略的细节不做的话粒子位置可能跑到上千reshape成聚类中心后距离所有样本都非常远适应度值巨大还浪费时间。3.3 FCM迭代细化与标签输出FCM迭代函数我按公式手写没有用工具箱现成函数。代码不复杂但有几个坑必须注意function [U, centers, J] fcm_iter(X, initCenter, m, maxIter, eps) C size(initCenter, 1); centers initCenter; N size(X, 1); U zeros(N, C); J zeros(maxIter, 1); for t 1:maxIter D pdist2(X, centers).^2; D(D 1e-12) 1e-12; % 防止除零 invD D .^ (-1/(m-1)); U invD ./ sum(invD, 2); centers (U.^m) * X ./ sum(U.^m, 1); J(t) sum(sum((U.^m) .* D)); if t 1 abs(J(t) - J(t-1)) eps J J(t); return; end end J J(end); end这段代码里最坑的就是距离矩阵D中如果有零计算D .^ (-1/(m-1))会出现无穷大导致隶属度全部变成NaN。我加了一个很小的下界1e-12确保数值稳定。另一个坑是当某个样本正好落在聚类中心上时正常数学上它对这个类的隶属度应该是1对其他类是0但数值计算需要这个防零处理。聚类结束后max(U, [], 2)得到每个样本归属度最大的类标签。但要注意聚类标签的顺序是不稳定的可能这次运行“类别1”对应的是上次的“类别3”。所以在业务输出时不能直接拿标签编号说事而是要根据每个类的聚类中心特征来命名。3.4 聚类结果的业务画像与可视化聚类跑完之后最重要的是看每一类用户长什么样。我一般先打印聚类中心矩阵再画二维对比图。centers_table array2table(centers, ... VariableNames, {AvgLoad,PeakRatio,ValleyRatio,LoadRate,NightRatio}); disp(centers_table); figure; bar(centers); legend({AvgLoad,PeakRatio,ValleyRatio,LoadRate,NightRatio}, ... Location,northeast); title(各类用户特征均值对比);注意这些聚类中心是在标准化后的空间里计算的直接看数值可以判断相对高低不能直接读“峰段占比是0.6”。最好在预处理时把标准化均值保存下来反归一化后再做业务报告。可视化方式我推荐两种。如果特征是5维以下用bar或parallelcoords足够清晰。如果要把典型日负荷曲线画出来那就要回到原始负荷数据把同一用户各时间点的功率按聚类标签分组再对每组求平均曲线。画出来之后IRL很容易看出“上班族”就是早晚两个峰“夜间型”就是晚上十点后持续走高。给聚类结果命名时不要用数学编号用“上班族”“夜间活力型”“全天平稳型”这种业务语言汇报时大家才听得懂。4. 参数调优、常见报错与避坑经验4.1 聚类个数C怎么定轮廓系数与Xie-Beni指数C值是最难直接拍脑袋定下来的参数。我的经验是分两步走第一步用小范围C分别跑PSO-FCM计算聚类质量指标第二步结合业务判断这个类数是否可解释。常用的指标有两个。轮廓系数用silhouette函数计算值越接近1越好Xie-Beni指标是FCM专用的紧凑性和分离度度量值越小说明簇内越紧凑、簇间越分离。计算流程可以写成for C 2:6 [centerC, ~] pso_fcm(X, C, m, 20, 50, 0.9, 0.5, 2, 2); [Uc, ~, ~] fcm_iter(X, centerC, m, 200, 1e-6); [~, labelC] max(Uc, [], 2); sC mean(silhouette(X, labelC, Euclidean)); disp([C, num2str(C), silhouette, num2str(sC)]); end但指标不是唯一的我遇到过C4时指标分数最高但每类之间的业务差异不够清晰的情况。这时候我会把C4和C5都画出来给业务同事看让他们选更容易解释的方案。算法服务业务而不是反过来。4.2 PSO早熟收敛和结果不稳定的处理办法PSO最容易出现的问题是早熟收敛也就是粒子群在没有找到全局最优时就全部集中到某个局部区域。现象是迭代曲线刚开始下降很快然后很长时间不动最终结果还不如随机初始化FCM跑几次取最优。这里我提供几个亲测有效的办法惯性权重线性递减从0.9到0.5避免后期速度太大来回震荡。限制最大速度把每个粒子的速度上限设为变量范围的10%到20%防止粒子飞得太远。粒子位置越界后重置为边界但速度要适当衰减否则粒子会在边界上反复弹跳。多次运行取最优PSO是随机算法单次结果有运气成分。一般跑5到10次记录全局最优适应度取最好的那个结果。如果做完这些还是不收敛我建议先做一次PCA降维看看特征到底有几个有效维度。有时候数据本身噪声大算法再强也救不回来。4.3 数据预处理阶段最容易踩的坑数据预处理占了这个项目一半以上的工作量也是最容易被低估的部分。我挑几个典型问题说。第一原始负荷数据里的尖峰。有些用户家里有快速大功率设备比如电磁炉、电热水器在采样点正好开机会产生一个特别大的功率值。如果不处理最大负荷特征会被这个单点污染直接带偏整个聚类。我的办法是用分位数截断比如把超过99.5分位数的值替换成99.5分位数。第二季节差异。某用户冬天用电暖器、夏天开空调如果只用某一个月的数据建模分类结果就带着季节性偏差。处理办法要么按月分别建模要么用全年平均特征但要在报告里说明定义。第三特征量纲。前面提过必须归一化这里再强调一次即使做了zscore如果有特征明显偏离正态分布最好先做一次对数变换再标准化。负荷数据通常右偏取对数后结构更清晰。第四缺失值不能一刀切。如果用户某天数据缺失多直接用均值填充会让他的特征缩向平均值聚类时容易变成“墙头草”。我会设置一个阈值缺失率超过20%的用户直接剔除缺失少的用前后时刻插值。4.4 Matlab运行报错定位清单写这个项目时我整理了几个常见的Matlab报错新手遇到不要慌。报错现象常见原因解决办法Error using pdist2输入矩阵维度不是N×D检查X和centers的列数是否一致矩阵维度不一致reshape时C*D与粒子的dim不一致检查生成位置时的dim是否等于C*D输出全是NaN距离矩阵有0或含有NaN特征加D(D1e-12)1e-12清洗输入数据运行时间过长PSO迭代次数太多FCM完整收敛PSO内部用少量FCM迭代只最后跑完整silhouette无法使用缺少统计工具箱换自写轮廓系数或使用Xie-Beni指标每次运行结果不一样随机初始化导致设rng固定随机种子或多轮取最优还有一个容易被忽略的问题Matlab函数文件名必须和函数名一致否则会报“找不到函数”。我把pso_fcm和fcm_iter分别存成pso_fcm.m和fcm_iter.m这是基础设施问题但确实卡过不少人。5. 把这个项目延伸到真实业务小技巧与扩展方向5.1 聚类只是开始业务标签才是抓手很多人在算法跑完后就觉得项目结束了这是最大的误区。聚类输出的标签编号对业务人员来说没有任何意义。我做完聚类后一定会做两件事第一给每类用户起一个业务名字第二统计每类用户的关键指标比如平均日电量、夜间用电占比、负荷率中位数形成一张画像表。只有到了这一步结果才能用于需求响应、台区负荷预测或者精准服务。举个例子。有一次聚类分出四类用户其中一类夜间电量占比极高、平均负荷低。业务同事看到后立刻想到“这可能是有储能设备的用户或者设置了低价谷段加热的热水器用户”。后来核查了一批样本果然这类用户要么家里有蓄热式电热水器要么参与了谷段电价套餐。这就是聚类分析的价值它不是直接告诉你用户家里有什么而是帮你圈定一个值得深挖的候选群体。5.2 从静态聚类走向动态用电行为追踪静态聚类只能反映一段时间的平均状态但用户行为是会变的。今天他是“上班族”可能因为换了工作变成“居家型”。要追踪这种迁移建议把时间窗滑动起来比如每隔一个月重新做一次PSO-FCM然后对比每个用户前后两次的最大隶属度类别是否变化。这样做可以提前识别用户状态改变对电力营销很有用。我在实际项目里的经验是不必每个月都重新跑PSO。第一次聚类时把聚类中心保存下来后面每个月只需要用历史聚类中心作为FCM初始中心拟合新数据再看隶属度变化。这样计算成本很低还能保持聚类标签语义的稳定性。PSO这种全局优化算法不是每次都要重跑只有觉得用户结构可能发生大变化时才重新跑一次完整流程。最后分享一个小技巧在pso_fcm里把每一代的最优适应度值都存下来画一条收敛曲线。这个曲线不只是用来发论文它能很直观地告诉你当前参数下算法到第几代基本稳定。如果迭代到30代就不动了后面70代就是浪费算力如果到100代还在明显下降那说明最大迭代次数设小了或者惯性权重衰减太快。我每次调参都会先看这条曲线比盯着最终指标高效得多。这个项目整体跑下来我的体会是PSO和FCM的搭配并不复杂真正决定成败的还是数据干净程度、特征设计和业务解读这三件事。