系统辨识实验报告.docx

《系统辨识实验报告.docx》由会员分享，可在线阅读，更多相关《系统辨识实验报告.docx（32页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、.word.自动化093 宋佳瑛 09051304系统辨识实睑报告实验一：系统辨识的经典方法实验目的：掌握系统的数学模型与系统的输入，输出信号之间 的关系，掌握经典辨识的实验测试方法和数据处理方法。熟悉 matlab/Simulink 环境。实验内容：1 .用系统阶跃响应法测试给定系统的数学模型。在系统没有噪声干扰的条件下通过测试系统的阶跃响应获得系 统的一阶加纯滞后或二阶加纯滞后模型，对模型进行验证。2 .在被辨识的系统加入噪声干扰,重复上述1的实验过程。1 .没有噪声 搭建对象带0.98遗忘因子的参数过度曲线和误差曲线根据上面两种方法所得到的误差曲线和参数过渡过程曲线，我 们可以看出来利用

2、最小二乘法得到的参数最终趋于稳定，为利 用带遗忘因子的最小二乘算法，曲线参数最终还是有小幅度震 荡。由此可以看出两种算法的一些特点与区别。最小二乘法：递推算法没获得一次新的观测数据就修正一次参 数估计值，随着时间的推移，便能获得满意的辨识结果。带遗忘因子的最小二乘法。其本质还是最小二乘法，只不过加 强新的数据提供的信息量，逐渐削弱老的数据，防止数据饱和。2 .参照index3 ,设计符合GLS和ELS的对象模型，改写参照 程序，实现相应的算法。.word.GLS与ELS的实现参照index3 ,设计符合GLS和ELS的对象模型，改写参照程序, 实现相应的算法。利用GLS算法求出参数如下:ala

3、2blb2-0.74110.08110.23640.1040即y(k)=0.7411y(kl)-0.0811y(k-2)+0.2364u(k-l)+0.1040u(k-2);利用ELS求解，得到如下的结果参数ala2blb2dld2ELS-0.7420.0810.2340.103004969绘制出ELS算法下参数的变化曲线以及误差曲线如下:实验程序：GLS.word.% 广义最识%cicn=2;N=799;% y=yl;% u = ul;y = UY(:,2);u = UY(:,l);% y=y3;% u = u3;Y=y(n+l:n+N);phil=-y(n:n + N-l);phi2 =

4、-y(n-l:n + N-2);phi3 = u(n:n + N-l);phi4=u(n-l:n + N-2);phi = phil phi2 phi3 phi4;theta=inv(phi*phi)*phi*Y%以上为最小二乘计 值%f0=10;10;while 1小二乘法辨小二乘法辨算Theta估计e(n+l:n+N,l)=y(n+l:n+N)-phi*theta;.word.omega(:/l)=-e(n:n+N-l);omega(:/2)=-e(n-l:n+N-2);f=inv(omega*omega)*omega,*e(n+l:n+N/l);if (f-fO)*(f-fO) 0.00

5、01breakendfO=f;%以上为最小二乘计算f估计值% for i=n + l:n + Nyy(i,D=y(i),y(i-i),y(i-2)*LfQ),f(2 比 uu(i,l)=u(i),u(il),u(i-2)*Lf 术2比 endyy(l,l)=y(l,l);uu(2/l) = u(2/l);%以上为数据滤波%phi(:,l)=-yy(n:n + N-l);phi(:/2)=-yy(n-l:n+N-2);phi(:/3) = uu(n:n + N-l);phi(:/4)=uu(n-l:n+N-2);theta=inv(phi*phi)*phi*yy(n+l:n + N);end t

6、heta.word.ELS%增广最小二乘的递推算法 %=产生均值为0 ,方差为1的高斯白噪声v(l)=0;v(2)=0;z= UY(:,2);M = UY(:,1);%递推求解P=100*eye(6); %估计方差Rstore=zeros(6z800);Pstore(:,l) = P(Ll),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6);Theta=zeros(6/401); %参数的估计值，存放中间过程估值Theta(:/1) = 3;3;3;3;3;3;K=10;10;10;10;10;10;for i=3:801h=-z(i-l);-z(i-2);M(i-l);M

7、(i-2);v(i-l);v(i-2);K=P*h*inv(h*P*h + l);Theta(:/i-l)=Theta(:/i-2) + K*(z(i)-h*Theta(:/i-2);v(i)=z(i)-h*Theta(:J-2);.word.P=(eye(6)-K*h)*P;Pstore(:,i-l) = P(l,l),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6); end%=dispC参数al、a2、bl. b2、dl、d2 估计结果：,) Theta(:,800)i = 1:800;figure plot(iJheta(l/:)zizTheta(2,:)zi,Th

8、eta(3/:)J/Theta(4,:)ziJheta( 5,:) /i/Theta(6/:)title。待估参数过渡过程)figure。)plot(i,Pstore(l,：),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Psto re(51:)J/Pstore(67:)titled估计方差变化过程).word.大作业程序第一题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;su

