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

丁 鋒,劉喜梅

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

在連載論文中，使用多新息辨識理論、遞階辨識原理、耦合辨識概念[1-6]，先后研究了信號模型、傳遞函數(shù)的參數(shù)估計[7-11]。最近，針對輸出誤差模型描述的動態(tài)隨機(jī)系統(tǒng):

研究了輔助模型(多新息)隨機(jī)梯度算法、輔助模型(多新息)遞推梯度算法、輔助模型(多新息)最小二乘算法等[12]，以及輔助模型遞階(多新息)梯度迭代算法、輔助模型遞階(多新息)最小二乘迭代算法等[13]。

有限脈沖響應(yīng)滑動平均模型(finite impulse response moving average model，F(xiàn)IR-MA模型)也稱為受控滑動平均模型(controlled moving average model，CMA模型)，它是一類特殊的方程誤差滑動平均模型?；瑒悠骄肼暩蓴_下的系統(tǒng)辨識方法稱為增廣辨識方法。

本研究針對有限脈沖響應(yīng)滑動平均(FIR-MA)系統(tǒng)，通過辨識模型分解，研究遞階增廣隨機(jī)梯度(HESG)算法、遞階多新息增廣隨機(jī)梯度(HMIESG)算法、遞階增廣梯度(HEG)算法、遞階多新息增廣梯度(HMI-EG)算法、遞階增廣最小二乘(HELS)算法、遞階多新息增廣最小二乘(HMIELS)算法等。文獻(xiàn)[14]提出了CMA系統(tǒng)的交互隨機(jī)梯度辨識算法。本研究提出的遞階辨識方法可以推廣用于其他線性和非線性隨機(jī)系統(tǒng)，以及信號模型的參數(shù)辨識[15-25]。

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

有限脈沖響應(yīng)滑動平均模型描述的系統(tǒng)稱為有限脈沖響應(yīng)滑動平均系統(tǒng)，受控滑動平均模型描述的系統(tǒng)稱為受控滑動平均系統(tǒng)，余同。

考慮下列有限脈沖響應(yīng)滑動平均模型(FIRMA模型)描述的隨機(jī)控制系統(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)輸入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。由式(5)和(2)可以寫為

定義虛擬輸出變量

由式(8)可得

式(9)、(6)是FIR-MA系統(tǒng)(2)的基于辨識模型分解得到的兩個子辨識模型(sub-identification model)，即遞階辨識模型。它們可以看作兩個子系統(tǒng)，一個包含了系統(tǒng)模型的參數(shù)向量b，一個包含了噪聲模型的參數(shù)向量d。辨識模型式(8)分解為子辨識模型式(9)、(6)的遞階結(jié)構(gòu)如圖1所示。由于分解，導(dǎo)致兩個子辨識模型間存在耦合變量(關(guān)聯(lián)項)b和d。因此在推導(dǎo)兩個子辨識算法時，需要根據(jù)遞階辨識原理，協(xié)調(diào)它們間的關(guān)聯(lián)項?；诜纸獾谋孀R方法也稱為遞階辨識方法。

圖1 辨識模型分解為子辨識模型的梯階結(jié)構(gòu)Fig．1 Hierarchical structure of decomposing the identification model

2 遞階增廣隨機(jī)梯度辨識算法

對于FIR-MA系統(tǒng)的遞階辨識模型式(9)和式(6)，定義兩個梯度準(zhǔn)則函數(shù)(gradient criterion function):

然而，算法式(10)～(13)不可實(shí)現(xiàn)，因?yàn)橛疫叞宋粗畔⑾蛄喀?t)，參數(shù)向量b和d，這里解決方法是應(yīng)用遞階辨識原理，細(xì)節(jié)如下。

從式(8)可得

上式中未知的ψ(t)，b和d分別用其估計和代替，則v(t)的估計(即殘差)可由式(15)計算:

式(10)～(13)中未知ψ(t)，d和b分別用其估計和代替，得到式(16)～(19)，聯(lián)立式(3)，(14)～(15)，便得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的遞階增廣隨機(jī)梯度算法(hierarchical extended stochatic gradient algorithm，HESG算法):

圖2 遞階增廣隨機(jī)梯度(HESG)算法計算流程Fig．2 Flowchart of the HESG algorithm

