我拿到一批居民智能电表的日负荷数据时第一反应就是跑一遍模糊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。比如聚类数c4、特征维度d8那么每一个粒子是一个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参数怎么定不是越大越好粒子群算法的参数选择有很强的经验性我的建议是从一组保守参数起步种群规模nPop30最大迭代次数MaxIter100c1c21.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%左右但耗时翻了三倍。这说明参数不是越大越好关键是找到收敛的阈值。对中小规模负荷聚类数据nPop30到40完全够用。4.2 模糊指数m和聚类数c两个最容易忽略的参数模糊指数m默认取2.0但值得多跑几组对比。m小于1.5时结果接近硬聚类样本容易被强行分到某一个类掩盖了用电行为的模糊性m大于3时隶属度过于扁平所有样本对每个类的隶属度都接近0.25聚类中心之间的差异被稀释类别就不清晰了。我建议在m1.8、2.0、2.2三档之间对比轮廓系数通常2.0附近表现最好。聚类数c的选择是一门独立的学问。常用的辅助判断指标有三个手肘法看目标函数J随c增大的下降拐点轮廓系数取最大值的cDBI指数取最小值的c。三者不一定指向同一个c需要结合业务解释性综合判断。我用轮廓系数筛出来c4最合适因为分到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]且c2、d4那么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轮之内快速下降然后趋于平缓如果到后期还在大幅跳动就得回去检查边界处理和随机种子是不是有问题了。