A Stochastic Analysis of Bike Sharing Systems
Shuang Tao Jamol Pender





      

A Stochastic Analysis of Bike Sharing Systems
     
  Shuang Tao 
 School of Operations Research and Information Engineering 
 Cornell University

 293 Rhodes Hall, Ithaca, NY 14853 
  st754@cornell.edu  
 
Jamol Pender 
 School of Operations Research and Information Engineering 
 Cornell University

 228 Rhodes Hall, Ithaca, NY 14853 
  jjp274@cornell.edu  
 
 
    December 30, 2023
==============================================================================================================================================================================================================================================================================================================================


 As more people move back into densely populated cities, bike sharing is emerging as an important mode of urban mobility.  In a typical bike sharing system, riders arrive at a station and take a bike if it is available. After retrieving a bike, they ride it for a while, then return it to a station near their final destinations.  Since space is limited in cities, each station has a finite capacity of docks, which cannot hold more bikes than its capacity.  In this paper, we study bike sharing systems with stations having a finite capacity.  By an appropriate scaling of our stochastic model, we prove a mean field limit and a central limit theorem for an empirical process of the number of stations with k bikes.  The mean field limit and the central limit theorem provides insight on the mean, variance, and sample path dynamics of large scale bike sharing systems.  We also leverage our results to estimate confidence intervals for various performance measures such as the proportion of empty stations, the proportion of full stations, and the number of bikes in circulation.  These performance measures have the potential to inform the operations and design of future bike sharing systems.      








§ INTRODUCTION


Bike sharing  is an emerging mode of eco friendly transportation that have launched in over 400 cities around the world (<cit.>).  In the United States, we are witnessing a transition where more people are deciding to live in large and densely populated cities.  As more people transition from the sprawling suburbs into densely populated cities, bike sharing programs will continue to grow in popularity since they provide easy transportation for citizens of these large cities.  As more people use these bike sharing systems(BSS), less people will drive motor vehicles on the road.  This reduction of vehicles on the road due to BSS has the potential to also reduce the growing traffic congestion in these growing cities.  BSS also promote healthy living as biking is a great form of exercise. They are environmentally friendly and they have the potential to reduce carbon emissions if operated correctly and efficiently, see for example <cit.> and their references within for more information on BSS.

For a typical system, riders simply arrive at a station and select a bike if there is one available for them to take.  If there is no bike available for the rider, the rider will leave the system.  Otherwise, the rider will take a bike and ride it for a while before returning it to another station near their final destination, if there is available space. If no space is available, the rider must find a nearby station to return the bike.  If there were an infinite supply of docks to store the bikes, then our model would be reduced to a network of infinite server queues, which is more tractable to analyze.  However, since the number of bike docks have finite capacity, the model become less tractable especially for systems with a large number of stations.  

Much of the complexity inherent in BSS lies in the scarcity of resources to move all riders around each city at all times of the day. Riders can encounter the scarcity of resources in two fundamental ways.  First, a rider can encounter insufficient resources by finding an empty station with no bikes when a rider needs one.  Secondly, a rider can find a station full with bikes when attempting to return a bike.  Thus, from a managerial point of view, having stations with no bikes or too many bikes are both problematic for riders.  There are several reasons why bike stations either have no or too many bikes.  One main reason why stations might have too many bikes or too few bikes is that the system is highly inhomogeneous.  Not only is the arrival rate a non-constant function of the time of day (see Figure <ref>), but also riders do not evenly distribute themselves amongst the available stations.  For example, many stations that are located in residential areas have fewer bikes available for riding as many riders take bikes to more commercial areas.  Another example that illustrates the inhomogeneous dynamics is that riders tend to take bikes from up-hill stations to go to down-hill stations, however, very few riders take bikes from down-hill stations to go up-hill.  Thus, as more bikes flow from residential to commercial areas during rush hours, or up-hill to down-hill stations, this causes the system to be more imbalanced over time. 

Since BSS are quite complex, researchers have been inspired to study these systems in great depth.  The subsequent analysis of BSS has generated many insights on these systems, especially for rebalancing the fleet of bikes.  Although there is a large community that studies these systems, few analytical models have been proposed, especially stochastic analytical models.  This is primarily because the stochastic models for BSS are often very complex and are rather intractable to analyze without making strong assumptions.  Nevertheless, the analysis of such stochastic and mathematical models could provide insights on the behavior of these BSS and how to manage them effectively. In fact, a deeper analysis of  stochastic models for such systems could help researchers understand the impact of different incentive algorithms for taking or returning bikes, which can lead to significant improvements in the overall system performance.   

Most BSS can be viewed as closed queueing networks. One of the first papers to model the bike sharing system as a closed queueing network is by  <cit.>. However, in the bike sharing context, the customers are replaced by bikes. The number of bikes, also known as the fleet size, is fixed and remains constant.  In this model, the bikes can go to two types of stations.  The first type of station is a single-server queue where the service times are the user inter-arrival times to this station. The second type of station is an infinite-server queue where the service times are the trip times on a route from station i to station j. The main drawback of the model by <cit.> is that it is based on infinite capacity queues.  This means that the model does not take into account the finite capacity of the stations and the related strategies of the users to return their bikes.  To overcome this major drawback, the model proposed in this paper allow the finite capacity at stations.  Thus, we are able to model the real system where customers are blocked from returning bikes to the stations that is nearest to their destination.  

In our model, we model the bike sharing system as state dependent M/M/1/K_i queueing networks.  When joining a saturated single-server queue, the user reattempts in another queue, after a time with the same distribution as the trip time, until he returns his bike.  Although this model seems to model the behavior of the network, it is not practical since it scales with the number of stations.  Thus, we follow an approach developed by <cit.> to study the bike sharing network's empirical process instead.  The empirical process still allows us to derive important performance measures of the original system, however, it scales with the maximum station capacity and not the number of stations, which is more practical for large networks like CitiBike.  







 

 §.§ Main Contributions of Paper


The contributions of this work can be summarized as follows:    


  * We construct a stochastic bike sharing queueing model that incorporates the finite capacity of stations. Since our model is difficult to analyze for a large number of stations, we propose to analyze an empirical process that describes the proportion of stations that have a certain number of bikes. 

  * We prove a mean field limit and a central limit theorem for our stochastic bike sharing empirical process, showing that the mean field limit and the variance of the empirical process can be described by a system of 1/2(K+4)(K+1) differential equations where K is the maximum station capacity.    

  * Using the mean field limit of the empirical process, we are able to approximate the mean proportion of empty and full stations.  Furthermore, with the central limit theorem of the empirical process, we are able to construct confidence intervals around the mean field limit for the same performance measures.

  * We compare the mean field limit and the central limit theorem to a simulation and show that the differential equations approximate the mean and variance of the empirical process extremely well. 
 






 §.§ Organization of Paper


The remainder of this paper is organized as follows. Section <ref> provides a brief history of bike sharing programs and a literature review on the research streams concerning BSS. Section <ref> introduces our bike sharing model and notation of the paper.  In Section <ref>, we derive amd prove the mean field limit of the empirical measure process of the distribution of stations with different bikes and utilization rate. In Section <ref>, we derive the diffusion limit and prove a functional CLT for our model. We show that the diffusion process is a centered Gaussian OU process and we also obtain a closed form expression of the diffusion limit process. In Section <ref>, we extend our analysis in Section <ref> to a broader case with non-uniform routing probabilities and capacities, and derive the mean field limit and diffusion limit in this extended case. In Section <ref>, we discuss the simulation results of our model, with both stationary and non-stationary arrival processes. We also give a comparison between real CitiBike empirical measure and simulated ones using our model to show how well our model is in capturing reality.  Finally, in Section <ref>, we give concluding remarks and provide some future directions of research that we intend to pursue later.





§ HISTORY AND LITERATURE REVIEW


The literature that focuses on the analysis and operations of BSS is increasing rapidly. Early research that studied the history of BSS includes <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>.  These papers provide a history of bike sharing and how it has evolved over time.  The beginning of bike sharing can be traced back to the first generation of white bikes (or free bikes) in Amsterdam, The Netherlands as early as 1965. However, this first generation of BSS failed due to a large amount of bike theft.  The launch of Bycyklen in Copenhagen in 1995 marked the second generation of BSS.  Bycyklen was the first BSS to implement docking stations and coin-deposit systems to unlock bikes.  These coin-deposit systems helped with the bike theft problem and thus made the BSS more reliable. However, even with these improvements, the Bycyklen could not eliminate bike thefts mainly due to the fact that customer still remained anonymous and there were no time limits on how long a customer could use a bike.  

The failures of the second generation of BSS inspired the present-day or the third generation of BSS by combining docks with information technology. These new systems incorporate information technology for bicycle reservations, pickups, and drop-offs.  This new technology has enabled many BSS to keep better track of bicycles and the users that use them, thus eliminating virtually all bike theft.  One example of a third generation BSS is the Paris’ bike sharing program Velib.  Velib was launched in Paris in July 2007 and has emerged as the most prominent example of a successful bike sharing program in the modern world. As a result of the success of Velib,  many cities like New York City (Citi Bike), Chicago (Divvy), and even Ithaca (Big Red Bikes) have implemented large-scale BSS and bike sharing has become a widely used form of transportation in these cities.  For the interested reader,  <cit.> provides a comprehensive survey of the vehicle/bike sharing literature.

Rebalancing is currently the biggest stream of research concerning BSS.  In rebalancing operations, there are two methods of rebalancing: (1) deploying a truck fleet or (2) providing user incentives, where deploying a truck fleet is often referred to as bike repositioning. Both methods involve static and dynamic cases.  Static repositioning usually is moving bikes during the night when traffic flow is low, while dynamic repositioning is moving bikes during the day based on current state of the system. Most research on this area focuses on the static case, partly because it is easier to model and also because the impact of repositioning is more important during the night (<cit.>). <cit.> is one of the first papers to study static repositioning of BSS, using mixed integer linear programming by maximizing customer demand satisfaction. <cit.> consider a similar problem, where a single truck repositions bikes to bring the inventory of each station to a predetermined value.  However, their objective is to minimize the routing cost as opposed to maximizing customer satisfaction. In the case of dynamic repositioning, <cit.> and <cit.> consider the case when the trucks respond in real time to the current state of the system.  However, <cit.> and <cit.> consider the situation where the time-dependent demand is known a priori and the rebalancing operations are computed in an off-line fashion. Yet, none of these papers really explore stochastic dynamics and they for the most part mainly exploit optimization techniques to tackle the problem.

Unlike the rebalancing literature, our paper falls into the performance analysis literature with an emphasis in supply analysis. We focus on analyzing the most salient performance measures such as the mean, variance, covariance and sample path dynamics of the bike distributions in a large-scale BSS.  There is not much literature that explores the fluctuations of BSS around the mean field limit.  In this paper, we prove a mean field limit and a functional central limit theorem under some smoothness conditions and show that the diffusion limit is characterized by a multi-dimensional Ornstein-Uhlenbeck (OU) process.  The functional central limit theorem not only gives us information about the sample path fluctuations of the queue length process, but it also allows us to construct approximate confidence intervals for various performance measures such as the proportion of empty stations, the proportion of full stations, and the mean number of bikes in circulation. Unlike the previous literature, our paper also considers non-stationary arrivals to stations, which is much more realistic given the user patterns we observe in the historical data from Citi Bike. In Figure <ref>, we plot the empirical mean of the number of trips (5 minute intervals) during the week for the time period Jan 1st-Dec 31st, 2015. We observe from Figure <ref> that the arrival rate is non-stationary and clearly reflects the morning rush and evening rush during the peak times. Thus, analyzing the non-stationary dynamics is crucial for understanding the impact of system inhomogeneity since it can provide useful guidelines for rebalancing operations.








§ BIKE-SHARING QUEUEING MODEL
 



In this section, we construct a Markovian bike sharing queueing model where customers can pick-up and drop-off bikes at each station if there is available capacity.  Figure <ref> provides an illustration of a typical Citi Bike station in New York City (NYC).  As one can see in Figure <ref>, the bikes are attached to docks and the number of docks is finite with roughly 40 docks.  Figure <ref> shows a map of Citi Bike stations, the nation's largest bike sharing program, with over 10,000 bikes and 600 stations across Manhattan, Brooklyn, Queens and Jersey City. Citi Bike was designed for quick, affordable and convenient trips, and has become an essential part of the transportation infrastructure in NYC.



 
Motivated by the Citi Bike bike sharing system, we consider a bikes sharing system with N stations and a fleet of M bikes in total. We assume that the arrival of customers to the stations are independent Poisson processes with rate λ_i for station i.  When a customer arrives at a station, if there is no bikes available, they will then leave the system and are immediately blocked and lost. Otherwise, the customer will take a bike and ride to station j with probability P_j.  We assume that the travel time of the rider is exponentially distributed with mean 1/μ, for every transition from one station to another.  Since we are concerned with finite capacity stations, we assume that station i has a bike capacity of K_i, which is assumed to be finite for all stations.   Thus, when a customer arrives at station j, if there are less than K_j bikes in this station, he returns his bike and leaves the system. If there are exactly K_j bikes (i.e. the station is full), the customer randomly chooses another station k with probability P_k and goes to that station to drop the bike off.  As before, it takes a time that is exponentially distributed with mean 1/μ. Finally, the customer rides like this again until he can return his bike to a station that is not full. 

Below, we provide Table <ref> for the reader's convenience so that they understand the notation that we will use throughout the paper.  




To avoid cumbersome notation, throughout Sections <ref>, <ref> and <ref>, we assume without loss of generality that the service rate μ is equal to 1 and that the routing probability from each station is uniform, i.e. P_i=1/N. We also assume that capacities across all stations are equal, i.e. K_i=K, for  i=1,...,N. With our notation in hand, we are ready to develop our stochastic model for our bike sharing network.  At first glance, these assumptions seem restrictive, however, we explain in Section <ref>  how our model can be extended seamlessly to more complex settings with non-uniform routing probabilities and capacities.  We should mention, however, the extension to inter-station transition probabilities i.e probabilities that depends on the departing station as well as the returning station are non-Markovian models.  In Section <ref> we discuss more in detail how one can extend the state space, by tracking the number of in-transit bikes coming from each station, to make the queueing model Markovian, however this extended model makes the problem even more high-dimensional and adds difficulty to the analysis.


