丁 鋒，劉喜梅

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

一些辨識新思想、新理論、新原理、新概念,如輔助模型辨識思想、多新息辨識理論、遞階辨識原理、耦合辨識概念、濾波辨識理念等[1-7]與傳統(tǒng)迭代辨識方法相結(jié)合,誕生出許多新辨識方法,如輔助模型梯度迭代辨識方法、多新息最小二乘迭代辨識方法、遞階牛頓迭代辨識方法等。從2017年開始在《青島科技大學(xué)學(xué)報》連載“信號建模”論文,從2018年開始連載“傳遞函數(shù)辨識”論文26篇[8-16]。最近的連載論文研究了有限脈沖響應(yīng)滑動平均系統(tǒng)的遞階遞推增廣參數(shù)辨識方法[9]和遞階增廣迭代參數(shù)辨識方法[10],以及濾波增廣遞推參數(shù)辨識方法[17]。本工作利用濾波辨識理念,研究有限脈沖響應(yīng)滑動平均系統(tǒng)的濾波增廣迭代辨識方法。

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è)階次n b和n d已知,記n＝n b＋n d。定義參數(shù)向量和信息向量:

定義不可測干擾滑動平均噪聲:

式(4)和(1)可以寫為

這里的目標是利用系統(tǒng)的觀測數(shù)據(jù){u(t),y(t)}或{y(t),φ(t)},基于數(shù)據(jù)濾波技術(shù),研究濾波迭代辨識方法,估計系統(tǒng)的參數(shù)向量b和d。

定義濾波輸出yf(t)和濾波輸入信息向量φf(t)如下:

它們可以進一步表示為

式(6)兩邊同除以D(z)得到

即

式(11)和(5)是FIR-MA系統(tǒng)的濾波辨識模型(filtered identification model)。它們分別包含了系統(tǒng)模型的參數(shù)向量b和噪聲模型的參數(shù)向量d。

2 濾波增廣梯度迭代方法

設(shè)L為數(shù)據(jù)長度。根據(jù)式(7),(11)和(5),定義堆積輸出向量Y(L)和堆積輸入信息矩陣Φ(L),堆積濾波輸出向量Yf(L)和堆積濾波輸入信息矩陣Φf(L),堆積噪聲向量W(L)和堆積噪聲矩陣Ψ(L)如下:

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

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

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

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

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

梯度迭代算法(14),(16)和(18)～(21)不可實現(xiàn),因為只有{y(t),φ(t)}是可得到的觀測數(shù)據(jù),而右邊堆積濾波輸出向量Yf(L),堆積濾波信息矩陣Φf(L),堆積噪聲矩陣Ψ(L),以及參數(shù)向量b都是未知的,故需要先估算這些未知量。

用噪聲v(t－i)的估計^v k(t－i)構(gòu)造ψ(t)的估計:

由式(8)可得

根據(jù)此式,用信息向量的估計(t)和第k次迭代的參數(shù)估計和代替上式中未知的ψ(t),b和d,可以得到白噪聲項v(t)的估計,即殘差:

用噪聲模型參數(shù)向量d的第k次迭代估計＝構(gòu)造多項式D(z)的迭代估計:

由于式(9)中D(z)是未知的,無法得到y(tǒng)f(t),故用其第k次迭代估計^D k(z)計算yf(t)的估計:

同理,由于式(10)中D(z)是未知的,無法得到φ(t)的濾波向量φf(t),故用其第k次迭代估計(z)計算φf(t)的估計:

根據(jù)Yf(L),Φf(L)和Ψ(L)的定義式,分別用yf(t),φf(t)和ψ(t)的估計(t)(t)和(t)來構(gòu)造它們的估計:

用(L)和(L)分別代替式(14)和(18)或(20)中未知Yf(L)和Φf(L),得到式(29)～(30)的梯度迭代估計,用(L)和前一次迭代估計分別代替式(16)和(19)或(21)中未知Ψ(L)和b,得到式(31)～(32)的梯度迭代估計,聯(lián)立式(12)～(13),(26)～(28),(2),(22),(24)～(25)和(23),就得到估計FIR-MA系統(tǒng)參數(shù)向量b和d的基于濾波的增廣梯度迭代辨識算法,簡稱濾波增廣梯度迭代算法(filtered extended gradient-based iterative algorithm,F-EGI算法):

