∎

11institutetext: Tahiya Salam 22institutetext: 22email<EMAIL_ADDRESS>33institutetext: M. Ani Hsieh 44institutetext: 44email<EMAIL_ADDRESS>
55institutetext: 1GRASP Laboratory
University of Pennsylvania, Philadelphia, USA

# Heterogeneous robot teams for modeling and prediction of multiscale
environmental processes

Tahiya Salam 1 M. Ani Hsieh 1

(Received: date / Accepted: date)

###### Abstract

This paper presents a framework to enable a team of heterogeneous mobile
robots to model and sense a multiscale system. We propose a coupled strategy,
where robots of one type collect high-fidelity measurements at a slow time
scale and robots of another type collect low-fidelity measurements at a fast
time scale, for the purpose of fusing measurements together. The multiscale
measurements are fused to create a model of a complex, nonlinear
spatiotemporal process. The model helps determine optimal sensing locations
and predict the evolution of the process. Key contributions are: i)
consolidation of multiple types of data into one cohesive model, ii) fast
determination of optimal sensing locations for mobile robots, and iii)
adaptation of models online for various monitoring scenarios. We illustrate
the proposed framework by modeling and predicting the evolution of an
artificial plasma cloud. We test our approach using physical marine robots
adaptively sampling a process in a water tank.

###### Keywords:

Environment Monitoring Heterogeneity Multi-Robot Systems Marine Robots

##  1 Introduction

Multi-robot systems have a long-standing history of being used to solve
problems that are distributed in time and space, with increased efficiency
(Parker, 2003). In the past, contributions to multi-robots teams have come
from swarm intelligence, optimal control and optimization, motion planning and
coordination, task allocation, distributed perception and estimation, decision
making, and learning. While many existing works consider teams composed of
identical robots, i.e., homogeneous teams, there is an increased interest in
the use of heterogeneous teams in recent years (Prorok et al, 2015, 2016;
Vasilijevic et al, 2015; Manjanna et al, 2018; Manderson et al, 2019; Maini et
al, 2018). Heterogeneity in multi-robot teams refers to the differing physical
or behavioral characteristics amongst agents. The impact of heterogeneity has
been studied in the context of cooperative sensing and coordination of multi-
agent robotic systems. Complex tasks such as search and rescue often benefit
or may require robots with different capabilities to coordinate and
collaborate to accomplish the task. Advantages of heterogeneous robot teams
over their homogeneous counterparts include the ability to introduce varying
degrees of cooperation amongst team members, sharing information from
different sources, faster exploration times, and increased robustness in the
event of failure. These advantages are most evident when deploying robot teams
to track dynamic and multiscale phenomena that varies in both space and time,
e.g., oil spills, forest fires, and air pollution. For example, more agile
aerial robots coupled with ground vehicles that can get closer to a target of
interest may provide better situational awareness in complex environments but
still maintain widespread surveillance of the area. While such heterogeneous
teams clearly bring significant advantages over homogeneous ones,
heterogeneity brings new challenges in coordination, decision-making,
planning, network connectivity, and information exchange. Furthermore, as
robots permeate all aspects of our lives, e.g., self-driving cars, package
delivery drones, autonomous environmental monitoring systems for wildlife
preservation, to name a few, multi-robot interactions will become more
ubiquitous. As such, the question of how best to coordinate these agents and
better leverage heterogeneity and diversity in capabilities for complex tasks
must be addressed.

In this work, we focus on the task of optimal sensing and online map merging
for a team of heterogeneous robots tracking a multiscale dynamic process. In
many environmental monitoring applications, the objective is to enable the
robot team to map and track dynamic processes such as fluid flows, animal and
insect swarms, forest fires, dispersion of airborne pollutants, and crowds in
complex environments (Dunbabin and Marques, 2012; Singh et al, 2010; Manderson
et al, 2019; Strogatz, 2001). These processes are often multiresolution in the
sense that they exhibit complicated patterns across various space-time scales.
The many disparate space-time scales makes it difficult to estimate the
process with a single, or even single type of, robot. A robot team with
heterogeneous mobility and sensing capabilities can be more adept at
collecting and fusing the multiresolution data needed to model and track these
multiscale dynamic processes. The robots can represent these processes as maps
and use the maps to adapt their sensing locations to improve them in an online
fashion. In the online model adaptation component, the previously derived
models of the dynamic process are updated with the acquisition of new data. In
multi-robot modeling, the environment is often assumed to be fixed. In these
static environments, the time scales tend be very long and changes in the
environment are mostly discrete, e.g. an office door is either open or closed.
However, this is not the case in modeling multiresolution dynamic processes,
as they are both varying in time and also exhibiting different time scales. In
these dynamic environments, the model must account for the added complexity of
parts of the model changing. Thus, an online model adaptation is essential for
quickly incorporating new process data to maintain model quality.

The core contribution of this work is addressing how to cohesively model an
unknown, complex, spatiotemporal environment using a team of heterogeneous
robots. We define this broad class of environments simply as dynamical systems
to allow us to use tools from dynamical systems theory to model the
environment. Such a comprehensive problem statement encompasses many scenarios
of interest for robotic applications. For example, this can include
incorporating various types of robots with diverse sensing and mobility
modalities to build a unified map of the environment, such as ground and
aerial vehicles jointly mapping temperature profiles or marine and aerial
vehicles jointly mapping ocean currents, as in Fig. 1. To that end, these
problems require tight coordination amongst robots to ensure the model is
accurate, consistent, and updated online.

At the same time, this work leverages the unique strengths of heterogeneous
multi-robot teams without explicitly formulating a task allocation problem
that would be computationally intractable (Korsah et al, 2013). This
alleviates the need to explicitly enumerate the requirements of the task and
match them to the traits of the robots. In the proposed strategy, the notion
of heterogeneity is encapsulated through fusion of multiscale data from
disparate sources into a cohesive, unified model. This work extends the
techniques for adaptive modeling and prediction through modal analysis
techniques to heterogeneous multi-robot systems.

Figure 1: Heterogeneous robots collecting different sensing information work
to create a cohesive model of a time varying environment. Aerial vehicles
collect low-fidelity sensor measurements, such as overhead images, over a wide
area, and marine vehicles collect high-fidelity sensor measurements, such as
current speeds, over a small area. Sensor measurements are unified into one
model for estimation and prediction of a time varying process

The proposed online mapping and task allocation strategy is rigorously
formulated, with computational complexity analysis provided. The strategy is
also accompanied by simulation studies using data of a multiresolution dynamic
processes. We present the development of the online task allocation and data
assimilation strategy. Experimental studies using mixed-reality teams,
consisting of simulated aerial vehicles and actual marine vehicles were
conducted using simulated data for a barium plasma cloud. This enables us to
successfully address and present experimental scenarios with truly multiscale
data of processes that exhibit complex, unpredictable behaviors.

To summarize, the technical novelty of this work lies in the introduction of a
novel, heterogeneous, multi-robot online modeling and estimation framework,
which determines optimal sensing locations from the estimated learned model of
a multiscale process representing the environment, in order to improve the
model accuracy. The remainder of the paper is organized as follows: Sect. 2
overviews the related works in heterogeneous multi-robot systems. Sect. 3
outlines the problem definition, the assumptions, and the mathematical
foundations for the algorithms used. The proposed approach and main
contributions are presented in Sect. 4. Computational analysis and simulation-
based evaluations are presented in Sect. 5. Experimental evaluation studies
are demonstrated in Sect. 6. Sect. 7 ends with conclusions and future research
directions.

##  2 Related works

In this section, we review the state-of-the-art in heterogeneous multi-robot
systems, with special attention to current solutions for multi-robot modeling
and data fusion with heterogeneous data.

###  2.1 Heterogeneous multi-robot systems

In the multi-robot literature, heterogeneity is often framed as a task
allocation problem (Jones et al, 2006; Rossi et al, 2009; Korsah et al, 2013;
Prorok et al, 2015, 2016; Bae et al, 2019; Liang et al, 2019; Notomista et al,
2019). Starting with the classification of multi-robot coordination by Gerkey
and Matarić (2004), followed by its categorization by Korsah et al (2013), and
reviewed by Khamis et al (2015), existing work has mainly focused on
coordination in the form of task-to-trait allocation. In this framework,
robots possess traits that allow them to complete task specifications. These
traits can include mobility capabilities, on-board sensors, and/or other
characteristics that allow the robots to make progress on or complete
specified tasks. Though the premise of traits has been used in prior works,
the terms were first formalized by Prorok et al (2015) for heterogeneous
teams. According to Prorok et al (2015), robots possess a set of traits that
allow them to accomplish tasks where these tasks can also be decomposed in
terms of traits. Such a discretization of both the robots’ capabilities and
the tasks’ requirements may not always be feasible, especially in situations
where the tasks cannot be easily decomposed and/or robots capabilities may not
exactly fulfill the tasks’ requirements. Even in situations where one
considers the distributions of traits rather than their binary existence, such
a formulation may ignore some of the original requirements of the task (Prorok
et al, 2015, 2016). For example, in this work the task is to reconstruct and
predict the evolution of the environment using heterogeneous robots. The
discretization of this task cannot be easily decomposed into the sensing
capabilities of the individual robots, as all of the collective information
must be leveraged to accomplish the task.

###  2.2 Multi-robot sensing and modeling

The question of modeling with multi-robot teams has been a topic of interest
in the robotics community for many years (Nashashibi et al, 1992).
Heterogeneity in sensing and modeling has been explored in the coverage
control literature. Coverage control encompasses a type of optimal sensor
allocation problem. Coverage control has a long-standing history of using
multi-robot teams to sense environmental features through a density function
(Cortés et al, 2004). The density function used in coverage control techniques
can be thought of as a model of the process of interest. For heterogeneous
sensing modalities, each sensing modality is represented as a different
density function and the optimal sensor allocation is solved as an
optimization problem over these combined density functions (Santos et al,
2018; Santos and Egerstedt, 2018; Sadeghi and Smith, 2019). In the cases where
the model is unknown a priori, robots perform a distributed estimation of
unknown density functions (Julian et al, 2012). However, these techniques may
not be suited for fast or unknown dynamics where the density function changes
rapidly in time or is unknown to the robots. Lastly, these techniques neglect
how to assimilate new sensor measurements from robots.

Gaussian Process (GP) modeling has been used in multi-robot systems for
modeling, determining optimal sensor placement, prediction, and path planning
(Krause and Guestrin, 2007; Krause et al, 2008; Singh et al, 2009, 2010; Xu et
al, 2011). GPs are a non-parameteric model that can be derived using data
(Rasmussen, 2004), which makes them extremely appealing for modeling the
environment using data collected by robots. These techniques have been
extended and applied to heterogeneous teams, where heterogeneity refers to
robots of different types (e.g. aerial and ground) or different sensing
modalities (e.g. water quality and sampling apparatus) (Manderson et al, 2019;
Manjanna et al, 2018). In these scenarios, robots are delegated specialized
roles that require fulfilling disparate tasks in parallel, such as a robot
building a GP model while the other samples. Requiring one of the robots to
build the GP model means the robot runs into the same problems associated with
standard robotic GP modeling applications. For one, determining the
appropriate kernel function that actually captures the complex spatiotemporal
relationships in the environment is a challenging task and still remains an
open question (Singh et al, 2010). Additionally, the size of the covariance
matrix increases as more measurements are collected making its inverse
expensive to compute. Some studies propose efficient methods for using only a
subset of observations instead of all observations to compute the covariance
matrix (Xu et al, 2011). However, these methods impose restrictions on the
covariance functions.

Modal analysis techniques have been used to study spatiotemporal systems and
extract the important physical features (Taira et al, 2017). Modal analysis
techniques for fluid flows, finance, video processing, and robotics have been
well studied (Mezić, 2005; Rowley et al, 2009; Schmid, 2010; Tu et al, 2013;
Jovanović et al, 2014; Brunton et al, 2016; Zhang et al, 2019; Matsumoto and
Indinger, 2017; Manohar et al, 2018, 2019; Salam and Hsieh, 2019). Dynamic
mode decomposition (DMD) is a modal analysis technique that decomposes time
series data into modes, where each mode has a characteristic frequency and
growth or decay rate (Schmid, 2010). DMD relies only on time series data and
does not require governing equations. This method combines advantages from two
powerful data analytic tools: power spectral analysis for temporal analysis
and principal components analysis (PCA) for spatial analysis. Additionally,
DMD can be used for prediction of future states, optimal sensor placement,
estimation, and control, which can be useful in many robotic applications
(Taira et al, 2017; Manohar et al, 2018). Despite DMD and similar techniques
being used for robotic applications, most uses of DMD in robotics have focused
on control and not modeling (Berger et al, 2015; Folkestad et al, 2020).

###  2.3 Data fusion for heterogeneous multi-robot teams

