用極大似然法進(jìn)行全參數(shù)估計(jì)

1、實(shí)用文檔北京工商大學(xué)系統(tǒng)建模與辨識(shí)課程上機(jī)實(shí)驗(yàn)報(bào)告（2016年秋季學(xué)期）專業(yè)名稱：控制工程上機(jī)題目：用極大似然法進(jìn)行參數(shù)估計(jì)專業(yè)班級(jí)：計(jì)研3班學(xué)生姓名： 王瑤 吳超學(xué) 號(hào)：10011316259 10011316260指導(dǎo)教師：劉翠玲2017 年1 月實(shí)驗(yàn)?zāi)康耐ㄟ^(guò)實(shí)驗(yàn)掌握極大似然法在系統(tǒng)參數(shù)辨識(shí)中的原理和應(yīng)用。實(shí)驗(yàn)原理1極大似然原理設(shè)有離散隨機(jī)過(guò)程VJ與未知參數(shù)日有關(guān)，假定已知概率分布密度 fM|e)。如果我們 得到n個(gè)獨(dú)立的觀測(cè)值 Vi,V2,M,則可得分布密度 f(Vi|e), fMS,，f(VnB)。 要求根據(jù)這些觀測(cè)值來(lái)估計(jì)未知參數(shù)e，估計(jì)的準(zhǔn)則是觀測(cè)值 vk的出現(xiàn)概率為最大。為此，定

2、義一個(gè)似然函數(shù)LM,V2,Vn|8) = f(Vi|8)f(V2|8fM|8)(1.1)上式的右邊是n個(gè)概率密度函數(shù)的連乘，似然函數(shù)L是H的函數(shù)。如果L達(dá)到極大值，VkA的出現(xiàn)概率為最大。因此，極大似然法的實(shí)質(zhì)就是求出使L達(dá)到極大值的日的估值8。為了A便于求9,對(duì)式(1.1)等號(hào)兩邊取對(duì)數(shù)，則把連乘變成連加，即 nln L = £ ln f (Vi 日)(1.2 )由于對(duì)數(shù)函數(shù)是單調(diào)遞增函數(shù)，當(dāng) 對(duì)8的偏導(dǎo)數(shù)，令偏導(dǎo)數(shù)為 0,可得ln LA解上式可得日的極大似然估計(jì)8 ML。i工L取極大值時(shí)，InL也同時(shí)取極大值。求式(1.2)=0(1.3 )2系統(tǒng)參數(shù)的極大似然估計(jì)Newton-R

3、aphson法實(shí)際上就是一種遞推算法，可以用于在線辨識(shí)。不過(guò)它是一種依每 L 次觀測(cè)數(shù)據(jù)遞推一次的算法，現(xiàn)在我們討論的是每觀測(cè)一次數(shù)據(jù)就遞推計(jì)算一次參數(shù)估計(jì)值 得算法。本質(zhì)上說(shuō)，它只是一種近似的極大似然法。設(shè)系統(tǒng)的差分方程為a(z")y(k) = b(z,)u(k)十：(k)(2.1)式中a(zJ) = 1 4z)'anzb(z 工)=b0 ' b1z 1 bnz因?yàn)楫a(chǎn)(k)是相關(guān)隨機(jī)向量，故(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的高斯分布白噪聲序列。多項(xiàng)式a(z,),b(z,)和c(z,)中的系數(shù)a1,a,b0,bn,c1,cn和序列e(k)的均方差仃都是未知參數(shù)。設(shè)待估參數(shù)日=aan bobn Gg】(23)并設(shè)y(k)的預(yù)測(cè)值為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)為預(yù)測(cè)誤差；1，bi, C；為ai,可，ci的估值。預(yù)測(cè)誤差可表示為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)因此預(yù)測(cè)誤差ie(k)滿足關(guān)系式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假定預(yù)測(cè)誤差e(k)服從均值為0的高斯分布，并設(shè)序列 Q(k)具有相同的方差 仃2。因 為Q(k)與c(z"), a(z)和b(z,)有關(guān)，所以仃2是被估參數(shù)日的函數(shù)。為了書寫方便， 把式(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 )或?qū)懗?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個(gè)方程式，把這 N個(gè)方程式寫成向量-矩陣形式式中-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)因?yàn)橐鸭俣?

8、L(k)是均值為0的高斯噪聲序列，高斯噪聲序列的概率密度函數(shù)為f =(27：C 212rejp-口(y -m) )2(2.14 )式中y為觀測(cè)值，er2和m為y的方差和均值，那么,1r 12八、(2.15 )f =Texp-0e(k)22 -1(2二二)2對(duì)于e(k)符合高斯噪聲序列的極大似然函數(shù)為L(zhǎng)(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對(duì)上式(2.17)等號(hào)兩邊取對(duì)數(shù)得，、，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 )或?qū)憺镹N 91 n N 9(2.19 )ln L(YN 二二)二一一ln2二 一一ln；-2 ' e2(k)222二百求帖1(丫；日,仃)對(duì)仃2的偏導(dǎo)數(shù)，令其等于 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越小越好，因?yàn)楫?dāng)方差 小22燈 取小時(shí)， e(k)取小，即殘差取小。因此希望minJN(2.20 )(2.21 )(2.22 )a2的估值取最(2.23 )因?yàn)槭?2.10 )可理解為預(yù)測(cè)模型，而 e(k)可看做預(yù)測(cè)誤差。因此使式(2.22)最小 就是使誤差的平方之和最小，即使對(duì)概率密度不作任何假設(shè)，這樣的準(zhǔn)則也是有意義的。 因 此可按J最小來(lái)求現(xiàn),a,b0,bn,。,cn的估計(jì)值。由于e(k)式參數(shù)a1,a,b0,bn,G,cn的線性函數(shù)，因此J是這些參數(shù)的二次型函數(shù)。求使ln L(Yn，。)最大的S，等價(jià)于在式