或

或

F-EGI算法(29)～(44)計算參數(shù)估計和的步驟如下。

1)初始化:令迭代次數(shù)k＝1,給定數(shù)據(jù)長度L和參數(shù)估計精度ε,置參數(shù)估計初值＝1nb/p0,＝1n d/p0,p0＝106。置(t)＝1/p0(t)＝1nb/p0,(t)＝1/p0,t＝1,2,…,L。

2)采集觀測數(shù)據(jù)u(t)和y(t),用式(38)構(gòu)造信息向量φ(t),t＝1,2,…,L。用式(33)構(gòu)造堆積輸出向量Y(L),用式(34)構(gòu)造堆積輸入信息矩陣Φ(L)。

3)用式(39)構(gòu)造噪聲向量(t),t＝1,2,…,L。用式(37)構(gòu)造堆積噪聲矩陣(L)。

4)根據(jù)式(32)計算步長μk(t),用式(31)刷新參數(shù)估計向量。

5)用式(40)計算濾波輸出(t),用式(41)計算濾波信息向量(t),t＝1,2,…,L。

6)用式(35)構(gòu)造堆積濾波輸出向量(L),用式(36)構(gòu)造堆積濾波信息矩陣(L)。

7)根據(jù)式(30)計算步長μf,k,用式(29)刷新參數(shù)估計向量。

8)用式(42)計算白噪聲的估計^v k(t),t＝1,2,…,L。

設(shè)M1(L)∈?L×L和M2(L)∈?L×L為非負定加權(quán)矩陣。將準則函數(shù)J1(b)和J2(d)修改為加權(quán)準則函數(shù):

使用梯度搜索,極小化準則函數(shù)J′1(b)和J′2(d)可以得到

或

或

那么式(33)～(50)構(gòu)成了加權(quán)濾波增廣梯度迭代算法(W-F-EGI算法),或稱為濾波加權(quán)增廣梯度迭代算法(F-W-EGI算法)。進一步,設(shè)0＜λi≤1為遺忘因子,取加權(quán)矩陣為

就得到遺忘因子濾波增廣梯度迭代算法(FF-F-EGI算法),或稱為濾波遺忘因子增廣梯度迭代算法(FFF-EGI算法)。

定義迭代新息向量:

和迭代殘差向量:

引理1對于F-EGI辨識算法(29)～(44),新息向量Ef,k(L)和E k(L)與殘差向量(L)和(L)滿足關(guān)系:

3 濾波增廣最小二乘迭代方法

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

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

因為只有{y(t),[φ[(t)}是可得到的觀測數(shù)據(jù),而堆積濾波輸出向量Yf(L),堆積濾波輸入信息矩陣Φf(L)和堆積噪聲矩陣Ψ(L)都是未知的,所以上兩式無法求解出參數(shù)估計。解決的辦法是用式(26)～(27)中的(L)和(L)分別代替式(55)中未知的Yf(L)和Φf(L),得到式(57)的最小二乘迭代估計(t),用式(28)中的(L)和前一次迭代估計(t)分別代替式(56)中未知的Ψ(L)和b,得到式(58)的最小二乘迭代估計,聯(lián)立式(35)～(44),就得到估計FIR-MA系統(tǒng)參數(shù)向量b和d的基于濾波的增廣最小二乘迭代辨識算法,簡稱濾波增廣最小二乘迭代算法(filtered extended least squares-based iterative algorithm,F-ELSI算法):

F-ELSI算法(57)～(70)計算參數(shù)估計和的步驟如下。

1)初始化:令迭代次數(shù)k＝1,給定數(shù)據(jù)長度L和參數(shù)估計精度ε,置參數(shù)估計初值＝1nb/p0或^b0為隨機向量,p0＝106。置(t)＝1/p0(t)＝1nb/p0(t)＝1/p0,t＝1,2,…,L。隨機數(shù)初值是為保證式(57)～(58)中的矩陣可逆。

2)采集觀測數(shù)據(jù)u(t)和y(t),用式(64)構(gòu)造信息向量φ(t),t＝1,2,…,L。用式(59)構(gòu)造堆積輸出向量Y(L),用式(60)構(gòu)造堆積輸入信息矩陣Φ(L)。

