# An Improved Iterated Local Search Algorithm for the Static Partial Repositioning Problem in Bike-Sharing System

Authors: Qiong Tang; Zhuo Fu; Dezhi Zhang; Meng Qiu; Minyi Li
Journal: Journal of Advanced Transportation
 (2020)

---

## Abstract

In this paper, a single-vehicle static partial repositioning problem (SPRP) is investigated, which distinguishes the user dissatisfaction generated by different stations. The overall objective of the SPRP is to minimize the weighted sum of the total operational time and the total absolute deviation from the target number of bikes at all stations. An iterated local search is developed to solve this problem. A novel loading and unloading quantity adjustment operator is proposed to further improve the quality of the solution. Experiments are conducted on a set of instances from 30 to 300 stations to demonstrate the effectiveness of the proposed customized solution algorithm as well as the adjustment operator. Using a small example, this paper also reveals that the unit penalty cost has an effect on the repositioning strategies.


---

## Body

## 1. Introduction

As urban traffic congestion and environmental pollution problems become more and more serious, bike-sharing systems (BSSs), as a means of sustainable transportation, can effectively solve the “last mile” problem of urban public transportation and provide benefits for users in terms of finances, health, and a low-carbon lifestyle [1]. However, when the BSS is in operation, the demand and supply of bikes in each station often become unbalanced: at some stations, users might not be able to find available bikes, and in some others, there are not enough lockers for users to return their bikes. As a result, user satisfaction decreases in such situations. As pointed out by Laporte et al. [2], several related problems arising at the strategic, tactical, and operational levels can be studied to improve the BSS. Among those, one of the well-known, yet challenging program is the bike rebalancing problem or bike repositioning problem (BRP). Solutions to BRP aim to determine the optimal vehicle routes and the loading and unloading quantities at stations in order to minimize the sum of travel cost and user dissatisfaction subject to operational constraints. As a common practice, the BSS operators redistribute bikes among unbalanced stations accordingly using a dedicated fleet of trucks.Existing research on BRP are mainly divided into two categories: dynamic and static. Dynamic BRP considers daytime operations where the user demand is changing over time, while static BRP focuses on nighttime operations in which the temporal user demand change can be ignored. Despite that many of the recent works have focused on dynamic settings, there are certainly practical advantages and importance of investigating the static repositioning problem. For instance, the static repositioning allows the repositioning truck to travel quickly in the city without aggravating the traffic congestion and parking problems. Hence, most existing works that take into account the next day bike inventory arrangement in BSS have focused on static settings [3]. In terms of the number of trucks used, static BRP studies are mainly divided into two categories: single and multiple vehicle problems. However, compared with multiple-vehicle BRP [3–21], it is more realistic to consider a single-vehicle BRP for a particular street in the city [22–32]. In this paper, we focuses on the single-vehicle static BRP.Often in the operation process of static BRP, considering various constraints of resource including operation time, capacity of the relocation trucks, and the inherent imbalance between supply and demand of the BSS (i.e., the number of surplus bikes is not equal to the number of deficit bikes at all stations), there is no way to achieve a complete or perfect balance among all stations [3]. Therefore, in practice, only a partial balance among some stations can be achieved. This problem is defined as the static partial repositioning problem (SPRP) [13]. In the SPRP, it is necessary to determine which station needed to be rebalanced, and the loading and unloading quantities at those stations, as well as the route of trucks to minimize the sum of travel cost and user dissatisfaction.Existing models to SPRP define user dissatisfaction either as the sum of the unmet demand (only consider stations with too few bikes) [4–9] or as the sum of the absolute differences between the current inventory and the target inventory at all stations (including stations with too many bikes and stations with too few bikes) [10, 11, 32]. They commonly ignore the different level of dissatisfaction caused by different stations. However, existing study demonstrates that user dissatisfaction caused by not satisfying one unit of demand may vary from station to station [12]. For instance, at locations with high density of bike stations or convenient transportation, the level of user dissatisfaction may be low, because users can easily walk to other nearby stations or use other modes of transportation. However, it can be significantly higher in stations where bike sharing is the only option available to users. Therefore, it is important to distinguish the level of user dissatisfaction generated by not satisfying one unit of demand at different stations. It is not very scientific to quantify the total dissatisfaction of all stations by simply using the sum of the number of bikes that are not satisfied at all stations.In this paper, we build an SPRP on the model proposed by Di Gaspero et al. [9] and Szeto et al. [32] and develop customized solution to the SPRP. Different from existing work, our proposed model will distinguish the impact of different stations on user dissatisfaction by defining different unit penalty costs associated with unsatisfied one bike at different stations. To the best of our knowledge, we are among the first to develop customized solutions to SPRP that could take into account the differentiation of user dissatisfaction at different stations.The main contribution of this paper on the SPRP literature can be summarized as follows:(1)
A comprehensive review of the static BRP literature is provided, including the complete and partial repositioning aspects(2)
A new SPRP is introduced, where the level of user dissatisfaction could be different at different stations(3)
A heuristic algorithm based on iterated local search is proposed to solve the problem(4)
A novel loading and unloading quantity adjustment operator is designed to further improve the quality of the solutionThe remainder of this paper is organized as follows: Section2 reviews the BRP literature. Section 3 describes the mathematical formulation of the problem as well as the nonlinear programming model. Then, in Section 4, we present the improved iterated local search algorithm to solve the model. We present experimental results in Section 5 and conclude this paper with discussion on potential future research in Section 6.


## 2. Literature Review

The bike-sharing system has become an increasingly studied problem in the last two decades. Si et al. [33] provided a complete and insightful picture of all 208 relevant articles about bike-sharing research published from 2010 to 2018. As mentioned before, this paper focuses on the static bike repositioning problem (BRP), and hence, this section will mainly analyze the current literature on the static BRP.The research on the static BRP can be traced back to the seminal paper of Benchimol et al. [22], in which the authors introduced the static single-vehicle BRP and proved it to be NP-hard. For the first time, Pal and Zhang [13] categorized the static BRP into complete and partial repositioning problem based on the rigorous of a repositioning that needs to be performed. Table 1 presents a summary of the literature on the static BRP in terms of the type of repositioning (complete or partial), the number of trucks used (|V|), and the objective function.Table 1
Summary of the static BRP in the literature.
ReferencesType of repositioning|V|Objective functionBenchimol et al. [22]Complete repositioningBruck et al. [23]Chemla et al. [24]Cruz et al. [25]Erdoğan et al. [26]Lahoorpoor et al. [27]1Minimize total travel costBulhões et al. [15]Dell’Amico et al. [16, 17]≥1Pal and Zhang [13]≥1Minimize the makespan of the rebalancing fleetWang and Szeto [14]≥1Minimize the total CO2 emissionsRaviv et al. [3]Forma et al. [19]Ho and Szeto [12]Partial repositioning≥1Minimize the weighted sum of the penalty cost and total travel timeHo and Szeto [28]1Minimize total penalty costTang et al. [29]1The upper-level model minimizes the total penalty cost; the lower-level model minimizes the travel timeDi Gaspero et al. [5, 8]Raidl et al. [6]Rainer-Harbach et al. [4, 7]≥1Minimize the weighted sum of the total absolute deviation from the target number of bikes, the total number of loading/unloading quantities, and the total travel timeDi Gaspero et al. [9]≥1Minimize the weighted sum of the total absolute deviation from the target number of bikes, and the total travel timeAlvarez-Valdes et al. [18]≥1Minimizing the overall cost of unsatisfied demandsErdoğan et al. [30]1Minimize total travel cost and handling costLi et al. [31]1Minimize the sum of travel, imbalance, substitution, and occupancy costsLiu et al. [10]≥1Minimize the weighted sum of the total unmet demand, the inconvenience of getting a bike and total operational timeSzeto et al. [32]1Minimize the weighted sum of the total unmet demand and the total operational timeSchuijbroek et al. [20]≥1Minimize maximum the sum of total travel, loading and unloading costsSzeto and Shui [20]≥1Minimize the positive deviation from the tolerance of total demand dissatisfaction first and then service time in form of the total service time or maximum service timeYou [21]≥1Minimize the sum of the total unmet demand and travel costsThis paper1Minimize the weighted sum of the total deviation from the target number of bikes and the total operational time


### 2.1. Complete Repositioning Problem

In the complete repositioning problem, achieving a perfect balance among stations is considered as a hard constraint. In this case, the ideal inventory of each station must be met after the repositioning operation. The objective of this problem is mostly to study how to arrange the order of stations’ visits so that the travel cost is minimized, see, e.g., [15–17, 22–27]. The exceptions are Pal and Zhang [13] and Wang and Szeto [14] where the objective is to minimize the makespan of the rebalancing fleet and the total CO2 emissions, respectively.Various solutions have been proposed to solve the static single-vehicle complete repositioning problem, for instances, Benchimol et al. [22] developed approximation algorithms, Chemla et al. [24] presented a branch-and-cut method embedded Tabu search, and Erdoğan et al. [26] proposed the first exact algorithm that consists of a branch-and-cut algorithm based on Benders combinatorial cuts. Later on, Cruz et al. [25] addressed the same problem by proposing an iterated local search algorithm. Compared with the algorithms proposed by Chemla et al. [24] and Erdoğan et al. [26], the iterated local search developed by Cruz et al. [25] has more advantages in term of solution time and quality. Bruck et al. [23] introduced static single-vehicle BRP with forbidden temporary operations and proposed three exact algorithms based on different mathematical formulations.There are also studies on the static multivehicle complete repositioning problems, for instance, Bulhões et al. [15] presented a branch-and-cut algorithm and an iterated local search. Dell’Amico et al. [16] presented branch-and-cut algorithm. Dell’Amico et al. [17] developed a destroy and repair metaheuristic algorithm, which makes use of effective constructive heuristic and of several local search procedures.It is worth mentioning that there are three BRP articles occurring in the BSS a bit different from most of those complete repositioning problems studied in this survey. Pal and Zhang [13] modeled a static BRP in a free-floating bike-sharing system and presented a hybrid nested large neighborhood search with a variable neighborhood descent algorithm. Wang and Szeto [14] introduced a static green BRP with broken bikes, and applied CPLEX to determine the routing and loading decisions. Lahoorpoor et al. [27] proposed a bottom-up spatial cluster-based method to solve the static BRP with consideration of the users’ behavior in the BSS.


### 2.2. Partial Repositioning Problem