1)初始化:令t＝1。置初值(0)＝1nb/p0，(0)＝1n d/p0，r1(0)＝1，r2(0)＝1，u(t－i)＝0，(t－i)＝1/p0，i＝1，2，…，max[n b，n d]，p0＝106。給定參數(shù)估計精度ε。

2)采集輸入輸出數(shù)據(jù)u(t)和y(t)。用式(21)構(gòu)造信息向量φ(t)，用式(22)構(gòu)造噪聲信息向量(t)。

3)用式(17)計算r1(t)，用式(19)計算r2(t)。

注1對于HESG算法式(16)～(24)，如果取式(16)中的r1(t)與式(18)中的r2(t)相等，且都為

r(t)＝r(t－1)＋，r(0)＝1，

那么HESG算法就退化為ESG算法。

式(16)～(24)是基于殘差的HESG算法，下面是基于新息的HESG算法:

3 遞階多新息增廣隨機(jī)梯度算法

設(shè)正整數(shù)p表示新息長度?；诙嘈孪⒈孀R理論[16]和HESG辨識算法(16)～(24)，將系統(tǒng)輸出y(t)，信息向量φ(t)和擴(kuò)展為堆積輸出向量Y(p，t)，堆積信息矩陣Φ(p，t)和:

將式(16)和式(18)中標(biāo)量新息(scalar innovation)

擴(kuò)展為新息向量(innovation vector)

注意到Y(jié)(1，t)＝y(tǒng)(t)，E(1，t)＝e(t)，Φ(1，t)＝φ(t)，式(16)和式(18)可以等價表達(dá)為

這是新息長度p＝1的兩階段“多新息”隨機(jī)梯度算法。將兩式中Φ(1，t)和E(1，t)里的1換為p，得到式(38)和式(41)，聯(lián)立式(34)～(37)，(17)和(19)～(24)，便得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的遞階多新息增廣隨機(jī)梯度算法(hierarchical multi-innovation extended stochastic gradient algorithm，HMI-ESG算法):

當(dāng)新息長度p＝1時，HMI-ESG算法退化為HESG算法(16)～(24)。

1)初始化:令t＝1，給定新息長度p和參數(shù)估計精度ε。置初值/p0，r1(0)＝1，r2(0)＝1，y(t－i)＝0，u(t－i)＝0，(t－i)＝1/p0，i＝1，2，…，p＋max[n b，n d]，p0＝106。

2)采集輸入輸出數(shù)據(jù)u(t)和y(t)。用式(47)和(48)構(gòu)造信息向量φ(t)和(t)，用式(44)～(46)構(gòu)造堆積輸出向量Y(p，t)，堆積信息矩陣Φ(p，t)和(p，t)。

3)用式(39)計算新息向量E(p，t)，用式(40)計算r1(t)。

5)用式(42)計算r2(t)，根據(jù)式(41)刷新參數(shù)估計向量(t)。用式(43)計算殘差(t)。

式(38)～(50)是基于殘差的HMI-ESG算法，下面是辨識FIR-MA系統(tǒng)參數(shù)向量b和d的基于新息的HMI-ESG算法:

隨機(jī)梯度算法(包括多新息隨機(jī)梯度算法、遞階隨機(jī)梯度算法、遞階增廣隨機(jī)梯度算法)的收斂速度很慢，為改善隨機(jī)梯度類算法的性能，可以引入收斂指數(shù)、遺忘因子等提高算法的性能。

在HMI-ESG算法中引入收斂指數(shù)、遺忘因子等，得到對應(yīng)的修正遞階多新息隨機(jī)梯度(M-HMIESG)算法、遺忘因子遞階多新息隨機(jī)梯度(FF-2SMISG)算法、遺忘因子修正遞階多新息隨機(jī)梯度(FF-M-HMI-ESG)算法等。

4 遞階遞推增廣梯度辨識算法

對于FIR-MA系統(tǒng)的遞階辨識模型式(9)和式(6)，定義準(zhǔn)則函數(shù):

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

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

分別求準(zhǔn)則函數(shù)J3(b)和J4(d)對參數(shù)向量b和d的梯度，得到

定義遞推關(guān)系:

于是，梯度可以表示為

