基于LSTM的光伏发电数据预测的实现

场景：

　　这几天在看数字孪生的论文，看到这篇《基于LSTM与迁移学习的光伏发电功率预测数字孪生模型》，打算复现一下，正好学习一下机器学习，之前没了解过。由于是小白，主要记录一下大概实现流程，之后有空再深入了解原理。所有python代码都在jupyter上运行。

解决：

　　读取csv文件：

# 读入文件
dataset = pd.read_csv('5-Site_1.csv')

　　看到数据内容为，其中参数为

　　时间； 当前平均相位；有功输送能量接收；有功功率； 风速；天气温度（摄氏度）；全球水平辐射；风向； 降雨量； 最大风速；气压； 冰雹堆积； 日射强度； 温度探头1； 温度探头2：

　　

 

 

 　　然后设定要预测的标签，这里选取半个小时之后，由于5min一条数据，即6条数据差。选取‘Active_Power’作为发电功率：

# 其中太阳辐射强度，温度，湿度，风速会影响光伏输出功率
# 有功功率为负数，说明是发电
# 倒序
dataset.sort_index(inplace=True)
# 预测半小时之后
pre_mins = 6
# 设置标签，标签即为目标，如预测的结果
dataset['label'] = dataset['Active_Power'].shift(-pre_mins)
# 选择列显示
dataset.loc[:,['timestamp','Active_Power','label']]

　　然后看一下：

 

 　　很明显看到有NaN值的存在，于是做数据预处理：空值处理，这里把空值填0：

# 数据预处理-1
# 为NaN的值填为0
dataset.fillna(0)

　　然后是相关性分析，并不是所有参数都对发电功率有影响的，理论上风速，日射强度，降雨量会有影响（风速主要影响云层移动），做相关性分析可以实现降维：

# 相关性分析\
X = dataset.loc[0:1000,'Wind_Speed']
Y = dataset.loc[0:1000,'Active_Power']
# print(X)
# print(Y)
# pandas包下的相关性分析
result_1 = X.corr(Y)
# print(result_1)
result_3 = dataset.corr()
result_3

　　看一下结果：

 

 　　我们选取和Active_Power相关性大于0.1的（这里的阈值是一个经验值）：

# 去除相关性小于0.1的
result_4 = result_3['Active_Power'][result_3['Active_Power'] > 0.1]
result_4 = result_4.drop('Active_Power')
result_4 = result_4.drop('Current_Phase_Average_Mean')
result_4 = result_4.drop('label')
print(type(result_4)) 
result_4

　　结果如下：

 

 　　也可以绘制热力图，直观的看相关性：

# 相关性绘图
# 热力图
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
# 显示汉字
mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False
result_5 = dataset.loc[:,['Active_Power','Global_Horizontal_Radiation','Temperature_Probe_1','Temperature_Probe_2']]
figure, ax = plt.subplots(figsize=(12, 12))
sns.heatmap(result_5.corr(), square=True, annot=True, ax=ax)

　　如图：

 

 　　接下来，还是数据处理。剔除一些异常数据。我是先绘制时间曲线，看那个参数异常数据明显，再进行剔除：

# 绘制相关数据的时间曲线，看是否有异常数据
fig = plt.figure(figsize=(20,8),dpi=80)
# 行，列，索引
plt.subplot(4,1,1)
X_1 = dataset.loc[0:1000,'timestamp']
Y_1 = dataset.loc[0:1000,'Active_Power']
plt.plot(X_1,Y_1)
plt.xlabel('时间')
plt.ylabel('输出功率')
plt.title('输出功率变化情况')

plt.subplot(4,1,2)
X_1 = dataset.loc[0:1000,'timestamp']
Y_1 = dataset.loc[0:1000,'Global_Horizontal_Radiation']
plt.plot(X_1,Y_1)
plt.xlabel('时间')
plt.ylabel('水平辐射')
plt.title('水平辐射变化情况')

plt.subplot(4,1,3)
X_1 = dataset.loc[0:1000,'timestamp']
Y_1 = dataset.loc[0:1000,'Temperature_Probe_1']
plt.plot(X_1,Y_1)
plt.xlabel('时间')
plt.ylabel('温度1')
plt.title('温度1变化情况')

