The systematic benchmarking of metaheuristic algorithms on continuous optimization problems has been a central activity of the computational intelligence community since the early 2000s. The IEEE CEC has played a pivotal role in this effort by publishing standardized benchmark suites at regular intervals. Liang et al. [14] introduced the CEC 2013 benchmark comprising 28 functions spanning unimodal, multimodal, and composition categories, establishing a widely adopted evaluation framework. The CEC 2014 benchmark, curated by Liang et al. [15], extended this to 30 functions with shifted and rotated variants to prevent algorithms from exploiting structural regularities. Awad et al. [16] further expanded the suite for CEC 2017 with 29 functions, including hybrid compositions of increasing complexity.

Among the algorithms evaluated on these benchmarks, DE has consistently demonstrated strong performance. Das and Suganthan [17] provided a comprehensive survey of DE variants, documenting the algorithm's success across hundreds of benchmark experiments. Tanabe and Fukunaga [18] proposed Success-History-based Adaptive Differential Evolution (SHADE), which adapts the scaling factor F and crossover rate CR using a historical memory mechanism, achieving top rankings in the CEC 2013 competition. More recently, Mirjalili and Lewis [19] introduced the Whale Optimization Algorithm (WOA), inspired by the bubble-net feeding behavior of humpback whales, and demonstrated competitive performance on 29 benchmark functions. Mohamed et al. [20] proposed LSHADE-cnEpSin, an enhanced adaptive DE variant that won the CEC 2017 competition, further cementing DE's position as a leading continuous optimizer.

PSO has also been the subject of extensive benchmarking. Clerc and Kennedy [21] introduced the constriction coefficient approach to PSO, providing theoretical convergence guarantees under simplified assumptions. Liang et al. [22] proposed Comprehensive Learning Particle Swarm Optimization (CLPSO), which diversifies the personal best information used to update particle velocities, achieving significant improvements on multimodal functions. The GWO algorithm of Mirjalili et al. [11] has attracted considerable attention since its introduction; comparative studies by Faris et al. [23] documented GWO's competitive performance on CEC benchmarks, particularly on unimodal and simple multimodal functions, while noting limitations on highly deceptive composition functions.

The application of metaheuristic algorithms to the solution of systems of nonlinear equations represents a natural and increasingly productive line of research. Grosan and Abraham [24] provided one of the earliest systematic studies, demonstrating that reformulating a system of n equations in n unknowns as an unconstrained optimization problem minimizing f(x) = Σi = 1n [fi(x)] enables metaheuristic algorithms to locate multiple roots simultaneously, an advantage over Newton-type methods that converge to a single root per invocation. Mo et al. [25] applied a conjugate direction PSO to systems arising in chemical engineering, achieving convergence on systems with up to 20 variables where the Levenberg-Marquardt method failed due to poor initial guesses.

Pourjafari and Mojallali [26] proposed a modified GA with adaptive crossover and mutation rates for solving nonlinear systems in interval arithmetic, successfully locating all roots within specified intervals. Hirsch et al. [27] applied continuous Greedy Randomized Adaptive Search Procedure (GRASP) to nonlinear systems with up to 50 variables, demonstrating that multi-start strategies can achieve high success rates even on systems with dense Jacobian structures. More recently, Nadimi-Shahraki et al. [28] applied an improved WOA to power flow equations with 30 buses, reporting a 96.4% success rate and mean residuals below 10⁻¹².

Parameter estimation in dynamical systems governed by ordinary and PDE is a classical inverse problem with applications in physics, biology, chemistry, and climate science. The problem is inherently challenging because the forward model (the numerical solution of the ODE/PDE) must be evaluated at each optimizer iteration, and the resulting objective landscape is typically non-convex and multimodal. Babaei [29] applied a hybrid GA-Nelder-Mead approach to estimate thermal diffusivity in the one-dimensional heat equation, achieving relative errors below 0.5% on synthetic data with 5% added noise. Several authors have explored Runge-Kutta-based optimization methods; notably, Ahmadianfar et al. [30] proposed the Runge-Kutta Optimizer (RUN) algorithm inspired by the slope-arithmetic structure of Runge-Kutta integration, demonstrating state-of-the-art performance on CEC 2014 benchmarks and engineering design problems.

