Code covered by the BSD License

# General simulated annealing algorithm

### Joachim Vandekerckhove (view profile)

• 1 file
• 4.15625

28 Mar 2006 (Updated )

Minimizes a function with the method of simulated annealing

anneal(loss, parent, options)
```function [minimum,fval] = anneal(loss, parent, options)
% ANNEAL  Minimizes a function with the method of simulated annealing
% (Kirkpatrick et al., 1983)
%
%  ANNEAL takes three input parameters, in this order:
%
%  LOSS is a function handle (anonymous function or inline) with a loss
%  function, which may be of any type, and needn't be continuous. It does,
%  however, need to return a single value.
%
%  PARENT is a vector with initial guess parameters. You must input an
%  initial guess.
%
%  OPTIONS is a structure with settings for the simulated annealing. If no
%  OPTIONS structure is provided, ANNEAL uses a default structure. OPTIONS
%  can contain any or all of the following fields (missing fields are
%  filled with default values):
%
%       Verbosity: Controls output to the screen.
%                  0 suppresses all output
%                  1 gives final report only [default]
%                  2 gives temperature changes and final report
%       Generator: Generates a new solution from an old one.
%                  Any function handle that takes a solution as input and
%                  gives a valid solution (i.e. some point in the solution
%                  space) as output.
%                  The default function generates a row vector which
%                  slightly differs from the input vector in one element:
%                  @(x) (x+(randperm(length(x))==length(x))*randn/100)
%                  Other examples of possible solution generators:
%                  @(x) (rand(3,1)): Picks a random point in the unit cube
%                  @(x) (ceil([9 5].*rand(2,1))): Picks a point in a 9-by-5
%                                                 discrete grid
%                  Note that if you use the default generator, ANNEAL only
%                  works on row vectors. For loss functions that operate on
%                  column vectors, use this generator instead of the
%                  default:
%                  @(x) (x(:)'+(randperm(length(x))==length(x))*randn/100)'
%        InitTemp: The initial temperature, can be any positive number.
%                  Default is 1.
%        StopTemp: Temperature at which to stop, can be any positive number
%                  smaller than InitTemp.
%                  Default is 1e-8.
%         StopVal: Value at which to stop immediately, can be any output of
%                  LOSS that is sufficiently low for you.
%                  Default is -Inf.
%       CoolSched: Generates a new temperature from the previous one.
%                  Any function handle that takes a scalar as input and
%                  returns a smaller but positive scalar as output.
%                  Default is @(T) (.8*T)
%      MaxConsRej: Maximum number of consecutive rejections, can be any
%                  positive number.
%                  Default is 1000.
%        MaxTries: Maximum number of tries within one temperature, can be
%                  any positive number.
%                  Default is 300.
%      MaxSuccess: Maximum number of successes within one temperature, can
%                  be any positive number.
%                  Default is 20.
%
%
%  Usage:
%     [MINIMUM,FVAL] = ANNEAL(LOSS,NEWSOL,[OPTIONS]);
%          MINIMUM is the solution which generated the smallest encountered
%          value when input into LOSS.
%          FVAL is the value of the LOSS function evaluated at MINIMUM.
%     OPTIONS = ANNEAL();
%          OPTIONS is the default options structure.
%
%
%  Example:
%     The so-called "six-hump camelback" function has several local minima
%     in the range -3<=x<=3 and -2<=y<=2. It has two global minima, namely
%     f(-0.0898,0.7126) = f(0.0898,-0.7126) = -1.0316. We can define and
%     minimise it as follows:
%          camel = @(x,y)(4-2.1*x.^2+x.^4/3).*x.^2+x.*y+4*(y.^2-1).*y.^2;
%          loss = @(p)camel(p(1),p(2));
%          [x f] = ANNEAL(loss,[0 0])
%     We get output:
%               Initial temperature:     	1
%               Final temperature:       	3.21388e-007
%               Consecutive rejections:  	1027
%               Number of function calls:	6220
%               Total final loss:        	-1.03163
%               x =
%                  -0.0899    0.7127
%               f =
%                  -1.0316
%     Which reasonably approximates the analytical global minimum (note
%     that due to randomness, your results will likely not be exactly the
%     same).

%  Reference:
%    Kirkpatrick, S., Gelatt, C.D., & Vecchi, M.P. (1983). Optimization by
%    Simulated Annealing. _Science, 220_, 671-680.

%   joachim.vandekerckhove@psy.kuleuven.be
%   \$Revision: v5 \$  \$Date: 2006/04/26 12:54:04 \$

def = struct(...
'CoolSched',@(T) (.8*T),...
'Generator',@(x) (x+(randperm(length(x))==length(x))*randn/100),...
'InitTemp',1,...
'MaxConsRej',1000,...
'MaxSuccess',20,...
'MaxTries',300,...
'StopTemp',1e-8,...
'StopVal',-Inf,...
'Verbosity',1);

% Check input
if ~nargin %user wants default options, give it and stop
minimum = def;
return
elseif nargin<2, %user gave only objective function, throw error
error('MATLAB:anneal:noParent','You need to input a first guess.');
elseif nargin<3, %user gave no options structure, use default
options=def;
else %user gave all input, check if options structure is complete
if ~isstruct(options)
'Input argument ''options'' is not a structure')
end
fs = {'CoolSched','Generator','InitTemp','MaxConsRej',...
'MaxSuccess','MaxTries','StopTemp','StopVal','Verbosity'};
for nm=1:length(fs)
if ~isfield(options,fs{nm}), options.(fs{nm}) = def.(fs{nm}); end
end
end

% main settings
newsol = options.Generator;      % neighborhood space function
Tinit = options.InitTemp;        % initial temp
minT = options.StopTemp;         % stopping temp
cool = options.CoolSched;        % annealing schedule
minF = options.StopVal;
max_consec_rejections = options.MaxConsRej;
max_try = options.MaxTries;
max_success = options.MaxSuccess;
report = options.Verbosity;
k = 1;                           % boltzmann constant

% counters etc
itry = 0;
success = 0;
finished = 0;
consec = 0;
T = Tinit;
initenergy = loss(parent);
oldenergy = initenergy;
total = 0;
if report==2, fprintf(1,'\n  T = %7.5f, loss = %10.5f\n',T,oldenergy); end

while ~finished;
itry = itry+1; % just an iteration counter
current = parent;

% % Stop / decrement T criteria
if itry >= max_try || success >= max_success;
if T < minT || consec >= max_consec_rejections;
finished = 1;
total = total + itry;
break;
else
T = cool(T);  % decrease T according to cooling schedule
if report==2, % output
fprintf(1,'  T = %7.5f, loss = %10.5f\n',T,oldenergy);
end
total = total + itry;
itry = 1;
success = 1;
end
end

newparam = newsol(current);
newenergy = loss(newparam);

if (newenergy < minF),
parent = newparam;
oldenergy = newenergy;
break
end

if (oldenergy-newenergy > 1e-6)
parent = newparam;
oldenergy = newenergy;
success = success+1;
consec = 0;
else
if (rand < exp( (oldenergy-newenergy)/(k*T) ));
parent = newparam;
oldenergy = newenergy;
success = success+1;
else
consec = consec+1;
end
end
end

minimum = parent;
fval = oldenergy;

if report;
fprintf(1, '\n  Initial temperature:     \t%g\n', Tinit);
fprintf(1, '  Final temperature:       \t%g\n', T);
fprintf(1, '  Consecutive rejections:  \t%i\n', consec);
fprintf(1, '  Number of function calls:\t%i\n', total);
fprintf(1, '  Total final loss:        \t%g\n', fval);
end

```