丁 鋒,徐 玲,籍 艷,劉喜梅
(1.江南大學 物聯(lián)網工程學院,江蘇 無錫 214122;2.青島科技大學 自動化與電子工程學院,山東 青島 266061)
方程誤差自回歸滑動平均模型(equation-error ARMA model,EEARMA模型)是方程誤差類模型的最一般形式。它的特例有方程誤差滑動平均模型(equation-err or moving average model,EEMA模型),方程誤差自回歸模型(equation-err or autoregressive model,EEAR模型)。EEARMA模型也稱為受控自回歸自回歸滑動平均模型(controlled autoregressive ARMA model,CARARMA模型),即受控自回歸ARMA模型或受控ARARMA模型,或ARARMAX模型(即帶外加輸入的ARARMA模型)。
系統(tǒng)辨識是研究建立系統(tǒng)數(shù)學模型的理論與方法,高精度的數(shù)學模型是實現(xiàn)高精度控制的基礎。最近出版的系統(tǒng)辨識學術專著[1-7]和連載論文[8-32]中介紹了許多系統(tǒng)辨識方法。最近的連載論文研究了輸出誤差模型的輔助模型遞階迭代參數(shù)估計方法[25]、有限脈沖響應滑動平均系統(tǒng)的遞階增廣迭代參數(shù)估計方法[27]、線性回歸系統(tǒng)的遞階迭代參數(shù)估計方法[29]、線性回歸系統(tǒng)的變間隔遞階迭代參數(shù)估計方法[31]等?!断到y(tǒng)辨識:迭代搜索原理與辨識方法》第2章研究了方程誤差滑動平均系統(tǒng)的增廣迭代辨識方法、方程誤差自回歸系統(tǒng)的廣義迭代辨識方法、方程誤差自回歸滑動平均系統(tǒng)的廣義增廣迭代辨識方法[5]。本工作基于遞階辨識原理,研究方程誤差自回歸滑動平均系統(tǒng)的遞階廣義增廣迭代辨識方法。有關遞階辨識原理與方法可參見專著[7]。
考慮下列方程誤差自回歸滑動平均模型(EEARMA模型)描述的隨機控制系統(tǒng),
其中{u(t)}和{y(t)}分別是系統(tǒng)的輸入和輸出序列,{v(t)}是零均值不可測白噪聲序列,A(z),B(z),C(z)和D(z)是后移算子z-1的多項式(z-1y(t)=y(tǒng)(t-1)或zy(t)=y(tǒng)(t+1)),定義如下,
定義干擾噪聲,
這里w(t)是一個ARMA噪聲,是有色噪聲。假設階次na,nb,nc和nd已知,記n=na+nb+nc+nd。設t≤0時,各變量的初值為零:u(t)=0,y(t)=0,w(t)=0,v(t)=0。
定義參數(shù)向量:
和信息向量:
注意到輸入信息向量φu(t)包含了系統(tǒng)輸入u(t-i),輸出信息向量φy(t)包含了系統(tǒng)輸出y(t-i),它們是已知的;噪聲向量ψw(t)包含了不可測相關噪聲項w(t-i),噪聲向量ψv(t)包含了不可測白噪聲項v(t-i),它們是未知的。
由式(2)可得
定義虛擬輸出變量:
由式(7)可得4個子辨識模型:
式(12)~(15)是EEARMA系統(tǒng)(1)的基于模型分解得到的4個子辨識模型,即遞階辨識模型。4個子辨識模型分別包含了系統(tǒng)模型的參數(shù)向量a和b,以及噪聲模型的參數(shù)向量c和d。
由于分解,導致4個子辨識模型間存在耦合變量(關聯(lián)項)b,c,d和a。因此在推導4個子辨識算法時,需要根據遞階辨識原理,協(xié)調它們之間的關聯(lián)項?;诜纸獾谋孀R方法也稱為遞階辨識方法。
基于遞階辨識模型(12)~(15),定義4個準則函數(shù):
令k=1,2,3,…是迭代變量和分別是參數(shù)向量a,b,c和d的第k次迭代估計,1/r1(L),1/r2(L),1/r3,k(L)和1/r4,k(L)為非負迭代步長(收斂因子),使用負梯度搜索,極小化J1(a),J2(b),J3(c)和J4(d),得到梯度迭代關系:
式(16)~(19)右邊虛擬輸出yi(t),噪聲向量ψw(t)和ψv(t)是未知的,所以無法計算出參數(shù)估計和。解決的方案是,先獲得yi(t),ψw(t)和ψv(t)的估計和,將辨識算法中未知量用他們的估計代替。定義式(8)~(11)中yi(t)涉及到未知參數(shù)向量b,c,d和a,以及未知噪聲向量ψw(t)和ψv(t),為此用噪聲w(t-i)和v(t-i)的估計和定義噪聲向量ψw(t)和ψv(t)的估計:
根據yi(t)的定義式(8)~(11),用參數(shù)向量b,c,d和a的前 一次迭代 估計和,以及噪聲向量ψw(t)和ψv(t)的估計和定義yi(t)的估計,即
式(16)~(19)中未知的yi(t),ψw(t)和ψv(t)用其估計和代替,得到式(28)~(31),聯(lián)立式(20)~(27)和(3)~(4),便得到估計參數(shù)向量a,b,c和d的遞階廣義增廣梯度迭代算法(hierarchical generalized extended gradientbased iterative al gorit h m,HGEGI算法):
定理1對HGEGI算法(28)~(42),如果r1(L),r2(L),r3,k(L)和r4,k(L)滿足
也可采用下列遞推算式計算:
1)初始化:給定數(shù)據長度L和參數(shù)估計精度ε,置參數(shù)估計,。
2)令k=1,采集觀測數(shù)據u(t)和y(t),置,…,L。用式(36)~(37)構造輸出信息向量φy(t)和輸入信息向量φu(t),t=1,2,…,L。根據式(43)和(44)計算r1(L)和r2(L)。
設t=l為用作辨識的數(shù)據起點,L為數(shù)據長度?;谶f階辨識模型(12)~(15),定義4個準則函數(shù):
定義協(xié)方差陣Pa(L),Pb(L),Pc(L)和Pd(L)如下:
令k=1,2,3,…是迭代變量,和分別是參數(shù)向量a,b,c和d在時刻t的第k次迭代估計。令準則函數(shù)J5(a),J6(b),J7(c)和J8(d)在點,和的梯度為零,得到
假設φy(t),φu(t),ψw(t)和ψv(t)是持續(xù)激勵的,從以上4式可以求得
使用式(8)~(11),可得
讀者可以討論引入加權因子和(或)遺忘因子的加權遞階廣義增廣最小二乘迭代算法(W-HGELSI算法)和遺忘因子遞階廣義增廣最小二乘迭代算法(FF-HGELSI算法)。
1)初始化:給定數(shù)據長度L(通常取L?n)和參數(shù)估計精度ε,置參數(shù)估計初值,。
2)令k=1。采集輸入輸出數(shù)據u(t)和y(t),置,…,l+L-1。用式(79)~(80)構造輸出信息向量φy(t)和輸入信息向量φu(t),t=l,l+1,…,l+L-1。
5)用式(68),(71),(74)和(77)計算增益向量La(t),Lb(t),Lc(t)和Ld(t),用式(69),(72),(75)和(78)計算協(xié)方差陣Pa(t),Pb(t),Pc(t)和Pd(t)。
6)如果t<L,t就增加1,轉到步驟5);否則用式(67)、(70)、(73)和(76)刷新參數(shù)估計向量,和。用式(85)構成系統(tǒng)的參數(shù)估計向量。
7)用式(83)和(84)計算噪聲的估計^wk(t)和,t=1,2,…,l+L-1。
基于遞階辨識模型(12)~(15),定義4個準則函數(shù):
令k=1,2,3,…是迭代變量,,和分別是參數(shù)向量a,b,c和d在時刻t的第k次迭代估計,μ1(t),μ2(t),μ3(t)和μ4(t)為非負迭代步長(收斂因子),使用負梯度搜索,極小化J9(a),J10(b),J11(c)和J12(d),得到梯度迭代關系:
因為虛擬輸出yi(t),噪聲向量ψw(t)和ψv(t)是未知的,所以以上幾式無法計算出參數(shù)估計,和。因此首先要獲得yi(t),ψw(t)和ψv(t)的估計和)。
由yi(t)的定義式(8)~(11)可知,它們涉及到未知參數(shù)向量b,c,d和a,以及未知噪聲向量ψw(t)和ψv(t)。為此用噪聲w(t-i)和v(t-i)的估計和定義噪聲向量ψw(t)和ψv(t)的估計:
由yi(t)的定義式(8)~(11)可得
據此,用參數(shù)向量b,c,d和a的估計,和,以及噪聲向量ψw(j)和ψv(j)的估計和定義估算yi(j)的輔助模型:
由式(6)有
據此,用參數(shù)向量a和b的估計和,及信息向量φy(j)和φu(j)定義噪聲w(j)的估計
由式(5)有
據此,用w(j)的估計^wk(j),參數(shù)向量c和d的估計和,以及信息向ψw(j)和ψv(j)的估計和定義噪聲v(j)的估計
式(86)~(89)中未知的yi(j),ψw(j)和ψv(j)用其估計和代替,得到式(98)~(101),聯(lián)立式(90)~(97)和(3)~(4),便得到估計參數(shù)向量a,b,c和d的遞階多新息廣義增廣梯度迭代算法(hierarchical multi-innovation generalized extended gradient-based iterative algorith m,HMI-GEGI算法):
定義新息
定理2對H MI-GEGI算法(98)~(112),如果收斂因子滿足
為簡化計算,步長也可取為
遞階多新息廣義增廣梯度迭代算法又稱移動數(shù)據窗遞階廣義增廣梯度迭代算法(MDW-HGEGI算法)。H MI-GEGI算法(98)~(112)和(117)~(120)計算參數(shù)估計的步驟如下。
1)初始化:令t=1,給定移動數(shù)據窗長度p(通常取p?n),參數(shù)估計精度ε,最大迭代次數(shù)kmax。置參數(shù)估計初值,。
2)令k=1,置。用式(106)~(107)構造輸出信息向量φy(t-j)和輸入信息向量φu(t-j),用式(108)~(109)構造噪聲向量j)和。
3)采集觀測數(shù)據u(t)和y(t)。用式(106)構造輸出信息向量φy(t),用式(107)構造輸入信息向量φu(t),根據式(117)和(118)計算步長μ1(t)和μ2(t)。
7)如果k<kmax,k就增加1,轉到步驟(4);否則執(zhí)行下一步。
令k=1,2,3,…是迭代變量,和分別是參數(shù)向量a,b,c和d在時刻t的第k次迭代估計。令準則函數(shù)J9(a),J10(b),J11(c)和J12(d)在點a=和的梯度為零,得到
假設φy(t),φu(t),ψw(t)和ψv(t)是持續(xù)激勵的,從以上4式可以求得
因為以上4式右邊包含了未知參數(shù)向量b,c,d和a,以及未知噪聲向量ψw(j)和ψv(j),所以無法計算出參數(shù)估計和。解決的辦法是,未知噪聲向量ψw(j)和ψv(j)用其估計和)代替,未知參數(shù)向量b,c,d和a用其前一次迭代估計和代替,代替后得到式(129)~(132),聯(lián)立式(106)~(112),便得到估計參數(shù)向量a,b,c和d的遞階多新息廣義增廣最小二乘迭代算法(hierarchical multi-innovation generalized extended least squares-based iterative al gorit h m,HMI-GELSI算法):
遞階多新息廣義增廣最小二乘迭代算法又稱移動數(shù)據窗遞階廣義增廣最小二乘迭代算法(MDWHGELSI算法)。H MI-GELSI算法(129)~(139)計算參數(shù)估計和)步驟如下。
1)初始化:給定移動數(shù)據窗長度p(通常取p?n),令t=l(l≥p),給定參數(shù)估計精度ε,最大迭代次數(shù)kmax。置參數(shù)估計初值。采集輸入輸出數(shù)據u(j)和y(j),j=1,2,…,l。置噪聲初值=隨機數(shù),)=隨機數(shù),j=1,2,…,l。令nm=p+max[na,nb,nc,nd]。
2)令k=1,用式(133)~(134)構造輸出信息向量φy(t-j)和輸入信息向量φu(t-j),用式(135)~(136)構造噪聲向量和,j=1,2,…,nm。
3)采集觀測數(shù)據u(t)和y(t)。用式(133)構造輸出信息向量φy(t),用式(134)構造輸入信息向量φu(t)。
6)如果k<kmax,k就增加1,轉到步驟4);否則執(zhí)行下一步。
遞階多新息廣義增廣最小二乘迭代算法(129)~(139)可以等價表示為
H MI-GELSI算法(140)~(155)計算參數(shù)估計和的步驟如下。
1)初始化:給定移動數(shù)據窗長度p(通常取p?n),令t=l(l≥p),給定參數(shù)估計精度ε,最大迭代次數(shù)kmax。置參數(shù)估計初值,。采集輸入輸出數(shù)據u(j)和y(j),置^w0(j)=隨機數(shù),=隨機數(shù),j=1,2,…,l。令nm=p+max[na,nb,nc,nd]。
2)令k=1,用式(149)~(150)構造輸出信息向量φy(t-j)和輸入信息向量φu(t-j),用式(151)~(152)構造噪聲向量和,2,…,nm,
3)采集觀測數(shù)據u(t)和y(t)。用式(149)構造輸出信息向量φy(t),用式(150)構造輸入信息向量φu(t)。用式(144)構造堆積輸出向量Y(p,t),用式(145)和(146)構造堆積輸出信息矩陣Φy(p,t)和堆積輸入信息矩陣Φu(p,t)。
7)如果k<kmax,k就增加1,轉到步驟(4);否則執(zhí)行下一步。
針對方程誤差自回歸滑動平均(EEARMA)系統(tǒng),研究和提出了遞階(多新息)廣義增廣梯度迭代辨識算法、遞階(多新息)廣義增廣最小二乘迭代辨識算法。盡管這些遞階迭代辨識方法是針對有色噪聲干擾下的標量方程誤差隨機系統(tǒng)提出的,但是其思想可以推廣到有色噪聲干擾下的線性和非線性多變量輸出誤差隨機系統(tǒng)中[33-35]。