3order_rls1.m

来自「最小二乘递推算法RLS」· M 代码 · 共 143 行

clear;%清理工作间变量
load ZS1.txt;
load X1.txt;
M=length(ZS1);%数据长度
%u=ZS1+randint(200,1,[-5,5]);
u=ZS1;%输入数据
z=X1;%输出数据
figure(1);%第1个图形
plot(ZS1);%输入数据显示
figure(2);
plot(X1);%输出数据显示
%RLS递推最小二乘辨识
L1=30;L2=12;L3=7;%记忆长度
a1=0;%一阶核个数
a2=0;%二阶核个数
a3=0;%三阶核个数
for i=1:1:L1
      a1=a1+1;
      for j=i:1:L2
          a2=a2+1;
          for k=j:1:L3
              a3=a3+1;
          end
      end
end
c_length=a1+a2+a3;%计算总共核项数
c00=ones(c_length,1);          
c0=c00*0.0001;%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(c_length,c_length);%直接给出初始状态P0，即一个充分大的实数单位矩阵,其项数为nxn，n为未知参数个数
E=0.000000005;%相对误差E=0.000000005
c=[c0,zeros(c_length,M-L1)];%被辨识参数矩阵的初始值及大小,行数为辨识核的个数，列数为数据长度
e=zeros(c_length,M-L1+1);%相对误差的初始值及大小，大小根据c矩阵变化
HL=zeros(M-L1+1,c_length);
for k=L1:M %开始求K 
    h1=zeros(c_length,1); %分三块写出h(k)
    i=1:1:a1;
        h1(i,1)=u(k-i+1);
  counter1=0;
  %for i=(a1+1):1:(a1+a2);
    for m=1:L2
        for n=m:L2
            if m==n 
               h1(a1+1+counter1,1)=u(k-m+1)*u(k-n+1);
            else h1(a1+1+counter1,1)=u(k-m+1)*u(k-n+1)*2;
            end
            counter1=counter1+1;
        end
    end
  %end
 counter2=0;
 %for   i=(a2+a1+1):1:c_length;
    for g1=1:L3
        for g2=g1:L3
            for g3=g2:L3
                if (g1==g2)&(g2==g3)&(g3==g1)
                   h1(a1+a2+1+counter2,1)=u(k-g1+1)*u(k-g2+1)*u(k-g3+1);
                else if (g1~=g2)&(g2~=g3)&(g3~=g1)
                        h1(a1+a2+1+counter2,1)=6*u(k-g1+1)*u(k-g2+1)*u(k-g3+1);
                    else h1(a1+a2+1+counter2,1)=3*u(k-g1+1)*u(k-g2+1)*u(k-g3+1);
                    end
                end 
            counter2=counter2+1;
            end
        end 
    end
 %end
  
    HL(k-L1+1,:)=h1';
    x=h1'*p0*h1+1; x1=inv(x); %开始求K(k)
    k1=p0*h1*x1;%求出K的值
    d1=z(k)-h1'*c0; c1=c0+k1*d1;%求被辨识参数c
    e1=c1-c0;%求参数当前值与上一次的值的差值
    e2=e1./c0;%求参数的相对变化
    e(:,k-L1+1)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列       
    c0=c1;%新获得的参数作为下一次递推的旧参数
    c(:,k-L1+1)=c1;%把辨识参数c 列向量加入辨识参数矩阵的最后一列 
    p1=p0-k1*k1'*[h1'*p0*h1+1];%求出 p(k)的值
    p0=p1;%给下次用
    %if e2<=E break;%若参数收敛满足要求，终止计算
    if all(e2(:))<=E break;
    end%小循环结束
end%大循环结束
c%显示被辨识参数
e%显示辨识结果的收敛情况
%分离各不同阶参数
 i=1:a1;
    c1_1(i)=c1(i,1);

 i=1:a2;
    c1_2(i)=c1(a1+i,1);

 i=1:a3;
    c1_3(i)=c1(a1+a2+i,1);
 
   % ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:);%分离各阶参数误差 
c1_1;
c1_2;
c1_3;
figure(3);%第2个图形
plot(c1_1);
figure(4);
plot(c1_2);
figure(5);
plot(c1_3);


%二阶三维图绘制
c1_2rec=zeros(m,n);
counter3=0;
for m=1:L2
    for n=m:L2
        c1_2rec(m,n)=c1_2(1+counter3);
        counter3=counter3+1;
    end
end

for m=1:L2
    for n=m:L2
        c1_2rec(n,m)=c1_2rec(m,n)
    end
end
figure(6);
surf(c1_2rec);

%二阶图对角线绘制
for m=1:L2
    cross1(m)=c1_2rec(m,m);
end
for m=1:L2
    cross2(m)=c1_2rec(m,L2-m+1);
end
figure(7);
subplot(2,1,1);
plot(cross1);
subplot(2,1,2);
plot(cross2);

%利用计算结果重构输出波形
figure(8);
plot(X1);
hold;
Z_result=HL*c1;
plot(30:200,Z_result(1:171));