We define 𝐗(t)=(X_1(t),⋯, X_N(t)), where X_i(t) is the number of bikes at station i at time t. Then X_i(t) is a continuous time Markov chain(CTMC), in particular a state dependent M/M/1/K_i queue. In this model, the rate of dropping off bikes at station i is equal to μ P_i(M-∑_k=1^NX_k(t))1{X_i(t)<K_i}, and the rate of retrieving bikes at station i is equal to λ_i1{X_i(t)>0}.  Using these rates, we can construct the functional forward equations for our stochastic bike sharing model.  This construction is given below in the following proposition:

 
 For any integrable  function f: ℤ_+^N→ℝ, the CTMC 𝐗(t) satisfy the following functional forward equation,
 
    𝔼[f(𝐗(t)) | 𝐗(0) = 𝐱]  ≡d/dt𝔼[f(𝐗(t)) | 𝐗(0) = 𝐱]  
       =   ∑_i=1^N𝔼[(f(𝐗(t)-1_i)-f(𝐗(t))λ_i1_{X_i(t)>0}]
           +∑_i=1^N𝔼[(f(𝐗(t)+1_i)-f(𝐗(t))μ P_i(M-∑_k=1^NX_k(t))1_{X_i(t)<K_i}] .
 
 
 
 

The proof is found in the Appendix.  

 

The time derivatives of the mean, variance, and covariance of 𝐗(t) are given by

    𝔼[X_i(t)]      =   𝔼[ μ P_i(M-∑_k=1^NX_k(t))1{X_i(t)<K_i}]  - λ_iℙ[ X_i(t) > 0   ] ,
     
    Var[X_i(t)]     =   𝔼[X_i^2(t)]-2𝔼[X_i(t)]𝔼[X_i(t)]
       =    2 Cov[X_i(t),μ P_i(M-∑_k=1^NX_k(t))1{X_i(t)<K_i}]  
       -    2 Cov[X_i(t),  λ_i1{ X_i(t) > 0 }]
        +   𝔼[ μ P_i(M-∑_k=1^NX_k(t))1{X_i(t)<K_i}] +  λ_iℙ[ X_i(t) > 0   ]  ,



    Cov[X_i(t), X_j(t)] 
       =   Cov[μ P_j(M-∑_k=1^NX_k(t))1_{X_j(t)<K_j}-λ_j1_{X_j(t)>0},X_i(t) ]
       +    Cov[μ P_i(M-∑_k=1^NX_k(t))1_{X_i(t)<K_i}-λ_i1_{X_i(t)>0},X_j(t) ],

for i,j=1,⋯,N and i≠ j.













  
   

 §.§ Intractability of Individual Stations Model

   
Although the functional forward equations given in Corollary <ref> describe the exact dynamics of the mean, variance and covariance of the bike sharing system, the system of differential equations are not closed.  This non-closure property of the functional forward equations in this model arises from the fact that the bike sharing system has finite capacity.  More importantly, it also implies that we need to know a priori the full distribution of the whole stochastic process 𝐗(t) in order to calculate the mean or variance or any moment for that matter. Work by <cit.> could yield useful and accurate closure approximations for making the system closed.  Moreover, withe exception of <cit.>, there  are no error bounds on the accuracy of various closure approximations.  Thus, it is not clear how well the closure approximations would perform over a variety of parameter settings.  Finally if we even wanted to solve these equations and knew the entire distribution of 𝐗(t) a priori, there are still O(N^2) differential equations(around 180,900 equations in the CitiBike case) that would need to be numerically integrated.  This is very computationally expensive and thus, we must take a different approach to analyze our bike sharing system. 
 
Moreover, if we want to analyze the limiting behavior of {X_i}_i=1^N as  CTMCs, as we let N go to ∞, the mean field limit would become infinite dimensional, which is quite complicated.  However, if we instead analyze an empirical measure process for 𝐗(t), we can use the finite capacity nature of the bike sharing system to our advantage and have a finite dimensional CTMC for the empircal measure process.  

   

 §.§ An Empirical Measure Model

 
Following the model of <cit.>, we construct an empirical measure process that counts the proportions of stations with n bikes and utilization r.  This empirical measure process is given below by the following equation: 

    Y_t^N(r,n)=1/N∑_i=1^N1{r_i^N=r,X_i^N(t)=n}.


We further define that Y_t^N(n)=∑_rY_t^N(r,n). By observing the empirical measure process, we notice that Y_t^N=(Y_t^N(0),⋯,Y_t^N(K)) ∈ [0,1]^K+1.  Thus, for our empirical measure process, we only need to solve O(K^2) differential equations for understanding the mean and variance dynamics of the bike sharing system, where K<<N. More importantly, the empirical measure will also allow us to obtain salient performance measures such as Y_t^N(0) (the proportion of stations with no bikes), Y_t^N(K) (the proportion of stations that are full of bikes), M - ∑^K_j=0 j · Y_t^N(j)N (the number of bikes in circulation), among others.

Conditioning on Y_t^N(r,n)=y(r,n) and given our assumptions that μ=1 and P_i=1/N, the transition rates of y are specified as follows:

When a customer arrives to a station with n bikes and relative utilization r to retrieve a bike, the
proportion of stations having n bikes goes down by 1/N, the proportion of stations having n-1 bikes goes up by 1/N, and the transition rate Q^N is 

    Q^N(y,y+1/N(1_(r,n-1)-1_(r,n)) )    =    
    y(r,n)λ_r N 1_n>0
       =    y(r,n)μ P_i/RN 1_n>0
       =    y(r,n)1/NR N1_n>0
       =   y(r,n)/rR^N_max1_n>0.

When a customer returns a bike to a station with n bikes and relative utilization r, the proportion of stations having n bikes goes down by 1/N, the proportion of stations having n+1 bikes goes up by 1/N, and the transition rate Q^N is 

    Q^N(y,y+1/N(1_(r,n+1)-1_(r,n)) )    =   
     y(r,n) ·μ·(M-∑_n'∑_r'n'y(r',n')N)1_n<K
       =   
     y(r,n)N(M/N-∑_n'∑_r'n'y(r',n'))1_n<K.

Similarly, we have the functional forward equations for Y_t^N(r,n).


 For any integrable function f: [0,1]^K+1→ℝ, Y_t^N(r)=(Y_t^N(r,0),⋯,Y_t^N(r,K)) satisfies the following functional forward equation,
 
    𝔼(f(Y_t^N(r)) | Y_0^N(r) = y_0(r)]
       =   ∑_n=0^K𝔼[(f(Y_t^N(r)+1/N(1_r,n-1-1_r,n))-f(Y_t^N(r)))Y^N_t(r,n)/rR_max1_n>0]
       +   ∑_n=0^K𝔼[(f(Y_t^N(r)+1/N(1_r,n+1-1_r,n))-f(Y_t^N(r)))Y^N_t(r,n)N(M/N..
       -   ..∑_n'∑_r'n'Y_t^N(r',n'))1_n<K]

 
 
  
  The proof is similar to the proof of Proposition <ref>.
  
 


















 
 
  The time derivative of the mean of Y_t^N(r,n) is given by 
 
    𝔼[Y_t^N(r,n)]     =   𝔼[ 1/rNR^N_max(Y^N_t(r,n+1)1_n<K-Y^N_t(r,n)1_n>0) ] 
           + 𝔼[ (M/N-∑_n'∑_r'n'Y_t^N(r',n'))  (Y^N_t(r,n-1)1_n>0-Y^N_t(r,n)1_n<K) ].

 for n=0,⋯,K. 
 Denote Σ_i,j^N(r,t)=Cov[Y_t^N(r,i),Y_t^N(r,j)]. When i=j, 
the time derivative of the variance term Var[Y_t^N(r,i)] is given by 
   
    Var[Y_t^N(r,i)]     =   𝔼[Y_t^N(r,i)^2]-2𝔼[Y_t^N(r,i)]𝔼[Y_t^N(r,i)]
       =   2/rNR^N_max(Σ_i,i+1^N(r,t)1_i<K-Σ_i,i^N(r,t)1_i>0)+2M/N(Σ^N_i,i-1(r,t)1_i>0 -Σ^N_i,i(r,t)1_i<K)
          - 2 ∑_n'∑_r'n'(Cov[Y_t^N(r,i),Y_t^N(r',n')Y^N_t(r,i-1)]1_i>0.
          .-Cov[Y_t^N(r,i),Y_t^N(r',n')Y^N_t(r,i) ]1_i<K)
          + 1/rN^2R^N_max(𝔼[Y^N_t(r,i+1) ]1_i<K+𝔼[Y^N_t(r,i)]1_i>0)
          +M/N^2(𝔼[Y^N_t(r,i-1)]1_i>0+𝔼[Y^N_t(r,i)]1_i<K)
           -1/N∑_n'∑_r'n'(𝔼[Y_t^N(r',n')Y^N_t(r,i-1)]1_i>0+𝔼[Y_t^N(r',n')Y^N_t(r,i) ]1_i<K).

  When |i-j|>1, the time derivative of the covariance term Cov[Y_t^N(r,i), Y_t^N(r,j)] is given by
   
    Cov[Y_t^N(r,i), Y_t^N(r,j)]  
       =   𝔼[Y_t^N(r,i)Y_t^N(r,j)]-𝔼[Y_t^N(r,i)]𝔼[Y_t^N(r,j)] -𝔼[Y_t^N(r,j)]𝔼[Y_t^N(r,i)]
       =   1/rNR^N_max[Σ_i+1,j^N(r,t)1_i<K+Σ_i,j+1^N(r,t)1_j<K-Σ_i,j^N(r,t)(1_j>0+1_i>0)]
          +M/N[Σ^N_i-1,j(r,t)1_i>0+Σ^N_i,j-1(r,t)1_j>0 -Σ^N_i,j(r,t)(1_i<K+1_j<K)]
          -  ∑_n'∑_r'n'(Cov[Y_t^N(r,j),Y_t^N(r',n')Y^N_t(r,i-1)]1_i>0-Cov[Y_t^N(r,j),Y_t^N(r',n')Y^N_t(r,i) ]1_i<K)
          -  ∑_n'∑_r'n'(Cov[Y_t^N(r,i),Y_t^N(r',n')Y^N_t(r,j-1)]1_j>0-Cov[Y_t^N(r,i),Y_t^N(r',n')Y^N_t(r,j) ]1_j<K).

and when j=i+1, the time derivative of the covariance term Cov[Y_t^N(r,i), Y_t^N(r,i+1)] is given by

    Cov[Y_t^N(r,i), Y_t^N(r,i+1)]  
       =   𝔼[Y_t^N(r,i)Y_t^N(r,i+1)]-𝔼[Y_t^N(r,i)]𝔼[Y_t^N(r,i+1)] -𝔼[Y_t^N(r,i+1)]𝔼[Y_t^N(r,i)]
       =   1/rNR^N_max[Σ_i+1,i+1^N(r,t)+Σ_i,i+2^N(r,t)1_i<K-1-Σ_i,i+1^N(r,t)(1+1_i>0)]
          +M/N[Σ^N_i-1,i+1(r,t)1_i>0+Σ^N_i,i(r,t) -Σ^N_i,i+1(r,t)(1+1_i<K-1)]
          -  ∑_n'∑_r'n'(Cov[Y_t^N(r,i+1),Y_t^N(r',n')Y^N_t(r,i-1)]1_i>0.
          .-Cov[Y_t^N(r,i+1),Y_t^N(r',n')Y^N_t(r,i) ])
          -  ∑_n'∑_r'n'(Cov[Y_t^N(r,i),Y_t^N(r',n')Y^N_t(r,i)]-Cov[Y_t^N(r,i),Y_t^N(r',n')Y^N_t(r,i+1) ]1_i<K-1)
          -𝔼[Y_t^N(r,i+1)]/rN^2R_max  -M/N^2𝔼[Y_t^N(r,i)]+1/N∑_n'∑_r'n'𝔼[Y_t^N(r',n')Y_t^N(r,i)] .


 
 
The time derivatives of 𝔼[Y_t^N(r,i)] and Var[Y_t^N(r,i)] come directly from applying Proposition <ref> with f(Y_t(r))=Y_t(r,i), Y_t^2(r,i) respectively. For the covariance term, we let f(Y_t^N(r))=f(Y^N_t(r,i),Y^N_t(r,j)) to prove the result.


























 
 


Although we have equations for the moments of the empirical process and individual stations, it is still difficult to analyze them and gain insights from them directly.  One reason is that even though we have reduced the dimensionality of the analysis significantly, we have not removed the non-closure property of the differential equations.  Thus, we need to develop a new approach that will allow us to get around this complication.  The method that we choose to use is asymptotic analysis and will be described in more details in the sequel.  However, before we get to the asymptotic analysis we state some technicalities about weak convergence.

 

 §.§ Preliminaries of Weak Convergence


Following <cit.>, we assume that all random variables in this paper are defined on a common probability space (Ω, ℱ, ℙ).  Moreover, for all positive integers k, we let 𝒟([0 , ∞), ℝ^k) be the space of right continuous functions with left limits (RCLL) in ℝ^k that have a time domain in [0, ∞).  As is usual, we endow the space 𝒟([0 , ∞), ℝ^k) with the usual Skorokhod J_1 topology, and let M^k be defined as the Borel σ-algebra associated with the J_1 topology. We also assume that all stochastic processes are measurable functions from our common probability space (Ω, ℱ, ℙ) into (𝒟([0 , ∞), ℝ^k), M^k).   Thus, if {ζ}^∞_n=1 is a sequence of stochastic processes, then the notation ζ^n →ζ implies that the probability measures that are induced by the ζ^n's on the space (𝒟([0 , ∞), ℝ^k), M^k) converge weakly to the probability measure on the space (𝒟([0 , ∞), ℝ^k), M^k) induced by ζ.  For any x ∈ (𝒟([0 , ∞), ℝ^k), M^k)  and any T > 0, we define 

    ||x||_T ≡sup_0 ≤ t ≤ T max_i = 1,2,...,k |x_i(t)|

and note that ζ^n converges almost surely to a continuous limit process ζ in the J_1 topology if and only if 

    ||ζ^n - ζ||_T → 0    a.s.

for every T > 0. 






§ MEAN FIELD LIMIT
 

In this section, we prove the mean field limit for our bike sharing model.  A mean field limit describes the large station dynamics of the bike sharing network over time.  Deriving the mean field limit allows us to obtain new insights on average system dynamics, when the demand for bikes and the number of stations are very large. Thus, we avoid the need to study an N-dimensional CTMC and compute its steady state distribution in this high dimensional setting. 

First, we state the important assumptions that will be used through out the paper, to ensure the existence of a mean field limit of our model.


 There exists a probability measure I(r,k)  on ]0,1] ×ℕ with finite support and Λ > 0 such that, as N tends to infinity, we have


  i) 1/N∑_i=1^N1 _(r_i^N,K_i^N)⇒ I(r,k),

  ii)NR_max^N→Λ^-1,

  iii)M/N→γ.





Now we state the main theorem in this section that proves the convergence of empirical process to its mean field limit.

Let |· | denote the Euclidean norm in ℝ^K+1. Under Assumption <ref>, suppose that Y_0^N y_0,  then we have for any ϵ>0 and t_0>0,

    lim_N→∞P(sup_t≤ t_0|Y_t^N-y_t|>ϵ)=0.

Here y_t=(y_t(0),⋯,y_t(K)), where y_t(k)=∫_0^1dy_t(r,k) for k=0,⋯,K. And y_t is the unique solution to the following differential equation starting at y_0

    y_t=b(y_t)


where b:[0,1]^K+1→ℝ^K+1 is a vector field satisfies

    b(y_t)   =   ∬_]0,1]× [0,...,K][ Λ/r(1_(r,n-1)-1_(r,n))1_n>0+(γ-∑_n∫_0^1 ndy_t(r,n) )(1_(r,n+1)-1_(r,n))1_n<K] dy_t(r,n),

or componentwise

    b(y_t)(0)=-∫_0^1(γ-∑_n∫_0^1 ndy_t(r,n) )dy_t(r,0)_return a bike to a no-bike station+∫_0^1Λ/rdy_t(r,1)_retrieve a bike from a 1-bike station,




    b(y_t)(k)=   ∫_0^1Λ/rdy_t(r,k+1)_retrieve a bike from a k+1-bike station-∫_0^1(Λ/r+γ-∑_n∫_0^1 ndy_t(r,n) )dy_t(r,k)_retrieve and return a bike to a k-bike station
       +∫_0^1(γ-∑_n∫_0^1 ndy_t(r,n) )dy_t(r,k-1)_return a bike to a k-1-bike station,


for k=1,...,K-1, and 

    b(y_t)(K)=-∫_0^1Λ/r dy_t(r,K)_retrieve a bike from a K-bike station+∫_0^1(γ-∑_n∫_0^1 ndy_t(r,n) )dy_t(r,K-1)_return a bike to a K-1-bike station.




A similar theorem is given in the paper of <cit.>, however, a proof is not given in their work.  Thus, to make our paper self contained, we provide a full proof of the mean field limit for the convenience of the reader as it is essential for our future results.  Our proof exploits Doob's inequality for martingales and Gronwall's lemma. Moreover, we use Proposition <ref>, Proposition <ref>, and Proposition <ref> in the proof, and they are stated after the proof of Theorem <ref>.


Since Y_t^N is a semi-martingale, we have the following decomposition of Y_t^N ,

    Y_t^N=Y_0^N_initial condition+M_t^N_martingale+∫_0^tβ(Y_s^N)_drift termds

