What is a simulated annealing algorithm?

Simulated annealing (SA) is a stochastic local search algorithm inspired by the annealing process in metallurgy. It is an effective and widely used technique for approximating the global optimum of various optimization problems. Similar to a local search algorithm, it begins with an initial solution to the problem and then iteratively explores neighboring solutions by making small modifications to the current solution. However, unlike hill climbing local search algorithm simulated annealing allows occasional worse transitions. This ability to make worse transitions helps overcome the issue of getting stuck in a local optimum.

Simulated annealing introduces the concept of a temperature parameter, which controls the probability of accepting worse solutions. During the iterative process, the algorithm generates a neighboring solution and evaluates its quality based on a cost or objective function.

If the new solution is better than the current solution, it is accepted as the new current solution. However, if the new solution is worse, it may still be accepted with a certain probability. The probability of accepting a worse solution is determined by a formula that depends on the temperature and the difference in cost between the new and current solutions. The formula is designed in such a way that the probability of accepting a worse solution is higher when the temperature is higher. As the algorithm progresses, the temperature gradually decreases. As iterations proceed, it becomes less likely to transition to a worse solution, and the process stabilizes.

Simulated annealing continues iterating and exploring the search space until either a stopping criterion is met (such as reaching a maximum number of iterations or reaching the final temperature) or no further improvements are observed. The best solution obtained during the iterations is an approximation of the optimal solution to the problem.

In this Answer, we will delve into the steps of the algorithm and examine its pseudocode. Following that, we will apply this technique to a well-known and challenging optimization problem which is the traveling salesman problem.

Detailed description

Here is the formal description of the problem:

Let $\mathcal{S}$ be the set of all possible solutions of a given optimization problem. Let $S \in \mathcal{S}$ be a specific solution. The cost function $f: \mathcal{S} \to \mathbb{R}$ assigns a real number as the cost value to each solution, the value $f(S)$ represents the cost of solution $S$ . Our goal is to minimize the cost function $f(S)$ over all possible solutions $\mathcal{S}$ . This can be formulated as to find minimum value: $\min \limits_{S \in \mathcal{S}} f(S)$ .

If we want to maximize $f(S)$ , we can consider the problem as a minimization problem by negating $f(S)$ function: $\max \limits_{S \in \mathcal{S}}f(S)=\min \limits_{S \in \mathcal{S}}(-f(S))$ .

Here are the steps of the algorithm:

Initialization: Start with an initial solution to the problem. This can be a random solution or a solution obtained using some method.
Define parameters: Set the initial temperature and cooling schedule. Typically, initial temperature $T_0=100$ and cooling are done by formula $T_{new} = \alpha T_{current}$ , where $\alpha \in [0.7, 0.99]$ is a cooling factor. Set number of iterations per temperature: typically, this number is between $100$ and $1000$ .
Iterative process: Perform iterations until a stopping criterion is met. This can be a maximum number of iterations or reaching a final temperature.
Generate a neighboring solution: Modify the current solution to generate a neighboring solution. The modification can be swapping variables or changing values.
Evaluate the neighboring solution: Calculate the cost or objective function value for the neighboring solution, indicating its quality. The cost value $f(S_{neighbor})$ should be evaluated for neighbor solution $S_{neighbor}$ .
Acceptance criterion: Determine whether to accept the neighboring solution as the new current solution. If the neighboring solution is better than the current solution, accept it: $S_{new}= S_{neighbor}$ if $f(S_{current}) > f(S_{neighbor})$ . If the neighboring solution is worse, accept it with a probability determined by the acceptance probability formula: $S_{new} = S_{neighbor}$ if $\exp\{{-(f(S_{neighbor})-f(S_{current}))}\} > r$ , where $r \in [0, 1]$ is some random number.
Update the temperature: Adjust the temperature according to the cooling schedule: $T_{new} = \alpha T_{current}$ .
Repeat: Go back to step four and continue the iterative process.
Termination: Stop the algorithm when the stopping criterion is met.
Output: Return the best solution (with the lowest cost value) obtained during the iterations. This solution approximates the global optimum solution.

