一種線性/非線性自回歸模型及其在建模和預(yù)測(cè)中的應(yīng)用

馬家欣 許飛云 黃 仁

(東南大學(xué)機(jī)械工程學(xué)院,南京 211189)

時(shí)間序列分析是現(xiàn)代系統(tǒng)辨識(shí)的重要方法之一,它在缺乏明確或者完全的系統(tǒng)輸入與輸出因果關(guān)系的情況下,將白噪聲看作系統(tǒng)輸入,從數(shù)理統(tǒng)計(jì)的角度揭示因果關(guān)系,最終求得系統(tǒng)的等價(jià)模型．該方法在各個(gè)領(lǐng)域都得到了廣泛的應(yīng)用[1-4]．

ARMA模型(包括AR模型和MA模型)是時(shí)序方法中最基本的時(shí)序模型,然而,它是在線性回歸模型的基礎(chǔ)上引申發(fā)展起來的,很難模擬工程實(shí)際中的非線性現(xiàn)象．而經(jīng)典的非線性時(shí)間序列模型,包括門限自回歸模型、雙線性模型、指數(shù)自回歸模型等,一般都是在特定的工程背景下提出的,通用性不夠．因此亟須提出一種適用于線性/非線性系統(tǒng)的時(shí)序模型．GNAR模型的提出和應(yīng)用,在形式上統(tǒng)一了不同背景下提出的非線性時(shí)序模型的表達(dá)式,并有機(jī)結(jié)合了線性模型[5-7]．另一方面,盡可能地利用更多的已知信息,獲取更多的系統(tǒng)特性,可進(jìn)一步提高建模精度．因此,在經(jīng)典時(shí)序模型的基礎(chǔ)上,把影響系統(tǒng)輸出的已知相關(guān)數(shù)據(jù)作為系統(tǒng)輸入,可得到帶有外部輸入的時(shí)序模型,如ARX,ARMAX等模型．這些模型不僅在理論上得到了驗(yàn)證,在工程實(shí)際中也得到了廣泛應(yīng)用[8-11]．本文針對(duì)系統(tǒng)輸入部分已知的特點(diǎn),在GNAR模型基礎(chǔ)上提出了GNARX模型,給出了模型的參數(shù)估計(jì)和結(jié)構(gòu)選取方法,并將該模型應(yīng)用于仿真和實(shí)際數(shù)據(jù),效果優(yōu)于其他模型．

1 GNARX模型表達(dá)式

GNAR模型建模時(shí),認(rèn)為系統(tǒng)輸入是零均值的白噪聲.如果系統(tǒng)已知單個(gè)輸入,記作序列{ut},系統(tǒng)的輸出序列記作{wt},白噪聲記作{at},那么GNAR模型就變形為單輸入的GNARX模型．

對(duì)于單輸入系統(tǒng),設(shè)t時(shí)刻輸出值為wt,表示為函數(shù)f的表達(dá)式,即

wt=f(wt-1,…,wt-∞,ut-τu,…,ut-τu-∞,

at-1,…,at-∞)+at

(1)

式中,f(·)為任意一般線性/非線性函數(shù);wt-i為t-i時(shí)刻的系統(tǒng)輸出觀測(cè)值,i=1,2,…;ut-τu-i為t-τu-i時(shí)刻的系統(tǒng)輸入,i=1,2,…;at-i為t-i時(shí)刻的白噪聲,i=1,2,…;τu為系統(tǒng)輸入ut的延遲．

對(duì)于一般線性平穩(wěn)系統(tǒng),式(1)中函數(shù)f表示ARMAX模型;而當(dāng)系統(tǒng)表現(xiàn)出非線性、非平穩(wěn)特征時(shí),由Weierstrass逼近原理[12]知,可用多項(xiàng)式逼近式(1)中的f函數(shù):

λi1,i2at-i1at-i2+γi1,i2wt-i1ut-τu+1-i2+ηi1,i2wt-i1at-i2+

πi1,i2ut-τu+1-i1at-i2)+…+at

(2)

式中,α1,α2,…,β1,β2,…,λ1,λ2,…是模型參數(shù)．

假設(shè)式(2)描述的系統(tǒng)具有零初始狀態(tài)(即wt=0,t≤0),且系統(tǒng)在初始時(shí)刻之前沒有輸入(即ut=0,t≤0)．實(shí)際建模中,模型階次p取有限值,令xt,i表示式(2)中的各i階項(xiàng),xt,i,1(i=1,2,…,p)表示i階項(xiàng)xt,i中的1次項(xiàng),即

