# A New Intelligent Rebalancing Management Method for Multiperiod and Multiobjective Bike-Sharing System Based on Machine Learning-Enabled Signal Processing Techniques

Authors: Jianhua Cao; Weixiang Xu
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

With the rapid development of information technology, the sharing economy has developed rapidly all over the world as a new mode of distributing business profit, among which the bike-sharing system (BSS) has become popular in many cities because of its low cost, convenience, and environmental protection. The application of the 5th generation mobile communication technology (5G) in BSS makes users to search the bikes more accurately and quickly and enables operators to spot noncompliant bike sharing as soon as possible, significantly improving the efficiency of bike-sharing management. However, one of the thorny issues for operators is the bike-sharing rebalancing problem (BRP). It is the key to improve the efficiency of rebalancing, reduce the rebalancing cost, and realize the sustainable development of BSS on how to excavate the huge amount of customer cycling data, respond quickly to customer demand, and use intelligence optimization algorithm to rebalance bikes among stations. However, most of the previous studies dealt with only one period BRP and rarely considered multiperiod issues. At the same time, most researches have focused on minimizing the total cost or time of rebalancing or customer dissatisfaction, but few have aimed at minimizing the rebalancing amount. In addition, the demand gap can reflect the real rental and returning requirements of customers over a certain period of time, which is rarely considered in solving BRP. First of all, this paper presents a multiperiod and multiobjective bike-sharing rebalancing problem (MMBRP). Secondly, a mathematical model is formulated with the objective of minimizing both the total rebalancing cost and amount. In order to solve MMBRP, an improved multiobjective backtracking search genetic algorithm (IMBSGA) is designed. Finally, the effectiveness and competitiveness of IMBSGA in solving MMBRP are verified by numerous experiments comparing with state-of-the-art algorithms.


---

## Body

## 1. Introduction

In the modern computer era, with the rapid development of information technology, a large number of new economic and business models have emerged. As innovations of profit distribution and consumption pattern based on information technology, sharing economy sprang up and developed rapidly all over the world [1]. The rapid development of the sharing economy around the world, supported by new technologies such as mobile Internet, big data, and cloud computing, has become an important way to improve resource allocation efficiency. The sharing economy is characterized by formalization, high efficiency, openness, and distribution. The 5th generation mobile communication technology (5G) is a new generation of broadband mobile communication technology with the characteristics of high speed, low delay, and large connection. It is the network infrastructure to realize man-machine interconnection. The development of 5G network and its application in the sharing economy are of great significance to improve the work efficiency and the service level of the operation platform and alleviate the contradiction between supply and demand.As the typical representative of the sharing economy, bike-sharing systems (BSS) have gained worldwide popularity and have maintained their competitiveness and success in the growing business trend since they were first introduced in Amsterdam, Netherlands, in 1965 [2]. The main reason is that in recent years, the rapid development of motorized transportation and private vehicles caused a great deal of exhaust gas emissions and serious traffic jams, especially during rush hours, causing serious threats and great inconvenience to the health of city residents [3]. Air pollution ranks fourth among major mortality risk factors globally, according to official statistics, causing nearly 6.75 million premature deaths and 213 million losses of healthy lives in 2019 (Institute, 2019). Traffic jams caused urban Americans to travel extra 8.8 billion hours and purchase extra 3.3 billion gallons of fuel for a congestion cost of $166 billion [4]. Compared to motorized transport, BSS provides an alternative to short-distance travel, effectively addresses the last mile travel problem, and significantly reduces traffic accidents and congestion [5].However, during BSS operation, one of the thorny issues for operators is that some stations in the system have no enough bikes to be rented while others lack empty docks to return bikes, due to the randomness and real-time characteristics of the consumer demand for bikes. For example, at morning peaks, bikes are cycled by users from residential areas near their homes to working areas, while at evening rush hours, users start riding in the opposite directions. As a result, it is difficult to find bikes at stations near residential areas during the morning peaks or empty docks near working areas, and the opposite condition happens during evening rush hours. It should also be noted that people often rent and return bikes at different stations, which also leads to an imbalance in the system. The imbalance of the system not only leads to the dissatisfaction of users and the decline of service quality but also increases the economic losses of operators and finally affects the normal operation and sustainable development of BSS.The common way to solve this problem is the operator uses trucks to transport bikes from surplus stations to insufficient stations, which is bike-sharing rebalancing (repositioning) problem (BRP). BRP is a NP-hard problem with different optimization objectives [6]. Most of the previous literatures have usually set a single objective, such as the total repositioning cost, the total time, and the user satisfaction. Chemla et al. [7], Cruz et al. [8], and Bruck et al. [9] regarded the repositioning process as one-commodity pickup and delivery problem. Chemla et al. established a BRP model allowing stations to be visited multiple times or serve as a buffer to minimize the rebalancing cost. Cruz et al. proposed an iterated local search- (ILS-) based heuristic to solve the problem while Bruck et al. proposed three exact algorithms to solve BRP. Bulhões et al. [10] introduced a static BRP using multiple vehicles and allowing to visit the stations multiple times and presented an integer programming formulation with the objective of minimizing cost, which implemented under a branch-and-cut scheme, in addition to an iterated local search metaheuristic. Schuijbroek et al. [11] proposed a new cluster-first route-second heuristic and the polynomial-sized clustering problem to minimize makespan to solve BRP. Lu et al. [12] aimed to minimum cost and Kadri et al. [13] to minimize the overall time, both of them presented a highly effective memetic algorithm to solve this NP-hard problem. Legros [14] developed an implementable decision-support tool to minimize the rate of arrival of unsatisfied users. Ho and Szeto [15] recommended adopting the total penalty value of each station as the cost and using tabu search algorithm to determine the route and amount of repositioning. Szeto and Shui [16] target to minimize the excess total demand dissatisfaction and consider two service times including the total service time and the maximum route duration of the fleet to examine a novel set of loading and unloading strategies.Some of the literatures design multiple objectives, but they are usually weighted, which is essentially a single objective. Kloimüllner et al. [17] developed a dynamic BRP model minimizing the weight of unfulfilled demands, deviation from the target fill levels, loading instructions, and the total driving time and then solved it using Grey and Pilot construction heuristic to perform variable neighborhood search (VNS). Forma et al. (2015) aimed to minimize a weighted sum of the expected number of unserved users and the total traveling distance and proposed a 3-step mathematical programming-based heuristic to solve the problem. Szeto et al. [18] presented a model with the goal of minimizing the weighted sum of unmet user demand and repositioning time, designed a chemical reaction optimization (CRO) to determine the truck route, and proposed a subroutine to determine the loading and unloading amount at each station. Li et al. [19] focused on a new BRP that took into account multiple types of bikes and proposed a hybrid genetic algorithm to solve the problem. Liu et al. [20] regarded minimizing the weighted sum of the inconvenience level of getting bikes from the system, the total unmet demand, and the total operational time as the optimization objective, which is solved by an enhanced version of CRO. However, multiobjective weights are difficult to determine and cannot meet the demand preferences of different decision-makers.The difficulty and complexity of BRP will increase dramatically if multiple objectives are to be achieved simultaneously, so few studies have focused on this issue. However, compared with the single-objective BRP, the multiobjective BRP is more suitable to the actual situation and helpful to solve the practical problems. S. Zhu and F. Zhu [21] formulated a multiobjective integer linear programming model with the objective of maximizing the accessibility and level of service, minimizing the number of intersections and total construction cost, and solved by augmented constraint method to generate a set of nondominated solutions. Huo et al. [22] considered unbalanced vehicle distributions and high relocation costs and proposed a data-driven optimization model with the consideration of demand uncertainty to improve BSS efficiency and overall profit. Jia et al. [23] built a multiobjective optimization model and solved it through a modified multiobjective artificial bee colony algorithm. Jia et al. [24] proposed a bike-sharing rebalancing problem with balance intervals, formulated a biobjective mixed-integer programming model with the aim of determining both the minimum cost route for a single capacitated vehicle and the maximum average rebalance utility, and presented a multistart multiobjective particle swarm optimization algorithm to solve the model.Recently, BBS expanded in scale and the complexity of BRP increases rapidly. For large-scale BRP, how to excavate the huge amount of customer cycling data, respond to customer demand quickly, and rebalance promptly with intelligence optimization algorithm are very important to improve the efficiency, reduce the cost, and realize the sustainable development of the system. These problems perplex operators of the system, because if not properly solved, operators may lose their business value, which also arouses the attention of many experts and scholars. And many achievements have been made in BRP research; however, there are still a lot of insufficiencies and potential problems existing.(i)
Most studies assumed that the loading and unloading amounts of all problem stations are equal; however, they tend to be inconsistent in reality [7, 25]. Therefore, researches based on this premise maybe arbitrary and biased(ii)
Few studies have utilized nonrebalancing stations for rebalancing; in fact, many of which have spare bikes or empty docks in addition to meet their own customer demand(iii)
Single-period BRP rebalances based on the demand of one period and sets the initial number of bikes of stations at the beginning of the period. Multiperiod BRP, which covers the constantly varying demand of stations over multiple periods, only sets the initial number of bikes at the first period. The reasonable initial number of bikes of stations can effectively reduce the rebalancing amount and delay the rebalancing interval during multiple periods. However, most of the literature uses station capacity percentage [27] or ratio of rental demand to returning demand [28] to determine the initial number of bikes, which cannot dynamically reflect the characteristics of customer demand in finer time granularity(iv)
Few previous literature aimed at minimizing the rebalancing amount, but the rebalancing amount not only influences the rebalancing time but also the rebalancing workload and hence the rebalancing costsIn this paper, we introduce a multiperiod and multiobjective bike-sharing rebalancing problem, referred to as MMBRP. First, we present it as a tuple of 10 elements that describe the characteristics of stations and customer demand, respectively. Second, stations are classified into nonrebalancing stations and problem stations according to the correlation between the demand gap and the number of bikes or empty docks of stations. And then, from nonrebalancing stations, we select some stations used for rebalancing, known as coordination stations; they are further divided into loading and unloading stations. Third, we introduce a multiobjective mathematic model which minimizes the total rebalancing cost and amount, simultaneously. To solve this problem, we used an improved multiobjective backtracking search genetic algorithm (IMBSGA), which include an initial number of bike setting method based on periods of demand gap and a coordination station selection and rebalancing amount determining strategy. Finally, the competitiveness and effectiveness of the proposed model and algorithm are verified by a large number of experiments.The contributions of this paper are as follows.(i)
Based on station classification, the paper selects some stations from nonrebalancing stations and defines them as coordination stations involved in rebalancing, which not only makes full use of the station resources but also greatly reduces the rebalancing amount(ii)
Taking into account the demand gap, the initial number of bikes at station is determined to reflect the customer demand in finer time granularity, which effectively extends the rebalance interval and greatly reduces the rebalance(iii)
A multiperiod and multiobjective formulation is proposed, which can optimize the total rebalancing cost and the total rebalancing amount at the same time(iv)
An improved multiobjective backtracking search genetic algorithm is presented to solve MMBRP. The algorithm improves selection I to update the historical population through the Pareto front of current population and retains historically superior information while synthesizing the better information of the current population. In addition, the algorithm also proposes a coordination station selection and rebalancing amount determining strategy, which can quickly determine the daily coordination stations and their rebalancing amounts, thus reducing the complexity of the algorithmThe reminder of this paper is as follows. Section2 is the problem description and formulation. Section 3 presents IMBSGA to solve the problem. Section 4 presents the experimental results and analysis. Section 5 presents the conclusions and prospects.