9、my=0;sumxy=0;sumx2=0;sumx2y=0;sumx3 = 0;sumx4=0;for k=l:nasumx=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);.word.sumx3=sumx3+x(k)A3;sumx4=sumx4+x(k)A4;endA=nazsumxzsumx2;sumxzsumx2/sumx3;sumx2zsumx3/su mx4;B=sumy;sumxy;sumx2y;C=AB ;第二题aa = 5;NNPP=

10、7;ts=2;RR=ones(7)+eye(7);UU=UY(15:21,l),;UY(14:20,l),;UY(13:19zl),;UY(12:18,l),;UYQ1:17,1),;UY(10:16,l);UY(9:15,l);YY=UY(8:14Z2);GG=(RR*U U*YY+4.4474)/aa*aa*(N N PP+l)*ts;plot(0:2:13,GG);hold onstem(0:2:13zGG；filled);grid on;title。脉冲响应序列。.word.求系统的脉冲传递函数0.2673g = -0.00360.29220.36830.17050.08200.048

11、9;a二口；for k=l:l:6a=a;g(k)fg(k+l);endG=a;g ,g(l);Gl=g(3:7fl);g(l);g(2);A=G(-G1);al=A(2)a2=A(l)D=l 0;al 1;C=g(l);g(2);B=D*C;bl=B(l)b2=B(2)zuixiao相关二步法z=UY(:,2);M = UY(;1);P=100*eye(4);.word.Theta=zeros(4,196);Theta(:,l) = 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;f

12、ori = 5:196h=-z(i-l);-z(i-2);M(i-l);M(i-2);hstar=M(i-l);M(i-2);M(i-3);M(i-4);K=P*hstar*inv(h*P*hstar+l);Theta(:,i)=Theta(:/i-l) + K*(z(i)-h,*Theta(:,i-l);P=(eye(4)-K*h)*P;endi = l:196;plot(i/Theta(l/:)/i/Theta(2,:)J/Theta(3,:),izTheta(4/:) title。相关二步法参数过度过程)；grid on最小二乘递推算法相关参数cic;lamt=l;z=UY(:,2);u

13、 = UY(:,l);.word.测试对象流程图实验结果为：搭建对象.word.c0=0.001 0.001 0.001 0.001;p0=10A4*eye(4z4);c=c0,zeros(4,198);for k=3:200hl=-z(k-l)rz(k-2),u(k-l),u(k-2),;x=hl*pO*hl+l*lamt;xl=inv(x);kl=p0*hl*xl;dl=z(k)-hl*c0;cl=c0+kl*dl;pl=l/lamt*(eye(4)-kl*hl)*p0;c(:,k-l)=cl;c0=cl;p0=pl;endal=c(l,:); a2=c(2,:); bl=c(3z:);

14、b2=c(4,:);i = l:199;plot(i,aLr,i,a2b,i,bl,y,i,b2,black);title。递推最小二乘法参数辨识,)；grid on带遗忘因子的最小二乘法将上述程序中的gmt=0.9 ;.word.辨识系统阶次方法一、利用残差最小辨识Data = UY;L = length(Data);Jn = zeros(lz10);t = zeros(lz10);rm = zeros(10z10);for n = l:l:10;N = L-n;FIA = zeros(N,2*n);du = zeros(n,l);dy = zeros(n + lzl);rl = 0;r0

15、= 0;for i = 1:Nfor I = l:n*2 if(l0)FIA(iJ) = Data(n + i-l+2zl); end end end.word.Y = Data(n+l:n+N,2);thita = inv(FIA*FIA)*FIA,*Y;JnO = 0;for k = n + l:n + Nfor j = l:ndu(j) = Data(k-jJ);dy(j) = Data(k+l-jz2);enddy(n + l) = Data(lz2);El(k) = l/thita(l:n),*dy-thita(n + l:2*n)*du;Jnl = JnO+El(k)A2;t(n)

