做居民用电行为分析这事我最初用的是传统FCM聚类。跑了几次实验结果一次一个样——有时候聚出来的几类曲线几乎重叠有时候某一类只剩一两条曲线拿给业务看根本解释不通。后来把粒子群算法PSO引进来给FCM做初始化优化情况才稳定下来。这篇博客就是完整记录我怎么用Matlab把PSO和FCM串起来对居民日负荷曲线做聚类分析的过程包括数据预处理、算法原理、代码实现和一堆实际踩坑的细节。适合正在做电力负荷分类、客户分群或者论文里需要做算法改进与对比的读者参考。1. 传统FCM聚类在居民用电行为分析上的明显短板1.1 FCM到底解决了什么问题它为什么能做用户分群模糊C均值聚类FCMFuzzy C-Means是K-Means的模糊版本。K-Means硬性地把每个样本划到某一类而FCM允许一个样本以不同隶属度属于多个类隶属度之和为1。在居民负荷曲线这个场景里这种软划分特别合适——因为一个家庭的用电行为本来就不可能是纯粹的单一定义白天可能有人在家办公晚上又全家人用电FCM能自然描述这种叠加属性。FCM的核心是下面这个目标函数[ J \sum_{i1}^{n}\sum_{j1}^{k} u_{ij}^{m} |x_i - c_j|^2 ]其中 (x_i) 是第i条日负荷曲线通常是96点或者24点的向量(c_j) 是第j个聚类中心(u_{ij}) 是样本i对类别j的隶属度m是模糊指数一般取2。算法通过交替迭代更新隶属度矩阵和聚类中心不断降低J值。J的物理意义就是所有样本到聚类中心的加权距离总和越小说明类内越紧凑。用电行为分析的典型做法就是把预处理后的日负荷曲线作为样本用FCM把几千户居民分成几类每一类代表一种典型用电模式。比如早上和晚上出现双峰的上班族家庭下午到凌晨持续用电的夜猫子户负荷曲线一直平平稳稳的老人家庭等等。有了这些分群后续的套餐推荐、需求响应、异常用电检测才能往下做。1.2 初始聚类中心敏感和局部最优问题FCM的本质是坐标下降法加梯度下降思想它对初始聚类中心的选择非常敏感。你先随机初始化一次跑出来的类中心是一组形态换个随机种子再初始化结果可能完全不一样。我试过同一份数据普通FCM连续跑20次轮廓系数最低0.21、最高0.53这个波动幅度根本没有工程可用性。问题出在目标函数J是一个非凸函数存在大量局部极小点。FCM的迭代过程是从初始中心出发沿着梯度方向滑到附近的谷底但谷底不一定是全局最优。尤其居民负荷数据有很强的相似性很多曲线形态差别不大类别之间边界模糊FCM稍不留神就会把两类的中心都初始化到同一片区域导致最后聚类结果失衡。1.3 为什么想到用粒子群算法来救场粒子群算法PSO是典型的群智能全局搜索方法。它模拟鸟群觅食行为每个粒子代表解空间里的一个候选解通过自身历史最优pbest和群体历史最优gbest来不断调整飞行速度与位置。PSO对目标函数是否连续可导没有要求也不怕陷入局部极值因为种群多样性保证了它有能力跳出局部坑。于是很自然的思路就出来了把FCM的初始聚类中心看成PSO里的粒子位置把FCM的目标函数J当作PSO的适应度函数。粒子群在全局范围内搜索到一组较好的聚类中心组合后再交给FCM去做精细迭代收敛。这样等于用PSO做全局粗寻优用FCM做局部精搜索两边优势互补。我当时试过直接完全用PSO替代FCM做聚类让每个粒子直接表示所有样本的隶属度矩阵计算量爆炸不说效果也差。而PSO只优化初始中心、FCM负责微调是很多论文里验证过的高性价比方案。2. 居民用电负荷数据怎么处理才能喂给聚类算法2.1 原始数据长什么样怎么变成96点的日负荷曲线大部分居民电表的数据粒度是15分钟一条一天下来是96个采集点形成一条96维的日负荷曲线。如果个别地区是30分钟粒度就是48点。我们这次用的是15分钟粒度所以目标是把每一户每一天整理成一条长度为96的向量。原始数据通常从营销系统导出来是一个三列表用户ID时间戳有功功率kW。我习惯先把数据读入表格然后按用户和时间戳进行透视变成每一行是一个用户一天每一列是96个点的矩阵。% 读取原始数据CSV格式UserID, TimeStamp, Power data readtable(load_data.csv); % 提取日期信息 data.Date datetime(data.TimeStamp, InputFormat, yyyy-MM-dd HH:mm:ss); data.Hour hour(data.Date); data.Minute minute(data.Date); data.PointIndex floor((data.Hour * 60 data.Minute) / 15) 1; % 1~96 % 透视成宽表每行一个用户一天每列一个采样点 loadMatrix unstack(data, Power, PointIndex, ... GroupingVariables, {UserID, Date}, AggregationFunction, mean);这里还有几个坑要提前说明。第一用户编号和数据量都很大时unstack可能会很慢建议先用groupsummary把同一个用户同一天同一个点的重复值先求平均再透视能快不少。第二缺失值处理不能直接删因为缺一个15分钟点整条曲线就断了。我一般用前后点的线性插值如果连续缺失超过2个小时8个点这条曲线直接丢掉因为要么是表计通信故障要么是电表停走插值出来的虚拟曲线会影响聚类真实性。2.2 归一化与特征选择不是把所有维度都扔进去96维数据直接进聚类算法理论上可以但实际效果并不好。一是维度高导致距离计算非常容易受到个别异常点影响二是计算量大粒子群每个粒子都要算96个维度的距离迭代100次谁都受不了。所以我做了两步处理。第一步是选特征。不必用满96个点我用的特征是日用电量一天总kWh、峰段电量占比早上8点到晚上10点电量除以总量、谷段电量占比晚上10点到次日早上8点、最大负荷出现时刻以及负荷率平均负荷除以最大负荷。这几个指标基本抓住了居民用电行为的核心区别。有的论文还会用信息熵或者波动率可加可不加。第二步是归一化。聚类算法依赖距离度量量纲不一致会让某些特征主导结果。比如日用电量是几十kWh负荷率只有0.2这种小数如果不归一化距离计算基本只看用电量。我用的是min-max归一化% dataFeatures: N x 特征数矩阵 minVals min(dataFeatures, [], 1); maxVals max(dataFeatures, [], 1); dataNorm (dataFeatures - minVals) ./ (maxVals - minVals eps);归一化之后所有特征都在0到1范围内。这个细节很重要后面聚类结果的稳定性和可解释性跟它直接挂钩。对于那些负荷分布特别不均匀的数据集比如少数大电量用户建议归一化改成Z-score避免极端值把所有正常用户挤到一个很小的空间里。2.3 数据清洗的经验那些僵尸用户和大坑做数据清洗时我踩过几个很典型的坑这里直接列出来全零曲线长期空置房、或者电表未装智能模块的家庭一天96个点全是0。如果你不做清洗FCM会专门分出一个全零类这一类毫无业务价值。我的处理办法是直接过滤掉全天最大负荷小于额定启动阈值的记录比如最大负荷低于0.1kW的视为无效。异常尖峰个别用户当天有一两个点特别大比如功率达到几十kW大概率是设备启动冲击或者数据跳变。我采用中位数滤波把超过该曲线中位数5倍的点替换为中位数。跨天用电边界如果你算的是日电量要注意跨天用电的归属。冬天晚开空调的用电高峰可能在零点之后单纯按自然日切割会把一段连续行为劈成两半。我后来用曲线相似度聚类的时候对这种问题没那么敏感因为找的是形态组但如果你用日电量做特征最好还是统一采用气温暖区间的连续日数据避开跨季节突变。清洗完成后最好再把数据可视化一下随便抽个几十户画一下日负荷曲线叠加图看看有没有明显异常。这一步虽然不生成算法但能帮你提前发现很多数据质量问题。3. PSO-FCM的Matlab实现从编码到迭代逐段拆解3.1 算法总体流程整个PSO-FCM聚类的流程分五步初始化种群把FCM的K个聚类中心拼接成一个向量作为粒子的位置向量。假设特征维度为D一个粒子的位置就是一个K×D维的向量。更新隶属度固定聚类中心按FCM公式计算每个样本到各类中心的隶属度。更新聚类中心按隶属度加权平均重新计算聚类中心这一步其实就是FCM的一次迭代。计算适应度用FCM的目标函数J作为适应度值J越小表示该粒子代表的聚类中心组合越好。粒子群迭代更新每个粒子的速度和位置重复步骤2-4直到满足最大迭代次数或gbest不再下降。这里有个关键点粒子群迭代的每一轮里我们并不是只做一次FCM更新而是做3-5次FCM内部迭代让这个粒子位置的潜力充分体现出来。如果只做一次适应度函数计算不够精准会干扰PSO的判断。3.2 关键代码适应度函数怎么写适应度函数是整个算法的核心。它接收一个粒子的位置向量把它还原成K个聚类中心然后执行FCM的更新公式返回最小化的目标函数值。function [fitness, centers, membership] fcmFitness(particle, X, K, m, innerIter) % particle: 1行, K*D列, 编码了K个聚类中心 % X: N行D列, 样本数据 % K: 聚类数 % m: 模糊指数, 通常为2 % innerIter: 粒子内部FCM迭代次数 [N, D] size(X); centers reshape(particle(1:K*D), K, D); % 边界保护防止中心出现NaN centers(isnan(centers)) min(X(:)); centers(isinf(centers)) max(X(:)); for iter 1:innerIter % 计算距离矩阵: N x K distMatrix zeros(N, K); for j 1:K diff X - centers(j, :); distMatrix(:, j) sum(diff.^2, 2); end distMatrix max(distMatrix, eps); % 避免除零 % 更新隶属度 invDist distMatrix .^ (-1/(m-1)); membership invDist ./ sum(invDist, 2); % 更新聚类中心 um membership .^ m; centers (um * X) ./ sum(um, 1); end % 计算目标函数 J 0; for i 1:N for j 1:K J J membership(i,j)^m * sum((X(i,:) - centers(j,:)).^2); end end fitness J; end注意这里我用了两层循环算距离效率不高但在Matlab里向量化写法有点绕为了可读性先保留。真实工程里可以用pdist2代替速度会快不少。后面我在加速技巧那一节专门讲。3.3 主程序流程与参数设置这是PSO主算法的框架。我把粒子位置初始化为从样本里随机抽K个点作为中心保证初始解有一定质量而不是全随机。惯性权重从0.9线性衰减到0.4让算法前期探索、后期收敛。% 参数设置 K 4; % 聚类数 N size(X, 1); % 样本数 D size(X, 2); % 特征维度 m 2; % 模糊指数 innerIter 5; % 粒子内部FCM迭代次数 swarmSize 30; % 粒子群规模 maxIter 50; % PSO最大迭代次数 c1 2.0; c2 2.0; % 学习因子 wMax 0.9; wMin 0.4; % 初始化粒子群每个粒子是K*D维向量 particles zeros(swarmSize, K*D); velocities zeros(swarmSize, K*D); for i 1:swarmSize idx randi(N, 1, K); % 随机抽K个样本作为初始中心 particles(i, :) reshape(X(idx, :), 1, []); velocities(i, :) 0.1 * randn(1, K*D); end % 初始化个体最优和全局最优 pbest particles; pbestFitness inf(swarmSize, 1); gbest zeros(1, K*D); gbestFitness inf; for iter 1:maxIter w wMax - (wMax - wMin) * iter / maxIter; for i 1:swarmSize [fitVal, ~, ~] fcmFitness(particles(i, :), X, K, m, innerIter); if fitVal pbestFitness(i) pbestFitness(i) fitVal; pbest(i, :) particles(i, :); end if fitVal gbestFitness gbestFitness fitVal; gbest particles(i, :); end end % 更新速度和位置 for i 1:swarmSize velocities(i, :) w * velocities(i, :) ... c1 * rand(1, K*D) .* (pbest(i, :) - particles(i, :)) ... c2 * rand(1, K*D) .* (gbest - particles(i, :)); particles(i, :) particles(i, :) velocities(i, :); % 边界约束把粒子位置限制在样本空间范围内 minBound min(X, [], 1); maxBound max(X, [], 1); particles(i, :) min(max(particles(i, :), minBound), maxBound); end % 每轮输出方便观察收敛情况 fprintf(Iter %d: gbestFitness %.4f\n, iter, gbestFitness); end % 最后用最优粒子跑一次完整FCM得到最终聚类结果 [~, finalCenters, finalMembership] fcmFitness(gbest, X, K, m, 50); [~, clusterLabel] max(finalMembership, [], 2);这里要解释一下为什么粒子的每一维度用样本的min和max做边界。因为聚类中心本来就是样本空间里的点如果PSO把中心推到了样本范围之外计算距离时会产生很大的空类区域最后聚出来可能只有两三类是有效的其他类别都是远离数据的空中心。边界约束是保证有效性的最低要求。你要是想给点余量可以把minBound和maxBound各扩大10%再夹逼。3.4 边界处理与加速技巧边界处理上面代码里已经体现了一版。但如果你直接用min(max())这种夹逼做法粒子容易被钉死在边界上速度归零丧失多样性。我后来改成了如果粒子越界则随机重置该维度的位置在该中点和边界之间的随机点同时速度取反效果比单纯夹逼好一点。不过设置上要复杂一些对于1152维K4、D288这种已经是够用了。另一个大坑是循环效率。我最初写的是三层循环粒子群循环 样本循环 聚类中心循环。跑一个30个粒子、50次迭代、5000个样本的实验要等一晚上。后来做了三件事用pdist2(X, centers)计算距离矩阵代替手写双层循环速度提升了大概20倍。隶属度更新整个向量化不需要对j循环。用parfor并行计算所有粒子的适应度。PSO每个粒子的适应度计算是相互独立的天然适合并行。但有两点要注意parfor里面不能访问工作区外的变量要用结构体或Table传参且rng的处理要小心不然每次并行结果不是完全一致。我在不追求可复现性的正式实验里用parfor在做稳定性验证时回到串行。4. 聚类数目K怎么定聚类效果怎么评判4.1 用轮廓系数和DB指数选KK的取值没有标准答案得用指标评估。我最常用的是轮廓系数Silhouette Coefficient和Davies-Bouldin指数。轮廓系数计算每个样本的类内平均距离a和最近邻类平均距离b指标是 (b-a)/max(a,b)取值范围-1到1越大越好。对所有样本取平均值就是整个聚类的轮廓系数。Matlab里直接用silhouette函数。DB指数是计算每一类的类内散度与两类中心距离的比值再取所有类对的最大值。DB越小说明类内紧凑、类间距离大。Matlab有现成函数但通常要自己实现个简短脚本。我通常的做法是先让K从2到8跑一遍PSO-FCM记录每个K的轮廓系数和DB值画一条曲线来选。但注意轮廓系数有时候会鼓励你选特别多的小类因为每个小类内都很紧凑。为了避免过度细分我还会结合业务可解释性——如果某个K值下有一类样本占比不足5%那这个K可能就是分得太碎了除非这一小类确实有明确的业务含义比如整夜开空调的极少数用户。4.2 和普通FCM、K-means的对比实验为了确认PSO-FCM不是自我安慰我拿同一套数据跑了普通FCM随机初始化跑20次取最优、K-meanskmeans函数重复20次和PSO-FCM30个粒子50次迭代。整体结果如下算法最优轮廓系数日电量类间标准差多次运行结果波动K-means0.483.1 kWh±0.06普通FCM0.523.5 kWh±0.15PSO-FCM0.614.2 kWh±0.02注意这里日电量类间标准差是我自己定义的一个指标表示各类日电量均值之间的离散程度越大说明分出来的类别在用电量上有明显差异。PSO-FCM在这个指标上高出不少说明它找出的聚类中心能够拉开类别之间的用电水准差距而不是把所有人都混在一堆。还有一个细节普通FCM跑20次取最优和PSO-FCM单跑一次的结果日常差别不大。但普通FCM单次随机初始化的结果可能很差而PSO-FCM即使设置不同的初始随机种子最后结果基本落在同一水平。对工程应用来说稳定本身就是巨大价值因为你不可能每次都跑20次然后挑最好的那只是实验室行为。4.3 三组典型居民用电行为画像聚类的结果最终要能落到业务上。用我们选的K4方案把某区域2500户居民分成了四类我挑其中三个典型类别说明第一类标准通勤家庭大约占38%。日负荷曲线呈现早高峰6点到8点和晚高峰18点到22点两个明显峰谷白天工作时段负荷很低。日用电量中等偏上负荷率低。这类居民对应上班族家庭早晚做饭、洗澡、开灯密集用电。业务上可以推荐峰谷分时电价套餐鼓励他们把部分洗衣、充电等大功率需求挪到谷段。第二类全天平稳型大约占22%。负荷曲线没有明显起伏白天略高晚上略低但降幅不大。结合日期信息看这类家庭工作日和周末的曲线差别很小大概率是白天有老人或全职主妇在家的多人口家庭。他们全天都有基础用电需求对电价调整的响应意愿相对低更适合保底套餐或电量包产品。第三类深夜活跃型大约占15%。晚上11点之后到凌晨2、3点负荷显著上升白天较低下午也有一个小高峰。典型画像包括夜班工作者、喜欢深夜娱乐的年轻人或者家里有燃气壁挂炉这类夜间保温设备。这类用户是需求响应的重要潜力户如果当地执行分时电价他们原本用电时段就在低谷不用调峰就能享受低电价但也意味着他们可能对电费不敏感需要从其他权益角度做运营。第四类占比很少主要是长期高耗能的特殊家庭这里不展开。5. 工程化落地最容易被忽视的几个细节5.1 初始随机种子对结果的影响以及怎么固定Matlab里rand、randi、rng控制随机数生成。PSO-FCM一旦涉及随机初始化两次运行结果就可能不同。科研要做对比实验时必须在脚本开头加rng(2024)这样一句确保每次跑出来的初始种群完全一样。但注意固定随机种子只是保证可复现不保证结果最好。我推荐的做法是正式实验前先不固定种子用不同种子跑5次看结果偏差如果偏差过大说明PSO参数有问题需要增加种群规模或迭代次数如果偏差很小再挑其中某次固定种子作为最终结果。有时候你会发现两次跑出来的gbest适应度接近但聚类标签差异很大特别是两类边界附近的样本归属不稳定。这时候应该看隶属度而不是只看标签隶属度在0.5附近的样本本来就可以归属两边不必强行分到某一类。5.2 距离度量选欧氏距离还是相关系数FCM默认用欧氏距离对幅值敏感。两条曲线形状完全一样但一户用电量是另一户的3倍欧氏距离会非常大导致它们被分到不同类。这在用电行为分析里可能是错的——我们通常更关心曲线的形态而不是绝对数值。一种改法是用Pearson相关系数转化为距离(d1-r)。相关系数对幅值变化不敏感能把早高峰晚高峰双峰形态相似的曲线聚到一起。但缺点是相关系数忽略电量水平建模意义会变弱。我的经验是在特征工程阶段保留日电量这个特征同时把曲线进行幅值归一化再用欧氏距离跑PSO-FCM效果比单独用相关系数好。具体做法是每条曲线先除以自身的日平均负荷变成归一化形态向量再和日电量、负荷率等标量特征拼接起来。这样既能抓形态又不丢电量的绝对信息。5.3 当数据规模上去之后矢量化和并行我一开始用三层循环5000个样本还能忍但跑到几万户就很痛苦。后面把距离计算改成矩阵运算% 计算所有样本到所有中心的距离矩阵 distMatrix pdist2(X, centers); % 更新隶属度 membership distMatrix .^ (-2/(m-1)); membership membership ./ sum(membership, 2);这样就消掉了一个维度上的循环。粒子群里的适应度计算用parfor并行的话每一步耗时几乎不会因为粒子数增多而线性增长。但要注意parfor传参的X是大矩阵每个worker都会保留一份拷贝内存会膨胀。我在32GB内存的机器上用distributed或者把X分割成几个块传给worker可以避免频繁内存交换。另外如果样本数量特别大比如十万户建议先对样本做一次小规模初步聚类或者在特征上做主成分分析降维到5-10维。PSO-FCM计算的复杂性主要在距离矩阵维度降下来之后计算量成倍减少而聚类结果通常并不会有太大损失。5.4 从聚类结果到运营建议的转化聚类本身不产生价值产生价值的是聚类之后的分析与策略。拿前面四类用户来说我最后的分析报告里包含这几项对比每一类的日负荷曲线与当地峰谷时段设置计算该类用户的谷段电量占比把谷段占比很低的用户标为削峰填谷潜力用户。计算每一类的负荷同时率。通勤家庭的晚高峰时间非常集中同时率高如果这一区域有很多这类家庭变压器晚高峰过载风险高可以考虑需求响应引导错峰。把聚类标签和用户基础信息户型面积、家庭人数、是否安装光伏做交叉分析寻找那些特征相似但用电行为异常的用户比如同样面积的房子突然从不做饭变成深夜高负荷可能存在用电设备新增或异常。这些转化工作不是算法层面的但却是实际项目里甲方最关心的东西。写论文的时候聚类结果展示到典型曲线图、类别占比、轮廓系数就足够了但做工程项目一定要往下走一步把聚类结果和为什么有用能做什么对应起来才有说服力。最后分享一个实操经验PSO-FCM调参时不要一上来就追求粒子群迭代300次、种群100个。先用小规模数据几百户快速把K和特征确定下来再放大到全量数据。我经历过在2500户上把PSO迭代次数从50调到200轮廓系数几乎没涨但耗时长了好几倍。对不同K值PSO参数甚至可以微调——K越大粒子编码维度越高搜索空间越大通常需要更大的种群或更多迭代。最实用的办法是设定一个早停条件如果gbest适应度连续10次没有下降超过0.1%就终止迭代能省下大量的无效计算。这个技巧在数据量大时尤其管用。
做居民用电行为分析这事我最初用的是传统FCM聚类。跑了几次实验结果一次一个样——有时候聚出来的几类曲线几乎重叠有时候某一类只剩一两条曲线拿给业务看根本解释不通。后来把粒子群算法PSO引进来给FCM做初始化优化情况才稳定下来。这篇博客就是完整记录我怎么用Matlab把PSO和FCM串起来对居民日负荷曲线做聚类分析的过程包括数据预处理、算法原理、代码实现和一堆实际踩坑的细节。适合正在做电力负荷分类、客户分群或者论文里需要做算法改进与对比的读者参考。1. 传统FCM聚类在居民用电行为分析上的明显短板1.1 FCM到底解决了什么问题它为什么能做用户分群模糊C均值聚类FCMFuzzy C-Means是K-Means的模糊版本。K-Means硬性地把每个样本划到某一类而FCM允许一个样本以不同隶属度属于多个类隶属度之和为1。在居民负荷曲线这个场景里这种软划分特别合适——因为一个家庭的用电行为本来就不可能是纯粹的单一定义白天可能有人在家办公晚上又全家人用电FCM能自然描述这种叠加属性。FCM的核心是下面这个目标函数[ J \sum_{i1}^{n}\sum_{j1}^{k} u_{ij}^{m} |x_i - c_j|^2 ]其中 (x_i) 是第i条日负荷曲线通常是96点或者24点的向量(c_j) 是第j个聚类中心(u_{ij}) 是样本i对类别j的隶属度m是模糊指数一般取2。算法通过交替迭代更新隶属度矩阵和聚类中心不断降低J值。J的物理意义就是所有样本到聚类中心的加权距离总和越小说明类内越紧凑。用电行为分析的典型做法就是把预处理后的日负荷曲线作为样本用FCM把几千户居民分成几类每一类代表一种典型用电模式。比如早上和晚上出现双峰的上班族家庭下午到凌晨持续用电的夜猫子户负荷曲线一直平平稳稳的老人家庭等等。有了这些分群后续的套餐推荐、需求响应、异常用电检测才能往下做。1.2 初始聚类中心敏感和局部最优问题FCM的本质是坐标下降法加梯度下降思想它对初始聚类中心的选择非常敏感。你先随机初始化一次跑出来的类中心是一组形态换个随机种子再初始化结果可能完全不一样。我试过同一份数据普通FCM连续跑20次轮廓系数最低0.21、最高0.53这个波动幅度根本没有工程可用性。问题出在目标函数J是一个非凸函数存在大量局部极小点。FCM的迭代过程是从初始中心出发沿着梯度方向滑到附近的谷底但谷底不一定是全局最优。尤其居民负荷数据有很强的相似性很多曲线形态差别不大类别之间边界模糊FCM稍不留神就会把两类的中心都初始化到同一片区域导致最后聚类结果失衡。1.3 为什么想到用粒子群算法来救场粒子群算法PSO是典型的群智能全局搜索方法。它模拟鸟群觅食行为每个粒子代表解空间里的一个候选解通过自身历史最优pbest和群体历史最优gbest来不断调整飞行速度与位置。PSO对目标函数是否连续可导没有要求也不怕陷入局部极值因为种群多样性保证了它有能力跳出局部坑。于是很自然的思路就出来了把FCM的初始聚类中心看成PSO里的粒子位置把FCM的目标函数J当作PSO的适应度函数。粒子群在全局范围内搜索到一组较好的聚类中心组合后再交给FCM去做精细迭代收敛。这样等于用PSO做全局粗寻优用FCM做局部精搜索两边优势互补。我当时试过直接完全用PSO替代FCM做聚类让每个粒子直接表示所有样本的隶属度矩阵计算量爆炸不说效果也差。而PSO只优化初始中心、FCM负责微调是很多论文里验证过的高性价比方案。2. 居民用电负荷数据怎么处理才能喂给聚类算法2.1 原始数据长什么样怎么变成96点的日负荷曲线大部分居民电表的数据粒度是15分钟一条一天下来是96个采集点形成一条96维的日负荷曲线。如果个别地区是30分钟粒度就是48点。我们这次用的是15分钟粒度所以目标是把每一户每一天整理成一条长度为96的向量。原始数据通常从营销系统导出来是一个三列表用户ID时间戳有功功率kW。我习惯先把数据读入表格然后按用户和时间戳进行透视变成每一行是一个用户一天每一列是96个点的矩阵。% 读取原始数据CSV格式UserID, TimeStamp, Power data readtable(load_data.csv); % 提取日期信息 data.Date datetime(data.TimeStamp, InputFormat, yyyy-MM-dd HH:mm:ss); data.Hour hour(data.Date); data.Minute minute(data.Date); data.PointIndex floor((data.Hour * 60 data.Minute) / 15) 1; % 1~96 % 透视成宽表每行一个用户一天每列一个采样点 loadMatrix unstack(data, Power, PointIndex, ... GroupingVariables, {UserID, Date}, AggregationFunction, mean);这里还有几个坑要提前说明。第一用户编号和数据量都很大时unstack可能会很慢建议先用groupsummary把同一个用户同一天同一个点的重复值先求平均再透视能快不少。第二缺失值处理不能直接删因为缺一个15分钟点整条曲线就断了。我一般用前后点的线性插值如果连续缺失超过2个小时8个点这条曲线直接丢掉因为要么是表计通信故障要么是电表停走插值出来的虚拟曲线会影响聚类真实性。2.2 归一化与特征选择不是把所有维度都扔进去96维数据直接进聚类算法理论上可以但实际效果并不好。一是维度高导致距离计算非常容易受到个别异常点影响二是计算量大粒子群每个粒子都要算96个维度的距离迭代100次谁都受不了。所以我做了两步处理。第一步是选特征。不必用满96个点我用的特征是日用电量一天总kWh、峰段电量占比早上8点到晚上10点电量除以总量、谷段电量占比晚上10点到次日早上8点、最大负荷出现时刻以及负荷率平均负荷除以最大负荷。这几个指标基本抓住了居民用电行为的核心区别。有的论文还会用信息熵或者波动率可加可不加。第二步是归一化。聚类算法依赖距离度量量纲不一致会让某些特征主导结果。比如日用电量是几十kWh负荷率只有0.2这种小数如果不归一化距离计算基本只看用电量。我用的是min-max归一化% dataFeatures: N x 特征数矩阵 minVals min(dataFeatures, [], 1); maxVals max(dataFeatures, [], 1); dataNorm (dataFeatures - minVals) ./ (maxVals - minVals eps);归一化之后所有特征都在0到1范围内。这个细节很重要后面聚类结果的稳定性和可解释性跟它直接挂钩。对于那些负荷分布特别不均匀的数据集比如少数大电量用户建议归一化改成Z-score避免极端值把所有正常用户挤到一个很小的空间里。2.3 数据清洗的经验那些僵尸用户和大坑做数据清洗时我踩过几个很典型的坑这里直接列出来全零曲线长期空置房、或者电表未装智能模块的家庭一天96个点全是0。如果你不做清洗FCM会专门分出一个全零类这一类毫无业务价值。我的处理办法是直接过滤掉全天最大负荷小于额定启动阈值的记录比如最大负荷低于0.1kW的视为无效。异常尖峰个别用户当天有一两个点特别大比如功率达到几十kW大概率是设备启动冲击或者数据跳变。我采用中位数滤波把超过该曲线中位数5倍的点替换为中位数。跨天用电边界如果你算的是日电量要注意跨天用电的归属。冬天晚开空调的用电高峰可能在零点之后单纯按自然日切割会把一段连续行为劈成两半。我后来用曲线相似度聚类的时候对这种问题没那么敏感因为找的是形态组但如果你用日电量做特征最好还是统一采用气温暖区间的连续日数据避开跨季节突变。清洗完成后最好再把数据可视化一下随便抽个几十户画一下日负荷曲线叠加图看看有没有明显异常。这一步虽然不生成算法但能帮你提前发现很多数据质量问题。3. PSO-FCM的Matlab实现从编码到迭代逐段拆解3.1 算法总体流程整个PSO-FCM聚类的流程分五步初始化种群把FCM的K个聚类中心拼接成一个向量作为粒子的位置向量。假设特征维度为D一个粒子的位置就是一个K×D维的向量。更新隶属度固定聚类中心按FCM公式计算每个样本到各类中心的隶属度。更新聚类中心按隶属度加权平均重新计算聚类中心这一步其实就是FCM的一次迭代。计算适应度用FCM的目标函数J作为适应度值J越小表示该粒子代表的聚类中心组合越好。粒子群迭代更新每个粒子的速度和位置重复步骤2-4直到满足最大迭代次数或gbest不再下降。这里有个关键点粒子群迭代的每一轮里我们并不是只做一次FCM更新而是做3-5次FCM内部迭代让这个粒子位置的潜力充分体现出来。如果只做一次适应度函数计算不够精准会干扰PSO的判断。3.2 关键代码适应度函数怎么写适应度函数是整个算法的核心。它接收一个粒子的位置向量把它还原成K个聚类中心然后执行FCM的更新公式返回最小化的目标函数值。function [fitness, centers, membership] fcmFitness(particle, X, K, m, innerIter) % particle: 1行, K*D列, 编码了K个聚类中心 % X: N行D列, 样本数据 % K: 聚类数 % m: 模糊指数, 通常为2 % innerIter: 粒子内部FCM迭代次数 [N, D] size(X); centers reshape(particle(1:K*D), K, D); % 边界保护防止中心出现NaN centers(isnan(centers)) min(X(:)); centers(isinf(centers)) max(X(:)); for iter 1:innerIter % 计算距离矩阵: N x K distMatrix zeros(N, K); for j 1:K diff X - centers(j, :); distMatrix(:, j) sum(diff.^2, 2); end distMatrix max(distMatrix, eps); % 避免除零 % 更新隶属度 invDist distMatrix .^ (-1/(m-1)); membership invDist ./ sum(invDist, 2); % 更新聚类中心 um membership .^ m; centers (um * X) ./ sum(um, 1); end % 计算目标函数 J 0; for i 1:N for j 1:K J J membership(i,j)^m * sum((X(i,:) - centers(j,:)).^2); end end fitness J; end注意这里我用了两层循环算距离效率不高但在Matlab里向量化写法有点绕为了可读性先保留。真实工程里可以用pdist2代替速度会快不少。后面我在加速技巧那一节专门讲。3.3 主程序流程与参数设置这是PSO主算法的框架。我把粒子位置初始化为从样本里随机抽K个点作为中心保证初始解有一定质量而不是全随机。惯性权重从0.9线性衰减到0.4让算法前期探索、后期收敛。% 参数设置 K 4; % 聚类数 N size(X, 1); % 样本数 D size(X, 2); % 特征维度 m 2; % 模糊指数 innerIter 5; % 粒子内部FCM迭代次数 swarmSize 30; % 粒子群规模 maxIter 50; % PSO最大迭代次数 c1 2.0; c2 2.0; % 学习因子 wMax 0.9; wMin 0.4; % 初始化粒子群每个粒子是K*D维向量 particles zeros(swarmSize, K*D); velocities zeros(swarmSize, K*D); for i 1:swarmSize idx randi(N, 1, K); % 随机抽K个样本作为初始中心 particles(i, :) reshape(X(idx, :), 1, []); velocities(i, :) 0.1 * randn(1, K*D); end % 初始化个体最优和全局最优 pbest particles; pbestFitness inf(swarmSize, 1); gbest zeros(1, K*D); gbestFitness inf; for iter 1:maxIter w wMax - (wMax - wMin) * iter / maxIter; for i 1:swarmSize [fitVal, ~, ~] fcmFitness(particles(i, :), X, K, m, innerIter); if fitVal pbestFitness(i) pbestFitness(i) fitVal; pbest(i, :) particles(i, :); end if fitVal gbestFitness gbestFitness fitVal; gbest particles(i, :); end end % 更新速度和位置 for i 1:swarmSize velocities(i, :) w * velocities(i, :) ... c1 * rand(1, K*D) .* (pbest(i, :) - particles(i, :)) ... c2 * rand(1, K*D) .* (gbest - particles(i, :)); particles(i, :) particles(i, :) velocities(i, :); % 边界约束把粒子位置限制在样本空间范围内 minBound min(X, [], 1); maxBound max(X, [], 1); particles(i, :) min(max(particles(i, :), minBound), maxBound); end % 每轮输出方便观察收敛情况 fprintf(Iter %d: gbestFitness %.4f\n, iter, gbestFitness); end % 最后用最优粒子跑一次完整FCM得到最终聚类结果 [~, finalCenters, finalMembership] fcmFitness(gbest, X, K, m, 50); [~, clusterLabel] max(finalMembership, [], 2);这里要解释一下为什么粒子的每一维度用样本的min和max做边界。因为聚类中心本来就是样本空间里的点如果PSO把中心推到了样本范围之外计算距离时会产生很大的空类区域最后聚出来可能只有两三类是有效的其他类别都是远离数据的空中心。边界约束是保证有效性的最低要求。你要是想给点余量可以把minBound和maxBound各扩大10%再夹逼。3.4 边界处理与加速技巧边界处理上面代码里已经体现了一版。但如果你直接用min(max())这种夹逼做法粒子容易被钉死在边界上速度归零丧失多样性。我后来改成了如果粒子越界则随机重置该维度的位置在该中点和边界之间的随机点同时速度取反效果比单纯夹逼好一点。不过设置上要复杂一些对于1152维K4、D288这种已经是够用了。另一个大坑是循环效率。我最初写的是三层循环粒子群循环 样本循环 聚类中心循环。跑一个30个粒子、50次迭代、5000个样本的实验要等一晚上。后来做了三件事用pdist2(X, centers)计算距离矩阵代替手写双层循环速度提升了大概20倍。隶属度更新整个向量化不需要对j循环。用parfor并行计算所有粒子的适应度。PSO每个粒子的适应度计算是相互独立的天然适合并行。但有两点要注意parfor里面不能访问工作区外的变量要用结构体或Table传参且rng的处理要小心不然每次并行结果不是完全一致。我在不追求可复现性的正式实验里用parfor在做稳定性验证时回到串行。4. 聚类数目K怎么定聚类效果怎么评判4.1 用轮廓系数和DB指数选KK的取值没有标准答案得用指标评估。我最常用的是轮廓系数Silhouette Coefficient和Davies-Bouldin指数。轮廓系数计算每个样本的类内平均距离a和最近邻类平均距离b指标是 (b-a)/max(a,b)取值范围-1到1越大越好。对所有样本取平均值就是整个聚类的轮廓系数。Matlab里直接用silhouette函数。DB指数是计算每一类的类内散度与两类中心距离的比值再取所有类对的最大值。DB越小说明类内紧凑、类间距离大。Matlab有现成函数但通常要自己实现个简短脚本。我通常的做法是先让K从2到8跑一遍PSO-FCM记录每个K的轮廓系数和DB值画一条曲线来选。但注意轮廓系数有时候会鼓励你选特别多的小类因为每个小类内都很紧凑。为了避免过度细分我还会结合业务可解释性——如果某个K值下有一类样本占比不足5%那这个K可能就是分得太碎了除非这一小类确实有明确的业务含义比如整夜开空调的极少数用户。4.2 和普通FCM、K-means的对比实验为了确认PSO-FCM不是自我安慰我拿同一套数据跑了普通FCM随机初始化跑20次取最优、K-meanskmeans函数重复20次和PSO-FCM30个粒子50次迭代。整体结果如下算法最优轮廓系数日电量类间标准差多次运行结果波动K-means0.483.1 kWh±0.06普通FCM0.523.5 kWh±0.15PSO-FCM0.614.2 kWh±0.02注意这里日电量类间标准差是我自己定义的一个指标表示各类日电量均值之间的离散程度越大说明分出来的类别在用电量上有明显差异。PSO-FCM在这个指标上高出不少说明它找出的聚类中心能够拉开类别之间的用电水准差距而不是把所有人都混在一堆。还有一个细节普通FCM跑20次取最优和PSO-FCM单跑一次的结果日常差别不大。但普通FCM单次随机初始化的结果可能很差而PSO-FCM即使设置不同的初始随机种子最后结果基本落在同一水平。对工程应用来说稳定本身就是巨大价值因为你不可能每次都跑20次然后挑最好的那只是实验室行为。4.3 三组典型居民用电行为画像聚类的结果最终要能落到业务上。用我们选的K4方案把某区域2500户居民分成了四类我挑其中三个典型类别说明第一类标准通勤家庭大约占38%。日负荷曲线呈现早高峰6点到8点和晚高峰18点到22点两个明显峰谷白天工作时段负荷很低。日用电量中等偏上负荷率低。这类居民对应上班族家庭早晚做饭、洗澡、开灯密集用电。业务上可以推荐峰谷分时电价套餐鼓励他们把部分洗衣、充电等大功率需求挪到谷段。第二类全天平稳型大约占22%。负荷曲线没有明显起伏白天略高晚上略低但降幅不大。结合日期信息看这类家庭工作日和周末的曲线差别很小大概率是白天有老人或全职主妇在家的多人口家庭。他们全天都有基础用电需求对电价调整的响应意愿相对低更适合保底套餐或电量包产品。第三类深夜活跃型大约占15%。晚上11点之后到凌晨2、3点负荷显著上升白天较低下午也有一个小高峰。典型画像包括夜班工作者、喜欢深夜娱乐的年轻人或者家里有燃气壁挂炉这类夜间保温设备。这类用户是需求响应的重要潜力户如果当地执行分时电价他们原本用电时段就在低谷不用调峰就能享受低电价但也意味着他们可能对电费不敏感需要从其他权益角度做运营。第四类占比很少主要是长期高耗能的特殊家庭这里不展开。5. 工程化落地最容易被忽视的几个细节5.1 初始随机种子对结果的影响以及怎么固定Matlab里rand、randi、rng控制随机数生成。PSO-FCM一旦涉及随机初始化两次运行结果就可能不同。科研要做对比实验时必须在脚本开头加rng(2024)这样一句确保每次跑出来的初始种群完全一样。但注意固定随机种子只是保证可复现不保证结果最好。我推荐的做法是正式实验前先不固定种子用不同种子跑5次看结果偏差如果偏差过大说明PSO参数有问题需要增加种群规模或迭代次数如果偏差很小再挑其中某次固定种子作为最终结果。有时候你会发现两次跑出来的gbest适应度接近但聚类标签差异很大特别是两类边界附近的样本归属不稳定。这时候应该看隶属度而不是只看标签隶属度在0.5附近的样本本来就可以归属两边不必强行分到某一类。5.2 距离度量选欧氏距离还是相关系数FCM默认用欧氏距离对幅值敏感。两条曲线形状完全一样但一户用电量是另一户的3倍欧氏距离会非常大导致它们被分到不同类。这在用电行为分析里可能是错的——我们通常更关心曲线的形态而不是绝对数值。一种改法是用Pearson相关系数转化为距离(d1-r)。相关系数对幅值变化不敏感能把早高峰晚高峰双峰形态相似的曲线聚到一起。但缺点是相关系数忽略电量水平建模意义会变弱。我的经验是在特征工程阶段保留日电量这个特征同时把曲线进行幅值归一化再用欧氏距离跑PSO-FCM效果比单独用相关系数好。具体做法是每条曲线先除以自身的日平均负荷变成归一化形态向量再和日电量、负荷率等标量特征拼接起来。这样既能抓形态又不丢电量的绝对信息。5.3 当数据规模上去之后矢量化和并行我一开始用三层循环5000个样本还能忍但跑到几万户就很痛苦。后面把距离计算改成矩阵运算% 计算所有样本到所有中心的距离矩阵 distMatrix pdist2(X, centers); % 更新隶属度 membership distMatrix .^ (-2/(m-1)); membership membership ./ sum(membership, 2);这样就消掉了一个维度上的循环。粒子群里的适应度计算用parfor并行的话每一步耗时几乎不会因为粒子数增多而线性增长。但要注意parfor传参的X是大矩阵每个worker都会保留一份拷贝内存会膨胀。我在32GB内存的机器上用distributed或者把X分割成几个块传给worker可以避免频繁内存交换。另外如果样本数量特别大比如十万户建议先对样本做一次小规模初步聚类或者在特征上做主成分分析降维到5-10维。PSO-FCM计算的复杂性主要在距离矩阵维度降下来之后计算量成倍减少而聚类结果通常并不会有太大损失。5.4 从聚类结果到运营建议的转化聚类本身不产生价值产生价值的是聚类之后的分析与策略。拿前面四类用户来说我最后的分析报告里包含这几项对比每一类的日负荷曲线与当地峰谷时段设置计算该类用户的谷段电量占比把谷段占比很低的用户标为削峰填谷潜力用户。计算每一类的负荷同时率。通勤家庭的晚高峰时间非常集中同时率高如果这一区域有很多这类家庭变压器晚高峰过载风险高可以考虑需求响应引导错峰。把聚类标签和用户基础信息户型面积、家庭人数、是否安装光伏做交叉分析寻找那些特征相似但用电行为异常的用户比如同样面积的房子突然从不做饭变成深夜高负荷可能存在用电设备新增或异常。这些转化工作不是算法层面的但却是实际项目里甲方最关心的东西。写论文的时候聚类结果展示到典型曲线图、类别占比、轮廓系数就足够了但做工程项目一定要往下走一步把聚类结果和为什么有用能做什么对应起来才有说服力。最后分享一个实操经验PSO-FCM调参时不要一上来就追求粒子群迭代300次、种群100个。先用小规模数据几百户快速把K和特征确定下来再放大到全量数据。我经历过在2500户上把PSO迭代次数从50调到200轮廓系数几乎没涨但耗时长了好几倍。对不同K值PSO参数甚至可以微调——K越大粒子编码维度越高搜索空间越大通常需要更大的种群或更多迭代。最实用的办法是设定一个早停条件如果gbest适应度连续10次没有下降超过0.1%就终止迭代能省下大量的无效计算。这个技巧在数据量大时尤其管用。
