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

北京工商大学系统建模与辨识?课程上机实验报告2021年秋季学期专业名称 ：控制工程上机题目:用极大似然法进行参数估计专业班级：计研3班学生姓名： 王瑶 吴超学 号：10011316259 10011316260指导教师:刘翠玲2021 年1 月实验目的通过实验掌握极大似然法在系统参数辨识中的原理和应用。实验原理1极大似然原理设有离散随机过程vj与未知参数e有关，假定概率分布密度 f(v)。如果我们 得到n个独立的观测值 Vi,V2, - ,Vn,那么可得分布密度 f(Vi|e), f (V2 0) , , f(VnB)。 要求根据这些观测值来估计未知参数e，估计的准那么是观测值 vk的出现概率为最大。为此，定义一个似然函数LM,V2：,Vnp)= f(Vi|e)f(V2ef(Vn|e)(1.1)上式的右边是n个概率密度函数的连乘，似然函数L是的函数。如果L到达极大值，Vk的出现概率为最大。因此，极大似然法的实质就是求出使l到达极大值的a的估值e。为了A便于求n,对式(1.1)等号两边取对数，那么把连乘变成连加，即nln L =乏 ln f (V 6)(1.2 )i =1由于对数函数是单调递增函数，当对8的偏导数，令偏导数为 0,可得:ln LA解上式可得臼的极大似然估计 日讯。L取极大值时，InL也同时取极大值。求式(1.2)=0(1.3 )2系统参数的极大似然估计Newton-Raphson法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每 L 次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数估计值 得算法。本质上说，它只是一种近似的极大似然法。设系统的差分方程为a(z) y(k) = b(z)u(k) + 土 (k)(2.1)式中a(zJ) =1 qz . , anzb(z ) = b0 b1z4 . bnz因为-(k)是相关随机向量，故(2.1 )可写成a(z*)y(k) =b(z*)u(k) +c(z*)8(k)(2.2 )式中c(zU)E(k)=顷k)(2.3 )c(z、)=1 +c1z+CiZJ(2.4 )&(k)是均值为0的高斯分布白噪声序列。多项式a(zW) , b(zW)和c(z)中的系数a1.,a,b0/L bn,c1/L cn和序列任)的均方差。都是未知参数。设待估参数e =aian bobn GG 】(2.5)并设y(k)的预测值为 AAAAAy(k) = -a1 y(k -1) v-an y(kn) b0 u(k)一bn u(kn)AA(2.6)G e(k 1) 丁cn e(k -n)式中e(k -i)为预测误差；ai , bi , ci为a,耳，q的估值。预测误差可表示为 n ne(k) = y(k) -y(k) = y(k) -aiy(k -i) 、bu(k-i)M ce(k T) =(1 a1z- anz)y(k) -(b0 b1z- bnzq)u(k) 1_2_n 、(ci z +c2z +cn z )e(k)(2.7)或者 1n1n(1 ci z Pn z )e( k) = (1 a z an z ) y( k) -_n、 (b + b1 z +bn z )u(k)(2.8)因此预测误差ie(k)满足关系式c(zJ1)e(ka(z_1)y(kb(z_1)u(k)(2.9 )式中 a( z ) = 1 a zan z八， A A d . A n b(z ) =b 灯 zz11nc(z)=1Gzcnz假定预测误差e(k)服从均值为0的高斯分布，并设序列e(k)具有相同的方差。2。因 为&Q与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 T)any(k - n) -b0u(k T) -nu(k T) -bnu(kn)c1e(k-1) 一cn(kn), k = n+1,n + 2,(2.11 )或写成 nnne(k)=y(k)+ aiy(ki)W Ku(ki) Ge(k i)(2.12)i 兰i=0i =1令k=n+1,n+2,n+N,可得e(k)的N个方程式，把这 N个方程式写成向量-矩阵形式一-y(n)- y(1)u(n+1)u(1)e(n)e(1) 1_y(n+1)-y(2)u(n+2) -u(2)e(n + 1)e(2)一y(n + N _1)-y(N)u(n十N广u(N)e(n 十 N _1)e(N)_：JN =B(k)是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为式中-y(n+1)y(n+2)n+1) |e(n+2)anbo(2.13 )_y(n+N)_e(n + N)因为已假定,11,、2、f =1 翎 (y - m)2 22二(2二)2式中y为观测值，。2和m为y的方差和均值，那么,1r 12,、f = exp-we(k)22(2二一-)2对于e(k)符合高斯噪声序列的极大似然函数为(2.14 )(2.15 )1exp(2二；2)21222,e (n 1) e (n 2) e (n N) 2。1exp(2二2)2T 、2 eNeN)2。(2.16 )L(K)=fexpq：W)(2二了2)2*对上式(2.17)等号两边取对数得,、，1.,1 t 、 N _ _ln L(Yn 、，；- ) =ln一 ln exo( - 一一)= 一一ln 2 -一 ln ；-2222(2二)22。2(2.17 )1_T _2 eN eN 2L(Yn 代。)=Le(n +1),e(n+2),，e(n + N)|B = f e(n+1)们 f e(n+2)忡fe(n + N)B(2.18 )或写为NN 91 n N 9(2.19 )ln L(Yn 已。)=一一ln2二 一一ln；-一 e2(k)222二 k求lnL(YN。,。)对。2的偏导数，令其等于 0,可得式中冬工一鼻土二如=0丁 二【nNe2(k)=1nNe2(k)二jN 5 N 2u i NJ =1e2(k)2静iO2越小越好，因为当方差 小22a 取小时， e(k)取小，即残差取小。因此希望U2=3minJN(2.20 )(2.21 )(2.22 )a2的估值取最(2.23 )因为式(2.10 )可理解为预测模型，而 e(k)可看做预测误差。因此使式(2.22)最小 就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准那么也是有意义的。 因 此可按J最小来求%,a,b0,aGL cn的估计值。由于e(k)式参数a1,a,b0，bn,G,cn的线性函数,因此J是这些参数的二次型函数。求使lnL(YN|e,b)最大的等价于在式(2.10 )的约束条件下求 #使J为最小。由于J对G是非线性的，因而求J的极小值问题并不好解，只能用迭代方法求解。求J极小值的常用迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤如下：(1) 确定初始的&值。对于 洗中的a,a, b0,bn可按模型e(k) =a(z ) y(k)b(z )u(k)(2.24)A用最小二乘法来求，而对于 日。中的G ,%可先假定一些值。(2) 计算预测误差Ae(k) =y(k) y(k)(2.25 )给出J =1孔2)2k并计算二2 七?)N ki(2.26 )JT-2J(3) 计算J的梯度 涕 和海赛矩阵 一-,有-;2JcQn N= e(k)k=n 1；:e(k)c9(2.27 )式中T(k) _ fe(k). ?e(k) fe(k). ?e(k)fe(k). ；e(k).:e(k)：a一七一：a1jan.:b0.:bnEgjcny(k) ay(k-1) any(k -n)-bou(k)-biu(k-1)一一bnu(k-n)- .:aice(k -1) -Cne(k - n)一勺一Cn2(2.28 )同理可得.:e(k).:airai(2.29 )毛(k)n=-u(k i) - L Cjj日fe(k - j)A(2.30 )业-e(ki)Cj 史3CijdCCi将式(2.29)移项化简，有y(k_i)=+ CjM Cj caij 1caij =0:a因为e(k - j) = e(k)zj(2.31 )(2.32 )(2.33 )(2.34 )(2.35 )由e(k j)求偏导，故fe(k - j) _ fe(k)z-j将(2.34 )代入(2.32 ),所以y(k-i)G Cj=? Cj=2 不三j =0： aij -0: a: ai j 心c(z) = 1 , qz , cnz所以得c(zLWi= y(k -i)(2.36 )同理可得2.30 和2.31 为根据2.36 构造公式c(z)业:bi-u(k -i)必=_e(ki)-：Ci伞)* ：(i 一 j) = yk (i 一 j) 一 j =y(k i)将其代入2.36 ,可得c(T(2.40(2.37 )(2.38 )(2.39 )消除cz可得:ai:ajfaV j1 ( 2.38 )式-：e(k)=fe(k -i j)_ fe(k -i)(2.42 )A%fe(k) _：e(k -ij) _：e(k -i 1)(2.43 )9e.:e(k) fe(k - i j)fe(k -i 1)0,同理可得2.37 和式2.29 、式2.30 和式2.31 均为差分方程，这些差分方程的初始条件为可通过求解这些差分方程，分别求出ek关于a：.,a,b0,bn,G,G的全部偏导数，而这些偏导数分别为yk , uk和ek的线性函数。下面求关于 e的二阶偏导数，即Tn N，二 e(k) 22：2e(k)(2.44 )A当e接近于真值e时，ek接近于0。在这种情况下，式2.44 等号右边第2项接近C J于0, T可近似表示为u2(2.45 )皂=;业迪T_了 k占1 F _ 七tt . .:2J _那么利用式(2.45 )计算J2比拟简单。H2按牛顿-拉卜森计算0的新估值 奋，有旨i =go |(旦)旦(2.46 )重复(2)至(4)的计算步骤，经过r次迭代计算之后可得 隹,近一步迭代计算可得(2.47 )如果(2.48 )E1 -r4:10； r那么可停止计算，否那么继续迭代计算。式(2.48 )说明，当残差方差的计算误差小于0.01%时就停止计算。这一方法即使在噪声比拟大的情况也能得到较好的估计值0。三实验内容设SIS。系统差分方程为y(k)a1y(k -1) a2y(k - 2) = bu(k -1) b2u(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 ELML 算法的程序流程图五代码实现程序如下：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.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 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.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.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.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.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 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.905 -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);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可知,估计值 a1=-0.1222 , a2=0.0365 , b1=1.7962 , b2= 0.0004 , C1= 0.0655 , C2= -0.0018a计估数参1.510.50-0.5-1-1.5-2 020406080100120140160180200k图1A参数波形变化情况b计估数参21.510.50-0.5-1-1.5020406080100120140160180200图2 B参数波形变化情况C计估数参图3C参数波形变化情况图10-11图10-12有162个数据，结果如下： thetael =-0.70950.04891.8688-0.9850-0.1768-0.0005可知，估计值 a=-0.7095 , a2=0.0489 , b1= 1.8688 , b2= -0.9850 , C1= -0.1768 , C2= -0.0005a计估数参k实用文档图1A参数波形变化情况图2 B参数波形变化情况b计估数参图3 C参数波形变化情况C计估数参六实验总结通过此次实验，对极大似然法的原理和方法进行了研究和对其matlab的仿真，可发现其离线辨识效果很好，但是其在线辨识的计算量很大， 实验数据不同，出现的效果也不一样， 递推的极大似然法呈现病态，所得参数估计取值很不稳定，以至于不能采用。实用文档