## 2. Problem Description and Formulation

### 2.1. Station Classification

In this section, we first describe the status of a bike-sharing station, which is the basis for station classification and can be defined as a tuple consisting of ten elements: <type, capacity, bike, empty dock, rent, return, demand gap, available loading amount, available unloading amount, and rebalance amount>.Type: whether or not involve in rebalance, including problem stations, nonrebalancing stations, and coordination stationsCapacity: station capacity, denoted by cBike: number of bikes at a station, denoted by bEmpty dock: number of empty docks at a station, denoted by l, where l=c−bRent: renting demand at a station, denoted by fReturn: returning demand at a station, denoted by gDemand gap: difference in the number of bikes rented from and returned to a station, denoted by dg, where dg=f−g. If dg is positive, it represents the minimum number of bikes to meet customer demand. If dg is negative, its absolute value represents the minimum number of empty docks to meet customer demandAvailable loading amount: number of bikes that can be loaded from a station in addition to meeting its own customer demand, denoted by alAvailable unloading amount: number of bikes that can be unloaded to a station in addition to meeting its own customer demand, denoted by aulRebalance amount: number of bikes loading to or unloading from a station in rebalancing, denoted by rIn addition to meeting their own customer demands, some nonrebalancing stations have excess bikes (al) or empty docks (aul) provided for rebalancing, which not only makes full use of station resources but also greatly reduces the number of bikes trucks are carrying, thus saving the total cost.According to the correlation between the demand gap dg and the bikesb or empty docks l of a station, stations can be classified into nonrebalancing and problem stations which can further be divided into loading and unloading stations. It is worth noting that problem stations are ones that cannot meet customer demand and need to be rebalanced. The station classification is shown in equation (1) below and can be described as follows.
(i)
Ifdg>0 and dg>b at a station, which indicates that the number of bikes of the station cannot meet the rental demand of its customer, it is called a loading station belonging to problem stations with the rebalancing amount r=dg−b(ii)
If0<dg<=b at a station, which stands that the number of bikes or empty docks can meet its customer demand, it is defined as a nonrebalancing station with the rebalancing amount r=0(iii)
Ifdg<0 and absdg>l, which denotes that the empty docks of the station cannot meet the returning demand of its customer, it is referred to as an unloading station belonging to problem stations with the rebalancing amount r=dg+l(iv)
Ifdg<0 and absdg<=l, it is also defined as a nonrebalancing station with the rebalancing amount r=0(1)dg>0dg>bunloadingstations,dg≤bnonrebalancingstations,dg<0absdg>lloadingstations,absdg≤lnonrebalancingstations.In practice, the total loading and unloading amount of all problem stations is usually unequal. The general solution to this problem is to use trucks carrying bikes when they pull out of the garage, which will significantly increase the rebalancing cost. As mentioned earlier, nonrebalancing station ordinarily often has spare bikes and empty docks to be used in rebalancing, which not only makes full use of station resources but also greatly reduces the number of bikes trucks are carrying, thus saving the total cost.In this paper, some stations are selected from nonrebalancing stations to participate into rebalance, which are defined as coordination stations. Now, it is necessary to determine the available loading amount al and available unloading amount aul and then further determine the rebalancing amountr of a coordination station according to the total rebalancing amount of all problem stations.
(i)
Ifdg>0 at a coordination station, it indicates that the station needs dg bikes to meet its customer demand; thus, the available unloading amount is b−dg and the available loading amount is l(ii)
Ifdg<0 at a coordination station, it stands that the station only needs dg empty docks to meet its customer demand; the available loading amount is dg+l and the available unloading amount is bFigure1(a) is an example of station classification. As can be seen, station #1 is a loading station with the rebalancing amount r=25−20; station #2 is an unloading station with the rebalancing amount r=−30+20; station #3 is a nonrebalancing station with the rebalancing amount r=0. Figure 1(b) is an example of coordination station, here following station #3 in Figure 1(a). #3-1 describes the status of a coordination station before rebalancing, and its available rebalancing amount al and aul are 15 and 20, respectively. The total rebalancing amount of the two problem stations is 5−10=−5, meaning that 5 bikes need to be unloaded to the coordination station to fill the loading and unloading amount gap of the two problem stations. The available loading amount of the coordination station is 20>5, and then, its rebalancing amount is set to +5 and the corresponding parameters being updated, as shown in #3-2.Figure 1
Station classification and coordination station rebalancing amount determination.
(a)
Station classification(b)
Coordination station example


### 2.2. Problem Description

MMBRP can be described as follows. Given a period of timeT, BSS includes a fleet of trucks K with the capacity of e, a set of nodes N=S,P which represent the station set S and a garage P where trucks are stored, stations having the capacity c and the number of bikes b. In the first period t (a basic period is often a day), the initial number of bikes bi1 of station i is determined taking into account the demand gaps of multiperiod firstly. Next, according to the correlation of the demand gap and the bikes or empty docks of station i, all stations can be divided into nonrebalancing station set Z and problem station set Q, and the rebalancing amount of problem station ritQ can be determined. Further, to fill up the difference of the loading and unloading amount of all problem stations, some stations are selected from nonrebalancing stations to participate into the rebalance, defined as coordination station set H. Moreover, the available rebalancing amount of each coordination station is calculated and further divided into available loading amount alit and available unloading amount aulit. And then, the rebalancing amount of coordination stations can be determined according to the rebalancing amount of all problem stations. Finally, the optimal path is figured out, through which trucks depart from the garage, go through all problem stations and coordination stations, loading bikes to or unloading bikes from them, eventually back to the garage. In the second period, based on the number of bikes at the end of the first period, station classification and path calculation will continue. In subsequent periods, the process of the second period is repeated until the end of all periods and multiobjective optimization is achieved, including the rebalancing costs (the sum of traveling costs, inventory costs, and fixed costs) and the rebalancing amount of all periods.


