A heuristic algorithm for a single vehicle static bike sharing rebalancing problem

The static bike rebalancing problem (SBRP) concerns the task of repositioning bikes among stations in self-service bike-sharing systems. This problem can be seen as a variant of the one-commodity pickup and delivery vehicle routing problem, where multiple visits are allowed to be performed at each station, i.e., the demand of a station is allowed to be split. Moreover, a vehicle may temporarily drop its load at a station, leaving it in excess or, alternatively, collect more bikes from a station (even all of them), thus leaving it in default. Both cases require further visits in order to meet the actual demands of such station. This paper deals with a particular case of the SBRP, in which only a single vehicle is available and the objective is to find a least-cost route that meets the demand of all stations and does not violate the minimum (zero) and maximum (vehicle capacity) load limits along the tour. Therefore, the number of bikes to be collected or delivered at each station should be appropriately determined in order to respect such constraints. We propose an iterated local search (ILS) based heuristic to solve the problem. The ILS algorithm was tested on 980 benchmark instances from the literature and the results obtained are quite competitive when compared to other existing methods. Moreover, our heuristic was capable of finding most of the known optimal solutions and also of improving the results on a number of open instances.


Introduction
The task of repositioning a commodity from one location to another is a well-known problem arising in different contexts such as logistics, transportation, and various disciplines, notably industrial engineering and operations management.A practical application arises in self-service bike sharing systems (BSS), which are becoming increasingly popular in recent years.Users rent bikes and return them at stations distributed over a region.A vehicle with limited load capacity then periodically collects and delivers bikes across different stations so as to rebalance the system.
Alternatives to the street traffic are important not only because of its impact in urban congestion, but also in the environment, commuting, and so on.The emerging worldwide BSS are proving to be an effective solution to mitigate the effects of traffic issues in large urban centers.Up to 2009, there were about 120 bike sharing programs around the world and, according to DeMaio (2009), they have a favorable impact on: decreasing traffic congestion, improving public health, and helping reducing the level of CO 2 emissions.One of the most famous systems is the Vélib' system in Paris, with 1800 stations and more than 20,000 bikes.
In such systems, each station has an inventory with a load capacity, an initial number of bikes, and consequently a number of free slots where users can return bikes to the system.Throughout the day, some stations may have no bike to be rented or free slots to store returned bikes.Therefore, an attempt to avoid this scenario, which is unpleasant for users, is to determine an initial acceptable number of bikes (and free slots) at each station.This task can be done based on demand history and peaks at each station.
The activity of repositioning bikes among stations on a regular basis is called rebalancing, and this is done by one or more vehicles that move bikes from one station to another in order to restore its inventory to the initial desired configuration.As per DeMaio (2009), good rebalancing systems are present in successful bike sharing programs, and since the vehicles move back and forth across an urban area, a vehicle routing optimization can be utilized.
The rebalancing is either static, performed when nearly no bikes are being used, or dynamic, which is done while the system is still in use.The static bike rebalancing problem (SBRP) is motivated by the fact that very few bikes are being used at night.Indeed, according to Chemla et al. (2013a) and Dell'Amico et al. (2014), there are many systems that are even closed during the night.
In this work we consider the single vehicle SBRP, which is clearly N P-hard, because it includes, among others, the classical traveling salesman problem (TSP) as a special case.The SBRP can be seen as a variant of the one-commodity pickup and delivery TSP (Hernández-Pérez and Salazar-González, 2004a,b), where multiple visits are allowed to be performed at each station, i.e., the demand of a station is allowed to be split.Moreover, a vehicle may arbitrarily drop its load at a station, leaving it in excess or, alternatively, collect more bikes from a station (even all of them), thus leaving it in default.Both cases require further visits in order to meet the actual demands of such station.This strategy of allowing a station to act as a buffer or a temporary depot is denoted as preemption or temporary operation (i.e., temporary pickup and temporary dropoff).Finally, visits to balanced stations are optional for the SBRP.
Salazar-González and Santos-Hernández (2015) considered a similar yet different problem, in which an upper limit is imposed on the number of visits to the customers and to the depot, and the single vehicle that performs the rebalancing is not forced to leave the depot with an empty load.
The works of Chemla et al. (2013a) and Erdogan et al. (2015) were the only ones to consider the same variant dealt in the present paper.Chemla et al. (2013a) proposed a mathematical formulation over an extended graph, where each station is replicated according to an upper bound on the number of visits.Due to its visible intractability, two relaxations were developed.The authors also presented among other contributions, a polynomial algorithm to compute optimal bike displacements for a given sequence of vertices, which is useful to determine if a route is feasible or not, as well as tabu search heuristics and a branch-and-cut algorithm that solves a relaxation of the problem.Erdogan et al. (2015) proposed