For biological models, Mendes and Kell [31] showed that evolutionary strategies outperform gradient-based methods for estimating parameters in the Lotka-Volterra predator-prey system, particularly when data are sparse or noisy. Rodriguez-Fernandez et al. [32] compared scatter search, DE, and Levenberg-Marquardt for parameter estimation in biochemical reaction networks, concluding that DE achieved the best balance between accuracy and computational cost. In the context of epidemiological models, Godio et al. [33] applied PSO to estimate SIR model parameters for COVID-19 in Italy, obtaining excellent fits to early-phase incidence data with R² values exceeding 0.998.

Combinatorial optimization problems, characterized by discrete, finite (but exponentially large) search spaces, pose a qualitatively different challenge for metaheuristic algorithms. The TSP, a canonical NP-hard problem, has served as the primary benchmark since the inception of the field. Dorigo and Gambardella [34] introduced Ant Colony Optimization (ACO) and demonstrated near-optimal performance on TSP instances with up to 200 cities. Helsgaun [35] proposed the Lin-Kernighan-Helsgaun (LKH) algorithm, a sophisticated local search heuristic that has produced optimal or near-optimal solutions for instances with over 10,000 cities. For GA-based approaches to TSP, Larranaga et al. [36] provided a comprehensive survey of representation schemes and crossover operators, establishing the Ordered Crossover (OX) and Partially Matched Crossover (PMX) as standard operators for permutation-encoded problems.

Graph coloring, another NP-hard problem, has been addressed by metaheuristic algorithms with increasing success. Malaguti et al. [37] proposed a GA with greedy partition crossover for the graph coloring problem, achieving competitive results on benchmark instances from the Discrete Mathematics and Theoretical Computer Science Center (DIMACS). Lu and Hao [38] introduced a memetic algorithm combining tabu search with crossover, producing state-of-the-art results on several DIMACS instances, including the challenging flat300 and le450 families. More recently, hybrid approaches combining metaheuristics with exact methods have shown promise; Moalic and Gondran [39] proposed a population-based approach with intensive local search that established new best-known solutions on several DIMACS instances.

Six metaheuristic algorithms were selected for this comparative study, chosen to represent the major algorithmic paradigms in the field: evolutionary computation (GA, DE), swarm intelligence (PSO, SSA), and physics- and mathematically-inspired optimization (GWO, SCA). All algorithms were implemented in MATLAB R2023b and executed on an Intel Core i9-13900K processor (24 cores, 5.8 GHz) with 64 GB DDR5 RAM under Ubuntu 22.04 LTS. The population size for all algorithms was fixed at N = 100 to ensure fair comparison, and the maximum number of function evaluations was set to 500,000 for all continuous optimization experiments. Each experiment was replicated 51 times with different random seeds to ensure statistical robustness.

GA: The GA implementation follows the canonical framework of Deb et al. [40]. Binary tournament selection is employed to select parent individuals for reproduction. The Simulated Binary Crossover (SBX) operator with distribution index ηc = 20 generates offspring, while polynomial mutation with distribution index ηm = 20 introduces perturbations. The crossover probability is set to pc = 0.9 and the mutation probability to pm = 1/D, where D is the problem dimension. Elitism is implemented by preserving the best individual across generations. For combinatorial problems such as TSP, the OX operator and the swap mutation are used with permutation encoding.

PSO: the PSO variant used in this study follows the recommendations of Shi and Eberhart [41] with inertia weight linearly decreasing from w = 0.9 to w = 0.4 over the course of the run. Cognitive and social acceleration coefficients are set to c1 = c2 = 2.0. Velocity clamping is applied with Vmax = 0.2 × (xmax − xmin) to prevent excessive velocity magnitudes that could cause particle divergence. The global-best topology is employed, in which all particles share the swarm's overall best position.

DE: the DE implementation uses the DE/rand/1/bin strategy, as recommended by Storn and Price [10] for general-purpose optimization. The scaling factor is set to F = 0.8 and the crossover rate to CR = 0.9. For each target vector xi, three distinct individuals xr1, xr2, and xr3 are randomly selected from the population, and the mutant vector is constructed as vi = xr1 + F · (xr2 − xr3). Binomial crossover then generates the trial vector ui, which replaces the target if it yields a superior fitness value.

