滤波辨识(9): 多变量CARARMA系统的滤波递阶广义增广参数辨识

丁 锋,万立娟,栾小丽,徐 玲,刘喜梅

(1.江南大学 物联网工程学院,江苏 无锡 214122;2.青岛科技大学 自动化与电子工程学院,山东 青岛 266061)

只有一个输入一个输出的系统称为单输入单输出系统,简称为标量系统。与标量系统相对应的是多变量系统。多变量系统有多个输入和(或)多个输出,也称为多输入多输出系统。它包括多输入单输出系统,单输入多输出系统。为了便于区分辨识方法的类别,如果在辨识中将输入向量和输出向量作为一个整体对待,这样的系统称为多变量系统。将包含了特征多项式的多变量系统称为类多变量系统。如果在辨识中使用了输入向量和输出向量的分量,例如按照输出的数目将多变量系统分解为多个子系统进行辨识,用这样的处理方法的多变量系统称为多输入多输出系统。由此分为多变量系统辨识方法,类多变量系统辨识方法,多输入多输出系统辨识方法三类。

与辅助模型辨识思想、多新息辨识理论、递阶辨识原理一样[1-7],滤波辨识理念也可以与梯度搜索、最小二乘搜索、牛顿搜索方法相结合,研究和提出不同的滤波递推辨识方法和滤波迭代辨识方法。《青岛科技大学学报(自然科学版)》上的连载论文研究了有限脉冲响应滑动平均系统、方程误差自回归系统、输出误差自回归滑动平均系统滤波递推辨识方法和滤波迭代辨识方法[8-14]。最近的连载论文研究了类多变量方程误差自回归滑动平均(M-EEARMA-like)系统,即类多变量受控自回归自回归滑动平均(M-CARARMA-like)系统的滤波递阶广义增广参数辨识方法[15],以及类多变量输出误差自回归滑动平均(M-OEARMA-like)系统,即类多变量Box-Jenkins系统的滤波辅助模型递阶广义增广递推参数辨识方法[16]。

本研究利用滤波辨识理念和递阶辨识原理,研究多变量受控自回归自回归滑动平均(M-CARARMA)系统(即多变量方程误差自回归滑动平均(MEEARMA)系统)的滤波递阶广义增广递推辨识方法。相关工作参见文献[2,7,17-27]。

1 多变量受控自回归ARMA 系统

考虑多变量受控自回归自回归滑动平均模型(multivariable controlled autoregressive autoregressive moving average model,M-CARARMA 模型),即多变量方程误差自回归滑动平均模型(multivariable equation-error autoregressive moving average model,M-EEARMA模型)。

描述的多变量系统:

这里考虑噪声w(t)是一个自回归滑动平均过程:

定义输出信息矩阵ψy(t),输入信息矩阵ψu(t)和滤波信息向量φ(t)如下:

这就是多变量受控自回归ARMA 系统(1)的滤波递阶辨识模型,它不仅包含了参数向量ρ,而且包含了参数矩阵θ。因此需要利用递阶辨识原理将其分解,从而研究递阶辨识方法。

对于多变量受控自回归ARMA 系统(1),可得到的只有系统输入输出数据{u(t),y(t)},在对应的滤波递阶辨识模型(11)中,信息矩阵ψy(t)中的输出向量y(t-i)是已知的,而滤波输出向量yf(t)是未知的,信息向量φ(t)中滤波输出向量yf(t)和滤波输入向量uf(t)都是未知的。这是辨识的困难所在。因此必须借助于滤波辨识理念和递阶辨识原理,利用系统的观测数据{u(t),y(t):t=1,2,3,…},通过对未知滤波输入向量uf(t),滤波输出向量yf(t)进行估算,研究M-CARARMA 系统的滤波递阶(多新息)广义增广随机梯度辨识方法、滤波递阶(多新息)广义增广递推梯度辨识方法、滤波递阶(多新息)广义增广最小二乘辨识方法,来辨识系统参数向量ρ和参数矩阵θ。

因为辨识模型涉及到未知滤波输入向量uf(t)和滤波输出向量yf(t),所以必须采用递推方法或迭代方法交替估算这些未知量,以及系统参数向量ρ和参数矩阵θ。下面先构造估算这些未知向量和未知矩阵的辅助模型。

2 估算未知向量和未知矩阵的辅助模型

