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, 3437].

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 [512, 1422], it is more realistic to consider a single vehicle BRP for a part of the city BSS [3, 4, 2330]. 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.

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 [35, 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) [810], LNS (large neighborhood search) [7, 11, 17], ILS (iterated local search) [4, 20], and GRASP (greedy randomized adaptive search procedure) [10].

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.

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.

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.

s0GenerateInitialSolution
sLocalSearch(s0)
while  termination criterion is not met  do
  Perturbation(s,history);
  LocalSearch();
  sAcceptanceCriterion(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), 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 x0, xn+1,and x0=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), d0=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.

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.

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)
rs; rs; 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
rPerturbation(r, Perb_LW);
end
Return  r
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 s0 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)
s0Initialization;
ss0; ConsIter_BL≔0;
while  ConsIter_BL<=MaxConsIter_BL  do
   RVND(s0, LS_BL);
   sILS_LM(,LS_LM, Perb_LM);
if TC(s )<TC(s),then
  ss;
 ConsIter _BL≔0;
else
 ConsIter _BL≔ ConsIter _BL +1;
end
  s0Perturbation(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 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 l1.

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

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)
s0Initialization; Tuba=;
ss0; ss0; 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.
sILS_LM(, LS_LM, Perb_LM);
  if TC(s)<TC(s),then
   ss;
   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 , 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 [email protected] 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 l2). 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 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.

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.

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

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