劉半藤，陳 唯，尹則高，孫 萍

(1. 常州大學(xué)信息科學(xué)與工程學(xué)院，江蘇 常州 213164；2. 中國(guó)海洋大學(xué)水利工程學(xué)院，山東 青島 266003；3. 浙江樹(shù)人大學(xué)信息科技學(xué)院，浙江 杭州 310015)

1 引言

時(shí)間序列預(yù)測(cè)指的是通過(guò)已知的歷史數(shù)據(jù)去預(yù)測(cè)未知的未來(lái)數(shù)據(jù)，已經(jīng)在醫(yī)學(xué)、金融、軍事等[1-3]領(lǐng)域中得到了廣泛的應(yīng)用，由于實(shí)際生活中的時(shí)間序列大多呈現(xiàn)非線性和不穩(wěn)定性，因此對(duì)于非線性和不穩(wěn)定時(shí)間序列的預(yù)測(cè)問(wèn)題一直備受各個(gè)領(lǐng)域中研究學(xué)者的關(guān)注[4]。

目前，對(duì)于非線性和不穩(wěn)定時(shí)間序列進(jìn)行預(yù)測(cè)的主要方法是采用回聲狀態(tài)網(wǎng)絡(luò)(Echo State Network，ESN)[5]。其特點(diǎn)就是采用被稱之為“儲(chǔ)備池”的大規(guī)模隨機(jī)稀疏連接網(wǎng)絡(luò)作為傳統(tǒng)神經(jīng)網(wǎng)絡(luò)的隱層，去處理非線性和不穩(wěn)定的時(shí)間序列，并且在訓(xùn)練過(guò)程中只需訓(xùn)練儲(chǔ)備池至輸出層的輸出權(quán)值，簡(jiǎn)化了網(wǎng)絡(luò)的訓(xùn)練過(guò)程，避免了傳統(tǒng)神經(jīng)網(wǎng)絡(luò)中存在的易于陷入局部最優(yōu)、訓(xùn)練算法復(fù)雜等問(wèn)題[6]。目前，一部分學(xué)者為了提高ESN的性能，紛紛開(kāi)始研究改進(jìn)ESN的方法，主要集中在儲(chǔ)備池拓?fù)浣Y(jié)構(gòu)優(yōu)化和輸出權(quán)值優(yōu)化兩個(gè)方面[7-8]。在儲(chǔ)備池拓?fù)浣Y(jié)構(gòu)方面，傳統(tǒng)ESN的儲(chǔ)備池中采用隨機(jī)網(wǎng)絡(luò)致使模型訓(xùn)練無(wú)目的性，為了解決此問(wèn)題，李菡[9]等學(xué)者提出采用同時(shí)具備隨機(jī)性和規(guī)則性的小世界回聲狀態(tài)網(wǎng)絡(luò)作為ESN的儲(chǔ)備池去預(yù)測(cè)非線性時(shí)間序列，提高了預(yù)測(cè)模型的適應(yīng)性和預(yù)測(cè)精度；但以上小世界回聲狀態(tài)網(wǎng)絡(luò)的節(jié)點(diǎn)連接屬于確定性連接，對(duì)具有時(shí)變性和模糊性的時(shí)間序列預(yù)測(cè)精度不高，因此，倫淑嫻[10]等學(xué)者提出采用改進(jìn)的小世界回聲狀態(tài)網(wǎng)絡(luò)作為ESN的儲(chǔ)備池，即儲(chǔ)備池網(wǎng)絡(luò)節(jié)點(diǎn)間的加邊概率根據(jù)節(jié)點(diǎn)間距離的負(fù)指數(shù)函數(shù)進(jìn)行修正，縮短了訓(xùn)練時(shí)間并提高了對(duì)具有時(shí)變性和模糊性的非線性時(shí)間序列的預(yù)測(cè)精度。

