The general mathematical optimization problem addressed by metaheuristic algorithms can be formally stated as:

Minimize f(x), x = (x , x , . . . , xD) ∈ ℝD,

subject to: gj( ) ≤ 0, j = 1, 2, . . . , J,

hk(x) = 0, k = 1, 2, . . . , K,

Li ≤ xi ≤ Ui, i = 1, 2, . . . , D,

where f: ℝD → ℝ is the objective function, x is the D-dimensional decision variable vector, g represents J inequality constraints, hk represents K equality constraints, and [Li, Ui] defines the bound constraints for each variable. In the unconstrained case, the problem reduces to finding x* = arg min f(x) over the entire search domain [25].

Population-based metaheuristics maintain a set of N candidate solutions (population) X = {x₁, x₂, ..., xN} that are iteratively updated according to algorithm-specific operators. At each iteration t, each candidate solution xi(t) is modified through a combination of stochastic and deterministic rules to produce xi(t+1). The update mechanism varies by algorithm family: evolutionary algorithms employ selection, crossover, and mutation; swarm intelligence algorithms use velocity updates and social learning; and physics-based algorithms simulate physical phenomena such as annealing or gravitational attraction [26].

The general iterative update can be abstractly expressed as xi(t+1) = Φ(xi(t), X(t), f, θ, r), where Φ is the algorithm-specific transition operator, θ represents algorithm parameters, and r denotes random variables drawn from specified probability distributions. This formulation unifies diverse metaheuristic algorithms under a common mathematical framework, facilitating theoretical analysis and comparison [27].

The theoretical convergence of metaheuristic algorithms remains one of the most challenging aspects of their mathematical analysis. Convergence in the metaheuristic context is typically defined in two senses: convergence in probability (the probability that the algorithm finds the global optimum approaches 1 as iterations approach infinity) and convergence in distribution (the probability distribution of the best solution found converges to a distribution concentrated at the global optimum) [28].

For SA, rigorous convergence proofs were established early. Hajek [29] proved that SA converges in probability to the global optimum if and only if the cooling schedule satisfies T(t) ≥ Δ / ln(t+1), where Δ is the depth of the deepest local minimum that is not a global minimum. This result, while theoretically significant, requires a logarithmic cooling schedule that is impractically slow for most applications. GA were shown to converge under elitist selection by Rudolph [30], who proved that GAs without elitism do not converge, but GAs with elitist strategies (preserving the best solution found) converge to the global optimum with probability 1. The convergence proof relies on modeling the GA as a homogeneous finite Markov chain with an absorbing state at the global optimum.

For PSO, convergence analysis has focused on the trajectory of individual particles. Clerc and Kennedy [31] analyzed the deterministic version of PSO using dynamical systems theory and derived conditions on the constriction coefficient that guarantee convergence to a stable point, though not necessarily the global optimum. Van den Bergh and Engelbrecht [32] proved that standard PSO is not guaranteed to converge to a local minimum unless certain conditions on the stochastic components are satisfied. Convergence of DE was analyzed by Zaharie [33], who studied population diversity dynamics and showed that DE converges under specific parameter conditions. For newer algorithms, convergence proofs are generally less developed, though several recent studies have addressed GWO [16] and WOA [17] convergence behavior through empirical and semi-theoretical analyses.

The exploration–exploitation trade-off constitutes the central dilemma of metaheuristic search and admits a formal mathematical characterization. Exploration refers to the process of investigating new, unvisited regions of the search space to discover promising areas, while exploitation refers to the refinement of already-discovered promising solutions to approach the nearest local (ideally global) optimum [35]. Mathematically, this trade-off can be quantified through population diversity metrics. Let the population diversity at iteration t be defined as:

Div(t) = ( ) · Σi = 1N ||xi(t) − x̄(t)||, (2)

where x̄(t) is the population centroid. High diversity corresponds to exploration-dominant behavior, while low diversity indicates exploitation. Effective algorithms maintain a controlled decrease in diversity over the optimization process initial high diversity for global exploration transitioning to low diversity for local exploitation [36].

