利用統計測試和機器學習分析和預測太陽能發電的性能測試和對比
本文將討論通過使用假設測試、特征工程、時間序列建模方法等從數據集中獲得有形價值的技術。我還將解決不同時間序列模型的數據泄漏和數據準備等問題,并且對常見的三種時間序列預測進行對比測試。
介紹
時間序列預測是一個經常被研究的話題,我們這里使用使用兩個太陽能電站的數據,研究其規律進行建模。首先將它們歸納為兩個問題來解決這些問題:
- 是否有可能識別出性能欠佳的太陽能組件?
- 是否可以預報兩天的太陽能發電量?
在繼續回答這些問題之前,讓我們先了解太陽能發電廠是如何發電的。

上圖描述了從太陽能電池板模塊到電網的發電過程。太陽能通過光電效應直接轉化為電能。當硅(太陽能電池板中最常見的半導體材料)等材料暴露在光線下時,光子(電磁能量的亞原子粒子)被吸收并釋放自由電子,從而產生直流電(DC)。使用逆變器,直流電被轉換成交流電(AC)并發送到電網,在那里它可以被分配到家庭。
數據
原始數據由每個太陽能發電廠的兩個逗號分隔值(CSV)文件組成。一份文件顯示了發電過程,另一份文件顯示了太陽能發電廠傳感器記錄的測量數據。每個太陽能發電廠的兩個數據集都被整理成一個pandas的df。
太陽能發電廠1號(SP1)和太陽能發電廠2號(SP2)的數據每15分鐘收集一次,從2020年5月15日到2020年6月18日。SP1和SP2數據集都包含相同的變量。
- Date Time -間隔15分鐘
- Ambient temperature-模塊周圍空氣的溫度
- Module temperature-模塊的溫度
- Irradiation——模塊上的輻射
- DC Power (kW) -直流
- AC Power (kW) -交流
- Daily Yield-每日發電量總和
- Total Yield-逆變器的累計產量
- Plant ID -太陽能電站的唯一標識
- Module ID -每個模塊的唯一標識
天氣傳感器用于記錄每個太陽能發電廠的環境溫度、組件溫度和輻射。
對于這個數據集直流功率將是因變量(目標變量)。我們目標是的試圖找到性能不佳的太陽能模塊。
兩個獨立的df用于分析和預測。唯一的區別是用于預測的數據被重新采樣為每小時的間隔,而用于分析的數據幀包含15分鐘的間隔。
首先我們刪除Plant ID,因為它對試圖回答上述問題沒有任何價值。Module ID也從預測數據集中刪除。表1和表2顯示了數據示例。


在繼續分析數據之前,我們對太陽能發電廠做了一些假設,包括:
- 數據采集儀器無故障
- 模塊定期清洗(忽略維護的影響)
- 兩個太陽能電站周圍都沒有遮擋問題
探索性數據分析(EDA)
對于數據科學的新手來說,EDA是通過繪圖可視化和執行統計測試來理解數據的關鍵一步。我們首先通過繪制SP1和SP2的DC和AC,可以觀察到每個太陽能發電廠的性能。

SP1顯示的直流功率比sp2高一個數量級。假設SP1采集的數據是正確的,用于記錄數據的儀器沒有故障,這就說明SP1中逆變器需要進行更深入的研究

通過按每個模塊的日頻率聚合AC和DC功率,圖3顯示了SP1中所有模塊的逆變器效率。根據領域內知識,太陽能逆變器的效率應該在93-96%之間。由于所有模塊的效率范圍為9.76% - 9.79%,這里可以說明需要調查逆變器的性能,以及是否需要更換。
由于SP1顯示了逆變器的問題,因此僅在SP2上進行了進一步的分析。
盡管這一小段分析是我們花了更多的時間對逆變器進行研究,但它并沒有回答確定太陽能模塊性能的主要問題。
由于SP2的逆變器正常工作,可以通過深入挖掘數據,來識別和調查任何異常情況。
圖4中顯示了模塊溫度和環境溫度之間的關系,并且有模塊溫度極高的情況。

這看起來似乎違反我們的認知,但是可以看到高溫對太陽能電池板的確有負面影響。當光子與太陽能電池內的電子接觸時,它們會釋放自由電子,但在更高的溫度下,更多的電子已經處于激發態,這降低了電池板可以產生的電壓,進而降低了效率。
考慮到這一現象,下面的圖5顯示了SP2的模塊溫度和直流功率(環境溫度低于模塊溫度的數據點和模塊運行數量較少的一天中的時間已經過過濾,以防止數據傾斜)。