在輸出權(quán)值方面，傳統(tǒng)ESN的輸出權(quán)值計(jì)算采用偽逆法，但偽逆法求解高維線性回歸時(shí)易出現(xiàn)共線性問(wèn)題[11]；為了解決此問(wèn)題，Wang[12-13]等學(xué)者提出采用Ridge回歸、Lasso回歸等線性回歸方法計(jì)算輸出權(quán)值，通過(guò)添加L2范數(shù)、L1范數(shù)解決共線性問(wèn)題，但Ridge回歸、Lasso回歸對(duì)偏大的輸出權(quán)值施加了更大的懲罰，屬于有偏估計(jì)，模型預(yù)測(cè)時(shí)容易出現(xiàn)過(guò)擬合問(wèn)題[14]。為了解決以上問(wèn)題，在計(jì)算輸出權(quán)值時(shí)需要采用漸近無(wú)偏正則化方法，提高預(yù)測(cè)模型的預(yù)測(cè)精度與泛化性能。常見(jiàn)的漸近無(wú)偏正則化方法有SCAD(Smoothly Clipped Absolute Deviation)正則化法[15]和MCP(Minimax Concave Penalty)正則化法[16]，目前，已有學(xué)者將SCAD正則化法成功應(yīng)用于小世界回聲狀態(tài)網(wǎng)絡(luò)的輸出權(quán)值優(yōu)化中[17]，提高了小世界回聲狀態(tài)網(wǎng)絡(luò)對(duì)非線性時(shí)間序列的預(yù)測(cè)精度。但采用MCP正則化法優(yōu)化小世界回聲狀態(tài)網(wǎng)絡(luò)輸出權(quán)值的方法還未曾提出，而且MCP正則化法的懲罰函數(shù)具有最小最大凸性，能對(duì)偏大或偏小的輸出權(quán)值做更恰當(dāng)?shù)膽土P，更適用于處理多維非線性數(shù)據(jù)[18-19]，因此，本文提出了一種基于MCP正則化SWESN的時(shí)間序列預(yù)測(cè)方法(Minimax Concave Penalty-Small World Echo State Network， MCP-SWESN)，提高小世界回聲狀態(tài)網(wǎng)絡(luò)對(duì)非線性時(shí)間序列的預(yù)測(cè)能力。

2 基于MCP正則化SWESN的時(shí)間序列預(yù)測(cè)方法

本文采用改進(jìn)的小世界網(wǎng)絡(luò)作為儲(chǔ)備池，得到小世界回聲狀態(tài)網(wǎng)絡(luò)(Small World Echo State Network，SWESN)。其拓?fù)浣Y(jié)構(gòu)如圖1所示。

圖1 小世界回聲狀態(tài)網(wǎng)絡(luò)拓?fù)鋱D

小世界回聲狀態(tài)網(wǎng)絡(luò)的狀態(tài)方程和輸出方程分別為

x(t)=f(Winu(t)+Wxx(t-1))

(1)

y(t)=xT(t)Wout

(2)

其中，u(t)∈RL、x(t)∈RM和y(t)∈R分別表示儲(chǔ)備池t時(shí)刻的輸入變量、狀態(tài)變量和輸出變量；激活函數(shù)f通常取雙曲正切tanh函數(shù)；Win∈RM×L、Wx∈RM×M和Wout∈RM分別為輸入權(quán)值矩陣、儲(chǔ)備池內(nèi)部權(quán)值矩陣和輸出權(quán)值矩陣。輸入權(quán)值矩陣Win隨機(jī)生成，確定后不再改變。

改進(jìn)的小世界網(wǎng)絡(luò)中儲(chǔ)備池內(nèi)部權(quán)值矩陣Wx通過(guò)建立加邊概率與節(jié)點(diǎn)間距離的函數(shù)關(guān)系式獲得，確定后不再改變。加邊概率P值隨著節(jié)點(diǎn)間的距離增大按指數(shù)方式遞減，即：

P=α×e(-β×d)

(3)

其中，P取值范圍為[0， 1]，表示節(jié)點(diǎn)間的連接權(quán)值。d表示節(jié)點(diǎn)間的歐氏距離，α用來(lái)調(diào)節(jié)距離靈敏度，β用來(lái)調(diào)整網(wǎng)絡(luò)的整體密度。

輸出權(quán)值矩陣Wout在訓(xùn)練時(shí)獲得，即最小化目標(biāo)函數(shù)對(duì)應(yīng)的Wout值，如式(4)所示，并通過(guò)最小二乘法求解獲得，如式(5)所示：

(4)

Wout=X?Y=(XTX)-1XTY

(5)

其中，(X，Y)是訓(xùn)練樣本，X?是X的偽逆。

正則化方法是在最小化目標(biāo)函數(shù)的基礎(chǔ)上添加懲罰函數(shù)，最小化添加懲罰項(xiàng)的目標(biāo)函數(shù)對(duì)應(yīng)的Wout估計(jì)值如式(6)所示：

(6)

其中，J表示變量個(gè)數(shù)，ρλ，γ表示罰函數(shù)。

