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>