系統(tǒng)辨識(shí)實(shí)驗(yàn)報(bào)告 (2)

1、.自動(dòng)化09-3 宋佳瑛 09051304系統(tǒng)辨識(shí)實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)一：系統(tǒng)辨識(shí)的經(jīng)典方法實(shí)驗(yàn)?zāi)康模赫莆障到y(tǒng)的數(shù)學(xué)模型與系統(tǒng)的輸入，輸出信號(hào)之間的關(guān)系，掌握經(jīng)典辨識(shí)的實(shí)驗(yàn)測(cè)試方法和數(shù)據(jù)處理方法。熟悉matlab/Simulink 環(huán)境。實(shí)驗(yàn)內(nèi)容：1.用系統(tǒng)階躍響應(yīng)法測(cè)試給定系統(tǒng)的數(shù)學(xué)模型。在系統(tǒng)沒有噪聲干擾的條件下通過測(cè)試系統(tǒng)的階躍響應(yīng)獲得系統(tǒng)的一階加純滯后或二階加純滯后模型，對(duì)模型進(jìn)行驗(yàn)證。2.在被辨識(shí)的系統(tǒng)加入噪聲干擾，重復(fù)上述1的實(shí)驗(yàn)過程。1. 沒有噪聲搭建對(duì)象 測(cè)試對(duì)象流程圖實(shí)驗(yàn)結(jié)果為：2、加入噪聲干擾搭建對(duì)象實(shí)驗(yàn)結(jié)果：加入噪聲干擾之后水箱輸出不平穩(wěn)，有波動(dòng)。實(shí)驗(yàn)二：相關(guān)分析法搭建對(duì)象：處理

2、程序：for i=1:15 m(i,:)=UY(32-i:46-i,1);endy=UY(31:45,2);gg=ones(15)+eye(15);g=1/(25*16*2)*gg*m*y;plot(g);hold on;stem(g);實(shí)驗(yàn)結(jié)果：相關(guān)分析法最小二乘法建模：二、三次實(shí)驗(yàn)本次實(shí)驗(yàn)要完成的內(nèi)容：1參照index2，設(shè)計(jì)對(duì)象，從workspace空間獲取數(shù)據(jù)，取二階，三階對(duì)象實(shí)現(xiàn)最小二乘法的一次完成算法和最小二乘法的遞推算法（LS and RLS）;2對(duì)設(shè)計(jì)好的對(duì)象，在時(shí)間為200-300之間，設(shè)計(jì)一個(gè)階躍擾動(dòng)，用最小二乘法和帶遺忘因子的最小二乘法實(shí)現(xiàn)，對(duì)這兩種算法的特點(diǎn)進(jìn)行說明；實(shí)

3、驗(yàn)內(nèi)容結(jié)果與程序代碼：利用LS和RLS得到的二階，三階參數(shù)算法階次A1A2A3B0B1B2B3LS二階-0.78420.1373-0.00360.56680.3157RLS二階-0.78240.1373-0.00360.56680.3157LS三階-0.4381-0.12280.0407-0.00780.56520.51060.1160RLS三階-0.4381-0.12280.0407-0.00780.56520.51060.1160以下給出RLS中的參數(shù)估計(jì)過程曲線和誤差曲線程序清單：LS（二階）：M=UY(:,1);z=UY(:,2);H=zeros(199,5);for i=1:199H