In the partial repositioning problem, achieving a perfect balance among stations is considered as a soft constraint; that is to say, it is not necessary to achieve the ideal inventory at each station after repositioning operation. As can be seen from the Table1, with the exceptions of Ho and Szeto [28], and Alvarez-Valdes et al. [18], the partial repositioning problems in the literature are modeled as multiobjective optimization problems. The goal of these problems is not only to pursue the minimum transportation cost by determining routing decisions, but also to minimize user dissatisfaction by determining loading and unloading decisions at each station. User dissatisfaction is commonly measured by the penalty costs [3, 12, 19, 28, 29], the absolute deviation from the target number of bikes [4–9], and the unmet demand [10, 11, 32].Raviv et al. [3] creatively defined the user dissatisfaction as penalty function which represents the expected number of shortages at any inventory level at each station. The authors addressed a multivehicle SPRP, and proposed a two-phase heuristic to first solve the routing subproblem then to solve the loading and unloading subproblem. Then, Forma et al. [19] presented a 3-step math heuristic, and later on, Ho and Szeto [12] developed a hybrid large neighborhood search algorithm (proven to be more effective than that proposed in [19]) to solve the same problem as in [3]. Ho and Szeto [28] presented an iterated Tabu search to solve a single-vehicle SPRP. Tang et al. [29] proposed a bilevel programming model by taking the outsourcing transportation mode into account and presented an iterated local search algorithm and Tabu search algorithm to solve the optimization problem.In [4], Rainer-Harbach et al. measured the user dissatisfaction based on the deviation from the target number of bikes, and addressed the multivehicle SPRP. They proposed a heuristic approach for the static BRP where the routing decisions are calculated by a variable neighborhood search metaheuristic and the loading decisions are computed by a helper algorithm. There are a list of other solutions proposed for this problem setting, for instance, in [5], Di Gaspero et al. introduced a combination of constraint programming and ant colony optimization approach; in [6], Raidl et al. provided an efficient method for determining optimal loading operations based on two maximum flow computations; in [7], Rainer-Harbach et al. proposed PILOT, a greedy randomized adaptive search procedure, along with a variable neighborhood search algorithm to solve this problem; and in [8], Di Gaspero et al. employed constraint programming to tackle this problem. Last but not least, in [9], Di Gaspero et al. described a multivehicle SPRP with the objective of minimizing both the weighted sum of the total absolute deviation from the target number of bikes, as well as the total travel time. They proposed a constraint programming approach that utilizes a smart branching strategy and large neighborhood search algorithm to solve the problem.For the first time, Szeto et al. [32] measured the user dissatisfaction based on the unmet demand, and addressed a single-vehicle SPRP by proposing an enhanced chemical reaction optimization algorithm. Liu et al. [10] studied a multivehicle SPRP for the free-floating bike-sharing system. They took into account the different convenience level of getting bikes, and proposed an enhanced chemical reaction optimization algorithm to solve the problem. You [11] studied a multivehicle SPRP under a minimum service requirement over a planning horizon, and proposed a two-phase heuristic algorithm to solve the problem based on linear programming.In this section, there are five articles differentiate themselves from the above other studies. Alvarez-Valdes et al. [18] addressed a static BRP to optimize the level of service quality by a two-stage algorithm. Erdoğan et al. [30] addressed a static BRP where the target inventory at each station must lie in a predetermined interval. They developed and implemented a Benders decomposition scheme and a branch-but-cut algorithm to solve the problem. Li et al. [31] investigated a static BRP with multiple bike types, and proposed a combined hybrid genetic algorithm. Schuijbroek et al. [20] unified dual-bounded service level constraints (with added inventory flexibility) and vehicle routing for static BRP and proposed a cluster-first route-second heuristic. Szeto and Shui [21] developed an enhanced artificial bee colony algorithm that incorporated a new set of optimal loading and unloading strategies to solve the static BRP to minimize first the positive deviation from the tolerance of the total demand dissatisfaction and then the total service time of all trucks.As demonstrated above, in complete repositioning problem, the loading and unloading quantities to a station are given in advance, as opposed to our partial repositioning problem where they are decision variable. The increase of decision variables also increased the computational complexity of solving partial positioning problems which is more in line with the operation of the actual BSS. When it comes to partial positioning problems, most researchers have focused their objective on minimizing the user dissatisfaction measured by the sum of unmet demand or deviations from the target fill levels at all stations [4–11, 32]. However, as mentioned previously, the different level of user dissatisfaction caused by unmet demands or deviations at different stations is another important aspect to take into account although it has been ignored in most of the literature. In this paper, we do not only pursue the overall minimization of operational costs and user dissatisfaction costs but also consider the differences in user dissatisfaction at different stations. As summary above, two works have solved the complete repositioning problem by iterated local search [23, 25], only one work has used iterated local search to solve the bilevel partial repositioning problem [29]. It is obvious that the iterated local search algorithm in the above three papers cannot effectively solve the partial repositioning problem studied in this paper. Considering the characteristics and solution quality of the problem, we improve the traditional iterated local search algorithm and propose a customized algorithm for solving the problem in this paper. Computational results show that the improved algorithm is effective.


## 2.1. Complete Repositioning Problem

In the complete repositioning problem, achieving a perfect balance among stations is considered as a hard constraint. In this case, the ideal inventory of each station must be met after the repositioning operation. The objective of this problem is mostly to study how to arrange the order of stations’ visits so that the travel cost is minimized, see, e.g., [15–17, 22–27]. The exceptions are Pal and Zhang [13] and Wang and Szeto [14] where the objective is to minimize the makespan of the rebalancing fleet and the total CO2 emissions, respectively.Various solutions have been proposed to solve the static single-vehicle complete repositioning problem, for instances, Benchimol et al. [22] developed approximation algorithms, Chemla et al. [24] presented a branch-and-cut method embedded Tabu search, and Erdoğan et al. [26] proposed the first exact algorithm that consists of a branch-and-cut algorithm based on Benders combinatorial cuts. Later on, Cruz et al. [25] addressed the same problem by proposing an iterated local search algorithm. Compared with the algorithms proposed by Chemla et al. [24] and Erdoğan et al. [26], the iterated local search developed by Cruz et al. [25] has more advantages in term of solution time and quality. Bruck et al. [23] introduced static single-vehicle BRP with forbidden temporary operations and proposed three exact algorithms based on different mathematical formulations.There are also studies on the static multivehicle complete repositioning problems, for instance, Bulhões et al. [15] presented a branch-and-cut algorithm and an iterated local search. Dell’Amico et al. [16] presented branch-and-cut algorithm. Dell’Amico et al. [17] developed a destroy and repair metaheuristic algorithm, which makes use of effective constructive heuristic and of several local search procedures.It is worth mentioning that there are three BRP articles occurring in the BSS a bit different from most of those complete repositioning problems studied in this survey. Pal and Zhang [13] modeled a static BRP in a free-floating bike-sharing system and presented a hybrid nested large neighborhood search with a variable neighborhood descent algorithm. Wang and Szeto [14] introduced a static green BRP with broken bikes, and applied CPLEX to determine the routing and loading decisions. Lahoorpoor et al. [27] proposed a bottom-up spatial cluster-based method to solve the static BRP with consideration of the users’ behavior in the BSS.


## 2.2. Partial Repositioning Problem

In the partial repositioning problem, achieving a perfect balance among stations is considered as a soft constraint; that is to say, it is not necessary to achieve the ideal inventory at each station after repositioning operation. As can be seen from the Table1, with the exceptions of Ho and Szeto [28], and Alvarez-Valdes et al. [18], the partial repositioning problems in the literature are modeled as multiobjective optimization problems. The goal of these problems is not only to pursue the minimum transportation cost by determining routing decisions, but also to minimize user dissatisfaction by determining loading and unloading decisions at each station. User dissatisfaction is commonly measured by the penalty costs [3, 12, 19, 28, 29], the absolute deviation from the target number of bikes [4–9], and the unmet demand [10, 11, 32].Raviv et al. [3] creatively defined the user dissatisfaction as penalty function which represents the expected number of shortages at any inventory level at each station. The authors addressed a multivehicle SPRP, and proposed a two-phase heuristic to first solve the routing subproblem then to solve the loading and unloading subproblem. Then, Forma et al. [19] presented a 3-step math heuristic, and later on, Ho and Szeto [12] developed a hybrid large neighborhood search algorithm (proven to be more effective than that proposed in [19]) to solve the same problem as in [3]. Ho and Szeto [28] presented an iterated Tabu search to solve a single-vehicle SPRP. Tang et al. [29] proposed a bilevel programming model by taking the outsourcing transportation mode into account and presented an iterated local search algorithm and Tabu search algorithm to solve the optimization problem.In [4], Rainer-Harbach et al. measured the user dissatisfaction based on the deviation from the target number of bikes, and addressed the multivehicle SPRP. They proposed a heuristic approach for the static BRP where the routing decisions are calculated by a variable neighborhood search metaheuristic and the loading decisions are computed by a helper algorithm. There are a list of other solutions proposed for this problem setting, for instance, in [5], Di Gaspero et al. introduced a combination of constraint programming and ant colony optimization approach; in [6], Raidl et al. provided an efficient method for determining optimal loading operations based on two maximum flow computations; in [7], Rainer-Harbach et al. proposed PILOT, a greedy randomized adaptive search procedure, along with a variable neighborhood search algorithm to solve this problem; and in [8], Di Gaspero et al. employed constraint programming to tackle this problem. Last but not least, in [9], Di Gaspero et al. described a multivehicle SPRP with the objective of minimizing both the weighted sum of the total absolute deviation from the target number of bikes, as well as the total travel time. They proposed a constraint programming approach that utilizes a smart branching strategy and large neighborhood search algorithm to solve the problem.For the first time, Szeto et al. [32] measured the user dissatisfaction based on the unmet demand, and addressed a single-vehicle SPRP by proposing an enhanced chemical reaction optimization algorithm. Liu et al. [10] studied a multivehicle SPRP for the free-floating bike-sharing system. They took into account the different convenience level of getting bikes, and proposed an enhanced chemical reaction optimization algorithm to solve the problem. You [11] studied a multivehicle SPRP under a minimum service requirement over a planning horizon, and proposed a two-phase heuristic algorithm to solve the problem based on linear programming.In this section, there are five articles differentiate themselves from the above other studies. Alvarez-Valdes et al. [18] addressed a static BRP to optimize the level of service quality by a two-stage algorithm. Erdoğan et al. [30] addressed a static BRP where the target inventory at each station must lie in a predetermined interval. They developed and implemented a Benders decomposition scheme and a branch-but-cut algorithm to solve the problem. Li et al. [31] investigated a static BRP with multiple bike types, and proposed a combined hybrid genetic algorithm. Schuijbroek et al. [20] unified dual-bounded service level constraints (with added inventory flexibility) and vehicle routing for static BRP and proposed a cluster-first route-second heuristic. Szeto and Shui [21] developed an enhanced artificial bee colony algorithm that incorporated a new set of optimal loading and unloading strategies to solve the static BRP to minimize first the positive deviation from the tolerance of the total demand dissatisfaction and then the total service time of all trucks.As demonstrated above, in complete repositioning problem, the loading and unloading quantities to a station are given in advance, as opposed to our partial repositioning problem where they are decision variable. The increase of decision variables also increased the computational complexity of solving partial positioning problems which is more in line with the operation of the actual BSS. When it comes to partial positioning problems, most researchers have focused their objective on minimizing the user dissatisfaction measured by the sum of unmet demand or deviations from the target fill levels at all stations [4–11, 32]. However, as mentioned previously, the different level of user dissatisfaction caused by unmet demands or deviations at different stations is another important aspect to take into account although it has been ignored in most of the literature. In this paper, we do not only pursue the overall minimization of operational costs and user dissatisfaction costs but also consider the differences in user dissatisfaction at different stations. As summary above, two works have solved the complete repositioning problem by iterated local search [23, 25], only one work has used iterated local search to solve the bilevel partial repositioning problem [29]. It is obvious that the iterated local search algorithm in the above three papers cannot effectively solve the partial repositioning problem studied in this paper. Considering the characteristics and solution quality of the problem, we improve the traditional iterated local search algorithm and propose a customized algorithm for solving the problem in this paper. Computational results show that the improved algorithm is effective.


## 3. Mathematical Formulation

The bike-sharing system studied in this paper comprises a depot, a set of stations, and a truck. Each station has a repositioning demand. If the repositioning demand at a station is greater than 0, the station is defined as a pickup station and can provide bikes to the delivery stations. The set of pickup stations is denoted byP. On the other hand, if the repositioning demand at a station is smaller than 0, then the station is defined as a delivery station and should be supplied with bikes from pickup stations. The set of delivery stations is denoted by D. If the demand at a station is equal to 0, then the station is called a balanced station. Balanced stations can be ignored in the partial repositioning problem as there is no visiting or loading/unloading operations needed.We consider the scenario that only a single truck with a given capacity is available. The truck collects bikes from pickup stations and transports them to delivery stations. It starts from and returns to the depot empty and operates within a given repositioning time constraint. Moreover, each station is allowed to be visited at most once, and not all stations need to be visited. The objective is to determine the route of the truck and the loading and unloading quantities at each visited station to minimize the weighted sum of the total operational time as well as the total absolute deviation from the target number of bikes at all stations.


### 3.1. Notations

#### 3.1.1. Sets

N: set of stationsN0: set of nodes, including the stations (indexed by i, i∈N) and the depot (indexed by 0)P: set of pickup stations, where P⊂ND: set of delivery stations, where D⊂NR: set of stations visited by the truck, where R⊂NU: set of stations not visited by the truck, where U⊂N


#### 3.1.2. Parameters

di: demand of station i at the beginning of the repositioning operation, if di>0, then i∈P; if di<0, then i∈Dtij: operational time that includes the travel time from station i to station j and the expected time required to load or unload bikes at station jQ: truck capacityT: repositioning timeM: a very large numberwi: the unit penalty cost associated with unsatisfied one bike at station i, wi∈0,1,∀i∈Nα: the weight associated with the truck’s total operational time


#### 3.1.3. Decision Variables

xij: binary variable that equal 1 if truck travels directly from station i to station j, and 0 otherwiseyij: the number of bikes onto the truck when it travels directly from station i to station jyiP: the number of bikes loaded on the truck at station i,yiP≥0,∀i∈P, if i∉P, then yiP=0yiD: the number of bikes unloaded from truck at station i,yiD≥0,∀i∈D, if i∉D, then yiD=0