GWO: the GWO algorithm of Mirjalili et al. [11] models the social hierarchy of grey wolves (alpha, beta, delta, omega) and their encircling-hunting behavior. The parameter a is linearly decreased from 2 to 0 over the course of iterations, controlling the transition from exploration to exploitation. The positions of the three best wolves (α, β, δ) guide the movement of the remaining search agents. No additional parameters require tuning beyond the population size and maximum iterations.

SSA: the SSA of Mirjalili et al. [12] distinguishes between leader salps (at the front of the chain) and follower salps. The key parameter c1 is computed as c1 = 2e − (4l/L)², where l is the current iteration, and L is the maximum number of iterations. This parameter controls the balance between exploration and exploitation, with large values of c1 encouraging exploration in early iterations.

SCA: the SCA of Mirjalili [13] updates search agent positions using sine and cosine functions. The parameter a is set to 2, and r1 is linearly decreased from a to 0 over iterations, determining whether the sine or cosine function is used and the magnitude of the position update.

The CEC 2014 benchmark suite [15] was adopted as the primary testbed for continuous function optimization. This suite comprises 30 functions organized into four categories: unimodal functions (F1–F3), simple multimodal functions (F4–F16), hybrid functions (F17–F22), and composition functions (F23–F30). For this study, 23 representative functions were selected spanning all four categories: F1–F3 (unimodal), F4– F16 (multimodal), and F17–F23 (composition). Each function is defined over the domain [−100, 100]D with a global minimum value shifted to f* = i × 100 for function Fi. All functions are shifted and rotated to prevent algorithms from exploiting coordinate-axis alignment.

Four dimensions were tested: D = 10, 30, 50, and 100. Each algorithm was run 51 times per function per dimension, yielding a total of 51 × 23 × 4 × 6 = 28,152 independent optimization runs for this problem set alone. The error metric f(x) − f ∗ was recorded at the termination of each run.

Eight systems of nonlinear equations were drawn from the mathematical and engineering literature, spanning dimensionalities from n = 2 to n = 100. Each system was reformulated as an unconstrained minimization problem by defining the objective function as f(x) = Σi = 1n [gi(x)]², where gi(x) = 0 represents the i-th equation. A solution was deemed to have been found when f(x) < 10 − 10. Each algorithm was allowed 500,000 function evaluations per run, with 51 independent runs per system.

Five parameter estimation problems involving ordinary and PDE were formulated. In each case, synthetic "observed" data were generated by solving the ODE/PDE with known true parameter values using a high-accuracy numerical solver (MATLAB's ode45 with absolute and relative tolerances of 10−12 for ODEs; a Crank-Nicolson finite difference scheme with Δx = 0.01 and Δt = 0.001 for PDEs). Gaussian noise with zero mean and standard deviation σ = 0.01 × max|y(t)| was added to simulate experimental measurement error. The objective function for each problem is the Mean Squared Error (MSE) between the observed data and the model output evaluated at the candidate parameter vector:

MSE(θ) = (1/N) Σj = 1N [yobs(tj) − ymodel(tj; θ)]². (1)

The five systems are: 1) Lotka-Volterra predator-prey model (dx/dt = αx − βxy, dy/dt = δxy − γy; four parameters: α, β, δ, γ), 2) SIR epidemic model (dS/dt = −βSI/N, dI/dt = βSI/N − γI, dR/dt = γI; three parameters: β, γ, I0), 3) Van der Pol oscillator (x″ − μ(1 − x²)x′ + x = 0; two parameters: μ, x(0)), 4) one-dimensional heat equation (∂u/ ∂t = κ ∂²u/ ∂x²; three parameters: thermal diffusivity κ, and two boundary condition coefficients), and 5) Burgers' equation (∂u/ ∂t + u ∂u/ ∂x = ν ∂²u/ ∂x²; two parameters: viscosity ν and initial amplitude A).

