untitled21.m

来自「最小二乘法系统模型辨识」· M 代码 · 共 71 行

load V10.txt;
u1=V10;
u=u1(1:100);
load T5.txt;
z1=T5;
z=z1(1:100);
N=5;
c3=eye(1,2*N);
h2=eye(1,2*N);
d=zeros(2*N,1);
%RLS递推最小二乘法
c0=c3';
p0=10^6*eye(2*N,2*N);
%E=0.00000000005;%相对误差E=0.00000000005 
c=[c0,zeros(2*N,14)];
e=zeros(2*N,15);
for k=N+1:100
    for i=1:N
    c3(1,i)=0.001;
    h2(1,i)=-z(k-i);
end
for j=N+1:2*N
    c3(1,j)=0.001;
    h2(1,j)=u(k-j+N);
end 
    h1=h2';
    x=h1'*p0*h1+1; 
    x1=inv(x); 
    k1=p0*h1*x1;%求出K值 
    m1=z(k)-h1'*c0; 
    c1=c0+k1*m1;%辨识参数c 
   % v(k)=z(k)-h1'*c1; 
  e1=c1-c0; 
    e2=e1./c0; 
   e(:,k)=e2; 
    c0=c1; 
    c(:,k)=c1;%把递推出的辨识参数C的列向量加入辨识参数矩阵 
    p1=p0-k1*h1'*p0;%find p(k) 
    p0=p1;%给下次使用 
%if e2<=E break;%若参数收敛满足要求，终止计算 
%end%小循环结束 
end%大循环结束
d=[d,c1];
d%显示被辨识参数
%a1=d(t,2);a2=d(2,2);a3=d(3,2);a4=d(4,2);b1=d(5,2);b2=d(6,2);b3=d(7,2);b4=d(8,2);
%y(k)=-a1*y(k-1)-a2*y(k-2)-a3*y(k-3)-a4*y(k-4)+b1*u(k-1)+b2*u(k-2)+b3*u(k-3)+b4*u(k-4);
for i=1:N
y(i)=z(i);%y(2)=z(2);y(3)=z(3);y(4)=z(4);
Y(i)=z(i);%Y(2)=z(2);Y(3)=z(3);Y(4)=z(4);
end
for k=N+1:100
    for i=1:N
    c3(1,i)=0.001;
    h2(1,i)=-y(k-i);
end
for j=N+1:2*N
    c3(1,j)=0.001;
    h2(1,j)=u(k-j+N);
end 
  y(k)=h2*d(:,2);
  Y(k)=y(k);
end
Y;
e3=Y'-z;
subplot(1,3,1);
plot(1:100,z,'r');
subplot(1,3,2);
plot(1:100,Y,'g');
subplot(1,3,3);
plot(1:100,e3,':');