††thanks: These authors contributed equally to this work
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>These authors contributed equally to this work
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
# Multi-time quantum process tomography of a superconducting qubit

Christina Giarmatzi  School of Computer Science, University of Technology
Sydney, Ultimo, Sydney, New South Wales 2007, Australia  ARC Centre of
Excellence for Engineered Quantum Systems, St. Lucia, Brisbane, Queensland
4072, Australia  Tyler Jones  ARC Centre of Excellence for Engineered Quantum
Systems, St. Lucia, Brisbane, Queensland 4072, Australia  School of Maths and
Physics, University of Queensland, St. Lucia, Brisbane, Queensland 4072,
Australia  Fiasqo, Brisbane, Queensland 4072, Australia  Alexei Gilchrist  ARC
Centre of Excellence for Engineered Quantum Systems, St. Lucia, Brisbane,
Queensland 4072, Australia  Department of Physics and Astronomy, Macquarie
University, Sydney, New South Wales 2122, Australia  Prasanna Pakkiam  ARC
Centre of Excellence for Engineered Quantum Systems, St. Lucia, Brisbane,
Queensland 4072, Australia  School of Maths and Physics, University of
Queensland, St. Lucia, Brisbane, Queensland 4072, Australia  Arkady Fedorov
ARC Centre of Excellence for Engineered Quantum Systems, St. Lucia, Brisbane,
Queensland 4072, Australia  School of Maths and Physics, University of
Queensland, St. Lucia, Brisbane, Queensland 4072, Australia  Fabio Costa
School of Maths and Physics, University of Queensland, St. Lucia, Brisbane,
Queensland 4072, Australia  Nordita, Stockholm University and KTH Royal
Institute of Technology, Stockholm, 106 91, Sweden

###### Abstract

Non-Markovian noise poses a formidable challenge to the scalability of quantum
devices, being both ubiquitous in current quantum hardware and notoriously
difficult to characterise. This challenge arises from the need for a full
reconstruction of a multi-time process, a task that has proven elusive in
previous efforts. In this work, we achieve the milestone of complete
tomographic characterisation of a multi-time quantum process on a
superconducting qubit by employing sequential measure-and-prepare operations
with an experimentally motivated post-processing technique, utilising both in-
house and cloud-based superconducting quantum processors. Employing the
process matrix formalism, we reveal intricate landscapes of non-Markovian
noise and provide evidence that components of the noise originate from quantum
sources. Our findings and techniques have significant implications for
advancing error-mitigation strategies and enhancing the scalability of quantum
devices.

††preprint: APS/123-QED

##  I Introduction

The characterisation of noise is a critical requirement for the advancement of
quantum technologies. Noise is present in all current quantum devices as they
interact with their environment; for many devices, this includes non-Markovian
noise due to temporal correlations across a multi-time process. While most
current techniques assume Markovian noise, this assumption fails in realistic
quantum devices and non-Markovian noise has been reported in state-of-the-art
quantum computing devices (such as IBM and Google) [1, 2, 3]. Past general
approaches to capture non-Markovian noise only provide access to two times of
a quantum process through dynamical maps [4, 5, 6, 7, 8, 9, 10]. Beyond this,
a number of techniques designed to reduce non-Markovian noise have been
developed [11, 12, 13, 14, 15], and have been shown to increase computational
accuracy [11], but all rely on ad-hoc methods, and a rigorous approach with
potential for scalability, is still lacking.

Recently, a formalism was developed that constructs a _process matrix_ [16,
17], which encodes all multi-time correlations in a quantum process. In this
formalism, non-Markovian noise is rigorously captured and can be further
investigated in terms of its nature (quantum or classical) and amount [18,
19], as well as memory length and location in the device [20] with a prospect
of scalability [21] and extending current characterisation methods to non-
Markovian regimes [22]. The approach has been used in distinct experimental
scenarios, successfully detecting non-Markovian noise [1, 23, 20, 24, 25, 26,
27]. However, full reconstruction of a multi-time process matrix requires
sequential informationally complete measure-and-prepare operations at every
time step. Previous tomography approaches did not have access to such
operations, and thus could only construct a “restricted” process matrix, out
of a limited set of experimental data [27, 28, 29, 30]. Process matrices also
have applications in quantum foundations, as a way to explore indefinite
causal order [31, 32, 33, 34, 35]; in this context, tomography of a process
matrix without causal order has been recently achieved in a quantum-optics
setup [36].

Here, we present the first full process tomography of a multi-time quantum
process implemented on a superconducting qubit and obtain measures of non-
Markovianity. We use devices from the labs of the University of Queensland and
the cloud capabilities of IBM Quantum [37], with varying parameters,
evaluating a total of ten three-time quantum processes. We measure general and
quantum non-Markovian noise in all cases and build a theoretical model to
compare our findings. This offers a robust method for characterising non-
Markovian noise, with potential for scalability, significantly advancing both
our theoretical understanding and the development of noise mitigation
techniques.

##  II Multi-time quantum process tomography

Consistent with our experimental setup, we consider a scenario where a quantum
system undergoes a sequence of three operations at specified times, see Fig 1.
The purpose of multi-time process tomography is to extract, from a finite set
of measurements, sufficient information that allows one to predict the joint
probability for any other sequence of measurements [38, 39, 30]. In the
presence of general system-environment correlations, knowledge of the density
matrix of the system at each time does not capture all multi-time
correlations. The mathematical object that encodes all needed information is
the process matrix—also known as quantum strategy [40], comb [38], or process
tensor [41]—which is a convenient representation of a quantum stochastic
process [42].

Figure 1: A three-time quantum process involving state preparation and two
gates acting on both systems.

Let us label $A,B,C$ the times at which the operations take place. The
operations are represented by general completely positive (CP), trace-non-
increasing maps from an input to an output Hilbert space so that, for example
at time $A$, the joint input-output space is
$\mathcal{H}^{A}=\mathcal{H}^{A_{I}}\otimes\mathcal{H}^{A_{O}}$. The map can
be expressed in Choi form [43, 44] (see Methods) as an operator $M^{A}_{a\mid
x}\in\mathcal{L}(\mathcal{H}^{A})$, $M^{A}_{a\mid x}\geq 0$, where $a$ denotes
the measurement outcome, $x$ denotes the setting (i.e., the choice of
operation), and $\mathcal{L}(\mathcal{H})$ the space of linear operators over
$\mathcal{H}$ [45]. For each setting choice $x$, a complete set of outcomes
defines an instrument [46], expressed as a set $\\{M^{A}_{a\mid x}\\}_{a}$
such that $\sum_{a}\operatorname{Tr}_{A_{O}}M^{A}_{a\mid
x}=\mathds{1}^{A_{I}}$. The process matrix is an operator on the joint space
of all inputs and outputs,
$W^{ABC}\in\mathcal{L}(\mathcal{H}^{A}\otimes\mathcal{H}^{B}\otimes\mathcal{H}^{C})$,
$W^{ABC}\geq 0$, which allows one to calculate arbitrary joint probabilities
for a sequence of outcomes $a,b,c$, given settings $x,y,z$, through the
generalised Born rule [47]

$p(a,b,c\mid x,y,z)=\operatorname{Tr}[W^{ABC}(M_{a\mid x}^{A}\otimes M_{b\mid
y}^{B}\otimes M_{c\mid z}^{C})].$ (1)

Formally, we can regard $W^{ABC}$ as a “density operator over time”, where
each time instant corresponds to a pair of subsystems, an input and an output
space. As output spaces encode the causal influence from an operation to
future ones, the final output space can always be discarded (hence, the final
time $C$ has only an input space). Furthermore in our scenario the first
operation is a strong state preparation, which decorrelates the system from
the environment allowing us to discard the first input space $A_{I}$. Thus for
our experiment the process matrix is $W^{ABC}\equiv W^{A_{O}B_{I}B_{O}C_{I}}$
corresponding to a four-qubit state.

By selecting informationally complete instruments at each time, and exploiting
the process-state isomorphism, one can reconstruct the process matrix using
similar techniques as for ordinary state tomography, namely one can invert Eq.
(1) and write $W$ as a function of the probabilities. A convenient choice for
the instruments is projective measurements $\\{\Pi_{a\mid x}\\}_{a}$
immediately followed by the preparation of a state $\rho_{x}$, which
correspond to operations of the form $M_{a\mid x}^{A}=\Pi_{a\mid
x}^{A_{I}}\otimes\rho_{x}^{A^{T}_{O}}$, where $A^{T}_{O}$ denotes
transposition on system $A_{O}$ and reflects our choice of convention,
consistent with Eq. (1). This choice makes the protocol equivalent to
tomography of a two-qubit channel that maps systems $A_{O}B_{O}$ to
$B_{I}C_{I}$, with $W^{ABC}$ as its Choi matrix. As each measurement projects
the system onto the measured state, at face value re-preparing an independent
state requires a conditional qubit rotation. The fast feed-forward needed for
this procedure was the main obstacle preventing full tomography in previous
work [27, 28, 29, 30]. As detailed in the next section, we overcome this
requirement through an appropriate post-processing scheme.

##  III Multi-time tomography protocol of a superconducting qubit

The experiment was conducted separately on two sets of superconducting
hardware; an in-house processor at the University of Queensland, which will be
referred to as uq, and a cloud processor from IBM Quantum with the moniker
ibm_perth.

The experiment conducted on the uq processor used two capacitively-coupled
flux-tunable transmons on a multi-qubit chip, mounted in a dilution
refrigerator with a base temperature of around 15 mK. Each qubit is coupled to
a dedicated resonator for readout. Flux lines are used in conjunction with a
magnetic coil fixed to the chip to control transition frequencies, providing a
degree of freedom by which qubit-qubit interaction strengths can be
controlled.

Figure 2: A visualisation of the pulse protocol used to generate tomography
data. The solid line on which pulses are applied represents the system qubit.
The dotted line represents the environment, with the two yellow boxes
indicating the system-environment interaction that multi-time process
tomography aims to characterise.

The experiments on the ibm_perth processor also utilised a two-transmon subset
of a larger multi-qubit chip, but the system is comprised of fixed-frequency
transmons. This removes flux control and thereby fixes the effective qubit
interaction strength. In both experiments, the system of interest is a single
qubit, while the other qubits are part of the environment (which also
comprises any other degrees of freedom interacting with the system, as well as
classical sources of noise).

To implement the three-time process tomography, we set the Pauli operators as
our observables $\mathcal{M}=\\{X,Y,Z\\}$ and their eigenstates as the input
states $\mathcal{B}=\\{X_{+},X_{-},Y_{+},Y_{-},Z_{+},Z_{-}\\}$. These form a
tomographically (over) complete basis. We also define a set of operations
within a basis gate set
$\mathcal{G}=\\{I,R_{X}(\pi/2),R_{X}(-\pi/2),R_{Y}(\pi/2),R_{Y}(-\pi/2),R_{X}(\pi)\\}$.
These operations are, by necessity, sufficient to prepare and measure all
combinations of basis states and observables. Multi-time process tomography
consists of temporally stacking two traditional process tomography protocols,
forming a (prepare $\rightarrow$ process $\rightarrow$ measure $\rightarrow$
prepare $\rightarrow$ process $\rightarrow$ measure) protocol, where
preparations are on spaces $A_{O},B_{O}$, while measurements are on
$B_{I},C_{I}$. In these experiments, no operation is performed between each
preparation and the subsequent measurement, so the idealised process is
composed of two identity channels from $A_{O}$ to $B_{I}$ and from $B_{O}$ to
$C_{I}$. Tomography aims to characterise departures from this idealisation,
representing noise in the process due to system-environment interaction. The
preparation states are drawn from $\mathcal{B}$ and measurement observables
from $\mathcal{M}$, resulting in
$\lvert\mathcal{B}\rvert^{2}\lvert\mathcal{M}\rvert^{2}=324$ unique protocols.
The general form of this protocol is illustrated in Fig. 2, where we perform
it on the system qubit (solid bottom line) and represent the environment with
the top dotted line. Together with Fig. 1, we can see the multi-time quantum
process and the times where we make the operations to gather the data.

To begin, all qubits were initialised in the ground state through a standard
long wait-time routine. The system qubit is then prepared into some basis
state $b_{0}\in\mathcal{B}$ using the appropriate rotation
$p_{0}\in\mathcal{G}$. This occurs at $A$. The system and environment are then
left to freely evolve for time $t_{1}$.

At $B$ we apply a projective measurement of the observable
$m_{0}\in\mathcal{M}$ on the system qubit, consisting of some rotation gate
$r_{0}\in\mathcal{G}$ to facilitate measurement of the correct axis, followed
by a readout signal. The reflected readout signal is then classified as 0 or 1
and the result is stored in memory. After the mid-circuit measurement, another
rotation $p_{1}\in\mathcal{G}$ is applied to prepare a new qubit state.
Conditional on the eigenstate collapse in the mid-circuit measurement, this
operation prepares the system qubit into one of the basis states
$b_{1},\tilde{b}_{1}\in\mathcal{B}$. For example, if a protocol measures
$m_{0}=X$ and applies rotation $p_{1}=I$, the prepared state is either
$b_{1}=X_{+}$ or $\tilde{b}_{1}=X_{-}$. The actual prepared state is then
identified by checking the result of the mid-circuit measurement. In this way,
we are able to obtain a complete set of preparations over all experimental
runs without the need for each prepared state to be chosen deterministically.
This simple but critical post-processing technique is designed to avoid the
need for dynamic, fast feed-forward control, which would otherwise be required
for conditional operations. Following these operations, the system and
environment are left to freely evolve for time $t_{2}$.

