非侵入式负荷检测与分解：电力数据挖掘新视角

电力数据挖掘

概述
案例背景
分析目标
分析过程
数据准备
- 数据探索
- 缺失值处理
属性构造
- 设备数据
- 周波数据
- 模型训练
性能度量
推荐阅读

在这里插入图片描述

主页传送门：📀 传送

概述

摘要：本案例将根据已收集到的电力数据，深度挖掘各电力设备的电流、电压和功率等情况，分析各电力设备的实际用电量，进而为电力公司制定电能能源策略提供一定的参考依据。更多详细内容请参考《Python数据挖掘：入门进阶与实用案例分析》一书。

案例背景

为了更好地监测用电设备的能耗情况，电力分项计量技术随之诞生。电力分项计量对于电力公司准确预测电力负荷、科学制定电网调度方案、提高电力系统稳定性和可靠性有着重要意义。对用户而言，电力分项计量可以帮助用户了解用电设备的使用情况，提高用户的节能意识，促进科学合理用电。

在这里插入图片描述

分析目标

本案例根据非侵入式负荷检测与分解的电力数据挖掘的背景和业务需求，需要实现的目标如下:

分析每个用电设备的运行属性。
构建设备判别属性库。
利用K最近邻模型，实现从整条线路中“分解”出每个用电设备的独立用电数据。

分析过程

在这里插入图片描述

数据准备

数据探索

在本案例的电力数据挖掘分析中，不会涉及操作记录数据。因此，此处主要获取设备数据、周波数据和谐波数据。在获取数据后，由于数据表较多，每个表的属性也较多，所以需要对数据进行数据探索分析。在数据探索过程中主要根据原始数据特点，对每个设备的不同属性对应的数据进行可视化，得到的部分结果如图1~图3所示。

无功功率和总无功功率：

在这里插入图片描述

电流轨迹：
在这里插入图片描述

电压轨迹：
在这里插入图片描述

可视化数据代码示例：

import pandas as pdimport matplotlib.pyplot as pltimport osfilename = os.listdir('../data/附件1')  # 得到文件夹下的所有文件名称n_filename = len(filename)  # 给各设备的数据添加操作信息，画出各属性轨迹图并保存def fun(a):save_name = ['YD1', 'YD10', 'YD11', 'YD2', 'YD3', 'YD4','YD5', 'YD6', 'YD7', 'YD8', 'YD9']plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号for i in range(a):Sb = pd.read_excel('../data/附件1/' + filename[i], '设备数据', index_col = None)Xb = pd.read_excel('../data/附件1/' + filename[i], '谐波数据', index_col = None)Zb = pd.read_excel('../data/附件1/' + filename[i], '周波数据', index_col = None)# 电流轨迹图plt.plot(Sb['IC'])plt.title(save_name[i] + '-IC')plt.ylabel('电流（0.001A）')plt.show()# 电压轨迹图lt.plot(Sb['UC'])plt.title(save_name[i] + '-UC')plt.ylabel('电压（0.1V）')plt.show()# 有功功率和总有功功率plt.plot(Sb[['PC', 'P']])plt.title(save_name[i] + '-P')plt.ylabel('有功功率（0.0001kW）')plt.show()# 无功功率和总无功功率plt.plot(Sb[['QC', 'Q']])plt.title(save_name[i] + '-Q')plt.ylabel('无功功率（0.0001kVar）')plt.show()# 功率因数和总功率因数plt.plot(Sb[['PFC', 'PF']])plt.title(save_name[i] + '-PF')plt.ylabel('功率因数（%）')plt.show()# 谐波电压plt.plot(Xb.loc[:, 'UC02':].T)plt.title(save_name[i] + '-谐波电压')plt.show()# 周波数据plt.plot(Zb.loc[:, 'IC001':].T)plt.title(save_name[i] + '-周波数据')plt.show()fun(n_filename)

缺失值处理

通过数据探索，发现数据中部分“time”属性存在缺失值，需要对这部分缺失值进行处理。由于每份数据中“time”属性的缺失时间段长不同，所以需要进行不同的处理。对于每个设备数据中具有较大缺失时间段的数据进行删除处理，对于具有较小缺失时间段的数据使用前一个值进行插补。

在进行缺失值处理之前，需要将训练数据中所有设备数据中的设备数据表、周波数据表、谐波数据表和操作记录表，以及测试数据中所有设备数据中的设备数据表、周波数据表和谐波数据表都提取出来，作为独立的数据文件，生成的部分文件如图4所示。

在这里插入图片描述

提取数据文件代码示例：