因为滤波递阶辨识模型中涉及未知滤波输入向量uf(t-i)和未知滤波输出向量yf(t-i),所以辨识方案需要采用递推算法或迭代算法来交替估计参数向量和参数矩阵,以及未知滤波向量。这种未知量用其估计代替的方法在辨识中是极其常用的。这种先进方法在其它科学中是极其少见的。

3 滤波递阶广义增广随机梯度辨识方法

基于滤波递阶辨识模型(11),定义准则函数

注意到信息矩阵ψy(t)和信息向量φ(t)中含有未知滤波输出向量yf(t-i)和未知滤波输入向量uf(t-i),所以式(17)～(19)无法实现。因此,式(17)～(19)中未知信息矩阵ψy(t)和信息向量φ(t)用其估计代替,得到式(20)～(23),联立式(12)～(14)和辅助模型(15)～(16),便得到辨识多变量受控自回归ARMA 系统(1)对应的滤波递阶辨识模型(11)参数向量ρ和参数矩阵θ的滤波递阶广义增广随机梯度算法(filtered hierarchical generalized extended stochastic gradient algorithm,F-HGESG算法)[20-21]:

从F-HGESG 辨识算法,可以得到一些特殊的滤波递阶随机梯度辨识算法。

1)当nγ=0和nδ=0,即γ(z)=1和δ(z)=1时,F-HGESG 辨识算法(20)～(30)退化为多变量受控自回归(M-CAR)系统的多变量随机梯度算法(multivariable stochastic gradient algorithm,MSG 算法)[2,7,17,20]。

2)当nγ=0,即γ(z)=1时,F-HGESG 辨识算法(20)～(30)退化为多变量受控自回归滑动平均(M-CARMA)系统的滤波递阶增广随机梯度算法(filtered hierarchical extended stochastic gradient algorithm,F-HESG 算法)[7,20]。

3)当nδ=0,即δ(z)=1时,F-HGESG 辨识算法(20)～(30)退化为多变量受控自回归自回归(M-CARAR)系统的滤波递阶广义随机梯度算法(filtered hierarchical generalized stochastic gradient algorithm,F-HGSG 算法)。

4)读者可以写出多变量受控自回归自回归滑动平均(M-CARARMA)系统的加权滤波递阶广义增广随机梯度(W-F-HGESG)算法或滤波加权递阶广义增广随机梯度(F-W-HGESG)算法、遗忘因子滤波递阶广义增广随机梯度(FF-F-HGESG)算法或滤波遗忘因子递阶广义增广随机梯度(F-FFHGESG)算法、加权遗忘因子滤波递阶广义增广随机梯度(W-FF-F-HGESG)算法或加权滤波遗忘因子递阶广义增广随机梯度(W-F-FF-HGESG)算法。

例如,为提高参数估计精度和收敛速率,可在F-HGESG 算法(20)～(30)中引入遗忘因子λ,即将式(20)修改为

或

就得到滤波遗忘因子递阶广义增广随机梯度算法(filtered forgetting factor HGESG algorithm,FFF-HGESG算法)。

在F-HGESG 辨识算法(20)～(30)中,e(t)∈ℝm为新息向量,定义残差向量

4 滤波递阶多新息广义增广随机梯度方法

设p为新息长度。根据多新息辨识理论[6,28],基于F-HGESG 算法(20)～(30),按照式(46)～(51)定义堆积输出向量Y(p,t),堆积输出矩阵Y1(p,t),堆积信息矩阵Φ(p,t),堆积信息向量Φ1(p,t),堆积信息矩阵Ψ(p,t)和Ψ1(p,t),将式(20)中信息矩阵扩展为堆积信息矩阵Ψ(p,t),新息向量e(t)扩展为一个大新息向量E(p,t),得到

联立式(36)～(40)和式(24)～(30),就得到辨识多变量受控自回归ARMA 系统(1)参数向量ρ和参数矩阵θ的滤波递阶多新息广义增广随机梯度算法(filtered hierarchical multi-innovation generalized extended stochastic gradient algorithm,F-HMIGESG 算法)[20]:

从F-HMI-GESG 辨识算法,可以得到一些特殊的滤波递阶多新息随机梯度辨识算法。

