# Dynamic Capacitated Arc Routing Problem in E-Bike Sharing System: A Monte Carlo Tree Search Approach

Authors: Shiqi Tan; Zhiheng Li; Na Xie
Journal: Journal of Advanced Transportation
 (2021)

---

## Abstract

This paper studies a dynamic capacitated arc routing problem for battery replacement in an e-bike sharing system, where workers replace batteries for underpowered e-bikes along street segments dynamically. The objective is to replace as many batteries as possible and minimize pickup failures. The temporal dependency of the routing decisions, the conflict of the workers’ operations, and the stochastic and dynamic nature of user demands all make this a difficult problem. To cope with these difficulties, a “Partition-First, Route-Second” bi-level solution framework is adopted to describe the problem in two different time scales. On the large time scale, a spatiotemporal partitioning method, which divides the road network into nonoverlapping subzones, is proposed to decompose the problem. On the small time scale, this paper models the routing decision process of individual worker as a Markov Decision Process. We adopt a lookahead policy that simulates future information and decisions over some horizons to evaluate the long-term influence of current feasible decisions. A Monte Carlo Tree Search algorithm is also used to improve the simulation efficiency. By performing numerical computation experiments on a test case study and comparing with some benchmarking policies, we demonstrate the effectiveness and efficiency of the suggested method.


---

## Body

## 1. Introduction

The “Mobile Internet” and “Sharing Economy” have greatly promoted the development of shared micromobility services, represented by free-floating bike sharing system, shared electric scooter system, and e-bike sharing system, which provide a more convenient and environmentally friendly “start/last mile” travel solution. The E-Bike Sharing System (EBSS) is a relatively new form of shared micromobility service, which could provide auxiliary power through the battery and help users reduce the physical energy consumption when they are riding against the wind, uphill, with loads, or over long distances and thus can significantly improve user experience. In China, the country with one of the world’s largest markets for electric bicycles, the shared micromobility service providers, such as Hellobike and Qingju Bike, have already launched EBSS in many cities, especially those with lots of slopes in the road network.Despite the improvement in user experience, EBSS faces more challenges in daily operation. A new challenge is how to replace the detachable batteries of underpowered e-bikes scattered throughout the city in a timely and cost-effective way. At present, the EBSS providers in China usually take the approach of setting up “virtual parking spots” as illustrated in Figure1. The virtual parking spots are areas where parking is allowed as determined by the built-in GPS positioning and are usually distributed along the street segments. Users who park e-bikes outside these areas will face punitive dispatch fees. In this way, it not only reduces the worker’s efforts to find an e-bike but also restrains the messy parking behavior.Figure 1
An illustration of the distribution of virtual parking spots. Users can easily query the surrounding and the nearest virtual parking spots through the EBSS platform.In the above context, since the underpowered e-bikes are mostly aggregated over the virtual parking spots along the street segments, the routing problem of battery replacement in EBSS can be modeled as a capacitated arc routing problem (CARP), just like the Chinese postman problem [1]. An optional way is to replace batteries at night when the user demands can be ignored, corresponding to a static version of CARP (SCARP). In the static version, the distribution of underpowered e-bikes is fixed and known in advance, and a set of minimum cost routes that service all underpowered e-bikes is expected to be found for a fleet of capacitated vehicles located at a given depot. However, there is a drawback with this passive static way: it cannot guarantee the number of e-bikes with sufficient power maintained at a high level throughout the next day, resulting in lower utilization rates of the system, especially during the evening peak.Therefore, this paper mainly considers a dynamic version of the CARP (DCARP) in EBSS: the battery replacement is performed during daytime, when the user’s participation must be considered and new underpowered e-bikes may stochastically appear after the fleet departs from depot. Although there is a large amount of studies related to the dynamic capacitated vehicle routing problem, which was surveyed by Pillac et al. [2] and Psaraftis et al. [3], direct studies on DCARP are relatively few. Some literatures [4–7] investigated the periodic CARP in certain scenarios and proposed various heuristic algorithms; some literatures [8–12] studied the dynamic case where the demand or deadheading cost of each arc is stochastic and proposed some repair operators or adaptive algorithms rather than just an initial static solution.However, the DCARP in the scenario of battery replacement in EBSS has the following difficulties: the first difficulty is the temporal dependency of the routing decisions. The current routing decisions for battery replacement will affect the success of subsequent pickup demands, which in turn affects the future state of e-bikes and the future routing decisions. This intricate temporal dependency results in a significant increase in the complexity of the problem. The second difficulty is the conflict of the workers’ operations. Since there is usually more than one worker, it is important to design a cooperation mechanism to avoid simultaneous operations on the same street segment by different workers and to improve the operation efficiency throughout the entire system. The third difficulty is how to anticipate the future influence of current decisions with stochastic and dynamic user demands. In the dynamic case, the goal is not only to maximize the number of batteries replaced but also to maximize the successful attempts to pick up e-bikes, and this requires anticipating potential future user demands when making current decisions.These difficulties make it hard to compute a fixed decision sequence in advance. To cope with these difficulties, this paper adopts a “Partition-First, Route-Second” bi-level solution framework to describe this problem in two different time scales, as shown in Figure2.Figure 2
The “Partition-First, Route-Second” bi-level solution framework.On the upper level, the road network is partitioned into subzones to decompose the problem into multiple smaller subproblems. This idea of partitioning streets into clusters has been widely used in the classical CARP and its variants [13–17]. In the setting of our problem, the road network partition provides constraints of feasible operation area and reward to each worker, preventing the potential conflict of workers’ operations. The partitioning also reduces the temporal dependence between different workers’ routing decisions to some extent, due to the short-distance nature of e-bike trips. To adapt the partition to the real-time state, a spatiotemporal partitioning method is proposed on a large time scale, where an initial partition is generated at the beginning of a day and will be adjusted dynamically with the real-time state of e-bikes based on some heuristic algorithms.On the lower level, the routing decision process of each worker is modeled as a Markov Decision Process (MDP) on a small time scale, which is small enough to capture the details of the user demands and worker operations. In the MDP, decisions should be made about which adjacent street segment to visit next when a worker has finished the task on the previous street segment. This paper adopts a lookahead policy with Monte Carlo Tree Search (MCTS) algorithm, which focuses on solving an approximation of the problem over some horizon to anticipate the future impacts [18, 19]. At a decision point of a worker, we model all feasible sequences of routing decisions for the worker over some horizon into a decision tree, then search and simulate it in the simulator based on the MCTS algorithm, which will allocate more simulation runs to more highly rewarding routes, and thus save computational resources in the simulation process.The main contributions of this paper are as follows. We model the routing problem of battery replacement in EBSS into a DCARP, based on the distribution of virtual parking spots. Then, we adopt a “Partition-First, Route-Second” solution framework to cope with the difficulties of this problem. We take both user satisfaction and operational efficiency into account and adopt a lookahead policy with MCTS algorithm to balance the relationship of exploration and exploitation during the simulation.The next parts of this paper are structured as follows. In Section2, we briefly review the related literature. In Section 3, we summarize the notations used in this paper and formulate the problem. In Section 4, a spatiotemporal partitioning method is proposed in detail. In Section 5, we introduce the lookahead policy with MCTS algorithm. In Section 6, we present the design of the system simulator. In Section 7, we analyze the results of numerical studies. Finally, we conclude the paper in Section 8.


## 2. Literature Review

As a convenient way for short-distance travel, bike sharing system has been expected to help mitigate traffic congestion, alleviate the growing environmental issues, and promote sustainable transportation development [20]. However, the development of the bike sharing systems has also undergone a long trial-and-error process. The first generation of bike sharing systems was dockless and operated as nonprofits but failed due to theft and damage [21]. Then, to facilitate maintenance, station-based bike sharing systems became the dominant form for a long time and went through stages such as coin-deposit systems, IT-based systems, and demand-responsive multimodal systems that made it easier for users to pick up or return bikes [22]. Besides, numerous studies [23–27] have investigated the bike allocation and reposition problems in station-based bike sharing systems to rebalance the system from nonhomogeneous asymmetric demand processes. However, the user base of the station-based bike sharing system has been limited, with the constraint that it is difficult to build dense stations due to the high cost and road space required.In recent years, with the support of emerging technologies such as mobile payment and electronic locking, the Internet-based free-floating bike sharing system has become the main mode, where users are allowed to pick up and return bikes at any location near their origin and destination. This increases the coverage of the system and attracts more users. However, with the rapid expansion of free-floating bike sharing system, a series of new problems have arisen, such as the recycling of faulty bikes [28] and the messy parking behaviors [29]. Therefore, the promotion of setting virtual parking spots in some areas is necessary [29]. The EBSS, which needs frequent replacement of batteries for underpowered e-bikes, has adopted the approach of setting up virtual parking spots in many cities in China. These virtual parking spots are quite different from the parking stations in traditional station-based bike sharing systems. The virtual parking spots do not require the construction of physical parking stations or docks and thus could be set up more densely throughout the city at a lower cost and in a more flexible way. However, to the best of our knowledge, there is no research on the problem of dynamic battery replacement in the EBSS with virtual parking spot setting so far.From another perspective, the dynamic battery replacement problem in EBSS is similar to the dynamic rebalancing problem in traditional station-based bike sharing systems in terms of user demand distribution and system operation mechanism. Thus, it would be helpful to review some of the dynamic rebalancing policies in station-based bike sharing systems. The rebalance problem is usually modeled as a pickup and delivery problem in the literature. Brinkmann et al. [23] proposed a short-term rebalancing policy, and according to this policy, the vehicle is sent to the nearest shared bike station where the safety buffer is breached. This policy is simple and computationally efficient but does not take future information into account. Ghosh et al. [24] simplified the relocation problem by aggregating the stations and trips and proposed a mixed-integer linear program on a rolling horizon. The high computational cost of this algorithm makes it difficult to meet the requirements of real-time applications. Brinkmann et al. [25] presented a dynamic lookahead policy to anticipate potential future demands over a predefined horizon by rollout simulation. However, this policy is less efficient when multiple vehicles should be considered and when the decision process occurs too frequently. Li et al. [26] used neural network to approximate the long-term value function of current decisions via a reinforcement learning approach. This deep reinforcement learning approach has the appeal of fast decision speed and the ability to consider long-term effects. However, it is still a challenging work to learn an accurate and stable value function neural network due to the stochastic exogenous information and complexity of the environment in EBSS.Taking the distribution of virtual parking spots into consideration, we model the battery replacement problem in EBSS as a dynamic capacitated arc routing problem. The capacitated arc routing problem is aimed to determine an optimal traversal of an arc subset of a graph by some capacitated vehicles [1] and has been widely studied in the past years, such as by converting to the capacitated vehicle routing problem [30, 31], or directly by exact optimization methods [32, 33], or by some lower bounding procedures [34–37] and heuristic methods [13, 38–41]. Among them, the “Partition-First, Route-Second” heuristic has been widely used in CARP and its variants and has been proved to be effective [13–17]. This paper adopts this heuristic idea to decompose the original problem into some smaller subproblems by partitioning the arcs into some subzones.The dynamic version of CARP considers the stochastic and dynamic demand or deadheading cost, which greatly increases the problem complexity. Christiansen et al. [8] considered a variant of the CARP with stochastic demands on each arc and proposed a branch-and-price algorithm. They incorporated the stochastic demand into the pricing problem. Laporte et al. [9] proposed an adaptive large-scale neighborhood search heuristic method for the CARP with stochastic demands. Tagmouti et al. [10, 11] studied the CARP with time-dependent service costs and proposed a variable neighborhood descent heuristic method. Mei [12] aimed to find a robust solution for all possible environments instead of the optimal solution for a specific environment for the dynamic CARP and designed a repair operator which could repair the solution to feasible dynamically. However, these methods may fail to anticipate the difficulties of the DCARP in EBSS, such as the temporal dependency of the routing decisions and the influence of the future demand trend. Thus, this paper proposes a lookahead policy based on MCTS algorithm for the DCARP in EBSS, which is an efficient algorithm for sequential decision problems [19].


## 3. Notation and Problem Formulation

In this section, we summarize the notations used in this paper in Section3.1, formulate the DCARP in Section 3.2, and introduce the decomposition of the decision process of multiple workers and model the routing decision process of an individual worker as MDP in Section 3.3.


### 3.1. Notation

