A Memetic Algorithm for the Cumulative Capacitated Vehicle Routing Problem Including Priority Indexes

: This paper studies the Cumulative Capacitated Vehicle Routing Problem, including Priority Indexes, a variant of the classical Capacitated Vehicle Routing Problem, which serves the customers according to a certain level of preference. This problem can be effectively implemented in commercial and public environments where customer service is essential, for instance, in the delivery of humanitarian aid or in waste collection systems. For this problem, we aim to minimize two objectives simultaneously, the total latency and the total tardiness of the system. A Mixed Integer formulation is developed and solved using the AUGMECON2 approach to obtain true efﬁcient Pareto fronts. However, as expected, the use of commercial software was able to solve only small instances, up to 15 customers. Therefore, two versions of a Memetic Algorithm with Random Keys (MA-RK) were developed to solve the problem. The computational results show that both algorithms provided good solutions, although the second version obtained denser and higher quality Pareto fronts. Later, both algorithms were used to solve larger instances (20–100 customers). The results were mixed in terms of quality but, in general, the MA-RK v2 consistently outperforms the ﬁrst version. The models and algorithms proposed in this research provide useful insights for the decision-making process and can be applied to solve a wide variety of business situations where economic, customer service, environmental, and social concerns are involved.


Introduction
In this paper, we study the biobjective variant of the Cumulative Capacitated Vehicle Routing Problem (CCVRP), a "customer-centric" variant of the classical Capacitated Vehicle Routing Problem (CVRP) [1] in which a fleet of k vehicles serves a set of customers by respecting their priority, defined as an index assigned to each node to indicate the preference to be served. Unlike the traditional VRP, which focuses on the impact that routing costs have on logistics and, in particular, in the transportation activities within the supply chain, the CCVRP rises as a particularization that covers broad objectives centered around service level issues. This problem is relevant in contexts where both customer satisfaction and company profits are prioritized. Due to the importance of sustainable business practices nowadays, there is a need to develop distribution strategies aimed at reducing the negative impact that transportation activities have on the environment. An important application can also be found in the context of emergency logistics, where the distribution of medical aid becomes crucial, results for instances up to 40 customers and extending the computational experimentation over the metaheuristics for instances up to 256 customers.
Considering the aforementioned discussion, to the best of our knowledge, no tailored approach has been proposed in the literature for the problem considering the risk-averse perspective. For this reason, in this paper ,we study the Biobjective Cumulative Capacitated Vehicle Routing Problem considering Profits (BCCVRP-Pr). To handle this new problem, we propose a mixed-integer formulation and a MA-RK procedure to efficiently deal with instances of reasonable size.
Elshaer and Awad [20] studied the diversity of variants of vehicle routing problems solved using metaheuristics, where eleven papers were identified using memetic algorithms. In a recent example, Li et al. [21] proposed a hybrid metaheuristic that combines a memetic algorithm, sequential variable neighborhood descent, and a revised 2-opt method, for the plug-in hybrid electric vehicle routing problem. In addition, Zhang et al. [22] developed a multiobjective memetic algorithm for the vehicle routing problem with time windows. To our knowledge, only Ngueveu et al. [2] applied memetic algorithms to the cumulative capacitated vehicle routing problem. They created two versions of their procedures to solve a single objective problem to minimize the sum of arrival times at the customers.
On the other side, memetic algorithms have been combined with random keys in some applications. One example is the hybridization of He et al. [23], where memetic algorithms were combined with a biased random key genetic algorithm and adaptive large neighbourhood search to solve a scheduling problem. Other applications to scheduling problems using memetic algorithms and random keys were proposed by Li and Yin [24] and Ghrayeb and Damodaran [25], among others. Although this combination of memetic algorithms and random keys have been used to solve sequencing problems, such as in the traveling salesman problem [26], this combination has not been used to solve complex routing problems such as the one presented in this paper.
The literature described above shows the novelty of the combination of memetic algorithms and random keys to solve complex multiobjective routing problems.
The remainder of this paper is organized as follows. Section 2 describes the proposed mathematical formulation, whereas Section 3 presents the algorithm developed. Section 4 reports the results obtained from the experimentation using a set of benchmark instances. Finally, Section 5 summarizes the major findings and proposes future research directions.

