系统辨识MATLAB

1、系统辨识作业方法一: 递推最小二乘法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=107*eye(4,4); %初始状态P0也采用直接取方式，取一个充分大的实数单位矩阵 E=0.00000005; %相对误差E参考值 c=c0,zeros(4,253); %被辨识参数矩阵的初始值及

2、大小 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列向量加入辨识参数矩阵的最后一列 p

3、1=p0-k1*h1'*p0; %求p1值 p0=p1; %把当前值给下次用 if norm(e1)<=E break; %若参数收敛满足要求，终止计算 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到254 plot(i,a1,'k',i,a2,'b',i,b1,'r',i,b2,'g') %画出a1，a2

4、，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&#

5、39;,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.00

6、15 b2 =-0.0014ea1 =-1.4483e-004 ea2 =3.5589e-004eb1 = -0.0020 eb2 = -0.00212、 似然法辨识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);

7、%赋初值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); b2(k)=o(4); d1(k)=o

8、(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

9、(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',

10、9;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('误差曲线')3、 增广法辨识z=xlsread('shuju',&#

11、39;B4:B257');u=xlsread('shuju','C4:C257');v=randn(1,254);c0=0.001 0.001 0.001 0.001 0.001 0.001 0.001' %直接给出被辨识参数的初始值,即一个充分小的实向量p0=106*eye(7,7); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.00000000005; %相对误差E=0.000000005c=c0,zeros(7,253); %被辨识参数矩阵的初始值及大小e=zeros(7,254); %相对误差的初始值及大小for k=3:2

12、54; %开始求K h1=-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); %K d1=z(k)-h1'*c0; c1=c0+k1*d1; %辨识参数c e1=c1-c0; e2=e1./c0; %求参数误差的相对变化 e(:,k)=e2; c0=c1; %给下一次用 c(:,k)=c1; %把递推出的辨识参数c 的列向量加入辨识参数矩阵 p1=p0-k1*k1'*h1'*p0*h1+1; %find p(k) p0=

13、p1; %给下次用 if e<=E break; %若收敛情况满足要求，终止计算 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的收

14、敛情况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('Parameter Identification with Recursive Least Squares Method') %标题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('Identification Precision') %标题figure(3); subplot(4,1,1); %画出第四个图形，第一个子图i=1:254;plot(i,zs(i),'r') %