# Hybrid Multi-agent Deep Reinforcement Learning for
Autonomous Mobility on Demand Systems

Tobias Enders<EMAIL_ADDRESS>
Technical University of Munich  Germany  James Harrison
<EMAIL_ADDRESS>
Google Research  Brain Team  California  USA  Marco Pavone<EMAIL_ADDRESS>
Stanford University  California  USA  Maximilian Schiffer<EMAIL_ADDRESS>
Technical University of Munich  Germany

###### Abstract

We consider the sequential decision-making problem of making proactive request
assignment and rejection decisions for a profit-maximizing operator of an
autonomous mobility on demand system. We formalize this problem as a Markov
decision process and propose a novel combination of multi-agent Soft Actor-
Critic and weighted bipartite matching to obtain an anticipative control
policy. Thereby, we factorize the operator’s otherwise intractable action
space, but still obtain a globally coordinated decision. Experiments based on
real-world taxi data show that our method outperforms state of the art
benchmarks with respect to performance, stability, and computational
tractability.

###### keywords:

hybrid learning and optimization, multi-agent learning, deep reinforcement
learning, autonomous mobility on demand

##  1 Introduction

Mobility on demand (MoD) systems, in which a fleet of free-floating vehicles
serves customers’ ad hoc requests for point-to-point transportation, have
transformed urban mobility in recent years. Companies like Uber, Lyft, and
DiDi, made MoD more accessible compared to taxi-based ride hailing services.
Autonomous vehicles will further transform MoD systems; besides much lower
prices, a major benefit of autonomous MoD (AMoD) is its improved potential for
advanced control strategies, as a central operator obtains full control over
the entire fleet. This transformation changes the fleet operator’s control
problem substantially: MoD operators focus primarily on revenue maximization,
as human drivers’ income is (almost) a fixed cost that dominates mileage-
dependent operational cost. Contrarily, AMoD operators focus on the
maximization of their operating profit, because operational costs dominate
their total cost balance. In this context, the central operator can leverage
its full knowledge about the system state and fleet control to make improved
(proactive) dispatching decisions, i.e., request to vehicle assignment and
rejection, to maximize its profit.

Since an operator does not have profound knowledge about future trip requests,
it faces an online decision-making problem in a stochastic environment. Hence,
it is promising to apply deep reinforcement learning (DRL) to this problem.
However, AMoD systems entail many vehicles and trip requests, such that an
operator’s action space is very large and possibly time-varying as the number
of requests changes over time. It is thus infeasible to apply off-the-shelf
single-agent DRL. To solve this problem, we propose a novel combination of a
multi-agent DRL algorithm with optimization-based centralized decision-making
through weighted bipartite matching. This hybrid algorithm combines the
advantages of multi-agent approaches, DRL, and combinatorial optimization.

###  1.1 Related Work

To keep this literature overview concise, we focus on literature for
controlling (autonomous) MoD systems in the following. For a review of multi-
agent DRL, we refer to [Gronauer and Diepold (2022)] and further elaborate on
how we build on the multi-agent DRL literature in Section 3.

Classical approaches for dispatching and explicit rebalancing decisions
focused on greedy or hand-crafted feature-based policies (Liao, 2003; Zhang et
al., 2017), queueing theoretical approaches (Zhang and Pavone, 2016), and
model predictive control (MPC) (Alonso-Mora et al., 2017).

Recently, many works applied DRL in the context of (autonomous) MoD, often
including or purely focusing on explicit rebalancing (e.g., Jiao et al., 2021;
Gammelli et al., 2021; Skordilis et al., 2022; Liang et al., 2022).
Contrarily, other works focused on DRL for non-myopic dispatching, which
entails an implicit rebalancing decision that avoids additional costs due to
empty driving. Early approaches (cf. Xu et al., 2018; Wang et al., 2018) were
shown to be inferior to at least one of the subsequent works: Li et al. (2019)
proposed a mean field multi-agent actor-critic algorithm. Tang et al. (2019)
used bipartite matching based on learned $V$-values. Zhou et al. (2019)
combined a multi-agent Deep Q-Network with minimization of the Kullback-
Leibler divergence (KL-divergence) between the vehicle and the request
distribution. Finally, Sadeghi Eshkevari et al. (2022) described how DiDi
recently rolled out DRL for dispatching in practice.

Since these works strive to improve the operations of today’s MoD systems,
they aim at maximizing the drivers’ revenue or the number of orders served,
rather than at maximizing the profit of an AMoD system. While Xu et al.
(2018); Wang et al. (2018); Tang et al. (2019); Sadeghi Eshkevari et al.
(2022); Liang et al. (2022) also use a combination of multi-agent DRL and
weighted matching, they all employ value-based algorithms. Contrarily, we use
an actor-critic algorithm, enabling more advanced strategies to mitigate
problems arising from multi-agent learning, in particular, decentralized
actors with centralized critics, see Section 3. By using Soft Actor-Critic
(SAC) (Haarnoja et al., 2018), we can nevertheless benefit from the improved
sample-efficiency of off-policy algorithms.

###  1.2 Contributions

To the best of our knowledge, we are the first to consider the problem of
making proactive dispatching decisions for a profit-maximizing AMoD system
operator with DRL. We propose a novel method that combines multi-agent SAC
with centralized final decision-making through weighted matching. We perform
experiments based on real-world data and, similar to related works, benchmark
our method against a greedy policy. In addition, we are the first to compare
our method against an MPC approach. We show that our method outperforms the
greedy policy on all instances by up to 5%. Moreover, we outperform the MPC
approach in most cases. Our DRL method shows a significantly more stable
performance across varying instances, while MPC may perform arbitrarily bad—in
single cases up to 60% worse than the greedy policy. Our code can be found at
https://github.com/tumBAIS/HybridMADRL-AMoD.

##  2 Problem Formulation: Markov Decision Process

We consider a profit-maximizing operator who centrally controls a fixed-size
fleet of vehicles to serve customer trip requests revealed over time within an
operating area. The operator can accept or reject requests and dispatches
accepted requests to vehicles. These decisions must be made in real-time and
immediately, i.e., the operator cannot defer requests to a later time step, as
customers are not willing to wait for feedback. If the operator accepts a
request, customers must be picked up within a known maximum waiting time
$\omega^{\text{max}}\in\mathbb{N}_{0}$ after the request was placed. We
formalize this control problem as a Markov decision process (MDP) as follows.

Preliminaries. We consider a discrete time horizon
$\mathcal{T}=\left\\{0,1,...,T\right\\}$. During one time step, multiple
requests can enter the system. The operator makes one decision per time step
for multiple requests simultaneously, which allows to optimize over a batch of
requests. We represent the operating area as a graph $G=(V,E)$ with weight
vectors
$\prescript{e}{}{\bm{w}}=\left(\prescript{e}{}{w}^{1},\prescript{e}{}{w}^{2}\right)\in\mathbb{R}_{>0}\times\mathbb{N}$,
denoting the distance ($\prescript{e}{}{w}^{1}$) of and the time steps
($\prescript{e}{}{w}^{2}$) to traverse an edge $e\in E$. The nodes of $G$ may
represent, e.g., the centers of zones into which the operating area is
divided.

States. We describe the system state at time $t\in\mathcal{T}$ by
$\bm{S}_{t}=\bigl{(}t,\left(\prescript{t}{}{\bm{r}}^{i}\right)_{i\in\left\\{1,...,R_{t}\right\\}}\\!,\bigl{(}\bm{k}_{t}^{j}\bigr{)}_{j\in\left\\{1,...,K\right\\}}\bigr{)}$,
with $R_{t}$ being the variable number of new requests
$\prescript{t}{}{\bm{r}}^{i}$, $i\in\\{1,...,R_{t}\\}$, at time step $t$, and
$K$ vehicles $\bm{k}_{t}^{j}$, $j\in\\{1,...,K\\}$. A request
$\bm{r}=\left(\omega,o,d\right)$ consists of a waiting time
$\omega\in\mathbb{N}_{0}\cup\emptyset$, an origin $o\in V$, and a destination
$d\in V\setminus\\{o\\}$; $\omega$ tracks the elapsed time from request
placement to pickup, where we set $\omega\leftarrow\emptyset$ at pickup. We
denote a vehicle by $\bm{k}=(v,\tau,\bm{r}^{1},\bm{r}^{2})$, with position
$v\in V$ and the number of time steps $\tau\in\mathbb{N}_{0}$ left to reach
this position. Here, $v$ can either be the current node if the vehicle idles
or the next node that will be reached if the vehicle travels. Furthermore,
slightly abusing notation, a vehicle can have at most two assigned requests
$\bm{r}^{1},\bm{r}^{2}$. Assigning more requests to one vehicle is
unreasonable for realistic trip lengths and maximum waiting times. We denote
the position of vehicle $\bm{k}_{t}^{j}$ by $\prescript{j\,}{}{\\!v}_{t}$ and
denote other components of the vehicle vector likewise.

Actions. The action space describing feasible decisions of the operator is

$\displaystyle\mathcal{A}\left(\bm{S}_{t}\right)=\biggl{\\{}\left(a_{t}^{1},...,a_{t}^{R_{t}}\right)\biggl{\lvert}\
$ $\displaystyle a_{t}^{i}=0\ \lor\
\left(a_{t}^{i}=j\in\left\\{1,...,K\right\\}\ \wedge\
\prescript{j}{}{\bm{r}}_{t}^{2}=\emptyset\right)\ \forall\
i\in\left\\{1,...,R_{t}\right\\},$
$\displaystyle\sum_{i=1}^{R_{t}}\mathds{1}\\!\left(a_{t}^{i}=j\right)\leq 1\ \
\forall j\in\left\\{1,...,K\right\\}\biggr{\\}}.$ (1)

The operator can take one decision $a_{t}^{i}$ per request
$\prescript{t}{}{\bm{r}}^{i}$, $i\in\left\\{1,...,R_{t}\right\\}$, either
rejecting it ($a_{t}^{i}=0$), which means that the request leaves the system,
or assigning it to vehicle $\bm{k}^{j}$ ($a_{t}^{i}=j$), which is only
possible if the vehicle does not already have two assigned requests, i.e., if
$\prescript{j}{}{\bm{r}}_{t}^{2}=\emptyset$ holds. The final condition in (1)
implies that at most one new request is assigned to each vehicle in each time
step, which is a realistic simplification facilitating the application of a
matching algorithm. The central operator’s action space size is of order
$(K+1)^{R_{t}}$.