Mathematical Formulation
In this section, we present the formal definition of the BCCVRP-Pr as well as its corresponding mathematical formulation. First, the following parameters and variables are considered:

Parameters
• N: Number of customers • K: Number of vehicles available The BCCVRP-Pr can be formally defined as an undirected graph G = (V, A) where V = {0, 1, 2, ..., N} denote the node set and A is the set containing all arcs. The node 0 corresponds to the depot and the rest of the nodes form the set of customers V = {1, 2, ..., N}. Each arc (i, j) ∈ A has an associated travel time c ij and each node j ∈ V has an associated demand d j . A heterogeneous fleet of K vehicles is available, each with capacity Q k , k ∈ {1, 2, ..., K}. It is assumed that the vehicle's overall capacity is enough to satisfy the total demand of all the clients. In addition, to consider the priority, a precedence matrix P is defined in which p ij = 1 represents that customer i should be serviced before customer j and p ij = 0 means that customer i can be served after the customer j.
The following additional parameters are considered: • Q max : Maximum capacity of any of the vehicles • M: Maximum travel distance allowed (the same for all vehicles) The goal is to find k tours that have only node 0 in the first position, and to cover all the nodes, while minimizing the sum of the latencies of the trips. The demand of all customers must be satisfied without exceeding the capacity of each vehicle. Customers should be served (preferably) according to their priority level, minimizing the total tardiness of the system.
For each path, the tardiness arises when a customer with lower priority is served before a customer with higher priority (even if both customers belong to different routes). In other words, the arrival time (t i ) for a customer with a lower priority index is less than the arrival time of a customer with a higher priority index (t j ). Qualitatively, tardiness is associated with customer dissatisfaction. However, because latency is estimated as a function of distance, then tardiness is obtained as the difference between their arrival times (when p ij = 1), I ij = t j − t i . Hence, the total tardiness of the system is computed, as shown in Equation (1): Figures 1 and 2 graphically exhibit the conflict among the values of latency and tardiness for an instance that contains 12 nodes. The number above each node denotes the customer priority index (higher values indicate higher priorities). In Figure 1, the minimum total tardiness was obtained subject to the minimum total latency. Correspondingly, Figure 2 indicates the optimal total latency while assuring the minimum total tardiness.  The formulation is based on the model presented in [10] for the classical CCVRP and following the single flow formulation proposed for the multiple traveling salesman problem [27].
Variables w k 0j = 1, if the vehicle k is assigned to a path from node 0 to customer j, 0, otherwise.
x ij = 1, if the arc (i, j) is in the path of a vehicle, 0, otherwise.
y ij = number of remaining nodes after the node i on a route if x ij = 1; 0, otherwise.
v k 0j = the sum of demands of the nodes after node 0 on the route k if w k 0j = 1; 0, otherwise. v ij = sum of demands of the nodes after node i on a route if x ij = 1; 0, otherwise. t k 0j = Arrival time of node j from node 0 on a route k if w k 0j = 1; 0, otherwise. t ij = Cumulative time at node j on a route if x ij = 1; 0, otherwise. I ij = Tardiness in the arrival time to node i if node j is served first(p ij = 1); 0, otherwise. The corresponding mathematical model is stated as follows: subject to: ∑ j∈V In this formulation, the objective functions in Equations (2) and (3) seek a trade-off between the objective of the total latency and the total tardiness of the system. The constraints in Equations (4) and (5) ensure that exactly K arcs leave and return to the depot. The constraints in Equation (6) establish that a customer node can be visited for exactly one vehicle coming from the depot. Additionally, the constraints in Equations (7) and (8) impose that exactly one arc enters and leaves for each node associated with a customer. The constraints in Equations (10)-(13) help to avoid sub-tours and allow to calculate the position of the customer j on its respective route. The constraints in Equations (12) and (13) force variables y 0j and y ij to be zero when w k 0j = 0 and x ij = 0, respectively. Regarding the capacity of the vehicles, the constraints in Equations (14) and (15) allow establishing lower and upper bounds for the cumulative demand of any route. These constraints are derived from a generalization of restrictions in Equations (10) and (11). The constraints in Equations (16) and (17)  and v ij to be zero when variables w k 0j = 0 and x ij = 0, respectively. The constraints in Equation (18) ensure that the demand at each node i is fulfilled and in conjunction with Equations (16) and (17)

