A biobjective capacitated vehicle routing problem using metaheuristic ILS and decomposition

Vehicle routing problems (VRPs) have usually been studied with a single objective function defined by the distances associated with the routing of vehicles. The central problem is to design a set of routes to meet the demands of customers at minimum cost. However, in real life, it is necessary to take into account other objective functions, such as social functions, which consider, for example, the drivers' workload balance. This has led to growth in both the formulation of multiobjective models and exact and approximate solution techniques. In this article, to verify the quality of the results, first, a mathematical model is proposed that takes into account both economic and work balance objectives simultaneously and is solved using an exact method based on the decomposition approach. This method is used to compare the accuracy of the proposed approximate method in test cases of medium mathematical complexity. Second, an approximate method based on the Iterated Local Search (ILS) metaheuristic and Decomposition (ILS/D) is proposed to solve the biobjective Capacitated VRP (bi-CVRP) using test cases of medium and high mathematical complexity. Finally, the nondominated sorting genetic algorithm (NSGA-II) approximate method is implemented to compare both medium-and high-complexity test cases with a benchmark. The obtained results show that ILS/D is a promising technique for solving VRPs with a multiobjective approach.


Introduction
The capacitated vehicle routing problem (CVRP) consists of the distribution of goods from a depot to a set of clients. Each client needs a certain amount of goods; for this, we have a fleet of vehicles, and a set of routes associated with a single vehicle must be designed in such a way that the demand of each of the customers is met taking into account the following operational restrictions: i) the vehicle capacity must not be exceeded, ii) in each proposed solution, a customer must be served by a single vehicle, iii) each of the vehicles in the fleet starts at the depot, delivers the respective demands, and returns to the depot. The objective function is to minimize the total distance traveled by designing a set of routes taking into account the constraints. The CVRP is an NP-hard combinatorial problem (Irnich, Toth, & Vigo, 2014). Vehicle routing problems (VRPs) have usually been studied with a single objective function; however, in real life, it is necessary to take into account other objective functions, such as social functions, which consider, for example, the drivers' workload balance.
In particular, in this paper, the biobjective CVRP is addressed, taking into account the economic and social objectives. The economic objective is to minimize the cost associated with the CVRP route design; the social objective is to balance the workload of each of the drivers (Matl, Hartl, & Vidal, 2017). The formulation of the social objective must satisfy several restrictions established by union contracts and company regulations (for example, work periods during the day, number and duration of breaks during service, maximum duration of driving periods and extra hours). In addition, the allocation must be equitable among the drivers, so it must take into account factors such as traffic conditions and the number of customers assigned to each vehicle (Lee & Ueng, 1999). Furthermore, in some cases, depending on the characteristics of the goods and services, the concept of equity is essential (Schwarze & Voß, 2013).
Several works, that introduce the concept of equity in the VRP, are highlighted. In (Matl et al., 2017), an exhaustive review and categorization of the literature related to VRPs and equity is performed, taking into account objective functions classified into five groups; the first is based on the CVRP with route balance (CVRPRB), the second adds time windows to this basic model, the third is called a "min-max VRP", the fourth includes complex variants of the VRP, such as periodic VRP or arc routing problems, and the fifth shows a wide selection of applications based on the VRP. In addition, (Matl et al., 2017) present a theoretical analysis in which in the case of VRPRB, there is no substantial discussion about the relative merits or limitations of the model; however, the decision to use the duration of the trip as the equity metric seems to be justified in practice.
Halvorsen-Weare and Savelsbergh (2016) present the only study, so far, to examine the effects of different equity functions in VRP solutions in a biobjective environment. Specifically, the authors consider the mixed capacitated general routing problem (MCGRP) as a variant of the CVRP and the capacitated arc routing problem (CARP). The biobjective MCGRP minimizes the cost and an equity objective based on four measures, range, maximum travel, absolute deviation with respect to the average length of the route, and absolute deviation with respect to the total length of the route. Some of the computational experiments suggest that the adoption of the balance measure in terms of the difference between the longest and the shortest route has good performance in several instances. In the same way, the number of Pareto front solutions can vary substantially depending on the equity function chosen. Schwarze and Voß (2013) formulate the model for the skill VRP taking into account the load balance and the reuse of resources to obtain a properly distributed workload for vehicles.
Others review different existing approaches in the literature that measure the load balance in VRPs (Matl et al., 2017;Schwarze & Voß, 2013). Basically, there are six different types, min-max is the simplest measure of existing inequality but has several disadvantages; for example, it does not distinguish between distributions with identical worst cases. Lexicographic min-max is an extension of the first case in which not only the worst case is minimized but also the second case (subject to the minimization of the first) and so on. Although it solves the problem of the previous case, it does not quantify the remaining differences. Range is the difference between the worst and the best case and provides more information. Mean absolute deviation (MAD) is affected by all the values in the solution and not only by the extreme values and is therefore "stronger" than the range. Standard deviation, as mentioned in (Matl et al., 2017), is possibly the best known statistical measure of dispersion and satisfies the Pigou-Dalton (PD) principle of transfers. However, its biggest disadvantage is the computational complexity. The Gini coefficient (GC) is used in economics and in inequality studies in general. The GC assumes values between 0 and 1, with values close to 0 corresponding to low inequality. However, its biggest disadvantage is its computational complexity.
Finally, several applications of the VRP have been presented, including the social dimension. For example, in (J. Q. Li, Borenstein, & Mirchandani, 2008), several social benefits are included in the programming model for solid waste collection in the city of Porto Alegre (Brazil). A model of location (recycling facilities) and routing of vehicles solved with a heuristic is formulated, delivering better programming in terms of the social aspects of the solid waste collection program. In (Lee & Ueng, 1999), the research focuses on workers, or more precisely on drivers. Given that employees are the most valuable resource of a successful business, a load-balancing factor is taken into account that indicates their well-being, improving the company's competitiveness. For the above, several considerations are assumed: the size of the fleet of vehicles, a linear relationship between time and distance traveled, the demand of each consumer, and the service time at the demand points.
With regard to the simultaneous formulation of the dimensions mentioned above for the VRP, several approaches have been taken into account, with an emphasis on the solution method. For example, in (Oyola & Løkketangen, 2014), the CVRP is solved taking into account costs and route balancing simultaneously. A GRASP-based heuristic is used to obtain the approximate Pareto front. Sapong et al. (2013) in contrast, present an exact mathematical model for this type of problem that the CVRP solves with a multiobjective approach based on the generation of columns, where one of the objectives is a linearized min-max function that is calculated by generating columns. Finally, Bertazzi et al. (2014) present an analysis that validates the use of the min-max objective function to solve the multiobjective CVRP problem, differentiating it from other approximate functions that measure the balance of routes as min-sum. The analysis involves the study of the worst case, identifying the artificial balance that can occur in some routes of the solution. In addition, the generation of the approximate Pareto front using min-sum and min-max is a promising research area.
The proposed method solves the multiobjective CVRP problem, which takes into account both costs and workload balance simultaneously (bi-CVRP), where the metaheuristic used is based on a trajectory called an Iterated Local Search (ILS), which includes the Variable Neighborhood Search (VNS). The ILS is composed of two stages that maintain diversification and intensification in the search space. The diversification stage is performed through a perturbation mechanism that makes it possible to explore promising regions in the solution space. The intensification stage is implemented with the VNS, which consists of specialized operators responsible for reducing the solution space by performing the search in nearby surroundings (neighborhoods) of the current solution. In the present work, the VNS is implemented by using two types of operators, interroute operators, which look for a better solution between two routes, and intra-route operators, which look for a better solution for only one route. Both operators are based on shift, swap, and 2-opt strategies.
The decomposition process is used to solve multiobjective problems that decompose a multiobjective optimization problem into a number of single-objective subproblems and optimize them simultaneously (Trivedi, Srinivasan, Sanyal, & Ghosh, 2017). In (Zhang & Li, 2007), the decomposition strategy is used within evolutionary algorithms, which means that this procedure has a lower computational complexity in each generation than the multiobjective genetic local search (MOGLS) and the nondominated sorting genetic algorithm (NSGA-II). The experimental results show that the multiobjective evolutionary algorithm based on decomposition (MOEA/D) with simple decomposition methods, outperforms MOGLS and NSGA-II on multiobjective 0-1 knapsack problems and continuous multiobjective optimization problems. It has been shown that MOEA/D, using objective normalization, can deal with disparately scaled objectives, and MOEA/D, with an advanced decomposition method, can generate a set of very evenly distributed solutions for three-objective test instances.
In this paper, the ILS metaheuristic and decomposition strategy (ILS/D) is used to solve the bi-CVRP; in addition, to make a fair comparison, we implement an exact method based on decomposition, and implement the well-known NSGA-II proposed by Deb et al. (2002). The methodology is validated with benchmark instances from the literature taken from (Augerat, Belenguer, Benavent, Corberán Naddef, & Rinaldi, 1995). For the analysis of obtained fronts, different performance metrics are used to determine both the quality of the solutions and the validation of the proposed method. From the results obtained, it can be concluded that the proposed method is competitive compared to the other algorithms, and it is possible to apply it to other routing problems, which should be of interest to the academic community.
The rest of the article is organized as follows: Section 2 presents the new solution method to solve the bi-CVRP using ILS/D. In Section 3, the results are analyzed using the instances from (Augerat et al., 1995) and compared with those obtained by the exact, and NSGA-II, approaches. Section 4 presents the conclusions, considerations, and guidelines for future works.

