广义最小二乘法

1、4.5 广义最小二乘法(GLS)GLSGeneralized Least Squares1 .基本原理广义最小二乘法的基本思想在于引入一个所谓成形滤波器(白化滤波器)，把相关噪声(k)转化成白噪声(k)。由方程(4-4)、(4-5),系统的差分方程可以表示为,11、a(z )y(k) b(z )u(k)(k)(4-114)式中 112na(z ) 1 a1z a2zanzb(z 1) b0 b1z 1b2z 2bnz n如果知道有色噪声序列(k)的相关性，则可以把 (k)看成白噪声通过线性系统后 所得的结果。这种线性系统通常称为成形滤波器，其差分方程为c(z 1) (k) d (z 1) (k

2、)(4-115)式中(k)是均值为零的白噪声序列，C(z 1)、d(z 1)是z 1的多项式。令 f(z1) C(z )1 f1z 1 f2z2 L L fmz m(4-116)d(z1)-11.、有 f(z ) (k)(k)或(k) 1- (k)(4-117)f (z )即(1 Lz1 f2z2 L Lfmzm) (k)(k)(4-118)或(k)f1 (k 1) f2 (k 2)fm (k m) (k) k n 1,L L ,n N(4119)这一噪声模型(自回归模型)的阶m, 一般事先是不知道的，实际经验表明，若指定为2或3,就可以获得令人满意的描述(k)的模型。把方程(4-119)看作

3、输入为零的差分方程，并由此式来写出N个方程。(n 1)fi (n) f2 (n 1)(n 2)G (n 1) f2 (n)fm (n 1 m)(n 1)fm (n 2 m)(n 2)(n N)f1 (n N写成向量矩阵形式为f(n 1)其中，f(n N)(n)(n 1)(n N 1)1) f2 (n N 2)f1(n 1)fm(n N)(n 1)(n)(n N 2)fm (n N m) (n N)(4-120)(n 1 m)(n 2 m)(n N m)，称作最小二乘格式。应用我(4-120)式所示的线性组合关系是辨识问题的基本表达形式们熟知的最小二乘法，可求出 f的估值。T 1 T? (4-1

4、21)1考虑在噪声模型(4-117)下，系统的差分方程。将(4-117)式得到的(k) (k), f(z )代入(4-114)式有111a(z )y(k) b(z )u(k)(k)(4f(z )-122)变换成a(z 1) f (z 1 )y(k) b(z 1) f (z 1)u(k)(k)(4-123)令f(z 1)y(k) y(k)(4-124)f(z1)u(k) u(k)(4-125)于是可得a(z1)y(k) b(z 1)u(k)(k)(4-126)y(k)aiy(k1)a2y(k2) Lay(kn)b0U(k)DU(k1) LbU(kn) (k)(4-127)在(4-126)或(4-

5、127)式中，(k)为不相关的随机序列(白噪声)，故可以用最小二乘法得到的无偏估计(即a1, an,bo,b1bn )。由此可见，广义最小二乘法(GLS)是建【3】1立在最小二乘法(LS)的基础之上的。基本最小二乘法只是广义最小二乘法在f(z1) 1时的特例。2 .计算步骤：广义最小二乘法的关键问题是如何用比较简便的方法找到成形滤波器的系数。其计算是逐次逼近法。下面以(4-121)和(4-123)式为依据来讨论。第一步：应用输入、输出数据 u(k)和y(k)(k 1,2 n N),按最初的模型11a y k b u k k求出的最小二乘估计?国1皆咛?1考an1)玲1)*1)T这个估值是不精确

6、的，它只是被估参数的一次近似。第二步：计算残差e(1)(),并拟合成形滤波器的模型：e()？(1)(1)y(k) b?(1)(1)u(k)(4-128)或e(1) ( ) y( ) aij1) y(k 1)sin1)y(k n) b(1)u(k) l?1)u(k 1)l?(1)u(k n)k n 1,n N(4-129)e(k)称为广义残差，用广义残差ek1)代替 (k)注意：残差e(k)是模型噪声，而 k 是系统噪声，可用最小乘法拟合出一个成型滤波器的模型。e(k)f1e(1)(k 1)fme(1)(k m) (k)利用(4-121)式得到f的最小二乘估值?(1)T(1)1(1)Tc(1)

