# On the eternal non-Markovianity of non-unital quantum channels

Shrikant Utagi<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Technology Madras, Chennai - 600036, India.  Subhashish Banerjee
<EMAIL_ADDRESS>Indian Institute of Technology Jodhpur, Jodhpur -
342037, India.  R. Srikanth<EMAIL_ADDRESS>Theoretical Sciences Division,
Poornaprajna Institute of Scientific Research, Bengaluru - 562164, India.

###### Abstract

The eternally non-Markovian Pauli channel is an example of a unital channel
characterized by a negative decay rate for all time $t>0$. Here we consider
the problem of constructing an analogous non-unital channel, and show in
particular that a $d$-dimensional generalized amplitude damping (GAD) channel
cannot be eternally non-Markovian when the non-Markovianity originates solely
from the non-unital part of the channel. We study specific ramifications of
this result for qubit GAD. Specifically, we construct a quasi-eternally non-
Markovian qubit GAD channel, characterized by a time $t^{\ast}>0$, such that
the channel is non-Markovian only and for all time $t>t^{\ast}$. We further
point out that our negative result for the qudit GAD channel, namely the
impossibility of the eternal non-Markovian property, does not hold for a
general qubit or higher-dimensional non-unital channel.

##  I Introduction

Open quantum systems is the study of the dynamics of the system of interest
taking into account the effect of its ambient environment and an important
aspect of this is the study of non-Markovian (NM) dynamics [1, 2, 3], which
has gained increasing attention recently in the quantum information theory
community [4, 5, 6, 7, 8, 9, 10]. A number of definitions of non-Markovianity
have been proposed namely, based on distinguishability [11], completely
positive (CP-) divisibility, fidelity [12], negativity of the decay rate [13,
14], quantum temporal correlations [15, 16] and a more general definition
based on deviation from semigroup structure [17], to name a few.

The notion of quantum non-Markovianity has always been debated, see for
instance Ref. [18, 19, 20, 21, 22, 7]. Its revived interest is due to the
advent of newer experiments with finer engineered baths, which highlight the
fact that there exists no universal definition or criterion which faithfully
characterizes, detects or measures quantum non-Markovianity. Most of the
developments focused on the two-time dynamical maps and the divisibility
criterion based on the two-time correlation functions in determining the
Markovian status of open system dynamics, See Ref. [4, 5, 6]. A stronger
notion of quantum Markovianity has recently been developed, that is based on
an operational point of view [20, 19, 22, 21, 18, 23] and eschews the use of
the concept of the not-completely positive (NCP) intermediate map. Moreover,
this approach meets the status of necessity and sufficiency of characterizing
non-Markovianity when multi-time correlations correlations are taken into
account. For the present purpose of this work, we restrict ourselves to the
divisibility criteria based on the two-time correlations and two-time
dynamical map formalism. In particular, we will be concerned with the eternal
CP-indivisibility and P-indivisibility of channels whose non-Markovianity
originates solely from the non-unital part, which we definite shortly.

An open system evolution is generally given by the famous Gorini-Kossakowski-
Lindblad (GKSL)-like master equation [24, 25]:

$\displaystyle\dot{\rho}=\mathcal{L}(t)[\rho]=-i[H,\rho]-\sum_{j}\gamma_{j}(t)\bigg{(}A_{j}(t)\rho
A_{j}^{\dagger}(t)-\frac{1}{2}\\{A_{j}^{\dagger}(t)A_{j}(t),\rho\\}\bigg{)},$
(1)

where $\gamma_{j}(t)$ are the canonical decay rates [14] and $A_{j}(t)$ are
the jump operators. Assuming that the system state $\rho_{S}(t=0)$ and
environment initially fixed in the state $\rho_{E}=\ket{e_{0}}\bra{e_{0}}$ is
a product i.e., $\rho_{SE}=\rho_{S}(0)\otimes\rho_{E}$, one can deduce a
completely positive trace preserving (CPTP) dynamical map on the system state
via the operator-sum form [25, 2]:

$\mathcal{E}(t)[\rho]=\sum_{j}\bra{e_{j}}U(\rho_{S}(0)\otimes\ket{e_{0}}\bra{e_{0}})U^{\dagger}\ket{e_{j}}=\sum_{j}E_{j}(t)\rho_{S}E_{j}^{\dagger}(t),$

where $K_{j}(t)\equiv\bra{e_{j}}U\ket{e_{0}}$ are called the Kraus operators,
where $U$ is joint system-environment unitary. The map $\mathcal{E}$ is in
fact the solution to the master equation
$\dot{\mathcal{E}}=\mathcal{L}\mathcal{E}$, which, for all $t\geq s\geq
t_{0}$, is given by
$\mathcal{E}(t,t_{0})=\mathcal{T}\exp\\{\int_{t_{0}}^{t}ds\mathcal{L}(s)\\}$,
where $\mathcal{T}$ is the time-ordering operator. However, it is well-known
that when the initial factorization assumption is relaxed, the ensuing
dynamical map need not be CP [26, 27, 28]. Moreover, if a map could be written
down, one must sacrifice either complete positivity or linearity. Therefore,
we stress that the dynamical maps defined throughout this work are obtained
based on the initial factorization assumption.

Now, we are in a position to define the concept of divisibility of a map.

A CPTP qubit map $\mathcal{E}(t_{2},t_{0})$ (a qubit quantum channel) is said
to be CP-divisible if for arbitrary intermediate time $t_{1}$ it can be
written as

$\displaystyle\mathcal{E}(t_{2},t_{0})=\mathcal{E}(t_{2},t_{1})\mathcal{E}(t_{1},t_{0}),$
(2)

such that $\mathcal{E}(t_{2},t_{1})$ is a CP-map. The Choi-Jamiolkowski matrix
[29] or B-matrix [30] for a map $\mathcal{E}(t_{2},t_{1})$ is given by

$\displaystyle\mathcal{\chi}=(\mathcal{E}(t_{2},t_{1})\otimes
I)[\ket{\Psi}\\!\bra{\Psi}],$ (3)

where $\ket{\Psi}=\frac{1}{\sqrt{2}}(\ket{00}+\ket{11})$ is a maximally
entangled state. When $\chi$ is negative, the channel
$\mathcal{E}(t_{2},t_{0})$ is CP-indivisible, in the sense that it is a CPTP
map that cannot be written as composition of CP maps as in Eq. (2). For
bijective maps, the above requirement is equivalent to the following condition
for CP-divisibility

$\displaystyle\frac{d}{dt}\|[\mathbb{I}_{d+1}\otimes\mathcal{E}(t)](\rho_{1}-\rho_{2})\|_{1}\leq
0$ (4)

being always satisfied for given any two density operators $\rho_{1}$ and
$\rho_{2}$ in the space of bounded operators of the bipartite system
$\mathcal{B}(H_{d+1}\otimes H_{d})$ [31]. Here, $\|A\|_{1}={\rm
Tr}(\sqrt{A^{\dagger}A})$ is the trace norm of an operator $A$. The channel
$\mathcal{E}(t_{2},t_{0})$ satisfies P-divisibility if

$\displaystyle\frac{d}{dt}\|\mathcal{E}(t)[\rho_{1}-\rho_{2}]\|_{1}\leq 0$ (5)

for any initial pair of states $\rho_{1}$ and $\rho_{2}$. Note that when the
condition (5) breaks down, the intermediate map $\mathcal{E}(t_{2},t_{1})$
would no longer be positive, in the sense that it outputs a negative state.
However, the full map $\mathcal{E}(t_{2},t_{0})$ is still CPTP but no longer
P-divisible.

The notion of an “eternally non-Markovian” (ENM) Pauli channels was proposed
in [14] via the canonical master equation given by

$\displaystyle\mathcal{L}(t)[\rho]=\sum_{j=1}^{3}\gamma_{j}(t)(\sigma_{j}\rho\sigma_{j}-\rho)$
(6)

where $\sigma_{j}$ are the Pauli operators, and the canonical decay rates
$\gamma_{j}(t)$ are given by $\gamma_{1}=\gamma_{2}=\frac{c}{2}$ and
$\gamma_{3}(t)=-\frac{c}{2}\tanh t$, where $c$ is some real constant. This
particular dynamics corresponds to CPTP map
$\mathcal{E}(t)\rho=\sum_{j=1}^{3}E_{j}\rho E_{j}^{\dagger}$, with the Kraus
operators given by

$E_{1}=\sqrt{2k}\sigma_{I},\leavevmode\nobreak\ \leavevmode\nobreak\
E_{2}=\sqrt{k}\sigma_{X};\leavevmode\nobreak\ \leavevmode\nobreak\
E_{3}=\sqrt{k}\sigma_{y}$ (7)

where

$k\equiv\frac{1}{4}(1-e^{-ct}).$ (8)

For Pauli channels, it is straightforward to arrive at Eq. (7) from Eq. (6)
from the so-called $H$-matrix method [32, 33]. This example was later shown to
arise under the mixing of two Pauli semigroups [34], providing an example of
the fact that the space of CP-divisible Pauli channels is not convex [35].
This is an instance in which the criterion Eq. (5) for witnessing non-
Markovianity fails [14]. Interestingly, in [36] the authors studied whole
families of eternal and quasi-eternal unital channels. ENM evolution arising
from a convex combination of Markovian semigroups has recently been extended
to the context of generalized Pauli channels in dimension $k$ [37], where it
is shown that $(d-1)^{2}$ out of $d^{2}-1$ decay rates can be always negative.

The examples of ENM evolution given above all have two broad restrictions: (1)
the concept of non-Markovianity here refers to CP-indivisibility, and not the
stronger condition of P-indivisibility; (2) they consider only unital–
specifically Pauli– channels. Here we wish to address the 2nd of these points.
So far as known to us, there are only a few works addressing similar questions
about non-unital channels. Here we undertake to do this, studying in detail
the generalized amplitude damping (GAD) channels [2, 38].

Before embarking on our results, we present a notion of non-unital non-
Markovianity as below. Any $d$-dimensional dynamical map may be given a matrix
representation as $\mathcal{E}(t)\rightarrow F(t)$ as

$F(t)=\left(\begin{array}[]{c|c}1&{\bm{0}}_{1\times(d^{2}-1)}\\\
\hline\cr\bm{\tau}&M\end{array}\right),$ (9)

where $\bm{\tau}\in\mathbb{R}^{d^{2}-1}$ and $M$ is a $d^{2}-1\times d^{2}-1$
real matrix. Given the state
$\rho=\frac{1}{d}(\mathbb{I}+\sum_{i=1}^{d^{2}-1}r_{i}G_{i}$), where $G_{i}$
are orthonormal basis given by the generalized Gell-Mann matrices with
$G_{0}=\frac{\mathbb{I}}{\sqrt{d}}$ and the Hermitian $G_{i}$ where
$i\in\\{1,.....d^{2}-1\\}$, then the map $F_{ij}(t):={\rm
Tr}(G_{i}\mathcal{E}[G_{j}])$, which transforms a Bloch vector $\bm{r}$ as

$\displaystyle\bm{r}\rightarrow\bm{r}^{\prime}(t)=M(t)\bm{r}(0)+\bm{\tau}(t).$
(10)

If $\bm{r_{1}}$ and $\bm{r_{2}}$ are two Bloch vectors corresponding to the
states $\rho_{1}$ and $\rho_{2}$, then from Eq. (9) and (10), the
transformation $\mathcal{E}(t)[\rho_{1}-\rho_{2}]$ corresponds to
$M(t)[\bm{r_{1}}-\bm{r_{2}}]$. Therefore, the trace distance fails to witness
non-Markovianity originating solely from $\bm{\tau}$. However, the full non-
Markovianity, in the sense of P-indivisibility is encoded in the map $F(t)$.
Interestingly, the divisibility property of $\mathcal{E}(t)$ carried over as
$F(t)=F(t,s)F(s)$ leads to a non-trivial relationship between unital and non-
unital parts of the channel:

$\displaystyle
M(t)=M(t,s)M(s)\quad\text{and}\quad\bm{\tau}(t)=\bm{\tau}(t,s)+M(t,s)\bm{\tau}(s).$
(11)

Therefore, the full P-indivisibility is captured by Eq. (11).

###### Definition 1

When non-Markovianity originates solely from the non-unital part, i.e., from
$\bm{\tau}(t)$, then we call it non-unital non-Markovianity and such a channel
will be called purely non-unital non-Markovian channel.

Given this notion, we ask: is eternal non-unital non-Markovianity, as per the
Definition (1), possible? We answer this in the negative for a class of non-
unital channels. However, we will later show that for a class of phase-
covariant channels and higher dimensional non-unital channels, eternal non-
Markovianity is still possible, nonetheless such a non-Markovianity originates
from the unital part of the channel.