Reformulation Using Epsilon Constraint
In this subsection, we describe the characteristics of the model and the proposed reformulation. A fundamental task in multiobjective optimization is to find Pareto-optimal solutions. As a biobjective approach, we decided to implement a multiobjective method of resolution to generate an exact front of efficient solutions.
In the mathematical model, we note that the objective functions are separable. In other words, each of them involves different decision variables. On the one hand, y ij allows estimating the total arrival time to the customers. On the other hand, I ij computes the total tardiness for the case in which customers having a minor priority are served earlier than customers with higher priority index (regardless if they are located in the same route or belong to different routes).
To clarify this, the particular structure of the biobjective problem proposed herein is described below: subject to :

Equations(4)-(24)
The aforementioned characteristics of the biobjective model are exploited by using an improved version of the ε-constraint method, named as AUGMECON2 [28], as a solution procedure. For every single routing decision in Equations (4)- (18), the minimum tardiness (min T) problem bounded by the constraints in Equations (19)-(24) is solved as principal objective function, transforming the latency function (L) into constraint. The results of the proposed method are shown in Section 4.

Metaheuristic Algorithm
This section is devoted to describe the metaheuristic approach capable of obtaining high quality solutions for small instances and able to deal with large size instances. The proposed method is based on a Memetic Algorithm (MA), a population procedure that has shown its effectiveness in solving sizeable combinatorial optimization problems by incorporating a local search procedure within a classical genetic algorithm. This procedure has been successfully applied for addressing the CCVRP [2], introducing efficient move evaluation procedures in operations O(1) for some particular neighborhood structures. MAs have been also employed in solving other routing problems such as split delivery vehicle routing problems [29], capacitated location routing problems [30,31], vehicle routing problems with time windows [32], school bus routing problems [33], and green and healthcare routing problems [34,35].

Proposed Memetic Algorithm
Holland [36] was the first to propose Genetic Algorithms (GAs) inspired on ideas of evolution theory. Due to their simple and yet effective search procedure, several papers (e.g., [37]) describe their successful implementations in vehicle routing problems. In particular, Memetic algorithms (MAs) belong to the class of evolutionary algorithms that intensify the search by including local search within a classical genetic algorithm framework. According to Moscato and Cota [38], MAs are intrinsically concerned with exploiting all available knowledge about the problem under study. Due to this, a random key mechanism is included during the construction procedure in order to enhance the performance of the procedure. In this work, the MA proposed is adapted from the NSGA-II, successfully implemented in [39], and consists of the following procedures: initialization, recombination (crossover and local search), and classification in fronts. In our construction procedure, we include a random key that helps to generate a diverse initial population of feasible solutions of size N. Subsequently, during a predetermined number of successive generations (iterations), an offspring (P t ) of N individuals is generated from P t−1 through recombination and local search mechanisms, involving members representing tentative solutions (high-quality or non-dominated solutions) and members representing diverse solutions. After this, the individuals who belong to the previous and current generations are evaluated and grouped into fronts, according to the level of non-domination, as explained in [40]. To obtain the resulting offspring population of size N, the individuals are inserted into the set, starting with the one that belongs to the front of non-dominated solutions (F 0 ). Algorithm 1 shows the pseudocode of the overall MA with Random Keys procedure.

