家庭用电预测——在 Keras 深度学习库中开发多变量时间序列预测的 LSTM 模型

[MoModel]

家庭用电预测

智能电表的兴起以及太阳能电池板等发电技术的广泛积累了大量的用电数据。该数据代表了功率的相关变量，这些变量又可用于建模甚至预测未来的电力消耗。像长期短期记忆网络（ LSTM ）这样的神经网络能够处理多个输入变量的问题。这在时间序列预测中具有很大的益处，而传统的线性方法难以适应多变量或多输入预测问题。

第一步 准备工作

在本教程中，您将了解如何在 Keras 深度学习库中开发多变量时间序列预测的 LSTM 模型。
完成本教程后，您将知道：

• 如何将原始数据集转换为可用于时间序列预测的类型
• 如何搭建解决多变量时间序列预测问题的 LSTM 模型
• 如何做出预测并将结果重新调整到原始单位

特征介绍
**

1. date: 日期格式为 dd/mm/yy
2. time: 时间格式为 hh:mm:ss
3. Global_active_power: 家庭消耗的总有功功率（千瓦），在交流电路中，电源在一个周期内发出瞬时功率的平均值(或负载电阻所消耗的功率)，称为"有功功率"
4. Global_reactive_power: 家庭消耗的总无功功率（千瓦），在具有电感或电容的电路中，在每半个周期内，把电源能量变成磁场(或电场)能量贮存起来，然后，再释放，又把贮存的磁场（或电场）能量再返回给电源，只是进行这种能量的交换，并没有真正消耗能量，我们把这个交换的功率值，称为" 无功功率"
5. voltage: 平均电压（伏特）
6. Global_intensity: 平均电流强度（安培)
7. sub_metering_1: 厨房的有功功率
8. sub_metering_2: 用于洗衣机等电器的有功功率
9. sub_metering_3: 空调热水器等电器的有功功率

这里我们使用的是时间序列预测模型，利用历史数据来预测之后的 Global_active_power。

第二步 导入数据并进行数据预处理

# 将 txt 文档读入并转换为 csv 文件格式
df = pd.read_csv(path, sep=';', 
                 parse_dates={'dt' : ['Date', 'Time']}, infer_datetime_format=True, 
                 low_memory=False, na_values=['nan','?'], index_col='dt')

上面代码执行了以下操作：
1 ）将 'Date' 和 'Time' 两列合并为 'dt'
2 ）将上面的数据转换为时间序列类型，将时间作为索引。

# 我们查看前 5 条数据
df.head()

我们可以看出 Global_active_power 大于 Global_reactive_power， voltage 基本稳定在 233 伏特。

# 了解数据的分布
df.describe()

我们可以通过上表了解数据的分布特征，比如均值和方差，还有最小值等等。

处理缺失值

在原始计量数据，特别是用户电量抽取过程中，发现存在缺失现象。若将这些值抛弃掉，会严重影响用电预测的结果。为了达到较好的建模效果，需要对缺失值进行处理。

# 找到所有有缺失值的列
total = df.isnull().sum().sort_values(ascending=False)
display(total)
# 用各列的均值填充缺失值
for j in range(0,7):        
        df.iloc[:,j] = df.iloc[:,j].fillna(df.iloc[:,j].mean())
# 查看是否还有缺失值
df.isnull().sum()

第三步 数据可视化

大部分真实的数据集都难以观察，因为它们有很多列变量，以及很多行数据。理解信息这方面大量都依赖于视觉。查看数据基本等价于了解数据。然而，基本上我们只能理解视觉上的二维或者三维数据，最好是二维。所以数据可视化能够帮助我们提高对数据的理解。

对数据集中呈现的结构和相关性进行观察，会让它们易于理解。一个准确的机器学习模型给出的预测，应当能够反映出数据集中所体现的结构和相关性。要明确一个模型给出的预测是否可信，对这些结构和相关性进行理解是首当其冲的。