Map merging techniques are methods for fusing information into a global map,
when multiple robots have been used to explore the environment and collect
information. Map merging focuses mainly on map matching, the problem of
finding the correspondences between maps, and map fusion, the problem of
merging the data from maps in the case where there is a known or estimated
correspondence between the maps. For a review of the state of art in map
merging techniques, we refer the interested reader to the work by Andersone
(2019). Most map merging techniques focus on homogeneous robots. Works
focusing on heterogeneous robots have mostly focused on map fusion for
occupancy grid maps of various scales (Topal et al, 2010; Park et al, 2016; Ma
et al, 2016; Shahbandi and Magnusson, 2019). However, there is no universal
solution to map merging, especially for heterogeneous robots, due to the
varying degrees of heterogeneity in format, sensing type, and scale of data.
As such, there are still many open challenges to merging of data and
developing coordination strategies for optimal sensing and modeling by
heterogeneous robot teams.

###  2.4 Contributions

The work presented in this paper avoids the formulation of heterogeneity as a
task allocation problem and instead considers multiple tasks that must be
solved sequentially using all available traits of the team. Unlike coverage
control and GP techniques, the proposed approach explicitly considers the
impact of heterogeneous data sources and online data assimilation in building
a single map. Likewise, map merging is still not well studied for robotic
frameworks leveraging heterogeneity. As mentioned previously, the core
contribution focuses on modeling and optimal sensing tasks and requires the
coupling of different resources from the heterogeneous robots, while
addressing the challenges of map merging and online updates.

##  3 Problem statement

The main contribution of this work is the synthesis of techniques from
nonlinear dynamical systems theory and robotics in the development of a
distributed estimation and sensing framework for environmental monitoring by
heterogeneous robot teams. The proposed framework enables us to model the
environment in principled ways while accounting for the constraints and
challenges posed by multi-robot systems. In the proposed framework,
heterogeneous robot capabilities are abstracted through the reduced-order
models created by the robots. For example, aerial vehicles generally have
wider fields of view, while marine vehicles are able to directly sample the
process. Thus, the heterogeneity of the robots is reflected through their
collection of multiresolution data. Aerial vehicles collect data with low
spatial resolution at higher frequency and marine vehicles collect data with
high spatial resolution at slower frequency.

##### Problem statement:

This work extends the techniques for adaptive modeling and prediction through
dimensionality reduction and learning techniques to heterogeneous multi-robot
systems. Accomplishing this requires combining data from different sources and
of varying resolutions into one cohesive model, while maintaining the
independence of data sources. In order to resolve the various temporal rate
and spatial resolutions, the proposed scheme ensures that lower resolution
data collected at faster rates are interpolated to the same resolution as the
higher resolution data collected at lower rates. The multi-resolution data is
aggregated in a centralized fashion. Model-learning techniques are then
employed on this resolved data to obtain a map capable of capturing the
spatiotemporal dynamics of the system, which is then used to guide the
allocation of the aerial and marine robots. The overall problem is comprised
of three main aspects: 1) model estimation and uncertainty quantification, 2)
uncertainty-reducing task assignment, and 3) online adaptation of models,
described in Sections 3.2, 3.3, and 3.4 respectively.

###  3.1 Background and assumptions

The proposed framework employs DMD for model estimation and uncertainty
quantification. For a more detailed exposition on DMD and its relationship to
dynamical systems theory for nonlinear differential equations, we refer the
interested reader to (Mezić, 2005; Rowley et al, 2009; Budišić et al, 2012).

Given $T+1$ sequential snapshots
$\\{\mathbf{x}(0),\mathbf{x}(1),\dots,\mathbf{x}(T)\\}$, let each
$\mathbf{x}(t)\in\mathbb{R}^{N}$ represents field values at $N$ spatial
locations at time $t$. Then, assume there is a linear dynamical relationship

$\mathbf{x}(t+1)=\mathbf{Ax}(t)$ (1)

between two consecutive snapshots. The matrix $\mathbf{A}$ is typically
extremely high-dimensional. DMD allows us to learn a reduced-order model of
the matrix $\mathbf{A}$ that captures the important spatial and temporal
characteristics of the data. Let

$\displaystyle\mathbf{X}=[\mathbf{x}(0)\quad\mathbf{x}(1)\quad\dots\quad\mathbf{x}(T-1)]$
(2)
$\displaystyle\mathbf{Y}=[\mathbf{x}(1)\quad\mathbf{x}(2)\quad\dots\quad\mathbf{x}(T)]$
(3)

be defined as the snapshot matrices, a matrix with its columns as snapshots of
collected data. Then, in matrix notation, Eq. (1) can then be written as
$\mathbf{Y}=\mathbf{A}\mathbf{X}$. Thus, DMD is an approximate
eigendecomposition of the operator

$\mathbf{A}=\mathbf{Y}\mathbf{X}^{\dagger},$ (4)

where † is the pseudoinverse operator.

While this is a linear approximation, many of the papers cited above mention
the idea that DMD is able to characterize nonlinear dynamics through an
analysis of some approximating linear system (Tu et al, 2013). Rowley et al
(2009) established the connection between DMD and nonlinear dynamical systems
and helped explain the validity of DMD when applied to nonlinear systems.

Given that the DMD operator $\mathbf{A}$ may be high dimensional and difficult
to compute, we use an efficient algorithm that computes the eigendecomposition
of $\mathbf{A}$ through a low-dimensional approximation (Schmid, 2010). For
this procedure, we compute a singular value decomposition (SVD) on
$\mathbf{X}$ such that

$\mathbf{X}=\mathbf{U}\bm{\Sigma}\mathbf{W}^{T}$ (5)

where $\bm{\Sigma}$ is an $r\times r$ diagonal matrix of non-zero singular
values $\\{\sigma_{0},\sigma_{1},\dots,\sigma_{r-1}\\}$ and $r$ is the rank of
the snapshot matrix. Then, combining Eqs. (4) and (5) gives
$\mathbf{Y}=\mathbf{A}\mathbf{U}\bm{\Sigma}\mathbf{W}^{T}$ and therefore an
approximation of $\mathbf{A}$ can be given by

$\displaystyle\hat{\mathbf{A}}=\mathbf{U}^{T}\mathbf{A}\mathbf{U}$ (6)
$\displaystyle\hat{\mathbf{A}}=\mathbf{U}^{T}\mathbf{Y}\left(\mathbf{X}\right)^{\dagger}\mathbf{U}=\mathbf{U}^{T}\mathbf{Y}\mathbf{W}\bm{\Sigma}^{-1}\mathbf{U}.$
(7)

Given the operator $\hat{\mathbf{A}}$, we can compute its eigendecomposition
as

$\hat{\mathbf{A}}\mathbf{V}=\mathbf{V}\mathbf{\Lambda}$ (8)

where $\mathbf{V}$ contains the eigenvectors of $\hat{\mathbf{A}}$ and
$\mathbf{\Lambda}$ is a diagonal matrix of the eigenvalues, $\lambda_{i}$, of
$\hat{\mathbf{A}}$.

The DMD modes can then be computed as

$\mathbf{\Phi}=\mathbf{Y}\mathbf{W}\bm{\Sigma}^{-1}\mathbf{V}$ (9)

where each column,
$\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\phi}_{i}$ of
$\mathbf{\Phi}$ is the DMD mode corresponding to the eigenvalue, $\lambda_{i}$
found in Eq. (8). The reconstruction of the data is then written as

$\hat{\mathbf{x}}(t)=\mathbf{\Phi}\mathbf{\Lambda}^{t}\bm{\alpha}$ (10)

where $\bm{\alpha}$ is computed using least-squares for
$\mathbf{x}(0)=\mathbf{\Phi}\bm{\alpha}$. Thus, the spatial modes are captured
by the vectors
$\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\phi}_{i}$ and
the temporal dynamics of the spatial modes are captured by $\lambda_{i}$. The
magnitude of the vector
$\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\phi}_{i}$
represents the spatial correlations between the set of locations, and the
magnitude and phase components of the eigenvalue $\lambda_{i}$ represents the
rate of growth/decay and frequency of oscillation of the corresponding spatial
mode
$\textpdfrender{TextRenderingMode=FillStroke,LineWidth=.35pt,}{\phi}_{i}$.

Given a workspace, $\mathcal{W}\subset\mathbb{R}^{d}$ for $d\in\\{2,3\\}$
discretized into $N$ points, each point can be denoted as
$p_{i}\in\mathbb{R}^{d}$, for $i\in I=\\{1,\dots,N\\}$, where $I$ is the index
set of the discretization. That is, each point $p_{i}$ can be uniquely
identified by an index $i$. The time varying process of interest $\mathcal{P}$
is observed over the workspace $\mathcal{W}$ and sensed over the times
$\\{0,\dots,\mathcal{T}\\}$. Each point $p_{i}$ has an associated weight,
$w_{i}$, corresponding to the informativeness of that location with respect to
mapping $\mathcal{P}$ and some specified statistical criterion. Robots are
able to collect various sensor measurements over a set of points $S_{i}$,
where $S_{i}\subseteq{I}$, depending on their sensing quality, sensing radius,
and location $p_{i}$ in the workspace. Let $m$ be the number of robots on a
team with one type of sensing modality, $m^{\prime}$ be the number of robots
with a different type of sensing modality, and the total number of robots be
$M=m+m^{\prime}$. Each robot $j$, for $j\in\\{1,\dots,M\\}$, collects
information at different time scales $T_{j}$, where
$T_{j}\subseteq\\{0,\dots,\mathcal{T}\\}$, depending on the robot’s sensing
capabilities. This means that each robot $j$ at location $p_{j}$ has its
unique sensing locations $S_{j}$ and collects data at the times $T_{j}$. In
this case, the sensing locations are non-overlapping by the allocation of the
robots but the data collection times may be overlapping.

###  3.2 Model estimation and error quantification

In order for the heterogeneous robots to build a model of the previously
unknown area, the proposed strategy needs to provide a framework for a
compact, meaningful representation of the environment. Hence, the
heterogeneous data collected from the robots should be unified to reconcile
the differences in spatial and temporal scales, build a cohesive model, and
extract the key features of the environment from the model. The multiscale,
multiresolution data collected should be combined as input to the algorithm
that constructs a model.

The environment is discretized into points within the workspace and
represented as a state $\mathbf{x}$, where the state changes in time. This
representation of the environment as a states lends to high-dimensional
systems. Given that the high-dimensionality of these systems is
computationally prohibitive, the model derived by the robots should have lower
dimensionality and elucidate the key spatial and temporal characteristics of
the process to allow for planning. Selecting a model that highlights the
meaningful features of the various spatial and temporal properties of the
dynamical system, such as DMD, allows for not only model reconstruction but
also model-informed planning of optimal sensing locations. In order to so, the
selected methodology for modeling has to support a measure of informativeness
with respect to sensing information, wherein collecting data in different
locations has measurably different effects on model quality. Additionally, the
model should reduce the error in reconstruction from limited sensing data
collected from mobile robots.

###  3.3 Uncertainty-reducing task assignment

The optimal sensor allocation problem is formulated as an optimization over
the developed map, specifically the spatial and temporal components captured
in $\mathbf{\Phi}$ and $\mathbf{\Lambda}$, as seen in Eq. (7) and (9). During
this phase, the robots are assigned to the optimal sensor locations that
lowers the map uncertainty. This couples the coordination of the robot team, a
task allocation problem, to the quality of completing the task, increasing
model quality.

Planning where mobile robots should go next to improve the model is essential
in the online planning framework. Given that the robots are operating in an
extremely high dimensional space and only have limited sensing radius and
capacity, this entails finding a small subset of optimal sensing locations
such that the uncertainty is reduced or the information gain is maximized.
Given informativeness of optimal sensing locations, determination of where
robots should go to collect the sensing measurements is NP-hard. However,
instead of task to trait assignment, as in the heterogeneity literature
(Korsah et al, 2013), this task allocation problem is restrained to the domain
of assigning sensing regions to robots, as in the spirit of coverage problems.

To begin, the model is unknown and the robots collect sensor measurements to
build an initial model of the environment. After an initial model is built,
the aim is to determine which locations of radius $k$ are expected to lead to
the most uncertainty reduction in estimating the model. These locations are
reevaluated and redetermined after more sensing data is collected and the
models are updated.

For every point $p_{i}\coloneqq(x,y)\in\mathcal{W}$, we can define a sensing
region, $S_{i}$, from that location $p_{i}$ as a collection of points,
$\\{p_{j}\coloneqq(x^{\prime},y^{\prime})\\}$, such that

$S_{i}\coloneqq\\{p_{j}\enskip|\enskip||(x-x^{\prime},y-y^{\prime})^{\top}||\leq
k\\}.$ (11)

The cumulative weight $W_{i}$ of each sensing regions $p_{i}$ is then
$W_{i}=\sum\limits_{p_{j}\in S_{i}}w_{j}$.

Define $L$ as a set containing the elements within the sets
$S_{1},S_{2},\dots,S_{m}$ that are arbitrary sensing regions, such that for
any two sets $S_{i}$ and $S_{j}$ with elements in $L$, $S_{i}\cap
S_{j}=\varnothing$. Thus, $L$ corresponds to all of the locations in $m$ non-
overlapping sensing regions. Note, there are multiple ways to satisfy the
definition of $L$. Intuitively, this means $L$ contains all of the points for
which sensing data is being collected over a specific allocation of the robots
to certain locations. This is shown in detail in Fig. 2.

(a)

(b)

(c)

