# A Bilevel Programming Model and Algorithm for the Static Bike Repositioning Problem

Authors: Qiong Tang; Zhuo Fu; Meng Qiu
Journal: Journal of Advanced Transportation
 (2019)

---

## Abstract

In this paper, by taking the outsourcing transportation mode into account, a bilevel programming model is proposed to formulate the static bike repositioning (SBR) problem, which can be used to determine the number of bikes loaded and unloaded at each station and the optimal truck routes in bike sharing systems (BSS). The upper-level BSS providers determine the optimal loading and unloading quantities at stations to minimize the total penalties. The lower-level truck owner pursues the minimum transportation route cost. An iterated local search and tabu search are developed to solve the model. Computational tests on a set of instances from 20 to 200 bikes demonstrate the effectiveness of the model and algorithms proposed, together with some insightful findings.


---

## Body

## 1. Introduction

Bike sharing systems (BSS) provide a carbon-free, green mode of transportation for urban users and reduce urban traffic congestion and air pollution, so they are very popular in urban environments at present [1]. One of the major problems faced by these BSS is the operational issue of repositioning bikes between different stations. Demand variability often causes the bikes in stations to become an excess of inventory or too little inventory to service new customers effectively. This not only affects the desired service level but also incurs spurious operational costs. [2]. In practice, aiming to increase customer satisfaction, BSS providers always need to transport bikes between stations by trucks to meet the customer demand. This relocation problem is called a bike repositioning problem (BRP), which determines the number of bikes loaded and unloaded at each station and optimal truck routes, subject to various constraints including truck capacity, station capacity, and operational time constraints. To reduce the operational costs, many BSS providers entrust the relocation task to third-party logistics companies (the truck owners). In such a real-life operation situation, the BSS providers determine the number of bikes loaded and unloaded at each station to maximize customer satisfaction with regard to the truck capacity constraint, and the truck owner determines the optimal truck routes to pursue the minimum transport scheduling cost based on the pickup and drop-off quantities at each station provided by the BSS providers. In the BRP under consideration, the BSS provider determines the pickup and drop-off quantities at each station subject to the station capacity and truck capacity constraints and seeks to optimize customer satisfaction, while the truck owner determines the truck routes subject to the truck capacity and operational time constraints and seeks to optimize the transportation cost. Once the BSS provider sets the pickup and drop-off quantities at each station, the truck owner reacts by providing the truck routes that optimize his objective function. In this realistic background, the BRP involves decision makers at two distinct levels with a hierarchical relationship. First, the BSS provider, at the upper level of the hierarchy, controls the loading and unloading quantities variables and is influenced by the resource constraints of the truck owner. Second, the truck owner, at the lower level of the hierarchy, controls the truck route variables and is limited by the decisions of the BSS provider. Such unique characteristics mean that the BRP can be modeled by using bilevel programming, which has been proposed in the literature as an appropriate model for hierarchical decision processes with two decision makers.To the best of our knowledge, only this paper addresses the problem of the balance between the BSS provider (e.g., Mobike) and truck owner (e.g., third-party logistics companies). In order to demonstrate a realistic BRP between the BSS provider and truck owner, in this paper, we develop a bilevel programming model to research the BRP, seeking to satisfy both of the decision makers, and to optimize customer satisfaction and the transportation cost between the upper- (BSS provider) and lower- (the truck owner) level decision makers.The remainder of this paper is organized as follows. In the next section, we review related work on the BRP. In Section3, we present the assumptions, notations, and bilevel programming model formulation. The solution algorithm is proposed in Section 4. Numerical examples are designed to demonstrate the efficiency of the proposed solution algorithm and the computational results are presented in Section 5. Finally, conclusions and further studies are set out in Section 6.


## 2. Literature Survey

The BRP has been constantly studied since the appearance of the seminal publication on it by Benchimol et al., in which the authors proved it to be NP-hard [3]. The BRP falls into a larger class of problems that is called the pickup and delivery vehicle routing problem (PDVRP), which generalizes the capacitated vehicle routing problem (CVRP) [13]. Furthermore, according to the classification introduced by Parragh et al., the BRP is a many-to-many pickup and delivery problem (M-M), where the origins and destinations of the requests are multiple [31]. When the BRP considers only one truck, the problem becomes the one-commodity pickup-and-delivery TSP introduced by Hernández-Pérez and Salazar-González [32].Generally, the models proposed for the BRP in the literature are classified into static and dynamic optimization problems according to the operational scenarios [33]. The static problem considers night-time scenarios when the usage rate of the systems is negligible, while the dynamic problem considers daytime scenarios that take real-time usage of the system into account. Only a few publications have tackled the dynamic BRP [1, 34–37].The publications on the static BRP are summarized in Table1. The static BRP is motivated by the fact that very few bikes are used at night [4]. Raviv et al. [33] pointed out that the static mode of operation benefits from the practical advantage that it allows the repositioning truck to travel swiftly around the city without aggravating the traffic congestion and parking problems and that static repositioning is useful for arranging the inventory of bikes in the BSS for the next day. In terms of the number of trucks used in static BRP publications, both single and multiple vehicle problems are studied. However, compared with multiple vehicle BRP [5–12, 14–22], it is more realistic to consider a single vehicle BRP for a part of the city BSS [3, 4, 23–30]. In this paper, our study focuses on the single vehicle static BRP, which needs to be solved once at the beginning of every night, based on the status of the BSS at that time and the demand forecast for the next day and when the city is divided into districts.Table 1
Summary of the static BRP in the literature.
References Number of trucks Objectives Pickup and dropoff quantities as decision variables Benchimol et al. [3] 1 Minimize total travel cost No  Cruz et al. [4] 1 Minimize route cost No  Lin and Chou [5] ≥1 Minimize total transportation time or distance No  Di Gaspero et al. [6, 7], Raidl et al.[8], Rainer-Harbach et al. [9, 10] ≥1 Minimize the weighted sum of the total absolute deviation from the target number of bikes, the total number of loading/unloading quantities, and the overall time required for all routes Yes  Di Gaspero et al. [11] ≥1 Minimize the weighted sum of the total absolute deviation from the target number of bikes, and the travel times Yes  Angeloudis et al. [12] ≥1 Minimize total travel time Yes  Dell’Amico et al. [13, 14] ≥1 Minimize total route cost No  Forma et al. [15] ≥1 Minimize the weighted sum of the penalty cost and total travel time Yes  Alvarez-Valdes et al. [16] ≥1 Minimize the weighted sum of the total service time and the coefficient of variations of the duration of all routes Yes  Ho and Szeto [17] ≥1 Minimize the weighted sum of the penalty cost and total travel time Yes  Schuijbroek et al. [18] ≥1 Minimize maximum tour length Yes  Pal and Zhang [19] ≥1 Minimize the makespan of the rebalancing fleet (maximum tour length) No  Bulhões et al. [20] ≥1 Minimize total route cost No  Szeto and Shui [21] ≥1 Minimize first the positive deviation from the tolerance of total demand dissatisfaction and then service time in the form of the total service time or maximum service time of all vehicles Yes  Liu et al. [22] ≥1 Minimize the weighted sum of the total number of unsatisfied customers, the inconvenience of getting a bike and total operational time Yes  Chemla et al. [23] 1 Minimize route cost No  Nair et al. [24] 1 Minimize total redistribution cost Yes  Erdoğan et al. [25] 1 Minimize total route and handling cost Yes  Ho and Szeto [26] 1 Minimize total penalty cost Yes  Erdoğan et al. [27] 1 Minimize route cost No  Kadri et al. [28] 1 Minimize the total waiting time of the imbalanced stations No  Li et al. [29] 1 Minimize the sum of travel, imbalance, substitution, and occupancy costs Yes  Szeto et al. [30] 1 Minimize the weighted sum of the total number of unsatisfied customers and the vehicle’s total operational time Yes  This paper 1 The upper-level model minimizes the total penalty cost; the lower-level model minimizes the travel cost YesIn terms of the problem objective functions addressed in the BRP, common measures are(1) transportation cost or other related measures such as transportation time or distance, and (2) customer dissatisfaction with the BSS, in the form of absolute deviation from the target optimal number of bikes, penalty cost, and the total number of unsatisfied customers or unmet demands. If all the demand of stations is a hard constraint to be satisfied, the objective of the BRP is only measured by the transportation cost [3–5, 13, 14, 20, 23, 27]. In reality, due to the uncertainty of demand, capacity limits of the trucks and stations, and repositioning time, not all the demand of stations can be satisfied. That is to say, the truck can only service a subset of unbalanced stations using limited repositioning resources, so it is necessary for the BSS provider to determine the loading instruction at all stations with the aim of achieving the highest possible customer satisfaction to maximize its revenue and profit. In this paper, similar to Raviv et al. [33], Ho and Szeto [17, 26], and Forma et al. [15], we measure customer dissatisfaction by a penalty function. Details of the penalty function can be found in Section 3.2.The challenges that arise in the static BRP necessitate efficient methodologies. Exact methods such as the branch-and-cut algorithm have been used to solve the static BRP. However, such methods are intractable for dealing with medium and large realistic BRP. We refer the interested reader to the exact methods papers by Erdoğan et al. [25, 27] and Kadri et al. [28]. As summarized in Table 2, several types of heuristics have been used to solve the BRP, for instance, TS (tabu search) [23, 26], GA (genetic algorithm) [29], CRO (chemical reaction optimization) [22, 30], ACO (ant colony optimization) [6], ABC (artificial bee colony algorithm) [21], VND (variable neighborhood descent) [9, 19], VNS (variable neighborhood search) [8–10], LNS (large neighborhood search) [7, 11, 17], ILS (iterated local search) [4, 20], and GRASP (greedy randomized adaptive search procedure) [10].Table 2
Summary of inexact solution methods for the static BRP in the literature.
References Solution methods Benchimol et al. [3] 9.5-approximation algorithm  Dell’Amico et al. [13] Destroy and repair algorithm  Cruz et al. [4] ILS  Lin and Chou [5] Classical heuristics  Di Gaspero et al. [6] Constraint programming with ACO  Di Gaspero et al. [7, 11] Constraint programming with LNS  Raidl et al. [8] VNS  Rainer-Harbach et al. [9] VNS with VND  Rainer-Harbach et al. [10] GRASP with VNS  Angeloudis et al. [12] Multistage heuristic  Forma et al. [15] 3-step math heuristic  Alvarez-Valdes et al. [16] Heuristic based on minimum cost flow  Ho and Szeto [17] Hybrid LNS  Schuijbroek et al. [18] Cluster-first route-second heuristic  Pal and Zhang [19] Nested LNS with VND  Bulhões et al. [20] Branch-and-cut with ILS  Szeto and Shui [21] Enhanced ABC  Liu et al. [22] Enhanced CRO  Chemla et al. [23] Branch-and-cut method with TS  Ho and Szeto [26] Iterated TS  Li et al. [29] Combined hybrid GA  Szeto et al. [30] Enhanced CRO  This paper ILS and TSTo deal efficiently with the proposed bilevel model, a heuristic based on ILS that iteratively solves the upper-level and lower-level models is developed. The proposed model and heuristic are evaluated using a test instance, and the results show the effectiveness of the algorithm.It appears that no previous work has addressed the BRP from this approach where the objectives of the different stakeholders are taken into consideration. In this paper, we intend to close a significant gap in the BRP literature and make a twofold contribution. First, the BRP is modeled as a mixed-integer, bilevel program, where the different objectives of the BSS provider and truck owner are considered. Second, we provide a mathematical formulation and heuristic algorithm to solve the introduced problem. While the bilevel model presented here is customized for the BRP under the situation that the BSS provider outsources the transportation work between bike stations to the third-party logistics companies, we hope that it will open a new line of research on the BRP.


## 3. Problem Definition and Model Description

In this section, first of all, we describe different aspects of the problem and then present its mathematical model. The basic idea of this bilevel model in the BRP is to design relocation transportation with the object of maximizing customer satisfaction with the BSS provider by considering the repositioning resources of third-party companies. We model this decision-making environment as a bilevel problem, so that, in this structure, the BSS provider is considered as the decision maker of the upper-level problem (ULP) and the truck owner as the decision maker of the lower-level problem (LLP) (see Figure1). The BSS provider as the ULP decision maker seeks the optimal loaded and unloaded quantities of bikes at stations. The ULP objective function minimizes the total penalties cost. This function is dependent on the penalty function. The truck owner as the LLP decision maker aims to minimize the travel costs.Figure 1
The conceptual structure of bilevel bike-sharing repositioning problem.The static BRP considers a group of stations and a depot where the truck is originally located at the depot. If the optimal inventory (i.e., the inventory with the lowest penalty cost) at a station is smaller than the initial inventory at that station, then the station is assumed as a pick-up station. The initial inventory at pick-up station is too large and the extra bikes need to be taken away. On the other hand, if the optimal inventory at a station is greater than the initial inventory at that station, then the station is assumed as a drop-off station [17, 26]. The initial inventory at drop-off station is too small and the insufficient bikes need to be replenished. The repositioning truck with constrained capacity transports bikes from pick-up stations to drop-off stations.


### 3.1. Assumptions and Notation

To facilitate the presentation of the essential ideas without loss of generality, the following basic assumptions are made in this paper:( 1 ) The outsourcing cost of BSS provider outsourcing relocation transportation to the third-party logistics companies is fixed.( 2 ) Each station is allowed to be visited at most once, as in the literature [26, 33], which means that the demand of a station is not allowed to be split. Similar to Forma et al. [15] and Ho and Szeto [17, 26], not all stations need to be visited, and initially balanced stations do not need to be visited.( 3 ) The depot is assumed to have no bikes.( 4 ) The load on the truck cannot be greater than the truck capacity.( 5 ) The truck starts from and returns to the depot and operates within a given time constraint.The notations used throughout this paper are listed in Table3 unless otherwise specified.Table 3
Notations for the bi-level model.
N Set of stations  N o Set of nodes, including the stations and the depot (indexed by 0)  i, j Indices of stations  s i 0 Initial inventory level at stationibefore the repositioning operation  s i Final inventory level at stationiafter the repositioning operation  f i ( s i ) Convex penalty function ofsi for stationi  s i I The target optimal inventory at stationi, argminsifi(si),i.e., the optimal number of bikes at station i  P Set of pick-up stations, i.e.,{i∈N∣si0>siI}  D Set of drop-off stations, i.e.,{i∈N∣si0<siI}  Q Truck capacity  c i Capacity of stationi  T Repositioning time  L Time required to load a bike from a station onto the truck  U Time required to unload a bike from the truck to a station  M A very large number  t i j Traveling time from stationito station j  x i j Binary variable that equals 1 if the truck travels directly from stationito stationj, and 0 otherwise  y i j Number of bikes on the truck when it travels directly from stationito stationj  y i P The loading quantities at stationi,  yiP≥0,∀i∈P,yiP=0,ifi∉P  y i D The unloading quantitiesat stationi,  yiD≥0,∀i∈D,yiD=0,ifi∉D  U The vector of the upper-level decision variables,U={yiD,yiP,si,∀i∈N}  L The vector of the lower-level decision variables,L={xij,yij,∀i∈N,j∈N,i≠j}


