simulated_annealing 0.1.4-nullsafety

Simulated annealing framework for Dart. Enables quickly setting up a simulator for finding the global minimum of multi-variate functions.

Simulated Annealing For Dart

Introduction

Simulated annealing (SA) is an algorithm aimed at finding the global minimum of a function E(x₀, x₁, ..., x_n) for a given region ω(x₀, x₁, ..., x_n). The function to be minimized can be interpreted as the system energy. In that case, the global minimum represents the ground state of the system.

The algorithm name was coined by Kirkpatrick et al. [1] and was derived from the process of annealing a metal alloy or glass. The first step of the annealing process consists of heating a solid material above a critical temperature. This allows its atoms to gain sufficient kinetic energy to be able to rearrange themselves. Then the temperature is decreased sufficiently slowly in order to minimize atomic lattice defects as the material solidifies.

Simulated annealing works by randomly selecting a new point in the neighbourhood of the current solution, evaluating the energy function, and deciding if the new solution is accepted or rejected.

If for a newly selected point the energy E is lower that the previous minimum energy E_min, the new solution is accepted: P(ΔE < 0, T) = 1, where ΔE = E - E_min.

Crucially, if ΔE > 0, the algorithm still accepts the new solution with probability: P(ΔE > 0, T) = e^{-ΔE/(k_B·T)}. Accepting up-hill moves provides a method of escaping from local energy minima.

The process is demonstrated in the animation above. The left figure shows a spherical 3D search space while the energy value is represented by colour. The figure on the right shows a projection of the energy function onto the x-y plane. Initially, random points are chosen from a large region encompasing the entire spherical search space. In the simulation shown above, intermediate solutions near the local minimum are followed by up-hill moves. As the temperature drops the search neighourhood is contracted and the solution converges to the global minimum.

Usage

To use this package include simulated_annealing as a dependency in your pubspec.yaml file.

The following steps are required to set up the SA algorithm.

1. Specify the search space ω.
2. Define an annealing schedule.
3. Define the system energy field.
4. Extend the class Simulator implementing the methods prepareLog() and recordLog() or create an instance of LoggingSimulator.
5. Start the simulated annealing process.

Click to show source code.

import 'dart:io';
import 'dart:math';

import 'package:list_operators/list_operators.dart';
import 'package:simulated_annealing/simulated_annealing.dart';

void main() async {

  // Defining a spherical space.
final radius = 2;
final x = FixedInterval(-radius, radius);
final y = ParametricInterval(
  () => -sqrt(pow(radius, 2) - pow(x.next(), 2)),
  () => sqrt(pow(radius, 2) - pow(x.next(), 2)),
);
final z = ParametricInterval(
  () => -sqrt(pow(radius, 2) - pow(y.next(), 2) - pow(x.next(), 2)),
  () => sqrt(pow(radius, 2) - pow(y.next(), 2) - pow(x.next(), 2)),
);
final dxMin = <num>[1e-6, 1e-6, 1e-6];
final space = SearchSpace([x, y, z], dxMin: [1e-6, 1e-6, 1e-6]);

// Defining an energy function.
final xGlobalMin = [0.5, 0.7, 0.8];
final xLocalMin = [-1.0, -1.0, -0.5];
num energy(List<num> x) {
  return 4.0 -
      4.0 * exp(-4 * xGlobalMin.distance(x)) -
      2.0 * exp(-6 * xLocalMin.distance(x));
}

// Constructing an instance of `EnergyField`.
final energyField = EnergyField(
  energy,
  space,
);
  // Constructing an instance of `LoggingSimulator`.
  final simulator = LoggingSimulator(energyField, exponentialSequence,
      iterations: 750, gammaStart: 0.7, gammaEnd: 0.05);

  print(await simulator.info);

  final xSol = await simulator.anneal((_) => 1, isRecursive: true);
  await File('../data/log.dat').writeAsString(simulator.rec.export());

  print('Solution: $xSol');
}

