Simulated Annealing the Swiss Army Knife of Global Optimization


Google+ Facebook Twitter LinkedIn Dzone Reddit Digg Blogger Hacker News Addthis

By Markus SprunckRevision: 1.0; Status: final; Last Content Change: Jul 13, 2013;
This article describes the use of Simulated Annealing for global optimization problems using the example of JavaThe application reads a caller-callee-graph, some initial cluster information and makes a cluster analysis. 

Some years ago I wrote the core functionality of this application in C++ and re-implemented it in pure Java, with some improvements. Actually, the new Java implementation is better in terms of speed than the old C++ implementation. The reason for this significant performance improvement is that now the global cost function is equated just as delta and not always in total.

What is Simulated Annealing?

In Wikipedia you may read the following description:

"Simulated annealing (SA) is a generic probabilistic metaheuristic for the global optimization problem of
 locating a good approximation to the global optimum of a given function in a large search space. [...]
 For certain problems, simulated annealing may be more efficient than exhaustive enumeration — provided
 that the goal is merely to find an acceptably good solution in a fixed amount of time, rather than the best
 possible solution." [1]


Simulated Annealing is no silver bullet and you can't solve all problems with this algorithm, but is a good choice in cases you have no reliable standard algorithm with good performance and an approximation is enough. Simulated Annealing is more or less a Swiss Army Knife that should be in your pocket.

How Simulated Annealing works?

The term Simulated Annealing comes from annealing of metal in physics. If you take a piece of metal (e.g. copper) and heat to very high temperatures, the diffusion of atoms within a solid material will be more likely. 

This diffusion is a stochastic process depending on the average kinetic energy of the atoms. The diffusion of atoms changes the rate of defects in the crystal structure. Sometimes the energy of the crystal is increases and sometimes decreased. A slow cooling will repair a lot of the misplaced atoms, as long the cooling is slow. As a result the metal will be softer than before heating. The global energy in the crystal is minimized and in principle the same happens in the simulated annealing algorithm.

In Simulated Annealing there are random changes to a model function which reduce or increase a global cost function. Starting with a high temperatures the likelihood of steps with increasing cost function are also high. With decreasing temperature also the likelihood of a step that increases the cost function will be smaller. 

At the end of the optimization just steps that reduce the cost function will be allowed (hill climber algorithm).
The slow cooling helps avoid to get trapped in local optimums of the cost function. 

Example Application for Path Optimization and Cluster Analysis

In Figure 1 the output of a small test graph is shown. 

Figure 1: Cluster Analysis with Simulated Annealing

The left diagram in the screen shot is the caller-callee-graph as imported from a file. The nodes are in more or less random order. In the middle you see the re-arranged graph, for this simulated annealing is used to minimize the path length of the links. Here the global length of all the paths between the nodes are minimized like in Traveling Salesman Problem [2].

This is done in the following method: 

private void evaluateNextStepPosition() {
// current value of the cost function
double currentCost = getGlobalCost(costFactor);
double bestCost = currentCost;
status.reset();
status.setTemperature(getCurrentTemperature());
status.incrementNumber();

for (long index = 1; index <= iterations; index++) {

final double currentCostOLD = currentCost;
currentCost = evaluateNodePositionCandidate(currentCost);
if (currentCost > bestCost) {
final double temp = (bestCost - currentCost) / status.getTemperature();
if (temp > -15.0 && Math.exp(temp) > Math.random()) {
bestCost = currentCost;
status.incrementWorse();
} else {
swapNodes();
status.incrementRejected();
currentCost = currentCostOLD;
}
} else if (currentCost < bestCost) {
bestCost = currentCost;
status.incrementBetter();
} else {
status.incrementConst();
}
}
}

The right diagram in figure 1 shows the nodes after the cluster analysis and with a re-arranged node distribution. The cluster analysis is done in the following method:

private void evaluateNextStepCluster() {
// current value of the cost function
double currentCost = getGlobalCost(costFactor);
double bestCost = currentCost;
status.reset();
status.setTemperature(getCurrentTemperature());
status.incrementNumber();

for (long index = 1; index <= iterations; index++) {
final double currentCostOLD = currentCost;
currentCost = evaluateNodeClusterCandidate(currentCost);
if (currentCost > bestCost) {
final double temp = (bestCost - currentCost) / status.getTemperature();
if (temp > -15.0 && Math.exp(temp) > Math.random()) {
bestCost = currentCost;
status.incrementWorse();
} else {
firstNode.setCluster(firstCluster);
status.incrementRejected();
currentCost = currentCostOLD;
}
} else if (currentCost < bestCost) {
bestCost = currentCost;
status.incrementBetter();
} else {
status.incrementConst();
}
}
}