xt,i,1={wt-1,wt-2,…,wt-nw,i,ut-τu,

ut-τu-1,…,ut-τu-nu,i+1}

(3)

則系統(tǒng)輸出wt可寫成如下形式:

(4)

式中,θ(i1),θ(i1,i2),…為模型參數(shù);xt,i,1(j)為向量xt,i,1中的第j個(gè)元素;nw,j(j=1,2,…,p)為輸出{wt}各階項(xiàng)的記憶步長(zhǎng);nu,j(j=1,2,…,p)為輸入{ut}各階項(xiàng)的記憶步長(zhǎng)．該模型簡(jiǎn)記為GNARX(p;τu;nw,1,nw,2,…,nw,p;nu,1,nu,2,…,nu,p)．

當(dāng)系統(tǒng)具有雙輸入{ut},{vt}時(shí),把式(4)推廣至雙輸入的GNARX模型,簡(jiǎn)記為GNARX(p;τu,τv;nw,1,nw,2,…,nw,p;nu,1,nu,2,…,nu,p;nv,1,nv,2,…,nv,p)．

將式(3)和(4)改寫成如下形式:

xt,i,1={wt-1,…,wt-nw,i,ut-τu,…,ut-τu-nu,i+1,

vt-τv,…,vt-τv-nv,i+1}

(5)

(6)

式中,nv,j(j=1,2,…,p)為輸入{vt}各階項(xiàng)的記憶步長(zhǎng);τv為系統(tǒng)輸入{vt}的延遲．

同理,式(6)可推廣到多輸入系統(tǒng),不再贅述．

2 GNARX模型的參數(shù)估計(jì)與結(jié)構(gòu)選取

GNARX模型的參數(shù)估計(jì)用最小二乘法很容易求得.以式(4)所示的單輸入GNARX模型為例,在時(shí)刻t,p階項(xiàng)記作xt,p,是一維向量形式．求取過程中用xt,p,i(i=1,2,…,p)表示p階項(xiàng)xt,p中的i次項(xiàng),也是一維向量,其中xt,p,i(j)是標(biāo)量,表示向量xt,p,i中的第j個(gè)元素．

xt,p,1={wt-1,wt-2,…,wt-nw,p,ut-τu,

ut-τu-1,…,ut-τu-nu,p+1}

xt,p,2={xt,p,1(1){xt,p,1(1)},

xt,p,1(2){xt,p,1(1),xt,p,1(2)},…,

xt,p,1(mp,1)xt,p,1}

?

xt,p,p={xt,p,p-1(1){xt,p,p-1(1)},

xt,p,p-1(2){xt,p,p-1(1),xt,p,p-1(2)},…,

xt,p,p-1(mp,p-1)xt,p,p-1}

(7)

令

xt,p=xt,p,p

(8)

記

xt={xt,1,xt,2,…,xt,p}
θ={θ(i1),θ(i1,i2),…}T

于是,t時(shí)刻到t+k時(shí)刻的方程式寫成

wt=xt·θ+at
wt+1=xt+1·θ+at+1
?
wt+k=xt+k·θ+at+k

θ參數(shù)最小二乘法估計(jì)為

(9)

其中

y={wt,wt+1,…,wt+k}

根據(jù)式(9)即可對(duì)單輸入GNARX模型進(jìn)行參數(shù)估計(jì)．同理,可獲得雙輸入、多輸入模型的參數(shù)估計(jì)方法．

GNARX模型結(jié)構(gòu)的辨識(shí),包括模型階次p的確定和記憶長(zhǎng)度的確定,這直接影響著模型的建模性能和預(yù)測(cè)性能及模型的復(fù)雜程度,同時(shí)也關(guān)系到計(jì)算速度、預(yù)報(bào)實(shí)時(shí)性等問題．

對(duì)于GNARX模型,定義單輸入模型i階記憶步長(zhǎng)si=nw,i+nu,i(i=1,2,…,p),雙輸入模型i階記憶步長(zhǎng)si=nw,i+nu,i+nv,i(i=1,2,…,p),多輸入模型i階記憶步長(zhǎng)可依此類推得到．則i階參數(shù)量ri為

式中,i=1,2,…,p．

則GNARX模型參數(shù)數(shù)量為

R=r1+r2+…+rp

修正AIC準(zhǔn)則函數(shù)[6]如下所示:

(10)

