Read More‎ > ‎

Simulated Annealing the Swiss Army Knife of Global Optimization

BY MARKUS SPRUNCK

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. 

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: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
	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:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
	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 Model.java is the main class with Simulated Annealing Algorithm

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
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
 1.1  Dec 17, 2014  Markus Sprunck   improve structure

Sponsored Link