机器学习实战之风功率预测
目前机器学习与气象数据的结合已经在实际生产中有了应用，比如风电场风功率预测、光伏功率预测和负荷预测。本文以风功率预测作为一个小栗子： 风功率预测是指以风电场的历史功率、历史风速、地形地貌、数值天气预报、风电机组运行状态等数据建立风电场输出功率的预测模型，以风速、功率或数值天气预报数据作为模型的输入，结合风电场机组的运行状态及运行工况，得到风电场未来的输出功率，预测时间尺度包括短期预测和超短期预测，目的是上报国家电网，利于国家电网调度。目前主流方案是结合数值天气预报和机器学习算法（LSTM、SVM等）对风功率进行时序预测，包含超短期预报（未来4个小时共16个时刻）和短期预报（未来三天共288个时刻）。 本文主要利用WRF的气象要素预报数据和LSTM算法进行风功率预测。
下面对建模过程进行介绍。
加载所需的Python库
import pandas as pd
from datetime import datetime
from math import sqrt
from numpy import concatenate
from matplotlib import pyplot
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import mean_squared_error
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
读取数据
风机轮毂高度处风速、风向等来自WRF模式预报，实际功率来自风场实测。
file_name='清洗完毕的风场数据.xls'
need=pd.DataFrame(pd.read_excel(file_name))#
need.columns
Index(['real power', 'wind speed', 'wind direction', 'pressure', 'temperature', 'humidity', 'air density'], dtype='object')
定义函数，将各列数据转换为监督学习形式
#参数:
#data: list or NumPy array.
#n_in: Number of lag observations as input (X).
#n_out: Number of observations as output (y).
#dropnan: Boolean whether or not to drop rows with NaN values.
#Returns:
#Pandas DataFrame of series framed for supervised learning.
def series_to_supervised(data, n_in=, n_out=, dropnan=True):
n_vars = if type(data) is list else data.shape[]
df = DataFrame(data)
cols, names = list(), list()
# input sequence (t-n, ... t-1)
for i in range(n_in, , -1):
cols.append(df.shift(i))
names += [('var%d(t-%d)' % (j+, i)) for j in range(n_vars)]
# forecast sequence (t, t+1, ... t+n)
for i in range(, n_out):
cols.append(df.shift(-i))
if i == :
names += [('var%d(t)' % (j+)) for j in range(n_vars)]
else:
names += [('var%d(t+%d)' % (j+, i)) for j in range(n_vars)]
# put it all together
agg = concat(cols, axis=)
agg.columns = names
# drop rows with NaN values
if dropnan:
agg.dropna(inplace=True)
return agg
数据归一化
values = need.values
# 确保数据为32位浮点型
values = values.astype('float32')
#开始归一化
scaler = MinMaxScaler(feature_range=(, ))
scaled = scaler.fit_transform(values)
将归一化后的数据使用series_to_supervised处理
n_in=3
n_out=2
n_feature=need.columns.size #特征数量，包含实际风功率
reframed = series_to_supervised(scaled,n_in, n_out)
print(reframed.column）Index(['var1(t-3)', 'var2(t-3)', 'var3(t-3)', 'var4(t-3)', 'var5(t-3)', 'var6(t-3)', 'var7(t-3)', 'var1(t-2)', 'var2(t-2)', 'var3(t-2)', 'var4(t-2)', 'var5(t-2)', 'var6(t-2)', 'var7(t-2)', 'var1(t-1)', 'var2(t-1)', 'var3(t-1)', 'var4(t-1)', 'var5(t-1)', 'var6(t-1)', 'var7(t-1)', 'var1(t)', 'var2(t)', 'var3(t)', 'var4(t)', 'var5(t)', 'var6(t)', 'var7(t)', 'var1(t+1)', 'var2(t+1)', 'var3(t+1)', 'var4(t+1)', 'var5(t+1)', 'var6(t+1)', 'var7(t+1)'], dtype='object')
根据实际需要，选取需要的变量，这里选取的预测目标为var1(t)，即预测当前t时刻的风功率.
reframed_selected=reframed [ ['var2(t-3)', 'var3(t-3)', 'var4(t-3)', 'var5(t-3)',
'var6(t-3)', 'var7(t-3)', 'var2(t-2)', 'var3(t-2)',
'var4(t-2)', 'var5(t-2)', 'var6(t-2)', 'var7(t-2)',
'var2(t-1)', 'var3(t-1)', 'var4(t-1)', 'var5(t-1)', 'var6(t-1)',
'var7(t-1)', 'var2(t)', 'var3(t)', 'var4(t)', 'var5(t)',
'var6(t)', 'var7(t)', 'var1(t)']]
分割数据
# 选取前90%的数据作为训练集，余下部分为测试集
values = reframed_selected.values #通过dataframs的values方法获取值
n_train_15min=int(values.shape[0]*0.9)# 行数*90%得到索引位置
train = values[:n_train_15min, :]
test = values[n_train_15min:, :]
# 再分别将测试集和训练集 按列分割为 test_X, test_y, tran_X, train_y
train_X, train_y = train[:, :-1], train[:, -1]# 倒数第一列为预测目标y
test_X, test_y = test[:, :-1], test[:, -1]
将输入变形为三维 [样本量, 时间步长, 特征量]
#将n_feature 减去1 是因为 tran_X 和test_X中没有包含var1，即没有包含实际功率这一特征
#步长为n_in+1=4 ，包含t-1、t-2、t-3、t 共四步，也就是说，本文用t-1、t-2、t-3、t 时刻的气象预报数据来预测t时刻的风功率
train_X = train_X.reshape(train_X.shape[], n_in+, n_feature-)
test_X = test_X.reshape(test_X.shape[], n_in+, n_feature-)
设计keras神经网络模型
model = Sequential()#创建序贯模型
#增添LSTM层，指定输入大小(train_X.shape[1], train_X.shape[2])
model.add(LSTM(64, input_shape=(train_X.shape[1], train_X.shape[2]))) #n_units=64
model.add(Dense(1)) # train_y ,test_y为一维，所以 Dense(1)
model.compile(loss='mse', optimizer='adam')#编译,指定损失函数和优化器
# 训练并拟合
history = model.fit(train_X, train_y, epochs=30, batch_size=300, validation_data=(test_X, test_y), verbose=2, shuffle=False)
#画出训练过程
pyplot.plot(history.history['loss'], label='train')
pyplot.plot(history.history['val_loss'], label='test')
pyplot.legend()
pyplot.show()
结果如下：
预测
训练好模型后进行预测，运行此代码块若报错，可直接运行 yhat = model.predict(test_X)
import tensorflow as tf
global graph,model
graph = tf.get_default_graph()
with graph.as_default():
yhat = model.predict(test_X)
将三维数据再转变为二维
test_X_2d= test_X.reshape((test_X.shape[],(n_in+)*(n_feature-)))
将预测出的yhat与 test_x_2d的后六列合并，也就是 'var2(t)', 'var3(t)','var4(t)', 'var5(t)','var6(t)', 'var7(t)'合并 ，这样inv_yhat就成了7维，与当初归一化时values变量维度一样，且特征对应，从而能够顺利反归一化。也就是说归一化与反归一化时候，列数要保持一致，如果不一致则会报错。
inv_yhat = concatenate((yhat, test_X_2d[:, -(n_feature-1):]), axis=)
inv_yhat = scaler.inverse_transform(inv_yhat)
inv_yhat = inv_yhat[:,]#挑选出预测出的风功率，即第一列
这里也对测试集进行反归一化，只不过y是真实值，不是预测值
test_y = test_y.reshape((len(test_y), 1))
inv_y = concatenate((test_y, test_X_2d[:, -(n_feature-1):]), axis=1)
inv_y = scaler.inverse_transform(inv_y) #w
inv_y = inv_y[:,0] #挑选出反归一化后真实的风功率
评估
#画出图像
#rmse = sqrt(mean_squared_error(inv_y, inv_yhat))
#print('Test RMSE: %.3f' % rmse)
pyplot.plot(inv_y,label='actual')
pyplot.plot(inv_yhat, label='model')
pyplot.legend()
pyplot.show()
此栗子只是一个线下的基线模型，大量参考了Jason Brownlee的博客开源代码，许多地方都可以改进，比如增加网络层数，改变网络结构（Bi-LSTM,CNN- LSTM）,添加dropout层,使用earlystop,使用动态学习率等。在调参时，推荐使用调参库hyperopt或者talos，本文通过hyperopt对epoches、n_units、batch_size三个参数调参，选择出最优组合，提高了2%的合格率。另外，特征工程也至关重要，本文并没有对原始数据进行特征工程处理。 在上线的过程中，若要及时预预报处多个风场的风功率，如何进行分布式运行等一系列问题需要去线上实验解决。
---------------------手动分割线------------------
keras的LSTM运行过程中可能会产生报错： TypeError: while_loop() got an unexpected keyword argument 'maximum_iterations' 解决方法 ： 使用 Keras 2.1.2， Tensorflow 1.4.1 版本, 可在window cmd中或者linux 命令行中输入： pip install keras==2.1.2 pip install tensorflow==1.4.0
本文代码主要参考了如下网络资源：
https://machinelearningmastery.com/multivariate-time-series-forecasting-lstms-keras/
https://cloud.tencent.com/developer/article/1041442
https://www.jqr.com/article/000438