1)当nγ=0和nδ=0时,F-HMI-GESG 辨识算法(41)～(58)退化为多变量受控自回归(M-CAR)系统的多变量多新息随机梯度算法(multivariable multi-innovation stochastic gradient algorithm,MMISG 算法)。

2)当nγ=0时,F-HMI-GESG 辨识算法(41)～(58)退化为多变量受控自回归滑动平均(M-CARMA)系统的滤波递阶多新息增广随机梯度算法(filtered hierarchical multi-innovation extended stochastic gradient algorithm,F-HMI-ESG 算法)。

3)当nδ=0时,F-HMI-GESG 辨识算法(41)～(58)退化为多变量受控自回归自回归(M-CARAR)系统的滤波递阶多新息广义随机梯度算法(filtered hierarchical multi-innovation generalized stochastic gradient algorithm,F-HMI-GSG 算法)。

4)读者可以写出多变量受控自回归自回归滑动平均(M-CARARMA)系统的加权滤波递阶多新息广义增广随机梯度(W-F-HMI-GESG)算法或滤波加权递阶多新息广义增广随机梯度(F-W-HMIGESG)算法、遗忘因子滤波递阶多新息广义增广随机梯度(FF-F-HMI-GESG)算法或滤波遗忘因子递阶多新息广义增广随机梯度(F-FF-HMI-GESG)算法、加权遗忘因子滤波递阶多新息广义增广随机梯度(W-FF-F-HMI-GESG)算法或加权滤波遗忘因子递阶多新息广义增广随机梯度(W-F-FF-HMIGESG)算法。

例如,为提高参数估计精度和收敛速率,可在F-HMI-GESG 算法(41)～(58)中引入遗忘因子λ,即将式(45)修改为

就得到滤波遗忘因子递阶多新息广义增广随机梯度算法(filtered forgetting factor HMI-GESG algorithm,F-FF-HMI-GESG 算法)。当然,一些梯度类算法都可以引入收敛指数、加权因子、遗忘因子等得到相应的辨识算法,不一一论述[2-3]。

3)用式(46)～(47)构造堆积输出向量Y(p,t)和堆积输出矩阵Y1(p,t),用式(48)～(49)构造堆积信息矩阵Φ(p,t)和堆积信息向量Φ1(p,t),用式(50)～(51)构造堆积信息矩阵Ψ(p,t)和Ψ1(p,t)。

5 滤波递阶广义增广递推梯度辨识方法

根据多变量受控自回归ARMA 系统(1)的滤波递阶辨识模型(11),定义准则函数

仿照递阶梯度算法的推导[5,7],使用梯度搜索,极小化准则函数J2(ρ,θ),可以得到下列梯度递推关系:

式(63)～(70)中的未知信息矩阵ψy(t)和信息向量φ(t)用其估计代替,联立式(24)～(30),便得到辨识多变量受控自回归ARMA 系统(1)参数向量ρ和参数矩阵θ的滤波递阶广义增广递推梯度算法(filtered hierarchical generalized extended recursive gradient algorithm,F-HGERG 算法),简称为滤波递阶广义增广梯度算法(filtered hierarchical generalized extended gradient algorithm,F-HGEG 算法)[13]:

从F-HGEG 辨识算法,可以得到一些特殊的滤波递阶递推梯度辨识算法。

1)当nγ=0和nδ=0时,F-HGEG 算法(71)～(85)退化为多变量受控自回归(M-CAR)系统的多变量递推梯度算法(multivariable recursive gradient algorithm,M-RG算法)[7]。

2)当nγ=0时,F-HGEG 算法(71)～(85)退化为多变量受控自回归滑动平均(M-CARMA)系统的滤波递阶增广梯度算法(filtered hierarchical extended gradient algorithm,F-HEG 算法)[7]。

3)当nδ=0时,F-HGEG 算法(71)～(85)退化为多变量受控自回归自回归(M-CARAR)系统的滤波递阶广义梯度算法(filtered hierarchical generalized gradient algorithm,F-HGG 算法)。

4)读者可以写出多变量受控自回归自回归滑动平均(M-CARARMA)系统的加权滤波递阶广义增广梯度(W-F-HGEG)算法或滤波加权递阶广义增广梯度(F-W-HGEG)算法、遗忘因子滤波递阶广义增广梯度(FF-F-HGEG)算法或滤波遗忘因子递阶广义增广梯度(F-FF-HGEG)算法、加权遗忘因子滤波递阶广义增广梯度(W-FF-F-HGEG)算法或加权滤波遗忘因子递阶广义增广梯度(W-F-FF-HGEG)算法。