### 2.3. Mathematical Model

The model assumes the following:(1)
The rental and returning demand at each station is known, and the amount does not vary with the number of bikes at station(2)
Rebalancing is done by the same type truck with a fixed capacity(3)
The distance between nodes of the road network is the shortest and the matrix is symmetricN: set of nodes, indexed by i=0,1,2,⋯,n, when i=0 is depot PS: set of stations, indexed by i=1,2,⋯,nQ: set of problem stations, Q∈SZ: set of nonrebalancing stations, Z∈SH: set of coordination stations, H∈ZK: set of trucks, indexed by k=0,1,2,⋯,zT: set of periods, indexed by t=0,1,2,⋯,pbit: number of bikes at the station i∈S before the rebalanceci: capacity of station i∈Sek: truckload k∈Kdgit: demand gap of station i∈Sdij: distance from node i∈N to node i∈NritQ: rebalancing amount of problem station i∈Qalitt: available loading amount of coordination station i∈Haulit: available unloading amount of coordination station i∈HC1: truck traveling costsC2: truck inventory costsC3: truck fixed costsThe decision variables are as follows:qi: number of demand gap periods of station i∈Sxijkt: a binary variable equals 1 if truck k travels from node i to node j and 0 otherwiseyit: binary variable equals 1 if the nonrebalancing station i∈Z participates in rebalancing and it is also the coordination station and 0 otherwiseritH: integer variable, the rebalancing amount of the coordination station i∈H and the absolute value is greater than 0vkt: binary variable equals 1 if truck k is used for rebalancing and 0 otherwisezijkt: a nonnegative integer variable, number of bikes carried by truck k from node i to node jThe mathematical model is as follows:(2)minC1∑t∈T∑k∈K∑i∈N∑j∈Ndijxijkt+C2∑t∈T∑k∈K∑i∈Sz0ikt+C3∑t∈T∑k∈Kvkt(3)min∑t∈T∑i∈Hrityit+∑i∈Qrit(4)s.t.bi1=α×ci+∑t=1qidgit,∀i∈S,(5)0≤bit≤ci,∀i∈S,t∈T,(6)bi,t+1=bi,t+gi,t+ri,t,∀ri,t∈ri,tQ,ri,tH,i∈Q∪H,t∈T,(7)alit=bit−git,∀git>0,i∈H,t∈T,ci−bit0+git,∀git<0,i∈H,t∈T,(8)aulit=ci−bit,∀git>0,i∈H,t∈T,bit,∀git<0,i∈H,t∈T,(9)−ulit≤rit≤lit,∀i∈H,t∈T,(10)∑k∈K∑j∈Sz0jktx0jkt+∑i∈QritQ+∑i∈HriHyit=0,t∈T,(11)∑j∈Sx0jkt=∑i∈Sxi0kt=1,∀k∈K,t∈T,(12)∑k∈K∑j∈Sz0jktx0jkt+∑i∈QritQargmini∈Zalit≥∑i∈Zyit,∀∑k∈K∑j∈Sz0jktx0jkt+∑i∈Qritq>0,(13)∑k∈K∑j∈Sz0jktx0jkt+∑i∈QritQargmaxi∈Zaulia≥∑i∈Zyit,∀∑k∈K∑j∈Sz0jktx0jkt+∑i∈Qritq<0,(14)xijkt=0,∀i,j∈N,i=j,k∈K,t∈T,(15)∑i∈Nxijkt=0,∀rjt=0,rjt∈rjtQ,rjtH,j∈Q∪H,k∈K,t∈T,1,∀rjt!=0,rjt∈rjtQ,rjtH,j∈Q∪H,k∈K,t∈T,(16)∑k∈Kvkt≥∑i∈Nxijkt=∑i∈Nxjikt,∀i≠j,j∈S,t∈T,(17)0≤∑i∈Nxijktzijkt+rjt≤e,∀rjt∈rjtQ,rjtH,j∈Q∪H,k∈K,t∈T,(18)∑i∈Nzjikt=∑i∈Nzijkt+rjt,∀i≠j,j∈S,k∈K,t∈T,(19)xijkt∈0,1,∀i∈N,j∈N,k∈K,t∈T,(20)yit=0,∀i∉H,i∈Z,t∈T,1,∀i∈H,i∈Z,t∈T,(21)rit∈Z,∀i∈H,t∈T,(22)zijkt∈N,∀i∈N,j∈N,k∈K,t∈T,(23)vkt∈0,1,∀k∈K,t∈T,(24)qi∈0,1⋯,p,∀i∈N.The objective function (2) is the total cost of rebalancing, including three costs: the first is the travel costs, the second is the truck inventory cost, and the last is the truck fixed cost. The objective function (3) is the total rebalancing amount of all periods. Constraint (4) is to determine initial number of bikes at station. Constraint (5) limits the data range for the initial number of bikes. Constraint (6) defines the number of bikes of station i at the beginning of period t. Constraints (7) and (8) define the available rebalancing amount of coordination stations. Constraint (9) is the constraint on the rebalancing amount of coordination stations. Constraint (10) ensures that the total rebalancing amount in the system is 0, that is, the sum of the number of bikes carried by truck plus the total rebalancing amount of problem stations plus the rebalancing amount of coordination stations is 0. Constraint (11) restricts the truck to return to its original depot. Constraints (12) and (13) limit the minimum number of coordination stations participating in rebalancing. Constraint (14) does not allow all intrastation visits. Constraint (15) forbids the truck to visit stations of empty demand (ri=0) and requires it visit all the stations with nonzero rebalancing amount. Constraint (16) does not allow temporary stops. Constraint (17) restricted the truck by its maximum capacity. Constraint (18) defines the flow constraint at each station. Constraints (19)–(24) define variables.


## 2.1. Station Classification

In this section, we first describe the status of a bike-sharing station, which is the basis for station classification and can be defined as a tuple consisting of ten elements: <type, capacity, bike, empty dock, rent, return, demand gap, available loading amount, available unloading amount, and rebalance amount>.Type: whether or not involve in rebalance, including problem stations, nonrebalancing stations, and coordination stationsCapacity: station capacity, denoted by cBike: number of bikes at a station, denoted by bEmpty dock: number of empty docks at a station, denoted by l, where l=c−bRent: renting demand at a station, denoted by fReturn: returning demand at a station, denoted by gDemand gap: difference in the number of bikes rented from and returned to a station, denoted by dg, where dg=f−g. If dg is positive, it represents the minimum number of bikes to meet customer demand. If dg is negative, its absolute value represents the minimum number of empty docks to meet customer demandAvailable loading amount: number of bikes that can be loaded from a station in addition to meeting its own customer demand, denoted by alAvailable unloading amount: number of bikes that can be unloaded to a station in addition to meeting its own customer demand, denoted by aulRebalance amount: number of bikes loading to or unloading from a station in rebalancing, denoted by rIn addition to meeting their own customer demands, some nonrebalancing stations have excess bikes (al) or empty docks (aul) provided for rebalancing, which not only makes full use of station resources but also greatly reduces the number of bikes trucks are carrying, thus saving the total cost.According to the correlation between the demand gap dg and the bikesb or empty docks l of a station, stations can be classified into nonrebalancing and problem stations which can further be divided into loading and unloading stations. It is worth noting that problem stations are ones that cannot meet customer demand and need to be rebalanced. The station classification is shown in equation (1) below and can be described as follows.
(i)
Ifdg>0 and dg>b at a station, which indicates that the number of bikes of the station cannot meet the rental demand of its customer, it is called a loading station belonging to problem stations with the rebalancing amount r=dg−b(ii)
If0<dg<=b at a station, which stands that the number of bikes or empty docks can meet its customer demand, it is defined as a nonrebalancing station with the rebalancing amount r=0(iii)
Ifdg<0 and absdg>l, which denotes that the empty docks of the station cannot meet the returning demand of its customer, it is referred to as an unloading station belonging to problem stations with the rebalancing amount r=dg+l(iv)
Ifdg<0 and absdg<=l, it is also defined as a nonrebalancing station with the rebalancing amount r=0(1)dg>0dg>bunloadingstations,dg≤bnonrebalancingstations,dg<0absdg>lloadingstations,absdg≤lnonrebalancingstations.In practice, the total loading and unloading amount of all problem stations is usually unequal. The general solution to this problem is to use trucks carrying bikes when they pull out of the garage, which will significantly increase the rebalancing cost. As mentioned earlier, nonrebalancing station ordinarily often has spare bikes and empty docks to be used in rebalancing, which not only makes full use of station resources but also greatly reduces the number of bikes trucks are carrying, thus saving the total cost.In this paper, some stations are selected from nonrebalancing stations to participate into rebalance, which are defined as coordination stations. Now, it is necessary to determine the available loading amount al and available unloading amount aul and then further determine the rebalancing amountr of a coordination station according to the total rebalancing amount of all problem stations.
(i)
Ifdg>0 at a coordination station, it indicates that the station needs dg bikes to meet its customer demand; thus, the available unloading amount is b−dg and the available loading amount is l(ii)
Ifdg<0 at a coordination station, it stands that the station only needs dg empty docks to meet its customer demand; the available loading amount is dg+l and the available unloading amount is bFigure1(a) is an example of station classification. As can be seen, station #1 is a loading station with the rebalancing amount r=25−20; station #2 is an unloading station with the rebalancing amount r=−30+20; station #3 is a nonrebalancing station with the rebalancing amount r=0. Figure 1(b) is an example of coordination station, here following station #3 in Figure 1(a). #3-1 describes the status of a coordination station before rebalancing, and its available rebalancing amount al and aul are 15 and 20, respectively. The total rebalancing amount of the two problem stations is 5−10=−5, meaning that 5 bikes need to be unloaded to the coordination station to fill the loading and unloading amount gap of the two problem stations. The available loading amount of the coordination station is 20>5, and then, its rebalancing amount is set to +5 and the corresponding parameters being updated, as shown in #3-2.Figure 1
Station classification and coordination station rebalancing amount determination.
(a)
Station classification(b)
Coordination station example