模型結(jié)構(gòu)的最終確定,需要綜合考慮模型的建模能力、預(yù)測(cè)能力和模型復(fù)雜程度等因素．在具體建模過程中,可根據(jù)現(xiàn)場(chǎng)要求的不同,調(diào)整式(10)中的重要度權(quán)系數(shù),來滿足實(shí)際建模需求．如模型用于濾波,則可增加建模重要度權(quán)系數(shù)α;若模型以預(yù)報(bào)為主,則可以增加預(yù)測(cè)重要度權(quán)系數(shù)β;而當(dāng)現(xiàn)場(chǎng)側(cè)重于建模實(shí)時(shí)性時(shí),則可以增加模型復(fù)雜重要度權(quán)系數(shù)γ．

3 GNARX模型的建模和預(yù)測(cè)應(yīng)用

GNARX模型適用范圍廣,對(duì)線性和非線性數(shù)據(jù)有著良好的建模預(yù)測(cè)能力．下文用GNARX模型對(duì)仿真數(shù)據(jù)和實(shí)際數(shù)據(jù)進(jìn)行建模和預(yù)測(cè),并與其他模型預(yù)測(cè)效果進(jìn)行對(duì)比．

為比較建模預(yù)測(cè)效果,定義如下衡量指標(biāo)：

建模均方誤差

(11)

預(yù)測(cè)均方誤差

(12)

建模平均絕對(duì)百分比誤差

(13)

預(yù)測(cè)平均絕對(duì)百分比誤差

(14)

3.1 對(duì)仿真數(shù)據(jù)的建模和預(yù)測(cè)

仿真數(shù)據(jù)模型具體形式如下:

ARX模型

xt-1.5xt-1+0.7xt-2=ut-3+0.5ut-4+at

(15)

ARMAX模型

xt-1.5xt-1+0.7xt-2=ut-1+0.5ut-2+

at-0.3at-1+0.2at-2-0.7at-3

(16)

NLARX模型

xt=F(xt-1,xt-2,ut-1,ut-2)+at

(17)

式中,at為白噪聲,服從正態(tài)分布N(0,0.1);ut-i(i=1,2,3,…)為系統(tǒng)輸入，服從平均分布U(0,1);xt-i(i=0,1,2,…)為系統(tǒng)輸出;F(·)取Matlab中S形網(wǎng)絡(luò)非線性估計(jì)函數(shù),單元數(shù)10,S形函數(shù)表達(dá)式為f(z)=1/(1+exp(-z))．

仿真數(shù)據(jù)選用前,先對(duì){xt}序列進(jìn)行歸一化處理,再舍去前50個(gè)過渡區(qū)數(shù)據(jù),取中間250個(gè)數(shù)據(jù)用于建模,后100個(gè)數(shù)據(jù)用于預(yù)測(cè)．分別用AR模型、GNAR模型、ARX模型、BP網(wǎng)絡(luò)和GNARX模型對(duì)上述各仿真數(shù)據(jù)進(jìn)行建模預(yù)測(cè)．其中時(shí)序模型是在合適的模型結(jié)構(gòu)選取范圍內(nèi)選取的,并用式(10)所示修正AIC判斷,重要度權(quán)系數(shù)嘗試選取α=1,β=1,γ=1,模型具體選取范圍見表1．

表1 針對(duì)不同數(shù)據(jù)建模時(shí)各模型結(jié)構(gòu)的選取范圍

表1中,BP神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)包括輸入層、隱層和輸出層3層,其中輸出層節(jié)點(diǎn)數(shù)mo為1,隱層節(jié)點(diǎn)數(shù)待定．對(duì)于網(wǎng)絡(luò)的輸入,可分為2種情況:

1) 僅以建模預(yù)測(cè)數(shù)據(jù)xt作為網(wǎng)絡(luò)的輸入,輸入樣本數(shù)據(jù)重合度記作r,輸入層節(jié)點(diǎn)數(shù)記作mi,網(wǎng)絡(luò)簡(jiǎn)記為BP1(mi;mh;mo;r)．

2) 數(shù)據(jù)xt,ut同時(shí)作為網(wǎng)絡(luò)的輸入,xt輸入維數(shù)為mx,ut輸入維數(shù)為mu,延遲為τu,輸入數(shù)據(jù)xt的重合度記作r,此時(shí)輸入層節(jié)點(diǎn)數(shù)為mx+mu,網(wǎng)絡(luò)簡(jiǎn)記為BP2(mx;mu,τu;mh;mo;r)．