16、 = abs(JnO-Jnl)/Jnl*(N-2*n-2)/2);JnO = Jnl;endJn(n) = JnO;for m = 0:1:9for m2 = n + l:l:L-mrl = rl+El(m2)*El(m2+m)/(L-m-n);rO = rO+El(m2)A2/(L-m-n);endrm(n/m + l) = rl/rO;.word.end endsubplot。,1,1)；plot(l:10,Jn,r);title。残差平方和jn曲线图)subplot,1,2);plot(l:l:10,t/g1);title(F检验结果图)；figure(2);plot(0:9,rm(l,

17、:)/g),hold on;plot(0:9,rm(2,:)b),hold on;plot(0:9,rm(3,:)/r);title。残差白色丁阶结果图);方法二、AIC方法订阶次Data = UY;L = length(Data);U = Data(:,l);Y = Data(:,2);for n = 1:1:5N = 201-nyd(l:N,l) = Y(n+l:n+N); for i = l:Nfor l = l:l:2*n.word.FIA(iJ) = -Y(n+i-l);elseFIA(iJ) = U(2*n + i-|);end endendthita = inv(FIAl*FIA

18、)*FIA,*yd;omiga = (yd - FI A*t h ita) *(yd - FI A*t h ita)/N;AIC(n) = N*log(omiga)+4*n;endplot(AIC/-);title(AI 订阶法)； xlabel(n); ylabel(AIC); grid on第三题 广义最小二乘法:n=2;N=799;.word.u=UY(:,l);y = UY(:,2);Y=y(n + l：n + N);phil=-y(n:n+N-l);phi2=-y(n-l:n+N-2);phi3=u(n+l:n+N);phi4=u(n:n + N-l);phi5=u(n-l:n+N-

19、2)phi = phil phi2 phi3 phi4 phi5;theta = inv(phi*phi)*phi*Y;f0=10;10;while 1e(n+l:n+Nzl)=y(n+l:n+N)-phi*theta;omega(:/l) = e(n:n + N-l);omega(:/2)=-e(n-l:n+N-2);f=inv(omega*omega)*omega*e(n+l:n+N/l);if (f-fO)*(f-fO) = 10A-6*eye(N)El = Zgll-glOL*Xsthita;%计算残差 E omiga(2:N/l) = -E1(1:N-1);omiga(3:Nz2)

20、= -El(l:N-2);D = omiga*M*omiga;Fx = inv(D)*omiga*M*Zgll;thitab = Fa*omiga*Fx;Xsthita = Xsthita - thitab;F = thitab - thitabO;thitabO = thitab;.word.enddispC夏氏修正法)Xsal = Xsthita(l),Xsa2 = Xsthita(2),Xsbl =Xsthita(3)zXsb2 = Xsthita(4)夏氏改良法Data = UY;%最小二乘法，取b0=0n = 2;%根据题目,本系统为2阶系统L = length(Data);%辨识所

21、需数据的总长度N = L-n;%构造测量矩阵的数据长度FIA = zeros(N,2*n);%构造测量矩阵for i=1:N%测量矩阵赋值for I = l:n*2if(ln)FIA(iJ) = Data(n+2+i-l,l);elseFIA(iJ) = -Data(n+i-l12);endendendY = Data(n + l:n + N,2);% 输出数据矩阵thita = inv(FIA*FIA)*FIA*Y;% 计算参数矩阵.word.实验结果:实验二：相关分析法搭建对象：.word.thitaLS = thita;% 最小二乘估值m = inv(FIA,*FIA)*FIA,;for

22、 i = 1:1:1000%构造 OmegaErr = Y-FIA*thitaLS;%计算残差 Eomiga(2:N,l) = -Err(l:N-l);omiga(3:N,2) = -Err(l:N-2);F = inv(omiga*omiga)*omiga,*Err;thita = thitaLS-m*omiga*F;enddispC使用夏氏改良法辨识结果为：1)Xgal = thita(l)zXga2 = thita(2)zXgbl = thita(3)zXgb2 = thita(4)第五题 用一般递推方法辨识z(k)=-1.5*z(k-l)+4.5*u(k-l)+2.8*u(k-2)A=

23、6;xO=l;M=255;for k=1:10000x2=A*xO;.word.xl=mod (x2zM);vl=xl/256;v(:zk)=(vl-0.5)*2;xO=xl;vO=vl;endnum=v;kl=k;yl=l;y2 = l;y3 = l;y4=0;for i=l:255;xl=xor(y3zy4);x2=yl;x3=y2;x4=y3;y(i)=y4；ify(i)0.5/u(i) = -l;else u(i) = l;endyl=xl;y2=x2;y3=x3;y4=x4;lamt=l;z(2)=0;z(l)=0;for k=3:15;z(k)=-1.5*z(k-l)+4.5*u(

24、k-l)+2.8*u(k-2)+0*num(Lk);.word.end%c0=0.001 0.001 0.001,;p0=10A3*eye(3/3)；c=c0,zeros(3,14);for k=3:15;hl = -z(k-lXu(k-l),u(k-2),;x=hl*pO*hl+l*lamt;xl=inv(x);kl=p0*hl*xl;dl=z(k)-hl*c0;cl=c0+kl*dl;pl=l/lamt*(eye(3)-kl*hl)*p0;e(:,k)=e2;c(:,k)=cl;c0=cl;p0=pl;end%al=c(lz:); bl=c(2z:); b2=c(3,1);i=l:15;p

25、lot(i,al/,r,/i/a2；b,/i/blz,k,/i/b2；g,)title。递推最小二乘参数辨识)用阻尼递推算法辨识为：.word.A=6;xO=l;M=255;for k=1:10000x2=A*xO;xl=mod (x2,M);vl=xl/256;v(:,k) = (vl-0.5)*2;x0=xl;v0=vl;endbeita=0.95;u=0.95;num=v;kl=k;yl=l;y2 = l;y3 = l;y4=0;for i = 1:255;xl=xor(y3,y4);x2=yl;x3=y2;x4=y3;y(i)=y4 ;ify(i)0.5zu(i) = -l;else