關(guān)系式(70)～(77)不可實(shí)現(xiàn)，因?yàn)橛疫叞宋粗畔⑾蛄喀?t)，d和b。解決的辦法是用它們的估計(t)(t－1)和(t－1)代替，聯(lián)立式(20)～(24)，便得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的遞階增廣梯度算法(hierarchical extended gradient algorithm，HEG算法):

這是基于殘差的HEG算法。將算法(78)～(90)中的式(86)和式(88)修改為就得到基于新息的HEG算法。

圖3 遞階增廣梯度(HEG)算法計算流程Fig．3 Flowchart of the HEG algorithm

1)初始化:取μ＝1，給定ε。令t＝1，置初值1/p0，i＝1，2，…，max[n b，n d]，p0＞0是一個大常數(shù)，如p0＝106。

2)采集輸入輸出數(shù)據(jù)u(t)和y(t)。用式(87)構(gòu)造信息向量φ(t)，用式(88)構(gòu)造信息向量(t)。

3)用式(79)計算r1(t)，用式(80)計算向量ξ1(t)，用式(81)計算矩陣R1(t)。用式(78)刷新參數(shù)估計向量(t)。

4)用式(83)計算r2(t)，用式(84)計算向量，用式(85)計算矩陣R2(t)。用式(82)刷新參數(shù)估計向量(t)。用式(86)計算殘差(t)。

5 遞階多新息增廣梯度辨識算法

設(shè)正整數(shù)p表示新息長度。對于FIR-MA系統(tǒng)的遞階辨識模型(9)和(6)，利用從t－p＋1到t的最新p組數(shù)據(jù)，定義堆積輸出向量Y(p，t)和Y1(p，t)，堆積噪聲向量W(p，t)，堆積信息矩陣Φ(p，t)和信息矩陣Ψ(p，t)如下:

定義兩個多新息準(zhǔn)則函數(shù):

定義堆積輸出信息向量Zi，t和堆積信息矩陣Ωi，t如下:

則準(zhǔn)則函數(shù)J5(b)和J6(d)可以寫為

求準(zhǔn)則函數(shù)J5(b)和J6(d)分別對參數(shù)向量b和d的一階偏導(dǎo)數(shù)，可以得到兩個準(zhǔn)則函數(shù)的梯度向量:

其中

設(shè)μ1(t)≥0和μ2(t)≥0為收斂因子。根據(jù)負(fù)梯度搜索，極小化準(zhǔn)則函數(shù)J5(b)和J6(d)，可以得到下列梯度遞推關(guān)系:

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

據(jù)此可取μi(t)＝1/r i(t)，r i(0)＝1，i＝1，2，且

由于式(95)～(102)右邊存在未知量Ψ(p，t)，ψ(t)，d和b，所以無法計算參數(shù)估計(t)和(t)，解決的辦法是應(yīng)用遞階辨識原理，那些未知量用其估計和(t－1)代替，聯(lián)立式(44)～(50)，便得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的遞階多新息增廣梯度算法(hierarchical multi-innovation extended gradient algorithm，HMI-EG算法):

這是基于殘差的HMI-EG算法。將算法(103)～(118)中的式(111)和式(116)修改為

就得到基于新息的HMI-EG算法。

1)初始化:給定新息長度p和參數(shù)估計精度ε。令t＝1，置初值max[n b，n d]，p0＞0是一個大常數(shù)，如p0＝106。

2)采集輸入輸出數(shù)據(jù)u(t)和y(t)。用式(115)構(gòu)造信息向量φ(t)，用式(116)構(gòu)造信息向量(t)。用式(112)～(114)構(gòu)造堆積輸出向量Y(p，t)，堆積信息矩陣Φ(p，t)和(p，t)。

3)用式(104)計算r1(t)，用式(105)計算向量ξ1(t)，用式(106)計算矩陣R1(t)。用式(103)刷新參數(shù)估計向量(t)。

4)用式(108)計算r2(t)，用式(109)計算向量ξ2(t)，用式(110)計算矩陣R2(t)。用式(107)刷新參數(shù)估計向量(t)。用式(111)計算殘差(t)。