Figure 2: Depiction of selecting $15$ candidate sensing regions from
workspace. a) The workspace $\mathcal{W}$ is the space over which the time
varying process $\mathcal{P}$ is observed. b) The process $\mathcal{P}$ can be
divided into candidate sensing regions, $S_{i}$, where for each set represents
a set of points that a robot would be able to sense if they were located in
the center of the region. c) In order to construct a valid set $L$ that
contains all of the points from one possible assignment of robots to regions,
overlapping sets $S_{5}$, $S_{8}$, $S_{13}$, $S_{18}$, and $S_{19}$ have been
removed. Note, that there are many possible sets $L$, such as the set of
points resulting from removing $S_{10}$ instead of $S_{8}$. The whitespace
$\mathcal{W}\setminus L$ is inferred using a model and the measurements in the
set $L$

The values of the process $\mathcal{P}$ can be inferred at all locations
without available sensor measurement, defined as the points in the set
$\mathcal{W}\,\setminus\,L$, for a specific $L$ given the choice of model and
an appropriate estimation scheme.

The optimization problem is then

$L^{*}\coloneqq\operatorname*{argmax}\limits_{L}\sum\limits_{S_{j}\subset
L}W_{j},$ (12)

as in finding the set of optimal locations $L^{*}$ such that the cumulative
weight over these $m$ sensing locations is maximized.

This is a combinatorial optimization problem for $m$ instances over the entire
workspace of dimension $N$, where $N$ is a very large number, with a
computational complexity of $\Theta(N^{m})$. For the systems that are being
modeled and studied, explicitly solving this combinatorial optimization
problem is infeasible. Thus, we need to devise a computationally efficient
method for determining the $m$ best sensing regions that approximates the
combinatorial optimization solution.

###  3.4 Online adaptation of models

In this work, the data can be assimilated at varying frequencies, allowing for
single updates, with the acquisition of data at a single time, or batch
updates, with the acquisition of data over some period of time. This
flexibility reduces connectivity requirements. For example, robots can explore
their assigned optimal sensing locations for some time without maintaining
connectivity with the rest of the robots, to return after some period of time,
share their acquired data, and update the model according to the online model
adaptation algorithm.

The assimilation of newly acquired data into the models after a sensing period
is a fundamental component of modeling a dynamic environment. In trying to
model a dynamic process, new sensing data must be assimilated into the
existing model to adequately capture the changing environment. Given the
dimensionality of the data and the distributed nature of collecting sensing
measurements, it is infeasible for all robots to keep the full time series
data and global models on-board. The techniques for online assimilation of the
data to the existing model may also depend on the duration of the process and
computational limitations of the vehicles. The key challenge in this framework
is determining which data and features are needed in updating and maintaining
models of the environment.

Different assimilation strategies may be used for different types of dynamics
and environments. Some environmental processes require long-term monitoring of
extremely high-dimensional spaces. Alternatively, other types of environmental
processes are quickly time varying processes with lower dimensionality. Thus,
various methodologies with diverse computational and algorithmic constraints
are required to efficiently adapt the models. This is addressed by selecting
different model adaptation strategies for the type of monitoring. For example,
the technique for extremely high-dimensional spaces has to address memory and
computational efficiency. Alternatively, the technique for quickly time
varying processes has to provide greater emphasis on newly acquired data.
Furthermore, the adaptation techniques are constrained by the choice of model,
as in the online data assimilation mechanisms must be consistent with the
modeling and optimal sensing strategies devised.

##  4 Proposed approach

Figure 3: Schematic of proposed framework with corresponding section numbers.
Noisy low fidelity measurements from aerial vehicles are upsampled and sparse
high-fidelity sensor measurements from marine vehicles are used to extrapolate
the full state. Estimates are combined and used as input to the modal
decomposition algorithm that provides decoupled temporal and spatial
characteristics. Characteristics are used independently to determine the
optimal sensing locations for data collection and eventual data assimilation

The proposed framework, demonstrated in Fig. 3, solves the model estimation
and uncertainty quantification through interpolation and extrapolation
techniques that allow for the combination of various data sources and outputs
temporal and spatial characteristics of the environment. The uncertainty-
reducing task assignment is solved efficiently by leveraging these
characteristics from the modal decomposition. Finally, the solution for online
adaptation of the model relies on the properties of the modes itself.

###  4.1 Reconciling spatial and temporal scales for model construction

Figure 4: Overview of data collected by heterogeneous multi-robot team used to
produce low- and high-fidelity estimates of field. Low-fidelity measurements
from the aerial vehicles are used to create low-fidelity estimates of the true
process as shown at $\text{time}=4$. High-fidelity sensor measurements, which
are not always available due to sensing and computational limits, are
collected at specific sensing locations by the marine vehicles, as shown by
the circular regions at $\text{time}=1$. A more detailed combined model is
constructed using aerial vehicle model estimates and sparse sensor
measurements from the marine vehicle. This model can be used for prediction
when there is no marine vehicle data available. Predictions and new sensor
measurements are assimilated into the model. The combined model is able to
capture more interesting patterns through the strategic placement of the
sensors, as seen at $\text{time}=5$

In order to determine the informativeness of the points in the workspace,
first we must define a method for unifying the measurements from disparate
sources. This requires combining data into one cohesive model, while
maintaining the independence of data sources. Though the data from different
sources is attempting to model and predict the same phenomena, they come at
various temporal and spatial resolutions. In order to resolve the various
temporal rates and spatial resolutions, we assure that the data collected at
faster rates but with lower spatial resolutions are interpolated to the same
resolution as the data collected at higher spatial resolution. This higher
spatial resolution data is being collected at a lower rate. Additionally, data
collected at higher spatial resolution but with sparse measurement, as in the
case where robots can densely sample but only at limited sensing location, is
estimated at high spatial resolution across the entire discretization of the
workspace. This is illustrated in Fig. 4.

(a)

(b)

Figure 5: Sampling of low-fidelity samples to higher spatial resolution. The
original field has measurements on a $96\times 384$ spatial grid. (a) At time
= $1.004$ secs, low-fidelity measurements on a $10\times 39$ grid are
upsampled to same spatial scale as high-fidelity measurement on a $48\times
192$ grid and closely agree with the true field. (b) At time = $97.22$ secs,
upsampling of low-fidelity measurements to same spatial scale as high-fidelity
measurement fail to capture relevant features of true field

In the model-dependent optimal sensor allocation component, the data collected
by robots is combined and then represented using DMD (Schmid, 2010). For
completeness, we review some of the key ideas behind DMD, first presented in
Sect. 3.1, again here. DMD learns a low-dimensional model that contains a set
of modes, where each mode is associated with a fixed oscillation frequency and
decay or growth rate, from the collected data. Using DMD allows for the
observation of oscillations in dynamic processes with both short and long time
scales (Tu et al, 2013). Additionally, the output model from DMD has
physically interpretable meaning and can be used to determine optimal sensing
locations in space (Mezić, 2005; Manohar et al, 2018, 2019). However, given
the inherent multiresolution nature of the data collected by the heterogeneous
robots, the data must first be resolved into a cohesive model. There are
various ways to do this, including sequential DMD, non-sequential DMD, (Tu et
al, 2013). In this paper, we use standard interpolation techniques to resolve
the data to a uniformly sampled model.

In Sects. 5 and 6, the focus is on marine vehicles sampling with high spatial
resolution at sparse locations, and aerial vehicles sampling at low spatial
resolution but with a faster temporal rate. While the proposed approach holds
true for this scenario, we argue that this framework is also applicable to
heterogeneous teams in general. The approach provides a model and prediction
scheme irrespective of the natures of the data sources provided.

We propose combining the measurements at varying temporal and spatial scales
of the process into one compact model that provides the dominant spatial and
temporal characteristics of the process. Initially, models of two types are
constructed: high-fidelity models with sparse measurements from marine robots,
as they collect precise measurements of the environment but only within
sensing regions, and low-fidelity models with abundant measurements from
aerial robots, as they survey large areas but cannot obtain precise samples.
High-fidelity models provide information for areas of the dynamic process to
allow for locally optimal estimates of the environment. Low-fidelity models
characterize high-level descriptions of the dynamic process. While we use
aerial and marine robots, more generally, multi-fidelity models encapsulate
heterogeneity from various sources, such as from robots with different
mobility or sensor types. Models from aerial robots are processed at a much
higher frequency, due to their increased computation power, while models from
marine robots are processed at lower frequency due to the fact that the data
comes in much slower and must be processed to estimate the full field. The
aerial vehicles maintain a coverage formation and upsample their low-fidelity
measurements to the same spatial resolution as the high-fidelity measurement
(Cortés et al, 2004).

Initially, for some time $T+1$, the low-fidelity measurements from the aerial
vehicles are collected. The aerial vehicles assume their positions using
coverage control and assume a uniform density function across the space. The
low-fidelity measurements are upsampled using bilinear interpolation, such
that spatial resolution of the interpolated data is the same as the candidate
high-fidelity data. In this scenario, standard interpolation techniques are
easy to implement, computationally efficient, and perform just as well as
other complicated learning techniques, such as neural networks. In our
studies, we use standard interpolation as there was not a statistical
difference in estimation when compared to other learning techniques and
standard interpolation does not require prior training and historical
knowledge. However, this method is general enough to use alternative
estimation techniques, so long as the dimensions are compatible with the high-
fidelity estimates. While there is obviously some relationship between the
various fidelity sensor measurements, we do not take this into account in the
upsampling as to account for the independence in sensor measurements from
disparate sensing sources. With this technique, we capture the coarse spatial
features of the environment but at a much higher frequency. Unsurprisingly,
for simple spatial features, this approach is successful, as shown in Fig. 5
LABEL:sub@interp-early. In these instances, the interpolated values closely
match the true values of the field. However, for more complex features, the
simple interpolation method is alone not sufficient Fig. 5 LABEL:sub@interp-
late. The interpolated do not match the true values of the field.

Given a dataset that is now characterized by a fast temporal rate and high
spatial resolution, we can use dimensionality reduction techniques to i) build
a cohesive model of the system that can be used for estimation of missing
sensor values, and ii) determining the relative importance of points in the
workspace to best estimate the dynamic process. Specifically, given a combined
data set $\mathbf{X}$, we can use the DMD analysis described in Sect. 3 to
compute a reduced order model. As given by the reconstruction equation, Eq.
(10), any state of the dynamical system at time $t$, $\mathbf{x}(t)$, can be
approximated using the DMD modes and eigenvalues as
$\hat{\mathbf{x}}(t)=\mathbf{\Phi}\mathbf{\Lambda}^{\top}\bm{\alpha}$. Note,
that for the DMD modes, $\mathbf{\Phi}\in\mathbb{R}^{N\times r}$, where $N$ is
the dimension of the discretization of the workspace. This means that the row
$i$ of the matrices $\mathbf{\Phi}$ and $\mathbf{\Lambda}$ corresponds to a
specific point $p_{i}$ from the discretization of the workspace. In Sect. 3,
we discussed how $\mathbf{\Phi}$ captures the key spatial modes and
$\mathbf{\Lambda}$ captures the temporal dynamics of the spatial modes. We
will use the spatial information in the approximation relationship, captured
by the matrix $\mathbf{\Phi}$, to reconstruct the system using the high-
fidelity sensor measurements from the marine robots. Similarly, this
relationship will also be leveraged in determining the optimal sensing
locations for the marine robots. This will be discussed in further detail in
the following section. Alternatively, we will use the temporal information
information, captured by the matrix $\mathbf{\Lambda}$ to construct a density
function for a coverage formation for the aerial vehicles (Cortés et al,
2004). As the components of $\mathbf{\Lambda}$ represent either the rate of
growth or decay and frequency of oscillation of the corresponding spatial
modes, we use this to emphasize areas of either high growth or high decay in
the density function in an attempt to capture some of the quickly time varying
features that may be missed by the marine robots.

###  4.2 Fast task allocation for optimal sensing locations

Given the spatial and temporal characteristics of the environment, we can use
these to inform us of the next sensing locations for the robots. For robots
sampling at select sensing regions and receiving high-fidelity data, the
optimization problem in Eq. (12) can be solved using the spatial
characteristics (Sect. 4.2.1). To determine optimal sensing locations for
robots receiving low-fidelity measurements over larger regions at a faster
rate we rely on standard coverage control techniques using the temporal
characteristics (Sect. 4.2.2).

####  4.2.1 Block QR pivoting using spatial characteristics

To begin, $m$ robots collecting high-fidelity sensing measurements are
initially deployed to random positions. Given that robots taking high spatial
resolution measurements can only be located at $m$ points, the locations
needed to be selected such that the informativeness of the sensing data
collected is maximized. To reiterate, using the definition of sensing regions
$S_{i}$ provided in Eq. (11), $L$ is a set containing the elements within the
sets $S_{1},S_{2},\dots,S_{m}$ that are arbitrary sensing regions, such that
for any two sets $S_{i}$ and $S_{j}$ with elements in $L$, $S_{i}\cap
S_{j}=\varnothing$. Let $\mathbf{x}_{L}$ contain the states only at the
locations specified by the set $L$, meaning $\mathbf{x}_{L}$ is just a vector
of sensor measurements. Furthermore, let $\mathbf{C}_{L}$ be the observation
matrix, such that $\forall p_{i}\in L$, ${\mathbf{C}_{L}}_{ii}=1$ and for all
other entries, ${\mathbf{C}_{L}}_{jk}=0$.

