余 楊，彭湘媛

（同濟(jì)大學(xué)控制科學(xué)與工程系，上海 200092）

液化天然氣（LNG）為當(dāng)今世界發(fā)展最快的燃料之一。據(jù)預(yù)測，世界上的LNG 需求到2030年會(huì)比之前高出兩倍［1-2］。近十幾年來，我國LNG 消費(fèi)的需求量飛速上升［3］，尤其是在東部地區(qū)和東南沿海等經(jīng)濟(jì)較為發(fā)達(dá)的區(qū)域，LNG在能源消費(fèi)總值中的占比不斷提升。根據(jù)中國可持續(xù)發(fā)展油氣資源戰(zhàn)略研究報(bào)告的預(yù)測，未來15年，我國LNG 需求量將出現(xiàn)指數(shù)增長，年均增速達(dá)10.8%，而LNG 生產(chǎn)量增長的速度卻遠(yuǎn)低于需求增長速度，年增長率僅為7.5%，依賴進(jìn)口的程度會(huì)越來越大［4］。

同時(shí)，LNG能源具有一定的不確定性，這也使整個(gè)產(chǎn)業(yè)面臨著一定的挑戰(zhàn)。要在市場中站穩(wěn)腳跟，實(shí)現(xiàn)利益最大化，需要提前嚴(yán)格把控每月的LNG 儲(chǔ)量，在不造成積壓的前提下盡量滿足市場需求。LNG銷售公司需要綜合考慮各種影響因素，對公司LNG 儲(chǔ)量做出相應(yīng)預(yù)測，提高經(jīng)濟(jì)效益，避免LNG 庫存儲(chǔ)量不夠、滿足不了市場需求，或是儲(chǔ)量過剩但銷售渠道不足［5］。

針對LNG 儲(chǔ)量預(yù)測問題，本文首先闡述了針對實(shí)際天然氣儲(chǔ)量及影響因素的一些預(yù)處理方法，對影響因素的相關(guān)性進(jìn)行了分析，運(yùn)用互信息理論對影響因素進(jìn)行篩選。然后利用LSTM 長短期記憶（Long Short Term Memory）神經(jīng)網(wǎng)絡(luò)［6］對LNG儲(chǔ)量進(jìn)行了預(yù)測并對模型進(jìn)行優(yōu)化，檢驗(yàn)了預(yù)測結(jié)果的準(zhǔn)確度和變化趨勢，最后與傳統(tǒng)預(yù)測方法——灰度預(yù)測［7］對比，驗(yàn)證LSTM模型的正確性和可行性。

1 基于互信息的相關(guān)因素分析

根據(jù)實(shí)際調(diào)研，影響LNG 儲(chǔ)量的因素有LNG 儲(chǔ)罐大小、運(yùn)營水平、市場價(jià)格、美元兌人民幣匯率、通航能力、碼頭能力和天氣等。這些因素對LNG 儲(chǔ)量的影響不盡相同。本文基于互信息的方法對這些因素進(jìn)行相關(guān)性分析［8-10］。

1.1 基于互信息的相關(guān)性分析

信息熵（Entropy）是信息的基本單位，用于描述離散變量是分散還是聚合。信息熵越大，代表變量的分布越離散，描述該變量需要的信息越多。離散型隨機(jī)變量X的信息熵定義如式（1）所示，即

其中，p（x）表示每種可能的取值x的概率，底數(shù)b可取不同的值，表示信息熵有不同的量綱。

互信息表示兩個(gè)或多個(gè)變量之間共享的信息量?；バ畔⒃酱?，變量之間的相關(guān)性越高。本文均為離散變量，對于兩個(gè)隨機(jī)離散變量X，Y之間的互信息定義如式（2）所示，即

其中，Xi代表變量X中的第i個(gè)值；Yj代表變量Y中的第j個(gè)值；P(Xi，Yj)表示變量X為第i個(gè)值與變量Y為第j個(gè)值時(shí)的聯(lián)合概率密度，P(X) 和P(Y) 為獨(dú)立密度。I(X；Y)越大，變量X包含關(guān)于Y的信息就越多。因此可用互信息來度量變量間的相關(guān)性。