Finally, at $C$ we measure our second observable $m_{1}\in\mathcal{M}$. In
all, the protocol is repeated 8192 times for each of 324 prepare-measure-
prepare-measure permutations, resulting in approximately $5.3\times 10^{6}$
individual shots per process characterisation experiment (one on uq and nine
at ibm_perth). This data is shaped into the number of realised counts we
obtain for each of the 324 permutations. An example of these permutations is
$(X_{+},X_{-},Y_{+},Z_{-})$ denoting that at $A$ we prepare the state $X_{+}$,
at $B$ we projectively measure $X_{-}$ and then re-prepare $Y_{+}$, and at $C$
we projectively measure $Z_{-}$.

After obtaining the data, we are ready to calculate the experimental process
matrix. A process matrix can written in the Pauli basis,
$W=\sum_{i,j,k,l}\alpha_{i,j,k,l}\sigma_{i}\otimes\sigma_{j}\otimes\sigma_{k}\otimes\sigma_{l}$,
where (i,j,k,l) = [0,3] denotes the Pauli basis
${\mathds{1},\sigma_{x},\sigma_{y},\sigma_{z}}$. From Eq. 1, we see that if we
apply Pauli operations at $A,B,C$, their expectation values will be the
coefficients
$\langle\sigma_{i},\sigma_{j},\sigma_{k},\sigma_{l}\rangle=\alpha_{i,j,k,l}$.
Hence, we can reconstruct the experimental process matrix, by calculating the
expectation values of the sigma terms from the experimental data

$W_{exp}=\frac{1}{16}\sum_{i,j,k,l}\langle\sigma_{i},\sigma_{j},\sigma_{k},\sigma_{l}\rangle\sigma_{i}\otimes\sigma_{j}\otimes\sigma_{k}\otimes\sigma_{l}.$
(2)

$W_{exp}$ is the process matrix we construct from the experimental data, but
it typically does not represent a physical one (it might be non-positive, or
fail the process matrix constraints). To obtain a physical process,
$W_{phys}$, we use SemiDefinite Programming (SDP), where we find a positive
matrix that minimises its distance (we used Frobenius distance) with the
$W_{exp}$. The SDP is the following:

variable $W_{phys}$ Hermitian Semidefinite minimize $\displaystyle\ \
\texttt{Frobenius Norm}(W_{exp}-W_{phys})$ (3) subject to $\displaystyle\ \
L_{V}(W_{phys})==W_{phys}$ $\displaystyle\ \ \text{Tr}(W_{phys})==4,$

where $L_{V}$ is an operator that projects a matrix to the space of valid
process matrices [32] and the normalisation constraint helps the SDP to bound
the solution. Using the $W_{phys}$ we report in the Results section our
findings on tomography and non-Markovianity.

Our tomography protocol differs from related works in that we implement full
multi-time process tomography by reconstructing the full process matrix.
Concretely, Xiang _et al._ implement projective mid-circuit measurements but
limit their analysis to an under-complete set of POVMs to construct a
‘restricted’ process matrix [28]. Similarly, White _et al._ restrict their
tomography to only sequences of unitary operations and final measurement, with
reference to the difficulty of mid-circuit measurement in superconducting
processors [30]. Related work by the same authors extends their analysis to
non-unital maps by introducing an ancillary qubit and entangling operations,
obtaining a separate class of restricted information known as classical
shadows [29].

##  IV Results: Tomography and non-Markovian noise

Figure 3 provides insight into the qualitative differences between two
experimental (physical) process matrices, $W_{phys}$, using data from the in-
house uq processor and from the ibm_perth processor. The process matrix from
ibm_perth is visually almost indistinguishable from a noiseless process with
identity channels between the times $A$, $B$ and $C$, while the uq process
maintains the qualitative structure of the identity channels with evident
noise. This is a reflection of the operational fidelity of the processors, in
conjunction with the significantly higher interaction strength between the
system and the rest of the qubits on the uq processor (which was controlled by
setting the appropriate qubit-qubit detuning).

Overall, we took tomographic data for one process on the uq processor and nine
on ibm_perth with varying $t_{1}$ and $t_{2}$. Variations in free evolution
times were explored after simulating each system with a simple one-qubit
environment model and finding that both general and quantum non-Markovianity
depend highly on these times, with a particular structure (see Methods). No
such structure was observed in the experiments we conducted. This indicates
that even though in both platforms the system qubit couples directly only to
one other qubit, such interaction does not account for all environmental
effects.

Figure 3: A visualisation of the absolute valued entries of the full multi-
time process matrix, constructed with tomographic data from (a) the uq
processor and (b) the ibm_perth processor.

Given the physical process matrices for each of the ten processes
characterised (broadly denoted as $W_{phys}$), we aim to detect and quantify
non-Markovianity. To this end, we calculate a measure of distance of each
$W_{phys}$ with its closest Markovian process matrix $W_{M}$. We construct
$W_{M}$ by discarding temporal correlations, combining the two channels [30,
48]

$\displaystyle W_{M}=T_{1}^{A_{O}B_{I}}\otimes T_{2}^{B_{O}C_{I}},$ (4)

where $T_{1}^{A_{O}B_{I}}=Tr_{B_{O}C_{I}}W_{phys}$ and
$T_{2}^{B_{O}C_{I}}=Tr_{A_{O}B_{I}}W_{phys}$. As a measure of distance, we
consider the Jensen-Shannon divergence (JSD) for two matrices $W_{1}$ and
$W_{2}$:

$JSD\left(W_{1},W_{2}\right)=H\left(\frac{W_{1}+W_{2}}{2}\right)-\frac{1}{2}\left(H\left(W_{1}\right)+H\left(W_{2}\right)\right),$
(5)

where $H(W)=-\operatorname{Tr}{W\log{W}}$ is the von Neumann entropy. The
square root of the JSD operates as a distance measure between states or
channels, and has utility in particular for quantifying non-Markovian dynamics
[49]. The results for $\sqrt{JSD}$ across the ten processes are reported in
Table 1. These indicates that non-Markovian effects were significant in both
platforms for all choices of parameters, but were more pernicious on the uq
processor. This is consistent with the appreciably stronger idle qubit-qubit
interaction on this chip, as seen in Table 2.

| Evolution times | General non-Mark. noise | Quantum non-Mark. noise
---|---|---|---
| ($t_{1},t_{2}$) ns |  ${\sqrt{JSD}}$ ${\times 10^{-1}}$ |  ${\mathcal{N}}$ ${\times 10^{-2}}$
uq | (97, 97) | $5.04\pm 0.07$ | $8.66\pm 0.86$
ibm_perth | (21, 21) | $1.77\pm 0.17$ | $1.61\pm 0.58$
| (21, 24) | $1.64\pm 0.17$ | $1.66\pm 0.61$
| (21, 28) | $1.81\pm 0.15$ | $2.22\pm 0.69$
| (24, 21) | $1.60\pm 0.16$ | $1.94\pm 0.64$
| (24, 24) | $1.85\pm 0.18$ | $1.71\pm 0.65$
| (24, 28) | $2.06\pm 0.14$ | $2.44\pm 0.68$
| (28, 21) | $1.96\pm 0.15$ | $2.00\pm 0.71$
| (28, 24) | $1.85\pm 0.14$ | $2.04\pm 0.67$
| (28, 28) | $1.76\pm 0.17$ | $1.83\pm 0.62$
Table 1: The square root Jensen-Shannon divergence ($\sqrt{JSD}$) and
negativity ($\mathcal{N}$) of the multi-time process matrices reconstructed
for each processor. Each experiment characterises two periods of free
evolution, which are denoted by labels $t_{1}$ and $t_{2}$, shown in
nanoseconds. The ibm_perth data consisted of nine different experimental runs
with slightly varying free evolution times, for the purposes of comparison to
a theoretical model.

A natural subsequent inquiry is whether non-Markovianity originates from a
quantum source; the value of distinguishing between classical or quantum
environments lies in the manifestly different dynamics they induce [50, 51,
52, 53]. The quantum nature of non-Markovianity in the process is captured by
quantum correlations across the evolution from $A$ to $B$ and from $B$ to $C$.
This maps to entanglement across the bipartition $A_{O}B_{I}|B_{O}C_{I}$ in
the four-qubit state representation of the process,
$\rho^{A_{O}B_{I}|B_{O}C_{I}}:=\frac{W^{ABC}}{\operatorname{Tr}W^{ABC}}$ . To
measure entanglement, we take the computationally inexpensive metric of
negativity,

$\mathcal{N}(\rho)=\left\lvert\sum_{\lambda_{i}<0}\lambda_{i}\right\rvert=\sum_{i}\frac{\lvert\lambda_{i}\rvert-\lambda_{i}}{2},$
(6)

where $\rho=\rho^{(A_{O}B_{I})^{T}|B_{O}C_{I}}$, where $T$ denotes partial
transposition (it can be on either of the two subsystems, $A_{O}B_{I}$ or
$B_{O}C_{I}$), and $\lambda_{i}$ are the eigenvalues of $\rho$ [54].

Quantum multi-time correlations are observable on each processor, but they are
markedly enhanced on the uq processor. This is again consistent with the
increased interaction strength observed between the system qubit and its
primary interacting ‘memory’ qubit on the uq processor as documented in Table
2. The negativities reported over the nine ibm_perth runs all lay within each
others’ credible intervals, indicating that the quantum non-Markovianity
present in each process was effectively invariant across the ($t_{1}$,
$t_{2}$) landscape explored. In all measurements, the reported uncertainty
denotes 95% credible intervals determined by bootstrapping techniques as
described in Methods, accounting for sampling noise.

Using the natural log in Eq. 5, the $\sqrt{JSD}$ metric is bounded between 0
and $\log{2}$ [55], lending scale to the magnitude of our findings which lie
on the order of $10^{-1}$. This is a strong signal that the multi-time
processes characterised on each device are easily distinguishable from the
most likely temporally uncorrelated candidates. Negativity serves as an
entanglement monotone and is therefore bounded from below by zero (completely
separable) [54]. The monotonic nature of negativity means that, in the case of
states, it is upper bounded by the negativity of the maximally entangled
state, which for four qubits is $\frac{3}{2}$, roughly two orders of magnitude
larger than in our measurements. This state does not correspond to a valid
process matrix (because the space of states is larger than that of process
matrices) but it gives a sense of measure of our results on negativity, which
are two orders of magnitude lower than the state maximum.

In summary, we observe and quantify non-Markovian dynamics and quantum
correlations with high confidence on both uq and ibm_perth processors,
reinforcing the importance of correlated noise in the design of mitigation
techniques. Our results, interpreted alongside the bounds provided,
demonstrate the practicality of our methods for characterising this noise.

##  V Conclusion

We report the first experimental full characterisation of a multi-time quantum
process on a superconducting qubit. We implement ten three-time processes
using an in-house superconducting processor and a cloud processor from IBM
Quantum. We characterise the processes using the process matrix formalism and
provide a measure of their non-Markovian noise; both general and stemming from
quantum correlations across time. We find that non-Markovian noise is present
in all cases, with a significant proportion originating from genuine quantum
correlations. A comparison with a simple system-environment model indicates
that the influence of nearest-neighbour qubits is not sufficient to account
for all the observed quantum non-Markovianity. Our results demonstrate the
capacity of our methods to fully describe multi-time processes and confidently
measure and characterise non-Markovian noise, laying the foundation for novel
error-mitigation techniques set to impact current and future quantum devices.

###### Acknowledgements.

We thank Xin He for providing support in preparing the tomography experiments
on the uq processor, and Gerardo Paz Silva for sharing relevant literature. We
acknowledge the use of IBM Quantum services for this work. The views expressed
are those of the authors and do not reflect the official policy or position of
IBM or the IBM Quantum team. CG acknowledges funding from the Sydney Quantum
Academy Fellowship and the UTS Chancellor’s Research Fellowship. FC was
supported by the Wallenberg Initiative on Networks and Quantum Information
(WINQ). This work was supported by the Australian Research Council (ARC)
Centre of Excellence for Quantum Engineered Systems grant (CE 110001013).
Nordita is supported in part by NordForsk. Finally, we acknowledge the
traditional owners of the land on which the University of Queensland, the
University of Technology Sydney, and Macquarie University are situated, the
Turrbal and Jagera people, the Gadigal people of the Eora Nation, and
Wallumattagal clan of the Dharug Nation.

## Methods

### Experimental devices

