Theoretical Analysis of Multi-Objective Evolutionary Algorithms on Integer Spaces with Local Optima

1 Introduction

Multi-objective optimization (MOO) problems arise in many real-world applications, where multiple conflicting objectives must be optimized simultaneously. The goal in MOO is to identify a set of Pareto-optimal solutions that represent trade-offs among conflicting objectives. Multi-objective evolutionary algorithms (MOEAs) have become one of the most widely used approaches for solving MOO problems [Deb, 2001; Zhou et al., 2019]. By maintaining a population of candidate solutions, MOEAs can approximate the Pareto front in a single run. Well-established MOEAs include the Global Simple Evolutionary Multi-objective Optimizer (GSEMO) [Giel, 2003], Non-dominated Sorting Genetic Algorithm II (NSGA-II) [Deb et al., 2002], S-Metric Selection Evolutionary Multi-objective Optimization Algorithm (SMSEMOA) [Beume et al., 2007], MOEA based on Decomposition (MOEA/D) [Zhang and Li, 2007], and SPEA2 [Zitzler et al., 2001]. These MOEAs have been successfully applied to many real-world domains, such as neural architecture search [Lu et al., 2019] and production scheduling [Ojstersek et al., 2020]. Despite their widespread use, theoretical study of MOEAs has lagged behind practical applications.

In recent years, theoretical analysis of MOEAs, particularly running time analysis, has attracted considerable attention, aiming to reveal, for example, which evolutionary mechanisms are provably effective and under which conditions they ensure efficient optimization [Auger and Doerr, 2011; Neumann and Witt, 2010; Zhou et al., 2019; Doerr and Neumann, 2020]. Early studies mainly focused on simple MOEAs (e.g., GSEMO), while recent work shifted to practical MOEAs, including MOEA/D [Li et al., 2015; Corus et al., 2025], NSGA-II [Zheng et al., 2022], SPEA2 [Ren et al., 2024], SMS-EMOA [Bian et al., 2025], and NSGA-III [Wietheger and Doerr, 2024a]. Beyond the algorithm framework, theoretical studies have investigated the influence of specific components of MOEAs, such as selection mechanisms [Doerr et al., 2025a; Ren et al., 2026], archiving [Bian et al., 2024; Ren et al., 2025; Ren et al., 2026], population update mechanisms [Bian et al., 2025], mutation operators [Doerr et al., 2025b], and crossover operators [Doerr and Qu, 2023; Dang et al., 2024a], thereby providing insights for efficient algorithm design. At the problem-setting level, theoreticians [Wietheger and Doerr, 2024b; Opris, 2025] have analyzed many-objective optimization problems (with more than two objectives), aiming to explain the increasing difficulty of MOEAs as the number of objectives grows.

However, these studies in MOEAs are primarily focused on synthetic pseudo-Boolean functions, such as ONEMINMAX [Giel and Lehre, 2010], LEADINGONESTRAILINGZEROES [Laumanns et al., 2004b], COUNTONESCOUNTZEROS [Laumanns et al., 2004b], and ONEJUMPZEROJUMP [Doerr and Zheng, 2021] benchmarks, all of which are obtained by extending classical single-objective functions [Liang and Li, 2025]. For instance, ONEMINMAX simultaneously maximizes the number of 1-bits (extended from ONEMAX [Droste et al., 2002]) and the number of 0-bits of a binary bit string; LEADINGONESTRAILINGZEROES simultaneously maximizes the number of leading 1-bits (extended from LEADINGONES [Rudolph, 1997]) and the number of trailing 0-bits; and ONEJUMPZEROJUMP generalizes the JUMP [Doerr and Neumann, 2020] function to a bi-objective setting. There also exist other benchmarks, such as ONETRAPZEROTRAP [Dang et al., 2024b], and multi-objective REALROYALROAD [Dang et al., 2023]. Beyond these synthetic benchmarks, only a few studies consider MOEAs on multi-objective combinatorial optimization problems, such as multi-objective knapsack [Laumanns et al., 2004a], minimum spanning tree [Cerf et al., 2023], and shortest path [Horoba, 2009].

Practical MOO problems, however, are not limited to pseudo-Boolean domains, which can lead to MOEAs exhibiting behaviors that differ substantially from those observed on synthetic problems [Liang and Li, 2026]. Recently, there has been progress in extending the running time analysis of MOEAs to integer spaces. To our knowledge, the only work in this vein is from [Doerr et al., 2025b] which analyzes a benchmark with an unbounded integer space (presented by [Rudolph, 2023]). This problem is defined over n-dimensional integer vectors and aims to minimize the distance of a solution to two target points. It induces a simple, unimodal fitness landscape in which the Pareto optimal solutions lie on the line connecting the two target points, and any move towards the Pareto set results in an improvement in fitness. Doerr et al. [2025b] analyzed the performance of (G)SEMO with different mutation operators, and showed that unit-step mutation can be slow when initial solutions are far from the Pareto front, whereas exponential-tail and powerlaw mutations yield consistently better performance across different problem parameters and starting points.