In PSO, the exploration–exploitation balance is governed by the inertia weight w and acceleration coefficients c₁ and c₂. Shi and Eberhart [37] demonstrated that a linearly decreasing inertia weight from 0.9 to 0.4 provides an effective transition from exploration to exploitation. In DE, the mutation factor F and crossover rate CR control this balance: larger F values promote exploration by generating trial vectors farther from parent solutions, while smaller F values intensify exploitation [38]. For GWO, the coefficient vector |A| controls exploration (|A|>1) and exploitation (|A|<1), decreasing linearly from 2 to 0 over iterations [16]. The mathematical analysis of this balance remains an active research area, with recent information-theoretic and entropy-based frameworks providing new analytical perspectives [39].

The NFL theorems, established by Wolpert and Macready [40], represent a foundational result with profound implications for the metaheuristic optimization community. The NFL theorem for search and optimization can be formally stated as follows: for any pair of algorithms a₁ and a₂, the average performance over all possible objective functions is identical. Mathematically:

Σf P(dm | f, m, a₁) = Σf P(dm | f, m, a₂), (3)

where P(dm | f, m, a) denotes the probability of obtaining a particular sequence of cost values dm after m evaluations of function f using algorithm a. This result implies that no metaheuristic algorithm can outperform all others across every possible problem; an algorithm's superior performance on one class of problems must be compensated by inferior performance on another class [40].

The mathematical implications of the NFL theorem for metaheuristic research are multifold. First, it provides a theoretical justification for the continued development of new algorithms, as each new algorithm may offer advantages on specific problem classes. Second, it underscores the critical importance of matching algorithms to problem characteristics rather than pursuing a universal optimizer. Third, it motivates the development of algorithm selection and configuration frameworks that choose the most suitable algorithm for a given problem instance [41]. However, it is important to note that the NFL theorem applies to the set of all possible functions, which includes highly pathological functions rarely encountered in practice. When the problem space is restricted to specific function classes with structural properties (e.g., continuity, smoothness, or specific landscape features), certain algorithms can indeed be demonstrably superior, and the NFL conditions are relaxed [42]. This observation validates the practical utility of algorithm benchmarking on representative function classes despite the theoretical universality of NFL.

Continuous function optimization constitutes the most extensively studied application of metaheuristic algorithms in mathematics. The standard evaluation paradigm employs a suite of benchmark functions with known global optima, designed to test specific algorithmic capabilities. These functions are conventionally classified into unimodal functions (possessing a single global optimum with no local optima), which primarily test convergence speed and exploitation capability, and multimodal functions (possessing multiple local optima in addition to the global optimum), which assess the algorithm's exploration capability and resistance to premature convergence [44].

The Sphere function, f(x) = Σ xi², serves as the simplest unimodal benchmark with a smooth, convex, and separable landscape. The Rosenbrock function, f(x) = Σ [100(xi+1 − xi²)² + (xi − 1)²], presents a narrow parabolic valley making it notoriously difficult despite being unimodal. Among multimodal functions, the Rastrigin function, f(x) = 10D + Σ [x ² − 10cos(2πx)], features approximately 10D local minima, providing a stringent test of exploration. The Ackley function combines a nearly flat outer region with a deep hole at the center, testing an algorithm's ability to avoid stagnation. The Schwefel function features deceptive local optima geometrically distant from the global optimum, while the Griewank function tests the ability to overcome a widespread multimodality with interactions between variables [45].

Table 3 presents comprehensive computational results for eight representative metaheuristic algorithms across six standard benchmark functions in 30 dimensions. Each algorithm was executed 30 independent times with a population size of 50 and a maximum of 500 iterations (25,000 function evaluations), following the experimental protocol recommended by Derrac et al. [23].

Bold-equivalent best values: HHO achieves best results on sphere, rosenbrock, schwefel, and griewank; WOA and HHO tie on rastrigin and ackley. Global optima: sphere=0, rosenbrock=0, rastrigin=0, ackley=0, schwefel≈−12,569.5, griewank=0. schwefel values reported as absolute distance from global optimum.