3)用式(65)構(gòu)造噪聲向量(t),t＝1,2,…,L。用式(63)構(gòu)造堆積噪聲矩陣(L)。

4)用式(58)刷新參數(shù)估計向量。

5)用式(66)計算濾波輸出(t),用式(67)計算濾波信息向量(t),t＝1,2,…,L。

6)用式(61)構(gòu)造堆積濾波輸出向量(L),用式(62)構(gòu)造堆積濾波信息矩陣(L)。

7)用式(57)刷新參數(shù)估計向量^b k。用式(68)計算殘差(t),t＝1,2,…,L。

極小化準則函數(shù)J′1(b)和J′2(d),令其對b和d的梯度為零,可以推導(dǎo)出

那么式(59)～(72)構(gòu)成了加權(quán)濾波增廣最小二乘迭代算法(W-F-ELSI算法),或稱為濾波加權(quán)增廣最小二乘迭代算法(F-W-ELSI算法)。進一步,設(shè)0＜λi≤1為遺忘因子,取加權(quán)矩陣為

就得到遺忘因子濾波增廣最小二乘迭代算法(FF-FELSI算法),或稱為濾波遺忘因子增廣最小二乘迭代算法(F-FF-ELSI算法)。

4 濾波多新息增廣梯度迭代方法

設(shè)p為數(shù)據(jù)窗長度(p?n)。根據(jù)式(7),(11)和(5),定義堆積輸出向量Y(p,t)和堆積信息矩陣Φ(p,t),堆積濾波輸出向量Yf(p,t)和堆積濾波輸入信息矩陣Φf(p,t),堆積噪聲向量W(p,t)和堆積噪聲矩陣Ψ(p,t)如下:

注意到W(p,t)＝Y(jié)(p,t)－Φ(p,t)b∈?p。根據(jù)辨識模型(11)和(5),定義關(guān)于參數(shù)向量b和d的二次準則函數(shù):

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

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

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

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

梯度迭代算法(75),(77)和(79)～(82)不可實現(xiàn),因為堆積濾波輸出向量Yf(p,t),堆積濾波信息矩陣Φf(p,t)和堆積噪聲矩陣Ψ(p,t)都是未知的,故需要先構(gòu)造這些未知變量的估計,進而在辨識算法中用它們的估計代替這些未知變量,細節(jié)如下。

用噪聲v(t－i)的估計(t－i)定義噪聲向量ψ(t)的估計:

由式(8)可得

用第k次迭代的參數(shù)估計(t)和(t),以及信息向量的估計(j)分別代替上式中未知的b和d,以及噪聲向量ψ(j),便得到白噪聲項v(j)的估計:

用噪聲模型參數(shù)向量d在時刻t的第k次迭代估計

構(gòu)造多項式D(z)在時刻t的第k次迭代估計:

由于D(z)是未知的,無法得到y(tǒng)f(t),故用其估計(t,z)計算yf(t)的估計:

同理,由于式(10)中D(z)是未知的,無法得到φ(t)的濾波向量φf(t),故用其代估計^D k(t,z)計算φf(t)的估計:

根據(jù)Yf(p,t),Φf(p,t)和Ψ(p,t)的定義式,分別用yf(t),φf(t)和ψ(t)的估計(t)(t)和(t)構(gòu)造Yf(p,t),Φf(p,t)和Ψ(p,t)的估計:

用Yf,k(p,t)和Φf,k(p,t)分別代替式(75)和(79)或(81)中未知的Yf(p,t)和Φf(p,t),得到式(90)～(92)的梯度迭代估計(t),用(p,t)和當前時刻t的前一次迭代估計^b k－1(t)分別代替式(77)和(80)或(82)中未知的Ψ(p,t)和b,得到式(93)～(95)的梯度迭代估計(t),聯(lián)立式(73)～(74),(87)～(89),(2),(83),(85)～(86)和(84)就得到估計FIR-MA系統(tǒng)參數(shù)向量b和d的基于濾波的多新息增廣梯度迭代辨識算法,簡稱濾波多新息增廣梯度迭代算法(filtered multi-innovation extended gradient-based iterative algorithm,F-MI-EGI算法):

