基于LSTM的光伏发电预测

本文是博主的《机器学习》课程的作业，很简单，没有涉及到什么十分高深的内容，最终的实现是调用了库。

我也是一个新手，希望大家多多包容~

数据我会后面上传，我直接上传原数据，处理数据的过程我感觉也是很有成就感的，大家自己处理数据试试吧！

数据下载地址：https://download.csdn.net/download/qq_37511129/11356993

数据处理

首先我们来看数据，训练集和测试集的标签并不统一，训练集中很多标签并没有出现在测试集中，所以我们把多余的特征删除。

剩下的标签有：

'Temp_Out','Out_Hum','Dew_Pt','Wind_Speed','Wind_Dir','Hi_Dir','Wind_Chill','Heat_Index','THW_Index','THSW_Index',’ Bar’,’Rain’,’Rain_Rate’,'Solar_Rad'。

通过整理，我们可以看出此次数据集总共有13个特征，给出的数据是每一个传感器在每一分钟记录的数据，总共是11个月的数据。数据量比较大，为了减少计算量，我们需要减少数据量，将数据集降维。

因为不是每一个特征都对实验结果有很大的影响，所以我们需要一个办法来确定相关程度，然后删除掉相关度不大的特征，实现数据降维。我选择的是皮尔森系数。

皮尔森相关系数是一种线性相关系数。皮尔森相关系数是用来反映两个变量线性相关程度的统计量。系数越大，说明相关性越强。

具体的计算结果我就不放出来了，皮尔森系数的计算方法大家也可以直接去网上查一下。

相关系数超过0.1的是前十个特征，所以我们就选取前十个特征进行分析。

我们来看数据分布。首先我定义了如下函数来读取数据，并绘制出各个数据分布。

def drow_pollution():
dataset = pd.read_csv('training_data_1.csv', header=0, index_col=(0,1))
values = dataset.values
# specify columns to plot
groups = [0, 1, 2, 3, 5, 6, 7, 8, 9]
i = 1
# plot each column
pyplot.figure(figsize=(10,10))
for group in groups:
pyplot.subplot(len(groups), 1, i)
pyplot.plot(values[:, group])
pyplot.title(dataset.columns[group], y=0.5, loc='right')
i += 1
pyplot.show()

绘制出的图像如下：

图表 1 原始数据分布图像

个分布图像中可以看出来，很多数据上存在着一些突变，这些突变应该是由于一些外部原因造成的异常值，我们需要处理掉这些异常值。

我主要选用的方法为3 法，先算出平均值，然后算出标准差，然后根据3来删除掉异常值，具体的代码如下：

import pandas as pd
dataset = pd.read_csv('training_data.csv', header=0, index_col=(0,1))
print(dataset.shape)
groups = [0, 1, 2, 3, 5, 6, 7, 8, 9]
pyplot.figure(figsize=(10,10))
for group in groups:
values = dataset.values
dataset=dataset[abs((values[:,group]-alues[:,group].mean())/values[:,group].std())<3]

删除掉异常值后的结果为：

图表 2 异常值处理过后的分布图像

值的分布图像和原分布图像进行对比我们可以明显看出之前的那些异常值已经处理掉了。

模型构建

从数据的结构来看，给出的数据是一个时间序列。所以我选择使用LSTM模型来实现数据预测的问题。

LSTM神经网络的特点是通过一个门来确定是完全保留和完全舍弃之前的信息。他可以将之前的时间序列也参考在模型之内。所以首先我们需要重新生成一个输入，我们需要定义一个函数，这个函数可以将之前n个时间序列都变成输入，具体实现如下：

def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
# convert series to supervised learning
n_vars = 1 if type(data) is list else data.shape[1]
df = pd.DataFrame(data)
cols, names = list(), list()
# input sequence (t-n, ... t-1)
for i in range(n_in, 0, -1):
cols.append(df.shift(i))
names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]
# forecast sequence (t, t+1, ... t+n)
for i in range(0, n_out):
cols.append(df.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)]
# put it all together
agg = pd.concat(cols, axis=1)
agg.columns = names
# drop rows with NaN values
if dropnan:
agg.dropna(inplace=True)
return agg

