The static bike sharing rebalancing problem with forbidden temporary operations

Authors papers to INFORMS journals by means of a style ﬁle template, which includes the journal title. However, use of a template does not certify that the paper has been accepted for publication in the named journal. INFORMS journal templates are for the exclusive purpose of submitting to an INFORMS journal and should not be used to distribute the papers in print or online or to submit the papers to another publication. This paper introduces and solves the static bike rebalancing problem with forbidden temporary operations. In this problem, one aims at ﬁnding a minimum cost route in which a vehicle performs a series of pickup and delivery operations, while satisfying demand and capacity constraints. In addition, a vehicle can visit stations multiple times but cannot use them to temporary store or provide bikes. Apart from bike rebalancing, the problem also models courier service transportation and repositioning of inventory between retail stores, where temporary operations are frequently disliked because they require additional manual work and service time. We present some theoretical results concerning problem complexity and worst case analysis, and then propose three exact algorithms based on diﬀerent mathematical formulations. Extensive computational results on instances involving up to 80 stations show that an exact algorithm based on a minimal extended network produces the best average results.


Introduction
In the static bike rebalancing problem with forbidden temporary operations (SBRP-FT), we are given a directed graph G = (V, A), where V = I ∪ {0} is the set of vertices, I = {1, ..., n} is the set of stations, 0 is the depot, and A = {(i, j) : i, j ∈ V, i = j} is the set of arcs. Each arc (i, j) ∈ A is associated with a traveling cost c ij , possibly asymmetric but satisfying the triangle inequality. Each station i ∈ V has a demand d i of bikes, and is said to be in excess if d i > 0, in default if d i < 0, or 3 transport of chemical products, where temporary dropoffs might require the presence of additional protective clothing for drivers and are allowed only at sites that meet the requirements established by regulations on risk mitigation (see Cefic-The European Chemical Industry Council 2013).
To our knowledge, this is the first work to directly approach the SBRP-FT. However, its closest variant, called the static bicycle rebalancing problem (SBRP), where temporary operations are allowed, has already been studied by Chemla, Meunier, and Wolfer Calvo (2013), Erdoǧan, Battarra, and Wolfler Calvo (2015), and Cruz et al. (2017). Salazar-González and Santos-Hernández (2015) studied a general one-commodity pickup and delivery problem with split demands, and described how their model could be modified to forbid temporary operations and impose zero load on the vehicle when leaving the depot. Another realistic application close to the SBRP-FT is the bike sharing rebalancing problem, studied by Dell'Amico et al. (2014), where multiple vehicles are used but multiple visits and demand splits are not allowed. Details about these works are provided in the following section. Note that all the methods that we propose could easily deal with the cases in which demand splits are not allowed (by simply imposing each vertex to be visited once), or a service time is imposed on each delivery (by adding the service time to the traversal time of the arcs).
By forbidding temporary pickups and deliveries, the complexity of operations at stations is reduced, but routing costs might be higher than in the SBRP. As an illustrative example, consider Figure   1-(a), which shows an optimal solution for the SBRP on a benchmark instance. Each arc is traversed exactly once and the flow of commodity passing through the arc is shown over the arc. Each station is associated univocally with a vertex and the station demand is shown nearby. In the first visit to station 11, coming from station 3, instead of performing a pickup operation the vehicle temporarily deliveries 2 bikes and then continues with an empty load to station 12. In the second visit, coming from station 5, the vehicle fulfills the demand of station 11 and recollects the 2 units delivered during the first visit, thus proceeding to station 17 with a load of 7 bikes. Figure 1- (b) presents an optimal solution for the same instance, but for the SBRP-FT. Note that forbidding temporary operations results in a small increase in routing costs, from 6012 to 6025, but it also substantially reduces the complexity of operations at station 11.
In this paper, we propose three exact algorithms for the SBRP-FT, all making use of a branchand-cut (B&C) framework, but with different emphases. The first two are based on the use of an aggregate formulation derived from the work of Chemla, Meunier, and Wolfer Calvo (2013), where an integer variable expresses the number of times the vehicle passes through an arc. The solution of this formulation may be infeasible for the SBRP-FT because temporary operations might be performed in stations visited twice or more. Thus, our first algorithm iteratively removes infeasible solutions, one at a time, in a process that we call branch-and-reject. Our second algorithm is a direct extension Bruck, Cruz, Iori, and Subramanian: The static bike sharing rebalancing problem with forbidden temporary operations 4 Article submitted to Transportation Science; manuscript no.  Example of optimal solutions for the SBRP and the SBRP-FT on a benchmark instance.
of the one presented in Erdoǧan, Battarra, and Wolfler Calvo (2015) and is built upon an expanded arc structure that enables one to only use binary variables. With this structure, constraints on paths can be used to easily forbid infeasibilities. Our third algorithm removes instead infeasibilities by duplicating vertices associated with stations visited multiple times, but attempts at duplicating as few vertices as possible. Furthermore, we also present some theoretical results on problem complexity Bruck, Cruz, Iori, and Subramanian: The static bike sharing rebalancing problem with forbidden temporary operations Article submitted to Transportation Science; manuscript no.

