#### Abstract

This research studies the capacitated green vehicle routing problem (CGVRP), which is an extension of the green vehicle routing problem (GVRP), characterized by the purpose of harmonizing environmental and economic costs by implementing effective routes to meet any environmental concerns while fulfilling customer demand. We formulate the mathematical model of the CGVRP and propose a simulated annealing (SA) heuristic for its solution in which the CGVRP is set up as a mixed integer linear program (MILP). The objective of the CGVRP is to minimize the total distance traveled by an alternative fuel vehicle (AFV). This research conducts a numerical experiment and sensitivity analysis. The results of the numerical experiment show that the SA algorithm is capable of obtaining good CGVRP solutions within a reasonable amount of time, and the sensitivity analysis demonstrates that the total distance is dependent on the number of customers and the vehicle driving range.

#### 1. Introduction

The increase in environmental consciousness has turned the green supply chain into a critical issue. Since 1990, many companies have been managing their supply chains in line with environmental regulations imposed by government sectors [1, 2]. Zhu and Sarkis [3] researched green supply chain management practices in Chinese manufacturing enterprises and stressed that globalization may result in more pressure to improve their environmental performance. Environmental issues affect supply chain decisions, such as sourcing of raw materials, location, transportation planning, and modal selection [2], and supply chain decisions determine both the network cost and environmental impact [4]. Therefore, a study considering the environmental impact on the supply chain is needed.

As a result of the negative impact on the environment, the most crucial green supply chain decision at the operation level is transportation planning [5]. The transportation sector plays a crucial role in producing emissions, particularly greenhouse gases (GHG). In 2009, the transportation sector contributed approximately 24% of all CO2 emissions from fossil fuel combustion, with road transportation comprising 16.7% of the overall GHG emissions [6]. Because transportation is the most important economic activity in the supply chain, transportation cost is one of the main components of the total logistics cost. Stock and Lambert [7] report that 60% of total logistics cost is in transportation. Bowersox et al. [8] noted that the most substantial element of logistical cost in a supply chain is due to the mode of transportation. Thus, effective and efficient transportation planning when considering the choice of vehicle type in vehicle routing can help mitigate these problems.

Some vehicle type characterizations include the body type, make/model, fuel type, vintage, and vehicle acquisition type [9]. The type of fuel a vehicle uses has a significant impact on the environment. There are two key driving factors when developing an alternative fuel vehicle (AFV), namely, energy security and GHG emissions reduction [10]. Therefore, AFVs play a crucial role in green transportation [11].

Arranging effective and efficient vehicle routing, usually referred to as the vehicle routing problem (VRP), is a way of mitigating the negative impact of transportation on the environment. Specifically, the green VRP (GVRP) is a variant of the VRP that focuses on solving the transportation problem in the green supply chain. The GVRP encompasses a broad and extensive class of problems concerning how to mitigate environmental issues as well as the design of a vehicle’s route. Research on the green VRP spans from the reduction of fuel consumption to designing an efficient vehicle route. However, research on designing the VRP to mitigate environmental issues using AFVs is still limited. Therefore, this present paper focuses on designing vehicle routes in the GVRP by considering AFVs as a way to mitigate environmental issues. Furthermore, by considering the vehicle capacity, this research extends that of Erdoğan and Miller-Hooks [12]. Hence, the concept of a capacitated GVRP (CGVRP) is introduced.

The CGVRP, a variant of the GVRP that is an NP-Hard, belongs to the same class. Thus, a simulated annealing (SA) heuristic approach is proposed. SA is one of the most common metaheuristic approaches and has been successfully used to solve the VRP and its variants [13–18]. Metaheuristics can consistently produce high-quality solutions, despite their disadvantage regarding computation times compared to heuristics [19]. The concept of SA was initially introduced by Metropolis et al. [20] and was first used by Kirkpatrick et al. [21] to solve a combinatorial optimization problem.