Transitions. We first describe the action-dependent transition from the pre-
decision to post-decision state. Then, we describe the transition from the
post-decision state to the next pre-decision state, which is independent of
the action and only determined by the system dynamics.

A reject decision has no impact on the state. When $a^{i}=j$, we add the
request to the vehicle state, i.e., if
$\prescript{j}{}{\bm{r}}^{1}=\emptyset$, then
$\prescript{j}{}{\bm{r}}^{1}\leftarrow\prescript{t}{}{\bm{r}}^{i}$, and
$\prescript{j}{}{\bm{r}}^{2}\leftarrow\prescript{t}{}{\bm{r}}^{i}$ otherwise.

The following transitions apply to all vehicles: if the vehicle picks up a
customer, i.e., if $\bm{r}^{1}\neq\emptyset\ \wedge\ \tau=0\ \wedge\
v=o\\!\left(\bm{r}^{1}\right)$, where $o\\!\left(\bm{r}^{1}\right)$ denotes
the origin of request $\bm{r}^{1}$, then
$\omega\\!\left(\bm{r}^{1}\right)\leftarrow\emptyset$. If the vehicle moves
between two nodes, i.e., if $\tau>0$, then $\tau\leftarrow\tau-1$. If the
vehicle is at a node but moves to serve a request, i.e., if $\tau=0\ \wedge\
\bm{r}^{1}\neq\emptyset$, then $v$ is replaced by the next node $v^{\prime}$
on the vehicle’s route to serve the request, going to origin or from origin to
destination, and $\tau\leftarrow\prescript{(v,v^{\prime})}{}{w}^{2}-1$. If a
vehicle drops off a customer before the next decision is made, i.e., if
$\bm{r}^{1}\neq\emptyset\ \wedge\ \omega\\!\left(\bm{r}^{1}\right)=\emptyset\
\wedge\ \tau=0\ \wedge\ v=d\\!\left(\bm{r}^{1}\right)$, we shift requests:
$\bm{r}^{1}\leftarrow\bm{r}^{2}$ and $\bm{r}^{2}\leftarrow\emptyset$. We
increment the waiting times $\omega\neq\emptyset$ of requests that have not
been picked up yet, i.e., $\omega\leftarrow\omega+1$, where $\omega$ refers to
$\omega\\!\left(\bm{r}^{1}\right)$ and/or $\omega\\!\left(\bm{r}^{2}\right)$.
Moreover, independent of the vehicles’ states, customers place new requests,
i.e.,
$\left(\prescript{t}{}{\bm{r}}^{i}\right)_{i\in\left\\{1,...,R_{t}\right\\}}$
is replaced by
$\left(\prescript{t+1}{}{\bm{r}}^{i}\right)_{i\in\left\\{1,...,R_{t+1}\right\\}}$.
Note that we do not know the underlying time-dependent probability
distribution which generates new requests, but we can simulate the resulting
requests by replaying historic data. We assume that the new requests arrive
independently of the state-action history, such that the Markov property
holds. Finally, $t\leftarrow t+1$.

Rewards. Since the operator maximizes its profit and fixed costs are
independent of the control problem, our reward function focuses on the
operating profit, which is the revenue from serving requests minus operational
costs, e.g., for fuel and maintenance. The operator obtains the revenue for a
request $\bm{r}$ when a vehicle picks up the request within the maximum
waiting time. The revenue is given by a function
$\mathrm{rev}(\bm{r})\in\mathbb{R}_{>0}$, representing the operator’s pricing
model. For improved readability, we express the profit components as functions
of the post-decision state $\bm{S}_{t^{+}}$ and write $t$ for $t^{+}$. Then,
the total revenue at time $t$ is

$\mathrm{Rev}(\bm{S}_{t})=\sum_{j=1}^{K}\,\mathds{1}\Bigl{(}\prescript{j}{}{\bm{r}}_{t}^{1}\neq\emptyset\
\wedge\ \prescript{j}{}{\tau}_{t}=0\ \wedge\
\prescript{j}{}{v}_{t}=o\\!\left(\prescript{j}{}{\bm{r}}_{t}^{1}\right)\
\wedge\
\omega\\!\left(\prescript{j}{}{\bm{r}}_{t}^{1}\right)\leq\omega^{\text{max}}\Bigr{)}\cdot\mathrm{rev}\\!\left(\prescript{j}{}{\bm{r}}_{t}^{1}\right).$

When a vehicle starts to move from $v$ to $v^{\prime}$, the operator incurs
operational costs $c\in\mathbb{R}_{>0}$ per distance unit, as commonly assumed
(see, e.g., Bösch et al., 2018). Thus, the total cost at time $t$ is

$\mathrm{Cost}(\bm{S}_{t})=c\cdot\sum_{j=1}^{K}\,\mathds{1}\Bigl{(}\prescript{j}{}{\tau}_{t}=0\
\wedge\
\prescript{j}{}{\bm{r}}_{t}^{1}\neq\emptyset\Bigr{)}\cdot\prescript{\left(\prescript{j}{}{v}_{t},\
\prescript{j}{}{v}_{t}^{\prime}\right)}{}{w}^{1}.$

The total profit at time $t^{+}$ is
$\mathrm{Profit}(\bm{S}_{t^{+}})=\mathrm{Rev}(\bm{S}_{t^{+}})-\mathrm{Cost}(\bm{S}_{t^{+}})$.
Note that $\bm{S}_{t^{+}}$ is a function of $\bm{S}_{t}$ (pre-decision) and
$\bm{a}_{t}\in\mathcal{A}\left(\bm{S}_{t}\right)$, such that we write
$\mathrm{Profit}(\bm{S}_{t^{+}})=\mathrm{Profit}(\bm{S}_{t},\bm{a}_{t})$.

The AMoD operator wants to find a policy
$\pi\left(\bm{a}_{t}\lvert\bm{S}_{t}\right)$ that maximizes the expected total
reward over all time steps, given the initial state $\bm{S}_{0}$:

$\mathrm{Profit}^{\ast}\\!\left(\bm{S}_{0}\right)=\max_{\pi}\,\mathbb{E}_{\left(\bm{S}_{t},\bm{a}_{t}\right)\sim\pi}\left[\,\left.\sum_{t=0}^{T-1}\mathrm{Profit}\left(\bm{S}_{t},\bm{a}_{t}\right)\,\right\rvert\,\bm{S}_{0}\,\right].$

To do so, we propose a hybrid DRL algorithm in the following section.

##  3 Method: Multi-agent Soft Actor-Critic with Global Matching

Analyzing our problem setting, we identify two key requirements to develop an
algorithm that constructs an effective control policy: first, it should
leverage information patterns that can be observed from historic trip data to
make non-myopic decisions. Second, it should be scalable to a realistic system
size to coordinate a large number of vehicles and requests. To account for the
second requirement, we formalize the centralized dispatching of vehicles to
requests as a bipartite matching problem (BMP). This BMP should be weighted to
allow for non-myopic dispatching decisions, anticipating the downstream impact
of decisions in a stochastic environment. The choice of weights heavily
impacts the policy’s performance. Therefore, we use DRL to parameterize these
weights, as it accounts for the downstream impact of decisions in stochastic
environments by design and allows to extract and use information from historic
data—thus, covering the first requirement.

However, single-agent DRL is not suitable for our problem setting, as the
central operator’s action space scales exponentially with the number of
vehicles and requests per time step and becomes intractable very quickly.
Thus, we leverage multi-agent DRL to factorize the action space at the price
of increased complexity, caused by having to coordinate the actions of
multiple DRL agents to finally take a centralized decision. Our hybrid
algorithm combines the advantages of multi-agent DRL with those of
combinatorial optimization: we use DRL agents as estimators to compute non-
myopic weights, serving as the input to a weighted bipartite matching
algorithm, which then makes a globally optimal and coordinated decision.

###  3.1 Overview

Figure 1 provides an overview of our method in which we leverage a DRL
algorithm to parameterize a weighted bipartite matching to take anticipatory
global dispatching decisions. To obtain a weight for each request-vehicle
combination, we consider each combination as one agent. We represent these
agents by an actor network, which we train using the SAC algorithm (Haarnoja
et al., 2018). To obtain the weights for our BMP, we post-process the actors’
outputs, such that from the perspective of the DRL agents and the computation
of policy parameter gradients, post-processing and matching are part of the
environment.

Figure 1: Overview of our method (gray text refers to parts which we use only
during training).

Since the DRL agents are a means to factorize the action space of the central
operator, rather than “real” individual agents, they can observe the global
system state. Although they should take cooperative decisions that eventually
benefit the central operator’s profit, the agents observe their own (egoistic)
rewards, not the global system reward, to avoid a credit assignment problem
(e.g., Agogino and Tumer, 2004) that otherwise occurs, in particular with many
agents. Here, we enforce coordination of the agents through the BMP and note
that varying the credit assignment scheme remains an interesting question for
future work.

Individual rewards imply a need for per-agent critic values. From one agent’s
perspective, the other agents are part of the environment. Since we train all
agents concurrently, their policies change simultaneously and the perceived
environment is non-stationary. To mitigate this, the critic gets the other
agents’ actions as an additional input (“centralized critic”), such that the
policy evaluation can explicitly account for other agents’ behavior (e.g.,
Lowe et al., 2017; Iqbal and Sha, 2019).

All agents represent a request-vehicle combination and are thus homogeneous.
Accordingly, they can share parameters and we need only one actor and one
critic network for all agents (cf. Iqbal and Sha, 2019). We can train those
centrally, i.e., the total parameter update is given by the sum of per-agent
updates. Still, the forward pass of the actor network is independent across
agents, allowing for decentralized and parallelized execution, which is
important for scalability.

Our method can handle a variable number of requests and thus a variable global
action space size, as we can use neural networks with parameter sharing for
any number of agents in parallel and the BMP does not require a fixed number
of requests. The same holds true for the vehicles, such that the system size
when testing may differ from training (see appendix).

We provide details on the individual components of our method in the
following.

###  3.2 Per-agent Post-processing and Global Matching

