系统辨识实验报告.doc

系统辨识实验报告自动化0903班 09051302 李姣实验一、系统辨识的经典方法系统的模块如图：（1）、对系统的传递函数进行辨识。对于一阶系统而言，未加入干扰信号时，其稳定值t0=20.0，h0=42.2040，加入干扰信号后其稳定值为t=40，h1=60.4937。现在分别取两个点为y1=30%对应的实际点为h1=42.2040+（60.4937-42.2040）*30%=47.6909；根据实际测试值，选取h1=47.8909，t1=20.6，对应的y1=（47.89*09-42.2040）/（60.4939-42.2040）=0.3109所以第一个点的取值为 y1=0.3109；t1=0.6；同理可得第二个点的数值为 y2=0.8033；t2=2.7；由公式 ：可得 T=1.6750；=0；由公式 可得 k=1.82899（2）、对传递函数进行检验下面对系统的辨识结果进行验证，用一个幅值为10的阶跃信号进行验证，程序如下：num=1.82899;den=1.675,1;t=0:0.1:10;y,x,t=step(num,den,t);plot(t,10*y)grid on;title(一阶系统模型的验证);xlabel(仿真时间);ylabel(系统的响应值);set(gca,xtick,0:0.5:10);set(gca,ytick,0:1:20);所得的仿真图形如下，实际系统加入测试信号后0.5s，从workspace中可发现系统的响应值为h=47.0929-42.2040=4.8889；验证是的对应仿真值为h=4.4720；其误差大小为：（4.8889-4.4720）/4.8889*100%=8.536%；同理，当仿真时间为3.8s时，h=16.3993；h=16.398；误差大小为： （16.3993-16.398）/16.3993*100%=0.08%；所以经过验证个，可以确定该辨识结果可以反应该系统的传递函数。实验二相关分析法aa=5;NNPP=15;ts=2;RR=ones(15)+eye(15);UU=UY(31:45,1);UY(30:44,1);UY(29:43,1);UY(28:42,1);UY(27:41,1); UY(26:40,1);UY(25:39,1);UY(24:38,1);UY(23:37,1);UY(22:36,1); UY(21:35,1);UY(20:34,1);UY(19:33,1);UY(18:32,1);UY(17:31,1);YY=UY(16:30,2);GG=(RR*UU*YY+4.4474)/aa*aa*(NNPP+1)*ts; plot(0:2:29,GG)hold onstem(0:2:29,GG,filled)（1） 最小二乘二阶一次完成算法HL=;N=15;m序列的周期n=2;系统的阶次for k=16:15+N a=-UY(n+k-1:-1:k,2) UY(n+k-1:-1:k,1); HL=HL;a; end ZL=UY(n+16:n+N+15,2);c1=HL*HL; c2=inv(c1); c3=HL*ZL; c=c2*c3; a1=c(1), a2=c(2), b1=c(3), b2=c(4)； 求得：a1=-0.7729;a2=0.1522;a3=0.5586;a4=0.3232;（2）最小二乘三阶的一次完成算法HL=;N=15;m序列的周期n=3;系统的阶次for k=16:15+N a=-UY(n+k-1:-1:k,2) UY(n+k-1:-1:k,1); HL=HL;a; end ZL=UY(n+16:n+N+15,2);c1=HL*HL; c2=inv(c1); c3=HL*ZL; c=c2*c3; a1=c(1), a2=c(2), a3=c(3), b1=c(4)；b2=c(5);b3=c(6); 求得相关系数为： a1=-0.1907; a2=-0.2461; a32=0.0392; b1=0.5609; b2=0.6509; b3=0.2252；（3）最小二乘法二阶一般递推算法%RLS递推最小z=UY(:,2);u=UY(:,1);c0=0.001 0.001 0.001 0.001；直接给出参数的初始值p0=104*eye(4,4); 直接给出初始状态P0，已给很大的单位实矩阵c=c0,zeros(4,198); 存储各部迭代后的参数值 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)=c1; c0=c1;更新C0 p0=p1;%更新C1enda1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); i=1:200;plot(i,a1,r,i,a2,b,i,b1,f,i,b2,k) title(递推最小二乘法参数辨识 )（3）最小二乘法三阶递推算法z=UY(:,2);u=UY(:,1);c0=0.001 0.001 0.001 0.001,0.001,0.001;p0=104*eye(6,6);c=c0,zeros(6,198);for k=4:200; h1=-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2),u(k-3); 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(6)-k1*h1)*p0; c(:,k-2)=c1; c0=c1;% p0=p1;%enda1=c(1,:); a2=c(2,:); a3=c(3,:); b1=c(4,:); b2=c(5,:);b3=c(6,:);i=1:199;plot(i,a1,r,i,a2,b,i,a3,black,i,b1,g,i,b2,y,i,b3,r) title(递推最小二乘参数辨识)（2）加阶跃扰动后的参数辨识带遗忘因子的最小二乘法clc;lamt=0.95;z=UY(:,2);u=UY(:,1);c0=0.001 0.001 0.001 0.001;p0=104*eye(4,4); c=c0,zeros(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(lamt=0.95的遗忘最小二乘法 );grid onlamt=1;z=UY(:,2);u=UY(:,1);c0=0.001 0.001 0.001 0.001;p0=104*eye(4,4); c=c0,zeros(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;%?C1enda1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); figure(2);i=1:199;plot(i,a1,r,i,a2,b,i,b1,y,i,b2,black) ;title(最小二乘法);grid on（3）搭建的对象为广义最小二乘法M=UY(1:800,1);v=randn(1,800);e=;e(1)=v(1);e(2)=v(2);for i=3:800e(i)=0*e(i-1)+0*e(i-2)+v(i);endz=UY(1:800,2);zf=;zf(1)=-1;zf(2)=0;for i=3:800zf(i)=z(i)-0*z(i-1)-0*z(i-2);enduf=;uf(1)=M(1);uf(2)=M(2);for i=3:800uf(i)=M(i)-0*M(i-1)-0*M(i-2);endP=100*eye(4); Theta=zeros(4,800);Theta(:,2)=3;3;3;3;K=10;10;10;10；PE=10*eye(2);ThetaE=zeros(2,800);ThetaE(:,2)=0.5;0.3;KE=10;10;for i=3:800h=-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;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;endi=1:800;figure(1);plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)title( )figure(3);plot(i,ThetaE(1,:),i,ThetaE(2,:);figure(3);plot(i,ThetaE(1,:),i,ThetaE(2,:);M=UY(1:800,1); v=randn(800);z=UY(1:800,2);P=100*eye(6);Theta=zeros(6,799);Theta(:,1)=3;3;3;3;3;3;K=10;10;10;10;10;10;for i=3:800h=-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;endi=1:799;figure(1)plot(i,Theta(1,:),i,Thet