### 3.2. Mathematical Formulation

The formulation for the SPRP addressed in this study is given as follows:(1)minα∑i∈N0∑j∈N0,i≠jtijxij+∑i∈Pwidi−yiP+∑i∈Dwi−di−yiD,which is subject to(2)∑j∈Nx0j=∑j∈Nxj0=1,(3)∑j∈Ny0j=∑j∈Nyj0=0,(4)∑j∈N0,j≠ixij≤1,∀i∈N,(5)∑j∈N0,j≠ixij=∑j∈N0,j≠ixji,∀i∈N,(6)∑i,j∈N0,i≠jtijxij≤T,(7)qj≥qi+1−M1−xij,∀i,j∈N0,i≠j,(8)∑j∈No,j≠iyij−∑j∈No,j≠iyji=yiP−yiD,∀i∈N,(9)yij≤Qxij,∀i,j∈N0,i≠j,(10)∑i∈NyiP−∑i∈NyiD=0,(11)yiD≤minmax−di,0⋅∑j∈N0,j≠ixij,∑j∈N0,j≠iyji,∀i∈N,(12)yiP≤minmaxdi,0⋅∑j∈N0,j≠ixij,Q−∑j∈N0,j≠iyji,∀i∈N,(13)qi≥0,∀i∈N0,(14)yiP≥0,∀i∈N,(15)yiD≥0,∀i∈N,(16)yij≥0,∀i,j∈N0,i≠j,(17)xij∈0,1,∀i,j∈N0,i≠j.Objective function (1) is defined as the weighted sum of the total operational time and the total weighted sum of deviation from the target number of bikes. Constraint (2) states that the truck departs from the depot only once and returns to the depot. Constraint (3) ensures that the truck starts from and returns to the depot empty. Constraints (4) ensure that the truck can visit a station at most once. Constraints (5) ensure that if the truck visits a station, it must leave that station. Constraint (6) ensures that total operational time of the truck does not exceed the repositioning time available. Constraints (7) are the subtour elimination constraint [34]. Constraints (8) require that the number of bikes loaded onto or unloaded from the truck at a given station equals the difference between the truck load before and after the station visit. Constraints (9) ensure that the load on the vehicle cannot be greater than the truck capacity. Constraint (10) ensures that bikes picked up from pickup stations by the truck are eventually delivered to delivery stations. Constraints (11) require that the unloading quantities at a station should not be greater than the number of bikes needed by that station and greater than the number of bikes on the truck when the vehicle reaches station. Constraints (12) require that the loading quantities at a station should not be greater than the number of bikes supplied by that station and greater than the available spaces on the truck when it arrives at the station. Constraints (13) ensure the auxiliary variable is nonnegative. Constraints (14)–(16) ensure that loading and unloading quantities at each station and the number of bikes on the truck are nonnegative integers. Constraints (17) define xij to be a binary variable.


## 3.1. Notations

### 3.1.1. Sets

N: set of stationsN0: set of nodes, including the stations (indexed by i, i∈N) and the depot (indexed by 0)P: set of pickup stations, where P⊂ND: set of delivery stations, where D⊂NR: set of stations visited by the truck, where R⊂NU: set of stations not visited by the truck, where U⊂N


### 3.1.2. Parameters

di: demand of station i at the beginning of the repositioning operation, if di>0, then i∈P; if di<0, then i∈Dtij: operational time that includes the travel time from station i to station j and the expected time required to load or unload bikes at station jQ: truck capacityT: repositioning timeM: a very large numberwi: the unit penalty cost associated with unsatisfied one bike at station i, wi∈0,1,∀i∈Nα: the weight associated with the truck’s total operational time


### 3.1.3. Decision Variables

xij: binary variable that equal 1 if truck travels directly from station i to station j, and 0 otherwiseyij: the number of bikes onto the truck when it travels directly from station i to station jyiP: the number of bikes loaded on the truck at station i,yiP≥0,∀i∈P, if i∉P, then yiP=0yiD: the number of bikes unloaded from truck at station i,yiD≥0,∀i∈D, if i∉D, then yiD=0


## 3.1.1. Sets

N: set of stationsN0: set of nodes, including the stations (indexed by i, i∈N) and the depot (indexed by 0)P: set of pickup stations, where P⊂ND: set of delivery stations, where D⊂NR: set of stations visited by the truck, where R⊂NU: set of stations not visited by the truck, where U⊂N


## 3.1.2. Parameters

di: demand of station i at the beginning of the repositioning operation, if di>0, then i∈P; if di<0, then i∈Dtij: operational time that includes the travel time from station i to station j and the expected time required to load or unload bikes at station jQ: truck capacityT: repositioning timeM: a very large numberwi: the unit penalty cost associated with unsatisfied one bike at station i, wi∈0,1,∀i∈Nα: the weight associated with the truck’s total operational time


## 3.1.3. Decision Variables

xij: binary variable that equal 1 if truck travels directly from station i to station j, and 0 otherwiseyij: the number of bikes onto the truck when it travels directly from station i to station jyiP: the number of bikes loaded on the truck at station i,yiP≥0,∀i∈P, if i∉P, then yiP=0yiD: the number of bikes unloaded from truck at station i,yiD≥0,∀i∈D, if i∉D, then yiD=0


## 3.2. Mathematical Formulation

The formulation for the SPRP addressed in this study is given as follows:(1)minα∑i∈N0∑j∈N0,i≠jtijxij+∑i∈Pwidi−yiP+∑i∈Dwi−di−yiD,which is subject to(2)∑j∈Nx0j=∑j∈Nxj0=1,(3)∑j∈Ny0j=∑j∈Nyj0=0,(4)∑j∈N0,j≠ixij≤1,∀i∈N,(5)∑j∈N0,j≠ixij=∑j∈N0,j≠ixji,∀i∈N,(6)∑i,j∈N0,i≠jtijxij≤T,(7)qj≥qi+1−M1−xij,∀i,j∈N0,i≠j,(8)∑j∈No,j≠iyij−∑j∈No,j≠iyji=yiP−yiD,∀i∈N,(9)yij≤Qxij,∀i,j∈N0,i≠j,(10)∑i∈NyiP−∑i∈NyiD=0,(11)yiD≤minmax−di,0⋅∑j∈N0,j≠ixij,∑j∈N0,j≠iyji,∀i∈N,(12)yiP≤minmaxdi,0⋅∑j∈N0,j≠ixij,Q−∑j∈N0,j≠iyji,∀i∈N,(13)qi≥0,∀i∈N0,(14)yiP≥0,∀i∈N,(15)yiD≥0,∀i∈N,(16)yij≥0,∀i,j∈N0,i≠j,(17)xij∈0,1,∀i,j∈N0,i≠j.Objective function (1) is defined as the weighted sum of the total operational time and the total weighted sum of deviation from the target number of bikes. Constraint (2) states that the truck departs from the depot only once and returns to the depot. Constraint (3) ensures that the truck starts from and returns to the depot empty. Constraints (4) ensure that the truck can visit a station at most once. Constraints (5) ensure that if the truck visits a station, it must leave that station. Constraint (6) ensures that total operational time of the truck does not exceed the repositioning time available. Constraints (7) are the subtour elimination constraint [34]. Constraints (8) require that the number of bikes loaded onto or unloaded from the truck at a given station equals the difference between the truck load before and after the station visit. Constraints (9) ensure that the load on the vehicle cannot be greater than the truck capacity. Constraint (10) ensures that bikes picked up from pickup stations by the truck are eventually delivered to delivery stations. Constraints (11) require that the unloading quantities at a station should not be greater than the number of bikes needed by that station and greater than the number of bikes on the truck when the vehicle reaches station. Constraints (12) require that the loading quantities at a station should not be greater than the number of bikes supplied by that station and greater than the available spaces on the truck when it arrives at the station. Constraints (13) ensure the auxiliary variable is nonnegative. Constraints (14)–(16) ensure that loading and unloading quantities at each station and the number of bikes on the truck are nonnegative integers. Constraints (17) define xij to be a binary variable.


## 4. Solution Method

Similar to other existing works (e.g., [15, 25, 29]), we use iterated local search (ILS) to solve the BRP problem. The ILS algorithm is a metaheuristic algorithm. It obtains the optimal solution through two mechanisms of local search (see Section 4.3) and perturbation search (see Section 4.4). In this paper, local search is executed through the randomized variable neighborhood decent (RVND) procedure. Algorithm 1 illustrates the ILSSPRP algorithm for the SPRP. According to Algorithm 1, after constructing the initial solution s0 by applying the initialization procedure (see Section 4.2), the procedure RVND (see Section 4.3) and perturbation search (see Section 4.4) are called repeatedly until the predefined termination criterion (the maximum consecutive iterations without improvement limit in our case) is satisfied. The local search including 6 neighborhood structures, such as Swap (N1), Insertion (N2), Delete (N3), Exchange (N4), 2-Opt (N5), and Reinsertion (N6) neighborhood structures, is executed using an RVND procedure, whereas the perturbation procedure performs Swap (P1), Insertion (P2), Exchange (P3), and Reinsertion (P4) neighborhood structures.Algorithm 1:ILSSPRP for SPRP.
s0:  = Initialization;s∗: = s0; s: = s0; ConsIter: = 0;whileConsIter ≤ MaxConsIter dos′ ≔ RVND (s, Nk); {1 ≤ k ≤ 6}s ≔ Adjustment (s′);if TC (s) < TC (s∗), thens∗: = s;ConsIter: = 0;elseConsIter: = ConsIter + 1;ends′′: = Perturbation search (s∗, Pm); {1 ≤ m ≤ 4}s: = Adjustment (s′′);EndReturns∗According to the characteristics of the SPRP, we improve the traditional ILS algorithm by adding a novel adjustment operator after RVND and perturbation search. Details of the adjustment operator can be found in Section4.5, and we also verify the impact of the proposed adjustment operator on the ILSSPRP in Section 5.2. Let s represent the current solution for the SPRP and s∗ represent the best solution that can be searched currently. For a solution s, TC (s) gives the objective value of solution s. We also denote ConsIter as the current consecutive iterations without improvement, and MaxConsIter as the maximum consecutive iterations without improvement, respectively.


### 4.1. Solution Representation

A solutionw is represented by w = (x, y, q) as shown in Figure 1, where x is a vector that stores the routing sequence of the truck (x0, x1, ..., xn, xn+1), where x0 = xn + 1 = 0 represents the depot, xhxh∈N,h∈1,…,n represents the station visited by the truck, the subscript h indicates the order of stations being visited by the truck. y is a vector that stores the demand that has been met at the station visited by the truck, where positive and negative values indicate the quantities of bikes collected and delivered, respectively. q is a vector used to store the cumulative load after the truck visits the station. According to equations (3) and (9), if a solution is feasible, then qh=0,∀h∈0,n+1 and 0≤qh≤Q,∀h∈0,1,…,n,n+1 where Q is the truck capacity.Figure 1
Solution representation.


### 4.2. Initialization

In this paper, the initial feasible solution is efficiently generated by a construction heuristic algorithm (see [12, 28, 29]). The basic construction heuristic algorithm was first proposed by Ho and Szeto [28] to solve the single-vehicle BRP, and Ho and Szeto [12] developed it to solve the multiple-vehicle BRP, and then Tang et al. [29] improved it to solve the bilevel BRP.The construction heuristic algorithm can be divided into two main steps:Step one, Sort: we first divide stations into two categories: pickup stations and delivery stations. Then, stations are sorted in descending order by the ranking criteria. Similar to Tang et al. [29], in this paper, the ranking criteria are widiQ, where diQ=minQ,di,i∈P;max−Q,di,i∈D.Step two, Insert: we first insert the pickup stations in turn into the route, if no pickup station can be inserted into the current route without exceeding the capacity limits; then, delivery stations are inserted in turn in the route. We repeat this process until no station can be inserted without violating the predetermined repositioning time constraint.


### 4.3. Local Search

After generating the initial solution, we improve this initial solution with local search technique. We consider two kinds of neighborhood structures: intersets and intraroute. The former performs moves between setR (the set of stations visited by the truck) and set U (the set of stations not visited by the truck), whereas the latter only considers moves within the route.


#### 4.3.1. Interset Neighborhood Structures

