Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 31 A Multi-Objective Solution of Green Vehicle Routing Problem Özgür KABADURMUŞ* 1 , Mehmet Serdar ERDOĞAN 1, Yiğitcan ÖZKAN 1 and Mertcan KÖSEOĞLU 1 1 Yasar University/International Logistics Management, Izmir, Turkey [Corresponding Author indicated by an asterisk *] Abstract— Distribution is one of the major sources of carbon emissions and this issue has been addressed by Green Vehicle Routing Problem (GVRP). This problem aims to fulfill the demand of a set of customers using a homogeneous fleet of Alternative Fuel Vehicles (AFV) originating from a single depot. The problem also includes a set of Alternative Fuel Stations (AFS) that can serve the AFVs. Since AFVs started to operate very recently, Alternative Fuel Stations servicing them are very few. Therefore, the driving span of the AFVs is very limited. This makes the routing decisions of AFVs more difficult. In this study, we formulated a multi-objective optimization model of Green Vehicle Routing Problem with two conflicting objective functions. While the first objective of our GVRP formulation aims to minimize total 𝑪𝑶 𝟐 emission, which is proportional to the distance, the second aims to minimize the maximum traveling time of all routes. To solve this multi-objective problem, we used 𝜺 -constraint method, a multi-objective optimization technique, and found the Pareto optimal solutions. The problem is formulated as a Mixed-Integer Linear Programming (MILP) model in IBM OPL CPLEX. To test our proposed method, we generated two hypothetical but realistic distribution cases in Izmir, Turkey. The first case study focuses on an inner-city distribution in Izmir, and the second case study involves a regional distribution in the Aegean Region of Turkey. We presented the Pareto optimal solutions and showed that there is a tradeoff between the maximum distribution time and carbon emissions. The results showed that routes become shorter, the number of generated routes (and therefore, vehicles) increases and vehicles visit a lower number of fuel stations as the maximum traveling time decreases. We also showed that as maximum traveling time decreases, the solution time significantly decreases. Index Terms— Green Vehicle Routing Problem, Alternative Fuel Vehicles, 𝜀 -Constraint, Multi-Objective Optimization, Pareto Optimality I. INTRODUCTION Factors such as technological developments, globalization and population growth have certain effects on the corporate strategy of companies. In the 1990s, the effects of globalization caused the competition to rise [1]. Therefore, firms were required to keep the price levels at a minimum. This situation forced many firms to lower internal costs in order to increase profits. The concept of Supply Chain Management was developed and studied for companies to minimize their operational costs. The main objective of Supply Chain Management is to provide the required products to the customers at the right time, at the right place, and at a low price. This goal can only be attained by effectively managing supply chain operations while minimizing operational costs. A basic supply chain would consist of suppliers, manufacturers, wholesalers (or distributors), and retailers. An illustration of a simple supply chain structure consisting of suppliers, manufacturers, distributors, retailers, and customers is given in Figure 1. Companies try to find the best way of improving their supply chain management by optimally planning their operations, such as strategic planning, demand planning, warehousing and logistics planning, manufacturing planning, inventory control, purchasing, and transportation planning. All these topics involve different decision processes. For example, in strategic planning, deciding the Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 32 location of factories and warehouses is the main focus, while in transportation planning finding the best the distribution plan among suppliers, manufacturers, retailers and end customers is the goal. Figure. 1. A typical supply chain Determination of a good transportation plan is difficult because the road infrastructure changes constantly due to the addition of new roads or the existence of roadworks. In addition, the variety of products has increased based on customer demands and it makes shipping activities even more difficult by forcing companies to find an optimal way of deliveries in order to minimize the distance, the fuel usage, 𝐶𝑂 2 emission, and therefore maximize their profits. Companies try to find more efficient ways to reach their customers from their depots and satisfy demand at the right time at the minimum cost. The Vehicle Routing Problem (VRP) addresses this issue. Vehicle Routing Problem is a generalization of Traveling Salesperson Problem (TSP) and it aims to reduce the transportation costs, avoid delivery delays, satisfy customer expectations, save fuel, and reduce environmental effects. Therefore, the VRP is very important for delivery costs and environmental effects. Figure 2 shows an example of VRP where the distributions to the customers are operated from a single depot using multiple routes (and vehicles). Figure. 2. An example of the VRP Green Vehicle Routing Problem (GVRP), first proposed by Erdogan and Miller-Hooks [2], focuses on the environmental aspects of the VRP. The use of vehicles creates pollution, wastage of fuel and traffic congestion. In the United States, 97% of energy for transportation involves the usage of gasoline. Throughout the world, new methods to reduce transportation pollution have been planned and implemented. For instance, municipal corporations, public enterprises, voluntary association, Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 33 and private companies are altering their traditional fleets of trucks to Alternative Fuel Vehicles (AFVs) that use non-petroleum fuel types such as biodiesel, electricity, and hydrogen. This is done mainly to minimize the negative impacts on the environment and adapt to the current environmental regulations. Green Vehicle Routing Problem includes a homogeneous fleet of AFVs and a set of Alternative Fuel Stations (AFSs) and aims to minimize total distance traveled by the vehicles. In this study, we developed a Multi-Objective Green Vehicle Problem (MOGVRP) by extending the GVRP. The problem has two conflicting objectives that are minimizing 𝐶𝑂 2 emission and the maximum traveling time of the routes. These two objectives are conflicting because as the maximum traveling time decreases, routes become shorter and the number of routes increases. Therefore, the total carbon emission, which is proportional to total distance, increases. We modeled the problem as a Mixed-Integer Linear Programming (MILP) model and used Epsilon Constraint Method to solve this multi-objective problem by using IBM OPL CPLEX. We also tested our model on two hypothetical but realistic case studies in Izmir, Turkey. The first case study focuses on an inner-city distribution in Izmir, and the second case study involves a regional distribution in the Aegean Region of Turkey. We obtained the distances and durations among locations (customers, alternative fuel stations, and the depot) by using Google Maps. We showed that there is a tradeoff between the maximum distribution time and carbon emissions and presented the Pareto optimal solutions. This study is organized as follows. In Section 2, the relevant VRP literature is summarized. In Section 3, the problem formulation is presented. Section 4 explains the solution methodologies used in our study. The results and managerial insights are provided in Section 5. The conclusion and future work are given in Section 6. II. LITERATURE REVIEW Green Vehicle Routing Problem was first proposed by [2] as an extension of the Vehicle Routing Problem and aims to minimize the total distance of deliveries to a set of customers using a homogenous fleet of alternative fuel vehicles that are originating from a single depot. It is important to point out that Alternative Fuel Stations are also added to the problem with the sole purpose of serving AFVs. Our study extends GVRP by adding another objective of minimizing the maximum travel time of the routes and therefore turning it to a multi-objective problem. In another GVRP study, [3] has designed a simulated annealing heuristic based exact solution to solve the GVRP. Pollution is an issue that has an impact on the environment; therefore, extensions of VRP relating to minimizing the negative effects of pollution were investigated in various studies. This line of research specifically names the problem as the Pollution-Routing Problem (PRP) and the major studies are [4], [5], [6], and [7]. The Pollution-Routing Problem is a recent extension of the Vehicle Routing Problem which decides the optimal routes for a set of vehicles to serve a number of customers while managing the speed of each vehicle on each route in order to minimize fuel consumption, emission and driver costs. Bektas and Laporte [4] introduced a more extensive objective function that accounts for the travel distance, the total amount of emissions, travel times, fuel consumption, and distribution costs. Demir et al. [5] presented an Adaptive Large Neighborhood Search (ALNS) algorithm, which integrates the classical ALNS scheme with a specially designed vehicle speed optimization algorithm for solving the PRP. A heterogeneous fleet of vehicles was added to the PRP by [6]. The main goal of [6] is to minimize the sum of vehicle fixed costs and routing costs where the latter includes the cost of fuel and 𝐶𝑂 2 emissions, and driver cost. Demir et al. [7] proposed a bi- objective PRP to minimize fuel consumption and total driving time. In their study, the combination of ALNS and a speed optimization procedure were used to solve the bi-objective PRP. There are a variety of extensions to the classical Vehicle Routing Problem in the literature, such as multi-depot, time-windows, and heterogeneous fleet. Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 34 The multi-depot vehicle routing problem (MDVRP) is the problem of delivering goods to a set of customers with a set of vehicles originating from several depots. It was initially presented by [8]. In another study, [9] applied the MDVRP in the healthcare sector using a heterogeneous fleet in which a set of heterogeneous vehicles (each belonging to a different facility) providing transportation services to patients. There are different vehicle types with different fleet sizes located in multiple healthcare facilities. Variable Neighborhood Search and Tabu Search algorithms were used to solve the problem in [9]. Another example is proposed by [10], which includes time windows of customers and a fleet of vehicles used for the delivery and installation of electronics. This is an important subject in the literature because transportation of installation services is a requirement for many companies that manufacture or sell products that are difficult to set up. Furthermore, [11] aimed to provide militaristic solutions for unmanned vehicle routing. A set of mixed ground and aerial vehicles at several depots traveling to certain targets were considered in their study and a branch-and-cut algorithm was developed. Crevier, Cordeau, and Laporte [12] presented a variant of the MDVRP by including intermediate depots to aid the vehicles to replenish. To solve the problem, a heuristic combining the Adaptive Memory Principle and a Tabu Search method was developed. Various heuristic solution methodologies presented by [13], [14], [15], and [16] to solve the MDVRP. The Vehicle Routing Problem with Time Windows (VRPTW) is a significant variation of VRP. In VRPTW, a vehicle must visit each customer within a specific time interval. The VRPTW was introduced to the literature by [17]. The model of this study aims to minimize route durations. [18] formulated the VRPTW as a set partitioning problem and solved the problem by using column generation. Furthermore, [19] formulated a Multi-Objective Genetic Algorithm for VRPTW. In their study, the first objective is to minimize the total distance while the second is to minimize the number of vehicles used. [20] proposed Vehicle Routing Problem with Soft Time Windows and Stochastic Travel Times where time windows of the customers are large. The objective function of the model considers both transportation cost and service cost and Tabu Search method was used to solve the model. Leung et al. [21] extended the VRPTW by considering simultaneous pick-up and delivery. In their study, customers require both pick-up and delivery services within the specified time windows. The main objective is to minimize the total distance traveled by the vehicles. VRP with the heterogeneous fleet (HVRP) includes multiple vehicle types with different capacities, fixed and variable costs. The main purpose of the HVRP is to minimize total variable routing cost and the vehicle fix costs. Angelelli and Mansini [22] proposed an HVRP with two-dimensional loading constraints. The problem considers two-dimensional loading configuration of products according to their sizes. A Simulated Annealing heuristic with a local search (SA_HLS) was used to solve the problem. Belfiore and Yoshizaki [23] extended the HVRP by considering time windows of customers and developed a Tabu Search algorithm to solve it. Furthermore, [24] developed an HVRP with time windows and split deliveries. They proposed a Scatter-Search (SS) approach to solving the model. Jair et al. [25] proposed the HVRP with time windows and multiple products where customers demand different types of products. In their study, an Ant Colony Optimization algorithm was developed to solve the problem. Furthermore, several multi-objective models were presented in the VRP literature ([7], [19], [26], [27]). For example, [26] proposed a hybrid meta-heuristic for multi-objective vehicle routing problems with time windows. Their model aims to minimize total distance traveled and the workload imbalance that is the ratio of distance traveled by a vehicle and load of the vehicle. The multi-objective optimization has also been used for militaristic purposes. Guerriero et al. [27] proposed a VRP with soft time windows for Autonomous Unmanned Aerial Vehicles (UAVs). The problem has three objectives that are the minimization of total distance traveled by the UAVs, minimization of the number of vehicles (UAVs) and maximization of customer satisfaction. Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 35 III. PROBLEM DEFINITION AND FORMULATION The aim of the GVRP is to find the optimal routes of a set of homogeneous AFV fleet that depart from a single depot in order to meet the customer demands. The fuel capacity of vehicles are set and Alternative Fuel Stations that vehicles can replenish their fuel are considered in the model. In addition, the maximum traveling time of a route is included in the model. In this study, we developed a Multi-Objective Green Vehicle Routing Problem (MOGVRP) that is an extension of the GVRP, which is originally proposed by [2]. The first objective of the problem is to minimize total 𝐶𝑂 2 emission that is proportional to total distance traveled by the vehicles. The second objective minimizes the maximum traveling time of all routes. The mathematical model is presented below. Sets: D Set of Depots C Set of Customers 𝐶 0 Set of Depot and Customers: 𝐷 ∪ 𝐶 F Set of Alternative Fuel Stations 𝐹 0 Set of Depot and Alternative Fuel Stations: 𝐷 ∪ 𝐹 V Set of All Nodes: 𝐷 ∪ 𝐶 ∪ 𝐹 E Set of Edges: (𝑣 𝑖 , 𝑣 𝑗 ): 𝑣 𝑖 , 𝑣 𝑗 𝜖 𝑉 Parameters: ∝ 𝐶𝑂 2 emission (in kg) per km 𝑑 𝑖𝑗 Distance from node i to node j: (𝑖 , 𝑗 ) 𝜖 𝐸 M A sufficiently large number r Vehicle fuel consumption rate (liters per km) Q Vehicle fuel tank capacity 𝑠 𝑖 Service time at node i: 𝑖 ∈ 𝑉 𝑡 𝑖𝑗 Travel time (hours) from node i to node j: (𝑖 , 𝑗 ) 𝜖 𝐸 Decision variables: 𝑥 𝑖𝑗 { 1 𝐼𝑓 𝑎 𝑣𝑒 ℎ𝑖𝑐𝑙𝑒 𝑡𝑟𝑎𝑣𝑒𝑙𝑠 𝑓𝑟𝑜𝑚 𝑖 𝑡𝑜 𝑗 0 𝑜𝑡 ℎ𝑒𝑟𝑤𝑖𝑠𝑒 } 𝑦 𝑗 The remaining fuel level before arriving to j. This variable is reset to Q upon visiting a fuel station or the depot 𝜏 𝑗 Arrival time of a vehicle at j (departure from the depot is assumed to be zero) 𝑇 𝑚𝑎𝑥 Maximum traveling time of all routes min ∑ ∝ 𝑑 𝑖𝑗 𝑥 𝑖𝑗 𝑖 ,𝑗 ∈𝑉 ′ 𝑖 ≠𝑗 (1) min 𝑇 𝑚𝑎𝑥 (2) 𝑠 . 𝑡 . Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 36 ∑ 𝑥 𝑖𝑗 = 1, ∀𝑖 ∈ 𝐼 𝑗 ∈𝑉 ′ 𝑗 ≠𝑖 (3) ∑ 𝑥 𝑖𝑗 𝑗 ∈𝑉 ′ 𝑗 ≠𝑖 ≤ 1, ∀𝑖 ∈ 𝐹 0 (4) ∑ 𝑥 𝑗𝑖 𝑖 ∈𝑉 ′ 𝑗 ≠𝑖 − ∑ 𝑥 𝑖𝑗 = 0, ∀𝑗 ∈ 𝑉 ′ 𝑖 ∈𝑉 ′ 𝑗 ≠𝑖 (5) ∑ 𝑥 0𝑗 𝑗 ∈𝑉 ′\{0} ≤ 𝑚 (6) ∑ 𝑥 𝑗 0 ≤ 𝑚 𝑗 ∈𝑉 ′\{0} (7) 𝜏 𝑗 ≥ 𝜏 𝑖 + (𝑡 𝑖𝑗 + 𝑝 𝑖 )𝑥 𝑖𝑗 − 𝑀 (1 − 𝑥 𝑖𝑗 ), ∀ 𝑖 ∈ 𝑉 ′ , 𝑗 ∈ 𝑉 ′ \{0} 𝑎𝑛𝑑 𝑖 ≠ 𝑗 (8) 0 ≤ 𝜏 0 ≤ 𝑇 𝑚𝑎𝑥 (9) 𝑡 0𝑗 ≤ 𝜏 𝑗 ≤ 𝑇 𝑚𝑎𝑥 − (𝑡 𝑗 0 + 𝑝 𝑗 ), ∀𝑗 ∈ 𝑉 ′ \{0} (10) 𝑦 𝑗 ≤ 𝑦 𝑖 − 𝑟 ∙ 𝑑 𝑖𝑗 𝑥 𝑖𝑗 + 𝑄 (1 − 𝑥 𝑖𝑗 ), ∀𝑗 ∈ 𝐼 𝑎𝑛𝑑 𝑖 ∈ 𝑉 ′ , 𝑖 ≠ 𝑗 (11) 𝑦 𝑗 = 𝑄 , ∀𝑗 ∈ 𝐹 0 (12) 𝑦 𝑗 ≥ 𝑚𝑖𝑛 {𝑟 ∙ 𝑑 𝑗 0 , 𝑟 ∙ (𝑑 𝑗𝑙 + 𝑑 𝑙 0 )}, ∀𝑗 ∈ 𝐼 , ∀𝑙 ∈ 𝐹 ′ (13) 𝑥 𝑖 ,𝑗 ∈ {0,1}, ∀𝑖 , 𝑗 ∈ 𝑉 (14) Objective Function (1) minimizes the total 𝐶𝑂 2 emission while objective function (2) minimizes the maximum traveling time of all routes. Constraints (3) ensures that a vehicle arrives at a customer exactly once, while Constraint (4) ensures that a vehicle departs from a customer exactly once. Constraints (5) is the flow balance constraint. Constraint (6) and (7) limit the number of available vehicles that may depart from and arrive at the depot. Constraint (8) calculates the time of arrival at each location. Constraint (9) sets the minimum departure time of the vehicles from the depot to zero and limits the arrival times of vehicles to the depot to 𝑇 𝑚𝑎𝑥 . Constraint (10) guarantees that the duration of each route does not exceed the maximum traveling time. Constraint (11) calculates the vehicle fuel level when it arrives at a customer. Constraint (12) guarantees that when a vehicle visits an AFS, the fuel tank of the vehicle becomes full. Constraint (13) ensures that vehicles have enough fuel to return the depot when they arrive at the last customer in the route. Lastly, Constraint (14) is the boundary constraint. IV. METHODOLOGY A. Multi-Objective Optimization Without loss of generality, a multi-objective optimization problem can be defined as a minimization problem of the form as shown in (14) through (16). Minimize 𝑓 (𝑥 ) = (𝑓 1 (𝑥 ), 𝑓 2 (𝑥 ) … … 𝑓 𝑛 (𝑥 )) (14) Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 37 Subject to; 𝑔 𝑖 (𝑥 ) ≤ 0 ∀𝑖 = 1, … . . , 𝑝 (15) ℎ 𝑖 (𝑥 ) = 0 ∀𝑖 = 1, … . . . , 𝑞 (16) Many real-world optimization problems involve multiple objectives. Since the objectives are generally conflicting, there is no way to find an optimal solution which is best for all objectives. Instead, there is a set of optimal solutions (i.e., Pareto optimal solutions) can be found. As in (14), a reduction in any objective function value can only be achieved by an increase in another objective function value. In the Pareto optimal set of solutions, no solution dominates another. By definition, solution A dominates solution B if at least one objective of “A” is better than “B” and the remaining objectives are at least as good as the ones of “B”. In this example, “A” is called a non-dominated solution if there is no other solution dominates “A”. Note that, the set of Pareto Optimal solutions contains only non-dominated solutions. The primary goals in the multi-objective optimization are: 1) Preserving non-dominated solutions, 2) Improving progress of the algorithm toward the Pareto front, 3) Maintaining solution diversification in the Pareto front, 4) Providing an abundant number of non-dominated solutions to the decision maker. Multi-objective optimization formulations are applied in various engineering optimization problems. Minimizing cost, emission, time, work-in-process and stocks and maximizing profit, performance and customer satisfaction are the examples of multiple conflicting objectives of many real-life engineering problems. In multi-objective optimization, there is no single optimal solution but a set of Pareto optimal solutions that are equally good. At this point, it is the decision maker’s choice to select one of the solutions from the Pareto Optimal Solutions set. In our model, we have two main objectives. Those are minimizing the maximum traveling time of all routes and minimizing the total 𝐶𝑂 2 emission. To obtain the true Pareto optimal solutions of the problem, we used 𝜀 -constraint method. B. 𝜀 -Constraint Method We applied 𝜀 -constraint method to analyze the trade-off between 𝐶𝑂 2 emission and maximum traveling time. 𝜀 -constraint method is one of the most common optimization methods in the multi- objective literature. Equations (17) and (18) explains the procedure of the 𝜀 -constraint method. For a multi-objective optimization problem, as shown in (17), x is the decision variable vector and 𝑓 1 (𝑥 ) through 𝑓 𝑛 (𝑥 ) are the n objectives of the problem that are to be minimized. Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 38 min(𝑓 1 (𝑥 ), 𝑓 2 (𝑥 ) … … 𝑓 𝑛 (𝑥 ) ) (17) 𝜀 -constraint method optimizes one of the objective functions by converting the other objectives into constraints as shown in (18). min 𝑓 1 (𝑥 ) Subject to; 𝑓 2 (𝑥 ) ≤ 𝑐 2 𝑓 2 (𝑥 ) ≤ 𝑐 2 (18) ...... ...... 𝑓 𝑛 (𝑥 ) ≤ 𝑐 𝑛 The method solves the problem iteratively by changing the right-hand-side (RHS) of each constrained objective, 𝑐 𝑖 , by a small 𝜀 value at each iteration. The iterations start from the upper bound of each objective and end at its lower bound. Thus, a different solution is obtained at each iteration. Various multi-objective approaches were developed in the literature such as NSGA-II, SPEA-II, MOEA/D and weighted sum method. We have chosen 𝜀 -constraint method because it has some advantages over the other methods. These advantages of the 𝜀 -constraint method are summarized as follows: 1- It gives the true Pareto Optimal solutions or yields a very close approximation, 2- Each iteration generates a different solution, 3- The number of obtained solutions can be pre-adjusted by changing the value of 𝜀 . In our problem, we have two objectives: (1) minimization of total carbon emissions, and (2) minimization of maximum traveling time of all routes. Equation (19) explains the procedure of the 𝜀 - constraint for our problem. The RHS value changes at each iteration t by an 𝜀 𝑡 value. Therefore, each iteration generates different 𝐶𝑂 2 emission and maximum traveling time values. min ∑ ∝ 𝑑 𝑖𝑗 𝑥 𝑖𝑗 𝑖 ,𝑗 ∈𝑉 ′ 𝑖 ≠𝑗 s.t. (19) 𝑇 𝑚𝑎𝑥 ≤ 𝑅𝐻𝑆 𝑇 𝑚𝑎𝑥 − 𝜀 𝑡 𝜀 𝑡 = 𝑡 ∗ 𝜀 V. RESULTS We created two different distribution scenarios: (1) City distribution (İzmir, Turkey), and (2) Regional distribution (Aegean region of Turkey). For both scenarios, we selected the locations of 13 different customers, 3 fuel stations and 1 depot using Google Maps. Tables 1 and 2 show the names and coordinates of the selected locations. After obtaining the coordinates of the locations, we generated the distance and traveling time matrices by using Google Maps Api in R programming language. Service times (loading and unloading times) at customers are randomly generated between 15 and 45 minutes, while the service times of alternative fuel stations are randomly generated between 5 and 10 minutes. The fuel capacity of vehicles is selected as 100 liters for İzmir city case, while it is selected as 500 liters for Aegean region case because we assumed that small- Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 39 sized vehicles are used in inner-city distribution. Lastly, the number of available vehicles is selected as 8. Table 1. The names and coordinates of the locations of Izmir city scenario Location No Location Name Latitude Longitude 0 Kemalpasa (Depot) 38.443019 27.382101 1 Urla 38.325408 26.766736 2 Torbali 38.15387 27.3613 3 Buca 38.384972 27.173453 4 Gaziemir 38.325188 27.12748 5 Alsancak 38.440644 27.154587 6 Menderes 38.251641 27.134917 7 Narlibahce 38.392787 27.008523 8 Menemen 38.609575 27.068201 9 Foca 38.662059 26.753541 10 Bornova 38.486544 27.212053 11 Camonu 38.102358 27.15145 12 Bademler 38.273696 26.829623 13 VillaKent 38.602069 26.919441 14 Torbali (Fuel Station) 38.209473 27.33527 15 Balcova (Fuel Station) 38.397213 27.066138 16 Menemen (Fuel Station) 38.682503 27.010836 Table 2. The names and coordinates of the location of Aegean region scenario Location No Location Name Latitude Longitude 0 Kemalpasa(Depot) 38.443019 27.382101 1 Akhisar 38.934484 27.845407 2 Ayvalik 39.263005 26.732903 Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 40 3 Demirci 39.047515 28.654982 4 Odemis 38.230259 27.971609 5 Aydin 37.849461 27.823558 6 Mugla 37.211026 28.362027 7 Denizli 37.774926 29.078012 8 Usak 38.675099 29.405645 9 Kutahya 39.419277 29.981497 10 Dursunbey 39.586409 28.628718 11 Afyon 38.768686 30.490718 12 Burdur 37.717731 30.281475 13 Susurluk 39.861124 28.153807 14 Civril (Fuel Station) 38.300789 29.735007 15 Sogutcuk (Fuel Station) 37.496792 28.109612 16 Kirkagac (Fuel Station) 39.273171 27.897975 The MILP model was coded in IBM CPLEX Optimization Studio and executed in an 8GB RAM, Intel® Core™ i7-4510U CPU @2.00 GHz processor Linux computer. 𝜀 -Constraint method is used to obtain Pareto Optimal solutions. Figure 6 shows the Pareto Optimal solutions of the İzmir city scenario. 28 non-dominated solutions were obtained. The solution with the minimum traveling time is at the one extreme and has 827.234 𝐶𝑂 2 kg emission and 3.7 hours maximum traveling time. The solution with the minimum total emission is at the other extreme and has 420.084 kg 𝐶𝑂 2 emission and 8.05 hours of maximum traveling time. Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 41 Figure. 6. Pareto optimal solutions for İzmir City scenario Figure 7 illustrates the optimal routes of these two extreme solutions of Izmir city scenario. On the map, Depot, Customers, and Alternative Fuel Stations are numbered. The names and coordinates of these locations can be found in Table 1. As shown in Figure 7, the number of routes increases and routes become shorter as the maximum traveling time decreases for İzmir city scenario. Thus, no alternative fuel station was visited for the scenario with the lowest 𝑇 𝑚𝑎𝑥 value (Figure 7a), however, five vehicles were used. In the scenario with the highest 𝑇 𝑚𝑎𝑥 value (Figure 7b), only two routes were used and the vehicle refuels at an alternative fuel station in the longest route. a) Lowest 𝑇 𝑚𝑎𝑥 and highest emission b) Highest 𝑇 𝑚𝑎𝑥 and lowest emission Figure 7. The optimal routes of two extreme solutions of Izmir city scenario Figure 8 shows the Pareto Optimal solutions of the Aegean region scenario. 15 non-dominated solutions were obtained in this scenario. The solution with the lowest 𝑇 𝑚𝑎𝑥 value generated the highest 𝐶𝑂 2 emission (3801.779 kg) and 10.7 hours of maximum traveling time. The solution with the highest 𝑇 𝑚𝑎𝑥 value yielded the lowest 𝐶𝑂 2 emission (2329.912 kg) and 16 hours of maximum traveling time. Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 42 Figure 8. Pareto optimal solutions for Aegean region scenario Figure 9 illustrates the optimal routes of these two extreme solutions of the Aegean region scenario. On the maps, Depot, Customers, and Alternative Fuel Stations are numbered. The names and coordinates of these locations can be found in Table 2. According to Figure 9, more routes (five) were used, routes become shorter, and no alternative fuel station was visited when maximum traveling time is the lowest (Figure 9a) in Aegean distribution scenario. Only three routes were used and one of the vehicles refueled at an alternative fuel station when the maximum traveling time is the highest (Figure 9b). a) Lowest 𝑇 𝑚𝑎𝑥 and highest emission b) Highest 𝑇 𝑚𝑎𝑥 and lowest emission Figure 9. Optimal routes of two extreme solutions of Aegean region scenario Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 43 As for the solution times, the mean and standard deviation of solution times are 36.668 and 64.819, while the minimum and maximum solution times are 1.129 and 340.833 for Izmir city scenario. The mean and standard deviation of solution times are 11.327 and 12.808, while the minimum and maximum solution times are 0.585 and 43.663 for Aegean region scenario. As 𝑇 𝑚𝑎𝑥 reduces, solution time also reduces for both scenarios. VI. CONCLUSION The Vehicle Routing Problem is one of the most popular combinatorial optimization problems. The aim of the VRP is to deliver goods to a set of customers with a fleet of vehicles departing from a depot at a minimum cost. Rather than minimizing cost, the aim of the VRP may also be minimizing emission, time, and distance or maximizing the service level. As an extension to the original VRP problem, the Green Vehicle Routing Problem focuses on the environmental aspects of the VRP. It considers a homogeneous fleet of Alternative Fuel Vehicles originating from a single depot and uses a set of Alternative Fuel Stations to refuel the vehicles. It aims to minimize the total distance of the routes. In this paper, we formulated a Multi-Objective Green Vehicle Routing Problem by extending the GVRP. Our multi-objective problem includes two objectives: (1) minimizing the maximum traveling time of all routes, and (2) minimizing the total 𝐶𝑂 2 emission, which is proportional to the total distance. These two objectives are conflicting because as the maximum traveling time decreases, the routes become shorter and more routes are used. Therefore, the total 𝐶𝑂 2 emission increases. We modeled the problem as a Mixed Integer Linear Programming model and coded in IBM OPL CPLEX. Since the problem is multi-objective, we applied 𝜀 -Constraint method to generate Pareto optimal solutions. The proposed method has been tested on two hypothetical but realistic case studies. The first case study considers a city distribution, whereas the second one addresses a regional distribution. The data of the problem are generated to reflect a real-life distribution operation. We determined 13 customers, three AFSs and one depot locations in both scenarios. Results show that 𝐶𝑂 2 emission decreases as maximum traveling time increases. In addition, routes become shorter, the number of generated routes (and therefore, vehicles) increases and vehicles visit a lower number of fuel stations as the maximum traveling time decreases. Also, as maximum traveling time decreases, the solution time significantly decreases. Furthermore, CPLEX solves the problem very fast for 17 locations. However, when the number of locations increases above 20, the solution time starts to increase dramatically. Thus, the future study may consider developing an effective heuristic algorithm to solve the large instances in a shorter time. In addition, extending the problem by considering different AFV types, the impact of vehicle load on emission, or multi-depot can be other directions for future research. REFERENCES 1. Dunning, J. H. (2014). The Globalization of Business (Routledge Revivals): The Challenge of the 1990s. London: Routledge. 2. Erdoğan, S. & Miller-Hooks, E. (2012). A green vehicle routing problem. Transportation Research Part E: Logistics and Transportation Review, 48(1), 100-114. 3. Koç, Ç. & Karaoglan, I. (2016). The green vehicle routing problem: A heuristic based exact solution approach. Applied Soft Computing, 39, 154-164. 4. Bektaş, T. & Laporte, G. (2011). The pollution-routing problem. Transportation Research Part B: Methodological, 45 (8), 1232-1250. 5. Demir, E., Bektaş, T. & Laporte, G. (2012). An adaptive large neighborhood search heuristic for the pollution-routing problem. European Journal of Operational Research, 223(2), 346-359. 6. Koç, Ç., Bektaş, T., Jabali, O. & Laporte, G. (2014). The fleet size and mix pollution-routing problem. Transportation Research Part B: Methodological, 70, 239-254. 7. Demir, E., Bektaş, T. & Laporte, G. (2014). The bi-objective pollution-routing problem. European Journal of Operational Research, 232(3), 464-478. 8. Laporte, G., Nobert, Y. & Taillefer, S. (1988). Solving a family of multi-depot vehicle routing and location-routing problems. Transportation Science, 22(3), 161-172. 9. Detti, P., Papalini, F. & de Lara, G. Z. M. (2017). A multi-depot dial-a-ride problem with heterogeneous vehicles and compatibility constraints in healthcare. Omega, 70, 1-14. Logistics & Sustainable Transport Vol. 10, No. 1, June 2019, 31-44 doi: 10.2478/jlst-2019-0003 44 10. Bae, H. & Moon, I. (2016). Multi-depot vehicle routing problem with time windows considering delivery and installation vehicles. Applied Mathematical Modelling, 40(13-14), 6536-6549. 11. Sundar, K. & Rathinam, S. (2017). Algorithms for heterogeneous, multiple depot, multiple unmanned vehicle path planning problems. Journal of Intelligent & Robotic Systems, 88(2-4), 513-526. 12. Crevier, B., Cordeau, J. F. & Laporte, G. (2007). The multi-depot vehicle routing problem with inter-depot routes. European Journal of Operational Research, 176(2), 756-773. 13. Cordeau, J.F. & Maischberger, M. (2012). A parallel iterated tabu search heuristic for vehicle routing problems. Computers & Operations Research, 39(9), 2033-2050. 14. Escobar, J. W., Linfati, R., Toth, P. & Baldoquin, M. G. (2014). A hybrid granular tabu search algorithm for the multi-depot vehicle routing problem. Journal of Heuristics, 20(5), 483-509. 15. Shimizu, Y. & Sakaguchi, T. (2014). A hierarchical hybrid meta-heuristic approach to coping with large practical multi- depot VRP. Industrial Engineering & Management Systems, 13(2), 163-171. 16. Subramanian, A., Uchoa, E. & Ochi, L. S. (2013). A hybrid algorithm for a class of vehicle routing problems. Computers & Operations Research, 40(10), 2519-2531. 17. Savelsbergh, M. W. P. (1992). The vehicle routing problem with time windows: Minimizing route duration. ORSA Journal on Computing, 4(2), 146-154. 18. Desrochers, M., Desrosiers, J. D. & Solomon, M. (1992). A new optimization algorithm for the vehicle routing problem with time windows. Operations Research, 40(2), 342-354. 19. Ombuki, B., Ross, B. J. & Hanshar, F. (2006). Multi-objective genetic algorithms for vehicle routing problem with time windows. Applied Intelligence, 24(1), 17-30. 20. Taş, D., Dellaert, N., Van Woensel, T. & De Kok, T. (2013). Vehicle routing problem with stochastic travel times including soft time windows and service costs. Computers & Operations Research, 40(1), 214-224. 21. Angelelli, E. & Mansini, R. (2002). The vehicle routing problem with time windows and simultaneous pick-up and delivery. In Klose, A., Speranza, M. G. & Van Wassenhove, L. N. (Eds.), Quantitative approaches to distribution logistics and supply chain management (pp. 249-267). Berlin: Springer. 22. Leung, S. C. H., Zhang, Z., Zhang, D., Hua, X. & Lim, M. K. (2013). A meta-heuristic algorithm for heterogeneous fleet vehicle routing problems with two-dimensional loading constraints. European Journal of Operational Research, 225(2), 199-210. 23. Jiang, J., Ng, K. M., Poh, K. L. & Teo, K. M. (2014). Vehicle routing problem with a heterogeneous fleet and time windows. Expert Systems with Applications, 41(8), 3748-3760. 24. Belfiore, P. & Yoshizaki, H. T. Y. (2013). Heuristic methods for the fleet size and mix vehicle routing problem with time windows and split deliveries. Computers & Industrial Engineering, 64(2), 589-601. 25. Jair, J., Paternina-Arboleda, C. D., Cantillo, V. & Montoya-Torres, J. R. (2013). A two-pheromone trail ant colony system— Tabu search approach for the heterogeneous vehicle routing problem with time windows and multiple products. Journal of Heuristics, 19(2), 233-252. 26. Baños, R. L., Ortega, J., Gil, C. N., Márquez, A. L. & De Toro, F. (2013). A hybrid meta-heuristic for multi-objective vehicle routing problems with time windows. Computers & Industrial Engineering, 65(2), 286-296. 27. Guerriero, F., Surace, R., Loscri, V. & Natalizio, E. (2014). A multi-objective approach for unmanned aerial vehicle routing problem with soft time windows constraints. Applied Mathematical Modelling, 38(3), 839-852. AUTHORS A. Özgür KABADURMUŞ, is an Assistant Professor at the Department of International Logistics Management, Yasar University (Izmir, Turkey). (e-mail: ozgur.kabadurmus@yasar.edu.tr) B. Mehmet Serdar ERDOĞAN, is a Research Assistant at the Department of International Logistics Management, Yasar University (Izmir, Turkey). (e-mail: mehmet.erdogan@yasar.edu.tr) C. Yiğitcan ÖZKAN, BA, is a student at the Department of International Logistics Management, Yaşar University, Izmir, Turkey (e-mail: yigitcan.ozkan@stu.yasar.edu.tr) D. Mertcan KÖSEOĞLU, MA, is a student at the Department of International Logisitics Management, Yaşar University, Izmir, Turkey (e-mail: mertcan.koseoglu@stu.yasar.edu.tr) Manuscript received by 23 February 2019. Published as submitted by the author(s).