傳遞函數(shù)辨識(20):有限脈沖響應(yīng)滑動平均系統(tǒng)的遞階增廣迭代參數(shù)估計

丁 鋒,劉喜梅

(1.江南大學(xué) 物聯(lián)網(wǎng)工程學(xué)院,江蘇 無錫214122;2.青島科技大學(xué) 自動化與電子工程學(xué)院,山東 青島266061)

輔助模型辨識思想、多新息辨識理論、遞階辨識原理、耦合辨識概念[1-6]是研究新型參數(shù)估計方法的有用工具,連載論文先后研究了信號模型、傳遞函數(shù)的參數(shù)估計[7-11]。最近,針對輸出誤差系統(tǒng),研究了輔助模型多新息隨機梯度算法、輔助模型多新息遞推梯度算法、輔助模型多新息最小二乘遞推算法等[12],以及輔助模型遞階多新息梯度迭代算法、輔助模型遞階多新息最小二乘迭代算法等[13];針對有限脈沖響應(yīng)滑動平均(FIR-MA)系統(tǒng),研究遞階增廣隨機梯度(HESG)算法、遞階多新息增廣隨機梯度(H MI-ESG)算法、遞階增廣遞推梯度(HERG)算法、遞階多新息增廣遞推梯度(H MI-ERG)算法、遞階增廣最小二乘(HELS)遞推算法、遞階多新息增廣最小二乘(H MI-ELS)遞推算法等[14]。本研究針對FIR-MA系統(tǒng),即受控滑動平均(CMA)系統(tǒng),研究遞階增廣梯度迭代(HEGI)算法、遞階多新息增廣梯度迭代(H MI-EGI)算法、遞階增廣最小二乘迭代(HELSI)算法、遞階多新息增廣最小二乘迭代(HMI-ELSI)算法。文獻[15]提出了CMA系統(tǒng)的交互隨機梯度辨識算法。本研究提出的遞階辨識方法可以推廣用于其他線性和非線性隨機系統(tǒng),以及信號模型的參數(shù)辨識[16-26]。

1 有限脈沖響應(yīng)滑動平均系統(tǒng)

考慮下列有限脈沖響應(yīng)滑動平均模型(FIRMA模型)描述的隨機控制系統(tǒng),

其中{u(t)}和{y(t)}分別是系統(tǒng)的輸入和輸出序列,{v(t)}是零均值不可測白噪聲序列,B(z)和D(z)是后移算子z－1的多項式(z－1y(t)＝y(tǒng)(t－1)或zy(t)＝y(tǒng)(t＋1)),定義如下:

定義參數(shù)向量和信息向量:

注意到輸入信息向量φ(t)包含了系統(tǒng)的可測輸入數(shù)據(jù)u(t－i),噪聲向量ψ(t)包含了未知噪聲項v(t－i)。

定義不可測干擾噪聲:

它是一個滑動平均過程。假設(shè)階次n b和n d已知,記n＝n b＋n d。設(shè)t≤0時,各變量的初值為零:u(t)＝0,y(t)＝0,w(t)＝0,v(t)＝0。式(3)和(1)可以寫為

定義虛擬輸出變量

由式(6)可得

式(7)和(4)是FIR-MA系統(tǒng)(1)的遞階辨識模型。它們可以看作兩個子系統(tǒng),一個包含了系統(tǒng)模型的參數(shù)向量b,一個包含了噪聲模型的參數(shù)向量d。由于分解,導(dǎo)致兩個子辨識模型間存在耦合變量(關(guān)聯(lián)項)b和d。因此在推導(dǎo)兩個子辨識算法時,需要根據(jù)遞階辨識原理,協(xié)調(diào)它們間的關(guān)聯(lián)項。

文獻[14]利用遞階辨識原理,研究了有限脈沖響應(yīng)滑動平均系統(tǒng)(1)的遞階增廣隨機梯度算法、遞階多新息增廣隨機梯度算法、遞階增廣遞推梯度算法、遞階多新息增廣遞推梯度算法、遞階增廣最小二乘算法、遞階多新息增廣最小二乘算法等。本工作研究FIR-MA系統(tǒng)的遞階迭代參數(shù)估計方法。

2 遞階增廣梯度迭代辨識方法

2.1 遞階增廣梯度迭代辨識算法

設(shè)L為數(shù)據(jù)長度。對于FIR-MA系統(tǒng)的遞階辨識模型(7)和(4),定義準(zhǔn)則函數(shù):

定義堆積輸出向量Y(L),堆積信息矩陣Φ(L),堆積虛擬輸出向量Y1(L),堆積噪聲向量W(L)和堆積噪聲矩陣Ψ(L)如下:

準(zhǔn)則函數(shù)J1(b)和J2(d)可以等價表示為

