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

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

2、。如果L达到极大值,VkA的出现概率为最大。因此，极大似然法的实质就是求出使L达到极大值的二的估值二。为A.i =1了便于求,对式（1.1）等号两边取对数，则把连乘变成连加，即(1.2)由于对数函数是单调递增函数,当L取极大值时,lnL也同时取极大值。求式（1.2）对，的偏导数，令偏导数为0,可得;In L(1.3)解上式可得二的极大似然估计ml。2系统参数的极大似然估计Newton-Raphson 法实际上就是一种递推算法 ，可以用于在线辨识。不过它是一种依 每L次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数 估计值得算法。本质上说，它只是一种近似的极大似然法。设

3、系统的差分方程为a(z)y(k)二b(z)u(k)(k)(2.1)式中11na(z) =1 aiz一 - . anz_ b(z 亠)=bo dz . bnZ因为| k)是相关随机向量，故(2.1)可写成a(z)y(k)二b(zJ)u(k) c(z) ；(k)(2.2)式中c(z) ；(k)二(k)(2.3)c(zJ) =1 qzCnZ(2.4)；(k)是均值为 0的高斯分布白噪声序列。多项式a(zJ)，z-1)和c(zJ)中的系数a1;.,a,b0L bn,q，cn和序列R(k)的均方差cr都是未知参数。设待估参数TV -an bobn Gcn 丨(2.5)并设y(k)的预测值为AAAAAy(

4、k) - -aj y(k -1) - an y(k -n) b0u(k)川川bn u(k - n)G e(k1)-cn e(kn)(2.6)式中e(k -i)为预测误差；J，bi，ci为aj，b，q的估值。预测误差可表示为A_ n 找n Ae(k)=y(k)-y(k) = y(k) - -E ajy(k-i)+W bu(k-i) +_ i 1i=0n 乱TA 1A nA A 1A n、 Ge(ki) = (1 a1 z -anz )y(k)(b。bz bn z )u(k)-y八 1 八2_n(c zC2 zcn z )e(k)( 2.7)或者AAAA(1 C1 zcnz)e(k) = (1 a

5、1zanz)y(k)-(bo + b1 z亠+bn z)u(k)(2.8)因此预测误差 (k)?满足关系式c(ze(k) =a(z)y(k) -b(zu(k)(2.9)式中a(zJ) =1 印 亠亠an z*AAAAb(z )=b0bzbn z假定预测误差e(k)服从均值为0的高斯分布，并设序列e(k)?具有相同的方差为g(k)与c(z), a(z)和b(z)有关，所以；2是被估参数v的函数。为了书写方便， 把式(2.9)写成c(z)e(k) = a(z)y(k) -b(z)u(k)(2.10)e(k) = y(k) a(k 一1)亠亠any(k - n) b0u(k -1) tiu(k -1

6、) -bnu(k -n) -Ge(k 一1) - -Cn(k - n), k = n 1,n 2,(2.11)或写成nnne(k)=y(k)+送 aiy(k-i)送 bu(ki)_送 qe(k i)( 2.12)ii zOi二(2.13)令k=n+1,n+2,n+N,可得e(k)的N个方程式，把这N个方程式写成向量-矩阵形式eN- YN_ ：,J N71式中Yn =y(n 1)|y(n+2)e(n 1)e(n +2)anb。丿(n+N)_f(n + N)_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)【eu(n N) u(N)

7、e(n N -1) e(N)L_y( n + N 1)-y(N)沁(k)，是均值为o的高斯噪声序列，高斯噪声序列的概率密度函数为1 1 2 f1 ex)2 (y - m)2 22b(2 二二)22 ，二和m为y的方差和均值，那么1 12f?exp 2e (k)2 22b(2 二二)2对于e(k)符合高斯噪声序列的极大似然函数为因为已假定(2.14)式中y为观测值(2.15 )LM 日,巧=Le(n +1),e(n+2),e(n + N)|8 = fe(n+1)日fe(n+2)|8fe(n + N)日122211 tNexp e (n 1) e(n 2) e (n N)Nexp(eNez)2 T

