Iterated local search multi-objective methodology for the green vehicle routing problem considering workload equity with a private fleet and a common carrier

Article history: Received May 5 2020 Received in Revised Format July 9 2020 Accepted August 19 2020 Available online August, 19 2020 A multi-objective methodology was proposed for solving the green vehicle routing problem with a private fleet and common carrier considering workload equity. The iterated local search metaheuristic, which is adapted to the solution of the problem with three objectives, was proposed as a solution method. A solution algorithm was divided into three stages. In the first, initial solutions were identified based on the savings heuristic. The second and third act together using the random variable neighbourhood search algorithm, which allows performing an intensification process and perturbance processes, giving the possibility of exploring new regions in the search space, which are proposed within the framework of optimizing the three objectives. According to the previous review of the state of the art, there is little related literature; through discussions with the productive sector, this problem is frequent due to increases in demand in certain seasons or a part of the maintenance vehicle fleet departing from service. The proposed methodology was verified using case studies from the literature, which were adapted to the problem of three objectives, obtaining consistent solutions. Where cases were not reported in the literature, these could be used as a reference in future research. © 2021 by the authors; licensee Growing Science, Canada


Introduction
The solution of the routing problem is of great interest to the productive sectors of a country since the costs associated with transporting goods or services is a significant part of their final value, which implies that it has a direct macroeconomic effect. The solution to this problem implies a positive impact on the economic variables of organizations. Transport companies are aware that the entire operation must include actions that add value to their activities to propose efficient, productive, and sustainable alternatives. According to Garetti et al. (2012), the concept of sustainability is based on three pillars, economic, environmental and social, as shown in fig. 1. In environmental matters, the transport industry generates significant amounts of annually, leading to an increase in environmental problems, such as the contribution to greenhouse gases (GHGs) that cause global warming. Currently, many companies seek to reduce these emissions to reduce their carbon footprint to sustainable levels in ecological terms. A final but not less important aspect is represented in the balance of the workload assigned to the drivers of the vehicle fleet to generate social welfare, which must be considered to comply with labour legislation or to guarantee restrictions specific to the routing problem that real applications demand.

Fig. 1. Pillars of sustainability
The problems related to routing vehicles, known as the vehicle routing problem (VRP), have their beginnings from the travelling salesman problem (TSP), first proposed by Flood (1956), that describes a traveller agent who must visit "n" cities before returning to the beginning point. The entire tour must be done while minimizing the cost. Dantzing et al. (1959) proposed a variation on the TPS problem that considers using "m" vehicles (replacing agents) with capacity restrictions. This problem is known as a capacitated vehicle routing problem (CVRP). Subsequently, new approaches to the problem were generated with all kinds of new variants, such as using time windows, restrictions on the distances travelled, and multiple deposits, among others. In general, in problems that are defined with three objectives to minimize, one of them corresponds to cost. It is not common for problems involving more than two objectives to be studied at the same time due to the complexity of their solution. According to the specialized literature review, the VRPPC is a problem that does not have studies considering multiple objectives, unlike the green vehicle route problem (GVRP) and the vehicle routing problem with route balancing (VRPBR), which are defined as bi-objective problems. The VRPPC is a problem that starts from the definition of the CVRP, but when the capacity of the vehicles is greater than the demand of the customers or there are vehicles in unavailability, using common carrier vehicles, which start their journey in the deposit and end on the last customer route, is considered. When other optimization criteria are considered, such as the impact on the environment due to operating vehicle fleets and workload balance, the problem is defined as a green vehicle route problem with private fleet and common carrier considering workload equity (GVRPPCWE).