Ten TSP instances were selected from the Traveling Salesman Problem Library (TSPLIB) repository [46], spanning sizes from 51 to 1,002 cities. For the TSP, each algorithm uses permutation encoding with problem-specific operators: OX and swap mutation for GA; random-key encoding with velocity-based reordering for PSO; random-key encoding with differential mutation for DE; and nearest-position rule-based updates for GWO, SSA, and SCA. Each algorithm was allotted 2,000,000 function evaluations for TSP instances, with 30 independent runs. Additionally, five graph-coloring instances from the DIMACS benchmark repository were selected, ranging in size from 25 to 450 vertices. The objective for graph coloring is to minimize the number of colors k such that no two adjacent vertices share a color.

Given the stochastic nature of metaheuristic algorithms, rigorous statistical analysis is essential for drawing reliable conclusions from experimental results. The statistical testing protocol adopted in this study follows the recommendations of Derrac et al. [47] for comparing multiple algorithms over multiple problems. The Friedman test [48] is a nonparametric statistical test for detecting differences among treatments (algorithms) across multiple blocks (problems). Under the null hypothesis H0, all algorithms perform equally, yielding identical mean ranks. The Friedman statistic χ²F follows a chi-squared distribution with k − 1 degrees of freedom, where k is the number of algorithms. The significance level is set to α = 0.05 throughout. When the Friedman test rejects H0, the Nemenyi post hoc test [49] is used for pairwise comparisons. Two algorithms are considered significantly different if their mean ranks differ by more than the Critical Difference (CD), CD = qα √(k(k + 1) / (6N)), where qα is the studentized range statistic, and N is the number of problems.

Additionally, Wilcoxon signed-rank tests with the Bonferroni correction are used for pairwise comparisons of algorithms across individual problem domains. The Bonferroni correction adjusts the significance level by dividing α by the number of pairwise comparisons (k(k−1)/2 = 15 for six algorithms), yielding an adjusted significance level of αadj = 0.05/15 ≈ = 0.0033.

Table 5 presents the mean error and standard deviation for each algorithm across all 23 CEC 2014 benchmark functions at dimension D = 30, over 51 independent runs. The best result for each function is highlighted. Functions F1–F3 are unimodal, F4–F16 are multimodal, and F17–F23 are composition functions.

The results in Table 5 reveal a clear hierarchy among the six algorithms on continuous benchmark functions at D = 30. DE achieves the best mean error on 19 of 23 functions, dominating all three unimodal functions (F1–F3) with errors approaching machine precision (order 10−22 to 10−28). On multimodal functions (F4– F16), DE maintains its lead, though the margins are narrower on functions with complex landscape features (F13, F14, F16). GWO emerges as the clear second-best algorithm, with errors typically 3–7 orders of magnitude above DE but substantially below the remaining algorithms. PSO ranks third overall, showing particular strength on unimodal functions where its global-best topology provides efficient convergence. The GA performs worst among the population-based methods on continuous functions, likely due to the disruptive nature of SBX crossover in high-dimensional continuous spaces. SCA consistently produces the worst results, suggesting that its sinusoidal position update mechanism lacks the adaptive precision needed for fine-grained convergence on shifted and rotated functions.

Composition functions (F17–F23) present a more nuanced picture. On F17 and F20, GWO outperforms DE, suggesting that the hierarchical encircling mechanism of GWO provides advantages on landscapes with multiple basins of attraction connected by saddle points. On F23, which features a highly deceptive multimodal landscape, all algorithms converge to similar mean errors, indicating that none can reliably escape the dominant local attractor.

Table 6 summarizes the scalability analysis by reporting the mean Friedman rank of each algorithm across all 23 functions at each dimensionality.

Several notable trends emerge from the scalability analysis: 1) DE's Friedman rank improves monotonically with dimensionality, from 1.52 at D = 10 to 1.35 at D = 100, confirming the well-documented scalability advantage of differential mutation in high-dimensional spaces, 2) GWO's rank degrades from 1.87 to 2.57 as dimensionality increases, suggesting that the three-leader encircling mechanism becomes less effective when the search space volume grows exponentially, and 3) SSA exhibits a non-monotonic pattern, improving from rank 3.78 at D = 10 to 3.52 at D = 100 (but with a peak degradation at D = 30), possibly reflecting the adaptive nature of the c1 parameter. The Friedman test rejects the null hypothesis of equal performance at all four dimensionalities (p < 0.001), confirming statistically significant differences among the algorithms. Convergence behavior analysis reveals further insights. On the Sphere function (F1) at D = 30, DE converges monotonically with the steepest descent rate, reaching error levels of 10−20 within approximately 200,000 function evaluations. GWO shows similarly rapid initial convergence but plateaus at error levels of 10−18– 10−21. PSO converges quickly to moderate accuracy (10−4) but slows considerably thereafter due to stagnation around the global best. On the Rastrigin function (F9), which features approximately 1030 local optima at D = 30, DE's convergence curve shows characteristic "staircase" behavior corresponding to successive escapes from local basins. GWO displays smoother convergence on F9, suggesting that the social hierarchy mechanism naturally preserves diversity.