In this paper, we present a class of integer-valued test functions that can incorporate diverse problem features, serving as a benchmark for the theoretical analysis of MOEAs on integer spaces. The main contributions are as follows:

• We introduce a class of multi-objective integer optimization functions, termed the L1DM-2SIP problem, which has an analyzable structure as well as practical relevance. This benchmark is inspired by the multi-point distance minimization problem (MPDMP) [K¨oppen and Yoshida, 2007; Ishibuchi et al., 2010; Fieldsend et al., 2021] in the continuous domain, which has been used to visually investigate the behavior of MOEAs. We instantiate a simple function of L1DM-2SIP that features local-optimal regions (Definition 5) and thus complements the current simple unimodal benchmark [Rudolph, 2023; Doerr et al., 2025b].

Compared to pseudo-Boolean benchmarks, the proposed benchmark exhibits a high degree of consistency between the solution and objective spaces. This consistency, arising from the L1 distance linking solution positions with dominance relations in objective space, makes the search process more intuitive and traceable in the solution space, especially when progressing from the local Pareto optimal set to the global Pareto optimal set.

• Note that the presented L1DM-2SIP problem is of practical relevance and can model real-world problems such as site location selection in urban planning, where the goal is to identify possible sites that simultaneously minimize the Manhattan distances to multiple facilities.
• For the considered L1DM-2SIP function, we prove in Theorem 1 that the widely-studied GSEMO with unitstep mutation [Doerr et al., 2025b] can be trapped in local-optimal regions and fails to identify the Pareto front. However, this difficulty can be overcome by introducing an ageing mechanism [Li et al., 2025] (as shown in Theorem 2) or employing heavy-tailed mutation operators, including exponential-tail and powerlaw mutations [Doerr et al., 2025b] (as shown in Theorem 3), which enables GSEMO escape local-optimal regions and successfully identify the Pareto front within a bounded expected running time.
• To show the extendability of the presented benchmark, we also provide an additional instantiation with a more complex landscape featuring numerous local optima (Definition 6). This instantiation more closely resembles the landscape characteristics of multi-modal problems commonly used in empirical research, such as those in the ZDT [Zitzler et al., 2000] and DTLZ [Deb et al., 2005] suites.

Finally, we also conduct empirical studies to complement the theoretical results. We hope this work is a step forward for the theoretical study of MOEAs on problems that are closely related to those commonly investigated in empirical research.

2 Proposed Benchmark L1DM-2SIP

In this section, we first introduce the basic concepts of multiobjective optimization, and then describe the bi-objective optimization problem in integer space studied by Doerr et al. [2025b]. Finally, we present a new integer-valued benchmark problem L1DM-2SIP, and provide one of its problem instances. In this paper, we use [p..q] to represent [p, q] ∩ Z for p ≤ q and p, q ∈ Z.

2.1 Multi-objective Optimization

Multi-objective optimization aims to simultaneously optimize two or more objective functions. In this paper, we focus on minimization. In multi-objective optimization, objectives are typically conflicting, meaning that there is no canonical complete order in the solution space X. To compare solutions, we use the domination relationship given in Definition 1. A solution is Pareto optimal if no other solution in X dominates it. The set of objective vectors corresponding to all Pareto optimal solutions forms the Pareto front (PF). The goal of multi-objective optimization is to find the Pareto front or a good approximation of it. In addition to global optimality, we also consider local optimality [Lapucci et al., 2023]. Given a neighborhood structure on X, a solution is locally Pareto-optimal if it is not dominated by any of its neighboring solutions. The definition of neighboring solutions studied in this paper is given by Definition 3. That is, two solutions are neighboring if and only if they differ by 1 in exactly one dimension. Note that a locally Pareto-optimal solution is not necessarily Pareto-optimal.

Definition 1. Let f = (f1, f2, . . . , fm) : X → Rm be the objective vector. For two solutions x and y ∈ X :

• x weakly dominates y (denoted as x ⪯ y) if for any 1 ≤ i ≤ m, fi(x) ≤ fi(y);
• x dominates y (denoted as x ≺ y) if x ⪯ y and fi(x) < fi(y) for some i;
• x and y are incomparable if neither x ⪯ y nor y ⪯ x.

Definition 2 (Local Pareto Optimality). Let N(x) ⊆ X denote the neighborhood of a solution x ∈ X. A solution x is called locally Pareto-optimal (with respect to N) if there exists no solution y ∈ N(x) such that y ≺ x.

Definition 3 (ℓ∞-Neighborhood). For x ∈ X ⊆ Zn , we have

2.2 Existing Integer-Space Benchmark

Doerr et al. [2025b] conducted the first running time analysis of (G)SEMO with an unbounded integer search space, using the problem benchmark (shown in Definition 4) presented by Rudolph [2023]. Their results showed that unit-step mutation is slow when starting far from the PF, while exponential-tail and power-law mutations perform better.

Definition 4 ([Rudolph, 2023]). Consider the bi-objective optimization problem fg : Zn → N20 with

where ∥ · ∥1 denotes the ℓ1-norm, z1 = (c, 0, ..., 0)T, and z2 = (-c, 0, ..., 0)T.

The Pareto set of fg is X∗ g = {(k, 0, . . . , 0)T ∈ Zn | k ∈ [-c..c]}, and the Pareto front of fg is F ∗ g = {(k, 2c - k)T ∈ N2 | k ∈ [0..2c]}, where |F ∗ | = 2c + 1. This problem is simple and exhibits a unimodal global fitness landscape, in which every point that moves towards the Pareto set achieves an improvement in fitness. Next, we will introduce a new benchmark with local optima.

2.3 Proposed Integer-Space Benchmark L1DM-2SIP

We consider a multi-objective optimization problem that simultaneously minimizes the L1-distance regarding a point to two sets of m points in an integer spac, called L1-distance minimization to two sets of integer points (L1DM-2SIP). Each set contains multiple candidate points, and typically each point in one set is paired with a corresponding point in the other set. One instance of the proposed L1DM-2SIP with two pairs of points is given as follows.

Definition 5 (L1DM-2SIP: instance I). Consider the problem of L1-distance minimization to two sets ({z1, u1} and {z2, u2}) of integer points: Zn → N20 with

where ∥ · ∥1 denotes the ℓ1-norm, z1 = (-c,-b, 0, ..., 0)T, z2 = (c,-b, 0, ..., 0)T, u1 = (-a, b, 0, ..., 0)T, u2 = (a, b, 0, ..., 0)T, a < c, and a + 2b - c > 1.

Under the conditions of a < c and a + 2b - c > 1, L1DM2SIP admits both Pareto-optimal and locally Pareto-optimal solutions. The set of locally Pareto-optimal solutions is given by XL = {(x1, -b, 0, . . . , 0)T ∈ Zn | x1 ∈ [-c..c]}, with |XL| = 2c + 1. The Pareto set of L1DM-2SIP is X∗ = {(x1, b, 0, . . . , 0)T ∈ Zn | x1 ∈ [-a..a]}, and the Pareto front is F ∗ = {(k, 2a-k)T ∈ N2 | k ∈ [0..2a]}, with |X∗ | = |F ∗ | = 2a + 1. Figure 1 shows an example of the above function with two decision variables, along with its fitness landscape under the Pareto dominance relation (i.e., Pareto landscape [Liang et al., 2024]).

In the left sub-figure of Figure 1, U1 = {(x1, b - c + a, 0, · · · , 0)T | x1 ∈ [-a..a]} and U2 = {(x1, b + c - a, 0, · · · , 0)T | x1 ∈ [-a..a]} (blue lines) have the same objective values as XLa = {(x1,-b, 0, . . . , 0)T ∈ Zn | x1 ∈ [-a..a]} (i.e., the set of solutions that satisfy x1 ∈ [-a..a] in XL and also a mapping of U1 ontoXL). The solutions in U3 = {(x1, x2, 0, · · · , 0)T | x1 ∈ [-c.. - a], x2 ∈ {b - c - x1, b + c + x1}} (the green line on the left) have the same objective values as (-c, b, 0, · · · , 0)T , and the solutions in U4 = {(x1, x2, 0, · · · , 0)T | x1 ∈ [a..c], x2 ∈ {b - c + x1, b + c-x1}} (the green line on the right) have the same objective values as (c, b, 0, · · · , 0)T . Thus, the solutions in the boundary U = U1 ∪ U2 ∪ U3 ∪ U4 (the boundary formed by the blue and green lines) and XL are incomparable or have the same objective values. The internal region enclosed by U is denoted as U + (the shaded area), which dominates the local Pareto optimal set XL. That is, any solution in U + is not dominated by the local Pareto optimal solutions in XL, while for any local Pareto optimal solution in XL, there exist solutions in U + dominating it.

It is worth noting that when c - a ≥ 3, two new sets of local Pareto optimal solutions XL1 and XL2 appear, where XL1 consists of all (x1, x2, . . . , xn)T ∈ Zn with -c ≤ x1 ≤ -a-3 and ⌈ a-c 2 ⌉+1 ≤ x2 ≤ ⌈ a-c 2 ⌉-x1-a-2, andXL2 consists of all (x1, x2, . . . , xn)T ∈ Zn with a + 3 ≤ x1 ≤ c and ⌈ a-c 2 ⌉ + 1 ≤ x2 ≤ ⌈ a-c 2 ⌉ +x1 -a- 2. However, XL1 and XL2 do not form a peak and will not affect the analysis of GSEMO, so the local Pareto optimal solutions mentioned later all refer to XL.

3 Running Time Analysis of GSEMO on L1DM-2SIP

In this section, we analyze the expected running time of GSEMO (an MOEA widely used in theoretical analysis of MOEAs [Neumann and Witt, 2010; Zhou et al., 2019]) with the same mutation operators as those used in [Doerr et al., 2025b] for solving the proposed integer-space benchmark L1DM-2SIP. As presented in Algorithm 1, GSEMO maintains a population of unfixed size to store all solutions that are not dominated by any other solution encountered so far. It starts with an initial solution randomly selected from the solution space. In each iteration, a solution x is randomly selected from the current population, and a new solution y is generated by applying mutation to x. If no solution in the current population dominates y, all individuals weakly dominated by y will be deleted from the population and y is added to the population.

For the mutation operator, we consider full-dimensional mutation, which independently modifies each dimension xi of a solution x with probability 1/n. That is, each dimension yi of the offspring solution y generated from x is

where Zi is an integer-valued random variable whose distribution depends on the chosen mutation operator.

Lemma 1 shows the maximum population size maintained by GSEMO during the optimization of L1DM-2SIP, which will be used in the proofs of Theorems 2 and 3. The proof of Lemma 1 proceeds by first dividing the solution space into different regions based on whether the solution is closer to z1 (z2) or u1 (u2), the integer points used in the definition of L1DM-2SIP. Then, in each region, the size of the nondominated solution set that can be retained in the population is determined based on the already established form of the objective functions. Since there are dominance relationships between the non-dominated solution sets in different regions, the size of the largest non-dominated solution set among them represents the maximum population size. The detailed proofs are provided in the Appendix due to space limitation.

Lemma 1. When using GSEMO to optimize L1DM-2SIP, the population size is at most 2c + 1.

In the following, we will first illustrate the risk of the GSEMO algorithm with unit-step mutation getting stuck in local Pareto optima when optimizing L1DM-2SIP, and then demonstrate that using an ageing strategy or heavy-tailed mutation operators can help escape from local Pareto optima and find the Pareto front.

3.1 Unit-stepMutation fails

We first consider GSEMO equipped with the basic unit-step mutation operator, and show in Theorem 1 that it runs the risk of getting stuck in local Pareto optima, preventing it from ever finding the Pareto front.

Unit-stepmutation. The random variable Zi in Eq. (1) is 1 with probability 1/2, and -1 otherwise.

Algorithm 1 Global Simple Evolutionary Multi-objective Optimizer (GSEMO) [Laumanns et al., 2004b]

• 1: Select a solution x0 uniformly fromX at random; 2: P ← {x0}; 3: while termination condition is not met do 4: Select a solution x from P uniformly at random ; 5: Apply mutation on x to generate y; 6: if ∄z ∈ P satisfies z ≺ y then 7: P ← (P \ {z ∈ P | y ⪯ z}) ∪ {y} 8: end if 9: end while

Theorem 1. For GSEMO with unit-step mutation optimizing L1DM-2SIP, the probability that the population P cannot find the Pareto front is greater than 0.

Proof. Suppose that the initial solution x0 randomly sampled from the solution space belongs to the local Pareto optimal set XL, which happens with some probability larger than 0. Then, the population will converge to the local optimal region XL, from which the algorithm cannot escape by plus/minus one on one dimension (using unit-step mutation). Thus, the Pareto front can never be found.

3.2 Ageing is Provably Helpful

Next, we show that for L1DM-2SIP, using an ageing strategy can help GSEMO escape from local Pareto optima. Ageing has several variants. A common version, i.e., static pure ageing, can automatically restart the algorithm after detecting stagnation in local optima [Horoba et al., 2009; Jansen and Zarges, 2011], and has been shown able to allow (µ+1)- EA to escape from local optima [Oliveto and Sudholt, 2014; Corus et al., 2019]. Here, we adopt the ageing strategy proposed by Li et al. [2025], which assigns each individual an age, increases it by one in each generation, and lets individuals whose age exceeds τ participate in survival selection. This strategy ensures that any newly generated individual can survive for at least τ generations, allowing inferior solutions to remain in the population long enough to reduce the risk of prematurely discarding potentially useful search directions. Meanwhile, the survival selection can keep non-dominated solutions encountered so far.

Algorithm 2 presents GSEMO with the ageing strategy. The age of the initial solution is set to 1. In each iteration, a solution x is randomly selected from the population and mutated. The new solution y has an age of 0. For individuals in the population who have reached age τ , only non-dominated solutions are retained, and among solutions with identical objective values, only one is randomly selected and kept. Individuals with an age less than τ are retained without further filtering. Finally, the new solution y is added to the population, and the age of all solutions is incremented by one.

In Theorem 2, we prove that the expected running time of GSEMO using the ageing strategy and unit-step mutation for finding the Pareto front of L1DM-2SIP can be upper bounded. Note that the running time is measured by the number of fitness evaluations.

Algorithm 2 GSEMO with ageing

• 1: Select a solution x0 ∈ X randomly, and set its age to 1;
• 2: P ← {x0};
• 3: while termination condition is not met do
• 4: Select a solution x from P uniformly at random;
• 5: Apply mutation on x to generate y and set y.age = 0;
• 6: Pτ ← {z ∈ P | z.age ≥ τ};
• 7: Ps = P \ Pτ ;
• 8: Pτ ← Retain the non-dominated solutions with distinct objective vectors in Pτ ;
• 9: P ←Pτ ∪ {y} ∪Ps;
• 10: Increase the age of all solutions in P by 1
• 11: end while

Theorem2. For GSEMOwith ageing and unit-step mutation optimizing L1DM-2SIP, the expected running time for finding the Pareto front is upper bounded by O((c+D(x0))(τ+c)n+ (2b-c+a)τ (τ +c)n(2en(τ +2c+1)/τ +1)2b-c+a ), where x0 is the initial solution, D(x0) = minv∈U∪XL ||x0 - v||1, and τ ≥ 1.

Proof. We divide the optimization process into five phases. The first phase aims to find a boundary solution in U = U1 ∪ U2 ∪ U3 ∪ U4 or a local Pareto optimal solution in XL. The second phase aims to find a boundary solution in U or a local Pareto optimal solution in XLa . The third phase is to find a solution in U + (i.e., the interior region surrounded by U). The fourth phase aims to find a global Pareto optimal solution in X∗ , and the last phase is to find the complete Pareto set X∗ .

According to Lemma 1, we know that the maximum size of the non-dominated solution set generated by GSEMO during the search process is 2c + 1. For GSEMO with ageing, since solutions with age less than τ are always retained, and only non-dominated solutions among those whose age reaches τ are kept, the maximum population size during the search is Pmax = τ + 2c + 1. During the analysis of the first two phases, we pessimistically assume that the solutions in U + are not found; otherwise, we directly go to the fourth phase.

Let disQ(x) = minv∈Q ||x-v||1 denote the distance from a point x in the solution space to a region Q. In the first phase, we consider D(P) = minx∈P disU∪XL (x), i.e., the minimum distance from the population P to the boundary U and the local Pareto optimal set XL. Because a solution will be deleted only if a new solution weakly dominating it appears, and the domination relationship between two solutions is consistent with the relationship between their values of disU∪XL (x), D(P) will not increase. To decrease D(P), it is sufficient to select a solution x′ = argminx∈P disU∪XL (x) and change the value of one of its first two dimensions in the correct direction. For example, if x′ = (c, -b-2, 0, . . . , 0)T , we only need to add one to its second dimension to obtain a solution that is closer to XL, and then D(P) will decrease by 1. For Algorithm 2, the probability of selecting x′ is at least 1/Pmax, the probability of changing the value in the corresponding dimension is (1/n) · (1 - 1/n)n-1 ≥ 1/(en), and the probability of unit-step mutation in the correct direction is 1/2. Thus, D(P) can decrease by 1 in each iteration with probability at least 1/(2enPmax). According to additive drift analysis (Lemma 4 in the Appendix) [He and Yao, 2001], the expected running time of the first phase is at most D(x0) · 2enPmax, as the initial distance D(P) = D(x0).

For the second phase, we assume that only solutions in XL \XLa are found after the first phase; otherwise, the goal of the second phase has reached. As the solutions in XL\XLa can only be weakly dominated by the solutions in U +, they will be always kept in the population. By selecting the solution x′ in XL \XLa which is closest to XLa , and only changing x′1 in the right direction by unit-step mutation, we can move closer to XLa . It is sufficient to perform such steps for c - a times. As such a step happens with probability at least 1/(2enPmax) in each iteration, the expected running time of the second phase is at most (c - a) · (2enPmax).

Now we are in the third phase, where some solutions in U ∪XLa have been found, and the goal is to reach the interior region U +. Consider starting from a solution u ∈ XLa in the current population. First, we select u for mutation and increase only its second dimension value u2, resulting in the offspring solution u′ . The probability of this event occurring is at least 1/(2enPmax). Next, we select u′ from the population and continue increasing its second dimension value. Since u′ is retained in the population for τ iterations, the probability that this event happens at least once in τ iterations is at least p = 1 - (1-1/(2enPmax))τ ≥ 1 - e- τ 2enPmax . Since e-x ≤ 1/(x + 1) for x ≥ -1, we have p ≥ 1 - 1 1+τ/(2enPmax) = τ 2enPmax+τ . This process of selecting a specific solution and increasing its second dimension value needs to be repeated 2b - c + a times to find a solution in U +, so the probability of successfully finding a solution in U + is at least 1 2enPmax · ( τ 2enPmax+τ )2b-c+a . Note that such one successful attempt requires (2b-c+a) · τ +1 iterations. If we consider starting from a boundary solution in U in the third phase, we can directly select it for mutation and change only one of its first two dimensions to enter U +, which happens with probability at least 1/(2enPmax) in one iteration. Thus, the expected running time of the third phase is O((2b - c + a)τnPmax · (2enPmax/τ + 1)2b-c+a ).

In the fourth phase, we have found one solution in U +, and aim to find a Pareto optimal solution in X∗ . For this purpose, we only need to select the solution nearest to X∗ for mutation and change one of its first two dimensions towards the direction of X∗ , which happens with probability at least 1/(2enPmax). After performing such a step for c - a - 1 times, X∗ will be reached. Thus, the expected running time of the fourth phase is at most (c - a - 1) · (2enPmax).

Finally, we are to find the whole Pareto optimal set X∗ . By selecting a proper Pareto optimal solution u ∈ X∗ ∩ P as the parent and only changing u1 in the right direction, we can obtain a new Pareto optimal solution in X∗ \ P. Such a behavior happens with probability at least 1/(2enPmax). As the size of X∗ is 2a+1, the expected running time of the last phase is at most 2a · (2enPmax).

Combining the above five phases, the total expected running time is D(x0)·2enPmax+(c-a)·(2enPmax)+O((2bc + a)τnPmax · (2enPmax/τ + 1)2b-c+a ) + (c - a - 1) · (2enPmax) + 2a · (2enPmax) = O((c + D(x0))(τ + c)n + (2b-c+a)τ (τ +c)n(2en(τ +2c+1)/τ +1)2b-c+a ), where the equality holds by Pmax = τ + 2c + 1.

The benefit of ageing has been demonstrated in [Li et al., 2025], which can accelerate SMS-EMOA in solving the mobjective pseudo-Boolean JUMP benchmark for all even m. Our results in Theorem 2 confirm its advantage in solving MOO problems in integer space. The ageing strategy plays a key role in the analysis of the third phase of the above proof. That is, it enables inferior solutions to survive in the population for a while, and helps the algorithm go across the inferior region to reach U + from the local Pareto optimal set XL.

3.3 Heavy-TailedMutations are Provably Helpful

In this section, we replace the unit-step mutation operator with two heavy-tailed mutation operators [Doerr et al., 2025b], i.e., exponential-tail [Rudolph, 1994] and powerlaw [Doerr et al., 2017] mutations, to allow multi-valued changes in each dimension of a solution, and analyze the expected running time of GSEMO on the L1DM-2SIP problem.

Exponential-tail mutation. The random variable Zi in Eq. (1) satisfies Pr(Zi = k) = α(2 - α)-1 (1 - α)|k| , where α ∈ (0, 1) is a decay parameter.

Power-law mutation. The random variable Zi in Eq. (1) satisfies Pr(Zi = k) = |k|-β /(2ζ(β)), where β ∈ (1, 2) is the power-law exponent and ζ(β) = ∑∞k=1 1/kβ is the Riemann zeta function.

We prove in Theorem 3 that GSEMO can find the whole Pareto front when using heavy-tailed mutations (either exponential-tail or power-law mutation). Unlike the unit-step mutation that only changes the value of each dimension by one, heavy-tailed mutations allow each dimension to change by a larger value. This property enables the algorithm to jump directly from an arbitrary state to a target state, which is crucial for escaping local Pareto optima. It would also be interesting to see whether GSEMO equipped with both ageing and heavy-tailed mutation [Corus et al., 2021] works, which we leave as a future study.

Theorem 3. Consider GSEMO optimizing L1DM-2SIP with the initial solution x0. Let D(x0) = minv∈U∪XL ||x0-v||1. The expected running time for finding the whole Pareto front satisfies the following bounds.

1. If the exponential-tail mutation with parameter α ∈ (0, 1) is used, then the expected running time is

1. If the power-law mutation with exponent β ∈ (1, 2) is used, then the expected running time is

The proof idea is similar to that of Theorem 2, i.e., dividing the optimization procedure into five phases and analyzing the expected running time of each phase. The main difference lies in that the different mutation operators are employed, and the largest population size is reduced from 2c+1+τ to 2c+1 without using the ageing strategy. In particular, in Phases 1 and 4, we apply Lemma 2 to bound the expected running time for reaching the target solution set U ∪ XL or X∗ . This lemma is developed from [Doerr et al., 2025b, Lemmas 14 and 18], with the main difference that their original analysis considers the point (0, 0, . . . , 0) as the search target, whereas in our setting the target is a set of solutions (i.e., U ∪ XL or X∗ ); the largest population size changes from 2a + 1 to 2c + 1. Phase 2 uses heavy-tailed mutation instead of unitstep mutation to find a solution in XLa from a solution in XL \ XLa . In Phase 3, we can use exponential-tail or powerlaw mutation to directly reach the interior region U + from a local Pareto optimal solution x in XLa within a single mutation step, e.g., by increasing x2 from -b to the range of [b - c + a + 1..b + c - a - 1], which happens with probability ∑2b+c-a-1 k=2b-c+a+1 α 2-α (1 - α)k by exponential-tail mutation or ∑2b+c-a-1 k=2b-c+a+1 k-β 2ζ(β) by power-law mutation. Finally, in Phase 5, Lemma 3 (adapted from [Doerr et al., 2025b, Lemmas 15 and 19] by replacing the largest population size 2a+1 in their analysis with 2c + 1 here) is used to bound the expected running time for extending from a single Pareto optimal solution to the entire Pareto front. Due to space limitation, the full proof of Theorem 3 is provided in the Appendix.

Lemma 2 (Adapted from [Doerr et al., 2025b]). Let Q ∈ {U ∪ XL,X∗ } be a target set, and let D(P) := minx∈P, g∈Q ∥x - g∥1 denote the minimum ℓ1-distance between two solutions in the population P and the target set Q. The expected running time of GSEMO for finding a solution in Q starting from P is bounded as follows:

1. If using the exponential-tail mutation with parameter α ∈ (0, 1), it is upper bounded by ((2c + 1)en/K)(nπ2 /(6α)+4αD(P )), where K is a constant.
2. If using the Power-law mutation with exponent β ∈ (1, 2), it is upper bounded by (2c + 1) · 2enζ(β)(21/(2-β) + 2(2 - β)(β - 1)-1 (D(P ))β-1 ).

Lemma 3 (Adapted from [Doerr et al., 2025b]). If a population P contains a Pareto optimal solution in X∗ , then the expected running time of GSEMO for finding the whole Pareto front F ∗ is bounded as follows:

1. If using the exponential-tail mutation with parameter α ∈ (0, 1), it is upper bounded by C(cnmax{ln(a + 1)/(aα), aα+ ln(a+1)}), where C is a large constant.
2. If using the Power-law mutation with exponent β ∈ (1, 2), it is upper bounded by 4 ln(2)enζ(β)(β - 1)(1 - 1.51-β )-1 (2a + 1)β (1 - 21-β )-2 + (2c + 1) · 2en ζ(β)(ln(a + 1) + 1).

4 Extension of L1DM-2SIP

In this section, we give a more general form of the benchmark problem L1DM-2SIP proposed in Section 2, which includes complex local optima regions.

Definition 6 (L1DM-2SIP). Consider the problem of L1distance minimization to two sets of integer points: Zn → N20 with

where hr(x) = mini∈[1..kz](∥x - zr,i∥1) and er(x) = minj∈[1..kv](∥x - vr,j∥1). Setting u1 = (-a, b, 0, ..., 0)T, u2 = (a, b, 0, ..., 0)T, z1,i = ((-1)i cγ(i), (-1)i bγ(i) + b, 0, . . . , 0)T, z2,i = ((-1)i+1 cγ(i), (-1)i bγ(i) + b, 0, . . . , 0)T, v1,j = ((-1)j eγ(j), (-1)j gγ(j) + b, 0, . . . , 0)T, and v2,j = ((-1)j eγ(j), (-1)j+1 gγ(j) + b, 0, . . . , 0)T, where γ(i) = 2 · ⌈i/2⌉-1, ck, bk, ek, and gk all increase gradually with k ≥ 1, and a < c1 < e1.

When c1 = c, b1 = 2b, and L1DM-2SIP depends only on the four points u1, u2, z1,1, and z2,1, the problem reduces to the instance considered in Sections 2 and 3. For the general form of L1DM-2SIP, under conditions of b1 > a + c1, e1 - a > g1 + 1, and for k ≥ 1, gγ(k) = cγ(k), ck + 2|bk+1 - bk| - ck+1 > 1, and gk + 2|ek - ek+1| - gk+1 > 1, L1DM-2SIP has the Pareto set X∗ = {(x1, b, 0, . . . , 0)T ∈ Zn | x1 ∈ [-a..a]}, and a series of local Pareto optimal solution sets XLh,i and XLe,j . The horizontal XLh,i = {(x1, (-1)i bγ(i) + b, 0, . . . , 0)T ∈ Zn | x1 ∈ [-cγ(i)..cγ(i)]}, and the vertical XLe,j consists of all ((-1)j eγ(j), x2, 0, . . . , 0)T ∈ Zn with x2 ∈ [-gγ(j) + b..gγ(j) + b]. According to Definition 2, the solution space may still contain a small number of discrete local Pareto optimal solutions, but these solutions will not form peaks and thus will not cause the algorithm to get stuck. Figure 2 gives an example of such a problem.

This problem can be instantiated to exhibit a more "interesting" structure by gradually increasing the distance between the point pairs from the global optimum towards outer regions, thereby forming a landscape with many local optima of progressively improving quality as they approach the global optimum. Such a landscape resembles those of many wellestablished functions like ZDT4 [Zitzler et al., 2000], DTLZ1 and DTLZ3 [Deb et al., 2005].

5 Experiments

In this section, we conduct experiments on the proposed integer-space function to complement the theoretical results. Table 1 presents the results of GSEMO solving L1DM2SIP with different number of fitness evaluations FEs ∈

{1×103 , 5×103 , 1×104 , 5×104 , 1×105 , 5×105 , 1×106 }, problem size n = 2, a = 3, b = 3, c = 4, and each dimension initialized within the range [-10, 10]. We present the average PF coverage and standard deviation over 50 runs for four variants of GSEMO: 1) GSEMO with unit-step mutation, 2) GSEMO with ageing (τ = 4) and unit-step mutation, 3) GSEMO with exponential-tail mutation (α = 0.5), and 4) GSEMO with power-law mutation (β = 1.5). The results show that increasing the number of fitness evaluations (FEs) does not improve the performance of GSEMO with unit-step mutation, indicating that the algorithmmay have encountered a bottleneck, specifically, getting stuck in local Pareto optima. Introducing an ageing strategy improved this situation, increasing the coverage of the PF at the same number of function evaluations. To further examine the effect of τ , we added results with τ = 2, 6, 8 in the Appendix, suggesting that τ cannot be set too small in our studied case. However, even with the ageing strategy, the results with unit-step mutation were still inferior to those with heavy-tailed mutation, because the latter allows for mutation strengths greater than 1.

6 Conclusion

In this paper, we take the first step towards the running time analysis of MOEAs on integer spaces with local optima. We begin by presenting a multi-objective benchmark defined on an integer space, characterized by an analyzable landscape and the presence of local optima. Next, we conduct a running time analysis on this benchmark and derive some theoretical results. Specifically, we prove that the widelystudied GSEMO with unit-step mutation may fail to find the Pareto front, whereas introducing an ageing mechanism or employing heavy-tailed (exponential-tail or power-law) mutation operators can enable GSEMO to escape local optima and identify the Pareto front within a bounded expected running time. Empirical studies also support our theoretical results. As a class of MOO problems with practical relevance, the proposed benchmark narrows the gap between the theoretical analyses and practice. We hope this work will inspire further research along this direction, e.g., considering other well-established MOEAs or the extension of L1DM-2SIP to broaden the theoretical understanding of search algorithms on integer spaces.

Acknowledgments

This work was supported by the National Science and Technology Major Project (2022ZD0116600), the NSFC (62276124), the Fundamental Research Funds for the Central Universities (14380020), the Fundamental and Interdisciplinary Disciplines Breakthrough Plan of the Ministry of Education of China (JYB2025XDXM118), and the "111 Center" (B26023). Chao Qian is the corresponding author.