### 3.2. Penalty Function

The penalty function used in this paper was first proposed by Raviv et al. [33]. Similar to Ho and Szeto [17, 26] and Forma et al. [15], we used it to define customer dissatisfaction. For any stationi, the penalty functionfi(si)is a function of the final inventory of stationi after the operation, which represents the anticipated shortage of stationi. The penalty cost of stationi can be calculated by the corresponding penalty functionfi(si), and the value range ofsiis 0,…,ci.Ifsi is out of the domain 0,…,ci, thenfi(si)=∞.As described by Forma et al. [15], all stations of the BSS were divided into three categories, namely, (1) residential stations where the demand of customers to rent bikes is large in the morning and that to return bikes is large in the afternoon; (2) business stations where the demand of customers to return bikes is large in the morning and that to rent bikes is large in the afternoon;  (3)  ‘‘touristic’’ stations where the demand of customers to rent or return bikes does not change much throughout the day. The optimal inventory level at stationi can be solved by argminsifisi.


### 3.3. Bilevel Model

In order to demonstrate the relationship between the BSS provider and the truck owner, a bilevel programming model modified on the basis of Ho and Szeto [26] is proposed. The BSS operator’s decision variables are si, yiP, and yiD. The upper-level model aims to minimize the total penalties cost. The lower-level model aims to minimize the transportation cost with the truck capacity constraint and the repositioning time constraint.


#### 3.3.1. The Upper-Level Model

The upper-level model aims to determine the optimalU at each station to minimize the total penalty costs and can be formulated as follows:(1)min∑i∈Nfisi(2)s.t.si=si0-yiP-yiD,∀i∈N(3)∑j∈No,j≠iyij-∑j∈No,j≠iyji=yiP-yiD,i∈N(4)yij≤Qxij,∀i,j∈No,i≠j(5)∑i∈PyiP-∑i∈DyiD=0(6)yiP≤si0,i∈P(7)yiD≤ci-si0,i∈D(8)yiP≥0,∀i∈P(9)yiD≥0,∀i∈N(10)si≥0,∀i∈N.The upper-level objective function (1) minimizes the total penalty costs generated by all stations. Constraint (2) defines the final inventory at stationiafter the operation as the initial inventory at stationiminus the loading quantities or plus the unloading quantities at stationi. Constraint (3) requires that the loading or unloading quantities at stationiare equal to the difference between the loads of the truck before and after the stationi. Constraint (4) specifies the upper bound of the truck’s load. Constraint (5) imposes that all bikes loaded from pick-up stations by the truck are finally unloaded at drop-off stations. Constraint (6) requires that the loading quantity at pick-up station should not exceed the initial inventory at that station. Constraint (7) requires that the unloading quantity at drop-off station should not exceed the residual capacity of that station. Constraints (8)-(10) make sure the variables are nonnegative. Decision variables xij can be solved from the lower-level model, while constraints (3) and (5) describe the relationship between the upper-level decision variables U and lower-level decision variables L.


#### 3.3.2. The Lower-Level Model

The lower-level model aims to determine the optimalL to minimize the transportation cost and can be formulated as follows:(11)min∑i∈N0∑j∈N0,i≠jtijxij(12)s.t.∑j∈Ny0j=∑j∈Nyj0=0(13)∑j∈No,j≠ixij≤1,i∈N(14)∑j∈No,j≠ixij-∑j∈No,j≠ixji=0,i∈N(15)∑i∈NLyiP+UyiD+∑i,j∈N0,i≠jtijxij≤T(16)qj≥qi+1-M1-xij,∀i,j∈No,i≠j(17)qi≥0,∀i∈N0(18)xij∈0,1,∀i,j∈N0,i≠j(19)yij≥0,∀i,j∈N0,i≠j.The lower-level objective function (11) minimizes the total travel time. Constraint (12) forces the fact that the truck leaves and returns to the depot empty. Constraint (13) imposes that every station can only to be visited once. Constraint (14) guarantees that it is impossible for the truck to stay at any station. Constraint (15) imposes that the total handling and travel time of the truck should not exceed the given repositioning time. Constraint (16) eliminates the subtours [38]. Constraint (17) ensures that the auxiliary variable is nonnegative. Constraint (18) declares that xij is a binary variable. Constraint (19) ensures that the load of the truck is nonnegative.Constraint (15) describes the relationship between the upper-level decision variables U and lower-level decision variables L. Note that, for a given upper-level solution U, solving the lower-level problems yields L. The routing variables are handled by the truck owner, while the loading instructions are determined by the BSS provider.


## 3.1. Assumptions and Notation

To facilitate the presentation of the essential ideas without loss of generality, the following basic assumptions are made in this paper:( 1 ) The outsourcing cost of BSS provider outsourcing relocation transportation to the third-party logistics companies is fixed.( 2 ) Each station is allowed to be visited at most once, as in the literature [26, 33], which means that the demand of a station is not allowed to be split. Similar to Forma et al. [15] and Ho and Szeto [17, 26], not all stations need to be visited, and initially balanced stations do not need to be visited.( 3 ) The depot is assumed to have no bikes.( 4 ) The load on the truck cannot be greater than the truck capacity.( 5 ) The truck starts from and returns to the depot and operates within a given time constraint.The notations used throughout this paper are listed in Table3 unless otherwise specified.Table 3
Notations for the bi-level model.
N Set of stations  N o Set of nodes, including the stations and the depot (indexed by 0)  i, j Indices of stations  s i 0 Initial inventory level at stationibefore the repositioning operation  s i Final inventory level at stationiafter the repositioning operation  f i ( s i ) Convex penalty function ofsi for stationi  s i I The target optimal inventory at stationi, argminsifi(si),i.e., the optimal number of bikes at station i  P Set of pick-up stations, i.e.,{i∈N∣si0>siI}  D Set of drop-off stations, i.e.,{i∈N∣si0<siI}  Q Truck capacity  c i Capacity of stationi  T Repositioning time  L Time required to load a bike from a station onto the truck  U Time required to unload a bike from the truck to a station  M A very large number  t i j Traveling time from stationito station j  x i j Binary variable that equals 1 if the truck travels directly from stationito stationj, and 0 otherwise  y i j Number of bikes on the truck when it travels directly from stationito stationj  y i P The loading quantities at stationi,  yiP≥0,∀i∈P,yiP=0,ifi∉P  y i D The unloading quantitiesat stationi,  yiD≥0,∀i∈D,yiD=0,ifi∉D  U The vector of the upper-level decision variables,U={yiD,yiP,si,∀i∈N}  L The vector of the lower-level decision variables,L={xij,yij,∀i∈N,j∈N,i≠j}


## 3.2. Penalty Function

The penalty function used in this paper was first proposed by Raviv et al. [33]. Similar to Ho and Szeto [17, 26] and Forma et al. [15], we used it to define customer dissatisfaction. For any stationi, the penalty functionfi(si)is a function of the final inventory of stationi after the operation, which represents the anticipated shortage of stationi. The penalty cost of stationi can be calculated by the corresponding penalty functionfi(si), and the value range ofsiis 0,…,ci.Ifsi is out of the domain 0,…,ci, thenfi(si)=∞.As described by Forma et al. [15], all stations of the BSS were divided into three categories, namely, (1) residential stations where the demand of customers to rent bikes is large in the morning and that to return bikes is large in the afternoon; (2) business stations where the demand of customers to return bikes is large in the morning and that to rent bikes is large in the afternoon;  (3)  ‘‘touristic’’ stations where the demand of customers to rent or return bikes does not change much throughout the day. The optimal inventory level at stationi can be solved by argminsifisi.


## 3.3. Bilevel Model

In order to demonstrate the relationship between the BSS provider and the truck owner, a bilevel programming model modified on the basis of Ho and Szeto [26] is proposed. The BSS operator’s decision variables are si, yiP, and yiD. The upper-level model aims to minimize the total penalties cost. The lower-level model aims to minimize the transportation cost with the truck capacity constraint and the repositioning time constraint.


### 3.3.1. The Upper-Level Model

The upper-level model aims to determine the optimalU at each station to minimize the total penalty costs and can be formulated as follows:(1)min∑i∈Nfisi(2)s.t.si=si0-yiP-yiD,∀i∈N(3)∑j∈No,j≠iyij-∑j∈No,j≠iyji=yiP-yiD,i∈N(4)yij≤Qxij,∀i,j∈No,i≠j(5)∑i∈PyiP-∑i∈DyiD=0(6)yiP≤si0,i∈P(7)yiD≤ci-si0,i∈D(8)yiP≥0,∀i∈P(9)yiD≥0,∀i∈N(10)si≥0,∀i∈N.The upper-level objective function (1) minimizes the total penalty costs generated by all stations. Constraint (2) defines the final inventory at stationiafter the operation as the initial inventory at stationiminus the loading quantities or plus the unloading quantities at stationi. Constraint (3) requires that the loading or unloading quantities at stationiare equal to the difference between the loads of the truck before and after the stationi. Constraint (4) specifies the upper bound of the truck’s load. Constraint (5) imposes that all bikes loaded from pick-up stations by the truck are finally unloaded at drop-off stations. Constraint (6) requires that the loading quantity at pick-up station should not exceed the initial inventory at that station. Constraint (7) requires that the unloading quantity at drop-off station should not exceed the residual capacity of that station. Constraints (8)-(10) make sure the variables are nonnegative. Decision variables xij can be solved from the lower-level model, while constraints (3) and (5) describe the relationship between the upper-level decision variables U and lower-level decision variables L.


### 3.3.2. The Lower-Level Model

The lower-level model aims to determine the optimalL to minimize the transportation cost and can be formulated as follows:(11)min∑i∈N0∑j∈N0,i≠jtijxij(12)s.t.∑j∈Ny0j=∑j∈Nyj0=0(13)∑j∈No,j≠ixij≤1,i∈N(14)∑j∈No,j≠ixij-∑j∈No,j≠ixji=0,i∈N(15)∑i∈NLyiP+UyiD+∑i,j∈N0,i≠jtijxij≤T(16)qj≥qi+1-M1-xij,∀i,j∈No,i≠j(17)qi≥0,∀i∈N0(18)xij∈0,1,∀i,j∈N0,i≠j(19)yij≥0,∀i,j∈N0,i≠j.The lower-level objective function (11) minimizes the total travel time. Constraint (12) forces the fact that the truck leaves and returns to the depot empty. Constraint (13) imposes that every station can only to be visited once. Constraint (14) guarantees that it is impossible for the truck to stay at any station. Constraint (15) imposes that the total handling and travel time of the truck should not exceed the given repositioning time. Constraint (16) eliminates the subtours [38]. Constraint (17) ensures that the auxiliary variable is nonnegative. Constraint (18) declares that xij is a binary variable. Constraint (19) ensures that the load of the truck is nonnegative.Constraint (15) describes the relationship between the upper-level decision variables U and lower-level decision variables L. Note that, for a given upper-level solution U, solving the lower-level problems yields L. The routing variables are handled by the truck owner, while the loading instructions are determined by the BSS provider.


## 3.3.1. The Upper-Level Model

The upper-level model aims to determine the optimalU at each station to minimize the total penalty costs and can be formulated as follows:(1)min∑i∈Nfisi(2)s.t.si=si0-yiP-yiD,∀i∈N(3)∑j∈No,j≠iyij-∑j∈No,j≠iyji=yiP-yiD,i∈N(4)yij≤Qxij,∀i,j∈No,i≠j(5)∑i∈PyiP-∑i∈DyiD=0(6)yiP≤si0,i∈P(7)yiD≤ci-si0,i∈D(8)yiP≥0,∀i∈P(9)yiD≥0,∀i∈N(10)si≥0,∀i∈N.The upper-level objective function (1) minimizes the total penalty costs generated by all stations. Constraint (2) defines the final inventory at stationiafter the operation as the initial inventory at stationiminus the loading quantities or plus the unloading quantities at stationi. Constraint (3) requires that the loading or unloading quantities at stationiare equal to the difference between the loads of the truck before and after the stationi. Constraint (4) specifies the upper bound of the truck’s load. Constraint (5) imposes that all bikes loaded from pick-up stations by the truck are finally unloaded at drop-off stations. Constraint (6) requires that the loading quantity at pick-up station should not exceed the initial inventory at that station. Constraint (7) requires that the unloading quantity at drop-off station should not exceed the residual capacity of that station. Constraints (8)-(10) make sure the variables are nonnegative. Decision variables xij can be solved from the lower-level model, while constraints (3) and (5) describe the relationship between the upper-level decision variables U and lower-level decision variables L.


## 3.3.2. The Lower-Level Model

The lower-level model aims to determine the optimalL to minimize the transportation cost and can be formulated as follows:(11)min∑i∈N0∑j∈N0,i≠jtijxij(12)s.t.∑j∈Ny0j=∑j∈Nyj0=0(13)∑j∈No,j≠ixij≤1,i∈N(14)∑j∈No,j≠ixij-∑j∈No,j≠ixji=0,i∈N(15)∑i∈NLyiP+UyiD+∑i,j∈N0,i≠jtijxij≤T(16)qj≥qi+1-M1-xij,∀i,j∈No,i≠j(17)qi≥0,∀i∈N0(18)xij∈0,1,∀i,j∈N0,i≠j(19)yij≥0,∀i,j∈N0,i≠j.The lower-level objective function (11) minimizes the total travel time. Constraint (12) forces the fact that the truck leaves and returns to the depot empty. Constraint (13) imposes that every station can only to be visited once. Constraint (14) guarantees that it is impossible for the truck to stay at any station. Constraint (15) imposes that the total handling and travel time of the truck should not exceed the given repositioning time. Constraint (16) eliminates the subtours [38]. Constraint (17) ensures that the auxiliary variable is nonnegative. Constraint (18) declares that xij is a binary variable. Constraint (19) ensures that the load of the truck is nonnegative.Constraint (15) describes the relationship between the upper-level decision variables U and lower-level decision variables L. Note that, for a given upper-level solution U, solving the lower-level problems yields L. The routing variables are handled by the truck owner, while the loading instructions are determined by the BSS provider.


## 4. Solution Method