Literature review for the VRPPC
The first VRPPC approaches were developed by Ball et al. (1983), for vehicles that must travel long distances to reach customers, so the solution should include the payment and use of subcontracted vehicles. Later Klincewicz et al. (1990), proposed a model that considers the area geographically divided into sectors with random daily demands with a single deposit and decide the best way to address each sector, either with private fleet vehicles or a common carrier. Chu (2005), considered that a fixed number of vehicles with a limited capacity, with a single deposit and customers, have known demands. Bolduc et al. (2007), proposed an improved heuristic called selection, routing and improvement (SRI). After the initial solutions are obtained, improvement is performed, which consists of applying a customer exchange through a procedure called interchange that was proposed by Osman (1993), to obtain the final solution. Bolduc et al. (2008) developed a metaheuristic that implements perturbance to improve exploration; this metaheuristic is called randomized construction improvementperturbation (RIP), which consists of building an initial solution with a random savings algorithm that has two subsequent stages of improvement. Potvin et al. (2011) proposed a metaheuristic based on the tabu search methodology with ejection chains in which the customers are ordered according to a weight to be assigned to the common carrier and routes are created for the remaining customers be attended to by the private fleet through a heuristic based on travel costs. The presented model allows infeasibility in the vehicle capacity in the objective function through a penalty parameter. Huang et al. (2011) proposed a heuristic that consists of relaxing the first constraint of its formulation called set coving problem (SCP), in which it incorporates Lagrange multipliers that must be updated in each iteration. Moon et al. (2012) present a problem considering subcontracted vehicles and drivers overtime called vehicle routing problem with time windows considering overtime and outsourcing vehicles (VRPTWOV). A hybrid genetic algorithm was proposed, which consisted of using the simulated annealing (SA) metaheuristic together with a genetic algorithm (GA), which act collaboratively to find solutions. Euchi et al. (2013) proposed a methodology that uses a hybrid iterative density evolutionary estimation algorithm (IDEA / 2-opt). This metaheuristic is a type of evolutionary algorithm in which an initial population of individuals is improved by applying stochastic operators. This algorithm uses a procedure called pressure selection towards diversity (Bosman, et al., 2002). Huijink et al. (2014) proposed a methodology based on an adaptive search in a variable neighbourhood (AVNS), which can be fast or slow. In both applications, perturbance operators such as cycling move, jump, shifting, create, destroy, split and bomb are used. A tabu search heuristic was used as the local search technique. Toro-Ocamp et al. (2016) proposed a mathematical model for the location and routing problem with private fleet and common carrier (CLRPPC), optimally guaranteeing the routes. The mathematical model presented in this work was carried out based on it. Euchi et al. (2013) proposed a tabu search algorithm for neighbourhoods with ejection chains to solve the VRPPC. This methodology creates a sophisticated neighbourhood structure that allows efficient search in the solution space, as it incorporates the efficiency of tabu search in the intensification and construction of multi-factor neighbourhood movements with ejection chains, which allow various degrees of feasibility generated by powerful exchange movements. Dabia (2019) proposed a new variant for the VRPPC, a cost model associated with the load assigned to the vehicles of the common carriers. This model contemplates the variation in cost due to the amount of cargo and integrates it into the strategic design of medium-term routing. Additionally, time windows were considered. They called this problem a rich vehicle routing problem with private fleet and common carrier (RVRPPC).