了解数据分布

我们可以使用 resample 函数使特征按不同单位进行聚合。例如:使用参数 'H' 调用此函数使时间索引的数据按小时聚合。
下面我们对 Global_active_power 按天进行聚合，并比较它的总和和平均值。 可以看出，重采样数据集的平均值和总和具有相似的结构。

# 对 Global_active_power 数据按天进行聚合，并比较均值和总和
df.Global_active_power.resample('D').sum().plot(title='Global_active_power resampled over day for sum') 
plt.tight_layout()
plt.show()   

df.Global_active_power.resample('D').mean().plot(title='Global_active_power resampled over day for mean', color='red') 
plt.tight_layout()
plt.show()

# 对 'Global_active_power' 按季度进行聚合
df['Global_active_power'].resample('Q').mean().plot(kind='bar')
plt.xticks(rotation=60)
plt.ylabel('Global_active_power')
plt.title('Global_active_power per quarter (averaged over quarter)')
plt.show()

# 对'Voltage' 按月进行聚合
df['Voltage'].resample('M').mean().plot(kind='bar', color='red')
plt.xticks(rotation=60)
plt.ylabel('Voltage')
plt.title('Voltage per quarter (summed over quarter)')
plt.show()

我们可以看出每个月电压平均值变化幅度非常小，基本保持稳定。

# 下面我们比较对不同特征以天进行重采样的数值
cols = [0, 1, 2, 3, 5, 6]
i = 1
groups=cols
# 统计以天进行重采样的平均值
values = df.resample('D').mean().values
# 对每个 column 进行绘图
plt.figure(figsize=(15, 10))
for group in groups:
    # 对每个特征添加子图
	plt.subplot(len(cols), 1, i)
    # 进行绘图
	plt.plot(values[:, group])
    # 添加标题
	plt.title(df.columns[group], y=0.75, loc='right')
    # 更换子图位置
	i += 1
plt.show()

# 下面我们看看 ‘ Global_active_power ‘ 数值分布情况
sns.distplot(df['Global_active_power']);

可以看出家庭消耗的总有功功率主要集中在 0-2kw 范围内

特征相关性分析

# 查看 'Global_intensity' 和 'Global_active_power' 特征之间的关系
data_returns = df.pct_change()
sns.jointplot(x='Global_intensity', y='Global_active_power', data=data_returns)  
plt.show()

# 查看 'Voltage' 和 'Global_active_power' 之间的关系
sns.jointplot(x='Voltage', y='Global_active_power', data=data_returns)  
plt.show()

从上面的两个图中可以看出 'Global_intensity' 和 'Global_active_power' 是线性相关的。 但 'Voltage' 和 'Global_active_power' 的相关性较低， 这是机器学习所要观察的。

# 对各特征按月进行聚合
plt.title('resampled over month',size=15)
sns.heatmap(df.resample('M').mean().corr(),linewidths=0.1,vmax=1.0, 
            square=True, cmap=plt.cm.RdBu, linecolor='white', annot=True)

# 对各特征按年进行聚合
plt.title('resampled over year',size=15)
sns.heatmap(df.resample('A').mean().corr(),linewidths=0.1,vmax=1.0, 
            square=True, cmap=plt.cm.RdBu, linecolor='white', annot=True)

从上面可以看出，采用重采样技术可以改变特征之间的相关性， 这对于特征工程非常重要。

第四步 多变量 LSTM 预测模型

在可以使用深度学习之前，必须将时间序列预测问题调整为监督学习问题，形成输入和输出序列对，利用前一时间的 Global_active_power 和其他特征预测当前时间的 Global_active_power。

# 下面我们对 ‘ Global_active_power' 向前移动一个单位
df['Global_active_power'].resample('h').mean().shift(1)
# 下面我们对 ‘ Global_active_power' 向后移动一个单位
df['Global_active_power'].resample('h').mean().shift(-1)