令k＝1,2,3,…,是迭代變量和?n d分別是參數(shù)向量b和d的第k次迭代估計,μ≥0和μk≥0是迭代步長(收斂因子),設(shè)λmax[X]為實對稱陣X的最大特征值。使用負(fù)梯度搜索,極小化J1(b)和J2(d),得到梯度迭代關(guān)系:

式(12)和(14)可以看作狀態(tài)為和的離散時間動態(tài)系統(tǒng),為了保證估計向量和^d k收斂,要求矩陣[I nb－μΦT(L)Φ(L)]和[I n d－μkΨT(L)Ψ(L)]的特征值均在單位圓內(nèi),且單位圓上沒有重特征值,即步長μ和μk必須滿足

故收斂因子μ和μk的保守選擇是

由于特征值計算很復(fù)雜,故收斂因子也可簡單更保守地取為

梯度迭代算法(11),(13)和(15)～(18)不可實現(xiàn),因為只有{y(t),φ(t)}是可得到的觀測數(shù)據(jù),即堆積信息矩陣Φ(L)是已知的,而堆積噪聲矩陣Ψ(L)是未知的,故需要構(gòu)造Ψ(L)的估計。

根據(jù)ψ(t)定義式的結(jié)構(gòu),用v(t－i)的估計構(gòu)造信息向量ψ(t)的估計:

從式(6)可得

其中未知的ψ(t),b和d分別用其估計和代替,則v(t)的迭代估計可以通過式(20)計算,

用ψ(t)的估計構(gòu)造Ψ(L)的估計:

基于遞階辨識原理和輔助模型辨識思想,用分別代替式(11),(13)和(15)～(18)右邊未知的Ψ(L),d和b,得到式(22)～(26)的梯度迭代估計和,聯(lián)立式(8)～(9),(21),(2),(19)～(20)就得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的遞階增廣梯度迭代算法(hierarchical extended gradient-based iterative algorith m,HEGI算法):

或

或

HEGI算法(22)～(35)計算量如表1所示(n＝n b＋n d),計算參數(shù)估計和的流程如圖1所示,計算參數(shù)估計的步驟如下。

表1 遞階增廣梯度迭代(HEGI)算法的計算量Table 1 The computational efficiency of the HEGI algorith m

圖1 遞階增廣梯度迭代(HEGI)算法計算流程Fig.1 Flo wchart of the HEGI algorit h m

1)初始化:假設(shè)在t≤0時所有變量取值均為零。令迭代次數(shù)k＝1,給定數(shù)據(jù)長度L和參數(shù)估計精度ε,置參數(shù)估計初值。置初值

2)采集觀測數(shù)據(jù)u(t)和y(t),用式(32)構(gòu)造信息向量φ(t),t＝1,2,…,L。

3)用式(29)構(gòu)造堆積輸出向量Y(L),用式(30)構(gòu)造堆積信息矩陣Φ(L),根據(jù)式(23)或(24)確定一個盡量大步長μ。

5)根據(jù)式(26)或(27)確定一個盡量大步長μk,用式(22)刷新參數(shù)估計向量,用式(25)刷新參數(shù)估計向量。

2.2 加權(quán)遞階增廣梯度迭代辨識算法

設(shè)L為數(shù)據(jù)長度,M(L)＝MT(L)∈?L×L為非負(fù)定加權(quán)矩陣(當(dāng)然可以使用兩個不同的加權(quán)矩陣)。對于FIR-MA系統(tǒng)的遞階辨識模型(7)和(4),將準(zhǔn)則函數(shù)J1(b)和J2(d)修改為加權(quán)準(zhǔn)則函數(shù):

仿照HEGI算法(22)～(35)的推導(dǎo),極小化J′1(b)和J′2(d),可以得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的加權(quán)遞階增廣梯度迭代算法(weighteg hierarchical extended gradient-based iterative algorith m,W-HEGI算法):

或

或

設(shè)η為遺忘因子,在W-HEGI辨識算法(36)～(49)中,取加權(quán)矩陣為

就得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的遺忘因子遞階增廣梯度迭代算法(f or getting factor HEGI al gorith m,FF-HEGI算法)。

3 遞階增廣最小二乘迭代方法

3.1 遞階增廣最小二乘迭代辨識算法

令準(zhǔn)則函數(shù)J1(b)和J2(d)分別對參數(shù)向量b和d的偏導(dǎo)數(shù)為零:

假設(shè)ΦT(L)Φ(L)和ΨT(L)Ψ(L)是可逆矩陣,那么從上兩式可以求得參數(shù)向量b和d的最小二乘估計:

因為在式(50)～(51)中,只有堆積輸出向量Y(L)和堆積信息矩陣Φ(L)是已知的,而堆積噪聲矩陣Ψ(L)是未知的,右邊還包含了未知參數(shù)向量d和b,故無法由上兩式求得參數(shù)估計。解決方案是基于遞階辨識原理和輔助模型辨識思想,用式(21)中的,以及前一次迭代參數(shù)估計和分別代替式(50)～(51)右邊未知的Ψ(L),d和b,得到式(52)～(53)的最小二乘迭代參數(shù)估計和,聯(lián)立式(28)～(35),就得到估計FIR-MA系統(tǒng)參數(shù)向量b和d的遞階增廣最小二乘迭代算法(hierarchical extended least squares-based iterative al gorit h m,HELSI算法):

HLSI算法(52)～(61)計算量如表2所示(n＝n b＋n d),計算參數(shù)估計和的流程如圖2所示,計算參數(shù)估計的步驟如下。

表2 遞階增廣最小二乘迭代(HELSI)算法的計算量Table 2 The co mputational efficiency of the HELSI algorith m

圖2 遞階增廣最小二乘迭代(HELSI)算法的計算流程Fig.2 Flowchart of the HELSI algorit h m

1)初始化:假設(shè)在t≤0時所有變量取值均為零。令迭代次數(shù)k＝1,給定數(shù)據(jù)長度L和參數(shù)估計精度ε,置參數(shù)估計初值和為隨機向量,或和置為隨機數(shù),t＝1,2,…,L。隨機數(shù)初值是為保證式(52)～(53)中的矩陣可逆。

2)采集觀測數(shù)據(jù)u(t)和y(t)。用式(58)構(gòu)造信息向量φ(t),t＝1,2,…,L。

3)用式(55)構(gòu)造堆積輸出向量Y(L),用式(56)構(gòu)造堆積信息矩陣Φ(L)。

5)用式(52)刷新參數(shù)估計向量,用式(53)刷新參數(shù)估計向量。

3.2 加權(quán)遞階增廣最小二乘迭代辨識算法

仿照HELSI算法(52)～(61)的推導(dǎo),極小化J′1(b)和J′2(d),可以得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的加權(quán)遞階增廣最小二乘迭代算法(weighted hierarchical extended least squaresbased iterative algorith m,W-HELSI算法):

設(shè)η為遺忘因子,在W-HELSI辨識算法(62)～(71)中,取加權(quán)矩陣為

就得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的遺忘因子遞階增廣最小二乘迭代算法(f or getting factor HELSI algorith m,FF-HELSI算法)。

4 遞階多新息增廣梯度迭代辨識方法

4.1 遞階多新息增廣梯度迭代辨識算法

設(shè)p為數(shù)據(jù)窗長度(p?n b＋n d)。根據(jù)式(7)和(4),定義堆積輸出向量Y(p,t)和堆積信息矩陣Φ(p,t),堆積虛擬輸出向量Y1(p,t),堆積噪聲向量W(p,t):和堆積噪聲矩陣Ψ(p,t)如下:

根據(jù)辨識模型(7)和(4),定義關(guān)于參數(shù)向量b和d的二次準(zhǔn)則函數(shù):

令k＝1,2,3,…,是迭代變量,μ(t)≥0和μk(t)≥0是迭代步長(收斂因子)和?n d分別是參數(shù)向量b和d在當(dāng)前時刻t的第k次迭代估計。使用負(fù)梯度搜索,極小化準(zhǔn)則函數(shù)J3(b)和J4(d),得到梯度迭代關(guān)系:

故收斂因子μ(t)和μk(t)的保守選擇是

由于特征值計算很復(fù)雜,故收斂因子也可簡單取為

梯度迭代算法(75),(77)和(79)～(82)不可實現(xiàn),因為右邊包含了未知的堆積噪聲矩陣Ψ(p,t),以及未知參數(shù)向量d和b。為此根據(jù)ψ(t)定義式的結(jié)構(gòu),用υ(t－i)的估計^υk－1(t－i)構(gòu)造信息向量ψ(t)的估計:

從式(6)可得

其中未知的ψ(j),b和d分別用其估計和代替,則υ(t)的估計可以通過式(84)計算,

或

或

式(92)可以用下列2式代替:

遞階多新息增廣梯度迭代算法又稱移動數(shù)據(jù)窗遞階增廣梯度迭代算法(MDW-HEGI算法)。H MIEGI算法(86)～(99)計算參數(shù)估計(t)和(t)的流程如圖3所示,計算參數(shù)估計的步驟如下。

圖3 遞階多新息增廣梯度迭代(HMI-EGI)算法的計算流程Fig.3 Flowchart of the HMI-EGI algorith m

