A Decision Support Tool for Urban Freight Transport Planning Based on a Multi-Objective Evolutionary Algorithm

We present an optimization procedure based on a hybrid version of an evolutionary multi-objective decision-making algorithm for its application in urban freight transportation planning problems. This tool is intended to solve the planning problems of a merchandise distribution firm that dispatches small volume fractional loads of fresh foods on daily schedules. The firm owns a network of distribution centers supplying a large number of small businesses in Buenos Aires and its surroundings. The recombination operator of the evolutionary algorithm used here has been designed specifically for this problem. It is intended to embody a strategy that takes into account constraints like temporary closeness, closeness time window and connectivity in order to improve its performance in the clustering phase. The representation allows incorporating specific information about the actual instances of the problem and uses adaptive control of the parameters in the calibration stage. The performance of the proposed optimizer was tested against the results obtained by two evolutionary algorithms, NSGA II and SPEA 2, widely used in similar problems. We use hypervolume as a measure of convergence and dispersion of Pareto fronts. The statistical analysis of the results obtained with the three algorithms uses the Wilcoxon rank sum test, which yields evidence that our procedure provides good results.


I. INTRODUCTION
Decision-making tools based on bio-inspired algorithms have been successfully used in logistics during the last decades. They have been continuously improved in the context of urban freight transport (UFT). The goal has always been increasing the efficiency and competitiveness of the firms, an objective usually hampered by the atomization of the sector and the complexity of logistic management at this stage of supply chains. A frequent issue involves taking into account in the decision-making process the needs of third parties since externalities over the relations with other agents may lead to quality and competitiveness losses in merchandise deliverance.
The associate editor coordinating the review of this manuscript and approving it for publication was Ching-Ter Chang .
We seek here to overcome those limitations by changing to a multi-objective cooperative objective approach, taking into account the interests of all the parties involved in the process, ranging from managers of distribution centers to the final customers. We proceed by developing a hybrid version of an evolutionary multi-objective algorithm addressing the problem of a firm delivering perishable fresh goods from several distribution centers, carrying relatively small fractional volumes to a large number of grocery stores in Buenos Aires and its satellite counties.

