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

1、实用文档北京工商大学系统建模与辨识课程上机实验报告（2016年秋季学期）专业名称：控制工程上机题目：用极大似然法进行参数估计专业班级：计研3班学生姓名： 王瑶 吴超学 号：10011316259 10011316260指导教师：刘翠玲2017 年1 月实验目的通过实验掌握极大似然法在系统参数辨识中的原理和应用。实验原理1极大似然原理设有离散随机过程VJ与未知参数日有关，假定已知概率分布密度 fM|e)。如果我们 得到n个独立的观测值 Vi,V2,M,则可得分布密度 f(Vi|e), fMS,，f(VnB)。 要求根据这些观测值来估计未知参数e，估计的准则是观测值 vk的出现概率为最大。为此，定

2、义一个似然函数LM,V2,Vn|8) = f(Vi|8)f(V2|8fM|8)(1.1)上式的右边是n个概率密度函数的连乘，似然函数L是H的函数。如果L达到极大值，VkA的出现概率为最大。因此，极大似然法的实质就是求出使L达到极大值的日的估值8。为了A便于求9,对式(1.1)等号两边取对数，则把连乘变成连加，即 nln L = £ ln f (Vi 日)(1.2 )由于对数函数是单调递增函数，当 对8的偏导数，令偏导数为 0,可得ln LA解上式可得日的极大似然估计8 ML。i工L取极大值时，InL也同时取极大值。求式(1.2)=0(1.3 )2系统参数的极大似然估计Newton-R

