基于多元LSTM模型的短时交通流量预测

正好这学期选修了一门交通流理论的课程，这里介绍一下我最后的课程论文，如何基于上一篇提到的在PeMS数据库下载到的数据，应用LSTM模型进行短时交通流量预测。

1. 数据准备

2019年9月1日至9月28日四周的交通流量数据，数据粒度是5分钟一统计，属性包含流率(veh/ 5 minutes),速度(mph),占有率(%),重车比例(%)。

2. 数据描述性统计

(1) 绘制流率、速度、密度三参数基本图

(2) 概率密度分布

绘制流率、速度的概率密度分布图和累积概率密度分布图。

(3) 时间序列变化

观察流量、重车分布、速度、占有率的时序变化。 以9月1日一天的数据为例：

所有四周的流量时间序列变化：

3. LSTM模型构建

以前三周数据训练，最后一周数据预测。时间步设为12，即前一个小时预测下一个5分钟。 话不多说，直接上代码。

import pandas as pd
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split
from keras.preprocessing.sequence import TimeseriesGenerator
from keras.models import Sequential
from keras.layers import LSTM, Dense, Activation
from sklearn.metrics import mean_absolute_error as mae
from sklearn.metrics import mean_squared_error as mse
from sklearn.metrics import r2_score
import matplotlib.pyplot as plt
import numpy as np

# 将序列转换为监督学习问题
def series_to_supervised(datasetX,datasetY,lag):
    print('x的属性数量',datasetX.shape[1])
    data_gen = TimeseriesGenerator(datasetX,datasetY, length=lag,
                                   batch_size=1)  # length是滑lag,sampling_rate是timesteps,batch_size是y的维度
    X_out, Y_out = data_gen[0]
    for j in range(0, len(data_gen)):
        batch = data_gen[j]
        x, y = batch
        X_out = np.concatenate((X_out, x))
        Y_out = np.concatenate((Y_out, y))
    X_out = np.delete(X_out, 0, axis=0)
    Y_out = np.delete(Y_out, 0)
    print('X shape',X_out.shape)
    print('Y shape', Y_out.shape)
    return X_out,Y_out

# 加载数据集
path='.\'
olddf=pd.read_excel(path+'PeMS数据.xlsx',sheet_name='Sheet1')
df=olddf.copy(deep=False)
df['Flow (Veh/5 Minutes)']=df['Flow (Veh/5 Minutes)'].map(lambda x:3*x)#转换为veh/h/ln

#生成lstmX,Y
input_attributes=['Flow (Veh/5 Minutes)','Occupancy (%)','Speed (mph)','Truck Prop (%)']
datasetX=df[input_attributes].values
datasetY=df['Flow (Veh/5 Minutes)']
scalerx = MinMaxScaler(feature_range=(0, 1))
scalery = MinMaxScaler(feature_range=(0, 1))
X_scaled = scalerx.fit_transform(datasetX)
Y_scaled= scalery.fit_transform(np.array(datasetY).reshape(-1, 1))
X,Y=series_to_supervised(X_scaled,Y_scaled,12)#12个数据点即1小时

#划分训练集与测试集
X_train,X_test,Y_train,Y_test= train_test_split(X,Y,test_size=0.25,random_state=0)
print('X_train shape',X_train.shape)
print('X_test shape',X_test.shape)
print('Y_train shape',Y_train.shape)
print('Y_test shape',Y_test.shape)


#构建模型
# input_dim是输入的train_x的最后一个维度，train_x的维度为(n_samples, time_steps, input_dim)
model = Sequential()
model.add(LSTM(30, input_shape=(X_train.shape[1], X_train.shape[2])))
model.add(Dense(1))
model.add(Activation('tanh'))
model.compile(loss='mse', optimizer='adam')

# 拟合神经网络模型
history = model.fit(X_train, Y_train, epochs=50, batch_size=64, validation_data=(X_train, Y_train), verbose=2,
                    shuffle=False)

# 绘制训练loss数据
fig1 = plt.figure()
ax1 = fig1.add_subplot(111)
ax1.plot(history.history['loss'], label='train')
ax1.plot(history.history['val_loss'], label='test')
ax1.set_ylabel('Loss')
ax1.set_xlabel('Epochs')
plt.legend()
plt.show()

# 做出预测
yhat = model.predict(X_test)
invY_predict=scalery.inverse_transform(np.array(yhat).reshape(-1,1))
invY_test=scalery.inverse_transform(np.array(Y_test).reshape(-1,1))

print('lstm结果')
print('r_square', r2_score(invY_test,invY_predict))
print('mae', mae(invY_test,invY_predict))
print('rmse', np.sqrt(mse(invY_test,invY_predict)))

fig2 = plt.figure(figsize=(10,7) , dpi=80)
ax2 = fig2.add_subplot(111)
ax2.plot(invY_test,color='blue', linestyle='-',label='Ground Truth')
ax2.plot(invY_predict,color='red', linestyle='--',label='LSTM Prediciton')
ax2.set_ylabel('Flow (veh/h/ln)')
ax2.set_xlabel('Date')
xlabel=[]
for i in range(22,30):
    xlabel.append('9/%d'%(i))
ax2.set_xticks(range(0,2304,288))
ax2.set_xticklabels(xlabel,rotation=90)
plt.legend()
plt.show()

#绘制9/28一天的图
fig3 = plt.figure(figsize=(10,7) , dpi=80)
ax3 = fig3.add_subplot(111)
ax3.plot(invY_test[-288:],color='green', linestyle='-',label='Ground Truth')
ax3.plot(invY_predict[-288:],color='purple', linestyle='--',label='LSTM Prediciton')
ax3.set_ylabel('Flow (veh/h/ln)')
ax3.set_xlabel('Time')
xlabel=[]
for i in range(0,25):
    xlabel.append('%d:00'%(i))
ax3.set_xticks(range(0,289,12))
ax3.set_xticklabels(xlabel,rotation=90)
plt.legend()
plt.show()

结果展示：

最后一周所有的预测：

9月28日结果： 预测的效果还可以，毕竟只设置了单隐层。