1、系統(tǒng)辨識大作業(yè)最小二乘法及其相關(guān)估值方法應(yīng)用學(xué)院：自動化學(xué)院專業(yè)：信息工程學(xué)號：姓名：馬志強(qiáng)日期：2010.11.14基于最小二乘法的多種系統(tǒng)辨識方法研究1. 最小二乘法的引出在系統(tǒng)辨識中用得最廣泛的估計(jì)方法是最小二乘法(LS)。設(shè)單輸入-單輸出線性定長系統(tǒng)的差分方程為xk+a1xk-1+ank-n=b0uk+bnuk-n，k=1，2，3，(5.1.1)式中: uk為隨機(jī)干擾；xk為理論上的輸出值。xk只有通過觀測才能得到，在觀測過程中往往附加有隨機(jī)干擾。xk的觀測值yk可表示為yk=xk+nk(5.1.2)式中: nk為隨機(jī)干擾。由式(5.1.2)得xk=yk-nk(5.1.3)將式(5.1

2、.3)帶入式(5.1.1)得yk+a1yk-1+anyk-n=b0uk+b1uk-1+bnuk-n+nk+i=1nai(k-i)(5.1.4)我們可能不知道nk的統(tǒng)計(jì)特性，在這種情況下，往往把nk看做均值為0的白噪聲。設(shè)(k)=n(k)+i=1nai(k-i)(5.1.5)則式(5.1.4)可寫成yk=-a1yk-1-a2yk-2-anyk-n+b0uk+b1uk-1+bnuk-n+(k)(5.1.6)在觀測uk時(shí)也有測量誤差，系統(tǒng)內(nèi)部也可能有噪聲，應(yīng)當(dāng)考慮它們的影響。因此假定(k)不僅包含了xk的測量誤差，而且包含了uk的測量誤差和系統(tǒng)內(nèi)部噪聲。假定(k)是不相關(guān)隨機(jī)序列(實(shí)際上(k)是相關(guān)

3、隨機(jī)序列)?，F(xiàn)分別測出n+N個(gè)隨機(jī)輸入值y1，y2，yn+N，u1，u2，u(n+N)，則可寫成N個(gè)方程，即yn+1=-a1yn-a2yn-1-any1+b0un+1+b1un+bnu1+(n+1)yn+2=-a1yn+1-a2yn-any2+b0un+2+b1un+1+bnu2+(n+2)yn+N=-a1yn+N-1-a2yn+N-2-anyN+b0un+N+b1un+N-1+bnuN+(n+N)上述N個(gè)方程可寫成向量-矩陣形式y(tǒng)(n+1)y(n+2)y(n+N)=-y(n)-y(n+1)-y(1)u(n+1)u(1)-y(2)u(n+2)u(2) -y(n+N-1)-y(N)u(n+N)u

4、(N)a1anb0bn+(n+1)(n+2)(n+3)(5.1.7)設(shè)y=y(n+1)y(n+2)y(n+N)，=a1anb0bn，=(n+1)(n+2)(n+N)=-y(n)-y(n+1)-y(1)u(n+1)u(1)-y(2)u(n+2)u(2) -y(n+N-1)-y(N)u(n+N)u(N)則式(5.1.7)可寫為y=+(5.1.8)式中：y為N維輸出向量；為N維噪聲向量；為(2n+1)維參數(shù)向量；為N(2n+1)測量矩陣。因此式(5.1.8)是一個(gè)含有(2n+1)個(gè)未知參數(shù)，由N個(gè)方程組成的聯(lián)立方程組。如果N2n+1，方程數(shù)少于未知數(shù)數(shù)目，則方程組的解是不定的，不能唯一地確定參數(shù)向量

5、。如果N=2n+1，方程組正好與未知數(shù)數(shù)目相等，當(dāng)噪聲=0時(shí)，就能準(zhǔn)確地解出=-1y(5.1.9)如果噪聲0，則=-1y-1(5.1.10)從上式可以看出噪聲對參數(shù)估計(jì)是有影響的，為了盡量較小噪聲對估值的影響。在給定輸出向量y和測量矩陣的條件下求系統(tǒng)參數(shù)的估值，這就是系統(tǒng)辨識問題?？捎米钚《朔▉砬蟮墓乐?，以下討論最小二乘法估計(jì)。2. 最小二乘法估計(jì)算法設(shè)表示 的最優(yōu)估值，y表示 y的最優(yōu)估值，則有y=(5.1.11)y=y(n+1)y(n+2)y(n+N)，=a1anb0bn寫出式(5.1.11)的某一行，則有yk=-a1yk-1-a2yk-2-anyk-n+b0uk+b1uk-1+bnuk