5
and worst case analysis, as well as an effective feasibility check procedure that is used by all our algorithms.
The remainder of this paper is organized as follows. Section 2 reviews the related literature. Section 3 presents the aggregate formulation and the branch-and-reject algorithm, whereas Section 4 discusses the complexity of the problem and introduces the theoretical results. Section 5 and 6 describe the binary arc formulation and the minimal extended approach, respectively. Section 7 contains the computational results and Section 8 concludes.

Related work
Bike sharing systems are an important tool for supporting sustainable mobility and reducing urban traffic and pollution. Because of this, in recent years there has been a growing trend in self-service bike sharing systems. According to the Bike Sharing World Map (see Meddin 2016), in August of 2016 there were at least 1005 cities worldwide with a bike sharing system in operation, and more than 300 others with programs under development. This accounts for, approximately, 1.4 million shared bikes being used on a daily basis.
The SBRP-FT belongs to the class of pickup and delivery problems (PDPs), for which we refer to, e.g., the recent surveys by Battarra, Cordeau, and Iori (2014) and Doerner and Salazar-González (2014). According to the classification scheme proposed by Berbeglia et al. (2007), the SBRP-FT is a many-to-many pickup and delivery problem (M-M-PDP), in which commodities may have multiple origins and destinations. In the literature, there are a few M-M-PDPs that are related to the SBRP-FT. A notable example is the one-commodity pickup and delivery traveling salesman problem (1-PDTSP), introduced by Hernández-Pérez and Salazar-González (2004a), where a single vehicle must perform a series of pickups and deliveries of a single commodity. The 1-PDTSP is different from the SBRP-FT because: (i) all customers must be visited exactly once and thus splits are not allowed; and (ii) the depot serves as a customer and therefore the route is note forced to start and end with an empty vehicle. Hernández-Pérez and Salazar-González (2004a) proposed a mathematical formulation for both the symmetric and the asymmetric versions of the 1-PDTSP, as well as a B&C algorithm.
This approach was later improved by Hernández-Pérez and Salazar-González (2007) A dynamic variant of the BRP was studied by Contardo, Morency, and Rousseau (2012), in which the demand of stations is not static and the fleet of vehicles is heterogeneous. The authors presented a mathematical formulation and proposed a solution approach based on a time-indexed formulation on a fixed time horizon that was solved using Dantzig-Wolfe and Benders' decompositions.
Static rebalancing of bikes was studied also by Raviv, Tzur, and Forma (2013), who aimed at minimizing both operational costs and users' dissatisfaction. The latter term was modeled as a convex function that considers the expected number of shortage events (users cannot rent a bike at an empty station or cannot return a bike at a full station). To solve the problem, they proposed two mathematical formulations, strengthened them with valid inequalities and dominance rules, and conducted tests on a variety of large instances based on real data, showing that small optimality gaps could be achieved within a reasonable time. Recently, Forma, Raviv, and Tzur (2015) solved the same problem by the use of a matching-based heuristic, obtaining good quality solutions for instances with up to 200 stations.
Another static rebalancing variant was addressed by Schuijbroek, Hampshire, and van Hoeve (2017). They combined two types of decisions, namely, determining service level requirements at each bike sharing station and designing low-cost vehicle routes to rebalance the inventory. They solved the resulting decision problem by means of a cluster-first route-second heuristic, in which a polynomial-size clustering subproblem is invoked to simultaneously determine service level feasibility and approximate routing costs. They tested their algorithm on realistic data from the cities of Boston (MA) and Washington (DC).
Recently, Salazar-González and Santos-Hernández (2015) proposed the split-demand onecommodity pickup and delivery traveling salesman problem (SD1PDTSP), which is a generalization of the 1-PDTSP. In the SD1PDTSP, the stations and the depot might be visited multiple times by a single vehicle. In practice, the depot is considered as a standard customer having its own demand Bruck, Cruz, Iori, and Subramanian: The static bike sharing rebalancing problem with forbidden temporary operations Article submitted to Transportation Science; manuscript no.

