Journal of Advanced Transportation
Volume 2019 (2019), Article ID 8641492, 19 pages
https://doi.org/10.1155/2019/8641492

Research Article

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

Qiong Tang,1,2,3 Zhuo Fu,1,3 and Meng Qiu1,3

¹School of Traffic and Transportation Engineering, Central South University, Changsha 410075, China
²College of Economics and Management, Hengyang Normal University, Hengyang 421002, China
³Smart Transport Key Laboratory of Hunan Province, Hunan, China

Correspondence should be addressed to Qiong Tang; qiongtang@hynu.edu.cn

Received 12 February 2019; Accepted 16 May 2019; Published 2 June 2019

Academic Editor: Juan C. Cano

Copyright © 2019 Qiong Tang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 Section 3, 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 Table 1. 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
Yes

In terms of the problem objective functions addressed in the BRP, common measures are transportation cost or other related measures such as transportation time or distance, and 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 TS

To 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 Figure 1). 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:

The outsourcing cost of BSS provider outsourcing relocation transportation to the third-party logistics companies is fixed.

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.

The depot is assumed to have no bikes.

The load on the truck cannot be greater than the truck capacity.

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 Table 3 unless otherwise specified.

Table 3: Notations for the bi-level model.

N
Set of stations

No
Set of nodes, including the stations and the depot (indexed by 0)

i, j
Indices of stations

Initial inventory level at station i before the repositioning operation

Final inventory level at station i after the repositioning operation

Convex penalty function of for station i

The target optimal inventory at station i, ,i.e., the optimal number of bikes at station i

P
Set of pick-up stations, i.e.,

D
Set of drop-off stations, i.e.,

Q
Truck capacity

Capacity of station i

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

Traveling time from station i to station j

Binary variable that equals 1 if the truck travels directly from station i to station j, and 0 otherwise

Number of bikes on the truck when it travels directly from station i to station j

The loading quantities at station i,  

The unloading quantitiesat station i,  

The vector of the upper-level decision variables,  

The vector of the lower-level decision variables,  

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 station i, the penalty function fi(si) is a function of the final inventory of station i after the operation, which represents the anticipated shortage of station i. The penalty cost of station i can be calculated by the corresponding penalty function fi(si), and the value range of si is . If si is out of the domain , then fi(si)= .

As described by Forma et al. [15], all stations of the BSS were divided into three categories, namely, 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; 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;    ‘‘touristic’’ stations where the demand of customers to rent or return bikes does not change much throughout the day. The optimal inventory level at station i can be solved by .

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 , , and . 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 optimal at each station to minimize the total penalty costs and can be formulated as follows:

The upper-level objective function (1) minimizes the total penalty costs generated by all stations. Constraint (2) defines the final inventory at station i after the operation as the initial inventory at station i minus the loading quantities or plus the unloading quantities at station i. Constraint (3) requires that the loading or unloading quantities at station i are equal to the difference between the loads of the truck before and after the station i. 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 can be solved from the lower-level model, while constraints (3) and (5) describe the relationship between the upper-level decision variables and lower-level decision variables .

3.3.2. The Lower-Level Model

The lower-level model aims to determine the optimal to minimize the transportation cost and can be formulated as follows:

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 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 and lower-level decision variables . Note that, for a given upper-level solution , solving the lower-level problems yields . 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 Table 2, 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 Sections 4.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 solution s, to yield a local optimum s∗. 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 a Perturbation, leading to a solution . After applying LocalSearch, 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 the LocalSearch; (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.

 s₀≔ GenerateInitialSolution

s ≔ LocalSearch(s₀)

while  termination criterion is not met  do

        ≔ Perturbation(s,history);

         ≔ LocalSearch();

     s≔ AcceptanceCriterion(s, s, history);

    end while

Return  s

Algorithm 1: ILS.

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), where x 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, and q is a vector to store the truck load during the route.

