1、第 8 章系統(tǒng)辨識 MATLAB 仿真實例MATLAB 是一套高性能數(shù)字計算和可視化軟件, 它集成概念設(shè)計 , 算法開發(fā) , 建模仿真 , 實時實現(xiàn)于一體 ,構(gòu)成了一個使用方便、界面友好的用戶環(huán)境,其強大的擴展功能為各領(lǐng)域的應(yīng)用提供了基礎(chǔ)。對于一個簡單的系統(tǒng) ,可以通過分析其過程的運動規(guī)律 ,應(yīng)用一些已知的定理和原理 ,建立數(shù)學(xué)模型 , 即所謂的 “白箱建模 ”。 但對于比較復(fù)雜的生產(chǎn)過程 ,該建模方法有很大的局限性。因為過程的輸入輸出信號一般總是可以測量的,而且過程的動態(tài)特性必然表現(xiàn)在這些輸入輸出數(shù)據(jù)中,那么就可以利用輸入輸出數(shù)據(jù)所提供的信息來建立過程的數(shù)學(xué)模型。這種建模方法就稱為系統(tǒng)辨識。

2、把辨識建模稱作 “黑箱建模 ”。本章列舉了在系統(tǒng)辨識中常用的 基本算法的 MA TLAB 程序。8.1 白噪聲的產(chǎn)生 如果在計算機上比較經(jīng)濟的產(chǎn)生統(tǒng)計上理想的各種不同分布的白噪聲序列，則對系統(tǒng)辨識仿真研究提供 了極大方便。為了簡單起見，常把各種不同分布的白噪聲序列稱為隨機數(shù)。從理論上講，只要有了一種具有 連續(xù)分布的隨機數(shù)，就可以通過函數(shù)變換的方法產(chǎn)生其他任意分布的隨機數(shù)。顯然，在具有連續(xù)分布的隨機 數(shù)中， 0,1）均勻分布的隨機數(shù)是最簡單、最基本的一種，有了0,1 ）均勻分布的隨機數(shù)，便可以產(chǎn)生其他任意分布的隨機數(shù)和白噪聲。本文首先產(chǎn)生 0,1 ）均勻分布的隨機序列，然后對程序稍加改動，將產(chǎn)生

3、0,1 ）均勻分布的隨機數(shù)統(tǒng)統(tǒng)減去 0.5，相當于將原隨機序列圖的橫坐標向上平移0.5，原隨機序列變成了 -0.5,0.5 ）的白噪聲，然后乘以存儲器 f 的系數(shù)，這里取法 f=2 ，則可得到產(chǎn)生 -1,1）均勻分布的白噪聲。8.1.1 乘同余法的介紹用如下遞推同余式產(chǎn)生正整數(shù)序列 xi xi Axi 1（mod M ）,i 1,2,3 2 是整數(shù)；若 k=3，則 mod8）或 k=8，則 mod256），且 A 應(yīng)取適中的值，如 3AM可以證明序列 i 是偽隨機數(shù)序列；同時還可以證明偽隨機序列 i 的循環(huán)周期達到最大值 2k 2。 將式8.1）和 8.2）合并為i A i 1 8.3） 式中