在有了这个序列转换的函数之后，我们就可以开始构造我们的数据了。首先我们要读入我们的training_data和test_data，运用一个库将数据进行了归一化处理。然后我设置的validation的比例是0.2也就是说training_data中百分之八十的数据是用来训练的，百分之二十的数据用来验证。

具体实现如下：

def train_test():
# 读取文件，包含头信息
training_data = pd.read_csv('training_data_1.csv')
test_data = pd.read_csv('test_data.csv')
# 删掉不用的字段
feature=['Temp_Out','Out_Hum','Dew_Pt','Wind_Speed','Wind_Dir','Hi_Dir','Wind_Chill','Heat_Index','THW_Index','THSW_Index','Solar_Rad']
training_data=training_data[feature]
test_data=test_data[feature]
ok_y = test_data['Solar_Rad']
# df转array
values_1=test_data.values
values = training_data.values
# 原始数据标准化
scaler = MinMaxScaler(feature_range=(0, 1))
scaled = scaler.fit_transform(values[:,:9])
scaled_1 = scaler.fit_transform(values_1[:,:9])
# 根据过去一小时预测当前
n_hours = 1
n_features = scaled.shape[1]
# 构造特征，过去三小时与当前数据集合
reframed = series_to_supervised(scaled, n_hours, 1)
reframed_1 = series_to_supervised(scaled_1, n_hours, 1)
values = reframed.values
values_1 = reframed_1.values
# 划分训练集与测试集
size = 0.8
num = int(values.shape[0]*size)
train = values
test = values_1
train_X, train_y = train[:num, :n_features*n_hours], train[:num, -1]
validation_X,validation_y=train[num:, :n_features*n_hours], train[num:, -1]
test_X, test_y = test[:, :n_features*n_hours], test[:, -1]
# reshape为 3D [samples, timesteps, features]，将n_hours看成n个独立的时间序列而不是一个整体的
train_X = train_X.reshape((train_X.shape[0], n_hours, n_features))
validation_X = validation_X.reshape((validation_X.shape[0], n_hours, n_features))
test_X = test_X.reshape((test_X.shape[0], n_hours, n_features))
return train_X,train_y,test_X,test_y,validation_X,validation_y,scaler,n_hours,n_features,ok_y

然后就是模型的训练了。模型训练我直接使用的keras库，首先是确定为一个Sequential模型，然后add LSTM层，Dense代表的是输出的维度，因为只有一个维度，所以设置为1.损失函数使用的是mae，优化器使用的是adam。

def fit_network_1(train_X,train_y,test_X,test_y,validation_X,validation_y,scaler,n_hours,n_features,ok_y):
model = Sequential()
model.add(LSTM(units=50, input_shape=(train_X.shape[1], train_X.shape[2])))
model.add(Dense(1))
#model.add(Activation('linear'))
model.compile(loss='mae', optimizer='adam')
# fit network
history = model.fit(train_X, train_y, epochs=20, batch_size=1440, validation_data=(validation_X,validation_y), verbose=2, shuffle=False)
# plot history
pyplot.plot(history.history['loss'], label='train')
pyplot.plot(history.history['val_loss'], label='test')
pyplot.legend()
pyplot.show()
# make a prediction
yhat = model.predict(test_X)
test_X = test_X.reshape((test_X.shape[0], n_hours*n_features))
# 将预测y与当前时间的x组合
inv_yhat0 = concatenate((test_X[:, -n_features:-1], yhat), axis=1)
inv_yhat1 = scaler.inverse_transform(inv_yhat0)
inv_yhat = inv_yhat1[:,-1]
# 将实际y与当前时间的x组合
test_y = test_y.reshape((len(test_y), 1))
inv_y0 = concatenate((test_X[:, -n_features:-1],test_y), axis=1)
print('inv_y0')
print(inv_y0.shape)
inv_y1 = scaler.inverse_transform(inv_y0)
inv_y = inv_y1[:,-1]
# calculate RMSE
rmse = sqrt(mean_squared_error(inv_y, inv_yhat))
print('Test RMSE: %.3f' % rmse)

运行的结果如下：

算出的RMSE为2.210（每一次的结果都会有一些偏差），但大体的结果基本稳定在这个值左右。