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

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

2、|8)=f(Vi|8)f(V2|8fM|8)(1.1)上式的右边是n个概率密度函数的连乘，似然函数L是H的函数。如果L达到极大值，VkA的出现概率为最大。因此，极大似然法的实质就是求出使L达到极大值的日的估值8。为了A便于求9,对式(1.1)等号两边取对数，则把连乘变成连加，即nlnL=lnf(Vi日)(1.2)i工L取极大值时，InL也同时取极大值。求式(1.2)=0(1.3)2系统参数的极大似然估计Newton-Raphson法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每L次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数估计值得算法。本质上说，它只是

3、一种近似的极大似然法。设系统的差分方程为a(z)y(k)=b(z,)u(k)十：(k)(2.1)式中a(zJ)=14z)anzb(z工)=b0b1z1bnz因为产(k)是相关随机向量，故(2.1)可写成a(z)y(k)=b(z)u(k)+c(z)8(k)(2.2)式中由于对数函数是单调递增函数，当对8的偏导数，令偏导数为0,可得lnLA解上式可得日的极大似然估计8ML。c(z,)Mk)=巴(k)(2.3)C(ZA)=1+c1z,+cnzj(2.4)s(k)是均值为0的高斯分布白噪声序列。多项式a(z,),b(z,)和c(z,)中的系数a1.,a,b0,Lbn,c1；-Lcn和序列e(k)的均方

