用极大似然法进行参数估计

北京工商大学系统辨识课程上机实验报告（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)上式的右边是n个概率密度函数的连乘，似然函数L是的函数。如果L达到极大值，V k的出现概率为最大。因此，极大似然法的实质就是求出使L达到极大值的的估值介。为了 便于求 ，对式(1.1)等号两边取对数，则把连乘变成连加，即1.2)In L = In f (V )ii =1由于对数函数是单调递增函数，当L取极大值时，InL也同时取极大值。求式(1.2) 对 的偏导数，令偏导数为0，可得d ln L = 01.3)解上式可得的极大似然估计ML 。dQ2 系统参数的极大似然估计New to n-Raphson法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每L 次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数估计值 得算法。本质上说，它只是一种近似的极大似然法。设系统的差分方程为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的高斯分布白噪声序列。多项式 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) 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) = 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.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)一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 兀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=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)最大的$，等价于在式(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) 计算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) 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) = -e( k - i)dci2.38)根据（2.36）构造公式c(z-1)為伙2 理=y伙-(i - j) - j = y(k - i)da j将其代入（2.36），可得2.39)c( z T)dek (i - j) =z -1) de(k)dadaji(2.40)消除c（z-1）可得de（k）dadaij同理可得（2.37）和（2.38）式de(k - i + j) _ de(k - i +1)da12.41)de(k) _ de(k - i + j)dbidbjde(k - i)db02.42)2.43)de(k)de(k - i + j)de(k - i +1)dcdcdcij1式（2.29）、式（2.30）和式（2.31）均为差分方程，这些差分方程的初始条件为 0可通过求解这些差分方程，分别求出e（k）关于a ,a,b，b ,c，c的全部偏导数，而这1,0n 1n些偏导数分别为y（k），u（k）和e（k）的线性函数。下面求关于0的二阶偏导数，即d 2 J50 2k _n+1de(k)de(k)-dbdJ +艺 e(k)k _ n +1d 2e(k)d0 22.44)当0接近于真值0时，e（k）接近于0。在这种情况下，式（2.44）等号右边第2项接近于0，J可近似表示为d0 2k_n+1d2 J _ 戈n de(k)1 de(k) lT d0 2d0 I d0 I则利用式(2.45)计算比较简单。(4) 按牛顿-拉卜森计算0的新估值$ 1，有0aj302.46)重复(2)至(4)的计算步骤，经过r次迭代计算之后可得近一步迭代计算可得0 =0 -r+1r(竺)一亘30 2302.47)如果AAQ 2Q 2r+1 r 10 - 42.48)AQ2r则可停止计算，否则继续迭代计算。式(2.48)表明，当残差方差的计算误差小于0.01% 时就停止计算。这一方法即使在噪声 比较大的情况也能得到较好的估计值0 。三 实验内容设 SISO 系统差分方程为y(k) + a y(k 1) + a y(k 2) = b u(k 1) + b u(k 2) + g (k)1 2 1 2其中极大似然估计法默认e(k)为：e(k)=C(z1)g(k)=g(k)+cg(k1)+c g(k2) 12辨识参数向量为0 = a1a2b1 b2 C c2T式中，g(k)为噪声方差各异的白噪声或有色噪声；u(k)和y(k)表示系统的输入输出变量。四 实验流程图1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890.0.433 -1.177 -0.390 -0.982 1.435 -0.119 -0.769 -0.899 0.882. -1.008 -0.844 0.628 -0.679 1.541 1.375 -0.984 -0.582 1.609.0.090 -0.813 -0.428 -0.848 -0.410 0.048 -1.099 -1.108 0.259.-1.627 -0.528 0.203 1.204 1.691 -1.235 -1.228 -1.267 0.309.0.043 0.043 1.461 1.585 0.552 -0.601 -0.319 0.744 0.829.-1.626 -0.127 -1.578 -0.822 1.469 -0.379 -0.212 0.178 0.493.-0.056 -0.1294 1.228 -1.606 -0.382 -0.229 0.313 -0.161 -0.810.-0.277 0.983 -0.288 0.846 1.325 0.723 0.713 0.643 0.463.0.786 1.161 0.850 -1.349 -0.596 1.512 0.795 -0.713 0.453.-1.604 0.889 -0.938 0.056 0.829 -0.981 -1.232 1.327 -0.681.0.114 -1.135 1.284 -1.201 0.758 0.590 -1.007 0.390 0.836.-1.52 -1.053 -0.083 0.619 0.840 -1.258 -0.354 0.629 -0.242.1.680 -1.236 -0.803 0.537 -1.100 1.417 -1.024 0.671 0.688.-0.123 -0.952 0.232 -0.793 -1.138 1.154 0.206 1.196 1.013.1.518 -0.553 -0.987 0.167 -1.445 0.630 1.255 0.311 -1.726.0.975 1.718 1.360 1.667 1.111 1.018 0.078 -1.665 -0.760.1.184 -0.614 0.994 -0.089 0.947 1.706 -0.395 1.222 -1.351.0.231 1.425 0.114 -0.689 -0.704 1.070 0.262 1.610 1.489.-1.602 0.020 -0.601 -0.020 -0.601 -0.235 1.245 1.226 -0.204.0.926 -1.297 %输入数据Y=0.086 2.210 0.486 -1.812 -3.705 -2.688 1.577 2.883 3.705 .1.642 0.805 -2.088 0.946 -0.039 1.984 -2.545 -1.727 -0.231.2.440 3.583 2.915 1.443 3.598 0.702 2.638 3.611 -0.168.1.732 0.666 2.377 -0.554 -2.088 2.698 0.189 -1.633 -2.010.1.716 -1.641 -1.885 1.061 -0.968 2.911 3.088 -1.629 -1.533.3.030 0.614 -1.483 -1.029 -1.948 -1.066 -0.113 -2.144 -2.626.0.134 -3.043 -1.341 0.338 2.702 3.813 -1.924 -2.813 -1.795.3.002 1.027 1.027 2.755 3.584 1.737 -0.837 -0.617 1.703.2.045 -2.886 -0.542 -2.991 -1.859 3.045 0.068 -0.375 0.451.1.036 0.153 -0.474 2.512 -2.681 -0.954 -0.307 0.628 -0.270.-0.277 0.983 -0.288 0.846 1.325 0.723 1.750 1.401 1.340.0.916 1.396 2.446 2.103 2.432 -1.486 3.031 2.373 -0.763.-0.752 -3.207 1.385 -1.642 -0.118 1.756 -1.613 -1.690 2.136.-1.136 -0.005 -2.210 2.331 -2.204 0.983 1.347 -1.691 0.595.1.809 -2.204 -2.330 -0.454 1.290 2.080 -1.990 -0.770 1.240. -0.252 3.137 -2.379 1.206 1.221 -1.977 2.471 -1.680 1.148.1.816 0.055 -1.856 0.269 -1.323 -2.486 1.958 0.823 2.481.2.209 3.167 -0.762 -2.225 -0.123 -2.786 1.026 2.843 1.071.-3.317 1.514 3.807 3.388 3.683 -1.935 -1.935 0.309 -3.390.-2.124 2.192 -0.855 -1.656 0.016 1.804 3.774 -0.059 2.371.-2.322 -0.032 2.632 0.565 -1.460 -1.839 1.917 0.865 3.180.3.261 -2.755 -0.536 -1.171 -0.905 -3.303 -0.834 2.490 3.039.0.134 1.901%输出数据na=2;nb=1;nc=2;d=1;nn=max(na,nc);L=size(Y,2);xiek=zeros(nc,1); %白噪声估计初值yfk=zeros(nn,1); %yf(k-i)ufk=zeros(nn,1); %uf(k-i)xiefk=zeros(nc,1); %vf(k-i)thetae_1=zeros(na+nb+1+nc,1); %参数估计初值P=eye(na+nb+1+nc);for k=3:L%构造向量phi=-Y(k-1);-Y(k-2);U(k-1);U(k-2);xiek; %组建 h ( k)xie=Y(k)-phi*thetae_1;phif=-yfk(1:na);ufk(d:d+nb);xiefk;%递推极大似然参数估计算法K=P*phif/(1+phif*P*phif);thetae(:,k)=thetae_1+K*xie;P=(eye(na+nb+1+nc)-K*phif)*P;yf=Y(k)-thetae(na+nb+2:na+nb+1+nc,k)*yfk(1:nc); %yf(k)uf=U(k)-thetae(na+nb+2:na+nb+1+nc,k)*ufk(1:nc); %uf(k)xief=xie-thetae(na+nb+2:na+nb+1+nc,k)*xiefk(1:nc); %vf(k)%更新数据thetae_1=thetae(:,k);for i=nc:-1:2xiek(i)=xiek(i-1);xiefk(i)=xiefk(i-1);endxiek(1)=xie;xiefk(1)=xief;for i=nn:-1:2yfk(i)=yfk(i-1);ufk(i)=ufk(i-1);endyfk(1)=yf;ufk(1)=uf;endthetae_1figure(1)plot(1:L,thetae(1:na,:);xlabel(k); ylabel (参数估计 a);legend(a_1,a_2); axis(0 L -2 2);figure(2)plot(1:L,thetae(na+1:na+nb+1,:);xlabel(k); ylabel (参数估计 b);legend(b_1,b_2); axis(0 L -1.5 2);figure(3)plot(1:L,thetae(na+nb+2:na+nb+nc+1,:);xlabel(k); ylabel (参数估计 c);legend(c_1,c_2); axis(0 L -2 2);六 实验结果实验结果如下所示。thetae_1 =-0.12220.03651.79620.00040.0655-0.0018可知，估计值a =-0.1222, a =0.0365, b =1.7962, b = 0.0004, c = 0.0655, c =-0.0018 1 2 1 2 1 2图 1 A 参数波形变化情况卽估数参卽估数参图2 B 参数波形变化情况图 3 C 参数波形变化情况图 10-11 图 10-12 有 162 个数据，结果如下：thetae_1 =-0.70950.04891.8688-0.9850-0.1768-0.0005可知，估计值 a =-0.7095，a =0.0489，b = 1.8688，b = -0.9850，c = -0.1768，c = -0.00051 2 1 2 1 2251卽估数参50500O52 .1-20608060O51-图 1 A 参数波形变化情况20 40 60 80 100 120 140 160 k0粹估数参505051图2 B 参数波形变化情况图 3 C 参数波形变化情况六 实验总结通过此次实验，对极大似然法的原理和方法进行了研究和对其mat lab的仿真，可发现 其离线辨识效果很好，但是其在线辨识的计算量很大，实验数据不同，出现的效果也不一样 递推的极大似然法呈现病态，所得参数估计取值很不稳定，以至于不能采用。