Constructive Procedure Based on Random Keys
The constructive procedure creates an initial population of feasible solutions based on generating a chain S a = {1, 2, · · · , n} and, for each customer, an auxiliary random key R a is used to encode the solution. Additionally, an empty set S p is used to save the temporary assignment of the customers to the routes.

Algorithm 1:
Memetic algorithm with random keys.

begin
it ← 0; Initialize a population (P 0 ) of σ chromosomes implementing the constructive procedure based on random keys; Sort P 0 in fronts following the non-dominated sorting approach; Generate an offspring population P it , of size N, from P it−1 , using selection, crossover and local-search operators; Combine parent and offspring population R it = P it ∩ P it−1 ; Sort population using the non-dominated sorting approach, identify fronts F j = (1, 2, · · · ), and calculate the crowding distance for each solution in F j ; Sort solutions in F j in decreasing order according to crowding distances, select the first N − |T it+1 | elements of F j and add it to T it+1 ; end end end Encoding mechanism: To encode the solution, a real number drawn randomly from [0, 1) is assigned to every single position in R a . Figure 3 depicts an example of this mechanism. Decoding mechanism: The decoding mechanism is applied based on the information of the random key R a . The R a chain is sorted in a non-decreasing order and their respective positions in S a are sorted correspondingly. As a result, a random ordered chain S a is obtained. Figure 4 exhibits the decoding mechanism. Assignment mechanism: In every iteration, the algorithm selects the corresponding jth customer in S a and systematically tries to insert it into into a temporary set S p in the first available position (procuring to maintain feasibility in the capacity of the vehicles). For instance, if in the first potential route, two customers have been previously assigned, the next open position will be the third one. In the case that the customer cannot be inserted in the route due to the capacity constraints, the insertion will be evaluated in the next available route. It is important to emphasize that, since the number of routes is given in advance, the construction procedure considers a parallel routing mechanism. In other words, it performs the evaluation of feasible insertions over all of the routes. If so, the algorithm continues by selecting the next customer at S a . Otherwise, the construction mechanism stops. If the algorithm reaches a feasible assignment, then S ← S p and the solution is inserted into the population Q t . Otherwise, the algorithm destroys the partial constructed solution in S p and generates a new random key R a (to sort R a ).
The entire constructive procedure finishes when all of the customers have been assigned into S p or after having a successive number of attempts without producing a feasible solution. When a feasible assignment is reached, the set S represents a feasible initial solution of routes. The customers already inserted in S are removed from S a . Figure 5 shows an illustrative example of a feasible assignment. After reaching a feasible assignment, the sequencing mechanism is applied to construct each route by respecting the order of insertion and, based on this, the corresponding calculations of the objective functions L (representing the total latency of the system) and T (total tardiness of the system, based on the priorities of the customers) are performed. Algorithm 2 shows the pseudocode for this algorithm.
Create an auxiliary random key chain R a with values [0, 1); Sort customers in S a in a non-decreasing order according with their corresponding random values in R a ; if The insertion of the jth customer is feasible to insert in route l then Insert the jth customer in the lth route S p ; Remove the assigned customer from S a ; j + +; f lag = 0; end else if l < K then l + +; end else Destroy the partial solution, S p ← ∅; Establish S a = {1, 2, · · · , n}; Create a new random key R a ; f lag = 1 ; end end Compute the values of total latency L and total tardiness T for the individual; end end