(i)
Swap-N1: randomly select a station i from the route (set R), then select another station jdi≤dj of the same type from set U, and exchange them.(ii)
Insertion-N2: randomly select a station from the set U to insert into the route (set R).(iii)
Delete-N3: randomly select a station from the route (set R) and delete it.


#### 4.3.2. Intraroute Neighborhood Structures

(i)
Exchange-N4: randomly select two stations from the route and exchange them.(ii)
2-Opt-N5: randomly select two stations from the route and exchange the access order of the stations between the two stations.(iii)
Reinsertion-N6: randomly select a station and reinsert it into the route.(iv)
The new neighborhood solutions obtained by applying above moves may no longer remain feasible. In order to ensure that a sufficient number of feasible solutions can be obtained, we selectively repair some infeasible solutions. Here, we repair infeasible solutions obtained by applyingN2 and N3 (details to be discussed in Sections 4.6.1 and 4.6.2). In terms of the other moves, only feasible solutions can be accepted as the new neighborhood solutions.


### 4.4. Perturbation Search

In order to avoid falling into local optimum and to realize global search, perturbation is used to constantly change the current local optimum solution. The following four kinds of moves are used to implement perturbations:(i)
Swap-P1: randomly select a station i from the route, then select another station j of the same type from set U, and exchange them.(ii)
Insertion-P2: randomly select a station from the set U and insert it into the route. For the infeasible neighborhood solution, when using the repair operator to repair, Insertion-P2 randomly selects a feasible solution after repair as the new neighborhood solution, while Insertion-N2 chooses the best feasible solution after repair as the new neighborhood solution.(iii)
Exchange-P3: randomly select two adjacent stations from the route and swap them with two other adjacent stations from the route.(iv)
Reinsertion-P4: randomly select two adjacent stations from the route, delete them, and reinsert them elsewhere on the route.(v)
Here, we repair infeasible solutions obtained by applyingP2 (details to be discussed in Section 4.6.3). For other moves P1, P3, and P4, only feasible solutions are accepted as the new neighborhood solutions.


### 4.5. Adjustment Operator

The basic idea of the adjustment operator is that a solution can be improved by properly adjusting the loading and unloading quantities between stations in a heuristically way [10, 12, 28, 32].Similar to the pickup-delivery cycle proposed by Lei and Ouyang [35], we first define a pickup-delivery pair as a sequence of consecutive visits that starts from a pickup station (preceded by the depot or delivery station) and ends at the corresponding delivery station (followed by the depot or pickup station). Figure 2 illustrates an example in which the truck visits two pickup-delivery pairs (with Q = 10). The first pickup-delivery pair contains one pickup station (station 1) and one delivery station (station 2), and the second pickup-delivery pair also contains one pickup station (station 4) and one delivery station (station 5). The truck starts from the depot with empty, picks up 5 bikes at station 1, and delivers directly to station 2, where 3 bikes are unloaded. Then, 6 bikes are further collected at station 4, delivers 8 bikes at station 5, and finally ends at the depot empty.Figure 2
Illustration of pickup-delivery pair.In this paper, we divide the adjustment operator into two categories: adjustment of the quantities of different types of stations in one pickup-delivery pair and adjustment of the quantities of the same type of stations in different pickup-delivery pairs. We first adjust the quantities of stations within each pickup-delivery pair, and then adjust the quantities of the same type of stations across different pickup-delivery pairs.


#### 4.5.1. Adjustment of the Quantities of Different Types of Stations in One Pickup-Delivery Pair

In this case, on the current route, we adjust the quantities of different type of stations within each pickup-delivery pair, until all pickup-delivery pairs are adjusted. For one pickup-delivery pair, we consider one pickup stationxi and one delivery station xj at a time (i and j indicate the order of stations being visited by the truck, respectively), with their loading and unloading quantities increased r at the same time, which is formulated by r=mindxi−yxiP,dxj−yxjD,Q−qi, subject to capacity and loading constraints. Then, the loading and unloading quantities on station xi and station xj are updated at the same time by yxiP=yxiP+r;yxjD=yxjD+r.As shown in Figure3(a), in the first pickup-delivery pair, pickup station 1 and delivery station 2 are selected (with Q = 10) to adjust the loading and unloading quantities. In this case, the increased loading/unloading quantities r=mind1−y1P,d2−y2D,Q−q1=min7−5,6−3,10−5=2. Then, y1P=y1P+2=5+2=7,y2D=y2D+2=3+2=5. In the second pickup-delivery pair, pickup station 4 and delivery station 5 are selected as shown in Figure 3(b), r = 1; then, y4P=y4P+1=6+1=7,y5D=y5D+1=8+1=9.Figure 3
Neighbor solutions after applying the adjustment operator in one pickup-delivery pair. (a) Adjustment of the loading and unloading quantities between stations in the first pickup-delivery pair. (b) Adjustment of the loading and unloading quantities between stations in the second pickup-delivery pair.
(a)(b)


#### 4.5.2. Adjustment of the Quantities of the Same Type of Stations in Different Pickup-Delivery Pairs

For two stations of the same type, adjusting the loading (for pickup stations) and unloading (for delivery stations) quantities between them can improve the current solution, under the condition of ensuring the feasibility of the solution after the transformation. Each time, two stations of the same type in different pickup-delivery pairs are selected, and the loading and unloading quantities of one station are increased by 1 unit, while that of the other station is reduced by 1 unit, resulting in a decrease in the value of the objective function (i.e., the weighted sum of the operational time and deviation from target). The procedure is repeated until there is no loading and unloading quantity adjustment between stations to improve the current solution.Consider a current feasible solution shown in Figure2. Figures 4(a) and 4(b) show an example where pickup station 1 (in the first pickup-delivery pair) and pickup station 4 (in the second pickup-delivery pair) are selected; then, we shift two unit quantity between them. For instance, in Figure 4(a), we increase the loading quantities at station 1 by two, and reduce the loading quantities at station 4 by two. We do the opposite, and it is shown in Figure 4(b). Under the two transformations, the new solution that could reduce the objective value will be accepted.Figure 4
Neighbor solutions after applying the adjustment operator between two pickup-delivery pairs. (a) Adjustment of the loading quantities between station 1 and station 4. (b) Adjustment of the loading quantities between station 1 and station 4. (c) Adjustment of the unloading quantities between station 2 and station 5. (d) Adjustment of the unloading quantities between station 2 and station 5.
(a)(b)(c)(d)Similarly, as shown in Figures4(c) and 4(d), delivery station 2 in the first pickup-delivery pair and delivery station 5 in the second pickup-delivery pair are selected; then, we shift two unit quantities between them. Figure 4(c) shows the operation where we increase the unloading quantities at station 2 by two, and decrease the unloading quantities at station 5 by two. And Figure 4(d) shows the opposite operation. Similarly, between the two newly generated solutions, we accept the one that could reduce the objective value as the current solution.


### 4.6. Repair Operator

There is a couple of constraints for a solution to be feasible. First of all, the sum of the total satisfied demand of stations visited by the truck on the route is 0 (according to equation (10)). For instance, a feasible solution is shown in Figure 5, suppose Q = 10, and the total satisfied demand of stations visited by the truck is 0, that is, 4−3 + 2 + 2−5 = 0. Therefore, if we want to ensure that the neighbor solution obtained by applying the Insertion and Delete is feasible, we must make sure that the sum of the total satisfied demand of stations visited by the truck is 0. Besides, we also need to make sure that the solution meets the capacity and loading constraints (according to equations (3) and (9)). In the following sections, we describe the repair procedure for infeasible solutions obtained by local search and perturbation search, respectively.Figure 5
A feasible solution.


#### 4.6.1. Repairing Infeasible Solutions after Applying N2

Given a feasible solution, if we insert a station into the solution, the obtain new neighbor solution is certainly not feasible. However, the feasibility of the new neighbor solution can be ensured by adjusting the loading/unloading quantity of a certain station. Assuming that the current satisfied demand of the inserted station isΔI, then it is necessary to reduce ΔI on the basis of the current satisfied demand of a certain station on the route, in order to balance out the sum of the total satisfied demand on the route to be 0. In this case, the new neighbor solution after applying the repair procedure is feasible if the capacity and loading constraints are satisfied. There might be more than one feasible adjustment scheme, under which the solution with the lowest objective value is selected as the new neighbor solution.As shown in Figure6, if we insert station 6 ΔI=y6P=2 into the current feasible solution (as shown in Figure 5), the total of all satisfied demands is 2 but not 0. Then, we need to choose one station from the route and reduce its satisfied demand by ΔI. As shown in Figure 7, there are three feasible neighbor solutions can be obtained by applying the repair procedure, and we choose the neighbor solution with the lowest objective value as the new feasible solution.Figure 6
An infeasible solution with station 6 inserted.Figure 7
Feasible solutions after applying repair operator. (a) Adjustment of the loading quantities of station 1. (b) Adjustment of the unloading quantities of station 2. (c) Adjustment the unloading quantities of station 5.
(a)(b)(c)


#### 4.6.2. Repairing Infeasible Solutions after Applying N3

Assuming that the current satisfied demand for the deleted station isΔD, then it is necessary to increase ΔD on the basis of the current satisfied demand of a certain station on the route, to ensure that the sum of the total satisfied demand on the route is 0. Similarly, if there are more than one feasible neighbor solutions, the solution with the lowest objective value is selected as the new neighbor solution.As shown in Figure8, if we delete station 4 ΔD=y4P=2 from the current feasible solution (Figure 5), the total of all satisfied demand is −2. Therefore, we need to choose one station from the route and increase its satisfied demand by ΔD. As shown in Figure 9, there are three feasible neighbor solutions can be obtained from the repair procedure, and we choose the neighbor solution with the lowest objective value as the new neighbor solution.Figure 8
An infeasible solution with station 4 deleted.Figure 9
Feasible solutions after applying repair operator. (a) Adjustment of the loading quantities of station 1. (b) Adjustment of the unloading quantities of station 2. (c) Adjustment of the unloading quantities of station 5.
(a)(b)(c)


#### 4.6.3. Repairing Infeasible Solutions after Applying P2

In this case, the repair operation would be the same as in Section4.6.1, the only difference is that when there is more than one feasible neighborhood solution generated, we randomly select one as the new neighbor solution.


## 4.1. Solution Representation

A solutionw is represented by w = (x, y, q) as shown in Figure 1, where x is a vector that stores the routing sequence of the truck (x0, x1, ..., xn, xn+1), where x0 = xn + 1 = 0 represents the depot, xhxh∈N,h∈1,…,n represents the station visited by the truck, the subscript h indicates the order of stations being visited by the truck. y is a vector that stores the demand that has been met at the station visited by the truck, where positive and negative values indicate the quantities of bikes collected and delivered, respectively. q is a vector used to store the cumulative load after the truck visits the station. According to equations (3) and (9), if a solution is feasible, then qh=0,∀h∈0,n+1 and 0≤qh≤Q,∀h∈0,1,…,n,n+1 where Q is the truck capacity.Figure 1
Solution representation.


## 4.2. Initialization

In this paper, the initial feasible solution is efficiently generated by a construction heuristic algorithm (see [12, 28, 29]). The basic construction heuristic algorithm was first proposed by Ho and Szeto [28] to solve the single-vehicle BRP, and Ho and Szeto [12] developed it to solve the multiple-vehicle BRP, and then Tang et al. [29] improved it to solve the bilevel BRP.The construction heuristic algorithm can be divided into two main steps:Step one, Sort: we first divide stations into two categories: pickup stations and delivery stations. Then, stations are sorted in descending order by the ranking criteria. Similar to Tang et al. [29], in this paper, the ranking criteria are widiQ, where diQ=minQ,di,i∈P;max−Q,di,i∈D.Step two, Insert: we first insert the pickup stations in turn into the route, if no pickup station can be inserted into the current route without exceeding the capacity limits; then, delivery stations are inserted in turn in the route. We repeat this process until no station can be inserted without violating the predetermined repositioning time constraint.


## 4.3. Local Search

After generating the initial solution, we improve this initial solution with local search technique. We consider two kinds of neighborhood structures: intersets and intraroute. The former performs moves between setR (the set of stations visited by the truck) and set U (the set of stations not visited by the truck), whereas the latter only considers moves within the route.


### 4.3.1. Interset Neighborhood Structures

(i)
Swap-N1: randomly select a station i from the route (set R), then select another station jdi≤dj of the same type from set U, and exchange them.(ii)
Insertion-N2: randomly select a station from the set U to insert into the route (set R).(iii)
Delete-N3: randomly select a station from the route (set R) and delete it.