For sensor measurements with additive Gaussian white noise,
$\eta\sim\mathcal{N}(0,\sigma^{2})$, the relationship between the full state,
$\mathbf{x}$, and the available sensor measurements, $\mathbf{x}_{L}$ are
defined as

$\displaystyle\mathbf{x}_{L}$ $\displaystyle=\mathbf{C}_{L}\mathbf{x}+\eta$
(13)
$\displaystyle=\mathbf{C}_{L}\mathbf{\Phi}\mathbf{\Lambda}^{\top}\bm{\alpha}+\eta.$
(14)

When the full state is unknown, as in the case with collecting high-fidelity
sensor measurements at sparse locations, the full state must be predicted from
the observed state. In this case, the matrix $\mathbf{C}$ is fixed by our
choice of $L$, the sensing locations determined optimal as in Eq. (12). Then,
the optimal least-squares estimate for the full state is determined by
approximating the unknown quantity
$\mathbf{a}(t)=\mathbf{\Lambda}^{\top}\bm{\alpha}$ in Eq. (10). The optimal
least-squares then gives us
$\mathbf{a}(t)\approx(\mathbf{C}_{L}\mathbf{\Phi})^{\dagger}\mathbf{x}_{L}=\mathbf{\hat{a}}(t)$.
Substituting $\mathbf{\hat{a}}(t)$ back into the DMD estimation from Eq. (10)
gives the reconstructed full state $\mathbf{\hat{x}}$ from sensor measurements
$\mathbf{x}_{L}$ as

$\mathbf{\hat{x}}\approx\mathbf{\Phi}(\mathbf{C}_{L}\mathbf{\Phi})^{\dagger}\mathbf{x}_{L}.$
(15)

This procedure is known as the gappy POD method (Everson and Sirovich, 1995)
and has been well-studied in using sparse measurements to reconstruct a
process of interest from basis functions (Willcox, 2006; Manohar et al, 2019).

There is the central question of how to select the set $L$, and consequently
the matrix $\mathbf{C}_{L}$, such that it fulfills the objective of minimizing
the least-squares approximation error between an estimated
$\mathbf{a}(t)-\mathbf{\hat{a}}(t)$. One metric for estimating the least-
squares approximation error is the error covariance, $\mathbf{K}$, (Joshi and
Boyd, 2009; Manohar et al, 2018), where

$\mathbf{K}=Var(\mathbf{a}-\mathbf{\hat{a}})=\sigma^{2}\left[\left(\mathbf{C}_{L}\mathbf{\Phi}\right)^{\dagger}\mathbf{C}_{L}\mathbf{\Phi}\right]^{-1}.$
(16)

The error covariance $\mathbf{K}$ describes the $\rho$-confidence ellipsoid,
$\varepsilon_{\rho}$, that contains the least squares error
$\mathbf{a}-\mathbf{\hat{a}}$ with probability $\rho$ (Joshi and Boyd, 2009).
Intuitively, this means that the ellipsoid $\varepsilon_{\rho}$ contains all
the vectors $\mathbf{\hat{a}}$ that could be $\mathbf{a}$ with confidence
$\rho$. In this sense, we are seeking to minimize the volume of the ellipsoid,
defined as vol($\varepsilon_{\rho}$) = $\delta_{\rho,r}$det$\mathbf{K}^{1/2}$.

When picking the $p_{i}\in L$ optimal point locations from which sensor
measurements can be obtained, with $\mathbf{C}_{L}$ determined by $L$, this
optimization problem can be formulated as

$\displaystyle
L^{\star}=\operatorname*{argmax}_{L}\text{det}\left[\left(\mathbf{C}_{L}\mathbf{\Phi}\right)^{\dagger}\left(\mathbf{C}_{L}\mathbf{\Phi}\right)\right].$
(17)

We first note that the minimization of the determinant of $\mathbf{K}$, which
minimizes the volume of the error ellipsoid, is equivalent maximization of the
determinant to the inverse
$\left[\left(\mathbf{C}_{L}\mathbf{\Phi}\right)^{\dagger}\left(\mathbf{C}_{L}\mathbf{\Phi}\right)\right]$.
The optimal solution to this requires a combinatorial search of
$\Theta(N^{m})$ over the space of high-dimension, $N$, and $m$ candidate
sensing locations. However, the solution to Eq. (17) can be found via an
extremely efficient greedy optimization method using a modified the pivoted QR
factorization of matrices (Manohar et al, 2018, 2019). The pivoted QR
factorization is used to maximize the volume of successive submatrices,
meaning the absolute value of the determinant is also being maximized.

We first describe the traditional QR factorization method used for optimal
design, and then discuss our novel modifications that allow the method to be
generalized to our purposes. The traditional QR factorization differs from our
method in two regards: i) it does not explicitly offer a weighting scheme to
evaluate the quality of individual measurements, and ii) it does not take into
account sensing radius. In general, optimal sensor placement methods focus on
sensing at a fixed number of disjoint points. Alternatively, in this approach
we have robots with a finite sensing radius over a set of points, as in Eq.
(11). Thus, the optimization problem is jointly finding the best set of points
over all available robots to collect sensor measurements, as in Eq. (12). In
the modifications presented, we address both of these issues in a
computationally efficient way.

First, we present the traditional QR algorithm for optimal design, as shown in
Alg. 1. Then, we will present our pivoted QR factorization algorithm with
modifications made for selecting $m$ sensor locations based on sensing radius,
$k$, as shown in Alg. 2. We will show that these modifications admit a greedy
solution that maximizes the absolute value of the determinant, in Eq. (17),
while taking into account the sensing constraints of the robots.

(a)

(b)

(c)

(d)

Figure 6: Block pivoted QR factorization for optimal sensor location. (a) The
locations in the environment have a informativeness associated with collected
sensing at that point, where darker colors correspond to greater
informativeness. This information is capture in a projection matrix, with
columns of the corresponding to individual sensing locations. (b) The
maximally informative sensing location is selected as the the solid line
circle according to the maximization weighting scheme that selects the
associated optimal columns. (c) The matrix columns are pivoted and weights of
the remaining locations are adjusted by removing its orthogonal projection.
(d) A new maximally informative sensing location is selected as the dash line
circle and the process is repeated

In the traditional QR algorithm, QR factorization with column pivoting gives
us the decomposition
$\mathbf{\Phi}\mathbf{\Phi}^{\top}\mathbf{C}^{\top}=\mathbf{QR}$. The pivoting
provides an approximate solution for the optimization problem in Eq. (17).
This procedure is known as submatrix volume maximization since the matrix
volume is the absolute value of the determinant and at each iteration of the
algorithm, the absolute determinant is greedily maximized. In the QR column
pivoting procedure, a pivot column with the maximal 2-norm is selected that
corresponds to the optimal sensing point with respect to minimizing the error
ellipsoid. The orthogonal projection of all of the remaining columns onto this
column is subtracted from the remaining columns.

The method of using blocks within the pivoted QR algorithm admits a valid
factorization of the desired $\mathbf{\Phi}\mathbf{\Phi}^{\top}$ that allows
us to greedily optimize the absolute value of the desired matrix for a
specified sensing radius $k$. This blocked pivoted QR is demonstrated in Fig.
6. In this scenario, we can perform our block version of the pivoted QR
factorization of $\mathbf{\Phi}\mathbf{\Phi}^{\top}$ since these values will
coincide with the groups of singular values of the
$\left[\left(\mathbf{C}\mathbf{\Phi}\right)^{\top}\mathbf{C}\mathbf{\Phi}\right]^{-1}$
(Manohar et al, 2019). Pivoted QR factorization of matrices decomposes a
matrix $\mathbf{M}$ into a unitary matrix $\mathbf{Q}$ and an upper-triangular
matrix $\mathbf{R}$ such that
$\mathbf{\Phi}\mathbf{\Phi}^{\top}\mathbf{C}^{\top}=\mathbf{QR}$. Using the
properties of determinants, we have that
$|\text{det}\mathbf{\Phi}\mathbf{\Phi}^{\top}\mathbf{C}^{\top}|=|\text{det}\mathbf{Q}||\text{det}\mathbf{R}|=\prod_{i}|\mathbf{R}_{ii}|$.
The QR pivoting enforces the diagonal dominant structure on the diagonal
entries of the matrix $\mathbf{R}$ such that
$\sigma_{i}^{2}=|\mathbf{R}_{ii}|^{2}\geq\sum_{j=i}^{c}|\mathbf{R}_{jc}|,\quad
1\leq i\leq c\leq T+1$. This then maximizes the absolute value of the
determinant.

The block pivoted QR factorization evaluates all of the sets of columns that
correspond to candidate, potentially non-adjacent sensing locations from the
matrix $\mathbf{\Phi}\mathbf{\Phi}^{\top}$. The candidate sensing locations,
$\mathcal{C}$, are determined by first enumerating the set of all points in
each possible circle of radius $k$, corresponding to the sensing radius,
within the workspace $\mathcal{W}$. These candidate circles within
$\mathcal{W}$ are mapped to the columns corresponding to spatial locations in
the matrix. The optimal sensing location is the set of columns within the
matrix $\mathbf{\Phi}\mathbf{\Phi}^{\top}$ that correspond to the spatial
regions $S^{*}$ where the sum of the column 2-norms $\sum_{p\in
S_{i}}{||(\mathbf{\Phi}\mathbf{\Phi}^{\top})_{p}||_{2}}$ is maximized. The
algorithm then selects this optimal sensing location from the matrix and
pivots it to the beginning of the matrix, removing its orthogonal projection
from the remaining columns of matrix $\mathbf{\Phi}\mathbf{\Phi}^{\top}$. The
candidate sensing locations are pruned such that any sensing locations with
points in the selected optimal sensing location are discarded. This process is
repeated the until all optimal sensing locations have been found.

Algorithm 1 Pivoted QR factorization algorithm for determining optimal point
measurements

1:procedure Pivoted QR Factorization($\mathbf{\Phi}\mathbf{\Phi}^{\top}$, $m$)

2: $P\leftarrow[]$ $\triangleright$ Set of selected point measurements

3: $W\leftarrow[]$ $\triangleright$ Weights corresponding to points

4: for $i=1,\dots,k$ do

5: $p=argmax_{j\not\in P}||(\mathbf{\Phi}\mathbf{\Phi}^{\top})_{j}||_{2}$
$\triangleright$ $j$-th column

6: $P\leftarrow[P\quad p]$

7: Determine Householder matrix $H$

8: $\mathbf{\Phi}\mathbf{\Phi}^{\top}\leftarrow
diag(I,H)\mathbf{\Phi}\mathbf{\Phi}^{\top}$ $\triangleright$ Remove orthogonal
projection of $(\mathbf{\Phi}\mathbf{\Phi}^{\top})_{s}$ from all columns

9: return $P,W$

Algorithm 2 Pivoted QR factorization algorithm for determining optimal sensor
locations and their corresponding weights

1:procedure Pivoted QR Factorization($\mathbf{\Phi}\mathbf{\Phi}^{\top}$, $m$,
$\mathcal{W}$, $k$)

2: $\mathcal{C}\leftarrow$ create candidate sensing location($\mathcal{W}$,
$r$)

3: $L^{*}\leftarrow[]$ $\triangleright$ Set of optimal sensing locations

4: $W\leftarrow[]$ $\triangleright$ Weights corresponding to points

5: for $i=1,\dots,m$ do

6: $S^{*}=argmax_{S_{i}\in S,S_{i}\not\in L^{*}}\sum_{p\in
S_{i}}{||(\mathbf{\Phi}\mathbf{\Phi}^{\top})_{p}||_{2}}$

7: $L^{*}\leftarrow[L^{*}\quad S^{*}]$

8: $\mathcal{C}\leftarrow\textsc{remove overlapping
regions}(\mathcal{W},r,S^{*})$

9:
$w=\frac{\left\lVert(\mathbf{\Phi}\mathbf{\Phi}^{\top})_{S^{*}}\right\rVert_{2}}{\sum_{i=0}^{n}\left\lVert\mathbf{\Phi}\mathbf{\Phi}^{\top}\right\rVert_{2}}$

10: $W\leftarrow[W\quad w]$

11: Determine Householder matrix $H$

12: $\mathbf{\Phi}\mathbf{\Phi}^{\top}\leftarrow
diag(I,H)\mathbf{\Phi}\mathbf{\Phi}^{\top}$ $\triangleright$ Remove orthogonal
projection of $(\mathbf{\Phi}\mathbf{\Phi}^{\top})_{S^{*}}$ from all columns

13: return $S,W$

####  4.2.2 Coverage control using temporal characteristics

For robots with low spatial resolution measurements collected over time, a
coverage control technique can be used to determine the placement of the
robots.

In coverage control, the workspace $\mathcal{W}$ is partitioned into Voronoi
regions, such that each robot $j$ located at point $p_{j}$ is sensing the
points in the Voronoi region $V_{j}$ closest to its position, as in

$V_{j}=\\{q\in\mathcal{W}\,|\,\left\lVert q-p_{j}\right\rVert\leq\left\lVert
q-p_{i}\right\rVert,i\neq j\\},$ (18)

where $\mathcal{I}$ is the index set of all the points in the workspace. Let
$R$ be the set of locations of all $m^{\prime}$ robots for this sensing
modality.

