function [VLs,S,VS]=estorderiv(y,u,nmax,L); //function [V,S,VS]=estorderiv(y,u,nmax,L); //function for model order estimation. The ordre is estimated using Instrumental Variable Method //with delayed inputs as the instrumental variables. The criterion has the following //forme: // 1 // V(n)=---Min{|Y(0) - Rn(n).theta|^2} n=1,2,3,...nmax // N //where // Y(0)= Zm.y(0) // Rn(n)=Zm.R(n) //inputs: //y ... vector of acquired system outputs //u ... vector of applied inputs //nmax ...maximum order to be considered //L ... size of instrumental matrix. Recommended value for L is 70, otherwise L should be >2nmax. //otputs: //V... normalized vector of the criterion (V=[V(0) V(1) V(2), ...]^T, where V(0)=1) //S... vector of penalty coefficients. Model complexity n is penalised. //VS... normalized vector of penalised criterion VS=V+S //programmed by: Hynek Prochazka // //written by: H. Prochazka, I.D. Landau //7th june 2002 sz=size(y); if(sz(2)~=1), y=y';end;//column sz=size(u); if(sz(2)~=1), u=u';end; VLs=[]; S=0; VS=[]; L=2*nmax; for nn=1:nmax, N=length(y)+1-L; //construction of R(nn) and Z Z=[]; uz=u; for k=1:2*nmax, Z=[Z uz(2*nmax-k+1:2*nmax-k+N-1)]; end; Rnny=[]; Rnnu=[]; for k=1:nn, yp=y(L-k+1:L-k+N-1); up=u(L-k+1:L-k+N-1); Rnny=[Rnny yp]; Rnnu=[Rnnu up]; end; y0=y(L+1:length(y)); //QR-triangularization [Q,R]=qr(Z); Z1=R(1:L,:); Zm=inv(Z1')*Z'; Y0=Zm*y0; R1nn=Zm*[Rnny Rnnu]; //theta computation theta_min=inv(R1nn'*R1nn)*R1nn'*Y0;//inv(R1nn)*Y0; criter=sum((Y0-R1nn*theta_min).^2)/N; VLs=[VLs;criter]; penal=0.4*nn*log10(N)/N; S=[S;penal]; end; V0=sum(Y0.^2)/N; VLs=[V0;VLs]; VLs=VLs/VLs(1); VS=VLs+S; VS=VS/VS(1);