The main notations and their values used in this paper are summarized in Table1.Table 1
Notations.
NotationMeaningModelNvNumber of intersections of the road networkviIntersection of the road network, fori∈1,2,…,Nvv0The depotVSet of vertexes, including the intersections and depot,V=v0,v1,v2,…vNvvi,vjStreet segment (edge) of the road networkvi,vj∈V,ESet of edges,E=vi,vj|vi,vj∈Vandi<jvi,vjDirected street segment (arc),vi,vj∈VASet of arcs,A=vi,vj|vi,vj∈Vandi≠jGGraph of the road network,G=V,E in the undirected representation or G=V,A in the directed representationTNumber of time steps in an operating day (1440)tTime step on the small time scale,t∈1,2,…,TΔtInterval between time steps (1 min)tbeginningThe beginning of the working hours (7:00 am)tendThe end of the working hours (10:00 pm)ψmaxMaximum range of e-bike battery (40 km)ψminMinimum battery range threshold for an e-bike to be available (3 km)τbServing time cost of replacing battery for an e-bike (1 minute)NebNumber of e-bikes in system (1000)NwNumber of workers/vehicles in system (4)wiWorker,i∈1,2,…,NwWSet of workers,W=w1,w2,…wNwBBattery capacity of each vehicle (30)cijDeadheading time cost of arcvi,vjτdTime cost of reloading a battery in the depot (0.5 minutes)ψ¯userExpectation of battery range consumption for e-bike trips (2.75 km)btijNumber of batteries replaced on arcvi,vj at time step t, t∈1,2,…,TdtijNumber of successful pickups on arcvi,vj at time step t, t∈1,2,…,TGiSubzone of the road network,i∈1,2,…,NwViVertex set in subzoneGi, i∈1,2,…,NwAiArc set in subzoneGi, i∈1,2,…,NwEiEdge set in subzoneGi, i∈1,2,…,NwpTime period on the large time scaleΔpTime interval between time periods (180 minutes)KmThe decision point set of workerwm, m∈1,2,…,NwtkmTime step for thekth decision point of worker wm, tkm∈1,2,…,T, m∈1,2,…,Nw, k∈KmηkmA vector of lengthA indicating the available e-bike number on the arcs in the subzone of worker wm at the kth decision point (the values of arcs in other subzones are masked to 0), m∈1,2,…,Nw, k∈KmφkmA vector of lengthA indicating the underpowered e-bike number on the arcs in the subzone of worker wm at the kth decision point (the values of arcs in other subzones are masked to 0), m∈1,2,…,Nw, k∈KmLkmA one-hot vector of lengthA indicating which arc worker wm is located on at the kth decision point, m∈1,2,…,Nw, k∈KmCkmA nonnegative number indicating the remaining battery number in the vehicle of workerwm at the kth decision point, m∈1,2,…,Nw, k∈KmSmState set of workerwm, m∈1,2,…,NwskmState of workerwm at the kth decision point, skm∈Sm, and skm=tkm,ηkm,φkm,Lkm,Ckm, m∈1,2,…,Nw, k∈KmXskmFeasible decision set under stateskmxkmDecision of workerwm at the kth decision point, xkm is a one-hot vector of length A indicating the arc to visit next, xkm∈Xskm, m∈1,2,…,Nw, k∈KmrkmInstant reward for workerwm at the kth decision point, m∈1,2,…,Nw, k∈KmπmPolicy for workerwm, m∈1,2,…,Nwπm∗Optimal policy for workerwm, m∈1,2,…,NwMethodmijThe link relationship of theith and the jth edge in E, i,j∈1,2,…,EMAdjacency matrix of edge set,M=mijE×EeijA 0-1 variable to indicate whether thejth edge in E belongs to the subset Ei, i∈1,2,…,Nw, j∈1,2,…,EeiA 0-1 vector to indicate the edges in subsetEi, ei=eij∀j∈1,2,…,E, i∈1,2,…,NwupjThe number of underpowered e-bikes on thejth edge at time period p, j∈1,2,…,EupThe vector of the number of underpowered e-bikes on each edge at time periodp, up=upj∀j∈1,2,…,EliThe length of theith edge (km), i∈1,2,…,ElThe vector of the length of each edge,l=li∀i∈1,2,…,EβTemperature in the simulated annealing algorithmβ0Initial temperature in the simulated annealing algorithm (100)αThe decay factor in temperature in the simulated annealing algorithm (0.9)MaxIterinnerInner iteration limit in the simulated annealing algorithm (10)MaxIterouterOuter iteration limit in the simulated annealing algorithm (20)IrSimulation runs for each feasible decision in the rollout policyImaxMaximum iteration number in the MCTS algorithmHLookahead simulation horizons (180 time steps)γA weighting factor in the UCB algorithm (5)λA time discount parameter (0.95)


### 3.2. Problem Formulation

Figure3 shows a simple scenario of EBSS, which consists of the road network, e-bikes scattered among the virtual parking spots or used on the road, and vehicles/workers replacing batteries for underpowered e-bikes.Figure 3
A sample scenario of EBSS.The road network is abstracted as a connected graphG=V,E, where V=v0,v1,v2,…vNv is the vertex set and E=vi,vj|vi,vj∈Vandi<j is the edge set. In the vertex set, v0 is the depot and vi|i=1,2,…,Nv is the set of intersections. Assume that each edge has two directions and the vehicle can only serve one direction at a time. Thus, an edge vi,vj can be split into two directed arcs: vi,vj and vj,vi, as shown by the edge v1,v2 in Figure 3. Then, the road network can also be represented as a directed graph: G=V,A, where A=vi,vj|vi,vj∈Vandi≠j is the arc set. The undirected graph representation is used in the road network partitioning and the directed graph representation is used in the routing decision process.The operating day is divided intoT discrete time steps: t∈1,2,…,T, according to a predefined time interval Δt (1 minute). The dynamic user demand is reflected in the stochastically generated pickup and return demands in each virtual parking spot at each time step, which leads to dynamic changes in the location and battery power level of e-bikes.The number of e-bikes in the system is supposed to be fixed atNeb, without new e-bikes being put into the system or removed from the system. The battery range of each e-bike is in the range of 0,ψmax (km). When the battery range is below a threshold ψmin (km), an e-bike is regarded as underpowered and unavailable to users. For the sake of simplicity, we assume that the time cost for a worker to replace a battery is a constant: τb (e.g., 1 minute).Assume there areNw workers (or vehicles) in the system, denoted by the set W=w1,w2,…,wNw. Each vehicle has the same battery capacity of B (e.g., 40). Suppose tbegin (e.g., 7:00 am) and tend (e.g., 10:00 pm) are the beginning and the end of working hours. During working hours, workers will continue to replace batteries for underpowered e-bikes. Once the available batteries are used up, workers need to return to the depot and reload batteries. Suppose vehicles travel with a certain speed that is independent of time, leading to a positive deadheading cost cij of each arc vi,vj∈A. Assume that the time cost of reloading a battery in the depot is τd (e.g., 0.5 minutes) and there are always enough fully charged batteries in the depot, so the time cost of reloading a vehicle’s battery in the depot is τd⋅B.Under the limited working hours and time cost constraints introduced above, on the one hand, we hope to route the vehicles in a way that replaces as many batteries as possible during a limited time; on the other hand, we also hope to prioritize services for those arcs that have the potential to reduce user pickup failures in the following time. Take workerw1 in Figure 3 for example; if there will be a user picking up e-bike on arc v4,v7 in the following time, it would be better for w1 to serve arc v4,v7 first instead of arc v6,v7. However, directly maximizing the sum of the number of replaced batteries and the number of successful pickups may be not appropriate. In fact, replacing one more battery has a higher expected reward than reducing one more pickup failure because a fully charged e-bike can usually be used multiple times and can reduce more potential pickup failures in the following days. For this consideration, the optimization objective function is designed as a weighted sum of these two factors, as shown in the following equation:(1)max∑t∈1,2,…,T∑vi,vj∈Aψmax−ψmin⋅btij+ψ¯user⋅dtij,where btij is the number of batteries replaced on arc vi,vj at time step t, dtij is the number of successful pickups on arc vi,vj at time step t, and ψ¯user indicates the expectation of battery power consumption for each e-bike trip. Therefore, ∑t∈1,2,…,T∑vi,vj∈Aψ¯user⋅dtij denotes the battery power consumed by the users, indicating the short-term rewards during current day, while ∑t∈1,2,…,T∑vi,vj∈Aψmax−ψmin⋅btij denotes the battery power replenished by the workers, indicating the long-term rewards of the system sustainability.


### 3.3. Decomposition of the Problem

For the above dynamic capacitated arc routing problem, it is hard to compute the optimal route for each worker in advance due to the stochastic and complex nature of the system. Therefore, we seek to formalize the problem as a sequential decision process in cooperative multiagent system. A decision point occurs when a worker completes the task on the current arc and a decision about which one of the adjacent arcs to visit next should be made. It should be noted that as the number of tasks and deadheading cost varies for each arc, the executing duration of each decision may be different and the decision processes between different workers are thus asynchronous. At each time step, only a subset of the workers may make a decision while others keep on executing the tasks on the current arcs.To simplify the problem, we assume that each worker is assigned a nonoverlapping subzone and is allowed to perform the battery replacement operation only within his subzone. As there areNw workers, the road network G=V,A will be partitioned into Nw subzones: G1,G2,…,GNw, satisfying the condition that Gi=Vi,Ai is strongly connected graph for i=1,2,…,Nw and A1∪A2∪…∪ANw=A, Ai∪Aj=0 for i≠j. In this case, the routing decisions of each worker have very limited effects on other workers, and we can decompose the cooperative multiagent decision problem in the whole road network into multiple mono-agent decision problems in each subzone approximatively.As shown in Figure4, the arc set A is partitioned into two subsets A1 and A2 for workers w1 and w2, and the decision process of each worker is modeled as a Markov Decision Process. The operating day is divided into discrete time steps t according to a small time interval Δt (e.g., 1 minute). At each time step, a worker can make a routing decision or continue with his current task. For example, the decision points for worker w1 in Figure 4 are time steps t+2, t+4, and t+6. We also divide the day into some discrete time periods p according to a larger time interval Δp (e.g., 180 minutes). At each time period (e.g., at p1 and p2 in Figure 4), the partition can be adjusted according to the real-time e-bike state.Figure 4
An illustration of the decomposition of decision processes.To model the routing decision process of each worker as a MDP, the following assumptions are made:(i)
After a worker leaves the depot, he follows the shortest path (i.e., calculated by the Dijkstra algorithm) to the nearest vertex of his subzone. During this process, the worker does not replace battery. Once the remaining batteries in the vehicle run out, the worker follows the shortest path back to the depot to reload the batteries and then goes to the previous location and continues his operation.(ii)
As new underpowered e-bikes may be generated on the arcs that have been served, multiple visits and multiple services to an arc are allowed. However, we suppose that temporary stops of vehicles are not allowed and that the vehicle’s power is always sufficient during working hours.Based on the above assumptions, the MDP for themth worker wm consists of the following components:A state:Sm is the feasible state set of worker wm. skm∈Sm is state of the kth decision point of worker wm. The state skm=tkm,ηkm,φkm,Lkm,Ckmincludes three main parts: the system time, the state of e-bikes, and the state of workers. The time step of the kth decision point of worker wm is denoted as tkm; the state of e-bikes includes ηkm and φkm, indicating the available e-bike number and the underpowered e-bike number on the arcs in the subzone of worker wm at the kth decision point, respectively. The state of worker includes the locations Lkm and the remaining battery number Ckm of worker wm at the kth decision point. A detailed description of the state information is given in Table 1.A decision:Xskm denotes the feasible decision set under state skm. Xskm is the intersection of the set of arcs adjacent to the current arc of wm and the set of arcs in the subzone of wm. For example, in Figure 4, v4,v0,v4,v1,v4,v7 is the set of arcs adjacent to v7,v4. However, v4,v1 belongs to different subzones from v7,v4. Therefore, when a worker is making a routing decision on arc v7,v4, the feasible decision set is v4,v0,v4,v7. The kth decision of worker wm is denoted as xkm∈Xskm, representing the next arc to visit.State transition function:Pmsk+1m|skm,xkm is the transition probability function for worker wm and describes the probability that the system will move to the next decision state sk+1m after making a decision xkm under the current decision state skm. The next decision state sk+1m will also be influenced by the stochastic exogenous information in the system (e.g., the stochastic user demand) and the operations of other workers.Reward function:rkm=Rmskm,xkm is the instant reward for worker wm at the kth decision point. Rmskm,xkm denotes the reward function for worker wm, reflecting the instant reward received after taking an decision xkm under state skm. rkm can be derived from the weighted sum of the number of successful pickups and the number of batteries replaced within the subzone of wm between time step tkm and tk+1m (if k is the last decision point of wm, then tk+1m is the last time step of the day), as shown in the following equation:(2)rkm=∑t=tkmtk+1m−1∑vi,vj∈Amψmax−ψmin⋅btij+ψ¯user⋅dtij.Since the subzones do not overlap, the cumulative reward in a subzone is largely related to the decisions of the worker in that subzone and less related to the decisions of other workers, and the sum of the cumulative rewards of all subzones is the total system reward, which is consistent with the optimization objective in equation (1).Policy:πm describes the policy to make a decision under current state for worker wm: xkm=πmskm. Suppose that s0m is the initial state of worker wm and Km is the decision point set of worker wm. For worker wm, our goal is to find an optimal policy πm∗ that maximizes his expected accumulation of rewards:(3)πm∗=argmaxπmE∑k∈KmRmskm,πmskm|s0m.