The primary device used was a superconducting quantum processor consisting of
five flux-tunable transmons, fabricated by Quantware [56]. System
characteristics of the transmon we performed multi-time tomography on
(‘System’) and the dominant interacting transmon (‘Memory’) are provided in
Table 2. These characteristics include T1 and T2, the relaxation and dephasing
lifetimes. The ‘System’ qubit only has one nearest-neighbour coupling to the
‘Memory’ qubit, which in turn is coupled to two other qubits. Charge lines and
flux lines to each qubit permit XY control and frequency tunability
respectively. A common feedline couples to individual dedicated resonators for
each transmon qubit.

The custom-built open-source library sqdtoolz [57] is used to orchestrate the
array of high-precision instruments required to execute control protocols on
our quantum processor. This includes the synchronised generation and timing of
microwave-frequency pulsed signals for control and measurement, a continuous
microwave-frequency pump signal for parametric amplifier operation and a
stable DC signal for frequency tuning purposes.

Characteristic | System qubit | Memory qubit
---|---|---
Symmetry point frequency | 5.97 GHz | 5.03 GHz
Operational frequency | 5.11 GHz | 5.03 GHz
Coupling to system | N/A | 11 MHz
T1 (relaxation time) | 14 $\mu$s | 10 $\mu$s
T2 (dephasing time) | 1.8 $\mu$s | 6 $\mu$s

Table 2: The characteristics of the major interacting transmon qubits in the
experimental processor. The system qubit is operating at a transition
frequency far from its symmetry point, which increases its sensitivity to flux
noise and in turn reduces T2.

Experiments are also performed on the cloud 7-qubit quantum processor
ibm_perth, as provided by IBM Quantum [37]. This system offers increased
coherence times and gate fidelities at the expense of frequency control.
Transmon characteristics can be found in Table 3.

Characteristic | System qubit | Memory qubit
---|---|---
Qubit frequency | 5.16 GHz | 4.98 GHz
Coupling to system | N/A | 3 MHz
T1 (relaxation time) | 207 $\mu$s | 102 $\mu$s
T2 (dephasing time) | 295 $\mu$s | 113 $\mu$s

Table 3: The characteristics of the major interacting transmon qubits in the
IBM cloud processor. For this system, we have no control over frequency and
therefore no control over interaction strength.

### Error estimation

The credible intervals reported in Table 1 were calculated using the non-
parametric resampling technique of bootstrapping. Bootstrapping a standard
deviation for some metric $\hat{\theta}(\hat{X})$ relies on the assumption
that the relative frequency of the measured $\hat{X}$ is representative of the
underlying distribution which $\hat{X}$ is drawn from. In our case, $\hat{X}$
is a measured set of shot probabilities, which may look something like
$\\{$00: 0.7773, 10: 0.1738, 01: 0.0421, 11: 0.0068$\\}$, where, 00 for
example corresponds to the probability of getting the plus state in both
measurements at $B$ and $C$. We resample 8000 shots from a multinomial
distribution with these underlying relative frequencies, over 1000 trials,
creating a set of bootstrapped results $\hat{X}_{i}$. The simulated
distribution of the metric of interest is then evaluated as
$\hat{\theta}(\hat{X}_{i})$ and twice the standard deviation of this
distribution is used to calculate the 95% credible interval of the reported
metric $\hat{\theta}(\hat{X})$. This technique is designed to account for shot
noise in the finite-sample estimates of observables, and does not account for
systematic error at the hardware level (gate miscalibration, readout
misclassification). Finally, we rely on bias correction in the bootstrapping
process; the sample mean of the metric $\hat{\theta}(\hat{X}_{i})$ was in
general not equal to the measured $\hat{\theta}(\hat{X})$ (sometimes outside
of the credible interval).

### Simulated process

Figure 4: A three-time quantum process involving state preparation and two
gates acting on both systems.

In both platforms used, the qubit we probe belongs to a larger multi-qubit
array. The particular qubits we choose each have one nearest-neighbour qubit,
with which they are directly coupled, as well as additional distant qubits
(two or more degrees of separation), with which interactions will be weaker.
We therefore take the simplest model of our environment to be that of the
single coupled qubit, as shown in Fig. 4. This interaction can be described by
the Hamiltonian [58]

$H=\frac{\omega_{1}}{2}(\mathds{1}\otimes\sigma_{z})+\frac{\omega_{2}}{2}(\sigma_{z}\otimes\mathds{1})+g_{12}(\sigma_{+}\otimes\sigma_{-}+\sigma_{-}\otimes\sigma_{+}),$
(7)

where $\omega_{1,2}$ are the resonance frequencies of the system and memory
qubit, respectively, and $\sigma_{+,-}$ are the creation and annihilation
operators, respectively. Then the free evolution of system-memory is expressed
as $U=\exp{-iHt}$, where $t$ is the time of the evolution.

To write the process matrix for our system, we start with writing the terms:
the initial state of the environment $\rho^{E_{1}}$, the channel
$T_{1}^{A_{O}E_{1}^{T}B_{I}E2}$ and the second channel
$T_{2}^{B_{O}E_{2}^{T}C_{I}E_{3}}$ (see Fig 4), where the superscript $T$ on a
system denotes partial transposition for that system. To obtain the process
matrix, we combine the terms and trace out all the environment systems

$W=\text{Tr}_{E_{1},E_{2},E_{3}}(\rho^{E_{1}}\otimes
T_{1}^{A_{O}E_{1}^{T}B_{I}E2}\otimes T_{2}^{B_{O}E_{2}^{T}C_{I}E_{3}}),$ (8)

where $\rho^{E_{1}}=|{0}\rangle$ because all qubits on the chip are
initialised to the ground state,
$T_{1}^{A_{O}E_{1}B_{I}E2}=[[U]]^{A_{O}E_{1}B_{I}E2}$ and
$T_{2}^{B_{O}E_{2}^{T}C_{I}E_{3}}=[[U]]^{B_{O}E_{2}C_{I}E_{3}}$, and $[[U]]$
is the Choi form of $U$, without the partial transposition.

_Choi matrix definition:_ The Choi matrix $M^{A_{I}A_{O}}\in{\cal L}({\cal
H}^{A_{I}}\otimes{\cal H}^{A_{O}})$, isomorphic to a CP map ${\cal
M}^{A}:{\cal L}({\cal H}^{A_{I}})\rightarrow{\cal L}({\cal H}^{A_{O}})$ is
defined as $M^{A_{I}A_{O}}:=[{\cal I}\otimes{\cal
M}({{|{\mathds{1}}\rangle\\!\rangle}\\!{\langle\\!\langle{\mathds{1}}|}}{\mathds{1}})]^{T}$,
where $\cal I$ is the identity map,
${|{\mathds{1}}\rangle\\!\rangle}=\sum_{j=1}^{d_{A_{I}}}|{jj}\rangle\in{\cal
H}^{A_{I}}\otimes{\cal H}^{A_{I}}$, $\\{|{j}\rangle\\}^{d_{A_{I}}}_{j=1}$ is
an orthonormal basis on ${\cal H}^{A_{I}}$ and ${T}$ denotes matrix
transposition in that basis and some basis of ${\cal H}^{A_{O}}$.

We see that the free evolution, U, depends on the parameters of the
Hamiltonian, $\omega_{1}$, $\omega_{2}$, $g_{1,2}$, and the time $t$. Hence,
our process matrix will also depend on these parameters,
$W(\omega_{1},\omega_{2},g_{1,2},t_{1},t_{2})$. To simulate the process on uq,
experimental parameters were measured to be as given in Table 2, and we scan
$t_{1}$ and $t_{2}$ from 85 to 115 ns, capturing the 97 ns process conducted
in the experiment. For the process on ibm_perth, the parameters were fixed at
the values available in Table 3, and we scan $t_{1}$ and $t_{2}$ between 17
and 33 ns, capturing the permutations of $t_{1}=[21.33,24.89,28.44]$ ns and
$t_{2}=[21.33,24.89,28.44]$ ns conducted in experiment. The ibm_perth process
simulation yielded the results in Fig. 5.

Figure 5: Simulated (a) $\sqrt{JSD}$ and (b) $\mathcal{N}$ values for the
ibm_perth processor over a range of free evolution times $t_{1}$ and $t_{2}$.

According to these results, the experimental runs where evolution times are
$t_{1}=t_{2}=24.89$ ns should maximise both metrics of non-Markovianity on
ibm_perth, while other settings ($t_{1}=21.33$ ns and $t_{2}=28.44$ ns) should
be close to fully suppressing non-Markovian dynamics. In practice, both
metrics were observed to be smeared across the free evolution times scanned.
This is not a particularly surprising result when contrasting the simplicity
of the coupled-qubit environment described in Eq. 7 to the complexity of
current-day quantum hardware, including but not limited to the extra qubits on
each chip.

|  | max$({\sqrt{JSD}})$ | max$({\mathcal{N}})$
---|---|---|---
uq | measured | 0.504 | 0.087
| simulated | 0.270 | 0.091
ibm_perth | measured | 0.206 | 0.024
| simulated | 0.031 | 0.001
Table 4: The maximum values observed for the square root Jensen-Shannon
divergence ($\sqrt{JSD}$) and negativity ($\mathcal{N}$) across the
experimentally measured (across the experimental runs, see Table 1) and
simulated (at the peak of a ($t_{1}$, $t_{2}$) oscillation, see Fig. 5) multi-
time processes.

A point of some interest is the comparison of the maximum $\sqrt{JSD}$ and
$\mathcal{N}$ predicted in simulation (across a reasonable ($t_{1}$, $t_{2}$)
range) to the values observed in experiment. A summary of the these are
provided in Table 4. The qualitative structure of the simulated values match
the results presented in this work, and are reasonable quantitatively on the
uq processor. Both metrics are heavily underestimated on the ibm_perth
processor. This is not an unusual result; the Hamiltonian presented in Eq 7 is
certainly missing several terms, both characterisable (e.g. other qubits on
the chip) and unknown. Hence, although the naive model is a reasonable first-
order approximation, our results show that we can probe more complex dynamics
than we can simulate, making our methods a powerful tool to detect non-
Markovian noise.

