劉怡然, 王東杰, 鄧雪峰, 劉振宇
(1. 山西農(nóng)業(yè)大學(xué) 信息與電氣工程學(xué)院, 山西 晉中 030801; 2. 中國農(nóng)業(yè)科學(xué)院 農(nóng)業(yè)信息研究所, 北京 100081)
近年來,生豬價(jià)格市場的異常波動(dòng)受到社會(huì)輿論的廣泛關(guān)注,及時(shí)準(zhǔn)確地掌握生豬價(jià)格的變化趨勢對確保農(nóng)產(chǎn)品市場平穩(wěn)健康發(fā)展有重要意義.如果能夠?qū)ιi價(jià)格進(jìn)行預(yù)測,決策者能夠合理地調(diào)控生豬市場價(jià)格,養(yǎng)殖戶也可據(jù)此調(diào)整圈舍中仔豬、種豬和育肥豬存欄量,因此,生豬價(jià)格預(yù)測成為研究者密切關(guān)注的領(lǐng)域.
早期的生豬價(jià)格預(yù)測模型結(jié)構(gòu)相對單一,模型泛化能力較差,預(yù)測精度難以保證[1].隨著計(jì)算機(jī)計(jì)算能力的提升,生豬價(jià)格預(yù)測模型的復(fù)雜程度越來越高,從單一的灰色系統(tǒng)模型、向量自回歸模型、神經(jīng)網(wǎng)絡(luò)模型等趨于日漸復(fù)雜的組合模型[2-3].LI Z. M.等[4]以生豬歷史價(jià)格為依據(jù),將遺傳算法優(yōu)化的混沌神經(jīng)網(wǎng)絡(luò)用于生豬價(jià)格的短期預(yù)測,能夠預(yù)測20 d的生豬價(jià)格.任青山等[5]將多元回歸分析和反向傳播神經(jīng)網(wǎng)絡(luò)(back propagation neural network, BPNN)結(jié)合,提高了生豬價(jià)格預(yù)測的準(zhǔn)確性和可靠性.然而,生豬價(jià)格存在偽周期現(xiàn)象,即生豬價(jià)格周期是重復(fù)出現(xiàn)的,但又不是完全重復(fù),采用固定序列長度對生豬價(jià)格進(jìn)行預(yù)測并不十分合適[6-7].為了解決這一問題,姜百臣等[8]引入集成經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition, EEMD)方法剖析“豬周期”價(jià)格波動(dòng)原理,并結(jié)合遺傳算法(GA)優(yōu)化的支持向量機(jī)(SVM)建立了生豬價(jià)格預(yù)測模型.LIU Y. R.等[9]設(shè)計(jì)了一種相似子序列搜索算法并將其與支持向量回歸模型組合,排除了生豬價(jià)格周期成分在時(shí)間軸上的彎曲和伸縮造成的干擾,解決生豬價(jià)格預(yù)測過程中序列周期長度不固定的問題.
隨著生豬價(jià)格及其影響因素?cái)?shù)據(jù)量的積累,可以采用深度學(xué)習(xí)的方法解決這一問題.LSTM神經(jīng)網(wǎng)絡(luò)使循環(huán)神經(jīng)網(wǎng)絡(luò)(recurrent neural network, RNN)自循環(huán)的權(quán)重視上下文而定,而不是固定的.由于LSTM積累的時(shí)間尺度可以動(dòng)態(tài)改變,非常適合解決生豬價(jià)格序列預(yù)測問題[10-11].
因此,針對生豬價(jià)格的偽周期現(xiàn)象,文中提出螢火蟲算法優(yōu)化長短時(shí)記憶神經(jīng)網(wǎng)絡(luò)(firefly algorithm-long short term memory, FA-LSTM)的生豬價(jià)格預(yù)測方法.首先采用螢火蟲算法搜索LSTM的參數(shù),再根據(jù)得到的最優(yōu)參數(shù)建立預(yù)測模型,將本算法與其他經(jīng)典機(jī)器學(xué)習(xí)預(yù)測算法進(jìn)行比較,具有更高的預(yù)測準(zhǔn)確度.
生豬價(jià)格序列屬于典型的非線性時(shí)間序列,且由于其周期變化的特點(diǎn),采用LSTM神經(jīng)網(wǎng)絡(luò)對其進(jìn)行預(yù)測,考慮到其關(guān)鍵參數(shù)對預(yù)測結(jié)果影響非常大,因此采用螢火蟲算法對LSTM模型進(jìn)行參數(shù)優(yōu)化.具體方法流程如圖1所示.
圖1 基于FA-LSTM的生豬價(jià)格預(yù)測方法總體流程
長短時(shí)記憶神經(jīng)網(wǎng)絡(luò)適用于解決非線性時(shí)間可變時(shí)間序列預(yù)測問題[12].LSTM的單元結(jié)構(gòu)如圖2所示,其關(guān)鍵是類似于“傳送帶”的單元狀態(tài)Ct(cell state),通過輸入門、遺忘門和輸出門的調(diào)節(jié),決定哪些信息可以送上“傳送帶”.解決了循環(huán)神經(jīng)網(wǎng)絡(luò)多次展開中導(dǎo)致的梯度爆炸和梯度消失問題,令其能夠處理長時(shí)間依賴[13-15].
圖2 LSTM的單元結(jié)構(gòu)
遺忘門決定上一時(shí)刻的單元狀態(tài)信息有多少可以保留到當(dāng)前時(shí)刻,則
ft=σ(Wf[ht-1,xt]+bf),
(1)
式中:ft表示遺忘門;ht-1為上一個(gè)單元狀態(tài)的輸出,它與當(dāng)前時(shí)刻的輸入xt共同拼接成一個(gè)輸入矩陣;Wf和bf分別表示遺忘門的權(quán)重矩陣和偏置;σ為sigmoid激活函數(shù).
it=σ(Wi[ht-1,xt]+bi),
(2)
(3)
(4)
式中:W和b分別表示相應(yīng)的權(quán)重矩陣和偏置.
根據(jù)當(dāng)前單元狀態(tài),通過輸出門ot決定當(dāng)前輸出信息.
ot=σ(Wo[ht-1,xt]+bo),
(5)
ht=ottanhCt.
(6)
進(jìn)行LSTM訓(xùn)練時(shí),采用均方誤差(mean squared error,MSE)作為損失函數(shù),其定義如下:
(7)
訓(xùn)練LSTM時(shí),由公式(1)-(6)得到單元狀態(tài)輸出后,由公式(7)可以得到存儲(chǔ)單元的誤差項(xiàng),根據(jù)誤差項(xiàng)更新權(quán)重梯度和權(quán)重,文中采用自適應(yīng)動(dòng)量估計(jì)算法(adaptive moment estimation, ADAM)更新權(quán)重[16].
LSTM神經(jīng)網(wǎng)絡(luò)需要調(diào)節(jié)的參數(shù)較多,主要包括隱含層數(shù)目、隱含層神經(jīng)元數(shù)、學(xué)習(xí)率、批次大小等參數(shù)等.因此,需要結(jié)合參數(shù)搜索算法得到合適的參數(shù)進(jìn)行網(wǎng)絡(luò)設(shè)置優(yōu)化.與粒子群優(yōu)化(partical swarm optimization, PSO)等算法相比,螢火蟲算法在解決嘈雜的非線性優(yōu)化問題方面有優(yōu)勢,成為一種有利的優(yōu)化工具[16-18].
螢火蟲算法中,搜索過程與螢火蟲的相對熒光亮度I和相互吸引度β有關(guān),其定義如公式(8)和公式(9)所示.發(fā)光亮的螢火蟲會(huì)吸引發(fā)光弱的螢火蟲向它移動(dòng),發(fā)光越亮代表其位置越好,最亮螢火蟲代表函數(shù)最優(yōu)解[16].
I=I0e-γrij,
(8)
(9)
式中:I0表示最亮亮度;γ表示光吸收系數(shù),通常取γ=1;r為兩只螢火蟲之間的距離;β0為其初始吸引度.每只螢火蟲將會(huì)朝著所有亮度比自己高的螢火蟲移動(dòng),其新的位置Xi+1如下:
Xi+1=Xi+β(r)(Xi-Xj)+αrand,
(10)
式中:Xi和Xj分別表示i,j2只螢火蟲的空間位置;α為步長;rand為隨機(jī)移動(dòng).
為了對生豬價(jià)格進(jìn)行預(yù)測,設(shè)計(jì)了一個(gè)具有2層的LSTM網(wǎng)絡(luò)結(jié)構(gòu),將螢火蟲算法作為目標(biāo)規(guī)劃算法輔助尋找LSTM的第1、2隱含層神經(jīng)元個(gè)數(shù)、學(xué)習(xí)率和批次大小.
構(gòu)建FA優(yōu)化LSTM的預(yù)測模型如圖3所示.預(yù)測模型的主要步驟如下: ① 初始化螢火蟲算法的參數(shù)和目標(biāo)參數(shù)值域,主要包括迭代次數(shù)、種群大小、光吸收系數(shù)γ、初始吸引度β0和步長α,目標(biāo)參數(shù)值域包括LSTM的第1、2隱含神經(jīng)元個(gè)數(shù)、學(xué)習(xí)率和批次大小的取值范圍.② 在目標(biāo)參數(shù)的取值范圍內(nèi),隨機(jī)初始化一個(gè)螢火蟲種群,即目標(biāo)參數(shù)組合.③ 將訓(xùn)練集數(shù)據(jù)分割為訓(xùn)練集和驗(yàn)證集,使用初始的目標(biāo)參數(shù)訓(xùn)練得到LSTM模型,并通過驗(yàn)證集數(shù)據(jù)計(jì)算初始適應(yīng)度集合,選取其中最佳者作為最佳適應(yīng)度.④ 判斷是否到達(dá)最大迭代次數(shù),若達(dá)到,則輸出有最優(yōu)適應(yīng)度相對應(yīng)的目標(biāo)參數(shù)組合,使用最優(yōu)目標(biāo)參數(shù)訓(xùn)練得到LSTM模型,若未達(dá)到,則進(jìn)行步驟⑤、⑥.⑤ 計(jì)算螢火蟲種群亮度,并使每一只螢火蟲向更亮的螢火蟲移動(dòng).⑥ 使用新的種群訓(xùn)練新的LSTM模型并計(jì)算新種群的適應(yīng)度,判斷是否有螢火蟲的適應(yīng)度大于最佳適應(yīng)度,無則直接返回步驟④,有則更新最佳適應(yīng)度,并記錄該適應(yīng)度對應(yīng)目標(biāo)參數(shù)組合,返回步驟④.
圖3 FA優(yōu)化LSTM的生豬價(jià)格預(yù)測模型建立流程圖
為了滿足LSTM模型訓(xùn)練的需要,采集中國種豬信息網(wǎng)、豬價(jià)格網(wǎng)和山東畜牧獸醫(yī)網(wǎng)等農(nóng)產(chǎn)品價(jià)格網(wǎng)站的生豬價(jià)格、豬肉價(jià)格、仔豬價(jià)格、玉米價(jià)格和豆粕價(jià)格作為試驗(yàn)數(shù)據(jù),相關(guān)數(shù)據(jù)采集情況如表1所示.
表1 數(shù)據(jù)采集情況描述
爬蟲程序采集的數(shù)據(jù)在量的單位、采樣頻率和精確度上存在差異,所以需要經(jīng)過預(yù)處理操作得到正確完整價(jià)格序列,以便進(jìn)行建模和預(yù)測.
1) 數(shù)據(jù)清洗.首先檢測缺失值和異常值,對前后時(shí)間間隔不大的缺失數(shù)據(jù)采用線性插值進(jìn)行了填補(bǔ),如式(11)所示,對出現(xiàn)的異常數(shù)據(jù)采用均值平滑法進(jìn)行處理,如式(12)所示.
(11)
(12)
式中:xa+i為a+i時(shí)刻的缺失值;xb為b時(shí)刻的異常值;xa、xa+j、xb+i、xb-i為分別為a、a+j、b+i、b-i時(shí)刻的有效數(shù)據(jù).
2) 數(shù)據(jù)變換.豬價(jià)格網(wǎng)的數(shù)據(jù)采集范圍是全國各縣,以省為單位進(jìn)行統(tǒng)計(jì),得到每日全國各省的平均值.豬價(jià)格網(wǎng)的生豬、豬肉和仔豬的價(jià)格采樣頻率是每日1次,然而豆粕價(jià)格和玉米價(jià)格并非如此,且為了與其他2個(gè)數(shù)據(jù)來源的數(shù)據(jù)保持一致,采用周平均值將豬價(jià)格網(wǎng)的價(jià)格序列變換為每周1個(gè)樣本.經(jīng)過數(shù)據(jù)變換后,將3個(gè)數(shù)據(jù)來源的樣本整合為一個(gè)數(shù)據(jù)集,得到樣本2 131個(gè),其基本情況描述如表2所示.
表2 數(shù)據(jù)集基本情況描述 元·kg-1
3) 數(shù)據(jù)規(guī)約.因各商品價(jià)格數(shù)量級有差異,采用式(13)對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理.
(13)
式中:xi表示原始值,i表示其序號(hào);xm表示序列平均值;xstd表示序列標(biāo)準(zhǔn)差.
4) 參數(shù)選擇.研究[19-20]表明,生豬、豬肉、仔豬、玉米和豆粕等相關(guān)因素的變動(dòng)都有可能影響生豬價(jià)格,然而并不是所有變量都能夠?qū)ιi價(jià)格預(yù)測產(chǎn)生積極影響.為了避免無關(guān)變量干擾預(yù)測精度,采用主成分分析法(principal component analysis, PCA)和相關(guān)系數(shù)分析對預(yù)測參數(shù)進(jìn)行選取[21-22].由主成分分析得到各個(gè)主成分的方差貢獻(xiàn)率如下:成分1為0.709 8;成分2為0.229 9;成分3為0.055 0;成分4為0.005 2,前2個(gè)主成分貢獻(xiàn)了約94%方差,因此選取這2個(gè)主成分得到成分矩陣如表3所示.各因素相關(guān)系數(shù)如下:生豬價(jià)格,1.000 0;豬肉價(jià)格,0.986 4;仔豬價(jià)格,0.894 7;玉米價(jià)格,0.255 6;豆粕價(jià)格,0.113 2.
表3 成分矩陣
生豬價(jià)格、豬肉價(jià)格和仔豬價(jià)格都對成分1貢獻(xiàn)較大,而成分1有70%的方差貢獻(xiàn)率;生豬價(jià)格與豬肉價(jià)格、仔豬價(jià)格的相關(guān)性都較高,而與玉米、豆粕的相關(guān)性較低.因此,選取生豬價(jià)格、豬肉價(jià)格和仔豬價(jià)格作為預(yù)測生豬價(jià)格的參數(shù).
假設(shè)有訓(xùn)練樣本集S={X,Y},X為樣本特征向量(矩陣),Y為標(biāo)簽,樣本集中任意一個(gè)樣本點(diǎn)i表示為(X(i),Y(i)).若訓(xùn)練集與標(biāo)簽之間有一個(gè)長度為l的時(shí)間步長,每個(gè)樣本有長度為m的時(shí)間窗口,則有如下表示的樣本點(diǎn)(X(i),Y(i)).
1≤i≤n-m-l+1,i∈N,
(14)
為了驗(yàn)證本方法在進(jìn)行生豬價(jià)格預(yù)測時(shí)的優(yōu)越性,將本模型產(chǎn)生的結(jié)果與BPNN模型和SVR模型等兩種經(jīng)典的淺層預(yù)測模型進(jìn)行了對比,同時(shí)也將本模型與另一種循環(huán)神經(jīng)網(wǎng)絡(luò)門控循環(huán)單元(gated recurrent unit, GRU)模型進(jìn)了對比[23-26].
試驗(yàn)平臺(tái)硬件配置為Intel(R) Core(TM)i5-7360U CPU@2.3GHz處理器,8 GB運(yùn)行內(nèi)存;試驗(yàn)所采用的軟件工具為python編程語言和keras框架.
參數(shù)優(yōu)化試驗(yàn)中,螢火蟲算法的迭代次數(shù)設(shè)置為200、種群大小為20、光吸收系數(shù)γ=1、初始吸引度β0=1,步長α=0.5[17];LSTM的第1、2隱含神經(jīng)元個(gè)數(shù)的搜索范圍設(shè)置在[1, 100],批次大小的取值范圍是[16, 64];LSTM采用Adam算法進(jìn)行自適應(yīng)優(yōu)化,取值范圍[0.001,0.1].預(yù)測試驗(yàn)中,時(shí)間步長l分別取值1、4和8,每次試驗(yàn)的結(jié)果都與單隱含層的LSTM、SVR、BPNN和GRU進(jìn)行對比,SVR、BPNN和GRU的參數(shù)也是由螢火蟲算法搜索得出的,時(shí)間窗口m=16.由于對生豬價(jià)格的預(yù)測屬于時(shí)序數(shù)據(jù)預(yù)測,訓(xùn)練集和測試集不宜隨機(jī)選取,故取數(shù)據(jù)集的前90%樣本作為訓(xùn)練集,后10%樣本作為測試集.
采用平均絕對誤差(mean absolute error,MAE)、均方根誤差(root mean squared error,RMSE)和確定系數(shù)(Rsquared,R2)評估模型的準(zhǔn)確性,計(jì)算公式如下:
(15)
(16)
(17)
3.2.1 螢火蟲算法參數(shù)搜索結(jié)果
時(shí)間步長為1時(shí),利用螢火蟲算法對LSTM、SVR、BPNN和GRU模型進(jìn)行參數(shù)搜索的適應(yīng)度函數(shù)值變化如圖4所示.FA-LSTM表示螢火蟲算法優(yōu)化的LSTM模型,其他同理.
圖4 適應(yīng)度函數(shù)值變化
由圖4可知,螢火蟲算法能在較少的迭代次數(shù)下迅速達(dá)到最佳適應(yīng)度,非常適合進(jìn)行參數(shù)搜索.文中采用的適應(yīng)度函數(shù)是均方根誤差的倒數(shù),從適應(yīng)度曲線的取值可以看出,LSTM模型有較小的訓(xùn)練誤差.
為了證明LSTM模型適合對生豬價(jià)格的預(yù)測,排除參數(shù)搜索算法不同或者不進(jìn)行參數(shù)搜索對模型預(yù)測準(zhǔn)確性的干擾,也采用螢火蟲算法對其他預(yù)測模型的參數(shù)進(jìn)行了搜索,表4為不同時(shí)間步長下各個(gè)模型所取得的參數(shù).
表4 不同時(shí)間步長下各個(gè)模型所取得的參數(shù)
3.2.2 對生豬價(jià)格的預(yù)測結(jié)果
時(shí)間步長分別取1、4和8時(shí), 3次生豬價(jià)格預(yù)測結(jié)果如圖5所示,可以看出FA優(yōu)化的LSTM模型對真實(shí)價(jià)格的擬合度較好.
圖5 3次生豬價(jià)格預(yù)測結(jié)果
時(shí)間步長為1、4和8時(shí),F(xiàn)A-LSTM模型、LSTM模型、FA-SVR模型、FA-BPNN模型和FA-GRU模型對生豬價(jià)格預(yù)測結(jié)果與真實(shí)生豬價(jià)格的平均絕對誤差、均方根誤差和確定系數(shù)如表5所示.
表5 3次生豬價(jià)格預(yù)測模型評價(jià)參數(shù)
由表5可知,在時(shí)間步長為1、4和8時(shí),F(xiàn)A優(yōu)化LSTM都能取得比較好的預(yù)測結(jié)果,尤其是在時(shí)間步長為1時(shí),F(xiàn)A優(yōu)化LSTM模型的確定系數(shù)為92.57%,說明了其對真實(shí)價(jià)格走勢的擬合度較高.當(dāng)時(shí)間步長為1時(shí),與不經(jīng)過螢火蟲算法優(yōu)化的LSTM模型(其參數(shù)設(shè)置與GRU模型相同)相比,平均絕對誤差和均方根誤差分別降低了20.03%和36.71%,確定系數(shù)增加了4.88%,當(dāng)時(shí)間步長取4、8時(shí),與不經(jīng)過螢火蟲算法搜尋參數(shù)的LSTM模型相比,平均絕對誤差和均方根誤差也都下降了,而確定系數(shù)也都增加了,說明參數(shù)的選擇對LSTM預(yù)測準(zhǔn)確率影響非常大,但是人工進(jìn)行參數(shù)搜索效率太低,螢火蟲算法大大提升了該模型的參數(shù)搜索效率.
由表5也可看出,在時(shí)間步長為1、4和8時(shí),各個(gè)模型的表現(xiàn)都較為不穩(wěn)定.FA優(yōu)化的SVR模型在時(shí)間步長取1時(shí),對生豬價(jià)格的預(yù)測表現(xiàn)比較好,僅次于FA優(yōu)化的LSTM模型,然而當(dāng)時(shí)間步長取值為8時(shí),它的表現(xiàn)又成為幾個(gè)模型中最差的.與之相反的是,F(xiàn)A優(yōu)化的GRU模型在時(shí)間步長取1時(shí)表現(xiàn)較差,僅優(yōu)于FA優(yōu)化的BPNN模型的表現(xiàn),但當(dāng)時(shí)間步長取8時(shí),其表現(xiàn)能夠和FA優(yōu)化的LSTM旗鼓相當(dāng).但FA優(yōu)化LSTM模型在3次試驗(yàn)中表現(xiàn)均優(yōu)于其他模型,說明FA優(yōu)化的LSTM的生豬價(jià)格預(yù)測模型具有較為理想的泛化性能,可以很好地?cái)M合非線性變化的生豬價(jià)格.
雖然LSTM每個(gè)時(shí)間步長和權(quán)重的計(jì)算復(fù)雜度為O(1),且其每個(gè)時(shí)間步上權(quán)重是共享的,但螢火蟲算法、粒子群算法、人工蜂群等一系列啟發(fā)式算法都是NP-hard問題,導(dǎo)致現(xiàn)有試驗(yàn)條件下參數(shù)搜索階段時(shí)間仍較長.其他相似的研究也存在這樣的問題,在試驗(yàn)的后期使用阿里云服務(wù)器ECS,縮短了程序執(zhí)行時(shí)間,所使用配置為Intel Xeon Platinum 8163處理器(8核,2.5 GHz),16 GB內(nèi)存.由于訓(xùn)練數(shù)據(jù)是相同的,可以從目標(biāo)參數(shù)搜索范圍和步長設(shè)置上對螢火蟲算法進(jìn)行分級,采用分布式的參數(shù)搜索方式也可以對螢火蟲算法進(jìn)行自適應(yīng)步長的改進(jìn).
1) 考慮到生豬價(jià)格序列在時(shí)間軸上的遲滯問題,提出了一種基于長短時(shí)記憶神經(jīng)網(wǎng)絡(luò)算法的生豬價(jià)格預(yù)測模型,首先對生豬價(jià)格數(shù)據(jù)進(jìn)行預(yù)處理操作,并采用螢火蟲算法對長短時(shí)記憶神經(jīng)網(wǎng)絡(luò)進(jìn)行了參數(shù)優(yōu)化,能夠快速找到合適的長短時(shí)記憶神經(jīng)網(wǎng)絡(luò)參數(shù),再利用經(jīng)過優(yōu)化的參數(shù)建立LSTM預(yù)測模型.
2) 所建立的基于長短時(shí)記憶神經(jīng)網(wǎng)絡(luò)的生豬價(jià)格預(yù)測模型,能夠預(yù)測未來1周、1月和2月的生豬價(jià)格,與傳統(tǒng)的淺層預(yù)測模型和未經(jīng)優(yōu)化的長短時(shí)記憶神經(jīng)網(wǎng)絡(luò)模型相比,提出的預(yù)測模型在不同時(shí)間間隔內(nèi)的預(yù)測結(jié)果都具有更高的預(yù)測準(zhǔn)確度.
3) 文中提出的生豬價(jià)格預(yù)測模型具有良好的預(yù)測性能和泛化能力,不僅能夠?yàn)橹贫ǚ€(wěn)定生豬市場的決策提供量化分析工具,也能為其他農(nóng)產(chǎn)品價(jià)格的預(yù)測提供參考.