anneal.m

matlab算法集 matlab算法集

function [x,ev,j] = anneal (x0,f0,f1,tol,eta,seed,gamma,mu,m,a,b,f)
%-----------------------------------------------------------------------
% Usage:       [x,ev,j] = anneal (x0,f0,f1,tol,eta,seed,gamma,mu,m,a,b,f);
%
% Description: Use the simulated annealing method to solve the following
%              n-dimensional constrained optimization problem:
%
%                  minimize:   f(x)
%                  subject to: a <= x <= b 
%         
%              The function anneal first performs a global search using 
%              a simulated annealing method, and then performs a local 
%              search using the function dfp.
%
% Inputs:      x0    = n by 1 vector containing initial guess
%              f0    = lower bound on f(x)
%              f1    = upper bound on f(x) (f1 > f0)
%              tol   = global error tolerance (tol >= 0)
%              eta   = local error tolerance (eta >= 0)
%              seed  = integer used to initialize random number
%                      generator (for a random seed use seed <= 0)  
%              gamma = localizaton parameter (0 < gamma < 1)  
%              mu    = penalty parameter (mu >= 0)
%              m     = maximum number of iterations (m >= 1)
%              a     = n by 1 vector of lower bounds on x
%              b     = n by 1 vector of upper bounds on x (b[k] >= a[k])
%              f     = string containing name of objective function: 
%                      f(x).  The form of f(x) is:
%
%                      function y = f(x)
%
%                      When f is called with n by 1 vector x, it must 
%                      return the value of the objective function f(x).
%
% Outputs:     x  = n by 1 solution vector
%              ev = number of scalar function evaluations
%              k  = number of iterations performed. Iterations for the
%                   initial global search continue until (f(x) - f0) <
%                   tol*(f1 - f0) or m iterations are performed.  
%                   Iterations for the final local search continue until
%                   m additional iterations are performed or the 
%                   following criterion is satisfied where h is the 
%                   step length and em is the machine epsilon:
%
%                      (||df(x)/dx|| < eta) or (h*||x|| < em)
%
%                   Normally, anneal is called with:
%
%                      0 < eta << tol << 1. */
%-----------------------------------------------------------------------
   
   global a_limit b_limit

% Initialize

   n = length(x0);
   chkvec (x0,1,'anneal');
   f1    = args (f1,f0,f1,3,'anneal');
   tol   = args (tol,0,tol,4,'anneal');
   eta   = args (eta,0,eta,5,'anneal');
   gamma = args (gamma,0,1,6,'anneal');
   mu    = args (mu,0,mu,7,'anneal');
   m     = args (m,1,m,8,'anneal');

   randinit (seed);
   x = x0;
   for i = 1 : n
      b(i) = args (b(i),a(i),b(i),10,'anneal');
   end
   a_limit = a;
   b_limit = b;
   eval = 0;                              % function evaluations 
   j = 0;                                 % iterations 
   g1 = zeros (n,1);                      % gradient
   x1 = zeros (n,1);                      % new x
   x2 = zeros (n,1);                      % optimal x so far
   alpha = log(2)/gamma;
   F0 = feval(f,x);
   eval = eval + 1;
   x2 = x;
   F2 = F0;
   Df = F0 - f0;
   d = max(b-a);
   delta = tol + 1;

% Find x

   hwbar = waitbar(0,'Computing Optimum: anneal');
   while (delta> tol) & (j < m)
      
% Generate a random step 

      waitbar (max(j/m,tol/delta))
      s = gamma*d*Df/(f1 - f0);
      for i = 1 : n
         x1(i) = x(i) + randg (1,1,0,s);
         x1(i) = min(b(i),max(a(i),x1(i)));
      end
      F1 = feval(f,x1);
      eval = eval + 1; 
      df = F1 - F0;

% Should we use it? 

      if df < 0                         % a descent
         x = x1;
         F0 = F1;
         if F1 < F2
            x2 = x1;
            F2 = F1;
         end
         j = j + 1;
      else
         p = exp(-alpha*df/Df);
         z = randu (1,1,0,1);
         if p > z
            x = x1;
            F0 = F1;
            j = j + 1;
         end
      end

% Check termination criterion 

      Df = F0 - f0;
      delta = Df/(f1 - f0);
   end
   close(hwbar)

% Perform a local search 

   x1 = x2;
   [x,ev,k] = penalty (x1,mu,tol,m,f,'','funannea');
   ev = ev + eval;
   j = j + k;