或

或

濾波多新息增廣梯度迭代算法又稱濾波移動數(shù)據(jù)窗增廣梯度迭代算法(filtered moving-data-window extended gradient-based iterative algorithm,FMDW-EGI算法)。

F-MI-EGI算法(90)～(107)計算參數(shù)估計(t)和(t)的步驟如下。

1)初始化:令t＝1,給定移動數(shù)據(jù)窗長度p(通常取p?n),參數(shù)估計精度ε,最大迭代次數(shù)kmax。置參數(shù)估計初值(t)＝1nb/p0,(t)＝1n d/p0,p0＝106。

2)令k＝1。置(j)＝1/p0(j)＝1/p0,(j)＝1nb/p0,j＝1,2,…,t。

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

4)用式(102)構(gòu)造噪聲向量(t),用式(100)構(gòu)造堆積噪聲矩陣(p,t)。

5)根據(jù)式(94)或(95)計算步長μk(t),用式(93)刷新參數(shù)估計向量(t)。

6)用式(103)計算濾波輸出(j),用式(104)計算濾波信息向量(j),j＝1,2,…,t。

7)用式(98)構(gòu)造堆積濾波輸出向量(p,t),用式(99)構(gòu)造堆積濾波信息矩陣(p,t)。

8)根據(jù)式(91)或(92)計算步長μf,k(t),用式(90)刷新參數(shù)估計向量(t)。

9)用式(105)刷新殘差(j),j＝1,2,…,t。

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

設(shè)Λ1(t)∈?p×p和Λ2(t)∈?p×p為非負定加權(quán)矩陣。將準則函數(shù)J3(b)和J4(d)修改為加權(quán)準則函數(shù):

使用梯度搜索,極小化準則函數(shù)J′3(b)和J′4(d)可以得到

或

或

那么式(96)～(113)構(gòu)成了加權(quán)濾波多新息增廣梯度迭代算法(weighted F-MI-EGI algorithm,W-FMI-EGI算法),或稱為濾波加權(quán)多新息增廣梯度迭代算法(F-W-MI-EGI算法)。

定義迭代新息向量:

和迭代殘差向量:

引理2對于F-MI-EGI辨識算法(90)～(107),新息向量Ef,k(p,t)和E k(p,t)與殘差向量(p,t)和(p,t)滿足關(guān)系:

5 濾波多新息增廣投影迭代方法

下面利用最速下降法推導(dǎo)濾波多新息增廣投影辨識方法。考慮從j＝t－p＋1到j(luò)＝t的數(shù)據(jù)窗里的p組數(shù)據(jù)。定義堆積噪聲向量:

基于FIR-MA系統(tǒng)的濾波辨識模型(11)和(5),定義滑動數(shù)據(jù)窗準則函數(shù):

令k＝1,2,3,…是迭代變量,μf,k(t)≥0和μk(t)≥0是迭代步長(收斂因子),設(shè)(t)∈?nb是參數(shù)向量b在當前時刻t的第k次迭代估計,是(t)∈?n d參數(shù)向量d在當前時刻t的第k次迭代估計。使用最速下降法(負梯度搜索),極小化準則函數(shù)J5(b)和J6(d),可以得到梯度遞推關(guān)系:

其中

下面求最佳步長(best step size)μf,k(t)。將b＝(t)代入準則函數(shù)J5(b)中,可得

極小化g(μf,k(t)),令g′(μf,k(t))＝0,可求得最佳步長為

同理可得

聯(lián)立式(118)～(123)和(96)～(107),就得到辨識FIR-MA系統(tǒng)參數(shù)向量b和d的基于濾波的多新息增廣投影迭代辨識算法,簡稱濾波多新息增廣投影迭代算法(filtered multi-innovation extended projection-based iterative algorithm,F-MI-EPI算法):