7
and capacity, and the vehicle is not forced to depart from or arrive at the depot empty. In addition, a maximum number of visits is imposed to each station, and split pickups/deliveries as well as temporary operations are allowed. The authors modeled the SD1PDTSP on an extended network where each station is associated to a number of vertices and each vertex is visited at most once.
They then developed a B&C algorithm and strengthened it by the use of Benders' decomposition and valid inequalities. They also discussed how to constraint their model to force empty load on vehicles leaving the depot and to forbid temporary operations. Their approach could thus be used to solve the SBRP-FT, although at the expense of the large number of variables required by the underlying extended formulation.
To our knowledge, Chemla, Meunier, and Wolfer Calvo (2013) were the first to study the SBRP, but under the name of single vehicle one-commodity capacitated pickup and delivery problem. The authors proposed a complete formulation to define the problem, presented two relaxations to provide lower bounds, and developed a tabu search algorithm to obtain feasible solutions of good quality.
The SBRP was recently heuristically solved by Cruz et al. (2017), who obtained improved results on the existing benchmark instances by means of an iterated local search algorithm. Erdoǧan, Battarra, and Wolfler Calvo (2015) studied the same problem addressed by Chemla, Meunier, and Wolfer Calvo (2013), using the name SBRP, and another problem variant that they called non-preemptive SBRP. In the non-preemptive version, a bike can be loaded on the vehicle and unloaded at a station at most once, in contrast with the SBRP, where bikes can be temporarily loaded and unloaded any number of times. For the solution of both problem variants they proposed an exact algorithm based on a binary arc formulation (that we resume in Section 5).
Note that the non-preemptive SBRP allows temporary operations. Indeed, the situation depicted in Figure 1-(a), although infeasible for the SBRP-FT, is feasible for both the preemptive and the non-preemptive versions of the SBRP studied by Erdoǧan, Battarra, and Wolfler Calvo (2015) (the additional bikes picked up and later delivered at vertex 11 are moved only once).
We finally mention the recent work by Bruck and Iori (2017), who presented exact algorithms to deal with the class of one-to-many-to-one single vehicle routing problems with pickups and deliveries, in which a vehicle is used to serve a set of customers requiring a pickup, a delivery, or both. A notable difficulty that is encountered when solving this class of problems is that each customer can be visited up to two times. In Section 6 we propose a non-trivial way of adapting the best algorithm in Bruck and Iori (2017) (minimal extended network algorithm) to the case where demands may have multiple origins and multiple destinations, and stations may be visited more than twice, so as to efficiently deal with the SBRP-FT.

Aggregate formulation and branch-and-reject algorithm
In this section, we introduce an effective lower bound that is based on the aggregate formulation by Chemla, Meunier, and Wolfer Calvo (2013), and then discuss a first algorithm that produces an exact SBRP-FT solution.

Aggregate formulation
As unbalanced stations can be visited an arbitrary number of times, the size of the route may grow exponentially with respect to the number of stations. In the SBRP-FT, as temporary operations are forbidden and the cost matrix satisfies the triangle inequality, a station is visited only to meet at least a part of its demand. In the following, we thus impose that each station i ∈ I might be visited at most β i = |d i | times, accounting for the situation in which at each visit a single bike is delivered/collected.
Furthermore, under the same assumptions, any arc connecting two stations having both positive or negative demand values is traversed at most once. In this sense, let u ij be an upper bound on the number of traversals on arc (i, j) ∈ A, defined as Let us also define for any set Hereafter, we assume that initially balanced stations are removed from the instances. By setting x ij as an integer variable specifying the number of times arc (i, j) ∈ A is traversed, we obtain the following aggregate formulation (AF): subject to The objective function (2) minimizes the total distance traveled. Constraints (3) and (4) define, respectively, the degree for the depot and the minimum degree for the stations, while constraints of the branching tree by a twofold procedure. First, we check the existence of subtours using a fast inspection procedure that works by simply traversing the arcs of the solution. For each subtour found, if any, a cut of type (6) is added to the model. In case none are found, we apply a standard max-flow procedure that looks for violations on the vehicle capacity. Constraints (7) define the domain of the variables. In the following, we denote G = (V, A) as the aggregate network, and a solution (z, x) to AF as an aggregate solution.

