Extremal Optimization (EO)

Contents

1. Introduction
2. Implementation of the algorithm
3. Test results

Introduction

Many real-world problems, particularly trading ones, are characterized by complex discrete objective function landscapes with multiple local extrema, discontinuities, and non-differentiable regions, making classical gradient-based methods inapplicable. Numerous metaheuristic algorithms have been developed to solve such problems, and each approach has its own advantages and disadvantages in balancing exploration and exploitation of the search space.

Extremal Optimization (EO) is a metaheuristic optimization algorithm inspired by the Bak-Sneppen model. The algorithm was developed by Stefan Boettcher and Allon Percus in 1999 as a method inspired by the concept of self-organized criticality, according to which complex systems naturally evolve toward a critical state where avalanche-like changes of different scales occur. A population-based variant of EO was developed for continuous optimization using iterative population-level updates.

Implementation of the algorithm

The Bak-Sneppen model was developed to describe the evolution of ecosystems and demonstrate the phenomenon of self-organized criticality. This is a simple coevolution model that shows how local interactions can lead to global critical events.

Principles of the Bak-Sneppen model

Ecosystem of species:

• N species are arranged in a chain (or on a lattice)
• Each species i is assigned a fitness value fi ∈ [0,1]
• Fitness represents the adaptability of a species to its environment.

Dynamics of evolution

At each time step:
1. Find the species with minimum fitness: fmin = min{fi}
2. Replace the fitness of this species and its neighbors with new random values
3. Repeat the process

Self-organized criticality

• The system self-organizes to a critical state
• The threshold value fc ≈ 0.67 occurs
• Species with fitness < fc die out in avalanche-like events
• Avalanche sizes follow a power law

Long periods of stagnation, when few changes are followed by sudden avalanche-like changes of varying scales in the absence of external control, the system self-organizes. Power-law changes affect the sizes of avalanches: P(s) ~ s^(-τ), where the lifetimes of species obey a power law in the absence of a characteristic scale. Changes in one species affect its neighbors. Global dynamics emerge from local interactions.

Boettcher and Perkus adapted the principles of the Bak-Sneppen model to create an optimization algorithm. Instead of improving the best, we improve the worst. Focusing on the worst elements is counterintuitive, but effective for breaking out of local optima. In the Bak-Sneppen model, evolution occurs through "punctuated equilibrium"—long periods of stasis interrupted by sudden changes. EO uses this principle: when only bad components are changed, and when changing one component dramatically improves the solution, improving one component may make other components "bad".

Extremal Optimization applies the principles of self-organized criticality — originally studied in ecosystem models—to optimization. The key innovation is the use of the natural tendency of complex systems to self-organize to automatically balance the exploration and exploitation of the search space. We will test this claim on practical problems. Let us now outline the algorithm in pseudocode.

Initialization:

1. Create a population out of popSize agents 
2. Initialize parameters
3. Create arrays for ranking agents and components

First run (population initialization):

1. For each agent in the population:
   • For each component (coordinates):
     • assign a random value within the acceptable range
     • apply discretization if necessary

Basic optimization loop:

For each agent in the population, execute:

1. Agent ranking:
   • create a list of all agents with their fitness
   • sort agents by fitness (worst first for maximization)
2. Select a target agent:
   • use the power-law distribution P(n) ∝ n^(-τ)
   • choose a rank according to this distribution
   • define the agent with the selected rank as the target one
3. Target agent's component ranking:
   • For each component, calculate its fitness:
     • fitness = 1 - (normalized deviation from the best known value)
   • Sort components by fitness (worst first)
4. Component selection and mutation:
   • use a power-law distribution to choose the component rank
   • replace the selected component with a new random value
   • check the boundaries and apply discretization

Results update:

1. Sort the entire population by fitness (the best at the beginning)
2. Update the global best solution if an improvement is found

Let's implement the algorithm in code. We will write a class that will represent the implementation of the EO algorithm, based on the optimization method using the principle of removing the worst elements to find solutions. The class inherits from the base class C_AO and additionally includes algorithm-specific parameters and methods. Constructor and destructor: initializing parameters and releasing resources respectively. Main parameters

• popSize — population size, number of agents in the population;
• tau — power-law distribution that determines the probabilities of selecting elements;
• greedyStart — proportion of agents with greedy initialization;
• eliteUpdate — proportion of agents participating in the update per iteration.

The params parameter array contains settings and copies the values of the original parameters for integration with the user interface or external control. Methods:

• SetParams — extract the current parameter values from the array and apply them;
• Init — initialize the population taking into account the ranges of parameters and the number of epochs;
• Moving — execute one step of the algorithm updating the state of the population;
• Revision — updates evaluations or prepare for the next stage. 

Internal structures:

• RankedComponent — structure for storing the ranking of agent components by their quality (Fitness);
• RankedAgent — structure for storing the ranking of the agents themselves according to their overall fitness.

Internal variables and parameters:

• compRanks and agentRanks — arrays of structures for storing sorted components and agents;
• tau, greedyStart, eliteUpdate control the behavior of the algorithm

Internal methods:

• ApplyExtremalOptimization — basic mechanism for selecting and updating elements based on their ranks, taking into account the power law distribution.
• CalculateComponentFitness — calculate the fitness of the agent component.
• SelectRankByPowerLaw — select the rank of a component or agent using a power law distribution, which allows focusing on the worst elements.

The class is intended for the application of the extreme optimization algorithm, especially in problems requiring the search for a high-quality solution by successively eliminating the least fit components in the population.

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

  //------------------------------------------------------------------
  ArrayResize (compRanks, coords);
  ArrayResize (agentRanks, popSize);

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

The initialization method of the EO class is responsible for preparing the initial conditions for further execution. It accepts parameters for the minimum and maximum value ranges for the variables, as well as step sizes for changing them, and the number of epochs (iterations) if needed.

The method first calls the standard initialization procedure, which configures basic parameters and prepares common data structures. If this procedure fails, initialization is aborted.

If the initialization is successful, the method resizes the arrays used to store the ranked components and agents. The size of the component array is set according to the number of coordinates or variables of the problem, and the array of agents is set according to the given population size. Once these steps are completed, the method returns the result of successful initialization.

//————————————————————————————————————————————————————————————————————
//--- Main loop of the algorithm
void C_AO_EO::Moving ()
{
  // Initial population setup
  if (!revision)
  {
    int greedyCount = (int)(popSize * greedyStart);

    for (int i = 0; i < popSize; i++)
    {
      // Random initialization for the rest
      for (int c = 0; c < coords; c++)
      {
        a [i].c [c] = u.RNDfromCI (rangeMin [c], rangeMax [c]);
        a [i].c [c] = u.SeInDiSp (a [i].c [c], rangeMin [c], rangeMax [c], rangeStep [c]);
      }
    }

    revision = true;
    return;
  }

  // Apply Extremal Optimization ---------------------------------
  ApplyExtremalOptimization ();
}
//————————————————————————————————————————————————————————————————————

The ApplyExtremalOptimization method implements the main loop of the EO algorithm for updating the agent population. Its goal is to improve solutions by modifying agent components based on EO principles. The number of agents to update per iteration is determined taking into account the eliteUpdate parameter (the proportion of agents to update) and the population size. Agents are sorted based on their fitness. The algorithm ranks agents from worst to best so that the algorithm is designed for maximization. Using a power-law probability distribution, the agent that will be subject to change is selected.

For the selected agent, its components are ranked. Each component is assigned a fitness value to evaluate its contribution to the agent fitness. Components are ranked from worst to best based on their "fitness". A power-law probability distribution is used to select one of the agent's components for modification. This is done to make it more likely that components with the worst rank are selected.

The selected agent component is replaced with a random value within the allowed range, which is specified by the rangeMin and rangeMax parameters. The resulting value is checked for validity to ensure that it complies with the specified boundaries and step sizes.

The method repeats these steps for each agent to be updated, implementing the key logic of the EO algorithm: selecting components with the worst performance and randomly replacing them in order to improve the fitness of solutions.

//————————————————————————————————————————————————————————————————————
//--- Apply Extremal Optimization
void C_AO_EO::ApplyExtremalOptimization ()
{
  // Number of agents to update in this iteration
  int numUpdates = MathMax (1, (int)(popSize * eliteUpdate));

  // Update selected agents based on EO principle
  //for (int update = 0; update < numUpdates; update++)
  for (int update = 0; update < popSize; update++)
  {
    // Step 1: Select an agent to modify
    // Use ranking by overall fitness
    int targetAgent;

    // Rank agents by fitness (from worst to best for maximization)
    for (int i = 0; i < popSize; i++)
    {
      agentRanks [i].idx = i;
      agentRanks [i].fitness = a [i].f;
    }

    // Sort (worst first to maximize)
    for (int i = 0; i < popSize - 1; i++)
    {
      for (int j = i + 1; j < popSize; j++)
      {
        if (agentRanks [i].fitness > agentRanks [j].fitness)
        {
          RankedAgent temp = agentRanks [i];
          agentRanks [i] = agentRanks [j];
          agentRanks [j] = temp;
        }
      }
    }

    // Select an agent according to the power-law distribution
    int rank    = SelectRankByPowerLaw (popSize);
    targetAgent = agentRanks [rank].idx;

    // Step 2: Rank the components of the selected agent
    for (int c = 0; c < coords; c++)
    {
      compRanks [c].agentIdx     = targetAgent;
      compRanks [c].componentIdx = c;
      compRanks [c].fitness      = CalculateComponentFitness (targetAgent, c);
    }

    // Sort components (worst first)
    for (int i = 0; i < coords - 1; i++)
    {
      for (int j = i + 1; j < coords; j++)
      {
        if (compRanks [i].fitness > compRanks [j].fitness)
        {
          RankedComponent temp = compRanks [i];
          compRanks [i] = compRanks [j];
          compRanks [j] = temp;
        }
      }
    }

    // Step 3: Select a component to change according to P(n) ∝ n^(-τ)
    int compRank = SelectRankByPowerLaw (coords);
    int compIdx  = compRanks [compRank].componentIdx;

    // Step 4: Replace the selected component with a new random value
    // This is the key principle of EO - unconditional replacement with random
    a [targetAgent].c [compIdx] = u.RNDfromCI (rangeMin [compIdx], rangeMax [compIdx]);

    // Check boundaries
    a [targetAgent].c [compIdx] = u.SeInDiSp (a [targetAgent].c [compIdx],
                                              rangeMin [compIdx],
                                              rangeMax [compIdx],
                                              rangeStep [compIdx]);
  }
}
//————————————————————————————————————————————————————————————————————

The CalculateComponentFitness method calculates the fitness value for a single component of the agent. This function plays an important role in the Extremal Optimization (EO) algorithm, evaluating the contribution of each component to the fitness of the solution. The algorithm uses a simple metric based on calculating "fitness" as the relative contribution of a component to the fitness.

Initially, "fitness" is initialized to zero. Next, the range of component values (range) is calculated. If this range is greater than zero, the normalized deviation is calculated. The deviation is calculated as the absolute value of the difference between the agent component value and the base value of cB, normalized by the range of component values. Fitness is calculated by subtracting "deviation" from 1.0. This makes "fitness" inversely proportional to the deviation: the smaller the deviation of a component from the target value, the higher its "fitness".

Thus, the function returns a fitness score for a particular component, which indicates how well the current value of the component matches the target value, which is necessary for further ranking of components.

//————————————————————————————————————————————————————————————————————
//--- Calculate the fitness component
double C_AO_EO::CalculateComponentFitness (int agentIdx, int componentIdx)
{
  // For the general optimization problem we use a simple metric
  // λi = relative contribution of the component to the fitness

  double fitness = 0.0;

  double range = rangeMax [componentIdx] - rangeMin [componentIdx];
  if (range > 0)
  {
    // Normalized deviation
    double deviation = MathAbs (a [agentIdx].c [componentIdx] - cB [componentIdx]) / range;
    fitness =1.0 - deviation;// Invert so that bigger = better
  }

  return fitness;
}
//————————————————————————————————————————————————————————————————————

The SelectRankByPowerLaw method implements the selection of the position (rank) of an element from a list based on a power probability distribution. This method is key to the EO algorithm, since it is through it that the selection of agents and components for modification is made. The method takes maxRank as input, which defines the maximum rank of the element.

The method is based on the principle of inverse transformation to obtain a random value, following the power law P(n) ∝ n^(-τ), where: n is the rank of the element; τ (tau) is the parameter determining the shape of the power-law distribution.

Depending on the value of the "tau" parameter, two scenarios are realized:

1. General case (τ ≠ 1.0):

   • a normalizing factor "norm" is calculated, which ensures the correct scaling of probabilities;
   • 'r' random number in the range [0, 1) is generated using u.RNDprobab();
   • the inverse transform formula is used to calculate the "rank" based on "r", "norm" and "τ";
   • the rank is limited to [0, maxRank - 1] to avoid array out-of-bounds.

2. Special case (τ = 1.0):

   • the normalizing factor "norm" is calculated for the case τ = 1.0;
   • a special formula is used to calculate the "rank";
   • the rank is limited within the range [0, maxRank - 1].

As a result, the method returns an integer "rank", which indicates the selected position of the element in the ranked list. This selection is probabilistic, that is, elements with lower rank (worst) have a higher probability of being selected, in accordance with the principles of EO.

//————————————————————————————————————————————————————————————————————
//--- Select a rank according to a power-law distribution
int C_AO_EO::SelectRankByPowerLaw (int maxRank)
{
  // P(n) ∝ n^(-τ), where n is a rank from 1 to maxRank
  // Use the inverse transformation method

  double r = u.RNDprobab ();

  if (tau != 1.0)
  {
    // General case: inverse transform for P(n) ∝ n^(-τ)
    double norm = (1.0 - MathPow (maxRank + 1.0, 1.0 - tau)) / (1.0 - tau);
    double x = r * norm;
    int rank = (int)MathPow ((1.0 - tau) * x + 1.0, 1.0 / (1.0 - tau)) - 1;

    if (rank >= maxRank) rank = maxRank - 1;
    if (rank < 0) rank = 0;

    return rank;
  }
  else
  {
    // Special case τ = 1: P(n) ∝ 1/n
    double norm = MathLog (maxRank + 1.0);
    int rank = (int)(MathExp (r * norm) - 1.0);

    if (rank >= maxRank) rank = maxRank - 1;
    if (rank < 0) rank = 0;

    return rank;
  }
}
//————————————————————————————————————————————————————————————————————

The "Revision" method is designed to update information about the best solutions found during the operation of the EO algorithm. Its task is to determine the current best solution and update global variables that store information about it. A temporary array of aT is created to store the sorted agents. The 'a' population of agents is sorted using the built-in u.Sorting function. The sorting is performed to maximize, that is, the best agents (those with the highest 'f' fitness function value) are placed at the beginning of the 'a' array (which is the main data structure) containing information about the agents, including the values of their 'c' components and the value of the 'f' fitness function.

After sorting, it is checked whether the 'f' fitness function value of the best agent (i.e. a[0]) is greater than the current best fB value (stores the fitness function value of the best solution found so far). If the current best agent is better than the global best solution, then the values of the best agent's components (a[0].c) are copied into the cB array, which likely stores the components of the best solution, and fB is updated to the value of the best agent's fitness function (a[0].f).

Thus, the Revision function updates the information about the best solution, preserving its components and the value of the fitness function, and sorts the population so that the best solutions are always available at the beginning of the array. 

//————————————————————————————————————————————————————————————————————
//--- Update the best solutions
void C_AO_EO::Revision ()
{
  // Sort the population for MAXIMIZATION
  static S_AO_Agent aT [];
  ArrayResize (aT, popSize);

  // Use the built-in sorting function
  u.Sorting (a, aT, popSize);

  // Update the global best solution
  if (a [0].f > fB)
  {
    ArrayCopy (cB, a [0].c, 0, 0, WHOLE_ARRAY);
    fB = a [0].f;
  }
}
//————————————————————————————————————————————————————————————————————

Test results

After testing, the algorithm did not show good results, so I decided to review all the original methods and make a modified version aimed at improving the convergence of the algorithm.

EO|Extremal Optimization (Boettcher-Percus)|50.0|1.4|0.5|0.3|
=============================
5 Hilly's; Func runs: 10000; result: 0.5146369337195529
25 Hilly's; Func runs: 10000; result: 0.29089804555433085
500 Hilly's; Func runs: 10000; result: 0.25192095557138877
=============================
5 Forest's; Func runs: 10000; result: 0.367128650332966
25 Forest's; Func runs: 10000; result: 0.19477408852361866
500 Forest's; Func runs: 10000; result: 0.15367465708144543
=============================
5 Megacity's; Func runs: 10000; result: 0.2584615384615384
25 Megacity's; Func runs: 10000; result: 0.12707692307692314
500 Megacity's; Func runs: 10000; result: 0.09413846153846227
=============================
All score: 2.25271 (25.03%)

Below is my own interpretation of this optimization method. We implement a new class that inherits the base optimization class. 

The constructor sets the initial parameters of the algorithm: population size, the number of worst solutions to be increased (improved), the probability of mutation, and the degrees of distribution for selection and mutation. The parameters are stored in an array of structures for further modification and customization.

The SetParams method updates the class's internal variables based on the values in the parameters array. Methods are declared for initializing the algorithm, performing movement steps (updating the state), and revision (reevaluates solutions or prepares the next iteration). The private part of the class stores variables that track the current epoch and the total number of epochs of the algorithm's execution. A method for mutating a specific component of a solution in a population is defined, which is a key part of the implemented algorithm. 

//————————————————————————————————————————————————————————————————————
class C_AO_EOm : public C_AO
{
  public: //----------------------------------------------------------
  ~C_AO_EOm () { }
  C_AO_EOm ()
  {
    ao_name = "EOm";
    ao_desc = "Extremal Optimization M";
    ao_link = "https://www.mql5.com/en/articles/18755";

    popSize        = 50;      // Population size
    popRaising     = 3;       // Boost the worst agents
    mutationRate   = 0.1;     // Mutation probability
    powCh          = 2.0;     // Selection power-law exponent
    powMut         = 8.0;     // Mutation power-law exponent 

    ArrayResize (params, 5);

    params [0].name = "popSize";        params [0].val = popSize;
    params [1].name = "popRaising";     params [1].val = popRaising;
    params [2].name = "mutationRate";   params [2].val = mutationRate;
    params [3].name = "powCh";          params [3].val = powCh;
    params [4].name = "powMut";         params [4].val = powMut;
  }

  void SetParams ()
  {
    popSize        = (int)params [0].val;
    popRaising     = (int)params [1].val;
    mutationRate   = params      [2].val;
    powCh          = params      [3].val;
    powMut         = params      [4].val;
  }

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

  void Moving   ();
  void Revision ();

  //------------------------------------------------------------------
  int    popRaising;       // Boost the worst agents
  double mutationRate;     // Mutation probability
  double powCh;            // Selection power-law exponent
  double powMut;           // Mutation power-law exponent

  private: //---------------------------------------------------------
  int    currentEpoch;     // current epoch
  int    totalEpochs;      // total number of epochs

  void MutateComponent (int agentIdx, int componentIdx);
};
//————————————————————————————————————————————————————————————————————

The Init method is responsible for initializing the algorithm. The function takes as input arrays with minimum, maximum values and step for optimization parameters, as well as the total number of epochs (cycles) of the algorithm operation.

Inside the method, the StandardInit function from the C_AO base class is first called. This function performs general initialization, memory allocation, or setting up basic parameters related to ranges and steps for optimization parameters.

If StandardInit fails, then the Init function also exits, returning 'false'. If StandardInit is successful, the method initializes the variables currentEpoch (current epoch) to 0 and totalEpochs (total number of epochs) to the value received in the epochsP parameter.

At the end, the method returns 'true', signaling successful initialization. Thus, the method prepares the algorithm for work by setting the initial parameters and epoch counters.

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

  //------------------------------------------------------------------
  currentEpoch = 0;
  totalEpochs  = epochsP;

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

The Moving method represents the main cycle of the extreme optimization algorithm. It is responsible for the evolution of the solution population in the process of searching for the optimal solution. Epoch increment: the currentEpoch counter is incremented. Population initialization (first epoch only): if the 'revision' flag (indicating the need to revise the decisions) is equal to 'false' (i.e. in the first cycle of the algorithm), the population is initialized.

For each "agent" (solution) in the population, for each "coordinate" (solution parameter), a random value is generated in a given range (rangeMin, rangeMax) for the current coordinate. The SeInDiSp discretization is used to ensure that the parameter values correspond to the specified steps. After initialization, "revision" is set to 'true' so that this block of code is executed only once. The method terminates. 

Applying extremal optimization in subsequent epochs: a temporary array aT is created to store new (mutated) solutions. For each agent in the population, the aT[i] structure is initialized. For each solution coordinate, a random number is generated and raised to the power of powCh. This is used to select the "parent".

The random number is scaled to select the parent index in the ind population. The mutation type is determined: if the random number is less than mutationRate (the mutation is enabled with mutationRate probability), then the mutation is applied according to the PowerDistribution power-law distribution law to the value of the parent coordinate. Otherwise, there is a "directed movement" towards the best solution with the addition of random noise. The new value is calculated as the sum of the parent coordinate value, a random number, and the difference between cB[c] (the best coordinate found) and the parent coordinate. The SeInDiSp discretization is used. After mutation, all coordinates from the temporary array aT are copied to the 'a' main population array.

Thus, the method implements an evolution cycle, where at each step (epoch) the solutions mutate (change), using the mechanisms of selection, random mutation and directed search, with the aim of improving their values (minimizing or maximizing the objective function) in order to ultimately find the optimal solution.

//————————————————————————————————————————————————————————————————————
//--- Main loop of the algorithm
void C_AO_EOm::Moving ()
{
  currentEpoch++;

  // Initial population setup
  if (!revision)
  {
    for (int i = 0; i < popSize; i++)
    {
      for (int c = 0; c < coords; c++)
      {
        a [i].c [c] = u.RNDfromCI (rangeMin [c], rangeMax [c]);
        a [i].c [c] = u.SeInDiSp (a [i].c [c], rangeMin [c], rangeMax [c], rangeStep [c]);
      }
    }

    revision = true;
    return;
  }

  //Apply Extremal Optimization---------------------------------------
  static S_AO_Agent aT []; ArrayResize (aT, popSize);

  for (int i = 0; i < popSize; i++)
  {
    aT [i].Init (coords);

    for (int c = 0; c < coords; c++)
    {
      double rnd = u.RNDprobab (); rnd = pow (rnd, powCh);
      int ind = (int)u.Scale (rnd, 0.0, 1.0, 0, popSize - 1);

      // Select the mutation type
      double mutType = u.RNDprobab ();

      if (mutType < mutationRate)
      {
        aT [i].c [c] = u.PowerDistribution (a [ind].c [c], rangeMin [c], rangeMax [c], powMut);
      }
      else
      {
        // Directed movement towards the better with noise
        aT [i].c [c] = a [ind].c [c] + u.RNDprobab () * (cB [c] - a [ind].c [c]);
      }

      // Check boundaries
      aT [i].c [c] = u.SeInDiSp (aT [i].c [c], rangeMin [c], rangeMax [c], rangeStep [c]);
    }
  }

  for (int i = 0; i < popSize; i++) ArrayCopy (a [i].c, aT [i].c);
}
//————————————————————————————————————————————————————————————————————

The method is responsible for updating information about the best solutions during the algorithm's operation. 

Sorting the population: A temporary array of aT is created to store a copy of the population. The population 'a' is sorted by the value of the objective function. The u.Sorting method is used for sorting. 

Update the global best solution: it is checked whether the value of the objective function "f" of the first element of the population (i.e. the best solution) is better than the current global best value of fB. If the current one is better, then cB is copied from the coordinates of the first element of a[0].c, and fB is updated with the value of f of the better solution. 

Worst solution update: fW is updated with the value of the function of the last element of the population. 

Raising the worst solutions: in a loop executed popRaising times, the worst solutions in the population are raised: for popRaising the worst agents (located at the end of the sorted array of 'a'), the value of the objective function 'f' is recalculated as a random value located between fW (the worst current result) and fB (the best current result). This is done to introduce variety or to shift the worst solutions towards the better. 

Resorting the population: population 'a' is re-sorted using the u.Sorting method. This is necessary to bring the population back into order after changes in the 'f' values.

//————————————————————————————————————————————————————————————————————
//--- Update the best solutions
void C_AO_EOm::Revision ()
{
  // Sort the population  --------------------------------------------
  static S_AO_Agent aT []; ArrayResize (aT, popSize);
  u.Sorting (a, aT, popSize);

  // Update the global best solution
  if (a [0].f > fB)
  {
    ArrayCopy (cB, a [0].c, 0, 0, WHOLE_ARRAY);
    fB = a [0].f;
  }

  fW = a [popSize - 1].f;

  //------------------------------------------------------------------
  for (int i = 0; i < popRaising; i++)
  {
    a [popSize - 1 - i].f = u.RNDfromCI (fW, fB);
  }

  u.Sorting (a, aT, popSize);
}
//————————————————————————————————————————————————————————————————————

Now the algorithm (in some sense inheriting the ideas of EO, albeit in a modified way), uses a power-law distribution to select donors (not necessarily bad ones).

a [popSize - 1 - i].f = u.RNDfromCI (fW, fB);

I have added a feature not found in EO descriptions: the algorithm essentially "resurrects" the worst agents, assigning them a random fitness in the range from the worst to the best in the population. PowerDistribution for mutation is an interpretation of the "power law distribution mutation". Directed movement towards the better (for mutType >= mutationRate ) — addition to accelerate convergence.

Now we can look at the results of the modified version. We can include the following results in our rating table.

EOm|Extremal Optimization Mod|50.0|3.0|0.1|2.0|8.0|
=============================
5 Hilly's; Func runs: 10000; result: 0.7616651321732368
25 Hilly's; Func runs: 10000; result: 0.7724295992586084
500 Hilly's; Func runs: 10000; result: 0.3174703398632668
=============================
5 Forest's; Func runs: 10000; result: 0.9999936711143977
25 Forest's; Func runs: 10000; result: 0.7675163269928252
500 Forest's; Func runs: 10000; result: 0.23527203643380376
=============================
5 Megacity's; Func runs: 10000; result: 0.7476923076923077
25 Megacity's; Func runs: 10000; result: 0.5396923076923076
500 Megacity's; Func runs: 10000; result: 0.14249230769230903
=============================
All score: 5.28422 (58.71%)

The visualization shows a large scatter of results for low-dimensional functions (green lines).

EOm on the Hilly test function

EOm on the Forest test function

EOm on the Megacity test function

After testing, the EOm algorithm ranks 12th in our ranking table.

#	AO	Description	Hilly			Hilly Final	Forest			Forest Final	Megacity (discrete)			Megacity Final	Final Result	% of MAX
			10 p (5 F)	50 p (25 F)	1000 p (500 F)		10 p (5 F)	50 p (25 F)	1000 p (500 F)		10 p (5 F)	50 p (25 F)	1000 p (500 F)
1	ANS	across neighbourhood 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 ptimization 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	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
9	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
10	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
11	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
12	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
13	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
14	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
15	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
16	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
17	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
18	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
19	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
20	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
21	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
22	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
23	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
24	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
25	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
26	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
27	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
28	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
29	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
30	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
31	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
32	(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
33	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
34	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
35	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
36	WOAm	wale 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
37	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
38	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
39	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
40	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
41	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
42	BFO-GA	bacterial foraging optimization - ga	0.89150	0.55111	0.31529	1.75790	0.96982	0.39612	0.06305	1.42899	0.72667	0.27500	0.03525	1.03692	4.224	46.93
43	SOA	simple optimization algorithm	0.91520	0.46976	0.27089	1.65585	0.89675	0.37401	0.16984	1.44060	0.69538	0.28031	0.10852	1.08422	4.181	46.45
44	ABHA	artificial bee hive algorithm	0.84131	0.54227	0.26304	1.64663	0.87858	0.47779	0.17181	1.52818	0.50923	0.33877	0.10397	0.95197	4.127	45.85
45	ACMO	atmospheric cloud model optimization	0.90321	0.48546	0.30403	1.69270	0.80268	0.37857	0.19178	1.37303	0.62308	0.24400	0.10795	0.97503	4.041	44.90
	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

Summary

The presented compact modified implementation of EOm demonstrates an interesting example in the field of metaheuristic optimization: a significant deviation from theoretical principles can lead to improved practical results. The algorithm, which ranked 12th among 45 population-based methods, is actually a hybrid that retains only the key idea of the power-law distribution from the original EO. 

"Resurrection" mechanism of the worst agents prevents premature convergence, maintain population diversity, and create additional opportunities for exploration. Simplification of the computational structure: abandoning component-wise evaluation reduces computational complexity, and working directly with solutions speeds up convergence.

The success of the modified version of EO confirms an important principle in the design of metaheuristics: the efficiency of an algorithm is determined not by its fidelity to the original idea, but by the balance between exploration and exploitation of the search space . The presented implementation, despite its departure from the classical principles of Extremal Optimization, demonstrates high convergence speed, competitive results, and ease of implementation and configuration.

This makes it a useful tool in the arsenal of population optimization methods, especially for problems where the balance between solution quality and computational cost is important.

Figure 1. Color gradation of algorithms according to the corresponding tests

Figure 2. Histogram of algorithm testing results (scale from 0 to 100, the higher the better, where 100 is the maximum possible theoretical result, in the archive there is a script for calculating the rating table)

EOm pros and cons:

Pros:

1. Simple implementation
2. Fast and efficient
3. Good results on discrete problems

Cons:

1. Large scatter of results on low-dimension functions
2. Average results on "smooth" low-dimensional problems

The article is accompanied by an archive with the current versions of the algorithm codes. The author of the article is not responsible for the absolute accuracy in the description of canonical algorithms. Changes have been made to many of them to improve search capabilities. The conclusions and judgments presented in the articles are based on the results of the experiments.

• github: https://github.com/JQSakaJoo/Population-optimization-algorithms-MQL5

Programs used in the article

#	Name	Type	Description
1	#C_AO.mqh	Include	Parent class of population optimization algorithms
2	#C_AO_enum.mqh	Include	Enumeration of population optimization algorithms
3	TestFunctions.mqh	Include	Library of test functions
4	TestStandFunctions.mqh	Include	Test stand function library
5	Utilities.mqh	Include	Library of auxiliary functions
6	CalculationTestResults.mqh	Include	Script for calculating results in the comparison table
7	Testing AOs.mq5	Script	The unified test stand for all population optimization algorithms
8	Simple use of population optimization algorithms.mq5	Script	A simple example of using population optimization algorithms without visualization
9	Test_AO_EOm.mq5	Script	EOm test stand