本文采用最小最大凹罰(MCP)作為懲罰函數(shù)，MCP罰函數(shù)在原點(diǎn)產(chǎn)生奇異值，能夠產(chǎn)生稀疏解。并且，在|θ|>γλ時(shí)，直接將變量置為零，滿足對(duì)變量θ的近似無(wú)偏估計(jì)，MCP罰函數(shù)如式(7)所示

(7)

其中，γ，λ為可調(diào)超參數(shù)(γ>2，λ>0)，采用遍歷手段獲得，θ為參數(shù)向量，本文中為輸出權(quán)值Wout。

(8)

out=arg min(‖Y-XWout‖2+

(9)

其中，D為Wout中非零元素個(gè)數(shù)，通過(guò)對(duì)式(9)重復(fù)執(zhí)行Ridge回歸解可得輸出權(quán)值估計(jì)

(10)

3 仿真研究

3.1 測(cè)試數(shù)據(jù)

為了驗(yàn)證本方法的有效性，分別用Lorenz混沌時(shí)間序列、Mackey-Glass混沌時(shí)間序列和實(shí)際的PM2.5濃度時(shí)間序列進(jìn)行測(cè)試。

1)Lorenz混沌時(shí)間序列生成方法如下

(11)

取a=10，b=28，c=8/3，x(0)=12，y(0)=2，z(0)=9，系統(tǒng)呈現(xiàn)混沌特性。利用四階龍格庫(kù)塔算法對(duì)Lorenz系統(tǒng)求得1000個(gè)時(shí)刻作為測(cè)試數(shù)據(jù)。

2)Mackey-Glass混沌時(shí)間序列生成方法如下

(12)

取a=0.2，b=0.1，c=10，τ=17，x(0)=1.2，t(0)=0，系統(tǒng)呈現(xiàn)混沌特性。利用四階龍格庫(kù)塔算法對(duì)Mackey-Glass系統(tǒng)求得1000個(gè)時(shí)刻作為測(cè)試數(shù)據(jù)。

3)實(shí)際時(shí)間序列采用北京市某年度的PM2.5濃度時(shí)間序列，數(shù)據(jù)來(lái)源為中國(guó)空氣質(zhì)量在線監(jiān)測(cè)分析平臺(tái)(https：∥www.aqistudy.cn/)。

3.2 仿真分析

通過(guò)將Lorenz混沌時(shí)間序列、Mackey-Glass混沌時(shí)間序列和實(shí)際PM2.5濃度時(shí)間序列歸一化、相空間重構(gòu)處理后進(jìn)行仿真測(cè)試和分析，分別取前500個(gè)時(shí)刻進(jìn)行訓(xùn)練，而后200個(gè)時(shí)刻進(jìn)行預(yù)測(cè)，測(cè)試結(jié)果進(jìn)行反歸一化。MCP-SWESN對(duì)Lorenz混沌時(shí)間序列、Mackey-Glass混沌時(shí)間序列和北京市某年度的PM2.5濃度時(shí)間序列的預(yù)測(cè)結(jié)果如圖2-圖4所示。

圖2 MCP-SWESN對(duì)Lorenz序列預(yù)測(cè)結(jié)果

圖3 MCP-SWESN對(duì)Mackey-Glass序列預(yù)測(cè)結(jié)果

圖4 MCP-SWESN對(duì)實(shí)際PM2.5濃度序列預(yù)測(cè)結(jié)果

圖2-圖4表明MCP-SWESN對(duì)Lorenz混沌時(shí)間序列、Mackey-Glass混沌時(shí)間序列和實(shí)際PM2.5濃度時(shí)間序列的預(yù)測(cè)均與實(shí)際曲線走勢(shì)一致，結(jié)果吻合度較好。

儲(chǔ)備池規(guī)模會(huì)影響ESN及其改進(jìn)方法的預(yù)測(cè)結(jié)果，因此為進(jìn)一步分析預(yù)測(cè)誤差，需要比較不同規(guī)模儲(chǔ)備池的ESN、SWESN、Ridge-SWESN、Lasso-SWESN、SCAD-SWESN和MCP-SWESN對(duì)Lorenz混沌時(shí)間序列、Mackey-Glass混沌時(shí)間序列和實(shí)際PM2.5濃度時(shí)間序列的預(yù)測(cè)結(jié)果，并采用標(biāo)準(zhǔn)均方根誤差(NRMSE)作為全部仿真預(yù)測(cè)的性能指標(biāo)

