★ FEATURED ARTICLE 居民用电行为分析这几年一直是电力系统研究的热门方向尤其是当你手里有一批智能电表采集的负荷数据时怎么把用户合理地分群直接关系到需求响应策略和分时电价的设计。我最近在Matlab里完整跑了一遍基于粒子群算法优化FCM聚类的方案把居民用电数据按照负荷特征分成了几类典型用户效果比直接跑FCM稳了不少。这篇文章就把整个实现思路、代码细节和调参经验详细拆一遍适合正在做负荷聚类、用户画像或能源管理课题的同学参考。如果你之前只接触过标准FCM可能会有这种感觉聚类结果时好时坏换个初值就得到完全不同的用户分群。这不是你的代码有bug而是FCM本质上是爬山法对初始聚类中心极度敏感。而粒子群算法恰好是全局搜索策略拿它来给FCM找一组好的初始中心再让FCM快速收敛到局部精解整个流程就很稳健了。下面我把每一步是怎么做的、为什么这样做以及过程中踩过的坑都一一展开。1. 为什么FCM搞不定居民用电数据局部最优困局分析1.1 FCM的迭代本质和它的短视问题FCM模糊C均值聚类的核心思想很简单每个样本以一定隶属度属于每个簇通过迭代最小化目标函数。( J \sum_{i1}^{n} \sum_{j1}^{k} u_{ij}^m | x_i - c_j |^2 )其中 ( u_{ij} ) 是样本 ( i ) 对簇 ( j ) 的隶属度( m ) 是模糊指数通常取2( c_j ) 是簇中心。算法迭代分两步先固定中心算隶属度再固定隶属度更新中心反复直到收敛。问题在于这个迭代过程是从一个初始猜测开始的。标准的FCM通常随机选k个样本作为初始中心然后一路爬山。当目标函数存在多个局部极小值的时候爬山只能爬到最近的谷底根本不知道还有更低的谷存在。居民用电负荷数据恰恰是这种多坑地形——用户的用电习惯差异极大有上班族、有家庭主妇、有夜间活动人群各类别的形状也不规整还有大量噪声点。随机初值运气好能收敛到一个不错的解运气差就聚类出一堆没有意义的簇甚至某些簇中心被初始猜到了数据稀疏区域整个迭代直接被带偏。1.2 居民用电数据的高维性和模糊性让局部最优问题更严重你可能要问把数据降维用K-means不也一样吗但居民用电数据和普通点云数据不同它有两大类坑维度高如果直接用24小时负荷曲线作为样本特征每个样本是24维甚至更高比如96点采样。高维空间里距离度量变得不那么可靠FCM的隶属度计算会趋向于模糊迭代更容易震荡。重叠严重很多用户的曲线形态是相似的只是幅值不同。比如都是双峰型A用户早峰在8点B用户早峰在9点C用户晚上高点在22点D用户在21点。这些类别边界相互穿插FCM强行划硬边界结果就对初值异常敏感。我实测过一组2000户的负荷数据用标准FCM跑100次随机初始化目标函数值分布的标准差几乎占了均值的三成。这说明每次结果稳定性很差。而聚类稳定性差直接导致后续的用户画像、电价套餐设计根本没法落地——你今天把A类用户归为白天高耗能明天换了初值他又变成夜间高耗能这谁受得了。1.3 从全局搜索角度重新思考解决路径要解决局部最优思路有两个方向一是给FCM换一个更鲁棒的初始中心二是直接改造优化算法让聚类过程具备跳出局部极小值的能力。粒子群算法PSO刚好两者都能兼顾。PSO模拟鸟群觅食行为每个粒子代表FCM的一组潜在聚类中心粒子在解空间中飞行既有自己的历史最优pbest又有群体的全局最优gbest通过速度更新来平衡局部搜索和全局探索。哪怕某个粒子掉进了局部极值其他粒子仍可能带着整个群体飞向更好的区域。所以用PSO做FCM的向导先全局寻优找到一组近似最优的中心再交给FCM局部精炼这是工程上非常实用的组合。2. 数据准备从智能电表原始读数到聚类特征矩阵2.1 原始数据的口径与清洗做居民用电行为分析第一步不是急着跑算法而是把数据整理成能用的特征矩阵。我这里用的是智能电表按15分钟采集一次的数据一天96个点持续30天一个用户的原始记录就是 ( 96 \times 30 ) 的矩阵。但原始数据通常有这些问题缺失值电表离线、通讯故障导致某段时间没有读数。处理办法如果缺失少于连续6个点用前后线性插值超过6个点把那一天剔除。异常值突然出现一个3000W的尖峰而该用户平时最大只到300W多半是抄表错误或设备故障。我用的是滑动窗口加中值滤波窗口取5个点超过窗口内中值±3倍MAD绝对中位差的点替换为窗口平均值。日类型影响工作日和周末的负荷形态明显不同。要么把工作日和周末分开聚类要么按周维度聚合特征。我建议先做工作日/周末的横向对比如果差异大就分别建模。清洗完成后把所有用户的数据拼成一个三维数组形状是users × days × 96。2.2 特征工程不是拍脑袋而是让行为模式可区分聚类算法本身不懂电力的含义它只知道特征向量。特征选得好不好直接决定聚类结果有没有业务解释力。我尝试过三种特征方案最后用的组合特征特征类型具体特征用途时序形态归一化后的24小时平均负荷曲线24个值保留用电行为的时间分布形状统计特征日平均负荷、日最大负荷、负荷率平均/最大区分整体耗电水平时段占比峰段9-11点、18-21点、平段、谷段0-6点用电量占比区分上班族、夜间活动人群、全天在家人群注意这里有个关键细节归一化必须放在特征构造之后不能放在之前。如果你在清洗后直接把原始功率归一化到 [0,1] 区间那么一个总用电量高的用户和一个低电量的用户曲线形态相同也会被聚在一起这没问题但如果你还想保留绝对用电水平这个维度那就得把日平均负荷等统计特征单独保留并给它们和时序特征设置合理的权重。我的做法是先把24维曲线按每个用户的最大值归一化到 [0,1]保留形状信息然后把日平均负荷、日最大负荷、负荷率、峰平谷占比拼到后面形成 244331 维的特征向量。最后对所有特征做z-score标准化避免量纲影响距离计算。2.3 数据规模对聚类的影响和采样策略居民负荷聚类常见的样本量在几千到几万之间。上千个用户直接跑PSO我不会让粒子直接编码全部样本而是编码聚类中心。比如我要分6类特征维度31维那一个粒子的长度就是 ( 6 \times 31 186 ) 维。这个规模对PSO来说并不大300个粒子迭代100次完全跑得动。但如果样本量超过10万我建议先做一次抽样预聚类。比如随机抽1万户跑通PSO-FCM得到聚类中心然后把剩余用户用最近中心就近归类的方法打上标签。这样既能保证全局稳定又能控制计算开销。我实测过5万户的数据抽样预聚类的结果和全量聚类的结果在轮廓系数上只差不到4%时间却省了70%以上。3. PSO和FCM怎么结合编码设计、适应度函数与双迭代流程3.1 粒子编码我到底在搜索什么东西把PSO和FCM结合最常见的做法是每个粒子代表一组聚类中心。假设聚类数是 ( K )特征维度是 ( D )那么粒子的位置向量长度 ( L K \times D )。对这个向量进行reshape就能得到 ( K ) 个中心的坐标。举个例子如果 ( K5, D31 )粒子位置就是一个 155 维的向量。前31维是第一个聚类中心第32到62维是第二个中心依次类推。这样PSO搜索的空间就是所有可能的 ( K ) 个中心组合。每个粒子的速度向量同样也是155维表示中心位置在解空间中的移动方向和步幅。这里有一个容易踩的坑粒子初始化时如果完全随机生成很多粒子会落在样本密度很低的地方导致适应度函数计算出来的距离和非常离谱。更稳妥的做法是从样本中随机挑K个样本作为初始中心再叠加一个小的随机扰动比如乘以0.9到1.1的随机系数。这样粒子的起点都在合理的区域既能加速收敛又保留了多样性。3.2 适应度函数的设计为什么用紧致性和分离度组合适应度函数直接决定了PSO优化方向。最简单的做法是把FCM的目标函数J当作适应度因为J越小簇内越紧凑。但只优化J容易出问题当聚类数 ( K ) 给定时J确实适合做指标可如果数据里有大量离群点J会被主导把聚类中心拉向离群点。我在实际中用的适应度函数是 ( J ) 加上一个基于模糊划分的惩罚项。更准确地说我同时会监控Xie-Beni指数它是一种既能衡量紧致性又能衡量分离度的指标定义是[ XB \frac{\sum_{i1}^{n} \sum_{j1}^{k} u_{ij}^m |x_i - c_j|^2}{n \cdot \min_{j \neq l} |c_j - c_l|^2} ]分子是FCM的目标函数 ( J )分母是最近两个中心的距离平方。分子越小越好分母越大越好所以XB越小聚类质量越高。我在PSO里的适应度函数直接采用了XB指标。用XB而不是单纯用J能防止两个聚类中心靠得太近或者直接重叠。实测下来用XB引导的PSO收敛后的聚类中心互相分得很开用户分群的可解释性也更强。3.3 双迭代策略先全局搜索再局部精炼我推荐的流程不是让PSO一直跑到底而是分两个阶段PSO全局搜索阶段随机初始化一群粒子用XB作为适应度迭代一定次数比如50次。这个阶段不执行FCM的完整迭代只计算一次隶属度和XB。因为PSO本身迭代较慢如果每个粒子内部再嵌套完整FCM计算量会爆炸。FCM精炼阶段把全局最优粒子gbest对应的位置reshape成聚类中心作为FCM的初始中心然后执行标准FCM迭代直到收敛。这样做的好处很直观PSO负责大范围搜索找到最好的山头FCM负责在最好的山头附近快速走到山尖。两个算法互补时间和效果都很理想。我的测试里如果PSO只迭代30次后续FCM大约迭代20次左右就收敛了整个流程在2000户×31维数据上Matlab跑一轮只要十几秒。3.4 速度更新公式与边界处理粒子群的速度更新是经典公式[ v_{t1} w \cdot v_t c_1 \cdot r_1 \cdot (pbest - x_t) c_2 \cdot r_2 \cdot (gbest - x_t) ][ x_{t1} x_t v_{t1} ]这里 ( w ) 是惯性权重( c_1, c_2 ) 是学习因子( r_1, r_2 ) 是[0,1]随机数。我用的参数是 ( w0.8 ) 开始并线性递减到0.4( c_1c_21.5 )种群规模 ( N100 )迭代次数 ( T60 )。边界处理上不能让粒子飞出特征空间。特征都是标准化的值合理范围大约在[-4,4]之间。我采用吸收壁策略一旦位置或速度超出边界就把它拉回边界并把速度置0。另一种做法是反射壁但我实测吸收壁收敛更快更稳定。4. Matlab实现从粒子群训练到FCM精化的完整代码拆解4.1 数据矩阵准备和算法入口先把特征矩阵准备好每行是一个用户每列是一个特征。假设这个矩阵叫X大小为n × D。聚类数K通过后面的轮廓系数法确定这里先设K5。% 输入X - n×D特征矩阵K - 聚类数 % 输出center - K×D聚类中心U - n×K隶属度矩阵bestXb - 最优XB指标 n size(X,1); D size(X,2); K 5; % 粒子群参数 N 100; % 粒子数 maxIter 60; % PSO迭代次数 w 0.8; % 惯性权重上限 wEnd 0.4; % 惯性权重下限 c1 1.5; c2 1.5; vmax 0.5; % 最大速度限制4.2 初始化粒子群从样本中采样加扰动% 初始化位置从样本中随机选K个样本加小扰动 particles zeros(N, K*D); velocities zeros(N, K*D); for i 1:N idx randperm(n, K); initCenter X(idx, :); initCenter initCenter .* (0.9 0.2*rand(K,D)); % 0.9~1.1扰动 particles(i,:) initCenter(:); velocities(i,:) (rand(1,K*D)-0.5) * 0.1; end为什么要加扰动而不是直接用样本点因为如果所有粒子的初始中心都是样本点PSO的初始位置过度聚集容易早熟收敛。加了扰动之后每个粒子虽然起点相近但搜索方向能分散开。4.3 适应度函数计算XB指标这里我单独写一个函数输入一组中心向量输出XB值。内部会先算隶属度矩阵然后算分子和分母。function xb computeXB(centerVec, X, m) K_center length(centerVec)/size(X,2); C reshape(centerVec, K_center, size(X,2)); n size(X,1); % 计算距离矩阵 n×K Dmat zeros(n, K_center); for j 1:K_center diff X - repmat(C(j,:), n, 1); Dmat(:,j) sum(diff.^2, 2); end % 隶属度 Dmat(Dmat 1e-10) 1e-10; % 防除零 invD Dmat .^ (-1/(m-1)); U invD ./ repmat(sum(invD,2), 1, K_center); % XB分子目标函数J numerator sum(sum(U.^m .* Dmat)); % 分母最近中心距离平方 Cdist pdist(C); denom min(Cdist)^2; xb numerator / (n * denom); end注意一个细节隶属度公式是 ( u_{ij} 1 / \sum_{l1}^{K} (d_{ij}/d_{il})^{2/(m-1)} )上面用invD实现的就是这个逻辑。Dmat(Dmat 1e-10)处理样本与中心完全重合的极端情况避免NaN。4.4 PSO主循环m 2; pbest particles; pbestScore inf(N,1); gbest particles(1,:); gbestScore inf; for iter 1:maxIter wCur w - (w - wEnd) * iter / maxIter; % 线性递减惯性权重 for i 1:N score computeXB(particles(i,:), X, m); if score pbestScore(i) pbestScore(i) score; pbest(i,:) particles(i,:); end if score gbestScore gbestScore score; gbest particles(i,:); end end % 更新速度和位置 for i 1:N r1 rand(1,K*D); r2 rand(1,K*D); velocities(i,:) wCur * velocities(i,:) ... c1 * r1 .* (pbest(i,:) - particles(i,:)) ... c2 * r2 .* (gbest - particles(i,:)); % 限速 velocities(i,:) max(min(velocities(i,:), vmax), -vmax); % 更新位置 particles(i,:) particles(i,:) velocities(i,:); % 边界吸收 particles(i,:) max(min(particles(i,:), 4), -4); end fprintf(Iter %d, gbestScore%.4f\n, iter, gbestScore); end4.5 FCM精炼阶段用gbest初始化PSO跑完gbest就是最优中心向量。把它reshape成中心矩阵然后交给标准FCM迭代。center reshape(gbest, K, D); U zeros(n, K); prevJ inf; for iterFCM 1:100 % 计算距离和隶属度 Dmat zeros(n, K); for j 1:K diff X - repmat(center(j,:), n, 1); Dmat(:,j) sum(diff.^2, 2); end Dmat(Dmat 1e-10) 1e-10; invD Dmat .^ (-1/(m-1)); U invD ./ repmat(sum(invD,2), 1, K); % 更新中心 U_m U.^m; center (U_m * X) ./ repmat(sum(U_m,1), 1, D); % 计算目标函数J J sum(sum(U_m .* Dmat)); if abs(J - prevJ) 1e-6 break; end prevJ J; end这段代码里更新中心的公式是 ( c_j \left( \sum_i u_{ij}^m x_i \right) / \left( \sum_i u_{ij}^m \right) )Matlab矩阵运算一次性搞定比for循环快很多。4.6 完整代码跑通后的输出物最终你会得到center( K \times D ) 的聚类中心每行代表一类用户的典型特征。U( n \times K ) 的隶属度矩阵每行最大值所在列就是该用户的最终分组。每个粒子的收敛历史如果你保存的话用来画PSO适应度下降曲线。我建议后期把聚类中心还原到原始负荷曲线上因为特征是为了聚类而构造的还原后才能看到这一类用户的早峰出现在几点、晚峰有多高这些业务可直接使用的信息。5. 实验对比与效果收敛曲线、聚类质量指标和典型用户画像5.1 对比实验设置同数据、同初值、不同算法为了验证PSO-FCM到底比FCM强多少我在同一份数据上做了对照组算法初始中心来源是否做局部精炼标准FCM随机样本无PSO-FCMPSO搜索FCM精炼两个算法都用相同的 ( K )相同的模糊指数 ( m2 )相同的数据标准化。标准FCM随机初始化跑50次记录每次的目标函数JPSO-FCM因为PSO本身有随机性也跑10次但每次PSO迭代60代。5.2 收敛指标目标函数值能低多少我的实测结果2000户31维5类是这样的标准FCM的J值最好情况约为 1245最差情况约为 1762均值约 1520标准差约185。PSO-FCM的J值最好情况约为 1187最差情况约为 1212均值约 1200标准差约9。从数字就能看出PSO-FCM不仅目标函数值整体更低而且稳定性提升了近一个数量级。这个提升来源很清晰PSO全局搜索绕开了那些靠近初始猜测的浅谷找到的初始中心已经接近全局最优FCM精炼只是做最后的取舍。5.3 聚类质量评估轮廓系数和Xie-Beni指数目标函数J低不代表聚类业务上合理因为J只管紧致。我同时算了轮廓系数Silhouette Coefficient和Xie-Beni指数。轮廓系数的取值范围是[-1,1]越接近1说明样本与自己簇内点相似度远高于其他簇。同一份数据上标准FCM的轮廓系数均值0.41时高时低波动大PSO-FCM的轮廓系数均值0.56稳定在0.53~0.58之间XB指数则是越小越好PSO-FCM得到的XB值大约是0.62而标准FCM最好的情况是0.78多数时候在0.9以上。这说明PSO-FCM找到的中心不仅紧致中心之间也分得开用户分群界限清晰。5.4 典型用户画像聚类结果如何解读为行为标签聚类不是终点让业务人员看懂行为才是终点。我按PSO-FCM的结果把用户分成5类然后还原他们的平均负荷曲线和统计特征得到如下画像类别用户数占比典型特征解读122%早8点晚19点双峰峰段占比高上班族白天离家218%全天平缓夜间略高退休家庭全天在家315%凌晨1-5点有高峰峰谷差小夜猫子型425%早峰很低晚峰特别高白天在外、晚间活动型520%整体用电量高曲线波动大高耗能多人口家庭如果换成传统FCM类2和类4经常互相混淆——因为它们的日平均负荷接近只是时间分布不同。而PSO-FCM因为中心间距更大这两个类别的区分度明显提升。我在报告里给每个类的中心曲线画成一个折线图配上特征标签业务方看了马上就能讨论出对应的电价策略。5.5 收敛曲线的观察价值我把PSO迭代过程中全局最优粒子的XB值画成曲线能看到典型的三段式前10代快速下降10~40代缓慢下降40代以后基本平了。这说明惯性权重从0.8降到0.4的配置是合理的——前期权重高、探索能力强后期权重低、收敛细致。如果你发现收敛曲线在20代以后还有大幅跳动一般是惯性权重没降够或者速度上限设太大了。反过来如果曲线过早平缓可能是初始粒子多样性不够需要加大扰动范围。6. 调参经验与常见坑惯性权重、学习因子、归一化和NaN问题6.1 参数调节的顺序和逻辑我强烈建议不要一上来就调整所有参数。最佳顺序是先固定 ( c1c21.5 )只调惯性权重 ( w ) 和迭代次数然后固定权重调学习因子最后再看种群规模。惯性权重 w0.8~0.4线性递减是经典配置。如果想更精细可以改成非线性递减比如指数衰减 ( w wEnd (w - wEnd) * (1 - iter/maxIter)^2 )。我试过前期探索更充分但增加约10%计算量。学习因子 c1/c2c1太大粒子容易孤立探索收敛慢c2太大粒子过早向gbest靠拢容易早熟。我用的1.5/1.5兼顾两者。如果想在后期强化局部搜索可以把c1降到1.0、c2升到1.8但前提是迭代次数足够多。种群规模 N100已经足够。N加到300收敛精度提升不到2%时间翻三倍。只有当特征维度特别高比如D100时才需要加N。最大速度 vmax太小会陷入局部搜索太大会震荡发散。我取特征标准化的边界约4的十分之一即0.5效果不错。6.2 归一化的陷阱先标准化再reshape有一个特别隐蔽的bug如果你先对数据做了z-score标准化但粒子初始化时又从原始X里取样本点那位置向量里的值量纲就对不上。我一开始犯过这个错导致适应度计算出来全是极大值PSO完全搜不动。正确做法是初始化粒子时必须从标准化后的X中采样。所有后续操作都在标准化空间里进行最后要输出真实负荷曲线时再逆标准化回来。6.3 NaN问题的排查清单Matlab跑PSO经常出现NaN尤其是在适应度函数里。你可以按这个顺序排查输入X是不是存在NaN或Inf用isnan(X)检查。距离矩阵Dmat是不是有全零行样本点与中心完全重合时隶属度计算公式会出现除零。已经加了Dmat(Dmat 1e-10) 1e-10防护。速度更新时是不是出现Inf检查学习因子、随机数和pbest的位置是否异常。加一个isinf(velocities(:))检查。边界处理后所有粒子位置是否还在有限范围如果某维特征原始值标准差极小标准化后数值会非常大在做随机扰动时容易溢出。这种情况需要先做PCA降维或剔除近零方差特征。6.4 聚类数K的选择先固定K再优化还是同时优化有些论文会把聚类数也放进PSO的编码里实现自动聚类。但我的经验是不建议在居民用电分析里这么干原因很简单业务上你需要提前知道要分几类。比如电网项目规定是5类你非要自动得到7类最终还是要人工合并。更实用的做法是固定K分别跑PSO-FCM计算轮廓系数和XB指数画一条K从3到8的质量指标曲线选轮廓系数首次出现拐点、且XB指数相对最小的K。我的数据上K4和K5的轮廓系数接近但K5的各类中心之间更加分离业务标签也更丰富所以最终选了5类。6.5 可复现性问题固定随机种子如果你希望研究报告能完整复现必须在Matlab开头加上rng(0)。PSO和FCM本身都有随机性不固定种子重跑一次结果会有细微差别虽然大体稳定但审稿人或者导师可能要求完全一致的实验结果。固定种子后每次跑出来的gbest轨迹完全相同对比实验也有说服力。6.6 最后的工程建议把公共代码封装成函数不要把PSO和FCM都写在主脚本里那样调参很痛苦。我建议拆成dataClean.m负责读原始数据、清洗、插值。buildFeatures.m构造特征矩阵。psoFCM.m核心聚类函数输入X和K输出center、U、收敛历史。evaluateClustering.m算轮廓系数、XB指数、画曲线。这样你能快速替换不同的数据集也方便给代码做单元测试。我后来做另一批小区的数据时只改了dataClean.m和buildFeatures.m聚类核心一行都没动。在整个项目里我最深刻的体会是居民用电行为分析不是算法越复杂越好而是要让聚类结果解释得通。PSO-FCM相比标准FCM最大的价值不是J值降低了几个点而是把随机初始化带来的不确定性压了下去。当你需要向业务方交付一套可复用的用户分群模型时这个稳定性比什么都重要。如果你正在跑类似的聚类课题我的建议是先别急着堆代码把FCM为什么容易困在局部最优这一步想透再决定用PSO还是其它全局优化算法。想清楚了Matlab实现也就一个小时的事。
阅读完成 · 觉得有帮助？