6、-n+k=-i=1naiy(k-i)+i=0nbiu(k-i)，k=n+1，n+2，n+N(5.1.12)設(shè)e(k)表示y(k)與y(k)之差，即ek= yk-yk=yk-yk=yki=1naiyk-i+i=0nbiuk-i=1+a1z-1+anz-nyk-b0+b1z-1+bnz-nuk=az-1yk-bz-1uk，k=n+1，n+2，n+N(5.1.13)式中az-1=1+a1z-1+anz-nbz-1=b0+b1z-1+bnz-nek成為殘差。把k=n+1，n+2，n+N分別代入式(5.1.13)可得殘差en+1，en+2，en+N。設(shè)e=en+1en+2en+NT則有ek= y-y=y

7、-(5.1.14)最小二乘估計(jì)要求殘差的平方和為最小，即按照指數(shù)函數(shù)J= eTe=(y-)T(y-)(5.1.15)為最小來確定估值。求J對的偏導(dǎo)數(shù)并令其等于0可得J=-2Ty-=0(5.1.16)T=Ty(5.1.17)由式(5.1.17)可得的最小二乘估計(jì)=(T)-1T y(5.1.18)3.遞推最小二乘法為了實(shí)現(xiàn)實(shí)時(shí)控制，必須采用遞推算法，這種辨識方法主要用于在線辨識。設(shè)已獲得的觀測數(shù)據(jù)長度為N，將式(5.1.8)中的y和分別用YN，N，N來代替，即YN=N+N(5.3.1)用N表示 的最小二乘估計(jì)，則N=(NTN)-1NT YN(5.3.2)設(shè)PN=（NTN）-1(5.3.5)于是N=

8、PNNTYN(5.3.6)如果再獲得1組新的觀測值un+N+1和y(n+N+1)，則又增加1個(gè)方程yN+1=N+1T+N+1(5.3.7)式中yN+1=yn+N+1，N+1=(n+N+1)N+1T=-yn+N -yN+1 un+N+1 u(N+1)將式(5.3.1)和式(5.3.7)合并，并寫成分塊矩陣形式，可得YNyN+1=NN+1T +NN+1(5.3.8)根據(jù)上式可得到新的參數(shù)估值N+1=NN+1TTNN+1T-1NN+1TTYNyN+1=PN+1NN+1TTYNyN+1=PN+1(NTYN+N+1yN+1)(5.3.9)式中PN+1=NN+1TTNN+1T-1=（NTN+N+1N+1T

9、）-1根據(jù)矩陣求逆引理可以求得遞推最小二乘法辨識公式N+1=N+KN+1(yN+1-N+1TN)(5.3.19)KN+1=PN+1N+1(1+N+1TPNN+1)-1(5.3.20)PN+1=PN-PNN+1(1+N+1TPNN+1)-1N+1TPN(5.3.21)由于進(jìn)行遞推計(jì)算需要給出N和PN初值P0和0，通過計(jì)算證明，可以取初值：0=0，P0=c2I，c是充分大的常數(shù)，I為2n+1(2n+1)單位矩陣，則經(jīng)過若干次遞推之后能夠得到較好的參數(shù)估計(jì)。3. 輔助變量法輔助變量法是一種可克服最小二乘有偏估計(jì)的一種方法，對于原辨識方程y=+(5.4.1)當(dāng)(k)是不相關(guān)隨機(jī)序列時(shí)，最小二乘法可以得

