Monday, March 6, 2017

Simulated Annealing in C#

After I explained genetic algorithms in a previous post, I today want to present another heuristic which is inspired by nature: Simulated annealing.
Prototype is the cooling of materials, in which the atoms move freely and a lot in the beginning, when temperature is high (in the beginning of the optimization we traverse big parts of the search space). As time goes on the movement of the atoms slows down, the system reaches a stable condition (towards the end of the search we only look for local improvements of our solution and end up in a local minimum).
This procedure is for example used for designing integrated circuits.
In contrast to genetic algorithms in simulated annealing we do not work with a pool of solutions, but just with a single one. A pool of solutions has its advantages, but slows down the process of course. The creation of a new solution can be compared with asexual reproduction: Based on our current solution we modify this to reach a new solution. For this we define a graph, which represents the neighbor relation of solutions, and traverse this.  The heavy movement of the atoms for high temperatures and the resulting scanning of many solutions we model with an acceptance probability: If the new solution is better than the current one, we always accept this as the new current solution. If the solution is worse though, we only accept it with a probability of e-Δ/T.
Here Δ is the difference of solution values (result of new solution - result of current solution) and T the current temperature. That means, the probability to accept a new solution decreases with the temperature and with the worsening amount of the solution value (just plug in some values to see that).
In every step of the algorithm the temperature is reduced by some scheme, here we choose the geometric cooling schedule: TNew = TOld * α.
In pseudocode the algorithm looks as follows:
Current_Solution = random solution
T = TStart
While T > TEnd:
   Temp = Randon neighbored solution of Current_Solution
   If Temp better than Current_Solution:
      Current_Solution = Temp
   Else:
      Current_Solution = Temp with probability mentioned above
   T = T * α. 

Let us now apply simulated annealing to the TSP - problem, which was already presented in the post about genetic algorithms. In our neighborhood graph two solutions are adjacent, if they can be created out of each other by swapping two cities in the roundtrip. As an example: The solution ABC is a neighbor of the solution BAC, CAB and ACB.
At the beginning of the algorithm we create a random initial solution in the function Initialize(). For this we try out as long as we find a valid solution, meaning one in which all cities are connected with each other. This is neccessary, since the algorithm can hardly improve from a bad solution with costs infinity by swapping just two cities.
In the function Go() we see the main loop of the algorithm. In every iteration with the function CreateRandomNeighbor() a random new solution is created based on the current one.
For this, as mentioned, two random cities are swapped.
The decision, whether the new solution is accepted, is made upon the criterion mentioend above.
As a starting temperature I here chose 1000, as an ending temperature 0.01, for α 0.999999.
With this parameters the algorithm reached for 100 random instances of size 20 a score of 4516 on average.
Soon I will write an own post about the comparison of genetic algorithms and simulated annealing and about the selection of parameters, but in a first comparison this gives simulted annealing the edge.
This also confirms my opinion, I personally prefer the elegant, easy but effective idea of simulated annealing, also the speed is superior.
For genetic algorithms it is further not clear, if the principle of sexual reproduction transfers to mathematical problems in a meaningful way.