This is work is organized as follows. In Section II we introduce a $d$-level
GAD channel for $d$-dimensional systems and show that it cannot be eternally
non-Markovian. In Section III, we further provide rigorous proof for the case
of qubit GAD channel and show that the utmost one can achieve in this
direction is a quasi-eternal non-Markovian GAD channel, of which we explicitly
construct a qubit example. In Section IV, we point out that this negative
result does not extend to non-unital channels in general, by providing an
example of a non-unital channel that is ENM. We present our conclusions in
Section V.

##  II Generalized amplitude damping and the impossibility of eternal non-
Markovianity

###  II.1 Qudit generalized amplitude damping

A well-known exemplar of a qubit non-unital channel is the GAD channel [2,
39]. In this work, generalizing the well-known $d$-level amplitude damping
channel [40, 41], we introduce a $d$-level GAD channel $\mathcal{E}_{G}$ for
$d$-dimensional systems given by

$\mathcal{E}_{G}=\sum_{l=0}^{d-1}p_{l}\mathcal{E}_{l},$ (12)

where $p_{l}$ is the probability with which the $l^{\rm th}$ damping channel
$\mathcal{E}_{l}$ acts on a $d$-dimensional system, along with the constraint
that $\sum_{l}p_{l}=1$. Here $\mathcal{E}_{l}$ damps an arbitrary input state
$\rho$ to the state $\ket{l}\\!\bra{l}$ with probability $\lambda$, and is
defined by the Kraus representation
$\mathcal{E}_{l}[\rho]=\sum_{j=0}^{d-1}E_{l,j}\rho E_{l,j}^{\dagger}$

with the Kraus operators in terms of damping parameter $\lambda$ as

$\displaystyle E_{l,l}$
$\displaystyle=\ket{l}\\!\bra{l}+\sqrt{1-\lambda}\sum_{\begin{subarray}{c}k=0\\\
k\neq l\end{subarray}}^{d-1}\ket{k}\\!\bra{k};$ $\displaystyle E_{l,j}$
$\displaystyle=\sqrt{\lambda}\ket{l}\\!\bra{j},\quad\quad(l\neq j).$ (13)

A key point about the damping channels is that they have a unique fixed point
$\mathfrak{F}$. For the amplitude damping channels $\mathcal{E}_{l}$ it is
$\ket{l}\\!\bra{l}$, and for the higher-dimensional GAD channel the fixed
point is $\mathfrak{F}\equiv\sum_{l=0}^{d-1}p_{l}\ket{l}\\!\bra{l}$.

We have the following result.

###### Theorem 1

A non-unital GAD channel in any finite dimension $d$ that is eternally CP-
indivisible in the non-unital part is impossible.

Proof: In the Markovian case, any initial state is driven monotonically
towards the fixed point $\mathfrak{F}$, fulfilling Eq. (4) for all pair of
states $\rho_{1}$ and $\rho_{2}$. If this CP-divisibility condition is
violated at time $t=0$, then there is a pair of states such that

$\displaystyle\frac{d}{dt}\|[\mathbb{I}\otimes\mathcal{E}(t)](\rho_{1}-\rho_{2})\|_{t=0}>0.$
(14)

With the idea of maximizing the l.h.s in Eq. (14), we let
$\rho_{2}\equiv\mathbb{I}\otimes\mathfrak{F}$ and
$\rho_{1}\equiv\mathbb{I}\otimes\rho_{S}$ for an arbitrary system state
$\rho_{S}$.

It can be shown that the action of the GAD channel Eq. (12) is described by

$\mathcal{E}(\rho_{S})[t]=\lambda(t)\mathfrak{F}+\\{(1-\lambda(t))\mathfrak{F}+\sqrt{1-\lambda(t)}\mathfrak{G}\\},$
(15)

where $\mathfrak{G}$ is a matrix consisting only of off-diagonal terms
dependent on the initial state. In the Markovian case, $\lambda(t)$ increases
monotonically from 0 to 1, and thus $\mathcal{E}(\rho_{S})[t]$ evolves
monotonically from the initial state to $\mathfrak{F}$. Note that this
corresponds manifestly to non-unital behavior. In specific,
$\mathcal{E}(\mathfrak{F})=\mathfrak{F}$ at all times $t$, as required by the
definition of a fixed point. We have from Eq. (15) that

$\|\mathcal{E}(\rho_{S})[t]-\mathfrak{F}\|=\sqrt{1-\lambda(t)}\|\mathfrak{G}\|.$
(16)

Given that GAD constitutes a CP map, $\lambda(t)\geq 0$ such that
$\lambda(0)=0$. Thus, we have that $\dot{\lambda}|_{t=0}\geq 0$. Then it
follows from Eq. (16) that

$\displaystyle\frac{d}{dt}\|\mathcal{E}({\rho}_{S})-\mathfrak{F}\|_{t=0}\leq
0,$ (17)

contrary to the implication of Eq. (14). $\square$ In the following, the above
result for the general case of qudits will be elucidated in the case of
qubits.

###  II.2 Case of qubits

In the case of qubits, we can illustrate the above results for maps in terms
of implications for the decoherence rates.

$\displaystyle E_{1}$ $\displaystyle=\sqrt{1-p}\left[\begin{array}[]{cc}1&0\\\
0&\sqrt{1-\lambda}\end{array}\right];\;E_{2}=\sqrt{1-p}\left[\begin{array}[]{cc}0&\sqrt{\lambda}\\\
0&0\end{array}\right];$ (22) $\displaystyle E_{3}$
$\displaystyle=\sqrt{p}\left[\begin{array}[]{cc}\sqrt{1-\lambda}&0\\\
0&1\end{array}\right];\;E_{4}=\sqrt{p}\left[\begin{array}[]{cc}0&0\\\
\sqrt{\lambda}&0\end{array}\right],$ (27)

where $\lambda,p\in[0,1]$.

A key point is that the form Eq. (27) constrains the Markovian behavior of the
channel, even without knowing details of the time dependence of the
parameters.

###### Lemma 1

The damping parameter $\lambda(t)$ of the GAD channel is such that

$\lim_{t\rightarrow 0+}\dot{\lambda}\geq 0,$ (28)

meaning that for sufficiently small time $t$, the variable $\lambda(t)$ cannot
decrease.

Proof: The functional form of Eq. (27) implies that

$\displaystyle\lambda(t=0)=0\;\;\text{and}\;\;\forall_{t\geq
0}\leavevmode\nobreak\ \lambda(t)\geq 0,$ (29)

where the first equation is a consequence of requiring that a channel must
reduce to the identity map $I$ at $t=0$, and the second that of requiring that
the map specified by Eq. (27) must be complete. Eq. (29) implies Eq. (28).
$\square$

It is well known that trace-distance can witness P-indivisibility. Our first
result is to show that GAD cannot manifest eternal P-indivisibility. To this
end, for the above map we have the affine representation

$\displaystyle\vec{r}^{\prime}(t)=M(t)\vec{r}(0)+\vec{\tau}(t)$ (30)

where $M$ is a $3\times 3$ real matrix which modifies the Bloch vector
$\vec{r}=(r_{1},r_{2},r_{3})^{T}$, and
$\vec{\tau}=(\tau_{1},\tau_{2},\tau_{3})^{T}$ is the shift vector representing
the displacement of Bloch vector towards the fixed point of the non-unital
channel under consideration. Now, for a GAD channel we have

$\displaystyle M$
$\displaystyle=\left(\begin{array}[]{ccc}\sqrt{1-\lambda}&0&0\\\
0&\sqrt{1-\lambda}&0\\\ 0&0&1-\lambda\\\
\end{array}\right);\quad\vec{\tau}=\left(\begin{array}[]{c}0\\\ 0\\\
(1-2p)\lambda\\\ \end{array}\right).$ (37)

The non-Markovianity can arise from the unital part (the matrix $M$) or the
non-unital part (the shift $\vec{\tau}$). For the most part in this paper, we
set the notation $p(t)=p$ and $\lambda(t)=\lambda$ for simplicity.

The block diagonal representation of a dynamical map provides a
straightforward way of connecting trace distance with the decay rates [14],
hence we make use of this feature here. Eq. (30) obeys the canonical master
equation in the Block diagonal representation

$\displaystyle\frac{d\vec{r}(t)}{dt}=D(t)\vec{r}(t)+\vec{\mu}(t)$ (38)

where $D(t)$ and $\vec{\mu}(t)$ are called damping matrix and drift vector
respectively. In relation with Eq. (30) we have
$\frac{dM(t)}{dt}=D(t)M(t)\quad\implies\quad M(t)=e^{\int_{t_{0}}^{t}D(s)ds}.$
The channel $\mathcal{E}$ is CP-divisible when $\text{Tr}[D(t)]<0$ and
P-divisible when $h_{\rm max}(D(t)+D(t)^{\rm T})<0$ [14], where $h_{\rm
max}(X)$ is the largest eigenvalue of the matrix $X$.

###### Theorem 2

A qubit GAD channel that is eternally non-Markovian as witnessed by the trace-
distance measure is impossible.

Proof: For the qubit GAD channel Eq. (27), we have:

$D(t)=\begin{pmatrix}\frac{-\dot{\lambda}}{2(1-\lambda)}&0&0\\\
0&\frac{-\dot{\lambda}}{2(1-\lambda)}&0\\\
0&0&\frac{-\dot{\lambda}}{1-\lambda}\end{pmatrix}.$ (39)

At time $t\rightarrow 0$, we have $h_{\rm max}(D(t)+D(t)^{\rm
T})=\frac{-\dot{\lambda}}{1-\lambda}$, which by virtue of Lemma 1 is negative.
$\square$

We also have that in the limit $t\rightarrow 0$,
$\textrm{Tr}[D(t)]=\frac{-2\dot{\lambda}}{1-\lambda}<0$ on account of Eq.
(28). Therefore, the Bloch volume measure [42, 14] also is unable to indicate
an ENM property.

One might suppose that even though the GAD channel cannot be eternally
P-indivisible, it could still be eternally CP-indivisible. Our next result is
that even this is ruled out, because the channel must start off as Markovian.

###### Theorem 3

A qubit GAD channel that is eternally CP-indivisible is impossible.

Proof: For the qubit GAD channel, the canonical master equation is

$\displaystyle\dot{\rho}(t)$ $\displaystyle=\mathcal{L}[\rho(t)]$
$\displaystyle=\gamma_{1}(t)[\sigma_{-}\rho(t)\sigma_{+}-\frac{1}{2}\\{\sigma_{+}\sigma_{-},\rho(t)\\}]$
$\displaystyle+\gamma_{2}(t)[\sigma_{+}\rho(t)\sigma_{-}-\frac{1}{2}\\{\sigma_{-}\sigma_{+},\rho(t)\\}]$
(40)

where $\sigma_{+}=\ket{1}\bra{0}$ and $\sigma_{-}=\ket{0}\bra{1}$ are the
raising and lowering operators. The canonical decay rates are given by:

$\displaystyle\gamma_{1}(t)$
$\displaystyle=\lambda\dot{p}+\frac{\dot{\lambda}p}{1-\lambda},$
$\displaystyle\gamma_{2}(t)$
$\displaystyle=\lambda\dot{q}+\frac{\dot{\lambda}q}{1-\lambda},$ (41)

where $p+q=1$. We thus have
$\gamma_{1}+\gamma_{2}=\frac{\dot{\lambda}}{1-\lambda}$ in as much as
$\dot{p}=-\dot{q}$. Eqs. (41) entail that

$\displaystyle\lim_{t\rightarrow 0+}\gamma_{1}(t)$
$\displaystyle=\dot{\lambda}p,$ $\displaystyle\lim_{t\rightarrow
0+}\gamma_{2}(t)$ $\displaystyle=\dot{\lambda}(1-p).$ (42)

It then follows by virtue of Eqs. (28) and (42), and noting that $p$ is a
probability, that $\lim_{t\rightarrow 0+}\gamma_{j}\geq 0$ ($j=1,2$).
$\square$

##  III Quasi-eternally non-Markovian GAD channel

Given the impossibility of an ENM qubit GAD channel, the utmost that can be
achieved in this direction is quasi-eternal non-Markovianity, whereby one of
the rates becomes negative after a finite time $t_{0}$, and remains negative
forever thereafter. The closest approximation to eternal non-Markovianity
would then be by letting $t_{0}\rightarrow 0$. An example of this sort is
discussed below. To this end, we introduce an eternal “unmixing” effect, i.e.,
having $\dot{p}(t)<0$ for all times. In that spirit, in Eq. (27), let

