Simulated Annealing
The simulated annealing (SA) algorithm was inspired by the physical annealing process of solids.
Annealing is a thermal process used to reduce the hardness and improve the workability of a materials.
Molding and forming a metal object leads to the creation of crystalline structures (a.k.a. “grains”) to form within the material.
When cooled quickly, these structures are smaller, leading to hard but brittle materials — more prone to crack under stress along the juncture lines between grains.
The annealing process increases the size of these crystal grains, making the material more ameaneable to stress. It consists of two steps:
- Heat the material above its recrystallization temperature but below melting temperature.
- Cool the material carefully such that new crystall grains can develop.
In the SA algorithm:
- The candidate solutions of a combinatorial optimization problem are the states of the physical system (i.e., the arrange of the crystall grains).
- The cost of a solution is equivalent to the energy of a state.
Setting up the resources used throughout the document:
Outline
- Initialization. Start with a feasible initial trial solution.
- Iteration. Use the move selection rule to select the next trial solution. (If none of the immediate neighbors of the current trial solution are accepted, the algorithm is terminated.)
- Check the temperature schedule. When the desired number of iterations have been performed at the current temperature \(T\), decrease \(T\) to the next value in the temperature schedule and resume performing iterations at this next value.
Acceptance criteria
Outline:
- \(Z_c\) = objective function value for the current trial solution,
- \(Z_n\) = objective function value for the current candidate to be the next trial solution,
- \(T\) (temperature) = a control parameter that measures the tendency to accept the current candidate to be the next trial solution if this candidate is not an improvement on the current trial solution.
Move selection rule (assuming MINIMIZATION): - If \(Z_n \leq Z_c\), always accept the candidate. - If \(Z_n > Z_c\), accept the candidate with the following probability \(e^x\) where \(x = \frac{Z_c - Z_n}{T}\) (Boltzmann distribution)
Hence:
\[ \begin{align} \text{Prob. acceptance} = \left\{ \begin{array}{cl} 1 & Z_n \leq Z_c, \\ e^{\frac{Z_c - Z_n}{T}} & Z_n > Z_c. \end{array} \right. \end{align} \]
In the following, an implementation example of the acceptance criteria:
Probability of acceptance vs. solution quality
Let us plot the probability of acceptance for a range of new candidate solutions \(Z_n = Z_c, \ Z_c+1,\ Z_c+2, \dots\) (i.e., worse than a current solution \(Z_c = 20\)) assuming temperature \(T=100\):
Probability of acceptance vs. temperature
The probability of acceptance increases the higher is the temperature \(T\). For example, for \(Z_c = 20\) and \(Z_n = 40\) (i.e., a bad quality solution) we have:
T=200: Probability of acceptance of Z_n = 40: 0.995
T=100: Probability of acceptance of Z_n = 40: 0.990
T= 20: Probability of acceptance of Z_n = 40: 0.951
T= 10: Probability of acceptance of Z_n = 40: 0.905
T= 5: Probability of acceptance of Z_n = 40: 0.819Below, we show the probability of acceptance for a range of candidate solutions for different values of \(T\). First, we generate random candidate solutions \(n\) with costs \(Z_n \geq 0\) randomly distributed around the value of the current solution \(Z_c\):
Notice that the probability of acceptance decreases sharpely for low temperatures:
Study cases
In the following, we present two study cases:
- Random city distribution: you can set the number of nodes that will be spread randomly in a \([0,1]\times[0,1]\) plane.
- 48 US capitals: a classic TSP instance comprised of the 48 US capitals with minimal tour length 33,523. You can find more classic TSP instances here.
Random city distribution
Set parameter Username
Set parameter LicenseID to value 2617855
Academic license - for non-commercial use only - expires 2026-02-04
Set parameter LazyConstraints to value 1
Gurobi Optimizer version 12.0.1 build v12.0.1rc0 (win64 - Windows 11.0 (26100.2))
CPU model: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Non-default parameters:
LazyConstraints 1
Optimize a model with 61 rows, 1830 columns and 3660 nonzeros
Model fingerprint: 0xbbe19efa
Variable types: 0 continuous, 1830 integer (1830 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [2e-03, 1e+00]
Bounds range [1e+00, 1e+00]
RHS range [2e+00, 2e+00]
Presolve time: 0.00s
Presolved: 61 rows, 1830 columns, 3660 nonzeros
Variable types: 0 continuous, 1830 integer (1830 binary)
Root relaxation: objective 5.100682e+00, 83 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 5.10068 0 6 - 5.10068 - - 0s
0 0 5.19704 0 16 - 5.19704 - - 0s
0 0 5.22478 0 16 - 5.22478 - - 0s
0 0 5.26602 0 8 - 5.26602 - - 0s
0 0 5.26885 0 8 - 5.26885 - - 0s
0 2 5.26885 0 8 - 5.26885 - - 0s
* 6 2 2 5.4200636 5.37706 0.79% 6.3 0s
Cutting planes:
Zero half: 9
Lazy constraints: 13
Explored 10 nodes (167 simplex iterations) in 0.17 seconds (0.03 work units)
Thread count was 8 (of 8 available processors)
Solution count 1: 5.42006
Optimal solution found (tolerance 1.00e-04)
Best objective 5.420063643805e+00, best bound 5.420063643805e+00, gap 0.0000%
User-callback calls 222, time in user-callback 0.04 sec
48 US capitals
[ 0 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
[(6734, 1453), (2233, 10), (5530, 1424), (401, 841), (3082, 1644), (7608, 4458), (7573, 3716), (7265, 1268), (6898, 1885), (1112, 2049), (5468, 2606), (5989, 2873), (4706, 2674), (4612, 2035), (6347, 2683), (6107, 669), (7611, 5184), (7462, 3590), (7732, 4723), (5900, 3561), (4483, 3369), (6101, 1110), (5199, 2182), (1633, 2809), (4307, 2322), (675, 1006), (7555, 4819), (7541, 3981), (3177, 756), (7352, 4506), (7545, 2801), (3245, 3305), (6426, 3173), (4608, 1198), (23, 2216), (7248, 3779), (7762, 4595), (7392, 2244), (3484, 2829), (6271, 2135), (4985, 140), (1916, 1569), (7280, 4899), (7509, 3239), (10, 2676), (6807, 2993), (5185, 3258), (3023, 1942)] 48
Set parameter LazyConstraints to value 1
Gurobi Optimizer version 12.0.1 build v12.0.1rc0 (win64 - Windows 11.0 (26100.2))
CPU model: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Non-default parameters:
LazyConstraints 1
Optimize a model with 48 rows, 1128 columns and 2256 nonzeros
Model fingerprint: 0xdc564e0d
Variable types: 0 continuous, 1128 integer (1128 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [1e+02, 8e+03]
Bounds range [1e+00, 1e+00]
RHS range [2e+00, 2e+00]
Presolve time: 0.00s
Presolved: 48 rows, 1128 columns, 2256 nonzeros
Variable types: 0 continuous, 1128 integer (1128 binary)
Root relaxation: objective 3.166926e+04, 78 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 31669.2574 0 14 - 31669.2574 - - 0s
H 0 0 38335.094169 31669.2574 17.4% - 0s
H 0 0 37725.521880 31669.2574 16.1% - 0s
H 0 0 37029.139102 31669.2574 14.5% - 0s
H 0 0 36754.544644 31669.2574 13.8% - 0s
0 0 32320.8108 0 6 36754.5446 32320.8108 12.1% - 0s
H 0 0 35949.039512 32320.8108 10.1% - 0s
H 0 0 35606.532587 32320.8108 9.23% - 0s
H 0 0 35541.897685 32320.8108 9.06% - 0s
0 0 33021.6741 0 - 35541.8977 33021.6741 7.09% - 0s
0 0 33035.1856 0 - 35541.8977 33035.1856 7.05% - 0s
0 0 33043.1273 0 - 35541.8977 33043.1273 7.03% - 0s
0 0 33081.8621 0 - 35541.8977 33081.8621 6.92% - 0s
0 0 33255.2634 0 20 35541.8977 33255.2634 6.43% - 0s
0 0 33255.2634 0 14 35541.8977 33255.2634 6.43% - 0s
H 0 0 34966.443166 33255.2634 4.89% - 0s
0 0 33299.9040 0 12 34966.4432 33299.9040 4.77% - 0s
0 0 33299.9040 0 12 34966.4432 33299.9040 4.77% - 0s
H 0 0 34047.150934 33299.9040 2.19% - 0s
0 0 33434.9079 0 22 34047.1509 33434.9079 1.80% - 0s
0 0 33434.9079 0 12 34047.1509 33434.9079 1.80% - 0s
0 0 33484.2144 0 - 34047.1509 33484.2144 1.65% - 0s
0 0 33508.6532 0 33 34047.1509 33508.6532 1.58% - 0s
H 0 0 33706.845991 33508.6532 0.59% - 0s
0 0 33511.1608 0 19 33706.8460 33511.1608 0.58% - 0s
* 0 0 0 33523.708507 33523.7085 0.00% - 0s
Cutting planes:
Gomory: 6
MIR: 1
Zero half: 9
Lazy constraints: 2
Explored 1 nodes (442 simplex iterations) in 0.36 seconds (0.05 work units)
Thread count was 8 (of 8 available processors)
Solution count 10: 33523.7 33706.8 34047.2 ... 37725.5
Optimal solution found (tolerance 1.00e-04)
Best objective 3.352370850744e+04, best bound 3.352370850744e+04, gap 0.0000%
User-callback calls 322, time in user-callback 0.12 sec
Tuning
Initial solution
We are going to compare the outcome of the SA algorithm starting from different initial solutions generated using:
- a random approach,
- the farthest addition algorithm,
- the neighborhood search algorithm (apply 2-Opt on a solution generated by the FA algorithm).
SA is executed starting from each of these initial routes:
The SA was able to substantially improve solutions 1 (random route), 2 (farthest addition), and 3 (neighborhood search) in the 48 US capitals TSP instance.
Initial temperature
Performance comparison starting from different temperatures:
We use the best initial temperature found to modify other SA parameters:
T=200 - Cost = 33,632.13Previously, the SA would spend 1,000 iterations at each temperature. Now, we test a range of different iterations:
Elapsed time for n. of iterations = 1000:- Runtime = 5.89s- Cost = 33,632.13
Elapsed time for n. of iterations = 3000:- Runtime = 16.59s- Cost = 33,632.13
Elapsed time for n. of iterations = 5000:- Runtime = 60.64s- Cost = 33,632.13
Elapsed time for n. of iterations = 10000:- Runtime = 61.76s- Cost = 33,632.13