## 3.1. Notation

The main notations and their values used in this paper are summarized in Table1.Table 1
Notations.
NotationMeaningModelNvNumber of intersections of the road networkviIntersection of the road network, fori∈1,2,…,Nvv0The depotVSet of vertexes, including the intersections and depot,V=v0,v1,v2,…vNvvi,vjStreet segment (edge) of the road networkvi,vj∈V,ESet of edges,E=vi,vj|vi,vj∈Vandi<jvi,vjDirected street segment (arc),vi,vj∈VASet of arcs,A=vi,vj|vi,vj∈Vandi≠jGGraph of the road network,G=V,E in the undirected representation or G=V,A in the directed representationTNumber of time steps in an operating day (1440)tTime step on the small time scale,t∈1,2,…,TΔtInterval between time steps (1 min)tbeginningThe beginning of the working hours (7:00 am)tendThe end of the working hours (10:00 pm)ψmaxMaximum range of e-bike battery (40 km)ψminMinimum battery range threshold for an e-bike to be available (3 km)τbServing time cost of replacing battery for an e-bike (1 minute)NebNumber of e-bikes in system (1000)NwNumber of workers/vehicles in system (4)wiWorker,i∈1,2,…,NwWSet of workers,W=w1,w2,…wNwBBattery capacity of each vehicle (30)cijDeadheading time cost of arcvi,vjτdTime cost of reloading a battery in the depot (0.5 minutes)ψ¯userExpectation of battery range consumption for e-bike trips (2.75 km)btijNumber of batteries replaced on arcvi,vj at time step t, t∈1,2,…,TdtijNumber of successful pickups on arcvi,vj at time step t, t∈1,2,…,TGiSubzone of the road network,i∈1,2,…,NwViVertex set in subzoneGi, i∈1,2,…,NwAiArc set in subzoneGi, i∈1,2,…,NwEiEdge set in subzoneGi, i∈1,2,…,NwpTime period on the large time scaleΔpTime interval between time periods (180 minutes)KmThe decision point set of workerwm, m∈1,2,…,NwtkmTime step for thekth decision point of worker wm, tkm∈1,2,…,T, m∈1,2,…,Nw, k∈KmηkmA vector of lengthA indicating the available e-bike number on the arcs in the subzone of worker wm at the kth decision point (the values of arcs in other subzones are masked to 0), m∈1,2,…,Nw, k∈KmφkmA vector of lengthA indicating the underpowered e-bike number on the arcs in the subzone of worker wm at the kth decision point (the values of arcs in other subzones are masked to 0), m∈1,2,…,Nw, k∈KmLkmA one-hot vector of lengthA indicating which arc worker wm is located on at the kth decision point, m∈1,2,…,Nw, k∈KmCkmA nonnegative number indicating the remaining battery number in the vehicle of workerwm at the kth decision point, m∈1,2,…,Nw, k∈KmSmState set of workerwm, m∈1,2,…,NwskmState of workerwm at the kth decision point, skm∈Sm, and skm=tkm,ηkm,φkm,Lkm,Ckm, m∈1,2,…,Nw, k∈KmXskmFeasible decision set under stateskmxkmDecision of workerwm at the kth decision point, xkm is a one-hot vector of length A indicating the arc to visit next, xkm∈Xskm, m∈1,2,…,Nw, k∈KmrkmInstant reward for workerwm at the kth decision point, m∈1,2,…,Nw, k∈KmπmPolicy for workerwm, m∈1,2,…,Nwπm∗Optimal policy for workerwm, m∈1,2,…,NwMethodmijThe link relationship of theith and the jth edge in E, i,j∈1,2,…,EMAdjacency matrix of edge set,M=mijE×EeijA 0-1 variable to indicate whether thejth edge in E belongs to the subset Ei, i∈1,2,…,Nw, j∈1,2,…,EeiA 0-1 vector to indicate the edges in subsetEi, ei=eij∀j∈1,2,…,E, i∈1,2,…,NwupjThe number of underpowered e-bikes on thejth edge at time period p, j∈1,2,…,EupThe vector of the number of underpowered e-bikes on each edge at time periodp, up=upj∀j∈1,2,…,EliThe length of theith edge (km), i∈1,2,…,ElThe vector of the length of each edge,l=li∀i∈1,2,…,EβTemperature in the simulated annealing algorithmβ0Initial temperature in the simulated annealing algorithm (100)αThe decay factor in temperature in the simulated annealing algorithm (0.9)MaxIterinnerInner iteration limit in the simulated annealing algorithm (10)MaxIterouterOuter iteration limit in the simulated annealing algorithm (20)IrSimulation runs for each feasible decision in the rollout policyImaxMaximum iteration number in the MCTS algorithmHLookahead simulation horizons (180 time steps)γA weighting factor in the UCB algorithm (5)λA time discount parameter (0.95)


## 3.2. Problem Formulation

Figure3 shows a simple scenario of EBSS, which consists of the road network, e-bikes scattered among the virtual parking spots or used on the road, and vehicles/workers replacing batteries for underpowered e-bikes.Figure 3
A sample scenario of EBSS.The road network is abstracted as a connected graphG=V,E, where V=v0,v1,v2,…vNv is the vertex set and E=vi,vj|vi,vj∈Vandi<j is the edge set. In the vertex set, v0 is the depot and vi|i=1,2,…,Nv is the set of intersections. Assume that each edge has two directions and the vehicle can only serve one direction at a time. Thus, an edge vi,vj can be split into two directed arcs: vi,vj and vj,vi, as shown by the edge v1,v2 in Figure 3. Then, the road network can also be represented as a directed graph: G=V,A, where A=vi,vj|vi,vj∈Vandi≠j is the arc set. The undirected graph representation is used in the road network partitioning and the directed graph representation is used in the routing decision process.The operating day is divided intoT discrete time steps: t∈1,2,…,T, according to a predefined time interval Δt (1 minute). The dynamic user demand is reflected in the stochastically generated pickup and return demands in each virtual parking spot at each time step, which leads to dynamic changes in the location and battery power level of e-bikes.The number of e-bikes in the system is supposed to be fixed atNeb, without new e-bikes being put into the system or removed from the system. The battery range of each e-bike is in the range of 0,ψmax (km). When the battery range is below a threshold ψmin (km), an e-bike is regarded as underpowered and unavailable to users. For the sake of simplicity, we assume that the time cost for a worker to replace a battery is a constant: τb (e.g., 1 minute).Assume there areNw workers (or vehicles) in the system, denoted by the set W=w1,w2,…,wNw. Each vehicle has the same battery capacity of B (e.g., 40). Suppose tbegin (e.g., 7:00 am) and tend (e.g., 10:00 pm) are the beginning and the end of working hours. During working hours, workers will continue to replace batteries for underpowered e-bikes. Once the available batteries are used up, workers need to return to the depot and reload batteries. Suppose vehicles travel with a certain speed that is independent of time, leading to a positive deadheading cost cij of each arc vi,vj∈A. Assume that the time cost of reloading a battery in the depot is τd (e.g., 0.5 minutes) and there are always enough fully charged batteries in the depot, so the time cost of reloading a vehicle’s battery in the depot is τd⋅B.Under the limited working hours and time cost constraints introduced above, on the one hand, we hope to route the vehicles in a way that replaces as many batteries as possible during a limited time; on the other hand, we also hope to prioritize services for those arcs that have the potential to reduce user pickup failures in the following time. Take workerw1 in Figure 3 for example; if there will be a user picking up e-bike on arc v4,v7 in the following time, it would be better for w1 to serve arc v4,v7 first instead of arc v6,v7. However, directly maximizing the sum of the number of replaced batteries and the number of successful pickups may be not appropriate. In fact, replacing one more battery has a higher expected reward than reducing one more pickup failure because a fully charged e-bike can usually be used multiple times and can reduce more potential pickup failures in the following days. For this consideration, the optimization objective function is designed as a weighted sum of these two factors, as shown in the following equation:(1)max∑t∈1,2,…,T∑vi,vj∈Aψmax−ψmin⋅btij+ψ¯user⋅dtij,where btij is the number of batteries replaced on arc vi,vj at time step t, dtij is the number of successful pickups on arc vi,vj at time step t, and ψ¯user indicates the expectation of battery power consumption for each e-bike trip. Therefore, ∑t∈1,2,…,T∑vi,vj∈Aψ¯user⋅dtij denotes the battery power consumed by the users, indicating the short-term rewards during current day, while ∑t∈1,2,…,T∑vi,vj∈Aψmax−ψmin⋅btij denotes the battery power replenished by the workers, indicating the long-term rewards of the system sustainability.


## 3.3. Decomposition of the Problem

For the above dynamic capacitated arc routing problem, it is hard to compute the optimal route for each worker in advance due to the stochastic and complex nature of the system. Therefore, we seek to formalize the problem as a sequential decision process in cooperative multiagent system. A decision point occurs when a worker completes the task on the current arc and a decision about which one of the adjacent arcs to visit next should be made. It should be noted that as the number of tasks and deadheading cost varies for each arc, the executing duration of each decision may be different and the decision processes between different workers are thus asynchronous. At each time step, only a subset of the workers may make a decision while others keep on executing the tasks on the current arcs.To simplify the problem, we assume that each worker is assigned a nonoverlapping subzone and is allowed to perform the battery replacement operation only within his subzone. As there areNw workers, the road network G=V,A will be partitioned into Nw subzones: G1,G2,…,GNw, satisfying the condition that Gi=Vi,Ai is strongly connected graph for i=1,2,…,Nw and A1∪A2∪…∪ANw=A, Ai∪Aj=0 for i≠j. In this case, the routing decisions of each worker have very limited effects on other workers, and we can decompose the cooperative multiagent decision problem in the whole road network into multiple mono-agent decision problems in each subzone approximatively.As shown in Figure4, the arc set A is partitioned into two subsets A1 and A2 for workers w1 and w2, and the decision process of each worker is modeled as a Markov Decision Process. The operating day is divided into discrete time steps t according to a small time interval Δt (e.g., 1 minute). At each time step, a worker can make a routing decision or continue with his current task. For example, the decision points for worker w1 in Figure 4 are time steps t+2, t+4, and t+6. We also divide the day into some discrete time periods p according to a larger time interval Δp (e.g., 180 minutes). At each time period (e.g., at p1 and p2 in Figure 4), the partition can be adjusted according to the real-time e-bike state.Figure 4
An illustration of the decomposition of decision processes.To model the routing decision process of each worker as a MDP, the following assumptions are made:(i)
After a worker leaves the depot, he follows the shortest path (i.e., calculated by the Dijkstra algorithm) to the nearest vertex of his subzone. During this process, the worker does not replace battery. Once the remaining batteries in the vehicle run out, the worker follows the shortest path back to the depot to reload the batteries and then goes to the previous location and continues his operation.(ii)
As new underpowered e-bikes may be generated on the arcs that have been served, multiple visits and multiple services to an arc are allowed. However, we suppose that temporary stops of vehicles are not allowed and that the vehicle’s power is always sufficient during working hours.Based on the above assumptions, the MDP for themth worker wm consists of the following components:A state:Sm is the feasible state set of worker wm. skm∈Sm is state of the kth decision point of worker wm. The state skm=tkm,ηkm,φkm,Lkm,Ckmincludes three main parts: the system time, the state of e-bikes, and the state of workers. The time step of the kth decision point of worker wm is denoted as tkm; the state of e-bikes includes ηkm and φkm, indicating the available e-bike number and the underpowered e-bike number on the arcs in the subzone of worker wm at the kth decision point, respectively. The state of worker includes the locations Lkm and the remaining battery number Ckm of worker wm at the kth decision point. A detailed description of the state information is given in Table 1.A decision:Xskm denotes the feasible decision set under state skm. Xskm is the intersection of the set of arcs adjacent to the current arc of wm and the set of arcs in the subzone of wm. For example, in Figure 4, v4,v0,v4,v1,v4,v7 is the set of arcs adjacent to v7,v4. However, v4,v1 belongs to different subzones from v7,v4. Therefore, when a worker is making a routing decision on arc v7,v4, the feasible decision set is v4,v0,v4,v7. The kth decision of worker wm is denoted as xkm∈Xskm, representing the next arc to visit.State transition function:Pmsk+1m|skm,xkm is the transition probability function for worker wm and describes the probability that the system will move to the next decision state sk+1m after making a decision xkm under the current decision state skm. The next decision state sk+1m will also be influenced by the stochastic exogenous information in the system (e.g., the stochastic user demand) and the operations of other workers.Reward function:rkm=Rmskm,xkm is the instant reward for worker wm at the kth decision point. Rmskm,xkm denotes the reward function for worker wm, reflecting the instant reward received after taking an decision xkm under state skm. rkm can be derived from the weighted sum of the number of successful pickups and the number of batteries replaced within the subzone of wm between time step tkm and tk+1m (if k is the last decision point of wm, then tk+1m is the last time step of the day), as shown in the following equation:(2)rkm=∑t=tkmtk+1m−1∑vi,vj∈Amψmax−ψmin⋅btij+ψ¯user⋅dtij.Since the subzones do not overlap, the cumulative reward in a subzone is largely related to the decisions of the worker in that subzone and less related to the decisions of other workers, and the sum of the cumulative rewards of all subzones is the total system reward, which is consistent with the optimization objective in equation (1).Policy:πm describes the policy to make a decision under current state for worker wm: xkm=πmskm. Suppose that s0m is the initial state of worker wm and Km is the decision point set of worker wm. For worker wm, our goal is to find an optimal policy πm∗ that maximizes his expected accumulation of rewards:(3)πm∗=argmaxπmE∑k∈KmRmskm,πmskm|s0m.


