Mathematical Problems in Engineering
Volume 2019 (2019), Article ID 2358258, 18 pages
https://doi.org/10.1155/2019/2358258

Research Article

A Simulated Annealing Heuristic for the Capacitated Green Vehicle Routing Problem

Nur Mayke Eka Normasari,1,2 Vincent F. Yu,¹ Candra Bachtiyar,¹ and Sukoyo³

¹Department of Industrial Management, National Taiwan University of Science and Technology, 43, Section 4, Keelung Road, Taipei 10607, Taiwan
²Department of Mechanical and Industrial Engineering, Gadjah Mada University, Jl. Grafika No. 2, Yogyakarta 55281, Indonesia
³Department of Industrial Engineering and Management, Institut Teknologi Bandung, Jl. Ganesha 10, Bandung, Jawa Barat, Indonesia

Correspondence should be addressed to Nur Mayke Eka Normasari; mayke@ugm.ac.id

Received 28 May 2018; Revised 31 October 2018; Accepted 19 December 2018; Published 13 January 2019

Academic Editor: Sitek Paweł

Copyright © 2019 Nur Mayke Eka Normasari et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

Table 1: Summary of GVRP research.

No.
Aspect of review
Kara et al. [22]
Kuo [14]
Xiao et al. [18]
Erdoğan and Miller-Hooks [12]
Schneider et al. [23]
This research (2014)

1
Topic
Energy minimization
Fuel consumption
Fuel consumption
Using AFV & AFS
Using electrical vehicle
Using AFV & AFS

2
Objective
Minimize energy and distance
Minimize fuel consumption
Minimize fuel cost
Minimize distance
Minimize distance
Minimize distance

3
Decision variable
Traveling decision from vertex i to j and the weight of the vehicle if it goes from i to j
Traveling decision from vertex i to j
Traveling decision from vertex i to j and the weight of the (cargo) vehicle if it goes from i to j
Traveling decision from vertex i to j, fuel level, time
Traveling decision from vertex i to j, time, remaining, cargo, battery capacity
Traveling decision from vertex i to j, fuel level, time, capacity

4
Method/ solver
CPLEX 8.0
Simulated Annealing
Simulated Annealing
MCWS & DBCA
VNS/TS
Simulated Annealing

5
Number of depot
Single
Single
Single
Single
Single, multiple
Single

6
Capacity of vehicle
Available
Available
Available
Not available
Available
Available

7
Demand
Available
Available
Available
Not available
Available
Available

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.

Figure 1: An example of the CGVRP solution.

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.

Figure 2: Solution representation.

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.

Table 2: SA neighborhood structure.

Neighborhood structure
Description

swap
Performed by selecting a random customer in a route and another customer from another route to swap their positions

insertion
Conducted by randomly choosing a customer in a route and another customer from a different route. Then, the first chosen customer is deleted and inserted into the front of the second chosen customer from another route

insert_AFS
An AFS node is inserted into the infeasible route to make it a feasible one. This is a way to remedy infeasible routes due to fuel capacity

delete_AFS
Performed when the removal of an AFS node can reduce the cost incurred in the system while preserving the feasibility of the route

4.4. Parameters Used

The SA heuristic uses five parameters, Iiter, T₀, TF, Nnon -improving, and α. Iiter denotes the number of iterations for which the search proceeds at a particular temperature, while T₀ represents the initial temperature, and TF represents the final temperature, below which the SA procedure is stopped. Nnon -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, Xbest, 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₀ 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₀ equals and (2) when no improvement is attained in consecutive temperature reductions. Figure 3 shows the flowchart of the proposed SA heuristic.

Figure 3: Flow chart of the simulated annealing algorithm for the CGVRP.

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.

Table 3: Scenarios for benchmark instances from Erdoğan and Miller-Hooks [12].

Scenario
Description
Details

S1
Uniform customer configuration impact
10 instances of 20 uniform customer locations and 3 fixed AFS locations

S2
Clustered customer configuration impact
10 instances of 20 clustered customer locations with 3 fixed AFS locations

S3
AFS configuration impact
10 instances, half from S1 and half from S2, with 6 random AFS locations