Solution SimulatedAnnealing(ProblemInstance problemInstance) 
{
    //Get initial random solution
    Solution currentSolution = InitSolution(problemInstance);
    //Initialize parameters
    double T = InitTemperature();
    const double coolingFactor = InitCoolingFactor();
    const int stepsPerT = InitStepsPerTemperature();
    Solution bestSolution = solution;
    while(!StopCriteria()) 
    {
        for(int i = 0; i < stepsPerT; ++i) 
        {
            //Generate neighbor solution
            Solution nextSolution = GetNeighbor(currentSolution);
            // Calculate difference between cost values of new and current solutions
            double delta = nextSolution.GetCostValue() - currentSolution.GetCostValue();
            // Acceptance criterion
            if(delta < 0 || exp(-delta / T) > GetRandom(0, 1))
            {
                currentSolution = nextSolution;
                if(currentSolution.GetCostValue() < bestSolution.GetCostValue())
                {
                    bestSolution = currentSolution;
                }
            }
        }
        T *= coolingFactor;
    }
    return solution;
}

Applications

We will focus on applications of simulated annealing algorithm for the traveling salesman problem (TSP):

"Given a list of cities and the distances between each pair, what is the shortest possible route that visits each city exactly once and returns to the origin city?"

The set $\mathcal{S}$ of all possible solutions can be described as the routes that pass through all cities and return to the origin. If $n$ is the number of cities, the solution space can be described as $(n-1)!$ circular permutations of the list of cities. Therefore, a solution $S = (S_1, ...., S_n), \text{where }S_i \in \{1, \dots, n\},$ is a circular permutation of $\{1, \dots, n\}$ that corresponds to one route. Let $S = (S_1, ...., S_n)$ be a solution, the cost function $f(S)$ is given by: $f(S) = \sum_{i=1}^{n}d(S_{i}, S_{i+1})$ , where $d(S_{i}, S_{i+1})$ is the distance between cities $S_{i}$ and $S_{i+1}$ and $S_{n+1}=S_1$ .