11、(2.10)的約束條件下求 S使J為最小。由于J對(duì)G是非線性的，因而求J的極小值問(wèn)題并不好解， 只能用迭代方法求解。 求J極小值的常用 迭代算法有拉格朗日乘子法和牛頓 -拉卜森法。下面介紹牛頓-拉卜森法。整個(gè)迭代計(jì)算步驟 如下：(1)確定初始的$0值。對(duì)于00中的a1,a, b0,bn可按模型e(k) =a(z )y(k) -b(z )u(k)(2.24)A用最小二乘法來(lái)求，而對(duì)于日0中的c1, cn可先假定一些值。(2)計(jì)算預(yù)測(cè)誤差A(yù)e(k) =y(k)-y(k)(2.25)給出1n N 2J e (k) 2k土 1并計(jì)算n N二2v e2(k)N書(2.26)J一,'2J(3)計(jì)算

12、J的梯度。9和海賽矩陣 -J ,有:2n N=e(k)k國(guó)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)移項(xiàng)化簡(jiǎn)，有y(k.i)=HCjieUcn Cj 網(wǎng)0 a j 1caij =0tai因?yàn)閑(k - j) = e(k)z-j(2.31 )(2.32 )(2.33 )(2.34 )(2.35 )由e(k j)求偏導(dǎo)，故：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 )根據(jù)(2.36 )構(gòu)造公式c(z1)邳 _(i _ j) = yk (i j) j = y(k -i)(2.39 )同將其代入(2.36 ),可得c(z嚴(yán)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,可通過(guò)求解這些差分方程，分別求出e(k)關(guān)于ai,.,a,b0,bn,G，cn的全部偏導(dǎo)數(shù)，而這些偏導(dǎo)數(shù)分別為y(k), u(k)和e(k)的線性函數(shù)。下面求關(guān)于T三 _ nC £e(k)kR1nN，二 e(k)kR 1：2e(k)26的二階偏導(dǎo)數(shù)，即(2.44 )A當(dāng)8接近于真值8時(shí)，e(k)接近于0。在這種情況下，式(2.44)等號(hào)右邊第2項(xiàng)接近于0,f2J-J可近似表不為3_曾陽(yáng)k)砥k) T 【J k n 1 F _ F(2.45 )tI 茜:2J

16、_ _則利用式(2.45 )計(jì)算J2-±比較簡(jiǎn)單。32按牛頓-拉卜森計(jì)算0的新估值01,有01 =窯1()- 1(2.46 )I的怎心重復(fù)(2)至(4)的計(jì)算步驟，經(jīng)過(guò)r次迭代計(jì)算之后可得 $，近一步迭代計(jì)算可得丁. I /2J 1F r 1 - 1 r -()一(2.47 )如果71 - ；- r4:：10(2.48 )則可停止計(jì)算，否則繼續(xù)迭代計(jì)算。式(2.48)表明，當(dāng)殘差方差的計(jì)算誤差小于0.01%時(shí)就停止計(jì)算。這一方法即使在噪聲比較大的情況也能得到較好的估計(jì)值0。三實(shí)驗(yàn)內(nèi)容設(shè)SISO系統(tǒng)差分方程為y(k) a1y(k-1) a2 y(k。2) = b1u(k。1) b2u(

17、k-2)(k)其中極大似然估計(jì)法默認(rèn)e(k)為：e(k) =C(z,)(k) = (k) g (k1) C2 (k -2)辨識(shí)參數(shù)向量為8=a1a2燈b2 c 1 c 2 T式中，？ k)為噪聲方差各異的白噪聲或有色噪聲;u(k)和y(k)表示系統(tǒng)的輸入輸出變量。四實(shí)驗(yàn)流程圖陽(yáng)s.s RM1.算法的程序流程圖五代碼實(shí)現(xiàn)程序如下：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 %輸入數(shù)據(jù)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% 輸出數(shù)據(jù)na=2;nb=1;nc=2;d=1;nn=max(na,nc);L=size(Y,2);xiek=zeros(nc,1); % 白噪聲估計(jì)初值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); % 參數(shù)估計(jì)初值P=eye(na+nb+1+nc);for k=3:L%構(gòu)造向量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;%遞推極大似然參數(shù)估計(jì)算法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)%更新數(shù)據(jù)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('參數(shù)彳計(jì) a');legend('a_1','a_2'); axis(0 L -2 2);figure(2)plot(1:L,thetae(