滤波辨识(3):方程误差自回归系统的滤波递推广义参数估计

丁 锋,刘海波,刘喜梅

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

梯度方法和牛顿方法是研究优化问题的基本方法。用于系统辨识就得到梯度辨识方法和牛顿辨识方法。对于线性参数系统,牛顿辨识方法退化为最小二乘辨识方法。

梯度方法和牛顿方法与辅助模型辨识思想、多新息辨识理论、递阶辨识原理、耦合辨识概念等相结合[1-7],诞生出许多辨识方法,如辅助模型随机梯度方法、辅助模型最小二乘方法、辅助模型牛顿辨识方法、多新息随机梯度辨识方法、多新息最小二乘辨识方法、多新息牛顿辨识方法、递阶随机梯度辨识方法、递阶最小二乘辨识方法、递阶牛顿辨识方法、耦合随机梯度辨识方法、耦合最小二乘辨识方法、耦合牛顿辨识方法等。

梯度方法、最小二乘方法、牛顿方法与滤波辨识理念相结合,诞生出滤波梯度辨识方法、滤波最小二乘辨识方法、滤波牛顿辨识方法等。《青岛科技大学学报(自然科学版)》的连载论文研究了线性回归系统、方程误差系统和输出误差系统的递阶递推辨识方法和递阶迭代辨识方法[8-16]。

最近的连载论文研究了有限脉冲响应滑动平均系统 的递阶递推增广参数辨识方法[9]和递阶增广迭代参数辨识方法[10],以及滤波增广递推参数辨识方法[17]和滤波增广迭代辨识方法[18]。本研究利用滤波辨识理念,研究方程误差自回归系统,即受控自回归自回归系统的滤波递推广义辨识方法。

基于数据滤波的辨识方法一般是针对有色噪声干扰的系统而言。其基本思想是,在系统方程两边同乘以一个滤波器(通常取为噪声模型传递函数的逆),对输入输出数据进行滤波,滤波后的系统是一个白噪声干扰的系统,而原系统的传递函数不变,然后使用滤波后的输入输出数据等进行辨识。因为滤波器(噪声模型的传递函数)是未知的,滤波后的输入输出数据也是未知的,实际中通常采用噪声模型传递函数的估计进行滤波,所以基于数据滤波的辨识方法只能以递推方式或迭代方式实现。

如果知道系统先验知识,滤波器也可以采用参数已知的固定模型,否则可采用动态变化模型,如即将讨论的用噪声传递函数的估计作为滤波器。

辨识中数据滤波与通信和信号处理中的滤波有本质不同。后者滤波是剔除被污染信号中的噪声,使得噪声在信号中的比重降低,如低通滤波器和高通滤波器。而滤波辨识方法只是改变系统模型的结构,不改变系统的输入-输出关系。总之,滤波辨识是将有色噪声干扰的系统模型变换为白噪声干扰的模型,即模型白色化(不是数据白色化),它能够提高参数辨识精度。

1 系统描述与滤波辨识模型

方程误差自回归模型(EEAR模型)又称为受控自回归自回归模型。考虑下列受控自回归自回归模型(controlled autoregressive autoregressive model,CARAR模型)描述的动态随机系统:

其中{u(t)}和{y(t)}分别是系统的输入和输出序列,{v(t)}是零均值、方差为σ2的随机白噪声序列,A(z),B(z)和C(z)是单位后移算子z－1的多项式:

设阶次na,nb和nc已知。记n:＝na＋nb＋nc。且设t≤0时,所有变量的初值均为零,如y(t)＝0,u(t)＝0,v(t)＝0。

置有关参数向量和信息向量如下:

定义中间变量:

取滤波器L(z):＝C(z)。定义滤波输入uf(t),滤波输出yf(t)分别为

定义滤波输出信息向量φaf(t)和滤波输入信息向量φbf(t)如下:

式(1)两边乘以L(z)得到

或

这个滤波后的模型是一个方程误差模型(即CAR模型)。将它写为向量形式为

将式(6)代入式(10)得到

其中

在辨识模型(11)中,参数向量θ包含了系统的所有参数ai,bi和ci。目标就是利用系统的输入输出数据u(t)和y(t),研究和提出估计参数向量θ的辨识方法。

对于辨识模型(11),由于多项式C(z)是未知的(即L(z)是未知的),故uf(t),yf(t),以及φaf(t)和φbf(t)都是未知的,且由它们构成的信息向量φ(t)也是未知的,故需要采用估计值代替未知变量来推导基于数据滤波的辨识方法。

2 滤波广义随机梯度辨识算法

对于EEAR系统的辨识模型(11),定义梯度准则函数:

这些递推关系式无法实现,因为式(13)～(14)右边信息向量φ(t)是未知的,所以需要先构造它们的估计。由噪声模型参数估计

构造多项式C(z)的估计:

根据定义式(7)～(8)的结构,用yf(t－i)和uf(ti)的估计和构造滤波输出信息向量φaf(t)和滤波输入信息向量φbf(t)的估计:

将式(13)～(14)中右边未知信息向量φ(t)用其估计代替,得到式(20)～(22),联立式(19),(17)～(18),(2)～(3),(15)～(16),可以得到估计EEAR系统参数向量θ的基于滤波的广义随机梯度辨识算法(filtering-based generalized stochastic gradient identification algorithm),简称为滤波广义随机梯度算法(filtered generalized stochastic gradient algorithm,F-GSG算法):

F-GSG算法(20)～(30)计算参数估计向量^θ(t)的步骤如下。

1)初始化:令t＝1。置初值,i)＝1/p0,i＝1,2,…,max[na,nb,nc],p0是一个大正数,如p0＝106。给定数据长度Le。

2)采集输入输出数据u(t)和y(t),用式(24)～(27)构造信息向量和φc(t),用式(23)构造信息向量。

3)用式(22)计算r(t),用式(21)计算新息e(t)。

6)如果t＜Le,t就增加1,转到步骤2);否则输出参数估计,终止递推计算过程。

在F-GSG辨识算法(20)～(30)中,e(t)＝为辨识新息。定义辨识残差。

3 滤波多新息广义随机梯度算法

为改进F-GSG算法的收敛速度,一种方法是借助于多新息辨识理论,通过扩展新息维数,推导出基于滤波的多新息广义随机梯度辨识算法。

设正整数p表示新息长度。基于F-GSG算法(20)～(30),将系统输出y(t)和信息向量扩展为堆积输出向量Y(p,t)和堆积信息矩阵(p,t):

将式(20)中标量新息e(t)＝y(t)－^φT(t)^θ(t－1)∈ℝ扩展为新息向量

这是新息长度p＝1的基于滤波的多新息广义随机梯度算法。将上式中和E(1,t)中的“1”换为p,得到式(34),联立式(31)～(33)和(22)～(30),就得到新息长度为p的估计EEAR系统参数向量θ的基于滤波的多新息广义随机梯度辨识算法(filtering-based multi-innovation generalized stochastic gradient identification algorithm),简称为滤波多新息广义随机梯度算法(filtered multi-innovation generalized stochastic gradient algorithm,FMI-GSG算法):

1)初始化:令t＝1,给定新息长度p。置初值,p0是一个大正数,如p0＝106。给定数据长度Le。

2)采集输入输出数据u(t)和y(t),用式(40)～(43)构造信息向量和φc(t),用式(39)构造信息向量。

3)用式(37)和(38)构造堆积输出向量Y(p,t)和堆积信息矩阵。

4)用式(35)计算新息向量E(p,t),用式(36)计算r(t)。

6)如果t＜Le,t就增加1,转到步骤2);否则输出参数估计,终止递推计算过程。

引理1对于F-MI-GSG辨识算法(34)～(46),新息向量E(p,t)与残差向量

满足关系:

4 滤波多新息广义投影辨识方法

下面利用最速下降法推导滤波多新息广义投影辨识方法。考虑从j＝t－p＋1到j＝t的数据窗里的p组数据。定义堆积输出向量Y(p,t)和堆积信息矩阵Φ(p,t)如下:

基于EEAR系统的滤波辨识模型(11),定义滑动数据窗准则函数:

假设步长(step size)为μ(t)。使用最速下降法(负梯度搜索),极小化准则函数J2(θ),可以得到梯度递推关系:

其中

下面求最佳步长(best step-size)μ(t)。将代入准则函数J2(θ)中,可得

由此可求得最佳步长为

由于堆积信息矩阵Φ(p,t)是未知的,所以算法(48)～(51)无法实现。解决的方法是用φ(t－i)的估计构造堆积信息矩阵Φ(p,t)的估计

式(49)～(51)中未知Φ(p,t)用其估计代替,得到式(53)～(55),联立式(48),(52)和(39)～(46),就得到辨识EEAR系统参数向量θ的基于滤波的多新息广义投影辨识算法(filtering-based multi-innovation generalized projection identification algorithm),简称为滤波多新息广义投影算法(filtered multi-innovation generalized projection algorithm,F-MI-GProj算法):

因为收敛因子的计算式比较复杂,故对其进行简化。由于对于任意实向量x和非负定对称矩阵Q,不等式成立,其中为矩阵Q的最大特征值,所以

于是,收敛因子可以保守取为

因为计算矩阵的迹(trace)比计算特征值简单,所以收敛因子可以更保守取为

因此,将式(53)修改为

或

则式(54)和(56)～(67)构成了简化的滤波多新息广义投影算法(F-MI-GProj算法)。

5 滤波递推广义梯度辨识算法

定义堆积输出向量Y(t)和堆积输入信息矩阵Φ(t)如下:

根据式(11),定义二次准则函数:

定义递推关系:

递推关系式(68)～(71)不可实现,因为右边包含了未知信息向量φ(t)。解决方案是用其估计代替,得到式(72)～(75),联立式(23)～(30),便得到辨识EEAR系统参数向量θ的基于滤波的递推广义梯度辨识算法(filtering-based recursive generalized gradient identification algorithm),简称为滤波递推广义梯度算法(filtered recursive generalized gradient algorithm,F-RGG算法):