Crossover Procedure with Local Search Strategies
The proposed crossover procedure consists of the combination of two solutions, A and B, to create a new solution C. A tournament selection operator is incorporated to diversify the creation of new solutions. After obtaining the new individual C, a selective local search mechanism can be applied to improve it.
The procedure receives the following inputs: the current population (Q), the updated population (R), and the number of children to create (N). Then, the mechanism starts by selecting two individuals from the current generation (P t ). The first individual will belong to the front F 0 , whereas the second one will be chosen at random from the entire population (generation). Subsequently, the customers of the routes for each solution are grouped into a single big chain by following the assignment order (starting from the first customer belonging to the first route and ending at the last customer of the latter one). As a result, the chains for the corresponding chromosomes A and B are obtained.
The creation of the new individual C is based on the information of the random keys (R a ) of each parent (chromosome). Then, a probability for inheriting is assigned to each parent. These probabilities are complementary. In other words, if the probability of P A = α is assigned for the first chromosome (A), then the second chromosomes will receive a probability of P B = 1 − α. Then, for every position to fill in the R a belonging to the customer, the roulette wheel is spun (RN) to determine if the element belonging to the R a for the first or second parent must be selected to insert in the child. If the value of RN ≤ P A , the ith element R a is included in the child. Otherwise, the element belonging to the second parent is selected. The mechanism stops when all of the positions have been evaluated. Figure 6 illustrates and example of this mechanism. To enhance the creation of reasonable quality solutions, the probability assigned to the parent belonging to F 0 is always greater than 50%. Since the decoding procedure is a simple mechanism, it might occur that different random keys lead to an identical solution. Once the Ra for a child solution C is obtained, its feasibility is evaluated by calling back the constructive procedure. If the resulting solution is feasible, then the total latency L C and tardiness T C objectives are computed. On the contrary, the child is discarded, and a new second parent is selected (preserving the first individual) from P t to restart the crossover procedure. Algorithm 3 depicts the pseudocode of the crossover mechanism. The LS procedure is based on local search strategies, applied to intensify the search in pursuit of finding local minima. This procedure consists of five different neighborhood structures arranged into two classes, namely intra-route and inter-route mechanisms, performing them iteratively. This strategy has proved to be successful for a mono-objective version of the CCVRP [10]. Below, we describe each type of move: • Intra-route swap. The procedure exchanges the positions of two customers belonging to the same route. For instance, if the customers to exchange belong to positions h and i, then arcs (h − 1, h), (h, h + 1), (i − 1, i) and (i, i + 1) are removed and replaced by arcs (h − 1, i), (i, h + 1), (i − 1, h) and (h, i + 1). It is important to remark that these movements do not affect feasibility in terms of capacity. • Intra-route reallocation. This mechanism deletes a customer from its current position and reinserts it into another position on the same route. • Inter-routes interchange. This strategy exchanges two customers belonging to different routes, as long as the move keeps feasibility (in terms of capacity). • Inter-routes reallocation. For a given customer, the operator searches for the best position of the customer to move in any of the routes. If the best-identified position is different from the current one, the movement is performed.
The two major strategies operate as follows: At first, the initial solutions are sent to the intra-local search procedure, where the intra-local search strategies are applied. Then, the local minimum is forwarded to perform inter-routes local search strategies. These procedures are iteratively implemented, while the current solution value L keeps improving. In each process, the reallocation movement is performed first, and the execution of the interchange movement next. The first improvement criterion (FI) is used. Algorithm 4 exhibits this process. As observed, this mutation procedure seeks to insert improved individuals to the next generation, although the mechanism does not guarantee that the chromosome selected can be deeply improved. Because the size of this subset is relatively small, it is always possible to find a chromosome to be improved.
This procedure differentiates the two versions of the MA. For the first version (MA-RK v1), all of the feasible individuals generated by the crossover mechanism are sent to the LS procedure (threshold = 1). For the second version (MA-RK v2), the local search mechanism is applied only to a certain percentage of individuals, expecting to accelerate the performance of the algorithm in terms of CPU time.