# 将xlsx文件转化为CSV文件import globimport pandas as pdimport mathdef file_transform(xls):print('共发现%s个xlsx文件' % len(glob.glob(xls)))print('正在处理............')for file in glob.glob(xls):  # 循环读取同文件夹下的xlsx文件combine1 = pd.read_excel(file, index_col=0, sheet_name=None)for key in combine1:combine1[key].to_csv('../tmp/' + file[8: -5] + key + '.csv', encoding='utf-8')print('处理完成')xls_list = ['../data/附件1/*.xlsx', '../data/附件2/*.xlsx']file_transform(xls_list[0])  # 处理训练数据file_transform(xls_list[1])  # 处理测试数据

提取数据文件完成后，对提取的数据文件进行缺失值处理，处理后生成的部分文件如图5所示。

在这里插入图片描述
缺失值处理代码示例：

# 对每个数据文件中较大缺失时间点数据进行删除处理，较小缺失时间点数据进行前值替补def missing_data(evi):print('共发现%s个CSV文件' % len(glob.glob(evi)))for j in glob.glob(evi):fr = pd.read_csv(j, header=0, encoding='gbk')fr['time'] = pd.to_datetime(fr['time'])helper = pd.DataFrame({'time': pd.date_range(fr['time'].min(), fr['time'].max(), freq='S')})fr = pd.merge(fr, helper, on='time', how='outer').sort_values('time')fr = fr.reset_index(drop=True)frame = pd.DataFrame()for g in range(0, len(list(fr['time'])) - 1):if math.isnan(fr.iloc[:, 1][g + 1]) and math.isnan(fr.iloc[:, 1][g]):continueelse:scop = pd.Series(fr.loc[g])frame = pd.concat([frame, scop], axis=1)frame = pd.DataFrame(frame.values.T, index=frame.columns, columns=frame.index)frames = frame.fillna(method='ffill')frames.to_csv(j[:-4] + '1.csv', index=False, encoding='utf-8')print('处理完成')evi_list = ['../tmp/附件1/*数据.csv', '../tmp/附件2/*数据.csv']missing_data(evi_list[0])  # 处理训练数据missing_data(evi_list[1])  # 处理测试数据

属性构造

虽然在数据准备过程中对属性进行了初步处理，但是引入的属性太多，而且这些属性之间存在重复的信息。为了保留重要的属性，建立精确、简单的模型，需要对原始属性进一步筛选与构造。

设备数据

在数据探索过程中发现，不同设备的无功功率、总无功功率、有功功率、总有功功率、功率因数和总功率因数差别很大，具有较高的区分度，故本案例选择无功功率、总无功功率、有功功率、总有功功率、功率因数和总功率因数作为设备数据的属性构建判别属性库。

处理好缺失值后，每个设备的数据都由一张表变为了多张表，所以需要将相同类型的数据表合并到一张表中，如将所有设备的设备数据表合并到一张表当中。同时，因为缺失值处理的其中一种方式是使用前一个值进行插补，所以产生了相同的记录，需要对重复出现的记录进行处理，处理后生成的数据表如表1所示。
在这里插入图片描述
合并且去重设备数据代码示例：

import globimport pandas as pdimport os# 合并11个设备数据及处理合并中重复的数据def combined_equipment(csv_name):# 合并print('共发现%s个CSV文件' % len(glob.glob(csv_name)))print('正在处理............')for i in glob.glob(csv_name):  # 循环读取同文件夹下的CSV文件fr = open(i, 'rb').read()file_path = os.path.split(i)with open(file_path[0] + '/device_combine.csv', 'ab') as f:f.write(fr)print('合并完毕！')# 去重df = pd.read_csv(file_path[0] + '/device_combine.csv', header=None, encoding='utf-8')datalist = df.drop_duplicates()datalist.to_csv(file_path[0] + '/device_combine.csv', index=False, header=0)print('去重完成')csv_list = ['../tmp/附件1/*设备数据1.csv', '../tmp/附件2/*设备数据1.csv']combined_equipment(csv_list[0])  # 处理训练数据combined_equipment(csv_list[1])  # 处理测试数据

周波数据

在数据探索过程中发现，周波数据中的电流随着时间的变化有较大的起伏，不同设备的周波数据中的电流绘制出来的折线图的起伏不尽相同，具有明显的差异，故本案例选择波峰和波谷作为周波数据的属性构建判别属性库。

由于原始的周波数据中并未存在电流的波峰和波谷两个属性，所以需要进行属性构建，构建生成的数据表如表2所示。

在这里插入图片描述
构建周波数据中的属性代码示例：