$\displaystyle p(t)=\frac{e^{-mt}}{n}\;\;;\;\;\lambda(t)=1-e^{-\nu t}$ (43)

where $m,n$ and $\nu$ are some real positive constants. From Eq. (41), we find
that

$\displaystyle\gamma_{1}(t)$
$\displaystyle=\frac{1}{n}e^{-mt}\left(m\left(e^{-\nu t}-1\right)+\nu\right)$
(44a) $\displaystyle\gamma_{2}(t)$ $\displaystyle=\nu-\gamma_{1}(t)$ (44b)

are the canonical decay rates, which are plotted in Figure (1) for the values
$m=3,n=2$ and $\nu=1$.

Decoherence rate $\gamma_{2}$ remains positive throughout. The initial time
when rate $\gamma_{1}$ turns negative is found using Eq. (44a) to be

$t^{\ast}=\frac{1}{\nu}\log\left(\frac{m}{m-\nu}\right).$ (45)

Note that since $m>\nu$, thus there is an initial interval given by $0\leq
t<t^{\ast}$, during which $\gamma_{1}$ is necessarily positive. In Eq. (44a),
the term $(m(e^{-\nu t}-1)+\nu)\rightarrow\nu-m$, in the limit
$t\rightarrow\infty$. Therefore the decay rate approaches $0^{-}$, provided
$m>\nu$. This provides the condition for the quasi-eternal non-Markovianity of
the channel. The symmetry between the two rates is a consequence of Eq. (44b),
and more generally, Eq. (41), and is reflected in Figure 1.

Figure 1: (Color online) The decay rates for the quasi-eternal non-Markovian
GAD channel described by the Lindblad rates Eq. (44), with values $m:=3,n:=2$
and $\nu:=1$. Whilst rate $\gamma_{2}$ is always positive, rate $\gamma_{1}$
is positive until time $t^{\ast}$, Eq. (45), before becoming negative, and
thereafter remains negative throughout thereafter. It is a consequence of Eq.
(45) that $t^{\ast}$ can be made arbitrarily close to 0, but not 0 itself, by
virtue of Theorem 1. The symmetry between the two decay rates is a consequence
of the fact, which follows from Eqs. (44a) and 44b), that
$\gamma_{1}(t)+\gamma_{2}(t)=\nu$.

This channel by design is P-divisible in the sense that the intermediate map
$F(t,s)$, according to Eq. (11), is positive at all times and thus will be
indicated to be non-Markovian neither by a distinguishability criterion such
as trace distance [43], nor by a criterion based on entropic quantities [44].
On the other hand, the occurrence of negative $\gamma_{1}$ implies that it is
CP-indivisible. The non-Markovianity of the channel may be quantified using a
measure due to Hall-Cresser-Li-Anderson [14], which is given by
$\mathcal{\xi}_{\rm\small HCLA}:=-\int_{t:\gamma(t)<0}\gamma(t)dt$. In the
context of this paper, $\gamma_{1}(t)<0$ for all $t>t^{\ast}$, hence we find
that

$\displaystyle\mathcal{\xi}_{\rm\small HCLA}$
$\displaystyle:=-\int_{{t^{\ast}}}^{\infty}\gamma_{1}(t)dt$
$\displaystyle=\frac{\nu}{(\nu+m)n}\left(-\frac{m}{\nu-m}\right)^{-\frac{\nu+m}{\nu}}.$
(46)

where $t^{\ast}$ is given in by Eq. (45). The HCLA measure is plotted in Fig.
(2).

Figure 2: (Color online) The HCLA measure of QENMGAD channel (27) against $m$.

##  IV Non-unital channels and eternal non-Markovianity

The impossibility of ENMity of GAD provokes the question of whether this is
generic for all qubit non-unital channels. Yet, a simple example suffices to
show that this is not the case. Consider the non-unital channel generated by
linearly combining the generators the above Pauli ENM channel and an amplitude
damping (AD) channel. An example would be:

$\displaystyle\mathcal{L}(\rho)$
$\displaystyle=[\sigma_{X}\rho\sigma_{X}-\rho]+[\sigma_{Y}\rho\sigma_{Y}-\rho]-\tanh(t)[\sigma_{Z}\rho\sigma_{Z}-\rho]$
$\displaystyle+\gamma(t)[\sigma_{-}\rho\sigma_{+}-\frac{1}{2}\\{\sigma_{+}\sigma_{-},\rho\\}],$
(47)

where $\gamma(t)$ is arbitrary decoherence rate for an AD channel. The action
of a unital channel leaves the identity operator unchanged, whence
$\dot{\rho}=\mathcal{L}_{\rm unital}(\rho)=0$ if $\rho\propto\mathbb{I}$. The
non-unitality of the channel described in Eq. (47) follows, noting that
$L(\mathbb{I})\neq 0$. Employing the notation
$\sigma_{\mp}\equiv\frac{1}{{2}}(\sigma_{X}\pm i\sigma_{Y})$, the above
equation can be rewritten as

$\displaystyle\mathcal{L}(\rho)=\sum_{j,k}d_{j,k}(\sigma_{j}\rho\sigma_{k}-\frac{1}{2}\\{\sigma_{k}\sigma_{j},\rho\\}),$
(48)

where the decoherence matrix ${\bf d}\equiv\\{d_{j,k}\\}$ [14] is given by

${\bf d}=\begin{pmatrix}1+\frac{\gamma}{2}&-i\frac{\gamma}{4}&0\\\
i\frac{\gamma}{4}&1+\frac{\gamma}{2}&0\\\ 0&0&-\tanh(t)\end{pmatrix}$ (49)

The eigenvalues of d, which correspond to the canonical rates, are found to be
$1+\gamma,1,-\tanh(t)$. The last eigenvalue obviously entails ENM-ity. Note
that this example can be readily extended by using the GAD channel instead of
the AD channel.

Ref. [45] points out an example of non-unital ENM channel belonging to the
class of phase covariant qubit channels $\mathcal{E}_{\rm pc}$, which are
characterized by the property that $e^{i\mathcal{E}\sigma_{z}}\mathcal{E}_{\rm
pc}[\rho]e^{-i\mathcal{E}\sigma_{z}}=\mathcal{E}_{\rm
pc}[e^{i\mathcal{E}\sigma_{z}}\rho e^{-i\mathcal{E}\sigma_{z}}]$. The
generator has the form $\mathcal{L}_{\rm
pc}\equiv\gamma_{z}(t)[\sigma_{z}\rho\sigma_{z}-\rho]+\gamma(t)[\sigma_{+}\rho\sigma_{-}-\frac{1}{2}\\{\sigma_{-}\sigma_{+},\rho\\}]+\gamma(t)[\sigma_{-}\rho\sigma_{+}-\frac{1}{2}\\{\sigma_{+}\sigma_{-},\rho\\}],$
where it can be arranged so that $\gamma_{z}<0$ for all $t>0$ [45].
Incidentally, this model was also considered in Ref. [46] in a different
context.

Although this settles the question regarding the existence of a purely non-
unital ENM channel, this also leads to the observation that in the above
mentioned instances, the channel is not purely non-unital, in that the
generator in either case can be expressed as a mixture that includes a unital
component. Moreover, it is this unital component that introduces the ENM-ity.

It could be asked whether this no-go result is specific to two-level systems.
It is straightforward to show that these results hold in general to arbitrary
finite dimensions, by taking a trivial extension of the qubit dynamics, by
embedding the qubit in two dimensions and padding the additional dimensions
with an identity operation. This point can be manifestly illustrated by
considering another simple way to produce an example of non-unital ENM-ity,
namely by combining the Pauli ENM channel with non-unital channels taking
recourse to a higher dimensional space. Consider a ququart (4-level quantum
system), wherein the “lower” two-dimensional subspace, spanned by $\ket{0}$
and $\ket{1}$, is subject to the above Pauli ENM channel, whilst its “higher”
two-dimensional subspace, spanned by $\ket{2}$ and $\ket{3}$, is subject to an
AD or GAD channel. In other words, the two parts of the ququart are locally
decohered, so that the sum of diagonals is preserved within each subspace and
there are no cross-terms in the noise operator across the two subspaces (e.g.,
$\ket{1}\bra{2}$, etc.). This is mathematically captured by requiring that the
generator for the ququart evolution is given by:

$\mathcal{L}_{4}=\mathcal{L}_{01}(t)\oplus\mathcal{L}_{23}(t),$ (50)

where the subscripts to the symbol $\mathcal{L}$ denote the subspaces acted
upon. Evidently, the ququart as a whole is subject to non-unital noise, which
is ENM by virtue of the evolution in first qubit subspace. The resulting map
is given is given by

$\displaystyle e^{\int_{0}^{t}\mathcal{L}_{4}(s)ds}[\rho]$
$\displaystyle=\begin{pmatrix}\mathcal{E}_{01}&\textbf{0}\\\
\textbf{0}&\mathcal{E}_{23}\\\ \end{pmatrix}[\rho]$
$\displaystyle=\begin{pmatrix}\mathcal{E}_{01}&\textbf{0}\\\
\textbf{0}&\mathbb{I}_{2}\\\
\end{pmatrix}\cdot\begin{pmatrix}\mathbb{I}_{2}&\textbf{0}\\\
\textbf{0}&\mathcal{E}_{23}\\\ \end{pmatrix}[\rho]$
$\displaystyle=\sum_{j=1}^{3}\sum_{l=1,2}M_{j}N_{l}\rho
N_{l}^{\dagger}M_{j}^{\dagger}$ (51)

where 0 represents a $2\times 2$ null matrix. Here the Kraus operators $M_{j}$
for the Pauli channel are obtained by padding an identity operator:
$M_{j}=(E_{j}\oplus\mathbb{I}_{2})$, where $E_{j}$ are given by Eq. (7), In
particular, $M_{1}=\sqrt{2k}\mathbb{I}_{4}$, and

$M_{2}=\begin{pmatrix}0&\sqrt{k}&0&0\\\ \sqrt{k}&0&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{pmatrix};\leavevmode\nobreak\ \leavevmode\nobreak\
M_{3}=\begin{pmatrix}0&-i\sqrt{k}&0&0\\\ i\sqrt{k}&0&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{pmatrix};$

The Kraus operators for the AD channel similarly are given by:
$N_{k}\propto(\mathbb{I}_{2}\oplus K_{k})$, where $K_{k}$ are the qubit Kraus
operators given in Eq. (27), with $p{:=}1$.

$N_{1}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&\sqrt{1-\lambda}\end{pmatrix};\leavevmode\nobreak\ \leavevmode\nobreak\
N_{2}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&0&\sqrt{\lambda}\\\
0&0&0&0\end{pmatrix}$

the ququart Kraus operators for the ENM and AD channels, respectively, where
$k$ is given by Eq. (8). These diverse examples examined in the Section
suggest that a purely non-unital channel (i.e., one lacking any unital
component) is impossible. Although we have only studied the qubit case, we
conjecture that it is true for non-unital channels in arbitrary finite
dimensions.

##  V Conclusions

The ENM Pauli channel [14] presented the rather counter-intuitive idea of a
quantum noisy channel that is recoherent for all $t>0$. Here we address the
question of the possibility of a non-unital channel that is ENM, and show that
it is impossible for the well-known class of GAD channels of qubits.

We point out that this impossibility of eternal non-Markovianity does not hold
for arbitrary non-unital channels. Specifically, as evident from examples
discussed above, one can construct an ENM non-unital channel by a linear
combination of the generators of an ENM unital channel and an arbitrary non-
unital one. However, we conjecture that for a non-unital channel with no such
unital component (in the above sense) to be ENM, is impossible.

A consequence of our above negative result is that the best that a qubit GAD
channel can approach the ENM property is to be non-Markovian for all times
above a non-zero threshold time $t^{\ast}$, i.e., strictly $t^{\ast}>0$. By
way of illustration, we provided an example of a “quasi-ENM” GAD channel,
which is non-Markovian for all time $t\geq t^{\ast}>0$.

## Acknowledgments

SU was supported by IIT Madras via the Institute Postdoctoral Fellowship. R.S.
acknowledges the support of Department of Science and Technology (DST), India,
Grant No.: MTR/2019/001516. RS and SB acknowledge support from
Interdisciplinary Cyber Physical Systems (ICPS) programme of the Department of
Science and Technology (DST), India, Grant No.: DST/ICPS/QuST/Theme-1/2019/6
and DST/ICPS/QuST/Theme-1/2019/13. SB also acknowledges support from the
Interdisciplinary Program (IDRP) on Quantum Information and Computation (QIC)
at IIT Jodhpur.