where
Y_0^N is the initial condition and M_t^N is a family of martingales.  Moreover, ∫_0^tβ(Y_s^N)ds is the integral of the drift term where the drift term is given by β: [0,1]^K+1→ℝ^K+1 or

    β(y)   =∑_x≠ y(x-y)Q(y,x)
       = ∑_n,r[1/rNR_max(1_(r,n-1)-1_(r,n))1_n>0+(M/N-∑_n'∑_r'n'y(r',n'))(1_(r,n+1)-1_(r,n))1_n<K]y(r,n).



We want to compare the empirical measure Y_t^N with the mean field limit y_t defined by

    y_t=y_0+∫_0^tb(y_s)ds.


The remaining of the proof of this theorem can be found in the Appendix.




For any stopping time T such that 𝔼(T)<∞, we have

    𝔼(sup_t≤ T|M_t^N|^2)≤ 4𝔼∫_0^Tα(Y_t^N)dt.
  

The proof is found in the Appendix.  




The drift function b(y) given in Equation (<ref>) is a Lipschitz function with respect to the Euclidean norm in ℝ^K+1. 

The proof is found in the Appendix.  




Under Assumption <ref>, we have for any ϵ>0 and s≥ 0,

    lim_N→∞ℙ(|β(Y_s^N)-b(Y_s^N)|>ϵ)= 0.


The proof is found in the Appendix.  





We have proved mean field limit for our bike sharing model.  Our analysis yields that as the number of stations goes towards infinity, we can solve a set of ordinary differential equations to obtain important performance measure information. The performance measures that we can approximate are the mean proportion of empty or saturated stations, and the average number of bikes in circulation.  Moreover, we can analyze how factors such as fleet size and capacity change the value of the performance measures.  

However, just knowing the mean field limit is not enough.  One reason is that we would like to know more about the stochastic variability of the system, i.e. the fluctuations around the mean field limit.  The mean field limit cannot explain the stochastic fluctuations of the BSS and therefore, we need to analyze the BSS in a different way. Thus, in the subsequent section we develop a functional central limit theorem for our bike sharing model, and explain why it is important for understanding stochastic fluctuations of bike sharing networks.






§ DIFFUSION LIMIT

In this section, we derive the diffusion limit of our stochastic empirical process bike sharing model.  Diffusion limits are critical for obtaining a deep understanding of the sample path behavior of stochastic processes.  One reason is that diffusion limits describe the fluctuations around the mean field limit and can help understand the variance or the asymptotic distribution of the stochastic process being analyzed.  We define our diffusion scaled bike sharing model by subtracting the mean field limit from the empirical measure process and rescaling it by √(N).  Thus, we obtain the following expression for the diffusion scaled bike sharing empirical process

    D_t^N=√(N)(Y_t^N-y_t) .


Unlike many other ride-sharing systems such as Lyft or Uber, bike sharing programs cannot use pricing as a mechanism for redistributing bikes to satisfy demand in real-time.  For this reason, it is essential to understand the dynamics and behavior of D_t^N.  D_t^N can be useful for describing the probability that the proportion of stations with i bikes exceeds a threshold i.e. ℙ(Y^N_t > x) for some x ∈ [0,1]^K+1.  It also describes this probability in a situation where there is no control or rebalancing of bikes in the system.  This knowledge of the uncontrolled system is especially important for newly-started bike sharing systems who are still in the process of gathering information about the system demand. The diffusion limit helps managers of BSS to understand the system dynamics and stability, which in turn helps them make short term and long term managerial decisions. It is also helpful in the case when the operators of the bike sharing system have no money for rebalancing the system to meet  real-time demand.  


Using the semi-martingale decomposition of Y^N_t given in Equation (<ref>), we can write a similar decomposition for D_t^N as follows:

    D_t^N   =√(N)(Y_0^N-y_0)+√(N)M_t^N+∫_0^t√(N)[β(Y_s^N)-b(y_s)]ds
       =D_0^N+√(N)M_t^N+∫_0^t√(N)[β(Y_s^N)-b(Y_s^N)]ds+∫_0^t√(N)[b(Y_s^N)-b(y_s)]ds.


Define 

    D_t=D_0+∫_0^tb'(y_s)D_s ds+M_t

where b'(y)=(∂ b(y)(i)/∂ y(j))_ij∈ℝ^(K+1)× (K+1) and M_t=(M_t(0),⋯,M_t(K))∈ℝ^K+1 is a real continuous centered Gaussian martingale, with Doob-Meyer brackets given by

    M(k) _t   =   ∫_0^t(b_+(y_s)(k)+b_-(y_s)(k))ds,
     M(k),M(k+1) _t   =   -∫_0^t[∫_0^1Λ/rdy_s(r,k+1). 
          + .∫_0^1(γ-∑_n∫_0^1ndy_s(r,n)) dy_s(r,k)]ds   for  k<K,
     M(k),M(j) _t   =   0  for  |k-j|>1.

Here b_+(y)=max(b(y),0) and b_-(y)=-min(b(y),0) denote the positive and the negative parts of function b(y) respectively.

Now we state the functional central limit theorem for the empirical measure process as follows,

Consider D_t^N in 𝔻(ℝ_+,ℝ^K+1) with the Skorokhod J_1 topology, and suppose that 


  1)
lim sup_N→∞√(N)(min_iλ_i^N-Λ)<∞,

  2)
lim sup_N→∞√(N)(M/N-γ)< ∞.

Then if D_0^N converges in distribution to D_0, then D_t^N converges to the unique OU process solving D_t=D_0+∫_0^tb'(y_s)D_s ds+M_t in distribution.







To prove Theorem <ref>, we take the following 4 steps,



  1).√(N)M_t^N is a family of martingales independent of D_0^N with Doob-Meyer brackets given by

    √(N)M^N(k)_t   =   ∫_0^t(β_+(Y_s^N)(k)+β_-(Y_s^N)(k))ds,
    √(N)M^N(k),√(N)M^N(k+1)_t   =   -∫_0^t∑_r[1/rNR_maxY_s^N(r,k+1).
          .+(M/N-∑_n'∑_r'n'Y_s^N(r',n'))Y_s^N(r,k)] ds   for  k<K,
    √(N)M^N(k),√(N)M^N(j)_t   =   0  for  |k-j|>1.



  2).For any T≥ 0, 
    lim sup_N→∞𝔼(|D_0^N|^2)<∞⇒lim sup_N→∞𝔼(sup_0≤ t≤ T|D_t^N|^2)<∞.


  3). If (D_0^N)_N=1^∞ is tight then (D^N)_N=1^∞ is tight and its limit points are continuous.

  4).	If D_0^N converges to D_0 in distribution, then D_t^N converges to the unique OU process solving D_t=D_0+∫_0^tb'(y_s)D_s ds+M_t in distribution.






√(N)M_t^N is a family of martingales independent of D_0^N with Doob-Meyer brackets given by

    √(N)M^N(k)_t   =   ∫_0^t(β_+(Y_s^N)(k)+β_-(Y_s^N)(k))ds,
    √(N)M^N(k),√(N)M^N(k+1)_t   =   -∫_0^t∑_r[1/rNR_maxY_s^N(r,k+1).
          .+(M/N-∑_n'∑_r'n'Y_s^N(r',n'))Y_s^N(r,k)] ds   for  k<K,
    √(N)M^N(k),√(N)M^N(j)_t   =   0  for  |k-j|>1.



The proof is found in the Appendix. 



For any s≥ 0,

    lim sup_N→∞√(N)|β(Y_s^N)-b(Y_s^N)|<∞.


The proof is found in the Appendix. 




For any T≥ 0, if 
    lim sup_N→∞𝔼(|D_0^N|^2 ) < ∞ ,
 then we have 
    lim sup_N→∞𝔼(sup_0≤ t≤ T|D_t^N|^2 ) < ∞ .





The proof is found in the Appendix. 



If (D_0^N)_N=1^∞ is tight then (D^N)_N=1^∞ is tight and its limit points are continuous.


The proof is found in the Appendix. 






b(y) is continuously differentiable with the derivatives 
∂ b(y(r,i))/∂ y(r,j) as follows,

    ∂ b(y(r,0))/∂ y(r,j)   =   j· y(r,0)+Λ/r1_{j=1}-(γ-∑_n=0^Kn(∑_r'y(r',n)))1_{j=0},
    ∂ b(y(r,k))/∂ y(r,j)   =   j·( y(k)-y(k-1))+Λ/r(1_{j=k+1}-1_{j=k})
          +(γ-∑_n=0^Kn(∑_r'y(r',n)))(1_{j=k-1}-1_{j=k})  for  0<k<K,
    ∂ b(y(r,K))/∂ y(r,j)   =   -j· y(K-1)-Λ/r1_{j=K}+(γ-∑_n=0^Kn(∑_r'y(r',n)))1_{j=K-1}.





The above equations can be obtained by directly taking derivatives to Equation (<ref>). We can see that the derivatives of b(y) are linear in y. Thus we can conclude that b(y) is continuously differentiable with respect to y.



By Theorem 4.1 in Chapter 7 of <cit.>,  it suffices to prove that the following condition holds

    sup_t≤ T|∫_0^t{√(N)[b(Y_s^N)-b(y_s)]-b'(y_s)D_s^N}ds| 0.

By Proposition <ref>, we know that b(y_t) is continuously differentiable with respect to y_t. By the mean value theorem, for every 0≤ s≤ t there exists a vector Z_s^N in between Y_s^N and y_s such that

    b(Y_s^N)-b(y_s)=b'(Z_s^N)(Y_s^N-y_s).

Therefore

    ∫_0^t{√(N)[b(Y_s^N)-b(y_s)]-b'(y_s)D_s^N}ds=∫_0^t[b'(Z_s^N)-b'(y_s)]D_s^Nds.

We know that 
    lim_N→∞sup_t≤ T|b'(Z_s^N)-b'(y_s)|=0  in probability

by the mean field limit convergence and the uniform continuity of b'.
By applying Chebyshev inequality we have   that D_s^N is bounded in probability, then by Lemma 5.6 in <cit.>,

    sup_t≤ T|∫_0^t{√(N)[b(Y_s^N)-b(y_s)]-b'(y_s)D_s^N}ds| 0.




The SDE (<ref>) has a unique solution 

    D_t=e^∫_0^tb'(y_s)dsD_0+∫_0^te^∫_s^tb'(y_u)dudM_s.


Define 𝒜(t)=b'(y_t), ℬ(t)=(d/dt M(i),M(j)_t)_ij, then the expectation E(D_t) is

    𝔼[D_t]=e^∫_0^t𝒜(s)ds𝔼[D_0],

 and the covariance matrix Σ (t)=Cov[D_t,D_t] is 
 
    Σ (t)=e^∫_0^t𝒜(s)dsΣ(0)e^∫_0^t𝒜^⊤(s)ds+∫_0^te^∫_s^t𝒜(u)duℬ(s)e^∫_s^t𝒜^⊤(u)duds.

Moreover, differentiation with respect to t yields

    d𝔼[D_t]/dt=𝒜(t)𝔼[D_t],



    dΣ(t)/dt=Σ(t)𝒜(t)^⊤+𝒜(t)Σ(t)+ℬ(t).




The proof is found in the Appendix. 























































































































§ EXTENSIONS
 
In Section <ref>, we assumed WLOG that the routing probabilities and capacities are uniform throughout all stations. If one views our statement as assumptions, then it seems like our model is limiting, however, we emphasize that these are not assumptions and that our analysis extends to the broader case of non-uniform routing probabilities and capacities.  Thus, in the non-uniform setting, we are still able to prove the same fluid and diffusion limits and more importantly reduce the dimensionality of the bike sharing network.  



 §.§ Extensions to non-uniform routing probabilities  and capacities



The first possible extension of our current model is to incorporate non-uniform	routing probabilities, and take into account the origin-destination pairs. That is, after a customer picks up a bike from station i, the probability that he will drop off at station j is equal to P_ij. However, to preserve the Markovian property of the queueing process, we also need to track the number of bikes in use which originated from station i, denoted as U^N_t(i), in addition to the empirical measure process Y^N_t. As the scale of the system N goes to infinity, the dimensionality of the new queueing process (Y^N_t, U^N_t) will also go to infinity, therefore we are back to the same infinite-dimensionality problem with X^N(t), which is illustrated more in details in Section <ref>, and lose the benefit of finite-dimensionality we get by studying the empirical measure Y^N_t.

However, we can still extend our model to non-uniform routing probabilities and capacities without considering the origin-destination pairs. In particular, we assume the probability that a bike being dropped off at station i is equals to P_i, and the capacity at station i is equal to K_i. 

We define relative routing probability at station i as 

    p_i=P_i/max_iP_i.


We consider the empirical measure process Y^N_t(r,n,p,k) that counts the proportion of stations that have n bikes, relative utilization r, routing probability p, and capacity k,

    Y^N_t(r,n,p,k)=1/N∑_i^N1{X_i^N(t)=n, r_i^N=r, p_i^N=p,K_i^N=k}.




Conditioning on Y_t^N(r,n,p,k)=y(r,n,p,k), the transition rates of y are specified as follows:

When a customer arrives to a station with n bikes,  relative utilization r, relative routing probability p, and capacity k, to retrieve a bike, the
proportion of stations having n bikes goes down by 1/N, the proportion of stations having n-1 bikes goes up by 1/N, and the transition rate Q^N is 

    Q^N(y,y+1/N(1_(r,n-1,p,k)-1_(r,n,p,k)) )    =    
    y(r,n,p,k)λ_r N 1_n>0
       =    y(r,n,p,k)μ P/RN 1_n>0
       =    y(r,n,p,k)pP^N_max/rR^N_maxN 1_n>0.

When a customer returns a bike to a station with n bikes ,  relative utilization r, relative routing probability p, and capacity k, to retrieve a bike,, the proportion of stations having n bikes goes down by 1/N, the proportion of stations having n+1 bikes goes up by 1/N, and the transition rate Q^N is 

    Q^N(y,y+1/N(1_(r,n+1,p,k)-1_(r,n,p,k)) )    =   
     y(r,n,p,k) N·μ· P(M-∑_n'∑_r',p',k'n'y(r',n',p',k')N)1_n<k
       =   
     y(r,n,p,k)N^2pP^N_max(M/N-∑_n'∑_r',p',k'n'y(r',n',p',k'))1_n<k.


We have the following functional forward equations for Y_t^N(r,n,p,k).


 For any integrable function f: [0,1]^k+1→ℝ, and 
Y_t^N(r,p,k)=(Y_t^N(r,0,p,k),⋯,Y_t^N(r,k,p,k)) satisfies the following functional forward equation,
 
    𝔼(f(Y_t^N(r,p,k)) | Y_0^N(r,p,k) = y_0(r,p,k)]
       =   ∑_n=0^k𝔼[(f(Y_t^N(r,p,k)+1/N(1_r,n-1,p,k-1_r,n,p,k))-f(Y_t^N(r,p,k)))Y^N_t(r,n,p,k)NpP^N_max/rR^N_max1_n>0]
       +   ∑_n=0^k𝔼[(f(Y_t^N(r,p,k)+1/N(1_r,n+1,p,k-1_r,n,p,k))-f(Y_t^N(r,p,k)))Y^N_t(r,n,p,k)N^2pP^N_max(M/N..
       -   ..∑_n'∑_r',p',k'n'Y_t^N(r',n',p',k'))1_n<k]

 
 
 


  §.§.§ Mean Field Limit

We first state the following assumptions that we use throughout this section,

 There exists a probability measure I(r,p,k)  on [0,1]× [0,1] ×ℕ and Λ > 0 such that, as N tends to infinity, we have


  i) 1/N∑_i=1^N1 _(r_i^N,p_i^N,K_i^N)⇒ I(r,p,k),

  ii)P^N_max/R_max^N→Λ, NP^N_max→𝒫, M/N→γ,

  iii)The set 𝒦=⋃_N=1^∞{K_i^N,i=1,...,N} is finite.





Now we state the main theorem in this section that proves the convergence of empirical process to its mean field limit.

Let |· | denote the Euclidean norm in ℝ^K_max+1. Under Assumption <ref>, suppose that Y_0^N y_0,  then we have for any ϵ>0 and t_0>0,

    lim_N→∞ℙ(sup_t≤ t_0|Y_t^N-y_t|>ϵ)=0.

Here y_t=(y_t(0),⋯,y_t(K_max)), where y_t(j)=∑_k∈𝒦∫_0^1∫ _0^1dy_t(r,j,p,k) for j=0,⋯,K_max. And y_t is the unique solution to the following differential equation starting at y_0

    y_t=b(y_t)


where b:[0,1]^K+1→ℝ^K_max+1 is a vector field satisfies

    b(y_t)   =   ∑_k∈𝒦∑_n=0^k∬_]0,1]× [0,1][ pΛ/r(1_(r,n-1,p,k)-1_(r,n,p,k))1_n>0+.
          .p𝒫(γ-∑_n∑_k≥ n∬_]0,1]×[0,1] ndy_t(r,n,p,k) )(1_(r,n+1,p,k)-1_(r,n,p,k))1_n<k] dy_t(r,n,p,k),

or componentwise


    b(y_t)(0)   =   -∑_k∈𝒦∬_]0,1]× [0,1]p𝒫(γ-∑_n∑_k∈𝒦∬_]0,1]× [0,1] ndy_t(r,n,p,k) )dy_t(r,0,p,k)_return a bike to a no-bike station
          +∑_k∈𝒦∬_]0,1]× [0,1]pΛ/rdy_t(r,1,p,k)_retrieve a bike from a 1-bike station,




    b(y_t)(j)=   ∑_k∈𝒦∬_]0,1]× [0,1]pΛ/rdy_t(r,j+1,p,k)_retrieve a bike from a j+1-bike station
        -∑_k∈𝒦∬_]0,1]× [0,1](pΛ/r+p𝒫(γ-∑_n∑_k∈𝒦∬_]0,1]× [0,1] ndy_t(r,n,p,k) ))dy_t(r,j,p,k)_retrieve and return a bike to a j-bike station
       +∑_k∈𝒦∬_]0,1]× [0,1]p𝒫(γ-∑_n∑_k∈𝒦∬_]0,1]× [0,1] ndy_t(r,n,p,k) )dy_t(r,j-1,p,k)_return a bike to a j-1-bike station,