## References

  * [1] J. Morris, F. A. Pollock, and K. Modi, Quantifying non-markovian memory in a superconducting quantum computer.
  * McEwen  _et al._ [2021] M. McEwen _et al._ , Removing leakage-induced correlated errors in superconducting quantum error correction, Nature Communications 12, 10.1038/s41467-021-21982-y (2021).
  * Young  _et al._ [2020] K. Young, S. Bartlett, R. J. Blume-Kohout, J. K. Gamble, D. Lobser, P. Maunz, E. Nielsen, T. J. Proctor, M. Revelle, and K. M. Rudinger,  _Diagnosing and Destroying Non-Markovian Noise_, Tech. Rep. (United States, 2020).
  * Piilo  _et al._ [2008] J. Piilo, S. Maniscalco, K. Härkönen, and K.-A. Suominen, Non-markovian quantum jumps, Phys. Rev. Lett. 100, 180402 (2008).
  * Breuer  _et al._ [2009] H.-P. Breuer, E.-M. Laine, and J. Piilo, Measure for the degree of non-markovian behavior of quantum processes in open systems, Phys. Rev. Lett. 103, 210401 (2009).
  * Rivas  _et al._ [2010] A. Rivas, S. F. Huelga, and M. B. Plenio, Entanglement and non-markovianity of quantum evolutions, Phys. Rev. Lett. 105, 050403 (2010).
  * Rivas  _et al._ [2014] Á. Rivas, S. F. Huelga, and M. B. Plenio, Quantum non-markovianity: characterization, quantification and detection, Rep. Prog. Phys. 77, 094001 (2014).
  * Breuer  _et al._ [2016] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Colloquium: Non-markovian dynamics in open quantum systems, Rev. Mod. Phys. 88, 021002 (2016).
  * Chruściński  _et al._ [2011] D. Chruściński, A. Kossakowski, and A. Rivas, Measures of non-markovianity: Divisibility versus backflow of information, Phys. Rev. A 83, 052128 (2011).
  * de Vega and Alonso [2017] I. de Vega and D. Alonso, Dynamics of non-markovian open quantum systems, Reviews of Modern Physics 89, 015001 (2017).
  * Takeda  _et al._ [2018] K. Takeda _et al._ , Optimized electrical control of a si/sige spin qubit in the presence of an induced frequency shift, npj Quantum Information 4, 54 (2018).
  * Freer  _et al._ [2017] S. Freer _et al._ , A single-atom quantum memory in silicon, Quantum Science and Technology 2, 015009 (2017).
  * Savytskyy  _et al._ [2023] R. Savytskyy _et al._ , An electrically driven single-atom “flip-flop”qubit, Science Advances 9, eadd9408 (2023).
  * Zwerver  _et al._ [2022] A. M. J. Zwerver _et al._ , Qubits made by advanced semiconductor manufacturing, Nature Electronics 5, 184 (2022).
  * Philips  _et al._ [2022] S. G. J. Philips _et al._ , Universal control of a six-qubit quantum processor in silicon, Nature 609, 919 (2022).
  * Oreshkov  _et al._ [2012] O. Oreshkov, F. Costa, and Č. Brukner, Quantum correlations with no causal order, Nat. Commun. 3, 1092 (2012), arXiv:1105.4464 [quant-ph] .
  * Oreshkov and Giarmatzi [2016] O. Oreshkov and C. Giarmatzi, Causal and causally separable processes, New J. Phys. 18, 093020 (2016).
  * Giarmatzi and Costa [2021] C. Giarmatzi and F. Costa, Witnessing quantum memory in non-Markovian processes, Quantum 5, 440 (2021).
  * Shrapnel  _et al._ [2018a] S. Shrapnel, F. Costa, and G. Milburn, Quantum markovianity as a supervised learning task, International Journal of Quantum Information 16, 1840010 (2018a).
  * Taranto  _et al._ [2021] P. Taranto, F. A. Pollock, and K. Modi, Non-markovian memory strength bounds quantum process recoverability, npj Quantum Information 7, 149 (2021).
  * Fux  _et al._ [2021] G. E. Fux, E. P. Butler, P. R. Eastham, B. W. Lovett, and J. Keeling, Efficient exploration of hamiltonian parameter space for optimal control of non-markovian open quantum systems, Phys. Rev. Lett. 126, 200401 (2021).
  * Figueroa-Romero  _et al._ [2021] P. Figueroa-Romero, K. Modi, R. J. Harris, T. M. Stace, and M.-H. Hsieh, Randomized benchmarking for non-markovian noise, PRX Quantum 2, 040351 (2021).
  * Milz  _et al._ [2018] S. Milz, F. A. Pollock, and K. Modi, Reconstructing non-markovian quantum dynamics with limited control, Phys. Rev. A 98, 012108 (2018).
  * Muñoz  _et al._ [2022] R. N. Muñoz, L. Frazer, G. Yuan, P. Mulvaney, F. A. Pollock, and K. Modi, Memory in quantum dot blinking, Phys. Rev. E 106, 014127 (2022).
  * Goswami  _et al._ [2021] K. Goswami, C. Giarmatzi, C. Monterola, S. Shrapnel, J. Romero, and F. Costa, Experimental characterization of a non-markovian quantum process, Phys. Rev. A 104, 022432 (2021).
  * Morris  _et al._ [2022] J. Morris, F. A. Pollock, and K. Modi, Quantifying non-markovian memory in a superconducting quantum computer, Open Systems & Information Dynamics 29, 10.1142/S123016122250007X (2022).
  * White  _et al._ [2020] G. A. L. White, C. D. Hill, F. A. Pollock, L. C. L. Hollenberg, and K. Modi, Demonstration of non-markovian process characterisation and control on a quantum processor, Nature Communications 11, 10.1038/s41467-020-20113-3 (2020).
  * Xiang  _et al._ [2021] L. Xiang _et al._ , Quantify the non-markovian process with intervening projections in a superconducting processor (2021).
  * White  _et al._ [2021] G. A. L. White, F. A. Pollock, L. C. L. Hollenberg, C. D. Hill, and K. Modi, From many-body to many-time physics (2021).
  * White  _et al._ [2022] G. White, F. Pollock, L. Hollenberg, K. Modi, and C. Hill, Non-markovian quantum process tomography, PRX Quantum 3, 10.1103/prxquantum.3.020344 (2022).
  * Chiribella  _et al._ [2013] G. Chiribella, G. M. D’Ariano, P. Perinotti, and B. Valiron, Quantum computations without definite causal structure, Phys. Rev. A 88, 022318 (2013), arXiv:0912.0195 [quant-ph] .
  * Araújo  _et al._ [2015] M. Araújo, C. Branciard, F. Costa, A. Feix, C. Giarmatzi, and Č. Brukner, Witnessing causal nonseparability, New J. Phys. 17, 102001 (2015).
  * Procopio  _et al._ [2015] L. M. Procopio, A. Moqanaki, M. Araújo, F. Costa, I. Alonso Calafell, E. G. Dowd, D. R. Hamel, L. A. Rozema, Č. Brukner, and P. Walther, Experimental superposition of orders of quantum gates, Nature Communications 6, 10.1038/ncomms8913 (2015).
  * Rubino  _et al._ [2017] G. Rubino, L. A. Rozema, A. Feix, M. Araújo, J. M. Zeuner, L. M. Procopio, Č. Brukner, and P. Walther, Experimental verification of an indefinite causal order, Science Advances 3, e1602589 (2017).
  * Goswami  _et al._ [2018] K. Goswami, C. Giarmatzi, M. Kewming, F. Costa, C. Branciard, J. Romero, and A. White, Indefinite causal order in a quantum switch, Physical Review Letters 121, 10.1103/physrevlett.121.090503 (2018).
  * Antesberger  _et al._ [2023] M. Antesberger, M. T. Quintino, P. Walther, and L. A. Rozema, Higher-order process matrix tomography of a passively-stable quantum switch (2023), arXiv:2305.19386 [quant-ph] .
  * IBM Quantum [2023] IBM Quantum,  (2023).
  * Chiribella  _et al._ [2009] G. Chiribella, G. M. D’Ariano, and P. Perinotti, Theoretical framework for quantum networks, Phys. Rev. A 80, 022339 (2009), arXiv:0904.4483 [quant-ph] .
  * Costa and Shrapnel [2016] F. Costa and S. Shrapnel, Quantum causal modelling, New J. Phys. 18, 063032 (2016).
  * Gutoski and Watrous [2006] G. Gutoski and J. Watrous, Toward a general theory of quantum games, in  _In Proceedings of 39th ACM STOC_ (2006) pp. 565–574, quant-ph/0611234 .
  * Pollock  _et al._ [2018] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, Non-markovian quantum processes: Complete framework and efficient characterization, Physical Review A 97, 012127 (2018).
  * Lindblad [1979] G. Lindblad, Non-markovian quantum stochastic processes and their entropy, Comm. Math. Phys. 65, 281 (1979).
  * Jamiołkowski [1972] A. Jamiołkowski, Linear transformations which preserve trace and positive semidefiniteness of operators, Rep. Math. Phys 3, 275 (1972).
  * Choi [1975] M.-D. Choi, Completely positive linear maps on complex matrices, Linear Algebra and its Applications 10, 285 (1975).
  * Heinosaari and Ziman [2011] T. Heinosaari and M. Ziman,  _The Mathematical Language of Quantum Theory: From Uncertainty to Entanglement_ (Cambridge University Press, 2011).
  * Davies and Lewis [1970] E. B. Davies and J. T. Lewis, An operational approach to quantum probability, Communications in Mathematical Physics 17, 239 (1970).
  * Shrapnel  _et al._ [2018b] S. Shrapnel, F. Costa, and G. Milburn, Updating the born rule, New Journal of Physics 20, 053010 (2018b).
  * Giarmatzi and Costa [2018] C. Giarmatzi and F. Costa, A quantum causal discovery algorithm, npj Quantum Information 4, 17 (2018).
  * Settimo  _et al._ [2022] F. Settimo, H.-P. Breuer, and B. Vacchini, Entropic and trace-distance-based measures of non-markovianity, Physical Review A 106, 10.1103/physreva.106.042212 (2022).
  * Jouzdani and Mucciolo [2014] P. Jouzdani and E. R. Mucciolo, Numerical evaluation of the fidelity error threshold for the surface code, Phys. Rev. A 90, 012315 (2014).
  * Paz-Silva  _et al._ [2017] G. A. Paz-Silva, L. M. Norris, and L. Viola, Multiqubit spectroscopy of gaussian quantum noise, Physical Review A 95, 10.1103/physreva.95.022121 (2017).
  * Beaudoin  _et al._ [2018] F. Beaudoin, L. M. Norris, and L. Viola, Ramsey interferometry in correlated quantum noise environments, Phys. Rev. A 98, 020102 (2018).
  * Paz-Silva  _et al._ [2019] G. A. Paz-Silva, L. M. Norris, F. Beaudoin, and L. Viola, Extending comb-based spectral estimation to multiaxis quantum noise, Physical Review A 100, 10.1103/physreva.100.042334 (2019).
  * Vidal and Werner [2002] G. Vidal and R. F. Werner, Computable measure of entanglement, Physical Review A 65, 10.1103/physreva.65.032314 (2002).
  * Lin [1991] J. Lin, Divergence measures based on the shannon entropy, IEEE Transactions on Information Theory 37, 145 (1991).
  * Quantware [2023] Quantware,  (2023).
  * Pakkiam  _et al._ [2023] P. Pakkiam, R. Beriwal, and J. Dearlove, Sqdtoolz (2023).
  * Blais  _et al._ [2007] A. Blais, J. Gambetta, A. Wallraff, D. I. Schuster, S. M. Girvin, M. H. Devoret, and R. J. Schoelkopf, Quantum-information processing with circuit quantum electrodynamics, Physical Review A 75, 10.1103/physreva.75.032329 (2007).


Sk-69 175 | – | $0.19\pm 0.26$ | – | $-0.20\pm 0.28$ | $-0.05\pm 0.19$ | – | – | $0.61\pm 0.18$ | – | –
Sk-69 246 | $0.26\pm 0.15$ | $0.14\pm 0.17$ | $0.04\pm 0.15$ | $-0.09\pm 0.14$ | $-0.08\pm 0.10$ | $0.18\pm 0.16$ | – | $0.11\pm 0.18$ | – | –
Sk-69 279 | $0.38\pm 0.16$ | $0.14\pm 0.15$ | $0.08\pm 0.10$ | $-0.11\pm 0.12$ | $0.10\pm 0.08$ | – | – | – | – | –
Sk-70 115 | $0.21\pm 0.14$ | $0.25\pm 0.15$ | $0.02\pm 0.12$ | $-0.09\pm 0.12$ | $-0.05\pm 0.09$ | $0.25\pm 0.14$ | – | $0.25\pm 0.10$ | – | $-0.23\pm 0.11$
Sk-70 79 | $0.05\pm 0.16$ | – | – | $-0.15\pm 0.20$ | $-0.05\pm 0.13$ | – | $0.36\pm 0.21$ | $-0.09\pm 0.19$ | – | –
Sk-71 45 | $0.26\pm 0.20$ | $0.27\pm 0.16$ | $0.04\pm 0.10$ | $-0.08\pm 0.10$ | $0.04\pm 0.08$ | – | – | – | – | $-0.14\pm 0.10$
Sk-71 50 | $-0.04\pm 0.23$ | $0.12\pm 0.12$ | $0.02\pm 0.09$ | $-0.12\pm 0.09$ | $-0.05\pm 0.07$ | – | – | $0.19\pm 0.15$ | $0.10\pm 0.20$ | –
Table 4: Same as Table 3, but for the SMC. ID | [Mg/Fe]nucl | [Si/Fe]nucl | [Cr/Fe]nucl | [Ni/Fe]nucl | [Zn/Fe]nucl | [Ti/Fe]nucl | [S/Fe]nucl | [Mn/Fe]nucl
---|---|---|---|---|---|---|---|---
AzV 18 | $0.38\pm 0.12$ | – | $0.10\pm 0.16$ | $-0.07\pm 0.14$ | $-0.07\pm 0.12$ | $0.38\pm 0.16$ | – | $-0.47\pm 0.12$
AzV 26 | $0.14\pm 0.09$ | $0.29\pm 0.09$ | $0.04\pm 0.10$ | $-0.05\pm 0.09$ | $0.03\pm 0.08$ | $0.23\pm 0.10$ | – | $-0.37\pm 0.08$
AzV 47 | $0.28\pm 0.15$ | $0.22\pm 0.10$ | $0.05\pm 0.12$ | $-0.08\pm 0.11$ | $0.02\pm 0.10$ | $0.28\pm 0.12$ | $0.35\pm 0.11$ | $-0.35\pm 0.10$
AzV 78 | $0.24\pm 0.19$ | – | $0.30\pm 0.21$ | $-0.08\pm 0.20$ | – | $0.31\pm 0.23$ | – | $-0.35\pm 0.18$
AzV 80 | $0.35\pm 0.16$ | – | $0.04\pm 0.18$ | $-0.06\pm 0.18$ | $-0.15\pm 0.13$ | $0.38\pm 0.20$ | – | $-0.43\pm 0.15$
AzV 95 | $0.19\pm 0.17$ | $0.32\pm 0.10$ | $0.07\pm 0.12$ | $-0.10\pm 0.11$ | $-0.03\pm 0.10$ | $0.30\pm 0.13$ | $0.31\pm 0.11$ | $-0.37\pm 0.10$
AzV 104 | $0.43\pm 0.23$ | $0.38\pm 0.24$ | $-0.04\pm 0.30$ | $-0.07\pm 0.28$ | $-0.07\pm 0.25$ | – | $0.32\pm 0.21$ | –
AzV 207 | – | $0.29\pm 0.41$ | $0.05\pm 0.46$ | – | $-0.15\pm 0.31$ | $0.34\pm 0.51$ | – | $-0.40\pm 0.39$
AzV 216 | $0.39\pm 0.12$ | – | $0.02\pm 0.13$ | $-0.03\pm 0.12$ | $-0.01\pm 0.11$ | – | – | –
AzV 229 | $0.32\pm 0.12$ | $0.39\pm 0.09$ | $0.11\pm 0.13$ | $0.10\pm 0.10$ | $-0.14\pm 0.10$ | $0.42\pm 0.11$ | $0.41\pm 0.08$ | $-0.29\pm 0.09$
AzV 242 | $0.34\pm 0.12$ | $0.36\pm 0.09$ | $-0.09\pm 0.23$ | $0.01\pm 0.10$ | $0.31\pm 0.19$ | $0.38\pm 0.10$ | $0.39\pm 0.09$ | $-0.42\pm 0.09$
AzV 321 | – | $0.32\pm 0.12$ | – | $0.15\pm 0.16$ | $0.02\pm 0.13$ | $0.41\pm 0.11$ | $0.42\pm 0.09$ | $-0.31\pm 0.15$
AzV 332 | – | $0.39\pm 0.10$ | $0.02\pm 0.11$ | $-0.14\pm 0.15$ | $0.02\pm 0.09$ | – | $0.41\pm 0.10$ | $-0.22\pm 0.10$
AzV 388 | $0.60\pm 0.16$ | $0.38\pm 0.12$ | – | – | $-0.08\pm 0.11$ | $0.45\pm 0.17$ | $0.44\pm 0.10$ | $-0.37\pm 0.14$
AzV 456 | – | $0.41\pm 0.23$ | – | – | $0.41\pm 0.20$ | $0.41\pm 0.32$ | – | $-0.25\pm 0.25$
AzV 476 | – | – | $0.03\pm 0.14$ | $-0.10\pm 0.15$ | $0.23\pm 0.19$ | $0.24\pm 0.15$ | – | $-0.33\pm 0.11$
Sk 190 | – | $0.61\pm 0.15$ | – | $0.09\pm 0.19$ | – | $0.49\pm 0.20$ | $0.47\pm 0.12$ | –
Sk 191 | – | $0.37\pm 0.15$ | – | – | – | $0.41\pm 0.20$ | $0.43\pm 0.14$ | $-0.29\pm 0.17$
AzV 238 | – | $0.35\pm 0.06$ | $0.09\pm 0.11$ | – | $0.07\pm 0.07$ | – | – | –
AzV 327 | – | $0.66\pm 0.17$ | $-1.20\pm 1.04$ | – | $0.15\pm 0.24$ | – | – | –
Sk 9 | – | $0.20\pm 2.10$ | $0.00\pm 2.10$ | – | $1.50\pm 2.10$ | – | – | –
Sk 116 | – | $0.30\pm 0.06$ | $0.00\pm 0.06$ | – | $0.03\pm 0.05$ | – | – | –

