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

1、北京工商大学系统辨识课程上机实验报告（2014年秋季学期）专业名称：控制工程上机题目：极大似然法进行参数估计专业班级：2015年1月实验目的通过实验掌握极大似然法在系统参数辨识中的原理和应用。二实验原理1极大似然原理设有离散随机过程Vk与未知参数 有关，假定已知概率分布密度 fM )。如果我们f(Vi ) , f(V2 )，,f(Vn )。要求根据这些观测值来估计未知参数为此，定义一个似然函数L(VlV ,Vn )上式的右边是n个概率密度函数的连乘,，估计的准则是观测值Vk的出现概率为最大。f(Vi )f(V2 ) f(Vn )(1.1 )似然函数L是 的函数。如果L达到极大值，Vk的出现概率

2、为最大。 因此，极大似然法的实质就是求出使L达到极大值的的估值。为了得到n个独立的观测值V1,V2,，Vn，则可得分布密度便于求 ，对式(1.1 )等号两边取对数，则把连乘变成连加，即nln L ln f (Vi )(1.2 )i 1由于对数函数是单调递增函数，当L取极大值时，lnL也同时取极大值。求式(1.2 )对 的偏导数，令偏导数为0，可得0(1.3 )解上式可得的极大似然估计ML。2系统参数的极大似然估计Newton-Raphson法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每 L 次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数估计值 得算法。

3、本质上说，它只是一种近似的极大似然法。设系统的差分方程为a(z 1)y(k) b(z 1)u(k)(k)(2.1 )式中a(z1)11na z an z1b(z ) b0dZ 1.bnZ n因为(k)是相关随机向量，故(2.1 )可写成a(z 1)y(k) b(z1)u(k)c(z 1) (k)(2.2 )式中c(z 1) (k)(k)(2.3 )c(z 1)1 c1z 1nqz(2.4 )(k)是均值为 0的咼斯分布白噪声序列。多项式a(z 1),b(z 1)和c(z 1)中的系数a. ,a,b°, bi,G,设待估参数G和序列 (k)的均方差都是未知参数。(2.5 )a1Tanb

4、0bnGcn并设y(k)的预测值为y(k)a1 y(k 1)an y(k n) b°u(k)bnu(k n)(2.6 )e(k)y(k)y(k)y(k)nai y(ki 1ni)bu(k i)i 0nqe(k i)i 1(1a1 z 1an z n)y(k) (b0b1 z 1bn z n)u(k)(1c1 z2c2 zCn zn)e(k)(2.7 )或者(11C1 zcn z n)e(k)=(1a1 z 1anz n)y(k)(b0b1 z 1bn z n)u(k)(2.8 )式中e(k i)为预测误差；ai , bi , ci为a, b , c的估值。预测误差可表示为因此预测误差