Literature review for the GVRP
For vehicle routing problems, the environmental impact related to emissions are considered an important factor in the increase in GHGs. Emissions depend on the type of vehicle used, speed, distance travelled and the load it carries. Palmer (2007) presented a model that considers speed, congestion and time windows, but did not consider vehicle load in their model. Maden et al. (2010) proposed a heuristic to minimize VRP travel time with time windows. This reduction achieved a savings of 7% in emissions for the cases studied. Kara et al. (2007) proposed a new cost function based on the distance and loading of vehicles for the PTRC. The energy cost (E) of vehicle travel from the depot to the customers and back to the depot is defined by a function: , = , + * , , where , is the vehicle load between nodes and , is the empty vehicle weight and , is the distance between nodes and and . The NSGA-II metaheuriscic, which is a genetic algorithm widely used for bi-objective problems, was implemented to solve the problem. Figliozzi (2010) analysed the impact of congestion and customer demand on emissions for vehicle routing with flexible time windows, the vehicle routing problem with soft time windows (VRPSTW). The IRCI algorithm, which consists of sequential construction of the routes followed by applying an improvement stage, was implemented for the solution. Kuo (2010) proposed a model for calculating fuel consumption for the time-dependent vehicle routing problem (TDVRP), considering speed and travel time. The simulated annealing (SA) algorithm was implemented for the solution. Fuel consumption is dependent on the travel time since, according to the departure time, the vehicle can increase or decrease its speed to arrive on time. Suzuki (2011) developed a solution for the travelling salesman problem with load vehicles and time windows (TSPTW). The mathematical approach was taken from that proposed by Wang et al. (2009), using gradients for arcs and an estimate of the load effect. An exhaustive enumeration technique that explores the entire solution space to find the optimal routes was used for instances between 5 and 10 customers. For instances of 15 customers, a metaheuristic called compressed annealing was used. Bektaş et al. (2011) presented the pollution-routing problem (PRP). The GHG emissions rate in grams per second is directly related to the fuel usage rate through the relationship = + , y . The problem was solved using an integer linear programming model. Ubeda (2011), modified the approach of the CVRP, the distance matrix by an environment matrix calculated as = × , where indicates the vehicle load on the arc and is the emissions factor. To solve the problem, a heuristic method was proposed that consists of solving the approximate environment matrix and then re-optimizing the solution with the complete saving matrix. Demir et al. (2012) developed an approach for fuel consumption and emissions based on proposals by Barth et al. (2005), Scora et al. (2006) and Barth et al. (2007). The fuel consumption rate was defined by = + / / , where is the mass ratio between fuel and air, is the friction factor of the engine, is the engine speed, is the engine displacement, and are constant and is the power per second in kilowatts. The problem was solved by metaheuristic adaptive large neighbourhood search. The fuel capacitated vehicle routing problem (FCVRP) was proposed by Xiao et al. (2012) considering the fuel consumption rate (FCR) that depends on the vehicle load. The FCR is formulated as = + + , where corresponds to the weight of the unladen vehicle and Q corresponds to loaded. and are regression coefficients. The simulated annealing algorithm was used to solve the FCVRP. Kucukoglu et al. (2013), in their study of the Green Capacitated Vehicle Routing Problem (GCVRP), proposed calculating the fuel consumption by linear regression as: = + ℎ + . The problem was solved by a model of integer linear programming. Toro et al. (2017) developed a model to calculate O emissions for multi-objective GCVRP, analysing the forces acting on the movement of the vehicle and assuming constant speeds. The proposed formulation was taken as a reference in this research for calculating gas emissions. The calculation is shown in Equation (1): where corresponds to the energy required to perform the displacement in an arc , . is related to the movement of the vehicle without load or empty and the second component, , is related to the movement of the vehicle with the load. Multi-objective GCVRP was solved using an exact model of mixed-integer linear programming.