S4
Station density impact
10 instances, half from 1 instance of S1 and half from 1 instance of S2,
with increasing the AFSs in increments of 2 from 2 to 10

Large 1
Transitioning to AFV
111 customers with 21 AFSs

Large 2
Increasing the number of customers impacted
Increasing customer number from 200 to 500

Large 3
Increasing AFS availability impact
Customers are fixed at 111, and AFSs are increased from 22 to 28 in increments of 2

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

Table 4: Performance comparison on small GVRP instances.

Sample
CPLEX
MCWS
DBCA
VNS/TS
SA

E S (miles)
N
C
N
C
T C
Diff (%)
T C
Diff (%)
T C
N
C
Time
Diff (%)
T C
N
C
Time
Diff

(miles)
(miles)
(miles)
(s)
(miles)
(s)
(%)

20c3sU1
1632.0
5
20
6
20
1818.4
11.4
1797.5
10.1
1797.5
6
20
41.4
10.1
1632.02
5
20
11.71
0.00

20c3sU2
1542.6
6
20
6
20
1614.2
4.6
1613.5
4.6
1574.8
6
20
38.4
2.1
1574.13
6
20
12.36
2.04

20c3sU3
1624.7
6
20
7
20
1969.6
21.2
1964.6
20.9
1704.5
6
20
38.4
4.9
1624.73
6
20
11.63
0.00

20c3sU4
1462.9
5
20
5
20
1508.4
3.1
1487.2
1.7
1482.0
5
20
39.0
1.3
1462.89
5
20
11.56
0.00

20c3sU5
1661.3
6
20
6
20
1752.7
5.5
1752.7
5.5
1689.4
6
20
40.2
1.7
1661.29
6
20
11.84
0.00

20c3sU6
1537.1
5
20
6
20
1668.2
8.5
1668.2
8.5
1618.7
6
20
40.2
5.3
1537.07
5
20
11.74
0.00

20c3sU7
1586.4
6
20
6
20
1730.5
9.1
1730.5
9.1
1713.7
6
20
38.4
8
1586.42
6
20
11.84
0.00

20c3sU8
1705.8
6
20
6
20
1718.7
0.8
1718.7
0.8
1706.5
6
20
40.2
0
1705.82
6
20
11.50
0.00

20c3sU9
1708.1
6
20
6
20
1714.4
0.4
1714.4
0.4
1708.8
6
20
39.6
0
1708.1
6
20
11.96
0.00

20c3sU10
1180.9
4
20
5
20
1309.5
10.9
1309.5
10.9
1181.3
4
20
38.4
0
1212.24
5
20
11.75
2.65

20c3sC1
1173.1
4
20
5
20
1300.6
10.9
1300.6
10.9
1173.6
4
20
37.2
0
1173.1
4
20
11.50
0.00

20c3sC2
1539.3
5
19
5
19
1553.5
0.9
1553.5
0.9
1540.0
5
19
34.8
0
1567.16
6
19
12.10
1.81

20c3sC3
879.9
3
12
4
12
1083.1
23.1
1083.1
23.1
880.2
3
12
15.0
0
879.862
3
12
9.78
0.00

20c3sC4
978.4
4
18
5
18
1135.9
16.1
1091.8
11.6
1059.4
4
18
31.8
8.3
978.428
4
18
11.59
0.00

20c3sC5
2189.8
7
19
7
19
2190.7
0.0
2190.7
0.0
2156.0
7
19
36.0
1.5
2138.69
7
19
12.64
2.33

20c3sC6
2332.4
7
17
9
17
2883.7
23.6
2883.7
23.6
2758.2
8
17
42.6
18.3
2332.4
7
17
12.10
0.00

20c3sC7
1055.0
3
6
5
6
1701.4
61.3
1701.4
61.3
1394.0
4
6
10.8
32.1
1054.98
3
6
8.28
0.00

20c3sC8
3139.7
9
18
10
18
3319.7
5.7
3319.7
5.7
3139.7
9
18
37.2
0
3235.57
10
18
12.96
3.05

20c3sC9
1643.6
5
19
6
19
1811.1
10.2
1811.1
10.2
1799.9
6
19
36.0
9.5
1687.43
6
19
11.94
2.67