The actor network parameterizes a categorical probability distribution over
the two actions that can be taken for a request-vehicle combination: reject or
accept. We post-process the actor output per agent to transform it into a per-
agent score, that we then use in the global weighted matching.

Algorithm 3.2 defines the post-processing. First, we mask infeasible actions
by setting the accept probability $p_{\text{a}}$ to zero if the vehicle
already has two assigned requests. Then, we sample a reject/assign decision
from the masked probability distribution; when testing, we instead take the
argmax of the probabilities. A reject decision ($\delta=0$) at the per-agent
level implies a request-to-vehicle reject decision at the global level, such
that we set the respective score to zero. For an accept decision, we use the
accept probability as score111At first sight, it might seem more intuitive to
let the actor parameterize a continuous distribution, from which we can take a
sample to directly obtain the score, and/or have separate outputs for the
reject/accept probability and the score. We tested both approaches, but
empirically observed that they perform worse than the variant described here.
If we choose one of those approaches, we cannot compute the terms in the loss
functions (see Section 3.3) which are an expectation w.r.t. the policy, i.e.,
$\pi_{\phi}(a\,\lvert\,s,i)^{T}\cdot(...)$, since the action space is not
(purely) categorical anymore. Then, we need to sample to estimate the
expectation, as in the version of SAC for continuous action spaces, which
increases the variance. We hypothesize that this harms the algorithm’s
performance and explains our empirical observation., such that a higher
probability leads to a higher score for the weighted matching. An accept
decision at the per-agent level does not always imply an accept decision at
the global level, as the matching might assign the request to a different
vehicle.

Per-agent post-processing $p_{\text{r}},p_{\text{a}}\in[0,1]\text{ s. t.
}p_{\text{r}}+p_{\text{a}}=1\ \ ;\ \ \bm{k}^{j}$ score $s$ (*[f]reject if
already two assigned requests)$\prescript{j}{}{\bm{r}}^{2}\neq\emptyset$
$p_{\text{r}}\leftarrow 1\ ,\ p_{\text{a}}\leftarrow 0$ training
$\delta\sim\mathrm{Categorical}\bigl{(}(p_{\text{r}},p_{\text{a}})\bigr{)}$
*sample $\delta\in\\{0,1\\}$ when training $\delta\leftarrow 0\text{ if
}p_{\text{r}}\geq p_{\text{a}},\ \delta\leftarrow 1\text{ if
}p_{\text{r}}<p_{\text{a}}$ *[f]argmax when testing $\delta=0$ $s\leftarrow 0$
*score is zero if rejected $s\leftarrow p_{\text{a}}$ *[f]score is accept
probability if accepted

We use all agents’ scores to create a bipartite graph, with vehicles and
requests as nodes, and edges between all vehicle and request nodes for which
we obtain per-agent accept decisions (i.e., $s>0$). The edges’ weights
correspond to the respective scores. We solve the resulting maximum weighted
BMP (formally defined in the appendix) using the Hungarian algorithm (Kuhn,
1955) to get a globally coordinated decision, where each request is assigned
at most once.

###  3.3 Multi-agent Soft Actor-Critic

SAC is an entropy-regularized, off-policy actor-critic algorithm. It trains a
stochastic policy $\pi\left(\bm{a}_{t}\lvert\bm{S}_{t}\right)$ with entropy
maximization, incentivizing exploration through a random policy:

$\pi^{\ast}=\arg\max_{\pi}\,\mathbb{E}_{\left(\bm{S}_{t},\bm{a}_{t}\right)\sim\pi}\left[\,\sum_{t=0}^{T-1}\mathrm{Profit}\left(\bm{S}_{t},\bm{a}_{t}\right)+\alpha\,H\\!\left(\pi\left(\,\cdot\,\lvert\,\bm{S}_{t}\right)\right)\right].$

The entropy of the policy is defined as
$H\\!\left(\pi\left(\,\cdot\,\lvert\,\bm{S}_{t}\right)\right)=-\mathbb{E}_{\bm{a}_{t}\sim\pi}\log\pi\left(\bm{a}_{t}\lvert\bm{S}_{t}\right)$
and the entropy coefficient $\alpha\in\mathbb{R}_{\geq 0}$ is a hyperparameter
that controls the exploitation/exploration trade-off. While SAC was originally
developed for continuous action spaces in Haarnoja et al. (2018), it can also
be applied to discrete actions (Christodoulou, 2019).

We chose to use an actor-critic algorithm to enable our multi-agent approach
of decentralized actors with centralized critics as explained in Section 3.1.
Since exploration is paramount for our problem setting, we use SAC, which lets
us explicitly tune how much the policy explores.

We parameterize the actor network with parameters $\phi$ and the critic with
$\theta$. SAC uses two critic networks, $Q\in\left\\{Q^{1},Q^{2}\right\\}$, as
well as corresponding target networks with parameters $\bar{\theta}$, which
are an exponential moving average of the primary parameters $\theta$. Based on
the multi-agent approach as described in Section 3.1, the loss function for
the actor with shared parameters is

$J_{\pi}(\phi)=\mathbb{E}_{s\sim
D}\left[\sum_{i}\pi_{\phi}(a\,\lvert\,s,i)^{T}\cdot\left(\alpha\log\pi_{\phi}(a\,\lvert\,s,i)-\min_{j\in\\{1,2\\}}\left\\{Q_{\theta}^{j}\left(a\,\lvert\,s,i,\bar{a}_{-i}\right)\right\\}\right)\right].$

Here, we use a simplified notation for improved readability: we denote a
transition by $\left(s,\bar{a},r,s^{\prime}\right)$, with global states
$s,s^{\prime}$, global action $\bar{a}$ (after the matching), and rewards $r$.
For agent $i$, $r_{i}$ denotes its reward, $\bar{a}_{-i}$ is the global action
except for agent $i$’s action, and $a$ is a per-agent action (reject/assign),
such that $\pi_{\phi}(a\,\lvert\,s,i)\in[0,1]^{2}$ and
$Q_{\theta}^{j}\left(a\,\lvert\,s,i,\bar{a}_{-i}\right)\in\mathbb{R}^{2}$. We
sample states (or transitions) from the replay buffer $D$ and denote the
discount factor by $\gamma$. For the actor loss, we do not sample the global
action from the replay buffer, but compute it based on the state $s$ and the
current policy, as in Iqbal and Sha (2019). For each of the two critics
$Q\in\left\\{Q^{1},Q^{2}\right\\}$, the loss function is

$\displaystyle J_{Q}(\theta)$
$\displaystyle=\mathbb{E}_{(s,\bar{a},r,s^{\prime})\sim
D}\left[\sum_{i}\frac{1}{2}\left(Q_{\theta}\left(a\,\lvert\,s,i,\bar{a}_{-i}\right)\Big{\rvert}_{\bar{a}_{i}}-y_{i}\right)^{2}\right]\text{,
with}$ $\displaystyle y_{i}$
$\displaystyle=r_{i}+\gamma\cdot\pi_{\phi}\left(a^{\prime}\,\big{\lvert}\,s^{\prime},i\right)^{T}\cdot\left(\min_{j\in\\{1,2\\}}\left\\{Q_{\bar{\theta}}^{j}\left(a^{\prime}\,\big{\lvert}\,s^{\prime},i,\bar{a}^{\prime}_{-i}\right)\right\\}-\alpha\log\pi_{\phi}\left(a^{\prime}\,\big{\lvert}\,s^{\prime},i\right)\right).$

Here, the notation $\big{\rvert}_{\bar{a}_{i}}$ means “evaluated at
$\bar{a}_{i}$”, i.e., of the two $Q$-values that we compute for the two
possible actions of agent $i$, we use the one corresponding to the global
decision $\bar{a}$. The term after $\gamma$ is the $V$-value estimate for
$s^{\prime}$ based on $\bar{\theta}$, for which we compute the next global
action $\bar{a}^{\prime}$ with the current policy. The number of requests and
thus the number of agents can change between subsequent time steps. This poses
a numerical problem for the critic loss computation, which requires the same
number of agents for $s$ and $s^{\prime}$. We solve this problem by amending
the requests in $s^{\prime}$ when saving a transition to the replay buffer and
provide details on this in the appendix.

The actor network obtains all vehicle states and requests for which a decision
must be made in the current time step as an input. To deal with these
(potentially) many inputs, we train a single request embedding and a single
vehicle embedding to encode all requests and vehicles, respectively. To
account for the variable number of requests and to let each agent focus on the
parts of the input that are important for this particular agent, we equip the
neural network with an attention mechanism (cf. Holler et al., 2019; Kullman
et al., 2022). Together with the request and vehicle embeddings for the agent
and additional features, we pass the context computed by the attention
mechanism to a sequence of feedforward layers. The critic network has the same
architecture, but receives the global action as an additional input. We remove
the action of the agent from this input, since the critic outputs $Q$-values
for both possible actions. Further details on the neural networks, e.g., a
formal description of the attention mechanism and hyperparameters, can be
found in the appendix.

##  4 Experiments

To validate our method, we perform experiments based on historic taxi data
that is publicly available for New York City (NYC TLC, 2015). We use a hexagon
grid for spatial discretization and consider two different instances: one with
11 small zones (approx. 500 meters distance between neighboring zones) and one
with 38 large zones (approx. 1 km distance), both in Manhattan. We consider
the time interval from 8:30 am to 9:30 am during morning rush hour as one
episode. Our data set contains data for 245 different dates in 2015, which we
split into 200 training dates, 25 validation dates, and 20 test dates. We use
a time step size of one minute and choose revenue and cost parameters such
that a vehicle that serves a customer without empty driving achieves an
operating profit margin of 10%. We consider different numbers of vehicles to
simulate different degrees of supply shortage; only cases with supply shortage
are interesting, as the operator can serve each request immediately in the
case of infinite supply, such that a myopic policy would be sufficient. For
additional details on the data set, system setup, and hyperparameters, we
refer to the appendix.

To benchmark our method, we compare its test performance against two
“classical”, non RL-based, algorithms: a greedy policy and an MPC approach.
The greedy policy accepts any request that can be served with a positive
profit (accounting for the cost for empty driving to the pickup location) and
rejects all others. It is a reasonable choice when there is no reliable
estimate of future requests. If we have such an estimate, it is promising to
apply MPC (see, e.g., Alonso-Mora et al., 2017). We adapt this approach to our
setting, using a request distribution estimate for mixed-integer-based
receding-horizon optimization. Details on both benchmarks can be found in the
appendix.

