user
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 them. 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 (Multiplication by kT is a placeholder, representing computing temperature as a function of 1-k/kmax):
* 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).

You will see that the Energy may grow to a local optimum, before decreasing to a global optimum.

**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:  158.1412 <== local optimum
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.
assistant
(formerly Perl 6)

**Translation of**: Go

```
# simulation setup
my \cities = 100;  # number of cities
my \kmax   = 1e6;  # iterations to run
my \kT     = 1;    # initial 'temperature'

die 'City count must be a perfect square.' if cities.sqrt != cities.sqrt.Int;

# locations of (up to) 8 neighbors, with grid size derived from number of cities
my \gs = cities.sqrt;
my \neighbors = [1, -1, gs, -gs, gs-1, gs+1, -(gs+1), -(gs-1)];

# matrix of distances between cities
my \D = [;];
for 0 ..^ cities² -> \j {
    my (\ab, \cd)       = (j/cities, j%cities)».Int;
    my (\a, \b, \c, \d) = (ab/gs, ab%gs, cd/gs, cd%gs)».Int;
    D[ab;cd] = sqrt (a - c)² + (b - d)²
}

sub T(\k, \kmax, \kT)      { (1 - k/kmax) × kT }                                 # temperature function, decreases to 0
sub P(\ΔE, \k, \kmax, \kT) { exp( -ΔE / T(k, kmax, kT)) }                        # probability to move from s to s_next
sub Es(\path)              { sum (D[ path[$_]; path[$_+1] ] for 0 ..^ +path-1) } # energy at s, to be minimized

# variation of E, from state s to state s_next
sub delta-E(\s, \u, \v) {
    my (\a,   \b,  \c,  \d) = D[s[u-1];s[u]], D[s[u+1];s[u]], D[s[v-1];s[v]], D[s[v+1];s[v]];
    my (\na, \nb, \nc, \nd) = D[s[u-1];s[v]], D[s[u+1];s[v]], D[s[v-1];s[u]], D[s[v+1];s[u]];
    if    v == u+1 { return (na + nd) - (a + d) }
    elsif u == v+1 { return (nc + nb) - (c + b) }
    else           { return (na + nb + nc + nd) - (a + b + c + d) }
}

# E(s0), initial state
my \s = @ = flat 0, (1 ..^ cities).pick(*), 0;
my \E-min-global = my \E-min = $ = Es(s);

for 0 ..^ kmax -> \k {
    printf "k:%8u  T:%4.1f  Es: %3.1f\n" , k, T(k, kmax, kT), Es(s)
            if k % (kmax/10) == 0;

    # valid candidate cities (exist, adjacent)
    my \cv = neighbors[(^8).roll] + s[ my \u = 1 + (^(cities-1)).roll ];
    next if cv ≤ 0 or cv ≥ cities or D[s[u];cv] > sqrt(2);

    my \v  = s[cv];
    my \ΔE = delta-E(s, u, v);
    if ΔE < 0 or P(ΔE, k, kmax, kT) ≥ rand { # always move if negative
        (s[u], s[v]) = (s[v], s[u]);
        E-min += ΔE;
        E-min-global = E-min if E-min < E-min-global;
    }
}

say "\nE(s_final): " ~ E-min-global.fmt('%.1f');
say "Path:\n" ~ s».fmt('%2d').rotor(20,:partial).join: "\n";

```

Output:

```
k:       0  T: 1.0  Es: 522.0
k:  100000  T: 0.9  Es: 185.3
k:  200000  T: 0.8  Es: 166.1
k:  300000  T: 0.7  Es: 174.7
k:  400000  T: 0.6  Es: 146.9
k:  500000  T: 0.5  Es: 140.2
k:  600000  T: 0.4  Es: 127.5
k:  700000  T: 0.3  Es: 115.9
k:  800000  T: 0.2  Es: 111.9
k:  900000  T: 0.1  Es: 109.4

E(s_final): 109.4
Path:
 0 10 20 30 40 50 60 84 85 86 96 97 87 88 98 99 89 79 78 77
67 68 69 59 58 57 56 66 76 95 94 93 92 91 90 80 70 81 82 83
73 72 71 62 63 64 74 75 65 55 54 53 52 61 51 41 31 21 22 32
42 43 44 45 46 35 34 24 23 33 25 15 16 26 36 47 37 27 17 18
28 38 48 49 39 29 19  9  8  7  6  5  4 14 13 12 11  2  3  1
 0
```