I Introduction
Logistics management involving flow of goods from points of origins to points of destinations, via delivery agents moving in a regular or complex network, is an interesting example of complex systems in which a large number of components interact with each other in strongly nonlinear ways[1][2]. While the rules of interaction are generally quite simple, the dynamics of such systems can be diverse and rather unpredictable, leading to universal emergent behaviours that can be exploited for system optimisation[3]. From a practical point of view, the management of such logistics system mainly involves efficient resource allocation in both spatial and temporal domains, especially when there is a need to respond to goods with origins and destinations generated on a real time and dynamical basis[4][5][6][7].
Taking the taxi system as an example, we have a system consisting of mainly two types of agents  commuters and taxis  interacting over the domain of an urban road network. The main task for the taxis is to deliver commuters from their origins to their respective destinations, and the available resource is the fleet of empty taxis at any moment of the day. The management of such resources is given by the dynamical allocation of the empty taxis over the entire road network, subject to factors such as the connectivity of the road network, the real time traffic conditions, and vehicle travelling speeds. In principle, such management can be effectively implemented by controlling the total number of taxis in the system, as well as the routing of individual taxis
[8]. The routing from origins to destinations for occupied taxis is generally straightforward, normally with the shortest path taken subject to certain constraints (e.g. traffic conditions or tolls, etc.). Highly nontrivial algorithms can be developed if ridesharing is involved, where occupied taxis may also be available to additional commuters based on judicious routematching algorithms[9]. Routing for empty taxis is nontrivial even without ridesharing, with several approaches in customer searching and demand predictions[10][11][12][13][14][15][16]. The optimal routing algorithm also depends on the available communication technologies. While in olden days roadside hailing is the only option for empty taxis to meet commuters, nowadays advanced realtime taxi booking via smartphones can be easily implemented, so that empty taxis have additional information of where the commuters are waiting and where their destinations are, before selecting an efficient routing strategy.There are also other similar logistics systems, consisting of agents interacting in analogous ways. Instead of commuters, picking up and delivering goods efficiently on a large scale is a common challenge for maritime port management, storage management and post delivery, etc. It is becoming more prevalent with the popularity of onlinetooffline services (e.g. grocery delivery) and the advances in automated technologies (e.g. autonomous vehicles and drones). A highly efficient and decentralised delivery network/system could be an important part of the overall logistics management in the near future, especially when the demand is generated in a realtime and less predictable manner. In general there can be many practical complications and constraints for such delivery systems[17], and sophisticated pickup and delivery algorithms could be developed especially in skip delivery problems (SDP) or split delivery vehicle routing problems (SDVRP)[18][19]. In this work, however, we use abstract models to simplify these complications, with the aim of finding universal features resulting from the dynamics of such systems that are unaffected by certain details systematically ignored in our analysis and simulations. Understanding such universal features could be useful in benchmarking and optimising these logistics systems with various resource allocations, routing and delivery algorithms.
Throughout this paper we denote the agents that are delivered from origins to destinations as , or agents; and the delivery agents that moves agents around as , or agents. While this paper does not focus on actual system optimisation strategies themselves, here we give a brief discussion on a number of optimisations the analysis in this paper can be useful for. The optimisation of the number of the agents, as well as the routing algorithm, can depend on a number of factors. In general, we would like to minimise the number of agents needed to save cost. In addition, a good routing algorithm should also minimise the average distance travelled by agents, which can be physically related to the “fuel cost” or the “maintenance cost”. For transportation systems, reducing the number of agents and the number of trips can also help reduce traffic congestions. If rewards are given for every agent delivered, we can also try to maximize the total or average earnings of agents. From the perspective of agents, the most important factor is the average waiting time; for every agent generated in the system, we would like an agent to pick it up as soon as possible. If pooling or ridesharing is involved, the trip time for a agent to be with an agent not only depends on the intrinsic origin and destination for agents, but also on possible detours for additional pickup/dropoff of the shared parties. It is thus also sensible to minimise the total travel time (which is the sum of the waiting time and the trip time), or the total cost of the trip (i.e. for taking a taxi). The overall utility function of the logistics system can be one or a combination of a number of factors mentioned above.
In this paper, we do not focus on how to optimize such logistics systems from the perspectives mentioned in the previous paragraph, as this will be discussed in detail elsewhere[20]. Instead, we give a welldefined mathematical formulation of the dynamics of such complex systems, followed by employing both analytical and numerical tools to understand and formulate universal and important features of the dynamics of such systems. These features will be useful for benchmarking system efficiency, predicting qualitatively different behaviours of the agents, as well as guiding specific optimisations of these logistics systems. Given the nonlinear interactions between agents over the complex network, the challenge is to solve such systems with limited analytical tools. Simplification of the mathematical formulation by ignoring some detailed agent behaviours, combined with macroscopic phenomenological reformulation of the problem where analytical treatment is possible, can lead to useful approximate solutions that captures the essential dynamics. It is, however, indispensable to have microscopic agent based simulations that are as realistic as possible, to evaluate and validate these approximate solutions. Such microscopic simulations are also very useful in exploring general behaviours of the system dynamics under various different scenarios.
This paper will be organised as follows: in Sec. II we formulate a network based logistics system with two types of agents in rigorous mathematical languages; in Sec. III we study mathematical models with simple settings but nontrivial behaviours, that capture essential features of the dynamics of the networkbased logistics systems with two types of agents; in Sec. IV we study numerically a particular example of the taxi systems, revealing various useful dynamical features that can be understood via the mathematical analysis in Sec. II; in Sec. V we focus on the two dynamical phases of the taxi systems, showing that even with the additional dynamics involving ridesharing, the formulation of different phases from formal and numerical analysis of the biagent logistics system can be useful for predicting complicated behavours affecting the quality of the taxi services; in Sec. VI we summerize our work and discuss about the outlooks.
Ii Mathematical Formulation
At a more abstract level, the logistics systems we focus on is a biagent system consisting of two types of agents, with nontrivial dynamics on a potentially complex network. We represent the network as a graph , where is the collection of nodes of the graph, and is the collection of directed, weighted edges connecting the nodes. The graph can also be fully represented by an adjacency matrix . For nodes , is nonzero if there is an edge directed from to . This is the case if and only if an agent can move from to without passing through any other nodes (so the node is directly connected from node). The matrix element gives the inverse of the time it takes for an agent to move from to . Thus can be timedependent, and if different agents have different velocities, each agent will have its own corresponding adjacency matrix. implies no edge between the two nodes, as the time it takes to travel from to is infinity.
We denote the two types of agents as and ; the corresponding rate of generation of the agents at node is given by and respectively, and both in principle can be time dependent. One or both agents can move from one node to another via edges, and when an agent and an agent (note the upper indices give the agent index, while the lower indices give the node index) meet at the same node, they “annihilate” each other in pairs, and leave the system. Let and be the number of agents of and respectively, the most important dynamics of such logistics systems is the time dependence of these two quantities, denoted by and .
The motion of the agents are described by an policy matrix , where is the index of the agents. For an agent at node
, its probability of moving towards a neighbouring node
at the next time step is given by ; thus we have , and if and only if . In principle, each agent can have its own policy matrix. If the motion of the agents are controlled in a centralized manner, then the policy matrix can differ from one agent to another based on the control algorithm, and can also be time dependent. On the other hand, for decentralized systems in which agents make their own decisions, all agents of the same type ( or ) can follow the same policy matrix. If that is the case, we omit the upper indices and simply denote the policy matrix as . It is also common for agents of the same type to follow the same policy matrix under certain conditions, and individualized policy matrices under other conditions (e.g. the implementation of booking policies in taxi systems, as we will discuss later on).The dynamics of this interacting system is thus completely defined by , and the annihilation rule between the two types of agents. We would like to emphasize here that for realistic systems, and may not be independent. As we shall see in Sec. IIA, the “annihilation” between agents from and can also correspond to the formation of “bound states” with nontrivial dynamics, with the “bound state” delivering the bound agent in it from its origin to its destination. At the destination the agent leaves the system, while the agent in the “bound state” reenters the system. If the latter is the only physical way for agents from to be generated in the system, then depends on and the origin/destination distribution of . Theoretically, it is useful to treat the dynamics of the “bound states” as hidden, and as dynamical and tunable. The characteristic dynamical behaviours can be systematically studied both analytically and numerically, as we will show in Sec. III and Sec. IV.
Iia Taxi System as a Special Case
In a taxi system, the network corresponds to the road network, and the two types of agents are the empty taxis () and the commuters (). The nodes are locations where commuters can board or alight the taxis, and the “annihilation” process corresponds to a commuter being picked up by an empty taxi when they meet at the same node. The commuters are stochastically generated at different nodes in the road network, with probability per time step either obtained from historical empirical data, or artificially synthesized. When a commuter is picked up by an empty taxi, both agents will disappear from the system. Empty taxis will reemerge at different nodes according to , corresponding to the roaming of the empty taxis and the spatiotemporal distribution of the destinations of the commuters, where commuters alight (but not reenter the system) and the taxis become empty again.
In this particular system, we treat commuters as nonmobile: once they are generated, they stay at the nodes with a trivial policy matrix . The empty taxis, on the other hand, are mobile with a nontrivial policy matrix. The efficiency of the taxi system depends strongly on how good this policy matrix is, corresponding to the strategy of empty taxis in anticipating where potential commuters will be in the road network.
The taxi system is also a typical example with nontrivial “hidden dynamics” of “bound states” representing an occupied taxis with a boarded commuter. In general, an occupied taxi will go along the shortest path from the origin to the destination of the passenger. Thus, part of depends on factors including and the distribution of the destinations, which we denote as . Physically, is the probability that a commuter boarding at node will alight at node , so we have . In our theoretical treatment in Sec. III, we ignore the dynamics of the occupied taxis, and treat the emergent as a phenomenological input. For microscopic agent based simulations in Sec. IV, on the other hand, the dynamics of both empty and occupied taxis are fully accounted for, and meaningful comparisons between theory and simulations can be made.
We also would like to comment here that the dynamics of the occupied taxis is generally quite straightforward if no ridesharing is involved, as is the case in the most part of this paper. Realtime adaptive ridesharing is a very important area of research, where the routing of occupied taxis and the route matching of different commuters can be highly nontrivial (see ref.[3] and the references therein). We will analyze this more complex system in details elsewhere; a number of theoretical and numerical results in this work will also lay foundation to a systematic characterization of taxis with dynamical ridesharing, as we will illustrate in Sec. V.