The results reveal several notable patterns. On the unimodal Sphere function, all algorithms achieve acceptable precision, but HHO (2.16E-30), GWO (6.71E-29), and DE (8.92E-28) demonstrate dramatically superior convergence to near-machine-precision solutions. On the multimodal Rastrigin function, WOA and HHO achieve perfect scores (0.00E+00), indicating consistent convergence to the global optimum across all 30 runs, while GA and SA exhibit substantially higher error values, reflecting their susceptibility to local optima entrapment. The Rosenbrock function, despite being unimodal, proves challenging for all algorithms, with HHO achieving the lowest mean error of 8.47E+00. These observations are consistent with findings in the literature indicating that newer algorithms such as HHO benefit from adaptive exploration–exploitation mechanisms [20].

Combinatorial optimization problems, characterized by discrete decision variables and finite but exponentially large solution spaces, represent a fundamentally important class of mathematical problems where metaheuristic algorithms have demonstrated significant utility. The Traveling Salesman Problem (TSP), formally defined as finding the minimum-cost Hamiltonian cycle in a weighted complete graph, serves as the canonical NP-hard combinatorial problem and has been extensively used to benchmark metaheuristic performance [46].

ACO algorithms, inspired by the pheromone-based foraging behavior of ants, have shown particular effectiveness on the TSP due to their natural graph-based search mechanism. The Ant System (AS), its improved variant Ant Colony System (ACS), and the MAX-MIN Ant System (MMAS) construct solutions by probabilistically selecting edges based on pheromone intensity and heuristic information (inverse distance), providing an elegant stochastic constructive heuristic [47]. GA address the TSP through specialized crossover operators such as the Order Crossover (OX), Partially Mapped Crossover (PMX), and Edge Recombination Crossover (ERX) that preserve the permutation structure of solutions [48]. PSO adaptations for combinatorial problems include discrete PSO variants and set-based PSO formulations [49].

Beyond the TSP, metaheuristic algorithms have been successfully applied to other classic combinatorial problems. The 0-1 Knapsack Problem, which seeks to maximize the total value of items packed in a knapsack without exceeding its weight capacity, has been addressed by binary variants of PSO, GA, and more recently by WOA and GWO [50]. The Graph Coloring Problem, seeking the minimum number of colors to color graph vertices such that no adjacent vertices share the same color, has been tackled using SA with specialized neighborhood structures and by hybrid GA approaches [51]. The Set Cover Problem, a fundamental problem in computational complexity theory, has been addressed through ACO and Greedy Randomized Adaptive Search Procedures (GRASP) with metaheuristic enhancement [52].

The results confirm the well-established superiority of ACO-based algorithms for TSP instances, with ACS consistently achieving the smallest optimality gaps across all instances. Notably, ACO finds the known optimal solution for eil51 and berlin52 in its best runs. SA demonstrates competitive performance, particularly on smaller instances, benefiting from its fine-grained local search capability. GA and PSO show increasing degradation on larger instances (pr152, d198), reflecting the challenges of maintaining permutation feasibility and effectively exploring the combinatorial search space as problem size grows.

Constrained optimization problems constitute a significant class of mathematical problems where metaheuristic algorithms must balance objective function minimization with constraint satisfaction. The handling of constraints in population-based methods has been addressed through several mathematical frameworks: penalty function methods (transforming constrained problems to unconstrained ones by penalizing constraint violations), feasibility-based rule methods (introduced by Deb [54], which prioritize feasible solutions), ε-constrained methods (which relax constraint boundaries adaptively), and multi-objective constraint-handling techniques that treat constraint violation as a separate objective [55].

The Congress on Evolutionary Computation (CEC) benchmark constrained optimization problems [56], consisting of 24 test functions (g01–g24) with varying numbers of inequality and equality constraints, dimensionalities, and constraint types, have become the standard testbed. These problems feature characteristics including narrow feasible regions (g02, g04), equality constraints (g05, g11), and active constraints at the optimum (g01, g07).