##  5 Results and Discussion

We provide plots illustrating the training process in the appendix. Figure 2
summarizes the performance of greedy, MPC, and our RL method on the test data
for all considered instances. On average, our RL method always outperforms the
greedy policy, by up to 5% over the 20 test dates. For individual dates, RL
outperforms greedy by up to 17%. It performs by at most 6% worse than greedy
for less than 20% of the individual dates. MPC is (substantially) worse than
greedy and RL in many cases, although it sometimes outperforms the RL method.
This means that MPC can provide a benefit in certain situations, but comes
with an unstable performance across instances, which limits its practical
applicability. In particular, MPC does not perform well in situations where
there is a large shortage of vehicles, which are handled well by our RL
method. Thus, our method provides a stable alternative, that always achieves
at least the greedy performance and outperforms it by a substantial margin in
many cases. Note that the order of magnitude of this performance improvement
is significant for our application area (cf. Sadeghi Eshkevari et al., 2022).
Given the large scale at which AMoD systems operate, the seemingly small
percentage improvements translate into significant monetary value for the
operator.

Figure 2: Test performance of MPC and our RL method compared to greedy (greedy
is 0%, values $<0\%$ indicate a performance worse than greedy). Each dot
represents one test date.

Figure 3 shows the performance of MPC and RL for the instance with 38 large
zones and different amounts of training data, i.e, different estimation
qualities of the request probability distribution. Since our problem setting
excludes fixed costs, additional resources, i.e., vehicles, are free of
charge. Thus, the greedy policy performs better with very few or many
vehicles, compared to instances with a medium number of vehicles. With few
vehicles, many requests are available for each vehicle, such that vehicles are
rarely idle or drive without a customer. With many vehicles, most requests can
be served quickly without empty driving because vehicles are usually
available. Consequently, with sufficient training data, for both MPC and our
RL approach, the performance gain vs. greedy is largest for a medium number of
vehicles. However, with few vehicles, MPC performs worse than greedy, as it is
not robust against mistakes when sampling future requests. Such errors have a
larger effect with fewer vehicles. With less training data, the performance
gain of our method decreases by about one percentage point, but it remains
reliably better than greedy. For all instances except the 250 vehicles case,
the performance loss is much larger for MPC; it is not always able to sustain
a performance better than greedy, even for instances where it outperforms
greedy with more training data. With 250 vehicles, there are many resources
free of charge, such that the mistakes made by MPC have such a small effect
that it is robust against a poor estimation quality. Based on these
observations, we conclude that our RL method is more robust against a poor
estimation quality due to insufficient training data than MPC. These results
might seem surprising, as RL is in general not very sample-efficient—although
SAC has better sample-efficiency than most policy gradient-based algorithms,
since it is an off-policy algorithm and uses a replay buffer. However, for our
problem setting, less training data does not mean that the RL agents must
learn from fewer samples, as the available training data can be replayed
multiple times in the simulated environment. The performance loss that we
observe for the RL method is more likely due to the decreasing diversity of
the training data to which the RL agents are exposed, leading to less
generalization.

Figure 3: Mean test performance of MPC and RL vs. greedy as a function of the
KL-divergence between the true request distribution and the distribution
estimated from the training data (with different amounts of training data,
resulting in different KL-divergence values). We use the request distribution
estimated from the real data for 38 large zones as the true distribution and
run the experiments with synthetic data sampled from this distribution.

Finally, a major advantage of our RL method over MPC is its shorter
computational time during execution. We can train the network parameters
offline in advance and easily scale the online execution, because the per-
agent actor computations are fast and straightforward to parallelize. Figure 4
shows that although we solve a combinatorial optimization problem in each time
step, the computational time of our method is very short even for large system
sizes. On the other hand, the computational time of MPC increases quickly for
large instances: for 3000 vehicles, the action computation for a single step
takes more than 30 seconds with MPC, while our RL method (including the
matching) takes less than 0.2 seconds. Thus, the practical application of MPC
at scale is greatly limited by its computational time, while our RL method can
be scaled to much larger system sizes.

Figure 4: Time to compute one action (mean over 60 steps, based on 38 large
zones, number of requests scaled with number of vehicles). MPC solve refers to
solving the mixed integer program (MIP), MPC total also includes the
generation of the MIP instance in each step.

##  6 Conclusion

We consider the dispatching problem of a profit-maximizing AMoD operator with
centralized control over a fleet of autonomous vehicles, who accepts (and
serves) or proactively rejects requests in real-time. To solve this problem,
we use a combination of multi-agent SAC with centralized final decision-making
through weighted matching. Our experiments based on real-world data show that
our method outperforms two strong benchmarks on most problem instances, that
it is stable across these instances and robust against a poor estimation of
the request distribution, and that it can be easily scaled to large system
sizes. In future work, we will investigate the use of global instead of ego
rewards. Furthermore, we will extend our framework to more complex use cases,
e.g., dispatching and charge scheduling.

##  Appendix A Method

In the following, we provide complementary details on our method. First, we
describe how we assign per-agent rewards. Second, we formally define the
weighted BMP. Third, we explain our solution to the dimensionality problem in
the critic loss computation. Fourth, we provide additional details on the
neural networks. Fifth, we give an overview of alternative approaches to post-
processing and weighted matching which we tested and discarded.

###  A.1 Per-agent Rewards

Agents observe their own (egoistic) rewards. When matching agents, we can
already compute the (potentially negative) profit that will result from the
decision to assign a certain request to a certain vehicle. This profit
consists of the revenue that will be obtained for the request, minus the
operational costs to drive from the vehicle’s position after it finished
serving its current request (if any) to the new request’s origin and the
operational costs to drive from there to the new request’s destination. If a
request is matched to a vehicle at the global level, the corresponding agent
immediately observes this profit as the reward. All agents for which the
request is not matched to the vehicle at the global level observe a reward of
zero. The sum over all these per-agent rewards equals the system reward,
although we virtually forward rewards in time. When we sample transitions from
the replay buffer, we normalize the sampled rewards by dividing them by the
standard deviation of all rewards currently stored in the replay buffer (cf.
Kurin et al., 2022).

###  A.2 Bipartite Matching Problem

At time step $t$, the weighted BMP is formally defined as

$\displaystyle\max_{x_{ij}}\ $
$\displaystyle\sum_{i=1}^{R_{t}}\sum_{j=1}^{K}s_{ij}\cdot x_{ij}$ s. t.
$\displaystyle\sum_{i=1}^{R_{t}}x_{ij}\leq 1\quad\text{for all}\ \
j\in\left\\{1,...,K\right\\},$ $\displaystyle\sum_{j=1}^{K}x_{ij}\leq
1\quad\text{for all}\ \ i\in\left\\{1,...,R_{t}\right\\},$

where $s_{ij}$ is the score computed for request $\prescript{t}{}{\bm{r}}^{i}$
and vehicle $\bm{k}^{j}$, and the decision variables are $x_{ij}=1$ if
$\prescript{t}{}{\bm{r}}^{i}$ is assigned to $\bm{k}^{j}$ and $x_{ij}=0$
otherwise.

###  A.3 Critic Loss Computation Despite Variable Number of Agents

The number of agents changes between time steps, since the number of requests
varies over time. However, for the TD-update of the $Q$-function, we need the
same number of agents at time $t$ and $t+1$. Thus, when adding transitions to
the replay buffer, we replace the requests state for $t+1$ by the requests
state for $t$ to obtain a matching number of agents. Given the small time step
size, it is strongly stochastic which requests occur at which time step, and
the underlying distribution can be assumed to be very similar for two
consecutive time steps, such that the requests at $t$ are in expectation a
similarly realistic sample for $t+1$ as the requests which we actually
observe.

###  A.4 Neural Networks

We first describe the features which we use as inputs to the neural networks.
We identify the zones which we use for spatial discretization by a horizontal
and a vertical index. All locations (e.g., the request origin) correspond to
one of these zones. We encode them by a vector with the horizontal and the
vertical index, each normalized to $[0,1]$. We also tried to use sinusoidal
positional encodings (cf. Vaswani et al., 2017) instead, but this did not
improve the performance of our method. A request encoding then consists of:

  * •

The encoding of the request origin

  * •

The encoding of the request destination

  * •

The distance from origin to destination on the graph $G$, normalized to
$[0,1]$

We encode a vehicle state by:

  * •

The encoding of the vehicle’s position, where we use the current position of
the vehicle if it does not have an assigned request or the destination of the
assigned request that will be served last

  * •

The number of time steps left to reach this position, normalized to $[0,1]$

  * •

The number of assigned requests, normalized to $[0,1]$

Apart from request and vehicle states, we use some additional features:

  * •

The current time step, normalized to $[0,1]$

  * •

A flag in $\\{0,1\\}$ indicating if the vehicle under consideration will be
able to serve the request under consideration within the maximum waiting time
if it is matched to the vehicle

  * •

The time steps to reach the position summed over all vehicles, normalized to
$[0,1]$, indicating how busy the fleet currently is

  * •

The number of requests placed since the current episode started, divided by
the count of requests that are placed on average until the current time step,
indicating how much demand was observed compared to an average episode