因为这里我们预测 ‘ Global_active_power' 不仅用到过去时间的 ‘ Global_active_power' 还会用到其他的特征，比如：'Voltage ‘，这时候我们把此类问题叫做多变量时间序列预测。下面我们展示将时间序列问题转换为监督学习问题的过程。

def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
	n_vars = 1 if type(data) is list else data.shape[1]
	dff = pd.DataFrame(data)
	cols, names = list(), list()
	# 输入序列(t-n, ... t-1)
	for i in range(n_in, 0, -1):
		cols.append(dff.shift(i))
		names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]
	# 预测序列 (t, t+1, ... t+n)
	for i in range(0, n_out):
		cols.append(dff.shift(-i))
		if i == 0:
			names += [('var%d(t)' % (j+1)) for j in range(n_vars)]
		else:
			names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)]
	# 组合起来
	agg = pd.concat(cols, axis=1)
	agg.columns = names
	# 丢掉 NaN
	if dropnan:
		agg.dropna(inplace=True)
	return agg
 

我们将数据以小时进行聚合，这样也可以减少计算时间，可以快速获得测试模型的结果。 我们以小时进行聚合（原始数据以分钟为单位）。这将把数据大小从 2075259 减少到 34589，但依然保持数据的整体结构。

# 将数据按小时聚合
df_resample = df.resample('h').mean() 
df_resample.shape

对特征进行归一化

数据归一化处理是数据挖掘的一项基础工作。不同指标往往具有不同的量纲，数值间的差别可能很大，不进行处理可能会影响数据分析的结果。为了消除指标间的量纲和取值范围差异的影响，需要进行标准化处理，将数据按照比例进行缩放，使之落入特定的区域，以便于进行综合分析。

同样我们必须对用户电量的各指标进行处理，这里我们用到最小最大规范化。

# 把所有特征进行规范化
values = df_resample.values 
# 特征归一化
scaler = MinMaxScaler(feature_range=(0, 1))
scaled = scaler.fit_transform(values)
# 转化为监督问题

# 转化为监督问题
reframed = series_to_supervised(scaled, 1, 1)
# 删除不需要的特征
reframed.drop(reframed.columns[[8,9,10,11,12,13]], axis=1, inplace=True)
print(reframed.head())

把数据集拆分成训练集和测试集

这里，我们将前三年的数据作为训练集，后一年的样本作为测试集，并将数据改为 3 维格式。

# 对样本集拆分成训练集和测试集
values = reframed.values
n_train_time = 365*24*3
train = values[:n_train_time, :]
test = values[n_train_time:, :]
train_X, train_y = train[:, :-1], train[:, -1]
test_X, test_y = test[:, :-1], test[:, -1]
# 把数据转换为 3 维
train_X = train_X.reshape((train_X.shape[0], 1, train_X.shape[1]))
test_X = test_X.reshape((test_X.shape[0], 1, test_X.shape[1]))
print(train_X.shape, train_y.shape, test_X.shape, test_y.shape)

搭建 LSTM 模型

模型架构
1 ） LSTM 在第一个可见层中有 100 个神经元
2 ）丢弃 20 ％，防止过拟合
3 ）输出层中 1 个神经元，用于预测 Global_active_power
4 ）使用平均绝对误差（ MAE ）损失函数和随机梯度下降的 Adam 优化器
5 ）该模型 epoch 为 20，批次大小为 70 

# 搭建网络模型
model = Sequential()
model.add(LSTM(100, input_shape=(train_X.shape[1], train_X.shape[2])))
model.add(Dropout(0.2))
model.add(Dense(1))
model.compile(loss='mean_squared_error', optimizer='adam')
# 训练网络
history = model.fit(train_X, train_y, epochs=20, batch_size=70, validation_data=(test_X, test_y), verbose=2, shuffle=False)
# 统计 loss 值
plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='upper right')
plt.show()
# 做出预测
yhat = model.predict(test_X)
test_X = test_X.reshape((test_X.shape[0], 7))

