简介面向电力负荷数据分析与挖掘的MATLAB实现包针对KMeans算法在用户用电行为聚类中的应用涵盖数据清洗、标准化、聚类执行及典型用电模式提取等环节。包内共14个文件以13个MATLAB脚本(.m)为主另含1个MATLAB数据文件(.mat)压缩包整体大小9.81MB。脚本除实现KMeans核心流程外还提供DBSCAN密度聚类、DBI、轮廓系数等聚类有效性评估代码便于对比不同聚类效果并确定合理类别数k。目前已有328人学习/浏览。通过该资源可快速搭建负荷聚类实验环境直接运行脚本查看用电数据分组结果理解从原始负荷序列到典型模式形成的完整链路各类脚本注释清晰、模块划分明确适合电力系统研究人员、数据分析初学者及课程设计使用者参考。1. kmean聚类在电力负荷分析里的活儿就是把用电数据按曲线形态归堆kmean聚类在电力负荷分析里的活儿说白了就是把海量用电数据按曲线形态归堆。供电公司运营监控平台里沉淀着大量96点负荷数据——每个计量点每天按15分钟一个点存下96个有功功率值。业务方想在几十万甚至上百万条历史记录里抽出“这栋楼工作日的典型用电形态”“这组台区夏季空调负荷典型日曲线”用于需求响应预案、配变容量校核和新用户接入方案比选。数据量大、形态杂、标签几乎为零靠人工枚举早晚高峰或平段特征很难兜住所有形态。K-means只需给定一个k就能把几万条日曲线按形状相近程度归成簇再从每簇中提取最接近质心的实际样本作为典型数据。这条链路从数据整理到结果落库一天内能跑通下文按完整流程逐步展开。2. K-means聚类原理与负荷样本的构造方式2.1 K-means为什么适合做负荷曲线聚类K-means本质上解决的是几何问题而非物理问题给定n条96维负荷曲线先初始化k个随机质心然后交替执行两步——把每条曲线分到距离最近的质心再用簇内全部样本的逐点均值更新质心反复迭代直到质心不再移动。算法收敛时极小化的是簇内平方误差和WCSS也就是sklearn里直接暴露的inertia_指标。它不关心一条曲线早上有没有峰只关心欧氏距离意义上的形状接近程度这种“不讲业务故事”的特性在负荷聚类选型时反而是最省心的。工程上优先选K-means而不是层次聚类或DBSCAN有三个现实原因。第一是规模日常一张负荷宽表少则十万行多则百万行K-means单次迭代复杂度约O(n×k×d)在sklearn的OpenMP并行下收敛很快层次聚类仅距离矩阵计算就是O(n²)样本超过五万条就很痛苦树状图几乎无法阅读。第二是簇形态排除表计故障和极端事件后日负荷曲线的类别边界相对清晰每个簇大体呈凸分布没有严重违背K-means对簇形状的假设。第三是产出物K-means的质心直接就是一条96点平均曲线天然符合“典型数据”的定义DBSCAN不给质心层次聚类还要自己再取均值。如果任务只是把上万条日曲线分成若干个有业务含义的类K-means是综合成本最低的起点。DBSCAN在高维曲线样本上需要调min_samples和eps对96维空间的密度定义极不直观调参成本远高于只调一个k。层次聚类虽然能输出可解释的树状图但合并过程不可逆个别异常日会一步步带偏整棵树的形态数据量上来后计算资源也扛不住。只有在业务目标变成“专门抓离群的高耗能时段”时我才会改用孤立森林或DBSCAN。聚类目标不同模型选型结论就不同没必要在常规典型数据提取场景里把工具搞复杂。2.2 负荷样本构造一行是一天一列是一个时刻点把原始电表冻结数据变成聚类矩阵是整个流程里最需要留神的一步。常见做法是每个计量点每天一行按15分钟颗粒度展开成96个列列名依次是p_00到p_95。下面这张宽表结构是后面所有处理的基础。字段说明示例值meter_id计量点编号M001234data_date数据日期2024-06-15p_00 ~ p_95当日96个15分钟冻结点功率(kW)12.5, 12.8, ...构造宽表时最容易出的问题是点序号错位。冻结数据时间标识从0:00开始、23:45结束共96点如果透视时不先排序或者把前一天23:45的点拼进来曲线会整体平移一个点位聚类结果里会出现一批看起来很陡的“伪尖峰簇”本质上只是数据质量问题。所以在pivot之前必须做点序校验确保每个计量点每天的点数完整。import pandas as pd # 原始表结构meter_id, ts(时间戳), power(kW) raw pd.read_parquet(load_15min.parquet) raw[data_date] raw[ts].dt.normalize() raw[point_no] raw[ts].dt.hour * 4 raw[ts].dt.minute // 15 # 每个计量点每天必须且仅有96个点过滤不完整日期 counts raw.groupby([meter_id, data_date]).size() valid_pairs counts[counts 96].index raw raw.set_index([meter_id, data_date]).loc[valid_pairs].reset_index() # 按顺序透视成宽表一行即某计量点某天的完整曲线 raw raw.sort_values([meter_id, data_date, point_no]) wide raw.pivot_table(index[meter_id, data_date], columnspoint_no, valuespower).reset_index() wide wide.rename(columns{c: fp_{c:02d} for c in range(96)})先算point_no再统计点数可以把采集系统重复补招造成的一天出现192个点这类问题暴露出来。point_no用整点编号乘4再加刻钟偏移适配15分钟颗粒度如果源表是5分钟颗粒度公式改成ts.dt.hour * 12 ts.dt.minute // 5列数变成288。透视后的宽表保存建议用parquet不用csv96列浮点数据用csv读写慢一倍还可能丢精度。后续所有清洗和聚类操作都以这个宽表为准。2.3 归一化用形状还是用数值直接决定典型曲线含义K-means默认用欧氏距离96个维度的量纲都是kW但不同计量点的负荷绝对值差异巨大。一台主变峰值8000kW一个商铺峰值只有200kW不归一化的话距离主要由数值大的样本主导小负荷曲线无论形态如何都会被揉进同一个簇。做这一步前必须先决定聚类关心的是“形状”还是“形状加量级”。通常有两套做法。第一套按每条曲线自身最大值归一化得到0到1的比值曲线只保留形状特征适合产业结构分析和典型日形态提取峰值为0的样本要加极小值兜底否则除零后整行变成NaN。第二套用z-score逐列标准化保留曲线间量级差异适合区域负荷叠加和配变容量规划。下面把两套做法封装成函数并同时返回还原系数为后面把典型曲线还原成实际功率做准备。import numpy as np from sklearn.preprocessing import StandardScaler from joblib import dump feat_cols [fp_{i:02d} for i in range(96)] def build_feature_matrix(wide_df, modeshape): x wide_df[feat_cols].fillna(0).values if mode shape: denom x.max(axis1, keepdimsTrue) denom[denom 0] 1e-6 return x / denom, {mode: shape, denom: denom.squeeze()} if mode zscore: scaler StandardScaler() x_norm scaler.fit_transform(x) return x_norm, {mode: zscore, mean: scaler.mean_, scale: scaler.scale_} x_norm, coef build_feature_matrix(wide, modeshape) dump(coef, norm_coef.joblib)denom保存每条样本当天的峰值还原时把归一化曲线乘回峰值即可。z-score模式返回mean和scale各96个值还原公式是逐点乘scale再加mean这两个向量必须和模型及聚类标签放在同一个归档目录。实际项目中最常见的问题就是只存了聚类模型忘了存归一化系数等要还原时再回原始表反查既慢又容易串计量点。另外注意fillna(0)只是兜底真正的大段缺失应该在清洗阶段处理掉不能在特征矩阵里靠填零蒙混过关。3. 用电数据清洗与样本特征构造先洗掉坏数据再谈聚类3.1 三类必坏数据缺失、跳变、与容量矛盾聚类算法对离群样本很敏感几条坏曲线就可能同时拉偏多个簇的质心。清洗阶段我按三类规则依次处理不混在一起做方便后面定位是哪种原因导致淘汰率异常高。第一类是缺失。前面按96点完整性筛过整行但还存在个别点位为NaN或0的情况前后点都有值。处理策略是缺失段不超过4个点且两端有连续有效值用线性插值缺失或连续全零超过8个小时即32个点整个样本丢弃。零值要区分冻结失败和真实空载规则上以“连续8个小时全零”作为无效日判断短时零值照常保留。第二类是跳变。有功功率在极短时间内出现远超正常波动的脉冲值多数来自冻结失败或抄表异常。常见做法是计算每条样本96点的相邻差分把差分绝对值超过该条曲线最大值30%的点记为跳变点一天内跳变点多于3个则整行剔除。跳变和真实冲击负荷如电焊机、轧机的区别在于持续性真实冲击会在一段时间内反复波动而采集跳变通常只是孤立一两个点。第三类是与容量矛盾。计量点档案里有运行容量负荷值超过容量1.5倍或出现明显负值需要标记处理。但负功率在分布式光伏场景里可能是真实的倒送电不能一刀切删掉要看用户档案里是否登记了光伏。只有没有分布式电源合同却持续出现负值时才判定为表计或接线问题。import numpy as np def max_zero_run(row): run best 0 for v in row: if v 0: run 1 best max(best, run) else: run 0 return best def clean_wide(wide, capacity_map, max_ratio1.5): out wide.copy() vals out[feat_cols].values # 连续零值超过32个点(8小时)视为无效日 zero_run np.apply_along_axis(max_zero_run, axis1, arrvals) out out[zero_run 32] # 与容量矛盾当日峰值超过容量*max_ratio则丢弃 cap out[meter_id].map(capacity_map) over_cap out[feat_cols].max(axis1) cap * max_ratio out out[~over_cap] # 跳变点差分超过当日最大值的30%且一天超过3个 diff np.abs(np.diff(vals, axis1)) spike_ratio diff / (vals.max(axis1, keepdimsTrue) 1e-6) spike_cnt (spike_ratio 0.3).sum(axis1) out out[spike_cnt 3] return outmax_zero_run用最朴素的循环统计一行里最长连续零值长度比字符串分割写法直白维护成本低。容量判定时max_ratio1.5是常用起始值有冲击负荷的工业用户需要单独调高到2.0。跳变点阈值30%也不是固定值如果清洗完聚类结果中仍然出现大量单点尖峰簇说明阈值太宽可以收到20%重新跑一遍。3.2 特征列构造星期、节假日和峰谷标签不放进聚类矩阵聚类矩阵里只放功率曲线本身不放星期几、是否节假日这类外部标签。原因很直接这些变量是类别型直接进欧氏距离计算会稀释曲线形状的贡献而且一旦放进去工作日和休息日被人为拉开距离聚类结果就变成在按日历分类而不是按负荷形态分类。日期类型更适合在聚类完成之后做结果解释比如统计每个簇里工作日、周末、节假日的占比判断某个簇是“工作日双峰型”还是“周末单峰型”。日期特征在样本构造阶段仍然有用它作为辅助列挂在宽表旁边不参与标准化和距离计算。下面是我通常保留的辅助列辅助列含义聚类中用途weekday星期几(0-6)聚类后分析簇内构成is_holiday是否法定节假日识别节假日专属形态season季节判断是否需要分季节聚类季节对负荷形态影响很大落地时常有两种做法全量样本一起聚类再看簇与季节的对应关系或者先按季节分桶每个季节单独聚类。后者适合夏季空调负荷和冬季采暖负荷差异极大的地区但要注意春秋季样本量偏少单独聚类时k要相应调小。我一般优先做全量聚类再根据簇内季节构成决定要不要拆模型因为拆分后每个子模型都要重新选k运维成本直接翻倍。3.3 清洗效果检验三点确认再进聚类清洗逻辑写完一定要先做输出检验别直接丢给KMeans。我会检查三件事清洗前后样本量变化和淘汰率、清洗后96列各自的极值和均值、随机抽取20条曲线做可视化侧写。第一件事能发现清洗规则是不是过猛比如一个月数据被删掉40%大概率是阈值设错了。第二件事检查归一化指标如果某列全部是常数标准化之后会出现除零问题。第三件事靠肉眼扫一遍曲线确认双峰、单峰、平直、晚高峰等典型形态仍然存在于数据里。快速统计代码用来打印淘汰率和列极值范围。before len(wide) wide_clean clean_wide(wide, capacity_map).reset_index(dropTrue) drop_rate (before - len(wide_clean)) / before print(f清洗淘汰率: {drop_rate:.2%}) print(每列最大值范围: f{wide_clean[feat_cols].max().min():.1f} ~ f{wide_clean[feat_cols].max().max():.1f} kW)淘汰率超过20%时建议把清洗规则拆开单独统计逐条看是缺失、跳变还是容量矛盾引起的不要笼统地调阈值。因为清洗后还要保持与x_norm矩阵的行号对齐这里提前执行了reset_index(dropTrue)后面取典型曲线时就不会出现索引错位。可视化抽检我一般随机抽20条直接画折线图如果所有曲线都平滑到“完美”也要警惕真实数据总有噪声清洗过度会把高峰时段原本存在的波动抹掉最后聚类典型曲线会比实际偏平。这个阶段的目标是去掉坏数据不是让数据变好看。4. k值选择与聚类执行的三个必调参数4.1 肘部法则选kinertia拐点怎么读K-means必须提前给定簇数k这是整个流程里最需要业务判断的参数。最常用的辅助工具是肘部法则对一批候选k值分别跑聚类记录每个k对应的inertia再画折线图找拐点也就是inertia下降幅度骤减的位置。拐点意味着再增加簇数对整体紧凑度的改善已经有限。from sklearn.cluster import KMeans import matplotlib.pyplot as plt inertia_list [] k_candidates range(2, 11) for k in k_candidates: model KMeans(n_clustersk, n_init10, random_state42) model.fit(x_norm) inertia_list.append(model.inertia_) plt.plot(list(k_candidates), inertia_list, markero) plt.xlabel(k) plt.ylabel(inertia) plt.grid(True) plt.savefig(elbow.png, dpi120)注意inertia曲线在k增大时必然单调下降因为簇越多每个簇越紧凑。曲线从陡降转为平缓的位置对应的k就是参考值。实际读取拐点时有相当强的主观成分尤其当负荷形态接近时曲线可能没有明显拐点。所以我不只依赖这条曲线而是同时看轮廓系数再与业务方确认可解释的类别数上限。4.2 轮廓系数辅助判断业务口径最后拍板轮廓系数衡量一个样本与自己簇内样本的相似度、以及和最近其他簇样本相似度的差值范围在-1到1之间越接近1说明簇内紧致且簇间分离好。sklearn的silhouette_score可以直接对归一化矩阵计算。from sklearn.metrics import silhouette_score sil_list [] for k in k_candidates: model KMeans(n_clustersk, n_init10, random_state42) labels model.fit_predict(x_norm) sil_list.append(silhouette_score(x_norm, labels)) best_k_sil k_candidates[int(np.argmax(sil_list))] print(f轮廓系数最大值在 k{best_k_sil}, score{max(sil_list):.3f})轮廓系数最大的k与肘部法则读出的k经常不一致前者偏向于让所有簇都紧凑且可分后者偏向于数据自然层级。实际落地时业务口径往往比统计指标更关键如果业务方明确需要“工作日型、周末型、季节性尖峰型”三到四类典型曲线那么即使指标建议6个簇也应该在4的粒度先出一版结果看质心曲线是否还能被清晰解释。反过来如果指标给k3但某一簇内部曲线形态差异极大说明k取小了需要调大。4.3 三个必调参数n_init、max_iter、random_stateK-means对初始质心敏感同一份数据不同随机种子可能得到不同结果。n_init控制用多少个随机初始质心方案各自跑完整轮最终保留inertia最小的那套结果sklearn默认值是10大数据量时降到5可以提速max_iter是单次运行最大迭代轮数默认300在大多数负荷数据上够用如果出现未收敛警告提高到600random_state固定随机种子保证离线调参和线上预测时同一批数据能复现完全一致的聚类结果。典型曲线归档时如果同一个k前后跑出两套标签下游业务对不上信任问题处理起来比想象中麻烦得多。参数常见值调整方向n_init10结果波动大时提高到20max_iter300未收敛警告时提高到600random_state42固定一个正整数不建议改动还有个容易被忽略的参数是tol它控制两次迭代质心移动量的阈值默认1e-4。负荷数据归一化后数值范围较小质心移动量本来就小如果发现迭代次数异常少可以检查是不是tol设太大导致提前终止。质心保存在cluster_centers_里注意那是归一化空间的坐标必须配合前面保存的coef才能还原成kW量级。final_k best_k_sil if best_k_sil 3 else 4 model KMeans(n_clustersfinal_k, n_init20, max_iter600, random_state42, tol1e-6) labels model.fit_predict(x_norm) np.save(cluster_labels.npy, labels) dump(model, kmeans_model.joblib)fit_predict一步完成拟合和打标签。n_init20在数据量百米万行时也就增加十几秒运行时间但能让结果稳定很多我一般调参阶段用10出正式归档结果时提高到20。tol1e-6是为了让质心真正收敛消除归一化后数值范围过小带来的提前停止风险。5. 典型曲线提取后怎么用校验、归档与预测冷启动5.1 选离质心最近的真实样本而不是直接用质心不直接用质心当典型数据原因很实际质心是均值会被簇内离群点拉偏画出来过于平滑而真实样本对应具体的用户和日期可以回去做现场核验也可以在报告里写清楚“典型曲线取自某用户某日”可信度完全不同。提取逻辑很直接取每个簇中到质心欧氏距离最小的那条实际曲线。from scipy.spatial.distance import cdist # 确保wide_clean与x_norm行对齐 wide_clean wide_clean.reset_index(dropTrue) x_norm, _ build_feature_matrix(wide_clean, modeshape) for c in range(final_k): idx np.where(labels c)[0] dists cdist(x_norm[idx], model.cluster_centers_[c].reshape(1, -1)).squeeze() nearest wide_clean.iloc[idx[int(np.argmin(dists))]] print(nearest[meter_id], nearest[data_date], nearest[feat_cols].sum().round(1), kWh)这里在归一化空间选最近点避免大负荷用户因为数值大总是被选中打印出的电量是实际值可以直接归档。一个簇需要多条典型曲线时取距离排序后的前三条覆盖该簇内的主要波动形态。5.2 校验、落库与冷启动应用典型曲线落库前做两项校验每个簇样本量占比不低于5%否则把该簇合并到距离最近的簇簇内样本与质心的平均距离如果超过全量均值的两倍就要考虑k是否偏小。归档时用长表结构包含日期、计量点、cluster_id和96点功率不要用宽表宽表对查询和新维度扩展都不友好。新装用户冷启动预测是最常见的落地场景没有历史负荷时把现有的一两天数据归一化后计算到各质心距离归入最近簇再用该簇典型曲线乘以容量倍率生成初始预测曲线。归簇判定同样要在归一化空间完成线上预处理逻辑必须与离线聚类完全一致包括同一套归一化系数和特征列顺序。这条流程跑通后典型数据提取就不只是分析报告而是真正接入日常负荷预测数据流水线的一个环节。本文还有配套的精品资源点击获取