Python数据分析案例52——基于SSA-LSTM的风速预测(麻雀优化)

案例背景

又要开始更新时间序列水论文的系列的方法了，前面基于各种不同神经网络层，还有注意力机制做了一些缝合模型。
其实论文里面用的多的可能是优化算法和模态分解，这两个我还没出专门的例子，这几天正好出一个优化算法的例子来做一个时间序列模型的缝合版。
想看更多的发论文用的模型可以参考我数据分析案例之前的文章，或者关注我后面的文章。

其实优化算法在python里面的生态不如MATLAB，现有的包很少，所以都是现写的。我自己也有优化算法专栏，以后有机会都写上去。本次的Python版的麻雀算法就是手写的，网上基本没有。

本次就简单点，使用优化算法里面表现较好的麻雀优化算法，优化算法我也做过一些测试，虽然都是各有优势，但是从通用性和整体表现来看，麻雀优化算法表现是较好的，那些什么混沌麻雀，自适应麻雀也差不多，可能在特殊的情况下表现会好一些。什么，你问问为什么不用粒子群，退火，遗传这种算法？emmmm，你自己去找些函数试试就知道他们比麻雀算法差多少了。。。

数据介绍

本次数据集有两个csv，一个桥面风速，一个气象站风速。

，

一般来说，桥面风速是好测量的，气象站的风速是被认为是真实的风速，所以我们当前的用气象站的风速作为y，之前的桥面风速和气象站风速作为X。

当然，需要本次案例演示数据和全部代码文件的同学还是可以参考：风速预测。

代码实现

导入包，深度学习的包有点多、

import os
import math
import time
import datetime
import random as rn
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams ['font.sans-serif'] ='SimHei'               #显示中文
plt.rcParams ['axes.unicode_minus']=False               #显示负号from sklearn.model_selection import train_test_split
from sklearn.preprocessing import MinMaxScaler,StandardScaler
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error,r2_scoreimport tensorflow as tf
import keras
from keras.layers import Layer
import keras.backend as K
from keras.models import Model, Sequential
from keras.layers import GRU, Dense,Conv1D, MaxPooling1D,GlobalMaxPooling1D,Embedding,Dropout,Flatten,SimpleRNN,LSTM
from keras.callbacks import EarlyStopping
#from tensorflow.keras import regularizers
#from keras.utils.np_utils import to_categorical
from tensorflow.keras  import optimizers

读取数据，展示前五行

data0=pd.concat([pd.read_excel('bridge2.xlsx').set_index('时间'),pd.read_excel('weather_station.xlsx').set_index('时间')],axis=1).sort_index().fillna(0)
data0.head()

一行代码读取两个文件，并且合并，代码风格还是简洁优雅的。

注意想换成自己的数据集的话，要预测的y放在最后一列。

构建训练集和测试集

时间序列的预测一套滑动窗口，构建的函数如下

def build_sequences(text, window_size=24):#text:list of capacityx, y = [],[]for i in range(len(text) - window_size):sequence = text[i:i+window_size]target = text[i+window_size]x.append(sequence)y.append(target)return np.array(x), np.array(y)def get_traintest(data,train_ratio=0.8,window_size=24):train_size=int(len(data)*train_ratio)train=data[:train_size]test=data[train_size-window_size:]X_train,y_train=build_sequences(train,window_size=window_size)X_test,y_test=build_sequences(test,window_size=window_size)return X_train,y_train[:,-1],X_test,y_test[:,-1]

对x和y进行标准化

data=data0.to_numpy()
scaler = MinMaxScaler() 
scaler = scaler.fit(data[:,:-1])
X=scaler.transform(data[:,:-1])   y_scaler = MinMaxScaler() 
y_scaler = y_scaler.fit(data[:,-1].reshape(-1,1))
y=y_scaler.transform(data[:,-1].reshape(-1,1))

划分训练集和测试集

train_ratio=0.8     #训练集比例   
window_size=64      #滑动窗口大小，即循环神经网络的时间步长
X_train,y_train,X_test,y_test=get_traintest(np.c_[X,y],window_size=window_size,train_ratio=train_ratio)
print(X_train.shape,y_train.shape,X_test.shape,y_test.shape)

