一种线性/非线性自回归模型及其在建模和预测中的应用

马家欣 许飞云 黄 仁

(东南大学机械工程学院,南京 211189)

时间序列分析是现代系统辨识的重要方法之一,它在缺乏明确或者完全的系统输入与输出因果关系的情况下,将白噪声看作系统输入,从数理统计的角度揭示因果关系,最终求得系统的等价模型．该方法在各个领域都得到了广泛的应用[1-4]．

ARMA模型(包括AR模型和MA模型)是时序方法中最基本的时序模型,然而,它是在线性回归模型的基础上引申发展起来的,很难模拟工程实际中的非线性现象．而经典的非线性时间序列模型,包括门限自回归模型、双线性模型、指数自回归模型等,一般都是在特定的工程背景下提出的,通用性不够．因此亟须提出一种适用于线性/非线性系统的时序模型．GNAR模型的提出和应用,在形式上统一了不同背景下提出的非线性时序模型的表达式,并有机结合了线性模型[5-7]．另一方面,尽可能地利用更多的已知信息,获取更多的系统特性,可进一步提高建模精度．因此,在经典时序模型的基础上,把影响系统输出的已知相关数据作为系统输入,可得到带有外部输入的时序模型,如ARX,ARMAX等模型．这些模型不仅在理论上得到了验证,在工程实际中也得到了广泛应用[8-11]．本文针对系统输入部分已知的特点,在GNAR模型基础上提出了GNARX模型,给出了模型的参数估计和结构选取方法,并将该模型应用于仿真和实际数据,效果优于其他模型．

1 GNARX模型表达式

GNAR模型建模时,认为系统输入是零均值的白噪声.如果系统已知单个输入,记作序列{ut},系统的输出序列记作{wt},白噪声记作{at},那么GNAR模型就变形为单输入的GNARX模型．

对于单输入系统,设t时刻输出值为wt,表示为函数f的表达式,即

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

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

(1)

式中,f(·)为任意一般线性/非线性函数;wt-i为t-i时刻的系统输出观测值,i=1,2,…;ut-τu-i为t-τu-i时刻的系统输入,i=1,2,…;at-i为t-i时刻的白噪声,i=1,2,…;τu为系统输入ut的延迟．

对于一般线性平稳系统,式(1)中函数f表示ARMAX模型;而当系统表现出非线性、非平稳特征时,由Weierstrass逼近原理[12]知,可用多项式逼近式(1)中的f函数:

λ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,…是模型参数．

假设式(2)描述的系统具有零初始状态(即wt=0,t≤0),且系统在初始时刻之前没有输入(即ut=0,t≤0)．实际建模中,模型阶次p取有限值,令xt,i表示式(2)中的各i阶项,xt,i,1(i=1,2,…,p)表示i阶项xt,i中的1次项,即

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

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

(3)

则系统输出wt可写成如下形式:

(4)

式中,θ(i1),θ(i1,i2),…为模型参数;xt,i,1(j)为向量xt,i,1中的第j个元素;nw,j(j=1,2,…,p)为输出{wt}各阶项的记忆步长;nu,j(j=1,2,…,p)为输入{ut}各阶项的记忆步长．该模型简记为GNARX(p;τu;nw,1,nw,2,…,nw,p;nu,1,nu,2,…,nu,p)．

当系统具有双输入{ut},{vt}时,把式(4)推广至双输入的GNARX模型,简记为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}各阶项的记忆步长;τv为系统输入{vt}的延迟．

同理,式(6)可推广到多输入系统,不再赘述．

2 GNARX模型的参数估计与结构选取

GNARX模型的参数估计用最小二乘法很容易求得.以式(4)所示的单输入GNARX模型为例,在时刻t,p阶项记作xt,p,是一维向量形式．求取过程中用xt,p,i(i=1,2,…,p)表示p阶项xt,p中的i次项,也是一维向量,其中xt,p,i(j)是标量,表示向量xt,p,i中的第j个元素．

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时刻到t+k时刻的方程式写成

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

θ参数最小二乘法估计为

(9)

其中

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

根据式(9)即可对单输入GNARX模型进行参数估计．同理,可获得双输入、多输入模型的参数估计方法．

GNARX模型结构的辨识,包括模型阶次p的确定和记忆长度的确定,这直接影响着模型的建模性能和预测性能及模型的复杂程度,同时也关系到计算速度、预报实时性等问题．

对于GNARX模型,定义单输入模型i阶记忆步长si=nw,i+nu,i(i=1,2,…,p),双输入模型i阶记忆步长si=nw,i+nu,i+nv,i(i=1,2,…,p),多输入模型i阶记忆步长可依此类推得到．则i阶参数量ri为

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

则GNARX模型参数数量为

R=r1+r2+…+rp

修正AIC准则函数[6]如下所示:

(10)

模型结构的最终确定,需要综合考虑模型的建模能力、预测能力和模型复杂程度等因素．在具体建模过程中,可根据现场要求的不同,调整式(10)中的重要度权系数,来满足实际建模需求．如模型用于滤波,则可增加建模重要度权系数α;若模型以预报为主,则可以增加预测重要度权系数β;而当现场侧重于建模实时性时,则可以增加模型复杂重要度权系数γ．

3 GNARX模型的建模和预测应用

GNARX模型适用范围广,对线性和非线性数据有着良好的建模预测能力．下文用GNARX模型对仿真数据和实际数据进行建模和预测,并与其他模型预测效果进行对比．

为比较建模预测效果,定义如下衡量指标：

