《用极大似然法进行参数估计46242》由会员分享，可在线阅读，更多相关《用极大似然法进行参数估计46242（17页珍藏版）》请在金锄头文库上搜索。
1、北京工商大学系统辨识课程上机实验报告（2014 年秋季学期）专业名称:控制工程上机题目:极大似然法进行参数估计专业班级:2015 年 1 月实验目的通过实验掌握极大似然法在系统参数辨识中的原理和应用。二 实验原理1 极大似然原理设有离散随机过程匕与未知参数有关，假定已知概率分布密度f(vk )。如果我们 得到n个独立的观测值V ,V,V，则可得分布密度f (V )，f (V )，f (V )。12,n12n要求根据这些观测值来估计未知参数，估计的准则是观测值V 的出现概率为最大。k 为此，定义一个似然函数L(V,V，,V ) = f (V)f (V )/(V )12n12n(11)上式的右边是
2、n个概率密度函数的连乘，似然函数L是的函数。如果L达到极大值，V k的出现概率为最大。因此，极大似然法的实质就是求出使L达到极大值的的估值宀。为了便于求 ，对式(1.1)等号两边取对数，则把连乘变成连加，即1.2)In L = In f (V )ii=1由于对数函数是单调递增函数，当L取极大值时，lnL也同时取极大值。求式(1.2) 对 的偏导数，令偏导数为0，可得皿=01.3)解上式可得。的极大似然估计ML 。62 系统参数的极大似然估计New to n-Raphson法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每L 次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递
3、推计算一次参数估计值 得算法。本质上说，它只是一种近似的极大似然法。设系统的差分方程为a(zt)y(k) = b(z-1)u(k) + (k)(2.1)式中a(z-1) = 1 + a z-1 +. + a z-n1nb(z-1) = b + b z-1 + + b z-n01n因为E (k)|是相关随机向量，故(2.1)可写成a(z-1)y(k) = b(z-1)u(k) + c(z-1)e (k)(2.2)式中c (z-1) (k) = (k)(2.3)c(z-1) = 1 + c z-1 + + c z-n(2.4)1n (k)是均值为0的高斯分布白噪声序列。多项式 a(z-1)，b(z
4、-1)和c(z-1)中的系数a ,a,b，b ,c，c和序列& (k)的均方差a都是未知参数。 1,0n 1n设待估参数2.5)9= aa bb cc】T1 n 0 n 1 n 并设y (k)的预测值为y (k) = - a y (k -1)a y (k - n) + bb u (k) hh bb u (k - n) +1nc e(k -1) hh c e(k - n)1n式中e(k -i)为预测误差；a.，b.，c.为a，b, iiiiie(k) = y(k) - y (k) = y (k)-2.6)或者 c e(k - i) = (1 +ii=1c 的估值。预测误差可表示为i工 a y (
5、k i) + 工 b u (k i) +ii=1ii=0a z-i hh a z-n) y (k) - (b + b z-1 hh b z-n )u (k)-1n01n(c z-i + c z-2 + c z-n )e(k)12n2.7)(1 + 勺 z-i hf Cn z-n )e(k) = (1 + ax z-1 hf an z-n) y(k)-(b + b z-i hf b z-n )u (k)01n因此预测误差(k卩满足关系式2.8)式中c (z -i)e(k) = a (z-1) y (k) - b (z -i)u(k)2.9)A /-gAAa(z-1) = 1 + a z-1 hh
6、 a z-n1nb (z -1) = b? + b z -1 hh b z-n01nAAAc( z -1) = 1 + c z -1 + H c z -n假定预测误差e(k)服从均值为0的高斯分布，并设序列b(k)具有相同的方差02。因 为4(k)与c(z-1)，a(z-1)和b(z-1)有关，所以a2是被估参数9的函数。为了书写方便， 把式(2.9)写成c(z-1)e(k) = a(z-1)y(k)-b(z-1)u(k)(2.10)e(k) = y(k) + a y(k-1)hha y(k-n)-b u(k-1)-bu(k-1)1n01b u(k n) c e(k 1) c (k n), k
7、 = n +1, n + 2, (2.11)n1n或写成e(k)= y(k)ha y(k -i)-bu(k-i)-ce(k-i)iiii=1i =0i=1令k=n+l,n+2,n+N,可得e(k)的N个方程式，把这N个方程式写成向量-矩阵形式2.12)式中e = Y 一 9N N N2.13)y(n +1)e(n +1)y(n + 2), eN =e(n + 2)_ y(n + N)_e(n + N)9Y=Nanb0- y (n)-y (n +1) - y(1) - y(2)u (n +1)u (1) u (n + 2)u (2)e(n) e( n +1) e(1) e(2)u (n + N)
8、u (N)e(n + N -1)e(N)-y (n + N -1)一 y (N)因为已假定b(k“是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为11exp (y m)2丄2c 2(2kc 2)22.14)式中y为观测值，c 2和m为y的方差和均值，那么11 expe2(k)丄2c 2(2兀c 2)2对于e(k)符合高斯噪声序列的极大似然函数为2.15)L(Yn 9, c) = Le(n +1), e(n + 2),e(n + N) 9 = f e(n +1) 9 f e(n + 2)p f e(n + N) 9 1exp-e 2(n +1) + e 2(n + 2) Hb e2(n +
9、 N) =1exp( 1 N2c 2N(2 兀c 2)2(2 兀c 2)2eTe )2J 2 NN2.16)1L(Y 9,c) =exp-n一2c 2(2兀c 2)2对上式(2.17)等号两边取对数得11lnL(Y 9,c) = ln+ lnexp(一 ereN一2c 2 N .(2kc 2)2(Y -9 )t (Y -9 人 科2.17)2-1eTe2c 2 N N2.18)或写为“NN1lnL(Y 9,c) = -一ln2兀一一Inc2 - n222c 2求lnL(Yn9,c)对c 2的偏导数，令其等于0,可得lnL(YN9，c)=-上 + 丄龙一e2(k) = 0Qc 22c 2 2c
10、4k=n+1龙e 2(k)k =n+12.19)2.20)式中2=N + e2(k)=N 2+ e2(k)=Nk =n+1k =n+12.21)J =艺 e 2(k)2k = n +1b 2越小越好，因为当方差b2最小时，e 2( k )最小，即残差最小。因此希望b 2的估值取最 小2.22)A 22.7b = mm JN因为式(2.10)可理解为预测模型，而e(k)可看做预测误差。因此使式(2.22)最小 就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。因 此可按J最小来求a ,a,b，b ,c的估计值。1,0n 1n由于e(k)式参数a ,a,b，b ,c，c
11、的线性函数，因此J是这些参数的二次型函数。1,0n 1n2.23)求使InL(YnP)最大的$，等价于在式(2.10)的约束条件下求0使J为最小。由于J对c 是非线性的，因而求 J 的极小值问题并不好解，只能用迭代方法求解。求 J 极小值的常用 i迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤 如下：确定初始的 0值。对于 0中的a1,., a，b0，bn可按模型e(k) = a( z t) y (k) - b( z-i)u(k)2.24)用最小二乘法来求，而对于0中的CCn可先假定一些值。(2)计算预测误差e( k) = y (k) - y (k)2.25
12、)给出二 2 艺e 2(k)k=n+1并计算2=N 岀e 2(k)k = n+12.26)式中Q 2 J(3)计算J的梯度60和海赛矩阵 ，有26J/?+N (k)6e(k)6060k=n+1de(k)=d0de(k)de(k) de(k)de(k) de(k) de(k)da dadbdbdcdc1n0n1nde(k)=dadaiidy(k) + a y(k -1) + a y(k -n)-b u(k)-bu(k -1)b u(k -n)-01c e(k -1)c e(k - n)1二 y(k -i) - c1de(k-1)de(k-2)cda2daiide( k - n)dai2.28)de(k)daiy(k - i)-工c de(k-j)jj=1dai2.29)同理可得讐=-u (k-i)上 cjij =1de(k j)dbi2.30)沁= -e(k-i)-Ycdcjdcij =1i将式(2.29)移项化简，有2.31)因为y (k - i)二型 +c de(k - j)二工c de(k 一 j) daj dajdaij Tij=0i2.32)e(k - j)二 e(k)z-j2.33)由e（k - j）求偏导，故Qe