如果式(126)和(129)的分母為零,就令(t)＝因為收斂因子的計算式比較復(fù)雜,故對其進行簡化。由于對于任意實向量x和非負定對稱矩陣Q,不等式

成立,其中λmax[Q]為矩陣Q的最大特征值。于是,收斂因子可以保守取為

因為計算矩陣的跡(trace)比計算特征值簡單,所以收斂因子可以更保守取為

因此,將式(124)和(127)合并修改為

或則式(125),(128)和(130)～(145)構(gòu)成了簡化的濾波多新息增廣投影迭代算法(F-MI-EPI算法),它等同于F-MI-EGI算法(90)～(107)。

6 濾波多新息增廣最小二乘迭代方法

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

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

因為堆積濾波輸出向量Yf(p,t),堆積濾波信息矩陣Φf(p,t)和堆積噪聲矩陣Ψ(p,t)都是未知的,所以解決的方法是用式(87)～(88)中的(p,t)和(p,t)分別代替式(146)中未知的Yf(p,t)和Φf(p,t),得到式(148)的最小二乘迭代估計(t),用式(89)中的(p,t)和前一次迭代參數(shù)估計(t)分別代替式(147)中未知Ψ(p,t)和b,得到式(149)的最小二乘迭代估計(t),聯(lián)立式(130)～(141),就得到估計FIR-MA系統(tǒng)參數(shù)向量b和d的基于濾波的多新息增廣最小二乘迭代辨識算法,簡稱濾波多新息增廣最小二乘迭代算法(filtered multi-innovation extended least squaresbased iterative algorithm,F-MI-ELSI算法):

濾波多新息增廣最小二乘迭代算法又稱濾波移動數(shù)據(jù)窗增廣最小二乘迭代算法(filtered movingdata-window extended least squares-based iterative algorithm,F-MDW-ELSI算法)。

F-MI-ELSI算法(148)～(161)計算參數(shù)估計和的步驟如下。

1)初始化:令t＝l,給定移動數(shù)據(jù)窗長度p(通常取p?n),參數(shù)估計精度ε,最大迭代次數(shù)kmax,(t)＝1nb/p0,p0＝106。采集觀測數(shù)據(jù){y(j),φ(j):j＝1,2,…,l－1}l≥p,l≥p這里為一個正整數(shù)(請讀者思考l的意義)。

2)令k＝1。置(j)為隨機數(shù)(j)為隨機數(shù)(j)為隨機向量,j＝1,2,…,t。隨機數(shù)初值是為保證式(148)和(149)中的矩陣可逆。

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

4)用式(156)構(gòu)造噪聲向量(t),用式(154)構(gòu)造堆積噪聲矩陣(p,t)。

5)用式(149)刷新參數(shù)估計向量(t)。

6)用式(157)計算濾波輸出(j),用式(158)計算濾波信息向量(j),j＝1,2,…,t。

7)用式(152)構(gòu)造堆積濾波輸出向量(p,t),用式(153)構(gòu)造堆積濾波信息矩陣(p,t)。

8)用式(148)刷新參數(shù)估計向量(t)。

9)用式(159)刷新殘差(j),j＝1,2,…,t。

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

極小化準則函數(shù)J′3(b)和J′4(d),令其對b和d的梯度為零,可以推導(dǎo)出

那么式(150)～(163)構(gòu)成了加權(quán)濾波多新息增廣最小二乘迭代算法(weighted F-MI-ELSI algorithm,W-F-MI-ELSI算法),或稱為濾波加權(quán)多新息增廣最小二乘迭代算法(F-W-MI-ELSI算法)。

7 結(jié) 語

基于系統(tǒng)的觀測輸入輸出數(shù)據(jù),利用迭代搜索和濾波辨識理念,研究了有限脈沖響應(yīng)滑動平均系統(tǒng)的濾波增廣迭代辨識方法,包括濾波(多新息)增廣梯度迭代辨識方法、濾波(多新息)增廣最小二乘迭代辨識方法,以及濾波加權(quán)和遺忘因子增廣迭代辨識方法等。提出的濾波增廣迭代辨識方法可以推廣到其它有色噪聲干擾下的多變量隨機系統(tǒng)和非線性隨機系統(tǒng)中[18-28]。