4、，初值為 0 x0c 。0 M c8.1.2 舉例利用乘同余法，選 R=2,A=6,k=8, M=2k 256 ,遞推 100 次，采用 MA TLAB 的仿真語言編程，產(chǎn)生 編程如下：A=6。 x0=1。 M=255 。 f=2。 N=100； %初始化；x0=1 。 M=255 。for k=1: N % 乘同余法遞推 100 次；x2=A*x0 。%分別用 x2和x0表示 xi+1和xi-1；x1=mod （x2,M 。 %取x2存儲器的數(shù)除以 M 的余數(shù)放 x1=(v1-0.5 *f 。 %將 v1中的數(shù) 表示行不變、列隨遞推循環(huán)次數(shù)變化；x0=x1 。% xi-1= xi；v0=v1

5、 。end%遞推 100 次結(jié)束；v2=v % 該語句后無 ；，實現(xiàn)矩陣存儲器 v中隨機數(shù)放在 v2中，且可直接顯示在 MA TLAB 的 window 中；k1=k。%grapher% 以下是繪圖程序；k=1:k1 。plot(k,v,k,v,r 。xlabel(k, ylabel(v 。tktle( (-1,+1 均勻分布的白噪聲 程序運行結(jié)果如圖 8.1所示。圖8.1采用MATLAB 產(chǎn)生的 (-1,+1均勻分布的白噪聲序列產(chǎn)生的 (-1， 1均勻分布的白噪聲序列 在程序運行結(jié)束后，產(chǎn)生的 (-1， 1均勻分布的白噪聲序列，直接從MATLAB 的 window 界面中 copy 出來如下

6、 v2 中每行存 6 個隨機數(shù))：v2 =-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.101

7、6-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500

8、-0.4844/ 220.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359-0.9531 -0.7188 0.6875 -0.83598.2 M 序列的產(chǎn)生 在進行系統(tǒng)辨識時，選用白噪聲作為辨識輸入信號可以保證獲得較好的辨識效果，但在項目上難以實 現(xiàn)。 M 序列是一種很好的辨識輸入信號，它具有近似白噪聲的性質(zhì)，不僅可以保證有較好的辨識效果，而且 項目上又易于實現(xiàn)。8.2.1 M 序列的產(chǎn)生方式M 序列是一種離散二位式隨機序列，所謂 “二位式 ”是指每個隨機變量只有兩種狀態(tài)。離散二位式隨機 序列是按照確定的方式產(chǎn)生的，實際上是一種確定序列?？捎枚嗉壘€性反饋移位寄

9、存器產(chǎn)生 M 序列。每級 移位寄存器由雙穩(wěn)態(tài)觸發(fā)器和門電路組成，稱為1位，分別以 0和 1來表示 2中狀態(tài)。當移位脈沖來到時，每位的內(nèi)容 0 或 1)移到下一位，最后 1 位編程如下：X1=1 。X2=0 。 X3=1 。 X4=0 。 %移位寄存器輸入 Xi 初 T態(tài) 。 % 異或運算if Y4=0U(i=-1 。elseU(i=Y4 。endendM=U%繪圖i1=i k=1:1:i1 。 plot(k,U,k,U,rx xlabel(k ylabel(M 序列 title(移位寄存器產(chǎn)生的 M 序列 程序運行結(jié)果如圖 8.2 所示。/ 22圖 8.2 軟件實現(xiàn)的移位寄存器產(chǎn)生的 M 序列

10、圖四級移位寄存器產(chǎn)生的M 序列M =Columns 1 through 10-1 1 -1 1 1 11-1-1Columns 11 through 201 -1 -1 1 1 -11-11Columns 21 through 301 1 -1 -1 -1 1-1-11Columns 31 through 40-1 1 -1 1 1 11-1-1Columns 41 through 501 -1 -1 1 1 -11-11Columns 51 through 601 1 -1 -1 -1 1-1-11-111-111i1 =608.3 最小二乘一次完成算法的產(chǎn)生 考慮仿真對象(8.4其中， v

11、( k )是服從正態(tài)分布的白噪z(k) 1.5z(k 1) 0.7z(k 2) u(k 1) 0.5u(k 2) v(k)聲 N(0,1 。輸入信號采用 4階M 序列，幅度為 1。選擇如下形式的辨識模型M 序列，則待辨識參數(shù) ?LS為?LS = ( H LH L ) 1H LzL。其中，被z(k) a1z(k 1) a2z(k 2) b1u(k 1) b2u(k 2) v(k) (8.5 設(shè)輸入信號的取值是從 k =1 到 k =16 的辨識參數(shù) ?LS 、觀測矩陣 zL、HL 的表達式為?LSa1a2b1b2zLz(3)z(4)z(16)HLz(2)z(3)z(15)z(1)z(2)z(14

12、)u(2)u(3)u(1)u(2)(8.6u(15)u(14)/ 22程序框圖如圖 8.3 所示。賦輸入信號初值 uMatlab 仿真程序如下： %二階系統(tǒng)的最小二乘一次完成算法辨識程序， u=-1,1,-1,1,1,1,1,-1,-1,-1定,1義,-1輸,-出1,觀1,測1值。的長 % z=zeros(1,16。 % 定義輸出觀測值的長度 for k=3:16z(k=1.5*z(k-1-0.7*z(k-2+u(k-1+0.5*u(k-2 end subplot(3,1,1 % 畫三行一列圖形窗口中的第一個圖形 stem(u %畫輸入信號 u的徑線根圖G據(jù)(形z式1 (3.27計算參數(shù) ?L

13、S subplot(3,1,2 % 畫三行一列圖形窗口中的第二個圖形 i=1:1:16 。 %橫坐標范圍從是?LS中1分到離1出6并，顯步示長出為被辨1識 參數(shù) a1, a2, b1, b2 plot(i,z % 圖形的橫坐標是采樣時刻 i, 縱坐標是輸出觀測值 z, 圖形格式為連續(xù)曲線 subplot(3,1,3 % 畫三行一列圖形窗口中的 stem(z,grid on % 畫出輸出觀測值 z 的徑線圖形，并顯示坐標網(wǎng)格 u,z %顯示輸入信號和輸出觀圖測8信.3號最小二乘一次完成算法程序框圖 %L=14 % 數(shù)據(jù)長度 并統(tǒng)計辨算系識統(tǒng)的的輸輸入出信值號為一個周期的M 序列畫出輸入和輸出觀測

14、值的圖形。 % 用理想輸出值作為觀測值圖 3.7 最小二乘遞推算法辨識三停機個圖形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 。 -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

15、 -z(13 u(14 u(13 。 -z(15 -z(14 u(15 u(14 % 給樣本矩陣 HL 賦值ZL=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 % 給 樣本矩陣 zL 賦值 %Calculating Parametersc1=HL*HL 。 c2=inv(c1 。 c3=HL*ZL 。 c=c2*c3 % 計算并顯示 ?LS%Display Parametersa1=c(1, a2=c(2, b1=c(3,b2=c(4 % 從 ?LS 中分離出并顯示 a1

16、、 a2、 b1、 b2 %End 程序運行結(jié)果：1u = -1 ，1， -1，1，1，1，1，-1，-1，-1，1，-1，-1，1，z = 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.5800，-2.5185HL =00 1.0000 -1.0000-0.50000 -1.0000 1.0000-0.2500-0.50001.0000-1.0000-0.5250-0.25001.00001.0000-2.1125-0.52501.00001.00

17、00-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.0000/ 220.9386-3.26701.0000-1.00003.19490.9386-1.00001.00004.63523.1949-1.0000-1.00006.21654.63521.0000-1.00005.58006.21651.00001.0000ZL = 0.5000 ，0.2500，0.5250， 2.1125，4.3012，6.4731，6.1988，3.2

18、670，-0.9386，-3.1949， -4.6352，- 6.2165， -5.5800， -2.5185 Tc = -1.5000 ，0.7000，1.0000，0.5000 Ta1 = -1.5000a2 = 0.7000b1 = 1.0000b2 =0.500010-1100-10REMENTS OF OUTPUTWAVE OF MEASU246810121416ST EM PLO T OF MEASUREMENTS OF OUTPUT100-1002468 10 1214 16圖 8.4 最小二乘一次完成算法仿真實例中輸入信號和輸出觀測值從仿真結(jié)果表 8.4 可以看出，因為所用的輸

19、出觀測值沒有任何噪聲成分，所以辨識結(jié)果也無任何誤差。8.4 最小二乘遞推算法的表產(chǎn)生8.1 最小二乘一次完成算法的辨識結(jié)果 考慮圖 8.5 參所示數(shù)的仿真a對1 象 ， a2b1b2真 值 -1.5 0.7 1.0 0.5圖中， v(k)是服估從計值N的 不 0.相7 關(guān) 隨 機N1.(0噪z 1聲 ) 。 0且.5G (z 1) Be(zk1 ) ， N (z 1) D ( z 1) ， A+( z 1)C (z 1)u(kAy(zk1 ) 1 1.z5(ak1z 1 0.7z 2 C(z 1)圖 8.5 最小二乘遞B(推z算1法) 辨1識.0z實例1 0.5z 2D(z 1) 1選擇圖 8

20、.5 所示的辨識模型。仿真對象選擇如下的模型結(jié)構(gòu) z(k) a1z(k 1) a2z(k 2) b1u(k 1) b2u(k 2) v(k) (8.7 其中， v(k)是服從正態(tài)分布的白噪聲 N構(gòu)造 h (k ；加權(quán)陣取單位陣 L I ；計算 K(k 、和 P(k ，計算各次參數(shù)辨識的相對誤差，精度滿足要求 后停機。最小二乘遞推算法辨識的 Malab 程序流程如圖 8.6 所示。圖 8.6 最小二乘遞推算法辨識的工作間M清al零ab6.0 程序流程圖/ 22給 M 序列的長度 L 和移位寄存器的輸入賦初第四個移位寄存器的輸出取反，并將幅值變?yōu)?0.03 得到辨識系統(tǒng)的輸入信號 下面給出具體程序

21、。%最小二乘遞推算法辨識程序 , clear % 清理工作間變量L=15。 % M 序列的周期 y1=1。y2=1。y3=1。y4=0 。 %四個移位寄存器的輸出初始值 for i=1:L 。 %開始循環(huán)，長度為 Lx1=xor(y3,y4 。 %第一個移位寄存器的輸入是第三個與第四個移位寄存器的輸出的“或”x2=y1 。 %第二個移位寄存器的輸入是第一個移位寄存器的輸出x3=y2 。 %第三個移位寄存器的輸入是第二個移位寄存器的輸出x4=y3 。 %第四個移位寄存器的輸入是第三個移位寄存器的輸出y(i=y4 。 %取出第四個移位寄存器的幅值為 0和 1 的輸出信號，即 M 序列if y(i0

22、.5,u(i=-0.03 。 % 如果 M 序列的值為 1, 辨識的輸入信號取 “-0.03 ”else u(i=0.03 。 %如果 M 序列的值為 0, 辨識的輸入信號取 “ 0.03 ”end % 小循環(huán)結(jié)束y1=x1 。y2=x2 。y3=x3 。y4=x4 。 %為下一次的輸入信號做準備end %大循環(huán)結(jié)束，產(chǎn)生輸入信號 ufigure(1 。 %第一個圖形 stem(u,grid on % 顯示出輸入信號徑線圖并給圖形加上網(wǎng)格 z(2=0。z(1=0。 %設(shè) z 的前兩個初始值為零 for k=3:15 。 %循環(huán)變量從 3 到 15z(k=1.5*z(k-1-0.7*z(k-2+

23、u(k-1+0.5*u(k-2 。 % 輸出采樣信號 end%RLS 遞推最小二乘辨識c0=0.001 0.001 0.001 0.001 。 %直接給出被辨識參數(shù)的初始值 ,即一個充分小的實向量 p0=106*eye(4,4 。 %直接給出初始狀態(tài) P0，即一個充分大的實數(shù)單位矩陣E=0.000000005 。 %取相對誤差 E=0.000000005 c=c0,zeros(4,14 。 %被辨識參數(shù)矩陣的初始值及大小 e=zeros(4,15。 % 相對誤差的初始值及大小 for k=3:15 。 %開始求 Kh1=-z(k-1,-z(k-2,u(k-1,u(k-2 。 x=h1*p0*h

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

25、1,: 。 a2=c(2,:。 b1=c(3,: 。 b2=c(4,: 。 ea1=e(1,:。 ea2=e(2,:。 eb1=e(3,:。 eb2=e(4,: 。 figure(2 。 %第二個圖形i=1:15 。 %橫坐標從 1 到 15 plot(i,a1,r,i,a2,:,i,b1,g,i,b2,: % 畫出 a1，a2，b1，b2 的各次辨識結(jié)果 title(Parameter Identification with Recursive Least Squares Method%圖形標題 figure(3 。 %第三個圖形 i=1:15 。 %橫坐標從 1 到 15 plot(i,e

26、a1,r,i,ea2,g,i,eb1,b,i,eb2,r: % 畫出 a1，a2，b1，b2 的各次辨識結(jié)果的收斂情況 title(Identification Precision % 圖形標題程序運行結(jié)果： c =0.0010 0 0.0010 -0.4984 -1.2328 -1.4951 -1.4962 -1.4991 -1.4998 -1.49990.0010 0 0.0010 0.0010 -0.2350 0.6913 0.6941 0.6990 0.6998 0.69990.0010 0 0.2509 1.2497 1.0665 1.0017 1.0020 1.0002 0.999

27、9 0.99980.0010 0 -0.2489 0.7500 0.5668 0.5020 0.5016 0.5008 0.5002 0.5002-1.5000 -1.5000-1.5000 -1.4999 -1.49990.7000 0.70000.7000 0.7000 -0.70000.9999 0.9999 0.99990.99990.99990.5000 0.5000 0.50000.50000.5000e =0 0 0 -499.4200 1.4734 0.2128 0.0007 0.0020 0.0004 0.00000 0 0 0 -235.9916 -3.9416 0.004

28、2 0.0070 0.0012 0.00010 0 249.8612 3.9816 -0.1466 -0.0607 0.0003 -0.0018 -0.0003 -0.00010 0 -249.8612 -4.0136 -0.2443 -0.1143 -0.0007 -0.0016 -0.0012 -0.00010.0001 0.0000 -0.0000 -0.0000 0.00000.0001 -0.0000 0.0000 0.0000 0.00000.0001 0.0000 0.0000 0.0000 -0.0000-0.0004 0.0000 -0.0000 0.0000 -0.0000

29、Fig.1 Input Signal圖 8.7 最小二乘遞推算法的參數(shù)辨識數(shù)/ 221.5-1.50Fig. 2 Parameter Ident5 10 ification with Recursive Least Squares Method15Convergence of b1Convergence of a2Convergence of b2Convergence of a1圖 8.8 最小二乘遞推算法的參數(shù)辨識結(jié)果3002001000-100-200-300-400-500Identification Precision0 5 10 15Fig.3 Identification Err

30、or圖 8.9 最小二乘遞推算法的參數(shù)辨識結(jié)果收斂情況 表 8.2 最小二乘遞推算法的辨識結(jié)果參數(shù) a1a2b1b2真值-1.50.71.00.5估計值-1.49990.70.99990.5000仿真結(jié)果表明，大約遞推到第十步時，參數(shù)辨識的結(jié)果基本達到穩(wěn)定狀態(tài)，即 a1=-1.4999, a2= 0.7000, b1=0.9999, b2=0.5000 。此時，參數(shù)的相對變化量E 0.000000005。從整個辨識過程來看，精度的要求直接影響辨識的速度。雖然最終的精度可以達到很小，但開始階段的相對誤差會很大，從圖 8.7(3 圖形中可看出， 參數(shù)的最大相對誤差會達到三為數(shù)廣義最小二乘算法的產(chǎn)生

31、 廣義最小二乘法所用的模型是：1 1 1A(z 1)Z(k) B(z 1)u(k)1 e( k) 8.9)c(z )z(k) a1z(k 1) a2z(k 2) b1u(k 1) b2u(k 2) e(k) e(k) c1 e(k 1) c2e(k 2) (k)其中我們選定模型參數(shù)： a1=1.5,a2=-0.7,b1=1,b2=0.5,c1=0,c2=0 廣義最小二乘遞推算法的計算步驟：?(0)(充分小的實向量 )Pf(0) a2I (a為充分大的數(shù) )1.給定初始條件 f?e (0) 0Pe (0) I/ 222.利用式 zf (k) C(z 11)z(k)，計算 zf (k)和uf (k

32、)； uf (k) C(z 1)z(k) f f3.利用式構(gòu)造 hf(k) ； a1, ,ana,b1, ,bnb hf(k) zf(k 1), , zf(k na ),uf (k 1), ,uf (k nb)?(k) ?(k 1) K f (k)zf (k) hf (k) ?(k 1)4.利用式 K f (k) Pf (k 1)hf(k) hf (k)Pf(k 1)hf (k) 1 1 遞推計算 ?(k)； Pf (k) I K f(k)hf (k)Pf(k 1)5.利用 e?(k) z(k) h (k) ?(k) 和 h(k) z(k 1), , z(k na ),u(k 1), ,u(k

33、 nb) ， 計算 e?(k) ；6.根據(jù) he(k) e?(k 1), , e?(k nc) 來構(gòu)造 he(k) ； ?e(k) ?e(k 1) K e (k) e?(k ) he(k) ?e(k 1)7.利用 Ke(k) Pe(k 1)he (k) he (k)Pe (k 1)he (k) 1 1 Pe(k) I K e(k)he(k)Pe(k 1)廣義最小二乘法程序代碼 clear % 清理工作間變量L=55。 % M 序列的周期y1=1。y2=1。y3=1。y4=0 。 %四個移位寄存器的輸出初始值for i=1:L 。 %開始循環(huán)，長度為 Lx1=xor(y3,y4 。 %第一個移位

34、寄存器的輸入是第三個與第四個移位寄存器的輸出的“或”x2=y1 。 %第二個移位寄存器的輸入是第一個移位寄存器的輸出x3=y2 。 %第三個移位寄存器的輸入是第二個移位寄存器的輸出x4=y3 。 %第四個移位寄存器的輸入是第三個移位寄存器的輸出y(i=y4 。 %取出第四個移位寄存器的幅值為 0和 1 的輸出信號，即 M 序列if y(i0.5,u(i=-1 。 %如果 M 序列的值為 1, 辨識的輸入信號取 “-1”else u(i=1 。 %如果 M 序列的值為 0, 辨識的輸入信號取 “ 1”end % 小循環(huán)結(jié)束y1=x1 。y2=x2 。y3=x3 。y4=x4 。 %為下一次的輸入

35、信號做準備end %大循環(huán)結(jié)束，產(chǎn)生輸入信號 u z(2=0。z(1=0。 %設(shè) z 的前四個初始值為零 for k=3:45 。 %循環(huán)變量從 5 到 45z(k=1.5*z(k-1-0.7*z(k-2+1*u(k-1+0.5*u(k-2 。 % 輸出采樣信號end%RGLS 廣義最小二乘辨識c0=0.001 0.001 0.001 0.001 。 %直接給出被辨識參數(shù)的初始值 ,即一個充分小的實向量 pf0=106*eye(4,4 。 %直接給出初始狀態(tài) P0，即一個充分大的實數(shù)單位矩陣 ce0=0.001 0.001 。pe0=eye(2,2。/ 22 c=c0,zeros(4,44 。

36、 %被辨識參數(shù)矩陣的初始值及大小 ce=ce0,zeros(2,44 。e=zeros(4,45。 % 相對誤差的初始值及大小 ee=zeros(2,45。for k=3:45 。 % 開始求 K zf(k=z(k+0*z(k-1+0*z(k-2 。 uf(k=u(k+0*u(k-1+0*u(k-2 。hf1=-zf(k-1,-zf(k-2,uf(k-1,uf(k-2 。 x=hf1*pf0*hf1+1 。 x1=inv(x 。 % 開始求 K(k k1=pf0*hf1*x1 。 %求出 K 的值d1=zf(k-hf1*c0 。 c1=c0+k1*d1 。 % 求被辨識參數(shù) ce1=c1-c0

37、 。 %求參數(shù)當前值與上一次的值的差值 e2=e1./c0。 % 求參數(shù)的相對變化e(:,k=e2 。 %把當前相對變化的列向量加入誤差矩陣的最后一列c0=c1 。 %新獲得的參數(shù)作為下一次遞推的舊參數(shù)c(:,k=c1 。 %把辨識參數(shù) c 列向量加入辨識參數(shù)矩陣的最后一列 pf1=pf0-k1*hf1*pf0 。 %求出 p(k 的值pf0=pf1 。 % 給下次用 h1=-z(k-1,-z(k-2,u(k-1,u(k-2 。ee(k=z(k-h1*c1 。 he1=-ee(k-1,-ee(k-2 。x=he1*pe0*he1+1 。 x1=inv(x 。k1=pe0*he1*x1 。 d1

38、=ee(k-he1*ce0 。ce1=ce0+k1*d1 。 pe1=pe0-k1*he1*pe0 。ce0=ce1。 ce(:,k=ce1 。pe0=pe1。end %大循環(huán)結(jié)束% 顯示被辨識參數(shù)及其誤差 (收斂 情況%分離參數(shù)a1=c(1,: 。 a2=c(2,: 。 b1=c(3,: 。 b2=c(4,: 。 c1=ce(1,: 。c2=ce(2,: 。 ea1=e(1,:。 ea2=e(2,:。 eb1=e(3,:。 eb2=e(4,: 。figure(2 。 %第二個圖形i=1:45 。 %橫坐標從 1到 25plot(i,a1,r,i,a2,:,i,b1,b,i,b2,:,i,c1

39、,b,i,c2,: % 畫出 a1， a2,b1， b2,c1,c2 的各次辨識結(jié)果 title( 參數(shù)變化曲線 % 圖形標題figure(3 。 %第三個圖形i=1:45 。 %橫坐標從 1到 25plot(i,ea1,r,i,ea2,k:,i,eb1,g,i,eb2,c: % 畫出 a1， a2，b1， b2,的各次辨識結(jié)果的收斂情況 title(誤差曲線 %圖形標題/ 22圖 8.10 廣義最小二乘法辨識參數(shù)變化曲線圖 8.11 廣義最小二乘法辨識誤差變化曲線結(jié)果分析：廣義最小二乘法的基本思想是基于對數(shù)據(jù)先進行一次濾波預(yù)處理，然后利用普通的最小二乘 法對濾波后的數(shù)據(jù)進行辨識。辨識結(jié)果 a

40、1=1.5016,a2=-0.7010,b1=1.0120,b2=0.4986,c1=0.0040,c2=0.0027 與最 小二乘法遞推算法相比較，可以發(fā)現(xiàn)它們的結(jié)果基本上是一致的。這是因為本例的仿真對象相當于=1 。這種對象利用最小二乘法已可獲得無偏一致估計，但廣義最小二乘法同時又能給出噪聲模型的參數(shù)估計。增廣最小二乘遞推算法的產(chǎn)生程序考慮圖 8.12 所示的仿真對象，圖中， v(k)是服從 N0,1 )分布的不相關(guān)隨機噪聲。且G(z 1)B(z 1)，A( z 1)N (z 1)D (z 1 ) ，C (z 1)A(z 1) 1 1.5a1z 1 0.7z 2 C(z 1) B(z 1)

41、 1.0z 1 0.5z 2D(z 1) 1 z 1 0.2z 2模型結(jié)構(gòu)選用如下形式：z(k)a1z(k1)a2(k2)b1u(k 1)b2u(k2) v(k)d1v(k 1)d2(v(k 2)增廣最小二乘遞推辨識程序流程如圖 8.13 所示。圖 8.13 增廣最小二乘遞推工算作法間辨清零識的 Malab6.0 程序流程圖Matlab 程序如下。 %增廣最小二乘遞推辨識程序， L=60。 %4 位移位寄存器產(chǎn)生的y1=1。y2=1。y3=1。y4=0 。 %4 for i=1:L 。M 序列的周期用 4 位移位寄存器產(chǎn)生幅值為 1 的輸入信個移位寄存器的輸出初始值產(chǎn)生隨機噪聲信號x1=xor

42、(y3,y4 。 % 第一個移位寄存器的輸入x2=y1 。 %第二個移位寄存器的輸入信畫號出輸入信號徑線圖形及隨機噪聲圖 x3=y2 。 %第三個移位寄存器的輸入信號x4=y3 。 %第四個移位寄存器的輸入信號產(chǎn)生輸出采樣信號y(i=y4 。 %第四個移位寄存器的輸出信號， M 序列, 幅值 0和1 ， 給被辨識參數(shù) 和 P 賦初值 if y(i0.5,u(i=-1 。 %M 序列的值為 1 時 ,辨識的輸入信號取 “-1 else u(i=1 。end%M 序列的值為 0 時,辨識的輸入信號取 “ 1”按照式 (3.80的第二式遞推計算 K(ky1=x1 。 y2=x2 end figure

43、(1 。 %畫。y3=x3 。y4=x4 。%為下一次的輸入信號準備按 照式(3.80的第一式 遞推計算第一個圖形% 畫第一個圖形的第一個按子照圖式subplot(2,1,1 。stem(u,grid on % 畫出 M 序列輸入信號 v=randn(1,60 。按照式 (3.80 的第三式遞推計算 P(k計算系統(tǒng)的實際輸出響應(yīng)及模型的響% 產(chǎn)生一組 60 個正態(tài)分布的隨機噪聲 %畫第一個圖形的第二按 % 畫出隨機噪聲信號號subplot(2,1,2 。 plot(v,grid on 。u ，v %顯示輸入信號和噪聲信 z=zeros(7,60 。 zs=zeros(7,60 。 zmN=ze

44、ros(7,60 。 出矩陣、不考慮噪聲時模型輸出矩陣、模型輸出矩陣的大小顯示被辨識參數(shù)、辨識精度、輸出采樣、系統(tǒng)實際輸出、模型輸 z(2=0 。 z(1=0 。 zs(2=0 。 zs(1=0 。 zm(2=0 。 zm(1=0 。 噪聲時系統(tǒng)輸出、不考慮噪聲時模型輸出、.67計算被辨識參數(shù)的相對變化量Y收斂滿足要s(7,60 。 %輸出采樣矩陣、不考慮噪聲時系統(tǒng)輸zmd(2=0 。 zmd(1=0 。 % 輸出采樣、不考慮型辨輸識參出初值的各次估計值及其誤,即一個充分小的實向量c0=0.001 0.001 0.001 0.001 0.001 0.001 0.001 。 % 直接給出被辨識參

45、數(shù)的初始值 p0=106*eye(7,7 。 %直接給出初始狀畫態(tài)出系統(tǒng)P0的，實即際一輸出個、充輸分出大采的樣實及數(shù)模型單輸位矩陣 E=0. 。 %相對誤差 E=0.000000005c=c0,zeros(7,14 。 %被辨識參數(shù)矩陣的初始值及大小 停機 e=zeros(7,15。 % 相對誤差的初始值及大小。 % 系統(tǒng)在 M 序列輸入下的輸出采for k=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/ 22樣信號h1=-z(k-1,-z(k-2,u(k-1,u(k-2,v(k,v(k-1,

46、v(k-2 。 % 為求 K(k 作準備 x=h1*p0*h1+1 。 x1=inv(x 。 k1=p0*h1*x1 。 %Kd1=z(k-h1*c0 。 c1=c0+k1*d1 。 % 辨識參數(shù) c zs(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ù)誤差的

47、相對變化e(:,k=e2 。c0=c1 。 %給下一次用c(:,k=c1 。 % 把遞推出的辨識參數(shù) c 的列向量加入辨識參數(shù)矩陣p1=p0-k1*k1*h1*p0*h1+1 。 %find p(kp0=p1 。 % 給下次用if e2 。 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=e(2,:。 eb1=e(3,: 。 eb2=e(4,: 。 %分離出 a1、 a2、 b1、 b2 的收斂情況 ed1=e(5

48、,:。 ed2=e(6,:。 ed3=e(7,: 。 %分離出 d1 、d2 、 d3的收斂情況 figure(2 。 %畫第二個圖形i=1:60 。plot(i,a1,r,i,a2,r:,i,b1,b,i,b2,b:,i,d1,g,i,d2,g:,i,d3,g+ % 畫出各個被辨識參數(shù) title(Parameter Identification with Recursive Least Squares Method % 標題 figure(3 。 i=1:60 。 %畫出第三個圖形plot(i,ea1,r,i,ea2,r:,i,eb1,b,i,eb2,b:,i,ed1,g,i,ed2,g:

49、,i,ed2,r+ % 畫出各個參數(shù)收斂情況 title(Identification Precision % 標題figure(4 。 subplot(4,1,1 。 % 畫出第四個圖形，第一個子圖i=1:60 。 plot(i,zs(i,r % 畫出被辨識系統(tǒng)在沒有噪聲情況下的實際輸出響應(yīng)subplot(4,1,2 。 i=1:60 。 plot(i,z(i,g % 第二個子圖，畫出被辨識系統(tǒng)的采樣輸出響應(yīng) subplot(4,1,3 。 i=1:60 。 plot(i,zm(i,b % 第三個子圖，畫出模型含有噪聲的輸出響應(yīng) subplot(4,1,4 。 i=1:60 。 plot(i

50、,zs(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 = 1.1418,1.5519,1.3836,-0.7581,0.4427,0.9111,-1.0741,0.2018,0.7629,-1.2882,-0.9530,0.7782,-0.0063,0.5245

51、, 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.2492/ 220.3692,0.1792,-0.0373,-1.6033,0.3394,-0.1311,0.4852,0.5988,-0.0860,0.3253,-0.3351,-

52、0.3224,-0.3824,-0.9534,0.2336,1.2352,-0.5785,-0.5015c =0.0010 0 0.0010 -0.2789-0.9025 -0.9390 -0.8620-1.1236 -1.5000-1.5000-1.50000.0010 0 0.00100.0010 -0.0734 -0.2661-0.1199 0.3259 0.7000 0.70000.70000.0010 0 0.05920.45600.53000.52930.94580.54021.00001.00001.00000.0010 00.05720.81860.85000.86391.11

53、571.29570.50000.50000.50000.0010 00.07960.73420.5417 0.5516 0.38470.24591.00001.00001.00000.0010 0-0.0894-0.5981-0.3420-0.4456-0.4044-0.8443-1.0000-1.0000 -1.00000.0010 0-0.0655-0.7795-0.8572-0.7412-0.4506-0.19740.20000.2000 0.2000e =0 0 0 -279.9426 2.2355 0.0404 -0.0820 0.3035 0.3350 0.0000 0.00000

54、 0 0 0 -74.3764 2.6263 -0.5495 -3.7189 1.1479 -0.0000 0.000000 58.2212 6.7007 0.1622 -0.0013 0.7868 -0.4289 0.8512 0.0000 0.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.0000 0.000000-90.35565.6936-0.42820.3029-0.09251.08800.

55、1844-0.0000 0.000000-66.475410.90590.0996-0.1354-0.3920-0.5619-2.01300.0000 0.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

56、,-2.3342,-2.9824,-3.9807,-5.5271,-6.6924,-8.7267,-6.5221,-2.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.093

57、2,-2.1863,-7.9830,-5.8500,-0.5991, -1.6810,2.3509,2.9533,0.3337,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

58、.1232,5.2248,1.9393,-4.2395,-3.3718,-1.9258, 6.9389,7.3995,1.2763 zm =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.982

59、4,-3.9807,-5.5271,-6.6924,-8.7267,-6.5221,-2.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.0478Fig.1(2) Stochastic Noises v圖 8.14 辨識系統(tǒng)在沒有噪聲情況下的實際輸出響應(yīng)/ 221.5b1&d1圖 8.15 被辨識系統(tǒng)的采樣輸出響應(yīng)100

60、0-500Error of a2Error of b2Error of a1Errorm of d1Error of b1Error of d2Error of d350000 1020 30 4050 60Fig.3 Identification Error圖 8.16 模型含有噪聲的輸出響應(yīng)100-10100-10100-10200-20SystemResponse Without Noises0102030405060SystemResponse With Noises0102030405060Model Response With Noises0102030405060Model Re