Solving Steiner problem with the help of steel annealing algorithm.

Introduction

First of all lets discus what is Steiner problem. We have a plane with n points. And we can put more m points anywhere we want. After that we should build minimum astringent tree. The problem in this task is to put m points in such way that will allow to get a smallest astringent tree.<o:p>

So why don’t we use a simple brute force algorithm? Because it will be to slow.<o:p> 

And what about hill climbing ? Unfortunately  the function that we are trying to optimize has a to much local minimums so hill climbing will stuck on one of then.<o:p>

Steel annealing algorithm is a one of the heuristic algorithms that is used to find a global minimum. The main Idea of anneal algorithm is the same like hill climbing algorithm: if new solution is better than the solution we have accept it. But anneal algorithm sometimes accept the solution which is not the best. Such “stupid” decisions allow escaping a local minimum trap.<o:p>

For example let’s take a look on the following function:

<v:shapetype id="_x0000_t75" coordsize="21600,21600" o:spt="75" o:preferrelative="t" path="m@4@5l@4@11@9@11@9@5xe" filled="f" stroked="f"> <v:stroke joinstyle="miter"> <v:formulas> <v:f eqn="if lineDrawn pixelLineWidth 0"> <v:f eqn="sum @0 1 0"> <v:f eqn="sum 0 0 @1"> <v:f eqn="prod @2 1 2"> <v:f eqn="prod @3 21600 pixelWidth"> <v:f eqn="prod @3 21600 pixelHeight"> <v:f eqn="sum @0 0 1"> <v:f eqn="prod @6 1 2"> <v:f eqn="prod @7 21600 pixelWidth"> <v:f eqn="sum @8 21600 0"> <v:f eqn="prod @7 21600 pixelHeight"> <v:f eqn="sum @10 21600 0"> <v:path o:extrusionok="f" gradientshapeok="t" o:connecttype="rect"> <o:lock v:ext="edit" aspectratio="t"> <v:shape id="Picture_x0020_1" o:spid="_x0000_i1026" type="#_x0000_t75" style="width: 467.25pt; height: 336pt; visibility: visible; "> <v:imagedata src="file:///C:\Users\Sergey\AppData\Local\Temp\msohtmlclip1\01\clip_image001.png"> <o:p>

If will start somewhere in a red point: 

<v:shape id="Picture_x0020_2" o:spid="_x0000_i1025" type="#_x0000_t75" style="width: 467.25pt; height: 336pt; visibility: visible; "> <v:imagedata src="file:///C:\Users\Sergey\AppData\Local\Temp\msohtmlclip1\01\clip_image002.png"> <o:p>

We can stuck in a point (0,0). In order to get to a real minimum we will need to move to the position which is larger than a local minimum.<o:p>

So in this article I am going to describe how did I implement the solution for the Steiner problem using annealing algorithm. As a programming language I used C#.<o:p> 

Implementation of  astringent tree length calculation. 

First of all lets define the function which we are going to minimize. The function we are going to minimize will be a sum of length of all edges from the astringent tree. In order to build a astringent tree we will use a Prim’s algorithm.<o:p>

The idea of Prim’s algorithm is rather simple:

1) Take any point from the tree and add it to the set of selected points.

2)  Add a most light edge that connects selected and unselected point. Add unselected point to the selected point set.

3) If the unselected points set is not empty go to 2. 

4) Finish. 

Here is a code that illustrates my implementation of that algorithm.<o:p> 

Pair FindTheSmallestEdgeAndInitializeSets(HashSet<int> unUsed, HashSet<int> used, int problemNumber)
<pre>{
    double bestValue = double.MaxValue;
    Pair result = new Pair();
    for (int i = 0; i < data.points[problemNumber].Length; ++i)
    {
        for (int j = 0; j < data.points[problemNumber].Length; ++j)
        {
            if (bestValue > distances_[i, j] && (i != j))
            {
                bestValue = distances_[i, j];
                result.first = i;
                result.second = j;
            }
        }
    }
    used.Add(result.first);
    used.Add(result.second);
    unUsed.Remove(result.first);
    unUsed.Remove(result.second);
    return result;
}
private Pair FindBestEdgeToAdd(HashSet<int> unUsed, HashSet<int> used, int problemNumber)
{
    double bestValue = double.MaxValue;
    Pair result = new Pair();
    foreach (int i in used)
    {
        foreach (int j in unUsed)
        {
            if (bestValue > distances_[i, j])
            {
                bestValue = distances_[i, j];
                result.first = i;
                result.second = j;
            }
        }
    }
    unUsed.Remove(result.second);
    used.Add(result.second);
    return result;
}
private void PrimeAlgorithm(int problemNumber)
{
    HashSet<int> usedPoints = new HashSet<int>();
    HashSet<int> unUsedPoints = new HashSet<int>();
    initializeSet(unUsedPoints, data.points[problemNumber].Length);
    int edgesCount = 0;
    edges = new Pair[data.points[problemNumber].Length - 1];
    edges[edgesCount++] = FindTheSmallestEdgeAndInitializeSets(unUsedPoints, usedPoints, problemNumber);
    while (unUsedPoints.Count != 0)
    {
        edges[edgesCount++] = FindBestEdgeToAdd(unUsedPoints, usedPoints, problemNumber); 
    }
} 
		A brief description of how to use the article or code. The
		class names, the methods and properties, any tricks or tips.