for j=1,...,K_max-1, and 

    b(y_t)(K_max)   =   -∑_k∈𝒦∬_]0,1]× [0,1]pΛ/r dy_t(r,K_max,p,k)_retrieve a bike from a K_max-bike station
          +∑_k∈𝒦∬_]0,1]× [0,1]p𝒫(γ-∑_n∑_k∈𝒦∬_]0,1]× [0,1] ndy_t(r,n,p,k) )dy_t(r,K_max-1,p,k)_return a bike to a K_max-1-bike station.



The only things we need to prove in this extension case are that b(y) is Lipschitz and that β(Y_t^N) b(Y_t^N) for any t≥ 0. The rest of the proof stays the same as in Section <ref>.





The drift function b(y) given in Equation (<ref>) is a Lipschitz function with respect to the Euclidean norm in ℝ^K_max+1. 

The proof is found in the Appendix.  




Under Assumption <ref>, we have for any ϵ>0 and s≥ 0,

    lim_N→∞ℙ(|β(Y_s^N)-b(Y_s^N)|>ϵ)= 0.


The proof is found in the Appendix.  




  §.§.§ Diffusion Limit

Now we state the functional central limit theorem for the empirical measure process in the extension case as follows,

Consider D_t^N in 𝔻(ℝ_+,ℝ^K_max+1) with the Skorokhod J_1 topology, and suppose that 


  1)
lim sup_N→∞√(N)(P_max^N/R_max^N-Λ)<∞,

  2)
lim sup_N→∞√(N)(M/N-γ)< ∞,

  3) lim sup_N→∞√(N)(NP_max^N-𝒫)< ∞.

Then if D_0^N converges in distribution to D_0, then D_t^N converges to the unique OU process solving D_t=D_0+∫_0^tb'(y_s)D_s ds+M_t in distribution.
Here b'(y)=(∂ b(y)(i)/∂ y(j))_ij∈ℝ^(K_max+1)× (K_max+1) and M_t=(M_t(0),⋯,M_t(K_max))∈ℝ^K_max+1 is a real continuous centered Gaussian martingale, with Doob-Meyer brackets given by

    M(k) _t   =   ∫_0^t(b_+(y_s)(k)+b_-(y_s)(k))ds,
     M(k),M(k+1) _t   =   -∫_0^t[∑_K∈𝒦∬_]0,1]× [0,1]pΛ/rdy_s(r,k+1,p,K). 
          + .∑_K∈𝒦∬_]0,1]× [0,1]p𝒫(γ-∑_n∑_K∈𝒦∬_]0,1]× [0,1]ndy_s(r,n,p,K)) dy_s(r,k,p,K)]ds 
            for  k<K_max,
     M(k),M(j) _t   =   0  for  |k-j|>1.

Here b_+(y)=max(b(y),0) and b_-(y)=-min(b(y),0) denote the positive and the negative parts of function b(y) respectively.


The only things we need to prove in this extension case are listed as propositions below. The rest of the proof stays the same as in Section <ref>.



√(N)M_t^N is a family of martingales independent of D_0^N with Doob-Meyer brackets given by

    M(k) _t   =   ∫_0^t(b_+(y_s)(k)+b_-(y_s)(k))ds,
     M(k),M(k+1) _t   =   -∫_0^t[∑_K∈𝒦∬_]0,1]× [0,1]pΛ/rdy_s(r,k+1,p,K). 
          + .∑_K∈𝒦∬_]0,1]× [0,1]p𝒫(γ-∑_n∑_K∈𝒦∬_]0,1]× [0,1]ndy_s(r,n,p,K)) dy_s(r,k,p,K)]ds 
            for  k<K_max,
     M(k),M(j) _t   =   0  for  |k-j|>1.





 For any s≥ 0, 

    lim sup_N→∞√(N)|β(Y_s^N-b(Y_s^N)|< ∞


The proof is found in the Appendix.  




For k<K_max

    lim_N→∞ℙ(sup_t≤ T|√(N)M^N(k),√(N)M^N(k+1)_t-  M(k),M(k+1)_t|>ϵ)=0


The proof is found in the Appendix.  







b(y) is continously differentiable with the derivatives 
∂ b(y(r,i,p,K))/∂ y(r,j,p,K) as follows,

    ∂ b(y(r,0,p,K))/∂ y(r,j,p,K)   =   p𝒫j· y(r,0,p,K)+pΛ/r1_{j=1}-p𝒫(γ-∑_n=0^K_maxn(∑_r',p',K'y(r',n,p',K')))1_{j=0},
    ∂ b(y(r,k,p,K))/∂ y(r,j,p,K)   =   p𝒫 j·( y(r,k,p,K)-y(r,k-1,p,K))+pΛ/r(1_{j=k+1}-1_{j=k})
          +p𝒫(γ-∑_n=0^K_maxn(∑_r',p',K'y(r',n,p',K')))(1_{j=k-1}-1_{j=k})  for  0<k<K_max,
    ∂ b(y(r,K,p,K))/∂ y(r,j,p,K)   =   -p𝒫 j· y(K_max-1)-pΛ/r1_{j=K_max}
          +p𝒫(γ-∑_n=0^K_maxn(∑_r',p',K'y(r',n,p',K')))1_{j=K_max-1}.





The above equations can be obtained by directly taking derivatives to Equation (<ref>). We can see that the derivatives of b(y) are linear in y. Thus we can conclude that b(y) is continuously differentiable with respect to y.






 §.§ Extensions to adding bike repositioning

We can also extend our model to incorporate repositioning of bikes, which is adopted by most bike-sharing companies to rebalance the network over time. To avoid cumbersome notations, we restrict ourselves to the uniform routing probabilities and capacities scenario for the model discussed in this section.

We may assume that "rebalancers" will arrive to the system according to a independent Poisson process with rate λ_RY^N_t(K)N, where λ_R>0 is a constant and Y^N_t(K) is the proportion of full stations at time t. They pick a full station uniformly at random among all full stations, and remove a certain number C of the bikes from that station. 

In terms of repositioning these bikes, we assume that C number of bikes are being put back to stations with bikes ranging from 0 to K-C. The probability of repositioning to a station with j bikes is equal to Y^N_t(j)g_j/∑_i=0^K-CY^N_t(i)g_i, where {g_i}_i are positive constants and g_0>>g_1>>g_2⋯>>g_K-C>0. This is a smoothed version of choosing the station with minimum number of bikes when repositioniong, and it preserves the Lipschitz property of the drift function, which is essential for deriving the fluid limit and the diffusion limit. We assume additionally that M<N(K-C) to avoid the case where the denominator ∑_i=0^K-CY^N_t(i)g_i can be 0. This can be easily shown through a proof by contradiction. For simplicity of the model, we assume repositioning time are exponentially distributed with rate μ_R. Finally, we use R^N_t to denote the number of bikes currently in repositioning by the rebalancers at time t.

Conditioning on (Y_t^N(r,n),R^N_t)=(y(r,n),R), , we have the following transition rate for the four types of events:

When a customer arrives to a station with n bikes,  relative utilization r, to retrieve a bikes, the proportion of stations having n bikes goes down by 1/N, the proportion of stations having n-1 bikes goes up by 1/N, and the transition rate Q^N is 

    Q^N((y,R),(y+1/N(1_(r,n-1)-1_(r,n)),R) )    =   y(r,n)/rR_max^N1_n>0.



When a customer returns a bike to a station with n bikes and relative utilization r, the proportion of stations having n bikes goes down by 1/N, the proportion of stations having n+1 bikes goes up by 1/N, and the transition rate Q^N is 

    Q^N((y,R),(y+1/N(1_(r,n+1)-1_(r,n)),R) )= y(r,n) ·μ·(M-∑_n'∑_r'n'y(r',n')N-R)1_n<K.



When a rebalancer arrives to a station with K bikes,  relative utilization r, to retrieve C bikes, the proportion of stations having K bikes goes down by 1/N, the proportion of stations having K-C bikes goes up by 1/N, and the transition rate Q^N is 

    Q^N((y,R),(y+1/N(1_(r,K-C)-1_(r,K)),R+C) )    =    
    y(r,K)λ_R N.


When a rebalancer arrives to a station with j bikes and relative utilization r, where 0≤ j≤ K-C, to put back C bikes, the proportion of stations having j bikes goes down by 1/N, the proportion of stations having j+C bikes goes up by 1/N, and the transition rate Q^N is 

    Q^N((y,R),(y+1/N(1_(r,j+C)-1_(r,j)),R-C) )=μ_R y(r,j)g_j/∑_i=0^K-Cy(r,i)g_i·R/C.



We have the following functional forward equations for (Y_t^N(r,n),R).


 For any integrable function f: [0,1]^K+1× [0,∞)→ℝ, 
Y_t^N(r)=(Y_t^N(r,0),⋯,Y_t^N(r,K)) satisfies the following functional forward equation,
 
    𝔼(f(Y_t^N(r),R) | Y_0^N(r) = y_0(r), R_0^N=R_0]
       =   ∑_n=0^K𝔼[(f(Y_t^N(r)+1/N(1_r,n-1-1_r,n),R_t^N)-f(Y_t^N(r),R_t^N))Y^N_t(r,n)/rR^N_max1_n>0]
       +   ∑_n=0^K𝔼[(f(Y_t^N(r)+1/N(1_r,n+1-1_r,n),R_t^N)-f(Y_t^N(r)),R_t^N)Y^N_t(r,n)N(M/N..
       -   ..∑_n'∑_r'n'Y_t^N(r',n')-R^N_t/N)1_n<K]
       +   𝔼[(f(Y_t^N(r)+1/N(1_r,K-C-1_r,K),R^N_t+C)-f(Y_t^N(r),R^N_t))Y^N_t(r,K)λ_R N]
       +   ∑_n=0^K-C𝔼[(f(Y_t^N(r)+1/N(1_r,n+C-1_r,n),R^N_t-C)-f(Y_t^N(r),R^N_t))μ_R Y^N_t(r,n)g_n/∑_i=0^K-CY^N_t(r,i)g_i·R^N_t/C].

 




  §.§.§ Mean Field Limit

Now we show the mean field limit result for the empirical process with rebalancing.

Let |·| denote the Euclidean norm in ℝ^K+2. Under Assumption <ref>, suppose that (Y_0^N,R_0^N/N) (y_0,r̅_0) ,  then we have for any ϵ>0 and t_0>0,

    lim_N→∞P(sup_t≤ t_0|(Y_t^N,R^N_t/N)-(y_t,r̅_t)|>ϵ)=0.

Here y_t=(y_t(0),⋯,y_t(K)), where y_t(k)=∫_0^1dy_t(r,k) for k=0,⋯,K. And (y_t,r̅_t) is the unique solution to the following differential equation starting at (y_0,r̅_0)

    y_t   =   b_1(y_t,r̅_t),
    r̅_t   =   b_2(y_t,r̅_t),

where b_1:[0,1]^K+1→ℝ^K+1 and b_2:[0,∞)→ℝ are vector fields satisfy

    b_1(y_t,r̅_t)
       =   ∬_]0,1]× [0,...,K][ Λ/r(1_(r,n-1)-1_(r,n))1_n>0+(γ-∑_n∫_0^1 ndy_t(r,n)-r̅_t )(1_(r,n+1)-1_(r,n))1_n<K] dy_t(r,n)
       +   ∫_]0,1]λ_R 1_(r,K-C-1_(r,K))dy_t(r,K)+∬_]0,1]× [0,...,K-C]μ_R g_n/∑_i=0^K-Cy_t(r,i)g_i·r̅_t/C (1_(r,n+C)-1_(r,n)) dy_t(r,n),

or componentwise

    b_1(y_t,r̅_t)(0)   =   -∫_0^1(γ-∑_n∫_0^1 ndy_t(r,n) -r̅_t )dy_t(r,0)_return a bike to a no-bike station+∫_0^1Λ/rdy_t(r,1)_retrieve a bike from a 1-bike station
          - ∫_0^1μ_R g_0/∑_i=0^K-Cy_t(r,i)g_i·r̅_t/Cdy_t(r,0)_rebalancer put back C bikes to a 0-bike station,


    b_1(y_t,r̅_t)(k)=   ∫_0^1Λ/rdy_t(r,k+1)_retrieve a bike from a k+1-bike station-∫_0^1(Λ/r+γ-∑_n∫_0^1 ndy_t(r,n) -r̅_t )dy_t(r,k)_retrieve and return a bike to a k-bike station
       +∫_0^1(γ-∑_n∫_0^1 ndy_t(r,n) -r̅_t )dy_t(r,k-1)_return a bike to a k-1-bike station-∫_0^1μ_R g_k/∑_i=0^K-Cy_t(r,i)g_i·r̅_t/Cdy_t(r,k)_rebalancer put back C bikes to a k-bike station,

for k=1,...,K-C-1, and 

    b_1(y_t,r̅_t)(k)=   ∫_0^1Λ/rdy_t(r,k+1)_retrieve a bike from a k+1-bike station-∫_0^1(Λ/r+γ-∑_n∫_0^1 ndy_t(r,n) -r̅_t )dy_t(r,k)_retrieve and return a bike to a k-bike station
       +∫_0^1(γ-∑_n∫_0^1 ndy_t(r,n) -r̅_t )dy_t(r,k-1)_return a bike to a k-1-bike station+∫_0^1μ_R g_k-C/∑_i=0^K-Cy_t(r,i)g_i·r̅_t/Cdy_t(r,k-C)_rebalancer put back C bikes to a k-C-bike station
        +∫_0^1 Cλ_R dy_t(r,K)1_k=K-C_rebalancers retrieve C bikes from a K-bike station,

for k=K-C,...,K-1, and 

    b_1(y_t,r̅_t)(K)   =   -∫_0^1Λ/r dy_t(r,K)_retrieve a bike from a K-bike station+∫_0^1(γ-∑_n∫_0^1 ndy_t(r,n) -r̅_t )dy_t(r,K-1)_return a bike to a K-1-bike station
          -∫_0^1 Cλ_R dy_t(r,K)1_k=K-C_rebalancers retrieve C bikes from a K-bike station,

and

    b_2(y_t,r̅_t)=∫_0^1 Cλ_R dy_t(r,K)_rebalancers retrieve C bikes from a K-bike station-∑_n=0^K-C∫_0^1μ_R g_n/∑_i=0^K-Cy_t(r,i)g_ir̅_t dy_t(r,n)_rebalancers put back C bikes to a n-bike station.



This result can be shown by similar techniques used in Section <ref>.




§ NUMERICAL EXAMPLES AND SIMULATION

In this section,  we confirm our theoretical results of our bike sharing model with a stochastic simulation.   We perform simulations with both stationary  and non-stationary arrival rates, which are discussed substantially in the following subsections. 


Figures <ref> and <ref> provide the user patterns of BSS from the historical data of Citi Bike. In Figure <ref>, we can see that the arrival  rate of users to CitiBike stations varies not only on the time of the day, but also on the day of the week. So in the subsequent numerical examples, we consider both stationary and non-stationary arrival rates to provide insights for the behavior of such systems under different demand and usage conditions.  However, unlike the arrival rate, in Figure <ref>, we observe that the mean travel duration of Citi Bike users does not vary significantly as a function of the time of day or the day of the week.  Therefore, in the subsequent numerical examples we assume a constant mean travel time 1/μ.









 §.§ Simulation Experiments with a Stationary Arrival Rate

 In this section, we provide the results of the simulation studies of our stochastic bike sharing model with a stationary arrival process.  Some of the common model parameters that we use in all of our simulation experiments are given as follows:  N=100, λ=1, μ=1. The number of sample paths we use is 50. Other parameters are specified in the illustration of each figure below.


  
Figure <ref>: Figure <ref> provides the simulation result of the different components of Y_t^N and their mean field limits y_t when K=10. Each station has 5 bikes at the begining of the simulation. The solid lines of Figure <ref> represent the different components of Y_t^N, and the dashed lines of the same color represent the corresponding mean field limit y_t. It is surprising to notice that even though we only have 100 stations in the simulation, the mean field limit provides a high quality approximation to the empirical measure on each component. Furthermore, since arrival rate is constant, the empirical measure converges to its equilibrium point as time tends towards infinity.



  
Figures <ref> - <ref>  : Figures <ref> - <ref> show the simulation results of Y_t^N with a 95% confidence interval, vs. the mean field limit y_t with two  standard deviations of the unscaled diffusion limit, when K=3. The initial distribution of bikes is given as Y_0^N=(0,0.5,0.5,0). In Figure <ref>, the black curve shows Y_t^N(0), the simulated dynamics (average of 50 sample paths) of the proportion of stations with no bikes. The red dashed line shows y_t(0), the mean field limit of the proportion of stations with no bikes over time, which is obtained by solving the system of ODEs derived in Theorem <ref>. The green curves show 2σ(Y_t^N(0)) or two standard deviations of Y_t^N(0) from the sample mean. The blue dashed lines show 2σ(D_t)/N or two standard deviations of the unscaled diffusion limit from the mean field limit y_t(0). Here σ(D_t) is obtained from solving the system of ODEs in Theorem <ref> Equation (<ref>).  We observe that the mean field limit and the diffusion limit provide high quality approximations for our stochastic bike sharing model.  The high quality of the approximations serve to provide additional evidence for using the mean field and diffusion limits we derived earlier.  Similarly, Figures <ref>, <ref>, and <ref> show the dynamics of the proportion of stations with 1, 2, and 3 bikes respectively, like Figure <ref> does for the proportion of stations with no bikes.



  
Figures <ref> :  Figure <ref> shows the simulation results of the variance of D_t^N vs. the numerical solution for the system of ODEs regarding the covariance matrix of D_t (see Equation (<ref>)) when K=3. Again, the approximation of the variance of diffusion limit to the variance of the actual diffusion process is pretty good. 




  
Figure <ref>:  Figure <ref> shows the number of bikes in circulation over time, which is denoted as 

    C_t^N≜ M-∑_j=0^Kj· Y_t^N(j)N.
 
