function [B,A]=foloe(y,u,na,nb,d,A0,Fin,lam1,lam0) // FOLOE is used to identify a discrete time model of a plant operating in // open-loop based on the output error method with filtered observation. // // [B,A]=foloe(y,u,na,nb,d,A0,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 // // A0 is the denominator of the initial model (fixed compensator) // // 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 nd= 1, error('A0 should be stable!'); end na0=length(A0)-1; np=max([na+1,nb+d,na0+1]); if ~exists('A0'), error('This routin needs more parameters!');end; if ~exists('Fin'), Fin=1000;end; if ~exists('lam1'), lam1=1;end; if ~exists('lam0'), lam0=1;end; if isempty(Fin), Fin=1000;end if isempty(lam1),lam1=1;end if isempty(lam0), lam0=1;end if (lam1>1 | 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 theta=zeros(nth,1); phi=zeros(2*np,1); phi_f=zeros(2*np,1); F=Fin*nth*eye(nth,nth); 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 for t=1:nd yhat=theta'*phi(j); e_apri=y(t)-yhat; e_apost=e_apri/(1+phi_f(j)'*F*phi_f(j)); theta=theta+F*phi_f(j)*e_apost; F=1/lam1*(F-((F*phi_f(j)*phi_f(j)'*F)/(lam1+phi_f(j)'*F*phi_f(j)))); lam1=lam0*lam1+1-lam0; yhat=theta'*phi(j); // Filtering phi by 1/A phi(i+1)=phi(i);phi(1)=-yhat;phi(np+1)=u(t); yuf=[phi(1) phi(np+1)]-[A0(2:na0+1) zeros(1,np-na0-1)]*[phi_f(1:np-1) phi_f(np+1:2*np-1)]; phi_f(i+1)=phi_f(i);phi_f(1)=yuf(1);phi_f(np+1)=yuf(2); end A=[1;theta(1:na)]'; B=[zeros(d+1,1);theta(na+1:na+nb)]';