Clément Legrand
email: clement.legrand4.etu@univ-lille.fr

Diego Cattaruzza
email: diego.cattaruzza@centralelille.fr

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

Marie-Eléonore Kessaci
email: marie-eleonore.kessaci@univ-lille.fr

Improving Neighborhood Exploration into MOEA/D Framework to Solve a Bi-Objective Routing Problem

Local search (LS) algorithms are efficient metaheuristics to solve combinatorial problems. The performance of LS highly depends on the neighborhood exploration of solutions. Many methods have been developed over the years to improve the efficiency of LS on different problems of operations research. In particular, the exploration strategy of the neighborhood and the exclusion of irrelevant neighboring solutions are design mechanisms that have to be carefully considered when tackling NP-hard optimization problems. A MOEA/D framework including an LS-based mutation and knowledge discovery mechanisms is the core algorithm used to solve a bi-objective vehicle routing problem with time windows (bVRPTW) where the total traveling cost and the total waiting time of drivers have to be minimized. We enhance the classical LS exploration strategy of the neighborhood from the literature of scheduling and propose new metrics based on customers distances and waiting times to reduce the neighborhood size. We conduct a deep analysis of the parameters to give a fine-tuning of the MOEA/D framework adapted to the LS variants and to the bVRPTW. Experiments show that the proposed neighborhood strategies lead to better performance on both Solomon's and Gehring and Homberger's benchmarks.

Introduction

With the increasing demand in terms of services (transportation, delivery, storage), comes an explosion of new challenges that have to be faced. Most of them require taking into account several conflicting aspects (ecological, economical, societal), each one characterized by an objective function to optimize. When at least two objectives are optimized simultaneously, we are facing a multi-criterion decision. In some cases, the decision maker provides a score of interest for each objective, which can lead to a reduction of the multi-objective problem to a single-objective problem where the different objectives are aggregated. Unfortunately, most of the time, the decision maker does not know precisely how the objectives interact with each other, which may lead to a bad strategy to solve the problem, and consequently to irrelevant results.

In the context of this paper, we are interested in the Vehicle Routing Problem with Time Windows (VRPTW). It can be considered a logistic problem, where vehicles have to serve customers within a precise time interval. This problem has been introduced decades ago, and many objectives have been investigated. 1 However, to the best of our knowledge, only a small subset of these objectives has received enough interest to have publicly available results. Thus, it is very hard to know if an algorithm or a method is competitive on the problem itself (independent of the objectives chosen) or if it is only relevant to the objectives studied. To overcome this lack in the literature, we propose to focus on a bi-objective VRPTW (bVRPTW) where one objective has been highly studied, that is the total transportation cost, and an objective that has rarely been studied, that is, the total waiting time incurred when drivers arrive before the opening of the time window [START_REF] Castro-Gutierrez | Nature of real-world multi-objective vehicle routing with evolutionary algorithms[END_REF][START_REF] Zhou | A local search-based multiobjective optimization algorithm for multiobjective vehicle routing problem with time windows[END_REF]. Using both of these objectives is interesting in many real-life situations, like food delivery, where the waiting time of a driver impacts the heat of the meals of the next customers, and consequently customer satisfaction. Another typical situation concerns the transportation of people, more precisely when a patient has a medical appointment, we do not want that he/she waits too much.

To solve the problem, we use MOEA/D, a Multi-Objective Evolutionary Algorithm based on Decomposition [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF] where the mutation step is replaced by a local search, and enhanced with Knowledge Discovery (KD) [START_REF] Legrand | Enhancing MOEA/D with learning: Application to routing problems with time windows[END_REF]. Local searches (LS) are known to be powerful algorithms used in evolutionary algorithms to improve their performance [START_REF] Knowles | Local-search and hybrid evolutionary algorithms for Pareto optimization[END_REF]. Indeed, LS are able to intensify the search by focusing on a specific region of the search space. They are based on neighborhood operators that link solutions together and a neighborhood exploration strategy defines how the neighbors are explored and when the exploration is stopped. Many LS have been developed, and many operators are available, for this problem, but most of them only consider a single objective which is the total transportation cost. Moreover, LS steps are quite time-consuming, which is why different strategies exist to speed up the search and reduce the time allocated to neighborhood exploration. Indeed, we can explore more or less the neighborhood according to the strategy used. With a first strategy, the first improving solution met is kept, allowing a fast exploration of the space but with a relatively long path to a local optimum. While a best strategy explores the whole neighborhood to apply the best improving move, allowing a smaller path but the exploration itself is more time-consuming. Since routing problems produce large neighborhoods, reduction techniques have been designed to avoid irrelevant moves. The most common one is probably the granular search introduced by [START_REF] Toth | The granular tabu search and its application to the vehicle-routing problem[END_REF]. It is based on the idea that two distant customers have a low chance to be served one right after the other in the same route. However, it requires a proper definition of the distance between customers, which may not be obvious when several objectives are considered and are not in the same range of values.

The contribution of the paper is to enhance the basic local search used in MOEA/D considering new neighborhood mechanisms. First, we present a new strategy to explore the neighborhood of bVRPTW solutions inspired by the state of the art for permutation flow-shop scheduling problem [START_REF] Ruiz | A simple and effective iterated greedy algorithm for the permutation flowshop scheduling problem[END_REF]. The neighborhood is then divided into subsets of neighbors. Contrary to the best strategy that explores the neighborhood entirely before making a decision or the first strategy that makes a decision for each visited neighbor, the proposed strategy, called first-best, explores a subset entirely before making a decision. Second, we define a new metric that considers not only the Euclidean distance between the customers but also their respective time windows to evaluate the waiting time. With such a metric, we hope to obtain a more accurate neighborhood, depending on the two objectives optimized. Two variants of the metric are proposed. One variant aggregates the cost and the waiting time with variable weights, while the other uses the same weights for the aggregation. These contributions are supported by a rigorous analysis of the different parameters used in the MOEA/D framework, through irace [START_REF] López-Ibáñez | The irace package: Iterated racing for automatic algorithm configuration[END_REF], to highlight the most influential ones. Therefore, each variant tested is then tuned with irace to keep the configuration the best adapted for each one. Experiments are conducted on Solomon's benchmark with instances of size 100, and on Gehring and Homberger's benchmark with instances of size 200 to evaluate the performance of each MOEA/D variant and define the best-suited neighborhood mechanisms.

The paper is structured as follows. Section 2 introduces the problem and motivates our choice of ob-jectives to optimize. Concepts related to multi-objective optimization are presented in Section 3. Our contribution concerning the local search, as well as the algorithms compared, are detailed in Section 4. The experimental setup (i.e., benchmarks, termination criterion, performance assessment, and machines used) is described in Section 5. Our experimental protocol, concerning the tuning and the evaluation of the performances, is presented in Section 6. The tuning is analyzed in Section 7, and the results obtained are discussed in Section 8. We conclude the paper in Section 9.

Description of the Problem

Vehicle Routing Problem with Time Windows (VRPTW)

The VRPTW [START_REF] Toth | Vehicle routing: problems, methods, and applications[END_REF] is defined on a graph G = (V, E), where V = {0, 1, . . . , N } is the set of vertices and E = {(i, j) | i, j ∈ V } is the set of arcs. It is possible to travel from i to j, incurring a travel cost c ij and a travel time t ij . Usually, the cost is computed as the Euclidean distance between the customers. Vertex 0 represents the depot where a fleet of K identical vehicles with limited capacity Q is based. Vertices 1, . . . , N represent the customers to be served, each one having a demand q i , a time window [a i , b i ] during which service must occur, and a service time s i estimating the required time to perform the delivery. Vehicles may arrive before a i . In that case, the driver has to wait until a i to accomplish service incurring a waiting time. Arriving later than b i is not allowed (the time windows are hard). We recall that a route r is an elementary cycle on G that contains the depot (that is vertex 0) and is expressed as a sequence of vertices r = (v 0 , v 1 , . . . , v |r| , v |r|+1 ) where v 0 = v |r|+1 = 0 and vertices v 1 , . . . , v |r| are all different. The cost c r of a route r is then given as the sum of traveling costs on arcs used to visit subsequent vertices, that is

|r| i=0 c vi,vi+1
. A solution x is represented as a set of K (possibly empty) routes, that is x = {r 1 , . . . , r K }, and its cost is expressed as:

f 1 (x) = K k=1 c rk (1)
The VRPTW calls for the determination of at most K routes such that the traveling cost is minimized and the following conditions are satisfied: (a) each route starts and ends at the depot, (b) each customer is visited by exactly one route, (c) the sum of the demands of the customers in any route does not exceed Q, (d) time windows are respected.

Objectives Related to the VRPTW

A wide range of objectives has already been investigated in the VRPTW context. This section provides a broad overview of the most common objectives found in the literature. Following the classification of [START_REF] Jozefowiez | Multi-objective vehicle routing problems[END_REF], these objectives are divided into three categories according to the component of the problem they are associated with the route, the node or arc activity, and the resources.

Objectives Associated to the Route

Among the objectives related to the route, the most common is surely minimizing the cost of the solutions generated [START_REF] Schneider | Designing granular solution methods for routing problems with time windows[END_REF]. For instance, the cost can be the distance traveled or the time required. It can also be something more complex, like CO 2 emissions. More generally, minimizing cost is linked to an economic criterion. When the total cost is ignored, the makespan (i.e., minimizing the length of the longest route) is an objective frequently minimized [START_REF] Zhou | A local search-based multiobjective optimization algorithm for multiobjective vehicle routing problem with time windows[END_REF]. These objectives are motivated by the applications of the problem. For example, if electric vehicles are used, the length of the tours should not exceed the autonomy of the vehicles.

Objectives Associated to Nodes or Arcs

Concerning the objectives related to the node or arc activity, they often involve time windows, since it is the most constraining part of the problem in general. When considering hard time windows (i.e., that cannot be violated), the driver's waiting time due to earliness can be optimized (Zhang et al., 2019). If we consider soft time windows (i.e., that can be violated with a penalty), the delay time of drivers is more likely to be optimized [START_REF] Zhou | A local search-based multiobjective optimization algorithm for multiobjective vehicle routing problem with time windows[END_REF]. More generally, the number of violated time windows can also be minimized [START_REF] Geiger | A computational study of genetic crossover operators for multi-objective vehicle routing problem with soft time windows[END_REF].