10、到參數(shù)向量的一致無偏估計(jì)。但是，在實(shí)際應(yīng)用中(k)往往是相關(guān)隨機(jī)序列。假定存在著一個(gè)(2n+1)N的矩陣Z滿足約束條件limN1NZT=0limN1NZT=Q(5.4.2)式中Q是非奇異的。用ZT乘以式(5.4.1)等號兩邊得ZTy=ZT+ZT (5.4.3)由上式得=(ZT)-1ZTy-(ZT)-1ZT(5.4.4)如果取IV=(ZT)-1ZTy(5.4.5)作為估值，則稱IV為輔助變量估值，矩陣Z成為輔助變量矩陣，Z中的元素稱為輔助變量。常用的輔助變量法有遞推輔助變量參數(shù)估計(jì)法，自適應(yīng)濾波法，純滯后等。4. 廣義最小二乘法廣義最小二乘法是能克服最小二乘法有偏估計(jì)的另一種方法，這種方法計(jì)算比

11、較復(fù)雜但效果比較好。下面直接介紹廣義最小二乘法的計(jì)算步驟：(1)應(yīng)用得到的輸入和輸出數(shù)據(jù)u(k)和yk(k=1，2，3，n+N)，按模型az-1yk=bz-1uk+k求出 的最小二乘估計(jì)(1)=a1(1)an(1)b0(1)bn(1)(2)計(jì)算殘差e1(k)e1k=a(1)z-1yk-b(1)z-1uk(3)用殘差e1k代替 k，計(jì)算f(1)=(1)T(1)-1(1)Te1(4)計(jì)算y1(k)和u1(k)y1k=yk+f11yk-1+fm1y(k-m)u1k=uk+f11uk-1+fm1u(k-m)(5)應(yīng)用得到的y1(k)和u1(k)按模型az-1y1k=bz-1u1k+(k)用最小二乘法重

12、新估計(jì)，得到的第2次估值(2)。然后按步驟(2)計(jì)算殘差e2(k)，按步驟(3)重新估計(jì)f，得到估值f 2。再按照步驟(4)計(jì)算y2(k)和u2k，按照步驟(5)求的第3次估值(3)。重復(fù)上述循環(huán)，之道的估值(i)收斂為止。5. 一種交替的廣義最小二乘法求解技術(shù)(夏式法)這種方法是夏天長提出來的，又稱夏式法。以上討論過的廣義最小二乘法的特點(diǎn)在于系統(tǒng)的輸入和輸出信號反復(fù)過濾。一下介紹的夏式法是一種交替的廣義最小二乘法求解技術(shù)，它不需要數(shù)據(jù)反復(fù)過濾，因而計(jì)算效率較高。這種方法可消去最小二乘估計(jì)中的偏差，而且由這種方法導(dǎo)出的計(jì)算方法也比較簡單?；谝陨系膸追N方法，有y=+=f+(5.7.1)因而有y

13、=+f+=f+(5.7.2)應(yīng)用最小二乘法可得到參數(shù)估值f=TT-1TTy(5.7.3)可以推出=(T)-1T y-(T)-1f(5.7.11)上式中的第1項(xiàng)是最小二乘估計(jì)LS，第2項(xiàng)是偏差項(xiàng)B，所以必須準(zhǔn)確計(jì)算B。為了準(zhǔn)確計(jì)算B，可采用迭代的方法。6. 專題解答設(shè)但輸入-單輸出系統(tǒng)的差分方程為yk=-a1yk-1-a2yk-2+b1uk-1+b2uk-2+(k)k=k+a1k-1+a2k-2取真實(shí)值T=a1a2b1b2=1.6420.7150.390.35，輸入數(shù)據(jù)如下所示ku(k)ku(k)ku(k)11.14711-0.958210.48520.201120.810221.6333-0.

14、78713-0.044230.0434-1.159140.947241.3265-1.05215-1.474251.70660.86616-0.71926-0.34071.15217-0.086270.89081.57318-1.099281.14490.626191.450291.177100.433201.15130-0.390用的真實(shí)值利用查分方程求出yk作為測量值，k為均值為0，方差為0.1，0.5的不相關(guān)隨機(jī)序列。(1) 用最小二乘法估計(jì)參數(shù)T=a1a2b1b2。(2) 用遞推最小二乘法估計(jì)T=a1a2b1b2。(3) 用輔助變量法估計(jì)參數(shù)T=a1a2b1b2。(4) 設(shè)k+2k=k

