企业官网建设流程全解析

做用户侧数据分析这几年，有一个问题几乎每周都会被问到：这一万多户居民，到底该怎么分群？看24小时负荷曲线，每一户都不一样，但业务上不可能一户一户去设计方案，必须归纳成几类典型行为。我最后用的方案，是粒子群算法优化Kmeans聚类，Matlab一套脚本跑通，既解决了标准Kmeans对初始中心太敏感的问题，又能把结果拿去做业务解读。这套思路适合刚接触聚类算法、想做用电行为分析的研究生，也适合电力营销、需求响应岗位的从业者直接参考。数据清洗、特征构造、PSO和改进Kmeans的融合方式、代码怎么拆、参数怎么调，我都会讲到，尽量给出一套能直接复现的流程。

1. 为什么居民用电行为分析绕不开聚类

1.1 业务视角：用电行为数据到底能给谁用

智能电表普及以后，电网侧能拿到的数据粒度已经非常细了。常见的是按15分钟或1小时冻结一次电量，一天下来就是96点或24点负荷曲线。数据量一多，难题就来了：这些数字本身不能直接指导业务，必须转化成“哪类用户、什么特征、适合什么策略”的结论。

比如同一个台区里，A用户晚上7点到10点负荷很高，B用户白天9点到17点负荷很高，两户的平均日用电量可能差不多，但它们在分时电价、需求响应、有序用电场景里的价值完全不同。如果用“平均负荷”这样的单一指标去描述，信息的损失太大。聚类就是用来解决这个问题的：把N个用户按负荷曲线特征分成K个簇，每个簇内部用户在用电模式上尽量相似，簇与簇之间尽量不同。

聚类结果一旦出来，用途非常直接。营销侧可以把“晚高峰型”用户单独挑出来，设计峰谷套餐；配网侧可以用不同簇的负荷曲线叠加，估算台区尖峰负荷；需求响应侧则可以针对“可转移负荷占比高”的簇做削峰邀约。所以这不是一个纯算法题，而是一个从数据到业务决策的完整链路。

1.2 算法视角：为什么这类表格型数据普遍用Kmeans

居民用电负荷数据经过特征提取后，形态上就是一个二维矩阵，行是用户，列是特征，比如24个时刻的功率值、峰谷电量占比等。在这种低维度、稠密矩阵上做无监督分群，Kmeans几乎是首选。

原因很直接。第一，Kmeans的时间复杂度接近O(N·K·D)，在几万用户、几十个特征的数据集上跑起来很快；第二，Matlab里内置了kmeans函数，还有silhouette等评价函数，工程实现成本很低；第三，Kmeans输出的簇中心本身就是“平均负荷曲线”，可以直接画出来看，可解释性很强。

但这不代表Kmeans可以直接用。我用标准Kmeans跑同一份居民负荷数据时，连续试了五六次，每次聚类结果都不一样，有时候轮廓系数高一些，有时候低不少。问题就出在Kmeans的迭代机制上：它对初始聚类中心的选择极其敏感，初始中心选得不好，算法就收敛到一个局部较优解，而不是全局理想解。

2. Kmeans的固有缺陷，以及PSO能帮什么忙

2.1 Kmeans迭代求解的本质与三大软肋

先看一下Kmeans在做什么。它本质上在最小化一个目标函数，也就是所有样本到它所属簇中心的距离平方和：

J = sum(i=1..K) sum(x in C_i) || x - u_i ||^2

其中u_i是第i个簇的中心，C_i是第i个簇的样本集合。算法流程很简单：先随机选K个初始中心，然后每个样本就近归类，再重新计算每个簇的均值作为新中心，重复两步直到中心不再明显变化。

这个过程有一个很直观的类比：一个人站在山区里想找最低点，Kmeans只会沿着当前能看到的下坡方向走，走到一个山谷底部就停了，但这个山谷未必是整片山区的最低点。初始中心不同，相当于把这个人放在不同位置，最后掉进的山谷可能完全不同。

具体到居民用电数据，问题集中在三点：

