1、實用標準文案系統辨識大作業學號:12051124精彩文檔班級:自動化1班張青姓名:信息與控制工程學院自動化系2015-07-11®k)，由g(k)，參照講第一題模仿index2，搭建對象，由相關分析法，獲得脈沖響應序列 義，獲得系統的脈沖傳遞函數 G(z)和傳遞函數G(s);應用最小二乘辨識,獲得脈沖響應序列6"同圖顯示兩種方法的辨識效果圖；應用相關最小二乘法，擬合對象的差分方程模型;(可以用辨識工具箱)辨識模型構建時變對象，用最小二乘法和帶遺忘因子的最小二乘法， 的參數，比較兩種方法的辨識效果差異答：根據index2搭建結構框圖:相關分析法：利用結構框圖得到 UY和tou

2、t矗屯山1fl Ba&ed國 Select cdtaFR UV '：2015<2 cJoublo12 11 345Jj15SO£25-34 J 55a二 13455-14.144-.3471554.773655旳275-2 眉02Vdrid'bit Ecitor - UY其中的 U就是題目中要求得出的m 口 XNameVslutSuv*201S lout<100*' n r X畫倒國斑M序列，根據結構框圖可知序列的周期是n4Np 212115。在comma nd window 中輸入下列指令，既可以得到 脈沖相應序列g(k)aa=5;N N