We use N(γ-∑_j=0^Kj· y_t(j)) to approximate 𝔼[C_t^N], the expectation of the number of bikes in circulation. We use 

    Var[M-∑_j=0^Kj· Y_t^N(j)N ] 
       =   Var[∑_j=0^Kj· Y_t^N(j)N ]
       =    N^2[∑_j=0^Kj^2 ·Var[Y^N_t(j)]+∑_i=0^K∑_j=i+1^K2ij·Cov[Y^N_t(i),Y^N_t(j)]]
       =    N^2[∑_j=0^Kj^2 ·Var[y_t(j)+1/√(N)D_t^N(j)].
          .+∑_i=0^K∑_j=i+1^K2ij·Cov[y_t(i)+1/√(N)D_t^N(i),y_t(j)+1/√(N)D_t^N(j)]]
       ≈    N[∑_j=0^Kj^2 ·Var[D_t(j)]+∑_i=0^K∑_j=i+1^K2ij·Cov[D_t(i),D_t(j)]]

to approximate Var[C^N_t], the variance of the number of bikes in circulation. The black curve shows the simulated result of C_t^N, the average number of bikes in circulation.  The red dashed line shows the mean field limit of the number of bikes in circulation over time, which is obtained by solving the ODEs derived in Theorem <ref>. The green curves show 2σ(C_t^N) or two standard deviations of the simulated (50 sample paths) number of bikes in circulation from its sample mean. The blue dashed lines show our approximation of  2σ(C_t^N)  using the diffusion limit (see Equation (<ref>)), which can be computed from solving the  ODEs derived in Theorem <ref>  Equation (<ref>). We also see that the mean field limit and the diffusion limit are good approximations for the actual dynamics of the average number of bikes in circulation. In this specific example, the number of bikes in circulation in equilibrium is roughly 60.





  
Figure <ref>:  Figure <ref> shows that the effect of increasing the arrival rate on the average number of bikes across stations. From Figure <ref>, we see that the average number of bikes decreases as arrival rate increases, and the intuition is that there will be fewer bikes at stations when more people arrive at the stations to retrieve bikes. From an operator's perspective, this plot illustrates that operators of BSS should be more aware of empty stations when usage of the bike sharing platform increases, and vice versa.



  
Figure <ref>: Figure <ref> shows the bike distribution in equilibrium for different values of γ when K=20 . We can see that in equilibrium, the shape of the bike distribution depends on the value of γ. The distribution is right-skewed when γ is small, and left-skewed when γ is large. In particular, the distribution is roughly symmetric when γ is equal to 11.








































 §.§ Simulation Experiments with a Non-stationary Arrival Rate

In this section, we provide the results of the simulation studies of our stochastic bike sharing model with a non-stationary arrival process.  We perform our simulation with the following parameters: N=100, λ(t)=1+0.5sin(t/2), μ=1. The number of sample paths we use in our simulation is 50. Other parameters are specified in the illustration of each figure below. For the time-varying arrival rate, we use a periodic function λ(t)=1+0.5sin(t/2) to mimic the patterns of morning rush and the afternoon-evening rush that we observe from the Citi Bike data (See Figure <ref>).


  
Figure <ref>:

 In Figure <ref>, we simulate the different components of Y_t^N and their subsequent mean field limits y_t when K=10. Each station has 5 bikes at the beginning of the simulation. The solid lines represent the simulation results of different components of Y_t^N, and the dashed lines of the same color represent the corresponding mean field limit y_t. We also add the purple dashed line to represent the time-varying arrival rate. We observe that the mean field limit provides an accurate approximation to the empirical measure on each component.  We also notice that the empirical measure is indeed affected by the non-stationary arrival rate and is also non-stationary. 



  
Figures <ref> - <ref>:

Figures <ref> - <ref> show the simulation results of Y^N_t with a 95% confidence interval, vs. the mean field limit y(t) with  two  standard deviations of the unscaled diffusion limit, when K=3. The initial distribution of bikes is given as Y_0^N=(0,0.5,0.5,0). In Figure <ref>, the black curve shows Y_t^N(0), the simulated proportion of stations with no bikes over time. To produce the simulation, we take an average of 50 independent sample paths. The red dashed line shows y_t(0), the mean field limit of the proportion of stations with no bikes over time, which is obtained from solving the system of ODEs in Theorem <ref>. The green curves show 2σ(Y_t^N(0)) or two standard deviations of Y_t^N(0), from the 50 sample paths in simulation. Lastly, the blue dashed lines show 2σ(D_t)/N or two standard deviations of the unscaled diffusion limit, which is obtained from solving the system of ODEs from Theorem <ref>  Equation (<ref>). We also add a purple dashed line to represent the time-varying arrival rate, with black dashed vertical lines to indicate the peaks and valleys of the arrival rate and the corresponding mean field limit.  Again, we observe that the mean field limit and the diffusion limit provide an accurate approximation to the simulation results.  This accuracy serves to validate the correctness of the mean field and diffusion limits we proved earlier. Moreover, with a small time lag, the simulation results of Y_t^N(0) are positively associated with movements of the arrival rate. Figures <ref>,  <ref>, and <ref> show the dynamics of the proportion of stations with 1, 2, and 3 bikes respectively. Unlike the case of zero bikes, we observe that the dynamics are negatively associated with movements of the arrival rate.  To observe these dynamics in one graph, Figure <ref> combines Figures <ref> - <ref> in one graph, but only keeps the mean field limit curves.   This allows the reader to visualize the dynamics of the mean field empirical measure under a time-varying arrival rate.



  
Figure <ref>  : Figure <ref> shows the simulation results of the variance of D_t^N vs. the numerical solution for the system of ODEs regarding the covariance matrix of D_t (see Equation (<ref>)) when K=3. The approximation of the variance of diffusion limit to the variance of the actual diffusion process is quite accurate. We also observe that as arrival rate increases, the variance of the proportion of stations with no bikes increases, while the proportion of stations with 1, 2, or 3 bikes decreases, and vice versa.



  
Figure <ref>  :  Figure <ref> shows the the average number of bikes in circulation over time, which is denoted as 

    C_t^N≜ M-∑_j=0^Kj· Y_t^N(j)N.
 
We use N(γ-∑_j=0^Kj· y_t(j)) to approximate 𝔼[C_t^N], the expectation of the number of bikes in circulation. We use 

    Var[M-∑_j=0^Kj· Y_t^N(j)N ] 
       =   Var[∑_j=0^Kj· Y_t^N(j)N ]
       =    N^2[∑_j=0^Kj^2 ·Var[Y^N_t(j)]+∑_i=0^K∑_j=i+1^K2ij·Cov[Y^N_t(i),Y^N_t(j)]]
       =    N^2[∑_j=0^Kj^2 ·Var[y_t(j)+1/√(N)D_t^N(j)].
          .+∑_i=0^K∑_j=i+1^K2ij·Cov[y_t(i)+1/√(N)D_t^N(i),y_t(j)+1/√(N)D_t^N(j)]]
       ≈    N[∑_j=0^Kj^2 ·Var[D_t(j)]+∑_i=0^K∑_j=i+1^K2ij·Cov[D_t(i),D_t(j)]]

to approximate Var[C^N_t], the variance of the number of bikes in circulation. The black curve shows the simulated result of the average number of bikes in circulation.  The red dashed line is the mean field limit of the number of bikes in circulation over time, from solving the system of ODEs obtained by Theorem <ref>. The green curves are 2σ(C_t^N), or two standard deviations from the mean of the number of bikes in circulation we obtain by simulating 50 independent sample paths of the stochastic bike sharing model. Lastly, the blue dashed lines show  our approximation of  2σ(C_t^N)  using the diffusion limit (see Equation (<ref>)), which can be computed from solving the system of ODEs from Theorem <ref> Equation (<ref>). We also add a purple dashed line to represent the time-varing arrival rate, with black dashed vertical lines to indicate the peaks and valleys of the arrival rate and the corresponding mean field limit. Again, we can see that the mean field limit and the diffusion limit provide a high quality approximation for the actual dynamics of the average number of bikes in circulation.  Moreover, we observe two important things.  First, the average number of bikes in circulation is also time varying and fluctuates between 40 bikes and 80 bikes.  Second, the average number of bikes in circulation lags slightly behind the arrival rate, which is a common phenomenon in non-stationary queues.  



  §.§.§ Additional Commentary on the Lag Effect


In the case where the arrival process is non-stationary, we observe a lag effect on the change of empirical measure in response to the change of the arrival rate.   Here we explain what the dynamics of the lag effect are in this non-stationary case.  We observe in Figure <ref> and Figure <ref> that after time 5, when the effect of the transient behavior is reduced,  the proportion of stations with no bikes increases as arrival rate increases. However, the proportion of stations with 3 bikes goes down as arrival rate increases. The intuition is that when more people starting picking up bikes at stations, we end up with more empty stations and less full stations.  Similarly, as the arrival rate decreases, the proportion of stations with no bikes decreases while the proportion of stations with 3 bikes increases. The intuition is that when less people starting picking up bikes at stations, we end up with less empty stations and more full stations.  

We also observe a lag between the peak(valley) of arrival and the peak(valley) of the proportion of stations with no bikes (Figure <ref>). There is also a lag between the peak(valley) of arrival and the valley(peak) of the proportion of stations with 3 bikes (Figure <ref>). This is because it takes time for the empirical process to respond to the change in the arrival rate.  In Figure <ref>, we plot the size of the lag for different values of the service rate μ.  We see that the lag effect has the most impact on the proportion of stations with one bike.  It also has an effect on the proportion of stations with no bikes, 2 bikes and 3 bikes, however, the time lag is much smaller than that of one bike.  We also find that the lag effect does not affect maximums and minimums the same way.  We observe in all of the plots, except for proportion of stations with one bike, that the lag effect is more pronounced for minimums than maximums.  This is especially true in the case of the proportion of stations having two bikes.  











































 §.§ Simulation Experiments with Heterogeneous Arrival Rates and Capacities

 In this section, we provide the results of the simulation studies of our stochastic bike sharing model with heterogeneous arrival rates and capacities. We also extend the scale of the system in our example to make it closer to reality. Some of the common model parameters that we use in all of our simulation experiments are given as follows:  N=800, μ=1, and for each station, K=20,30,40, or 50 with equal probabilities. The number of sample paths we use is 30. Other parameters are specified in the illustration of each figure below.
 


  
Figure <ref> :  Figure <ref> shows the distribution of station capacities of Citi Bike (December 2017). We observe that most station capacities range from  20 to 50. Thus we set the parameter K=20,30,40, or 50 with equal probabilities. 



  
Figure <ref> : Figure <ref> shows the average empirical distribution of CitiBike in a day (September 25th, 2017), with one standard deviation. This is exactly the empirical measure y we try to analyze in our bike sharing model. Studying the empirical measure and the its fluctuation over time is of great importance to understanding the system dynamics, providing prediction and guidance for rebalancing and strategic design.



  
Figure <ref> : In Figure <ref>, we simulate the real Citi Bike network (New York City area) and show the empirical measure results from the simulation. The network setup is exactly the same as CitiBike over the month of September, 2017, where we have N=695 number of stations with capacities being the same as Citi Bike. We use Fourier regression to provide a fit for arrival functions for each station, and use them as heterogeneous time-varying arrival rates in our bike sharing model. We used a constant travel time rate μ=3.75, which comes from taking average of the trip duration data. We run the simulation for a 24-hour period and plot the average empirical measure from the simulation. 



  
Figure <ref> : Figure <ref> shows the difference in average empirical measure between simulation and real CitiBike data. We can see that the simulated empirical measure have at most 0.02 difference from the real CitiBike data. We can also observe that stations with 10-28 bikes have a positive difference while stations with less than 10 bikes have a negative difference. This is because rebalancing happens in the real CitiBike network, where bikes will be moved from stations with more bikes to stations with less bikes during peak hours to make the system more balanced, something we are not able to capture in our model. This causes the result from our model to have a larger number of  stations with more bikes, and a smaller number of stations with less stations, compared to reality. Besides, we also didn't consider the effect of information on CitiBike customers in our model, such as smartphone app that tells people the number of bikes and docks at stations in real time. In reality, a lot of people rely on the app to find stations for picking up and dropping off, which drives people to stations that have more bikes or docks. This would cause the same effect in making simulation results slightly different from reality. Overall, we can conclude from the simulation results that the bike sharing model proposed in our paper is doing well in capturing the real-life situation at CitiBike, which shows that our mean field limit and diffusion limit results will be of great value to the understanding the bike sharing systems in real life.



  
Figure <ref> :This figure shows the empirical measure from simulation and its mean field limit in the heterogeneous arrivals and capacities case where arrivals are stationary. The arrival rates are set up to be λ=0.25,0.5,0.75, or 1 with equal probabilities. In this figure we showed the the empirical measure and its corresponding mean field limit with 95% confidence interval at a given time point. The blue bars show the empirical measure from simulation, which is an average of 30 sample paths. The red bars show the corresponding mean field limit we get from solving the ODEs in (<ref>). The error bars show the 95% confidence interval from the simulation results. We can see that the mean field limit is fitting well with empirical measure, which shows that our model is able to capture heterogeneous arrivals and capacities very well, and it can be easily adapted to more complex models, such as those with non-uniform routing probabilities, as shown in Section <ref>.



  
Figure <ref> : This figure shows the evolution of different components of the empirical measure over time in the heteogeneous arrivals and capacities setting with the same parameters setup as in Figure <ref>.  The solid lines represents proportion of stations with k bikes where k=5, 10, 15, 20, 25. The dashed lines represents their corresponding mean field limits. Again the mean field limits are fitting well with the empirical measure.



  
Figure <ref> :This figure shows the empirical measure from simulation and its mean field limit in the heterogeneous arrivals and capacities case where arrivals are non-stationary. The arrival rates are set up to be λ(t)=4λ_0(1+sin(2t)) where λ_0=0.25,0.5,0.75,1 with equal probabilities. In this figure we showed the the empirical measure and its corresponding mean field limit with 95% confidence interval at a given time point. The blue bars show the empirical measure from simulation, which is an average of 30 sample paths. The red bars show the corresponding mean field limit we get from solving the ODEs in (<ref>). The error bars show the 95% confidence interval from the simulation results. We can see that the mean field limit is fitting well with empirical measure, which shows that our model is able to capture heterogeneous arrivals that are also non-stationary very well, and it can be easily adapted to more complex models, such as those with non-uniform routing probabilities, as shown in Section <ref>.




  
Figure <ref> : This figure shows the evolution of different components of the empirical measure over time in the heterogeneous arrivals and capacities setting with the same parameters setup as in Figure <ref>.  The solid lines represents proportion of stations with k bikes where k=0, 10, 20, 30, 40, 50. The dashed lines represents their corresponding mean field limits. We see that the empirical measure is also non-stationary and that the mean field limits are fitting well with the empirical measure.



  
Computational cost:
Another major benefit of our model compared to just using simulation to study the bike sharing system is that our model is extremely computationally inexpensive, given that it only involves numerically solving ODEs. The computational time for running the simulation example in Figure <ref>, with 10 sample paths for 8 unit times is 13.5 hours, while getting the mean field limit for the same example only takes less than 30 seconds.
 
































§ APPLIED AND PRACTICAL VALUE OF OUR WORK

The results presented in this paper have great value to the operations of bike sharing systems. First of all, not only is the empirical measure itself an important performance measure to the bike sharing systems, but it also allows us to obtain salient performance measures such as y_t(0) (the proportion of stations with no bikes), y_t(K) (the proportion of stations that are full), M - ∑^K_j=0 j · y_t(j)N (the number of bikes in circulation), among others.  The empirical measure approach is a significant reduction in computational complexity when compared to the full stochastic model.  