引理3 对于F-HGEG 算法(71)～(85),定义新息向量E(t)和新息矩阵E1(t)如下:

它们有下列关系:

6 滤波递阶多新息广义增广递推梯度方法

从F-HMI-GEG 辨识算法,可以得到一些特殊的滤波递阶多新息梯度辨识算法。

1)当nγ=0 和nδ=0 时,F-HMI-GEG 算法(90)～(110)退化为多变量受控自回归(M-CAR)系统的多变量多新息递推梯度算法(multivariable multi-innovation recursive gradient algorithm,MMIRG 算法)[7]。

2)当nγ=0时,F-HMI-GEG 算法(90)～(110)退化为多变量受控自回归滑动平均(M-CARMA)系统的滤波递阶多新息增广梯度算法(filtered hierarchical multi-innovation extended gradient algorithm,F-HMI-EG 算法)[7]。

3)当nδ=0时,F-HMI-GEG 算法(90)～(110)退化为多变量受控自回归自回归(M-CARAR)系统的滤波递阶多新息广义梯度算法(filtered hierarchical multi-innovation generalized gradient algorithm,F-HMI-GG 算法)。

4)读者可以写出多变量受控自回归自回归滑动平均(M-CARARMA)系统的加权滤波递阶多新息广义增广递推梯度(W-F-HMI-GERG)算法或滤波加权递阶多新息广义增广梯度(F-W-HMI-GEG)算法、遗忘因子滤波递阶多新息广义增广递推梯度(FF-F-HMI-GERG)算法或滤波遗忘因子递阶多新息广义增广梯度(F-FF-HMI-GEG)算法、加权遗忘因子滤波递阶多新息广义增广递推梯度(W-FF-FHMI-GERG)算法或加权滤波遗忘因子递阶多新息广义增广梯度(W-F-FF-HMI-GEG)算法。

3)用式(98)～(99)构造堆积输出向量Y(p,t)和堆积输出矩阵Y1(p,t),用式(100)～(101)构造堆积信息矩阵Φ(p,t)和堆积信息向量Φ1(p,t)。用式(102)～(103)构造堆积信息矩阵Ψ(p,t)和Ψ1(p,t)。

4)用式(91)计算r1(t),用式(92)计算向量ξ1(t),用式(93)计算矩阵R1(t),用式(90)刷新参数估计向量

5)用式(95)计算r2(t),用式(96)计算向量ξ2(t),用式(97)计算矩阵R2(t),用式(94)刷新参数估计矩阵

7 滤波递阶广义增广最小二乘辨识方法

根据滤波递阶辨识原理,令J2(ρ,θ)分别对ρ和θ的偏导数为0,可推导出最小二乘递推关系[3]:

式(111)～(116)中的未知信息矩阵ψy(t)和信息向量φ(t)用其估计代替,联立式(24)～(30),便得到辨识多变量受控自回归ARMA 系统(1)参数向量ρ和参数矩阵θ的滤波递阶广义增广最小二乘算法(filtered hierarchical generalized extended least squares algorithm,F-HGELS 算法)[18-20]:

从F-HGELS辨识算法,可以得到一些特殊的滤波递阶最小二乘辨识算法。

1)当nγ=0和nδ=0时,F-HGELS算法(117)～(130)退化为多变量受控自回归(M-CAR)系统的多变量递推最小二乘算法(multivariable recursive least squares algorithm,M-RLS算法)[2,3,7,18-20]。

2)当nγ=0时,F-HGELS算法(117)～(130)退化为多变量受控自回归滑动平均(M-CARMA)系统的滤波递阶增广最小二乘算法(filtered hierarchical extended least squares algorithm,F-HELS算法)[7,20]。

3)当nδ=0时,F-HGELS 算法(117)～(130)退化为M-CARAR 系统的滤波递阶广义最小二乘算法(filtered hierarchical generalized least squares algorithm,F-HGLS算法)。