4、差仃都是未知参数。设待估参数日=aanbobnGg】(2.5)并设y(k)的预测值为AAAAAy(k)=-a1y(k-1)一any(k-n)b0u(k)-一bnu(k-n)AAGe(k-1);.cne(k-n)式中e(ki)为预测误差；1，bi,I为ai,可，ci的估值。预测误差可表示为nne(k)=y(k)-y(k)=y(k)-、aiy(k-i)xtu(k-i)_i1i=0ce(k-i)=(1a1zdanz)y(k)-(b0b1zJ-bnz)u(k)-T1_2_n、(CIz+c2z+cnz)e(k)(2.7)或者1n1n(1CIz-cnz)e(k)=(1a1z-anz)y(k)-1(b0+b

5、1z+bnz)u(k)(2.8)因此预测误差ie(k)满足关系式c(zJ1)e(k)=a(z_1)y(k)-b(z_1)u(k)(2.9)式中a(z)=1a1zanz八1八八一A.b(z)=bnz，bnz11nC(Z)=1GZgz假定预测误差e(k)服从均值为0的高斯分布，并设序列e(k)具有相同的方差仃2。因为Q(k)与c(z),2(2)和吊2,)有关，所以仃2是被估参数6的函数。为了书写方便，把式(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(k-1)一一bnu(kn)c1e(k-1)一c

6、n(kn),k=n+1,n+2,(2.11)或写成nnne(k)=y(k)十aiy(ki)bu(ki)cie(k-i)(2.12)Ti=0i=1令k=n+1,n+2,n+N,可得e(k)的N个方程式，把这N个方程式写成向量-矩阵形式(2.6)式中y为观测值，er2和m为y的方差和均值，那么,1r12八、f=rexp-72e(k)22-1(2二二)2对于e(k)符合高斯噪声序列的极大似然函数为L(YNd。)=Le(n+1),e(n+2),，e(n+N帆=fe(n+1)8fe(n+2帆fe(n+N)刃(2.18)或写为NN91nN9lnL(YN二二)二一一ln2二一一ln；-2e2(k)式中-y(

7、n+1)1y(n+2)y(n+N)N-y(n)-y(n+1)eN-YN-N71Hn+1)|e(n+2)e(n+N)-y(1)-y(2)y(n+N1)-y(N)(2.13)anbou(n1)u(1)u(n2)u(2)u(nN)u(N)e(n)e(n1)e(1)【e(2)e(nN-1)e(N)因为已假定L(k)是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为f=(27：C212exo(y-m)2(2.14)(2.15)1Nexp(2二二2)2122_2 一2e2(n1)e2(n2)e2(nN)2二1N；exp(2二二2)212二2T、或L(Y,exp(N)2o2-1(2二二2)2对上式(2.

8、17)等号两边取对数得，1.,1T、N_N12lnL(YN)=lnNlnexo(2eNeN)=一一ln2一一ln；-2222(2二二)2(2.16)(2.17)12二2TeNeN(2.19)222二百求帖1(丫；日,仃)对仃2的偏导数，令其等于0,可得因为式(2.10)可理解为预测模型，而e(k)可看做预测误差。因此使式(2.22)最小就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。因此可按J最小来求现,a,b0,bn,。,cn的估计值。由于e(k)式参数a1,a,b0，bn,G,cn的线性函数，因此J是这些参数的二次型函数。求使lnL(YN，。)最大的S，等价于

9、在式(2.10)的约束条件下求$使J为最小。由于J对G是非线性的，因而求J的极小值问题并不好解，只能用迭代方法求解。求J极小值的常用迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤如下：(1)确定初始的$0值。对于60中的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)给出1nN2Je(k)2k土1并计算nN二2e2(k)N书(2.26)J一,2J式中W(YNJ:二马12(八03二2二22二工121nN2；

10、二二一e(k)Nkn1nN=N2k/(k)1nN2J=e2(k)2kzn1仃2越小越好，因为当方差小22灯取小时，e(k)取小，即残差取小。因此希望minJN(2.20)(2.21)(2.22)a2的估值取最(2.23)(3)计算J的梯度。9和海赛矩阵,有:2式中T-e(k)二；:e(k).?e(k)fe(k).Fe(k)?e(k),鼻ajan.:b0.:bnjcny(k)a1y(k-1)-any(k-n)-bou(k)biu(k-1)1“-bnu(k-n)-血ce(k-1)-Cne(k-n)同理可得由e(kj)求偏导，故：e(k-j)_；:e(k)z-j汨汨nN=e(k)k=n1fe(k)c

11、9(2.27).:e(k)”)一03。肆自一-卫川H 汨函(2.28)；:e(k)ja；ny(k-J).二Cjjife(k-j)(2.29)；：e(k)：:bin=-u(k-i)-，Cjjife(kj)(2.30):：e(k)n沦(k-j)二-e(k7Cj：Cijd二G将式(2.29)移项化简，有y(k.i)=H1,星厂Cj网 0ji蜘j=o二ai因为e(k-j)=e(k)z-j(2.31)(2.32)(2.33)(2.34)(2.35)将(2.34)代入(2.32),所以ygi)CjJ)Cjk:四zTj=o-a；j-o、a-a；j=oC(zA=1C|Z,Cnz所以得c(z)ek)=y(k-i

12、)(2.36).:3i同理可得(2.30)和(2.31)为c(z，)史k)=u(ki)(2.37)由c(zJ1)-(-k)=-e(k-i)(2.38)根据(2.36)构造公式c(z1)邳邳_(i_j)=yk(ij)j=y(k-i)(2.39)同将其代入(2.36),可得c(z严严k-(i-j)=c(z)四(2.40)二aj二ai消除c(z)可得沟(k)；:e(k-ij)沦(k-i1)：:aiFj：为1同理可得(2.37)和(2.38)式fe(k)_；:e(k-ij)_；:e(k-i)能由fb。：:e(k)_fe(kij)_：:e(k-i1)Xi论沁式(2.29)、式(2.30)和式(2.31)

13、均为差分方程，这些差分方程的初始条件为0,可通过求解这些差分方程，分别求出e(k)关于ai：.,a,b0,bn,G,g的全部偏导数，而这A当8接近于真值8时，e(k)接近于0。在这种情况下，式(2.44)等号右边第2项接近于0,f2J-J可近似表不为(2.41)(2.42)(2.43)些偏导数分别为y(k),u(k)和e(k)的线性函数。下面求关于T三_nCe(k)23nN，二e(k)kR1： 2e(k)26的二阶偏导数，即(2.44)3_曾阳k)阳k)【Jk更1F_F(2.45)tI茜:2J_则利用式(2.45)计算J2比较简单。:2按牛顿-拉卜森计算0的新估值01,有01=00-1()-1

14、(2.46)I的怎心重复(2)至(4)的计算步骤，经过r次迭代计算之后可得$，近一步迭代计算可得丁.I/2JjFr1-1r-()一(2.47)如果71-；-r4：10(2.48)则可停止计算，否则继续迭代计算。式(2.48)表明，当残差方差的计算误差小于0.01%时就停止计算。这一方法即使在噪声比较大的情况也能得到较好的估计值0。三实验内容设SISO系统差分方程为y(k)a1y(k-1)a2y(k。2)=b1u(k。1)b2u(k-2)(k)其中极大似然估计法默认e(k)为：e(k)=C(z,)(k)=(k)G(k1)C2(k-2)辨识参数向量为8=a1a2bib2c1c2T式中，？k)为噪声