### 4.3.2. Intraroute Neighborhood Structures

(i)
Exchange-N4: randomly select two stations from the route and exchange them.(ii)
2-Opt-N5: randomly select two stations from the route and exchange the access order of the stations between the two stations.(iii)
Reinsertion-N6: randomly select a station and reinsert it into the route.(iv)
The new neighborhood solutions obtained by applying above moves may no longer remain feasible. In order to ensure that a sufficient number of feasible solutions can be obtained, we selectively repair some infeasible solutions. Here, we repair infeasible solutions obtained by applyingN2 and N3 (details to be discussed in Sections 4.6.1 and 4.6.2). In terms of the other moves, only feasible solutions can be accepted as the new neighborhood solutions.


## 4.3.1. Interset Neighborhood Structures

(i)
Swap-N1: randomly select a station i from the route (set R), then select another station jdi≤dj of the same type from set U, and exchange them.(ii)
Insertion-N2: randomly select a station from the set U to insert into the route (set R).(iii)
Delete-N3: randomly select a station from the route (set R) and delete it.


## 4.3.2. Intraroute Neighborhood Structures

(i)
Exchange-N4: randomly select two stations from the route and exchange them.(ii)
2-Opt-N5: randomly select two stations from the route and exchange the access order of the stations between the two stations.(iii)
Reinsertion-N6: randomly select a station and reinsert it into the route.(iv)
The new neighborhood solutions obtained by applying above moves may no longer remain feasible. In order to ensure that a sufficient number of feasible solutions can be obtained, we selectively repair some infeasible solutions. Here, we repair infeasible solutions obtained by applyingN2 and N3 (details to be discussed in Sections 4.6.1 and 4.6.2). In terms of the other moves, only feasible solutions can be accepted as the new neighborhood solutions.


## 4.4. Perturbation Search

In order to avoid falling into local optimum and to realize global search, perturbation is used to constantly change the current local optimum solution. The following four kinds of moves are used to implement perturbations:(i)
Swap-P1: randomly select a station i from the route, then select another station j of the same type from set U, and exchange them.(ii)
Insertion-P2: randomly select a station from the set U and insert it into the route. For the infeasible neighborhood solution, when using the repair operator to repair, Insertion-P2 randomly selects a feasible solution after repair as the new neighborhood solution, while Insertion-N2 chooses the best feasible solution after repair as the new neighborhood solution.(iii)
Exchange-P3: randomly select two adjacent stations from the route and swap them with two other adjacent stations from the route.(iv)
Reinsertion-P4: randomly select two adjacent stations from the route, delete them, and reinsert them elsewhere on the route.(v)
Here, we repair infeasible solutions obtained by applyingP2 (details to be discussed in Section 4.6.3). For other moves P1, P3, and P4, only feasible solutions are accepted as the new neighborhood solutions.


## 4.5. Adjustment Operator

The basic idea of the adjustment operator is that a solution can be improved by properly adjusting the loading and unloading quantities between stations in a heuristically way [10, 12, 28, 32].Similar to the pickup-delivery cycle proposed by Lei and Ouyang [35], we first define a pickup-delivery pair as a sequence of consecutive visits that starts from a pickup station (preceded by the depot or delivery station) and ends at the corresponding delivery station (followed by the depot or pickup station). Figure 2 illustrates an example in which the truck visits two pickup-delivery pairs (with Q = 10). The first pickup-delivery pair contains one pickup station (station 1) and one delivery station (station 2), and the second pickup-delivery pair also contains one pickup station (station 4) and one delivery station (station 5). The truck starts from the depot with empty, picks up 5 bikes at station 1, and delivers directly to station 2, where 3 bikes are unloaded. Then, 6 bikes are further collected at station 4, delivers 8 bikes at station 5, and finally ends at the depot empty.Figure 2
Illustration of pickup-delivery pair.In this paper, we divide the adjustment operator into two categories: adjustment of the quantities of different types of stations in one pickup-delivery pair and adjustment of the quantities of the same type of stations in different pickup-delivery pairs. We first adjust the quantities of stations within each pickup-delivery pair, and then adjust the quantities of the same type of stations across different pickup-delivery pairs.


### 4.5.1. Adjustment of the Quantities of Different Types of Stations in One Pickup-Delivery Pair

In this case, on the current route, we adjust the quantities of different type of stations within each pickup-delivery pair, until all pickup-delivery pairs are adjusted. For one pickup-delivery pair, we consider one pickup stationxi and one delivery station xj at a time (i and j indicate the order of stations being visited by the truck, respectively), with their loading and unloading quantities increased r at the same time, which is formulated by r=mindxi−yxiP,dxj−yxjD,Q−qi, subject to capacity and loading constraints. Then, the loading and unloading quantities on station xi and station xj are updated at the same time by yxiP=yxiP+r;yxjD=yxjD+r.As shown in Figure3(a), in the first pickup-delivery pair, pickup station 1 and delivery station 2 are selected (with Q = 10) to adjust the loading and unloading quantities. In this case, the increased loading/unloading quantities r=mind1−y1P,d2−y2D,Q−q1=min7−5,6−3,10−5=2. Then, y1P=y1P+2=5+2=7,y2D=y2D+2=3+2=5. In the second pickup-delivery pair, pickup station 4 and delivery station 5 are selected as shown in Figure 3(b), r = 1; then, y4P=y4P+1=6+1=7,y5D=y5D+1=8+1=9.Figure 3
Neighbor solutions after applying the adjustment operator in one pickup-delivery pair. (a) Adjustment of the loading and unloading quantities between stations in the first pickup-delivery pair. (b) Adjustment of the loading and unloading quantities between stations in the second pickup-delivery pair.
(a)(b)


### 4.5.2. Adjustment of the Quantities of the Same Type of Stations in Different Pickup-Delivery Pairs

For two stations of the same type, adjusting the loading (for pickup stations) and unloading (for delivery stations) quantities between them can improve the current solution, under the condition of ensuring the feasibility of the solution after the transformation. Each time, two stations of the same type in different pickup-delivery pairs are selected, and the loading and unloading quantities of one station are increased by 1 unit, while that of the other station is reduced by 1 unit, resulting in a decrease in the value of the objective function (i.e., the weighted sum of the operational time and deviation from target). The procedure is repeated until there is no loading and unloading quantity adjustment between stations to improve the current solution.Consider a current feasible solution shown in Figure2. Figures 4(a) and 4(b) show an example where pickup station 1 (in the first pickup-delivery pair) and pickup station 4 (in the second pickup-delivery pair) are selected; then, we shift two unit quantity between them. For instance, in Figure 4(a), we increase the loading quantities at station 1 by two, and reduce the loading quantities at station 4 by two. We do the opposite, and it is shown in Figure 4(b). Under the two transformations, the new solution that could reduce the objective value will be accepted.Figure 4
Neighbor solutions after applying the adjustment operator between two pickup-delivery pairs. (a) Adjustment of the loading quantities between station 1 and station 4. (b) Adjustment of the loading quantities between station 1 and station 4. (c) Adjustment of the unloading quantities between station 2 and station 5. (d) Adjustment of the unloading quantities between station 2 and station 5.
(a)(b)(c)(d)Similarly, as shown in Figures4(c) and 4(d), delivery station 2 in the first pickup-delivery pair and delivery station 5 in the second pickup-delivery pair are selected; then, we shift two unit quantities between them. Figure 4(c) shows the operation where we increase the unloading quantities at station 2 by two, and decrease the unloading quantities at station 5 by two. And Figure 4(d) shows the opposite operation. Similarly, between the two newly generated solutions, we accept the one that could reduce the objective value as the current solution.


## 4.5.1. Adjustment of the Quantities of Different Types of Stations in One Pickup-Delivery Pair

In this case, on the current route, we adjust the quantities of different type of stations within each pickup-delivery pair, until all pickup-delivery pairs are adjusted. For one pickup-delivery pair, we consider one pickup stationxi and one delivery station xj at a time (i and j indicate the order of stations being visited by the truck, respectively), with their loading and unloading quantities increased r at the same time, which is formulated by r=mindxi−yxiP,dxj−yxjD,Q−qi, subject to capacity and loading constraints. Then, the loading and unloading quantities on station xi and station xj are updated at the same time by yxiP=yxiP+r;yxjD=yxjD+r.As shown in Figure3(a), in the first pickup-delivery pair, pickup station 1 and delivery station 2 are selected (with Q = 10) to adjust the loading and unloading quantities. In this case, the increased loading/unloading quantities r=mind1−y1P,d2−y2D,Q−q1=min7−5,6−3,10−5=2. Then, y1P=y1P+2=5+2=7,y2D=y2D+2=3+2=5. In the second pickup-delivery pair, pickup station 4 and delivery station 5 are selected as shown in Figure 3(b), r = 1; then, y4P=y4P+1=6+1=7,y5D=y5D+1=8+1=9.Figure 3
Neighbor solutions after applying the adjustment operator in one pickup-delivery pair. (a) Adjustment of the loading and unloading quantities between stations in the first pickup-delivery pair. (b) Adjustment of the loading and unloading quantities between stations in the second pickup-delivery pair.
(a)(b)


## 4.5.2. Adjustment of the Quantities of the Same Type of Stations in Different Pickup-Delivery Pairs

For two stations of the same type, adjusting the loading (for pickup stations) and unloading (for delivery stations) quantities between them can improve the current solution, under the condition of ensuring the feasibility of the solution after the transformation. Each time, two stations of the same type in different pickup-delivery pairs are selected, and the loading and unloading quantities of one station are increased by 1 unit, while that of the other station is reduced by 1 unit, resulting in a decrease in the value of the objective function (i.e., the weighted sum of the operational time and deviation from target). The procedure is repeated until there is no loading and unloading quantity adjustment between stations to improve the current solution.Consider a current feasible solution shown in Figure2. Figures 4(a) and 4(b) show an example where pickup station 1 (in the first pickup-delivery pair) and pickup station 4 (in the second pickup-delivery pair) are selected; then, we shift two unit quantity between them. For instance, in Figure 4(a), we increase the loading quantities at station 1 by two, and reduce the loading quantities at station 4 by two. We do the opposite, and it is shown in Figure 4(b). Under the two transformations, the new solution that could reduce the objective value will be accepted.Figure 4
Neighbor solutions after applying the adjustment operator between two pickup-delivery pairs. (a) Adjustment of the loading quantities between station 1 and station 4. (b) Adjustment of the loading quantities between station 1 and station 4. (c) Adjustment of the unloading quantities between station 2 and station 5. (d) Adjustment of the unloading quantities between station 2 and station 5.
(a)(b)(c)(d)Similarly, as shown in Figures4(c) and 4(d), delivery station 2 in the first pickup-delivery pair and delivery station 5 in the second pickup-delivery pair are selected; then, we shift two unit quantities between them. Figure 4(c) shows the operation where we increase the unloading quantities at station 2 by two, and decrease the unloading quantities at station 5 by two. And Figure 4(d) shows the opposite operation. Similarly, between the two newly generated solutions, we accept the one that could reduce the objective value as the current solution.


## 4.6. Repair Operator

There is a couple of constraints for a solution to be feasible. First of all, the sum of the total satisfied demand of stations visited by the truck on the route is 0 (according to equation (10)). For instance, a feasible solution is shown in Figure 5, suppose Q = 10, and the total satisfied demand of stations visited by the truck is 0, that is, 4−3 + 2 + 2−5 = 0. Therefore, if we want to ensure that the neighbor solution obtained by applying the Insertion and Delete is feasible, we must make sure that the sum of the total satisfied demand of stations visited by the truck is 0. Besides, we also need to make sure that the solution meets the capacity and loading constraints (according to equations (3) and (9)). In the following sections, we describe the repair procedure for infeasible solutions obtained by local search and perturbation search, respectively.Figure 5
A feasible solution.


### 4.6.1. Repairing Infeasible Solutions after Applying N2