Let n be the total number of stations (excluding the starting and ending depot) visited by the truck. The stations visited by the truck are denoted as xh, , and the repositioning truck is assumed to start and end repositioning at the depot (denoted by x₀, xn+1,and x₀=xn+1=0).

The elements of d are defined as . A positive integer represents a loading operation, while a negative integer indicates an unloading operation. According to the constraints (12), d₀=dn+1=0.

The last elements of q are defined as , where denotes the loading/unloading quantity at station xh. In other words, qh gives the value of the load on the truck after visiting station xh 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:

A feasible solution for the repositioning operation performed by the truck with the capacity of 10 bikes is shown in Figure 2. 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 are rather than used by Ho and Szeto [26], where . The main reason for improvement is that if the loading or unloading quantities at station i exceed Q, 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 Section 4.6) to the new solution feasible with the Best Choose Criterion (see Section 4.6).

Swap: The move replaces station i in the solution with another station j which is not in the current solution. The station i and station j 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 station j 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 Section 4.6) to make the neighbor solution feasible with the Best 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 Section 4.6) to the new solution feasible with the Random 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 Section 4.6) to make the neighbor solution feasible with the Random 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 Section 4.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

Perturbation 
mechanism
Perb_LM
Rand Exchange
Perb_BL
Rand Insert

Rand Relocation
Rand Delete

4.5.1. ILS for Lower Level

Algorithm 2 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 solution s 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, let r represent the current route, r represent the best route that can be searched for at present, RC(r) represent the travel cost of current route r, ConsIter_LM represent the current consecutive iterations without improvement, and MaxConsIter_LM represent the maximum consecutive iterations without improvement.

Function ILS_LM(s, LS_LM, Perb_LM, MaxConsIter_LM)

r≔s; r≔s; ConsIter _LM=0;

while  ConsIter_LM<=MaxConsIter_LM  do

 ≔RVND(r, LS_LW);

 if  RC()<RC(r),then

    r≔ ;

    ConsIter _LM≔0;

 else

 ConsIter _LM≔ ConsIter _LM +1;

 end

 r≔Perturbation(r, Perb_LW);

end

Return  r

Algorithm 2: ILS for lower-level model.

4.5.2. ILS for Bilevel

Algorithm 3 provides an outline of the ILS for the bilevel model. According to Algorithm 3, after constructing the initial solution s₀ 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, let s 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 solution s, ConsIter_BL represent the current consecutive iterations without improvement, and MaxConsIter_BL represent the maximum consecutive iterations without improvement.

Function  ILS_BL(LS_LM, Perb_LM, LS_BL, Perb_BL, MaxConsIter_BL)

s₀≔Initialization;

s≔s₀; ConsIter_BL≔0;

while  ConsIter_BL<=MaxConsIter_BL  do

          ≔RVND(s₀, LS_BL);

         s≔ILS_LM(,LS_LM, Perb_LM);

 if TC(s )<TC(s),then

     s≔ s;

    ConsIter _BL≔0;

 else

    ConsIter _BL≔ ConsIter _BL +1;

 end

        s₀≔Perturbation(s, Perb_BL);

end

Return  s

Algorithm 3: ILS for bi-level model.

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 solution w identified to the best solution in a subset of the neighborhood N(w) of w. 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 as l₁.

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 as l₂.

For a swap move, if station i in the current solution is swapped with another station j which is not in the current solution, it means that this swap move integrated removing station i and inserting station j. When generating a neighborhood by a swapping move, if station i exists in the tabu list TABU1 or station j 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 solution is generated from all feasible neighborhood operators (LS_BL, see Section 4.4.2) of the current solution s. Then, the optimal low-level solution is obtained through ILS_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, let s 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 solution s, TC(s) represent the total penalties cost of current solution s, 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.

Function  TS_BL(LS_LW, Perb_LW, LS_BL)

s₀≔Initialization; Tuba=;

s≔s₀; s≔s₀; 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 solution 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()<TC(s)),