As shown in Table2, several papers have been dedicated to the application of heuristic algorithms in solving the BRP arising in a single-level programming model. As pointed out by Moore and Bard [39], the mixed integer bilevel programming models are NP-hard. In this section, we focus on the applications of ILS and TS for solving the bilevel BRP. The proposed general heuristic is mostly based on the ILS framework. Before describing the solution method, a brief outline of this metaheuristic is first provided. Then, we describe the solution representation, the initialization procedure, the neighborhood operators, and the perturbation operators that are used in all heuristic algorithms.


### 4.1. Basic Description of ILS

ILS was introduced by Lourenço et al. in 2001 and iteratively alternates between local search (intensification, see Sections4.4.1 and 4.4.2) and perturbation (diversification, see Sections 4.4.3 and 4.4.4) mechanisms with a view to finding high-quality solutions until a given stopping criterion is satisfied [40]. The local search consists of a set of neighborhood moves performed on solutions, to yield a local optimums∗. The main drawback of local descent is that it may get trapped in local optima that are significantly worse than the global optimum. Much like simulated annealing, ILS escapes from local optima by applying perturbations to drive the search towards a new region of the solution space. The basic ILS structure is given by Lourenço et al. [41], as shown in Algorithm 1. Starting with a local minimum s∗, we apply aPerturbation, leading to a solution s′. After applyingLocalSearch, we find a new local minimum s∗′ that may be better than s∗. To implement an ILS algorithm, four procedures should be specified: (i)GenerateInitialSolution, where an initial feasible solution is constructed; (ii)LocalSearch, which improves the solution initially obtained; (iii)Perturbation, where a new starter point is generated through a perturbation of the solution returned by theLocalSearch; (iv)Acceptance Criterion, that determines from which solution the search should be accepted as the new current solution. The number of iterations performed depends on the predefined termination criterion.Algorithm 1:ILS.
s0≔ GenerateInitialSolution s ≔ LocalSearch(s0) while  termination criterion is not met  do s′≔ Perturbation(s,history); s′ ≔ LocalSearch(s′); s≔ AcceptanceCriterion(s, s∗′, history); end while Return  s∗


### 4.2. Solution Representation

A solution to the bilevel BRP with a single truck considered in this work can be represented as a sequence of visits to stations, starting and ending at the depot, along with the number of bikes loaded or unloaded at each visit.A solutionω is represented by ω= (x,d,q), wherex is a vector to store the truck route;d is a vector to store the loading/unloading quantity at the station visited by the truck, andq is a vector to store the truck load during the route.Letn be the total number of stations (excluding the starting and ending depot) visited by the truck. The stations visited by the truck are denoted asxh, xh∈N,h∈[1,n],and the repositioning truck is assumed to start and end repositioning at the depot (denoted byx0,xn+1,andx0=xn+1=0).The elements ofd are defined as dh=yxhP-yxhD,∀h∈[0,n+1]. A positive integer represents a loading operation, while a negative integer indicates an unloading operation. According to the constraints (12),d0=dn+1=0.The last elements ofqare defined as qh=∑k=0hdk,∀h∈[0,n+1], where dh denotes the loading/unloading quantity at stationxh. In other words,qh gives the value of the load on the truck after visiting stationxh of the route. The truck leaves the depot with no load; then, according to Hernández-Pérez and Salazar-González [32], a route is feasible if the following inequality and equality are satisfied:(20)maxk=0n+1qk-mink=0n+1qk≤Qmink=0n+1qk=0.A feasible solution for the repositioning operation performed by the truck with the capacity of 10 bikes is shown in Figure2. In this example,n = 8, and the truck starts and ends the repositioning at depot 0 empty.Figure 2
Solution representation.


### 4.3. Initialization

In this paper, the initial feasible solution is quickly generated by a construction heuristic algorithm based on the penalty cost. The construction heuristic algorithm was first proposed by Ho and Szeto [26] to solve the single vehicle BRP, and then Ho and Szeto [17] used it to solve the multiple vehicle BRP.The approximate flow of the construction heuristic algorithm can be divided into two steps. Setp one, sort all stations in descending order by dividing them into pick-up and drop-off stations. The ranking criteria arefi(si0)-fi(siQ) rather than fi(si0)-fi(siI) used by Ho and Szeto [26], where siQ=si0-minQ,si0-siI,i∈P;si0-max(-Q,si0-siI),i∈D. The main reason for improvement is that if the loading or unloading quantities at stationi exceedQ, there is no way to meet the excess demand. Step two, assign insert in turn; first assign the pick-up stations; if no pick-up station can be inserted into the current route without exceeding the capacity limits, then drop-off stations are inserted. The other detailed algorithm flow can be found in Ho and Szeto [26].


### 4.4. Local Search and Perturbation Mechanism

The neighborhood structure of local search is one of the most important components that determine the size of the search space and the quality of the final solution. The perturbation procedure is a diversification scheme that helps the search process to escape from local optima. In this paper, the local search in our proposed ILS is executed by a randomized variable neighborhood descent (RVND) procedure. RVND consists of systematically examining different types of neighborhoods in a random ordering of the neighborhoods. All possible moves of each neighborhood are examined and the search resumes from the best improving neighbor. If a neighborhood is not capable of finding an improved solution, the procedure then selects another one at random. The procedure ends when all neighborhoods fail to refine the current solution.


#### 4.4.1. Local Neighborhood Structures for Lower-Level Model

Focusing on the problem of determining the optimal truck route when the pickup and drop-off quantities are given by the upper-level model, the following neighborhood structures for the lower-level model were implemented.Reinsertion (σ): one subsequence σ is removed and reinserted in another position on the route. The size of σ was limited to 1–3, thus resulting in three different neighborhoods.Exchange (σ): one subsequence σ is interchanged with another subsequence σ. The size of σ was restricted to 1 and 2, thus leading to two distinct neighborhoods (disregarding symmetries).2-opt: Two nonadjacent arcs are removed from the route and then two new ones are inserted into the route. In other words, a subsequence of the tour is reversed.


#### 4.4.2. Local Neighborhood Structures for Bilevel Model

Insert: The move inserts a station which is not in the current solution into the current solution and then performs the repairing heuristic algorithm (see Section4.6) to the new solution feasible with theBest Choose Criterion (see Section 4.6).Swap: The move replaces stationi in the solution with another stationj which is not in the current solution. The stationi and stationj must be of the same type, only the station number is exchanged, and the loading and unloading quantity are not changed. This exchange is feasible if the inventory level of stationj satisfies constraints (6) and (7), and the resultant traveling time does not violate constraint (15).Delete: The move deletes a station from the current solution and then performs the repairing heuristic algorithm (see Section4.6) to make the neighbor solution feasible with theBest Choose Criterion (see Section 4.6).


#### 4.4.3. Perturbation Mechanism for Lower-Level Model

The perturbation procedure randomly selects one of the following mechanisms:Rand Reinsertion(σ): Reinsertion (σ) moves are applied at random and we limited σ to be at most 3.Rand Exchange: Exchange(σ) moves are applied at random. In our case, we limited σ to be at most 2.


#### 4.4.4. Perturbation Mechanism for Bilevel Model

Rand Insert: The move inserts a station which is not in the current solution into the current solution and then performs the repairing heuristic algorithm (see Section4.6) to the new solution feasible with theRandom Choose Criterion (see Section 4.6).Rand Delete: The move deletes a station from the current solution and then performs the repairing heuristic algorithm (see Section4.6) to make the neighbor solution feasible with theRandom Choose Criterion (see Section 4.6).


### 4.5. Solution Algorithm for Bilevel Model

In this study, the bilevel BRP is solved as follows. With a given inventory level U, the truck route L can be found by solving the lower-level model using ILS (see Section4.5.1). With the truck route L solved from the lower-level model, the objective function of the upper-level model can be evaluated and a new set of U can be generated. The above procedure will then be repeated until it is converged to the optimal solution of the bilevel BRP model. We develop ILS to solve the lower-level model and propose ILS and TS to solve the bilevel model. In this section, the ILS procedure for the lower-level model (see Section 4.5.1) will first be proposed, and then the ILS (see Section 4.5.2) and TS (see Section 4.5.3) procedure for the bilevel model will be described in turn. The local search and perturbation neighborhoods performed in ILS and TS are summarized in Table 4.Table 4
Local search and perturbation mechanism.
Algorithm Neighborhood ILS_LM ILS_BL TS_BL Local search LS_LM 2-opt LS_BL Insert LS_BL Insert Relocation Swap Swap Exchange Remove Remove  Perturbationmechanism Perb_LM Rand Exchange Perb_BL Rand Insert Rand Relocation Rand Delete


#### 4.5.1. ILS for Lower Level

Algorithm2 illustrates the flowchart of ILS for the low-level model. The ILS heuristic performs a procedure that iteratively calls the local search and perturbation moves until MaxConsIter_LM consecutive iterations without improvement have been executed. The best solution of the current iteration is always the one chosen to be perturbed. With the initial solutions given by the upper-level model, the local search which consists of 2-opt, Relocation, and Exchange neighborhoods, as LS_LM, is executed using a RVND procedure (see Section 4.4.1), whereas the perturbation procedure, as Perb_LM, performs Random Relocation and Exchange neighborhoods (see Section 4.4.2). In this procedure, letr represent the current route,r∗ represent the best route that can be searched for at present,RC(r) represent the travel cost of current router, ConsIter_LM represent the current consecutive iterations without improvement, and MaxConsIter_LM represent the maximum consecutive iterations without improvement.Algorithm 2:ILS for lower-level model.
Function ILS_LM(s, LS_LM, Perb_LM, MaxConsIter_LM) r≔s; r∗≔s; ConsIter _LM=0; while  ConsIter_LM<=MaxConsIter_LM  do r′≔RVND(r, LS_LW); if  RC(r′)<RC(r∗),then r∗≔ r′; ConsIter _LM≔0; else ConsIter _LM≔ ConsIter _LM +1; end r≔Perturbation(r∗, Perb_LW); end Return  r∗


#### 4.5.2. ILS for Bilevel

Algorithm3 provides an outline of the ILS for the bilevel model. According to Algorithm 3, after constructing the initial solutions0 by applying the initialization procedure (see Section 4.3), the procedures RVND, ILS_LM (see Section 4.5.1), and Perturbation (see Section 4.4.4) are called repeatedly until the predefined termination criterion is satisfied. The best current solution is chosen to be perturbed by perturbation moves. In this procedure, lets represent the current solution for the bilevel model,s∗ represent the best solution that can be searched for at present,TC(s) represent the total penalties cost of current solutions, ConsIter_BL represent the current consecutive iterations without improvement, and MaxConsIter_BL represent the maximum consecutive iterations without improvement.Algorithm 3:ILS for bi-level model.
Function  ILS_BL(LS_LM, Perb_LM, LS_BL, Perb_BL, MaxConsIter_BL) s 0≔Initialization; s ∗≔s0; ConsIter_BL≔0; while  ConsIter_BL<=MaxConsIter_BL  do s′ ≔RVND(s0, LS_BL); s≔ILS_LM(s′,LS_LM, Perb_LM); ifTC(s )<TC(s∗),then s∗≔ s; ConsIter _BL≔0; else ConsIter _BL≔ ConsIter _BL +1; end s0≔Perturbation(s∗, Perb_BL); end Return  s∗


#### 4.5.3. TS for Bilevel

The TS algorithm was introduced by Glover in 1986 and is an extension to the local search algorithm. It has been successfully applied to tackle a wide variety of routing problems, such as the traveling salesman problem, capacity vehicle routing problems, pickup and delivery problem, BRP [23, 26], and other routing problems. At each iteration, TS performs an exploration of the solution space by moving from a solutionw identified to the best solution in a subset of the neighborhoodN(w) ofw. To avoid being trapped in the local optima and prevent cycling on the same solution, TS uses an explicit memory structure (i.e., tabu list) to record the search track. If the move is in the tabu list, the move is forbidden. Sometimes tabu lists may prevent excellent solutions being chosen. The aspiration criterion provides the opportunity to overcome this obstacle. That is, the aspiration criterion enables superior neighboring solutions to be chosen, even though the solution is obtained by a banned move in the tabu list. We utilize two independent tabu lists, TABU1 and TABU2, as follows for the removal and insertion moves, respectively.TABU1: for a removal move, the tabu list records the station number which is removed from the current solution. When generating a neighborhood, if the station selected to be removed from the current solution exists in the tabu list, the move is marked as tabu. The length of the tabu list is denoted asl1.TABU2: for an insertion move, the tabu list records the station number which is inserted into the current solution. When generating a neighborhood, if the station selected to be inserted into the current solution exists in the tabu list, the move is marked as tabu. The length of the tabu list is denoted asl2.For a swap move, if stationi in the current solution is swapped with another stationj which is not in the current solution, it means that this swap move integrated removing stationiand inserting stationj. When generating a neighborhood by a swapping move, if stationi exists in the tabu list TABU1 or stationj exists in the tabu list TABU2, the swap move is marked as tabu.In fact, after generating an initial solution, the following two steps are repeated until the stopping criterion is met. The stopping criterion of TS is achieving a predetermined number of iterations with no improvement (MaxConsIter_TS).Step 1.
Generating a new solution: a temporary solutions′ is generated from all feasible neighborhood operators (LS_BL, see Section 4.4.2) of the current solutions. Then, the optimal low-level solution s″ is obtained throughILS_lower (see Section 4.5.1).Step 2.
Evaluating the new solution: in this step, the objective value of the new solution is compared with that of the best solution found. If the new solution improves the best solution, it is accepted and the best solution will be updated; otherwise, the new solution is accepted only if its move is not in the tabu list.In this TS_BL procedure, lets represent the current solution to the bilevel model,s∗ represent the best solution that can be searched for at present, tabu represent the two tabu lists,N(s) represent the collection of the neighborhood solutions generated by the current solutions, TC(s) represent the total penalties cost of current solutions,TC(s∗) represent the minimum total penalties cost so far, ConsIter_TS represent the current consecutive iterations without improvement, and MaxConsIter_TS represent the maximum consecutive iterations without improvement. The detailed pseudocode of TS for the bilevel model is presented  in Algorithm 4.Algorithm 4:TS for bi-level model.
Function  TS_BL(LS_LW, Perb_LW, LS_BL) s 0≔Initialization; Tuba=; s ∗≔s0; s≔s0; ConsIter_TS≔0; while  ConsIter_TS<=MaxConsIter_TS  do N(s) ≔ LocalSearch(s, LS_BL);(N(s) consists of all feasible neighbor solutions obtained by applying the three local search moves). Select a solutions′ in the set N(s) that minimizes TS(s), and the solution is either a non-tabu solution or a solution better than the best one so far (i.e.,TC(s′)<TC(s∗)), and update the tabu. s≔ILS_LM(s′, LS_LM, Perb_LM); ifTC(s)<TC(s∗),then s∗≔ s; ConsIter_TS≔0; else ConsIter_TS≔ ConsIter_TS+1; end end Return  s∗