15、，用廣義最小二乘法估計(jì)參數(shù)T=a1a2b1b2。(5) 用夏式法估計(jì)參數(shù)T=a1a2b1b2。(6) 詳細(xì)分析和比較所獲得的參數(shù)辨識結(jié)果，并說明上述參數(shù)便是方法的優(yōu)缺點(diǎn)。根據(jù)題目要求的解法，利用Matlab編程實(shí)現(xiàn)系統(tǒng)辨識的估值利用最小二乘法估計(jì)的結(jié)果如下：最小二乘法方差a1a2b1b20.00011.62800.70280.39713.4491.63160.70590.39470.34941.63540.71200.39180.34631.63620.70820.39700.35271.63600.71650.39060.35341.62890.70460.39080.34370.0011.

16、15430.67660.40640.33041.55770.63710.38680.32491.60500.68600.37370.32441.60600.68160.35830.31671.61950.70300.39070.33661.56700.65720.37520.31400.011.35380.50100.34860.17090.89560.16370.42370.06971.00080.20360.42360.12681.34030.47070.38260.26151.05740.22890.36820.13171.12310.29630.35920.15060.11.14240

17、.27100.32840.22161.02550.17360.37660.18440.88960.11090.38930.08130.81820.11140.42980.09230.81000.01530.41220.13520.77150.13110.47140.06400.50.97510.10170.02710.07920.89380.07400.34840.16340.19270.01970.37620.15210.55060.03920.65100.05840.75600.04940.33720.14370.94590.13770.38150.1853部分程序運(yùn)行結(jié)果遞推最小二乘法方

18、差a1a2b1b20.00011.67540.67870.52070.39241.30760.29000.00750.32841.51460.69631.14010.16391.57330.77820.41490.71101.16020.47530.67360.31591.20910.31920.52770.02750.0011.47670.40400.26790.55121.62590.75940.32530.37801.53930.47570.12680.43461.15480.17000.19260.82510.88580.07600.33850.04061.41290.31270.09

19、920.83800.011.34850.34450.31940.37101.16390.32960.78130.21621.99461.23231.48520.03041.39240.35430.33190.45721.39820.36080.77730.31521.63460.72290.57800.39470.11.56240.71320.44220.41121.73350.71520.08440.63991.47630.53660.32550.31161.44770.34890.22180.22651.62160.70820.65950.42751.51050.40000.01130.2

20、2130.51.79270.94110.27300.34711.55560.88770.59720.12171.78681.25381.12480.21001.57330.74340.35890.13871.31930.60841.29710.30291.59590.53860.01410.6947部分程序運(yùn)行結(jié)果：輔助變量法方差a1a2b1b20.00011.77990.85880.41470.40461.30760.29000.00750.32841.67350.75780.40600.34841.58120.65460.37710.36111.66570.74690.37720.3561

21、1.52810.65090.36450.33020.0011.62950.67750.40820.35011.64250.73050.39370.35431.55950.60520.35630.33621.41450.49250.40210.28731.63710.72700.36810.34181.35390.47330.39060.24890.011.34510.47000.38220.26801.36570.48930.43490.24961.37020.50090.43880.26111.18840.37070.34570.15211.36360.53300.41350.20641.3

22、1580.48360.45620.21690.11.55450.61670.41040.35121.59000.66480.40980.36641.66100.70290.40010.34951.51040.60140.39680.32561.56200.64960.39500.30971.44180.57060.41830.27830.51.49520.57040.37690.34831.55920.65410.43870.33301.36370.53020.38650.19161.55430.63240.32080.29291.53850.59940.36610.35761.35110.4

23、8390.34120.2208部分程序運(yùn)行結(jié)果：廣義最小二乘法方差a1a2b1b20.00011.64510.71820.38950.35031.64680.72050.38660.34921.63080.70380.38520.34571.64670.71880.38980.35441.64110.71560.38980.34771.64410.71850.38960.34880.0011.70270.76780.40590.37611.64770.71660.39730.36431.63910.71130.38360.35041.65860.71660.37910.36721.68310.

24、74720.39940.37311.59640.68500.36550.32310.011.68260.74690.38590.36001.72450.76960.35650.36421.65770.71860.39100.38201.66560.72630.34250.34131.69030.73920.36330.38381.69420.74390.38340.38800.10.89960.15150.44350.06041.03800.30440.34100.04641.65710.65690.45890.42591.21800.45840.31520.10300.19570.09740

