Multi-Objective Model and Variable Neighborhood Search Algorithms for the Joint Maintenance Scheduling and Workforce Routing Problem

: This paper addresses a problem faced by maintenance service providers: performing maintenance activities at the right time on geographically distributed machines subjected to random failures. This problem requires determining for each technician the sequence of maintenance operations to perform to minimize the total expected costs while ensuring a high level of machine availability. To date, research in this area has dealt with routing and maintenance schedules separately. This study aims to determine the optimal maintenance and routing plan simultaneously. A new bi-objective mathematical model that integrates both routing and maintenance considerations is proposed for time-based preventive maintenance. The ﬁrst objective is to minimize the travel cost related to technicians’ routing. The second objective can either minimize the total preventive and corrective maintenance cost or the failure cost. New general variable neighborhood search (GVNS) and variable neighborhood descent (VND) algorithms based on the Pareto dominance concept are proposed and performed over newly generated instances. The efﬁciency of our approach is demonstrated through several experiments. Compared to the commercial solver and existing multi-objective VND and GVNS, these new algorithms obtain highly competitive results on both mono-objective and bi-objective variants.


Introduction
Large-scale systems are essential to today's industry, such as aeronautics, railways, telecommunications, etc. The management of such systems is often complicated and requires an adapted maintenance strategy. Maintenance is indeed a vital primary service in industries, especially when failures are likely to impact personnel safety and cause important environmental damages. It can also have a major impact on profitability.
Maintenance aims to retain equipment in its operating condition or restore it to a previous state enabling its required functions. To perform maintenance activities, it is crucial to effectively manage a crew of operators and sequence their visits over a planning horizon. Two main categories of maintenance can be distinguished [1]: (1) Corrective Maintenance (CM) that is performed after the failure. It includes actions such as repair or replacement. (2) Preventive Maintenance (PM), which occurs before the failure and intends to reduce the risks of unforeseen breakdowns. It can itself be divided into two categories: 1.
Time-Based Maintenance (TBM) that includes the periodic replacement and the age replacement policy (ARP). 2. Condition-Based Maintenance (CBM) that is based on the inspection of the unit or system before the intervention.
In the first case of TBM, the preventive replacement is performed at regular time intervals. For the ARP, the unit is replaced depending on its age or its remaining useful life [1]. CBM is based on the inspection of the unit or system before the intervention. It focuses on predicting the health condition of the system using information collected from sensors.
The companies mainly outsource their maintenance operations to a service provider. This maintenance provider is required to ensure good quality maintenance services at the lowest overall cost. Figuring out the right time and manner to carry out the maintenance are therefore major concerns. The planning of maintenance operations on geographically dispersed machines subjected to random failures is, therefore, a complex problem faced by service providers. In the industry, machines are exposed to random failures that can lead to significant penalties. The opportune maintenance times are sometimes defined by the equipment manufacturer, but more often, the companies entrust their maintenance providers to define the maintenance schedule. The service provider tries to reduce the number of maintenance operations as much as possible, which can be accomplished by maximizing the periods between two successive visits to the same machine. At the same time, he seeks to avoid the huge costs resulting from carrying out corrective maintenance operations on broken machines. The challenge is therefore finding the optimal number of visits to perform maintenance operations. Mathematical maintenance policies are used to develop an effective preventive maintenance plan. The objective of such policies is to design a maintenance schedule including both preventive replacement and corrective replacement. The maintenance models dealing with these policies are probabilistic due to the uncertainty of the mechanism that causes failure. When machines are geographically distributed, it is necessary to sequence the visits and determine the best routing-maintenance policy. Cost, availability, reliability, and maintainability are the benchmarks for evaluating maintenance decisions, similar to cost, quality, and time objectives in logistics. The maintenance company would be interested in finding the best compromise between services and maintenance costs on one hand and operational and transport costs on the other hand. The objective is to jointly consider the technical aspects of maintenance and the organizational aspects of operations management. To consider these two aspects, many works in the literature consider both maintenance and routing problems [2][3][4][5]. To handle this optimization problem, it is necessary to simultaneously investigate the maintenance scheduling and vehicle routing problem (VRP). The vehicle routing problem with time windows (VRPTW) has been widely used in the literature as the main variant of VRP for planning and scheduling maintenance operations [4][5][6][7].The problem addressed in this paper consists of determining the best routing maintenance policy when planning operators' schedules to perform maintenance operations on geographically distributed machines subjected to random failures.
Herein, the main contributions of this work are summarized as follows.
(1) A new bi-objective model is proposed, aiming to minimize both maintenance and transport costs in the case of time-based preventive maintenance. The first objective minimizes the total travel cost related to technicians' routing. The second objective can either minimize the total preventive and corrective maintenance cost or the failure cost. A penalty cost is incurred when the maintenance activities are performed after the optimal time interval. The present paper comes with the novelty of introducing a nonlinear and uncertain failure cost that uses information from equipment degradation. Moreover, the failure and maintenance costs adopted have not been used in a multiobjective study with the routing objective. We also include a continuous-time for the last restoration in the second objective of the model. (2) Variable Neighborhood Descent and General Variable Neighborhood Search variants have been designed and implemented to solve the mono-objective and the bi-objective problems. The proposed multi-objective VND and GVNS algorithms are based on the Pareto dominance strategy and start with a unique initial solution. The algorithm MOVND/P is designed to be an intensified local search component of the GVNS algorithms, whereas MOVND/PI and MOGVNS/P are intended to solve the multiobjective problem. They use a proposed multi-objective best improvement strategy called MOBI/P and new adaptations of VNS mechanisms in the multi-objective context. (3) Extensive experiments demonstrated that the proposed algorithms outperformed the literature algorithms. We also test some other GVNS variants for comparison reasons. An analysis is conducted to clearly show the differences between the performance of the algorithms proposed and the literature algorithms. It permits measuring the influence of the proposed mechanisms to improve the results.
This study will be presented as follows: Section 2 sets a summary of the corresponding literature. Section 3 delineates the problem and sets out its mathematical formulation. A description of the proposed solution approaches and their implementation details are given in Section 4. Section 5 then presents experimental results. At last, Section 6 lays out concluding observations and perspectives for upcoming research.