The coverage control technique considers minimizing the function

$\mathcal{H}(R)=\sum\limits_{j=1}^{m^{\prime}}\int_{V_{j}}f\left(\left\lVert
q-p_{j}\right\rVert\right)\phi(q)dq.$ (19)

In coverage control, the density function $\phi$ is a measure of information
or probability that some event takes place. In this scenario, we can use the
temporal characteristics $\mathbf{\Lambda}$ as the density function $\phi$ in
Eq. (19) and use standard implements for coverage control techniques to
determine the optimal locations for these robots. The minimization of the
function Eq. (19) can be solved using centroidal Voronoi partitions. These
partitions correspond to the optimal partition of the space according to the
density function such that the best coverage is achieved. We refer the
interested reader to the work by Cortés et al (2004) for full details on
coverage control methods and implementation. The optimal sensing locations of
the vehicles collecting low-fidelity sensor measurements are determined by
partitioning the region according to the centers of the centroidal Voronoi
partitions.

###  4.3 Online adaptation of spatial and temporal features

We proposed two methods to allow for the online adaptation of the DMD modes
and eigenvalues. First, we propose a general method that works on large
dimensional systems and provides a robust approximation of the DMD operator.
However, this framework requires the retention of all of the snapshots
collected. The second method presented allows for fast computations and does
not require the retention of all of the snapshots collected, but imposes
restrictions on the workspace. The latter framework is amenable to long-term
monitoring of environmental phenomena.

In both procedures, we want to incorporate updates to the eigendecomposition
of the operator $\mathbf{A}$ after the acquisition of new data
$\mathbf{x_{T+1}}$, $\mathbf{x_{T+2}},\dots\mathbf{x_{T+1+\tau}}$. Define
$\mathbf{X_{new}}$ and $\mathbf{Y_{new}}$ as

$\displaystyle\mathbf{X_{new}}=[\mathbf{x_{T+1}}\,\mathbf{x_{T+2}}\,\dots\,\mathbf{x_{T+\tau}}]$
(20)
$\displaystyle\mathbf{Y_{new}}=[\mathbf{x_{T+2}}\,\mathbf{x_{T+3}}\,\dots\,\mathbf{x_{T+1+\tau}}]$
(21)

with the combined data in matrices

$\displaystyle\mathbf{X^{\prime}}=[\mathbf{X}\quad\mathbf{X_{new}}]$ (22)
$\displaystyle\mathbf{Y^{\prime}}=[\mathbf{Y}\quad\mathbf{Y_{new}}].$ (23)

Then, for the online adaptation, the desired decomposition is for the operator
$\mathbf{A}^{\prime}$ for the equation
$\mathbf{Y}^{\prime}=\mathbf{A}^{\prime}\mathbf{X}^{\prime}$.

####  4.3.1 Generalized adaptation for varying spatial and temporal scale
modeling

This method allows for adapting the eigenmodes and eigenvectors of
$\mathbf{A}$ to $\mathbf{A}^{\prime}$ in the presence of new data collected
from a variety of systems. This method is referred to as the “general” method.
While this method works with different spatial and temporal scales, it
requires storing the previous SVD of the data matrix $\mathbf{X}$ and all of
the snapshots of the data. However, as in the case of the traditional DMD
algorithm, this online adaptation scheme alleviates the need for directly
computing $\mathbf{A}^{\prime}$ and instead computes the eigendecomposition of
an approximation of $\mathbf{A}^{\prime}$. We generalize the techniques in
(Matsumoto and Indinger, 2017), which allow for single updates at a time,
using the results from (Brand, 2002) that allow for updating SVD computations
with missing or incomplete data. We can formulate a new matrix
$\mathbf{X^{\prime}}$ in terms of the SVD of the previous data matrix
$\mathbf{X}$ from Eq. (5) and the new data $\mathbf{X_{new}}$.

$\begin{split}\mathbf{X^{\prime}}&=[\mathbf{X}\quad\mathbf{X_{new}}]=[\mathbf{U}\bm{\Sigma}\mathbf{W}^{\top}\quad\mathbf{X_{new}}]\\\
&=[\mathbf{U}\quad\mathbf{J}]\begin{pmatrix}\bm{\Sigma}&\mathbf{L}\\\
\mathbf{0}&\mathbf{P}\end{pmatrix}\begin{pmatrix}\mathbf{W}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}\end{pmatrix}^{T}\\\
&=[\mathbf{U}(\mathbf{I}-\mathbf{UU}^{\top})/\mathbf{P}]\begin{pmatrix}\bm{\Sigma}&\mathbf{U}^{\top}\mathbf{X_{new}}\\\
\mathbf{0}&\mathbf{P}\end{pmatrix}\begin{pmatrix}\mathbf{W}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}\end{pmatrix}^{T}\end{split}$ (24)

where $\mathbf{L}=\mathbf{U}^{\top}\mathbf{X_{new}}$ is the projection of the
new data onto the orthogonal basis $\mathbf{U}$,
$\mathbf{P}=\mathbf{J}^{\top}\mathbf{H}$,
$\mathbf{H}=\mathbf{U}(\mathbf{I}-\mathbf{UU}^{\top})\mathbf{X_{new}}=\mathbf{X_{new}}-\mathbf{UL}$
is the component of the new data orthogonal to the subspace spanned by
$\mathbf{U}$, $\mathbf{J}$ is an orthogonal basis of $\mathbf{H}$ from the QR
decomposition $\mathbf{JR}=\mathbf{H}$, and
$\mathbf{K}=\mathbf{J}^{\top}\mathbf{H}$ be the projection of the new data
onto the subspace orthogonal to $\mathbf{H}$.

Denote an intermediate matrix $\mathbf{Z}$ as

$\mathbf{Z}=\begin{pmatrix}\bm{\Sigma}&\mathbf{U}^{\top}\mathbf{X_{new}}\\\
\mathbf{0}&\mathbf{P}\end{pmatrix}.$ (25)

We can then take the SVD of
$\mathbf{Z}=\mathbf{U^{\prime}}\bm{\Sigma^{\prime}}\mathbf{W^{\prime T}}$

Finally, we can update $\mathbf{U}$, $\bm{\Sigma}$, $\mathbf{W}$.

$\displaystyle\mathbf{U}\leftarrow[\mathbf{U}\quad\mathbf{J}]\mathbf{U^{\prime}}$
$\displaystyle\bm{\Sigma}\leftarrow\bm{\Sigma^{\prime}}$
$\displaystyle\mathbf{W}\leftarrow\begin{pmatrix}\mathbf{W}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}\end{pmatrix}\mathbf{W^{\prime}}$

With the updated $\mathbf{U}$, $\bm{\Sigma}$, $\mathbf{W}$, we can proceed
with Eq. (7) using the updated SVD matrices and the updated
$\mathbf{Y^{\prime}}$. This allows us to compute a low-dimensional
approximation of the $\mathbf{A}^{\prime}$ operator.

####  4.3.2 Fast adaptation for long-term modeling of environmental phenomena
with coarse spatial scale

Alternatively, this method allows for direct computation of the operator
$\mathbf{A}^{\prime}$ on data collected over an extended period of time.
However, the spatial dimension of the environment being sampled should be
small (less than 200). While the previous $\mathbf{A}$ matrix needs to be
retained for updating, the snapshots do not. Instead, a compact, lower
dimensional representation of the previous data is stored and updated. This
method is referred to as the “long-term” method. We extend and generalize the
results of (Zhang et al, 2019), where the authors allow for estimating updates
to the new $\mathbf{A}^{\prime}$ matrix after the acquisition of a single new
training datum. In this work, we allow for simultaneous batch updates to
compute the new $\mathbf{A}^{\prime}$ matrix.

Assuming that $\mathbf{X}$ is a matrix with full row rank, we can write the
pseudoinverse of the matrix as

$(\mathbf{X})^{\dagger}=(\mathbf{X})^{\top}(\mathbf{X}\mathbf{X}^{\top})^{-1}$
(26)

Thus, $\mathbf{A}$ can be written in terms of

$\displaystyle\mathbf{A}=\mathbf{Y}(\mathbf{X})^{\top}(\mathbf{X}\mathbf{X}^{\top})^{-1}=\mathbf{Q}\mathbf{S}$
(27) $\displaystyle\mathbf{Q}=\mathbf{Y}(\mathbf{X})^{\top}$ (28)
$\displaystyle\mathbf{S}=(\mathbf{X}\mathbf{X}^{\top})^{-1}.$ (29)

In order to calculate the updated $\mathbf{A}$, we need $\mathbf{A^{\prime}}$
= $\mathbf{Q^{\prime}}\mathbf{S^{\prime}}$, where we can write
$\mathbf{Q^{\prime}}$ and $\mathbf{S^{\prime}}$ in terms of $\mathbf{Q}$,
$\mathbf{S}$, and the new training data in $\mathbf{X_{new}}$ and
$\mathbf{Y_{new}}$.

$\displaystyle\mathbf{Q^{\prime}}$
$\displaystyle=[\mathbf{Y}\quad\mathbf{Y_{new}}][\mathbf{X}\quad\mathbf{X_{new}}]^{\top}$
$\displaystyle=\mathbf{Y}\mathbf{X}+\mathbf{Y_{new}}\mathbf{X_{new}^{\top}}$
$\displaystyle=\mathbf{Q}\mathbf{S}+\mathbf{Y_{new}}\mathbf{X_{new}^{\top}}$
(30) $\displaystyle(\mathbf{S^{\prime}})^{-1}$
$\displaystyle=[\mathbf{X}\quad\mathbf{X_{new}}][\mathbf{X}\quad\mathbf{X_{new}}]^{\top}$
$\displaystyle=\mathbf{X}\mathbf{X}^{\top}++\mathbf{X_{new}}\mathbf{X_{new}}^{\top}$
$\displaystyle=\mathbf{S^{-1}}+\mathbf{X_{new}}\mathbf{X_{new}}^{\top}$ (31)

The new $\mathbf{A^{\prime}}$ matrix is then

$\displaystyle\mathbf{A^{\prime}}$
$\displaystyle=\mathbf{Q^{\prime}}\mathbf{S^{\prime}}$
$\displaystyle=(\mathbf{Q}\mathbf{S}+\mathbf{Y_{new}}\mathbf{X_{new}^{\top}})(\mathbf{S}^{-1}+\mathbf{X_{new}}\mathbf{X_{new}^{\top}})^{-1}.$
(32)

From Eq. (32), we see that $\mathbf{Q^{\prime}}$ is easily computed using the
previously computed matrices $\mathbf{Q}$ and $\mathbf{S}$ and the newly
acquired data. However, we still need to compute
$\mathbf{S^{\prime}}=(\mathbf{S^{-1}}+\mathbf{X_{new}}\mathbf{X_{new}^{\top}})^{-1}$.
To do so, we use the Woodbury matrix identity given by

$(\mathbf{C}+\mathbf{DEF})^{-1}=\mathbf{C}^{-1}-\mathbf{C}^{-1}\mathbf{D}(\mathbf{E}^{-1}+\mathbf{FC}^{-1}\mathbf{D})^{-1}\mathbf{F}\mathbf{C}^{-1}.$
(33)

Let
$\mathbf{C}=\mathbf{S}^{-1},\mathbf{D}=\mathbf{X_{new}},\mathbf{E}=\mathbf{I},\mathbf{F}=\mathbf{X_{new}^{\top}}$
in Eq. (33), we can rewrite $\mathbf{S^{\prime}}$ as

$\displaystyle\mathbf{S^{\prime}}$
$\displaystyle=(\mathbf{S}^{-1}+\mathbf{X_{new}}\mathbf{X_{new}^{\top}})^{-1}$
$\displaystyle=\mathbf{S}-\mathbf{S}\mathbf{X_{new}}(\mathbf{I}+\mathbf{X_{new}^{\top}}\mathbf{S}\mathbf{X_{new}})^{-1}\mathbf{X_{new}^{\top}}\mathbf{S}.$
(34)

Let
$\mathbf{\Gamma}=(\mathbf{I}+\mathbf{X_{new}^{\top}}\mathbf{S}\mathbf{X_{new}})^{-1}$,
combining Eqs. (32) and (4.3.2), we can arrive at a simplified expression for
$\mathbf{A^{\prime}}$:

$\mathbf{A^{\prime}}=\mathbf{A}+(\mathbf{Y_{new}}-\mathbf{A}\mathbf{X_{new}})\mathbf{\Gamma}\mathbf{X_{new}^{\top}}\mathbf{S}.$
(35)

In this scenario, we also impose a weighting scheme to place greater
importance on training data that has been collected more recently, which is
helpful in the context of long-term monitoring where early training data may
not be as useful. To weight samples with varying importance, we apply a
forgetting factor $\gamma$ to past training data, receiving new expressions
for $\mathbf{Q^{\prime}}$ and $\mathbf{S^{\prime}}$ given by

$\displaystyle\mathbf{Q^{\prime}}=\gamma\mathbf{Q}+\mathbf{Y_{new}}\mathbf{X_{new}^{\top}}$
(36)
$\displaystyle(\mathbf{S^{\prime}})^{-1}=\gamma\mathbf{S^{-1}}+\mathbf{X_{new}^{\top}}\mathbf{X_{new}}.$
(37)

