Dr. Dobb's is part of the Informa Tech Division of Informa PLC

This site is operated by a business or businesses owned by Informa PLC and all copyright resides with them. Informa PLC's registered office is 5 Howick Place, London SW1P 1WG. Registered in England and Wales. Number 8860726.


Channels ▼
RSS

Design

Simulated Annealing: A Heuristic Optimization Algorithm


Sep01: Algorithm Alley

Girish is a scientist at the Tata Research Development and Design Centre (TRDDC) in Pune, India. His research interests are in AI and mathematical logic. He can be contacted at [email protected].


There are many practical problems that need to find a lowest cost solution that minimizes some given cost function. For example, in the famous traveling salesman problem (TSP), given a map connecting N cities through a road network, the task is to find a shortest length tour starting and ending at a given city and where the salesman visits every other city exactly once. The obvious cost here is the total distance traveled in a tour. Of course, there may be additional factors that contribute to evaluating the desirability of a proposed tour; for example, the shortest time to be spent in cold regions. As another example, the placement of components on a circuit board is expected to minimize the total length of the wire used and also the number of crossings.

These so-called "combinatorial optimization" problems have some common characteristics. There is a state-space in which the search for an optimal solution needs to be conducted. Each state (or configuration) in such a state-space is a candidate solution. For example, the state- space for TSP consists of all possible tours in the given graph and the state-space for the circuit placement problem consists of all possible placements of given components in a given MXN grid, subject to connectivity in the circuit diagram. The state-space is often finite, though very large, and is characterized in terms of a small number of parameters where each parameter has a finite domain of possible values. There are domain specific, often symbolic, constraints that exclude certain states as possible solutions. The cost can also include many complex factors and dependencies. To give some idea of the state-space size, consider the problem of placing three components on a 100X50 grid; each component can be placed on a square in this grid and no two components can occupy the same square. One possible state (that is, a placement) is to place the components in squares given by the coordinates (10,10), (15,15), and (20,20). Clearly, the total number of possible placement for three components=5000X4999X998=124925010000≈125 billion; this is the size of the state-space.

It can be proved that many of these combinatorial optimization problems are intractable; that is, there are no known efficient systematic search algorithms that find the optimal solution for all input instances of the problem. Fortunately, if the optimal solution is too hard to find, a solution that is sufficiently close to it is often acceptable in practice. A common approach is then to define heuristic search algorithms that efficiently find such near optimal solutions most of the time. The heuristics refer to domain-specific knowledge or expertise to control and expedite the search. Among the many heuristic search algorithms, I will outline one rather simple approach (and present a C implementation skeleton for it) called "hill climbing," and use it to introduce some useful terms that describe common limitations of search algorithms. I then present a well-known probabilistic local search technique called "simulated annealing" (SA) that mimics the natural process of the slow cooling of liquids, which leads to a solid that has the lowest energy. Finally, I provide a general C implementation of the SA algorithm and illustrate it by applying it to the problem of producing an aesthetic drawing of a given graph.

Hill Climbing

Hill climbing is a greedy nearest neighbor approach to search the state-space to locate a lowest cost solution. It uses domain knowledge in the form of a distance function that returns an estimate of how close a given state is to the global minimum. This function is provided separately by the user for each problem domain; for example, the distance function for the placement problem is different from that for TSP. Listing One is a C program skeleton of the hill-climbing search algorithm, which is actually for a specific version called the "steepest ascent" hill climbing or "gradient search." The problem-specific definitions of the C structure STATE and the assumed functions need to be provided for using this code.

It is easy to see the limitations of the hill-climbing algorithm for locating the global maximum. Hill climbing often reaches a local minimum — that is, a state that has a lower cost than all its neighbors, but which is not a global minimum. It may also reach a plateau — a state with neighboring states that all have the same cost and hence no next state can be chosen for further exploration. The hill-climbing algorithm can be modified to include simple strategies to deal with these problems; backtracking to some earlier state, taking a jump to some random state (when a plateau is encountered), moving in several directions simultaneously, and so on. Nevertheless, in general, hill climbing has been found to be limited in practical combinatorial optimization problems. Simulated annealing (SA) is another local search algorithm (among a variety of other techniques) that overcomes some of the limitations of hill climbing and often succeeds in efficiently locating good low-cost solutions. (For more information, see "Simulated Annealing," by Michael P. McClauphlin, DDJ, September 1989.)

