|
| 1 | +--- |
| 2 | +tags: |
| 3 | + - Original |
| 4 | +--- |
| 5 | + |
| 6 | +# Simulated Annealing |
| 7 | + |
| 8 | +**Simulated Annealing (SA)** is a randomized algorithm, which approximates the global optimum of a function. It's called a randomized algorithm, because it employs a certain amount of randomness in its search and thus its output can vary for the same input. |
| 9 | + |
| 10 | +## The Problem |
| 11 | + |
| 12 | +We are given a function $ E(s) $, which calculates the potential of the state $ s $. We are tasked with finding the state $ s_{best} $ at which $ E(s) $ is minimized. SA is suited for problems where the states are discrete and $ E(s) $ has multiple local minima. We'll take the example of the [Travelling Salesman Problem (TSP)](https://en.wikipedia.org/wiki/Travelling_salesman_problem). |
| 13 | + |
| 14 | +### Travelling Salesman Problem (TSP) |
| 15 | + |
| 16 | +You are given a set of nodes in 2 dimensional space. Each node is characterised by its $ x $ and $ y $ coordinates. Your task is to find the an ordering of the nodes, which will minimise the distance to be travelled when visiting these nodes in that order. |
| 17 | + |
| 18 | +### State |
| 19 | + |
| 20 | +State space is the collection of all possible values that can be taken by the independent variables. |
| 21 | +A State is a unique point in the state space of the problem. In the case of TSP, all possible paths that we can take to visit all the nodes is the state space, and any single one of these paths can be considered as a State. |
| 22 | + |
| 23 | +### Neighbouring State |
| 24 | + |
| 25 | +It is a State in the State space which is close to the previous State. This usually means that we can obtain the neighbouring State from the original State using a simple transform. In the case of the Travelling salesman problem, a neighbouring state is obtained by randomly choosing 2 nodes, and swapping their positions in the current state. |
| 26 | + |
| 27 | +### E(s) |
| 28 | + |
| 29 | +$ E(s) $ is the function which needs to be minimised (or maximised). It maps every state to a real number. In the case of TSP, $ E(s) $ returns the distance of travelling one full circle in the order of nodes in the state. |
| 30 | + |
| 31 | +## Approach |
| 32 | + |
| 33 | +We start of with a random state $ s $. In every step, we choose a neighbouring state $ s_{next} $ of the current state $ s $. If $ E(s_{next}) < E(s) $, then we update $ s = s_{next} $. Otherwise, we use a probability acceptance function $ P(E(s),E(s_{next}),T) $ which decides whether we should move to $ s_{next} $ or stay at s. T here is the temperature, which is initially set to a high value and decays slowly with every step. The higher the temperature, the more likely it is to move to s_next. |
| 34 | +At the same time we also keep a track of the best state $ s_{best} $ across all iterations. Proceed till convergence or time runs out. |
| 35 | + |
| 36 | + |
| 37 | +## Probability Acceptance Function |
| 38 | +$$ |
| 39 | +P(E,E_{next},T) = |
| 40 | + \begin{cases} |
| 41 | + \text{True} &\quad\text{if } \exp{\frac{-(E_{next}-E)}{T}} \ge random(0,1) \\ |
| 42 | + \text{False} &\quad\text{otherwise}\\ |
| 43 | + \end{cases} |
| 44 | +
|
| 45 | +$$ |
| 46 | +This function takes in the current state, the next state and the Temperature , returning a boolean value, which tells our search whether it should move to $ s_{next} $ or stay at s. Note that for $E_{next} < E$ , this function will always return True. |
| 47 | + |
| 48 | +```cpp |
| 49 | +bool P(double E,double E_next,double T){ |
| 50 | + double e = 2.71828; |
| 51 | + double prob = (double)rand()/RAND_MAX; // Generate a random number between 0 and 1 |
| 52 | + if(pow(e,-(E_next-E)/T) > prob) return true; |
| 53 | + else return false; |
| 54 | +} |
| 55 | +``` |
| 56 | +## Code Template |
| 57 | +
|
| 58 | +```cpp |
| 59 | +class state{ |
| 60 | + public: |
| 61 | + state(){ |
| 62 | + // Generate the initial state |
| 63 | + } |
| 64 | + state next(){ |
| 65 | + state next; |
| 66 | + next = s; |
| 67 | + // Make changes to the state "next" and then return it |
| 68 | + return next; |
| 69 | + } |
| 70 | + double E(){ |
| 71 | + // implement the cost function here |
| 72 | + }; |
| 73 | +}; |
| 74 | +
|
| 75 | +int main(){ |
| 76 | + double T = 1000; // Initial temperature |
| 77 | + double u = 0.99; // decay rate |
| 78 | + state s = state(); |
| 79 | + state best = s; |
| 80 | + double E = s.E(); |
| 81 | + double E_next; |
| 82 | + double E_best = E; |
| 83 | + while (T > 1){ |
| 84 | + state next = s.next(); |
| 85 | + E_next = next.E(); |
| 86 | + if(P(E,E_next,T)){ |
| 87 | + s = next; |
| 88 | + if(E_next < E_best){ |
| 89 | + best = s; |
| 90 | + E_best = E_next; |
| 91 | + } |
| 92 | + } |
| 93 | + E = E_next; |
| 94 | + T *= u; |
| 95 | + } |
| 96 | + cout << E_best << "\n"; |
| 97 | + return 0; |
| 98 | +} |
| 99 | +``` |
| 100 | +## How to use: |
| 101 | +Fill in the state class functions as appropriate. If you are trying to find a global maxima and not a minima, ensure that the $ E() $ function returns negative of the function you are maximising and finally print out $ -E_{best} $. Set the below parameters as per your need. |
| 102 | + |
| 103 | +### Parameters |
| 104 | +- T : Temperature. Set it to a higher value if you want the search to run for a longer time |
| 105 | +- u : Decay. Decides the rate of cooling. A slower cooling rate (larger value of u) usually gives better results. Ensure $u < 1$. |
| 106 | + |
| 107 | +The number of iterations the loop will run for is given by the expression |
| 108 | +$$ |
| 109 | +N = \lceil -\log_{u}{T} \rceil |
| 110 | +$$ |
| 111 | + |
| 112 | +To see the effect of decay rate on solution results, run simulated annealing for decay rates 0.95 , 0.97 and 0.99 and see the difference. |
| 113 | + |
| 114 | +### Example State class for TSP |
| 115 | +```cpp |
| 116 | +class state{ |
| 117 | + public: |
| 118 | + vector<pair<int,int>> points; |
| 119 | + |
| 120 | + state(){ // Initial random order of points |
| 121 | + points = {{0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0}}; |
| 122 | + } |
| 123 | + state next(){ // picks 2 random indices and swaps them |
| 124 | + state s_next; |
| 125 | + s_next.points = points; |
| 126 | + int a = ((rand()*points.size())/RAND_MAX); |
| 127 | + int b = ((rand()*points.size())/RAND_MAX); |
| 128 | + pair<int,int> t = s_next.points[a]; |
| 129 | + s_next.points[a] = s_next.points[b]; |
| 130 | + s_next.points[b] = t; |
| 131 | + return s_next; |
| 132 | + } |
| 133 | + |
| 134 | + double euclidean(pair<int,int> a, pair<int,int> b){ // return euclidean distance between 2 points |
| 135 | + return pow(pow((a.first-b.first),2)+pow((a.second-b.second),2),0.5); |
| 136 | + } |
| 137 | + double E(){ // calculates the round cost of travelling one full circle. |
| 138 | + double dist = 0; |
| 139 | + bool first = true; |
| 140 | + int n = points.size(); |
| 141 | + for(int i = 0;i < n; i++){ |
| 142 | + dist += euclidean(points[i],points[(i+1)%n]); |
| 143 | + } |
| 144 | + return dist; |
| 145 | + }; |
| 146 | +}; |
| 147 | +``` |
| 148 | + |
| 149 | +## Extra Modifications to the Algorithm: |
| 150 | + |
| 151 | +- Add a time based exit condition to the while loop to prevent TLE |
| 152 | +- You can replace the e value in the Probability Acceptance function to any real number > 1. For a given $ E_{next} - E > 0 $, a higher e value reduces the chance of accepting that state and a smaller e value, increases it. |
| 153 | + |
| 154 | + |
| 155 | +## Problems |
| 156 | + |
| 157 | +- https://usaco.org/index.php?page=viewproblem2&cpid=698 |
| 158 | +- https://codeforces.com/contest/1556/problem/H |
| 159 | +- https://atcoder.jp/contests/intro-heuristics/tasks/intro_heuristics_a |
0 commit comments