The Accessibility Vehicle Routing Problem

Accepted Manuscript The Accessibility Vehicle Routing Problem O.J. Ibarra-Rojas, L. Hernandez, L. Ozuna PII:

S0959-6526(17)32553-2

DOI:

10.1016/j.jclepro.2017.10.249

Reference:

JCLP 11041

To appear in:

Journal of Cleaner Production

Received Date: 28 August 2017 Revised Date:

5 October 2017

Accepted Date: 22 October 2017

Please cite this article as: Ibarra-Rojas OJ, Hernandez L, Ozuna L, The Accessibility Vehicle Routing Problem, Journal of Cleaner Production (2017), doi: 10.1016/j.jclepro.2017.10.249. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

The Accessibility Vehicle Routing Problem O.J. Ibarra-Rojasa,∗, L. Hernandezb , L. Ozunac a Universidad

Aut´ onoma de Nuevo Le´ on, Faclutad de Ciencias F´ısico-Matem´ aticas Aut´ onoma de Nuevo Le´ on, Facultad de Ciencias Qu´ımicas c Universidad Aut´ onoma de Nuevo Le´ on, Facultad de Ingenier´ıa Mec´ anica y El´ ectrica

RI PT

b Universidad

Abstract

D

M AN U

SC

In a distribution process where the demand relates to essential products or services, is important to consider the access for people to fulfill their needs. In particular, for land use and urban transportation planning, accessibility relates to appropriately allocating opportunities to satisfy a demand or provide a service considering the cost of mobility. Measuring accessibility is a challenging task, indeed, it depends on the context of the study and has not been properly considered in the definition of vehicle routing problems, which are commonly used to represent distribution processes. In the study reported here, we addressed a vehicle routing problem to optimize accessibility based on six indicators: the number of zones with access to opportunities with delimited mobility, the number of zones covered by the route, the cost of travel, the distance to the nearest opportunity, the number of opportunities, and geographical disaggregation. We defined a mixed-integer linear formulation for the proposed problem that we used to show the potential benefits of our approach compared with a maximum coverage vehicle routing problem for small instances. In turn, we designed an iterated local search algorithm and analyzed its efficiency according to a benchmark of randomly generated instances. Numerical results show that we obtain high-quality solutions for acceptable computational times.

3 4 5 6 7 8 9 10 11 12 13 14 15

Given the world’s ever-increasing human population, the sustainability of planning and of the implementation of processes, whether in manufacturing, public services, or land use, among other fields, is important to guarantee equilibrium between populations and their ecosystems. To attain the desired characteristics of such activities, different indicators have been considered. Juwana et al. (2012), for example, have reviewed indicator-based water sustainability assessments, in which access to water resources for people was an important indicator. In general, improving accessibility means guaranteeing more and better-allocated opportunities for people to obtain a service or meet a demand without relying as much on mobility and is therefore undoubtedly important from a socioeconomic perspective. Indeed, measures of accessibility have been applied to land use and transportation analysis, among other studies on socioeconomic factors. However, perspectives on accessibility differ depending on the context, and to the best of our knowledge, no consensus on the definition of the term or best measures for it exist. According to Geurs and van Wee (2004), any measure of accessibility should consider the impedance, or cost, imposed upon an individual to cover the distance between an origin and a destination, as well as the number and locations

EP

2

1. Introduction

AC C

1

TE

Keywords: Accessibility, vehicle routing problem, mixed-integer programming, iterated local search

∗ Corresponding

author Email addresses: [email protected] (O.J. Ibarra-Rojas), [email protected] (L. Hernandez), [email protected] (L. Ozuna)

Preprint submitted to Elsevier

October 24, 2017

ACCEPTED MANUSCRIPT

41

2. Related literature

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

42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

SC

22

M AN U

21

D

20

TE

19

Measures of accessibility measure commonly consider different elements or indicators based on the context. For example, Liu et al. (2015) analyzed the impact of different factors on park visitation among citizens, including the accessibility of green areas, defined as a weighted sum of the indicators of the number of parks within a given radius, the walking time to the nearest park, and the number of turns along the shortest path to the nearest park. More recently, S´anchez-Garc´ıa et al. (2017) proposed a method based on geographical information systems to evaluate the location of wood-fired power plants in relation to the availability of biomass feedstock, both geographically and in terms of accessibility. Their indicators related to accessibility used to evaluate the locations of wood-fired plants were the distances from biomass feedstock to major roads and the varying levels of accessibility to direct the supply of wood fuel based on the interaction of slopes in the terrain, road networks, population centers, government restrictions in logging operations, and legal limitations in protected areas such as national parks, natural parks, sites of community importance, and special protection areas. As also done by Liu et al. (2015), their indicators were combined in a measure of cost for the location of wood-fired plants. We consider accessibility to optimize distribution processes that can be represented with a VRP, particularly one with coverage indicators as part of an arguably global measure of accessibility, which, to the best of our knowledge, has not been previously addressed by researchers. Vehicle routing problems with route length constraints are formulated as orienteering problems and maximal covering tour problems. The classic orienteering problem considers an objective function of

EP

18

AC C

17

RI PT

40

of opportunities along the way, ideally in light of the spatial distribution of the services that meet the demands for those opportunities. Handy and Niemeier (1997) have posited that four interrelated elements should be considered by any measure of accessibility: the degree of spatial disaggregation, the definition of origin and destinations, travel impedance, and the attractiveness of a zone due to the provided services in the area. As Makr´ı and Folkesson (1999) have argued, a function of travel impedance regarding accessibility can be defined in terms of the distance or time estimated in light of straight-line distance, network distance, network models simulating travel demand, or field surveys of actual driving times, if not a combination of those elements. In short, the greater the travel time or distance to reach a destination, the less accessible that destination is. Even when meeting demands is critical in distribution processes—for instance, in supplying humanitarian relief—and accessibility is therefore crucial, previous studies have not consider the concept (see reviews of humanitarian relief of Altay and Green III, 2006; de la Torre et al., 2012). In response, this study propose an accessibility measure based on the definition of accessibility articulated by Bertolini et al. (2005) which is related to amount and diversity of places that can be reached within a given travel time and/or cost. This measure is integrated in a Vehicle Routing Problem with route length constraints. Our major contributions are the definition of a new variant for the VRP including accessibility, a mathematical formulation for the proposed problem, a solution algorithm capable of obtaining good quality solutions in acceptable computational times, and a comparative analysis that exhibits the relevance of considering the accessibility concept in routing. The structure of this paper is as follows. Section 2 presents related literature to our study. Section 3 introduces our Accessibility Vehicle Routing Problem (denoted as AVRP) as well as its mathematical formulation. Section 4 introduces an Iterated Local Search to obtain high quality solutions for large instances of the AVRP. Section 5 presents a comparison of AVRP and a Maximal Covering Vehicle Routing problem, as well as an analysis of the efficiency of our ILS algorithm. Finally, Section 6 presents our conclusion and further research.

16

2

ACCEPTED MANUSCRIPT

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

RI PT

66

SC

65

M AN U

64

D

63

TE

62

EP

61