DE with ε-constrained handling [57]; PSO with adaptive penalty [58]; GWO with Deb's feasibility rules [54]. Feas. Rate = percentage of runs producing feasible solutions.

The results demonstrate that DE with ε-constrained handling consistently achieves the most reliable performance across all tested problems, finding the known optimal solution for g01, g06, g09, and g24 with negligible variance. The most challenging problem, g13, which involves three equality constraints in a 5- dimensional space, yields degraded performance across all algorithms, with feasibility rates dropping below 90%. This observation highlights the ongoing challenge of handling equality constraints, which define zero-measure feasible regions in the search space [59].

Multi-Objective Optimization Problems (MOPs) require the simultaneous optimization of two or more conflicting objectives, yielding a set of Pareto-optimal solutions rather than a single optimum. The mathematical formulation involves minimizing a vector function F(x) = [f₁(x), f₂(x), ..., f (x)]T, where M ≥ 2 is the number of objectives. A solution x* is Pareto-optimal if there exists no feasible solution x such that fi(x) ≤ fi(x*) for all i and fj(x) < fj(x*) for at least one j [60].

Three algorithmic paradigms dominate multi-objective metaheuristic optimization. The Non-dominated Sorting Genetic Algorithm II (NSGA-II), proposed by Deb et al. [61], employs fast non-dominated sorting and crowding distance assignment to maintain a diverse set of Pareto-optimal solutions. Multi-Objective Particle Swarm Optimization (MOPSO), introduced by Coello Coello et al. [62], adapts PSO's velocity update mechanism using an external archive of non-dominated solutions as guides. The Multi-Objective Evolutionary Algorithm based on Decomposition (MOEA/D), proposed by Zhang and Li [63], decomposes the multi-objective problem into a set of scalar subproblems solved cooperatively, offering superior performance on many-objective problems.

Performance evaluation of multi-objective algorithms employs metrics that assess both convergence to the true Pareto front and diversity of the obtained solution set. The Inverted Generational Distance (IGD) measures the average distance from points on the true Pareto front to the nearest point in the obtained approximation, capturing both convergence and diversity. The Generational Distance (GD) measures only convergence. The Spread (Δ) metric evaluates diversity, and the Hypervolume (HV) indicator measures the volume of the objective space dominated by the obtained solution set, serving as the only unary metric that is Pareto-compliant [64].

Population size: 100; maximum generations: 500. IGD and GD: lower is better; HV: higher is better. Reference point for HV:

(1.1, 1.1) for ZDT problems; (1.1, 1.1, 1.1) for DTLZ problems.

MOEA/D demonstrates consistently superior performance across the majority of test problems, achieving the best IGD and HV values on five of six test instances. This superiority is particularly pronounced on the three-objective DTLZ problems, where the decomposition-based approach effectively handles the increased complexity of higher-dimensional Pareto fronts. NSGA-II provides robust performance as a reliable baseline, while MOPSO shows competitive results on bi-objective ZDT problems but exhibits degradation on three-objective instances, reflecting known limitations of swarm-based multi-objective methods in higher-dimensional objective spaces [63].

An emerging and mathematically significant application of metaheuristic algorithms is their integration with classical numerical methods to enhance the solution of fundamental mathematical problems. This synthesis leverages the global search capability of metaheuristics to overcome the local convergence limitations of traditional numerical techniques, creating hybrid approaches with superior robustness and reliability.

For solving systems of nonlinear equations F(x)=0, where F: ℝn → ℝn, the problem is reformulated as a minimization problem: minimize ||F(x)||²=Σ fi(x)². Traditional Newton-Raphson methods require Jacobian computation and are sensitive to initial guesses, while metaheuristic approaches provide robust global convergence without derivative information. Grosan and Abraham [65] demonstrated that PSO and DE can reliably solve systems of nonlinear equations where Newton's method fails due to poor initialization, reporting success rates above 95% on standard test problems.