Estimating the Value of the System Boltzmann Constant

The Boltzmann constant k_B relates the system temperature with the kinetic energy of particles in a gas.

In the context of SA, k_B relates the system temperature with the probability of accepting a solution:

P(ΔE > 0, T) = e^{-ΔE/(k_B·T)}    P(ΔE <0, T) = 1.0.

The expression above ensures that the acceptance probability decreases with decreasing temperature (for ΔE > 0). As such, the temperature is a parameter that controls the probability of up-hill moves.

Most authors set k_B = 1 and scale the temperature to control the solution acceptance probability. I find it more practical to use an independent temperature scale with the lowest value T_end and calculate the system dependent constant k_B (see section Algorithm Tuning).

An estimate for the average scale of the variation of the energy function ΔE can be obtained by sampling the energy function E at random points in the search space ω and calculating the sample standard deviation σ_E [3].

For continuous problems, the size of the search region around the current solution is gradually contracted to ω_end in order to generate a solution with the required precision.

The constant k_B is fixed such that: P(ΔE_end, T_end) = e^{-ΔE_end/(k_B·T_end)} = γ_end, where T_end is the final annealing temperature.

The initial temperature is then set such that the initial acceptance probability is: P(ΔE_start,T_start) = e^{-ΔE_start/(k_B·T_start)} = γ_start.

Algorithm Tuning

For discrete problems it can be shown that by selecting a sufficiently high initial factor k_B·T the algorithm converges to the global minimum if temperture decreases on a logarithmic scale (slow cooling schedule) and the number of inner iterations (Markov chain length) is sufficiently high [2].

Practical implementations of the SA algorithm aim to generate an acceptable solution with minimal computational effort. For such fast cooling schedules, algorithm convergence to the global minimum is not strictly guaranteed. In that sense, SA is a heuristic approach and some degree of trial and error is required to determine which annealing schedule works best for a given problem.

The behaviour of the annealing simulator can be tuned using the following optional parameters of the class Simulator:

• tEnd: The final annealing temperature with default value 1e-4. It is arbitrary as it is scaled with k_B.
• gammaStart: Initial acceptance probability with default value 0.7. Useful values for γ_start are in the range of (0.7, 0.9). If γ_start is too low, up-hill moves are unlikely (potentially) preventing the SA algorithm from escaping a local miniumum. If γ_start is set close to 1.0 the algorithm will accept too many up-hill moves at high temperatures wasting computational time and delaying convergence.
• gammaEnd: Final acceptance probability. Towards the end of the annealing process one assumes that the solution has converged towards the global minimum and up-hill moves should be restricted. For this reason γ_end has default value 0.05.
• deltaEnergyStart: A critical SA parameter used to estimate the initial temperature. If ΔE_start is too large the algorithm will oscillate wildy between random points and will most likely not converge towards an acceptable solution. On the other hand, if ΔE_start is too small up-hill moves are unlikely and the solution most likely converges towards a local minimum or a point situated in a plateau-shaped region.
• deltaEnergyEnd: Typical energy variation ΔE if the current position is perturbed within the minimum search neighbourhood ω_end. It is used to calculate k_B.
• iterations: Determines the number of temperature steps in the annealing schedule.

Selecting an annealing schedule.

It is recommended to start with a higher number of outer iterations (number of entries in the sequence of temperatures) and log quantities like the current system energy, temperature, and the intermediate solutions.

The figure below shows a typical SA log where the x-coordinate of the solution (green dots) converges asymptotically to 0.5. The graph is discussed in more detail here.

The number of inner iterations (performed while the temperature is kept constant) is also referred to as Markov chain length and is determined by a function with typedef MarkovChainLength see method anneal.

For fast cooling schedules convergence to an acceptable solution can be improved by increasing the number of inner iterations.

Examples

Further information can be found in the folder example. The following topics are covered:

• search space,
• annealing schedule,
• system energy and logging simulator.

Features and bugs

Please file feature requests and bugs at the issue tracker.