### 4.6. Repairing Infeasible Solutions

Sometimes the neighbor solution may no longer remain feasible with reference to the loading and capacity constraints after implementing the Insert and Delete neighborhood moves. Hence, the solution can be repaired by adjusting the pick-up/drop-off quantities at one of the stations on the new route visited by the truck.Letω= (x,d,q) be the current solution, set x=(x0,x1,…,xn,xn+1),x0=xn+1=0, let S(x) be the set of stations in the route, S(x)={x1,…,xn}. Let σ=(σ1,…,σnσ) be a subsequence of a routex of solutionw, let djσ be the loading and unloading quantities at the station σj, let qjσ be the cumulative load at station σj, and then qjσ=∑i=1jdiσ.Letqminσ=min{0,q1σ,…,qnσσ},qmaxσ=max{0,q1σ,…,qnσσ}; that is to say, qminσis the minimum cumulative load and qmaxσ is the maximum cumulative load. According to Bulhões et al. [20], we can obtain the minimum and maximum flow of bikes allowed to enter the subsequence σ so as to ensure route feasibility. The minimum flow of bikes allowed to enter so as to ensure feasibility is lminσ=-qminσ; the maximum flow of bikes allowed to enter so as to ensure feasibility is lmaxσ=Q-qmaxσ.Delete: If we delete the station xh from the current solution, the loading and unloading quantities of the station xh in the current solution are dh, set ΔD=dh. After performing the delete move, let the new route of the new solution be xD=(x0,x1,…,xh-1,xh,…,xn+1)=(x0D,x1D,…xnDD,xnD+1D),nD=n-1.LetS(xD)=x1D,…,xnDD be the set of stations in the route. According to constraints (3) and (5), the total quantities of the route xD reduce to ΔD; in order to maintain the quantity balance, we need to increase ΔDat another station in S(xD). First, we put all feasible stations xkD∈S(xD),k∈[1,nD] that can adjust the quantity to obtain a feasible solution by increasing ΔD into the set FD,FD⊆S(xD). Then, we choose one station from FD according to the predetermined rules. If we choose station x∗D,x∗D∈FD, then set dx∗D=dx∗D+ΔD. The characteristics of the stations in SF must first be analyzed, and then the rules to choose the station to generate the feasible neighbor solution will be proposed.∀ k ∈ [ 1 , n D ], xD=(x0D,x1D,…xnDD,xnD+1D),nD=n, let σD=(xkD,…,xnDD) be a subsequence of a route xD. The station xkDcan be put into FD only if all the following conditions are satisfied:(21)ifxkD∈P,anddxkD+ΔD>0;dxkD+ΔD≤sxkD0;qksum+ΔD∈lminσD,lmaxσD.ifxkD∈D,anddxkD+ΔD<0;dxkD+ΔD≤cxkD-sxkD0;qksum+ΔD∈lminσD,lmaxσD.There are two rules to choose the station by adjusting the quantities to generate the feasible new solution:Best Choose Criterion andRandom Choose Criterion. The Best Choose Criterion is to choose the stationi from set FD with the least increasing penalty cost. Stationi is determined from argmini∈FD{f(siI-ΔD)-f(siI)}. The Random Choose Criterion is to randomly choose the stationi from set FD.Insert: If we insert a station m∈N∖S(x) into the current route between xj and xj+1, the loading and unloading quantities of the new stationm are dm, and let ΔI=dm. After performing the insert move, set the new route of the new solution as xI=(x0,x1,…xj,h,xj+1,…,xn,xn+1)=(x0I,x1I,…,xnII,xnI+1I),nI=n+1.SetS(xI)={x1I,…,xnII} as the set of stations in the new route xI. According to constraints (3) and (5), the total quantities of the route xI increase ΔI, and in order to maintain the quantity balance, we need to decrease ΔD at another station in S(xI). First, we put all feasible stations xkI∈S(xI),k∈[1,j]∪[j+1,nI] in S(xI), which can adjust the quantity to obtain a feasible solution through decreasing ΔI into the set FI,FI⊆S(xD). Then, we choose one station from FI according to theBest/Random Choose Criterion. The Best Choose Criterion is to choose the stationi from set FI with the least increasing penalty cost. The stationi is determined from argmini∈FI{f(siI+ΔI)-f(siI)}. The Random Choose Criterion is to randomly choose the stationi from set FI.If we choose stationx∗I,x∗I∈FI, then set dx∗I=dx∗I-ΔI. Obviously, inserting a station into the route will need more time for traveling. In addition to satisfying this truck capacity constraint, a feasible neighbor solution after performing the insert move also needs to obey the time constraint of the truck.∀ k ∈ [ 1 , j ] ∪ [ j + 2 , n I ], xI=(x0I,x1I,…,xnII,xnI+1I),nI=n+1, let σI=(xkI,…,xnII) be a subsequence of a route xI.The station xkIcan be put into FI only if all the following conditions are satisfied:(22)ifxkI∈P,anddxkI-ΔI>0;dxkI-ΔI≤sxkD0;qksum-ΔI∈lminσI,lmaxσITimexD+∑j∈Sx,j≠kd+dxkD-ΔD≤T.ifxkI∈D,anddxkI-ΔI<0;dxkI-ΔI≤cxkI-sxkI0;qksum-ΔI∈lminσI,lmaxσITimexD+∑j∈Sx,j≠kd+dxkD-ΔD≤T.


## 4.1. Basic Description of ILS

ILS was introduced by Lourenço et al. in 2001 and iteratively alternates between local search (intensification, see Sections4.4.1 and 4.4.2) and perturbation (diversification, see Sections 4.4.3 and 4.4.4) mechanisms with a view to finding high-quality solutions until a given stopping criterion is satisfied [40]. The local search consists of a set of neighborhood moves performed on solutions, to yield a local optimums∗. The main drawback of local descent is that it may get trapped in local optima that are significantly worse than the global optimum. Much like simulated annealing, ILS escapes from local optima by applying perturbations to drive the search towards a new region of the solution space. The basic ILS structure is given by Lourenço et al. [41], as shown in Algorithm 1. Starting with a local minimum s∗, we apply aPerturbation, leading to a solution s′. After applyingLocalSearch, we find a new local minimum s∗′ that may be better than s∗. To implement an ILS algorithm, four procedures should be specified: (i)GenerateInitialSolution, where an initial feasible solution is constructed; (ii)LocalSearch, which improves the solution initially obtained; (iii)Perturbation, where a new starter point is generated through a perturbation of the solution returned by theLocalSearch; (iv)Acceptance Criterion, that determines from which solution the search should be accepted as the new current solution. The number of iterations performed depends on the predefined termination criterion.Algorithm 1:ILS.
s0≔ GenerateInitialSolution s ≔ LocalSearch(s0) while  termination criterion is not met  do s′≔ Perturbation(s,history); s′ ≔ LocalSearch(s′); s≔ AcceptanceCriterion(s, s∗′, history); end while Return  s∗


## 4.2. Solution Representation

A solution to the bilevel BRP with a single truck considered in this work can be represented as a sequence of visits to stations, starting and ending at the depot, along with the number of bikes loaded or unloaded at each visit.A solutionω is represented by ω= (x,d,q), wherex is a vector to store the truck route;d is a vector to store the loading/unloading quantity at the station visited by the truck, andq is a vector to store the truck load during the route.Letn be the total number of stations (excluding the starting and ending depot) visited by the truck. The stations visited by the truck are denoted asxh, xh∈N,h∈[1,n],and the repositioning truck is assumed to start and end repositioning at the depot (denoted byx0,xn+1,andx0=xn+1=0).The elements ofd are defined as dh=yxhP-yxhD,∀h∈[0,n+1]. A positive integer represents a loading operation, while a negative integer indicates an unloading operation. According to the constraints (12),d0=dn+1=0.The last elements ofqare defined as qh=∑k=0hdk,∀h∈[0,n+1], where dh denotes the loading/unloading quantity at stationxh. In other words,qh gives the value of the load on the truck after visiting stationxh of the route. The truck leaves the depot with no load; then, according to Hernández-Pérez and Salazar-González [32], a route is feasible if the following inequality and equality are satisfied:(20)maxk=0n+1qk-mink=0n+1qk≤Qmink=0n+1qk=0.A feasible solution for the repositioning operation performed by the truck with the capacity of 10 bikes is shown in Figure2. In this example,n = 8, and the truck starts and ends the repositioning at depot 0 empty.Figure 2
Solution representation.


## 4.3. Initialization

In this paper, the initial feasible solution is quickly generated by a construction heuristic algorithm based on the penalty cost. The construction heuristic algorithm was first proposed by Ho and Szeto [26] to solve the single vehicle BRP, and then Ho and Szeto [17] used it to solve the multiple vehicle BRP.The approximate flow of the construction heuristic algorithm can be divided into two steps. Setp one, sort all stations in descending order by dividing them into pick-up and drop-off stations. The ranking criteria arefi(si0)-fi(siQ) rather than fi(si0)-fi(siI) used by Ho and Szeto [26], where siQ=si0-minQ,si0-siI,i∈P;si0-max(-Q,si0-siI),i∈D. The main reason for improvement is that if the loading or unloading quantities at stationi exceedQ, there is no way to meet the excess demand. Step two, assign insert in turn; first assign the pick-up stations; if no pick-up station can be inserted into the current route without exceeding the capacity limits, then drop-off stations are inserted. The other detailed algorithm flow can be found in Ho and Szeto [26].


## 4.4. Local Search and Perturbation Mechanism

The neighborhood structure of local search is one of the most important components that determine the size of the search space and the quality of the final solution. The perturbation procedure is a diversification scheme that helps the search process to escape from local optima. In this paper, the local search in our proposed ILS is executed by a randomized variable neighborhood descent (RVND) procedure. RVND consists of systematically examining different types of neighborhoods in a random ordering of the neighborhoods. All possible moves of each neighborhood are examined and the search resumes from the best improving neighbor. If a neighborhood is not capable of finding an improved solution, the procedure then selects another one at random. The procedure ends when all neighborhoods fail to refine the current solution.


### 4.4.1. Local Neighborhood Structures for Lower-Level Model

Focusing on the problem of determining the optimal truck route when the pickup and drop-off quantities are given by the upper-level model, the following neighborhood structures for the lower-level model were implemented.Reinsertion (σ): one subsequence σ is removed and reinserted in another position on the route. The size of σ was limited to 1–3, thus resulting in three different neighborhoods.Exchange (σ): one subsequence σ is interchanged with another subsequence σ. The size of σ was restricted to 1 and 2, thus leading to two distinct neighborhoods (disregarding symmetries).2-opt: Two nonadjacent arcs are removed from the route and then two new ones are inserted into the route. In other words, a subsequence of the tour is reversed.


### 4.4.2. Local Neighborhood Structures for Bilevel Model

Insert: The move inserts a station which is not in the current solution into the current solution and then performs the repairing heuristic algorithm (see Section4.6) to the new solution feasible with theBest Choose Criterion (see Section 4.6).Swap: The move replaces stationi in the solution with another stationj which is not in the current solution. The stationi and stationj must be of the same type, only the station number is exchanged, and the loading and unloading quantity are not changed. This exchange is feasible if the inventory level of stationj satisfies constraints (6) and (7), and the resultant traveling time does not violate constraint (15).Delete: The move deletes a station from the current solution and then performs the repairing heuristic algorithm (see Section4.6) to make the neighbor solution feasible with theBest Choose Criterion (see Section 4.6).


### 4.4.3. Perturbation Mechanism for Lower-Level Model

The perturbation procedure randomly selects one of the following mechanisms:Rand Reinsertion(σ): Reinsertion (σ) moves are applied at random and we limited σ to be at most 3.Rand Exchange: Exchange(σ) moves are applied at random. In our case, we limited σ to be at most 2.


### 4.4.4. Perturbation Mechanism for Bilevel Model

Rand Insert: The move inserts a station which is not in the current solution into the current solution and then performs the repairing heuristic algorithm (see Section4.6) to the new solution feasible with theRandom Choose Criterion (see Section 4.6).Rand Delete: The move deletes a station from the current solution and then performs the repairing heuristic algorithm (see Section4.6) to make the neighbor solution feasible with theRandom Choose Criterion (see Section 4.6).


## 4.4.1. Local Neighborhood Structures for Lower-Level Model

Focusing on the problem of determining the optimal truck route when the pickup and drop-off quantities are given by the upper-level model, the following neighborhood structures for the lower-level model were implemented.Reinsertion (σ): one subsequence σ is removed and reinserted in another position on the route. The size of σ was limited to 1–3, thus resulting in three different neighborhoods.Exchange (σ): one subsequence σ is interchanged with another subsequence σ. The size of σ was restricted to 1 and 2, thus leading to two distinct neighborhoods (disregarding symmetries).2-opt: Two nonadjacent arcs are removed from the route and then two new ones are inserted into the route. In other words, a subsequence of the tour is reversed.


## 4.4.2. Local Neighborhood Structures for Bilevel Model

Insert: The move inserts a station which is not in the current solution into the current solution and then performs the repairing heuristic algorithm (see Section4.6) to the new solution feasible with theBest Choose Criterion (see Section 4.6).Swap: The move replaces stationi in the solution with another stationj which is not in the current solution. The stationi and stationj must be of the same type, only the station number is exchanged, and the loading and unloading quantity are not changed. This exchange is feasible if the inventory level of stationj satisfies constraints (6) and (7), and the resultant traveling time does not violate constraint (15).Delete: The move deletes a station from the current solution and then performs the repairing heuristic algorithm (see Section4.6) to make the neighbor solution feasible with theBest Choose Criterion (see Section 4.6).


## 4.4.3. Perturbation Mechanism for Lower-Level Model

The perturbation procedure randomly selects one of the following mechanisms:Rand Reinsertion(σ): Reinsertion (σ) moves are applied at random and we limited σ to be at most 3.Rand Exchange: Exchange(σ) moves are applied at random. In our case, we limited σ to be at most 2.


## 4.4.4. Perturbation Mechanism for Bilevel Model

Rand Insert: The move inserts a station which is not in the current solution into the current solution and then performs the repairing heuristic algorithm (see Section4.6) to the new solution feasible with theRandom Choose Criterion (see Section 4.6).Rand Delete: The move deletes a station from the current solution and then performs the repairing heuristic algorithm (see Section4.6) to make the neighbor solution feasible with theRandom Choose Criterion (see Section 4.6).


