ftest.m

遗忘因子法辨识系统参数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 遗忘因子最小二乘递推算法(RFF)
% wjx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 三阶系统的模型
% z(k)=a1*z(k-1)+a2*z(k-2)+a3*z(k-3)+b1*u(k-1)+b2*u(k-2)+b3*u(k-3)+v(k)
% 辨识出参数 theta=[a1,a2,a3,b1,b2,b3]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 清理工作间变量
clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 生成输入信号9级的M序列
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 置初始值
c=[1 0 1 0 1 1 0 1 1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% M序列的周期
len=length(c);
T=2^len-1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 循环运算
for i=1:T;
    x1=xor(c(8),c(9)); % 模2运算，也即二进制的加法运算，xor()表示异或
    x2=c(1); 
    x3=c(2); 
    x4=c(3);
    x5=c(4);
    x6=c(5);
    x7=c(6);
    x8=c(7);
    x9=c(8);
    c(i)=c(9); % 输出M序列
    u(i)=c(i);
    c(1)=x1;c(2)=x2;c(3)=x3;c(4)=x4;c(5)=x5;c(6)=x6;c(7)=x7;c(8)=x8;c(9)=x9;   
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 系统的阶次
N=3;
% z的前三个初始值置零
z(1)=0;z(2)=0;z(3)=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 输出数据(或含噪声数据)
v=rand(T,1);C=input('噪声系数C= ');
% v=whnoise(T);
for k=N+1:T;    
    %z(k)=1.4*z(k-1)-0.6*z(k-2)+0.8*z(k-3)+u(k-1)+0.7*u(k-2)+0.9*u(k-3); 
    z(k)=1.4*z(k-1)-0.6*z(k-2)+0.8*z(k-3)+u(k-1)+0.7*u(k-2)+0.9*u(k-3)+C*v(k);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 遗忘因子最小二乘递推算法(RFF)
L=input('数据长度L= ');
m=input('遗忘因子0<m<=1,m= ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for N=1:4;
% 置参数theta的初始值(很小的实向量)
theta0=10^(-5)*diag(eye(2*N,2*N));
% 置协方差阵P的初始值(a^2*I,a充分大)
p0=10^7*eye(2*N,2*N);
% 相对误差
epsilon=10^(-8);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 递推算法
for k=5:L+4;
    switch N
        case 1
            h1=[-z(k-1),u(k-1)]';
        case 2
            h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';
        case 3
            h1=[-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2),u(k-3)]';
        case 4
            h1=[-z(k-1),-z(k-2),-z(k-3),-z(k-4),u(k-1),u(k-2),u(k-3),u(k-4)]';
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 求k时刻的增益阵K(k)
    k1=p0*h1*inv(h1'*p0*h1+m);   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 求k时刻的协方差阵P(k)
    p1=(p0-k1*h1'*p0)/m; 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % 求k时刻的参数theta(k)
    t1=z(k)-h1'*theta0; 
    theta1=theta0+k1*t1;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %求损失函数J
    g0=0;
    g1=m*(g0+t1^2/(h1'*p1*h1+m));
    g0=g1;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    % 递推算法的停机标准
    e=(theta1-theta0)./theta0; 
    if abs(e)<epsilon 
        theta=theta1;
    else
        theta0=theta1;
        p0=p1;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %保存三阶系统估计参数
    if  N==3;%本程序中估计参数很不稳定，随L,r取值不同，变化很大（不作为参数辨识的依据），但能有时能正确辨识系统阶次。
        par(:,k-4)=theta1;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       
end %for k=5:L+4
%G(N)=g1+(g1==0)*eps
G(N)=g1+eps%防止出现除O的情况
end%for N=1:4
switch L
    case 100
        ta=3.09;
    case 300
        ta=3.03;
    case 500
        ta=3.01;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%绘制三阶系统参数变化图（计算不稳定，仅提供同前程序对比）
x=5:(L+4);
plot(x,par(1,:),x,par(2,:),x,par(3,:),x,par(4,:),x,par(5,:),x,par(6,:))
legend( 'a1','a2','a3','b1','b2','b3',4)
%axis([0 20 -2 2])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:3
    
    f(k)=0.5*(G(k)-G(k+1))*(L-2*k-2)/G(k+1);
end
for k=1:2
   n0=(f(k)>ta)&(f(k+1)<=ta);
    if n0>0;
       n0=k+1
   else
       disp('error')
   end
end