Branch-and-reject algorithm
Formulation AF does not necessarily solve the SBRP-FT to optimality. Instead, it just provides a valid relaxation for the problem. This happens because AF does not consider the evolution on the number of bikes at each visit to a station. As a result, it might produce aggregate solutions with temporary operations, and others where there are no temporary operations but the associated bike displacements are still infeasible. An example of the second case, adapted from a similar example in Chemla, Meunier, and Wolfer Calvo (2013), is shown in Figure 2. In the figure, I = {1, 2, 3, 4, 5}, each depicted arc is traversed just once, and the vehicle load is reported over each arc. It is easy to see that, for Q ≥ 2, this aggregate solution satisfies all constraints (3)- (7), and is thus feasible for AF. However, it is infeasible for the SBRP-FT, because when the vehicle arrives at station 2, starting in 0 and passing through 1, it has only 1 bike to deliver as station 5 has not yet been visited. It is worth mentioning that this kind of infeasibilities may only exist in aggregate solutions and when there is at least one vertex visited more than once. Because the x variables are integer (not just binary), it is not possible to separate the standard infeasible path constraints, nor the enhanced combinatorial Benders' cuts (see, e.g., Codato and Fischetti 2006) to forbid infeasible aggregate solutions.
A straightforward approach to overcome this issue and define a complete problem formulation could be achieved by solving the SBRP-FT under a so-called extended network, obtained by duplicating Article submitted to Transportation Science; manuscript no.
each station i ∈ I into β i vertices, and imposing each vertex to be visited at most once. By doing so, we could keep track of the evolution of the number of bikes at each station during the successive visits, thus building extended solutions where the aforementioned infeasibilites for the SBRP-FT are eliminated. This was indeed the approach implemented by Salazar-González and Santos-Hernández (2015) to solve the SD1PDTSP, and one of the approaches implemented by Bruck and Iori (2017) to solve the single vehicle routing problem with deliveries and selective pickups (a one-to-many-to-one PDP in which deliveries are mandatory, pickups are optional but generate a revenue if performed, and customers are visited at most twice). However, these two works showed that formulations on the extended network become easily intractable, even for small size instances, due to the large number of duplicate vertices and the even larger number of possible connections between the vertices. Indeed, the maximum number of visits to a vertex is 2 in the tests in Bruck and Iori (2017), and 3 in the tests in Salazar-González and Santos-Hernández (2015). We thus disregarded this approach from our tests as typical β i values of our instances can be much larger than 2 or 3.
Our algorithms rely instead on the idea of building feasible SBRP-FT solutions starting from feasible AF solutions (or from reformulations/generalizations of AF). To this aim, following the same notation adopted for the SBRP by Chemla, Meunier, and Wolfer Calvo (2013), we say that an aggregate solution (z,x) to AF induces a feasible SBRP-FT solution if it is possible to find a route that travels exactlyx ij times on each arch (i, j), has costz, and satisfies all SBRP-FT constraints.
The following feasibility check is solved a number of times by our algorithms.
Definition 1. Given an aggregate solution (z,x) to AF, the disaggregation check is to determine whether or not (z,x) induces a feasible solution for the SBRP-FT.
In the following, we call disaggregate solution a feasible SBRP-FT solution that is produced by a disaggregation check. Our first solution approach, denoted as branch-and-reject (B&R) algorithm, consists in solving AF under the B&C framework provided by a mixed integer linear programming (MILP) solver, but solving the disaggregation check for any integer solution that is encountered during the process. AF solutions inducing feasible disaggregate solutions are kept and used to possibly update the incumbent solution. AF solutions for which the disaggregation check returned answer "no" are instead reported to be infeasible and hence disregarded. More in detail, there are two possible sources for an AF infeasible solution: either (i) it has been generated by the general purpose heuristics of the MILP solver, or (ii) it originates from a current integral node of the search tree. In case (i) the solution is simply rejected, whereas in case (ii) the corresponding node is fathomed without adding explicit constraints. Note that, when an integer point is declared infeasible, the branching process within a MILP solver continues, and thus B&R implicitly explores all integer points in the polyhedron and returns an optimal solution.
Details on the method that we have implemented to solve the disaggregation check are given in Section 6.3 and the B&R is computationally evaluated in Section 7.
Bruck, Cruz, Iori, and Subramanian: The static bike sharing rebalancing problem with forbidden temporary operations Article submitted to Transportation Science; manuscript no.

Complexity and worst case analysis
In this section, we discuss the computational complexity of the disaggregation check and study the worst-case performance of some important PDP problems.

Complexity of the disaggregation check
We prove the complexity of the disaggregation check by means of two consecutive results. For the sake of conciseness, all proofs of this paper are included in the appendix.
Property 1. The disaggregation check is N P-complete.
Proof. In the appendix.
This first result derives from a similar one that has been proposed for the SBRP by Chemla, Meunier, and Wolfer Calvo (2013), but it is based on a slightly different proof that allows to better understand our second and stronger result.
Property 2. The disaggregation check is strongly N P-complete.
Proof. In the appendix.

Worst case analysis
In the well-known vehicle routing problem (VRP), a fleet of vehicles must visit a set of customers, each exactly once, so as to satisfy demands and minimize costs. Archetti, Savelsbergh, and Speranza (2006) studied the relaxation in which each customer can be visited more than once (a.k.a. split-delivery VRP), thus allowing demand splits. In particular, they showed that z(VRP)/z(VRP with splits) ≤ 2, where z(·) denotes the optimal solution value, and that the bound is tight. This proves that allowing splits may halve the VRP solution cost.
A related study was conducted by Nowak, Ergun, and White (2008) on the area of one-to-one PDPs (refer again to the notation by Berbeglia et al. 2007). The authors focused on the pickup and delivery VRP (PDVRP), a VRP generalization in which each customer demand consists of transporting a load from a given pickup vertex to a given destination vertex, and conjectured that z(PDVRP)/z(PDVRP with splits) ≤ 2. To the best of our knowledge this conjecture still holds.
We investigated similar properties in the area of M-M PDPs, where, in contrast to the VRP and the PDVRP, demands may have multiple origins and destinations. We could not find a relationship between the SBRP (that allows splits and temporary operations) and the SBRP-FT (that only allows splits), neither between the SBRP and the SBRP variant in which only splits are forbidden. However, we can point out the relationship that exists between the SBRP and the SBRP variant in which both splits and temporary operations are forbidden, and show that the same result applies to the SD1PDTSP studied by Salazar-González and Santos-Hernández (2015).
Property 3. Forbidding both the demand splitting and the use of temporary operations in the SBRP may lead to an arbitrarily large increase in the solution cost, and the same holds for the SD1PDTSP.