1)初始化:令t＝1,置初值,,…,max[na,nb,nc],p0是一个大正数,如p0＝106。给定参数估计精度ε。

2)采集输入输出数据u(t)和y(t),用式(77)～(80)构造信息向量和φc(t),用式(76)构造信息向量。

3)用式(73)计算r(t),用式(74)计算向量ξ(t),用式(75)计算矩阵R(t)。用式(72)刷新参数估计向量)。

6 滤波多新息递推广义梯度算法

设正整数p表示新息长度。定义堆积输出向量Y(p,t)和堆积输入信息矩阵Φ(p,t)如下:

基于F-RGG算法(72)～(83),将式(73)～(75)中系统输出y(t)和输入信息向量扩展为堆积输出向量Y(p,t)和堆积信息矩阵),得到式(87)～(89),联立式(72),(84)～(85)和(76)～(83),便得到辨识EEAR系统参数向量θ的基于滤波的多新息递推广义梯度辨识算法(filtering-based multi-innovation recursive generalized gradient identification algorithm),简称为滤波多新息递推广义梯度算法(filtered multi-innovation recursive generalized gradient algorithm,F-MI-RGG算法):

1)初始化:令t＝1, 给定新息长度p和精度指标ε。置初值,,p0是一个大正数,如p0＝106。

2)采集输入输出数据u(t)和y(t),用式(93)～(96)构造信息向量和φc(t),用式(92)构造信息向量。

3)用式(90)和(91)构造堆积输出向量Y(p,t)和堆积信息矩阵。

4)用式(87)计算r(t),用式(88)计算向量ξ(t),用式(89)计算矩阵R(t)。用式(86)刷新参数估计向量。

7 滤波递推广义最小二乘算法

参照递推最小二乘算法的推导,极小化J3(θ),可以得到下列最小二乘递推关系:

这些递推关系无法实现,因为式(100)～(102)右边的信息向量φ(t)是未知的,在辨识算法中使用它的估计代替,得到式(103)～(106),联立式(23)～(30),便得到辨识EEAR系统参数向量θ的基于滤波的递推广义最小二乘辨识算法(filteringbased recursive generalized least squares identification algorithm),简称为滤波递推广义最小二乘算法(filtered recursive generalized least squares algorithm,F-RGLS算法):

1)初始化:令t＝1。置初值,,,p0是一个大正数,如p0＝106。给定参数估计精度ε。

2)采集输入输出数据u(t)和y(t),用式(108)～(111)构造信息向量和φc(t),用式(107)构造信息向量。

3)用式(104)计算新息e(t),用式(105)和(106)计算增益向量L(t)和协方差阵P(t)。根据式(103)刷新参数估计向量。

8 滤波多新息广义最小二乘算法

基于F-RGLS算法(103)～(114),将系统输出y(t)和信息向量分别扩展为堆积输出向量Y(p,t)和堆积信息矩阵:

得到式(118)～(121),联立式(115)～(117)和(107)～(114),便得到辨识EEAR系统参数向量θ的基于滤波的多新息递推广义最小二乘辨识算法(filtering-based multi-innovation recursive generalized least squares identification algorithm),简称为滤波多新息广义最小二乘算法(filtered multi-innovation generalized least squares algorithm,F-MIGLS算法):

F-MI-GLS算法(118)～(131)计算参数估计向量的步骤如下。

1)初始化:令t＝1,给定新息长度p和精度指标ε。置初值,…,max[na,nb,nc],p0是一个大正数,如p0＝106。

2)采集输入输出数据u(t)和y(t),用式(125)～(128)构造信息向量和φc(t),用式(124)构造信息向量。

3)用式(122)和(123)构造堆积输出向量Y(p,t)和堆积信息矩阵。

4)用式(119)计算新息向量E(p,t),用式(120)和(121)计算增益矩阵L(t)和协方差阵P(t)。根据式(118)刷新参数估计向量。

多新息辨识方法包含多新息递推辨识方法和多新息迭代辨识方法。对于多新息递推辨识方法,常常省略递推两个字,将多新息递推最小二乘辨识方法简称为多新息最小二乘辨识方法。如,滤波多新息递推广义最小二乘辨识算法简称为滤波多新息广义最小二乘算法,滤波多新息递推广义最小二乘辨识算法简称为滤波多新息广义最小二乘算法。

9 结 语

针对方程误差自回归(EEAR)系统,研究和提出了滤波(多新息)广义随机梯度辨识方法、滤波多新息广义投影辨识方法、滤波(多新息)递推广义梯度辨识方法、滤波(多新息)广义最小二乘辨识方法。尽管这些滤波广义辨识方法是针对自回归噪声干扰下的方程误差自回归随机系统提出的,但是其思想可以推广到有色噪声干扰下的线性和非线性多变量随机系统中[19-31]。