通過式（2）計(jì)算LNG 儲(chǔ)量各影響因素與因變量的互信息，將互信息超過一定閾值（T1）的變量篩選出來，作為冗余過濾器的輸入集合。經(jīng)過計(jì)算，各變量與LNG儲(chǔ)量的相關(guān)度見表1。

表1 LNG影響因素相關(guān)度排名Tab.1 The rank of the correlation degree of LNG influencing facton

計(jì)算出各影響因素的相關(guān)性后，經(jīng)過反復(fù)實(shí)驗(yàn)，將相關(guān)性過濾器的閾值T1定位為0.050000，即相關(guān)性小于0.050000的影響因素被濾除，不再進(jìn)入冗余過濾器。因此，氣溫（℃）這一影響因素因?yàn)橄嚓P(guān)度僅為0.002439而被濾除。

1.2 基于互信息的冗余過濾器

一個(gè)特征值關(guān)于目標(biāo)變量的信息內(nèi)容可以被其他特征值表示出來時(shí)，這個(gè)特征值稱為冗余特征值。冗余特征不僅會(huì)惡化預(yù)測性能，而且還會(huì)影響算法學(xué)習(xí)的速度。因此，用于預(yù)測數(shù)據(jù)必須去除冗余特征值。本文提出了一種基于互信息的冗余性分析方法。算法的具體步驟如下：

（1）所有候選特征值進(jìn)行線性歸一化處理，范圍為［0，1］。歸一化的具體表達(dá)式如式（3）所示，即

其中，Xmax，Xmin分別代表變量X中的最大值和最小值，Y為歸一化后的結(jié)果。

（2）計(jì)算待選輸入與目標(biāo)變量相關(guān)性。第i個(gè)待選輸入Xi對于目標(biāo)變量Y的相關(guān)性記為D(Xi)，計(jì)算公式如式（2）。

（3）篩選出相關(guān)度超過預(yù)先設(shè)定的閾值T1的m個(gè)特征值，并放入相關(guān)性選擇特征集SRelevance中，用作冗余性分析。

（4）將SRelevance中相關(guān)性排名第一的特征值添加到SFinal中，將冗余過濾計(jì)數(shù)器的i變?yōu)?。

（5）SRelevance中被選中的第i個(gè)特征值的冗余度記為R(i)，計(jì)算公式如（4）所示，即

其中，Xs代表已經(jīng)通過冗余過濾器而得到的特征量，R（i）表示Xi與SFinal中特征量的最大冗余值。

（6）如果R（i）大于提前設(shè)定的閾值T2，特征值Xi冗余，將其濾除。否則Xi添加到SFinal，冗余過濾器的計(jì)數(shù)器加1。

（7）如果i達(dá)到m+1，結(jié)束特征值選擇。SFinal即為選擇后的結(jié)果，準(zhǔn)備輸入給LSTM 神經(jīng)網(wǎng)絡(luò)進(jìn)行預(yù)測訓(xùn)練。

經(jīng)過冗余濾波器篩選，可以看到這五個(gè)變量中通航能力的冗余性較大，即含有較多的重復(fù)信息。另一方面，因包含的信息越多，神經(jīng)網(wǎng)絡(luò)訓(xùn)練所得的結(jié)果就越準(zhǔn)確，因此也不能將變量減少過多。經(jīng)過反復(fù)實(shí)驗(yàn)，將冗余過濾器的閾值T2設(shè)定為0.80000。冗余度大于T2的變量，即通航能力將被濾除，無需進(jìn)入后續(xù)LSTM神經(jīng)網(wǎng)絡(luò)預(yù)測。

2 LSTM 神經(jīng)網(wǎng)絡(luò)預(yù)測模型

