目錄
1?數據處理
1.1 導入庫文件
1.2 導入數據集
1.3 缺失值分析
2 VMD經驗模態分解
3?構造訓練數據
4 LSTM模型訓練
5 預測
1?數據處理
1.1 導入庫文件
import time
import datetime
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sampen import sampen2 # sampen庫用于計算樣本熵
from vmdpy import VMD # VMD分解庫import tensorflow as tf
from sklearn.cluster import KMeans
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, mean_absolute_percentage_error
from sklearn.preprocessing import MinMaxScaler
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, Activation, Dropout, LSTM, GRU
from tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping# 忽略警告信息
import warnings
warnings.filterwarnings('ignore') 1.2 導入數據集
實驗數據集采用數據集8：新疆光伏風電數據集（下載鏈接），數據集包括組件溫度(℃) 、溫度(°)?? ?氣壓(hPa)、濕度(%)、總輻射(W/m2)、直射輻射(W/m2)、散射輻射(W/m2)、實際發電功率(mw)特征，時間間隔15min。對數據進行可視化：
# 導入數據
data_raw = pd.read_excel("E:\\課題\\08數據集\\新疆風電光伏數據\\光伏2019.xlsx")
data_rawfrom itertools import cycle
# 可視化數據
def visualize_data(data, row, col):cycol = cycle('bgrcmk')cols = list(data.columns)fig, axes = plt.subplots(row, col, figsize=(16, 4))fig.tight_layout()if row == 1 and col == 1: # 處理只有1行1列的情況axes = [axes] # 轉換為列表，方便統一處理for i, ax in enumerate(axes.flat):if i < len(cols):ax.plot(data.iloc[:,i], c=next(cycol))ax.set_title(cols[i])else:ax.axis('off') # 如果數據列數小于子圖數量，關閉多余的子圖plt.subplots_adjust(hspace=0.6)plt.show()visualize_data(data_raw.iloc[:,1:], 2, 4)
?單獨查看部分功率數據，發現有較強的規律性。
?因為只是單變量預測，只選取實際發電功率(mw)數據進行實驗：
1.3 缺失值分析
首先查看數據的信息，發現并沒有缺失值
data_raw.info()
?進一步統計缺失值
data_raw.isnull().sum()
2 VMD經驗模態分解
使用VMD將目標信號分解成若干個模態，進一步可視化分解結果
# VMD分解函數
# signal: 輸入信號
# alpha: 正則化參數
# tau: 時間尺度參數
# K: 分量數量
# DC: 是否包括直流分量
# init: 初始化方法
# tol: 收斂容限
# n_ite: 最大迭代次數
def vmd_decompose(series=None, alpha=2000, tau=0, K=7, DC=0, init=1, tol=1e-7, draw=True): # 得到 VMD 分解后的各個分量、分解后的信號和頻率imfs_vmd, imfs_hat, omega = VMD(series, alpha, tau, K, DC, init, tol) # 將 VMD 分解分量轉換為 DataFrame, 并重命名df_vmd = pd.DataFrame(imfs_vmd.T)df_vmd.columns = ['imf'+str(i) for i in range(K)]return df_vmddf_vmd = vmd_decompose(data_raw['實際發電功率(mw)']) # 對 df_raw_data['AQI'] 進行 VMD 分解，并將結果賦值給 df_vmd
# 繪制 df_vmd 的數據，以子圖形式顯示每個分量
ax = df_vmd.plot(title='VMD Decomposition', figsize=(16,8), subplots=True,fontsize=16)
for a in ax:a.legend(loc='upper right',prop={'size': 14})plt.subplots_adjust(hspace=0.5)
將原始數據和分解后的模態合并
df_vmd['sum'] = data_raw['實際發電功率(mw)'] # 將 data_raw['實際發電功率(mw)']添加到 df_vmd 中的 'sum' 列
?這里利用VMD-LSTM進行預測的思路是通過VMD將原始功率分解為多個變量，然后將分解變量作為輸入特征，將原始出力功率作為標簽，將單變量轉為多變量進行預測。
3?構造訓練數據
構造訓練數據，也是真正預測未來的關鍵。首先設置預測的timesteps時間步、predict_steps預測的步長（預測的步長應該比總的預測步長小），length總的預測步長，參數可以根據需要更改。
timesteps = 96*5 #構造x，為96*5個數據,表示每次用前96*5個數據作為一段
predict_steps = 96 #構造y，為96個數據，表示用后96個數據作為一段
length = 96 #預測多步，預測96個數據
feature_num = 7 #特征的數量通過前5天的timesteps數據預測后一天的數據predict_steps個，需要對數據集進行滾動劃分（也就是前timesteps行的特征和后predict_steps行的標簽訓練，后面預測時就可通過timesteps行特征預測未來的predict_steps個標簽）。因為是多變量，特征和標簽分開劃分，不然后面歸一化會有信息泄露的問題。
# 構造數據集，用于真正預測未來數據
# 整體的思路也就是，前面通過前timesteps個數據訓練后面的predict_steps個未來數據
# 預測時取出前timesteps個數據預測未來的predict_steps個未來數據。
def create_dataset(datasetx,datasety,timesteps=36,predict_size=6):datax=[]#構造xdatay=[]#構造yfor each in range(len(datasetx)-timesteps - predict_steps):x = datasetx[each:each+timesteps]y = datasety[each+timesteps:each+timesteps+predict_steps]datax.append(x)datay.append(y)return datax, datay
數據處理前，需要對數據進行歸一化，按照上面的方法劃分數據，這里返回劃分的數據和歸一化模型，函數的定義如下：
# 數據歸一化操作
def data_scaler(datax,datay):# 數據歸一化操作scaler1 = MinMaxScaler(feature_range=(0,1))scaler2 = MinMaxScaler(feature_range=(0,1))datax = scaler1.fit_transform(datax)datay = scaler2.fit_transform(datay)# 用前面的數據進行訓練，留最后的數據進行預測trainx, trainy = create_dataset(datax[:-timesteps-predict_steps,:],datay[:-timesteps-predict_steps,0],timesteps, predict_steps)trainx = np.array(trainx)trainy = np.array(trainy)return trainx, trainy, scaler1, scaler2然后對數據按照上面的函數進行劃分和歸一化。通過前5天的96*5數據預測后一天的數據96個，需要對數據集進行滾動劃分（也就是前96*5行的特征和后96行的標簽訓練，后面預測時就可通過96*5行特征預測未來的96個標簽）
datax = df_vmd[:,:-1]
datay = df_vmd[:,-1].reshape(df_vmd.shape[0],1)
trainx, trainy, scaler1, scaler2 = data_scaler(datax, datay)
4 LSTM模型訓練
首先搭建模型的常規操作，然后使用訓練數據trainx和trainy進行訓練，進行50個epochs的訓練，每個batch包含128個樣本（建議使用GPU進行訓練）。預測并計算誤差，訓練好將模型保存，并進行可視化，將這些步驟封裝為函數。
# # 創建lSTM模型
def LSTM_model_train(trainx, trainy):# 調用GPU加速gpus = tf.config.experimental.list_physical_devices(device_type='GPU')for gpu in gpus:tf.config.experimental.set_memory_growth(gpu, True)# LSTM網絡構建 start_time = datetime.datetime.now()model = Sequential()model.add(LSTM(128, input_shape=(timesteps, feature_num), return_sequences=True))model.add(Dropout(0.5))model.add(LSTM(128, return_sequences=True))model.add(LSTM(64, return_sequences=False))model.add(Dense(predict_steps))model.compile(loss="mean_squared_error", optimizer="adam")# 模型訓練model.fit(trainx, trainy, epochs=50, batch_size=128)end_time = datetime.datetime.now()running_time = end_time - start_time# 保存模型model.save('vmd_lstm_model.h5')# 返回構建好的模型return modelymodel = LSTM_model_train(trainx, trainy)5 預測
首先加載訓練好后的模型
# 加載模型
from tensorflow.keras.models import load_model
model = load_model('vmd_lstm_model.h5')準備好需要預測的數據，訓練時保留了6天的數據，將前5天的數據作為輸入預測，將預測的結果和最后一天的真實值進行比較。
y_true = datay[-timesteps-predict_steps:-timesteps]
x_pred = datax[-timesteps:]預測并計算誤差，并進行可視化，將這些步驟封裝為函數。
# 預測并計算誤差和可視化
def predict_and_plot(x, y_true, model, scaler, timesteps):# 變換輸入x格式，適應LSTM模型predict_x = np.reshape(x, (1, timesteps, feature_num)) # 預測predict_y = model.predict(predict_x)predict_y = scaler.inverse_transform(predict_y)y_predict = []y_predict.extend(predict_y[0])# 計算誤差r2 = r2_score(y_true, y_predict)rmse = mean_squared_error(y_true, y_predict, squared=False)mae = mean_absolute_error(y_true, y_predict)mape = mean_absolute_percentage_error(y_true, y_predict)print("r2: %.2f\nrmse: %.2f\nmae: %.2f\nmape: %.2f" % (r2, rmse, mae, mape))# 預測結果可視化cycol = cycle('bgrcmk')plt.figure(dpi=100, figsize=(14, 5))plt.plot(y_true, c=next(cycol), markevery=5)plt.plot(y_predict, c=next(cycol), markevery=5)plt.legend(['y_true', 'y_predict'])plt.xlabel('時間')plt.ylabel('功率(kW)')plt.show()return y_predicty_predict_nowork = predict_and_plot(x_pred, y_true, model, scaler2, timesteps)最后得到可視化結果，發下可視化結果并不是太好，可以通過調參和數據處理進一步提升模型預測效果。