zone coverage based on scores for zones visited and penalties for zones not visited by routes. The latter problem gives a single score for each zone and imposes a hard constraint for route time and length (see review of Vansteenwegen et al., 2011), which in turn allows it consider a weighted objective function for the number of zones on the route. A variant of the problem is the team orienteering problem, in which the goal is to determine a number of paths P , each limited by a maximum route time such that the total collected score is maximized (Tang and Miller-Hooks, 2005). Although those problems can also consider time windows (Tricoire et al., 2010; Vansteenwegen et al., 2009), the objective function is a sum score for visited zones, which is a rather simple objective compared to elements that should be considered by a measure of accessibility. A problem more pertinent to our study is the generalized orienteering problem, in which the total score is defined by nonlinear functions of vertices visited, for a certain combination of vertices can produce a higher or lower score than the sum of individual scores (Geem et al., 2005; Wang et al., 1995, 2008). Some studies on the generalized orienteering problem have focused on sightseeing activities; for example, Wang et al. (2008) defined multiple attributes for each touristic attraction, and the score was not a single value but consisted of multiple attributes considered in a weighted objective function based on the preferences of the different attributes. The authors proposed a genetic algorithm to solve large instances of the problem. The maximal covering tour problem is also related to the concept of accessibility. This problem was proposed by Current and Schilling (1994), who assumed that people can travel from their origin zone to other zones to obtain a service or meet a demand. That problem also involves a coverage radius for each visited zone, although the tour is limited to a maximum number of visits. As such, the goal is to minimize the sum of distances between zones not covered and their nearest coverage radius. In a sense, Current and Schilling (1994) considered indicators of accessibility, for their model can maximize the number of zones covered, as well as guarantee nearby opportunities. However, their model does not consider attributes such as spatial disaggregation or the number of opportunities. More recently, Flores-Garza et al. (2015) presented a generalization of the MCTP that considers a multivehicle case and assumes mandatory coverage for zones impossible to cover by the routes and zones whose coverage is optional. The goal of their model is to minimize the arrival time at all locations visited considering a limited route length, and to obtain good-quality solutions, the authors proposed a greedy randomized adaptive search procedure heuristic. In another study, Nolz et al. (2010) examined a multiobjective coverage tour problem in the context of humanitarian relief, with the proposed objectives of a weighted sum of distances from zones not covered to their nearest coverage radius, in which weights represent the number of inhabitants in each zone; the maximum visit time at each zone along the route in consideration of service time; and route duration. Since their problem is intractable for commercial solvers, the authors implemented a hybrid method based on genetic algorithms and a variable neighborhood search to obtain quality solutions for large instances. Although facility location problems also define coverage objective functions, as with VRPs, the common considerations are radius of coverage in terms of an acceptable service distance or the minimum distance from the farthest clients to service points. Constraints of the problems include single or multiple sources for the satisfaction of demands, limited capacity, mandatory coverage for specific clients, mandatory constraints upon closeness, the allowance or prohibition of partial coverage in terms of demand, and hierarchical location (Farahani et al., 2012). To the best of our knowledge, our study is the first to apply the concept of accessibility to a VRP.

AC C

60

3

ACCEPTED MANUSCRIPT

102

3. Accessibility Vehicle Routing Problem

103

3.1. Assumptions and input

111 112 113 114 115 116

117 118

119

120 121

122 123 124 125 126 127 128 129

RI PT

110

SC

109

M AN U

108

• A zone i is directly covered if it is in a route and indirectly covered if not but nevertheless it is in neighborhood N (j) of directly covered zone j. • If neither directly nor indirectly covered, then zone i is not covered. • A zone j represents an opportunity for zone i not covered only if zone j is directly covered and j ∈ A(i). Figure 1 illustrates coverage sets N (i), accessibility sets A(i), and the classification of zones using a graph. Two routes depart from one origin (represented with a square node of the graph). Sets N (i) are represented by black dashed lines and, for example, N (5) = {6}. Sets A(i) are represented by red doted lines, for example, A(7) = {4, 5, 6}. The solution shows that zone 6 is indirectly covered by zone 5 and that zone 7 is not covered but presents two opportunities to have its demand met (i.e., by zones 4 and 5). Zones 8, 9, and 10 are not covered, and because they have only one opportunity to have the demand met (i.e., in zone 1), they can be considered to constitute a small cluster of zones not covered—a situation that we would prefer to avoid.

D

107

TE

106

EP

105

We assume a set of zones with unitary demand and the number of vehicles with unlimited capacity seeking to cover those zones, as well as known travel times for each pair of zones connected by a direct road (independent of the vehicles). Moreover, the zones should be visited by at most one vehicle, it is possible to cover zones indirectly if they are near other zones covered by the route (proposed by Current and Schilling, 1994)—that is, we consider coverage neighborhoods—and a zone not covered could have the demand met by another zone within a predefined limit of mobility distance (i.e., we consider accessibility neighborhoods). Now, we introduce our notation, I is the set of zones, in which 0 represents the origin of the routes and K the number of homogeneous vehicles. The travel time from zone i to zone j is represented as tij for each pair of zones i and j connected with a direct road, and T is the time limit of the routes. N (i) is the set of zones covered indirectly if zone i is visited by a route, and A(i) is the set of accessible zones for zone i from which zone i can have the demand met. Besides the previous sets and parameters, we define the following classification of zones.

AC C

104

7 6

5

4 8

o 3

10 1 9

2

Figure 1: Example of two routes in a network with 10 zones that illustrate coverage and accessibility neighborhoods.

4

ACCEPTED MANUSCRIPT

133 134 135 136 137 138

Any solution for a vehicle routing problem with constraints of limited length of routes length involves sets I1 and I2 of zones directly covered and not covered, respectively. In the example in Figure 1, those sets are I1 = {1, 2, 3, 4, 5} and I2 = {7, 8, 9, 10}, whereas zones indirectly covered are defined by I − (I1 ∪ I2 ∪ {0}) = {6}. Based on sets I1 and I2 , we define different indicators of accessibility to merge into a single measure. First, we are interested in the number of zones with access to the demand either being directly covered or having at least one opportunity within their mobility radius. Then, we define the following indicator aai taking the value of 1 if the zone has access to obtain a service or meet a demand, and 0 otherwise.

139 140 141 142 143

1 0

147 148 149 150

 

P

f (tij )

j∈I1 ∩A(i)

P

f (tij )

j∈A(i)

if i ∈ / I2 ,

(2)

if i ∈ I2 .

D

TE

146

   1

Note that ati takes a value within (0,1), in which 0 and 1 represent the worst and best accessibility values, respectively. In particular, ati takes the value of 1 if zone i is covered or if not covered but all of its accessible zones in A(i) are directly covered. However, if ati takes the value of 1, it does not differentiate the number of opportunities of zone i compared with other zones i0 for which ati0 is also 1, and the same happens with the distance to opportunities. Thus, we propose accessibility indicator aoi = |I1 ∩ A(i)| to consider the number of opportunities for each zone i not covered, as well as the indicator ani = min{tij }, which represents the distance from a zone i not covered to the nearest directly

EP

145

(1)

Other than indicator aai , we are interested in the number ac of covered zones, which is simply defined as ac = |I1 |. We also define our function of travel impedance f (tij ) to access from zone i to zone j. That type of function is commonly defined as the inverse power of the distance or time from one zone to another (Makr´ı and Folkesson, 1999). Thus, we propose the impedance function f (tij ) = t1ij to define the following accessibility indicator ati for each zone i in terms of travel impedance.

ati =

144

if i ∈ I1 or j ∈ I1 for some j ∈ A(i), otherwise .

SC

aai =