Objectives Associated to Resources

The last category of objectives concerns resources. The main resources encountered in the literature are vehicles and goods. The minimization of the number of vehicles often appears in classical benchmarks and can be interpreted economically, in that fewer vehicles means less monetary investment. More precisely, this objective is often minimized first in Solomon's and Gehring and Homberger's benchmarks [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF][START_REF] Gehring | A parallel hybrid evolutionary metaheuristic for the vehicle routing problem with time windows[END_REF], since buying a truck or hiring a driver is the most expensive part. Finally, some objectives try to erase disparities between tours, to bring a fairness aspect into the problem. To define a balancing objective it is necessary to define a route's workload, which can be expressed as the number of customers visited or the number of goods delivered [START_REF] Melián-Batista | A bi-objective vehicle routing problem with time windows: A real case in tenerife[END_REF][START_REF] Baños | A hybrid meta-heuristic for multi-objective vehicle routing problems with time windows[END_REF].

Minimizing the Waiting Time of Drivers

When drivers arrive before the opening of a time window they must wait until the opening time. It increases the time of the route for the driver and may incur satisfaction issues. The total waiting time of drivers is an objective that has not been highly studied in the literature. In the last decade, [START_REF] Castro-Gutierrez | Nature of real-world multi-objective vehicle routing with evolutionary algorithms[END_REF] studied how this objective was correlated to four other common objectives in the VRPTW context. In particular, they showed that in Solomon's instances, the total traveled distance and the total waiting time were weakly correlated in clustered instances with tight time windows and in random instances with wide time windows. In addition, the objectives are in harmony (i.e., the minimization of one tends to also decrease the other) on clustered instances with wide time windows, whereas they are conflicting (i.e., the minimization of one tends to increase the other) on random instances with tight time windows. The total waiting time due to early arrival has also been studied by [START_REF] Zhou | A local search-based multiobjective optimization algorithm for multiobjective vehicle routing problem with time windows[END_REF], where a five-objective VRPTW is solved with an objective-wise local search. Moreover, it seems interesting to study together two continuous objectives (here the total traveled distance and the total waiting time) instead of using a discrete objective, like the number of vehicles. Indeed, with two continuous objectives, the fronts generated by the algorithm contain, generally, much more non-dominated solutions, that are relevant when using knowledge discovery methods.

We formalize the notion of waiting time as follows. The waiting time W i at a customer i is given as the maximum between 0 and the difference between the opening of the TW a i and the arrival time T i at location i, that is W i = max{0, a i -T i }. Note that each route r = (v 0 , v 1 , . . . , v |r| , v |r|+1 ) is associated with a feasible (i.e., consistent with traveling times and TWs) arrival time vector T r = (T v0 , T v1 , . . . , T v|r| , T v|r|+1 ) and the total waiting time W r (T r ) on route r, with respect to T r is given by W r (T r ) = |r| i=1 W vi . Thus the total waiting time of a solution x = {r 1 , . . . , r K } on a graph G, given a time arrival vector for each route in the solution, i.e.,T x = (T r1 , . . . , T rK ), is given by the following formula:

f 2 (x, T x ) = K k=1 W rk (T rk ) (2)
In the remainder of the paper, the problem considered is a Bi-Objective VRPTW (bVRPTW), where functions f 1 (Equation 1) and f 2 (Equation 2) are simultaneously minimized, hence forming a multi-objective problem. Key concepts to understand how to solve a multi-objective problem are given in Section 3.

MOEA/D Framework

Multi-Objective Optimization

A Multi-objective Combinatorial Optimization Problem (MoCOP) is formalized as follows [START_REF] Coello | Multi-objective combinatorial optimization: Problematic and context[END_REF]:

(M oCOP ) = Optimize F (x) = (f 1 (x), f 2 (x), . . . , f n (x)) s.t. x ∈ D, (3)
where n is the number of objectives (n ≥ 2), x is the vector of decision variables, D is the (discrete) set of feasible solutions and each objective function f i (x) has to be optimized (i.e.,minimized or maximized). In multi-objective optimization, the objective function F defines a so-called objective space denoted by Z. For each solution x ∈ D there exists a point in Z defined by F (x).

A dominance criterion is defined to compare solutions together: a solution x dominates a solution y, in a minimization context, if and only if for all i ∈ [1 . . . n], f i (x) ≤ f i (y) and there exists j ∈ [1 . . . n] such that f j (x) < f j (y). The dominance relation introduces a partial order on the solution: there are indeed solutions that cannot be compared among each other.

Then a set of non-dominated solutions is called a Pareto front. A feasible solution x * ∈ D is called Pareto optimal if and only if there is no solution x ∈ D such that x dominates x * . Resolving a MoCOP involves finding all the Pareto optimal solutions which form the Pareto optimal set. The true Pareto front of the problem is the image of the Pareto optimal set by the objective function.

Over the years, many metaheuristics based on local search techniques or using evolutionary algorithms [START_REF] Blot | Survey and unification of local search techniques in metaheuristics for multi-objective combinatorial optimisation[END_REF] have been designed to solve multi-objective problems. Moreover, many tools [START_REF] Riquelme | Performance metrics in multi-objective optimization[END_REF] have been developed to assess and compare the performance of multi-objective algorithms. When the results are compared to a Pareto front of reference, it is common to use the generational distance (GD) metric or the epsilon metric. The GD measures the cumulative (Euclidean) distance from each solution of the front generated to the closest solution in the reference front. The epsilon metric finds the smallest value ϵ, so that, if we multiply the objective vector of all the solutions obtained by ϵ, then the new front strictly dominates the front of reference. However, we do not have access to reference fronts for the bVRPTW studied, which makes it difficult to use these metrics. That is why, in this paper, we use the unary hypervolume (uHV) [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF]. It is a metric defined relatively to a reference point Z ref .

This indicator evaluates the accuracy, diversity, and cardinality of the front, and it is the only indicator with this capability. Moreover, it can be used without knowing the true Pareto front of the problem. It reflects the volume covered by the members of a non-dominated set of solutions. Thus, the larger the hypervolume, the better the set of solutions.

Reference MOEA/D

MOEA/D [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF] is a genetic algorithm that approximates the true Pareto front by decomposing the multi-objective problem into M scalar objective subproblems. The scalarization is obtained by weighting each of the n objectives f k with a weight w k ∈ [0, 1], such that n k=1 w k = 1. Thus the fitness of a solution x for the subproblem i, defined with a weight w i = (w i 1 , w i 2 , . . . , w i n ), is:

f (x|w i ) = n k=1 w i k • f k (x) (4)
During an iteration, MOEA/D generates a new solution for the i-th subproblem by breeding solutions obtained from close (notion defined hereafter) subproblems. The neighborhood, of size m, of a weight vector w i is defined as the set of its m closest (by means of the Euclidean distance) weight vectors among {w 1 , . . . , w M }. Then the neighborhood N m (i) of the i-th subproblem consists of the m subproblems defined with a weight vector belonging to the neighborhood of w i . Note that each subproblem is associated with its best solution found during the execution.

At the start of MOEA/D, M uniformly distributed weight vectors are determined and then it proceeds as described in Algorithm 1. For each subproblem, a random initial solution is generated and evaluated, and all these solutions form together the population used in the algorithm. Note that, the size of the population does not change during the execution of the algorithm, and each solution of the population is associated with exactly one subproblem. Each subproblem's neighborhood (of size m) is also computed. When optimizing subproblem i, a random pair of solutions is selected from its neighborhood. The Partially Mapped crossover (PMX), developed by [START_REF] Goldberg | Alleles, loci, and the traveling salesman problem[END_REF], is applied with probability p pmx , and only one solution is randomly kept. Then, a Local Search (LS), described in Section 4, is applied with probability p ls . Indeed, the mutation is frequently replaced by an LS [START_REF] Knowles | Local-search and hybrid evolutionary algorithms for Pareto optimization[END_REF] in genetic algorithms. Finally, the resulting solution is added to the set S of solutions generated during the iteration, and a few neighbors (at most two in our case) of the subproblem i are updated. When all subproblems have been considered, S is merged with the archive A, containing the non-dominated solutions generated during the execution of the algorithm. It is important to use an external archive since the population may not be able to track all the non-dominated solutions generated due to the discretization of the objective space. If memory issues are encountered, it is common to bound the size of A, but here, with the bVRPTW, it was not necessary. If the termination criterion is reached, the solutions of A are returned, otherwise, a new iteration is started. In the remaining of the paper, this algorithm is our reference algorithm, called R M OEA/D .

Hybrid MOEA/D

Description of the Hybrid Framework

Hybridizing machine learning methods and metaheuristics has become quite common to solve combinatorial problems. In [START_REF] Legrand | Improving moea/d with knowledge discovery. application to a bi-objective routing problem[END_REF], MOEA/D is hybridized with a knowledge discovery framework, where the notion of knowledge groups is described. The fitness space is divided into k G regions each representing a knowledge group. Each group is characterized by a vector, the reader is referred to the original paper for further details. The region of each group is defined using the subproblems of MOEA/D, in particular, if a solution belongs to the region of a knowledge group, then the knowledge extracted from the solution is added to that group. For this paper, we consider that each subproblem is associated with one group, being the closest one, considering the Euclidean distance between the aggregation of the subproblem and the Algorithm 1: Reference MOEA/D Framework (R M OEA/D ).

Input: M weight vectors w 1 , . . . , w M . Output: The external archive A /* Initialisation */ 1 (A, S) ← (∅, ∅) 2 P ← random initial population (x i for the i-th subproblem)

3 for i ∈ {1, . . . , M } do 4 N (i) ← indexes of the m closest weight vectors to w i 5 Obj i ← {f k (x i ) | 1 ≤ k ≤ n} /* Core of the algorithm */ 6 while stopping criterion not satisfied do 7 for i ∈ {1, . . . , M } do 8 (i 1 , i 2 ) ← Select (N (i)) 9 x ← PMX(x i1 , x i2 ) x ← LS(x) S ← S ∪ {x} updateNeighbors(P, N (i), x) A ← updateArchive(A, S) S ← ∅ return A
vector of each group. Moreover, the framework relies on two mechanisms, being the knowledge extraction (K ext ) and the knowledge injection (K inj ). Here, both mechanisms follow an intensification strategy, meaning that for the extraction (resp. injection) the knowledge extracted (resp. injected) belongs to one group only. In the context of the bVRPTW, K ext and K inj are described in Section 3.3.2.