Accounting for the weighting scheme, we modify the computation of
$\mathbf{S^{\prime}}$ as

$\mathbf{S^{\prime}}=\frac{1}{\gamma}*(\mathbf{S}-\mathbf{S}\mathbf{X_{new}}\mathbf{\Gamma}\mathbf{X_{new}^{\top}}\mathbf{S}).$
(38)

Thus, we can updated the operator $\mathbf{A}$ explicitly and use an
eigenvalue decomposition to derive the DMD modes.

##  5 Analysis of online algorithms

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7: Comparison of various online DMD methods for computing eigenvalues. The data is collected in a $20\times 20$ environment over $1000$ time steps from $0$ to $10$ seconds, where the online algorithms are initialized with data from $400$ time steps at $4$ seconds and updated every $10$ time steps. Various noise regimes are studied for additive white Gaussian noise with low, medium, and high variance. (a)-(c) Estimation of real part of eigenvalue of dynamical system, where true value is $-1$ where variance is $0.01$ in (a), $0.04$ in (b), and $0.1$ in (c). (d) - (f) Estimation of imaginary part of eigenvalue of dynamical system, where true value is $0$ where variance is $0.01$ in (d), $0.04$ in (e), and $0.1$ in (f) Figure 8: Comparison of various weights using online data assimilation with DMD methods with batch DMD algorithm Table 1: Comparison of Execution Times for Online DMD Algorithms on Various Environments Environmental Parameters |  Time Dimension $T$ $\tau$ | $500$ | $500$ | $1000$ | $2000$
---|---|---|---|---|---
$10\times 10$ | $10\times 10$ | $20\times 20$ | $20\times 10$
100 | 100 | 400 | 200
10 | 100 | 100 | 100
Online DMD Method (s) |  Batch General Long-term | 5.038 | 1.376 | 3.715 | 14.77
0.557 | 0.975 | 1.710 | 3.343
0.304 | 0.646 | 2.969 | 1.387
With Eigen. Computation (s) |  Batch General Long-term | 13.59 | 2.756 | 37.52 | 28.82
0.537 | 1.950 | 3.411 | 4.266
0.413 | 1.319 | 29.61 | 17.03

In this section, we evaluate the proposed online algorithms, as described in
Sect. 4.3, in various environments with different time scales and resolutions.
The two proposed approaches for online adaptation of the DMD modes and
eigenvalues, general adaptation for varying spatial and temporal scale
modeling and fast adaptation for long-term modeling of environmental phenomena
with coarse spatial scale, have their respective advantages and drawbacks. As
a basis for comparison, we will use the naive batch algorithm. At each update
period $\tau$, the batch algorithm computes an estimate $\mathbf{A}_{batch}$
with $\mathbf{Y}^{\prime}=\mathbf{A}_{batch}\mathbf{X}^{\prime}$ using the
pseudo-inverse operator. Then, an eigendecomposition is performed on
$\mathbf{A}_{batch}$ to estimate the DMD modes and eigenvalues. For low-
dimensional systems and short time scales, the performance of the systems are
similar. For an environment with dimension $20\times 20$ and a complex
spatiotemporal process simulated $1000$ time steps over $10$ seconds, the
dominant eigenvalue is such that $\mathrm{Im}(\lambda_{\text{DMD}})=0$ and
$\mathrm{Re}(\lambda_{\text{DMD}})=-1$. The data is simulated and injected
with Gaussian random noise $\mathcal{N}(0,\sigma^{2})$ for $\sigma^{2}=0.01$,
$0.04$, and $0.1$. The online DMD methods are first initialized with data
collected over $400$ time steps and updated every $10$ time steps. For each
$\sigma^{2}$ value, 10 trials are collected and the data is aggregated to
analyze their descriptive statistics. Fig. 7 shows the comparison of the naive
batch algorithm, the general adaptation method, and long-term adaptation
method. The boxplots represent the median, upper and lower quantiles, and
outliers for each online DMD method under low, medium, and high noise regimes.
While the online methods do not converge exactly to the eigenvalue, they
provide close, fast approximations when compared to the batch method. After
the initial DMD model is computed, convergence for all three methods occurs
after several updates, around $t=5s$. For all the three methods, the greater
the noise, the greater the spread of the estimates of the eigenvalue across
trials. The general method provides closer estimates of the real parts of the
eigenvalue, while the long-term method provides closer estimates of the
imaginary parts of the eigenvalue. The long-term methods exhibits greater
variance in its estimates compared to the general method.

The general adaptation scheme allows for computation of eigenmodes using the
low-rank operator $\hat{\mathbf{A}}$, as in Eq. (7). This method is
computationally efficient in a wide range of scenarios, but does require
storing the previous SVD of the data matrix $\mathbf{X}$ and all of the
snapshots of the data in order to compute the low-rank operator. There are
certain conditions under which this operator is unable to produce accurate
estimates of the operator $\hat{\mathbf{A}}$. As presented by Tu et al (2013),
the eigendecomposition of the operator $\hat{\mathbf{A}}$ produces the
eigenvalues and eigenvectors of $\mathbf{A}$ only when $\mathbf{Y}$ lies in
the span of $\mathbf{X}$. Otherwise, the eigenvalues and eigenvectors of the
low-dimensional operator are useless. While the authors suggest this is easily
fulfilled in dynamical systems, in reality, this is not a condition that is
known a priori and may lead to extremely poor estimates of the dynamical
system. For example, if a system is slowly varying between time steps where
$\mathbf{Y}$ is simply a time-shifted $\mathbf{X}$, it may very well be that
$\mathbf{Y}$ could be written as a linear combination of the matrix
$\mathbf{X}$ and therefore $\mathbf{Y}$ would lie in the span of $\mathbf{X}$.
In general, for nonlinear systems this should not be common as this would
suggest the matrix $\mathbf{Y}$ is not a simple linear combination of the
matrix $\mathbf{X}$. However, in practice, some downsampling in the time
domain may be necessary to use the algorithm.

The long-term modeling scheme directly computes the $\mathbf{A^{\prime}}$
operator and allows for the inclusion of a weighting parameter. The weighting
parameter, $\gamma$, is advantageous for disregarding the previous data that
may not capture the relevant dynamics, as shown in Fig. 8. Allowing for a
forgetting factor improves the estimate of the DMD spatial and temporal modes.
While this work was proposed by (Zhang et al, 2019), there does not exist, to
the best of our knowledge, a systematic way for picking $\gamma$ to ensure the
performance of the system. It has not been shown whether under certain
conditions, values of $\gamma$ allow for a convergence to the true DMD
operator. Though computing the operator $\mathbf{A^{\prime}}$ is advantageous
in terms of model accuracy, it comes with some caveats. First, this method
requires initialization with full rank matrices. This means that in order for
the DMD computation to be initialized at least $N$ data points will have to
have been collected for the stability of the algorithm. Otherwise, the matrix
will be ill-conditioned and the model will be erroneous. For high-dimensional
systems, this constraint is hard to fulfill. For example, for a $100\times
100$ environment, $N$ would equal $10000$, and thus, at least, $10000$ time
steps of data would be required for the algorithm. This may not be feasible in
all applications. It is for this reason that in Sect. 6, we were unable to
test the algorithm experimentally on a high-dimensional system projected onto
the water tank where we only had $\sim 100$ time steps of data. Second, there
is additional computational time associated with the construction of the
matrix and then its eigendecomposition, instead of directly computing its
eigenvalues and eigenvectors.

Table 1 demonstrates the execution times of the general and long-term online
adaptation methods compared with the batch algorithm. All of the timing
calculations were performed on the same synthetic data but with various
discretizations in time and dimension. The system is defined as

$f(x,y,t)=senh(x)senh(y)*1.9i^{-t},$

where $senh(z)=(e^{z}+e^{-z})/2$. The system is a hyperbolic sine with damped
oscillation. The data is simulated by using the equation and injecting
Gaussian random noise $\mathcal{N}(0,0.4)$ for each trial. Additionally, the
parameters $T$, the initial batch size of the data used to compute the DMD
model, and the parameter $\tau$, the batch size of the assimilated data, are
varied across environments. Note, the parameter $T$ depends on the size of the
dimension for the long-term algorithm, as mentioned earlier. For our
implementation of the online adaptation algorithms, we used the
eigendecomposition algorithms from the scipy Python library, and modified the
DMD functions from the PyDMD (Demo et al, 2018) library. All computations were
performed on a Dell Optiplex 9020 with a 3.40 GHz Intel i7-4770 CPU.

The batch algorithm and long-term algorithm compute estimates of the operator
$\mathbf{A}$ as specified by their respective algorithms. The general
algorithm, however, directly computes the eigenvalues and eigenmodes. For the
comparison of the online DMD methods with eigendecompositions, the batch and
long-term algorithms must perform an eigenvalue decompositions, while the
general algorithm merges the eigenvectors into a DMD basis. Thus, the batch
algorithm must incrementally perform both the pseudo-inverse operation and
eigendecomposition with increasingly larger matrices. Similarly, the long-term
algorithm must perform the eigendecomposition on increasingly larger matrices.

For small environments and short time scales, as in the case of the $10\times
10$ environments, there is negligible difference between the general and long-
term algorithms. Both of them outperform the batch algorithm, as they do not
have to perform a pseudoinverse and subsequent matrix decomposition operator
during each assimilation step. In this scenario, even though the long-term
algorithm requires explicit eigendecomposition of the estimated
$\mathbf{A^{\prime}}$ matrix, it avoids the repeated SVD computations required
of the general algorithm. This results in the faster performance of the long-
term algorithm compared to the general algorithm for environments with small
dimensions and short duration. The computations associated with the
eigendecomposition take less time than the various low-dimension matrix
operations for environments for small dimensions. For larger environments with
longer durations, such as the $20\times 20$ with $1000$ time steps, the
general adaptation method outperforms long-term adaptation method both with
and without the SVD computation. The online algorithm is not well-suited for
large environments. For smaller environments ($\sim 200$ spatial points) being
monitored over a long duration, such as the $20\times 10$ grid with $2000$
time steps, the long-term algorithm is faster than the general algorithm
without eigendecomposition. It is worth noting that the long-term algorithm
does not require storage of any past snapshots, while the general algorithm
does, which in the case of long-term monitoring can become prohibitively
large. For the general algorithm used for monitoring of high-dimensional
system in 6.1, this can be up to $3MB$, while for the long-term monitoring it
is only several hundred kilobytes.

##  6 Mixed reality experiments

###  6.1 Simulation of proposed approach on barium cloud data

(a)

(b)

(c)

Figure 9: Depiction of coverage control technique and sensing quality for
aerial vehicles. (a) Aerial vehicles are randomly assigned sensing to begin.
Locations are assigned to robots for coverage based on their proximity to the
robot. (b) The true field can be sensed at the spatial locations shown. (c)
The robots each take one sensing region, where the regions are split by the
black line. The further away a location within a robot’s sensing region is
from the robot’s position, the greater the noise in the data associated with
that position. Thus, quality of the sensor measurements deteriorates further
away from the position of the robots, as seen by the differences in colors
from robots’ sensing measurements and the true field

(a)

(b)

Figure 10: Depiction of assignment to sensing regions for marine vehicles. (a) The marine vehicles sense the true process at a greater spatial resolution than the aerial vehicles. If the marine vehicles had full coverage of the process, the resulting process would appear as such. (b) However, marine vehicles have limited sensing radius and do not have sensing data available for spatial locations outside of their regions. The circular regions represent the data being sensed by each of the robots, whereas the white space represents regions for which sensing data is unavailable. An estimate of the true process in the white space is constructed using the limited sensing data from robots and a model Table 2: Comparison of Mean-Squared Averaged Over Entire Spatial Field and Time for Various Simulated Artificial Plasma Clouds and Algorithm Parameters Environment | Data Sampling Rate |  |
---|---|---|---
Time | Space | Noise | MSE
AV | MV | AV | MV |  |
Cloud 1 | 1 | 5 | 12 | 4 | 0.15 | 1.014
1 | 5 | 9 | 3 | 0.12 | 0.9885
1 | 5 | 6 | 2 | 0.12 | 0.9709
Cloud 2 | 1 | 5 | 12 | 4 | 0.15 | 0.0968
1 | 5 | 6 | 2 | 0.12 | 0.0947
1 | 2 | 9 | 3 | 0.12 | 0.1235
Cloud 10 | 1 | 5 | 12 | 4 | 0.15 | 3.078
1 | 2 | 9 | 3 | 0.12 | 3.166

(a)

(b)

Figure 11: Error over time comparing models from heterogeneous data and
homogeneous data. Mean-squared error over entire spatial grid is calculated
over the time series for reconstruction from heterogeneous data, homogeneous
data from aerial vehicles, and homogeneous data from marine vehicles. (a) The
simulated environment, while nonlinear, is relatively simple compared to the
other tested environments. The model with heterogeneous outperforms both of
the homogeneous models for the duration of the process. (b) For a more
complicated process, the estimates from the aerial vehicle are better to
begin. However, as the process continues, the heterogeneous model estimates
result in lower error.

