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

用极大似然法进行参数估计 《系统辨识》上机实验 报告 软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载 北京工商大学 (2014年秋季学期) 专业名称 : 控制工程 上机题目 : 极大似然法进行参数估计 专业班级 : 2015年 1 月 一 实验目的 通过实验掌握极大似然法在系统参数辨识中的原理和应用。 《系统辨识》上机实验报告 二 实验原理 1 极大似然原理 ,设有离散随机过程与未知参数有关，假定已知概率分布密度。如果我们f(V,){V}kk 得到n个独立的观测值…，则可得分布密度，，…,。f(V,)f(V,)f(V,)V,V,V12,12nn ,要求根据这些观测值来估计未知参数，估计的准则是观测值{}的出现概率为最大。{V}k 为此，定义一个似然函数 L(V,V,？,V,),f(V,)f(V,)？f(V,)12n12n (1.1) ,上式的右边是n个概率密度函数的连乘，似然函数L是的函数。如果L达到极大值，{V}k , ,,的出现概率为最大。因此，极大似然法的实质就是求出使L达到极大值的的估值。为了 , ,便于求，对式(1.1)等号两边取对数，则把连乘变成连加，即 n (1.2) lnL,lnf(V,),i,1i 由于对数函数是单调递增函数，当L取极大值时，lnL也同时取极大值。求式(1.2),对的偏导数，令偏导数为0，可得 ,lnL,0,, (1.3) , ,,ML解上式可得的极大似然估计。 2 系统参数的极大似然估计 Newton-Raphson法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每L次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数估计值得算法。本质上说，它只是一种近似的极大似然法。 设系统的差分方程为 ,1,1 (2.1) a(z)y(k),b(z)u(k)，,(k) 式中 ,,,11nazazaz()1...,，，，1n ,,,11nbzbbzbz()...,，，，01n ,(k)因为是相关随机向量，故(2.1)可写成 ,1,1,1 (2.2) a(z)y(k),b(z)u(k)，c(z),(k) 式中 ,1 (2.3) c(z),(k),,(k) ,1,1,n c(z),1，cz，？，cz (2.4) 1n,1,1,1,(k)是均值为0的高斯分布白噪声序列。多项式a(z)，b(z)和c(z)中的系数 《系统辨识》上机实验报告 和序列的均方差都是未知参数。 a,a,b,？b,c,？c,{,(k)}1,？0n1n 设待估参数 T, c？c (2.5) b？b,,[a？a1n0n1n 并设的预测值为 y(k) ,,,,, y(k),,ay(k,1),？,ay(k,n)，bu(k)，？，bu(k,n)，1n0n ,, (2.6) ce(k,1)，？，ce(k,n)1n ,,, 式中为预测误差;，，为，，的估值。预测误差可表示为 abce(k,i)abciiiiii nn,,,， e(k),y(k),y(k),y(k),,ay(k,i)，bu(k,i)，,,ii,,1,0，ii n,,,,,,，,1,n,1,n ce(k,i),(1，az，？，az)y(k),(b，bz，？，bz)u(k),1n01n,i,，i,1 ,,,,1,2,n (2.7) (cz，cz，？，cz)e(k)12n 或者 ,,,,,1,n,1,n = (1，cz，？，cz)e(k)(1，az，？，az)y(k),1n1n ,,,,1,n (2.8) (b，bz，？，bz)u(k)01n 因此预测误差满足关系式 ,,e(k) ,,,,1,1,1 (2.9) c(z)e(k),a(z)y(k),b(z)u(k)式中 ,,,,1,1,na(z),1，az，？，az 1n ,,,,,1,1,nb(z),b，bz，？，bz 01n ,,,,1,1,nc(z),1，cz，？，cz 1n2,,,假定预测误差服从均值为0的高斯分布，并设序列e(k)具有相同的方差。因e(k) ,,,2,1,1,1,c(z)a(z)b(z),,为e(k)与，和有关，所以是被估参数的函数。为了书写方便，,把式(2.9)写成 ,1,1,1 (2.10) c(z)e(k),a(z)y(k),b(z)u(k) e(k),y(k)，ay(k,1)，？，ay(k,n),bu(k,1),bu(k,1),？, 1n01 bu(k,n),ce(k,1),？,c(k,n),k,n，1,n，2,？ (2.11) n1n 或写成 nnn (2.12) e(k),y(k)，ay(k,i),bu(k,i),ce(k,i),,,iii,1,0,1iii e(k)令k=n+1,n+2,…,n+N,可得的N个方程式，把这N个方程式写成向量-矩阵形式 e,Y,,, (2.13) NNN 式中 《系统辨识》上机实验报告 a，，1,,？y(n，1)e(n，1)，，，，,,,,,,,,ay(n，2)e(n，2)n,,,,,, , , Y,e,,,NNb,,,,？？0,,,,,,,,？y(n，N)e(n，N)，，，，,,b,,n，， ,y(n)？？？e(n)e(1),y(1)u(n，1)u(1)，， ,,,y(n，1)？？？e(n，1),y(2)u(n，2)e(2)u(2),, ,,N,,？？？？？？ ,,？？？,y(N)u(n，N)u(N)e(n，N,1)e(N),y(n，N,1)，， ,,e(k)因为已假定是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为 112f,exp[,(y,m)]12,222(2),, (2.14) 2,式中y为观测值，和m为y的方差和均值，那么 112f,exp[,e(k)] (2.15) 122,22(2),, 对于符合高斯噪声序列的极大似然函数为 e(k) ,,,,,,L(Y,),L[e(n，1),e(n，2),？,e(n，N)],f[e(n，1)]f[e(n，2)]？f[e(n，N)]N 1111222T,exp{,[e(n，1)，e(n，2)，？，e(n，N)]},exp(,ee)NNNN22,,222222(2)(2),,,, (2.16) 或 T,,1()()Y,,Y,,NN,,(,)exp[] LY,, (2.17) NN2,222(2,,) 对上式(2.17)等号两边取对数得 111NN2TT,,,,ln(,)lnlnexp()ln2lnLY,，,ee,,,,ee NNNNNN22,2222,22(2),, (2.18) 或写为 n，NNN122,,,, (2.19) lnL(Y,),,ln2,ln,e(k),N2,2221k,n，2,lnL(Y,,,)求对的偏导数，令其等于0，可得 N n，N,,,lnL(Y,)N12N ,,，e(k),0 (2.20) ,224,,,,221k,n， 则 《系统辨识》上机实验报告 ,n，Nn，N1212222 (2.21) ,,ek,ek,J()(),,NNN211k,n，k,n， 式中 n，N12 (2.22) J,e(k),21k,n，2222越小越好，因为当方差最小时，最小，即残差最小。因此希望的估值取最,,,e(k) 小 2,2 (2.23) ,,minJN 因为式(2.10)可理解为预测模型，而e(k)可看做预测误差。因此使式(2.22)最小就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。因此可按J最小来求的估计值。 a,a,b,？b,c,？c1,？0n1n 由于e(k)式参数的线性函数，因此J是这些参数的二次型函数。a,a,b,？b,c,？c1,？0n1n ,,求使lnL(Y,,,)最大的，等价于在式(2.10)的约束条件下求使J为最小。由于J对,,N 是非线性的，因而求J的极小值问题并不好解，只能用迭代方法求解。求J极小值的常用ci 迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤如下: ,, (1)确定初始的值。对于中的a,a,b,？b可按模型 ,,001,？0n ,,,1,1 (2.24) e(k),a(z)y(k),b(z)u(k) ,c,？c1n,0用最小二乘法来求，而对于中的可先假定一些值。 (2)计算预测误差 , e(k),y(k),y(k) (2.25) 给出 n，N12J,e(k),21k,n， 并计算 ,n，N122,,e(k),N1k,n， (2.26) ,J2,J,,(3)计算J的梯度 和海赛矩阵 ,有 2,, n，N,,Je(k), (2.27) e(k),,,,,1k,n， 式中 《系统辨识》上机实验报告 T ，，e(k)e(k)e(k),,,,e(k),e(k),e(k),e(k) ？,？？,,,c,c,b,b,aa,,,10n，1nn， ,e(k), ,[y(k)，ay(k,1)，？，ay(k,n),bu(k),bu(k,1),？,bu(k,n),1n01n,a,aii ce(k,1),？,ce(k,n)]1n ,e(k,1),e(k,2),e(k,n) ,y(k,i),c,c,,c (2.28) ？n12,a,a,aiii即 n,ek,ek,j()() ,yk,i,c (2.29) (),j,a,a,1jii 同理可得 n,ek,ek,j()(),,uk,i,c (2.30) (),j,b,b,1jii n,ek,ek,j()(),,ek,i,c (2.31) (),j,c,cj,1ii将式(2.29)移项化简，有 nn,ek,ek,j,ek,j()()()yk,i,，c,c (2.32) (),,jj,a,a,a,1,0jjiii因为 ,j (2.33) e(k,j),e(k)z 由求偏导，故 e(k,j) j,,ek,j,ekz()(), (2.34) ,a,aii 将(2.34)代入(2.32)，所以 ,jnnn,,,,ekjekzek()()(),j,,,,ykicccz (2.35) (),,,jjj,,,aaa000j,j,j,iii ,1,1,nc(z),1，cz，？，cz 1n所以得 ,e(k),1c(z),y(k,i) (2.36) ,ai 《系统辨识》上机实验报告 同理可得(2.30)和(2.31)为 ,e(k),1 (2.37) c(z),,u(k,i),bi ,e(k),1 (2.38) c(z),,e(k,i),ci 根据(2.36)构造公式 ,e[k,(i,j)],1c(z),y[k,(i,j),j],y(k,i) (2.39) ,aj 将其代入(2.36)，可得 ,ek,i,j,ek[()](),1,1cz,cz()() (2.40) ,a,aji ,1消除可得 c(z) ,,,，,,，e(k)e(kij)e(ki1),, (2.41) ,,,aaaij1 同理可得(2.37)和(2.38)式 ,,,，,,e(k)e(kij)e(ki),, (2.42) ,,,bbbij0 ,,,，,,，e(k)e(kij)e(ki1),, (2.43) ,,,cccij1 式(2.29)、式(2.30)和式(2.31)均为差分方程，这些差分方程的初始条件为0， a,a,b,？b,c,？c可通过求解这些差分方程，分别求出e(k)关于的全部偏导数，而这1,？0n1n些偏导数分别为，和的线性函数。下面求关于的二阶偏导数，即 {y(k)}{u(k)}{e(k)}, T22n，Nn，N,,,,Je(k)e(k)e(k)，，,， (2.44) e(k),,22,,,,,,,,,,，，11k,n，k,n， , ,,当接近于真值时，e(k)接近于0。在这种情况下，式(2.44)等号右边第2项接近 2,J于0，可近似表示为 2,, T2，nN,J,e(k),e(k)，， (2.45) ,,2,,,,,,,,，，1,，kn 《系统辨识》上机实验报告 2,J则利用式(2.45)计算比较简单。 2,, , (4)按牛顿-拉卜森计算的新估值，有 ,,1 2,,，，,,JJ,1,, (2.46) ,,(),10,,2,,,,,，，,0 , ,r重复(2)至(4)的计算步骤，经过r次迭代计算之后可得，近一步迭代计算可得 2,,，，,,JJ,1,, (2.47) ,,(),rr，1,,2,,,,,，，,r如果 ,,22,,,r，1r,4,10 (2.48) ,2,r 则可停止计算，否则继续迭代计算。 0.01%时就停止计算。这一方法即使在噪声式(2.48)表明，当残差方差的计算误差小于 , ,比较大的情况也能得到较好的估计值。 三 实验内容 设SISO系统差分方程为 y(k)，ay(k,1)，ay(k,2),bu(k,1)，bu(k,2)，,(k)1212其中极大似然估计法默认为: ek() ,1 ekCzkkckck()()()()(1)(2),,，,，,,,,,12 T,b[aab辨识参数向量为 , c c] 122121 式中，为噪声方差各异的白噪声或有色噪声;u(k)和y(k)表示系统的输入输出变量。 ,(k) 《系统辨识》上机实验报告 四 实验流程图 五 代码实现 程序如下: U=[1.147 0.201 -0.787 -1.584 -1.052 0.866 1.152 1.573 0.626... 0.433 -0.958 0.810 -0.044 0.947 -1.474 -0.719 -0.086 1.099... 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:2 xiek(i)=xiek(i-1); xiefk(i)=xiefk(i-1); end xiek(1)=xie; xiefk(1)=xief; for i=nn:-1:2 yfk(i)=yfk(i-1); ufk(i)=ufk(i-1); end yfk(1)=yf; ufk(1)=uf; end thetae_1 figure(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.1222 0.0365 1.7962 0.0004 0.0655 -0.0018 可知，估计值a=-0.1222，a=0.0365，b=1.7962，b= 0.0004，c= 0.0655，c= -0.0018 121212 2 a11.5a2 1 0.5 0 参数估计a-0.5 -1 -1.5 -2 020406080100120140160180200k 图1 A参数波形变化情况 2 b1 b1.52 1 0.5 0参数估计b -0.5 -1 -1.5 020406080100120140160180200k 《系统辨识》上机实验报告 图2 B参数波形变化情况 2 c1 1.5c2 1 0.5 0 参数估计c-0.5 -1 -1.5 -2 020406080100120140160180200k 图3 C参数波形变化情况 图10-11 图10-12有162个数据，结果如下: thetae_1 = -0.7095 0.0489 1.8688 -0.9850 -0.1768 -0.0005 可知，估计值a=-0.7095，a=0.0489，b= 1.8688，b= -0.9850，c= -0.1768，c= -0.0005 121212 2 a11.5a2 1 0.5 0 参数估计a-0.5 -1 -1.5 -2 020406080100120140160k 《系统辨识》上机实验报告 图1 A参数波形变化情况 2 b1 b1.521 0.5 0参数估计b -0.5 -1 -1.5 020406080100120140160k 图2 B参数波形变化情况 2 c11.5c21 0.5 0 参数估计c-0.5 -1 -1.5 -2 020406080100120140160k 图3 C参数波形变化情况 六 实验 总结 初级经济法重点总结下载党员个人总结TXt高中句型全总结.doc高中句型全总结.doc理论力学知识点总结pdf 通过此次实验，对极大似然法的原理和方法进行了研究和对其matlab的仿真，可发现 现的效果也不一样，其离线辨识效果很好，但是其在线辨识的计算量很大，实验数据不同，出 递推的极大似然法呈现病态，所得参数估计取值很不稳定，以至于不能采用。