Simulated Annealing

It is well known in physics that when a hot liquid is cooled slowly (annealed) it reaches a crystalline form. The atoms (or molecules) in a crystal are arranged in a regular and ordered configuration; this state has the lowest (potential+kinetic) energy and consequently, highest stability. On the other hand, when the hot liquid is cooled rapidly, the end result is an amorphous solid — this state has a somewhat higher energy. Qualitatively speaking, this happens because when the liquid is cooled slowly, the atoms have time to rearrange themselves until they achieve a state of lowest energy at that temperature — the so-called "thermal equilibrium." In a thermal equilibrium at temperature T (in K), the system achieves a state of energy E with a probability p(E)e-E/kT where k=1.3804X10-16 erg per degree is a constant of nature, called the "Boltzmann" constant. An interesting thing to note about this Boltzmann distribution is that it lets the system be in a higher energy state, even in a thermal equilibrium, although with a small chance. This chance decreases with temperature T, so that the lower the temperature, the higher the chance that the system is in the lowest energy state.

"Optimization by Simulated Annealing," by S. Kirkpatrick et al. (Science, 1983), identified the similarities between the annealing process and a general search for optimal solution in a combinatorial optimization problem. Essentially, the probability, at a temperature T, with which the search algorithm moves from a state s1 of energy E1 to another state s2 of energy E2 is given by e-(E1-E2)/kT. Clearly, whenever E2<E1 the system will definitely move to s2 but when E2>E1, the state change is probabilistic.

For applying search techniques to combinatorial optimization problems, you need to:

  • Clearly define the information associated with each state.
  • Identify an initial state, which is often generated randomly.

  • Define a function that (often probabilistically) returns a state in the neighborhood of a given state.

  • Define a Boolean function that returns True if the given state is a solution (that is, a minimum cost state) and False, otherwise.

