Dendrites with a closed set of end points

G Model

ARTICLE IN PRESS

ASOC 3307 1–11

Applied Soft Computing xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Applied Soft Computing journal homepage: www.elsevier.com/locate/asoc

The green vehicle routing problem: A heuristic based exact solution approach

1

2

3

Q1

4 5

C¸a˘grı Koc¸ a,∗ , Ismail Karaoglan b a b

CIRRELT and Canada Research Chair in Distribution Management, HEC Montréal, Montréal, Canada Department of Industrial Engineering, Selc¸uk University, Konya, Turkey

6

7 20

a r t i c l e

i n f o

a b s t r a c t

8 9 10 11 12 13

Article history: Received 22 March 2015 Received in revised form 4 August 2015 Accepted 22 October 2015 Available online xxx

14

19

Keywords: The green vehicle routing problem Mixed integer programming formulation Branch-and-cut algorithm Simulated annealing

21

1. Introduction

15 16 17 18

22Q2 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

This paper develops a simulated annealing heuristic based exact solution approach to solve the green vehicle routing problem (G-VRP) which extends the classical vehicle routing problem by considering a limited driving range of vehicles in conjunction with limited refueling infrastructure. The problem particularly arises for companies and agencies that employ a fleet of alternative energy powered vehicles on transportation systems for urban areas or for goods distribution. Exact algorithm is based on the branchand-cut algorithm which combines several valid inequalities derived from the literature to improve lower bounds and introduces a heuristic algorithm based on simulated annealing to obtain upper bounds. Solution approach is evaluated in terms of the number of test instances solved to optimality, bound quality and computation time to reach the best solution of the various test problems. Computational results show that 22 of 40 instances with 20 customers can be solved optimally within reasonable computation time. © 2015 Elsevier B.V. All rights reserved.

Today’s competitive economic environment requires strategic and operational decisions for companies in order to optimize and manage their logistic processes more efficiently. One of the most important operational decision concerns the design of vehicle routes since it offers great potential to reduce the costs and to improve the service quality. The classical vehicle routing problem (VRP) aims at routing a fleet of vehicles on a given network to serve a set of customers under specified supply and demand related constraints. Minimizing the total distance traveled by all vehicles or minimizing the overall travel cost are the typical objectives of the VRP and usually the cost is computed as a linear function of distance. Since its introduction by Dantzig and Ramser [1], the VRP and its variants have been studied extensively by researchers. Many heuristics have been developed in recent years for several variants of the VRP (see [2–5]). For a recent coverage of the state-of-the-art models and solution algorithms, the reader is referred to the survey by Cordeau et al. [6], Golden et al. [7] and Laporte [8], and to the books by Golden et al. [7] and Toth and Vigo [9].

∗ Corresponding author. Tel.: +1 514 343 6111; fax: +1 514 343 6111. E-mail addresses: [email protected] (C¸. Koc¸), [email protected] (I. Karaoglan).

The classical VRP assumes that the vehicle fuel tank capacity is unlimited and the fuel amount in the tank is always sufficient to serve all customers in any possible route. However, in reallife, vehicles need to refuel their tanks to continue and complete their tour. This situation is frequently encountered in the case of companies or agencies having alternative energy powered fleets (i.e., natural gas, electricity, ethanol) in which routes have to be planned taking additional difficulties associated with the limited refueling infrastructure into account. The alternative energy powered fleet operations (or more generally green logistics concept) have emerged as one of the latest extensions of the VRP literature in recent years. For example, many studies suggest that there are several opportunities for reducing carbon dioxide (CO2 ) emissions by extending the traditional VRP objectives to account for wider environmental and social impacts rather than just economic costs [10–12]. These studies are motivated by the activities of the transportation industry which has significant negative impacts on the environment, economy and human health. These impacts include increased resource consumption, toxic effects on ecosystems and humans, increased air and noise pollution, and also the climatic effects induced by greenhouse gas (GHG) emissions. GHG and in particular CO2 emissions are the most concerning ones as they have direct effects on human health, e.g., pollution, and indirect ones, e.g., climate change. Growing concerns about such hazardous effects of transportation on the environment call for revised planning approaches to road transportation by explicitly accounting for

http://dx.doi.org/10.1016/j.asoc.2015.10.064 1568-4946/© 2015 Elsevier B.V. All rights reserved.

Please cite this article in press as: C¸. Koc¸, I. Karaoglan, The green vehicle routing problem: A heuristic based exact solution approach, Appl. Soft Comput. J. (2015), http://dx.doi.org/10.1016/j.asoc.2015.10.064

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

G Model

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

ASOC 3307 1–11

ARTICLE IN PRESS

2

C¸. Koc¸, I. Karaoglan / Applied Soft Computing xxx (2015) xxx–xxx

such negative impacts. Juan et al. [13] studied the fleet size and mix vehicle routing problem with multiple driving ranges in which the total distance that each vehicle type can travel is limited. This problem arises in the routing of electric and hybrid-electric vehicles which can only travel limited distances due to the limited capacity of their batteries. A mathematical model is formulated and a multi-round heuristic is developed. The method is based on a biased randomized algorithm which can be used alone to create alternative fleet choices whenever the feasibility of the prespecified fleet configuration is not guaranteed. A set of benchmark instances were created to analyze how distance-based costs increase when considering “greener” fleet configurations. The method performed well on all benchmark instances and many different alternative solutions offer competitive distance-based costs while using fewer long- or medium-range vehicles than normally required. More recently, Koc¸ et al. [14] studied the fleet size and mix pollution-routing problem where the objective is a linear combination of vehicle, fixed cost, fuel cost and CO2 emissions, and driver cost. The authors formally defined the problem, presented a mathematical model and developed a hybrid evolutionary metaheuristic. For a further coverage of green issues at the operational level the reader is referred to the book chapter of Eglese and Bektas¸ [15] and to the surveys of Demir et al. [16] and Lin et al. [17]. Erdogan and Miller-Hooks [18] introduced the green vehicle routing problem (G-VRP). The problem design least-cost delivery routes from a depot to a set of geographically scattered customers within a pre-specified time limit and without exceeding the vehicle’s driving range that depends on fuel tank capacity to minimize the total distance traveled and/or total cost. In the G-VRP, vehicles have limited fuel tank capacity and are allowed to refuel when needed. The vehicles may be refueled at a limited number of fuel stations (FSs) which are available in the service area and at the depot node. In practice, the G-VRP is encountered, particularly, when the vehicle fleet includes the alternative fuel vehicles (AFVs). The authors proposed a mathematical formulation, and developed two construction heuristics; the Modified Clarke and Wright Savings heuristic and the density-based clustering algorithm, and a customized improvement technique. Another problem that is closely related to the G-VRP is the VRP with satellite facilities (VRP-SF) in which replenishment of a vehicle is allowed from another facility different from the depot. Bard et al. [19] have formulated the VRP-SF as a mixed integer programming (MIP) formulation with capacity and tour duration limitation constraints. Vehicles with capacity limitations have the option to stop at satellite facilities to reload in order to serve customer demand at the nodes. In their formulation, dummy nodes are introduced to represent multiple stops at intermediate depots for reloading vehicles with goods for delivery process. The authors have also developed an exact solution procedure (B&C algorithm) for solving VRP-SF to optimality. Regarding the complexity of the problem, different heuristic/meta-heuristic approaches have been also proposed to solve larger VRP-SF instances [20–22]. Frade et al. [23] studied the location of electric-vehicle charging stations in the city of Lisbon. The authors used a maximal covering model, and the aim to minimize the level of service and the number of charging stations. Chen et al. [24] later studied the electric vehicle charging stations location problem where a parking-based assignment method is presented for the city of Seattle. Nie and Ghamami [25] developed a conceptual optimization model to investigate travel of electric vehicles along a long corridor. Cavadas et al. [26] later described a method to locate the electric charging stations in the city of Coimbra where the aim is to maximize the number of electric vehicles under a fixed budget for building the stations. More recently, Schneider et al. [27] introduced an extended version of the G-VRP by considering an electric vehicle fleet with time windows, recharging at stations and limited vehicle load