Literature Review
The majority of existing studies tackling routing and maintenance problems dealt separately with each of them. Two principal research streams can be distinguished in the literature: First stream: The routing is used to schedule maintenance operations. Second stream: Maintenance and routing are viewed as an integrated problem by considering both their specific features.
The first stream of research includes workforce scheduling problems that are particularly useful in reality and affect many organizations. Indeed, the workforce cost constitutes one of the highest costs in any organization. The major problems addressed in the first stream are technician routing and scheduling problem (TRSP) [8], technician and task scheduling problem (TTSP) [7], service technician routing and scheduling problem (STRSP) [6], and workforce scheduling and geographically distributed asset maintenance problems (GDMP) [2]. Cordeau et al. [7] proposed a mixed-integer linear program and a solution approach for the ROADEF 2007 challenge. This challenge tackled the specific features of the TTSP in a large telecommunications company. A constructive heuristic was proposed to generate a feasible solution by defining teams and assigning tasks. An adaptive large neighborhood search metaheuristic is then used to solve the problem in the improvement phase. Each technician is specialized in different tasks with different skill levels. He is able to execute tasks requiring lower levels than his as well. Skills requirements and a time window characterize each maintenance task. An outsourcing budget, as well as task priorities, are considered. Kovacs et al. [6] introduced the service technician routing and scheduling problem. Their model minimized routing and outsourcing costs by considering skills and team-building constraints. They used an adaptive large neighborhood search algorithm to solve the problem. They also proposed to select destroy-repair operator pairs and adapted the adaptive mechanism consequently. Pillac et al. [8] dealt with the dynamic technicians routing and scheduling problem (DTRSP). The assignment of technicians with adequate skills, spare parts, and tools to tasks is realized to minimize the overall cost. To handle more requests, technicians can replenish tools and spares at the depot at any time. The authors proposed a re-optimization approach for the periodic problem. This approach relies on a parallel adaptive neighborhood search (RpALNS) algorithm, which generates a new route each time a request arrives. The suggested parallelization scheme distributes the computational effort among the different processors. Their metaheuristic has achieved the same performances as the state-of-the-art results for the dynamic vehicle routing problem with time windows. Çakırgil et al. [9] studied the multi-skill workforce scheduling and routing problem for an electricity company. Different skills requirements and priorities characterize the maintenance operations, and teams of technicians need to be formed to perform them. A mixed-integer programming model was proposed to complete higher priority tasks earlier and minimize total costs. The authors proposed a two-stage matheuristic based on the variable neighborhood search to solve the bi-objective problem.
The second stream of research is concerned with the combination of maintenance and routing characteristics. Workforce scheduling is not the only problem dealt with in this case since other maintenance problems are taken into account. Reliability analysis can be included to design solutions that assign technicians at the right time to perform maintenance operations. Workforce costs and maintenance costs are both dealt with in this case. Lopez-Santana et al. [4] proposed a mathematical model called the combined maintenance and routing (CMR) and a two-phase procedure to solve it. In the first phase, a maintenance model is solved to determine the optimal times to perform preventive maintenance operations, their frequency, and their time windows while minimizing the total expected maintenance cost. The output data of the maintenance model is then considered in a second phase as the input data of the routing model that schedules maintenance operations for each technician. The maintenance model is again solved using the updated start times of preventive operations obtained by the routing model to connect the two problems. This procedure is repeated until meeting the stopping criterion. The nonlinear and convex maintenance cost of the objective function is approximated using a piecewise linear function, and the problem is solved for small instances using a commercial solver. Jbili et al. [3] modeled an integrated strategy of vehicle routing and maintenance. They considered vehicles used in transcontinental transportation that are subject to random failures on the road. These failures have to be repaired, which may take random durations that induce delays. A policy consisting of replacing the critical component when arriving at selected customers is adopted. A mathematical model is proposed to determine the optimal routing and sequence of PM actions simultaneously. The objective is to minimize the total expected cost per unit time. It considers the reliability of the vehicle, maintenance (PM operations and minimal repairs), transportation costs, and finally, maintenance durations and penalties incurred by late arrivals. A genetic algorithm is then proposed to solve large instances. Chen et al. [2] studied the maintenance of gully pots or storm drains. They modeled a multi-period VRP, which considers the risk impact of gully pot failure and its failure behavior each day. The risk impact is estimated using meteorological information. A risk-driven analysis is adopted to evaluate maintenance actions. They focused on two factors: parked cars and gully pots status information. The latter has been proven to be the dominant factor that may negatively affect the scheduling of maintenance actions. Rashidnejad et al. [5] presented a bi-objective multi-period model of preventive maintenance planning in geographically dispersed systems through prognostic information and remaining useful life (RUL) called the integrated vehicle routing problem with time windows and maintenance scheduling (IVRPTW-MS). The first objective minimizes the total cost composed of the performing maintenance cost, the expected failure cost, and the travel cost, whereas the second objective minimizes the unavailability of the assets. The authors used fixed costs and did not include corrective maintenance. They used a non-dominated sorting genetic algorithm II (NSGA-II) to solve this NP-hard problem.
There is a relatively small amount of research combining preventive and corrective maintenance strategies since most of the papers examined dealt only with preventive scheduling, apart from [2,4]. Moreover, they did not consider the uncertainty of the breakdowns and did not deal with a multi-objective perspective. Papers that considered random breakdowns among those examined are [3,4]. The multi-objective formulation is considered in [5,9]. Analyzing the previous problems emphasizes the need to propose solutions to the real problem that includes the following realistic features: large scale, uncertainty, the combination of both corrective and preventive maintenance, and finally, routing of technicians in a multi-objective perspective. To the best of our knowledge, no previous study investigated these features together. These aspects make the NP-hard optimization problem even more challenging. Therefore, heuristic-based approaches are necessary for large instances since exact methods are practically limited. Due to their success in solving many mono-objective combinatorial problems, VND and GVNS algorithms were extended in several studies to deal with multi-objective optimization problems. However, most of these extensions do not question the utility of using properties working well for evolutionary algorithms such as dealing with a population of solutions, aggregation, etc. They also do not propose a specific local search strategy for multi-objective optimization apart from Paquette et al. [10]. Instead, most literature algorithms apply mono-objective local searches by dealing with each objective separately [11] or use previously proposed local search strategies [12,13] such as Pareto local search [13].
The present paper comes with the novelty of exploring the combination of vehicle routing and maintenance problems in a unique bi-objective model and proposes novel VND and GVNS algorithms to solve it. This paper proposes a bi-objective model for the joint maintenance and routing problem. The first objective minimizes the total travel cost and the penalty cost. The second objective can be minimizing the total preventive and corrective maintenance cost. This cost was proposed by Lopez-Santana et al. [4] in their single objective model. The paper also introduces a nonlinear and uncertain failure cost that uses information from equipment degradation, with the routing model as a second objective. The maintenance hypothesis of renewal theory is considered on both costs through a continuous time for the last restoration. Penalties for late arrivals are also associated with the second considered objective. Moreover, the failure and maintenance costs adopted have not been used in a multi-objective study with the routing objective. The workforce cost and the maintenance cost represent the highest costs in a plant. Failing to optimize these costs together can lead to a serious loss for the manufacturing company and reduce its profitability. The maintenance and transport are support processes often outsourced due to their importance. They are directly linked for a service provider to the production of its service. This model simultaneously optimizes these costs in an organization, making it appropriate for real-world situations. Maintenance service providers can use this model to provide the best services to their customers whenever the maintenance strategy adopted by those customers is time-based maintenance, including preventive and corrective maintenance. The previous research gaps of the problem in the literature are therefore filled. To solve the integrated maintenance and routing problem efficiently, multi-objective new adaptations of the variable neighborhood descent and the general variable neighborhood search algorithms are proposed. This paper describes the design of the improvement method, the new best improvement strategy proposed MOBI/P, the acceptance criterion, the stopping criterion, and the approach adopted to reduce the computational time. The paper finally presents an analysis to demonstrate the efficiency and novelty of the proposed mechanisms compared to the literature. This work differs from the literature since it proposes novel adaptations of VND and GVNS algorithms to solve multi-objective problems. Indeed, it does not use the existing local search strategies but a new one called MOBI/P. Moreover, the VND algorithm, which is to be used as an intensification algorithm, is designed to be far faster than the GVNS algorithm. In addition, these algorithms include more design features to lead to better solutions than the literature.