Problem description
Let n be the number of stations, V = {1, ..., n} be the set of vertices representing their locations (station 0 represents the depot), and A be the set of arcs in a complete and directed graph G = (V ∪{0}, A).For each arc a (i,j) ∈ A, there is a cost c a , satisfying the triangular inequality (c (i,j) + c (j,k) ≥ c (i,k) , ∀i, j, k ∈ V ).
For each i ∈ V , let p i ∈ Z be the amount of bikes initially stored, p i ∈ Z be the number of bikes requested by i after the service is performed, and d i = p i − p i be the total demand.When d i > 0 and d i < 0, we assume that i ∈ V is a delivery and a pickup station, respectively.A station i ∈ V may have no demand (p i = p i ) and in this case the visit becomes optional.Each station i ∈ V has a capacity q i ∈ Z and the depot is assumed to have no bikes, i.e., q 0 = p 0 = p 0 = 0. Finally, let Q ∈ Z be the vehicle capacity.
The objective is to find a least-cost route that starts and ends at the depot, visits each station with non-zero demand at least once, meets the demands of all stations (i.e., the initial load p i is changed to the target demand p i , ∀i ∈ V ), and does not violate the minimum (zero) and maximum (Q) load limits.Therefore, the number of bikes to be collected or delivered at each visit to a station should be appropriately determined in order to respect such constraints.
Finally, stations may serve to perform temporary operations (preemption), either as a temporary depot or a temporary buffer, i.e., supply more bikes than their initial demand or hold more bikes (without exceeding its inventory load capacity), and in both cases have their demand satisfied in later visits.
Figure 1 shows a graphical representation of an optimal solution for the benchmark instance n20q10D (n = 20 and Q = 10).The nodes are distributed according to the spatial coordinates of the stations.
The positive and negative values next to the nodes are the number of bikes collected and delivered, respectively.The arcs and their associated values represent the vehicle traveling to the next station in the sequence and the vehicle load, respectively.For example, the vehicle delivers 2 bikes in the first visit to station 12, collects 10 at station 10, returns to 12 to deliver 6 more (meeting the demand of 8) and then travels to station 14 with a load of 4 bikes.3 Proposed heuristic ILS iteratively alternates between local search (intensification) and perturbation (diversification) mechanisms with a view of finding high quality solutions.In our case, we embed a variable neighborhood descent (VND) (Mladenović and Hansen, 1997) based procedure in the local search phase of the metaheuristic.As in previous works (e.g., Penna et al., 2013, andSilva et al., 2015), the neighborhoods of our algorithm are examined in a random ordering during the search (RVND).
As opposed to most of the former ILS implementations cited in Section 1, infeasible solutions are temporary accepted after the application of perturbation moves, not only for the sake of diversification, but also as an attempt to escape from local optimal solutions.This modification was crucial for the favorable performance of the heuristic when dealing with the single vehicle SBRP considered here, which appears to be more challenging to solve than other VRPs where ILS was successfully applied to obtain high quality solutions by only considering the feasible search space.
The proposed hybrid heuristic, called ILS SBRP , combines multiple restarts, local search, perturbations mechanisms, and a repair phase.Figure 2 illustrates the flowchart of ILS SBRP .For each of the I R restarts, a feasible initial solution is generated using a simple greedy randomized constructive algorithm (see Section 3.3).Next, local search, perturbation and repair procedures are successively applied until the stopping criterion is met, that is, when the number of consecutive attempts to escape from a local optimal solution reaches I ILS trials.Because perturbation moves are allowed to produce infeasible solutions, we implemented a procedure called AddUnbalanced (Chemla et al., 2013a) (see Section 3.2) with the aim of repairing such solutions.Nevertheless, there is no guarantee that a solution will be feasible after applying this procedure.When an infeasible solution is not totally repaired and the local search does not find a move that leads to a feasible solution, then the infeasible solution is disregarded.Note that perturbation is always applied over the best solution of the current multi-start iteration.Finally, ILS SBRP returns the best solution found among all restarts.