capacity. The authors formulated the electric vehicle routing problem with time windows and recharging stations (E-VRPTW) as a MIP and employed a hybrid heuristic solution procedure that combines a variable neighborhood search algorithm with a Tabu Search heuristic. The authors have also improved the MIP formulation of Erdogan and Miller-Hooks [18] and used it to solve the set of small G-VRP benchmark instances with CPLEX. The algorithm is tested on a benchmark instances derived from the literature and on a newly generated benchmark instances. The results have shown that the proposed heuristic performed well and the hybridization mechanism had positive impacts on the solution quality. Schneider et al. [28] later studied the vehicle routing problems with intermediate stops (VRP-IS) in which stopping requirements at intermediate facilities may include replenishment/disposal and refueling/recharging stops. The authors have developed an adaptive variable neighborhood search algorithm to solve the VRP-IS instances. The algorithm is tested on several related VRP instances (G-VRP, E-VRPTW, etc.) in the literature and on the new benchmark instances. The authors reported that the proposed algorithm shows a satisfactory performance compared to the methods from the literature and is able to obtain numerous new best solutions. Felipe et al. [29] proposed constructive and local search heuristics within a simulated annealing framework to solve a variant of the G-VRP which considers multiple technologies and partial recharges. The authors tested their solution method on a newly generated benchmark instances, which indicate the efficiency of the algorithm. This brief review shows that only several heuristic algorithms have already solved the G-VRP, which does not guarantee the optimality. We believe there exists merit for the development of a new solution approach based on simulated annealing heuristic and branch-and-cut (B&C) algorithm, which is capable of optimally solving the G-VRP. This is the main motivation of this paper. The contributions of this paper are as follows: we introduce an efficient and powerful new solution approach to solve the G-VRP. We develop a new mathematical formulation having fewer variable and constraints without network augmentation, and adapts a set of valid inequalities to strengthen the linear programming relaxation of the formulation. We propose a simulated annealing heuristic to improve initial solution and upper bounds found during the search process of the solution approach. The rest of this paper is organized as follows. The problem description and a brief review of the literature are given in Section Q3 2. Section 2.1 shows the formulation of Erdogan and Miller-Hooks [18]. Section 2.2 presents the proposed mixed integer programming formulation, and subsequently Section 3 introduces valid inequalities. The proposed solution approach is described in Section 4. Finally, Section 5 presents the result of computational experiment results, followed by conclusions in Section 6.

2. Problem definition and formulations The G-VRP is defined on a complete directed graph G = (N, A), where N = “0” ∪ Nc ∪ F is the set of nodes, “0” corresponds to the depot, Nc = {n1 , n2 , . . ., nc } is the customer nodes, F = {nc+1 , nc+2 , . . ., nc+s } is the fuel stations, and A = {(ni , nj ):ni , nj , ∈ N} is the set of arcs that connect nodes in N. An unlimited number of homogeneous vehicle fleet is available at the depot to serve customers with fuel tank capacity Q (gallons) and fuel consumption rate r (gallons per mile). Each vehicle travels on the graph with constant speed sp (miles per hour). Each arc (ni , nj ) ∈ A is associated with a nonnegative distance dij , travel time tij associated with distance (tij = dij /sp) and also the triangular inequality holds (dij + djk ≥ dik ). Note that our formulation is also valid when the triangular inequality does not hold. Each node is associated with a service time pi , which represents the service

Please cite this article in press as: C¸. Koc¸, I. Karaoglan, The green vehicle routing problem: A heuristic based exact solution approach, Appl. Soft Comput. J. (2015), http://dx.doi.org/10.1016/j.asoc.2015.10.064

133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180

181

182 183 184 185 186 187 188 189 190 191 192 193 194 195

G Model

ARTICLE IN PRESS

ASOC 3307 1–11

C¸. Koc¸, I. Karaoglan / Applied Soft Computing xxx (2015) xxx–xxx 196 197 198 199 200 201 202 203

208 209 210 211

231

In this section, we present the formulation for the G-VRP in the literature. This formulation is originally developed by Erdogan and Miller-Hooks [18] (abbreviated as FEMH ) and improved by Schneider et al. [27]. In FEMH , graph G is augmented as G = (N , A ) with a set of s  dummy nodes, F = nc+s+1 , nc+s+2 , . . ., nc+s+s , one for each potential visit to a fuel station (FS) or depot serving as a FS (i.e., N  =  N ∪ F and A = {(ni , nj ) : ni , nj ∈ N  }). Based on this augmentation, the decision variables and FEMH are presented as follows: Decision variables: xij Binary variable equal to 1 if a vehicle travels from node i to j, 0 otherwise ∀i, j ∈ N . fi Fuel level variable specifying the remaining fuel level upon arrival to node i and it is reset to Q at each refueling station node i and at the depot (∀ i ∈ N ).  i Time variable specifying the arrival time of a vehicle at node i, initialized to zero upon departure from the depot (∀ j ∈ N ). Model (FEMH ):

232

min

214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230



dij xij

(1)

(i,j) ∈ N  ,i = / j 233

subject to

234

s.t.



j ∈ N  ,i

 235

j ∈ N  ,i 236

i ∈ N  ,i

i ∈ N  ,i

= / j

N  \ 0”

N  \Nc



(3)

xij = 0 ∀j ∈ N 

(4)

= / j

(5)

) xi0 ≤ m

238

(6)

N  \ 0”

)







j ≥ i + tij + pi xij − T max 1 − xij



241



(2)

xoj ≤ m



240



xji −

237