## 4. Spatiotemporal Partitioning Method

As introduced in the previous sections, on the upper level of the solution framework, we decompose the decision processes of multiple workers by partitioning the road network into subzones. Therefore, the way the road network is partitioned has a great impact on the allocation of labor force and operation efficiency. In this section, we first propose the objective function of the road network partitioning in Section4.1 and then introduce the algorithm of initial partition generation and partition adjustment at each time period in Sections 4.2 and 4.3.


### 4.1. Objective

The objective of road network partitioning should not only take into account the static topology of the road network but also make the amount of tasks, i.e., the number of underpowered e-bikes, in each subzone as consistent as possible. This is to avoid situations where some workers are overworked while others are underworked. We use an edge (i.e., a street segment) as the partition unit, so the two arcs of an edge will always be divided into the same subzone. Therefore, the edge setE will be partitioned into Nw (i.e., the number of workers) subsets: E1,E2,…,ENw, satisfying that E1∪E2∪…∪ENw=E, Ei∩Ej=0 for i≠j, and the subgraph induced by Ei is connected graph for i=1,2,…,Nw. If the subgraph Gi induced by above Ei is represented in the form of a directed graph: Gi=Vi,Ai, it is obvious that Gi=Vi,Ai satisfies the constraints for the subzone described in Section 3.3.Since the graph is partitioned on edges instead of vertexes, here we use the adjacent matrixM=mijE×E to represent the link relationship of edges in set E. mij=1 indicates that the ith and jth edges in E are linked; otherwise, mij=0. An edge is linked to another edge when and only when they have common vertexes. Suppose an edge is not linked to itself. For example, in Figure 3, edge v5,v8 is linked to edges v0,v5, v2,v5, v6,v8, and v8,v9.We chose three indicators to evaluate the road network partition: the first indicator is that the number of underpowered e-bikes in each subzone should be as similar as possible, which facilitates the efficiency of battery replacement and makes the workers’ income (which is usually related to how many batteries they replaced) as even as possible. The second indicator is that the sum of the lengths of the edges in each subzone should be as similar as possible to make each subzone similar in size. The third indicator is that the links between subzones should be as few as possible, which ensures that the street segments in each subzone are as closely connected as possible.The 0-1 vectorsei=eij∀j∈1,2,…,E for i=1,2,…,Nw are the decision variables for the partitioning, indicating the edges belonging to subset Ei. eij=1 indicates the jth edge belonging to the subset Ei; otherwise, eij=0. Also, we use vector up=upi∀i∈1,2,…,E to represent the real-time underpowered e-bike number on each edge at time period p and use l=li∀i∈1,2,…,E to represent the length of each edge. Let g represent a partition solution, i.e., a realization of the decision variables e1,e2,…,eNw. At a new time period p, the following objective of the partition is designed to be optimized:(4)mine1,e2,…,eNwTaskCutg=ϕ1Nw∑i=1NweiTup1Tup−1Nw2+ϕ2Nw∑i=1NweiTl1Tl−1Nw2+ϕ3∑i=1NweiTM1−ei1TM1,(5)s.t.∑i=1Nweij=1,∀j∈1,2,…,E,(6)eij∈0,1,∀i∈1,2,…,Nw,∀j∈1,2,…,E,(7)ConnectedMei,ei=1,∀i∈1,2,…,Nw,where eiTup denotes the underpowered e-bike number in Ei at time period p, eiTl denotes the total length of edges of Ei, eiTM1−ei denotes the link number between Ei and the complement of Ei in E, and ϕ1, ϕ2, and ϕ3 are the weight coefficients of the three indicators, respectively. Therefore, in equation (4), the first term represents the variance of the underpowered e-bike number (after normalization) in each subzone, the second term is the variance of the total length of edges (after normalization) in each subzone, the third term is the ratio of the link number between subzones and the total link number in the road network.Equation (5) ensures that each edge can only be divided into exactly one subzone, and equation (6) ensures the decision variables are 0-1 integers.In equation (7), Mei,ei is a submatrix of M which consists of the jth row and column of M where eij=1, representing the adjacent matrix of edges in set Ei. The ConnectedMei,ei function is used to determine whether the subgraph induced by the adjacent matrix Mei,ei is connected (e.g., according to the depth-first search algorithm). ConnectedMei,ei=1 indicates the subgraph induced by matrix Mei,ei is connected; otherwise, ConnectedMei,ei=0.Since the amount of tasks in each subzone depends on the real-time state of e-bikes and may change over time, it may not be appropriate to partition the road network in a fixed way throughout the day. However, the frequent changes to the road network partition may cause confusion to workers. Therefore, in the proposed spatiotemporal partitioning method, an initial partition solution is first generated at the beginning of the operating day and then will be adjusted dynamically in the large time scale.


### 4.2. Initial Partition Generation

When generating the initial road network partition, we first ignore the dynamical e-bike state and only consider the static topology of the road network to ensure that the initial partition is the same at the beginning of each day. We adopt the solution of the optimization problem (8) with constraints (5) and (6) as the initial road network partition.(8)mine1,e2,…,eNwRatioCutg=12∑i=1NweiTM1−ei1Tei.TheRatioCut problem in equation (8) is a classic graph-cutting problem and attempts to minimize the links between subzones while making the number of edges in each subzone as close as possible. The number of edges in a subzone is often correlated with the number of tasks and the total length of the edges in this subzone. Various algorithms to this problem have been proposed in the past and we adopt the solution based on spectral clustering algorithm, which can be referred in the literature [42, 43]. If there is a subgraph in the initial partition that is not connected, the nonmaximal connected components of this subgraph will be assigned to other adjacent subgraphs to ensure that all subgraphs in the initial partition are connected.


### 4.3. Partition Adjustment Algorithm

For a new time periodp and an existing road network partition g of the previous time period (for period p0, the previous road network partition is the initial partition), we first determine whether the existing partition needs to be adjusted according to the real-time e-bike state, and if so, we adjust the existing partition through minimizing the objective in equation (4).In the last two decades, various algorithms for community detection or graph partition have been proposed, including a large part of optimization-based methods, i.e., the modularity optimization method [44–47]. With reference to literature [46], a simulated annealing heuristic algorithm is adopted to adjust the road network partition. The simulated annealing algorithm is an effective heuristic algorithm that can effectively avoid falling into local minima and eventually converge to the global optimum. The process is shown in Algorithm 1.Algorithm 1: Simulated annealing-based partition adjustment algorithm.
Input:up at time period p, the edge length l, existing partition g, the initial temperature β0, MaxIterinner, MaxIterouter, α∈0,1, TaskCutminOutput: new road network partition gbestIfTaskCutg≤TaskCutmin:Keep the existing partitiong as gbestElse:Set temperatureβ=β0Repeat:Repeat:Create a new partitiong' from g using task balance operatorCalculateΔtask=TaskCutg′−TaskCutgAcceptg=g′ with probability of min1,exp−Δtask/βCreate a new partitiong' from g using length balance operatorCalculateΔtask=TaskCutg′−TaskCutg and accept g=g′ with probability of min1,exp−Δtask/βUntil iteration limit MaxIterinner reachedCreate a new partitiong' from g using local search operatorAcceptg=g′ if TaskCutg′<TaskCutgSetβ=αΔβUsegbest to record the best partition solution foundUntil iteration limit MaxIterouter reachedOutput the new road network partition:gbestIn Algorithm1, we use three operators to generate new solutions, i.e., the task balance operator, the length balance operator, and the local search operator. The task balance operator consists of the following steps as illustrated in Figure 5:Step 1: for the existing partitiong, the TaskCutg is calculated according to equation (4), as shown in Figure 5(a).Step 2: select two adjacent subzones randomly according to their difference in underpowered e-bike number. Here we ensure that the adjacent subzones with larger difference in underpowered e-bike number have a higher probability of being selected (e.g., through Boltzmann exploration strategy). Suppose the subzonesE2 and E3 in Figure 5(b) were selected. As e2Tup>e3Tup, we select the edges in E2 that are linked to the edges in E3 as the removable edge set of E2, i.e., the edges v0,v3 and v2,v5.Step 3: we randomly select an edge from the removable edge set ofE2 and transfer it to E3, such as the edge v2,v5. It should be noted that the edge selected here cannot change the connectivity of E2 after being transferred, i.e., satisfying the constraint in equation (7). Then, a new partition g′ is created as shown in Figure 5(c).Figure 5
An illustration of the task balance operator. The first number on each edge denotes the number of underpowered e-bikes of the edge at time periodp, and the second number on each edge denotes the length of the edge (km). Here, ϕ1,ϕ2, and ϕ3 are set to 100, 100, and 1, respectively. (a) Existing partition. (b) Select subzones. (c) Transfer edge.
(a)(b)(c)Thelength balance operator and the local search operator are the same as the task balance operator except for step 2 above. In the length balance operator, we select two adjacent subzones for edge transfer randomly according to their difference in the sum of edge lengths. In the local search operator, we select two adjacent subzones for the edge transfer completely at random.We define theTaskCutmin in Algorithm 1 as the minimum TaskCut value at the previous time period. This suggests that no change to the existing partition is needed as long as the variance of the underpowered e-bike number in each subzone at the current period is not greater than the previous time period. For period p0, we set TaskCutmin=0.


## 4.1. Objective