The framework of the hybridization is shown in Algorithm 2. Note that, if lines 3, 11, and 16 are removed, then the algorithm becomes the original MOEA/D described in Section 3.2. At line 3, the procedure createGroups creates all the knowledge groups. At line 11, the injection procedure K inj is applied to the current solution x, with probability p inj , using the knowledge stored in the closest group (intensification strategy). At line 16, the extraction procedure K ext is used to extract the knowledge from the set of solutions generated during the iteration. Then it updates the closest group (intensification strategy) of the subproblem being optimized.

Knowledge Mechanisms for the bVRPTW

In this section, the K inj and K ext learning mechanisms are defined to suit the bVRPTW. The mechanisms are inspired by [START_REF] Arnold | PILS: Exploring high-order neighborhoods by pattern mining and injection[END_REF] who developed PILS, an optimization method relying on frequent patterns from high-quality solutions to explore vast neighborhoods. PILS has already been integrated into the Hybrid Genetic Search (HGS) of [START_REF] Vidal | Multiobjective memetic algorithm based on adaptive local search chains for vehicle routing problem with time windows[END_REF] and the Guided Local Search (GLS) of [START_REF] Arnold | Knowledge-guided local search for the vehicle routing problem[END_REF] to solve the Capacitated Vehicle Routing Problem (CVRP) with one single objective and provided good results. Usually, patterns are defined as sequences of consecutive customers. Those with a size between 2 and size p , a user-defined parameter, are extracted from generated solutions by K ext . Since the depot belongs to all routes, it is not considered inside patterns. In particular, a route r = (0, v 1 , . . . , v |r| , 0), contains max(|r|-k +1, 0) patterns of size k. The extracted patterns are added to the corresponding group, where their frequency is incremented. K inj injects a set of N Inj patterns in the solution x provided. Note Algorithm 2: Hybrid MOEA/D Framework (H M OEA/D ).

Input: M weight vectors w 1 , . . . , w M . The number k G of knowledge groups.

Output: The external archive

A /* Initialisation */ 1 (A, S) ← (∅, ∅) 2 P ← random initial population (x i for the i-th subproblem) 3 G ← createGroups(k G ) 4 for i ∈ {1, . . . , M } do 5 N (i) ← indexes of the m closest weight vectors to w i 6 Obj i ← {f j (x i ) | 1 ≤ j ≤ n} /* Core of the algorithm */ 7 while stopping criterion not satisfied do 8 for i ∈ {1, . . . , M } do 9 (i 1 , i 2 ) ← Select (N (i)) x ← PMX(x i1 , x i2 ) x ← K inj (x, G, i) x ← LS(x) S ← S ∪ {x} updateNeighbors(P, N (i), x) A ← updateArchive(A, S) G ← K ext (G, S) S ← ∅ return A
that, a pattern is kept in the solution only when the solution is improved. To select a pattern, first, its size is randomly chosen among {2, . . . , size p }. Hence the selection is not biased towards smaller, more numerous, patterns. Then the pattern is randomly chosen among the N F requent most frequent patterns of the same size, which is also a parameter of the algorithm. When all the N Inj patterns have been selected, they are injected one by one according to the following steps. To inject a pattern, arcs connecting it are removed and the nodes are reconnected to form a feasible solution. The reconnection is done optimally, but reversed patterns are not considered due to time windows.

Neighborhood Exploration Mechanisms for Local Search

A solution to the problem is represented as a permutation of the customers, and it is evaluated with the split algorithm provided by [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF]. For this study, we consider the same neighborhood as defined by [START_REF] Schneider | Designing granular solution methods for routing problems with time windows[END_REF]. The applied operators are: Relocate, Swap, and 2-opt * . These operators are commonly used in local search algorithms for routing problems since they are basic operators and they are able to produce a large neighborhood. The relocate operator moves one customer to another location (possibly on the same route). The swap operator exchanges two customers. The 2-opt * operator generalizes the 2-opt (that is an exchange of two arcs in the same route), by involving different routes. The construction of the neighborhood is similar to the one used in the Randomized Variable Neighborhood Descent (see, e.g., [START_REF] Subramanian | A hybrid algorithm for a class of vehicle routing problems[END_REF], where the order of the neighborhood operators is kept during descent but shuffled each time the LS is re-applied.

The exploration strategy of the neighborhood (of an operator) is discussed in Section 4.1. Concerning the granular search, used to reduce the neighborhood, different metrics are proposed in Section 4.2 to evaluate how close customers are.

Exploration of the Neighborhood

In this section, more details are given on how the neighborhood of a solution is explored. In local search, the two classical neighborhood exploration strategies are the best and the first improvement methods [START_REF] Hoos | Stochastic local search: Foundations and applications[END_REF]. Both of them start with the current solution and try to find an improving solution in their neighborhood that is defined by an operator. While the best improvement generates all the neighboring solutions of the current solution to find the neighbor with the best improvement, the first improvement evaluates the neighboring solutions one by one and stops as soon as it finds one improving neighbor. Obviously, the best strategy guarantees the best possible progress at each iteration of the local search but the time allotted to this increases with the size of the problem. On the other hand, the first strategy allows a fast exploration but many improving neighbors are set aside. In routing problems, the most commonly used neighborhood exploration strategy is the classical best strategy [START_REF] Subramanian | A hybrid algorithm for a class of vehicle routing problems[END_REF][START_REF] Schneider | Designing granular solution methods for routing problems with time windows[END_REF][START_REF] Accorsi | A fast and scalable heuristic for the solution of large-scale capacitated vehicle routing problems[END_REF], where the best neighboring solution found by the application of one operator on the current solution is selected. The operators presented above (Relocate, Swap, and 2-opt * ) provide a number of possible neighbors that depends on the number of customers. Therefore, the neighborhood of a solution increases with the size of the problem. When considering large problems, it may be necessary to speed up each neighborhood exploration. However, the first strategy does not consider the best move for a considered customer. A compromise between these two classical strategies called first-best has been proposed in the literature for scheduling problems, particularly for permutation flow-shop [START_REF] Ruiz | A simple and effective iterated greedy algorithm for the permutation flowshop scheduling problem[END_REF]. Algorithm 3 gives the pseudo-code of the first-best procedure. The procedure requires a neighborhood operator (e.g., Swap, Relocate, or 2-opt * ), and the solution x which undergoes the LS. For the given operator each customer (l.5) is tentatively moved to its best location until an improving move is found (l.8). When it is the case, the customer is moved, and a new iteration starts. The authorized moves are generated by the procedure generate moves (l.6). The algorithm stops when the whole neighborhood is explored without finding any improving move.

Neighborhood Metric

Another way to accelerate neighborhood exploration is the reduction of the neighborhood size. In routing problems, it is known that most of the arcs in good solutions involve close customers. Therefore, it may be interesting to consider the closeness between customers when we generate the neighboring solutions of a solution. The problem is to define this closeness. Indeed, for the classical VRPTW, when minimizing the total traveled distance, the Euclidean distance between customers seems to be enough to quantify the closeness between customers. However, when two (or more) objectives are considered, it is not as easy to define a relevant metric. Once a metric is defined, a natural way to reduce the neighborhood of an operator is to consider moves including the δ nearest customers for the metric defined. Roughly speaking, this is the concept behind the granular search that was introduced by [START_REF] Toth | The granular tabu search and its application to the vehicle-routing problem[END_REF]. Having a method that restricts the neighborhood to moves that seem relevant, is interesting to spare time and resources during the LS.

For our study, we compare three different metrics. The first metric, called d 1 , is the classical metric x ′ ← best solution obtained by applying a move from moves

8 if f (x ′ ) < f (x) then 9 x ← x ′ 10 improve ← T rue 11 break 12 return
x used in single-objective routing problems: the closeness between two customers is simply evaluated as the Euclidean distance between them. The second metric, d 2 , is an aggregation between two costs, each one related to one of the objectives. The first cost represents the Euclidean distance between the customers (d 1 in fact), and it is related to the total traveled distance objective. The second cost evaluates the waiting time incurred by going from one customer to another one (not symmetric in general due to the TW). This cost is related to the total waiting time of drivers. More precisely, each subproblem generated in MOEA/D, with weight vector w = (w 1 , w 2 ), has its own metric that is defined as:

d w 2 (u, v) = w 1 • distance(u, v) -M inDistance M axDistance -M inDistance + w 2 • W T (u, v) -M inW T M axW T -M inW T (5)
where u and v are two customers, M axDistance (resp. M inDistance) is the largest (resp. smallest) distance between two customers, and M axW T (resp. M inW T ) is the largest (resp. smallest) value of WT between two customers (as defined by Equation 6). These values are used to normalize the distances and waiting times, and they are computed at the beginning of the execution since they only depend on the characteristics of the instance. The value W T (u, v) is the waiting time incurred by going to v from u.

If [a u , b u ] (resp. [a v , b v ])
is the time window of customer u (resp. v), s u the service time of customer u and t uv the traveling time from u to v, then W T is expressed as follows:

W T (u, v) = max(0, a v -(a u + s u + t uv )) (6)
The third and final metric, d 3 is a particular case of d 2 where the weights of the two costs are the same (i.e., 0.5 for both). This metric is interesting because it is independent of the aggregations defined in MOEA/D, and it can therefore be used in another framework, without defining any aggregations.

The two strategies presented in Section 4.1, and the three metrics defined in this section, lead to six variants, being f bd 1 , f bd 2 , f bd 3 , bd 1 , bd 2 and bd 3 , where f b (resp. b) stands for first-best (resp. best) strategy and d X is the metric. These six variants are generalized through the H M OEA/D framework (see Algorithm 2), where the exploration strategy and the metric are parameters of the LS. These variants are also compared with the reference algorithm, R M OEA/D , described in Section 3.2.

Experimental Setup

Benchmarks

Solomon's benchmark [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF] is a set of VRPTW instances, frequently used in the literature to evaluate the performance of multi-objective algorithms [START_REF] Ghoseiri | Multi-objective vehicle routing problem with time windows using goal programming and genetic algorithm[END_REF][START_REF] Qi | A decomposition based memetic algorithm for multi-objective vehicle routing problem with time windows[END_REF][START_REF] Moradi | The new optimization algorithm for the vehicle routing problem with time windows using multi-objective discrete learnable evolution model[END_REF]. Thus, we use this benchmark to evaluate the performance of all seven variants. The benchmark contains 56 instances, each with 100 customers, divided into three categories according to the type of generation used: R (random), C (clustered), or RC (random-clustered). The R category (23 instances) contains instances where customers are randomly located in a 100 × 100 grid, while the instances of category C (17 instances) contain clusters of customers. The category RC (16 instances) contains instances where half the customers are randomly located and the other half is clustered. Each category is divided into two classes, either 1 or 2, according to the width of time windows. Instances of class 1 have tighter time windows than instances of class 2, meaning that instances 1 are more constrained. For our experiments, the RC instances are discarded, to highlight the differences in the tuning between random instances and clustered instances.

We also consider bigger instances from the benchmark of Gehring and Homberger [START_REF] Gehring | A parallel hybrid evolutionary metaheuristic for the vehicle routing problem with time windows[END_REF]. This set extends the Solomon one and considers a larger number of customers. Indeed, it contains instances of size 200, 400, 600, 800, and 1000. These instances are very similar to the instances of Solomon since there are the three categories R, C, and RC, and for each category, two classes (1 or 2), depending on the width of the time windows. Contrarily to Solomon's set, each category of instance is equally represented (each one contains 10 instances with tight time windows (class 1) and 10 instances with wide time windows (class 2)) For our experiments, we considered the instances of size 200, without the RC category.

Termination Criterion and Performance Assessment

Given the size N of an instance (i.e., the number of customers), we set the termination criterion of all the variants to 9 * N 2 /125 seconds. It means that the algorithms are executed during 720 seconds (i.e., 12 minutes) on instances of size 100 and 2880 seconds (i.e., 48 minutes) on instances of size 200. During our preliminary studies, it was enough to obtain accurate and robust results.

The performance of the solution sets is evaluated using the unary hypervolume [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF]) (uHV), which quantifies the volume of the region dominated by the solutions in the set. Since the true Pareto front of the problem is unknown, metrics that depend on it (such as epsilon or GD) cannot be utilized. The normalization points for the objectives of each solution are obtained through experiments and are updated dynamically when a new point is discovered. Known optimal values concerning the total traveled distance are available at CVRPLIB. To compute the uHV we normalize the values of the objectives with the normalization points obtained and then we use the point (1.001, 1.001) as a reference. The experiments are run on two computers "Intel(R) Xeon(R) CPU E5-2687W v4 @ 3.00GHz", with 24 cores each. The variants have been implemented using the jMetalPy framework [START_REF] Benitez-Hidalgo | jmetalpy: A python framework for multi-objective optimization with metaheuristics[END_REF]. The code is available at https://github.com/Clegrandlixon/data itor2023.