Solution Method
To solve the biobjective CVRP, it is proposed initially to implement an exact method, which will serve as a reference to carry out a comparison with the results obtained using both the proposed method and the NSGA-II metaheuristic.

Exact Method
To solve the problem in an exact way, a mixed integer linear model (MILP) was defined in which the economic objective function consists of minimizing the costs associated with the design of the routes. To define the social objective function, we took into account what was defined in (Lozano, Gonz, Rodriguez-tello, & Lacomme, 2016), where seven objective functions related to the route balance were studied, and it was found that the most frequently used ones were max-min (also called range), which minimizes the difference between the maximum and minimum length of routes; min-max, which minimizes the maximum length path; and var, which minimizes the variance of the lengths of the routes. In a recent study, Matl et al. (2017) analyzed the objective functions mean absolute deviation (MAD), defined as the average absolute difference between the length of each route and the average of the route lengths, and the Gini coefficient, which is one of the most commonly used measures for inequality studies. In (Galindres-Guancha et al., 2018), the quadratic objective function var was used, showing satisfactory results with a high computational cost.
To solve the biobjective CVRP in an exact way, we used the traditional cost function and the linear objective function minmax, which seeks to minimize the maximum path length of the solution, as shown in (1).
ij c and ij x correspond to the cost and existence of the arc between nodes i and j. i  corresponds to the length of the route i.
The optimization process will use the decomposition concept, whose objective is to generate the Pareto front. Decomposition is a widely used strategy in multiobjective optimization problems (Zhang & Li, 2007), and consists of dividing the multiobjective problem into N scalar optimization subproblems.This method transforms multiple objectives into an aggregated objective function by multiplying each objective function by a weighting factor and summing up all weighted objective functions. In addition, to find a well-distributed Pareto front, for example, in (Samanlioglu, 2013), the weighted Tchebycheff aggregation method is used. In this paper, the weighted Tchebycheff method is also used, and equally spaced vectors  are generated using the method implemented in (Zhang & Li, 2007). According to (Marler & Arora, 2004), this formulation provides a necessary condition for Pareto optimality. The formulation that includes the weighted Tchebycheff aggregation is shown in Eq. (4): The parameters are generated, and Eq. (5) is modified by obtaining the minimum value of each objective function.