数据可视化

y_test = y_scaler.inverse_transform(y_test.reshape(-1,1))
test_size=int(len(data)*(1-train_ratio))
plt.figure(figsize=(10,5),dpi=256)
plt.plot(data0.index[:-test_size],data0.iloc[:,-1].iloc[:-test_size],label='Train',color='#FA9905')
plt.plot(data0.index[-test_size:],data0.iloc[:,-1].iloc[-(test_size):],label='Test',color='#FB8498',linestyle='dashed')
plt.legend()
plt.ylabel('Predict Series',fontsize=16)
plt.xlabel('Time',fontsize=16)
plt.show()

训练函数的准备

下面继续自定义函数，评价指标

def set_my_seed():os.environ['PYTHONHASHSEED'] = '0'np.random.seed(1)rn.seed(12345)tf.random.set_seed(123)def evaluation(y_test, y_predict):mae = mean_absolute_error(y_test, y_predict)mse = mean_squared_error(y_test, y_predict)rmse = np.sqrt(mean_squared_error(y_test, y_predict))mape=(abs(y_predict -y_test)/ y_test).mean()#r_2=r2_score(y_test, y_predict)return mse, rmse, mae, mape #r_2

我们使用回归问题常用的mse, rmse, mae, mape作为预测效果的评价指标。

自定义注意力机制的类

class AttentionLayer(Layer):    #自定义注意力层def __init__(self, **kwargs):super(AttentionLayer, self).__init__(**kwargs)def build(self, input_shape):self.W = self.add_weight(name='attention_weight',shape=(input_shape[-1], input_shape[-1]),initializer='random_normal',trainable=True)self.b = self.add_weight(name='attention_bias',shape=(input_shape[1], input_shape[-1]),initializer='zeros',trainable=True)super(AttentionLayer, self).build(input_shape)def call(self, x):# Applying a simpler attention mechanisme = K.tanh(K.dot(x, self.W) + self.b)a = K.softmax(e, axis=1)output = x * areturn outputdef compute_output_shape(self, input_shape):return input_shape

自定义模型的构建

def build_model(X_train,mode='LSTM',hidden_dim=[32,16]):set_my_seed()model = Sequential()if mode=='MLP':model.add(Dense(hidden_dim[0],activation='relu',input_shape=(X_train.shape[-2],X_train.shape[-1])))model.add(Flatten())model.add(Dense(hidden_dim[1],activation='relu'))elif mode=='LSTM':# LSTMmodel.add(LSTM(hidden_dim[0],return_sequences=True, input_shape=(X_train.shape[-2],X_train.shape[-1])))#model.add(LSTM(hidden_dim[1]))#model.add(Flatten())#model.add(Dense(hidden_dim[1], activation='relu'))elif mode=='GRU':#GRUmodel.add(GRU(hidden_dim[0],return_sequences=True, input_shape=(X_train.shape[-2],X_train.shape[-1])))model.add(GRU(hidden_dim[1]))elif mode == 'Attention-LSTM':model.add(LSTM(hidden_dim[0], return_sequences=True, input_shape=(X_train.shape[-2], X_train.shape[-1])))model.add(AttentionLayer())        #model.add(LSTM(hidden_dim[1]))#, return_sequences=Falsemodel.add(Flatten())model.add(Dense(hidden_dim[1], activation='relu'))#model.add(Dense(4, activation='relu'))elif mode=='SSA-LSTM':# LSTMmodel.add(LSTM(hidden_dim[0],input_shape=(X_train.shape[-2],X_train.shape[-1])))#return_sequences=True, model.add(Dense(hidden_dim[1], activation='relu'))model.add(Dense(1))model.compile(optimizer='Adam', loss='mse',metrics=[tf.keras.metrics.RootMeanSquaredError(),"mape","mae"])return model

自定义画损失图函数和预测对比函数