The objective of road network partitioning should not only take into account the static topology of the road network but also make the amount of tasks, i.e., the number of underpowered e-bikes, in each subzone as consistent as possible. This is to avoid situations where some workers are overworked while others are underworked. We use an edge (i.e., a street segment) as the partition unit, so the two arcs of an edge will always be divided into the same subzone. Therefore, the edge setE will be partitioned into Nw (i.e., the number of workers) subsets: E1,E2,…,ENw, satisfying that E1∪E2∪…∪ENw=E, Ei∩Ej=0 for i≠j, and the subgraph induced by Ei is connected graph for i=1,2,…,Nw. If the subgraph Gi induced by above Ei is represented in the form of a directed graph: Gi=Vi,Ai, it is obvious that Gi=Vi,Ai satisfies the constraints for the subzone described in Section 3.3.Since the graph is partitioned on edges instead of vertexes, here we use the adjacent matrixM=mijE×E to represent the link relationship of edges in set E. mij=1 indicates that the ith and jth edges in E are linked; otherwise, mij=0. An edge is linked to another edge when and only when they have common vertexes. Suppose an edge is not linked to itself. For example, in Figure 3, edge v5,v8 is linked to edges v0,v5, v2,v5, v6,v8, and v8,v9.We chose three indicators to evaluate the road network partition: the first indicator is that the number of underpowered e-bikes in each subzone should be as similar as possible, which facilitates the efficiency of battery replacement and makes the workers’ income (which is usually related to how many batteries they replaced) as even as possible. The second indicator is that the sum of the lengths of the edges in each subzone should be as similar as possible to make each subzone similar in size. The third indicator is that the links between subzones should be as few as possible, which ensures that the street segments in each subzone are as closely connected as possible.The 0-1 vectorsei=eij∀j∈1,2,…,E for i=1,2,…,Nw are the decision variables for the partitioning, indicating the edges belonging to subset Ei. eij=1 indicates the jth edge belonging to the subset Ei; otherwise, eij=0. Also, we use vector up=upi∀i∈1,2,…,E to represent the real-time underpowered e-bike number on each edge at time period p and use l=li∀i∈1,2,…,E to represent the length of each edge. Let g represent a partition solution, i.e., a realization of the decision variables e1,e2,…,eNw. At a new time period p, the following objective of the partition is designed to be optimized:(4)mine1,e2,…,eNwTaskCutg=ϕ1Nw∑i=1NweiTup1Tup−1Nw2+ϕ2Nw∑i=1NweiTl1Tl−1Nw2+ϕ3∑i=1NweiTM1−ei1TM1,(5)s.t.∑i=1Nweij=1,∀j∈1,2,…,E,(6)eij∈0,1,∀i∈1,2,…,Nw,∀j∈1,2,…,E,(7)ConnectedMei,ei=1,∀i∈1,2,…,Nw,where eiTup denotes the underpowered e-bike number in Ei at time period p, eiTl denotes the total length of edges of Ei, eiTM1−ei denotes the link number between Ei and the complement of Ei in E, and ϕ1, ϕ2, and ϕ3 are the weight coefficients of the three indicators, respectively. Therefore, in equation (4), the first term represents the variance of the underpowered e-bike number (after normalization) in each subzone, the second term is the variance of the total length of edges (after normalization) in each subzone, the third term is the ratio of the link number between subzones and the total link number in the road network.Equation (5) ensures that each edge can only be divided into exactly one subzone, and equation (6) ensures the decision variables are 0-1 integers.In equation (7), Mei,ei is a submatrix of M which consists of the jth row and column of M where eij=1, representing the adjacent matrix of edges in set Ei. The ConnectedMei,ei function is used to determine whether the subgraph induced by the adjacent matrix Mei,ei is connected (e.g., according to the depth-first search algorithm). ConnectedMei,ei=1 indicates the subgraph induced by matrix Mei,ei is connected; otherwise, ConnectedMei,ei=0.Since the amount of tasks in each subzone depends on the real-time state of e-bikes and may change over time, it may not be appropriate to partition the road network in a fixed way throughout the day. However, the frequent changes to the road network partition may cause confusion to workers. Therefore, in the proposed spatiotemporal partitioning method, an initial partition solution is first generated at the beginning of the operating day and then will be adjusted dynamically in the large time scale.


## 4.2. Initial Partition Generation

When generating the initial road network partition, we first ignore the dynamical e-bike state and only consider the static topology of the road network to ensure that the initial partition is the same at the beginning of each day. We adopt the solution of the optimization problem (8) with constraints (5) and (6) as the initial road network partition.(8)mine1,e2,…,eNwRatioCutg=12∑i=1NweiTM1−ei1Tei.TheRatioCut problem in equation (8) is a classic graph-cutting problem and attempts to minimize the links between subzones while making the number of edges in each subzone as close as possible. The number of edges in a subzone is often correlated with the number of tasks and the total length of the edges in this subzone. Various algorithms to this problem have been proposed in the past and we adopt the solution based on spectral clustering algorithm, which can be referred in the literature [42, 43]. If there is a subgraph in the initial partition that is not connected, the nonmaximal connected components of this subgraph will be assigned to other adjacent subgraphs to ensure that all subgraphs in the initial partition are connected.


## 4.3. Partition Adjustment Algorithm

For a new time periodp and an existing road network partition g of the previous time period (for period p0, the previous road network partition is the initial partition), we first determine whether the existing partition needs to be adjusted according to the real-time e-bike state, and if so, we adjust the existing partition through minimizing the objective in equation (4).In the last two decades, various algorithms for community detection or graph partition have been proposed, including a large part of optimization-based methods, i.e., the modularity optimization method [44–47]. With reference to literature [46], a simulated annealing heuristic algorithm is adopted to adjust the road network partition. The simulated annealing algorithm is an effective heuristic algorithm that can effectively avoid falling into local minima and eventually converge to the global optimum. The process is shown in Algorithm 1.Algorithm 1: Simulated annealing-based partition adjustment algorithm.
Input:up at time period p, the edge length l, existing partition g, the initial temperature β0, MaxIterinner, MaxIterouter, α∈0,1, TaskCutminOutput: new road network partition gbestIfTaskCutg≤TaskCutmin:Keep the existing partitiong as gbestElse:Set temperatureβ=β0Repeat:Repeat:Create a new partitiong' from g using task balance operatorCalculateΔtask=TaskCutg′−TaskCutgAcceptg=g′ with probability of min1,exp−Δtask/βCreate a new partitiong' from g using length balance operatorCalculateΔtask=TaskCutg′−TaskCutg and accept g=g′ with probability of min1,exp−Δtask/βUntil iteration limit MaxIterinner reachedCreate a new partitiong' from g using local search operatorAcceptg=g′ if TaskCutg′<TaskCutgSetβ=αΔβUsegbest to record the best partition solution foundUntil iteration limit MaxIterouter reachedOutput the new road network partition:gbestIn Algorithm1, we use three operators to generate new solutions, i.e., the task balance operator, the length balance operator, and the local search operator. The task balance operator consists of the following steps as illustrated in Figure 5:Step 1: for the existing partitiong, the TaskCutg is calculated according to equation (4), as shown in Figure 5(a).Step 2: select two adjacent subzones randomly according to their difference in underpowered e-bike number. Here we ensure that the adjacent subzones with larger difference in underpowered e-bike number have a higher probability of being selected (e.g., through Boltzmann exploration strategy). Suppose the subzonesE2 and E3 in Figure 5(b) were selected. As e2Tup>e3Tup, we select the edges in E2 that are linked to the edges in E3 as the removable edge set of E2, i.e., the edges v0,v3 and v2,v5.Step 3: we randomly select an edge from the removable edge set ofE2 and transfer it to E3, such as the edge v2,v5. It should be noted that the edge selected here cannot change the connectivity of E2 after being transferred, i.e., satisfying the constraint in equation (7). Then, a new partition g′ is created as shown in Figure 5(c).Figure 5
An illustration of the task balance operator. The first number on each edge denotes the number of underpowered e-bikes of the edge at time periodp, and the second number on each edge denotes the length of the edge (km). Here, ϕ1,ϕ2, and ϕ3 are set to 100, 100, and 1, respectively. (a) Existing partition. (b) Select subzones. (c) Transfer edge.
(a)(b)(c)Thelength balance operator and the local search operator are the same as the task balance operator except for step 2 above. In the length balance operator, we select two adjacent subzones for edge transfer randomly according to their difference in the sum of edge lengths. In the local search operator, we select two adjacent subzones for the edge transfer completely at random.We define theTaskCutmin in Algorithm 1 as the minimum TaskCut value at the previous time period. This suggests that no change to the existing partition is needed as long as the variance of the underpowered e-bike number in each subzone at the current period is not greater than the previous time period. For period p0, we set TaskCutmin=0.


## 5. Lookahead Policy with MCTS Algorithm

After partitioning the road network into subzones, on the lower level of the solution framework, we model the routing decision process of each worker in their subzone as a mono-agent MDP and adopt a lookahead policy with MCTS algorithm; the lookahead policy based on the rollout heuristic is introduced in Section5.1. A MCTS algorithm that can improve the simulation efficiency is proposed in Section 5.2. A heuristic pruning rule is added to the MCTS algorithm in Section 5.3.


### 5.1. Lookahead Policy

This paper takes an arc as the unit of the routing decision process as the two arcs of an edge generally cannot be served at the same time. We useQk∗skm,xkm to denote the optimal Q-value of the state-decision pair, representing the expected rewards obtained along the optimal policy after executing decision xkm under the state skm. Thus, the optimal decision can be deduced by the following equation:(9)xkm=argmaxx∈XskmQk∗skm,x.However, it is difficult to obtain the optimalQ-value in our problem. The difficulty lies not only in the complexity of the problem but also in the stochastic and time-dependent user demands.This paper adopts a lookahead policy for our problem. The rollout policy is one of the lookahead policies [18, 48] and will estimate the approximate value of a feasible decision from the simulated experience generated in a simulator. The simulator will generate the stochastic pickup and return demands, according to a priori knowledge or the prediction of the realistic environment. The rollout policy is suitable for problems where the state space is very large.The process of the lookahead policy with rollout heuristic is shown in Figure6. For the kth decision point of worker wm, we first simulate from each feasible decision x for a predefined simulation run Ir and then use the Q-value calculated by equation (10) to approximate the optimal Q-value in equation (9).(10)Qkskm,x=1Ir∑j=1Ir∑i=kkmaxjλi−k⋅rijm∀x∈Xskm,where λ is a time discount factor to reduce the impact of rewards in future because the uncertainty gets higher further into the future. As there are a large number of time steps in the operating day, we truncate the rollout trajectories and set the simulation horizon to H time steps. H is set to 180 according to preliminary tests. kmaxj is the terminal decision point within the horizon H in the jth simulation trajectory. rijm is the reward for the ith decision point of worker wm in the jth simulation trajectory. For the decision point i=k+1,k+2,…,kmaxj, the decision is made following the base heuristic polikers also follow the base heuristic policy when making decisions during the simulation.Algorithm 2: Base heuristic policy
Input: state skm=tkm,ηkm,φkm,Lkm,Ckm, current road network partition gOutput: decision xkm and calculate the feasible decision set Xskm based on the current worker’s location Lkm and the road network partition g;(1)
Calculate the ratio of the underpowered e-bike number and available e-bike numberudpx on each feasible arc x∈Xskm based on the e-bike state ηkm,φkm.(2)
If the available e-bike number is 0, we replace it with a sufficiently small value.(3)
Select a decisionxkm with Boltzmann exploration Pxkm|skm=expudpxkm/∑x∈Xskmexpudpx.Figure 6
A lookahead policy with rollout heuristic.


### 5.2. MCTS Algorithm

The rollout policy spends the same simulation runs for each feasible decision, resulting in a lack of exploitation of the decisions that appear to be promising and much simulation for some decisions that may be clearly inferior to others [48]. Therefore, this paper adopts the MCTS sampling algorithm, a variant of the rollout policy, that can better the trade-off between exploration and exploitation when allocating simulation resources.When we are making a routing decision, we are also attempting to find an arc belonging to the optimal route. Constructing a decision tree of all feasible decisions is equivalent to traversing all possible routes that a worker can take in future.In the decision tree, a node represents an arc. The root node (denoted asnode1) represents the current arc where the worker is located at. The child nodes of a given node indicate the arcs that can be expanded from the given arc. All children of the root node form the feasible decision set. The nodes between the root node and a leaf node form a possible route: ℒ=node1,node2,…. It is worth noting that in the simulation horizon, if the worker runs out of batteries, he will return to the depot immediately to reload the batteries and then go to his subzone to continue the operation. The arcs passed during this process will not be recorded in the decision tree.On the basis of the rollout policy, the MCTS algorithm will maintain a decision tree during the simulation. Each node in the decision tree has two attributes:ni denotes the number of times that nodei is visited and scorei denotes the value of nodei. The root node and its child nodes are initialized with ni=0,scorei=0.MCTS is an iterative algorithm, where the maximum iteration number is supposed to be set toImax. Each iteration contains the following four main parts: tree traversal, node expansion, random simulation, and backpropagation:(i)
Tree traversal: tree traversal starts from the root node and selects the child nodes layer by layer until the current node is a leaf node. One criterion for selecting child node is the UCB algorithm [49]:(11)UCBnodei=scorei+γ2lnNni,ifni>0,∞,ifni=0,wherenodeiis a child node, N denotes the total visited times of the parent node, and γ is a weighting factor. The child node with the largest UCB value will be selected. The simulator will simulate the user demand and the operation of workers after a child node is selected. When the worker completes the task, a new child node will be selected until the current node is a leaf node.(ii)
Node expansion: if the current leaf node has not been visited, then go to the next step; if the current leaf node has been visited, then expand its children to the decision tree and set the visited time and the value of these child nodes to 0. A child node is selected at random for the next step.(iii)
Random simulation: random simulation is performed from the current leaf node until the simulation horizonH according to the base heuristic policy in Algorithm 2. Suppose the worker is wm and the current decision point is k. Until the end of the simulation, a route ℒ of wm will be obtained. Suppose the terminal decision point in the simulation is kmax. The Scoreℒ in equation (12) is used to evaluate route ℒ:(12)Scoreℒ=∑i=kkmaxλi−k⋅rim.(iv)
Backpropagation: after getting theScoreℒ, the nodes of the decision tree that are also in route ℒ will be updated according to equations (13) and (14):(13)ni=ni+1∀nodeiinL,(14)scorei=ni−1⋅scorei+Scoreℒni∀nodeiinL.After performing the backpropagation, the decision tree is updated. Then, it returns to the tree traversal in step (1) and continues iteration until the iteration number reachesImax. Finally, the child node of the root node with the highest value will be selected as the next arc to visit.


