untitled3.m

广义最小二乘辨识的matlab实现

%%%%%%%%%%%%%广义最小二乘辨识的matlab实现

clear;   
N=225;   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
%        以下利用同余法生成白噪的过程              %   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
A=6;   
x_0=1;   
M=255;   
f=2;   
NN=450;    %共生成450个噪声值   
   
for i=1:NN   
 x_A=A*x_0;   
 xi=mod(x_A,M);   
 vi=xi/256;   
 v(i)=(vi-0.5)*f;   
 x_0=xi;   
end   
%plot(v);   
noise=v(1:N);              %全程噪声v，前期噪声noise，由广义最小二乘法辨识   
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
%        以下利用System Identification Toolbox     %   
%       工具箱提供系统辨识的输入信号函数idinput    %   
%                  生成M序列的过程                 %   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
u=idinput([15 1 30],'prbs',[0 1],[-1 1]);        %u采用15拍M序列   
u=u';  
input=u(1:N);             %全程输入u，前期输入input，由广义最小二乘法辨识  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%      以下利用Model Predictive Control Toolbox    %  
%         工具箱建立数学模型并模拟输入的过程       %  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
output_0=zeros(1,NN);                      %定义输出观测值的长度450  
n=2;                                         %A和B的最高阶次  
A=[1 -1.5 0.7];  
B=[0 2.31 0.5];  
p=2;                                         %C的阶数  
C=[1 0.125 0.312];                           %设定ai bi ci  
m=idpoly(A,B,C,1,1);                         %建立数学模型  
output_0=sim(m,[u' v']);                    %得到输出值  
output_0=output_0';   
output=output_0(1:N);                    %全程输出output_0，前期输出output   
%plot(output_0);   
m   
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
%        以下开始利用广义最小二乘法辨识参数        %   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    for i=1:n   
        x(:,i)=-output((n+1-i):(N-i));   
    end   
    for i=(n+1):(2*n+1)   
        x(:,i)=input((2*n+2-i):(N+n+1-i));   
    end                                         %得到X   
   
    y=output((n+1):N);                         %得到y   
   
    ab=inv(x'*x)*x'*y';                         %第一次最小二乘  
  
    a=ab(1:n);a=[1;a];  
    b=ab((n+1):(2*n+1));                        %分别提取A B  
  
    delta=ab;  
    ab_1=ab;                    
     
while norm(delta)>1E-6                                      %norm：计算范式  
  
        ek=(sim(idpoly(1,a'),output')-sim(idpoly(1,b'),input'))';        %残差模型    
         
        e=(ek((n+1):N))';  
                           
        om=zeros(N-n,p);  
        for i=1:p  
            om(:,i)=-ek((n+1-i):(N-i));        %残差矩阵om  
        end  
  
        c=[1;inv(om'*om)*om'*e];                %第二次最小二乘得c  
         
        uk=(sim(idpoly(1,c'),input'))';   
        yk=(sim(idpoly(1,c'),output'))';  
  
         xx=zeros(N-n,(2*n+1));  
         y=(yk((n+1):N))';                          %新的y   
   
         for i=1:n                             %得新的X   
             xx(:,i)=-yk((n+1-i):(N-i));   
         end   
         for i=(n+1):(2*n+1)   
             xx(:,i)=uk((2*n+2-i):(N+n+1-i));   
         end   
   
    ab=inv(xx'*xx)*xx'*y;                         %%第三次最小二乘，产生新的ab,提取到a,b   
    a=ab(1:n);a=[1;a];   
    b=ab((n+1):(2*n+1));   
   
    delta=ab-ab_1;                           %产生新的delta   
    ab_1=ab;                                %记录本次计算的ab   
end   
   
y1=(sim(idpoly(a',b',c'),[input' noise']))';     %广义最小二乘估计得到的拟和值   
subplot(211);   
P1=plot(1:225,output,'-r',1:225,y1,':b');                                %真实输出值   
legend(P1,'真实输出值','拟和输出值');         
title('广义最小二乘估计');                 
hold on   
a   
b   
c   
norm(y1-output)     %计算估计的损失