一是初始中心敏感。如果随机抽出来的K个用户恰好都集中在某一种用电模式里，比如全是低电量用户，那么高电量用户所在的簇可能没人去“占领”，最后聚类结果会严重失衡。

二是容易陷入局部最优。Kmeans每一步都在下降目标函数，但每一步都是贪心的局部调整，缺少跳出局部谷底的能力。

三是K值还需要单独确定。这个不是Kmeans本身能解决的，通常靠肘部法则、轮廓系数或业务经验去判断。

2.2 粒子群算法为什么适合做Kmeans的前端优化

粒子群算法（PSO）是模拟鸟群觅食行为的群体智能优化算法。它不要求目标函数可导，也不要求问题有连续解析解，只需要能对每个候选解算出一个适应度值，就能在这个解空间里搜索。

PSO的核心描述很简单：一群粒子在解空间里飞行，每个粒子都有一个位置（代表候选解）和速度（代表移动方向和步长）。每轮迭代，粒子会根据自己的历史最优位置pbest和整个群体的历史最优位置gbest来更新速度，再更新位置。速度更新公式是：

v(t+1) = w·v(t) + c1·r1·(pbest - x(t)) + c2·r2·(gbest - x(t)) x(t+1) = x(t) + v(t+1)

这里w是惯性权重，控制粒子保持原来运动趋势的程度；c1和c2是学习因子，分别控制个体认知和社会认知的权重；r1和r2是两个介于0和1之间的随机数。

把这个机制对比Kmeans的“贪心下山”，PSO的特点是每个粒子不仅仅在局部下降，它还会被“拉向”自己曾经到过的好位置和群体发现的好位置。所以在搜素前期，粒子比较分散，相当于一群人在山区不同位置同时找低洼处；搜素后期，大家逐渐靠拢到目前发现的最低点附近，再在这个区域精细搜索。

这里有一个很关键的思路转换：要把PSO用在Kmeans上，必须把“聚类问题”转换成“优化问题”。聚类中心u1, u2, ..., uK本身就是一组参数，我们把它们拼接成一个向量，就是PSO里一个粒子的位置。比如特征维度是D，簇数是K，那么一个粒子的维度就是K×D。对这个粒子解码，还原成K个簇中心，然后对所有样本计算距离平方和，得到的J值就是适应度。J越小，说明这组簇中心越好。

2.3 两种融合路径：质心初始化 vs 质心精调

实际做实验时，PSO和Kmeans的融合有两条路径可选，概念上容易混，分开说清楚。

第一条路径，是把PSO当作“找初始中心”的工具。先用PSO在K×D维空间里搜索一组比较好的簇中心，然后把这组中心作为Kmeans的初始质心，再让Kmeans继续迭代收敛。这样做的好处是，Kmeans保留了它成熟高效的局部聚类能力，PSO则负责提供一个更合理的起点。

第二条路径，是完全用PSO代替Kmeans的迭代。也就是在PSO每轮解出当前最优粒子后，直接把所有样本按最近距离原则划分到对应簇中，不再调用Kmeans做局部迭代。这种做法理论上更符合“优化”的语义，但在实际数据上往往没必要。因为Kmeans本身在给定较好初始中心后，局部收敛非常快，再让它多走几步并不会破坏PSO的搜索结果，反而能把边界样本的分簇质量修得更细一点。

我在代码里采用的是第一条路径。整体框架是：数据标准化后，用PSO搜索最优质心位置，再把最优质心喂给Matlab内置kmeans函数作为Start参数，最后输出标签和轮廓系数。这样既稳定，又符合大多数论文和业务报告里“PSO优化Kmeans”的表述。

3. 从原始计量数据到可聚类特征：数据准备与预处理

3.1 数据清洗和特征怎么构造

原始数据不能直接进算法。智能电表采回来的数据常见问题有：某天采集失败导致整行缺失、个别时点出现跳变（比如负数、超过表计倍率的极大值）、用户搬家导致一段时间无数据等等。我的处理习惯是：