##  Appendix A On the use of [Zn/Fe] as dust tracer

The main processes that shape the abundance patterns in the neutral ISM are
visualized in Figure 3. The different metals are used to determine the
chemical properties of the ISM (see Sect. B). The $x$-axis there refers to a
general refractory index, which is a way to measure how heavily each element
is typically depleted into dust. In our analysis, we use the $B2_{X}$
coefficients for this purpose, which are determined by De Cia et al. (2016)
and Konstantopoulou et al. (2022) preferentially using [Zn/Fe] to trace the
amount of dust in the neutral ISM. The fact that the observed abundances in
the neutral ISM overall line up linearly in the observed abundance patterns
both in the Milky Way (De Cia et al. 2021) and the Magellanic Clouds (Figs. 4
and 5) indicates that the values of $B2_{X}$ are, to a first approximation,
reliable. The metallicities resulting from the use of the $B2_{X}$
coefficients (i.e. the ’relative method’ based on the relative abundances of
metals) do overall agree with those resulting from an analogous but
independent method that does not make preferential use of [Zn/Fe] (De Cia et
al. 2021), which is based on $F*$ to trace the strength of dust depletion
(Jenkins 2009).

The reasons for the choice of [Zn/Fe] to trace dust is that Zn and Fe are
depleted into dust grains in very different ways, they are often well
observed, and their nucleosynthesis overall follow each other in the
metallicity regime of interest. This last point is crucial, specially given
that stellar observations of [Zn/Fe] can show large scatter and deviations
from solar values.

The overall distribution of [Zn/Fe] with [Fe/H] in stars in the Milky Way is
rather flat and near solar, in the metallicity range $-2.0<$ [Fe/H] $<+0.5$
and strongly deviates from solar at lower and higher metallicities (e.g.
Sneden et al. 1991; Primas et al. 2000; Nissen et al. 2007; Saito et al.
2009). Recently, Sitnova et al. (2022) study the distribution of [Zn/Fe] with
[Fe/H] in Milky Way stars (in the halo, thick and thin disk) in the
metallicity regime $-2.0<$ [Fe/H] $<+0.2$, using non-LTE measurements, and
find [Zn/Fe] enhancement of up to 0.2 dex in low-metallicity stars (in the
halo) and under-abundance of down to $-0.2$ dex in high-metallicity stars (in
the thin disk). We explore this possibility further in Sect. 5. Stellar
abundances of stars younger than a few Gyrs in the thin disk of the Milky Way
show also a large scatter in [Zn/Fe], and values down to $-0.2$ dex at solar
and super-solar metallicities (Delgado Mena et al. 2019). The stellar [Zn/Fe]
in halo stars in the Milky Way is close to zero for metallicities in the range
$-2.0<$ [Fe/H] $<-1.0$ and increases to higher values (by 0.1–02 dex) at lower
metallicities, and up to $+0.5$ dex at [Fe/H] $\sim-3$ (Nissen et al. 2007).
In the bulge of the Milky Way, a gas-poor environment with no star-formation
activity, the values of [Zn/Fe] go down to $-0.6$ dex at super-solar
metallicities (Barbuy et al. 2015; Duffau et al. 2017; da Silveira et al.
2018). In nearby dwarf spheroidal galaxies stellar measurements of [Zn/Fe]
also show a large scatter, with value down to about $-1$ dex (Skúladóttir et
al. 2017). Dwarf spheroidals are fundamentally different from DLAs because
they do not currently have gas and should have lost their gas content by
$z\sim 2.3$ (Salvadori & Ferrara 2012).

In the neutral ISM in local and distant galaxies there is no evidence so far
for large intrinsic (nucleosynthetic) deviations of [Zn/Fe]. The existence of
depletion sequences with scatter $<0.2$ dex (De Cia et al. 2016;
Konstantopoulou et al. 2022) limits the possible [Zn/Fe] deviations due to
nucleosynthesis. We speculate that any Zn deviations due to nucleosynthesis
should be lower than $0.2$ dex. In particular, in DLAs we observe values down
to [Zn/Fe] $=0$, and never strongly negative values. In this regime, dust
depletion is negligible. When [Zn/Fe] $=0$, DLAs show enhancement of $\alpha$
elements of about 0.3 dex, and under-abundance of Mn by about 0.3 dex (De Cia
et al. 2016), similar to the [$X$/Fe] expected from observations in the Local
Group and chemical evolution of galaxies (see Sect. 4.2). In addition, [Zn/Fe]
correlates well with [Si/Ti], [O/Si], and other dust tracers. This excludes
strong deviations of Zn from Fe due to nucleosynthesis. More details about why
[Zn/Fe] can be considered a reliable dust tracer, to a first approximation,
are discussed in De Cia (2018). Future studies of dust depletion based on
relative abundances may be able to expand the framework and avoid any
preferential use of [Zn/Fe] as a dust tracer.

##  Appendix B Choosing different sets of metals in the analysis of the
abundance patterns

We compare the results obtained from different choices of metals as reference
for the linear fit to the abundance patterns. First, we use all available
metals except Mn, and derive the metallicity $[{\rm M/H}]_{\rm
tot,all\,metals}$ and strength of dust depletion $[{\rm Zn/Fe}]_{\rm
fit,all\,metals}$. This basic approach already shows $\alpha$-element
enhancement and Mn under-abundance in several systems. The choice of the
exclusion of Mn is driven by the Mn under-abundance observed in the ISM of the
Magellanic Clouds (Jenkins & Wallerstein 2017; Konstantopoulou et al. 2022).
Indeed, including Mn in the fit to the abundance patterns produces sometimes
non-physical negative slopes for cases when not many relatively volatile
elements are available. Second, we consider the fit to all metals but
excluding all the $\alpha$ elements and Mn, to avoid misinterpreting the
effects of dust, nucleosynthesis, and metallicity. We denote the resulting
metallicity and strength of depletion as $[{\rm M/H}]_{\rm tot,non-\alpha}$
and $[{\rm Zn/Fe}]_{\rm fit,non-\alpha}$, respectively, and use them as
reference. In the three cases where neither Zn nor other non-$\alpha$ volatile
elements are constrained (namely, SK-67 14, SK 190 and SK 191), we consider
the fit to all metals including a posteriori an additional 0.32 dex999The 0.32
dex is chosen here to be in-between the typical values of $\alpha$-element
enhancement in the Milky Way ($\sim 0.35$ dex McWilliam 1997) and DLAs ($\sim
0.3$ dex De Cia et al. 2016). in the uncertainty on their total metallicity to
account for the paucity of data and potential nucleosynthesis effects causing
over- or underestimate of the metallicity. Third, we include only the $\alpha$
elements, with the exception of O, and linearly fit the abundance patterns
only for systems that have a constrained measurement of Ti, to ensure a large
enough dynamical range in the refractory index. This is possible for a limited
number of systems (the golden sample, with larger symbols in Figs. 15 and
16.). We denote this resulting metallicity as $[{\rm M/H}]_{\rm tot,\alpha}$
and the strength of dust depletion as $[{\rm Zn/Fe}]_{\rm fit,\alpha}$. Oxygen
is observed in only three LMC systems, and its exclusion from the fit is
driven by the potential presence of ISM mixture at least in one of them, which
we discuss in Section 4.5.

The reference values $[{\rm Zn/Fe}]_{\rm fit}$ and $[{\rm M/H}]_{\rm tot}$ are
those from the optimal approach, which is described in Sect. 3. Here we
compare these with the results from the same analysis but using three
different sets of metals: only the $\alpha$ elements, only the non-$\alpha$
elements, and all the metals. Figure 15 shows a comparison between the
strength of depletion $[{\rm Zn/Fe}]_{\rm fit}$ (i.e. the slope of the linear
fits to the abundance patterns) for the optimal case with the case where only
the non-$\alpha$ elements were used (top panel), or all metals (mid panel), or
the observed [Zn/Fe] (bottom panel). The optimal $[{\rm Zn/Fe}]_{\rm fit}$
correlates reasonably well with all the three quantities, with notable
exceptions.

Figure 15: Comparison between the overall strength of the dust depletion
[Zn/Fe]fit from the fit to the abundance patterns for the optimal choice of
metals (see Sect. 3) with the observed [Zn/Fe] (bottom panel), the fit to all
metals (middle panel), and the fit to only the non-$\alpha$ elements. Orange
triangles show the measurements for the LMC and blue diamonds those for the
SMC. The dotted curve is the $x=y$ line. Larger symbols highlight the golden
sample. Figure 16: Comparison between the total metallicity [M/H]tot from the
fit to the abundance patterns for the optimal choice of metals (see Sect. 3)
with the observed [Zn/H] (bottom panel), the fit to all metals (third panel),
and the fit to only the non-$\alpha$ elements (second panel) and the fit to
only the $\alpha$ elements (top panel). Orange triangles show the measurements
for the LMC and blue diamonds those for the SMC. The dotted curve is the $x=y$
line. The dashed line on the top panel is not a fit to the data, but marks the
$x=y+0.35$ dex line, the enhancement in the abundances of $\alpha$ elements.
Larger symbols highlight the golden sample.

In a similar way, Figure 16 shows a comparison between the total metallicity
$[{\rm M/H}]_{\rm tot}$ (i.e. the normalization of the linear fits to the
abundance patterns) for the optimal case and the case where only the $\alpha$
elements were used (top panel), or only the non-$\alpha$ elements (second
panel), or all metals (third panel), or the observed [Zn/H] (bottom panel).
Clearly, the observed [Zn/H] underestimates the actual metallicity, and this
is likely because of dust depletion. The values of $[{\rm M/H}]_{\rm
tot,non-\alpha}$ and $[{\rm M/H}]_{\rm tot,all\,metals}$ correlate well with
the optimal total metallicity, however with some deviations. The values of
$[{\rm M/H}]_{\rm tot,\alpha}$ are mostly higher than the optimal metallicity,
and this is likely due to $\alpha$-element enhancement. The dashed purple line
on the top panel of Fig. 16 is not a fit to the data, but the $x=y$ line
increased by 0.35 dex, as a proxy for the overall level of $\alpha$-element
enhancement in the Milky Way (e.g. Lambert 1987). Overall, because of the
likely $\alpha$-element enhancement, it is not optimal to choose derive the
metallicity from the fit to the abundance patterns of the $\alpha$ elements
only. Indeed, we refer to the total metallicity $[{\rm M/H}]_{\rm tot}$ as the
metallicity of the Fe-peak elements, after dust correction.


*[enumerate,1]label=0) tee

# Sample-Based Output-Feedback Navigation with Bearing Measurements

