前言

与其他机器学习算法不同，卷积神经网络能够从序列数据中自动学习特征，支持多变量数据，并可直接输出用于多步预测的向量。一维CNN已被证明可以很好地执行，甚至在具有挑战性的序列预测问题上也能达到最新的结果。

计划用两篇文章介绍如何开发 1D CNN 进行多步时间序列预测。主要内容如下：

如何为单变量数据开发多步时间序列预测的CNN模型；
如何为多变量数据开发多通道多步时间序列预测的CNN模型；
如何为多变量数据开发多头多步时间序列预测的CNN模型。

本文介绍如何为单变量数据开发多步时间序列预测的CNN模型。

代码环境：

python 3.7.6
tensorflow 2.1.0
keras 2.3.1

代码在 jupyter notebook编写。完整代码部分对比较难理解的地方做了注释，应该很容易理解。

适用于多时间步预测的CNN模型

数据处理和模型评估部分之前的文章已经介绍过了，本文不再赘述。之前的文章：传送门

卷积神经网络（CNN）可用于多步时间序列预测。CNN也可用于递归或直接预测策略，在该策略中，模型进行一步预测，并为每个要预测的时间步开发一个模型。或者，CNN可以用来预测整个输出序列，作为整个矢量的一步预测。这是前馈神经网络的普遍优势。使用CNN的一个重要好处是：支持多个1D输入以进行预测。如果多步输出序列是一个以上输入序列的函数，则此功能很有用。这可以使用两种不同的模型配置来实现。

多个输入通道（Multiple Input Channels）：每个输入序列作为一个单独的通道读取，可类比图像的不同颜色通道（如RGB）。
多个输入头（Multiple Input Heads）：每个输入序列由不同的CNN子模型读取，并且模型的内部表征在解释之前先进行组合再进行预测。按个人理解，命名为子模型异构。

本文介绍三种业务场景下的CNN模型：

单变量输入数据的多步时间序列预测CNN模型；
多变量多输入通道的多步时间序列预测CNN模型；
多变量子模型异构的多步时间序列预测CNN模型；

这些模型会在家庭用电量数据集上进行演示。通过之前的文章可知，如果一个模型比一个总RMSE约为465千瓦的七天预测朴素模型的RMSE小，就可以认为该模型是可用的。本文示例并没有调整超参数以获得最佳性能，所选择的结构和超参数都经过了小的尝试和误差。考虑到模型的随机性，最好对模型进行多次评估，并输出测试数据集的平均性能。

1 单变量多步预测 CNN 模型

1.1 业务需求

开发卷积神经网络进行多步时间序列预测，该预测使用上篇文章中转换好的日功耗单变量序列。业务需求：给定之前几天的总日耗电量，预测下一个标准周（周天开始，周六结束）的日耗电量。

1.2 1D CNN 模型

一维CNN模型要求输入数据的shape为：[样本，时间步长，特征]（[samples, timesteps, features]）。一个样本（sample）包含一周7天的日总有功功率，即滑动窗口的宽度为7；特征只有1个（原数据集有七个特征，不包含日期和时间），即7天的日总有功功率序列。训练数据集有159周的数据（将原数据集的每分钟采样一次的数据，重采样成每天采样数据，然后把前三年的作为训练集，最后一年的数据作为测试集，截取的都是完整周，丢弃不完整周的日数据。相关细节👓请看这里），因此训练集的shape为：[159,7,1]。

将原数据集（每分钟采样一次）重采样成每天的数据，新数据的shape为：

(1442, 8)

这种格式的数据将使用之前的标准周来预测下一个标准周。问题是对于神经网络来说，159个训练样本是远远不够的。创建更多训练数据的一种方法是在训练期间做如下更改：根据之前的7天预测未来7天，而不考虑是否为标准周，可以理解滑动窗口的步长由7变为1，相当于增加了样本数（samples），特征数（features）不变，即行数增加，列数不变。而不是，相当于人为地在原数据集每两个样本之间又增加了6个新样本。再说详细一点就是，现在标准周预测数据相当于只有159行，2列，第一行是前一周的日总有功功率及期望预测输出（下一周的数据），第二行是之后一周的日总有功功率及期望输出预测（下下周的数据）…以此类推，最后行数（样本数变为159×7行）。原来的训练样本是使用当前一个标准周（周天开始，周六结束）的数据预测下一个标准周；使用下一个标准周预测下下个标准周…以此类推。而现在是放开标准周的限制来预测，比如从周一开始，周天结束（预测的总天数不变，同样是七天）；周二开始，周一结束；周三开始，周二结束…以此类推。细说就是，由通过之前一周中，从周一开始到周天结束的非标准周的预测结果，相当于添加一行（一个样本）；接着添加从周二开始到下周周一结束的非标准周预测结果，又相当于添加了一行…以此类推。测试集不做改变，仍然使用标准周进行预测。经过以上分析，训练数据的shape由 [159,7,1] 变为 [159*7,7,1]。其实最关键的还是段首说的那句话：滑动窗口的滑动步长由7变为1，相当于增加了样本数，特征数不变。画个图说明一下数据集扩充的思路：
在这里插入图片描述
代码实现：