Problem Definition and Formulation
We consider a set M of machines that are located in dispersed customers' sites. They are subjected to random breakdowns due to the failure of a critical component. For each machine and at regular time intervals, preventive maintenance (PM) interventions are scheduled. A PM operation has a cost Cpm i and lasts a duration T pm i . Each preventive operation must be performed within its corresponding time window. However, if it is not possible due to the high workload of technicians, late arrivals are then allowed. Hence, technicians can arrive at an operation i following the end of its time window b i . A penalty cost p ik is incurred in this case. On the other hand, a team of technicians has to wait until the earliest time a i to start an operation i. Teams of technicians K are available to perform both preventive and corrective maintenance. They perform corrective maintenance (CM) operations when machines unexpectedly break down. In this case, the machine stays in a failure state until the arrival of the teams of technicians for a certain time W i . The unit waiting cost for each operation is Cw i . A CM operation has a cost Ccm i and lasts a duration Tcm i . A minimized number of vehicles among the available ones must be determined and used for all operations to reduce transport costs. The Node 0 is considered the departure point of maintenance teams, and the Node n + 1 is considered the final destination. We denote by O the set of all nodes of the PM operations related to all machines, and by V, the set of all locations (operations nodes plus the depots). Each team of technicians has to perform the maintenance operations at the customers' sites. Thus, the model could be defined as a directed complete graph G = (V, A), where A = {(i, j), i ∈ V o , j ∈ V d , i = j} with V o and V d as the sets of origin and destination vertices. The following assumptions are considered: • All maintenance teams have the same skills and qualifications to do the maintenance operations. • It is considered that the random variable of the time to failure for each machine follows a Weibull distribution whose shape and scale parameters are, respectively, β m and σ m , m ∈ M. Any other distribution resulting from the historical data of machine failures can be applied. • We suppose β m > 1, m ∈ M to deal with the third part of the bathtub curve. This part features the wear-out life of the machine when the failure rate is an increasing function of age or usage. Indeed, the wear-out is the phenomenon accelerating the risk of failure over time. Each PM operation must be performed in its associated time window. However, if it is not possible, we allow a team of technicians to arrive following the end of its time window. In this case, a penalty cost is added to the total cost. • The failure of a machine is generally due to the failure of its critical components. • Failures of individual machines are statistically independent. • A machine can have several maintenance operations over the planning horizon. • In case of failure, before the beginning of the next PM operation, the customer must wait for the team of technicians to perform a CM operation instead. The team is not rescheduled based on this new situation. • The travel times are deterministic and satisfy the triangle inequality.

Notation of the Joint Maintenance Scheduling and Workforce Routing Problem
The notations below are employed throughout the paper.

The Description of the General Approach
The maintenance model is solved to determine the optimal time φ i for each maintenance operation i that minimizes the total maintenance cost. The optimal date φ i to perform a PM operation i is used to determine the durations d i for the PM operations, the total number of operations n and a time window interval [a i , b i ] that minimizes the total maintenance cost. The above steps have been adopted by Lopez-Santana et al. [4]. We then solve the defined joint maintenance and routing model that minimizes the routing and the maintenance costs or the routing and the failure costs. Penalties for late arrival are considered in both cases. In the first case, the maintenance cost is minimized considering the workforce routing constraints that incorporate maintenance requirements (maintenance time windows, etc.). Integrating technical maintenance requirements with transport management is the first merit of the approach. In the second case, the failure cost is minimized to reduce the risk of failures. However, the optimal time chosen must be within a time window that minimizes the total maintenance cost. This latter case considers three objective requirements simultaneously while the minimization of maintenance cost is included in time windows constraints. Expressing some objectives of optimization problems as constraints reduces the problem complexity as adopted in [14]. The second merit of our approach is integrating several maintenance strategies. Indeed, the chosen optimal time reduces the risk of failures, the maintenance cost, and the travel cost (risk-based maintenance, etc.). The outputs of the problem are shown in Figure 1. The workforce routing constraints influence the time to perform maintenance operations, the time of the last restoration that depends on that latter in the combined model, and the waiting time. All the variables are interconnected. Therefore, it is essential to model an integrated problem and to solve it using multi-objective optimization.

The Maintenance Model
The decision model to determine the optimal interval between PM operations is referred to as the maintenance model. It is used when the aim of performing maintenance activities is to minimize the total related maintenance costs of preventive and corrective maintenance operations. The optimal period δ m to perform a PM operation is determined for each machine m with the frequency n m of PM operations in the planning horizon.
In the following, we define the main terms employed in the maintenance decision model. The probability of failure in the interval [0, δ m ] of machine m, where P m is the probability's notation, is: The model generally applied by researchers in the literature to settle optimal times to carry out PM operations in the case of periodic preventive maintenance is defined in [15]. Given a machine m, we seek the optimal time δ * m to execute PM operations that minimizes the total maintenance cost CM m (δ m ): The terms E[CM m (δ m )] and E[T m (δ m )] refer, respectively, to the total expected cost of a cycle and the expected cycle length for a machine m.
Lopez-Santana et al. [4] propose an extended maintenance cost which includes the waiting cost expression in addition to the PM and CM durations. These service times cannot be negligible when the systems under study are large-scale.
M m (δ m ) represents the expected time to failure of machine m assuming the failure happens before δ m .
According to Lopez-Santana [4], the waiting time W i of an operation i undertaken by the teams of technicians k is the difference between the technicians arrival time to the PM task θ i,k and the mean time to failure M m (θ i,k ): The values of θ i,k are needed to compute the waiting times. They are also the outputs of the routing model. Consequently, the number of W i obtained is equal to the number of PM operations to be performed for the machine m in the planning horizon. Since the maintenance model takes only one value, Lopez-Santana et al. [4] consider the average value of all W i , where i refers to the PM operations related to the machine m to update the waiting time W m and iterate the procedure. However, the values of θ i,k are unknown before solving the routing model. Therefore, we suppose a PM operation will be scheduled at δ m . This approximation has also been used by [16].
The waiting time is therefore defined as follows: The maintenance cost per time unit for a machine m at the time δ m is finally equal to: This nonlinear equation without constraints is solved to obtain the optimal period for each machine m, δ * m = argmin CM m (δ m ) . The frequency of a PM operation on the planning horizon H can be obtained as follows: The frequency of a PM task is the number of times it is realized in the horizon. Considering a frequency n m of a machine m, the PM operations need to be performed The execution date φ mo of an operation o corresponding to machine m is equal to φ i which represents the execution date of an operation i and is obtained as follows: The time window [a i , b i ] of PM tasks associated with a machine are obtained using [a m , b m ] the time window of the PM period of machine m on the first cycle. [a m , b m ] is determined using the percentage of time of postponing or preempting a PM task that we can tolerate. The cost stays relatively low and close to the minimal value CM m (δ * m ) in this interval. [a i , b i ] can be expressed as follows:

The Joint Maintenance Scheduling and Workforce Routing Model
As mentioned previously, a machine has to undergo a number n m of operations within the time horizon. They are indexed by i ∈ {1, . . . , n m }. In this case, the parameters related to operations i related to the same machine m are equal to the parameters of this machine m. F i , σ i and β i are, respectively, equal to F m , σ m and β m . This way, the values of the operations parameters in the routing model are obtained from the maintenance parameters. The mathematical model of the integrated maintenance scheduling and routing problem can be formulated as follows.
The first objective function f 1 regards the transport of technicians. The second objective function f l deals with maintenance. It equals either f 2 or f 3 : The objective function f 1 (23) minimizes the total travel cost related to technicians routing as well as the penalty cost of not respecting the operations time windows upper bounds. The objective function f 2 (24) minimizes the failure cost, the waiting cost, and the penalty cost. Minimizing the machines probability of failure is equivalent to maximizing machines reliability. The objective function f 3 (25) minimizes the total expected maintenance cost and the penalty cost incurred for arriving after the deadline b i associated to PM operations. The penalty term is intended to ensure that the time windows are respected. Therefore, it is considerably small compared to the first term of the objectives whether it is the routing, failure, or maintenance cost and does not increase the correlation between the objectives. The constraints (15) and (16) indicate that each PM operation has to be executed exactly once by one team of technicians. Constraints (17) ensure that the entry of a team of technicians to node i is mandatorily followed by their leaving. Constraints (18) make sub-tours impossible. The purpose of constraints (19) is ensuring that each PM operation is carried out within its time window. Note that we assume a 0 = 0 and b n+1 represents the maximum time of arrival to the depot. It is permitted to arrive after the deadline b i . A penalty p i is then calculated. Constraints (20) determine the number of vehicles needed and that minimizes the total costs. Constraints (21) ensure that there are different teams of technicians in the different tours. They also ensure that every tour starts at the depot. Finally, the constraints (22) impose domain conditions on the variables.