Solution representation
A solution for the single vehicle SBRP considered in the present work can be represented as a sequence of visits to stations, starting and ending at the depot, along with the amount of bikes collected or delivered at each visit.
Three vectors are used as data structures to store: (i) the route, where the first and last element are fixed at 0, i.e., the depot; (ii) the operation performed by the vehicle during a visit, where negative and positive values indicate the amount of bikes delivered and collected, respectively; and (iii) the vehicle load during the route.
As in Chemla et al. (2013a), a flow network is used to check in polynomial time whether or not a solution is feasible, with respect to bike displacements and vehicle capacity, given a sequence of vertices representing visits to stations.A detailed explanation can be found in Appendix A.
We also use another data structure which consists of a key-value map composed by n + 1 elements that store the number of visits performed at each station.This is useful, for example, to check whether a solution includes all stations with non-zero demand.Note that information hold in (ii) is extracted from the computed bike displacements when solving the max-flow problem (see Appendix A).From such, it is possible to derive, in linear time, the vehicle loads in (iii) by the adding or subtracting the bike displacements at each visit.

Repairing infeasible solutions
As already mentioned, infeasible solutions are allowed after perturbations.We therefore implemented a procedure called AddUnbalanced (Chemla et al., 2013a), which tries to repair a solution by adding stations to the route.More precisely, both the most unbalanced station in excess and in default, i and j, respectively, are selected and three moves are possible: (i) adding arcs (j, i) and (i, j) after the existing visit to j; (ii) adding arcs (i, j) and (j, i) after the existing visit to i; (iii) if both i and j are not in the sequence, adding (i, j) at the end of the sequence, before returning to the depot.
For example, let us consider a scenario where stations i = 12 and j = 14 are the most unbalanced.
More precisely i has initially 20 bikes and a demand of −10, i.e., a pickup station, while j is initially holding 3 out of 10 (target) bikes, i.e., a delivery station with demand 7.An infeasible solution is presented in Figure 3a, where the referred stations are not balanced, that is, their demands are not met, since only 4 bikes were collected in station 12 and 4 bikes were delivered at station 14. Figure 3b shows a modified solution, where after the addition of arcs (14, 12) and (12, 14), a new and feasible configuration of bike displacements were determined by means of the maximum flow check (see Appendix A).We can observe that the second visit complements the first one, meeting the demand of both stations: the vehicle deliveries 1 bike at station 14, collects the remaining 7 at station 12, now balanced, and finally meets the demand of station 14 by delivering 6 more bikes.
It is worth emphasizing that the AddUnbalanced procedure does not necessary lead to a feasible solution.However, in general, experimental results showed that such procedure has a high level of success in fully repairing infeasible solutions.

Constructive Procedure
The pseudocode of the greedy randomized constructive procedure is presented in Alg. 1.The algorithm stores and maintains a list of open vertices (OV ) corresponding to stations whose demands are still not fully met.Stations without demand are also included in this list.In order to ensure a level of diversity during the process of generating an initial solution, OV is sorted in random order (line 4).
The algorithm follows a greedy procedure by selecting the first vertex to be inserted at the end of However, it may come to a point where no station can be fully served in a single visit, either because the vehicle has not enough bikes to deliver, or the residual capacity is not sufficient to collect the required bikes at once.Hence, a split becomes necessary.The second part of the constructive procedure (lines 13-17) iterates over OV searching for a station whose demand maximizes the number of bikes that can be delivered or collected.Ties are broken according to the nearest insertion criterion.The station demand and vehicle load are updated after the insertion.Next, the algorithm restarts from line 5 and the entire insertion procedure is repeated until OV becomes empty.Note that the generated initial solution is always feasible.

