系統(tǒng)辨識實驗

實驗一：系統(tǒng)辨識的經(jīng)典方法實驗?zāi)康模赫莆障到y(tǒng)的數(shù)學(xué)模型與系統(tǒng)的輸入、輸出信號之間的關(guān)系，掌握經(jīng)典辨識的實驗測試方法和數(shù)據(jù)處理方法。熟悉MATLAB/Simulink環(huán)境。實驗內(nèi)容：用階躍響應(yīng)發(fā)測試給定系統(tǒng)的數(shù)學(xué)模型在系統(tǒng)沒有噪聲干擾的條件下通過測試系統(tǒng)的階躍響應(yīng)獲得系統(tǒng)的一階加純滯后或者二階加純滯后模型，對模型進(jìn)行驗證。在被辨識系統(tǒng)中加入噪聲干擾，重復(fù)上述1的實驗過程。實驗方法：在MATLAB環(huán)境下用simulink構(gòu)造測試環(huán)境，被測試的模型為水槽液位控制對象。利用非線性水槽模型（tank）可以搭建單水槽系統(tǒng)的模型，也可以搭建多水槽系統(tǒng)的模型，多水槽模型可以是高低放置，也可以并排放置。圖1.1二階水槽測試流程圖（1）一階環(huán)節(jié)傳函

地不網(wǎng)DW3』尊wuilllh：g■><-le_y(8)_27.1K—2.71|QtAfaT-￡?■|ftL.Lr.,.也,曰&3|￡Ji，ryt|呵。心-—EiJF"JD辜事孤才口3w>*zAuA地不網(wǎng)DW3』尊wuilllh：g■><-le_y(8)_27.1K—2.71|QtAfaT-￡?■|ftL.Lr.,.也,曰&3|￡Ji，ryt|呵。心-—EiJF"JD辜事孤才口3w>*z利用兩點法運(yùn)用公式得到t=2匕-1=2X1.1-2.3—-0.1T=2(,—<)-2x(2.3-1.1)=2.4所以確定Ke-Ts2.71xeo.isG所以確定Ke-Ts2.71xeo.isG(s)=Ts+12.4s+1(2)二階環(huán)節(jié)傳函確定K=Ay=衛(wèi)-變-2.66AuAu10找到^(<)-0.4,v(t2)-0.8的點，對應(yīng)得到匕-4.0，七-9」；運(yùn)用公式fT+T=(t+1).216=(4.0+9.1)/2.16=6.06《1221TT/(T+T)2n1.74tIt-0.55=1.74x4.0/9.1-0.55=0.21'12121!2求的所以確定T=2.1211T=0.909K_2.66(〈s+1)(7^+1)—(2.121s+1)(0.909s+1)思考題實驗二相關(guān)分析法實驗?zāi)康模赫莆障嚓P(guān)分析法測試系統(tǒng)數(shù)學(xué)模型的過程和方法。掌握應(yīng)用移位計數(shù)器設(shè)計M序列信號發(fā)生器的方法。實驗內(nèi)容：設(shè)計并實現(xiàn)PRBS偽隨機(jī)序列信號（M序列）發(fā)生器；應(yīng)用相關(guān)分析法測試給定系統(tǒng)的數(shù)學(xué)模型。使用偽隨機(jī)序列信號測試系統(tǒng)在有噪聲情況下的輸出，用相關(guān)分析法辨識系統(tǒng)的數(shù)學(xué)模型。模型驗證。實驗方法：測試對象可以與實驗一相同。應(yīng)用MATLAB/Simulink設(shè)計PRBS偽隨機(jī)信號（M序列）發(fā)生器，可以用JK觸發(fā)器構(gòu)造的移位計數(shù)器實現(xiàn)，也可以用程序?qū)崿F(xiàn)。對于n級M序列信號發(fā)生器反饋半加器的兩個輸入端，一個取自末級，一個曲子k級（見下表），輸出反饋到第一個觸發(fā)器的輸入端。n2345678910k11or21or32or31or51or3or4or62or3or44or53or7根據(jù)被辨識模型的特點，選取合適的時鐘周期、移位計數(shù)器的級數(shù)和信號幅值。M序列信號的設(shè)計測試對象

程序aa=5;NNPP=15;ts=2;RR=ones(15)+eye(15);fori=15:-1:1UU(16-i,:)=UY(16+i:30+i,1)';endYY=[UY(31:45,2)];GG=RR*UU*YY/[aa*aa*(NNPP+1)*ts];plot(0:2:29,GG)holdonstem(0:2:29,GG,'filled')結(jié)果：工作點的脈沖響應(yīng)序列