## References

  * [1] Heinz-Peter Breuer and Francesco Petruccione.  The theory of open quantum systems.  Oxford University Press, 2002.
  * [2] Michael A Nielsen and Isaac Chuang.  Quantum computation and quantum information.  Cambridge University, 2000.
  * [3] Subhashish Banerjee.  Open quantum systems.  Springer Nature Singapore Pte Ltd., 2018.
  * [4] Angel Rivas, Susana F Huelga, and Martin B Plenio.  Quantum non-markovianity: characterization, quantification and detection.  Rep. Prog. Phys, 77(9):094001, 2014.
  * [5] Heinz-Peter Breuer, Elsi-Mari Laine, Jyrki Piilo, and Bassano Vacchini.  Colloquium: Non-markovian dynamics in open quantum systems.  Rev. Mod. Phys, 88(2):021002, 2016.
  * [6] Inés de Vega and Daniel Alonso.  Dynamics of non-markovian open quantum systems.  Rev. Mod. Phys., 89:015001, Jan 2017.
  * [7] U. Shrikant and Prabha Mandayam.  Quantum non-markovianity: Overview and recent developments.  Frontiers in Quantum Science and Technology, 2, 2023.
  * [8] N Pradeep Kumar, Subhashish Banerjee, R Srikanth, Vinayak Jagadish, and Francesco Petruccione.  Non-markovian evolution: a quantum walk perspective.  Open Systems & Information Dynamics, 25(03):1850014, 2018.
  * [9] Javid Naikoo, Subhashish Banerjee, and CM Chandrashekar.  Non-markovian channel from the reduced dynamics of a coin in a quantum walk.  Physical Review A, 102(6):062209, 2020.
  * [10] KG Paulson, Ekta Panwar, Subhashish Banerjee, and R Srikanth.  Hierarchy of quantum correlations under non-markovian dynamics.  Quantum Information Processing, 20(4):1–26, 2021.
  * [11] Heinz-Peter Breuer, Elsi-Mari Laine, and Jyrki Piilo.  Measure for the degree of non-markovian behavior of quantum processes in open systems.  Phys. Rev. Lett, 103(21):210401, 2009.
  * [12] AK Rajagopal, AR Usha Devi, and RW Rendell.  Kraus representation of quantum evolution and fidelity as manifestations of markovian and non-markovian forms.  Phys. Rev. A, 82(4):042107, 2010.
  * [13] Ángel Rivas, Susana F Huelga, and Martin B Plenio.  Entanglement and non-markovianity of quantum evolutions.  Phys. Rev. Lett, 105(5):050403, 2010.
  * [14] Michael J. W. Hall, James D. Cresser, Li Li, and Erika Andersson.  Canonical form of master equations and characterization of non-markovianity.  Phys. Rev. A, 89:042120, Apr 2014.
  * [15] Shin-Liang Chen, Neill Lambert, Che-Ming Li, Adam Miranowicz, Yueh-Nan Chen, and Franco Nori.  Quantifying non-markovianity with temporal steering.  Physical review letters, 116(2):020503, 2016.
  * [16] Shrikant Utagi.  Quantum causal correlations and non-markovianity of quantum evolution.  Physics Letters A, 386:126983, 2021.
  * [17] Shrikant Utagi, R Srikanth, and Subhashish Banerjee.  Temporal self-similarity of quantum dynamical maps as a concept of memorylessness.  Scientific Reports, 10(1):1–10, 2020.
  * [18] Li Li, Michael J.W. Hall, and Howard M. Wiseman.  Concepts of quantum non-markovianity: A hierarchy.  Physics Reports, 759:1 – 51, 2018.
  * [19] Felix A. Pollock, César Rodríguez-Rosario, Thomas Frauenheim, Mauro Paternostro, and Kavan Modi.  Non-markovian quantum processes: Complete framework and efficient characterization.  Phys. Rev. A, 97:012127, Jan 2018.
  * [20] Felix A. Pollock, César Rodríguez-Rosario, Thomas Frauenheim, Mauro Paternostro, and Kavan Modi.  Operational markov condition for quantum processes.  Phys. Rev. Lett., 120:040405, Jan 2018.
  * [21] Adrián A. Budini.  Quantum non-markovian processes break conditional past-future independence.  Phys. Rev. Lett., 121:240401, Dec 2018.
  * [22] Simon Milz, MS Kim, Felix A Pollock, and Kavan Modi.  Completely positive divisibility does not mean markovianity.  Physical review letters, 123(4):040401, 2019.
  * [23] Adrián A. Budini.  Conditional past-future correlation induced by non-markovian dephasing reservoirs.  Phys. Rev. A, 99:052125, May 2019.
  * [24] G. Lindblad.  Commun. Math. Phys., 48:119, 1976.
  * [25] E. C. G. Sudarshan, P. M. Mathews, and Jayaseetha Rau.  Stochastic Dynamics of Quantum-Mechanical Systems.  Phys. Rev., 121(3):920–924, February 1961.
  * [26] Philip Pechukas.  Reduced dynamics need not be completely positive.  Physical review letters, 73(8):1060, 1994.
  * [27] Robert Alicki.  Comment on “reduced dynamics need not be completely positive”.  Physical review letters, 75(16):3020, 1995.
  * [28] Philip Pechukas.  Pechukas replies.  Physical review letters, 75(16):3021, 1995.
  * [29] Man-Duen Choi.  Completely positive linear maps on complex matrices.  Linear algebra and its applications, 10(3):285–290, 1975.
  * [30] E. C. G. Sudarshan, P. M. Mathews, and Jayaseetha Rau.  Stochastic dynamics of quantum-mechanical systems.  Phys. Rev., 121:920–924, Feb 1961.
  * [31] Dariusz Chruściński, Ángel Rivas, and Erling Størmer.  Divisibility and information flow notions of quantum markovianity for noninvertible dynamical maps.  Physical review letters, 121(8):080407, 2018.
  * [32] Dariusz Chruściński and Filip A Wudarski.  Non-markovian random unitary qubit dynamics.  Physics Letters A, 377(21-22):1425–1429, 2013.
  * [33] Dariusz Chruściński and Filip A Wudarski.  Non-markovianity degree for random unitary evolution.  Physical Review A, 91(1):012104, 2015.
  * [34] Filip A Wudarski and Dariusz Chruściński.  Markovian semigroup from non-markovian evolutions.  Physical Review A, 93(4):042120, 2016.
  * [35] Michael Marc Wolf, J Eisert, TS Cubitt, and J Ignacio Cirac.  Assessing non-markovian quantum dynamics.  Physical review letters, 101(15):150402, 2008.
  * [36] Nina Megier, Dariusz Chruściński, Jyrki Piilo, and Walter T Strunz.  Eternal non-markovianity: from random unitary to markov chain realisations.  Scientific reports, 7(1):1–11, 2017.
  * [37] Katarzyna Siudzińska and Dariusz Chruściński.  Quantum evolution with a large number of negative decoherence rates.  Journal of Physics A: Mathematical and Theoretical, 53(37):375305, 2020.
  * [38] R Srikanth and Subhashish Banerjee.  Squeezed generalized amplitude damping channel.  Phys. Rev. A, 77(1):012318, 2008.
  * [39] Subhashish Banerjee, R Srikanth, CM Chandrashekar, and Pranaw Rungta.  Symmetry-noise interplay in a quantum walk on an n-cycle.  Phys. Rev. A, 78(5):052316, 2008.
  * [40] Arijit Dutta, Junghee Ryu, Wiesław Laskowski, and Marek Żukowski.  Entanglement criteria for noise resistance of two-qudit states.  Physics Letters A, 380(27-28):2191–2199, 2016.
  * [41] Alejandro Fonseca.  High-dimensional quantum teleportation under noisy environments.  Physical Review A, 100(6):062311, 2019.
  * [42] Salvatore Lorenzo, Francesco Plastina, and Mauro Paternostro.  Geometrical characterization of non-markovianity.  Physical Review A, 88(2):020102, 2013.
  * [43] Jing Liu, Xiao-Ming Lu, and Xiaoguang Wang.  Nonunital non-markovianity of quantum dynamics.  Phys. Rev. A, 87:042103, Apr 2013.
  * [44] Nina Megier, Andrea Smirne, and Bassano Vacchini.  Entropic bounds on information backflow.  Physical Review Letters, 127(3):030401, 2021.
  * [45] Sergey N Filippov, AN Glinov, and Leevi Leppäjärvi.  Phase covariant qubit dynamics and divisibility.  Lobachevskii Journal of Mathematics, 41(4):617–630, 2020.
  * [46] Bassano Vacchini, Andrea Smirne, Elsi-Mari Laine, Jyrki Piilo, and Heinz-Peter Breuer.  Markovianity and non-markovianity in quantum and classical systems.  New Journal of Physics, 13(9):093004, 2011.


# Multi-Agent Terraforming:
Efficient Multi-Agent Path Finding via Environment Manipulation

David Vainshtein, Kiril Solovey, Oren Salzman

###### Abstract

Multi-agent pathfinding (MAPF) is concerned with planning collision-free paths
for a team of agents from their start to goal locations in an environment
cluttered with obstacles. Typical approaches for MAPF consider the locations
of obstacles as being _fixed_ , which limits their effectiveness in automated
warehouses, where obstacles (representing pods or shelves) can be _moved_ out
of the way by agents (representing robots) to relieve bottlenecks and
introduce shorter routes. In this work we initiate the study of MAPF with
movable obstacles. In particular, we introduce a new extension of MAPF, which
we call Terraforming MAPF (tMAPF), where some agents are responsible for
moving obstacles to clear the way for other agents. Solving tMAPF is extremely
challenging as it requires reasoning not only about collisions between agents,
but also _where_ and _when_ obstacles should be moved. We present extensions
of two state-of-the-art algorithms, CBS and PBS, in order to tackle tMAPF, and
demonstrate that they can consistently outperform the best solution possible
under a static-obstacle setting.

##  1 Introduction

> The impediment to action advances action.
> What stands in the way becomes the way.
>
> Marcus Aurelius

Multi-agent path finding (MAPF) is a popular algorithmic framework that
captures complex tasks involving mobile agents that need to plan individual
routes while avoiding collisions during plan execution (Stern et al. 2019;
Salzman and Stern 2020). This abstraction has been successfully applied to a
variety of settings (see, e.g., Wurman, D’Andrea, and Mountz (2008); Belov et
al. (2020); Li et al. (2021a); Choudhury et al. (2021)). However, in some
cases this formulation may not be expressive enough to fully capture the
underlying task, which can lead to suboptimal performance.

This is especially true in the context of automated warehouses where we are
given a stream of tasks, and the goal is to maximize the system’s throughput.
In this setting, formulated as a _lifelong_ MAPF (L-MAPF) problem and
typically solved via a sequence of MAPF queries (Ma et al. 2017; Liu et al.
2019; Li et al. 2021b), inventory pods that hold goods are manipulated by a
large team of mobile agents (or robots): agents pick up pods, carry them to
designated dropoff locations where goods are manually removed from the pods
(to be packaged for customers); each pod is then carried back by a robot to a
(possibly different) storage location (Wurman, D’Andrea, and Mountz 2008).
When applied to this setting, existing variants of MAPF and L-MAPF tend to
impose the following limiting and artificial constraint: pods that are not
currently carried to or from a dropoff location are modeled as _static
obstacles_ , which cannot be moved. Thus, those approaches overlook the fact
that pods can be manipulated to clear the way for agents and reduce the travel
time or distance of agents in the system.

(a) MAPF problem

(b) Static obstacles

(c) Terraforming

Figure 1:  Comparing MAPF and tMAPF for a toy problem with agents
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}a^{t}_{1}},{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}a^{t}_{2}},{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}a_{3}}$
(circles), and a row of obstacles (grey squares). (1(a)) MAPF problem
assigning agents to their goal locations (squares). Here, agent
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}a_{3}}$ is
a task agent (denoted as
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}a_{3}^{t}}$)
with no task assigned to him. (1(b)) MAPF solution where agents must avoid
collisions with obstacles and each other. Here, agent
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}a_{3}}$
cannot help the other two agents as all obstacles are static. (1(c)) tMAPF
solution where
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}a_{3}}$ is
a mover agent (denoted as
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}a_{1}^{m}}$)
that creates a shortcut by clearing a movable obstacle $o_{3}^{m}$. Video
https://bit.ly/3ImgfAw

To bridge the gap between existing MAPF formulations and the type of problems
they are intended to tackle in the real world, we explore the implications of
allowing agents the extra flexibility of _manipulating the environment by
moving obstacles_ (e.g., dynamically relocating pod’s locations in warehouse
applications). To this end, we introduce a new MAPF variant which we term
_“Terraforming MAPF”_ (or tMAPF in short)111 The term “Terraforming”, which
originated from SciFi literature and was recently used in the context of space
exploration, is the process of deliberately altering the environment of a
planet to make it habitable. . In tMAPF, formally defined in Sec. 3, the input
consist of the agents’ start and goal location and a set of obstacles (as in
MAPF), as well as a specification of the obstacles whose position can be
changed while answering the tMAPF query (see Fig. 1 and details below).

