lst_sq.m

用MATLAB实现遗传算法中的好多问题

%function [x,OPTIONS] = lst_sq(FUN,x,OPTIONS,GRADFUN,time_out,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,P18,P19,P20,P21)
% See leastsq.m for usage. Modified to time-out after "time_out" seconds.
%
%       X=LST_SQ('FUN',X,OPTIONS,GRADFUN,time_out,P1,P2,..) allows
%       coefficients, P1, P2, ... to be passed directly to FUN:
%       F=FUN(X,P1,P2,...). Empty arguments ([]) are ignored.

% Author             Date     Predecessor   Modification
% ======             ====     ===========   ============
% B.McKay M.Willis   25/5/95  leastsq.m     use etime instead of tic/toc
% J. Elsey           31/7/95  lst_sq.m      increase the number of allowable inputs
% B.McKay	     6/10/95     -          returns latest x in case of time out (rather than initial)
% 
% Function Calls:   
% ===============    
%
% Last Modification: 6/10/95
% =================
%

function [x,OPTIONS] = lst_sq(FUN,x,OPTIONS,GRADFUN,time_out,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,P18,P19,P20,P21)

global optimea optimeb optimec;

% ------------Initialization----------------
t1=clock;

x_initial = x;
XOUT=x(:);
[nvars]=length(XOUT);

evalstr = [FUN];
if ~any(FUN<48)
        evalstr=[evalstr, '(x'];
        
        optimes=clock;
        
        for i=1:nargin - 4
                evalstr = [evalstr,',P',num2str(i)];
        end
        
        optimef=clock;
        optimea=optimea+etime(optimef,optimes);
        evalstr = [evalstr, ')'];
end


if nargin < 3, OPTIONS=[]; end
if nargin < 4, GRADFUN=[]; end

if length(GRADFUN)
        evalstr2 = [GRADFUN];
        if ~any(GRADFUN<48)
                evalstr2 = [evalstr2, '(x'];
                
                optimes=clock;
                
                for i=1:nargin - 4
                    evalstr2 = [evalstr2,',P',num2str(i)];
                end
                
                optimef=clock;
                optimeb=optimeb+etime(optimef,optimes);
                
                evalstr2 = [evalstr2, ')'];
        end
end

f = eval(evalstr);
f=f(:);
nfun=length(f);
GRAD=zeros(length(XOUT),nfun);