3、PP=15;ts=2;RR=o nes(15)+eye(15);for i=15:-1:1UU(16-i,:)=UY(16+i:30+i,1)'endYY=UY(31:45,2);GG=RR*UU*YY/aa*aa*(NN PP +1)*ts;plot(0:2:29,GG)hold onstem(0:2:29,GG,'filled')Grid;title('脈沖序列 g( ©')最小二乘法建模的響應序列04030 250.2D 16D.0&10由于是二階水箱系統，可以假設系統的傳遞函數為bobiSG（S）1 a1S2，已知g(),求b0

4、,b1,a1, a2g()求得的各階矩,已知G( S)的結構，用長除法求得G(s)的s展開式，其系數等于從然后求G(s)的參數。得到結果:a1 =-1.1561 a2 =0.4283 bO =-0.0028 b1 =0.2961在comma nd window中輸入下列指令得到傳遞函數:>> nu>= 0, £56L -0. 0fl23 ;(leii= (0. 4233 -1. 1551 1 ；3y=-t'£ (nuiij den)Transfer function:0. 2961 s - 0. 002SC. 42Q3 S 2 - U 150 3

5、+ 1最小二乘一次算法相關參數%最小二乘法一次完成算法M=UY(:,1);z=UY(:,2);H=zeros(100,4);for i=1:100H(i,1)=-z(i+1);H(i,2)=-z(i);H(i,3)=M(i+1);H(i,4)=M(i);endEstimated nv(H'*H)*H'*(z(3:102) %結束得到相關系數為：Estimate =-0.78660.13880.57070.3115帶遺忘因子最小二乘法:%帶遺忘因子最小二乘法程序M=UY(:,1);z=UY(:,2);P=1000*eye(5);%Theta=zeros(5,200); %Thet

6、a(:,1)=0;0;0;0;0;K=zeros(4,400);%K=10;10;10;10;10;lamda=0.99;%遺忘因數 for i=3:201 h=-z(i-1);-z(i-2);M(i);M(i-1);M(i-2);K=P *h* in v(h'* P*h+lamda);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2);P=(eye(5)-K*h')* P/lamda;endi=1:200;figure plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta

7、(4,:),i,Theta(5,:) title（'帶遺忘因子最小二乘法'） grid%結束Estimate 可由仿真圖得出，可知兩種方法參數確定十分接近。兩種方法的辨識效果差異，用最小二乘法辨識出的圖形更采用相關分析法和最小二乘法辨識出的脈沖響應圖形可看出 像脈沖響應，他在最后無偏差，衰減為零，其可實現無偏差估計。而相關分析法圖形雖與相 關分析的相近，但它最后有偏差。帶遺忘因子的遞推最小二乘辨識的參數平均值可隨著實際 參數變化而突變。但他對噪聲比較敏感，只是參數圍繞參數實際值上下波動，而辨識出參數 的平均值接近實際值，而且比其他方法更加接近，并且可以防止數據飽和。第二題設SI

8、SO系統差分方程為(40 分)y(k)=aiy(k 1) a2y(k2)b1u(k 1) b2u(k 2)(k)辨識參數向量為=a1 a2lbb2 T ,輸入輸出數據詳見數據文件uy01.txt uy03.txt 。(k)為噪聲方差各異的白噪聲或有色噪聲。試求解:用最小二乘及遞推最小二乘法估計 用輔助變量法及其遞推算法估計 用廣義最小二乘法及其遞推算法估計 用夏氏偏差修正法、夏氏改良法及其遞推算法估計 用增廣矩陣法估計6)用極大似然法估計分析噪聲 (k)特性;答：1.基本最小二乘法：%首先將TXT中的數據裝載在 Matlab中:clear;clc;Wbs = 0;while Wbs<1d

9、isp('從下列文件中選擇所需加載文件： file name = inp ut('輸入所需打開文件名uy1.txt uy2.txt uy3.txt')endWbs,Message = fopen( file name,'。;:','s');打開文件;Data = fscan f(Wbs,'%g %g',2 in f);%生成數據矩陣Data = Data'%最小二乘法，取 bO=On = 2;%根據題目，本系統為2階系統L = len gth(Data);%辨識所需數據的總長度N = L-n;% 構造測量矩陣的數據

10、長度FIA = zeros(N,2* n);%構造測量矩陣for i = 1:N%測量矩陣賦值for l = 1:n*2if(l> n)FIA(i,l) = Data( n+2+i-l,1);elseFIA(i,l) = -Data( n+i-1,2);endend endY = Data( n+1: n+N,2);%輸出數據矩陣thita = in v(FIA'*FIA)*FIA'*Y;%計算參數矩陣')disp('使用最小二乘算法所得結果如下：%辨識參數數據輸出如下： lsa1 = thita(1),lsa2 = thita (2),lsb1 = th

11、ita(3),lsb2 = thita(4)實驗結果:Uy1.txtuy2.txtuy3.txtIsal =Isal =Isa】二i&5&-i. 1579-klsa2 =lsa2 =0.S51?0.S5130. 3639Isbl =Isbl =Isbl =0. 42500, 73440. 7451lsb2 =lsb2 =lsb2 =0. 56190. 57736 55432.遞推最小二乘法:%首先將TXT中的數據裝載在Matlab中:clear;clc;Wbs = 0;while Wbs<1disp('從下列文件中選擇所需加載文件： file name = inp

12、 ut('輸入所需打開文件名uy1.txt uy2.txt uy3.txt')endWbs,Message = fopen( file name,'r');%：','s');打開文件;Data = fscan f(Wbs,'%g %g',2 in f);%生成數據矩陣Data = Data'%遞推最小二乘法n = 2;%根據題目已知該系統階數為L = len gth(Data);%為辨識參數向量賦初值，這里均取為0.001for a = 1:1: n*2thita0(a,1) = 0.001;end%直接給出矩陣P

13、n的初始狀態P0P0 = 109*eye (n *2, n*2);thita = thita0,zeros (n *2,L);%需辨識參數矩陣的初值及大小for i = n+1:1:L %這里i從第n+1個數開始for m = 1:n*2 % 輸入矩陣賦值if(m<=n)FIA1(m,1) = -Data(i-m,2);endelseendendFIA1(m,1) = Data(i-m+2,1);X = FIA1'*P0*FIA1+1;%弓I入中間變量 X為K遞推公式中的求逆項,1維K1 = P 0*FIA1/X;D1 = Data(i,2)-FIA1'*thita0;P

14、1 = P0-K1*K1'*X;% 求出 P(K)矩陣P0 = P1;%作為下次迭代初值thita1 = thita0+K1*D1;%估值補償公式，thita1為求被辨識的參數向量thita0 = thita1;%所得的參數向量作為下次迭代初值disp('遞推最小二乘法所得結果如下:')顯示被辨識參數Dta仁thita1(1),Dta2=thita1(2),Dtb仁thital (3) ,Dtb2=thita1(4) %實驗結果:Dtal =It al 二Dt al =-U 1555-L 4579-1.4&3tl =Dt aZ =551'9'0.

15、8639Etb 1 =Etbl =Dthl =0. 42500. 73410* 7451Etb2 =Etb2 =0. 5019'0. 57730.5E43Uy1.txtuy2.txtuy3.txt3.輔助變量法：%首先將TXT中的數據裝載在 Matlab中:clear;clc;Wbs = 0;while Wbs<1disp('從下列文件中選擇所需加載文件： file name = inp ut('輸入所需打開文件名uy1.txt uy2.txt uy3.txt')endWbs,Message = fopen( file name,'。;:'

16、,'s');打開文件;Data = fscan f(Wbs,'%g %g',2 in f);%生成數據矩陣Data = Data'%使用遞推輔助變量法對 2階系統進行辨識，取 b0=0n = 2;L = len gth(Data);N = Ln;%分別取出輸入、輸出數據u,yU = Data(1:L,1);Y = Data(1:L,2);%首先使用最小二乘法，粗略得到初始的被辨識參數向量fzOL = -Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2);Y1 = Data(3:L,2);%構造輸出向量，也用于輔助變量法的輸出向量得到

17、初始參數向量thita = in v(fzOL'*fzOL)*fzOL'*Y1;%獲得初值，以下使用輔助變量法 for i=1:400Yf = fzOL*thita;% 輔助模型輸出%構造輔助變量矩陣，改變新的矩陣中的輸出量，輸入量不變Zfz = fzOL;Zfz(2:N,1) = -Yf(1:(N-1),1);Zfz(3:N,2) = -Yf(1:(N-2),1);thita = in v(Zfz'*fzOL)*Zfz'*Y1;%求取被估計參數的輔助變量估值enddisp('使用輔助變量法，得到如下辨識結果:')Fza1=thita(1),Fz

18、a2=thita(2),Fzb1=thita(3),Fzb2=thita(4)實驗結果:Uy1.txtuy2.txtuy3.txtFcal =Fzal =fsal 二-I. 475&-1. 4479-L 4662=0,8 翎30.essT0, B620Fzbl =Fzbl =Fzbl 二0. 43C40.73C20. 7441Fib2 =r3b2 =0. 55S90. 53210. S5544.遞推輔助變量法：%首先將TXT中的數據裝載在 Matlab中:uy1.txt uy2.txt uy3.txt'):','s');打開文件;clear;clc;Wb

19、s = 0;while Wbs<1disp('從下列文件中選擇所需加載文件： file name = inp ut('輸入所需打開文件名Wbs,Message = fopen( file name,'。;生成數據矩陣endData = fscan f(Wbs,'%g %g',2 in f);%Data = Data'%使用遞推輔助變量法對2階系統進行辨識，取 b0=0 n = 2;L = len gth(Data);%分別取出輸入、輸出數據u,yU = Data(1:L,1);Y = Data(1:L,2);N = 498;n0 = 100

20、;%取前100個數據利用最小二乘辨識FIA = -Y(2: n0-1),-Y(1: n0-2),U(2: n0-1),U(1: n0-2);Y1 = Data(3:n 0,2);%構造輸出向量，也用于輔助變量法的輸出向量thita = in v(FIA'*FIA)*FIA'*Y1;%得到初始參數向量Z =卜丫(2: n0-1) -Y(1: n0-2) U(2: n0-1) U(1: n0-2);thita = in v(Z'*FIA)*Z'*Y1;P0 = in v(Z'*FIA);%后面的數據計算，使用遞推輔助變量法 for n = n0+1:NY11

21、 = Y(1); Y(2); Z*thita;U( n-1) U(n-2);Zn = -Y11( n-1) -Y11( n-2)Z = Z; Zn;kesy = -Y( n-1); -Y( n-2); U( n-1); U( n-2);K = P0*Z n'/(1+kesy'* P0*Zn');P1 = P0 - K*kesy'* P0;P0 = P1;thita = thita + K*(Y (n) - kesy'*thita);enddisp('遞推輔助變量法辨識結果：');Dfa1=thita(1),Dfa2=thita (2) ,

22、Dfb1=thita (3),Dfb2=thita (4)實驗結果：Uy1.txtuy2.txtuy3.txt=HhI =Dfal 二-1,4671-1. 46bT-1.4617Df&2 =Dfa2 =0. 3624Q. B6250.B5&3DfbL =Dfbl =Dihl =0. 42910. 74650.7400Dfb2 =Dfb2 =n；b2 =a 55530. 56B50.55385.廣義最小二乘法：%首先將TXT中的數據裝載在 Matlab中:clear;clc;Wbs = 0;while Wbs<1disp('從下列文件中選擇所需加載文件： file

23、name = inp ut('輸入所需打開文件名uy1.txt uy2.txt uy3.txt')endWbs,Message = fopen( file name,'。;:','s');打開文件;Data = fscan f(Wbs,'%g %g',2 in f);%生成數據矩陣Data = Data'%使用廣義最小二乘法，辨識2階系統的參數n = 2;L = len gth(Data);N = L-n;for m=1:1:300%分別取出輸入、輸出數據 u,yU = Data(1:L,1); Y = Data(1:L,

24、2);%首先使用最小二乘法，求系統參數向量初值gyOL = -Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2);%LSZgy1 = Data(3:L,2);%輸出矩陣 y測量矩陣thita = in v(gyOL'*gyOL)*gyOL'*Zgy1;%參數向量初值Gya1=thita(1);Gya2=thita (2);Gyb1=thita (3) ;Gyb2=thita(4);%得到初值后，以下使用廣義最小二乘法進行辨識，由殘差代替噪聲誤差E(1) = Y(1);%k=1時，殘差的值E(2) = Y(2)+Gya1*Y(1)-Gyb1*U(1);%k=

25、2 for i = 3:N+n%殘差公式如下：時，殘差的值E(i)= Y(i)+Gya1*Y(i-1)+Gya2*Y(i-2)-Gyb1*U(i-1)-Gyb2*U(i-2);endE1 = E(n+1: n+N)：omiga(1:N,1) = -E( n:n+N-1)'omiga(1:N,2) = -E( n-1: n+N-2)'%fl表示由殘差代替噪聲項，而辨識得到的濾波器的估計參數fl=in v(omiga'*omiga)*omiga'*E1;%得到fl后，計算經過濾波的u , y數據；其中，k=1,2時:Data(1,1)=U(1);Data(1,2)=

26、Y(1);Data(2,2)=Y(2)+fl(1)*Y(1);Data(2,1)=U(2)+fl(1)*U(1);%k>=3 時:for k=3:N+nData(k,1) = U(k)+fl(1)*U(k-1)+fl (2) *U(k-2);Data(k,2) = Y(k)+fl(1)*Y(k-1)+fl(2)*Y(k-2);endenddisp('由廣義最小二乘法得到的辨識結果如下:')Gya仁thita(1),Gya2=thita (2),Gyb仁thita (3) ,Gyb2=thita (4)實驗結果:實驗結果:tyal =Gy al =G-yal =-1.469

27、4-1.-I.4030=Gy a 2 =G-ya2 =0. 36470.8&e&0.S607G-ybl =Gybl =&ybl =0* 42360, 71150.7425tyb2 =&yb2 =Gyb2 =Q 55490. 55030.5553Uy1.txtuy2.txtuy3.txt6.遞推廣義最小二乘法%首先將TXT中的數據裝載在 Matlab中:clear;clc;Wbs = 0;while Wbs<1disp('從下列文件中選擇所需加載文件：file name = inp ut('輸入所需打開文件名Wbs,Message = fop

28、en( file name,'。;uy1.txt uy2.txt uy3.txt'):','s');打開文件;endData = fscan f(Wbs,'%g %g',2 in f);%生成數據矩陣Data = Data'u1=Data(:,1);y1=Data(:,2);n=2;N=500-n;n0=50;%用前50組數據利用lsm估計出初值%構造矩陣phiphi(:,1)=-y1( n:n+n 0-1);phi(:,2)=-y1( n-1: n+n 0-2);phi(:,3)=u1( n:n+n0-1);phi(:,4)=u

29、1( n-1: n+n 0-2);y(:,1)=y1( n+1: n+n0);seta=(i nv(p hi'* phi)* phi'*y; seta0=seta;f0=0 0'p10=i nv( phi'* phi);p20=10A6*eye(2,2);for i=n 0:N-2y(1:i,1)=y1( n+1: n+i);u(1:i,1)=u1( n+1: n+i);%計算殘差e=y-p hi*seta0;%f估計omega=-e( n+i-2) -e(i+1-1)'K2=p20*omega* in v(1+omega'* p20*omega

30、);f1=f0+K2*(e( n+i-2)-omega'*f0);p2=p20-K2*omega'* p20;p20=p2;f0=f1;%濾波yy(1:i,1)=y1( n:n+i-1);yy(1:i,2)=y1( n-1: n+i-2);uu(1:i,1)=u1( n:(n+i-1);uu(1:i,2)=u1( n-1):( n+i-2);yg=y+yy*f1;%yg(k)=y(k)+f(1)*y(k-1)+f(2)*y(k-2)ug=u+uu*f1; % seta估計%ug(k)=u(k)+f(1)*u(k-1)+f(2)*u(k-2)ksai=-yg( n+i-2) -y

31、g( n+i-3) ug( n+i-2) ug( n+i-3)'K仁p 10*ksai* in v(1+ksai'* p10*ksai);% K(N+1)% P(N+1)值p1= p10-K1*ksar* p10;p10=p1;seta仁seta0+K1*(y1( n+i+1)-ksai'*seta0); seta0=seta1;% seta(N+1)值%構造新phi陣phi(1:i+1,1)=-y1( n:n+i-1+1);phi(1:i+1,2)=-y1( n-1: n+i-2+1);p hi(1:i+1,3)=u1( n:n+i-1+1);phi(1:i+1,4)

32、=u1( n-1: n+i-2+1);enddisp（'遞推廣義最小二乘法辨識結果：'）;a仁seta0(1),a2=seta0(2),b仁seta0(3),b2=seta0 (4)實驗結果：a1 =a1 =a1 =-0.7976-0.7872-0.8052a2 =a2 =a2 =-0.1941-0.2300-0.1857b1 =b1 =b1 =-0.00730.15770.1552b2 =b2 =b2 =0.36830.50900.4912Uy1.txtuy2.txtuy3.txt7.夏氏偏差修正法：%首先將TXT中的數據裝載在 Matlab中:clear;clc;Wbs =

33、 0;while Wbs<1disp('從下列文件中選擇所需加載文件： file name = inp ut('輸入所需打開文件名uy1.txt uy2.txt uy3.txt')endWbs,Message = fopen( file name,'。;:','s');打開文件;Data = fscan f(Wbs,'%g %g',2 in f);%生成數據矩陣Data = Data'%使用夏氏修正法，對2階系統進行參數辨識n = 2;L = len gth(Data);N = Ln;U = Data(:,1

34、);Y = Data(:,2);%首先使用最小二乘法，求系統參數向量初值glOL =-Y(2:L- 1),-Y(1:L-2),U(2:L-1),U(1:L-2);Zgl1 = Data(3:L,2);Sgl1 = glOL'*glOL;Sgl2=i nv(Sgl1);Sgl3=glO*Zgl1;Xsthita = Sgl2*Sgl3;%計算參數矩陣%得到初值后，以下使用夏氏修正法進行參數辨識:thitabO = 0;%設偏差項的偏差初值為0Fa = Sgl2*glOL'M = eye(N)-glOL*Sgl2*glOL'F = eye(N);if(F>=l0A-6

35、*eye(N)E1 = Zgl1-glOL*Xsthita;%計算殘差 Eomiga(2:N,1) = -E1(1:N-1);omiga(3:N,2) = -E1(1:N-2);D = omiga'*M*omiga;Fx = in v(D)*omiga'*M*Zgl1;thitab = Fa*omiga*Fx;Xsthita = Xsthita - thitab;F = thitab - thitabO;thitabO = thitab;enddisp('夏氏修正法')Xsa1 = Xsthita(1),Xsa2 = Xsthita(2),Xsb1 = Xsth

36、ita (3),Xsb2 = Xsthita(4)實驗結果:Uy1.txtuy2.txtuy3.txtXsa1 =-1.4789Xsa2 =0.8659Xsb1 =0.4259Xsb2 =0.5553Xsa1 =-1.4680Xsa2 =0.8548Xsb1 =0.7392Xsb2 =0.5712Xsa1 =-1.4651Xsa2 =0.8603Xsb1 =0.7423Xsb2 =0.55628.夏氏改良法：%首先將TXT中的數據裝載在 Matlab中:clear;clc;Wbs = 0;while Wbs<1disp('從下列文件中選擇所需加載文件： file name = i

37、np ut('輸入所需打開文件名uy1.txt uy2.txt uy3.txt')endWbs,Message = fopen( file name,'。;:','s');打開文件;Data = fscan f(Wbs,'%g %g',2 in f);%生成數據矩陣Data = Data'%最小二乘法，取b0=0n = 2;%根據題目，本系統為2階系統L = len gth(Data);%辨識所需數據的總長度N = L-n;% 構造測量矩陣的數據長度FIA = zeros(N,2* n);%構造測量矩陣%測量矩陣賦值for

38、 i = 1:Nfor l = 1:n*2if(l> n)FIA(i,l) = Data( n+2+i-l,1);endendelseendFIA(i,l) = -Data( n+i-1,2);Y = Data( n+1: n+N,2);%輸出數據矩陣thita = in v(FIA'*FIA)*FIA'*Y;%計算參數矩陣thitaLS = thita;%最小二乘估值m = in v(FIA'*FIA)*FIA'for i = 1:1:1000% 構造 OmegaErr = Y-FIA*thitaLS;% 計算殘差 Eomiga(2:N,1) = -Er

39、r(1:N-1);omiga(3:N,2) = -Err(1:N-2);F = in v(omiga'*omiga)*omiga'*Err;thita = thitaLS-m*omiga*F;enddisp('使用夏氏改良法辨識結果為:')Xga1 = thita(1),Xga2 = thita(2),Xgb1 = thita(3),Xgb2 = thita(4)實驗結果:Uy1.txtuy2.txtuy3.txtXga1 =Xga1 =Xga1 =-1.4779-1.4674-1.4651Xga2 =Xga2 =Xga2 =0.86540.85440.8603

40、Xgb1 =0.4258Xgb1 =0.7390Xgb1 =0.7424Xgb2 =0.5557Xgb2 =0.5716Xgb2 =0.5562實用標準文案9.夏氏遞推法：%首先將TXT中的數據裝載在 Matlab中: clear;clc;Wbs = 0;while Wbs<1enddisp('從下列文件中選擇所需加載文件:file name = inp ut('輸入所需打開文件名Wbs,Message = fopen( file name,'。;Data = fscan f(Wbs,'%g %g',2 in f);%Data = Data'

41、;%這里選用450組數據進行遞推迭代n = 2;N = 450-n;U = Data(:,1);Y = Data(:,2);%初始賦值如下:beta = 0 0 0 0 0 0' P0= 1010*eye(6,6);E1 = 0;E2 = 0;%循環迭代for i = 3:NE1uy1.txt uy2.txt uy3.txt')：','s');打開文件;生成數據矩陣Y(n+i)-beta(1)*Y( n+i-1)-beta (2)*Y( n+i-2)-beta(3)*Y( n+i-1)-beta (4) *U( n+i-2);E2Y(n+i-1)-bet

42、a(1)*Y( n+i-2)-beta (2) *Y( n+i-3)-beta(3)*U( n+i-2)-beta (4) *U( n+i-3);FIA = -Y( n+i) -Y(n+i-1) U(n+i) U(n+i-1) -E1 -E2'K = P 0*FIA*i nv(1+FIA'* P0 *FIA);P1 = P0-K*FIA'* P0;P0 = P1;betal = beta+K*(Y( n+i+1)-FIA*beta);beta = betal;enddisp('取450組數據進行遞推夏氏辨識的結果為:');Dxa1=beta(1),Dxa

43、2=beta (2),Dxb1=beta(3),Dxb2=beta (4)實驗結果:Dxa1 =-1.4067Dxa1 =-1.5289Dxa1 =-1.4919Uy1.txtuy2.txtuy3.txtDxa2 =Dxa2 =Dxa2 =0.80490.94510.8810Dxb1 =精彩文檔0.4300Dxb1 =Dxb1 =Dxb2 =0.56500.74930.7385Dxb2 =Dxb2 =0.54320.5410實用標準文案精彩文檔10.增廣矩陣法：%首先將TXT中的數據裝載在 Matlab中:clear;clc;Wbs = 0;while Wbs<1disp('從下

44、列文件中選擇所需加載文件： file name = inp ut('輸入所需打開文件名uy1.txt uy2.txt uy3.txt')endWbs,Message = fopen( file name,'。;：','s');打開文件;Data = fscan f(Wbs,'%g %g',2 in f);%生成數據矩陣Data = Data'%使用增廣矩陣法進行 2階系統的參數辨識n = 2;L = len gth(Data);U = Data(1:L,1);Y = Data(1:L,2);for i = 1:3*n %首

45、先給辨識參數向量賦初值，這里取thita0為0thitaO(i,1) = 0;endP0 = 1010*eye(3* n,3* n);%然后設定初始狀態 PO為充分大的對角陣ERRO = 0.00000000001;%thita = thita0,zeros(3* n,L);%zs = zeros(L,1);% 構造初始噪聲矩陣收斂情況指標，相對誤差 ERR0設定參數向量的初值及大小hn = 0,0,0,0,0,0'for j = n+1:1:Lzs(j-1) = Y(j-1)-h rTthita。;hn = -Y(j-1),-Y(j-2),U(j-1),U(j-2),zs(j-1),

46、zs(j-2)'計算修正矩陣KK = P 0*h n*i nv(h rf*P 0*h n+1);%thita0 = thita0+K*Y(j)-h n'*thita0;P0 = P 0-K*h rT P0;%計算狀態矩陣 P(N)thita(:,j) = thita0;E = abs(thita(:,j-1)-thita0)./thita0);if E<=ERR0 break;%根據設定的誤差限，判斷何時跳出endendZgma1=thita0(1),Zgma2=thita0(2),Zgmb1=thita0(3),Zgmb2=thita0(4),Zgmc1=thita0(

47、5),Zgmc2=thita 0(6)Uy1.txtuy2.txtuy3.txtZgmai = -1.4683Zgma2 =0.8626Zgmbi =0.4246Zgmai = -1.4629Zgma2 =0.8520Zgmbi =0.7396Zgma1 = -1.4573Zgma2 =0.8543Zgmb1 =0.7427Zgmb2 =0.5540Zgmb2 =0.5751Zgmb2 =0.5563Zgmcl =-0.9205Zgmcl =0.0282Zgmcl =0.5183Zgmc2 =0.0436Zgmc2 =-0.2689Zgmc2 =0.280611.極大似然法：%首先將TXT中的

48、數據裝載在 Matlab中:clear;clc;Wbs = 0;while Wbs<1disp('從下列文件中選擇所需加載文件： file name = inp ut('輸入所需打開文件名uy1.txt uy2.txt uy3.txt')endWbs,Message = fopen( file name,'。;:','s');打開文件;Data = fscan f(Wbs,'%g %g',2 in f);%生成數據矩陣Data = Data'n = 2;L = len gth(Data);N = Ln;U =

49、 Data(:,1);Y = Data(:,2);mthita = zeros(3* n,1);%給被辨識參數矩陣賦初始值%用最小二乘法求初始的被辨識參數矩陣OL = -Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2);Z1 = Data(3:L,2);thita = in v(OL'*OL)*OL'*Z1;mthita(1:2* n,1) = thita;%得到出之后，使用牛頓一拉布森極大似然法編程E = zeros(L,1);J = 0;sgma0 = 1;Ea1 = zeros(L,1);Ea2 = zeros(L,1);Eb1 = zeros(L

50、,1);Eb2 = zeros(L,1);Ec1 = zeros(L,1);Ec2 = zeros(L,1);%給J的梯度和Hassian矩陣賦初始值Ethita = zeros(3* n,1);fd = zeros(3* n,1);Hassia n = zeros(3* n,3* n);for p = 1:500Ja1mthita(1);Ja2mthita(2);Jb1mthita(3);Jb2mthita(4);Jc1mthita(5);Jc2 = mthita(6);for i = n+1:Ly2(i) = -Ja1*Y(i-1)-Ja2*Y(i-2)+Jb1*U(i-1)+Jb2*U(i-2)+Jc1*E(i-1)+Jc2*E(i-2);E(i) = Y(i)-y2(i);J = 0.5*(E(i)A2);sgma = (1/N)*(E(i)2);endfor j = n+1:1:LEa1(j)Y(j-1)-Jc1,Jc2*Ea1(j-1),Ea1(j-2)'Ea2(j)Y(j-2)-Jc1,Jc2*Ea2(j-1),Ea2(j-2)'Eb1(j)-U(j