先按用户ID和日期排序，对单条缺失时点做线性插值；如果某用户一天的完整采集点数不够80%，或者连续多日没有数据，直接剔除这个用户；对超过每日平均负荷3倍以上且持续多个时点的异常值，用前后两天的同一时点值替换。

清洗之后，面临的第二个问题是特征怎么选。最直观的做法是把一天的24小时平均功率作为24个特征。这个方式信息最全，聚出来的簇可以直接画平均日负荷曲线，业务上非常好解释。缺点是特征维度高，PSO粒子维度K×D会变大，计算量上升。

我建议在24维负荷曲线基础上，根据业务目标再补几类关键特征：

这些特征本质上都是对24点曲线的压缩，能降低维度，也方便后期解释。如果原始数据没有峰谷时段定义，可以按本地分时电价时段来切，或者先用固定规则：峰时段取8点到22点，谷时段取22点到次日8点。

3.2 标准化和K值估计

特征构造完成后，标准化这一步不能省。Kmeans基于欧氏距离计算相似度，如果某个特征的量纲特别大，距离就会被它主导。比如日均负荷是几千瓦，峰谷差率是0到1之间的小数，如果不处理，几千瓦的特征在距离计算里会直接把其他特征压死。

Matlab里直接调用zscore函数即可，把每个特征变成均值为0、标准差为1的序列。标准化后的矩阵记为X，每一行一个用户，每一列一个特征，后续聚类和PSO优化都用X。

K值的选择可以用两种方式交叉验证。一种是肘部法则，画出K从2到10时目标函数J的变化曲线，找拐点；另一种是轮廓系数，取轮廓系数最大的K值。但在实际业务里，我一般不会盲信这两个指标，而是把业务约束也放进来。比如电网项目往往希望分3到5类，因为太多类业务执行不过来，太少类又区分不出差异。多数情况下K=4是一个比较合理的起点，后期再根据聚类质量和业务反馈微调。

4. Matlab代码实现：PSO-Kmeans从建模到出图

4.1 主脚本框架与参数设置

下面这段代码是我整理出来的核心主脚本。为了方便讲解，假设已经清洗完数据，得到了一个N×24的负荷矩阵loadData，每一行对应一个用户的24小时平均功率。

% PSO_Kmeans_demo.m % 基于粒子群算法优化Kmeans聚类的居民用电行为分析 % 输入: loadData N×24 原始日负荷曲线 % 输出: idx 聚类标签 % silAvg 平均轮廓系数 rng(1); % 固定随机种子，保证可复现 loadData = fillmissing(loadData, 'linear'); % 线性插值补缺失 X = zscore(loadData); % 标准化，消除量纲影响 N = size(X, 1); D = size(X, 2); K = 4; % 聚类数，可由肘部法则/业务确定 % PSO参数 nP = 30; % 粒子数 maxIter = 60; % 迭代代数 wStart = 0.9; % 惯性权重初始值 wEnd = 0.4; % 惯性权重结束值 c1 = 1.5; % 个体学习因子 c2 = 1.5; % 社会学习因子 dim = K * D; % 粒子维度: K个簇中心

PSO参数不是随便拍的。粒子数一般取20到40，特征是24维、簇数取4的时候，粒子维度是96，30个粒子在96维空间里搜索，覆盖能力已经够用。迭代代数取50到80，太多收益递减，太少容易搜索不充分。w从0.9线性降到0.4是为了实现“前期大范围探索、后期精细收敛”，这是PSO调参里比较通用的做法。

初始化粒子位置时，我不用纯随机数，而是从样本里随机抽K个用户，把这K个用户的特征向量拼接成一个粒子。这样初始质心一定落在真实数据范围内，不会出现一开始就跑到特征空间外部、距离完全失真的情况。

pos = zeros(nP, dim); vel = zeros(nP, dim); pbestPos = zeros(nP, dim); pbestVal = inf(nP, 1); for p = 1:nP idxSample = randperm(N, K); % 随机抽K个用户 initCenter = X(idxSample, :); % K×D pos(p, :) = initCenter(:)'; % 拼接成1×(K*D) end gbestPos = pos(1, :); gbestVal = inf;

