系统辨识实验报告 (2)

1、.自动化09-3 宋佳瑛 09051304系统辨识实验报告实验一：系统辨识的经典方法实验目的：掌握系统的数学模型与系统的输入，输出信号之间的关系，掌握经典辨识的实验测试方法和数据处理方法。熟悉matlab/Simulink 环境。实验内容：1.用系统阶跃响应法测试给定系统的数学模型。在系统没有噪声干扰的条件下通过测试系统的阶跃响应获得系统的一阶加纯滞后或二阶加纯滞后模型，对模型进行验证。2.在被辨识的系统加入噪声干扰，重复上述1的实验过程。1. 没有噪声搭建对象 测试对象流程图实验结果为：2、加入噪声干扰搭建对象实验结果：加入噪声干扰之后水箱输出不平稳，有波动。实验二：相关分析法搭建对象：处理

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);实验结果：相关分析法最小二乘法建模：二、三次实验本次实验要完成的内容：1参照index2，设计对象，从workspace空间获取数据，取二阶，三阶对象实现最小二乘法的一次完成算法和最小二乘法的递推算法（LS and RLS）;2对设计好的对象，在时间为200-300之间，设计一个阶跃扰动，用最小二乘法和带遗忘因子的最小二乘法实现，对这两种算法的特点进行说明；实

3、验内容结果与程序代码：利用LS和RLS得到的二阶，三阶参数算法阶次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中的参数估计过程曲线和误差曲线程序清单：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); %估计方差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); %参数的估计值，存放中间过程估值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(待估参数过渡过程)figure(2)plot(i,P

6、store(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:)title(估计方差变化过程)同理可以写出三阶的LS以及RLS算法，此处略去。2.在t=250出加入一个阶跃扰动，令扰动值为0.05利用RLS和带遗忘因子的RLS计算结果如下表。算法遗忘因子A1A2A3B0B1B2B3RLS1-0.4791-0.09170.0349-0.00680.56590.48820.1033RLS0.98-0.5498-0.03900.0270-0.00690.56280.45150.0777RLS的参数变化过程曲线和误差曲线如下：带0.9

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

8、S的对象模型，改写参照程序，实现相应的算法。利用GLS算法求出参数如下：a1a2b1b2-0.74110.08110.23640.1040即y(k)=0.7411y(k1)-0.0811y(k-2)+0.2364u(k-1)+0.1040u(k-2);利用ELS求解，得到如下的结果参数a1a2b1b2d1d2ELS-0.74240.08190.23460.103900绘制出ELS算法下参数的变化曲线以及误差曲线如下：实验程序：GLS%广义最小二乘法辨识%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 %以上为最小二乘计算Theta估计值%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; %以上为最小二乘计算f估计值%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);%以上为数据滤波%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、增广最小二乘的递推算法%=产生均值为0，方差为1 的高斯白噪声=v(1)=0;v(2)=0;z=UY(:,2);M=UY(:,1);%递推求解P=100*eye(6); %估计方差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); %参数的估计值，存放中间过程估值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(参数a1、a2、b1、b2、d1、d2 估计结果：)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(待估参数过渡过程)figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:),i,Pstore(6,:)title(估计方差变化过程)大作业程序第一题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(脉冲响应序列)求系统的脉冲传递函数 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 相关二步法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(相关二步法参数过度过程);grid on最小二乘递推算法相关参数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(递推最小二乘法参数辨识 );gr

19、id on带遗忘因子的最小二乘法将上述程序中的lamt=0.9；辨识系统阶次方法一、利用残差最小辨识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检验结果图);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(残差白色丁阶结果图);方法二、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;%计算残差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;%根据题目，本系统为2阶系统L = length(Data);%辨识所需数据的总长度N = L-n;%构造测量矩阵的数据长度FIA = zeros(N,2*n);%构造测量矩阵for i = 1:N %测量矩阵赋值 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);%输出数据矩阵thita = inv(FIA*FIA)*FIA*Y;%计算参数矩阵thitaLS = thita;%最小二乘估值m = inv(FIA*FIA)*FIA;for i = 1:1:1000%构造Omega Err = Y-FIA*thitaLS;%计算残差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(使用夏氏改良法辨识结果为：)Xga1 = thita(1),Xga2 = thita(2),Xgb1 = thita(3),Xgb2 = thita(4)第五题用一般递推方法辨识 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+