Given a feasible solution, if we insert a station into the solution, the obtain new neighbor solution is certainly not feasible. However, the feasibility of the new neighbor solution can be ensured by adjusting the loading/unloading quantity of a certain station. Assuming that the current satisfied demand of the inserted station isΔI, then it is necessary to reduce ΔI on the basis of the current satisfied demand of a certain station on the route, in order to balance out the sum of the total satisfied demand on the route to be 0. In this case, the new neighbor solution after applying the repair procedure is feasible if the capacity and loading constraints are satisfied. There might be more than one feasible adjustment scheme, under which the solution with the lowest objective value is selected as the new neighbor solution.As shown in Figure6, if we insert station 6 ΔI=y6P=2 into the current feasible solution (as shown in Figure 5), the total of all satisfied demands is 2 but not 0. Then, we need to choose one station from the route and reduce its satisfied demand by ΔI. As shown in Figure 7, there are three feasible neighbor solutions can be obtained by applying the repair procedure, and we choose the neighbor solution with the lowest objective value as the new feasible solution.Figure 6
An infeasible solution with station 6 inserted.Figure 7
Feasible solutions after applying repair operator. (a) Adjustment of the loading quantities of station 1. (b) Adjustment of the unloading quantities of station 2. (c) Adjustment the unloading quantities of station 5.
(a)(b)(c)


### 4.6.2. Repairing Infeasible Solutions after Applying N3

Assuming that the current satisfied demand for the deleted station isΔD, then it is necessary to increase ΔD on the basis of the current satisfied demand of a certain station on the route, to ensure that the sum of the total satisfied demand on the route is 0. Similarly, if there are more than one feasible neighbor solutions, the solution with the lowest objective value is selected as the new neighbor solution.As shown in Figure8, if we delete station 4 ΔD=y4P=2 from the current feasible solution (Figure 5), the total of all satisfied demand is −2. Therefore, we need to choose one station from the route and increase its satisfied demand by ΔD. As shown in Figure 9, there are three feasible neighbor solutions can be obtained from the repair procedure, and we choose the neighbor solution with the lowest objective value as the new neighbor solution.Figure 8
An infeasible solution with station 4 deleted.Figure 9
Feasible solutions after applying repair operator. (a) Adjustment of the loading quantities of station 1. (b) Adjustment of the unloading quantities of station 2. (c) Adjustment of the unloading quantities of station 5.
(a)(b)(c)


### 4.6.3. Repairing Infeasible Solutions after Applying P2

In this case, the repair operation would be the same as in Section4.6.1, the only difference is that when there is more than one feasible neighborhood solution generated, we randomly select one as the new neighbor solution.


## 4.6.1. Repairing Infeasible Solutions after Applying N2

Given a feasible solution, if we insert a station into the solution, the obtain new neighbor solution is certainly not feasible. However, the feasibility of the new neighbor solution can be ensured by adjusting the loading/unloading quantity of a certain station. Assuming that the current satisfied demand of the inserted station isΔI, then it is necessary to reduce ΔI on the basis of the current satisfied demand of a certain station on the route, in order to balance out the sum of the total satisfied demand on the route to be 0. In this case, the new neighbor solution after applying the repair procedure is feasible if the capacity and loading constraints are satisfied. There might be more than one feasible adjustment scheme, under which the solution with the lowest objective value is selected as the new neighbor solution.As shown in Figure6, if we insert station 6 ΔI=y6P=2 into the current feasible solution (as shown in Figure 5), the total of all satisfied demands is 2 but not 0. Then, we need to choose one station from the route and reduce its satisfied demand by ΔI. As shown in Figure 7, there are three feasible neighbor solutions can be obtained by applying the repair procedure, and we choose the neighbor solution with the lowest objective value as the new feasible solution.Figure 6
An infeasible solution with station 6 inserted.Figure 7
Feasible solutions after applying repair operator. (a) Adjustment of the loading quantities of station 1. (b) Adjustment of the unloading quantities of station 2. (c) Adjustment the unloading quantities of station 5.
(a)(b)(c)


## 4.6.2. Repairing Infeasible Solutions after Applying N3

Assuming that the current satisfied demand for the deleted station isΔD, then it is necessary to increase ΔD on the basis of the current satisfied demand of a certain station on the route, to ensure that the sum of the total satisfied demand on the route is 0. Similarly, if there are more than one feasible neighbor solutions, the solution with the lowest objective value is selected as the new neighbor solution.As shown in Figure8, if we delete station 4 ΔD=y4P=2 from the current feasible solution (Figure 5), the total of all satisfied demand is −2. Therefore, we need to choose one station from the route and increase its satisfied demand by ΔD. As shown in Figure 9, there are three feasible neighbor solutions can be obtained from the repair procedure, and we choose the neighbor solution with the lowest objective value as the new neighbor solution.Figure 8
An infeasible solution with station 4 deleted.Figure 9
Feasible solutions after applying repair operator. (a) Adjustment of the loading quantities of station 1. (b) Adjustment of the unloading quantities of station 2. (c) Adjustment of the unloading quantities of station 5.
(a)(b)(c)


## 4.6.3. Repairing Infeasible Solutions after Applying P2

In this case, the repair operation would be the same as in Section4.6.1, the only difference is that when there is more than one feasible neighborhood solution generated, we randomly select one as the new neighbor solution.


## 5. Computational Results

In this section, we present experiments conducted to verify the effectiveness of the proposed ILSSPRP algorithm as well as the proposed adjustment operator. We compare the performance of our proposed solution with the optimal solution of Lingo 18. We also analyze through experiments the influence of the change of weight coefficient α on the proposed ILSSPRP and the change of unit penalty cost wi on the repositioning strategies.We follow the instances proposed by Rain-Harbach et al. [7] to conduct the experiments, where the size of the instances is ranging from 30 to 300 stations. For each instance, we select the first seed of the 30 independent seeds. The initial objective value, the total pickup, and total drop-off quantities of each instances are listed in Table 2. We consider instances with 30 and 60 stations as the small-size instances; instances contain 120 and 180 stations as the medium-size instances; and instances with 240 and 300 stations as the large-size instances. In this study, unless otherwise specified, we follow the setting in Szeto et al. [32], α = 0.00001. Three repositioning time constraints (in seconds) were used, T = 14400 s, T = 21600 s, and T = 28800 s. Two truck capacities (in terms of the number of bikes) were considered: Q = 10 and Q = 20, respectively. The unit penalty cost wi associated with one unsatisfied bike at station i was randomly generated in the interval (0, 1], ∀i∈N.Table 2
Characterization of the instances.
|N|Initial objective valueTotal pickup quantitiesTotal drop-off quantities3099.210710760170.3179179120330.7354354180606.4557557240760.4785785300852.2883883The experiments are conducted in Matlab, on a computer with Intel Core i5-4590 CPU @ 3.30 GHz and 4 GB RAM. We run each instance for 20 times, the best and average of the solutions and the average computing time of those 20 runs are used to evaluate the performance of the ILSSPRP.The main parameter of ILSSPRP to be calibrated is the maximum consecutive iterations without improvement (MaxConsIter). Here, we set its value as a function of the instance, MaxConsIter = max{100, |N|}, similar to Bulhões et al. [15].


### 5.1. Comparison between Lingo and ILSSPRP

In terms of the quality of the solution and the calculation efficiency, Tables3–5 show the comparison results obtained by Lingo and ILSSPRP under different repositioning times (T = 14400 s, T = 21600 s, and T = 28800 s), truck capacity (Q = 10 and Q = 20), and number of stations (|N|). GapBest and GapAvg (%) indicate the performance of the ILSSPRP relative to that of Lingo based on the best and average objective values, respectively. The CPU (in seconds) indicates the average computing time of the ILSSPRP or the computing time for Lingo, respectively.Table 3
Result comparison of Lingo and ILSSPRP (T = 14400 s).
InstanceLingoILSSPRPGapBest (%)GapAvg (%)|N|QOptimal valueCPUBest valueAverage valueCPU301048.644023348.644048.64404.270.000.002043.5440360041.438641.43864.58−5.08−5.086010117.14203600117.6428117.74245.330.430.5120113.44403600111.6422112.98244.86−1.61−0.4112010273.74403600271.2440271.91016.65−0.92−0.6720296.81803600256.9386257.01926.74−15.52−15.4818010554.14403600538.4374538.72098.02−2.92−2.8620542.04403600528.2434528.84308.11−2.61−2.5024010690.0386690.03868.4720671.9338671.93388.9130010796.3374796.90978.9220789.1392792.123610.63Avg3179.137.12−3.53−3.31Table 4
Result comparison of Lingo and ILSSPRP (T = 21600 s).
InstanceLingoILSSPRPGapBest (%)GapAvg (%)|N|QOptimal valueCPUBest valueAverage valueCPU301028.8148360028.517828.70408.32−1.04−0.392023.5148360021.213021.73989.42−10.85−8.16601094.2160360093.213093.72758.12−1.08−0.522084.4150360081.722683.947412.92−3.29−0.5612010249.31603600243.5154244.153115.46−2.38−2.1120229.41543600222.4112224.314214.45−3.15−2.2718010513.21603600506.4052506.405211.77−1.34−1.3420594.46103600493.0178494.476122.22−20.58−20.2224010656.0130656.462217.6920629.4130630.375732.6930010764.4004764.400517.2220743.1136744.145046.51Avg360018.07−5.46−4.45Table 5
Result comparison of Lingo and ILSSPRP (T = 28800 s).
InstanceLingoILSSPRPGapBest (%)GapAvg (%)|N|QOptimal valueCPUBest valueAverage valueCPU301020.0888360014.286215.344614.00−40.62−30.922013.586236008.68509.192818.10−56.43−47.79601079.5874360074.780875.396716.96−6.43−5.562065.7868360063.588664.140719.10−3.46−2.5712010223.88803600220.4964221.820524.22−1.54−0.9320243.31663600193.2868195.197833.08−25.88−24.6518010593.44323600475.4784476.172132.82−24.81−24.6320588.86603600454.5760454.992036.97−29.54−29.4224010623.5814625.904041.9120582.9796583.620756.0530010719.6832720.362153.9820711.2862712.261557.28Avg360033.71−23.59−20.81From Tables3–5, we can see that Lingo is not able to obtain any feasible solutions for the large-size problems (where |N| = 240 and |N| = 300). It also generally fails to obtain any global optimal solutions within 3600 seconds even for small-size problems. In contrast, for the various size of problems from small to large, the average computing time of the ILSSPRP varies from 7.12 s to 33.71 s.As shown in the Tables3–5, almost all of the values of GapBest and GapAvg are negative (with only one exception for the results obtained under |N| = 60, Q = 10, T = 14400 s). The absolute values of the average GapBest and GapAvg values obtained by the ILSSPRP for the scenarios with long repositioning time are larger than those with short repositioning time. GapBest and GapAvg are, respectively, −3.53 and −3.31 for T = 14400 s, and −5.46 and −4.45 for T = 21600 s, but −23.59 and −20.81 for T = 28800 s. This implies that the ILSSPRP performs better than Lingo when solving the SPRP with a longer repositioning time.Therefore, we can draw the conclusion that the proposed ILSSPRP can obtain better feasible solutions much faster than Lingo.


### 5.2. The Effectiveness of Introducing Adjustment Operator into ILSSPRP