II. LITERATURE REVIEW
Herbert A. Simon pioneered the view of decision-making as an iterative process in which rationality is bounded by the inherent features of the decision-maker and the context in which the decision has to be made [1]. Simon insisted that this process can be enhanced with the help of computational tools, providing rational information to human decision-makers, that will be able to interpret it in the light of their knowledge and beliefs [2]. This does not mean that the entire process should become automatized. On the contrary, the computational contribution consists in providing data and computer processing to help in the decision process [3], being human judgement the source of the final decision.
The study of decision-making in UFT indicates that agents make decisions evaluating simultaneously different objectives, usually conflicting ones [4]. Salazar, Carrasquero and Galván [5], took Simon's insight, developing a decision-aid tool for this framework. The ensuing process proceeds by stages.
The first stage involves the construction of the decisionmaking model. That is, to state analytically a multi-objective optimization model and to choose a solution technique. A second stage involves the actual process of searching solutions. The human decision maker intervenes here by expressing her preferences with respect to the alternative solutions (see [6]- [12] and [13]). Finally, in the last stage, the decisionmaker chooses the final solution from among several alternatives.
The initial characterization of the decision problem requires revising the different approaches to the UFT vehicle routing and truck loading problems. These approaches can be classified according to the characteristics of the distribution network. Miguel, Frutos and Tohmé [14] present such taxonomy. Considering that classification, we find that our problem of interest has features proper of certain variants of the Vehicle Routing Problem (VRP), namely the Capacitated Vehicle Routing Problem with Time Windows (CVRP-TW), the Multi Depot Vehicle Routing Problem (MD-CVRP) and the Time Dependent Vehicle Routing Problem (TD-CVRP). Figure 1 highlights the main components of our own problem (in the gray box).
The Capacitated Vehicle Routing Problem (CVRP) starts by considering the existence of a certain number of clients to be supplied with a given volume of merchandise from a single depot. To carry out the distribution there are a certain number of vehicles with given capacities. Each vehicle visits exactly once each of its assigned clients and each client is visited by a single vehicle. The sum of the demands of the clients assigned to each vehicle should not exceed its capacity. The objective is to determine the sequence of deliveries minimizing the total cost of distribution, which is assumed to be proportional to the distances traveled [15]- [22]. The Multi Depot Vehicle Routing Problem (MD-CVRP), is a variant of the CVRP which assumes multiple depots or storage sites with different locations [23]- [25]. The Vehicle Routing Problem with Time Windows (CVRP-TW), is the variant that assumes the existence of allowable time intervals for the deliverance to each client. The CVRP-TW, in turn, can adopt different variants, since time windows may exit for different elements of the service, be it the time of arrival at the clients' stores, the working time of drivers, the hours of activity of depots, etc. [26]- [33]. Finally, the Time Dependent Vehicle Routing Problem (TD-CVRP), takes into account the variations of traffic density in different areas, causing fluctuations in the speed of the vehicles, independently of the distances traveled. Because of this, we focus on the minimization of traveling times and not the distance traveled by the vehicles [34].
Gayialis, Konstantakopoulos and Tatsiopoulos in a review of the literature on the Vehicle Routing Problem (VRP) [35], found that 92% of the contributions take into account capacity constraints (C-VRP); 39% of the works analyzed consider the time windows for delivery (VRP-TW); 12% minimizes the travel distance (TD-VRP) while 18% analyzes deliveries up from different centers (MD-VRP). According to [36] and [37], from all the studies of UFT problems with multiple storage sites, only 10% consider multiple objectives, and from these only a few assume simultaneously time windows and capacity constraints (MD-VRP-TW), but without considering temporal dependencies [38] and [39]. On the other hand, few articles consider simultaneously time dependences, time windows and capacity constraints (TD-CVRP-TW) [40] but without assuming multiple storage constraints. We compose these different formulations in a single framework. We call the resulting overall problem to be addressed in this paper the Multi-depot Time Dependent Capacitated Vehicle Routing Problem with Time Windows (MD-TDCVRP-TW).
With respect to the optimization criteria, according again Vega-Mejía et al. [37], among the multi-objective approaches, 52% minimize the total distance traveled, while the minimization of travelling time is an objective only in 13% of those contributions. Other objectives, as pointed out by those and other reviewers, are the minimization of the number of vehicles, as a way of reducing the sub-use of the vehicles [32] and [39]; the load balance in the vehicles, be it respect to the total time on each route, the amount of merchandise on each route, or the number of clients served by route [39]. This goal is particularly interesting from the point of view of the use of human resources [41]. Other objectives in the literature are the minimization of risks, present in the delivery of hazardous materials [23]; the minimization of CO 2 emissions [40] or the maximization of client satisfaction, derived from the satisfaction of the time constraints posed by them to minimize tardiness or earliness.
With respect to the selection of the solution method, since the problem is NP-hard ( [42]- [44]), exact methods can only solve relatively small instances of the aforementioned types. It is natural then, that most of the methods applied in the literature are based on the use of meta-heuristics (80% of the contributions in the last decades), and pure heuristics (in the order of 15%), while all other methods are used only in 5% of the contributions.
The meta-heuristics that have been applied in the UFT context are ant colony optimization [30], artificial bee colony [34], biogeography-based algorithms [23], the firefly algorithm [19], evolutionary algorithms [15], [17], [18], [24], [25], [31], [21], taboo search [33], [38] and simulated annealing [29], [40], among others. From the articles that address multi-objective versions of VRP with more than 20 nodes (clients), in more than 80% of them use, as said, meta-heuristics, and among them, in particular Evolutionary Multi-Objective Algorithms (MOEAs) [36]. The successes in the use of MOEAs for the treatment of problems related to ours lead us to adopt this approach for our analysis. A summary of CVRP literature and its variants, relevant for our analysis, is presented in Table XIV (in the Appendix).