4)读者可以写出多变量受控自回归自回归滑动平均(M-CARARMA)系统的加权滤波递阶广义增广最小二乘(W-F-HGELS)算法或滤波加权递阶广义增广最小二乘(F-W-HGELS)算法、遗忘因子滤波递阶广义增广最小二乘(FF-F-HGELS)算法或滤波遗忘因子递阶广义增广最小二乘(F-FF-HGELS)算法、加权遗忘因子滤波递阶广义增广最小二乘(W-FF-F-HGELS)算法或加权滤波遗忘因子递阶广义增广最小二乘(W-F-FF-HGELS)算法。

3)用式(118)计算新息向量e(t),用式(119)计算增益矩阵L1(t),用式(120)计算协方差矩阵P1(t),用式(117)刷新参数估计向量,用式(122)计算增益矩阵L2(t),用式(123)计算协方差矩阵P2(t),用式(121)刷新参数估计矩阵

在F-HGELS辨识算法(117)～(130)中,e(t)∈ℝm为新息向量,定义残差向量

引理4对于F-HGELS辨识算法(117)～(130),新息向量e(t)与残差向量有下列关系:

8 滤波递阶多新息广义增广最小二乘方法

设p为新息长度。根据多新息辨识理论[6,28],基于F-HGELS算法(117)～(130),按照式(146)～(151)定义堆积输出向量Y(p,t),堆积输出矩阵Y1(p,t),堆积输入信息矩阵Φ(p,t),堆积输入信息向量Φ1(p,t),堆积信息矩阵Ψ(p,t)和Ψ1(p,t),将式(117)中新息向量e(t)扩展为一个大新息向量E(p,t),得到

将式(121)中新息向量e(t)扩展为新息矩阵E1(p,t),得到

从F-HMI-GELS辨识算法,可以得到一些特殊的滤波递阶多新息最小二乘辨识算法。

1)当nγ=0 和nδ=0 时,HMI-GELS 算法(138)～(158)退化为M-CAR 系统的多变量多新息最小二乘算法(multivariable multi-innovation least squares algorithm,M-MILS算法)[7,20]。

2)当nγ=0 时,F-HMI-GELS 算法(138)～(158)退化为M-CARMA 系统的滤波递阶多新息增广最小二乘算法(filtered hierarchical multi-innovation extended least squares algorithm,F-HMI-ELS算法)[7,20]。

3)当nδ=0 时,F-HMI-GELS 算法(138)～(158)退化为M-CARAR 系统的滤波递阶多新息广义最小二乘算法(filtered hierarchical multi-innovation generalized least squares algorithm,F-HMIGLS算法)。

4)读者可以写出多变量受控自回归自回归滑动平均(M-CARARMA)系统的加权滤波递阶多新息广义增广最小二乘(W-F-HMI-GELS)算法或滤波加权递阶多新息广义增广最小二乘(F-W-HMIGELS)算法、遗忘因子滤波递阶多新息广义增广最小二乘(FF-F-HMI-GELS)算法或滤波遗忘因子递阶多新息广义增广最小二乘(F-FF-HMI-GELS)算法、加权遗忘因子滤波递阶多新息广义增广最小二乘(W-FF-F-HMI-GELS)算法或加权滤波遗忘因子递阶多新息广义增广最小二乘(W-F-FF-HMIGELS)算法。

3)用式(146)～(147)构造堆积输出向量Y(p,t)和堆积输出矩阵Y1(p,t),用式(148)～(149)构造堆积信息矩阵Φ(p,t)和堆积信息向量Φ1(p,t)。用式(150)～(151)构造堆积信息矩阵Ψ(p,t)和Ψ1(p,t)。

4)用式(139)计算新息向量E(p,t),用式(140)计算增益矩阵L1(t),用式(141)计算协方差矩P1(t),用式(138)刷新参数估计向量

5)用式(143)计算新息矩阵E1(p,t),用式(144)计算增益矩阵L2(t),用式(145)计算协方差矩阵P2(t),用式(142)刷新参数估计矩阵

9 结语

针对多变量受控自回归自回归滑动平均(MCARARMA)系统,利用滤波辨识理念和递阶辨识原理,提出了滤波递阶(多新息)广义增广随机梯度辨识方法、滤波递阶(多新息)广义增广递推梯度辨识方法、滤波递阶(多新息)广义增广最小二乘辨识方法。这些滤波递阶广义增广辨识方法可以推广到其它有色噪声干扰下的线性和非线性多变量随机系统中[30-44]。