3、aphson法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每 L 次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数估计值 得算法。本质上说，它只是一种近似的极大似然法。设系统的差分方程为a(z")y(k) = b(z,)u(k)十：(k)(2.1)式中a(zJ) = 1 4z)'anzb(z 工)=b0 ' b1z 1 bnz因为产(k)是相关随机向量，故(2.1 )可写成a(z")y(k) =b(z")u(k) +c(z")8(k)(2.2 )式中c(z')Mk)=巴(k)(2.3 )C(Z

4、A) =1 +c1z，+ +cnz,(2.4 )s(k)是均值为 0的高斯分布白噪声序列。多项式a(z,),b(z,)和c(z,)中的系数a1,a,b0,bn,c1,cn和序列e(k)的均方差仃都是未知参数。设待估参数日=aan bobn Gg】(23)并设y(k)的预测值为AAAAAy(k) = -a1 y(k -1)-an y(k -n) b0 u(k) -一bn u(k - n)AA(2.6)c1 e(k -1) cn e(k -n)式中e(k i)为预测误差；1，bi, C；为ai,可，ci的估值。预测误差可表示为nne(k) = y(k) -y(k) = y(k) - -、aiy(k

5、 -i) x bu(k-i) _ i 1i=0二：Ge(k T) =(1 aiz, -：-anz")y(k)7bo biz- -bnz“)u(k) T1_2_n 、(ci z +c2z +cn z )e(k)(2.7)或者 1n1n(1 ci z -»cn z )e(k) = (1 a1 z -» an z ) y(k) -1_n、(bo + bz + +bnz )u(k)(2.8)因此预测误差ie(k)满足关系式c(zJL)e(k) =a(z_1)y(k) -b(z_L)u(k)(2.9 )式中a( z ) = 1 a1 zan z八1 八八一 A.b(z

6、9;) =bo 1blz，bnz"11nc(z)=1gzcnz假定预测误差e(k)服从均值为0的高斯分布，并设序列 Q(k)具有相同的方差 仃2。因 为Q(k)与c(z"), a(z)和b(z,)有关，所以仃2是被估参数日的函数。为了书写方便， 把式(2.9)写成c(z")e(k) =a(z,)y(k) b(z')u(k)(2.10 )e(k) = y(k) a1y(k -1)any(k - n) b0u(k -1) bu* -1) 一 一bnu(kn)c1e(k-1) 一cn(kn), k = n+1,n + 2,(2.11 )或写成 nnne(k)=y

7、(k)十£ aiy(ki)£ bu(ki)£ cie(k -i)(2.12)Ti=0i =1令k=n+1,n+2,n+N,可得e(k)的N个方程式，把这 N个方程式写成向量-矩阵形式式中-y(n+1) 1y(n+2)y(n +N)<T>'N-y(n)-y(n+1)eN- YN- N71Hn+1) Ie(n+2)&n + N)-y(1)-y(2)y(n + N 1)-y(N)(2.13 )anbou(n 1)u(1)u(n 2)u(2)u(n N) u(N)e(n)e(n 1)e(1) 1 e(2)e(n N -1) e(N)因为已假定

8、L(k)是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为f =(27：C 212rejp-口(y -m) )2(2.14 )式中y为观测值，er2和m为y的方差和均值，那么,1r 12八、(2.15 )f =Texp-0e(k)22 -1(2二二)2对于e(k)符合高斯噪声序列的极大似然函数为L(Yn d。)=Le(n +1),e(n+2),，e(n + N)e = f e(n+1)8 f e(n+2帆fe(n + N)刃1 Nexp (2 二二2)2122_2一2e2(n 1) e2(n 2)e2(n N)2。1 Kexp( (2二二 2)212二2T 、 eNeN )或L(2,

9、87;exp(2 o2-1(2二二 2)2对上式(2.17)等号两边取对数得，、，1.,1 T 、 N , c N ,2ln L(Yn ' ,- ) =lnN ln exo( 2 eNeN) = 一 一 ln 2 一 一 ln；-2222(2二二)2(2.16 )(2.17 )12二2T eN eN(2.18 )或写为NN 91 n N 9(2.19 )ln L(YN 二二)二一一ln2二 一一ln；-2 ' e2(k)222二百求帖1(丫；日,仃)对仃2的偏导数，令其等于 0,可得式中W零，:马2(八03二2二2 2二工 i21 nN 2；二二一 e (k)N k n 1n

10、N市2kl1e2(k)1 nN 2J = ' e2(k)2kzn i仃2越小越好，因为当方差 小22灯 取小时， e(k)取小，即残差取小。因此希望minJN(2.20 )(2.21 )(2.22 )a2的估值取最(2.23 )因为式(2.10 )可理解为预测模型，而 e(k)可看做预测误差。因此使式(2.22)最小 就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。 因 此可按J最小来求现,a,b0,bn,。,cn的估计值。由于e(k)式参数a1,a,b0,bn,G,cn的线性函数，因此J是这些参数的二次型函数。求使ln L(Yn，。)最大的S，等价于在式

11、(2.10)的约束条件下求 S使J为最小。由于J对G是非线性的，因而求J的极小值问题并不好解， 只能用迭代方法求解。 求J极小值的常用 迭代算法有拉格朗日乘子法和牛顿 -拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤 如下：(1)确定初始的$0值。对于00中的a1,a, b0,bn可按模型e(k) =a(z )y(k) -b(z )u(k)(2.24)A用最小二乘法来求，而对于日0中的c1, cn可先假定一些值。(2)计算预测误差Ae(k) =y(k)-y(k)(2.25)给出1n N 2J e (k) 2k土 1并计算n N二2v e2(k)N书(2.26)J一,'2J(3)计算

12、J的梯度。9和海赛矩阵 -J ,有:2n N=e(k)k国1fe(k)c9(2.27 )式中T-e(k)二；:e(k).汨(k)fe(k). Fe(k) ?e(k), 鼻.:e(k).:aia.:a1jan.:b0.:bnjc1：:cny(k) a1y(k-1) -any(k-n)-boU(k)biU(k-1)1“- -bnu(k-n)-血ce(k -1) -Cne(k - n)H泪泪(2.28 )；:e(k)jainy(k -J).二 Cjj 1论(k-j)Fai(2.29 )同理可得；：e(k)：:bin=-u(k -i) -，Cjj 1fe(k - j)(2.30 )治(k)n 沦(k-

13、j)二-gk-i cj：Cijd二 Ci将式(2.29)移项化简，有y(k.i)=HCjieUcn Cj 网0 a j 1caij =0tai因为e(k - j) = e(k)z-j(2.31 )(2.32 )(2.33 )(2.34 )(2.35 )由e(k j)求偏导，故：e(k - j) _ ：:e(k)z-j汨 汨将(2.34 )代入(2.32 ),所以y(k-i) Cjkcn Cjk :四zTj =o- ai j -o、a- ai j =oC(z J) = 1 ' C|Z ,'cnz所以得c(z)£ekl = y(k i)(2.36 )同理可得(2.30 )

14、和(2.31 )为c(z,)回i = u(ki)(2.37 )由c(z,)e(k) = -e(k-i)(2.38 )根据(2.36 )构造公式c(z1)邳 _(i _ j) = yk (i j) j = y(k -i)(2.39 )同将其代入(2.36 ),可得c(z严k-(i-j)=c(z)(2.40)二aj二 ai消除c(z)可得(2.41 )(2.42 )(2.43 )沟(k)；:e(k -i j)沦(k -i 1)：:aiFj：:ai同理可得(2.37 )和(2.38 )式fe(k) _ ；:e(k -i j) _ ；:e(k - i)/由：:bo：:e(k) _ fe(k i j)

15、_ ：:e(k -i 1).:Ci论沁式(2.29)、式(2.30)和式(2.31 )均为差分方程，这些差分方程的初始条件为 0,可通过求解这些差分方程，分别求出e(k)关于ai,.,a,b0,bn,G，cn的全部偏导数，而这些偏导数分别为y(k), u(k)和e(k)的线性函数。下面求关于T三 _ nC £e(k)kR1nN，二 e(k)kR 1：2e(k)26的二阶偏导数，即(2.44 )A当8接近于真值8时，e(k)接近于0。在这种情况下，式(2.44)等号右边第2项接近于0,f2J-J可近似表不为3_曾阳k)砥k) T 【J k n 1 F _ F(2.45 )tI 茜:2J

16、_ _则利用式(2.45 )计算J2-±比较简单。32按牛顿-拉卜森计算0的新估值01,有01 =窑1()- 1(2.46 )I的怎心重复(2)至(4)的计算步骤，经过r次迭代计算之后可得 $，近一步迭代计算可得丁. I /2J 1F r 1 - 1 r -()一(2.47 )如果71 - ；- r4:：10(2.48 )则可停止计算，否则继续迭代计算。式(2.48)表明，当残差方差的计算误差小于0.01%时就停止计算。这一方法即使在噪声比较大的情况也能得到较好的估计值0。三实验内容设SISO系统差分方程为y(k) a1y(k-1) a2 y(k。2) = b1u(k。1) b2u(

17、k-2)(k)其中极大似然估计法默认e(k)为：e(k) =C(z,)(k) = (k) g (k1) C2 (k -2)辨识参数向量为8=a1a2灯b2 c 1 c 2 T式中，？ k)为噪声方差各异的白噪声或有色噪声;u(k)和y(k)表示系统的输入输出变量。四实验流程图阳s.s RM1.算法的程序流程图五代码实现程序如下：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

18、 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

19、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 -

20、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.2421.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.0131.518 -0.553 -0.987 0.167 -1.445 0.630 1.255

21、0.311 -1.7260.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.3510.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.2040.926 -1.297 %输入数据Y=0.086 2.210 0.486 -1.812 -3.705 -2.688 1

22、.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.168.1.732 0.666 2.377 -0.554 -2.088 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

23、.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.703.2.045 -2.886 -0.542 -2.991 -1.859 3.045 0.068 -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.340.0.916 1.3

24、96 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.5951.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

25、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.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.90

26、5 -3.303 -0.834 2.490 3.0390.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);xie

27、k; %组建 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-thet

28、ae(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(