In this study, the adjustment operator (presented in Section4.5), which considers the characteristics of the SPRP, is introduced to improve the solution quality obtained by the proposed ILSSPRP. In this section, we study through experiments the effectiveness of the proposed adjustment operator by comparing the performance of the ILSSPRP with and without adjustment operator using the various instances described in Section 5.1.Tables6–8 present the results obtained by the proposed ILSSPRP algorithm with and without the adjustment operator, under different repositioning times, respectively. CPU denotes the average computing time in seconds. GapBest and GapAvg here are the deviations of the best/average of the objective values from 20 runs obtained by the ILSSPRP with the adjustment operator from the corresponding optimal solution obtained by ILSSPRP without the adjustment operator, respectively. From Tables 6–8, we can see ILSSPRP without the adjustment operator resulting in a shorter computation time than the ILSSPRP with the adjustment operator.Table 6
Result comparison of ILSSPRP with and without the adjustment operator (T = 14400 s).
InstanceILSSPRP without adjustmentILSSPRPGapBest (%)GapAvg (%)|N|QBest valueAverage valueCPUBest valueAverage valueCPU301051.644051.94402.0648.644048.64404.27−6.17−6.782048.338648.90392.2341.438641.43864.58−16.65−18.026010118.1410118.14101.49117.6428117.74245.33−0.42−0.3420125.8314125.83152.11111.6422112.98244.86−12.71−11.3712010277.6434278.06502.20271.2440271.91016.65−2.36−2.2620262.0434263.84932.42256.9386257.01926.74−1.99−2.6618010551.5374552.83163.86538.4374538.72098.02−2.43−2.6220536.1428537.96504.44528.2434528.84308.11−1.50−1.7224010691.0392691.03923.92690.0386690.03868.47−0.15−0.1520675.6338675.63384.77671.9338671.93388.91−0.55−0.5530010797.3254797.32544.09796.3374796.90978.92−0.12−0.0520794.0284794.02844.48789.1392792.123610.63−0.62−0.24Avg3.177.12−3.81−3.90Table 7
Result comparison of ILSSPRP with and without the adjustment operator (T = 21600 s).
InstanceILSSPRP without adjustmentILSSPRPGapBest (%)GapAvg (%)|N|QBest valueAverage valueCPUBest valueAverage valueCPU301032.911233.46013.3128.517828.70408.32−15.41−16.572027.304628.12554.3821.213021.73989.42−28.72−29.376010103.0052103.87583.5193.213093.72758.12−10.51−10.832095.403495.40444.5981.722683.947412.92−16.74−13.6512010252.6992252.99564.13243.5154244.153115.46−3.77−3.6220234.0136235.92095.58222.4112224.314214.45−5.22−5.1718010509.4076509.40786.83506.4052506.405211.77−0.59−0.5920495.8064496.653710.66493.0178494.476122.22−0.57−0.4424010659.7034659.70347.88656.0130656.462217.69−0.56−0.4920643.8136644.399612.00629.4130630.375732.69−2.29−2.2230010764.5106765.76747.95764.4004764.400517.22−0.01−0.1820757.1046759.441615.69743.1136744.145046.51−1.88−2.06Avg7.2118.07−7.19−7.10Table 8
Result comparison of ILSSPRP with and without the adjustment operator (T = 28800 s).
InstanceILSSPRP without adjustmentILSSPRPGapBest (%)GapAvg (%)|N|QBest valueAverage valueCPUBest valueAverage valueCPU301017.885019.33785.2414.286215.344614.00−25.19−26.022012.277813.60115.498.68509.192818.10−41.37−47.95601078.479079.55034.8274.780875.396716.96−4.95−5.512067.873669.47465.7463.588664.140719.10−6.74−8.3212010227.0844229.59967.66220.4964221.820524.22−2.99−3.5120202.2880204.22009.13193.2868195.197833.08−4.66−4.6218010477.6772477.867110.48475.4784476.172132.82−0.46−0.3620463.6724465.063012.46454.5760454.992036.97−2.00−2.2124010624.5778626.352413.92623.5814625.904041.91−0.16−0.0720604.0832606.060414.53582.9796583.620756.05−3.62−3.8430010732.4820734.762615.30719.6832720.362153.98−1.78−2.0020723.5718725.361217.74711.2862712.261557.28−1.73−1.84Avg10.2133.71−7.97−8.85However, the GapBest and GapAvg values for all discussed instances are negative, which implies that the adjustment operator is necessary to improve the quality of the solution. As expected, when T increases, the absolute average value of GapBest and GapAvg increases. This also indicates that the ILSSPRP with the adjustment operator is more effective than the ILSSPRP without the adjustment operator to solve the SPRP with long repositioning time. Therefore, it can be concluded that it is beneficial to incorporate adjustment operator into the ILSSPRP to solve the SPRP.


### 5.3. Sensitivity Analysis onα

In this set of experiments, we analyze the sensitivity of the performance of the proposed approach to the values ofα. Note that increasing the value of α means that we are increasing the impact of operational time on the whole weighted objective function. When we set α = 0, it means that the objective function only considers user dissatisfaction.We test the sensitivity ofα with a small-size instance consisting of 30 stations (|N| = 30) and consider the three repositioning time constraints (T = 14400 s, T = 21600 s and T = 28800 s) and two truck capacity constraints (Q = 10 and Q = 20).Different values ofα from 0 to 0.001 were considered, and for each value, we perform 20 runs of the experiment. The results obtained with Lingo 18 and the ILSSPRP were compared in Table 9. CPU (in seconds) is the computing time for Lingo or the average computing time of 20 runs for the ILSSPRP. GapBest and GapAvg are the deviations of the best/average of the objective values from 20 runs obtained by the ILSSPRP from the corresponding optimal solution obtained by Lingo, respectively.Table 9
Sensitivity analysis onα.
InstanceLingoILSSPRPGapBest (%)GapAvg (%)αTQOptimal valueCPUBest valueAverage valueCPU0144001048.500026148.500048.50002.950.000.002043.2000360040.600041.13002.98−6.40−5.03216001029.8000360029.500030.02003.29−1.020.732022.3000360020.500021.20005.05−8.78−5.19288001014.4000360013.000013.97505.79−10.77−3.042012.700036009.100010.04447.20−39.56−26.440.00001144001048.644023348.644048.64404.270.000.002043.5440360041.438641.43864.58−5.08−5.08216001028.8148360028.517828.70408.32−1.04−0.392023.5148360021.213021.73989.42−10.85−8.16288001020.0888360014.286215.344614.00−40.62−30.922013.586236008.68509.192818.10−56.43−47.790.0001144001049.940018149.940049.94002.610.000.002045.1400360041.986042.86512.85−7.51−5.31216001032.3360360028.342030.35279.39−14.09−6.532023.8300360022.936023.300911.71−3.90−2.27288001018.4800360016.626016.96099.78−11.15−8.962016.3140360010.832011.742016.54−50.61−38.940.001144001062.90009662.900062.90003.350.000.002057.5400360054.460055.39644.31−5.66−3.87216001050.6800360047.200048.82407.33−7.37−3.802044.0000360041.220042.470710.74−6.74−3.60288001042.4000360039.580040.573311.91−7.12−4.502036.3000360036.240036.255022.24−0.17−0.12Avg3032.138.28−12.29−8.72As shown in Table9, almost all of the values of GapBest and GapAvg are negative (with only one exception for the results obtained under α = 0, T = 21600 s, Q = 10). The average values of GapBest and GapAvg are −12.29 and −8.72, respectively. This then indicates that the proposed ILSSPRP can always obtain better solutions than Lingo under different values of α, repositioning time and truck capacity. The average computing time of the ILSSPRP is 8.28 s, which is significantly less than that of Lingo (3032.13 s). Therefore, we can conclude that the ILSSPRP works well under different values of α.


### 5.4. Sensitivity Analysis onwi

In this section, a small example is constructed to investigate the sensitivity of the proposed model to the parameterwi. We consider 2 scenarios and investigate the repositioning strategies, respectively:(i)
Scenario 1: we solve the model without taking into account the different values ofwi(ii)
Scenario 2: we consider different values ofwi and solve the model accordinglyWe select the first six stations in the instance of |N| = 30 mentioned in Section 5.1, and set T = 7200 s and Q = 10. The tested values of the parameters including di and wi are shown in Table 10. Note that when di > 0, station i is a pickup station; otherwise, if di < 0, station i is a delivery station. The results are displayed in Table 11, including the value of the objective function, the route of the truck, and the loading/unloading quantities at each station.Table 10
Parameter setting.
Station parameter123456di96−6−68−5wi0.80.60.80.60.20.7Table 11
The value of objective function and the route of the truck in different scenarios.
ScenarioThe value of objective functionRoute and loading/unloading quantities at stations118.06480 ⟶ 2(+5) ⟶ 6(−5) ⟶ 1(+6) ⟶ 4(−6) ⟶ 027.67080 ⟶ 1(+9) ⟶ 3(−6) ⟶ 2(+2) ⟶ 6(−5) ⟶ 0As shown in Table11, the experimental results of scenario 1 and scenario 2 are significantly different. For pickup station 1 and 2, in scenario 2, the value of w1 is relatively large. This means that the level of user dissatisfaction (for one bike) of station 1 is higher than that of station 2. Therefore, in order to reduce the total cost, it is necessary to increase the loading quantity at station 1 to meet the demand of station 1 as much as possible. The loading quantity of station 1 increases from 6 bikes in scenario 1 (where we disregard the different values of wi) to 9 bikes in scenario 2 (d1 = 9, the maximum loading quantity of station 1 is 9 bikes); similarly, the loading quantity of station 2 decreases from 5 bikes in scenario 1 to 2 bikes in scenario 2.Regarding the delivery stations 3, 4, and 6, in scenario 2, becausew3>w6>w4, the priorities of meeting the demand are first given to station 3, followed by station 6, and last to station 4. From the results, we can see that the unloading quantity of station 3 increases from 0 bike in scenario 1 to 6 bikes in scenario 2 (d3 = −6, the maximum unloading quantity of station 3 is 6 bikes), while the unloading quantity of station 4 decreases from 6 bikes in scenario 1 to 0 bike in scenario 2.From the above results, it can be seen that different values ofwi can directly impact the loading/unloading quantities of repositioning truck at each station, as well as the route of the truck, i.e., the repositioning strategy of the BSS operators. Therefore, it can be concluded that it is necessary to take distinguishing the user dissatisfaction generated by different stations into account when modeling the bike repositioning problem, which is the key factor to the optimization of the bike repositioning.


## 5.1. Comparison between Lingo and ILSSPRP

In terms of the quality of the solution and the calculation efficiency, Tables3–5 show the comparison results obtained by Lingo and ILSSPRP under different repositioning times (T = 14400 s, T = 21600 s, and T = 28800 s), truck capacity (Q = 10 and Q = 20), and number of stations (|N|). GapBest and GapAvg (%) indicate the performance of the ILSSPRP relative to that of Lingo based on the best and average objective values, respectively. The CPU (in seconds) indicates the average computing time of the ILSSPRP or the computing time for Lingo, respectively.Table 3
Result comparison of Lingo and ILSSPRP (T = 14400 s).
InstanceLingoILSSPRPGapBest (%)GapAvg (%)|N|QOptimal valueCPUBest valueAverage valueCPU301048.644023348.644048.64404.270.000.002043.5440360041.438641.43864.58−5.08−5.086010117.14203600117.6428117.74245.330.430.5120113.44403600111.6422112.98244.86−1.61−0.4112010273.74403600271.2440271.91016.65−0.92−0.6720296.81803600256.9386257.01926.74−15.52−15.4818010554.14403600538.4374538.72098.02−2.92−2.8620542.04403600528.2434528.84308.11−2.61−2.5024010690.0386690.03868.4720671.9338671.93388.9130010796.3374796.90978.9220789.1392792.123610.63Avg3179.137.12−3.53−3.31Table 4
Result comparison of Lingo and ILSSPRP (T = 21600 s).
InstanceLingoILSSPRPGapBest (%)GapAvg (%)|N|QOptimal valueCPUBest valueAverage valueCPU301028.8148360028.517828.70408.32−1.04−0.392023.5148360021.213021.73989.42−10.85−8.16601094.2160360093.213093.72758.12−1.08−0.522084.4150360081.722683.947412.92−3.29−0.5612010249.31603600243.5154244.153115.46−2.38−2.1120229.41543600222.4112224.314214.45−3.15−2.2718010513.21603600506.4052506.405211.77−1.34−1.3420594.46103600493.0178494.476122.22−20.58−20.2224010656.0130656.462217.6920629.4130630.375732.6930010764.4004764.400517.2220743.1136744.145046.51Avg360018.07−5.46−4.45Table 5
Result comparison of Lingo and ILSSPRP (T = 28800 s).
InstanceLingoILSSPRPGapBest (%)GapAvg (%)|N|QOptimal valueCPUBest valueAverage valueCPU301020.0888360014.286215.344614.00−40.62−30.922013.586236008.68509.192818.10−56.43−47.79601079.5874360074.780875.396716.96−6.43−5.562065.7868360063.588664.140719.10−3.46−2.5712010223.88803600220.4964221.820524.22−1.54−0.9320243.31663600193.2868195.197833.08−25.88−24.6518010593.44323600475.4784476.172132.82−24.81−24.6320588.86603600454.5760454.992036.97−29.54−29.4224010623.5814625.904041.9120582.9796583.620756.0530010719.6832720.362153.9820711.2862712.261557.28Avg360033.71−23.59−20.81From Tables3–5, we can see that Lingo is not able to obtain any feasible solutions for the large-size problems (where |N| = 240 and |N| = 300). It also generally fails to obtain any global optimal solutions within 3600 seconds even for small-size problems. In contrast, for the various size of problems from small to large, the average computing time of the ILSSPRP varies from 7.12 s to 33.71 s.As shown in the Tables3–5, almost all of the values of GapBest and GapAvg are negative (with only one exception for the results obtained under |N| = 60, Q = 10, T = 14400 s). The absolute values of the average GapBest and GapAvg values obtained by the ILSSPRP for the scenarios with long repositioning time are larger than those with short repositioning time. GapBest and GapAvg are, respectively, −3.53 and −3.31 for T = 14400 s, and −5.46 and −4.45 for T = 21600 s, but −23.59 and −20.81 for T = 28800 s. This implies that the ILSSPRP performs better than Lingo when solving the SPRP with a longer repositioning time.Therefore, we can draw the conclusion that the proposed ILSSPRP can obtain better feasible solutions much faster than Lingo.


