The Biobjective Bike-Sharing Rebalancing Problem with Balance Intervals: A Multistart Multiobjective Particle Swarm Optimization Algorithm

(e bike-sharing system (BSS), as a sustainable way to deal with the “last mile” problem of mass transit systems, is increasingly popular in recent years. Despite its success, the BSS tends to suffer from the mismatch of bike supply and user demand. BSS operators have to transfer bikes from surplus stations to deficit stations to redistribute them among stations by means of trucks. In this paper, we deal with the bike-sharing rebalancing problem with balance intervals (BRP-BIs), which is a variant of the static bike-sharing rebalancing problem. In this problem, the equilibrium of station is characterized by a balance interval instead of a balance point in the literature. We formulate the BRP-BI as a biobjective mixed-integer programming model with the aim of determining both theminimum cost route for a single capacitated vehicle and the maximum average rebalance utility, an index for the balanced degree of station.(en, a multistart multiobjective particle swarm optimization (MS-MOPSO) algorithm is proposed to solve the model such that the Pareto optimal solutions can be derived. (e proposed algorithm is extended with crossover operator and variable neighbourhood search to enhance its exploratory capability. Compared with Hybrid NSGA-II andMOPSO, the computational experimental results demonstrate that our MS-MOPSO can obtain Pareto optimal solutions with higher quality.


Introduction
Nowadays, with the increasing concerns regarding global warming, countries around the world have begun to save energy and reduce carbon emissions. In order to achieve the low-carbon and sustainable transport, bike-sharing systems (BSSs) are increasingly popular, as they can reduce travel costs, road congestion, and greenhouse gas emissions [1]. e first BSS was introduced in Amsterdam in 1965, and then its popularity increased as a new short-distance transport mode in the last two decades across the globe [2,3].
Despite its success, the BSS tends to suffer from the mismatch of bike supply and user demand. is happens when a user arrives at a station that has no bike available or that has no parking slot available. As can be seen from Figure 1, there are evident morning and evening peaks in the bike-sharing usage on workdays. During the morning rush hour, users pick up bikes in residential areas, make a oneway point to point transportation, and return them in working areas. is morning peak leads to excessive bikes in working areas as well as insufficient bikes in residential areas, depicted in Figure 2, and vice versa during the evening rush hour. is mismatch results in poor service level and low customer satisfaction, for the BSS has little chance to complete self-balance in a short time.
To improve the service level, it is necessary for the BSS operator to realize the rebalance between bike supply and user demand by transporting bikes from surplus stations to deficient stations with the help of trucks, referred to as the bike-sharing rebalancing problem (BRP) [4]. BRPs can be divided into static problems and dynamic ones according to whether user activity is negligible or not.
In this paper, we focus on the static bike-sharing rebalancing problem with balance intervals (BRP-BIs), where the equilibrium of station is characterized by balance interval instead of balance point in the literatures. And then, we formulate this problem as a biobjective mixed-integer programming model with objectives of minimizing the total travel distance as well as maximizing the average rebalance utility, an index for the balanced degree at each station. e definition and description of balance interval and rebalance utility are described in Section 3.
BRP-BI, as a variant of one-commodity pickup-anddelivery traveling salesman problem (1-PDTSP) [5], is a complex problem that includes decisions regarding the route that the vehicle should follow and the number of bikes that should be loaded or unloaded at each station. It is more complicated than traveling salesman problem, a well-known NP-hard problem, since the numbers of bikes loaded or unloaded are decision variables. Besides, BRP-BI is based on the balance interval, and it has two conflicting objectives. It is therefore more complex and more difficult to be solved exactly. As the problem size increases, the calculation time will increase exponentially. us, the real-life BRP-BI instances cannot be solved exactly within short computation time.
To handle this biobjective NP-hard problem in a runtime acceptable to BSS operators, we develop the multistart multiobjective particle swarm optimization (MS-MOPSO) algorithm. MOPSO has been found to be successful to solve a wide variety of multiobjective optimization problems, but until now, it has not been extended to deal with biobjective BRP-BI. MS-MOPSO is thus proposed to obtain the Pareto optimal solutions for this biobjective problem. To avoid falling into the local optimum and to enhance the exploratory capability, crossover operator and variable neighborhood search are incorporated to enable MS-MOPSO to achieve balanced exploration and exploitation capabilities. e main contributions of this paper lie on the following: (i) To introduce a new and practical concept of balance interval, which is different from balance point in the literatures, and then, BRP-BI is presented (ii) To construct a biobjective mixed-integer programming model to formulate this problem (iii) To develop a MS-MOPSO enhanced by the crossover operator and variable neighborhood search to deal with this biobjective model e remainder of this paper is organized as follows: Section 2 summarizes the literature review on the related problem. Section 3 formally proposes the concept of balance interval and rebalance utility and then the biobjective mixedinteger programming model for BRP-BI. Section 4 develops the MS-MOPSO to obtain the nondominated optimal solutions for BRP-BI. Section 5 contains the computational experiments and sensitivity analysis. Finally, Section 6 draws some conclusions about this work.

Literature Review
During the past decade, BRPs have attracted the considerable attention in the literatures. Generally speaking, they can be roughly classified into two categories: static and dynamic, depending on whether the user activity of a BSS is negligible or not. If the user activity is low and has no effect on the rebalancing operations (such as during the midnight), it is static BRP (SBRP). Otherwise, it is dynamic BRP (DBRP) [6]. In the following, we will discuss the literatures related to the SBRP and DBRP, respectively.

Static Bike-Sharing Rebalancing
Problem. Many approaches to solve the SBRP revolve around optimization models that maximize profit or minimize travel costs. Dell'Amico et al. [7] provide four mixed-integer programming models with the objective of minimizing total costs, where a fleet of capacitated vehicles is employed to redistribute the bikes. e same objective is presented by Zhang et al. [8], who propose a hybrid discrete PSO algorithm to solve it. Erdogan et al. [5] and Cruz et al. [9] solve the same model with a single vehicle, where multiple visits are allowed at each station. Duan et al. [10] focus on a vehicle routing problem in the rebalancing operations with the objective of minimizing the total driving distance, and a greedy algorithm is proposed. Shi et al. [11] propose an improved PSO algorithm to minimize the total costs including the fixed costs and the variable transport costs. In addition, some researchers have taken time into consideration in the SBRP models. Casazza et al. [12] incorporate time into the constraints such that each route does not exceed time limitation.
In the study of Bulhões et al. [3], service time is limited and the objective is to minimize the travel time. All the problems above need to fully meet the rebalance demands, but other researchers are trying to relax these constraints. Papazek et al. [13] set the primary goal to minimize the absolute deviation between target and final fill levels for all stations. Gaspero et al. [14] solve the problem with the aim of minimizing the weighted sum of the total travel time and the total absolute deviation from the target number of bikes, while Raviv and Kolka [15] take it as a penalty cost. Szeto et al. [16] and Espegren et al. [17] also consider the same factors, but the objective is to minimize the weighted sum of unmet customer demand and operational time. Meanwhile, Ho and Szeto [18] define the objective as the sum of the penalty cost incurred at each station and the weighted total travel time. Raviv et al. [19] define the goal of the problem as minimizing the total cost of the system consisting of the weighted sum of the penalties incurred at all stations and the total operating costs. Ho and Szeto [20] present an iterated tabu search heuristic process to solve the SBRP with the goal of minimizing the total penalties incurred at all stations. Li et al. [21] consider solving a multitype bike rebalancing problem, which is formulated as a mixed-integer programming problem, aiming at minimizing the total costs consisting of the travel cost, the penalties due to unmet demands, and the penalties associated with the substitution and occupancy strategies. Broken bikes are considered by Wang and Szeto [22] with the objective of minimizing the total CO 2 emissions of all vehicles. Usama et al. [23] also consider replacing faulty bikes in the system with the following objectives: user dissatisfaction to unmet demand, vehicle routing costs including handling, loading, unloading, and traveling cost, and the waiting time of repositioning vehicles.

Dynamic
Bike-Sharing Rebalancing Problem. Due to the participations of users, there is a great deal of uncertainties in the BSS during the daytime and the SBRP models have failed to meet the actual dynamic rebalance demands. erefore, some scholars have turned their attention to the DBRP under uncertainty circumstances.
Some researchers deal with uncertainty through predictions. Alvarez-Valdes et al. [2] estimate the unsatisfied demand and solve the problem with the objective of minimizing the global expected dissatisfaction for each possible number of bicycles at the beginning of the period. And then they use these estimates to guide redistribution algorithms. Zhang et al. [24] propose a dynamic bicycle repositioning methodology that considers inventory-level forecasting, user arrivals forecasting, bicycle repositioning, and vehicle routing in a unified manner. A multicommodity time-space network flow model is presented with the objective of minimizing the total vehicle travel costs plus the expected user dissatisfaction in the system. Caggiani et al. [25] propose a new comprehensive methodology that starts from the prediction of the number and position of bikes and ends with a relocation decision support system, aiming at achieving a high degree of user satisfaction and keeping the vehicle repositioning costs as low as possible.
In addition, other researchers embody uncertainty through stochastic demand based on scenarios. In the study of Yan et al. [26], the objective is to maximize the expected profit of the system, that is, the expected total bike rental income on the demand arcs minus the expected total bike holding cost on the holding arcs over all scenarios. Dell'Amico et al. [27] develop the stochastic programming model with the objective of minimizing the travel cost and the penalty costs for unfulfilled demands. Maggioni et al. [28] propose the two-stage and multistage stochastic optimization models to determine the optimal number of bikes to assign to each station at the beginning of the service. e objective is to minimize the sum of the procurement cost for the assigned bikes, the expected stock-out cost for shortage, the expected time-waste cost for overflow, and the expected transshipment cost for repositioning bikes at the end of the service. Yuan et al. [29] present a mixed-integer programming model with the objective of minimizing the daily costs including the capital cost and the expected operating cost.
Furthermore, there are still some studies dealing directly with dynamic factors. Legros [30] focuses on dynamic repositioning strategy in a BSS with the objective of minimizing the rate of arrival of unsatisfied users and solves it by a Markov decision process approach. You [31] develops a constrained mathematical model to deal with a multivehicle BRP, aiming to develop dynamic decisions to minimize the Complexity 3 sum of the transportation costs and unmet costs under service-level constraints over a planning horizon. Cao et al. [32] propose a station location problem under a zoning scheme in consideration of multiperiod demand, and the objective of this problem is to minimize the fixed construction costs and variable operational costs. Warrington and Ruchti [33] present a novel nested-flow formulation of the dynamic problem to minimize a generalized cost function, which includes penalty costs for unserved customer demands, costs for the movements of vehicles, and costs of loading and unloading. Regue and Recker [34] propose a framework to solve the DBRP based on four models: a demand forecasting model, a station inventory model, a redistribution needs model, and a vehicle-routing model. eir approach is proactive instead of reactive, as bike repositioning occurs before inefficiencies are observed.

Summary.
Most studies directly give out the unbalanced demands in their BRP models without involving the concept of equilibrium of station. According to the equilibrium of station, model objective, and solution approach, Table 1 summarizes the literature considering the equilibrium of station.
In terms of equilibrium of station, the majority of these studies regard it as a balance point, but only two studies [5,35] begin to pay attention to the concept of interval. Obviously, the balance point is a special case of the balance interval. e advantage of the balance interval is that it can decrease the number of rebalancing operations, thereby reducing the complexity of the BRP model. In the study of Erdogan et al. [5], demand interval is proposed to refer to an inventory of bicycles lying between given lower and upper bounds. In their BRP model, the bounding constraints must be satisfied. Kadri et al. [35] propose a known threshold levels indicating the minimum and maximum number of bicycles required for each station, and the rebalance demand is used only to meet the boundary of the interval. Different from their research, we define the equilibrium of station as the balance interval, and the inventory of bikes at a station can be beyond this interval. If so, the rebalance utility of the station will decrease. In addition, one objective of our biobjective BRP-BI is to maximize the average rebalance utility of all stations. And another objective is to minimize the vehicle travel distance.
Regarding the model objective and the solution approach, the monoobjective BRPs as well as their corresponding exact or heuristic solution approaches have been studied in many recent works. But the multiobjective counterpart has not widely been considered yet, and only two papers start to go beyond the scope of monoobjective. Mohammadi and Shirouyehzad [36] devise a multiobjective multiperiod BRP to minimize the slacks and surpluses in stations as well as the rebalancing costs. Szeto and Shui [37] investigate a BRP with bilevel objective. e first priority objective is to minimize the total demand dissatisfaction, and the second priority objective is to minimize the service time. However, no multiobjective method has been proposed to solve these problems. For real-life biobjective BRP-BI, it is more realistic and more useful to calculate the Pareto optimal solutions instead of just one solution. erefore, we will design a multiobjective solution approach, MS-MOPSO, to solve our biobjective BRP-BI. is will expand the existing research direction of BRP.

Mathematical Formulation
BRP-BI studied in this paper is defined on a complete graph is the set of arcs, and d ij is the distance between vertex i and vertex j.

Balance Interval.
In most literatures, such as the study of Dell'Amico et al. [27], the equilibrium of a BSS station is defined to a prespecified balance point, namely, the target number of bikes in station. Each station must be brought from its current situation to its equilibrium situation to increase the service level of the BSS as well as users' satisfaction. Suppose the equilibrium of station j is b j and the current inventory of bikes is l 1 j . If l 1 j � b j , station j is defined as balanced and its rebalance demand q j � 0. Otherwise, station j is unbalanced, and its rebalance demand To avoid the defects caused by the balance point, the equilibrium of station is redefined in this paper as balance interval, that is, an interval of inventory of bikes, depicted in Figure 4. By introducing threshold parameter θ, let (1 − θ)b j and (1 + θ)b j be the lower and upper bounds of balance interval of station j. Unlike the balance point, if , station j is defined as balanced, and its rebalance demand q j � 0. In such a case, station j can be removed from set N, and then N ′ � N/ j and some bikes need to be picked up from station j. e balance interval can reduce the number of rebalancing operations, thereby reducing the complexity of BRP-BI. Furthermore, the balance point is a special case of our balance interval when θ � 0. Obviously, the different value of θ will result in different sizes of balance interval and then different complexity of BRP-BI. e discussion on the value of θ will be carried out later in Section 5.

Rebalance Utility.
Each station must carry both enough bikes to meet users' rental demands and enough vacant slots to satisfy users' return demands. Since the station has a fixed capacity, a larger quantity of bikes implies a smaller quantity of vacant slots. erefore, the measurement of station service capability includes the bike rental capability as well as the bike return capability, depicted in Figure 5. According to the number of bikes at station j, there are three intervals: In the left interval, when the station is empty, bicycles cannot be rented and all slots are available to return, so the rental capability is 0 and the return capability is 1. With the increase in the number of bikes, the rental capability increases, but the return capability is still 1 since there are still enough vacant slots to meet users' return requests. When the number of bikes increases to (1 − θ)b j , suppose the bikes available are enough for rental, the rental capability becomes 1. erefore, the station service capability in the left interval depends on the rental capability.
In the middle interval, namely, the balance interval, the station has both enough bikes and enough vacant slots, so the station service capability remains at 1.
In the right interval, as the number of bikes increases, the rental capability remains at 1, but the return capability decreases because there are fewer and fewer vacant slots. When the station is full, on the contrary to empty, the rental capability is still 1 and the return capability is 0. Similarly, the station service capability in the right interval depends on the return capability.
In summary, the station service capability depends on the rental capability in the left interval [0, , and depends on the return capability in the right interval In our study, we propose the rebalance utility function to evaluate the service capability of a station after rebalancing operations, whatever the rebalance demand has fully met or not.
Suppose l 2 j is the inventory of bikes hold by station j after the rebalancing operations, the rebalance utility function μ(l 2 j ) is defined as follows:

References
Equilibrium of station Model objective Solution approach Erdogan et al. [5] Interval Single Benders decomposition, branch-and-cut Dell'amico et al. [7] Point Single Branch-and-cut Cruz et al. [9] Point Single Iterated local search Duan et al. [10] Point Single Greedy algorithm Ho and Szeto [20] Point Single Iterated tabu search heuristic Dell'Amico et al. [27] Point Single L-shaped, branch-and-cut You [31] Point Single Two-phase heuristic Kadri et al. [35] Interval Single Branch-and-bound Mohammadi et al. [36] Point Multiple Lingo Szeto and Shui [37] Point Bilevel Enhanced artificial bee colony algorithm is study Interval For the purpose of promoting the station service capability and improving user satisfaction, maximizing the average rebalance utility is chosen as one objective of our biobjective BRP-BI.

Problem Formulations.
In BRP-BI, a vehicle with capacity Q is used for the rebalancing operations. It leaves the depot with q 1 inventory of bikes, performs a tour visiting each station exactly once to pick up s + j or deliver s − j inventory of bikes at station j ∈ N ′ (n ′ is the number of stations in set N ′ ), and returns to the depot with q 2 inventory of bikes. If the vehicle travels from vertex i to vertex j, the binary decision variable x ij � 1; otherwise, x ij � 0. Another decision variable f ij represents the quantity of bikes flowing on the arc (i, j) ∈ A. One objective is to maximize the average rebalance utility, and the other one is to minimize the total travel distance.
Referring to the model of Dell'Amico et al. [27], the biobjective mixed-integer programming model for BRP-BI can be written as follows: i∈S i∈S Objective function (2) minimizes the total travel distance, and objective function (3) maximizes the average rebalance utility. Obviously, these two objectives cannot be optimized at the same time. Equation (4) imposes that each station has to be visited exactly once. Equation (5) ensures flow conservation, and constraint (6) eliminates subtours. Constraint (7) guarantees that the vehicle capacity is not

Return capability
Rental capability Station service capability (rebalance utility) 1 6 Complexity exceeded throughout its tour. Constraints (8) and (9) specify that the difference between the flows entering and leaving a station is equal to its rebalance demand. Constraint (10) enforces that a station has either delivery rebalance demand or pickup rebalance demand and cannot have both. Constraint (11) guarantees that the bikes picked up at a station cannot exceed the remaining capacity of the vehicle. Constraint (12) enforces that the bikes delivered to a station cannot exceed the inventory hold by the vehicle. Furthermore, constraint (13) restricts the decision variables that are binary.

Solution Method
In recent years, many efficient metaheuristics have been proposed to solve various multiobjective optimization problems, such as EM-BBO [38], NSGA-II [39], MOEA/D [40], and MOPSO [41]. Compared with single-objective methods, multiobjective algorithms can generate a set of nondominated optimal solutions, from which suitable schemes can be selected for different actual scenarios. Due to the advantages of simplicity, low-cost implementation, and fast convergence, MOPSO has become a well-known algorithm for optimizing practical multiobjective problems. To address our biobjective BRP-BI, MS-MOPSO, an improved version of the MOPSO proposed by Coello et al. [42], is presented in this section. In order to achieve balanced exploration and exploitation capabilities, we have introduced variable neighborhood search [43] and multistart strategy [44] to enhance our method. e former can avoid falling into local optima and expand exploration capabilities. e latter can help find higher quality Pareto optimal solutions. e procedure and key components of our algorithm are described in detail below.

Procedure of the Algorithm.
e procedure of our MS-MOPSO is demonstrated in Figure 6. Unlike single-objective PSO, there are two different criteria to determine the personal optimum pbest and the global optimum gbest. We first select the unbalanced stations from the BSS and then use the fitness evaluation function to obtain the nondominated particles. All the nondominated particles are gathered into an external repository (REP), from which every particle chooses its gbest. For pbest, domination-based and probabilistic rules are utilized. To avoid being trapped in local optima, we have introduced the crossover operator of genetic algorithm to update the population and the variable neighborhood search to improve the exploratory capability of our algorithm. Finally, when the number of runs for multistart strategy (ms) reaches the specified number, we will get the nondominated optimal solutions for biobjective BRP-BI.

Initial Solutions Generation.
e stations with inventory of bikes in their balance intervals can be regarded as in equilibrium and do not need to be rebalanced. For example, there are 15 stations, numbered 1-15, whose current inventories of bikes are illustrated in Table 2.
Suppose the target inventory of stations is 50 and θ � 0.2, the balance interval is [40,60]. It can be seen from  Complexity can be removed from the BSS and will not be visited in the rebalancing operations. In brief, the vehicle only needs to visit the unbalanced stations, and this will be helpful to reduce the scale of the problem as well as the total travel distance.
A feasible solution of the BRP-BI is made up of two parts: one is the vehicle route and the other is the rebalance demand. In the following, we will discuss how to construct the initial vehicle route and how to calculate the rebalance demand for each unbalanced station.

Vehicle Route Construction.
A vehicle route can be represented as (0, i 1 , i 2 , . . . , i n , 0), where 0 is the depot and i 1 , i 2 , . . . , i n is a sequence of unbalanced stations that need to be rebalanced. By random permutation for i 1 , i 2 , . . . , i n , a typical "vehicle route" is displayed in the first line of Table 3, where depot 0, as the start and end points, is omitted. Another vehicle route can surely be obtained by changing the permutation of stations. In other words, different permutations of stations correspond to different vehicle routes.

Rebalance Demand Calculation.
e notations used in this section are set as follows: inv i1 and inv i2 are the inventories in station i before and after rebalancing operations, respectively; ideal_inv is the target inventory; lb/ub is the lower/upper boundary of the balance interval; q i is the load on the vehicle after visiting station i; and Q is the capacity of vehicle. e pseudocode for calculation of rebalance demand of station i (i � 1, 2, . . ., n) is demonstrated in Algorithm 1. As mentioned above, the balanced stations have already been eliminated from the rebalancing operations, thus inv i1 is either less than lb or more than ub.
If inv i1 < lb, there is a delivery rebalance demand at station i, and the unloading quantities depend on the current load on vehicle. If the load is enough, station i will be replenished to the target inventory ideal_inv, which is commonly in the middle of the balance interval. Otherwise, the load is insufficient to satisfy the delivery demand, and station i will get all the load on vehicle. In such a case, if inv i2 < lb, the rebalance utility of station i will decrease.
On the contrary, if inv i1 > ub, there is a pickup rebalance demand at station i, and the loading quantities depend on the free space of vehicle. If the free space is enough, the extra inventory will be loaded into the vehicle to make the inventory of station i be ideal_inv; otherwise, the vehicle will be filled. Similarly, if inv i2 > ub, its rebalance utility will also decrease.
In Table 3, the "inventories before rebalancing" of the unbalanced stations also comes from Table 2. Given the target inventory ideal_inv � 50, the "rebalance demands" can be obtained by "inventories before rebalancing" minus ideal_inv, where "+" means the pickup rebalance demand and "− " means the delivery rebalance demand.
Suppose the vehicle's capacity is 50 and its initial load is zero, the "vehicle's load after rebalancing" and "inventories after rebalancing" have to be calculated iteratively by using Algorithm 1. In the first two stations, the rebalancing operations perform well, and the vehicle picks up 48 bikes at station 1 and unloads 18 bikes at station 2. When the vehicle arrives at station 5, its load is 30 bikes, but the delivery demand of station 5 is 31 bikes. According to Algorithm 1, station 5 is replenished by only 30 bikes and is not fully satisfied. is can avoid increasing the travel distance because the vehicle must return back to the depot to pick up one more bike if fully meeting the demands is required. In other words, we have trade-off between the travel distance and the rebalance utility. e similar analysis can be applied to the remaining stations.
Moreover, the "quantities satisfied" are the differences between "inventories before rebalancing" and "inventories after rebalancing." By comparing the "rebalance demands" and "quantities satisfied," the grey columns demonstrate that station 5, 10, and 15 have insufficient bikes and they are not fully satisfied. However, the rebalancing operations are satisfying because "inventories after rebalancing" for all stations is in the balance interval [40,60]. is also illustrates the superiority of our balance interval.
In summary, to generate the initial solution population, different vehicle routes can be constructed by random permutation, and the rebalance demands can be calculated by using Algorithm 1.

Fitness Function and Personal Best Particle.
Particles in PSO are sharing information and moving towards both their own personal best particle pbest and the global best particle gbest. But in multiobjective optimization, "nondominated" takes the place of "best." And there is more than one criterion to determine the nondominated particles in our MS-MOPSO. In this section, we will discuss how to choose the pbest and then in the following, gbest will be involved.
In the biobjective BRP-BI, the dominance relationship of the particles is obtained by comparing their fitness functions originated from the objective functions (2) and (3). e fitness functions used in our algorithm are specified as follows: Based on these fitness functions, the Pareto dominance in our algorithm is defined as follows: for two particles z1 Otherwise, there is no dominance relationship between z1 and z2. Particle z * is said to be a nondominated particle, or Pareto particle, if and only if f 1 (z * ) ≤ f 1 (z) and f 2 (z * ) < f 2 (z), or f 1 (z * ) < f 1 (z) and f 2 (z * ) ≤ f 2 (z) for any particle z. If z * is a nondominated particle, (f 1 (z * ), f 2 (z * )) is named as a nondominated point and all these points constitute the Pareto front.
In an initial solution, for each particle, pbest is itself. In the process of solving, if the current particle can dominate its pbest, pbest is updated to this particle; if pbest dominates the current particle, pbest is not changed; otherwise, the current particle and pbest have no dominance relationship, and one of them is randomly selected as pbest.

REP and Global Best Particle.
In our MS-MOPSO, REP is used to keep a historical record of all nondominated solutions found during the search process. At the beginning of the search, the nondominated particles in the initial solution are all added to REP. Along the search process, the current nondominated particles found at each iteration are compared one by one with all the solutions in REP. If one current particle is dominated by an individual of REP, it is discarded. Otherwise, if all individuals of REP cannot dominate this current particle, such particle can be appended to REP. If there are solutions of REP dominated by the new member, those individuals are removed from REP. If REP has reached its maximum capacity, we use the adaptive grid method proposed by Coello et al. [42] to remove extra solutions and keep the solutions in REP well distributed.
For the selection of gbest from REP, we refer to the method of Dabiri et al. [41]. e solution space is divided into multiple equal grids, and the grids with fewer particles have a higher chance of being selected. If there is more than one particle in the grid selected, one of them is taken randomly as gbest.

Particle Crossover Operator.
e traditional position and velocity update formulas cannot directly apply to renew the particles in MS-MOPSO. erefore, the crossover operator of the genetic algorithm is introduced to update the particles. And our particle crossover operator is defined as follows: where X i (t + 1) and X i (t) are the solutions of the i th particle in the current iteration (t) and the next iteration (t + 1), respectively; pbest i (t) is the personal optimal solution of the i th particle, and gbest(t) is the global optimal solution; and symbol ⊗ is on behalf of the crossover operator. An example is taken to illustrate the framework of our crossover operator. As shown in Figure 7(a), X i (t) and its pbest i are initially selected as the parents P 1 and P 2 , and then the vehicle routes are updated by exchanging their segments of chromosome selected randomly and the corresponding rebalance demands are computed using Algorithm 1. After crossover operations, two offspring particles (O 1 and O 2 ) are generated, displayed in Figure 7(b). ey are compared by their fitness function values computed by formulas (14) and (15), and the nondominated one is selected to perform the subsequent operations. If there is no dominance relationship, one of them is randomly selected. In this example, O 2 is selected as a new parent P 1 ′ and gbest is the other parent P 2 ′ . In the same way, P 1 ′ and P 2 ′ are compared to produce their offsprings O 1 ′ and O 2 ′ , demonstrated in Figure 7(c). Similarly, we take the nondominated individual or a random one (if there is no dominance relationship) as the updated particle X i (t + 1).

Particle Mutation Operator.
For the purpose of increasing the diversity of particles, expanding the search range, and jumping out of local optimum, the particle mutation operator enhanced by variable neighborhood Complexity 9 search is proposed to update the particles. ere are three neighborhood structures used in our algorithm: NS1 : one station is selected randomly and removed from its original position, and then it is inserted into another random position, depicted in Figure 8(a). NS2 : two stations are selected randomly, and their positions are swapped, displayed in Figure 8(b). NS3 : two stations are selected randomly, and the sequence between them is arranged in the reverse order, shown in Figure 8(c).
In our mutation operator, the three structures above are performed randomly. After the mutation operations, if the new particle dominates the original, the new particle will replace the original; if the original dominates the new particle, the new particle will be discarded. Otherwise, if there is no dominance relationship between these two particles, then one of them is randomly chosen.

Results and Discussion
We have implemented the MS-MOPSO algorithm using MATLAB R2014b, and all computational experiments have been carried out on a notebook with Intel Core i5-8265U CPU and 8 GB RAM. e test instances for numerical analysis are generated randomly.

Test Instances and Parameters.
e performance of MS-MOPSO is evaluated by 12 test instances (3 types, 4 instances per type), displayed in Table 4. ese three types are small size (n � 20), medium size (n � 60), and large size (n � 100). In these instances, the stations' coordinates are randomly generated in the range of longitude 31.15°∼31.25°and latitude 121.33°∼121.43°. e depot was assumed to be located at the center of this area, and its coordinate is (31.20°, 121.38°). e station capacity is 100, and its current quantity of bikes is randomly drawn from a range [0, 100]. A single vehicle with capacity of 50 is used to serve the unbalanced stations selected in advance. Unless otherwise specified, we set θ � 0.2 in all instances; thus, the balance interval is [40,60]. e numbers of unbalanced stations are also shown in Table 4.
In the following section, the number of runs for multistart strategy is 3 for small-size instances, 4 for mediumsize instances, and 5 for large-size instances. e parameters of MS-MOPSO and MOPSO are set as follows: particle population is 100, REP size is 100, and mutation rate is 1. And HNSGA-II runs using a population size of 100, a crossover rate of 0.7, a mutation rate of 1, and a mutation percentage of 0.4. For the small-size, medium-size, and large-size instances, the maximal iterations are set to 300, 800, and 1800, respectively.
erefore, multistart strategy of our algorithm is conducive to finding the better Pareto optimal solutions.

Comparisons with HNSGA-II and MOPSO.
To validate the efficiency and efficacy of our algorithm, we compare it with HNSGA-II and MOPSO, using two criteria: the quality of solutions and the running time.
First, we compare the three algorithms by the quality of solutions. Figure 10 displays the graphical results of the Pareto optimal solutions in all instances. For the small-size instances, it is noticeable that the difference is not large. e Pareto optimal solutions for instances 1 and 4 are almost the same. But in instances 2 and 3, the quality of solutions obtained by MS-MOPSO and HNSGA-II is still slightly better than that obtained by MOPSO. For the medium-size and large-size instances, the performance of MS-MOPSO is obviously better than that of MOPSO in all instances. In instances 6, 7, and 11, MS-MOPSO and HNSGA-II have similar performance. But in instances 5, 8, 9, 10, and 12, MS-MOPSO performs better. Especially in instances 5 and 10, it achieves complete superiority. ough some solutions produced by MS-MOPSO are slightly worse than HNSGA-II, the performance of our MS-MOPSO is satisfying if we consider that HNSGA-II is normally considered good enough to solve most multiobjective optimization problems. en we compare the three algorithms by running time. We run 10 times for each instance and each algorithm, and the average running time is shown in Figure 11. Although the solutions of MOPSO are not the best, it has complete advantages in running time in all instances. Obviously, MS-MOPSO takes longer than MOPSO due to the multistart strategy. However, in the small-size instances, MS-MOPSO has a shorter running time than HNSGA-II. In the mediumsize instances, the running time of MS-MOPSO is slightly longer than that of HNSGA-II. In the large-size instances, MS-MOPSO requires nearly 3000 seconds, which is also acceptable to actual BSS operators because the application scenario of BRP-BI is at night and there is enough time to calculate a better solution. Furthermore, our experiments are conducted on a notebook. For a more powerful computer, the running time of MS-MOPSO will be greatly reduced.

Statistical Analysis of Performance Metrics.
To evaluate the convergence and the diversity of multiobjective algorithms, several performance metrics are proposed, such as hypervolume (HV), inverse generation distance (IGD), and spacing metrics (SM) [45,46]. HV is a comprehensive index used to evaluate both convergence and diversity. e bigger the HV is, the better the nondominated solution set is. IGD is the average distance between the Pareto optimal front and the nondominated solution set obtained by the algorithm. e smaller the IGD is, the better the diversity and the convergence of the nondominated solution set are. SM measures how well distributed the nondominated solutions are. e smaller the SM is, the better the solutions spread along the obtained front are.
All the metrics above require a true Pareto front. Unfortunately, the true Pareto front of our biobjective BRP-BI is unknown since it is impossible to be exactly solved. To overcome this drawback, Baños et al. [45] propose a way to construct a quasi-optimal Pareto front by executing all algorithms over a very large number of iterations. After that, all the nondominated solutions are combined, discarding the dominated solutions and maintaining the nondominated ones. In our study, the same method is used to construct the quasi-optimal Pareto front for our biobjective BRP-BI. Each algorithm has been executed for 3 hours for each instance, and the composed quasi-optimal Pareto front is used to evaluate the performance of multiobjective algorithms.
To provide reliable statistics, each algorithm is run 10 times for instance 4 (small size), instance 8 (medium size), and instance 12 (large size). To calculate HV, we set the reference point as twice the objective f1 and half the objective f2, where f1 and f2 are the average value of objective values in the quasi-optimal Pareto front. e statistical results of HV, IGD, and SM obtained by using MS-MOPSO, HNSGA-II, and MOPSO are displayed in Tables 5-7, respectively. In all the three tables, the best, worst, average, and variance of the three metrics are given.   To sum up, compared with HNSGA-II and MOPSO, MS-MOPSO performs better in convergence, diversity, and distribution. Besides, MS-MOPSO is more stable. MS-MOPSO is an improved version of MOPSO and keeps the advantages of simplicity, diversity, and fast convergence. By combining the multistart strategy as well as variable neighborhood search to enhance its exploratory capability, MS-MOPSO has achieved a well-balanced exploration-exploitation capability. erefore, our algorithm can generate higher quality Pareto optimal solutions than HNSGA-II and MOPSO in the biobjective BRP-BI.

Discussion on the Value of θ.
Parameter θ directly affects the size of the balance interval and then the number of unbalanced stations that needs to be served. In other words, the change of θ will result in the change of the balance interval, which can affect whether the stations need to be rebalanced or not. Generally speaking, the larger θ is, the less the unbalanced stations are. In this section, the sensitivity analysis of θ is analyzed using the instance 8, a medium-size instance of 60 stations. Different values of θ are from 0.3 to 0 in decrement of 0.02. Table 8 displays the balance interval, the number of unbalanced stations, the unbalanced ratio, and two objectives for different θ. With the decrease of θ, the balance interval is getting narrower and narrower. us, both the number of unbalanced stations and unbalanced ratio are larger and larger. e objective g 1 (f 1 ), i.e., the total travel distance, generally increases with the decrease of θ. Especially when θ � 0, the total travel distance is the highest because the balance interval shrinks to a point; that is, the station is in equilibrium only when its inventory is just 50. Obviously, there is a fluctuation. When the value of θ is reduced from 0.22 to 0.2, the number of unbalanced stations increases from 42 to 43, but the value of g 1 (f 1 ) decreases from 80.6227 to 79.1583. is is because that the biobjective BRP-BI has both pickup demand and delivery demand, so the sequence of stations visited by the vehicle will affect the travel distance. Due to the fixed capacity of the vehicle, when the load is almost full, the vehicle may have to visit a remote delivery station. If the added station is near the current station and has a delivery demand, the travel distance can be  reduced. Conversely, when the load is almost empty, there may exist a station with pickup demand, which can also shorten the travel distance. As for the objective g 2 (− f 2 ), with the decrease of θ, the average rebalance utility also keeps the uptrend with tiny fluctuations. It can be observed that θ � 0.2 may be an appropriate value for instance 8. For other instances, the sensitivity analysis can be performed by solving a series of biobjective BRP-BI problems with different θ. According to the computation results, the suitable value of θ can be determined. However, this method cannot be widely used in real world since it is time-consuming. For the real-life BSS operators, a more realistic method is that they can use big data methods to find the appropriate value of θ from historical data.

Conclusions
We have introduced, modeled, and solved the biobjective BRP-BI, a multiobjective variant of BRP encountered in the BSS operations. We first extend the station's equilibrium from balance point to balance interval and then formulate the BRP-BI as a biobjective mixed-integer programming with the objectives of minimizing the total travel distance as well as maximizing the average rebalance utility. Furthermore, the MS-MOPSO algorithm is enhanced by the crossover operator, and variable neighborhood search is proposed to solve this problem. To illustrate the efficiency and efficacy of our algorithm, it has been tested on three types of twelve instances. e computational results reveal the huge advantage of multistart strategy and MS-MOPSO. e statistical analysis results of three performance metrics (i.e., hypervolume, inverse generation distance, and spacing metrics) also confirm the same conclusion: MS-MOPSO achieves better nondominated optimal solutions in terms of convergence, diversity, and distribution, outperforming those obtained by HNSGA-II and MOPSO. Furthermore, MS-MOPSO has better stability. e sensitivity analysis results of parameter θ indicate that as θ decreases, both objectives keep uptrend with fluctuations. e BSS decision maker can use the nondominated optimal solutions of the biobjective BRP-BI in two ways. One method is to choose an optimal rebalanced scheme from the Pareto optimal set in the actual BSS operation. To this end, specific criteria of the decision maker must be introduced. Without these criteria, all solutions in the Pareto optimal set are not comparable. After setting the preferences for objectives, the decision maker can determine an optimal scheme from the Pareto optimal set. e other method is based on the above method. e nondominated optimal solutions allow the decision maker to balance the two conflicting objectives more flexibly. By weighing these objectives, the decision maker can quickly adjust the optimal scheme to adapt to complex and dynamic practical conditions. For example, if the total travel distance must be reduced, how much the rebalance utility can be sacrificed? Or, if the rebalance utility must be greater than a given value, what is the minimum travel distance required? From the nondominated optimal solutions, the decision maker can freely determine which optimal scheme suits the current situation. Compared with the monoobjective solution method, this is an obvious advantage of our MS-MOPSO.
Future research can extend our model and algorithm to solve the BRP-BI with multiple vehicles. In addition, dynamic or stochastic BRP with balance intervals is also a good research direction.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.