A New Approach to the Bi-objective Green Vehicle Routing Problem: Optimization in Newspaper Distribution

The purpose of this work is to present a methodology to provide a solution to a Bi-objective Green Vehicle Routing Problem (BGVRP). The methodology, illustrated using a case study (newspaper distribution problem) and instances from the literature, was divided into three stages: Stage 1, data treatment; Stage 2, metaheuristic approaches (hybrid or non-hybrid), used comparatively, and, Stage 3, analysis of the results, with a comparison of the algorithms

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 (Braekers, Ramaekers, & Nieuwenhuyse, 2016).
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 (Erdelic & Caric, 2019).
Meanwhile, the field of Multi-Objective Optimization (MOO) addresses the process of optimizing two or more conflicting objectives simultaneously, subject to constraints (Kumar et al., 2016). 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 (Ehrgott, Wang, Raith, &Van Houtte, 2012 andSteiner, Datta, Steiner Neto, Scarpin, &Figueira, 2015).
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 (Lin, Choy, Ho, Chung, & Lam, 2014;Norouzi, Sadegh-Amalnick, & Tavakkoli-Moghaddam, 2017;Gong et al., 2018).
The aim of this work is to present a methodology to solve a BGVRP, illustrated using a newspaper distribution problem and instances from the literature. In this context, a multiobjective approach could aid the complex decision-making process in order to achieve the intended objectives. which begins at 21h on Fridays, with the deliveries finalized around 8h 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 8,500 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 2h, 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 N or NM (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 N and NM. 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.

Literature Review
In this section, some concepts related to the Green VRP are addressed, in addition to some correlated works.

Green VRP
To address the theme of the Green VRP, the classification proposed by Lin et al. (2014) 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.

Correlated Works
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. Govindan et al. (2014) presented the two-echelon location-routing problem with time-windows (2E-LRPTW), applied in a case study to address perishable food distribution. The authors used multi-objective particle swarm optimization (MOPSO), adapted multi-objective variable neighborhood search (AMOVNS) and the multiobjective hybrid approach.
Through multi-objective route based fuel consumption vehicle routing problems (MRFCVRPs), Psychas, Marinaki, Marinakis, and Migdalas (2017) used the Pareto Front (PF), applied to parallel multi-start non-dominated sorting differential evolution algorithms (PMS-NSDEs) 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 (2017). The authors compared two metaheuristic procedures, NSGA-ӀӀ and MOPSO, with the former presenting the best performance.

Theoretical Foundation
This section contains a brief presentation of the main concepts and algorithms that provide the foundation for the present study.

4.1.Vehicle Routing Problem
In the literature, some mathematical approaches to solving the VRP can be found, such as those proposed by Fisher and Jaikumar (1981) and Subramanian, Penna, Ochi and Souza (2013). According to Lin et al. (2014), the VRP is a problem with very wide-ranging applications in the fields of distribution and logistics management, such as newspaper distribution (Ferreira, Steiner, & Guersola, 2017), and numerous other possible applications.

4.2.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 .

4.3.Tabu Search
According to Ferreira et al. (2017), Glover (1986) 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
Multi-Objective Optimization (MOO) is a process that optimizes two or more conflicting objectives simultaneously, subject to constraints (Demir, Bektas, & Laporte, 2014). When addressing the problem of minimization, the mathematical model is defined as in (1) (2) and (3) define the space of the decision variables in R n 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 (Wang et al., 2020).

4.5.NSGA II Algorithm
The Genetic Algorithm (GA), developed by Holland (1975) 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 (2002), 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-nondominated-sorting (FNS) procedure, in which the classification of the population at different levels through the use of Pareto dominance occurs.

4.6.MOPSO Algorithm
According to Abad et al. (2018), Multi-objective particle swarm optimization (MOPSO), proposed by Coello (2000), 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 (2001) can be used, as shown in (8), where wmax and wmin address, respectively, the maximum and minimum amount of inertia.

4.7.Utility Theory
According to Poterba and Rotemberg (2018), 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.

4.8.Environmental Considerations
Here, the environmental considerations related to GHG emissions are addressed. According to Carvalho (2011), 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.

Proposed Methodology
The proposed methodology to resolve the BGVRP for the study consisted of three stages: Evaluation of the case study and instances from the literature.

5.1.Stage 1: Data "Treatment"
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 contain 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.5, will be applied to the BGVRP in the cases study in question and the instances from the literature. 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.

MOPSO
The MOPSO algorithm, presented in Section 4.6, will be applied to the BGVRP in the case in question and the instances from the literature. The parameters and procedures for the case study will be the same as those presented in Section 5.2.1 (NSGA-II). Furthermore, for the calculation of speed and position, Equations (6) and (7)  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.
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.

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 was described in Section 4.2, and the NSGA-II in Section 4.5. In the literature, correlated works can be found, such as that of Wang et al. (2018), 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.1 (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 was described in Section 4, the TS in Section 4.3 and the NSGA-II in Section 4.5. The parameters and procedures are the same as those presented in 5.2.3 (CWNSGA-II).
The differential of the CWTSNSGA-II is in the loop (Figure 4), as it follows 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. Source: The authors.

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.
where N is the number of elements in the set of non-dominated solutions of each algorithm; 1 2 + 2 2 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.

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 instances from the literature, Matlab R2016b software was used, installed in a notebook with an Intel® Core™ 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 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.    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 form the sides of the motorcycle) for distribution.       Source: The authors.
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.  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.  197 200 195 199 VRP 4 199 197 199 196 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.Instances fom the Literature
The algorithms presented in this work (  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.

Conclusions
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). In the CWNSGA-II algorithm (Figure 3), the initial solution was obtained through the application of the C&W with a TSP, making swap moves until a number of individuals for the initial population was generated and then adapting the TSP to the VRP. C&W was also used to refine the solutions generated after the crossover, more specifically in each "petal" of the VRP solutions. The differential of the CWTSNSGA-II is in the loop shown in Figure 4, which follows the same procedures as the CWNSGA-II and a certain number of TS moves: it calculates fitness, excludes duplicates, updates the PF, CD and generates a new population. Moreover, the C&W resulted in an improvement in the "petals" of the VRP. The TS moves are swaps on the PF and in the solution with the highest CD value to expand the number of solutions close to the PF and diversify isolated solutions. Therefore, the algorithm improves its population considerably at each iteration.
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 instances from the literature, 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, always seeking to improve results and achieve scientific advances through new techniques, using the expanded version of the database, with approximately 8,500 delivery points (newspaper subscribers).
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.