12
Article submitted to Transportation Science; manuscript no.
Proof. In the appendix.
It is worth mentioning that, as often happens on combinatorial problems, there is quite a large gap between a proven worst case and what happens on average in practical cases. In fact, Cruz et al. (2017) have empirically shown that the difference in the solution cost, when one forbids temporary operations in the SBRP, is relatively small. They conducted experiments with and without temporary operations, and found that the increase in the solution cost was always smaller than 4% in their tests. Erdoǧan, Battarra, and Wolfler Calvo (2015) proposed an exact algorithm to solve the SBRP and further adapted it to the non-preemptive SBRP. In this section we present a direct extension of their algorithm which is capable of solving the SBRP-FT.

Binary arc formulation
As mentioned in Section 3, the main drawback that prevents AF from exactly solving the SBRP-FT is the fact of not considering the evolution of the number of bikes at stations after each visit.
Provided that the triangle inequality holds, we can use (1) to transform each integer variable x ij in AF into a set of binary variables as where y k ij is a binary variable assuming value 1 if arc (i, j) is traversed at least k times by the vehicle. Equations (8) can be intuitively seen as a representation of the value of x ij as a binary number, where each variable y k ij gives the k-th bit of the number (see, e.g., Vanderbeck and Wolsey 2010). For the sake of simplicity, let γ ij = ⌊log 2 (u ij )⌋. We now create a multigraph in which each arc (i, j) ∈ A is replaced by a set of arcs (i, j, k) with k = 0, 1, . . . , γ ij . Let R be the set of routes that are infeasible for the SBRP-FT and A ′ (r) be the subset of arcs in the multigraph used by route r ∈ R. We are now ready to present the following binary arc formulation (BinArc).
Bruck, Cruz, Iori, and Subramanian: The static bike sharing rebalancing problem with forbidden temporary operations Article submitted to Transportation Science; manuscript no.

13
The objective function (9) minimizes the total distance traveled. Constraints (10)-(12) impose the degrees for the depot and the stations. Constraints (13) specify that each arrival at a vertex is followed by a departure, while (14) are capacity constraints for the vehicle. Infeasible routes are forbidden by constraints (15), whereas constraints (16) define the domain of the y variables. Note that constraints (12) are not strictly necessary, but are used to strengthen the formulation, given that γ ij k=0 2 k y k ij might be greater than u ij . In practice, the model made by (9)-(14) and (16) is a reformulation of AF, and thus a valid relaxation, but the inclusion of (15), made possible by the use of the binary expansion, allows to obtain a complete SBRP-FT representation. To separate (14), we first use equation (8) to aggregate the y variables, and then we apply a standard max-flow procedure. To separate (15), we solve the disaggregation check using the method described later in Section 6.3. In short, formulation (9)-(16) is solved by a B&C procedure, called BinArc algorithm, that separates constraints (14) and (15) only at integer nodes of the branching tree.

The minimal extended network approach
In this section, we introduce a formulation capable of modeling any state of the problem network, from the aggregate network, where each station is associated with a single vertex, to the extended network, where each station i is associated with β i vertices. We then use the formulation as a basis of an iterative B&C algorithm.

General formulation
Let G ′′ = (V ′′ , A ′′ ) be a complete and directed graph, where V ′′ = ∪ i∈I V i ∪ {0} is the set of vertices, V i = {i 1 , . . . , iβ i } is the subset of vertices associated with station i ∈ I, andβ i is a parameter satisfyinḡ β i ≤ β i for each i ∈ I. In addition, let I k specify the station associated with vertex k ∈ V ′′ . To represent the evolution that our iterative algorithm imposes on this network, set V ′′ is partitioned as As an illustrative example consider the vertex set with depot and 4 stations depicted in Figure 3-(a). This set contains only aggregate vertices that belong to V a and can be visited multiple times. Suppose that vertex 2 is duplicated by our iterative procedure, as shown in Figure 3-(b). The initial aggregate vertex is replaced by a first extended vertex belonging to V f , and a partially aggregate vertex belonging to V p . Vertices in V f are visited at most once, while those in V p may be visited multiple times. Also, a partially aggregate vertex may only be visited when all other vertices associated with the same station are also visited. Suppose now that a second duplication of vertex 2 occurs, as shown in Figure 3-(c). This results in an extended vertex, belonging to V e , which Article submitted to Transportation Science; manuscript no. may be visited at most once. To better understand the difference between V e and V f , consider two vertices, k ∈ V f and l ∈ V e , having I k = I l = i ∈ I. The only difference between k and l is that the former can only be visited when the latter is also visited. Additionally, in case d i = 0, vertex k must be visited exactly once. Subsequent duplications of the same vertex, not explicitly reported in the figure, would maintain exactly one vertex of type f and one vertex of type p, and increase only the number of vertices of type e.

