Русский
preview
Butterfly Optimization Algorithm (BOA)

Butterfly Optimization Algorithm (BOA)

MetaTrader 5Trading |
195 0
Andrey Dik
Andrey Dik

Table of Contents

  1. Introduction
  2. Algorithm Implementation
  3. Test Results
  4. Conclusions


Introduction

In this article, we will examine one of the newer optimization algorithms inspired by nature. Butterflies are among the most amazing creatures, and their lives are inextricably linked to searching: for food, for a mate, and for a place to lay their eggs. Unlike many other insects, butterflies have a unique system of chemical communication based on the detection and release of aromatic substances — pheromones. This very feature formed the basis of the Butterfly Optimization Algorithm (BOA), proposed by Indian researchers S. Arora and S. Singh in 2019.

In nature, butterflies use specialized chemoreceptors located on their antennae to detect chemical signals over astonishingly long distances. Males of some species can detect female pheromones from several kilometers away. The intensity of a perceived fragrance depends on two factors: the strength of the fragrance source and the distance to it. According to Stevens' law, which describes the psychophysical perception of stimuli in living organisms, the perceived magnitude of a sensation is related to the physical intensity of the stimulus by a power law. The authors of the algorithm formalized this principle as the fragrance equation: f = c·I^a, where f is the perceived fragrance magnitude, "I" is the stimulus intensity (related to the quality of the food source), "c" is the sensory modality (the butterfly’s ability to perceive odors), and the exponent "a" determines the nature of the relationship between perception and intensity.

The foraging behavior of butterflies can be divided into two phases. When a butterfly detects the strong fragrance of a flower rich in nectar, it flies directly toward the source — this is global search, guided by the best-known resource. However, if there is no clear leader or the butterfly is in a cloud of mixed odors from several flowers, it makes local movements between the nearest sources, exploring the surrounding area — this is local search. The probability of choosing between global and local search is determined by the switching parameter p.

An important feature of chemical communication in butterflies is that the fragrance fades with distance. Pheromone molecules disperse in the air, are absorbed by obstacles, and are broken down by ultraviolet light. The parameter "a" in the fragrance formula models precisely this absorption effect: when "a" is close to zero, the fragrance spreads with virtually no loss and can be detected from any point in space, which facilitates global exploration; when "a" is close to one, the fragrance rapidly fades, and the butterflies orient themselves primarily toward the nearest sources, which reinforces the local exploitation of found solutions.

Another characteristic feature of butterfly behavior is that they themselves are sources of fragrance. Each butterfly releases pheromones that attract other butterflies. The better the spot it finds (more nectar, better conditions), the stronger its pheromone signal. In the algorithm, this is modeled by the relationship between the stimulus intensity "I" and the value of the objective function: butterflies with higher fitness emit a stronger fragrance and have a greater influence on the population's movement.

Thus, the swarm of virtual butterflies in the BOA algorithm constitutes a self-organizing system in which each individual serves simultaneously as both a seeker and a beacon for the others.



Algorithm Implementation

While implementing the BOA butterfly optimization algorithm, I ran into a serious problem: the original formulas from the research paper caused the algorithm to behave incorrectly. Instead of converging toward the optimum, the butterfly population systematically collapsed toward the origin, and on multidimensional functions, the algorithm performed worse than a random search. A detailed mathematical analysis revealed the cause of this phenomenon.

In the original paper, the global search formula is written as x_new = x + (r²·g* - x)·f, where x is the butterfly's current position, g* is the best found solution, r is a random number in [0,1], and f is the butterfly’s fragrance. At first glance, the formula seems logical: the butterfly should move toward the best solution. However, on closer inspection of the expression in parentheses, (r²·g* - x), the problem becomes obvious.

Let's consider a specific example: suppose the best solution is g* = 10, the current position is x = 5, the random number is r = 0.7 (so r² = 0.49), and the fragrance is f = 0.5. Substituting the values, we get: x_new = 5 + (0.49·10 - 5)·0.5 = 5 + (4.9 - 5)·0.5 = 5 + (-0.1)·0.5 = 4.95. The butterfly moved not toward the optimum g* = 10, but in the opposite direction — toward zero! This happens because, when r < 1, the expression r²·g* always yields a value less than g*, and the difference (r²·g* - x) points not toward g*, but toward the point r²·g*, which is closer to zero.

A similar problem exists in the local search formula x_new = x + (r²·x_j - x_k)·f. The expression (r²·x_j - x_k) does not create a meaningful direction between the two butterflies j and k. Instead, it generates a vector that is systematically shifted toward the origin because of the factor r² < 1 multiplying x_j.

It is interesting to note that the original formulas can yield acceptable results on standard test functions (Sphere, Rastrigin, Ackley, and others), whose global optimum lies precisely at the point (0, 0, ..., 0). In this case, the shift toward the origin coincides by chance with the movement toward the optimum, which masks an error in the algorithm's logic. However, on real optimization problems, where the optimum may be located at an arbitrary point in the search space, the algorithm proves to be completely ineffective.

To remedy the situation, it was necessary to change the placement of parentheses in the formulas while retaining all the original components of the algorithm. The adjusted global search formula is x_new = x + r²·(g* - x)·f. Now the expression (g* - x) represents a vector pointing from the current position toward the best solution, which corresponds to the biological metaphor: a butterfly flies toward the fragrance of the best flower. The multiplier r² scales the step length, and the fragrance f determines the intensity of the motion. Let's check using the same example: x_new = 5 + 0.49·(10 - 5)·0.5 = 5 + 0.49·5·0.5 = 5 + 1.225 = 6.225. Now the butterfly moves correctly from x = 5 toward g* = 10.