4.2 适应度函数与PSO主循环

适应度函数是整个优化过程的核心。给定一个粒子位置，先reshape成K×D的质心矩阵，然后计算所有样本到最近质心的距离平方和。这个值就是PSO要最小化的目标。

function [fit, label] = psoKmeansFitness(X, pos, K, D) % 把粒子位置解码成K个聚类中心，返回距离平方和与最近中心标签 center = reshape(pos, K, D); N = size(X, 1); label = zeros(N, 1); distSum = 0; for i = 1:N dist = sum((X(i, :) - center).^2, 2); % K×1 [minD, minIdx] = min(dist); label(i) = minIdx; distSum = distSum + minD; end fit = distSum; end

这个函数写法比较直观，但注意循环遍历所有样本，如果用户量到几十万，单次适应度计算会偏慢。单机实验几万用户以内问题不大，后续做大规模数据时可以用矩阵化重写，后面我会说到。

PSO主循环的逻辑是：先遍历所有粒子计算适应度，更新每个粒子的个体最优和全局最优；再按速度公式更新粒子的速度和位置；位置更新后要做边界截断，防止质心飞到特征空间外面。

xMin = min(X); % 各特征最小值 xMax = max(X); % 各特征最大值 for iter = 1:maxIter w = wStart - (wStart - wEnd) * iter / maxIter; % 评估当前所有粒子 for p = 1:nP [fit, ~] = psoKmeansFitness(X, pos(p, :), K, D); if fit < pbestVal(p) pbestVal(p) = fit; pbestPos(p, :) = pos(p, :); end if fit < gbestVal gbestVal = fit; gbestPos = pos(p, :); end end % 更新速度和位置 for p = 1:nP r1 = rand(1, dim); r2 = rand(1, dim); vel(p, :) = w * vel(p, :) ... + c1 * r1 .* (pbestPos(p, :) - pos(p, :)) ... + c2 * r2 .* (gbestPos - pos(p, :)); pos(p, :) = pos(p, :) + vel(p, :); % 边界截断 pos(p, :) = max(pos(p, :), repmat(xMin, 1, K)); pos(p, :) = min(pos(p, :), repmat(xMax, 1, K)); end end

迭代结束以后，gbestPos就包含了PSO找到的最优质心。注意optimum位置可能落在特征空间的边界上，如果经常贴边，通常说明K值或者特征构造有问题，值得回头查数据。

4.3 把最优中心交给Kmeans并可视化

得到gbestPos之后，我习惯再交给kmeans函数跑一次。这样做的原因是，PSO搜索到的是“质心放在哪里整体距离更小”，而Kmeans在给定这个起点后，还能通过局部迭代把质心调整到更精细的位置。

bestCenter = reshape(gbestPos, K, D); % 把PSO找到的中心作为Kmeans初始质心 idx = kmeans(X, K, 'Start', bestCenter, ... 'MaxIter', 300, ... 'Replicates', 1, ... 'Distance', 'sqeuclidean'); % 轮廓系数 sil = silhouette(X, idx); silAvg = mean(sil); fprintf('平均轮廓系数: %.4f\n', silAvg);

这一步有一个容易踩的坑：kmeans函数里的'Start'参数在传入矩阵时，矩阵每一行就是一个初始中心，所以传入的必须是K×D矩阵，而不是拼接后的向量。因此必须先reshape。

结果可视化我一般分两张图。第一张是聚类标签下的PCA降维散点图，适合快速看簇的重叠程度；第二张是每一类用户的平均日负荷曲线，适合业务解读。

% 图1: PCA降维后看聚类分布 [coeff, score] = pca(X); figure; gscatter(score(:,1), score(:,2), idx); xlabel('PC1'); ylabel('PC2'); title('PCA降维后的聚类分布'); % 图2: 每类的平均日负荷曲线 figure; for k = 1:K subplot(2, 2, k); plot(mean(loadData(idx == k, :), 1), 'LineWidth', 1.5); xlabel('小时'); ylabel('平均功率/kW'); title(['第' num2str(k) '类用户平均日负荷曲线']); grid on; end