## 4.5. Solution Algorithm for Bilevel Model

In this study, the bilevel BRP is solved as follows. With a given inventory level U, the truck route L can be found by solving the lower-level model using ILS (see Section4.5.1). With the truck route L solved from the lower-level model, the objective function of the upper-level model can be evaluated and a new set of U can be generated. The above procedure will then be repeated until it is converged to the optimal solution of the bilevel BRP model. We develop ILS to solve the lower-level model and propose ILS and TS to solve the bilevel model. In this section, the ILS procedure for the lower-level model (see Section 4.5.1) will first be proposed, and then the ILS (see Section 4.5.2) and TS (see Section 4.5.3) procedure for the bilevel model will be described in turn. The local search and perturbation neighborhoods performed in ILS and TS are summarized in Table 4.Table 4
Local search and perturbation mechanism.
Algorithm Neighborhood ILS_LM ILS_BL TS_BL Local search LS_LM 2-opt LS_BL Insert LS_BL Insert Relocation Swap Swap Exchange Remove Remove  Perturbationmechanism Perb_LM Rand Exchange Perb_BL Rand Insert Rand Relocation Rand Delete


### 4.5.1. ILS for Lower Level

Algorithm2 illustrates the flowchart of ILS for the low-level model. The ILS heuristic performs a procedure that iteratively calls the local search and perturbation moves until MaxConsIter_LM consecutive iterations without improvement have been executed. The best solution of the current iteration is always the one chosen to be perturbed. With the initial solutions given by the upper-level model, the local search which consists of 2-opt, Relocation, and Exchange neighborhoods, as LS_LM, is executed using a RVND procedure (see Section 4.4.1), whereas the perturbation procedure, as Perb_LM, performs Random Relocation and Exchange neighborhoods (see Section 4.4.2). In this procedure, letr represent the current route,r∗ represent the best route that can be searched for at present,RC(r) represent the travel cost of current router, ConsIter_LM represent the current consecutive iterations without improvement, and MaxConsIter_LM represent the maximum consecutive iterations without improvement.Algorithm 2:ILS for lower-level model.
Function ILS_LM(s, LS_LM, Perb_LM, MaxConsIter_LM) r≔s; r∗≔s; ConsIter _LM=0; while  ConsIter_LM<=MaxConsIter_LM  do r′≔RVND(r, LS_LW); if  RC(r′)<RC(r∗),then r∗≔ r′; ConsIter _LM≔0; else ConsIter _LM≔ ConsIter _LM +1; end r≔Perturbation(r∗, Perb_LW); end Return  r∗


### 4.5.2. ILS for Bilevel

Algorithm3 provides an outline of the ILS for the bilevel model. According to Algorithm 3, after constructing the initial solutions0 by applying the initialization procedure (see Section 4.3), the procedures RVND, ILS_LM (see Section 4.5.1), and Perturbation (see Section 4.4.4) are called repeatedly until the predefined termination criterion is satisfied. The best current solution is chosen to be perturbed by perturbation moves. In this procedure, lets represent the current solution for the bilevel model,s∗ represent the best solution that can be searched for at present,TC(s) represent the total penalties cost of current solutions, ConsIter_BL represent the current consecutive iterations without improvement, and MaxConsIter_BL represent the maximum consecutive iterations without improvement.Algorithm 3:ILS for bi-level model.
Function  ILS_BL(LS_LM, Perb_LM, LS_BL, Perb_BL, MaxConsIter_BL) s 0≔Initialization; s ∗≔s0; ConsIter_BL≔0; while  ConsIter_BL<=MaxConsIter_BL  do s′ ≔RVND(s0, LS_BL); s≔ILS_LM(s′,LS_LM, Perb_LM); ifTC(s )<TC(s∗),then s∗≔ s; ConsIter _BL≔0; else ConsIter _BL≔ ConsIter _BL +1; end s0≔Perturbation(s∗, Perb_BL); end Return  s∗


### 4.5.3. TS for Bilevel

The TS algorithm was introduced by Glover in 1986 and is an extension to the local search algorithm. It has been successfully applied to tackle a wide variety of routing problems, such as the traveling salesman problem, capacity vehicle routing problems, pickup and delivery problem, BRP [23, 26], and other routing problems. At each iteration, TS performs an exploration of the solution space by moving from a solutionw identified to the best solution in a subset of the neighborhoodN(w) ofw. To avoid being trapped in the local optima and prevent cycling on the same solution, TS uses an explicit memory structure (i.e., tabu list) to record the search track. If the move is in the tabu list, the move is forbidden. Sometimes tabu lists may prevent excellent solutions being chosen. The aspiration criterion provides the opportunity to overcome this obstacle. That is, the aspiration criterion enables superior neighboring solutions to be chosen, even though the solution is obtained by a banned move in the tabu list. We utilize two independent tabu lists, TABU1 and TABU2, as follows for the removal and insertion moves, respectively.TABU1: for a removal move, the tabu list records the station number which is removed from the current solution. When generating a neighborhood, if the station selected to be removed from the current solution exists in the tabu list, the move is marked as tabu. The length of the tabu list is denoted asl1.TABU2: for an insertion move, the tabu list records the station number which is inserted into the current solution. When generating a neighborhood, if the station selected to be inserted into the current solution exists in the tabu list, the move is marked as tabu. The length of the tabu list is denoted asl2.For a swap move, if stationi in the current solution is swapped with another stationj which is not in the current solution, it means that this swap move integrated removing stationiand inserting stationj. When generating a neighborhood by a swapping move, if stationi exists in the tabu list TABU1 or stationj exists in the tabu list TABU2, the swap move is marked as tabu.In fact, after generating an initial solution, the following two steps are repeated until the stopping criterion is met. The stopping criterion of TS is achieving a predetermined number of iterations with no improvement (MaxConsIter_TS).Step 1.
Generating a new solution: a temporary solutions′ is generated from all feasible neighborhood operators (LS_BL, see Section 4.4.2) of the current solutions. Then, the optimal low-level solution s″ is obtained throughILS_lower (see Section 4.5.1).Step 2.
Evaluating the new solution: in this step, the objective value of the new solution is compared with that of the best solution found. If the new solution improves the best solution, it is accepted and the best solution will be updated; otherwise, the new solution is accepted only if its move is not in the tabu list.In this TS_BL procedure, lets represent the current solution to the bilevel model,s∗ represent the best solution that can be searched for at present, tabu represent the two tabu lists,N(s) represent the collection of the neighborhood solutions generated by the current solutions, TC(s) represent the total penalties cost of current solutions,TC(s∗) represent the minimum total penalties cost so far, ConsIter_TS represent the current consecutive iterations without improvement, and MaxConsIter_TS represent the maximum consecutive iterations without improvement. The detailed pseudocode of TS for the bilevel model is presented  in Algorithm 4.Algorithm 4:TS for bi-level model.
Function  TS_BL(LS_LW, Perb_LW, LS_BL) s 0≔Initialization; Tuba=; s ∗≔s0; s≔s0; ConsIter_TS≔0; while  ConsIter_TS<=MaxConsIter_TS  do N(s) ≔ LocalSearch(s, LS_BL);(N(s) consists of all feasible neighbor solutions obtained by applying the three local search moves). Select a solutions′ in the set N(s) that minimizes TS(s), and the solution is either a non-tabu solution or a solution better than the best one so far (i.e.,TC(s′)<TC(s∗)), and update the tabu. s≔ILS_LM(s′, LS_LM, Perb_LM); ifTC(s)<TC(s∗),then s∗≔ s; ConsIter_TS≔0; else ConsIter_TS≔ ConsIter_TS+1; end end Return  s∗


## 4.5.1. ILS for Lower Level

Algorithm2 illustrates the flowchart of ILS for the low-level model. The ILS heuristic performs a procedure that iteratively calls the local search and perturbation moves until MaxConsIter_LM consecutive iterations without improvement have been executed. The best solution of the current iteration is always the one chosen to be perturbed. With the initial solutions given by the upper-level model, the local search which consists of 2-opt, Relocation, and Exchange neighborhoods, as LS_LM, is executed using a RVND procedure (see Section 4.4.1), whereas the perturbation procedure, as Perb_LM, performs Random Relocation and Exchange neighborhoods (see Section 4.4.2). In this procedure, letr represent the current route,r∗ represent the best route that can be searched for at present,RC(r) represent the travel cost of current router, ConsIter_LM represent the current consecutive iterations without improvement, and MaxConsIter_LM represent the maximum consecutive iterations without improvement.Algorithm 2:ILS for lower-level model.
Function ILS_LM(s, LS_LM, Perb_LM, MaxConsIter_LM) r≔s; r∗≔s; ConsIter _LM=0; while  ConsIter_LM<=MaxConsIter_LM  do r′≔RVND(r, LS_LW); if  RC(r′)<RC(r∗),then r∗≔ r′; ConsIter _LM≔0; else ConsIter _LM≔ ConsIter _LM +1; end r≔Perturbation(r∗, Perb_LW); end Return  r∗


## 4.5.2. ILS for Bilevel

Algorithm3 provides an outline of the ILS for the bilevel model. According to Algorithm 3, after constructing the initial solutions0 by applying the initialization procedure (see Section 4.3), the procedures RVND, ILS_LM (see Section 4.5.1), and Perturbation (see Section 4.4.4) are called repeatedly until the predefined termination criterion is satisfied. The best current solution is chosen to be perturbed by perturbation moves. In this procedure, lets represent the current solution for the bilevel model,s∗ represent the best solution that can be searched for at present,TC(s) represent the total penalties cost of current solutions, ConsIter_BL represent the current consecutive iterations without improvement, and MaxConsIter_BL represent the maximum consecutive iterations without improvement.Algorithm 3:ILS for bi-level model.
Function  ILS_BL(LS_LM, Perb_LM, LS_BL, Perb_BL, MaxConsIter_BL) s 0≔Initialization; s ∗≔s0; ConsIter_BL≔0; while  ConsIter_BL<=MaxConsIter_BL  do s′ ≔RVND(s0, LS_BL); s≔ILS_LM(s′,LS_LM, Perb_LM); ifTC(s )<TC(s∗),then s∗≔ s; ConsIter _BL≔0; else ConsIter _BL≔ ConsIter _BL +1; end s0≔Perturbation(s∗, Perb_BL); end Return  s∗


## 4.5.3. TS for Bilevel

The TS algorithm was introduced by Glover in 1986 and is an extension to the local search algorithm. It has been successfully applied to tackle a wide variety of routing problems, such as the traveling salesman problem, capacity vehicle routing problems, pickup and delivery problem, BRP [23, 26], and other routing problems. At each iteration, TS performs an exploration of the solution space by moving from a solutionw identified to the best solution in a subset of the neighborhoodN(w) ofw. To avoid being trapped in the local optima and prevent cycling on the same solution, TS uses an explicit memory structure (i.e., tabu list) to record the search track. If the move is in the tabu list, the move is forbidden. Sometimes tabu lists may prevent excellent solutions being chosen. The aspiration criterion provides the opportunity to overcome this obstacle. That is, the aspiration criterion enables superior neighboring solutions to be chosen, even though the solution is obtained by a banned move in the tabu list. We utilize two independent tabu lists, TABU1 and TABU2, as follows for the removal and insertion moves, respectively.TABU1: for a removal move, the tabu list records the station number which is removed from the current solution. When generating a neighborhood, if the station selected to be removed from the current solution exists in the tabu list, the move is marked as tabu. The length of the tabu list is denoted asl1.TABU2: for an insertion move, the tabu list records the station number which is inserted into the current solution. When generating a neighborhood, if the station selected to be inserted into the current solution exists in the tabu list, the move is marked as tabu. The length of the tabu list is denoted asl2.For a swap move, if stationi in the current solution is swapped with another stationj which is not in the current solution, it means that this swap move integrated removing stationiand inserting stationj. When generating a neighborhood by a swapping move, if stationi exists in the tabu list TABU1 or stationj exists in the tabu list TABU2, the swap move is marked as tabu.In fact, after generating an initial solution, the following two steps are repeated until the stopping criterion is met. The stopping criterion of TS is achieving a predetermined number of iterations with no improvement (MaxConsIter_TS).Step 1.
Generating a new solution: a temporary solutions′ is generated from all feasible neighborhood operators (LS_BL, see Section 4.4.2) of the current solutions. Then, the optimal low-level solution s″ is obtained throughILS_lower (see Section 4.5.1).Step 2.
Evaluating the new solution: in this step, the objective value of the new solution is compared with that of the best solution found. If the new solution improves the best solution, it is accepted and the best solution will be updated; otherwise, the new solution is accepted only if its move is not in the tabu list.In this TS_BL procedure, lets represent the current solution to the bilevel model,s∗ represent the best solution that can be searched for at present, tabu represent the two tabu lists,N(s) represent the collection of the neighborhood solutions generated by the current solutions, TC(s) represent the total penalties cost of current solutions,TC(s∗) represent the minimum total penalties cost so far, ConsIter_TS represent the current consecutive iterations without improvement, and MaxConsIter_TS represent the maximum consecutive iterations without improvement. The detailed pseudocode of TS for the bilevel model is presented  in Algorithm 4.Algorithm 4:TS for bi-level model.
Function  TS_BL(LS_LW, Perb_LW, LS_BL) s 0≔Initialization; Tuba=; s ∗≔s0; s≔s0; ConsIter_TS≔0; while  ConsIter_TS<=MaxConsIter_TS  do N(s) ≔ LocalSearch(s, LS_BL);(N(s) consists of all feasible neighbor solutions obtained by applying the three local search moves). Select a solutions′ in the set N(s) that minimizes TS(s), and the solution is either a non-tabu solution or a solution better than the best one so far (i.e.,TC(s′)<TC(s∗)), and update the tabu. s≔ILS_LM(s′, LS_LM, Perb_LM); ifTC(s)<TC(s∗),then s∗≔ s; ConsIter_TS≔0; else ConsIter_TS≔ ConsIter_TS+1; end end Return  s∗


## 4.6. Repairing Infeasible Solutions