The remainder of the paper is organized as follows. Section 2 presents the literature review. Section 3 presents the problem definition. Section 4 describes the proposed simulated annealing heuristic for the CGVRP. Section 5 discusses the numerical experiments. Finally, Section 6 presents conclusions and the future direction.

#### 2. Literature Review

Transportation is an essential sector in connecting one location to another for the delivery or pickup of freight and/or passengers. However, it comes with hazardous and fatal emissions that severely impact global health. Therefore, “greening” road transportation has become a top-ranking necessity [24].

Transportation with a low negative impact on human health and the environment is defined as green transportation [25], and through realizing the importance and the capability of transportation to promote improvements in a company, the number of companies that are considering “greening” as a strategy is growing [26]. Reducing emissions can be initiated by an efficiently designed VRP. In recent years, the GVRP has considered special issues related to the VRP and the environment. Lin et al. [27] defined that the objective of GVRP is to be able to harmonize the economic and environmental costs by applying effective routes to meet financial indices and environmental concerns. Therefore, the vehicle route should be designed to satisfy green transportation requirements by reducing consumption levels and consequently reducing CO2 emissions from road transportation to minimize the total cost [24]. The GVRP research stream is divided into two groups. One group considers minimizing the fuel consumption [14, 22] and another considers the use of AFVs & alternative fuel stations (AFSs) [12].

According to [27], other notable research streams sharing the concept of the GVRP on minimizing the negative impact on the environment are the Pollution-Routing Problem (PRP) and the VRP in reverse logistics. Kara et al. [22] began studies on energy (i.e., fuel consumption) minimization by considering that transportation cost is affected by load and distance. They formulated the problem as an integer linear program and solved it exactly using CPLEX 8.0. Kuo [14] considered the minimization of total fuel consumption for the time-dependent vehicle routing problem (TDVRP). The transportation speed was used to calculate the fuel consumption in this model. Xiao et al. [18] also focused on minimizing fuel consumption but used a load-dependent function where the fuel consumption rate (FCR) is linearly associated with the vehicle’s load. A simulated annealing heuristic was used to solve the problem. Erdoğan and Miller-Hooks [12] considered the use of alternative fuel in the problem of recharging or refueling of the vehicle, the alternative fuel powered vehicle (AFV) in particular. They used a Modified Clarke and Wright Savings (MCWS) heuristic and the Density-Based Clustering Algorithm (DBCA) to solve the problem. Schneider et al. [23] examined the VRP time window when employing an electrical vehicle. They presented a hybrid heuristic that combines a variable neighborhood search (VNS) algorithm with a tabu search (TS) heuristic to solve their problem. The GVRP is one benchmark instance that demonstrates the performance of the VNS/TS heuristic. The VNS/TS heuristic outperforms both heuristic methods proposed by [12].

Fuel oil being used as the primary source of energy has, in fact, a side effect on the environment. Therefore, the use of alternative energy is intended to reduce the adverse effects of fossil fuels. In fact, there are many AFVs that are marketed and sold to consumers around the world. USA has 13,776 alternative fuel stations. Due to environmental concerns, an increasing oil demand, and rising oil production costs, efforts for cleaner alternative fuels and advanced propulsion technologies are the driving factors for developing AFVs. Furthermore, AFVs hold the potential to solve many environmental challenges that relate to emissions caused by transportation [28–30]. Table 1 lists the comparative research related to the GVRP from [12, 14, 18, 22, 23].

Due to the capacity limitation in each vehicle, some studies look at capacity as an extension of the VRP. Lin et al. [15] defined the CVRP as a variant of the VRP in which the vehicles have the maximal loading capacity to serve a set of customers with known demands in order to minimize the total cost. As defined by Breedan [17], the demand of the standard VRP is deterministic and known. In light of the properties of cost in the matrix distance, the CVRP can be further partitioned into the symmetric CVRP (SCVRP) and asymmetric CVRP (ACVRP) [31].

#### 3. Problem Definition

