# Canonical symplectic particle-in-cell method for long-term large-scale
simulations of the Vlasov-Maxwell system

Hong Qin  School of Nuclear Science and Technology and Department of Modern
Physics, University of Science and Technology of China, Hefei, Anhui 230026,
China  Plasma Physics Laboratory, Princeton University, Princeton, NJ 08543
Jian Liu  School of Nuclear Science and Technology and Department of Modern
Physics, University of Science and Technology of China, Hefei, Anhui 230026,
China  Key Laboratory of Geospace Environment, CAS, Hefei, Anhui 230026, China
Jianyuan Xiao  School of Nuclear Science and Technology and Department of
Modern Physics, University of Science and Technology of China, Hefei, Anhui
230026, China  Key Laboratory of Geospace Environment, CAS, Hefei, Anhui
230026, China  Ruili Zhang  School of Nuclear Science and Technology and
Department of Modern Physics, University of Science and Technology of China,
Hefei, Anhui 230026, China  Key Laboratory of Geospace Environment, CAS,
Hefei, Anhui 230026, China  Yang He  School of Nuclear Science and Technology
and Department of Modern Physics, University of Science and Technology of
China, Hefei, Anhui 230026, China  Key Laboratory of Geospace Environment,
CAS, Hefei, Anhui 230026, China  Yulei Wang  School of Nuclear Science and
Technology and Department of Modern Physics, University of Science and
Technology of China, Hefei, Anhui 230026, China  Key Laboratory of Geospace
Environment, CAS, Hefei, Anhui 230026, China  Yajuan Sun  LSEC, Academy of
Mathematics and Systems Science, Chinese Academy of Sciences, P.O. Box 2719,
Beijing 100190, China  Joshua W. Burby  Plasma Physics Laboratory, Princeton
University, Princeton, NJ 08543  Leland Ellison  Plasma Physics Laboratory,
Princeton University, Princeton, NJ 08543  Yao Zhou  Plasma Physics
Laboratory, Princeton University, Princeton, NJ 08543

###### Abstract

Particle-in-Cell (PIC) simulation is the most important numerical tool in
plasma physics. However, its long-term accuracy has not been established. To
overcome this difficulty, we developed a canonical symplectic PIC method for
the Vlasov-Maxwell system by discretizing its canonical Poisson bracket. A
fast local algorithm to solve the symplectic implicit time advance is
discovered without root searching or global matrix inversion, enabling
applications of the proposed method to very large-scale plasma simulations
with many, e.g., $10^{9}$, degrees of freedom. The long-term accuracy and
fidelity of the algorithm enables us to numerically confirm Mouhot and
Villani’s theory and conjecture on nonlinear Landau damping over several
orders of magnitude using the PIC method, and to calculate the nonlinear
evolution of the reflectivity during the mode conversion process from
extraordinary waves to Bernstein waves.

Particle-in-Cell Simulations, Vlasov-Maxwell Equations, Canonical Symplectic
Algorithm

###### pacs:

52.65.Rr, 52.25.Dg

In modern plasma physics, numerically solving the Vlasov-Maxwell (VM)
equations using the Particle-In-Cell (PIC) method has become the most
important tool Birdsall91 ; Hockney88  for theoretical studies in the last
half century. Many innovative algorithms, such as the Boris scheme for
advancing particles Boris70 ; Boris77  and Villasenor-Buneman’s charge-
conserving deposition scheme Villasenor92 , have been developed and
successfully applied. Recently, new geometric numerical methodology has been
adopted for PIC simulations. This exciting trend begins with the discovery of
symplectic algorithms for Hamiltonian equations that govern charged particle
dynamics Ruth83 ; Feng85 ; Feng86 ; Feng10 ; Forest90 ; Channell90 ; Candy91 ;
Marsden01-357 ; Hairer02 ; Qin08-PRL ; Qin09-PoP . The Boris algorithm was
discovered to be volume-preserving Qin13-084503 ; Zhang15-pop  and high-order
volume-preserving methods have been found He15 . In addition, the Vlasov-
Maxwell system Squire12 ; Xiao13 ; Kraus14 ; Shadwick14 ; Xiao15  and the
Vlasov-Poisson system Evstatiev13  have been discretized from a variational
symplectic perspective that preserves symplectic structures and exhibits
excellent long-term accuracy and fidelity.

In this letter, we develop a new canonical symplectic PIC method for solving
the VM equations by discretizing its canonical Poisson bracket Marsden82 . The
distribution function $f$ is first discretized in phase space through the
Klimontovich representation by a finite number of Lagrangian sampling points
$(\bm{\mathbf{X}}_{i},\mathbf{P}_{i})$ $(i=1,...,N)$, where $\mathbf{X}_{i}$
and $\mathbf{P}_{i}$ are the position and canonical momentum of the $i$-th
particle, and $N$ is the total number of sampling points. The electromagnetic
field is discretized point-wise on a given spatial grid, and the Hamiltonian
functional is expressed as a function of the sampling points and the
discretized electromagnetic field. This procedure generates a finite-
dimensional Hamiltonian system with a canonical symplectic structure. The
number of degrees of freedom for the discrete system is $D=3N+3M,$ where $M$
denotes the total number of the discrete grid-points.

In general, for a Hamiltonian function whose momentum dependence and position
dependence are not separable, it is not possible to make symplectic
integration algorithms explicit Feng10 ; Hairer02 . For the discrete
Hamiltonian system developed here for the VM equations, the dimension of
system is usually very large, and root searching algorithms required by
implicit methods are too time-consuming to be practical. However, we
discovered that if the symplectic Euler algorithm Feng10  is applied to the
discrete VM Hamiltonian system at hand, the implicit time advance can be
carried out as inexpensively as an explicit method by just inverting a
$3\times 3$ matrix for every particle separately. The resulting canonical
symplectic PIC method for the VM system inherits all the good numerical
features of canonical symplectic algorithms, such as the long-term bound on
energy-momentum error. Being symplectic means that the numerical solution
satisfies $D(2D-1)$ constraints as the exact solution does. Since $D$ is a
large number, the symplectic condition is much stronger than a few constraints
on global energy and momentum. The symplectic condition is almost as strong as
imposing local conservation everywhere in phase space.

Two examples of application are given. In the first example, we simulate the
dynamics of nonlinear Landau damping. It also serves as a test of the
algorithm. The discrete VM Hamiltonian system for this study has more than
$2.69\times 10^{8}$ degrees of freedom. The damping rate from the numerical
results agrees exactly with the theoretical result. Furthermore, long-term
simulations reveal that the phase mixing dynamics in velocity space is the
physical mechanism of the nonlinear Landau damping, as recently proved by
Mouhot and Villani Mouhot11 ; Villani14  for the Vlasov-Poisson system and
conjectured by Villani Villani14-i  for the Vlasov-Maxwell system. In the
second application, we study the nonlinear mode conversion process from
extraordinary modes to Bernstein modes (X-B mode conversion) in an
inhomogeneous hot plasma. Simulations show that nonlinear mode excitations and
self-interaction of the Bernstein waves significantly modify the reflectivity
and conversion rate. It is the long-term accuracy and fidelity of the
canonical symplectic PIC algorithm that enables us to numerically confirm
Mouhot and Villani’s theory and conjecture over several orders of magnitude
using the PIC method, and to calculate the nonlinear evolution of the
reflectivity during the X-B mode conversion.

We start from the canonical Poission bracket and Hamiltonian for the Vlasov-
Maxwell equations Marsden82 ,

$\displaystyle\\{F,G\\}\equiv\int f\left\\{\frac{\delta F}{\delta
f},\frac{\delta G}{\delta
f}\right\\}{}_{\mathbf{xp}}d\mathbf{x}d\mathbf{p}+\int\left(\frac{\delta
F}{\delta\mathbf{A}}\frac{\delta G}{\delta\mathbf{Y}}-\frac{\delta
G}{\delta\mathbf{A}}\frac{\delta F}{\delta\mathbf{Y}}\right)d\mathbf{x}.$ (1)
$\displaystyle
H(f,\mathbf{A},\mathbf{Y})=\frac{1}{2}\int(\mathbf{p}-\mathbf{A})^{2}fd\mathbf{x}d\mathbf{p}+\frac{1}{2}\int\left[\mathbf{Y}^{2}+(\nabla\times\mathbf{A})^{2}\right]d\mathbf{x}.$
(2)

Here, $F$, $G$, and the Hamiltonian $H$ are functionals of the distribution
function $f$, vector potential $\mathbf{A}$, and
$\mathbf{Y}\equiv\partial\mathbf{A}/\partial t$. The bracket
$\\{h,g\\}_{\mathbf{xp}}$ inside the first term on the right hand side of Eq.
(1) is the canonical Poisson bracket for functions $h$ and $g$ of canonical
phase space $(\mathbf{x},\mathbf{p})$. The temporal gauge, i.e., $\phi=0$, has
been explicitly chosen for this Poisson bracket to be valid. The Poisson
bracket defined in Eq. (1) can be formally derived from the point of view of
co-adjoint orbit theory Marsden82 , and can be used to derive the non-
canonical Morrison-Marsden-Weinstein bracket in the
$(f,\mathbf{E},\mathbf{B})$ and $(\mathbf{x},\mathbf{v})$ coordinates
Morrison80 ; Weinstein81 ; Marsden82 . First, we discretize the distribution
function using the Klimontovich representation

$f(\mathbf{x},\mathbf{p},t)=\sum_{i=1}^{N}\delta\left(\mathbf{x}-\mathbf{X}_{i}\right)\delta\left(\mathbf{p}-\mathbf{P}_{i}\right),$
(3)

where $(\mathbf{X}_{i},\mathbf{P}_{i})$ $(i=1,...,N)$ are particles’
coordinates in phase space. Under this discretization, it can be shown that

$\displaystyle\frac{\delta F}{\delta\mathbf{X}_{i}}$ $\displaystyle=$
$\displaystyle\int\delta\left(\mathbf{x}-\mathbf{X}_{i}\right)\delta\left(\mathbf{p}-\mathbf{P}_{i}\right)\frac{\partial}{\partial\mathbf{x}}\left(\frac{\delta
F}{\delta f}\right)d\mathbf{x}d\mathbf{p},$ (4) $\displaystyle\frac{\delta
F}{\delta\mathbf{P}_{i}}$ $\displaystyle=$
$\displaystyle\int\delta\left(\mathbf{x}-\mathbf{X}_{i}\right)\delta\left(\mathbf{p}-\mathbf{P}_{i}\right)\frac{\partial}{\partial\mathbf{p}}\left(\frac{\delta
F}{\delta f}\right)d\mathbf{x}d\mathbf{p},$ (5)

from which we obtain

$\frac{\delta F}{\delta\mathbf{X}_{i}}\frac{\delta
G}{\delta\mathbf{P}_{i}}-\frac{\delta G}{\delta\mathbf{X}_{i}}\frac{\delta
F}{\delta\mathbf{P}_{i}}=\int\delta\left(\mathbf{x}-\mathbf{X}_{i}\right)\delta\left(\mathbf{p}-\mathbf{P}_{i}\right)\left\\{\frac{\delta
F}{\delta f},\frac{\delta G}{\delta
f}\right\\}{}_{\mathbf{xp}}d\mathbf{x}d\mathbf{p}.$ (6)

It then follows that the first term on the right-hand side of Eq. (1) is

$\int f\left\\{\frac{\delta F}{\delta f},\frac{\delta G}{\delta
f}\right\\}{}_{\mathbf{xp}}d\mathbf{x}d\mathbf{p}=\sum_{i=1}^{N}\left(\frac{\delta
F}{\delta\mathbf{X}_{i}}\frac{\delta G}{\delta\mathbf{P}_{i}}-\frac{\delta
G}{\delta\mathbf{X}_{i}}\frac{\delta F}{\delta\mathbf{P}_{i}}\right).$ (7)

