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_score
 
import 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 capacity
    x, 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 mechanism
        e = K.tanh(K.dot(x, self.W) + self.b)
        a = K.softmax(e, axis=1)
        output = x * a
        return output
 
    def 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':
        # LSTM
        model.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':
        #GRU
        model.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=False
        model.add(Flatten())
        model.add(Dense(hidden_dim[1], activation='relu'))
        #model.add(Dense(4, activation='relu'))
    elif mode=='SSA-LSTM':
        # LSTM
        model.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}',:]=s
    s=[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.T
           
    return 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() - 1
            X_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] = GbestScore
    
    return 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.index
    plt.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数据分析领域的代码文章~(需要定制类似的代码可私信)