Approximated Method
To solve the problem in an approximate way, we propose a novel method based on the metaheuristic ILS, and again the concept of decomposition (ILS/D), to generate the front. In this paper, an adaptation of the ILS and reduced variable neighborhood descending techniques (ILS-RVND) shown in (Silva, Subramanian, Vidal, & Ochi, 2012) is made, which consists of performing a local search starting from an initial solution to improve this solution, and then applying a perturbation mechanism to obtain a new solution. This process is performed in a certain number of iterations. The local search is based on the construction of neighborhoods whose solutions are accepted according to a criterion. This method has been shown to be efficient in this type of problem (Subramanian, Ochi, & Cabral, 2008). Decomposition is used as in the previous case; each subproblem corresponds to the initial solution used in the ILS-RVND.
The N subproblems are optimized using the MOEA/D framework, as shown in (Zhang & Li, 2007). In addition, two modifications to the previous framework were taken into account to increase its efficiency: i) MOEA/D-GR (Global Replacement) is used to carry out a more adequate update in the neighboring subproblems   B i . In this scheme, we first find a more appropriate neighborhood to be replaced by the candidate solution i x . The details are shown in (Wang, Zhang, Gong, & Zhou, 2014). ii) MOEA/D-GRA (generalized resource allocation) is used to dedicate more resources to the subproblems that are more likely to improve. In each iteration, computational resources are assigned to some subproblems according to the vector PoI. Specifically, the resource allocation strategy, Online, is used to define a utility function based on the relative improvement in the last T  iterations. The details are shown in (Zhou & Zhang, 2016).
The details of the proposed ILS/D are shown in Algorithm 1.
Line 1 of Algorithm 1 sets a vector of N elements weighted and equally spaced, where N is the number of subproblems. Line 2 defines the size of the neighborhood for each subproblem; usually, T is a fraction of the total number of subproblems . In lines 9 and 10, the criteria for stopping the algorithm when a maximum number of iterations is reached are defined. In lines 11 to 19, the ILS algorithm is executed for each subproblem, taking into account that the subproblem with the highest probability of improvement is optimized (line 12). In line 13, the improved solution ' i x for subproblem i is obtained using the ILS algorithm (Subramanian et al., 2008 shown in the strategy proposed by (Wang et al., 2014), to maintain the balance between convergence and diversity in the optimization process. To find the most appropriate subproblem j for ' i x , the following expression is used:

Non-Dominated Sorting Genetic Algorithm (NSGA-II)
The NSGA-II was proposed by ( Deb et al., 2002). It sorts a population of size N according to the level of nondomination and the crowding distance. The NSGA-II uses selection, diversification, and mutation operators to create half of the population following the selection of the best solutions. For most problems, the results show that the NSGA-II can find diverse solutions and good convergence close to the optimal Pareto front in comparison to the multiobjective evolutionary algorithms (MOEAs). Details of the NSGA-II are shown in (Gallego Rendón, Toro Ocampo, & Escobar Zuluaga, 2015).
The NSGA-II has been used to solve the VRP with a multiobjective approach. For example, (Jozefowiez, Semet, & Talbi, 2005) address a biobjective VRP in which the total length of routes is minimized as well as the balance of routes, that is, the difference between the maximal route length and the minimal route length. In addition, to improve the efficiency of the NSGA-II, parallelization and elitist diversification strategies have been added, indicating a strict improvement of the generated Pareto set. Finally, in (Sun, Liang, Zhang, & Wang, 2017), a memetic algorithm based on NSGA-II is used to solve the biobjective vehicle routing problem with route balance (VRPB), which considers the total length of routes and the balance issue among different routes, minimizing the maximal route length.

Performance Metrics
To measure the quality of the front, we used the most common metrics in multiobjective optimization that address the indicators of both convergence and diversity (K. Li, Wang, Zhang, & Ishibuchi, 2018) (Riquelme, Von Lücken, & Barán, 2015):