Experimental Protocol

Tuning

The reference MOEA/D (R M OEA/D ) uses the following parameters: M , the number of subproblems considered, and m the size of the neighborhood of each subproblem. The probabilities associated with each mechanism are p pmx for the crossover and p ls for the LS. The granularity parameter δ is used to reduce the neighborhood during LS (cf. Section 4.2). The exploration strategy S (either best or first-best) and the metric d (either d 1 , d 2 , or d 3 ) are also parameters of MOEA/D, and the choice is left to the tuning. The hybrid MOEA/D (H M OEA/D ) has, in addition, the following parameters: the number of knowledge groups created k G , the probability of injection p inj , the maximal size size p of the patterns extracted, and the number N Inj of patterns injected, chosen among the N F requent most frequent patterns. According to a preliminary study [START_REF] Legrand | Improving moea/d with knowledge discovery. application to a bi-objective routing problem[END_REF] and existing works [START_REF] Arnold | PILS: Exploring high-order neighborhoods by pattern mining and injection[END_REF], we set m = 1/4 × M and N F requent = 5 × N Inj . We propose a different range of values for the eight remaining parameters (cf. Table 1) to define the configuration space in irace.

The tuning is decomposed into three main steps. A first tuning is performed on all the variants, with the whole configuration space defined above. A budget of 2000 is granted to irace, spread between 8 iterations. We generated a set of 20 instances (of size 100) per category to perform the tuning, following the guidelines provided by [START_REF] Uchoa | New benchmark instances for the capacitated vehicle routing problem[END_REF]. The generated instances mimic Solomon's instances. Moreover, we separate the tuning of the algorithms on instances R and C, since they do not have the same structure. The tuning of the original MOEA/D allows seeing which neighborhood exploration strategy and metric are better when not performing learning steps. In addition, the tuning of the hybrid framework grants a set of good configurations that can be used to analyze the influence of parameters. Moreover, we perform an additional tuning, where the exploration strategy and the metric are not set in case of learning (i.e., tuning of H M OEA/D directly), allowing us to compare which strategy and metric are better when the learning mechanism is used. See Section 7.1 for more details.

Secondly, we analyze the elite configurations returned for the hybrid variants. Each configuration returned is a local optimum for irace and thus should be the best configuration (statistically) in its close neighborhood. For each configuration, we generate all the possible configurations distant by 1 (with the hamming distance). More precisely, only one parameter is changed at a time (the possible values are the values available in Table 1), while keeping the values of the other parameters unchanged. In order to check whether the configurations returned by irace are good or not, we evaluate the generated configurations on the tuning instances (like before), with the same seed to be fair with the methodology of irace. The hypervolumes obtained are compared in Section 7.2.

In the third part, we use the results of the experiment above to adapt and reduce the configuration space. With the reduction of the space, irace should return more accurate results. Moreover, the number of aggregations (i.e., the number of subproblems) generated is set to 40. By setting this parameter, it is easier to see the impact of the number of groups on the performances of the algorithms. We perform the new tuning of the six variants f bd 1 , f bd 2 , f bd 3 , bd 1 , bd 2 , and bd 3 on the new space. The budget allocated to irace is reduced accordingly to 1000. The instances are those used for the first tuning. The configurations returned are discussed in Section 7.3.

Evaluation of Performances

Through our experiments, we investigate the influence of both the strategy of exploration and the metric used during the LS. The algorithms are evaluated with the best configuration returned by irace at the end of Table 1: Parameter's space given to irace. The granularity and the number of patterns injected are expressed as a percentage of the size of the instance. The parameter representing the number of groups is defined as a percentage of the number of subproblems (when this parameter is equal to 1, only 1 group is created, since 0 groups is not a valid value). The space contains 388 800 configurations. In the first experiment, we use Solomon's instances of size 100 (categories R and C), which represent a total of 40 instances. Each variant is executed 30 times on the test instances. For each algorithm, the k-th run of an instance is executed with the same seed, being 10(k -1), allowing a fair comparison. We recall that the termination criterion is set to 720 seconds (i.e., 12 minutes) for all variants. For each category of instance (either R or C), we compute the average uHV obtained over the 30 runs. We perform a Friedman test on the average uHV, to know if all algorithms are equivalent, and if it is not the case, we apply a pairwise Wilcoxon test with the Bonferroni correction to know which algorithms are statistically better.

Parameter

A second experiment is performed on Gehring and Homberger's instances of size 200 (on categories R and C too). We keep the same configurations as before (parameters are scaled accordingly, in particular, the granularity and the number of patterns injected). Again, we perform 30 runs on each instance, and each algorithm is executed with the same seeds. The termination criterion is now set to 2880 seconds (i.e., 48 minutes). We perform the same statistical tests to compare the average uHV obtained.

Tuning and Parameters

Preliminary Study

For this set of experiments, we recall, that we performed a tuning of all considered algorithms with irace. The configuration space is given in Table 1. The range of each parameter is chosen to explore various parts of the space. For more details about the tuning, see Section 6.1. The tuning performed by irace on C (resp. R) instances led to the elite configurations stored in Table 2 (resp. Table 3). Only the best elite configuration for each variant is reported (i.e., the configuration of rank 1). The reference algorithm, R M OEA/D , seems more efficient with the first-best exploration strategy on both categories of instances. The same conclusion is reached for the hybrid one (H M OEA/D ). Moreover, the metric d 3 seems more appropriate to C instances, and d 2 to R instances, according to the results obtained for R M OEA/D and H M OEA/D . In fact, C instances are, in general, more simple to solve. Indeed, the closeness between customers is structurally integrated into the instance, thus a simplified metric, that is d 3 , where the waiting time between customers discriminates easily close customers, seems enough to obtain good results. The metric d 2 is more relevant for R instances, less structured, since the weights between the distance and the waiting time vary according to the aggregation of the subproblem solved. Hence, d 2 discriminates customers according to each aggregated objective. Note that, in neither case the metric d 1 is chosen, meaning that d 1 does not seem adapted to solve these instances in a bi-objective context, as expected.

Concerning the other parameters, we remark that the LS is, in most cases, applied with probability 0.10, which is relatively low, but it is a costly step in terms of execution time. In addition, the granularity remains high (with a value of the 50% or 75% closest customers considered during the search), and the choice of the metric does not seem to have a great impact on it. Moreover, the LS is much more impacting at the beginning of the algorithm, when solutions are bad, thus it may not be interesting to use it too frequently. On the other hand, the injection is applied with a much higher probability (either 0.75 or 1.00). Indeed, the injection operator needs fewer resources and is useful (in terms of diversification and intensification) throughout the execution of the algorithms. During the injection, at least 50 patterns are tentatively injected, and the maximal size of the extracted patterns is almost always set to 5 (surprisingly, it was set to 8 only for the f bd 2 variant on C instances). The value of 5 was the value chosen by the authors of PILS [START_REF] Arnold | PILS: Exploring high-order neighborhoods by pattern mining and injection[END_REF]. The crossover operator is applied with a probability below 0.50 when using the KD mechanism, i.e., in the hybrid variants (H M OEA/D , f bd 1 , f bd 2 , f bd 3 , bd 1 , bd 2 , and bd 3 ), whereas it is applied in every iteration with the reference MOEA/D. Indeed, in R M OEA/D , the crossover is the only operator that introduces diversity in the solutions, necessary to escape local optima. It is not the case with the other variants, due to the presence of the injection operator. The number of aggregations used is coherent with the values used in general (around 40 aggregations). Only R M OEA/D requires 80 aggregations on R instances to generate better and diverse solutions. Finally, the number of knowledge groups used is very variable, but low values are not used. In particular, strictly more than one groups are created. To be more precise, at least 4 groups are used for each variant. Our main hypothesis concerning the use of many groups is that the structure between two solutions (in terms of patterns extracted) in the objective space may vary significantly, even between solutions with close objective values. As a consequence, many groups are created to somehow reduce the structural gap between solutions in the same group.