Next, we describe the actor network. The request and vehicle encodings are
used as the input for the request and vehicle embeddings, respectively. Both
the request and the vehicle embedding are a feedforward layer with 32 units
and ReLU activation. We denote the resulting embedding vectors by
$\bm{e}_{r^{i}}$, $i\in\left\\{1,...,R_{t}\right\\}$, and $\bm{e}_{k^{j}}$,
$j\in\left\\{1,...,K\right\\}$. The attention mechanism computes a global
context, which is the concatenation of a requests context and a vehicles
context. The requests context is computed as
$\bm{c}_{r}=\sum_{i=1}^{R_{t}}\beta_{r^{i}}\cdot\bm{e}_{r^{i}}$ with
$\beta_{r^{i}}=\sigma\left(\bm{w}_{r}\cdot\tanh\left(W_{r}\cdot\bm{e}_{r^{i}}\right)\right)\in\mathbb{R}$,
where $\sigma$ is the sigmoid activation function and
$\bm{w}_{r}\in\mathbb{R}^{256}$ as well as $W_{r}\in\mathbb{R}^{256\times 32}$
are trainable parameters (the parameters are the same across all requests).
The vehicles context is computed as
$\bm{c}_{k}=\sum_{i=1}^{K}\beta_{k^{i}}\cdot\bm{e}_{k^{i}}$ with
$\beta_{k^{i}}=\sigma\left(\bm{w}_{k}\cdot\tanh\left(W_{k}\cdot\bm{e}_{k^{i}}\right)\right)\in\mathbb{R}$
and trainable parameters $\bm{w}_{k}\in\mathbb{R}^{128}$,
$W_{k}\in\mathbb{R}^{128\times 32}$ (the parameters are the same across all
vehicles). Thereby, we obtain a fixed size global representation of variable
sized inputs. The structure of the embeddings and the attention mechanism is
similar to Holler et al. (2019); Kullman et al. (2022). We use the global
context as well as the request and vehicle embedding corresponding to the
agent under consideration together with the aforementioned additional inputs
as the input to a sequence of feedforward layers. We use five layers with ReLU
activation and 1024, 512, 128, 32, and 8 units, respectively, followed by the
output layer of size 2 with softmax activation. For all layers, we use L2
regularization with a regularization coefficient of $10^{-4}$. We chose the
specific architecture (number and size of layers) through hyperparameter
tuning.

The critic networks are identical to the actor network except for the
following differences: For all request and vehicle states that do not
correspond to the agent under consideration, we add the action to their
encodings. For the requests, this is a flag in $\\{0,1\\}$ indicating if the
request was rejected or accepted. For the vehicles, this is the origin and
destination of the request that was newly assigned to the vehicle (zeros if no
request was newly assigned to the vehicle). Moreover, we use the request and
vehicle states of the agent under consideration as inputs to the feedforward
layer of size 1024, but not for the embedding and attention layers. Finally,
there is no activation function in the output layer.

###  A.5 Alternatives to Post-processing and Weighted Matching

We tested two alternatives to our weighted matching approach to obtain a
global decision from the per-agent actor network outputs: generating a global
probability distribution and non-weighted matching.

For the global probability distribution, we use the per-agent reject/accept
probabilities to construct the joined probability distribution over all
feasible global decisions. We then sample, by taking the argmax when testing,
from the global distribution to obtain a global decision. While this approach
removes the need for a combinatorial optimization algorithm and gives
promising results on very small instances, it becomes intractable quickly due
to the exponentially increasing size of the global action space.

For the non-weighted matching, we post-process the per-agent actor network
outputs as described in Algorithm 3.2, but do not use scores to construct a
weighted bipartite graph. Instead, we create a non-weighted bipartite graph
based on the per-agent reject/assign decisions $\delta$ and use a maximum
(non-weighted) bipartite matching algorithm to obtain a global decision. While
this approach is scalable, it performs worse than the weighted matching
variant.

##  Appendix B Experiments

In the following, we provide complementary information on our experiments.
First, we give more details on the system setup and how we pre-process the
real-world taxi data set. Second, we state the hyperparameter values used in
our experiments. Third, we provide details on the two benchmark algorithms.

###  B.1 Data Set and System Setup

We use yellow taxi trip records from the year 2015 and exclude weekends and
holidays. We assume that requests are placed at the time reported as the
pickup time in the data set. Besides, we use the pickup and dropoff
longitude/latitude from the data set and keep only trips for which pickup and
dropoff coordinates are located on the main island of Manhattan. We discretize
space with a hexagon grid as shown in Figure 5. We assign a pickup and dropoff
zone to each request based on the shortest distance from the
longitude/latitude coordinates and remove trips that start and end in the same
zone.

Figure 5: Hexagon grid laid over Manhattan for spatial discretization. The
operating areas which we consider are marked in green.

Each zone is represented by a node in the graph $G$, that contains edges only
between nodes that represent neighboring zones. The distance between
neighboring zones is 459 meters and 917 meters for the small and large zones,
respectively, and we assume a travel time of two and five time steps based on
a realistic average driving speed. When vehicles travel between non-
neighboring zones, they take the shortest route on $G$.

We consider the two operating areas depicted in Figure 5, i.e., we consider
only requests that have a pickup and dropoff location within the green area.
For the 38 large zones, we downscale the trip data by a factor of 20, i.e., we
use only every 20th request for our simulation, to have a system size suitable
for the hardware that we used for the experiments. On average, this results in
360 requests per episode for the 11 small zones instance, with up to 23
requests in a single time step. For the 38 large zones instance, we observe
828 requests per episode on average, with up to 20 requests in a single time
step. Note that the mean trip distance is larger for the 38 large zones
instance, such that the number of vehicles required to serve a certain number
of requests is larger than for 11 small zones.

We assume a maximum waiting time of five minutes. To achieve an operating
profit margin of 10% when a request is served without empty driving to the
pickup location, we set the revenue to 5.00 USD per km and the operational
costs to 4.50 USD per km. Note that these numbers might be considered to be
unrealistic, but can be scaled to a different level without any effect on the
control problem and our results, since we report all results relative to the
greedy performance.

With the eight to 80 vehicles which we consider for the 11 small zones, the
greedy policy serves 27% to 78% of the requests. For the 38 large zones with
50 to 250 vehicles, the greedy policy serves 30% to 76% of the requests.

###  B.2 Hyperparameters

We train for 200,000 steps, update the network parameters every 20 steps, and
test the performance of the current policy on the validation data every 2,880
steps (48 episodes). During the first 20,000 steps, we collect experience with
a random policy and do not update the network parameters.

For the critic loss, we use the Huber loss with a delta value of 10 instead of
the squared error. Moreover, we use gradient clipping with a clipping ratio of
10 for actor and critic gradients. We use the Adam optimizer with a learning
rate of $3\cdot 10^{-4}$. We sample batches of size 128 from a replay buffer
with maximum size 100,000. We set the discount factor to 0.9 since this gives
a better and particularly more stable performance than other values which we
tested. For the update of the target critic parameters we use an exponential
moving average with smoothing factor $5\cdot 10^{-3}$. We tune the entropy
coefficient individually per instance and use values between 0.35 and 1.30
across the experiments reported in this paper.

We repeat each training run with multiple random seeds and use the model with
the best validation performance across runs to test the performance of our
method on the test data set. Results reported throughout this paper correspond
to these test results. For the 11 small zones instance, we use five random
seeds, while we use three random seeds for the 38 large zones instance.

The MPC results are based on an average over multiple random seeds. For the
results in Figure 2, we use five random seeds, while we use three random seeds
for the results in Figure 3.

###  B.3 Benchmarks

We benchmark our algorithm against two well-known policies that are as
follows.

Greedy. The greedy policy considers requests in their arrival order. Whenever
there is at least one vehicle that will be able to serve the request within
the maximum waiting time and with a positive profit, the greedy policy accepts
the request. If there is no such vehicle, the policy rejects the request. The
profit calculation takes into account the revenue from serving the request and
the cost to drive from the request’s origin to its destination, as well as the
cost to drive from the destination of the request that the vehicle currently
serves (the position of the vehicle if it is idle) to the origin of the new
request. If there is more than one vehicle which fulfills these conditions,
the greedy policy assigns the request to the vehicle that will be closest to
the request origin once it has finished its current job.

MPC. We adapt the approach by Alonso-Mora et al. (2017) to our problem
setting. First, for each 15 minute interval, we estimate a probability
distribution over origin-destination pairs from the training data, i.e., we
obtain the probability that a new request shall be picked up at this origin
and dropped off at this destination based on frequentist statistics. Here, we
use Laplace smoothing to mitigate a bias from potentially sparse data. In
addition, we estimate the number of requests that can be expected. We then use
those estimates for online decision-making as follows: in each time step, we
observe the new (real) requests. In addition, for some sampling horizon, we
sample the expected number of requests from the estimated probability
distribution (virtual requests). With the current vehicle states, real and
virtual requests, we solve an offline optimization problem, maximizing for
total profit assuming perfect information over the sampling horizon, with
mixed integer programming. The solution gives a decision for the real
requests, which we use as the action for the current time step. We repeat this
process in a receding horizon fashion. We set the sampling horizon to five
minutes, which turned out to yield the best performance.

##  Appendix C Complementary Results

In the following, we provide additional result plots. Figure 6 shows the
validation reward to illustrate the training process for two exemplary
instances: 11 small zones with 24 vehicles and 38 large zones with 100
vehicles. Figure 7 shows the test performance of a policy trained with some
original number of agents in a system with an increased number of agents. The
performance does not deteriorate, illustrating that our method is non-
parametric w.r.t. the number of agents.

Figure 6: Validation reward along the training process for two exemplary
instances. The green line is the mean over random seeds, the shaded area
depicts the minimum and maximum values over random seeds. Plots in the right
column are zoomed in versions of the plots in the left column to make the part
of the training process after a reward of zero is reached better visible.
Figure 7: Average performance of our method vs. greedy over the 20 test dates
for the 38 large zones instance. We train the RL agents with 50, 100, 150, and
250 vehicles on the original system. The left plot shows the test performance
on the original system. The right plot shows the test performance of the same
policies, without additional training, for an increased system size, with
twice as many vehicles and requests (downscaled by factor 10 instead of 20),
i.e., we use the policy trained with 50 vehicles for the 100 vehicles on the
larger system, and likewise for the other vehicle numbers.

We thank Gerhard Hiermann for his help with the implementation of the MIP for
the MPC-based benchmark. The work of Tobias Enders was supported by a
fellowship within the IFI program of the German Academic Exchange Service
(DAAD). This work was supported by the German Research Foundation (DFG) under
grant no. 449261765.