We suggest two algorithms to solve the tMAPF problem which are described in
Section 5. The first, which we call TF-CBS, is based on the celebrated CBS
algorithm (Sharon et al. 2015) which we review in Sec. 4. TF-CBS is complete
and is guaranteed to produce a cost-optimal path. The second algorithm, which
we call TF-PBS, offers computational efficiency by trading completeness and
optimality guarantees in favor of rapidly attaining high-quality solutions. In
our evaluation, described in Sec. 6, we demonstrate how both algorithms
consistently out-perform the optimal solution produced by classical MAPF
algorithms on several different metrics.

In this work we concentrate on the algorithmic implication of terraforming, or
environment manipulation, in the context of MAPF. However, our ultimate goal
is to apply terraforming to the L-MAPF problem. As we will see, this is
extremely challenging which is why in this work we limit ourselves to a
simplified tMAPF setting. We discuss the research challenges of moving from
terraforming in the context of MAPF to L-MAPF in Sec. 7.

##  2 Related work

A variety of approaches were developed to solve the MAPF problem and its many
variants using algorithmic tools such as network flow (Yu and LaValle 2012),
satisfiability (Surynek et al. 2016), Answer Set Programming (Erdem et al.
2013) and search-based methods (Barer et al. 2014; Boyarski et al. 2015;
Sharon et al. 2015; Li et al. 2019).

In this work we adapt approaches from the latter group. Of specific relevance
to our work are Conflict-Based Search (CBS) (Sharon et al. 2015) which is used
as an algorithmic building block in many state-of-the-art MAPF solvers (see,
e.g. Greshler et al. (2021)) and Priority-Based Search (PBS) (Ma et al. 2019)
which is commonly used when solving the L-MAPF problem (Ma et al. 2017; Li et
al. 2021b). As both are used to develop our new algorithms, we elaborate on
both algorithms in Sec. 4.

Arguably, the most closely-related work to our new problem formulation is by
Bellusci, Basilico, and Amigoni (2020) who introduce a new variant of the MAPF
problem, called Configurable MAPF (C-MAPF). Here, the structure of the
environment is configurable (given a set of constraints) and the problem calls
for computing an environment configuration as well as a solution to a MAPF
query. Note that a solution to the C-MAPF problem is a _fixed_ environment (as
well as a set of paths dictating where each agent should go). In contrast, in
tMAPF the environment _dynamically changes_ as agents execute their paths.

Finally, we mention the works of Hauser (2013, 2014) that focus on computing
paths that minimize the number of obstacles to remove from a given
environment. However, these works are formulated for a single-agent and even
then, the agent is not tasked with removing said obstacles. Also related is
the work on Multi-Robot Clutter Removal (MRCR) (Tang, Han, and Yu 2020). Here,
a group of agents are tasked with clearing obstacles from a given environment
while avoiding collisions with each other and any remaining obstacles.

##  3 Problem Formulation

In this section we provide a formal definition of the tMAPF problem. In
preparation, we first formally define the (classical) MAPF problem.

###  3.1 Multi-Agent Path Finding (MAPF)

A MAPF problem is a tuple

$\langle\mathcal{G},\mathcal{A},\mathcal{O},\mathcal{V}_{\text{start}},\mathcal{V}_{\text{goal}}\rangle,$

where $\mathcal{G}=(\mathcal{V},\mathcal{E})$ is the environment graph,
$\mathcal{A}=\\{a_{1},\ldots,a_{n}\\}$ is the agent set,
$\mathcal{O}=\\{o_{1},\ldots,o_{\ell}\\}\subset\mathcal{V}$ is the obstacle
set, $\mathcal{V}_{\text{start}}=\\{s_{1},\ldots,s_{n}\\}\subset\mathcal{V}$
is the set of agents’ start vertices and
$\mathcal{V}_{\text{goal}}=\\{g_{1},\ldots,g_{n}\\}\subset\mathcal{V}$ is the
set of agents’ goal vertices. We now elaborate on those ingredients.

We assume that the graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ on which the
agents $\mathcal{A}$ operate is undirected and includes self edges to all the
vertices to simulate agent wait actions, i.e., $(v,v)\in\mathcal{E}$ for every
vertex $v\in\mathcal{V}$ and when an agent moves along such an edge we will
say that it waits in place. In typical MAPF formulations, _static obstacles_ ,
which block certain agent positions, are implicitly encoded via the graph
$\mathcal{G}$, where vertices representing static obstacle locations are
removed. In our setting, in preparation to the tMAPF problem where some agents
can share locations with obstacles, it will be convenient to explicitly
account for obstacles. In particular, the set of obstacles $\mathcal{O}$
denotes graph vertices that are blocked, and which the agents cannot visit.

#### States.

A MAPF _state_ $S=(v_{1},\ldots,v_{m})$ is a vector where
$v_{i}\in\mathcal{V}$ represents the location of agent $a_{i}$. We say that
state $S$ is _valid_ if the following conditions are met:

  * S1

No two agents share the same location, i.e., $v_{i}\neq v_{j}$ for any two
agents $a_{i}\neq a_{j}$.

  * S2

Agents do not collide with obstacles, i.e., $v_{i}\not\in\mathcal{O}$ for any
agent $a_{i}$.

#### Transitions.

A _transition_ between two valid states $S=(v_{1},\ldots,v_{m})$ and
$S^{\prime}=(v^{\prime}_{1},\ldots,v^{\prime}_{m})$ is _valid_ if the
following conditions are met:

  * T1

Agents move along edges, i.e., $\forall
a_{i}\in\mathcal{A},~{}(v_{i},v^{\prime}_{i})\in\mathcal{E}$.

  * T2

Agents do not swap locations over the same edge, i.e., $\forall
a_{i},a_{j}\in\mathcal{A}~{}s.t.~{}i\neq j$ it holds that $v_{i}\neq
v^{\prime}_{j}$ or $v^{\prime}_{i}\neq v_{j}$.

#### Solution.

A _solution_ to the above MAPF problem is a sequence of states
$\pi=(S^{1},\ldots,S^{k})$, such that each state
$S^{i}:=(v_{1}^{i},\ldots,v_{n}^{i})$ is valid, each transition from $S^{i}$
to $S^{i+1}$ is valid for any $i$, $S^{1}=\mathcal{V}_{\text{start}}$ and
$S^{k}=\mathcal{V}_{\text{goal}}$.

#### Solution cost in MAPF.

To define the cost of a solution $\pi$ we first define a
$\text{cost}(\pi,a_{i})$ for agent $a_{i}\in\mathcal{A}$ to be the earliest
arrival time to its goal after which $a_{i}$ does not change its location.
Namely, $\text{cost}(\pi,a_{i})$ is the smallest $j$ s.t.
$\forall\tau\in[j,k],v^{\tau}=g_{i}$. The cost of a solution is then defined
as $\text{Cost}(\pi):=\sum_{i}\text{cost}(\pi,a_{i})$ and is commonly referred
to as the _sum of costs_ (we explain in Sec. 3.3 why we use this cost function
and not other commonly-used ones such as the makespan (Stern et al. 2019)).

###  3.2 Terraforming MAPF

A tMAPF problem, which generalizes the MAPF problem (and is visualized in Fig.
1), is a tuple

$\langle\mathcal{G},\mathcal{A}=\mathcal{A}^{t}\cup\mathcal{A}^{m},\mathcal{O}=\mathcal{O}^{s}\cup\mathcal{O}^{m},\\\
\mathcal{V}_{\text{start}}=\mathcal{V}_{\text{start}}^{t}\cup\mathcal{V}_{\text{start}}^{ma}\cup\mathcal{V}_{\text{start}}^{mo},\mathcal{V}_{\text{goal}}^{t}\rangle.$

Here, $\mathcal{G}=(\mathcal{V},\mathcal{E})$ is a graph defined as in the
MAPF setting. However, now there are two types of agents: $n_{t}$ _task
agents_ $\mathcal{A}^{t}$, and $n_{m}$ _mover agents_ $\mathcal{A}^{m}$.
Similar to the MAPF setting, task agents
$\mathcal{A}^{t}=\\{a_{1}^{t},\ldots,a^{t}_{n_{t}}\\}$ are forbidden from
moving to either of the two types of obstacle vertices (see below) and have
designated start vertices
$\mathcal{V}_{\text{start}}^{t}=\\{s_{1}^{t},\ldots,s^{t}_{n_{t}}\\}\subset\mathcal{V}$
and goal vertices
$\mathcal{V}_{\text{goal}}^{t}=\\{g_{1}^{t},\ldots,g^{t}_{n_{t}}\\}\subset\mathcal{V}$.
In contrast, mover agents
$\mathcal{A}^{m}=\\{a_{1}^{m},\ldots,a^{m}_{n_{m}}\\}$ can share locations
with obstacles (to simulate robots going underneath pods in automated
warehouses), and only have designated start vertices
$\mathcal{V}_{\text{start}}^{ma}=\\{s_{1}^{ma},\ldots,s^{ma}_{n_{m}}\\}\subset\mathcal{V}$
without predefined goal vertices.222As a convention, superscript ‘t’
corresponds to task agents, superscript ‘s’ corresponds to static obstacles
while superscript ‘m’ corresponds to either mover agents or movable obstacles.
When we need to distinguish between the latter two we use superscripts ‘ma’
and ‘mo’ to refer to mover agents and movable obstacles, respectively.

In tMAPF we have two types of obstacles: $\ell_{s}$ _static obstacles_
$\mathcal{O}^{s}$, and $\ell_{m}$ _movable obstacles_ $\mathcal{O}^{m}$ whose
location can change via a mover agents (to be explained shortly). Static
obstacles $\mathcal{O}^{s}=\\{o_{1}^{s},\ldots,o^{s}_{\ell_{s}}\\}$ are
associated with their vertex locations, i.e.,
$\mathcal{O}^{s}\subset\mathcal{V}$, as in the MAPF setting. Every movable
obstacles $o^{m}_{i}\in\mathcal{O}^{m}$ has a unique start position
$s_{1}^{mo}\in\mathcal{V}_{\text{start}}^{mo}\subset\mathcal{V}$ that also
serves as its goal position. The location of a movable obstacle can change
throughout the execution via the mover agents.

#### States.

A tMAPF state encodes the locations of the two types of agents, as well as the
locations of the movable obstacles. In particular, a state
$S=(v_{1},\ldots,v_{n_{t}},u_{1},\ldots,u_{n_{m}},w_{1},\ldots,w_{\ell_{m}})$
is a vector of vertices, where $v_{i}\in V$ represents the location of task
agent $a^{t}_{i}\in\mathcal{A}^{t}$, $u_{j}\in V$ represents the location of a
mover agent $a^{m}_{j}\in\mathcal{A}^{m}$, and $w_{k}\in V$ represents the
location of the movable obstacles $o^{m}_{k}\in\mathcal{O}^{m}$. A state $S$
is _valid_ if the following conditions are met:

  * S1’

No two agents share the same vertex, i.e., $v_{i}\neq v_{j}$ for any
$a^{t}_{i}\neq a^{t}_{j}$, $u_{i}\neq u_{j}$ for any $a^{m}_{i}\neq
a^{m}_{j}$, and $v_{i}\neq u_{j}$ for any
$a^{t}_{i}\in\mathcal{A},a^{m}_{j}\in\mathcal{A}^{m}$.

  * S2’

No two obstacles share the same vertex, i.e., $w_{i}\neq w_{j}$ and
$w_{i},w_{j}\not\in\mathcal{O}^{s}$ for any $o^{m}_{i}\neq
o^{m}_{j}\in\mathcal{O}^{m}$.

  * S3’

No task agents collide with obstacles, i.e., $v_{i}\not\in\mathcal{O}^{m}$ and
$v_{i}\neq w_{j}$ for any $a_{i}\in\mathcal{A}$ and
$o^{m}_{j}\in\mathcal{O}^{m}$.

#### Transitions.

A _transition_ between two valid states
$S=(v_{1},\ldots,v_{n_{t}},u_{1},\ldots,u_{n_{m}},w_{1},\ldots,w_{\ell_{m}})$
and
$S^{\prime}=(v^{\prime}_{1},\ldots,v^{\prime}_{n_{t}},u^{\prime}_{1},\ldots,u^{\prime}_{n_{m}},w^{\prime}_{1},\ldots,w^{\prime}_{\ell_{m}})$
is _valid_ if the following conditions are met:

  * T1’