## 5.2. The Effectiveness of Introducing Adjustment Operator into ILSSPRP

In this study, the adjustment operator (presented in Section4.5), which considers the characteristics of the SPRP, is introduced to improve the solution quality obtained by the proposed ILSSPRP. In this section, we study through experiments the effectiveness of the proposed adjustment operator by comparing the performance of the ILSSPRP with and without adjustment operator using the various instances described in Section 5.1.Tables6–8 present the results obtained by the proposed ILSSPRP algorithm with and without the adjustment operator, under different repositioning times, respectively. CPU denotes the average computing time in seconds. GapBest and GapAvg here are the deviations of the best/average of the objective values from 20 runs obtained by the ILSSPRP with the adjustment operator from the corresponding optimal solution obtained by ILSSPRP without the adjustment operator, respectively. From Tables 6–8, we can see ILSSPRP without the adjustment operator resulting in a shorter computation time than the ILSSPRP with the adjustment operator.Table 6
Result comparison of ILSSPRP with and without the adjustment operator (T = 14400 s).
InstanceILSSPRP without adjustmentILSSPRPGapBest (%)GapAvg (%)|N|QBest valueAverage valueCPUBest valueAverage valueCPU301051.644051.94402.0648.644048.64404.27−6.17−6.782048.338648.90392.2341.438641.43864.58−16.65−18.026010118.1410118.14101.49117.6428117.74245.33−0.42−0.3420125.8314125.83152.11111.6422112.98244.86−12.71−11.3712010277.6434278.06502.20271.2440271.91016.65−2.36−2.2620262.0434263.84932.42256.9386257.01926.74−1.99−2.6618010551.5374552.83163.86538.4374538.72098.02−2.43−2.6220536.1428537.96504.44528.2434528.84308.11−1.50−1.7224010691.0392691.03923.92690.0386690.03868.47−0.15−0.1520675.6338675.63384.77671.9338671.93388.91−0.55−0.5530010797.3254797.32544.09796.3374796.90978.92−0.12−0.0520794.0284794.02844.48789.1392792.123610.63−0.62−0.24Avg3.177.12−3.81−3.90Table 7
Result comparison of ILSSPRP with and without the adjustment operator (T = 21600 s).
InstanceILSSPRP without adjustmentILSSPRPGapBest (%)GapAvg (%)|N|QBest valueAverage valueCPUBest valueAverage valueCPU301032.911233.46013.3128.517828.70408.32−15.41−16.572027.304628.12554.3821.213021.73989.42−28.72−29.376010103.0052103.87583.5193.213093.72758.12−10.51−10.832095.403495.40444.5981.722683.947412.92−16.74−13.6512010252.6992252.99564.13243.5154244.153115.46−3.77−3.6220234.0136235.92095.58222.4112224.314214.45−5.22−5.1718010509.4076509.40786.83506.4052506.405211.77−0.59−0.5920495.8064496.653710.66493.0178494.476122.22−0.57−0.4424010659.7034659.70347.88656.0130656.462217.69−0.56−0.4920643.8136644.399612.00629.4130630.375732.69−2.29−2.2230010764.5106765.76747.95764.4004764.400517.22−0.01−0.1820757.1046759.441615.69743.1136744.145046.51−1.88−2.06Avg7.2118.07−7.19−7.10Table 8
Result comparison of ILSSPRP with and without the adjustment operator (T = 28800 s).
InstanceILSSPRP without adjustmentILSSPRPGapBest (%)GapAvg (%)|N|QBest valueAverage valueCPUBest valueAverage valueCPU301017.885019.33785.2414.286215.344614.00−25.19−26.022012.277813.60115.498.68509.192818.10−41.37−47.95601078.479079.55034.8274.780875.396716.96−4.95−5.512067.873669.47465.7463.588664.140719.10−6.74−8.3212010227.0844229.59967.66220.4964221.820524.22−2.99−3.5120202.2880204.22009.13193.2868195.197833.08−4.66−4.6218010477.6772477.867110.48475.4784476.172132.82−0.46−0.3620463.6724465.063012.46454.5760454.992036.97−2.00−2.2124010624.5778626.352413.92623.5814625.904041.91−0.16−0.0720604.0832606.060414.53582.9796583.620756.05−3.62−3.8430010732.4820734.762615.30719.6832720.362153.98−1.78−2.0020723.5718725.361217.74711.2862712.261557.28−1.73−1.84Avg10.2133.71−7.97−8.85However, the GapBest and GapAvg values for all discussed instances are negative, which implies that the adjustment operator is necessary to improve the quality of the solution. As expected, when T increases, the absolute average value of GapBest and GapAvg increases. This also indicates that the ILSSPRP with the adjustment operator is more effective than the ILSSPRP without the adjustment operator to solve the SPRP with long repositioning time. Therefore, it can be concluded that it is beneficial to incorporate adjustment operator into the ILSSPRP to solve the SPRP.


## 5.3. Sensitivity Analysis onα

In this set of experiments, we analyze the sensitivity of the performance of the proposed approach to the values ofα. Note that increasing the value of α means that we are increasing the impact of operational time on the whole weighted objective function. When we set α = 0, it means that the objective function only considers user dissatisfaction.We test the sensitivity ofα with a small-size instance consisting of 30 stations (|N| = 30) and consider the three repositioning time constraints (T = 14400 s, T = 21600 s and T = 28800 s) and two truck capacity constraints (Q = 10 and Q = 20).Different values ofα from 0 to 0.001 were considered, and for each value, we perform 20 runs of the experiment. The results obtained with Lingo 18 and the ILSSPRP were compared in Table 9. CPU (in seconds) is the computing time for Lingo or the average computing time of 20 runs for the ILSSPRP. GapBest and GapAvg are the deviations of the best/average of the objective values from 20 runs obtained by the ILSSPRP from the corresponding optimal solution obtained by Lingo, respectively.Table 9
Sensitivity analysis onα.
InstanceLingoILSSPRPGapBest (%)GapAvg (%)αTQOptimal valueCPUBest valueAverage valueCPU0144001048.500026148.500048.50002.950.000.002043.2000360040.600041.13002.98−6.40−5.03216001029.8000360029.500030.02003.29−1.020.732022.3000360020.500021.20005.05−8.78−5.19288001014.4000360013.000013.97505.79−10.77−3.042012.700036009.100010.04447.20−39.56−26.440.00001144001048.644023348.644048.64404.270.000.002043.5440360041.438641.43864.58−5.08−5.08216001028.8148360028.517828.70408.32−1.04−0.392023.5148360021.213021.73989.42−10.85−8.16288001020.0888360014.286215.344614.00−40.62−30.922013.586236008.68509.192818.10−56.43−47.790.0001144001049.940018149.940049.94002.610.000.002045.1400360041.986042.86512.85−7.51−5.31216001032.3360360028.342030.35279.39−14.09−6.532023.8300360022.936023.300911.71−3.90−2.27288001018.4800360016.626016.96099.78−11.15−8.962016.3140360010.832011.742016.54−50.61−38.940.001144001062.90009662.900062.90003.350.000.002057.5400360054.460055.39644.31−5.66−3.87216001050.6800360047.200048.82407.33−7.37−3.802044.0000360041.220042.470710.74−6.74−3.60288001042.4000360039.580040.573311.91−7.12−4.502036.3000360036.240036.255022.24−0.17−0.12Avg3032.138.28−12.29−8.72As shown in Table9, almost all of the values of GapBest and GapAvg are negative (with only one exception for the results obtained under α = 0, T = 21600 s, Q = 10). The average values of GapBest and GapAvg are −12.29 and −8.72, respectively. This then indicates that the proposed ILSSPRP can always obtain better solutions than Lingo under different values of α, repositioning time and truck capacity. The average computing time of the ILSSPRP is 8.28 s, which is significantly less than that of Lingo (3032.13 s). Therefore, we can conclude that the ILSSPRP works well under different values of α.


## 5.4. Sensitivity Analysis onwi

In this section, a small example is constructed to investigate the sensitivity of the proposed model to the parameterwi. We consider 2 scenarios and investigate the repositioning strategies, respectively:(i)
Scenario 1: we solve the model without taking into account the different values ofwi(ii)
Scenario 2: we consider different values ofwi and solve the model accordinglyWe select the first six stations in the instance of |N| = 30 mentioned in Section 5.1, and set T = 7200 s and Q = 10. The tested values of the parameters including di and wi are shown in Table 10. Note that when di > 0, station i is a pickup station; otherwise, if di < 0, station i is a delivery station. The results are displayed in Table 11, including the value of the objective function, the route of the truck, and the loading/unloading quantities at each station.Table 10
Parameter setting.
Station parameter123456di96−6−68−5wi0.80.60.80.60.20.7Table 11
The value of objective function and the route of the truck in different scenarios.
ScenarioThe value of objective functionRoute and loading/unloading quantities at stations118.06480 ⟶ 2(+5) ⟶ 6(−5) ⟶ 1(+6) ⟶ 4(−6) ⟶ 027.67080 ⟶ 1(+9) ⟶ 3(−6) ⟶ 2(+2) ⟶ 6(−5) ⟶ 0As shown in Table11, the experimental results of scenario 1 and scenario 2 are significantly different. For pickup station 1 and 2, in scenario 2, the value of w1 is relatively large. This means that the level of user dissatisfaction (for one bike) of station 1 is higher than that of station 2. Therefore, in order to reduce the total cost, it is necessary to increase the loading quantity at station 1 to meet the demand of station 1 as much as possible. The loading quantity of station 1 increases from 6 bikes in scenario 1 (where we disregard the different values of wi) to 9 bikes in scenario 2 (d1 = 9, the maximum loading quantity of station 1 is 9 bikes); similarly, the loading quantity of station 2 decreases from 5 bikes in scenario 1 to 2 bikes in scenario 2.Regarding the delivery stations 3, 4, and 6, in scenario 2, becausew3>w6>w4, the priorities of meeting the demand are first given to station 3, followed by station 6, and last to station 4. From the results, we can see that the unloading quantity of station 3 increases from 0 bike in scenario 1 to 6 bikes in scenario 2 (d3 = −6, the maximum unloading quantity of station 3 is 6 bikes), while the unloading quantity of station 4 decreases from 6 bikes in scenario 1 to 0 bike in scenario 2.From the above results, it can be seen that different values ofwi can directly impact the loading/unloading quantities of repositioning truck at each station, as well as the route of the truck, i.e., the repositioning strategy of the BSS operators. Therefore, it can be concluded that it is necessary to take distinguishing the user dissatisfaction generated by different stations into account when modeling the bike repositioning problem, which is the key factor to the optimization of the bike repositioning.


## 6. Conclusions

In this paper, we study the static partial repositioning problem where the level of user dissatisfaction differs in different stations. We propose an iterated local search algorithm ILSSPRP to find solutions to the problem as well as an adjustment operator to further improve the obtained solutions. We conduct a comprehensive set of experiments to study the performance of the proposed approach in terms of solution quality and computational time, compared with Lingo. The experimental results show that the proposed ILSSPRP almost always outperforms Lingo from small- to large-size problems, as well as under different relative importance set to the two objectives (i.e., operational time and deviation from the target number of bikes). We also use a small example to illustrate that the unit penalty cost wi has an effect on the repositioning route and the loading/unloading quantities at each station.There are potential directions for future research. For instance, unbalance demands in the bike-sharing system will not only increase user dissatisfaction, but also aggravate urban traffic congestion. Therefore, in addition to considering economic benefits, we will also need to consider environmental benefits in the bike repositioning problem. How to introduce environmental benefits into the objective function would be an interesting topic for future research. Besides, we are also interested in developing efficient and effective algorithms to verify feasibility of a given solution.



---

*Source: 3040567-2020-11-17.xml*