-0.050工作點的脈沖序列響應(yīng)0.40.350.30.250.20.150.10.0551015202530實驗三最小二乘法參數(shù)估計-0.050工作點的脈沖序列響應(yīng)0.40.350.30.250.20.150.10.0551015202530實驗?zāi)康模赫莆兆钚《朔ǖ幕驹砗退惴ㄒ约霸谙到y(tǒng)辨識中的應(yīng)用，編制系統(tǒng)辨識的數(shù)據(jù)處理、最小二乘法程序一次辨識程序和帶遺忘因子的遞推最小二乘法程序。實驗內(nèi)容：1.參照index2，設(shè)計對象，從workspace空間獲取數(shù)據(jù)，取二階，三階對象實現(xiàn)最小二乘法的一次完成算法和最小二乘法的遞推算法(LSandRLS);從workspace工作空間中獲取數(shù)數(shù)據(jù)，用一次完成算法實現(xiàn)最小二乘，LS程序如下：M=UY(:,1);z=UY(:,2);H=zeros(140,5);fori=1:140H(i,1)=-z(i+11);H(i,2)=-z(i+10);H(i,3)=M(i+12);H(i,4)=M(i+11);H(i,5)=M(i+10);endEstimate=inv(H'*H)*H'*(z(13:152))估計結(jié)果Estimate=-0.78540.1387-0.00200.56940.3136RLS程序：%最小二乘的遞推算法M=UY(:,1);z=UY(:,2);%遞推求解P=100*eye(4);%估計方差Pstore=zeros(4,140);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)]';Theta=zeros(4,140);%參數(shù)的估計值，存放中間過程估值Theta(:,1)=[3;3;3;3];%K=zeros(4,400);%增益矩陣K=[10;10;10;10];fori=3:142h=[-z(i-1);-z(i-2);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(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)]';endi=1:141;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:))title('待估參數(shù)過渡過程')figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:))title('估計方差變化過程')估計結(jié)果:

32.521.510.50-0.5-1-1.5待估參數(shù)過渡過程50100150估計方差變化過程10090807060504030201032.521.510.50-0.5-1-1.5待估參數(shù)過渡過程50100150估計方差變化過程100908070605040302010501001502.LS程序估計結(jié)果：Estimate=-1.22050.24930.05230.61730.1588帶遺忘因子的最小二乘法實現(xiàn)lamt=1;z=UY(:,2);u=UY(:,1);c0=[0.0010.0010.0010.001]';p0=10A4*eye(4,4);c=[c0,zeros(4,198)];fork=3:200h1=[-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;c(:,199)結(jié)果ans=-1.32370.35160.60770.11353.參照index3，設(shè)計符合GLS和ELS的對象模型，改寫參照程序，實現(xiàn)相應(yīng)的算法。GLS程序：v=randn(1,400);e=[];e(1)=v(1);e(2)=v(2);fori=3:400e(i)=0*e(i-1)+0*e(i-2)+v(i);endM=UY(:,1);z=UY(:,2);%變換后的觀測序列zf=[]；zf(1)=-1；zf(2)=0;fori=3:142zf(i)=z(i)-0*z(i-1)-0*z(i-2);end%變換后的輸入序列uf=[];uf(1)=M(1);uf(2)=M(2);fori=3:142uf(i)=M(i)-0*M(i-1)-0*M(i-2);end%賦初值P=100*eye(4);%估計方差Pstore=zeros(4,142);Pstore(:,2)=[P(1,1),P(2,2),P(3,3),P(4,4)]';Theta=zeros(4,140);%參數(shù)的估計值，存放中間過程估值Theta(:,2)=[3;3;3;3];K=[10;10;10;10];%增益PE=10*eye(2);ThetaE=zeros(2,140);ThetaE(:,2)=[0.5;0.3];KE=[10;10];fori=3:142h=[-zf(i-1);-zf(i-2);uf(i-1);uf(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i)=Theta(:,i-1)+K*(z(i)-h'*Theta(:,i-1));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)]';he=[-e(i-1);-e(i-2)];KE=PE*he*inv(1+he'*PE*he);ThetaE(:,i)=ThetaE(:,i-1)+KE*(e(i)-he'*ThetaE(:,i-1));PE=(eye(2)-KE*he')*PE;end%=====================輸出結(jié)果及作圖=========================disp('參數(shù)a1a2b1b2的估計結(jié)果：’)Theta(:,140)disp('噪聲傳遞系數(shù)c1c2的估計結(jié)果：’)

ThetaE(:,141)i=1:142;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:))title('待估參數(shù)過渡過程')figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:))title('估計方差變化過程')figure(3)plot(i,ThetaE(1,:),i,ThetaE(2,:));辨識結(jié)果：a1a2b1b2ans=-0.74050.08230.23730.1052c1c2ans=-0.05750.049032.521.510.50-0.5待估參數(shù)過渡過程-1

