系统辨识MATLAB.docVIP

系统辨识作业 u=NUM; [NUM1]=xlsread(shuju,C4:C257); z=NUM1; N=length(u); c0=[0.001,0.001,0.001,0.001]; %直接给出被辨识参数的初始值，取一个充分小的实向量 p0=10^7*eye(4,4); %初始状态P0也采用直接取方式，取一个充分大的实数单位矩阵 E=0 %相对误差E参考值 c=[c0,zeros(4,253)]; %被辨识参数矩阵的初始值及大小 e=zeros(4,254); %相对误差的初始值及大小 %开始递推运算 for k=3:N; h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]; %求h1 x=h1*p0*h1+1; x1=inv(x); k1=p0*h1*x1;%求k1 c1=c0+k1*(z(k)-h1*c0); %求c e1=(c1-c0)./c0;%求参数的相对变化 e(:,k)=e1; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1; %新获得的参数作为下一次递推的旧参数 c(:,k)=c1; %把当前所辨识参数的c1列向量加入辨识参数矩阵的最后一列 p1=p0-k1*h1*p0; %求p1值 p0=p1; %把当前值给下次用 if norm(e1)=E break; %若参数收敛满足要求，终止计算 end end %分离参数 a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:); ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:); figure(1); %画第1个图形 i=1:254; %横坐标从1到254 plot(i,a1,k,i,a2,b,i,b1,r,i,b2,g) %画出a1，a2，b1，b2的各次辨识结果 legend(a1,a2,b1,b2); xlabel(k); ylabel(辨识参数); %标注纵轴变量 title(最小二乘各次递推参数估计值) %图形标题 a1=c(1,254) a2=c(2,254) b1=c(3,254) b2=c(4,254) figure(2); %画第2个图形 i=1:254; %横坐标从1到254 plot(i,ea1,k,i,ea2,b,i,eb1,r,i,eb2,g); %画出a1，a2，b1，b2的各次辨识结果的收敛情况 legend(ea1,ea2,eb1,eb2); xlabel(k) %标注横轴变量 ylabel(参数误差) %标注纵轴变量 title(参数的误差收敛情况) %图形标题 ea1=e(1,254) ea2=e(2,254) eb1=e(3,254) eb2=e(4,254) 最小二乘递推辨识结果： a1 = -0.6418 a2 =-0.3451 b1 = 0.0015 b2 =-0.0014 ea1 =-1.4483e-004 ea2 =3.5589e-004 eb1 = -0.0020 eb2 = -0.0021 似然法辨识 v=randn(254,1); %产生正态分布随机数 V=0; %计算噪声方差 for i=1:254 V=V+v(i)*v(i); end V1=V/254; N=xlsread(shuju.xls) ; %取数据作为输入输出矩阵 A=N([1 : 254],[1,2,3]); x=A(:,2); y=A(:,3); u=x; z=y;%得到输入输出数据 o1=0.001*ones(6,1);p0=eye(6,6); %赋初值 zf(1)=0.1;zf(2)=0.1;vf(2)=0.1;vf(1)=0.1;uf(2)=0.1;uf(1)=0.1; %迭代计算参数值和误差值 for k=3:254 h=[-z(k-1);-z(k-2);u(k-1);u(k-2);v(k-1);v(k-2)]; hf=h; K=p0*hf*inv(hf*p0*hf+1); p=[eye(6,6)-K*hf]*p0; v(k)=z(k)-h*o1; o=o1+K*v(k) ; p0=p; o1=o; a1(k)=o(1); a2(k)=o(2); b1(k)=o(3)