plt.subplot(4,1,4)
X_1 = dataset.loc[0:1000,'timestamp']
Y_1 = dataset.loc[0:1000,'Temperature_Probe_2']
plt.plot(X_1,Y_1)
plt.xlabel('时间')
plt.ylabel('温度2')
plt.title('温度2变化情况')

plt.show()

　　时间曲线：

 

 　　此时我注意到，在夜晚时间中，是不用考虑数据功率的，于是我绘制一天的时间变化，选取白天06:00--18:00的数据：

# 显示一天中输出功率的大概时间
# 行，列，索引
plt.figure(figsize=(20,10),dpi=80)
X_2 = dataset.loc[0 + 3:288 + 3 :12,'timestamp']
Y_2 = dataset.loc[0 + 3:288 + 3 :12,'Active_Power']

plt.xticks(rotation = 45)
plt.plot(X_2,Y_2)
plt.xlabel('时间')
plt.ylabel('水平辐射')
plt.title('一天内水平辐射变化情况')
plt.show()

# 由图可知，晚上输出功率为0，选取每天06:00--18:00的数据
from datetime import datetime
dataset['datetime'] = dataset['timestamp'].apply(lambda x:datetime.strptime(x,"%Y-%m-%d %H:%M:%S").hour)
filred_data = dataset[(dataset['datetime'] >= 6) & (dataset['datetime'] <= 18 )]
filred_data = filred_data.fillna(0)

　　一天内的曲线图：

 

 　　最后去除不相关的数据，去除夜晚数据看一下：

# 去除不相关的列
filred_data = filred_data.drop('datetime',axis = 1)
filred_data = filred_data.drop('Active_Energy_Delivered_Received',axis = 1)
filred_data = filred_data.drop('Current_Phase_Average_Mean',axis = 1)
filred_data = filred_data.drop('Wind_Speed',axis = 1)
filred_data = filred_data.drop('Weather_Temperature_Celsius',axis = 1)
filred_data = filred_data.drop('Wind_Direction',axis = 1)
filred_data = filred_data.drop('Weather_Daily_Rainfall',axis = 1)
filred_data = filred_data.drop('Max_Wind_Speed',axis = 1)
filred_data = filred_data.drop('Air_Pressure',axis = 1)
filred_data = filred_data.drop('Hail_Accumulation',axis = 1)
filred_data = filred_data.drop('Pyranometer_1',axis = 1)
filred_data
# 调整序号
filred_data = filred_data.reset_index(drop=True)
filred_data

　　如下：

 

 　　这个时候再来查找异常值，异常值的定义有很多方式：3a法，Z_score法，箱线法等。这里采用3a法：

# 用3sigma法找出超出（u-3a,u+3a）的值
def three_sigma(ser):
    # 计算平均值
    mean_data = ser.mean()
    # 计算标准差
    std_data = ser.std();
    rule = (mean_data - 3 * std_data > ser) | (mean_data + 3 * std_data < ser)
    # 生成一个从0开始，到ser长度-1结束的连续索引，再根据rule列表中的True值，直接保留所有为True的索引，也就是异常值的行索引
    index = np.arange(ser.shape[0])[rule]
    # 获取异常值
    outliers = ser.iloc[index]
    return outliers

# 异常值处理
mean_data = filred_data['Temperature_Probe_1'].mean();
mean_data

　　异常值的处理有很多方式，插值替代，均值替代等。这里我选用的是均值替代：

# 由图可知，对温度1的数据进行异常值处理:均值替代
# 经常用的还有插值法
error_data = three_sigma(filred_data['Temperature_Probe_1'])
true_data = filred_data['Temperature_Probe_1']
print(error_data)
for i, v in error_data.items():
    true_data[i] = mean_data

plt.figure(figsize=(20,10),dpi=80)
X_1 = filred_data.loc[0:1000,'timestamp']
Y_1 = true_data[0:1001]
plt.plot(X_1,Y_1)
plt.xlabel('时间')
plt.ylabel('温度1')
plt.title('修正后温度1变化情况')
plt.show()