In this section, we show that using our proposed algorithm,the performance of
the heterogeneous system is better than that of a homogeneous system and that
the optimal sensing locations derived from our proposed algorithm impact the
overall quality of the model in complex environments. A simulated environment
is created using the density of an artificial plasma cloud in the near-Earth
space environment. The density of this artificial plasma cloud was derived in
simulation, using the models presented by Zalesak et al (1987). These
processes do not simply convect, but rather can become unstable and undergo
complex nonlinear evolution. When used as an environmental processes of
interest, the evolution of the density of the artificial plasma cloud allows
us to study the ability of the robots to track an unstable, complex nonlinear
phenomena.

We compare various clouds on different algorithm parameters and study their
average over time. The clouds are computed using computation fluid dynamics
models of artificial plasma clouds that serve as representation for barium
clouds in the atmosphere. These models were initialized under various
parameters from the equations in Zalesak et al (1987). To begin, the simulated
data is collected and interpolated by the simulated aerial vehicles for a
certain predetermined number of time steps. This is used to create an initial
model using DMD. The marine robots’ positions are randomly initialized. Then,
the sensing data from the marine robots is used to estimate the full field and
combined into the model. The spatial and temporal characteristics are
extracted from the model. Based on this initial information, the robots are
sent to new optimal sensing locations. This is done using the methods
described in Sections 4.2.1 and 4.2.2 for the marine vehicles, and aerial
vehicles, respectively. The data from the aerial vehicle is collected for the
next predetermined number of time steps, at which point the data is combined
with the high-fidelity sensing data and estimates from the marine vehicle. All
of this data is assimilated into the original model using the generalized
adaptation techniques. At this point, there is new spatial and temporal
characteristics extracted. The process is repeated for the whole time series,
at which point the algorithm is terminated.

First, we assume two aerial vehicles have full coverage and five marine
vehicles each have a sensing radius of $15$ units. It is assumed aerial
vehicles and marine vehicles are collecting sensor measurements using the
paradigms shown in Fig. 9 and Fig. 10, respectively. The mean-squared error
over the entire spatial grid calculated over the time series as $mse$ is
calculated for field estimate $\bm{\hat{x}(t)}$ using
$mse=\frac{1}{N}\sum_{i=1}^{N}\left(\bm{\hat{x}}_{i}(t)-\bm{{x}}_{i}(t)\right)^{2}$.
In this case, the model estimates at each location $i$, $\bm{\hat{x}}_{i}(t)$
is compared against the value of the density from the simulated data,
$\bm{{x}}_{i}(t)$. This is then averaged over all time steps. The results from
these simulations are shown in Table 2. The data sampling rates are varied for
the aerial vehicle (AV) and marine vehicles (MV). The data sampling rate is an
number integer number, $n$, indicates that every $n^{th}$ element was taken
from the full simulated data in either the time or space domain. The time
sampling rate for the marine vehicle also indicates how after the model is
updated. For example for the MV time sampling rate $5$, the model is updated
every $5$ time steps, where the marine vehicle only has data at the $5^{th}$
time step. The same noise is injected into both the aerial vehicle and marine
vehicle sensor data.

Using the results from the Cloud $1$ simulation, we see higher noise results
in worse estimates, as expected. We also see that greater spatial resolution
leads to better results. However, for faster update times, seen in the results
from Cloud $2$, where the MV time sampling rate is $2$, we observe worse
results from an increase in time sampling. This is due to the fact that early
models of the process from few data points are poor and these errors from the
early propagate throughout the duration of the process. We also observe very
low errors in Cloud $2$, when compares to Clouds $1$ and $10$. This is due to
the relative complexity of Cloud $2$ to Clouds $1$ and $10$. Cloud $2$
exhibits much simpler behavior throughout the process.

We also observe the model estimate over the duration of the full process and
compare them to models created using just homogeneous data. The estimates of
the field are computing by using data from just the aerial vehicles, data from
just the marine vehicles, and data combined from both types of vehicles. The
mean-squared error is computed over the entire spatial grid for each time
step, as shown in Fig. 11. This corresponds to the simulation in Table 2 for
row $1$ in the environment using Cloud $2$ in (a) and row $1$ in the
environment using Cloud $2$ in (b). The estimates using the data from the
aerial vehicles are computed using standard interpolating techniques. The
estimates using the data from the marine vehicles are computed by constructing
a DMD model from the model, as in Eqs. (8) and (9), and using prediction
techniques for times when sensing data is unavailable, as in Eq. (10). The
estimates using data combined from both types of vehicles are computed as
described above. Note, that for the DMD model computed using homogeneous
marine vehicle data, we assume the homogeneous data is collected from the
optimal sensing locations generated by our model. Using random locations
results in extremely high errors. The homogeneous data from the aerial
vehicles is noisy and collected at a much lower spatial resolution then the
true process. As the process becomes more complex, the inclusion of multiple
types of data allows the proposed approach to outperform either of the other
estimations. This is especially evident for Fig. 11 (b), where towards the end
of the process the optimal sensing locations are critical in improving the
model performance for the model using heterogeneous data.

###  6.2 Experimental validation in water tank test bed

In the following, we present results to show our method successfully uses a
heterogeneous team of robots to model and infer the properties of a complex
spatiotemporal phenomena. We evaluate this claim for the general adaptation
strategy presented above and also present results on a test bed of real robot
to validate the use of the development of these algorithms in physical
environments.

(a)

(b)

Figure 12: Experimental setup used for testing algorithms. (a) Micro-
autonomous surface vehicles collecting high-fidelity sensor measurements. (b)
Indoor $4.5m\times 3.0m\times 1.2m$ water tank equipped with motion capture
system for multi-robot experiments

To model and predict properties of a complex, nonlinear process, we use a
fleet of three real micro-autonomous surface vehicles (mASVs), one simulated
mASV, and two simulated aerial vehicles. The mASVs are differential drive and
include a micro-controller board, an XBee radio module, and an inertial
measurement unit. The XBee radio modules are used to communicate the
information collected onboard. The testbed tank is $4.5m\times 3.0m\times
1.2m$ and is equipped within a motion capture system used for localization.
The mASVS are all localized using this motion capture system. This
experimental setup is shown in Fig. 12.

We ran experiments on $100$ seconds of data to evaluate our general online
adaptation strategy in Sect. 4.3 and using the environmental setup described
in Sect. 6.1. The process we are tracking is the density of an artificial
plasma cloud in the near-Earth space environment, derived in the same way as
the simulated environment described in Sect. 6.1 and simulated using
computational fluid dynamics models based on the equations in Zalesak et al
(1987). This process is projected onto the tank as an image, where the pixel
values of the image are mapped to concentration values of the density of the
artificial plasma cloud for the aerial vehicle. We use simulated sensing data
as noisy sparse measurements from the concentration data for marine vehicles.
In simulation, travel time is not explicitly accounted for in the
determination of optimal sensing locations. Instead, it is assumed that the
robots will be able to travel to these locations before the next adaptation
time. This is validated in these experiments, where robots are in fact able to
determine optimal sensing locations, assign robots in the team to locations
based on minimizing cumulative time, travel to optimal sensing locations, and
collect sensor measurements all before the next adaptation time.

We evaluate the full framework with the generalized adaptation framework. The
environment is originally defined over a $384\times 768$ workspace. In this
environment, we have four marine robots and two aerial vehicles. The workspace
is down-sampled to a $10\times 20$ grid for the low-fidelity data. This low-
fidelity data is simulated as on-board sensing data for the aerial vehicles.
Aerial vehicles are assumed to have full views of their sensing region as
determined by a coverage control method, where the quality of their sensing
data diminishes at points further away from their current sensing location, as
in Fig. 9. This is simulated by calculating the desired position of the aerial
vehicles from the coverage control technique and adding Gaussian random noise
to the true sensing data as a function of the distance from the desired
position. This data is then interpolated to the dimensions of the high-
fidelity data collected by the marine robots which is on a $48\times 96$ grid.
The marine robots have a sensing radius of $10$, meaning that within the
$48\times 96$ grid, a marine robot can collect sensing information over all of
the grid points laying in a circle with radius $10$ centered at its position,
as in Fig. 10. The high-fidelity sensing data from the marine vehicles is
combined into a single vector $\mathbf{x}_{L}$ and used to extrapolate the
data from the $48\times 96$ grid using Eq. (15), based off of Gappy POD
methods (Everson and Sirovich, 1995; Manohar et al, 2018). While the aerial
vehicles are able to collect sensor measurements across the full time series,
marine vehicles are only able to collect data periodically at regular
intervals.

(a)

(b)

(c)

(d)

Figure 13: Depiction of field reconstruction using adapted DMD models and
sensor measurements from robots. (a) True field at full spatial resolution as
obtained by artificial plasma cloud simulation at time step $40$. (b)
Reconstructed of field using heterogeneous data at time step $40$, with marine
robot locations depicted by black circles. (c) True field at full spatial
resolution as obtained by artificial plasma cloud simulation at time step
$40$. (d) Reconstruction of field using heterogeneous data at time step $40$,
with updated marine robot locations

An example of the robots tracking an artificial plasma cloud is shown in Fig.
13. In the beginning of the experiment, as in Fig. 13 (a), the robots are able
to create an reconstruction of the field, shown by Fig. (b). It should be
noted, the robots are taking sensor measurements in the areas of the process
such that they are not collecting varied information, as seen in Fig. (b).
However, as the experiment continues and the robots are able to collect more
data, the model is improved, such that the optimal sensing locations
correspond to the areas of interesting dynamics in the process, as evident in
Fig. 13 (c) and (d). It is evident that the reconstruction is able to capture
the complex, nonlinear dynamics of the process, shown in Fig. 13 (d). More
specifically, we see robots collected sensor measurements at or near the
spatial locations that correspond to the complex phenomena in the process. The
robots in the mixed reality experiment are able to successfully track a
complex, nonlinear environmental process using their heterogeneous multiscale,
multiresolution data.

##  7 Conclusion

In this work, we contribute to the understanding of using heterogeneous robot
systems for fusion of information to the end of modeling complex and
multiscale nonlinear phenomena. We consider the problem of leveraging
multiresolution sensor information from a heterogeneous robot team to model
and predict the evolution of a spatiotemporal process. Using these models, we
infer the optimal sensing locations using computationally efficient methods.
We show that the collection of this information and adaptation of the models
with measurements from disparate sources results in accurate reconstructions
of the estimated field. The frameworks presented allows for a decoupling of
the temporal and spatial modes apparent in the data. This decoupling is then
used within a task allocation framework for the various types of robots.
Instead of relying on the standard task-trait allocation approaches typically
used in heterogeneous robotic framework, this approach leverages the unique
strengths of the robots to jointly complete a task. Additionally, the
framework proposes two distinct ways for assimilating online data. In the case
of data assimilation for long-term modeling, the forgetting factor is useful
in scenarios where the environmental process is quickly time varying. However,
the precise nature of the forgetting factor, specifically what the optimal
value is for a specific process, is not well understood and requires greater
investigation (Zhang et al, 2019). The reliance of this framework on DMD and
its inherent capability for future state-estimation, as in Eq. (10), may allow
for the use of other optimal sensing techniques that account for not only the
current dynamics of the system but also the future dynamics of the system.
Much of the literature in the intersection of optimal sensing and robotics
rely on the use of GPs. However, in these scenarios, the spatial dimension of
the system, $N$ is extremely high ($\sim 10,000$ points). For modeling time
series data, GPs require knowledge of the kernel and the complexity of the
algorithm would be $\mathcal{O}^{N^{3}}$, which is prohibitively large in this
scenario. This is contrast to standard DMD methods that have few assumptions
and a computational complexity of $\mathcal{O}^{NTk}$ (Erichson et al, 2019).
This work presents and experimentally validates a framework to combine
heterogeneous sensing and mobility capabilities to learn, optimally sense with
respect to, and adapt a model.