def sliding_window(train, sw_width=7, in_start=0):'''该函数实现窗口宽度为7、滑动步长为1的滑动窗口截取序列数据'''data = train.reshape((train.shape[0] * train.shape[1], train.shape[2])) # 将以周为单位的样本展平为以天为单位的序列X, y = [], []for _ in range(len(data)):in_end = in_start + sw_widthout_end = in_end + sw_width# 保证截取样本完整，最大元素索引不超过原序列索引，则截取数据；否则丢弃该样本if out_end < len(data):# 训练数据以滑动步长1截取train_seq = data[in_start:in_end, 0]train_seq = train_seq.reshape((len(train_seq), 1))X.append(trian_seq)y.append(data[in_end:out_end], 0)in_start += 1return np.array(X), np.array(y)

在训练集上运行该函数，将159个样本转换为1099个样本。转换后，数据集的形状（shape）为 X=[1099, 7, 1] 和 y=[1099, 7]。接下来，可以在训练数据上定义和拟合CNN模型。这个多步时间序列预测问题是一个自回归问题（只有一个特征，日总有功功率）。

较少的特征和数据意味着只需要一个简单的模型即能满足业务需求。定义参数配置为 filters=16 和 kernel_size=3 的卷积层。这意味着7天的输入序列将通过卷积操作一次读取3个时间步长的值，该操作将执行16次。池化层把这些特征映射（feature maps）缩小1/4，之后将内部表示压缩为一个长向量。然后，全连接的层对其进行解释，最后在输出层输出序列预测结果。

使用均方误差损失函数（mse），因为它与我们选择的RMSE误差度量很好地匹配。使用随机梯度下降算法Adam，epochs=20，batch_size=4。该算法的小批量和随机性意味着同一个模型在每次训练时将学习输入到输出的略有不同的映射。这意味着评估模型时，结果可能会有所不同。可以尝试多次运行模型并计算模型性能的平均值。

1.3 完整代码

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt# 设置中文显示
plt.rcParams['font.sans-serif'] = ['Microsoft JhengHei']
plt.rcParams['axes.unicode_minus'] = Falseimport math
import sklearn.metrics as skm
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, Flatten
from tensorflow.keras.layers import Conv1D,MaxPooling1Ddef split_dataset(data):'''该函数实现以周为单位切分训练数据和测试数据'''# data为按天的耗电量统计数据，shape为(1442, 8)# 测试集取最后一年的46周（322天）数据，剩下的159周（1113天）数据为训练集，以下的切片实现此功能。train, test = data[1:-328], data[-328:-6]train = np.array(np.split(train, len(train)/7)) # 将数据划分为按周为单位的数据test = np.array(np.split(test, len(test)/7))return train, testdef evaluate_forecasts(actual, predicted):'''该函数实现根据预期值评估一个或多个周预测损失思路：统计所有单日预测的 RMSE'''scores = list()for i in range(actual.shape[1]):mse = skm.mean_squared_error(actual[:, i], predicted[:, i])rmse = math.sqrt(mse)scores.append(rmse)s = 0 # 计算总的 RMSEfor row in range(actual.shape[0]):for col in range(actual.shape[1]):s += (actual[row, col] - predicted[row, col]) ** 2score = math.sqrt(s / (actual.shape[0] * actual.shape[1]))print('actual.shape[0]:{}, actual.shape[1]:{}'.format(actual.shape[0], actual.shape[1]))return score, scoresdef summarize_scores(name, score, scores):s_scores = ', '.join(['%.1f' % s for s in scores])print('%s: [%.3f] %s\n' % (name, score, s_scores))def sliding_window(train, sw_width=7, n_out=7, in_start=0):'''该函数实现窗口宽度为7、滑动步长为1的滑动窗口截取序列数据'''data = train.reshape((train.shape[0] * train.shape[1], train.shape[2])) # 将以周为单位的样本展平为以天为单位的序列X, y = [], []for _ in range(len(data)):in_end = in_start + sw_widthout_end = in_end + n_out# 保证截取样本完整，最大元素索引不超过原序列索引，则截取数据；否则丢弃该样本if out_end < len(data):# 训练数据以滑动步长1截取train_seq = data[in_start:in_end, 0]train_seq = train_seq.reshape((len(train_seq), 1))X.append(train_seq)y.append(data[in_end:out_end, 0])in_start += 1return np.array(X), np.array(y)def cnn_model(train, sw_width, in_start=0, verbose_set=0, epochs_num=20, batch_size_set=4):'''该函数定义 1D CNN 模型'''train_x, train_y = sliding_window(train, sw_width, in_start=0)n_timesteps, n_features, n_outputs = train_x.shape[1], train_x.shape[2], train_y.shape[1]model = Sequential()model.add(Conv1D(filters=16, kernel_size=3, activation='relu', input_shape=(n_timesteps, n_features)))model.add(MaxPooling1D(pool_size=2))model.add(Flatten())model.add(Dense(10, activation='relu'))model.add(Dense(units=n_outputs))model.compile(loss='mse', optimizer='adam', metrics=['accuracy'])print(model.summary())model.fit(train_x, train_y,epochs=epochs_num, batch_size=batch_size_set, verbose=verbose_set)return modeldef forecast(model, pred_seq, sw_width):'''该函数实现对输入数据的预测'''data = np.array(pred_seq)data = data.reshape((data.shape[0]*data.shape[1], data.shape[2]))input_x = data[-sw_width:, 0] # 获取输入数据的最后一周的数据input_x = input_x.reshape((1, len(input_x), 1)) # 重塑形状[1, sw_width, 1]yhat = model.predict(input_x, verbose=0) # 预测下周数据yhat = yhat[0] # 获取预测向量return yhatdef evaluate_model(model, train, test, sd_width):'''该函数实现模型评估'''history_fore = [x for x in train]predictions = list() # 用于保存每周的前向验证结果；for i in range(len(test)):yhat_sequence = forecast(model, history_fore, sd_width) # 预测下周的数据predictions.append(yhat_sequence) # 保存预测结果history_fore.append(test[i, :]) # 得到真实的观察结果并添加到历史中以预测下周predictions = np.array(predictions) # 评估一周中每天的预测结果score, scores = evaluate_forecasts(test[:, :, 0], predictions)return score, scoresdef model_plot(score, scores, days, name):'''该函数实现绘制RMSE曲线图'''plt.figure(figsize=(8,6), dpi=150)plt.plot(days, scores, marker='o', label=name)plt.grid(linestyle='--', alpha=0.5)plt.ylabel(r'$RMSE$', size=15)plt.title('CNN 模型预测结果',  size=18)plt.legend()plt.show()def main_run(dataset, sw_width, days, name, in_start, verbose, epochs, batch_size):'''主函数：数据处理、模型训练流程'''# 划分训练集和测试集train, test = split_dataset(dataset.values)# 训练模型model = cnn_model(train, sw_width, in_start, verbose_set=0, epochs_num=20, batch_size_set=4)# 计算RMSEscore, scores = evaluate_model(model, train, test, sw_width)# 打印分数summarize_scores(name, score, scores)# 绘图model_plot(score, scores, days, name)if __name__ == '__main__':dataset = pd.read_csv('household_power_consumption_days.csv', header=0, infer_datetime_format=True, engine='c',parse_dates=['datetime'], index_col=['datetime'])days = ['sun', 'mon', 'tue', 'wed', 'thr', 'fri', 'sat']name = 'cnn'sliding_window_width=7input_sequence_start=0epochs_num=20batch_size_set=4verbose_set=0main_run(dataset, sliding_window_width, days, name, input_sequence_start,verbose_set, epochs_num, batch_size_set)