### 5.3.k-Cycle Pruning Rule

The rules for expanding the child nodes were introduced earlier, i.e., expanding the arcs that are adjacent to the parent node and in the same subzone. However, following this relatively loose rule may generate some apparently low-value routes, such as circular routes. Take Figure3 for example; the circular route v7,v4,v4,v7,v7,v4,v4,v7,… may yield low reward because of repeated visits to the same arcs during a short time. To avoid such situation, we propose a k-cycle pruning rule.Fornodei at the ith level of the decision tree (suppose the root node node1 is the first level), we denote the route from the root node to it by ℒ=node1,node2,…,nodei−1,nodei. If a child node of nodei is in the set of nodei+1−mink,i,nodei+2−mink,i,…,nodei−1,nodei, this child node will not be expanded (unless there is no other expandable child node). k is an empirical constant and is tuned artificially. The routes with circles of length less than or equal to k will be pruned from the decision tree. For example, if k is set to 2, the decision tree from the arc v7,v4 in Figure 3 will be cut as shown in Figure 7. In the experiment, k is tuned to 6 according to preliminary tests.Figure 7
An example of thek-cycle pruning rule.


## 5.1. Lookahead Policy

This paper takes an arc as the unit of the routing decision process as the two arcs of an edge generally cannot be served at the same time. We useQk∗skm,xkm to denote the optimal Q-value of the state-decision pair, representing the expected rewards obtained along the optimal policy after executing decision xkm under the state skm. Thus, the optimal decision can be deduced by the following equation:(9)xkm=argmaxx∈XskmQk∗skm,x.However, it is difficult to obtain the optimalQ-value in our problem. The difficulty lies not only in the complexity of the problem but also in the stochastic and time-dependent user demands.This paper adopts a lookahead policy for our problem. The rollout policy is one of the lookahead policies [18, 48] and will estimate the approximate value of a feasible decision from the simulated experience generated in a simulator. The simulator will generate the stochastic pickup and return demands, according to a priori knowledge or the prediction of the realistic environment. The rollout policy is suitable for problems where the state space is very large.The process of the lookahead policy with rollout heuristic is shown in Figure6. For the kth decision point of worker wm, we first simulate from each feasible decision x for a predefined simulation run Ir and then use the Q-value calculated by equation (10) to approximate the optimal Q-value in equation (9).(10)Qkskm,x=1Ir∑j=1Ir∑i=kkmaxjλi−k⋅rijm∀x∈Xskm,where λ is a time discount factor to reduce the impact of rewards in future because the uncertainty gets higher further into the future. As there are a large number of time steps in the operating day, we truncate the rollout trajectories and set the simulation horizon to H time steps. H is set to 180 according to preliminary tests. kmaxj is the terminal decision point within the horizon H in the jth simulation trajectory. rijm is the reward for the ith decision point of worker wm in the jth simulation trajectory. For the decision point i=k+1,k+2,…,kmaxj, the decision is made following the base heuristic polikers also follow the base heuristic policy when making decisions during the simulation.Algorithm 2: Base heuristic policy
Input: state skm=tkm,ηkm,φkm,Lkm,Ckm, current road network partition gOutput: decision xkm and calculate the feasible decision set Xskm based on the current worker’s location Lkm and the road network partition g;(1)
Calculate the ratio of the underpowered e-bike number and available e-bike numberudpx on each feasible arc x∈Xskm based on the e-bike state ηkm,φkm.(2)
If the available e-bike number is 0, we replace it with a sufficiently small value.(3)
Select a decisionxkm with Boltzmann exploration Pxkm|skm=expudpxkm/∑x∈Xskmexpudpx.Figure 6
A lookahead policy with rollout heuristic.


## 5.2. MCTS Algorithm

The rollout policy spends the same simulation runs for each feasible decision, resulting in a lack of exploitation of the decisions that appear to be promising and much simulation for some decisions that may be clearly inferior to others [48]. Therefore, this paper adopts the MCTS sampling algorithm, a variant of the rollout policy, that can better the trade-off between exploration and exploitation when allocating simulation resources.When we are making a routing decision, we are also attempting to find an arc belonging to the optimal route. Constructing a decision tree of all feasible decisions is equivalent to traversing all possible routes that a worker can take in future.In the decision tree, a node represents an arc. The root node (denoted asnode1) represents the current arc where the worker is located at. The child nodes of a given node indicate the arcs that can be expanded from the given arc. All children of the root node form the feasible decision set. The nodes between the root node and a leaf node form a possible route: ℒ=node1,node2,…. It is worth noting that in the simulation horizon, if the worker runs out of batteries, he will return to the depot immediately to reload the batteries and then go to his subzone to continue the operation. The arcs passed during this process will not be recorded in the decision tree.On the basis of the rollout policy, the MCTS algorithm will maintain a decision tree during the simulation. Each node in the decision tree has two attributes:ni denotes the number of times that nodei is visited and scorei denotes the value of nodei. The root node and its child nodes are initialized with ni=0,scorei=0.MCTS is an iterative algorithm, where the maximum iteration number is supposed to be set toImax. Each iteration contains the following four main parts: tree traversal, node expansion, random simulation, and backpropagation:(i)
Tree traversal: tree traversal starts from the root node and selects the child nodes layer by layer until the current node is a leaf node. One criterion for selecting child node is the UCB algorithm [49]:(11)UCBnodei=scorei+γ2lnNni,ifni>0,∞,ifni=0,wherenodeiis a child node, N denotes the total visited times of the parent node, and γ is a weighting factor. The child node with the largest UCB value will be selected. The simulator will simulate the user demand and the operation of workers after a child node is selected. When the worker completes the task, a new child node will be selected until the current node is a leaf node.(ii)
Node expansion: if the current leaf node has not been visited, then go to the next step; if the current leaf node has been visited, then expand its children to the decision tree and set the visited time and the value of these child nodes to 0. A child node is selected at random for the next step.(iii)
Random simulation: random simulation is performed from the current leaf node until the simulation horizonH according to the base heuristic policy in Algorithm 2. Suppose the worker is wm and the current decision point is k. Until the end of the simulation, a route ℒ of wm will be obtained. Suppose the terminal decision point in the simulation is kmax. The Scoreℒ in equation (12) is used to evaluate route ℒ:(12)Scoreℒ=∑i=kkmaxλi−k⋅rim.(iv)
Backpropagation: after getting theScoreℒ, the nodes of the decision tree that are also in route ℒ will be updated according to equations (13) and (14):(13)ni=ni+1∀nodeiinL,(14)scorei=ni−1⋅scorei+Scoreℒni∀nodeiinL.After performing the backpropagation, the decision tree is updated. Then, it returns to the tree traversal in step (1) and continues iteration until the iteration number reachesImax. Finally, the child node of the root node with the highest value will be selected as the next arc to visit.


## 5.3.k-Cycle Pruning Rule

The rules for expanding the child nodes were introduced earlier, i.e., expanding the arcs that are adjacent to the parent node and in the same subzone. However, following this relatively loose rule may generate some apparently low-value routes, such as circular routes. Take Figure3 for example; the circular route v7,v4,v4,v7,v7,v4,v4,v7,… may yield low reward because of repeated visits to the same arcs during a short time. To avoid such situation, we propose a k-cycle pruning rule.Fornodei at the ith level of the decision tree (suppose the root node node1 is the first level), we denote the route from the root node to it by ℒ=node1,node2,…,nodei−1,nodei. If a child node of nodei is in the set of nodei+1−mink,i,nodei+2−mink,i,…,nodei−1,nodei, this child node will not be expanded (unless there is no other expandable child node). k is an empirical constant and is tuned artificially. The routes with circles of length less than or equal to k will be pruned from the decision tree. For example, if k is set to 2, the decision tree from the arc v7,v4 in Figure 3 will be cut as shown in Figure 7. In the experiment, k is tuned to 6 according to preliminary tests.Figure 7
An example of thek-cycle pruning rule.


## 6. System Simulator

Referring to literature [50], when there is no real e-bike trip data or no condition to test the proposed policy in reality, we turn our focus to simulation, which is an integral part of most studies of dynamic planning in complex systems. A system simulator is designed to capture the basic dynamics of the interactions of road network, e-bikes, workers, and users, which is shown in Figure 8.Figure 8
The simulation flowchart.In the “road network partition” module, an initial road network partition will be generated at the beginning of the operating day and will be adjusted every new period.An e-bike in the EBSS has three attributes: location, battery power level, and usage state. It is assumed that users only pick up or return e-bikes in the virtual parking spots. Therefore, e-bikes are only parked within the virtual parking spots or used on the road. The configuration at the beginning of the simulation includes the number of e-bikes in each virtual parking spot and the initial battery power level of each e-bike. In each time step, the attributes of the e-bike will be changed by two factors: the worker’s operation and the user’s usage.In the “worker’s operation” module, the simulator will update the worker’s operations one by one. Suppose a worker is operating on arc vi,vj; the deadheading cost cij is used as the expected travel time for the worker on arc vi,vj. At every time step, if there is an underpowered e-bike on arc vi,vj, the worker will spend τb minutes to replace the battery for this e-bike. Also, the battery power level of the e-bike will be updated to ψmax accordingly. If there is no underpowered e-bike, the expected travel time will be reduced by one time step. If the worker’s remaining travel time on the current arc becomes 0, a decision point occurs and the worker will select an adjacent arc to visit next according to the proposed policy.In the “pickup demand” module, the simulator will generate a batch of new orders every time step. The order contains the following attributes: original virtual parking spot, destination, virtual parking spot, distance, and expected travel time.The origins of the orders are first generated stochastically according to a certain distribution. Spatially, the virtual parking spots are aggregated based on the arc, i.e., virtual parking spots on the same arc have the same pickup demand distribution. Temporally, the pickup demand at each time step is assumed to conform to a certain Poisson distribution. Given the origin of an order, the destination is generated from a Rayleigh distribution:fx=x/σ2exp−x2/2σ2, where x is the Manhattan distance from the destination to the origin and σ is a parameter that is set to 3 empirically. Given the origin and the destination, the shortest path between the origin and destination in the road network is calculated as the distance (e.g., calculated by the Dijkstra algorithm or Floyd algorithm). The user’s speed is assumed to be generated from a normal distribution: N12km/h,5km/h2, empirically. Therefore, the expected travel time can be estimated and will be approximated as an integer number of time steps.For each new order, the simulator will first try to dispatch it to an e-bike. If there is an available e-bike in the original virtual parking spot or neighboring virtual parking spots on the same arc, the dispatch will be successful and the order will be tracked by the simulator. The usage state of the dispatched e-bike will be changed to be in use; else, the dispatch is failed and the order will be deleted.In the “return demand” module, at each time step, the simulator will subtract a time step from the expected travel time of the orders in progress. If the expected travel time of an order becomes zero, this order is finished and will be deleted from the simulator. The attributes of the e-bike dispatched to this order will be updated: the location will be changed to the destination, the battery power level (represented in range) will be subtracted by the distance of this order, and the state will be changed to be not in use.Finally, if the time step reachesT, the simulation is ended and the simulation data will be analyzed.


## 7. Numerical Studies

In this section, the results of the numerical study on a test case are shown. In Section7.1, the test case is presented; Section 7.2 introduces the benchmark policies; a comparison of the performance of different policies is presented in Section 7.3; and more insights into the experiment are given in Sections 7.4–7.6. Some parameters of the model, as well as some tuned parameters in the algorithms, are given in Table 1.


