Crystal Structure Algorithm (CryStAl)
739
0
Contents
Introduction
I got my hands on an optimization algorithm to examine, and working on it revealed the "pitfalls" of implementing a promising idea. Meet the Crystal Structure Algorithm (CryStAl) — a metaheuristic optimization algorithm inspired by the physical process of crystal structure formation. It was proposed in 2021 by the research group of A. Kaveh and M. Kooshkebaghi and published in an article by Siamak Talatahari. CryStAl simulates the process of crystal formation, where the solutions are atoms, and the search for the optimum involves arranging these atoms into a stable structure. The algorithm uses a basis and a lattice as analogs of the spatial distribution of solutions, symmetry and regularity to form a balanced population, and energy criteria to filter out unstable solutions.
Implementation of the Algorithm
Let's examine the operating concept of the CryStAl algorithm through detailed pseudocode for all of its stages:
Input Parameters of the Algorithm
Crystal population size — the number of crystals in the population (usually 30)
Number of random crystals — the number of random crystals used to calculate the mean value of Fc (usually 3)
Problem dimensionality — the number of variables to be optimized
Maximum number of iterations — how many times the algorithm will repeat the optimization cycle
Search bounds — the lower and upper values for each variable
STAGE 1: POPULATION INITIALIZATION
Creation of Initial Crystals
For each crystal in the population and for each coordinate (variable) of that crystal:
- Formula 3 — Random initialization: New coordinate = Lower bound + Random number × (Upper bound - Lower bound), where the random number is uniformly distributed between 0 and 1
Calculate the fitness function value for this crystal and store the crystal's personal records
Initialization of global records
Find the crystal with the maximum fitness function value
Find the crystal with the minimum fitness function value
STAGE 2: MAIN OPTIMIZATION LOOP
Repeat the following steps from the first iteration up to the maximum number of iterations:
STEP 1: Calculating the mean position of all crystals (Cr_main)
For each coordinate: reset the sum for that coordinate to zero
For each crystal in the population: add the value of that coordinate for this crystal to the sum
Divide the sum by the number of crystals to obtain the mean value
Result: a vector of mean values for all coordinates has been obtained
STEP 2: Determining the best crystal (Cr_b)
Set the best-crystal index to the first crystal
For each subsequent crystal, from the second to the last:
If the fitness of this crystal is greater than that of the current best crystal: update the best-crystal index to the index of this crystal
Result: the index of the crystal with the best fitness function value is known
STEP 3: Updating the positions of all crystals
For each crystal in the population:
Stage 3.1: Saving the current state
Save the crystal's current position as the old position
Stage 3.2: Selecting a position update strategy
Randomly select one of the four position update strategies:
- Strategy 0: Simple cubic update
- Strategy 1: Cubic update with the best crystal
- Strategy 2: Cubic update with the mean of random crystals
- Strategy 3: Combined update with the best crystal and the mean position
Step 3.3: Applying the selected strategy
IF STRATEGY 0 is selected — Simple cubic update:
For each crystal coordinate: generate a random number r between 0 and 1
- Formula 4 — New coordinate = Old coordinate + r × Mean coordinate of all crystals
ELSE IF STRATEGY 1 is selected — With the best crystal:
For each crystal coordinate: generate two random numbers, r1 and r2, between 0 and 1
- Formula 5 — New coordinate = Old coordinate + r1 × Mean coordinate of all crystals + r2 × Coordinate of the best crystal
ELSE IF STRATEGY 2 is selected — With the mean of random crystals:
Substep 2.1: Calculating the mean of random crystals (Fc)
For each coordinate: reset the sum for that coordinate to zero
Repeat the specified number of times (usually 3): select a random crystal from the population (excluding the one currently being processed)
For each coordinate: add the value of this coordinate for the randomly selected crystal to the sum
For each coordinate: divide the sum by the number of selected random crystals
Result: the mean of several randomly selected crystals has been obtained
Substep 2.2: Position update
For each crystal coordinate: generate two random numbers, r1 and r2, between 0 and 1
- Formula 6 — New coordinate = Old coordinate + r1 × Mean coordinate of all crystals + r2 × Mean coordinate of random crystals
ELSE, STRATEGY 3 is selected — Combined:
Substep 3.1: Calculation of the mean of random crystals (Fc)
(Performed in exactly the same way as in strategy 2)
Substep 3.2: Position update
- For each crystal coordinate: generate three random numbers r1, r2, and r3 between 0 and 1
- Formula 7 — New coordinate = Old coordinate + r1 × Mean coordinate of all crystals + r2 × Coordinate of the best crystal + r3 × Mean coordinate of random crystals
Step 3.4: Handling boundary violations (Boundary Handling)
For each coordinate of the crystal's new position: if the coordinate is less than the lower bound OR greater than the upper bound: the crystal went out of bounds — a correction is needed
- Generate a random number between 0 and 1
- Recalculate the coordinate within the valid bounds: Coordinate = Lower bound + Random number × (Upper bound - Lower bound)
Update the crystal's position to the newly calculated position
STEP 4: Evaluation of new positions and updating records
For each crystal in the population:
Step 4.1: Calculating fitness
Calculate the fitness function value for the crystal's new position
Step 4.2: Updating the crystal's personal records
If the new fitness value is better than the personal best fitness value: update the personal best fitness value to the current value. Save the current position as the personal best position.
If the new fitness value is worse than the personal worst fitness value: update the personal worst fitness value to the current value. Save the current position as the personal worst position.
Step 4.3: Updating Global Records
If this crystal's fitness value is better than the global best fitness value: Update the global best fitness value to this crystal's fitness value. Save this crystal's position as the global best position.
If this crystal's fitness value is worse than the global worst fitness value: Update the global worst fitness value to this crystal's fitness value. Save this crystal's position as the global worst position.
STEP 5: Proceed to the next iteration
Proceed to the next iteration of the main optimization loop
End of the main optimization loop
STAGE 3: ALGORITHM TERMINATION
After completing all iterations, return the result:
- Global best position (solution to the optimization problem)
- Global best fitness value
The diagram in the figure below shows the main loop, a visualization of the four strategies, and an example of a search space with crystals.
Figure 1. Flowchart of the CryStAl Algorithm
The CryStAl algorithm strives to maintain a balance between two important aspects of optimization:
Exploration — the exploration of new regions of the search space: it is provided by strategy 0 through movement toward the mean position, enhanced by strategy 2 through information from random crystals, and helps prevent getting stuck in local optima.
Exploitation — an in-depth search near good solutions: it is implemented by strategy 1 through movement toward the best crystal, accelerates convergence to the optimum, and enables effective use of the information found.
A combined approach:
- Strategy 3 combines all three components
- Randomly selecting strategies at each iteration ensures adaptability
- Each strategy has an equal 25% probability of being selected
Role of key components:
Mean position (Cr_main) — represents the center of mass of the population, directs the search toward the central regions of the explored space, and is updated at each iteration based on the current positions of all crystals.
Best crystal (Cr_b) — contains the best solution found so far, attracts other crystals to a promising region, and enhances exploitation of the good solutions found.
Mean of random crystals (Fc) — represents local information from neighbors, adds stochasticity and diversity to the search, and helps avoid premature convergence.
Now that we have taken a detailed look at how the algorithm works, let's move on to implementing it in code. We will write a class named C_AO_CryStAl that inherits from the C_AO class and is intended to implement the Crystal Structure Algorithm.
Constructor:- Two parameters are initialized: popSize (the size of the crystal population, default 30) and numRandomCrystals (the number of random crystals used for one of the "Fc" calculations, default 3).
- These parameters are added to a structure that provides information about the algorithm's parameters.
- This method allows the values of the popSize and numRandomCrystals parameters to be updated using data from an external parameter structure.
- A check is performed to ensure that numRandomCrystals is always at least 1 and no more than popSize minus 1, in order to avoid incorrect calculations.
- Init () — a method for initializing the algorithm, which requires information about the ranges and steps of the input variables, as well as the number of epochs;
- Moving () — a method that implements the step for moving crystals in the algorithm;
- Revision () — a method responsible for revising or updating the population's state.
- numRandomCrystals — the number of random crystals used to calculate Fc.
- meanPosition — an array that stores the mean position of all crystals in the population;
- bestCrystalIdx — the index of the best crystal in the current population;
- isFirstIteration — a Boolean flag indicating whether the current iteration is the first one after initialization;
- CalculateMeanPosition () — a method for calculating the mean position of the crystals;
- CalculateFc () — a method for calculating the value of Fc, which is a criterion that depends on the selected crystals;
- SelectRandomCrystal () — a method for selecting a random crystal from the population, with the option to exclude a specified crystal.
Overall, the C_AO_CryStAl class models an optimization algorithm based on the behavior of crystals. It operates on a population of "crystals," each of which has its own position. The algorithm includes steps for moving crystals, calculating their quality metrics, and determining the best solution.
//———————————————————————————————————————————————————————————————————— class C_AO_CryStAl : public C_AO { public: //---------------------------------------------------------- ~C_AO_CryStAl () { } C_AO_CryStAl () { ao_name = "CryStAl"; ao_desc = "Crystal Structure Algorithm"; ao_link = "https://www.mql5.com/ru/articles/19899"; popSize = 30; // crystal population size numRandomCrystals = 3; // number of random crystals for Fc ArrayResize (params, 2); params [0].name = "popSize"; params [0].val = popSize; params [1].name = "numRandomCrystals"; params [1].val = numRandomCrystals; } void SetParams () { popSize = (int)params [0].val; numRandomCrystals = (int)params [1].val; if (numRandomCrystals < 1) numRandomCrystals = 1; if (numRandomCrystals > popSize - 1) numRandomCrystals = popSize - 1; } bool Init (const double &rangeMinP [], const double &rangeMaxP [], const double &rangeStepP [], const int epochsP); void Moving (); void Revision (); //------------------------------------------------------------------ int numRandomCrystals; // number of random crystals for calculating Fc private: //--------------------------------------------------------- double meanPosition []; // the mean position of all crystals (Cr_main) int bestCrystalIdx; // best crystal index bool isFirstIteration; // flag for the first iteration after initialization void CalculateMeanPosition (); void CalculateFc (int currentCrystal, double &fcArray []); int SelectRandomCrystal (int excludeCrystal); }; //————————————————————————————————————————————————————————————————————
The Init method is designed to initialize the Crystal Structure Algorithm before it begins running.
First, the StandardInit method is called, which performs basic preparatory steps common to all algorithms derived from C_AO. If this basic step fails, the Init method also fails. An array named meanPosition (mean position) is created. The size of this array is determined by the "coords" variable, which stores the number of dimensions. All elements of this array are initialized to zero. This means that, at the start of the algorithm, the mean position of all crystals is assumed to be zero.
Setting the initial values for the flags:- bestCrystalIdx (best crystal index) is set to 0. This implies that, at the start, the first crystal in the population is considered the best.
- isFirstIteration (first iteration flag) is set to 'true'. This indicates that the next iteration to be executed will be the first one following full initialization.
//———————————————————————————————————————————————————————————————————— //--- Initialization bool C_AO_CryStAl::Init (const double &rangeMinP [], const double &rangeMaxP [], const double &rangeStepP [], const int epochsP) { if (!StandardInit (rangeMinP, rangeMaxP, rangeStepP)) return false; //------------------------------------------------------------------ // Array initialization ArrayResize (meanPosition, coords); ArrayInitialize (meanPosition, 0.0); bestCrystalIdx = 0; isFirstIteration = true; return true; } //————————————————————————————————————————————————————————————————————
The CalculateMeanPosition method is designed to calculate the mean position of all "crystals" in the population. First, all elements of the meanPosition array (which stores the mean for each coordinate) are reset. This is done so that we start with a clean slate before we begin summing the current positions of the crystals.
The method then iterates over each crystal in the population. For each crystal, it also iterates through all of its coordinates (defined by "coords"). The value of each coordinate of the current crystal is added to the corresponding element of the meanPosition array. Thus, meanPosition gradually accumulates the sum of the values of each coordinate across all crystals.
Once the summation is complete, the method iterates through all the coordinates again. Each value accumulated in meanPosition is divided by the total number of crystals, popSize. This operation converts the sum of the coordinates into their mean value. As a result of executing this method, the meanPosition array will contain the mean value of each coordinate across all crystals in the population.
//———————————————————————————————————————————————————————————————————— //--- Calculating the mean position of all crystals (Cr_main) void C_AO_CryStAl::CalculateMeanPosition () { // Reset the mean values to zero ArrayInitialize (meanPosition, 0.0); // Sum the current positions of all crystals for (int i = 0; i < popSize; i++) { for (int c = 0; c < coords; c++) { meanPosition [c] += a [i].c [c]; } } // Divide by the number of crystals to obtain the mean for (int c = 0; c < coords; c++) { meanPosition [c] /= (double)popSize; } } //————————————————————————————————————————————————————————————————————
The CalculateFc method is designed to calculate the so-called "mean value of randomly selected crystals". It takes as input the index of the current crystal, currentCrystal, and an array, fcArray, into which the result will be written.
At the start, all elements of the fcArray array are set to zero. This ensures that we start the summation from scratch for each new Fc position. The method then enters a loop that repeats numRandomCrystals times. In each iteration of this loop:
- a random crystal index, randomCrystal, is selected. It is important to note that this selection is made in such a way as to exclude currentCrystal itself (i.e., the crystal for which we are currently calculating Fc);
- The position (coordinates) of this randomly selected crystal is added to the corresponding elements of the fcArray array. Thus, this array accumulates the sum of the coordinates of all the selected random crystals.
As a result, the fcArray array will contain the coordinates of the "average" crystal, obtained by averaging the positions of several randomly selected crystals from the population (excluding the current one).
//———————————————————————————————————————————————————————————————————— //--- Calculation of Fc (average value of randomly selected crystals) void C_AO_CryStAl::CalculateFc (int currentCrystal, double &fcArray []) { // Initialize the array with zeros ArrayInitialize (fcArray, 0.0); // Select numRandomCrystals random crystals (excluding the current one) for (int i = 0; i < numRandomCrystals; i++) { int randomCrystal = SelectRandomCrystal (currentCrystal); // Sum the positions of the selected crystals for (int c = 0; c < coords; c++) { fcArray [c] += a [randomCrystal].c [c]; } } // Calculate the mean value for (int c = 0; c < coords; c++) { fcArray [c] /= (double)numRandomCrystals; } } //————————————————————————————————————————————————————————————————————
The SelectRandomCrystal method is designed to select a random "crystal" from the entire population, but with one important condition: it must not select the crystal specified as excludeCrystal.
The method enters an infinite loop that will continue until the exit condition is met. Inside the loop, a random numeric index is generated, chosen from a range that covers all possible "crystals" in the population (from 0 to popSize minus 1).
After a random index is generated, a check is performed:
- if the generated index "selected" is equal to the index of the crystal we want to exclude, "excludeCrystal", then the loop continues and a new random index is generated;
- if the generated index is not equal to excludeCrystal, then the loop exit condition is satisfied.
//———————————————————————————————————————————————————————————————————— //--- Selecting a random crystal (excluding the specified one) int C_AO_CryStAl::SelectRandomCrystal (int excludeCrystal) { int selected; do selected = u.RNDminusOne (popSize); while (selected == excludeCrystal); return selected; } //————————————————————————————————————————————————————————————————————
The Moving method is the main operational loop of the algorithm and is responsible for the evolution of the population of "crystals" at each iteration. If the algorithm is run for the first time, when "revision" is 'false', it initializes the positions of all crystals in the population. For each crystal and each of its coordinates, a random value between 0 and 1 is generated. This random value is then scaled to the specified working range (rangeMin–rangeMax) for that coordinate. Discretization is used to map the position to the allowed step size rangeStep. After initialization, the "revision" flag is set to 'true' so that subsequent calls to the method perform an update rather than an initialization.
After the first iteration, the mean position of all crystals in the population is calculated first. The algorithm finds the crystal with the best fitness function value "f" in the current iteration and stores its index.
For each crystal in the population (iterating from 0 to popSize):
- the crystal's current position is copied to temporary storage, which will be used as the "old position" (Cr_old) for calculations;
- one of several available position update strategies is selected at random (there are four in total, numbered from 0 to 3);
- for two of the four strategies (strategies 2 and 3), the "mean of randomly selected crystals" (Fc) is calculated, which will also be used in the update;
- depending on the selected strategy, the new crystal position Cr_new is calculated based on its old position, the mean position of the population, the position of the best crystal, and Fc. Each strategy uses its own combination of these factors and random numbers (r, r1, r2, r3) to apply updates;
- After the new position is calculated, it is checked to ensure it falls within the specified working ranges (rangeMin–rangeMax) for each coordinate. If the new position goes out of bounds, it is adjusted by generating a new random position within the allowed bounds;
- The discretization procedure is applied to the position again using the u.SeInDiSp() function.
The Moving method iteratively moves "crystals" within the search space, using various strategies based on current information about the population (mean position, best crystal) and random factors.
//———————————————————————————————————————————————————————————————————— //--- Main algorithm step void C_AO_CryStAl::Moving () { // Initial population setup (formula 3) if (!revision) { // Random initialization of the crystal population for (int i = 0; i < popSize; i++) { for (int j = 0; j < coords; j++) { double xi = u.RNDfromCI (0.0, 1.0); a [i].c [j] = rangeMin [j] + xi * (rangeMax [j] - rangeMin [j]); a [i].c [j] = u.SeInDiSp (a [i].c [j], rangeMin [j], rangeMax [j], rangeStep [j]); } } revision = true; return; } // CRITICAL: Calculate the mean position and find the best crystal BEFORE updating CalculateMeanPosition (); // Find the best crystal for the current iteration bestCrystalIdx = 0; for (int i = 1; i < popSize; i++) { if (a [i].f > a [bestCrystalIdx].f) { bestCrystalIdx = i; } } // Main position update loop (formulas 4–7) for (int i = 0; i < popSize; i++) { // Save the current position as the old one (Cr_old) ArrayCopy (a [i].cP, a [i].c, 0, 0, coords); // Select a random strategy (0–3) int strategy = u.RNDminusOne (5); // Prepare the array for Fc (only for strategies 2 and 3) double fc []; ArrayResize (fc, coords); if (strategy == 2 || strategy == 3) { CalculateFc (i, fc); } // Apply the selected position update strategy switch (strategy) { case 0: // Strategy 1: Simple Cubic (formula 4) // Cr_new = Cr_old + r * Cr_main for (int c = 0; c < coords; c++) { double r = u.RNDfromCI (0.0, 1.0); a [i].c [c] = a [i].cP [c] + r * meanPosition [c]; } break; case 1: // Strategy 2: Cubic with the best crystals (formula 5) // Cr_new = Cr_old + r1 * Cr_main + r2 * Cr_b for (int c = 0; c < coords; c++) { double r1 = u.RNDfromCI (0.0, 1.0); double r2 = u.RNDfromCI (0.0, 1.0); a [i].c [c] = a [i].cP [c] + r1 * meanPosition [c] + r2 * a [bestCrystalIdx].c [c]; } break; case 2: // Strategy 3: Cubic update with the mean crystals (formula 6) // Cr_new = Cr_old + r1 * Cr_main + r2 * Fc for (int c = 0; c < coords; c++) { double r1 = u.RNDfromCI (0.0, 1.0); double r2 = u.RNDfromCI (0.0, 1.0); a [i].c [c] = a [i].cP [c] + r1 * meanPosition [c] + r2 * fc [c]; } break; case 3: // Strategy 4: Cubic update with the best and mean crystals (formula 7) // Cr_new = Cr_old + r1 * Cr_main + r2 * Cr_b + r3 * Fc for (int c = 0; c < coords; c++) { double r1 = u.RNDfromCI (0.0, 1.0); double r2 = u.RNDfromCI (0.0, 1.0); double r3 = u.RNDfromCI (0.0, 1.0); a [i].c [c] = a [i].cP [c] + r1 * meanPosition [c] + r2 * a [bestCrystalIdx].c [c] + r3 * fc [c]; } break; } // Checking and correcting bounds for a new position for (int c = 0; c < coords; c++) { // If the position went out of bounds, apply boundary handling if (a [i].c [c] < rangeMin [c] || a [i].c [c] > rangeMax [c]) { // Generate a new random position within the valid range double xi = u.RNDfromCI (0.0, 1.0); a [i].c [c] = rangeMin [c] + xi * (rangeMax [c] - rangeMin [c]); } // Apply discretization a [i].c [c] = u.SeInDiSp (a [i].c [c], rangeMin [c], rangeMax [c], rangeStep [c]); } } } //————————————————————————————————————————————————————————————————————
The Revision method is responsible for updating and saving the best solutions found, both at the level of each individual "crystal" (personal best solution) and at the level of the entire population (global best solution).
The method examines each "crystal" in the population to assess its current state. For each crystal, its current fitness value "f" is compared with its best-ever fitness value fB. If the current fitness value is better than the previous best, the personal best fitness value fB is updated to the current value, and the crystal position corresponding to this new best value is copied to its personal best solution cB.
Alongside updating the personal best solutions, the algorithm tracks the best solution found by the entire population so far. The crystal's current fitness value "f" is also compared with the population's overall best fitness value fB. If the crystal's current fitness value is better than the current global best, then the global best fitness value fB is updated with the current value, and the crystal's position corresponding to this new global best value is copied to the global best solution cB.
Upon completion of the Revision method, the variables fB and cB store the highest fitness values and their corresponding coordinates found over the entire run of the algorithm. Each crystal is also assigned its own best-found solution.
//———————————————————————————————————————————————————————————————————— //--- Updating and checking results void C_AO_CryStAl::Revision () { // Updating the best and worst solutions for each crystal for (int i = 0; i < popSize; i++) { // Update the personal best solution if (a [i].f > a [i].fB) { a [i].fB = a [i].f; ArrayCopy (a [i].cB, a [i].c, 0, 0, coords); } // Update the global best solution if (a [i].f > fB) { fB = a [i].f; ArrayCopy (cB, a [i].c, 0, 0, coords); } } } //————————————————————————————————————————————————————————————————————
Test Results
According to the test results, the algorithm achieves very modest performance metrics. The result itself — 26.33% — is noteworthy: it falls within the random-search score, random walk, in our ranking table. So, what does that mean? We'll come back to that later.
CryStAl|Crystal Structure Algorithm|30.0|3.0|
=============================
5 Hilly's; Func runs: 10000; result: 0.49130791883895686
25 Hilly's; Func runs: 10000; result: 0.3254668138312991
500 Hilly's; Func runs: 10000; result: 0.25762861219210875
=============================
5 Forest's; Func runs: 10000; result: 0.3770837459944563
25 Forest's; Func runs: 10000; result: 0.21650897859329796
500 Forest's; Func runs: 10000; result: 0.158890697466887
=============================
5 Megacity's; Func runs: 10000; result: 0.2938461538461538
25 Megacity's; Func runs: 10000; result: 0.1498461538461539
500 Megacity's; Func runs: 10000; result: 0.09889230769230854
=============================
All score: 2.36947 (26.33%)
These results did not leave me indifferent, so I decided to make several adjustments and improve the algorithm's convergence. Its interaction logic turned out to be not entirely effective, which made it possible to revise some illogical methods while preserving the general idea. Let's see how it turned out in the modified version.
Position Update Formulation (Strategies 0–3): the modified version uses a different update style, which looks like adding an offset to the current position a[i].c[c]. The formulas look different, although the strategies use the same concepts. Adding a certain proportion of the difference between the current position and the mean position results in a significant change in the mathematical implementation.
Additional Random Perturbation: an additional condition is introduced: "if (u.RNDprobab() < 0.1)". If this condition is met with a 10% probability, the crystal's position is updated using a[i].c[c] = u.PowerDistribution(cB[c], rangeMin[c], rangeMax[c], 20), which reduces the probability of the algorithm getting stuck. The remaining strategies are applied with a 90% probability.
Next, regarding saving the old position a [i].c P: the original version saves the old position in a [i].cP before calculating the new position. In the modified version, this position is neither saved nor used; the update is performed directly from the current position a [i].c [c].
Boundary Handling: in the original implementation, after applying the position update strategy, the method iterates through all coordinates, and if a new position went out of bounds, it generates a completely new random position within the valid range. In C_AO_CryStAlm, the update strategy is applied first, and then only discretization is performed using the u.SeInDiSp method.
//———————————————————————————————————————————————————————————————————— //--- Main step of the algorithm void C_AO_CryStAlm::Moving () { // Initial population initialization (formula 3) if (!revision) { // Random initialization of the crystal population for (int i = 0; i < popSize; i++) { for (int j = 0; j < coords; j++) { double xi = u.RNDfromCI (0.0, 1.0); a [i].c [j] = rangeMin [j] + xi * (rangeMax [j] - rangeMin [j]); a [i].c [j] = u.SeInDiSp (a [i].c [j], rangeMin [j], rangeMax [j], rangeStep [j]); } } revision = true; return; } // CRITICAL: Calculate the mean position and find the best crystal BEFORE updating CalculateMeanPosition (); // Main position update loop (formulas 4–7) for (int i = 0; i < popSize; i++) { // Select a random strategy (0–3) int strategy = u.RNDminusOne (5); // Prepare the array for Fc (only for strategies 2 and 3) double fc []; ArrayResize (fc, coords); if (strategy == 2 || strategy == 3) { CalculateFc (i, fc); } for (int c = 0; c < coords; c++) { if (u.RNDprobab () < 0.1) { a [i].c [c] = u.PowerDistribution (cB [c], rangeMin [c], rangeMax [c], 20); } else { // Apply the selected position update strategy switch (strategy) { case 0: // Strategy 1: Simple cubic (formula 4) { // Cr_new = Cr_old + r * Cr_main double r = u.RNDfromCI (0.0, 1.0); a [i].c [c] = a [i].c [c] + r * (meanPosition [c] - a [i].c [c]); } break; case 1: // Strategy 2: Cubic with the best crystals (formula 5) { // Cr_new = Cr_old + r1 * Cr_main + r2 * Cr_b double r1 = u.RNDfromCI (0.0, 1.0); double r2 = u.RNDfromCI (0.0, 1.0); a [i].c [c] = a [i].c [c] + r1 * (meanPosition [c] - a [i].c [c]) + r2 * (cB [c] - a [i].c [c]); } break; case 2: // Strategy 3: Cubic with the mean crystals (formula 6) { // Cr_new = Cr_old + r1 * Cr_main + r2 * Fc double r1 = u.RNDfromCI (0.0, 1.0); double r2 = u.RNDfromCI (0.0, 1.0); a [i].c [c] = a [i].c [c] - r1 * (meanPosition [c] - a [i].c [c]) + r2 * (fc [c] - a [i].c [c]); } break; case 3: // Strategy 4: Cubic with the best and mean crystals (formula 7) { // Cr_new = Cr_old + r1 * Cr_main + r2 * Cr_b + r3 * Fc double r1 = u.RNDfromCI (0.0, 1.0); double r2 = u.RNDfromCI (0.0, 1.0); double r3 = u.RNDfromCI (0.0, 1.0); a [i].c [c] = a [i].c [c] + r1 * (meanPosition [c] - a [i].c [c]) + r2 * (cB [c] - a [i].c [c]) + r3 * (fc [c] - a [i].c [c]); } break; } } // Apply discretization a [i].c [c] = u.SeInDiSp (a [i].c [c], rangeMin [c], rangeMax [c], rangeStep [c]); } } } //————————————————————————————————————————————————————————————————————
Initialization has not changed. There are no differences in the CalculateMeanPosition method between the versions for C_AO_CryStAl and C_AO_CryStAlm, or in the CalculateFc and SelectRandomCrystal methods. The Revision method is implemented identically in both versions, C_AO_CryStAl and C_AO_CryStAlm.
Let's take a look at the new results below: they confirm the validity of my theoretical assumptions. Now it is to find out why the original version turns a directed search into a random walk with drift.
CryStAlm|Crystal Structure Algorithm M|30.0|3.0|
=============================
5 Hilly's; Func runs: 10000; result: 0.8065665859751701
25 Hilly's; Func runs: 10000; result: 0.5614723387106613
500 Hilly's; Func runs: 10000; result: 0.2909726108642525
=============================
5 Forest's; Func runs: 10000; result: 0.7522883470073202
25 Forest's; Func runs: 10000; result: 0.48674341114821484
500 Forest's; Func runs: 10000; result: 0.2101468006131327
=============================
5 Megacity's; Func runs: 10000; result: 0.48461538461538467
25 Megacity's; Func runs: 10000; result: 0.2701538461538461
500 Megacity's; Func runs: 10000; result: 0.11110769230769342
=============================
All score: 3.97407 (44.16%)
When we compare the two versions of the Crystal Structure Algorithm, it may seem at first glance that the changes are minimal. Both methods initialize the population in the same way, compute the mean position, and use four position update strategies. However, behind this superficial similarity lie fundamental differences that lead to an almost twofold difference in performance. Let's break down each change to understand why the modified version is so much more effective.
A fundamental error: absolute addition instead of directed movement. To understand the key problem with the original algorithm, let’s consider a one-dimensional optimization problem: finding the maximum of a function in the range from 0 to 100. The mean position of the population (the point toward which the crystals are attracted) is at point 60. We have three crystals in different positions: the first at point 10, the second at point 60, and the third at point 90. The random coefficient r is the same for all of them: 0.5.
Original formula: x_new = x_old + r × mean
- Crystal 1: 10 + 0.5 × 60 = 10 + 30 = 40 (shift of +30)
- Crystal 2: 60 + 0.5 × 60 = 60 + 30 = 90 (shift of +30)
- Crystal 3: 90 + 0.5 × 60 = 90 + 30 = 120 (shift of +30, OUT OF BOUNDS!)
Did you notice the problem? ALL THREE CRYSTALS MOVED BY +30, no matter where they were! The first crystal was far to the left of the target and moved to the right — that is fine. The second one was already at the target and was pushed to the right — that is absurd. The third one was to the right of the target; it was pushed even farther to the right and went out of bounds — a disaster.
This is not "movement toward the target." This is an identical push applied to everyone in one direction, as if a strong wind had blown and shifted everyone by the same distance. Moreover, the direction and strength of this "wind" are determined by the absolute coordinate values of the target, not by the relative positions of the crystals.
Modified formula: x_new = x_old + r × (mean - x_old)
- Crystal 1: 10 + 0.5 × (60 - 10) = 10 + 25 = 35 (moved 50% of the way toward the target)
- Crystal 2: 60 + 0.5 × (60 - 60) = 60 + 0 = 60 (already at the target, stayed in place)
- Crystal 3: 90 + 0.5 × (60 - 90) = 90 - 15 = 75 (traveled 50% of the way BACK to the target)
Do you see the difference? Each crystal moves toward its target, in proportion to its distance from it. Those who are far away take a big step. Those who are close take a small one. Those who are already there stand still. Whoever is beyond the target moves back. This is exactly what is needed for optimization.
Problem 1: Blow-up with large coordinates. If the variables being optimized have a range, say, from 0 to 1000, and the mean position of the population happens to be around 600–700, then a number on the order of 600–700 (multiplied by a random coefficient from 0 to 1) is added to EVERY crystal. Even a crystal located at position 650 (almost exactly at the target!) will get a displacement of +300...+600 and jump to the 1000 boundary or beyond. The algorithm turns into a random number generator at the boundaries of the search space.
Problem 2: Dependence on the coordinate scale. If you have a problem with variables of different scales (for example, x ranging from 0 to 10 and y ranging from 0 to 1000), then the displacements along the y-axis will be 100 times greater than those along the x-axis, simply because of the difference in scale! This has nothing to do with how far the crystal is from the optimum along each coordinate.
Problem 3: All crystals move in the same direction. When you add the same vector to all the crystals (albeit with different random coefficients), the entire population shifts as a single unit in one direction. This is not “searching for the optimum”; it is “population drift.” Imagine a flock of birds that is supposed to be looking for food, but instead all the birds fly in unison toward the northeast, with slight random deviations. This is not search; it is migration.
Problem 4: Constant reinitialization. Because of the enormous displacements, the crystals are constantly going out of bounds. The algorithm has to reinitialize them at random locations. This means that **any accumulated information about the direction of movement is lost**. The crystal was moving toward a promising area, went out of bounds, and was teleported to a random location — its entire movement history was erased. Essentially, the algorithm reverts to a random search every few iterations.
Memory of the best solution: local vs. global. The second important change concerns which "best solution" is used in the formulas. The modified version keeps a global record book. When someone reaches an altitude of 3,000 meters, it is recorded and never forgotten. Even if all the explorers later went back down, the algorithm remembers: "The best spot we ever found was right there, at an altitude of 3,000." And it is precisely this point that the other explorers are attracted to. This is particularly important in optimization problems, where the population may temporarily deteriorate. For example, you applied a mutation, and all the solutions became worse. In the original version, the algorithm will start following inferior solutions. In the modified version, it will continue to stay on course toward the truly best solution found over the entire run.
Adding intensive local exploitation. The third key change concerns a new mechanism that is not present in the original. In the modified version, each crystal coordinate is generated with a 10% probability in a special way around the global best solution, using the "PowerDistribution" function with a parameter value of 20. The original algorithm works as a pure metaheuristic: the crystals move throughout the space according to their strategies. The modified version adds a local search component: every tenth coordinate jumps into the neighborhood of the best solution found. A value of 20 in PowerDistribution indicates a very high concentration — roughly as if you were looking for your keys within a half-meter radius of where you remember last seeing them. This is not a random point in this region, but a point that is very likely to be very close to the center.
Simplified boundary handling. The fourth change concerns what happens when a crystal attempts to go out of bounds in the search space. In the modified version, the coordinate is simply passed to the u.SeInDiSp function, which clips the value to the bounds and applies discretization. Why was this change possible? Because, thanks to the correct movement formulas, the crystals in the modified version very rarely go out of bounds
In the original version, because of absolute vector addition, the crystals constantly go out of bounds, making reinitialization necessary. But this reinitialization comes at a high cost: first, it requires the generation of random numbers and additional computations; second, it turns directed motion into random movement. The crystal might have been moving in a promising direction, but it went out of bounds and ended up in a completely random place, losing its entire movement trajectory.
Memory usage optimization. The fifth change may seem technically minor, but it demonstrates an understanding of what the algorithm truly needs to function. In the original version, before updating the crystal's position, its current position is copied into a special array called cP (previous position), and this copy is then used in the formulas. For example, the new coordinate is calculated as the old coordinate from the cP array plus something else. The modified version does not include this copying. The formula simply takes the current value of the coordinate, calculates the change, and stores the result back in the same variable.
Why did the original need a copy? Because the formula looked like "x_new = x_old + something." To calculate the right-hand side, x_old had to be saved before making any changes. In the modified version, the formula has been rearranged so that it can work directly with the current value: "x = x + change." For a population of 30 crystals in 10-dimensional space, this saves 300 data copy operations on each iteration. It may seem like a minor detail, but when an algorithm runs for thousands of iterations, these small details add up.
In the visualization, we can observe how the modified version forms “crystalline shapes” during operation; the spread of values for functions of low and medium dimensionality is also noticeable.
CryStAlm on the Hilly test function
CryStAlm on the Forest test function
CryStAlm on the Forest test function
CryStAlm on the standard Ackley test function
In the ranking table, the CryStAlm algorithm is shown for reference; it fell just short of being assigned a place.
| 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) | 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 | DOAdingom | dingo_optimization_algorithm_M | 0.47968 | 0.45367 | 0.46369 | 1.39704 | 0.94145 | 0.87909 | 0.91454 | 2.73508 | 0.78615 | 0.86061 | 0.84805 | 2.49481 | 6.627 | 73.63 |
| 2 | 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 |
| 3 | 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 |
| 4 | 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 |
| 5 | (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 |
| 6 | 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 |
| 7 | 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 |
| 8 | 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 |
| 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 | 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 |
| 26 | 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 |
| 27 | 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 |
| 28 | 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 |
| 29 | 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 |
| 30 | 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 |
| 31 | 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 |
| 32 | 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 |
| 33 | 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 |
| 34 | (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 |
| 35 | 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 |
| 36 | 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 |
| 37 | 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 |
| 38 | 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 |
| 39 | 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 |
| 40 | 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 |
| 41 | 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 |
| 42 | 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 |
| 43 | 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 |
| 44 | DA_duelist | duelist_algorithm | 0.92782 | 0.53778 | 0.27792 | 1.74352 | 0.86957 | 0.47536 | 0.18193 | 1.52686 | 0.62153 | 0.33569 | 0.11715 | 1.07437 | 4.345 | 48.28 |
| 45 | 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 |
| CryStAlm | crystal_structure_algorithm_M | 0.80656 | 0.56147 | 0.29097 | 1.65900 | 0.75228 | 0.48674 | 0.21014 | 1.44916 | 0.48461 | 0.27015 | 0.11111 | 0.86587 | 3.974 | 44.16 | |
| 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
A beautiful idea, ruined by its implementation. The Crystal Structure Algorithm was introduced in 2021 with an appealing concept: using the principles of crystal structure formation for global optimization. The authors proposed a population-based algorithm in which "crystals" move under the influence of the population's mean position, the best solution found, and randomly selected neighbors. Four movement strategies, an elegant mathematical model. On paper, it all looked convincing.
The problem is that there is a critical stage between the conceptual idea and a working implementation — the translation of a physical metaphor into specific mathematical operations. And this is exactly where the original CryStAl fails. Formulas that look reasonable in the paper turn guided search into chaotic wandering in practice. This is not just a story about "finding a bug and fixing it"; it is a fundamental lesson that, in the development of metaheuristics, the concept and the mathematical implementation must be in harmony; otherwise, the algorithm will not work as intended, or it will not work at all.
The gap between intent and execution. Let's look at what the authors WANTED to implement in the original CryStAl. The idea is clear: crystals should be attracted to three key points — the mean position of the entire population (center of mass, global information), the best crystal found so far (exploitation), and the mean of random neighbors (local information, diversity). The four strategies combine these attractions in different ways. In theory, this creates a balance between exploring new areas of the search space and deepening the search in promising regions. It is a beautiful physical analogy: atoms interact with their neighbors and strive to achieve an equilibrium configuration. The authors WANTED to implement "attraction toward the target," but WROTE "vector addition." These are different operations!
A systemic error in thinking: formulas instead of geometry. The root of the problem runs deeper than a single incorrect formula. This is a cognitive error in algorithm design: focusing on formulas instead of understanding what those formulas do in the geometry of the search space. The authors of the original CryStAl wrote: "the crystal moves toward the mean position, formula: x + r × mean." They focused on the algebraic expression without visualizing what would happen to the crystals in different parts of the search space. The right approach requires asking the following questions:
- What should happen to a crystal that is already at the target point? Answer: it should stay there or take a very small step. The original formula flings it away from the target. The modified version leaves it in place (when x_old = mean, we get x + 0 × (mean - x) = x).
- What should happen to a crystal that is beyond the target? Answer: it should move back toward the target. The original formula pushes it farther in the same direction. The modified version reverses direction (a negative difference results in backward motion).
- Should the step size depend on the absolute values of the coordinates? Answer: no, it should depend on the distance to the target. The original depends on absolute values. The modified version depends on distances.
- What happens when variables have different scales? Original: the displacements along different coordinates are not proportional to their relative distances to the target. Modified version: displacements are automatically scaled according to relative distances.
The answers to these questions immediately show that the formula x + r × target is geometrically incorrect for modeling attraction. By contrast, the formula x + r × (target - x) is geometrically correct: it is interpolation — moving a fraction r of the way toward the target — with the natural property of decelerating as the point approaches the target.
No amount of formulas can save a flawed foundation. It may seem as though CryStAl’s four strategies offset each other’s shortcomings. Strategy 0 with a single target point does not work very well, but should not Strategy 3 with three target points be better? In practice, the opposite happens: stacking flawed operations only makes the problem worse.
Strategy 3 formula: x_new = x_old + r1 × mean + r2 × best + r3 × fc, where fc is the mean of random crystals. Now we add THREE large numbers (the absolute coordinates of the three target points). Let's take an example: variables range from 0 to 1000, the crystal is at [640, 430, 770], mean = [650, 420, 780], best = [645, 425, 775], fc = [655, 415, 785], and the coefficients are r1 = r2 = r3 = 0.4.
First coordinate: 640 + 0.4 × 650 + 0.4 × 645 + 0.4 × 655 = 640 + 260 + 258 + 262 = 1420 — even further out of bounds!
Strategy 3 was supposed to be the most complex and effective one (balancing global information, local information, and exploitation), but in practice it generates the largest outliers because it sums even more absolute values. The problem isn't the number of components, but rather that the basic operation (absolute addition) is fundamentally incorrect.
From random search to directed optimization. What does the original CryStAl actually do in practice? After several iterations, the following mode of operation is established: at each step, the population receives huge displacements, most crystals go out of bounds, and their coordinates are randomly reinitialized. In fact, every 2–3 iterations, the algorithm generates a new, nearly random population, retaining only a small fraction of the coordinates from the previous generation.
The modified CryStAlm works in a fundamentally different way. The crystals actually move toward their target points rather than moving randomly. The trajectories are smooth and controlled. Boundary violations are rare. The population methodically converges toward promising regions while maintaining diversity through random strategies and a repulsion mechanism. Each iteration accumulates information about the function landscape and uses it for subsequent steps, without erasing the previous state through random regeneration.
A methodological lesson for researchers. The story of CryStAl is not just an example of a failed implementation. This is a systemic problem in the development of metaheuristic algorithms.
A broader context: an epidemic of weak algorithms. CryStAl is not unique. In the literature on metaheuristics over the past 15 years, there have been hundreds of “new” algorithms inspired by various metaphors, and many of them suffer from the same problem: a beautiful idea, a non-rigorous implementation, testing on a limited set of problems, and publication with claims of competitiveness. Then other researchers use these algorithms as baselines, compare their own work against them, and create a chain of comparisons with fundamentally weak foundations.
Modifying CryStAl into CryStAlm gives us material for an in-depth analysis. Minimal changes to the code (rewrite the motion formulas, add PowerDistribution, use cB instead of bestCrystalIdx), but fundamental changes in behavior.
A comparison of CryStAl and CryStAlm offers several critical lessons for the research community in the field of metaheuristic optimization:
Lesson 1: A concept is not the same as its implementation. A beautiful algorithm idea is just the beginning. Translating an idea into mathematical operations requires a deep understanding of the geometry of the search space. Formulas should do what you want them to do, not just look plausible on paper. Check not only "what the formula calculates," but also "what it does with agents in different situations."
Lesson 2: Abstraction hides errors. When you write “the crystal is attracted to the mean position,” it is easy to overlook the fact that the formula implements a parallel shift rather than attraction. Abstract thinking ("moving toward a target") must be accompanied by concrete verification ("exactly how the coordinates will change"). Always check formulas using specific numerical examples.
Lesson 3: Standard tests are not enough. Testing on a CEC benchmark with normalized variables may mask critical issues that become apparent at realistic scales or in mixed ranges. Include tests with "tricky" conditions: very large ranges, variables with different scales, and narrow optima. Visualize the algorithm's behavior, not just the results.
Lesson 4: A multitude of mechanisms cannot compensate for a flawed foundation. If the basic operation (for example, the movement formula) is incorrect, adding complex layers on top will not help. First, make sure the foundation is solid, then add improvements. One good mechanism is better than ten bad ones.
Lesson 5: Comparisons require checking the baseline. If your new algorithm outperforms the existing one, ask yourself: does that baseline even work? Quick check: compare the baseline against random search. If the difference is small, you are comparing against a dummy baseline, and your result is meaningless.
Lesson 6: Implementation details are critical. The difference between an “almost correct” formula and a “correct” one can be enormous in terms of performance. x + r × mean and x + r × (mean - x) differ by just a few characters in the code, but the quality of the results differs by a factor of two. The details matter. Do not dismiss them as "minor technicalities."
CryStAl demonstrates that metaheuristics is not merely a collection of metaphors and formulas. This is an engineering discipline in which understanding how operations affect the dynamics of search within the solution space is just as important as the original concept. A beautiful idea, implemented carelessly, results in an algorithm that behaves like a random search. The same idea, implemented with an understanding of geometry and dynamics, results in a working tool.
To move the field forward, it is not enough to simply come up with new metaphors. We need to create algorithms that actually do what they claim to do — explore the space in a directed manner, accumulate information, and converge to optima. And this requires not just writing down formulas, but deeply analyzing whether those formulas correctly implement the intended concept.
Figure 2. Color gradation of algorithms for the corresponding tests
Figure 3. Histogram of algorithm test results (on a scale from 0 to 100; the higher the value, the better, where 100 is the maximum possible theoretical result; the archive contains a script for calculating the rating table)
Pros and cons of the CryStAlm algorithm:
Pros:
- A relatively simple implementation.
- Fast.
- Consistently average results (yellow, with no major failures on individual tests).
Cons:
- Spread of values for functions of low and medium dimensionality.
An archive containing the latest versions of the algorithm code is attached to the article. The author of this article is not responsible for the absolute accuracy of the descriptions of the canonical algorithms; changes have been made to many of them to improve search capabilities. The conclusions and findings presented in the articles are based on the results of the experiments conducted.
Programs used in the article
| # | Name | Type | Description |
|---|---|---|---|
| 1 | #C_AO.mqh | include file | Parent class for population-based optimization algorithms |
| 2 | #C_AO_enum.mqh | include file | Enumeration of population-based optimization algorithms |
| 3 | TestFunctions.mqh | include file | Test function library |
| 4 | TestStandFunctions.mqh | include file | Test bench function library |
| 5 | Utilities.mqh | include file | Utility function library |
| 6 | CalculationTestResults.mqh | include file | Script for calculating results for a 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_AO_CryStAlm.mq5 | Script | Test bench for CryStAlm |
Translated from Russian by MetaQuotes Ltd.
Original article: https://www.mql5.com/ru/articles/19899
Warning: All rights to these materials are reserved by MetaQuotes Ltd. Copying or reprinting of these materials in whole or in part is prohibited.
This article was written by a user of the site and reflects their personal views. MetaQuotes Ltd is not responsible for the accuracy of the information presented, nor for any consequences resulting from the use of the solutions, strategies or recommendations described.
Implementing and Benchmarking Bag-of-SFA-Symbols (BOSS) Against Dynamic Time Warping (DTW)
Building a Divergence System (Part III): The Adaptive SuperTrend EA
Building a Volume-Based Liquidity Heatmap Indicator in MQL5
Streaming MetaTrader 5 Trade Events to a Local HTTP Server Using WinINet in MQL5
- Free trading apps
- Over 8,000 signals for copying
- Economic news for exploring financial markets
You agree to website policy and terms of use