Sometimes the neighbor solution may no longer remain feasible with reference to the loading and capacity constraints after implementing the Insert and Delete neighborhood moves. Hence, the solution can be repaired by adjusting the pick-up/drop-off quantities at one of the stations on the new route visited by the truck.Letω= (x,d,q) be the current solution, set x=(x0,x1,…,xn,xn+1),x0=xn+1=0, let S(x) be the set of stations in the route, S(x)={x1,…,xn}. Let σ=(σ1,…,σnσ) be a subsequence of a routex of solutionw, let djσ be the loading and unloading quantities at the station σj, let qjσ be the cumulative load at station σj, and then qjσ=∑i=1jdiσ.Letqminσ=min{0,q1σ,…,qnσσ},qmaxσ=max{0,q1σ,…,qnσσ}; that is to say, qminσis the minimum cumulative load and qmaxσ is the maximum cumulative load. According to Bulhões et al. [20], we can obtain the minimum and maximum flow of bikes allowed to enter the subsequence σ so as to ensure route feasibility. The minimum flow of bikes allowed to enter so as to ensure feasibility is lminσ=-qminσ; the maximum flow of bikes allowed to enter so as to ensure feasibility is lmaxσ=Q-qmaxσ.Delete: If we delete the station xh from the current solution, the loading and unloading quantities of the station xh in the current solution are dh, set ΔD=dh. After performing the delete move, let the new route of the new solution be xD=(x0,x1,…,xh-1,xh,…,xn+1)=(x0D,x1D,…xnDD,xnD+1D),nD=n-1.LetS(xD)=x1D,…,xnDD be the set of stations in the route. According to constraints (3) and (5), the total quantities of the route xD reduce to ΔD; in order to maintain the quantity balance, we need to increase ΔDat another station in S(xD). First, we put all feasible stations xkD∈S(xD),k∈[1,nD] that can adjust the quantity to obtain a feasible solution by increasing ΔD into the set FD,FD⊆S(xD). Then, we choose one station from FD according to the predetermined rules. If we choose station x∗D,x∗D∈FD, then set dx∗D=dx∗D+ΔD. The characteristics of the stations in SF must first be analyzed, and then the rules to choose the station to generate the feasible neighbor solution will be proposed.∀ k ∈ [ 1 , n D ], xD=(x0D,x1D,…xnDD,xnD+1D),nD=n, let σD=(xkD,…,xnDD) be a subsequence of a route xD. The station xkDcan be put into FD only if all the following conditions are satisfied:(21)ifxkD∈P,anddxkD+ΔD>0;dxkD+ΔD≤sxkD0;qksum+ΔD∈lminσD,lmaxσD.ifxkD∈D,anddxkD+ΔD<0;dxkD+ΔD≤cxkD-sxkD0;qksum+ΔD∈lminσD,lmaxσD.There are two rules to choose the station by adjusting the quantities to generate the feasible new solution:Best Choose Criterion andRandom Choose Criterion. The Best Choose Criterion is to choose the stationi from set FD with the least increasing penalty cost. Stationi is determined from argmini∈FD{f(siI-ΔD)-f(siI)}. The Random Choose Criterion is to randomly choose the stationi from set FD.Insert: If we insert a station m∈N∖S(x) into the current route between xj and xj+1, the loading and unloading quantities of the new stationm are dm, and let ΔI=dm. After performing the insert move, set the new route of the new solution as xI=(x0,x1,…xj,h,xj+1,…,xn,xn+1)=(x0I,x1I,…,xnII,xnI+1I),nI=n+1.SetS(xI)={x1I,…,xnII} as the set of stations in the new route xI. According to constraints (3) and (5), the total quantities of the route xI increase ΔI, and in order to maintain the quantity balance, we need to decrease ΔD at another station in S(xI). First, we put all feasible stations xkI∈S(xI),k∈[1,j]∪[j+1,nI] in S(xI), which can adjust the quantity to obtain a feasible solution through decreasing ΔI into the set FI,FI⊆S(xD). Then, we choose one station from FI according to theBest/Random Choose Criterion. The Best Choose Criterion is to choose the stationi from set FI with the least increasing penalty cost. The stationi is determined from argmini∈FI{f(siI+ΔI)-f(siI)}. The Random Choose Criterion is to randomly choose the stationi from set FI.If we choose stationx∗I,x∗I∈FI, then set dx∗I=dx∗I-ΔI. Obviously, inserting a station into the route will need more time for traveling. In addition to satisfying this truck capacity constraint, a feasible neighbor solution after performing the insert move also needs to obey the time constraint of the truck.∀ k ∈ [ 1 , j ] ∪ [ j + 2 , n I ], xI=(x0I,x1I,…,xnII,xnI+1I),nI=n+1, let σI=(xkI,…,xnII) be a subsequence of a route xI.The station xkIcan be put into FI only if all the following conditions are satisfied:(22)ifxkI∈P,anddxkI-ΔI>0;dxkI-ΔI≤sxkD0;qksum-ΔI∈lminσI,lmaxσITimexD+∑j∈Sx,j≠kd+dxkD-ΔD≤T.ifxkI∈D,anddxkI-ΔI<0;dxkI-ΔI≤cxkI-sxkI0;qksum-ΔI∈lminσI,lmaxσITimexD+∑j∈Sx,j≠kd+dxkD-ΔD≤T.


## 5. Computational Results

The heuristic algorithm presented in Section4 was coded in MATLAB, and all computational experiments were carried out on a Lenovo notebook with an Intel Core CPU@2.5 Hz. To tune our proposed heuristic and illustrate its efficiency and efficacy, sets available at http://www.eng.tau.ac.il/~talraviv/Publications/3step%20data.zip were used. These sets contain instances ranging from 20 to 200 stations, which were first used by Forma et al. [15]. Two vehicle capacities were used,Q=10 andQ=20, two repositioning times (in seconds) were used,T=14400 andT=21600, and two times for loading and unloading time (in seconds) were considered asL=U=30,L=U=60. Other input data (i.e., si0, ci, and fi(si)) were taken from Forma et al. [15]. The algorithm was run 20 times for each instance; the best and average of the solutions and the average computing time of 20 runs were used to evaluate the performance of the solution algorithm.The proposed algorithm ILS_BL relies mainly on two parameters, namely, the maximum consecutive iterations without improvement in the ILS_LM (MaxConsIter_LM) and the maximum consecutive iterations without improvement in the ILS_BL (MaxConsIter_BL). The TS_BL relies on four parameters, namely, the maximum consecutive iterations without improvement in the ILS_LM (MaxConsIter_LM), the maximum consecutive iterations without improvement in the TS_BL (MaxConsIter_TS), and the lengths of the two tabu lists (l1andl2). The value of the parameters MaxConsIter_LM and MaxConsIter_BL in ILS, and the parameter MaxConsIter_TS in TS were defined as a function of N (the number of all stations) andLen (the number of stations visited by the repositioning truck). Similar to Ho and Szeto [26] set l1=[(T/3600)log10(N),where[•] is the nearest integer function, and l2=[1.5l1]. After some preliminary experiments, the results from the algorithm were obtained by setting MaxConsIter_LM=[1.5∗Len], MaxConsIter_BL =[1.5∗Len], and MaxConsIter_TS=[1.5∗Len] when N<=100; MaxConsIter_LM=[2∗Len], MaxConsIter_BL =[2∗Len], and MaxConsIter_TS=[2∗Len] when N>100.


### 5.1. Effects of Different Initial Inventory Level on Solving the Bilevel BRP

As pointed out by Forma et al. [15], the initial inventory levels at the stations influence the amount of work necessary to reach their ideal inventory levels. The further the initial inventory level of a station is from its ideal level, the more necessary the work is. That is, the number of bicycles to load/unload to/from the station is higher. Therefore, similar to Forma et al. [15], we first test the performance of our solution algorithm based on three different cases of initial inventories called ‘‘light,” ‘‘real,’’ and ‘‘heavy,’’ referring to different workloads. For the ‘‘light’’ case, we generated initial inventories from a normal distribution with a mean that corresponds to the ideal inventory and a standard deviation that equals 0.2ci. In the ‘‘heavy’’ case, the standard deviation is equal to 0.6ci. The ‘‘real’’ case corresponds to real initial inventories that were observed in the system at midnight on a random day.Table5 provides some characteristics of the instances considered in this experimental study. In the first column in Table 5, N denotes the number of stations. The columns entitled “Light level,” “Real level,” and “Heavy level” denote the three different workloads of each instance. “Initial cost” represents the expected number of shortage events, assuming that no repositioning is done. The last column entitled “Ideal cost” presents the expected number of shortage events, assuming that the number of bikes at each station before the demand commences is optimal (i.e., ideal, ignoring the repositioning effort costs and constraints). The column entitled “Δ” represents the difference between the “Initial cost” and the “Ideal cost” of the corresponding workload. The value of Δ becomes larger when N increases. This means that the amount of repositioning work increases as the number of stations increases.Table 5
Characterization of the numerical instances.
N Initial cost  Ideal cost Light level Δ Real level Δ Heavy level Δ 20 159.9751 24.6652 159.4013 24.0914 251.0328 115.7229 135.3099  40 304.3676 40.6803 317.6306 53.9433 467.3612 203.6739 263.6873  60 458.1260 60.9716 462.8977 65.7433 707.0593 309.9049 397.1544  80 594.52439 68.0321 601.1442 74.6524 870.232 343.7424 526.4918  100 734.5029 86.7962 737.1299 89.4232 1048.1935 400.4868 647.7067  120 898.2554 118.6006 884.1184 104.4636 1286.6367 506.9819 779.6548  140 1057.5950 140.8475 1042.6900 125.9425 1510.8252 594.0777 916.7475  160 1203.5684 157.8831 1205.3941 159.7088 1699.9332 654.2479 1045.6853  180 1378.5382 195.3434 1388.4729 205.2781 1986.7625 803.5677 1183.1948  200 1566.6917 227.475 1548.7675 219.5508 2300.6803 971.4636 1329.2167Table6 reports the results obtained by settingQ=20,U=L=60,T=14400, which are under different initial inventory levels of stations. “Best” are the best objective values (zupper) obtained by the solution algorithm, “Avg” are the average objective values (zupper), and “Cpu” is the average computing time in seconds of the 20 runs. In order to measure the performance of the solution algorithm for all the different instances ranging from 20 to 200 stations, we consider the percentage of completion repositioning (PCR), which can be defined as follows: (23)PercentageofCompletionRepositioning=Initialcost-zupperInitialcost-Idealcost×100%.The value of the PCR represents customer satisfaction; the higher the value, the greater the degree of customer satisfaction. For the upper-level decision makers, there is no way to make the PCR 100%, due to the limited resources of the lower-level operator, such as the repositioning time and the capacity of the truck. Obviously, the value of PCR increases with the increase of the number of bikes picked up or dropped off at the stations visited. PCRB is the best percentage of repositioning completion and PCRA is the average percentage of repositioning completion, which are defined as follows:(24)PCRB=Initialcost-zupperBestInitialcost-Idealcost×100%,zupperBestisthebestobjectivevalueofupperlevelPCRA=Initialcost-zupperAvgInitialcost-Idealcost×100%,zupperAvgistheaverageobjectivevalueofupperlevel.When the workload is at the real level, the runtime is smaller than that for the other two levels. The runtime for the light level (3.8320) is smaller than that of the heavy level (4.6087). This is because demand that is too large and too small will increase the chance of having infeasible solutions. The results of PCRB and PCRA shown in Table 6 indicate the opposite trend as the runtime. As expected, in all the different workloads, the PCRB and PCRA decrease as the number of stations N increases. When the workload is at a real level, the PCRB and PCRA are higher than those of the other two levels. The PCRB and PCRA for the light level are higher than those of the heavy level. However, the runtime does not increase with increasing N, as shown in Table 6. The runtime decreases as N increases. This may be because an increase in N prompts the BSS providers to easily choose stations with a large number of shortages to minimize the penalty function, which leads to the truck needing more time to load/unload bikes. The larger value of loading/unloading time leaves little time for the truck to travel and the number of stations the truck visits is smaller, thus shrinking the solution space.Table 6
Results comparison of different workloads (Q=20,L=U=60, T=14400).
N Light level Real level Heavy level Best Avg Cpu P C R B (%) P C R A(%) Best Avg Cpu P C R B (%) P C R A(%) Best Avg Cpu P C R B (%) P C R A (%) 20 135.3294 135.3517 6.0980 99.9209 99.8305 138.1661 138.1661 5.9453 88.1443 88.1443 154.5146 154.5998 7.3875 83.4030 83.3310  40 270.2953 271.4314 5.7781 83.7563 80.9635 271.2084 271.824 4.0347 86.0574 84.9162 356.7850 357.7206 6.8385 54.2908 53.8314  60 421.6972 422.4598 4.4906 59.7472 58.4964 415.5485 418.1919 3.3344 72.0213 68.0005 593.7192 594.7404 5.5156 36.5725 36.2430  80 557.5409 558.8367 5.4187 54.3611 52.4564 553.2614 554.6676 4.8453 64.1410 62.2573 759.5554 761.2930 3.7969 32.1982 31.6927  100 699.3996 700.4135 3.2047 40.4434 39.2752 694.8413 696.0872 2.6360 47.2904 45.8971 936.6365 939.3733 5.8656 27.8554 27.1720  120 854.5804 856.4121 3.7313 36.8253 35.2809 839.5803 842.9962 3.0297 42.6350 39.3651 1181.2395 1184.1772 4.3281 20.7891 20.2097  140 1011.7422 1012.9844 3.3391 32.5549 31.6730 998.4752 999.0857 2.2234 35.1071 34.6224 1402.5349 1405.9053 3.0641 18.2283 17.6610  160 1157.7237 1159.6947 2.2297 29.0371 27.7887 1153.7849 1154.4173 1.4031 32.3146 31.9186 1585.4935 1592.0343 3.6047 17.4918 16.4921  180 1328.6067 1331.4690 2.0203 25.5609 24.0956 1335.9792 1337.0632 1.8750 25.5720 25.0439 1878.4544 1881.6833 2.5594 13.4784 13.0766  200 1504.3961 1505.5343 2.0094 22.9896 22.4892 1495.6188 1497.4245 1.8703 24.2079 23.3855 2190.1434 2192.0859 3.1266 11.3784 11.1784  Avg   3.8320 48.5197 47.2349   3.1197 51.7491 50.3551   4.6087 31.5686 31.0888


### 5.2. Effectiveness of the ILS_BL Compared with TS_BL