and update the tabu.

    s≔ILS_LM(, LS_LM, Perb_LM);

        if TC(s)<TC(s),then

         s≔ s;

         ConsIter_TS≔0;

        else

         ConsIter_TS≔ ConsIter_TS+1;

        end

  end

    Return  s

Algorithm 4: TS for bi-level model.

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 , let be the set of stations in the route, . Let be a subsequence of a route x of solution w, let be the loading and unloading quantities at the station , let be the cumulative load at station , and then .

Let ,; that is to say, is the minimum cumulative load and 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 ; the maximum flow of bikes allowed to enter so as to ensure feasibility is .

Delete: If we delete the station from the current solution, the loading and unloading quantities of the station in the current solution are , set . After performing the delete move, let the new route of the new solution be .

Let be the set of stations in the route. According to constraints (3) and (5), the total quantities of the route reduce to ; in order to maintain the quantity balance, we need to increase at another station in . First, we put all feasible stations that can adjust the quantity to obtain a feasible solution by increasing into the set . Then, we choose one station from according to the predetermined rules. If we choose station , then set . The characteristics of the stations in must first be analyzed, and then the rules to choose the station to generate the feasible neighbor solution will be proposed.

, , let be a subsequence of a route . The station can be put into only if all the following conditions are satisfied:

There are two rules to choose the station by adjusting the quantities to generate the feasible new solution: Best Choose Criterion and Random Choose Criterion. The Best Choose Criterion is to choose the station i from set with the least increasing penalty cost. Station i is determined from . The Random Choose Criterion is to randomly choose the station i from set .

Insert: If we insert a station into the current route between and , the loading and unloading quantities of the new station m are , and let . After performing the insert move, set the new route of the new solution as .

Set as the set of stations in the new route . According to constraints (3) and (5), the total quantities of the route increase , and in order to maintain the quantity balance, we need to decrease at another station in . First, we put all feasible stations in , which can adjust the quantity to obtain a feasible solution through decreasing into the set . Then, we choose one station from according to the Best/Random Choose Criterion. The Best Choose Criterion is to choose the station i from set with the least increasing penalty cost. The station i is determined from . The Random Choose Criterion is to randomly choose the station i from set .

If we choose station , then set . 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.

, , let be a subsequence of a route .The station can be put into only if all the following conditions are satisfied:

5. Computational Results

The heuristic algorithm presented in Section 4 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 and Q=20, two repositioning times (in seconds) were used, T=14400 and T=21600, and two times for loading and unloading time (in seconds) were considered as L=U=30, L=U=60. Other input data (i.e., , , and ) 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 (l1 and l₂). 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 (the number of all stations) and Len (the number of stations visited by the repositioning truck). Similar to Ho and Szeto [26] set is the nearest integer function, and . After some preliminary experiments, the results from the algorithm were obtained by setting MaxConsIter_LM=[1.5Len], MaxConsIter_BL =[1.5Len], and MaxConsIter_TS=[1.5Len] when <=100; MaxConsIter_LM=[2Len], MaxConsIter_BL =[2Len], and MaxConsIter_TS=[2Len] when >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.

Table 5 provides some characteristics of the instances considered in this experimental study. In the first column in Table 5, 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 increases. This means that the amount of repositioning work increases as the number of stations increases.

Table 5: Characterization of the numerical instances.

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.2167

Table 6 reports the results obtained by setting Q=20,U=L=60,T=14400, which are under different initial inventory levels of stations. “Best” are the best objective values () obtained by the solution algorithm, “Avg” are the average objective values (), 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:

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. is the best percentage of repositioning completion and is the average percentage of repositioning completion, which are defined as follows:

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 and shown in Table 6 indicate the opposite trend as the runtime. As expected, in all the different workloads, the and decrease as the number of stations increases. When the workload is at a real level, the and are higher than those of the other two levels. The and for the light level are higher than those of the heavy level. However, the runtime does not increase with increasing , as shown in Table 6. The runtime decreases as increases. This may be because an increase in 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).