The development of transportation planning is essential to mitigating the negative impact of transportation on the environment. The CGVRP is an extension of the GVRP that considers the importance of vehicle capacity. Previously, Zhang et al. [32] proposed the CGVRP without considering the maximum time constraint in their model. In this paper, we develop a mixed integer linear program (MILP) model for the CGVRP that focuses on determining vehicle routes by considering the utilization of AFVs, vehicle capacity, and the maximum time constraint.

The CGVRP is expected to generate tours of at most* m* vehicles with the objective of minimizing the total distance traveled. Each vehicle makes one single tour, and each vehicle starts and ends at the depot, visiting a subset of customer vertices in between and, if needed, an AFS. The total demand in each tour (route) should not exceed the vehicle capacity. In addition, the total time spent for each tour (route) should not exceed the maximum route duration, .

Adopting the model from that of Erdoğan and Miller-Hooks [12], the CGVRP is an undirected graph with , where* V* is a vertex set that includes a set of customers , depot and a set of AFS with . Moreover, , where . The set of the edge is , which connects the vertices of* V*.

The CGVRP model considers a deterministic and static problem. It assumes that the vehicles have a constant travel speed and the refueling stations have unlimited capacity; hence, the vehicle tank is filled up to full capacity after visiting an AFS. An AFS can be visited more than once, so to allow for multiple visits to an AFS, a set of dummy AFSs is needed. The set of dummy AFSs is represented by , and then the set of vertices becomes . The technique of adding dummy vertices was introduced by Bard et al. [33]. The CGVRP formulation is as follows.

*Sets * Depot Set of customers, Set of depot and customers, Set of AFSs, Set of dummy AFSs created from Set of dummy AFSs and the depot, Set of vertices, Set of vertices with dummy vertices,

*Parameters * Fuel capacity of vehicle (gallons) Rate of vehicle fuel consumption (gallons per mile) Maximum route time Travel time from vertices* i* to* j* Service time at customer vertex* i *if (refueling time at the AFS vertex* i* if ) Demand at vertex* i*, set to 0 if Vehicle capacity *m* Maximum number of vehicles

*Decision Variables * 1 if there is travel from vertex* i* to* j*, 0 otherwise Remaining tank fuel capacity at vertex* j* Arrival time at vertex* j* Remaining capacity of vehicle at vertex* j*

*Mathematical Programming Model *

The objective function (1) minimizes the total distance of the routes. Constraint set (2) guarantees that each customer is visited exactly once. Constraint set (3) guarantees that each AFS is visited at most once. Constraint set (4) is the flow conservation constraint. Constraint sets (5) and (6) ensure that at most* m* vehicles are used. The arrival times at vertex* j* are shown in constraint set (7). Constraint sets (8) and (9) ensure the limit for route duration. The remaining fuel level at vertex* j* is computed in constraint sets (10) and (11). Constraint set (12) ensures that there is sufficient fuel up to the depot. Constraint set (13) shows the remaining vehicle capacity at vertex* j*. Constraint set (14) guarantees that the amount loaded at vertex* i* should not exceed the vehicle capacity. Binary and integrality are guaranteed with constraint set (15).

#### 4. Simulated Annealing for the CGVRP

We use an SA heuristic for the CGVRP. The SA is a local search heuristic that has the capability to escape from local optima by accepting a worse solution with a certain probability function [20]. In the following subsection, we discuss the proposed SA in detail, including the solution representation, initial solution, neighborhood, parameter setting, and the proposed SA procedure.

##### 4.1. Solution Representation

The CGVRP network involves a depot, a set of customer nodes, a set of AFS nodes, and a set of dummy AFS nodes. Figure 1 illustrates an example of the CGVRP solution. In this example, the AFVs have a fuel capacity of 60 gallons and a consumption rate of 0.2. The solution of this example generates 6 routes that are assigned in 6 AFVs. As shown in Figure 1, the first vehicle is assigned to visit customers 5 and 9. The vehicle requires more than 60 gallons to reach customer 5; thus, it decides to visit the AFS. The remaining vehicles are assigned in the same way as the first vehicle.