Literature Review for the VRPEW
According to the literature review by Matl et al. (2017), the vehicle routing problem with equity (VRPEW) is divided into five groups. The first focuses on solving multi-objective problems with route balancing. This problem is known as the vehicle routing problem with route balancing (VRPRB). In the second group, an extension of the first was proposed using time windows. The third group implements a single objective as min-max VRP, where equity is the aim to optimize. The fourth group includes other types of equity implementations for the VRP. The fifth group contains work done as applications to the problem of equity. In this investigation, the techniques referenced in the first and third group will be used to find the balance in the routes. The VRPRB was first presented by Jozefowiez et al. (2002) and treats the problem as bi-objective, developing a parallel and hybrid algorithm for its solution, through a combination of genetic algorithms and tabu search. Based on previous work, Jozefowiez et al. (2005) performed an adaptation for the NGSA-II metaheuristics based on parallel search and improvement in population density management. This work was followed in Jozefowiez et al. (2007), where a solution algorithm called pareto aiming for objective, with a hybrid NGSA-II algorithm, tabu searching and objectives programming techniques, was presented. Lacomme et al. (2015) proposed a multi-objective genetic algorithm that is improved using constructive heuristics for generating the initial population and including a local search process.  presented a population algorithm based on local search to solve the bi-objective problem of minimizing route lengths and balancing them. The randomized algorithm from Clarke et al. (1964) was used for the initial solution. Later on,  extended this work with the application of pareto ant colony optimization (PACO), which improves the deficiencies of the local search approaches in the information flow among solution sets. Borgulya (2008) presented a multi-objective evolutionary algorithm. The algorithm uses an explicit collective memory method, extended virtual loser (EVL), an evolutionary algorithm in which recombination is omitted in favour of an adaptation of the mutation operator. Jozefowiez et al. (2009) proposed a metaheuristic based on an evolutionary algorithm with classical multi-objective operators and to improve efficiency, new operators were implemented for evaluating individuals in the population and descent selection. To favour search diversification, an elitist selection mechanism was implemented. Oyola et al. (2014) proposed a formulation for VRPRB where the balance of the routes points to optimizing equity measuring the range (length to the largest). Tabu search and GRASP-ASP based heuristics were used to solve the problem. Reiter et al. (2012) researched the capacitated vehicle routing problem (CVRP) with two objectives, (1) minimization of the total trip cost and (2) minimization the route length. They presented a solution formulation that is an extension of that presented by Laporte et al. (1984) and Achuthan et al. (1996) for the distance-constrained CVRP (DCVRP), which consists of applying an adaptive method of constraint-based on a hybrid model of a branch-and-cut and two genetic algorithms (NSGA-II and GA). Sarpong et al. (2013) proposed a generator column for bi-objective vehicle routing problems with min-max. One of the objectives was defined as a function of minimums and maximums that is optimized through generating columns that have a strategy to make them more effective. One of the most recent articles on the bi-objective VRPRB was presented by Lacomme et al. (2015). Their solutions were generated from an acceptance criterion to find a set of non-dominated solutions for VRPRB through the TSP (indirect solutions); this step was achieved by applying a division algorithm or splitting algorithm that considers both optimization criteria. To improve the solutions, a method based on a path relinking algorithm is used.

Formulation of the GVRPPCWE
The strategy for constructing routes that solve this routing problem is based on Granada et al.'s (Granada, et al., 2019) proposal. Its methodology allows constructing topologies that eliminate generating sub tours because a set of restrictions is proposed focused on obtaining valid solutions. The topologies formed are graphs with nodes of lesser degree equal to two. In this problem we have the following considerations: i.
There is a single deposit, where the private vehicles will return after visiting the last customer, the common carrier should not return to the deposit. ii.
Deposit capacity is unlimited or enough for the demand. iii.
Vehicles have limited and uniform capacity, including common carriers. This causes the problem to be homogeneous. iv.
The set of customers have a known demand and the sum is greater than the capacity of the private fleet . v.
The network associated with the problem consists of two complete graphs. The first is associated with the arcs travelled by the private fleet's routes and the second with the common carrier ones. vi.
Each customer will be visited only once with any of the vehicles, either private fleet or common carrier. vii.
Vehicles on common carrier routes may only carry out one journey and their journey ends at the last customer served.
O emissions are obtained directly by calculating vehicle fuel consumption. The model used in this work was the one implemented in Toro et al. (2017), which is defined through the nature of the forces that intervene in the movement of vehicles. To determine the equity criterion, and are defined as the distance travelled by vehicles and , as the set of routes where , ∈ . The quantifiable value of equity was obtained from the difference between the largest and shortest distance of all routes (Galindres, et al., 2017). This criterion is called range and is described in the following equation as: The GVRPPCWE was obtained based on the formulation proposed in Castaneda et al. (2019) and was defined as a mixedinteger linear problem. The problem is represented as a complete graph of the form = , , where represents a set with all the vertices of the graph (deposit and customers) and the set of arcs connecting the customers and deposit. The set = {1,2,3, … } contains the customers' nodes (deposit will have node zero assigned). Variables are defined as: Sets: : customers. : deposits (In this case it is unique) : nodes, = ∪ .
Decision binary variables: : cost associated with travel between nodes − . : variable that represents the activation of the arcs travelled by the private fleet. : variable that represents the arc of return to the deposit. : variable that represents the activation of the arcs travelled by the common carrier.
Parameters: : a distance between customers and . α: a constant function to the inclination, empty vehicle weight, energy for transitory speed, the resistances of wheels and wind and the internal losses of the vehicle. γ: a constant function of the slope of the terrain and the resistance of the wheels.
: penalty for using the common carrier. : a flow of goods between nodes − for both fleets.
The functions for cost, environmental impact and equity are defined as:

= (5)
Therefore, the model is stated as follows: The description of the equations that make up the constraints can be found in detail in Castaneda et al. (2019). Eq.  [S0] ← Perturb (S1);

Solution method
To solve the GVRPPCEW, the ILS metaheuristic was used, which implements the generation of subsets in the solution space around local optimal solutions, in which local search procedures will be executed applying appropriate neighbourhood schemes. Diversification procedures were applied to leave the local optimum, allowing access to other solution spaces Rendon et al. (2015). The characteristics of the random variable neighbourhood search (RVNS) were based on the random use of search operators in the so-called intra-route and inter-route neighbourhoods. The algorithm comprises the following steps: (1) generate the local initial solution, (2) local search, (3) perturbation, (4) dominance and (5) acceptance criteria. The pseudocode is presented in Algorithm 1.

Initial Solution
After generating the initial solution from the savings algorithm, the solutions generated by the ILS throughout its execution are stored for later use as initial solutions. The selected solution was assigned a probability of being perturbation, allowing diversification of the starting points Algorithm 1 shows the pseudo-code for selecting the ILS local initial solution. For constructing the heuristics, the concepts presented by Chu (2005) and the heuristics of Clarke et al. (1964), called the savings algorithm, were used. This heuristic is not implemented as a multi-objective. In this search for initial solutions that minimize only the cost, the other two objectives are minimized through the ILS local search, since trying to incorporate them into an initial heuristic would be of little use since the local search of the multi-objective algorithm can more efficiently reach solutions that minimize the objectives simultaneously. The implementation of the heuristic for the first initial solution is divided into four stages: (1) construction of routes for the private fleet and selecting customers for the subcontracted routes, which is carried out in two phases. (2) inserts. This stage also runs in two phases. (3) application of the classic savings algorithm for selected customers for the common carrier fleet. The insertion of customers to the routes is done in parallel; this allows insertions in all the routes created in advance of the algorithm. (4) improvement stage.

Random variable neighbourhood search (RVNS)
This function groups the local search, with the characteristic of neighbourhood construction chosen randomly from a defined set = , , … , . These neighbourhoods are built through inter-route exchange structures, where each neighbour corresponds to a feasible solution, formed by exchanging customers among routes. RVNS scans neighbours and stores the best solution. The procedure is described in Algorithm 3.
In line (4), the inter-router operator is randomly selected, which is applied at line (5). Line (7) contains the improvement conditional, which indicates that all the objective functions of the solutions found must be improved by the inter-route operators. This conditional makes each solution of the pareto be improved by another that strictly dominates it. The conditional of line (13) only verifies improvement in the objective function of cost since, because customer movements are on the same route, emissions decrease proportionally and the equity criterion does not vary significantly. Operators that do not improve the solution are eliminated in line (14). After intra-route improvement, inter-router operators are restarted to start the search for the best solution again. In line (25), intra-route operators that do not improve the solution are eliminated. The algorithm ends when there are no inter-route operators available.