Similar derivation of the discretized bracket can be found in the context of
Hamiltonian description of vortex fluid Morrison81-1783 ; Morrison81-1788 .

To discretize the second term on the right-hand side of Eq. (1), we first
discretize the fields $\mathbf{A}(t)$ and $\mathbf{Y}(t)$ on a Eulerian
spatial grid as

$\mathbf{A}(\mathbf{x},t)=\sum_{J=1}^{M}\mathbf{A}_{\mathbf{J}}(t)\Psi(\mathbf{x}-\mathbf{x}_{J})\,,\quad\mathbf{Y}(\mathbf{x},t)=\sum_{J=1}^{M}\mathbf{Y}_{J}(t)\Psi(\mathbf{x}-\mathbf{x}_{J})\,,$
(8)

where the discrete fields $\mathbf{A}_{J}(t)$ and $\mathbf{Y}_{J}(t)$ are the
fields evaluated on the grid-point $\mathbf{x}_{J}$. The subscript $J$ is the
index of the grid-point, and $M$ is the total number of the grid-points. Here,
$\Psi(\mathbf{x}-\mathbf{x}_{J})$ is the step function,

$\Psi(\mathbf{x}-\mathbf{x}_{J})=\left\\{\begin{array}[]{lc}1,&|x-x_{J}|<\frac{\Delta
x}{2},|y-y_{J}|<\frac{\Delta y}{2},|z-z_{J}|<\frac{\Delta z}{2}~{},\\\
0,&\textrm{elsewhere}~{}.\end{array}\right.$ (9)

Under this discretization of $\mathbf{A}(t)$ and $\mathbf{Y}(t)$, we have

$\frac{\delta
F}{\delta\mathbf{A}}=\sum_{J=1}^{M}\frac{\delta\mathbf{A}_{J}}{\delta\mathbf{A}}\frac{\partial
F}{\partial\mathbf{A}_{J}}=\sum_{J=1}^{M}\frac{1}{\Delta
V}\Psi(\mathbf{x}-\mathbf{x}_{J})\frac{\partial
F}{\partial\mathbf{A}_{J}}\thinspace,$ (10)

where $\Delta V$ is the volume of each cell, which is taken to be a constant
in the present study. For the second equal sign in Eq. (10), use has been made
of the fact that

$\frac{\delta\mathbf{A}_{J}}{\delta\mathbf{A}}=\frac{1}{\Delta
V}\Psi(\mathbf{x}-\mathbf{x}_{J})\thinspace.$ (11)

The discretization of the second term on the right-hand side of Eq. (1) is
thus __

$\int\left(\frac{\delta F}{\delta\mathbf{A}}\frac{\delta
G}{\delta\mathbf{Y}}-\frac{\delta G}{\delta\mathbf{A}}\frac{\delta
F}{\delta\mathbf{Y}}\right)d\mathbf{x}=\sum_{J=1}^{M}\left(\frac{\partial
F}{\partial\mathbf{A}_{J}}\frac{\partial
G}{\partial\mathbf{Y}_{J}}-\frac{\partial
G}{\partial\mathbf{A}_{J}}\frac{\partial
F}{\partial\mathbf{Y}_{J}}\right)\frac{1}{\Delta V}\thinspace.$ (12)

Finally, the discrete Poisson bracket for the VM system is

$\qquad\left\\{F,G\right\\}=\sum_{i=1}^{N}\left(\frac{\delta
F}{\delta\mathbf{X}_{i}}\frac{\delta G}{\delta\mathbf{P}_{i}}-\frac{\delta
G}{\delta\mathbf{X}_{i}}\frac{\delta
F}{\delta\mathbf{P}_{i}}\right)+\sum_{J=1}^{M}\left(\frac{\partial
F}{\partial\mathbf{A}_{J}}\frac{\partial
G}{\partial\mathbf{Y}_{J}}-\frac{\partial
G}{\partial\mathbf{A}_{J}}\frac{\partial
F}{\partial\mathbf{Y}_{J}}\right)\frac{1}{\Delta V}\thinspace$ (13)

for functions $F$ and $G$ of the particles $(\mathbf{X}_{i},\mathbf{P}_{i})$
and the discretized field $(\mathbf{A}_{J},\mathbf{Y}_{J})$.

Next, we need to express the Hamiltonian functional given by Eq. (2) in terms
of $(\mathbf{X}_{i},\mathbf{P}_{i})$ and $(\mathbf{A}_{J},\mathbf{Y}_{J})$.
The particles’ total kinetic energy is the sum of each particle’s kinetic
energy. The vector potential at a particle’s position can be interpolated from
$\mathbf{A}_{J}(t)$ as

$\mathbf{A}(\mathbf{X}_{j},t)=\sum_{J=1}^{M}\mathbf{A}_{J}(t)W(\mathbf{X}_{j}-\mathbf{x}_{J})\,,$
(14)

where $W(\mathbf{\mathbf{X}}_{j}-\mathbf{x}_{J})$ is a chosen interpolation
function. Note that $W(\mathbf{\mathbf{X}}_{j}-\mathbf{x}_{J})$ is not
necessarily the same as the step function $\Psi(\mathbf{x}-\mathbf{x}_{J})$ in
Eq. (8). This is of course allowed as long as the consistency condition is
satisfied, i.e., the continuous limit is recovered when the grid-size goes to
zero. The Hamiltonian then becomes

$\displaystyle\tilde{H}(\mathbf{X}_{i},\mathbf{P}_{i},\mathbf{A}_{J},\mathbf{Y}_{J})=\frac{1}{2}\sum_{i=1}^{N}\left(\mathbf{P}_{i}^{2}-2\mathbf{P}_{i}\cdot\sum_{J=1}^{M}\mathbf{A}_{J}W(\mathbf{X}_{i}-\mathbf{x}_{J})\right.$
(15)
$\displaystyle\left.\sum_{J,L=1}^{M}\mathbf{A}_{J}\cdot\mathbf{A}_{L}W(\mathbf{X}_{i}-\mathbf{x}_{J})W(\mathbf{X}_{i}-\mathbf{x}_{L})\right)+\frac{1}{2}\sum_{J=1}^{M}\left[\mathbf{Y}_{J}^{2}+\left(\nabla_{d}\times\mathbf{A}\right)_{J}^{2}\right]\Delta
V,$

where $\left(\nabla_{d}\times\mathbf{A}\right)_{J}$ is the discrete curl
operator acting on the discrete vector potential evaluated at the $J$-th grid-
point. Finally, the discrete Hamiltonian (15) and discrete Poisson structure
(13) form a canonical symplectic discretization of the original continuous
Vlasov-Maxwell system. The ordinary differential equations for the canonical
system are

$\displaystyle\dot{\mathbf{X}}_{i}$ $\displaystyle=$
$\displaystyle\left\\{\mathbf{X}_{i},\tilde{H}\right\\}_{d}=\mathbf{P}_{i}-\sum_{J=1}^{M}\mathbf{A}_{J}W(\mathbf{X}_{i}-\mathbf{x}_{J})\,,$
(16) $\displaystyle\dot{\mathbf{A}}_{J}$ $\displaystyle=$
$\displaystyle\left\\{\mathbf{A}_{J},\tilde{H}\right\\}_{d}=\mathbf{Y}_{J}\,,$
(17) $\displaystyle\dot{\mathbf{P}}_{i}$ $\displaystyle=$
$\displaystyle\left\\{\mathbf{P}_{i},\tilde{H}\right\\}_{d}=\sum_{J=1}^{M}\left(\mathbf{P}_{i}\cdot\mathbf{A}_{J}\right)\nabla
W(\mathbf{X}_{i}-\mathbf{x}_{J})-$
$\displaystyle\sum_{J,L=1}^{M}\left(\mathbf{A}_{J}\cdot\mathbf{A}_{L}\right)W(\mathbf{X}_{i}-\mathbf{x}_{J})\nabla
W(\mathbf{X}_{i}-\mathbf{x}_{L})\,,$ (18) $\displaystyle\dot{\mathbf{Y}}_{J}$
$\displaystyle=$
$\displaystyle\left\\{\mathbf{Y}_{J},\tilde{H}\right\\}_{d}=\sum_{i=1}^{N}\mathbf{P}_{i}W(\mathbf{X}_{i}-\mathbf{x}_{J})\frac{1}{\Delta
V}-$
$\displaystyle\sum_{i=1}^{N}\sum_{L=1}^{M}\mathbf{A}_{L}W(\mathbf{X}_{i}-\mathbf{x}_{J})W(\mathbf{X}_{i}-\mathbf{x}_{L})\frac{1}{\Delta
V}-\left(\nabla_{d}^{T}\times\nabla_{d}\times\mathbf{A}\right)_{J}\,.$ (19)

This equation system consists of $6(M+N)$ equations describing the dynamics of
$N$ particles and fields on $M$ discrete grid-points. The last term in Eq.
(19) is defined to be

$\left(\nabla_{d}^{T}\times\nabla_{d}\times\mathbf{A}\right)_{J}\equiv\frac{1}{2}\frac{\partial}{\partial
A_{J}}\left[\sum_{L=1}^{M}\left(\nabla_{d}\times\mathbf{A}\right)_{L}^{2}\right].$
(20)

The notation of $\nabla_{d}^{T}\times\nabla_{d}\times\mathbf{A}$ indicates
that the right hand side of Eq. (20) can be viewed as the discretized
$\nabla\times\nabla\times\mathbf{A}$ for a well-chosen discrete curl operator
$\nabla_{d}$. To wit, we note that the term $\nabla\times{\bf A}$ in the
Hamiltonian is discretized using the step function $\Psi({\bf x}-{\bf
x}_{J})$,

$\nabla\times{\bf A}=\sum_{J=1}^{M}\left(\nabla_{d}\times{\bf
A}\right)_{J}\Psi({\bf x}-{\bf x}_{J}).$ (21)

As an example, we define the discrete curl operator
$(\nabla_{d}\times\mathbf{A})_{J}=\left(\nabla\times\mathbf{A}\right)_{i,j,k}$
to be

$\left(\nabla_{d}\times\mathbf{A}\right)_{J}\equiv\left(\begin{array}[]{c}\frac{A_{i,j,k}^{3}-A_{i,j-1,k}^{3}}{\Delta
y}-\frac{A_{i,j,k}^{2}-A_{i,j,k-1}^{2}}{\Delta z}\\\
\frac{A_{i,j,k}^{1}-A_{i,j,k-1}^{1}}{\Delta
z}-\frac{A_{i,j,k}^{3}-A_{i-1,j,k}^{3}}{\Delta{\bf x}}\\\
\frac{A_{i,j,k}^{2}-A_{i-1,j,k}^{2}}{\Delta{\bf
x}}-\frac{A_{i,j,k}^{1}-A_{i,,j-1,k}^{1}}{\Delta y}\end{array}\right),$ (22)

which can be written as a linear operator on the space of $3M$-vectors as

$\nabla_{d}\times\mathbf{A}=\left(\begin{array}[]{c}\left(\nabla_{d}\times\mathbf{A}\right)_{1}\\\
\vdots\\\
\left(\nabla_{d}\times\mathbf{A}\right)_{M}\end{array}\right)=\Gamma\left(\begin{array}[]{c}\mathbf{A}_{1}\\\
\vdots\\\ \mathbf{A}_{M}\end{array}\right)\thinspace.$ (23)

Here, $\Gamma$ is a $3M\times 3M$ sparse matrix specifying the discrete curl
operator. The partial derivative with respect to $\mathbf{A}_{J}$ can be
expressed as

$\frac{1}{2}\frac{\partial}{\partial\mathbf{A}_{J}}\left[\sum_{L=1}^{M}(\nabla_{d}\times\mathbf{A})_{L}^{2}\right]=\left[\Gamma^{T}\Gamma\left(\begin{array}[]{c}\mathbf{A}_{1}\\\
\vdots\\\
\mathbf{A}_{M}\end{array}\right)\right]_{J}\equiv\left(\nabla_{d}^{T}\times\nabla_{d}\times\mathbf{A}\right)_{J},$
(24)

where $\Gamma^{T}$ is the transposition of $\Gamma$. Obviously, the notation
in Eq. (20) or (24) is meaningful for any linear discrete curl operator
$\nabla_{d}$.

It is clear from Eqs. (16)-(19) that the particles and fields interact through
the interpolation function $W(\mathbf{X}_{i}-\mathbf{x}_{J})$. The function
$W(\mathbf{X}_{i}-\mathbf{x}_{L})/\Delta V$ distributes particles’ charge over
the grid-points as if they are “charged clouds” with finite-size Birdsall91 .

Once the canonical symplectic structure is given, canonical symplectic
algorithms can be readily constructed using well-developed methods Ruth83 ;
Feng85 ; Feng86 ; Feng10 ; Channell90 ; Forest90 ; Candy91 ; Marsden01-357 .
For a reason soon to be clear, we adopt the semi-explicit symplectic Euler
method for time advance. The symplectic Euler method for a generic canonical
Hamiltonian system is

$\displaystyle p^{n+1}$ $\displaystyle=$ $\displaystyle p^{n}-\Delta
t\frac{\partial H}{\partial q}(p^{n+1},q^{n}),$ (25) $\displaystyle q^{n+1}$
$\displaystyle=$ $\displaystyle q^{n}+\Delta t\frac{\partial H}{\partial
p}(p^{n+1},q^{n}).$ (26)

where $\Delta t$ is the time-step, and the superscript $n$ in $p^{n}$ and
$q^{n}$ denotes that they are the value at the $n$-th time step. It is
implicit for $p$, but explicit for $q$. Making use of this algorithm, the
iteration rules for Eqs. (16) -(19) are

$\displaystyle\frac{\mathbf{X}_{i}^{n+1}-\mathbf{X}_{i}^{n}}{\Delta
t}=\mathbf{P}_{i}^{n+1}-\sum_{J=1}^{M}\mathbf{A}_{J}^{n}W(\mathbf{X}_{i}^{n}-\mathbf{x}_{J})\,,$
(27) $\displaystyle\frac{\mathbf{A}_{J}^{n+1}-\mathbf{A}_{J}^{n}}{\Delta
t}=\mathbf{Y}_{J}^{n+1}\,,$
$\displaystyle\frac{\mathbf{P}_{i}^{n+1}-\mathbf{P}_{i}^{n}}{\Delta
t}=\sum_{J=1}^{M}\left(\mathbf{P}_{i}^{n+1}\cdot\mathbf{A}_{J}^{n}\right)\nabla
W(\mathbf{X}_{i}^{n}-\mathbf{x}_{J})$
$\displaystyle\qquad-\sum_{J,L=1}^{M}\left(\mathbf{A}_{J}^{n}\cdot\mathbf{A}_{L}^{n}\right)W(\mathbf{X}_{i}^{n}-\mathbf{x}_{J})\nabla
W(\mathbf{X}_{i}^{n}-\mathbf{x}_{L})\,,$
$\displaystyle\frac{\mathbf{Y}_{J}^{n+1}-\mathbf{Y}_{J}^{n}}{\Delta
t}=\sum_{i=1}^{N}\mathbf{P}_{i}^{n+1}W(\mathbf{X}_{i}^{n}-\mathbf{x}_{J})\frac{1}{\Delta
V}$
$\displaystyle\qquad-\sum_{i=1}^{N}\sum_{L=1}^{M}\mathbf{A}_{L}^{n}W(\mathbf{X}_{i}^{n}-\mathbf{x}_{J})W(\mathbf{X}_{i}^{n}-\mathbf{x}_{L})\frac{1}{\Delta
V}-\left(\nabla_{d}^{T}\times\nabla_{d}\times\mathbf{A}^{n}\,\right)_{J}.$
(30)

These difference equations furnish a canonical symplectic PIC method for the
Vlasov-Maxwell equations.

As discussed above, symplectic algorithms for a Hamiltonian system with non-
separable momentum and position dependence are implicit in general. This is
indeed the case for the difference equations (27)-(30), because the right-hand
sides of Eqs. (27)-(30) depend on values of the $(n+1)$-th time-step. However,
they are semi-explicit, because Eqs. (27), (Canonical symplectic particle-in-
cell method for long-term large-scale simulations of the Vlasov-Maxwell
system), and (30) are explicit for $\mathbf{X}_{i}^{n+1}$,
$\mathbf{A}_{J}^{n+1}$, and $\mathbf{Y}_{J}^{n+1},$ respectively. Another good
property of the system is that the only implicit equation (Canonical
symplectic particle-in-cell method for long-term large-scale simulations of
the Vlasov-Maxwell system) is linear in terms of $\mathbf{P}_{i}^{n+1}$, and
it is only implicit for each particle, i.e., Eq. (Canonical symplectic
particle-in-cell method for long-term large-scale simulations of the Vlasov-
Maxwell system) does not couple $\mathbf{P}_{i}^{n+1}$ and
$\mathbf{P}_{k}^{n+1}$ when $i\neq k.$ Therefore, the system can be solved
without root searching iterations as follows. We first solve the linear
equation (Canonical symplectic particle-in-cell method for long-term large-
scale simulations of the Vlasov-Maxwell system) for $\mathbf{P}_{i}^{n+1}$ for
every index $i$ separately, which amounts to inverting a $3\times 3$ matrix
for every $i$. Then $\mathbf{X}_{i}^{n+1}$ and $\mathbf{Y}_{J}^{n+1}$ are
advanced explicitly according to Eqs. (27) and (30), and the last step is to
advance $\mathbf{A}_{J}^{n+1}$, also explicitly, according to Eq. (Canonical
symplectic particle-in-cell method for long-term large-scale simulations of
the Vlasov-Maxwell system).

The preservation of the symplectic structure exerts $D(2D-1)$ constraints on
the numerical solution. Because $D$ is a large number, preservation of
symplectic structure is a very strong constraint and significantly reduces the
errors of numerical solutions. We can also appreciate this advantage from the
viewpoint of symplectic capacity, which is defined on any open set of the
phase space. Symplectic maps preserves symplectic capacity Hofer94 , and in
principle there are infinite constraints that symplectic algorithms can
satisfy as the continuous systems do.

We now apply this canonical symplectic PIC scheme to simulate the nonlinear
Landau damping process. This study also serves as a test of the algorithm.
Previously, similar study and test have been performed for other algorithms,
e.g., the Eulerian algorithms for the Vlasov-Poisson system Zhou01 ; Kraus14 .
The ions are treated as a uniform positively charged background, and the
dynamics of electrons are simulated. The electron density is
$n_{e}=1.149\times 10^{16}/m^{3}$, and the thermal velocity of electrons is
$v_{T}=0.1c$, where $c$ is the light velocity in vacuum. The three-dimensional
computational region is divided into $896\times 1\times 1$ cubic cells. The
size of grid is chosen to be $\Delta x=2.4355\times 10^{-4}m$, the time step
$\Delta t=\Delta x/2c$. The interpolation function is chosen to be 8th order,
i.e.,

$\displaystyle W(\mathbf{x})$ $\displaystyle=$ $\displaystyle W_{1}(x/\Delta
x)W_{1}(y/\Delta x)W_{1}(z/\Delta x)~{},$ (31) $\displaystyle W_{1}(q)$
$\displaystyle=$ $\displaystyle\left\\{\begin{array}[]{lc}0,&q>2~{},\\\
\frac{15}{1024}q^{8}-\frac{15}{128}q^{7}+\frac{49}{128}q^{6}-\frac{21}{32}q^{5}+\frac{35}{64}q^{4}-q+1,&1<q\leq
2~{},\\\
-\frac{15}{1024}q^{8}-\frac{15}{128}q^{7}+\frac{7}{16}q^{6}-\frac{21}{32}q^{5}+\frac{175}{256}q^{4}-\frac{105}{128}q^{2}+\frac{337}{512},&0<q\leq
1~{},\\\
-\frac{15}{1024}q^{8}+\frac{15}{128}q^{7}+\frac{7}{16}q^{6}+\frac{21}{32}q^{5}+\frac{175}{256}q^{4}-\frac{105}{128}q^{2}+\frac{337}{512},&-1<q\leq
0~{},\\\
\frac{15}{1024}q^{8}+\frac{15}{128}q^{7}+\frac{49}{128}q^{6}+\frac{21}{32}q^{5}+\frac{35}{64}q^{4}+q+1,&-2<q\leq-1~{},\\\
0,&q<-2~{}.\end{array}\right.$ (38)

It can be proved that the kernel function $W$ is 3rd order continuous in the
whole space. According to our performance benchmarks, this 8th order kernel is
about 30% more computationally costly than a 2nd order kernel, which is
acceptable since a higher order continuous kernel gives more numerical
fidelity. Initially, $10^{5}$ sampling points of electrons are distributed in
each cell. The total number of particles is $N=8.96\times 10^{7}$, and the
number of degrees of freedom is $D=2.69\times 10^{8}.$ The initial electrical
field perturbation is $\mathbf{E}_{1}=E_{1}\cos\left(kx\right)\mathbf{e}_{x},$
where the wave number is $k=2\pi/224\Delta x$, and the amplitude of the
perturbation electric field is $E_{1}=9.103\times 10^{4}V/m.$ The simulations
are performed for 80000 time steps, during which a complete picture of the
nonlinear Landau damping is revealed. As expected for symplectic algorithms,
the numerical error of energy does not increase with time and is bounded
within 1% for all time. The theoretical damping rate calculated from the
dispersion relation is $\omega_{i}=-1.3926\times 10^{9}/s$, and the
theoretical real frequency is $\omega_{r}=9.116\times 10^{9}/s$. In Fig. 1,
the slope of the green line is the theoretical damping rate, and the blue
curve is the evolution of the electrical field observed in the simulation. It
is evident that the simulation and theory agree perfectly. After
$t=30/\omega_{r}$, the energy of the wave drops below the level of numerical
noise, and the damping process stops. The evolution of the electron
distribution function is plotted in Fig. 2, which clearly demonstrates the
mechanism of phase mixing in velocity space. We observe in Fig. 2 that the
wave-number in velocity space increases with time, which results in a decrease
in density perturbation and thus attenuation of the electrical field. More
importantly, this mechanism of phase mixing is the dominant physics for the
entire nonlinear evolution of the Landau damping, as proved by Mouhot and
Villani Mouhot11 ; Villani14  recently for the electrostatic Vlasov-Poisson
system. In addition, our simulation is electromagnetic and it shows that this
physical picture of nonlinear Landau damping is also valid for the
electromagnetic Vlasov-Maxwell system, as Villani conjectured Villani14-i . It
is the long-term accuracy and fidelity of the canonical symplectic PIC
algorithm that enables the confirmation of Mouhot and Villani’s theory and
conjecture over several orders of magnitude.

Figure 1: Perturbed electrical field as a function of time. The slope of the
green line is the theoretical damping rate. Figure 2: Electron distribution
function in velocity space as a function of time. Different colors denote the
amplitude of the perturbation. The mechanism of phase mixing in velocity space
is clearly demonstrated. The wave-number in velocity space increases with
time, which results in a decrease in density perturbation and thus attenuation
of the electrical field.

Even though Mouhot and Villani rigorously proved that when the initial
perturbation amplitude is small enough, the electrical perturbation will decay
to zero, plasma physicists have known this fact since 1960s Armstrong67 ;
Canosa74 . It has also been known that when the initial electrical
perturbation is large enough, the perturbation will bounce back after initial
phase of damping Armstrong67 ; Canosa74 . One such case simulated is plotted
in Fig. 3. For this case, the amplitude of the initial electrical field is
0.494MV/m and the wave number is $k=2\pi/272\Delta x$. Figure 4 shows the
bounce-time as a function of the initial amplitude of the electrical field
obtained in simulations. The physics demonstrated in our PIC simulations
agrees with that obtained from Eulerian solvers Zhou01 , except that our
simulations are carried out for the full Vlasov-Maxwell system.

Figure 3: In nonlinear Landau damping, the perturbation will bounce back after
initial phase of damping if the initial perturbation is large enough. The
amplitude of the initial electric field in this case is 0.494MV/m.  Figure 4:
Relationship between the initial amplitude $E_{1}$ and the bounce time $t_{b}$
in nonlinear Landau damping.

In the second application, the 1D nonlinear mode conversion of extraordinary
waves to electron Bernstein waves (X-B mode conversion) in a inhomogeneous hot
plasma is simulated for a long time. The plasma density profile in the
simulation is

$n_{e}(x)=n_{0}\left\\{\begin{array}[]{lc}\exp\left[-\left(\frac{x/\Delta
x-n_{r}-320}{0.4n_{r}}\right)^{2}\right]~{},&0<\frac{x}{\Delta
x}\leq\left(n_{r}+320\right)~{},\\\
1~{},&\left(n_{r}+320\right)<\frac{x}{\Delta x}\leq
1300~{},\end{array}\right.$ (39)

where $n_{0}=2.3\times 10^{19}~{}\mathrm{m}^{-3},$ $n_{r}=380,$ and $\Delta
x=2.773\times 10^{-5}$m is the grid size. The thickness of the plasma boundary
is $n_{r}\Delta x$. The electron temperature is $T_{e}=57.6\mathrm{{eV}}$, and
constant external magnetic field is $\mathbf{B}=B_{0}\mathbf{e}_{z}$ with
$B_{0}=0.6\mathrm{{T}}$. The simulation domain is a $1584\times 1\times 1$
cubic mesh. At both boundaries in the $x$-direction, the Mur’s absorbing
condition are used, and periodic boundary conditions are adopted in the $y$\-
and $z$-directions. The time step is chosen to be $\Delta t=\Delta x/2c$. At
the left boundary, a source is placed to excite an electromagnetic
perturbation at $\omega=0.0145/\Delta t$ with
$\mathbf{E}_{1}=E_{1}\mathbf{e}_{y}$ and $E_{1}=900\mathrm{{kV}}$. As
illustrated in Fig. 5, the extraordinary wave excited at the left boundary
first propagates to the region of cutoff-resonance-cutoff near $x/\Delta
x=500.$ The wave is then partially reflected back, and partially converted to
electron Bernstein waves Xiao15 . The reflectivity evolution is plotted in
Fig. 6. Nonlinear excitations and self-interactions of the Bernstein modes
dominate the long time dynamics of the mode conversion process. As a
consequence, the reflectivity of the incident wave evolves nonlinearly as
shown in Fig. 6. For this set of chosen parameters, the incident wave is
completely reflected at the later time of the process.

Figure 5: The space-time dependence of $E_{y}(t,x)$ during the nonlinear X-B
mode conversion. Figure 6: The evolution of reflectivity during the nonlinear
X-B mode conversion.

In conclusion, we have developed a canonical symplectic particle-in-cell
simulation method for the Vlasov-Maxwell system by discretizing its canonical
Poisson bracket. In phase space, the distribution function is discretized by
the Klimontovich representation using Lagrangian markers, and the
electromagnetic field is discretized point-wise on a spatial grid. The
resulting canonical Hamiltonian system with a large number of degrees of
freedom is integrated by the symplectic Euler method, whose difference
equations can be solved inexpensively by inverting a $3\times 3$ matrix
locally for every particle. Implicit root searching and global matrix
inversion are avoided entirely. This technique makes large-scale applications
of the developed canonical symplectic PIC method possible. To suppress
numerical noise caused by the coarse sampling, smoothing functions for
sampling points can also be conveniently implemented in the canonical
symplectic PIC algorithm. By incorporating the smoothing functions into the
Hamiltonian functional before the discretization, we are able to rein in all
the benefits of smoothing functions without destroying the canonical
symplectic structure. Progress in this and other directions will be reported
in future publications.

###### Acknowledgements.

This research is supported by ITER-China Program (2015GB111003, 2014GB124005,
2013GB111000), JSPS-NRF-NSFC A3 Foresight Program in the field of Plasma
Physics (NSFC-11261140328), the CAS Program for Interdisciplinary
Collaboration Team, the Geo-Algorithmic Plasma Simulator (GAPS) project, and
the the U.S. Department of Energy (DE-AC02-09CH11466).

## References

  * (1) C. K. Birdsall and A. B. Langdon.  Plasma Physics via Computer Simulation.  Adam Hilger, Bristol.
  * (2) R. W. Hockney and J. W. Eastwood.  Computer Simulation Using Particles.  Institute of Physics Publishing, Bristol, 1988.
  * (3) J. Boris.  In Proceedings of the Fourth Conference on Numerical Simulation of Plasmas, page 3. Naval Research Laboratory, Washington D. C., 1970.
  * (4) J. P. Boris.  Dynamic stabilization of the imploding shell rayleigh-taylor instability.  Comments Plasma Phys. Controlled Fusion, 3:1, 1977.
  * (5) J. Villasenor and O. Buneman.  Rigorous charge conservation for local electromagnetic field solvers.  Computer Physics Communications, 69:306, 1992.
  * (6) R. D. Ruth.  A canonical integration technique.  IEEE Trans. Nucl. Sci, 30:2669, 1983.
  * (7) K. Feng.  On difference schemes and sympletic geometry.  In K. Feng, editor, the Proceedings of 1984 Beijing Symposium on Differential Geometry and Differential Equations, pages 42–58. Science Press, 1985.
  * (8) K. Feng.  Difference schemes for Hamiltonian formalism and symplectic geometry.  J. Comput. Maths., 4:279–289, 1986.
  * (9) Kang Feng and Mengzhao Qin.  Symplectic Geometric Algorithms for Hamiltonian Systems.  Springer-Verlag, 2010.
  * (10) E. Forest and R. D. Ruth.  4th-order symplectic integration.  Physica D, 43:105–117, 1990.
  * (11) P. J. Channell and C. Scovel.  Symplectic integration of Hamiltonian systems.  Nonlinearity, 3:231–259, 1990.
  * (12) J. Candy and W. Rozmus.  A symplectic integration algorithm for separable Hamiltonian functions.  Journal of Computational Physics, 92:230–256, 1991.
  * (13) J. E. Marsden and M. West.  Dicrete mechanics and variational integrators.  Acta Numerica, 10:357–514, 2001.
  * (14) E. Hairer, C. Lubich, and G. Wanner.  Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations.  Springer, New York, 2002.
  * (15) H. Qin and X. Guan.  A variational symplectic integrator for the guiding center motion of charged particles for long time simulations in general magnetic fields.  Physical Review Letters, 100:035006, 2008.
  * (16) H. Qin, X. Guan, and W. M. Tang.  Variational symplectic algorithm for guiding center dynamics and its application in tokamak geometry.  Physics of Plasmas, 16:042510, 2009.
  * (17) H. Qin, S. Zhang, J. Xiao, J. Liu, Y. Sun, and W. M. Tang.  Why is boris algorithm so good?  Physics of Plasmas, 20:084503, 2013.
  * (18) R. Zhang, J. Liu, H. Qin, Y. Wang, Y. He, and Y. Sun.  Volume-preserving algorithm for secular relativistic dynamics of charged particles.  Physics of Plasmas, 22:044501, 2015.
  * (19) Y. He, Y. Sun, J. Liu, and H. Qin.  Volume-preserving algorithms for charged particle dynamics.  Journal of Computational Physics, 281:135, 2015.
  * (20) J. Squire, H. Qin, and W. M. Tang.  Geometric integration of the Vlasov-Maxwell system with a variational particle-in-cell scheme.  Physics of Plasmas, 19:084501, 2012.
  * (21) J. Y. Xiao, J. Liu, H. Qin, and Z. Yu.  A variational multi-symplectic particle-in-cell algorithm with smoothing functions for the Vlasov-Maxwell system.  Physics of Plasmas, 20:102517, 2013.
  * (22) M. Kraus.  Variational integrators in plasma physics.  PhD thesis, Technical University of Munich, 2014.
  * (23) B. A. Shadwick, A. B. Stamm, and E. G. Evstatiev.  Variational formulation of macro-particle plasma simulation algorithms.  Physics of Plasmas, 21:055708, 2014.
  * (24) J. Xiao, J. Liu, H. Qin, Z. Yu, and N. Xiang.  Variational symplectic particle-in-cel simulation of nonlinear mode conversion from extraordinary waves to Bernstein waves.  Physics of Plasmas, page in press, 2015.
  * (25) E. Evstatiev and B. Shadwick.  Variational formulation of particle algorithms for kinetic plasma simulations.  J. Comput. Phys., 245:376, 2013.
  * (26) J. E. Marsden and A. Weinstein.  The Hamiltonian structure of the Maxwell-Vlasov equations.  Physica, 4D:394, 1982.
  * (27) C. Mouhot and C. Villani.  On landau damping.  Acta Mathematica, 207:29, 2011.
  * (28) C. Villani.  Particle systems and nonlinear landau damping.  Physics of Plasmas, 21:030901, 2014.
  * (29) C. Villani.  Plasmatalks: Ron Davidson interviews Cedric Villani, 2014.
  * (30) P. J. Morrison.  The Maxwell-Vlasov equations as a continuous Hamiltonian system.  Physics Letters, 80A:383, 1980.
  * (31) A. Weinstein and P. J. Morrison.  Comments on: The Maxwell-Vlasov equaitons as a continuous Hamiltonian system.  Physics Letters, 86A:235, 1981.
  * (32) P. J. Morrison.  Hamiltonian field description of two dimensional vortex fluids and guiding center plasmas.  Technical Report PPPL-1783, Princeton Plasma Physics Laboratory, 1981\.
  * (33) P. J. Morrison.  Hamiltonian field description of the one-dimensional poisson-vlasov equations.  Technical Report PPPL-1788, Princeton Plasma Physics Laboratory, 1981\.
  * (34) H. Hofer and E. Zehnder.  Symplectic invariants and Hamiltonian dynamics.  Birkhauser Verlag, 1994.
  * (35) T. Zhou, Y. Guo, and C. W. Shu.  Numerical study on Landau damping.  Phsica D, 157:322, 2001.
  * (36) T. P. Armstrong.  Numerical studies of the nonlinear Vlasov equation.  Physics of Fluids, 6:1269, 1967.
  * (37) J. Canosa and J. Gazdog.  Threshold conditions for electron trapping by nonlinear waves.  Physics of Fluids, 17:2030, 1974.


# An Evaluation of ChatGPT-4’s
Qualitative Spatial Reasoning Capabilities in RCC-8

A G Cohn

###### Abstract

Qualitative Spatial Reasoning (QSR) is well explored area of Commonsense
Reasoning and has multiple applications ranging from Geographical Information
Systems to Robotics and Computer Vision. Recently many claims have been made
for the capabilities of Large Language Models (LLMs). In this paper we
investigate the extent to which one particular LLM can perform classical
qualitative spatial reasoning tasks on the mereotopological calculus, RCC-8.

## Introduction

Qualitative Spatial Reasoning (QSR111We may use QSR as shorthand for both
Qualitative Spatial Reasoning and Qualitative Spatial Representation; context
should usually make clear which is intended.) (Cohn and Renz 2008; Chen et al.
2015; Cohn and Hazarika 2001) is a well developed field which is concerned
with the representation of qualitative spatial information and reasoning with
it. In natural language, spatial information is usually represented
qualitatively (using prepositions such as _on, in, left of, part of, under,
touching, …)_ and many calculi have been developed to represent such
information. There are calculi for mereological relations (such as RCC-5
(Jonsson and Drakengren 1997)), mereotopological relations (such as RCC-8
(Randell, Cui, and Cohn 1992; Cohn et al. 1997)), directions (such as OPRA
(Moratz 2006)), size (Gerevini and Renz 2002) for example as well as calculi
combining two different aspects of spatial information, such as the Rectangle
Algebra (Guesgen 1989; Mukerjee and Joe 1990) which can represent both
mereotopological information as well as directional. What is common to all
these calculi is that they consist of a set of _jointly exhaustive and
pairwise disjoint_ (JEPD) _base_ relations. For example, RCC-8 contains eight
JEPD _base_ relations, illustrated in 2D in Fig. 1.

Figure 1: The eight relations of the RCC-8 calculus illustrated in 2D.

_Large Language Models_ (LLMs) (Devlin et al. 2019; Brown et al. 2020), such
as ChatGPT-4 (Roumeliotis and Tselikas 2023) are a recent example of so called
_Foundation Models_ which have been trained on very large textual corpora in
order to generate text in response to a prompt. This is not the place to
survey this burgeoning field, but we note that many claims have been made for
the power and apparent intelligent behaviour that these models can display. In
particular their performance on some benchmarks may lead one to believe that
they possess, at least to some degree, the ability to perform commonsense
reasoning. Spatial reasoning is usually regarded as one core aspect of common
sense so it is natural to ask whether LLMs can reason about qualitative
spatial information. This is the question that we address here.

In earlier work (Cohn and Hernandez-Orallo 2023) we use extended dialogues
with an LLM to try to map the boundaries of spatial commonsense in some LLMs,
addressing a variety spatial challenges, and examining not only the response
given but also the explanation/justification of the response, but did not
specifically focus on existing QSRs, though some questions were asked which do
correspond to particular reasoning steps in an existing QSR. Here we focus on
one specific QSR and ask the question as to what extent an LLM can perform
reasoning in that calculus, and conduct a more exhaustive evaluation, but
looking at the ability to perform compositions between relations and also to
reason about the conceptual neighbourhood diagram of the calculus. Weaknesses
in the reasoning powers of LLMs have previously been noted (e.g. (Cai, Chang,
and Han 2023)) so one might not expect LLMs to perform well in this regard.
But on the other hand, there are a large number of papers about QSR in the
literature and these are likely to have formed part of the training corpus of
an LLM, and thus might facilitate correctly responding to prompts – though the
information concerning the actual reasoning steps are often given in tables
(in particular _composition tables_ – see below) and thus might be hard for
LLM training procedures to process well.

There are now many LLMs in the literature. Some of these are open source and
are explicit about the training corpus; others are closed and give no specific
information about the training, or the precise corpus, such as the GPT family
of LLMs. Nevertheless since we observed previously (Cohn and Hernandez-Orallo
2023) that ChatGPT-4 and GPT4 were the most performant for spatial reasoning,
we use ChatGPT-4 as the LLM with which we perform our experiments. In each of
the experiments below, an initial prompt gave the problem setting and the task
to be performed. Subsequent prompts in the conversation probed about one
specific inference (e.g. one cell in a composition table). Each experiment was
a separate conversation and the version number of ChatGPT-4 was given as “May
24” (presumably, the version of May 24 2023). The conversations can be found
at the following location: https://tinyurl.com/qr23sup.

## Compositional Reasoning

The most researched form of reasoning with Qualitative Spatial Relations is
that of composition: i.e. given a two facts R1($x,y$), and R2($y,z$), then
what relations are possible between $x$ and $z$? In general more than one base
relation is possible, so the answer may take the form of a disjunction of
relations, as can be seen in the composition table for RCC-8 in Fig. 2. Given
the ubiquity of RCC-8 in the QSR literature, in this paper we focus on the
abilities of an LLM to reason with RCC-8.

Figure 2: The RCC-8 Composition Table (Cohn et al. 1997)

.

### Experiment 1: Compositional Reasoning in RCC-8

In this experiment the standard names of the eight relations were given as
well as a sentence stating that the relations were part of the RCC-8 calculus:

> Consider the following set of eight pairwise disjoint and mutually
> exhaustive binary spatial relations. These relations form part of the well
> known RCC-8 qualitative spatial reasoning calculus. DC(x,y) means that x and
> y are disconnected and share no spatial parts. EC(x,y) means that x and y
> touch at a boundary but do not share any interior parts. PO(x,y) means that
> x and y share a spatial part, but neither is part of the other. TPP(x,y)
> means that x is part of y and touches y’s boundary. NTPP(x,y) means that x
> is part of y but does not touch y’s boundary. TPPi(x,y) is the same as
> TPP(y,x). NTPPi(x,y) is the same as NTPP(y,x). Finally, EQ(x,y) means that x
> and y are coincident. I will now ask you a series of questions about these
> relations. There may be more than one possible relation, in which case name
> all of the possible answers. If all eight relations are possible, just say
> ALL.

As an example of a prompt and ChatGPT-4’s response, consider the first one:

> Prompt: _If DC(x,y) and DC(y,z) then what are the possible relationships
> between x and z._ Response:  If DC(x,y) and DC(y,z), then x and z can either
> be disconnected (DC) or externally connected (EC). This is because the
> disconnection between x and y, and y and z, only guarantees that neither x
> nor z is connected to y. It does not preclude the possibility that x and z
> themselves might be disconnected or touching at a boundary. However, the
> other relations like part-of (TPP, NTPP) and their inverses (TPPi, NTPPi),
> partial overlap (PO), and equality (EQ) are not possible because they would
> imply some level of connection between x and z through y, which contradicts
> our initial conditions. So, the possible relationships between x and z are:
> DC(x,z): x and z are disconnected. EC(x,z): x and z are externally
> connected.

Notice that in fact the composition of DC with DC is the union of all eight
RCC-8 relations, so the reasoning is incorrect here.

The results of the 7*7=49 prompts issued (compositions involving EQ were
omitted since they should be trivial) are summarised in Fig. 3. There is a
total of 85 relations correctly predicted (those in black), 197 relations
which are correctly not predicted (those not in the table). Thus there are a
total of 282 correctly predicted (both present and missing) relations, i.e.
71.94% of the total possible (7*7*8 = 392). 11 of the cells (out of 49) are
completely correct. This is clearly much better than chance but far from a
correct composition table. Fig. 4 shows the absolute number and relative
percentages of the 8 relations in the proportion that they are correctly or
incorrectly predicted. As can be seen, DC is never missed as a correct
prediction, and PO is never incorrectly predicted as present. DC has the
highest accuracy which is perhaps not surprising since DC might be regarded as
intuitively easier to reason about (and this fits with the predictions about
human preferred models in the section below). Including both correctly
predicted and correctly missing relationships, all eight relations have
broadly similar accuracies.

Figure 3: The Composition Table for RCC-8 produced by ChatGPT-4. The entry in
each cell uses the following coding: D (DC), E(EC), P(PO), T(TPP), N(NTPP),
t(TPPi), n(NTPPi), Q(EQ). Black means that relation is correctly predicted (85
times), red means that relation is incorrectly predicted (61 times), blue
means that the relation was incorrectly not predicted (49 times).

.

Figure 4: Relation statistics for the Composition Table for RCC-8 produced by
ChatGPT-4. The upper chart shows the absolute number of relations, and the
lower the relative percentage for each relation.

.

In order to test whether the result was influenced by prior knowledge of RCC-8
gained as part of its training, we also performed the same experiment, but
with all the relation names prefixed by an X to disguise the connection to
RCC-8. The prompt was the same as above except for the change of relation
names and the omission of the second sentence. The results are given in Fig. 5
while Fig. 6 shows the absolute number and relative percentages of the 8
relations in the proportion that they are correctly or incorrectly predicted.
As can be seen, DC again is never missed as a correct prediction, and EC is
only missed twice; again PO is never incorrectly predicted as present. As
before, DC, EC, and POhave the highest accuracies, along with EQ, but EQ is
never predicted as present correctly, only incorrectly. The overall average of
correctly predicted relations (present and missing) drops from 71.94% in the
non-anonymous case above to 67.09% so there is some loss of performance though
whether is due to the anonymisation of the relations or the stochasticity of
ChatGPT-4 is not clear.

Figure 5: The Composition Table for RCC-8 produced by ChatGPT-4 using
anonymised relations. The entry in each cell uses the following coding: D
(DC), E (EC), P (PO), T (TPP), N (NTPP), t (TPPi), n (NTPPi), Q (EQ). Black
means that relation is correctly predicted (95 times), red means that relation
is incorrectly predicted (57 times), blue means that the relation was
incorrectly not predicted (72 times).

.

Figure 6: Relation statistics for the Composition Table for RCC-8 with
anonymised relations produced by ChatGPT-4. The upper chart shows the absolute
number of relations, and the lower the relative percentage for each relation.

.

### Experiment 2: Preferred Compositions in RCC-8

As noted above, in general a composition of two relations will yield more than
one possible base relations, but it turns out that humans tend to have a
“preferred” relation. For example, Ragni et al (2007) report on experiments
performed on native German speakers and native Mongolian speakers for RCC-8.
In their experiments the relations were described, but the human subjects were
not allowed to draw possible configurations, so the setting is essentially
equivalent to an LLM setting.

Given that humans may struggle to see all the possible relations222The fact
that some humans may struggle to compute the composition table does not stop
it being a valid question to see if an LLM can determine the correct entries.,
determining whether there is agreement about the most preferred is good
question to ask. It turns out that there is good agreement in general across
and within the two cultures, with the the percentage of people agreeing with
the same preferred relation ranging from 30% to 87.5% (a random choice would
yield 12.5% on average since there are eight relations to choose from). (They
did not query cases where the composition yields a unique relation, nor did
they consider EQ as one of the two relations as this should be a trivial
task.) This agreement is perhaps surprising since the two languages are
linguistically very different. Ragni et al (2007) do report some differences
though – for example although both language speakers preferred DC whenever it
was consistent, Mongolians preferred PO over NTPPi whereas for Germans the
converse was true. Both cultures only chose EQ when composing a relation with
its inverse (e.g. TPP with TPPi).

The theory of _preferred mental models_ (Knauff, Rauh, and Schlieder 1995)
states that people construct the simplest (computationally cheapest) model
consistent with the premises. Their experiments showed that humans prefer
models with the smallest overlapping complexity which explains the preference
for DC noted above.

Given the difficulty reported in Experiment 1 in ChatGPT-4 correctly inferring
all possible relations in a composition, asking this more specific query seems
a natural task to try; this gives some measure of whether the cognitive
preferences of humans correspond at all to the choices made by ChatGPT-4.

The initial prompt given to ChatGPT-4 was identical to the initial prompt
above when querying about composition. In its responses, sometimes ChatGPT-4
gave a clear preference for one relation, for example:

> Prompt: _If DC(x,y) and DC(y,z) then what is your preferred relationship
> between x and z?_ Response: _The relations DC(x,y) and DC(y,z) mean that
> region x is disconnected from region y and region y is disconnected from
> region z. These relationships don’t provide any information about the
> spatial relationship between x and z. x and z could be disconnected, they
> could touch at a boundary, or they could partially overlap, among other
> possibilities. However, in absence of any other information and considering
> spatial continuity, I would say the most likely relation is DC(x, z),
> assuming x and z are also disconnected. Please note that this is just a
> statistical assumption; it’s possible for x and z to have any of the other
> relations defined in RCC-8._

In other cases ChatGPT-4 justified its choice by saying it was the “cautious”
choice, or the “safest choice”. The humans in the (Ragni, Tseden, and Knauff
2007) experiment were not asked to justify their choices, but the claim made
by the authors of that paper noted above is that choices were based on
computational/simplicity considerations – which differs from the reasons
claimed here by ChatGPT-4.

The initial prompt to ChatGPT-4 given was:

> Consider the following set of eight pairwise disjoint and mutually
> exhaustive binary spatial relations, between two one piece spatial 2D
> regions. These relations form part of the well known RCC-8 qualitative
> spatial reasoning calculus. DC(x,y) means that x and y are disconnected and
> share no spatial parts. EC(x,y) means that x and y touch at a boundary but
> do not share any interior parts. PO(x,y) means that x and y share a spatial
> part, but neither is part of the other. TPP(x,y) means that x is part of y
> and touches y’s boundary. NTPP(x,y) means that x is part of y but does not
> touch y’s boundary. TPPi (x,y) is the same as TPP(y,x). NTPPi(x,y) is the
> same as NTPP(y,x). Finally, EQ(x,y) means that x and y are coincident. I
> will now ask you a series of questions about these relations. It is possible
> that in some cases there could be more than one relation that might hold; in
> these cases just give me your single most preferred relationship.

Figure 7: The table of preferred compositions for RCC-8 produced by ChatGPT-4.
Green text means it agreed with the human most likely model. Orange text means
it agreed with the human most likely model but failed to notify that it was
the only possible model. Red means it chose an impossible relation. Blue means
it chose a possible but not preferred model. A green background means that it
agreed with preferred relation of one of the language groups but not the
overall average.

.

As can be seen in Fig. 7, ChatGPT-4 only agreed with the average human on
20/49 (40.82%) though in a further three cases it agreed with one of the
language groups (twice Mongolian, once German) giving a total of 23/49
(46.93%). In seven cases it actually predicted an impossible relation as its
preferred composition333Sometimes the humans in the (Ragni, Tseden, and Knauff
2007) experiment also predicted an impossible relation as can be seen in Table
2 of their publication., and in all cases where the composition is unique it
failed to note that. In all the other cases, 18/49 (36.63%) its preferred
relation was possible but not preferred by either language group or overall.

As can be seen, sometimes ChatGPT-4 made a logically inconsistent prediction,
for example with EC/NTPPi; looking at the ChatGPT-4’s response in this case we
can see that this is because it seems to misunderstand the meaning of NTPPi:

> _If EC(x,y) and NTPPi(y,z), this means that region x touches the boundary of
> region y but does not overlap with it, and region y is a part of region z
> but does not touch z’s boundary, with the perspective of y being contained
> in z. In this scenario, it’s likely that x and z are either disconnected
> (DC) or touch at the boundary (EC). Since y, which is part of z but does not
> touch z’s boundary, is touching x at the boundary, there’s a possibility for
> indirect spatial contact between x and z. Without additional context, it’s
> safer to assume that x and z touch at the boundary (EC). However, this is
> just an assumption, and the actual relationship between x and z could be
> different depending on the specific spatial arrangement._

In a different case, TPPi/TPPi, the initial understanding of the relations is
correct inferring that $z$ is part of $x$, but then ChatGPT-4 makes a wrong
step in the final part of its reasoning, and claims that TPP($x,y)$:

> _If TPPi(x,y) and TPPi(y,z), this means that region y is a part of region x
> and touches x’s boundary, and region z is a part of region y and touches y’s
> boundary. Given these conditions, it’s likely that z is a part of x and
> touches x’s boundary (TPP). This is because z, which is a part of y and
> touches y’s boundary, is also indirectly a part of x through y, which is a
> part of x and touches x’s boundary. Therefore, my preferred relationship
> between x and z would be TPP(x,z), assuming that z is a part of x and
> touches x’s boundary. Please note that this is an assumption and the actual
> relationship between x and z could be different depending on the specific
> spatial arrangement._

ChatGPT-4 makes exactly the same error with TPPi/NTPPi. For the wrong NTPPi/EC
case, the the description of its reasoning contains a more fundamental error:

> _…This is because y, which is a part of x and does not touch x’s boundary,
> is also only touching z’s boundary, so x and z would not have any overlap…_

A similar more fundamental error is made in the cases of NTPPi / TPP, whilst
in the final two wrong cases, ChatGPT-4 misunderstands the second relation as
its inverse.

In order to test whether the result was influenced by prior knowledge of RCC-8
gained as part of its training, we also performed the same experiment, but
with all the relation names prefixed by an X to disguise the connection to
RCC-8. The prompt was the same as above except for the change of relation
names and the omission of the second sentence. The results are given in Fig.
8. ChatGPT-4 only agreed with the average human on 24/49 (48.99%) though in a
further five cases it agreed with one of the language groups (twice times
Mongolian, twice German), and in one case the Mongolian preference was equal
between the overall choice and the choice ChatGPT-4 made giving a an agreement
with the German speakers of 26/49 (53.06%) and 27/49 (55.1%) with the
Mongolian speakers. In three cases (6.12%) it actually predicted an impossible
relation as its preferred composition (this also happened sometimes in the
human experiment), and in two cases (out of 11) where t he composition is
unique it failed to note that. In all the other cases, 20/49 (60.82%) its
preferred relation was possible but not preferred by either language group or
overall. The overall agreement with humans in this anonymised case is actually
significantly higher than the non-anonymised experiment above; it is not clear
if there is a reason for this and if so what.

Figure 8: The table of preferred compositions for RCC-8 produced by ChatGPT-4
with anonymised relations. The relation names are coloured in the same way as
in Table 7.

.

## Experiment 3: Spatial Continuity

Continuity networks were introduced by Randell and Cohn(1989) to represent the
set of possible ‘next’ relations that might obtain, assuming that motion is
continuous and any transformations in the shape and/or size of an object are
also continuous. Subsequently these have been termed “conceptual
neighbourhoods” (Freksa 1992). Fortunately, the abbreviation for both of these
terms is CN, which we shall use henceforth. The CN for RCC-8 is depicted in
Fig. 1. The initial prompt was the same as for the composition table prompt
above, but with two additional sentences at the end: _In all these questions,
assume that all motion is continuous, and any deformations of the spatial
regions are also continuous. Also please describe the motion or changes of
size or shape that would be required for the change of relation to take
place_. Each subsequent prompt also emphasises the need for the relation to be
an immediate next relation. As an example prompt/response pair consider the EC
case:

> Prompt: If EC(a,b) holds, then what are the next possible relations that
> might hold as the two regions translate or change their shape or size? I am
> interested in the immediate next relation, so no intermediate relation must
> hold. Response: If EC(a,b) holds, which means that the two regions a and b
> are just touching at a boundary but do not share any interior parts, then
> the immediate next possible relations as the two regions translate or change
> their shape or size could be either DC(a,b) or PO(a,b). DC(a,b): This would
> occur if the regions move away from each other and are no longer in contact
> at all, becoming fully disconnected. PO(a,b): This would occur if one region
> begins to overlap with the other, but neither region is completely contained
> within the other. They would then be partially overlapping, sharing some but
> not all of their spatial parts. The transitions to these relations would
> involve either moving the regions apart or moving them further into each
> other. In the case of transition to PO(a,b), there could also be some
> deformation of the regions if they change shape to allow for the partial
> overlap.

An example of a faulty prediction is from PO:

> _If PO(a,b) holds, indicating that regions a and b partially overlap but
> neither is fully contained within the other, the immediate next relations
> could be EC(a,b), DC(a,b), TPP(a,b), TPPi(a,b), or EQ(a,b). EC(a,b): This
> would occur if the overlapping parts of a and b are reduced to just boundary
> parts, meaning the regions are touching but no longer overlapping. DC(a,b):
> This would occur if regions a and b are moved or deformed such that they no
> longer overlap or touch at all. TPP(a,b) or TPPi(a,b): This would occur if
> one region becomes fully contained within the other, either by moving into
> it or shrinking. TPP(a,b) means a is part of b and touches b’s boundary.
> TPPi(a,b) means b is part of a and touches a’s boundary. EQ(a,b): This could
> occur if the regions are reshaped or resized such that they coincide
> exactly, becoming the same region in space. The transitions to these
> relations would involve either moving the regions, resizing them, or
> reshaping them in some way._

Here, ChatGPT-4 has mistakenly said that DC is a next relation, failing to
notice that it has to pass through EC on the way.

As can be seen, the predictions are mostly correct with just five errors:
three where a link is falsely predicted and two where a link is missing. The
two missing links are both from EQ (to NTPP and to NTPPi). Interestingly the
links in the reverse direction are correctly predicted, so ChatGPT-4 is not
able to reason that a link in one direction implies the link in the other
direction should also be present. Similarly whilst the absence of a link from
DC to PO is correctly predicted, the inverse case is not, and is also the case
for PO to/from NTPP and NTPPi.

Figure 9: The Continuity Table for RCC8 produced by ChatGPT-4. An ‘x’ means
that the relation in that column is predicted as an immediate neighbour of the
relation in that row. An empty box means that the relation is not predicted as
an immediate neighbour. Green means that the prediction was correct and red
that it was incorrect. The leading diagonal is white since a relation is not a
next relation of itself.

.

In order to test whether the result was influenced by prior knowledge of RCC-8
gained as part of its training, we also performed the same experiment, but
with all the relation names prefixed by an X to disguise the connection to
RCC-8. The prompt was the same as above except for the change of relation
names and the omission of the second sentence. The results are given in Fig.
10. There are 3 incorrectly predicted links, 3 missing links, 19 correctly
predicted links and 31 correct missing links, giving an accuracy of 50/56
(89.2%). This is slightly worse than the case above. There is one more missing
link but the missing links are all different in the two cases. Although there
are the same number of wrong links, only one of these is in common (PO to DC).
Overall the results are broadly similar and may be due to the stochastic
nature of ChatGPT-4’s responses, suggesting that either the disguise was not
very effective, or that prior training did not really affect the response and
it was able to reason from ‘first principles’ (if not always correctly) in
response to each prompt.

Figure 10: The Continuity Table for RCC-8 produced by ChatGPT-4 using
disguised relation names. The meaning of the colouring is the same as in Fig.
9.

.

## Concluding Remarks and Future Work

This investigation has supported the widely-held view that LLMs can struggle
to do reasoning tasks444 Bender et al (2021) have observed that LLMs might be
regarded just as “stochastic parrots” and thus it is not suprising that
precise, logically correct deductive reasoning is challenging for an LLM.. In
the case of Experiment 1, in which ChatGPT-4 was asked to compute the entire
composition table for RCC-8, this is a non trivial task even for humans, so it
is perhaps not surprising that ChatGPT-4 did not achieve 100% accuracy – the
scores of 71.94% (and 67.09% for the anonymised relatins) are clearly much
better than chance and do suggest a reasonable facility to perform such
computations. A detailed analysis of the actual conversations in the
supplementary material shows that sometimes ChatGPT-4 does appear able to do
some interesting (qualitative) spatial reasoning, but often fails, sometimes
making elementary mistakes. It also shows inconsistency in being able to
reason correctly about a relation but not its inverse. It also sometimes
confuses a relation with its inverse. It is possible that fine tuning,
explicit chain-of-thought prompting, or more carefully engineered prompts
might improve performance; however, given the stochastic nature of LLMs it
seems unlikely that the results would be as good as logical reasoning (the
experiment on preferred relations is of course not strictly a logical
reasoning exercise, except for the requirement not to predict spatially
impossible relations).

There are a variety of avenues for further work which present themselves.
Other calculi could be experimented with – for example the coarser calculus
RCC-5, or calculi for reasoning about direction or size (Cohn and Renz 2008).
Other LLMs could be evaluated – though since new LLMs and new LLM versions are
continually being released, this is a challenge with no definite stopping
point. Tracking the change in performance of a particular LLM across releases
would also be of interest – though in the case of closed LLMs such as
ChatGPT-4 where the owners have the right to harvest user conversations and
use them for future training, it will not be clear if any improvement is the
result of leakage from the previous conversation or more general performance
improvement555However, note that no feedback was given to ChatGPT-4 as to
whether the proffered response was correct or not.. It has already been
observed (Cohn and Hernandez-Orallo 2023) that different LLMs have different
strengths – determining which LLMs are better at which spatial reasoning tasks
would also be worth of future investigation. The overall conclusion that LLMs
in general struggle with more complex spatial reasoning tasks is likely to
remain the case, at least for the foreseeable future. In the API version of
GPT, different temperatures could be tried, and multiple runs with averages
computed. Different prompts and prompting strategies could be tried, though
arguably since QSR has always been viewed as a form of commonsense reasoning,
it should not be necessary to devise specific prompts to elicit commonsense
behaviour.

It is not clear how successful the anonymisation was – in one case I mistyped
an X relation and it was able to suggest the intended relation name,
suggesting that it has the ability to dissect relation names; thus more
sophisticated anonymisation might be tried. In earlier work (Cohn and
Hernandez-Orallo 2023) we had already done some limited experimentation asking
an LLM to reason about spatial relations in a real world context rather than
the purely abstract setting used in the experiments in this paper – it would
be interesting to conduct more extensive tests LLMs doing compositional
reasoning in a more realistic setting, and similarly for the continuity
experiment.

Experiment 2 above already investigated how LLM performance compared to human
performance to a limited extent but further investigation would be worthwhile,
including a head-to-head comparison rather than simply taking a result from
the literature originally intended to investigate a different question.
Another interesting avenue for further work will be to explore the use of
multimodal FMs – when humans perform spatial reasoning tasks including the
challenge of building a composition table, it is natural to use pencil and
paper to sketch diagrams and possible scenarios – investigating whether a
multi-modal FM with such abilities (including the ability to analyse its own
drawings) would be of great interest to the spatial reasoning community.

As mentioned above, another possible avenue of research is to investigate
different prompting strategies, including k-shot (Dang et al. 2022), chain-of-
thought(Wei et al. 2022) and tree-of-thought(Yao et al. 2023) strategies. Not
doing so was deliberate in this paper as I was interested in exploring in how
the “vanilla” LLM would perform. Whilst for specific downstream tasks, fine-
tuning or employing specific prompting strategies may reasonable, there is an
argument to be made that for commonsense reasoning, this is not a reasonable
strategy since the task is a general one rather than a specific downstream
task.

## Data statement

All the conversations with ChatGPT-4 that support the summary tables in this
paper can be found at http://tinyurl.com/qr23sup.

## Acknowledgments

This work was supported by: The Alan Turing Institute; the Economic and Social
Research Council (ESRC) under grant ES/W003473/1; the Turing’s Defence and
Security programme through a partnership with the UK government in accordance
with the framework agreement between GCHQ and The Alan Turing Institute.

## References

  * Bender et al. (2021) Bender, E. M.; Gebru, T.; McMillan-Major, A.; and Shmitchell, S. 2021.  On the Dangers of Stochastic Parrots: Can Language Models Be Too Big?  In _Proceedings of the 2021 ACM conference on fairness, accountability, and transparency_ , 610–623.
  * Brown et al. (2020) Brown, T.; Mann, B.; Ryder, N.; Subbiah, M.; Kaplan, J. D.; Dhariwal, P.; Neelakantan, A.; Shyam, P.; Sastry, G.; Askell, A.; et al. 2020.  Language models are few-shot learners.  _Advances in Neural Information Processing Systems_ , 33: 1877–1901.
  * Cai, Chang, and Han (2023) Cai, Z.; Chang, B.; and Han, W. 2023.  Human-in-the-Loop through Chain-of-Thought.  _arXiv preprint arXiv:2306.07932_.
  * Chen et al. (2015) Chen, J.; Cohn, A. G.; Liu, D.; Wang, S.; Ouyang, J.; and Yu, Q. 2015.  A survey of qualitative spatial representations.  _The Knowledge Engineering Review_ , 30(1): 106–136.
  * Cohn et al. (1997) Cohn, A. G.; Bennett, B.; Gooday, J.; and Gotts, N. M. 1997.  Qualitative Spatial Representation and Reasoning with the Region Connection Calculus.  _Geoinformatica_ , 1(3): 275–316.
  * Cohn and Hazarika (2001) Cohn, A. G.; and Hazarika, S. M. 2001.  Qualitative spatial representation and reasoning: An overview.  _Fundamenta Informaticae_ , 46(1-2): 1–29.
  * Cohn and Hernandez-Orallo (2023) Cohn, A. G.; and Hernandez-Orallo, J. 2023.  Dialectical language model evaluation: An initial appraisal of the commonsense spatial reasoning abilities of LLMs.  _arXiv preprint arXiv:2304.11164_.
  * Cohn and Renz (2008) Cohn, A. G.; and Renz, J. 2008.  _Handbook of Knowledge Representation_ , chapter Qualitative Spatial Representation and Reasoning.  Elsevier.
  * Dang et al. (2022) Dang, H.; Mecke, L.; Lehmann, F.; Goller, S.; and Buschek, D. 2022.  How to Prompt? Opportunities and Challenges of Zero- and Few-Shot Learning for Human-AI Interaction in Creative Applications of Generative Models.  arXiv:2209.01390.
  * Devlin et al. (2019) Devlin, J.; Chang, M.-W.; Lee, K.; and Toutanova, K. 2019.  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , 4171–4186. Minneapolis, Minnesota: Association for Computational Linguistics.
  * Freksa (1992) Freksa, C. 1992.  Temporal reasoning based on semi-intervals.  _Artificial intelligence_ , 54(1-2): 199–227.
  * Gerevini and Renz (2002) Gerevini, A.; and Renz, J. 2002.  Combining topological and size information for spatial reasoning.  _Artificial Intelligence_ , 137(1-2): 1–42.
  * Guesgen (1989) Guesgen, H. W. 1989.  Spatial reasoning based on Allen’s temporal logic, Technical Report TR-89-049.  Technical report, International Computer Science Institute Berkeley.
  * Jonsson and Drakengren (1997) Jonsson, P.; and Drakengren, T. 1997.  A Complete Classification of Tractability in RCC-5.  _J. Artif. Intell. Res._ , 6: 211–221.
  * Knauff, Rauh, and Schlieder (1995) Knauff, M.; Rauh, R.; and Schlieder, C. 1995.  Preferred mental models in qualitative spatial reasoning: A cognitive assessment of Allen’s calculus.  In _Proceedings of the seventeenth annual conference of the Cognitive Science Society_ , 200–205.
  * Moratz (2006) Moratz, R. 2006.  Representing relative direction as a binary relation of oriented points.  In _ECAI_ , volume 6, 407–411.
  * Mukerjee and Joe (1990) Mukerjee, A.; and Joe, G. 1990.  A Qualitative Model for Space.  In Shrobe, H. E.; Dietterich, T. G.; and Swartout, W. R., eds., _Proceedings of the 8th National Conference on Artificial Intelligence. Boston, Massachusetts, USA, July 29 - August 3, 1990, 2 Volumes_ , 721–727. AAAI Press / The MIT Press.
  * Ragni, Tseden, and Knauff (2007) Ragni, M.; Tseden, B.; and Knauff, M. 2007.  Cross-cultural similarities in topological reasoning.  In _Spatial Information Theory: 8th International Conference, COSIT 2007_ , 32–46. Springer.
  * Randell and Cohn (1989) Randell, D.; and Cohn, A. 1989.  Modelling topological and metrical properties in physical processes.  In _Proceedings of the First International Conference on Principles of Knowledge Representation and Reasoning_ , 357–368.
  * Randell, Cui, and Cohn (1992) Randell, D.; Cui, Z.; and Cohn, A. G. 1992.  A Spatial Logic Based on Regions and Connection.  In _3rd International Conference on Knowledge Representation and Reasoning, 1992_ , 165–176.
  * Roumeliotis and Tselikas (2023) Roumeliotis, K. I.; and Tselikas, N. D. 2023.  ChatGPT and Open-AI Models: A Preliminary Review.  _Future Internet_ , 15(6): 192.
  * Wei et al. (2022) Wei, J.; Wang, X.; Schuurmans, D.; Bosma, M.; Chi, E.; Le, Q.; and Zhou, D. 2022\.  Chain of thought prompting elicits reasoning in large language models.  _arXiv preprint arXiv:2201.11903_.
  * Yao et al. (2023) Yao, S.; Yu, D.; Zhao, J.; Shafran, I.; Griffiths, T. L.; Cao, Y.; and Narasimhan, K. 2023.  Tree of thoughts: Deliberate problem solving with large language models.  _arXiv preprint arXiv:2305.10601_.


# Real-Time Escape Route Generation in Low Visibility Environments using
Reinforcement Learning

Author Names Omitted for Anonymous Review. Paper-ID [add your ID here]  Hari
Srikanth  Mira Loma High School
Sacramento, California 95761
Telephone: (916) 693-0389
Email<EMAIL_ADDRESS>
###### Abstract

Structure fires are responsible for the majority of fire-related deaths
nationwide. In order to assist with the rapid evacuation of trapped people,
this paper proposes the use of a system that determines optimal search paths
for firefighters and exit paths for civilians in real time based on
environmental measurements. Through the use of a LiDAR mapping system
evaluated and verified by a trust range derived from sonar and smoke
concentration data, a proposed solution to low visibility mapping is tested.
These independent point clouds are then used to create distinct maps, which
are merged through the use of a RANSAC based alignment methodology and
simplified into a visibility graph. Temperature and humidity data are then
used to label each node with a danger score, creating an environment tensor.
After demonstrating how a Linear Function Approximation based Natural Policy
Gradient RL methodology outperforms more complex competitors with respect to
robustness and speed, this paper outlines two systems (savior and refugee)
that process the environment tensor to create safe rescue and escape routes,
respectively.

##  I Motivation

Structure fires pose a serious danger to many communities across the world,
and this problem has only grown over the past decade, with over a 33% increase
in deaths from fire since 2012 (Hall and Evarts [4]). Of these, structure
fires are most prominent, constituting 79% of civilian deaths and 86% of
civilian injuries. These casualties stem from people being trapped within the
building for an extended period of time, resulting in exposure to collapsing
structures, high heat/flames, and toxic smoke. In addition, the unpredictable
nature of fires reduces the efficacy of predetermined escape strategies. In
order to determine safe rescue and escape routes, a live mapping of the
environment with the measurement of aspects of interest is necessary. However,
current methodologies are insufficient to handle this task, due to an
inability to map rapidly in densely obstructed environments, as well as slow
onboard machine learning.

##  II Background

The first task for the system is to generate a map of the environment in real
time. A well researched method for this is Simultaneous Localization and
Mapping (SLAM). SLAM is a process where the agent pose is determined within an
environment which is (simultaneously) being mapped. While SLAM methodologies
have been well tested in regular environments, low visibility environments
clouded with smoke or dust pose a problem for most systems, as they utilize
optical rangefinding systems. In order to develop a resilient mapping system,
various solutions have been proposed:

  * •

SLAM with Visual and Thermal Imaging Cameras (Brunner et al. [1]): Utilizes
thermal imaging in order to counterbalance visual camera obfuscation. However,
increased reliance on thermal data resulted in a decrease in localization
accuracy.

  * •

SLAM with laser range finders and 94 GHz Frequency Modulated Continuous Wave
Radar (Gerardo-Castro and Peynot [3]): Proposed a system with mm-wave radar to
penetrate dense smoke that laser scanners cannot. Very precise, but some false
negatives remained. In addition, the custom antenna system is very costly and
locks mobility.

  * •

SLAM with laser range finders and sonar sensor array (Machado Santos et al.
[6]): Uses a fuzzy logic system Couceiro et al. [2] to decide between sensor
inputs, reducing the risk of false negatives.

Regardless of the specific mapping methodology, the map resolution and
effective range drastically decrease as the adversarial noise parameter
increases. Given these limitations, the use of Multi-Robot Systems (MRS) is
necessary to optimize mapping speed and robustness, with a Random Sample
Consensus (RANSAC) based map unification to ensure live processing does not
become too computationally expensive (Lázaro et al. [5]). The second aspect of
the problem is the determination of optimal rescue or escape paths. In order
to do this, the situation can be modeled through a Markov Decision Process and
solved using Reinforcement Learning. Reinforcement Learning (RL) is a paradigm
where an agent seeks to maximize the reward it gains through refining its
policy. At each timestep t, the agent observes the environmental state and
according to some policy $\pi$ it takes some action. This action changes the
environmental state and returns some reward, and this is used to retrain the
policy. A well researched RL method is Actor Critic, which utilizes an actor
that updates the policy and a critic which evaluates said policy for fast
optimization. While industry standard algorithms Trust Region Policy
Optimization (Schulman et al. [8]) and Proximal Policy Optimization Schulman
et al. [9] utilize complex neural networks to estimate the value function, my
prior research on reinforcement learning determined that a Linear Function
Approximation based Natural Policy Gradient algorithm (LFA-NPG) would
determine the optimal policy with much less iterations than either TRPO and
PPO, in low dimensional standard and sparse reward RL benchmarks. In addition,
my robustness analysis determined that LFA-NPG was significantly more noise
resistant to adversarial noise than TRPO and PPO, maintaining identical
performance across data sampled within 20% of the true value. Although the
determination of navigation paths appears to be optimized for operations
research methodology, the unique nature of each fire does not lend itself to
the determination of the heuristics often employed in these algorithms.

##  III Design

The main necessities of a real world robotics system are as follows:

  * •

Speed: Data must be sampled at a high rate, in order for the system to react
appropriately

  * •

Robustness: Algorithms must be resistant to adversarial noise, an inherent
factor when live systems are considered

In addition to these two criteria, the search and rescue application
motivating this paper requires the system to have some measure of its own
accuracy. Fig. 1 details the proposed system, taking these criteria into
consideration.

Figure 1: System Organization

###  III-A Data Acquisition & Processing

Given that each second the system takes to find an exit path is a second of
exposure to hazardous conditions, the key priority of my solution is speed.
Fast navigation, fast mapping, fast processing, and fast path determination
take precedence over exact optimization or high resolution. In order to map at
this rate, the system will utilize a fleet of autonomous drones. Each drone
will be outfitted with a mapping module (consisting of a LiDAR rangefinder and
four sonar scanners) and a sensor module (to collect temperature and humidity
data, as well as determine a quantitative adversarial noise parameter). Using
the temperature and humidity data, the map will then be populated with scores
quantifying the danger of each area of the building. From here, the LFA-NPG
agent will determine the optimal paths, both for rescue (for the firefighters)
and for escape (for the inhabitants). The base hardware for each device is a
light quadrotor with ROS on a RPi running ubuntu. Using 3D printed parts to
keep the weight as low as possible, the drone will have a target thrust/weight
of 2:1. Each RPi interfaces with a Navio 2 flight control unit running an
ardupilot APM autopilot. ROS features a variety of nodes to run SLAM
calculations, such as HectorSLAM, GMapping, KartoSLAM, CoreSLAM, and LagoSLAM.
Of these, GMapping is the most optimal for mapping in a smoky environment, as
it remains robust and has a large amount of support.

The raw data that the mapping module collects from each drone is laser scan
information from the LiDAR, four sonar distance measurements, and a measure of
adversarial noise. The core of the mapping system is the LiDAR data: even
though its accuracy is compromised in smoky environments, it is the most high-
resolution sensor available within the price range. The sonar data comes from
four sensors at 90 degrees to one another. The sonar data is then transformed
to a laser scan frame to directly evaluate it in comparison to the LiDAR data.
In order to determine which sensor data should be used (LiDAR for high
resolution vs Sonar for higher accuracy) when constructing the map, the device
incorporates a fuzzy logic system, through a measurement of confidence called
the Trust Range. The trust range is determined by the difference between laser
and sonar data and the adversarial noise parameter. In a design intended for
real world use, the adversarial noise parameter in the use case of a fire
would be determined by a particulate matter sensor. However, the logistics of
producing real smoke multiple times for testing is a safety concern, and as
such the effect of smoke is simulated through the use of a smoke machine.
Because smoke machines do not produce particulate smoke, but rather vaporize
alcohol-based fog solution, the noise parameter is instead determined through
the use of an alcohol sensor (no difference in algorithm, just uses a
different function and hardware when calculating). After the trust range is
determined, the points outside the trust range are eliminated from the data
and points within are added to an aggregate point cloud. This aggregate point
cloud includes data from each drone exploring the system. The data from each
point cloud is then passed into GMapping, creating a map. The set of maps from
each drone is then fed into the processing module.

The processing module begins by merging each of the individual maps together,
using the algorithm outlined below:

  1. 1.

Use a correlative scan matcher to identify matching edges of each scan

  2. 2.

Add each edge solution to a pool of candidates

  3. 3.

For any two maps, determine a translation between them such that one candidate
edge is identical (zero error).

  4. 4.

Determine the errors of all other candidate edges. Depending on if the error
of the other edges is high or low, it is possible to determine inliers and
outliers

  5. 5.

Once the correct translation has been determined, merge the two maps into a
global map. This global map can then be merged with the next map, which can
repeat until all maps have been integrated into the global map.

After the global map has been created, the pose of each drone is passed to the
flight controller, completing the localization aspect of SLAM. While this map
could just be exported as is, more processing and simplification is necessary
to determine exit paths. First, the merged map is smoothed into a visibility
graph, which maps complex environmental features to a significantly more
condensed set of traversable nodes. Each node is then assigned a danger score
based on temperature and humidity data. For nodes where no temperature data
has been collected, data from nearby nodes is weighted and used to estimate
the conditions at the location. The complete environment tensor with
traversable nodes and dangers is then ready to be fed into the RL training
models.

###  III-B RL Data Analysis

LFA-NPG is an policy-space RL Paradigm that utilizes the Natural Actor Critic
(Peters and Schaal [7]) Architecture, optimizing the policy through natural
gradient descent while estimating the value function (how well a given policy
will produce reward) through Linear Function Approximation. LFA-NPG was
evaluated on standard RL benchmarks, Cartpole, and Acrobot (a sparse-reward
environment, which is more common in real world robotics systems). Its
performance is as shown in Fig. 2.

Figure 2: LFA-NPG Robustness Analysis

As shown, LFA-NPG converges to the optimal policy even as the noise parameter
(represented by the true state sampling error passed through to the model)
increases up to 20%. This is due to its simple value estimation method: in
contrast to the complex neural networks present in other forms of RL, LFA-
NPG’s simpler value estimation methodology enables it to be more resistant to
fluctuations in input. In addition to robustness, LFA-NPG also is much faster
than competing algorithms in low dimensionality applications (such as CartPole
and Acrobot). It reaches a similar level of reward compared to PPO and TRPO,
with a much faster and more logarithmic convergence, ideal for high speed use
cases. This is evidenced in Fig. 3.

Figure 3: LFA-NPG Model Convergence w.r.t Iterations & Time

Considering that the data that can be collected during a live system activity
is limited in scope, it will tend to be sufficiently low dimensionality enough
for LFA-NPG to function at an optimal level. Moreover, RL methodology has been
shown to be competitive with leading operations research algorithms currently
used for traversal problems, and since the structure of an MDP is more
cohesive with raw sensor inputs than the heuristic information necessary for
OR algorithms, it is uniquely positioned as the optimal real time path planner
for real world robotics.

In the presented use case, the environment tensor can be used to inform both
the optimal strategy for a firefighter and the optimal strategy for a trapped
civilian. As such, two RL models will train simultaneously using the same base
state space, with different reward functions. The savior system determines the
strategy for the firefighter: Given the locations of possible entryways into
the building, the system then determines the optimal paths to navigate to any
point within the building, in the form of a tree. The refugee system is used
to determine the escape routes. The starting nodes are all lethal areas
(determined by danger score) and have a reward function dependent on duration
and level of exposure to hazardous conditions. These two objectives are each
used to formulate distinct MDPs, which LFA-NPG based solvers optimize through
the use of the actor critic architecture.

In order to determine the feasibility of this technology in real world
applications, a minimum of 2 drones must be utilized (to evaluate the efficacy
of the merged map) in a simulated indoor fire. This indoor fire will be
simulated through the use of a smoke machine, which the glycol sensor uses to
determine the danger parameter. As it is difficult to simulate the high
temperatures of a fire, the environment tensor will instead be fed a possible
matrix of danger scores within the building.

##  IV Conclusion

This work considers the relevant Search and Rescue use case for autonomous
robot systems for the traversal of a burning building. In addition to
discussing the efficacy of current low visibility perception systems, this
work proposes a novel perception system that makes use of a multi agent
mapping array. Each individual agent prioritizes certain flows of data through
a trust range, determined by the measure of adversarial noise, which is
additionally used to evaluate the safety of each location in the map. It goes
on to discuss a possible way in which this data can be analyzed, namely for
the generation of escape/rescue routes that maximize safety. In order to
accomplish this, the paper recommends the use of a linear function
approximation based natural policy gradient reinforcement learning
methodology, demonstrating its high speed and strong resistance to adversarial
noise in sufficiently low dimensional systems.

## References

  * Brunner et al. [2013] Christopher Brunner, Thierry Peynot, Teresa Vidal-Calleja, and James Underwood.  Selective combination of visual and thermal imaging for resilient localization in adverse conditions: Day and night, smoke and fire.  _Journal of Field Robotics_ , 30(4):641–666, 2013\.  doi: https://doi.org/10.1002/rob.21464.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/rob.21464.
  * Couceiro et al. [2012] Micael S. Couceiro, J.A. Tenreiro Machado, Rui P. Rocha, and Nuno M.F. Ferreira.  A fuzzified systematic adjustment of the robotic darwinian pso.  _Robotics and Autonomous Systems_ , 60(12):1625–1639, 2012.  ISSN 0921-8890.  doi: https://doi.org/10.1016/j.robot.2012.09.021.  URL https://www.sciencedirect.com/science/article/pii/S0921889012001753.
  * Gerardo-Castro and Peynot [2012] Marcos Paul Gerardo-Castro and Thierry Peynot.  Laser-to-radar sensing redundancy for resilient perception in adverse environmental conditions.  01 2012.  URL https://www.researchgate.net/publication/264545306_Laser-to-radar_sensing_redundancy_for_resilient_perception_in_adverse_environmental_conditions.
  * Hall and Evarts [2021] Shelby Hall and Ben Evarts.  Fire loss in the united states during 2021, 2021.  URL https://www.nfpa.org/~/media/fd0144a044c84fc5baf90c05c04890b7.ashx#:~:text=In%202021%2C%20local%20fire%20departments,14%2C700%20reported%20civilian%20fire%20injuries.
  * Lázaro et al. [2013] M. T. Lázaro, L. M. Paz, P. Piniés, J. A. Castellanos, and G. Grisetti.  Multi-robot slam using condensed measurements.  In _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , pages 1069–1076, Nov 2013.  doi: 10.1109/IROS.2013.6696483.
  * Machado Santos et al. [2014] João Machado Santos, Micael S. Couceiro, David Portugal, and Rui P. Rocha.  Fusing sonars and lrf data to perform slam in reduced visibility scenarios.  In _2014 IEEE International Conference on Autonomous Robot Systems and Competitions (ICARSC)_ , pages 116–121, 2014.  doi: 10.1109/ICARSC.2014.6849772.
  * Peters and Schaal [2008] Jan Peters and Stefan Schaal.  Natural actor-critic.  _Neurocomputing_ , 71(7):1180–1190, 2008.  ISSN 0925-2312.  doi: https://doi.org/10.1016/j.neucom.2007.11.026.  URL https://www.sciencedirect.com/science/article/pii/S0925231208000532.  Progress in Modeling, Theory, and Application of Computational Intelligenc.
  * Schulman et al. [2017a] John Schulman, Sergey Levine, Philipp Moritz, Michael I. Jordan, and Pieter Abbeel.  Trust region policy optimization, 2017a.  URL https://arxiv.org/abs/1502.05477.
  * Schulman et al. [2017b] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov.  Proximal policy optimization algorithms, 2017b.