To demonstrate the effectiveness of ILS_BL, comparisons of the performance of the ILS_BL and TS_BL were conducted. Table7 reports the results obtained under the “real level” workload of each station withQ=20,L=U=60, andT=14400. GapB (%) and GapA (%) indicate the performance of the ILS_BL relative to that of TS_BL based on the best and the average PCR, respectively. The value of Gap implies the relative improvement (or deterioration) on the PCR of ILS_BL compared to that of the TS_BL. If it is positive, the performance of the ILS_BL is better than the TS_BL. From Table 7, all GapB and GapA values are positive, which means that the quality of the ILS_BL solution is better than that of the TS_BL. On average, the improvement was 3.1434%, with TS_BL requiring slightly more computing time. As expected, under a given repositioning time and truck capacity, when the size N becomes larger, the values of PCRB and PCRA become smaller. When N is 80, the values of GapB and GapA take the maximum value. This indicates that the ILS_BL is more effective than the TS_BL to solve the medium-size BRP.Table 7
Results comparison of ILS_BL and TS_BL (Q=20,L=U=60, T=14400).
N TS_BL ILS_BL Best Avg Cpu P C R B (%) P C R A (%) Best Avg Cpu P C R B (%) P C R A (%) G a p B G a p A 20 138.4328 138.4328 7.4031 87.0373 87.0373 138.1661 138.1661 5.9453 88.1443 88.1443 1.1070 1.1070  40 271.7736 272.4146 5.4016 85.0096 83.8213 271.2084 271.8240 4.0347 86.0574 84.9162 1.0478 1.0949  60 419.9921 422.2647 4.8875 65.2623 61.8055 415.5485 418.1919 3.3344 72.0213 68.0005 6.7590 6.1950  80 559.8575 560.9618 5.1563 55.3053 53.8260 553.2614 554.6676 4.8453 64.1410 62.2573 8.8357 8.4313  100 698.6127 699.7694 3.0641 43.0729 41.7794 694.8413 696.0872 2.6360 47.2904 45.8971 4.2175 4.1177  120 843.4800 846.9870 3.4813 38.9020 35.5448 839.5803 842.9962 3.0297 42.6350 39.3651 3.7330 3.8203  140 1002.6325 1003.3431 2.3906 31.8062 31.2420 998.4752 999.0857 2.2234 35.1071 34.6224 3.3009 3.3804  160 1155.8360 1156.7478 1.6354 31.0303 30.4594 1153.7849 1154.4173 1.4031 32.3146 31.9186 1.2843 1.4592  180 1337.7577 1339.0367 2.0235 24.7056 24.0825 1335.9792 1337.0632 1.8750 25.5720 25.0439 0.8664 0.9614  200 1498.2623 1499.3265 2.2625 23.0039 22.5192 1495.6188 1497.4245 1.8703 24.2079 23.3855 1.2040 0.8663  Avg   3.7706 48.5135 47.2117   3.1197 51.7491 50.3551 3.2356 3.1434


### 5.3. Effects of DifferentQ,U,L, andT for the Upper-Level Cost

In this section, we are interested in evaluating the impact of the truck capacity, the loading and unloading time per bike, and the repositioning time.Table8 shows the results comparison betweenQ=20 andQ=10 underU=L=60,T=14400 with sizes ranging from 20 to 200 stations. GapB and GapA indicate the performance ofQ=20 relative to that ofQ=10 based on the best and the average PCR, respectively. All negative GapB and GapA values indicate that as the capacity of the repositioning truck decreases, the objective function value increases, and on average the deterioration was 2.4245%. The Cpu of the ILS_BL required to solve the two instances are 3.1197 and 4.0147, respectively, which indicates that ILS_BL needs a little more time to obtain solutions as the truck capacity becomes smaller. This may be because a decrease in truck capacity increases the chance of having infeasible solutions.Table 8
Results comparison of different vehicle capacitiesQ(U=L=60, T=14400).
N Q=20 Q=10 Best Avg Cpu P C R B (%) P C R A (%) Best Avg Cpu P C R B (%) P C R A (%) G a p B G a p A 20 138.1661 138.1661 5.9453 88.1443 88.1443 138.1661 138.1661 6.9360 88.1443 88.1443 0 0  40 271.2084 271.824 4.0347 86.0574 84.9162 273.2792 274.2790 5.6156 82.2186 80.3651 -3.8388 -4.5511  60 415.5485 418.1919 3.3344 72.0213 68.0005 419.4120 419.7309 3.9047 66.1447 65.6596 -5.8766 -2.3409  80 553.2614 554.6676 4.8453 64.1410 62.2573 556.7589 558.0334 4.9063 59.4560 57.7487 -4.6850 -4.5086  100 694.8413 696.0872 2.6360 47.2904 45.8971 694.9719 697.3861 3.4297 47.1444 44.4446 -0.1460 -1.4525  120 839.5803 842.9962 3.0297 42.6350 39.3651 841.0893 843.126 4.3875 41.1905 39.2408 -1.4445 -0.1243  140 998.4752 999.0857 2.2234 35.1071 34.6224 1001.7933 1004.6089 3.2047 32.4725 30.2369 -2.6346 -4.3855  160 1153.7849 1154.4173 1.4031 32.3146 31.9186 1160.7397 1162.2074 1.7469 27.9599 27.0409 -4.3547 -4.8777  180 1335.9792 1337.0632 1.8750 25.5720 25.0439 1336.9554 1340.2226 2.7328 25.0964 23.5048 -0.4756 -1.5391  200 1495.6188 1497.4245 1.8703 24.2079 23.3855 1497.3002 1498.4470 3.2828 23.4421 22.9198 -0.7658 -0.4657  Avg   3.1197 51.7491 50.3551   4.0147 49.3269 47.9306 -2.4222 -2.4245Table9 shows the results comparison between loading/unloading timeL=U=60 and loading/unloading timeL=U=30 underQ=20,T=14400, with sizes ranging from 20 to 200 stations. GapB and GapA indicate the performance ofL=U=60 relative to that ofL=U=30 based on the best and the average PCR, respectively. All positive GapB and GapA values indicate that as the loading/unloading time increases, the objective function value increases, and on average the improvement was 8.1178%. The Cpu of the ILS_BL required to solve the two instances is 3.1197 and 4.2709, respectively, which indicates that a decrease in loading/unloading time leads to a longer computation time. This may be explained by the fact that if the loading/unloading time of bikes is shortened, there is more time for the truck to visit more stations, resulting in an expansion of the solution space.Table 9
Results comparison of different unloading and loading service times (Q=20, T=14400).
N L = U =60 L = U =30 Best Avg Cpu P C R B (%) P C R A (%) Best Avg Cpu P C R B (%) P C R A (%) G a p B G a p A 20 138.1661 138.1661 5.9453 88.1443 88.1443 138.1661 138.1661 6.8984 88.1443 88.1443 0 0  40 271.2084 271.824 4.0347 86.0574 84.9162 264.9843 265.3725 5.8828 97.5956 96.8760 11.5382 11.9598  60 415.5485 418.1919 3.3344 72.0213 68.0005 409.9184 411.0709 4.7121 80.5851 78.8321 8.5638 10.8316  80 553.2614 554.6676 4.8453 64.1410 62.2573 545.789 547.3909 5.9654 74.1506 72.0048 10.0096 9.7475  100 694.8413 696.0872 2.6360 47.2904 45.8971 683.9921 684.5961 3.7500 59.4228 58.7474 12.1324 12.8503  120 839.5803 842.9962 3.0297 42.6350 39.3651 832.4907 833.3217 3.7250 49.4217 48.6262 6.7867 9.2611  140 998.4752 999.0857 2.2234 35.1071 34.6224 987.7048 988.9890 4.0281 43.6590 42.6393 8.5519 8.0169  160 1153.7849 1154.4173 1.4031 32.3146 31.9186 1141.8162 1142.4030 3.6500 39.8086 39.4411 7.4940 7.5225  180 1335.9792 1337.0632 1.8750 25.5720 25.0439 1322.5730 1324.8880 2.0813 32.1027 30.9749 6.5307 5.9310  200 1495.6188 1497.4245 1.8703 24.2079 23.3855 1484.8130 1486.3200 2.0156 29.1297 28.4432 4.9218 5.0577  Avg   3.1197 51.7491 50.3551   4.2709 59.4020 58.4729 7.6529 8.1178Table10 shows the results comparison between repositioning timeT=14400 andT=21600 underU=L=60,Q=20 with sizes ranging from 20 to 200 stations. GapB and GapA indicate the performance ofT=14400 relative to that ofT=21600 based on the best and the average PCR, respectively. All positive GapB and GapA values indicate that as the repositioning time increases, the objective function value decreases, and on average the deterioration was 14.4469%. The Cpu of the ILS_BL required to solve the two instances is 3.1197 and 4.4563, respectively, which indicates that an increase in repositioning time leads to a longer computation time. This is expected because the truck can visit more stations with the increased repositioning time, resulting in an expansion of the solution space.Table 10
Results comparison of different repositioning times (Q=20, U=L=60).
N T=14400 T=21600 Best Avg Cpu P C R B (%) P C R A (%) Best Avg Cpu P C R B (%) P C R A (%) G a p B G a p A 20 138.1661 138.1661 5.9453 88.1443 88.1443 138.1661 138.1661 7.1078 88.1443 88.1443 0 0  40 271.2084 271.824 4.0347 86.0574 84.9162 263.8787 264.0471 6.1563 99.6452 99.3330 13.5878 14.4168  60 415.5485 418.1919 3.3344 72.0213 68.0005 404.8966 405.6153 4.8651 88.2236 87.1304 16.2023 19.1299  80 553.2614 554.6676 4.8453 64.1410 62.2573 541.4491 542.3618 6.3177 79.9641 78.7415 15.8231 16.4842  100 694.8413 696.0872 2.6360 47.2904 45.8971 674.1251 675.9928 4.1979 70.4569 68.3683 23.1665 22.4712  120 839.5803 842.9962 3.0297 42.6350 39.3651 822.3390 823.2306 4.5764 59.1396 58.2861 16.5046 18.9210  140 998.4752 999.0857 2.2234 35.1071 34.6224 976.7542 977.8183 3.6290 52.3539 51.5090 17.2468 16.8866  160 1153.7849 1154.4173 1.4031 32.3146 31.9186 1131.2477 1132.9226 2.2207 46.4260 45.3773 14.1114 13.4587  180 1335.9792 1337.0632 1.8750 25.5720 25.0439 1309.3586 1310.8218 2.8672 38.5401 37.8273 12.9681 12.7834  200 1495.6188 1497.4245 1.8703 24.2079 23.3855 1474.0601 1475.6522 2.6250 34.0274 33.3022 9.8195 9.9167  Avg   3.1197 51.7491 50.3551   4.4563 65.6921 64.8019 13.9430 14.4469


## 5.1. Effects of Different Initial Inventory Level on Solving the Bilevel BRP

As pointed out by Forma et al. [15], the initial inventory levels at the stations influence the amount of work necessary to reach their ideal inventory levels. The further the initial inventory level of a station is from its ideal level, the more necessary the work is. That is, the number of bicycles to load/unload to/from the station is higher. Therefore, similar to Forma et al. [15], we first test the performance of our solution algorithm based on three different cases of initial inventories called ‘‘light,” ‘‘real,’’ and ‘‘heavy,’’ referring to different workloads. For the ‘‘light’’ case, we generated initial inventories from a normal distribution with a mean that corresponds to the ideal inventory and a standard deviation that equals 0.2ci. In the ‘‘heavy’’ case, the standard deviation is equal to 0.6ci. The ‘‘real’’ case corresponds to real initial inventories that were observed in the system at midnight on a random day.Table5 provides some characteristics of the instances considered in this experimental study. In the first column in Table 5, N denotes the number of stations. The columns entitled “Light level,” “Real level,” and “Heavy level” denote the three different workloads of each instance. “Initial cost” represents the expected number of shortage events, assuming that no repositioning is done. The last column entitled “Ideal cost” presents the expected number of shortage events, assuming that the number of bikes at each station before the demand commences is optimal (i.e., ideal, ignoring the repositioning effort costs and constraints). The column entitled “Δ” represents the difference between the “Initial cost” and the “Ideal cost” of the corresponding workload. The value of Δ becomes larger when N increases. This means that the amount of repositioning work increases as the number of stations increases.Table 5
Characterization of the numerical instances.
N Initial cost  Ideal cost Light level Δ Real level Δ Heavy level Δ 20 159.9751 24.6652 159.4013 24.0914 251.0328 115.7229 135.3099  40 304.3676 40.6803 317.6306 53.9433 467.3612 203.6739 263.6873  60 458.1260 60.9716 462.8977 65.7433 707.0593 309.9049 397.1544  80 594.52439 68.0321 601.1442 74.6524 870.232 343.7424 526.4918  100 734.5029 86.7962 737.1299 89.4232 1048.1935 400.4868 647.7067  120 898.2554 118.6006 884.1184 104.4636 1286.6367 506.9819 779.6548  140 1057.5950 140.8475 1042.6900 125.9425 1510.8252 594.0777 916.7475  160 1203.5684 157.8831 1205.3941 159.7088 1699.9332 654.2479 1045.6853  180 1378.5382 195.3434 1388.4729 205.2781 1986.7625 803.5677 1183.1948  200 1566.6917 227.475 1548.7675 219.5508 2300.6803 971.4636 1329.2167Table6 reports the results obtained by settingQ=20,U=L=60,T=14400, which are under different initial inventory levels of stations. “Best” are the best objective values (zupper) obtained by the solution algorithm, “Avg” are the average objective values (zupper), and “Cpu” is the average computing time in seconds of the 20 runs. In order to measure the performance of the solution algorithm for all the different instances ranging from 20 to 200 stations, we consider the percentage of completion repositioning (PCR), which can be defined as follows: (23)PercentageofCompletionRepositioning=Initialcost-zupperInitialcost-Idealcost×100%.The value of the PCR represents customer satisfaction; the higher the value, the greater the degree of customer satisfaction. For the upper-level decision makers, there is no way to make the PCR 100%, due to the limited resources of the lower-level operator, such as the repositioning time and the capacity of the truck. Obviously, the value of PCR increases with the increase of the number of bikes picked up or dropped off at the stations visited. PCRB is the best percentage of repositioning completion and PCRA is the average percentage of repositioning completion, which are defined as follows:(24)PCRB=Initialcost-zupperBestInitialcost-Idealcost×100%,zupperBestisthebestobjectivevalueofupperlevelPCRA=Initialcost-zupperAvgInitialcost-Idealcost×100%,zupperAvgistheaverageobjectivevalueofupperlevel.When the workload is at the real level, the runtime is smaller than that for the other two levels. The runtime for the light level (3.8320) is smaller than that of the heavy level (4.6087). This is because demand that is too large and too small will increase the chance of having infeasible solutions. The results of PCRB and PCRA shown in Table 6 indicate the opposite trend as the runtime. As expected, in all the different workloads, the PCRB and PCRA decrease as the number of stations N increases. When the workload is at a real level, the PCRB and PCRA are higher than those of the other two levels. The PCRB and PCRA for the light level are higher than those of the heavy level. However, the runtime does not increase with increasing N, as shown in Table 6. The runtime decreases as N increases. This may be because an increase in N prompts the BSS providers to easily choose stations with a large number of shortages to minimize the penalty function, which leads to the truck needing more time to load/unload bikes. The larger value of loading/unloading time leaves little time for the truck to travel and the number of stations the truck visits is smaller, thus shrinking the solution space.Table 6
Results comparison of different workloads (Q=20,L=U=60, T=14400).
N Light level Real level Heavy level Best Avg Cpu P C R B (%) P C R A(%) Best Avg Cpu P C R B (%) P C R A(%) Best Avg Cpu P C R B (%) P C R A (%) 20 135.3294 135.3517 6.0980 99.9209 99.8305 138.1661 138.1661 5.9453 88.1443 88.1443 154.5146 154.5998 7.3875 83.4030 83.3310  40 270.2953 271.4314 5.7781 83.7563 80.9635 271.2084 271.824 4.0347 86.0574 84.9162 356.7850 357.7206 6.8385 54.2908 53.8314  60 421.6972 422.4598 4.4906 59.7472 58.4964 415.5485 418.1919 3.3344 72.0213 68.0005 593.7192 594.7404 5.5156 36.5725 36.2430  80 557.5409 558.8367 5.4187 54.3611 52.4564 553.2614 554.6676 4.8453 64.1410 62.2573 759.5554 761.2930 3.7969 32.1982 31.6927  100 699.3996 700.4135 3.2047 40.4434 39.2752 694.8413 696.0872 2.6360 47.2904 45.8971 936.6365 939.3733 5.8656 27.8554 27.1720  120 854.5804 856.4121 3.7313 36.8253 35.2809 839.5803 842.9962 3.0297 42.6350 39.3651 1181.2395 1184.1772 4.3281 20.7891 20.2097  140 1011.7422 1012.9844 3.3391 32.5549 31.6730 998.4752 999.0857 2.2234 35.1071 34.6224 1402.5349 1405.9053 3.0641 18.2283 17.6610  160 1157.7237 1159.6947 2.2297 29.0371 27.7887 1153.7849 1154.4173 1.4031 32.3146 31.9186 1585.4935 1592.0343 3.6047 17.4918 16.4921  180 1328.6067 1331.4690 2.0203 25.5609 24.0956 1335.9792 1337.0632 1.8750 25.5720 25.0439 1878.4544 1881.6833 2.5594 13.4784 13.0766  200 1504.3961 1505.5343 2.0094 22.9896 22.4892 1495.6188 1497.4245 1.8703 24.2079 23.3855 2190.1434 2192.0859 3.1266 11.3784 11.1784  Avg   3.8320 48.5197 47.2349   3.1197 51.7491 50.3551   4.6087 31.5686 31.0888