The CGVRP solution is represented by a permutation of* n* customers denoted by the set . These numbers represent either customer nodes or a combination of customer and AFS nodes. Figure 2 shows the solution representation of the example. In Figure 2, Route 1 illustrates the route of the first vehicle. It starts at the depot then visits the AFS, then customers 5 and 9, and then returns to the depot.

##### 4.2. Initial Solution

The initial solution is constructed using the Nearest Neighbor (NN) method [34]. Starting from the depot, the vehicle visits the nearest unvisited customer vertex. The process continues until all customers are visited, subject to fuel capacity, route duration, and vehicle capacity. This allows for the vehicle to visit the nearest AFS when needed. Finally, it returns to the depot.

##### 4.3. Neighborhood

We use the standard SA with a random neighborhood structure that features movement procedures, including swap and insertion. These moves are commonly embedded in the SA heuristic or other metaheuristics frameworks [35]. In our proposed SA, we use* swap*,* insertion*,* insert_AFS*, and* delete_AFS*. These SA neighborhood structures are described in Table 2. The probabilities of selecting any of the four neighborhood moves are predetermined to be equal; that is, each move has a probability of 0.25 in every step of the algorithm.

##### 4.4. Parameters Used

The SA heuristic uses five parameters,* I*_{iter},* T*_{0},* T*_{F},* N*_{non -improving}, and *α*.* I*_{iter} denotes the number of iterations for which the search proceeds at a particular temperature, while* T*_{0} represents the initial temperature, and* T*_{F} represents the final temperature, below which the SA procedure is stopped.* N*_{non -improving} is the maximum allowable number of reductions in temperature, during which the best objective function value is not improved. Finally, *α* is the coefficient controlling the cooling schedule from Kirkpatrick et al. [21].

##### 4.5. SA Procedure

The SA heuristic consists of two parts: initial phase and improvement phase. The initial phase employs the Nearest Neighbor (NN) algorithm to create an initial feasible solution. After finishing the initial phase, the algorithm continues to the improvement phase. The improvement phase tries to improve the solution from the initial phase using the neighborhood structures and local search.