其中數(shù)據(jù)重合度針對(duì)輸入樣本{xt}而言,設(shè)第1組輸入數(shù)據(jù)為{xt,xt-1,…,xt-mx},第2組輸入數(shù)據(jù)為{xt-k,xt-k-1,…,xt-k-mx},把k記作輸入樣本{xt}的平移量,則數(shù)據(jù)重合度定義為

(18)

網(wǎng)絡(luò)訓(xùn)練時(shí),最大迭代次數(shù)設(shè)為1.0×104,目標(biāo)誤差為1.0×10-3,在一定范圍內(nèi)改變網(wǎng)絡(luò)待定系數(shù)(如表中mi,mh和r等),最終求取相對(duì)較優(yōu)的網(wǎng)絡(luò)模型．

針對(duì)上述3組仿真數(shù)據(jù)的建模預(yù)測(cè),各模型結(jié)構(gòu)選取結(jié)果見表2,表中還給出了各模型建模預(yù)測(cè)的各項(xiàng)誤差指標(biāo).通過對(duì)比建模預(yù)測(cè)效果可看出,GNARX模型建模預(yù)測(cè)誤差最小,其預(yù)測(cè)效果見圖1～圖3．

3.2 對(duì)振動(dòng)位移采樣電流數(shù)據(jù)的建模和預(yù)測(cè)

振動(dòng)位移采樣電流數(shù)據(jù)取自文獻(xiàn)[4]．用GNARX模型對(duì)振動(dòng)位移采樣電流數(shù)據(jù)(共150個(gè))建模并預(yù)測(cè)時(shí),首先對(duì)振動(dòng)位移{x1t}和動(dòng)態(tài)切削力{x2t}分別進(jìn)行歸一化處理,然后對(duì){x1t}取前120個(gè)數(shù)據(jù)用于建模,后30個(gè)數(shù)據(jù)用于預(yù)測(cè)．因?yàn)檎駝?dòng)位移與動(dòng)態(tài)切削力存在一定關(guān)系,那么對(duì)應(yīng)的采樣電流數(shù)據(jù)也必然存在某種聯(lián)系,所以本文用動(dòng)態(tài)切削力采樣電流數(shù)據(jù){x2t}作為GNARX模型的輸入．各模型選取范圍見表1,模型選取結(jié)果及建模預(yù)測(cè)誤差列于表3中(α=1,β=1,γ=1)，預(yù)效果見圖4．

通過分析GNARX模型應(yīng)用結(jié)果可知:

1) 在明確系統(tǒng)輸入的情況下,ARX和GNARX模型比AR和GNAR模型更為精確,建模預(yù)測(cè)效果更好．

2) GNARX模型應(yīng)用于仿真數(shù)據(jù)(包括線性和非線性模型數(shù)據(jù)),其建模和預(yù)測(cè)精度明顯高于其他模型,體現(xiàn)了GNARX模型良好的建模預(yù)測(cè)能力．

表2 各模型對(duì)仿真數(shù)據(jù)建模預(yù)測(cè)的誤差對(duì)比

圖1 GNARX模型對(duì)ARX數(shù)據(jù)的預(yù)測(cè)效果

圖2 GNARX模型對(duì)ARMAX數(shù)據(jù)的預(yù)測(cè)效果

圖3 GNARX模型對(duì)NLARX數(shù)據(jù)的預(yù)測(cè)效果

表3 各模型對(duì)振動(dòng)位移采樣電流建模預(yù)測(cè)的誤差對(duì)比

圖4 GNARX模型對(duì)振動(dòng)位移采樣電流數(shù)據(jù)預(yù)測(cè)效果圖

3) GNARX模型應(yīng)用于振動(dòng)位移采樣電流數(shù)據(jù),其建模預(yù)測(cè)效果同樣優(yōu)于其他模型,說明GNARX模型適用于實(shí)際數(shù)據(jù)建模,具有工程應(yīng)用可行性．

4) GNARX模型應(yīng)用于振動(dòng)位移采樣電流數(shù)據(jù)的優(yōu)越性并不明顯,這可能是模型輸入選取不當(dāng)所造成的,說明動(dòng)態(tài)切削力和振動(dòng)位移2組數(shù)據(jù)間并沒有明顯的輸入輸出關(guān)系．

4 結(jié)語