Moreover, it is the first time that a diffusion limit of the empirical measure of a inhomogeneous bike sharing systems is derived, which gives great insights to the fluctuation of the systems performance over time from a queueing and risk management perspective. It is a much more computationally efficient way to study the system behaviors in strategic planning stage, where given the design of the system and parameters such as arrival rate, travel time distribution and fleet size, you can easily get the proportion of problematic stations (empty and full stations) over time, which is a key measure that we want to minimize in a bike sharing system. More importantly, the diffusion limit of the empirical measure provides a refinement to the mean field limit, which gives us a better understanding to system fluctuations over time. By using the diffusion limit, we are able to build confidence intervals for the proportion of problematic stations. This helps us design a BSS with low blocking experiences, not just in expectation, but with high probability.  This is especially important for managers of BSSs who want to control the dynamics of bike stations and reduce the volatility of station fluctuations.

Another benefit of our work is that deriving the mean field limit and diffusion limit provides a way to formulate optimization problems associated with BSS. For example, the mean field and diffusion limit for the proportion of problematic stations can be used as objective functions that we try to minimize given system parameters. The expectation and variance of the number of bikes in circulation can computed through the mean field limit and diffusion limit of the empirical process, and therefore can be used to determine optimal fleet size.  The current literature only uses mean field limits and our diffusion limits can be used to determine optimal fleet size under a more complex stochastic setting instead of a deterministic setting.  

Our analysis also benefits the rebalancing of bike sharing systems. It can be used to provide short-term prediction of the empirical measure of the system, which helps the operators of BSS identify when key measures such as proportion of problematic stations, or number of bikes in circulation, will go beyond a threshold and act beforehand. Different from traditional data analysis methods that predict patterns of BSS solely using history data, our method analyzes the system behaviors from a more fundamental way, one that does not heavily rely on data.  Most importantly, as peak hours only last a few hours, decisions for rebalancing need to be made fast. In this case our method is much more computationally efficient and effective, than just doing simulations, which tends to be very slow and intractable for large scale systems like CitiBike.  

Another major benefit of using the empirical measure approach is that if CitiBike chooses to add a station, then in the empirical measure approach the dimensionality will only increase if the number of docks at the new station is larger than all of the rest.  However, in the individual station model, it will automatically increase by one.  Even though systems like CitiBike are large, they continue to add stations and increase the complexity of simulating the system.  The empirical measure approach that we advocate in this work does not get worse when the management chooses to add stations.


From a broader perspective, the framework of mean field and diffusion limits we established in this paper provide an  effective and efficient way to analyze different problems associated with BSS, such as, designing reasonable architecture of a BSS, finding a better path scheduling, improving inventory management, redistributing the bikes among stations or clusters, price optimization, application of intelligent information technologies and so forth.  Our work serves as the initial step to exploring these important problems facing BSS.   





§ CONCLUSION


In this paper, we construct a bike sharing queueing model that incorporates the finite capacity of stations. Since our model is intractable to analyze directly, especially for a large number of stations, we propose to analyze the limiting dynamics of an empirical process that describes the proportion of stations that have a certain number of bikes.  We prove a mean field limit and a functional central limit theorem for our stochastic bike sharing empirical process, showing that the mean field limit and the variance of the empirical process can be described by a system of 1/2(K+4)(K+1) differential equations where K is the maximum station capacity.  We compare the mean field limit and the functional central limit theorem with simulation and show that the differential equations approximate the mean and variance of the empirical process extremely well. 

There are many directions for future work.  The first direction would be to generalize the arrival and service distribution to follow general distributions.  As Figure <ref> shows, the trip durations are not exponential and are closer to a lognormal distribution.  An extension to general distributions would aid in showing how the non-exponential distributions affect the dynamics of the empirical process.  Recent work by <cit.> provides a Poisson process representation of phase type distributions and Markovian arrival processes.  This work might be useful in deriving new limit theorems for the queueing process with non-renewal arrival and service processes.  

In the non-stationary context, it is not only important to understand the dynamics of the mean field limit, but also it is important to know various properties of the mean field limit. For example, it would be informative to know the size of the amplitude and the frequency of the mean field limit when the arrival rate is periodic. One way to analyze the amplitude and the frequency is to exploit methods from non-linear dynamics like Lindstedt's method and the two-variable expansion method in <cit.>.

Lastly, it is also interesting to consider a spatial model of arrivals to the bike sharing network.  In this case, we would consider customers arriving to the system via a spatial Poisson process and customers would choose among the nearest stations to retrieve a bike.  This spatial process can model the real choices that riders make and would model the real spatial dynamics of bike sharing networks.  We intend to pursue these extensions in future work.   





§ ACKNOWLEDGMENTS

Jamol Pender gratefully acknowledges the support of National Science Foundation (NSF) for Jamol Pender's Career Award CMMI # 1751975.