III. CHARACTERIZATION OF THE PROBLEM
As indicated previously, we focus on the case faced by the head of logistics of a fresh foods distribution company that works on daily schedules, delivering its merchandise with a fleet of trucks to wholesale sellers in four regional marketplaces, which have contracts with two hundred grocery stores in the Buenos Aires metropolitan area (GBA). Tables 10 and 11 of the Appendix present detailed information about the distribution centers and the clients.
These goods are distributed in urban settings. Several of the delivery routes are congested during working hours, and thus the average speed of the vehicles on them are lower than in other routes. We consider then t ij , the time that takes to go from site i site j on a given route, independently of the distance traversed.
The load unit consists of a homogeneous cardboard box weighting 20 kg, with a length of 50 cm, width 32 cm and height 30 cm. These boxes can be piled up, ensuring a compact use of space in the truck. The vehicle fleet consists of homogeneous trucks with capacity and speed that agree with normative rulings for urban transportation. We assume here that each of them can load 8.000 kg (400 load units).
Each client prepares a request detailing the number of load units and the time window at which it is able to receive it. This means that no delays can be admitted and, if a truck arrives earlier it has to wait to carry out the deliverance. The service times at the facilities of the clients involves the times from arrival to departure, which depend in particular from the intermediate activities of parking, unloading, signing documentation, etc.
With respect to the objectives, we consider the agent in charge of designing the distribution plan seeks to minimize the distribution times (a proxy for the costs in urban contexts) complying with the requirements of the clients; another goal is to balance the workload among the different vehicles. In

summary, the Multi-objective Multi-Depot Time Dependent Capacitated Vehicle Routing Problem with Time Windows (MO-MDTDCVRP-TW), can be defined as follows:
Find the Pareto-optimal solutions (i.e. the Pareto front) obtained by simultaneously minimizing the total time of distribution and balancing the workload of the different vehicles, while satisfying the requests of many clients at the committed delivery times, departing from several distribution centers, using a fleet of vehicles with homogeneous load capacities and restricted working hours, taking into account the traffic flow density on the respective routes. Upon the determination of the Pareto front, the decisionmaker picks one solution on it.