And here is the code:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace SA_TSP
{
    class Program
    {
        static void Main(string[] args)
        {
            TSPInstance Instance = CreateRandomTSPInstance(20);
            SimulatedAnnealing SA = new SimulatedAnnealing();
            SA.Go(Instance, true, 1000, 0.01, 0.999999);
        }

        // Create a random TSP instance, connect two cities with probability 70 % and choose a length between 0 and 1000.
        public static TSPInstance CreateRandomTSPInstance(int n)
        {
            TSPInstance Instance = new TSPInstance(n);
            Random Rnd = new Random();
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    if (Rnd.Next(100) < 70)
                    {
                        Instance.AddConnection(i, j, Rnd.Next(1000));
                    }
                }
            }
            return Instance;
        }

    
    }

    public class TSPInstance
    {
        public int n;
        public int[][] Connections;

        public TSPInstance(int _n)
        {
            n = _n;
            Connections = new int[n][];
            for (int i = 0; i < n; i++)
            {
                Connections[i] = new int[n];
                for (int j = 0; j < n; j++)
                {
                    Connections[i][j] = int.MaxValue;
                }
            }
        }

        public void AddConnection(int start, int end, int distance)
        {
            Connections[start][end] = distance;
        }
    }

    public class SimulatedAnnealing
    {
        public class Individual : ICloneable
        {
            public double Fitness;
            public int[] Tour;

            public Individual(int length)
            {
                Tour = new int[length];
            }

            // creates a deep copy of an individual
            public object Clone()
            {
                Individual Cloned = new Individual(Tour.Length);

                for (int i = 0; i < Tour.Length; i++)
                {
                    Cloned.Tour[i] = Tour[i];
                }
                Cloned.Fitness = Fitness;
                return Cloned;
            }
        }

        Individual CurrentSolution;
        TSPInstance Instance;

        public void Initialize()
        {
            // fill the initial population, for each individual create a random permutation of the cities

            Random Rnd = new Random();
            double Fitness = int.MaxValue;

            // create random permutations, as long as a valid tour is found
            while (Fitness == int.MaxValue)
            {
                List<int> Cities = new List<int>();
                for (int j = 0; j < Instance.n; j++)
                {
                    Cities.Add(j);
                }
                int Counter = 0;
                while (Cities.Count > 0)
                {
                    int Index = Rnd.Next(Cities.Count);
                    CurrentSolution.Tour[Counter++] = Cities[Index];
                    Cities.RemoveAt(Index);
                }
                CurrentSolution.Fitness = CalculateFitness(CurrentSolution);
                Fitness = CurrentSolution.Fitness;
            }
        }

        public Individual CreateRandomNeighbor()
        {
            // switch to random cities in the tour
            Random Rnd = new Random();
            Individual Cloned = (Individual)CurrentSolution.Clone();
            int End = Rnd.Next(Instance.n);
            int Start = Rnd.Next(Instance.n);
            int Temp = Cloned.Tour[Start];
            Cloned.Tour[Start] = Cloned.Tour[End];
            Cloned.Tour[End] = Temp;
            Cloned.Fitness = CalculateFitness(Cloned);
            return Cloned;
        }

        public Individual Go(TSPInstance inst, Boolean print, double tstart, double tend, double alpha)
        {
            Instance = inst;
            double T = tstart;
            CurrentSolution = new Individual(inst.n);

            Initialize();

            double LastPrint = tstart;

            // repeat as long as the temperature is over some threshold
            while (T > tend)
            {
                Individual NewSolution = CreateRandomNeighbor();
                if (NewSolution.Fitness < CurrentSolution.Fitness)
                    CurrentSolution = NewSolution;
                else
                {
                    double AcceptanceProbability = Math.Exp(-(NewSolution.Fitness - CurrentSolution.Fitness) / T);
                    Random Rnd = new Random();
                    if (Rnd.NextDouble() <= AcceptanceProbability)
                        CurrentSolution = NewSolution;
                }
                T = T * alpha;
                if (print && (LastPrint - T > 50))
                {
                    Console.WriteLine("Temperature: " + T.ToString());
                    Console.WriteLine("Current Solution: " + CurrentSolution.Fitness.ToString());
                    LastPrint = T;
                }
            }

            if (print)
                Console.WriteLine("Best: " + CurrentSolution.Fitness.ToString());

            return CurrentSolution;
        }

        public int CalculateFitness(Individual ind)
        {
            // sum over entire tour and return costs
            int Costs = 0;
            for (int i = 0; i < ind.Tour.Length; i++)
            {
                if (i == ind.Tour.Length - 1)
                {
                    if (Instance.Connections[ind.Tour[i]][ind.Tour[0]] == int.MaxValue)
                        return int.MaxValue;
                    Costs += Instance.Connections[ind.Tour[i]][ind.Tour[0]];
                }
                else
                {
                    if (Instance.Connections[ind.Tour[i]][ind.Tour[i + 1]] == int.MaxValue)
                        return int.MaxValue;
                    Costs += Instance.Connections[ind.Tour[i]][ind.Tour[i + 1]];
                }
            }
            return Costs;
        }
    }
}