Annealing algorithm implementation.<o:p>

As it was mentioned before we are not always take the best result. Let’s understand when do we select a bad solution. The name “anneal algorithm” was taken because in this algorithm we use the idea of cooling of steel. When the metal is cooling it’s particles is cooling to, but sometimes the temperature of a single particle is increasing. The higher is the metal temperature the higher is the probability for the particle to become hotter. <o:p>

So we will have some initial temperature. After each iteration the temperature will be less and less. So the probability of acceptance of a bad result will be lover.<o:p>

Here is my implementation of annealing algorithm:

private void Annealing(int problemNumber)
{
    double alpha = 0.99;
    double temperature = 100.0;
    double initialTemperature = temperature;
    double epsilon = 0.2;
    double deltha;
    double accepted = 0;
    double rejected = 0;
    double multiplicator = 0.05;
    best = GetTreeLength(problemNumber);
    Point[] bestPoints = new Point[data.points[problemNumber].Length];
    Pair[] bestEdges = new Pair[edges.Length];
    Array.Copy(data.points[problemNumber], bestPoints, bestPoints.Length);
    Array.Copy(edges, bestEdges, edges.Length);
    Random rand = new Random();
    
    while (temperature > epsilon)
    {
        double prevLength = GetTreeLength(problemNumber);
        double[,] offsets = new double[2, mutablePointsNumber]; 
        for (int i = 0; i < mutablePointsNumber; ++i)
        {
            offsets[0,i] = rand.NextDouble() * multiplicator - 0.5 * multiplicator;
            offsets[1,i] = rand.NextDouble() * multiplicator - 0.5 * multiplicator;
            int currentIndex = data.points[problemNumber].Length - i - 1;
            ModifyPoint(currentIndex, problemNumber, offsets[0, i], offsets[1, i]);
        }
            
            
        double currentLength = GetTreeLength(problemNumber);
        double proba = rand.NextDouble();
        deltha = currentLength - prevLength;
        double e = Math.Exp((-deltha * 13 * initialTemperature) / temperature);
        if (proba < e || deltha < 0)
        {
            accepted++;
            if (currentLength < best)
            {
                best = currentLength;
                Array.Copy(data.points[problemNumber], bestPoints, bestPoints.Length);
                Array.Copy(edges, bestEdges, edges.Length);
            }
        }
        else
        {
            rejected++;
            for (int i = 0; i < mutablePointsNumber; ++i)
            {
                int currentIndex = data.points[problemNumber].Length - i - 1;
                ModifyPoint(currentIndex, problemNumber, -offsets[0, i], -offsets[1, i]);
            }
        }
        
        chart1.Series[0].Points.AddY(currentLength);
        if (accepted > mutablePointsNumber * 2  || rejected > 4 * mutablePointsNumber )
        {
            accepted = 0;
            rejected = 0;
            temperature *= alpha;
        }
    }
    Array.Copy(bestPoints, data.points[problemNumber], bestPoints.Length);
    Array.Copy(bestEdges, edges, edges.Length);
} 

Results.<o:p>

To illustrate results I had used the following graphics<o:p> 

<v:shapetype id="_x0000_t75" coordsize="21600,21600" o:spt="75" o:preferrelative="t" path="m@4@5l@4@11@9@11@9@5xe" filled="f" stroked="f"> <v:stroke joinstyle="miter"> <v:formulas> <v:f eqn="if lineDrawn pixelLineWidth 0"> <v:f eqn="sum @0 1 0"> <v:f eqn="sum 0 0 @1"> <v:f eqn="prod @2 1 2"> <v:f eqn="prod @3 21600 pixelWidth"> <v:f eqn="prod @3 21600 pixelHeight"> <v:f eqn="sum @0 0 1"> <v:f eqn="prod @6 1 2"> <v:f eqn="prod @7 21600 pixelWidth"> <v:f eqn="sum @8 21600 0"> <v:f eqn="prod @7 21600 pixelHeight"> <v:f eqn="sum @10 21600 0"> <v:path o:extrusionok="f" gradientshapeok="t" o:connecttype="rect"> <o:lock v:ext="edit" aspectratio="t"> <v:shape id="Picture_x0020_5" o:spid="_x0000_i1026" type="#_x0000_t75" style="width: 468pt; height: 369pt; visibility: visible; "> <v:imagedata src="file:///C:\Users\Sergey\AppData\Local\Temp\msohtmlclip1\01\clip_image001.png"> <o:p>

This one illustrates how do the length changes during the algorithm iterations. And in the right bottom corner we can see the difference between initial length and a result length.<o:p>

To illustrate the positions of the points I made a following visualization. <o:p>

<v:shape id="Picture_x0020_6" o:spid="_x0000_i1025" type="#_x0000_t75" style="width: 468pt; height: 369pt; visibility: visible; "> <v:imagedata src="file:///C:\Users\Sergey\AppData\Local\Temp\msohtmlclip1\01\clip_image002.png"> <o:p>

Code and test data of the program is Download source.<o:p>