## 2.2. Problem Description

MMBRP can be described as follows. Given a period of timeT, BSS includes a fleet of trucks K with the capacity of e, a set of nodes N=S,P which represent the station set S and a garage P where trucks are stored, stations having the capacity c and the number of bikes b. In the first period t (a basic period is often a day), the initial number of bikes bi1 of station i is determined taking into account the demand gaps of multiperiod firstly. Next, according to the correlation of the demand gap and the bikes or empty docks of station i, all stations can be divided into nonrebalancing station set Z and problem station set Q, and the rebalancing amount of problem station ritQ can be determined. Further, to fill up the difference of the loading and unloading amount of all problem stations, some stations are selected from nonrebalancing stations to participate into the rebalance, defined as coordination station set H. Moreover, the available rebalancing amount of each coordination station is calculated and further divided into available loading amount alit and available unloading amount aulit. And then, the rebalancing amount of coordination stations can be determined according to the rebalancing amount of all problem stations. Finally, the optimal path is figured out, through which trucks depart from the garage, go through all problem stations and coordination stations, loading bikes to or unloading bikes from them, eventually back to the garage. In the second period, based on the number of bikes at the end of the first period, station classification and path calculation will continue. In subsequent periods, the process of the second period is repeated until the end of all periods and multiobjective optimization is achieved, including the rebalancing costs (the sum of traveling costs, inventory costs, and fixed costs) and the rebalancing amount of all periods.


## 2.3. Mathematical Model

The model assumes the following:(1)
The rental and returning demand at each station is known, and the amount does not vary with the number of bikes at station(2)
Rebalancing is done by the same type truck with a fixed capacity(3)
The distance between nodes of the road network is the shortest and the matrix is symmetricN: set of nodes, indexed by i=0,1,2,⋯,n, when i=0 is depot PS: set of stations, indexed by i=1,2,⋯,nQ: set of problem stations, Q∈SZ: set of nonrebalancing stations, Z∈SH: set of coordination stations, H∈ZK: set of trucks, indexed by k=0,1,2,⋯,zT: set of periods, indexed by t=0,1,2,⋯,pbit: number of bikes at the station i∈S before the rebalanceci: capacity of station i∈Sek: truckload k∈Kdgit: demand gap of station i∈Sdij: distance from node i∈N to node i∈NritQ: rebalancing amount of problem station i∈Qalitt: available loading amount of coordination station i∈Haulit: available unloading amount of coordination station i∈HC1: truck traveling costsC2: truck inventory costsC3: truck fixed costsThe decision variables are as follows:qi: number of demand gap periods of station i∈Sxijkt: a binary variable equals 1 if truck k travels from node i to node j and 0 otherwiseyit: binary variable equals 1 if the nonrebalancing station i∈Z participates in rebalancing and it is also the coordination station and 0 otherwiseritH: integer variable, the rebalancing amount of the coordination station i∈H and the absolute value is greater than 0vkt: binary variable equals 1 if truck k is used for rebalancing and 0 otherwisezijkt: a nonnegative integer variable, number of bikes carried by truck k from node i to node jThe mathematical model is as follows:(2)minC1∑t∈T∑k∈K∑i∈N∑j∈Ndijxijkt+C2∑t∈T∑k∈K∑i∈Sz0ikt+C3∑t∈T∑k∈Kvkt(3)min∑t∈T∑i∈Hrityit+∑i∈Qrit(4)s.t.bi1=α×ci+∑t=1qidgit,∀i∈S,(5)0≤bit≤ci,∀i∈S,t∈T,(6)bi,t+1=bi,t+gi,t+ri,t,∀ri,t∈ri,tQ,ri,tH,i∈Q∪H,t∈T,(7)alit=bit−git,∀git>0,i∈H,t∈T,ci−bit0+git,∀git<0,i∈H,t∈T,(8)aulit=ci−bit,∀git>0,i∈H,t∈T,bit,∀git<0,i∈H,t∈T,(9)−ulit≤rit≤lit,∀i∈H,t∈T,(10)∑k∈K∑j∈Sz0jktx0jkt+∑i∈QritQ+∑i∈HriHyit=0,t∈T,(11)∑j∈Sx0jkt=∑i∈Sxi0kt=1,∀k∈K,t∈T,(12)∑k∈K∑j∈Sz0jktx0jkt+∑i∈QritQargmini∈Zalit≥∑i∈Zyit,∀∑k∈K∑j∈Sz0jktx0jkt+∑i∈Qritq>0,(13)∑k∈K∑j∈Sz0jktx0jkt+∑i∈QritQargmaxi∈Zaulia≥∑i∈Zyit,∀∑k∈K∑j∈Sz0jktx0jkt+∑i∈Qritq<0,(14)xijkt=0,∀i,j∈N,i=j,k∈K,t∈T,(15)∑i∈Nxijkt=0,∀rjt=0,rjt∈rjtQ,rjtH,j∈Q∪H,k∈K,t∈T,1,∀rjt!=0,rjt∈rjtQ,rjtH,j∈Q∪H,k∈K,t∈T,(16)∑k∈Kvkt≥∑i∈Nxijkt=∑i∈Nxjikt,∀i≠j,j∈S,t∈T,(17)0≤∑i∈Nxijktzijkt+rjt≤e,∀rjt∈rjtQ,rjtH,j∈Q∪H,k∈K,t∈T,(18)∑i∈Nzjikt=∑i∈Nzijkt+rjt,∀i≠j,j∈S,k∈K,t∈T,(19)xijkt∈0,1,∀i∈N,j∈N,k∈K,t∈T,(20)yit=0,∀i∉H,i∈Z,t∈T,1,∀i∈H,i∈Z,t∈T,(21)rit∈Z,∀i∈H,t∈T,(22)zijkt∈N,∀i∈N,j∈N,k∈K,t∈T,(23)vkt∈0,1,∀k∈K,t∈T,(24)qi∈0,1⋯,p,∀i∈N.The objective function (2) is the total cost of rebalancing, including three costs: the first is the travel costs, the second is the truck inventory cost, and the last is the truck fixed cost. The objective function (3) is the total rebalancing amount of all periods. Constraint (4) is to determine initial number of bikes at station. Constraint (5) limits the data range for the initial number of bikes. Constraint (6) defines the number of bikes of station i at the beginning of period t. Constraints (7) and (8) define the available rebalancing amount of coordination stations. Constraint (9) is the constraint on the rebalancing amount of coordination stations. Constraint (10) ensures that the total rebalancing amount in the system is 0, that is, the sum of the number of bikes carried by truck plus the total rebalancing amount of problem stations plus the rebalancing amount of coordination stations is 0. Constraint (11) restricts the truck to return to its original depot. Constraints (12) and (13) limit the minimum number of coordination stations participating in rebalancing. Constraint (14) does not allow all intrastation visits. Constraint (15) forbids the truck to visit stations of empty demand (ri=0) and requires it visit all the stations with nonzero rebalancing amount. Constraint (16) does not allow temporary stops. Constraint (17) restricted the truck by its maximum capacity. Constraint (18) defines the flow constraint at each station. Constraints (19)–(24) define variables.


