bracket.m

matlab算法集 matlab算法集

function [a,b,c,err,ev] = bracket (h,x,d,mu,f,p,q,dm)
%-----------------------------------------------------------------------
% Usage:       [a,b,c,err,k] = bracket (h,x,d,mu,f,p,q,dm)
%
% Description:  Construct a three-point pattern {a,b,c} such that
%               a < b < c, F(b) < F(a), and F(b) <= F(c) where:
%
%                   F(y) = f(x+y*d) + mu*(P(x+y*d) + Q(x+y*d))
%
%                The interval [a,c] then contains the minimum of a
%                unimodal function F(y).  Here P(x) is a penalty
%                function associated with the r equality constraints,
%                p(x) = 0, and Q(x) is a penalty function associated
%                with the s inequality constraints, q(x) >= 0.        
%
% Inputs:        h  = initial estimate of step size (h > 0)
%                x  = n by 1 start vector
%                d  = n by 1 search direction vector
%                mu = penalty paramter (mu >= 0)
%                f  = name of objective function: minimize: f(x)
%                p  = name of equality constraint function: p(x) = 0
%                q  = name of inequality constraint function: q(x) >= 0
%
%                     The forms of f, p, and q are:
%
%                     function y = f(x)
%                     function u = p(x)
%                     function v = q(x)
%
%                     When f is called with n by 1 vector x, it must
%                     return the scalar y = f(x).  When p is called with 
%                     n by 1 vector x,it must compute r by 1 vector 
%                     u = p(x). When q is called with n by 1 vector x, 
%                     it must compute s by 1 vector v = q(x).
%               
%                dm = optional display mode.  If present,
%                     intermediate results are displayed.
% 
% Outputs	 a   = lower limit of three-point pattern
%                b   = interior point of three-point pattern
%                c   = upper limit of three-point pattern
%                err = an error code.  If *err = 0, the search
%                      was successful.  Otherwise, a descent
%                      direction was NOT found using the minimum 
%                      step length.
%                ev  = number of scalar function evaluations 
%-----------------------------------------------------------------------

% Initialize

   h  = args (h,eps,h,1,'bracket');   		
   chkvec (x,2,'bracket');
   chkvec (d,3,'bracket');
   mu = args (mu,0,mu,4,'bracket');
 
   display = nargin > 7;
   ev = 2;		                  
   n = length (x);
   u = zeros(n,1);
   hmin = eps*max([eps,norm(x,inf)]);
   hmax = 1/(eps*eps);
   if display
      fprintf ('\nbracket range: %g to %g',hmin,hmax);
   end
   err = 0;
   a = 0;
   b = h;
   c = 2*h;
   F0 = fmu (f,p,q,x,mu);
   u = x + b*d;
   F1 = fmu (f,p,q,u,mu);
   
% Find b by halving candidate interval 

   while (F1 >= F0) & (b > hmin) 
      c = b;
      b = b/2; 	
      F2 = F1;
      u = x + b*d;  
      F1 = fmu(f,p,q,u,mu);
      ev = ev  + 1;
      if display
         fprintf ('\n f(%g) = %g, f(%g) = %g, f(%g) = %g',...
                  a,F0,b,F1,c,F2);  
      end
   end

% Check for problems finding b 

   if F1 >= F0 
      err = 1;
      if display	
         fprintf ('\nUnable to find a descent in bracket using a step');
         fprintf ('\nof length %g. Initial step length = %g. ',b,h);
         fprintf ('\nf(%g) = %g, f(%g) = %g, f(%g) = %g',...
                   a,F0,b,F1,c,F2);  
         wait
      end   
      return 
   end
      
% Find c by doubling candidate interval 

   if ev <= 2 
      u = x + c*d;  
      F2 = fmu (f,p,q,u,mu);
      ev = ev + 1;
      if display
         fprintf ('\n  f(%g) = %g, f(%g) = %g, f(%g) = %g',...
                   a,F0,b,F1,c,F2);  
      end

      while (F2 < F1) & (c < hmax)
         a = b;
         b = c;
         c = 2*c;
         F0 = F1;
         F1 = F2;
         u = x + c*d;  
         F2 = fmu (f,p,q,u,mu);
         ev = ev + 1;
         if display
            fprintf ('\n  f(%g) = %g, f(%g) = %g, f(%g) = %g',...
                      a,F0,b,F1,c,F2);  
         end
      end
   end
      
% Check for problems finding c 

   if display
      if (F2 < F1) 
         fprintf ('\nUnable to reach minimum in bracket using a step');
         fprintf ('\nof length %g. Initial step length = %g. ',c,h);
         fprintf ('\nf(%g) = %g, f(%g) = %g, f(%g) = %g',...
                   a,F0,b,F1,c,F2);  
         wait 
       end
   end         

% Finalize 

   if display
      wait
   end
%-----------------------------------------------------------------------