More generally, we observe some differences between the tuning performed for C instances and the tuning performed for R instances. For example, the probability of injection tends to be higher for C instances, while the number of injected patterns tends to be higher for R instances. It validates our choice to tune separately the algorithms on instances R and C. Moreover, many configurations are very similar. Indeed most of them only differ by two or three parameters or use close parameters. In the following, we investigate the influence of the parameters on a subset of elite configurations, in order to see, if the configurations returned are truly local optima in the configuration space provided, and if the configurations can be improved by reducing the space.

Influence of Parameters on Elite Configurations

Given an elite configuration from Table 2 (resp. Table 3), the set of its neighbors obtained by modifying one parameter (e.g. the number of aggregations) is evaluated on the 20 generated instances of category C (resp. R). The same seed is used for all the evaluations, which is similar to a run performed by irace. The results are reported in Figure 1a (resp. Figure 1b). The (blue) boxplot on the row Initial represents the performances of the elite configuration as returned by irace. In other words, it is the reference boxplot. Each boxplot on the rows below represents the performance of the configurations obtained when the corresponding parameter is modified. We remark that only a few parameters have a noticeable impact on the configuration. The most impacting parameters are the strategy (S) used during the LS, the probability of the PMX, and the The strategy only takes two values, and according to the results obtained during the tuning, it seems coherent that changing the best strategy to the first-best strategy improves the results obtained, while the contrary deteriorates the results. Note that, here, the injection process is activated (all elite configurations have a nonnegative probability of injection).

The detailed results for the LS probability are shown in Figure 2. For each possible value of the parameter, the hypervolumes obtained with the corresponding configuration are represented. If nothing really interesting is observed when the strategy is best, with the first-best one it seems better to consider lower values (between 0.10 and 0.25).

Similarly, the uHV tends to increase when the probability of PMX increases until a plateau is reached after the value of 0.75. The influence of the crossover is even more important with the best strategy.

Concerning the probability of injection, the average uHV increases, along with the probability. In particular, when the injection is disabled (probability 0.0), we retrieve the fact that first-best strategy performs better than best strategy when no learning is used. Moreover, the injection should be applied with high probability (above 0.75) in general. For the other parameters, there are almost no variations of the uHV

Elite fbd1

Elite fbd2

Elite fbd3 Elite bd1 Elite bd2 Elite bd3 0.7 0.8 0.9 0.7 0.8 0.9 0.7 0.8 0.9 0.7 0.8 0.9 0.7 0.8 0.9 0.7 0.8 0.9

Elite bd1

Elite bd2 Elite bd3 Elite fbd1 Elite fbd2 Elite fbd3 0.7 0.8 0.9 0.7 0.8 0.9 0.7 0.8 0.9 0.7 0.8 0.9 0.7 0.8 0.9 0.7 0.8 0.9 fronts between two instances can be very different.

Reduced Configuration Space and Final Tuning

As explained in the former section, some parameters do not have much influence on the hypervolumes returned. That is why we decided to set the number of patterns injected to 100% of the size of the instance (i.e., 100 patterns are tentatively injected on instances of size 100). Note that, it was the value chosen by [START_REF] Arnold | PILS: Exploring high-order neighborhoods by pattern mining and injection[END_REF] in their own experiments. Moreover, we decided to set the number of aggregations to 40, in order to see if the number of groups used can be more precisely tuned during the tuning. According to what has been said in the previous section, we reduced the range of values for the probability of PMX to (0.25, 0.50, 0.75), and for the probability of injection to (0.50, 0.75, 1.00), being the most interesting values. Concerning the probability of LS, we changed the range to (0.10, 0.20, 0.30), in order to give more adapted parameters to irace. Finally, the granularity is also reduced to (25,50,75), since extreme values never belong to elite configurations. All the changes are reported in Table 4.

Using the reduced space of configuration, we performed a new tuning, only on the six variants f bd 1 , f bd 2 , f bd 3 , bd 1 , bd 2 , bd 3 . Indeed, R M OEA/D already had a small configuration space (it has fewer parameters to fix), and the elite configurations returned for H M OEA/D were already in the reduced configuration space defined. The budget allocated to irace is reduced from 2000 to 1000. The best elite configurations are 5. At first, we notice that the differences between the two categories of instance are more important, considering the granularity parameter (δ), and the maximum size of patterns extracted (size p ). Indeed, the granularity is always set to 50 on C instances, which seems coherent, since these instances contain clusters of customers. The granularity tends to be slightly higher on R instances, except for the variant bd 3 , which has the lowest granularity. The maximum size of the patterns is now 8 for most of the variants on R instances, meaning that bigger patterns are interesting when instances are not well structured. For C instances, the maximum size of the pattern remains coherent, since smaller patterns are enough to perform interesting moves.

Another major change concerns the number of knowledge groups used during the execution of the algorithms. The variant being set, more groups are generated on C instances than on R instances. In addition, more groups are generated with strategy best than with strategy first-best.

Concerning the other parameters, the probability of PMX is set to 0.50 for all variants, the probability of LS is set to 0.10 (except for the variant f bd 2 where it is 0.20), and the probability of injection remains similar to what we saw during the first tuning (see Table2 and Table 3).

Experimental Results and Discussion

The first set of experiments is performed on Solomon's instances of size 100 belonging to categories R and C. According to the category, we use the corresponding elite configuration in Table 5. The average uHV obtained through the 30 runs, for each category of instances, is shown in Figure 3 (plot above). The results obtained with the first-best strategy are slightly more robust, with smaller boxplots, than those obtained with the best one. Moreover, the variant f bd 2 always seems better than the other variants, especially on R instances. Among the six variants, it is the variant bd 3 that returned the worst results.

Table 6 contains the average ranking of the different algorithms on all classes of instances. In this table, the ranks of statistically best algorithms on a class of instances are put in bold. In particular, the variants f bd 2 and f bd 3 are statistically equivalent on C1 instances and are better than the other variants. On C2 instances, the variants f bd 1 , f bd 2 , f bd 3 , and bd 1 are equivalent. It is not surprising since C2 instances are the easiest ones. Moreover, it means that the strategy best with the metric d 1 , which is usually used in the literature, returns good results when the problem can be modeled by one objective. In addition, all the hybrid MOEA/Ds are better than the reference R M OEA/D , which has very low uHV. In fact, R M OEA/D gets easily stuck in local optima, in particular when there must be only one solution on the final front. On R instances, statistically, the best overall algorithm is f bd 2 . However, f bd 1 is equivalent to f bd 2 on R1 instances. As a consequence, the strategy first-best with the metric d 2 , seems to be the best combination, considering all the classes of instances of size 100. This result has already been observed since these parameters were selected for the configuration returned by irace for H M OEA/D , i.e. when the exploration strategy and the metric were left as parameters to be tuned (see Table 3). More precisely, the results obtained for all C instances (except C103 and C104) show that the Pareto front contains only one solution, minimizing both objectives. In all the other instances the optimal value of the total traveled distance is not reached by our algorithms, however on some instances (e.g. R101 and R102) it is very close. Concerning the total waiting time, it is always possible to create a solution without waiting time, and our algorithms find them. In addition, we measured the correlation between the objectives on the reference fronts obtained. The reference fronts used, and the associated correlations are available at dataItor2023. Without doing a rigorous analysis of the objective space of the instances like Castro-Gutierrez et al. ( 2011), we can not say more than some instances should be avoided when studying multi-objective optimization. Table 7 and Table 8 show more detailed results about the convergence of the algorithms f bd 2 , bd 2 , and R M OEA/D . We can remark that, as expected, the final uHV obtained by f bd 2 is higher than the one obtained by bd 2 , and R M OEA/D . Moreover, f bd 2 is able to reach 80% of the uHV of the reference front significantly faster than R M OEA/D and bd 2 on many instances. When the algorithm does not reach 80% of the uHV at the end of the execution, the maximal time allocated is used instead. It shows that using the learning mechanism with the first-best strategy speeds up the convergence process.

The second set of experiments is performed on Gehring and Homberger's instances of size 200. Again, we kept only the instances of categories R and C. The configurations used are the same as the configurations used for the instances of Solomon, except that the granularity is doubled, as well as the number of patterns injected. Indeed, both of these parameters are directly linked to the size of the instance. In Figure 3, we clearly observe that for the benchmark of size 200, the variants using the first-best strategy are much more robust than the variants using the best one. Moreover, when the exploration strategy is set, there is almost no change when the metric is modified. One hypothesis should be that the given values for the granularity are too high for these instances. Considering the Table 6, the variant f bd 3 is statistically better on instances C1 (with tight time windows), while the variants f bd 2 and f bd 3 are statistically equivalent on instances C2 (with wide time windows). Moreover, these two algorithms remain better than all the others. On R1 instances, we do not reach the same conclusion. Indeed, the variant f bd 1 is statistically better now. However, on R2 instances, the three variants f bd 1 , f bd 2 , and f bd 3 are statistically equivalent and the results returned by these three variants are very close. It is important to notice that with a specific tuning on instances of size 200, we may have reached different conclusions. That being said, strategy first-best is always better than strategy best. Concerning the metrics, we finally observe that they have a relatively low impact on instances of a bigger size. We meet again the result obtained in Section 7.2, where the metrics did not seem to have an impact on the elite configurations. The reference fronts, the correlations, and additional convergence results are publicly available at dataItor2023.

Conclusion