The local search formula was adjusted in the same way: x_new = x + r²·(x_j - x_k)·f. The expression (x_j - x_k) now represents the difference vector between two random butterflies, setting the direction of local exploration. This results in random wandering within the current population distribution, which corresponds to the butterflies' local search for food within a limited area.

Mathematical notation used in scientific articles may contain ambiguities or errors, and simply copying equations without understanding their geometric meaning can lead to an algorithm that does not work. When implementing any optimization method, it is necessary to analyze the physical or biological meaning of each operation and verify that the algorithm operates correctly on problems with a known optimum located at different points in the search space.

BOA_Illustration

Illustration of the BOA Algorithm in Action

The illustration shows a clearing at night with flowers and butterflies. The best flower (g*) is a large yellow one with a golden halo of fragrance and concentric waves of fragrance propagation. Three butterflies (orange, blue, and purple) are flying toward it along dotted orange paths — global search. The other two (yellow and turquoise) explore the area around the second flower in a zigzag pattern along green trajectories — a local search.

Let's write pseudocode for the BOA algorithm.

INPUT:

N — population size
c — sensory modality, [0, 1]
aStart — initial absorption exponent
p — probability of switching between global and local search
T — total number of epochs
D — dimension of the search space
min[d], max[d] — search bounds for each coordinate d

OUTPUT:
g* — best found position
fB — best objective function value

//-----------------------------------------------------------------------------
// Initialization
//-----------------------------------------------------------------------------
1. aCurrent ← aStart
2. fB ← −∞
3. FOR EACH butterfly i = 1..N:
4. fragrance[i] ← c
5. intensity[i] ← 0.5
6. FOR EACH coordinate d = 1..D:
7. x[i][d] ← min[d] + rand(0,1) × (max[d] − min[d])
8. f[i] ← FitnessFunction(x[i])


//-----------------------------------------------------------------------------
// Main loop
//-----------------------------------------------------------------------------
9. FOR EACH epoch t = 1..T:

//--- Update the global best solution ---
10. FOR EACH butterfly i = 1..N:
11. IF f[i] &gt; fB:
12. fB ← f[i]
13. g* ← x[i]

//--- Stimulus intensity normalization ---
14. fMin ← min(f[1..N])
15. fMax ← max(f[1..N])
16. range ← fMax − fMin
17. FOR EACH butterfly i = 1..N:
18. IF range &lt; ε:
19. intensity[i] ← 0.5
20. ELSE:
21. intensity[i] ← 0.1 + 0.9 × (f[i] − fMin) / range

//--- Update the absorption exponent ---
22. aCurrent ← aStart + (t / T) × (1.0 − aStart)
23. IF aCurrent &gt; 1.0: aCurrent ← 1.0

//--- Calculate the fragrance (Formula 1) ---
24. FOR EACH butterfly i = 1..N:
25. fragrance[i] ← c × intensity[i] ^ aCurrent

//--- Moving butterflies ---
26. FOR EACH butterfly i = 1..N:
27. rnd ← rand(0,1)

28. IF rnd < p:
//--- Global search (Formula 2) ---
29. FOR EACH coordinate d = 1..D:
30. r ← rand(0,1)
31. r² ← r × r
32. direction ← g*[d] − x[i][d]
33. x[i][d] ← x[i][d] + r² × direction × fragrance[i]

34. ELSE:
//--- Local search (Formula 3) ---
35. j ← a random butterfly from [1..N]
36. k ← a random butterfly from [1..N], k ≠ j
37. FOR EACH coordinate d = 1..D:
38. r ← rand(0,1)
39. r² ← r × r
40. direction ← x[j][d] − x[k][d]
41. x[i][d] ← x[i][d] + r² × direction × fragrance[i]

//--- Boundary check ---
42. FOR EACH coordinate d = 1..D:
43. IF x[i][d] < min[d]: x[i][d] ← min[d]
44. IF x[i][d] > max[d]: x[i][d] ← max[d]

//--- Calculating fitness ---
45. f[i] ← FitnessFunction(x[i])

46. RETURN g*, fB

Let's move on to implementing the BOA algorithm.

The C_AO_BOA class inherits from the base class C_AO and implements the Butterfly Optimization Algorithm. The constructor specifies the algorithm's name, its description, and a link to the article, as well as the default parameter values: population size, sensory modality, initial absorption exponent, and the probability of switching between global and local search.

The "params" array contains four elements, each of which stores the name and value of the corresponding parameter. The public fields "sensorModC," "aStart," and "switchP" are accessible from outside the class for directly setting their values before initialization. The private section of the class contains the current absorption exponent "aCurrent," the total number of epochs "epochs," the current epoch counter "epochNow," as well as two arrays: "fragrance" for storing the calculated fragrance of each butterfly and "intensity" for storing the normalized stimulus intensity. The SetParams method reads values from the "params" array and assigns them to the corresponding class fields.

//————————————————————————————————————————————————————————————————————
class C_AO_BOA : public C_AO
{
  public:
  ~C_AO_BOA () { }

  C_AO_BOA ()
  {
    ao_name = "BOA";
    ao_desc = "Butterfly Optimization Algorithm";
    ao_link = "https://www.mql5.com/en/articles/21209";

    popSize    = 50;
    sensorModC = 0.9;
    aStart     = 0.5;
    switchP    = 0.8;

    ArrayResize (params, 4);
    params [0].name = "popSize";    params [0].val = popSize;
    params [1].name = "sensorModC"; params [1].val = sensorModC; // sensory modality [0, 1]
    params [2].name = "aStart";     params [2].val = aStart;     // initial absorption exponent
    params [3].name = "switchP";    params [3].val = switchP;    // switch probability
  }