issueloadData还在前面定义过，所以在subplot里能用。如果原始数据不是24点而是96点，把横轴刻度改成15分钟间隔即可。

5. 实验效果怎么看：稳定性、轮廓系数与业务解读

5.1 评价一个聚类结果，不能只看目标函数

我试过很多次，单独看目标函数J会有一个问题：PSO把J降得很低，但聚类结果未必是业务上最舒服的。因为J本质上衡量的是“紧凑程度”，没有充分考虑到簇之间的分离程度。所以做效果对比时，我至少看三个角度。

第一个是稳定性。同样一份数据，把随机种子换几个，跑多次，看簇中心和标签是否变化。标准Kmeans如果随机初始化几次，结果经常忽大忽小，有时甚至某个簇只有几个人。PSO-Kmeans在这点上改善非常明显，因为PSO的搜索不是纯随机撒点，而是有记忆、有方向的搜索，最终质心落点更稳定。

第二个是轮廓系数。Matlab里直接算，值域在-1到1之间，越接近1说明样本离自己簇的中心越近、离邻近簇越远。我个人经验是，居民用电行为数据能做到0.2到0.35就算不错，毕竟用户行为本身是连续过渡的，没有绝对清晰的边界。如果有人报告说轮廓系数0.8，大概率是特征或者K值选择有问题，比如把完全重复的数据当成了特征。

第三个是运行时间。PSO-Kmeans增加了适应度计算这一步，时间开销肯定比单次Kmeans大，但一般不会失控。拿1万户、24维特征、30个粒子、60代迭代来说，跑一趟通常在几十秒到两三分钟，取决于机器和循环写法。如果数据到了几十万用户，就需要考虑矩阵化加速或者特征降维。

我自己的经验是，把普通Kmeans（单次）、Kmeans重复10次取最优、PSO-Kmeans放在同一份数据上对比时，差距主要体现在波动性上。普通Kmeans不同随机种子跑出来的轮廓系数可能从0.18到0.26之间乱跳；重复10次取最优会稳定在较高值，但用户耗时要乘10；PSO-Kmeans的轮廓系数我一般能稳定在0.25左右，而且耗时通常低于重复10次Kmeans的总时间。

5.2 聚类结果的行为画像怎么解读

聚类算法只是中间产品，真正交付给业务的是“每类用户是谁、有什么特征、该怎么运营”。对居民用户来说，K=4是比较常见的分法，我在实际项目里经常会看到下面几类画像。

第一类是晚高峰主导型。平均负荷曲线在18点到22点有明显高峰，白天低谷偏低。这类用户大多是上班族，傍晚回家后集中用电，包括做饭、洗浴、娱乐负荷。峰谷差率通常比较大，是分时电价和削峰填谷的重点目标。

第二类是全天平稳型。曲线整体偏高且平坦，负载率很高，夜间负荷也不低。这类用户可能是家里有老人小孩常驻，或者有持续运行的设备。它们用电量基数大，参与需求响应的空间也大，但响应意愿不一定高，需要考虑不影响正常生活为前提的互动策略。

第三类是白天工作型。工作日白天负荷明显，周末和夜间下降。这类用户很多是家庭作坊、小店、或在家办公人群。它们对电价比较敏感，更可能响应白天的激励政策。

第四类是夜间用电型。曲线在22点到次日6点之间偏高。常见原因包括电动汽车充电、蓄热式电暖器等。这类用户在谷电阶段用能比例高，适合谷电套餐，也是电网消纳新能源、填谷的主力。

分完类以后，我一般会把每类用户占总用户数的比例、户均电量、峰谷差率做一张透视表，再结合空间位置画到地图上，这样业务部门可以直接拿去用。

6. 调参、踩坑和可复现性建议

6.1 PSO参数调参经验

PSO的参数不算多，但每个参数对结果都有影响，分享一下我的调参手感。

