《用极大似然法进行参数估计》由会员分享，可在线阅读，更多相关《用极大似然法进行参数估计（21页珍藏版）》请在金锄头文库上搜索。
1、北京工商大学系统辨识课程上机实验报告（2014 年秋季学期）专业名称 ：控制工程上机题目:极大似然法进行参数估计专业班级:2015 年 1 月实验目的通过实验掌握极大似然法在系统参数辨识中的原理和应用。二 实验原理1 极大似然原理设有离散随机过程V 与未知参数有关，假定已知概率分布密度f (V )。如果我们 kk得到n个独立的观测值V ,V,V，则可得分布密度f (V )，f (V )，f (V )。12,n12n要求根据这些观测值来估计未知参数，估计的准则是观测值V 的出现概率为最大。k 为此，定义一个似然函数L(V,V，,V ) = f(V)f (V ”f(V )门“12n12n(1.1)
2、上式的右边是n个概率密度函数的连乘，似然函数L是的函数。如果L达到极大值，V k的出现概率为最大。因此，极大似然法的实质就是求出使L达到极大值的的估值介。为了 便于求 ，对式(1.1)等号两边取对数，则把连乘变成连加，即1.2)In L = In f (V )ii =1由于对数函数是单调递增函数，当L取极大值时，lnL也同时取极大值。求式(1.2) 对 的偏导数，令偏导数为0，可得d ln L = 01.3)解上式可得的极大似然估计ML 。dQ2 系统参数的极大似然估计New to n-Raphson法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每L 次观测数据递推一次的算法，现在
3、我们讨论的是每观测一次数据就递推计算一次参数估计值 得算法。本质上说，它只是一种近似的极大似然法。设系统的差分方程为a(z-1)y(k) = b(z-1)u(k) + g (k)(2.1)式中a(z-1)=1+a z-1 + .+ a z-n 1n b(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)8 (k) = (k)(2.3)c(z-1) = 1 + c z-1 + + c z-n(2.4)1n (k)是均值为0的高斯分布白噪声序列。
4、多项式 a(z-1)，b(z-1)和c(z-1)中的系数a ,a,b，b ,c，c和序列s (k)的均方差a都是未知参数。 1,0 n 1 n设待估参数9= a a b b c c 】T1 n 0 n 1 n 并设y (k)的预测值为y(k)= -a y(k -1)-a y(k -n)+b u(k)+b u(k -n)+ 1nc e(k -1) + c e(k - n)1n式中e(k -i)为预测误差；a.，b.，c.为a，b, iiiiic 的估值。预测误差可表示为 i工a y(k i) + 工bu(k i) +ii i=1i=0e(k)二 y(k) - y(k)二 y(k)n z -n)
5、y(k) - (b 0+哲 z-1 + + b ” z -n )u (k)-工 ce(k - i) = (1+a z-1 + ai1i=1(c z-1 + c z-2 + b c z-n )e(k)12n或者(1 + z-i + C” z-n)e(k) = (1 +z-1 + a” z-n)y(k)-(b + b z-1 + + b z-n )u(k) 01n因此预测误差(k )】满足关系式c (z j)e(k) = a (z -1) y (k) - b (z -1)u (k)2.5)2.6)2.7)2.8)2.9)式中a(z-1) = 1 + a z-1 + a z-n1nb(z-1) =
6、b + b z-1 + b z-n01nc(z-1) = 1 + c z-1 + c z-n假定预测误差e(k)服从均值为0的高斯分布，并设序列e(k)具有相同的方差a2。因 为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) + a y(k-n)-b u(k-1)-bu(k-1)1n01b u(k n) c e(k 1) c (k n), k = n +1, n + 2, (2.11)n1n或写成2
7、.12)e(k) = y(k) + 工a y(k -i) -工bu(k -i) -工ii i=1i =0令k=n+l,n+2,n+N,可得e(k)的N个方程式，c e(k - i )ii=1把这N个方程式写成向量-矩阵形式e = Y 一 9N N NY=Ny (n +1) y (n + 2)，eN =e(n +1)e(n + 2)y (n + N)e(n + N)9 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)u (N)e(n + N 1) e(N)y(n + N 1)一
8、y(N)因为已假定 e(k) 是均值为0 的高斯噪声序列，高斯噪声序列的概率密度函数为f =-exp (y m)212c 2(2kc 2)22.14)式中y为观测值，c 2和m为y的方差和均值，那么exp1 (2 兀c 2)212c 2e 2(k)2.15)对于e(k)符合高斯噪声序列的极大似然函数为L(Yn 9, c) = Le(n +1), e(n + 2),e(n + N) 9 = f e(n +1) 9 f e(n + 2) 9 f e(n + N) 9 1expe2(n +1) + e2(n + 2) Hb e2 (n + N) =1exp( 1 N2c 2N(2 兀c 2)2(2
9、兀c 2)2eTe )2c 2 N N2.16)(Y 9 )t (Y 9 人 氓L(Y 9 ,c) =exp nNN2c 2(2 兀c 2)2对上式(2.17)等号两边取对数得ln L(Y 9 ,c) = ln1+ ln exp( eTeNN2c 2 N .(2 兀c 2)22或写为InL(Y 9,c) = Nln2兀-Nn2求ln L(Yn9 ,c)对c 2的偏导数，a ln L(Y 9 ,c)N=ac2Nc 1ln c 2 22c 2k=n+1令其等于 0，可得N 1 n Ne2(k) = 0+ 2c 22c 4k=n+l2.17)eTe2c 2 N N2.18)2.19)2.21)人2=
10、N 迅e 2(k)=N 2 迅e 2(k)=NJ式中k n+1k n+12.22)J =丈 n e2(k)2kn+1A b2N2越小越好，因为当方差Q2最小时，e2(k)最小，即残差最小。因此希望b 2的估值取最 小2.23)因为式(2.10)可理解为预测模型，而e(k)可看做预测误差。因此使式(2.22)最小 就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。因 此可按J最小来求a ,a,b，b ,c，c的估计值。1,0n 1n由于e(k)式参数a ,a,b，b ,c，c的线性函数，因此J是这些参数的二次型函数。1,0n 1n求使InL(YnQ)最大的$，等价于在
11、式(2.10)的约束条件下求$使J为最小。由于J对c.是非线性的，因而求J的极小值问题并不好解，只能用迭代方法求解。求J极小值的常用 i迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤 如下：确定初始的$0值。对于$0中的Ja耳，bn可按模型e(k) a (z -1) y (k) - b( z -1)u (k)(2.24)用最小二乘法来求，而对于$ 0中的ci Cn可先假定一些值。(2) 计算预测误差2.25)2.26)2.27)e(k)= y(k)-y(k)给出J 艺 e2(k)2k=n+1并计算b 2 刃 e 2(k)Nk=n+1dJ6 2 J(3) 计算
12、J的梯度凤 和海赛矩阵，有6$ 26J6$k = n +1de(k)de(k) de(k) de( k)de(k) de(k) de(k)dadadbdbdcdc1n0n1nde(k) _dadaiidy(k) + a y(k 1) + a y(k n) b u(k) bu(k 1)b u(k n)01ce(k 1)c e(k n)1_ y(k i) cde(k 1)de(k 2)cda2idaide(k 一 n)dai2.28)de(k)daiy(k-i) Yc de(k j)jj_1dai2.29)同理可得讐 _u(k i)丄cjij _1de(k j)dbi2.30)沁 _e(k i)
13、Yc de(k j) dcj dcij _1i将式(2.29)移项化简，有2.31)因为y(k i) _ 凹 + Yc de(k j) _Yc de(k j) dajdajdaij _1ij _0i2.32)e(k j)二 e(k) z- j2.33)由e（k j）求偏导，故de(k 一 j) _ de(k)zjdai将（2.34）代入（2.32），所以dai2.34)c z - j daji j _ 02.35)y(k -i) _Yc de(k j) _Yc de(k)z一j _沁工 jj_0c(z-i) _ 1 + c z1 + + c zn1n所以得(八 de(k)c(z-1)_ y(k i)dai同理可得（2.30）和（2.31）为c(z i)i2.37)c( z-1) e(k) =