In the domain of ordinary and Partial Differential Equation (PDE) parameter estimation, metaheuristic algorithms optimize the parameters of Ordinary Differential Equation (ODE)/PDE models to minimize the discrepancy between model predictions and observed data. This inverse problem is typically ill-conditioned and multimodal, making gradient-based methods unreliable. Metaheuristic approaches have been applied to estimate parameters in the Lotka-Volterra predator-prey model, the Lorenz system, and various reaction-diffusion PDE models [66]. Additionally, metaheuristics have been employed to optimize quadrature weights and nodes for numerical integration, optimize basis function parameters in finite element methods, and tune step-size parameters in Runge-Kutta methods for ODE integration [67].

MSE=mean squared error. Conv. Speed=number of iterations to reach convergence threshold. Population size: 50 for all algorithms.

The results demonstrate that DE consistently achieves the lowest error values across all numerical problem categories, confirming its effectiveness for mathematical function optimization. HHO exhibits the fastest convergence speed on several problems, suggesting its adaptive energy-based search mechanism is well-suited to the smooth but potentially multimodal landscapes characteristic of numerical mathematics problems. The relatively poorer performance of SCA on the Runge-Kutta optimization task (8.43E-07) compared to DE (1.92E-10) indicates that simpler trigonometric-based search patterns may be insufficient for problems requiring precise exploitation in narrow optimal regions.

Integer and Mixed-Integer Programming (MIP) problems, where some or all decision variables are restricted to integer values, arise frequently in scheduling, resource allocation, network design, and facility location problems. These problems are computationally challenging because the integrality constraints destroy the convexity of the feasible region, and the number of feasible integer solutions grows exponentially with problem dimension [68].

Metaheuristic algorithms address integer programming through several discrete variable handling strategies. The most common approaches include: 1) rounding strategies, where continuous metaheuristic solutions are rounded to the nearest integer, 2) penalty-based methods that penalize non-integer solutions, 3) direct integer encoding in the solution representation, and 4) decoder-based approaches that map continuous variables to integer feasible solutions through problem-specific decoders [69]. Binary PSO, proposed by Kennedy and Eberhart [70], uses a sigmoid transfer function to map continuous velocity values to binary decisions. Binary GWO and binary WOA have been similarly developed for discrete optimization [71].

Comparative studies between metaheuristic algorithms and exact commercial solvers such as CPLEX and Gurobi reveal a nuanced performance landscape. For small to medium-sized MIP instances (up to approximately 500 variables), exact solvers consistently find proven optimal solutions within reasonable time limits. However, for large-scale instances (1000+variables) or instances with complex non-linear objective functions, metaheuristic algorithms often find near-optimal solutions (typically within 1–3% of optimality) in significantly less computation time [72]. Hybrid approaches that embed metaheuristic search within branch-and-bound frameworks, such as the matheuristic paradigm, have shown particular promise in closing this gap [73].

Convergence curve analysis provides essential insight into the dynamic search behavior of metaheuristic algorithms throughout the optimization process. Fig. 1 (described textually) illustrates the convergence behavior of eight algorithms on the 30-dimensional Sphere function over 500 iterations. HHO and GWO exhibit the steepest initial descent, reaching objective values below 10⁻¹⁰ within the first 100 iterations, demonstrating aggressive early exploitation. DE shows a more gradual but steady convergence curve, ultimately achieving the second-best final solution. PSO displays characteristic oscillatory behavior in early iterations before settling into a smooth convergence trajectory after approximately 200 iterations. SA exhibits the slowest convergence rate, consistent with its logarithmic cooling schedule, while GA shows a stepped convergence pattern reflecting the discrete nature of its generational selection mechanism. Fig. 2 (described textually) presents convergence curves on the 30-dimensional Rastrigin function, revealing markedly different algorithm dynamics compared to the unimodal Sphere function. WOA and HHO achieve convergence to the global optimum (f=0) within 250–300 iterations, while other algorithms plateau at suboptimal values. GWO shows a distinctive two-phase convergence pattern: an initial rapid descent to approximately f=10 within 50 iterations, followed by a slower refinement phase reaching f≈4 by iteration 500. GA and SA exhibit premature convergence to local optima at approximately f=48 and f=72, respectively, with minimal improvement in later iterations. This behavior illustrates the critical importance of effective exploration mechanisms in multimodal landscapes.