## 3. An Improved Multiobjective Backtracking Search Genetic Algorithm

At present, machine learning methods have been widely used in traffic flow prediction and station clustering in the field of bike-sharing research and achieved many performances. In this section, we develop an improved multiobjective backtracking search genetic algorithm (IMBSGA) to solve MMBRP. The main reason for adopting this algorithm is that BRP itself is a NP-hard problem and genetic algorithms (GA) are an effective method to solve such problem. In addition, backtracking search algorithm (BSA) [26] can retain the historical information for a certain period of time, thus increasing population diversity and avoiding falling into local optimum. On this basis, we improve the algorithm, first of all, coding the number of demand gap periods and then calculating the initial number of bikes of stations. Further, we propose a strategy to determine the coordination stations and their rebalancing amount. And then, we update historical population based on the Pareto front in the first selection stage (selection I) of BSA, which ensures the optimal solutions of current population can be recorded in historical populations.The algorithm flow is as shown in Figure2.Figure 2
Algorithm flow.


### 3.1. Coding

In this paper, each cell on a chromosome represents the cumulative demand gap periods at a station. Figure3 is an example of a chromosome structure consisting of six stations, each location representing a station and the figure denoting the cumulative demand gap periods at the station. Therefore, the cumulative demand gap periods of the six stations are 3, 4, 2, 7, 5, and 3, respectively. For instance, 3 is the cumulative demand gap periods of the first station, similarly as the following figure.Figure 3
An example of a chromosome structure.


### 3.2. Calculating the Initial Number of Bikes Based on the Periods of Demand Gap

In this section, we propose a method for calculating the initial number of bikes considering the demand gap periods, which can reflect the trend of rental and returning demand of users for a certain period of time. According to the cumulative demand gaps of multiple periods, the initial number of bikes can be determined to extent the rebalance interval and reduce the rebalance amount. The formula is as equation (4). The initial number of bikes bi1 comprises two parts. The first part is determined based on the percentage capacity of station i [27]. And then, the number of demand gap periods is obtained by coding and the sum of demand gaps periods of station i during period zi is calculated, which form the second part.


### 3.3. Determining Strategy of Coordination Station Selection and Rebalancing Amount

This paper innovatively proposes the use of coordination stations to solve MMBRP; however, the problem station set is determined according to the initial number of bikes and customer demand of the period. As the rebalance proceeds, the problem station set of the next period will change with the set of the period before, which becomes a dynamic problem. That is, only determining the initial number of bikes and the rebalancing amount of each problem and coordination station of the first period can the problem station for the second period be identified and then coordination stations be selected. By this means, the problem station set and coordination station set of all periods can be obtained. Therefore, in order to reduce the complexity of the algorithm, we propose a coordination station selection and rebalancing amount determining strategy, which can directly determine the coordination stations and their rebalancing amount only based on the initial number of bikes in the first period and then obtain the problem stations and coordination stations of all periods.The process of strategy is as follows.Step 1.
Determine the problem station setQ based on the initial number of bikes in the first period and calculate the total rebalancing amount of problem stations R=∑i∈Qri.Step 2.
Each problem station selects the two closest nonrebalancing stations to join the coordination station setH and computes the available rebalance amount mi.Step 3.
Select the coordination stationi closest to the problem station; if mi<=R, the rebalancing amount ri=mi; otherwise, ri=R and R=R−ri.Step 4.
IfR=0, the strategy terminates. Otherwise, repeat Step 3.


### 3.4. Initialization

Initialization is the first step in IMBSGA, and the first population is generated during this process. The value of each cell in the chromosome is randomly generated within the range1,m, of which m is the maximum value of periods. When M chromosomes are generated, the initialization end and current population P generate, while replicating P to generate history population OldP.


### 3.5. Selection I (Population Selection Mechanism Based on Pareto Front)

At this stage, the algorithm determines the historical population by selection strategy, which will lead the population evolution in the selection step. Before determining the historical population, the algorithm uses equation (24) to determine whether to update the population in this iteration, where a and b are random numbers generated from an even distribution at 0,1. If the historical population of this iteration is updated, the algorithm used equation (26) to update OldP, which first places the nondominated solution set in the current population P into OldP and then remakes the Pareto selection to produce OldP of N size and changes the order of the individuals in it at random. The mechanism ensures that the BSA can retain historically superior information while synthesizing the better information of the current population.
(25)OldP=UpdateP,a<b,OldP,otherwise,(26)UpdateP=ParetoSelectOldP+nondominatedsolutionsetofP.


### 3.6. Crossover and Mutation

Crossover operation, which combines some characteristics of different parent individuals, is the main way to generate new individuals. In this paper, we adopt a uniform crossover strategy of intersecting parent generations, one of which is current population, and the other is historical population. It exchanges the crossover point on each parent according to probabilities, thus creating two new individuals. The general steps of this strategy are as follows:(1)
Randomly select two chromosomes as parents(2)
Swap each crossover point in the parents according to probabilityThe process of uniform crossing is shown in Figure4.Figure 4
The process of uniform crossing.This paper uses a uniform mutation method, which used a random number in the range of1,m to replace the original gene value at each gene in the individual with a small probability.


### 3.7. Selection II

Different from traditional BSA, all the solutions are based on nondominated sorting technique in IMBSGA. So we combine the generated populationA with the current population B and select the next generation from A+B according to the Pareto hierarchy.


## 3.1. Coding

In this paper, each cell on a chromosome represents the cumulative demand gap periods at a station. Figure3 is an example of a chromosome structure consisting of six stations, each location representing a station and the figure denoting the cumulative demand gap periods at the station. Therefore, the cumulative demand gap periods of the six stations are 3, 4, 2, 7, 5, and 3, respectively. For instance, 3 is the cumulative demand gap periods of the first station, similarly as the following figure.Figure 3
An example of a chromosome structure.


## 3.2. Calculating the Initial Number of Bikes Based on the Periods of Demand Gap

In this section, we propose a method for calculating the initial number of bikes considering the demand gap periods, which can reflect the trend of rental and returning demand of users for a certain period of time. According to the cumulative demand gaps of multiple periods, the initial number of bikes can be determined to extent the rebalance interval and reduce the rebalance amount. The formula is as equation (4). The initial number of bikes bi1 comprises two parts. The first part is determined based on the percentage capacity of station i [27]. And then, the number of demand gap periods is obtained by coding and the sum of demand gaps periods of station i during period zi is calculated, which form the second part.


## 3.3. Determining Strategy of Coordination Station Selection and Rebalancing Amount

This paper innovatively proposes the use of coordination stations to solve MMBRP; however, the problem station set is determined according to the initial number of bikes and customer demand of the period. As the rebalance proceeds, the problem station set of the next period will change with the set of the period before, which becomes a dynamic problem. That is, only determining the initial number of bikes and the rebalancing amount of each problem and coordination station of the first period can the problem station for the second period be identified and then coordination stations be selected. By this means, the problem station set and coordination station set of all periods can be obtained. Therefore, in order to reduce the complexity of the algorithm, we propose a coordination station selection and rebalancing amount determining strategy, which can directly determine the coordination stations and their rebalancing amount only based on the initial number of bikes in the first period and then obtain the problem stations and coordination stations of all periods.The process of strategy is as follows.Step 1.
Determine the problem station setQ based on the initial number of bikes in the first period and calculate the total rebalancing amount of problem stations R=∑i∈Qri.Step 2.
Each problem station selects the two closest nonrebalancing stations to join the coordination station setH and computes the available rebalance amount mi.Step 3.
Select the coordination stationi closest to the problem station; if mi<=R, the rebalancing amount ri=mi; otherwise, ri=R and R=R−ri.Step 4.
IfR=0, the strategy terminates. Otherwise, repeat Step 3.


## 3.4. Initialization

Initialization is the first step in IMBSGA, and the first population is generated during this process. The value of each cell in the chromosome is randomly generated within the range1,m, of which m is the maximum value of periods. When M chromosomes are generated, the initialization end and current population P generate, while replicating P to generate history population OldP.