5、 e(k)满足关系式(2.9)G e(k 1)Cn e(k n)c(z 1)e(k)a(z 1)y(k)b(z 1)u(k)a(z 1)1a-i zb(z1)b0bi zc(z 1)1G znan znbnZnCn Z式中为 e(k)与c(z ), a(z )和b(z )有关，所以是被估参数的函数。为了书写方便，把式(2.9)写成c(z 1)e(k) a(z 1)y(k) b(z1)u(k)(2.10 )e(k) y(k) a(k 1)any(kn) b°u(k 1)biu(k 1)bnu(k n) &e(k 1)Cn(kn), k n 1,n2,(2.11 )或写成ne(k

6、)y(k)aiy(k i)nbiu(kni)cie(ki)(2.12 )假定预测误差e(k)服从均值为o的高斯分布，并设序列e(k)具有相同的方差 2。因i 0i 1N个方程式写成向量-矩阵形式i 1令k=n+1,n+2,n+N,可得e(k)的N个方程式，把这(2.13 )eNYnn式中Yny(ny(n2),eNe(n 2)5ab0N)e(n N)bny(n)y(1)u(n1)u(1)e(n)e(1)y(n1)y(2)u(n2)u(2)e(n1)e(2)y(nN 1)y(N)u(nN)u(N)e(nN 1)e(N)ai1)e(n1)y(n高斯噪声序列的概率密度函数为0的高斯噪声序列，因为已假定

7、e(k)是均值为, 11 exo$(y m)22 2 2(2 )22和m为y的方差和均值，那么1 1 2 ?exp 2 e (k) (22)2对于e(k)符合高斯噪声序列的极大似然函数为式中y为观测值,L(Yn , ) Le(n 1),e( n 2), exp(2 2尸i2(n 1)L(Yn ,(2.14)(2.15)对上式(2.17)(2 2)7等号两边取对数得12)7lnL(YN , ) ln (2或写为ln L(Yn ,)Nln22求In L(Yn ,)对2的偏导数,lnLM ,),e(ne2(n(YnIn exp(N) fe(n2)e2(n1) fe(n 2)1FTfe(n N)T(Y

8、n2 eNeN)Nln21 n22k令其等于 0,可得n Ne2(k)n 1N)-(2exp(eNeN)2 2(2.16)e2(k)1(2.17)Nln21_ T2 eNeN(2.18)(2.19)(2.20)式中n N2en 1(k)1nNe2(k)1JN 2k n 1N(2.21 )2越小越好，因为当方差 小1 n N2 k n2最小时，e2(k)i2e (k)最小，即残差最小。因此希望(2.22 )2的估值取最(2.23 )2 2 min JN因为式(2.10 )可理解为预测模型，而 e(k)可看做预测误差。因此使式(2.22 )最小 就是使误差的平方之和最小， 即使对概率密度不作任何假

9、设， 这样的准则也是有意义的。 因 此可按J最小来求a1, ,a,b0, bn,c1, cn的估计值。由于e(k)式参数a1, ,a,bo, bn,c, cn的线性函数，因此J是这些参数的二次型函数。求使ln L(Yn ,)最大的，等价于在式(2.10 )的约束条件下求使J为最小。由于J对C|是非线性的，因而求J的极小值问题并不好解，只能用迭代方法求解。 求J极小值的常用迭代算法有拉格朗日乘子法和牛顿 -拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤 如下：(1)确定初始的°值。对于°中的a1, ,a,b0,bn可按模型(2.24 )e(k)a(z 1)y(k) b(z

10、 1)u(k)用最小二乘法来求，而对于0中的c Cn可先假定一些值。(2)计算预测误差给出并计算式中e(k) y(k)1 n N2ky(k)(2.25 )e2(k)n 1n N2e1(k)(2.26 )(3)计算J的梯度和海赛矩阵n N- e(k)k n 1e(k)(2.27 )e(k)e(k)e(k)e(k)e(k) e(k)Te(k)aianbobnC1Cne(k)y(k)aiy(k1)any(kn)bou(k) 0u(k 1)bnu(k n)aaiqe(k 1)Cne(kn)y(k i) se(k3 C2-e(k 2)C e(k n)Cn（2.28）aiaiai即e(k)y(kni)e(

11、kCjj)（2.29）aij iai同理可得e(k)u(kni)e(k Cjj)（2.30）bijibe(k)e(kni)e(k Cjj)（2.31）Cij1Ci将式（2.29 ）移项化简，有y(k i)-e(k)ainj 1e(k j) ne(k j)Cj5aij 0ai（2.32）因为e(kj)e(k)z j（2.33）由e(kj）求偏导，故e(kj)e(k)z j（2.34）aiai将（2.34 ）代入（2.32 ）,所以y(k i)n c e(k cjj 0aj) nijc(zc e(k)z jcje(k) njcjz ai j 0Cnz n（2.35）01) 1aqz1所以得c(z1

12、) e(k)y(ki)（2.36）ai同理可得（2.30 ）和（2.31）为c(z 1)e(k)biu(ki)(2.37 )c(z 1)e(k)Cie(ki)(2.38 )根据（2.36 ）构造公式c(z1)非(ij)yk (ij)j y(k i)(2.39 )aj将其代入（2.36），可得c(z1)耿(ij) c(z1)-e(k)(2.40)ajai消除c（z 1）可得e(k)e(ki j)e(ki 1)(2.41 )aiaja1同理可得（2.37 ）和（2.38）式e(k)e(ki j)e(ki)(2.42 )bibjb0e(k)e(ki j)e(ki 1)(2.43 )Cicj式（2.2

13、9 ）、式（2.30）:和式（2.31 ）均为差分方程，这些差分方程的初始条件为0，些偏导数分别为y（k），u（k）和e（k）的线性函数。下面求关于的二阶偏导数，即2J2TnN e(k) e(k) nNe(k) 2e(k)e(k)2-k n 1(2.44)当接近于真值 时,e（k）接近于0。在这种情况下，式（2.44）等号右边第2项接近可通过求解这些差分方程，分别求出e（k）关于ai ,a,b0,bn,G,cn的全部偏导数，而这于0, £可近似表示为(2.45)nN e(k) e(k) T则利用式(2.45 )计算2比较简单。按牛顿-拉卜森计算的新估值1，有(2)1(2.46)如果重

14、复(2)至(4)的计算步骤，经过r次迭代计算之后可得r，近一步迭代计算可得I)12 /(2.47)10 4(2.48)则可停止计算，否则继续迭代计算。式(2.48 )表明，当残差方差的计算误差小于0.01%时就停止计算。这一方法即使在噪声比较大的情况也能得到较好的估计值三实验内容设siso系统差分方程为y(k) a(k1)a2y(k 2)du(k 1) b2u(k 2)(k)其中极大似然估计法默认 e(k)为:e(k) C(z1) (k)(k)G (k 1) C2 (k 2)辨识参数向量为=a1a2Rb2 C1 C2T式中，(k)为噪声方差各异的白噪声或有色噪声；u(k)和y(k)表示系统的输

15、入输出变量。四实验流程图五代码实现程序如下：U=1.147 0.201 -0.787 -1.584 -1.052 0.866 1.152 1.573 0.6260.433 -0.958 0.810 -0.044 0.947 -1.474 -0.719 -0.086 1.0991.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.8900.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.

16、582 1.6090.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.3090.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.81

17、0-0.277 0.983 -0.288 0.846 1.325 0.723 0.713 0.643 0.4630.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.6810.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.2421.680 -1.236 -0.

18、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.0131.518 -0.553 -0.987 0.167 -1.445 0.630 1.255 0.311 -1.7260.975 1.718 1.360 1.667 1.111 1.018 0.078 -1.665 -0.7601.184 -0.614 0.994 -0.089 0.947 1.706 -0.395 1.222 -1.3510.231 1.425 0.114 -0.689 -0.704 1

19、.070 0.262 1.610 1.489-1.602 0.020 -0.601 -0.020 -0.601 -0.235 1.245 1.226 -0.2040.926 -1.297 %输入数据Y=0.086 2.210 0.486 -1.812 -3.705 -2.688 1.577 2.883 3.7051.642 0.805 -2.088 0.946 -0.039 1.984 -2.545 -1.727 -0.2312.440 3.583 2.915 1.443 3.598 0.702 2.638 3.611 -0.1681.732 0.666 2.377 -0.554 -2.088

20、 2.698 0.189 -1.633 -2.0101.716 -1.641 -1.885 1.061 -0.968 2.911 3.088 -1.629 -1.5333.030 0.614 -1.483 -1.029 -1.948 -1.066 -0.113 -2.144 -2.6260.134 -3.043 -1.341 0.338 2.702 3.813 -1.924 -2.813 -1.7953.002 1.027 1.027 2.755 3.584 1.737 -0.837 -0.617 1.7032.045 -2.886 -0.542 -2.991 -1.859 3.045 0.0

21、68 -0.375 0.4511.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.3400.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

22、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.1481.816 0.055 -1.856 0.269 -1.323 -2.486 1.958 0.823 2.4812.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

23、.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.1803.261 -2.755 -0.536 -1.171 -0.905 -3.303 -0.834 2.490 3.0390.134 1.901%输岀数据n a=2; nb=1; nc=2;d=1;nn=max (na,n c);L=size(Y,2);xiek=zeros( nc,1); %白噪声估计初值yfk=zeros (nn ,1); %yf(k-i)ufk=zeros( nn ,1); %uf

24、(k-i) xiefk=zeros (n c,1); %vf(k-i) thetae_ 仁 zeros( na+nb+1+ nc,1); %参数估计初值P=eye( na+n b+1+ nc);for k=3:L组建h（ k）%构造向量phi=-Y(k-1);-Y(k-2);U(k-1);U(k-2);xiek; %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+n b

25、+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+n b+2:na+n b+1+nc,k)'*xiefk(1: nc); %vf(k)%更新数据thetae_1=thetae(:,k);for i=n c:-1:2xiek(i)=xiek(i-1);xiefk(i)=xiefk(i-1);endxiek(1)

26、=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');lege nd('a_1','a_2'); axis(0 L -2 2);figure(2)plot(1:L,thetae( na+1: na+n b+1,:); xlabel('k'); ylabel('参数估计 b');lege nd('b_1','b_2'); axis(