在圖5中,紅線表示平均溫度。這里可以看到有一個明確的臨界點和直流電源停滯的跡象。在~52°C開始平穩。為了找到性能次優的太陽能模塊,所有顯示模塊溫度超過52°C的行都被刪除。
下面的圖6顯示了SP2中每個模塊在一天中的直流功率。這樣就基本符合了預期,午間發電量較大。但是還有個問題,在運行高峰時期,發電量較低。我們很難總結造成這種情況的原因,因為當天的天氣條件可能很差,或者SP2可能需要進行日常的維護等等。
圖6中也有低性能模塊的跡象。它們可以被識別為圖上偏離最近群集的模塊(單個數據點)。

為了確定哪些模塊表現不佳,我們可以進行統計測試,同時將每個模塊的性能與其他模塊進行比較,從而確定性能。
每隔15分鐘,不同模塊的直流電源在同一時間的分布是正態分布,通過假設檢驗可以確定哪些模塊表現不佳。計數是指模塊落在99.9%置信區間之外且p值< 0.001的次數。
圖7按降序顯示了每個模塊在統計上顯著低于同期其他模塊的次數。

從圖7中可以清楚地看出,模塊' Quc1TzYxW2pYoWX '是有問題的。這些信息可以提供給SP2的相關工作人員,調查原因。
建模
下面我們開始使用三種不同的時間序列算法:SARIMA、XGBoost和CNN-LSTM,進行建模并比較
對于所有三個模型,都使用預測下一個數據點進行預測。Walk-forward驗證是一種用于時間序列建模的技術,因為隨著時間的推移,預測會變得不那么準確,因此更實用的方法是在實際數據可用時,用實際數據重新訓練模型。
在建模之前需要更詳細地研究數據。圖8顯示了SP2數據集中所有特征的相關熱圖。熱圖顯示了因變量直流功率,與模塊溫度、輻照和環境溫度的強相關性。這些特征可能在預測中發揮重要作用。
在下面的熱圖中,交流功率顯示皮爾森相關系數為1。為了防止數據泄漏問題,我們將直流功率從數據中刪除。

SARIMA
季節自回歸綜合移動平均(SARIMA)是一種單變量時間序列預測方法。由于目標變量顯示出24小時循環周期的跡象,SARIMA是一個有效的建模選項,因為它考慮了季節影響。這可以從下面的季節分解圖中觀察到。

SARIMA算法要求數據是平穩的。有多種方法來檢驗數據是否平穩,例如統計檢驗(增強迪基-福勒檢驗),匯總統計(比較數據的不同部分的均值/方差)和可視化分析數據。在建模之前進行多次測試是很重要的。
增強迪基-富勒(ADF)檢驗是一種“單位根檢驗”,用于確定時間序列是否平穩。從根本上說,這是一個統計顯著性檢驗,其中存在一個零假設和替代假設,并根據得出的p值得出結論。
零假設:時間序列數據是非平穩的。
替代假設:時間序列數據是平穩的。
在我們的例子中,如果p值≤0.05,我們可以拒絕原假設,并確認數據沒有單位根。
from statsmodels.tsa.stattools import adfuller
result = adfuller(plant2_dcpower.values)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
print('\t%s: %.3f' % (key, value))

從ADF檢驗來看,p值為0.000553,< 0.05。根據這一統計數據,可以認為該數據是穩定的。然而,查看圖2(最上面的圖),有明顯的季節性跡象(對于被認為是平穩的時間序列數據,不應該有季節性和趨勢的跡象),這說明數據是非平穩的。因此,運行多個測試非常重要。
為了用SARIMA對因變量建模,時間序列需要是平穩的。如圖9(第一個和第三個圖)所示,直流電有明顯的季節性跡象。取第一個差值[t-(t-1)]去除季節性成分,如圖10所示,因為它看起來類似于正態分布。數據現在是平穩的,適用于SARIMA算法。

SARIMA的超參數包括p(自回歸階數)、d(差階數)、q(移動平均階數)、p(季節自回歸階數)、d(季節差階數)、q(季節移動平均階數)、m(季節周期的時間步長)、trend(確定性趨勢)。