Notations  Physical meaning 
Rate of generation of agent at node  
Rate of generation of agent at node  
The OD matrix, the matrix element gives the probability  
for the commuter generated at node i to have  
destination at node  
Weighted, directed adjacency matrix for the network,  
the matrix element gives the inverse of travel time across the edge  
Policy matrix, the matrix element gives the probability  
for the agent at node to move to node 
Iii Theoretical Modelling and Analysis
Given the adjacency matrix of the network , the spatiotemporal distribution of the rate of generation of the two agents given by , and the policy matrix governing the movements of the agents, we would like to analyse both the equilibrium and nonequilibrium dynamics of the agents in the system. In particular, one can calculate many useful quantities (e.g. average waiting time, average travel distances, etc.) from and , the number of the two types of agents as a function of time.
It is useful to first look at the simple model of a single node with fixed probabilities and , where for each time step is the probability of generating a agent, while is the probability of generating an agent. The latter is equivalent to removing a agent at the node if and only if there is one present, since a agent and an agent annihilate each other. We also assume that agents of type can accumulate at the node, while agents of type do not; thus if an agent is generated at the node with no agents present, this agent will be removed at the end of the time step. This corresponds to the case that in the full network, agents do not move, while agents roams in the network with some policy matrix.
One should also note that given our system is discrete in time domain, both and (as well as and in the rest of this paper) are the average number of agents of and generated in , where is the time interval between consecutive time steps. Thus and scale with , and there is thus no constraint that and cannot be larger than . For or , implementations in numerical simulations can also be carried out straightforwardly by rescaling , while keeping the physically relevant probability densities constant. Alternatively, for small systems one can employ the Gillespie algorithm[24]
to simulate the continuoustime Markov chain directly.
Iiia The Single Node Model
We now look at the simple model with only one node, with the associated and both satisfying , which can always be achieved by rescaling . Similar problems with the constraint that has been studied before[21]. Our analysis here is more general, and we can view the dynamics of the accumulated waiting time as a Markov chain with state space . More concretely, we have the number of agents at the node at time t, . We thus have the initial condition and the recursive relation as follows:
(1) 
where we have defined , the probability of having agents waiting at the node at time step . Of course, we have the initial condition
. By computing the eigenvalues and eigenvectors associated with the transition probability matrix, we can derive the following explicit form of
:(2) 
where
(3) 
and
(4) 
Furthermore, by performing asymptotic analysis for
, we find that(5) 
where is the statistical average of the number of waiting agents at the node at time
. This signals a secondorder phase transition at the line
.An important quantity in this system is the average waiting time of the agents, which is averaged over the total number of agents being generated over the simulation time . Its formal expression is given by the following:
(6) 
To carry out nonasymptotic analysis, we write down the explicit expression for
(7) 
with given by Eq. (2). In particular, for the case , by summing the series in explicitly, we get
(8) 
In Fig.(1) we plotted Eq.(8) for different values of and , in the regime that . As one can see clearly, for small and , we have slow convergence of to its equilibrium value. This will be reflected in the numerical simulations in Sec. IV, and could be important for certain logistics systems in actual practice.
IiiB The Noninteracting Limit
The analysis above is completely general, as any system can be rescaled to satisfy by redefining the time step interval . However, either physically or numerically, we can also encounter situations in which at smallest possible both and can be quite large (as we would see in Sec. IV). In such cases, the dynamics of the system can be understood more straightforwardly in the noninteracting limit.
The analysis is nontrivial when because queuing of
agents will invariably occur at some, if not most, of the time steps. Heuristically speaking, queuing
agents interact with each other, as not every one of them will see a pickup rate of . Such interaction requires detailed analytic treatment as we presented in Sec. IIIA.Let us now look at a single agent generated at the node. Assuming no further agents are generated, the probability of this agent surviving for exactly time steps is given by . The average waiting time of this agent is thus given by
(9) 
Eq.(9) is obviously only valid for , and in the case that agents are continuously generated with probability at each time step, the average waiting time for all agent is no longer given by Eq.(9). This is because when there is queuing at the node, the effective probability of the generation of agent, , for the queued agents are reduced, due to this rather nontrivial interaction between existing agents at the node.
The noninteracting limit can be realized with and . In this case, each time step there is always one agent at the node. Given that , agents will not interact with each other as there will be no queuing. In principle, for any we can rescale , the time interval between two consecutive time steps, to . We thus map to the case of , Eq. (9) can be applied for all agent in the system, with the average waiting time as “zero”. Note that physically this means the average waiting time is smaller than , which will not be resolved in numerical simulations in fixed . Thus Eq.(9) is only useful for large rate of generation of agents; otherwise for very small rate, needs to be large for the noninteracting limit to apply, and the conclusions from the noninteracting limit is not very useful. For example, if is one hour, only knowing that the average waiting time is less than one hour from the noninteracting limit cannot help us optimise it.
This noninteracting limit is thus usually the case for realistic logistics systems where is generally large even with small time resolution. The dynamics of such systems can be more straightforwardly understood in the formalism presented here. The simplification comes from the fact that the detailed dynamics of the system within the time resolution are ignored. Such dynamics are still fully captured in the more general analysis in Sec. IIIA, though they may not be physically relevant when the noninteracting limit is applicable. As we will see in Sec. IV, for taxi systems with booking policies, the noninteracting limit generally applies. In contrast, with stochastic policies corresponding to roadside hailing, the more general analysis is needed to explain certain features of the taxi dynamics.
IiiC The Network Effects
The generalisation of the single node mode to the entire network is straightforward if we know the spatial distribution of the rate of generation of agent, (where is the node index). It is, however, rather nontrivial to determine this distribution, which depends on many factors of the transport system. There are two parts contributing to : one part comes from the intrinsic generation of new agents, the other part comes from moving of agents from one node to another. For the latter, it depends both on the underlying network and the policy matrix governing the movements of agents. Microscopically, the equation governing the dynamics of is given as follows:
(10)  
Here can be viewed as the number of empty taxis at the node at time . The policy matrix and the adjacency matrix are both given and static; is the source term describing the intrinsic generation of new agents at the beginning of the time step, while is the sink term describing the intrinsic removal of existing agents at the end of the time step. Both and depends on the detailed dynamics of the logistics system. Eq.(10) is not equivalent to diffusion or random walk on graphs[22] where the edge weights are interpreted as transition probabilities from one node to another. Here, the agents have to travel from one node to another, and the edge weights are inversely proportional to the time it takes for the agent to move between neighbouring nodes. While on average we would expect Eq.(10) can be effectively described as some diffusion processes, exactly solving Eq.(10) is difficult. It does allow us to numerically study the time evolution of without doing the full scale agent based modelling, which we will explore in Sec. IV.
Iv Microscopic Numerical Simulations
We now carry out large scale agent based numerical simulations of more realistic systems, to understand their dynamics based on the theories developed in the previous sections. Our large scale simulation platform consists of (1) an underlying network, (2) the agents of with their origins and destinations denoted by the nodes in the network, and (3) the agents of that move along the edges of the network from one node to another. The positions and status of all agents are updated at every time step, which we can nominally represent as one second. The agents are generated stochastically according to ; thus at each time step, a “dice” is thrown for each of the nodes to determine if a new agent is generated at this particular node. The generation of agents can be either done similarly with , or via any specific “hidden dynamics” that we can now fully capture with microscopic simulations. For example, an agent can roam the network according to some predetermined policy matrix (the number of time steps for an agent to move from node to node is given by ), bind with an agent if they happen to be at the same node at a specific time step, delivering it to its destination node and become a free agent again.
As a concrete example, we carry out detailed microscopic simulation of the taxi systems as described in Sec. IIA, and compare the simulation results with the analysis in Sec. III. We focus on the cases where the road network
and the commuter generation probability distribution
are both timeindependent. The empty taxi generation probability is time dependent and governed by the dynamics of taxis picking up commuters and dropping them off at the destinations, in addition to empty taxis roaming about the road network. We thus require the input of an origin/destination matrix and the policy matrix of the empty taxis as introduced in Sec. IIA, which we also take to be time independent. For each simulation, the total number of taxis (the sum of empty taxis and occupied taxis) is also conserved and time independent. We would like to emphasize that for occupied taxis (which is not treated as agents in our theoretical formulation), they move from the origin of the boarded commuter to the respective destination via the shortest path computed from the road network using the Dijkstra’s algorithm.Iva Stochastic Policies and Booking Policies
Given a specific road network and the commuter demand patterns (as encoded in ), we would expect the dynamics of the taxi system to be strongly dependent on the movement of the empty taxis governed by . Two types of the policies are considered in this work. The first type is the completely decentralised stochastic policies, where is the same for every empty taxi looking for commuters. This type corresponds to the traditional roadside hailing of taxis, in which empty taxis do not know where exactly commuters are waiting. At each node, gives the probability distribution of which neighbouring node the empty taxi will move to. The distribution itself can be random or highly nontrivial. Instead of being completely random, experienced taxi drivers know where to look for commuters with a more effective .
The second type is based on realtime booking of taxis, nowadays quite common via smartphone Apps or telephone calls. In this case, when a commuter is generated at any node in the network, a certain algorithm is implemented to assign a nearby empty taxi to come pick the commuter up. The assigned empty taxi becomes booked, which will go from its current position to the location of the commuter via the shortest path. In this way, the policy matrix is different for different booked taxis (which are still counted as agents of type ), which we formally denote as . In this work, we use a simple assignment algorithm, in which the nearest empty within a specific range from the commuter will be assigned. For empty taxis that are not booked, they still move with a decentralized policy , and in our simulations it is a completely random policy. Thus empty taxis roam the road network with a random walk, before being booked and move directly towards the origin of the booking commuter.
IvB Spatial Generation of Empty Taxis
The spatiotemporal generation of empty taxis, denoted as , is highly important from resource allocation point of view. Given that the total number of taxis is fixed, depends strongly on , the road network structure , and the policy matrix . An accurate prediction of the empty taxi generation pattern is difficult to solve analytically from Eq.(10); we thus employ a practical approach for most of the cases, and study the pattern of numerically. Decentralized movement (i.e. roadside hailing) and taxi booking will be considered separately; and since we only look at cases where are static, we will only focus on the spatial distribution of after the system reaches equilibrium.
For the case of roadside hailing, has two contributions at the node: , the empty taxi generated when an occupied taxi delivered a commuter at its destination at node ; , coming from an empty taxi moving to node from a neighbouring node. In particular, is equivalent to the source term in Eq.(10
). In the special case where the destination is uniformly distributed, so that
, the inverse of the total number of nodes , is uniform and independent of . When the number of taxis is small so that in the long time limit all nodes have queuing commuters, is negligible, so that is also uniformly distributed. This corresponds to the case in Eq.(10) where at any time step. When the number of taxis is large and is sufficiently large, both and are only uniform if there is full translational symmetry of the road network and the policy matrix: is constant and independent of , gives a regular lattice with all edges of the same weight, and , where is the number of neighbours of the node. In particular, if is an arbitrary network, does not generally equilibrate as one can see in Fig. (2).In the case that is nonuniform, on the other hand, is dominated by the spatial distribution of , and only has a subleading effect. It is useful to look at the dependence of on , which we plotted in Fig.(3). While depends on only, also strongly depends on , in addition to . In particular, from a heuristic point of view a good will lead to higher if is larger: if the node has a higher probability of generating commuters per time step, the effective arrival rate of empty taxis at this node should also be higher. In Fig.(3a) the two plots are from a random walk implemented with , and a more intelligent
from the recursive value model enhanced with reinforcement learning
[23]. One can clearly see for the latter, increases appreciably with , indicating more efficient policy matrix for commuter seeking.With taxibooking, we look at the case where the commitment of booking is binding: an empty taxi that is booked can only pick up the respective commuter, and a commuter who books a taxi will only board the assigned taxi. In this way, only has contribution from the booked empty taxis arriving at the locations of their respective commuters. The approximately linear relationship between and is evident from the numerical simulations, as one can see from Fig.(3b). One should also note that in contrast to the stochastic policies where the interaction of the commuters and empty taxis are strictly local (confining to the same node), with taxibooking such interaction is intrinsically nonlocal, depending on the range of booking . The range of booking is defined as such: only empty taxis within a radius of R from the commuter has a chance to be booked. We will discuss about such distinction about the locality of interaction in Sec. V.
IvC Critical Taxi Number
Since each node in the network has a welldefined rate of generation of (for commuters) and (for empty taxis), we can use the results from the models in Sec. III to analyze if there will be queuing for each node, and the overall average waiting time of the commuters of the entire system. In particular, from Sec. III we know that for each node there is a secondorder phase transition at . Heuristically speaking, in the long run there will be accumulation of commuters at the node if , and no accumulation if . For the entire road network, there will be no commuter accumulation if for each node; otherwise if we assume commuters will wait for taxis no matter how long it takes, the average waiting time of the commuters will diverge. This generic argument should apply to a wide variety of taxi systems, including various forms of roadside hailing and taxibooking. One should note that in this case, average waiting time of the commuters is equivalent to the average number of commuters in the road network, over all the time steps. Formally, the average waiting time is given by
(11) 
where , is the total number of time steps for simulation, and is the number of commuters (the agents) in the entire road network at time step . For a single node it reduces to Eq.(7).
Given the complicated dependence of on various aspects of the taxi system (especially the details of the road network and the policy matrix), we explore numerically how the average waiting time of the taxi system depends on the number of taxis (sum of empty and occupied taxis) in the simulation. We expect should increase with if all other aspects are fixed. Thus for small , we expect for majority of nodes, while for large we have for majority of the nodes.
For stochastic policies, we let empty taxis roam randomly in the road network. Different stochastic policies lead to qualitatively similar results, which we will discuss in details elsewhere[20]. It is clear from the numerical simulation that when the number of taxis is large enough, given a spatially uniform rate of generation (here we look at the simple case of , where is the total number of nodes of the network, and one commuter is generated in the entire road network per time step), eventually for every node . In this regime, the average waiting time no longer grows with the simulation time, as there is no queuing of commuters at any node. Conversely, as decreases with decreasing , queuing will eventually occur at nodes with . In this case, the average waiting time will increase with , as one can see in Fig.(4). If both and are uniform, a sharp transition can be observed in the limit of large , as one can see the trend in Fig.(4 a). For nonuniform and and finite , there will be a much smoother transition area where queuing only occurs at some nodes. Thus for actual finite simulations, the case with nonuniform will lead to a much smoother dependence of on , as one can see in Fig.(4 b).
While the dynamics of the stochastic policies can be generally understood via the analysis in Sec. IIIA, the booking policies in most cases correspond to the noninteracting limit as we discussed in Sec. IIIB. The general idea of the booking policy is that when a commuter is generated, an optimal taxi is chosen and given the information about the location of this commuters. This taxi then becomes booked and heads directly towards the respective commuter. In contrast to the stochastic policies where the interaction between taxis and commuters is localized (the empty taxi has no information about the location of the waiting commuters until it is at the same node of the waiting commuters), such interaction is global with the booking policies, when the information of the location of the waiting commuter is made available to the booked empty taxi that potentially can be quite far away, depending on the range of the booking policy, .
To understand that, let us look at the simple example, where the nearest empty taxi (if available) at a particular time step is booked for all commuters generated at that time step, no matter how far away this empty taxi is (i.e. ). The trip time for the booked taxi to arrive at the waiting commuter is bounded given the finite size of the network; its statistical average is welldefined and can be computed from the road network and the density of the empty taxis. Apart from this contribution to the average waiting time, we can thus treat the entire network as a single “node”, with the rate of generation of commuters and empty taxis given by and respectively. The noninteracting limit thus applies if (since in all simulations we take as one second between each time step), which is generally the case for any realistically large road network. In Fig.(4 c,d), we take . We thus observe a much sharper transition even for small in Fig.(4 c), and the dynamics is not affected by the distribution of because of the global nature of the booking policies. The transition occurs at , and when there is no queuing at any nodes (every commuter generated can book an empty taxi almost right away). Such qualitative behaviours should apply to other booking policies, i.e. the more realistic cases with finite .
There is thus a welldefined critical number of taxis in the limit of , such that when there are no queuing at any node, and when there are queuing of commuters at some nodes in the road network. For the booking policies with , occurs when , and a rather sharp transition occurs at as one can see from Fig.(4c,d), where . For booking policies with finite range , or with stochastic policies, needs to be large enough so that for , for all . For finite , the transition at may not be sharp because of the slow convergence of the taxi dynamics as one can see from Eq.(8) and Fig.(1). When or are nonuniform, we also expect only a portion of nodes have queuing commuters for a range of . Nevertheless, the transition will be sharp in the limit , and is still welldefined. In particular, for stochastic policies with random walk, as one can see from Fig.(4a), with the same road network and the same total rate of generation .
Interestingly, is smaller in the case of the stochastic policies, as compared to the booking policy with no range constraints. In general, quantitatively depends on all details of the taxi system, and it is an important quantity characterising the efficiency of the taxi system. We also expect this critical number of delivering agent to be a universal concept the dynamics of the biagent logistics system in an arbitrary network. It serves as the boundary of two distinctive dynamical phases of such logistics system, as we will illustrate below in Sec. V.
V Dynamical Phases of the Taxi System
The qualitatively different dependence of the average waiting time on the simulation time for and allows us to define two distinct dynamical phases of the taxi system, with serving as the phase boundary. We denote the phase with the oversaturated phase, where the demand for taxis exceeds the supply, leading to queuing of commuters at some of the nodes in the road network. On the other hand, the phase with is denoted as the undersaturated phase, where the supply of taxis exceeds the demand. In this case, no nodes have queuing commuters in the long run, and there is an excess of empty taxis running around looking for commuters. From the logistics point of view, given a specific demand and a specific road network, a supply of is optimal, in the sense that all taxis are efficiently utilised, and all commuters can be picked up within reasonable amount of waiting time. When the supply of taxis decreases from , the average waiting time of the commuters increases rapidly, suggesting a severe degrading of the service quality of the taxi system. In contrast, when the supply of taxis increases from , the average waiting time of the commuters decreases insignificantly, indicating very small marginal gain. On the other hand, the cost of supplying taxis increases linearly (in the form of fuel cost and/or compensation to drivers). It is thus undesirable to deviate too much from on both sides, and it is important to know that given a particular taxi system, what is the optimal number of taxis of this system.
Implications in optimisation– Once we have as much details as possible (e.g. traffic conditions, commuter loading and unloading time, etc.), comprehensive numerical simulation can lead to accurate prediction of given a specific demand, as well as the corresponding average waiting time at . This has important implications on how we should optimise the average waiting time. While conventionally, the optimal number of taxis in a city is determined by making sure that the average waiting time is below a certain chosen benchmark (e.g. ). We now know that if , we can easily increase the number of taxis in the oversaturated phase to reduce the average waiting time to . On the other hand, if , we need to increase the number of taxis dramatically in the undersaturated phase to achieve , which may not be the best option. In the latter case, we should focus on improving the empty taxi routing policies or booking algorithms, or other approaches to shift , which could be a more economical way in reducing the average waiting time, in contrast to the increase in the total number of taxis in a city.
Va Tuning of
The critical taxi number most directly depends on the rate of generation of commuters . When the total demand increases, more taxis are needed if all other factors are kept constant. With the same total demand given by , different spatial distribution (given by the dependence of on , the node index), the destination distribution (given by the OD matrix ), the underlying road network (given by the adjacency matrix ), and the manoeuvring strategies of empty taxis (given by the policy matrix ) can also lead to variation of . In most cases, the quantitative analysis can only be done numerically, as the effective strongly depends on those factors in highly nontrivial ways.
In this part we analyse the dependence of the dimensionless regularized waiting time on the number of taxis in the system. In the undersaturated phase , we have ; while a second order phase transition occurs at , and in the oversaturated phase with , . Such phase transition is most easily detected numerically with the booking policy. This is because as we have shown in Sec. IIIB, equilibrium can be quickly reached for and in the noninteracting limit. Different booking policies can affect , as we show here in Fig.(5).
While all factors mentioned above can influence with the booking policies, they do not do so independently. We identify two factors that fundamentally affects the critical number of taxis, when the total average demand given by is kept constant. The first factor is the booking policy . We show in Fig.(5a) that reducing the range of booking (so that only empties within a distance of from a commuter can be booked, and the nearest empty taxi, if available, will be booked) generally reduces , making the booking process more efficient in terms of the number of taxis needed for the entire taxi system. Intuitively, this is because if an empty taxi far away from the commuter is booked, it is committed and unable to take other commuters along its way to pick up the designated commuter, leading to suboptimal matching of available taxis and waiting commuters. The tradeoff is that in the undersaturated phase, the average waiting time increases with decreasing (hard to distinguish from Fig.(5a) due to the scale but the change is unambiguous). Thus when there is an oversupply of empty taxis roaming on the streets, we can reduce average waiting time by increasing the range of booking, provided such alteration does not push the taxi system from the undersaturated phase to the oversaturated phase.
The other fundamental factor is the average distances between origindestination pairs generated in the system. This factor depends collectively on and . The greater the average distances, the longer it takes on average for occupied taxis to send commuters from their origins to destinations, which is the hidden part of the taxi dynamics. If the total number of taxis is fixed, heuristically this will lead to lower effective on average, resulting in larger . This is indeed the case as one can see from Fig.(5b), there is an almost linear relationship between and the average trip distance for various different cases. Both artificial lattice road network and empirical islandwide road network of the city state of Singapore are included in the plots. For the Singapore road network, we also compare the synthetic and empirical commuter demand patterns (encoded in and ). We thus have strong numerical evidence that the linear relationship is quite general for different types of road networks and taxi demand patterns.
With booking policies, the movement of booked and occupied taxis are determined by the origins or destinations of respective commuters via the shortest path algorithm. For empty taxis that are not yet booked, their movement is completely decentralised, and in Fig.(5) we choose the simple random walk policy. More intelligent policies for empty, unbooked taxis can also be implemented; they do not alter the qualitative features illustrated in this work, and we will discuss in more details elsewhere. On the other hand, if we take the limit of the range of booking , the booking policy is reduced to the stochastic policy, where all empty taxis will only pick up commuters if they are at the same node of the waiting commuters. We thus expect many of the features on the dependence of on various factors of the taxi system to be similar with the stochastic policies.
In general for realistic large road network, if physically each time step in the simulation represents one second, is very small for almost all nodes. Empirically for the city state of Singapore, on a typical day there are around six commuters generated for every second for the entire road network of nodes, so that [3]. Mathematically, is only welldefined in the limit of the simulation time , when the taxi system reaches equilibrium for and at each node. As we have shown in Sec. III, convergence to asymptotic values of average waiting time can be slow even for a single node, when both and are small. Thus, theoretical investigation of the dependence of requires much longer simulation time, and extrapolation to infinite simulation time needs to be performed, for stochastic policies. Indeed, for a wide number of taxi systems, the behaviours of agrees with that with the booking policies in the limit of the booking range .
VB Phase characteristics in the context of ridesharing
The phase transition in this biagent logistics system has wider applications in more complicated and realistic systems. We illustrate this with the example of the taxi systems in which realtime, dynamic ridesharing between commuters can be implemented. In this system, if the occupant of the taxi and the waiting commuter are both willing to share their rides, and if sharing of their trips will result in detours within the acceptable limits, ridesharing will be implemented. Strictly speaking, there are four types of active agents in this system: empty taxis, occupied taxis with with room to share and the passenger open to ridesharing, commuters open to ridesharing, and commuters not open to ridesharing. For occupied taxis with no room to share or with passengers not open to ridesharing, they are still part of the “hidden dynamics” that is only responsible for the effective generation of other active agents. When no commuters are willing to share, the system reduces to the original biagent system we studied in details in this paper.
Given the nature of realtime, dynamic ridesharing, only booking policies are natural for the available taxis to pick up commuters. An important parameter to tune is , the percentage of commuters generated in the system who are open to ridesharing. This is a phenomenological parameter describing the degree of acceptance of ridesharing for a particular society, and may have important cultural, social and economical implications. When , the theoretical analysis in Sec. III does not strictly apply. Heuristically, however, for each node there is still an effective from both empty taxis and taxis open to ridesharing, and an effective from both types of commuters. This is confirmed numerically, and both effective and now also depend on , the tuning of which can lead to interesting phase transitions with significant impact on the dynamics of the taxi system.
Large scale numerical simulations are carried out with a specific route matching algorithm that ensures when ridesharing occurs, the extra time involved in the detours will not exceed five minutes for any shared parties. We will not go into the details of the simulation, which can be found in Ref[3]. We analyse the dependence of the average waiting time (the time that commuters spend waiting at the nodes) and the average trip time (the time commuters spend inside the taxi) on , the percentage of commuters open to ridesharing. Note that even when a commuter is open to ridesharing, his/her trip may not be shared with another commuter. This is because at the time such a commuter makes a booking request, there may not be suitable partners to share the trip that satisfies the constraint of the detour time.
As shown in Fig.(6), in such systems we also have very welldefined , which decreases with increasing . This confirms the intuition that ridesharing makes the taxi system more efficient. In the oversaturated phase, the average waiting time decreases rapidly with increasing ; while in the undersaturated phase, the decrease is much more marginal. In both phases, on the other hand, the average trip time increases with . This is because with more commuters open to ridesharing, more detours will be executed, leading to longer trip time. Interestingly, when the number of taxis is large enough, increasing can lead to a phase transition, as one can see from Fig.(6) and Fig.(7c,d). This phase transition leads to a significant drop both in terms of the average waiting time and the average trip time. This is because the phase transition physically corresponds to the transition from the case with no empty taxis roaming, to the case with empty taxis roaming. With the latter case, waiting time will be significantly shortened with the booking policy as we have shown in Sec. IV. With empty taxis roaming, the number of shared trips will also reduce significantly (as in the oversaturated phase with no empty taxis, occupied taxis are the only options for many commuters), leading to fewer detours and lower average trip time.
Vi Discussions and Conclusions
In conclusion, this paper mainly focuses on theoretical studies of the general dynamics of logistics systems consisting of two types of agents, which we can generically labelled as the “delivery” agents ( agents) and the “goods” agents ( agents). Both agents are generated in the system on a realtime, dynamical basis, and the tasks of the agents are to find agent and deliver them to their respective destinations. We use the taxi system as the main example, whereby the agents are empty taxis, and the agents are commuters. The interaction between the two types of agents are formulated in precise mathematical languages, and the domain of the interaction is over a potentially complex network (e.g. the road network for taxis and commuters). We constructed simple models that capture the essential features of such logistics systems, and solved them analytically, revealing a second order phase transition when the system reaches equilibrium. Large scale agentbased modelling is also carried out for realistic settings and with additional details, revealing a number interesting emergent behaviours that can be well explained with our simple models and the analytical results.
One of the main results is the identification of a second order transition from the oversaturated phase (where supply of agents exceeds that of the agents) and the undersaturated phase (vice versa). The phase boundary is both mathematically welldefined and easily observable in more realistic, finite systems via agentbased simulations. The dynamical behaviours of the logistics systems are qualitatively different in these two phases, and we studied in details on how the phase boundary can be affected by various components of the logistics system. Moreover, the phase boundary gives a proper definition of an optimal number of agents in the system, from the system efficiency point of view. In contrast to some similar definitions in the literature, this optimal number depends entirely on the dynamics over the spatiotemporal distribution of the rate of generation of and agents. It can be welldefined without explicit reference to any specific utility functions we construct over the logistics system. Given its generality, we conjecture that sharp phase transitions can be detected with most of the commonly used utility functions across the same . We illustrated this with the taxi systems allowing ridesharing, where we considered average waiting time and travel time for the commuters, and average trip time which can corresponds to the fuel costs for the drivers.
Going forward, we will give a more detailed analysis on how to apply the theoretical tools and conceptual development in this work for benchmarking various optimisation schemes for this class of logistics systems, especially when considering different utility functions from customers’, operators’ and policymakers’ perspectives. Complex logistics systems as the ones discussed in the work can have many components to tune, in order for certain aspects of performances to improve. The theoretical treatment in this work allows us to understand both intuitively and quantitatively that given a specific situation, tuning which component may lead to greater marginal gain, as we have illustrated in the text right before Sec. VA. The existence of the oversaturated phase and the undersaturated phase could be prevalent in many systems similar to the ones we are focusing here. We propose that a thorough analysis of the phase boundary dependence on various factors in the system, with a combination of numerical and analytical tools, should be the first step in optimising dynamical spatiotemporal resource allocation.
Acknowledgment
This research was partially supported by Singapore ASTAR SERC “Complex Systems” Research Programme grant 1224504056.
References
 [1] Gerardo Berbeglia, JeanFran ois Cordeau and Gilbert Laporte, “Dynamic pickup and delivery problems”, European Journal of Operational Research 202 (2010) 8?15.
 [2] Francesco Ferrucci and Stefan Bock, “Realtime control of express pickup and delivery processes in a dynamic environment”, Transportation Research Part B: Methodological, 63, 114 (2014).
 [3] Bo Yang, Shen Ren, Erika Fille Legara, Zengxiang Li, Edward Y.X. Ong, Louis Lin and Christopher Monterola, “Phase Transition in Taxi Dynamics and Impact of Ridesharing”, arXiv: 1801.00462.
 [4] Ilgaz Sungur, Yingtao Ren, Fernando Ordonez, Maged Dessouky and Hongsheng Zhong, “A Model and Algorithm for the Courier Delivery Problem with Uncertainty”, Transportation Science, 44, 193?205 (2010).
 [5] Oguz Solyali, JeanFrancois Cordeau and Gilbert Laporte, “Robust Inventory Routing Under Demand Uncertainty”, Transportation Science, 46, 327340 (2011).
 [6] Arthur Flajolet, Sebastien Blandin and Patrick Jaillet, “Robust Adaptive Routing Under Uncertainty”, Operations Research, 66, 210229 (2017).
 [7] John Gunnar Carlsson and Erick Delage, “Robust Partitioning for Stochastic Multivehicle Routing”, Operations Research, 61, 727744 (2013).
 [8] Burak Eksioglu, Arif Volkan Vural and Arnold Reisman, “The vehicle routing problem: A taxonomic review”, Computers and Industrial Engineering, 57, 1472 (2009).
 [9] Furuhata, F., Dessouky, M., Ordóez, F., Brunet, M.E., Wang, X., Koening, S., “Ridesharing: the stateoftheart and future directions”, Transp. Res. Part B 57 28–46 (2013).
 [10] R.C.P. Wong, W.Y. Szeto and S.C. Wong, “Bilevel decisions of vacant taxi drivers traveling towards taxi stands in customersearch: Modeling methodology and policy implications”, Transport Policy, 33, 73 (2014).
 [11] L. MoreiraMatias, J. Gama, M. Ferreira and L. Damas, “A predictive model for the passenger demand on a taxi network”, International IEEE Conference on Intelligent Transportation Systems (ITSC), DOI: 10.1109/ITSC.2012.6338680.
 [12] R.C.P. Wong, W.Y. Szeto, S.C. Wong and H. Yang, “Modelling multiperiod customersearching behaviour of taxi drivers”, Transportmetrica B, 2, 40 (2014).
 [13] R. Bai, J. Li, J.A.D. Atkin and G. Kendall, “A novel approach to independent taxi scheduling problem based on stable matching”, Journal of the Operational Research Society, 65, 1501 (2014).
 [14] M. Nourinejad and M. Ramezani, “Developing a largescale taxi dispatching system for urban networks”, IEEE 19th International Conference on Intelligent Transportation Systems (ITSC), 2016, DOI: 10.1109/ITSC.2016.7795592