Figure 3 Example of a vertex duplication
Similarly to AF, let x ij be a set of integer variables specifying the number of times arc (i, j) in A ′′ is traversed, and let y i be a set of binary variables that define whether or not vertex i ∈ V ′′ \ {0} is visited. In addition, let g i be an unconstrained variable that reports the quantity of demand that has been collected/delivered during the visit to vertex i ∈ V ′′ \ {0}. We are now ready to introduce the following general formulation (GF), which is the stepping stone of our algorithm: subject to x(δ + (0)) = 1 (18) The objective function (17)  Constraints (24) impose the order for the usage of the duplicates, stating that a vertex i k+1 can only be visited if its predecessor i k is also visited. The demand d i of each station must be served, as stated by constraints (25), and the vehicle capacity Q must not be exceeded, as imposed by (26) and (27).
Note that constraints (26) and (27) also serve the purpose of preventing subtours and are separated using the same procedure applied to constraints (6). The max function in the right hand side of (26) induces a nonlinear term, which is required because − i∈S g i could be 0 or lower. We deal with this by simply checking at each call to the separation procedure which term is the highest, and then using it as the right hand side of the constraint that is actually added to the model. Constraints (28) and (29) specify lower and upper bounds for the g i variables, which should be equal to 0 when y i = 0, lower than 0 when y i = 1 and d I i < 0, and greater than 0 when y i = 1 and d I i < 0. It is worth mentioning that constraints similar to (28) This formulation has some interesting properties. In case each station is associated with a single (aggregate) vertex, we have V ′′ = V a ∪ {0} and V p = V f = V e = ∅, which corresponds to the scenario adopted for AF. Moreover, g j = d I j for each j ∈ V ′′ , and consequently constraints (26) and (27) are equivalent to (6), and constraints (28) and (29) can be removed, as they become redundant. This leads to: Observation 1. In case each station is associated with a single vertex in G ′′ , GF is equivalent to AF. Now recall that the shortcomings that might render an aggregate solution infeasible may only exist when there is at least one vertex being visited multiple times. By associating each station i ∈ I with β i duplicates, we obtain the extended network, where each vertex is visited at most once. Therefore, solving GF under this network always produces a feasible solution, and hence: Observation 2. In case each station i ∈ I is associated with β i duplicates in G ′′ , an optimal solution of GF is also an optimal solution for SBRP-FT.
Article submitted to Transportation Science; manuscript no.