输出：

Model: "sequential_1"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
conv1d_1 (Conv1D)            (None, 5, 16)             64        
_________________________________________________________________
max_pooling1d_1 (MaxPooling1 (None, 2, 16)             0         
_________________________________________________________________
flatten_1 (Flatten)          (None, 32)                0         
_________________________________________________________________
dense_2 (Dense)              (None, 10)                330       
_________________________________________________________________
dense_3 (Dense)              (None, 7)                 77        
=================================================================
Total params: 471
Trainable params: 471
Non-trainable params: 0
_________________________________________________________________
None
actual.shape[0]:46, actual.shape[1]:7
cnn: [400.800] 429.3, 391.4, 346.6, 395.3, 404.1, 320.0, 494.9

评估模型，打印一周预测的总体RMSE，以及一周中每天的每天RMSE。我们可以看到，在这种情况下，该模型比朴素模型性能更好，实现了400千瓦左右的整体RMSE，低于朴素模型实现的465千瓦。
在这里插入图片描述

每日RMSE图。该图显示，周二和周五的总有功功率比其他时间更容易预测，按标准周预测时，周六是最难预测的一天。

可以把滑动窗口宽度设置为14 重新实验：

sliding_window_width = 14
main_run(dataset, sliding_window_width, days, name, input_sequence_start,verbose_set, epochs_num, batch_size_set)

输出：

Model: "sequential_2"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
conv1d_2 (Conv1D)            (None, 12, 16)            64        
_________________________________________________________________
max_pooling1d_2 (MaxPooling1 (None, 6, 16)             0         
_________________________________________________________________
flatten_2 (Flatten)          (None, 96)                0         
_________________________________________________________________
dense_4 (Dense)              (None, 10)                970       
_________________________________________________________________
dense_5 (Dense)              (None, 7)                 77        
=================================================================
Total params: 1,111
Trainable params: 1,111
Non-trainable params: 0
_________________________________________________________________
None
actual.shape[0]:46, actual.shape[1]:7
cnn: [383.272] 390.3, 392.4, 348.3, 386.5, 374.8, 311.2, 462.5