7、e(4-131 )其中f?？L le(1)e(1) (n 1)e(1)(n2)e(n N)e(n)Me(1)(n N 1)e(1)(n m 1)Me(1)(n m N)第三步：应用而所得的成形滤波器，对输入输出数据滤波由(4-124) (4-125)y(1) (k) F(1)u(1) (k) F1)y(k)1)u(k)(4-132)y(k) y(k) U(k) u(k)f?(1) y(k 1)?(1)u(k1)f?(1)y(kd1%m)(4-133)m)第四步：求出参数的第二次估值？按模型重新估计a(z-1)y(1)(k)b(z 1)u (1) (k)(k),得其最小=乘估值？(2),并转入第

8、二步。以下将重复以上步骤：由72),按步骤2计算残差e(2) k按步骤3计算y(k)、U(k);按步骤求的第3次估值？。，估值?(2);进行多次循环，直到的第i次估值？收敛为止。循环程序的收敛性判断:iHm?(i)(z1) 1(4-134)这意味着残差e(k)已经白噪声化，数据不需要继续滤波的一个良好估计。广义最小二乘辩识算法的程序框图:输入 u(k).y(k) n N组图4.12广义最小二乘辩识算法的程序框图(k)(k)3 .算例：将 GLS算法与LS算法进行比较有一单输入一单输出系统y(k)a1y(k 1) a2y(k 2) b1u(k 1)的真值为Ta a2 b 卜0.5 0.5 1.0

9、输入u (k)是具有零均值和单位方差的独立高斯随机变量序列；(k)的成形滤波器模型为：f (z 1) (k) (1.0 0.85z 1) (k)即f0 1.0, f1 0.85( G较大是为了使残差强烈相关)2(k)是具有方差0.64的零均值白噪声。辩识结果：利用了 N=300的输入输出数据，用广义最小二乘法( GLS)进行迭代计算。每次迭代计算，都计算出残差的均方误差:T2 ef,计算结果绘于图 X中。N4.优缺点：优：能够克服当存在有色噪声干扰时，基本最小二乘估计的有偏性，估计效果较好，在 实际中得到较好的应用。缺： 计算量大，每个循环要调用两次最小二乘法及一次数据滤波， 求差分方程（4-

10、123）的参数估值，是一个非线性最优化问题，不一定总能保证算法对最优解的收敛性。广义最小二乘法本质上是一种逐次逼近法。对于循环程序的收敛性还没有给出证明。 GLS算法的最小二乘指标函数 J中可能存在一个以上局部极小值，（特别在信噪比不大时，J可能是多举的）。GLS方法的估计结果往往取决于所选用参数的初始估值。参数 估计初值应选得尽量接近优参数。在没有验前信息的情况下，最小二乘估值被认为是最好的初始条件。 广义最小二乘法的收敛速度不是很高。（图4.11即是一个说明）。4.6递推广义最小二乘法(RGLS )RGLS-Recursive Generalized Least Squares首先回顾整批

11、GLS辨识方法中的有关算式：f(z1) (k)(k)(4-117)a(z 1)f (z 1)y(k) b(z 1) f (z 1)u(k)(k)(4-123)1、y(k)f (z )y(k)(4-124)1、 u (k)f (z )u(k)(4-125)1 _1 _a(z )y(k) b(z )u(k)(k)(4-126)(4-126)式的形式和(4-4)、(4-5)式完全一样，相应有Y -Y(4-33) (4-135)为建立递推计算公式，当观测数据长度为N时，把(4-135)式写成：YNN N N(4-136)广义最小二乘法的递推计算过程可分成二部分：(1)按递推最小二乘法(RLS),随着N

12、的增大，不断计算 ?N (逐步接近于无偏)和?N (逐步使噪声白化)；(2)在递推过程中，%和？N是时变的，则过滤信号 U(k).y(k)及残差e(k)是由时变系统产生，要不断计算 U(k).y(k)及e(k)。因而，递推广义最小二乘法由两组普通的递推最小二乘法组成，它们是通过滤波算法联系起来的：令Y(k) ?(z1)Y(k)U(k) ?(z 1)u(k) y n 1 y(n N 1) "J 1 y(n N) y(n N 1)L L y(N 1)u(n N 1)L L U(N 1) e(k) ?(z1)y(k) b(z 1)u(k)N 1 e(n N) e(n N 1)L L e(n N m 1) eN 1 e(n N 1)为由e(n)到e(n N m)所组成的矩阵块，则有递推公式：? K (f N 1 N KN 1(yN 1(4-138)(4-139)(4-140)P(2)T P(2)Pn(21PN2)PN F.PN(4-143)1 N 1PNN 1(4-141 )(4-142)L PN2) T 1(4-144)RGLS辩识算法的程序框图(见参考文献13之P197)。递