8、22 y(2二)2(2二)2(2.16 )L(Yn 咕)= expS 一 V (药一八) (2砂 2)2对上式(2.17)等号两边取对数得1ln L(Yn 十)=ln (2.17 )1 TNN-nr ln exo(2eNeN)ln 2ln；二2 厅222(2二)21eT e2 eN eN2-(2.18 )或写为ln L(YnL)= 一号In 2 二Nln；2 字 Ig22；k 出 1(2.19 )求In L(YnR；)对二2的偏导数，：ln L(Yn 叮)令其等于0，可得1 nN 2 e2(k) =02；- k=B 1(2.20)式中2ek =n 1(k)n N2k：e2(krJ(2.21 )

9、2二越小越好，因为当方差小n N e2(k)2 1、二最小时，e2(k)最小，即残差最小。因此希望(2.22 )2-的估值取最(2.23)2 2min JN(2.22)最小就，这样的准则也是有意义的。因因为式(2.10)可理解为预测模型，而e(k)可看做预测误差。因此使式 是使误差的平方之和最小，即使对概率密度不作任何假设 此可按J最小来求&,a,b,bn,C1 / Cn的估计值。由于e(k)式参数a：.,a,g，bnqL q的线性函数，因此J是这些参数的二次型函数。求使In L(Yn6)最大的日，等价于在式(2.10)的约束条件下求 日使J为最小。由于J对Ci是非线性的，因而求J的极小值问题

10、并不好解，只能用迭代方法求解。求J极小值的常用迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤如下(1)确定初始的二0值。对于二0中的a1；a,b,bn可按模型e(k) =a(zy(k) -b(z)u(k)(2.24)用最小二乘法来求，而对于二0中的C1Cn可先假定一些值(2)计算预测误差(2.25)e(k) =y(k) -y(k)给出J f e2(k)2 k出i并计算A 1 仆(2.26)宀 NO2(k)J计算J的梯度和海赛矩阵/J胡nN：e(k)八e(k)警k=n 1门(2.27)式中T 空(k) _ 比(k) 血(k)晁(k) 空(k) 5e(k)6e(

11、k)1 c9 |_ ca-itan和0cbncC|ccn-bnu(k - n) _y(k) a-y(k -1)any(k -n)-bu(k)-bu(k -1)-:a:a.c1e(k T)-cne(k n)二 y(k _i) _g.:e(k -1)：e(k-2)C2 -：e(k -n)(2.28)：e(k - j)：:e(k)-a,同理可得ie(k).:bi55 ).:bi(2.30)-:e(k)-：Ci-e(k(2.31)将式（2.29）移项化简，有y(k-i)acaizcfj)因为由e（k - j）求偏导，故.:ai八 CjJJj z0口(2.32)e(k - j) = e(k)zJ(2.3

12、3)e(k-j).:e(k)z.ai(2.34)将（2.34 ）代入（2.32）,所以y(k -i)八 cj出-：e(k-j)八j=0；:e(k)z_lCj.ai所以得同理可得（2.30）和.:e(k) n.aiCjZ_j=0(2.35)c(z)c(z：e(k)-:ai(2.31)为c(ze(k)-bic(ze(k)-：Ci二1C|Z亠 亠cnzy(k -i)-u(k -i)-e(k - i)-n(2.36)(2.37)(2.38)C(z4) 二 yk -(i - j) - j = y(k -i)(2.39)将其代入(2.36),可得陀：卅一(1)胡z)兰巴6_ai消除c(z)可得；：e(k)

13、：e(k - i j)；：e(k -i 1).:ai.:aj同理可得(2.37)和(2.38)式fe(k)fe(k-i j) ：e(k-i)bi_bj；:b0：e(k)：e(k - i j)：e(k -i 1)Ci. C j: Ci(2.40 )(2.41 )(2.42 )(2.43)式(2.29 )、式(2.30 )和式(2.31 )均为差分方程，这些差分方程的初始条件为0,可通过求解这些差分方程，分别求出e(k)关于3|；.,a,b0,bn,c1/ cn的全部偏导数，而这些偏导数分别为y(k)，u(k)和e(k)的线性函数。下面求关于二的二阶偏导数，即T：2J _ : N ：e(k) ：e

14、(k) - 2k=n 1n：N二 e(k)k =n 1；：2e(k)(2.44 )A当二接近于真值时，e(k)接近于0。在这种情况下，式(2.44)等号右边第2项接近I手可近似表示为吹 Wk)”e(k)T(2.45 )则利用式(2.45)计算比较简单。按牛顿-拉卜森计算二的新估值X，有重复(2)至(4)的计算步骤，经过r次迭代计算之后可得,近一步迭代计算可得(2.47)：2 J 1如果(2.48)则可停止计算，否则继续迭代计算。式(2.48)表明，当残差方差的计算误差小于0.01%时就停止计算。这一方法即使在噪声A比较大的情况也能得到较好的估计值r。三实验内容设SISO系统差分方程为y(k)

15、a(k -1) a2y(k - 2) = 0u(k -1) b2u(k - 2)(k)其中极大似然估计法默认e(k)为：e(k) =C(z)(k)二(k) c, (k-1) c2 (k -2)辨识参数向量为 71= a,a2 bb2 ci c2T式中，(k)为噪声方差各异的白噪声或有色噪声；u(k)和y(k)表示系统的输入输出变量四实验流程图岡卑.5 BLM1-漳法的程序诜理图五代码实现程序如下：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.

16、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.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

17、.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.810-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.8

18、89 -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.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.16

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

20、.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 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.0

21、66 -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.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

22、.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.1

23、48 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

24、.755 -0.536 -1.171 -0.905 -3.303 -0.834 2.490 3.039 0.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(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

25、%构造向量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+n b+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+

26、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_仁thetae(:,k);for i=n c:-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); lege nd(a_1,a_2); axis(0 L -2 2); figure(2)plot(1:L,thetae (n a+1: na+