Stochastic Pareto local search : Pareto neighbourhood exploration and perturbation strategies

Pareto local search (PLS) methods are local search algorithms for multi-objective combinatorial optimization problems based on the Pareto dominance criterion. PLS explores the Pareto neighbourhood of a set of non-dominated solutions until it reaches a local optimal Pareto front. In this paper, we discuss and analyse three different Pareto neighbourhood exploration strategies: best, ﬁrst, and neutral improvement. Furthermore, we introduce a deactivation mechanism that restarts PLS from an archive of solutions rather than from a single solution in order to avoid the exploration of already explored regions. To escape from a local optimal so-lution set we apply stochastic perturbation strategies, leading to stochastic Pareto local search algorithms (SPLS). We consider two perturbation strategies: mutation and path-guided mutation. While the former is unbiased, the latter is biased towards preserving common substructures between 2 solutions. We apply SPLS on a set of large, correlated bi-objective quadratic assignment problems (bQAPs) and observe that SPLS signiﬁcantly outperforms multi-start PLS. We investigate the reason of this performance gain by studying the ﬁtness landscape structure of the bQAPs using random walks. The best performing method uses the stochastic perturbation algorithms, the ﬁrst improvement Pareto neigborhood exploration and the deactivation technique.


Introduction
Pareto local search algorithms (PLS) (Paquete et al. 2004;Angel et al. 2004;Basseur 2006) are local-or neighbourhood-search techniques for solving multiobjective combinatorial optimization problems.PLS applies an exploration strategy to iteratively move from a set of current solutions to a neighbouring set that improves upon the current solution.The algorithm stops in a Pareto local optimal set: a set of solutions that have no improving solutions in their neighbourhood.Because the local search gets stuck in a local optimal set, the search needs to be restarted from different regions of the search space.Multi-start PLS restarts the search from randomly generated points in the search space.However, this does not exploit information present in previously found solutions.
In the literature some variants of PLS algorithms can be found, each having small differences in their neighbourhood exploration technique.Paquete et al.'s PLS (Paquete et al. 2004) stores only those solutions in the archive that are Pareto optimal.Basseur's PLS (Basseur 2006) explores all solutions in the archive exactly once, even though they might have become dominated by newly included solutions.Both PLSs explore the entire neighbourhood of solutions in the archive.Thus, these PLSs use a best improvement exploration strategy (Hansen and Mladenović 2006).Angel et al.'s PLS (Angel et al. 2004) uses dynasearch and dynamic programming to explore the neighbourhood of the bi-objective travelling salesman problem.Liefooghe et al. (2011) experimentally compare various instances of PLS algorithms using different parameters and settings on bi-objective instances of travelling salesman and scheduling problems.
Variations of PLS were proposed to improve PLS performance (Hamacher and Ruhe 1994).Two phase PLS (Lust and Teghem 2010;Dubois-Lacoste et al. 2011a) first approximates the Pareto local optimal set with solutions generated with singleobjective solvers.In the second phase, the algorithm uses PLS (Paquete et al. 2004;Angel et al. 2004) to generate solutions not found in the first phase.Guided PLS (Alsheddy and Tsang 2009) uses a variant of PLS that prioritizes the objectives of the search space to uniformly spread the solutions over the Pareto local optimal set.Guided PLS is somehow the inverse of the two phase PLS because it first performs PLS and then tries to escape from a Pareto local optimal set by adjusting the weights of different objectives.
Iterated PLS (Drugan and Thierens 2010) uses stochastic perturbation operators to escape from a Pareto local optimal set.Iterated PLS is the multi-objective equivalent of iterated local search (Lourenco et al. 2003) a technique for single objective optimization that restarts LS from solutions generated by mutating the currently found local optima.In addition to mutation SPLS can generate restarting solutions using recombination operators because there are multiple solutions in the Pareto local optimal set that can be used as parent solutions.The Pareto local optimal set can be used directly as the population or alternatively, a separate elitist population can be used to generate new solutions (Lopez-Ibanez et al. 2006).

Main contributions of the paper
We propose several techniques to ameliorate the performance of Pareto local search algorithms.Our variation of PLS uses a technique-called deactivation-that restarts PLS from a newly generated solution and a set of incomparable solutions that were generated with previously restarted PLS runs.In Sect.2, we show that PLS using deactivation can avoid the repetitive exploration of (non-promising) regions.
PLS algorithms usually explore the entire neighbourhood of a solution (Paquete et al. 2004).We show here though that first improvement neighbourhood exploration strategies can be more efficient.A similar observation has recently been made in (Liefooghe et al. 2011).Pareto neighbourhood exploration strategies-or improvement strategies-are the equivalent of single-objective improvement strategies (Hansen and Mladenović 2006) for multi-objective spaces.We view Pareto improvement as a relationship between an initial solution, a solution from the neighbourhood of the initial solution, and all the other solutions of the current Pareto set.In Sect.3, we design and analyse three Pareto neighbourhood exploration strategies.We show that these Pareto improvement strategies are well performing.We prove that PLS using our Pareto neighbourhood exploration strategies are proper, sound and complete.We also discuss the relationship between our Pareto improvement strategies and another, often used, technique to search in multiobjective search spaces, namely the hypervolume unary indicator (Zitzler et al. 2003;Bringmann and Friedrich 2010).The (Pareto) improvement strategies are independent of the landscape they are applied on, thus these results can straightforwardly be applied to other combinatorial problems.
To escape from a local optimal solution set we apply stochastic perturbation strategies, leading to stochastic Pareto local search algorithms (SPLS).We consider two perturbation strategies: mutation and path-guided mutation.While the former is unbiased, the latter is biased towards preserving common substructures between two solutions.In Sect. 4 we discuss a SPLS that uses mutation and path-guided mutation as perturbation operators.Since path-guided mutation can be seen as a specific form of recombination we refer to this algorithm as Genetic PLS (GPLS).The structure of the neighbourhood and the perturbation operators depend on the representation and particularities of the problem they are tested on.In Sect.5, we experimentally investigate the efficiency of the proposed mechanisms (i.e., the Pareto improvement strategies, deactivation, genetic operators) for PLSs on instances of the bi-objective QAP (bQAP) (Knowles and Corne 2003).The solutions of bQAPs are represented by permutations of facilities to different locations.
The landscape analysis in Sect.6 records the quality of the restarting solutions with random walks (RW).We extend the scope of the RW method to study the relationship between the distance to a Pareto local optimal set and quality measures, like the hypervolume unary indicator, or mechanisms like deactivation.RW is a logical choice for analysing the behaviour of GPLS in relationship with the genetic operators because the solutions generated with these operators are all accepted as restarting points for PLS (Merz and Freisleben 2000).RW reveals correlations between the structure of the search space and the genetic operator used.Section 7 concludes the paper.

Stochastic Pareto local search (SPLS)
We define stochastic Pareto local search (SPLS) algorithms as a combination of stochastic perturbation operators and Pareto local search (PLS).In Sect.2.1, we first introduce a slight generalization of PLS (Paquete et al. 2004;Angel et al. 2004;Basseur 2006) to allow the design of more efficient SPLS algorithms.In Sect.2.3, we present the genetic PLS algorithms that use genetic operators to restart PLS.

Notation and definitions
Let us consider the multi-objective function f = (f 1 , . . ., f m ), where f : S → O. S is the set of solutions in a countable solution space.The image of the solution set S under f is denoted with the objective space O = f(S) = {y | ∃x ∈ S : y = f(x)}.We consider the single objective space as a special case of the multi-objective space with m = 1.A bi-objective space has m = 2.
Let a, b ∈ S be two solutions and f(a), f(b) ∈ O be their images in objective space, called objective vectors.To define an optimization problem, we first consider an order binary relationship, , between any two elements in the objective space O. Without loss of generality, only minimization is considered.In Table 1, we present the relationships between objective vectors and sets of objective vectors used in this paper (Zitzler et al. 2003).
A Pareto set is a set of non-dominated solutions.The set of Pareto optimal solutions, or Pareto optimal set, A ⊂ S, is Min(S, ≺) = {a ∈ A ⊂ S | ¬∃x ∈ S : f(x) ≺ f(a)}.A Pareto local optimal set is the set of non-dominated solutions A alg generated with a search algorithm alg, Min(S alg , ≺) = {a ∈ A alg ⊂ S alg ⊂ S | ¬∃x ∈ S alg : f(x) ≺ f(a)}, where S alg ⊂ S is the subset of solutions explored with the search algorithm.A maximal Pareto local set is a Pareto local optimal set containing solutions for which the solutions in their neighbourhood are either dominated or in the Pareto local optimal set.

Pareto local search
Algorithm 1 presents the pseudo-code of our PLS algorithm.For ease of discussion, we consider an archive of unlimited size that can store all the solutions in the Pareto set.For very large Pareto sets, for instance for real coded problems where unlimited  size Pareto sets are unrealistic, different archiving methods can be applied (Knowles and Corne 2004;Deb et al. 2005;Lopez-Ibanez et al. 2011).
PLS(A, I) has two input parameters.A is a Pareto set containing at least one solution s with visited flag set to false.The Pareto improvement strategy, I(s, A \ {s}, N , , f), explores the neighbourhood N of a solution s using the Pareto set A \ {s} and the domination relationship over the function f.The output of I(•) is a Pareto set, A .The neighbourhood function N : S → P(S) has as input a solution from S and returns the set of neighbours for that solution.N (s) depends on the problem (e.g., multi-objective QAPs) and the solution s it is applied on.I's algorithm does not depend on the landscape; it either exhaustively explores the entire neighbourhood and selects all non-dominated solutions in it, or stops when the first improvement to the initial set of solutions is found.
Each iteration of PLS, a solution s with the visited flag set to false, s.visited = false, is randomly chosen from the current Pareto set (or archive) A. The solutions in the neighbourhood of s are evaluated with I.The function merge merges M ≥ 2 Pareto sets into a new Pareto set, where Note that merge also removes the dominated solutions from the input Pareto set.The newly added solutions from the neighbourhood have the visited flag set on "false".The search continues until there are no unvisited solutions in A .Thus, the search stops when there are no "new" incomparable or dominating solutions found.Proposition 1 Let PLS be as in Algorithm 1, with I the best Pareto improvement strategy.Consider two initial Pareto sets A 1 = {s} and A 2 = {s, q 1 , . . .q k } for PLS.Consider the same sequence π N (s) of visiting the solutions in N (s).Then, the number of solutions from N (s) visited with PLS(A 1 , I) is smaller or equal with the number of solutions N (s) visited with PLS(A 2 , I).The set of solutions in N (s) visited only with PLS(A 1 , I) are dominated by at least one solution from {q 1 , . . .q k }.

Deactivation technique
Proof Without loss of generality, let's assume that k = 1.Then, A 2 = {s, q 1 }.Consider the sequence of solutions {s 1 , . . ., s n } ∈ π N (s) and these solutions' relationships with s and q 1 .The only case when the two PLSs are behaving differently is when ∃s i that is dominated by q 1 and incomparable with s.With PLS(A 1 , I), s i is added to A .With PLS(A 2 , I), s i is not added to A .We conclude that there are less solutions from π N (s) that are evaluated with PLS using A 2 than when A 1 is used, and the extra solutions evaluated with A 1 are dominated by q 1 .
Example 1 shows situations where the deactivation techniques is useful.
Example 1 In Fig. 1, a solution is denoted with a point.A circle or a part of a circle denotes the neighbourhood of a solution.A Pareto set is denoted with a curve.In Fig. 1, the initial solution s has in its neighbourhood N (s) two incomparable solutions s 1 and s 2 , f(s 1 ) f(s 2 ), f(s) f(s 2 ) and f(s) f(s 1 ), where s 1 , s 2 ∈ N (s).Let s 3 be in the neighbourhood of s 1 , s 3 ∈ N (s 1 ), such that s 3 dominates s, s 1 and s 2 , f(s 3 ) ≺ f(s) and f(s 3 ) ≺ f(s 1 ) and f(s 3 ) ≺ f(s 2 ).Let s 4 and s 5 be in the neighbourhood of s 2 , s 4 , s 5 ∈ N (s 2 ), such that s 4 dominates all other solutions, and s 5 is incomparable with s 3 and dominates s 1 , s 2 and s.
Let's consider PLS, PLS(A, I), from Algorithm 1 with I the best Pareto improvement strategy, and different Pareto sets A.
In Fig. 1(a), we consider A = {s} and two equal probable sequences to visit the neighbourhood of s, {s 1 , s 2 } and {s 2 , s 1 }.When the neighbours of s are visited in the A ← ∅ while Stopping criterion, T , is NOT met do Generate s uniform randomly A ← Deactiv(s, A) A ← merge(A, PLS(A , I)) end while return A sequence {s 1 , s 2 }, the PLS's output is A = {s 3 } because s 3 ∈ N (s 1 ) is dominating s 2 that is deleted from the archive A .When the neighbours of s are visited in the sequence {s 2 , s 1 }, the PLS's output is A = {s 4 } because s 3 ∈ N (s 1 ) is added to the archive A after adding s 4 ∈ N (s 2 ) and s 4 dominates s 3 .
In Fig. 1(b), q and q are two solutions that are incomparable with s, and q dominates s 2 and is incomparable with s 1 , and q dominates s 1 and is incomparable with s 2 .When A = {s, q}, the PLS's output is A = {s 3 } because s 2 is dominated by q and is deleted from the archive A .When A = {s, q }, then s 1 is dominated by q and is deleted from the archive A , and the output is A = {s 4 }.
Note that both PLSs can delete s 3 whose neighbourhood contains s 4 that dominates all the other solutions.In the first case, the output of PLS depends on the sequence in which solutions in the neighbourhood of s are evaluated.In the second case, the output of PLS depends on the initial Pareto set.
In Algorithms 3 and 4, the input Pareto set A of PLS is equal to the output of the deactivation method.

The multi-restart PLS algorithm (MPLS)
MPLS is a straightforward algorithm to escape from a Pareto local optimal set.It restarts multiple times from uniform randomly chosen initial solutions.Algorithm 3 presents the pseudo-code for MPLS.The input parameters for MPLS are a Pareto Algorithm 4 Genetic PLS GPLS(I, α, T 1 , T ) A ← MPLS(I, T 1 ) while Stopping criterion, T , is NOT met do Select s uniform randomly from A if α > U(0, 1) or |A| < 2 then s ← Mutation(s) else Select s 1 = s from A s ← Recombine(s, s 1 ) end if A ← Deactiv(s , A) A ← merge(A, PLS(A , I)) end while return A improvement strategy I and a stopping criterion T .The stopping criterion can be, for example, a maximum number of restarted PLSs.The Pareto set A is initialized as the empty set, A ← ∅.Until T is met, a solution s is uniform randomly sampled from the search space.The dominating solutions in A by s are deactivated, A ← Deactiv(s, A).The output Pareto set A is the input Pareto set for PLS, PLS(A , I).The Pareto set A is updated by merging it with the output of PLS.
The difference between this multi-start PLS and Paquete's multi-start PLS is: (i) the use of the deactivation technique, and (ii) the use of different Pareto improvement strategies I. Paquete et al.'s PLS starts to explore the neighborhood of a solution with an empty archive.MLS cannot really profit from an exploration archive, because the restarting solution s is uniform randomly generated, it is improbable that there are solutions in A beside s.That is, all the solutions in the Pareto local optimum set A are probably dominating a randomly generated solution s.The deactivation technique will lead to an empty archive.This hypothesis is confirmed in Sect.6, where MPLS from Algorithm 3 and Paquete's MPLS have similar performance.Therefore, deactivation is not useful unless fit restarting solutions are generated, for example using genetic operators on solutions from the Pareto local optimum set A. We show that the use of this exploration archive is very useful for GPLS.

Restarting PLS using genetic operators
To improve upon multi-restart PLS's performance, stochastic Pareto local search aims to escape from Pareto local optimal sets by stochastic perturbation operators that preserve partial information of the perturbed solutions.In case of genetic PLS (GPLS), the perturbation operators are a combination of mutation and recombination.Algorithm 4 gives the pseudo-code of GPLS.
First, as initialization phase, an MPLS runs in order to construct an initial Pareto set A with diverse solutions from different basins of attraction.The stopping criterion T 1 is the number of times N > 0 MPLS is restarted from uniform randomly generated solutions.The type of operator, i.e. mutation or recombination, is selected proportionally with a fixed probability α ∈ [0, 1].A new individual s is generated with mutation, Mutation(s), or with recombination, Recombination(s, s 1 ), where s 1 is another individual from A, s = s 1 .The solutions in A are deactivated with Deactiv(s , A) and the resulting Pareto set A is the initial Pareto set for PLS.PLS is restarted from s .The archive A is updated with the output of PLS, A .This process is repeated until a stopping criterion, T , is met.T depends on the problem, e.g., mQAPs, and is presented in the next section.The GPLS algorithm returns a Pareto local optimal set.
When α = 1 the individuals are generated only with mutation.The advantage of mutation-generated restarts is that the local search is restarted from nearby areas of the landscape.When α = 0.5 the individuals are generated 50 % with mutation and 50 % with recombination.Recombination generates solutions at larger distance from the current solution but still aims to preserve partial information of the two parent solutions.
Note that the deactivation method does not have any effect if there are no solutions in the Pareto set A that are incomparable with the solution s.We assume that the individuals generated with mutation or recombination from solutions in the current Pareto set are "good" solutions that are incomparable with some solutions in the Pareto set.

Pareto neighbourhood exploration strategies
The performance of LS algorithms depends on the choice of improvement strategy (Hansen and Mladenović 2006;Liefooghe et al. 2011).Until recently, best improvement algorithms were commonly used for the neighbourhood exploration strategies in Pareto local search algorithms.The first non-dominating neighbour (Liefooghe et al. 2011) stops after the first non-dominating solution found.The first dominating neighbour (Aguirre and Tanaka 2005;Liefooghe et al. 2011;Dubois-Lacoste et al. 2011b) stops after the first dominating solution.
In this section, we introduce a new concept for neighbourhood exploration strategies that compares the solutions in a neighbourhood N (s) not only to s but also to the solutions from the Pareto set A. In Algorithm 1, A is a Pareto set of solutions incomparable with s that were previously generated in a PLS run.We show that usage of A can: (i) improve the quality of the outputted Pareto set, and (ii) return a smaller Pareto set.We compare our Pareto improvement strategies with another very important method to search through the multi-objective space, the hypervolume unary indicator.
Because we consider one type of neighbourhood function, N , one type of order relationship , and one function f, we will often use I(s, A) instead of I(s, A, N , , f) to shorten the notation.

Best Pareto improvement, I B
Like the best improvement for single objective spaces, the best Pareto improvement strategy explores all the solutions in the neighbourhood of the input solution s.The pseudo-code for I B is given in Algorithm 5.
The best Pareto improvement strategy from Algorithm 5, I B (s, A), evaluates all the individuals in the neighbourhood of s, N (s).The outputted set A is empty if there are no solutions in N (s) that are incomparable or dominating all solutions in {s} ∪ A. Otherwise, A contains solutions from N (s) that are dominating or incomparable with {s} ∪ A and all other solutions from N (s).The visited flag of newly added solutions in A is set to false.
Note that the output A does not contain solutions from the initial set {s} ∪ A. This initial set is eliminated in the last line of the algorithm I B , A ← A \ ({s} ∪ A).In the next proposition, we show that the difference between the outputs of the standard best improvement strategy, I B (s, ∅), and I B (s, A) is the size of the outputted Pareto set A .The output of I B (s, ∅) could be larger than I B (s, A), where the differences are solutions that are dominated by at least one solution from A. We also show that merging the resulting Pareto set of the two I B with A leads to the same result.
Proposition 2 Let I B be as in Definition 2.Then, between the Pareto sets I B (s, A) and I B (s, ∅) are the following relationships: Proof The first relationship holds because there can be solutions in N (s) that are incomparable with s but dominated by a solution from A. Then, the output of I B (s, ∅) contains such a solution, but I B (s, A) does not include it.In all other cases, the output of the two improvement strategies are the same.The second relationship, the equality for the best Pareto improvement, holds because all the solutions in N are explored regardless of the initial Pareto set.

Neutral Pareto improvement strategy, I N
The pseudo-code for the neutral Pareto improvement I N is presented in Algorithm 6.I N stops the exploration of a neighbourhood N (s) when the first solution that is non-dominated by the initial solution s and non-dominated by the Pareto set A is found.
Return A end if end for return ∅ Definition 3 Let I N (s, A) be the neutral Pareto improvement strategy from Algorithm 6 and π N (s) denote a random permutation of N (s).I N stops after the first solution a ∈ π N is found such that a is incomparable with s and incomparable with all solutions from A. If a exists, the outputted Pareto set is A ← {a}, and a.visited ← false.Otherwise, A ← ∅.
Again, the initial solutions {s} ∪ A are not contained in the output A of I N .In the following proposition, we show that, even after merging the output of I N with A, I N (s, A) is at least as good-meaning equal or dominating-with I N (s, ∅).
Proposition 3 Let I N (s, ∅) be the neutral Pareto improvement strategy from Definition 3.Then, the outputted Pareto sets of I N (s, A) and I N (s, ∅) are in the following relationships Proof Let s and s ∈ π N be two solutions such that: (i) s is the first and s is the second solution in π N (s) that are incomparable with s, and (ii) s is dominated by at least a solution in A, and (iii) s is incomparable or dominating all solutions in A. Then, I N (s, ∅) stops at s and I N (s, A) selects s .For the rest of the cases, the output of I N (s, A) and I N (s, ∅) are equal.We conclude that the output of I N (s, A) is nondominated by the output of I N (s, ∅).Otherwise, the output of This means that there are cases where neutral Pareto improvement is outperforming first non-dominated solution (Liefooghe et al. 2011) but there are no situations where I N is outperformed by the first non-dominated solution.

First Pareto improvement strategy, I F
The pseudo-code for I F is given in Algorithm 7. In I F , the exploration of the neighbourhood stops when the first solution that dominates the initial solution is found.s, A) be the first Pareto improvement strategy from Algorithm 7 and let π N (s) denote a random permutation of N (s).I F stops at the first solution a ∈ N (s) that dominates s.If ∃a, the output A contains a and all the solutions from π N (s) that precede a and that are: (i) incomparable with a and (ii) dominating or are incomparable with all solutions from {s} ∪ A and all the other solutions from π N (s) that precede a. Otherwise, if ¬∃a, then A contains all the solutions from N (s) that are: (i) incomparable with s, (ii) incomparable or dominate all solutions in A, and (iii) incomparable or dominating all other solutions in N (s).The visited flag of all solutions included in A is set on false.
Like for the other two improvement strategies, from the output of I F are deleted the initial set of solutions {s} ∪ A. Similarly with I B , the Pareto set of I F (s, A) is smaller or equal with the output of I F (s, ∅), where the extra solutions are dominated by A. After merging the outputs of I F with A, the resulting Pareto sets are equal.
Proposition 4 Let I F be as in Definition 4.Then, we have Proof Let us assume that ∃a ∈ N (s) such that f(a) ≺ f(s).Then, a is non-dominated by the currently explored solutions for π N (s) and all solutions from A. Compared with I F (s, A), I F (s, ∅) returns solutions that are incomparable with a but dominated by a solution from A. Assume that ¬∃a ∈ N (s), f(a) ≺ f(s), then all the solutions in N (s) are evaluated regardless of the initial Pareto set.Again, both strategies return the same result.
The first improvement strategy (Liefooghe et al. 2011) is slightly different from I F (s, ∅).With the first improvement all non-dominated solutions with s are proposed for integration with A that includes solutions dominated by a or by solutions in A.
With I F (s, ∅), the solutions that are dominated by a are already deleted from the output.This makes a difference in computational time for large neighbourhoods with lots of incomparable solutions.

Pareto improvement strategies vs searching with the hypervolume unary indicator
The use of the initial Pareto set A allows us to compare the Pareto improvement strategies with a single-objective first improvement exploration strategy using the hypervolume unary indicator, a technique currently often used for multi-objective optimization (Zitzler et al. 2003;Bringmann and Friedrich 2010).This observation is important because some popular algorithms use the hypervolume indicator.For example, Emmerich et al. (2005), Beume et al. (2007), Bader and Zitzler (2011) use the hypervolume to select a subset of solutions from a larger set or to discard a solution from a set.The hypervolume unary indicator function, H, measures in the objective space the hypervolume contained between a reference point, r, and the Pareto set A. The hypervolume H(r; A) is maximized, meaning that the larger the hypervolume is, the "better" its set A is considered.For a given reference point the hypervolume indicator of A is by design equal or smaller than the hypervolume indicator of ({s} ∪ A) (Bader 2010).Algorithm 8 gives the pseudo-code for the hypervolume-based first improvement strategy.
Definition 5 Let the hypervolume H(r; A) act on the reference point r and a Pareto set A. The hypervolume-based first improvement strategy from Algorithm 8, I H (s, A, N , ≥, H), evaluates the hypervolume indicators of each of the sets merge({s }, {s} ∪ A), where s ∈ π N , a random permutation of N (s).I H stops after the first solution a ∈ π N for which H(r; A) < H(r; merge({a}, A)).If ∃a, the output is A ← {a}.Otherwise, A ← ∅.
We show the equivalence of I N with I H for any reference point r.Then, Proof We split the proof in two parts.First, we show that I N ⊆ I H , and then we show that I H ⊆ I N .I N stops at the first generated solution s 1 ∈ N (s) for which f(s 1 ) f(s), and ¬∃s ∈ A, f(s for which the hypervolume indicator of H(r; merge({s 1 }, {s} ∩ A)) is larger than H(r; {s} ∩ A).Then, I H selects s 1 and stops.I H stops at the first generated solution s 2 ∈ N (s) such that H(r; {s} ∪ A) < H(r; merge({s 2 }, {s} ∪ A)).This means that {s 2 } is either: (i) incomparable with s and all other points in A, (ii) s 2 dominates s.In both cases, I N selects s 2 .We conclude that To summarize, the output of local search that uses the hypervolume-based first improvement indicator is equivalent with the output of PLS with the neutral Pareto improvement strategy, but the computational cost of using the hypervolume indicator is much higher than the cost of running PLS with the neutral Pareto neighbourhood exploration.

Comparing Pareto neighbourhood exploration strategies
The behaviour of the three Pareto neighbourhood exploration strategies-I B , I F , and I N -are compared by means of an example.
Example 2 In Fig. 2, we denote the solutions with points and crosses and Pareto sets with line segments.Consider five solutions in N (s): (i) s 1 and s 2 are dominating s and are incomparable with each other, (ii) s 3 dominates s but is dominated by s 1 and s 2 , (iii) s 4 is incomparable with s but dominated by s 1 , s 2 and s 3 , (iv) s 5 is dominated by all other solutions.In Fig. 2(b), additionally to the solutions from Fig. 2(a), the initial Pareto set is A = {p}, where p is incomparable with both s 1 and s 2 .
Both I B (s, ∅) in Fig. 2(a) and I B (s, {p}) in Fig. 2(b) return the dominating solutions s 1 and s 2 , A = {s 1 , s 2 }.I F (s, ∅) and I F (s, {p}) select one of the solutions that dominates s: s 1 , s 2 or s 3 .Then, the output is A = {s 1 }, or A = {s 2 }, or A = {s 3 }.I N (s, ∅) selects one of the four solutions that dominates or is incomparable with s: s 1 , s 2 , s 3 , and s 4 .I N (s, {p}) selects one of the three solutions that dominates or is incomparable with {s, p}: s 1 , s 2 , and s 3 .In Fig. 2(b), s 4 is not selected because it is dominated by solutions in A. Then, I N (s, ∅) is dominated by I N (s, {p}).s 5 is never selected.
To study the differences between the three Pareto neighbourhood exploration strategies I B , I N and I F , we identify the conditions under which these strategies return the same result.The following statements result directly from their definitions.Let π N be a permutation of N (s) and let |I(s, A)| be the number of solutions evaluated in the algorithm I(s, A).The output of I B (s, A) and I F (s, A) are the same Pareto set if and only if one of the following conditions is met: 1. there exists no solution s ∈ N (s) such that s dominates s, 2. there exists exactly one solution s ∈ N (s) that: (i) dominates s, (ii) dominates all solutions in N (s) \ {s }, and (iii) precedes in π N (s) , if exists, the other solutions that dominate s but are dominated by s .
The output of I N (s, A) is the same as the output of I F (s, A) if and only if one of the following conditions hold: 3. there are in N (s) only solutions that are dominated by at least a solution in {s}∪A, 4. there exists exactly one solution s ∈ N (s) that: (i) is incomparable or dominating all solutions in {s} ∩ A, and (ii) all other solutions in N (s) \ {s } are dominated by at least a solution in {s} ∩ {A} The output of I N (s, A) is the same as the output of I B (s, A) if and only if the above conditions 3 and 4 hold.It is interesting to note that the conditions under which the outputs of I F and I B are equal include the conditions for which the outputs of I F and I N are equal, which, at their turn are equal with the conditions for I B and I N to be equal.
In the following, we compare the number and the quality of the solutions returned when the three Pareto neighbourhood exploration strategies are used.According to the above proposition, there is a trade-off between the number of evaluated solutions and the quality of the newly added solutions.The advantage of I F and I N is that they evaluate only a part of the neighbourhood, whereas I B evaluates the entire neighbourhood.The disadvantage of I F and I N is that the solutions selected, in general, are dominated by the solutions selected with I B .The solution(s) generated with I F are non-dominated by the solutions selected with I N .Note that the time complexity of evaluating a solution with any of the three Pareto neighbourhood exploration strategies presented before is linear with the number of objectives, m, and the number of solutions in the neighbourhood of a solution s, N (s), and the number of solutions in the initial Pareto set A.
3.6 Soundness and completeness of PLS using Pareto improvement strategies Paquete et al. (2007) show that their PLS is always stopping and, it returns a maximal Pareto local optimal set.Our version of PLS from Algorithm 1 differs in certain aspects from Paquete's PLS: (i) the multiple Pareto neighbourhood exploration strategies , I, and (ii) the starting Pareto set A. In Paquete's PLS, the Pareto set A is reduced to a single, randomly generated solution and there is a single improvement strategy, the best (Pareto) improvement strategy.We show that PLS using one of the three Pareto improvement strategies and an initial Pareto set has the local optimality property-that is, it always stops in a Pareto local optimal set.For simplicity of the proof, we assume that the initial Pareto set A contains only solution, s, with the visited flag set on true.This proof follows the proof of Theorems 5.3 and 5.4 from (Paquete et al. 2007), but with different improvement strategies for PLS.
Theorem 1 In Algorithm 1, PLS(A, I B ) and PLS(A, I F ) always stop in a maximal Pareto local optimal set.PLS(A, I N ) stops in a subset of a maximal Pareto local optimal set.Proof In the first part of this proof, we show that no matter in which sequence solutions are visited, the output contains always a set of solutions that dominates the other solutions.Let us assume that there are two solutions, s and s in A , such that f(s ) ≺ f(s ) and both solutions are evaluated with the same PLS.If s is visited first, it is added to A .When s is visited, s is discarded from A .If s is visited first, it is added to A .When s is visited, it is not added to the archive because s is dominated by s .We conclude that s is always added to A regardless of the succession in which solutions are visited.
Second, we show that PLS using I B and I F are always stopping in a maximal Pareto local optimal set.Thus, there are no solutions in the neighbourhoods of the solutions from the Pareto local optimal set that are incomparable with all the solutions in this set and are not added to the Pareto set.This property is already proven by Paquete et al. (2007) for best improvement strategies where all the neighbourhoods from all the solutions in A are explored entirely.The property also holds for I F because, like for I B , all the incomparable solutions from the neighbourhood of a solution in the Pareto local optimal set are visited.
Consider I N and a solution s ∈ A .Consider there are two solutions s and s in N (s) that are incomparable with s and all other solutions in A but that are not in the neighbourhood of each other.If s is visited before s , then s is added to A but s is discarded and s.visited set to true and thus s is not visited again.We conclude that PLS using I N stops in a subset of a maximal Pareto local optimum set.
Third, we show that PLS is always terminating.We show that a solution s eliminated from the current archive A because it is dominated by a newly added solution s , where f(s ) ≺ f(s), will not be added a second time into the archive A .Suppose that the same solution s is compared with solutions in A a second time.
A contains solutions that are dominating or incomparable with s .Thus, A contains solutions that are dominating s, and s is not included in this archive.If there are no other solutions that dominate the solutions from A , the PLS algorithm from Algorithm 1 stops by design.
Note that all three strategies add the dominating solutions, keep the incomparable solutions and delete the dominated solutions from the current Pareto set A .The difference between these three Pareto neighbourhood exploration strategies is their stopping criterion.I N stops the earliest, since it only needs to find a new incomparable solution.I B stops after all solutions are visited, and I F stops somewhere in-between I N and I B .

Genetic PLS (GPLS) for multi-objective quadratic assignment problems (mQAPs)
In order to design the neighbourhood, N , the stopping criterion, T , and genetic operators for GPLS, we need to specify a type of multi-objective combinatorial optimization problem.In this section, we design instances of GPLS for mQAPs.The neighbourhood and the stopping criterion are applicable for all kinds of permutation problems.For an optimal performance, the genetic operators are designed for each problem they are applied on.
Multi-objective QAPs Single and multi-objective QAPs are NP-hard combinatorial optimization problems that model many real-world problems (i.e., scheduling, vehicle routing, etc.).Intuitively, QAPs can be described as the (optimal) assignment of n facilities to n locations.A distance is specified between each pair of locations, and for each pair of facilities the amount of materials (or flows) transported between these facilities is given.The goal is to find the assignment of facilities to locations that minimizes the sum of the products between distances and flows.We consider multi-objective QAPs introduced by Knowles and Corne (2003).These mQAPs have for each objective different flow matrices and a single distance matrix.The flow matrices are correlated with some correlation ρ.Let us consider n facilities, a set Π(n) of all permutations of {1, 2, . . ., n} and the n × n distance matrix D = (d ij ), where d ij is the distance between location i and location j .We assume an m objective space, and m flow matrices B k = (b k ij ), each with n × n elements, where b k ij represents the flow from facility i to facility j in the k-th objective.The goal is to minimize for all objectives the set of cost functions , where π(•) is a permutation from Π(n).It takes quadratic time to evaluate this function.
The neighbourhood N (•) Multi-objective QAPs are permutation problems.A suitable neighbourhood operator for mQAPs is the q-exchange operator that swaps the position of q facilities.For example, the 2-exchange swapping operator swaps the position of two different facilities.This operator is attractive because of its linear time to compute the change in the cost function with the condition that all matrices D and B k are symmetrical (Paquete and Stützle 2006).We have The size of the neighbourhood increases quadratically with the number of facilities.
Stopping criterion, T The stopping criterion for GPLS is chosen to fairly compare its performance with MPLS.The distance between two solutions is defined as the minimum number of exchanges needed to obtain one solution from the other solution (Schiavinotto and Stützle 2007).The distance between a solution and the solution obtained with q-exchange mutation is q − 1.The search in GPLS is halted when the same number of swaps is executed as with MPLS.The difference in cost function, c k , for two solutions with distance q is linear in the number of facilities and the number of exchanges.Counting the number of swaps is equivalent to counting the number of fitness evaluations.

Genetic operators for GPLS on mQAPs
In this section, we present two perturbation operators for mQAPs (Drugan and Thierens 2010).In Algorithm 4, the Mutation algorithm is the parametrized mutation, and the path-guided mutation has the role of Recombination.We use cycles (Schiavinotto and Stützle 2007) to design operators that generate children at a given distance from their parent.These types of operators allow us to analyse the behaviour of GPLS on the tested problem instances.
The parametrized mutation, or the q-exchange mutation, uniform randomly selects, without replacement, q > 2 distinct locations in a solution s, {l 1 , . . ., l q }.To generate a new solution, these locations are exchanged from left to right or from right to left with equal probability to not bias the generation of individuals.When positions are exchanged from right to left, a position l i takes the value of its right neighbour l i+1 .Then, in this sequence, b ← l 1 , l 1 ← l 2 and so on until l q−1 ← l q and l q ← b, where b is a buffer variable.Note that s and the resulting solution form a cycle of size q that contains the mutated positions {l 1 , . . ., l q } and are q swaps apart.When PLS uses the 2-exchange operator to generate a neighbourhood, the mutation operator should exchange at least 3 facilities to escape from the region of attraction of a Pareto local optimal set.
The path-guided mutation, or the q-exchange path-mutation, uses two solutions, s and s , uniform randomly selected without replacement from the current Pareto set such that the distance between them is at least q.A child s is generated by copying the first parent s , and mutate s.The second parent s , is used to construct the path between the two solutions.The set of cycles common for the two solutions, s and s , are identified.A cycle is a minimal subset of locations such that the set of their facilities is the same in both parents, s and s .It is possible to switch a subset from one parent to the other one while keeping a valid permutation.
A cycle, c, is randomly chosen.For q − 1 times, choose at random a position i in the cycle c from solution s, where s[i] = s [j ] and i = j .Exchange the values of s[i] and s [j ].With this swap, the distance between s and its first parent, s , is increased with 1 and the distance between s and the second parent, s , is decreased with 1.If the size of c is smaller or equal than q, a second cycle is chosen.This process of randomly selecting a cycle and swap locations is repeated until the distance between s and s is becoming q.Then distance between s and s is − q + 1.If there are no parent solutions at distance larger or equal with q, we generate a solution with parametrized mutation.
The path-guided mutation is respectful: if the two parents have the same facility on a position i then their child will also have the same facility on the i-th position.Consequently, to be able to generate any solution in the search space we need to alternate path-guided mutation and mutation operators.The path-guided mutation operator resembles the path relinking operator (Jaszkiewicz and Zielniewicz 2009) because it generates solutions on the path between two solutions.
Cycle recombination (Poon and Carter 1995) uses cycles to generate children, but, unlike path-guided mutation, it exchanges an entire block of cycles between parents.Partially mapped crossover (PMX) (Goldberg and Lingle 1985) is the most popular recombination operator for the (multi-objective) travelling salesman problem (TSP).Some positions are uniform randomly chosen in the two parents and information is exchanged between them such that the generated children are valid permutations.Unlike path-guided mutation, PMX creates, or connects, cycles rather than respecting them.In Sect.6, we show that operators that do not respect cycles function worse than the operators proposed here that do respect them.

The expected number of solutions compared with a Pareto improvement strategy
In this section we derive a formula to compute the expected number of solutions generated with the three Pareto neighbourhood exploration strategies.We show that the number of solutions generated with I N and I F are in expectation much smaller than with I B .Let us consider a neighbourhood N (s) generated with the 2-exchange operator that has: z-the number of solutions that dominate s, z The best Pareto improvement strategy evaluates all the solutions in the neighbourhood.The number of solutions evaluated with The probability of evaluating a certain number of solutions with I F follows a hyper-geometrical distribution describing the number of successes in a finite sequence where samples are drawn without replacement.
Lemma 1 Let I F and N (s) be as before.Let z be the number of solutions in N (s) that dominate s.The probability to find, after exploring exactly i solutions from N (s), the first solution that dominates s is Proof I F stops at the first improvement in the sequence of solutions generated from N (s).In total there are z improvements in this sequence of C n 2 solutions.Let successes be associated with improvements.The number of improvements, and thus success, before I F stops is 1.Using the definition of the hyper-geometrical distribution, we assume that there are z successes and C n 2 − z failures.The probability of finding the first improvement after i trials from the z total improvements is the hypergeometrical value divided by i, z 1 For the neutral Pareto improvement strategy, the probability to discover in i steps the first solution that dominates or is incomparable with s is Proposition 7 Let I B , I F , I N be as before.Let N be the neighbourhood generated with 2-exchange operator, and let z be the number of solutions that dominate s, v the number of solutions that are incomparable with s, and π N (s) a uniform randomly generated permutation of N (s) as before.
The expected number of evaluated neighbours with I B is where E{|g(•)|} is the expected value of the function g.The expected number of evaluated neighbours with I F is The expected number of evaluated neighbours with I N is Proof The number of evaluated neighbours is always the same for I B , C n 2 .Thus, To calculate the expected evaluated solutions for I F and I N , we use the definition of the expected value for the function g(•), E{g} = x g(x) • P (x).
Here, the probability P (x = (i)) is, as in Eq. 2, the probability that a Pareto improvement strategy generates i solutions before the first improvement on s is met.The sum is over all i = {1, . . ., C 2 n − z}, where C 2 n − z is the maximum number of solutions that can be evaluated before the first improvement on s is found.The function g(i) = i as usually.Then, the expected number of evaluated solutions for first Pareto improvement strategy is To proof the expected value for I N , we follow the same line of reasoning.Some properties of these expected values are given in the following proposition.
, from Eqs. 2 and 1, decrease.Thus, E{|I F |} and E{|I N |} decrease.Proving that when v increases, E{|I N |} decreases follows the same line of reasoning."Flat" multi-objective landscapes An interesting, but logical, property is that I N and I F return less solutions in a flat landscape (i.e.landscape with many incomparable solutions) than I B .Single objective spaces with many equal values are called flat and they are known as inappropriate for heuristics in general.Multi-objective spaces with many incomparable values are more common than flat single objective spaces.In general, all the subsets of elements {x i , . . ., x j } ∈ S such that f ) form a Pareto set.For "flat" spaces, or spaces with a large "flat" region, when PLS uses I B , the entire region is explored and added to the current Pareto set.That leads to computational and storage problems for PLS.In this case, a Pareto improvement strategy, like I F and I N , that finds a single or few solutions from the space is less demanding in memory and run-time than I B .

Experimental results
In this section, we experimentally compare six different (stochastic) Pareto local search algorithms.The first three are multi-start PLS algorithms (MPLS) each using a different Pareto neighbourhood exploration strategy: best improvement I B , first improvement I F , and neutral improvement I N .The other three are SPLS algorithms that combine the above three neighbourhood exploration strategies with stochastic perturbation that consists of unbiased mutation and path-guided mutation, each applied with equal probability.All algorithms generate their neighbourhood N with the 2-exchange operator.We denote the six algorithms as follows: 1. bMPLS uses MPLS(I B , T ), where T means that PLS is restarted for M = 100 times; 2. fMPLS uses MPLS(I F , T ), where T means that PLS is restarted until the same number of swaps are generated as with bMPLS; 3. nMPLS uses MPLS(I N , T ), with T as in fMPLS; 4. bGPLS uses GPLS(I B , 0.5, T 1 , T 2 ), where T 1 randomly restarts PLS N = 10 times and T 2 means that PLS is restarted using genetic operators until the same number of swaps are generated as with bMPLS; 5. fGPLS uses GPLS(I F , 0.5, T 1 , T 2 ), with T 1 and T 2 as in bGPLS; 6. nGPLS uses GPLS(I N , 0.5, T 1 , T 2 ), with T 1 and T 2 as in bGPLS.
We generate restarting solutions with both mutation and path-guided mutation with an exchange rate q, uniform randomly selected between 3 and n/3 .These exchange rates are chosen based on the empirical observations from Drugan andThierens (2010, 2011).Also the landscape analysis in the paper shows that these exchange rates are the most effective.For example, for n = 25, we have q = {3, . . ., 8}, and for n = 100, we have q = {3, . . ., 33}.Each algorithm independently runs 50 times.
It should be noted that we have performed preliminary experiments with unstructured and structured with negative and zero correlation bQAPs.As Paquete et al. noticed, on negative unstructured bQAP instances more complicated LS-based heuristics and the multi-restart PLS have approximatively the same performance.The Pareto sets of these bQAPs are very large, for instance for ρ = 0 and n = 75 the Pareto optimal set typically found has about 2000 solutions.This can be explained by the fact that when there is no correlation, many solutions have a different ordering for the two objectives.Similarly when there is negative correlation (ρ = −0.75), the sizes of the Pareto optimal sets typically found are huge.In this case, the ordering of each solution is opposite in the two objectives.For such large Pareto optimal sets the search with PLS using I B therefore becomes too expensive.Knowles and Corne (2004) acknowledge that the structured QAPs are easier than the unstructured bQAPs for local search.
To summarize, we consider that the complexity is highest when there are a large number of attractors containing Pareto optimal solutions, compared to few attractors with a large number of Pareto optimal solutions.Section 5.1 looks at the performance of the six algorithms, Sect.5.2 discusses their exploration properties, while Sect.5.3 focuses on their dynamical behaviour.

Comparing the performance of stochastic PLS algorithms
Pareto sets are typically compared by the hypervolume unary indicator and/or attainment functions (Zitzler et al. 2003).The unary attainment function (Fonseca and Fleming 1996) (EAF) gives the probability of attaining each point (independently) in the objective space.Contour surfaces through certain probabilities can then be drawn.To compute the hypervolume unary indicator, we use the PISA package (Bleuler et al. 2003).For comparison purposes we first normalize the outputs of all six algorithms, and the corresponding reference set, if it exists, from http://eden.dei.uc.pt/~paquete/qap/,A * (n, ρ), where n = {0.25,0.5, 0.75} and ρ as before.For each instance, the reference set A * is obtained with a two phase tabu algorithm that is run much longer than the MPLSs are run.There are algorithms that outperform this tabu algorithm in some cases (Lopez-Ibanez et al. 2004).This normalization script assigns to the smallest point(s) in an objective the value 1 and to the highest point(s) the value 2. All the other points are scaled to a value between 1 and 2 in both objectives.The normalized outputs are the inputs for the hypervolume calculation.The reference point is (2.1, 2.1) which has the largest value in both objectives.Table 2(A) shows the average and the standard deviation of the hypervolume unary indicators of the six SPLSs on twelve testing bQAPs.Table 2(B) compares the six SPLS algorithms: (i) the mean hypervolume for the three values of ρ and four numbers of facilities, n, (ii) the number of wins of fGPLS, bGPLS and fMPLS over the six SPLSs, and (iii) the rank of each SPLS algorithms based on the number of wins.For the same bQAP instance, we measure if the difference between two measurements of Table 2(A) is statistical significant using the Wilcoxon non-parametric two samples test, where p < 0.05. Figure 4 presents five EAFs at 2 %, 25 %, 50 %, 75 % and 100% for two algorithms, fMPLS and fGPLS, on two bQAPs, bQAP(75, 0.5) and bQAP(75, 0.75).Overall, fGPLS is the best tested SPLS and nMPLS is the worst algorithm.
All GPLS algorithms -that are fGPLS, nGPLS, bGPLS-outperform all the multirestart PLS-that are bMPLS, fMPLS and nMPLS-on the twelve tested bQAPs.In Fig. 4, the EAFs contours are closer to the lowest point (1, 1), and thus better, for fGPLS than for fMPLS.Furthermore, for fGPLS the best EAF contour at 2 % is Fig. 4 Attainment surfaces at 2 %, 25 %, 50 %, 75 % and 100 % (lines from bottom left to top right) of the normalized outputs.The two tested SPLSs are: (i) fMPLS on the left and (ii) fGPLS on the right.The two tested bQAPs are: (i) bQAP(75, 0.5) top and (ii) bQAP(75, 0.75) bottom.These EAFs are compared with the best known Pareto set for the two bQAPs: (i) A * (75, 0.5) top and (ii) A * (75, 0.75) bottom.The objective values found at 2 % EAFs correspond to the best Pareto set found over 50 runs and those found at 50 % EAFs are the median outcome better than 1 % attainment surface for A * from 100 runs on both bQAP instances.The objective values of 2 % EAFs are better for A * than for fMPLS.Thus, the two perturbation operators of GPLS increase the efficiency of the SPLS algorithms.
Table 2(B) shows that the SPLS algorithms that use perturbation operators have the highest performance, rank 1, 2, and 3.The difference between GPLS and its corresponding MPLS with the same Pareto neighbourhood exploration strategy increases when the correlation ρ and the number of facilities n increase.For each SPLS, the hypervolume indicator decreases with the increase of ρ and n.The smallest, statistically the same, hypervolume values are registered with the three MPLSs on bQAP(100, 0.75).On the same bQAP instance, the smallest difference between the hypervolume of GPLS and MPLS is 0.41 when I B and I N are used.The largest hy-pervolume difference, 0.45, is between fGPLS and fMPLS.Note that, on bQAP(100, 0.75) the hypervolume indicators of GPLSs are approximatively two times larger than the hypervolume indicators of MPLSs.On the contrary, on bQAP(25, 0.25) all SPLS algorithms perform about the same, the exception is nMPLS, which has a lower performance.We conclude that the difference between the performance of GPLS and MPLS increases with increasing problem complexity-that is the increase in n and ρ.
In Fig. 4, on bQAP(75, 0.5), attainment surfaces, EAFs, are clustered together, whereas on bQAP(75, 0.75) they are spread.Attainment surfaces of fGPLS are more spread than the attainment functions of fMPLS.In Table 2, the clustered EAFs correspond to a small variance of hypervolume values whereas the spread EAFs correspond to large variances of the hypervolume indicators.On the other side, on bQAP(75, 0.5), in Fig. 4(top), each attainment surface is reasonably well spread over the two objectives.On bQAP(75, 0.75), in Fig. 4(bottom), the attainment surfaces are rather narrow and have a small amount of points.For large correlations, there are few incomparable solutions and a large difference between the worst and the best objective values.For small correlations, there is a large amount of solutions with different precedence in the two objectives and a smaller difference between the worst and the best objectives.
The SPLS algorithms that use I F have a better or equal performance than the corresponding SPLSs that use I B .In Table 2(B), both the best genetic PLS, rank 1, and the best multi-restart PLS, rank 4, use I F .The SPLSs that use I F or I B have a better or equal performance than the corresponding SPLSs that use I N .For bQAPs with low correlation, ρ = 0.25, there are no statistical differences between algorithms using I B and I F .bQAPs with high correlation, ρ = 0.75, give the largest difference between SPLSs that use I F and I B .For bQAPs with ρ = {0.25,0.5}, SPLSs that use I N are the algorithms with the worst performance.The bad performance of SPLSs using I N is due to the long walk on the Pareto set of incomparable solutions instead of searching for dominating solutions that greatly increase the hypervolume.For ρ = 0.75, SPLSs that use I N have similar performance with SPLSs that use I B .The small number of solutions in the Pareto sets of bQAPs with ρ = 0.75, restricts the time I N spends in a neighbourhood of a solution.
The overall conclusion is that the SPLS algorithm with the highest performance is fGPLS using the first improvement Pareto neighbourhood exploration I F and both perturbation methods, mutation and path-guided mutation.The bQAP, where the difference between fGPLS and the other SPLSs is maximal, has the largest number of facilities and the largest correlation, bQAP(100, 0.75).In the next section, we give more insights on the SPLSs' behaviours on the tested bQAPs.

Landscape characteristics of bQAPs explored with SPLSs
Table 3 shows properties of bQAPs when explored with SPLSs.Table 3(A) shows the number of times PLS is restarted in a SPLS run, #PLS, where all the SPLS algorithms have the same run-time, measured in number of swaps.Table 3( B) shows that the number of calls of the Pareto improvement strategy, #I is proportional with n and inverse proportional with ρ.In Table 3(C), the number of solutions in the output Pareto local optimal set of a PLS, |A|, is proportional with n and inverse proportional with ρ.  6) 80 ( 4) 571 (1) 116 ( 5) 563 ( 2) 126 ( 4) bQAP (C) The number of solutions in the Pareto local optimal set |A| 25 0.25 28 There is a natural relationship between Table 3(A), (B) and (C).In Sect.4.2, we showed that the expected number of solutions generated with I B is much larger than the expected number of solutions generated with I F , which, at its turn, is much larger than the solutions generated with I N .The more often the Pareto improvement strategy I is called, the less times SPLS has to restart PLS.Thus, #I and the number of solutions generated with one of the Pareto neighbourhood exploration strategies are inversely correlated with #PLS, while #I is correlated with |A|.Logically, the larger the size of a Pareto local optimal set, the larger the number of I calls that are necessary to explore the neighbourhood.For all algorithms, the largest #I is obtained for bQAP(100, 0.25) which has the largest fronts |A|.
The least number of times a SPLS restarts PLS is 100 for bMPLS, closely followed by fMPLS.Similar with bMPLS, Paquete and Stützle (2006) restarts PLS 100 times.bGPLS and fGPLS restart up to 3-3.5 times more often.The largest number of times PLS is restarted, #PLS, is obtained with nGPLS, which can be up to 40 times larger than 100.When I B is used, #PLS is independent of n and ρ because the expected number of solutions visited with I B is also constant with n. #PLS of SPLSs using I N varies greatly with n and ρ because the expected number of solutions evaluated with I N varies greatly with n and the number of incomparable and dominating solutions in the neighbourhood.As expected, #PLS when I F is used has values between #PLS using I B and I N .
Interesting, for I N , the smaller ρ is, the larger #PLS is, whereas for I F , the larger ρ is, the larger #PLS is.We corroborate this observation with the observation that bQAP with lower ρ have more incomparable solutions, thus larger Pareto local optimal sets.The lower the correlation ρ, the lower the number of solutions in the neighbourhood that are visited with the I N exploration strategy, and the more often the SPLS using I N has to restart.On the other side, the SPLS using the I F exploration strategy runs longer for low ρ because of the large number of incomparable solutions in the neighbourhood and the smaller number of dominating solutions that would stop I F .
#I and |A| depend on the type of SPLS, GPLSs have smaller Pareto local optimal sets and smaller #I than MPLSs.The SPLSs using I N have a much smaller |A| than the SPLSs using I B and I F .This is caused by the numerous incomparable solutions in Pareto local optimal sets of bQAPs with low correlation, ρ = 0.25.Recall that I N stops after adding a single incomparable solution to A, whereas I F and I B add a set of solutions to A.
To conclude, GPLSs spend less time in reaching a Pareto local optimal set than MPLSs do.SPLSs using I N restart most often and they have the worst performance because I N prematurely stops the evaluation of incomparable solutions on large Pareto sets.SPLSs using I F have the best performance and are restarted more often than SPLSs using I B .

Dynamical behaviour of SPLSs on bQAPs
In this section, we empirically analyse the dynamical behaviour of SPLSs using the relationship between Pareto local optimal sets generated at consecutive iterations of PLS.Let us consider any two consecutive iterations of PLS.The Pareto local optimal Fig. 5 The correlation of fMPLS on the left and of fGPLS on the right on bQAP(75, 0.75): (top) between their distance in objective space of consecutively restarted PLSs and their distance in solution space, (bottom) between the hypervolume of restarted PLSs and their average distance in solution space set of the first PLS is denoted as the parent Pareto local optimal set, or shortly the parent Pareto set.The output of the second PLS is denoted as the offspring Pareto local optimal set, or shortly the offspring Pareto set.The distance in the objective space, obj , between any parent and its offspring Pareto set is calculated as the average Euclidean distance between any two solutions, one from the parent and one from its child.The distance in the solution space, sol, between a parent and its offspring Pareto set is calculated as the average of the minimum number of swaps between any two solutions of the two Pareto local optimal sets.For this experiment, Pareto local optimal sets of a single random run of SPLSs are considered.
Figure 5(a) and (b) show, on bQAP(75, 0.75), the relationship between obj and sol. Figure 5(c) and (d) present, on bQAP(75, 0.75), the relationship between the distance in solution space, sol and the hypervolume of the offspring Pareto set.Note the resemblance between the figures for the same SPLS.For fMPLS, solutions   4 generalizes the results from Fig. 5.For MPLSs, the distance in solution space, sol, is proportional with n and is independent of ρ and I.For GPLSs, sol is proportional with n and inverse proportional with ρ. sol has the smallest values for GPLSs using I B , and the largest values for GPLSs using I N .
sol for GPLSs is smaller than for MPLSs.obj is proportional with n and inverse proportional with ρ.For GPLS, a large distance is objective space obj implies also a large distance in solution space sol.obj has the smallest values for SPLSs using I B , and the largest values for SPLSs using I N .
To conclude, genetic PLSs perform better than multi-restart PLSs by exploiting the structure of the bQAP search space, a structure caused by the positive correlation between the objectives.As shown above, the distance in solution space is correlated with the distance in objective space between consecutive Pareto local optimal sets and with the hypervolume indicator.

Random walks for landscape analysis of bQAPs using SPLSs
The previous section showed that the best performance was obtained with the SPLS using genetic perturbation operators and the first Pareto improvement neighbourhood exploration strategy.We presented some experimental results where GPLSs exploit the correlation in bQAPs.In this section, we deepen our analysis using random walks to investigate the effects of various mechanisms of stochastic PLSs in isolation.A random walk (RW) is a sequence of steps generated with a random process on the path between two "good" solutions in the search space using a genetic operator.
Section 6.1 shows the usefulness of the genetic operators that use cycles.The efficiency of deactivation for different correlation values of bQAPs is studied in Sect.6.2.A novel landscape analysis with RW for multi-objective spaces is given in Sect.6.3.The properties of different Pareto neighbourhood exploration strategies are investigated with RW in Sect.6.4.
Methodology Let s and s be two solutions selected uniform randomly without replacement from the best known Pareto local optimal set A * of a bQAP instance.Let (s, s ) be the distance between the two solutions and let s 1 be a solution generated on the path between two solutions s and s .In this experiment, (s, s ) = 70.For each value of q ∈ {1, 2, . . ., (s, s )}, 100 solutions are generated by swapping q pairs of facilities with each genetic operator.
Using random walks, we study the relationship between the quality of the solutions generated on the path at a certain distance from s and the quality of the Pareto local optimal set of PLSs started from these solutions.We consider two initial Pareto sets A for PLS(A, I): (a) A = {s 1 }, which does not consider deactivation, and (b) A = Deactiv(s 1 , A * ) which considers deactivation.By definition, Deactiv(s 1 , A * ) contains s 1 and all the solutions in A * that are incomparable with s 1 .For the ease of the discussion, we have set the visited flag of the solutions from A * included in Deactiv(s 1 , A * ) to true.When no deactivation is used, the exploratory archive is empty and solutions are only compared to the newly generated solution s.Later, the newly found solutions are merged with the complete archive, so that only non-dominated solutions remain.Another alternative would be to restart PLS from a set merge({s 1 }, A * ).Preliminary experiments showed a very poor performance of such an algorithm because the probability that merge({s 1 }, A * ) = A * is very small.Most of the time, such PLS is stuck in the same, early found, Pareto local optimum set.
The quality of a Pareto local optimal set is measured with the hypervolume unary indicator computed as indicated in Sect. 5.The objectives of all generated Pareto local optimal sets, including the objectives in A * , are normalized, each objective having values between 1 and 2. The reference solution for the hypervolume unary indicator is, as before, (2.1, 2.1).A large hypervolume value means a better Pareto local optimal set.A PLS is restarted for each of the 100 points generated at certain distances on the path.

Cycle vs non-cycle based genetic operators for bQAPs
It is important to recall that the perturbation operators we used in Sect.4.1, mutation and path-guided mutation, preserve the cycles in the parents.To motivate our choice, we show here that these operators are more efficient than perturbation operators that do not respect cycles.The efficiency is measured with the hypervolume unary indicator and the Euclidean distance in objective space.In Table 2, the largest difference between GPLSs and MPLSs are for bQAPs with ρ = 0.75.We focus here on this test problem .
We compare four operators.The two operators that respect cycles are: (i) the q-exchange mutation, mut, and (ii) the q-exchange path-mutation, path.We also consider two operators, mutation mut R and path-guided mutation path R , that do not respect cycles.mut R swaps q uniform randomly chosen pairs of facilities, with replacement.path R swaps q pairs of facilities such that, for each swap: (i) the first location, i, is uniform randomly chosen from a parent s, and (ii) the second location, j , is chosen to have the same value in s like the value in the i-th location in s.
In Fig. 6, the x-axis represents the distance in solution space between the initial solution s and the solution s 1 , (s, s 1 ) generated with one of the four operators.Note that for the operators that do respect cycles, mut and path, we have (s, s 1 ) = q.For operators that do not respect cycles, mut R and path R , we have (s, s 1 ) ≤ q.To see this, one should note that by randomly sampling positions for pairing, some positions might be swapped two or more times.In this case the distance between the generated and the initial solution does not increase.In Fig. 6(a) and (b), the y-axis denotes the Euclidean distance between the objectives of s and s 1 , (f (s), f (s 1 )).In Fig. 6(c) and (d), the y-axis denotes the hypervolume unary indicator for a Pareto local optimal set generated by PLS(Deactiv(s 1 , A * ), I B ) from Algorithm 1.The hypervolume indicator value for A * is shown at the borders in Fig. 6(c) and (d).
In Fig. 6(a) and (c), mut generates solutions that are closer in objective space to the initial solution, s, than the solutions generated with mut R .The difference between the two operators, in both objective space and hypervolume indicators, is diminished for large (s, s 1 ), i.e.50 to 70 swaps.Note that for a large amount of swaps, s 1 is actually randomly generated.The largest difference between the two operators is for 20 to 30 swaps.This indicates a better exploration of the structure of the problem for mut than with mut R .Wilcoxon non-parametric two-side tests showed that these differences in hypervolume are significant.We conclude that mut R is worse than mut.
In Fig. 6(b) and (d), path generates solutions that are closer in objective space to s and s 1 than solutions generated with path R .Similarly with mutation, the largest Fig. 6 Random walk (RW) between two solutions s and s from A * (75, 0.75) for four genetic operators (a) mut and mut R and (b) path and path R .The top part shows the Euclidean distance in objective space between the initial solution s and its offspring s 1 on the path towards s , (f(s), f(s )).The bottom part shows the hypervolume of a Pareto local optimal set of PLSs that are started from the solutions generated from s 1 on the path between two solutions, s and s , PLS(Deactiv(s 1 , A * ), I B ) using (c) mut and mut R , and (d) path and path R difference between path and path R is between 20 and 30 swaps which coincides with the largest and most significant difference between the corresponding hypervolume indicators.Note that, for path R , there are no solutions s 1 generated at large distances from s, (s, s 1 ) > 40.
Let us compare mutation mut and path-guided mutation path.The solutions that are far-away from the first parent in path are, by definition, close to the second parent and thus bounded.For (s, s 1 ) > 35, the solutions generated with path have lower Euclidean distances in objective space than the solutions generated with mut.The corresponding hypervolume unary indicator for path is significantly larger than for mut.Thus, the path-guided mutation is better than mutation in exploring the structure of the landscape.Using the same arguments, path R is better than mut R .
To conclude, the operator that best explores the structure of the search space is path which preserves both the cycles and walks on the path between two good solutions.It is interesting to note the relationship between the hypervolume indicator of a Pareto Fig. 7 The Pareto local optimal set of PLSs started from s 1 on the path between two solutions in A * (75, 0.5): (a) using the deactivation technique PLS(Deactiv(s 1 , A * ), I B ), (b) without deactivation PLS(s 1 , I B ). (c) The number of solutions in A * (75, 0.5) that are incomparable with the solutions generated with mut and path on the path between s and s local optimal set of PLSs restarted from a solution s 1 of the path and the distance from the initial solution s to s 1 .This is a particularity of bQAP's landscape captured also in Fig. 5.It basically means that solutions that are correlated in solution space are also correlated in objective space.That is why genetic PLSs are so much better than the multi-restart PLS for highly correlated bQAPs.

The use of deactivation in bQAPs
In this section, we measure the effect of the deactivation technique on PLS's performance.In Fig. 7(a) and (b), we compare the hypervolume indicators of PLSs restarted from s 1 on the path between two solutions s and s from A * (75, 0.5).In Fig. 7(c), we count the solutions from A * (75, 0.5) that are incomparable with the generated solution s 1 .
The hypervolume indicator for PLSs restarted with solutions generated with path and using the deactivation technique are the largest for (s, s 1 ) ∈ {1, . . ., 15} ∪ {55, . . ., 69}.For the same distances, there is a positive number of solutions, maximum 18 from a total of 24, that are incomparable with s 1 .For path, the increase in the quality of Pareto local optimal sets generated with PLS restarted in s 1 is reflected in the increased number of incomparable solutions in A * .For mut, when s 1 is such that (s, s 1 ) ∈ {1, . . ., 15}, the number of incomparable solutions is as large as for solutions generated with path.For a large number of swaps, the number of incomparable solutions in A * is 0 indicating a rapid decrease in quality of the generated solutions.Note that for solutions generated with mut the deactivation technique does not increase the quality of the output of the restarted PLSs.
Note the difference in the hypervolume indicator of bMPLS from Table 2(A) and the hypervolume indicators for PLS restarted at a large distance from s from a solution s 1 generated with mut from Fig. 6(c) and (d), and Fig. 7(a) and (b).In this section, to compute the hypervolume indicator, we normalize all the Pareto local optimal sets of PLSs restarted on the path and A * of a bQAP.In the previous section, we normalized all the Pareto local optimal sets of the tested SPLSs.Some of the solutions in these Pareto sets are dominating solutions in A * , and some of these solutions are much worse than the Pareto local optimal sets from the experiments in this section.
For bQAP(75, 0.75), for both mut and path, there are few solutions, at most 2 from a total number of solutions of 11, that are incomparable with the current solution s 1 in A * (75, 0.75), where (s, s 1 ) ∈ {1, 2}.For the rest of solutions s 1 on the path, (s, s 1 ) > 2, there are no incomparable solutions in A * .
In conclusion, for ρ = 0.75, the size of Pareto local optimal sets are small and the deactivation mechanism does not work.For lower correlations, ρ = 0.5, the size of Pareto local optimal sets are larger and the deactivation selects several incomparable solutions, which improves results a lot with the path-guided mutation operator.

The properties of PLSs restarted from RW solutions
In this subsection, we investigate the properties of PLSs that are started from the solutions on the path between s and s , s, s ∈ A * , using one of the genetic operators mut and path.Consider PLS(Deactiv(s 1 , A * ), I B ) and bQAP(75, 0.5) as before.
Let us measure how many times in a PLS run the Pareto improvement strategy I B is called.The number of swaps in the PLS run is the number of calls of I B multiplied with the number of solutions in a neighbourhood, i.e.C n 2 with n the number of facilities.Figure 8 (a) shows the relationship between (s, s 1 ) and #I B , the number of calls of I B in the PLS started from s 1 .For mut, (s, s 1 ) is proportional with #I B .For path, the largest #I B is obtained for (s, s 1 ) = (s, s )/2.In Table 3(B), for bQAP(0.75,0.5), #I B for bMPLS is 457 that is approximately the same as #I B for PLSs restarted using mut at a large distance from s.The number of restarts I B for bGPLS is 67 which corresponds to small (s, s 1 ) ∈ {10, . . ., 20}.Recall that the exchange rates for mutation and path-guided mutation for GPLSs on bQAP(75, 0.5) are uniform randomly selected from {3, . . ., 38}.The number of restarts in bGPLS cannot be directly related with #I B from this experiment.
Let us compare I B from Fig. 8  In Fig. 8(b), the size of a Pareto local optimal set, |A| between 10 and 50 solutions, is increasing with the increase of (s, s 1 ).Again, for a large (s, s 1 ), |A| is larger for solutions generated with mut than for solution generated with path.This relationship is found back in Table 3(C) for bQAP(75, 0.5) where |A| for bMPLS is 35 and for bGPLS is 22.This also means that larger archives do not always correspond to better Pareto local optimal sets.
To conclude, the results obtained with RW are closely related to the results in Sect.5, meaning that RW can be used in the context of SPLS to study the performance of these algorithms.

RW landscape analysis for different Pareto neighbourhood exploration strategies
In this section we compare the properties of the three Pareto neighbourhood exploration strategies, I B , I F and I N , on a bQAP instance bQAP(75, 0.75).We show that the basins of attraction are explored differently by the three Pareto neighbourhood exploration strategies.
In Fig. 9(a), if s 1 is generated with mut such that (s, s 1 ) ≥ 30, PLS using I B escapes the basin of attraction of s.When path is used, the probability to find s , the second parent of s 1 , in PLS's output is non-zero if the distance between s 1 and s , (s 1 , s ), is lower than 25.In Fig. 9(b) and (c), PLSs that use I F and I N explore less well the basins of attraction of s and s than I B does.The probability to return in s with I F is non-zero for (s, s 1 ) < 20 and it is smaller than with I B .This is because the search is stopped before all the solutions in the neighbourhood are explored.As expected, PLSs using I F explore better the basins of attraction than PLSs using I N .
To conclude, note that the "optimal" number of exchanges in the two genetic operators is always a trade-off between the probability of leaving the basin of attraction, and the probability of generating a new local optimum that is closest to-and thus Fig. 9 The probability that a Pareto local optimal set of PLSs started on the path between two solutions in A * (75, 0.75) contains the parent solutions.PLS uses one of the following Pareto neighbourhood exploration strategies: (a) I B , (b) I F , and (c) I N most correlated with-the current Pareto local optimal set.For a small number of exchanges, the solutions generated are "fit" in the sense that the objective values approach the objective value of the initial solution and the hypervolume of the Pareto local optimal sets of PLS restarted from these solutions is high.The smaller the number of exchanges, the higher the probability that solutions belong to the basin of attraction of the initial solution.

Conclusions and discussion
We introduced stochastic PLS (SPLS), a class of algorithms that combine Pareto local search with stochastic perturbation operators to efficiently explore multi-objective search problems.The paper focuses on three issues: (i) Pareto neighbourhood exploration strategies to explore the Pareto neighbourhood of a solution, (ii) stochastic perturbation operators, mutation and path-guided mutation, to restart the Pareto local search, and (iii) the deactivation method to run PLS from perturbed solutions in the context of previous runs of Pareto local search.
We proposed a new concept for Pareto neighbourhood exploration strategies that compares the solutions in a neighbourhood with: (a) the current solution from a Pareto set, and (b) all the other solutions in that Pareto set which were previously generated with Pareto local search.We show that our Pareto improvement strategies, especially neutral and first Pareto improvement, improve the performance of similar strategies from the literature (Liefooghe et al. 2011).Interestingly, we showed that the use of the neutral Pareto improvement strategy with PLS is equivalent to the use of the first improvement strategy with hypervolume based local search.The hypervolume unary indicator is widely used to transform a multi-objective problem into a single-objective problem (Beume et al. 2009).However, our experiments indicate that SPLSs that use the neutral Pareto improvement-or equivalently, hypervolume based local search using the first improvement exploration strategy-perform the worst especially for bQAPs with a large number of incomparable solutions because of the lack of exploration around the evaluated solutions.SPLSs that use the first Pareto improvement strategy have the highest performance.
We presented genetic operators, mutation and path-guided mutation, specifically designed for multi-objective QAPs.We used random walks to show that our genetic operators that use cycles efficiently exploit this search space.We also added explanations on why some algorithms work better than the others.Overall, our experimental results showed that genetic PLSs outperform multi-restart PLSs for bi-objective QAPs that have a structured search space that is exploitable with genetic operatorsin our case positive correlation between objectives.
In the experimental section, we explained the difference in behaviour of six different (S)PLS algorithms using some characteristics of the landscape.Random walks were successfully used to analyse the fitness landscape and PLS's behaviour.We showed that the deactivation method improves the performance of PLSs for bQAPs with large Pareto local optimal sets, thus with a large amount of incomparable solutions.To make the results in this paper comparable with the results from Paquete's study, we restarted bMPLS the same number of times as Paquete did.Some of our SPLS algorithms performed better than Paquete's algorithm, better even than the tabu search algorithm which they run much longer in order to generate a reference set.The best performance with SPLS is obtained using the first Pareto improvement strategy, which stops the neighbourhood exploration when a dominating solution is found, and using stochastic perturbation operators which exploit the particularities of the search space.All the novel techniques proposed in this paper contributed to the efficiency of the SPLS algorithm.
To conclude, we can state that stochastic Pareto local search (SPLS) is a powerful framework for multi-objective, combinatorial optimization problems.Although we only showed experimental results on a single problem type (mQAPs), one can safely predict that SPLS will also be successful for other types of problems, just like stochastic local search (SLS) algorithms are successfully applied in a very large range of single-objective combinatorial optimization problems.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.

Algorithm 1
a) f(b) f(a) f(b) and f(b) f(a) A B A B and B A Non-dominated by f(a) f(b) f(a) f(b) or f(a) f(b) A B A B or A B Pareto local search PLS(A, I) Require: ∃s ∈ A, s.visited = false A ← A while ∃s ∈ A , with s.visited = false do A ← I(s, A \ {s}, N , , f) for all s ∈ A do s .visited← false end for A ← merge(A , A ) Set s.visited ← true end while return A We use a Pareto set A to restart PLS in Algorithm 1 as opposed to a single initial solution as Paquete et al.'s PLS does.Definition 1 Let's consider the deactivation method from Algorithm 2, Deac- tiv(s, A), with A a Pareto set and s a solution.The output is an archive A I containing s and the solutions of A that are incomparable with s.The differences between PLS from Algorithm 1 and Paquete et al.'s PLS are: (i) we use a Pareto set A to restart PLS, and (ii) we propose several instances for the Pareto Algorithm 2 Deactivation Deactiv(s, A) Require: A: a Pareto set A I ← {s} for all s ∈ A do if f(s) f(s ) then A I ← merge({s }, A I ) end if end for return A I improvement algorithm I. PLS restarting from a single solution s and using the best improvement neighbourhood exploration strategy resembles the most Paquete et al.'s PLS.The following proposition shows that deactivation can decrease the number of visited solutions, the extra solutions visited with PLS without deactivation that are dominated by solutions in the deactivation set.

Fig. 2
Fig. 2 Consider N (s) = {s 1 , s 2 , s 3 , s 4 , s 5 }.We compare the three Pareto neighbourhood exploration strategies (a) with an empty initial Pareto set or (b) with the initial Pareto set A = {p}

Proposition 6
Let I B , I N and I F be Pareto neighbourhood exploration strategies, and π N (s) the same permutation of the evaluated neighbourhood as before.Let us consider |I(s, A)| the number of solutions evaluated with I.Then, I N (s, A) ≤ I F (s, A) ≤ I B (s, A) and I B (s, A) I F (s, A) I N (s, A) Proof I B generates all the neighbours of s thus it generates always at least as many individuals as the other two exploration strategies, I F and I N .This also means that I B always finds the solution, if it exists, that dominates {s} ∩ A and all the other solutions N (s).Thus I B (s, A) I F (s, A) and |I F (s, A)| ≤ |I B (s, A)|.We have I F (s, A) I N (s, A) because I F stops in a dominating solution for s, if it exists, whereas I N stops in any non-dominated solution for s and A. For a permutation π N (s) , |I N (s, A)| ≤ |I F (s, A)| because I N is always stopping earlier than I F unless one of the conditions 3 or 4 hold in which case |I N (s, A)| = |I F (s, A)|.

Fig. 3
Fig. 3 The expected number of evaluated solutions in a neighborhood when (a) I B and (b) I F are used.The number of facilities n = {20, . . ., 40} and the number of dominating solutions z = {1, 5, 10, 15}

Proposition 8
If n increases, then each of the expected values E{|I N |}, E{|I F |}, E{|I B |} also increases.If z increases, then E{|I F |} and E{|I N |} decrease.If v increases, then E{|I N |} decreases.Proof The first property follows directly from the definition of expected value.If z increases, then the terms C

Figure 3
Figure 3 shows the expected number of evaluated solutions in a neighbourhood with (a) best Pareto improvement and (b) first Pareto improvement for several values of the number of facilities, n, and the number of dominating solutions, z.The expected number of evaluated solutions in a neighbourhood is always largest for best Pareto improvement, E{|I B |}.Even if only one solution in the neighbourhood is dominating the current solution, z = 1, E{|I F |} is half of E{|I B |}.As stated in (a) with the hypervolume indicators from Fig. 7(a).A small number of calls of I B corresponds with high hypervolume unary indicators.A large #I B corresponds with low hypervolume unary indicator.This means that

Fig. 8
Fig. 8 Consider the Pareto local optimal set of PLSs started from s 1 on the path between two solutions in A * (75, 0.5), PLS(Deactiv(s 1 , A * ), I B ).(a) The number of times I B is called, #I B .(b) The number of solutions in the Pareto local optimal sets, |A|

Table 1
Relations between objective vectors and sets of objective vectors considered in this paper

Table 2
The performance of the six SPLS algorithms on 12 bQAP instances

Table 4
The distance between parent and offspring in (A) solution space, sol, and (B) objective space, obj Pareto sets are at a large distance in solution space and they are uncorrelated with each other.For fGPLS, solutions in parent and offspring Pareto set are correlated.A small distance in solution space corresponds to a small distance in the objective space and a high hypervolume indicator for the offspring Pareto set.A large distance in solution space corresponds to a large distance in objective space and a lower hypervolume indicator for the offspring.In Table 4, we calculate: (i) obj in Table 4(A), and (ii) sol in Table 4(B).As expected, Table