rgls.m

关于系统辨识的一些仿真程序

clc
clear
close all
var_v=1;
J=25;
u=sqrt(var_v)*randn(1,J);
z=zeros(1,J);
zf=zeros(1,J);
v=0.05*sqrt(var_v)*randn(1,J);%均值为0,方差为var_v的白噪声信号
ab=[-3.808 5.434 -3.445 0.8187 1 -3.935 5.806 -3.807 0.9362 1 0 0 0 0]';
c_1=[1 0 0 0 0];
for k=1:4
    uf(k)=u(k);
end
for k=5:J; 
    e(k)=c_1*[v(k);v(k-1);v(k-2);v(k-3);v(k-4)];
    z(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.806*u(k-2)-3.807*u(k-3)+0.9362*u(k-4)+e(k);
    zf(k)=c_1*[z(k);z(k-1);z(k-2);z(k-3);z(k-4)];
    uf(k)=c_1*[u(k);u(k-1);u(k-2);u(k-3);u(k-4)];
end
c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001]';
ce0=[0 0 0 0]';
p0=10^6*eye(9,9);
pe0=eye(4,4);
E=0.0000000001;
c=[c0,zeros(9,J-1)];%被辨识参数矩阵的初始值及大小
ce=[ce0,zeros(4,J-1)];
e=zeros(9,J);%相对误差的初始值及大小
for k=5:J 
    h1=[-zf(k-1),-zf(k-2),-zf(k-3),-zf(k-4),uf(k),uf(k-1),uf(k-2),uf(k-3),uf(k-4)]';
    h2=[-z(k-1),-z(k-2),-z(k-3),-z(k-4),u(k),u(k-1),u(k-2),u(k-3),u(k-4)]';
    x=h1'*p0*h1+1; x1=inv(x); 
    k1=p0*h1*x1;
    d1=zf(k)-h1'*c0; c1=c0+k1*d1;
    
    ez(k)=z(k)-h2'*c1;
    he=[-ez(k-1),-ez(k-2),-ez(k-3),-ez(k-4)]';
    xe=he'*pe0*he+1; xe1=inv(xe); 
    ke1=pe0*he*xe1;
    de1=ez(k)-he'*ce0; ce1=ce0+ke1*de1;
    
    e1=c1-c0;
    e2=e1./c0;
    e(:,k)=e2;        
    c0=c1;
    ce0=ce1;
    c(:,k)=c1; 
    ce(:,k)=ce1; 
    p1=p0-k1*k1'*[h1'*p0*h1+1];
    pe1=pe0-ke1*ke1'*[he'*pe0*he+1];
    p0=p1;
    pe0=pe1;
    if e2<=E break;
    end
end
c,e,ce,
a1=c(1,:); a2=c(2,:);a3=c(3,:);a4=c(4,:); b1=c(5,:); b2=c(6,:);b3=c(7,:);b4=c(8,:);b5=c(9,:);
m1=ce(1,:); m2=ce(2,:);m3=ce(3,:);m4=ce(4,:);
ea1=e(1,:); ea2=e(2,:);ea3=e(3,:);ea4=e(4,:); eb1=e(5,:); eb2=e(6,:);eb3=e(7,:);eb4=e(8,:);eb5=e(9,:); 

figure(1)
i=1:J;
plot(i,a1,'r',i,a2,'b',i,a3,'g',i,a4,'y',i,b1,'r:',i,b2,'b:',i,b3,'g:',i,b4,'y:',i,b5,':',i,m1,'r:',i,m2,'b:',i,m3,'g:',i,m4,'y:') 
title('Parameter Identification with RGLS Method')%图形标题
figure(2)
i=1:J; 
plot(i,ea1,'r',i,ea2,'g',i,ea3,'b',i,ea4,'y',i,eb1,'r:',i,eb2,'b:',i,eb3,'g:',i,eb4,'y:',i,eb5,':') 
title('RGLS Identification Precision') %图形标题

ze=zeros(1,J);
ee=zeros(1,J);
for k=5:J 
    z(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.806*u(k-2)-3.807*u(k-3)+0.9362*u(k-4)+v(k);
    v(k)=m1(k)*e(k)+m2(k)*e(k-1)+m3(k)*e(k-2)+m4(k)*e(k-3);
    ze(k)=-a1(k)*ze(k-1)-a2(k)*ze(k-2)-a3(k)*ze(k-3)-a4(k)*ze(k-4)+b1(k)*u(k)+b2(k)*u(k-1)+b3(k)*u(k-2)+b4(k)*u(k-3)+b5(k)*u(k-4)+e(k);
    ee(k)=abs((ze(k)-z(k))/z(k));
end
figure(3)
i=1:J; 
plot(i,ee,'r')
title('RGLS算法模型误差')