圖11顯示了自相關(ACF)、部分自相關(PACF)和季節性ACF/PACF圖。ACF圖顯示了時間序列與其延遲版本之間的相關性。PACF顯示了時間序列與其滯后版本之間的直接相關性。藍色陰影區域表示置信區間。SACF和SPACF可以通過從原始數據中取季節差(m)來計算,在本例中為24,因為在ACF圖中有一個明顯的24小時的季節效應。
根據我們的直覺,超參數的起點可以從ACF和PACF圖中推導出來。如ACF和PACF均呈逐漸下降的趨勢,即自回歸階數(p)和移動平均階數(q)均大于0。p和p可以通過分別觀察PCF和SPCF圖,并計算滯后值不顯著之前具有統計學顯著性的滯后數來確定。同樣,q和q可以在ACF和SACF圖中找到。
差階(d)可以通過使數據平穩的差的數量來確定。季節差異階數(D)是根據從時間序列中去除季節性成分所需的差異數來估計的。
這些超參數選擇可以看這篇文章:https://arauto.readthedocs.io/en/latest/how_to_choose_terms.html
也可以采用網格搜索方法進行超參數優化,根據最小均方誤差(MSE)選擇最優超參數,包括p = 2, d = 0, q = 4, p = 2, d = 1, q = 6, m = 24, trend = ' n '(無趨勢)。
from time import time
from sklearn.metrics import mean_squared_error
from statsmodels.tsa.statespace.sarimax import SARIMAX
configg = [(2, 1, 4), (2, 1, 6, 24), 'n']
def train_test_split(data, test_len=48):
"""
Split data into training and testing.
"""
train, test = data[:-test_len], data[-test_len:]
return train, test
def sarima_model(data, cfg, test_len, i):
"""
SARIMA model which outputs prediction and model.
"""
order, s_order, t = cfg[0], cfg[1], cfg[2]
model = SARIMAX(data, order=order, seasonal_order=s_order, trend=t,
enforce_stationarity=False, enfore_invertibility=False)
model_fit = model.fit(disp=False)
yhat = model_fit.predict(len(data))
if i + 1 == test_len:
return yhat, model_fit
else:
return yhat
def walk_forward_val(data, cfg):
"""
A walk forward validation technique used for time series data. Takes current value of x_test and predicts
value. x_test is then fed back into history for the next prediction.
"""
train, test = train_test_split(data)
pred = []
history = [i for i in train]
test_len = len(test)
for i in range(test_len):
if i + 1 == test_len:
yhat, s_model = sarima_model(history, cfg, test_len, i)
pred.append(yhat)
mse = mean_squared_error(test, pred)
return pred, mse, s_model
else:
yhat = sarima_model(history, cfg, test_len, i)
pred.append(yhat)
history.append(test[i])
pass
if __name__ == '__main__':
start_time = time()
sarima_pred_plant2, sarima_mse, s_model = walk_forward_val(plant2_dcpower, configg)
time_len = time() - start_time
print(f'SARIMA runtime: {round(time_len/60,2)} mins')

圖12顯示了SARIMA模型的預測值與SP2 2天內記錄的直流功率的比較。
為了分析模型的性能,圖13顯示了模型診斷。相關圖顯示在第一個滯后后幾乎沒有相關性,下面的直方圖顯示在平均值為零附近的正態分布。由此我們可以說模型無法從數據中收集到進一步的信息。

XGBoost
XGBoost (eXtreme Gradient Boosting)是一種梯度增強決策樹算法。它使用集成方法,其中添加新的決策樹模型來修改現有的決策樹分數。與SARIMA不同的是,XGBoost是一種多元機器學習算法,這意味著該模型可以采用多特征來提高模型性能。
我們采用特征工程提高模型精度。還創建了3個附加特性,其中包括AC和DC功率的滯后版本,分別為S1_AC_POWER和S1_DC_POWER,以及通過交流功率除以直流功率的總體效率EFF。并將AC_POWER和MODULE_TEMPERATURE從數據中刪除。圖14通過增益(使用一個特征的分割的平均增益)和權重(一個特征在樹中出現的次數)顯示了特征的重要性級別。