## 3.5. Selection I (Population Selection Mechanism Based on Pareto Front)

At this stage, the algorithm determines the historical population by selection strategy, which will lead the population evolution in the selection step. Before determining the historical population, the algorithm uses equation (24) to determine whether to update the population in this iteration, where a and b are random numbers generated from an even distribution at 0,1. If the historical population of this iteration is updated, the algorithm used equation (26) to update OldP, which first places the nondominated solution set in the current population P into OldP and then remakes the Pareto selection to produce OldP of N size and changes the order of the individuals in it at random. The mechanism ensures that the BSA can retain historically superior information while synthesizing the better information of the current population.
(25)OldP=UpdateP,a<b,OldP,otherwise,(26)UpdateP=ParetoSelectOldP+nondominatedsolutionsetofP.


## 3.6. Crossover and Mutation

Crossover operation, which combines some characteristics of different parent individuals, is the main way to generate new individuals. In this paper, we adopt a uniform crossover strategy of intersecting parent generations, one of which is current population, and the other is historical population. It exchanges the crossover point on each parent according to probabilities, thus creating two new individuals. The general steps of this strategy are as follows:(1)
Randomly select two chromosomes as parents(2)
Swap each crossover point in the parents according to probabilityThe process of uniform crossing is shown in Figure4.Figure 4
The process of uniform crossing.This paper uses a uniform mutation method, which used a random number in the range of1,m to replace the original gene value at each gene in the individual with a small probability.


## 3.7. Selection II

Different from traditional BSA, all the solutions are based on nondominated sorting technique in IMBSGA. So we combine the generated populationA with the current population B and select the next generation from A+B according to the Pareto hierarchy.


## 4. Experimental Results and Analysis

The information management system is the core component of BSS, which is responsible for user information management, real-time data monitoring, and other functions. Note that the information management for bike sharing refers to the terminal information monitoring of bike-sharing stations, including station layout, road conditions, the distribution, and bike-sharing renting and returning information. In this paper, signal processing technology is used for data collection from Citi Bike System in New York, which can be accessed from here:https://s3.amazonaws.com/tripdata/201903-citibike-tripdata.csv.zip.Citi Bike is the first bike-sharing program in New York, which was founded in September 2011 after winning a bid of the bike-sharing program organized by the New York Department of Transportation. Citi Bike is a typical large-scale bike-sharing system with an initial deployment of 6000 bikes and 300 stations. Since then, the scale of Citi Bike system continues to expand and it has grown to 12,000 bikes and 773 stations as of March 2019.The data must be cleaned and processed before they can be used for analysis and experiment. Due to the less missing and good versatility of the data after preprocessing, many experts and scholars apply it to the study of large-scale BRP problems. In this paper, the location, capacity, and customer demand of 770 stations from March 1 to 3, 2019, have been obtained.All the methods are programmed with Python 3.6 and implemented on an Intel(R) Core (TM) i9-9900K CPU @ 3.60 GHz 64 G computer equipped with Ubuntu operating system. The case parameters of this paper are as follows: truck routing costsC1=0.001, truck inventory cost C2=0.5, and truck capacity is 45, corresponding to fixed truck costs C3=20. The number of trucks K=1. The parameters of the IMBSGA are as follows: the population scale CN=50, the cyclic iteration maximum maxcycle=50, cross rate cr=0.8, and mutation rate mr=0.2.


### 4.1. Comparison with Other Setting Methods of Initial Number of Bikes

In order to verify the validity of the method proposed in this paper, the previous research results are selected for comparison, including setting the initial number of bikes according to the station capacity percentage (SCP) [27] and ratio of rental demand to return demand (RRDRD) [28]; the results are shown in Table 1.Table 1
The total rebalancing amount and cost of IMBSGA, SCP, and RRDRD.
MethodThe total rebalancing amountThe total costIMBSGA (cost minimization)106116.9IMBSGA (rebalancing amount minimization)96126.46SCP (50%)769213.664RRDRD565170.59As can be seen from Table1, the total rebalancing amount or the total cost obtained using IMBSGA is much lower than that of other methods, and the results calculated by RRDRD are better than those of SCP. These results indicate that setting the initial number of bikes of stations considering the demand characteristics of customer can effectively reduce the total rebalancing amount and cost, and the utilization of demand gap is more effective.To further verify the effectiveness of IMBSGA, we compare the number of problem stations (NPS), the number of coordination stations (NCS), the rebalancing amount (RA), and the rebalancing cost (RC) obtained by the three methods over three days (the basic study period is one day), as shown in Table2.Table 2
The NPS, NCS, RA, and RC of IMBSGA, SCP, and RRDRD.
DayIMBSGASCPRRDRDNPSNCSRARCNPSNCSRARCNPSNCSRARC1733732.122912264.071569351.122541535.29411519283.48261113967.5531274459.077439455106.12522933391.92Sum241496126.4613763769253.679346565210.59NPS is an important index of BRP. In general, the fewer NPS there is, the smaller the RA and therefore the lower the NCS. In addition, the fewer the total number of stations participating in rebalancing NPS+NCS there is, the lower the RC. Table2 shows that the NPS in IMBSGA, SCP, and RRDRD are 7, 22, and 15 on the first day, respectively, and the result obtained by IMBSGA is far lower than the results obtained by the other two methods, the same as the indices of NCS and RA. Because the total number of stations that participated in rebalancing on the first day in IMBSGA is only 7+3, the total rebalancing cost is greatly reduced compared to SCP and RRDRD. Further, the NPS, NCS, RA, and RC calculated by IMBSGA of the second and third days are also far below than those of SCP and RRDRD. From the view of the sum of three days, the NPS in IMBSGA is only 24, which is much smaller than that in the other methods; the index in SCP and RRDRD is 137 and 93, respectively, the same as the remaining three indices. All indices in IMBSGA are significantly lower than those in SCP and RRDRD, which indicates that the superiority of IMBSGA is more obvious. Moreover, it illustrates that the initial number of bikes obtained by the proposed method in this paper can significantly extend the rebalancing interval and decrease the NPS per day, thus reducing the workload of rebalancing.


### 4.2. Comparison with Other Well-Known Multiobjective Optimization Algorithms

The Nondominated Sorting Genetic Algorithm II (NSGA-II) is an efficient multiobjective genetic algorithm that uses nondominated sorting and crowding distance to solve Pareto optimal solutions. It is also one of the well-regarded variants of GA, and its chromosome rank is based on their level domination. Compared to the earlier variants of GA (NSGA), NSGA-II achieves a “fast nondominating sorting technique,” a “crowding distance technique,” and a parameterless niching operator, greatly improving the complexity and accuracy of the algorithm [29, 30].SPEA2 is an evolutionary method for multicriteria optimization proposed by Eckart Zitzler and Lothar Thiele in 1999, also known as enhanced Pareto genetic algorithm. In contrast to the original Strong Pareto Evolutionary Algorithm (SPEA), the improvements made in SPEA2 are related to the increased detail in the calculation of individual aptitude, the size of the external file, and the parent selection process. First, the aptitude calculation differences in that for every solution, both dominating and dominated individuals are considered. Additionally, a new parameter called “density” is added that establishes the proximity of every individual to its neighbors. Finally, parent selection is conducted from the external file only, which maintains a permanent size throughout algorithm execution until the maximum number of generations (T) has elapsed or another termination criterion is reached [31].


#### 4.2.1. The Performance of Pareto Front

In this paper, we study a multiobjective optimization problem (MOP), whose solution is to find a compromise solution that approximately optimizes multiple objectives simultaneously, namely, the Pareto optimal solution set (PS). The surface in the target space formed by the objective function corresponding to all Pareto solutions in the PS is called the Pareto front (PF).Table3 shows the nondominant solution set of 50 population using IMBSGA. As can be seen from the experimental results, the rebalancing amount is distributed between 96 and 106 and the total cost distributes between 116.9 and 126.46, among which the smallest rebalancing amount solution is [96, 126.46] and the lowest total cost solution is [106, 116.9].Table 3
Nondominant solution set of IMBSGA.
IDThe total rebalancing amountThe total rebalancing cost1103117.932106116.41396126.464105116.95104117.456105116.48798121.31899120.69The performances of IMBSGA, NSGA-II, and SPEA2 [31] were compared based on the analysis of Pareto front, and the nondominant solution set of these three algorithms is shown in Figure 5. The results show that the nondominated solution set distribution on the Pareto front of IMBSGA is the most evenly distributed, suggesting that IMBSGA is superior to other algorithms. First, the historical population of IMBSGA increases the ability to exploit the solution space and maintains historical information of the population, thus ensuring the diversity of population for future generations. In addition, IMBSGA also extends the search scope of solution space to avoid falling into local optimum and instability in the process of evolution. Finally, the improved selection I makes it possible for historical population to record nondominant solution set of current population, thus enhancing the iterative speed of algorithm so that better solutions can be found more quickly.Figure 5
Nondominant solution set of IMBSGA, NSGA-II, and SPEA2.


