function [wlossf,ulossf,wrni,urni,wyhat,uyhat]=olvalid(B,A,C,y,u) % OLVALID is used for validation of plant models identified in open loop. % % [wlossf,ulossf,wrni,urni,wyhat,uyhat]=olvalid(B,A,C,y,u) % % The estimation of the auto-correlation of the % 1-step prediction error and the estimated output. % % The estimation of the cross-correlation function between the % output error prediction and the corresponding prediction error. % % The upper bound for the cross-correlation and autocorrelation estimations is given % for a confidence interval of %97. % % The plant output y and the estimated ones (1-step predictor and output error predictor) % are compared and the loss functions are computed. % % B and A are respectively the numerator (including the pure time delay) and denominator of the plant % model to be validated. % % C is the numerator of the noise model. If no noise model is avalaible C must be set to 1. % % y and u are the column vectors of output and input signals in % open loop operation. % % % wyhat : Open loop 1-step estimated output % wlossf : Loss function (the sum of (y-wyhat)^2 divided by number of data) % % uyhat : Open loop output error prediction % ulossf : Loss function (the sum of (y-uyhat)^2 divided by number of data) % % wrni : Auto-correlations % urni : Uncorrelations % % written by: G. Zito, I.D. Landau % Last modified on November 12, 2004 na=length(A); nB=length(B); d=-1;i=1; while B(i)==0,d=d+1;i=i+1;end if d==-1, error('B should contain at least one zero!');end if d==nB-1, error('B should be a nonzero vector!');end nb=nB-d-1; if isempty(C) nC=0; else nC=length(C); end nd=min(length(u),length(y)); % number of data nc=max(max(nC,max(na-1,nb)),4); th_lim=num2str(2.17/sqrt(nd)); disp([' ']) disp('Whiteness Test') [wyhat,err_pred]=modelsim(B,A,C,y,u,d,nd); [wlossf,wrni]=correl(err_pred,err_pred,nc,nd) disp(['Theoretical upper limit : ' th_lim]) disp([' ']) if (max(wrni)>0.15) beep; disp('WARNING!'); disp('The model does not pass the whiteness test'); disp([' ']) end %%%%%%%%%%%%%%%%%%%%% figure; x1=1:nc; xx=0;yy=0; for i=1:nc,xx=[xx x1(i) x1(i) x1(i)];yy=[yy 0 wrni(i) 0];end xx=[xx nc+1];yy=[yy 0]; plot(xx,yy,0:nc+1,2.17/sqrt(nd)*ones(1,length(0:nc+1)),0:nc+1,0.15*ones(1,length(0:nc+1)),'-') hold on xlabel('i'); ylabel('WRN(i)'); title('Whiteness test for the OL identification') legend('Correlation terms',strcat('Theoretical limit: ',th_lim),'Practical limit: 0.15',0) %%%%%%%%%%%%%%%%%%%%% disp([' ']) disp('Uncorrelation Test') disp([' ']) sys_hat=filt(B,A,1); uyhat=lsim(sys_hat,u); err_dec=uyhat-y; ulossf=err_dec'*err_dec/nd [normr,urni]=correl(err_dec,uyhat,nc,nd); urni disp(['Theoretical upper limit : ' th_lim]) disp([' ']) r0=abs(err_dec'*uyhat/nd); % Warning if max(URNI) > 0.15 if (max(urni)>0.15) beep; disp('WARNING!'); disp('The model does not pass the uncorrelation test'); disp([' ']) end %%%%%%%%%%%%%%%%%%%%% figure; x1=1:nc; xx=0;yy=0; for i=1:nc,xx=[xx x1(i) x1(i) x1(i)];yy=[yy 0 urni(i) 0];end xx=[xx nc+1];yy=[yy 0]; plot(xx,yy,0:nc+1,2.17/sqrt(nd)*ones(1,length(0:nc+1)),0:nc+1,0.15*ones(1,length(0:nc+1)),'-') hold on xlabel('i'); ylabel('URN(i)'); title('Uncorrelation test for the OL identification') legend('Correlation terms',strcat('Theoretical limit: ',th_lim),'Practical limit: 0.15',0) %%%%%%%%%%%%%%%%%%%%% % Warning if ( max(WRNI) and max(URNI) ) > 0.15 if (max(wrni)>0.15 & max(urni)>0.15) beep; disp('WARNING!'); disp('The model has to be improved'); end function [y_pred,err_pred]=modelsim(B,A,C,y,u,d,nd) na=length(A)-1; nb=length(B)-1-d; nc=length(C)-1; nth=na+nb+nc; np=max([na+1,nb+d,nc+1]); i=[1:np-1 np+1:2*np-1 2*np+1:3*np-1]; %shift index for vector phi j=[1:na np+d+1:np+d+nb 2*np+1:2*np+nc]; %phi index for theta theta=[A(2:end) B(2+d:end) C(2:end)]'; phi=zeros(3*np,1); for t=1:nd y_pred(t)=theta'*phi(j); err_pred(t)=y(t)-y_pred(t); phi(i+1)=phi(i);phi(1)=-y(t);phi(np+1)=u(t);phi(2*np+1)=err_pred(t); end err_pred=err_pred'; function [normr,rni] = correl(x,y,nc,nd) switch nargin case 4 case 3 case 2 nc=4; case 1 y=x; nc=4; case 0 error('This function needs at least one argument') otherwise error('Too many arguments') end normr=sqrt((y'*y)*(x'*x))/nd; for i=1:nc ri(i)=x(1+i:nd)'*y(1:nd-i)/(nd-i); rni(i)=abs(ri(i)/normr); end