通過網格搜索確定建模使用的超參數,結果為:*learning rate = 0.01, number of estimators = 1200, subsample = 0.8, colsample by tree = 1, colsample by level = 1, min child weight = 20 and max depth = 10
我們使用MinMaxScaler將訓練數據縮放到0到1之間(也可以試驗其他縮放器,如log-transform和standard-scaler,這取決于數據的分布)。通過將所有自變量向后移動一段時間,將數據轉換為監督學習數據集。
import numpy as np
import pandas as pd
import xgboost as xgb
from sklearn.preprocessing import MinMaxScaler
from time import time
def train_test_split(df, test_len=48):
"""
split data into training and testing.
"""
train, test = df[:-test_len], df[-test_len:]
return train, test
def data_to_supervised(df, shift_by=1, target_var='DC_POWER'):
"""
Convert data into a supervised learning problem.
"""
target = df[target_var][shift_by:].values
dep = df.drop(target_var, axis=1).shift(-shift_by).dropna().values
data = np.column_stack((dep, target))
return data
def xgb_forecast(train, x_test):
"""
XGBOOST model which outputs prediction and model.
"""
x_train, y_train = train[:,:-1], train[:,-1]
xgb_model = xgb.XGBRegressor(learning_rate=0.01, n_estimators=1500, subsample=0.8,
colsample_bytree=1, colsample_bylevel=1,
min_child_weight=20, max_depth=14, objective='reg:squarederror')
xgb_model.fit(x_train, y_train)
yhat = xgb_model.predict([x_test])
return yhat[0], xgb_model
def walk_forward_validation(df):
"""
A walk forward validation approach by scaling the data and changing into a supervised learning problem.
"""
preds = []
train, test = train_test_split(df)
scaler = MinMaxScaler(feature_range=(0,1))
train_scaled = scaler.fit_transform(train)
test_scaled = scaler.transform(test)
train_scaled_df = pd.DataFrame(train_scaled, columns = train.columns, index=train.index)
test_scaled_df = pd.DataFrame(test_scaled, columns = test.columns, index=test.index)
train_scaled_sup, test_scaled_sup = data_to_supervised(train_scaled_df), data_to_supervised(test_scaled_df)
history = np.array([x for x in train_scaled_sup])
for i in range(len(test_scaled_sup)):
test_x, test_y = test_scaled_sup[i][:-1], test_scaled_sup[i][-1]
yhat, xgb_model = xgb_forecast(history, test_x)
preds.append(yhat)
np.append(history,[test_scaled_sup[i]], axis=0)
pred_array = test_scaled_df.drop("DC_POWER", axis=1).to_numpy()
pred_num = np.array([pred])
pred_array = np.concatenate((pred_array, pred_num.T), axis=1)
result = scaler.inverse_transform(pred_array)
return result, test, xgb_model
if __name__ == '__main__':
start_time = time()
xgb_pred, actual, xgb_model = walk_forward_validation(dropped_df_cat)
time_len = time() - start_time
print(f'XGBOOST runtime: {round(time_len/60,2)} mins')
圖15顯示了XGBoost模型的預測值與SP2 2天內記錄的直流功率的比較。

CNN-LSTM
CNN-LSTM (convolutional Neural Network Long - Short-Term Memory)是兩種神經網絡模型的混合模型。CNN是一種前饋神經網絡,在圖像處理和自然語言處理方面表現出了良好的性能。它還可以有效地應用于時間序列數據的預測。LSTM是一種序列到序列的神經網絡模型,旨在解決長期存在的梯度爆炸/消失問題,使用內部存儲系統,允許它在輸入序列上積累狀態。
在本例中,使用CNN-LSTM作為編碼器-解碼器體系結構。由于CNN不直接支持序列輸入,所以我們通過1D CNN讀取序列輸入并自動學習重要特征。然后LSTM進行解碼。與XGBoost模型類似,使用scikitlearn的MinMaxScaler使用相同的數據并進行縮放,但范圍在-1到1之間。對于CNN-LSTM,需要將數據重新整理為所需的結構:[samples, subsequences, timesteps, features],以便可以將其作為輸入傳遞給模型。
由于我們希望為每個子序列重用相同的CNN模型,因此使用timedidistributedwrapper對每個輸入子序列應用一次整個模型。在下面的圖16中可以看到最終模型中使用的不同層的模型摘要。

在將數據分解為訓練數據和測試數據之后,將訓練數據分解為訓練數據和驗證數據集。在所有訓練數據(包括驗證數據)的每次迭代之后,模型可以進一步使用這一點來評估模型的性能。
學習曲線是深度學習中使用的一個很好的診斷工具,它顯示了模型在每個階段之后的表現。下面的圖17顯示了模型如何從數據中學習,并顯示了驗證數據與訓練數據的收斂。這是良好模特訓練的標志。