The inverted generational distance (IGD)
IGD shows the deviation of the smallest distances between the points of the optimal front with respect to the approximate front (Jiang, Ong, Zhang, & Feng, 2014), defined as shown in (6).
where P is the approximate front, S is the optimal front, i d is the minimum distance between each point of P and the closest solutions in S, and 2 q  .

The hypervolume (HV) indicator
HV classified in (Jiang et al., 2014) as a metric that measures convergence and diversity simultaneously, was proposed by (Zitzler & Thiele, 1999) and is shown in (7). In particular, the closer the solutions of S are to the true Pareto fronts, the larger the value of HV. At the same time, a higher HV could also indicate that solutions of S are scattered more evenly in the objective space. However, the downside of HV lies in the high computational complexity.

HV(S)
where i v is the hypercube formed by point i of front S and a reference point (nadir point).

Generational distance (GD)
GD is among the indicators commonly used in MOEAs to measure convergence with respect to the optimal front, according to (Jiang et al., 2014), and is defined as shown in (8).
where P is the approximate front and S is the optimal front; i d is the minimum distance from each point of the optimal front to each closest point of the approximate front, and 2 q  .

Coverage (C)
C is described by (Zitzler & Thiele, 1999), and the equation is defined as shown in (9). C is used to measure two sets of points (i.e. 1 S and 2 S ). The value 1 2 ( , ) 1 C S S  means that all points in 2 S are dominated by 1 S .

Computational Results
The instances used in this work correspond to those proposed in (Augerat et al., 1995). To analyze the results obtained, three scenarios were constructed. The first is used to compare the results of the exact method with a similar one available in the literature. The second uses a medium-sized instance to compare the results of both the ILS/D and the NSGA-II with those produced by the exact method. The third uses a large instance to compare the performance of the ILS/D with that of the NSGA-II, and it summarizes the results for the instances proposed in (Augerat et al., 1995), showing the performance of ILS/D and NSGA-II algorithms using the HV and coverage metrics to compare them.
First, the exact method was implemented using the CPLEX solver and the AMPL language. Due to the high computational effort in larger instances, it was only used with the instance of medium complexity, A-n37-k5, and, in addition, the results obtained in this instance were used for comparison with those obtained by the ILS/D and the NSGA-II metaheuristics. This instance is also analyzed in (Reiter & Gutjahr, 2012), in which the biobjective CVRP problem is solved using a hybrid technique based on the   constraint method and the multiobjective NSGA-II algorithm. Here, we used a similar model and the same objective functions, minimization of the total costs and minimization of the maximal length of the route (min-max).
Here, it should be noted that the extreme solutions for the A-n37-k5 instance of the Pareto front obtained using the exact method and shown in Fig. 1. are the same as those obtained in (Reiter & Gutjahr, 2012).
Both the proposed ILS/D and NSGA-II algorithms were coded in C++, and instances of medium and high mathematical complexity were run for 36 and 79 clients, respectively, taken from the benchmark presented by (Augerat et al., 1995).
Second, to validate the quality of the proposed ILS/D and NSGA-II algorithms, instance A-n37-k5 was used. This instance is initially solved with an exact method that will serve as a reference for comparative purposes, and whose front is presented in Fig. 2 in red. In addition, 20 independent runs of both ILS/D and NSGA-II metaheuristics were performed. The fronts with the best values of the HV indicator obtained by ILS/D and NSGA-II were compared with the front obtained by the exact method. The fronts presented in Fig. 2 marked with blue and green colors, correspond to the fronts for the proposed ILS/D and NSGA-II methods, respectively.  As shown in Fig. 3, for the HV indicator, it can be seen that both algorithms generate solutions with similar medians. ILS/D has more variability than NSGA-II but after the 75th percentile, the ILS/D solutions are better than the NSGA-II solutions on average. For the indicators IGD and GD, shown in Figs. 4 and 5, respectively, it is more evident that the solutions generated by ILS/D are more convenient in magnitude and variability. Therefore, ILS/D has greater reliability and robustness for both indicators.
Third, the fronts with the highest values of HV obtained for the larger instance A-n80-k10 of the group presented in (Augerat et al., 1995) are shown in Fig. 6. The front highlighted in blue corresponds to the one obtained with the ILS/D algorithm. The front highlighted in green corresponds to the front obtained with the NSGA-II algorithm. The statistical information of the ILS/D and NSGA-II algorithms for the A-n80-k10 instance is presented in Fig. 7 and Fig. 8. 20 independent runs were made with each algorithm. It should be noted, that for the calculation of IGD and GD metrics, the utopic point was taken as the reference front. As shown in Fig. 7, for the HV indicator, similar values are observed for both ILS/D and NSGA-II; however, it is observed that ILS/D presents less variability. In Fig. 8, when comparing the algorithms with regard to the GD indicator, ILS/D shows better behavior, generating solutions of smaller magnitude and variability. According to what was observed with the study of instances of medium and high mathematical complexity, and when they are analyzed with the ILS/D and NSGA-II algorithms, it follows that the ILS/D algorithm shows better performance in the HV and GD indicators. Finally, Table 1 shows the performance of both the proposed ILS/D and NSGA-II. The comparison uses the HV and coverage metrics. The first column shows the instance name and then shows the maximum, average and standard deviation obtained values for the HV metric. Next, it shows the obtained value to the coverage metric. In this column, s1 and s2 correspond to the fronts with the highest value for HV generated by the ILS/D and NSGA-II algorithms, respectively. The last column of Table 1 shows the run average time. The Avg HV column of Table 1 highlights the values for HV where ILS/D is better than NSGA-II.