Search methodology of the multi-objective ILS algorithm
ILS is a search algorithm with a high degree of downward trend because it focuses its search on the location of a specific solution and through the construction of neighbourhood structures, it can intensively explore solutions that are found nearby.
In the search methodology for problems with several objectives, the factor that indicates whether one solution is better than another is indicated by dominance but this could not be used in the local search because, by the definition given for dominance, it accepts solutions that are not strictly better; the search might not converge towards the optimal solutions and get caught in search loops. For this reason, the algorithm is defined by the condition: one solution is better than another in local search when it is strictly better on all objectives simultaneously. This is done by the instruction in line (7) of Algorithm 3. This implies that the ILS algorithm performs a local search for each solution with inter-router operators, which obtain solutions only in their space or region of dominance. In this way, the algorithm improves the pareto front. To easily visualize the search space, Fig. 2 shows the space where the solutions obtained by the ILS algorithm are located in a problem with two objectives. For three or more objectives, descent vectors point to the origin for minimization problems on all their objectives.

Tabu movement
The implementation of an algorithm like ILS involves making many movements between routes and verifying the value of the objective functions. For this reason, it is necessary to avoid performing repetitive movements and calculations that do not improve objective functions. Therefore, the tabu search metaheuristic, which implements a short-term memory causing a part of the movements to be temporarily prohibited, was used as a reference. In this work, tabu movements are defined as those that have already been verified by intra-route exchange operators and do not improve the objective function. Because the ILS applies a perturbation procedure to some or all the routes, the modified routes will be enabled again, removing this restriction.

Exchange operators
Interchange operators are structures in charge of creating neighbourhoods through exchanging customers between routes. The choice of each structure is made randomly. The exchange structures are described as: Inter-route operators: the shift operator transfers one or more customers from route to route . The swap operator performs exchanges between a pair of routes and . Unlike the shift operator, it ensures a simultaneous exchange of customers. The K-shift operator consists of transferring consecutive customers from route to the end of route . Intra-route operators: the Or-opt operator selects a customer to be removed from his current position to be inserted in any other position on the route. The Or-opt(n) operator has adjacent nodes change position on the route. The 2-Opt operator executes a movement that consists of removing two non-adjacent arcs and entering them again in a crossed way. Exchange performs a position permutation movement between two nodes.

Perturbation
The perturbation mechanism allows exploration in different spaces of the solution and consists of applying intra-route movements randomly. This is called diversification. These movements allow leaving local optimums and reaching other points in the solution space. The disturbance movements were applied by randomly selecting multiple inter-route interchange operators. The perturbation operators are described as multiple-shift (1,1), a non-simultaneous but random exchange operator of two customers between a pair of routes, and multiple-swap (1,1), unlike the other operator, these perturbance exchanges must be carried out simultaneously, which means that the transferred customers exchange positions, so both movements must be feasible.

Perturbance movement restriction
The perturbance procedures can mean that the solution used to perform local searches is not in relatively close neighbourhoods; this affects the performance of the ILS algorithm since increasing diversification loses search intensification. To balance these aspects, a distance criterion that limits the disturbance movements was implemented. This criterion consists of generating a distance from a route reference point and a candidate node to transfer or exchange. For this, the centroid of each route was determined. The reference distance was calculated between the centroid of the route and the deposit. It is defined that in the centroid of each route there will be a radius that draws a circumference, the nodes of other routes that are inside this circumference are called perturbing to the route, that is, they can be transferred or exchanged.

Dominance
For construction solutions in algorithms with > 1 objectives, the definition of dominance, which indicates whether one solution is better than another, is necessary. When one solution is better compared to another, it is called a dominated solution, the opposite case is non-dominated. The concept of dominance was defined by Deb (2014) as follows: One solution dominates another solution when the following conditions are true:  Solution is not worse than solution on any of the objectives.  Solution is strictly better than solution on at least one of the objectives. If any of the conditions is not met, the solution does not dominate solution .
For the selection of the non-dominated, each iteration of the ILS algorithm generates a new candidate solution to enter the pareto front. To determine if a solution obtained can enter the pareto, the following procedures were used:  For the candidate solution to enter the pareto, its region of dominance is verified. There should be no solution in the pareto that dominates it. If it is true, the solution can enter the pareto; otherwise, it is discarded.