20c3sC10
2583.4
8
15
8
15
2648.8
2.5
2644.1
2.4
2583.4
8
15
27.0
0
2571.24
9
15
11.56
0.47

S1_2i6s
1567.4
6
20
6
20
1614.2
3.0
1614.2
3.0
1578.1
6
20
42.6
0.7
1567.39
6
20
12.95
0.00

S1_4i6s
1396.7
5
20
5
20
1561.3
11.8
1541.5
10.4
1397.3
5
20
45.0
0
1396.69
5
20
12.05
0.00

S1_6i6s
1519.7
5
20
6
20
1616.2
6.4
1616.2
6.4
1560.5
5
20
43.8
2.7
1519.7
5
20
12.18
0.00

S1_8i6s
1623.6
6
20
6
20
1902.5
17.2
1882.5
16.0
1692.3
6
20
44.4
4.2
1623.57
6
20
12.31
0.00

S1_10i6s
1173.0
4
20
5
20
1309.5
11.6
1309.5
11.6
1173.5
4
20
42.6
0
1173.04
4
20
12.41
0.00

S2_2i6s
1633.2
6
20
6
20
1645.8
0.8
1645.8
0.8
1633.1
6
20
45.0
0
1645.1
6
20
12.28
0.73

S2_4i6s
1498.5
5
19
6
19
1505.1
0.4
1505.1
0.4
1533.0
5
19
52.8
2.3
1504.42
6
19
12.26
0.40

S2_6i6s
2626.0
8
20
10
20
3115.1
18.6
3115.1
18.6
2431.3
7
20
46.8
7.4
2626.46
8
20
12.97
0.02

S2_8i6s
2158.4
7
16
9
16
2722.6
26.1
2722.6
26.1
2158.4
7
16
34.2
0
2174.85
7
16
11.41
0.76

S2_10i6s
1486.1
6
16
6
16
1995.6
34.3
1995.6
34.3
1958.5
6
17
36.6
31.8
1511.38
6
16
11.46
1.70

S1_4i2s
1581.6
6
20
6
20
1582.2
0.0
1582.2
0.0
1582.2
6
20
37.8
0
1581.56
6
20
12.23
0.00

S1_4i4s
1459.5
5
20
6
20
1580.5
8.3
1580.5
8.3
1460.1
5
20
40.8
0
1459.49
5
20
12.09
0.00

S1_4i6s
1396.7
5
20
5
20
1561.3
11.8
1541.5
10.4
1397.3
5
20
45.0
0
1396.69
5
20
11.99
0.00

S1_4i8s
1396.7
5
20
6
20
1561.3
11.8
1561.3
11.8
1397.3
6
20
49.2
0
1396.69
5
20
12.33
0.00

S1_4i10s
1395.5
5
20
5
20
1536.0
10.1
1529.7
9.6
1396.0
5
20
51.0
0
1395.16
5
20
12.00
0.02

S2_4i2s
978.4
4
18
5
18
1135.9
16.1
1117.3
14.2
1059.4
4
18
30.6
8.3
978.428
4
18
11.46
0.00

S2_4i4s
1365.0
5
19
6
19
1522.7
11.6
1522.7
11.6
1446.1
5
19
36.0
5.9
1365.08
5
19
12.78
0.01

S2_4i6s
1353.0
5
20
6
20
1786.2
32.0
1730.5
27.9
1434.1
5
20
41.4
6
1353.03
5
20
12.16
0.00

S2_4i8s
1353.0
5
20
6
20
1786.2
32.0
1786.2
32.0
1434.1
5
20
45.0
6
1353.03
5
20
12.13
0.00

S2_4i10s
1353.1
5
20
6
20
1783.6
31.8
1729.5
27.8
1413.1
5
20
46.8
4.4
1353.02
5
20
12.11
0.01

Average

13.1

12.6

39.0
4.57

11.90
0.47

E S = exact solution, N = number of tours, C = customers served, T C = total cost, Diff = absolute difference with CPLEX.

Figure 4: Performance comparison on small GVRP instances in % average Diff.

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.

