Traffic congestion in large urban centers is a cause for concern because it results in economic losses, lower productivity at the global level and a negative impact on the environment. A basic logistic management issue is defining configurations of distribution to meet demand points economically, offering increasingly better qualified levels of service in terms of ability to meet demand in progressively shorter time intervals .
The Vehicle Routing Problem (VRP), proposed by Dantzig and Ramser , is of great importance in terms of effective logistic distribution . In this context, we can find in the literature a number of works that address the optimization of fuel consumption , CO2 emissions  , both of these cases simultaneously  and other approaches  . To complement the VRP, the Green VRP analyzes routing activities while taking environmental aspects into consideration   .
In recent decades, there has been an increase in the number of publications related to the VRP. At the same time, computational power has evolved in terms of performance, contributing to and stimulating the growth in research and development in this field. Even so, it remains difficult to find an exact solution to large problems, and it is often necessary to resort to heuristic and metaheuristic procedures .
Heuristic and metaheuristic algorithms are alternative procedures with regard to exact methods because of the latter’s high computational cost. Heuristic methods seek to achieve results to problems quickly, but do not guarantee that the solutions they provide are optimal solutions. On the other hand, metaheuristic methods do not usually “cling” to local optima in the search for global optimization .
Meanwhile, the field of Multi-Objective Optimization (MOO) addresses the process of optimizing two or more conflicting objectives simultaneously, subject to constraints . Thus, the objective of a MOO problem is to find non-dominated solutions and quantify the trade-offs in satisfaction between the different established objectives  .
The proposal for the Bi-objective Green Vehicle Routing Problem (BGVRP) presented here is to combine the VRP, MOO techniques and environmental considerations to generate economically and environmentally feasible routing. In this context, many researchers have adapted their solution procedures to attend to issues related to factors such as cost, travel time, CO2 emissions, fuel consumption, reverse logistics and environmental impacts   .
The aim of this work is to present a methodology to solve a BGVRP, illustrated using a newspaper distribution problem and literature instances. In this context, a multi-objective approach could aid the complex decision-making process in order to achieve the intended objectives.
The main contribution of this work lies in the methodology proposed for the BGVRP and in the approach to the case study of newspaper distribution, with a view to minimizing the distance traveled (and the consequent minimization of fuel consumption and CO2 emissions), in addition to the homogenization of demand, which is achieved by minimizing the difference between the demand for each route and the average demand considering all the routes, in terms of absolute value. Furthermore, the asymmetry of distances/real travel times is also considered here instead of the usual symmetrical Euclidean distances. Environmental considerations are included in all these aspects, treated as a mechanism for conscious decision making.
The work is structured into seven sections. Following this introduction, Section 2 contains a description of the problem to be studied, describing the methodology used. Section 3 contains correlated works, while in Section 4 the theoretical foundations for the development of the study are presented. The methodology and data treatment for the case study are addressed in Section 5. Section 6 contains the results, with the conclusion of the study presented in Section 7.
2. Characterization of the Problem
The present case study, used to illustrate the methodology adopted here to solve a BGVRP, refers to the process of distributing newspapers and magazines by a printer located in the municipality of Curitiba, in Paraná State, Brazil.
Currently, the printer delivers the newspaper to five municipalities belonging to Curitiba and its metropolitan region (Curitiba, Colombo, Almirante Tamandaré, Pinhais and São José de Pinhais), all of which are in Paraná State, in Brazil. The distribution of the newspaper is divided between the North and South regions. The average delivery is to 18,000 subscribers and 1400 copies to individual sales outlets (including newsagents, bakeries and markets).
The printer is responsible for printing and distributing newspapers and magazines for the company that hires it. The Newspaper or Newspaper with Magazine are delivered to customers at their homes. The printer receives the “artwork” for the newspaper and informed of the amount to be produced every Thursday. On Thursday, if there is time, or on Friday mornings, the newspapers and magazines are printed by 14 h and forwarded to the distribution sector. On the same day, at the distribution sector, the orders for the newspapers and magazines are organized. The process ends with the physical distribution, which begins at 21 h on Fridays, with the deliveries finalized around 8 h on Saturday mornings.
Based on the delivery list sent for the week, the distribution department is responsible for producing labels to put on all the items and dividing the sectors that, for the South region, includes the south of Curitiba and the municipalities of Pinhais and São José dos Pinhais. It should be highlighted that the “sectors” are clusters of “balanced” demand points (with approximately the same number of delivery points). This region has around 8500 subscribers and is served by 40 deliverymen.
The present study addresses only one of these 89 sectors in the South region, located in the neighborhood of Água Verde, which contains 479 addresses with 791 items, which are currently delivered by one deliveryman who uses a motorcycle. This sector is supported by a mobile depot (a van) where the newspapers are stored, and the deliveryman can “reload” his newspapers. To serve the sector in question, the deliveryman needs to reload three times, generating four routes. Total delivery time in the sector (four routes) averages 2 h, varying according to the deliveryman’s experience. It should be pointed out that each deliveryman works in around three sectors every weekend.
To simplify the database with the 479 delivery points, the printer established “48 macro points” (or fragments of streets) as a reference. This form of data treatment will be explained in more detail in Section 4.
The current solution used by the printer, i.e., the four routes of the sector, distances traveled, journey times and used capacity of the vehicles, is presented in Table 1.
Some considerations of the problem with regard to the printer should be highlighted. The employees work shifts to ensure the printing and distribution of the newspapers. There is no established operational sequence for the delivery of the newspapers. In other words, the deliveryman is responsible for the delivery sequence in his sector. There is no time window or time set by customers to receive their Newspaper or Newspaper with Magazine (every customer has agreed that delivery will be made by Saturday morning). Customers can request a new delivery if the newspaper is not delivered or is not in good condition. No limit has been set on the vehicle’s volume capacity, although the ideal limit would be up to 210 volumes. There have been situations when the deliveryman clearly loaded his motorcycle with more newspapers than it was equipped to carry. No maximum limit is set on how long it takes to cover the routes. With specific regard to the route in this study, the deliveryman is obliged to deliver only the Newspaper and Newspaper with Magazine. The company has the autonomy to use the routes in the sector as it sees fit. The journeys are mostly asymmetric because in the municipality of Curitiba there are many one-way streets.
3. Literature Review
In this section, some concepts related to the Green VRP are addressed, in addition to some correlated works.
Table 1. Routes of the printer’s current solution (newspaper distribution).
To address the theme of the Green VRP, the classification proposed by Lin et al. , was used. Therefore, the Green VRP was classified as: VRPFuel, Pollution Routing Problem (PRP) and VRP in Reverse Logistics (VRPRL).
VRPFuel addresses the optimization of energy consumption in transport and reduced fuel consumption. The PRP is used to choose a vehicle routing scheme with lower pollution levels, particularly the reduction of Greenhouse Gases (GHG). It can also include broader goals that reflect environmental costs. The VRPRL concentrates on aspects of reverse logistic distribution, which in turn is divided into four subcategories: Selective Pickups with Pricing, Waste Collection, End-of-life Goods Collection, and Simultaneous Distribution and Collection.
In the literature, we can highlight some works on the proposed theme in recent years. Those that were most instrumental in inspiring the development of the present study are described below, ordered by year of publication. Through multi-objective route based fuel consumption vehicle routing problems, Psychas, Marinaki, Marinakis, and Migdalas  used the Pareto Front (PF), applied to parallel multi-start non-dominated sorting differential evolution algorithms to minimize the cost of the route, travel time and fuel consumption. A design of a Forward/Reverse Logistics Network with Environmental Considerations was proposed by Rabbani, Saravi, and Farrokhi-Asl . The authors compared two metaheuristic procedures, NSGA-ӀӀ and MOPSO, with the former presenting the best performance.
Gupta et al. , presented a shortest path evolutionary algorithm applied to the green multi-objective multi-attribute VRP, considering time windows, pickup and delivery, heterogeneous fleets of vehicles and traffic congestions. Tests were performed in the literature and some real case studies were conducted in Singapore. Hassanzadeh and Rasti-Barzoki  presented a bi-objective mathematical model, taking into consideration production programming and vehicle routing in order to minimize the consumption of resources and delayed deliveries. The model can also reduce CO2 emissions. For this purpose, a new Genetic Algorithm was developed, the variable neighborhood search based non-dominated sorting genetic algorithm II.
With the green open location-routing problem, a mixed integer linear programming problem, Toro et al.  used the ε-constraint algorithm with a view to minimize costs, fuel consumption and CO2 emissions. Finally, the authors realized a need for other approaches to calculate fuel consumption.
Abad et al.  presented a general application for the bi-objective model for pickup and delivery PRP with integration and consolidation shipments in a crossdocking system. As a solution procedure, the non-dominated ranking genetic algorithm, NSGA-II and MOPSO were used. The algorithms were treated with a view to minimizing costs and fuel consumption. The main contribution lies in the use of cross-docking, which can be applied during loading, sorting or unloading. Gong et al.  treated the new closed-loop supply chain logistics network of the vehicle routing problem with simultaneous pickups and deliveries dominated by the remanufacturer, where they presented the bee evolutionary algorithm guiding nondominated sorting genetic algorithm II to solve it with three goals in mind: to minimize fuel consumption, minimize waiting time and minimize the distance to delivery. The problem was applied in the context of logistic distribution and generic tests were conducted. The effectiveness of the proposed algorithm was verified in comparison with the traditional NSGA-II.
Liu, Shan, Zhang, and Zhang  presented the green vehicle routing optimization problem. For this purpose, the authors used the Hybrid Quantum Immune Algorithm. Wang et al.  presented Recycling Vehicle Routing Optimization in Two-Echelon Reverse Logistics Networks, and to solve it the authors used an algorithm based on the Clarke and Wright savings method, along with the NSGA-II. In this context of reverse logistics, applied to a case in China, it was possible to minimize costs and the number of vehicles.
Erdelic and Caric  proposed a study on the Electric Vehicle Routing Problem, a variation on VRPFuel. Finally, Wang, Leng, Wang, Li, and Zhao  presented a hyper-heuristic approach to the Location-Routing Problem of Cold Chain Logistics, considering fuel consumption. The authors sought to minimize the total cost (including the cost of fuel consumption) and minimize vehicle distribution time.
4. Theoretical Foundation
This section contains a brief presentation of the main concepts and algorithms that provide the foundation for the present study.
Vehicle Routing Problem: In the literature, some mathematical approaches to solving the VRP can be found, such as those proposed by Fisher and Jaikumar  and Subramanian, Penna, Ochi and Souza . According to Lin et al. , the VRP is a problem with very wide-ranging applications in the fields of distribution and logistics management, such as newspaper distribution , and numerous other possible applications.
Clark and Wright Savings Algorithm (C&W Savings): The Clarke and Wright algorithm, proposed in 1964, is also known as the C&W Savings Heuristic or C&W Savings. This is one of the classic heuristic algorithms used for routing construction .
Tabu Search: According to Ferreira et al. , Glover  is considered the father of the Tabu Search (TS) metaheuristic due to the author’s numerous publications in this context. TS explores solutions in the solution space beyond the local optima. The method follows the premise of local search, where it generates moves at each iteration through neighborhood solutions. The move transforms one solution into another and the neighborhood is a set of solutions that enable a set of possible moves.
Multi-Objective Optimization (MOO): Is a process that optimizes two or more conflicting objectives simultaneously, subject to constraints . When addressing the problem of minimization, the mathematical model is defined as in (1) to (5), below.
where x: n-dimension decision variable vector; Z and f(x): correspond to the k-dimensional objective vector; g(x): is the set of m-dimensional inequality constraints; h(x) is the set of p-dimensional inequality constraints. In (1), the space Z = f(x) addresses the image of X, thus referred to as the feasible region for the objective function space. The constraints presented in (2) and (3) define the space of the decision variables in Rn of the feasible region X and any point of x Î X as a feasible solution. In (4) there is the decision variable domain and, finally, in (5), the objective function domain .
NSGA II Algorithm: The Genetic Algorithm (GA), developed by Holland  in the nineteen sixties, has been widely used in applications with multi-objectives. The Non-dominated Sorting Genetic Algorithm II (NSGA-II), proposed by Deb, Pratap, Aguarwal and Meyarivan , is a multi-objective GA that classifies solutions through the concept of Pareto dominance. The differential of this algorithm occurs through “Elitism”, which guarantees the preservation of good solutions in the search process, and also through the application of the Fast-non-dominated-sorting (FNS) procedure, in which the classification of the population at different levels through the use of Pareto dominance occurs.
MOPSO Algorithm: According to Abad et al. , multi-objective particle swarm optimization (MOPSO), proposed by Coello , has proved to be highly efficient and is widely used in the literature to solve optimization problems. The three steps of this algorithm are initialization, unloaded sorting and crowding distance (CD). Equations (6) and (7), below, show the speed and updating of the particle position.
where, r1 and r2 are uniform random numbers; C1 and C2 are constants that specify the acceleration of particles to pBest and gBest, and w is the inertial weight. To calculate the inertial weight, the linear distribution proposed by Kennedy, Eberhart, and Shi  can be used, as shown in Equation (8), where wmax and wmin address, respectively, the maximum and minimum amount of inertia.
Utility Theory: According to Poterba and Rotemberg , the fundamental property of Utility Theory is that the Monetary Utility Function (UF) of a Decision Maker (DM) is that it is indifferent between two alternatives if they yield the same expected utility.
Therefore, U(x) = 0 can be configured for the lowest utility value and U(x) = 1 for the highest utility value. All other utilities must lie within this interval. This procedure facilitates the definition of the relative utility of all the results from the worst to the best. Finally, when using a UF in a decision analysis process, it must be adapted to the preferences and values of the DM.
Environmental Considerations: Here, the environmental considerations related to GHG emissions are addressed. According to Carvalho , it is possible to observe the average CO2 emissions from urban transport through the different modes of transport: subway 3.16; Bus 1.28; Car 0.19; Motorcycle 0.07; Heavy vehicles 1.28. The values presented, in Carbon emissions per kilometer (kg of CO2/km), provide the input parameters (constants) to measure the amount of GHG proposed in this work.
The proposed methodology to resolve the BGVRP for the study consisted of three stages: Stage 1, Data “Treatment” of the case study (Newspaper Distribution): Macro points and Collect the times and, distances. Stage 2, Metaheuristic approaches for the BGVRP: NSGA-II, MOPSO, proposed hybrids algorithms, CWNSGA-II and CWTSNSGA-II. Stage 3, Analysis of the Results: Analysis of the Pareto front for the case study; Stratification of the opportunities obtained from the case study solution; Selection of the best route through Utility Function for the case study; Analysis of the instances selected from the literature and; Evaluation of the case study and literature instances.
5.1. Stage 1: Data “Treatment”
The R software was used along with Google Maps to collect the times and distances between the 48 macro points. The site that was used makes available some modes of transport as a reference, namely car, public transport, walking, bicycle and plane. For the present study, the “car” was used as a reference for the data collection in the software as it is closest to the reality of the activity, given that the route is covered by a motorcycle (as already mentioned in Section 2).
For the “treatment” of the case study data, the 479 delivery points (demand points) were organized in such a way as to consider fragments of the streets as a reference. Thus, the database was reduced to “48 macro points” (or “fragments of streets”) for delivery, which contains the 479 original points. These 479 points were “joined” to form the 48 macro points.
The printer provided the proposed “fragmentation” of the streets, as it uses them to support the development of the delivery itineraries of the newspapers and, therefore, made use of this proposal. Therefore, the first delivery address will be referred to as a “macro point” and the other delivery points of the respective fragments will be the delivery points belonging to the macro point in question.
Therefore, the Objective Functions (OFs) used in this case study to optimize the newspaper distribution process are, as already mentioned, “OF1: minimization of CO2 emissions” and “OF2: absolute minimization of the demands of each route in relation to the average demand”.
5.2. Stage 2: Metaheuristic Approaches
Here, a description is given of the multi-objective metaheuristic procedures used in the study.
NSGA-II: The NSGA-II algorithm, presented in Section 4, will be applied to the BGVRP in the cases study in question and the literature instances. For the case study, the parameters used, after several preliminary experiments, were: population = 100; vehicle capacity = 210 volumes of newspaper; number of iterations = 100; percentage of crossover = 0.5 and percentage of mutation = 0.1.
The formation of the route in the initial solution is random. To make the crossover, first of all, the VRP was transformed into a Traveling Salesman Problem (TSP), and then a random cut was made of one point between the second and second-last gene. The chromosome was then adapted through mutation to contain all the addresses on the route as a TSP. Finally, the TSP was “transformed” into a VRP making the “distribution” of the chromosome according to the “maximum travel time” along with “vehicle capacity” to validate the routes. Thus, the creation of a new route (chromosome) will be feasible. A flowchart of the implemented NSGA-II algorithm is shown in Figure 1.
Figure 1. Flowchart of the implemented NSGA-II algorithm.
MOPSO: The MOPSO algorithm, presented in Section 4, will be applied to the BGVRP in the case in question and the literature instances. The parameters and procedures for the case study will be the same as those presented in Section 5.2 (NSGA-II). Furthermore, for the calculation of speed and position, Equations (6) and (7) were selected, in accordance with the recommendations from the literature : acceleration constant or cognition rate C1 and C2 = 1; inertial weight W through Equation (8); Wmax = 1 and Wmin = 0.2.
It should be observed that the positions (routes) were simplified for the TSP and, following the calculation procedures of Equations (6) and (7), were readapted for the VRP. It was established that the velocity vector can contain up to N changes, where N is the number of addresses to be visited, in addition to the depot.
As the pBest particle, the evolution of the particle itself was considered. The gBest particle is the random choice of a particle from the repository. In the repository, the updated PF at each iteration is found. It should be highlighted that all the routes were validated by the capacity of the vehicle (motorcycle). A flowchart of the MOPSO algorithm is shown in Figure 2.
Figure 2. Flowchart of the implemented MOPSO algorithm.
CWNSGA-II Hybrid Algorithm: The proposed procedure of the Hybrid Algorithm of Clarke and Wright’s Savings and the Non-dominated Sorting Genetic Algorithm II (CWNSGA-II) is presented here. The C&W algorithm and the NSGA-II were described in Section 4. In the literature, correlated works can be found, such as that of Wang et al. , with the authors using the C&W as a procedure in the generation of the initial population. Here, the C&W was applied and a swap operation was used to achieve the established number of the population to form the initial population. The parameters and procedures are the same as those presented in 5.2 (NSGA-II) in the solution of the proposed case study. The C&W will be applied to refine the solutions generated after the crossover, more specifically, at each “petal” of the VRP solutions. A flowchart of the CWNSGA-II algorithm is presented in Figure 3.
CWTSNSGA-II Hybrid Algorithm: The proposed procedure of the Hybrid Algorithm of Clarke & Wright’s Savings, Tabu Search and the Non-dominated Sorting Genetic Algorithm II (CWTSNSGA-II) is presented. The C&W algorithm, TS and NSGA-II were described in Section 4. The parameters and procedures are the same as those presented in 5.2 (CWNSGA-II).
The differential of the CWTSNSGA-II is in the loop (Figure 4), as it follows
Figure 3. Flowchart of the CWNSGA-II algorithm.
Figure 4. Flowchart of the implemented CWTSNSGA-II algorithm in (a) Stages in the iterations of the implemented CWTSNSGA-II algorithm; (b) Initial population of the iteration; (c) Population after crossover (children); (d) Population after Tabu Search moves on the Pareto front; (e) Population after Tabu Search moves in the diversification of individuals with high CD; (f) New Population.
the same procedures as the NSGA-II, in addition to a number of TS moves, calculates fitness, removes duplicate solutions, updates the PF and CD and generates a new population. Furthermore, the petals of the VRP are improved by the C&W. The TS moves are swaps on the PF and in the highest value solution of CD, to increase the number of solutions close to the PF and identify isolated solutions. Thus, the algorithm considerably improves its population at each iteration. A flowchart of the CWTSNSGA-II algorithm and commented illustration of the stages in the iterations of this procedure is presented in Figure 4.
5.3. Stage 3: Analysis of the Results
To analyze the results, the PF resulting from each algorithm was initially considered, in other words, the solution with the greatest hypervolume of three simulations. Then, with the best solution from each algorithm, they were evaluated and compared according to their hypervolume and other evaluation metrics described below, and the most adequate solution was selected. Finally, the best route was selected according to the Utility Function.
The metrics used are related to quality, processing time, number of solutions, mean ideal distance  and hypervolume .
Quality Metric: Quality Metric or percentage of non-dominated solutions. Here, the PFs resulting from each algorithm are compared and evaluated in terms of dominance between the presented solutions. Finally, the percentage of non-dominated solutions for each algorithm is presented.
Processing time: The processing time of each algorithm in seconds (s);
Number of solutions: Number of solutions presented on the PF;
Mean Ideal Distance (MID): The MID uses an imaginary ideal point as a reference to calculate the mean Euclidean distances. This point is obtained using the minimum value of each objective function among all the solutions from the PF. Thus, the best algorithm will have the lowest MID value. The calculation procedure is presented in Equation (9).
where N is the number of elements in the set of non-dominated solutions of each algorithm; is the Euclidian distance between the ideal point and the non-dominated solution.
Hypervolume: The calculation of the hypervolume allows an estimation of the proximity of the solution in relation to the optimal PF. As a procedure, it calculates the volume of the covered region between the points of the solutions on the PF and a determined reference point. Thus, the sum of the hypercubes formed by the non-dominated solutions is calculated. This procedure is used when the optimal PF is not known. It helps to measure the quality of a solution compared with other solutions and, thus, the higher the hypervolume, the better the solution will be.
6. Results and Discussion
In this section, the results that were obtained using the techniques proposed in Section 5 are presented, analyzed and discussed. To test the results of the case study, as well as the literature instances, Matlab R2016b software was used, installed in a notebook with an Intel® CoreTM i7-7500U processor CPU@2.7GHz 2.90GHZ, and 16.00 GB of RAM, with a 64-bit operational system.
6.1. Case Study: Newspaper Distribution
Here, the results are presented for each of the techniques used (NSGA-II; MOPSO; CWNSGA-II and CWTSNSGA-II) for the newspaper distribution. In each technique, four rounds were used, varying the initial populations in each one, followed by the hypervolume to select the best solution. The initial population, final population and Pareto front of the best performance of each technique are presented in Figure 5.
As shown in Figure 5, the final population of the NSGA-II, CWNSGA-II and
Figure 5. Pareto front with the initial and final population (left) and Pareto front (right) using the NSGA-II (a), MOPSO (b), CWNSGA-II (c) and CWTSNSGA-II (d) techniques.
CWTSNSGA-II algorithms showed greater convergence than the final population of the MOPSO algorithm.
Figure 6 contains a graph with the Pareto fronts obtained from the techniques used for the solution to this case study. The feasible reference point selected by the authors is also presented, point = (5; 45), for the calculation of the hypervolume and solution for the Utility Function, point = (1.7; 5.0), the calculation of which will be presented later in the study.
As shown in Figure 6, the PF resulting from the CWNSGA-II and CWTSNSGA-II algorithms, located on the left, were closer and presented good results. The PF resulting from the NSGA-II and MOPSO algorithms were more distant and had worse results.
Figure 7 presents the result of the evaluation metrics. The main result to highlight is that the hypervolume of the CWNSGA-II algorithm was slightly better (greater) than that of the CWTSNSGA-II, although both algorithms presented the same number of solutions and same percentage of non-dominated solutions. The distancing between the solutions of the CWTSNSGA-II (Figure 6) resulted in a higher MID value, showing that it is unfeasible compared with the other solutions. The MID value for the CWNSGA-II (4.9) algorithm was slightly lower than that of the NSGA-II (5.0).
To complement these results, the result of the evaluation of the Pareto fronts was measured by hypervolume. It should be noted that the higher the value of the hypervolume, the better the set of solutions will be. Therefore, the calculation of the hypervolume reinforces the dominance of the solution of the CWNSGA-II (143.9) algorithm over that of the CWTSNSGA-II (139.2), NSGA-II (59.3) and MOPSO (17.0), respectively.
Figure 8 presents the calculations for the application of the Utility Function
Figure 6. Analysis of the Pareto fronts of the employed techniques.
Figure 7. Result of the evaluation metrics in the case study solutions.
Figure 8. Utility Function applied on the dominant PF (of CWNSGA-II). The first frame shows the value of the objective functions and the second frame, the utility function values.
on the dominant PF resulting from the CWNSGA-II algorithm. For this procedure, the decision maker was asked to assign weights to each of the objective functions. Thus, 70% importance was assigned to the first objective function and 30% to the second objective function. Finally, Solution 4, with coordinates (1.74; 5.00), presented the greatest Utility Function, and this was the proposed solution.
A comparison of the current solution and proposed solution for the printer is presented in Table 2. The optimized solution used the same number of routes (motorcycles) and was the best in terms of all the requirements. The highlight was the “balancing” of demand, i.e., a better balance with regard to the number of newspapers loaded into the bags (bags hanging from the sides of the motorcycle) for distribution.
Table 2 shows that CO2 emission (OF1) was optimized by 19.9% and difference in demand (OF2) by 87.5%. These are the main parameters measured in this study.
Figure 9 shows in graphic form a comparison of the current solution (left) and the optimized solution (right), and Table 3 shows a description of the routes of the sector of the optimized solution for the newspaper distribution.
As shown in Table 3, although more delivery points have been allocated to tours 2 and 4 than to tours 1 and 3, there was greater homogenization of demand between the routes.
Table 4 contains information obtained through the objective functions, more specifically, the routes of the sector. These were obtained from the Pareto front,
Table 2. Comparison of current and proposed solutions.
Table 3. Description of the routes of the proposed sector.
Table 4. Information obtained through the objective functions using the CWNSGA-II technique.
Figure 9. Comparison of the (a) current and (b) optimized solutions.
CO2 emission (OF1), distance traveled, fuel consumption, travel time, number of vehicles used and difference in demand (OF2).
It should also be highlighted in Table 4 that the 5 VRPs are the solutions (global) resulting from the PF of the CWNSGA-II technique for the sector in question. Each VRP addresses a set of four routes for the newspaper distribution.
Table 5 shows the stratification of the assigned capacity in each route of the sector (VRP tours) obtained by the Pareto front.
Table 5. Stratification of assigned capacity in each tour of the sector of solutions using the CWNSGA-II technique.
As shown in Table 5, the proposed optimized solutions (VRPs) have a much more balanced assigned capacity than the printer’s current solution.
6.2. Literature Instances
The algorithms presented in this work (NSGA-II; MOPSO; CWNSGA-II and CWTSNSGA-II) were applied to a set of literature Instances. For the case of the BGVRP, a set of 27 instances of the CVRP (more specifically, Augerat’s Set A), with the dimension varying from 32 to 80 cities (considering the depot), standard vehicle capacity of 100 units, and 5 to 10 routes in the solutions, available at https://neo.lcc.uma.es/vrp/vrp-instances/.
To evaluate the algorithms presented in this work through literature Instances, the same parameters of Hassanzadeh and Rasti-Barzoki  were used regarding the number of iterations, population and/or particle size, percentage of crossover and percentage of mutation; and Abad et al.  for the cognition rate. The assigned parameters are: Iterations 100; Population/Particle 60; Crossover 0.7; Mutation 0.3; C1 and C2 is 1.
The OFs that were used are the same as in the case study, i.e., “minimization of CO2 emissions” and “absolute minimization of the demands of each route in relation to the average demand”. To evaluate the instances, three simulations were performed and the response with the greatest hypervolume was selected in addition to the other evaluative metrics, all of which are presented in 5.3.
The results of the tests with the literature Instances are synthesized in Figure 10. The CWTSNSGA-II algorithm had the best performance in terms of rate of dominance, MID and Hypervolume, albeit with the worst computational time. The CWNSGA-II algorithm presented the highest number of solutions. In general, the CWNSGA-II and CWTSNSGA-II presented better results than the NSGA-II and MOPSO.
As also shown in Figure 10, the use of C&W to improve the routes led to the performance of the CWNSGA-II being almost double the computational time of the NSGA-II. Moreover, the systematized complexity in the CWTSNSGA-II algorithm resulted in it being slower than the others, although it had the best results in the tests.
Figure 10. Result of the literature instances.
This work presents a new approach to the BVRP that can be used in logistic distribution problems in general. The use of the methodology is illustrated here through a case study involving newspaper distribution. It should be highlighted that to collect the data for the case study, the asymmetry of distances traveled was considered and, consequently, the travel times, instead of the usual Euclidean (symmetric) distances.
Use was then made of the hybrid algorithms of the NSGA-II with C&W Savings (CWNSGA-II), as well as NSGA-II with C&W and TS (CWTSNSGA-II). Comparative tests of the algorithms were performed through hypervolume (the bigger, the better) and other evaluation metrics. Afterwards, the utility function was used to choose the best solutions, supported by the weights assigned by the decision maker.
To define the macro points of the case study, the “direction of the street” was considered rather than the mean point normally used by researchers. Because of this and the other viewpoints that were addressed, the deliveryman can carry out his tasks with fewer adaptations when making deliveries. The CWNSGA-II algorithms proved to be superior to the others. The optimization level that was achieved was 19.9% for OF1 (minimization of CO2 emissions), and consequently the same percentage of minimization for total distance, and 87.5% for OF2, (minimization of difference in demand).
In this context, multi-objective optimization was helpful in the conflict between minimizing CO2 emission and minimizing difference in demand. Thus, regarding OF1 it was possible to stratify distance travelled, fuel consumption and travel time. CO2 emission represents the environmental concerns that transformed this problem into a BGVRP. OF2, difference in demand, contributes to the homogenization of occupied vehicle capacity. The environmental concerns, currently with high values, are treated here as a mechanism for conscious decision making. This is why OF1 had a greater weight.
Finally, with regard to the literature Instances, CWNSGA-II and CWTSNSGA-II algorithms presented better results than the others.
There remain many opportunities with a view to exploring applications and techniques regarding works related to the BGVRP. In this context, the intention is to develop more scenarios, developing new multi-objective metaheuristic algorithms in addition to hybrid versions. The theoretical aspects presented here bring advances in the use of routing techniques and in the multi-objective approach, through a very original methodology. The use of two objectives simultaneously is more realistic, as well as result in lower operating costs and better organization of activities in the company will bring greater satisfaction to users.
This study was financed in part by PUCPR and the Coordination for the Improvement of Higher Education Personnel, Brazil (CAPES; 1st author) and by the National Council for Scientific and Technological Development, Brazil (CNPq; 2nd author).
 Erdelic, T. and Caric, T. (2019) A Survey on the Electric Vehicle Routing Problem: Variants and Solution Approaches. Journal of Advanced Transportation, 2019, Article ID: 5075671.
 Validi, S., Bhattacharya, A. and Byrne, P.J. (2015) A Solution Method for a Two-Layer Sustainable Supply Chain Distribution Model. Computers & Operations Research, 54, 204-217.
 Abad, H.K.E.A., Vahdani, B., Sharifi, M. and Etebari, F. (2018) A Bi-Objective Model for Pickup and Delivery Pollution-Routing Problem with Integration and Consolidation Shipments in Cross-Docking System. Journal of Cleaner Production, 193, 784-801.
 Ebrahimi, S.B. (2018) A Stochastic Multi-Objective Location-Allocation-Routing Problem for Tire Supply Chain Considering Sustainability Aspects and Quantity Discounts. Journal of Cleaner Production, 198, 704-720.
 Fathollahi-Fard, A.M., Hajiaghaei-Keshteli, M.H.K. and Tavakkoli-Moghaddam, R. (2018) A Bi-Objective Green Home Health Care Routing Problem. Journal of Cleaner Production, 200, 423-443.
 Amer, H., Salman, N., Hawes, M., Chaqfeh, M., Mihaylova, L. and Mayfield, M. (2016) An Improved Simulated Annealing Technique for Enhanced Mobility in Smart Cities. Sensors (Switzerland), 16, 1-23.
 Ghezavati, V.R. and Beigi, M. (2016) Solving a Bi-Objective Mathematical Model for Location-Routing Problem with Time Windows in Multi-Echelon Reverse Logistics Using Metaheuristic Procedure. Journal of Industrial Engineering International, 12, 469-483.
 Fu, P., Li, H., Wang, X., Luo, J., Zhan, S.L. and Zuo, C. (2017) Multiobjective Location Model Design Based on Government Subsidy in the Recycling of CDW. Mathematical Problems in Engineering, 2017, Article ID: 9081628.
 Sawik, B., Faulin, J. and Pérez-Bernabeu, E. (2017) A Multicriteria Analysis for the Green VRP: A Case Discussion for the Distribution Problem of a Spanish Retailer. Transportation Research Procedia, 22, 305-313.
 Toro, E.M., Franco, J.F., Echeverri, M.G. and Guimaraes, F.G. (2017) A Multi-Objective Model for the Green Capacitated Location-Routing Problem Considering Environmental Impact. Computers and Industrial Engineering, 110, 114-125.
 Soleimani, H., Chaharlang, Y. and Ghaderi, H. (2018) Collection and Distribution of Returned-Remanufactured Products in a Vehicle Routing Problem with Pickup and Delivery Considering Sustainable and Green Criteria. Journal of Cleaner Production, 172, 960-970.
 Braekers, K., Ramaekers, K. and Nieuwenhuyse, I.V. (2016) The Vehicle Routing Problem: State of the Art Classification and Review. Computers & Industrial Engineering, 99, 300-313.
 Kumar, R.S., Kondapaneni, K., Dixit, V., Goswami, A., Thakur, L.S. and Tiwari, M.K. (2016) Multi-Objective Modeling of Production and Pollution Routing Problem with Time Window: A Self-Learning Particle Swarm Optimization Approach. Computers and Industrial Engineering, 99, 29-40.
 Ehrgott, M., Wang, J.Y.T., Raith, A. and Van Houtte, C. (2012) A Bi-Objective Cyclist Route Choice Model. Transportation Research Part A: Policy and Practice, 46, 652-663.
 Steiner, M.T.A., Datta, D., Steiner Neto, P.J., Scarpin, C.T. and Figueira, J.R. (2015) Multi-Objective Optimization in Partitioning the Healthcare System of Parana State in Brazil. Omega, 52, 53-64.
 Lin, C., Choy, K.L., Ho, G.T.S., Chung, S.H. and Lam, H.Y. (2014) Survey of Green Vehicle Routing Problem: Past and Future Trends. Expert Systems with Applications, 41, 1118-1138.
 Norouzi, N., Sadegh-Amalnick, M. and Tavakkoli-Moghaddam, R. (2017) Modified Particle Swarm Optimization in a Time-Dependent Vehicle Routing Problem: Minimizing Fuel Consumption. Optimization Letters, 11, 121-134.
 Gong, X., Deng, Q., Gong, X., Zhang, L., Wang, H. and Xie, H. (2018) A Bee Evolutionary Algorithm for Multiobjective Vehicle Routing Problem with Simultaneous Pickup and Delivery. Mathematical Problems in Engineering, 2018, Article ID: 2571380.
 Psychas, I.D., Marinaki, M., Marinakis, Y. and Migdalas, A. (2017) Non-Dominated Sorting differential Evolution Algorithm for the Minimization of Route Based Fuel Consumption Multiobjective Vehicle Routing Problems. Energy Systems, 8, 785-814.
 Rabbani, M., Saravi, N.A. and Farrokhi-Asl, H. (2017) Design of a Forward/Reverse Logistics Network with Environmental Considerations. International Journal of Supply and Operations Management, 4, 115-132.
 Gupta A., Heng C.K., Ong Y.S., Tan P.S. and Zhang A.N. (2017) A Generic Framework for Multi-Criteria Decision Support in Eco-Friendly Urban Logistics Systems. Expert Systems with Applications, 71, 288-300.
 Hassanzadeh, A. and Rasti-Barzoki, M. (2017) Minimizing Total Resource Consumption and Total Tardiness Penalty in a Resource Allocation Supply Chain Scheduling and Vehicle Routing Problem. Applied Soft Computing, 58, 307-323.
 Liu, X.H., Shan, M.Y., Zhang, R.L. and Zhang, L.H. (2018) Green Vehicle Routing Optimization Based on Carbon Emission and Multiobjective Hybrid Quantum Immune Algorithm. Mathematical Problems in Engineering, 2018, Article ID: 8961505.
 Wang, Y., Peng, S., Assogba, K., Liu, Y., Wang, H., Xu, M. and Wang, Y. (2018) Implementation of Cooperation for Recycling Vehicle Routing Optimization in Two-Echelon Reverse Logistics Networks. Sustainability, 10, 1358.
 Wang, Z., Leng, L., Wang, S., Li, G. and Zhao, Y. (2020) A Hyperheuristic Approach for Location-Routing Problem of Cold Chain Logistics Considering Fuel Consumption. Computational Intelligence and Neuroscience, 2020, Article ID: 8395754.
 Subramanian, A., Penna, P.H.V., Ochi, L.S. and Souza, M.J.F. (2013) Um algoritmo heurístico baseado em iterated local search para problemas de roteamento de veículos. In: Lopes, H.S., Rodrigues, L.C.A. and Steiner, M.T.A., Eds., Meta-Heurísticas em Pesquisa Operacional, Ed. Omnipax, Curitiba, 165-180.
 Ferreira, J.C., Steiner, M.T.A. and Guersola, M.S. (2017) A Vehicle Routing Problem Solved through Some Metaheuristics Procedures: A Case Study. IEEE Latin America Transactions, 15, 943-949.
 Demir, E., Bektas, T. and Laporte, G. (2014) A Review of Recent Research on Green Road Freight Transportation. European Journal of Operational Research, 237, 775-793.
 Deb, K., Pratap, A., Aguarwal, S. and Meyarivan, T. (2002) A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation, 6, 182-197.
 Carvalho, C.H.R. (2011) 1606: Texto para discussao. Emissoes relativas de poluentes do transporte motorizado de passageiros nos grandes centros urbanos brasileiros. Instituto de Pesquisa Economica Aplicada (IPEA), Brasília.