系统辨识报告.docx

1、系统辨识报告实验一：系统辨识的经典方法实验目的：掌握系统的数学模型与系统的输入，输出信号之间的关系，掌握经典辨识的实验测试方法和数据处理方法。熟悉matlab/Simulink 环境。实验内容：1.用系统阶跃响应法测试给定系统的数学模型。在系统没有噪声干扰的条件下通过测试系统的阶跃响应获得系统的一阶加纯滞后或二阶加纯滞后模型，对模型进行验证。2.在被辨识的系统加入噪声干扰，重复上述1的实验过程。实验内容：利用非线性水槽模型搭建单水槽系统模型，如下图由上图及其计算的H：H为一二元数组，分别表示第一个、第二个水箱的液位。二实验方法：运用阶跃响应法：第一个水箱的参数辨识（一阶）：一阶惯性环节的传递函

2、数为：其中： 将H归一化： 在H中查得时对应T=2.3故模型为 第二个水箱的参数辨识（二阶）：二阶系统的传递函数为：其中：在H中可得： 故有：由公式可求出：故： 实验二：相关分析法搭建对象：处理程序：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空间获取数据，取二阶，三阶对象实现最小二乘法的一

3、次完成算法和最小二乘法的递推算法（LS and RLS）;2对设计好的对象，在时间为200-300之间，设计一个阶跃扰动，用最小二乘法和带遗忘因子的最小二乘法实现，对这两种算法的特点进行说明；实验三 最小二乘法参数估计一 实验内容（1）参照index2，设计对象，从workspace空间获取数据，取二阶，三阶对象实现最小二乘法的一次完成算法和最小二乘法的递推算法（LS and RLS）;（2）对设计好的对象，在时间为200-300之间，设计一个阶跃扰动，用最小二乘法和带遗忘因子的最小二乘法实现，对这两种算法的特点进行说明（3）参照index3，设计符合GLS和ELS的对象模型，改写参照程序，实

