function [B,A]=udrls(y,u,na,nb,d,Fin,lam1,lam0) % RLS is used to identify a discrete time model of a plant operating in % open-loop based on the recursive least squares method and the UD factorization % for updating the adaption gain matrix (F). % % [B,A]=udrls(y,u,na,nb,d,Fin,lam1,lam0) % % y and u are the column vectors containing respectively the output and the % excitation signal. % % na, nb are the order of the polynomials A,B and d is the pure time delay % % % Fin is the initial gain F0=Fin*(na+nb)*eye(na+nb) (Fin=1000 by default) % % lam1 and lam0 make different adaptation algoritms as follows: % % lam1=1;lam0=1 :decreasing gain (default algorithm) % 0.952, error('This routine is only for SISO systems'),end [nl,nc]=size(u); if nc>2, error('This routine is only for SISO systems'),end if (na<0 | nb<0 | d<0), error('The order of A,B and d should not be negative!!'),end nd=min(length(u),length(y)); % number of data nth=na+nb; if nd1 | lam0>1), error('lam1 and lam0 should be less than 1');end if (lam1<0.95 | lam0<0.95), disp ('warning :lam1 and lam0 are normally greater than 0.95');end i=[1:np-1 np+1:2*np-1]; %shift index for vector phi j=[1:na np+d+1:np+d+nb]; %phi index for theta theta=zeros(nth,1); phi=zeros(2*np,1); U=eye(nth); D=Fin*nth*eye(nth); F=U*D*U'; for t=1:nd yhat=theta'*phi(j); e_apri=y(t)-yhat; Beta(1)=1; Delta(1)=lam1; V=U'*phi(j); G=D*V; for J=1:nth Beta(J+1)=Beta(J)+ V(J)*G(J); Delta(J+1)=Delta(J)+ V(J)*G(J); D(J,J)=Delta(J)*D(J,J)/(Delta(J+1)*lam1); Gamma(J)=G(J); M(J)=-V(J)/Delta(J); for I=1:J-1 Ut=U(I,J); U(I,J)=U(I,J)+Gamma(I)*M(J); Gamma(I)=Gamma(I)+Ut*Gamma(J); end end for II=1:nth Gamma(II)=Gamma(II)/Beta(nth+1); end e_apost=e_apri/Beta(nth+1); theta=theta+F*phi(j)*e_apost; % or : theta=theta+Gamma'*e_apri; F=U*D*U'; lam1=lam0*lam1+1-lam0; phi(i+1)=phi(i);phi(1)=-y(t);phi(np+1)=u(t); end A=[1;theta(1:na)]'; B=[zeros(d+1,1);theta(na+1:na+nb)]';