#### 4.2.2. Spread Metric

The spread metric is used to measure the extent of the diffusion and distribution of solutions in approximate concentrations, and it is measured by the following:(27)spread=df+dl+∑i=1N−1di−d¯df+dl+N−1×d¯,whereN is the scale of the nondominated set, di is the Euclidean distance between consecutive solutions, d¯ is the average distance of di, and df and dl are the Euclidean distances between the extreme solutions and the boundary solutions of the nondominated solution set. The more uniformly the solution set is distributed, the closer ∑i=1N−1di−d¯ will close to zero; the wider the distribution of the solution set, the closer df and dl will close to zero. Thus, the smaller the spread, the better the distribution of the nondominated solution set.Using three algorithms of IMBSGA, NSGA-II, and SPEA2, the tests were independently conducted in 30, 35, 40, 45, and 50 iterations. Statistical performances such as the mean and standard deviation (std) of spread under different iterations are shown in Table4. Obviously, the mean of NASG-II is the worst, far greater than the other two algorithms. The metric of SPEA2 is not stable, and the std is higher than the other two algorithms, but still better than NASG-II in diversity. The mean and std of IMBSGA are much smaller than those of NASG-II and SPEA2, indicating that IMBSGA has better diversity and more stable performance.Table 4
Spread of IMBSGA, NSGA-II, and SPEA2.
IterationsIMBSGANASG-IISPEA2MeanstdMeanstdMeanstd300.652870.022761.143990.030940.987640.04992350.639210.0184311.13740.021131.016870.05973400.651790.020941. 123880.030820.988470.05251450.643670.022421.133820.033030.971450.03618500.636990.0194441.130260.023021.006230.04666


#### 4.2.3. Inverted Generational Distance (IGD)

This metric calculates the average Euclidean distance from the solution in the Pareto front to their closest solution in the approximation set, which can be calculated as follows:(28)IGD=∑i=1Ndi∗2N.di∗ is the minimum Euclidean distance from the i-th solution to the Pareto front. The lower the IGD value, the better the metric, indicating that the solution set has good convergence and uniformity.With 50 iterations, 10 independent experiments were carried out under three algorithms of IMBSGA, NSGA-II, and SPEA2, and the statistical properties of IGD such as mean and standard deviation (std) are shown in Table5.Table 5
IGD of IMBSGA, NSGA-II, and SPEA2.
IMBSGANASG-IISPEA2MeanstdMeanstdMeanstd0.249090.048620.258630.024570.281860.03294As can be seen from Table5, the mean and std of IGD obtained by IMBSGA are smaller than those obtained by NSGA-II and SPEA2, suggesting that while improving diversity, convergence can be improved along with diversity.


## 4.1. Comparison with Other Setting Methods of Initial Number of Bikes

In order to verify the validity of the method proposed in this paper, the previous research results are selected for comparison, including setting the initial number of bikes according to the station capacity percentage (SCP) [27] and ratio of rental demand to return demand (RRDRD) [28]; the results are shown in Table 1.Table 1
The total rebalancing amount and cost of IMBSGA, SCP, and RRDRD.
MethodThe total rebalancing amountThe total costIMBSGA (cost minimization)106116.9IMBSGA (rebalancing amount minimization)96126.46SCP (50%)769213.664RRDRD565170.59As can be seen from Table1, the total rebalancing amount or the total cost obtained using IMBSGA is much lower than that of other methods, and the results calculated by RRDRD are better than those of SCP. These results indicate that setting the initial number of bikes of stations considering the demand characteristics of customer can effectively reduce the total rebalancing amount and cost, and the utilization of demand gap is more effective.To further verify the effectiveness of IMBSGA, we compare the number of problem stations (NPS), the number of coordination stations (NCS), the rebalancing amount (RA), and the rebalancing cost (RC) obtained by the three methods over three days (the basic study period is one day), as shown in Table2.Table 2
The NPS, NCS, RA, and RC of IMBSGA, SCP, and RRDRD.
DayIMBSGASCPRRDRDNPSNCSRARCNPSNCSRARCNPSNCSRARC1733732.122912264.071569351.122541535.29411519283.48261113967.5531274459.077439455106.12522933391.92Sum241496126.4613763769253.679346565210.59NPS is an important index of BRP. In general, the fewer NPS there is, the smaller the RA and therefore the lower the NCS. In addition, the fewer the total number of stations participating in rebalancing NPS+NCS there is, the lower the RC. Table2 shows that the NPS in IMBSGA, SCP, and RRDRD are 7, 22, and 15 on the first day, respectively, and the result obtained by IMBSGA is far lower than the results obtained by the other two methods, the same as the indices of NCS and RA. Because the total number of stations that participated in rebalancing on the first day in IMBSGA is only 7+3, the total rebalancing cost is greatly reduced compared to SCP and RRDRD. Further, the NPS, NCS, RA, and RC calculated by IMBSGA of the second and third days are also far below than those of SCP and RRDRD. From the view of the sum of three days, the NPS in IMBSGA is only 24, which is much smaller than that in the other methods; the index in SCP and RRDRD is 137 and 93, respectively, the same as the remaining three indices. All indices in IMBSGA are significantly lower than those in SCP and RRDRD, which indicates that the superiority of IMBSGA is more obvious. Moreover, it illustrates that the initial number of bikes obtained by the proposed method in this paper can significantly extend the rebalancing interval and decrease the NPS per day, thus reducing the workload of rebalancing.


## 4.2. Comparison with Other Well-Known Multiobjective Optimization Algorithms

The Nondominated Sorting Genetic Algorithm II (NSGA-II) is an efficient multiobjective genetic algorithm that uses nondominated sorting and crowding distance to solve Pareto optimal solutions. It is also one of the well-regarded variants of GA, and its chromosome rank is based on their level domination. Compared to the earlier variants of GA (NSGA), NSGA-II achieves a “fast nondominating sorting technique,” a “crowding distance technique,” and a parameterless niching operator, greatly improving the complexity and accuracy of the algorithm [29, 30].SPEA2 is an evolutionary method for multicriteria optimization proposed by Eckart Zitzler and Lothar Thiele in 1999, also known as enhanced Pareto genetic algorithm. In contrast to the original Strong Pareto Evolutionary Algorithm (SPEA), the improvements made in SPEA2 are related to the increased detail in the calculation of individual aptitude, the size of the external file, and the parent selection process. First, the aptitude calculation differences in that for every solution, both dominating and dominated individuals are considered. Additionally, a new parameter called “density” is added that establishes the proximity of every individual to its neighbors. Finally, parent selection is conducted from the external file only, which maintains a permanent size throughout algorithm execution until the maximum number of generations (T) has elapsed or another termination criterion is reached [31].


### 4.2.1. The Performance of Pareto Front

In this paper, we study a multiobjective optimization problem (MOP), whose solution is to find a compromise solution that approximately optimizes multiple objectives simultaneously, namely, the Pareto optimal solution set (PS). The surface in the target space formed by the objective function corresponding to all Pareto solutions in the PS is called the Pareto front (PF).Table3 shows the nondominant solution set of 50 population using IMBSGA. As can be seen from the experimental results, the rebalancing amount is distributed between 96 and 106 and the total cost distributes between 116.9 and 126.46, among which the smallest rebalancing amount solution is [96, 126.46] and the lowest total cost solution is [106, 116.9].Table 3
Nondominant solution set of IMBSGA.
IDThe total rebalancing amountThe total rebalancing cost1103117.932106116.41396126.464105116.95104117.456105116.48798121.31899120.69The performances of IMBSGA, NSGA-II, and SPEA2 [31] were compared based on the analysis of Pareto front, and the nondominant solution set of these three algorithms is shown in Figure 5. The results show that the nondominated solution set distribution on the Pareto front of IMBSGA is the most evenly distributed, suggesting that IMBSGA is superior to other algorithms. First, the historical population of IMBSGA increases the ability to exploit the solution space and maintains historical information of the population, thus ensuring the diversity of population for future generations. In addition, IMBSGA also extends the search scope of solution space to avoid falling into local optimum and instability in the process of evolution. Finally, the improved selection I makes it possible for historical population to record nondominant solution set of current population, thus enhancing the iterative speed of algorithm so that better solutions can be found more quickly.Figure 5
Nondominant solution set of IMBSGA, NSGA-II, and SPEA2.