-11.510.50-0.550100150ELS：待估參數(shù)過渡過程-11.510.50-0.550100150程序：%增廣最小二乘的遞推算法M=UY(:,1);z=UY(:,2);%遞推求解P=100*eye(6);%估計方差Pstore=zeros(6,140);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)];Theta=zeros(6,140);%參數(shù)的估計值，存放中間過程估值Theta(:,1)=[3;3;3;3;3;3];%K=zeros(4,400);%增益矩陣K=[10;10;10;10;10;10];fori=3:141h=[-z(i-1);-z(i-2);M(i-1);M(i-2);v(i-1);v(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(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估計結(jié)果：’)Theta(:,140)i=1:140;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:),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('估計方差變化過程')辨識結(jié)果：ans=-0.74120.08110.23640.10420.00020.0001543210-1-2待估參數(shù)過渡過程20406080100120140估計方差變化過程1009080706050403020100020406080100120140附錄：大作業(yè)程序：第一題：z=[1.012.033.024.0156.027.038.049.0310]';M=[9.64.11.30.40.050.10.71.73.89.1]';fori=1:10H(i,1)=z(i)*z(i);H(i,2)=z(i);H(i,3)=1;endEstimate=inv(H'*H)*H'*M第二題：求脈沖響應(yīng)序列：GG=[-0.00280.29290.38180.26830.15990.07550.03780.02610.02080.01820.01060.01310.00950.0096-0.0009];h=zeros(2,2);a=zeros(2,1);b=zeros(2,1);c=GG(1,3:4)';fori=1:-1:0h(2-i,:)=GG(1,2-i:1:3-i)';enda=inv(h)*cb=[1,0;a(1),1]*[GG(1)GG(2)]'RLS:M=UY(:,1);z=UY(:,2);P=100*eye(5);%估計方差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ù)的估計值，存放中間過程估值Theta(:,1)=[0.1;0.1;0.1;0.1;0.1];%K=zeros(4,400);%增益矩陣K=[10;10;10;10;10];fori=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,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:))title('估計方差變化過程')Theta(:,200)帶遺忘因子：lamt=0.95;z=UY(:,2);u=UY(:,1);c0=[0.0010.0010.0010.001]';p0=10"4*eye(4,4);c=[c0,zeros(4,198)];fork=3:200h1=[—z(kT),—z(k—2),u(kT),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');gridonc(:,199)AIC定階：%%%%%%%%%%%%%%%階次辨識(AKaike)%%%%%%%%%%%y1=UY(:,2);u1=UY(:,1);forn=1:5N=200—n;Y1=y1(n+1:n+N,1);fori=1:nphi1(1:N,i)=-y1(n+1-i:n+N-i,1);phi1(1:N,n+i)=u1(n+1-i:n+N-i,1);endtheta1=inv(phi1'*phi1)*phi1'*Y1;e1(n+1:n+N,1)=y1(n+1:n+N)-phi1*theta1;fori=1:npsai1(1:N,i)=-y1(n+1-i:n+N-i,1);psai1(1:N,n+i)=u1(n+1-i:n+N-i,1);psai1(1:N,2*n+i)=e1(n+1-i:n+N-i,1);endthetaa1=inv(psai1'*psai1)*psai1'*Y1;E1(n+1:n+N,1)=y1(n+1:n+N)-psai1*thetaa1;sigma1(n,1)=E1'*E1/N;m(n)=n;phi1=0;psai1=0;B1(n,1)=N*log(sigma1(n,1))+3*n;endplot(m,B1,'r-o');gridontitle('數(shù)據(jù)1AKaike階次辨識')；殘差定階：%%%%%%%%%%%%%%階次辨識(殘差最小)％%%%%%%%%%%u=UY(:,1);y=UY(:,2);forn=1:5%對各階次計算JnN=200-n;Y=y(n+1:n+N,1);forj=1:Nfork=1:nphi(j,k)=-y(n+j-k);phi(j,n+k)=u(n+j-k);endendtheta=inv(phi'*phi)*phi'*Y;e(n+1:n+N,1)=y(n+1:n+N)-phi*theta;J(n,1)=e'*e;%計算Jnphi=0;m(n)=n;endplot(m,J,，r-o，);title('數(shù)據(jù)1階次曲線圖，)；holdongridon第三題GLSv=randn(1,400);e=[];e(1)=v(1);e(2)=v(2);fori=3:400e(i)=0*e(i-1)+0*e(i-2)+v(i);endM=UY(:,1);z=UY(:,2);%變換后的觀測序列zf=[];zf(1)=-1;zf(2)=0;fori=3:142zf(i)=z(i)-0*z(i-1)-0*z(i-2);end%變換后的輸入序列uf=[];uf⑴=M(1);uf(2)=M(2);fori=3:142uf(i)=M(i)-0*M(i-1)-0*M(i-2);end%賦初值P=100*eye(4);%估計方差Pstore=zeros(4,142);Pstore(:,2)=[P(1,1),P(2,2),P(3,3),P(4,4)]';Theta=zeros(4,140);%參數(shù)的估計值，存放中間過程估值Theta(:,2)=[3;3;3;3];K=[10;10;10;10];%增益PE=10*eye(2);ThetaE=zeros(2,140);ThetaE(:,2)=[0.5;0.3];KE=[10;10];fori=3:142h=[-zf(i-1);-zf(i-2);uf(i-1);uf(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i)=Theta(:,i-1)+K*(z(i)-h'*Theta(:,i-1));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)]';he=[-e(i-1);-e(i-2)];KE=PE*he*inv(1+he'*PE*he);ThetaE(:,i)=ThetaE(:,i-1)+KE*(e(i)-he'*ThetaE(:,i-1));PE=(eye(2)-KE*he')*PE;end%======輸出結(jié)果及作圖disp('參數(shù)a1a2b1b2的估計結(jié)果：’)Theta(:,140)disp('噪聲傳遞系數(shù)c1c2的估計結(jié)果：’)ThetaE(:,141)i=1:142;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:))title('待估參數(shù)過渡過程')figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:))title('估計方差變化過程')figure(3)plot(i,ThetaE(1,:),i,ThetaE(2,:));ELS：M=UY(:,1);z=UY(:,2);%遞推求解P=100*eye(6);%估計方差Pstore=zeros(6,140);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)];Theta=zeros(6,140);%參數(shù)的估計值，存放中間過程估值Theta(:,1)=[3;3;3;3;3;3];%K=zeros(4,400);%增益矩陣K=[10;10;10;10;10;10];fori=3:141h=[-z(i-1);-z(i-2);M(i-1);M(i-2);v(i-1);v(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(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估計結(jié)果：’)Theta(:,140)i=1:140;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:),i,Theta(6,:))title('待估參數(shù)過渡過程')Iegend('a1','a2','b1','b2','d1','d2');figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:),i,Pstore(6,:))title('估計方差變化過程')第四題：夏氏法：Data=UY;n=2;L=length(Data);N=L-n;U=Data(:,1);Y=Data(:,2);OL=[-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2)];Z1=Data(3:L,2);S1=OL'*OL;S2=inv(S1);S3=OL'*Z1;Xsthita=S2*S3;E1=Z1-OL*Xsthita;E2(1,1)=0.001;E2(2:L-1,1)=E1(1:L-2);fn=[-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2),-E2(1:L-2),-E2(2:L-1)];n1=fn,*fn;thita0=inv(n1)*fn'*Z1;P0=10"9*eye(n*3,n*3);fori=n+1:1:Lform=1:1:n*3if(m<=2)fm(m,1)=-Data(i-m,2);elseif(m>2&&m<=4)fm(m,1)=Data(i-m+2,1);elsefm(m,1)=-E2(i-m+4,1);endendX=fm,*P0*fm+1;K1=P0*fm/X;D1=Data(i,2)-fm'*thita0;P1=P0-K1*K1'*[fm'*P0*fm+1];P0=P1;%X^IaIA，I^u，u360^thita1=thita0+K1*D1;thita0=thita1;endXsa1=Xsthita(1),Xsa2=Xsthita(2),Xsb1=Xsthita(3),Xsb2=Xsthita(4)夏氏改良：Data=UY;n=2;L=length(Data);N=L-n;FIA=zeros(N,2*n);fori=1:Nforl=1:n*2if(l>n)FIA(i,l)=Data(n+2+i-l,1);elseFIA(i,l)=-Data(n+i-l,2);endendend=Data(n+1:n+N,2);thita=inv(FIA'*FIA)*FIA'*Y;thitaLS=thita;m=inv(FIA'*FIA)*FIA';fori=1:1:1000%11O!OmegaErr=Y-FIA*thitaLS;omiga(2:N,1)=-Err(1:N-1);omiga(3:N,2)=-Err(1:N-2);F=inv(omiga'*omiga)*omiga'*Err;thita=thitaLS-m*omiga*F;enddisp('日FAiA^LAA%,??土建，屈1苗afa￡°')Xga1=thita(1),Xga2=thita(2),Xgb1=thita(3),Xgb2=thita(4)夏氏修正：Data=UY;n=2;L=