模型階次辨識.

1、盛曉婷 模型階次辨識 模型階次辨識 盛曉婷0908110618、二、三、四、五、 利用行列式比估計模型的階次 2利用殘差的方差估計模型的階次 3利用Akaike準則估計模型的階次4利用最終預報誤差準則估計模型的階次5根據(jù)Hankel矩陣的秩估計模型的階次 7利用行列式比估計模型的階次 9附錄1附錄2利用殘差的方差估計模型的階次 11附錄3利用Akaike準則估計模型的階次13附錄4利用最終預報誤差準則估計模型的階次15附錄5利用Hankel矩陣的秩估計模型的階次 17 1盛曉婷模型階次辨識一、利用行列式比估計模型的階次 例13.2考慮如下方針對 象 z ( k ) - 1.630 z ( k

2、- 1) + 0.638 z ( k - 2) = 0.01u ( k ) + 0.0004u (k - 1) - 0.0024u (k - 1) + v (k )其中，v (k )是均值為零，方差為 cv2( 0 (T v2 )的不相關A* ? L n +1 n +1 ? A隨機噪聲；輸入信號u (k )采用M序列。數(shù)據(jù)長度L=1000.定義行列式比DR ( n =* a det H * (n) det H * (n + 1) a a 其中，1 *t *?* a a a ? h ( n) = L H n H n ? ? a ? H * (n + 1) = 1 H a*T H a* ?L n

3、+1 n +1 ? a ?az (n - 1) . z (1)|? z ( n)a ? a * z ( n)z (2)| ? z (n + 1)屮=n ? |? a a?z (n + L - 1) z (n+ L -2) . z ( L)| ? a u (n) u (n - 1) . u (1) ? ? a a ? u (n + 1)u(n) . u (2) ? ? a a u (n + L - 1) u (n+ L - 2) . u ( L)? ? a a * a a 當n從 1 開始逐一增加，如若DR ( n)較DR (n - 1)有顯著增加，則這時的n認為已比較接* a近過程模型的真 實

4、階次，即應取 n0 = n。a dr* (1) = 1.6837 DR* (2) = 3.1289 DR* (3) = 3.3043 當 n = 2時，行列式比顯著增加，故模型的階次可定為2階。Matlab程序見附錄1. a2盛曉婷 模型階次辨識 利用行列式比估計模型階次 3.5 3 2.5行列式比2 1.5 1 0.5 0 1 1.2 1.4 1.6 1.8 2階次2.2 2.4 2.6 2.8 3二、利用殘差的方差估計模型的階次 例 13.4 考慮如下的仿真對象 z ( k ) - 1.5 z ( k - 1) + 0.7 z ( k - 2) = u ( k - 1) + 0.5u (

5、k -2) + v (k )其中，輸入u ( k )為M序列，噪聲v ( k )為均值為0,方差為1的白 噪聲。最小二乘意義下，系數(shù)的辨識值為 ？A ? ? a1 ? ?A ? A a2 A 0 ? :a ?= (H at H a ) -1 H at Zn0 ? ? n n n ? b1 ? ? a ? ?b2 ? ? ?hT ？ ？- Z (0)-Z (-1) u (0) u (1 - n)? ? Z (1) ? ? T ? ? ? ? ?u (1) u (2- n)? ?= ?-Z (1)-Z (0) ? Z ？ ， H = ? h 其中，Z L = L ? . ? ? . ?.? ? ?

6、 ? ? hT (1002) ? ?- Z ( L - 1) - Z ( L - 2) u ( L - 1)u ( L - n)?Z(L) ? ? a其中，n表示模型的階次估計值，n0為過程的真實階次。3盛曉婷模型階次辨識模型的輸出殘差為a a z n = zn0 - H a 0 n殘差的方差 為V1 ( n) = a a 1 at a zn * zn L a a a n從1開始逐一增加，若L充分大，則 V1 (n)將隨n遞減，減到某一固定值時，對應的n可認為是過程模型的階次。由 下圖可見，過程模型的階次是 2.利用殘差的方差估計模型的階次 2.8 2.6 2.4 2.2 殘差方差2 1.8

7、1.6 1.4 1.2 1 1 1.5 2 2.5階次3 3.5 4 Matlab程序見附錄2.三、利用Akaike準則估計模型的階次 考慮如下的仿真對象z ( k ) - 1.5 z ( k - 1) + 0.7 z ( k -2) = u ( k - 1) + 0.5u ( k - 2) + v (k )其中，輸入 u ( k )為 M 序列，噪聲 v ( k )為 均值為0,方差為1的白噪聲。最小二乘意義下，系數(shù)的辨識值為 4盛曉婷 模型階次辨識 T A A A A ? A 9 ML = a1 a 2 b1 b 2? = ( H AT H A ) -1 H ATZ n0 ? nn n ?

8、2 a a 1 ? a 小=(zn - H n9 ML )T ( zn - H n ?9 ML )? L ? hT ？ ？ - Z (0) - Z (-1) u (0) u (1 - n) ? ? Z (1) ? ? T ? ? ? ? ? u (1) u (2 - n) ? ? = ? - Z (1) - Z (0) ? Z ？ ， H = ? h 其中，Z L = L ? . ? . ? . ? ? ? ? ? ? ? ? hT (1002) ? ? - Z ( L - 1) - Z ( L - 2) u ( L - 1)u ( L - n) ? ? ? Z ( L) ? ? ? a其中，

9、n表示模型的階次估計值，n0為過程的真 實階次。AIC 準則為 AIC (n) = L log( c v ) + 4 n a a a 2 a a分別計算 AIC ( n) 值，找到使其最小的n作為模型的階次，對應的n可n從1開始逐一增加，認 為是過程模型的階次。由下圖可見，過程模型的階次是2. AIC判斷模型階次1000 800 600 AIC 400 200 0 -200 1 1.5 2 2.5階次 3 3.5 4 Matlab程序見附錄 3.四、 利用最終預報誤差準則估計模型的階次考慮如下的仿真對象5盛曉婷模型階次辨識 z ( k ) - 1.5 z ( k - 1) + 0.7 z (

10、k - 2) = u ( k - 1) + 0.5u ( k - 2)+ v (k )其中，輸入u ( k )為M序列，噪聲v ( k )為均值為0,方差為1的白噪 聲。最小二乘意義下，系數(shù)的辨識值為 ？A ? ? a1 ? ?A ? A a2 A 9?去?=(H at H a ) -1 H atZ n0 ? n n n ? b1 ? ? a ? ? b2 ? ? ?hT ？ ？ - Z (0)-Z (-1) u (0) u (1-n) ? ?Z (1) ? ? T ? ? ? ? ? u (1) u (2 -n)? ?= ?- Z (1)-Z (0) ? Z ？，H = ?h 其中，Z L

11、= L ? . ? ? . ?.? ? ? ? ? hT (1002)? - Z (L - 1) - Z ( L - 2) u ( L - 1) u ( L - n)?Z( L) ? ? a其中，n表示模型的階次估計值，n0為過程的真實階次。模型的輸出殘差為 a a z n = zn0 - H a 9 n殘差的方差為 c = a 2 1 - at a zn * zn L L + (n a + nb ) L - (n a + nb ) a a a a 最終預報誤差準貝q為 fPE (n a , nb ) = a a a a 其 a 2 中，n a、b為模型階次估計值，此處設 n b已知，且nb

12、= 2。n a a a a a a n從1 開始逐一增加，當FPE (n a , nb )最小，對應的n可認為是過程模型的階次。由下 圖可見，過程模型的階次是2. 6盛曉婷 模型階次辨識 利用最終預報誤差準則估計模型的階次3 2.8 2.6 2.4 2.2 FPE2 1.8 1.6 1.4 1.2 1 1 1.5 2 2.5階次 3 3.5 4 Matlab程序見附錄 4.五、根據(jù) Hankel 矩 陣的秩估計模型的階次 例13.1已知過程的脈沖響應序列如下所示： g=1 0.8 0.650.54 0.46 0.39 0.35 0.31 0.280.260.240.230.220.21 0.2

13、0 0.19 0.19 0.18 0.18 0.180.17 0.17 0.17 0.16 0.16 0.150.150.150.150.140.14 0.14 0.13 0.13 0.13 0.13 0.120.12 0.12 0.12 0.12 0.11 0.110.110.110.110.110.11; 1 4國 det H (2, k ) 46 k =1D2 = = -81.8700 1 44 刀 det H (3, k ) 44 k =1 1 44刀 det H (3, k ) 44 k =1 = 32.6610D3 = 1 42 det H (4, k )刀 42 k =1 1 4

14、2 刀 det H (4,-862426031=1 40刀 det H (5, k ) 40 k =1 7 盛曉婷模型階次辨識g (k + 1) ? g (k ) ? g (k + 1) g (k + 2)其中H (l , k ) = ? ? g (k + l - 1) g (k + l ) . g (k + l - 1)? . g (k + l ) ? ? ? . g (k +21 - 2) ?顯然，D3 = max，所以過程模型的介詞可定為 3階。Matlab程序見附 錄5. 8盛曉婷 模型階次辨識 附錄1利用行列式比估計模型的階次 close all clearcic %= % 利用行列

15、式比估計 模型的階次%見過程辨識，方崇智，清華大學出版社，P337,例13.2 % %=產(chǎn)生 M 序列作為輸入= x=0 1 0 1 1 0 1 1 1; %initial value n=1003; %n 為脈沖數(shù)目 M=; % 存放 M 序列 for i=1:n temp=xor(x,x(9); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp; end % 產(chǎn)生高 斯白噪聲 v=randn(1,1003); z=; z(1)=-1; z(2)=0; L=1000; for i=3:L+3 z(i)=1.630*z(i- 1)-0.638

16、*z(i-2)+0.0004*M(i-1)-0.0024*M(i-2)+v(i); end % n=1 for i=1:L H1(i,1)=z(i); H1(i,2)=M(i); end A=H1*H1/L; 9 盛曉婷 模型階次辨識 % n=2 for i=1:L H2(i,1)=z(i+1); H2(i,2)=z(i); H2(i,3)=M(i+1); H2(i,4)=M(i); end B=H2*H2/L; %n=3 for i=1:L H3(i,1)=z(i+2); H3(i,2)=z(i+1); H3(i,3)=z(i); H3(i,4)=M(i+2); H3(i,5)=M(i+1)

17、; H3(i,6)=M(i); end C=H3*H3/L; %n=4 for i=1:L H4(i,1)=z(i+3); H4(i,2)=z(i+2); H4(i,3)=z(i+1); H4(i,4)=z(i); H4(i,5)=M(i+3); H4(i,6)=M(i+2); H4(i,7)=M(i+1); H4(i,8)=M(i); end D=H4*H4/L;DR(1)=det(A)/det(B); DR(2)=det(B)/det(C); DR(3)=det(C)/det(D); i=1:3; stem(i,DR) title(利用行列式比估計模型階次)xlabel(階次)ylabel

18、(行列式比)10盛曉婷 模型階次辨識 附錄2利用殘差的方差估計模型的階次%利用殘差的方差估計模型的階次方崇智，過程辨識，清華大學出版社；P345,例13.4 %= close all clearclc %=產(chǎn)生 M 序列作為輸入= x=0 1 0 1 1 0 1 11; %initial value n=1003; %n 為脈沖數(shù)目 M=; % 存放 M 序列 for i=1:n temp=xor(x(4),x(9); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp; end % 產(chǎn)生高 斯白噪聲 v=randn(1,1004); z=;

19、z(1)=-1; z(2)=0; L=1000; for i=3:L+4 z(i)=1.5*z(i- 1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i); end % 模型階次 n=1 for i=1:L H1(i,1)=z(i); H1(i,2)=M(i); end estimate=i nv(H1*H1)*H1*z(2:1001); e=z(2:1001)-H1*estimate; V1=e*e/L; 11 盛曉婷 模型階次辨識 % 模型階次 n=2 for i=1:L H2(i,1)=z(i+1); H2(i,2)=z(i); H2(i,3)=M(i+1); H2(

20、i,4)=M(i); end estimate=i nv(H2*H2)*H2*z(3:1002); e=z(3:1002)-H2*estimate; V2=e*e/L; % 模型階次 n=3 for i=1:L H3(i,1)=z(i+2); H3(i,2)=z(i+1); H3(i,3)=z(i); H3(i,4)=M(i+2); H3(i,5)=M(i+1); H3(i,6)=M(i); end estimatednv(H3*H3)*H3*z(4:1003); e=z(4:1003)-H3*estimate; V3=e*e/L; % 模型 階次 n=4 for i=1:L H4(i,1)=

21、z(i+3); H4(i,2)=z(i+2); H4(i,3)=z(i+1); H4(i,4)=z(i); H4(i,5)=M(i+3); H4(i,6)=M(i+2); H4(i,7)=M(i+1); H4(i,8)=M(i); end estimate=i nv(H4*H4)*H4*z(5:1004); e=z(5:1004)-H4*estimate; V4=e*e/L; plot(1:4,V1 V2 V3 V4) title(利用殘差的方差估計模型的階次)xlabel(階次) ylabel(殘差方差)12盛曉婷模型階次辨識 附錄3利用Akaike準則估計模型的階次 利用Akaike準則

22、估計模型的階次方崇智，過程辨識，清華大學出版社；P348,例13.6 %= close all clearclc %=產(chǎn)生 M 序列作為輸入= x=0 1 0 1 1 0 1 11; %initial value n=1003; %n 為脈沖數(shù)目 M=; % 存放 M 序列 for i=1:n temp=xor(x(4),x(9); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp; end % 產(chǎn)生高 斯白噪聲 v=randn(1,1004); z=; z(1)=-1; z(2)=0; L=1000; for i=3:L+4 z(i)=1.

23、5*z(i- 1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i); end % 模型階次 n=1 for i=1:L H1(i,1)=z(i); H1(i,2)=M(i); end estimate1=i nv(H1*H1)*H1*z(2:1001); D仁(z(2:1001)- H1*estimate1)*(z(2:1001)-H1*estimate1)/L; AIC1= L*log(D1)+4*1; 13 盛曉婷 模型階次辨識 % 模型階次 n=2 for i=1:L H2(i,1)=z(i+1); H2(i,2)=z(i); H2(i,3)=M(i+1); H2(i

24、,4)=M(i); end estimate2=i nv(H2*H2)*H2*z(3:1002); D2=(z(3:1002)-H2*estimate2)*(z(3:1002)-H2*estimate2)/L;AIC2=L*log(D2)+4*2; % 模型階次 n=3 for i=1:L H3(i,1)=z(i+2); H3(i,2)=z(i+1); H3(i,3)=z(i); H3(i,4)=M(i+2); H3(i,5)=M(i+1); H3(i,6)=M(i); endestimate3=i nv(H3*H3)*H3*z(4:1003); D3=(z(4:1003)-H3*estima

25、te3)*(z(4:1003)- H3*estimate3)/L; AIC3=L*log(D3)+4*3; % 模型階次 n=4 for i=1:L H4(i,1)=z(i+3); H4(i,2)=z(i+2); H4(i,3)=z(i+1); H4(i,4)=z(i); H4(i,5)=M(i+3); H4(i,6)=M(i+2); H4(i,7)=M(i+1); H4(i,8)=M(i); end estimate4=i nv(H4*H4)*H4*z(5:1004); D4=(z(5:1004)-H4*estimate4)*(z(5:1004)-H4*estimate4)/L; AIC4=

26、L*log(D4)+4*4; plot(1:4,AIC1 AIC2 AIC3 AIC4) title(AIC判斷模型階次)xlabel(階次)ylabel(AIC) 14盛曉婷 模型階次辨識 附錄4利用最終預報誤差準則估計模型的階次%利用最終預報誤差準則估計模型的階次方崇智，過程辨識，清華大學出版社；P355 %= close all clearclc %=產(chǎn)生 M 序列作為輸入= x=0 1 0 1 1 0 1 11; %initial value n=1003; %n 為脈沖數(shù)目 M=; % 存放 M 序列 for i=1:n temp=xor(x(4),x(9); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp; end % 產(chǎn)生高 斯白噪聲 v=randn(1,1004); z=; z(1)=-1