為了解決時(shí)序數(shù)列預(yù)測的需求問題，循環(huán)神經(jīng)網(wǎng)絡(luò)RNN 被提出，其獨(dú)特的網(wǎng)絡(luò)結(jié)構(gòu)為處理上下文存在關(guān)聯(lián)性的內(nèi)容提供了極大的方便，解決了BP神經(jīng)網(wǎng)絡(luò)等傳統(tǒng)神經(jīng)網(wǎng)絡(luò)存在的效率不高、準(zhǔn)確性較低等問題。傳統(tǒng)的RNN 神經(jīng)網(wǎng)絡(luò)會(huì)在層與層之間傳遞信息時(shí)出現(xiàn)爆炸或梯度消失的情況，并且每次都會(huì)采用獨(dú)立的數(shù)據(jù)向量，沒有內(nèi)存幫助其處理需要記憶的任務(wù)，因此LSTM長短期記憶預(yù)測模型被提出。

依據(jù)相關(guān)性信息的篩選，本文使用了5 類輸入特征，包括運(yùn)營水平（車）、市場價(jià)格（元/噸）、美元兌人民幣匯率、碼頭能力（船）以及LNG 歷史儲(chǔ)量（萬噸），通過不同輸入數(shù)據(jù)篩選出對LNG 儲(chǔ)量預(yù)測較為合適的特征進(jìn)行分析。

2.1 輸入數(shù)據(jù)處理

（1）異常值處理。

在構(gòu)造數(shù)據(jù)集過程中難免會(huì)遇到異常值問題，由于相關(guān)性特征構(gòu)造時(shí)需要對所有的數(shù)據(jù)進(jìn)行簡化操作，最終主要的時(shí)間序列以篩選后的數(shù)據(jù)為準(zhǔn)，其他相關(guān)的多余數(shù)據(jù)只能做刪除處理；同時(shí)當(dāng)碼頭因突發(fā)情況無法運(yùn)行時(shí)，需要對其進(jìn)行異常值填充處理，對于因政治原因無法通行的情況直接將其填充為0。

（2）數(shù)據(jù)歸一化處理。

由于影響因素間各個(gè)數(shù)據(jù)存在量級(jí)上的差距，有的特征如市場價(jià)格達(dá)到了千以上的數(shù)量級(jí)，而有的特征如碼頭能力、通航能力只有1-10 數(shù)量級(jí)，在神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)時(shí)會(huì)偏向于將數(shù)值較大的特征認(rèn)為是主要特征而忽略了細(xì)小的特征，然而對LNG 銷售產(chǎn)業(yè)而言，碼頭能力的關(guān)注度卻經(jīng)常大于市場價(jià)格。為了不對預(yù)測結(jié)果造成偏差，需要對各個(gè)特征進(jìn)行歸一化處理，將數(shù)值縮放到［0，1］區(qū)間，同時(shí)也是為了便于神經(jīng)網(wǎng)絡(luò)進(jìn)行學(xué)習(xí)，歸一化方程如式（5）所示，即

2.2 模型結(jié)構(gòu)設(shè)計(jì)

