时间序列预测 — CNN-LSTM-Attention实现多变量负荷预测(Tensorflow)：多变量滚动

???专栏链接：https://blog.csdn.net/qq_41921826/category_12495091.html

专栏内容

??所有文章提供源代码、数据集、效果可视化

??文章多次上领域内容榜、每日必看榜单、全站综合热榜

?

?

?

?

?

?

?

时间序列预测存在的问题

??现有的大量方法没有真正的预测未来值，只是用历史数据做验证

??利用时间序列分解算法存在信息泄露的问题：有人用emd+lstm对时间序列进行预测，是否存在原理上的问题? - 知乎

---

目录

1?数据处理

1.1 导入库文件

1.2 导入数据集

?1.3 缺失值分析

?2 构造训练数据

3 CNN_LSTM-Attention模型训练

3.1 CNN_LSTM模型?

3.2 搭建Attention模型

3.3 搭建CNN_LSTM-Attention模型

4 CNN_LSTM-Attention模型预测

4.1 分量预测

4.2 可视化

---

1?数据处理

1.1 导入库文件

import scipy
import pandas as pd
import numpy as np
import math
import datetime
from matplotlib import pyplot as plt

# 导入深度学习框架tensorflow
import tensorflow as tf    
from tensorflow import keras 
from tensorflow.keras import Sequential, layers, callbacks
from tensorflow.keras.layers import Input, Reshape,Conv2D, MaxPooling2D, LSTM, Dense, Dropout, Flatten, Reshape, TimeDistributed
from keras import backend as K

from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, mean_absolute_percentage_error 

plt.rcParams['font.sans-serif'] = ['SimHei']     # 显示中文
plt.rcParams['axes.unicode_minus'] = False  # 显示负号
plt.rcParams.update({'font.size':18})  #统一字体字号

1.2 导入数据集

实验数据集采用数据集6：澳大利亚电力负荷与价格预测数据（下载链接），包括数据集包括日期、小时、干球温度、露点温度、湿球温度、湿度、电价、电力负荷特征，时间间隔30min。

# 导入数据
data_raw = pd.read_excel("澳大利亚电力负荷与价格预测数据.xlsx")
data_raw = data_raw[-365*24*6-49:-1].reset_index(drop=True)
data_raw

?对数据进行可视化

from itertools import cycle
# 可视化数据
def visualize_data(data, row, col):
    cycol = cycle('bgrcmk')
    cols = list(data.columns)
    fig, axes = plt.subplots(row, col, figsize=(16, 4))
    fig.tight_layout()
    if row == 1 and col == 1:  # 处理只有1行1列的情况
        axes = [axes]  # 转换为列表，方便统一处理
    for i, ax in enumerate(axes.flat):
        if i < len(cols):
            ax.plot(data.iloc[:,i], c=next(cycol))
            ax.set_title(cols[i])
        else:
            ax.axis('off')  # 如果数据列数小于子图数量，关闭多余的子图
    plt.subplots_adjust(hspace=0.6)
    plt.show()

visualize_data(data_raw.iloc[:,2:], 2, 3)

单独查看部分负荷数据。

# 预测结果可视化
plt.figure(dpi=100, figsize=(14, 4))
plt.plot(data_raw['电力负荷'], markevery=5)
plt.xlabel('时间')
plt.ylabel('负荷')
plt.show()

??

?1.3 缺失值分析

首先查看数据的信息，发现并没有缺失值

data_raw.info()

?进一步统计缺失值

data_raw.isnull().sum()

?2 构造训练数据

构造数据前先将数据变为数值类型

data = data_raw.iloc[:-48*8,2:].values

构造训练数据，也是真正预测未来的关键。首先设置预测的timesteps时间步、predict_steps预测的步长（预测的步长应该比总的预测步长小），length总的预测步长，参数可以根据需要更改。

timesteps = 48*7 #构造x，为48*7个数据,表示每次用前48*7个数据作为一段
predict_steps = 4 #构造y，为4个数据，表示用后4个数据作为一段
length = 48  #预测多步，预测48个数据
feature_num = 6 #特征的数量