Table 5: Performance comparison for large GVRP instances.

Instance
MCWS
DBCA
VNS/TS
SA

T C (miles)
N
C S
Time (s)
T C (miles)
N
C S
Time (s)
T C (miles)
N
C S
Time (s)
T C (miles)
N
C S
Time (s)

111c21s
5626.64
20
109
NA
5626.64
20
109
NA
4797.15
17
109
1305.6
4648.24
18
109
41.40

200c21s
10428.6
35
190
NA
10413.59
36
191
NA
8963.46
35
192
4599
8847.24
35
192
75.38

300c21s
14242.6
49
281
NA
14229.92
49
282
NA
12594.77
46
283
10934
12334.6
49
283
123.96

350c21s
16471.5
57
329
NA
16460.3
57
329
NA
14323.02
51
329
13922
14326.6
56
329
128.26

400c21s
19472.1
67
378
NA
19099.04
66
373
NA
16850.21
61
378
18307
16561
69
378
154.04

450c21s
21854.2
75
424
NA
21854.19
75
424
NA
18521.23
68
424
31531
18497.2
76
424
176.69

500c21s
24527.5
84
471
NA
24517.08
84
471
NA
21170.9
76
471
21361
20887.8
83
471
203.42

111c22s
5610.57
20
109
NA

4802.16
17
109
1413.6
4646.8
18
109
41.72

111c24s
5412.48
20
109
NA

4786.96
17
109
1314
4679.73
18
109
41.24

111c26s
5408.38
20
109
NA

4778.62
17
109
1507.2
4671.94
18
109
40.79

111c28s
5331.93
20
109
NA

4799.15
17
109
1450.2
4650.9
18
109
39.45

Average
12216.9

NA
16028.68

NA
10580.69

9785.84
10432.00

96.94

T C = total cost, N = number of tours, C S = customers served.