def plot_loss(hist,imfname=''):plt.subplots(1,4,figsize=(16,2))for i,key in enumerate(hist.history.keys()):n=int(str('14')+str(i+1))plt.subplot(n)plt.plot(hist.history[key], 'k', label=f'Training {key}')plt.title(f'{imfname} Training {key}')plt.xlabel('Epochs')plt.ylabel(key)plt.legend()plt.tight_layout()plt.show()
def plot_fit(y_test, y_pred):plt.figure(figsize=(4,2))plt.plot(y_test, color="red", label="actual")plt.plot(y_pred, color="blue", label="predict")plt.title(f"拟合值和真实值对比")plt.xlabel("Time")plt.ylabel('power')plt.legend()plt.show()

可能有的小伙伴觉得看不懂了，没关系，我都是高度的封装，不需要知道每个函数里面的细节，大概知道他们是做什么的就行。因为下面要把他们全部打包为训练函数，改一下参数就可以使用不同的模型，很方便，

df_eval_all=pd.DataFrame(columns=['MSE','RMSE','MAE','MAPE'])
df_preds_all=pd.DataFrame()
def train_fuc(mode='LSTM',batch_size=64,epochs=30,hidden_dim=[32,16],verbose=0,show_loss=True,show_fit=True):#构建模型s = time.time()set_my_seed()model=build_model(X_train=X_train,mode=mode,hidden_dim=hidden_dim)earlystop = EarlyStopping(monitor='loss', min_delta=0, patience=5)hist=model.fit(X_train, y_train,batch_size=batch_size,epochs=epochs,callbacks=[earlystop],verbose=verbose)if show_loss:plot_loss(hist)#预测y_pred = model.predict(X_test)y_pred = y_scaler.inverse_transform(y_pred)#print(f'真实y的形状：{y_test.shape},预测y的形状：{y_pred.shape}')if show_fit:plot_fit(y_test, y_pred)e=time.time()print(f"运行时间为{round(e-s,3)}")df_preds_all[mode]=y_pred.reshape(-1,)s=list(evaluation(y_test, y_pred))df_eval_all.loc[f'{mode}',:]=ss=[round(i,3) for i in s]print(f'{mode}的预测效果为：MSE:{s[0]},RMSE:{s[1]},MAE:{s[2]},MAPE:{s[3]}')print("=======================================运行结束==========================================")return s[0]

所有的函数都准备完了，下面初始化参数，开始准备训练模型

window_size=64
batch_size=64
epochs=30
hidden_dim=[32,16]verbose=0
show_fit=True
show_loss=True
mode='LSTM'  #MLP,GRU

MLP模型训练

train_fuc(mode='MLP',batch_size=batch_size,epochs=epochs,hidden_dim=hidden_dim,verbose=1)

可以看到这个训练函数运行完后，可以清晰的看到每个训练轮的损失，损失的变化图，预测的效果对比图，还有评价指标的计算结果。

换模型也很便捷，只需要该mode这一个参数就行。

GRU模型训练

修改一下mode就行，其他参数你可以改也可以不改

train_fuc(mode='GRU',batch_size=batch_size,epochs=epochs,hidden_dim=hidden_dim,verbose=1)

LSTM训练

train_fuc(mode='LSTM',batch_size=batch_size,epochs=epochs,hidden_dim=hidden_dim,verbose=1)

Attention-LSTM模型的训练

train_fuc(mode='Attention-LSTM',batch_size=batch_size,epochs=epochs,hidden_dim=hidden_dim,verbose=1)

好像加了注意力机制的效果只变好了一点点。

麻雀搜索优化算法(SSA)

这里直接写上SSA的源代码，python版本的网上几乎是没有的