Agents and movable obstacles move along graph edges, i.e., $\forall
i~{}(S[i],S^{\prime}[i])\in\mathcal{E}$ where $S[i]$ ($S^{\prime}[i]$) denotes
the $i$’th element of $S$ ($S^{\prime}$).

  * T2’

Agents and movable obstacles do not swap locations using the same edge, i.e.,
$\forall i\neq j~{}S[i]\neq S^{\prime}[j]$ or $S[j]\neq S^{\prime}[i]$.

  * T3’

Movable obstacles can only move via a mover agent, i.e., if $w_{i}\neq
w^{\prime}_{i}$, for some movable obstacle $o^{m}_{i}\in\mathcal{O}^{m}$, then
there exists a mover $a^{m}_{j}\in\mathcal{A}^{m}$ such that $w_{i}=u_{j}$ and
$w^{\prime}_{i}=u^{\prime}_{j}$.

#### Solution.

A _solution_ to the above tMAPF problem is a sequence of states
$\pi=(S^{1},\ldots,S^{k})$, such that each state
$S^{i}:=(v_{1}^{i},\ldots,v_{n_{t}}^{i},u_{1}^{i},\ldots,u_{n_{m}}^{i},w_{1}^{i},\ldots,w_{\ell_{m}}^{i})$
is valid, each transition from $S^{i}$ to $S^{i+1}$ is valid for any $i$,
$S^{1}=(s^{t}_{1},\ldots,s^{t}_{n_{t}},s^{ma}_{1},\ldots,s^{ma}_{n_{m}},s^{mo}_{1},\ldots,s^{mo}_{\ell_{m}})$,
$v^{k}_{i}=g^{t}_{i}$ for any task agent $a^{t}_{i}\in\mathcal{A}^{t}$ and
$w^{k}_{i}=s^{mo}_{i}$ for any movable obstacle $o^{m}_{i}\in\mathcal{O}^{m}$.

#### Simplifying assumptions.

In this work we impose several simplifying assumptions on our tMAPF problem:

  * A1

A mover agent can move at most one movable obstacle.

  * A2

The number of movable obstacles is equal to the number of mover agents.

In Sec. 7 we discuss how these simplifying assumptions can be lifted in future
work.

#### Solution cost in tMAPF.

We consider several cost metrics for a given solution
$\pi=(S^{1},\ldots,S^{k})$. In preparation, we define different cost functions
for the different agents and for the movable obstacles. Similar to the MAPF
setting, for a task agent $a^{t}_{i}\in\mathcal{A}^{t}$, we define
$\text{cost}^{t}(\pi,a^{t}_{i})$ to be the smallest $j$ s.t.
$\forall\tau\in[j,k],v_{i}^{\tau}=g_{i}$. We define the total task-agent cost
as $\text{Cost}^{t}(\pi)=\sum_{i}\text{cost}^{t}(\pi,a^{t}_{i})$.

For a mover agent $a^{m}_{i}\in\mathcal{A}^{m}$, denote
$\text{cost}^{ma}_{\bcancel{w}}(\pi,a^{m}_{i})$ to be the number of steps to
reach the movable obstacle it is going to move without accounting for wait
actions, i.e., the cost of a wait action is zero. The total cost for mover
agents to reach the movable obstacles is
$\text{Cost}^{ma}_{\bcancel{w}}(\pi)=\sum_{i}\text{cost}^{ma}_{\bcancel{w}}(\pi,a^{m}_{i})$.
Additionally, define $\text{cost}^{mo}_{\bcancel{w}}(\pi,a^{m}_{i})$ to be the
number of steps a movable obstacle takes without wait actions and define the
total cost for movable obstacles to be
$\text{Cost}^{mo}_{\bcancel{w}}(\pi)=\sum_{i}\text{cost}^{mo}_{\bcancel{w}}(\pi,a^{m}_{i})$.

In this work we will consider the following two cost functions accounting for
motions taken by both types of agents as well as the moveable obstacles.

$\text{Cost}_{1}(\pi)=\text{Cost}^{t}(\pi)+\text{Cost}^{mo}_{\bcancel{w}}(\pi).$
(1)
$\text{Cost}_{2}(\pi)=\text{Cost}^{t}(\pi)+\text{Cost}^{ma}_{\bcancel{w}}(\pi)+\text{Cost}^{mo}_{\bcancel{w}}(\pi).$
(2)

In both cost functions we treat task agents just as in the MAPF setting and
the difference arises in how we treat mover agents and movable obstacles. In
addition, in both cost functions we do not account for wait actions (see
discussion in Sec. 3.3) incurred by the mover agents and movable obstacles.
Intuitively, in Eq. (1) we ignore the cost of reaching a movable obstacle by a
mover agent and only account for the “work” required to move movable obstacle.
In Eq. (2) we add the motions required by a mover agent.

###  3.3 Discussion

When formalizing the tMAPF problem, there are many subtle-yet-important
variants one can consider. E.g., “what cost function to use?” and “are mover
agents allowed to move under static obstacles?”. In our formulation, we made
sure to have our tMAPF formulation a _generalization_ of the standard MAPF
formulation while at the same time serve as a stepping stone to our ultimate
goal of applying terraforming to the L-MAPF problem. Indeed, in Sec. 7 we
discuss what are the steps and challenges required to reach this goal.

Specifically, in the MAPF setting we chose to use the sum of costs as our cost
function as it extends naturally to the lifelong setting where we wish to
maximize a system’s throughput. Moving to the terraforming version of L-MAPF
in warehouse applications, we envision that there will be no pre-allocation of
the agents to two distinct groups of task and mover agents. These two groups
will form naturally where agents en-route to a pick an item (obstacle in our
formulation) will serve as mover agents while agents already carrying said
items will serve as task agents. This motivated us to (i) allow mover agents
to move under static obstacles and (ii) to propose cost functions that focus
on the work done to move movable obstacles and not necessarily to reach them
as we assume that this will be done by agents en-route to reaching a goal.
Thus, Eq. (1) may be seen as a _lower-bound_ on the cost to move obstacles as
we do not account for the steps taken to reach it and assume that the mover
agent would have passed next to the movable obstacle. Similarly, Eq. (2) may
be seen as an _upper-bound_ on the cost to move obstacles as we assume that at
least some of the steps taken to reach it would have been carried out
regardless.

##  4 Algorithmic Background

Before we present our approaches for tMAPF, we describe in this section two
algorithmic building blocks. Namely, the CBS and PBS algorithms for
(classical) MAPF.

###  4.1 Conflict-Based Search

Conflict-Based Search (CBS) is a popular approach for the (classical) MAPF
problem, which is both complete and optimal. We now provide an overview of
CBS. Specifically, we describe a recent variant by Li et al. (2019) that uses
both positive and negative constraints (to be explained shortly) as it was
shown to have better runtime both empirically and when looking at a worst-case
complexity analysis (Gordon, Filmus, and Salzman 2021). We refer the reader to
(Sharon et al. 2015) and (Li et al. 2019) for the full description.

CBS maintains constraints between agents which are used to resolve conflicts.
On the high-level search, CBS explores a _constraint tree_ (CT), where a given
node $N$ of CT encodes a set of constraints $C_{N}$ on the locations of agents
in time and space. In particular, CBS includes positive constraints denoted by
$(+,a_{i},v,\tau)$ which specify that an agent $a_{i}$ must visit vertex $v$
at time step $\tau$, as well as negative constraints denoted by
$(-,a_{i},v^{\prime},\tau^{\prime})$ which prohibits agent $a_{i}$ from
visiting vertex $v^{\prime}$ at time step $\tau^{\prime}$. Similar constraints
are imposed on edges. In addition to the constraints, each CT node maintains
single-agent paths that represent the current MAPF solution (possibly
containing conflicts) as well as the cost of the solution.

CBS starts the high-level search with the tree root whose constraint set is
empty, and assigns to each agent its shortest path, while avoiding static
obstacles but ignoring interactions between agents. Whenever CBS expands a
node $N$, it invokes a low-level search to compute a new set of paths that
abide by the constraint set $C_{N}$ (see details below). If a collision
between agents, e.g., $a_{i}$ and $a_{j}$, at a vertex $v$ (or an edge) at
time step $\tau$ is encountered in the new paths, CBS generates two child CT
nodes $N_{1},N_{2}$ with the updated constraints
$C_{N_{1}}=C_{N}\cup\\{(+,a_{i},v,\tau)\\},C_{N_{2}}=C_{N}\cup\\{(-,a_{i},v,\tau)\\}$,
respectively. That is, the node $N_{1}$ includes constraints to enforce that
$a_{i}$ visits $v$ at time $\tau$ which implicitly adds the constraints
$(-,a_{\ell},v,\tau)$ for $\ell\neq i$. The node $N_{2}$ encodes the opposite
situation with respect to $a_{i}$ forcing it to avoid $v$ at time $\tau$. The
high-level search chooses to expand at each iteration a CT node with the
lowest cost. The high-level search terminates when a valid solution is found
at some node $N$, or when no more nodes for expansion remain, in which case,
CBS declares failure.

The low-level search of CBS proceeds as follows: For a given CT node $N$, an
A* (Hart, Nilsson, and Raphael 1968) search is invoked for a particular agent
$a_{i}$ that violates the node’s constraints. Importantly, the search
simultaneously explores agent positions in time and space, and while doing so
ensures that the constraints $C_{N}$ are satisfied. Once a path for the agent
is found, the set of solution paths for $N$ is updated.

###  4.2 Priority-Based Search

Priority-based search (PBS) is a recent approach for MAPF. It forgoes the
completeness and optimality guarantees for the sake of computational
efficiency. We now provide an overview of PBS, and refer the reader to (Ma et
al. 2019) for the full description. At its core, PBS has some resemblance to
CBS, in the sense that it is a hierarchical approach with high and low level
search. However, unlike CBS which maintains space-time constraints between the
agents in the high-level search, PBS maintains priorities between agents.

On the high-level search, PBS explores a _priority tree_ (PT), where a given
node $N$ of PT encodes a (partial) priority set $P_{N}=\\{a_{h}\prec
a_{i},~{}a_{j}\prec a_{l}$, …}. A priority $a_{i}\prec a_{j}$ means that agent
$a_{i}$ has precedence over agent $a_{j}$ whenever a low-level search is
invoked (see below). In addition to the ordering, each PT node maintains
single-agent paths that represent the current MAPF solution (possibly
containing conflicts). PBS starts the high-level search with the tree root
whose priority set is empty, and assigns to each agent its shortest path.
Whenever PBS expands a node $N$, it invokes a low-level search to compute a
new set of paths which abide by the priority set $P_{N}$. If a collision
between agents, e.g., $a_{i}$ and $a_{j}$, is encountered in the new paths,
PBS generates two child PT nodes $N_{1},N_{2}$ with the updated priority sets
$P_{N_{1}}=P_{N}\cup\\{a_{i}\prec a_{j}\\},P_{N_{2}}=P_{N}\cup\\{a_{j}\prec
a_{i}\\}$, respectively. The high-level search chooses to expand at each
iteration a PT node in a depth-first search manner. The high-level search
terminates when a valid solution is a found at some node $N$, or when no more
nodes for expansion remain, in which case, PBS declares failure.

The low-level search of PBS proceeds in the following manner. For a given PT
node $N$, PBS performs a topological sort of the agents according to $P_{N}$
from high priority to low, and plans individual-agent paths based on the
ordering. For a given topological ordering
$(a^{\prime}_{1},\ldots,a^{\prime}_{k^{\prime}})\subset\mathcal{A}$, for some
$1\leq k^{\prime}\leq k$, the low-level iterates over the $k^{\prime}$ agents
in the topological ordering, and updates their paths such that they do not
collide with any higher-priority agents. (Note that agents that do not appear
on this list maintain their original plans.) It then checks whether collisions
occur between all the agents combined.

##  5 Algorithmic framework

In this section we present our algorithmic contributions for tackling the
tMAPF problem. We first describe a complete and optimal approach that is based
on CBS, and then proceed to a faster but incomplete PBS-based method.

###  5.1 A CBS-based approach for tMAPF