i∈(

= / j

xij ≤ 1 ∀i ∈



j∈(

xij = 1 ∀i ∈ Nc

= / j



239

 

∀i ∈ N  ; ∀j ∈ N  \ 0

; i= / j

 (7)

242 243



Objective function (1) minimizes the total distance traveled by the vehicles. Constraint sets (2)–(4) are known as degree and flow conservation constraints. While constraint sets (2) and (3) ensure that each customer must be visited exactly once, and each FS (and associated dummy nodes) will have at most one successor node, constraint set (4) guarantees that entering and leaving arcs to each node are equal. Constraint sets (5) and (6) ensure that at most m vehicles are used to serve customers. Value m should be set as the number of available vehicles in case of the number of vehicle restrictions, otherwise these constraints should be removed. Constraint sets (7)–(9) both determine the arrival time at each node and ensure that each vehicle returns to the depot no later than Tmax . Constraint set (10) is associated with a vehicle’s fuel level based on node sequence and type. Constraint sets (7) and (10) are also eliminates the illegal routes (such as cycle without depot node). Constraint set (11) replenishes the vehicle fuel tank level to Q when it arrives at any FS. Constraint set (12) guarantees that there will be enough remaining fuel to return to the depot by the way of a FS from any customer location en route. Finally, constraint set (13) is known as integrality constraints.

2.1. The formulation of Erdogan and Miller-Hooks [18]

207



∀j ∈ N  \ 0

• • • •

213

206



t0j ≤ j ≤ T max − tj0 + pj

212

205



time at the customer node and refueling time at fuel station node. The depot can serve as a fuel station and all fuel stations have unlimited capacities. There is no limitation on the number of stops for refueling and fuel tank of vehicle is assumed to be full after leaving a fuel station (FS). There exists a time horizon (Tmax ) which establishes the duration of a workday. The problem is to determine the corresponding vehicle routes so as to minimize the total cost subject to the following assumptions: Each vehicle is used for at most one route. Each route starts and ends at the depot. Each customer is served by exactly one vehicle. Fuel level at the vehicle’s tank must be greater than or equal to the fuel consumption between any two nodes. • The amount of fuel in a vehicle’s tank is sufficient to be able to visit between any pair of nodes. • The duration of the route assigned to a vehicle cannot exceed Tmax .

204

3

0 ≤ 0 ≤ T max



fj ≤ fi − r · dij xij + Q 1 − xij fj = Q

∀j ∈ N  \Nc

fj ≥ r · dji xji xij ∈







∀j ∈ Nc ; ∀i ∈ N  ; i = / j



∀i ∈ (F ∪ F) ; ∀j ∈ Nc

0, 1

∀ (i, j) ∈ N



(9)

244

(10)

245

(11)

246

(12)

247

(13)

248

2.2. A new formulation for the G-VRP

dijk = dik + dkj

∀i, j ∈ N; ∀k ∈ F; i = / j

249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268

269

The G-VRP formulation described in previous section requires the augmentation of a network with dummy nodes in order to represent multiple visits to an alternative fuel station. In the worst case, the number of dummy nodes for each FS visit can be equal to the number of customers (i.e., |N|) which increase the total number of nodes in the augmented network from (|Nc | + |F|+1) to (|Nc | + |Nc ||F|+1). Although, it is acceptable for small test instances, this method requires quite more solution time for larger instances since the number of binary variables in a formulation increases with the additional dummy nodes. As it is known, the size of enumeration tree in any MIP solver increases exponentially based on the number of binary variables in a problem. Furthermore, the memory requirement of any commercial software increases with the size of the matrix dimensioned with the number of variables and constraints. Thus, it is expected that the formulation obtained by augmentation of a network has poorer performance than that of one with no augmentation in terms of linear programming (LP) relaxation, CPU time, etc. In this paper, the G-VRP is formulated without network augmentation. In our formulation, Miller–Tucker–Zemlin (MTZ) capacity and subtour elimination constraints for the VRP are adapted for the GVRP, which were proposed by Miller et al. [30] for the traveling salesman problem (TSP). Kulkarni and Bhave [31] adapted them to the capacitated VRP. Desrochers and Laporte [32] lifted them and introduced new bounding constraints on additional variables. Kara et al. [33] corrected the lifted version of these constraints, and finally, exact meanings to auxiliary variables were given by Kara [34]. For the sake of simplicity, we add the service time (pi ) to the succeeding arc (i . e . , tij = pi + (dij /sp)) and introduced two new  parameters which are indirect distance dijk and indirect travel time  t ijk calculated as follows: 

(8)





(14)

Please cite this article in press as: C¸. Koc¸, I. Karaoglan, The green vehicle routing problem: A heuristic based exact solution approach, Appl. Soft Comput. J. (2015), http://dx.doi.org/10.1016/j.asoc.2015.10.064

270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301

302

G Model ASOC 3307 1–11

ARTICLE IN PRESS

4

C¸. Koc¸, I. Karaoglan / Applied Soft Computing xxx (2015) xxx–xxx





303



∀i, j ∈ N; ∀k ∈ F; i = / j

t ijk = pi + pk + dijk /sp

(15)



304 305 306 307 308 309 310 311 312 313 314 315 316 317 318

319

where dijk defines the distance between node i and j through fuel  station k and t ijk determines the total duration which includes the service time, refueling time and en route time on arcs (i, k) and (k, j). Based on these definitions, the decision variables of the G-VRP are presented as follows: Decision variables: xij Binary variable equal to 1 if a vehicle travels from node i to j directly, 0 otherwise (∀ i, j ∈ (0 ∪ Nc )). yijk Binary variable equal to 1 if a vehicle travels from node i to j through fuel station k, 0 otherwise (∀ i, j ∈ (0 ∪ Nc ) ; ∀ k ∈ F). fi Fuel level variable specifying the remaining tank fuel level upon arrival to node i (∀ i ∈ Nc ).  i Time variable specifying the arrival time of a vehicle at node i (∀ i ∈ Nc ). Model (FKK ):



min

dij xij +

i,j ∈ (0∪Nc )

320

321

323

xij − xji



+





(17)



= 0 ∀i ∈ (0 ∪ Nc )



3. Valid inequalities for the G-VRP











≤ M1ijk − tij − t ijk



j ≥ t0j x0j +



∀i, j ∈ Nc ; ∀k (0 ∪ F)

(19)



(t 0jk y0jk ) ∀j ∈ Nc



j ≤ Tmax − Tmax − t0j x0j −

i ≤ Tmax − ti0 xi0 −





(20)



Tmax − tˆ0jk y0jk

∀j ∈ Nc

(21)

k∈F

tˆi0k yi0k



(22) j ∈ NC

k∈F 330 331 332



333 334



f j − fi + M2ij xij + M2ij − r · dij − r · dji xji ≤ M2ij − r · dij

fj ≤ Q − r · d0j x0j −

∀i, j ∈ Nc ; i = / j

 

(23)

r · dkj yijk



∀j ∈ Nc

(24)

i ∈ (0∪Nc ) k ∈ F

335

fi ≥ r · di0 xi0 +

 

r · dik yijk



∀i ∈ Nc

(25)

j ∈ (0∪Nc ) k ∈ F 336

i , fi ≥ 0 ∀i ∈ Nc

337

xij ∈

338

yijk ∈





0, 1





0, 1

(26)

∀i, j ∈ (0 ∪ Nc )

(27)

∀i, j ∈ (0 ∪ Nc ) ; k ∈ F.

(28)

340

In the given formulation, two constants (M1ij and M2ij ) are used for all (i, j) pairs and are evaluated as follows:

341

M1ijk = T max + tij + tˆijk − ti0 − t0j

339



342

M2ij = Q + r · dij − min r · dkj k∈F



∀i, j ∈ Nc ; k ∈ F



− min r · dik k∈F



∀i, j ∈ Nc

344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366

(29) (30)

367

The success of any exact solution procedure mainly depends on the performance of the upper and lower bounding procedures. Rapid improvements on these bounds shorten the solution procedure. In this section, we describe the valid inequalities derived from the VRP literature to improve the lower bounds. Valid inequalities are one of the most practical ways to strengthen the linear relaxations of the formulations. These inequalities eliminate some fractional solutions from the solution space in such a way that a stronger lower bound can be obtained for the problem. The first valid inequality which bounds below the number of routes originating from the depot is given as follows:



∀i ∈ Nc

343

(18)

+ M1ijk − tij yijk + M1ijk − tij − t ijk − t jik yjik

k∈F

329

yijk − yjik





328

 



326 327

yijk ⎠ = 1 ∀j ∈ Nc

i − j + M1ijk − t ijk xij + M1ijk − t ijk − tij − tji xji



325

(16)

k ∈ (0∪F)



324

dijk yijk

k ∈ (0∪F)

i ∈ (0∪Nc ) 322



⎞



⎝xij +

i ∈ (0∪Nc )

 



i,j ∈ (0∪Nc )k ∈ (0∪F)

⎛



s.t.



Objective function (16) minimizes the total distance traveled by the vehicles. Constraint sets (17) and (18) are known as degree constraints. While constraint set (17) ensures that each customer must be visited exactly once, constraint set (18) guarantees that entering and leaving arcs to each node are equal. The arrival time of a vehicle to each customer is determined by constraint set (19). Constraint set (19) along with constraint sets (20)–(22) guarantee that each vehicle must return to the depot before time Tmax . Constraint sets (20) and (21) both bound the decision variables ( i ) and specify the arrival time of a vehicle at the first customer in any route. Constraint set (22) ensures that the arrival time of a vehicle at the last customer must allow the vehicle to return to the depot before time Tmax . Constraint set (23) determines the fuel level at each customer node based on the distance and the vehicle’s fuel consumption rate. Constraint set (24) specifies the vehicle’s fuel level at the first customer node of any route or first customer node after refueling. Constraint set (25) guarantees that vehicle’s fuel level is enough to return to the depot or to reach a FS. Constraint sets (26)–(28) are known as non-negativity and integrality constraints, respectively. It should be noted that whereas FKK has O(|NC |2 |F|) binary variables, O(|NC |) continuous variables and O(|NC |2 |F|) constraints, FEMH has O(|NC |2 |F|2 ) binary variables, O(|NC | |F|) continuous variables and O(|NC |2 |F|2 ) constraints.

 x0j +



368 369 370 371 372 373 374 375 376 377 378

y0jk

≥ rG-VRP (NC )

(31)

379

k∈F

where rG-VRP (NC ) is the minimum number of vehicles. Similar bounding constraint has been used by Achuthan et al. [35] for the capacitated VRP and [36,37] for the location routing problem with simultaneous pickup and delivery. Although, the minimum number of vehicles required can easily be evaluated for vehicle routing problems ⎡ having ⎤customer load and vehicle capacity (i.e., rVRP (NC ) = ⎢ ⎢



⎢i ∈ NC

pi /Cap⎥ ⎥: where pi is

⎥

the demand of customer i, and Cap is the vehicle capacity), it is not as easy as for routing problems without load restrictions. In this paper, we use the following property of any feasible route for G-VRP to obtain rG-VRP (NC ).

380 381 382 383 384 385

386

387 388 389

390 391

Property 1 us consider a feasible route P=  Let v0 , v1 , v2 , . . ., vs−1 , vs , vs+1 with arbitrary s customers (as illustrated in Fig. 1) where v0 and vs+1 represent the depot node "0". The total en route time of P, called as tP , is evaluated as tP = tv0 ,v1 + tv1 ,v2 + · · · + tvs−1 ,vs + tvs ,vs+1

(32)

Please cite this article in press as: C¸. Koc¸, I. Karaoglan, The green vehicle routing problem: A heuristic based exact solution approach, Appl. Soft Comput. J. (2015), http://dx.doi.org/10.1016/j.asoc.2015.10.064

392 393 394 395 396

397

G Model

ARTICLE IN PRESS

ASOC 3307 1–11

C¸. Koc¸, I. Karaoglan / Applied Soft Computing xxx (2015) xxx–xxx

5

knowledge, this evaluation scheme is the first for the VRP without customer demands. Constraint (31) is a special case of the following exponentialsize constraints derived from capacity and subtour elimination constraints of the VRP [38].







xij +



 

≤ S  − rG-VRP (S)

yijk

 

∀S ⊂ NC , S  > 2

(35)

where rG-VRP (S) is calculated as in constraint (31). The constraint set (35) guarantees that the number of vehicles visiting a set of customers is not smaller than the corresponding lower bound. Another exponential-size inequality for the G-VRP is based on the generalized large multistar (GLM) inequalities which have been originally proposed for the VRP [39,40]. We have adapted the GLM for the G-VRP as follows:



 

xij +

i ∈ S j ∈ (N\S)

Fig. 1. The total en route time evaluation for any feasible solution.

399

or equivalently tP =

+··· +

400

401

402

tv0 ,v1 + 2

t

+

2

t

vs−1 ,vs + tvs ,vs+1



2

+

tv0 ,v1 + 2

s  

+

tvs ,vs+1

(33)

2

404 405 406 407 408 409

tvi−1 ,vi + tvi ,vi+1 2

+

tvs ,vs+1 2

Table 1 Procedure to obtain the minimum number of vehicles for a set of nodes. Procedure: Evaluation of rG-VRP (P) Input: Node set (P) not including the depot node Output: rG-VRP (P) Step 1. Initialize parameters t¯ P ← 0, NumDepLink ← 0,  ← ∅, //Route time evaluation without depot node Step 2. Perform the following steps to all nodes  (v) in P Step 2.1. Find the first (¯vcur ) and second vcur

nearest node to (vcur ).

v¯ cur ← argmini ∈ P (ti,vcur ) vcur ← argmini ∈ P:i =/ v¯ cur (tvcur ,i ) Step 2.2. Evaluate the effect of vcur to t¯P tv¯ cur ,vcur +t vcur ,vcur 2

t¯ P ← t¯P + //Route time evaluation considering depot node Step 3. Find the best unlinked customer node (i*) in P causing least route time increment because of linking to depot



∗

t

t0,i −

i,vi



(34)



t

t0i∗ −

i∗ ,vi∗

414

415

416 417 418 419 420 421 422

423

⎛

 

 − → tx ij xij + ← txji xji

ty ijk yijk + ← tyjik yjik



424

∀S ⊂ NC , S = / ∅

(36)

425

− → tx ij =

ti,0 t

i,j ← txji =

 − → ty ijk = ← tyjik =

426

if node j is depot + tj,0

otherwise

t0,i

if node j is depot

t0,j + tj,i

otherwise

tˆi0

if node j is depot

ˆ tijk + tj,0 otherwise tˆ0i

if node j is depot

t0,j + tˆjik

otherwise

∀i, j ∈ Nc

∀i, j ∈ Nc ; k ∈ F

(37)

427

(38)

428

The constraints set (36) ensure that the number of leaving arcs (or equivalently the number of routes) from a set of customers (S) must satisfy the required number of routes. These inequalities are useful when t¯S /Tmax is close to the next integer. 4. Solution approach This section presents a detailed description of our heuristic based exact solution approach for the G-VRP. The method implements a combination of cutting planes and implicit enumeration to solve any combinatorial optimization problem. The basic idea is the identification of the violated inequalities that are valid throughout the enumeration tree. Therefore, at each step of the solution method, violated inequalities which are identified by solving the separation problem are added to the formulation and the corresponding linear program is reoptimized.

429 430 431 432

433

434 435 436 437 438 439 440 441 442

2

Step 4. Evaluate the effect of k to t¯P t¯ P ← t¯P +

yijk

i ∈ S j ∈ (N\S)



In the best case, the preceding (vi−1 ) and the succeeding node (vi+1 ) of each customer (vi ) on any route are the first and the second nearest nodes to vi in P, respectively. Therefore, the total en route time evaluation without sequence information in P using the first and the second nearest nodes, called t¯P , gives a lower bound on tP .䊐 According to Property 1, a polynomial time procedure is used to obtain rG-VRP (NC ) which is illustrated in Table 1. To the best of our

i ← argmini ∈ P:i/∈



where t¯S is given as in Eq. (34) for the customer set S and



413

k∈F

i=1 403



 − →

2

412

k∈F

⎝t¯S +

Tmax

 t  v0 ,v1 + tv1 ,v2 v1 ,v2 + tv2 ,v3

and more formally tP =

1

≥

411

k∈F

(i,j) ∈ S

398

410



4.1. Simulated annealing

443

2

NumDepLink ← NumDepLink + 1 r G−VRP (P) ← t¯P /Tmax ←  ∪ i * Step 5. if NumDepLink < (2 * rG-VRP (P)) then go to Step 3 else report rG-VRP (P) and stop

At each iteration of the simulated annealing (SA) heuristic, neighbors of the current solution are generated by using all of the moving strategies. The best one among them is chosen as a new solution (Snew ) for the problem. In recent years, SA has received increased attention for solving several routing problems [41]. Banos

Please cite this article in press as: C¸. Koc¸, I. Karaoglan, The green vehicle routing problem: A heuristic based exact solution approach, Appl. Soft Comput. J. (2015), http://dx.doi.org/10.1016/j.asoc.2015.10.064

444 445 446 447 448

G Model ASOC 3307 1–11

ARTICLE IN PRESS

6

C¸. Koc¸, I. Karaoglan / Applied Soft Computing xxx (2015) xxx–xxx

483

et al. [42] developed a SA heuristic which is based on a multiobjective procedure to solve the VRP with Time Windows. The SA algorithm implements a candidate list strategy in generating the neighbors since searching whole neighborhood of the current solution by a moving strategy is a very time consuming process. According to this strategy, each moving step randomly generates a subset of neighbors, LS, in order to satisfy the G-VRP constraints. Then, this neighbors are gathered in a pool to select the new solution (Snew ). If the new solution is better than the current solution then it is accepted as the current  solution.  Otherwise, it is accepted with probability of exp −s/Titr as the current solution. s is the relative percent deviation of the quality of the new solution from the current solution which is calculated by [(f(Snew ) − f(Scur ))/f(Scur )] * 100. In each iteration of the SA, temperature (Titr ) is reduced using a geometric cooling schedule, i.e., Titr ← ˛Titr−1 where ˛ is the cooling rate (0 < ˛ < 1). At the beginning of the SA, the temperature is set to initial temperature (i.e., Titr ← T0 ) and the SA stops when the temperature reaches the final temperature (Tf ). SA heuristics use a set of correlated parameters and configuration decisions. In our implementation, we initially used the parameters suggested by [36,37] for the SA, where an extensive meta-calibration procedure was applied to generate effective parameter values. But we have conducted several experiments to further fine-tune these parameters. After several preliminary experiments, following parameter values are implemented in the SA. The initial temperature (T0 ) is set 665 in which an inferior solution (inferior by 70% relative to current solution) is accepted with a probability of 0.90, the final temperature (Tf ) is chosen 0.15 such that a solution which is inferior by 1% relative to current solution is accepted with a probability of 0.001. The number of neighbors generated by each moving strategy is set |NC | for Swap operator, |F| for Merge, FSAdd and FSDrop operators. Furthermore, cooling rate is selected as 0.95. Because of the stochastic nature of the SA, it is run 5 times for each instance with different random seeds.

484

4.2. B&C algorithm

449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482

Table 2 The branch-and-cut algorithm for the G-VRP. Algorithm: Branch-and-cut for the G-VRP Input: G-VRP problem data Output: Optimal solution (Sbest ) for the G-VRP Step 1.

Step 2.

Step 3. Step 4.

Step 5. Step 6.

Step 7.

Step 8.

Step 9. 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512

The steps of the proposed B&C algorithm for the G-VRP are presented in Table 2, where LP is the LP-relaxation of the MIP model, St∗ is the optimal solution of LP at a particular node t, S0 is the initial feasible solution, ˚ is the set of unexplored nodes of enumeration tree. The B&C is initialized by Step 1 in which the initial LP model is constructed, an initial feasible solution S0 is obtained, and S0 is accepted as the best feasible solution. The procedure to obtain an initial solution in our algorithm can be summarized as follows: first, we implement an extension of the well-known Clarke and Wright heuristic [43]. In this heuristic, a simple route for each customer from the depot is constructed. Then, all pairs of routes (R1 , R2 ) in the solution are examined in terms of cost savings obtained by their merge and the merge is implemented by combining the pair of nodes providing the largest cost savings. This strategy is repeated until no feasible combination is found. After an initial feasible solution (S0 ) is generated, a classical simulated annealing (SA) algorithm is implemented to obtain an initial solution (S0 ) for the B&C algorithm. SA, which stems from the simulation of the annealing of solids, is a stochastic search technique that is able to escape from local optima using a probability function [44]. It starts with an initial feasible solution (S0 ) and improves it until a stopping criterion is met. We employ the following moving strategies, two of which are well-known in the VRP literature and the others are developed for the G-VRP. Merge: Two routes (R1 and R2 ) are selected randomly, R2 is located to each best position in R1 , and the best feasible position is selected. An example of this operator is illustrated in Fig. 2. In this

Step 10.

Construct initial linear relaxation of the MIP model (LP). Obtain the initial solution (S0 ) Evaluate objective function value of initial solution (f(S0 )) Set Sbest ← S0 , f(Sbest ) ← f(S0 ) Add the root node to enumeration tree () and set the root node as current node (t) if (=∅) or (termination criterion is met) then report Sbest and stop else select a node (t) from enumeration tree (t ∈ ) end if Solve the corresponding LP and obtain the optimal solution (St∗ ). if St∗ is infeasible then ← \ t, go to Step 2 end if if f (St∗ ) > f (Sbest ) then ← \ t and go to Step 2 if St∗ is integer then if f (St∗ ) < f (Sbest ) then Sbest ← St∗ , f (Sbest ) ← f (St∗ ) end if ← \ t, go to Step 2 end if if (St∗ is not integer) then Obtain  feas  feasible solution from current fractional solution St if f (Stfeas ) < f (Sbest ) then Sbest ← Stfeas , f (Sbest ) ← f (Stfeas ) end if end if if a tailing off exists then go to Step 10 else solve separation problem to identify violated inequalities and put them cut pool (K) end if if k = / 0 then add violated inequalities to the LP, go to Step 3 end if Select a non-integer decision variable according to branching rule, create two new nodes and add them to ˚ and go to Step 2

figure, routes (0-6-2-5-7-1-0) and (0-3-4-0) are merged into new one as (0-6-2-5-7-3-4-1-0). Swap: Two randomly selected customers, which are in the same route or in two different routes, are exchanged. Swap operator is illustrated in Fig. 3. In this figure, customers 3 and 6 are swapped as illustrated in Fig. 5(a) and (b). FSAdd: Two successive customers are selected randomly and a FS visit is added between them if there is no previous FS visit. Although this operator may increase the routing cost, better results can be obtained at the successive iterations by using other operators

Fig. 2. A sample for Merge operator.

Please cite this article in press as: C¸. Koc¸, I. Karaoglan, The green vehicle routing problem: A heuristic based exact solution approach, Appl. Soft Comput. J. (2015), http://dx.doi.org/10.1016/j.asoc.2015.10.064

513 514 515 516 517 518 519 520 521 522

G Model

ARTICLE IN PRESS

ASOC 3307 1–11

C¸. Koc¸, I. Karaoglan / Applied Soft Computing xxx (2015) xxx–xxx

7

Fig. 5. A sample for FSDrop operator. Fig. 3. A sample for Swap operator.

∗ ) to the last customer i* in the route R continues largest xij∗ (or yijk as long as addition is feasible. This step is repeated until all cus∗ is selected as the largest value, tomers are assigned to routes. If yijk then a FS (k*) visit is included between the customers i* and j* (i.e., 0, . . ., i*, k*, j*, . . ., 0). Note that a simple route is constructed for customer j who is not assigned to any routes since nodes i with ∗ , y∗ > 0 are assigned previously to another route. After xij∗ , xji∗ , yijk jik

generating a feasible solution (St ), an improved solution (Stfeas ) is obtained by implementing SA algorithm explained in Step 1. To reduce the computation time spent by this procedure, SA is applied to all nodes up to depth level ten and all nodes on every tenth depth level in enumeration tree. Step 8 determines the valid inequalities, if a tailing off does not exist. In our B&C algorithm, inequality (31) is imposed directly to the model at the beginning of the algorithm (Step 1) since there is only one constraint in (31). However, we need a procedure to obtain constraints from the sets (35) and (36). For this purpose, we introduce a greedy constructive heuristic procedure to separate these inequalities. In each iteration of this procedure, a customer is selected randomly as a seed node and customer set S is initialized with this seed. Then, a new customer node (m*) that minimizing 35 } for (35) and the slack of the constraint (i.e. m∗ ← argmin {slackm

Fig. 4. A sample for FSAdd operator.

523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539

(i.e., merge or swap) as it is not possible for the current iteration due to the fuel tank capacity limitations. An illustrative example for this operator can be seen in Fig. 4. In this figure, a FS visit between customers 3 and 1 is added. After this insertion it may be possible to merge routes (0-2-3-6-1-0) and (0-4-5-0) into new one as (0-23-6-1-4-5-0) for the next iteration. FSDrop: A randomly selected FS visit is removed if possible. This operator is illustrated in Fig. 5. In this figure, FS visit between customers 2 and 3 are removed and a new arc is added to these customers. Steps 2–6 are the classical steps of the B&C algorithm. Step 2 specifies the termination criterion and selects an unexplored node having the smallest objective function value from enumeration tree for additional processing. Then, the next step solves the corresponding LP model and obtains the (probably fractional) optimal solution (St∗ ). Steps 4 and 5 prune the current node t, if St∗ is infeasible or its objective function value is

36 } for (36)) is selected and S is expanded with m∗ ← argmin {slackm

36 slackm =

 xij∗ +

i ∈ (S∪m)j ∈ (N\(S∪m))

⎛

⎜ ⎜ ⎜ − 1/Tmax ⎜t¯(S∪m) + ⎜ ⎝



m* (i.e. S ← S ∪ m*). The definitions of the slack values are given as follows for (35) and (36), respectively.





35 slackm = S ∪ m − rG-VRP (S ∪ m) −





xij∗ +





∗ yijk

(39)

557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574

575

576 577

578

579

∗ yijk

580

k∈F

i ∈ (S ∪ m)

556

k∈F





555

m ∈ (NC \S)

(i,j) ∈ (S∪m)



554

m ∈ (NC \S)

540



553

⎞  − → ∗ tx ij xij + ← txji xji∗ +

 − →

∗ ty ijk yijk

k∈F

+←

∗ tyjik yjik

⎟  ⎟ ⎟ ⎟ ⎟ ⎠

(40)

581

j ∈ (N\ (S ∪ m))

541 542 543 544 545 546 547 548 549 550 551 552

worse than that of the best feasible solution found so far (Sbest ), respectively. Step 6 prunes the current node if St∗ is integer, and then Sbest is updated if a better integer solution is found. In B&C literature, it has been observed that optimal fractional solutions may include some information about optimal integer solutions, and it is highly possible to obtain good quality integer solutions from fractional solutions by using a heuristic and/or exact solution procedure. Based on this idea, we obtain a feasible solution from St∗ by using the following heuristic procedure in Step 7. ∗ be an optimal fractional solution of any node on enuLet xij∗ and yijk meration tree. Firstly, a route R is started by selecting customer i* ∗ (or y∗ ) and then appending a new customer j* with with largest x0i 0ik

where rG-VRP (S ∪ m) and t¯(S∪m) is defined as previously. When the customer m* is selected, the validity of the current constraint is checked. If the violation is occurred for the set S is added to the cut pool. This procedure is repeated until all violated cuts have been identified. These cuts are then imposed to the LP model and reoptimized in Step 9. Finally, Step 10 creates two new nodes by applying branching rule, if a tailing off exists or new violated inequalities are not identified. 5. Computational results This section presents the computational results of numerical experiments. Firstly, we give brief information about the test

Please cite this article in press as: C¸. Koc¸, I. Karaoglan, The green vehicle routing problem: A heuristic based exact solution approach, Appl. Soft Comput. J. (2015), http://dx.doi.org/10.1016/j.asoc.2015.10.064

582 583 584 585 586 587 588 589

590

591 592

G Model ASOC 3307 1–11

ARTICLE IN PRESS

8

C¸. Koc¸, I. Karaoglan / Applied Soft Computing xxx (2015) xxx–xxx

Table 3 Linear programming relaxation of formulations on the test instances. Scenario 1

Scenario 2

Scenario 3

EMH

KK

Instance

EMH

KK

Instance

EMH

KK

Instance

EMH

KK

20c3sU1 20c3sU2 20c3sU3 20c3sU4 20c3sU5 20c3sU6 20c3sU7 20c3sU8 20c3sU9 20c3sU10

61.38 56.48 57.49 51.00 53.90 51.91 61.87 52.86 66.51 41.23

38.05 26.22 40.33 29.30 30.70 31.37 35.54 35.77 40.59 13.90

20c3sC1 20c3sC2 20c3sC3 20c3sC4 20c3sC5 20c3sC6 20c3sC7 20c3sC8 20c3sC9 20c3sC10

62.53 67.32 60.39 61.14 81.71 75.85 88.87 76.45 70.75 76.62

32.60 52.08 31.31 32.44 48.67 60.89 47.10 51.70 60.24 66.98

S1 S1 S1 S1 S1 S2 S2 S2 S2 S2

61.66 54.95 57.40 55.24 41.63 62.94 69.45 70.21 57.95 56.79

26.37 25.02 28.81 35.23 13.32 49.48 52.06 41.70 38.44 29.98

S1 S1 S1 S1 S1 S2 S2 S2 S2 S2

55.01 51.36 54.95 56.69 56.85 59.47 57.07 68.19 68.60 68.60

33.78 28.24 25.02 25.02 24.95 32.44 35.93 34.90 34.90 34.90

Average

55.46

32.17

72.16

48.40

58.82

34.04

59.68

31.01

598

problems. Then, we investigate the effectiveness of our formulation (FKK ) and Erdogan and Miller-Hooks’ formulation (FEMH ) with respect to linear programming relaxations. Finally, the performance of the proposed B&C algorithm is analyzed. All experiments were conducted on Intel Xeon 3.16 GHz equipped with 8 GB RAM computer (the operating system is Windows 7-x64).

599

5.1. Test problems

593 594 595 596 597

620

The algorithm is tested using the benchmark set generated by Erdogan and Miller-Hooks [18]. Each instance involves 20 customers and FSs varying from 2 to 10 (with increments of two). Each customer and FS node are located in a grid of 330 by 300 miles and the depot is located near the center of the grid. To investigate the effect of the location configuration of both customers and stations and the number of FSs, four different scenarios have been considered. In Scenario-1 (2) customers are located uniformly (clustered) and three FSs are fixed between the depot and the grid boundaries in westerly, northerly and southeasterly directions. In Scenario-3, while half of the test instances are selected from Scenario-1, others are selected from Scenario-2, and for each of them 3 additional FSs are located randomly. In Scenario-4, while half of the test instances are created from one instance of Scenario-1, others are created from one instance of Scenario-2, and the number of FSs is increased from 2 to 10 by increments of two for each instance. In the original test set, some instances are infeasible because of time and refueling restrictions, such as some customers cannot be served in a given maximum route duration or cannot be reached without more than one refueling stop. Therefore, these customers are identified and removed from the test instances.

621

5.2. Comparison of MIP formulations

600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619

622 623

Scenario 4

Instance

In this section, we compare two MIP formulations with respect to LP relaxation quality. The LP/MIP solver IBM ILOG CPLEX 12.5

2i6s 4i6s 6i6s 8i6s 10i6s 2i6s 4i6s 6i6s 8i6s 10i6s

4i2s 4i4s 4i6s 4i8s 4i10s 4i2s 4i4s 4i6s 4i8s 4i10s

is used as LP solver. To analyze the computational results, we use lower bound gap which is the percentage gap between the objective function value of the LP relaxation (ZLP ) corresponding to a particular formulation and the optimal or the best found integer feasible solution (Z*) obtained by proposed B&C algorithm. The percentage gap values are calculated as  = 100 * ((Z* − ZLP )/Z*). Note that, these gap values report how close the corresponding LP bounds to the optimal/best known solution. It is well-known that the quality of the bounds is one of the critical issues in reducing execution time of exact solution procedure. A tighter bound can help the exact solution procedure to get optimal solutions in a shorter computational time. Table 3 presents the computational results for each test problem. In Table 3, each triple of columns reports the name of the instance, the percentage gap values for the Erdogan and MillerHooks’ formulation (EMH ) and our formulation (KK ) for each scenario, respectively. As shown, our formulation provides better results than [18] formulation for all instances. Average percentage gap values on all scenarios are obtained as 36.41% for our formulation while this value is 61.53% for [18] formulation. These results indicate that it is plausible to use the proposed formulation in B&C algorithm since tight bounds are one of the most crucial components for exact solution procedures.

5.3. Results for the heuristic algorithm

Scenario 2

In this section, we investigate the effects of heuristic solution procedures (i.e., extended Clarke and Wright heuristic (ECW) and simulated annealing (SA) algorithm). Table 4 presents the computational results for each scenario, respectively. In this table, successive three columns give instance name, average percentage gap values for ECW and SA heuristics at the beginning of B&C, respectively. The gap values are evaluated as 100 * ((Zheur − UB)/Zheur ), where Zheur is the feasible solution value

Scenario 3

Scenario 4

Instance

GapECW

GapSA

Instance

GapECW

GapSA

Instance

GapECW

GapSA

Instance

GapECW

GapSA

20c3sU1 20c3sU2 20c3sU3 20c3sU4 20c3sU5 20c3sU6 20c3sU7 20c3sU8 20c3sU9 20c3sU10

17.17 6.56 15.56 15.28 14.76 5.85 2.51 22.74 0.56 10.85

0.87 0.00 0.30 7.95 0.00 1.90 0.98 3.06 0.00 10.85

20c3sC1 20c3sC2 20c3sC3 20c3sC4 20c3sC5 20c3sC6 20c3sC7 20c3sC8 20c3sC9 20c3sC10

14.81 2.67 23.05 4.76 9.04 10.61 43.86 9.49 3.26 9.97

12.01 0.16 2.98 4.76 1.61 7.68 43.86 0.00 0.00 2.19

S1 S1 S1 S1 S1 S2 S2 S2 S2 S2

6.34 22.27 6.76 23.77 11.59 1.18 46.59 25.62 0.00 17.85

0.65 13.79 3.56 4.91 11.59 0.75 46.59 0.20 0.00 4.61

S1 S1 S1 S1 S1 S2 S2 S2 S2 S2

7.98 17.01 22.27 22.27 22.38 4.76 3.49 22.87 22.87 22.80

0.00 1.61 13.79 0.03 13.34 4.76 1.45 2.26 7.11 7.05

Average

11.18

2.59

13.15

7.52

16.20

8.67

16.87

5.14

2i6s 4i6s 6i6s 8i6s 10i6s 2i6s 4i6s 6i6s 8i6s 10i6s

625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646

647

Table 4 Computational results of the heuristic solution procedures. Scenario 1

624

4i2s 4i4s 4i6s 4i8s 4i10s 4i2s 4i4s 4i6s 4i8s 4i10s

Please cite this article in press as: C¸. Koc¸, I. Karaoglan, The green vehicle routing problem: A heuristic based exact solution approach, Appl. Soft Comput. J. (2015), http://dx.doi.org/10.1016/j.asoc.2015.10.064

648 649 650 651 652 653 654 655

G Model

ARTICLE IN PRESS

ASOC 3307 1–11

C¸. Koc¸, I. Karaoglan / Applied Soft Computing xxx (2015) xxx–xxx

9

Table 5 Computational results for the test problems of Scenario-1. Instance

20c3sU1 20c3sU2 20c3sU3 20c3sU4 20c3sU5 20c3sU6 20c3sU7 20c3sU8 20c3sU9 20c3sU10

#Cust

20 20 20 20 20 20 20 20 20 20

#FS

3 3 3 3 3 3 3 3 3 3

Average

UB

LP

Root

B&C

Gap

CPU

Gap

CPU

#(35)

#(36)

Gap

CPU

#(35)

#(36)

1797.49 1574.78 1704.48 1482.00 1689.37 1618.65 1713.66 1706.50 1708.82 1181.31

38.0 26.2 40.3 29.3 30.7 31.4 35.5 35.8 40.6 13.9

0.0 0.0 0.0 0.0 0.1 0.0 0.0 0.0 0.0 0.0

15.6 11.6 28.0 18.3 13.7 14.9 15.0 9.0 20.3 0.8

0.6 0.4 0.1 0.1 0.6 0.2 0.3 0.4 0.5 0.1

70 76 38 24 57 48 61 63 64 27

257 114 40 32 247 77 134 199 182 25

0.0 4.7 0.0 9.5 0.0 2.7 9.8 0.0 8.3 0.0

172.3 3600.0 1789.0 3600.0 2165.5 3600.0 3600.0 1601.3 3600.0 2.3

3918 136,128 126,733 543,111 76,617 206,853 144,715 42,187 157,942 27

12,445 351,763 192,425 915,252 199,557 593,611 411,704 182,116 349,572 25

1617.71

32.2

0.0

14.7

0.3

53

131

3.5

2373.0

143,823

320,847

Table 6 Computational results for the test problems of Scenario-2. Instance

20c3sC1 20c3sC2 20c3sC3 20c3sC4 20c3sC5 20c3sC6 20c3sC7 20c3sC8 20c3sC9 20c3sC10

#Cust

20 19 12 18 19 17 6 18 19 15

#FS

3 3 3 3 3 3 3 3 3 3

Average

656 657 658 659 660 661 662 663 664 665 666

667

668 669 670 671

UB

LP

Root

B&C

Gap

CPU

Gap

CPU

#(35)

#(36)

Gap

CPU

#(35)

#(36)

1173.57 1539.97 880.20 1059.35 2156.01 2758.17 1393.99 3139.72 1799.94 2583.42

32.6 52.1 31.3 32.4 48.7 60.9 47.1 51.7 60.2 67.0

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

21.2 26.2 9.3 15.3 37.9 42.5 26.2 25.6 39.7 48.0

0.1 0.3 0.1 0.1 0.2 0.1 0.0 0.2 0.3 0.2

18 51 36 30 39 15 3 29 52 58

20 95 60 24 117 31 1 42 149 101

7.6 0.0 0.0 11.0 0.0 0.0 0.0 0.0 0.0 0.0

3600.0 1164.5 25.4 3600.0 2246.4 61.6 0.1 53.7 113.9 2067.5

388,839 21,531 1527 237,796 79,900 2160 3 860 2075 53,172

396,159 44,110 1793 485,153 324,194 8332 1 748 5153 222,441

1848.43

48.4

0.0

29.2

0.2

33

64

1.9

1293.3

78,786

148,808

for corresponding heuristic solution value and UB is optimal/best found solution value. Please note that, the computation times of these heuristics are less than 1 s. Table 4 shows that average percentage gap values are 14.35% and 5.98% for ECW and SA heuristic, respectively. For ECW, the average gap values are 11.18%, 13.15%, 16.20% and 16.87% for scenario 1, 2, 3 and 4, respectively. For SA heuristic, the average gap values are 2.59%, 7.52%, 8.67% and 5.14% for scenario 1, 2, 3 and 4, respectively. According to these results, although the solutions obtained by ECW heuristic are far from the optimal/best found solution proposed SA heuristic improves this solution within short computation time.

5.4. Results for the solution approach In this section, we analyze the performance of our B&C algorithm in terms of the number of test instances solved to optimality, computation time, and bound qualities of several components (i.e., LP relaxation, root node of enumeration tree, and overall

B&C algorithm). Bound quality is defined as the gap between corresponding lower bound (LB) and the optimal/best found solution (UB) and is evaluated as 100 * ((UB − LB)/UB). We implement our B&C algorithm using the SCIP v1.2.0 framework [45] and [46] as underlying LP solver, the default parameters were used, and a time limit of 1 h is set for each instance. Tables 5–8 summarize the computational results for each scenario, respectively. In these tables, the first four columns display the name of instance, the number of customers, the number of FSs and upper bounds obtained by the proposed B&C algorithm. Bold and italic entries in UB columns indicate that the corresponding problem solved to optimality in one hour. Successive columns give information about LP relaxation, root node of the enumeration tree and overall B&C algorithm. While the columns labeled Gap report the percentage gap, the columns labeled CPU show the solution time. The columns labeled #(35) and #(36) display the number of added valid inequalities (35) and (36), respectively. In these tables, it can be seen that the average lower bound percentage gap values of LP relaxation from FKK are obtained as 32.2%,

Table 7 Computational results for the test problems of Scenario-3. Instance

#Cust

S1 S1 S1 S1 S1 S2 S2 S2 S2 S2

20 20 20 20 20 20 19 20 16 16

2i6s 4i6s 6i6s 8i6s 10i6s 2i6s 4i6s 6i6s 8i6s 10i6s

Average

#FS

6 6 6 6 6 6 6 6 6 6

UB

LP

Root

B&C

Gap

CPU

Gap

CPU

#(35)

#(36)

Gap

CPU

#(35)

#(36)

1578.12 1397.27 1560.49 1692.32 1173.48 1633.10 1505.07 2431.33 2158.35 1585.46

26.4 25.0 28.8 35.2 13.3 49.5 38.4 52.1 41.7 30.0

0.0 0.0 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0

11.5 13.6 7.9 8.1 0.1 11.8 22.8 8.8 16.7 7.8

0.5 0.1 0.5 0.6 0.1 0.3 0.3 0.2 0.1 1.5

64 21 66 58 26 47 55 35 30 69

150 33 168 202 28 85 76 88 59 174

0.0 4.3 0.0 0.0 0.0 0.0 9.5 0.0 0.0 0.0

1626.6 3600.0 523.5 817.8 1.9 66.0 3600.0 1801.0 3.2 1.3

51,655 249,452 10,447 12,168 29 853 98,062 44,320 41 69

112,355 484,415 25,838 52,278 30 2824 157,808 121,929 140 174

1671.50

34.0

0.0

10.9

0.4

47

106

1.4

1204.1

46,710

95,779

Please cite this article in press as: C¸. Koc¸, I. Karaoglan, The green vehicle routing problem: A heuristic based exact solution approach, Appl. Soft Comput. J. (2015), http://dx.doi.org/10.1016/j.asoc.2015.10.064

672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690

G Model ASOC 3307 1–11

ARTICLE IN PRESS

10

C¸. Koc¸, I. Karaoglan / Applied Soft Computing xxx (2015) xxx–xxx

Table 8 Computational results for the test problems of Scenario-41. Instance

S1 S1 S1 S1 S1 S2 S2 S2 S2 S2

4i2s 4i4s 4i6s 4i8s 4i10s 4i2s 4i4s 4i6s 4i8s 4i10s

#Cust

20 20 20 20 20 18 19 20 20 20

Average

#FS

2 4 6 8 10 2 4 6 8 10

UB

LP

Root CPU

Gap

CPU

#(35)

#(36)

Gap

CPU

#(35)

#(36)

1582.21 1460.09 1397.27 1397.27 1396.02 1059.35 1446.08 1434.14 1434.14 1434.13

33.8 28.2 25.0 25.0 24.9 32.4 35.9 34.9 34.9 34.9

0.0 0.1 0.1 0.1 0.1 0.0 0.0 0.0 0.1 0.1

23.7 17.3 13.6 13.6 13.3 14.5 19.7 10.5 10.6 10.2

0.1 0.2 0.1 0.2 0.2 0.1 0.3 0.3 0.5 0.5

18 28 21 20 32 37 67 54 92 72

24 32 33 28 43 31 65 79 89 109

6.2 3.8 4.3 0.0 5.0 8.2 6.2 6.8 6.6 6.7

3600.0 3600.0 3600.0 2133.6 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0

652,305 398,654 249,514 144,411 146,158 182,052 93,266 104,626 61,921 49,624

1,020,534 632,741 484,537 191,305 265,565 381,171 146,004 152,519 87,563 66,906

1404.07

31.0

0.0

14.7

0.2

44

53

5.4

3453.4

208,253

342,885

730

48.4%, 34.0% and 31.0% for each scenario, sequentially. However, at the root node, these values are reduced to 14.7%, 29.2%, 10.9% and 14.7% without significant increase in solution time (i.e., less than 1 s). This result indicates that valid inequalities (35) and (36) are very efficient in improving the solution quality. Tables 5–8 also show that 22 out of 40 instances are solved to optimality using the proposed B&C algorithm. For the other 18 test instances that is not solved in optimally, the same results found as in [27,28] which reported the best G-VRP results in the literature. Average percentage gap values are obtained as 3.5%, 1.9%, 1.4% and 5.4% with the average solution time of 2373.0, 1293.3, 1204.1 and 3453.4 s for each scenario, respectively. Tables 5 and 6 indicate that test problems where customers are uniformly located (i.e., Scenario-2) are more difficult to solve than clustered test problems (i.e., Scenario-1). Average percentage gap values are obtained as 3.5% and 1.9% with the average solution time of 2373.0 and 1293.3 s for uniformly and clustered located test problems, respectively. Table 7 reports the results obtained from the test problems in Scenario-3 which are obtained by increasing the number of FSs from three to six on the selected test problems from Scenarios 1 and 2. As is seen in this table, the average gap value and CPU time are found as 1.4% and 1204.1 s while these values are 2.8% and 1935.1 when the original problems for each row are considered (see Tables 5 and 6), meaning that the problem is easy to solve when the number of FSs increases. In another experiment, we analyze the effect of the number of FSs for the same customer configuration. Table 8 shows the results obtained from the test problems in Scenario-4. These results are obtained by increasing the number of FSs from two to ten by increments of two on the selected test problems from Scenarios 1 and 2 (20c3sU4 and 20c3sC4, respectively). The results clearly show that the objective function value for the test problems derived from 20c3sU4 is improved (i.e., decreased from 1582.21 to 1396.02) when the number of FSs increases. However, there is some fluctuation on the objective function value for the test problems derived from 20c3sC4 since the some customers become feasible as the newly added fuel stations and the number of customers are increased. Therefore, the objective function value increases when the number of customers increases.

731

6. Conclusion

691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729

732 733 734 735 736 737 738

B&C

Gap

We have developed a heuristic based solution approach to solve the green vehicle routing problem (G-VRP). We have presented a mixed integer programming formulation and developed a B&C algorithm for the exact solution of the G-VRP. We have adapted several valid inequalities from literature to solve the G-VRP. We have also developed several separation algorithms for these inequalities and a simulated annealing heuristic to improve the quality of

solutions. The solution approach has evaluated in terms of optimality gap and computation time to reach the best solution on the various test problems derived from the literature. Our results clearly indicated that the proposed mathematical formulation in the B&C algorithm is effectively able to tighten the bounds. The results also show that the SA heuristic improves this solution within short computation time, however the solutions obtained by ECW heuristic are far from the optimal/best found solution. We realized that the G-VRP is easy to solve when the number of fuel stations increases. Our results show that optimal solutions for 22 out of 40 test instances with 20 customers are obtained within reasonable computation time. Furthermore, the computational results indicate that the simulated annealing heuristic is able to improve the performance of the B&C algorithm in finding high quality solutions. Concerning the future research directions, the proposed MIP formulation and B&C algorithm can be modified by considering more realistic aspects of the G-VRP including customer demands, time windows, multiple depots, etc. The proposed B&C algorithm implements a straightforward branching scheme (i.e., branching most fractional variable), and three set of valid inequalities. Further researches may also address how to affect different branching schemes and valid inequalities to the G-VRP. CO2 emissions or other greenhouse gas emission can be embedded into the objective function. Uncited references [47,48]. Acknowledgements The authors are thankful to the SCIP team for providing us the source code of branch-and-cut-and-price framework. The authors are also grateful to “IBM Academic Initiative” program to provide us with CPLEX callable library. Thanks are due to referees for their valuable comments. References [1] G.B. Dantzig, J. Ramser, The truck dispatching problem, Manag. Sci. 6 (1959) 80–88. [2] A.A. Juan, J. Faulin, R. Ruiz, B. Barrios, S. Caballé, The SR-GCWS hybrid algorithm for solving the capacitated vehicle routing problem, Appl. Soft Comput. 10 (2010) 215–224. [3] M.R. Khouadjia, B. Sarasola, E. Alba, L. Jourdan, E.G. Talbi, A comparative study between dynamic adapted PSO and VNS for the vehicle routing problem with dynamic requests, Appl. Soft Comput. 12 (2012) 1426–1439. [4] Y. Marinakis, G.R. Iordanidou, M. Marinaki, Particle swarm optimization for the vehicle routing problem with stochastic demands, Appl. Soft Comput. 13 (2013) 1693–1704.

Please cite this article in press as: C¸. Koc¸, I. Karaoglan, The green vehicle routing problem: A heuristic based exact solution approach, Appl. Soft Comput. J. (2015), http://dx.doi.org/10.1016/j.asoc.2015.10.064

739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763

Q4

764

765

766

767 768 769 770 771

772

773 774 775 776 777 778 779 780 781 782 783

G Model ASOC 3307 1–11

ARTICLE IN PRESS C¸. Koc¸, I. Karaoglan / Applied Soft Computing xxx (2015) xxx–xxx

784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836

[5] Z. Ursani, D. Essam, D. Cornforth, R. Stocker, Localized genetic algorithm for vehicle routing problem with time windows, Appl. Soft Comput. 11 (2011) 5375–5390. [6] J.-F. Cordeau, G. Laporte, M.W.P. Savelsbergh, D. Vigo, Vehicle routing, in: C. Barnhart, G. Laporte (Eds.), Transportation, Handbooks in Operations Research and Management Science, 1st ed., Elsevier Science Publishers, Amsterdam, 2007, pp. 367–428. [7] B.L. Golden, S. Raghavan, E.A. Wasil, The Vehicle Routing Problem: Latest Advances and New Challenges, Springer Science + Business Media, LLC, New York, 2008. [8] G. Laporte, Fifty years of vehicle routing, Transp. Sci. 43 (2009) 408–416. [9] P. Toth, D. Vigo, Vehicle Routing: Problems, Methods, and Applications, MOSSIAM Series on Optimization, Philadelphia, 2014. [10] A.C. McKinnon, M.I. Piecyk, Measurement of CO2 emissions from road freight transport: a review of UK experience, Energy Policy 37 (2009) 3733–3742. [11] A. Sbihi, R.W. Eglese, Combinatorial optimization and green logistics, Ann Oper Res 5 (2007) 99–116. [12] T. Bektas¸, G. Laporte, The pollution-routing problem, Transp. Res. B: Methodol. 45 (2011) 1232–1250. [13] A.A. Juan, J. Goentzel, T. Bektas¸, Routing fleets with multiple driving ranges: is it possible to use greener fleet configurations? Appl. Soft Comput. 21 (2014) 84–94. [14] C¸. Koc¸, T. Bektas¸, O. Jabali, G. Laporte, The fleet size and mix pollution-routing problem, Transp. Res. B: Methodol. 70 (2014) 239–254. [15] R. Eglese, T. Bektas¸, Green vehicle routing, in: P. Toth, D. Vigo (Eds.), Vehicle Routing: Problems, Methods, and Applications, MOS-SIAM Series on Optimization, Philadelphia, 2014, pp. 437–458. [16] E. Demir, T. Bektas¸, G. Laporte, A review of recent research on green road freight transportation, Eur. J. Oper. Res. 237 (2014) 775–793. [17] C. Lin, K.L. Choy, G.T.S. Ho, S.H. Chung, H.Y. Lam, Survey of green vehicle routing problem: past and future trends, Expert Syst. Appl. 41 (2014) 1118–1138. [18] S. Erdogan, E. Miller-Hooks, A green vehicle routing problem, Transp. Res. E: Logist. Transp. Rev. 48 (2012) 100–114. [19] J.F. Bard, L. Huang, M. Dror, P. Jaillet, A branch and cut algorithm for the VRP with satellite facilities, IIE Trans. 30 (1998) 821–834. [20] Y. Chan, S.F. Baker, The multiple depot, multiple traveling salesmen facility-location problem: vehicle range, service frequency, and heuristic implementations, Math. Comput. Model. 41 (2005) 1035–1053. [21] B. Crevier, J. Cordeau, G. Laporte, The multi-depot vehicle routing problem with inter-depot routes, Eur. J. Oper. Res. 176 (2007) 756–773. [22] C.D. Tarantilis, E.E. Zachariadis, C.T. Kiranoudis, A hybrid guided local search for the vehicle-routing problem with intermediate replenishment facilities, Inf. J. Comput. 20 (2008) 154–168. [23] I. Frade, A. Ribeiro, G. Gonc¸alves, A.P. Antunes, Optimal location of charging stations for electric vehicles in a neighborhood in Lisbon, Portugal, Transp. Res. Rec. J. Transp. Res. Board 2252 (1) (2011) 91–98. [24] T.D. Chen, K.M. Kockelman, M. Khan, Electric vehicle charging station location problem: a parking-based assignment method for Seattle, in: Transportation Research Board 92nd Annual Meeting (No. 13-1254), 2013. [25] Y.M. Nie, M.A. Ghamami, Corridor-centric approach to planning electric vehicle charging infrastructure, Transp. Res. B: Methodol. 57 (2013) 172–190. [26] J. Cavadas, G. Correia, J. Gouveia, Electric vehicles charging network planning, in: Computer-based Modelling and Optimization in Transportation, Springer International Publishing, 2014, pp. 85–100.

11

[27] M. Schneider, A. Stenger, D. Goeke, The electric vehicle-routing problem with Q5 time windows and recharging stations, Transp. Sci. (2014) (in press). [28] M. Schneider, A. Stenger, J. Hof, An Adaptive VNS Algorithm for Vehicle Routing Problems with Intermediate Stops, Technical Report LPIS-01/2014, 2014. ˜ G. Righini, G. Tirado, A heuristic approach for the green [29] Á. Felipe, M.T. Ortuno, vehicle routing problem with multiple technologies and partial recharges, Transp. Res. E: Logist. Transp. Rev. 71 (2014) 111–128. [30] C. Miller, A. Tucker, R. Zemlin, Integer programming formulations and traveling salesman problems, J. ACM 7 (1960) 326–329. [31] R.V. Kulkarni, P.R. Bhave, Integer programming formulations of vehicle routing problems, Eur. J. Oper. Res. 20 (1985) 58–67. [32] M. Desrochers, G. Laporte, Improvements and extensions to the Miller–Tucker–Zemlin subtour elimination constraints, Oper. Res. Lett. 10 (1991) 27–36. [33] I. Kara, G. Laporte, T. Bektas, A note on the lifted Miller–Tucker–Zemlin subtour elimination constraints for the capacitated vehicle routing problem, Eur. J. Oper. Res. 158 (2004) 793–795. [34] I. Kara, Two Indexed Polynomial Size Formulations for Vehicle Routing Problems, Bas¸kent University, Ankara/Turkey, 2008. [35] N.R. Achuthan, L. Caccetta, S.P. Hill, An improved branch-and-cut algorithm for the capacitated vehicle routing problem, Transp. Sci. 37 (2003) 153–169. [36] I. Karaoglan, F. Altiparmak, I. Kara, B. Dengiz, A branch and cut algorithm for the location-routing problem with simultaneous pickup and delivery, Eur. J. Oper. Res. 211 (2011) 318–332. [37] I. Karaoglan, F. Altiparmak, I. Kara, B. Dengiz, The location-routing problem with simultaneous pickup and delivery: formulations and a heuristic approach, Omega 40 (2012) 465–477. [38] G. Laporte, Y. Nobert, P. Pelletier, Hamiltonian location problems, Eur. J. Oper. Res. 12 (1983) 82–89. [39] L. Gouveia, A result on projection for the vehicle routing problem, Eur. J. Oper. Res. 85 (1995) 610–624. [40] A.N. Letchford, J.-J. Salazar-González, Projection results for vehicle routing, Math. Program. 105 (2006) 251–274. [41] K.A. Dowsland, J.M. Thompson, Simulated annealing, in: Handbook of Natural Computing, Springer, Berlin, 2012, pp. 1623–1655. [42] R. Banos, J. Ortega, C. Gil, A. Fernandez, F. De Toro, A simulated annealingbased parallel multi-objective approach to vehicle routing problems with time windows, Expert Syst. Appl. 40 (2013) 1696–1707. [43] G. Clarke, J.W. Wright, Scheduling of vehicles from a central depot to a number of delivery points, Oper. Res. 12 (1964) 568–581. [44] B. Suman, P. Kumar, A survey of simulated annealing as a tool for single and multiobjective optimization, J. Oper. Res. Soc. 57 (2005) 1143–1160. [45] T. Achterberg, SCIP – A Framework to Integrate Constraint and Mixed Integer Programming, Zuse Institute, Berlin, 2004. [46] IBM-ILOG, CPLEX 12.2 http://www-01.ibm.com/software/integration/ optimization/cplex-optimizer/. [47] R. Baldacci, E. Hadjiconstantinou, A. Mingozzi, An exact algorithm for the capacitated vehicle routing problem based on a two-commodity network flow formulation, Oper. Res. 52 (2004) 723–738. [48] T. Bektas¸, G. Erdo˘gan, S. Røpke, Formulations and branch-and-cut algorithms for the generalized vehicle routing problem, Transp. Sci. 45 (2011) 299–316.

Please cite this article in press as: C¸. Koc¸, I. Karaoglan, The green vehicle routing problem: A heuristic based exact solution approach, Appl. Soft Comput. J. (2015), http://dx.doi.org/10.1016/j.asoc.2015.10.064

837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890