26、u(i)=l;endyl=xl;y2=x2;y3=x3;y4=x4;.word.endz(2)=0;z(l)=0;for k=3:200;z(k) = -1.5*z(k-l)+4.5*u(k)+2.8*u(k-2)+0*num(l/k); endc0=0.001 0.001 0.001;p=zeros(3/3)；m=0.95;p0=10A3*eye(3z3)；c=cOzzeros(3z14);I=eye(3,3)；for k=3:15;fai=-z(k-l) u(k-l) u(k-2);p=m*I-beita*m*I+beita*inv(pO)+fai*fai;pl=inv(p);c(:/k)

27、=c(:/k-l)+beita*m*pl*c(:/k-l)-c(:/k-2)+pl*fai*z(k)- fai*c(:,k-l);p0=pl;endal=c(l；:); bl=c(2,:); b2=c(3,:);.word.ploFaLrXbLkjbZg） %title（阻尼递推最小二乘参数辨识）grid on.word.m(i,:)=UY(32-i:46-i,l); endy=UY(31:45,2);gg=ones(15)+eye(15);g = l/(25*16*2)*gg*m*y;plot(g);hold on;stem(g);实验结果： 相关分析法.word.最小二乘法建模:.word

28、.二.三次实验本次实验要完成的内容：1 .参照index2 ,设计对象,从workspace空间获取数据，取 二阶，三阶对象实现最小二乘法的一次完成算法和最小二乘法 的递推算法(LS and RLS );.对设计好的对象,在时间为200-300之间，设计一个阶跃 扰动，用最小二乘法和带遗忘因子的最小二乘法实现，对这两 种算法的特点进行说明；实验内容结果与程序代码：利用LS和RLS得到的二阶，三阶参数算法阶次A1A2A3B0B1B2B3LS二阶-0.70.13-0.00.560.31842730366857RLS二阶-0.70.13-0.00.560.31824730366857LS三阶-0.4

29、-0.10.04-0.00.560.510.1138122807078520660RLS三阶-0.4-0.10.04-0.00.560.510.1138122807078520660.word.以下给出RLS中的参数估计过程曲线和误差曲线(二阶)：M = UY(:rl);z=UY(:,2);H=zeros(199/5);for i=l:199H(i,l)=-z(i+l);H(i,2) = -z(i);H(i,3) = M(i+2);H(i,4) = M(i + l);H(L5) = M(i);endEstimate=inv(H*H)*H*(z(3:201)RLS （二阶）:.word.cicM

30、 = UY(:,1);z=UY(:,2);P=100*eye(5); %估计方差Pstore=zeros(5,200);Pstore(:,l)=P(Ll),P(2,2),P(3,3),P(4,4),P(5,5)r；Theta=zeros(5/200); %参数的估计值，存放中间过程估值 Theta(:zl) = 0;0;0;0;0;K=10;10;10;10;10;10;10;for i=3:201h=-z(i-l);-z(i-2);M(i);M(i-l);M(i-2);K=P*h*inv(h*P*h + l);Theta(:/i-l)=Theta(:/i-2) + K*(z(i)-h*The

31、ta(:J-2);P=(eye(5)-K*h)*P;Pstore(:,i-1)=PQ,1),P(2,2),P(3,3),P(4,4),P(5,5)1； end i=l:200;figure(l)plot(i/Theta(l/:)/izTheta(2z:),izTheta(3/:zi/Theta(4,:)zi,Theta( 5,：)titled待估参数过渡过程).word.figure(2)plot(i/Pstore(lz:),izPstore(2/:,i/Pstore(3/:/i/Pstore(4,:)zi/Psto re(5,:)titled估计方差变化过程)同理可以写出三阶的LS以及RLS算法,此处略去。2.在t=250出加入一个阶跃扰动，令扰动值为0.05 利用RLS和带遗忘因子的RLS计算结果如下表。算法遗忘 因子A1A2A3B0B1B2B3RLS1-0.4-0.00.03-0.00.560.480.1079191749068598233RLS0.98-0.5-0.00.02-0.00.560.450.0749839070069281577RLS的参数变化过程曲线和误差曲线如下: