Nicolas Szczepanski
email: nicolas.szczepanski@univ-lille.fr

Lucien Mousin
email: lucien.mousin@univ-catholille.fr

Nadarajen Veerapen
email: nadarajen.veerapen@univ-lille.fr

Laetitia Jourdan
email: laetitia.jourdan@univ-lille.fr

Automatic Configuration of Multi-Thread Local Search: Preliminary Results on Bi-objective TSP

Keywords: parallel computing, automatic configuration, multi-objective local search

In solving combinatorial problems, the advent of multi-core machines has led to the development of new parallel methods offering higher performance as well as better solutions. However, finding the best parallel algorithm on such architectures for a given problem and programming these methods are still challenging tasks today. As a matter of fact, only a few parallel solvers have been designed so far to tackle multi-objective combinatorial optimisation problems. Therefore this paper first proposes a new highly parametric multi-thread and multiobjective local search algorithm dedicated to tackling optimisation problems. Specifically, this new parallel solver incorporates and combines most of methods available in the literature, such as single-walk and multi-walk parallel local search, which can be independent of each other or cooperative by sharing some solutions. In addition, this solver is also capable of making some innovative hybridisations by mixing both single-walk and multiwalk approaches.

The major problem is that it then becomes very difficult, if not impossible, for a human expert to determine the best configurations. This obstacle is due to the fact that the number of parameters in parallel methods is excessively too large. Indeed, all configuration possibilities include the combination of two kinds of parameters. The first ones are the parameters of the sequential algorithms within the parallel solver for multiwalk-based methods. The second ones are those controlling the main parallel components such as hybridisation, diversification or choice of communications. Fortunately, automatic algorithm configuration allows this problem to be taken into account. Thus, we use a configurator especially designed for multi-objective problems called MO-ParamILS in order to expose some best parallel configurations for the bi-objective travelling salesman problem.

I. INTRODUCTION

For NP-Hard multi-objective combinatorial optimisation problems, using an approximation (or incomplete) algorithm, such as metaheuristic, is one of the most popular optimization methods. Multi-objective local search (MOLS) algorithms are metaheuristics able to deal with multi-objective problems by This work has been achieved within the framework of CPER Data project. CPER Data is co-financed by European Union with the financial support of European Regional Development Fund (ERDF), French State and the French Region of Hauts-de-France.

handling the optimisation of several criteria simultaneously. More specifically, MOLS algorithms progressively improve the objective values of one or more solutions of a given optimisation problem through a local search algorithm in order to obtain high quality solutions. In the literature, MOLS algorithms have been used to tackle combinatorial optimisation problems such as the travelling salesman problem [START_REF] Liefooghe | On dominance-based multiobjective local search: design, implementation and experimental analysis on scheduling and traveling salesman problems[END_REF], the quadratic assignment problem [START_REF] Dubois-Lacoste | Anytime pareto local search[END_REF], and the permutation flowshop scheduling problem (PFSP) [START_REF] Blot | Neutral but a Winner! How Neutrality Helps Multiobjective Local Search Algorithms[END_REF], [START_REF] Dubois-Lacoste | Anytime pareto local search[END_REF], [START_REF] Liefooghe | On dominance-based multiobjective local search: design, implementation and experimental analysis on scheduling and traveling salesman problems[END_REF]. Due to the complexity of these problems, improving MOLS algorithms is still a difficult and tedious task.

In order to improve the quality of the solutions within a reasonable amount of time, a basic approach from a hardware perspective consists in running existing solvers on more powerful processing units. This sounds interesting since Moore's Law states that the number of transistors in a processor will double every two years. Nevertheless, Moore's Law may not hold forever. There are physical limits to the ability to continually improve the overall processing power of computers [START_REF] Waldrop | The chips are down for Moore's law[END_REF]. Gordon Moore himself expects it will not hold beyond 2025. Thus, instead of increasing clock rates, processor manufacturers prefer to arrange multiple processors onto the same chip at the expense of some loss of performance (caused by synchronisations, atomic operation, load balancing, and data consistency) and more difficult programming. This is referred to as multi-core architecture.

Leveraging the power of multi-core architectures has had some success for solving hard combinatorial optimisation problems in both mono-objective and multi-objective contexts [START_REF] Codognet | Parallel local search[END_REF]. In the context of this paper, we are interested in parallel methods that focus on the local search algorithm. Parallel local search approaches distinguish between single-walk and multiple-walk methods. On the one hand, single-walk methods (fine-grained parallelism) are limited to the neighborhood of the current solution in the local search algorithm [START_REF] Van Luong | Local search algorithms on graphics processing units. a case study: The permutation perceptron problem[END_REF]. On the other hand, multi-walk methods (coarse-grained parallelism) consist in developing concurrent explorations of the search space by using distinct local search algorithms [START_REF] Yazdani | Flexible jobshop scheduling with parallel variable neighborhood search algorithm[END_REF]. Let us also point out that multi-walk methods can be either independent or cooperative, with some communications between concurrent processes [START_REF] Gabriel Crainic | Cooperative parallel variable neighborhood search for the p-median[END_REF]. In addition, other so-called hybrid methods are based on both single-walk and multi-walk at the same time in order to improve performance [START_REF] Arbelaez | A gpu implementation of parallel constraint-based local search[END_REF]. Nevertheless, there is little work comparing these methods, and it is still difficult to find the best parallel method for a given problem. Moreover, even if it seems natural to design parallel solvers to take advantage of technological advances, only a few parallel solvers have been designed so far for multi-objective problems. In fact, multi-threaded applications are generally hard to implement.

To overcome these disadvantages, we first propose a new highly parametric multi-thread MOLS algorithm which can instantiate the majority of known parallel methods with or without communications. However, as the number possible components and parameters increases, so does the difficult to analyze all possible combinations of strategies in order to find the best configurations.

To address this problem, the recent research field of automatic algorithm configuration (AAC) is an efficient and increasingly popular way to find the best configurations for metaheuristics. Indeed, today some algorithms containing an exponential number of configurations can use an AAC configurator in order to improve solver performances on a homogeneous set of instances representing a combinatorial problem. Thus, the advantages of using an AAC on a highly parametric algorithm have already been exploited in the past for single-objective SLS algorithms [START_REF] Marmion | Automatic design of hybrid stochastic local search algorithms[END_REF] but also on the MOLS algorithms [START_REF] Blot | Automatic Configuration of Multi-Objective Local Search Algorithms for Permutation Problems[END_REF]. In fact, this paper is a continuation of existing work [START_REF] Blot | Automatic Configuration of Multi-Objective Local Search Algorithms for Permutation Problems[END_REF] which aims to move from sequential to parallel computing by exploiting multi-thread architectures.

Thus, our second contribution is to use an automatic algorithm configuration called MO-ParamILS [START_REF] Blot | MO-ParamILS: A Multi-objective Automatic Algorithm Configuration Framework[END_REF] to automatically search for the best configurations of our parallel MOLS algorithm. To our knowledge, there is no work that uses such a tool for a multi-thread MOLS algorithm. On top of that, we propose to use the same AAC configurator to also find good configurations of the sequential MOLS algorithm in order to use them as parameters inside of a portfolio mode within our parallel MOLS algorithm.

The paper is organized as follows. In the next section, basic notions about multi-objective, as well as sequential and parallel, MOLS algorithms are provided. Section III describes the concept of automatic configuration (AC) and enumerates the different strategies implemented in the multi-thread MOLS algorithm. We explain the experimental protocol and present and discuss results in Section IV. We then conclude and provide perspectives in Section V.

II. MULTI-OBJECTIVE LOCAL SEARCH ALGORITHMS

A. Multi-Objective Optimisation (MOO)

Multi-objective optimisation aims to simultaneously improve several criteria (objective functions) that directly affect the solution quality of a given problem. More precisely, a MOO problem consists in determining the set of solutions in which objective functions f i (x) reach their optimal values. We may consider minimization without loss of generally:

argmin x∈D (f 1 (X 1 ), f 2 (X 2 ), . . . , f n (X n )) (1)
where n is the number of objectives (n ≥ 2) and D is the set of feasible solutions where a solution x is represented by the vector of k decision variables x = (x 1 , x 2 , . . . , x k ) of the given problem. Moreover, the sets X n ⊆ x represent several vectors of decision variables which can be different depending on the associated objective function. Note that mixed MOO problems, consisting of both some objective functions to maximize and minimize, can easily be transformed into minimization MOO problems by changing their sign, f ′ i (x) = -f i (x). The concept of Pareto dominance is used to distinguish solutions according to criteria of multi-objective combinatorial problems. A solution s 1 dominates another solution s 2 if, and only if, (i) s 1 is better than or equal to s 2 for all criteria, and (ii) s 1 is strictly better than s 2 for at least one criterion. A set of non-dominated solutions {s 1 , s 2 , . . . , s m }, i.e, in which there is no couple (s i , s j ) (i ̸ = j) such that s i dominates s j is called a Pareto set, a Pareto front, or an archive in the context of multi-objective local search algorithms. Solving a MOO problem consists in finding the Pareto optimal set S * ⊂ D, i.e, the best Pareto set in which there is no other feasible solution x ′ ∈ D that dominates any x ∈ S * .

To assess the Pareto set's quality, various indicators have been proposed. Hypervolume (HV) is one of the most broadly used performance indicators in the literature on multi-objective optimisation [START_REF] Von Lücken | performance metrics in multi-objective optimization[END_REF]. Assuming normalised objective values in [0, 1], unary hypervolume measures the volume between a given Pareto set of solutions and the point (1, 1). While HV is primarily a convergence indicator, it also captures information about the diversity of the set of solutions. Another indicator is a variant of ∆ spread [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: Nsga-ii[END_REF], used to capture the distributional properties of a Pareto set. Given a Pareto set S, ordered according to the first criterion, we define

∆ ′ := |S|-1 i=1 |di-d| (|S|-1)• d ,
where d denotes the average over the Euclidean distances d i for i ∈ [1, |S| -1] between adjacent solutions on the ordered set S. This indicator has to be minimised; it takes small values for large Pareto sets with evenly distributed solutions, and takes values close to or greater than 1 for Pareto sets with few or unevenly distributed solutions.

Since many combinatorial optimisation problems are hard to solve due to their NP-hard complexity, one of efficient solving methods is the use of approximation algorithm in order to obtain high-quality solutions in a reasonable amount of time. These metaheuristics may be classified either as natureinspired (evolutionary, genetic, ant colony, . . . ) algorithms or local search-based (tabu search, iterated local search, variable neighborhood search, . . . ) algorithms.

B. Multi-Objective Local Search (MOLS)

Local search methods explore the search space by iteratively making small changes to a single solution. The procedure generally starts from a good initial solution constructed heuristically using, for example, a hill climbing algorithm. The set of all possible moves for a given solution is called a neighborhood. Next, a better solution than the current one is selected from the neighborhood to become the new current solution. However, the neighborhood may not contain any improving solution. The current solution is then called a local optimum. To remedy this problem, a local search algorithm tries to escape from this local optimum by moving to another area of the search space. Iterative local search (ILS) deals with local optima by restarting from a new random solution or by performing a sufficiently large perturbation on the current solution.

Multi-objective local search (MOLS) algorithms are most often based on Pareto local search (PLS) [START_REF] Paquete | Pareto Local Optimum Sets in the Biobjective Traveling Salesman Problem: An Experimental Study[END_REF], that gradually improves a Pareto set. In the literature, numerous extensions to PLS have emerged, such as the iterated PLS [START_REF] Drugan | Stochastic pareto local search: Pareto neighbourhood exploration and perturbation strategies[END_REF], the stochastic PLS [START_REF] Madalina | Path-guided mutation for stochastic pareto local search algorithms[END_REF], the anytime PLS [START_REF] Dubois-Lacoste | Anytime pareto local search[END_REF] and the dominancebased multi-objective local search (DMLS) [START_REF] Liefooghe | On dominance-based multiobjective local search: design, implementation and experimental analysis on scheduling and traveling salesman problems[END_REF]. As opposed to single-objective local search, improving only one solution, DMLS maintains and improves multiple candidate and nondominated solutions in an archive. The sequential iterated MOLS algorithm used as base in this paper is more recent and proposes a new generalisation of MOLS that exploits many strategies available in the literature [START_REF] Blot | Automatic Configuration of Multi-Objective Local Search Algorithms for Permutation Problems[END_REF]. This algorithm starts by creating an archive of two solutions and then improve them by iteratively executing four distinct components: selection, exploration, archive and perturbation.

At each iteration, the first step is to select the solutions from the archive that will be explored (selection phase). In this phase, the select-size parameter allows the control of the number of solutions to be selected while select-strat represents the way to select them. The latter can be to select all solutions (all) or a certain number of solutions chosen either randomly (rand), or according to the time spent in the archive (newest and oldest). Thereafter the exploration phase seeks new candidate solutions to add to the archive by successively exploring the neighborhoods of the selected solutions. The explor-strat parameter determines which solutions should be added: either the improving ones (all_imp and imp), or the non-dominated ones (ndom) or the both improving and non-dominated ones (all and imp_ndom). Furthermore, the exploration is terminated when a given number explor-size of added solutions is reached except when the value of explor-strat is all_imp or all. In this case, there is no limit to the number of solutions added in the archive. In addition, the explor-ref parameter controls the comparison of neighbors for the improving and non-dominated criteria. This is carried out with respect to either the current solution (sol) or all the solutions in the archive (arch). Once the new solutions have been added to the archive, filtering is performed to keep only the non-dominated ones. Next, if the limit archive-size is exceeded then the archive is reduced by deleting uniformly at random some solutions (archive phase). This ensures the search space is kept within check for the next iteration. The last phase of the iteration, called perturbation, is then applied in order to diversify the search thanks to the parameter perturb-strat. Either a restart is performed by generating new random solutions (restart), or some kick moves are performed on all solutions of the archive (kick_all) or only on some of them (kick) thanks to the parameter perturb-size. Finally, the number of moves in a solution is carried out randomly according to the parameter perturb-strength. In the rest of this paper, we intend to exploit all of these parameters in a new Parallel Multi-Objective Local Search (PMOLS) algorithm.

C. Parallel Multi-Objective Local Search (PMOLS)

Due to the large number of algorithmic components that do not contain any dependent tasks between them (some task A depends on B if it needs the result of B to be executed), local search methods naturally present various forms of parallelism. Parallel Multi-Objective Local Search (PMOLS) algorithms can thus take advantage of numerous opportunities such as fine-grained or coarse-grained parallelism, either with or without communication, . . . . Moreover, PMOLS can be used on various types of parallel architectures such as multi-cores, GPUs or cloud computing. Figure 1 presents the two main methods of parallel local search.

Single-walk methods, sometimes called neighborhood decomposition strategies, consist in using parallelism inside a single search process and exploit the independence between the tasks representing the neighborhood exploration of the local search. Indeed, this exploration can be easily divided into several parts where each part is computed by a thread (Figure 1a). In theory, this method can achieve, in some cases, linear speedup thanks to the independence between these parts. Nevertheless, as for all fine-grained parallel methods, an overload due to thread synchronisation can have a significantly detrimental impact on the speedup. Indeed, parallel processes need to be synchronized in order to choose the most promising neighbor to explore the next solution. Thus, this method can have very poor performance if the calculation cost of evaluations during the exploration is not high enough. However, single-walk methods have proved to be effective in solving different optimisation problems on both multi-thread [START_REF] Codognet | Parallel local search[END_REF] and GPU architectures [START_REF] Van Luong | Gpu computing for parallel local search metaheuristics[END_REF]. Indeed, as the evaluation is the same algorithm but with different data (the solutions) and without dependencies, GPUs are typically dedicated to single-walk methods.

Contrary to single-walk, multi-walk methods are simply the execution of several and distinct local search algorithms, i.e, a portfolio of solvers launched in parallel (Figure 1b). Nevertheless, since the search space is not divided, several processes can explore the same solutions (as illustrated by the threads T 2 and T 3 on the figure). Conversely, too much diversification between the processes prevents local search from focusing on the best solutions. In the field of parallel computing of combinatorial problems, this problem is known as the diversification versus intensification dilemma between search processes. Multi-walk methods can be divided into two categories: independent and cooperative approaches. This last category adds a communication mechanism allowing to exchange information in order to improve either the diversification or the intensification between search processes [START_REF] Codognet | Parallel local search[END_REF].

NEIGHBORHOOD 1 NEIGHBORHOOD 2 SEARCH SPACE T 1 T 2 T 3 T 1 T 2 T 3 (a) Single-walk parallelisms SEARCH SPACE T 1 T 1 T 2 T 2 T 2 T 3 T 3 T 3 T 3 (b) Multi-walk parallelisms
Fig. 1: Different parallel trajectories (arrows) of three threads (T 1 , T 2 and T 3 ) for the single-walk and multiple-walk methods.

The dots are solutions in the search space. Among them, unfilled circles represent the initial solutions of each thread.

In the literature, some cooperative multi-walk approaches are popular. For instance, there are some works on the metaheuristic parallelization which consist in creating multiple solver entities and sharing some information like elite solutions with the master node, which maintains a centralized shared memory for the best solvers [START_REF] Gendreau | Parallel tabu search for real-time vehicle routing and dispatching[END_REF]. Other parallel metaheuristic solvers prefer to exchange their configurations every k iterations and each solver instance decides whether it adopts a received configuration or continues its current search process [START_REF] Cordeau | A parallel iterated tabu search heuristic for vehicle routing problems[END_REF], [START_REF] Machado | Parallel Performance of Declarative Programming using a PGAS Model[END_REF]. More recently, [START_REF] Caniou | Large-scale parallelism for constraint-based local search: the costas array case study[END_REF] presents a parallel adaptive search in which only integer values are exchanged between entities in order to trigger or not restart procedures. Other works are agent-based approaches and are composed of a set of agents that perform specific tasks (search agent, intensifying agent and diversifying agent). Agents are often provided from others domains like machine learning techniques or genetic algorithms [START_REF] Talbi | Cosearch: A parallel cooperative metaheuristic[END_REF]. Finally, ParadisEO-PEO is a framework dedicated to multi-objective constrained combinatorial optimisation and supports different levels of parallel metaheuristics, from neighborhood decomposition (single-walk) to independent and cooperative multiple-walk [START_REF] Talbi | Designing cellular networks using a parallel hybrid metaheuristic on the computational grid[END_REF].

In the paper, we have implemented a hybrid approach mixing both single-walk and multi-walk methods. Nevertheless, to the best of our knowledge, there is less work dealing with this hybridisation. Among them, Arbelaez and Codognet [START_REF] Arbelaez | A gpu implementation of parallel constraint-based local search[END_REF] take advantage of both GPU and CPU architectures at the same time by executing multiple instances of the adaptive search solver and by performing the evaluation of large neighborhoods in parallel.

III. AUTOMATIC CONFIGURATION OF PARALLEL MULTI-OBJECTIVE LOCAL SEARCH A. Multi-Objective Automatic Configuration

The Automatic Algorithm Configuration (AAC) research field provides tools called configurators allowing to determine the best configurations among all valid combinations of parameter values of a target algorithm. Indeed, since this configuration space is exponentially large, in most situations, finding the best configurations manually cannot be accomplished, even by experienced individuals. Thus AAC is becoming an essential ingredient in the design of powerful solvers for challenging optimisation problems. However, most of the existing work on AAC focuses on configuration procedures that optimize a sequential and single-objective algorithm. In contrast, in this paper, we deal with a more challenging AAC problem which aims to find the best configurations of a target algorithm combining both parallel computing and multi-objective optimisation.

On the one hand, to optimize the target algorithm, most AAC configurators in the literature (irace [START_REF] López-Ibáñez | The irace package: Iterated racing for automatic algorithm configuration[END_REF], ParamILS [START_REF] Hutter | Paramils: An automatic algorithm configuration framework[END_REF], or SMAC [START_REF] Hutter | Sequential model-based optimization for general algorithm configuration[END_REF]) use a single indicator usually representing the resolution time or the solution quality. But more recently, as part of the work of Blot et al. [START_REF] Blot | Automatic Configuration of Multi-Objective Local Search Algorithms for Permutation Problems[END_REF], multiple performance indicators dedicated to multi-objective optimisation problems have shown their effectiveness through the MO-ParamILS configurator [START_REF] Blot | MO-ParamILS: A Multi-objective Automatic Algorithm Configuration Framework[END_REF]. Moreover, the authors of [START_REF] Blot | Automatically Configuring Multiobjective Local Search Using Multi-objective Optimisation[END_REF] have achieved better performance with MO-ParamILS rather than with a single-objective configurator that performs an aggregation of two objectives. Let us also note that another tool called SPRINT-Race, based on multi-objective model racing, was used to identify the Pareto optimal parameter settings of Ant Colony Optimization algorithms [START_REF] Zhang | Anagnostopoulos. Sprint multi-objective model racing[END_REF].

On the other hand, AAC development and applications remain largely focused on sequential rather than parallel approaches. In parallel computing, some AAC configurators has been employed in order to determine the best parallel configurations [START_REF] Lindauer | Automatic construction of parallel portfolios via algorithm configuration[END_REF]. More precisely, Automatic Construction of Parallel Portfolios (ACPP) involves optimizing some sequential parameters of each solver into the portfolio. For instance, Lindauer et al. [START_REF] Lindauer | Automatic construction of parallel portfolios via algorithm configuration[END_REF] expose parameterized SAT solvers that produce significantly better SAT solvers than state-of-the-art parallel solvers built by human experts.

Finally, we have not encountered any work on the application of a configurator to a target algorithm that is both parallel and multi-objective, and this paper addresses this topic.

B. Parameters of Parallel MOLS

The parameters are classified into two parts: parameters coming from the sequential MOLS algorithm and those dedicated to the parallel MOLS program. The first parameters are those controlling the search of a single solver, i.e. a single MOLS algorithm. Among them, we find the totality of parameters used by Blot et al. [START_REF] Blot | Automatic Configuration of Multi-Objective Local Search Algorithms for Permutation Problems[END_REF] presented in section II-B such as explor-strat or archive-size. In addition, the other parameters deal with the main parallel components such as the single-walk, the multi-walk, the hybridisation, or the choice of communications.

Nevertheless, some parameters depend on other parameters coming from either sequential or parallel MOLS algorithms. For example, the parameter explor-size has to be used only when the value of explor-strat is imp, imp_ndom or ndom. In addition, the parallel single-walk method is inherently incompatible with certain MOLS methods. In fact, in order to have enough tasks with sufficient computational load in a single-walk method, the values of parameters select-size and explor-size have to be large enough.

To alleviate this problem, we carry out two different experiments: one for the single-walk and hybrid methods, and one for the multi-walk methods. Thus, each of these experiments has different parameter values as exposed in the table I. Respectively according to the single-walk and hybrid experiment and the multi-walk experiment, the total number of configurations is to 5460 and 10920 for the sequential part and is to 7500 and 7560 for the parallel part.

The three columns of Table I represent, respectively, the name of the parameter and the parameter values, first for the single-walk and the hybrid experiment, and then for the multiwalk experiment. The lines of the table summarize firstly in a first group the sequential parameters used, whereas the second group exposes the new parameters implemented in our parallel MOLS algorithm.

hybridisation. This parameter is composed of two numbers: the first is defined as the number of threads representing the portfolio of the multi-walk method, noted n, where each thread is a MOLS algorithm, while the second is the number of threads dedicated to the single-walk methods for each MOLS algorithm of the portfolio, noted m. Thus the total number of threads used in our PMOLS algorithm is the product of these two numbers n × m. For example, the values ⟨4, 5⟩ represent a PMOLS algorithm of n = 4 MOLS algorithms in parallel, and each of them divides the computational load of its neighborhood exploration across m = 5 threads thanks to the single-walk method, making a total of 20 threads. Two extreme cases of this parameter should also be noted: the pure single-walk method and the pure multi-walk method without hybridisation and respectively represented by the values ⟨1, 20⟩ and ⟨20, 1⟩. configuration-n. This set of parameters is used to choose the values of sequential parameters (i.e. configurations) for each of n MOLS algorithms in the portfolio among x available sequential configurations. In the example of Table I and in our experiments, x = 5 in order to limit the exponential explosion of an excessively large number of sequential configurations. By definition, these parameters depend on the value of parameter hybridisation, i.e., configuration-n is taken into account (i.e. used as parameter) only if the value of hybridisation is equal to ⟨n, m⟩ (i.e, contains n MOLS algorithms in the portfolio). For example, when hybridisation is ⟨10, 2⟩, only configuration-10 is used as parameter, and other configuration-n such as n ̸ = 10 are not possible parameters. A possible value of the parameter configuration-n is a set of distinct numbers between 1 and 5 where each number represents a sequential configuration among the x available sequential configurations. Moreover, the values of configuration-n are all combinations of size 1, 2, 3, 4 or 5 of the five available sequential configurations, respectively for n equal to 1, 2, 4, 10 or 20. For instance, the values of configuration-4 are of size 3 and some of them can be, for example, ⟨1, 2, 3⟩ or ⟨2, 4, 5⟩.

In addition, we remove also those that are equivalent (⟨1, 2, 3⟩ and ⟨2, 3, 1⟩). This allows to browse through a wide range of sequential configuration combinations for the portfolio without adding unnecessary or equivalent configurations. When the number of MOLS algorithms n is greater than the number of available sequential configurations x, then our parallel algorithm repeats the distribution in the order of configurations given by parameter values of configuration-n. For example, for hybridisation ⟨10, 2⟩, the chosen parameter is configuration-10, and if this parameter is ⟨1, 2, 2, 5⟩, then the sequential configurations of 10 MOLS algorithms are the configurations 1, 2, 2, 5, 1, 2, 2, 5, 1 and 2.

diversification. This parameter specifies the mechanism of the diversification of each MOLS in the portfolio. Either only their initial solutions are distinct (solution) or only their neighborhood explorations are performed in different orders (neighborhood) or both.

communication and communication-send. These parameters control the amount of exchanged solutions at the end of each neighborhood exploration. The communication can be turned off thanks to the value none of communication. In the other cases, solutions are received either by one or all MOLS algorithms or half of them. In addition, the parameter communication-send is the number of solutions to send to others. Note also that the solutions sent are among the best from the current archive according to the Pareto-dominance.

IV. EXPERIMENTS

A. Case Study

The Symmetric Travelling Salesman Problem (TSP) can be modelled as an undirected weighted and complete graph G where the vertices denote the cities whereas the edges (pairs of non-ordered vertices) correspond to distances between cities. Given a TSP instance G, the goal is to determine the shortest possible route that visits each city (vertex) and returns to the origin city, i.e, finding a minimum-weight Hamiltonian cycle in G. The decision version is an NP-complete problem and a possible solution (not necessarily the best route) is represented by a permutation of cities, i.e., a Hamiltonian G-cycle.

The experiments in this paper focus on the bi-objective symmetric TSP (bTSP) problem, which is defined as the TSP but with the peculiarity of having two weights to minimize per edge instead of one. In other words, the aim is to minimize the total distance covered by a round trip according to each of the two objectives. Note also that these two objectives are not correlated since they were computed by combining two independently generated distance matrices computed using Euclidean distance between cities randomly placed on a twodimension grid.

In order to fairly test the final parallel configurations found by the configurator on other different instances, we have two disjoint sets of distinct bTSP instances: training instances and test instances. Furthermore in each of these sets, we classify the instances according to the number of cities which are: 100, 300 and 500 cities.

The training instances are new instances generated by the original generator from the DIMACS challenge and contain 30 instances for each number of cities representing a total of 435 possible combinations for the bTSP problem per city size. In contrast, test instances taken from the literature are used to measure the performance of bTSP algorithms [START_REF] Liefooghe | On dominance-based multiobjective local search: design, implementation and experimental analysis on scheduling and traveling salesman problems[END_REF], [START_REF] Paquete | Pareto Local Optimum Sets in the Biobjective Traveling Salesman Problem: An Experimental Study[END_REF]. This set consists of 6 instances for each number of cities totaling 15 different pairwise independent combinations per city size.

B. Protocol

In order to configure our parallel MOLS without aggregating hypervolume and spread into a single performance indicator, we have made the choice to use the recent multiobjective AAC configurator MO-ParamILS, itself based on a MOLS algorithm. Thanks to this feature, MO-ParamILS is, in effect, fully capable of exploiting the multi-objective aspect. The usual MO-ParamILS configuration protocol uses three phases: training, validation and test.

In the training phase, we run the MO-ParamILS configurator on the training instances many times with different random seeds in order to find diverse and various Pareto fronts representing some distinct good configurations. In order to complete this phase within a reasonable period of time, note that MO-ParamILS does not test the configurations on all instances but only on a subset of them. Next, to reduce the number of configurations by taking only the best of them, a second phase allowing a fairer comparison of configurations is necessary. This is achieved by the validation phase which consists in executing the configurations found by the training phase on the complete set of training instances. The result thus obtained is in the form of a single Pareto front of some non-dominated configurations. To finish, the Pareto set of configurations of the validation phase is evaluated on the test instances in order to check the quality of these configurations on other different instances.

However, this protocol usually used for sequential algorithms has to be adapted to the parallel features of our PMOLS program. We have therefore chosen to divide the protocol into two steps: one for the sequential MOLS and one for the parallel MOLS.

The first step is composed of a training phase followed by a validation phase and aims to find five good configurations of the sequential MOLS algorithm (among the parameters of the first part of the table I) intended for the parameter configuration-n of the parallel MOLS algorithm. These five configurations are first chosen in the Pareto front of the validation phase and then among other dominated configurations if this Pareto front is out of stock. The second step consists of finding the best configurations of the parallel MOLS algorithm by using the five good configurations of the previous step. However, due to limited resources, the training and validation phases have been revised by dividing by four the number of training instances representing a total of 30 instances.

The experiments have been conducted on a cluster of four computers each containing 20 Intel XEON cores running at 2.2 GHz and with 500 GiB RAM. The cluster is equipped with an Ethernet controller at 1 GiB/s. The sequential MOLS algorithms have been rewritten from [START_REF] Blot | Automatic Configuration of Multi-Objective Local Search Algorithms for Permutation Problems[END_REF] into a new library called MH-builder. This library now takes advantage of the object-oriented programming in order to be fully modular and thus help in the building of many MOLS algorithms. The parallel methods and the communications have been implemented thanks to the pFactory library available at https://github.com/crillab/pfactory [START_REF] Audemard | Pfactory: A generic library for designing parallel solvers[END_REF].

For each MOLS or PMOLS algorithm, the base times are 20, 60 and 100 seconds, respectively, for 100, 300 and 500 cities. Note that the training phases are consist of 10 MO-ParamILS runs achieving each 100 MOLS algorithm. The total time of the validation phases depends of the number of configurations found in the training phases. Thus, for 500 cities and the single walk and hybrid experiment, the sequential step has lasted about 9 days: 1 day for the training phase (1000 × 100 seconds) and 8 days for the validation phase (15 configurations ×455 instances ×100 seconds). For the sake of comparison, an exhaustive approach on the training instances will last 15.06 years (10920 configurations ×435 instances ×100 seconds). Next, the parallel step has taken 2 days: about 1 day for the training phase (1000 × 100 seconds), 1 day for the validation phase [START_REF] Marmion | Automatic design of hybrid stochastic local search algorithms[END_REF] for the test phase. Once again, note that an exhaustive approach on the training instances would have lasted 0.71 years (7500 configurations ×30 instances ×100 seconds). In the same way, each other experiment takes between 4 and 15 days of computation.

C. Results and Discussion

Let us recall that we have performed two experiments, the hybrid one, and the multi-walk one. Each of these experiments follows the protocol previously indicated, consisting of a sequential step and a parallel step. Figure 2 shows the results of the test phases of the parallel step (the triangles and squares) and, for the sake of comparison, the results of the five sequential configurations (the circles) used in the associated parallel methods on the same test instances with identical run times (20, 60 and 100 seconds, respectively, for 100, 300 and 500 cities). Each point represents the averages of 1-hypervolume and spread values of each configuration for all test instances. Thus, the closer the points are to the origin (with small values), the better their hypervolume and spread. Note that the comparisons between the sequential and parallel methods are performed only w.r.t. the hypervolume and spread values, and not on the speedup in terms of time.

On the one hand, the multi-walk experiment is represented by the first three sub-figures 2a, 2b and 2c. We can see that the parallel methods are better than the sequential ones for both the hypervolume and spread values. On the other hand, the hybrid experiment (sub-figures 2d, 2e and 2f) exposes greater discrepancies between sequential and parallel methods, notably for the hypervolume values. Furthermore, for all experiments, the larger the city size, the greater the gap. For the single-walk methods, this is due to the fact that they are more efficient when the computation load of the function evaluation is higher, i.e, for 500 cities. More generally, the parallel performances are also less strong for 100 cities because sequential methods are already very close to optimality. By comparing both experiments, we can observe that the multi-walk methods are the best, even though some hybrid configurations are similar in terms of hypervolume for 300 cities. We conclude that both hybrid and multi-walk methods are both useful and powerful.

Concerning sequential parameters, the imp-ndom, ndom and all-imp strategies with different kick sizes and select strategies are the most commonly used. The archive size is often set to 1000. We also found explor-ref, which controls what to compare neighbors against when evaluating improvement and non-domination, to be equal to either the archive or the solution. For parallel parameters, communications are useful because they are used in most of the best parallel configurations by often sending one solution to some other solvers of the portfolio. The parameter diversification is never set to solution and often set to both (i.e. a combination of distinct initial solutions and performing neighborhood exploration in different orders). Concerning hybridisation, the best parallel configurations use a portfolio of 4 solvers where each solver performs a singlewalk of 5 threads for 100 cities (value 4_5). For the other city sizes, the parameter hybridisation is often also equal to 4_5 but sometimes to 10_2 or 2_10. We can see also two pure single-walk methods (1_20) in the configurations of 500 cities. Therefore, this shows that all parallel components may be useful, within a given context, and they can improve the parallel MOLS algorithm performances when configured correctly.

V. CONCLUSION AND PERSPECTIVE

This paper deals with, and obtains encouraging preliminary results, on a problem that had never been tackled until now: how to find the best configurations of a local search algorithm that is both parallel and multi-objective. As a first step, we have implemented a new highly parametric multithread MOLS dedicated to tackling optimisation problems. In addition, this solver is also capable of making some innovative hybridisations by mixing both single-walk and cooperative multi-walk approaches. Next we have exposed some best parallel configurations by using MO-ParamILS thanks to two distinct experiments. In the near future, we plan to build upon this work and explore various perspectives. Firstly, we will test other sequential configurations for the portfolio of the parallel MOLS algorithm because, as we can see from our results, it is not necessarily the best sequential configurations that form the best parallel solver. Secondly, we are also going to add the possibility for solvers to communicate their search space in order to improve their diversification potential. Finally, the performance of our parallel MOLS algorithm should also be investigated on other problems, such as the rule mining problem [START_REF] Sara | Multi-objective automatic algorithm configuration for the classification problem of imbalanced data[END_REF]. Indeed, because the evaluation steps of this problem are more expensive, we believe that hybrid methods would be more suitable.

Fig. 2 :

 2 Fig.2: Results of the single-walk and hybrid experiment and the multi-walk experiment on 100, 300 and 500 cities.

TABLE I :

 I Sequential and parallel parameters.

  configurations ×30 instances ×100 seconds) and about 1 hour (2 configurations ×15 instances ×100 seconds)

	1					1.2
														1.8						PARETO FRONT
						1.1														OTHERS
	0.9													1.6						SEQUENTIALS
						1
	0.8													1.4
	0.7 ∆ SPREAD					0.7 0.8 0.9 ∆ SPREAD							1 1.2 ∆ SPREAD
	0.6
						0.6							0.8
	0.5
						0.5							0.6
	0.4					0.4							0.4
	0.1	0.12	0.14	0.16	0.18	0.2	0.08	0.09	0.1	0.11	0.12	0.13	0.14	0.05	0.06	0.07	0.08	0.09	0.1	0.11	0.12	0.13	0.14
			1 -HYPERVOLUME					1 -HYPERVOLUME						1 -HYPERVOLUME
		(a) Multi-walk -100 cities				(b) Multi-walk -300 cities					(c) Multi-walk -500 cities
	1					1.2
														1.8
						1.1
	0.9
														1.6
						1
	0.8													1.4
	0.7 ∆ SPREAD					0.7 0.8 0.9 ∆ SPREAD							1 1.2 ∆ SPREAD
	0.6
						0.6							0.8
	0.5
						0.5							0.6
	0.4					0.4							0.4
	0.1	0.12	0.14	0.16	0.18	0.2	0.08	0.09	0.1	0.11	0.12	0.13	0.14	0.05	0.06	0.07	0.08	0.09	0.1	0.11	0.12	0.13	0.14
			1 -HYPERVOLUME					1 -HYPERVOLUME						1 -HYPERVOLUME
		(d) Hybridisation -100 cities			(e) Hybridisation -300 cities				(f) Hybridisation -500 cities


Kenneth N Reid

Jingpeng Li

Alexander Brownlee

Mathias Kern
email: mathias.kern@bt.com

Gilbert Owusu
email: gilbert.owusu@bt.com

A Hybrid Metaheuristic Approach to a Real World Employee Scheduling Problem

Keywords: CCS CONCEPTS, Computing methodologies → Randomized search, • Applied computing → Industry and manufacturing, Evolutionary Ruin and Stochastic Recreate, Metaheuristics, Employee Scheduling, Variable Neighbourhood Search

Employee scheduling problems are of critical importance to large businesses. These problems are hard to solve due to large numbers of conflicting constraints. While many approaches address a subset of these constraints, there is no single approach for simultaneously addressing all of them. We hybridise 'Evolutionary Ruin & Stochastic Recreate' and 'Variable Neighbourhood Search' metaheuristics to solve a real world instance of the employee scheduling problem to near optimality. We compare this with Simulated Annealing, exploring the algorithm configuration space using the irace software package to ensure fair comparison. The hybrid algorithm generates schedules that reduce unmet demand by over 28% compared to the baseline. All data used, where possible, is either directly from the real world engineer scheduling operation of around 25,000 employees, or synthesised from a related distribution where data is unavailable.

INTRODUCTION

In this paper we present a solution to a real world employee scheduling problem, using a hybridization of Evolutionary Ruin & Stochastic Recreate (ER&SR) with Variable Neighbourhood Search (VNS) and an evolutionary Monte-Carlo acceptance criterion (EMCAC).

The problem specification and associated data were provided by a FSO (Field Service Operations company) and specific problem details are intentionally omitted. This collaboration with the FSO ensured the precise design of constraints and fitness criteria. Since this paper deals with a real world problem, the algorithm's effectiveness is empirically grounded.

Employee scheduling is a problem where employees are assigned to one or more shifts in order to meet a number of objectives (e.g. improving employee efficiency in their roles and meeting demand). Employee scheduling is a top concern for management in improving organisational efficiency [START_REF] Enz | Key issues of concern in the lodging industry: what worries managers[END_REF]: failure to optimize shift allocation for employees can mean a loss in profits, customer satisfaction and failure to meet demand.

ER&SR and VNS have been successful when used individually to produce satisfiable schedules in the domain of personalised scheduling [START_REF] Reid | Variable neighbourhood search: A case study for a highly-constrained workforce scheduling problem[END_REF][START_REF] Reid | Shift scheduling and employee rostering: An evolutionary ruin & recreate solution[END_REF]. Our proposed approach in this paper compared three combinations of VNS and ER&SR, and uses EMCAC to control the end criterion. The three approaches are as follows: ER&SR followed sequentially by VNS (i.e. each algorithm runs only once), an ER&SR and VNS hybrid (VNS within ER&SR) and a third approach which combines an internal VNS and a sequential VNS with ER&SR. This is detailed further in Section 4.

The irace package [START_REF] López-Ibáñez | The irace package: Iterated racing for automatic algorithm configuration[END_REF] was used to automatically configure the parameters for the various ER&SR and VNS settings. This is to ensure comparisons were fair with the same time budgets for each configuration, excluding the handwritten baseline.

The contributions of this paper are as follows:

(1) Propose the new hybridised algorithm of ER&SR and VNS.

(2) Analysis of this algorithm's application to a real world employee scheduling problem. (3) A systematic parameter search by irace to find a suitable configuration. (4) A comparison of the hybrid algorithm with four other configurations and the improvement relative to the baseline.

The problem domain is similar to that considered in a previous work [START_REF] Reid | Shift scheduling and employee rostering: An evolutionary ruin & recreate solution[END_REF], in which both the shift scheduling and employee rostering problems with personalised schedules are tackled. A personalised schedule is unlike a Roster Pattern (RP) based approach, in that shift start times, end times, shift lengths are variable day to day. In this paper we utilise a RP based approach where a set of shift patterns are input parameters and cannot be modified, instead the selected pattern per employee, and start week, are modified.

The base problem requires the assignment of a skill for each employee, for each shift. The additional constraints are expressed in the form of employee RPs, which requires the additional weekwise assignment of a RP an employee. The search space cannot be exhaustively explored in polynomial time and similar employee scheduling problems have been considered NP-Hard [START_REF] Heimerl | Scheduling and staffing multiple projects with a multiskilled workforce[END_REF]. Metaheuristic solutions were used to solve this problem to acceptable quality in reasonable business time frames, as defined by the FSO. We propose a hybrid solution in anticipation of jointly benefiting from the exploratory nature of ER&SR and the exploitative local search of VNS.

Employee rostering is often used interchangeably with employee scheduling and workforce rostering. In this paper we use the term 'employee rostering' to mean the allocation of an employee to a shift or set of shifts which are preset and immutable. In this paper we also use the term 'employee' and 'engineer' interchangeably as the real world problem pertains to engineers.

This paper is structured as follows: Section 1 completes with a discussion on related works in subsection 2. In Section 3 the problem description is provided as are the benefits of this approach. Section 4 details the solution implementation, using the theoretical framework provided by Li et al [START_REF] Li | Search with evolutionary ruin and stochastic rebuild: A theoretic framework and a case study on exam timetabling[END_REF] as a basis, and the unique elements of our implementation are described. The results of this approach are analysed in Section 5 and Section 6 presents conclusions.

RELATED WORK

In the field of employee rostering there is much motivation for obtaining high quality employee rosters. Managing the shifts of a large number of employees is a complex task -optimisation of this process is a necessity for efficient use of resources. The field of employee rostering has been extensively surveyed: one of the first surveys of the state of the art was produced by Baker [START_REF] Baker | Workforce allocation in cyclical scheduling problems: A survey[END_REF] in 1976. In 2004 Ernst et al produced a review of applications, methods and models [START_REF] Ernst | Staff scheduling and rostering: A review of applications, methods and models[END_REF], which not only successfully highlighted some of the most influential research in the field, but also correctly predicted future trends. Such trends include employee scheduling in airlines, which has since been further examined by several reports [START_REF] Benlic | Breakout local search for the multiobjective gate allocation problem[END_REF] [4], and the need for more robust frameworks to tackle the complexity of these problems.

The extensive 2013 review 'Personnel Scheduling: A Literature Review' by Van Den Bergh et al [START_REF] Van Den Bergh | Personnel scheduling: A literature review[END_REF], gives a thorough exploration of the field, in which papers are categorised by criteria such as personnel characteristics, decision types, shift flexibility, coverage constraints (both hard and soft), whether or not under-and overstaffing is allowed, by skills, and more. There have been several other notable domain specific literature reviews, such as home care routing and scheduling [START_REF] Fikar | Home health care routing and scheduling: A review[END_REF] and physician scheduling [START_REF] Erhard | State of the art in physician scheduling[END_REF].

The term 'metaheuristic' was first coined by Glover [START_REF] Glover | Future paths for integer programming and links to artificial intelligence[END_REF] in the celebrated paper introducing Tabu search. Metaheuristics are highlevel procedures for exploring complex solution spaces, designed to solve the problem of entrapment within local optima [START_REF] Gendreau | Handbook of metaheuristics[END_REF]. Metaheuristics are 'higher level' heuristics in that they define generic solution mechanisms: for example, a metaheuristic framework for Simulated Annealing can be customised for different problem domains (e.g. Travelling Salesperson, Vehicle Routing, Scheduling etc.) via the provision of domain-specific heuristics for solution neighborhood and solution quality. Popular metaheuristics include Simulated Annealing [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], VNS [START_REF] Mladenović | Variable neighborhood search[END_REF], tabu search [START_REF] Glover | Future paths for integer programming and links to artificial intelligence[END_REF], genetic algorithms, memetic algorithms [START_REF] El-Yaakoubi | Tabu search and memetic algorithms for a real scheduling and routing problem[END_REF] and many others.

Scheduling vast numbers of employees while taking demand, location, skill sets, preferences and contract types into consideration is a complex task which cannot be completed manually for a realworld workforce. It is often impossible to exhaustively search for an optimal solution in a timely fashion [START_REF] Burke | Search methodologies[END_REF]. As such other approaches have been developed by operational engineers, computer scientists and employee managers to generate optimal or close-to-optimal solutions, e.g. using evolutionary computation approaches [START_REF] Pinedo | Scheduling in the service industries: An overview[END_REF]. Improving labour costs by only a few percent could prove very beneficial, since this is often the main expenditure for companies [START_REF] Van Den Bergh | Personnel scheduling: A literature review[END_REF].

Rather than the expensive (and often impossible) alternative of analysing every potential roster, an evolutionary algorithm can sample the solution space for a result which is often at least close-tooptimal, within a reasonable time. VNS is a methodology which is also commonly used in optimisation problems. In a nurse rostering problem example [START_REF] Rahimian | A hybrid integer programming and variable neighbourhood search algorithm to solve nurse rostering problems[END_REF], VNS is hybridized with integer programming. In this approach a greedy heuristic is employed to tackle the initial problem and create a base solution, which is successfully improved upon iteratively with VNS deep-embedded with integer programming. Bruecker et al in a 2015 state-of-the-art review [START_REF] De Bruecker | Workforce planning incorporating skills: State of the art[END_REF] on workforce planning incorporating skills concluded that more real world application is essential to advance the field. The authors also found that mixed integer linear programming appears to be the most popular mathematical programming approach, while custom approaches are the most popular in heuristic methodologies.

PROBLEM DESCRIPTION

In this section we describe the problem as defined by collaboration with the FSO. We then state the Hard Constraints (HCs) and Soft Constraints (SCs), giving context for a quality measure on feasible solutions.

Similarly complex employee scheduling problems in the literature are NP-Hard [START_REF] Heimerl | Scheduling and staffing multiple projects with a multiskilled workforce[END_REF], so the goal of this research is to find a solution meeting all HCs and as many SCs as possible within a time span deemed acceptable by the FSO. This time span is within 2 minutes on a medium specification workbook issued to employees within the FSO. For testing purposes run time is greatly increased, as described in Section 5.

The problem and data are both real world and provided by an FSO company. The information in the paper has been anonymised and generalised for security and data protection reasons.

The problem we tackle in this paper schedules up to 200 employees, varying depending on input data region, from a pool of around 25,000 employees across the UK. There are therefore hundreds of possible sets of input data from various localities. There are a variety of local differences which change frequently due to variation in demand forecasting or problem region. Differences can include different temporal patterns in demand, different skills required, different skills available from engineers and a different number of engineers available. This problem has the following unique combination of traits as defined in collaboration with the FSO:

(1) Due to the real-world data being provided by a UK based FSO, UK Employment Law must be adhered to. This includes maximum working hours allowed per week, number of rest hours in between shifts, etc. (2) The employees described in the data have a variety of skills.

Demand is provided for the problem in skill-based manner. Demand may fluctuate greatly day to day, week to week. The objective of meeting demand includes not simply providing manpower to cover number of employees required per shift, but skills required per shift, and per day, while still meeting HCs. The term "shift" refers to a time period within some 24 hour period. A shift can have a single employee allocated to it. Shifts have the attributes of start time, end time, length, demand covered (in terms of an integer array with potential values of 0 or 1 per hour) and skill covered.

A solution is deemed feasible if the associated schedule does not violate any legal, contractual or ethical guideline defined by the FSO, each of which is represented by a HC. Such guidelines prevent, for example, working illegal number of hours per week or working shifts that violate an employee's contract. As this problem works from RPs provided by the FSO, most HCs that would normally be considered in employee scheduling are unnecessary, as every RP meets all legal and contractual obligations. HCs therefore cannot be violated by nature of the encoding of the problem: the algorithm will only generate and test solutions which meet all HCs. Thus the concern of this paper is only to minimise the cost of SCs, or inversely improve the quality (fitness) of schedules.

• HC1: Employees must have exactly one RP.

• HC2: An optionally specified maximum number of RP changes can occur. • HC3: An optionally specified maximum number of starting week changes can occur. • HC4: If RP preferences are provided, only RPs preferred by each employee can be explored. HC1 is hardcoded, an employee cannot be given more than or less than one RP outside of a single iteration of ER&SR or VNS (meaning that employees do have RPs removed, but they are then given another RP, or potentially the same one again).

HC2 and HC3 can optionally have the maximum number of changes set to no maximum, zero (meaning no changes of this type allowed) or a specific number. An example could be maximum of 10 RP changes allowed. When 10 have occurred, only those modified RPs can be changed. If one of these is returned to the original pattern, then the current number of changes is reduced to 9, and another employee's pattern could potentially be changed. In this paper HC2 and HC3 are set to infinite, meaning they are always satisfied.

HC4 can be set to on or off. If set to off, any employee could potentially be given any RP. This is generally not realistic in many regions as certain RPs are created for special circumstances, specific contractual requirements, etc. It would also mean employees could be given the incorrect number of working hours per week. However, if the available RPs are all verified to be interchangeable, then this is not an issue.

For most regions, employees have a number of preferred RPs, and HC4 would be toggled on. The preferred patterns will always meet contractual obligations. These patterns are agreed upon by both employee and regional managers. Every employee has at least one preferred RP, which is their current pattern. There is no maximum number of preferred RPs, other than the number available in the pool, which is dependent on how many have been created by managers in the past. In the region used for this paper there are around 400 potential RPs employees could have in their preferred patterns list.

There are no hard constraints regarding starting week optimisation, so any employee can still have rotas modified, as well as skills worked every shift even if they cannot change RP directly. For example, if an employee has a rotating rest day off (Monday on week 1, Tuesday on week 2, ..., Saturday on week 6) then changing the starting week can drastically affect potential skills available on week days across a number of months.

A solution's quality is judged with SCs. SCs consider elements that are important for the quality of life for engineers, for meeting demand and for meeting customer satisfaction. These are important, but do not break any laws or contractual requirements if not met. The following SCs are considered to judge the quality of a solution:

• SC1: Shifts should cover the maximum demand requirement they possibly can. • SC2: Employees should be assigned to shifts that require a skill listed in their preferences. • SC3: Allocations of employees to shifts should favour employee skill preferences in descending order. • SC4: Variation in the percentage of demand met per hour, day and skill should be minimised.

To be clear: all engineers have the skills required for all work, however employees have a list of "preferred" skills. The reasons for the skills being preferred are varied, from personal enjoyment, to proficiency or managerial preference. As such, SC2 measures whether employees are working a skill listend, SC3 judges quality based on level of preference.

Combined weighted constraint violations are used to judge solution fitness. Specific weightings are provided by the FSO; in this paper all weightings are equal (25% per SC). In real world usage the weightings are dependent on locale specific priorities, equal weighting was judged to be acceptable for testing within this paper. As an example of why this is useful, if demand is not being met to a satisfiable level, SC1 could be set to 100%, and the others set to 0%.

IMPLEMENTATION 4.1 Evolutionary Ruin & Stochastic Recreate

ER&SR [START_REF] Li | Search with evolutionary ruin and stochastic rebuild: A theoretic framework and a case study on exam timetabling[END_REF] is a variant of the original Ruin & Recreate method [START_REF] Schrimpf | Record breaking optimization results using the ruin and recreate principle[END_REF] by the addition of evolutionary features, selection and mutation, as the ruining strategies. In this section we describe the differences from previous work on ER&SR. ER&SR was used to solve not only an employee rostering but also a shift scheduling problem, without VNS in previous works [START_REF] Reid | Shift scheduling and employee rostering: An evolutionary ruin & recreate solution[END_REF], however in this instance ER&SR is used in hybrid with VNS to solve a real world employee rostering problem, not a shift scheduling problem.

Our implementation was designed to be modular, with a VNS component which can be called at any stage of the algorithm. This allows us to consider a variety of approaches and how to best exploit the local search advantages from VNS in conjunction with the more exploratory ER&SR. As such we implemented the following hybridised algorithms:

Listing 1: ER&SR Pseudocode (i) ER&SR with VNS running sequentially after ER&SR stopping criterion (referred to as 'ER&SR with VNS Sequential') (ii) ER&SR with VNS running after accepting a modified schedule (referred to as 'ER&SR with VNS Internal') (iii) ER&SR with both the above (referred to as 'ER&SR with both VNS internal and sequential') Pseudocode of ER&SR with VNS can be seen in Listing 1. Lines 19 and 25 reference VNS, of which is described more fully in subsection 4.2, with VNS pseudocode available in Listing 2.

Our previous paper [START_REF] Reid | Shift scheduling and employee rostering: An evolutionary ruin & recreate solution[END_REF] worked with a personalised scheduling problem as opposed to RPs, we introduced a number of changes. The same phases exist and perform the same functions with modifications, and the VNS hybridisation has since been implemented.

(1) The initial set up phase initialises the software and reads the current hand-written solution in place by the FSO. (2) The EMCAC control mechanism is initialised with the parameterised inputs, controlling the number of times ER&SR loops (but not the sequential VNS element). (3) Phase 1 -solution decomposition no longer considers HCs. This is due to RPs which ensure contractual and legal requirements cannot be breached. More precisely, the HCs considered in this paper can be temporarily breached during the evolutionary ruin phase, but are all then restored during stochastic recreate. SCs are considered in the solution decomposition, and the fitness evaluations can be completed by considering only the potential changes from the current state (referred to by [START_REF] Li | Search with evolutionary ruin and stochastic rebuild: A theoretic framework and a case study on exam timetabling[END_REF] as componential fitness). This is necessary as all four SCs consider the percentage of all individual shifts meeting critera. The solution decomposition serves the purpose of analysing fitness during ER&SR run-time as well as after VNS swaps have occurred to decide whether to keep or revert changes. (4) Phase 2 -evolutionary ruin is mostly unchanged. Instead of ruining individual shifts, both the selection and mutation elements now ruin employee RPs, by giving a null allocation.

Selection will destroy elements with the following probability:

P(s) = r * x 1 ,
where P(s) is the probability of acceptance, r is a random number between 0 and 1, and x 1 is a parameter to limit the stochasticity of selection. If P is higher than the fitness of the employee allocated to their current RP, then the employee has their RP removed for a new one. After selection is complete, mutation considers all non-affected employees, and randomly removes employees from RP. This adds an additional level of stochasticity to the search. Fitness is not considered during mutation.

Mutation can be overly destructive and hinder exploitation, so careful parameterisation is important. The chance of ruination is calculated as follows:

P(b) = r * x 2 ,
where P(b) is the probability, r is a random number between 0 and 1, x 2 is the parameter used to increase or decrease probability. All of this phase is likely to break all hard constraints before the schedule is reconstructed and made feasible in the next phase. ( 5) Phase 3 -stochastic recreate no longer finds random shifts to allocate employees to, but instead allocates employees to random RPs and starting weeks. Employees to be allocated a new shift pattern are chosen randomly. This is so the search space can be explored if a maximum bound is set via HC2 or HC3. These HCs are considered with every RP allocation to ensure the maximum bound is not breached. If HC4 is toggled on, only RPs the employee has selected may be considered for allocation. Otherwise, any RP can be allocated. ( 6) Phase 4 -solution acceptance will probabilistically accept the changes. This is unchanged from the above mentioned previous work, and still uses the EMCAC. The acceptance probability is calculated as follows:

P = exp((o -n)/t) > r ,
where P is the probability of accepting a worsened state, o is the original solution energy cost (modelled as inverse fitness), n is the energy cost of the new solution, t is current temperature and r is a uniformly generated random number between 0 and 1.

Variable Neighbourhood Search

After the ER&SR algorithm has accepted a new solution, VNS is called (if parameters indicate to do so). The VNS method chosen is basic VNS [START_REF] Gendreau | Handbook of metaheuristics[END_REF]. While more complex varieties of VNS are available (such as other hybrid solutions [START_REF] Rahimian | A hybrid integer programming and variable neighbourhood search algorithm to solve nurse rostering problems[END_REF], or with adaptive memory improvements [START_REF] Simeonova | The heterogeneous fleet vehicle routing problem with light loads and overtime: Formulation and population variable neighbourhood search with adaptive memory[END_REF]), ER&SR already provides a feasible close-tooptimal solution. Therefore VNS is used to search nearby neighbourhoods for minor improvements rather than building from the baseline.

Assuming iterations left is above zero, VNS first finds the fitness of the current schedule. Each employee defines a neighbourhood. A random employee is selected, and so is their RP. If HC4 is enabled, only RPs the employee has a preference for are considered for the next step, otherwise all RPs are considered. A random RP from the selected list is chosen. Checking if no HCs are failed by a potential swap, then the swap goes ahead. Then fitness of this new schedule

Order of Events

After the internal VNS has completed, ER&SR repeats until the completion criterion is met. VNS can then run after ER&SR has completed, if parameterised. This would allow the best solution found by ER&SR to be exploited further by VNS. In this paper we apply VNS either a) during ER&SR (shortly after a solution has been accepted) for a number of parameterised iterations, b) after the ER&SR iterations have all completed, or c) both. This means this algorithm can be run with ER&SR only, ER&SR with VNS Internal, ER&SR with VNS Sequential, and ER&SR with both.

RESULTS

In this section we discuss results from tests run on an Intel i5 2.49GHz CPU. The configurations tested were: The baseline is also shown in the results for comparison (the baseline is the fitness computed from the manually generated solution in place at the FSO). The baseline solution is always feasible.

The effectiveness of metaheuristics are largely dependent on well-tuned parameters [START_REF] Burke | Search methodologies[END_REF]. This is particularly important in this instance as there are many parameters and ranges to select from. Thus the parameters for our testing were chosen by the irace parameter configuration tool [START_REF] López-Ibáñez | The irace package: Iterated racing for automatic algorithm configuration[END_REF], such automated parameter tuning being preferable to onerous manual configuration [START_REF] Hoos | Programming by optimization[END_REF]. During irace testing, the testing of elite configurations was enabled, meaning that only the best configurations found during the run of irace, the final best, will be used in the testing phase. The first elimination test value was set to the default of five. The default statistical test configuration for irace was used (F-test).

Each configuration had 30 iterations within irace. The parameters tuned by irace are listed in Table 1; temperature, cooling rate, selection rate, mutation phase 1 rate, mutation phase 2 rate, VNS internal iterations, and VNS sequential iterations. Not all of the parameters were tuned for every algorithm configuration: e.g. cooling rate does not apply to 'VNS only' or the baseline; vnsintr only applies when VNS Internal is enabled. Note that Table 1 shows the parameters for irace for a single test (specifically ER&SR with VNS internally), as an example.

These dependencies are reflected in Table 2, the tuned configurations found by irace. Here, 'N/A' is used when a variable is not tuned. The cooling rate was limited to a maximum of 0.5 due to previous observations during testing finding that too few iterations caused by a higher cooling rate caused poorer results. Similarly the VNS internal number of iterations was capped at 100, as higher values caused a dramatic increase in run-time.

After the parameters were tuned with irace, the algorithms were then each run, with tuned parameters, thirty times. The mean values for these tests can be seen in Table 3. A box-plot showing the range of fitness values can be viewed in Fig. 2. The results shown in this paper reference a single region and around 180 employees being scheduled. HC2 and HC3 have no maximum bound, allowing any number of changes -which both increases the complexity of the problem but also increases the chance of finding better soft constraint compliant allocations. The demand forecast data used for testing is provided by the FSO is calculated by an internal research team for the region the employees tested work.

The SA metaheuristic, given by Config F, has the most variability in terms of fitness provided by the soft constraints. Across all of the configurations that use ER&SR, interestingly there is no solution generated which is worse than the baseline, but that is likely due to the high number of iterations present in the run configuration parameters provided by irace (both directly by VNS iterations and indirectly via temperature & cooling rate values).

A non-parametric ANOVA test (Kruskal-Wallis) was applied to the results and returned a p-value of p < 2.2e -16 . The datasets were confirmed visually as normally distributed in histogram plots as well as in Q-Q plots, which is a plot of the quantiles of the two distributions. The Q-Q plot of Configuration A can be viewed at Figure 3. We can reject the null hypothesis that the results for each configuration are drawn from the same distribution.

There is clear indication that ER&SR provides better schedules, and more so when using the VNS module at any tested stage, than a real world problem provided by a Field Service Operations company solving the employee rostering problem. The novelty of this approach lies in the improvement of the solution compared to standalone ER&SR. A further novelty is application of hyper-parameter tuning to ER&SR (via the irace algorithm configuration tool). This solution also provides an improved schedule for a real world application when compared to the approach currently in place.

/A N/A N/A E N/A N/A N/A N/A N/A N/A 1963 N/A F 9189272 0.004 N/A N/A N/A N/A N/A 1 G N/A N/A N/A N/A N/A N/A N/A N/A
There is also scope for further research on ER&SR, in particular when hybridised with industry solvers such as CPLEX, which would include a mathematical formulation of the problem. There is also potential in replacing the stochastic rebuild section of ER&SR with a different mechanism, for example with an exact method with lowered bounds, other metaheuristics or matheuristics. While in this paper the number of VNS iterations are auto configured by irace, an interesting future development may be a reverse EM-CAC to control the depth of the search of the VNS, so that more exploitative search can take place towards the end of run-time. This may be particularly useful as it was found during experimentation that the ruin method is very destructive early in run-time, negating any benefits of VNS. Additionally, future work should provide random seed trial search dynamics across methods tested to improve replicability of these stochastic algorithms.

DATA ACCESS STATEMENT

The raw employee data nor roster pattern intellectual property cannot be shared due to GDPR and intellectual property laws.

  i t n e s s i s p a r a m e t e r i s e d w e i g h t e d SCs 5 f = g e t F i t n e s s ( ) 6 / / n e i g h b o u r h o o d i s a l l o w e d RPs f o r c o n s i d e r a l l o w e d RPs f o r a . 14 i f ( HC4isTurnedOn ( ) ) 15 s ' = R P F r o m P r e f e r e n c e s ( a ) 16 / / c o n s i d e r a l l RPs .

Fig 1

 1 provides a visual overview of the hybrid algorithm. In the diagram, the internal VNS (to the left) is shorthand for the same process to the right of the diagram. Both the internal VNS and the sequential VNS are only active if parameterised to run.

(

  A) ER&SR with both VNS Internal and VNS Sequential. (B) ER&SR with VNS Internal. (C) ER&SR with VNS Sequential. (D) ER&SR only. (E) VNS only. (F) Simulated Annealing only. (G) Baseline

Figure 1 :

 1 Figure 1: Diagrammatic Overview of ER&SR & VNS Hybrid

Figure 2 :Figure 3 :

 23 Figure 2: Box-Plot of Average Fitness (A-F) and Baseline Fitness (G)

Table 1 :

 1 Parameters To Be Tuned by irace

	Config	Temp	Cooling Rate	Selection Rate	Mutation P1 Rate	Mutation P2 Rate	VNS Internal	VNS Sequential	Simulated Annealing
	A	116621.6	0.045	0.434	0.041	0.353	35.951	743.145	N/A
	B	33904.68	0.221	0.707	0.045	0.597	81.053	N/A	N/A
	C	711429.8	0.024	0.301	0.791	0.486	N/A	103	N/A
	D	706984.6	0.007	0.881	0.721	0.696	N

Table 2 :

 2 Tuned Parameters the schedules created by SA, or by VNS alone. The results show that the method this paper provides is an improved alternative to current RPs / starting weeks employees are working on the baseline solution for this dataset and demand forecast. The schedules are improved in terms of average soft constraint fitness, which is how the quality of a solution is judged. View Table4for comparison of demand hours met, which shows unmet demand is reduced by over 28% from Config A compared to the baseline of Config G.The results also show that this algorithm can provide a small but relatively quick improvement upon standard ER&SR, where Config ER&SR VNS Internal VNS Sequential Simulated Annealing Average Fitness

	A	1	1	1	0	0.55993
	B	1	1	0	0	0.54985
	C	1	0	1	0	0.54227
	D	1	0	0	0	0.51884
	E	0	0	1	0	0.49834
	F	0	0	0	1	0.44840
	G	0	0	0	0	0.42987

Table 3 :

 3 Average Fitness Per Run Type

Table 4 :

 4 Demand Met Per Config

	Config	Original Demand	Demand Hours Met	Demand Hours Not Met
	A	80305	55426	24879
	B	80305	55198	25107
	C	80305	55180	25125
	D	80305	54699	25606
	E	80305	52980	27325
	F	80305	51512	28793
	G	80305	45373	34932

6 CONCLUSIONS

This paper presents a hybrid algorithm combining Evolutionary Ruin & Stochastic Recreate (ER&SR) with Variable Neighbourhood Search (VNS) and compares it to a Simulated Annealing (SA) approach and a first-improvement VNS approach. We apply this to

ACKNOWLEDGEMENTS

The authors of this paper extend gratitude to those involved and the following funders: British Telecommunications Plc and the UK's EPSRC DAASE project (grant no. EP/J017515/1). J. Swan acknowledges the support of the EU H2020 SAFIRE project (Ref. 723634). Additional thanks to Dr Paul McMenemy for additional support.