MEN algorithm
Our third exact algorithm, called the minimal extended network (MEN) algorithm, generalizes the best algorithm in Bruck and Iori (2017) to deal with the case where demands have multiple origins and destinations and stations may be visited more than twice. The MEN algorithm is based on the following two observations. Firstly, most often than not, in optimal solutions for the SBRP-FT stations are visited no more than 2 or 3 times (this behavior was already hinted by Salazar-González and Santos-Hernández (2015) for the related SD1PDTSP). Secondly, by performing extensive computational experiments, we have verified that quite often an optimal solution of AF corresponds to an optimal solution for the SBRP-FT.
Algorithm MEN starts by solving GF over the pure aggregate network, where each station is associated with a single vertex. Formulation GF is solved by using a B&C algorithm in which constraints (26) and (27) are separated with the max-flow procedure described in Chemla, Meunier, and Wolfer Calvo (2013). We opted to separate those constraints only at integer nodes. A disaggregation check is then performed on the solution obtained by the B&C to determine whether the solution is feasible for the SBRP-FT (details are provided in Section 6.3). If feasible, then MEN terminates with an optimal SBRP-FT solution. Otherwise, we inspect the solution and determine the numberβ i of visits to i.
Then, a new iteration solves GF under the modified graph. The procedure iterates until the optimal solution found is feasible for the SBRP-FT.
Recall that AF provides a valid relaxation of the SBRP-FT. Because of this fact and of Observation 1, algorithm MEN provides at the first iteration a valid lower bound. Then, in the worst case scenario, it iterates for i∈I β i times, performing a single duplication at each iteration, possibly improving the lower bound value, and culminating in the extended network in the last iteration, where, because of Observation 2, it is guaranteed to produce a feasible SBRP-FT solution. Consequently, the solution value would be both a lower and an upper bound, and hence an optimum, so we can conclude that: Observation 3. Algorithm MEN converges to the optimal solution of the SBRP-FT.
According to Observation 3, algorithm MEN might be as time consuming (or even more) as directly solving the extended network. However, computational experiments indicate this is not the case and, in practice, only a few iterations are usually necessary to converge to an optimal SBRP-FT solution.
Furthermore, to speed up convergence, at each iteration of MEN a list of the 10 best incumbent solutions found is kept. In case the solution found is infeasible, this list is checked as an attempt of finding a feasible upper bound. The incumbent upper bound is then used as a warm start for the next iteration. During computational experiments we observed that, for some instances, the first iteration of MEN produces an infeasible aggregate solutionx, but the reverse solutionx, havingx ij =x ji for It is worth mentioning that an extended network approach was already used by Salazar-González and Santos-Hernández (2015) to model the related SD1PDTSP, so the most fruitful contribution of algorithm MEN consists in iteratively solving the problem, starting from the simplest aggregate network and then iteratively extending the graph through additional vertices until an optimum is found. To this regard, note that there is no guarantee that after the execution of algorithm MEN the resulting graph will be of overall minimum cardinality. In fact, finding a graph of minimal cardinality that could allow MEN to produce an optimal solution appears to be as hard as solving the SBRP-FT itself (and also finding non-trivial bounds on such cardinality is not easy).
Even duplicating just one vertex at a time, instead of duplicating all vertices i withβ i > 1, has no guarantee to lead to a minimum cardinality.
It is also worth mentioning that the MEN algorithm can be easily extended to deal with the case of rebalancing by means of a fleet of multiple vehicles. This can be obtained by simply replacing "1" with the number of available vehicles in constraint (18). Note, indeed, that because we forbid temporary operations there is no transshipment, and this in turn prevents synchronization issues that could require additional constraints.
In addition, defineG = (Ṽ ,Ã) as the support graph obtained by duplicating each vertex i ∈V >1 intō β i duplicates, and including all vertices inV 1 directly intoṼ . Let σ(i) specify the vertex in V ′′ that is associated with i ∈Ṽ . We check whether or not (z,x,ḡ) leads to a disaggregate feasible solution for the SBRP-FT by solving a formulation, called duplicate and inspect (DIF) that uses the same variables x and g as in GF, but is defined overG as follows: (DIF) min 0 Article submitted to Transportation Science; manuscript no.
Constraints (35) ensure that the search is performed for a disaggregate solution of the same valuez of the original aggregate one, whereas constraints (36) and (37)

Computational experiments
Our algorithms have been coded in C++ and the computational experiments have been run on a PC with an Intel Core i7-3770 3.40 GHz with 8 GB of RAM. A time limit of 1 hour has been given to each execution. We have used CPLEX 12.6.2 with default options and single thread to solve the formulations and to implement the B&C algorithms.
We first tested our algorithms on a set of instances obtained by slightly modifying those proposed by Chemla, Meunier, and Wolfer Calvo (2013) for the SBRP. The results on these instances are summarized in Section 7.1. We then selected our two best solution methods and used them to solve

19
to a station), there were no significant efficiency issues with the duplicate-and-inspect procedure.
However, this was not the case for both the B&R and the BinArc algorithms, as the check is called every time an incumbent solution is found, and the first solutions found are likely to be of poor quality (possibly containing several stations that are visited multiple times). We overcame this issue by providing an initial solution of good quality. In view of this, we have adapted to the SBRP-FT the iterated local search (ILS) algorithm proposed by Cruz et al. (2017) for the SBRP, and used it to provide an initial solution (called UB 0 ) for both the B&R and the BinArc algorithm. This has been obtained by using the MIPstart procedure provided by CPLEX. The two resulting algorithms are called B&R+MIPstart and BinArc+MIPstart below. In addition, we tested two variants for MEN: a first one, called MEN+MIPstart, in which UB 0 is used as a MIPstart when solving the GF models; and a second one, simply called MEN, in which UB 0 is only used to provide an output solution in case MEN did not find any lower cost solution.

Computational results on SBRP benchmark instances
We have slightly modified the set of benchmark instances proposed by Chemla, Meunier, and Wolfer Calvo (2013) for the SBRP and later used by Erdoǧan, Battarra, and Wolfler Calvo (2015), by imposing triangle inequality on the cost matrix  Two interesting comments can be made when comparing MEN+MIPstart with MEN. On one side, it is interesting to observe that providing an initial solution for the MEN algorithm did not improve its overall efficiency. In fact, for several instances the opposite behavior was verified. This might be explained by the observation that GF usually finds feasible upper bounds quite fast, but takes longer to improve the lower bound, and providing a warm start might change the branching decisions taken by CPLEX. The fact that performing modifications in the initial conditions of a MILP solver may lead to significant changes in the resulting performance has already been noticed by Fischetti and Monaci (2014) and Lodi and Tramontani (2013), among others. On the other side, there are instances for which MEN alone cannot provide any feasible solution. That happened for 2 instances with n = 40, 5 instances with n = 50, and 10 instances with n = 60 (we refer to the appendix for details). In these cases, the use of a good initial heuristic is very important.  Table 2 it can be also verified that the larger the capacity, the easier the instance. In fact, instances having Q larger than 35 are not very meaningful, as they are easily solved by both