This paper addresses the bi-objective vehicle routing problem with time windows (bVRPTW), a routing problem where both the total traveling cost and the total waiting time of drivers have to be minimized. We chose to solve the bVRPTW with a hybrid MOEA/D [START_REF] Legrand | Improving moea/d with knowledge discovery. application to a bi-objective routing problem[END_REF] where knowledge discovery mechanisms are used to extract patterns, sequences of customers here, and inject them into solutions. Indeed, MOEA/D is known to obtain pretty good performance when solving a bi-objective combinatorial optimization problem. This hybrid MOEA/D framework is modular and both strategical components and parameter values can be fixed. In the mutation step of MOEA/D, local search algorithms are widely used to 2.3 2.3 1.4 1.0 1.9 1.4 2.6 1.7 fbd 3 2.2 2.9 2.5 2.7 1.1 1.7 2.4 2.6 bd 1

3.6 3.3 5.1 4.7 4.0 5.3 6.1 6.5 bd 2 4.9 4.2 3.9 5.0 5.3 5.6 5.0 6.5 bd 3 4.4 5.9 6.1 7.0 5.7 6.5 4.3 5.0 intensify the search in specific regions and find the best local optimum solutions. While routing problems generally used a local search with a best improvement strategy [START_REF] Hoos | Stochastic local search: Foundations and applications[END_REF] and with a distance metric to reduce the neighborhood size [START_REF] Schneider | Designing granular solution methods for routing problems with time windows[END_REF], we adapted from scheduling the first-best improvement strategy and enhanced the distance metric with the waiting time information of customers. Experiments have been conducted on two benchmarks of the literature. We showed that the first-best strategy was better adapted than the best strategy, in most cases, to solve bVRPTW instances. These results are supported by a rigorous tuning of the hybrid MOEA/D and a deep analysis of the influence of its parameters. It emerges from this analysis that only a few parameters are really impacting, like the probability of injection, the probability of crossover, and the probability of local search. Moreover, it is unnecessary to consider the interval [0, 1] entirely, but only a subset containing the most interesting values. The other parameters should either evolve during the execution of the algorithm, in order to adapt to the instances solved (like the granularity, the number of groups, or the number of patterns injected), or be directly set (like the number of aggregations and the maximum size of patterns extracted). This study also allowed us to assess the configurations returned by irace, in particular on the choice of the exploration strategy and the distance metric. Nevertheless, the granularity value returned is higher than expected. Since it does not highly impact the performances of the algorithms, we advise using a smaller value (like 0.25 or 0.50 instead of 0.75) to reduce the memory allocated during the execution. Moreover, the analysis of the neighborhood of elite configurations shows that the configuration space i.e., the values of parameters, was

Algorithm 3 :

 3 The f irst-best procedure. Input: A solution x and a neighborhood operator N Output: A local optimum 1 improve ← T rue 2 while improve do 3 improve ← F alse 4 customers ← shuffle([1 . . . N ]) 5 for customer ∈ customers do 6 moves ← generate moves(customer, N ) 7

Fig. 2 :

 2 Fig. 2: Influence of the LS probability on instances of (a) category C and (b) category R. The red dot represents the mean uHV.

Fig. 3 :

 3 Fig. 3: Results obtained on the different categories of instances, of (left) Solomon's benchmark for size 100, and (right) Gehring and Homberger's benchmark for size 200.

Table 2 :

 2 Best elite configurations returned by irace when tuning on C instances. The symbol "-" reflects the absence of the parameter during the tuning.Parameter R M OEA/D H M OEA/D fbd 1 fbd 2 fbd 3 bd 1 bd 2 bd 3

	M	40	40	20	40	40	20	40	40
	δ (%)	50	75	50	75	50	75	50	50
	p pmx	1.00	0.50	0.50 0.25 0.25 0.50 0.10 0.25
	p ls	0.10	0.10	0.10 0.10 0.10 0.10 0.25 0.10
	S	first-best	first-best	-	-	-	-	-	-
	d	d 3	d 3	-	-	-	-	-	-
	k G (%)	-	50	25	10	25	100	10	50
	size p	-	5	5	8	5	5	5	5
	N Inj (%)	-	100	50	50	75	50	50	75
	p inj	-	0.75	1.00 1.00 1.00 1.00 1.00 1.00

Table 3 :

 3 Best elite configurations returned by irace when tuning on R instances. The symbol "-" reflects the absence of the parameter during the tuning.Parameter R M OEA/D H M OEA/D fbd 1 fbd 2 fbd 3 bd 1 bd 2 bd 3

	M	80	40	20	60	40	40	40	20
	δ (%)	75	75	50	75	75	50	50	50
	p pmx	1.00	0.50	0.50 0.50 0.25 0.50 0.25 0.50
	p ls	0.10	0.10	0.25 0.10 0.10 0.10 0.10 0.10
	S	first-best	first-best	-	-	-	-	-	-
	d	d 2	d 2	-	-	-	-	-	-
	k G (%)	-	100	25	100	75	100	10	100
	size p	-	5	5	5	5	5	5	5
	N Inj (%)	-	100	75	75	75	75	75	75
	p inj	-	0.75	1.00 1.00 0.75 0.75 0.75 0.75
	probability of injection.

Table 4 :

 4 Reduction of the configuration space. The parameters with a modified range are underlined in red. The space now contains 1 458 different configurations.

	Parameter

Table 5 :

 5 Best elite configurations, from the reduced configuration space, returned by irace on instances of (a) category C, and (b) category R.Param. fbd 1 fbd 2 fbd 3 bd 1 bd 2 bd 3

			(a) Category C.					(b) Category R.
								fbd 1 fbd 2 fbd 3 bd 1 bd 2 bd 3
	δ (%)	50	50	50	50	50	50	75	75	75	50	50	25
	p pmx	0.50 0.50 0.50 0.50 0.50 0.50	0.50 0.50 0.50 0.50 0.50 0.50
	p ls	0.10 0.10 0.10 0.10 0.10 0.10	0.10 0.20 0.10 0.10 0.10 0.10
	k G (%) 50	50	100 100 100 100	50	25	50	100 100	50
	size p	5	5	5	5	2	2	8	8	8	8	8	2
	p inj	0.75 1.00 1.00 1.00 0.75 0.75	1.00 1.00 0.75 0.75 1.00 0.75
	reported in Table

Table 6 :

 6 Average ranking of all the algorithms on the different categories of instances. Bold ranks are statistically equivalent.

	Size	100		200
	Category	C	R	C	R
	Class	C1 C2 R1 R2 C1 C2 R1 R2
	R M OEA/D 7.0 7.0 6.8 5.3 7.0 4.8 6.8 4.0
	fbd 1	3.7 2.4 2.3 2.3 3.0 2.9 1.1 1.7
	fbd 2

Table 8 :

 8 Detailed results for algorithms f bd 2 , bd 2 , and R M OEA/D on R instances. From left to right: the average size of the front, the average uHV, and the average time to reach 80% of the reference uHV.

		fbd 2		bd 2		R M OEA/D
	Inst.	|F | uHV Time (s) |F | uHV Time (s) |F | uHV Time (s)
	R101 72.2 0.915	65.0	70.6 0.911	229.7	47.7 0.907	73.4
	R102 45.6 0.923	137.8	45.0 0.911	296.8	27.1 0.876	232.7
	R103 23.7 0.937	277.8	23.4 0.899	457.0	20.3 0.724	721.1
	R104 6.0 0.915	341.2	6.3 0.881	533.4	6.4 0.499	721.2
	R105 12.3 0.941	166.2	13.9 0.906	299.6	7.7 0.862	438.6
	R106 6.2 0.894	258.8	8.0 0.878	406.5	5.5 0.763	700.1
	R107 6.0 0.851	496.2	5.1 0.828	602.9	6.8 0.524	720.9
	R108 1.8 0.880	358.2	1.5 0.851	532.8	2.5 0.389	720.9
	R109 1.2 0.915	259.5	1.4 0.877	450.8	1.5 0.727	680.3
	R110 1.1 0.843	368.9	1.3 0.849	496.6	1.3 0.637	717.1
	R111 1.9 0.850	423.1	1.9 0.809	583.4	2.4 0.517	720.7
	R112 1.0 0.845	480.2	1.0 0.804	559.9	1.0 0.387	720.4
	R201 52.1 0.854	230.3	52.2 0.843	376.4	33.6 0.814	525.9
	R202 44.7 0.857	311.6	42.1 0.849	457.3	25.8 0.824	515.1
	R203 34.5 0.859	338.4	31.2 0.819	600.0	21.3 0.797	613.9
	R204 13.2 0.806	548.7	11.9 0.752	679.7	8.1 0.719	706.6
	R205 20.6 0.839	413.7	22.3 0.833	509.6	14.4 0.798	637.2
	R206 22.5 0.867	349.4	22.9 0.839	534.9	13.0 0.821	581.8
	R207 19.4 0.828	467.1	16.4 0.805	622.5	10.7 0.776	680.6
	R208 5.5 0.808	541.6	6.5 0.790	605.5	3.8 0.720	718.8
	R209 12.0 0.808	543.2	12.6 0.776	648.0	10.0 0.755	712.1
	R210 20.4 0.855	368.3	20.5 0.830	555.0	12.3 0.791	649.3
	R211 1.0 0.786	548.7	1.0 0.768	618.9	1.0 0.732	673.5

Acknowledgements

This work has been supported by the French National Research Agency through the AI PhD@Lille program (ANR-20-THIA-0014). This support is gratefully acknowledged.

Elite bd1

Elite bd2

Elite bd3 when the value of the parameter changes. In theory, it means that we could choose any possible value for the remaining parameters, without changing the uHV too much. Moreover, it also explains why the set of elite configurations returned by irace for one variant, contains in general very different values for those parameters. In particular, considering the number of groups, no value seems more interesting than another one to obtain good results. It is very dependent on the other parameters (and in particular, the number of aggregations). We think, that this parameter should also depend on the instance considered since the Pareto well defined since a single parameter cannot significantly change the performance obtained on the tuning instances.

In future works, we will investigate the adaptive side of the algorithm by automatically updating some parameters during the execution of the algorithm. Besides, we will tackle different variants of the bVRPTW to validate or not our conclusions when the objectives are modified.


Lucas Terrei
email: lucas.terrei@univ-lorraine.fr

Davood Zeinali

Zoubir Acem

Véronique Marchetti

Paul Lardet

Pascal Boulet

Gilles ; Parent

Zeinali, Davood; LEMTA Lucas ; Lemta

Véronique; CSTB Lemta Marchetti