Table 7 presents the performance of all six algorithms on the eight nonlinear systems, reporting the success rate (the percentage of runs achieving f(x∗) < 10−10), the mean objective value at termination, the mean iteration count, the mean computation time, and the number of distinct solutions found.

GWO demonstrates the strongest overall performance on nonlinear equation solving, achieving the highest success rate on all eight systems and the lowest mean objective values. The overall success rate of GWO across all systems is 97.8%, substantially exceeding DE (93.4%), PSO (88.0%), GA (81.7%), SSA (78.4%), and SCA (65.4%). GWO's superiority is particularly pronounced in higher-dimensional systems: on the 100- variable Large Sparse System, GWO achieves a 94.1% success rate, compared to 67.3% for GA and 41.2% for SCA. This advantage likely stems from GWO's three-leader hierarchy (α, β, δ), which maintains three simultaneous reference points for guiding the search and an effective strategy for navigating the narrow valleys characteristic of nonlinear equation residual landscapes.

A notable finding is the number of distinct solutions found on the Robot Kinematics system, which has 16 known solutions. GWO locates 15 of the 16 solutions across the 51 runs, while DE finds 13, PSO finds 11, and GA finds only 9. This result demonstrates GWO's ability to explore different basins of attraction across multiple runs, a property valuable for practical applications where identifying all solutions of a nonlinear system is desirable.

Table 8 shows the accuracy of the best solution found by each algorithm on a representative solution of the Robot Kinematics system. Both DE and GWO recover the true parameter values to six decimal places, with GWO achieving a residual of 8.91 × 10−18 compared to DE's 4.56 × 10−16. GA and PSO achieve 4–5 decimal place accuracy, while SCA shows the lowest accuracy with errors in the fourth decimal place and a residual of 3.45 × 10−6.

Table 9 presents the parameter estimation accuracy for all five ODE/PDE models. For each model and algorithm, the table reports the estimated parameter values, the true values, the relative error, the MSE, and the coefficient of determination (R²) between the model output with estimated parameters and the synthetic observed data.

DE achieves the best overall parameter estimation accuracy across all five models, with a mean R² of 0.99994 and maximum relative errors consistently below 0.02%. GWO performs as the close second with a mean R² of 0.99981, followed by PSO (mean R² = 0.99952), GA (mean R² = 0.99834), SSA (mean R² = 0.99761), and SCA (mean R² = 0.99423). The superior performance of DE on this problem class can be attributed to its differential mutation mechanism, which generates trial vectors with self-adaptive step sizes that naturally refine toward the true parameter values as the population converges.

The Lotka-Volterra system provides a particularly instructive case study. The predator-prey dynamics generate oscillatory trajectories that are highly sensitive to parameter values. Small perturbations in α, β, δ, or γ can qualitatively change the phase portrait. Despite this sensitivity, DE recovers all four parameters with relative errors below 0.003%, producing predicted trajectories that are visually indistinguishable from the synthetic observed data. The predicted-vs-observed plots show that DE-estimated parameters capture both the amplitude and phase of the population oscillations with exceptional fidelity, even in the presence of 1% measurement noise.