import numpy as np
import random
import copy
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D''' 种群初始化函数 '''
def initial(pop, dim, ub, lb):X = np.zeros([pop, dim])for i in range(pop):for j in range(dim):X[i, j] = random.random()*(ub[j] - lb[j]) + lb[j]return X,lb,ub'''边界检查函数'''
def BorderCheck(X,ub,lb,pop,dim):for i in range(pop):for j in range(dim):if X[i,j]>ub[j]:X[i,j] = ub[j]elif X[i,j]<lb[j]:X[i,j] = lb[j]return X'''计算适应度函数'''
def CaculateFitness(X,fun):pop = X.shape[0]fitness = np.zeros([pop, 1])for i in range(pop):fitness[i] = fun(X[i, :])return fitness'''适应度排序'''
def SortFitness(Fit):fitness = np.sort(Fit, axis=0)index = np.argsort(Fit, axis=0)return fitness,index'''根据适应度对位置进行排序'''
def SortPosition(X,index):Xnew = np.zeros(X.shape)for i in range(X.shape[0]):Xnew[i,:] = X[index[i],:]return Xnew'''麻雀发现者更新'''
def PDUpdate(X,PDNumber,ST,Max_iter,dim):X_new  = copy.copy(X)R2 = random.random()for j in range(PDNumber):if R2<ST:X_new[j,:] = X[j,:]*np.exp(-j/(np.random.random()*Max_iter))else:X_new[j,:] = X[j,:] + np.random.randn()*np.ones([1,dim])return X_new'''麻雀加入者更新'''            
def JDUpdate(X,PDNumber,pop,dim):X_new = copy.copy(X)for j in range(PDNumber+1,pop):if j>(pop - PDNumber)/2 + PDNumber:X_new[j,:]= np.random.randn()*np.exp((X[-1,:] - X[j,:])/j**2)else:#产生-1，1的随机数A = np.ones([dim,1])for a in range(dim):if(random.random()>0.5):A[a]=-1       AA = np.dot(A,np.linalg.inv(np.dot(A.T,A)))X_new[j,:]= X[1,:] + np.abs(X[j,:] - X[1,:])*AA.Treturn X_new                    '''危险更新'''   
def SDUpdate(X,pop,SDNumber,fitness,BestF):X_new = copy.copy(X)Temp = range(pop)RandIndex = random.sample(Temp, pop)SDchooseIndex = RandIndex[0:SDNumber]for j in range(SDNumber):if fitness[SDchooseIndex[j]]>BestF:X_new[SDchooseIndex[j],:] = X[0,:] + np.random.randn()*np.abs(X[SDchooseIndex[j],:] - X[1,:])elif fitness[SDchooseIndex[j]] == BestF:K = 2*random.random() - 1X_new[SDchooseIndex[j],:] = X[SDchooseIndex[j],:] + K*(np.abs( X[SDchooseIndex[j],:] - X[-1,:])/(fitness[SDchooseIndex[j]] - fitness[-1] + 10E-8))return X_new'''麻雀搜索算法'''
def SSA(pop,dim,lb,ub,Max_iter,fun):ST = 0.6 #预警值PD = 0.7 #发现者的比列，剩下的是加入者SD = 0.2 #意识到有危险麻雀的比重PDNumber = int(pop*PD) #发现者数量SDNumber = int(pop*SD) #意识到有危险麻雀数量X,lb,ub = initial(pop, dim, ub, lb) #初始化种群fitness = CaculateFitness(X,fun) #计算适应度值fitness,sortIndex = SortFitness(fitness) #对适应度值排序X = SortPosition(X,sortIndex) #种群排序GbestScore = copy.copy(fitness[0])GbestPositon = np.zeros([1,dim])GbestPositon[0,:] = copy.copy(X[0,:])Curve = np.zeros([Max_iter,1])for i in range(Max_iter):BestF = fitness[0]X = PDUpdate(X,PDNumber,ST,Max_iter,dim)#发现者更新X = JDUpdate(X,PDNumber,pop,dim) #加入者更新X = SDUpdate(X,pop,SDNumber,fitness,BestF) #危险更新X = BorderCheck(X,ub,lb,pop,dim) #边界检测fitness = CaculateFitness(X,fun) #计算适应度值fitness,sortIndex = SortFitness(fitness) #对适应度值排序X = SortPosition(X,sortIndex) #种群排序if(fitness[0]<=GbestScore): #更新全局最优GbestScore = copy.copy(fitness[0])GbestPositon[0,:] = copy.copy(X[0,:])Curve[i] = GbestScorereturn GbestScore,GbestPositon,Curve

这个模块可以写外面，从工程的角度来看放在一个py里面导入是最合理的。但是我们这是案例，考虑到简洁和易学性，所以我们都放在一个文件里面了。