### 7.1. Test Case

The simulated scenario of the test case study is shown in Figure9, covering an area of 6∗6=36km2, with 81 vertexes and 272 arcs. Assume that vertex 41 is also the location of the depot (i.e., vertex v0). To make the scenario as realistic as possible, we set three types of districts, indicated by red, blue, and green, respectively. Suppose the red arcs represent the central business district, the blue arcs represent the residential districts, and the green arcs represent the connective district.Figure 9
Road network of the test case study.Arcs in different colored districts are supposed to have different temporal trends of user demand. For example, in residential districts, the pickup demand level for e-bikes is higher in the morning rush than in the evening rush, while in central business districts, the pickup demand level for e-bikes is higher in the evening rush than in the morning rush. The distribution of user demand is assumed to be the same for the time steps during an hour. The relative pickup demand level for virtual parking spots in the three different districts during each hour is shown in Figure10.Figure 10
Relative pickup demand level for virtual parking spots in each district.Given the relative pickup demand level and the total pickup demand number, the specific pickup demands can be generated for each virtual parking spot at each time step. The total pickup demand number in a day has three levels to verify the validity of the proposed policy: 2000 pickup demands/day, 4000 pickup demands/day, and 8000 pickup demands/day. Given the origin of an order, the destination, distance, and the expected travel time are generated from the rules in Section6.There are three lengths of arcs in the road network: 0.5 km, 1 km, and 2 km. Assume that there are 1–3 virtual parking spots randomly distributed in each arc, so the density of virtual parking spots is greater in the residential and central business districts than in the connecting districts because most arcs in the residential and central business districts are shorter than the arcs in the connecting districts. This is usually in line with the reality.The total number of e-bikes in the system (i.e.,Neb) is fixed at 1000. In the initial configuration, the e-bikes are randomly distributed in each virtual parking spot and the battery level is uniformly distributed in the range 0,ψmax, where ψmax is set to 40 km.Assume that there are 4 workers/vehicles and the battery capacity of each vehicle is 30. The speed of the worker when operating in his subzone is assumed to be constant at 10 km/h in the residential districts and central business districts and 12 km/h in the connective district. The speed of the worker when traveling between the depot and his subzone is assumed to be constant at 15 km/h. From this speed parameter, the deadheading cost of each arc can be calculated. The time cost to replace a battery is set to 1 minute, and the time cost to reload a battery at the depot is set to 0.5 minutes.We simulated for 7 days to evaluate each policy. The working hours are from 7 am to 10 pm. Suppose it is possible to make the initial settings of the system generated from the same distribution for every day to test the effectiveness of the dynamic operations.


### 7.2. Benchmarks

To verify the effectiveness of the suggested policy, the following benchmark policies are designed.Random policy: the random policy selects the feasible decisions with equal probability.ϵ−greedy policy: without considering future information, ϵ−greedy policy selects an arc to visit next based on the ratio of the underpowered e-bike number and available e-bike number on each feasible arc (if the available e-bike number is 0, we replace it with a sufficiently small value). Suppose the number of feasible arcs is n. The feasible arc with the maximum ratio is selected with the probability of 1−ϵ+ϵ/n; any other feasible arc is selected with the probability of ϵ/n. Here ϵ is tuned to 0.3.Rollout policy: the lookahead policy with rollout heuristic (denoted as the rollout policy) introduced in Section5.1 is also used as a benchmark. The simulation horizons and the evaluation way of the rollout policy are consistent with the MCTS algorithm.


### 7.3. Results

We compare the suggested lookahead policy with MCTS algorithm (denoted as the MCTS policy) with the benchmarks at multiple demand levels.The simulation run for each feasible decisionIr is tuned to 16 in the rollout policy. The maximum iteration Imax is set to the number of feasible decisions multiplied by 16 in the MCTS policy to ensure that the total simulation run at each decision point is the same for the rollout policy and the MCTS policy.Three demand levels of 2000 pickup demand/day, 4000 pickup demands/day, and 8000 pickup demands/day are tested. The average score of each policy is calculated from the simulation trajectories by equation (1). The result is shown in Figure 11.Figure 11
Evaluation of different policies at three demand levels. (a) Successful pickup number. (b) Number of replaced batteries. (c) Score.
(a)(b)(c)The result shows that the rollout policy and the MCTS policy outperformed the myopic random policy andϵ−greedy policy significantly. Except at the demand level of 2000 pickup demands/day, the MCTS policy outperformed the rollout policy. Notably, at the demand level of 2000 pickup demands/day, the MCTS policy performed slightly worse than the rollout policy and the ϵ−greedy policy performed slightly worse than the random policy. This phenomenon occurred because at a low demand level (with the same e-bike number and worker number), the amount of tasks for each worker is relatively low, resulting in a relatively abundant labor force for these tasks and a low sensitivity of the results to different policies under the assumption that temporary stops of vehicles are not allowed during operating.The number of replaced batteries varies more while the successful pickup number varies less under different policies, resulting in a higher impact of the number of replaced batteries on the final score. In fact, we found that the impact on the success rate of user pickups through battery replacement operation alone is limited because the EBSS also has other problems that can affect the success of user pickups, such as the unbalanced distribution of e-bike supply and user demand, which could not be solved by battery replacement operation.


### 7.4. The Partitioning

For a given scenario, the road network partition way will be adjusted every time period. Take the test case with a demand level of 4000 pickup demands/day for example; the road network partition during some working hours of a day is shown in Figure12.Figure 12
An example of the road network partition result. The blue, yellow, red, and green arcs denote the subzone of workersw1,w2,w3, and w4, respectively. (a) 8 am. (b) 12 am. (c) 6 pm. (d) 10 pm.
(a)(b)(c)(d)In Figure12, the edges with the same color belong to the same subzone. From the result, it could be seen that the spatiotemporal partitioning method is well able to partition the road network according to the topological characteristics and can adjust over the temporal workload characteristics. The results show that each subzone is relatively close in size and the partitioning results in each period are relatively close, which ensures that the change of the subzone does not cause too much confusion to the workers.The average number of batteries replaced by each worker at different demand levels is shown in Table2. In Table 2, the average number of initial underpowered e-bikes in each subzone at the beginning of the operating day is notated as “initial_underpowered,” the average number of new emerged underpowered e-bikes in each subzone during the operating day is notated as “new_underpowered,” and the average number of batteries replaced by each worker is notated as “battery_replaced.” The % completion_rate is defined as the ratio of the “battery_replaced” and the sum of the “initial_underpowered” and “new_underpowered” in the following equation:(15)%completion_rate=battery_replacedinitial_underpowered+new_underpowered∗100%.Table 2
Results of the batteries replaced by each worker.
Workerinitial_underpowerednew_underpoweredbattery_replaced% completion_rate(a) 2000 pickup demands/dayw113.840.647.487.1w225.645.462.487.9w313.434.841.886.7w420.049.459.285.3(b) 4000 pickup demands/dayw114.886.883.682.3w225.294.697.281.1w314.263.863.881.8w422.698.494.478.0(c) 8000 pickup demands/dayw114.0144.2120.075.9w229.0168.8145.073.3w314.0127.0106.075.2w420.4195.6155.471.9As the number of user demands increases, the number of batteries replaced by each worker also increases, but the % completion_rate of each worker decreases. For a given demand level, the number of batteries replaced by each worker is relatively close, which proves the effectiveness of the partitioning algorithm. From the results, we can see that the batteries replaced by workersw2 and w4 are relatively more than the workers w1 and w3. In fact, this is because in the subzone where the road network is sparse and the e-bikes are more sparsely distributed (such as the subzone of w1 and w3 in Figure 12), the worker has to spend more time running on the road, resulting in a lower number of batteries replaced in the same operating time.From the results, we also can see that at each demand level, the % completion_rate of different workers is relatively close, which reflects that the partition is reasonable and allows the labor force to be appropriately distributed. Also, we can find that the worker number is sufficient at the demand level of 2000 pickup demands/day and the demand level of 4000 pickup demands/day because the batteries replaced are more than the new appeared underpowered batteries, but more workers are needed at the demand level of 8000 pickup demands/day as the batteries replaced by each worker are all less than the new appeared underpowered batteries.


### 7.5. Dynamic and Static Operations

To verify the necessity of dynamic operation, we compared the lost user number under dynamic battery replacement operation and under static battery replacement operation. In the static case, no operation is applied to the e-bikes during the daytime, but the initial setting of a day is set to the same as the dynamic case under the assumption that there is a static battery replacement operation at night.The lost users are divided into two categories: lost users because there is no e-bike in the virtual parking spots, denoted as “unbalance lost,” and lost users because the e-bikes in the virtual parking spots are underpowered, denoted as “underpowered lost.”The average number of lost users under static and dynamic operations at three demand levels is shown in Table3. From the results, we can see that the dynamic operation could reduce the total lost users compared to the static case. In addition, the result shows that dynamic battery replacement mainly reduced the number of the “underpowered lost” users but may exacerbate the unbalance of the e-bike supply, resulting in more “unbalance lost” users. One reason for the high “unbalance lost” number in the dynamic operation case is that for an arc with underpowered e-bikes, in the static operation case, the subsequent failed pickups will be recorded as “underpowered lost”; however, in the dynamic operation case, the underpowered e-bikes may be replaced battery by the worker in time and will be picked up by users who arrive first, and the subsequent failed pickups on this arc will be recorded as “unbalance lost.” Therefore, both the battery replacement operation and rebalancing operation are essential for the EBSS.Table 3
Number of lost users under dynamic and static operation.
Demand levelLost unbalancedLost underpoweredLost total(a) Static2000203.4233.8437.24000536.8842.41379.280001117.83037.04154.8(b) Dynamic (MCTS policy)2000266.083.4349.44000759.6305.61065.280001878.61144.43023.0The average number of lost users reduced by the dynamic operation (compared with the static operation) at three demand levels is 87.8, 314.0, and 1131.8, respectively, which suggests that at the low demand level, such as 2000 pickup demands/day, it may not be cost effective to perform dynamic operations or the number of workers can be reduced appropriately due to the low ratio of user demand over e-bike supply.


### 7.6. Simulation Runs

In the MCTS policy and rollout policy, the number of simulation runs is a key factor for estimating the value of feasible decisions. In general, the more the simulation runs, the more accurate the estimation, but it also requires more computational cost.As the MCTS policy is able to allocate more simulation runs to more highly rewarding routes, it requires less computational costs to achieve the same effect as the rollout policy. We compared the results of these two policies for different number of simulation runs at the demand level of 8000 pickup demand/day, as shown in Table4.Table 4
Results under different simulation runs.
Simulation runs/decisionPickupReplaceScore(a) Rollout44958.8519.032839.784958.2521.232919.4164948.6527.233115.0324962.8532.433346.5(b) MCTS44961.4512.832617.484969.2526.433142.1164972.6530.033347.4325021.0541.233832.2Four numbers of simulation runs are tested: 4, 8, 16, and 32 simulation runs per feasible decision. Take the 8 simulation runs for example; if the number of feasible decisions of a decision point is 3, then the total simulation runs at that decision point is 24. The result shows that the more the simulation runs, the higher the final score in both the rollout policy and the MCTS policy. The results of experiments with 8, 16, and 32 simulation runs all show that the MCTS policy performs better than the rollout policy, indicating that the MCTS policy requires only half the simulation runs to obtain the same or even better results than the rollout policy. However, in the experiment with 4 simulation runs, the MCTS policy performs worse than the rollout policy. One reason for this result is that the simulation runs for some feasible decisions may be so inadequate as to cause great deviation in the MCTS policy when the total number of simulations runs is relatively low.


## 7.1. Test Case

