搜索
系统辨识作业方法一:递推最小二乘法[NUM]=xlsread('shuju','B4:B257');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.00000005;%相对误差E参考值c=[c0,zeros(4,253)];%被辨识参数矩阵的初始值及大小e=zeros(4,254);%相对误差的初始值及大小%开始递推运算fork=3:N;h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';%求h1x=h1'*p0*h1+1;x1=inv(x);k1=p0*h1*x1;%求k1c1=c0+k1*(z(k)-h1'*c0);%求ce1=(c1-c0)./c0;%求参数的相对变化e(:,k)=e1;%把当前相对变化的列向量加入误差矩阵的最后一列c0=c1;%新获得的参数作为下一次递推的旧参数c(:,k)=c1;%把当前所辨识参数的c1列向量加入辨识参数矩阵的最后一列p1=p0-k1*h1'*p0;%求p1值p0=p1;%把当前值给下次用ifnorm(e1)=Ebreak;%若参数收敛满足要求，终止计算endend%分离参数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到254plot(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到254plot(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.6418a2=-0.3451b1=0.0015b2=-0.0014ea1=-1.4483e-004ea2=3.5589e-004eb1=-0.0020eb2=-0.0021二、似然法辨识v=randn(254,1);%产生正态分布随机数V=0;%计算噪声方差fori=1:254V=V+v(i)*v(i);endV1=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;%迭代计算参数值和误差值fork=3:254h=[-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);b2(k)=o(4);d1(k)=o(5);d2(k)=o(6);e1(k)=abs(a1(k)+1.2);e2(k)=abs(a2(k)-0.6);e3(k)=abs(b1(k)-1.0);e4(k)=abs(b2(k)-0.5);e5(k)=abs(d1(k)+1.0);e6(k)=abs(d2(k)-0.2);zf(k)=z(k)-d1(k)*zf(k-1)-d2(k)*zf(k-2);uf(k)=u(k)-d1(k)*uf(k-1)-d2(k)*uf(k-2);vf(k)=v(k)-d1(k)*vf(k-1)-d2(k)*vf(k-2);hf=[-zf(k-1);-zf(k-2);uf(k-1);uf(k-2);vf(k-1);vf(k-2)];end%数%绘图figure(1)k=1:254;plot(k,a1,'k:',k,a2,'b',k,b1,'r',k,b2,'m:',k,d1,'g',k,d2,'k');xlabel('k')ylabel('参数辨识')legend('a1=-1.2,','a2=0.6','b1=1.0','b2=0.5','d1=-1.0','d2=0.2');%图标炷title('递推极大使然RML');figure(2)k=1:254;plot(k,e1,'k',k,e2,'b',k,e3,'r',k,e4,'m',k,e5,'g',k,e6,'k');xlabel('k');ylabel('误差过度过程');title('误差曲线')三、增广法辨识z=xlsread('shuju','B4:B257');u=xlsread('shuju','C4:C257');v=randn(1,254);c0=[0.0010.0010.0010.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(7,7);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.00000000005;%相对误差E=0.000000005c=[c0,zeros(7,253)];%被辨识参数矩阵的初始值及大小e=zeros(7,254);%相对误差的初始值及大小fork=3:254;%开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)]';%为求K(k)作准备x=h1'*p0*h1+1;k1=p0*h1*(inv(x));%Kd1=z(k)-h1'*c0;c1=c0+k1*d1;%辨识参数ce1=c1-c0;e2=e1./c0;%求参数误差的相对变化e(:,k)=e2;c0=c1;%给下一次用c(:,k)=c1;%把递推出的辨识参数c的列向量加入辨识参数矩阵p1=p0-k1*k1'*[h1'*p0*h1+1];%findp(k)p0=p1;%给下次用ife=Ebreak;%若收敛情况满足要求，终止计算end%判断结束end%循环结束%分离变量a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);%分离出a1、a2、b1、b2d1=c(5,:);d2=c(6,:);d3=c(7,:);%分离出d1、d2、d3ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:);%分离出a1、a2、b1、b2的收敛情况ed1=e(5,:);ed2=e(6,:);ed3=e(7,:);%分离出d1、d2、d3的收敛情况figure(1);%画第二个图形i=1:254;plot(i,a1,'r',i,a2,'r:',i,b1,'b',i,b2,'b:',i,d1,'g',i,d2,'g:',i,d3,'g+')%画出各个被辨识参数title('ParameterIdentificationwithRecursiveLeastSquaresMethod')%标题figure(2);i=1:254;%画出第三个图形plot(i,ea1,'r',i,ea2,'r:',i,eb1,'b',i,eb2,'b:',i,ed1,'g',i,ed2,'g:',i,ed2,'r+')%画出各个参数收敛情况title('IdentificationPrecision')%标题figure(3);subplot(4,1,1);%画出第四个图形，第一个子图i=1:254;plot(i,zs(i),'r')%画出被辨识系统在没有噪声情况下的实际输出响应subplot(4,1,2);i=1:254;plot(i,z(i),'g')%第二个子图，画出被辨识系统的采样输出响应subplot(4,1,3);i=1:254;plot(i,zm(i),'b')%第三个子图，画出模型含有噪声的输出响应subplot(4,1,4);i=1:254;plot(i,zs(i),'b')%第四个子图，画出模型去除噪声后的输出响应