  void SetParams ()
  {
    popSize    = (int)params [0].val;
    sensorModC = params      [1].val;
    aStart     = params      [2].val;
    switchP    = params      [3].val;
  }

  bool Init (const double &rangeMinP  [],
             const double &rangeMaxP  [],
             const double &rangeStepP [],
             const int     epochsP = 0);

  void Moving   ();
  void Revision ();

  //------------------------------------------------------------------
  double sensorModC;  // sensory modality c (sensory modality)
  double aStart;      // initial absorption exponent
  double switchP;     // switch probability p (switch probability)

  private: //—————————————————————————————————————————————————————————
  double aCurrent;    // current absorption exponent
  int    epochs;      // total number of epochs
  int    epochNow;    // the current epoch

  double fragrance []; // the fragrance of each butterfly
  double intensity []; // stimulus intensity (normalized fitness)

  void CalculateFragrance ();
  void NormalizeIntensity ();
};
//————————————————————————————————————————————————————————————————————

The Init method takes arrays of the minimum and maximum bounds of the search space, an array of discretization steps, and the total number of epochs. First, the StandardInit method of the base class is called, which performs standard initialization: allocating memory for the population, copying the search bounds and steps, and determining the number of coordinates. If standard initialization fails, the method returns 'false'.

Next, the total number of epochs is stored, the current epoch counter is reset to zero, and the current absorption exponent is set to the initial value "aStart". Next, memory is allocated for the fragrance and intensity arrays of size "popSize," and each element of the fragrance array is assigned the initial value "sensorModC", while each element of the intensity array is assigned a value of 0.5. This means that, before the first fitness evaluation, all butterflies are considered equally attractive, and their initial fragrance is determined solely by the sensory modality.

//————————————————————————————————————————————————————————————————————
bool C_AO_BOA::Init (const double &rangeMinP  [],
                     const double &rangeMaxP  [],
                     const double &rangeStepP [],
                     const int     epochsP = 0)
{
  if (!StandardInit (rangeMinP, rangeMaxP, rangeStepP)) return false;

  //------------------------------------------------------------------
  epochs   = epochsP;
  epochNow = 0;
  aCurrent = aStart;

  //------------------------------------------------------------------
  // Initialize the fragrance and intensity arrays
  //------------------------------------------------------------------
  ArrayResize (fragrance, popSize);
  ArrayResize (intensity, popSize);

  for (int i = 0; i < popSize; i++)
  {
    fragrance [i] = sensorModC;
    intensity [i] = 0.5;
  }

  return true;
}
//————————————————————————————————————————————————————————————————————

The Moving method implements the movement of all butterflies in the current iteration. At the beginning, the "epochNow" counter is incremented. If the "revision" flag is set to 'false' — which indicates the first iteration — each butterfly is assigned a random position: for each coordinate, a random number is generated from a uniform distribution on the interval from zero to one, and the position is calculated as a linear interpolation between the minimum and maximum bounds, followed by discretization using the SeInDiSp function. After initialization, the "revision" flag is set to 'true', and the method completes.

In all subsequent iterations, the CalculateFragrance method is called first; it calculates the fragrance of each butterfly based on the current stimulus intensity and absorption exponent. Next, for each butterfly "i", a random number is generated and compared with the switching probability "switchP". If the random number is less than "switchP", a global search is performed: for each coordinate "d," a random number "r" is generated, its square "r²" is calculated, the direction is determined as the difference between the coordinate of the global best solution "cB" and the butterfly’s current coordinate, and then the position is updated by adding the product of "r²", the direction, and the fragrance of that butterfly. Thus, the butterfly moves toward the best found solution, with the step size proportional to the distance to the target, a random multiplier, and the fragrance strength.

If the random number is greater than or equal to "switchP", a local search is performed. Two different butterflies, "j" and "k", are randomly selected from the population. For each coordinate "d", "r²" is generated in the same way, the direction is calculated as the difference between the coordinates of butterflies "j" and "k", and the position is updated by adding the product of "r²", the direction, and the fragrance. This mechanism creates a random walk whose scale is determined by the current spread of the population.

Additionally, after a local search, with a probability of 0.2, one randomly selected coordinate of the butterfly is replaced with a value generated from a normal distribution centered at the corresponding coordinate of the best solution "cB" and bounded by the limits of the search space. This technique is not part of the original algorithm and serves to prevent the population from getting stuck in local extrema by introducing controlled noise near the best-known solution.

Finally, for each coordinate of each butterfly, the SeInDiSp function is called, which both limits the value to the permissible range and performs discretization with the specified step size.

//————————————————————————————————————————————————————————————————————
void C_AO_BOA::Moving ()
{
  epochNow++;

  //------------------------------------------------------------------
  // First iteration—random population initialization
  // The positions of the butterflies are generated randomly within the search space
  //------------------------------------------------------------------
  if (!revision)
  {
    for (int i = 0; i < popSize; i++)
    {
      for (int d = 0; d < coords; d++)
      {
        double rnd = u.RNDfromCI (0.0, 1.0);
        a [i].c [d] = rangeMin [d] + rnd * (rangeMax [d] - rangeMin [d]);
        a [i].c [d] = u.SeInDiSp (a [i].c [d], rangeMin [d], rangeMax [d], rangeStep [d]);
      }
    }
    revision = true;
    return;
  }

  //------------------------------------------------------------------
  // Calculate the fragrance for each butterfly (formula 1)
  // f = c * I^a
  //------------------------------------------------------------------
  CalculateFragrance ();

  //------------------------------------------------------------------
  // Main cycle—movement of butterflies
  //------------------------------------------------------------------
  for (int i = 0; i < popSize; i++)
  {
    if (u.RNDprobab () < switchP)
    {
      //==============================================================
      // GLOBAL SEARCH (formula 2, revised)
      // x_i(t+1) = x_i(t) + r² × (g - x_i(t)) × f_i
      // The butterfly moves toward the global best solution g
      //
      // Logic: the direction (g - x) points from the current position toward the best one,
      // the step is scaled by r² (randomness) and f (fragrance/fitness)
      //==============================================================
      for (int d = 0; d < coords; d++)
      {
        double r = u.RNDprobab ();
        double r2 = r * r;

        // Direction from the current position to the best position
        double direction = cB [d] - a [i].c [d];

        // Step: r² × direction × fragrance
        double step = r2 * direction * fragrance [i];

        a [i].c [d] = a [i].c [d] + step;
      }
    }
    else
    {
      //==============================================================
      // LOCAL SEARCH (formula 3, adjusted)
      // x_i(t+1) = x_i(t) + r² × (x_j(t) - x_k(t)) × f_i
      // Random walk: the direction is determined by the difference between
      // two random butterflies, j and k
      //
      // Logic: (x_j - x_k) yields a random direction in space,
      // determined by the current distribution of the population
      //==============================================================

      // Select two random butterflies, j and k
      int j = u.RNDminusOne (popSize);
      int k = u.RNDminusOne (popSize);

      // Ensure that j and k are different
      while (k == j) k = u.RNDminusOne (popSize);

      for (int d = 0; d < coords; d++)
      {
        double r = u.RNDprobab ();
        double r2 = r * r;

        // A random direction between two butterflies
        double direction = a [j].c [d] - a [k].c [d];

        // Step: r² × direction × fragrance
        double step = r2 * direction * fragrance [i];

        a [i].c [d] = a [i].c [d] + step;
      }

      //==============================================================
      // Adding noise to the coordinates of the best solution
      // to prevent getting stuck in local extrema
      // This is not present in the original BOA
      //==============================================================
      if (u.RNDprobab () < 0.2)
      {
        int ind = u.RNDminusOne (coords);
        a [i].c [ind] = u.GaussDistribution (cB [ind], rangeMin [ind], rangeMax [ind], 1);
      }
    }

    //----------------------------------------------------------------
    // Boundary checking and discretization
    //----------------------------------------------------------------
    for (int d = 0; d < coords; d++)
    {
      a [i].c [d] = u.SeInDiSp (a [i].c [d], rangeMin [d], rangeMax [d], rangeStep [d]);
    }
  }
}
//————————————————————————————————————————————————————————————————————

The Revision method performs three operations after all butterflies have been moved and their fitness has been calculated. The first operation is to update the global best solution: for each butterfly "i", the system checks whether its fitness value exceeds the current best value "fB", and if so, "fB" is updated and the butterfly's position is copied into the array of best coordinates "cB". The second operation is to call the "NormalizeIntensity" method, which recalculates the stimulus intensity for each butterfly based on the current fitness values.

The third operation is to update the absorption exponent "aCurrent" by linearly interpolating from the initial value "aStart" to one, based on the ratio of the current epoch to the total number of epochs. If the calculated value exceeds one, it is capped at one. The increase in the absorption exponent from small values toward one means that, at the start of optimization, the fragrance of all butterflies is approximately the same and the search proceeds uniformly, whereas by the end, butterflies with high fitness receive a significantly stronger fragrance, which enhances the exploitation of the promising regions that have been found.

//————————————————————————————————————————————————————————————————————
void C_AO_BOA::Revision ()
{
  //------------------------------------------------------------------
  // 1. Update the global best solution
  //------------------------------------------------------------------
  for (int i = 0; i < popSize; i++)
  {
    if (a [i].f > fB)
    {
      fB = a [i].f;
      ArrayCopy (cB, a [i].c, 0, 0, coords);
    }
  }

  //------------------------------------------------------------------
  // 2. Normalize the stimulus intensity (I)
  //    Intensity I is proportional to fitness f(x)
  //------------------------------------------------------------------
  NormalizeIntensity ();

  //------------------------------------------------------------------
  // 3. Update the absorption exponent a
  //    According to the pseudocode: "Update the value of a"
  //    a increases linearly from aStart to 1.0
  //    This gradually reduces the influence of intensity on the fragrance
  //------------------------------------------------------------------
  if (epochs > 0)
  {
    aCurrent = aStart + (double)epochNow / (double)epochs * (1.0 - aStart);
    if (aCurrent > 1.0) aCurrent = 1.0;
  }
}
//————————————————————————————————————————————————————————————————————

The NormalizeIntensity method is responsible for normalizing fitness values to a uniform stimulus intensity scale. First, the minimum and maximum fitness values are determined among all butterflies in the current population. If the difference between the maximum and minimum fitness values is negligible — less than 10⁻¹⁰, which means that all fitness values are practically equal — all butterflies are assigned the same intensity of 0.5.

Otherwise, for each butterfly, the intensity is calculated using the formula for a linear mapping of the fitness value to the range from 0.1 to 1.0: I = 0.1 + 0.9 × (f − fMin) / (fMax − fMin). The butterfly with the worst fitness receives an intensity of 0.1, and the one with the best fitness receives an intensity of 1.0. The lower bound of 0.1 was intentionally chosen to be nonzero: even the worst butterflies must have a nonzero fragrance in order to retain the ability to participate in the search process and avoid becoming completely passive agents.

//————————————————————————————————————————————————————————————————————
// Normalization of stimulus intensity
// According to the article: "Stimulus Intensity I at x is determined by f(x)"
// I is proportional to fitness
//————————————————————————————————————————————————————————————————————
void C_AO_BOA::NormalizeIntensity ()
{
  //------------------------------------------------------------------
  // Find the minimum and maximum fitness
  //------------------------------------------------------------------
  double minF = a [0].f;
  double maxF = a [0].f;

  for (int i = 1; i < popSize; i++)
  {
    if (a [i].f < minF) minF = a [i].f;
    if (a [i].f > maxF) maxF = a [i].f;
  }

  //------------------------------------------------------------------
  // Normalize the intensity to the range [0.1, 1.0]
  // The lower bound of 0.1 prevents the worst butterflies from having too little fragrance
  //------------------------------------------------------------------
  double range = maxF - minF;

  if (range < 1e-10)
  {
    // All fitness values are approximately equal
    for (int i = 0; i < popSize; i++)
    {
      intensity [i] = 0.5;
    }
  }
  else
  {
    for (int i = 0; i < popSize; i++)
    {
      intensity [i] = 0.1 + 0.9 * (a [i].f - minF) / range;
    }
  }
}
//————————————————————————————————————————————————————————————————————

The CalculateFragrance method calculates the fragrance of each butterfly using the algorithm's main formula: fragrance = c × I^a, where "c" is the sensory modality "sensorModC," "I" is the normalized stimulus intensity "intensity", and "a" is the current absorption exponent "aCurrent." This formula is based on Stevens' law from psychophysics, which describes the relationship between the physical intensity of a stimulus and its subjective perception.

When the value of the parameter "a" is small — as is typical at the start of optimization — raising the intensity to a low power yields values close to one for all butterflies, regardless of their fitness, and the fragrance of all butterflies turns out to be approximately equal to the sensory modality "c", which ensures uniform exploration of the space. As "a" approaches one, differences in intensity become increasingly apparent in the fragrance: butterflies with high fitness have a fragrance close to "c", while those with low fitness have a significantly weaker fragrance, which directs the search toward the most promising areas.

//————————————————————————————————————————————————————————————————————
// Calculate the fragrance using Formula 1: f = c * I^a
// where:
// c - sensory modality (determines the base step size)
// I - stimulus intensity (normalized fitness)
// a - absorption exponent
//
// For small a (beginning): I^a → 1 for all I, fragrance ≈ c (uniform search)
// For large a (end): I^a → I, the best butterflies have a stronger fragrance
//————————————————————————————————————————————————————————————————————
void C_AO_BOA::CalculateFragrance ()
{
  for (int i = 0; i < popSize; i++)
  {
    // Formula 1: fragrance = c * I^a
    fragrance [i] = sensorModC * MathPow (intensity [i], aCurrent);
  }
}
//————————————————————————————————————————————————————————————————————


Test Results

It is important to note that without the small block of code below, there is insufficient diversity in the butterfly population. This code snippet implements a stochastic mutation mechanism that is not present in the original description of the BOA algorithm and was added to compensate for its structural weakness. The problem is that, in the authors' original concept, the local search depends entirely on the difference between the positions of two random butterflies (x_j − x_k). As the population converges, the distances between butterflies decrease, the differences tend toward zero, and the steps of the local search become microscopic. The population loses the ability to move away from the vicinity of the current best solution, even if it is a local optimum.

The added mechanism works as follows: with a probability of 0.2, that is, for approximately every fifth butterfly during local search, one coordinate is randomly selected from all available ones, and its value is replaced with a number generated from a normal distribution centered at the corresponding coordinate of the global best solution "cB". The sigma parameter is equal to 1, and the values are limited to the valid coordinate range. This creates a controlled perturbation: the new value will, with high probability, be close to the best solution, but with some probability, it may be far enough away to push the butterfly into a new region of the search space.

It is important to note that only one coordinate is subject to mutation, not the entire position vector. This preserves most of the accumulated information about the promising region while introducing diversity along a single dimension. It is also important that the distribution is centered on the coordinate of the best found solution, rather than on the butterfly’s current position — thus, the mutation is not a blind random walk, but a directed exploration of the neighborhood of the best found point.

//==============================================================
      // Adding noise to the coordinates of the best solution
      // to prevent getting stuck in local extrema
      // this is not present in the original BOA
      //==============================================================
      if (u.RNDprobab () < 0.2)
      {
        int ind = u.RNDminusOne (coords);
        a [i].c [ind] = u.GaussDistribution (cB [ind], rangeMin [ind], rangeMax [ind], 1);
      }
    }

    //----------------------------------------------------------------

Here are the results without this perturbation.

BOA|Butterfly Optimization Algorithm|50.0|0.9|0.5|0.8|
=============================
5 Hilly's; Func runs: 10000; result: 0.6635729348082547
25 Hilly's; Func runs: 10000; result: 0.3933240668190092
500 Hilly's; Func runs: 10000; result: 0.2670951838030907
=============================
5 Forest's; Func runs: 10000; result: 0.5030842340559138
25 Forest's; Func runs: 10000; result: 0.31195194496509265
500 Forest's; Func runs: 10000; result: 0.1922137682213402
=============================
5 Megacity's; Func runs: 10000; result: 0.32461538461538464
25 Megacity's; Func runs: 10000; result: 0.16953846153846158
500 Megacity's; Func runs: 10000; result: 0.10381538461538553
=============================
Overall score: 2.92921 (32.55%)

And these are the results we get when we add noise to the coordinates of the best solution.

BOA|Butterfly Optimization Algorithm|50.0|0.9|0.5|0.8|
=============================
5 Hilly's; Func runs: 10000; result: 0.7293075035168879
25 Hilly's; Func runs: 10000; result: 0.5005498140056686
500 Hilly's; Func runs: 10000; result: 0.2754406215727097
=============================
5 Forest's; Func runs: 10000; result: 0.94839177748601
25 Forest's; Func runs: 10000; result: 0.45575888873544235
500 Forest's; Func runs: 10000; result: 0.19904714403024765
=============================
5 Megacity's; Func runs: 10000; result: 0.5
25 Megacity's; Func runs: 10000; result: 0.24523076923076922
500 Megacity's; Func runs: 10000; result: 0.10603076923077022
=============================
Overall score: 3.95976 (44.00%)

A visualization of how the algorithm works on our test functions, as well as on the standard test functions available in the program for free selection and testing.

Hilly

BOA on the Hilly test function

Forest

BOA on the Forest test function

Megacity

BOA on the Megacity test function

Ackley

BOA on the standard Ackley function

Shaffer

BOA on the standard Shaffer function

In the ranking table of the best population-based optimization methods, the BOA algorithm is included for reference based on the test results.

No. AO Description Hilly Hilly
Final
Forest Forest
Final
Megacity (discrete) Megacity
Final
Final
Result
% of
MAX
10 p (5 F) 50 p (25 F) 1,000 p (500 F) 10 p (5 F) 50 p (25 F) 1,000 p (500 F) 10 p (5 F) 50 p (25 F) 1,000 p (500 F)
1 ANS across-neighborhood search 0.94948 0.84776 0.43857 2.23581 1.00000 0.92334 0.39988 2.32323 0.70923 0.63477 0.23091 1.57491 6.134 68.15
2 CLA code lock algorithm (joo) 0.95345 0.87107 0.37590 2.20042 0.98942 0.91709 0.31642 2.22294 0.79692 0.69385 0.19303 1.68380 6.107 67.86
3 AMOm animal migration optimization M 0.90358 0.84317 0.46284 2.20959 0.99001 0.92436 0.46598 2.38034 0.56769 0.59132 0.23773 1.39675 5.987 66.52
4 (P+O)ES (P+O) evolution strategies 0.92256 0.88101 0.40021 2.20379 0.97750 0.87490 0.31945 2.17185 0.67385 0.62985 0.18634 1.49003 5.866 65.17
5 CTA comet tail algorithm (joo) 0.95346 0.86319 0.27770 2.09435 0.99794 0.85740 0.33949 2.19484 0.88769 0.56431 0.10512 1.55712 5.846 64.96
6 TETA time evolution travel algorithm (joo) 0.91362 0.82349 0.31990 2.05701 0.97096 0.89532 0.29324 2.15952 0.73462 0.68569 0.16021 1.58052 5.797 64.41
7 SDSm stochastic diffusion search M 0.93066 0.85445 0.39476 2.17988 0.99983 0.89244 0.19619 2.08846 0.72333 0.61100 0.10670 1.44103 5.709 63.44
8 ECBO enhanced_colliding_bodies_optimization 0.93479 0.75747 0.32471 2.01697 0.97436 0.77446 0.23037 1.97919 0.88923 0.58061 0.15224 1.62208 5.618 62.43
9 BOAm billiards optimization algorithm M 0.95757 0.82599 0.25235 2.03590 1.00000 0.90036 0.30502 2.20538 0.73538 0.52523 0.09563 1.35625 5.598 62.19
10 AAm archery algorithm M 0.91744 0.70876 0.42160 2.04780 0.92527 0.75802 0.35328 2.03657 0.67385 0.55200 0.23738 1.46323 5.548 61.64
11 ESG evolution of social groups (joo) 0.99906 0.79654 0.35056 2.14616 1.00000 0.82863 0.13102 1.95965 0.82333 0.55300 0.04725 1.42358 5.529 61.44
12 SIA simulated isotropic annealing (joo) 0.95784 0.84264 0.41465 2.21513 0.98239 0.79586 0.20507 1.98332 0.68667 0.49300 0.09053 1.27020 5.469 60.76
13 EOm extremal_optimization_M 0.76166 0.77242 0.31747 1.85155 0.99999 0.76751 0.23527 2.00277 0.74769 0.53969 0.14249 1.42987 5.284 58.71
14 BBO biogeography-based optimization 0.94912 0.69456 0.35031 1.99399 0.93820 0.67365 0.25682 1.86867 0.74615 0.48277 0.17369 1.40261 5.265 58.50
15 ACS artificial cooperative search 0.75547 0.74744 0.30407 1.80698 1.00000 0.88861 0.22413 2.11274 0.69077 0.48185 0.13322 1.30583 5.226 58.06
16 DA dialectical algorithm 0.86183 0.70033 0.33724 1.89940 0.98163 0.72772 0.28718 1.99653 0.70308 0.45292 0.16367 1.31967 5.216 57.95
17 BHAm black hole algorithm M 0.75236 0.76675 0.34583 1.86493 0.93593 0.80152 0.27177 2.00923 0.65077 0.51646 0.15472 1.32195 5.196 57.73
18 ASO anarchy society optimization 0.84872 0.74646 0.31465 1.90983 0.96148 0.79150 0.23803 1.99101 0.57077 0.54062 0.16614 1.27752 5.178 57.54
19 RFO royal flush optimization (joo) 0.83361 0.73742 0.34629 1.91733 0.89424 0.73824 0.24098 1.87346 0.63154 0.50292 0.16421 1.29867 5.089 56.55
20 AOSm atomic orbital search M 0.80232 0.70449 0.31021 1.81702 0.85660 0.69451 0.21996 1.77107 0.74615 0.52862 0.14358 1.41835 5.006 55.63
21 TSEA turtle shell evolution algorithm (joo) 0.96798 0.64480 0.29672 1.90949 0.99449 0.61981 0.22708 1.84139 0.69077 0.42646 0.13598 1.25322 5.004 55.60
22 BSA backtracking_search_algorithm 0.97309 0.54534 0.29098 1.80941 0.99999 0.58543 0.21747 1.80289 0.84769 0.36953 0.12978 1.34700 4.959 55.10
23 DE differential evolution 0.95044 0.61674 0.30308 1.87026 0.95317 0.78896 0.16652 1.90865 0.78667 0.36033 0.02953 1.17653 4.955 55.06
24 SRA successful restaurateur algorithm (joo) 0.96883 0.63455 0.29217 1.89555 0.94637 0.55506 0.19124 1.69267 0.74923 0.44031 0.12526 1.31480 4.903 54.48
25 BO bonobo_optimizer 0.77565 0.63805 0.32908 1.74278 0.88088 0.76344 0.25573 1.90005 0.61077 0.49846 0.14246 1.25169 4.895 54.38
26 CRO chemical reaction optimization 0.94629 0.66112 0.29853 1.90593 0.87906 0.58422 0.21146 1.67473 0.75846 0.42646 0.12686 1.31178 4.892 54.36
27 BIO blood inheritance optimization (joo) 0.81568 0.65336 0.30877 1.77781 0.89937 0.65319 0.21760 1.77016 0.67846 0.47631 0.13902 1.29378 4.842 53.80
28 DOA dream_optimization_algorithm 0.85556 0.70085 0.37280 1.92921 0.73421 0.48905 0.24147 1.46473 0.77231 0.47354 0.18561 1.43146 4.825 53.62
29 BSA bird swarm algorithm 0.89306 0.64900 0.26250 1.80455 0.92420 0.71121 0.24939 1.88479 0.69385 0.32615 0.10012 1.12012 4.809 53.44
30 DEA dolphin_echolocation_algorithm 0.75995 0.67572 0.34171 1.77738 0.89582 0.64223 0.23941 1.77746 0.61538 0.44031 0.15115 1.20684 4.762 52.91
31 HS harmony search 0.86509 0.68782 0.32527 1.87818 0.99999 0.68002 0.09590 1.77592 0.62000 0.42267 0.05458 1.09725 4.751 52.79
32 SSG saplings sowing and growing 0.77839 0.64925 0.39543 1.82308 0.85973 0.62467 0.17429 1.65869 0.64667 0.44133 0.10598 1.19398 4.676 51.95
33 BCOm bacterial chemotaxis optimization M 0.75953 0.62268 0.31483 1.69704 0.89378 0.61339 0.22542 1.73259 0.65385 0.42092 0.14435 1.21912 4.649 51.65
34 ABO african buffalo optimization 0.83337 0.62247 0.29964 1.75548 0.92170 0.58618 0.19723 1.70511 0.61000 0.43154 0.13225 1.17378 4.634 51.49
35 (PO)ES (PO) evolution strategies 0.79025 0.62647 0.42935 1.84606 0.87616 0.60943 0.19591 1.68151 0.59000 0.37933 0.11322 1.08255 4.610 51.22
36 FBA fractal-based algorithm 0.79000 0.65134 0.28965 1.73099 0.87158 0.56823 0.18877 1.62858 0.61077 0.46062 0.12398 1.19537 4.555 50.61
37 TSm tabu search M 0.87795 0.61431 0.29104 1.78330 0.92885 0.51844 0.19054 1.63783 0.61077 0.38215 0.12157 1.11449 4.536 50.40
38 BSO brain storm optimization 0.93736 0.57616 0.29688 1.81041 0.93131 0.55866 0.23537 1.72534 0.55231 0.29077 0.11914 0.96222 4.498 49.98
39 WOAm whale optimization algorithm M 0.84521 0.56298 0.26263 1.67081 0.93100 0.52278 0.16365 1.61743 0.66308 0.41138 0.11357 1.18803 4.476 49.74
40 AEFA artificial electric field algorithm 0.87700 0.61753 0.25235 1.74688 0.92729 0.72698 0.18064 1.83490 0.66615 0.11631 0.09508 0.87754 4.459 49.55
41 AEO artificial ecosystem-based optimization algorithm 0.91380 0.46713 0.26470 1.64563 0.90223 0.43705 0.21400 1.55327 0.66154 0.30800 0.28563 1.25517 4.454 49.49
42 CAm camel algorithm M 0.78684 0.56042 0.35133 1.69859 0.82772 0.56041 0.24336 1.63149 0.64846 0.33092 0.13418 1.11356 4.444 49.37
43 ACOm ant colony optimization M 0.88190 0.66127 0.30377 1.84693 0.85873 0.58680 0.15051 1.59604 0.59667 0.37333 0.02472 0.99472 4.438 49.31
44 CMAES covariance_matrix_adaptation_evolution_strategy 0.76258 0.72089 0.00000 1.48347 0.82056 0.79616 0.00000 1.61672 0.75846 0.49077 0.00000 1.24923 4.349 48.33
45 BOA butterfly_optimization_algorithm 0.72930 0.50054 0.27544 1.50528 0.94839 0.45575 0.19904 1.60318 0.50000 0.24523 0.10603 0.85126 3.960 44.00
RW random walk 0.48754 0.32159 0.25781 1.06694 0.37554 0.21944 0.15877 0.75375 0.27969 0.14917 0.09847 0.52734 2.348 26.09


Conclusions

The BOA algorithm is based on the idea that each butterfly is both a seeker and a source of fragrance, the strength of which is determined by the quality of the solution found. The fragrance formula, based on Stevens’ law from psychophysics, ensures a smooth transition from uniform exploration of space in early iterations to targeted exploitation of promising areas in later ones.

While working on the algorithm, an error was discovered and corrected in the original movement formulas published by the authors. Incorrect parenthesization in the global and local search expressions led to a systematic shift of the population toward the origin instead of movement toward the optimum. This error went unnoticed because standard test functions typically have their optimum precisely at the origin, and the flaw in the formulas happened to coincide with the direction toward the solution. After the parentheses were corrected, the algorithm worked correctly on problems with an arbitrarily located optimum in the search space.

Nevertheless, even with the corrected formulas, the algorithm exhibited a tendency toward premature convergence. As the population converges, the differences between the butterflies’ positions decrease, the local search steps become negligibly small, and the population loses the ability to leave the neighborhood of the current best solution. To compensate for this shortcoming, a stochastic mutation mechanism — which was absent from the original description — was added: with a small probability, one of the butterfly’s coordinates is replaced with a value generated from a normal distribution centered on the best found solution. This allowed us to improve the algorithm performance from 32.5% to 44% on the test function set.

Despite the improvement, the result achieved is not sufficient for the algorithm to enter the ranking table. The main reason lies in the architectural limitations of the approach itself. The butterfly movement mechanism relies on only two modes: a global mode, directed toward the single best solution, and a local mode, determined by the difference between the positions of two random agents. Both modes generate steps scaled by the square of a random number r², with a mean value of approximately 0.33; when combined with a fragrance factor of about 0.5, this results in a typical step size of about 12% of the distance to the target. Such dynamics are too conservative for effectively exploring multidimensional spaces with many local optima.

Overall, the butterfly optimization algorithm should be regarded as a viable, though not outstanding, metaheuristic. Its strength lies in its ease of implementation. A weakness remains the lack of diversity in search strategies and the tendency for steps to decay as the population converges.

tab

Figure 2. Color gradation of algorithms for the corresponding tests

chart

Figure 3. Histogram of algorithm test results (on a scale from 0 to 100; the higher the score, the better, where 100 is the maximum possible theoretical result; the archive contains a script for calculating the ranking table)


Pros and Cons of the BOA Algorithm:

Pros:

  1. Decent performance on Forest-type problems (smooth problems with "sharp" extrema).

Cons:

  1. Prone to getting stuck.

An archive containing the latest versions of the algorithm source code is attached to the article. The author of this article does not guarantee absolute accuracy in the description of the canonical algorithms; many of them have been modified to improve their search capabilities. The conclusions and assessments presented in the articles are based on the results of the experiments conducted.

Updates to the unified test bench: fixed the algorithm selection enumerator (corrected the selection of the following algorithms):

  • AO_CryStAlm (Crystal Structure Algorithm M)
  • AO_CoSO (Community of Scientist Optimization)



Code files used in this article

# Name Type Description
1 #C_AO.mqh
Included file
Parent class of population-based optimization algorithms
2 #C_AO_enum.mqh
Included file
Enumeration of population-based optimization algorithms
3 TestFunctions.mqh
Included file
Test function library
4
TestStandFunctions.mqh
Included file
Test bench function library
5
Utilities.mqh
Included file
Utility functions library
6
CalculationTestResults.mqh
Included file
Script for calculating results for the comparison table
7
Testing AOs.mq5
Script Unified test bench for all population-based optimization algorithms
8
Simple Use of Population Optimization Algorithms.mq5
Script
A simple example of using population-based optimization algorithms without visualization
9
Test_BOA.mq5
Script Test bench for BOA


Translated from Russian by MetaQuotes Ltd.
Original article: https://www.mql5.com/ru/articles/21209

Attached files |
BOA.zip (340.31 KB)
Features of Custom Indicators Creation Features of Custom Indicators Creation
Creation of Custom Indicators in the MetaTrader trading system has a number of features.
Neural Networks in Trading: The Temporal Query Model (TQNet) Neural Networks in Trading: The Temporal Query Model (TQNet)
The TQNet framework opens up new possibilities for modeling and forecasting financial time series by combining modularity, flexibility, and high performance. The article explores the possibility of implementing complex mechanisms for handling global correlations, including advanced parameter initialization methods.
Features of Experts Advisors Features of Experts Advisors
Creation of expert advisors in the MetaTrader trading system has a number of features.
Developing a Multi-Currency Expert Advisor (Part 30): From Trading Strategy to Launching a Multi-Currency Expert Advisor Developing a Multi-Currency Expert Advisor (Part 30): From Trading Strategy to Launching a Multi-Currency Expert Advisor
The article outlines the complete process of creating a multi-currency Expert Advisor using the Adwizard library for MetaTrader 5: from setting up the environment for creating optimization projects to obtaining the final multi-currency Expert Advisors, which combine multiple instances of a simple trading strategy. We will walk through setting up the necessary input parameters, conventions for convenient file names, and launching three instances of the final Expert Advisors on different trading accounts with different parameters.