Optimization Strategies for Dockless Bike Sharing Systems via two Algorithms of Closed Queuing Networks

: The dockless bike sharing system (DBSS) has been globally adopted as a sustainable transportation system. Due to the robustness and tractability of the closed queuing network (CQN), it is a well-behaved method to model DBSSs. In this paper, we view DBSSs as CQNs and use the mean value analysis (MVA) algorithm to calculate a small size DBSS and the ﬂow equivalent server (FES) algorithm to calculate the larger size DBSS. This is the ﬁrst time that the FES algorithm is used to study the DBSS, by which the CQN can be divided into different subnetworks. A parking region and its downlink roads are viewed as a subnetwork, so the computation of CQN is reduced greatly. Based on the computation results of the two algorithms, we propose two optimization functions for determining the optimal ﬂeet size and repositioning ﬂow, respectively. At last, we provide numerical experiments to verify the two algorithms and illustrate the optimal ﬂeet size and repositioning ﬂow. This computation framework can also be used to analyze other on-demand transportation networks.


Introduction
Bike sharing systems (BSSs) are a type of smart transportation. Nowadays, there are primarily two types of BSSs. The first one, the station-based bike sharing system (SBBSS), consists of stations distributed across the city. The capacity of a station is limited. An arriving user rents a bike from a station and returns it to a destination station if there is a vacant dock. Otherwise, the user must find another station to return the bike. SBBSSs have been carried out in over 2900 cities around the world ( [1]). Another type of BSS, with no fixed stations for parking bikes, is the dockless bike sharing system (DBSS). Bikes can be parked everywhere allowable without capacity limit (refer to Shi et al. [2]). The DBSSs first arose in 2014 and now they have been deployed in over 250 cities in the world. Since the difference of capacity and parking positions, the modeling and computation of the SBBSS and DBSS are their distinguishing features. In this paper, we provide research for the DBSS.
Since the BSS is large scale and dynamic, it is vital to determine an optimal fleet size. The optimal fleet size is a tradeoff between the earnings from users and maintaining cost of the BSS. Bike availability is defined to be the probability that the parking region is not empty. It is an essential task to understand the relationship between the bike availability of each parking region and the fleet size. If the fleet is in surplus, it is a waste of resource and many bikes will pile up at the roadside and cause traffic congestion. It is a severe problem in China, notably for Mobike and ofo (two popular dockless

Literature Review
The research of fleet sizing is a vital topic of BSS. Some work has provided different ideas to determine the optimal fleet size. Kőchel et al. [4] provided a joint simulation optimization idea of fleet size and bike repositioning strategy and proposed an iterative means to obtain the fleet size. Fanti et al. [5] viewed an electric car sharing system as a CQN to determine the optimal fleet size and developed an optimization problem to maximize the system revenue aiming to find the optimal fleet size. Zhai et al. [6] proposed a sparse matrix construction solution method for determining the fleet size of the DBSS. Bazan et al. [7] studied mobility-on-demand networks, which are similar to BSS. They presented a simulation model combined with optimization function and queuing network and computed the optimal fleet size. These papers provide different methods to find the optimal fleet size. However, in order to consider the stochasticity fully, we choose the CQN to model the DBSS. We determine the optimal fleet size based on an optimization function, in which the parameters are obtained by computing the CQN.
Another important measure of BSS is repositioning management. For the static repositioning, Liu et al. [8] studied a DBSS with multiple heterogeneous vehicles, depots, and visits. Their conclusion was that the enhanced chemical reaction optimization (CRO) had better solutions than the preliminary CRO. Schuijbroek et al. [9] proposed a mixed-integer programming based on decomposing the multi-vehicle repositioning matters into single-vehicle problems. They also provided a heuristic of cluster-first route-second to mitigate its running time. For other static repositioning research, readers can refer to Szeto et al. [10], Bulhões et al. [11], Jia et al. [12], and Tang et al. [13]. Static repositioning typically can not adequately capture the surges in user demand. Thus dynamic repositioning needs to be studied sufficiently. For the dynamic repositioning, Shu et al. [14] found that repositioning execution can get an additional 15-20% trips. Caggiani et al. [15] carried out dynamic bike repositioning at a constant gap time. The aim is to get a high user satisfaction and a low repositioning cost. Legros et al. [16] developed an implementation decision-support tool to decide the number of bikes and which station to carry out the repositioning. They also proved that there is an optimal inventory level. Mellou and Jaillet [17] proposed a novel mixed-integer programming formulation to solve the dynamic repositioning problem and provided a linear programming model to capture the bike flows from all trips. The authors of these studies computed the repositioning flows requiring a huge number of trucks, but in this paper we relax this requirement by setting the availability of parking regions to a certain level. We also provide a probabilistic guarantee for meeting any given service level.
CQN is an appropriate method to model BSS. George et al. [18] modeled a vehicle rental system with fixed stations as a CQN. They analyzed the asymptotic behavior of vehicle availability and provided an optimization function for determining the optimal fleet size. Li et al. [19] analyzed an SBBSS as a virtual CQN, in which parking regions and roads were considered as "virtual nodes". This was the first time that road nodes were analyzed in detail in the CQN. Li et al. [20] extended the CQN model of Li et al. [19] from three points: (i) from strong connections between stations to an irreducible graph of stations; (ii) from Poisson arrival of users to Markovian arrival processes; (iii) from exponential distributions of ridding-bike times on the roads to phase-type distributions. Furthermore, Li et al. [21] developed a fluid and diffusion approximation of a multiclass CQN to analyze the SBBSS. In this paper, we extend the computation of CQN by applying the MVA and FES algorithms. Further, some other important examples include Zhang et al. [22], Fricker et al. [23], Mizuno et al. [24], Calafiore et al. [25], Samet et al. [26], Iglesias et al. [27], and Vishkaei et al. [28].

Model Description
In this section, we give a detailed model description for the DBSS.
(1) Parking regions: In the DBSS, parking regions are viewed as nodes. Since the geographical locations and surrounding traffic environment may be different, we determine that the parking regions are different. We assume that the parking space of each parking region is sufficient so that bikes can always be parked upon arrival at a parking region immediately, that is, having no blocking effects. We assume that there are N parking regions distributed in the DBSS and denote the set of parking region nodes by S = {1, 2, . . . , N}.
(2) Roads: Parking regions are connected by roads. We express the set of roads as I. Since there are origin-direction links of parking regions, we view the roads between parking regions are directed. Road i → j and Road j → i are different roads for 1 ≤ i, j ≤ N. We view Road i → j as a node with finite servers and assume that the travel time from Parking region i to Parking region j is exponential with mean random variable 1/µ ij .
We denote the set of downlink parking regions of Parking region i by Θ i . Similarly, we denote the set of uplink parking regions of Parking region i by ∆ i .
To express the set of parking regions starting from Parking region i for 1 ≤ i ≤ N, we write The number of directed roads in the set R (i) at most is N − 1. The number of elements in the set . We assume that there is an irreducible path in the DBSS. The irreducibility is guaranteed by a suitable road construction with R (i) for 1 ≤ i ≤ N. The set of road nodes is denoted by I = {N + 1, N + 2, . . . , M}. For the detailed expression of irreducible path, readers can refer to Li et al. [20].
(3) User arrival processes: We assume that the arrivals of users at Parking region i is a Poisson process and set the mean arrival rate as λ i > 0 for 1 ≤ i ≤ N.
(4) Bike using processes: When a user arrives at Parking region i with at least one bike, he rents a bike and rides it on the Road i → j with probability p ij for j ∈ Θ i and ∑ j∈Θ i p ij = 1 for each i = 1, 2, . . . , N. The service times on the Road i → j are i.i.d. and exponential with mean rate µ ij > 0.
(5) The departure discipline: The users will leave the system with two cases: (i) when a user arrives at a parking region that has no bikes, he will leave the system directly. (b) When a user finishes his trip and returns the bike to a parking region successfully, he leaves the system immediately.
All the random variables above are independent of each other. At each parking region, bikes form a queue while waiting for users to rent. Thus, from the perspective of bikes, we consider the parking regions as single-server (SS) nodes and the roads as finite-server (FS) nodes (the number of the servers on the road is less than or equal to the total number of bikes in the DBSS). Parking region i can be viewed as a server with service rate λ i and follows First Come First Served (FCFS) scheduling for 1 ≤ i ≤ N. The service time of a parking region is the interval time of customer arrival. A bike leaves an SS node and moves to an FS node. Parking region nodes are viewed as M/M/1 queuing systems, while road nodes are modeled as M/M/n queuing systems. Figure 1 illustrates a DBSS with three parking regions.

Closed Queuing Network
All the bikes are just circulating among parking regions and roads in the DBSS, and the number of bikes in a DBSS is fixed. From bikes' perspective, we model the DBSS as a closed queuing network.
First, we let n i (t) denote the number of bikes in Node i at time t and Q (t) = (n 1 (t) , n 2 (t) , · · · , n N+M (t)). Thus {Q (t) : t ≥ 0} is a Markov process. When a user rents a bike at Parking region i, he selects a destination of Parking region j with probability p ij , where p ij > 0, p ii = 0, and ∑ j∈Θ i p ij = 1. For each i ∈ I, we define u (i) and d (i) to be the uplink and downlink nodes of i, which denote the origin and destination nodes of i, respectively. The elements of routing matrix P are written as otherwise.
The service time of each node is exponential. For i ∈ S, µ i = λ i . For Node i → j, µ ij = n ij µ ij where n ij is the number of bikes on Node i → j.
We set n i as the number of bikes in Node i. To study the CQN, the service rates, the routing matrix and relative arrival rates are needed to be determined first. For a given fleet size K, the state space of the where n = (n 1 , n 2 , . . . n N , n N+1 , n N+2 , . . . , n N+M ). The relative arrival rate e i is written as Our network consists of finite-server nodes and single-server nodes with exponential service time, so this model is a BCMP network, which has a product form solution (See Baskett et al. [29]). The product form solution is given by in which Note that n ∈ Ω and For the detailed proof of the product form solution, readers can refer to Baskett et al. [29]. For the DBSS, there are two vital performance measures: (1) Expected revenue E charged by the bikes on roads (2) Service level A i measured by the mean bike availability Expected revenue E is the utmost concern of operators. The value A i is an acknowledged measure of the quality of service. Incentive measures can be implemented based the results of A i . For example, measures include guiding customers to rent bikes at parking regions with high relative utilization and return bikes to parking regions with lower utilization.

MVA Algorithm of the Closed Queuing Network
Since the computation of normalization constant G is expensive, we apply the MVA algorithm (see Reiser and Lavenberg [30]) to compute the CQN by avoiding the calculating of the normalization constant. Through the MVA algorithm, we can calculate the mean waiting times T i (k) and the mean queue lengths K i (k) at each Node i of a CQN, where k = 1, . . . , K.
Based on the the MVA algorithm, we can get the mean response time at the ith node Based on T i (k), the overall throughput is given by and the throughput of each node can be written as in which e i is determined by Equation (1). Based on T i (k) and λ i (k), the mean number of bikes at the ith node is written as We provide the detailed computation by MVA in Algorithm 1.

Algorithm 1
Mean value analysis (MVA) algorithm for closed queuing network.

Step 2. Computation for the relative arrival rates
Let e 1 = 1. Compute relative arrival rates e j for j = 2, 3, . . . , N + M, based on equation .
Get the normalization constant based on equation G (k) = G(k−1) λ(k) . Step 4.3 Compute the mean number of bikes at each node based on equation K i (k) = e i λ (k) T i (k).
Step 5. Computation for performance measures E i and A i in the closed queuing network.
We are interested in the availability of bikes at each parking region (i.e., the probability that there is at least one bike at the parking region). Set T i (0) = K i (0) = 0. By iteration, we can get For small networks, this algorithm provides a method to obtain the optimal fleet size straightly. However, for a larger network, the computation is too expensive. Thus, we need to find another approach to reduce the computation.

FES Algorithm of the Closed Queuing Network
MVA algorithm can be computed directly without calculating the normalizing constant, but it needs high storage requirement and expensive computation. So we try to find another approach to reduce the computation. We use the FES algorithm to compute the DBSS.
In the closed queuing network of DBSS, since bikes are served in parking regions as well as on roads, we divide the nodes into SS nodes and FS nodes. However, when we consider the two performance measures: expected revenue E i and available rate A i , we only consider the results of SS nodes. Since the difference between parking region nodes and road nodes, we try to simplify the computation by partitioning them.
FES method is based on Norton's theorem, refer to Chandy et al. [31] and Akyildiz et al. [32]. Based on the FES method, we partition the closed queuing network into N disjoint subnetworks, in which Node i and node i → j are a subnetwork for i = 1, . . . , N, j ∈ Θ i . We write the subnetwork as SN 1 , SN 2 , . . . , SN N . Analyzing one of the subnetworks, we need to short-circuit all the other nodes that do not belong to the subnetwork.
First, we need to determine the throughput λ SN j (k) and the normalization constant G j (k) for 1 ≤ j ≤ N and 1 ≤ k ≤ K. The arrival rates of nodes in the subnetwork are the same as the original network. Consider each subnetwork as an FES node and construct an equivalent reduced network. The dependent service rates µ ej of the FES Node j are identical with the throughputs in the corresponding jth subnetwork, which is written as µ ej (k) = λ SN j (k) , for j = 1, . . . , N and k = 1, . . . , K.
We provide the detailed FES in Algorithm 2. In Algorithm 2, the time and storage requirement are reduced considerably. Algorithm 2 Flow equivalent server (FES) algorithm for closed queuing network.

Step 1. Partition the network into N subnetworks.
Partition the network into subnetworks SN 1 , SN 2 , . . . , SN N , in which Node i and Node i → j constitute SN i subnetwork for i = 1, . . . , N, j ∈ Θ i .
Compute the normalization constant of the reduced network by convoluting the N normalizing vectors, based on the equation

Performance Analysis
In this section, we provide the results of a performance analysis of DBSSs. Based on the MVA and FES algorithms, we get the results of the expected number of bikes and the availability for each parking region. Then, we propose two optimization functions to determine the optimal fleet size and the repositioning flow based on the results.

Optimal Fleet Size
Since the optimal fleet size is a tradeoff between the earnings got from users and maintaining cost of the fleet, we present an optimization function to determine the optimal fleet size. Let r i be the revenue obtained from a bike being in use at road nodes per-unit-time. The BSS will cause an unavailability cost c j whenever bikes park at parking regions idly. Then the objective is to maximize the profit in which E i is the expected queue length at Node i for i ∈ I and A j is bike availability at Node j for j ∈ S in steady state. Based on Shanthikumar and Yao [33] and George et al. [18], we get that for a single-class CQN with FCFS exponential servers, both the throughput at each parking region and the system throughput are non-decreasing concave functions of fleet size. The values {A j (K), j ∈ S} and {E i (K), i ∈ I} are all non-decreasing concave functions, therefore F (K) is concave.
To solve Equation (7), we can calculate the MVA algorithm iteratively for each K ∈ {1, 2, ...} and compute the associated profit. Since the concavity of F (K), the method will terminate once profit decreases from K to K + 1.

Repositioning Flow
Since bikes will inevitably accumulate at one or more parking regions in the DBSS. To ease this problem, we introduce a set of repositioning flow to rebalance the system, i.e., moving bikes to parking regions that lack bikes. We use y ij to denote the rate of bikes repositioning flow from Parking region i to Parking region j. We discuss how to achieve a target service level with different bike resources.
We assume that Y i = ∑ j∈Θ i y ij . The availability of Parking region i is given by If we set the availability A i of Parking region i equal to 1, it needs a huge number of bikes idly waiting at each parking region to meet the availability (refer to Li et al. [34]). Thus, we set the objective as a service level min{A i } ≥ σ ∈ (0, 1). Therefore, we can realize the certain level of availability and release stressful demand of the fleet size. Then the repositioning cost minimization problem can be written by where b j is the cost of maintaining bikes and b j = ∑ j∈Θ i b ij p ij .
We set ξ i = e i /λ i , ξ i = max {e i /λ i }, γ i = e i / (λ i + X i ), and γ i = max {e i / (λ i + X i )} for i ∈ S. Since A i = γ i /γ i , the optimal equation can be rewritten as The optimal solution is

Numerical Experiments
For the closed queuing network, the MVA algorithm can be used to compute a small size DBSS and the FES algorithm can be used to calculate larger size DBSS. We used a small size system to verify the computational accuracy of the MVA and FES algorithms. Then we provide the results of the computation of optimal fleet size and repositioning flow. For all the numerical experiments, we implemented the computation in Matlab on a laptop.

Comparison of Algorithm 1 and Algorithm 2
We studied a small size DBSS with three parking regions and the three parking regions are connected to each other and so there are six road nodes. We describe the DBSS as a CQN and set the fleet size as 45. The physical structure of the DBSS is depicted in Figure 2. The routing matrix P is given by Based on the routing matrix P and Equation (1), we can obtain relative arrival rates. We list the service rates and the relative arrival rates in Table 1. Using the FES algorithm, we divided the nine nodes into three subnetworks. i.e., Node 1, Node 4 and Node 6 are subnetwork 1; Node 2, Node 5, and Node 8 are subnetwork 2; Node 3, Node 7, and Node 9 are subnetwork 3. We computed the DBSS by Algorithm 1 and 2 respectively and verified that the results are consistent. We provide the results in Table 2.

Optimal Fleet Size and Repositioning Flow
(1) The three parking region DBSS Let the hourly rental fee r i for i ∈ I be $2 and the maintenance cost c i per-bike per-hour for i ∈ S be $0.20. By Algorithm 2 and Equation (7), we compute the three parking region DBSS depicted in Figure 2. The variation of profits corresponding to fleet size is shown in Figure 3. We can find that the optimal fleet size is 60. This method considers the differences between parking regions and loose structures of the routing matrix. We can find that the profit first increases and then decreases as the fleet size K increases. This result offers operators quantitative basis of the fleet size design. When setting min i∈S {A i } ≥ 0.9, we get the optimal repositioning strategy to add rebalancing flows 1 → 4 with a rate of 0.0641, 1 → 6 with a rate of 0.0547, and 3 → 8 with a rate of 0.1431. The min i∈S {A i } ≥ 0.9. The computation of repositioning flow considers the dynamics of DBSS, which is close to the actual operation.
(2) A 50 parking region DBSS By Algorithm 2, we computed a DBSS with 50 parking regions. We partitioned the DBSS into 50 subnetworks and depict the optimal fleet size in Figure 4. The optimal repositioning flows are listed in Table 3, in which y i,j is the flow from Parking region i to Road i → j. Algorithm 2 extends the computation of CQN of DBSS greatly. Since there are actual service processes on the road in the DBSS, when we analyze the system we must study the service rates on the road nodes. It is well known that the number of road nodes is enormous in a DBSS, which will cause the state space explosion and lead to the computation of normalization constant too expensive. By partitioning the system into subnetworks, the number of nodes is reduced significantly, which realizes the performance computation of larger size DBSS. Thus the result provides operators with valuable quantitative reference.

Conclusions
In this paper, we provide a closed queuing network of dockless bike sharing systems. To avoid the expensive computation of the normalization constant, we use the MVA algorithm and the FES algorithm to analyze the small size and larger size DBSS, respectively. The FES algorithm is proposed to compute the DBSS for the first time. We let a parking region and its downlink roads as a subnetwork and partition the DBSS into different subnetworks. The FES algorithm reduces the number of nodes greatly, and thus the computation of the CQN is lessened significantly. This paper opens a new avenue to study the CQN of DBSS. Based on the computation results of the CQN, we determine the optimal fleet size and repositioning flow by two optimization functions. The computation of repositioning flow adequately considers the dynamics of bikes during the repositioning process, which is close to the actual DBSS. According to our analysis, the method can be used as guidelines for designing and controlling the fleet and arrange the repositioning flow to achieve a high level of service.
For future work, the CQN model presented in this paper can be used to analyze other on-demand transportation networks. Another upcoming issue is that we can consider the unusable bikes in the DBSS. Since the failure rates of bikes are high, unusable bikes account for a large proportion in the DBSS. Therefore when determining the optimal fleet and the repositioning flow, the failure rate of bikes can be introduced in.

Conflicts of Interest:
The authors declare no conflict of interest.