北京工商大学系统建模与辨识课程上机实验报告(2016 年秋季学期)专业名称:控制工程上机题目:用极大似然法进行参数估计专业班级 :计研3班学生姓名:王瑶 吴超学 号：1001131625910011316260指导教师:刘翠玲2017 年 1 月一 实验目的通过实验掌握极大似然法在系统参数辨识中的原理和应用。二 实验原理1 极大似然原理设有离散随机过程V与未知参数0有关，假定已知概率分布密度f (Vkp )。如果我们得到n个独立的观测 值v ,v,v,贝y可得分布密度f (V 0)，f (V p)，f (V p)。要求根据这些观测值来估计未知参数0，估 12,n12n计的准则是观测值V 的出现概率为最大。为此，定义一个似然函数k1.1)L(V,V，,V 0) = f (V0)f (V 0)/(V 0)12n12n上式的右边是n个概率密度函数的连乘，似然函数L是0的函数。如果L达到极大值，W的出现概率为最大。因此，极大似然法的实质就是求出使L达到极大值的的估值。为了便于求*，对式(1.1)等号两边取对数, 则把连乘变成连加，即InL = Inf (V )(1.2)ii=1由于对数函数是单调递增函数，当L取极大值时，lnL也同时取极大值。求式(1.2)对的偏导数，令偏 导数为 0，可得d In L1.3)=0 d解上式可得的极大似然估计 ml。2 系统参数的极大似然估计Newton-Raphson 法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每 L 次观测数据递推一 次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数估计值得算法。本质上说，它只是一种近似的 极大似然法。设系统的差分方程为a(z-1)y(k)二 b(z-1)u(k) + g (k)(2.1)式中 因为g (k)|是相关随机向量，故(2.1)可写成a(z-1)y(k) = b(z-1)u(k) + c(z-1)e (k)(2.2)式中c (z-1) (k) = g (k)(2.3)c(z-1) = 1 + c z-1 + + c z-n(2.4)1n (k)是均值为0的高斯分布白噪声序列。多项式a(z-1)，b(z-1)和c(z-1)中的系数a ,a,b，b ,c，c和序1,0n 1n列 (k)的均方差b都是未知参数。设待估参数= a a b b c -c 】T(2.5)1 n 0 n 1 n并设y (k)的预测值为c e(k-1) + + c e(k-n)(2.6)1n式中e(k -i)为预测误差；ai，b，为a，b，c的估值。预测误差可表示为iiiiii(c z-1 + cz-2 + + c z-n)e(k)(2.7)12n或者(1 + c z -1 + c” z -n )e(k) = (1 + a z -1 + a” z -n) y (k)-(b0+b1 z-1 + +b z-n)u(k)(2.8)01n因此预测误差越(k)1满足关系式c( z -1)e(k) = a (z-1) y (k) 一 b( z -1)u (k)(2.9)式中假定预测误差e(k)服从均值为0的高斯分布，并设序列4(k)具有相同的方差Q2。因为4(k)与c(z-1),a(z-1)和b(z-1)有关，所以c2是被估参数9的函数。为了书写方便，把式(2.9)写成2.10)2.11)c( z-1)e(k) = a (z-1) y (k) 一 b( z-1)u (k)b u(k 一 n) c e(k 一 1)c (k 一 n), k = n +1, n + 2,n1或写成e(k) = y(k) + 工a y(k - i) -工bu(k -i) -工ce(k -i) iiii=1i=0/ =1令k=n+l,n+2,n+N,可得e(k)的N个方程式，把这N个方程式写成向量-矩阵形式e 二 Y 一 9 2.12)2.13)式中Y=Ny (n +1) y (n + 2),e=e(n +1)e(n + 2)y (n + )e(n + )9anb0b因为已假定(k，是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为11exp (y m)2丄2c 22.14)2.15)(2 兀c 2)2式中y为观测值，c 2和m为y的方差和均值，那么11 exp 一 e-(k)丄2c 2(2 兀c 2)2对于e(k)符合高斯噪声序列的极大似然函数为2.16)1L(Y 9 ,c) =expNN(2兀c 2)2(Y 9 )t(1 9)2g 22.17)一N心-Nlnc2丄eTe222c 2 n N2.18)对上式(2.17)等号两边取对数得1 1lnL(Y 9,c) = ln+ lnexp(一 ereN2c 2 (2 兀c 2)2或写为,nlnL(Y 9,c)二ln2兀 N2求lnL(Y p,c)对c2的偏导数，令其等于0,可得 d ln L(Y 9, c) N=dc 2訓 c 2 一 -L-刘e 2(k)k =n+12c2 + -k 刘e-(k) = k = n +12.19)式中2 二 N 艺e 2(k)二 N 2 艺e 2(k)二 Nk=n+lk=n+l2.21)J 二艺e2(k)2 k n+12越小越好，因为当方差Q2最小时，e2(k)最小，即残差最小。因此希望Q2的估值取最小2.22) min J2.23)因为式(2.10)可理解为预测模型，而e(k)可看做预测误差。因此使式(2.22)最小就是使误差的平方之 和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。因此可按J最小来求a , a, b，b , c,c1,0n 1n的估计值。由于e(k)式参数件，a，叮弋,V的线性函数，因此J是这些参数的二次型函数。求使ln L(Yn P，6最大的$，等价于在式(2.10)的约束条件下求使J为最小。由于J对c.是非线性的，因而求J的极小值问题并i不好解，只能用迭代方法求解。求 J 极小值的常用迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿- 拉卜森法。整个迭代计算步骤如下：确定初始的0 值。对于0 中的j a，件，bn可按模型e(k) a( z-i) y (k) - b( z-)u (k)2.24)用最小二乘法来求，而对于0 0中的C1，Cn可先假定一些值。(2) 计算预测误差e(k) y(k)- y(k)给出并计算 2 e 2(k)Nk n +1J d 2 J(3) 计算J的梯度60和海赛矩阵,有30 22.25)2.26)3J30 e(k)k n +13e(k)302.27)3e(k -n) c 一n3ai式中3e(k-1)3e(k-2)y(k - i) - c- c -13a23aii同理可得沁=y (k - i)上 c jdajdaijTi沁 一u(ki)丄c de(巳dbj dbij Ti2.29)2.30)de( k)dei_e(k i)上e de(k - j)j dcj=1i2.31)将式（2.29）移项化简，有2.32)de(k) * 工de(k - j) _ 牙de(k - j)dajdaj daij=1ij=0i因为e(k - j) = e(k)z- j2.33)由e（k - j）求偏导，故de(k j) _ de(k) z- jii2.34)将（2.34）代入（2.32），所以y(k-i) = Ze Tj)=e de(k)z-j _ de(k)瓦j=0c z- jjj=02.35)所以得(1A de(k) e(z-1)_ y(k -1)dai 同理可得(2.30)和(2.31)为e(z-1)空2 _-u(k - i) dbi2.36)2.37)e( z -1)de(k)de_ -e(k - i)2.38)根据（2.36）构造公式e (z -1)j_ yk - (i - j) - j _ y(k - i)2.39)消除c(z-1)可得同理可得(2.37)和(2.38)式c( z-i)a=c( z j)de(k)da(2.40)de(k)de(k 一 i + j)de(k i +1)( 2.41)daidajda1de( k)de( k 一 i + j)de(k 一i)( 2.42)dbidbjdb0de(k)de(k 一 i + j)de(k i +1)( 2.43)dcidcjdc1式(2.29)、式(2.30)和式(2.31)均为差分方程，这些差分方程的初始条件为 0，可通过求解这些差分方程，分别求出e(k)关于a ,a,b，b ,c，c的全部偏导数，而这些偏导数分别为y(k), u(k)和e(k)1,0n 1n的线性函数。下面求关于e的二阶偏导数，即d 2 Jde2艺k=n+1de(k) r de(k) de L de J +忧 e(k)k=n+1d2e(k)de22.44)nd 2 J当e接近于真值e时，e(k)接近于0。在这种情况下，式(2.44)等号右边第2项接近于0，可近似表de 2示为d2J yn de(k)1 de(k) lT de 2de de _|k n +1( 2.45)则利用式(245)计算迹比较简单。(4)按牛顿-拉卜森计算e的新估值$ 1，有en1=en(巴)-1 -(de 2)dJdeen2.46)重复(2)