Light level
Real level
Heavy level

Best
Avg
Cpu
(%)
(%)
Best
Avg
Cpu
(%)
(%)
Best
Avg
Cpu
(%)
(%)

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. Table 7 reports the results obtained under the “real level” workload of each station with Q=20, L=U=60, and T=14400. (%) and (%) 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 and 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 becomes larger, the values of and become smaller. When is 80, the values of and 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).

TS_BL
ILS_BL

Best
Avg
Cpu
(%)
(%)
Best
Avg
Cpu
(%)
(%)

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 Different Q, U, L, and T 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.

Table 8 shows the results comparison between Q=20 and Q=10 under U=L=60, T=14400 with sizes ranging from 20 to 200 stations. and indicate the performance of Q=20 relative to that of Q=10 based on the best and the average PCR, respectively. All negative and 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 capacities Q (U=L=60, T=14400).

=20
=10

Best
Avg
Cpu
(%)
(%)
Best
Avg
Cpu
(%)
(%)

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.4245

Table 9 shows the results comparison between loading/unloading time L=U=60 and loading/unloading time L=U=30 under Q=20, T=14400, with sizes ranging from 20 to 200 stations. and indicate the performance of L=U=60 relative to that of L=U=30 based on the best and the average PCR, respectively. All positive and 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).

=60
=30

Best
Avg
Cpu
(%)
(%)
Best
Avg
Cpu
(%)
(%)

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.1178

Table 10 shows the results comparison between repositioning time T=14400 and T=21600 under U=L=60, Q=20 with sizes ranging from 20 to 200 stations. and indicate the performance of T=14400 relative to that of T=21600 based on the best and the average PCR, respectively. All positive and 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).

=14400
=21600

Best
Avg
Cpu
(%)
(%)
Best
Avg
Cpu
(%)
(%)

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.

Data Availability

The data used in this study can be accessed from the web site http://www.eng.tau.ac.il/~talraviv/Publications/3step%20data.zip.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors acknowledge the National Natural Science Foundation of China (Grant no. 71271220) and the Hengyang Social Science Foundation Project (Grant no. 2017D031).

References

L. Caggiani, R. Camporeale, M. Ottomanelli, and W. Y. Szeto, “A modeling framework for the dynamic management of free-floating bike-sharing systems,” Transportation Research Part C: Emerging Technologies, vol. 87, pp. 159–182, 2018.

Z. Haider, A. Nikolaev, J. E. Kang, and C. Kwon, “Inventory rebalancing through pricing in public bike sharing systems,” European Journal of Operational Research, vol. 270, no. 1, pp. 103–117, 2018.

M. Benchimol, P. Benchimol, B. Chappert et al., “Balancing the stations of a self service “bike hire” system,” RAIRO Operations Research, vol. 45, no. 1, pp. 37–61, 2011.

F. Cruz, A. Subramanian, B. P. Bruck, and M. Iori, “A heuristic algorithm for a single vehicle static bike sharing rebalancing problem,” Computers & Operations Research, vol. 79, pp. 19–33, 2017.

J. H. Lin and T. C. Chou, “A geo-aware and VRP-based public bicycle redistribution system,” International Journal of Vehicular Technology, vol. 2012, 2012.

L. Di Gaspero, A. Rendl, and T. Urli, “A hybrid ACO+CP for balancing bicycle sharing systems,” in Hybrid Metaheuristics, vol. 7919 of Lecture Notes in Computer Science, pp. 198–212, Springer, Berlin, Heidelberg, Germany, 2013.

L. Di Gaspero, A. Rendl, and T. Urli, “Balancing bike sharing systems with constraint programming,” Constraints. An International Journal, vol. 21, no. 2, pp. 318–348, 2016.

G. R. Raidl, B. Hu, M. Rainer-Harbach, and P. Papazek, “Balancing bicycle sharing systems: improving a VNS by efficiently determining optimal loading operations,” in Hybrid Metaheuristics, vol. 7919 of Lecture Notes in Computer Science, pp. 130–143, Springer, Berlin, Heidelberg, Germany, 2013.