## 5.2. Effectiveness of the ILS_BL Compared with TS_BL

To demonstrate the effectiveness of ILS_BL, comparisons of the performance of the ILS_BL and TS_BL were conducted. Table7 reports the results obtained under the “real level” workload of each station withQ=20,L=U=60, andT=14400. GapB (%) and GapA (%) indicate the performance of the ILS_BL relative to that of TS_BL based on the best and the average PCR, respectively. The value of Gap implies the relative improvement (or deterioration) on the PCR of ILS_BL compared to that of the TS_BL. If it is positive, the performance of the ILS_BL is better than the TS_BL. From Table 7, all GapB and GapA values are positive, which means that the quality of the ILS_BL solution is better than that of the TS_BL. On average, the improvement was 3.1434%, with TS_BL requiring slightly more computing time. As expected, under a given repositioning time and truck capacity, when the size N becomes larger, the values of PCRB and PCRA become smaller. When N is 80, the values of GapB and GapA take the maximum value. This indicates that the ILS_BL is more effective than the TS_BL to solve the medium-size BRP.Table 7
Results comparison of ILS_BL and TS_BL (Q=20,L=U=60, T=14400).
N TS_BL ILS_BL Best Avg Cpu P C R B (%) P C R A (%) Best Avg Cpu P C R B (%) P C R A (%) G a p B G a p A 20 138.4328 138.4328 7.4031 87.0373 87.0373 138.1661 138.1661 5.9453 88.1443 88.1443 1.1070 1.1070  40 271.7736 272.4146 5.4016 85.0096 83.8213 271.2084 271.8240 4.0347 86.0574 84.9162 1.0478 1.0949  60 419.9921 422.2647 4.8875 65.2623 61.8055 415.5485 418.1919 3.3344 72.0213 68.0005 6.7590 6.1950  80 559.8575 560.9618 5.1563 55.3053 53.8260 553.2614 554.6676 4.8453 64.1410 62.2573 8.8357 8.4313  100 698.6127 699.7694 3.0641 43.0729 41.7794 694.8413 696.0872 2.6360 47.2904 45.8971 4.2175 4.1177  120 843.4800 846.9870 3.4813 38.9020 35.5448 839.5803 842.9962 3.0297 42.6350 39.3651 3.7330 3.8203  140 1002.6325 1003.3431 2.3906 31.8062 31.2420 998.4752 999.0857 2.2234 35.1071 34.6224 3.3009 3.3804  160 1155.8360 1156.7478 1.6354 31.0303 30.4594 1153.7849 1154.4173 1.4031 32.3146 31.9186 1.2843 1.4592  180 1337.7577 1339.0367 2.0235 24.7056 24.0825 1335.9792 1337.0632 1.8750 25.5720 25.0439 0.8664 0.9614  200 1498.2623 1499.3265 2.2625 23.0039 22.5192 1495.6188 1497.4245 1.8703 24.2079 23.3855 1.2040 0.8663  Avg   3.7706 48.5135 47.2117   3.1197 51.7491 50.3551 3.2356 3.1434


## 5.3. Effects of DifferentQ,U,L, andT for the Upper-Level Cost

In this section, we are interested in evaluating the impact of the truck capacity, the loading and unloading time per bike, and the repositioning time.Table8 shows the results comparison betweenQ=20 andQ=10 underU=L=60,T=14400 with sizes ranging from 20 to 200 stations. GapB and GapA indicate the performance ofQ=20 relative to that ofQ=10 based on the best and the average PCR, respectively. All negative GapB and GapA values indicate that as the capacity of the repositioning truck decreases, the objective function value increases, and on average the deterioration was 2.4245%. The Cpu of the ILS_BL required to solve the two instances are 3.1197 and 4.0147, respectively, which indicates that ILS_BL needs a little more time to obtain solutions as the truck capacity becomes smaller. This may be because a decrease in truck capacity increases the chance of having infeasible solutions.Table 8
Results comparison of different vehicle capacitiesQ(U=L=60, T=14400).
N Q=20 Q=10 Best Avg Cpu P C R B (%) P C R A (%) Best Avg Cpu P C R B (%) P C R A (%) G a p B G a p A 20 138.1661 138.1661 5.9453 88.1443 88.1443 138.1661 138.1661 6.9360 88.1443 88.1443 0 0  40 271.2084 271.824 4.0347 86.0574 84.9162 273.2792 274.2790 5.6156 82.2186 80.3651 -3.8388 -4.5511  60 415.5485 418.1919 3.3344 72.0213 68.0005 419.4120 419.7309 3.9047 66.1447 65.6596 -5.8766 -2.3409  80 553.2614 554.6676 4.8453 64.1410 62.2573 556.7589 558.0334 4.9063 59.4560 57.7487 -4.6850 -4.5086  100 694.8413 696.0872 2.6360 47.2904 45.8971 694.9719 697.3861 3.4297 47.1444 44.4446 -0.1460 -1.4525  120 839.5803 842.9962 3.0297 42.6350 39.3651 841.0893 843.126 4.3875 41.1905 39.2408 -1.4445 -0.1243  140 998.4752 999.0857 2.2234 35.1071 34.6224 1001.7933 1004.6089 3.2047 32.4725 30.2369 -2.6346 -4.3855  160 1153.7849 1154.4173 1.4031 32.3146 31.9186 1160.7397 1162.2074 1.7469 27.9599 27.0409 -4.3547 -4.8777  180 1335.9792 1337.0632 1.8750 25.5720 25.0439 1336.9554 1340.2226 2.7328 25.0964 23.5048 -0.4756 -1.5391  200 1495.6188 1497.4245 1.8703 24.2079 23.3855 1497.3002 1498.4470 3.2828 23.4421 22.9198 -0.7658 -0.4657  Avg   3.1197 51.7491 50.3551   4.0147 49.3269 47.9306 -2.4222 -2.4245Table9 shows the results comparison between loading/unloading timeL=U=60 and loading/unloading timeL=U=30 underQ=20,T=14400, with sizes ranging from 20 to 200 stations. GapB and GapA indicate the performance ofL=U=60 relative to that ofL=U=30 based on the best and the average PCR, respectively. All positive GapB and GapA values indicate that as the loading/unloading time increases, the objective function value increases, and on average the improvement was 8.1178%. The Cpu of the ILS_BL required to solve the two instances is 3.1197 and 4.2709, respectively, which indicates that a decrease in loading/unloading time leads to a longer computation time. This may be explained by the fact that if the loading/unloading time of bikes is shortened, there is more time for the truck to visit more stations, resulting in an expansion of the solution space.Table 9
Results comparison of different unloading and loading service times (Q=20, T=14400).
N L = U =60 L = U =30 Best Avg Cpu P C R B (%) P C R A (%) Best Avg Cpu P C R B (%) P C R A (%) G a p B G a p A 20 138.1661 138.1661 5.9453 88.1443 88.1443 138.1661 138.1661 6.8984 88.1443 88.1443 0 0  40 271.2084 271.824 4.0347 86.0574 84.9162 264.9843 265.3725 5.8828 97.5956 96.8760 11.5382 11.9598  60 415.5485 418.1919 3.3344 72.0213 68.0005 409.9184 411.0709 4.7121 80.5851 78.8321 8.5638 10.8316  80 553.2614 554.6676 4.8453 64.1410 62.2573 545.789 547.3909 5.9654 74.1506 72.0048 10.0096 9.7475  100 694.8413 696.0872 2.6360 47.2904 45.8971 683.9921 684.5961 3.7500 59.4228 58.7474 12.1324 12.8503  120 839.5803 842.9962 3.0297 42.6350 39.3651 832.4907 833.3217 3.7250 49.4217 48.6262 6.7867 9.2611  140 998.4752 999.0857 2.2234 35.1071 34.6224 987.7048 988.9890 4.0281 43.6590 42.6393 8.5519 8.0169  160 1153.7849 1154.4173 1.4031 32.3146 31.9186 1141.8162 1142.4030 3.6500 39.8086 39.4411 7.4940 7.5225  180 1335.9792 1337.0632 1.8750 25.5720 25.0439 1322.5730 1324.8880 2.0813 32.1027 30.9749 6.5307 5.9310  200 1495.6188 1497.4245 1.8703 24.2079 23.3855 1484.8130 1486.3200 2.0156 29.1297 28.4432 4.9218 5.0577  Avg   3.1197 51.7491 50.3551   4.2709 59.4020 58.4729 7.6529 8.1178Table10 shows the results comparison between repositioning timeT=14400 andT=21600 underU=L=60,Q=20 with sizes ranging from 20 to 200 stations. GapB and GapA indicate the performance ofT=14400 relative to that ofT=21600 based on the best and the average PCR, respectively. All positive GapB and GapA values indicate that as the repositioning time increases, the objective function value decreases, and on average the deterioration was 14.4469%. The Cpu of the ILS_BL required to solve the two instances is 3.1197 and 4.4563, respectively, which indicates that an increase in repositioning time leads to a longer computation time. This is expected because the truck can visit more stations with the increased repositioning time, resulting in an expansion of the solution space.Table 10
Results comparison of different repositioning times (Q=20, U=L=60).
N T=14400 T=21600 Best Avg Cpu P C R B (%) P C R A (%) Best Avg Cpu P C R B (%) P C R A (%) G a p B G a p A 20 138.1661 138.1661 5.9453 88.1443 88.1443 138.1661 138.1661 7.1078 88.1443 88.1443 0 0  40 271.2084 271.824 4.0347 86.0574 84.9162 263.8787 264.0471 6.1563 99.6452 99.3330 13.5878 14.4168  60 415.5485 418.1919 3.3344 72.0213 68.0005 404.8966 405.6153 4.8651 88.2236 87.1304 16.2023 19.1299  80 553.2614 554.6676 4.8453 64.1410 62.2573 541.4491 542.3618 6.3177 79.9641 78.7415 15.8231 16.4842  100 694.8413 696.0872 2.6360 47.2904 45.8971 674.1251 675.9928 4.1979 70.4569 68.3683 23.1665 22.4712  120 839.5803 842.9962 3.0297 42.6350 39.3651 822.3390 823.2306 4.5764 59.1396 58.2861 16.5046 18.9210  140 998.4752 999.0857 2.2234 35.1071 34.6224 976.7542 977.8183 3.6290 52.3539 51.5090 17.2468 16.8866  160 1153.7849 1154.4173 1.4031 32.3146 31.9186 1131.2477 1132.9226 2.2207 46.4260 45.3773 14.1114 13.4587  180 1335.9792 1337.0632 1.8750 25.5720 25.0439 1309.3586 1310.8218 2.8672 38.5401 37.8273 12.9681 12.7834  200 1495.6188 1497.4245 1.8703 24.2079 23.3855 1474.0601 1475.6522 2.6250 34.0274 33.3022 9.8195 9.9167  Avg   3.1197 51.7491 50.3551   4.4563 65.6921 64.8019 13.9430 14.4469


## 6. Conclusions

The bike repositioning problem plays an important role in bike sharing systems. In this paper, the bilevel model we have presented offers a new research direction and field for studying BRP under the situation that BSS providers outsource the transportation between bike stations to the third-party logistics companies (the truck owners). In this bilevel model, the upper-level decision maker is the BSS provider whose goal is to minimize the total penalties cost caused at all stations, while the truck owner as the lower-level decision maker aims to minimize the travel cost. We proposed the ILS and TS to solve the bilevel model.Computational experiments are set up to illustrate the performance of the proposed algorithm and the properties of the problem. In particular, the results show that ILS outperforms TS in its solution quality and convergence speed. Moreover, the results also demonstrate that the size of the repositioning truck, the value of loading/unloading time for a single bike, and the length of repositioning time can significantly affect the total penalties cost of BSS provider.This study can be extended in several directions in the future: first, multiple repositioning trucks used by the third-party logistics companies can be considered to make this study more practical. Second, this problem will be more practical and flexible if some assumptions are relaxed. For example, it will also be more practical to relax the assumption that a station will be visited by a vehicle more than once. Third, it would be interesting to investigate the impact of temporary loading/unloading bikes at balanced stations.



---

*Source: 8641492-2019-06-02.xml*