## References

  * Agogino and Tumer (2004) Adrian K. Agogino and Kagan Tumer.  Unifying temporal and structural credit assignment problems.  In _Proceedings of the Third International Joint Conference on Autonomous Agents and Multiagent Systems_ , 2004.
  * Alonso-Mora et al. (2017) Javier Alonso-Mora, Alex Wallar, and Daniela Rus.  Predictive routing for autonomous mobility-on-demand systems with ride-sharing.  In _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2017.
  * Bösch et al. (2018) Patrick M. Bösch, Felix Becker, Henrik Becker, and Kay W. Axhausen.  Cost-based analysis of autonomous mobility services.  _Transport Policy_ , 64:76–91, 2018.
  * Christodoulou (2019) Petros Christodoulou.  Soft actor-critic for discrete action settings.  _arXiv preprint arXiv:1910.07207_ , 2019.
  * Gammelli et al. (2021) Daniele Gammelli, Kaidi Yang, James Harrison, Filipe Rodrigues, Francisco C. Pereira, and Marco Pavone.  Graph neural network reinforcement learning for autonomous mobility-on-demand systems.  In _Proceedings of the 60th IEEE Conference on Decision and Control (CDC)_ , 2021.
  * Gronauer and Diepold (2022) Sven Gronauer and Klaus Diepold.  Multi-agent deep reinforcement learning: a survey.  _Artificial Intelligence Review_ , 55(2):895–943, 2022.
  * Haarnoja et al. (2018) Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine.  Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor.  In _Proceedings of the 35th International Conference on Machine Learning_ , 2018.
  * Holler et al. (2019) John Holler, Risto Vuorio, Zhiwei Qin, Xiaocheng Tang, Yan Jiao, Tiancheng Jin, Satinder Singh, Chenxi Wang, and Jieping Ye.  Deep reinforcement learning for multi-driver vehicle dispatching and repositioning problem.  In _2019 IEEE International Conference on Data Mining (ICDM)_ , 2019\.
  * Iqbal and Sha (2019) Shariq Iqbal and Fei Sha.  Actor-attention-critic for multi-agent reinforcement learning.  In _Proceedings of the 36th International Conference on Machine Learning_ , 2019.
  * Jiao et al. (2021) Yan Jiao, Xiaocheng Tang, Zhiwei Qin, Shuaiji Li, Fan Zhang, Hongtu Zhu, and Jieping Ye.  Real-world ride-hailing vehicle repositioning using deep reinforcement learning.  _Transportation Research Part C: Emerging Technologies_ , 130, 2021\.
  * Kuhn (1955) Harold W. Kuhn.  The Hungarian method for the assignment problem.  _Naval Research Logistics Quarterly_ , 2(1–2):83–97, 1955.
  * Kullman et al. (2022) Nicholas D. Kullman, Martin Cousineau, Justin C. Goodson, and Jorge E. Mendoza.  Dynamic ride-hailing with electric vehicles.  _Transportation Science_ , 56(3):775–794, 2022\.
  * Kurin et al. (2022) Vitaly Kurin, Alessandro De Palma, Ilya Kostrikov, Shimon Whiteson, and M. Pawan Kumar.  In defense of the unitary scalarization for deep multi-task learning.  _arXiv preprint arXiv:2201.04122_ , 2022.
  * Li et al. (2019) Minne Li, Zhiwei Qin, Yan Jiao, Yaodong Yang, Jun Wang, Chenxi Wang, Guobin Wu, and Jieping Ye.  Efficient ridesharing order dispatching with mean field multi-agent reinforcement learning.  In _The World Wide Web Conference_ , 2019.
  * Liang et al. (2022) Enming Liang, Kexin Wen, William H. K. Lam, Agachai Sumalee, and Renxin Zhong.  An integrated reinforcement learning and centralized programming approach for online taxi dispatching.  _IEEE Transactions on Neural Networks and Learning Systems_ , 33(9):4742–4756, 2022.
  * Liao (2003) Ziqi Liao.  Real-time taxi dispatching using global positioning systems.  _Communications of the ACM_ , 46(5):81–83, 2003\.
  * Lowe et al. (2017) Ryan Lowe, Yi Wu, Aviv Tamar, Jean Harb, Pieter Abbeel, and Igor Mordatch.  Multi-agent actor-critic for mixed cooperative-competitive environments.  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , 2017.
  * NYC TLC (2015) NYC TLC.  Trip record data.  Online, 2015.  URL https://www.nyc.gov/site/tlc/about/tlc-trip-record-data.page.
  * Sadeghi Eshkevari et al. (2022) Soheil Sadeghi Eshkevari, Xiaocheng Tang, Zhiwei Qin, Jinhan Mei, Cheng Zhang, Qianying Meng, and Jia Xu.  Reinforcement learning in the wild: Scalable RL dispatching algorithm deployed in ridehailing marketplace.  In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , 2022.
  * Skordilis et al. (2022) Erotokritos Skordilis, Yi Hou, Charles Tripp, Matthew Moniot, Peter Graf, and David Biagioni.  A modular and transferable reinforcement learning framework for the fleet rebalancing problem.  _IEEE Transactions on Intelligent Transportation Systems_ , 23(8):11903–11916, 2022.
  * Tang et al. (2019) Xiaocheng Tang, Zhiwei Qin, Fan Zhang, Zhaodong Wang, Zhe Xu, Yintai Ma, Hongtu Zhu, and Jieping Ye.  A deep value-network based approach for multi-driver order dispatching.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2019.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In _Advances in Neural Information Processing Systems_ , 2017.
  * Wang et al. (2018) Zhaodong Wang, Zhiwei Qin, Xiaocheng Tang, Jieping Ye, and Hongtu Zhu.  Deep reinforcement learning with knowledge transfer for online rides order dispatching.  In _2018 IEEE International Conference on Data Mining (ICDM)_ , 2018\.
  * Xu et al. (2018) Zhe Xu, Zhixin Li, Qingwen Guan, Dingshui Zhang, Qiang Li, Junxiao Nan, Chunyang Liu, Wei Bian, and Jieping Ye.  Large-scale order dispatch in on-demand ride-hailing platforms: A learning and planning approach.  In _Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2018.
  * Zhang et al. (2017) Lingyu Zhang, Tao Hu, Yue Min, Guobin Wu, Junying Zhang, Pengcheng Feng, Pinghua Gong, and Jieping Ye.  A taxi order dispatch model based on combinatorial optimization.  In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2017.
  * Zhang and Pavone (2016) Rick Zhang and Marco Pavone.  Control of robotic mobility-on-demand systems: A queueing-theoretical perspective.  _International Journal of Robotics Research_ , 35(1–3):186–203, 2016.
  * Zhou et al. (2019) Ming Zhou, Jiarui Jin, Weinan Zhang, Zhiwei Qin, Yan Jiao, Chenxi Wang, Guobin Wu, Yong Yu, and Jieping Ye.  Multi-agent reinforcement learning for order-dispatching via order-vehicle distribution matching.  In _Proceedings of the 28th ACM International Conference on Information and Knowledge Management_ , 2019.


# Enhancement of a silicon waveguide single photon source by temporal
multiplexing

Jeremy C. Adcock<EMAIL_ADDRESS>Davide Bacco  Yunhong Ding
<EMAIL_ADDRESS>
(October 2021)

##  Abstract

Efficient generation of single photons is one of the key challenges of
building photonic quantum technology, such as quantum computers and long-
distance quantum networks. Photon source multiplexing—where successful pair
generation is heralded by the detection of one of the photons, and its partner
is routed to a single mode output—has long been known to offer a concrete
solution, with output probability tending toward unity as loss is reduced.
Here, we present a temporally multiplexed integrated single photon source
based on a silicon waveguide and a low-loss fibre switch and loop
architecture, which achieves enhancement of the single photon output
probability of $4.5\pm 0.5$, while retaining $g^{(2)}(0)=0.01$.

Fig. 1: Experimental set up and temporal-loop multiplexing dynamics. a.
Schematic of our experiment. Photon pairs are probabilistically generated in
each time bin via SFWM in a silicon waveguide. Detection of one of the pair
causes the other to be switched into a fibre loop, routing the signal photon
to the output time bin and boosting the output probability. b. A pulse train
defines $m$ single temporally independent single photon sources. c. Pulse
picking scheme and switching logic for temporal-loop multiplexing

##  1\. Introduction

Photonics is poised to play a key role in future quantum technologies, both in
quantum computers and the quantum networks that will connect them. In the last
decade, progress in photonic quantum information processing has accelerated
dramatically, with historic challenges such as achieving near-unity photon
purity recently being overcome1, 2. This has enabled demonstrations of quantum
advantage with photons3, 4, as well as the generation of large entangled
states, with 12 photon Greenberger–Horne–Zeilinger entanglement5 and
reconfigurable resource state generators of up to 8 qubits demonstrated6, 7.
However, these demonstrations are limited by photon generation efficiency,
with the mutliphoton data rate proportional to $Rp^{n}\eta^{n}$, for a $n$
photon system with probability $p$ of photon generation, transmission $\eta$
and clockspeed $R$. While $\eta$ can be improved by classical photonic
engineering, nonlinear sources have an intrinsically limited $p$, as to avoid
multiphoton emission, which scales as $O(p^{2})$, and remains unchanged since
the seminal quantum photonics experiments of the 1980s8, 9.

Today, there are two promising approaches to achieving a near-deterministic
single photon source. The first is based on solid-state emitters, such as
quantum dots 10, 11, NV-center defects in diamond12, or single molecules13,
14. To reach $p=1$, these systems must be engineered to deterministically emit
their single photons into a single guided mode. Over the last two decades, the
state of the art efficiency has risen to $p=57$% to optical fibre11 ($R=76$
MHz). Emission of near-indistinguishable photons from array of solid-state
emitters also remains challenging, as the wavelength and linewidth of the
emitted photons are subject to variation based on the precise arrangement of
emitter’s atomic environment (for example, a quantum dot’s shape and
interaction with the surrounding crystal lattice). Furthermore, while emission
at standard telecommunications wavelengths has been demonstrated15, telecom
quantum dots are less mature, with purity and efficiency metrics lagging
behind their counterparts in the near infra-red.

The second approach, known as photon source multiplexing, uses an array of
probabilistic nonlinear photon pair sources16, 17, and is the topic of this
paper. Here, the probability that at least one source among the array
generates a photon pair, $p_{h}$, tends to unity as the number of sources,
$m$, is increased. Upon pair generation, one of the photons is detected to
herald the presence of its partner, which can be dynamically switched to a
single output mode. In this way, if any source fires, a single photon is
routed to the output. Any loss experienced by the signal photon will reduce
the photon output probability, therefore the switching network is required to
have minimal loss, with unity probability reached with net transmission
$\eta\rightarrow 1$ and $m\rightarrow\infty$.