The Novelty of the Proposed Model
Several multi-objective or bi-objective models have been proposed for the maintenance scheduling and workforce routing problem. However, adopting the failure and maintenance costs as the second objective with the routing cost has never been proposed to date in the literature to the best of our knowledge. The second objective function of the model can either be minimizing the total preventive and corrective maintenance cost and the penalty cost (25) or minimizing the failure cost with the penalty cost (24). In both costs related to maintenance, the time of the last restoration previous to the i-th PM operation performed by the team of technicians k, ρ i,k is used to verify the hypothesis of the renewal theory. Including this time in the objectives functions is a novel aspect of the model. The definition of the failure cost itself and its possible use is another novelty of the proposed model.

The Failure Cost:
The second objective function f 2 minimizes the failure cost and the penalty cost. The aim is to maximize machines reliability. The uncertainty aspect is integrated into the failure cost function by incorporating the probability of failure and the reliability of each machine. In the second objective function, we have used the probability of failure, F i (θ i,k − ρ i,k ), for each operation i ∈ O to utilize information from equipment degradation and hence consider failure risks in technicians' assignment to tasks. It includes direct costs (failure cost) and indirect costs (production losses) illustrated by the waiting time. During that time, the failed machines were out of order and were not used by the organization for production activities. Losses are, therefore, supported by the company.
When we multiply it by the waiting cost per unit time, we obtain a cost that can be interpreted as the production loss incurred by the organization due to this machine failure. The failure cost is beneficial in industries where breakdowns are hazardous and influence safety. This term is nonlinear and includes a random variable. It is equal, in the case of the Weibull distribution, to: The Maintenance Cost: The objective function f 3 minimizes the total expected maintenance cost with the penalty cost. The maintenance cost balances both preventive and corrective maintenance to find the least cost considering the two strategies. The Task Duration: The duration of the maintenance operation i is probabilistic. It is simplified in the model as proposed by [4] to be deterministic. Its real expression proposed by [4] is:

Proposed Multi-Objective Algorithms
The proposed mathematical model is a mixed integer nonlinear program (MINLP) which cannot be solved by commercial solvers in a reasonable time when applied to largescale instances. VRPTW is NP-hard [5], and the classical maintenance scheduling problem using operations research techniques is NP-hard as well [5]. Naturally, the combined problem is an NP-hard problem.
The following section presents an adaptation of variable neighborhood descent (VND) and general variable neighborhood search (GVNS) to the multi-objective context to solve the proposed bi-objective combined maintenance and routing problem. At first, the method adopted is explained, and multi-objective notions used in this paper are defined; then the solution representation is shown. Next, the functions used by the algorithms and neighborhood structures are presented. Finally, the proposed multi-objective algorithms are described, as well as their novelty compared to the literature. The proposed algorithms extend single objective variable neighborhood descent (VND) and general variable neighborhood search (GVNS) proposed by Mladenovic and Hansen [17] to solve multiobjective problems. The algorithm MOVND/P is designed to be an intensification local search component of the GVNS algorithms, whereas MOVND/PI and MOGVNS/P are intended to solve the multi-objective problem. They use the Pareto dominance strategy and start with a unique initial solution. They also incorporate novel proposed strategies. Additionally, they all use a novel multi-objective best improvement strategy called MOBI/P. The other GVNS variants aim to measure the impact of the inclusion of a decomposition strategy on the algorithms. MOGVNS/D uses a weighted sum method instead of the Pareto dominance concept and a population of solutions, whereas MOGVNS/CDP tests the use of a population of solutions with MOGVNS/P.

Pareto Optimality
Pareto Optimal Dominance: "The solution s dominates another solution s " is denoted by

Solution Representation and Constraints Handling
A solution is represented as a set of K tours, where each tour is a permutation of PM operations. All the constraints considered are hard apart from verifying the upper bounds of the time windows. We allow the arrival of a team to an operation after the latest permitted time defined by its time window. We therefore accept only feasible solutions. Figure 2 shows the solution representation that was used in our implementation and the results for one instance when minimizing only the first objective (the routing cost). The solution that minimizes this cost is composed of two routes (two vehicle are therefore needed). Each route is composed of three tasks in the specified order. Section 3.4 explicitly indicates how the three objective costs are evaluated and how the start time is calculated θ ik , i ∈ O, k ∈ K. For instance, five machines are being considered. Each machine has one operation, apart from machine three that has two operations: three and four. The renewal time ρ ik for all operations apart from four is therefore zero. The renewal time ρ 42 for task four related to machine three is in route two. It is equal to the start time of task three associated with the same machine, and that precedes task four. The data of this specified instance is detailed in [18].

Initial Population Generation
We use a weighted aggregation of the two objectives when generating an initial solution. The weighted sum method is used in MOGVNS/D to evaluate and compare solutions. In contrast, the other proposed algorithms (MOVND/P, MOVND/PI, MOGVNS/P, and MOGVNS/CDP) use a multi-objective evaluation and can therefore start with any weightless solution. The weighted sum method is a linear combination of weights. The evaluated function using this method is: The algorithms start with initial solutions constructed using the best insertion heuristic. This latter calculates the minimum insertion cost for each operation that has to be inserted in the solution. The insertion cost of an operation i represents the difference between the cost solution with and without operation i. At each iteration of the heuristic, the operation with the minimum insertion cost is inserted at its best position in the current solution. The process stops when all the operations are inserted. The proposed algorithms, MOGVNS/CDP and MOGVNS/D deal with a population of solutions. Each population individual is generated using the best insertion heuristic and specific weighted sum method. The detailed procedure is described in Algorithm 1. In decomposition approaches, each population individual is associated with a subproblem. It is necessary to decompose the multi-objective problem effectively to reach the maximum parts of the Pareto front. Subproblems are defined in our case by firstly generating the weights. These weights (or search directions) are chosen to explore different directions. The best insertion heuristic is then used to construct each population's individual using these weights.

Local Search Procedure and Neighborhood Structures
Our VNS and VND algorithms rely on the set of the neighborhood structures used and particularly on the sequence of their execution. The neighborhood structures are classified and then applied from the best to the least performing [17]. This work adopts the following sequence of neighborhood structures: the swap, the insert, the 2-opt*, and the 2-opt operators. We use approximately the same order as the literature for VRP problems [19,20]. The results of our previous work [18] show that with well-chosen operators, the GVNS algorithm is an effective method to solve the single objective version of our problem. Its general structure can be applied to improve other methods performance [21]. Preliminary tests were realized using a steepest descent heuristic (best improvement local search) to verify this order. The semi-randomization and the randomization of the operators have been tested in both shaking and VND phases. They worsen the solutions obtained. In this study, we have used the classical BA strategy for the mono-objective GVNS and a novel proposed BA strategy, the MOBI/P, for the multi-objective algorithms to obtain a complete Pareto front. In our Pareto-based algorithms, a solution s is considered better than an other solution s if it dominates it or if it is incomparable to it. In other words, if after applying a move operator f (s ) ≺ f (s) or f (s ) ⊀ f (s) and f (s) ⊀ f (s ), the solution s is considered as a new efficient solution and the set A is updated by means of the addSolution method. The intra-route moves are used on a single route, while inter-route moves intervene on multiple routes. We examine in the neighborhood exploration all possible positions for each operator. The three first neighborhood structures are inter-routes, and the 2-opt is an intra-route operator. The insert and 2-opt* procedure may modify the number of operations in each route as well as their order. Here's a description of the different possible moves: The swap move exchanges two operations either in the same route or between different routes. The insert move deletes an operation from its position and inserts it in another position that can be in the same route or in a different route. The 2-opt* operator generates a neighboring solution by removing arcs (i, i + 1) and (j, j + 1) belonging to two distinct routes and reconnecting arcs (i, j + 1) and (j, i + 1). It withdraws two edges from a route and replaces them with two other edges to form new routes. This operator is therefore inter-route. The 2-opt operator removes arcs (i, i + 1) and (j, j + 1) from a route and links arcs (i, j) and (i + 1, j + 1) in the same route. This operator is the same as a reverse operator that reverses the elements between i and j + 1. This operator is intra-route since it intervenes on arcs belonging to the same route.