We present an extension of the CBS algorithm called TF-CBS, to solve the tMAPF
problem. The main idea behind TF-CBS is to associate each mover agent
$a^{m}_{i}\in\mathcal{A}^{m}$ with a particular movable obstacle
$o^{m}_{j}\in\mathcal{O}^{m}$ and treat those two elements as one entity in
both levels of the search. In other words, $a^{m}_{i}$ can be thought of as
agent that has two physical interpretations, that of the location of the
actual agent, and the location of the obstacle $o^{m}_{j}$ it is required to
move. The assignment of a mover to a movable obstacle is done once in the
beginning of the execution of TF-CBS. In particular, we use a greedy heuristic
where every mover agent $a^{m}_{i}\in\mathcal{A}^{m}$ is assigned to the
closest movable obstacle $o^{m}_{j}\in\mathcal{O}^{m}$ according to the
distance over $\mathcal{G}$, while ignoring obstacle locations, from
$s^{ma}_{i}$ to $s^{mo}_{j}$. If an agent’s closest obstacle is already taken
by another agent, then the next-closest obstacle is assigned (and so on). We
leave the study of stronger assignment approaches, such as the Hungarian
method (Kuhn 1955), for future work. For the remainder of this section, we
assume that each mover agent $a^{m}_{i}\in\mathcal{A}^{m}$ is associated with
a specific movable obstacle that, w.l.o.g, is denote by $o^{m}_{i}$.

Given that we treat $a^{m}_{i}$ and $o^{m}_{i}$ as one entity, we only impose
in the high-level search of TF-CBS constraints on the agents of
$\mathcal{A}^{t}$ and $\mathcal{A}^{m}$ but not on the movable obstacles.
Next, we describe the constraints we use in TF-CBS, which generalize those of
CBS. We first consider constraints involving task agents. For a given
$a^{t}_{i}\in\mathcal{A}^{t}$ the positive and negative constraints
$(+,a^{t}_{i},v,t)$ and $(-,a^{t}_{i},v,t)$ are defined exactly as in the
setting of CBS, i.e., $a^{t}_{i}$ should or should not visit $v$ at time $t$,
respectively.

Next, consider a mover agent $a^{m}_{i}$ and its assigned obstacle
$o^{m}_{i}$. The positive constraint $(+,a^{m}_{i},v,t)$ should be interpreted
as the agent should visit $v$ at time $t$, with or without $o^{m}_{i}$. The
negative constraint $(-,a^{m}_{i},v,t)$ simply means that $a^{m}_{i}$ should
not be in $v$ at time $t$, with or without $o^{m}_{i}$.

Now we introduce a new type of constraint to prevent collisions with movable
obstacles before they are reached by the designated mover, which is not
covered via the previous constraints. In particular, consider a movable
obstacle $o^{m}_{i}$ and suppose that we have a lower bound $t$ on the time of
arrival of agent $a^{m}_{i}$ to $s^{mo}_{i}$. In such a case we would like to
inform task agents to avoid getting to $s^{mo}_{i}$ before time $t$, as they
surely cannot reach this location since it is blocked and cannot be moved
until after timestep $t$. Similarly, we would like to prevent mover agents
that currently carry obstacles from reaching this vertex as well. Thus, we
introduce the timed constraint $(s^{mo}_{i},t)$ requiring that any agent
reaching $s^{mo}_{i}$ before time $t$ is a mover agent that is not carrying an
obstacle while passing through this vertex.

Before wrapping up the description of the high-level search, we mention that
TF-CBS maintains for every CT node its cost computed either using Eq. (1) or
Eq. (2). TF-CBS determines the next CT node to expand in a best-first search
manner according to the node cost.

We now proceed to describe the low-level search of TF-CBS (visualized in Fig.
1). For a task agent $a^{t}_{i}\in\mathcal{A}^{t}$ the search proceeds in a
manner similar to CBS using A* and while adhering to the constraints of the
current high-level node. Note that during this process the agent needs to
avoid collisions with static obstacles, but collision avoidance with movable
obstacles is enforced through positive constraints over mover agents, e.g.,
$(+,a^{m}_{i},v,\tau)$, or timed constraints of the form $(s^{mo}_{i},\tau)$.

When a low-level search for an agent $a^{m}_{i}$ is invoked, its search
conceptually consists of two parts: (i) $a^{m}_{i}$ moves from its start
location $s^{ma}_{i}$ to the start location $s^{mo}_{i}$ of $o^{m}_{i}$, and
(ii) $a^{m}_{i}$ moves from $s^{mo}_{i}$ to $s^{mo}_{i}$, while carrying the
obstacle $o^{m}_{j}$. Note that in part (i) the agent $a^{m}_{i}$ is free to
visit vertices occupied by static or movable obstacles while abiding by the
constraints, whereas in part (ii) the search only permits $a^{m}_{i}$ to visit
vertices that abide by the constraints and do not include obstacles. After the
expansion of a node, new constraints are added to the child nodes, and their
costs are updated, as in CBS. In addition, a new lower bound $t$ for the
arrival of a mover agent $a^{i}_{m}$ to its movable obstacle is computed, and
its corresponding timed constraint is updated.

#### Theoretical properties and computational complexity

For a given assignment of mover agents to movable obstacles, TF-CBS derives
the properties of CBS (see (Sharon et al. 2015)) and is both complete and
optimal (proof omitted). Analyzing the running time is somewhat more
complicated. Gordon, Filmus, and Salzman (2021) show that the number of CT
node expansions for CBS may be as high as
$O\left((e|V|)^{|\mathcal{A}|C^{*}}\right)$, with $C^{*}$ being the cost of
the optimal solution. While our problem is not exactly the same, in our
setting we need to account for the agents $\mathcal{A}^{m}$, which incur an
exponential price in the worst case. An exact analysis is left for future
work.

###  5.2 A PBS-based approach for tMAPF

Figure 2:  Illustration of the adverse side-effect of transitive ordering.
Starting with a priority ordering $a^{t}_{1}\prec a^{t}_{2}$ at a PT node, the
low-level planner computes a path that causes agent $a^{t}_{1}$ to collide
with mover $a^{m}_{1}$ (at $s_{3}^{mo}$) while $a^{t}_{2}$ waits a single
timestep. Since obstacle displacement cost by $a^{m}_{1}$ is greater than the
path savings of $a^{t}_{1}$ (yet not accounting for the fact that $a^{t}_{2}$
would also benefit from $a^{m}_{1}$ clearing the way), the collision between
$a^{t}_{1}$ and $a^{m}_{1}$ would impose a priority $a^{m}_{1}\prec
a^{t}_{1}$. Then, transitive ordering implies a total ordering
$\\{a^{m}_{1}\prec a^{t}_{1}\prec a^{t}_{2}\\}$ that causes both $a^{t}_{1}$
and $a^{t}_{2}$ to avoid colliding with the path of mover $a^{m}_{1}$, never
utilizing the shortcut.

We describe a PBS-based approach for tMAPF called TF-PBS. Similarly to TF-CBS,
in TF-PBS we associate each mover agent $a^{m}_{i}\in\mathcal{A}^{m}$ with a
particular movable obstacle $o^{m}_{i}\in\mathcal{O}^{m}$ and treat those two
elements as one entity. On the high-level we explore a priority tree (PT)
whose nodes encode priority sets between task and mover agents, as in PBS. The
low-level search for a given agent proceeds in a manner similar to TF-CBS in
that for a task agent it searches for a path from its start to goal, and for a
mover agent $a_{i}^{m}$ it searches for a path from its start location
$s_{i}^{ma}$ to its obstacle’s location $s_{i}^{mo}$, and finishes at
$s_{i}^{mo}$.

We do, however, make a departure from the way priorities are treated in the
high-level search within PBS to accommodate the special circumstances of
tMAPF: Recall that when solving the MAPF problem, PBS enforces transitive
ordering of priorities induced by a topological sort. For instance, given a
priority set $\\{a_{3}\prec a_{1},~{}a_{1}\prec a_{2}\\}$, for some node $N$
of the PT tree and some three agents
$a_{1},a_{2},a_{3}\in\mathcal{A}^{t}\cup\mathcal{A}^{m}$, a topological sort
would yield the total order $a_{3}\prec a_{1}\prec a_{2}$, which implies that
the low-level search would compute a path for $a_{1}$ while avoiding collision
with $a_{3}$, and a path for $a_{2}$ while avoiding both $a_{3}$ and $a_{1}$.

Such an approach can lead to poor utilization of the movers’ capabilities of
clearing shortcuts in tMAPF. Suppose that agents $a_{1},a_{2}$ are task agents
$a^{t}_{1},a^{t}_{2}$ and $a_{3}$ is a mover $a^{m}_{1}$ whose shortest path
in the low-level search is such that it would not move once it reaches its
movable obstacle $o^{m}_{3}$. Considering that $a^{m}_{1}\prec a^{t}_{1}$,
agent $a^{t}_{1}$ would opt for a detour around a passage blocked by
$o^{m}_{3}$. Due to transitive ordering, $a^{t}_{2}$ would also avoid the
blocked passage, even if going through it would significantly improve its path
cost. This phenomenon is extended to the descendants of the current node, and
could be exacerbated with a few more agents. In more general terms,
$a^{m}_{1}$ chooses to block a shortcut, even though the shortcut can
potentially serve multiple agents and the aggregate path cost savings across
agents moving through the shortcut can offset the cost of obstacle
displacement.

In contrast, by reasoning about the exact priorities given in the priority set
(rather than transitive ordering) we can avoid such a situation. For instance,
if $a^{t}_{2}$ plans while avoiding collision only with $a^{t}_{1}$ (due to
the constraint $a^{t}_{1}\prec a^{t}_{2}$) it may choose to go through the
shortcut. This would cause a collision in the current PT node between
$a^{m}_{1}$ and $a^{t}_{2}$, which could be resolved in offspring PT nodes.
Fig. 2 depicts the problem with transitive ordering and illustrates the
phenomenon from the previous paragraph.

Thus, within TF-PBS we employ a combined approach that utilizes direct
prioritization whenever agents collide, and a greedy depth-first expansion
scheme that steers the high-level search towards nodes with lower node cost in
tie-breaks. Namely, the high-level search produces two child nodes and expands
first the child with lower cost among the two. It continues exploring the
descendants until a solution is found or until it reaches a node where the
low-level search deems that an agent cannot reach its goal due to higher-
ranking agents blocking the way.

(a)

(b)

Figure 3:  Warehouses used in our empirical evaluation. (3(a)) and (3(b))
depict the Large and Small environments, with $80$ and $50$ task agents
(colorful dots), respectively. Here rows of shelves (gray rectangles) form
long narrow aisles and goals (colorful rectangles) are selected from nearby
workstations around the periphery of each map.

##  6 Evaluation

A typical warehouse presents long rows of shelves that form narrow aisles,
with workstations located around the perimeter of the map (as illustrated in
Fig. 3). In autonomous warehouses, longer aisles allow for greater storage
capacity, but also lead to constrained environments that quickly become
congested as more agents are introduced. Hence, we evaluate our approach using
maps inspired by autonomous warehouses with intentionally long aisles, and
assess the impact of terraforming on measures of solution quality, node
expansions, and success-rate. We conduct our experiments on two warehouse-like
maps333Code and data will be made public upon publication.:

  * •

Small (Felner et al. 2018; Li et al. 2019) of size $24\times 47$.

  * •

Large (Li et al. 2020) of size $32\times 75$.

For each map we vary the number of task-agents $|A^{t}|$, and for every
combination of map and $|A^{t}|$ we generate $10$ scenarios where the agents’
start vertices are uniformly distributed. The agents’ goal vertices are
randomly selected to be either (1) around workstations (at the perimeter of
the map) or (2) empty vertices across the entire map (there is a 50/50 chance
to choose from (1) or (2)). In this manner, we obtain a flow of agents both to
and from workstations.

For tMAPF, recall that we make a simplifying assumption that our input also
includes a set of movable obstacles $\mathcal{O}^{m}$ as well as a set of
mover agents $\mathcal{A}^{m}$, and require that
$|\mathcal{A}^{m}|=|\mathcal{O}^{m}|$. Across all experiments, a set of $20$
and $42$ movable obstacles are selected for the Small and Large maps
respectively, situated in the middle of every row of shelves. When comparing
tMAPF with MAPF, the static environment treats all obstacles as static and
solves strictly for task agents $\mathcal{A}^{t}$ (omitting mover agents
$\mathcal{A}^{m}$).

We implemented the algorithms in Python and tested them on an Ubuntu machine
with 4GB RAM and a 2.7GHz Intel i7 CPU.

Figure 4:  Solution quality comparison between MAPF-optimal solutions, PBS
solutions, and tMAPF solutions relative to the _baseline_ solution (dashed
line).

#### Solution quality.

We report solution quality as the _sum of costs_ , as shown in Figure 4. The
horizontal axis specifies the Small and Large warehouse maps, denoted
Small-$n$ and Large-$n$, where $n=|\mathcal{A}^{t}|$. The vertical axis
measures solution suboptimality relative to a lower-bound MAPF solution of
ideal single-agent shortest paths, called _baseline_. The baseline (dashed
line) does not account for inter-agent collisions and interactions, which
regularly degrade from the quality of the optimal solution. Fig. 4 shows the
optimal MAPF solution (obtained with CBS) as being above the baseline due to
delays and congestion.