The algorithm starts by setting the current temperature and randomly generating an initial solution* X*. The current best solution,* X*_{best}, is set to be* X*, and the best objective function of* X*, denoted by , is set to be Obj(*X*). In each iteration, a new solution* Y* is obtained from a predefined neighborhood of the current solution* X*. After that, the objective function values of* X* and* Y* are evaluated using the formula, Δ = obj(*Y*)-obj(*X*).* Y* is better than* X* if Δ is negative, in which case* Y* will replace* X*. Otherwise,* Y* will replace* X* with a probability equal to exp (-Δ/*T*). Furthermore, the best solution obtained is set to be and .

The current temperature* T*_{0} is decreased after iterations according to the formula . Furthermore, is improved by a local search procedure. All possible* swap* moves to are used in the local search procedure. The best solution obtained from the local search is used to replace if it is better than .

The algorithm will terminate if one of the terminating conditions is reached. In this algorithm, there are two terminating conditions: (1) when the current temperature* T*_{0} equals and (2) when no improvement is attained in consecutive temperature reductions. Figure 3 shows the flowchart of the proposed SA heuristic.

#### 5. Numerical Experiment

A numerical experiment is conducted to demonstrate the problem and evaluate the performance of the proposed SA algorithm. We extend the instances from that of Erdoğan and Miller-Hooks [12] and solve the instances using CPLEX. The SA heuristic is coded in C++, and the experiment is run on a desktop computer using an Intel Core 2 Duo E8400 processor, running on a 32-bit platform with a 3.00 GHz processor speed using 4 GB of RAM under the Windows 7 operating system.

##### 5.1. Parameter Setting

We performed a parameter setting to determine the best parameter set to be used for the proposed SA algorithm.

The benchmark instances used in the parameter setting are the small instances taken from Erdoğan and Miller-Hooks [12]. These include scenarios denoted as S1, S2, S3, and S4, as shown in Table 3. A sample instance is randomly selected from the problem set of each of the four scenarios for preliminary testing. Then, a two-level () factorial design is conducted. The low level and high level values are obtained in a pilot experiment.

The SA parameters initially used for the pilot experiment are as follows: =50000; =1000; =0.92; =1; and =30. The next step is to perform the sensitivity analysis using the one-at-a-time (OAT) method, proposed by Morris [36], to identify the low and high level by comparing the SA results to the CPLEX result. The input dataset for the SA is generated five times for each problem to choose the best two values (average), namely, the high level and low level values. A factorial design experiment is also conducted to identify the best parameter setting for each SA parameter. The lowest difference in the heuristic result and CPLEX result is chosen as the best parameter. These parameters are used for all scenarios: =70000; =900; =0.92; =0,1; and =400.

##### 5.2. Performance of SA: Benchmark Instances

To evaluate the solution quality of the SA on the CGVRP, we evaluate the SA performance on the GVRP instance proposed by Erdoğan and Miller-Hooks [12] because there is currently no benchmark instance similar to the CGVRP. The benchmark method used to verify the proposed SA is from Erdoğan and Miller-Hooks [12] and Schneider et al. [23]. Using a hybrid of variable neighborhood search and tabu search (VNS/TS), Schneider et al. [23] employed the benchmark instance from Erdoğan and Miller-Hooks [12] to evaluate the quality of their proposed algorithm in the GVRP.

The percentage difference is computed between the SA and CPLEX results to assess the algorithm quality as follows: . The lower the value of is, the better the SA performance will be. The same conditions also work for the DBCA Diff column, MCWS Diff column, and VNS/TS Diff column.

In comparing the SA performance with the benchmark methods, namely, the MCWS heuristic, DBCA heuristic, and VNS/TS, we conducted numerical experiments for the same scenarios. Table 3 shows the instances used.

Table 4 and Figure 4 show the best solution obtained in five runs of the proposed SA for the small instance. It also shows the SA’s superior performance at an average* Diff* of 0.47% compared to that of the MCWS at 13.1%, the DBCA at 12.6%, and the VNS/TS at 4.57%.

For a large instance, Table 5 shows the SA’s performance at 10432 (in distance), which is the lowest compared to the MCWS at 12216.9, the DBCA at 16028.68, and the VNS/TS at 10580.69. Figure 5 shows that the SA has the lowest total cost in miles. Therefore, the SA outperforms the MCWS, DBCA, and VNS/TS in both average* Diff *and distance.

##### 5.3. Experiment on the CGVRP

The experiment on the GVRP is used to evaluate the solution obtained through the proposed SA in the CGVRP by using the same parameter settings that were evaluated in the GVRP instance. Because this research is the first to study the CGVRP, there is no benchmark instance to evaluate the solution of the problem. Therefore, we design an extension instance based on Erdoğan and Miller-Hooks [12], both for small and large instance problems.

In the CGVRP experiment, using the data for the small instance from [12], we extend the following: the customer demand is generated by a random number around 1-10, the vehicle capacity is 500, the fuel tank capacity is 60 gallons, the fuel consumption needed per mile per hour (mph) is 0.2 gallons, the total tour duration for each vehicle is 11 hours, the customer service time is 30 minutes, and the AFS service time is 15 minutes. The customer and AFS locations are similar to those in [12].

In the exact solution using CPLEX, some small instances are relaxed by determining the number of routes in each instance. The exact solutions (in the bold figure) are the optimal solutions obtained. To evaluate the quality of the proposed SA algorithm for the CGVRP problem, we calculated the difference relative to that of the CPLEX. Table 6 shows the SA’s average* Diff* result at 0.48%. Figure 6 shows the total cost obtained by the SA compared to the exact solution. The SA obtains the same solution as the exact solution in 60% of the total instances.

We conducted the numerical experiment for the large instance of the CGVRP. The results between the capacitated and uncapacitated vehicles in the GVRP are evaluated for the Large1 scenario, as explained in Table 3. The noncapacitated vehicle is represented by using a huge vehicle capacity (6000), and therefore the vehicle capacity constraint will not restrict the solution.

The result in Table 7 shows that the total cost of using the low capacity (volume capacity at 25) is 6161.8, which is 32.5% higher than uncapacitated GVRP. This is caused by the increasing number of tours needed when the vehicle capacity is limited. In this case, the number of tours for the CGVRP increases by 38.8% versus the uncapacitated GVRP. This result shows that the decision in choosing an appropriate capacity to fulfill customer demand is crucial.

For the CGVRP instance 111c21s, which has a huge vehicle capacity, the total cost obtained is similar to the total cost of the GVRP as shown in Table 5. Table 5 shows that the SA outperforms the MCWS, DBCA, and VNS/TS. Therefore, a vehicle capacity of 6,000 is used for the sensitivity analysis.

##### 5.4. Sensitivity Analysis

A sensitivity analysis is used to investigate the capability of the SA in solving large instances. It consists of four scenarios as described in Table 8. Tables 8–11 present the sensitivity analysis results.

Table 9 shows the impact on the total distance when the vehicle capacity varies. The distance traveled is reduced by 3.4% and 3.5% with an increase in the vehicle capacity from 1000 to 2000 and 2000 to 3000, respectively. However, increasing the vehicle capacity above 3000 does not improve the solution.

Table 10 shows that the number of unserved customers increases as the number of customers increases from 200 to 500. However, the result shows that the proposed algorithm can handle large-size problems and can find the lowest cost solution compared to the previous algorithms, as shown in Table 5.

Table 11 shows the impact of increasing the AFS availability. The AFS increases from 22 to 28 (in increments of 2) by keeping the basic number of AFSs. Increasing the AFS availability from 22 to 24 reduces the total travel distance by 1%. However, increasing the number of AFSs to more than 24 does not improve the solution.

Table 12 shows that by increasing the driving range from 200 to 500 miles (in increments of 100 miles), the number of customers served also increases. Moreover, all customers are served when using 500 miles. A summary of the sensitivity analysis result on total cost can be seen in Figure 7.

**(a)**

**(b)**

**(c)**

**(d)**

#### 6. Conclusion

This research formulates the CGVRP model as an extension of the GVRP using the simulated annealing (SA) algorithm to solve it. The proposed algorithm, simulated annealing, can be a reliable method for the GVRP solution. The algorithm shows superiority over the previous methods implemented by Erdoğan and Miller-Hooks [12].

In the CGVRP, the SA algorithm performs well both in the small and large instances. On average, the SA heuristic evaluates the fitness value 110 times to find a solution. Moreover, the sensitivity analysis based on a large problem instance shows that the number of customers influences the minimum number of vehicles that need to be used. Increasing the vehicle capacity from 1000 to 2000 and the vehicle driving range from 200 to 300 influences the total travel distance. The SA algorithm in this research can be used to help companies in evaluating the possible number of customers that can be served, the fleet size needed, and the travel distance requirement as a limitation of the fueling stops and driving range.

Future research directions may focus on considering heterogeneous fleets, different driving range limitations, fuel prices, or use of different fuel sources, including the type of fuel and the amount of fuel.

#### Data Availability

The basic GVRP instances data used to support the findings of this study were supplied by Sevgi Erdoğan under license and so cannot be made freely available. Requests for access to these data should be made to Sevgi Erdoğan (serdogan@umd.edu; Tel.: +1 301 405 2046; fax: +1 301 405 2585).

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was funded by Badan Penerbit dan Publikasi (BPP) of Gadjah Mada University and the Ministry of Science and Technology of the Republic of China (Taiwan) under grant MOST 106-2410-H-011-002-MY3. This research is based on Mr. Candra Bachtiyar’s thesis.