Updating Non-Dominated Solutions Set
The addSolution method described in Algorithm 2 is used to update the set of nondominated solutions in the archive A [12,22]. This method uses Pareto dominance to evaluate solutions.  The multi-objective variable neighborhood descent based on Pareto dominance is an adaptation of the variable neighborhood descent algorithm (VND) to solve our multiobjective problems. It is called MOVND/P. This algorithm has been designed primarily as a multi-objective local search component and an intensification algorithm. Its goal is to improve the performance of other algorithms. This algorithm is based on the general principle of exploring several neighborhoods when there is still an improvement and is inspired by the single objective VND.
Algorithm 3 describes the proposed MOVND/P algorithm. At each iteration and while the archive is still improving (steps 6-28), the current neighborhood of s is entirely explored with a novel multi-objective best improvement strategy called MOBI/P. This MOBI/P procedure tests if each neighbor of s is non-dominated with the best solution found so far in the neighborhood of s. This best solution changes during the search. The MOBI/P strategy is very impactful since it allows the search to be more efficient and diversified and enables the algorithm to converge rapidly. All non-dominated solutions by the best solution found so far in the neighborhood of s are then stored in the set P (step 11).
Each point in P generated with the MOBI/P procedure that is not in the archive A is evaluated to enter this latter or not (steps [13][14][15][16][17][18][19]. The solution x is included in the archive A if it is non-dominated by s and replaces the current solution s. This replacement of the current solution s by a solution that dominates it or that is incomparable to it is another new feature used. It can be noticed in line 15 and is directly inspired by a single-objective local search.
The counterAISecond measures if some solutions have been added to the archive from the previous iteration to the next iteration. It is incremented in line 17. The aim is to continue running the algorithm while there is still an improvement. The counterAISecond precisely measures if there is an improvement compared to the counter of the previous iteration counterAISecondPrevIt. There is, however, a stopping condition on the number of iterations iter max to avoid large computational time.
We define a new multi-objective neighborhood change procedure. In our case, we say that a neighborhood N l improves the archive A if all the solutions returned in P are nondominated by the solutions in A. An improvement means that all new points have been added to the archive A. This is recorded using the counter counterAI, which is incremented at each improvement.
A new characteristic was also incorporated in the neighborhood change procedure. In our algorithm, we stay in the improving neighborhood if all the solutions in P are nondominated by s but also when counter did not reach iterC max . The last part of the condition is required to avoid not exploring the other neighborhoods when the first neighborhood is constantly improving. We can obtain an efficient mix between staying in the best neighborhood and exploring the different neighborhoods. The multi-objective variable neighborhood descent improved based on Pareto dominance (MOVND/PI) presented in Algorithm 4 is an improvement of our Algorithm 3 presented in Section 4.6. The improved algorithm can be used to solve multi-objective problems but is not intended to be a local search component such as the previous algorithm. We allow it, therefore, to consume more computational time through a more sophisticated selection criterion. The new solution s for the following iteration is randomly chosen from the set of non-dominated solutions A among solutions that have not been previously explored in the search. The aim of this criterion is to avoid the intensification of a point already visited and exploited. The algorithm, therefore, does not include the line 15 of the Algorithm 3 to select a solution for the next iteration. The only criteria here to continue running the algorithm is the maximum number of iterations. Algorithm 5 describes the GVNS-based algorithm proposed. At each iteration of MOGVNS/P and while the archive is still improving (steps 7-33), a shaking procedure is applied on the current solution s to obtain s . The shaking procedure selects a random solution s from the current kth neighborhood of the current solution s (s ∈ N k (s)). The aim of this phase is to diversify the search process. This perturbation is applied nS times.
An intensification step is then applied to the perturbed solution s using the MOVND/P described in Section 4.6. The non-dominated solutions obtained by MOVND/P are stored and returned in the set P V ND (step 16). All the solutions of the P V ND set are compared to s to update the archive A using the function addSolution. The improvement of the archive A while exploring P V ND is recorded using counterAI and counterAISecond from one iteration to the other such as in the MOVND/P algorithm.
We stay in the same neighborhood if all solutions in P V ND are added to the archive (steps 25-29). The multi-objective neighborhood change procedure is the same as defined in MOVND/P and MOVND/PI.
In this algorithm, and differently from MOVND/P, the non-dominated solutions in the archive A are used to restart the current solution s. Indeed, the new solution s for the next iteration is randomly chosen from the archive A among solutions that have not been already explored.
The step-by-step procedure of MOGVNS/P can be summarized as follows. First, an initial solution is constructed using the best insertion heuristic and is added to the archive of efficient solutions A. The main procedure is then repeated while there is still an improvement of the archive A and while not reaching a maximum number of iterations iter max . It is composed of the following steps: the exploration of the neighborhood struc-tures while k ≤ k max and a random selection in the archive A for the next iteration of the algorithm. The neighborhood exploration phase is composed of four steps. It starts with a shaking procedure to perturb the current solution, followed by an intensification phase using the MOVND/P algorithm to generate the Pareto front P V ND . Next, a test is applied to verify that the points in P V ND are non-dominated by s to enter the set of efficient solutions A or not. Finally, the decision to stay in the improving neighborhood or explore other neighborhoods is made. Counters are used throughout the algorithms to define stopping criteria, including counters that record the improvement of the archive A.  The single objective VND applies a weighted sum function to evaluate and compare solutions. Given a solution s, the weighted sum function associated to a solution s and weight vector (w 1 , w l ) is calculated as follows: g(s, w 1 , w l ) = w 1 × f 1 (s) + w l × f l (s) with l = 2, 3. The initial population is generated using the procedure described in Algorithm 1. Compared to the previous ones based on Pareto dominance, the main drawback of this method is that the Pareto front returned is an approximation of the supported efficient solutions. It can also be time-consuming since we start from a population of solutions. On the other hand, this algorithm is easily convertible to the single objective GVNS algorithm since it is based on aggregation.

Novelty among Acceptance Strategies in Multi-Objective Optimization
Several authors have proposed and used several adaptations for the first-accept (FA) and best-accept (BA) in the multi-objective optimization. Therefore, we analyze and review the existing ones and the proposed method hereafter.

First-Accept based on Pareto dominance:
This procedure consists of accepting the first solution s that is non-dominated by s. It has been used by Cota et al. [12]. Best-Accept based on Pareto dominance PLS: Paquete et al. [10] proposed Pareto local search to apply the BA strategy in the multi-objective context. All the neighborhood of s is explored, and the new solution is accepted if it is not dominated by any solution in the set of efficient solutions. Best-Accept based on Pareto dominance MOBI/P: We use in this paper a novel best-accept strategy for multi-objective optimization called Multi-objective Best Improvement strategy based on Pareto dominance (MOBI/P). This procedure tests if each neighbor of s is nondominated with the best solution found so far in the neighborhood of s. This best solution changes during the search. If the new solution s dominates this best solution, it replaces it, and so on until exploring all the neighborhoods and retaining a unique last best solution that is the most converging. Choosing the last best solution permits a translation in the Pareto front, ensuring diversification if the new solution s is incomparable with the best solution or more convergence if it dominates it. Moreover, it is faster compared to the PLS strategy since we avoid comparing each solution in the neighborhood to each solution in the archive as in PLS. Best-Accept for each objective separately: This procedure consists of applying a single objective local search but for each objective separately. Duarte et al. [11] use VND-1 (respectively, VND-2) to improve the value of f1 (respectively, f2) regardless of the value of f2 (respectively, f1) in a bi-objective problem. The final local optimum is then tested to enter the archive or not. Combination between the FA and BA strategy in PLS: Dubois-Lacoste et al. [13] propose an alternative to accept only dominating solutions to speed up the search. They suggest switching to the criteria adopted in PLS of accepting non-dominated solutions if such solutions are no longer found. The authors also proposed an alternative strategy. It consists of using the first-improvement technique to converge first to a good approximation of the Pareto front until all solutions in the archive are flagged as visited. Afterward, one can move to the PLS best improvement strategy to complete the archive with the remaining neighbors.

Novelty in the Design of the Multi-Objective Algorithms
We consider an improvement if all the solutions explored have been added to the archive. This procedure is different than the definition of improvement proposed by [11], where improving means at least one new point has been added to the archive A.
Likewise, we stay in the improving neighborhood if all the solutions explored are non-dominated by s but also when a defined counter does not reach a maximum number of iterations. This last condition is required to avoid not exploring the other neighborhoods when the first neighborhood is constantly improving. In the literature, if there is an improvement, the exploration continues in the same neighborhood structure [11]. Queiroz and Mundim [23] propose to change the neighborhood systematically at each iteration in an adapted template from [11] to reduce the computational time.
The proposed algorithms use the MOBI/P procedure above as a multi-objective best improvement strategy, whereas there are several different FA and BA accept strategies in the literature.

Computational Experiments
This section presents computational experiments carried out to measure the performance of the proposed algorithms. We present the results of the mono-objective problem using the GVNS algorithm and the results of the multi-objective problem using the proposed algorithms MOVND/P, MOVND/PI, and MOGVNS/P. The results of the other proposed variants, MOGVNS/CDP and MOGVNS/D, are also shown to demonstrate the performance of the three aforementioned algorithms. In addition, CPLEX results are reported for the mono-objective problem. Heuristic Pareto fronts are also provided for comparison to our multi-objective problem. The maintenance model without constraints was solved using Python 3.6. The joint maintenance scheduling and workforce routing model was then solved using C++. This work has been compiled with GCC 7.4 in a Linux environment. To solve the MILP using an exact method, the concert technology library of CPLEX 12.10.0 version with default settings in C++ has been used. Experiments were conducted using the CALCULCO computing platform in an AMD EPYC 7702 with 2CPU, 2 gigahertz, and 1 core was dedicated to each instance. All methods are run using the same machine to avoid bias.

Instances Description
We use two sets of instances. The first one, denoted by Re, is inspired by industrial reality. In the Re instances, large and short maintenance durations are considered (Tcm and T pm go from 0.5 to 48 h). We set the shape parameter of the Weibull distribution β > 1 to consider the wear-out period of the machine's life. These machines are, therefore, old. Short (27 km) and long distances (up to 500 km) are considered. The speed is fixed to 60 km/h, as adopted in real-world scenarios. The horizon is set to H = 101 h for n = 12 and H = 80 h for n = 6. The penalty cost is fixed to c = 10. This cost is high enough to produce the desired effect of penalizing the objective functions whenever multiplied by the penalty variables. The depot time window is the interval between the lower bound a 0 = 0 and the upper bound b n+1 . For this class of instances, b n+1 = H. These data are shown in [18]. There are six machines that may need more than one maintenance task in the planning horizon in the class of instances Re. The travel time and distance between the same machine operations are naturally zero. The number of available vehicles is initially fixed to 4 for n = 12 and 2 for n = 6. We note that in some obtained solutions, we can use fewer vehicles than those available. The second set is derived from the well-known Solomon benchmark. These benchmark problems are available on Solomon's web page http://web.cba.neu.edu/ msolomon/problems.htm. We used the first 10, 25, 50 machines of Solomon's classes of instances. The number of operations used in these instances varies from 23 to 91 operations, and the number of vehicles ranges from 5 to 35. This latter can also be reduced during the search.
Solomon's instance are classified in three groups: • R: randomly distributed locations. • C: geographically clustered locations. • RC: partially randomly distributed and partially clustered locations.
There are six Solomon's classes with different locations coordinates: C101, C201, R101, R201, RC101, RC201. The unitary speed is adopted. For each of these six classes, we consider the horizon H = 100 h with a latest arrival time at the depot b n+1 = 200 h and H = 200 h with b n+1 = 400 h. The opening time is a 0 = 0. The maintenance parameters were generated as described in [4]: T pm i ∼ U c [5,10], Tcm i ∼ U c [15,30], [10,20], and σ ∼ t 0i * U c [2,5]. Only the parameter β ∼ U c [2,6] is generated differently compared to the one described in [4] so that it can be closer to reality. The continuous uniform distribution U c has been used to generate the random information. In the instances sets, the percentage of time windows tolerance is set to tol = 7%. It is inspired by the values usually used in large-scale industries that consider a restricted time window. Another significant time windows tolerance is considered tol = 30%.

Parameter Setting and Performance Metrics
To set the algorithm parameters, we have run several preliminary tests. We noticed that a high value of the shaking parameter nS and a high number of iterations have a negative impact on the computational time. The retained parameters for the maximum number of iterations is iter max = max(1, E[n/8]) for the routing objective (23) and the maintenance objective (25), where the symbol E stands for the whole experiments parts of the integer portion of the number. For the failure objective (24), this parameter is iter max = max(1, E[n/4]). The diversification parameter nS in the shaking phase is set to three. We compare the results of the CPLEX solver and our mono-objective GVNS algorithm for the routing objective function (23). In the CPLEX columns, we report the objective value, the CPU time, and whether the solution is optimal, feasible, or not found after 96 h of execution. In the GVNS algorithm, we indicate, respectively, the maximum, minimum, and average values for five runs. We also provide the corresponding maximum, minimum, and CPU time.
The pseudo-code of mono-objective GVNS and VND are shown in [18]. The monoobjective GVNS can directly be obtained with MOGVNS/D when the initial population is reduced to one individual. The gap between the objective value of our metaheuristic and the CPLEX solution is calculated as follows: This gap represents the percent decrease of the cost objective when the GVNS algorithm is used in comparison to the use of the CPLEX solver. It is the case whenever the CPLEX solver only found a feasible penalized solution, and the gap is different from zero.
The percent decrease of CPU time when using the GVNS algorithm compared to the CPLEX solver is given by:

Test Results
The results of the GVNS and CPLEX solver for the routing objective are represented in Table 1. Bold values are the minimum objective values obtained by either the GVNS or the solver. Whenever a number in the improvement column is bold, there is an improvement. The quality of the initial solution obtained with the best insertion heuristic considerably reduces the CPU time since we start from a reduced objective compared to random initial solutions. The impact of good initial solutions on the CPU time is illustrated in [18]. Indeed, providing a well-chosen starting solution positively influences the execution time and helps to improve the final solution quality rapidly. Optimizing the failure cost requires more iterations to reach high-quality solutions than optimizing the routing cost or the maintenance cost.
For small instances (less than 26 operations) of Solomon's class C, both CPLEX and GVNS were able to solve the problem optimally. However, the GVNS algorithm consumed less CPU time than CPLEX for most of the instances. When the number of operations increases to more than 34, some optimums were obtained for class C, R, or RC instances. For these instances, CPLEX returns mostly either feasible solutions or optimal solutions but after a very long CPU time (more than 48 h). In many instances, CPLEX fails to find optimal solutions after one week. In most cases, the GVNS algorithm improves the objective value of feasible solutions obtained by CPLEX. The Gap of GVNS for most instances is 0%. The algorithm is also very robust since the difference between the maximum, minimum, and average values is very small, almost negligible.
For the 47 tested instances, CPLEX returned either feasible or optimal values for 29 instances and failed to obtain any feasible solution for the remaining 18 instances, even after a week of execution. For the 29 instances for which we obtained either feasible or optimal solutions using CPLEX, the average gap or objective improvement of GVNS over CPLEX is 18.47%. The average improvement of CPU is 32.81%. When removing the four Re instances for which the number of operations is inferior to 12, we obtain an average objective improvement in favor of GVNS equal to 21.42 % and an average improvement of the CPU is 77.96%. The gap value is 0 for 17 instances among 29 instances solved. Therefore, the proposed GVNS found the optimums in these instances. The GVNS algorithm improves the CPLEX results on 11 instances out of the 29 instances for which CPLEX returns either feasible or optimal solutions. However, it obtains a slightly worse value for only one instance. The average gap for the 12 instances where the gap is different from 0 is 44.61%. This means that when the gap is different from zero (meaning only a feasible solution is found), GVNS finds a better objective value than CPLEX with a percentage of 44.61%. The average percent decrease of CPU time when using GVNS compared to CPLEX is 85.73% in this case.
Our GVNS algorithm outperforms CPLEX in terms of running time, the longest CPU time being only 3 h (10,906.1 s) compared to the maximum CPU time reached by CPLEX of 95.53 h (343,932 s). In mono-objective optimization, it is simple to evaluate the quality of a solution. It is a more difficult task in multi-objective optimization since the output is represented by sets of trade-off solutions, potentially incomparable in terms of Pareto dominance. Consequently, we use several indicators to measure the quality of an approximation of the Pareto front involving several criteria such as solution quality, computational effort, robustness, and other factors. The indicators assess solution quality and computational effort. The CPU indicator evaluates the time needed for the algorithm to return a final Pareto front. The quality indicators that we used are the number of non-dominated solutions and the hypervolume indicator to measure coverage and assess the front returned convergence.

Numerical Results
The hypervolume (HV) indicator is a metric that measures the size of the space covered. It assesses the space enclosed by all the solutions of the objective space. The HV of an estimated Pareto front is the sum of the hypercubes that each set of solutions contains. It shows the distribution of solutions along the Pareto front. The larger the HV indicator is, the better the Pareto front is.
The reference point ( f The hypervolume is the sum of the areas formulated as follows: where l = 2,3, and p is the number of solutions s i in the Pareto front A. Our bi-objective algorithms based on VND and GVNS frameworks were compared to those of [11]. Algorithms in [11] were run for a longer time to have a complete Pareto front and meaningful results. For a fair comparison, all methods were run on the same machine in the computing platform.
The improvement of our algorithms over the literature algorithm is evaluated using a percent increase of the HV indicator and the number of non-dominated points when using the proposed algorithms PA compared to the comparison algorithms CA of [11]. It is also evaluated using the percent decrease of the CPU time of the proposed algorithms PA over the comparison algorithms CA.

Test Results
Tests are realized over small, medium, and large instances. Small instances include 6 and 12 PM operations. The number of operations varies between 23 and 34 for mediumsized instances and between 52 and 73 for large-sized instances. The indicators assessed are the hypervolume (HV) to measure the convergence and coverage of the solutions in the Pareto front obtained, the number of non-dominated points (NDP), and the CPU time in seconds. Tables A1 and A2 show the results of the proposed variable neighborhood descent algorithms, MOVND/P, MOVND/PI, and the results of the MOVND suggested by [11]. Furthermore, Table A1 details the results for the maintenance cost as a second objective, whereas  Table A2 shows the results when the failure cost is the second objective of the problem.

Comparison with the Literature
Some results that were not displayed in detail are also reported here for the tested instances: the average percent improvement of HV (IHV), NDP (INDP), and CPU (ICPU). They are summarized in Table 2. There is an improvement for bold numbers. When the maintenance cost is the second objective considered, both MOVND/P and MOVND/PI have better average values than the MOVND of [11] on all the considered quality indicators. Indeed, the average percent improvement of MOVND/P compared to MOVND [11] of the hypervolume (IHV) and the number of non-dominated points (INDP) is, respectively, 7.82% and 2.79%. These two indicators improved slightly; however, the CPU time decreased considerably by 81.56%. The second proposed algorithm, MOVND/PI, considerably outperforms the MOVND of the literature [11] for the 46 instances tested. Indeed, the hypervolume and the number of non-dominated points increased, respectively, by an average of I HV = 104.12% and I NDP = 108.45%, which is a considerable improvement. However, the CPU decreased by only 4.73%. Despite the minor improvement for the CPU time, the other coverage and convergence indicators are considerably improved. We also measure this improvement for the proposed VND algorithms over MOVND of the literature [11] for the instances without the four small instances Re. For the 42 derived from Solomon's instances, the improvement of the indicators HV, NDP, and CPU of MOVND/P is, respectively, I HV = 2.32%, I NDP = −3.26%, and ICPU = 88.81%. The improvement of the same indicators for MOVND/PI becomes I HV = 106.38%, I NDP = 105.20%, and ICPU = 28.87%. Table 2. Comparisonbetween the proposed algorithms and the literature algorithms.

Improvement
The Maintenance Objective The Failure Objective

IHV (%) INDP (%) ICPU (%) IHV (%) INDP (%) ICPU (%)
MOVND/P over MOVND of [11]  We consider the failure cost as a second objective with the routing one. The average percent improvement of the indicators of MOVND/P in comparison to MOVND of the literature [11] is considerable, and the values of HV, NDP, and CPU indicators equal, respectively, 85.71%, 67.74%, and 80.79%. The average percent improvement of HV for MOVND/PI compared to MOVND of the literature is 303.78%. The number of nondominated points NDP increased by 304.90%. However, the improvement of the CPU is −10.21%. When the four small instances Re are deleted to consider only instances generated from Solomon's set, the improvement of the indicators stays considerable for both algorithms. The values of HV, NDP, and CPU indicators improved by, respectively, 67.16%, 67.11%, and 88, 81% for MOVND/P and by 324.52%, 325.54%, and 14.01% for MOVND/PI. The quality indicators are improved, and the CPU time of our algorithms is considerably less.
We can conclude that all the indicators were improved for both failure cost and maintenance cost as the second objective. There is a slight improvement for the HV and NDP indicator for MOVND/P for the maintenance cost and a considerable improvement in the CPU time that has decreased. On the other hand, MOVND/PI improves the HV and NDP indicators considerably with only a tiny improvement in the CPU time. The results of the failure cost as a second objective demonstrated a considerable improvement for all the quality and time indicators for both MOVND/P and MOVND/PI algorithms. Table A3 illustrates the results of the comparison between our proposed MOGVNS/P and MOGVNS of the literature suggested by [11]. When the maintenance cost is the second objective, the average percent improvement of HV when using MOVGVNS/P is 52.29%. The number of non-dominated points NDP increased by an average of 52.15%. The computational time CPU of MOVGVNS/P decreased by ICPU = 41.69%. When the failure cost is the second objective, the average percent improvement of HV of MOVGVNS compared to MOGVNS of the literature is 91.47%. The number of non-dominated points NDP increased by 92.11%. The running time is reduced by 47.74%. All the indicators have been improved. We can conclude that the MOGVNS/P algorithm considerably outperforms the literature algorithm in all the assessed indicators for both second objectives.
We fixed a time limit to one week to have all the instances. We have therefore tested the MOGVNS of the literature for fewer instances since it consumes more time.

Comparison between the Proposed Algorithms
Tables A4 and A5 present the results of the three indicators measuring the quality of the solutions obtained and the time to get those solutions for the proposed MOGVNS/P, MOGVNS/CDP, and MOGVNS/D algorithms. The last study aims to compare the performance of the algorithms presented and discuss their difference. When the maintenance cost is the second objective, and for the 46 tested instances, MOGVNS/P outperforms in average MOGVNS/CDP in the three indicators HV, number of NDP, and CPU with, respectively, 2.42%, 1.92%, and 24.73%. MOGVNS/P is better on average than the decomposition-based algorithm MOGVNS/D in the indicators HV, the number of NDP, and CPU with, respectively, 574.41%, 573.14%, and 19.92%. The MOGVNS/P improves MOVND/PI on the indicators HV and NDP by about 25.87% and 26.90%. However, MOGVNS/P consumes more time (−706.20%) than MOVND/PI. In the opposite direction, the MOVND/PI improves MOGVNS/P by −0.38%, −0.23%, and 65.74%. Therefore, we can conclude that MOVND/PI has similar performance for the HV and NDP indicators and is far faster than MOGVNS/P. When the failure cost is the second objective and for the 46 instances tested, MOGVNS/P outperforms on average MOGVNS/CDP in the three indicators HV, the number of NDP, and CPU with, respectively, −2.91%, −1.68%, and 18.24%. The MOGVNS/P is better than MOGVNS/D in the indicators HV, NDP, and CPU with 792.71%, 679.11%, and 13.03%. The MOGVNS/P improves MOVND/PI on the three indicators HV, NDP, and CPU by about 22.02%, 21.75%, and −387.50%. The HV and NDP indicators are improved; however here again, MOGVNS/P consumes more time than MOVND/PI when the failure cost is considered. In the opposite sense, the MOVND/PI improves MOGVNS/P by −7.55%, −7.41%, and 70.87%. We can conclude from this study that MOGVNS/P is better in all indicators compared to MOGVNS/CDP and MOGVNS/D. Moreover, MOVND/PI has approximately the same performance as MOGVNS/P, but it is considerably faster.
The convergence of the proposed algorithms for all the instances tested can be discussed from the HV results of Table 2. The improvement for the HV indicator of MOGVNS/P over MOGVNS/CDP is slight. The two algorithms therefore have similar convergence and coverage performance. They also have approximately the same number of non-dominated points NDP. The mono-objective GVNS, whose results are presented in Section 5.2.2 obtains a solution with the optimal cost found by the solver or a solution with a better objective value than the feasible solution returned by the solver. Therefore, the GVNS algorithm converges. The algorithm MOGVNS/D uses the same components (neighborhood structures, etc.) as the GVNS algorithm and an aggregated function to evaluate solutions. It therefore has the same convergence as the other GVNS algorithms. However, the set of efficient solutions returned includes only supported solutions, which explains its weak performance in the HV and NDP indicators compared to MOGVNS/P. Figure 3 illustrates for a given instance of 52 operations that MOGVNS/D has a good convergence but returned an incomplete Pareto front. MOVND/P is designed to be an intensification algorithm. It is a main component of MOGVNS/P and MOGVNS/CDP. It therefore has less convergence than the latter two. MOGVNS/P outperforms MOVND/PI in the HV and NDP indicators. MOGVNS/P also takes more time. Therefore, the convergence and the coverage of MOGVNS/P are better than those of MOVND/PI. The MOGVNS/P algorithm also returns more non-dominated points than MOVND/PI. This is shown in Figure 3. We can conclude that all the algorithms converge; however, the convergence of the GVNS algorithms is better than the convergence of VND-based algorithms. All the GVNS algorithms proposed in this paper have the same convergence as MOGVNS/P but differ only in terms of coverage.

Discussion of the Results
In conclusion, the proposed algorithms MOVND/P, MOVND/PI, and MOGVNS/P have a better performance compared to the MOVND, and MOGVNS of [11]. One reason can be that the VND used by [11] first tries to improve each objective separately independently from the other using a single objective evaluation. As a result, the algorithms are not totally based on Pareto dominance. The method we used to fully explore the neighborhoods MOBI/P is also very impactful in improving solution quality and reducing computational time. We also force the exploration of all the neighborhoods with specific stopping criteria, while in [11], staying in the same neighborhood is necessary when it still improves. MOGVNS/P has a similar performance as MOGVNS/CDP. However, the latter consumes more time since MOGVNS/CDP is the same as MOGVNS/P applied to a population of solutions instead of one initial solution. The shaking procedure is also a better perturbation than the method we used to generate the initial population and is sufficient to diversify the search. The importance of the shaking procedure also appears when comparing the performance between MOGVNS/P and MOVND/PI. MOGVNS/P outperforms for the quality indicators MOVND/PI. MOGVNS/P is far better than MOGVNS/D for all quality indicators since the latter algorithm returns only supported solutions. All methods have a good convergence towards a good Pareto front. The methods differ in the coverage and convergence indicator, the number of non-dominated solutions, and the computational time. Our best algorithms are MOGVNS/P and MOVND/PI. Figure 3 illustrates the Pareto front obtained on the instance C201/25/18/100/200/0.30. For this specific instance, MOGVNS/P and MOGVNS/CDP have good performance over all other algorithms. MOVND/P and MOVND/PI are better than the MOVND proposed in the literature. Moreover, MOVND/PI is also better than MOGVNS in the literature. MOGVNS/D returns an incomplete Pareto front but has a good convergence. This instance can represent the results of the performance of all the algorithms on average. The proposed MOGVNS/P algorithm has better convergence compared to MOVND/PI. The average percent results obtained in Table 2 show the same performance. However, the results of MOVND/PI, on average, are close to those of MOGVNS/P, although it is less converging.

Conclusions
In this paper, we address a joint maintenance and workforce routing problem. We propose a novel bi-objective mathematical model that aims to minimize both maintenance and transport costs in the case of time-based preventive maintenance. The set of machines is supposed to be geographically distributed and subject to non-deterministic failures. The joint maintenance and routing problem's objective is to simultaneously determine the optimal times to perform preventive maintenance operations on each machine and find the optimal sequence of these operations that simultaneously minimizes the maintenance and routing cost.
There are many contributions to this paper. We first propose a nonlinear stochastic failure cost that uses information from equipment degradation. It includes direct costs (failure cost) and indirect costs (production losses) illustrated by the waiting time. This objective is particularly valuable for industries where failures seriously influence personnel safety and cause environmental damage. We also investigate another maintenance cost previously proposed in the literature that aims to balance both preventive and corrective costs related to maintenance operations. For the first time, a bi-objective approach is proposed to deal with this problem, where we associate either the failure cost or the maintenance cost with the routing cost. Both costs also consider the time of the last restoration. The proposed model considers maintenance operations time windows, penalties related to late arrival, and maintenance costs under uncertainty which are interesting features of real industrial problems. New adaptations of variable neighborhood descent and general variable neighborhood search frameworks called respectively MOVND/P, MOVND/PI, and MOGVNS/P are proposed to deal with combinatorial multi-objective optimization problems, and this problem, particularly. To obtain high-quality solutions, we describe how to design the improvement method, the acceptance criterion, the stopping criterion, and the neighborhood change procedure. These algorithms are based on the Pareto dominance concept and use a new multi-objective best improvement strategy called MOBI/P. We test a pure decomposition approach. The resultant algorithm, MOGVNS/D, decomposes the problem into several scalar subproblems and uses the weighted sum method to evaluate the solutions. This algorithm is easily convertible to a mono-objective general variable neighborhood search. Finally, we test the use of a population of solutions with MOGVNS/P. The resulting variant is called MOVNS/CDP.
Several numerical experiments have been performed to validate our proposals. First, the mathematical model and the proposed algorithms were evaluated on randomly generated instances. For the single objective variant of the problem, the computational results for the linear routing objective indicate that the GVNS significantly outperforms the results of the commercial solver CPLEX in terms of solution quality and CPU time. The obtained results also demonstrate that MOVND/P, MOVND/PI, and MOGVNS/P outperform ex-isting MOVND and MOGVNS in the literature for all indicators measuring convergence and coverage in less computational time. Compared to the other proposed MOGVNS variants, MOGVNS/P is significantly better than MOGVNS/D for all the quality and time indicators since the latter algorithm returns only supported solutions. MOGVNS/P slightly outperforms MOGVNS/CDP in less computational time. Moreover, MOVND/PI and MOGVNS/P have almost similar performances with a time advantage in favor of MOVND/PI. Future research will include exploring the integration of other maintenance strategies with the routing problem to reduce maintenance costs and improve the quality of maintenance services. It will also be promising to test the hybridization of local search and population-based methods to solve the problem.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.