Global Index (short | long) | Local contents | Local Index (short | long)
[eta,it_inf]=pemss(z,etainit,index,maxiter,tol,lim,maxsize,T,...
PEMSS Minimizes prediction error criteria
TH = pemss(Z,TH0)
TH: The result is returned in this matrix. See HELP THETA
Z: The datamatrix Z = [Y U], where each column of Y and U
corresponds to a particular output and input signal.
TH0: This argument provides the model structure,
the initial conditions, and the quadratic norm.
(See MS2TH, MF2TH and THINIT)
With TH = pemss(Z,TH0,INDEX) only the adjustable parameters
with indices listed in the row vector INDEX will be estimated
Remaining parameters will be fixed to their values in ETA0
Optional arguments can be reached by
TH = pemss(Z,TH0,INDEX,MAXITER,TOL,LIM,MAXSIZE,T)
See Help AUXVAR for an explanation of the these arguments.
In the m-file NUDERSTEP the stepsize to be used for numerical
differentiation is set.
| This function calls | This function is called by |
|---|---|
function [eta,it_inf]=pemss(z,etainit,index,maxiter,tol,lim,maxsize,T,...
display,gui)
% L. Ljung 10-2-90,7-25-94
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 2.5 $ $Date: 1998/05/22 01:18:24 $
global XIDiter
stopp=0;
[Ncap,nz]=size(z);
maxdef=idmsize(Ncap);
if nargin<10, gui=0;end
if nargin<9, display=0;end
if nargin<8,T=etainit(1,2);end
if nargin<7,maxsize=maxdef;end
if nargin<6,lim=1.7;end
if nargin<5,tol=0.01;end
if nargin<4,maxiter=10;end
if nargin<3,index=[];end
[nr,nc]=size(etainit);
nd=etainit(1,5);
if nd==0,error('This model structure does not have any adjustable parameters!'),end
if any(imag(z)),if lim>0,
disp(['Lim set to 0; No robustification' ...
' for complex, multioutput data'])
lim=0;
end,end
if isempty(index),index=1:nd;end
if max(index)>nd
msg=str2mat('You have specified a parameter index number that is larger',...
'than the number of adjustable parameters in the model');
if gui,errordlg(msg);iduipoin(2);end
error(msg)
end
sspm=getmfth(etainit);
if any(etainit(2,8)==[2 3]), Tmod=-abs(T);else Tmod=abs(T);end
[th,P,lambda]=th2par(etainit);
[modarg,rarg,carg]=getargth(etainit);
n=nd;
[A,B,C,D,K,X0]=feval(sspm,th,Tmod,modarg);[ny,nx]=size(C);
[nx,nu]=size(B);[ny,nx]=size(C);
rowmax=max(length(index)*ny,nx+nz);
Mloop=floor(maxsize/rowmax);
if (nu+ny)~=nz, error('Incorrect number of inputs and outputs specified in data vector'),end
ei=eig(A-K*C);
if max(abs(ei))>1.1,error('Initial guess predictor unstable!'),end
for kc=1:Mloop:Ncap
jj=(kc:min(Ncap,kc-1+Mloop));
if jj(length(jj))<Ncap,jjz=[jj,jj(length(jj))+1];else jjz=jj;end
xh=ltitr(A-K*C,[K B-K*D],z(jjz,:),X0);
yh(jj,:)=(C*xh(1:length(jj),:).'+[zeros(ny,ny) D]*z(jj,:).').';
[nxhr,nxhc]=size(xh);X0=xh(nxhr,:).';
end
e=z(:,1:ny)-yh;
V=det(e'*e/length(e));
g=ones(n,1);
%** Determine limit for robustification of criterion ***
if lim~=0, lim=median(abs(e-ones(Ncap,1)*median(e)))*lim/0.7;end
if lim==0
el=e;
else
[ne,me]=size(e);
ll=ones(ne,1)*lim;la=abs(e)+eps*ll;el=e.*sqrt(min(la,ll)./la);
clear ll,clear la
end
%if lim~=0, lim=median(abs(e-ones(Ncap,1)*median(e)))*lim/0.7;end
%if lim==0, lim=1000*ones(1,ny)*V;end
%el=min(e,ones(Ncap,1)*lim); el=max(el,-ones(Ncap,1)*lim);
V=real(det(el'*el/length(e)));
clear e, l=0; st=0;
% ** The minimization loop **
testnorm=tol+1;
while [testnorm>tol l<maxiter st~=1]
l=l+1;
if gui,drawnow,display=get(XIDiter(2),'value');end
% * Compute the Gauss-Newton direction *
if norm(lambda-eye(ny))==0,lam=eye(ny);else ...
lam=inv(sqrtm(lambda));end
[g,R,ng]=gnns(z,th,sspm,Tmod,modarg,lam,lim,maxsize,index);
testnorm=abs(ng'*g/Ncap*100); % No normalization wrt V, since that
% is implicit in the lam-normalization
% * Search along the g-direction for a lower value of V *
[Vn,th1,st,lambda]=sess(V,z,th,sspm,Tmod,modarg,lam,g,lim,Mloop,display);
if st==1,
[Vn,th1,st,lambda]=...
sess(V,z,th,sspm,Tmod,modarg,lam,ng/trace(R)*length(R),lim,Mloop,display);
end
% * Display the current values along with the previous ones *
t=zeros(n,3); t(:,1)=th1.';t(:,2)=th.';t(:,3)=g;
if display
%home,clc
disp([' ITERATION # ' int2str(l)])
disp(['Current loss: ' num2str(Vn) ' Previous loss: ' num2str(V)])
disp(['Current th prev. th gn-dir. '])
disp(t)
disp(['Norm of gn-vector: ' num2str(norm(g))])
disp(['Expected improvement: ',num2str(testnorm),' %'])
disp(['Achieved improvement: ',num2str((V-Vn)*100/V),' %'])
if st==1
disp('No improvement of the criterion possible along the search direction')
disp('Iterations therefore terminated')
end
end %if display
if gui
set(XIDiter(4),'string',int2str(l));
set(XIDiter(5),'string',num2str(Vn));
set(XIDiter(9),'string',['(',num2str(V),')']);
set(XIDiter(6),'string',num2str(testnorm));
drawnow
stopp=get(XIDiter(7),'userdata');
end
th=th1;impr=Vn-V;V=Vn;
if stopp,break,end
end
it_inf=[l,impr,testnorm];
% *** Build up the ETA-matrix ***
eta = etainit;
P1=zeros(nd,nd);
P1(index,index)=inv(R);
P=P1;
eta(1,2)=sign(etainit(1,2))*abs(T);
eta(3,1:nd)=th;
e=pe(z,eta);
lambda=e'*e/length(e);
eta(1,1)=det(lambda);eta(2,1)=eta(1,1)*(1+nd/Ncap)/(1-nd/Ncap);
eta(4:nd+3,1:nd)=P;
eta(4+nd+etainit(1,6):3+ny+nd+etainit(1,6),1:ny)=lambda;
ti=fix(clock); ti(1)=ti(1)/100;
eta(2,2:6)=ti(1:5); eta(2,7)=23;
if eta(2,8)==3,eta(2,7)=33;end