通过前timesteps行历史数据预测后面predict_steps个数据，需要对数据集进行滚动划分（也就是前timesteps行的数据和后predict_steps行的数据训练，后面预测时就可通过timesteps行数据预测未来的predict_steps行数据）。这里需要注意的是，因为是多变量滚动预测多变量，特征就是标签（例如，前5行[干球温度、露点温度、湿球温度、电价、电力负荷]预测第6行[干球温度、露点温度、湿球温度、电价、电力负荷]，划分数据集时，就用前5行当做train_x，第6行作为train_y，此时的train_y有多列，而不是只有1列）。

# 构造数据集，用于真正预测未来数据
# 整体的思路也就是，前面通过前timesteps个数据训练后面的predict_steps个未来数据
# 预测时取出前timesteps个数据预测未来的predict_steps个未来数据。
def create_dataset(datasetx, datasety=None, timesteps=96*7, predict_size=12):
    datax = []  # 构造x
    datay = []  # 构造y
    for each in range(len(datasetx) - timesteps - predict_size):
        x = datasetx[each:each + timesteps]
        # 判断是否是单变量分解还是多变量分解
        if datasety is not None:
            y = datasety[each + timesteps:each + timesteps + predict_size]
        else:
            y = datasetx[each + timesteps:each + timesteps + predict_size]
        datax.append(x)
        datay.append(y)
    return datax, datay

??数据处理前，需要对数据进行归一化，按照上面的方法划分数据，这里返回划分的数据和归一化模型（单变量和多变量的归一化不同，多变量归一化需要将X和Y分开归一化，不然会出现信息泄露的问题），此时的归一化是相当于是单变量归一化，函数的定义如下：

# 数据归一化操作
def data_scaler(datax, datay=None, timesteps=36, predict_steps=6):
    # 数据归一化操作
    scaler1 = MinMaxScaler(feature_range=(0, 1))   
    datax = scaler1.fit_transform(datax)
    # 用前面的数据进行训练，留最后的数据进行预测
    # 判断是否是单变量分解还是多变量分解
    if datay is not None:
        scaler2 = MinMaxScaler(feature_range=(0, 1))
        datay = scaler2.fit_transform(datay)
        trainx, trainy = create_dataset(datax, datay, timesteps, predict_steps)
        trainx = np.array(trainx)
        trainy = np.array(trainy)
        return trainx, trainy, scaler1, scaler2
    else:
        trainx, trainy = create_dataset(datax, timesteps=timesteps, predict_size=predict_steps)
        trainx = np.array(trainx)
        trainy = np.array(trainy)
        return trainx, trainy, scaler1, None

然后分解的数据进行划分和归一化。

trainx, trainy, scalerx, scalery = data_scaler(data.reshape(-1, 1), timesteps=timesteps, predict_steps=predict_steps)

3 CNN_LSTM-Attention模型训练

3.1 CNN_LSTM模型?

CNN-LSTM 是一种结合了 CNN 特征提取能力与 LSTM 对时间序列长期记忆能力的混合神经网络。

CNN 主要由四个层级组成， 分别为输入层、 卷积层、 激活层（Relu 函数）和池化层。 每一层都会将数据处理之后送到下一层， 其中最重要的是卷积层， 这个层级起到的作用是将特征数据进行卷积计算， 将计算好的结果传到激活层， 激活函数对数据进行筛选。最后一层是 LSTM 层， 这一层是根据 CNN 处理后的特征数据，对其模型进行进一步的维度修偏， 权重修正等工作， 为下一步输出精度较高的预测值做好准备， 在 LSTM 训练的过程中， 由于其神经网络内部包括了输入、 遗忘和输出门， 通常的做法是通过增减遗忘门和输入门的个数， 来控制算法的精度。
?

来源：基于改进的 CNN-LSTM 短期风功率预测方法研究

对于输入到 CNN-LSTM 的数据，首先，经过 CNN 的卷积层对局部特征进行提取，将提取后的特征向量传递到池化层进行特征向量的下采样和数据体量的压缩。然后，将经过卷积层和池化层处理后的特征向量经过一个扁平层转化成一维向量输入到 LSTM 中， 每一层 LSTM 后加一个随机失活层以防止模型过拟合。

3.2 搭建Attention模型

参考文章：https://www.cnblogs.com/jiangxinyang/p/9367497.html

(1) Attention思想