Computational Results
This section is devoted to reporting the computational experiments conducted to assess the efficiency of the proposed approach. First, we provide the set of instances used to perform the tests, as well as the characteristics of the computational equipment used. Secondly, we present the parameter setting for our versions of the MA. Finally, the experimental results for both the mathematical formulation and the metaheuristic procedure are displayed, accompanied by the respective discussion.

Test Instances
The instances used to conduct the experimentation were adapted from the ones used in the literature to evaluate the multi depot VRP with heterogeneous fleet: Koulaeian et al. [41] (Kou15), Chunyu and Xiaobo [42] (CaX10), Gillett and Jhonson [43] (GaJ76-7-GaJ76-12), and Augerat et al. [44] (Pn16k8 and Pn23k8). Even though there are some instances proposed by Talliard [45] and Li et al. [46] for the classical Heterogeneous Fleet Vehicle Routing Problem, we decided to use these instances since some of them provide a reasonable size in the number of customers for assessing the model.
The size of the instances ranges from 12 to 100 nodes and from 2 to 10 vehicles. The generated problem instances are characterized by the following criteria: (i) number of customers; (ii) number of vehicles; (iii) coordinates (x,y) for all locations (including the depot); (iv) demand of each node; and (v) priority index for each node. Since the original instances consider multiple depots and do not consider any preference index between the customers to be served, we selected the first depot as the single-origin, and we included a priority parameter by assigning a numerical index within {1, √ n }, based on a uniform distribution. The customers having the highest value of the index represent the ones that should be first served (highest level of priority).
Since the modified instances are set to deal with a different problem, and to facilitate the report of the results, we decided to rename them using the nomenclature "FNO-x", where x denotes a consecutive number assigned according to the rank assigned to the instance (in terms of the number of nodes, following a non-increasing order). For example, the instance Kou15 is the one with the lowest number of customers (12); therefore, it was renamed as FNO1. The remaining instances based on Pn16k8, CaX10, Pn23k8, GaJ76-7, GaJ76-8, GaJ76-9, GaJ76-11, and GaJ76-12 were renamed as FNO2, FNO3, FNO4, FNO5, FNO6, FNO7, FNO8, and FNO9 respectively. In the case of the instance GaJ76-10, it was named FNO10 because it has the largest size in the number of customers. These new instances are available by request.
All of the experiments were conducted using a PC Intel Core i7 @2.30 GHz with 16 GB of RAM Memory under Windows 10. The formulation was modeled using AMPL and solved using Gurobi 9.0. For each instance, we established a time limit of 7200 s (2 h). In the case of the MA-RK, both versions were coded using the C++ language. In the next subsection, the parameter tuning is presented.

Parameters Setting
In the case of the MA-RK v1, the values for the parameters corresponding to the size of the population N, the threshold value β, and the maximum number of generations D were set as 1000, 0.1, and 100, respectively. In the case of MA-RK v2, we set a threshold of 0.4. These values were obtained after performing a preliminary analysis over a subset of instances randomly selected. In addition, several tests with different number of iterations were conducted, finding that, for all the analyzed instances, the MA-RK stops improving when it reaches 80% of the maximum number of iterations, depending on the instance. Lastly, to evaluate consistency, each instance was executed 10 times, and the best front obtained is reported. The next sections show the numerical results computed over the test instances.

