zgls.m

用得到的已知输入输出数据

function [IdentResult]=ZGLS(DATA,ORDER,Lamda,C)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 增广最小二乘程序                   %
% 系统模型A(q)y(t)=B(q)u(t)+C(q)e(t) %
% Y(k)=-SidaA(1)*Y(k-1)-...          %
%      -SidaA(Na)*Y(k-Na)            %
%      +SidaB(1)*U(k-1)+...          %
%      +SidaB(Nb)*U(k-Nb)            %
%      +SidaC(1)*E(k-1)+...          %
%      +SidaC(Nc)*E(k-Nc)            %
%      +E(k)                         %
% 辨识数据为列向量Y,U                %
% [DATA]=[Y,U],是必须输入的          %
% A,B,C阶次分别为Na,Nb,Nc            %
% [ORDER]=[Na,Nb,Nc],默认值为[2,2,2] %
% Lamda为遗忘因子,默认值为1          %
% C为初值常数,默认值为1000000        %
% IdentResult:辨识结果               %
% IdentResult.Sida:辨识系数          %
% IdentResult.AllSida:系数记录       %
% IdentResult.DataLength:数据长度    %
% IdentResult.StartPoint:辨识始点    %
% IdentResult.Jmin:准则函数          %
% IdentResult.Error:辨识残差         %
% IdentResult.Resdue:残差方差        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<4
    C=1000000;
end
if nargin<3
    Lamda=1;
end
if nargin<2
    ORDER=[2,2,2];
end
if nargin<1
    error('NO INPUT DATA');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y=DATA(:,1);         %输出数据
U=DATA(:,2);         %输入数据
Na=ORDER(1);         %A的阶次
Nb=ORDER(2);         %B的阶次
Nc=ORDER(3);         %C的阶次
DATALength=size(DATA,1);
iStart=max(Na,Nb)+1; %起始点
if DATALength<iStart
    error('TOO LITTLE DATA SAMPLED');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%辨识过程
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%递推最小二乘，带遗忘因子Lamda
Sida=zeros(Na+Nb+Nc,1); %初始系数
E=zeros(iStart-1,1);    %E初值
P=C*eye(Na+Nb+Nc);      %P矩阵
AllSida=[];             %系数记录
Jmin=0;                 %准则函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%下面是辨识算法
for i=iStart:DATALength
    temp1=[];temp2=[];temp3=[];
    for j=1:Na
        temp1=[temp1;-Y(i-j)];
    end
    for j=1:Nb
        temp2=[temp2;U(i-j)];
    end
    for j=1:Nc
        if i-j>0
            temp3=[temp3;E(i-j)];
        else 
            temp3=[temp3;0];
        end
    end
    Fi=[temp1;temp2;temp3];
    E(i)=Y(i)-Fi'*Sida;
    Jmin=Lamda*(Jmin+E(i)*E(i)/(Fi'*P*Fi+Lamda));
    L=P*Fi/(Lamda+Fi'*P*Fi);
    P=(eye(Na+Nb+Nc)-L*Fi')*P/Lamda;
    Sida=Sida+L*E(i);
    AllSida=[AllSida,Sida];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算残差
L=DATALength-iStart+1;       %真值数据长度
YReal=Y(iStart:DATALength);  %真实值
YHat=zeros(L,1);             %估计值
for i=1:Na
    YHat=YHat-Sida(i)*Y(iStart-i:DATALength-i);
end
for i=1:Nb
    YHat=YHat+Sida(i+Na)*U(iStart-i:DATALength-i);
end
Error=YReal-YHat;            %残差向量
Resdue=Error'*Error/L;       %残差方差
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%返回结果
IdentResult=struct('DATA',DATA,...
    'ORDER',ORDER,...
    'DataLength',DATALength,...
    'StartPoint',iStart,...
    'Sida',Sida,...
    'AllSida',AllSida,...
    'Jmin',Jmin,...
    'Error',Error,...
    'Resdue',Resdue);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Designed by Ling Qing  2002.01.16  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%