我们可以看出模型的收敛速度很快。

第五步 模型评估

预测模型对训练集进行预测而得到的准确率并不能很好地反映预测模型对未来的性能，为了有效判断一
个预测模型的性能表现，需要一组没有参加预测模型建立的数据集，并在该模型上评价预测模型的准确率，这组独立的数据集叫做测试集。在测试集进行预测并且评估，我们怎样对模型进行性能衡量？

回归问题的评价指标:通常用相对／绝对误差，平均绝对误差，均方误差，均方根误差等指标来衡量,分类问题的评价指标：准确率，精确率，召回率，ROC 曲线，混淆矩阵。

我们将预测与测试数据集相结合，并反演缩放。
以预测值和实际值为原始尺度，我们可以计算模型的误差分数。 在这种情况下，我们计算出与变量本身相同的单位产生误差的均方根误差（ RMSE)。

#对预测值进行反演缩放
inv_yhat = np.concatenate((yhat, test_X[:, -6:]), axis=1)
inv_yhat = scaler.inverse_transform(inv_yhat)
inv_yhat = inv_yhat[:,0]
# 对真实值进行反演缩放
test_y = test_y.reshape((len(test_y), 1))
inv_y = np.concatenate((test_y, test_X[:, -6:]), axis=1)
inv_y = scaler.inverse_transform(inv_y)
inv_y = inv_y[:,0]
# 计算 RMSE
rmse = np.sqrt(mean_squared_error(inv_y, inv_yhat))
print('Test RMSE: %.3f' % rmse)

## 我们比较下 200 小时真实值和预测值 
aa=[x for x in range(200)]
plt.plot(aa, inv_y[:200], marker='.', label="actual")
plt.plot(aa, inv_yhat[:200], 'r', label="prediction")
plt.ylabel('Global_active_power', size=15)
plt.xlabel('Time step', size=15)
plt.legend(fontsize=15)
plt.show()

第六步 思考，如何进一步的改进模型？

能不能进一步的改进模型呢？下面提出了一些可能的改进模型的方案，有兴趣的话可以试一试哦。
1.在缺失值处理中利用其他的插值方法
2.使用复杂的模型
3.调整 epoch 和 batch_size

第七步 总结

今天我们一起制作了一个家庭用户用电预测的应用，大家可以在项目源码地址 fork 这个项目
http://www.momodel.cn:8899/explore/5cde0ed11afd94371e5697ff?type=app

我们首先对数据进行预处理，处理缺失值；接着进行数据可视化，了解数据的结构和相关性；然后搭建 LSTM 模型，其中最为关键的是将问题转化为监督学习问题；最后我们对模型进行评估，并提出了优化模型的建议。

使用我们的模型同样也可以预测温湿度和股价等等，只需要略加修改就行，来做出自己的应用吧。

---

参考资料：
https://wenku.baidu.com/view/3973baa6951ea76e58fafab069dc5022aaea46b9.html
https://blog.csdn.net/sinat_22510827/article/details/80996937
https://blog.csdn.net/weixin_40651515/article/details/83895707
https://www.jianshu.com/p/bebf8ca6a946

---

Mo（网址：momodel.cn）是一个支持 Python 的人工智能在线建模平台，能帮助你快速开发、训练并部署模型。

---

Mo 人工智能俱乐部 是由网站的研发与产品设计团队发起、致力于降低人工智能开发与使用门槛的俱乐部。团队具备大数据处理分析、可视化与数据建模经验，已承担多领域智能项目，具备从底层到前端的全线设计开发能力。主要研究方向为大数据管理分析与人工智能技术，并以此来促进数据驱动的科学研究。
目前俱乐部每周六在杭州举办以机器学习为主题的线下技术沙龙活动，不定期进行论文分享与学术交流。希望能汇聚来自各行各业对人工智能感兴趣的朋友，不断交流共同成长，推动人工智能民主化、应用普及化。