(

RI PT

131 132

3.2. Accessibility measure

M AN U

130

j∈I1

152

153 154 155

covered zone j—a distance that we seek to minimize. Lastly, we defined an accessibility indicator in terms of the spatial disaggregation, as = min {tij }, which we maximize in order to avoid a situation in i,j∈I2

which zones not covered were near each other. Our accessibility measure A(I1 , I2 ) for sets I1 , and I2 was thus defined as the following weighted sum of indicators of accessibility, for which parameters wk represent the weights.

AC C

151

A(I1 , I2 ) = w1

X

aai + w2 ac + w3

i∈I

X i∈I

ati + w4

X

i∈I2

ani + w5

X

aoi + w6 as

(3)

i∈I2

158

Based on the above, our AVRP determines which zones of a set of given zones I have to be visited by a number of routes K with a limited length such that the accessibility measure A(I1 , I2 ) based on coverage, number and location of opportunities, and spatial disaggregation is maximized.

159

3.3. Mixed-integer linear program for the AVRP

156 157

160

To define our mathematical formulation, we define the following decision variables.

5

ACCEPTED MANUSCRIPT

zi

vij

aai

otherwise. if zone i is covered either directly o indirectly by a route, otherwise. if zone i is directly covered by a route, otherwise.

RI PT

yi

if a route visits zone j after visiting zone i,

if zone i is not covered and j is the nearest zone with zj = 1, otherwise.

SC

xij

 1 = 0  1 = 0  1 = 0  1 = 0  1 = 0

if zone i is covered or has opportunities within is mobility radius, otherwise.

M AN U

ati : indicator based on travel impedance and location of opportunities for zone i. ani : indicator for distance from zone i not covered to the nearest zone j with zj = 1. aoi : indicator representing the number of opportunities for zone i not covered. as : indicator based on spatial disaggregation of zones not covered. ui : route time before visiting zone i with a vehicle (u0 = 0). Then, our mixed-integer linear program is as follows, in which parameters wi define the weights of accessibility indicators whereas I 0 is defined as I 0 = I − {0}. max

X

EP

subject to:

AC C

163

(w1 aai + w2 yi + w3 ati + w4 ani + w5 aoi ) + w6 as

TE

i∈I 0

D

161 162

6

(4)

ACCEPTED MANUSCRIPT

X

x0i =

i∈I 0

(6)

X

xij =

∀i ∈ I 0

xji

j∈I

RI PT

j∈I

j∈I

(5)

xi0 ≤ K

i∈I 0

X

xi0

i∈I 0

X

X

X

xji ≤ 1

ui + tij ≤ uj + T (1 − xij )

j:i∈N (j)

yi ≤ zi +

X

zj

j:i∈N (j)

172 173 174 175

(13)

∀i ∈ I 0 , j ∈ A(i)

(14)

(18)

∀i ∈ I 0

(19)

∀i ∈ I 0

(20)

as ≤ tij + M2 (yi + yj )

∀i ∈ I , j ∈ I , j 6= i

(21)

∀i ∈ I, j ∈ I : j 6= i

(22)

j∈A(i)

zj

D

EP

X

∀i ∈ I

TE

vij tij

0

xij , yi , zi , vij , aai ∈ {0, 1}; ati ∈ (0, 1); aoi , ap ∈ Z+ ∪ {0}; ui , ani , as ≥ 0,

171

∀i ∈ I 0

∀i ∈ I 0

X

yi

aoi ≤ |A(i)|(1 − yi )

170

(12)

(17)

aoi ≤

169

∀i ∈ I 0

∀i ∈ I 0

vij

j∈I 0 :j6=i

167

(11)

(16)

X

j∈A(i) P zj tij j∈A(i) P 1 + tij j∈A(i)

ani :=

168

∀i ∈ I 0

∀i ∈ I 0

ati ≤

166

(10)

(15)

vij = 1 − yi

aai := yi +

165

(9)

0

0

j∈I 0 :j6=i

164

∀i ∈ I, j ∈ I 0 , j 6= i

AC C

vij ≤ zj X

zj ≤ M1 yi

M AN U

X

(8)

∀i ∈ I

j∈I

zi +

∀i ∈ I 0

SC

ui ≤ T X zi := xij

(7)

0

Objective function (4) represents our accessibility measure. Constraints (5) and (6) limit the number of vehicles, whereas (7) and (8) are flow balance restrictions. The inequalities (9) define the accumulated travel time ui before the visit to each zone i on a route whereas the inequalities (10) guarantee that each variable ui for zone i ∈ I 0 is less than or equal to the maximum route time T . Notice that we do not consider the travel time of vehicles returning to the depot. The definition of zi is given by (11). Constraints (12) force variable yi to be 1 when zone i is either directly or indirectly covered (M1 can be set to |N (i)| + 1) and inequalities (13) force variable yi to be 0 when zone i is not covered. Constraints (14) and (15) guarantee that each zone not covered has only one nearest opportunity j. Equations (16) guarantee that aai takes the value of 1 if the zone is either covered (yi = 1) or has at least one opportunity j within its mobility radius A(i), and 0 otherwise. The inequalities (17) and the objective function (4) define the value of indicator ati . Notice that by maximizing our accessibility measure, we try to define near opportunities for each zone i not covered, i.e., to activate variables zj with small values of tij when 7

ACCEPTED MANUSCRIPT

176 177 178 179

yi = 0. Equation (18) defines the variables ani , whereas inequalities (19) and (20), along with objective function (4), define aoi as the number of opportunities for zone i if not covered and 0 if covered. Lastly, Constraint (21), along with objective function (4), defines the minimum distance between zones not covered, in which M2 can be set to max{tij }. (i,j)

183

4. Iterated Local Search for the AVRP

186 187 188 189 190

.

s⇤

TE

Objective value

191

A simple, quick tool that can be improved by increasing the quality of each of its components, an ILS relies on evading local optimality by way of perturbation (i.e., search diversification) while increasing the quality of the local search (Louren¸co et al., 2010). Given an initial incumbent solution s, which can be a local optimum, we implement a perturbation to obtain an intermediate feasible solution s0 . Subsequently, our local searches could find a new local optimum s0∗ . If the new solution meets an acceptance criterion (typically, the objective function), then the incumbent solution would be s0∗ ; otherwise, s would remain the incumbent. We iterate until the algorithm reaches a maximum number of iterations (these steps are exhibit by Figure 2). Now, we detail all the components of our algorithm.

SC

185

M AN U

184

.

s0⇤

Pe

rtu

rb

at io

D

181

RI PT

182

By using commercial solvers, our proposed formulation is intractable since we only obtain optimal solutions for small instances. We therefore develop an iterated local search (ILS) metaheuristic to obtain high-quality feasible solutions for large instances of the AVRP.

180

n

.

s0

Solution space

193 194 195 196 197 198 199 200 201 202 203 204 205

4.1. Constructive algorithm

Our proposed constructive algorithm is based on the savings methods proposed by Clarke and Wright (1964). In particular, we define a list LS of pairs (i, j) where zone i is the first or the last zone in a route k whereas zone j is not visited by any route. This list LS is ordered decreasingly in terms of the saving parameter sij = t0i + t0j − λtij , where λ is a route shape parameter (see Yellow, 1970). Moreover, given a solution s and a pair (i, j) ∈ LS, we define a insertion procedure InsertSavings(s, i, j) which determines if zone j is inserted at the beginning or at the end of the route k visiting i, as well as the direction of the route, which is necessary since we do not consider the travel time of vehicles returning to the depot. Figure 3 exhibits an example of solution s consisting of a route visiting zones i and i0 , and a zone j to be inserted in that route. Cases (a)–(d) show the four solutions that are explored with our proposed operator InsertSavings(s, i, j), from which we choose the feasible solution with the minimal route length. We also define an alternative insertion procedure to explore each position of all routes. Given a solution s, and a zone i not visited by any route, the operator Insert(s, i) inserts zone i in the position

AC C

192

EP

Figure 2: Schematic view of an iterated local search algorithm.

8

ACCEPTED MANUSCRIPT Solution s :

j

i0

i

(a)

(b)

j

(c)

j

i

i0

(d)

j

i0

i

j

i

i0

RI PT

i

i0

Figure 3: Solution s and solutions (a)–(b) that are explored when implementing the operator InsertSavings(s, i, j) for a specific pair (i, j) ∈ LS.

226

4.2. Local search algorithms

214 215 216 217 218 219 220 221 222 223 224

227 228 229

230 231

232 233

234 235 236

M AN U

212 213

D

211

TE

209 210

EP

208

AC C

207

SC

225

of a route k that yields the minimal increment of route time, for which we consider both feasible and infeasible insertions along our ILS. We recall that our accessibility measure depends on the zones visited by routes. Then, we define a set L of different ordered lists to use in the insertion of zones leading to different solutions. In particular, L has three lists, where the first list order zones i ∈ I 0 decreasingly in terms of the value of |{j ∈ I 0 : i ∈ A(j), yj = 0}| but breaking ties by allocating first in L the zone i with maximal cardinality |N (i)|; the second list order zones i ∈ I 0 decreasingly in terms of the value of |N (i)| but breaking ties by allocating first in L the zone i with maximal cardinality |{j ∈ I 0 : i ∈ A(j), yj = 0}|; and finally, the third list order zones i ∈ I 0 decreasingly in terms of the value of |{j ∈ I : i ∈ A(j), yj = 0}| but breaking ties by allocating first in L the zone i with maximal value of t0i . Based on the above, we define a constructive algorithm Constructor(λ, L) consisting of three phases detailed in Algorithm 1. Similarly to the classic savings method, Phase I of our constructive algorithm defines K feasible routes visiting a single zone that is chosen from the ordered list L. In, Phase II, we implement our operator InsertSavings(s, i, j) for solution s and pairs (i, j) ∈ LS, such that, all pairs (i, j) ∈ LS where zone j has no access to opportunities (i.e., m = 0) are explored first, followed by pairs (i, j) ∈ LS with zone j not covered but with at least one opportunities (i.e., m = 1) and finally, pairs (i, j) ∈ LS with zone j indirectly covered (i.e., m = 2). Next, Phase III implement feasible operators Insert(s, i) based on the ordered list L. To obtain our initial solution, algorithm Constructor(λ, L) is implemented for different pairs (λ, L) and then, we choose the solution with the best value of the objective function.

206

We propose several greedy local searches based on the following operators, some of them have been successfully implemented to solve different SVRPs (e.g., Aras et al., 2011; Palomo-Mart´ınez et al., 2017; Wang et al., 2017). • Exchange(s, k, i, j): inserts zone j at the position of zone i on route k and removes zone i of that route only if that movement is feasible. • Reorder(s, k, i, j, m): cuts the segment of route k starting from zone i and ending in zone j and reinserts it after zone m, for which we consider both feasible and infeasible moves. • Swap(s, k, i, k 0 , i0 ): removes zone i from route k and inserts it at the position of route k 0 that yields the minimal increment on the route length. Similarly, we remove i0 from route k 0 and insert it in route k, for which we consider feasible swaps.

9

ACCEPTED MANUSCRIPT

Algorithm 1: Constructor(λ, L):

RI PT

Initialize empty solution s and decision variables equal to 0; Create savings list LS based on λ, T imek = 0, and k = 1; Phase I: while k ≤ K do Choose first zone i in list L with aai = 0 and t0i ≤ T ; Insert i at the beginning of the route, i.e., update decision variables x0i = xi0 = aai = yi = zi = 1, aaj = yj = 1 for each j ∈ N (i), and aaj = 1 for each j ∈ I : i ∈ A(j), and T imek = T imek + t0i ; Update list L = L − ({i} ∪ N (i)) and update k = k + 1;

M AN U

SC

Phase II: for m = 0, 1, 2 do Take first pair (i, j) ∈ LS such that aaj + yj + zj = m; s ← InsertSavings(s, i, j); if zj = 1 then Update route k (the one visiting i), and list LS according to the operator InsertSavings(s, i, j); Update decision variables aaj = yj = zj = 1, aaj 0 = yj 0 = 1 for each j 0 ∈ N (j), and aaj 0 = 1 for each j 0 ∈ I : j ∈ A(j 0 ), and list L = L − ({j} ∪ N (j)); else LS = LS − {(i, j)};

TE

D

Phase III: while T imek < T for some route k or L 6= ∅ do Take first zone i ∈ L; s ← Insert(s, i) considering only feasible movements; if zi = 1 then Update route k according to the operator Insert(s, i); Update decision variables aai = yi = zi = 1, aaj = yj = 1 for each j ∈ N (i), and aaj = 1 for each j ∈ I : i ∈ A(j);

238

239 240

241 242 243 244

245 246 247 248 249 250

• P osition(s, k, i, k 0 ): removes zone i from route k and insert it at the position of route k 0 that yields the minimal increment on the length of that route, for which we consider feasible moves. • Cut(s, k, i, k 0 , i0 ): cuts the route segment starting from i (i0 , resp) to the last zone on route k (k 0 , resp) and reinsert it in route k 0 (k, resp), for which we consider feasible operations.

AC C

237

EP

Update list L = L − {i};

• F easibility(s): creates a list Lk of zones visited for each infeasible route k and order it decreasingly based on the value of the objective function if zone i is removed. Next, for each infeasible route k, selects txhe first zone i in Lk then, removes i from route k and updates Lk = Lk − {i} and route k until that route reaches feasibility. Subsequently, we define LocalSearchI, a greedy algorithm that implements the operator Insert() iteratively over an ordered list L of zones not covered and indirectly covered, details of which appear in Algorithm 2. LocalSearchI implements only feasible insertions and a necessary reordering procedure because, if zone i is inserted, then yi takes the value of 1 and {j ∈ I : i ∈ A(j), yj = 0} is updated, which affects the order of list L. We also propose LocalSearchIF described by Algorithm 3, which consists of inserting a zone along a 10

ACCEPTED MANUSCRIPT

252

route, implementing feasibility recovery to obtain a feasible solution, and updating the incumbent when an improvement in the objective function is identified.

SC

251

RI PT

Algorithm 2: LocalSearchI(s): Create list L of indirectly covered and not covered zones ordered increasingly in terms of the value of yi (zones not covered first) but breaking ties by allocating first in L the zone i with the maximal cardinality of set {j ∈ I : i ∈ A(j), yj = 0} (zones not covered for which i could be an opportunity); while L 6= ∅ do Take first element i of L; Implement feasible Insertion(s, i); Update zones: L = L − {i};

M AN U

Algorithm 3: LocalSearchIF (s):

TE

D

improve = true; while improve do improve = false; for each route k and each zone i with zi = 0 do s0 ← Insert(s, i) but considering only route k; s0 ← F easibility(s0 ); if s0 improves the objective value then s ← s0 ; improve =true; else s0 ← s;

255

EP

254

Additionally, we define algorithm LocalSearchE, a greedy local search that implements the operator Exchange() in order to reduce route lengths and increase the value of the objective function by considering different zones on the routes. Details of the local search are described by Algorithm 4.

AC C

253

Algorithm 4: LocalSearchE(s): improve = true; while improve do improve = false; for each route k and each pair of zones (i, j) with yi = 1 and yj = 0 do s0 ← Exchange(s, k, i, j); if s0 improves the objective value then s ← s0 ; improve =true; else s0 ← s;

11

ACCEPTED MANUSCRIPT

258 259 260 261

Algorithm 5: LocalSearchS(s):

M AN U

improve = true; while improve do improve = false; for each route k and zone i in route k do for each route k 0 6= k and zone i0 in route k 0 do s0 ← Swap(s, k, i, k 0 , i0 ); if s0 diminish the sum of route lengths then s ← s0 ; improve = true; else s0 ← s;

262

RI PT

257

Finally, we propose local searches to diminish the route length so the value of the objective function can be increased by implementing insertion procedures. We first define LocalSearchS, which implements operator Swap() for each 4-tuple (s, k, i, k 0 , i0 ) and updates the solution if it diminishes the sum of the lengths of routes k; details appear in Algorithm 5. Analogously to LocalSearchS, we define two greedy algorithms, LocalSearchP and LocalSearchC, to implement the operators P osition() and Cut() for tuples (k, i, k 0 ) and (k, i, k 0 , i0 ), respectively.

SC

256

4.3. Perturbation procedure and acceptance criterion

270

4.4. Proposed Iterated Local Search

267 268

TE

266

EP

264 265

D

269

The perturbation procedure of an ILS must be an operator irreversible by local searches. The intensity of this perturbation bias the algorithm to diversification or intensification. We propose the operator P erturbation(s), which consists of implementing the operators Reorder(s, k, i, j, m) for randomly selected 4-tuple (k, i, j, m) and F easibility(s) to obtain a feasible solution. Having examined infeasible solutions and recovering feasibility in the local searches and the perturbation procedure, we define the acceptance criterion Acceptance(s, s0 ) based only on feasibility and the value of the objective function, for which we select the feasible solution with the highest value of accessibility.

263

277

5. Experimental stage

272 273 274 275

278 279 280

AC C

276

To obtain high-quality solutions for the AVRP, we define—in a preliminary stage—the order of local searches to be sequentially implemented in our ILS. The Algorithm 6 shows the detailed steps where we can notice that the local searches are used in three phases: improving value of objective function with LocalSearchIF and LocalSearchE; diminishing route length with LocalSearchS, LocalSearchP , and LocalSearchC; and improving value of objective function by implementing LocalSearchI. Our ILS stops when a number IterLimit of iterations without any improvement in the objective function is reached.

271

In this section, we present a benchmark of randomly generated instances. Moreover, we show that our metaheuristic outperforms the solver of CPLEX for large instances and that we obtain solutions congruent with the concept of accessibility.

12

ACCEPTED MANUSCRIPT

Algorithm 6: ILS():

281

282 283 284 285 286

M AN U

SC

RI PT

Define s as the best solution of Constructor(λ, L) among all pairs (λ, L); s ← LocalSearchIF (s) and s ← LocalSearchE(s); s ← LocalSearchS(s), s ← LocalSearchP (s), and s ← LocalSearchC(s); s ← LocalSearchI(s); Iter = 0; while Iter < IterLimit do s0 ← P erturbation(s); s0 ← LocalSearchIF (s0 ) and s0 ← LocalSearchE(s0 ); s0 ← LocalSearchS(s0 ), s0 ← LocalSearchP (s0 ), and s0 ← LocalSearchC(s0 ) ; s0 ← LocalSearchI(s0 ); s ← Acceptance(s, s0 ); if s0 6= s then Iter = Iter + 1;

5.1. Instances

We design different types of instances based on the number of zones |I|, the number of routes K, the maximum length T for each route, the size of the grid in which zones are randomly generated, a coverage radius rN to define set N (i) for each i ∈ I 0 (i.e., j ∈ N (i) if tij ≤ rN ), and the accessibility radius rA to define A(i) for each i ∈ I 0 (i.e., j ∈ A(i) if tij ≤ rA ). Table 1 shows the different instance sizes and for each of which we randomly generate 30 instances.

D

Table 1: Instances types for the AVRP.

|I| 15 16 18 19 60 61 120 121

287 288 289 290 291 292 293 294 295 296 297

AC C

EP

TE

Instance type A1 A2 B1 B2 C1 C2 D1 D2

K 2 2 2 2 5 5 12 12

T 60 60 80 80 120 120 180 180

size of the grid 60 × 60 60 × 60 60 × 60 60 × 60 120 × 120 120 × 120 240 × 240 240 × 240

rN 8 8 8 8 8 8 12 12

rA 16 16 16 16 18 18 30 30

For instances A1, B1, C1, and D1, all zones are randomly generated within the grid. In the case of instances types A2, B2, C2, and D2, the depot is randomly generated within the grid and then, |I|−1 3 zones are randomly generated within each one of the three regions in the grid layout illustrated by Figure 4. Notice that the area of the three regions is the same and by using the latter grid layout, we are able of representing scenarios were may exist groups of zones. Finally, all instances satisfy A(i) 6= ∅ for all i ∈ I 0 . We executed all experiments on a MacPro with 2×2.4 GHz quad core Intel Xeon and 20 GB of RAM. We solved MILPs by using CPLEX 16.2 with default settings and a stopping criterion of 3 h of computational time. After a preliminary experimental stage with the AVRP MILP, we fixed weight parameters as w1 = 1000, w2 = 250, w3 = 100, w4 = −1, w5 = 10, and w6 = 0.01. For our constructor algorithm, we use λ within {−30, −29.9, . . . , 29.9, 30} and three ordered lists in L, leading to explore 1800 sets of K routes to find our best initial solution.

13

ACCEPTED MANUSCRIPT (a) grid layout

(b) example 117 85 11494

82

99

105 96

106

region of zones

118

104 87

119 90

91

93 100 113

95 103

101 115

108 107

11297 89

98 110 120 84 102

86

88

81 83 116 111

92 109

empty 20 36

region

75

11 4 23

28

48 18 2139 3 13

region of zones

62 41

14 10

65 42

77

52 73

7

5 25 1

37

24

1230 32 19

17

31 26

60

2 40 9

22 8

6

46

54 50

72

67

6445

7658

27 33 16

74 80 49

61

35

0

44 55

43 66 78

RI PT

region of zones

70 56 4757 69 79

15

63

51

6871

59

34 29 38

53

299 300 301 302

303 304 305 306

5.2. AVRP versus MCVRP

To highlight the potential benefits of our proposed approach, we compare our AVRP with the VRP maximizing the weighted sum of the number of covered zones and the sum of distances from zones not covered to the nearest zone on the route, as proposed by Current and Schilling (1994). We refer to the P P latter problem as the MCVRP. The objective function of MCVRP is defined as w2 yi + w4 ani , in

since we could therefore consider only indicator ati for zones not covered by routes. Our proposed problem yields slightly worst results for indicators considered in the MCVRP; see negative values in columns for P P indicators yi and ani . However, improvement in the number of opportunities and the accessibility

D

309

i∈I 0

311 312 313

i∈I 0

indicator based on the impedance function for zones not covered are quite large (i.e., up to 66.68%). With the AVRP, the magnitude of the average improvement for the number of zones not covered with access to opportunities is larger than the magnitude of the average decrement of the number of covered zones for both types of instances.

EP

310

i∈I 0

TE

308

i∈I 0

which parameters w2 and w4 are defined similarly as set for the AVRP. Table 2 shows the average improvement of accessibility indicators solved with the AVRP over the MCVRP for 60 small instances of types A1 and A2 (the ones we can solve to optimality with commercial P solver). We focus on summation ati (1 − yi ) to analyze the accessibility indicators based on impedance i∈I 0

307

M AN U

298

SC

Figure 4: (a) Grid layout for instance types A2, B2, C2, and D2; and (b) an example of instance of type D2 generated using that grid layout. 1

AC C

Table 2: Average improvement for sum of indicators of solving AVRP over Current and Schilling

Instance type A1 A2

314 315 316 317

318

P

i∈I 0

aai

4.42% 1.95%

Average improvement of AVRP over MCVRP P P P P yi ati (1 − yi ) ani aoi i∈I 0

i∈I 0

-3.84% -1.50%

66.68% 53.69%

320

i∈I 0

-0.49% -0.68%

64.13% 28.45%

as 11.06% 4.36%

For a clearer presentation of results in Table 2, we present Figures 5 which shows the comparison for each accessibility indicator among all small instances, where we exhibits the improvements described by the previous analysis. For example, it can be seen that for 75% of the instances of AVRP we obtain the P same value of coverage indicator yi and there are small losses for the rest of the instances while we i∈I 0 P obtain significant improvements for indicator ati (1 − yi ) for 41.66% of the instances by solving the i∈I 0

319

i∈I 0

AVRP instead of solving the MCVRP. Figure 6 shows a graph of the solutions obtained by the AVRP and MCVRP for instance 7 of type A1 14

ACCEPTED MANUSCRIPT

2

14

14

3

14

14

13

4

14

13

12

5

14

14

11

6

14

14

10

7

14

11

8

14

14

9

14

13

10 AVRP

14

MCVRP 14

1 11

9 14

9 12

2 12

14 13

14 13

3 13

14 9

13 9

4 14

11 14

12 14

15 5

14 9

13 9

16 6

12 13

11 13

9

17 7

14 10

14 11

8

18 8

11 14

11 13

19 9

13 10

13 11

14

9 8

14 13 12 11 10

7

X

8 2.1841 14

22 212

14 10

1

23 313

98

5

8 10 3.7442

24 414

14 2.1232 11

14 0.4483 11

25 515

14 9 3.8337

14 10 3.8337

26 616

12 8

1

9 112

27 717

14 3.4694 12

14 012

1

28 818

14 10

3

12 311

0

29 919

14 9

4

30 1020

14 AVRP 11

3

14 MCVRP311

311 1121

15 8785.6929

14 78 1.6929 9

322 1222

15 11 10

4

15 11 3.7651 11

333 1323

15 7564.6316

15 56 2.0602 7

344 1424

15 342.2237 12

15 36 1.9226 12

15 15 57 57 311 5 11 366 15 15 9 2.7414 9 1.1921 1626 10 10 377 15 15 43 101 1727 2 1.4 10 10 15 15 388 41 41 1828 1.7037 11 3.7037 12 15 15 399 49 41 2.4573 1929 12 3.2399 12 AVRP 15 MCVRP 14 40 10 38 2.65 36 2.65 2030 11 11 5 5 15 15 41 68 58 111 2131 2.0632 10 4.5971 11 1 1 15 15 42 43 39 122 2232 4 115 15 14 10 43 15 15 133 60 63 1 2333 1 12 12 3 1 44 15 15 144 33 33 211 2434 2 11 7 7 5 15 14 45 62 52 15 2535 2.565 12 2.565 12 1 1 15 15 46 115 109 166 1.2293 2636 13 1.2293 13 6 0 47 15 15 177 27 21 2.9549 2737 3.4286 12 12 4 4 48 15 15 188 39 33 2838 3 013 13 7 3 13 13 49 62 58 199 2939 2 AVRP 13 MCVRP213 10 6 4 14 14 50 37 37 20 311 3040 3 11 9 9 1 9 3 11 51 13 10 21 74 73 1.5052 3141 3.4205 13 13 2 ** ** 11 9 12 52 15 15 22 48 57 32423 0 013 13 4 4 13 6 4 15 15 53 138 138 23 3343 1.5416 161.5416 13 11 11 4 14 4 3 14 14 54 24 24 24 2.364 3444 2.364 9 9 5 13 13 6 4 15 55 15 15 25 371.6663 37 1.6663 3545 11** 12 6 ** 8 5 16 15 15 56 66 66 26 36467 2 214 14 26 9 17 4 2 15 14 57 48 48 27 3747 1.1367 261.1367 26 8 12 12 18 5 3 58 14 14 28 33 50 210 3848 2 7 9 13 13 7 5 19 59 15 15 29 290.8673 29 0.8673 3949 10 11 11 17 28 5 5 20 15 15 60 41 41 30 4050 2.7663 1.7083 21 14 11 11 11 21 9 4 57 52 31 1.0095 4151 231.5385 25 12 8 10 22 1 32 07 0 0.398400000000001 4252 130.398400000000001 18 12 13 13 2 2 23 33 392.3482 39 4353 2.3482 14 12 12 14 14 4 4 24 47 47 34 4454 1.0556 1.0556 16 16 15 12 12 25 5 5 33 33 35 1.1392 4555 3.6154 10 16 112 11 3 3 26 36 21 0.481 21 0.481 4656 17 46 34 12 12 321 5 4 27 37 432.1318 43 4718 2.1318 20 20 12 13 57 28 3 0 19 17 38 322 4819 1.2615 11 11 1.2615 58 11 11 29 2 2 18 18 39 1.1551 4920 341.1551 34 59 124 12 323 4 30 40 56 51 2.4231 5021 3 2.4231 11 13 13 60 9 4 31 27 25 41 5122 2.9173 012 324 22 61 32 0 0 25 25 42 1 0.7826 1 5223 0.7826 33 4 4 325 43 47 47 148 5324 1 48 7 7 34 44 261.1332 26 1.1332 5425 30 30 3 3 35 52 41 45 13 13 5526 3.2383 1.7299 36 5 5 10 10 46 1.107 5 1.107 5 5627 37 3 3 47 31 31 28 0.290699999999999 23 11 57 3 5 4 38 48 221.8157 22 56 56 5829 1.8157 3 3 39 55 55 49 37 37 5930 1.7795 1.4223 40 5 3 59 59 50 9 9 6031 2 2 3 3 41 51 106 170 32 ** ** 2 2 42 52 27 27 33 11 11

1

ati (1

14 11

14 9 1.1452

MCVRP

i2I 0

60

AVRP

yi )

MCVRP

6 5 4 3 2

X

1

AVRP

ani

i2I 0

180 135 90 45 0

X

i2I 0

1

60

aoi

AVRP

14 12 10 8 6 4 2 0

as

60

MCVRP

MCVRP

1

60

AVRP

MCVRP

60 45

AC C

Instancia

AVRP

15

20AVRP 10

12 110

60

yi

i2I 0

21 111

10 13 2.1841

1

RI PT

X

SC

14 MCVRP 11

MCVRP

15

355 1525

Instancia

AVRP

aai

i2I 0

M AN U

Instancia

X

13

D

Instancia

MCVRP 13

EP

Instancia

AVRP 1

TE

Instancia

30 15 0

1

60

Figure 5: Comparison of access indicators for solutions of instances A1 and A2 obtained by MCVRP and AVRP. Higher is P better except for indicator ani . i∈I 0

where green nodes, black nodes, pink nodes, and red nodes represent directly covered zones, indirectly covered zones, zones not covered with opportunities, and zones not covered without opportunities, respectively. Black arcs represent route segments, red dashed arcs represent opportunities, and gray dashed arcs represent missed opportunities. It can be seen that zones 3, 5, and 14 have no access to opportunities in the case of MCVRP (Figure 6-(a)) leading to geographical disaggregation, whereas all zones not covered 15

ACCEPTED MANUSCRIPT

326

had opportunities in the case of AVRP, even when there were more zones not covered (6-(b)). 5

(a) MCVRP

5

(b) AVRP 3

3

14

14

6 8

0

8

4

4 11 2

11 2 13

1 10 12

RI PT

6 0

13

1 10 12

7

7

9

9

328 329 330 331

Figure 7 shows the solutions obtained by the AVRP and MCVRP for instance 27 of type A2. where it can be seen that our proposed AVRP improves all accessibility indicators except for the number of covered zones. In general, the chief difference between our AVRP with the MCVRP is that we obtained 1 1 more, better-allocated opportunities for inhabitants in zones not covered by routes, as consistent with the concept of accessibility. (a) MCVRP

12

11 13 14

15

12

11 13 14

3

1

0 5

D

2

9

8 10

TE

4

(b) AVRP

15

3 0 5

M AN U

327

SC

Figure 6: Comparison of solutions obtained by MCVRP and AVRP for instance 7 of type A1.

2

1 9

8 10

4

7

7 6

EP

6

Figure 7: Comparison of solutions obtained by MCVRP and AVRP for instance 27 of type A2

5.3. Results of implementing our ILS on generated instances

AC C

332

334

In this section, we compared the CPLEX’s solver and our ILS for small instances and we analyzed 1 1 the behavior of the proposed ILS on large instances.

335

5.3.1. Comparison of CPLEX and ILS

333

336 337 338 339 340 341 342 343

To compare our ILS with the solver of CPLEX 12.6, we focused on computational time in seconds and gaps computed as the relative difference of the best dual bound obtained by CPLEX and the best feasible solution obtained by each method. Table 3 shows the average gap and time obtained by CPLEX and the ILS for instance types A1, A2, B1, and B2. Although CPLEX often outperforms heuristic algorithms for small instances in complex problems, since the average gap obtained by the ILS for instance types A1 and A2 was less than 0.8%, we obtained optimal or nearly optimal solutions for the smallest instances generated in the experimental stage in less than 0.5 s of average computational time. In particular, we found 21 and 18 optimal solutions for types

16

ACCEPTED MANUSCRIPT

Table 3: Comparison of CPLEX and ILS on small instances.

346 347 348 349 350 351 352

A1 and A2, respectively. For types B1 and B2, the difference between the average gap obtained by our ILS and CPLEX was 0.36% at most. In particular, we found three and six optimal solutions for instances of type B1 and type B2, respectively. Average computational times of our ILS were less than 1 s, whereas CPLEX reached the time limit for most instances of those types. We also compared CPLEX and ILS in 10 instances of type C1. Table 4 shows the gap and time (in seconds) for each instance. Feasible solutions obtained by CPLEX in 3 h of computational time were of significantly lower quality than those obtained with the ILS for all instances, for a difference of 10.40% in the average gap. That result demonstrates the usefulness of our ILS as a tool to obtain quality solutions for large instances of the AVRP in acceptable computational times.

SC

345

M AN U

344

gap 0.71% 0.40% 1.83% 1.55%

ILS time (secs) 0.39 0.44 0.61 0.69

RI PT

Instance type A1 A2 B1 B2

CPLEX gap time (secs) 0% 94.59 0% 130.20 1.47% 9703.46 1.38% 7597.24

Table 4: Comparison of CPLEX and ILS on ten instances of type C1.

353 354 355 356 357 358 359 360 361 362 363

AC C

EP

TE

D

Instance type C1 1 C1 2 C1 3 C1 4 C1 5 C1 6 C1 7 C1 8 C1 9 C1 10 Average

CPLEX gap 24.36% 20.07% 3.49% 17.94% 17.13% 7.31% 12.84% 19.31% 12.41% 15.42% 15.03%

gap 8.82% 4.27% 2.66% 4.46% 2.51% 2.86% 3.76% 4.05% 3.86% 4.61% 4.19%

ILS time (secs) 7.06 8.38 8.3 7.55 8.48 8.24 7.84 9.04 3.86 8.24 8.18

5.3.2. ILS on generated instances To showcase the behavior of our ILS on the generated instances, Table 5 presents the average computational time and average improvement of our ILS over the constructive algorithm among all instance types. Figure 8 shows a box plot of the improvements summarized by Table 5. In this Figure, black dots represent the average improvement, black line within the box shows the median, the box is delimited by the first and third quartiles, the dashed lines reach lower and upper limits, and the blank circles are atypical values. All types showed average improvements of at least 0.60%, whereas average computational times were less than 1 min. In particular, for types B1 and B2, the maximum route length allowed visiting almost all zones, even in the constructive phase, which yielded less improvement when the ILS was implemented.

17

ACCEPTED MANUSCRIPT

Table 5: Average improvement of our ILS over the initial solution obtained by Algorithm 1.

improvement of ILS over Constructor

time(secs)

A1 A2 B1 B2 C1 C2 D1 D2

2.28% 2.98% 0.60% 0.65% 3.57% 2.42% 2.54% 1.25%

0.38 0.44 0.60 0.69 8.09 8.89 50.68 46.80

RI PT

Instance type

20

●

● ●

M AN U

15

●

● ●

10

Improvement (%)

SC

●

●

●

5

●

●

●

●

●

TE

0

D

●

A2

● ●

●

●

B1

B2

C1

C2

D1

D2

Instance

EP

A1

●

Figure 8: Box plot for the improvements of ILS over Constructor() for all instances types of the AVRP.

365 366 367 368 369 370 371 372 373 374 375 376 377

Even if improvements shown in Figure 8 appear small for some instances, they may represent significant changes in values for the different indicators of our accessibility measure given the value of the weight parameters of the objective function. For example, Figure 9 shows the case of instance 14 of type C1, for which we obtained the greatest improvement of that instance type (12.06%). The improvement of our ILS over the constructor algorithm occurred primarily due to the increment in the number of zones not covered that got access to opportunities (zones 22, 29, 46, and 58) and the number of zones covered. Moreover, the average value of ati based on travel impedance among all zones not covered was 0.89 for the solution obtained by the ILS compared with the value of 0.46 for the solution obtained with the constructor algorithm. In short, there were more, closer opportunities for those zones. Figure 10 shows a better example of instance 26 of type C2 with an improvement of 7.09%. In that case, there were also more zones not covered with access to opportunities (i.e., zone 9), but the more significant improvement was due to 13 more covered zones in the solution obtained by the ILS. Moreover, there were not groups of zones not covered—that is, there was greater separation between those zones— which constitutes an important improvement in the indicator of spatial disaggregation (250%).

AC C

364

18

ACCEPTED MANUSCRIPT (a) Constructor 18

48 50

31

(b) ILS 44

33

59

18

48 50

0

31

59

41

11

44

33

0 41

11

22

22

13 12 29

25

34 27

24 2 15

54

56 13 12 29

19 9

20

28

42

37 14

46

57

6

30

8

39 5352

3 58

6

30

8

3 58

4043 49 38

1

51

23

39 5352

4043 49

7 21

55

35

17

17

7

21

SC

57

45 32

16

37 14 46

54 19

45 32

16

10

4

47

9

20

28

42

25

34 27

24 2 15

4

47

36

26

10

RI PT

36

26 56

5

55

35

38

1

51

23

5

M AN U

Figure 9: Instance 14 of type C1 which corresponds to the largest improvement on the objective function of ILS over Constructor() for that instances type. (a) Constructor 1

46

50

4859 60

42

45 47

46

53

54 55 58 57

44

43

49

51

D

9

27

10

18 2

7

6

19

17

15 4 3 12

45 47

13

33

24

21 26 34

54 55 58 57

44 51

43

49 52

9

27

39

37 16 18

11

25 40

10 35

18 2

32 22 7

31

2329 28

39

14

20

38 0

30

36

EP

11

5

2329

28

TE

37

14

53

56

52

18

50

41

56

16

1

4859 60

42

41

20

(b) ILS

5

17

38 0

30

19

25 40

36

35 6

15 4 3 12

33 13

24 21 26 34

32 22 31

378 379 380 381 382 383 384 385 386 387 388

AC C

Figure 10: Instance 26 of type C2 which corresponds to the largest improvement on the objective function of ILS over Constructor() for that instances type.

1 of D1 and D2, Figure 11 shows the case of instance 1 In case of larger instances 28 of type D1, which demonstrated an improvement of 6.06%. Similarly to cases shown in Figures 9 and 10, the improvements emerged due to more zones covered and not covered with opportunities for access. That instance represents a scenario in which the maximum route length is not enough to avoid spatial disaggregation, since a large concentration of zones not covered appears in the low area of Figure 11, which could represent the limit of the distribution process. However, we increased the number of opportunities and the value of accessibility indicator based on the impedance function for zones not covered in that delimited area. Finally, we analyzed instance 24 of type D2 (Figure 12), which represented the largest improvement (i.e., 5.07%) for that instance type. In that case, 5.07% consisted of an additional zone not covered with access to opportunities and additional covered zones, among other changes, for the rest of the accessibility indicators.

19

ACCEPTED MANUSCRIPT (a) Constructor

35

55 27

77 57

14 85

65 11

108 45

11892

83 87

32

998

101

21 13 53

84 70

61 6344

83

87

54 5 98

73

109

111 68 998

24

72 1

31

46

107

56

88

113

84 70

59

54 5 98

33 112 4 74 9 26 100 95 105 110 75 114

61 6344

73

109

28

72

1

30

67

94

39

37 34

80

101

21 13 53

79

81

71

115

3

28

15

58

23

82

19 29 2

6

59

30

67

11892

93

47

90

104

102

20 106 7

32

33 112 4 74 9 26 100 95 105 110 75 114

68

108 45

56

113

55

10322 48 116 36 41 117 49

16

64 38 91 76 25 18

17

42

35

14 85

65 11

88

80

111

43 66

10

27

77 57

107

37 34

115

3

79

39

52

60

71

58

23

81

89

19 29 2

6

24

15

102

20 106 7

82

40

1278

119 062 96 50 69 51 86

97

93

47

90

104

10322 48 116 36 41 117 49

16

64 38 91 76 25 18

17

42

60

43 66

10

RI PT

52

SC

89

40

1278

119 062 96 50 69 51 86

97

(b) ILS

31

46

94

M AN U

Figure 11: Instance 28 of type D1 which corresponds to the largest improvement on the objective function of ILS over Constructor() for that instances type. (a) Constructor 1 105

111

92

109 104

112 83 118 82 95 98 115 114 103

87 96 91

93 99

102 107

18 36 219 10 26

3 12 6 16 11 15 23 30

4 13

31

3314

35

1

39 21 34 29 38 227 5

17 40 37 8 9 25

113 101 86

115 114 103

89 116 84108 117

102 107

45 71

42

70

5749 59 78 65

EP

24

94

0

43 64

77 68 44 58

109 104

112 83 118 82 95 98

D

27

28

88 110 85

1 111

92

106 9781 119 120

56

52 4867 41 66 74

75

20

54 46 69

TE

20 32

100

105

90

(b) ILS

27

28

32 18 36 219 10 26

50 80

55 73 79 62 76 53 6063 47 51 61

87 96 91

93 99

100

88 110 85

94

90 106 9781 119 120 113 101 86 89 116 84108 117

3 12 6 16 11 15 23 30

4 13

45 71 42

5749 59 78

31

3314 24 35

72

1

39 21 34 29 38 227 5

70

65

17 40

77 68 44 58

37 8 9 25 0

43 64

54 46 69

56 52 4867 41 66 74 75

50 80

55 73 79 62 76 53 6063 47 51 61

72

AC C

Figure 12: Instance 24 of type D2 which corresponds to the largest improvement on the objective function of ILS over Constructor() for that instances type.

391

Based on the numerical 1results in this section, our iterated local search1 is a tool that can be implemented on large scenarios to improve the accessibility of a distribution processes, in order to guarantee more opportunities requiring less mobility for inhabitants of zones not covered.

392

6. Conclusions and further research

389 390

393 394 395 396 397

Improving accessibility consists of guaranteeing more opportunities for people to obtain services or have demands met without having to rely as much on mobility. Despite the different definitions of accessibility for different contexts, common elements of measures of accessibility include the number of opportunities to obtain the service, the specific cost of traveling, and geographical disaggregation. In our study, we defined an accessibility measure as a weighted sum of indicators for those elements. 20

ACCEPTED MANUSCRIPT

423

Acknowledgments

403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421

SC

402

M AN U

401

D

400

TE

399

RI PT

422

We optimized this measure in our proposed new variant of the SVRP, which we called as AVRP, that considered a radius of coverage for zones along routes and a radius of mobility for zones not covered by any route. We find out that classic coverage metrics for VRPs do not consider all the elements of accessibility defined in our study. To exhibit the latter, we compared a maximal covering vehicle routing problem— focused only on the number of covered zones and the distance to the nearest opportunity—with our AVRP where we obtained more and better-allocated opportunities to fulfill a demand or receive a service and we diminish geographical disaggregation, which is more desirable according to an inclusive philosophy consistent with the definition of accessibility. In particular, the solutions obtained with our approach yielded improvements between 1.95% and 66.68% for indicators included only in the AVRP, while we obtain losses less than 3.84% for indicators optimized in the MCVRP. Our proposed mathematical formulation is intractable by implementing commercial solver, even for small instances. In turn, we designed an ILS capable of obtaining good-quality solutions for large instances of the AVRP in less than 1 min of computational time. As a result, we generated nearly optimal solutions—at most, of 1.83% of the average gap—for small instances types. Furthermore, by implementing our proposed local searches, we also significantly improved elements of accessibility in the objective function over initial solutions generated by our constructive algorithm. A further research area is including the demand, as well as capacities for zones, in order to determine an assignment of opportunities for inhabitants in zones not covered by routes. Moreover, the proposed measure of accessibility can be adapted for facility location problems in cases in which demand is critical for people (e.g., distribution of bare necessities and humanitarian relief) and a delimited mobility for inhabitants is considered. Besides, a comparison of different accessibility measures is of interest as well as the practical implementation of such approaches in real scenarios. Finally, considering stochastic parameters for the problem (for example, Allahviranloo et al., 2014) is also of interest since a potential implementation is in the context of humanitarian relief where uncertainty is present in the input data.

398

426

Appendix: Notation list

427

Sets

429

430

AC C

428

EP

425

This work was supported by the Programa para el Desarrollo Profesional Docente (PRODEP), Grants 511-6/17-7538 and DSA/103.5/16/14538.

424

I: set of zones.

I1 : set of directly covered zones. I2 : set of zones not covered by routes.

431

K: number of homogeneous vehicles.

432

N (i): set of zones covered indirectly by i if that zone is visited by any route.

433

A(i): set of accessible zones for zone i.

21

ACCEPTED MANUSCRIPT

434

Parameters tij : travel time from zone i to zone j.

436

T : time limit of the routes.

437

f (tij ) impedance function for traveling from zone i to zone j.

438

wk : weight parameter.

439

s: solution for the AVRP.

440

sij : savings parameter for pair of zones (i, j).

441

λ: route shape parameter.

442

L: set of ordered lists of zones.

443

L: ordered list of zones which denotes an element of L.

444

rN : coverage radius to define the set N (i).

445

rA : accessibility radius to define the set A(i).

446

LS: List of pairs of zones (i, j) ordered accordingly savings parameters sij .

451

452 453 454 455

456 457 458

459 460

461 462

463 464

465 466

SC

M AN U

D

450

Altay, N., Green III, W., 2006. OR/MS research in disaster operations management. European Journal of Operational Research 175, 475 – 493.

TE

449

Allahviranloo, M., Chow, J.Y., Recker, W.W., 2014. Selective vehicle routing problems under uncertainty without recourse. Transportation Research Part E: Logistics and Transportation Review 62, 68 – 88.

Aras, N., Aksen, D., Tekin, M.T., 2011. Selective multi-depot vehicle routing problem with pricing. Transportation Research Part C: Emerging Technologies 19, 866 – 884. Freight Transportation and Logistics (selected papers from ODYSSEUS 2009 - the 4th International Workshop on Freight Transportation and Logistics).

EP

448

References

Bertolini, L., le Clercq, F., Kapoen, L., 2005. Sustainable accessibility: a conceptual framework to integrate transport and land use plan-making. two test-applications in the netherlands and a reflection on the way forward. Transport Policy 12, 207 – 220.

AC C

447

RI PT

435

Clarke, G., Wright, J.W., 1964. Scheduling of vehicles from a central depot to a number of delivery points. Operations Research 12, 568–581. Current, J.R., Schilling, D.A., 1994. The median tour and maximal covering tour problems: Formulations and heuristics. European Journal of Operational Research 73, 114 – 126. Farahani, R., Asgari, N., Heidari, N., Hosseininia, M., Goh, M., 2012. Covering problems in facility location: A review. Computers and Industrial Engineering 62, 368–407. Flores-Garza, D.A., Salazar-Aguilar, M.A., Ngueveu, S.U., Laporte, G., 2015. The multi-vehicle cumulative covering tour problem. Annals of Operations Research , 1–20.

22

ACCEPTED MANUSCRIPT

474

475 476

477 478

479 480 481

482 483

484 485 486

487 488 489

490 491

492 493

494 495

496 497

498 499

500 501

502 503

RI PT

473

Handy, S., Niemeier, D., 1997. Measuring accessibility: An exploration of issues and alternatives. Environment and Planning A 29, 1175–1194. Juwana, I., Muttil, N., Perera, B., 2012. Indicator-based water sustainability assessment - a review. Science of The Total Environment 438, 357 – 371. Liu, H., Li, F., Xu, L., Han, B., 2015. The impact of socio-demographic, environmental, and individual factors on urban park visitation in beijing, china. Journal of Cleaner Production in press.

SC

472

Louren¸co, H.R., Martin, O.C., St¨ utzle, T., 2010. Iterated Local Search: Framework and Applications. Springer US, Boston, MA. pp. 363–397. Makr´ı, M., Folkesson, C., 1999. Accessibility measures for analyses of land use and travelling with geographical information systems, in: Proceedings of the 2nd KFB-Research Conference, Lund, Sweden. pp. 1–17.

M AN U

471

Geurs, K., van Wee, B., 2004. Accessibility evaluation of land-use and transport strategies: Review and research directions. Journal of Transport Geography 12, 127–140.

Nolz, P., Doerner, K., Gutjahr, W., Hartl, R., 2010. A Bi-objective Metaheuristic for Disaster Relief Operation Planning. Springer Berlin Heidelberg, Berlin, Heidelberg. pp. 167–187. Palomo-Mart´ınez, P., Salazar-Aguilar, M., Laporte, G., Langevin, A., 2017. A hybrid variable neighborhood search for the orienteering problem with mandatory visits and exclusionary constraints. Computers and Operations Research 78, 408–419.

D

470

S´anchez-Garc´ıa, S., Athanassiadis, D., Mart´ınez-Alonso, C., Tolosana, E., Majada, J., Canga, E., 2017. A gis methodology for optimal location of a wood-fired power plant: Quantification of available woodfuel, supply chain costs and ghg emissions. Journal of Cleaner Production 157, 201–212.

TE

469

Tang, H., Miller-Hooks, E., 2005. A tabu search heuristic for the team orienteering problem. Computers and Operations Research 32, 1379 – 1407.

EP

468

Geem, Z., Tseng, C., Park, Y., 2005. Harmony Search for Generalized Orienteering Problem: Best Touring in China. Springer Berlin Heidelberg, Berlin, Heidelberg. pp. 741–750.

de la Torre, L., Dolinskaya, I., Smilowitz, K., 2012. Disaster relief routing: Integrating research and practice. Socio-Economic Planning Sciences 46, 88 – 97.

AC C

467

Tricoire, F., Romauch, M., Doerner, K., Hartl, R., 2010. Heuristics for the multi-period orienteering problem with multiple time windows. Computers and Operations Research 37, 351 – 367. Vansteenwegen, P., Souffriau, W., Berghe, G., Oudheusden, D.V., 2009. Iterated local search for the team orienteering problem with time windows. Computers and Operations Research 36, 3281 – 3290. Vansteenwegen, P., Souffriau, W., Oudheusdena, D.V., 2011. The orienteering problem: a survey. European Journal of Operational Research 209, 1–10. Wang, Q., Sun, X., Golden, B., Jia, J., 1995. Using artificial neural networks to solve the orienteering problem. Annals of Operations Research 61, 111–120. Wang, X., Golden, B., Wasil, E., 2008. Using a Genetic Algorithm to Solve the Generalized Orienteering Problem. Springer US, Boston, MA. pp. 263–274.

23

ACCEPTED MANUSCRIPT

M AN U

SC

RI PT

Yellow, P., 1970. A computational modification to the savings method of vehicle scheduling. Operational Research Quarterly (1970-1977) 21, 281–283.

D

507

TE

506

EP

505

Wang, Y., Ma, X., Li, Z., Liu, Y., Xu, M., Wang, Y., 2017. Profit distribution in collaborative multiple centers vehicle routing problem. Journal of Cleaner Production 144, 203 – 219.

AC C

504

24