The PDE parameter estimation problems (heat equation, Burgers' equation) are computationally more expensive due to the need to solve the PDE numerically at each function evaluation. Despite this increased cost, DE and GWO achieve R² values above 0.999 on both PDE problems. The Burgers' equation, which features nonlinear advection that can generate sharp gradient regions (shock-like structures), poses the greatest challenge: SCA achieves only R² = 0.99423 on this problem, reflecting the difficulty of navigating the complex MSE landscape associated with shock-forming solutions.

Table 10 presents the TSP results across all 10 TSPLIB instances. For each algorithm, the best tour length found across 30 independent runs and the mean tour length are reported, along with the gap to the known optimal solution.

The results confirm the clear superiority of GA over other methods for combinatorial optimization problems. GA achieves the best tour length on 10 of 10 TSP instances, with a mean gap-to-optimal of only 0.87% across best solutions. The second-best performer is GWO (mean gap 1.98%), followed by PSO (2.89%), DE (3.78%), SSA (4.23%), and SCA (7.12%). The relatively weak performance of DE on TSP is notable given its dominance on continuous problems; this reflects the fundamental mismatch between DE's differential mutation operator (designed for continuous vector spaces) and the discrete permutation space of TSP, even with random-key encoding.

As the instance size increases, the performance gap between GA and the other algorithms widens. On eil51 (51 cities), GA's best tour is only 0.23% above optimal, while on pr1002 (1,002 cities), the gap increases to 1.70%. For comparison, SCA's gap grows from 3.52% on eil51 to 7.71% on pr1002, indicating substantially worse scalability on discrete problems.

On graph coloring, GA again dominates, finding the chromatic number on 3 of 5 instances (queen5_5, le450_5a, le450_15c) and using only one extra color on the remaining two (DSJC250.5, flat300_28). GWO is the second-best performer, matching GA on queen5_5 and le450_5a. SCA performs worst, using 1–6 extra colors beyond the chromatic number, reflecting its limited ability to handle discrete constraint satisfaction. The computation times reveal that GA's superior solution quality does not come at a disproportionate computational cost; in fact, GA achieves the best time on 4 of 5 instances, likely due to its faster convergence enabled by effective crossover operators for partition-based encoding.

Table 12 presents the overall Friedman rankings computed separately for each problem domain and aggregated across all domains.

The Friedman test applied to the overall rankings yields a test statistic of χ²F = 31.47 with 5 degrees of freedom, corresponding to p < 0.001. It decisively rejects the null hypothesis that all six algorithms perform equally across all problem domains. GWO achieves the best overall rank (2.05), reflecting its strong and consistent performance across all domains, particularly its rank-1 position in nonlinear equation solving. DE ranks second overall (2.43), with the best rank on CEC benchmarks and ODE/PDE estimation. PSO occupies a solid third position (2.76), while GA (3.04) benefits from its dominant performance on combinatorial problems, despite weaker continuous-optimization performance. SSA (3.90) and SCA (5.62) occupy the bottom positions, with SCA consistently ranked last on 3 of 5 domains.

The Nemenyi post hoc analysis reveals a clear pattern of statistically significant differences. At the α = 0.05 significance level, the following pairwise differences are significant: DE vs. SSA (p = 0.0034), DE vs. SCA (p < 0.001), GWO vs. SSA (p = 0.0012), GWO vs. SCA (p < 0.001), PSO vs. SCA (p < 0.001), GA vs. SCA (p < 0.001), and SSA vs. SCA (p = 0.0023). Importantly, the differences among the top three algorithms (DE, GWO, PSO) are not statistically significant, nor is the difference between DE and GA. It suggests that DE, GWO, and PSO form a "top tier" of algorithms whose overall cross-domain performance is statistically indistinguishable, while SCA is statistically inferior to all other algorithms.

The CD for the Nemenyi test at α = 0.05 with k = 6 algorithms and N = 5 problem domains is CD = 2.85. It means that two algorithms must differ by at least 2.85 rank positions to be declared significantly different. Examining the overall ranks GWO (2.05), DE (2.43), PSO (2.76), GA (3.04), SSA (3.90), and SCA (5.62), only the pairs involving SCA exceed this threshold against the top three algorithms. The rank difference between GWO and SSA (1.85) approaches but does not exceed the CD, resulting in a marginally non-significant comparison that would likely become significant with additional problem instances.

The dominance of DE on CEC benchmark functions and ODE/PDE parameter estimation problems can be attributed to two key mechanisms: 1) the differential mutation operator vi = xr1 + F · (xr2 − xr3) automatically adapts the step size and direction to the local geometry of the fitness landscape in the early stages, when the population is dispersed, the difference vectors xr2 − xr3 are large, promoting exploration; as the population converges, the difference vectors shrink, enabling fine-grained exploitation, and 2) the greedy selection mechanism (replacing the target vector only if the trial vector is superior) ensures monotonic improvement of each population member, preventing the premature convergence to suboptimal regions that can affect PSO and GA. The improvement of DE's Friedman rank from 1.52 at D = 10 to 1.35 at D = 100 suggests that the differential mutation mechanism scales well with dimensionality, likely because the number of useful difference vectors grows with the diversity of the population in high-dimensional spaces.