將影響系統(tǒng)輸入的相關(guān)數(shù)據(jù)作為外部輸入,在GNAR模型的基礎(chǔ)上,提出了帶有外部輸入的線性/非線性自回歸模型——GNARX模型,并給出了其一般表達(dá)式及其最小二乘參數(shù)估計(jì)方法．運(yùn)用一種綜合考慮建模誤差、預(yù)測(cè)誤差及模型復(fù)雜度的修正AIC準(zhǔn)則,確定GNARX模型結(jié)構(gòu),該準(zhǔn)則通過改變重要度權(quán)系數(shù)來適應(yīng)不同的現(xiàn)場(chǎng)需求,不斷試驗(yàn)與修改,尋求最優(yōu)模型．將該模型應(yīng)用于仿真和實(shí)際數(shù)據(jù)，結(jié)果表明，GNARX模型精度高,通用性好,無論是對(duì)ARX,ARMAX,NLARX等模型的仿真數(shù)據(jù),還是對(duì)實(shí)際工程數(shù)據(jù),GNARX模型都表現(xiàn)出了很好的建模預(yù)測(cè)能力,效果優(yōu)于傳統(tǒng)的時(shí)序模型及BP網(wǎng)絡(luò)．

由于非線性時(shí)間序列的分析與處理要比線性平穩(wěn)時(shí)序復(fù)雜得多,且GNARX模型結(jié)構(gòu)相對(duì)冗長(zhǎng),在模型結(jié)構(gòu)的確定、適用性檢驗(yàn)等方面尚無統(tǒng)一的方法和成熟的理論研究,因此,該模型的理論研究和實(shí)際應(yīng)用都有待進(jìn)一步探討．

)

[1]Box G E P,Jenkins G M,Reinsel G C.Timeseriesanalysis:forecastingandcontrol[M].4th ed.Hoboken,USA: John Wiley & Sons Inc.,2008.

[2]Flores J J,Graff M,Rodriguez H.Evolutive design of ARMA and ANN models for time series forecasting [J].RenewableEnergy,2012,44: 225-230.

[3]Hansen P R,Huang Z,Shek H H.Realized GARCH: a joint model for returns and realized measures of volatility [J].JournalofAppliedEconometrics,2012,27(6): 877-906.

[4]楊叔子,吳雅,軒建平.時(shí)間序列分析的工程應(yīng)用:下冊(cè)[M].2版.武漢:華中科技大學(xué)出版社,2007.

[5]陳茹雯,黃仁,張志勝,等.基于數(shù)學(xué)模型的視覺測(cè)量系統(tǒng)圖形畸變校正方法[J].機(jī)械工程學(xué)報(bào),2009,45(7):243-248.

Chen Ruwen,Huang Ren,Zhang Zhisheng,et al.Distortion correction method based on mathematic model in machine vision measurement system[J].JournalofMechanicalEngineering,2009,45(7): 243-248.(in Chinese)

[6]Huang Ren,Xu Feiyun,Chen Ruwen.General expression for linear and nonlinear time series models[J].FrontiersofMechanicalEngineeringinChina,2009,4(1): 15-24.

[7]陳茹雯,黃仁,史金飛,等.線性/非線性時(shí)間序列模型一般表達(dá)式及其工程應(yīng)用[J].東南大學(xué)學(xué)報(bào):自然科學(xué)版,2008,38(6):1077-1080.

Chen Ruwen,Huang Ren,Shi Jinfei,et al.General expression for linear and nonlinear time series model and Its engineering application[J].JournalofSoutheastUniversity:NaturalScienceEdition,2008,38(6): 1077-1080.(in Chinese)

[8]Muhannad Z,Yusoff Z M,Rahiman M H F,et al.Modeling of steam distillation pot with ARX model [C]//IEEE8thInternationalColloquiumonSignalProcessingandItsApplications.Melaka,Malaysia,2012: 194-198.

[9]Sanandaji B M,Vincent T L,Wakin M B,et al.Compressive system identification of LTI and LTV ARX models[C]//50thIEEEConferenceonDecisionandControlandEuropeanControlConference.Orlando,FL,USA,2011: 791-798.

[10]Raniman M H F,Taib M N,Salleh Y M.Black box modeling of steam distillation essential oil extraction system using ARMAX structure[C]//InternationalConferenceonIntelligentandAdvancedSystems.Kuala Lumpur,Malaysia,2007: 1059-1062.

[11]Stefanoiu D,Seraficeanu C,Culita J.Identification of MIMO-ARMAX models for glycemia and sodium ions tests through Particle Swarm Optimization[C]//IEEEInternationalConferenceonControlandAutomation.Christchurch,New Zealand,2009: 643-650.

[12]Giardina C R,Chirlian P M.Proof of Weierstrass approximation theorem using band-limited functions[J].ProceedingsoftheIEEE,1973,61(4): 512.