Next, we compare the solution cost obtained by TF-PBS for the tMAPF problem.
Recall that here there are additional mover agents ($20$ for Small and $42$
for Large) to be considered by the search. Fig. 4 shows both the upper-bound
cost (Eq. (2)) and the lower-bound cost (Eq. (1)). Recall that the lower-bound
reflects the path cost of regular agents and of mover agents as they carry
their assigned obstacle, whereas the upper-bound also accounts for movers’
path cost en-route to their obstacle. The results suggest that terraforming
has the capacity to _outperform_ the optimal solution available for a static
environment by alleviating bottlenecks and long detours. Interestingly, the
figure also demonstrates the potential to even surpass the baseline of a given
scenario, thanks to the ability of mover agents in tMAPF to create shortcuts.

#### Success rate and node expansions.

Success rate is measured by the number of scenarios solved within a $5$ minute
timeout, and use the number of (high-level) expanded nodes as a proxy for
average runtime necessary for each approach to obtain a solution. Figure 5
shows the success rate and expanded node count for the various combinations of
maps and number of task agents $|\mathcal{A}^{t}|$. We re-iterate that for
tMAPF, additional mover agents are introduced ($20$ for Small and $42$ for
Large). As a result, the total agent count $|\mathcal{A}|$ for tMAPF is
greater, which makes it more computationally challenging to solve than MAPF.
As expected from a congested warehouse environment, we see a rapid
deterioration in success rate as more agents are introduced. This is
especially pronounced in TF-CBS than CBS. In terms of the number of high-level
node expansions, we see a steep rise in CBS, and to a greater extent, in TF-
CBS.

Figure 5:  Comparison of success rate (top) and expanded nodes (bottom)
between TF-CBS, CBS, TF-PBS and PBS.

An encouraging trend is evident with TF-PBS, which is more robust in terms of
success rate and node expansions relative to the other approaches. The
warehouse environment presents narrow corridors that are prone to head-on
collisions between two or more agents. This is where the adaptive agent-
priority assignment of TF-PBS has a demonstrable advantage. When two agents
collide, TF-PBS imposes a priority ordering that prevents future collisions
between them. In this manner, local collisions between agents elicit lasting
constraints on the search space that partially eliminate unnecessary
explorations of the Priority Tree (PT), thus empirically reducing the number
of node expansions necessary to reach a solution.

##  7 Discussion and Future Work

In this work we explored the potential that terraforming has for MAPF-like
problems. As demonstrated in our evaluation, the result is a form of emergent
collaboration, in which agents create shortcuts and reduce the overall cost
beyond what can be achieved in the same static environment.

To harness the full potential of tMAPF, we envision its application to L-MAPF
where agents en-route to collect an item can serve as mover agents, hence
moving one or more movable objects with minimal overhead. This blurs the
simplifying assumption we made in our problem formulation where there is a
clear distinction between mover and task agents, and motivates the cost
functions we introduced for tMAPF. This also requires lifting assumptions A1
and A2.

The main challenge we see is how to dynamically assign movable obstacles to
agents that are not currently carrying items. We consider two alternatives: In
the first alternative, agents not carrying items plan to their original goal
and are assigned movable obstacles that lie on or near their already-planned
path. This will probably incur little overhead for the agent but may result in
an assignment to a movable obstacle whose displacement offers only a minor
benefit. In the second alternative, task agents that identify that moving a
movable obstacle will dramatically reduce their path length will “request”
that the obstacle be moved. This will trigger an assignment of that specific
obstacle to the closest agent not carrying an item. This can be seen as the
complementary case to the first alternative where moving an obstacle may incur
a large overhead for the agent but results in an assignment to a movable
obstacle that can significantly improve costs if moved.

Finally, we foresee applications of our work beyond our motivating example of
MAPF in warehouses. For example, consider sortation centers (see (Kou et al.
2020) and visualization in Fig. 6) where agents need to a reach a dropoff
station, obtain a parcel for delivery, and then deliver the parcel to a
sorting bin. We envision the sorting bins as having a mechanism that can be
used to automatically cover them (possibly incurring time). This will allow
agents to pass over covered bins hence reducing their path length. In the
context of our work, these bins are the movable obstacles considered in this
work. There are several challenges that need to be addressed (both mechanical
as well as algorithmic), but this showcases the potential power of
terraforming for increasing throughput of automated logistic centers.

Figure 6:  Sortation center. Figure adapted from %**** ̵tf˙mapf˙socs˙22.tex
̵Line ̵1050 ̵****https://tinyurl.com/bdj8n8kx.

## References

  * Barer et al. (2014) Barer, M.; Sharon, G.; Stern, R.; and Felner, A. 2014.  Suboptimal Variants of the Conflict-Based Search Algorithm for the Multi-Agent Pathfinding Problem.  In _European Conference on Artificial Intelligence (ECAI)_ , volume 263, 961–962.
  * Bellusci, Basilico, and Amigoni (2020) Bellusci, M.; Basilico, N.; and Amigoni, F. 2020.  Multi-Agent Path Finding in Configurable Environments.  In _Autonomous Agents and MultiAgent Systems (AAMAS)_ , 159–167.
  * Belov et al. (2020) Belov, G.; Du, W.; de la Banda, M. G.; Harabor, D.; Koenig, S.; and Wei, X. 2020\.  From Multi-Agent Pathfinding to 3D Pipe Routing.  In Harabor, D.; and Vallati, M., eds., _Symposium on Combinatorial Search (SoCS)_ , 11–19.
  * Boyarski et al. (2015) Boyarski, E.; Felner, A.; Stern, R.; Sharon, G.; Tolpin, D.; Betzalel, O.; and Shimony, S. E. 2015.  ICBS: Improved Conflict-Based Search Algorithm for Multi-Agent Pathfinding.  In _International Joint Conferences on Artificial Intelligence (IJCAI)_ , 740–746.
  * Choudhury et al. (2021) Choudhury, S.; Solovey, K.; Kochenderfer, M. J.; and Pavone, M. 2021.  Efficient Large-Scale Multi-Drone Delivery using Transit Networks.  _Journal of Artificial Intelligence Research_ , 70: 757–788.
  * Erdem et al. (2013) Erdem, E.; Kisa, D. G.; Oztok, U.; and Schüller, P. 2013.  A general formal framework for pathfinding problems with multiple agents.  In _Association for the Advancement of Artificial Intelligence (AAAI)_.
  * Felner et al. (2018) Felner, A.; Li, J.; Boyarski, E.; Ma, H.; Cohen, L.; Kumar, T. S.; and Koenig, S. 2018.  Adding heuristics to conflict-based search for multi-agent path finding.  In _International Conference on Automated Planning and Scheduling (ICAPS)_ , volume 28.
  * Gordon, Filmus, and Salzman (2021) Gordon, O.; Filmus, Y.; and Salzman, O. 2021.  Revisiting the Complexity Analysis of Conflict-Based Search: New Computational Techniques and Improved Bounds.  In _Symposium on Combinatorial Search (SoCS)_ , 64–72.
  * Greshler et al. (2021) Greshler, N.; Gordon, O.; Salzman, O.; and Shimkin, N. 2021.  Cooperative Multi-Agent Path Finding: Beyond Path Planning and Avoidance.  In _Symposium on Multi-Robot and Multi-Agent Systems (MRS)_ , 20–28.
  * Hart, Nilsson, and Raphael (1968) Hart, P. E.; Nilsson, N. J.; and Raphael, B. 1968.  A Formal Basis for the Heuristic Determination of Minimum Cost Paths.  _IEEE Transactions on Systems Science and Cybernetics_ , 4(2): 100–107.
  * Hauser (2014) Hauser, K. 2014.  The minimum constraint removal problem with three robotics applications.  _The International Journal of Robotics Research_ , 33(1): 5–17.
  * Hauser (2013) Hauser, K. K. 2013.  Minimum constraint displacement motion planning.  In _Robotics: science and systems_.
  * Kou et al. (2020) Kou, N. M.; Peng, C.; Ma, H.; Kumar, T. K. S.; and Koenig, S. 2020.  Idle Time Optimization for Target Assignment and Path Finding in Sortation Centers.  In _Association for the Advancement of Artificial Intelligence (AAAI)_ , 9925–9932.
  * Kuhn (1955) Kuhn, H. W. 1955.  The Hungarian method for the assignment problem.  _Naval research logistics quarterly_ , 2(1-2): 83–97.
  * Li et al. (2021a) Li, J.; Chen, Z.; Zheng, Y.; Chan, S.; Harabor, D.; Stuckey, P. J.; Ma, H.; and Koenig, S. 2021a.  Scalable Rail Planning and Replanning: Winning the 2020 Flatland Challenge.  In _International Conference on Automated Planning and Scheduling (ICAPS)_ , 477–485.
  * Li et al. (2020) Li, J.; Gange, G.; Harabor, D.; Stuckey, P. J.; Ma, H.; and Koenig, S. 2020.  New techniques for pairwise symmetry breaking in multi-agent path finding.  In _International Conference on Automated Planning and Scheduling (ICAPS)_ , volume 30, 193–201.
  * Li et al. (2019) Li, J.; Harabor, D.; Stuckey, P. J.; Felner, A.; Ma, H.; and Koenig, S. 2019.  Disjoint splitting for multi-agent path finding with conflict-based search.  In _International Conference on Automated Planning and Scheduling (ICAPS)_ , volume 29, 279–283.
  * Li et al. (2021b) Li, J.; Tinka, A.; Kiesel, S.; Durham, J.; Kumar, S.; and Koenig, S. 2021b.  Lifelong Multi-Agent Path Finding in Large-Scale Warehouses.  In _Association for the Advancement of Artificial Intelligence (AAAI)_.
  * Liu et al. (2019) Liu, M.; Ma, H.; Li, J.; and Koenig, S. 2019.  Task and path planning for multi-agent pickup and delivery.  In _Autonomous Agents and MultiAgent Systems (AAMAS)_.
  * Ma et al. (2019) Ma, H.; Harabor, D.; Stuckey, P. J.; Li, J.; and Koenig, S. 2019.  Searching with consistent prioritization for multi-agent path finding.  In _Association for the Advancement of Artificial Intelligence (AAAI)_ , volume 33, 7643–7650.
  * Ma et al. (2017) Ma, H.; Li, J.; Kumar, T. K. S.; and Koenig, S. 2017.  Lifelong multi-agent path finding for online pickup and delivery tasks.  In _Autonomous Agents and MultiAgent Systems (AAMAS)_ , 837–845.
  * Salzman and Stern (2020) Salzman, O.; and Stern, R. 2020.  Research Challenges and Opportunities in Multi-Agent Path Finding and Multi-Agent Pickup and Delivery Problems.  In _Autonomous Agents and MultiAgent Systems (AAMAS)_ , 1711–1715.
  * Sharon et al. (2015) Sharon, G.; Stern, R.; Felner, A.; and Sturtevant, N. R. 2015.  Conflict-based search for optimal multi-agent pathfinding.  _Artificial Intelligence_ , 219: 40–66.
  * Stern et al. (2019) Stern, R.; Sturtevant, N.; Felner, A.; Koenig, S.; Ma, H.; Walker, T.; Li, J.; Atzmon, D.; Cohen, L.; Kumar, S.; et al. 2019.  Multi-Agent Pathfinding: Definitions, Variants, and Benchmarks.  In _Symposium on Combinatorial Search (SoCS)_.
  * Surynek et al. (2016) Surynek, P.; Felner, A.; Stern, R.; and Boyarski, E. 2016.  Efficient SAT approach to multi-agent path finding under the sum of costs objective.  In _European Conference on Artificial Intelligence (ECAI)_ , 810–818. IOS Press.
  * Tang, Han, and Yu (2020) Tang, W. N.; Han, S. D.; and Yu, J. 2020.  Computing High-Quality Clutter Removal Solutions for Multiple Robots.  In _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 7963–7970.
  * Wurman, D’Andrea, and Mountz (2008) Wurman, P. R.; D’Andrea, R.; and Mountz, M. 2008.  Coordinating Hundreds of Cooperative, Autonomous Vehicles in Warehouses.  _AI Magazine_ , 29(1): 9.
  * Yu and LaValle (2012) Yu, J.; and LaValle, S. M. 2012.  Multi-agent Path Planning and Network Flow.  In _Workshop on the Algorithmic Foundations of Robotics (WAFR)_ , volume 86, 157–173. Springer.