建模均方误差

(11)

预测均方误差

(12)

建模平均绝对百分比误差

(13)

预测平均绝对百分比误差

(14)

3.1 对仿真数据的建模和预测

仿真数据模型具体形式如下:

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为白噪声,服从正态分布N(0,0.1);ut-i(i=1,2,3,…)为系统输入，服从平均分布U(0,1);xt-i(i=0,1,2,…)为系统输出;F(·)取Matlab中S形网络非线性估计函数,单元数10,S形函数表达式为f(z)=1/(1+exp(-z))．

仿真数据选用前,先对{xt}序列进行归一化处理,再舍去前50个过渡区数据,取中间250个数据用于建模,后100个数据用于预测．分别用AR模型、GNAR模型、ARX模型、BP网络和GNARX模型对上述各仿真数据进行建模预测．其中时序模型是在合适的模型结构选取范围内选取的,并用式(10)所示修正AIC判断,重要度权系数尝试选取α=1,β=1,γ=1,模型具体选取范围见表1．

表1 针对不同数据建模时各模型结构的选取范围

表1中,BP神经网络的结构包括输入层、隐层和输出层3层,其中输出层节点数mo为1,隐层节点数待定．对于网络的输入,可分为2种情况:

1) 仅以建模预测数据xt作为网络的输入,输入样本数据重合度记作r,输入层节点数记作mi,网络简记为BP1(mi;mh;mo;r)．

2) 数据xt,ut同时作为网络的输入,xt输入维数为mx,ut输入维数为mu,延迟为τu,输入数据xt的重合度记作r,此时输入层节点数为mx+mu,网络简记为BP2(mx;mu,τu;mh;mo;r)．

其中数据重合度针对输入样本{xt}而言,设第1组输入数据为{xt,xt-1,…,xt-mx},第2组输入数据为{xt-k,xt-k-1,…,xt-k-mx},把k记作输入样本{xt}的平移量,则数据重合度定义为

(18)

网络训练时,最大迭代次数设为1.0×104,目标误差为1.0×10-3,在一定范围内改变网络待定系数(如表中mi,mh和r等),最终求取相对较优的网络模型．

针对上述3组仿真数据的建模预测,各模型结构选取结果见表2,表中还给出了各模型建模预测的各项误差指标.通过对比建模预测效果可看出,GNARX模型建模预测误差最小,其预测效果见图1～图3．

3.2 对振动位移采样电流数据的建模和预测

振动位移采样电流数据取自文献[4]．用GNARX模型对振动位移采样电流数据(共150个)建模并预测时,首先对振动位移{x1t}和动态切削力{x2t}分别进行归一化处理,然后对{x1t}取前120个数据用于建模,后30个数据用于预测．因为振动位移与动态切削力存在一定关系,那么对应的采样电流数据也必然存在某种联系,所以本文用动态切削力采样电流数据{x2t}作为GNARX模型的输入．各模型选取范围见表1,模型选取结果及建模预测误差列于表3中(α=1,β=1,γ=1)，预效果见图4．

通过分析GNARX模型应用结果可知:

1) 在明确系统输入的情况下,ARX和GNARX模型比AR和GNAR模型更为精确,建模预测效果更好．

2) GNARX模型应用于仿真数据(包括线性和非线性模型数据),其建模和预测精度明显高于其他模型,体现了GNARX模型良好的建模预测能力．

表2 各模型对仿真数据建模预测的误差对比

图1 GNARX模型对ARX数据的预测效果

图2 GNARX模型对ARMAX数据的预测效果

图3 GNARX模型对NLARX数据的预测效果

表3 各模型对振动位移采样电流建模预测的误差对比

图4 GNARX模型对振动位移采样电流数据预测效果图

3) GNARX模型应用于振动位移采样电流数据,其建模预测效果同样优于其他模型,说明GNARX模型适用于实际数据建模,具有工程应用可行性．

4) GNARX模型应用于振动位移采样电流数据的优越性并不明显,这可能是模型输入选取不当所造成的,说明动态切削力和振动位移2组数据间并没有明显的输入输出关系．

4 结语

将影响系统输入的相关数据作为外部输入,在GNAR模型的基础上,提出了带有外部输入的线性/非线性自回归模型——GNARX模型,并给出了其一般表达式及其最小二乘参数估计方法．运用一种综合考虑建模误差、预测误差及模型复杂度的修正AIC准则,确定GNARX模型结构,该准则通过改变重要度权系数来适应不同的现场需求,不断试验与修改,寻求最优模型．将该模型应用于仿真和实际数据，结果表明，GNARX模型精度高,通用性好,无论是对ARX,ARMAX,NLARX等模型的仿真数据,还是对实际工程数据,GNARX模型都表现出了很好的建模预测能力,效果优于传统的时序模型及BP网络．

由于非线性时间序列的分析与处理要比线性平稳时序复杂得多,且GNARX模型结构相对冗长,在模型结构的确定、适用性检验等方面尚无统一的方法和成熟的理论研究,因此,该模型的理论研究和实际应用都有待进一步探讨．

)

[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]杨叔子,吴雅,轩建平.时间序列分析的工程应用:下册[M].2版.武汉:华中科技大学出版社,2007.

[5]陈茹雯,黄仁,张志胜,等.基于数学模型的视觉测量系统图形畸变校正方法[J].机械工程学报,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]陈茹雯,黄仁,史金飞,等.线性/非线性时间序列模型一般表达式及其工程应用[J].东南大学学报:自然科学版,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.