The basic algorithm for both completly differnt tasks is alway the same. For a number of iterations there happens a random change in the model. 

    currentCost = evaluateNodeClusterCandidate(currentCost);

Then the new cost function is equated. Depending on the temperature the new change is accepted even if the
currentCost is larger than bestCost. If the cost function doesn't change nothing happens. The as the temperature decreases the likelihood of decreasing steps gets smaller because of the term

    (temp > -15.0) && (Math.exp(temp) > Math.random())

In file 1 you may see the complete Model class and on GitHub (see link below) you can find the complete Eclipse project with all supporting and IO classes.

// File #1: Model.java - Main Class with Simulated Annealing Aglorithm 
package com.sw_engineering_candies.yasa.model;

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

import org.apache.log4j.Logger;

import com.sw_engineering_candies.yasa.io.Status;

public final class Model implements Parameter {

public static final String DEFAULT_CLUSTER_NAME = "-";

private static final String STRING_FORMAT_LONG = "%10d";

private static final String STRING_FORMAT_FLOAT = "%f";

/** standard logger (see log4j.properties file for details) */
private static final Logger LOGGER = Logger.getLogger(Model.class);

/** list of all nodes in the model; used for random access by number */
public static final int DEFAULT_SIZE_NODE_NUMBER = 2500;
private final List<Node> nodes = new ArrayList<Node>(DEFAULT_SIZE_NODE_NUMBER);

/** list of all links in the model */
private static final int DEFAULT_SIZE_LINK_NUMBER = 10000;
private final List<Link> links = new ArrayList<Link>(DEFAULT_SIZE_LINK_NUMBER);

/** list of all clusters in the model; used for random access by name */
public static final int DEFAULT_SIZE_CLUSTER_NUMBER = 100;
private final List<Cluster> clusters = new ArrayList<Cluster>(DEFAULT_SIZE_CLUSTER_NUMBER);

/** stores the current status during the optimization */
private final Status status = new Status();

/** dependent from the size of the model and the type of optimization */
private long iterations = 0L;

/** steps of optimization, for each step the temperature will be changed */
private long steps = 0L;

/** factor for the evaluation of the temperature */
private double decay = 0.0d;

/** start temperature */
private double temperature = 0.0d;

/** factor to scale the cost functions to 100% */
private double costFactor = 1.0d;

/** maximal number of columns of the grid */
private int columns = 0;

/** store the old cluster during an iteration to undo the change */
private Cluster firstCluster, oldCluster = null;

/** store the first node candidate during an iteration to undo the change */
private Node firstNode = null;

/** store the second node candidate during an iteration to undo the change */
private Node secondNode = null;

/** defines the type of the cost function */
private boolean optimizePosition = true;

private double getGlobalCost(final double d) {
long totalCosts = 0;
if (optimizePosition) {
for (final Node node : nodes) {
totalCosts += getPositionDeltaCost(node);
}
} else {
for (final Node node : nodes) {
if (!node.isClusterNode()) {
totalCosts += getClusterDeltaCost(node);
}
}
}
return totalCosts / d;
}

public long getClusterDeltaCost(final Node node) {
long totalCosts = 0;
for (final Link link : node.getLinks()) {
totalCosts += getLinkCost(link, node);
}
return totalCosts;
}

public long getPositionDeltaCost(final Node node) {
long cost = 0L;
for (final Link link : node.getLinks()) {
final Node target = link.getCallee();
final long target1 = (node.getColumn()                      - target.getColumn()) * (node.getColumn()                             - target.getColumn()) + (node.getRow() - target.getRow())
    * (node.getRow() - target.getRow());
cost += target1 == 1 ? 0 : link.isClusterLink() ? target1 : 2 * target1;

final Node source = link.getCaller();
final long source1 = (node.getColumn()                             - source.getColumn()) * (node.getColumn()                             - source.getColumn()) + (node.getRow() - source.getRow())
        * (node.getRow() - source.getRow());
cost += source1 == 1 ? 0 : link.isClusterLink() ? source1 : 2 * source1;
}
return cost;
}

private static long getLinkCost(final Link link, final Node node) {

if (link.getCaller() == node) {
return node.getCluster() != link.getCallee().getCluster() ? COST_FACTOR_LINKS : 0;
} else {
return node.getCluster() != link.getCaller().getCluster() ? COST_FACTOR_LINKS : 0;
}
}

private double evaluateNodeClusterCandidate(final double currentCost) {
do {
final int index = (int) Math.round(Math.floor(Math.random() * nodes.size()));
firstNode = nodes.get(index);
} while (!firstNode.isFrozen());

firstCluster = firstNode.getCluster();
final double oldDeltaCostFunction = getClusterDeltaCost(firstNode);

// create the new state
firstNode.setCluster(evaluateRandomCluster(firstCluster));

// return the cost function
final double newDeltaCostFunction = getClusterDeltaCost(firstNode);
return currentCost + (newDeltaCostFunction - oldDeltaCostFunction) / costFactor;
}

private double evaluateNodePositionCandidate(final double currentCost) {
// remember the old state
int index = (int) Math.round(Math.floor(Math.random() * nodes.size()));
firstNode = nodes.get(index);

index = (int) Math.round(Math.floor(Math.random() * nodes.size()));
secondNode = nodes.get(index);

final double oldDeltaCostFunction;
oldDeltaCostFunction = getPositionDeltaCost(firstNode) + getPositionDeltaCost(secondNode);

swapNodes();

// return the cost function
final double newDeltaCostFunction;
newDeltaCostFunction = getPositionDeltaCost(firstNode) + getPositionDeltaCost(secondNode);

return currentCost + (newDeltaCostFunction - oldDeltaCostFunction) / costFactor;
}

private void evaluateNextStepPosition() {
// current value of the cost function
double currentCost = getGlobalCost(costFactor);
double bestCost = currentCost;
status.reset();
status.setTemperature(getCurrentTemperature());
status.incrementNumber();

for (long index = 1; index <= iterations; index++) {

final double currentCostOLD = currentCost;
currentCost = evaluateNodePositionCandidate(currentCost);
if (currentCost > bestCost) {
final double temp = (bestCost - currentCost) / status.getTemperature();
if (temp > -15.0 && Math.exp(temp) > Math.random()) {
bestCost = currentCost;
status.incrementWorse();
} else {
swapNodes();
status.incrementRejected();
currentCost = currentCostOLD;
}
} else if (currentCost < bestCost) {
bestCost = currentCost;
status.incrementBetter();
} else {
status.incrementConst();
}
}
}

private void swapNodes() {
final long xTemp = firstNode.getColumn();
final long yTemp = firstNode.getRow();
firstNode.setColumn(secondNode.getColumn());
firstNode.setRow(secondNode.getRow());
secondNode.setColumn(xTemp);
secondNode.setRow(yTemp);
}

private void evaluateNextStepCluster() {
// current value of the cost function
double currentCost = getGlobalCost(costFactor);
double bestCost = currentCost;
status.reset();
status.setTemperature(getCurrentTemperature());
status.incrementNumber();

for (long index = 1; index <= iterations; index++) {
final double currentCostOLD = currentCost;
currentCost = evaluateNodeClusterCandidate(currentCost);
if (currentCost > bestCost) {
final double temp = (bestCost - currentCost) / status.getTemperature();
if (temp > -15.0 && Math.exp(temp) > Math.random()) {
bestCost = currentCost;
status.incrementWorse();
} else {
firstNode.setCluster(firstCluster);
status.incrementRejected();
currentCost = currentCostOLD;
}
} else if (currentCost < bestCost) {
bestCost = currentCost;
status.incrementBetter();
} else {
status.incrementConst();
}
}
}

public double getTemperature() {
return temperature;
}

private double getCurrentTemperature() {
return temperature * Math.pow(1 - (double) status.getStep() / (double) steps, decay);
}

private Cluster evaluateRandomCluster(final Cluster cluster) {
final int index = 1 + (int) Math.round(Math.floor(Math.random() * (clusters.size() - 1)));
final Cluster newCluster = clusters.get(index);
if (cluster != newCluster) {
return newCluster;
} else {
return evaluateRandomCluster(cluster);
}
}

public long getLinkCount() {
return links.size();
}

public long getClusterCount() {
return clusters.size();
}

public List<Node> getNodes() {
return nodes;
}

public List<Link> getLinks() {
return links;
}

public List<Cluster> getClusters() {
return clusters;
}

public int getColumns() {
return columns;
}

public void setColumns(final int c) {
columns = c;
}

public long getNodeCount() {
return nodes.size();
}

private void outputStatus() {
final String s1 = String.format(" %3d/%3d ", status.getStep(), steps);
final String s2 = String.format("%.2e", status.getTemperature());
final String s3 = String.format(STRING_FORMAT_LONG, status.getNoChange());
final String s4 = String.format(STRING_FORMAT_LONG, status.getBetter());
final String s5 = String.format(STRING_FORMAT_LONG, status.getWorse());
final String s6 = String.format(STRING_FORMAT_LONG, status.getRejected());
final double globalCost = getGlobalCost(costFactor);
final String s7 = String.format(" %9.2f%s", globalCost * 100, "%");

if (status.getStep() != 0) {
LOGGER.info(s1 + s2 + s3 + s4 + s5 + s6 + s7);
} else {
LOGGER.info(s1 + " - - - - -" + s7);
}
}

private static void outputStatusHeader() {
LOGGER.info("PROGRESS TEMPERATUR CONST BETTER WORSE REJECTED COST");
LOGGER.info("");
}

public void run() {
costFactor = getGlobalCost(1.0);

if (LOGGER.isInfoEnabled()) {
LOGGER.info("");
LOGGER.info("-----------------------------------------------------------------------");
LOGGER.info("optimize position = " + optimizePosition);
LOGGER.info("decay = " + String.format(STRING_FORMAT_FLOAT, decay));
LOGGER.info("temperature = " + String.format(STRING_FORMAT_FLOAT, temperature));
LOGGER.info("iterations = " + iterations);
LOGGER.info("steps = " + steps);
LOGGER.info("inital cost = " + costFactor);
LOGGER.info("");
}

// start optimization
status.setStep(0);
outputStatusHeader();
outputStatus();
while (status.getStep() < steps) {
if (optimizePosition) {
evaluateNextStepPosition();
} else {
evaluateNextStepCluster();
}
outputStatus();
if (status.getBetter() == 0) {
break;
}
}

LOGGER.info("-----------------------------------------------------------------------");
LOGGER.info("");

}

public void initParameters(final boolean position) {
this.optimizePosition = position;
iterations = ITTERATIONS_PER_NODE * nodes.size();
decay = RUN_DECAY + (nodes.size() > 0 ? Math.log(nodes.size()) : 0);
temperature = TEMPERATURE;
steps = STEPS;
oldCluster = clusters.get(0);
firstCluster = oldCluster;
}

public void initNodePostion() {
setColumns((int) Math.sqrt(getNodeCount()) - 1);
int x = 0;
int y = 0;
final List<Node> nodes = getNodes();
for (final Node node : nodes) {
node.setColumn(x);
node.setRow(y);
if (x == getColumns()) {
y++;
x = 0;
} else {
x++;
}
}
}

public long getIterations() {
return iterations;
}

public void setIterations(final long i) {
this.iterations = i;
}

public long getSteps() {
return steps;
}

public void setSteps(final long s) {
this.steps = s;
}

public double getDecay() {
return decay;
}

public void setDecay(final double d) {
this.decay = d;
}

public void setTemperature(final double t) {
this.temperature = t;
}

public Status getStatus() {
return status;
}

public void reconnectClusterNodes() {
LOGGER.info(String.format("reconnect cluster-nodes"));

// find all cluster nodes for later use
final Map<String, Node> nodesMap = new HashMap<String, Node>(Model.DEFAULT_SIZE_NODE_NUMBER);
for (final Node node : getNodes()) {
if (node.isClusterNode()) {
nodesMap.put(node.getName(), node);
}
}

// Create new cluster links
for (final Node node : getNodes()) {
if ("-".equals(node.getInitialClusterName())) {
final Cluster cluster = node.getCluster();
if (!cluster.getName().equals(DEFAULT_CLUSTER_NAME)) {
final Node clusterNode = nodesMap.get("C@" + cluster.getName());
if (clusterNode != null && !node.isClusterNode()) {
links.add(0, new Link(clusterNode, node, true));
}
}
}
}
}
}

Recommendations
  • The most time consuming step in simulated annealing is the calculation of the global cost function, so always try to find a way to just calculate the delta for each new step.

Change History

 Revision  Date  Author  Description
 1.0  Jul 13, 2013  Markus Sprunck  first version

Google+ Comments

You may press the +1 button to share and/or comment