Mahroo Bahreinian1, Marc Mitjans2, Roy Xing3 and Roberto Tron4 This work was
supported by ONR MURI N00014-19-1-2571 “Neuro-Autonomy: Neuroscience-Inspired
Perception, Navigation, and Spatial Awareness”1Mahroo Bahreinian is with
Division of Systems Engineering at Boston University, Boston, MA, 02215 USA.
Email<EMAIL_ADDRESS>Mitjans is with Department of Mechanical
Engineering at Boston University, Boston, MA, 02215 USA. M. Mitjans is
additionally supported by ”la Caixa” Foundation fellowship
LCF/BQ/AA18/11680117<EMAIL_ADDRESS>Roy Xing is with Department of
Mechanical Engineering at Boston University, Boston, MA, 02215 USA.
<EMAIL_ADDRESS>Tron is with Department of Mechanical and System
<EMAIL_ADDRESS>
###### Abstract

We consider the problem of sample-based feedback-based motion planning from
_bearing_ (direction-only) measurements. We build on our previous work that
defines a cell decomposition of the environment using RRT∗, and finds an
output feedback controller to navigate through each cell toward a goal
location using duality, Control Lyapunov and Barrier Functions (CLF, CBF), and
Linear Programming. In this paper, we propose a novel strategy that uses
relative bearing measurements with respect to a set of landmarks in the
environment, as opposed to full relative displacements. The main advantage is
then that the measurements can be obtained using a simple monocular camera. We
test the proposed algorithm in the simulation, and then in an experimental
environment to evaluate the performance of our approach with respect to
practical issues such as mismatches in the dynamical model of the robot, and
measurements acquired with a camera with a limited field of view.

##  I INTRODUCTION

Motion planning is a major research area in the context of mobile robots, as
it deals with the problem of finding a path from an initial state toward a
goal state while avoiding collisions [1, 2, 3]. Traditional path planning
methods focus on finding _single nominal paths_ in a given _known map_ , and
the majority of them makes the implicit assumption that the agent possesses a
lower-level _state feedback_ controller for following such nominal path in the
face of external disturbances and imperfect models. In this paper, we instead
use the alternative approach of synthesizing a set of output-feedback
controllers instead of single a nominal path. Focusing on controllers allows
us to directly consider the measurements (outputs) available to the agent,
instead of assuming full state knowledge (i.e., a perfect localization in the
environment); the main advantage of this approach is that it is intrinsically
robust to disturbances and limited changes in the environment. The main
challenge, however, is that we need to explicitly take into account the type
of measurements that are available. For instance, in this paper we consider
the case of a mobile robot that can use a monocular camera to detect a set of
landmarks (e.g., objects) in the environment. Using image-based information,
the robot can compute the bearing (i.e., the direction vector) with respect to
the landmarks. The goal of this work is then to synthesize controllers that
can use this type of measurements for feedback-based path planning.

_Previous works_. A well known class of techniques for motion planning is
sampling-based methods. Algorithms such as Probabilistic RoadMaps (PRM) [4],
Rapidly exploring Random Trees (RRT∗, [5, 6]) and asymptotically optimal
Rapidly Exploring Random Tree (RRT∗, [7]), have become popular in the last few
years due to their good practical performance, and their probabilistic
completeness [8, 6, 7]; there have also been extensions considering perception
uncertainty [9]. However, these algorithms only provide nominal paths, and
assume that a separate low-level controller exists to generate collision-free
trajectories at run time.

The concept of the sample-based output feedback controller is based on [10,
11] which, the output feedback controller takes as input relative displacement
measurements with respect to a set of landmarks. This algorithm is based on
solving a sequence of robust min-max Linear Programming (LP) problems on the
elements of a convex cell decomposition of the environment and using linear
Control Lyapunov Function (CLF) and Control Barrier Function (CBF)
constraints, to provide stability and safety guarantees, respectively.
Additionally, it extends planning beyond the discrete nominal paths, as
feedback controllers can correct deviations from such paths, and are robust to
discrepancies between the planning and real environment maps. The authors in
[10] assumes that a polyhedral convex cell decomposition of the environment is
available, which greatly reduces its applicability. The authors in [11]
extended upon [10] and introduces a representation of the environment that is
obtained via a modified version of the optimal Rapidly-exploring Random Trees
(RRT∗), combined with a convex cell decomposition of the environment.

In order to observe the environment, modern robots can be equipped with
different kinds of sensors, such as Inertial Measurement Units [12], visual
sensors [13], sonar sensors[14], and contact sensors [15], etc. Among these,
monocular cameras are nearly ubiquitous, given their favorable trade-offs with
respect to size, power consumption, and the richness of information in the
measurements; a peculiarity of this sensor is that, due to perspective
projection, it can provide relatively precise _bearing_ information (the
direction of an object or point in the environment) but not, from a single
image alone, _depth_ information. Many different techniques have been
introduced for monocular navigation and control. An example is visual
servoing, where image-based computations are used to find a control law that
guides the robot toward a home location using bearing alone [16, 17, 18, 19],
or by concurrently estimating depths [20, 21, 22].

_Proposed approach and contributions_. The goal of the present work is to
joint the two threads of work mentioned above: apply the synthesis method of
[10, 11] but with controllers that use as input a set of bearing measurements
between the robot and a set of landmarks. Since we consider a planning
problem, we assume that a map of the environment is available, and that the
positions of the landmarks are known and fixed in the reference frame of this
map; furthermore, we assume that a global compass direction is available (so
that the orientation, but not the position of the robot in the map is known).
The main contribution of this work is then to provide a way to synthesize
bearing-based feedback controllers that can be used to solve the path planning
problem by extending the method of [10, 11] via a dynamic rescaling of the
measurements. This contribution is tested in both simulations and experiments
with real robots.

##  II BACKGROUND

In this section, we review the system dynamics, bearing direction
measurements, CLF and CBF constraints, and the RRT∗ method in the context of
our proposed work.

###  II-A System Dynamics

Consider a driftless control affine system

$\dot{x}=Bu,$ (1)

where $x\in\mathcal{X}\subset n$ denotes the state, $u\in\mathcal{U}\subset
m{}$ is the system input, and $B\in n\times m$ define the linear dynamics of
the system. Note that We assume that the system (1) is controllable, that
$\mathcal{X}$ and $\mathcal{U}$ are polytopic,

$\displaystyle\mathcal{X}=\\{x\mid A_{x}x\leq b_{x}\\},$
$\displaystyle\mathcal{U}=\\{u\mid A_{u}u\leq b_{u}\\},$ (2)

and that $0\in\mathcal{U}$. The assumptions above are more restrictive than
those in, e.g., [11], but are necessary to allow our solution to the bearing-
based control problem.

Note that, in our case, $\mathcal{X}$ will be a convex cell centered around a
sample in the tree (Sec. II-D).

###  II-B Bearing direction measurements

In this work, we assume global compass direction is available (the bearing
directions can be compared in the same frame of reference) and the robot has
only access to bearing direction measurements (see Fig. 1a). The bearing
direction measurements are defined as

$\beta_{i}=d_{i}(x)^{-1}(l_{i}-x),\;\;\;\;i=\\{1,\ldots,N\\},$ (3)

where $N$ is the number of _landmarks_ and $d_{i}$ is the relative
displacement measurement between the position of the robot and landmark
$l_{i}$

$d_{i}(x)=\lVert l_{i}-x\rVert.$ (4)

The location of landmarks ${l}_{i}$ is assumed to be fixed and known at the
planning time. Additionally, we assume at the implementation time, the robot
can measure the bearing direction $\beta_{i}$ between its position $x$ and a
set of landmarks $l_{i}$.

(a)

(b)

(c)

Figure 1: a) Bearing direction measurements. Landmarks are shown by blue
circles and the robot is shown by a black circle. We assume the robot measures
the bearing direct to landmarks and they are shown by blue arrows. The
position of landmarks is known and bearing direction measurements between
landmarks are shown by purple, green, and orange arrows. b) As bearing
measurements are unit vectors, the robot sees all the landmarks within 1 unit
from itself. c) modify bearing measurements by assuming the landmark $l_{2}$
is fixed.

###  II-C Optimal Rapidly-Exploring Random Tree (RRT∗)

The RRT∗ sampling-based algorithm is typically used for single-query path
planning, but that can also be used (as in this paper) to build a tree
$\mathcal{T}=(\mathcal{V},\mathcal{E})$ representing the free configuration
space starting from a given root node. In this work, we build a tree in the
free configuration space by implementing the RRT∗ algorithm, and the only
modification we made to the original RRT∗ algorithm is that we store random
samples that were found to be in collision with an obstacle in the list
$\mathcal{V}_{\text{collision}}$ instead of discarding them; this list is then
returned and used to define the CBF constraints in our algorithm (see Sec.
II-G). RRT∗ is guaranteed to be asymptotically complete and optimal, although
these guarantees do not hold with a finite number of samples.

###  II-D Convex Decomposition of the Environment by Simplified RRT∗ Trees

We start with a tree $\mathcal{T}=(\mathcal{V},\mathcal{E})$ generated by the
traditional RRT∗ algorithm from Sec. II-C. Since the number of samples is
finite, the generated tree is not optimal although it has a large number of
nodes. We simplify the tree to reduce the number of nodes (while keeping all
the samples that are in collision with obstacles) by following the simplified-
RRT∗ algorithm in [11] and denote it as
$\mathcal{T}=(\mathcal{V}_{s},\mathcal{E}_{s})$.

Note that as a consequence of the simplifying steps above, it is still
possible to connect any sample that was discarded from the original RRT∗ to
the simplified tree with a straight line, suggests that the simplified tree
will be a good road-map representation [23] of the free configuration space
reachable from the root (up to the effective resolution given by the original
sampling).

Given the simplified tree $\mathcal{T}_{s}=(\mathcal{V}_{s},\mathcal{E}_{s})$,
for each node $i\in\mathcal{V}_{s}$ in the tree, we define a convex cell
$\mathcal{X}_{ij}$ similar to [11] such that the boundaries of
$\mathcal{X}_{ij}$ are defined as the bisectors hyper-plane between node $i$
and other nodes in the tree except node $j$ which node $j$ is the parent of
node $i$. The polyhedron $\mathcal{X}_{ij}$ is similar to a Voronoi region
[2](See Fig. 2 for an example). Note that $\mathcal{X}_{ij}$ contains all the
points that are closest to $i$ than other vertices in $\mathcal{T}_{s}$, but
it also includes the parent $j$.

Figure 2: The Voronoi-like region $\mathcal{X}_{ij}$ for edge $(i,j)$, and the
corresponding CBF constraints; green points and arrows: other vertices and
edges in the tree $\mathcal{T}_{s}$; red points: collision samples in
$\mathcal{V}_{\text{collision}}$ that are inside $\mathcal{X}_{ij}$; black
line: boundaries of the configuration space; dashed lines: boundaries of the
convex set $\mathcal{X}_{ij}$; red dotted lines: CBF constraints $h_{i}$.

The environment is implicitly represented by a simplified tree graph
$\mathcal{T}_{s}=(\mathcal{V}_{s},\mathcal{E}_{s})$ via sampling. According to
the simplified tree, the environment is decomposed the environment into a set
of convex cells.

Note that the landmarks $\hat{l}_{i}$, from the point of view of our
algorithms, can be arbitrary as long as there is at least two landmark visible
from any point $x$ in the free configuration space. The landmarks do not need
to be chosen from the samples of the RRT∗ algorithm, or from the obstacles.

###  II-E Control Lyapunov and Barrier Functions (CLF, CBF)

In this section, we review the CLF and CBF constraints, which are differential
inequalities that ensure stability and safety (set invariance) of a control
signal $u$ with respect to the dynamics (1). First, it is necessary to review
the following.

###### Definition 1

Given a sufficiently smooth function $h:n\to$ and a vector field $f:n\to n$,
we use the notation $L_{f}h=\nabla_{x}h^{\mathrm{T}}f$ to denote the Lie
derivative of $h$ along $f$, where $\nabla h$ represents the gradient field of
$h$.

The Lie derivative of a differentiable function $h(x)$ for the dynamics (1)
with respect to the vector field $B$ is defined as $\mathcal{L}_{B}h(x)$.
Applying this definition to (1) we obtain

$\mathcal{L}_{B}h(x)=\frac{\partial h(x)}{\partial x}Bu.$ (5)

In this work, we assume that Lie derivatives of $h(x)$ of the first order are
sufficient [24] (i.e., $h(x)$ has relative degree $1$ with respect to the
dynamics (1)); however, the result could be extended to the higher degrees, as
discussed in [10].

We now pass to the definition of the differential constraints. Consider a
continuously differentiable function $V(x):\mathcal{X}\to$, $V(x)\geq 0$ for
all $x\in\mathcal{X}$, with $V(x)=0$ for some $x\in\mathcal{X}$.

###### Definition 2

The function $V(x)$ is a Control Lyapunov Function (CLF) with respect to (1)
if there exists positive constants $c_{1},c_{2},c_{v}$ and control inputs
$u\in\mathcal{U}$ such that