#include <iostream>
#include <vector>
#include <algorithm>
#include <random>
#include <cmath>
using namespace std;
class TSPSolution {
private:
    vector<size_t> cities_; // Order in which cities are visited
    double cost_; // Total cost of the solution
    vector<vector<double>> distances_; // Matrix distances between cities
public:
    // Construct a random solution to the TSP problem
    TSPSolution(const vector<vector<double>>& distances):
    cost_(0.0),
    distances_(distances)
    {
        cities_.resize(distances_.size());
        for(size_t i = 0; i < distances_.size(); ++i) 
        {
            cities_[i] = i;
        }
        // Shuffle indices except the first one, as we assume that routes start from the city with index 0
        shuffle(cities_.begin() + 1, cities_.end(), mt19937(random_device{}()));
        CalculateCost();
    }
    // Get the cost of the solution
    double GetCost() const 
    {
        return cost_;
    }
    // Get the order in which cities are visited
    const vector<size_t>& GetCities() const 
    {
        return cities_;
    }
    // Calculate the cost of the solution given a matrix of distances between cities
    void CalculateCost() 
    {
        cost_ = 0.0;
        for(size_t i = 0; i < cities_.size() - 1; ++i) 
        {
            cost_ += distances_[cities_[i]][cities_[i+1]];
        }
        cost_ += distances_[cities_[cities_.size()-1]][cities_[0]];
    }
    // Swap two cities in the route and evaluate the cost
    void SwapCities(size_t i, size_t j) 
    {
        size_t size  = cities_.size();
        size_t prevI = (i + size - 1) % size;
        size_t nextI = (i + size + 1) % size;
        size_t prevJ = (j + size - 1) % size;
        size_t nextJ = (j + size + 1) % size;
        cost_ -= distances_[cities_[prevI]][cities_[i]];
        cost_ -= distances_[cities_[i]][cities_[nextI]];
        cost_ -= distances_[cities_[prevJ]][cities_[j]];
        cost_ -= distances_[cities_[j]][cities_[nextJ]];
        swap(cities_[i], cities_[j]);
        cost_ += distances_[cities_[prevI]][cities_[i]];
        cost_ += distances_[cities_[i]][cities_[nextI]];
        cost_ += distances_[cities_[prevJ]][cities_[j]];
        cost_ += distances_[cities_[j]][cities_[nextJ]];
    }
};
struct SAParameters {
    double initialTemperature;
    double finalTemperature;
    double coolingFactor;
    int iterationsPerTemperature;
};
// Function to generate a random integer within a range
size_t RandomInt(size_t min, size_t max) 
{
    static random_device rd;
    static mt19937 generator{rd()};
    uniform_int_distribution<size_t> distribution{min, max};
    return distribution(generator);
}
TSPSolution SimulatedAnnealing(const vector<vector<double>>& distances, const SAParameters& parameters) 
{
    // Create a random initial solution
    TSPSolution currentSolution(distances);
    TSPSolution bestSolution = currentSolution;
    double currentTemperature = parameters.initialTemperature;
    while(currentTemperature > parameters.finalTemperature) 
    {
        for(int i = 0; i < parameters.iterationsPerTemperature; ++i) 
        {
            // Generate a new neighbor solution by swapping two random cities
            TSPSolution newSolution(currentSolution);
            size_t cityIndex1 = RandomInt(1, distances.size() - 1);
            size_t cityIndex2 = RandomInt(1, distances.size() - 1);
            newSolution.SwapCities(cityIndex1, cityIndex2);
            // Calculate difference between costs of new and current solutions
            double delta = newSolution.GetCost() - currentSolution.GetCost();
            // Accept the new solution if it has a lower cost or according to the acceptance probability
            if (delta < 0 || exp(-delta / currentTemperature) > RandomInt(0, 1000) / 1000.0) 
            {
                currentSolution = newSolution;
                if(currentSolution.GetCost() < bestSolution.GetCost()) 
                {
                    bestSolution = newSolution;
                }
            }
        }
        // Cool down
        currentTemperature *= parameters.coolingFactor;
    }
    return bestSolution;
}
int main() 
{
    // Define the distances between cities
    vector<vector<double>> distances = {
        {0, 3, 7, 2, 9},
        {3, 0, 4, 1, 8},
        {7, 4, 0, 6, 5},
        {2, 1, 6, 0, 4},
        {9, 8, 5, 4, 0}
    };
    // Set parameters for simulated annealing algorithm
    SAParameters parameters = {
        .initialTemperature = 1000.0,
        .finalTemperature = 0.1,
        .coolingFactor = 0.9,
        .iterationsPerTemperature = 1000
    };
    // Run simulated annealing algorithm
    TSPSolution solution = SimulatedAnnealing(distances, parameters);
    cout << "Shortest circular route found: ";
    const vector<size_t>& route = solution.GetCities();
    for (size_t city : route) 
    {
        cout << city + 1 << " ";
    }
    cout << route[0] + 1 << endl;
    cout << "Total cost: " << solution.GetCost() << endl;
    return 0;
}

The TSPSolution class represents a route consisting of a sequence of cities, where the order is a circular permutation of indices $\{0, ..., n-1\}$ starting from index 0. The SAParameters structure defines the parameters that control the simulated annealing algorithm. These parameters include initial temperature, cooling factor, number of iterations per temperature, and final temperature. These parameters determine the behavior and convergence of the algorithm.

The SimulatedAnnealing() function is the implementation of the simulated annealing algorithm. It accepts the distances between cities and the specific control parameters through the SAParameters structure. In the main() function algorithm is executed on an example with five cities. The optimal solution/route for this specific problem instance is returned.

In conclusion, simulated annealing is a powerful optimization technique particularly useful in complex problems with multiple local optima. By leveraging randomness and the concept of temperature, it explores the solution space effectively, allowing it to escape local optima and potentially find high-quality solutions.

Free Resources