Function Repository Resource:

AntColonyOptimization

Source Notebook

Solve a combinatorial optimization problem using an Ant Colony Optimization (ACO) algorithm

Contributed by: Lars Ulke-Winter

ResourceFunction["AntColonyOptimization"][array,fitness]

returns the minimum found result for the function fitness with possible values for each variable given in array.

ResourceFunction["AntColonyOptimization"][array,fitness,All]

returns the intermediate optimization step results and optimization parameters as a summary Association.

Details and Options

ResourceFunction["AntColonyOptimization"] is a global stochastic optimization method for the solution of combinatorial arrangement problems. It imitates the swarm communication behaviour of ants during foraging (Stigmergic Optimization).
The optimization is based on the exchange of information between individuals (randomly generated solutions) using a virtual pheromone trail for the information exchange.
ResourceFunction["AntColonyOptimization"] finds a solution of vector x, where the ith element of x is chosen from the values of the ith list in array, such that fitness[x] has the smallest positive value found.
The input array takes a List of lists with available values at each row. The lists need not have uniform length.
The fitness function to be minimized must return a positive scalar value. Therefore, the standard method of negating the objective function to convert a minimization problem into a maximization problem doesn’t work in this ACO implementation.
Limitation: This implementation does not consider dependencies of design variables among each other. The assignment of the design variables must be independent of the settings of the others. The definition of the search space is therefore possible using only a given array and does therefore not require a given tree or graph structure of the search space. Optimization problems such as the traveling salesman cannot be encoded in this way.
Optimization parameters given as options include:
"PopulationSize"500number of ants
"MaxIterations"200maximum iterations
"FinishAtPheromoneThreshold"0.95stop criterion
"PheromoneFunction"(1/#)&pheromone concentration to be deposited depending on fitness
"WeatheringFactor"0.05fraction of decreasing pheromone concentration (selection probability) per iteration
"ProbabilitySelectBestPath"0.02proportion of generated solution alternatives that currently prefer the best solution
The computation will stop if the selection probability for a variable state exceeds the value of "FinishAtPheromoneThreshold".
"PheromoneFunction" is a function whose result marks the selection probability for the next step. By default, larger fitness results receive lower marking values than smaller fitness results.
Optimization parameters control the runtime behavior and solution quality as a tradeoff between exploitation and exploration. These can be modified depending on the target function and the number of variables, e.g. by trial and error.
The option "ProbabilitySelectBestPath" takes values between zero and one. Smaller values use probabilistic selecting resulting in slower convergence but exploring more possibilities.
The current minimal selection probability for a variable state is monitored during the iterations. The smaller it is, the greater the ability to explore new solutions.
The output of ResourceFunction["AntColonyOptimization"] is a List representing the best found solution. Due to the stochastic search method, this optimum does not necessarily have to be global, thus better solutions may exist.
When the third argument is All, the output is an Association with keys:
"Best"variable allocation which minimizes fitness
"Iterations"for each iteration step: {arguments which minimize fitness, normalized pheromone matrix, and current minimal selection probability for a variable state}
"NumberOfFunctionCalls"number of calls made to the fitness calls

Examples

Basic Examples (2) 

Create an array of possible assignments:

In[1]:=
SeedRandom[314];
possibValues = Table[RandomSample[Range[10]], {i, Range[10]}];
possibValues // Grid
Out[3]=

Define an objective function to minimize:

In[4]:=
fit[x_List] := Total[x]

Search for minimum:

In[5]:=
ResourceFunction["AntColonyOptimization"][possibValues, fit]
Out[5]=

Return the complete iteration process using the third argument of All:

In[6]:=
SeedRandom[314];
possibValues = Table[RandomSample[Range[10]], {i, Range[10]}];
fit[x_List] := Total[x]
In[7]:=
result = ResourceFunction["AntColonyOptimization"][possibValues, fit, All];
result["Best"]
Out[5]=

Convergence behavior:

In[8]:=
result["NumberOfFunctionCalls"]
Out[8]=
In[9]:=
Column@{ListLinePlot[fit /@ result["Iterations"][[All, 1]], PlotLabel -> "Fitness function", GridLines -> Automatic],
  ListLinePlot[result["Iterations"][[All, 3]], PlotLabel -> "Min. pheromone concentration", GridLines -> Automatic]}
Out[9]=
In[10]:=
ListAnimate[
 ArrayPlot[#, Mesh -> All, PlotLabel -> "Pheromone concentration"] & /@
   Table[result["Iterations"][[i, 2]], {i, Length@result["Iterations"]}]]
Out[10]=

Options (4) 

The above example with default options shown:

In[11]:=
SeedRandom[314];
In[12]:=
ResourceFunction["AntColonyOptimization"][
  Table[RandomSample[Range[10]], {i, Range[10]}], Total[#] &, All,
  "PopulationSize" -> 500, "MaxIterations" -> 200, "FinishAtPheromoneThreshold" -> 0.95, "WeatheringFactor" -> 0.05, "PheromoneFunction" -> (1/# &), "ProbabilitySelectBestPath" -> 0.02];
In[13]:=
%["Best"]
Out[13]=
In[14]:=
%%["NumberOfFunctionCalls"]
Out[14]=

FinishAtPheromoneThreshold (1) 

Terminate if a quarter of individuals have the same fitness value:

In[15]:=
SeedRandom[314];
In[16]:=
ResourceFunction["AntColonyOptimization"][
  Table[RandomSample[Range[10]], {i, Range[10]}], Total[#] &, All,
  "FinishAtPheromoneThreshold" -> 0.15];
In[17]:=
%["Best"]
Out[17]=
In[18]:=
%%["NumberOfFunctionCalls"]
Out[18]=

ProbabilitySelectBestPath (1) 

Have 50% of the individuals select the current shortest path (low exploration):

In[19]:=
SeedRandom[314];
In[20]:=
ResourceFunction["AntColonyOptimization"][
  Table[RandomSample[Range[10]], {i, Range[10]}], Total[#] &, All,
  "ProbabilitySelectBestPath" -> 0.5];
In[21]:=
%["Best"]
Out[21]=
In[22]:=
%%["NumberOfFunctionCalls"]
Out[22]=

PheromoneFunction (1) 

Use a lower pheromone marking at each step depending on individual fitness value (low convergence speed, high exploration):

In[23]:=
SeedRandom[314];
In[24]:=
ResourceFunction["AntColonyOptimization"][
  Table[RandomSample[Range[10]], {i, Range[10]}], Total[#] &, All,
  "PheromoneFunction" -> (1/(20 #) &)];
In[25]:=
%["Best"]
Out[25]=
In[26]:=
%%["NumberOfFunctionCalls"]
Out[26]=

WeatheringFactor (1) 

Reduce the selection probability for less-used paths more slowly:

In[27]:=
SeedRandom[314];
In[28]:=
ResourceFunction["AntColonyOptimization"][
  Table[RandomSample[Range[10]], {i, Range[10]}], Total[#] &, All,
  "WeatheringFactor" -> 0.001];
In[29]:=
%["Best"]
Out[29]=
In[30]:=
%%["NumberOfFunctionCalls"]
Out[30]=

Applications (5) 

How much money is in your pocket (knapsack problem)? (5) 

Total given mass:

In[31]:=
mWallet = Quantity[92., "Grams"];

Values for individual coins:

In[32]:=
moneyValue = <|"1¢" -> Quantity[0.01, "Euros"], "2¢" -> Quantity[0.02, "Euros"], "5¢" -> Quantity[0.05, "Euros"], "10¢" -> Quantity[0.1, "Euros"], "20¢" -> Quantity[0.2, "Euros"], "50¢" -> Quantity[0.5, "Euros"], "1€" -> Quantity[1., "Euros"], "2€" -> Quantity[2., "Euros"]|>;
In[33]:=
moneyMass = <|"1¢m" -> Quantity[2.3, "Grams"], "2¢m" -> Quantity[3.06, "Grams"], "5¢m" -> Quantity[3.92, "Grams"],
    "10¢m" -> Quantity[4.1, "Grams"], "20¢m" -> Quantity[5.74, "Grams"], "50¢m" -> Quantity[7.8, "Grams"], "1€m" -> Quantity[7.5, "Grams"], "2€m" -> Quantity[8.5, "Grams"]|>;

Fitness function (error function). In this problem the difference between the target coin weight and the calculated coin weight is minimized:

In[34]:=
massError[coins_List] := Module[{eps = 10.^-6, error},
  error = QuantityMagnitude[Abs@((Values[moneyMass] . coins) - mWallet), "Grams"];
  If[(Total@coins) == 0,
   10000.,
   If[error <= eps, eps, error]]]

Matrix of possible coins in wallet:

In[35]:=
possibCoins = Values@<|"1¢Range" -> Range[0, 5], "2¢Range" -> Range[0, 5], "5¢Range" -> Range[0, 5], "10¢Range" -> Range[0, 8], "20¢Range" -> Range[0, 7], "50¢Range" -> Range[0, 7], "1€Range" -> Range[0, 7], "2€Range" -> Range[0, 7]|>;
possibCoins // Grid
Out[8]=

Three optimization runs:

In[36]:=
Table[Values[
    moneyMass] . (result = (ResourceFunction["AntColonyOptimization"][
       possibCoins, massError, "PopulationSize" -> 1000, "FinishAtPheromoneThreshold" -> 0.99, "WeatheringFactor" -> 0.01, "MaxIterations" -> 30])) -> Values[moneyValue] . result,
 {3}]
Out[36]=

Properties and Relations (1) 

The fitness function must return a positive value. For maximization problems, fitness must be adjusted:

In[37]:=
ResourceFunction["AntColonyOptimization"][
 {
  {4, 2, 3, 1, 5, 6},
  {3, 1, 6, 4},
  {6, 3, 4, 1, 2, 5},
  {6},
  {5, 6},
  {2, 6, 1}
  },
 (Total[1/#] &)]
Out[37]=

Neat Examples (7) 

Minimize the mass of a rod structure (7) 

Optimize a structure geometry:

In[38]:=
nodes = {{1, {0, 0}}, {2, {0, 20}}, {3, {20, 0}}, {4, {20, 20}}, {5, {
    40, 0}}, {6, {40, 10}}, {7, {60, 0}}, {8, {60, 10}}, {9, {80, 0}}};
In[39]:=
rodElements = {{1, {1, 2, 
Subscript[a, 1]}}, {2, {1, 3, 
Subscript[a, 2]}}, {3, {2, 4, 
Subscript[a, 3]}}, {4, {4, 3, 
Subscript[a, 4]}}, {5, {1, 4, 
Subscript[a, 5]}}, {6, {2, 3, 
Subscript[a, 6]}}, {7, {3, 5, 
Subscript[a, 7]}}, {8, {5, 6, 
Subscript[a, 8]}}, {9, {4, 6, 
Subscript[a, 9]}}, {10, {3, 6, 
Subscript[a, 10]}}, {11, {4, 5, 
Subscript[a, 11]}}, {12, {5, 7, 
Subscript[a, 12]}}, {13, {7, 8, 
Subscript[a, 13]}}, {14, {6, 8, 
Subscript[a, 14]}}, {15, {5, 8, 
Subscript[a, 15]}}, {16, {6, 7, 
Subscript[a, 16]}}, {17, {7, 9, 
Subscript[a, 17]}}, {18, {8, 9, 
Subscript[a, 18]}}};
In[40]:=
structUndeformed = Graphics[
Sequence[{
ReplaceAll[rodElements, {
Pattern[no, 
Blank[Integer]], 
Pattern[prop, 
Blank[List]]} :> {Black, Dashed, 
AbsoluteThickness[
ReplaceAll[
Part[prop, 3], 
Table[Subscript[a, i] -> 1, {i, 18}]]], 
Line[{
Part[nodes, 
Part[prop, 1], 2], 
Part[nodes, 
Part[prop, 2], 2]}]}], 
ReplaceAll[nodes, {
Pattern[no, 
Blank[Integer]], 
Pattern[coord, 
Blank[List]]} -> {Gray, 
PointSize[0.02], 
Point[coord]}]}, AspectRatio -> Automatic]];
In[41]:=
Show[
Sequence[structUndeformed, 
Graphics[
Arrow[{{80, 10}, {80, 0}}]], AspectRatio -> Automatic]]
Out[41]=
In[42]:=
force = 1000;

Restrictions:

In[43]:=
rodSectionsValid = Table[10 i, {i, 1, 4}];
In[44]:=
maxDeform = 12.0; (* at the point of force application in force direction *)
In[45]:=
((Length@
  rodSectionsValid)^Length@rodElements) (* number of configurations, search space *)
Out[45]=

Define helper functions:

In[46]:=
addSubMatrix[orgM_?MatrixQ, subM_?MatrixQ, atX_Integer, atY_Integer] := Module[{M = orgM, M1 = subM, rows = Length[subM], columns = Length[
Part[subM, 1]], i, j}, For[i = 1, i <= rows, For[j = 1, j <= columns, Part[M, atX - 1 + i, atY - 1 + j] = Part[
         M, atX - 1 + i, atY - 1 + j] + Part[M1, i, j]; Increment[
       j]]; Increment[i]]; M]
In[47]:=
generateStiffnesMatrix[
  section_List] := Module[{stiff, rodelements, rodProp, rodLength, kSub, allStiff, c, s}, stiff = 10000.; rodelements = {{1, {1, 2, stiff, 
Part[section, 1]}}, {2, {1, 3, stiff, 
Part[section, 2]}}, {3, {2, 4, stiff, 
Part[section, 3]}}, {4, {4, 3, stiff, 
Part[section, 4]}}, {5, {1, 4, stiff, 
Part[section, 5]}}, {6, {2, 3, stiff, 
Part[section, 6]}}, {7, {3, 5, stiff, 
Part[section, 7]}}, {8, {5, 6, stiff, 
Part[section, 8]}}, {9, {4, 6, stiff, 
Part[section, 9]}}, {10, {3, 6, stiff, 
Part[section, 10]}}, {11, {4, 5, stiff, 
Part[section, 11]}}, {12, {5, 7, stiff, 
Part[section, 12]}}, {13, {7, 8, stiff, 
Part[section, 13]}}, {14, {6, 8, stiff, 
Part[section, 14]}}, {15, {5, 8, stiff, 
Part[section, 15]}}, {16, {6, 7, stiff, 
Part[section, 16]}}, {17, {7, 9, stiff, 
Part[section, 17]}}, {18, {8, 9, stiff, 
Part[section, 18]}}}; rodProp = ReplaceAll[rodelements, {
Pattern[no, 
Blank[Integer]], 
Pattern[prop, 
Blank[List]]} :> {no, {c = (Part[nodes, 
Part[prop, 2], 2, 1] - Part[nodes, 
Part[prop, 1], 2, 1])/(rodLength = Sqrt[(Part[nodes, 
Part[prop, 2], 2, 1] - Part[nodes, 
Part[prop, 1], 2, 1])^2 + (Part[nodes, 
Part[prop, 2], 2, 2] - Part[nodes, 
Part[prop, 1], 2, 2])^2]), s = (Part[nodes, 
Part[prop, 2], 2, 2] - Part[nodes, 
Part[prop, 1], 2, 2])/rodLength, Part[prop, 3] (Part[prop, 4]/rodLength)}}]; kSub = ReplaceAll[
    rodProp, {
Pattern[no, 
Blank[Integer]], 
Pattern[prop, 
Blank[List]]} :> {no, Part[prop, 3]
         ReplaceAll[{{c1^2, c1 s1}, {c1 s1, s1^2}}, {c1 -> Part[
           prop, 1], s1 -> Part[prop, 2]}]}]; allStiff = Table[
Table[0, {i, 1, 2 Length[nodes]}], {j, 1, 2 Length[nodes]}]; For[
   i = 1, i <= Length[rodelements], allStiff = addSubMatrix[allStiff, 
Part[kSub, i, 2], Part[rodelements, i, 2, 1] 2 - 1, Part[rodelements, i, 2, 1] 2 - 1]; allStiff = addSubMatrix[
      allStiff, -Part[kSub, i, 2], Part[rodelements, i, 2, 1] 2 - 1, Part[rodelements, i, 2, 2] 2 - 1]; allStiff = addSubMatrix[
      allStiff, -Part[kSub, i, 2], Part[rodelements, i, 2, 2] 2 - 1, Part[rodelements, i, 2, 1] 2 - 1]; allStiff = addSubMatrix[
      allStiff, 
Part[kSub, i, 2], Part[rodelements, i, 2, 2] 2 - 1, Part[rodelements, i, 2, 2] 2 - 1]; Increment[i]]; allStiff]
In[48]:=
computeMass[
  section_List] := Module[{length = {20, 20, 20, 20, 20 Sqrt[2], 20 Sqrt[2], 20, 10, 10 Sqrt[5], 10 Sqrt[5], 20 Sqrt[2], 20, 10, 20, 10 Sqrt[5], 10 Sqrt[5], 20, 10 Sqrt[5]}}, 
Dot[length, section]]
In[49]:=
computeDeformation[
  section_List] := Module[{sol, allStiff, uAll, unKnown, allForces, freac1, freac3, freac4, u2, u5, u6, u7, u8, u9, u10, u11, u12, u13,
    u14, u15, u16, u17, u18}, allStiff = generateStiffnesMatrix[
    section]; uAll = {0, u2, 0, 0, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16, u17, u18}; unKnown = {freac1, freac3, freac4,
      u2, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16, u17, u18}; allForces = {freac1, 0, freac3, freac4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -force}; sol = Flatten[
Solve[Dot[allStiff, uAll] == allForces, unKnown]]]

The fitness function uses a penalty term to enforce constraints:

In[50]:=
objective[sections_] := With[{uMax = maxDeform, u18 = N@Abs@(computeDeformation[sections] // Last // Values), m = N@computeMass[sections]},
   If[u18 > uMax, m + 250 (u18 - uMax)^2, m]
   ];

Search for minimum mass structure according to the displacement constraints:

In[51]:=
result = ResourceFunction["AntColonyOptimization"][
   ConstantArray[rodSectionsValid, Length@rodElements],
   objective, "PopulationSize" -> 800, "MaxIterations" -> 20, "FinishAtPheromoneThreshold" -> 0.80];

Show the result:

In[52]:=
Thread[{
Sequence[
Subscript[a, 1], 
Subscript[a, 2], 
Subscript[a, 3], 
Subscript[a, 4], 
Subscript[a, 5], 
Subscript[a, 6], 
Subscript[a, 7], 
Subscript[a, 8], 
Subscript[a, 9], 
Subscript[a, 10], 
Subscript[a, 11], 
Subscript[a, 12], 
Subscript[a, 13], 
Subscript[a, 14], 
Subscript[a, 15], 
Subscript[a, 16], 
Subscript[a, 17], 
Subscript[a, 18]]} -> result] (* optimal sections *)
Out[52]=
In[53]:=
computeDeformation[result] // Values // Last (* restriktion *)
Out[53]=
In[54]:=
computeMass[result] // N
Out[54]=

Plot the result with deformed structure:

In[55]:=
solUAll = Partition[{0, Subscript[u, 2], 0, 0, Subscript[u, 5], Subscript[u, 6], Subscript[u, 7], Subscript[u, 8], Subscript[u, 9], Subscript[
     u, 10], Subscript[u, 11], Subscript[u, 12], Subscript[u, 13], Subscript[u, 14], Subscript[u, 15], Subscript[u, 16], Subscript[
     u, 17], Subscript[u, 18]} /. Thread[{Subscript[u, 2], Subscript[u, 5], Subscript[u, 6], Subscript[u, 7], Subscript[u, 8], Subscript[u, 9], Subscript[u,
        10], Subscript[u, 11], Subscript[u, 12], Subscript[u, 13], Subscript[u, 14], Subscript[u, 15], Subscript[u, 16], Subscript[u, 17], Subscript[u, 18]} -> Values[computeDeformation[result]][[4 ;; -1]]], 2];
In[56]:=
defNodes = Table[{i, nodes[[i, 2]] + solUAll[[i]]}, {i, Length[nodes]}];
In[57]:=
structDeformed = Graphics[{rodElements /. {no_Integer, prop_List} :> {Blue, AbsoluteThickness[
        prop[[3]] /. Table[Subscript[a, i] -> result[[i]]/10, {i, 18}]], Line[{defNodes[[prop[[1]], 2]], defNodes[[prop[[2]], 2]]}]}, defNodes /. {no_Integer, coord_List} -> {Hue[0], PointSize[0.02`],
        Point[coord]}}, AspectRatio -> Automatic];
In[58]:=
Show[structUndeformed, structDeformed, Graphics[Arrow[{{80, 10}, {80, 0}}]], AspectRatio -> Automatic]
Out[58]=

Publisher

Lars Ulke-Winter

Version History

  • 1.0.0 – 13 February 2020

Source Metadata

Related Resources

Author Notes

Detailed explanation what is being optimized in the first example:

The array (called possibValues in the Basic Examples section) shows the valid variable assignments in each row. Each row stands for one variable with its corresponding possible values. In this example you have 10 variables, each of which can have 10 different states. This results in 1010 possible variable assignments (the size of the search space in this example).

The optimization task is to find the shortest path (the minimum) through array (with its 1010 different permissible configurations) with regard to the defined function fit. For the sake of simplicity, the total value of the individual variables was simply chosen as the objective function in this example, fit=Total[#]&. Thus, for example, a path as marked in red here has a corresponding path length of Total[{3,9,1,6,10,9,7,3,6,8}]62:

The optimal (shortest possible) path to be found by the AntColonyOptimization algorithm after some iterations is the following, which achieves Total[{1,1,1,1,1,1,1,1,1,1}]10 as the absolute minimum of the search space:

This example is rather trivial and only serves to outline the basic procedure and introduce the basic parameters, (such as array and fitness) that are necessary to formulate the optimization problem. Slightly more complex application examples are demonstrated in the various Example subsections above.

A brief summary:

array (a list of lists) defines the search space, where the number of rows corresponds to the number of variables and the number of columns corresponds to the possible individual states of each variable. The number of possible variable states can generally be different for each variable.
The function fitness assigns each valid variable configuration an associated quantity, the corresponding path length. This allows comparison of the different variable configurations.
Within the given search space (encoded in array) and the given virtual path length (defined by fitness), AntColonyOptimization tries to find the shortest path (the minimum) through the search space. It returns the related variable assignments.

License Information