Conclusions and Future Work
A multiobjective algorithm was implemented that solves the bi-CVRP considering the cost of the routes and their balance. In this algorithm, the decomposition concept was used. For the solution of the subproblems, the ILS was used. The results are initially compared with those obtained using an exact method and the well-known NSGA-II metaheuristic for an instance of medium mathematical complexity. Finally, the ILS/D is compared with the NSGA II, using instances of medium and high mathematical complexity, concluding a better performance of the ILS/D as compared to the NSGA II.
An exact model was proposed to solve the bi-CVRP problem that seeks to optimize, in addition to the economic costs, the workload of the drivers by balancing the length of the routes. The exact model corresponds to a multiobjective MILP formulation; to solve it, the decomposition strategy was used, which transforms a multiobjective problem into N singleobjective subproblems. The front obtained is used to validate the proposed methodology.
An algorithm based on decomposition taken from MOEA/D and ILS was implemented. The first decomposes the multiobjective problem into N mono-objective subproblems, and the second solves each of these subproblems. The algorithm implemented showed great efficiency in generating fronts, whose fundamental characteristics are diversity between solutions and convergence; that is, the performance metrics used in multiobjective optimization, such as HV and IGD, show an adequate distribution and convergence of the obtained fronts, with which an adequate behavior of the proposed algorithm in the solution of this problem is verified.
When using instances of 32 to 80 clients from the Augerat database, it is concluded, that the ILS/D algorithm that has been proposed, presents better performance when compared with the NSGAII algorithm for the analysis; the indicators HV and coverage, obtained the following percentages: Max HV 50.0%, Avg HV 63.6%, Coverage (S2, S1) 77.3%, which confirm the best performance of ILS/D.
As the main research contribution of this work, a methodology that uses ILS metaheuristics combined with the concept of decomposition is proposed. The conclusions show that the method is competitive compared to the traditional NSGA-II, especially in larger cases. Additionally, the use of an exact method in solving the problem is used as a reference to measure the degree of precision of the implemented metaheuristics.
In future work, the possibility of applying parallel processing is raised, because the decomposition strategy breaks down the multiobjective problem into N independent subproblems.