Simulated annealing
Quoted from the Wikipedia page : Simulated annealing (SA) is a probabilistic technique for approximating the global optimum of a given function. Simulated annealing interprets slow cooling as a slow decrease in the probability of temporarily accepting worse solutions as it explores the solution space.
Pseudo code from Wikipedia
Notations :
T : temperature. Decreases to 0.
s : a system state
E(s) : Energy at s. The function we want to minimize
∆E : variation of E, from state s to state s_next
P(∆E , T) : Probability to move from s to s_next.
if ( ∆E < 0 ) P = 1
else P = exp ( - ∆E / T) . Decreases as T → 0
Pseudo-code:
Let s = s0 -- initial state
For k = 0 through kmax (exclusive):
T ← temperature(k , kmax)
Pick a random neighbour state , s_next ← neighbour(s)
∆E ← E(s) - E(s_next)
If P(∆E , T) ≥ random(0, 1), move to the new state:
s ← s_next
Output: the final state s
Problem statement
We want to apply SA to the travelling salesman problem. There are 100 cities, numbered 0 to 99, located on a plane, at integer coordinates i,j : 0 <= i,j < 10 . The city at (i,j) has number 10*i + j. The cities are all connected : the graph is complete : you can go from one city to any other city in one step.
The salesman wants to start from city 0, visit all cities, each one time, and go back to city 0. The travel cost between two cities is the euclidian distance between there cities. The total travel cost is the total path length.
A path s is a sequence (0 a b ...z 0) where (a b ..z) is a permutation of the numbers (1 2 .. 99). The path length = E(s) is the sum d(0,a) + d(a,b) + ... + d(z,0) , where d(u,v) is the distance between two cities. Naturally, we want to minimize E(s).
Definition : The neighbours of a city are the closest cities at distance 1 horizontally/vertically, or √2 diagonally. A corner city (0,9,90,99) has 3 neighbours. A center city has 8 neighbours.
Distances between cities d ( 0, 7) → 7 d ( 0, 99) → 12.7279 d ( 23, 78) → 7.0711 d ( 33, 44) → 1.4142 // sqrt(2)
Task
Apply SA to the travelling salesman problem, using the following set of parameters/functions :
- kT = 1
- temperature (k, kmax) = kT * (1 - k/kmax)
- neighbour (s) : Pick a random city u > 0 . Pick a random neighbour city v > 0 of u , among u's 8 (max) neighbours on the grid. Swap u and v in s . This gives the new state s_next.
- kmax = 1000_000
- s0 = a random permutation
For k = 0 to kmax by step kmax/10 , display k, T, E(s). Display the final state s_final, and E(s_final).
Illustrated example Temperature charts
Numerical example
kT = 1 E(s0) = 529.9158 k: 0 T: 1 Es: 529.9158 k: 100000 T: 0.9 Es: 201.1726 k: 200000 T: 0.8 Es: 178.1723 k: 300000 T: 0.7 Es: 154.7069 k: 400000 T: 0.6 Es: 148.1412 k: 500000 T: 0.5 Es: 133.856 k: 600000 T: 0.4 Es: 129.5684 k: 700000 T: 0.3 Es: 112.6919 k: 800000 T: 0.2 Es: 105.799 k: 900000 T: 0.1 Es: 102.8284 k: 1000000 T: 0 Es: 102.2426 E(s_final) = 102.2426 Path s_final = ( 0 10 11 21 31 20 30 40 50 60 70 80 90 91 81 71 73 83 84 74 64 54 55 65 75 76 66 67 77 78 68 58 48 47 57 56 46 36 37 27 26 16 15 5 6 7 17 18 8 9 19 29 28 38 39 49 59 69 79 89 99 98 88 87 97 96 86 85 95 94 93 92 82 72 62 61 51 41 42 52 63 53 43 32 22 12 13 23 33 34 44 45 35 25 24 14 4 3 2 1 0)
Extra credit
Tune the parameters kT, kmax, or use different temperature() and/or neighbour() functions to demonstrate a quicker convergence, or a better optimum.
EchoLisp
<lang scheme> (lib 'math)
- distances
(define (d ci cj) (distance (% ci 10) (quotient ci 10) (% cj 10) (quotient cj 10))) (define _dists (build-vector 10000 (lambda (ij) (d (quotient ij 100) (% ij 100))))) (define-syntax-rule (dist ci cj) [_dists (+ ci (* 100 cj))])
- E(s) = length(path)
(define (Es path) (define lpath (vector->list path)) (for/sum ((ci lpath) (cj (rest lpath))) (dist ci cj)))
- temperature() function
(define (T k kmax kT) (* kT (- 1 (// k kmax))))
- |
- alternative temperature()
- must be decreasing with k increasing and → 0
(define (T k kmax kT) (* kT (- 1 (sin (* PI/2 (// k kmax)))))) |#
- ∆E = Es_new - Es_old > 0
- probability to move if ∆E > 0, → 0 when T → 0 (frozen state)
(define (P ∆E k kmax kT) (exp (// (- ∆E ) (T k kmax kT))))
- ∆E from path ( .. a u b .. c v d ..) to (.. a v b ... c u d ..)
- ∆E before swapping (u,v)
- Quicker than Es(s_next) - Es(s)
(define (dE s u v)
- old
(define a (dist [s (1- u)] [s u])) (define b (dist [s (1+ u)] [s u])) (define c (dist [s (1- v)] [s v])) (define d (dist [s (1+ v)] [s v]))
- new
(define na (dist [s (1- u)] [s v])) (define nb (dist [s (1+ u)] [s v])) (define nc (dist [s (1- v)] [s u])) (define nd (dist [s (1+ v)] [s u]))
(cond ((= v (1+ u)) (- (+ na nd) (+ a d))) ((= u (1+ v)) (- (+ nc nb) (+ c b))) (else (- (+ na nb nc nd) (+ a b c d)))))
- all 8 neighbours
(define dirs #(1 -1 10 -10 9 11 -11 -9))
(define (sa kmax (kT 10)) (define s (list->vector (cons 0 (append (shuffle (range 1 100)) 0)))) (printf "E(s0) %d" (Es s)) ;; random starter (define Emin (Es s)) ;; E0
(for ((k kmax)) (when (zero? (% k (/ kmax 10))) (printf "k: %10d T: %8.4d Es: %8.4d" k (T k kmax kT) (Es s)) )
(define u (1+ (random 99))) ;; city index 1 99 (define cv (+ [s u] [dirs (random 8)])) ;; city number #:continue (or (> cv 99) (<= cv 0)) #:continue (> (dist [s u] cv) 5) ;; check true neighbour (eg 0 9) (define v (vector-index cv s 1)) ;; city index
(define ∆e (dE s u v)) (when (or (< ∆e 0) ;; always move if negative (>= (P ∆e k kmax kT) (random))) (vector-swap! s u v) (+= Emin ∆e))
;; (assert (= (round Emin) (round (Es s)))) ) ;; for
(printf "k: %10d T: %8.4d Es: %8.4d" kmax (T (1- kmax) kmax kT) (Es s)) (s-plot s 0) (printf "E(s_final) %d" Emin) (writeln 'Path s)) </lang>
- Output:
(sa 1000000 1) E(s0) 501.0909 k: 0 T: 1 Es: 501.0909 k: 100000 T: 0.9 Es: 167.3632 k: 200000 T: 0.8 Es: 160.7791 k: 300000 T: 0.7 Es: 166.8746 k: 400000 T: 0.6 Es: 142.579 k: 500000 T: 0.5 Es: 131.0657 k: 600000 T: 0.4 Es: 116.9214 k: 700000 T: 0.3 Es: 110.8569 k: 800000 T: 0.2 Es: 103.3137 k: 900000 T: 0.1 Es: 102.4853 k: 1000000 T: 0 Es: 102.4853 E(s_final) 102.4853 Path #( 0 10 20 30 40 50 60 70 71 61 62 53 63 64 54 44 45 55 65 74 84 83 73 72 82 81 80 90 91 92 93 94 95 85 75 76 86 96 97 98 99 88 89 79 69 59 49 48 47 57 58 68 78 87 77 67 66 56 46 36 35 25 24 34 33 32 43 42 52 51 41 31 21 11 12 22 23 13 14 15 16 17 26 27 37 38 39 29 28 18 19 9 8 7 6 5 4 3 2 1 0)
J
Implementation:
<lang J>dist=: +/&.:*:@:-"1/~10 10#:i.100
satsp=:4 :0
kT=. 1
pathcost=. [: +/ 2 {&y@<\ 0 , ] , 0:
neighbors=. 0 (0}"1) y e. 1 2{/:~~.,y
s=. (?~#y)-.0
d=. pathcost s
step=. x%10
for_k. 1+i.x do.
T=. kT*1-k%x
u=. ({~ ?@#)s
v=. ({~ ?@#)I.u{neighbors
sk=. (?0 do.
s=.sk
d=.dk
end.
if. 0=step|k do.
echo k,T,d
end.
end.
0,s,0
)</lang>
Notes:
We do not do anything special for negative deltaE because the exponential will be greater than 1 for that case and that will always be greater than our random number from the range 0..1.
Also, while we leave connection distances (and, thus, number of cities) as a parameter, some other aspects of this problem made more sense when included in the implementation:
We leave city 0 out of our data structure, since it can't appear in the middle of our path. But we bring it back in when computing path distance.
Neighbors are any city which have one of the two closest non-zero distances from the current city (and specifically excluding city 0, since that is anchored as our start and end city).
Sample run:
<lang J> 1e6 satsp dist 100000 0.9 184.641 200000 0.8 174.6 300000 0.7 146.691 400000 0.6 150.49 500000 0.5 135.62 600000 0.4 126.064 700000 0.3 113.819 800000 0.2 104.728 900000 0.1 103.314 1e6 0 103.071 0 10 11 21 20 30 31 32 33 34 35 45 55 65 66 67 77 87 88 78 68 59 69 79 89 99 98 97 96 86 76 75 85 95 94 84 83 93 92 91 90 80 70 71 81 82 72 62 52 51 61 60 50 40 41 42 43 53 63 73 74 64 54 44 46 36 37 27 16 6 7 8 9 19 29 39 49 58 57 56 47 48 38 28 18 17 26 25 15 5 4 14 24 23 22 12 13 3 2 1 0</lang>