Figure 5: Performance comparison for large GVRP instances in total cost (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.

Table 6: SA performance for small CGVRP instances.

Sample
CPLEX
SA

Exact solution
Number of tours
Customers served
Total cost
Number of tours
Customers served
Time
Absolute diff CPLEX

(miles)
(miles)
(s)
(%)

20c3sUCAP1
1632.01
5
20
1632.02
5
20
12.98
0.00

20c3sUCAP2
1542.61
6
20
1574.13
6
20
13.67
2.04

20c3sUCAP3
1624.73
6
20
1624.73
6
20
12.51
0.00

20c3sUCAP4
1462.89
5
20
1462.89
5
20
11.99
0.00

20c3sUCAP5
1661.29
6
20
1661.29
6
20
11.88
0.00

20c3sUCAP6
1537.07
5
20
1537.07
5
20
12.84
0.00

20c3sUCAP7
1586.42
6
20
1586.42
6
20
12.65
0.00

20c3sUCAP8
1705.82
6
20
1705.82
6
20
12.14
0.00

20c3sUCAP9
1708.1
6
20
1708.1
6
20
12.87
0.00

20c3sUCAP10
1180.86
4
20
1212.24
5
20
12.63
2.66

20c3sCCAP1
1173.09
4
20
1173.1
4
20
12.68
0.00

20c3sCCAP2
1539.34
5
19
1567.16
6
19
12.84
1.81

20c3sCCAP3
879.863
3
12
879.862
3
12
9.76
0.00

20c3sCCAP4
978.428
4
18
978.428
4
18
12.17
0.00

20c3sCCAP5
2189.77
7
19
2138.69
7
19
12.90
2.33

20c3sCCAP6
2332.4
7
17
2332.4
7
17
11.83
0.00

20c3sCCAP7
1054.98
3
6
1054.98
3
6
8.00
0.00

20c3sCCAP8
3139.92
9
18
3117.55
9
18
13.59
0.71

20c3sCCAP9
1643.55
5
19
1687.43
6
19
12.27
2.67

20c3sCCAP10
2583.42
8
15
2571.24
9
15
12.44
0.47

S1_2i6sCAP
1567.39
6
20
1567.39
6
20
12.52
0.00

S1_4i6sCAP
1396.69
5
20
1396.69
5
20
13.42
0.00

S1_6i6sCAP
1519.7
5
20
1519.7
5
20
12.93
0.00

S1_8i6sCAP
1623.57
6
20
1623.57
6
20
12.68
0.00

S1_10i6sCAP
1173.04
4
20
1173.04
4
20
13.05
0.00

S2_2i6sCAP
1633.17
6
20
1645.88
6
20
12.73
0.78

S2_4i6sCAP
1498.5
5
19
1504.42
6
19
12.74
0.40

S2_6i6sCAP
2625.97
8
20
2626.46
8
20
13.54
0.02

S2_8i6sCAP
2158.35
7
16
2174.85
7
16
11.26
0.76

S2_10i6sCAP
1486.11
6
16
1552.78
6
16
12.35
4.49

S1_4i2sCAP
1581.56
6
20
1581.56
6
20
12.74
0.00

S1_4i4sCAP
1460.09
5
20
1459.49
5
20
12.75
0.04

S1_4i6sCAP
1396.69
5
20
1396.69
5
20
12.98
0.00

S1_4i8sCAP
1396.69
5
20
1396.69
5
20
12.74
0.00

S1_4i10sCAP
1395.46
5
20
1395.16
5
20
12.98
0.02

S2_4i2sCAP
978.428
4
18
978.428
4
18
11.96
0.00

S2_4i4sCAP
1364.98
5
19
1365.08
5
19
12.88
0.01

S2_4i6sCAP
1353.03
5
20
1353.03
5
20
12.22
0.00

S2_4i8sCAP
1353.03
5
20
1353.03
5
20
12.65
0.00

S2_4i10sCAP
1353.13
5
20
1353.02
5
20
12.30
0.01

Average

12.45
0.48

a indicates a nonoptimal solution because the model is relaxed by determining the number of tours.

Figure 6: SA performance for small CGVRP instances in total cost (distance).

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.

Table 7: SA performance for large CGVRP instances.

Scenario
Instance
SA

Total cost
Time
Number of tours
Customers served
Vehicle capacity

(miles)
(s)

Large1 capacitated
111c21s
6161.8
40.475
25
109
25

Large1 uncapacitated
111c21s
4648.96
37.568
18
109
6000

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 8: Scenarios for sensitivity analysis.

Scenario
Description
Details

CGVRP A1
Vehicle capacity impact
Using 111 customers, 21 AFSs and vehicle capacity increase from 1,000 to 6,000 in 1,000 increments

CGVRP A2
Increasing the number of customers impact
Increase the number of customers from 200 to 500 in increments of 100; add random customer locations; keep AFS locations fixed and vehicle capacity at 6,000

CGVRP A3
Increasing AFS availability impact
Similar to GVVRP A1, with adding AFS locations, it
increased in increments of 2 from 22 to 28 and vehicle capacity at 6000

CGVRP A4
Driving range limit impact
Similar to GVVRP A1, driving range increases in 100 mile increments from 200 miles to 500 miles and vehicle capacity at 6,000

Table 9: Scenario CGVRP A1: effect of vehicle capacity.

Scenario
Instance
SA

Total cost
Time
Number of tours
Customers served
Vehicle capacity

(miles)
(s)

CGVRP A1
111c21s
5041.78
39.659
23
109
1000

111c21s
4874.67
39.562
20
109
2000

111c21s
4710.18
39.164
19
109
3000

111c21s
4692.74
39.485
19
109
4000

111c21s
4673.66
38.823
18
109
5000

111c21s
4648.96
37.568
18
109
6000

Table 10: Scenario CGVRP A2: effect of increasing the number of customers.

Scenario
Instance
SA

Total cost (miles)
Time (s)
Number of tours
Customers served

CGVRP A2
200c21s
8786.25
71.596
35
192

300c21s
12691.1
112.058
46
283

400c21s
16680.3
156.797
67
378

500c21s
21126.8
200.556
84
471

Table 11: Scenario CGVRP A3: effect of the number of AFSs.

Scenario
Instance
Number of AFS
SA

Total distance (miles)
Time (s)
Number of tours
Customers served

CGVRP A3
111c22s
22
4711.73
37.171
19
109

111c24s
24
4651.46
38.043
19
109

111c26s
26
4650.12
38.28
19
109

111c28s
28
4650.47
37.84
19
109

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.

Table 12: Scenario CGVRP A4: effect of the vehicle driving range.

Scenario
Instance
SA

Total cost
Time (s)
Number of tours
Customers served
Driving range (miles)

(miles)

CGVRP A4
111c21s
4614.46
43.24
28
100
200

111c21s
4798.05
38.31
18
109
300

111c21s
4597.18
40.524
17
109
400

111c21s
5342.77
39.31
18
111
500

(a)

(b)

(c)

(d)

Figure 7: Impact on total cost by (a) vehicle capacity, (b) number of customers, (c) number of AFSs, and (d) driving range.

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.

References

S. K. Srivastava, “Green supply-chain management: a state-of-the-art literature review,” International Journal of Management Reviews, vol. 9, no. 1, pp. 53–80, 2007.

H.-J. Wu and S. C. Dunn, “Environmentally responsible logistics systems,” International Journal of Physical Distribution & Logistics Management, vol. 25, no. 2, pp. 20–38, 1995.

Q. Zhu and J. Sarkis, “Relationships between operational practices and performance among early adopters of green supply chain management practices in Chinese manufacturing enterprises,” Journal of Operations Management, vol. 22, no. 3, pp. 265–289, 2004.

J. Q. Frota Neto, J. M. Bloemhof-Ruwaard, J. A. E. E. van Nunen, and E. van Heck, “Designing and evaluating sustainable logistics networks,” International Journal of Production Economics, vol. 111, no. 2, pp. 195–208, 2008.

S. Ubeda, F. J. Arcelus, and J. Faulin, “Green logistics at eroski: a case study,” International Journal of Production Economics, vol. 131, no. 1, pp. 44–51, 2011.

I. T. Forum, “Reducing Transport Greenhouse Gas Emissions,” Trends data: OECD/ITF, 2010.

J. R. Stock and D. Lambert, Strategic Logistic Management, New York: Mc Graw Hill, 4th edition, 2001.

D. J. Bowersox, D. J. Closs, and M. B. Coope, Supply Chain Logistics, McGraw-Hill/Irwin, 3rd edition, 2010.

C. R. Bhat, S. Sen, and N. Eluru, “The impact of demographics, built environment attributes, vehicle characteristics, and gasoline prices on household vehicle holdings and use,” Transportation Research Part B: Methodological, vol. 43, no. 1, pp. 1–18, 2009.

X. Ou, X. Zhang, and S. Chang, “Scenario analysis on alternative fuel/vehicle for China's future road transport: Life-cycle energy demand and GHG emissions,” Energy Policy, vol. 38, no. 8, pp. 3943–3956, 2010.

M. Melaina and J. Bremson, “Refueling availability for alternative fuel vehicle markets: Sufficient urban station coverage,” Energy Policy, vol. 36, no. 8, pp. 3233–3241, 2008.

S. Erdoĝan and E. Miller-Hooks, “A green vehicle routing problem,” Transportation Research Part E: Logistics and Transportation Review, vol. 48, no. 1, pp. 100–114, 2012.

A. S. Alfa, S. S. Heragu, and M. Chen, “A 3-opt based simulated annealing algorithm for vehicle routing problems,” Computers & Industrial Engineering, vol. 21, no. 1–4, pp. 635–639, 1991.

Y. Kuo, “Using simulated annealing to minimize fuel consumption for the time-dependent vehicle routing problem,” Computers & Industrial Engineering, vol. 59, no. 1, pp. 157–165, 2010.

S.-W. Lin, K.-C. Ying, Z.-J. Lee, and F.-H. Hsi, “Applying simulated annealing approach for capacitated vehicle routing problems,” in Proceedings of the Systems, Man and Cybernetics, 2006. SMC'06. IEEE International Conference, vol. 1, pp. 639–644, 2006.

S.-W. Lin, V. F. Yu, and C.-C. Lu, “A simulated annealing heuristic for the truck and trailer routing problem with time windows,” Expert Systems with Applications, vol. 38, no. 12, pp. 15244–15252, 2011.

A. Van Breedam, “Improvement heuristics for the vehicle routing problem based on simulated annealing,” European Journal of Operational Research, vol. 86, no. 3, pp. 480–490, 1995.

Y. Xiao, Q. Zhao, I. Kaku, and Y. Xu, “Development of a fuel consumption optimization model for the capacitated vehicle routing problem,” Computers & Operations Research, vol. 39, no. 7, pp. 1419–1431, 2012.

J.-F. Cordeau, M. Gendreau, G. Laporte, J.-Y. Potvin, and F. Semet, “A guide to vehicle routing heuristics,” Journal of the Operational Research Society, vol. 53, no. 5, pp. 512–522, 2002.

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” The Journal of Chemical Physics, vol. 21, no. 6, pp. 1087–1092, 1953.

S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science, vol. 220, no. 4598, pp. 671–680, 1983.

I. Kara, B. Y. Kara, and M. K. Yetis, “Energy minimizing vehicle routing problem,” Combinatorial Optimization and Applications, pp. 62–71, 2007.

M. Schneider, A. Stenger, and D. Goeke, “The electric vehicle-routing problem with time windows and recharging stations,” Transportation Science, vol. 48, no. 4, pp. 500–520, 2014.

K. Salimifard, H. Shahbandarzadeh, and R. Raeesi, “Green transportation and the role of operation research,” in Proceedings of the 2012 International Conference on Traffic and Transportation Engineering, Singapore, 2012.

M. Björklund, “Influence from the business environment on environmental purchasing—drivers and hinders of purchasing green transportation services,” Journal of Purchasing and Supply Management, vol. 17, no. 1, pp. 11–22, 2011.

H. Min and I. Kim, “Green supply chain research: Past, present, and future,” Logistics Research, vol. 4, no. 1-2, pp. 39–47, 2012.

C. Lin, K. L. Choy, G. T. S. Ho, S. H. Chung, and H. Y. Lam, “Survey of green vehicle routing problem: past and future trends,” Expert Systems with Applications, vol. 41, no. 4, pp. 1118–1138, 2014.

A. R. Campbell, T. Ryley, and R. Thring, “Identifying the early adopters of alternative fuel vehicles: A case study of Birmingham, United Kingdom,” Transportation Research Part A: Policy and Practice, vol. 46, no. 8, pp. 1318–1327, 2012.

E. Graham-Rowe, B. Gardner, C. Abraham et al., “Mainstream consumers driving plug-in battery-electric and plug-in hybrid electric cars: A qualitative analysis of responses and evaluations,” Transportation Research Part A: Policy and Practice, vol. 46, no. 1, pp. 140–153, 2012.

B. van Bree, G. P. J. Verbong, and G. J. Kramer, “A multi-level perspective on the introduction of hydrogen and battery-electric vehicles,” Technological Forecasting & Social Change, vol. 77, no. 4, pp. 529–540, 2010.

P. Toth and D. Vigo, The Vehicle Routing Problem, Siam, 2002.

S. Zhang, Y. Gajpal, and S. S. Appadoo, “A meta-heuristic for capacitated green vehicle routing problem,” Annals of Operations Research, pp. 1–19, 2017.

J. F. Bard, L. Huang, M. Dror, and P. Jaillet, “A branch and cut algorithm for the VRP with satellite facilities,” Institute of Industrial Engineers (IIE). IIE Transactions, vol. 30, no. 9, pp. 821–834, 1998.

P. C. Pop, I. Zelina, V. Lupse, C. P. Sitar, and C. Chira, “Heuristic algorithms for solving the generalized vehicle routing problem,” International Journal of Computers, Communications & Control, vol. 6, no. 1, pp. 158–165, 2011.

V. F. Yu, S.-W. Lin, W. Lee, and C.-J. Ting, “A simulated annealing heuristic for the capacitated location routing problem,” Computers & Industrial Engineering, vol. 58, no. 2, pp. 288–299, 2010.

M. D. Morris, “Factorial sampling plans for preliminary computational experiments,” Technometrics, vol. 33, no. 2, pp. 161–174, 1991.