25、.21640.10351.68590.72320.41440.44040.51.40530.67030.36720.13721.57830.61480.35120.41731.02330.06580.20350.40601.59730.81030.35440.21311.61910.67090.31110.48661.08040.22100.75000.5185部分程序運(yùn)行結(jié)果：夏式法方差a1a2b1b20.00011.63670.70900.39100.35111.64510.71840.38760.34721.36180.21590.82050.58511.62910.70750.3911

26、0.34641.61730.69290.39550.34391.63070.70660.39120.34490.0011.53140.64070.40160.30921.56320.65750.42040.33191.61720.69460.39250.32581.56700.64890.40160.32331.53260.62430.39950.30701.62510.69460.36210.33600.010.88150.10890.37340.04401.37630.49720.41210.26541.50270.59520.38540.32771.34210.42680.39780.2

27、7901.24540.30240.36690.25251.33610.47410.49990.29000.10.66390.10740.68470.08271.14060.31540.42590.23670.17910.29160.94540.21261.20680.38270.74600.32200.64120.07040.58020.08771.03840.28780.74930.37930.50.93300.07520.97410.58820.62910.16860.46790.78690.59860.27920.26050.20080.44410.80330.74421.76670.6

28、2160.14400.58210.49060.65760.09991.00020.0263部分程序運(yùn)行結(jié)果：結(jié)論：通過編程計(jì)算，獲得在噪聲方差比較小的情況下，各種方法所獲得的估值比較理想，但隨著噪聲方差的增大，估值的偏差隨之增大，橫向比較看來夏式法與廣義最小二乘法能夠更好地還原參數(shù)值，當(dāng)觀測值足夠多時(shí)，各種方法都能很好地反映參數(shù)真實(shí)值。Matlab源程序：%最小二乘估計(jì)% clear u= 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719

29、-0.086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177 -0.390; n=normrnd(0, sqrt(0.1), 1, 31); z=zeros(1,30);for k=3:31 z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715*n(k-2);end h0=-z(2) -z(1) u(2) u(1); HLT=h0,zeros(4,28);for k=3:30 h1=-z(k

30、) -z(k-1) u(k) u(k-1); HLT(:,k-1)=h1;end HL=HLT; y=z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16);z(17);z(18);z(19);z(20);z(21);z(22);z(23);z(24);z(25);z(26);z(27);z(28);z(29);z(30);z(31);%求出FAI c1=HL*HL; c2=inv(c1); c3=HL*y; c=c2*c3;%display(方差=0.1時(shí)，最小二乘法估計(jì)辨識參數(shù)如下：); a1

31、=c(1); a2=c(2); b1=c(3); b2=c(4); clear%遞推最小二乘法估計(jì) u= 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.985 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.144 1.177 -0.390; z(2)=0; z(1)=0;n=normrnd(0, sqrt(0.1), 1, 31);for k=3:3

32、1 z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715*n(k-2);end c0=0.001 0.001 0.001 0.001; %直接給出被辨識參數(shù)的初始值,即一個(gè)充分小的實(shí)向量 p0=106*eye(4,4); %直接給出初始狀態(tài)P0，即一個(gè)充分大的實(shí)數(shù)單位矩陣 E=0.; %取相對誤差E=0. c=c0,zeros(4,30); %被辨識參數(shù)矩陣的初始值及大小 e=zeros(4,30); %相對誤差的初始值及大小for k=3:30; %開始求K h1=-z(k-1),-z(

33、k-2),u(k-1),u(k-2); x=h1*p0*h1+1; x1=inv(x); %開始求K(k) k1=p0*h1*x1;%求出K的值 d1=z(k)-h1*c0; c1=c0+k1*d1; %求被辨識參數(shù)c e1=c1-c0; %求參數(shù)當(dāng)前值與上一次的值的差值 e2=e1./c0; %求參數(shù)的相對變化 e(:,k)=e2; %把當(dāng)前相對變化的列向量加入誤差矩陣的最后一列 c0=c1; %新獲得的參數(shù)作為下一次遞推的舊參數(shù) c(:,k)=c1; %把辨識參數(shù)c 列向量加入辨識參數(shù)矩陣的最后一列 p1=p0-k1*k1*h1*p0*h1+1; %求出 p(k)的值 p0=p1; %給下