Cstb Boulet

Lemta

Experimental study of spruce wood reaction to fire in single burning item test

Keywords: Single Burning Item, Cone calorimeter, Wood degradation Single Burning Item, cone calorimeter, wood, degradation 22

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION 23

In Europe, building materials are subjected to standardized tests which allow classification according 24 to their reaction to fire. [START_REF] Lemta | The single burning item (sbi) test method-a decade of development 357 and plans for the near future[END_REF] One of the most widely used tests in reaction to fire is the Single Burning Item 25 (SBI). [START_REF] Standard | 13823: 2013 reaction to fire tests for building products-building products excluding 359 floorings exposed to the thermal attack by a single burning item[END_REF] The SBI test is an intermediate-scale test between bench-scale tests such as Fire Propagation Ap-26 paratus (FPA) test and full-scale test such as ISO Room Corner test. [START_REF] Iso | International standard-fire tests-full-scale room test for surface products[END_REF] This test simulates a single burning For the normative test, the two panels exposed to the burner do not have the same dimension. There 106 are one short panel (50 cm × 150 cm) and one large panel (100 cm × 150 cm). This arrangement makes 107 it possible to observe the horizontal spread of the flame during the test. Nevertheless, it was shown by 108 Zeinali et al. 7 that the flame spread on the long MDF panel never exceeds 50 cm, i.e. less than the short 109 panel width. So, for this study, two identical panels having a 50 cm width were used. A similar configura-110 tion was also chosen by Chaudhari et al. on PMMA panels. 17 Additionally, in the present configuration, 111 the remainder of the wall on the large panel side was completed by a 50 cm wide calcium silicate panel 112 to remain in the same configuration to the original SBI test and ensure an air supply similar to the tests 113 carried out in previous works. For the sake of clarity, the labels "large panel" and "short panel" denomi-114 nations will be still used in the paper as shown in Fig. 1. The spruce panels were exposed to burner HRR 115 values of 15, 20, 30 and 50 kW for exposure times of 15, 20 and 30 minutes.

116

Temperatures inside the wood were measured with K-type thermocouples during the tests. However, de-117 pending on how the thermocouples are implanted in the wood sample, the measured temperature can vary 118 by up to 400 • C for some experimental conditions, as shown in. [START_REF] Terrei | In-depth wood temperature 394 measurement using embedded thin wire thermocouples in cone calorimeter tests[END_REF][START_REF] Fahrni | Correct temperature measurements in fire exposed 400 wood[END_REF][START_REF] Fahrni | Investigation of different temperature measure-403 ment designs and installations in timber members as low conductive material[END_REF] A precise measurement method 119 was developed recently, [START_REF] Terrei | In-depth wood temperature 394 measurement using embedded thin wire thermocouples in cone calorimeter tests[END_REF] consisting in embedding very thin (0.1 mm diameter) wire thermocouples in- side the sample, along the isotherms parallel to the exposed surface. This method was validated at the 121 small scale with cone calorimeter tests and is implemented here, at the SBI scale. For that, the wood panel 122 was cut perpendicularly to the exposed surface, and then square-grooves with 0.2 mm depth and 0.2 mm 123 width were precisely machined at the desired locations. Then, wire K-type thermocouples composed of 124 alumel and chromel wires were welded by autogenous electric welding and put in the square-grooves.

125

Finally, wood panel was glued with Melamine-Urea-Formaldehyde (MUF). More details concerning the 126 protocol of the thermocouple implantation in the wood sample are presented in. [START_REF] Terrei | In-depth wood temperature 394 measurement using embedded thin wire thermocouples in cone calorimeter tests[END_REF] At the scale of SBI, 127 implanting the thermocouple wires in this way takes four hours per sample from grooving to bonding 128 and thermocouple welding. The temperature was only studied for the short panel. Figure 2 suggested by Reszka and Torero, [START_REF] Reszka | In-depth temperature measurements in wood exposed to intense radiant 406 energy[END_REF] the sheathed thermocouples were spaced by 20 mm in order to limit their mutual disturbance.

138

The cone calorimeter 139 As introduced by previous studies, [START_REF] Hansen | Prediction of time to flashover in the iso 9705 room corner test based 378 on cone calorimeter test results[END_REF][START_REF] Delichatsios | Application of upward flame spread for the prediction of sbi and iso room 380 corner (and parallel wall) experiments and classification[END_REF][START_REF] Hakkarainen | Application of a one-dimensional thermal flame spread model 383 on predicting the rate of heat release in the sbi test[END_REF] comparisons between the cone calorimeter [START_REF] Iso | Reaction-to-fire tests -heat release, smoke production and mass loss rate -part 1: 363 Heat release rate (cone calorimeter method) and smoke production rate (dynamic measurement)[END_REF] and the SBI 2 140 are very useful for the understanding of the fire behavior of materials. The HRR evolution or the time-141 to-ignition are among the most often used data to identify the correlation between the two scales. In the 142 present work, the char depth positions in the two scales were compared over time. Following the standard 143 ISO 5660-1, the wood samples in the cone calorimeter tests were wrapped with two layers of aluminum 144 foil, and the distance between the heater and the tested samples was 25 mm. For these experiments, the 145 cone heater and wood samples were oriented vertically like for SBI tests, without a pilot flame or spark.

146

The spruce wood material of the cone samples (10 cm × 10 cm × 5 cm) was identical to that tested at (Gardon type from Captec). This heat flux sensor was placed in holes with a diameter equal to that of 158 the fluxmeter (i.e., 2.5 cm), previously drilled on the short panel at various locations. Figure 3 shows the 159 locations of the gauge and the corresponding mapping of the measured heat fluxes for different burner 160 powers. The free holes were blocked with plugs from the same materials as the panel to prevent any gas 161 flow through them. The heat flux was measured on the small panel at a total of 22 points with locations 162 ranging from 15 to 105 cm along the height and from 0 to 40 cm along the width. For each burner 163 power and measurement point, the heat flux was recorded for 2 minutes. When the burner is ignited, it is The heat flux maps illustrated in Fig. 3 show that the heat flux received at the wall increases locally 169 with the power of the burner. For a 50 kW burner power, the measured heat flux is at least 30 kW.m -2 170 over a large area (20 cm wide and 85 cm high). Beyond 105 cm in height and 30 cm in width, the heat flux was not measured (except in position 22) but the heat flux does not exceed 5 kW.m -2 there. These The measurements indicate that the heat flux received by the wooden panel is systematically higher 182 than that received by the incombustible panel for the same burner powers (see Fig. 4). For a 15 kW burner 183 power, the heat flux is increased by 85 % for the wood panel, compared to the inert panel, but only 3 % for 184 the 50 kW burner power. This shows that the difference tends to decrease as the heating power increases.

185