Local search
Initial and perturbed solutions are possibly improved by means of an RVND based procedure during the local search.RVND consists of systematically examining different types of neighborhoods in a random ordering.In particular, if the best neighbor consists of an improving move, then the search may continue from any of the existing neighborhoods (including the last one to be explored) at random.Otherwise, a different neighborhood other than those that did not succeed in finding an improved move is randomly selected.The procedure ends when all neighborhoods fail to refine the current solution.Only feasible moves are accepted.
The following six neighborhood structures were implemented.

Working Paper
Algorithm 1 Initial Solution Constructive Procedure OV ← List sorted in random order with all stations where d i = 0 + random ones with d i = 0 5: repeat 6: inserted ← f alse 7: for all i ∈ OV do 8: Solution ← Solution ∪ i update Q 20: until OV = ∅ 21: return Solution 22: end GenerateInitialSolution.
• Reinsertion -N (1) : A station is removed and then reinserted in another position of the sequence.
• Or-opt2 -N (2) : Two consecutive stations are removed and then inserted in another position.
• Or-opt3 -N (3) : Three consecutive stations are removed and then inserted in another position.
• 2-opt -N (4) : Two non-adjacent arcs are removed from the sequence and then two new ones are inserted.In other words, a subsequence of the tour is reversed.
• Suppression -N (6) : Given a sequence L = i 0 , i 1 , ..., i k , a suppression list is composed of visits to stations i j , ∀j ∈ {1, . . ., k − 1}, such that p ij = p ij (zero demand) or p ij = p ij and i j is visited more than once in the tour.The best move, if any, consists in selecting one station to be removed from L so that the solution cost is minimized and the resulting new sequence L is feasible.For example, Figure 4b shows the removal of an additional visit to station 2, thus modifying the subsequence 2, 6, 2, 9, 0 to 2, 6, 9, 0. This neighborhood was originally proposed by Chemla et al. (2013a), but the authors considered all stations.
The first five are well-known TSP neighborhood structures, while the last is a problem-specific neighborhood.Figure 4a depicts an initial solution and Figures 4b to 4g illustrate modified solutions that were obtained after changing the previous one by means of one of the neighborhoods described above.
For example, Figure 4d shows a solution in which a 2-opt move was applied over the solution shown in Figure 4c.For ease of presentation, values of pickup/delivery operations as well as the vehicle load are omitted.

Perturbation mechanisms
One of the four mechanisms described below is selected at random whenever the algorithm enters the perturbation phase.• AddBuffer -P (1) : An additional visit to a station is included, expecting to act as buffer, using the cheapest insertion criterion.Unrouted stations are inserted twice using the same criterion (Chemla et al., 2013a).
• AddStations -P (2) : This perturbation mechanism generalizes the previous one in the sense of allowing multiple visits to be added in the solution, but with a different insertion criterion.More precisely, an additional visit (or two, in the case of unrouted stations) to up to three random stations are included towards the end of the route.Here we only consider stations that are visited at most once.Adjacent visits to the same station are forbidden.
• Double-Bridge -P (3) : Introduced by Martin et al. (1991) for the TSP, this perturbation consists of a permutation of two subsequences.As a result, four arcs are removed and four new ones are added so as to generate a new sequence.
• Suppression -P (4) : A suppression move (see Section 3.4) is applied at random, but in this case the resulting modified sequence is allowed to be infeasible.Figure 5 shows an example of perturbations applied over a (supposedly) local optimal solution.In Figure 5d, a Double-Bridge move is applied by interchanging subsequence 6, 4, 1 with subsequence 7, 9.
Figure 5b shows the AddBuffer perturbation, when an additional visit to station 7 is performed expecting it to act as a buffer.In Figure 5c, the perturbation AddStations is applied by adding two random visits: one to station 8 and another one to station 6.

Computational experiments
The ILS SBRP algorithm was coded in C++ (g++ 4.6.4)and the computational tests were carried on an Intel R Core TM i7-3770 with 3.40 GHz and 16 GB of RAM running Ubuntu 14.04.Only a a single thread was used during the experiments.