M. Rainer-Harbach, P. Papazek, B. Hu, and G. R. Raidl, “Balancing bicycle sharing systems: a variable neighborhood search approach,” in Evolutionary Computation in Combinatorial Optimization, vol. 7832 of Lecture Notes in Comput. Sci., pp. 121–132, Springer, Heidelberg, Germany, 2013.

M. Rainer-Harbach, P. Papazek, G. . Raidl, B. Hu, and C. Kloimüllner, “PILOT, GRASP, and VNS approaches for the static balancing of bicycle sharing systems,” Journal of Global Optimization, vol. 63, no. 3, pp. 597–629, 2015.

L. Di Gaspero, A. Rendl, and T. Urli, “Constraint-based approaches for balancing bike sharing systems,” in Principles and Practice of Constraint Programming, vol. 8124 of Lecture Notes in Computer Science, pp. 758–773, Springer, Berlin Heidelberg, Germany, 2013.

P. Angeloudis, J. Hu, and M. G. H. Bell, “A strategic repositioning algorithm for bicycle-sharing schemes,” Transportmetrica A: Transport Science, vol. 10, no. 8, pp. 759–774, 2014.

M. Dell'Amico, M. Iori, S. Novellani, and T. Stützle, “A destroy and repair algorithm for the bike sharing rebalancing problem,” Computers & Operations Research, vol. 71, pp. 149–162, 2016.

M. Dell'Amico, E. Hadjicostantinou, M. Iori, and S. Novellani, “The bike sharing rebalancing problem: Mathematical formulations and benchmark instances,” OMEGA - The International Journal of Management Science, vol. 45, pp. 7–19, 2014.

I. A. Forma, T. Raviv, and M. Tzur, “A 3-step math heuristic for the static repositioning problem in bike-sharing systems,” Transportation Research Part B: Methodological, vol. 71, pp. 230–247, 2015.

R. Alvarez-Valdes, J. M. Belenguer, E. Benavent et al., “Optimizing the level of service quality of a bike-sharing system,” OMEGA - The International Journal of Management Science, vol. 62, pp. 163–175, 2016.

S. C. Ho and W. Y. Szeto, “A hybrid large neighborhood search for the static multi-vehicle bike-repositioning problem,” Transportation Research Part B: Methodological, vol. 95, pp. 340–363, 2017.

J. Schuijbroek, R. C. Hampshire, and W.-J. van Hoeve, “Inventory rebalancing and vehicle routing in bike sharing systems,” European Journal of Operational Research, vol. 257, no. 3, pp. 992–1004, 2017, http://www.sciencedirect.com/science/article/pii/S0377221716306658.

A. Pal and Y. Zhang, “Free-floating bike sharing: Solving real-life large-scale static rebalancing problems,” Transportation Research Part C: Emerging Technologies, vol. 80, pp. 92–116, 2017.

T. Bulhões, A. Subramanian, G. Erdoğan, and G. Laporte, “The static bike relocation problem with multiple vehicles and visits,” European Journal of Operational Research, vol. 264, no. 2, pp. 508–523, 2018.

W. Y. Szeto and C. S. Shui, “Exact loading and unloading strategies for the static multi-vehicle bike repositioning problem,” Transportation Research Part B: Methodological, vol. 109, pp. 176–211, 2018.

Y. Liu, W. Y. Szeto, and S. C. Ho, “A static free-floating bike repositioning problem with multiple heterogeneous vehicles, multiple depots, and multiple visits,” Transportation Research Part C: Emerging Technologies, vol. 92, pp. 208–242, 2018.

D. Chemla, F. Meunier, and R. W. Calvo, “Bike sharing systems: solving the static rebalancing problem,” Discrete Optimization, vol. 10, no. 2, pp. 120–146, 2013.