If the solution entered the pareto, verification of its lower area is performed. All existing solutions that are in this region will be removed from the pareto because they will be dominated by the new solution that enters.
Algorithm 4 shows the procedure for selecting non-dominated solutions.

Pareto Metrics
In multi-objective problem solving, pareto fronts or frontiers are generated. This set of solutions creates difficulty in comparing with each other. Hence the need to use indicators that allow comparison of pareto fronts and that allow determining when one is better than the other. It is important to clarify that pareto metrics were not used as part of the optimization algorithm in the improvement principle in this work. These were used only to compare the results obtained. The metrics used in this work were obtained from the implementation by Tian et al. (2017). Pareto metrics were not used as incumbent on metaheuristics. These were used only to make the comparison between pareto fronts obtained after executing the ILS algorithm. For all metrics, a reference pareto front or at least one reference point is required on which the best solutions are determined. In this work, the reference point was defined as the Nadir point (Deb, 2014). This is defined from the maximum values of the set of pareto solutions obtained. Hypervolume is one of the indicators most commonly used as a pareto metric (Bradstreet, 2011). A pareto front is considered as points in the solution space; the hypervolume is the n-dimensional space contained among another front and a reference. In this implementation, the higher value hypervolume was a better metric since this value indicates how large the volume is between the front and reference point, which is of lower quality in this case (Nadir point). Spacing is a metric obtain from Schott (1995); it indicates how spaced the solutions of a Pareto front are, so it is defined as a density metric. The value is calculated by measuring the distance of each solution from its neighbours. The spread metric (Deb, et al., 2002), is used as a distribution and diversification parameter of the pareto front. Values close or equal to zero indicate that the solutions are continuous and evenly distributed among the extreme solutions.

Parallel processing
In this implementation, parallel processing is an important tool for building pareto fronts in less computational time. In this way, it is possible to use all the available processing capacity for constructing the pareto. Parallelism in processing consists of generating simultaneous and independent execution processes from the same initial solution, generating various pareto fronts that are unified to obtain a single front. To perform the unification of the various paretos, it is necessary to apply a global dominance procedure, which consists of unifying all the paretos, applying the dominance function. In the end, there will be only one pareto compound of non-dominated solutions. The density of the pareto front operates as a function of the depth of the defined local search. The greater the number of iterations defined for local search and the number of parallel processes, the denser the front. Fig. 3 shows the scheme of parallel processing for multi-objective ILS.

Analysis of results
The algorithm was implemented in C++ (Microsoft Visual C++ 2015) and Matlab®, running on a computer with a 64-bit Microsoft Windows® operating system and an Intel® Xeon™ CPU ES-2683 v4 @ 2.10Ghz, 16 processors and one memory 64GB RAM. For the implementation parameters, the formulation proposed in Penna et al. (2013) was used as a reference. The parameters used were: = 100 * where: : defines the number of local iterations. : establishes the relationship between the number of local iterations and the number of vehicles in the problem. For each iteration, this value will be randomly selected within the domain [1,20].
: number of vehicles from the private fleet and common carriers. : depth factor of the local search. For dense pareto fronts, it is defined as = 4, but for instances of high computational complexity (More than 50 customers), it is defined as = 2 to ensure convergence in reasonable times. : maximum number of iterations for the ILS.