IV. THE MODEL
The MO-MDTDCVRP-TW can be represented as a graph G = (V , E), where V is the set of nodes and E the set of edges. V can be partitioned in two subsets, that of storing sites and that of retailers, Each client i ∈ V R has a demand d i , a service time ts i > 0 and time window [o i , c i ]. These windows are of the ''hard'' type, meaning that if the vehicle reaches i before o i , it will not be received and will have to wait until o i . In turn, if it reaches node i after time c i , it will be unable to deliver the request, making the program unfeasible. Besides, to ensure feasibility, vehicle s cannot exceed r s , defined as the maximal time of operation allowed in a day for that truck.
For each deposit l ∈ V L we assume ts l = 0, implying that the vehicle is already loaded at the start of the program. The fleet is S l = {1, 2, . . . , K l , consisting of K l vehicles. Each s ∈ S l has a load capacity of Q = 400 load units.
Each edge (i, j) ∈ E, has assigned a time tr ij , including t ij , the time required to go from node i to node j on a given route, plus the service time at the destination node, ts j .

A. BINARY VARIABLES
We use the following binary variables: x sl ij , which equals 1 if vehicle s goes from i to j, departing from storage site l and 0, otherwise.

C. OBJECTIVE FUNCTIONS
As discussed before, the problem consists of the simultaneous minimization of two functions: Function f 1 represents the total time spent on all routes as well as the costs induced by early arrivals (i.e. waiting times): Function f 2 , represents the standard deviation of work times between vehicles. Its minimization generates balanced workloads on all routes.
These constraints indicate that the number of vehicles used on a route cannot exceed the size of the fleet at the depot at its origin (|S l | = K l , ∀l ∈ V L ).
The following constraints indicate that the total load cannot exceed the capacity of each vehicle: The following ensure that each vehicle starts and ends at the same depot: These constraints preserve the flow. That is, if a vehicle s reaches client r, it has to depart next from there: A vehicle cannot go from a depot to another: (7) indicates that the load on a truck traversing edge (i, j) should not exceed the capacity of the vehicle: 0 ≤ w sl ij ≤ Q · x sl ij ∀i, j ∈ V , ∀s ∈ S, ∀l ∈ V L (7) (7) In turn (8) indicates that the merchandise unloaded at r must be equal to the demand of that client: These constraints establish that the total time spent by a vehicle on a route cannot exceed the total time allowed for the route: The following constraints that delays at the depots cannot be allowed: (11) indicates that if a client j is served by s starting of a depot l, after serving client i; then the arrival time at j must be later than the departure time from i: The following constraints ensure that the deliveries verify the time windows: Finally, the following define the range of values of the variables:

V. SOLUTION METHOD
The algorithm we use to run the optimization process is a hybrid variant of the elitist non-dominated sorting genetic algorithm (NSGA-II), developed by Deb, Agrawal, Pratap and Meyarivan [45]. NSGA-II is flexible enough to admit a representation and a crossover operation specifically defined for our problem. The knowledge of the decision maker contributes to guide the algorithm towards the best solutions. Here, in particular, we replace the Pareto dominance originally assumed in NSGA-II by the criterion of g-dominance [6]. We call our variant, HY-NSGA-II. Table 1 presents its pseudocode.

A. REPRESENTATION
We apply a path representation, based on the permutation of integers with a single chromosome consisting of four genomes for the calibration stage and three for the optimization stage. The first genome contains the sequence of clients; genomes 2 and 3 contain information about genome 1. In turn, genome 4 contains information about the parameters to be evolved in the calibration stage.
We present, as an example, the case of an individual representing the potential solution of an instance with 20 clients, In this example G1 has twenty places, each one corresponding to one of the clients; each permutation of G1 is a sequence of visits. G2 has twelve entries, indicating where the routes are to be distinguished in G1 for each depot (we assume here that each depot has the same amount of potential routes assigned to it, K ). So, for instance, in G2 the first route of depot 3 ends at 15 while the second one ends at 19, meaning that the second route from depot 3 covers the clients at positions 16 to 19 of genome 1, i.e. 13,9,20,10 . Figure 2 shows the scheme of the chromosome. The zeros of G2 represent empty routes (non-contracted services) at each depot. The three places in G3 indicate the amount of active routes per depot. Finally G4, represents the parameters to be calibrated.

B. RECOMBINATION OPERATOR
We use an operator which we call ERX-MD, inspired in the edge recombination operator of Whitley, Starkweather and Fuquay [46]. The offspring is obtained by combining edges present in the chromosomes of the parents. This method has the better performances on representations based on integer permutations [47].
The improvements introduced in the progeny stems from the feasible edges of the parents, obtained according to the following steps: 1st. Feasible edge map: a list of feasible edges for each cluster corresponding to each depot is obtained. 2nd. Reclustering: the nodes in different clusters are regrouped, seeking to maximize the forward connections in each new cluster. 3rd Tour construction: the progeny is obtained from feasible edge map and the forward connections, using a variant of the insertion heuristic of [48] as follows: A node j becomes a child of i if it yields the lowest of the costs ϑ obtained by weighing these three measures: tr ij : Temporary Closeness: it is obtained as the average speed of traveling the edge between nodes i and j. δ j : Connectivity: computes the number of forward nodes from j. ϕ j : Closeness Time Window: gives priority to a node j with the closest time window, i.e. the urgency to insert it in a route.
The next node is obtained by comparing the values ϑ ij = α 1 · tr ij + α 2 · δ j + α 3 · ϕ j , of the candidate j nodes connected to i. Figure 3 presents the schematics of this procedure.    Figure 4 show a reclustering of two parent with ten nodes each and two clusters (depot 1 in dark gray and depot 2 in light gray), assuming a single route per cluster.
In the feasible edge map of each cluster we can see that node 6 belongs to both clusters. In the reclustering, node 6 has to be assigned to cluster 2, given that it has a higher degree of connectivity than in cluster 1. If, instead, node 6 were assigned to cluster 1, it would become isolated after node 5 is assigned to cluster 2 (see Figure 5).  . Sub-graphs of clusters 1 and 2 according to how node 6 is assigned. Figure 6 illustrates the graphs of the maps of the feasible edges, assigning node 6 to clusters 1 or 2. In this figure we can see that by assigning node 6 to cluster 1 it reduces the feasible edges in cluster 2, leaving it only with two feasible edges.
In summary, ERX-MD assigns first clients to depots and then it determines the routes, defining how the clients of each depot have to be visited. That is, it takes the information of the parents and regroups the over-assigned nodes to improve the connectivity in the subgraph of each cluster. Afterwards, the routes are determined without any further corrections.
This strategy solves, in a relatively cheap way, the problem of determining an offspring satisfying the constraints of this kind of problems. This is particularly relevant, considering the usual costs of eliminating unfeasible individuals.
With respect to the mutation operator, the relevant question is to introduce diversity by adapting a standard insertion operator. A gene from G1 is inserted in a feasible position of the same genome of a single individual, chosen at random. The change can have either inter-route, intra-route or interdepot effects.

C. INTRODUCING PARTIAL PREFERENCES (G-DOMINANCE)
In the context of our problem, the decision-maker 1 knows approximately well the zone of the objective space that is desirable according to the goals of all the parties. This information is used to guide the search towards feasible and efficient solutions. We use the concept of g-dominance to incorporate preferences up from a reference point g that determines the alternatives that are more or less preferred than g [6].   We show an initial population that evolves towards the region of interest of the decision maker, denoted FP * . So, instead of constructing the entire Pareto front, this method quickly leads to the region of interest.
To implement g-dominance penalties are applied on solutions outside the preferred regions. The corresponding pseudo-code is as follows: Therefore, g-dominance penalizes the solutions that are dominated by the given reference point g. In other words, Figure 8 describes the process of decision making followed daily. It starts by the consolidation of the data compiled from different sources, first with respect to origins or destinations, as well those generated in the process of Enterprise Resource Planning(ERP) as for instance the inventory available in the depots, the amounts of merchandise to deliver to each client, the average stop delivering time or the time windows delivery constraints, etc. Secondly, data provided by the Geographic Information System (GIS) on the road network, or traffictimes is also consolidated. On the basis of this information, the g-dominance procedure is applied once a reference point is defined.

D. THE DECISION PROCESS
This instance of a decision process is then solved by the application of our proposed HY-NSGA-II, yielding the portion of interest in the Pareto front (FP * ). By visualizing the resulting solutions the decision maker will select an alternative, upon which she will generate departures, route assignments, time tables, metrics of performance, amount of vehicles used, average speed on the road, etc. If this resulting information is not satisfactory, the decision maker can select another solution from the portion of interest or even introduce a new reference point to start from scratch the search. This can be repeated until a solution acceptable for the decision maker is found.
The reports generated by the chosen solution are translated into the documents necessary for the execution of the distribution schedule (e.g. the roadmap for the drivers) and for the accounting and traceability activities associated to the delivery of goods.

VI. COMPUTER EXPERIMENTS
Given the features of the problem analyzed here, we lack of published benchmarks to which to compare the pros and cons of our algorithm. Therefore, we proceed to validate and test its performance on a real world instance, comparing the results to those obtained from the application of Deb, Agrawal, Pratap and Meyarivan's [45] version of NSGAII without our variations and Zitzler, Laumanns and Thiele's SPEA2 [49].
In order to generate significant comparisons we incorporate in the benchmarks the g-dominance strategy as well as the same mutation operator. We use a quantitative measure of both the convergence to the Pareto front and the distribution of solutions on that front, namely the Hypervolume Indicator (or S metric) [50]. This choice facilitates the evaluation of the relative performance of multi-objective optimizers, in particular when the actual Pareto front of the problem is unknown [50].   At the calibration stage we used a self-adapting strategy of parameter control, according to which the parameters evolve with the individuals, being codified as part of the chromosomes (see the Representation subsection) and subject to the same rules of variation and selection. The strategy of selecting successful parameters takes simultaneously into account the results from the three optimizers, as to avoid biases in favor of any of them. We obtained thus optimal parameters to run the experiments, with a crossing probability of 0.8, a probability of mutation of 0.01, a maximum of 1000 generations and a population of 500 individuals.

A. GRAPHICAL REPRESENTATION OF THE SIMULATIONS
We start by running a simulation of the three algorithms, assuming that the regions of interest are at the extremes and the middle of the objectives space. We took the following as the reference points for g-dominance: G1 = (290, 250), G2 = (300, 180) and G3 = (345, 145). Figure 9 presents the results for HY-NSGA-II, NSGA-II and SPEA2.
Not knowing the actual Pareto front makes it impossible to compare the simulations with the actual front. But it is still possible to compare them with those obtained with the benchmark algorithms. In this sense, HY-NSGA-II seems to perform acceptably well.
We distinguished in the graphs the alternatives chosen by the decision maker. In Figure 10 we highlight those solutions (A, B and C, respectively) in the space of decisions. We transformed the geographical coordinates of the nodes representing the clients into Cartesian ones, as x ∈ [0, 1] and y ∈ [0, 1] while the edges (of a corresponding color) indicate the prescription of the chosen solution.
The next table specifies solution B, chosen for reference point G2: In this solution, between 48 and 58 clients are served per distribution center, using 4 or 5 vehicles in each route. Around a 70% of the capacity of the vehicle is used to transport the merchandise.

B. ANALYSIS OF THE RESULTS OF THE COMPUTER EXPERIMENTS
We evaluate the Hypervolume Indicator (H) for the solutions chosen by the decision maker, assuming a given reference point. We consider a point Q, with f 1 = 1.000, f 2 = 2.000, dominated by all the solutions generated by the three   algorithms. We consider three measures based on H, namely its average in 30 runs, for each algorithm, denotedH ; the standard deviations, H σ ; and the maximum value attained, H best . We obtained the following results: These results show that our algorithm yields a higher average hypervolume, achieving a better degree of convergence than the benchmarks. This means that HY-NSGA-II ensures, according to this measure, better approximations to the relevant region of the (unknown) Pareto front.
On the other hand, the variability, measured by the standard deviation is similar and sometimes worse than that of the other two algorithms.     Another relevant consequence is that our algorithm performs better for reference points that demand more equilibrium between the tasks. Tables 4, 5 and 6 indicate that the average hypervolume obtained with HY-NSGA-II is higher than that of NSGA-II and SPEA2, for reference points leaning towards the right, as G3. Figure 11 present the Boxplots capturing graphically the localization, variability and degree of asymmetry of the Hypervolume Indicator on the three reference points analyzed in Tables 4, 5 and 6.
This provides a graphical confirmation of the assessment already made on the basis of the statistical  information in Tables 4 to 6. The median of the H indicator for HY-NSGA-II is above that of the two other algorithms in all the reference points, but the improvement shows more clearly in regions requiring better balances of the loads, as in point G3.
We provide a further assessment of the algorithm, using the Wilcoxon Rank Sum Test [51]. This is a non-parametric test based on pairwise comparisons between the hypervolumes generated by the algorithms. Tables 7, 8 and 9 presents the results for the three reference points, at a significance level of 5%, with the null hypothesis: the algorithms yield the same median.
The results allow us to reject the null hypothesis and conclude that the median value of the hypervolume under HY-NSGA-II is higher than with the other algorithms.

VII. CONCLUSION
We presented a multi-objective optimization procedure for Urban Freight Transport planning. Starting from an actual problem of scheduling the delivery of merchandise to a network of clients, aligned on different routes departing from   several distribution centers, we modeled this situation as an optimization problem in which different objectives can be satisfied.
Among those objectives we considered the timing of the deliveries and the balance of loads on the vehicles. In both cases the idea is to reduce the costs involved in the distribution of goods. To solve the problem we took into consideration different constraints, ranging from the layout of the network of clients to the driving and parking regulations near the delivery points.
We called this problem MO-TDMDCVRPTW (Multiobjective Time Dependent Multi-Depot Capacitated Vehicle Routing Problem with Time Windows). Given its intractability, we chose to search for solutions applying a Multiobjective Evolutionary Algorithm (MOEA), a hybrid of the Non-dominated Sorting Genetic Algorithm (NSGA-II). We added two improvements, the first consisting in the incorporation of specific knowledge of the problem in the chromosomes of the four genomes on which run the evolutionary process. The second improvement is the use of the recombination procedure ERX-MD, which we specifically designed for this problem. It ensures that solutions satisfy the constraints without requiring costly repairs.
We also use the g-dominance strategy, which allows the introduction of the preferences of the decision maker, in order to lead the search towards regions of interest for her. The experience of the decision maker becomes then an integral part of the whole procedure.
Since the literature does not present previous results on this issue, the only way to assess the quality of our proposal is by comparison with other, more established, procedures. Using data from a real world case, we addressed it using our algorithm, HY-NSGA-II, as well as NSGA-II (in its original form, without our improvements) and SPEA2.
Our computer experiments allow us to conclude that our algorithm HY-NSGA-II is more efficient both in the convergence towards the actual Pareto front and in the distribution of solutions over the front for each of the three reference points. With respect to the variability, represented by the standard deviation of the hypervolume indicators, the results are similar to those obtained with the other two algorithms. We can see that the improvements obtained by HY-NSGA-II are larger for the reference points that correspond to more balanced schedules.
This means that HY-NSGA-II is significantly better than NSGA-II and SPEA2, which constitute, in turn, the best known algorithms for this kind of combinatorial optimization problems.
Future results involve the incorporation of further details of the real world instances of the problem. In particular, we intend to increase the number of objectives, capturing relevant aspects of decision making in logistics. Another extension involves applying other tools of computational intelligence, as fuzzy logic or neural networks, which may contribute to find even better solutions.

APPENDIX
See Tables 10-14. FERNANDO TOHMÉ received the bachelor's degree in mathematics and the Ph.D. degree in economics.
He was a former Fulbright Scholar, and held visiting positions at U.C. Berkeley, Washington University in St. Louis, and Endicott College. He is currently a Principal Researcher of CON-ICET (National Research Council of Argentina) and also a Full Professor with the Department of Economics, Universidad Nacional del Sur, Bahía Blanca, Argentina. His research has focused on decision problems, game theory, and optimization in socio-economic settings. He has published in the Theory and Decision, the Mathematics of Social Sciences, Omega, the Artificial Intelligence, the Mathematical and Computational Modeling, and the Annals of Operations Research, among others.
MÁXIMO MÉNDEZ BABEY received the M.Sc. degree in industrial engineering from the School of Industrial Engineering of Las Palmas, Spain, and the Ph.D. degree in computer sciences from the Universidad de Las Palmas de Gran Canaria (ULPGC), Spain.
He is currently an Associate Professor with the Department of Informatics and Systems, Universidad de Las Palmas de Gran Canaria, and is a member of the Instituto de Sistemas Inteligentes y Aplicaciones Numéricas en Ingeniería (SIANI) of the ULPGC as a Research Fellow. He has published more than 30 technical articles in international journals, conferences, and book chapters. His research interests include genetic algorithms, metaheuristics for multi-objective optimization, decision analysis, machine learning and real applications in mechanical engineering, renewable energies, transportation planning, and weather prediction.