34、次用 if e2=E break; %如果參數(shù)收斂情況滿足要求，終止計(jì)算 endend%display(方差為0.0001遞推最小二乘法辨識后的結(jié)果是：); a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:);%display(a1,a2,b1,b2經(jīng)過遞推最小二乘法辨識的結(jié)果是：);for i=3:31; if(c(1,i)=0) q1=c(1,i-1); break; endendfor i=3:31; if(c(2,i)=0) q2=c(2,i-1); break; endendfor i=3:31; if(c(3,i)=0) q3=c(3,i-1);

35、break; endendfor i=3:31; if(c(4,i)=0) q4=c(4,i-1); break; endend a1=q1; a2=q2; b1=q3 ; b2=q4;%clear%輔助變量遞推最小二乘法估計(jì) na=2; nb=2; siitt=1.642 0.715 0.39 0.35; siit0=0.001*eye(na+nb,1); p=106*eye(na+nb); siit(:,1)=siit0; y(2)=0;y(1)=0; x(1)=0;x(2)=0; j=0; u= 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1

36、.573 0.626 0.433 -0.985 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.144 1.177 -0.390; n=normrnd(0, sqrt(0.01), 1, 31);for k=3:31; h=-y(k-1),-y(k-2),u(k-1),u(k-2); y(k)=h*siitt+n(k)+1.642*n(k-1)+0.715*n(k-2); hx=-x(k-1),-x(k-2),u(k-1),u(k-2)

37、; kk=p*hx/(h*p*hx+1); p=eye(na+nb)-kk*h*p; siit(:,k-1)=siit0+kk*y(k)-h*siit0; x(k)=hx*siit(:,k-1); j=j+(y(k)-h*1.642 0.715 0.39 0.35)2; e=max(abs(siit(:,k-1)-siit0)./siit0); ess(:,k-2)=siit(:,k-1)-siitt; siit0=siit(:,k-1); end a1=siit0(1); a2=siit0(2); b1=siit0(3); b2=siit0(4); clear%廣義最小二乘估計(jì) clear;

38、 nn = normrnd(0,sqrt(0.5),1,31); uk=1.147 0.201 -0.787 -1.589 -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 1.144 1.177 -0.390; yk(1)=0; yk(2)=0;for i=1:29; yk(i+2)=-1.642*yk(i+1)-0.715*yk(i)+0.39*uk(

39、i+1)+0.35*uk(i)+nn(i+2)+1.642*nn(i+1)+0.715*nn(i);end;for i=1:29; A(i,:)=-yk(i+1) -yk(i) uk(i+1) uk(i);end siit=inv(A*A)*A*(yk(3:31)+nn(2:30); e(1)=yk(1); e(2)=yk(2)+siit(1)*yk(1)-siit(3)*uk(1);for i=3:31; e(i)=yk(i)+siit(1)*yk(i-1)+siit(2)*yk(i-2)-siit(3)*uk(i-1)-siit(4)*uk(i-2);endfor i=1:29; fai(

40、i,:)=-e(i+1) -e(i);end f=inv(fai*fai)*fai*e(3:31);for i=3:31; yk(i)=yk(i)+f(1)*yk(i-1)+f(2)*yk(i-2);end yk(2)=yk(2)+f(1)*yk(1);for i=3:30; uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);end uk(2)=uk(2)+f(1)*uk(1);for j=1:30for i=1:29; A(i,:)=-yk(i+1) -yk(i) uk(i+1) uk(i);end siit=inv(A*A)*A*yk(3:31); e(1)=yk

41、(1); e(2)=yk(2)+siit(1)*(yk(1)-siit(3)*uk(1);for i=3:31; e(i)=yk(i)+siit(1)*(yk(i-1)+siit(2)*(yk(i-2)-siit(3)*uk(i-1)-siit(4)*uk(i-2);endfor i=1:29; fai(i,:)=-e(i+1) -e(i);end f=inv(fai*fai)*fai*e(3:31); k1(j)=f(1); k2(j)=f(2);for i=3:31; yk(i)=yk(i)+f(1)*(yk(i-1)+f(2)*(yk(i-2); endyk(2)=yk(2)+f(1)*yk(1);for i=3:30 uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);end uk(2)=uk(2)+f(1)*uk(1);end siit;%用夏氏偏差修正法估計(jì)參數(shù) clear; u= 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.985 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.144 1.1