Experimental Results
The first set of experiments aims at evaluating the performance of the formulation concerning optimality (optimally solved instances), and the effectiveness of the MA-RK comparing the results with those obtained by the resolution of the model. We first present the results for Instances FNO1 and FNO2 (up to 15 nodes).
The following metrics were used to compare the performance of the exact and approximation procedures: • Number of points on the front (the larger, the better) • CPU time (in seconds) (the shorter, the better) • k-distance [47] (the smaller, the better) • Hypervolume [48] (the larger, the better) • The coverage of the fronts [48] computed of one front over another, denoted by c(X', X") (the higher, the better) These metrics have shown their successful implementation in biobjective VRPs [19,39]. The number of non-dominated points measures the ability of each method to find efficient fronts. Table 1 summarizes these results for Instances FNO1 and FNO2. Figures 7 and 8        A point to highlight in Figures 7 and 8 and Table 1 is that, for both instances, the Pareto front obtained by Gurobi is densely crowded. Additionally, notice that both versions of MA-RK performed differently over the solved instances. In particular, for Instance FNO1, the second version of the algorithm produced a front that is closer to the optimal front obtained by CPLEX. On the contrary, for Instance FNO2, both algorithms produced fronts near to the optimal front and, in particular, the MA-RK v1 produced a more dense front than the MA-RK v2. Regarding the computational time, the summarized results are presented in Table 2. In this case, the exact method required around 1 h for solving the FNO1 instance (12 nodes), whereas the time required to solve the FNO2 (15 nodes) instance was almost four times as long. In particular, both versions of the metaheuristic required less than 1 s to obtain the solutions. This fact, in conjunction with the metric of the quantity of non-dominated points, supports the evidence that, in general, the MA-RK algorithm performed very well. Regarding the density of the fronts, Table 3 shows the average k-distance value of all points on the efficient frontiers for each instance, with k = 3. Specifically, the MA-RK produced fronts with more density than AUGMECON2. In particular, for the FNO2 instance, it can be seen that the NSGA-v1 obtained the minimum values for the maximum and average distances, while the MA-RK v2 obtained a denser front for Instance FNO1. To verify the efficiency of the MA-RK (in both versions), we used the hypervolume metric. Table 4 displays the obtained results. Again, the MA-RK v2 provided a better value of hypervolume for Instance FNO1, while the MA-RK v1 performed better on Instance FNO2. Finally, we used the coverage measure (considering only the strict domination). Tables 5 and 6 exhibit the results. In these tables, a value of C(X', X") equal to 1 means that all points in the estimated efficient frontier X" are strictly dominated by points in the estimated efficient frontier X'. As expected, the exact method dominates both algorithms entirely in terms of the space covered. Regarding the metaheuristic procedures, for Instance FNO1, the MA-RK v2 was able to generate points that dominate almost 77.78% over the ones generated by the MA-RK v1. For Instance FNO2, the front of the MA-RK v2 dominates 33.33% of the points generated by the MA-RK v1 which, in turn, dominates 25% of the points generated by the MA-RK v2. We observe that, for Instances FNO1 and FNO2, when the minimum value of latency is obtained, the maximum amount of total tardiness rises to 3.2 times the cost of latency, which might translate to a higher level of customer dissatisfaction. On the contrary, when the minimum value of total tardiness is reached, the overall latency of the system rises up to 1.4 times the optimal (minimum) value. In this case, the increment of latency translates into a significant increase in the cost and, therefore, in a reduction of profit. In addition, it is obvious that the prioritization of the customers generates an unbalance in their demand, especially for the case of having routes where relatively few customers can have significantly high amounts of demand compared to the rest.
Another aspect to highlight is that the decision-making process can be seen from two perspectives: (1) savings in tardiness costs can represent up to 72% of the total costs; and (2) savings in latency produce savings for up to 28% of the total costs. In other words, according to the objective function, if a preference must be defined a priori, tardiness must be more important than latency.

Experimental Results for Larger Instances
The experimentation considering large size instances was conducted in both versions of the algorithm. The complementary experimentation involves instances of up to 100 nodes. Tables 7-12 display the results of our computational experimentation.
In Table 7, Column 1 displays the name of the instance. Columns 2 and 3 indicate the size in terms of the number of nodes and the number of routes. Columns 4 and 5 report the number of non-dominated solutions obtained by each algorithm. For Tables 9 and 10, Column 1 shows the name of the instance, whereas Columns 2 and 3 report the value of the corresponding algorithms over the evaluated metric. Specifically, for Table 11, two columns are used to indicate the maximum and average k-distances for each procedure.
Following the same sequence used in the previous section, the first metric to compare is the quantity of nondominated points. Table 7 reports the number of non-dominated points obtained by each algorithm (Pareto front). According to the information there, MA-RK v2 was able to obtain a higher quantity of non-dominated points. In some instances, the number of points reported almost doubled the amount of the ones obtained by MA-RK v1.This can be explained by the fact that MA-RK v2 generates more diverse solutions, since the process of intensification is selective. In Table 8, the minimum and maximum values for each objective are shown. According to the information obtained, for most of the instances, the MA-RK v2 reports better values for the total latency. In addition, the MA-RK v2 produces better values of tardiness for most of the instances. In summary, the selective version of the MA-RK clearly outperforms the MA-RK v1.      Due to this, the metric of the execution time was evaluated to verify if any of the versions performs more quickly. Table 9 displays the elapsed CPU time for the best execution. As expected, the required time increases as the size of the instances increased. However, both versions of the metaheuristic were working within a similar computational performance range.
As for the third metric, the hypervolume, the results of the algorithm are displayed in Table 10. There, we do not have enough evidence to confirm that any algorithm outperforms the other. What can be confirmed is that MA-RK v2 produced higher values of hypervolume for seven out of eight instances. However, for the instance where the MA-RK v1 obtained better values, the difference against the NSGA was small. Regarding the density of the frontiers, Table 11 shows the results obtained. From these results, it can be noticed that, in most of the cases, the MA-RK v2 produced lower values for the maximum distances (more compactness). However, the MA-RK v1 algorithm produced better values for the average distances. Finally, Table 12 reports the values obtained for the set coverage metric. The first column refers to the name of the instance, and the rest show comparisons in coverage between the algorithms. Again, MA-RK v2 performed better than MA-RK v1, by dominating the entire front provided by MA-RK v1. This confirms that the selective version of the MA-RK clearly dominates.  In summary, we conclude that, even when both metaheuristics provide good results in a reasonable computational time, the MA-RK v2 consistently outperforms the non-selective MA-RK version.

Conclusions and Future Work
This study addressed the biobjective Cumulative Capacitated Vehicle Routing Problem. This problem mainly arises in commercial contexts such as the delivery of perishable goods, in which there are differentiated based on priority indexes. In the case of a pooled transportation service, it might help to estimate the trade-off between delivering the orders in the same sequence as customers board the vehicle and the minimum arrival time of the system. For this problem, a mixed-integer programming formulation and two metaheuristic algorithms were developed. A commercial optimization software was able to solve the model for small size instances, whereas the algorithms showed their effectiveness by providing feasible results in a reasonable amount of computational time.
The algorithms showed their efficiency by providing good quality fronts for the small size instances. Additionally, for larger instances, both algorithms provided good values for the multiobjective metrics evaluated. Although none of the algorithms outperformed each other, the MA-RK v2 obtained fronts with a higher quantity of points, more density, and more coverage of the sets. However, the MA-RK v1 was slightly faster. One intuition is that MA-RK v1 is more intensive in the local search, stagnating in local optima, while MA-RK v2 maintains the diversity, allowing to escape from local optima and populating the Pareto-fronts. However, more computational experiments are needed to clarify this effect.
In summary, all procedures provided a positive contribution to a more sustainable balance between economic and customer service objectives. Our results provide useful insights for business applications in terms of considering customer satisfaction and gaining a valuable sustainable advantage given that the reduction in the traveled time translates into a reduction of CO 2 emissions. Future research lines include the design of routes using congested environments with travel speed variation. This fact can be addressed either by modifying the objective function to include the variability in the travel time during the day or using a risk aversion approach, by adding a profit to each node associated with the order of visit. In addition, considerations involving time windows as priority metrics, and demand uncertainty can be worth consideration, as well as factors such as labor costs or balance of the total traveled distance among routes, which seem to dominate the overall cost.