Computational results on realistic benchmark instances
We also tested our best algorithms on a set of real-world instances that were first proposed by the case in which multiple homogeneous vehicles can be used and each station is visited once. We adapted the instances to the SBRP-FT by assuming that a single vehicle is used, removing initially balanced stations, and allowing the unbalanced stations to be visited more than once. In addition, to deal with the fact that, originally, these instances are not balanced (i.e., i∈I d i = 0), we added a dummy station at the same coordinates of the depot and set its demand as − i∈I d i . The vehicle capacity was kept as in Dell'Amico et al. (2014), who attempted up to 3 different values for each graph. Note that we limited our study to the 56 instances having at most 80 stations, because this is the computational limit of our exact methods and also because it is unrealistic to consider that a single vehicle can visit more stations in a single day. with Buenos Aires have a small number of stations, but seem harder to solve than some larger instances. This is probably due to the fact that the demands of stations in Buenos Aires are very high, and the vehicle capacity is very tight, and this makes the problems more difficult to solve. On the other hand, instances associated with Guadalajara have a larger number of stations but very low demand values, and are easier to solve. In Table 4 we aggregate the results by value of Q. Also here one can notice that small vehicle capacities tend to increase the solution difficulty. We also performed an additional test by attempting Q=1000, and could solve all instances with both MEN+MIPstart and MEN in at most a few minutes. We also decided to perform additional tests that follow the observation by Liu et al. (2016), who pointed out that in many practical applications the demand distribution is usually geographically Bruck, Cruz, Iori, and Subramanian: The static bike sharing rebalancing problem with forbidden temporary operations Article submitted to Transportation Science; manuscript no.

23
unbalanced. For example, as people usually leave their residencies in the morning, residential areas are expected to have a shortage of bikes in this moment of the day, while commercial regions have a surplus. Similarly, in the late afternoon this scenario is usually reversed. To tackle this type of scenarios, we created a new set of instances where the demand of stations is distributed unevenly. These instances were generated by modifying the demands of the instances proposed by Chemla, Meunier, and Wolfer Calvo (2013) using the following procedure. For each instance, we first divide the geographic space into four quadrants as in the two-dimensional Cartesian system. Next, in a new set of 25 instances with realistic spatial distribution of demand. Note that the procedure we adopted tends to generate demands with values that are much higher than those of the other instances that we tested. The resulting instances are expected to be very challenging because split deliveries are forced by the fact that some demands are larger than the vehicle capacity.
The outcome of the computational experiments performed on this new set of instances is presented in Table 5. Detailed results are reported in the appendix. As expected, the instances are more challenging to solve. This can be noted by the reduced quality of the initial upper bound provided by the heuristic, which now never manages to reach the value of a proven optimal solution, and by the reduced scalability of the MEN algorithm, which now can solve to optimality only instances having 40 stations or less. The positive aspect is that the developed algorithms can directly tackle instances with demands that are higher than the vehicle capacity, showing a good adaptability to a wide range of problems.

Conclusions
In this paper, we have introduced the static bike rebalancing problem with forbidden temporary operations (SBRP-FT), a pickup and delivery problem that arises in bike rebalancing, courier service transportation, and repositioning of inventory cases in which temporary dropoff/pickup of products is forbidden. We have presented some theoretical results related to computational complexity and worst case analysis and proposed three exact algorithms. The first one consists of a branch-and-reject (B&R) algorithm that builds upon an incomplete integer programming formulation. The second one is an integer programming model, denoted as BinArc formulation, that is based on a binary network expansion. The third one is based on the so-called minimum extended network (MEN) algorithm, which iteratively enlarges a relaxed formulation by duplicating as few vertices as possible.
We have first performed extensive computational experiments with the proposed exact methods on a set of 450 benchmark instances. Although the BinArc formulation is usually outperformed by the other two approaches, it is capable of proving the optimality of 278 solutions, while B&R and MEN (using its best configuration) obtain 376 and 402 optimal solutions, respectively. The good behavior of the MEN algorithm was also proven by a second set of tests that we performed on realistic instances taken from the bike sharing literature.
Future work includes the application of the MEN algorithm to other routing problems with split demands, and the development of hybrid algorithms by combining heuristic and exact approaches to find improvements for the SBRP-FT open instances. Such instances are mainly those with small capacity values and in practice they tend to be harder to solve because the feasible solutions usually contain several multiple visits to the stations.