import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import MinMaxScaler
import keras
from keras.models import Sequential
from keras.layers.convolutional import Conv1D, MaxPooling1D
from keras.layers import LSTM, TimeDistributed, RepeatVector, Dense, Flatten
from keras.optimizers import Adam
n_steps = 1
subseq = 1
def train_test_split(df, test_len=48):
"""
Split data in training and testing. Use 48 hours as testing.
"""
train, test = df[:-test_len], df[-test_len:]
return train, test
def split_data(sequences, n_steps):
"""
Preprocess data returning two arrays.
"""
x, y = [], []
for i in range(len(sequences)):
end_x = i + n_steps
if end_x > len(sequences):
break
x.append(sequences[i:end_x, :-1])
y.append(sequences[end_x-1, -1])
return np.array(x), np.array(y)
def CNN_LSTM(x, y, x_val, y_val):
"""
CNN-LSTM model.
"""
model = Sequential()
model.add(TimeDistributed(Conv1D(filters=14, kernel_size=1, activation="sigmoid",
input_shape=(None, x.shape[2], x.shape[3]))))
model.add(TimeDistributed(MaxPooling1D(pool_size=1)))
model.add(TimeDistributed(Flatten()))
model.add(LSTM(21, activation="tanh", return_sequences=True))
model.add(LSTM(14, activation="tanh", return_sequences=True))
model.add(LSTM(7, activation="tanh"))
model.add(Dense(3, activation="sigmoid"))
model.add(Dense(1))
model.compile(optimizer=Adam(learning_rate=0.001), loss="mse", metrics=['mse'])
history = model.fit(x, y, epochs=250, batch_size=36,
verbose=0, validation_data=(x_val, y_val))
return model, history
# split and resahpe data
train, test = train_test_split(dropped_df_cat)
train_x = train.drop(columns="DC_POWER", axis=1).to_numpy()
train_y = train["DC_POWER"].to_numpy().reshape(len(train), 1)
test_x = test.drop(columns="DC_POWER", axis=1).to_numpy()
test_y = test["DC_POWER"].to_numpy().reshape(len(test), 1)
#scale data
scaler_x = MinMaxScaler(feature_range=(-1,1))
scaler_y = MinMaxScaler(feature_range=(-1,1))
train_x = scaler_x.fit_transform(train_x)
train_y = scaler_y.fit_transform(train_y)
test_x = scaler_x.transform(test_x)
test_y = scaler_y.transform(test_y)
# shape data into CNN-LSTM format [samples, subsequences, timesteps, features] ORIGINAL
train_data_np = np.hstack((train_x, train_y))
x, y = split_data(train_data_np, n_steps)
x_subseq = x.reshape(x.shape[0], subseq, x.shape[1], x.shape[2])
# create validation set
x_val, y_val = x_subseq[-24:], y[-24:]
x_train, y_train = x_subseq[:-24], y[:-24]
n_features = x.shape[2]
actual = scaler_y.inverse_transform(test_y)
# run CNN-LSTM model
if __name__ == '__main__':
start_time = time()
model, history = CNN_LSTM(x_train, y_train, x_val, y_val)
prediction = []
for i in range(len(test_x)):
test_input = test_x[i].reshape(1, subseq, n_steps, n_features)
yhat = model.predict(test_input, verbose=0)
yhat_IT = scaler_y.inverse_transform(yhat)
prediction.append(yhat_IT[0][0])
time_len = time() - start_time
mse = mean_squared_error(actual.flatten(), prediction)
print(f'CNN-LSTM runtime: {round(time_len/60,2)} mins')
print(f"CNN-LSTM MSE: {round(mse,2)}")
圖18顯示了CNN-LSTM模型的預測值與SP2 2天內記錄的直流功率的對比。

由于CNN-LSTM的隨機性,該模型運行10次,并記錄一個平均MSE值作為最終值,以判斷模型的性能。圖19顯示了為所有模型運行記錄的mse的范圍。

結果對比
下表顯示了每個模型的MSE (CNN-LSTM的平均MSE)和每個模型的運行時間(以分鐘為單位)。

從表中可以看出,XGBoost的MSE最低、運行時第二快,并且與所有其他模型相比具有最佳性能。由于該模型顯示了一個可以接受的每小時預測的運行時,它可以成為幫助運營經理決策過程的強大工具。
總結
在本文中我們分析了SP1和SP2,確定SP1性能較低。所以對SP2的進一步調查顯示,并且查看了SP2中那些模塊性能可能有問題,并使用假設檢驗來計算每個模塊在統計上明顯表現不佳的次數,' Quc1TzYxW2pYoWX '模塊顯示了約850次低性能計數。
我們使用數據訓練三個模型:SARIMA、XGBoost和CNN-LSTM。SARIMA表現最差,XGBOOST表現最好,MSE為16.9,運行時間為1.43 min。所以可以說XGBoost在表格數據中還是最優先得選擇。