OLDX=XOUT;
MATX=zeros(3,1);
MATL=[f'*f;0;0];
OLDF=f'*f;
FIRSTF=f'*f;
[OLDX,OLDF,OPTIONS]=lsint(XOUT,f,OPTIONS);
EstSum=0.5;
GradFactor=0; 

CHG = 1e-7*abs(XOUT)+1e-7*ones(nvars,1);

OPTIONS(10)=1;
status=-1;

while status~=1
if etime(clock,t1) > time_out; x(:)=OLDX; return; end

% Work Out Gradients
        if ~length(GRADFUN) | OPTIONS(9)
                OLDF=f;
                
                CHG = sign(CHG+eps).*min(max(abs(CHG),OPTIONS(16)),OPTIONS(17));
                
                optimes=clock;
                
                for gcnt=1:nvars
                    if XOUT==[] | any(isnan(XOUT(:)))==1 | any(XOUT(:)==inf)==1 | any(...
                       XOUT(:)==-inf)==1,
                       x=x_initial;return, end;
                    
                    temp = XOUT(gcnt);
                    XOUT(gcnt) = temp +CHG(gcnt);
                    x(:) = XOUT;
                    f(:) = eval(evalstr);
                    GRAD(gcnt,:)=(f-OLDF)'/(CHG(gcnt));
                    XOUT(gcnt) = temp;
                end
                
                optimef=clock;
                optimec=optimec+etime(optimef,optimes);
                
                f = OLDF;
                OPTIONS(10)=OPTIONS(10)+nvars;
% Gradient check
                if OPTIONS(9) == 1
                    GRADFD = GRAD;
                    x(:)=XOUT; GRAD = eval(evalstr2);
                    graderr(GRADFD, GRAD, evalstr2);
                    OPTIONS(9) = 0;
                end
        else
                x(:) = XOUT;
                OPTIONS(11)=OPTIONS(11)+1;
                GRAD = eval(evalstr2);
        end
        % Try to set difference to 1e-8 for next iteration
        if nfun==1
                CHG = nfun*1e-8./GRAD;
        else
                CHG = nfun*1e-8./sum(abs(GRAD)')';
        end

        gradf = 2*GRAD*f;
        fnew = f'*f;
%---------------Initialization of Search Direction------------------
        if status==-1
                if GRAD==[] | any(isnan(GRAD(:)))==1 | any(GRAD(:)==inf)==1 | any(...
                   GRAD(:)==-inf)==1,
                x=x_initial;return, end;
                
                if cond(GRAD)>1e8
                    SD=-(GRAD*GRAD'+(norm(GRAD)+1)*(eye(nvars)))\(GRAD*f);
                    if OPTIONS(5)==0, GradFactor=norm(GRAD)+1; end
                    how='COND';
                else
%               SD=GRAD'\(GRAD'*X-F)-X;
                    SD=-(GRAD*GRAD'+GradFactor*(eye(nvars)))\(GRAD*f);
                end
                FIRSTF=fnew;
                OLDG=GRAD;
                GDOLD=gradf'*SD;
                OPTIONS(18)=1;
                if OPTIONS(1)>0
                    disp([sprintf('%5.0f %12.3g %12.3g ',OPTIONS(10),fnew,OPTIONS(18)),sprintf('%12.3g  ',GDOLD)]);
                end
                XOUT=XOUT+OPTIONS(18)*SD;
                if OPTIONS(5)==0
                    newf=GRAD'*SD+f;
                    GradFactor=newf'*newf;
                    SD=-(GRAD*GRAD'+GradFactor*(eye(nvars)))\(GRAD*f); 
                end
                newf=GRAD'*SD+f;
                XOUT=XOUT+OPTIONS(18)*SD;
                EstSum=newf'*newf;
                status=0;
                if OPTIONS(7)==0; PCNT=1; end
                
        else
%-------------Direction Update------------------
                gdnew=gradf'*SD;
                if OPTIONS(1)>0, 
                    num=[sprintf('%5.0f %12.3g %12.3g ',OPTIONS(10),fnew,OPTIONS(18)),sprintf('%12.3g  ',gdnew)];
                end
                if gdnew>0 & fnew>FIRSTF

% Case 1: New function is bigger than last and gradient w.r.t. SD -ve
% ... interpolate. 
                    how='inter';
                    [stepsize]=cubici1(fnew,FIRSTF,gdnew,GDOLD,OPTIONS(18));
                    OPTIONS(18)=0.9*stepsize;
                elseif fnew<FIRSTF

%  New function less than old fun. and OK for updating 
%         .... update and calculate new direction. 
                    [newstep,fbest] =cubici3(fnew,FIRSTF,gdnew,GDOLD,OPTIONS(18));
                    if fbest>fnew,fbest=0.9*fnew; end 
                    if gdnew<0
                        how='incstep';
                        if newstep<OPTIONS(18),  newstep=(2*OPTIONS(18)+1e-4); how=[how,'IF']; end
                        OPTIONS(18)=abs(newstep);
                    else 
                        if OPTIONS(18)>0.9
                            how='int_step';
                            OPTIONS(18)=min([1,abs(newstep)]);
                        end
                    end
% SET DIRECTION.
% Gauss-Newton Method    
                    temp=1;
                    if OPTIONS(5)==1 
                        if OPTIONS(18)>1e-8 & cond(GRAD)<1e8
                            SD=GRAD'\(GRAD'*XOUT-f)-XOUT;
                            if SD'*gradf>eps,how='ERROR- GN not descent direction',  end
                            temp=0;
                        else
                            if OPTIONS(1) > 0
                                disp('Conditioning of Gradient Poor - Switching To LM method')
                            end
                            how='CHG2LM';
                            OPTIONS(5)=0;
                            OPTIONS(18)=abs(OPTIONS(18));                               
                        end
                    end
                    
                    if (temp)      
% Levenberg_marquardt Method N.B. EstSum is the estimated sum of squares.
%                                 GradFactor is the value of lambda.
% Estimated Residual:
                        if EstSum>fbest
                            GradFactor=GradFactor/(1+OPTIONS(18));
                        else
                            GradFactor=GradFactor+(fbest-EstSum)/(OPTIONS(18)+eps);
                        end
                        SD=-(GRAD*GRAD'+GradFactor*(eye(nvars)))\(GRAD*f); 
                        OPTIONS(18)=1; 
                        estf=GRAD'*SD+f;
                        EstSum=estf'*estf;
                        if OPTIONS(1)>0, num=[num,sprintf('%12.3g ',GradFactor)]; end
                    end
                    gdnew=gradf'*SD;

                    OLDX=XOUT;
% Save Variables
                    FIRSTF=fnew;
                    OLDG=gradf;
                    GDOLD=gdnew;        

                    % If quadratic interpolation set PCNT
                    if OPTIONS(7)==0, PCNT=1; MATX=zeros(3,1);  MATL(1)=fnew; end
                else 
% Halve Step-length
                    how='Red_Step';
                    if fnew==FIRSTF,
                        if OPTIONS(1)>0,
                            disp('No improvement in search direction: Terminating')
                        end
                        status=1;
                    else
                        OPTIONS(18)=OPTIONS(18)/8;
                        if OPTIONS(18)<1e-8
                            OPTIONS(18)=-OPTIONS(18);
                        end
                    end
                end
                XOUT=OLDX+OPTIONS(18)*SD;
                if OPTIONS(1)>0,disp([num,how]),end

        end %----------End of Direction Update-------------------
        if OPTIONS(7)==0, PCNT=1; MATX=zeros(3,1);  MATL(1)=fnew; end
% Check Termination 
        if max(abs(SD))< OPTIONS(2) & (gradf'*SD) < OPTIONS(3) & max(abs(gradf)) < 10*(OPTIONS(3)+OPTIONS(2))
                if OPTIONS(1) > 0
                    disp('Optimization Terminated Successfully')  
                end
                status=1; 
        elseif OPTIONS(10)>OPTIONS(14)
                disp('maximum number of iterations has been exceeded');
                if OPTIONS(1)>0
                    disp('Increase OPTIONS(14)')
                end
                status=1;
        else

% Line search using mixed polynomial interpolation and extrapolation.
                if PCNT~=0
                    while PCNT > 0
                    if etime(clock,t1) > time_out; x(:)=OLDX;; return; end
                    
                        x(:) = XOUT; f(:)  = eval(evalstr); OPTIONS(10)=OPTIONS(10)+1;
                        fnew = f'*f;
                        if fnew<OLDF'*OLDF, OX = XOUT; OLDF=f; end
                        [PCNT,MATL,MATX,steplen,fnew,how]=searchq(PCNT,fnew,OLDX,MATL,MATX,SD,GDOLD,OPTIONS(18),how);
                        OPTIONS(18)=steplen;
                        XOUT=OLDX+steplen*SD;
                        if fnew==FIRSTF,  PCNT=0; end
                    end
                    XOUT = OX;
                    f=OLDF;
                else
                    x(:)=XOUT; f(:) = eval(evalstr); OPTIONS(10)=OPTIONS(10)+1;
                end
        end
end
OPTIONS(8) = fnew;
XOUT=OLDX;
x(:)=XOUT;