# 求取周波数据中电流的波峰和波谷作为属性参数import globimport pandas as pdfrom sklearn.cluster import KMeansimport osdef cycle(cycle_file):for file in glob.glob(cycle_file):cycle_YD = pd.read_csv(file, header=0, encoding='utf-8')cycle_YD1 = cycle_YD.iloc[:, 0:128]models = []for types in range(0, len(cycle_YD1)):model = KMeans(n_clusters=2, random_state=10)model.fit(pd.DataFrame(cycle_YD1.iloc[types, 1:]))  # 除时间以外的所有列models.append(model)# 相同状态间平稳求均值mean = pd.DataFrame()for model in models:r = pd.DataFrame(model.cluster_centers_, )  # 找出聚类中心r = r.sort_values(axis=0, ascending=True, by=[0])mean = pd.concat([mean, r.reset_index(drop=True)], axis=1)mean = pd.DataFrame(mean.values.T, index=mean.columns, columns=mean.index)mean.columns = ['波谷', '波峰']mean.index = list(cycle_YD['time'])mean.to_csv(file[:-9] + '波谷波峰.csv', index=False, encoding='gbk ')cycle_file = ['../tmp/附件1/*周波数据1.csv', '../tmp/附件2/*周波数据1.csv']cycle(cycle_file[0])  # 处理训练数据cycle(cycle_file[1])  # 处理测试数据# 合并周波的波峰波谷文件def merge_cycle(cycles_file):means = pd.DataFrame()for files in glob.glob(cycles_file):mean0 = pd.read_csv(files, header=0, encoding='gbk')means = pd.concat([means, mean0])file_path = os.path.split(glob.glob(cycles_file)[0])means.to_csv(file_path[0] + '/zuhe.csv', index=False, encoding='gbk')print('合并完成')cycles_file = ['../tmp/附件1/*波谷波峰.csv', '../tmp/附件2/*波谷波峰.csv']merge_cycle(cycles_file[0])  # 训练数据merge_cycle(cycles_file[1])  # 测试数据

模型训练

在判别设备种类时，选择K最近邻模型进行判别，利用属性构建而成的属性库训练模型，然后利用训练好的模型对设备1和设备2进行判别。构建判别模型并对设备种类进行判别，如代码清单6所示。

建立判别模型并对设备种类进行判别代码示例：

import globimport pandas as pdfrom sklearn import neighborsimport pickleimport os# 模型训练def model(test_files, test_devices):# 训练集zuhe = pd.read_csv('../tmp/附件1/zuhe.csv', header=0, encoding='gbk')device_combine = pd.read_csv('../tmp/附件1/device_combine.csv', header=0, encoding='gbk')train = pd.concat([zuhe, device_combine], axis=1)train.index = train['time'].tolist()  # 把“time”列设为索引train = train.drop(['PC', 'QC', 'PFC', 'time'], axis=1)train.to_csv('../tmp/' + 'train.csv', index=False, encoding='gbk')# 测试集for test_file, test_device in zip(test_files, test_devices):test_bofeng = pd.read_csv(test_file, header=0, encoding='gbk')test_devi = pd.read_csv(test_device, header=0, encoding='gbk')test = pd.concat([test_bofeng, test_devi], axis=1)test.index = test['time'].tolist()  # 把“time”列设为索引test = test.drop(['PC', 'QC', 'PFC', 'time'], axis=1)# K最近邻clf = neighbors.KNeighborsClassifier(n_neighbors=6, algorithm='auto')clf.fit(train.drop(['label'], axis=1), train['label'])predicted = clf.predict(test.drop(['label'], axis=1))predicted = pd.DataFrame(predicted)file_path = os.path.split(test_file)[1]test.to_csv('../tmp/' + file_path[:3] + 'test.csv', encoding='gbk')predicted.to_csv('../tmp/' + file_path[:3] + 'predicted.csv', index=False, encoding='gbk')with open('../tmp/' + file_path[:3] + 'model.pkl', 'ab') as pickle_file:pickle.dump(clf, pickle_file)print(clf)model(glob.glob('../tmp/附件2/*波谷波峰.csv'),glob.glob('../tmp/附件2/*设备数据1.csv'))

性能度量

根据代码清单6的设备判别结果，对模型进行模型评估，得到的结果如下
在这里插入图片描述
混淆矩阵如图7所示：

ROC曲线如图8所示：

模型评估代码示例：