Instances
The benchmark instances used to test the proposed algorithm are those suggested by Hernández-Pérez and Salazar-González (2004a), which were originally created for the one-commodity pickup and delivery traveling salesman problem.The benchmark contains instances ranging from 20 to 500 customers (stations), and vehicle capacities ranging from 10 to 1000.For each pair of problem size and vehicle capacity, there are 10 instances named from A to J and, for each vertex i, there is a demand d i ∈ [−10, 10].Chemla et al. (2013a) and Erdogan et al. (2015) only reported results for a subset of instances of the referred benchmark.Therefore, in order to compare our results with theirs, we tested ILS SBRP for all instances considered in at least one of the two works (see Section 4.4).Furthermore, to compute the initial and final targets as well as the load capacity for each station, the same procedures adopted by such authors were employed: for each vertex i, p i = α × 10, p i = α × (10 + d i ), q i = α × 20, where α is an input parameter, and experiments were conducted with α = 1 and α = 3.In order to properly compare our results with those in Chemla et al. (2013a) and Erdogan et al. (2015), we adopted their same convention of rounding down the values of the cost matrix to the nearest integer (floor), although we noticed that this can cause slight violations of the triangle inequality.

Impact of the perturbation mechanisms
In this section we are interested in evaluating the impact of the perturbation mechanisms described in Section 3.5, that is, AddBuffer (P (1) ), AddStations (P (2) ), Double-Bridge (P (3) ), and Suppression (P (4) ).
In view of this, we selected a subset of 30 challenging instances for performing the experiments.These instances were chosen according to the largest gap values with respect to the lower bounds reported in Erdogan et al. (2015).We ran ILS SBRP 10 times for each of the 30 instances considering all possible combinations of perturbations.For this testing we arbitrarily adopted I R = 10 and I ILS = n.