1)初始化:假設(shè)在t≤0時所有變量取值均為零。令t＝1,給定移動數(shù)據(jù)窗長度p(通常取p?n),參數(shù)估計精度ε,最大迭代次數(shù)kmax。置參數(shù)估計初值

2)令k＝1,置n d。采集觀測數(shù)據(jù)u(t)和y(t)。用式(96)構(gòu)造信息向量φ(t),用式(93)構(gòu)造堆積輸出向量Y(p,t),用式(94)構(gòu)造堆積信息矩陣Φ(p,t)。根據(jù)式(87)或(88)確定一個盡量大步長μ(t)。

4)根據(jù)式(90)或(91)確定一個盡量大步長μk(t),用式(86)刷新參數(shù)估計向量,用式(89)刷新參數(shù)估計向量

6)如果k＜kmax,k就增加1,轉(zhuǎn)到步驟3);否則執(zhí)行下一步。

4.2 加權(quán)遞階多新息增廣梯度迭代辨識算法

設(shè)p為新息長度為非負(fù)定加權(quán)矩陣(當(dāng)然可以使用兩個不同的加權(quán)矩陣)。對于FIR-MA系統(tǒng)的遞階辨識模型(7)和(4),將準(zhǔn)則函數(shù)J3(b)和J4(d)修改為加權(quán)準(zhǔn)則函數(shù):

仿照H MI-EGI算法(86)～(99)的推導(dǎo),極小化J′3(b)和J′4(d),便得到估計FIR-MA系統(tǒng)參數(shù)向量b和d的加權(quán)遞階多新息增廣梯度迭代算法(weighted H MI-EGI algorit h m,W-HMI-EGI算法):

或

5 遞階多新息增廣最小二乘迭代方法

5.1 遞階多新息增廣最小二乘迭代辨識算法

由于最小二乘迭代算法涉及矩陣逆,所以至少要求p≥n b＋n d,以及信息向量是持續(xù)激勵的。令準(zhǔn)則函數(shù)J3(b)和J4(d)分別對參數(shù)向量b和d的偏導(dǎo)數(shù)為零:

假設(shè)ΨT(p,t)Ψ(p,t)是可逆矩陣,那么從兩上式可以求得參數(shù)向量b和d的最小二乘估計:

遞階多新息增廣最小二乘迭代算法又稱移動數(shù)據(jù)窗遞階增廣最小二乘迭代算法(MDW-HELSI algorit h m,MDW-HELSI算法)。H MI-ELSI算法(118)～(127)計算參數(shù)估計(t)和(t)的流程如圖4所示,計算參數(shù)估計的步驟如下。

1)初始化:假設(shè)在t≤0,時所有變量取值均為零。令t＝l,給定移動數(shù)據(jù)窗長度p(通常取p?n),參數(shù)估計精度ε,最大迭代次數(shù)kmax,置參數(shù)估計初值,采集觀測數(shù)據(jù){u(j),y(j):j＝1,2,…,l－1},這里l≥p為一個正整數(shù)(請思考l的意義)。

2)令k＝1。置為隨機數(shù),j＝1,2,…,L。隨機數(shù)初值是為保證式(119)中的矩陣可逆。

3)采集觀測數(shù)據(jù)u(t)和y(t)。用式(124)構(gòu)造信息向量φ(t),用式(121)構(gòu)造堆積輸出向量Y(p,t),用式(122)構(gòu)造堆積信息矩陣Φ(p,t)。

圖4 遞階多新息增廣最小二乘迭代(HMI-ELSI)算法的計算流程Fig.4 Flowchart of the HMI-ELSI algorith m

4)用式(125)構(gòu)造(t),用式(123)構(gòu)造堆積噪聲矩陣(p,t)。

7)如果k＜kmax,k就增加1,轉(zhuǎn)到步驟4);否則執(zhí)行下一步。

5.2 加權(quán)遞階多新息增廣最小二乘迭代辨識算法

仿照H MI-ELSI算法(118)～(127)的推導(dǎo),極小化J′3(b)和J′4(d),便得到估計FIR-MA系統(tǒng)參數(shù)向量b和d的加權(quán)遞階多新息增廣最小二乘迭代算法(weighted H MI-ELSI algorith m,W-HMIELSI算法):

6 結(jié) 語

針對最簡單的有色噪聲干擾的有限脈沖響應(yīng)滑動平均模型,基于遞階辨識原理,通過極小化兩個最小二乘準(zhǔn)則函數(shù),研究了遞階(多新息)增廣梯度迭代辨識算法和遞階(多新息)增廣最小二乘迭代辨識算法等。提出的遞階迭代方法可以推廣研究大規(guī)模線性隨機系統(tǒng)和復(fù)雜非線性隨機系統(tǒng)的迭代辨識,其干擾可以是有色噪聲。