4、(i,1)=-z(i+1);H(i,2)=-z(i);H(i,3)=M(i+2);H(i,4)=M(i+1);H(i,5)=M(i);endEstimate=inv(H*H)*H*(z(3:201)RLS（二階）：clcM=UY(:,1);z=UY(:,2); P=100*eye(5); %估計(jì)方差Pstore=zeros(5,200);Pstore(:,1)=P(1,1),P(2,2),P(3,3),P(4,4),P(5,5);Theta=zeros(5,200); %參數(shù)的估計(jì)值，存放中間過程估值Theta(:,1)=0;0;0;0;0;K=10;10;10;10;10;10;10;for

5、 i=3:201h=-z(i-1);-z(i-2);M(i);M(i-1);M(i-2);K=P*h*inv(h*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h*Theta(:,i-2);P=(eye(5)-K*h)*P;Pstore(:,i-1)=P(1,1),P(2,2),P(3,3),P(4,4),P(5,5);endi=1:200;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:)title(待估參數(shù)過渡過程)figure(2)plot(i,P

6、store(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:)title(估計(jì)方差變化過程)同理可以寫出三階的LS以及RLS算法，此處略去。2.在t=250出加入一個(gè)階躍擾動(dòng)，令擾動(dòng)值為0.05利用RLS和帶遺忘因子的RLS計(jì)算結(jié)果如下表。算法遺忘因子A1A2A3B0B1B2B3RLS1-0.4791-0.09170.0349-0.00680.56590.48820.1033RLS0.98-0.5498-0.03900.0270-0.00690.56280.45150.0777RLS的參數(shù)變化過程曲線和誤差曲線如下：帶0.9

7、8遺忘因子的參數(shù)過度曲線和誤差曲線根據(jù)上面兩種方法所得到的誤差曲線和參數(shù)過渡過程曲線，我們可以看出來利用最小二乘法得到的參數(shù)最終趨于穩(wěn)定，為利用帶遺忘因子的最小二乘算法，曲線參數(shù)最終還是有小幅度震蕩。由此可以看出兩種算法的一些特點(diǎn)與區(qū)別。最小二乘法：遞推算法沒獲得一次新的觀測(cè)數(shù)據(jù)就修正一次參數(shù)估計(jì)值，隨著時(shí)間的推移，便能獲得滿意的辨識(shí)結(jié)果。帶遺忘因子的最小二乘法。其本質(zhì)還是最小二乘法，只不過加強(qiáng)新的數(shù)據(jù)提供的信息量，逐漸削弱老的數(shù)據(jù)，防止數(shù)據(jù)飽和。2. 參照index3，設(shè)計(jì)符合GLS和ELS的對(duì)象模型，改寫參照程序，實(shí)現(xiàn)相應(yīng)的算法。GLS 與ELS的實(shí)現(xiàn)參照index3，設(shè)計(jì)符合GLS和EL

8、S的對(duì)象模型，改寫參照程序，實(shí)現(xiàn)相應(yīng)的算法。利用GLS算法求出參數(shù)如下：a1a2b1b2-0.74110.08110.23640.1040即y(k)=0.7411y(k1)-0.0811y(k-2)+0.2364u(k-1)+0.1040u(k-2);利用ELS求解，得到如下的結(jié)果參數(shù)a1a2b1b2d1d2ELS-0.74240.08190.23460.103900繪制出ELS算法下參數(shù)的變化曲線以及誤差曲線如下：實(shí)驗(yàn)程序：GLS%廣義最小二乘法辨識(shí)%clcn=2;N=799;% y=y1;% u=u1;y=UY(:,2);u=UY(:,1);% y=y3;% u=u3;Y=y(n+1:n+

9、N);phi1=-y(n:n+N-1);phi2=-y(n-1:n+N-2);phi3=u(n:n+N-1);phi4=u(n-1:n+N-2);phi=phi1 phi2 phi3 phi4;theta=inv(phi*phi)*phi*Y %以上為最小二乘計(jì)算Theta估計(jì)值%f0=10;10;while 1e(n+1:n+N,1)=y(n+1:n+N)-phi*theta;omega(:,1)=-e(n:n+N-1);omega(:,2)=-e(n-1:n+N-2);f=inv(omega*omega)*omega*e(n+1:n+N,1);if (f-f0)*(f-f0)0.0001

10、breakendf0=f; %以上為最小二乘計(jì)算f估計(jì)值%for i=n+1:n+Nyy(i,1)=y(i),y(i-1),y(i-2)*1,f(1),f(2);uu(i,1)=u(i),u(i-1),u(i-2)*1,f(1),f(2);endyy(1,1)=y(1,1);uu(2,1)=u(2,1);%以上為數(shù)據(jù)濾波%phi(:,1)=-yy(n:n+N-1);phi(:,2)=-yy(n-1:n+N-2);phi(:,3)=uu(n:n+N-1);phi(:,4)=uu(n-1:n+N-2);theta=inv(phi*phi)*phi*yy(n+1:n+N);endfthetaELS%

11、增廣最小二乘的遞推算法%=產(chǎn)生均值為0，方差為1 的高斯白噪聲=v(1)=0;v(2)=0;z=UY(:,2);M=UY(:,1);%遞推求解P=100*eye(6); %估計(jì)方差Pstore=zeros(6,800);Pstore(:,1)=P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6);Theta=zeros(6,401); %參數(shù)的估計(jì)值，存放中間過程估值Theta(:,1)=3;3;3;3;3;3;K=10;10;10;10;10;10;for i=3:801h=-z(i-1);-z(i-2);M(i-1);M(i-2);v(i-1);v(i-2);

12、K=P*h*inv(h*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h*Theta(:,i-2);v(i)=z(i)-h*Theta(:,i-2);P=(eye(6)-K*h)*P;Pstore(:,i-1)=P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6);end%=disp(參數(shù)a1、a2、b1、b2、d1、d2 估計(jì)結(jié)果：)Theta(:,800)i=1:800;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:)