Parameter tuning
The main ILS SBRP parameters to be calibrated are the number of restarts (I R ) and the maximum number of consecutive ILS iterations without improvement over the current local optimal solution (I ILS ).Here we set I R = 10, as in Silva et al. (2015), where the authors put forward a multi-start ILS that was capable of obtaining state-of-the-art results for the split-delivery VRP.
In previous works, such as those mentioned in Section 1, the value of I ILS was tuned based on the instance size.In VRPs (e.g., Penna et al. (2013) I min is more important for small size instances and its role is to prevent low values for I ILS , which in this case may lead to an insufficient number of ILS iterations required for obtaining high quality (or even optimal) solutions.We then tested several values for I min , more specifically, 100, 120, 140, 160, and 180.
For each of them, we ran ILS SBRP 10 times for all instances containing 20 and 30 stations.The average results obtained suggest that I min = 160 seems to provide a good compromise between solution quality and CPU time, since the algorithm managed to find almost all best known solutions in a relatively short amount of time when using this value.Therefore, we set I ILS = max{160, 4 × n}., 15, 20, 25, 30, 35, 40, 45, 1000}. Chemla et al. (2013a ran their experiments on an AMD Athlon 5600+ 2.8 GHz with 16 GB of RAM, while Erdogan et al. (2015) performed their testing on an Intel i7 3.60 GHz and 8 GB of RAM.On the one hand, because the hardware performance of the first appears to be quite inferior to the second, as well as to our intel i7 3.40 GHz, we decided to estimate an approximation factor based on the single thread rating values reported in https://www.cpubenchmark.net/compare.php?cmp[]=86&cmp[]=896, so as to better compare the runtime performance of the methods.According to the referred website, the AMD Athlon 5600+ 2.8 GHz is roughly 2.43 times slower than our processor.We thus report the original CPU time values of Chemla et al. (2013a) divided by a factor of 2.43.On the other hand, since the machine used in Erdogan et al. (2015) is rather equivalent to ours, perhaps even slightly faster, we decided to consider the original runtime values reported by the authors.

Results for instances with up to 60 stations
The aggregate average results for instances containing 20, 30, 40, 50, and 60 stations are reported in Tables 2 and 3, where Instance group denotes the set of 10 instances of a particular group (for example, group n20q10 contains 10 instances with n = 20 and Q = 10); UB1 Gap (%), UB2 Gap (%), and Gap (%) correspond to the gap between UB1, UB2, and the upper bound found by Erdogan et al. (2015), respectively, and the lower bound reported in From Table 2, it can be observed that the quality of the solutions found by ILS SBRP , as well as those obtained by the algorithm of Erdogan et al. (2015), are visibly superior than the ones determined by the tabu searches of Chemla et al. (2013a), especially TS1.Such superiority becomes even more prominent for α = 3, as presented in Table 3.Also, the average CPU times spent by ILS SBRP to find or improve the best solutions reported by Chemla et al. (2013a) (UB2) are rather small in most cases, sometimes only a matter of relatively few seconds, except for the instance group n40q10 when α = 3, where the proposed algorithm required more CPU time.
In addition, assuming the same values for the demands, the smaller the vehicle capacity, the larger the relative number of visits.This increases the size of the tour, thus affecting the number of operations performed during the local search, and possibly the number of ILS iterations, as more moves are required to be evaluated.However, this was somewhat expected since the CPU effort of the algorithms of Erdogan et al. (2015) and of Chemla et al. (2013a) tend to increase exponentially with increasing values of n.
Finally, Table 4 shows a summary of the best solutions found by the proposed algorithm, where #Opt  the proposed heuristic was capable of significantly improving the best known solution of all instances.
In addition, when considering α = 1, ILS SBRP required, on average, at most 8 seconds to find or improve the best results of Chemla et al. (2013a).In some cases, such as those involving Q ≥ 30, our algorithm spent, on average, only a fraction of a second to achieve a superior solution than the best one from the literature.For α = 3, more time was required, on average, to accomplish the same purpose, but mostly for Q = 10.

Concluding Remarks
In this work we proposed a hybrid ILS algorithm that was especially designed to solve a challenging single-vehicle SBRP variant.Extensive computational experiments were conducted on 980 instances from the literature ranging from 20 to 100 stations.The results were compared with those reported in Chemla et al. (2013a) and Erdogan et al. (2015).For the 900 instances containing up to 60 stations, the proposed heuristic, called ILS SBRP , was capable of finding 796 out of 823 known optimal solutions (97%) and improving the result of 61 out of 77 open instances (79%).Our algorithm only failed to be at least equal to the best known solution in 27 instances (3%).In addition, the average gap of the average solutions found by ILS SBRP and the lower bound reported in Erdogan et al. (2015), for each instance group, was always smaller than 0.7%, thus ratifying the robustness of our heuristic.As for the 80 instances involving 100 stations, ILS SBRP outperformed the best heuristics available for the problem by considerably improving the best known solution for all instances.
Future work may include the development of an enhanced procedure to verify whether or not the solution is feasible.Currently, this is the most time consuming part of the algorithm, where we use a relatively costly max-flow based procedure (Chemla et al., 2013a) for performing this task.Hence, any improvement on this procedure could possibly lead to an improvement on the CPU time.Also, other type of hybridizations could be experimented by combining, for example, efficient exact algorithms with the heuristic suggested in this work.

Appendix A Checking feasibility
Let L = i 1 , i 2 , ..., i k be a sequence of vertices, where i 1 = i k = 0.A directed graph can be built using p i , p i , and q i for each i in the sequence, as follows: • Let s be the source of the flow network, and for each vertex i representing the first occurrence of each station in the sequence let us define a set of arcs u i with capacity p i ; • Let t be the sink of the flow network, and for each i representing the last occurrence of each station in the sequence let us define a set of arcs w i with capacity p i ; • For each j = 2, ..., k − 1 let us define an arc b j,j+1 with capacity Q; and • If a station i is visited more than once, let us define an arc d e,e+1 with capacity q i , between the eth and (e + 1)th visits to i.
By computing an s-t maximum flow, one can find optimal bike displacements along the sequence L.
For each station i, let us define pi and p i , respectively, as the resulting s-t flow on arcs u i and w i .Flow on arcs b j,j+1 indicates the number of bikes from i j to i j+1 and flow on arcs d e,e+1 denotes the quantity of bikes remaining in a station i after the eth visit and before the (e + 1)th visit.Figure 10 depicts a flow network for a sequence L = 0, 1, 4, 2, 3, 5, 2, 4, 1, 0, where s and t correspond to the depot.

Appendix B Detailed results
This appendix presents the detailed results found by ILS SBRB , as well as those obtained by the exact algorithm of Erdogan et al. (2015), for the instances with up to 60 stations.

Figure 1 :
Figure 1: Representation of optimal solution with value 5989 for instance n20q10D

Figure 3 :
Figure 3: Handling an infeasible solution by considering additional visits to unbalanced stations

Figure 4 :
Figure 4: Example regarding the application of neighborhood structures

Figure 5 :
Figure 5: Example regarding the application of perturbation mechanisms ; Silva et al. (2015); Subramanian (2012), and Vidal et al. (2015)), this parameter is usually set as a function of the number of customers and vehicles.In TSP-like (or single machine scheduling) problems (e.g., Silva et al., 2012; Subramanian and Battarra, 2013; and Subramanian et al., 2014), I ILS was set only as a function of the number of customers (or jobs).We decided to use the same rationale as in Silva et al. (2012), by setting I ILS = max{I min , β × n}, where I min and β are input parameters.For the latter we set β = 4, as in Subramanian et al. (2014).Note that