深度学习里的Attention model其实模拟的是人脑的注意力模型，举个例子来说，当我们观赏一幅画时，虽然我们可以看到整幅画的全貌，但是在我们深入仔细地观察时，其实眼睛聚焦的就只有很小的一块，这个时候人的大脑主要关注在这一小块图案上，也就是说这个时候人脑对整幅图的关注并不是均衡的，是有一定的权重区分的。这就是深度学习里的Attention Model的核心思想。

(2) Encoder-Decoder框架

所谓encoder-decoder模型，又叫做编码-解码模型。这是一种应用于seq2seq问题的模型。seq2seq问题简单的说，就是根据一个输入序列x，来生成另一个输出序列y。Encoder-Decoder模型中的编码，就是将输入序列转化成一个固定长度的向量；解码，就是将之前生成的固定向量再转化成输出序列。

Encoder-Decoder（编码-解码）是深度学习中非常常见的一个模型框架，准确的说，Encoder-Decoder并不是一个具体的模型，而是一类框架。Encoder和Decoder部分可以是任意的文字，语音，图像，视频数据，模型可以采用CNN，RNN，BiRNN、LSTM、GRU等等。所以基于Encoder-Decoder，我们可以设计出各种各样的应用算法。

Encoder-Decoder框架可以看作是一种文本处理领域的研究模式，应用场景异常广泛，下图是文本处理领域里常用的Encoder-Decoder框架最抽象的一种表示：

(3) Attention模型?

????????在Encoder-Decoder框架中，在预测每一个yi时对应的语义编码c都是一样的，也就意味着序列X中点对输出Y中的每一个点的影响都是相同的。这样就会产生两个弊端：一是语义向量无法完全表示整个序列的信息，再者就是先输入的内容携带的信息会被后输入的信息稀释掉，或者说，被覆盖了。输入序列越长，这个现象就越严重。这就使得在解码的时候一开始就没有获得输入序列足够的信息， 那么解码的准确度自然也就要打个折扣了。

　　为了解决上面的弊端，就需要用到我们的Attention Model（注意力模型）来解决该问题。在机器翻译的时候，让生成词不是只能关注全局的语义编码向量c，而是增加了一个“注意力范围”，表示接下来输出词时候要重点关注输入序列中的哪些部分，然后根据关注的区域来产生下一个输出。模型结构如下：

关于模型的更多介绍可以查阅相关文献，下面给出Attention的代码

# CNN_LSTM_Attention模型
from keras.layers import multiply, Permute, RepeatVector, Multiply, Lambda
from keras.models import Model

def attention_function(inputs, single_attention_vector=False):    
    # 获取 inputs 的时间步数和特征维度   
    TimeSteps = K.int_shape(inputs)[1]
    input_dim = K.int_shape(inputs)[2]
 
    a = Permute((2, 1))(inputs)   #将 inputs 的维度进行转置，维度顺序变为 (特征维度, 时间步维度)  
    a = Dense(TimeSteps, activation='softmax')(a)  #经过全连接层
    
     # 如果为 True，单一注意力操作
    if single_attention_vector:
        a = Lambda(lambda x: K.mean(x, axis=1))(a)  #对第二个维度进行求平均，得到单一注意力向量
        a = RepeatVector(input_dim)(a)   # 将单一注意力向量进行复制，使其与 inputs 的维度一致  
         
    a_probs = Permute((2, 1))(a)    # 再次将注意力权重进行转置，维度顺序变为 (时间步维度, 特征维度)  
    output_attention_mul = Multiply()([inputs, a_probs])  # 使用 Multiply 层将 inputs 和注意力权重进行元素级乘法操作
    return output_attention_mul

3.3 搭建CNN_LSTM-Attention模型

首先搭建模型的常规操作，然后使用训练数据trainx和trainy进行训练，进行50个epochs的训练，每个batch包含64个样本（建议使用GPU进行训练，增加epochs）。