粒子数nP：特征维度越高，粒子数要多一点。24维特征、K=4时，30个粒子够用；如果特征升到40维以上，建议加到50。粒子数再多，时间开销增长明显，但精度收益会边际递减。

惯性权重w：我习惯用0.9到0.4线性递减。开始阶段w大，粒子速度快，负责全局搜索；后期w小，粒子速度慢，负责局部精细搜索。固定w也可以，但容易在前期收敛太快或后期震荡，所以动态w通常更稳。

学习因子c1、c2：一般取1.5到2.0，并且让c1和c2相等。c2如果明显大于c1，会让所有粒子过早被gbest吸引，种群多样性下降，容易早熟；c1过大则粒子各自为政，收敛变慢。遇到结果发散的情况，先把c1、c2都调到1.5，再观察。

边界处理：粒子位置越界后，简单粗暴地截断到边界。注意这只是硬约束，更平滑的做法是让粒子在越界方向上的速度分量归零，但代码复杂度会上去。样本量不大时，直接截断问题不大。

6.2 Matlab实现中的常见坑

第一个坑是维度拼接顺序。initCenter(:)'是按列取元素，reshape成K×D时也是按列填充。如果中间不小心用了initCenter(:)或者转置，质心矩阵的数据布局会完全错乱，聚类结果自然不对。调试时可以用一个很小的测试矩阵跑一遍适应度函数，打印reshape后的center确认一下。

第二个坑是标准化对象搞错。特征工程阶段，标准化后的X进聚类；但画平均日负荷曲线时要用原始的loadData。如果把标准化的曲线拿去做均值，画出来的形状虽然一样，但Y轴含义变成了“标准差倍数”，业务人员根本看不懂。

第三个坑是K值选错。PSO可以把初始质心优化得很好，但如果K本身不符合数据结构特点，结果依然没有业务价值。比如把K设成2，所有用户被硬分成“高用电”和“低用电”，虽然看起来轮廓系数还行，但行为细节全丢；K设成7，会出现两个空簇或者两个几乎重叠的簇。建议先在标准化后的数据上跑一遍肘部法则，再结合业务确定K。

第四个坑是随机种子。PSO、Kmeans、PCA里都有随机成分，科研和报告里只要强调复现性，就一定要在脚本开头加rng(1)之类的固定种子。不然今天跑出来四类，明天跑出来四类但用户编号完全变了，后面做案例分析会很被动。

第五个坑是样本量变大后的性能问题。前面那个适应度函数为了好懂用了for循环，1万用户没问题，但如果是10万用户，每次适应度计算就要做10万次距离计算，60代乘以30个粒子，总计算压力就上来了。这时候建议重写适应度函数，把距离计算改成矩阵运算：计算N×K距离矩阵，再对行求最小值和索引，速度能快一个数量级。更激进的做法是先PCA降到5到8维再做PSO，粒子维度一下子降下来，收敛速度会快很多。

6.3 从K=4开始跑通第一版

最后给一个实操建议：第一次跑的时候，不要急着做复杂的特征工程。直接拿24点负荷曲线标准化，K固定为4，PSO参数用上面代码里的默认值，先看聚类结果和平均日负荷曲线是否合理。如果分出来的某类曲线和业务常识完全对不上，先查数据清洗，再查K值，然后才是调PSO参数。

我曾经在真实数据上遇到过一个典型问题：有一类用户的平均曲线在凌晨5点突然蹿高，后来查数据发现是热水器定时加热，这类用户聚在一起其实是“凌晨加热型”，并不是数据错误。这就是聚类的价值——它会把你不曾预期但有真实业务逻辑的规律挖出来。

这套流程跑顺之后，后面要做的扩展无非是三件事：一是把适应度函数向量化，支撑更大规模的用户量；二是尝试用轮廓系数、DBI指标代替距离平方和作为适应度函数，看哪种目标更符合业务预期；三是把聚类结果与用户档案、电价方案做关联分析，让“分群”真正变成“运营动作”。在居民用电行为分析这个方向上，算法从来不是终点，把聚类结果用起来才是终点。