Fully integrated photon source multiplexing is a clear path to scaling today’s
quantum photonic integrated circuits6, 7, but remains a challenge due to the
lack of low-loss high-speed switching technology on integrated quantum
photonic platforms18. Furthermore, modern proposals for photonic quantum
computing architectures19, 20, 21 utilise multiplexing extensively—conditional
switching and measurement are crucial to both discrete and continuous variable
proposed quantum computer architectures and for reaching loss- and error-
tolerant thresholds22, 20, 19.

Today, chip- and fibre-based multiplexing systems have achieved enhancement
factors in the $1$–$3$ range, with $m\leq 4$ sources utilised23, 24, 25, 26,
27, though enhancement factors of up to $28$ have been realised with up to
$m=40$ time bins in bulk optics, benefiting from extremely low-loss Pockels
cells28 based switches. Repetition rates remain in the $R=0.5$–$15$ MHz
range23, 28 across these demonstrations. These relatively slow repetition
rates are limited by available switching technology, while spontaneous photon
generation has been demonstrated at clockspeeds up to 10 GHz29, offering a
potential increase of three orders of magnitude. Indeed, for fully-integrated
temporal multiplexing, GHz-speed switching is a requirement in order to reduce
loss and footprint. Achieving these high clockspeeds may be possible with
time- and wave-length division multiplexing of the pump30, while further
reductions can be achieved by multiplexing forward- and backward propagating
optical modes of each nonlinear source31.

In this paper, we present a temporally multiplexed single photon source based
on a low-loss fibre switch and loop architecture. Our resource efficient
design multiplexes $m=11$ temporally independent (time bin) nonlinear sources,
uses a single heralding detector, and achieves an enhancement heralded single
photon emission of $4.5\pm 0.5$.

##  2\. Multiplexed source

Photon pair sources can be multiplexed across any set of orthogonal photonic
modes. For example, any combination of spatial24, 25, temporal28, 25, 26, or
frequency30, 32 modes can be multiplexed, given the modes can be dynamically
switched with high fidelity and low loss. The loss of the mode conversion
(switching network) determines the maximum possible enhancement and maximum
number of effective sources, while the output probability is limited to the
transmission of the output photon through the switching network to its
application, $\eta_{s}$, known as the heralding efficiency. Since any optical
loss reduces output efficiency (as it also does with solid-state single photon
sources), a truly deterministic source $p=1$ is infeasible. However, loss-
tolerant bounds for quantum computing and repeater protocols exists in the
$>90$% range20, 19, 33.

Logarithmic tree structures provide the optimum enhancement for a given per-
pass switching loss, with each photon traversing the (lossy) switch
$\lceil\mathrm{log}_{2}(m)\rceil$ times34. These can be implemented in space
with a physical logarithimic tree structure, or with $\mathrm{log}_{2}(m)$
delay lines with delays of length $2^{n}\tau$. However, these topologies
require $O(m)$ switches and detectors, and thus inflate system size and cost
dramatically, and feature increased propagation loss in routing. Our temporal,
loop-based architecture16, 28, instead only requires a single switch and
detector, at the expense that the photons will instead traverse the switch
$(m+1)/2$ times on average, and achieves similar performance for small numbers
of bins.

In this work, we multiplex $m=11$ temporally distinct sources within an output
clock cycle, with photons generated by spontaneous four-wave mixing (SFWM) in
a $1.2$ cm long silicon waveguide (Fig. 1). Here, bright pump pulses from a
pulsed laser (1586 nm, 1 ps pulse duration, Calmar), are filtered with square-
shaped WDM filters (0.44 nm, Opneti), before being injected into the waveguide
via a grating coupler35. Pairs of photons in the pump are then converted into
pairs of photons at energy-conserving wavelengths symmetrically around the
pump by the $\chi^{(3)}$ nonlinearity of the silicon waveguide. The light is
then extracted and the pump removed with further WDM filters at $1591.26$ and
$1581.18$ nm, defining the herald and single photon output channels
respectively, before detection in superconducting nanowire single photon
detectors (SNSPDs, Quantum Opus), whose output is processed by a time-tagger
(qutools quTAG) to recover photon rate and arrival time data. The switch is an
unpolarised, $1$ MHz rate, $2\times 2$ fibre optical switch (Photonwares
NanoSpeed, proprietary), featuring an $8$ ns rise/fall time, sub-20 dB cross
talk, and an $80$ ns short pulse capability. A high-voltage electrical driving
system ($5$ V trigger) is also supplied.

To define a variable number of time-bins for multiplexing, pulses are picked
from the $16$ MHz repetition rate laser using a variable optical attenuator
(VOA). Here, $m$ pulses, spaced by $\tau=125$ ns, are selected in a 2 $\upmu$s
window, reducing the effective repetition rate to $8$ MHz. This defines the
output clock rate of the multiplexed source, $R=500$ kHz, as well as the
switch fibre loop length, which must delay the photons by $\tau$ per round
trip ($25$ m fibre). The switch and spliced fibre loop have a measured loss of
$1$ dB per round trip, which combined with the source output loss determines a
maximum photon probability enhancement of $4.5$, which is achieved in the
limit of low base pair generation probability, i.e. average photon number
$\mu\approx p\approx 0.01$.

An FPGA state machine, which is clocked by the laser, controls both the
multiplexing detect-switch-release logic, as well as generating the the pulse
picking signal. A fibre delay of $\delta=200$ ns is included in the signal
path before the fibre switch to compensate for the time taken to route the
heralding photon to the detector and the resulting electronic signal back to
the FPGA, which responds by activating the switch on the same clock cycle.

While the switch’s speed determines the output clockspeed, $R$, the pulse
response time of the switch determines the minimum size for individual time
bins, i.e. halving the pulse response time would allow twice as many time
bins, $m$, to be utilised for the same $R$. Meanwhile, the loss in the switch
determines the maximum attainable enhancement of the heralded single photon
probability, with increased the number of bins resulting in diminishing
returns toward a fixed limit corresponding to the net photon output
transmission $\eta_{s}$ (Fig. 2). Our system is capable of multiplexing up to
16 sources without sacrificing overall output clockspeed, $R$, though to
realise the potential of more than $m=11$ bins would require a switch with
increased transmission (improvements beyond $m=11$ were insignificant). In
principle, the FPGA state machine supports up to 64 multiplexed sources, at
the cost of reducing to $R$ to $125$ kHz.

Fig. 2: Photon probability enhancement by temporal-loop multiplexing. a.
Heralded photon probability and rate increase with the number of multiplexed
sources (time bins). The contribution from each source decreases exponentially
as photons in the $m^{\mathrm{th}}$ time bin must traverse the switch $m$
times, with $-1$ dB (79%) transmission per trip. b. The probability and rate
of detecting herald photons increases with number of multiplexed sources.
Solid lines are a calibrated numerical model28, using independently measured
loss and brightness parameters.

### 2.1. Results

We calibrate our base SFWM source to emit with a mean photon number of
$\mu=CR/(S_{h}S_{s})=0.009$ using the coincidence ($C$) to singles ($S$) ratio
in the low power limit, where $g^{(2)}(0)\approx\mu$ and find the single
photon transmissions $\eta_{s}=0.113$, $\eta_{h}=0.145$, including one pass
through the switch and loop for the signal photon. This low power is required
to suppress multiphoton noise below $1$%—required both for this calibration
and for scalable quantum photonic protocols—and furthermore enables us to
measure maximum attainable enhancement of the single photon probability.

We measure the singles and coincidence rate for $m=\\{1,\ldots,11\\}$ bins and
observe the contribution of each temporally independent photon pair source.
The contribution of each bin diminishes logarithmically, as photons generated
in bin $k$ must traverse the switch $k$ times and the loop $k-1$ times. Using
$m=11$ multiplexed time bin sources, we achieve an enhancement, $E$, of
$4.5\pm 0.5$, of the heralded photon detection probability from
$p\cdot\eta_{s}\eta_{h}=0.00014$ to $p_{m}=0.00063=Ep\cdot\eta_{s}\eta_{h}$.
Here, $p\approx\mu=0.009$ is the base generation probability in the source. We
also match this with a model using our measured loss and brightness parameters
for $m=\\{1,\dots,11\\}$ time bins28 (Fig. 2).

For a mean photon number $\mu=0.18$, as used in ref. 28, our model predicts an
output probability of $p_{m}=1.4$%, with an enhancement of $3.8$ at $m=11$
bins—though here the multiphoton emission rate is unacceptably high for most
quantum applications. Because we multiplex the first-generated photon pair,
increasing $\mu$ results in higher probability that photon pairs are generated
in earlier time bins, where the output photon experiences more loss, reducing
the relative enhancement.

Due to imperfect matching between the photon storage loop and the laser
repetition period, $\Delta_{\tau}=1.7$ ns, the output photons, which have a
coherence length $T_{c}\approx 5$ ps, are temporally distinguishable if they
were generated in different time bins. However, our simple and inexpensive
design lends itself to quantum key distribution systems where photon
indistinguishability is not required, but increasing photon probability boosts
secret key rate. Finally, our system comprises entirely of guided-mode optics,
and is therefore stable and miniaturisable with appropriate packaging of the
photonic integrated circuit. Thanks to this, the system is also modular. For
example, a high spectral purity fibre- or waveguide-based source2, 36, or more
efficient switch, could be swapped in to the system.

### 2.2. Outlook

Our device’s output clock speed, $R$, is limited to $500$ kHz by its driving
electronics (Photonwares). However, with spontaneous photon sources supporting
up to at least $R=10$ GHz29, up to three orders of magnitude improvement in
photon rate is possible. Future work, using an integrated switch, such as low-
loss thin-film lithium niobate (TFLN)37 will enable on-chip photon routing at
potentially at tens of GHz switching speeds with CMOS compatible voltages38,
enabling multiplexing within the recovery time of commercial SNSPDs. Here,
increasing the clock speed $R$, reduces propagation losses, and decreases
effective switch loss, resulting in higher attainable output efficiencies.
Furthermore, group delay can be precisely calculated and accounted for with
fine-tuning achieved via temperature stabilisation, for indistinguishable
photon generation. Hybrid silicon/TFLN integrated photonics39, achieved via
wafer bonding techniques, should enable photon sources and multiplexing switch
networks to be fabricated on a single chip, enabling orders of magnitude
improvements in efficiency and rate.