　　看一下修正后温度一的曲线，温度一应该是光伏面板表面温度：

 

 　　然后可以开始准备训练数据了。首先是标准化，主要是为了避免不同量纲的影响，可以采用0-1归一化：

# 标准化
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
# 注意左闭右开
# 归一化是为了避免不同量纲的影响
X_3 = scaler.fit_transform(filred_data.iloc[:,1:5])

　　然后划分训练集（80%），测试集（20%）以及其y值（就是预测的结果）：

# 0.8为训练集，0.2为测试集
train_size=int(len(X_3)*0.8)
test_size=int(len(X_3))-train_size
train_data = X_3[0:train_size,0:4]
# y值是对应的预测值
y_train_data = list(filred_data.iloc[0:train_size,5])
test_data = X_3[train_size:,0:4]
y_test_data = list(filred_data.iloc[train_size:,5])

　　然后做输入输出，输入是一个时间序列，这里取5，即5个连续时刻的参数作为输入，输出一个半小时后的预测值：

# LSTM模型构建,多输入单输出
# 取一定时间内的数据作为输入，预测该时间后一段时间的label，最后不断缩小差值loss
# 记忆的天数
men_days = 5

from collections import deque;
deq = deque(maxlen= men_days)

X_4 = []
for i in train_data:
    deq.append(list(i))
    if len(deq) == men_days:
        X_4.append(list((deq)))
X_Test = []
for i in test_data:
    deq.append(list(i))
    if len(deq) == men_days:
        X_Test.append(list((deq)))
        
X_4 = np.array(X_4)
X_Test = np.array(X_Test)

　　对齐一下array的长度：

y_train_data = np.array(y_train_data[:-men_days+1])
y_test_data = np.array(y_test_data[:])

　　开始创建并训练LSTM:

# 创建并训练LSTM
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM,Dense,Dropout

model = Sequential()
# 第一层，激活函数选择‘relu’
model.add(LSTM(units=32, input_shape=(X_4.shape[1], X_4.shape[2])))
model.add(Dropout(0.5))  
# 全连接层
model.add(Dense(1))  

# 模型编译
model.compile(loss='mae', optimizer='adam')  

# validation_data是测试的数据集
# epochs相当于迭代次数
history = model.fit(X_4,y_train_data,batch_size=20,epochs=20,validation_data=(X_Test, y_test_data), verbose=2,shuffle=False)
# val_mape是测试集上的数据
# mape是训练集的数据
# val_loss越小越好,尽量在0.02左右

　　这里看一下输出：

　　

 

 　　可以绘制损失图，这个值应该越小越好：

# 绘制损失图
plt.plot(history.history['loss'], label='train')
plt.plot(history.history['val_loss'], label='test')

plt.ylabel('loss', fontsize='10')
plt.xlabel('epoch', fontsize='10')
plt.legend()
plt.show()

　　如图：

 

 　　对比一下预测数据和实际数据，绘制对比图：

# 预测数据
from numpy import concatenate
print(X_Test)
y_predict = model.predict(X_Test)

# 绘制预测结果图
x_predict =  filred_data.loc[train_size:,'timestamp']

plt.figure(figsize=(20,10),dpi=80)
plt.plot(x_predict[0:1000],y_predict[0:1000],color='black', label='predict')

plt.plot(x_predict[0:1000],y_test_data[0:1000],color='red', label='true')

plt.xlabel('时间')
plt.ylabel('预测值')
plt.title('实际值和预测值的对比图')
plt.legend()
plt.show()

　　如下：

 

 　　由上可知，其实很难做到实时预测，应该预测本身带有滞后性。总之，还是了解了很多的，之后还要看一下强化学习的东西。

参考：

　　《基于LSTM与迁移学习的光伏发电功率预测数字孪生模型_史凯钰》

　　https://blog.csdn.net/qq_37511129/article/details/95983968

　　https://blog.csdn.net/akun1213/article/details/122676867

　　https://www.bilibili.com/video/BV1U341127k5?p=6&spm_id_from=pageDriver

　　https://blog.csdn.net/Jruo911/article/details/124182164

　　https://blog.csdn.net/heianzhongjinhua/article/details/122895038

　　数据集来自：https://dkasolarcentre.com.au/download?location=alice-springs