4、现相应的.算法。二实验源程序 内容（1）：（a）% 最小二乘法的一次完成算法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);endEstimate=inv(H*H)*H*(z(3:102)Estimate = -0.7867 0.1388 0.0007 0.57080.3116（b）%最小二乘的递推算法M=UY(:,1);z=UY(:,2);P=1000*eye(5); % Pstore=zeros(5,200);% Pstore(:,1)=P

5、(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=zeros(4,400); %K=10;10;10;10;10;for 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

6、:200;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:)title()grid 内容（2）：(b)%带遗忘因子的最小二乘法M=UY(:,1);z=UY(:,2);P=1000*eye(5); %Theta=zeros(5,200); %Theta(:,1)=0;0;0;0;0;K=zeros(4,400); %K=10;10;10;10;10;lamda=0.99;%for i=3:201h=-z(i-1);-z(i-2);M(i);M(i-1);M(i-2);K=P*h*inv(h*

7、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(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:)title()grid内容（3）：(a)%广义最小二乘法辨识GLS:clcn=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:n+N

8、-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 breakendf0=f;%f%for i=n+1:n+Nyy(i,1)=y(i),y(i-1),y(i-2)*1,f(1)

9、,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);endftheta(b)%增广最小二乘法的递推算法（EM）：u=UY(:,1);y=UY(:,2);n=2;N=799;thitak=ones(7,1)*0.001;p1=eye(7,7)*(

10、1.0e6); p2=zeros(7,7);K=zeros(7,1); e=zeros(801,1); I=eye(7,7);for i=3:801phi1=-y(i-1);phi2=-y(i-2);phi3=u(i);phi4=u(i-1);phi5=u(i-2);phi6=e(i-1);phi7=e(i-2);phi=phi1 phi2 phi3 phi4 phi5 phi6 phi7; K=p1*phi/(phi*p1*phi+1); p2=(I-K*phi)*p1; thita(:,i)=thitak+K*(y(i,1)-phi*thitak);p1=p2;thitak=thita(:

11、,i);e(i)=y(i)-phi*thitak;endthitakthitak = -0.7426 0.0821 0.0000 0.2364 0.1038 0.0010 0.0010三.实验结果内容（1）（b）最小二乘的递推算法的实验结果 内容（2）（a）用最小二乘法得出的辨识结果(b)用带遗忘因子的最小二乘法得到的结果当=0.99时：当=0.98时：当=0.97时： 当=0.95时：当=0.93时：当=0.92时：内容（3）：（a）广义最小二乘法辨识GLS得到的结果:f = -0.3506 0.2873theta = -0.7411 0.0809 0.23660.1041 则系统的差分方程

12、为(b)增广最小二乘法的递推算法（EM）得到的结果：四实验结果分析比较以上六个图可以知道：加入扰动时，最小二乘法的递推算法不能将扰动引起的干扰消除，其辨识出来的参数不能迅速复原；而渐消记忆的递推算法在扰动加入后，能较快复原，且随着遗忘因子的减小，其效果越明显；五附录：系统辨识大作业源程序清单作业一：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;X0=ones(10,1);A=X.2,X,X0;E=inv(A*A)*A*Y;YE=A*E;Dif=abs(Y)-ab

13、s(YE);plot(X,Y,-rs,X,YE,-gs,X,Dif,.);legend(观测值Y,拟合值YE,差值Y-YE);xlabel(x);ylabel(y);title(作业一 y=ax2+bx+c参数求解);作业二：%利用相关法求解系统的脉冲响应。Np=15;a=5;delta=2;C=4*Np+1;for i=1:Np C=C-1; for j=1:Np U(i,j)=UY(C+j-1,1); end Y(i)=UY(i+4*Np-1,2);endR=ones(15)+eye(15);G=R*U*Y/(a2*(Np+1)*delta);t=0:2:2*Np-2;stem(t,G,g

14、s)hold onplot(t,G,-)%利用脉冲相应求解系统G（z）n=2;%系统阶次是2ga=G(2+i) G(3+i);G(3+i) G(4+i);y=-G(4+i:5+i);az=inv(ga)*y;r=1 0 0;az(2) 1 0;az(1) az(2) 1;gb=G(1+i:3+i);bz=r*gb;Gz=tf(bz,1 az,2); %利用脉冲相应求解系统G（s）as=ga(-G(1:2);x=roots(as(2) as(1) 1);s=log(x);cs=1 1;x(1) x(2)G(1:2);den=conv(1 -s(1),1 -s(2);num=cs(1)*cs(2)

15、;Gs=tf(num,den) %利用相关最小二乘法求系统的参数 clc z=UY(:,2);M=UY(:,1);%递推求解P=100*eye(4); %估计方差Pstore=zeros(4,301);Pstore(:,1)=P(1,1),P(2,2),P(3,3),P(4,4);Theta=zeros(4,301); %参数的估计值，存放中间过程估计值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:301h=-z(i-1);-z(i-2);M(i

16、-1);M(i-2);hstar=M(i-1);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;Pstore(:,i-1)=P(1,1),P(2,2),P(3,3),P(4,4);end%=结果输出及作图=disp(参数 a1 a2 b1 b2的估计结果：)Theta(:,301)i=1:301;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Th

17、eta(4,:)title(参数过渡过程)figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)title(估计方差变化过程) %最小二乘递推算法(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=zer

18、os(4,400); %增益矩阵K=10;10;10;10;10;10;10;for 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,

19、Theta(5,:)title(待估参数过渡过程)figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:)title(估计方差变化过程)Theta%带遗忘因子的RLS算法clcM=UY(:,1);z=UY(:,2);larmda=0.99;%遗忘因子 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); % 参数的估计值，存放中间

20、过程估值Theta(:,1)=0;0;0;0;0;K=10;10;10;10;10;10;10;for i=3:201h=-z(i-1);-z(i-2);M(i);M(i-1);M(i-2);K=P*h*inv(h*P*h+larmda);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h*Theta(:,i-2);P=(eye(5)-K*h)*P/larmda;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

21、(3,:),i,Theta(4,:),i,Theta(5,:)title(待估参数过渡过程)figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:)title(估计方差变化过程)Theta%阶次辨识(残差最小)%clcy=UY(:,2);u=UY(:,1);%subplot(2,2,1);for n=1:5 %对各阶次计算Jn N=301-n; Y=y(n+1:n+N,1); for j=1:N for k=1:n phi(j,k)=-y(n+j-k); phi(j,n+k)=u(n+

22、j-k); end end theta=inv(phi*phi)*phi*Y; e(n+1:n+N,1)=y(n+1:n+N)-phi*theta; J(n,1)=e*e; %计算Jn phi=0; m(n)=n;endplot(m,J,r-o);title(数据阶次曲线图);hold ongrid on%*阶次计算AIC*%Data=UY;%按Akkaike信息准则定阶L = length(Data);U = Data(:,1);Y = Data(:,2);for n = 1:1:5 N = 301-n yd(1:N,1) = Y(n+1:n+N); for i=1:N for l=1:1:

23、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;endplot(AIC,-);title(AIC随阶次的变化曲线);xlabel(n);ylabel(AIC);作业三:广义最小二乘法（GLS）：%广义最小二乘法辨识clcn=2;N=799;% y=y1;% u=u1;y=UY(:,2);u=UY(:,1);

24、% y=y3;% u=u3;Y=y(n+1:n+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)

25、;if (f-f0)*(f-f0)0.0001 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*y

26、y(n+1:n+N);endfthetaELS：%增广最小二乘的递推算法%=产生均值为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

27、-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);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,:),i,Theta(6,:)title(待估参数过渡过程)figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Ps