The outstanding performance of GWO on nonlinear systems (97.8% overall success rate, rank 1.00) merits particular discussion. The social hierarchy mechanism in which the three best solutions (α, β, δ) simultaneously guide the movement of all other search agents provides a natural multi-reference-point search strategy that is well-suited to the narrow, steep valleys that characterize the residual landscapes of nonlinear systems. When the residual f(x) = Σ[gi(x)]² approaches zero near a solution, the landscape becomes increasingly ravine-like with extreme conditioning, and GWO's encircling mechanism can track these narrow features more effectively than DE's purely stochastic mutation. Additionally, GWO's parameter-free nature (beyond population size and maximum iterations) eliminates the risk of parameter misspecification that can degrade DE or PSO performance.

The clear superiority of GA on TSP and graph coloring problems (mean gap 0.87% on TSP, Friedman rank 1.00 on TSP and 1.20 on graph coloring) reflects the fundamental advantage of evolutionary operators that are natively designed for discrete, permutation-based, or partition-based representations. The OX operator preserves relative order information from parent tours, producing offspring that inherit high-quality sub-tours while introducing sufficient diversity through recombination. In contrast, the continuous-to-discrete encoding schemes required by PSO, DE, GWO, SSA, and SCA (e.g., random-key encoding and the nearest-position rule) introduce representational overhead that degrades solution quality, particularly on larger instances where the permutation space is vast.

The scalability analysis on CEC 2014 benchmarks reveals divergent trends among the algorithms. DE's performance improves relative to competitors as dimensionality increases, consistent with the theoretical analysis of Zaharie [50], which showed that the variance of the differential mutation operator scales appropriately with the search space volume. Conversely, GWO's relative performance degrades with dimensionality, which may be attributable to the "curse of dimensionality" affecting the encircling mechanism: in high dimensions, the distance between the three leaders and the remaining agents increases, reducing the effectiveness of position updates. SCA exhibits the most inconsistent scalability behavior, likely because the sinusoidal update equation lacks the adaptive mechanisms needed to adjust exploration intensity across scales.

While a formal comparison with classical methods was not conducted in this study, the results support the following qualitative assessment: for well-behaved problems (smooth, unimodal, differentiable), classical methods such as Newton-Raphson and conjugate Gradient will almost certainly outperform metaheuristics in both convergence speed and solution accuracy. However, the metaheuristic algorithms tested here, particularly DE and GWO, demonstrate robust performance across problem types and dimensions without requiring gradient information, analytic Jacobians, or careful initialization. For the 100-variable nonlinear system, for example, Newton-Raphson would require a reliable initial guess within the basin of attraction, which is essentially impossible to obtain a priori. GWO's 94.1% success rate on this system, achieved from random initialization, represents a qualitatively superior practical outcome.

Several limitations of this study should be acknowledged: 1) all algorithms were tested with fixed, canonical parameter settings; adaptive variants (e.g., SHADE, CLPSO, adaptive GWO) would likely yield superior performance but would confound the comparison by conflating algorithmic structure with parameter adaptation strategy, 2) the ODE/PDE parameter estimation experiments used synthetic data with controlled noise; real-world experimental data may present additional challenges including model mismatch, heteroscedastic noise, and missing observations, 3) the study did not include hybrid algorithms that combine metaheuristic global search with local search refinement (e.g., memetic algorithms), which are known to produce superior results on many problem types, and 4) the combinatorial optimization experiments used standard metaheuristic implementations without specialized local search operators (e.g., 2-opt, 3-opt for TSP), which would substantially improve the performance of all algorithms on TSP instances.