(13)

其中，Q表示預(yù)測(cè)時(shí)間序列長(zhǎng)度，Yd(t)表示目標(biāo)值，Y(t)表示預(yù)測(cè)值。avg(Yd(t))表示目標(biāo)值均值。

采用上述ESN及其改進(jìn)的六種方法分別對(duì)三組時(shí)間序列進(jìn)行仿真，重復(fù)30次，誤差結(jié)果取平均值。預(yù)測(cè)Lorenz混沌時(shí)間序列、Mackey-Glass混沌時(shí)間序列和實(shí)際PM2.5濃度時(shí)間序列的訓(xùn)練NRMSE、測(cè)試NRMSE結(jié)果分別如圖5-圖10所示。選擇最合適的儲(chǔ)備池規(guī)模所對(duì)應(yīng)的預(yù)測(cè)結(jié)果，結(jié)果如表1-表3所示。

圖5 不同方法的訓(xùn)練誤差結(jié)果(Lorenz)

圖6 不同方法的測(cè)試誤差結(jié)果(Lorenz)

表1 Lorenz混沌時(shí)間序列預(yù)測(cè)結(jié)果對(duì)比

圖7 不同方法的訓(xùn)練誤差結(jié)果(Mackey-Glass)

圖8 不同方法的測(cè)試誤差結(jié)果(Mackey-Glass)

表2 MG混沌時(shí)間序列預(yù)測(cè)結(jié)果對(duì)比

圖9 不同方法的訓(xùn)練誤差結(jié)果(PM2.5)

圖10 不同方法的訓(xùn)練誤差結(jié)果(PM2.5)

表3 實(shí)際PM2.5濃度時(shí)間序列預(yù)測(cè)結(jié)果對(duì)比

上述圖表展示了不同儲(chǔ)備池規(guī)模下，ESN、SWESN、Ridge-SWESN、Lasso-SWESN、SCAD-SWESN和MCP-SWESN六種預(yù)測(cè)方法對(duì)Lorenz混沌時(shí)間序列、Mackey-Glass混沌時(shí)間序列和實(shí)際PM2.5濃度時(shí)間序列進(jìn)行預(yù)測(cè)的訓(xùn)練NRMSE和測(cè)試NRMSE情況。MCP-SWESN方法在對(duì)上述三種時(shí)間序列進(jìn)行預(yù)測(cè)時(shí)測(cè)試NRMSE在六種模型中最小，預(yù)測(cè)精度最高；相比于ESN、SWESN和SCAD-SWESN方法，MCP-SWESN方法的訓(xùn)練NRMSE和測(cè)試NRMSE均較小，說(shuō)明MCP-SWESN預(yù)測(cè)方法的非線性擬合能力有了一定的提升；同時(shí)，相比于Ridge-SWESN、Lasso-SWESN預(yù)測(cè)方法，MCP-SWESN方法的訓(xùn)練NRMSE和測(cè)試NRMSE相近，說(shuō)明MCP-SWESN方法很好的解決了Ridge-SWESN和Lasso-SWESN預(yù)測(cè)模型出現(xiàn)的過(guò)擬合問(wèn)題。

4 結(jié)語(yǔ)

本文提出一種基于MCP正則化SWESN的時(shí)間序列預(yù)測(cè)方法，采用改進(jìn)的小世界回聲狀態(tài)網(wǎng)絡(luò)，縮短模型的訓(xùn)練時(shí)間并提高預(yù)測(cè)精度；采用MCP正則化方法優(yōu)化改進(jìn)的小世界回聲狀態(tài)網(wǎng)絡(luò)的輸出權(quán)值，解決常規(guī)回歸方法計(jì)算權(quán)值時(shí)出現(xiàn)的過(guò)擬合問(wèn)題，提高預(yù)測(cè)模型的非線性擬合能力。最后，仿真實(shí)現(xiàn)ESN、SWESN、Ridge-SWESN、Lasso-SWESN、SCAD-SWESN和MCP-SWESN六種預(yù)測(cè)方法對(duì)Lorenz混沌時(shí)間序列、Mackey-Glass混沌時(shí)間序列和北京市某年度的PM2.5濃度時(shí)間序列的預(yù)測(cè)，結(jié)果顯示：基于MCP-SWESN的時(shí)間序列預(yù)測(cè)方法具有更強(qiáng)的預(yù)測(cè)能力。