[15]
J. Ke, H. Zheng, H. Yang and X. Chen, “Shortterm forecasting of passenger demand under ondemand ride services: A spatiotemporal deep learning approach”, Transportation Research Part C,
85, 591 (2017).  [16] J. Long, W.Y. Szeto, J. Du and R.C.P. Wong, “A dynamic taxi traffic assignment model: A twolevel continuum transportation system approach”, Transportation Research Part B, 100, 222 (2017).
 [17] Hang Xu, ZhiLong Chen, Srinivas Rajagopal and Sundar Arunapuram, “Solving a Practical Pickup and Delivery Problem”, Transportation Science, 37, 347364 (2003).
 [18] Gabor Nagy, Niaz A. Wassan, M. Grazia Speranza and Claudia Archetti, “The Vehicle Routing Problem with Divisible Deliveries and Pickups”, Transportation Science, 49, 271294 (2013).
 [19] C. Archetti, R. Mansini and M. G. Speranza, “Complexity and Reducibility of the Skip Delivery Problem”, Transportation Science, 39, 182187 (2005).
 [20] Bo Yang et.al. work in progress.
 [21] Norman T.J. Bailey, “On Queuing Processes with Bulk Services”, Journal of the Royal Statistical Society. Series B (Methodological), 16, 8087 (1954).
 [22] Daniel A. Spielman. Spectral Graph theory. Combinatorial Scientific Computing. Chapman and Hall/CRC Press (2011).
 [23] Bo Yang and Qianxiao Li, “Turnbyturn Intelligent Manoeuvring of Driverless Taxis: A Recursive Value Model Enhanced by Reinforcement Learning”, IEEE Intelligent Vehicles Symposium 2018, in press.
 [24] Daniel T Gillespie, “Exact stochastic simulation of coupled chemical reactions”, The journal of physical chemistry, 81(25) 2340–2361 (1977).
Comments
There are no comments yet.