In addition, for applying SA to combinatorial optimization problems, you need to define analogues of temperature and cost. A common approach is to define energy as a function that returns the cost for a given state; the definition of the energy function includes all relevant factors that affect the cost. The Boltzmann constant k is often taken to be 1. The notion of temperature is often purely artificial (does not correspond to any real aspect of the problem) and is often defined to vary from a given maximum value Tmax to a given minimum value Tmin. A number of cooling or annealing strategies are used to systematically drop the temperature from Tmax until it reaches Tmin; the most common is the "geometric annealing schedule," where the next temperature is obtained from the current temperature by the geometric series formula Tnext=gXTcurrent, where g is a given constant between 0 and 1. For instance, when Tmax=1000, Tmin=10 and g=0.5, the sequence of temperatures obtained using geometric annealing schedules will be 1000, 500, 250, 125, 62.5, 31.25, 15.625, and 7.8125, which is less than Tmin; at this point the schedule terminates. Figure 1 is a high-level description of the SA algorithm. The termination condition is often T<Tmin; but additional factors can be included here. Also, since the expression e-(E-E')/kT need not be always between 0 and 1, it has to be scaled appropriately so that comparison with the output of a standard random number can be made.

Listings Two and Three are a C implementation of the SA algorithm. To apply this C implementation to a specific combinatorial problem, you need to provide some additional problem-specific details. I will illustrate this process by applying this C implementation to a specific problem.

Aesthetic Graph Drawing Using SA

A number of practical applications involving combinatorial optimization have used SA. Here I describe an interesting problem of automatically generating an aesthetic drawing of a given graph on an MXN grid (say, on screen), as discussed in "Drawing Graphs Nicely Using Simulated Annealing," by R. Davidson and D. Harel (ACM Transaction Graphics, October 1996). The task here is to take a given graph of n vertices and e edges (representing connectivity between vertices) and place the n vertices on an MXN grid such that when the adjacent vertices are connected by straight line edges, the resulting graph looks aesthetically neat. Of course, there are many criteria for what constitutes an aesthetic drawing of a graph: Vertices are evenly distributed in the grid, not too many are near the center or near the borders, edge lengths are not too small nor too long, edge crossings are minimized, and so on.

Additional C source code needed for applying the SA algorithm to the graph- drawing problem is available electronically (see "Resource Center," page 5). The C structure GRAPH contains the connectivity of the given graph. A state (as stored in the C structure STATE) represents a placement of all the n vertices of a given graph on the MXN grid; that is, a state associates a position, in the form of an x- and a y-coordinate, with each vertex in a given graph. The RandomState function generates a random placement for all the vertices and is used to generate the initial state. The NeighbourState function generates a neighbor state (for the given state) by randomly selecting a vertex and moving it by one square (up, down, right, left, up-right, up-left, down-right, or down-left).

The Energy function computes the cost associated with a given placement of a given graph. Several factors contribute to the calculation of the energy E(s) for a state s; I've implemented only the three, as in the formula for E(s) (Figure 2).

In Figure 2, n and e are the number of vertices and edges in the given graph, dij is the Euclidean distance between vertex i and vertex j (in the given placement), ri, li, ti, and bi are the distances of the ith vertex from the right, left, top, and bottom borders, respectively, dk is the length of the kth edge and 1, 2, and 3 are the parameters used to control the relative importance of the three factors. Contribution of the first factor increases as the vertices move closer to each other and thus SA will try to prevent vertices from coming too close to each other. Contribution of the second factor to E(s) increases as a vertex moves closer to the borders of the MXN grid and thus SA will try to prevent vertices from coming too close to the borders. Contribution of the third factor to E(s) increases as lengths of edges increase and thus SA will try to prevent edges that are too long. Davidson and Harel suggest more factors that contribute to E(s) and those can be easily added. Increasing 1 relative to other parameters causes SA to prefer pictures with smaller distances between vertices. Increasing 2 relative to 1 causes SA to push the vertices towards the center and decreasing it causes SA to use more drawing space near the borders. Increasing 3 causes SA to prefer shorter edges. When either dij=0 or when either of ri, li, ti, or bi is 0, then you assume that the corresponding energy takes a large value denoted in the program by the constant MAX_ENERGY.

The InitGraph function fills up a global variable with the input graph that is to be placed. The SimAnneal function implements the SA algorithm. Tmax and Tmin are encoded in the constants MAX_TEMP and MIN_TEMP, respectively.

As an example of the application of the SA algorithm, Figure 3 is a simple graph and its placements as generated by the SA algorithm on a 30X20 grid. Note that the algorithm did not generate a hexagonal placement (possibly because it did not attempt to prevent edge crossings); however, the generated placement does look symmetrical in a different way. By adding more factors to the energy calculations and by selecting appropriate values for the parameters (for 1, 2, and 3), the SA algorithm can be used to generate pleasing placements for more complex graphs; see Davidson and Harel. Even further complex factors that characterize aesthetic drawing of a graph can also be added; for example, for balancing the area occupied by the convex hull of the graph and the remaining blank area.

Conclusion

A number of other heuristic algorithms have been invented for combinatorial optimization problems — AO* search, constraint satisfaction, means-end analysis, genetic algorithms, and so on. These algorithms yield acceptable near-optimal solutions for a variety of problems. But there is no known universal algorithm that works efficiently for all problems. We examined a particular probabilistic heuristic optimization algorithm, called "simulated annealing." SA mimics the natural process of slow cooling of liquids that leads to a solid form that has the lowest energy.

A number of possible enhancements have been attempted for SA. The SA algorithm has been parallelized in various ways to obtain speedup benefits on multiprocessor systems. Experimenting with different types of annealing schedules is another possible direction.

SA itself is not always suitable for all combinatorial optimization problems. Empirically it has been observed that for SA to work satisfactorily, the cost function should not contain narrow and steep valleys. The acceptable near minimal solutions should not be within such narrow and steep valleys. The energy (cost) function should change smoothly upon changes in states. It should also be efficiently computable, since it is invoked a large number of times. SA usually finds a near optimal solution, but not the global minimum itself if it has a low probability of being found. Finding suitable temperature schedules and best values of control factors for SA often needs careful experimentation. In spite of these problems, SA is being successfully used in a large number of practical problems, including VLSI chip design and layout, floor planning, flow shop scheduling, channel routing, graph drawing, image processing, coding theory, graph coloring and partitioning, satisfiability, and so on.

Acknowledgments

I gratefully acknowledge the support of Professor Mathai Joseph. I also wish to thank Dr. Manasee Palshikar for her patience and strength.

DDJ

Listing One

// max. no. of neighboring states for a given state
#define MAX_NEIGHBOUR   100     

// external functions assumed (specific to each problem domain)
extern double Cost(STATE *pS);  // returns the cost of given state
extern int Solution(STATE *pS);  // 1 if given state is a solution; else 0 
// generate neighbor states of given state; return no. of neighbor states
extern int GetNeighbours(STATE *pS, STATE Next[MAX_NEIGHBOUR]); 

// hopefully find and return the lowest cost state
// STATE is a problem-specific representation of the structure of a state 
// pS0 and pS are pointers to given initial state and low cost state found. 
int HillClimbing(STATE *pS0, STATE *pS)
{
   STATE Next[MAX_NEIGHBOUR];
   int i, n, index;
   double c0, c, c1;

   if ( Solution(pS0) ) // found
   {
       CopyState(pS,pS0);
       return(1);
   }
   CopyState(pS, pS0); // initialize
   c0 = Cost(pS0);     // cost of initial state
   while ( (n = GetNeighbours(pS,Next)) > 0 ) // get neighbours of s
   {
       index = 0; // index (in Next) of lowest cost neighbour
       c = Cost(&Next[0]);
       for(i = 0; i < n; i++) // do for each neighbour state
       {
       if ( Solution(&Next[i]) ) // found
       {
           CopyState(pS,&Next[i]);
           return(1);
       }
       if ( (c1 = Cost(&Next[i])) < c )
       {
           index = i;
           c = c1;
       }
       } // end for 
       if ( c < c0 ) // found a lower cost neighbour
       {
       CopyState(pS, &Next[index]);
       c0 = c;
       }
       else // reached local maximum
       break;
   } // end while
   return(0); // global solution not found; S contains low cost state found 
}

Back to Article

Listing Two

/* Module: Simulated Annealing
 * Programmer: Girish Keshav Palshikar
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "simann.h"
#include "sim_app.h"    // problem-specific information needed by SA

// function prototypes 
int NextTemp(TEMP *pt);       
BOOL Boltzmann(ENERGY e1, ENERGY e2, TEMP t);
BOOL CoolEnough(TEMP t); 

// problem-specific functions assumed
extern BOOL NeighbourState(STATE *pS, STATE *pNew);
extern BOOL IsSolution(STATE *pS); // TRUE if state is an acceptable solution
extern ENERGY Energy(STATE *pS); // compute & return cost/energy of state
extern int CopyState(STATE *pDest, STATE *pSrc);

// s0 is the initial state and t0 is the initial temperature
// Limit is the no. of attempts to explore neighbourhoods at a given temp. 
BOOL SimAnneal(STATE *pInit, TEMP t0, int Limit, STATE *pDest)
{
   STATE s;  // current state
   TEMP t;   // current temp.
   ENERGY e, e_new; // energy of current and new state
   BOOL success = FALSE;
   int i;

   CopyState(&s, pInit);  // make current state = initial state
   e = Energy(&s); // compute current energy
   t = t0;
   CopyState(pDest, pInit);  // copy initial state into pDest
   
   while ( !CoolEnough(t) ) 
   {
       if ( IsSolution(&s) ) // done
       {
       success = TRUE;
       break;
       }
       for ( i = 0, success = FALSE; i < Limit; i++ ) 
       {
       if ( !NeighbourState(&s, pDest) ) // generate new state in pDest
           break;   // no more states in the neighbourhood 
       if ( IsSolution(pDest) ) // done
       {
           success = TRUE;
           break;
       }
       e_new = Energy(pDest);
       printf("temp = %lf i = %d e = %lf e_new = %lf\n", t, i, e, e_new); 
       if ( Boltzmann(e, e_new, t) ) // make pDest new current state?
       {    
           CopyState(&s, pDest);
           e = e_new;
       }
       }  // end for
       NextTemp(&t);  // get the next lower temperature
   } // end while          
   CopyState(pDest, &s); // copy last state into pDest
   return(success);
}
BOOL Boltzmann(ENERGY e1, ENERGY e2, TEMP t)
{
   double x;
   if ( e2 < e1 ) return(TRUE);
   x = exp( (-(e2 - e1)) / t );
   if ( x <= 1.0 && rand() < x ) return(TRUE);
   if ( rand() < (x - floor(x)) ) // x > 1, so check decimal fraction of x
       return(TRUE); 
   return(FALSE);
}
// implement cooling schedule; return the next lower temp
int NextTemp(TEMP *pt)       
{
    TEMP t1 = (*pt);
    (*pt) = t1 * TEMP_FACTOR;
    return(0);
}
// true if given temp < a fixed threshold; stop condition for SA 
BOOL CoolEnough(TEMP t) 
{ 
   if ( t < MIN_TEMP ) return(TRUE);
   return(FALSE);
}

Back to Article

Listing Three

/* Module: Simulated Annealing
 * File: simann.h
 * Programmer: Girish Keshav Palshikar
*/
#ifndef SIMANN_H
#define SIMANN_H

typedef double TEMP;
typedef double ENERGY;
typedef enum { FALSE = 0, TRUE = 1 } BOOL;

#define MAX_TEMP           5000.0
#define MIN_TEMP           1.0
#define BOLTZMANN_CONSTANT 1.0
#define TEMP_FACTOR        0.5

#endif

Back to Article


Related Reading


More Insights






Currently we allow the following HTML tags in comments:

Single tags

These tags can be used alone and don't need an ending tag.

<br> Defines a single line break

<hr> Defines a horizontal line

Matching tags

These require an ending tag - e.g. <i>italic text</i>

<a> Defines an anchor

<b> Defines bold text

<big> Defines big text

<blockquote> Defines a long quotation

<caption> Defines a table caption

<cite> Defines a citation

<code> Defines computer code text

<em> Defines emphasized text

<fieldset> Defines a border around elements in a form

<h1> This is heading 1

<h2> This is heading 2

<h3> This is heading 3

<h4> This is heading 4

<h5> This is heading 5

<h6> This is heading 6

<i> Defines italic text

<p> Defines a paragraph

<pre> Defines preformatted text

<q> Defines a short quotation

<samp> Defines sample computer code text

<small> Defines small text

<span> Defines a section in a document

<s> Defines strikethrough text

<strike> Defines strikethrough text

<strong> Defines strong text

<sub> Defines subscripted text

<sup> Defines superscripted text

<u> Defines underlined text

Dr. Dobb's encourages readers to engage in spirited, healthy debate, including taking us to task. However, Dr. Dobb's moderates all comments posted to our site, and reserves the right to modify or remove any content that it determines to be derogatory, offensive, inflammatory, vulgar, irrelevant/off-topic, racist or obvious marketing or spam. Dr. Dobb's further reserves the right to disable the profile of any commenter participating in said activities.

 
Disqus Tips To upload an avatar photo, first complete your Disqus profile. | View the list of supported HTML tags you can use to style comments. | Please read our commenting policy.