def CNN_LSTM_Attention_train(trainx, trainy, timesteps, feature_num, predict_steps):
    # 调用GPU加速
    gpus = tf.config.experimental.list_physical_devices(device_type='GPU')
    for gpu in gpus:
        tf.config.experimental.set_memory_growth(gpu, True)
    
    #搭建cnn模型
    inputs = Input(shape=(timesteps, feature_num))
    reshaped = Reshape((timesteps, feature_num, 1))(inputs)
    conv2d = Conv2D(filters=64, kernel_size=3, strides=1, padding="same", activation="relu")(reshaped)
    maxpool = MaxPooling2D(pool_size=2, strides=1, padding="same")(conv2d)
    dropout = Dropout(0.3)(maxpool)
    reshape2 = Reshape((timesteps, -1))(dropout)
    
    #搭建atttention模型
    attention_out = attention_function(reshape2)
    
    #搭建lstm模型
    lstm1 = LSTM(128, return_sequences=True, dropout=0.2)(attention_out)
    lstm2 = LSTM(128, return_sequences=False, dropout=0.2)(lstm1)  
    repeat_vector = RepeatVector(predict_steps)(lstm2) 
    outputs = TimeDistributed(Dense(feature_num))(repeat_vector)
    model = Model(inputs=inputs, outputs=outputs)
    model.compile(loss="mean_squared_error", optimizer="adam", metrics=['accuracy'])
    print(model.summary())
    model.fit(trainx, trainy, epochs=50, batch_size=128)
 
    return model

然后进行训练，将训练的模型、损失和训练时间保存。

# 模型训练
model = CNN_LSTM_Attention_train(trainx, trainy, timesteps, feature_num, predict_steps)
# 模型保存
model.save('cnn_lstm_attention.h5')

4 CNN_LSTM-Attention模型预测

4.1 分量预测

下面介绍文章中最重要，也是真正没有未来特征的情况下预测未来标签的方法。整体的思路也就是取出预测前48*7行数据预测未来的4行个数据，然后将4行数据添加进历史数据，再预测4行数据，滚动预测。因为每次只预测4行数据，但是我要预测48个数据，所以采用的就是循环预测12次的思路。

# #滚动predict
# #因为每次只能预测4行数据，但是我要预测48行数据，所以采用的就是循环预测的思路。
# #每次预测的4行数据，添加到数据集中充当预测x，然后在预测新的4行y，再添加到预测x列表中，如此往复,最终预测出48行。
def predict_using_LSTM(model, data, timesteps, predict_steps, length, feature_num, scaler):
    # 初始化预测输入和输出
    predict_xlist = np.array(data).reshape(1, timesteps, feature_num) 
    predict_y = np.array([]).reshape(0, feature_num)

    while len(predict_y) < length:
        # 从最新的predict_xlist取出timesteps个数据，预测新的predict_steps个数据
        predictx = predict_xlist[:,-timesteps:,:]

        # 预测新值
        lstm_predict = model.predict(predictx)
        
        # 将新预测出来的predict_steps个数据，加入predict_xlist列表，用于下次预测
        predict_xlist = np.concatenate((predict_xlist, lstm_predict), axis=1)
        
        # 预测的结果y，每次预测的predict_steps个数据，添加进去，直到预测length个为止
        lstm_predict = scaler.inverse_transform(lstm_predict.reshape(-1, feature_num))
        predict_y = np.concatenate((predict_y, lstm_predict), axis=0)
        
    return predict_y

然后对数据进行预测，得到预测结果。

from tensorflow.keras.models import load_model

# 加载模型
model = load_model('CNN_LSTM_Attention.h5')
pre_x = scalerx.transform(data[-48*8:-48])
y_true = data[-48:, -1]
# 预测
y_predict = predict_using_LSTM(model, pre_x, timesteps, predict_steps, length, feature_num, scalerx)

4.2 可视化

对预测的结果进行可视化并计算误差。

# 预测并计算误差和可视化
def error_and_plot(y_true,y_predict):
    # 计算误差
    r2 = r2_score(y_true, y_predict)
    rmse = mean_squared_error(y_true, y_predict, squared=False)
    mae = mean_absolute_error(y_true, y_predict)
    mape = mean_absolute_percentage_error(y_true, y_predict)
    print("r2: %.2f\nrmse: %.2f\nmae: %.2f\nmape: %.2f" % (r2, rmse, mae, mape))
    
    # 预测结果可视化
    cycol = cycle('bgrcmk')
    plt.figure(dpi=100, figsize=(14, 5))
    plt.plot(y_true, c=next(cycol), markevery=5)
    plt.plot(y_predict, c=next(cycol), markevery=5)
    plt.legend(['y_true', 'y_predict'])
    plt.xlabel('时间')
    plt.ylabel('功率(kW)')
    plt.show()   
    
    return 0

error_and_plot(y_true,y_predict[:,-1])

?