import globimport pandas as pdimport matplotlib.pyplot as pltimport seaborn as snsfrom sklearn import metricsfrom sklearn.preprocessing import label_binarizeimport osimport pickle# 模型评估def model_evaluation(model_file, test_csv, predicted_csv):for clf, test, predicted in zip(model_file, test_csv, predicted_csv):with open(clf, 'rb') as pickle_file:clf = pickle.load(pickle_file)test = pd.read_csv(test, header=0, encoding='gbk')predicted = pd.read_csv(predicted, header=0, encoding='gbk')test.columns = ['time', '波谷', '波峰', 'IC', 'UC', 'P', 'Q', 'PF', 'label']print('模型分类准确度：', clf.score(test.drop(['label', 'time'], axis=1), test['label']))print('模型评估报告：\n', metrics.classification_report(test['label'], predicted))confusion_matrix0 = metrics.confusion_matrix(test['label'], predicted)confusion_matrix = pd.DataFrame(confusion_matrix0)class_names = list(set(test['label']))tick_marks = range(len(class_names))sns.heatmap(confusion_matrix, annot=True, cmap='YlGnBu', fmt='g')plt.xticks(tick_marks, class_names)plt.yticks(tick_marks, class_names)plt.tight_layout()plt.title('混淆矩阵')plt.ylabel('真实标签')plt.xlabel('预测标签')plt.show()y_binarize = label_binarize(test['label'], classes=class_names)predicted = label_binarize(predicted, classes=class_names)fpr, tpr, thresholds = metrics.roc_curve(y_binarize.ravel(), predicted.ravel())auc = metrics.auc(fpr, tpr)print('计算auc：', auc)  # 绘图plt.figure(figsize=(8, 4))lw = 2plt.plot(fpr, tpr, label='area = %0.2f' % auc)plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')plt.fill_between(fpr, tpr, alpha=0.2, color='b')plt.xlim([0.0, 1.0])plt.ylim([0.0, 1.05])plt.xlabel('1-特异性')plt.ylabel('灵敏度')plt.title('ROC曲线')plt.legend(loc='lower right')plt.show()model_evaluation(glob.glob('../tmp/*model.pkl'),glob.glob('../tmp/*test.csv'),glob.glob('../tmp/*predicted.csv'))

根据分析目标，需要计算实时用电量。实时用电量计算的是瞬时的用电器的电流、电压和时间的乘积，公式如下。
在这里插入图片描述
其中，为实时用电量，单位是0.001kWh。为功率，单位为W。

实时用电量计算，得到的实时用电量如表3所示。

在这里插入图片描述
计算实时用电量

# 计算实时用电量并输出状态表def cw(test_csv, predicted_csv, test_devices):for test, predicted, test_device in zip(test_csv, predicted_csv, test_devices):# 划分预测出的时刻表test = pd.read_csv(test, header=0, encoding='gbk')test.columns = ['time', '波谷', '波峰', 'IC', 'UC', 'P', 'Q', 'PF', 'label']test['time'] = pd.to_datetime(test['time'])test.index = test['time']predicteds = pd.read_csv(predicted, header=0, encoding='gbk')predicteds.columns = ['label']indexes = []class_names = list(set(test['label']))for j in class_names:index = list(predicteds.index[predicteds['label'] == j])indexes.append(index)# 取出首位序号及时间点from itertools import groupby  # 连续数字dif_indexs = []time_indexes = []info_lists = pd.DataFrame()for y, z in zip(indexes, class_names):dif_index = []fun = lambda x: x[1] - x[0]for k, g in groupby(enumerate(y), fun):dif_list = [j for i, j in g]  # 连续数字的列表if len(dif_list) > 1:scop = min(dif_list)  # 选取连续数字范围中的第一个else:scop = dif_list[0   ]dif_index.append(scop)time_index = list(test.iloc[dif_index, :].index)time_indexes.append(time_index)info_list = pd.DataFrame({'时间': time_index, 'model_设备状态': [z] * len(time_index)})dif_indexs.append(dif_index)info_lists = pd.concat([info_lists, info_list])# 计算实时用电量并保存状态表test_devi = pd.read_csv(test_device, header=0, encoding='gbk')test_devi['time'] = pd.to_datetime(test_devi['time'])test_devi['实时用电量'] = test_devi['P'] * 100 / 3600info_lists = info_lists.merge(test_devi[['time', '实时用电量']],how='inner', left_on='时间', right_on='time')info_lists = info_lists.sort_values(by=['时间'], ascending=True)info_lists = info_lists.drop(['time'], axis=1)file_path = os.path.split(test_device)[1]info_lists.to_csv('../tmp/' + file_path[:3] + '状态表.csv', index=False, encoding='gbk')print(info_lists)cw(glob.glob('../tmp/*test.csv'),glob.glob('../tmp/*predicted.csv'),glob.glob('../tmp/附件2/*设备数据1.csv'))