优化算法定义完成后，定义目标函数

#import SSA
def fobj(X):s=train_fuc(mode='SSA-LSTM',batch_size=int(X[0]),epochs=int(X[1]),hidden_dim=[int(X[2]),int(X[3])],verbose=0,show_loss=False,show_fit=False)return s

进行优化算法的训练：

GbestScore1,GbestPositon1,Curve1 = SSA(pop=2,dim=4,lb=[8,20,30,12],ub=[40,40,80,42],Max_iter=2,fun=fobj)

我这里由于时间问题，我种群数量pop只用了2个，一般是30个，迭代次数一般是100-200次，我就改了2次，因为新电脑的TensorFlow不支持GPU加速，算的太慢了.......就没去搜索那么多次，就填了个较小的数字做演示好了。

打印最优的参数解和最优的适应度值

print('最优适应度值：',GbestScore1)
GbestPositon1=[int(i)for i in GbestPositon1[0]]
print('最优解为：',GbestPositon1)

带入最优解：

train_fuc(mode='SSA-LSTM',batch_size=GbestPositon1[0],epochs=GbestPositon1[1],hidden_dim=[GbestPositon1[2],GbestPositon1[3]],show_loss=True,show_fit=True,verbose=1)

虽然没搜索几次，但是这个效果还是不错的。

查看评价指标对比

好了，所有的模型都训练和预测了，评价指标都算完了，我们当然想对比了，我前面写训练函数都已经留了一手，预测的结果和效果都存下来了，和我一样一步步运行下来的可以直接查看预测效果。

df_eval_all

还是不直观，画图吧

bar_width = 0.4
colors=['c', 'orange','g', 'tomato','b', 'm', 'y', 'lime', 'k','orange','pink','grey','tan']
fig, ax = plt.subplots(2,2,figsize=(8,6),dpi=128)
for i,col in enumerate(df_eval_all.columns):n=int(str('22')+str(i+1))plt.subplot(n)df_col=df_eval_all[col]m =np.arange(len(df_col))plt.bar(x=m,height=df_col.to_numpy(),width=bar_width,color=colors)##plt.xlabel('Methods',fontsize=12)names=df_col.indexplt.xticks(range(len(df_col)),names,fontsize=10)plt.xticks(rotation=40)plt.ylabel(col,fontsize=14)plt.tight_layout() 
#plt.savefig('柱状图.jpg',dpi=512)
plt.show()

可以清楚地看见，SSA-lstm的效果最好，其次是GRU，然后是LSTM和attention-lstm。

所以说优化算法还是有效的，

继续画图他们的预测效果对比图：

plt.figure(figsize=(10,5),dpi=256)
for i,col in enumerate(df_preds_all[['MLP','GRU','LSTM','Attention-LSTM','SSA-LSTM']].columns):plt.plot(data0.index[-test_size-1:],df_preds_all[col],label=col) # ,color=colors[i]plt.plot(data0.index[-test_size-1:],y_test.reshape(-1,),label='Actual',color='k',linestyle=':',lw=2)
plt.legend()
plt.ylabel('wind',fontsize=16)
plt.xlabel('Date',fontsize=16)
#plt.savefig('点估计线对比.jpg',dpi=256)
plt.show()

也可以从这个图清楚的看到预测效果对比

总结

在这个案例里面的，SSA-LSTM效果好于GRU好于LSTM和attention-LSTM，说明优化算的效果是可以的，当然同学们还有时间可以用SSA-GRU，SSA-attention-LSTM都去试试，，看谁的效果好。模型修改就该buildmodel这个函数，很简单的，搭积木，要什么层就写什么层的名字就行。
画个数据也是很容易实验的。

这是优化算法+神经网络的方法啦，修改不同的优化算法就用自己自定义的算法替换就行，我后面的优化算法的专栏可能也会更新的，最近也有粉丝问问能不能出一个VMD或者CEEMDAN这些模态分解的对比，有时间我都写出来，可以关注我后面的文章。

创作不易，看官觉得写得还不错的话点个关注和赞吧，本人会持续更新python数据分析领域的代码文章~(需要定制类似的代码可私信)