The heralded photon probability is held back by our system’s relatively low
heralding efficiency, and moderate switch loss. On-chip propagation and
grating coupler losses comprise the bulk of the loss (around $-6$ dB), with
the remaining loss occurring in the filters ($-0.6$ dB), detectors ($-1$ dB),
and fibre connections (around $-0.8$ dB), resulting in net transmissions
(including detection) of $\eta_{h}=0.145$, $\eta_{s}=0.113$, including one
pass through the switch for the signal photon. This results in a base heralded
photon output rate of $X_{b}=Rp\eta_{s}\eta_{h}=69\pm 7$ Hz, which is
increased to $X_{m}=312\pm 7$ Hz when $m=11$ time bins are multiplexed.

In our system, a zero-loss switch would yield $p_{m}=\eta_{s}^{\prime}$,
$X_{m}=71.5$ kHz as $m\rightarrow\infty$ where $\eta_{s}^{\prime}=0.143$ is
the transmission of the filtered photon without the switch. Here, any increase
in $R$, limited by the switch, provides a factor increase in multiplexed
output rate $X_{m}$. Decreasing the chip coupling losses, for example by
utilising modern, sub-dB grating couplers would enable this system to achieve
heralding efficiencies of greater than $\eta_{s}^{\prime}=0.6$ to the output
fibre, which gives an upper bound of multiplexed source output probability for
a lossless switch, $p_{m}=\eta_{s}^{\prime}$, $X_{m}=300$ kHz.

##  3\. Conclusion

Efficient single photon sources are a key challenge in today’s proposed
photonic quantum computer architectures19, 21 and quantum networks40, 20,
which utilise feed-forward and dynamic quantum circuits as a core feature. In
this work we have boosted the heralded photon emission of a silicon waveguide
photon pair source by dynamically detecting and routing photon pairs generated
in up to $11$ multiplexed time bins. Our system achieved a factor $4.5\pm 0.5$
enhancement of the base probability rate, with negligible change to the
$g^{(2)}(0)=0.01$. This work demonstrates the increased efficiency of
multiplexed quantum photonic hardware via a compact, guided-mode approach, on
the path to developing future fully-integrated devices for large-scale quantum
technology.

##  4\. Acknowledgements

The authors would like to thank Kjeld Dalgaard for designing the FPGA state
machine which controlled the experiment. We acknowledge funding from Villum
Fonden Young Investigator project QUANPIC (ref. 00025298) and Danish National
Research Foundation Center of Excellence, SPOC (ref. DNRF123).

##  5\. Contributions

JCA conducted the experiment, analysed the data and wrote the manuscript, JCA,
DB and YD designed the experiment, and YD managed the project.

## References

  * 1 Graffitti, F., Kundys, D., Reid, D. T., Brańczyk, A. M. & Fedrizzi, A.  Pure down-conversion photons through sub-coherence-length domain engineering.  _Quantum Science and Technology_ 2, 035001 (2017).
  * 2 Paesani, S. _et al._ Near-ideal spontaneous photon sources in silicon quantum photonics.  _Nature Commun._ 11, 2505 (2020).
  * 3 Zhong, H.-S. _et al._ Quantum computational advantage using photons.  _Science_ 370, 1460–1463 (2020).
  * 4 Zhong, H.-S. _et al._ Phase-programmable gaussian boson sampling using stimulated squeezed light.  _arXiv preprint arXiv:2106.15534_ (2021).
  * 5 Zhong, H.-S. _et al._ 12-photon entanglement and scalable scattershot boson sampling with optimal entangled-photon pairs from parametric down-conversion.  _Physical review letters_ 121, 250505 (2018).
  * 6 Adcock, J. C., Vigliar, C., Santagati, R., Silverstone, J. W. & Thompson, M. G.  Programmable four-photon graph states on a silicon chip.  _Nature communications_ 10, 1–6 (2019).
  * 7 Vigliar, C. _et al._ Error-protected qubits in a silicon photonic chip.  _Nature Physics_ 1–7 (2021).
  * 8 Shih, Y. & Alley, C. O.  New type of einstein-podolsky-rosen-bohm experiment using pairs of light quanta produced by optical parametric down conversion.  _Physical Review Letters_ 61, 2921 (1988).
  * 9 Wu, L.-A., Kimble, H., Hall, J. & Wu, H.  Generation of squeezed states by parametric down conversion.  _Physical review letters_ 57, 2520 (1986).
  * 10 Somaschi, N. _et al._ Near-optimal single-photon sources in the solid state.  _Nature Photonics_ 10, 340–345 (2016).
  * 11 Tomm, N. _et al._ A bright and fast source of coherent single photons.  _Nature Nanotechnology_ 16, 399–403 (2021).
  * 12 Wan, N. H. _et al._ Large-scale integration of artificial atoms in hybrid photonic circuits.  _Nature_ 583, 226–231 (2020).
  * 13 Khasminskaya, S. _et al._ Fully integrated quantum photonic circuit with an electrically driven light source.  _Nature Photonics_ 10, 727 (2016).
  * 14 Toninelli, C. _et al._ Single organic molecules for photonic quantum technologies.  _Nature Materials_ 1–14 (2021).
  * 15 Müller, T. _et al._ A quantum light-emitting diode for the standard telecom window around 1,550 nm.  _Nature communications_ 9, 1–6 (2018).
  * 16 Pittman, T. B., Jacobs, B. C. & Franson, J. D.  Single photons on pseudodemand from stored parametric down-conversion.  _Physical Review A_ 66, 042303 (2002).
  * 17 Migdall, A. L., Branning, D. & Castelletto, S.  Tailoring single-photon and multiphoton probabilities of a single-photon on-demand source.  _Physical Review A_ 66, 053805 (2002).
  * 18 Adcock, J. _et al._ Advances in quantum machine learning.  _arXiv preprint arXiv:1512.02900_ (2015).
  * 19 Bartolucci, S. _et al._ Fusion-based quantum computation.  _arXiv preprint arXiv:2101.09310_ (2021).
  * 20 Bartolucci, S. _et al._ Creation of entangled photonic states using linear optics.  _arXiv preprint arXiv:2106.13825_ (2021).
  * 21 Bourassa, J. E. _et al._ Blueprint for a scalable photonic fault-tolerant quantum computer.  _Quantum_ 5, 392 (2021).
  * 22 Gimeno-Segovia, M., Shadbolt, P., Browne, D. E. & Rudolph, T.  From three-photon greenberger-horne-zeilinger states to ballistic universal quantum computation.  _Physical review letters_ 115, 020502 (2015).
  * 23 Ma, X.-s., Zotter, S., Kofler, J., Jennewein, T. & Zeilinger, A.  Experimental generation of single photons via active multiplexing.  _Physical Review A_ 83, 043814 (2011).
  * 24 Collins, M. J. _et al._ Integrated spatial multiplexing of heralded single-photon sources.  _Nature communications_ 4 (2013).
  * 25 Mendoza, G. J. _et al._ Active temporal and spatial multiplexing of photons.  _Optica_ 3, 127–132 (2016).
  * 26 Xiong, C. _et al._ Active temporal multiplexing of indistinguishable heralded single photons.  _Nature communications_ 7 (2016).
  * 27 Francis-Jones, R. J., Hoggarth, R. A. & Mosley, P. J.  All-fiber multiplexed source of high-purity single photons.  _Optica_ 3, 1270–1273 (2016).
  * 28 Kaneda, F. & Kwiat, P. G.  High-efficiency single-photon generation via large-scale active time multiplexing.  _Science advances_ 5, eaaw8586 (2019).
  * 29 Zhang, X. _et al._ High repetition rate correlated photon pair generation in integrated silicon nanowires.  _Applied Optics_ 56, 8420–8424 (2017).
  * 30 Zhang, X. _et al._ Enhancing the heralded single-photon rate from a silicon nanowire by time and wavelength division multiplexing pump pulses.  _Optics letters_ 40, 2489–2492 (2015).
  * 31 Xiong, C. _et al._ Bidirectional multiplexing of heralded single photons from a silicon chip.  _Optics letters_ 38, 5176–5179 (2013).
  * 32 Joshi, C., Farsi, A., Clemmen, S., Ramelow, S. & Gaeta, A. L.  Frequency multiplexing for quasi-deterministic heralded single-photon sources.  _Nature communications_ 9, 1–8 (2018).
  * 33 Morley-Short, S., Gimeno-Segovia, M., Rudolph, T. & Cable, H.  Loss-tolerant teleportation on large stabilizer states.  _Quantum Science and Technology_ (2018).
  * 34 Bonneau, D., Mendoza, G. J., O’Brien, J. L. & Thompson, M. G.  Effect of loss on multiplexed single-photon sources.  _New Journal of Physics_ 17, 043057 (2015).
  * 35 Ding, Y., Peucheret, C., Ou, H. & Yvind, K.  Fully etched apodized grating coupler on the SOI platform with $-0.58$ dB coupling efficiency.  _Optics letters_ 39, 5348–5350 (2014).
  * 36 Lugani, J., Francis-Jones, R. J., Boutari, J. & Walmsley, I. A.  Spectrally pure single photons at telecommunications wavelengths using commercial birefringent optical fiber.  _Optics express_ 28, 5147–5163 (2020).
  * 37 He, M. _et al._ High-performance hybrid silicon and lithium niobate mach–zehnder modulators for 100 gbit s- 1 and beyond.  _Nature Photonics_ 13, 359–364 (2019).
  * 38 Wang, C. _et al._ Integrated lithium niobate electro-optic modulators operating at CMOS-compatible voltages.  _Nature_ 562, 101–104 (2018).
  * 39 He, M. _et al._ High-performance hybrid silicon and lithium niobate mach–zehnder modulators for 100 Gbit $\mathrm{s}^{-1}$ and beyond.  _Nature Photonics_ 13, 359–364 (2019).
  * 40 Azuma, K., Tamaki, K. & Lo, H.-K.  All-photonic quantum repeaters.  _Nature communications_ 6, 6787 (2015).