15、方差各异的白噪声或有色噪声;u(k)和y(k)表示系统的输入输出变量。四实验流程图阳s.sRM1.算法的程序流程图五代码实现程序如下：U=1,1470.201-0.787-1.584-1.0520.8661.1521.5730.6260.433-0.9580.810-0.0440.947-1.474-0.719-0.0861.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.433-1.177-0.390-0.9821.435-0.119-0.769-0.8990.882-1.008-0.8440.628-0.6791.5411.375-0.

16、984-0.5821.6090.090-0.813-0.428-0.848-0.4100.048-1.099-1.1080.259-1.627-0.5280.2031.2041.691-1.235-1.228-1.2670.3090.0430.0431.4611.5850.552-0.601-0.3190.7440.829-1.626-0.127-1.578-0.8221.469-0.379-0.2120.1780.493-0.056-0.12941.228-1.606-0.382-0.2290.313-0.161-0.810-0.2770.983-0.2880.8461.3250.7230.

17、7130.6430.4630.7861.1610.850-1.349-0.5961.5120.795-0.7130.453-1.6040.889-0.9380.0560.829-0.981-1.2321.327-0.6810.114-1.1351.284-1.2010.7580.590-1.0070.3900.836-1.52-1.053-0.0830.6190.840-1.258-0.3540.629-0.2421.680-1.236-0.8030.537-1.1001.417-1.0240.6710.688-0.123-0.9520.232-0.793-1.1381.1540.2061.1

18、961.0131.518-0.553-0.9870.167-1.4450.6301.2550.311-1.7260.9751.7181.3601.6671.1111.0180.078-1.665-0.760.1.184-0.6140.994-0.0890.9471.706-0.3951.222-1.3510.2311.4250.114-0.689-0.7041.0700.2621.6101.489.-1.6020.020-0.601-0.020-0.601-0.2351.2451.226-0.2040.926-1.297%输入数据Y=0.0862.2100.486-1.812-3.705-2.

19、6881.5772.8833.7051.6420.805-2.0880.946-0.0391.984-2.545-1.727-0.2312.4403.5832.9151.4433.5980.7022.6383.611-320.6662.377-0.554-2.0882.6980.189-1.633-2.0101.716-1.641-1.8851.061-0.9682.9113.088-1.629-1.5333.0300.614-1.483-1.029-1.948-1.066-0.113-2.144-2.6260.134-3.043-1.3410.3382.7023.813-1

20、.924-2.813-1.7953.0021.0271.0272.7553.5841.737-0.837-0.6171.703.2.045-2.886-0.542-2.991-1.8593.0450.068-0.3750.4511.0360.153-0.4742.512-2.681-0.954-0.3070.628-0.270- 0.2770.983-0.2880.8461.3250.7231.7501.4011.340.0.9161.3962.4462.1032.432-1.4863.0312.373-0.763.- 0.752-3.2071.385-1.642-0.1181.756-1.6

21、13-1.6902.136- 1.136-0.005-2.2102.331-2.2040.9831.347-1.6910.5951.809-2.204-2.330-0.4541.2902.080-1.990-0.7701.240-0.2523.137-2.3791.2061.221-1.9772.471-1.6801.1481.8160.055-1.8560.269-1.323-2.4861.9580.8232.4812.2093.167-0.762-2.225-0.123-2.7861.0262.8431.071- 3.3171.5143.8073.3883.683-1.935-1.9350

22、.309-3.390- 2.1242.192-0.855-1.6560.0161.8043.774-0.0592.371- 2.322-0.0322.6320.565-1.460-1.8391.9170.8653.1803.261-2.755-0.536-1.171-0.905-3.303-0.8342.4903.0390.1341.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=

23、zeros(nc,1);%vf(k-i)thetae_1=zeros(na+nb+1+nc,1);%参数估计初值P=eye(na+nb+1+nc);fork=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

24、+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);fori=nc:-1:2xiek(i)=xiek(i-1);xiefk(i)=xiefk(i-1);endxiek(1)=xie;xiefk(1)=xief;fori=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(0L-22);figure(2)plot(1:L,thetae(