Indeed, for a low power of the burner, the heat flux received on the panels is mainly radiative since the 186 burner flame is too small to stick to the panel surface (more precisely at the fluxmeter's location). For For higher burner powers, the flame of the burner is large enough to stick to the panel surface so a 192 convective flux coming from the burner flame is also present. In this case, the pyrolysis gases feed the 193 burner flame but its thickness does increase any further. Moreover, the combustible panels in this case also 194 experience the so-called "blowing effect", i.e., reduction of convective heat transfer to the surface due to 195 the opposite flow of pyrolyzate mass. [START_REF] Zeinali | Experimental study of 366 corner fires part i: Inert panel tests[END_REF][START_REF] Drysdale | An introduction to fire dynamics[END_REF][START_REF] Quintiere | Fundamentals of fire phenomena[END_REF] Accordingly, as the burning rate of the panels increases with 196 higher burner powers, the convective heat transfer to the surface decrease, just as reported in. 7,[START_REF] Zeinali | Flame spread and fire behavior in a corner configuration[END_REF] This As illustrated in Fig. 5, it is possible to define three zones of degradation: (A) the zone next to the For a given exposure time, the propagation of the char front at the surface increases with the burner 219 power. There is also an asymmetry between the degradation patterns of the short panel and those of the 220 large panel. This observation is also confirmed by the char depth, which generally shows that the short 221 panel is more degraded than the large one. This observation concerning the in-depth temperatures and 222 heat fluxes was also observed by Zeinali et al. [START_REF] Zeinali | Experimental study of 366 corner fires part i: Inert panel tests[END_REF] and Zhang et al., [START_REF] Zhang | Assessment of fire dynamics simulator for heat flux 368 and flame heights predictions from fires in sbi tests[END_REF] and is expected to be due to the air 223 supply pattern in the setup which favors flame spread on the short panel. [START_REF] Zeinali | Flame spread and fire behavior in a corner configuration[END_REF] Indeed, this phenomenon 224 leads to higher temperatures of 100 • C 5 as well as 1 to 5 kW.m -2 larger heat fluxes [START_REF] Zhang | Assessment of fire dynamics simulator for heat flux 368 and flame heights predictions from fires in sbi tests[END_REF] on the short side. embedded thermocouples temperature measurements). [START_REF] Terrei | In-depth wood temperature 394 measurement using embedded thin wire thermocouples in cone calorimeter tests[END_REF][START_REF] Tran | Burning rate of solid wood measured in a heat release rate calorimeter[END_REF] The results are presented in Fig. 11. The in-depth char front positions determined via direct measurements are in good agreement with 259 those estimated by the 300 • C isotherm measurements made using the embedded thermocouples, just 260 as observed by Terrei et al. [START_REF] Terrei | In-depth wood temperature 394 measurement using embedded thin wire thermocouples in cone calorimeter tests[END_REF] The char front reached the depths of 10.5 and 18 mm after the exposure 261 times of 15 and 30 min, respectively. It appears that the propagation of the char front slows down with 262 time according to a power law, as showed in Fig. 11b. The curve shows that the charring rate is not 263 constant during the experiment. It starts at 1.5 mm.min -1 and decreases rapidly to 0.7 mm.min -1 after 264 10 minutes of exposure. Subsequently, the charring rate becomes more or less constant, staying between 265 0.6 and 0.7 mm.min -1 , and such values comply well with those obtained in previous studies performed 266 with spruce wood through cone calorimeter tests [START_REF] Terrei | In-depth wood temperature 394 measurement using embedded thin wire thermocouples in cone calorimeter tests[END_REF] and standard fire resistance tests. 25, 26 267 The in-depth char front propagation was measured for the entire width of the short panel as well as the 268 large panel by cutting the panels at every 5 cm height after each test and making direct measurements regarding the different thermal exposure on the small and large panel. Indeed, the average char front 276 depths in the small panel are always slightly higher than those in the large panel, indicating that the short 277 panel receives a slightly stronger thermal exposure, and this holds true for the entire panel width. The 278 char front contours in Fig. 14 are based on experimental measurement (see protocol in Fig. 12). Thus, In Fig. 15a, where the reference of heat flux is based on inert panels, the char depth positions in SBI The SBI was used to study the fire behavior of spruce wood at various burner powers and exposure 318 times, i.e. beyond those defined in the standard of SBI. [START_REF] Standard | 13823: 2013 reaction to fire tests for building products-building products excluding 359 floorings exposed to the thermal attack by a single burning item[END_REF] Correspondingly, a total of 22 tests were per- an incombustible panel. For the highest burner powers, the total heat flux measured over the combustible 326 panels was predominantly due to the burner power, rather than the burning of the panels themself. This For the charring rate, a quasi-steady-state plateau between 0.6 and 0.7 mm.min -1 was found while the 340 propagation of the charring front on the exposed surface increased according to the exposure time and the 341 burner power. In addition, it was found that the surface area of the highly degraded region varies linearly 342 with the burner power for a fixed exposure time.

343

The comparison of the char front depth in the cone calorimeter and the SBI setup showed a good agree-

  bottom of a room corner (like a burning trashcan) using a propane burner to test construction products and to classify them according to their fire reaction. The performance is evaluated over a period 29 of 20 minutes with visual and non-visual observations. During the experiment, visual observations are benefit of such an implantation.

Figure 1 :

 1 Figure 1: SBI experimental setup.

  shows the 129 position of the thermocouples in the wood sample. The thermocouple junctions (red points in Fig. 2) are 130 placed at 9 cm from the corner line and 45 cm from the bottom of the panel, at a location where the panel 131 should be severely affected by the burner flame.

  132

Figure 2 : 6 F

 26 Figure 2: Diagrams showing the in-depth thermocouple measurement points (red circles) inside the wood panel. The top view of the panel is shown to the left, while the side view of the measurement panel is shown to the right.

  147 intermediate scale. The incident heat fluxes for the cone calorimeter tests were chosen in order to match 148 the SBI heat fluxes measured both with the wood and the incombustible panels.

  presents the results of the SBI tests conducted with incombustible calcium silicate 151 panels (average density: 870 ±50 kg.m -3 and 12 mm thickness), followed by the test results of spruce 152 wood panels in SBI and the cone calorimeter, discussing the heat fluxes, in-depth temperatures, char front 153 on the exposed surface and char front depth in the samples. 154 Total incident heat fluxes on incombustible panels in SBI 155 First, an inert calcium silicate panel was exposed to the burner powers equal to 15, 20, 30 and 50 kW, 156 while the total incident heat fluxes were measured with a water-cooled (around 12 • C) heat flux gauge 157

  for approximately 30 seconds before reaching a steady state. As a result, the average 165 of the flux value and its standard deviation were only calculated for the last 90 seconds. The acquisition 166 frequency being 5 Hz, the average was taken over approximately 365 measurements. Figure3presents a 167 map of the heat fluxes obtained for each burner power.

Figure 3 :

 3 Figure 3: Contour plots of total heat fluxes obtained using burner powers of 15, 20, 30 and 50 kW on the incombustible short panel. The measurements were made using a water-cooled heat fluxmeter at the locations shown on the diagram to the right.

172Figure 4 :

 4 Figure 4: Comparison of the total incident heat fluxes measured on an incombustible panel (squares) and that with a spruce wood panel (triangles) as a function of the burner power.

187a

  combustible panel, the burning of the panel produces an additional convective heat flux on the panel 188 surface, as the pyrolysis gases feed and widen the burner flame. Therefore, the heat fluxes received by 189 the combustible panels are significantly larger than those received by the inert panels. This explains the 190 large differences in the irradiance of the two types of surface for the 15 and 20 kW tests.191

  the convergence of the heat fluxes of combustible and incombustible panel tests as the burner 198 power increases at the position where the heat flux was measured. 199 Surface char front propagation on the wood panels 200 After each test, the wood panels were extinguished with water in order to stop the wood combustion 201 and degradation. The inspection of the samples after the tests shows that the power of the burner plays 202 an important role in the degradation of the wood. Figure 5 presents spruce wood samples after tests at 203 different burner powers for an exposure time of 15 minutes. For the same exposure time, the quantity of 204 char increases with the imposed power. 205 (a) 15 kW. (b) 20 kW. (c) 30 kW. (d) 50 kW.

Figure 5 :

 5 Figure 5: Photos of test specimens after 15 minutes of testing with burner powers ranging from 15 to 50 kW.

Figure 6 :

 6 Figure 6: Degradation zones of wood samples for tests conducted with a 20 kW burner power. The error bars indicate the magnitude of one standard deviation across the tests.

Figure 6

 6 Figure 6 also indicates the standard deviation in the positioning of the boundaries of zones A, B and

  30 minutes exposure time.

Figure 7 :

 7 Figure 7: Evolution of the surface char front in zone A for different burner powers with exposure times of 15 min (left) and 30 min (right). The error bars indicate the magnitude of one standard deviation across three tests.

  225

Figure 8 Figure 8 :In-depth panel temperatures 232 Figure 9 Figure 9 :Figure 10 : 15 F

 88232991015 Figure8shows how the area of zone A increases as a function of the imposed burner power.226

  In-depth charring rate.

Figure 11 :

 11 Figure 11: Comparison of the in-depth char front propagation (left) and the deduced charring rate (right), determined either on direct measurements of the char layer thickness (square symbols) or estimated by the 300 • C isotherm measurements made using the embedded thermocouples (circle symbols), with a burner power of 30 kW.

269

  Fig. 13 and 14 in the form of contour plots for tests with burner powers of 10, 20, 30, and 50 kW and 271

Figure 12 :

 12 Figure 12: Steps for cutting the burning wood panel and measuring the location of the char front.

Figure 13

 13 Figure 13 presents a comparison of the spatial profiles over the width of the char front depths on the

Figure 13 :

 13 Figure13: Comparison of the spatial profile over the width of the char front depths in the small and large panels after 30 minutes exposure to a 30 kW burner power.

279Figure 14 :

 14 Figure 14: Contour plots of char front depths on a spruce wood panel (left) and contour plots of total heat fluxes on an incombustible panel (right) for burner powers of 15, 20, 30 and 50 kW.

  Comparison based on heat fluxes measured with the inert panels. Comparison based on heat fluxes measured with the wood panels.

Figure 15 :

 15 Figure 15: Comparison of the char front position between the cone and the SBI as a function of time, for different heat fluxes measured at position 9 shown in Fig. 3. The error bars indicate the magnitude of deviation across three repeatability tests, and the symbols indicate the average values.

303

  and cone tests are in good agreement for the highest heat flux (showing differences less than 2 mm), 304 although the two tests involve very different scales. Nevertheless, for the lowest heat fluxes (20 and 28 305 kW.m -2 ), the char front positions are no more comparable since for these burner powers the heat flux 306 measured on a wood panel differs significantly from the heat flux measured on an inert panel, as shown 307 in Fig. 4. This explains the extended char depth in the SBI tests. In Fig. 15b where the reference of heat 308 flux is based on wood panels, the comparison of the char front positions obtained at the SBI scale and at 309 the cone scale shows better agreement between the two scales. The deviations for each heat flux are less

319

  formed with measurements including char depth, heat flux, in-depth temperatures with very thin thermo-320 couples (0.1 mm diameter) embedded inside wood panels as well as position of the char on the exposed 321 surface. The heat fluxes measured with inert panels provided data about the distribution of total inci-322 dent heat fluxes over the panels for different burner powers. Accordingly, it was quantified how the heat 323 flux distribution over the surface of the panels is non-uniform, with heat fluxes ranging from 1.5 to 50 324 kW.m -2 . The measured heat fluxes with wood panel were up to 85 % higher than the heat fluxes with 325

327 22 331

 22 causes the heat fluxes to converge to those observed with incombustible panels. Indeed for these burner 328 powers, the flame has already reached its maximum thickness and thus the heat flux cannot increase 329 further. Moreover, high burner powers induce strong opposite flow of pyrolyzate mass over the surface, 330 causing a so-called 'blowing effect' that leads to convective cooling. As demonstrated in some recent studies, the in-depth temperature can be underestimated when the ther-332 mocouple is not fixed in an appropriate way. Very thin thermocouples implanted along the isotherms 333 were used in this study to quantify this underestimation and to accurately measure the in-depth temper-334 atures. As a result, it was found that the thermocouple wires implanted along the isotherms are faster to 335 detect the rise of temperatures compared to the sheathed thermocouples placed in holes drilled from the 336 back side, with up to 75 • C difference in the temperatures obtained using the two methods. Moreover, 337 the position of the char front deduced from the 300 • C isotherm was found to comply very well with the 338 direct measurements.

  339

  higher by nearly 100 % at an exposure time of 15 min, albeit the difference decreases to approximately 50[7] D. Zeinali, S. Verstockt, T. Beji, G. Maragkos, J. Degroote, and B. Merci, "Experimental study of 371 corner fires part ii: Flame spread over mdf panels," Combustion and Flame, vol. 189, pp. 491-505,

	372
	2018.
	F o r F o r F o r
	P e P e P e
	e r e r e r
	R e R e R e
	v v v
	i e i e i e
	w w w
	346
	22

344

ment for the higher heat fluxes. However, for the lower heat fluxes measured with an inert panel, there is 345 a significant difference between the char front depths in SBI and cone tests, namely those of the SBI are

https://mc.manuscriptcentral.com/jfsciences Journal of Fire Sciences

https://mc.manuscriptcentral.com/jfsciences

than 2 mm between the SBI and the cone tests. These comparisons show that it would be more accurate 311 to consider the heat flux with the wood panel for studying the scaling effect. However, results also show 312 that for heat fluxes measured on the inert panel greater than 43 kW.m -2 , it is possible to consider that the 313 positions of the char front will be similar between the two scales. At lower heat fluxes, the modification