## References

  * Parker (2003) Parker LE (2003) The Effect of Heterogeneity in Teams of 100+ Mobile Robots. Tech. rep.
  * Prorok et al (2015) Prorok A, Hsieh MA, Kumar V (2015) Fast redistribution of a swarm of heterogeneous robots. In: EAI International Conference on Bio-inspired Information and Communications Technologies (BICT), DOI 10.4108/eai.3-12-2015.2262349
  * Prorok et al (2016) Prorok A, Hsieh MA, Kumar V (2016) Formalizing the impact of diversity on performance in a heterogeneous swarm of robots. In: Proceedings - IEEE International Conference on Robotics and Automation, IEEE, vol 2016-June, pp 5364–5371, DOI 10.1109/ICRA.2016.7487748
  * Vasilijevic et al (2015) Vasilijevic A, Calado P, Lopez-Castejon F, Hayes D, Stilinovic N, Nad D, Mandic F, Dias P, Gomes J, Molina JC, Guerrero A, Gilabert J, Miskovic N, Vukic Z, Sousa J, Georgiou G (2015) Heterogeneous robotic system for underwater oil spill survey. _MTS/IEEE OCEANS 2015 - Genova: Discovering Sustainable Ocean Energy for a New World_ pp 1–7, DOI 10.1109/OCEANS-Genova.2015.7271492
  * Manjanna et al (2018) Manjanna S, Li AQ, Smith RN, Rekleitis I, Dudek G (2018) Heterogeneous Multi-Robot System for Exploration and Strategic Water Sampling. _Proceedings - IEEE International Conference on Robotics and Automation_ pp 4873–4880, DOI 10.1109/ICRA.2018.8460759
  * Manderson et al (2019) Manderson T, Manjanna S, Dudek G (2019) Heterogeneous Robot Teams for Informative Sampling. _Workshop on Informative Path Planning and Adaptive Sampling at Robotics Science and Systems_ URL http://arxiv.org/abs/1906.07208
  * Maini et al (2018) Maini P, Gupta G, Tokekar P, Sujit P (2018) Visibility-Based Monitoring of a Path Using a Heterogeneous Robot Team. In: IEEE International Conference on Intelligent Robots and Systems, pp 3765–3770, DOI 10.1109/IROS.2018.8593960
  * Dunbabin and Marques (2012) Dunbabin M, Marques L (2012) Robots for environmental monitoring: Significant advancements and applications. _IEEE Robotics and Automation Magazine_ 19(1):24–39, DOI 10.1109/MRA.2011.2181683
  * Singh et al (2010) Singh A, Ramos F, Durrant Whyte H, Kaiser WJ, Whyte HD, Kaiser WJ (2010) Modeling and Decision Making in Spatio-Temporal Processes for Environmental Surveillance Amarjeet. _Proceedings - IEEE International Conference on Robotics and Automation_ pp 5490–5497, DOI 10.1109/ROBOT.2010.5509934
  * Strogatz (2001) Strogatz S (2001) _Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering, Second Edition (Studies in Nonlinearity)_ , vol 32. Westview Press, DOI 10.5860/choice.32-0994
  * Korsah et al (2013) Korsah GA, Stentz A, Dias MB (2013) A comprehensive taxonomy for multi-robot task allocation. _International Journal of Robotics Research_ 32(12):1495–1512, DOI 10.1177/0278364913496484
  * Jones et al (2006) Jones EG, Browning B, Dias MB, Argall B, Veloso M, Stentz A (2006) Dynamically formed heterogeneous robot teams performing tightly-coordinated tasks. _Proceedings - IEEE International Conference on Robotics and Automation_ 2006(May):570–575, DOI 10.1109/ROBOT.2006.1641771
  * Rossi et al (2009) Rossi C, Aldama L, Barrientos A (2009) Simultaneous task subdivision and allocation for teams of heterogeneous robots. _Proceedings - IEEE International Conference on Robotics and Automation_ pp 946–951, DOI 10.1109/ROBOT.2009.5152299
  * Bae et al (2019) Bae J, Lee J, Chung W (2019) A heuristic for task allocation and routing of heterogeneous robots while minimizing maximum travel cost. _Proceedings \- IEEE International Conference on Robotics and Automation_ 2019-May:4531–4537, DOI 10.1109/ICRA.2019.8794257
  * Liang et al (2019) Liang CD, Wang L, Yao XY, Liu ZW, Ge MF (2019) Multi-target tracking of networked heterogeneous collaborative robots in task space. _Nonlinear Dynamics_ 97(2):1159–1173, DOI 10.1007/s11071-019-05038-x, URL https://doi.org/10.1007/s11071-019-05038-x
  * Notomista et al (2019) Notomista G, Mayya S, Hutchinson S, Egerstedt M (2019) An optimal task allocation strategy for heterogeneous multi-robot systems. _2019 18th European Control Conference, ECC 2019_ pp 2071–2076, DOI 10.23919/ECC.2019.8795895
  * Gerkey and Matarić (2004) Gerkey BP, Matarić MJ (2004) A formal analysis and taxonomy of task allocation in multi-robot systems. _International Journal of Robotics Research_ 23(9):939–954, DOI 10.1177/0278364904045564
  * Khamis et al (2015) Khamis A, Hussein A, Elmogy A (2015) Multi-robot task allocation: A review of the state-of-the-art. In: Cooperative Robots and Sensor Networks 2015, Springer, pp 31–51
  * Nashashibi et al (1992) Nashashibi F, Devy M, Fillatreau P (1992) Indoor scene terrain modeling using multiple range images for autonomous mobile robots. In: Proceedings - IEEE International Conference on Robotics and Automation, pp 40–46
  * Cortés et al (2004) Cortés J, Martínez S, Karatas T, Bullo F, Member S (2004) Coverage Control for Mobile Sensing Networks. _IEEE Transactions on Robotics and Automation_ 20(2):243–255, DOI 10.1109/TRA.2004.824698
  * Santos et al (2018) Santos M, Diaz-Mercado Y, Egerstedt M (2018) Coverage Control for Multirobot Teams with Heterogeneous Sensing Capabilities. _IEEE Robotics and Automation Letters_ 3(2):919–925, DOI 10.1109/LRA.2018.2792698
  * Santos and Egerstedt (2018) Santos M, Egerstedt M (2018) Coverage Control for Multi-Robot Teams with Heterogeneous Sensing Capabilities Using Limited Communications. _IEEE International Conference on Intelligent Robots and Systems_ pp 5313–5319, DOI 10.1109/IROS.2018.8594056
  * Sadeghi and Smith (2019) Sadeghi A, Smith SL (2019) Coverage control for multiple event types with heterogeneous robots. _Proceedings - IEEE International Conference on Robotics and Automation_ 2019-May:3377–3383, DOI 10.1109/ICRA.2019.8793639
  * Julian et al (2012) Julian BJ, Angermann M, Schwager M, Rus D (2012) Distributed robotic sensor networks: An information-theoretic approach. _International Journal of Robotics Research_ 31(10):1134–1154, DOI 10.1177/0278364912452675
  * Krause and Guestrin (2007) Krause A, Guestrin C (2007) Nonmyopic active learning of Gaussian processes: An exploration- exploitation approach. _ACM International Conference Proceeding Series_ 227:449–456, DOI 10.1145/1273496.1273553
  * Krause et al (2008) Krause A, Singh A, Guestrin C (2008) Near-optimal sensor placements in Gaussian processes: Theory, efficient algorithms and empirical studies. _Journal of Machine Learning Research_ 9:235–284, DOI 10.1145/1390681.1390689
  * Singh et al (2009) Singh A, Krause A, Guestrin C, Kaiser WJ (2009) Efficient Informative Sensing using Multiple Robots. _Journal of Artificial Intelligence Research_ 34:707–755, DOI 10.1613/jair.2674
  * Xu et al (2011) Xu Y, Choi J, Oh S (2011) Mobile sensor network navigation using Gaussian processes with truncated observations. _IEEE Transactions on Robotics_ 27(6):1118–1131, DOI 10.1109/TRO.2011.2162766
  * Rasmussen (2004) Rasmussen CE (2004) Gaussian Processes in machine learning. _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ 3176(Advanced Lectures on Machine Learning):63–71, DOI 10.1007/978-3-540-28650-9_4
  * Taira et al (2017) Taira K, Brunton SL, Dawson ST, Rowley CW, Colonius T, McKeon BJ, Schmidt OT, Gordeyev S, Theofilis V, Ukeiley LS (2017) Modal analysis of fluid flows: An overview. _AIAA Journal_ 55(12):4013–4041, DOI 10.2514/1.J056060
  * Mezić (2005) Mezić I (2005) Spectral properties of dynamical systems, model reduction and decompositions. _Nonlinear Dynamics_ pp 309–325, DOI 10.1007/s11071-005-2824-x
  * Rowley et al (2009) Rowley CW, Mezi I, Bagheri S, Schlatter P, Henningson DS (2009) Spectral analysis of nonlinear flows. _Journal of Fluid Mechanics_ 641:115–127, DOI 10.1017/S0022112009992059
  * Schmid (2010) Schmid PJ (2010) Dynamic mode decomposition of numerical and experimental data. _Journal of Fluid Mechanics_ 656:5–28, DOI 10.1017/S0022112010001217
  * Tu et al (2013) Tu JH, Rowley CW, Luchtenburg DM, Brunton SL, Kutz JN (2013) On dynamic mode decomposition - theory and applications. _Journal of computational dynamics_ (September):1–30, DOI 10.3934/jcd.2014.1.391
  * Jovanović et al (2014) Jovanović MR, Schmid PJ, Nichols JW (2014) Sparsity-promoting dynamic mode decomposition. _Physics of Fluids_ 26(2), DOI 10.1063/1.4863670
  * Brunton et al (2016) Brunton BW, Johnson LA, Ojemann JG, Kutz JN (2016) Extracting spatial-temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition. _Journal of Neuroscience Methods_ 258:1–15, DOI 10.1016/j.jneumeth.2015.10.010
  * Zhang et al (2019) Zhang H, Rowley CW, Deem EA, Cattafesta LN (2019) Online dynamic mode decomposition for time-varying systems. _SIAM Journal on Applied Dynamical Systems_ 18(3):1586–1609, DOI 10.1137/18M1192329, URL http://arxiv.org/abs/1707.02876
  * Matsumoto and Indinger (2017) Matsumoto D, Indinger T (2017) On-the-fly algorithm for Dynamic Mode Decomposition using Incremental Singular Value Decomposition and Total Least Squares, URL http://arxiv.org/abs/1703.11004
  * Manohar et al (2018) Manohar K, Brunton BW, Kutz JN, Brunton SL (2018) Data-Driven Sparse Sensor Placement for Reconstruction: Demonstrating the Benefits of Exploiting Known Patterns. _IEEE Control Systems_ 38(3):63–86, DOI 10.1109/MCS.2018.2810460
  * Manohar et al (2019) Manohar K, Kaiser E, Brunton SL, Kutz JN (2019) Optimized Sampling for Multiscale Dynamics. _Multiscale Modeling & Simulation_ 17(1):117–136, DOI 10.1137/17m1162366
  * Salam and Hsieh (2019) Salam T, Hsieh MA (2019) Adaptive Sampling and Reduced-Order Modeling of Dynamic Processes by Robot Teams. _IEEE Robotics and Automation Letters_ 4(2):477–484, DOI 10.1109/LRA.2019.2891475
  * Berger et al (2015) Berger E, Sastuba M, Vogt D, Jung B, Ben Amor H (2015) Estimation of perturbations in robotic behavior using dynamic mode decomposition. _Advanced Robotics_ DOI 10.1080/01691864.2014.981292
  * Folkestad et al (2020) Folkestad C, Pastor D, Mezic I, Mohr R, Fonoberova M, Burdick J (2020) Extended Dynamic Mode Decomposition with Learned Koopman Eigenfunctions for Prediction and Control. In: Proceedings of the American Control Conference, DOI 10.23919/ACC45564.2020.9147729
  * Andersone (2019) Andersone I (2019) Heterogeneous map merging: State of the art. _Robotics_ 8(3):1–29, DOI 10.3390/robotics8030074
  * Topal et al (2010) Topal S, Erkmen I, Erkmen AM (2010) A Novel Map Merging Methodology for Multi-Robot Systems. _World Congress on Engineering and Computer Science_ I
  * Park et al (2016) Park J, Sinclair AJ, Sherrill RE, Doucette EA, Curtis JW (2016) Map merging of rotated, corrupted, and different scale maps using rectangular features. In: Proceedings of the IEEE/ION Position, Location and Navigation Symposium, PLANS 2016, DOI 10.1109/PLANS.2016.7479743
  * Ma et al (2016) Ma L, Zhu J, Zhu L, Du S, Cui J (2016) Merging grid maps of different resolutions by scaling registration. _Robotica_ DOI 10.1017/S0263574715000168
  * Shahbandi and Magnusson (2019) Shahbandi SG, Magnusson M (2019) 2D map alignment with region decomposition. _Autonomous Robots_ DOI 10.1007/s10514-018-9785-7
  * Budišić et al (2012) Budišić M, Mohr R, Mezić I (2012) Applied Koopmanism. _Chaos_ 22(4), DOI 10.1063/1.4772195
  * Everson and Sirovich (1995) Everson R, Sirovich L (1995) Karhunen–Loève procedure for gappy data. _Journal of the Optical Society of America_ 12(8):1657, DOI 10.1364/JOSAA.12.001657
  * Willcox (2006) Willcox K (2006) Unsteady flow sensing and estimation via the gappy proper orthogonal decomposition. _Computers and Fluids_ 35(2):208–226, DOI 10.1016/j.compfluid.2004.11.006
  * Joshi and Boyd (2009) Joshi S, Boyd S (2009) Sensor selection via convex optimization. _IEEE Transactions on Signal Processing_ DOI 10.1109/TSP.2008.2007095
  * Brand (2002) Brand M (2002) Incremental singular value decomposition (SVD) of incomplete data. Tech. rep., Mitsubishi Electric Research Labs, URL http://www.merl.com
  * Demo et al (2018) Demo N, Tezzele M, Rozza G (2018) PyDMD: Python Dynamic Mode Decomposition. _The Journal of Open Source Software_ 3(22):530, DOI https://doi.org/10.21105/joss.00530
  * Zalesak et al (1987) Zalesak ST, Drake JF, Huba JD (1987) Dynamics of three dimensional ionospheric plasma clouds. _Physical Review Letters_ 58(3):278–281, DOI 10.1103/PhysRevLett.58.278
  * Erichson et al (2019) Erichson NB, Brunton SL, Kutz JN (2019) Compressed dynamic mode decomposition for background modeling. _Journal of Real-Time Image Processing_ DOI 10.1007/s11554-016-0655-2