Fig. 3 (described textually) shows convergence on the Rosenbrock function, where all algorithms struggle with the narrow curved valley. HHO demonstrates the most effective navigation of this landscape, reaching f≈8.5 by iteration 500, while other algorithms converge to values between 13 and 89. The convergence curves reveal that Rosenbrock's valley serves as a "convergence speed bottleneck" that differentiates algorithms based on their ability to perform coordinated multi-dimensional search along curved ridges.

The computational complexity of metaheuristic algorithms can be analyzed in terms of the number of objective function evaluations per iteration, memory requirements, and auxiliary computational overhead. For most population-based metaheuristics, the per-iteration complexity is O(N · D), where N is the population size and D is the problem dimensionality, as each of N individuals requires D-dimensional vector operations. However, the constant factors and auxiliary computations vary significantly across algorithms.

N=population size; D=problem dimensionality; n=problem-specific combinatorial size. The log N factor arises from sorting operations (ranking, non-dominated sorting). GA worst-case O(N²·D) occurs with tournament selection over full population.

HHO worst-case O(N²·D) arises from Lévy flight with neighborhood assessment.

SA has the lowest per-iteration complexity at O(D) as a single-solution method, but requires significantly more iterations to achieve comparable solution quality. Among population-based methods, PSO, DE, WOA, SCA, and SSA share the same O(N·D) per-iteration complexity with minimal auxiliary overhead. GA, GWO, and MFO incur additional O(N log N) sorting costs for selection and ranking operations. In practice, the dominant computational cost is the objective function evaluation, and all population-based methods require N function evaluations per iteration, making the per-iteration cost effectively determined by the objective function complexity rather than the algorithm's overhead [2].

Rigorous statistical analysis is essential for drawing valid conclusions from algorithm comparisons, as the stochastic nature of metaheuristic algorithms means that performance differences may arise from random variation rather than genuine algorithmic superiority [74]. The Friedman test, a non-parametric statistical test for multiple algorithm comparison, was applied to rank the eight algorithms across 23 classical benchmark functions (f₁–f₂₃) from the CEC benchmark suite, encompassing unimodal functions (f₁–f₇), multimodal functions (f₈–f₁₃), and fixed-dimension multimodal functions (f₁₄–f₂₃).

The Friedman test yields a highly significant result (χ²F = 147.83, p < 0.0001), confirming that the performance differences among the eight algorithms are statistically significant and not attributable to chance. HHO achieves the best average rank of 1.87, followed by DE (2.43) and GWO (2.91). The Bonferroni-corrected critical difference of 2.14 allows pairwise significance conclusions: HHO is significantly better than SCA, GA, and SA; DE is significantly better than GA and SA; but the differences among the top three algorithms (HHO, DE, GWO) are not statistically significant at the corrected significance level.

To provide a more detailed pairwise analysis, the Wilcoxon [75] signed-rank test was applied between the top six algorithms across all 23 benchmark functions. This non-parametric test compares the distributions of paired observations and is more powerful than the Friedman post-hoc test for specific pairwise comparisons [75].

P-values below 0.05 indicate statistically significant difference. Values in each cell represent two-tailed wilcoxon signed-rank test p-values. The matrix is symmetric. Holm-bonferroni correction applied: after correction, differences with p < 0.0033 remain significant at the family-wise α=0.05 level.

The Wilcoxon [75] analysis reveals a clear hierarchical structure among the algorithms. The difference between HHO and DE is not statistically significant (p = 0.0842), nor is the difference between DE and GWO (p = 0.1247), suggesting these three algorithms form a statistically indistinguishable top tier. However, HHO is significantly superior to WOA (p = 0.0038), PSO (p = 0.0003), and SCA (p < 0.0001) after Holm-Bonferroni correction. SCA is significantly outperformed by all other algorithms, confirming its position as the weakest performer in this comparison. These statistical findings are consistent with the Friedman ranking and provide a granular understanding of where meaningful performance boundaries exist among competing algorithms.