Click here to Skip to main content
15,867,453 members
Articles / General Programming / Algorithms

Solving Steiner problem with the help of steel annealing algorithm.

Rate me:
Please Sign up or sign in to vote.
4.88/5 (12 votes)
12 Aug 2012CPOL3 min read 25.1K   48   21   12
Steiner problem solution.

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> 

C#
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> 

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


Written By
Software Developer (Senior) Oracle
United States United States
This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.

Comments and Discussions

 
QuestionImages in the article Pin
Member 884486021-Jun-13 5:46
Member 884486021-Jun-13 5:46 
QuestionSource download broken? Pin
Kean-W6-Nov-12 8:28
Kean-W6-Nov-12 8:28 
Questionlitle question Pin
SergeyAB6-Aug-12 19:56
SergeyAB6-Aug-12 19:56 
AnswerRe: litle question Pin
Sergey Vystoropskiy6-Aug-12 20:31
Sergey Vystoropskiy6-Aug-12 20:31 
GeneralRe: litle question Pin
SergeyAB6-Aug-12 22:42
SergeyAB6-Aug-12 22:42 
GeneralRe: litle question Pin
Sergey Vystoropskiy6-Aug-12 23:27
Sergey Vystoropskiy6-Aug-12 23:27 
GeneralMy vote of 5 Pin
Kanasz Robert5-Aug-12 9:51
professionalKanasz Robert5-Aug-12 9:51 
outstanding. This article is well written and explained. I hope you will write more articles like that.
You've got my 5

Good luck
Robert
GeneralRe: My vote of 5 Pin
Sergey Vystoropskiy5-Aug-12 10:01
Sergey Vystoropskiy5-Aug-12 10:01 
GeneralWell done! Pin
Espen Harlinn5-Aug-12 2:41
professionalEspen Harlinn5-Aug-12 2:41 
GeneralRe: Well done! Pin
Sergey Vystoropskiy5-Aug-12 3:24
Sergey Vystoropskiy5-Aug-12 3:24 
GeneralRe: Well done! Pin
Espen Harlinn5-Aug-12 3:35
professionalEspen Harlinn5-Aug-12 3:35 
GeneralRe: Well done! Pin
Sergey Vystoropskiy5-Aug-12 3:40
Sergey Vystoropskiy5-Aug-12 3:40 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.