Wolfram Research

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.
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" 500 number of ants
"MaxIterations" 200 maximum iterations
"FinishAtPheromoneThreshold" 0.95 stop criterion
"PheromoneFunction" (1/#) & pheromone concentration to be deposited depending on fitness
"WeatheringFactor" 0.05 fraction of decreasing pheromone concentration (selection probability) per iteration
"ProbabilitySelectBestPath" 0.02 proportion 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

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

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

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

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

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

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

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

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

Pay attention, fitness function must return a positive value. For maximization problems, the 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

Minimize the mass of a rod structure

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]=

Resource History

Source Metadata

Related Resources

Author Notes

Detailed explanation what is being optimized in this example: The array possibValues in this simple example 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 them can have 10 different states. This results in 1010 possible variable assignments (size of the search space in this example). The optimization task is to find the shortest path (minimum) through this 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 follows: with Total[{1,1,1,1,1,1,1,1,1,1}]10 (absolute minimum of the search space). This example is rather trivial and has only the purpose to outline the basic procedure and to introduce the basic parameters, such as array and fitness function, which are necessary to formulate the optimization problem. Slightly more complex application examples follow below.

A brief summary : -The array (list of lists) defines the search space (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 quality, the corresponding path length. This allows to compare the different variable configurations. - Within the given search space (encoded in array) and the given virtual path length defined by fitness, AntColonyOptimization try to find the shortest path (minimum) through the search space and return its related variable assignments .

License Information