R. Nair, E. Miller-Hooks, R. C. Hampshire, and A. Bušić, “Large-scale vehicle sharing systems: analysis of Vélib',” International Journal of Sustainable Transportation, vol. 7, no. 1, pp. 85–106, 2012.

G. Erdoğan, G. Laporte, and R. Wolfler Calvo, “The static bicycle relocation problem with demand intervals,” European Journal of Operational Research, vol. 238, no. 2, pp. 451–457, 2014.

S. C. Ho and W. Y. Szeto, “Solving a static repositioning problem in bike-sharing systems using iterated tabu search,” Transportation Research Part E: Logistics and Transportation Review, vol. 69, pp. 180–198, 2014.

G. Erdoğan, M. Battarra, and R. Wolfler Calvo, “An exact algorithm for the static rebalancing problem arising in bicycle sharing systems,” European Journal of Operational Research, vol. 245, no. 3, pp. 667–679, 2015.

A. A. Kadri, I. Kacem, and K. Labadi, “A branch-and-bound algorithm for solving the static rebalancing problem in bicycle-sharing systems,” Computers & Industrial Engineering, vol. 95, pp. 41–52, 2016.

Y. Li, W. Y. Szeto, J. Long, and C. S. Shui, “A multiple type bike repositioning problem,” Transportation Research Part B: Methodological, vol. 90, pp. 263–278, 2016.

W. Y. Szeto, Y. Liu, and S. C. Ho, “Chemical reaction optimization for solving a static bike repositioning problem,” Transportation Research Part D: Transport and Environment, vol. 47, pp. 104–135, 2016.

S. N. Parragh, K. F. Doerner, and R. F. Hartl, “A survey on pickup and delivery problems,” Journal für Betriebswirtschaft, vol. 58, no. 2, pp. 81–117, 2008.

H. Hernández-Pérez and J. Salazar-González, “Heuristics for the one-commodity pickup-and-delivery traveling salesman problem,” Transportation Science, vol. 38, no. 2, pp. 245–255, 2004.

T. Raviv, M. Tzur, and I. A. Forma, “Static repositioning in a bike-sharing system: models and solution approaches,” EURO Journal on Transportation and Logistics, vol. 2, no. 3, pp. 187–229, 2013.

C. Contardo, C. Morency, and L.-M. Rousseau, “Balancing a dynamic public bike-sharing system,” Technical Report, CIRRELT-2012-09, 2012.

L. Caggiani and M. Ottomanelli, “A Modular Soft Computing based Method for Vehicles Repositioning in Bike-sharing Systems,” Procedia - Social and Behavioral Sciences, vol. 54, pp. 675–684, 2012.

D. Zhang, C. Yu, J. Desai, H. Y. K. Lau, and S. Srivathsan, “A time-space network flow approach to dynamic repositioning in bicycle sharing systems,” Transportation Research Part B: Methodological, vol. 103, pp. 188–207, 2017.

C. S. Shui and W. Y. Szeto, “Dynamic green bike repositioning problem – A hybrid rolling horizon artificial bee colony algorithm approach,” Transportation Research Part D: Transport and Environment, vol. 60, pp. 119–136, 2018.

C. E. Miller, A. W. Tucker, and R. A. Zemlin, “Integer programming formulation of traveling salesman problems,” Journal of the ACM, vol. 7, pp. 326–329, 1960.

J. T. Moore and J. F. Bard, “The mixed integer linear bilevel programming problem,” Operations Research, vol. 38, no. 5, pp. 911–921, 1990.

H. R. Lourenço, O. C. Martin, and T. Stützle, “A beginner’s introduction to iterated local search,” in Proceedings of the 4th Metaheuristics International Conference, pp. 16–20, Porto, Portugal, 2001.

H. Lourenço, O. Martin, and T. Stützle, Iterated Local Search, F. Glover and G. Kochenberger, Eds., Handbook of Metaheuristics, Boston, MA, USA, 2003.