在HMI-EG算法(103)～(118)中引入加權(quán)因子(加權(quán)矩陣)、遺忘因子，便得到加權(quán)遞階多新息增廣梯度(W-HMI-EG)算法、遺忘因子遞階多新息增廣梯度(FF-HMI-EG)算法、加權(quán)遺忘因子遞階多新息增廣梯度(W-FF-HMI-EG)算法。

6 遞階增廣最小二乘辨識算法

參考遞推最小二乘算法的推導(dǎo)，極小化準(zhǔn)則函數(shù)J3(b)和J4(d)，可得到以下最小二乘遞推關(guān)系:

然而，同一個問題出現(xiàn):算法(121)～(129)不可實(shí)現(xiàn)，因?yàn)橛疫叞宋粗畔⑾蛄喀?t)，參數(shù)向量d和b。解決的辦法是用它們的估計(t)(t－1)和(t－1)代替，得到式(130)～(135)，聯(lián)立式(20)～(24)，便得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的遞階增廣最小二乘算法(hierarchical extended least squares algorithm，HELS算法):

這個算法也稱為基于殘差的HELS算法，讀者可研究它的收斂性[3]。下面給出辨識FIR-MA系統(tǒng)參數(shù)向量b和d的基于新息的遞階增廣最小二乘算法(HELS算法):

圖4 遞階增廣最小二乘(HELS)算法計算流程Fig．4 Flowchart of the HELS algorithm

1)初始化:令t＝1。置初值max[n b，n d]，p0＝106。給定精度ε。

2)采集輸入輸出數(shù)據(jù)u(t)和y(t)。用式(137)構(gòu)造信息向量φ(t)，用式(138)構(gòu)造噪聲信息向量(t)。

3)用式(131)和式(134)計算增益向量L1(t)和L2(t)，用式(132)和式(135)計算協(xié)方差陣P1(t)和P2(t)。

7 遞階多新息增廣最小二乘算法

設(shè)p為新息長度。根據(jù)多新息辨識理論，基于HELS算法(130)～(140)，將輸出y(t)，信息向量φ(t)和(t)擴(kuò)展為堆積輸出向量Y(p，t)，堆積信息矩陣Φ(p，t)和(p，t):

將式(130)和式(133)中標(biāo)量新息(scalar innovation)

擴(kuò)展為新息向量

便得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的遞階多新息增廣最小二乘算法(hierarchical multi-innovation extended least squares algorithm，HMI-ELS算法):

HMI-ELS算法(152)～(166)計算參數(shù)估計向量(t)和(t)的步驟如下。

1)初始化:給定新息長度p和估計精度ε。令t＝1。置初值

2)采集輸入輸出數(shù)據(jù)u(t)和y(t)。用式(163)構(gòu)造信息向量φ(t)，用式(164)構(gòu)造噪聲信息向量(t)。用式(160)～(162)構(gòu)造堆積輸出向量Y(p，t)，堆積信息矩陣Φ(p，t)和(p，t)。

3)用式(153)計算新息向量E(p，t)，用式(154)和式(157)計算增益向量L1(t)和L2(t)，用式(155)和(158)計算協(xié)方差陣P1(t)和P2(t)。

式(152)～(166)是基于殘差的HMI-ELS算法，下面是辨識FIR-MA系統(tǒng)參數(shù)向量b和d的基于新息的遞階多新息增廣最小二乘算法(HMI-ELS算法):

8 結(jié) 語

針對最簡單的有色噪聲干擾的有限脈沖響應(yīng)滑動平均模型，利用系統(tǒng)的觀測輸入輸出數(shù)據(jù)，基于遞階辨識原理，對辨識模型進(jìn)行分解，得到兩個遞階辨識子模型，一個包含系統(tǒng)模型的參數(shù)，一個包含噪聲模型的參數(shù)，推導(dǎo)出遞階增廣隨機(jī)梯度辨識算法和遞階多新息增廣隨機(jī)梯度辨識算法，遞階增廣梯度辨識算法和遞階多新息增廣梯度辨識算法，以及遞階增廣最小二乘辨識算法和遞階多新息增廣最小二乘辨識算法等。提出的分解辨識思想可以推廣到大規(guī)模線性隨機(jī)系統(tǒng)和復(fù)雜非線性隨機(jī)系統(tǒng)的辨識，其干擾可以是有色噪聲。