The simulated scenario of the test case study is shown in Figure9, covering an area of 6∗6=36km2, with 81 vertexes and 272 arcs. Assume that vertex 41 is also the location of the depot (i.e., vertex v0). To make the scenario as realistic as possible, we set three types of districts, indicated by red, blue, and green, respectively. Suppose the red arcs represent the central business district, the blue arcs represent the residential districts, and the green arcs represent the connective district.Figure 9
Road network of the test case study.Arcs in different colored districts are supposed to have different temporal trends of user demand. For example, in residential districts, the pickup demand level for e-bikes is higher in the morning rush than in the evening rush, while in central business districts, the pickup demand level for e-bikes is higher in the evening rush than in the morning rush. The distribution of user demand is assumed to be the same for the time steps during an hour. The relative pickup demand level for virtual parking spots in the three different districts during each hour is shown in Figure10.Figure 10
Relative pickup demand level for virtual parking spots in each district.Given the relative pickup demand level and the total pickup demand number, the specific pickup demands can be generated for each virtual parking spot at each time step. The total pickup demand number in a day has three levels to verify the validity of the proposed policy: 2000 pickup demands/day, 4000 pickup demands/day, and 8000 pickup demands/day. Given the origin of an order, the destination, distance, and the expected travel time are generated from the rules in Section6.There are three lengths of arcs in the road network: 0.5 km, 1 km, and 2 km. Assume that there are 1–3 virtual parking spots randomly distributed in each arc, so the density of virtual parking spots is greater in the residential and central business districts than in the connecting districts because most arcs in the residential and central business districts are shorter than the arcs in the connecting districts. This is usually in line with the reality.The total number of e-bikes in the system (i.e.,Neb) is fixed at 1000. In the initial configuration, the e-bikes are randomly distributed in each virtual parking spot and the battery level is uniformly distributed in the range 0,ψmax, where ψmax is set to 40 km.Assume that there are 4 workers/vehicles and the battery capacity of each vehicle is 30. The speed of the worker when operating in his subzone is assumed to be constant at 10 km/h in the residential districts and central business districts and 12 km/h in the connective district. The speed of the worker when traveling between the depot and his subzone is assumed to be constant at 15 km/h. From this speed parameter, the deadheading cost of each arc can be calculated. The time cost to replace a battery is set to 1 minute, and the time cost to reload a battery at the depot is set to 0.5 minutes.We simulated for 7 days to evaluate each policy. The working hours are from 7 am to 10 pm. Suppose it is possible to make the initial settings of the system generated from the same distribution for every day to test the effectiveness of the dynamic operations.


## 7.2. Benchmarks

To verify the effectiveness of the suggested policy, the following benchmark policies are designed.Random policy: the random policy selects the feasible decisions with equal probability.ϵ−greedy policy: without considering future information, ϵ−greedy policy selects an arc to visit next based on the ratio of the underpowered e-bike number and available e-bike number on each feasible arc (if the available e-bike number is 0, we replace it with a sufficiently small value). Suppose the number of feasible arcs is n. The feasible arc with the maximum ratio is selected with the probability of 1−ϵ+ϵ/n; any other feasible arc is selected with the probability of ϵ/n. Here ϵ is tuned to 0.3.Rollout policy: the lookahead policy with rollout heuristic (denoted as the rollout policy) introduced in Section5.1 is also used as a benchmark. The simulation horizons and the evaluation way of the rollout policy are consistent with the MCTS algorithm.


## 7.3. Results

We compare the suggested lookahead policy with MCTS algorithm (denoted as the MCTS policy) with the benchmarks at multiple demand levels.The simulation run for each feasible decisionIr is tuned to 16 in the rollout policy. The maximum iteration Imax is set to the number of feasible decisions multiplied by 16 in the MCTS policy to ensure that the total simulation run at each decision point is the same for the rollout policy and the MCTS policy.Three demand levels of 2000 pickup demand/day, 4000 pickup demands/day, and 8000 pickup demands/day are tested. The average score of each policy is calculated from the simulation trajectories by equation (1). The result is shown in Figure 11.Figure 11
Evaluation of different policies at three demand levels. (a) Successful pickup number. (b) Number of replaced batteries. (c) Score.
(a)(b)(c)The result shows that the rollout policy and the MCTS policy outperformed the myopic random policy andϵ−greedy policy significantly. Except at the demand level of 2000 pickup demands/day, the MCTS policy outperformed the rollout policy. Notably, at the demand level of 2000 pickup demands/day, the MCTS policy performed slightly worse than the rollout policy and the ϵ−greedy policy performed slightly worse than the random policy. This phenomenon occurred because at a low demand level (with the same e-bike number and worker number), the amount of tasks for each worker is relatively low, resulting in a relatively abundant labor force for these tasks and a low sensitivity of the results to different policies under the assumption that temporary stops of vehicles are not allowed during operating.The number of replaced batteries varies more while the successful pickup number varies less under different policies, resulting in a higher impact of the number of replaced batteries on the final score. In fact, we found that the impact on the success rate of user pickups through battery replacement operation alone is limited because the EBSS also has other problems that can affect the success of user pickups, such as the unbalanced distribution of e-bike supply and user demand, which could not be solved by battery replacement operation.


## 7.4. The Partitioning

For a given scenario, the road network partition way will be adjusted every time period. Take the test case with a demand level of 4000 pickup demands/day for example; the road network partition during some working hours of a day is shown in Figure12.Figure 12
An example of the road network partition result. The blue, yellow, red, and green arcs denote the subzone of workersw1,w2,w3, and w4, respectively. (a) 8 am. (b) 12 am. (c) 6 pm. (d) 10 pm.
(a)(b)(c)(d)In Figure12, the edges with the same color belong to the same subzone. From the result, it could be seen that the spatiotemporal partitioning method is well able to partition the road network according to the topological characteristics and can adjust over the temporal workload characteristics. The results show that each subzone is relatively close in size and the partitioning results in each period are relatively close, which ensures that the change of the subzone does not cause too much confusion to the workers.The average number of batteries replaced by each worker at different demand levels is shown in Table2. In Table 2, the average number of initial underpowered e-bikes in each subzone at the beginning of the operating day is notated as “initial_underpowered,” the average number of new emerged underpowered e-bikes in each subzone during the operating day is notated as “new_underpowered,” and the average number of batteries replaced by each worker is notated as “battery_replaced.” The % completion_rate is defined as the ratio of the “battery_replaced” and the sum of the “initial_underpowered” and “new_underpowered” in the following equation:(15)%completion_rate=battery_replacedinitial_underpowered+new_underpowered∗100%.Table 2
Results of the batteries replaced by each worker.
Workerinitial_underpowerednew_underpoweredbattery_replaced% completion_rate(a) 2000 pickup demands/dayw113.840.647.487.1w225.645.462.487.9w313.434.841.886.7w420.049.459.285.3(b) 4000 pickup demands/dayw114.886.883.682.3w225.294.697.281.1w314.263.863.881.8w422.698.494.478.0(c) 8000 pickup demands/dayw114.0144.2120.075.9w229.0168.8145.073.3w314.0127.0106.075.2w420.4195.6155.471.9As the number of user demands increases, the number of batteries replaced by each worker also increases, but the % completion_rate of each worker decreases. For a given demand level, the number of batteries replaced by each worker is relatively close, which proves the effectiveness of the partitioning algorithm. From the results, we can see that the batteries replaced by workersw2 and w4 are relatively more than the workers w1 and w3. In fact, this is because in the subzone where the road network is sparse and the e-bikes are more sparsely distributed (such as the subzone of w1 and w3 in Figure 12), the worker has to spend more time running on the road, resulting in a lower number of batteries replaced in the same operating time.From the results, we also can see that at each demand level, the % completion_rate of different workers is relatively close, which reflects that the partition is reasonable and allows the labor force to be appropriately distributed. Also, we can find that the worker number is sufficient at the demand level of 2000 pickup demands/day and the demand level of 4000 pickup demands/day because the batteries replaced are more than the new appeared underpowered batteries, but more workers are needed at the demand level of 8000 pickup demands/day as the batteries replaced by each worker are all less than the new appeared underpowered batteries.


## 7.5. Dynamic and Static Operations

To verify the necessity of dynamic operation, we compared the lost user number under dynamic battery replacement operation and under static battery replacement operation. In the static case, no operation is applied to the e-bikes during the daytime, but the initial setting of a day is set to the same as the dynamic case under the assumption that there is a static battery replacement operation at night.The lost users are divided into two categories: lost users because there is no e-bike in the virtual parking spots, denoted as “unbalance lost,” and lost users because the e-bikes in the virtual parking spots are underpowered, denoted as “underpowered lost.”The average number of lost users under static and dynamic operations at three demand levels is shown in Table3. From the results, we can see that the dynamic operation could reduce the total lost users compared to the static case. In addition, the result shows that dynamic battery replacement mainly reduced the number of the “underpowered lost” users but may exacerbate the unbalance of the e-bike supply, resulting in more “unbalance lost” users. One reason for the high “unbalance lost” number in the dynamic operation case is that for an arc with underpowered e-bikes, in the static operation case, the subsequent failed pickups will be recorded as “underpowered lost”; however, in the dynamic operation case, the underpowered e-bikes may be replaced battery by the worker in time and will be picked up by users who arrive first, and the subsequent failed pickups on this arc will be recorded as “unbalance lost.” Therefore, both the battery replacement operation and rebalancing operation are essential for the EBSS.Table 3
Number of lost users under dynamic and static operation.
Demand levelLost unbalancedLost underpoweredLost total(a) Static2000203.4233.8437.24000536.8842.41379.280001117.83037.04154.8(b) Dynamic (MCTS policy)2000266.083.4349.44000759.6305.61065.280001878.61144.43023.0The average number of lost users reduced by the dynamic operation (compared with the static operation) at three demand levels is 87.8, 314.0, and 1131.8, respectively, which suggests that at the low demand level, such as 2000 pickup demands/day, it may not be cost effective to perform dynamic operations or the number of workers can be reduced appropriately due to the low ratio of user demand over e-bike supply.


## 7.6. Simulation Runs

In the MCTS policy and rollout policy, the number of simulation runs is a key factor for estimating the value of feasible decisions. In general, the more the simulation runs, the more accurate the estimation, but it also requires more computational cost.As the MCTS policy is able to allocate more simulation runs to more highly rewarding routes, it requires less computational costs to achieve the same effect as the rollout policy. We compared the results of these two policies for different number of simulation runs at the demand level of 8000 pickup demand/day, as shown in Table4.Table 4
Results under different simulation runs.
Simulation runs/decisionPickupReplaceScore(a) Rollout44958.8519.032839.784958.2521.232919.4164948.6527.233115.0324962.8532.433346.5(b) MCTS44961.4512.832617.484969.2526.433142.1164972.6530.033347.4325021.0541.233832.2Four numbers of simulation runs are tested: 4, 8, 16, and 32 simulation runs per feasible decision. Take the 8 simulation runs for example; if the number of feasible decisions of a decision point is 3, then the total simulation runs at that decision point is 24. The result shows that the more the simulation runs, the higher the final score in both the rollout policy and the MCTS policy. The results of experiments with 8, 16, and 32 simulation runs all show that the MCTS policy performs better than the rollout policy, indicating that the MCTS policy requires only half the simulation runs to obtain the same or even better results than the rollout policy. However, in the experiment with 4 simulation runs, the MCTS policy performs worse than the rollout policy. One reason for this result is that the simulation runs for some feasible decisions may be so inadequate as to cause great deviation in the MCTS policy when the total number of simulations runs is relatively low.


## 8. Conclusion

In this paper, we have modeled the battery replacement in an EBSS as a DCARP according to the distribution characteristics of the virtual parking spots. The objective is to maximize a weighted sum of the number of successful pickup demands and the number of batteries replaced. A “Partition-First, Route-Second” solution framework was proposed for the scenario with multiple workers. A spatiotemporal partitioning method was proposed to partition the road network into subzones on a large time scale, and the routing decision process of each worker was modeled as MDP on a small time scale. A lookahead policy with MCTS sampling algorithm was proposed to anticipate the potential future influence of the current decisions. The numerical study results show that the spatiotemporal partitioning method was able to avoid the conflict between workers’ operations and allocate the labor force reasonably. The MCTS policy proposed in this paper outperformed the myopic policies and could reduce the simulation runs significantly with no degradation in effectiveness compared to the rollout policy.The limitations of the application scenarios of the suggested policy also need to be noted. Like the rollout policy and other lookahead policies, the MCTS policy relies on the forecast of future exogenous information, such as the user demand trend. As analyzed in literature [25], if the impact of noise of user demand becomes too significant, the anticipation for future demand is less possible, and a policy only based on the current state may yield better results than the lookahead policy.Future research for the EBSS may focus on the case where battery replacement operation and the rebalancing operation are performed simultaneously in the system. The unbalance between supply and demand in bike sharing systems has been studied extensively and also exists in the EBSS. Our experiments show that only replacing the battery dynamically may have limited improvement in the success rate of user pickup. As the e-bike has not only location attributes but also battery power level attributes, the rebalancing problem of an EBSS may be more complicated and needs further study.



---

*Source: 9665340-2021-10-23.xml*