www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xregoptmgr/nlleastsq.m
function F= nlleastsq(F,~) %NLLEASTSQ general nonlinear least squares for parameter estimation % % F= nlleastsq(F,m) % Copyright 2000-2009 The MathWorks, Inc. and Ford Global Technologies, Inc. % function handles to run and setup functions F.RunFcn = @i_FitModel; % Options which user can change F= AddOption(F,'GradObj','off','on|off','',0); F= AddOption(F,'Jacobian','on','on|off','',0); F= AddOption(F,'LargeScale','on','on|off','',0); F= AddOption(F,'JacobPattern',[],'sparse','',0); F= AddOption(F,'TypicalX',[],'numeric','',0); F= AddOption(F,'PrecondBandWidth',0,{'int',[0,Inf]},'',0); % These options can be change by the optimisation manger's GUI setp F= AddOption(F,'Display','none','none|iter|final','',1); F= AddOption(F,'MaxFunEvals',1000,{'int',[1 Inf]},'Max. Function Evaluations',1); F= AddOption(F,'MaxIter',100,{'int',[1 Inf]},'Max. Iterations',1); F= AddOption(F,'TolFun',1e-6,{'numeric',[1e-12 Inf]},'Function Tolerance',1); F= AddOption(F,'TolX',1e-6,{'numeric',[1e-12 Inf]},'Variable Tolerance',1); F= AddOption(F,'MaxError',10,{'numeric',[0 Inf]},'Max. Error',1); F= AddOption(F,'ScaleResiduals',1,'boolean','Scale Residuals',1); F.name= 'Nonlinear least squares'; % function [UpdatedContextObj,costFcn,OK,varargout] = optimRunFcnTmpl(omgr,ContextObj,InitVals,varargin) function [m,cost,OK,xf,res,J] = i_FitModel(m,F,x0,x,y,varargin) % get model constraints % do we want a constraints manager class to do this if size(y,1)<length(x0) % protect cost = NaN; OK = false; xf=[]; res = []; J = []; return end [LB,UB,A,c,nlfunc]= getConstraints(F); % Put Optimisation ToolBox parameters and values into a cell matrix otbxParam= [{F.foptions(1:end-2).Param} ; {F.foptions(1:end-2).Value}]; if get(F,'ScaleResiduals') c0= sqrt( ssenonlin(x0,F,1,m,x,y,[],varargin{:}) ); if c0==0 % try the norm of y data c0= sqrt(y'*y); if c0==0 % use 1 c0=1; end end % c0 = c0/sqrt(size(y,1)); x0=x0(:); sx= abs(x0); sx(sx < max(sx)*1e-8 )= 1; sx= 10.^( floor(log10(sx)) ) ; if max(sx)/min(sx)>1e3 % scale if range of parameters too great x0= x0./sx; % scale constraints if ~isempty(A) nv=length(x0); A= A*spdiags(sx,0,nv,nv); end if ~isempty(LB) LB= LB./sx; end if ~isempty(UB) UB= UB./sx; end else sx=[]; end else c0=1; sx= []; end if isempty(A) && isempty(c) && isempty(nlfunc) % use lsqnonlin fopts= optimset(optimset('lsqnonlin'),... otbxParam{:}); if ~isempty(LB) || ~isempty(UB) fopts= optimset(fopts,'largescale','on'); end [xf,cost,res,flag]= lsqnonlin(@lsqerrors,x0,LB,UB,fopts,F,c0,m,x,y,sx,varargin{:}); else % use fmincon if ~isempty(nlfunc); nlfunc= @i_nlconstraints; end fopts= optimset(optimset('fmincon'),... otbxParam{:},... 'largescale','off',... 'Algorithm','active-set',... 'MaxSQPIter',1000,... 'HessPattern',ones(size(x0))); [xf,cost,flag]= fmincon(@ssenonlin,x0,A,c,[],[],LB,UB,nlfunc,fopts,F,c0,m,x,y,sx,varargin{:}); end if nargout>4 [res,J,OK]= lsqerrors(xf,F,1,m,x,y,sx,varargin{:}); else OK=1; end if ~isempty(sx) % scale final results xf= sx.*xf; end % final unscaled cost cost= cost*c0^2; OK= all(OK) & flag >= 0; % Sum of squares error function [cost,g]= ssenonlin(B,F,c0,m,x,y,sx,varargin) if nargout==1 res= lsqerrors(B,F,c0,m,x,y,sx,varargin{:}); else [res,J]= lsqerrors(B,F,c0,m,x,y,sx,varargin{:}); g= 2*res'*J; if ~all(isfinite(g)) ind= any(isfinite(J),2) | isfinite(res); g= 2*res(ind)'*J(ind,:); end end cost = res'*res; % residual function function [res,J,ok]= lsqerrors(B,F,c0,m,x,y,sx,varargin) if ~isempty(sx) B= sx.*B; end if nargout==1 res= feval(F.costFunc,B,m,x,y,varargin{:}); else % evaluate residuals and gradients [res,J]= feval(F.costFunc,B,m,x,y,varargin{:}); J= J/c0; if ~isempty(sx) n=size(J,2); J= J*sparse(1:n,1:n,sx,n,n); end end res = res/c0; % deal with nonfinite residuals ok= isfinite(res); if ~all(ok) res(~ok)= get(F,'MaxError'); end % deal with complex residuals if ~isreal(res) ir= abs(imag(res))>1e-6; ok= ok & ir; res(ir)= get(F,'MaxError'); res= real(res); end drawnow % Nonlinear Constraints function [c,ceq,gc,geq]= i_nlconstraints(B,F,~,m,x,y,sx,varargin) if ~isempty(sx) % scale parameters B= sx.*B; end ceq=[]; geq=[]; if nargout>2 [c,gc]=feval(F.Constraints.nlcon,B,m,x,y,varargin{:}); if ~isempty(sx) n = length(sx); gc= sparse(1:n,1:n,1./sx,n,n)*gc; end else c=feval(F.Constraints.nlcon,B,m,x,y,varargin{:}); end