§ APPENDIX
 


 The time derivative of the expectation 𝔼[f(𝐗(t))] can be derived by the following discretization method. Taking the expectation on f(𝐗(t+Δ)) conditioned on 𝐗(0)=𝐱 for some small Δ>0, we have
 
    𝔼[f(𝐗(t+Δ)|𝐗(0)=𝐱] 
       =   ∑_i=1^Nf(𝐗(t))[1-λ_iΔ1_{X_i(t)>0}-μ P_i(M-∑_k=1^NX_k(t))Δ1_{X_i(t)<K_i}] 
           +∑_i=1^N[f(𝐗(t)-1_i)λ_iΔ1_{X_i(t)>0}+ f(𝐗(t)+1_i)μ P_i(M-∑_k=1^NX_k(t))Δ1_{X_i(t)<K_i}]
          +o(Δ).

 Then 
 
    𝔼[f(𝐗(t+Δ)|𝐗(0)=𝐱]- f(𝐗(t))
       =   ∑_i=1^N[f(𝐗(t)-1_i)-f(𝐗(t)]λ_iΔ1_{X_i(t)>0}
           +∑_i=1^N[f(𝐗(t)+1_j)-f(𝐗(t)]μ P_i(M-∑_k=1^NX_k(t))Δ1_{X_i(t)<K_i}
          +o(Δ).

 By dividing by Δ and taking the expectation on both sides of Equation(<ref>) we get
 
    𝔼[f(𝐗(t+Δ))]-𝔼[f(𝐗(t))]/Δ
       =   ∑_i=1^N𝔼[(f(𝐗(t)-1_i)-f(𝐗(t))λ_i1_{X_i(t)>0}]
           +∑_i=1^N𝔼[(f(𝐗(t)+1_i)-f(𝐗(t))μ P_i(M-∑_k=1^NX_k(t))1_{X_i(t)<K_i}]
          +o(Δ)/Δ .

 Taking Δ→ 0 yields
 
    d/dt𝔼[f(𝐗(t))]
       =   ∑_i=1^N𝔼[(f(𝐗(t)-1_i)-f(𝐗(t))λ_i1_{X_i(t)>0}]
           +∑_i=1^N𝔼[(f(𝐗(t)+1_i)-f(𝐗(t))μ P_i(M-∑_k=1^NX_k(t))1_{X_i(t)<K_i}].

 Let f(𝐗(t))=X_i(t) for i=1,⋯,N, we have the following functional forward equations to each component of 𝐗(t),
 
    𝔼[f(X_i(t)) | X_i(0) = x_i]  ≡d/dt𝔼[f(X_i(t)) | X_i(0) = x_i]  
       =   𝔼[ ( f(X_i(t) +1) - f(X_i(t)) ) ·(μ P_i(M-∑_k=1^NX_k(t))1{X_i(t)<K_i}) ]  
          +  𝔼[ ( f(X_i(t) -1) - f(X_i(t)) ) ·( λ_i1{ X_i(t) > 0 }) ] .
 




Let |·| denote the Euclidean norm in ℝ^K+1, then

    |Y_t^N-y_t|   =   |Y_0^N+M_t^N+∫_0^tβ
    (Y_s^N)ds-y_0-∫_0^tb(y_s)ds|
        =   |Y_0^N-y_0+M_s^N+∫_0^t(β
    (Y_s^N)-b(Y_s^N))ds+∫_0^t(b(Y_s^N)-b(y_s))ds|.

Now define the random function f^N(t)=sup_s≤ t|Y_s^N-y_s|, we have

    f^N(t)≤ |Y_0^N-y_0|+sup_s≤ t|M_s^N|+∫_0^t|β(Y_s^N)-b(Y_s^N)|ds+∫_0^t|b(Y_s^N)-b(y_s)|ds.

By Proposition <ref>, b(y) is Lipschitz with respect to Euclidean norm. Let L be the Lipschitz constant of b(y), then

    f^N(t)   ≤    |Y_0^N-y_0|+sup_s≤ t|M_s^N|+∫_0^t|β(Y_s^N)-b(Y_s^N)|ds+∫_0^t|b(Y_s^N)-b(y_s)|ds 
       ≤    |Y_0^N-y_0|+sup_s≤ t|M_s^N|+∫_0^t|β(Y_s^N)-b(Y_s^N)|ds+L∫_0^t|Y_s^N-y_s|ds 
       ≤    |Y_0^N-y_0|+sup_s≤ t|M_s^N|+∫_0^t|β(Y_s^N)-b(Y_s^N)|ds+L∫_0^tf^N(s)ds.


By Gronwall's lemma (See <cit.>), 

    f^N(t) ≤(|Y_0^N-y_0|+sup_s≤ t|M_s^N|+∫_0^t|β(Y_s^N)-b(Y_s^N)|ds)e^Lt.


Now  to bound f^N(t) term by term, we define the function α: [0,1]^K+1→ℝ^K+1 as

    α(y)   =   ∑_x≠ y|x-y|^2Q(y,x)
       =   1/N∑_n∑_r[1/rNR_max(1_(r,n-1)+1_(r,n))1_n>0.
          +
    . (M/N-∑_n'∑_r'n'y(r',n'))(1_(r,n+1)+1_(r,n))1_n<K] · y(r,n)

and consider  the following four sets 

    Ω_0    =   {|Y_0^N-y_0|≤δ}, 
    Ω_1    =   {∫_0^t_0|β(Y_s^N)-b(Y_s^N)|ds≤δ}, 
    Ω_2    =   {∫_0^t_0α(Y_t^N)dt ≤ A(N)t_0 }, 
    Ω_3    =   {sup_t≤ t_0|M_t^N|≤δ} ,

where δ=ϵ e^-Lt_0/3. Here the set Ω_1 is to bound the initial condition, the set Ω_2 is to bound the drift term β and the limit of drift term b, and  the sets Ω_2,Ω_3 are to bound the martingale M_t^N.

Therefore on the event Ω_0∩Ω_1∩Ω_3,

    f^N(t_0)≤ 3δ e^Lt_0=ϵ.




Since lim_N→∞M/N=γ and lim_N→∞NR_max=1/Λ, we can
choose  N large enough such that

    M/N≤ 2γ,   NR_max≥1/2Λ.

And by the proof of Proposition <ref>, there exists C>0 such that lim_N→∞r_i^N≥Λ/C. See the proof of Proposition <ref> in the Appendix for details. 
Thus

    α(y)    ≤   1/N∑_n∑_r(2Λ/r· 2+2γ· 2 )  · y(r,n)
       ≤   1/N( 4Λ/Λ/C+4γ)∑_n∑_r  y(r,n)
       ≤   4/N(C+γ).

Consider the stopping time 
    T=t_0∧inf{t≥ 0:∫ _0^tα(Y_s^N)ds>A(N)t_0}.

By Proposition <ref>,

    𝔼(sup_t≤ T|M_t^N|^2)≤ 4𝔼∫_0^Tα(Y_t^N)dt≤ 4A(N)t_0.

On Ω_2, we have T =t_0, so 
Ω_2 ∩Ω_3^c⊂{sup_t≤ T|M^N_t|>δ}. By
Chebyshev's inequality, 

    ℙ(Ω_2 ∩Ω_3^c)≤ℙ(sup_t≤ T|M^N_t|>δ)≤𝔼(sup_t≤ T|M_t^N|^2)/δ^2≤ 4A(N)t_0/δ^2.

Thus, by Equation (<ref>), we have the following result,

    ℙ(sup_t≤ t_0|Y_t^N-y_t|>ϵ)   ≤ℙ(Ω_0^c∪Ω_1^c∪Ω_3^c)
       ≤ℙ(Ω_2 ∩Ω_3^c)+ℙ(Ω_0^c∪Ω_1^c∪Ω_2^c)
       ≤ 4A(N)t_0/δ^2+ℙ(Ω_0^c∪Ω_1^c∪Ω_2^c)
       =36A(N)t_0 e^2Lt_0/ϵ^2+ℙ(Ω_0^c∪Ω_1^c∪Ω_2^c).


Let A(N)=4(C+γ)/N, then Ω_2^c=∅.
And since Y_0^N y_0,   lim_N→∞ℙ(Ω_2^c)=0. Therefore we have 
    lim_N→∞ℙ(sup_t≤ t_0|Y_t^N-y_t|>ϵ)=lim_N→∞ℙ(Ω_1^c) .

By Proposition <ref>, lim_N→∞ℙ(Ω_1^c)=0.
Thus, we proved the final result

    lim_N→∞ℙ(sup_t≤ t_0|Y_t^N-y_t|>ϵ)=0.





<ref>[Bounding Martingale]
For any stopping time T such that 𝔼(T)<∞, we have

    𝔼(sup_t≤ T|M_t^N|^2)≤ 4𝔼∫_0^Tα(Y_t^N)dt.
  







Let μ̃ be the jump measure of Y_t^N, and ν be its compensator, defined on (0,∞)× [0,1]  by 

    μ̃=∑_t:Y_t^N≠ Y_t-^Nδ(t,Y_t^N),       ν(dt,B)=Q(Y_t-^N,B)dt   ∀ B∈ℬ([0,1]).

Let Ỹ_m,t^N be the jump chain of Y_t^N, J_m be the jump time, then we have for any t∈ [0,∞), J_n≤ t<J_n+1 for some n≥ 0.
The martingale M_t^N can be written as

    M_t^N   =   Y^N_t-Y_0^N-∫_0^tβ(Y_s^N)ds 
       =   ∑_m=0^n-1(Ỹ^N_m+1-Ỹ^N_m)-∫_0^t∫_0^1(y-Y_s-^N)Q(Y_s-^N,dy)ds
       =   ∫_0^t∫_0^1(y-Y_s-^N)μ̃(ds,dy)-∫_0^t∫_0^1(y-Y_s-^N)ν(ds,dy)
       =   ∫_0^t∫_0^1(y-Y_s-^N)(μ̃-ν)(ds,dy).

Note the following identity

    (M_t^N)^2   =   2∫_0^t∫_0^1M_s-(y-Y_s-^N)(μ̃-ν)(ds,dy)+∫_0^t∫_0^1(y-Y_s-^N)^2μ̃(ds,dy).

This can be established by verifying that the jumps of the left and right hand sides agree, and that their derivatives agree between jump times. Then we can write

    (M_t^N)^2=N_t^N+∫_0^tα(Y_t^N)ds

where

    N_t^N=∫_0^t∫_0^1H(s,y)(μ̃-ν)(ds,dy),

and

    H(s,y)=2M_s-(y-Y_s-^N)+(y-Y_s-^N)^2.

Consider the previsible process

    H_2(t,y)=H(t,y)1_{t≤ T∧ T_n}

where T_n=inf{t≥ 0: β(Y_t^N)>n}∧ n.

Then

    N_T∧ T_n=∫_0^∞∫_0^1H_2(t,y)(μ̃-ν)(dt,dy),

and

    𝔼∫_0^∞∫_0^1|H_2(s,y)|ν(ds,dy)   =   𝔼∫_0^T∧ T_n∫_0^1|2M_s-(y-Y_s-^N)+(y-Y_s-^N)^2|ν(ds,dy) 
       =    𝔼∫_0^T∧ T_n(2|M_t^N|β(Y_t^N)+α(Y_t^N))dt 
       ≤   𝔼∫_0^T∧ T_n2(2+n^2)ndt+𝔼∫_0^Tα(Y_t^N)dt 
       ≤   2n^4+4n^2+4(C+γ)/N𝔼(T)<∞.

By Theorem 8.4 in <cit.>, we have that N^T∧ T_n is a martingale. Replace t by T∧ T_n in Equation (<ref>) and take expectation to obtain

    𝔼(|M_T∧ T_n|^2)=𝔼(N_T∧ T_n)+𝔼∫_0^T∧ T_nα(Y_t^N)dt.

Since N^T∧ T_n is a martingale, 

    𝔼(N_T∧ T_n)=𝔼(N_0)=0.

Apply Doob's L^2-inequality to the martingale M^T∧ T_n to obtain

    𝔼(sup_t≤ T∧ T_n|M_t|^2)   ≤    4𝔼(|M_T∧ T_n|^2)
       =    4𝔼∫_0^T∧ T_nα(Y_t^N)dt.






<ref>[Asymptotic Drift is Lipschitz]
The drift function b(y) given in Equation (<ref>) is a Lipschitz function with respect to the Euclidean norm in ℝ^K+1. 



Assume that max_i(λ_i^N)≤ C
<∞ for all N, 
then 
    r_i^N=R_i^N/R_max=1/Nλ_i^N/max_i(1/Nλ_i^N)=min_iλ_i^N/λ_i^N.

By Assumption <ref>, 
    NR_max=1/min_iλ_i^N→1/Λ>0.

Therefore 
    lim_N→∞r_i^N≥lim_N→∞Λ/max_iλ_i^N≥Λ/C.

Thus the integral that defines function b should start from Λ/C instead of 0, i.e.

    b(y)=∬_[Λ/C,1]× [0,...,K][ Λ/r(1_(r,n-1)-1_(r,n))1_n>0+(γ-∑_n∫ ndy(r,n) )(1_(r,n+1)-1_(r,n))1_n<K] dy(r,n).

Now consider y,ỹ∈ [0,1]^K+1,

    |b(y)-b(ỹ)|≤ 2(Λ/Λ/C+γ)|y-ỹ|=2(C+γ)|y-ỹ|≜ L|y-ỹ|

where |·| denotes the Euclidean norm.



<ref>[Asymptotic Drift is Lipschitz]
The drift function b(y) given in Equation (<ref>) is a Lipschitz function with respect to the Euclidean norm in ℝ^K_max+1. 



Assume that max_i(λ_i^N)≤ C
<∞ for all N, 
then 
    p_i^N/r_i^N=P_i^N/P_max^N/R_i^N/R^N_max=λ_i^NR^N_max/P_max^N

By Assumption <ref>, 
    R^N_max/P^N_max→1/Λ>0.

Therefore 
    lim sup_N→∞p_i^N/r_i^N≤ C/Λ.

Thus the function to be integrated in b is bounded,

    b(y)   =   ∑_k∈𝒦∬_[0,1]× [0,1][ pΛ/r(1_(r,n-1,p,k)-1_(r,n,p,k))1_n>0.
          +.p𝒫(γ-∑_n∑_k∈𝒦∬_[0,1]× [0,1] ndy(r,n,p,k) )(1_(r,n+1,p,k)-1_(r,n,p,k))1_n<k] dy(r,n,p,k).


Now consider y,ỹ∈ [0,1]^K_max+1,

    |b(y)-b(ỹ)|≤ 2(C/Λ·Λ+𝒫γ)|y-ỹ|=2(C+𝒫γ)|y-ỹ|≜ L|y-ỹ|

where |·| denotes the Euclidean norm.




<ref>[Drift is Asymptotically Close to a Lipschitz Drift]
Under Assumption <ref>, we have for any ϵ>0 and s≥ 0,

    lim_N→∞ℙ(|β(Y_s^N)-b(Y_s^N)|>ϵ)= 0.



 
    |β(Y_s)-b(Y_s)|^2
    ≤   2^2∑_n|∑_r1/rNR_maxY_s(r,n)-∫_0^1Λ/rdY_s(r,n)|^2+2^2∑_n|∑_rM/NY_s(r,n)-∫_0^1γ dY_s(r,n)|^2
       +2^2∑_n|∑_r(∑_n'∑_r'n'Y_s(r',n'))Y_s(r,n)-∫_0^1(∑_n'∫_0^1n'dY_s(r',n'))dY_s(r,n)|^2.
 
 It suffices to show each term goes to zero as N→∞.

Since Y_s^N is a discrete random variable, we have for each n,

    ∫_0^1f(r)dY_s(r,n)=∑_rf(r)Y_s(r,n)
  holds for any function f.

Then

    |∑_r1/rNR_maxY_s(r,n)-∫_0^1Λ/rdY_s(r,n)|
    ≤   |∑_r(1/rNR_max-Λ/r) Y_s(r,n)|+|∑_rΛ/rY_s(r,n)-∫_0^1Λ/rdY_s(r,n)|
    
    =    |1/NR_max-Λ|∑_r1/rY_s(r,n)
    ≤   |1/NR_max-Λ| C/Λ∑_rY_s(r,n)
    ≤   |1/NR_max-Λ|C/Λ→ 0.

Similarly,
		
    |∑_rM/NY_s(r,n)-∫_0^1γ dY_s(r,n)|
    ≤   |∑_r(M/N-s) Y_s(r,n)|+|∑_rγ Y_s(r,n)-∫_0^1γ dY_s(r,n)|
    ≤   |M/N-γ|→ 0.
	
The last term is zero since Y_s^N is discrete and 

    ∫_0^1f(r)dY_s(r,n)=∑_rf(r)Y_s(r,n)
 for any function f.




<ref>[Drift is Asymptotically Close to a Lipschitz Drift]
Under Assumption <ref>, we have for any ϵ>0 and s≥ 0,

    lim_N→∞ℙ(|β(Y_s^N)-b(Y_s^N)|>ϵ)= 0.



 
    |β(Y_s)-b(Y_s)|^2
       ≤    2^2∑_n|∑_r,p,kpP^N_max/rR^N_maxY_s(r,n,p,k)-∑_k∈𝒦∬_[0,1]× [0,1]pΛ/rdY_s(r,n,p,k)|^2
          +2^2∑_n|∑_r,p,kpNP^N_maxM/NY_s(r,n,p,k)-∑_k∈𝒦∬_[0,1]× [0,1]p𝒫γ dY_s(r,n,p,k)|^2
          +2^2∑_n|∑_r,p,k(∑_n'∑_r',p',k'n'Y_s(r',n',p',k'))Y_s(r,n,p,k).
          -.∑_k∈𝒦∬_[0,1]× [0,1](∑_n'∑_k∈𝒦∬_[0,1]× [0,1]n'dY_s(r',n',p',k'))dY_s(r,n,p,k)|^2.

It suffices to show each term goes to zero as N→∞.

Since Y_s^N is a discrete random variable, we have for each n,

    ∑_k∈𝒦∬_[0,1]× [0,1]f(r,p,k)dY_s(r,n,p,k)=∑_k∈𝒦∬_[0,1]× [0,1]f(r,p,k)Y_s(r,n,p,k)
  holds for any function f.

Then

    |∑_r,p,kpP^N_max/rR^N_maxY_s(r,n,p,k)-∑_k∈𝒦∬_[0,1]× [0,1]Λ/rdY_s(r,n,p,k)|
    ≤   |∑_r,p,k(pP^N_max/rR^N_max-pΛ/r) Y_s(r,n,p,k)|+|∑_r,p,kpΛ/rY_s(r,n,p,k)-∑_k∈𝒦∬_[0,1]× [0,1]pΛ/rdY_s(r,n,p,k)|
    
    =    |P^N_max/R^N_max-Λ|∑_r,p,kp/rY_s(r,n,p,k)
    ≤   |P^N_max/R^N_max-Λ| C/Λ∑_r,p,kY_s(r,n,p,k)
    ≤   |P^N_max/R^N_max-Λ|C/Λ→ 0.

Similarly,

    |∑_r,p,kpNP^N_maxM/NY_s(r,n,p,k)-∑_k∈𝒦∬_[0,1]× [0,1]p𝒫γ dY_s(r,n,p,k)|
    ≤   |∑_r,p,k(NP^N_maxM/N-𝒫γ) pY_s(r,n,p,k)|+|∑_r,p,kp𝒫γ Y_s(r,n)-∑_k∈𝒦∬_[0,1]× [0,1]p𝒫γ dY_s(r,n,p,k)|
    ≤   |NP^N_maxM/N-𝒫γ|→ 0.
	
The last term is zero since Y_s^N is discrete and 

    ∑_k∈𝒦∬_[0,1]× [0,1]f(r,p,k)dY_s(r,n,p,k)=∑_k∈𝒦∬_[0,1]× [0,1]f(r,p,k)Y_s(r,n,p,k)
 for any function f.



By Dynkin's formula,

    √(N)M^N(k)_t=   ∫_0^tN∑_x≠ Y_s^N|x(k)-Y_s^N(k)|^2 Q(Y_s^N,x)ds
    
    =   N∫_0^tα(Y_s^N)(k)ds
    
    =   ∫_0^t∑_r[1/rNR_max(Y_s^N(r,k+1)1_k<K+Y_s^N(r,k)1_k>0).
       +.(M/N-∑_n'∑_r'n'Y_s^N(r',n'))(Y_s^N(r,k)1_k<K+Y_s^N(r,k-1)1_k>0)]ds
    
    =   ∫_0^t(β_+(Y_s^N)(k)+β_-(Y_s^N)(k))ds.


To compute  √(N)M^N(k),√(N)M^N(k+1)_t for k<K, since

    M^N(k)+M^N(k+1)_t
    
    =   ∫_0^t∑_x≠ Y_s^N|x(k)+x(k+1)-Y_s^N(k)-Y_s^N(k+1)|^2 Q(Y_s^N,x)ds
    
    =   1/N∫_0^t∑_r[1/rNR_max(Y_s^N(r,k+2)1_k<K-1+Y_s^N(r,k)1_k>0).
       +.(M/N-∑_n'∑_r'n'Y_s^N(r',n'))(Y_s^N(r,k+1)1_k<K-1+Y_s^N(r,k-1)1_k>0)]ds.

We have that

    √(N)M^N(k),√(N)M^N(k+1)_t
    
    =   N/2[ M^N(k)+M^N(k+1)_t- M^N(k)_t- M^N(k+1)_t]
    
    =   1/2∫_0^t∑_r[1/rNR_max(Y_s^N(r,k+2)1_k<K-1+Y_s^N(r,k)1_k>0).
       +.(M/N-∑_n'∑_r'n'Y_s^N(r',n'))(Y_s^N(r,k+1)1_k<K-1+Y_s^N(r,k-1)1_k>0)]ds
       -1/2∫_0^t(β_+(Y_s^N)(k)+β_+(Y_s^N)(k+1)+β_-(Y_s^N)(k)+β_-(Y_s^N)(k+1))ds
    
    =   -∫_0^t∑_r[1/rNR_maxY_s^N(r,k+1)+(M/N-∑_n'∑_r'n'Y_s^N(r',n'))Y_s^N(r,k)] ds.

When |k-j|>1, M^N(k) and M^N(j) are independent, thus 

    √(N)M^N(k),√(N)M^N(j)_t=0.



<ref>
For any s≥ 0,

    lim sup_N→∞√(N)|β(Y_s^N)-b(Y_s^N)|<∞.



 
    |β(Y_s)-b(Y_s)|
    ≤   2∑_n|∑_r1/rNR_maxY_s(r,n)-∫_0^1Λ/rdY_s(r,n)|+2∑_n|∑_rM/NY_s(r,n)-∫_0^1γ dY_s(r,n)|
       +2∑_n|∑_r(∑_n'∑_r'n'Y_s(r',n'))Y_s(r,n)-∫_0^1(∑_n'∫_0^1n'dY_s(r',n'))dY_s(r,n)|.
 
 By Equation (<ref>) and (<ref>), 
 
    |β(Y_s)-b(Y_s)|≤ 2(K+1)(|1/NR_max-Λ|C/Λ+|M/N-γ|).

 By the assumptions in Theorem <ref>, we have
 
    lim sup_N→∞√(N)(min_iλ_i^N-Λ)<∞,  lim sup_N→∞√(N)(M/N-γ)< ∞.

Thus 

    lim sup_N→∞√(N)|β(Y_s^N)-b(Y_s^N)|   ≤   lim sup_N→∞2(K+1)√(N)(|1/NR_max-Λ|C/Λ+|M/N-γ|) 
       <   ∞.





By Proposition <ref>, √(N)|β(Y_s^N)-b(Y_s^N)|=O(1), then

    |D_t^N|   ≤ |D_0^N|+√(N)|M_t^N|+O(1)t+∫_0^t√(N) |b(Y_s^N)-b(y_s)|ds
       ≤|D_0^N|+√(N)|M_t^N|+O(1)t+∫_0^t√(N)L|Y_s^N-y_s|ds
       =|D_0^N|+√(N)|M_t^N|+O(1)t+∫_0^tL|D_s^N|ds.

By Gronwall's Lemma,

    sup_0≤ t≤ T|D_t^N|≤ e^LT(|D_0^N|+O(1)T+sup_0≤ t≤ T|√(N)M_t^N|).

Then

    lim sup_N→∞𝔼(sup_0≤ t ≤ T|D_t^N|^2)≤ e^2LT[lim sup_N→∞𝔼(|D_0^N|)+O(1)T+lim sup_N→∞𝔼(sup_0≤ t≤ T√(N)|M_t^N|)]^2.

By Jensen's inequality and Proposition <ref>, we have that 

    [𝔼(sup_0≤ t≤ T√(N)|M_t^N|)]^2≤ N𝔼(sup_0≤ t≤ T|M_t^N|^2)≤ 4NA(N)T,

and that A(N)=O(1/N). Therefore

    lim sup_N→∞𝔼(sup_0≤ t≤ T√(N)|M_t^N| )<∞.

Together with our assumption lim sup_N→∞𝔼(|D_0^N|^2)<∞, we have

    lim sup_N→∞𝔼(sup_0≤ t ≤ T|D_t^N|^2 )<∞.




To prove the tightness of (D^N)_N=1^∞ and the continuity of the limit points, we only need to show that the following two tightness conditions hold for each T>0 and ϵ>0,


  (i) 

    lim_K→∞lim sup_N→∞ℙ(sup_0≤ t≤ T|D_t^N|>K )=0,


  (ii) 

    lim_δ→ 0lim sup_N→∞ℙ(w(D^N,δ,T)≥ϵ)=0


where for x∈𝔻^d,

    w(x,δ,T)=sup{sup_u,v∈[t,t+δ]|x(u)-x(v)|:0≤ t≤ t+δ≤ T}.

By Lemma <ref>, there exists C_0>0 such that

    lim_K→∞lim sup_N→∞ℙ(sup_0≤ t≤ T|D_t^N|>K )    ≤   lim_K→∞lim sup_N→∞𝔼(sup_0≤ t≤ T|D_t^N|^2 )/K^2
       ≤   lim_K→∞C_0/K^2
       =   0,

which proves condition (i).

For condition (ii), we have that


    D^N_u - D^N_v   =   √(N)· ( M^N_u - M^N_v)_first term + ∫^u_v√(N)( β(Y^N_z) - b(Y^N_z)  ) dz_second term
       +   ∫^u_v√(N)( b(Y^N_z) - b(y_z)  ) dz _third term

for any 0<t≤ u<v≤ t+δ≤ T.  Now it suffices to show that each of the three terms of D^N_u - D^N_v satisfies condition (ii).  In what follows, we will show that each of the three terms satisfies condition (ii) to complete the proof of tightness. 

 For the first term √(N)· ( M^N_u - M^N_v), we would like to show that the limiting sample path of √(N)M_t^N is a continuous Brownian motion,  by using the martingale central limit theorem.
 
Similar to the proof of Proposition <ref>, we can show that 

    sup_t≤ T|β_+(Y_t^N)-b_+(Y_t^N)|0,   sup_t≤ T|β_-(Y_t^N)-b_-(Y_t^N)|0.

And by the proof of Proposition <ref>, b_+(y), b_-(y) are also Lipschitz with constant L, then by the fact that the composition of Lipschitz functions are also Lipschitz,

    max{sup_t≤ T|b_+(Y_t^N)-b_+(y_t)|,sup_t≤ T|b_-(Y_t^N)-b_-(y_t)|}≤ Lsup_t≤ T|Y_t^N-y_t|.

By Theorem <ref>, 

    sup_t≤ T|Y_t^N-y_t| 0.

Thus combining Equations (<ref>), (<ref>) and (<ref>), we have 

    lim_N→∞ℙ(sup_t≤ T|√(N)M^N(k)_t-  M(k)_t|>ϵ)
       =   lim_N→∞ℙ(sup_t≤ T|∫_0^t(β_+(Y_s^N)(k)+β_-(Y_s^N)(k)- b_+(y_s)(k)-b_-(y_s)(k))ds|>ϵ)
       ≤   lim_N→∞ℙ(sup_t≤ TT|β_+(Y_t^N)(k)- b_+(Y^N_t)(k)|>ϵ/4)+lim_N→∞ℙ(sup_t≤ TT|b_+(Y_t^N)(k)- b_+(y_t)(k)|>ϵ/4) 
       +   lim_N→∞ℙ(sup_t≤ TT|β_-(Y_t^N)(k)- b_-(Y^N_t)(k)|>ϵ/4)+lim_N→∞ℙ(sup_t≤ TT|b_-(Y_t^N)(k)- b_-(y_t)(k)|>ϵ/4)
       ≤   lim_N→∞ℙ(sup_t≤ TT|β_+(Y_t^N)(k)- b_+(Y^N_t)(k)|>ϵ/4)+2lim_N→∞ℙ(sup_t≤ TLT|Y_t^N- y_t|>ϵ/4 ) 
        +   lim_N→∞ℙ(sup_t≤ TT|β_-(Y_t^N)(k)- b_-(Y^N_t)(k)|>ϵ/4)
       =    0,

for any ϵ>0 and 0≤ k≤ K. This result implies that

    sup_t≤ T|√(N)M^N(k)_t-  M(k)_t| 0.


For the adjacent terms, we have

    lim_N→∞ℙ(sup_t≤ T|√(N)M^N(k),√(N)M^N(k+1)_t-  M(k),M(k+1)_t|>ϵ)
       =   lim_N→∞ℙ(sup_t≤ T|∫_0^t[∑_r(1/rNR_maxY_s^N(r,k+1)+(M/N-∑_n'∑_r'n'Y_s^N(r',n'))Y_s^N(r,k)). ..
       -   ...(∫_0^1Λ/rdy_s(r,k+1)+∫_0^1(γ-∑_n∫_0^1ndy_s(r,n)) dy_s(r,k))]ds|>ϵ)
       ≤   lim_N→∞ℙ(sup_t≤ TT|∑_r1/rNR_maxY_t^N(r,k+1)-∫_0^1Λ/rdY_t(r,k+1) |>ϵ/3)
       +   lim_N→∞ℙ(sup_t≤ TT|(M/N-∑_n',r'n'Y_t^N(r',n'))Y_t^N(r,k)- ∫_0^1(γ-∑_n∫_0^1ndY_t(r,n)) dY_t(r,k)|>ϵ/3)
       +   lim_N→∞ℙ(sup_t≤ T2LT|Y_t^N- y_t|>ϵ/3 )
       ≤   lim_N→∞ℙ(sup_t≤ TT|C/Λ NR_max-Λ/Λ /C|∑_rY_s(r,k+1) >ϵ/3)
       +   lim_N→∞ℙ(sup_t≤ TT|M/N-γ| ∑_rY_s(r,k)>ϵ/3)
       ≤   lim_N→∞ℙ(sup_t≤ TCT/Λ|1/ NR_max-Λ| >ϵ/3)+lim_N→∞ℙ(sup_t≤ TT|M/N-γ| >ϵ/3)
       =   0,

which implies

    sup_t≤ T|√(N)M^N(k),√(N)M^N(k+1)_t-  M(k), M(k+1)_t| 0.

Since for all |i-j|>1, 

    √(N)M^N(i),√(N)M^N(j)_t = M(i), M(j)_t=0.

We can conclude that 

    sup_t≤ T|√(N)M^N(i),√(N)M^N(j)_t-  M(i), M(j)_t| 0.

for all 0≤ i,j≤ K.

We also know that the jump size of Y^N_t is 1/N, therefore

    lim_N→∞𝔼[sup_0<t≤ T|√(N)M^N_t-√(N)M^N_t-| ]=lim_N→∞𝔼[sup_0<t≤ T|√(N)Y^N_t-√(N)Y^N_t-| ]=0.

By Theorem 1.4 in Chapter 7 of <cit.>,  √(N)M^N_t converges to the continuous Brownian motion M_t in distribution in 𝔻(ℝ_+,ℝ^K+1). By Prohorov's theorem, (√(N)M^N)_N=1^∞ is tight. This automatically implies the tightness condition ii).
Thus

    lim_δ→ 0lim sup_N→∞ℙ(sup_u,v∈ [t,t+δ]| √(N)(M_u^N-M_v^N)|>ϵ)
       ≤   lim_δ→ 0lim sup_N→∞ℙ(sup_u,v∈ [t,t+δ]| √(N)M_u^N-M_u|>ϵ/3)+lim_δ→ 0lim sup_N→∞ℙ(sup_u,v∈ [t,t+δ]|√(N)M_v^N-M_v|>ϵ/3)
          +lim_δ→ 0lim sup_N→∞ℙ(sup_u,v∈ [t,t+δ]| M_u-M_v|>ϵ/3)
       ≤   lim_δ→ 0ℙ(sup_u,v∈ [t,t+δ]| M_u-M_v|>ϵ/3)
       =    0



For the second term ∫^u_v√(N)( β(Y^N_z) - b(Y^N_z)  ) dz, we have by Proposition <ref> that the quantity √(N)( β(Y^N_z) - b(Y^N_z)  ) is bounded for any value of z∈ [0,T].  Therefore, there exists some constant C_1 that does not depend on N such that

    sup_z∈ [0,T]√(N)| β(Y^N_z) - b(Y^N_z)  |≤ C_1.

Then

    lim_δ→ 0lim_N→∞ℙ(sup_u,v∈ [0,T],|u-v|≤δ∫^u_v√(N)| β(Y^N_z) - b(Y^N_z)  | dz > ϵ)
       ≤   lim_δ→ 0lim_N→∞ℙ(δsup_z∈ [0,T]√(N)| β(Y^N_z) - b(Y^N_z)  |  > ϵ)
       ≤   lim_δ→ 0ℙ(δ C_1  > ϵ)
       =    0.

 Thus, we have proved the oscillation bound for the second term.  
 
 Finally for the third term we have that 

    ∫^u_v√(N)| b(Y^N_z) - b(y_z)  | dz    ≤   ∫^u_v√(N) L| Y^N_z - y_z| dz
       =   ∫^u_v L ·|D^N_z| dz 
       ≤    Lδsup_t∈ [0,T]|D^N_t|.

By Lemma <ref>,  

    lim_δ→ 0lim_N→∞ℙ(sup_u,v∈ [0,T],|u-v|≤δ∫^u_v√(N)| b(Y^N_z) - b(y_z)  | dz>ϵ)
       ≤   lim_δ→ 0lim_N→∞ℙ(Lδsup_t∈ [0,T]|D^N_t|>ϵ)
       ≤   lim_δ→ 0lim_N→∞𝔼(sup_t∈ [0,T]|D^N_t|^2)/(ϵ/L δ)^2
       ≤   lim_δ→ 0C_0(Lδ)^2/ϵ^2
       =    0,

 which implies that the oscillation bound holds for the third term.  






To prove the existence and uniqueness of the SDE (<ref>), we show the following two conditions hold:
There exists a constant H>0 such that


  (1)Lipschitz condition: for any D,D̃∈ℝ^K+1, any t∈ (t_0,T), 

    |b'(y_t)D-b'(y_t)D̃|≤ H|D-D̃|.


  (2)Linear growth condition: for any D∈ℝ^K+1, any t∈ (t_0,T),

    |b_+(y_t)+b_-(y_t)|≤ H(1+|D|),    |b'(y_t)D|≤ H(1+|D|).



By Proposition <ref>, b(y) is Lipschitz, and by our assumption in Theorem <ref>, b(y) is also continuously differentiable, thus b'(y) is bounded. Then conditions (1)  and (2) follow, which prove that there exist a unique solution to the SDE (<ref>).

Take expectation on both sides of Equation (<ref>), since 𝔼[ ∫_0^te^∫_s^tb'(y_u)dudM_s]=0, we have 
    𝔼[D_t]=e^∫_0^t𝒜(s)ds𝔼[D_0].

Therefore

    D_t-𝔼[D_t]   =   e^∫_0^t𝒜(s)ds(D_0-𝔼[D_0])+∫_0^te^∫_s^t𝒜(u)dudM_s,


    Σ(t)   =   E[(D_t-𝔼[D_t])(D_t-𝔼[D_t])^⊤]
       =   e^∫_0^t𝒜(s)dsE[(D_0-𝔼[D_0])(D_0-𝔼[D_0])^⊤] (e^∫_0^t𝒜(s)ds)^⊤
       +   (∫_0^te^∫_s^t𝒜(u)dudM_s)(∫_0^te^∫_s^t𝒜(u)dudM_s)^⊤
       =   e^∫_0^t𝒜(s)dsΣ(0)e^∫_0^t𝒜^⊤(s)ds+∫_0^te^∫_s^t𝒜(u)duℬ(s)e^∫_s^t𝒜^⊤(u)duds.



<ref>
For any s≥ 0,

    lim sup_N→∞√(N)|β(Y_s^N)-b(Y_s^N)|<∞.



 
    |β(Y_s)-b(Y_s)|
    ≤   2∑_n|∑_r,p,kpP^N_max/rR^N_maxY_s(r,n,p,k)-∑_k∈𝒦∬_[0,1]× [0,1]pΛ/rdY_s(r,n,p,k)|
       +2∑_n|∑_r,p,kpNP^N_maxM/NY_s(r,n,p,k)-∑_k∈𝒦∬_[0,1]× [0,1]p𝒫γ dY_s(r,n,p,k)|
       +2∑_n|∑_r,p,k(∑_n'∑_r',p',k'n'Y_s(r',n',p',k'))Y_s(r,n,p,k).
       -.∑_k∈𝒦∬_[0,1]× [0,1](∑_n'∑_k∈𝒦∬_[0,1]× [0,1]n'dY_s(r',n',p',k'))dY_s(r,n,p,k)|.
 
 By Equation (<ref>) and (<ref>), 
 
    |β(Y_s)-b(Y_s)|≤ 2(K_max+1)(|P_max/R_max-Λ|C/Λ+|NP_maxM/N-𝒫γ|).

 By the assumptions in Theorem <ref>, we have
 
    lim sup_N→∞√(N)(P_max^N/R_max^N-Λ)<∞, lim sup_N→∞√(N)(M/N-γ)< ∞, lim sup_N→∞√(N)(NP_max^N-𝒫)< ∞.

Thus 

    lim sup_N→∞√(N)|β(Y_s^N)-b(Y_s^N)|   ≤   lim sup_N→∞2(K_max+1)√(N)(|P_max/R_max-Λ|C/Λ+|NP_maxM/N-𝒫γ|) 
       <   ∞.





<ref>
For k<K_max,

    lim_N→∞ℙ(sup_t≤ T|√(N)M^N(k),√(N)M^N(k+1)_t-  M(k),M(k+1)_t|>ϵ)=0.


For the adjacent terms, we have

    lim_N→∞ℙ(sup_t≤ T|√(N)M^N(k),√(N)M^N(k+1)_t-  M(k),M(k+1)_t|>ϵ)
       =   lim_N→∞ℙ(sup_t≤ T|∫_0^t[∑_r,p,K∈𝒦(pP_max/rR_maxY_s^N(r,k+1,p,K)+(pNP_maxM/N-∑_n'∑_r',p',K'n'Y_s^N(r',n',p',K')). ...
          .. . . Y_s^N(r,k,p,K))-(∑_K∈𝒦∬_[0,1]× [0,1]pΛ/rdY_s(r,k+1,p,K). ...
       +    .. . . ∑_K∈𝒦∬_]0,1]× [0,1]p𝒫(γ-∑_n∑_K∈𝒦∬_]0,1]× [0,1]ndy_s(r,n,p,K)) dy_s(r,k,p,K))]ds|>ϵ)
       ≤   lim_N→∞ℙ(sup_t≤ TT|∑_r,p,KpP_max/rR_maxY_t^N(r,k+1,P,K)-∑_K∈𝒦∬_[0,1]× [0,1]pΛ/rdY_s(r,k+1,p,K) |>ϵ/3)
       +   lim_N→∞ℙ(sup_t≤ TT|∑_r,p,K∈𝒦(pNP_maxM/N-∑_n'∑_r',p',K'n'Y_s^N(r',n',p',K'))Y_t^N(r,k,p,K)..
          .. - ∑_K∈𝒦∬_]0,1]× [0,1]p𝒫(γ-∑_n∑_K∈𝒦∬_]0,1]× [0,1]ndy_s(r,n,p,K)) dy_s(r,k,p,K)|>ϵ/3)
       +   lim_N→∞ℙ(sup_t≤ T2LT|Y_t^N- y_t|>ϵ/3 )
       ≤   lim_N→∞ℙ(sup_t≤ TT|CP_max/Λ R_max-Λ/Λ /C|∑_r,p,KY_s(r,k+1,p,K) >ϵ/3)
       +   lim_N→∞ℙ(sup_t≤ TT| NP_maxM/N-𝒫γ| ∑_r,p,KpY_s(r,k,p,K)>ϵ/3)
       ≤   lim_N→∞ℙ(sup_t≤ TCT/Λ|P_max/ R_max-Λ| >ϵ/3)+lim_N→∞ℙ(sup_t≤ TT|NP_maxM/N-𝒫γ| >ϵ/3)
       =   0,

which implies

    sup_t≤ T|√(N)M^N(k),√(N)M^N(k+1)_t-  M(k), M(k+1)_t| 0.





<ref>


b(y) is continously differentiable with the derivatives 
∂/∂ y(j)b(y)(i) as follows,

    ∂/∂ y(j)b(y)(0)   =   j· y(0)+1/y(j)∫_0^1Λ/rdy(r,1)1_{j=1}-(γ-∑_n=0^Kny(n))1_{j=0},
    ∂/∂ y(j)b(y)(k)   =   j·( y(k)-y(k-1))+1/y(j)(∫_0^1Λ/rdy(r,k+1)1_{j=k+1}-∫_0^1Λ/rdy(r,k)1_{j=k})
          +(γ-∑_n=0^Kny(n))(1_{j=k-1}-1_{j=k})  for  0<k<K,
    ∂/∂ y(j)b(y)(K)   =   -j· y(K-1)-1/y(j)∫_0^1Λ/rdy(r,K)1_{j=K}+(γ-∑_n=0^Kny(n))1_{j=K-1}.




Since y_t(k)=∫_0^1dy_t(r,k) for k=0,⋯, K, we have

    b(y_t)(k)=   ∫_0^1Λ/rdy_t(r,k+1)1_k<K_retrieve a bike from a k+1-bike station-∫_0^1(Λ/r1_k>0+(γ-∑_n∫_0^1 ndy_t(r,n) )1_k<K)dy_t(r,k)_retrieve and return a bike to a k-bike station
       +∫_0^1(γ-∑_n∫_0^1 ndy_t(r,n) )dy_t(r,k-1)1_k>0_return a bike to a k-1-bike station
    
    =    y_t(k+1)1_k<K∫_0^1Λ/rdy_t(r,k+1)/y_t(k+1)-y_t(k)1_k>0∫_0^1Λ/rdy_t(r,k)/y_t(k)
       +(γ-∑_n ny_t(n) )(y_t(k-1)1_k>0-y_t(k)1_k<K).

By condition i) in Assumption <ref>, 

    1/N∑_i=1^N1 _(r_i^N,K_i^N)⇒ I(r,k).
 As mentioned in the begining of Section <ref>,  we assume without loss of generality that K_i^N=K. This is valid because when {K_i^N}_i  are different, we can still carry out the proof of the mean field and diffusion limits by simply replacing r with (r,K). Thus, in this case the limiting measure I(r,k)=I(r,K) is one dimensional.

By our definition, r_i^N, the relative utilization at station i,  only denpends on the station number i. It does not depend on the number of bikes at that station. This means that the two components of the empirical measure y(r,n), the relative utilization r and the number of bikes n,  are independent. Then

    ∫_0^1Λ/rdy_t(r,k)/y_t(k)=∫_0^1Λ/rdI(r,K),

which no longer depends on y_t(k). Therefore we have the derivatives of b(y_t)(k) as follows,

    ∂/∂ y_t(j)b(y_t)(k)   =   j·( y_t(k)1_k<K-y_t(k-1)1_k>0)+1/y_t(j)(∫_0^1Λ/rdy_t(r,k+1)1_{j=k+1}.
          .-∫_0^1Λ/rdy_t(r,k)1_{j=k>0})+(γ-∑_n=0^Kny_t(n))(1_{j=k-1}-1_{j=k<K}) ,

for 0≤ k,j ≤ K. 

We can also write the derivatives of b(y) seperately as follows,

    ∂/∂ y(j)b(y)(0)   =   j· y(0)+1/y(j)∫_0^1Λ/rdy(r,1)1_{j=1}-(γ-∑_n=0^Kny(n))1_{j=0},
    ∂/∂ y(j)b(y)(k)   =   j·( y(k)-y(k-1))+1/y(j)(∫_0^1Λ/rdy(r,k+1)1_{j=k+1}-∫_0^1Λ/rdy(r,k)1_{j=k})
          +(γ-∑_n=0^Kny(n))(1_{j=k-1}-1_{j=k})  for  0<k<K,
    ∂/∂ y(j)b(y)(K)   =   -j· y(K-1)-1/y(j)∫_0^1Λ/rdy(r,K)1_{j=K}+(γ-∑_n=0^Kny(n))1_{j=K-1}.

We can see that the derivatives of b(y) are linear in y=(y(0),⋯,y(K)). Thus we can conclude that b(y) is continuously differentiable with respect to y.




plainnat