### 4.2.2. Spread Metric

The spread metric is used to measure the extent of the diffusion and distribution of solutions in approximate concentrations, and it is measured by the following:(27)spread=df+dl+∑i=1N−1di−d¯df+dl+N−1×d¯,whereN is the scale of the nondominated set, di is the Euclidean distance between consecutive solutions, d¯ is the average distance of di, and df and dl are the Euclidean distances between the extreme solutions and the boundary solutions of the nondominated solution set. The more uniformly the solution set is distributed, the closer ∑i=1N−1di−d¯ will close to zero; the wider the distribution of the solution set, the closer df and dl will close to zero. Thus, the smaller the spread, the better the distribution of the nondominated solution set.Using three algorithms of IMBSGA, NSGA-II, and SPEA2, the tests were independently conducted in 30, 35, 40, 45, and 50 iterations. Statistical performances such as the mean and standard deviation (std) of spread under different iterations are shown in Table4. Obviously, the mean of NASG-II is the worst, far greater than the other two algorithms. The metric of SPEA2 is not stable, and the std is higher than the other two algorithms, but still better than NASG-II in diversity. The mean and std of IMBSGA are much smaller than those of NASG-II and SPEA2, indicating that IMBSGA has better diversity and more stable performance.Table 4
Spread of IMBSGA, NSGA-II, and SPEA2.
IterationsIMBSGANASG-IISPEA2MeanstdMeanstdMeanstd300.652870.022761.143990.030940.987640.04992350.639210.0184311.13740.021131.016870.05973400.651790.020941. 123880.030820.988470.05251450.643670.022421.133820.033030.971450.03618500.636990.0194441.130260.023021.006230.04666


### 4.2.3. Inverted Generational Distance (IGD)

This metric calculates the average Euclidean distance from the solution in the Pareto front to their closest solution in the approximation set, which can be calculated as follows:(28)IGD=∑i=1Ndi∗2N.di∗ is the minimum Euclidean distance from the i-th solution to the Pareto front. The lower the IGD value, the better the metric, indicating that the solution set has good convergence and uniformity.With 50 iterations, 10 independent experiments were carried out under three algorithms of IMBSGA, NSGA-II, and SPEA2, and the statistical properties of IGD such as mean and standard deviation (std) are shown in Table5.Table 5
IGD of IMBSGA, NSGA-II, and SPEA2.
IMBSGANASG-IISPEA2MeanstdMeanstdMeanstd0.249090.048620.258630.024570.281860.03294As can be seen from Table5, the mean and std of IGD obtained by IMBSGA are smaller than those obtained by NSGA-II and SPEA2, suggesting that while improving diversity, convergence can be improved along with diversity.


## 4.2.1. The Performance of Pareto Front

In this paper, we study a multiobjective optimization problem (MOP), whose solution is to find a compromise solution that approximately optimizes multiple objectives simultaneously, namely, the Pareto optimal solution set (PS). The surface in the target space formed by the objective function corresponding to all Pareto solutions in the PS is called the Pareto front (PF).Table3 shows the nondominant solution set of 50 population using IMBSGA. As can be seen from the experimental results, the rebalancing amount is distributed between 96 and 106 and the total cost distributes between 116.9 and 126.46, among which the smallest rebalancing amount solution is [96, 126.46] and the lowest total cost solution is [106, 116.9].Table 3
Nondominant solution set of IMBSGA.
IDThe total rebalancing amountThe total rebalancing cost1103117.932106116.41396126.464105116.95104117.456105116.48798121.31899120.69The performances of IMBSGA, NSGA-II, and SPEA2 [31] were compared based on the analysis of Pareto front, and the nondominant solution set of these three algorithms is shown in Figure 5. The results show that the nondominated solution set distribution on the Pareto front of IMBSGA is the most evenly distributed, suggesting that IMBSGA is superior to other algorithms. First, the historical population of IMBSGA increases the ability to exploit the solution space and maintains historical information of the population, thus ensuring the diversity of population for future generations. In addition, IMBSGA also extends the search scope of solution space to avoid falling into local optimum and instability in the process of evolution. Finally, the improved selection I makes it possible for historical population to record nondominant solution set of current population, thus enhancing the iterative speed of algorithm so that better solutions can be found more quickly.Figure 5
Nondominant solution set of IMBSGA, NSGA-II, and SPEA2.


## 4.2.2. Spread Metric

The spread metric is used to measure the extent of the diffusion and distribution of solutions in approximate concentrations, and it is measured by the following:(27)spread=df+dl+∑i=1N−1di−d¯df+dl+N−1×d¯,whereN is the scale of the nondominated set, di is the Euclidean distance between consecutive solutions, d¯ is the average distance of di, and df and dl are the Euclidean distances between the extreme solutions and the boundary solutions of the nondominated solution set. The more uniformly the solution set is distributed, the closer ∑i=1N−1di−d¯ will close to zero; the wider the distribution of the solution set, the closer df and dl will close to zero. Thus, the smaller the spread, the better the distribution of the nondominated solution set.Using three algorithms of IMBSGA, NSGA-II, and SPEA2, the tests were independently conducted in 30, 35, 40, 45, and 50 iterations. Statistical performances such as the mean and standard deviation (std) of spread under different iterations are shown in Table4. Obviously, the mean of NASG-II is the worst, far greater than the other two algorithms. The metric of SPEA2 is not stable, and the std is higher than the other two algorithms, but still better than NASG-II in diversity. The mean and std of IMBSGA are much smaller than those of NASG-II and SPEA2, indicating that IMBSGA has better diversity and more stable performance.Table 4
Spread of IMBSGA, NSGA-II, and SPEA2.
IterationsIMBSGANASG-IISPEA2MeanstdMeanstdMeanstd300.652870.022761.143990.030940.987640.04992350.639210.0184311.13740.021131.016870.05973400.651790.020941. 123880.030820.988470.05251450.643670.022421.133820.033030.971450.03618500.636990.0194441.130260.023021.006230.04666


## 4.2.3. Inverted Generational Distance (IGD)

This metric calculates the average Euclidean distance from the solution in the Pareto front to their closest solution in the approximation set, which can be calculated as follows:(28)IGD=∑i=1Ndi∗2N.di∗ is the minimum Euclidean distance from the i-th solution to the Pareto front. The lower the IGD value, the better the metric, indicating that the solution set has good convergence and uniformity.With 50 iterations, 10 independent experiments were carried out under three algorithms of IMBSGA, NSGA-II, and SPEA2, and the statistical properties of IGD such as mean and standard deviation (std) are shown in Table5.Table 5
IGD of IMBSGA, NSGA-II, and SPEA2.
IMBSGANASG-IISPEA2MeanstdMeanstdMeanstd0.249090.048620.258630.024570.281860.03294As can be seen from Table5, the mean and std of IGD obtained by IMBSGA are smaller than those obtained by NSGA-II and SPEA2, suggesting that while improving diversity, convergence can be improved along with diversity.


## 5. Conclusions and Future Work

In this paper, a multiperiod and multiobjective bike-sharing rebalancing problem (MMBRP) is proposed first, which considers how to solve the multiobjective rebalancing problem under multiple periods. Then, a mathematical model is established to minimize the total rebalancing cost and amount, simultaneously. Aiming at the problem, an improved multiobjective backtracking search genetic algorithm (IMBSGA) is designed. In the algorithm, the historical population is updated by selecting the nondominated solution set of the current population in selection I, and an initial number of bike setting method based on the periods of demand gap and a coordination station selection and rebalancing amount determining strategy are proposed. Finally, the proposed method is compared with two other well-known multiobjective algorithms, NSGA-II and SPEA2, by conducting a real-world case study. In the initial number of bike setting, the methods of other papers are used for comparison. The experimental results show that this algorithm is competitive in solving MMBRP, effectively reducing the total rebalancing cost and amount in multiple periods.At present, the selection of coordination stations and the determination of rebalancing amount are still obtained through the proposed strategy. In the future, we plan to use dynamic programming methods to solve the problem. At the same time, the demands of user are known, and we plan to forecast the demand to solve the MMBRP with uncertain demand.



---

*Source: 1556467-2022-01-18.xml*