ILS
SBRP was executed 10 times for each instance with a time limit of 1 hour.The results found by our algorithm are compared with the upper bounds determined by two versions of the tabu search heuristic ofChemla et al. (2013a).The first one (TS1) starts from an initial solution generated by a greedy procedure, while the second version (TS2) receives the solution produced by their branch-and-cut algorithm (over a relaxation of the problem) as initial solution.Detailed results of TS1 and TS2 are available at http://cermics.enpc.fr/~meuniefr/SVOCPDP.html.According to the authors, it should be noted that UB1 is the best solution found by TS1 and UB2 is the best solution found considering both TS1 and TS2.A comparison is also performed with the lower and upper bounds obtained by the exact Branch-and-cut algorithm ofErdogan et al. (2015).Regarding the benchmark instances, Chemla et al. (2013a) considered n ∈ {20, 40, 60, 100} and Q ∈ {10, 30, 45, 1000}, whereas Erdogan et al. (2015) considered n ∈ {20, 30, 40, 50, 60} and Q ∈ {10

Figures 6 QFigure 7 :
Figures 6 and 7 show how the average gaps and CPU times of each method vary according to the value of Q.We omit the results of TS1 because the associated gaps are quite inferior when compared to those obtained by the other algorithms.While the average gaps of ILS SBRP tend to be larger for very small values of Q, the average CPU time decreases as the value of Q increases, both for α = 1 and α = 3.A similar behavior regarding the CPU time performance can be observed for the heuristic ofChemla et al. (2013a), as opposed to the algorithm ofErdogan et al. (2015), which does not seem to have a consistent pattern when considering this aspect.

Figure 10 :
Figure 10: Flow network for feasibility checking

Table 1 :
Impact of different combinations of the perturbation mechanisms

Table 2 :
Erdogan et al. (2015)verage solution and the best solution, respectively, found by ILS SBRP over the 10 runs with respect to the lower bounds inErdogan et al. (2015); Avg.Time (s) is the average CPU time in seconds spent by ILS SBRP to completion over the 10 runs; Avg.TT UB2 (s) denotes the average time over the 10 runs to find or improve the best heuristic solution found in Chemla et al. (2013a) (UB2); and Avg.NV is the average number of visits of the final solutions found by ILS SBRP .Detailed results are reported in Appendix B, including the best solution found when considering all experiments (Best of all exp.).Aggregate average results per instance group for n ∈ {20, 30, 40, 50, 60} and α = 1 Erdogan et al. (2015); Time (s), UB1 Time (s), and UB2 Time (s) indicate, respectively, the CPU time in seconds spent byErdogan et al. (2015), TS1, and TS2, where the last two are scaled to our processor as mentioned above; Avg.Gap (%) and Best Gap (%)

Table 4 :
Chemla et al. (2013a)mance of the best solutions aggregated by n Tables5 and 6illustrate the detailed results found by our algorithm and the tabu searches ofChemla et al. (2013a)for every instance containing 100 stations.In this case, the gaps are computed with respect to the lower bound reported inChemla et al. (2013a).The results obtained show that ILS SBRP clearly outperforms the methods from the literature, both in terms of solution quality and CPU time.In general,