13、,i,Theta(6,:)title(待估參數(shù)過渡過程)figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:),i,Pstore(6,:)title(估計(jì)方差變化過程)大作業(yè)程序第一題x=1.01 2.03 3.02 4.01 5 6.02 7.03 8.04 9.03 10;y=9.6 4.1 1.3 0.4 0.05 0.1 0.7 1.7 3.8 9.1;ma,na=size(x);sumx=0;sumy=0;sumxy=0;sumx2=0;sumx2y=0;sumx3=0;s

14、umx4=0;for k=1:na sumx=sumx+x(k); sumy=sumy+y(k) sumxy=sumxy+x(k)*y(k); sumx2=sumx2+x(k)*x(k); sumx2y=sumx2y+x(k)*x(k)*y(k); sumx3=sumx3+x(k)3; sumx4=sumx4+x(k)4;endA=na,sumx,sumx2;sumx,sumx2,sumx3;sumx2,sumx3,sumx4;B=sumy;sumxy;sumx2y;C=AB；第二題aa=5;NNPP=7;ts=2;RR=ones(7)+eye(7);UU=UY(15:21,1);UY(14:

15、20,1);UY(13:19,1);UY(12:18,1);UY(11:17,1); UY(10:16,1);UY(9:15,1);YY=UY(8:14,2);GG=(RR*UU*YY+4.4474)/aa*aa*(NNPP+1)*ts; plot(0:2:13,GG);hold onstem(0:2:13,GG,filled);grid on;title(脈沖響應(yīng)序列)求系統(tǒng)的脈沖傳遞函數(shù) g= -0.0036 0.2922 0.3683 0.2673 0.1705 0.0820 0.0489; a=;for k=1:1:6a=a;g(k),g(k+1);endG=a;g(7),g(1);G

16、1=g(3:7,1);g(1);g(2);A=G(-G1);a1=A(2)a2=A(1)D=1 0;a1 1;C=g(1);g(2);B=D*C;b1=B(1)b2=B(2)zuixiao 相關(guān)二步法z=UY(:,2);M=UY(:,1);P=100*eye(4); Theta=zeros(4,196); Theta(:,1)=3;3;3;3;Theta(:,2)=3;3;3;3;Theta(:,3)=3;3;3;3;Theta(:,4)=3;3;3;3;K=10;10;10;10;for i=5:196h=-z(i-1);-z(i-2);M(i-1);M(i-2);hstar=M(i-1);

17、M(i-2);M(i-3);M(i-4); K=P*hstar*inv(h*P*hstar+1);Theta(:,i)=Theta(:,i-1)+K*(z(i)-h*Theta(:,i-1);P=(eye(4)-K*h)*P;endi=1:196;plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)title(相關(guān)二步法參數(shù)過度過程);grid on最小二乘遞推算法相關(guān)參數(shù)clc;lamt=1;z=UY(:,2);u=UY(:,1);c0=0.001 0.001 0.001 0.001;p0=104*eye(4,4); c=c0,ze

18、ros(4,198); for k=3:200 h1=-z(k-1),-z(k-2),u(k-1),u(k-2); x=h1*p0*h1+1*lamt; x1=inv(x); k1=p0*h1*x1; d1=z(k)-h1*c0; c1=c0+k1*d1; p1=1/lamt*(eye(4)-k1*h1)*p0; c(:,k-1)=c1; c0=c1;p0=p1; enda1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); i=1:199;plot(i,a1,r,i,a2,b,i,b1,y,i,b2,black) ;title(遞推最小二乘法參數(shù)辨識(shí) );gr

19、id on帶遺忘因子的最小二乘法將上述程序中的lamt=0.9；辨識(shí)系統(tǒng)階次方法一、利用殘差最小辨識(shí)Data=UY;L = length(Data);Jn = zeros(1,10);t = zeros(1,10);rm = zeros(10,10);for n=1:1:10; N = L-n; FIA = zeros(N,2*n); du = zeros(n,1); dy = zeros(n+1,1); r1 = 0;r0 = 0;for i = 1:N for l = 1:n*2 if(l0) FIA(i,l) = Data(n+i-l+2,1); end endendY = Data(n

20、+1:n+N,2);thita = inv(FIA*FIA)*FIA*Y; Jn0 = 0; for k = n+1:n+N for j = 1:n du(j) = Data(k-j,1); dy(j) = Data(k+1-j,2); end dy(n+1) = Data(1,2); E1(k) = 1,thita(1:n)*dy-thita(n+1:2*n)*du; Jn1 = Jn0+E1(k)2; t(n) = abs(Jn0-Jn1)/Jn1*(N-2*n-2)/2); Jn0 = Jn1; end Jn(n) = Jn0; for m = 0:1:9 for m2 = n+1:1:

21、L-m r1 = r1+E1(m2)*E1(m2+m)/(L-m-n); r0 = r0+E1(m2)2/(L-m-n); end rm(n,m+1) = r1/r0; endendsubplot(2,1,1);plot(1:10,Jn,r);title(殘差平方和-jn曲線圖);subplot(2,1,2);plot(1:1:10,t,g);title(F檢驗(yàn)結(jié)果圖);figure(2);plot(0:9,rm(1,:),g),hold on;plot(0:9,rm(2,:),b),hold on;plot(0:9,rm(3,:),r);title(殘差白色丁階結(jié)果圖);方法二、AIC方法訂

22、階次Data = UY;L = length(Data);U = Data(:,1);Y = Data(:,2);for n = 1:1:5 N = 201-n yd(1:N,1) = Y(n+1:n+N); for i=1:N for l=1:1:2*n if(l=n) FIA(i,l) = -Y(n+i-l); else FIA(i,l) = U(2*n+i-l); end end end thita = inv(FIA*FIA)*FIA*yd; omiga = (yd-FIA*thita)*(yd-FIA*thita)/N; AIC(n) = N*log(omiga)+4*n;endpl

23、ot(AIC,-);title(AI訂階法);xlabel(n);ylabel(AIC);grid on第三題廣義最小二乘法：n=2;N=799;u=UY(:,1);y=UY(:,2);Y=y(n+1:n+N);phi1=-y(n:n+N-1);phi2=-y(n-1:n+N-2);phi3=u(n+1:n+N);phi4=u(n:n+N-1);phi5=u(n-1:n+N-2)phi=phi1 phi2 phi3 phi4 phi5;theta=inv(phi*phi)*phi*Y ;f0=10;10;while 1e(n+1:n+N,1)=y(n+1:n+N)-phi*theta;omeg

24、a(:,1)=-e(n:n+N-1);omega(:,2)=-e(n-1:n+N-2);f=inv(omega*omega)*omega*e(n+1:n+N,1);if (f-f0)*(f-f0)=10-6*eye(N) E1 = Zgl1-glOL*Xsthita;%計(jì)算殘差E omiga(2:N,1) = -E1(1:N-1); omiga(3:N,2) = -E1(1:N-2); D = omiga*M*omiga; Fx = inv(D)*omiga*M*Zgl1; thitab = Fa*omiga*Fx; Xsthita = Xsthita - thitab; F = thitab

25、 - thitab0; thitab0 = thitab;enddisp(夏氏修正法)Xsa1 = Xsthita(1),Xsa2 = Xsthita(2),Xsb1 = Xsthita(3),Xsb2 = Xsthita(4)夏氏改良法Data=UY;%最小二乘法，取b0=0n = 2;%根據(jù)題目，本系統(tǒng)為2階系統(tǒng)L = length(Data);%辨識(shí)所需數(shù)據(jù)的總長度N = L-n;%構(gòu)造測(cè)量矩陣的數(shù)據(jù)長度FIA = zeros(N,2*n);%構(gòu)造測(cè)量矩陣for i = 1:N %測(cè)量矩陣賦值 for l = 1:n*2 if(ln) FIA(i,l) = Data(n+2+i-l,1)

26、; else FIA(i,l) = -Data(n+i-l,2); end endendY = Data(n+1:n+N,2);%輸出數(shù)據(jù)矩陣thita = inv(FIA*FIA)*FIA*Y;%計(jì)算參數(shù)矩陣thitaLS = thita;%最小二乘估值m = inv(FIA*FIA)*FIA;for i = 1:1:1000%構(gòu)造Omega Err = Y-FIA*thitaLS;%計(jì)算殘差E omiga(2:N,1) = -Err(1:N-1); omiga(3:N,2) = -Err(1:N-2); F = inv(omiga*omiga)*omiga*Err; thita = thi

27、taLS-m*omiga*F;enddisp(使用夏氏改良法辨識(shí)結(jié)果為：)Xga1 = thita(1),Xga2 = thita(2),Xgb1 = thita(3),Xgb2 = thita(4)第五題用一般遞推方法辨識(shí) z(k)=-1.5*z(k-1)+4.5*u(k-1)+2.8*u(k-2)A=6;x0=1;M=255;for k=1:10000 x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=(v1-0.5)*2; x0=x1; v0=v1;endnum=v;k1=k;y1=1;y2=1;y3=1;y4=0;for i=1:255; x1=xor(y3,y4); x2=y1; x3=y2; x4=y3; y(i)=y4; if y(i)0.5,u(i)=-1; else u(i)=1; end y1=x1;y2=x2;y3=x3;y4=x4; lamt=1;z(2)=0;z(1)=0;for k=3:15; z(k)=-1.5*z(k-1)+4.5*u(k-1)+2.8*u(k-2)+0*num(1,k); end% c0=0.001 0.001 0.001;p0=103*eye(3,3);c=c0,zeros(3,14);for k=3:15; h1=-z(k-1),u(k-1),u(k-2); x=h1*p0*h1+