$\displaystyle\mathcal{L}_{B}V(x)u+c_{v}V(x)\leq 0,\forall x\in\mathcal{X}.$
(6)

Furthermore, (6) implies that $\lim_{t\to\infty}V(x(t))=0$.

Consider a continuously differentiable function $h(x):\mathcal{X}\to$ which
defines a safe set $\mathcal{C}$ such that

$\displaystyle\mathcal{C}$ $\displaystyle=\\{x\in n|\;h(x)\geq 0\\},$ (7)
$\displaystyle\partial\mathcal{C}$ $\displaystyle=\\{x\in n|\;h(x)=0\\},$
$\displaystyle{Int}(\mathcal{C})$ $\displaystyle=\\{x\in n|\;h(x)>0\\}.$

In our setting, the set $\mathcal{C}$ will represent a convex local
approximation of the free configuration space (in the sense that
$x\in\mathcal{C}$ does not contain any sample that was found to be in
collision). We say that the set $\mathcal{C}$ is _forward invariant_ (also
said _positive invariant_ [25]) if $x(t_{0})\in\mathcal{C}$ implies
$x(t)\in\mathcal{C}$, for all $t\geq 0$ [26].

######  Definition 3 (CBF, [27])

The function $h(x)$ is a _Control Barrier Function_ with respect to (1) if
there exists a positive constant $c_{h}$, control inputs $u\in\mathcal{U}$,
and a set $\mathcal{C}$ such that

$\mathcal{L}_{B}h(x)u+c_{h}h(x)\geq 0,\forall x\in\mathcal{C}.$ (8)

Furthermore, (8) implies that the set $\mathcal{C}$ is forward invariant.

###  II-F Stability by CLF

To stabilize the navigation along an edges of a tree, we define the Lyapunov
function $V_{ij}(x)$ as

$V_{ij}(x)=z_{ij}^{T}(x-x_{j}),$ (9)

where $z_{ij}=\frac{x_{j}-x_{i}}{\lVert x_{j}-x_{i}\rVert}$ is the unit vector
of _exit direction_ for edge $(i,j)$, $x_{j}\in\text{\emph{exit face}}$ is the
position of the parent of node $i$, and $V_{ij}(x)$ reaches its minimum
$V(x)=0$ at $x_{j}$. Note that the Lyapunov function represents, up to a
constant, the distance $d(x,x_{j})$ between the current system position and
the exit face. By Definition 2, $V_{ij}(x)$ is a CLF.

###  II-G Safety by CBF

In this section, we define barrier functions $h_{ij}(x)$ that define a cone
representing a local convex approximation of the free space between $i$ and
$j$, in the sense that it excludes all samples in $\mathcal{V}$ that are on
the way from $i$ to $j$. In particular, we use the following steps:

  1. 1.

Define set $\mathcal{O}_{ij}\subset\mathcal{V}_{\text{collision}}$ whose are
in $\mathcal{X}_{ij}$.

  2. 2.

From the set $\mathcal{O}_{ij}$, we choose two samples $\\{o_{u},o_{d}\\}$ on
two sides of the edge $(i,j)$ such that they are closet samples in
$\mathcal{O}_{ij}$ to edge $(i,j)$.

  3. 3.

We write the equations of two lines passing through $\\{j,o_{u}\\}$ and
$\\{j,o_{d}\\}$ in a matrix form using $A_{h_{i}}\in 2\times 2$, $b_{h}\in 2$
to define the invariant set

$\mathcal{C}_{ij}=\\{x\in 2:A_{h_{ij}}\mathbf{x}+b_{h_{i}}>0\\}.$ (10)

The corresponding CBF is then defined as

$h_{ij}(x)=A_{h_{ij}}x+b_{h_{ij}}.$ (11)

An example of the set $\mathcal{C}_{ij}$ is shown in Fig. 2. Note that the
region $\mathcal{C}_{ij}$ might not include the entire cell
$\mathcal{X}_{ij}$. However, the controller will be designed to satisfy the
CBF and CLF constraints over the entire cell $\mathcal{X}_{ij}$ [11]).

###  II-H Output Feedback Controller

We assume that the robot can only measure the relative displacements between
the robot’s position $x$ and the landmarks in the environment, which
corresponds to the output function

$\mathcal{Y}=(L-x\mathbf{1}^{\mathrm{T}})^{\vee}=L^{\vee}-\mathcal{I}x=\operatorname{stack}{(l_{i}-x)},$
(12)

where $L\in n\times n_{l}$ is a matrix of landmark locations,
$i=1,\ldots,n_{l}$ that $n_{l}$ is the number of landmarks, $A^{\vee}$
represents the vectorized version of a matrix $A$,
$\mathcal{I}=\mathbf{1}_{nl}\otimes I_{n}$, and $\otimes$ is the Kronecker
product. We define the feedback controller as

$u_{ij}=K_{ij}\mathcal{Y},$ (13)

where $K_{ij}\in m\times nn_{l}$ are the feedback gains that need to be found
for each cell $\mathcal{X}_{ij}$.

In (13), the input to the controller is the relative displacement measurements
between the position of the robot and the set of landmarks in the environment.
The controller in (13) is a weighted linear combination of the measured
displacements $y$. The goal is to design $u$ such that the system is driven
toward the exit direction $z_{ij}$ while avoiding obstacles. Assume the
controller is in the form of (13), and the output function is (12), following
the approach of [10, 11], and using the CLF-CBF constraints reviewed in Sec.
II, we encode our goal in the following feasibility problem:

$\displaystyle\textrm{find}\;\;{K_{ij}}$ (14) subject to:
$\displaystyle\text{CBF:}\;-(\mathcal{L}_{B}h_{ij}(x)u+c^{\mathrm{T}}_{h}h_{ij}(x))\leq
0,$
$\displaystyle\text{CLF:}\;\;\mathcal{L}_{B}V_{ij}(x)u+c^{\mathrm{T}}_{v}{V_{ij}}(x)\leq
0,$ $\displaystyle u\in\mathcal{U},\;\;\forall
x\in\mathcal{X}_{ij},\;\;(i,j)\in\mathcal{E}.$

The constraints in (14) need to be satisfied for all $x$ in the region
$\mathcal{X}_{ij}$, i.e., the same control gains should satisfy the CLF-CBF
constraints at every point in the region. In [11], this problem is solved by
rewriting (14) using a min-max formulation and Linear Programming (LP) method.
Starting from a point $x\in\mathcal{X}_{ij}$, $u_{ij}$ drives the robot toward
$x_{j}$, the robot switches its controller to $u_{jq}$ when $x$ reaches the
exit face of $\mathcal{X}_{ij}$.

##  III FEEDBACK CONTROL PLANNING WITH BEARING MEASUREMENTS

As mentioned in the introduction, in this paper we consider a robot equipped
with a monocular camera that, in general, does not provide the depth of a
target object in the image, and instead measures the corresponding relative
bearing (Sec. II-B). In this section, we show that in this case, it is
possible to still use the control synthesis method of Sec. II-H after
rescaling the bearing measurements, such that they are similar to the ideal
displacement measurements. We divide the section into three parts. First, we
give details on the rescaling procedure; then, we show that the resulting
bearing controller still solves the path planning problem, albeit with
modified CLF and CBF conditions; finally, we discuss how to handle the
practical problem of a camera with a limited field of view.

###  III-A Rescaling of the Bearing Measurements

For simplicity, in this section we consider only 2-D environments (however, an
extension to the 3-D case is possible with minor modifications). Without loss
of generality, we assume that the robot position $x$ is equal to the zero
vector (if not, all the computations below are still valid after introducing
opportune shifts by $x$).

The bearings $\\{\beta_{i}\\}$ can be identified with points on a unitary
circle centered at the origin (Fig. 1a).

For the current cell (and hence its associated controller) we pick a _fixed_
landmark $f$ among those available. We then define $\tilde{l}_{f}$ as the
point on the circle corresponding to $\beta_{f}$. Our goal is to rescale all
the other bearings $i\neq f$ such that they are similar to the full
displacements.

For this purpose, we define _inter-landmark_ bearing directions between
landmark $f$ and every other landmarks $i\neq f$ as

$\beta_{i}^{f}=d_{i}(l_{f})^{-1}(l_{i}-l_{f});$ (15)

note that these bearings can be pre-computed from the known locations of the
landmarks in the map.

For each landmark $i\neq f$, let $P_{i}^{f}$ be the line passing through the
fixed landmark $\tilde{l}_{f}$ with direction $\beta_{i}^{f}$, and let $P_{i}$
be the line passing through the robot position $x$ with direction $\beta_{i}$.
We then define the scaled landmark position $\tilde{l}_{i}$ as the
intersection between $P_{i}^{f}$ and $P_{i}$.

###### Lemma 1

Let $s(x)=d_{f}(x)$ be the distance between landmark $f$ and the robot; then,
we have that $l_{i}-x=s(x)(\tilde{l}_{i}-x)$ for each landmark $i$.

###### Proof:

The triangles $l_{f},x,l_{i}$ and $\tilde{l}_{f},x,\tilde{l}_{i}$ are similar
since they have identical internal angles. Moreover,
$\lVert\tilde{l}_{f}-x\rVert=1$ by construction, the ratio between the
segments $l_{f},x$ and $\tilde{l}_{f},x$ is equal to $s(x)$. Combining these
two facts, we have that the ratio between the segments $l_{f},x$ and
$\tilde{l}_{f},x$ is also $s(x)$; the claim then follows. ∎