ILS Search Stability Verification
To verify the stability of the algorithm concerning the solutions, the central limit theorem was used as a reference. In the process, = 30 executions were performed for instance A-n39-k6, obtaining solutions with a low coefficient of variation, which allows confirming that the algorithm is stable and descending. Table 1 shows the values obtained for the runs. The values for standardized bias and standardized kurtosis are in a range from [-2,2], indicating that the data associated with the pareto metrics do not have a significant deviation from normality. Fig. 3 shows a boxplot graph of the results of the normalized metrics with an outlier for time. Instances from the specialized literature for the classic CVRP problem described in Uchoa (2017) were used, classified as complex for each objective, where denotes the number of vehicles, number of customers and considering the problem symmetric. In the instances, the value of corresponds to the number of private fleet vehicles, calculated as = 0.8 / , where is the value of total customer demand and is the capacity of the vehicles. The vehicle fleet is homogeneous and the common carriers' routes have no restriction on their maximum vehicle limit. For the comparison of the pareto fronts found in each execution, the hypervolume metric was used to rank them from best to worst by their highest value, since this is larger as the pareto moves away from the reference Nadir point. This metric uses normalized values. The use of several metrics to classify the solutions by their quality could be inappropriate since the metrics are defined on different geometric concepts that imply that in some cases they are not in classification if the reference is a point and not a pareto. In the case of spacing, this metric is used only to rate the density of the front, not its quality. Table 2 shows the results obtained for instances defined as low mathematical complexity, with 50 customers and up to 10 vehicles. For these instances, 10 executions were carried out for each and the best solution was selected. The implemented instances with high mathematical complexity, with a greater number of customers than the previous instances, had 200 customers and 16 vehicles. For these instances, 5 executions of each were run to obtain results in reasonable times. Table 3 presents the results obtained. The pareto fronts for the solutions obtained are shown in the Appendix.

Conclusions
This investigation was proposed to solve the GVRPPCWE. The proposed metaheuristic uses the basic concepts of the ILS, intensification through local search and diversification over perturbance. The construction of solutions for the multi-objective problem was achieved by selecting the solutions obtained using dominance criteria and applying parallel processing. The verification of the results was carried out using pareto metrics, which allow comparison among fronts. The results obtained have no reference in the literature. The methodology that was implemented from the ILS algorithm starts from the construction of a first initial solution with the savings heuristic and the feedback of solutions that allow obtaining various starting points. This metaheuristic works in a descending method, finding the best solutions that improve the three objectives simultaneously, with a local search for solutions that are performed through the RVNS, selecting the solutions that are strictly better ( Ψ * < Ψ ), avoiding loops search and ensuring their convergence. Additionally, for local search, tabu movements were used that increase efficiency by not allowing repetition of movements that do not improve objective functions. Due to the high complexity of the problem, which implements open and closed graphs for the fleets and the minimization of three objectives, it was necessary to implement parallel processing for constructing the Pareto fronts. This methodology allows the use of all the computing capacity to obtain the solutions, drastically reducing the execution time and making it possible to obtain solutions in instances of up to 200 customers in reasonably short times. To verify the stability of the algorithm, an instance was selected to be executed 30 times according to the central limit theorem, obtaining that the results show stable behaviour and short variation and showing results with a normal probability distribution for the time of execution, hypervolume, spacing and spread of the fronts.
For instances of low mathematical complexity, defined as instances of up to 50 customers and 10 vehicles, the fronts are obtained with times in a range of fewer than 5119 seconds. For instances of high mathematical complexity, defined as up to 200 customers with 16 vehicles, an exponential jump in time was obtained for instances of more than 101 customers and 4 vehicles. These times are considered very low for the complexity associated with using two types of fleets and a multiobjective problem. This was accomplished by implementing an efficient search algorithm. The hypervolume metric was implemented to compare the different fronts obtained, using the Nadir point obtained from the combination of all the fronts as a reference. Other metrics are not used for pareto quality verification since the true (optimal) pareto is not available and these may not match the hypervolume because they are defined geometrically differently. The spacing and spread metrics are indicators of the quality in the construction of the front. The hypervolume values can be used as a reference for future research to improve the results. The spacing metric gives a reference for how far apart the solutions contained in the fronts are and is a value that characterizes the fronts obtained. According to the results obtained for the spread metric, the fronts obtained are not complete between the extreme points, since the methodology does not focus on finding solutions by steps or intervals. The range of the results is mostly ≤1, indicating that the fronts are of medium uniform distribution with some distant points.
In future works, since there is no specialized literature in the direction of this research, the results can be used as a reference to be compared with the implementation in other metaheuristics and exact models.