LSTM 網(wǎng)絡(luò)模型結(jié)構(gòu)主要是“輸入層-隱藏層-輸出層”，模型的隱藏層搭建通常運(yùn)用Keras 學(xué)習(xí)框架中的Dense 層、PReLU 層、LSTM 層和Dropout 層。在輸出層方面，LNG 儲(chǔ)量回歸預(yù)測利用Dense 層對預(yù)測結(jié)果進(jìn)行回歸，輸出值即為預(yù)測值。輸入層的神經(jīng)元個(gè)數(shù)由輸入?yún)?shù)所決定，輸出層結(jié)構(gòu)根據(jù)預(yù)測類型不同進(jìn)行區(qū)分，因?yàn)殡[藏層結(jié)構(gòu)目前還沒有具體的理論指導(dǎo)，因此嘗試使用Dense層與LSTM 層結(jié)合進(jìn)行多種模式探索。本文所用模型的網(wǎng)絡(luò)結(jié)構(gòu)如圖1 所示，其中隱藏層的層數(shù)需要經(jīng)過測試才能最終得出。

圖1 LSTM 模型結(jié)構(gòu)Fig.1 The structure of LSTM

網(wǎng)絡(luò)模型可以分為以Dense 為主的結(jié)構(gòu)、以LSTM為主的結(jié)構(gòu)以及兩者混合結(jié)構(gòu)。隱含層中根據(jù)需要添加Dropout層防止過擬合。加入Dropout就是為了隨機(jī)刪減某些數(shù)據(jù)，讓模型有更多種嘗試的可能性，從而提高模型的泛化能力。其中選擇PReLU 作為激活函數(shù)。神經(jīng)元個(gè)數(shù)關(guān)乎到模型學(xué)習(xí)速度的快慢，但是過多和過少都對訓(xùn)練沒有幫助，因此將使用試湊法進(jìn)行測試。

2.3 模型初步預(yù)測結(jié)果

模型初始時(shí)使用常用損失函數(shù)MSE（Mean Squared Error，均方誤差），優(yōu)化器為SGD，批量數(shù)初始值設(shè)置為10，訓(xùn)練次數(shù)為2000 次。仿真結(jié)果如圖2 所示。圖2（a）表明，預(yù)測結(jié)果與實(shí)際值仍有一定的誤差，尤其是在起伏較大的位置，預(yù)測值與實(shí)際值偏離較為嚴(yán)重。此外，預(yù)測趨勢應(yīng)該保持基本正確，即后一個(gè)月預(yù)測值較前一個(gè)月的儲(chǔ)量是增加還是減少，應(yīng)與實(shí)際情況保持一致。這樣可以增加企業(yè)運(yùn)營的穩(wěn)定性，使其更有能力對未來市場進(jìn)行把控。如果趨勢出現(xiàn)了較大偏差，同樣會(huì)引起較大利益損失。因此，本文對趨勢的變化進(jìn)行了作圖檢驗(yàn)（用后一個(gè)月儲(chǔ)量數(shù)值減去前一個(gè)月儲(chǔ)量數(shù)值），結(jié)果如圖2（b）所示。由圖2（b）可以看出，預(yù)測趨勢與實(shí)際變化趨勢存在嚴(yán)重的不一致。

圖2 初始LSTM模型（沒有優(yōu)化）的預(yù)測結(jié)果Fig.2 The prediction of LNG reserves by initial LSTM model

3 LSTM預(yù)測模型的優(yōu)化

如前文所述，LSTM 模型對于預(yù)測結(jié)果并不達(dá)標(biāo)，因此本文又在多個(gè)方面對其進(jìn)行了優(yōu)化，使得最終的相對誤差絕大多數(shù)控制在5% 以內(nèi)，即使存在個(gè)別異常值，也不超過10%。

3.1 修改Batch Size

Batch Size 為一次訓(xùn)練所選取的樣本數(shù)。在沒有使用Batch Size 之前，意味著網(wǎng)絡(luò)在訓(xùn)練時(shí)，是一次性地把所有的數(shù)據(jù)（整個(gè)數(shù)據(jù)庫）輸入網(wǎng)絡(luò)中，然后計(jì)算其梯度進(jìn)行反向傳播，由于在計(jì)算梯度時(shí)使用了整個(gè)數(shù)據(jù)庫，所以計(jì)算得到的梯度方向更為準(zhǔn)確。但在這種情況下，計(jì)算得到不同梯度值差別巨大，難以使用一個(gè)全局的學(xué)習(xí)率。在樣本數(shù)很小的數(shù)據(jù)庫中，沒有Batch Size 是可以接受的，而且效果也很好。但是一旦是稍大些的數(shù)據(jù)庫，一次性地把所有數(shù)據(jù)輸入網(wǎng)絡(luò)，肯定會(huì)引起內(nèi)存的爆炸，所以這時(shí)就需要一個(gè)Batch Size。

Batch Size 的數(shù)值關(guān)乎模型的優(yōu)化效率和速度，同時(shí)直接影響到GPU 內(nèi)存的使用情況，假如GPU 內(nèi)存不大或是數(shù)據(jù)集較小，該數(shù)值最好相應(yīng)設(shè)置得小一些。設(shè)置合適的Batch Size會(huì)使梯度準(zhǔn)確，預(yù)測效率和準(zhǔn)確率同時(shí)也得到提高。Batch Size 的優(yōu)化結(jié)果如表2所示。

由表2 可見，改變訓(xùn)練樣本數(shù)，預(yù)測效果會(huì)發(fā)生較為明顯的改變，從原定的Batch Size 為10 開始逐漸減小，預(yù)測曲線更加貼近真實(shí)值，相對誤差減小，準(zhǔn)確度提升。當(dāng)Batch Size 為4 時(shí)，模型達(dá)到飽和狀態(tài)，之后繼續(xù)減小訓(xùn)練樣本數(shù)，誤差將會(huì)再次增大。因此將Batch Size的最終值定為4。

表2 Batch size 優(yōu)化預(yù)測結(jié)果統(tǒng)計(jì)表Tab.2 The statistics of Batch size optimization effect

3.2 修改損失函數(shù)

通過損失函數(shù)可以看到模型的優(yōu)劣，為研究提供了優(yōu)化的方向，但是每個(gè)模型的損失函數(shù)并不是通用的［11］。損失函數(shù)的選取依賴于參數(shù)的數(shù)量、異常值、機(jī)器學(xué)習(xí)算法、梯度下降的效率、導(dǎo)數(shù)求取的難易和預(yù)測的置信度等若干方面。原損失函數(shù)為均方誤差MSE，它是預(yù)測值與目標(biāo)值之間差值的平方和，公式如式（6）所示，即

這一步優(yōu)化中將損失函數(shù)由MSE變換為平均絕對誤差（Mean Absolute Error）。MAE是另一種常用的回歸損失函數(shù)，它是目標(biāo)值與預(yù)測值之差絕對值的和，表示預(yù)測值的平均誤差幅度，而不需要考慮誤差的方向，范圍是0到無窮，計(jì)算公式如式（7）所示，即

一般情況下，利用均方誤差更容易求解，但平均絕對誤差則對于異常值更穩(wěn)健。由于均方誤差（MSE）在浮動(dòng)較大的位置損失通常高于平均絕對誤差（MAE），它會(huì)給這樣的值賦予更大的權(quán)重，而模型會(huì)試圖減小異常值造成的誤差，導(dǎo)致模型整體表現(xiàn)下降。所以當(dāng)訓(xùn)練數(shù)據(jù)中含有較多的異常值或起伏較大時(shí)，平均絕對誤差（MAE）更為適合。當(dāng)對所有觀測值進(jìn)行處理時(shí)，如果利用MSE進(jìn)行優(yōu)化會(huì)得到所有觀測的均值，而使用MAE則能得到所有觀測的中值。與均值相比，中值對于異常值的魯棒性更好，這就意味著平均絕對誤差對于異常值有著比均方誤差更好的魯棒性。MSE和MAE的預(yù)測結(jié)果如表3所示。

表3 損失函數(shù)優(yōu)化結(jié)果統(tǒng)計(jì)表Tab.3 The statistics of loss function optimization effect

由表3 可以看出，將損失函數(shù)由MSE 調(diào)整為MAE后，由于原始數(shù)據(jù)本身起伏較大，更適用于平均絕對誤差，預(yù)測效果得到較大程度改善。但與目標(biāo)值5% 的誤差仍然有一定距離，因此還需要繼續(xù)優(yōu)化。

3.3 修改優(yōu)化器

其實(shí)機(jī)器學(xué)習(xí)訓(xùn)練過程的本質(zhì)就是使損失最小化，故定義損失函數(shù)后，優(yōu)化器會(huì)在其中發(fā)揮重要作用。在深度學(xué)習(xí)中，通常對于梯度進(jìn)行優(yōu)化，優(yōu)化的目標(biāo)就是優(yōu)化網(wǎng)絡(luò)模型里的參數(shù)θ（θ是一個(gè)集合，θ1、θ2、θ3，…）。本文將原優(yōu)化器SGD 修改為RMSProp。理由是SGD 在隨機(jī)選擇梯度的同時(shí)會(huì)引入噪聲，使得權(quán)值更新的方向不一定正確，即更新比較頻繁，會(huì)造成損失函數(shù)有嚴(yán)重的震蕩。此外，SGD 也不能單獨(dú)克服局部最優(yōu)解的問題。

相比之下，RMSProp算法給每一個(gè)權(quán)值一個(gè)變量，MeanSquare（w，t）用來記錄第t次更新步長時(shí)前t-1 次的梯度平方的平均值。然后再用第t次的梯度除以前t-1 次的梯度平方的平均值，得到學(xué)習(xí)步長的更新比例。迭代更新公式如式（8-9）［13］，即

通過RMSprop，可以調(diào)整不同維度上的步長，加快收斂速度。根據(jù)比例會(huì)得到新的學(xué)習(xí)步長。如果當(dāng)前得到的梯度為負(fù)，學(xué)習(xí)步長就會(huì)減小一點(diǎn)；如果當(dāng)前得到的梯度為正，學(xué)習(xí)步長就會(huì)增大一點(diǎn)。 因此RMSProp算法步長的更新更加緩和。優(yōu)化器修改前后的預(yù)測結(jié)果如表4所示。

表4 優(yōu)化器優(yōu)化結(jié)果統(tǒng)計(jì)表Tab.4 The statistics of optimizer optimization

由表4 可知，將優(yōu)化器改為RMSProp 后，預(yù)測結(jié)果達(dá)到了小于5%的目標(biāo)，實(shí)現(xiàn)了相對可靠的預(yù)測。

3.4 綜合優(yōu)化

最后將三種優(yōu)化方式同時(shí)作用于LSTM 模型，得到的最終預(yù)測效果如圖3所示。

由圖3 可知，經(jīng)過優(yōu)化后預(yù)測效果得到較大改善，擬合度明顯上升。除去第33-35 個(gè)月實(shí)際儲(chǔ)量有較大幅度上升、預(yù)測數(shù)據(jù)沒有上升到相應(yīng)高度以外，其余月份中兩條曲線幾乎完全重合。計(jì)算所得殘差全部控制在0.6萬噸以內(nèi)、相對誤差基本控制在5%，只有極個(gè)別值臨近10%。作為預(yù)測模型，這樣的誤差完全可以接受，即已經(jīng)達(dá)到了一個(gè)較為滿意的預(yù)測水準(zhǔn)。

圖3 完全優(yōu)化后LSTM預(yù)測結(jié)果圖Fig.3 The prediction of the completely optimized LSTM

將改進(jìn)前后的預(yù)測結(jié)果進(jìn)行量化對比，如表5所示。

表5 優(yōu)化前后預(yù)測結(jié)果對比Tab.5 The prediction comparison of initial and optimized LSTM

由表5 可以看出，預(yù)測效果較之前改善了很多。尤其是當(dāng)后一個(gè)月比前一個(gè)月出現(xiàn)較大變化時(shí)，預(yù)測值不僅能與實(shí)際值保持同增同減，增減的幅度也近似相同。第5、16、26 個(gè)月這三個(gè)下降值較大的月份，變化程度甚至完全相同。在所有29 個(gè)月中，僅有3 個(gè)月的趨勢與實(shí)際值不同，即預(yù)測趨勢的準(zhǔn)確值達(dá)到了89.66%，因?yàn)槌霈F(xiàn)偏差的個(gè)別數(shù)值即使趨勢不同，但相差也很小，可以忽略，也就是說達(dá)到了企業(yè)的基本運(yùn)營要求。

4 總結(jié)

LNG儲(chǔ)量受到匯率、政策、國家基礎(chǔ)建設(shè)等諸多因素的影響。如何從已有信息中分析出有效信息輔助未來儲(chǔ)量的預(yù)測，是每個(gè)投資者都想要尋找的辦法。本文通過互信息的方法分析了LNG 影響因素的相關(guān)性和冗余性，建立了基于LSTM 神經(jīng)網(wǎng)絡(luò)的LNG 預(yù)測模型，并根據(jù)預(yù)測結(jié)果對LSTM 預(yù)測模型進(jìn)行了優(yōu)化，使其預(yù)測結(jié)果控制在較小的誤差范圍內(nèi)，從而為LNG 儲(chǔ)量計(jì)劃提供指導(dǎo)。