Our proposed solution is then to compute the scale displacements
$\tilde{\mathcal{Y}}=\operatorname{stack}(\\{\tilde{l}_{i}-x)$, which are then
used with the pre-computed controller

$\tilde{u}_{ij}(x)=K_{ij}\tilde{\mathcal{Y}}.$ (16)

###  III-B Analysis of the Bearing Controller

The following lemma shows that the original displacement-based controller
$u_{ij}$ and our proposed bearing-based controller $\tilde{u}_{ij}$ are
essentially equivalent from the point of view of path planning.

###### Proposition 1

Assume $s(x)$ is uniformly upper bounded (i.e., $s(x)<\infty$ for all
$x\in\mathcal{X}_{ij}$). The controllers $u_{ij}$ in (13) and $\tilde{u}_{ij}$
in (16) produce the same paths (but traced, in general, with different speeds)
for the drifless system (1) when started from the same initial condition.

###### Proof:

Let $x$ and $\tilde{x}$ be the trajectories of the system under $u_{ij}$ and
$\tilde{u}_{ij}$, respectively. Since both the dynamics and the controllers
are linear, we have that $\dot{x}=s\dot{\tilde{x}}$ when evaluated at the same
location. This implies that the two curves $x$ and $\tilde{x}$ are the same up
to a reparametrization of the velocity. ∎

In fact, we can also relate the new controller to the conditions in the
synthesis problem (14).

###### Proposition 2

Assume that $s_{\min}<s(x)\leq s_{\max}$, and that
$\tilde{u}=K_{ij}s\mathcal{Y}\in\mathcal{U}$ for all $x\in\mathcal{X}_{ij}$.
Then $K_{ij}$ is a feasible solution for (14) with the modified CLF and CBF
conditions:

$\displaystyle-(\mathcal{L}_{B}h_{ij}\tilde{u}+\tilde{c}_{h}h_{ij})$
$\displaystyle\leq 0,$ (17)
$\displaystyle\mathcal{L}_{B}V_{ij}\tilde{u}+\tilde{c}_{v}V_{ij}$
$\displaystyle\leq 0,$ (18)

where $\tilde{c}_{h}=s_{\min}^{-1}c_{h}$ and
$\tilde{c}_{v}=s_{\max}^{-1}c_{v}$.

###### Proof:

The claim follows by dividing the original CBF and CLF conditions by $s(x)$,
and then using the bounds $s_{\min}$, $s_{\max}$. ∎

Note that the bounds on $s(x)$ translate to bounds on the distance between the
cell $\mathcal{X}_{ij}$ and the landmarks $l_{i}$

Taken together, Propositions 1 and 2 show that the controller $K_{ij}$ found
by assuming a displacement-based controller can be used also for the bearing-
based case, although the speed of the resulting trajectories might be more
aggressive.

###  III-C Control With the Limited Field of View

In the formulation Sec. 14 it is implicitly assumed that the controller has
access to all the landmarks measurements at all times. However, in practice, a
robot will only be able to detect a subset of the landmarks due to a limited
field of view or environment occlusions (see Fig. 3b). To tackle this issue,
we show in this section that the controller $\tilde{u}$ in (16) can be
designed using all available landmarks (as in the preceding section), but then
implemented using two landmarks or more.

The main idea is to use the two available landmarks to estimate the missing
landmarks $\tilde{l}_{i}$ through the following steps (see also Fig. 3a):

  1. 1.

Choose one landmark $f$ as the fixed landmark, and denote as $f^{\prime}$ the
other landmark (in the figure $f=2$ and $f^{\prime}=1$). Compute
$\beta_{i}^{f}$ and $\beta_{i}^{f^{\prime}}$ for all
$i\notin\\{f,f^{\prime}\\}$.

  2. 2.

Following the same steps as in Sec. III-A, compute $\tilde{l}_{f}$ (from the
bearing $\beta_{f}$) and $\tilde{l}_{f}^{\prime}$ (via intersection).

  3. 3.

Let $P_{i}^{f}$ and $P_{i}^{f^{\prime}}$ be the lines passing through $l_{f}$
and $l_{f^{\prime}}$ in the directions $\beta_{i}^{f}$ and
$\beta_{i}^{f^{\prime}}$, respectively.

  4. 4.

Compute the non-visible modified landmark $\tilde{l}_{i}$ from the
intersection of $P_{i}^{f}$ and $P_{i}^{f^{\prime}}$.

After finding the modified bearings for all landmarks, the controller is
implemented as in Sec. III-A.

(a)

(b)

Figure 3: Implement controller with limited field of view. In Fig. 3a,
landmark $l_{3}$ is not visible from the robot camera. In Fig. 3b, the
modified observation $l_{3}$ is computed using the relative bearing
measurements between landmarks.

##  IV SIMULATION AND EXPERIMENTAL RESULTS

To assess the effectiveness of the proposed algorithm, we run a set of
validations using both MATLAB simulations and experiments using ROS on a
Jackal robot by Clearpath Robotics [28]. While the optimization problem
guarantees exponential convergence of the robot to the stabilization point, in
these experiments the velocity control input $u$ has been normalized to
achieve constant velocities along the edges of the tree.

(a) Configuration space

(b) Generated RRT*

(c) Simplified tree

(d) Simulation

Figure 4: Fig. 4a shows the configuration space of a kitchen-like environment,
where obstacles are represented by gray color. In Fig. 4b the root for RRT∗ is
located at the origin, red dots show the samples in collision with the
obstacles, and the tree generated by RRT∗ is plotted in green. Fig. 4c depicts
the corresponding simplified tree. In Fig. 4d, the agent is initialized at
different start points and moves toward the goal point.

###  IV-A MATLAB Simulation

The simulated MATLAB environment is presented in Fig. 4, where the obstacles
are represented by gray color. For RRT∗, we set the maximum number of
iterations in RRT∗ to $1500$, we choose $\eta=50$, and we place the root of
the tree at the origin of the environment (bottom left corner). The generated
tree from RRT∗ and its simplified form are shown in Fig. 4b and Fig. 4c
respectively. Then, we compute a controller for each edge of the simplified
tree as described in Sec. II-H. Fig. 4d shows the resulting trajectories from
four initial positions. In all cases, the robot reaches the desired goal
location by applying the sequence of controllers found during planning.

(a) Clearpath Robotics Jackal

(b) AprilTag

Figure 5: The Jackal robot used for the experiments is shown in Fig. 5a. We
use AprilTags (Fig. 5b) as the landmarks for the algorithm. Figure 6:
Experimental Environment

###  IV-B Clearpath Robotics Jackal Experiment

We used the Jackal robot [28] to test our algorithm in a lab environment. A
bird’s-eye view of the environment is presented in Fig. 6. We equipped the
Jackal with an Intel RealSense D455 camera [29]. The Jackal’s trajectory is
tracked by the OptiTrack motion capture system (with 44 infra-red cameras),
and we use AprilTags [30] with unique codes as the landmark fiducials, which
are placed in the environment at known positions and orientations with respect
to the inertial reference frame of OptiTrack. The data was managed by the
Robot Operating System (ROS) [31].

The experimental implementation of our algorithm requires taking into account
several practical considerations. First, due to the limited field of view
inherent in the RealSense camera, which prevents the Jackal from detecting all
required fiducials, we use the proposed method in Sec. III-C to compute the
controller based on two fiducials detected by the camera at each time instant.
Second, all bearing measurements described in previous sections are assumed to
be rotationally aligned with the global inertial reference frame. To
accommodate for this assumption, we compute the global rotation of the Jackal
with respect to the world reference frame as

${}^{W}\mathbf{R}_{J}={}^{W}\mathbf{R}_{AT}\left({}^{J}\mathbf{R}_{AT}\right)^{-1},$
(19)

where ${}^{J}\mathbf{R}_{AT}$ is the measured rotation of the AprilTag with
respect to the Jackal. Then, a normalized bearing ${}^{J}\mathbf{t}_{AT-J}$
detected by the Jackal can be aligned with the global reference frame as
${}^{W}\mathbf{t}_{AT-J}={}^{W}\mathbf{R}_{J}{}^{J}\mathbf{t}_{AT-J}$.

Finally, previous sections assumed a linear dynamical model for the robot,
while the Jackal has a unicycle dynamics. We map the original 2D input $u$ to
a linear velocity $u_{x}$ along the $x$ axis of the robot and an angular
velocity $\omega_{z}$ around the $z$ axis using a standard low-level
controller:

$\displaystyle
u_{x}=\small{\dfrac{\alpha}{\left\|u\right\|}\begin{bmatrix}\cos{\varphi}\\\
\sin{\varphi}\end{bmatrix}}^{\mathrm{T}}u,$
$\displaystyle\omega_{z}=\small{\frac{\beta}{\left\|u\right\|}\begin{bmatrix}0\\\
0\\\ 1\end{bmatrix}^{\mathrm{T}}\left(\begin{bmatrix}\cos{\varphi}\\\
\sin{\varphi}\\\ 0\end{bmatrix}\times\begin{bmatrix}u\\\
0\end{bmatrix}\right)}$ (20)

where $\varphi$ is the instantaneous yaw rotation of the robot with respect to
the world reference frame, which is extracted from ${}^{C}\mathbf{R}_{AT}$ and
${}^{W}\mathbf{R}_{AT}$, and $\times$ represent the 3-D cross product;
$\alpha$ and $\beta$ are user-defined scalar gains, 0.1 and 0.5 respectively.

We moved the Jackal across 4 different paths in the environment. The Jackal
followed the edges of the simplified RRT∗ tree and reached the expected goal
for all starting positions, despite the fact that the measurements were
obtained with vision alone, and despite different dynamics of the robot. The
video of the experiments is available at
https://github.com/Mahrooo/NavigationwithBearingMeasurements.

##  V CONCLUSIONS AND FUTURE WORKS

In this work, we proposed a novel approach to synthesize a set of out-put
feedback controllers on elements of a cell decomposition of the environment
that is based on the tree generated by a simplified version of sampling-based
RRT∗ method. We assumed the robot is equipped with a monocular camera and
there exists a global compass. We introduced a new approach to modify the
bearing measurements such that they form a uniform scaling of the set of
landmarks and defined the output feedback controllers take as input the
modified bearing measurements. Then, we computed the output-feedback
controllers by implementing the robust Linear Programming that satisfied the
stability of the system by the CLF constraint and satisfied the safety of the
system by CBF constraint. We validated our approach on both the simulation and
the experimental environment and represented the robustness of our approach to
the mismatches in the dynamical model of the robot, and to the measurements
acquired with a camera with a limited field of view. For the future, we plan
to prove to study the robustness of our algorithm to the deformation of the
environment. In addition, we plan to find the output feedback controller when
there exists an uncertainty in the bearing measurements.

## References

  * [1] J.-P. Laumond _et al._ , _Robot motion planning and control_. Springer, 1998, vol. 229.
  * [2] J.-C. Latombe, _Robot motion planning_. Springer Science & Business Media, 2012, vol. 124.
  * [3] H. M. Choset, S. Hutchinson, K. M. Lynch, G. Kantor, W. Burgard, L. E. Kavraki, and S. Thrun, _Principles of robot motion: theory, algorithms, and implementation_. MIT press, 2005.
  * [4] L. E. Kavraki, P. Svestka, J.-C. Latombe, and M. H. Overmars, “Probabilistic roadmaps for path planning in high-dimensional configuration spaces,” _IEEE transactions on Robotics and Automation_ , vol. 12, no. 4, pp. 566–580, 1996.
  * [5] S. M. LaValle, “Rapidly-exploring random trees: A new tool for path planning,” Iowa State University, Tech. Rep., 1998.
  * [6] S. M. LaValle and J. J. Kuffner Jr, “Randomized kinodynamic planning,” _The international journal of robotics research_ , vol. 20, no. 5, pp. 378–400, 2001.
  * [7] S. Karaman and E. Frazzoli, “Sampling-based algorithms for optimal motion planning,” _The international journal of robotics research_ , vol. 30, no. 7, pp. 846–894, 2011.
  * [8] S. M. LaValle, _Planning algorithms_. Cambridge university press, 2006.
  * [9] V. Renganathan, I. Shames, and T. H. Summers, “Towards integrated perception and motion planning with distributionally robust risk constraints,” _arXiv preprint arXiv:2002.02928_ , 2020.
  * [10] M. Bahreinian, E. Aasi, and R. Tron, “Robust planning and control for polygonal environments via linear programming,” _2020 IEEE American Control Conference (ACC)_ , 2020.
  * [11] M. Bahreinian, M. Mitjans, and R. Tron, “Robust sample-based output-feedback path planning,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, pp. 5780–5787.
  * [12] H. Zhao and Z. Wang, “Motion measurement using inertial sensors, ultrasonic sensors, and magnetometers with extended kalman filter for data fusion,” _IEEE Sensors Journal_ , vol. 12, no. 5, pp. 943–953, 2011.
  * [13] N. P. Papanikolopoulos, P. K. Khosla, and T. Kanade, “Visual tracking of a moving target by a camera mounted on a robot: A combination of control and vision,” _IEEE transactions on robotics and automation_ , vol. 9, no. 1, pp. 14–35, 1993.
  * [14] T. Yata, A. Ohya, and S. Yuta, “A fast and accurate sonar-ring sensor for a mobile robot,” in _Proceedings 1999 IEEE International Conference on Robotics and Automation (Cat. No. 99CH36288C)_ , vol. 1. IEEE, 1999, pp. 630–636.
  * [15] V. Lumelsky and K. Sun, “A unified methodology for motion planning with uncertainty for 2d and 3d two-link robot arm manipulators,” _The International journal of robotics research_ , vol. 9, no. 5, pp. 89–104, 1990\.
  * [16] D. Lambrinos, R. Möller, R. Pfeifer, and R. Wehner, “Landmark navigation without snapshots: the average landmark vector model,” in _Proc. Neurobiol. Conf. Göttingen_ , 1998.
  * [17] A. A. Argyros, K. E. Bekris, and S. C. Orphanoudakis, “Robot homing based on corner tracking in a sequence of panoramic images,” in _Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. CVPR 2001_ , vol. 2. IEEE, 2001, pp. II–II.
  * [18] M. Liu, C. Pradalier, Q. Chen, and R. Siegwart, “A bearing-only 2d/3d-homing method under a visual servoing framework,” in _2010 IEEE International Conference on Robotics and Automation_. IEEE, 2010, pp. 4062–4067.
  * [19] R. Tron and K. Daniilidis, “An optimization approach to bearing-only visual homing with applications to a 2-d unicycle model,” in _2014 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2014, pp. 4235–4242.
  * [20] A. Pentland, T. Darrell, M. Turk, and W. Huang, “A simple, real-time range camera,” in _1989 IEEE Computer Society Conference on Computer Vision and Pattern Recognition_. IEEE Computer Society, 1989, pp. 256–257.
  * [21] J. Ens and P. Lawrence, “An investigation of methods for determining depth from focus,” _IEEE Transactions on pattern analysis and machine intelligence_ , vol. 15, no. 2, pp. 97–108, 1993.
  * [22] J. Michels, A. Saxena, and A. Y. Ng, “High speed obstacle avoidance using monocular vision and reinforcement learning,” in _Proceedings of the 22nd international conference on Machine learning_ , 2005, pp. 593–600.
  * [23] H. M. Choset, K. M. Lynch, S. Hutchinson, G. Kantor, W. Burgard, L. Kavraki, and S. Thrun, _Principles of robot motion: theory, algorithms, and implementation_. MIT press, 2005.
  * [24] A. Isidori, _Nonlinear control systems_. Springer Science & Business Media, 1995.
  * [25] F. Borrelli, A. Bemporad, and M. Morari, _Predictive control for linear and hybrid systems_. Cambridge University Press, 2017.
  * [26] X. Xu, P. Tabuada, J. W. Grizzle, and A. D. Ames, “Robustness of control barrier functions for safety critical control,” _IFAC-PapersOnLine_ , vol. 48, no. 27, pp. 54–61, 2015.
  * [27] Q. Nguyen and K. Sreenath, “Exponential control barrier functions for enforcing high relative-degree safety-critical constraints,” in _2016 American Control Conference (ACC)_. IEEE, 2016, pp. 322–328.
  * [28] Clearpath, “Jackal unmanned ground vehicle,” https://clearpathrobotics.com/jackal-small-unmanned-ground-vehicle/, accessed: 2021-02-26.
  * [29] intel, “Intel realsense d455,” https://www.intelrealsense.com/depth-camera-d455/.
  * [30] J. Wang and E. Olson, “AprilTag 2: Efficient and robust fiducial detection,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, oct 2016, pp. 4193–4198.
  * [31] Stanford Artificial Intelligence Laboratory et al., “Robotic operating system, melodic morenia,” https://www.ros.org.