Matlab系統(tǒng)辨識仿真CH3,CH4,CH6程序注釋與剖析

1、第3章程序及注釋例3.3考慮仿真對象z(k)-1.5z(k-1)0.7z(k-2)，u(k-1)0.5u(k-2)v(k)(3.41)其中，v(k)是服從正態(tài)分布的白噪聲N(0,1)。輸入信號采用4階M序列，幅度為1。選擇如下形式的辨識模型(3.42)z(k)+az(k一1)+az(k一2)，bu(k一1)+bu(k一2)+v(k)1212設(shè)輸入信號的取值是從k=1到k=16的M序列,則待辨識參數(shù)0LS其中，被辨識參數(shù)0、觀測矩陣zL、HL的表達(dá)為0=(HtH)-iHtzLS式為LS/X0，a1a2，z，zz(4)，H，LSb1b2L_z(16)L-z(15)-z(2)-z(3)-z(1)-z

2、(2)-z(14u(2)u(3)u(1)u(2)(3.43)u(15)u(14)程序框圖如圖3.2所示。Matlab6.0仿真程序如下：圖3.2最小二乘一次完成算法程序框圖%二階系統(tǒng)的最小二乘一次完成算法辨識程序，在光盤中的文件名：FLch3LSeg1.mu=-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1;%系統(tǒng)辨識的輸入信號為一個周期的M序列z=zeros(1,16);%定義輸出觀測值的長度fork=3:16z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%用理想輸出值作為觀測值endsubplot(3,1,1)%畫三行一列

3、圖形窗口中的第一個圖形stem(u)%畫輸入信號u的徑線圖形subplot(3,1,2)%畫三行一列圖形窗口中的第二個圖形i=1:1:16;%橫坐標(biāo)范圍是1到16，步長為1plot(i,z)%圖形的橫坐標(biāo)是采樣時刻i,縱坐標(biāo)是輸出觀測值z,圖形格式為連續(xù)曲線subplot(3,1,3)%畫三行一列圖形窗口中的第三個圖形stem(z),gridon%畫出輸出觀測值z的徑線圖形，并顯示坐標(biāo)網(wǎng)格u,z%顯示輸入信號和輸出觀測信號%L=14%數(shù)據(jù)長度HL=-z(2)-z(1)u(2)u(1);-z(3)-z(2)u(3)u(2);-z(4)-z(3)u(4)u(3);-z(5)-z(4)u(5)u(4

4、);-z(6)-z(5)u(6)u(5);-z(7)-z(6)u(7)u(6);-z(8)-z(7)u(8)u(7);-z(9)-z(8)u(9)u(8);-z(10)-z(9)u(10)u(9);-z(11)-z(10)u(11)u(10);-z(12)-z(11)u(12)u(11);-z(13)-z(12)u(13)u(12);-z(14)-z(13)u(14)u(13);-z(15)-z(14)u(15)u(14)%給樣本矩陣HL賦值LZL=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(

5、16)%給樣本矩陣zL賦值%CalculatingParameterscl=HL*HL;c2=inv(cl);c3=HL*ZL;c=c2*c3%計算并顯示gLS%DisplayParametersal=c(l),a2=c(2),bl=c(3),b2=c(4)%從g中分離出并顯示a1、a2、化、LS121b2%End程序運行結(jié)果：u=-1，1，-1，1，1，1，1，-1，-1，-1，1，-1，-1，1，1z=0，0，0.5000，0.2500，0.5250，2.1125，4.3012，6.4731，6.1988，3.2670，-0.9386，-3.1949，-4.6352，6.2165，-5.5

6、800，-2.5185001.0000-1.0000-0.50000-1.00001.0000-0.2500-0.50001.0000-1.0000-0.5250-0.25001.00001.0000-2.1125-0.52501.00001.0000-4.3012-2.11251.00001.0000-6.4731-4.3012-1.00001.0000-6.1988-6.4731-1.0000-1.0000-3.2670-6.1988-1.0000-1.00000.9386-3.26701.0000-1.00003.19490.9386-1.00001.00004.63523.1949-1

7、.0000-1.00006.21654.63521.0000-1.00005.58006.21651.00001.0000ZL=0.5000，0.2500，0.5250，2.1125，4.3012，HL=6.4731，6.1988，3.2670，-0.9386，-3.1949，-4.6352，-6.2165，-5.5800，-2.5185Tc=-1.5000，0.7000，1.0000，0.5000TV=54.3000,61.8000,72.4000,88.7000,118.6000,194.0000a1=-1.5000a2=0.7000b1=1.0000b2=0.5000圖3.3最小二乘一次

8、完成算法仿真實例中輸入信號和輸出觀測值10-1100-10100-10從仿真結(jié)果表3.1可以看出，由于所用的輸出觀測值沒有任何噪聲成分，所以辨識結(jié)果也無任何誤差。表3.1最小二乘一次完成算法的辨識結(jié)果參數(shù)a102b真值-1.50.71.00.5估計值-1.507100.5例3.4根據(jù)熱力學(xué)原理，對給定質(zhì)量的氣體，體積V與壓力P之間的關(guān)系為PVa，其中a和卩為待定參數(shù)。經(jīng)實驗獲得如下一批數(shù)據(jù)，V的單位為立方英寸，P的單位為帕每平方英寸。V54.361.872.488.7118.6194.0P61.249.537.628.419.210.1試用最小二乘一次完成算法確定參數(shù)a和卩。首先要寫出系統(tǒng)的最

9、小二乘表達(dá)式。為此，把體積V與壓力P之間的關(guān)系PVa，改為對數(shù)關(guān)系，即，logP=alogV+log卩。此式與式(3.14)，z(k)，h(k)0+e(k)，對比可得：z(k)，logP，h(k)，logV1，0，alogF。例3.4的Matlab6.0程序如下。%實際壓力系統(tǒng)的最小二乘辨識程序，在光盤中的文件名：FLch3LSeg2.mclear%工作間清零V=54.3,61.8,72.4,88.7,118.6,194.0，P=61.2,49.5,37.6,28.4,19.2,10.1%賦初值并顯示V、P%logP=-alpha*logV+logbeita=-logV,1alpha,log(

10、beita)=HL*sita%注釋P、V之間的關(guān)系fori=1:6;Z(i)=log(P(i);%循環(huán)變量的取值為從1到6，系統(tǒng)的采樣輸出賦值End%循環(huán)結(jié)束ZL=Z%zL賦值HL=-log(V(1),1;-log(V(2),1;-log(V(3),1;-log(V(4),1;-log(V(5),1;-log(V(6),1%Hl賦值%CalculatingParametersc1=HL*HL;c2=inv(c1);c3=HL*ZL;c4=c2*c3%計算被辨識參數(shù)的值%SeparationofParametersalpha=c4(1)%a為c4的第一個元素beita=exp(c4(2)%為以自

11、然數(shù)為底的c4的第二個元素的指數(shù)程序運行結(jié)果：P=61.2000,49.5000,37.6000,28.4000,19.2000,10.1000ZL=4.1141,3.9020,3.6270,3.3464,2.9549,2.3125HL=-3.99451.00001-4.12391.0000-4.28221.0000-4.48531.0000-4.77581.0000-5.26791.0000丿c41.404219.6786alpha=1.4042beita=1.5972e+004仿真結(jié)果表明，用最小二乘一次完成算法可以迅速辨識出系統(tǒng)參數(shù)，即a=1.4042,=1.5972e+004。例3.5

12、考慮圖3.6所示的仿真對象，圖中，v(k)是服從N(0,1)分布的不相關(guān)隨機噪聲。且G(z-1)B(z-1)，N(z-1)D(z-1)，A(z-1)C(z-1)A(zT)1-1.5azT+0.7z2C(z_1)0.5,u(i)=-0.03;%如果M序列的值為1,辨識的輸入信號取“-0.03”elseu(i)=0.03;%如果M序列的值為0,辨識的輸入信號取“0.03”end%小循環(huán)結(jié)束y1=x1;y2=x2;y3=x3;y4=x4;%為下一次的輸入信號做準(zhǔn)備end%大循環(huán)結(jié)束，產(chǎn)生輸入信號ufigure(1);%第一個圖形stem(u),gridon%顯示出輸入信號徑線圖并給圖形加上網(wǎng)格z(2

13、)=0;z(l)=0;%設(shè)z的前兩個初始值為零fork=3:15;%循環(huán)變量從3到15z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%輸出采樣信號end圖3.7最小二乘遞推算法辨識的Malab6.0程序流程圖%RLS遞推最小二乘辨識c0=0.0010.0010.0010.001;%直接給出被辨識參數(shù)的初始值,即一個充分小的實向量p0=10W*eye(4,4);%直接給出初始狀態(tài)P0,即一個充分大的實數(shù)單位矩陣E=0.000000005;%取相對誤差E=0.000000005c=c0,zeros(4,14);%被辨識參數(shù)矩陣的初始值及大小e=zeros

14、(4,15);%相對誤差的初始值及大小fork=3:15;%開始求Kh1=-z(k-1),-z(k-2),u(k-1),u(k-2);x=h1*p0*h1+1;x1=inv(x);%開始求K(k)kl=p0*hl*xl;%求出K的值dl=z(k)-hl*c0;cl=c0+kl*dl;%求被辨識參數(shù)ce1=c1-c0;%求參數(shù)當(dāng)前值與上一次的值的差值e2=e1./c0;%求參數(shù)的相對變化e(:,k)=e2;%把當(dāng)前相對變化的列向量加入誤差矩陣的最后一列c0=c1;%新獲得的參數(shù)作為下一次遞推的舊參數(shù)c(:,k)=c1;%把辨識參數(shù)c列向量加入辨識參數(shù)矩陣的最后一列pl=p0-kl*kl*hl*p

15、0*hl+l;%求出p(k)的值p0=p1;%給下次用ife20.001000.0010-0.4984-1.2328-1.4951-1.4962-1.4991-1.4998-1.49990.001000.00100.0010-0.23500.69130.69410.69900.69980.69990.001000.25091.24971.06651.00171.00201.00020.99990.9998.0.00100-0.24890.75000.56680.50200.50160.50080.50020.50020.03Fig.1InputSignal圖3.8最小二乘遞推算法的參數(shù)辨識仿真

16、-1.5000-1.5000-1.5000-1.4999-1.49990.70000.70000.70000.7000-0.70000.99990.99990.99990.99990.99990.50000.50000.50000.50000.5000e=000-499.42001.47340000-235.991600249.86123.9816-0.146600-249.8612-4.0136-0.24430.00010.0000-0.0000-0.00000.21280.00070.00200.00040.0000-3.94160.00420.00700.00120.0001-0.060

17、70.0003-0.0018-0.0003-0.0001-0.1143-0.0007-0.0016-0.0012-0.00010.00000.0001-0.00000.00000.00000.00000.00010.00000.00000.0000-0.0000-0.00040.0000-0.00000.0000-0.0000表3.2最小二乘遞推算法的辨識結(jié)果參數(shù)aabb真值-1.50.71.00.5估計值-1.49990.70.99990.5000仿真結(jié)果表明，大約遞推到第十步時，參數(shù)辨識的結(jié)果基本達(dá)到穩(wěn)定狀態(tài)，即。=-1.4999,a2=0.7000,=0.9999,b2=0.5000o此

18、時，參數(shù)的相對變化量E0.000000005。從整個辨識過程來看，精度的要求直接影響辨識的速度。雖然最終的精度可以達(dá)到很小，但開始階段的相對誤差會很大，從圖3.8(3)圖形中可看出，參數(shù)的最大相對誤差會達(dá)到三為數(shù)例3.6考慮圖3.10所示的仿真對象，圖中，v(k)是服從N(0,1)分布的不相關(guān)隨機噪聲。u(k)圖3.10增廣最小二乘法辨識實例結(jié)構(gòu)圖D(z-1)，C(z-1)G(z-1),B(z一1)N(z-1)A(z-1)A(z-1)，1-1.5a1z-1+0.7z-2，C(z-1)B(z-1)，1.0z-1+0.5z-2D(z-1)，1-z-1+0.2z-2圖3.11增廣最小二乘遞推算法辨識

19、的Malab6.0程序流程圖模型結(jié)構(gòu)選用如下形式：z(k)+aZ(k一1)+a?(k-2)，b、u(k一1)+buk一2)+v(k)+d、v(k一1)+d?(v(k一2)增廣最小二乘辨識程序流程如圖3.11所示Matlab6.0程序如下。%增廣最小二乘辨識程序，在光盤中的文件名：FLch3ELSeg4.mclearL=60;%4位移位寄存器產(chǎn)生的M序列的周期y1=1;y2=1;y3=1;y4=0;%4個移位寄存器的輸出初始值fori=1:L;x1=xor(y3,y4);%第一個移位寄存器的輸入信號x2=y1;%第二個移位寄存器的輸入信號x3=y2;%第三個移位寄存器的輸入信號x4=y3;%第四

20、個移位寄存器的輸入信號y(i)=y4;%第四個移位寄存器的輸出信號，M序列,幅值0和1，ify(i)0.5,u(i)=-l;%M序列的值為1時，辨識的輸入信號取“-1”elseu(i)=1;%M序列的值為0時,辨識的輸入信號取“1”endy1=x1;y2=x2;y3=x3;y4=x4;%為下一次的輸入信號準(zhǔn)備endfigure(1);%畫第一個圖形subplot(2,1,1);%畫第一個圖形的第一個子圖stem(u),gridon%畫出M序列輸入信號v=randn(1,60);%產(chǎn)生一組60個正態(tài)分布的隨機噪聲subplot(2,1,2);%畫第一個圖形的第二個子圖plot(v),gridon

21、;%畫出隨機噪聲信號u，v%顯示輸入信號和噪聲信號z=zeros(7,60);zs=zeros(7,60);zm=zeros(7,60);zmd=zeros(7,60);%輸出采樣矩陣、不考慮噪聲時系統(tǒng)輸出矩陣、不考慮噪聲時模型輸出矩陣、模型輸出矩陣的大小z(2)=0;z(1)=0;zs(2)=0;zs(1)=0;zm(2)=0;zm(1)=0;zmd(2)=0;zmd(1)=0;%輸出采樣、不考慮噪聲時系統(tǒng)輸出、不考慮噪聲時模型輸出、模型輸出的初值c0=0.0010.0010.0010.0010.0010.0010.001;%直接給出被辨識參數(shù)的初始值,即一個充分小的實向量p0=10W*ey

22、e(7,7);%直接給出初始狀態(tài)P0,即一個充分大的實數(shù)單位矩陣E=0.00000000005;%相對誤差E=0.000000005c=c0,zeros(7,14);%被辨識參數(shù)矩陣的初始值及大小e=zeros(7,15);%相對誤差的初始值及大小fork=3:60;%開始求Kz(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);%系統(tǒng)在M序列輸入下的輸出采樣信號hl=-z(k-l),-z(k-2),u(k-l),u(k-2),v(k),v(k-l),v(k-2);%為求K(k)作準(zhǔn)備x=h1*p0*h1+1;x

23、1=inv(x);k1=p0*h1*x1;%Kd1=z(k)-h1*c0;c1=c0+k1*d1;%辨識參數(shù)czs(k)=-1.5*z(k-1)+0.7*z(k-2)+u(k-1)+0.5*u(k-2);%系統(tǒng)在M序列的輸入下的輸出響應(yīng)zm(k)=-z(k-1),-z(k-2),u(k-1),u(k-2)*c1(1);c1(2);c1(3);c1(4);%模型在M序列的輸入下的輸出響應(yīng)zmd(k)=h1*c1;%模型在M序列的輸入下的輸出響應(yīng)e1=c1-c0;e2=e1./c0;%求參數(shù)誤差的相對變化e(:,k)=e2;c0=c1;%給下一次用c(:,k)=c1;%把遞推出的辨識參數(shù)c的列向量

24、加入辨識參數(shù)矩陣p1=p0-k1*k1*h1*p0*h1+1;%findp(k)p0=p1;%給下次用ife2u=1,-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,1,-1,-1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-11.1418,1.5519,1.3836,-0.7581,0.4427,0.9111,-1.0741,0.2018,0.7629,-1.2882,-0.9530,0.7782,-0.0

25、063,0.5245,1.3643,0.4820,-0.7871,0.7520,-0.1669,-0.8162,2.0941,0.0802,-0.9373,0.6357,1.6820,0.5936,0.7902,0.1053,-0.1586,0.8709,-0.1948,0.0755,-0.5266,-0.6855,-0.2684,-1.1883,0.2486,-0.1025,-0.0410,-2.2476,-0.5108,0.24920.3692,0.1792,-0.0373,-1.6033,0.3394,-0.1311,0.4852,0.5988,-0.0860,0.3253,-0.335

26、1,-0.3224,-0.3824,-0.9534,0.2336,1.2352,-0.5785,-0.5015c=0.001000.0010-0.2789-0.9025-0.9390-0.8620-1.1236-1.5000-1.5000-1.50000.001000.00100.0010-0.0734-0.2661-0.11990.32590.70000.70000.70000.001000.05920.45600.53000.52930.94580.54021.00001.00001.00000.001000.05720.81860.85000.86391.11571.29570.5000

27、0.50000.50000.001000.07960.73420.54170.55160.38470.24591.00001.00001.00000.00100-0.0894-0.5981-0.3420-0.4456-0.4044-0.8443-1.0000-1.0000-1.0000000100-0.0655-0.7795-0.8572-0.7412-0.4506-0.19740.20000.20000.2000V000-279.94262.23550.0404-0.08200.30350.33500.00000.00000000-74.37642.6263-0.5495-3.71891.1

28、479-0.00000.00000058.22126.70070.1622-0.00130.7868-0.42890.85120.00000.000000-58.2212-15.30520.03840.01630.29150.1613-0.6141-0.0000-0.000000-80.5566-10.2283-0.26210.0183-0.3027-0.36073.06630.00000.000000-90.35565.6936-0.42820.3029-0.09251.08800.1844-0.00000.0000L00-66.475410.90590.0996-0.1354-0.3920

29、-0.5619-2.01300.00000.0000z=0,0,-0.4400,-3.9913,-5.7014,-6.9415,-7.8178,-3.9097,1.4543,2.4075,3.5810,6.6598,6.0079,3.5364,2.4376,-0.0964,-2.3472,-2.3178,-4.4100,-6.9915,-6.0233,-5.8181,-3.6094,1.7476,5.5068,6.5756,8.0416,6.3933,2.3550,0.6078,-2.3342,-2.9824,-3.9807,-5.5271,-6.6924,-8.7267,-6.5221,-2

30、.5583,2.1343,2.3062,4.1939,6.4870,6.3125,3.2878,0.8703,-3.0262,-2.7133,-3.2428,-3.7807,-4.8137,-6.6618,-5.5921,-2.9025,1.1385,3.1125,3.7363,6.0362,6.7499,2.6324,0.0478zs=0,0,-0.5000,-0.8401,4.1789,4.2583,6.9212,8.3677,1.8920,-5.4182,-2.0932,-2.1863,-7.9830,-5.8500,-0.5991,-1.6810,2.3509,2.9533,0.333

31、7,3.4925,5.9002,4.6409,6.0108,2.8415,-5.6480,-6.5369,-4.5087,-7.9595,-5.4608,1.4428,0.2369,4.4268,2.3396,2.3834,4.0042,4.6696,8.9054,5.1745,0.7719,-5.4923,-1.4653,-3.1765,-7.2947,-6.4279,-0.0129,0.4961,5.6486,1.4516,1.4649,1.9010,3.0741,7.1232,5.2248,1.9393,-4.2395,-3.3718,-1.9258,6.9389,-7.3995,1.2

32、763zm=0,0,-0.4400,-3.9913,-5.7014,-6.9415,-7.8178,-3.9097,1.4543,2.4075,3.5810,6.6598,6.0079,3.5364,2.4376,-0.0964,-2.3472,-2.3178,-4.4100,-6.9915,-6.0233,-5.8181,-3.6094,1.7476,5.5068,6.5756,8.0416,6.3933,2.3550,0.6078,-2.3342,-2.9824,-3.9807,-5.5271,-6.6924,-8.7267,-6.5221,-2.5583,2.1343,2.3062,4.

33、1939,6.4870,6.3125,3.2878,0.8703,-3.0262,-2.7133,-3.2428,-3.7807,-4.8137,-6.6618,-5.5921,-2.9025,1.1385,3.1125,3.7363,6.0362,6.7499,2.6324,0.0478v=0102030405060Fig.1(2)StochasticNoisesvFig.2ParameterIdentificationWithRecursiveLeastSquaresMethodFig.3IdentificationError0102030405060Fig.4Responses圖3.12

34、增廣最小二乘遞推算法辨識仿真10.50-0.5-1420-2-41.510.50-0.5-1-1.510005000-500100-10100-10100-10200-20v=表3.3增廣最小二乘遞推算法的辨識結(jié)果參數(shù)aia2b1bdd2d3真值-1.50.71.00.51.0-1.00.2估計值-1.50.71.00.51.0-1.00.2仿真結(jié)果表明，遞推到第9步時，參數(shù)辨識的結(jié)果基本達(dá)到穩(wěn)定狀態(tài)，即a1=-1.5000,a2=0.7000,b1=1.0000,b2=0.5000,d1=1.0000,d=-1.0000,d=0.2000。此時，辨識參數(shù)的相對變化量E0.000000005。

35、23與最小二乘遞推算法相比，增廣最小二乘遞推算法雖然考慮了噪聲模型，同樣具有速度快、辨識結(jié)果精確的特點。第4章程序與注釋第4章例：考慮圖4.2表示的線性時不變SISO系統(tǒng)。根據(jù)卷積定理，系統(tǒng)的輸出y（k）與輸入序列,（k_1）,（k_2）,.,（k-n）的關(guān)系可表示成y（k）Xg,（k-i）（4.40）gi1其中，g,g,g組成系統(tǒng)的脈沖響應(yīng)。系統(tǒng)在偽隨機碼輸入作用下12N的輸出響應(yīng)如表4.1所示。u(k)G(s)y(k)-圖4.2單輸入單輸出系統(tǒng)表4.1輸入輸出數(shù)據(jù)k012345678910111213141516171819.u(k)-1-1-1-1111-111-1-11-11-1-1-

36、1-11.y(k)0-2-6-7-7-3573-153-5-31-11-5-7-7.利用這些數(shù)據(jù)辨識系統(tǒng)的脈沖響應(yīng)，當(dāng)N取3時，有y(k)二g卩(ki)ii=1令(4.41)h(k)=卩(k-1),卩(k-2),卩(k-3)g=g，g，g匕123若系統(tǒng)的真實脈沖響應(yīng)記作g,則有(4.42)y(k)=h”k)g0根據(jù)式(4.16),脈沖響應(yīng)估計值直伙)的遞推算式為g(k+1)=g(k)+R(k)h(k)y(k)，ht(k)g(k)(4.43)(4.44)其中，脈沖響應(yīng)估計值的初始值取g=0，權(quán)矩陣R(k)取如下形式R*(k)=diaga(k)，力(k)，力(k)N123A(k)h2(k)iii=

37、1A(k)=1,A(k)=12117A3(k)=了(4.45a)(4.45b)確定性問題梯度校正遞推算法辨識的流程如圖4.3所示。下面給出Malab6.0程序。%梯度校正參數(shù)辨識程序，在光盤中的文件名：FLch4GAeg1.mclearU=-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,-1,-1,1;y=0,-2,-6,-7,-7,-3,5,7,3,-1,5,3,-5,-3,1,-1,1,-5,-7,-7;%畫出u和y圖形figure(1),subplot(2,1,1),stem(u),subplot(2,1,2),stem(y),holdonk=1:

38、20；plot(k,y)%給出初始值h1=-1,0,0;h2=-1,-1,0;g=0,0,0;I=1,0,0;0,1/2,0;0,0,1/4;h=h1,h2,zeros(3,16);%計算樣本數(shù)據(jù)h(k)fork=3:18h(:,k)=u(k),u(k-1),u(k-2);end%計算權(quán)矩陣R(k)和g的估計值fork=1:18a=h(l,k)A2+(h(2,k)A2)/2+(h(3,kF2)/4;al=l/a;R=al*I;%按照式(4.45a)和(4.45b)計算權(quán)矩陣g(:,k+1)=g(:,k)+R*h(:,k)*(y(k+1)-h(:,k)*g(:,k);%按照式(4.44)計算脈沖

39、響應(yīng)的估計值end%繪圖g1=g(1,:);g2=g(2,:);g3=g(3,:);figure(2);k=1:19;subplot(121);plot(k,g1,r,k,g2,g,k,g3,b),gridon%計算模型輸出ym及系統(tǒng)輸出與模型輸出之間的誤差Eyfork=1:18ym(k)=h(:,k)*g(:,k);Ey(k)=y(k+1)-ym(k);endk=1:18;subplot(122);plot(k,Ey),gridong,ym,Ey%顯示脈沖響應(yīng)估計值、模型輸出及系統(tǒng)輸出與模型輸出之0.99360.99720.99570.99630.99740.99740.99740.9972

40、0.9980間的誤差figure(3);x=0:1:3;y=0,g(1,18),g(2,18),g(3,18);xi=linspace(0,3);yi=interp1(x,y,xi,cubic);plot(x,y,o,xi,yi,m),gridon%畫出脈沖響應(yīng)估計值及其三次插值曲線圖4.3確定性梯度校正辨識的Malab6.0程序流程圖程序執(zhí)行結(jié)果：g=02.00004.66675.23815.23811.53742.14382.23931.96581.9559001.33331.61901.61903.46943.77263.82043.95713.9621L0000.14290.14291

41、.06800.91640.94031.00871.00622.00631.99181.99801.99531.99951.99951.99952.00012.00323.98733.99463.99773.99903.99693.99693.99693.99723.9988ym=0,-2.0000,-6.0000,-7.0000,3.4762,3.9388,6.8328,2.5213,-0.9826,4.9118,2.9746,-4.9891,-2.9953,1.0073,-1.0000,1.0000,-4.9990,-6.9945Ey=-2.0000,-4.0000,-1.0000,0.00

42、00,-6.4762,1.0612,0.1672,0.4787,-0.0174,0.0882,0.0254,0.0109,-0.0047,-0.0073,-0.0000,0,-0.0010,-0.005505101520Fig.1OutputErrorBetweenSystemandim!V.OdeljI=j210-1-2-3-4-5-62010Fig.2(2)-70110.99360.99720.99570.99630.99740.99740.99740.99720.9980-1-1005101520Fig.1(2)O20esponse-1010Fig.2(1)DiscreteImpulse

43、RutputSignal20-7010Fig.2(2)OutputErrorFig.3Impulsereponsecurve圖4.5確定性問題梯度校正參數(shù)辨識仿真（被辨識參數(shù)的個數(shù)為5）圖4.4表明，遞推到第10步時，被辨識參數(shù)基本上達(dá)到了穩(wěn)定狀態(tài)，即g12,g24，g3l；此時系統(tǒng)的輸出與模型的輸出誤差也基本達(dá)到穩(wěn)定狀態(tài)，即Ey0。由于被辨識參數(shù)的個數(shù)較少，遞推校正算法的收斂性比較好，也就是說，輸入輸出的觀測數(shù)據(jù)量已足夠。但從圖4.4Figure3中可看出，由于被辨識參數(shù)的個數(shù)較少，它還不足以充分顯示全部的系統(tǒng)脈沖響應(yīng)。0.99360.99720.99570.99630.99740.9974

44、0.99740.99720.9980為了充分地顯示出系統(tǒng)的脈沖響應(yīng)，可以把被辨識參數(shù)的個數(shù)增加到5個。圖4.5給出了被辨識參數(shù)的個數(shù)為5時的辨識結(jié)果。圖4.5基本顯示出了系統(tǒng)的全部脈沖響應(yīng)過程。但在圖4.5的辨識中，還是利用上面給出的20對輸入輸出數(shù)據(jù)。通過圖4.4與圖4.5的比較可以看出，在觀測數(shù)據(jù)量一定的情況下，隨著被辨識參數(shù)的個數(shù)的增加，梯度校正辨識算法的收斂性變差，這是由于觀測數(shù)據(jù)不足，遞推步數(shù)有限造成的。所以，隨著被辨識參數(shù)的個數(shù)的增加，若要保持梯度校正辨識算法的收斂效果，需要同步增加觀測數(shù)據(jù)量。第6章程序與注釋例6.1考慮圖6.1所示的仿真對象，圖中，v(k)是均值為零、方差為2的

45、服從正態(tài)分布的不相關(guān)隨機噪聲。v+圖6.1貝葉斯辨識結(jié)構(gòu)圖G(z-1),B(z一1)N(z-1),D(z1)A(z-1)C(z-1)A(z-1)，l-1.5az-1+0.7z-2，C(z-1)B(z-1)，1.0z-1+0.5z-2D(z-1)，1-z-1+0.2z-2模型結(jié)構(gòu)選用如下形式。z(k)+az(k-1)+a?(k一2)，b、u(k-1)+一2)+v(k)+d、v(k-1)+d?(v(k一2)Bayes辨識算法程序流程如圖6.2所示，Matlab6.0程序如下。%Bayes辨識程序，在光盤中的文件名：FLch6BAegl.m。clearL=60;%四位移位寄存器產(chǎn)生的M序列的長度y1

46、=1;y2=1;y3=1;y4=0;%四個移位寄存器的輸出初始值fori=1:L;x1=xor(y3,y4);%第一個移位寄存器的輸入信號x2=y1;%第二個移位寄存器的輸入信號x3=y2;%第三個移位寄存器的輸入信號x4=y3;%第四個移位寄存器的輸入信號y(i)=y4;%第四個移位寄存器的輸出信號ify(i)0.5,u(i)=-l;%M序列的值為,1時，辨識的輸入信號取“-1”elseu(i)=1;%M序列的值為0時，辨識的輸入信號取“1”endy1=x1;y2=x2;y3=x3;y4=x4;%為下一次的輸入信號作準(zhǔn)備endfigure(1);subplot(2,1,1);stem(u),

47、gridon%畫第一個圖形的第一個子圖M序列輸入信號v=randn(1,60);subplot(2,1,2);%產(chǎn)生一組長度為60的正態(tài)分布隨機噪plot(v),gridon;%畫第一個圖形窗口的第二個子圖隨機噪聲信號R=corrcoef(u,v);%計算輸入信號與隨機噪聲信號的相關(guān)系數(shù)r=R(1,2)%取出并顯示互相關(guān)系數(shù)rv=std(v)*std(v)%計算并顯示隨機噪聲的方差u，v%顯示輸入信號和噪聲信號z(2)=0;z(1)=0;zs(2)=0;zs(1)=0;%輸出采樣、系統(tǒng)實際輸出、模型輸出賦初值zm(2)=0;zm(1)=0;%模型輸出賦初值%Bayes辨識c0=0.0010.0

48、010.0010.0010.0010.0010.001;%直接給出被辨識參數(shù)的初始值,即一個充分小的實向量p0=10W*eye(7,7);%直接給出初始狀態(tài)P0,即一個充分大的實數(shù)單位矩陣E=0.00000000005;%取相對誤差E=0.000000005為停機準(zhǔn)則圖6.2貝葉斯辨識的Malab6.0程序流程圖c=c0,zeros(7,14);%被辨識參數(shù)矩陣的初始值及大小e=zeros(7,15);%相對誤差的初始值及大小fork=3:20;%開始求Kz(k)=-1.5*z(k-1)+0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);%

49、系統(tǒng)在M序列輸入下的輸出采樣信號hl=-z(k-l),-z(k-2),u(k-l),u(k-2),v(k),v(k-l),v(k-2);%為求K(k)作準(zhǔn)備x=h1*p0*h1+rv;x1=inv(x);k1=p0*h1*x1;%Kdl=Z(k)-hl*C0；%開始求被辨識參數(shù)cc1=c0+k1*d1;%辨識參數(shù)czs(k)=-1.5*z(k-l)+0.7*z(k-2)+u(k-l)+0.5*u(k-2);%系統(tǒng)在M序列輸入下的輸出響應(yīng)zm(k)=-z(k-1),-z(k-2),u(k-1),u(k-2)*c1(1);c1(2);c1(3);c1(4);%模型的輸出響應(yīng)e1=c1-c0;e2=

50、e1./c0;%求參數(shù)的相對變化e(:,k)=e2;c0=c1;%給下一次用c(:,k)=c1;%把辨識參數(shù)c列向量加入辨識參數(shù)矩陣p1=p0-k1*k1*h1*p0*h1+1;%findp(k)p0=p1;%給下次用ife2=Ebreak;%若收斂情況滿足要求，終止計算endendc,e,z,zs,zm%顯示被辨識參數(shù)、收斂情況、輸出采樣值、實際輸出值、模型輸出值%分離賦值a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);%分離出a1、a2、b1、b2d1=c(5,:);d2=c(6,:);d3=c(7,:);%分離出d1、d2、d3ea1=e(1,:);ea2

51、=e(2,:);eb1=e(3,:);eb2=e(4,:);%分離出a1、a2、b1、b2的收斂情況ed1=e(5,:);ed2=e(6,:);ed3=e(7,:);%分離出d1、d2、d3的收斂情況figure(2);%畫第二個圖形i=1:20;plot(i,a1,r,i,a2,r:,i,b1,b,i,b2,b:,i,d1,g,i,d2,g:,i,d3,g+)%畫出被辨識參數(shù)的各次估計值title(ParameterIdentificationwithRecursiveLeastSquaresMethod)%標(biāo)題figure(3);%畫出第三個圖形i=1:20;plot(i,ea1,r,i,

52、ea2,r:,i,eb1,b,i,eb2,b:,i,ed1,g,i,ed2,g:,i,ed2,r+)%畫出各個參數(shù)收斂情況title(IdentificationError)%標(biāo)題%Response%響應(yīng)figure(4);subplot(4,1,1);%畫出第四個圖形中的第一個子圖i=1:20;plot(i,zs(i),r)%畫出被辨識系統(tǒng)在沒有噪聲情況下的實際輸出響應(yīng)subplot(4,1,2);%畫出第四個圖形中的第二個子圖i=1:20;plot(i,z(i),g)%畫出被辨識系統(tǒng)輸出的采樣值subplot(4,1,3);%畫出第四個圖形中的第三個子圖i=1:20;plot(i,zm(i

53、),b)%畫出模型含有噪聲的輸出響應(yīng)subplot(4,1,4);%畫出第四個圖形中的第四個子圖i=1:20;plot(i,zmy(i),b)%畫出模型去除噪聲后的輸出響應(yīng)程序執(zhí)行結(jié)果如下：u=1,-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,1,-1,-1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1v=-0.4326,-1.6656,0.1253,0.2877,-1.1465,1.1909,1.189

54、2,-0.0376,0.3273,0.1746,-0.1867,0.7258,-0.5883,2.1832,TOC o 1-5 h z1-55555575575757-0.1364,0.1139,1.0668,0.0593,-0.0956,-0.8323,0.2944,-1.3362,0.7143,1.6236,-0.6918,0.8580,1.2540,-1.5937,55555555775557-1.4410,0.5711,-0.3999,0.6900,0.8156,0.7119,1.2902,0.6686,1.1908,-1.2025,-0.0198,-0.1567,-1.6041,0.2573,-1.0565,1.4151,-0.8051,0.5287,0.2193,-0.9219,-2.1707,-0.0592,-1.0106,0.6145,0.5077,1.6924,0.5913,-0.6436,555557555557570.3803,-1.0091,-0.0195,-0.0482rv=0.8991r=0.0591c=0.001000.0010-0.0004-0.1558-0.9692-1.2967-1.3339-1.5000-1.5000-1.5000-1.50000.001000.00