Mobile Sensor Path Planning for Kalman Filter Spatiotemporal Estimation

The estimation of spatiotemporal data from limited sensor measurements is a required task across many scientific disciplines. In this paper, we consider the use of mobile sensors for estimating spatiotemporal data via Kalman filtering. The sensor selection problem, which aims to optimize the placement of sensors, leverages innovations in greedy algorithms and low-rank subspace projection to provide model-free, data-driven estimates. Alternatively, Kalman filter estimation balances model-based information and sparsely observed measurements to collectively make better estimation with limited sensors. It is especially important with mobile sensors to utilize historical measurements. We show that mobile sensing along dynamic trajectories can achieve the equivalent performance of a larger number of stationary sensors, with performance gains related to three distinct timescales: (i) the timescale of the spatiotemporal dynamics, (ii) the velocity of the sensors, and (iii) the rate of sampling. Taken together, these timescales strongly influence how well-conditioned the estimation task is. We draw connections between the Kalman filter performance and the observability of the state space model and propose a greedy path planning algorithm based on minimizing the condition number of the observability matrix. This approach has better scalability and computational efficiency compared to previous works. Through a series of examples of increasing complexity, we show that mobile sensing along our paths improves Kalman filter performance in terms of better limiting estimation and faster convergence. Moreover, it is particularly effective for spatiotemporal data that contain spatially localized structures, whose features are captured along dynamic trajectories.


Introduction and Related Work
Many scientific disciplines require the estimation of spatio-temporal data from limited, point-source sensor measurements for the purpose of characterization, forecasting, reconstructing, and/or controlling a given system.The goal of optimal sensor placement, also known as sensor selection, is to find the optimal locations in the state space to place only a few sensors so as to achieve the best performance in one or more of the above listed metrics.The combinatorial optimization problem of sensor selection is NP-hard, so most algorithms aim to find sub-optimal solutions by leveraging greedy searches and low-rank subspace representations of the system in order to efficiently find a near-optimal solution [1].Greedy methods [2], [3] are computationally efficient and include QR decomposition [4] with column pivoting [1], [5], [6], (Q)DEIM [7]- [9], and Gappy-POD [10]- [12], all of which take advantage of the sub-modularity, or near-submodularity, of criteria such as the trace, spectral norm, condition number, determinant and/or its low-rank projection basis.Greedy searches can also be modified to include cost constraints in the sensor placement problem [5], [13].Greedy solutions can also be further refined, such as through a genetic algorithm [14].Other objectives, such as the reconstruction error [15] and the observability matrix [16], can also be used for sensor selection.Statistical methods using Gaussian process models [17], [18] also are effective in leveraging entropy or mutual information as the main objective for optimization.And more recently, shallow decoder networks can be trained within the context of greedy algorithms [19], [20].
In contrast to instantaneous estimation from sensor measurements, Kalman filtering provides a recursive method that estimates based on collective information from prior knowledge of the dynamical model and a time-history of the sensor measurements [21], [22].In the sensor placement problem, it is often required to have the number of sensors to be at least the same or more than the latent rank of the system in order to be able to capture enough information for reconstruction [1].But with Kalman filter estimation, fewer sensors can be used to achieve the same performance given that the system is observable with these sensor measurements [22].Commonly, the Kalman filter sensor selection (KFSS) problem studies the objective based on a posteriori error covariance, which is a metric in Kalman filtering for how much the estimates deviate from the truth.The metric can be considered within an observation period [23], but it is more commonly taken to the limit at the infinite-time horizon when the full convergence of Kalman filtering is reached.Although optimization over the trace of the error covariance matrix, which represents the estimation MSE, does not have a constant-factor polynomial-time approximation [24]- [26], greedy methods are still near-optimal [27].
The diversity of mathematical methods highlighted above for optimal sensor placement typically focus on stationary, point sensors.However, in many applications, sensors can be mobile, in which case sensors are allowed to freely move in the measurement space while collecting measurements along the way.The problem concerning the design of trajectories or paths of sensors is called the sensor path planning problem.In the field of engineering and robotics, path planning problem has been long considered for the purposes of navigation as well as estimation in a dynamical environment [28]- [32].The task of tracking and estimating a flow field has often been tackled by constructing a simplified, restricted problem that focuses on a network of sensors with a simple formation for efficient paramerization and optimization [33]- [37].Different control laws for the path of the sensors are considered for different tasks, including a simple circular or elliptical control [33], gradient climbing control [35], control along level curves [36], or control based on smoothed particle hydrodynamics [38].Lynch et al. [39] propose a decentralized mobile network to collectively estimate environmental functions through communication networks, while the sensors move according to a gradient control law that maximizes information.Shriwastav et al. [40] built a trajectory by connecting a cost-efficient path among optimal sensor placement locations under POD based reconstruction.For many of these work, the emphasis is on modeling and control of the sensor positions.Sensor scheduling [41], [42] is a similar problem that concerns a schedule of densely placed sensors.Unlike the path planning problem, the sensors do not move in the scheduling problem, although it still can be formulated and solved as a special case of the path planning problem.
While many consider the sensor path in an infinite-time horizon, theoretical studies [43]- [45] show that the optimal infinite-time schedule is independent of the initial error covariance and can be approximated arbitrarily closely by a periodic schedule.This provides a mathematical foundation for studying problems that consider the planning of a periodic sensor trajectory for spatio-temporal estimation with Kalman filtering.Lan and Schwager [46], [47] approach the periodic path planning problem with a rapidly exploring random cycles (RRC) method that constructs and evaluates cycles found by randomly exploring the state space using a tree structure; Chen et al. [48] utilize deep reinforcement learning instead as a learnable deterministic method for finding cycles.The problem extends to multiple sensors that do not have a set network formation, each following its individual path.These works are most closely related to the problem considered here.They approach the combinatorial optimization with a randomized or active search method, first searching for possible cycles, then evaluating their costs.By assumption, the sensors move to a different location at each discrete time step based on the trajectory found in this way.
In this paper, we consider the use of mobile sensors to improve the performance of estimating spatiotemporal data with Kalman filtering, where we focus on planning a periodic sensor trajectory that optimizes estimation.We consider the condition number of the observability matrix of the model as a metric for the Kalman filter estimation.The study of observability is not new and has been discussed previously in different sensor problems [16], [34], [49]- [51].In particular, Manohar et al. [49] present a balanced model reduction for sensor and actuator selection through observability and controllability in a linear quadratic Gaussian (LQG) controller setting.We build on these ideas, developing an optimization for the path planning of mobile sensors with the objective of dynamic Kalman filter estimation.We identify three distinct timescales related to Kalman filter design and estimation with mobile sensors: (i) the timescale of the spatio-temporal dynamics, (ii) the velocity of the sensors, and (iii) the rate of sampling.We propose an approach for greedy selection based on the empirical observability matrix for path planning, and leverage low-rank representation of the system to promote efficient computation complexity.Figure 1 shows how our overall strategy leverages lowrank approximations in order to determine trajectories in the spatio-temporal fields of interest.Compared with previous works, our approach does not restrict the formation of the sensor network nor the shape of Figure 1: Overview of proposed approach to sensor path planning for dynamic estimation.The panels are divided into two main steps for estimating spatio-temporal data under a Kalman filter setting.The top panel shows the construction of a low-rank representation of the data as the prior model for Kalman filter through DMD.The DMD modes and eigenvalues make up a linear dynamical model in a reduced dimension and a projection back to the original diemsnion.The dimension of the observability matrix is also reduced by the low-rank representation for efficient computation.The bottom panel illustrates the greedy path finding algorithm that optimizes the observability matrix along the path and improves Kalman filter estimation performance.It leverages a greedy row selection on the projected full observability matrix.Conceptually, at each time step, based on the historical selection of sensor locations, the sensors are led to the next valid locations within a velocity constraint.the trajectory and builds the path by leveraging a low-rank system representation and greedy optimization.Our approach provides a scalable and efficient periodic path planning procedure for multi-sensor and highdimensional problems.We conduct a series of experiments on synthetic data, the Kuramoto-Sivashinsky system, and sea surface temperature data to show that mobile sensing improves Kalman filter performance in terms of better limiting estimation and faster convergence.

Problem Formulation and Background Methods
The mathematical formulation of the sensor selection problem considers a discrete-time linear system model where x t ∈ R n , A ∈ R n×n , and w t ∈ R n is the system disturbance following a Gaussian distribution with zero mean and a covariance matrix 0 ≺ Q ∈ R n×n .The measurements from k sensors are of the form where y t ∈ R k , C t ∈ R k×n , and v t ∈ R k is the measurement noise following a Gaussian distribution with zero mean and a covariance matrix 0 ≺ R ∈ R k×k .Directly measuring in the state space, we write the matrix C t as a selection matrix made up of standard unit vectors as columns.We further assume that the measurement noise are independent and identical across sensors with variance ρ, so the covariance matrix for v t is R = ρI.In a time-invariant system with a stationary sensing scenario at fixed locations, then C t = C.We denote a sensor trajectory σ = {σ 1 , σ 2 , ...} of k sensors.σ t ⊆ [n], |σ t | = k is a set containing sensor locations at time t, which determine the selection matrix C t = C(σ t ) responsible for collecting measurements along the trajectory.In general, σ can extend to an infinite-time horizon.Zhang et al. [43], [44] show that any infinite-time trajectory can be approximated by a periodic trajectory.Therefore, we focus on a periodic trajectory of fixed cycle rather than a trajectory over infinite-time.In practice, periodic trajectories also make sense since many systems contain some periodic or quasi-periodic characteristics.Furthermore, it is often favorable to plan a trajectory such that the sensor can return to a specified location periodically for maintenance and sensor recharging.Then, we write σ = {σ 1 , σ 2 , ..., σ l } to be a periodic trajectory of length l, so that σ l+1 = σ 1 , σ l+2 = σ 2 , and so on.
In Sec.2.1, we discuss the use of low-rank representation of the system for sparse sampling.We then introduce observability of the system in Sec.2.2 and relate it to Kalman filter estimation performance in Sec.2.3.Finally, we give attention to three key timescales in the Kalman filter model design in Sec.2.4.

Reduced-Order Model and Sparse Sampling
In order to promote efficient computation and better model representation for sparse sampling, we consider a reduced-order model (ROM).Specifically, we consider a system with a low-rank linear representation, where z t ∈ R m (m < n) is the internal low-rank dynamics state, Λ ∈ R m×m is the low-rank dynamical system, and Ψ ∈ R n×m is the linear projection basis.Measurements y t are collected in the original highdimensional state space.One can define a projection basis to be a universal basis for compressed sensing, or a tailored POD basis for a data-driven approach [1].However, such basis does not necessarily project to a proper low-rank dynamical system.To find a low-rank representation, suppose that the dynamics A is known, we can take a spectral decomposition of A and truncate the eigenvalues and eigenvectors to a low-rank representation.Alternatively, a data-driven approach is to find a close estimation of the model from the data by using dynamic mode decomposition (DMD) and its many variants [52]- [56] that can be useful in sparse sensing [57].DMD modes constitute the linear projection from high-dimensional data to the low-rank representation.The DMD eigenvalues form a diagonal dynamics matrix for the low-rank system.
ROMs are commonly utilized in the stationary sensor placement problem [1], [9], [12].Assuming no disturbance and noise, the measurements can be expressed as y t = CΨz t .Then, we can obtain x t through a simple linear reconstruction via the Moore-Penrose pseudoinverse, xt = Ψẑ t = Ψ(CΨ) † y t .It is clear that the reconstruction depends on the conditioning of matrix CΨ.Given a tailored basis, Q-DEIM [9] uses QR factorization with column pivoting (QRcp) to greedily find near optimal selections.At each step, QRcp selects a new pivot column with the largest norm and removes the orthogonal projections onto the pivot column from the remaining columns.Controlling the condition number by maximizing the matrix volume, QRcp enforces a diagonal dominance structure through column pivoting and expands the sub-matrix volume.In a more recent work, GappyPOD+E [12] extends the Q-DEIM method to an "oversampling" case where the sample/selection size is larger than the basis rank to improve stability.Based on the theory of random sampling in GappyPOD [11], it is a deterministic method that utilizes a lower bound for the smallest eigenvalue of the submatrix to continue sensor selection over model rank.

Observability
Observability is concerned with the possibility of finding the states of the system from the observations.A time-varying system of the form x t+1 = Ax t , y t = C t x t , or a pair (A, C t ), is observable at time t if the system state can be determined from the observations in [t, τ ] for some τ > t [58].The system is said to be observable if it is true for all time.Observability of a system is examined through the observability Gramian.In our discrete-time system, it is equivalent to study the observability matrix, The system is observable if and only if the observability matrix has full (column) rank.When all states are measured, C t = I, the full observability matrix is . In the reduced-order model representation, the projected full observability matrix is ) .
In a time invariant system where C t = C is fixed in time, it may need multiple sensors or a long period in time to achieve full rank of the observability matrix.For example, for a fully measured system, O is trivially full rank and the system states can be determined immediately at each time step.But for sparse sensing on a state space of large dimension, observability of the system is harder to achieve.By considering mobile sensing, it opens up possibilities to generate better observability with limited sensors.
A fully observable system is necessary for an accurate estimation using sparse measurements.In particular, it allows Kalman filter estimation to converge to steady-state values on an infinite-time horizon.We discuss in more detail the connection between observability and Kalman filter estimation in the following section.

Kalman Filter
We use a Kalman filtering (KF) estimator for spatiotemporal estimation on a linear model.Under the assumption of Gaussian noise, KF is known to be the best linear estimator for minimizing mean squared error [21].Kalman filter algorithms combine the information from the prior knowledge of the system and the observed measurements over time to find an optimal estimate of the system.Let Σ t denote the error covariance matrix at time t in the Kalman filter estimation.By definition, its trace is the expected squared estimation error at time t.The error covariance follows a recurrence relation In the time-invariant case, when t → ∞, the limiting error covariance satisfies Σ = AΣA * −AΣC * (CΣC * + R) −1 CΣA * + Q, which is known as the discrete algebraic Riccati equation (DARE).When the system is observable, the error covariance is guaranteed to converge to a limit or a limit cycle in a periodic schedule [43].
The limiting error covariance is a common metric to evaluate a KF model.However, finding this limit by solving the DARE is difficult and computationally expensive since it does not have a closed-form solution.Therefore, knowing that observability is a necessary condition for KF estimation performance, we further show how the conditioning of the observability matrix drives the KF estimation.We first relate the limiting expected squared error from the KF estimation to the conditioning of C in a time-invariant model with full-rank measurements.We then show that any model with sparse measurements or periodic trajectory can be reformulated at a larger time step to a time-invariant representation with full-rank measurements.And the reformulated selection matrix is the same as the observability matrix in the original form.
The expected squared error is represented as the trace of the error covariance matrix, whose limit is the solution of a DARE.Since the DARE does not have a closed-form solution, we consider an upper and a lower bound for the trace of the solution.While various works have derived different bounds on the DARE solution [59], [60], we utilize the following results for our analysis: We then have bounds: , where We can easily show that the lower bound is monotonically decreasing with λ 1 (C * R −1 C), and the upper bound is monotonically decreasing with The condition is usually satisfied as well for a stable system in general when the disturbance covariance does not have a heavily dominant eigenvalue.
Since we consider the model with independent and identical measurement noise, R = rI, so , where σ i (C) is the i-th largest singular value of C. Therefore, in a time-invariant model where C is full rank, we can minimize the condition number κ(C) = σ1(C) σn(C) in order to achieve lower squared error.
However, in most scenarios the system model is more complicated.When using limited sensors, the measurements C will not be full rank.In the mobile sensor with periodic trajectory scenario where C t depends on time, the system is not time-invariant.We can show a reformulation of these models to a timeinvariant representation in which C is full rank.Then, the above result applies to these models as well.Consider the general model x t+1 = Ax t + w t , y t = C t x t + v t .Let k = nt be a larger time step where n is the dimension of the state space or multiples of the sensor trajectory period.Then we can follow [63] and rewrite the system as follows: In the reformation, Ĉ is time-invariant, and it is exactly the observability matrix O in the original form.If the system is observable, O has full rank, and so does Ĉ.By representing a time-variant system of mobile sensors in a time-invariant form, we can draw the same conclusion as the time-invariant system that the condition number of the observability matrix bounds the limiting error covariance matrix of KF estimation.Thus, lowering the condition number of the observability matrix leads to better KF estimation performance.

Kalman Filter Design Factors
In the mobile sensor scenario, besides planning the trajectory of the sensors, we should also consider in the model design the following key factors: the system Nyquist rate, discrete sampling rate, and sensor speed.These three timescales relate to the conditioning of the observability matrix of the system and the performance of Kalman filter estimation.Although not a definite guide, the following provides useful heuristics for estimation performance based on these timescales.
The Nyquist rate represents the internal time scale of the continuous-time dynamics.It is defined to be twice the highest frequency of the spatiotemporal dynamics.The discretization of the continuous-time system is considered good if it samples faster than the Nyquist rate.We believe the same applies to mobile sensing with Kalman filter estimation.At least one measurement should be collected within Nyquist rate to capture the highest frequency feature of the system at the most relevant location.
The sampling rate refers to the rate at which the measurements are collected.It also represents the time step of the discrete model.Faster sampling rate above Nyquist adds more measurements in a fixed time interval.In the stationary sensor setting, this improves stability of the estimation.With mobile sensors, faster sampling rate further adds more information since the measurement locations change.This leads to better system observability and KF estimation until the statistical independence of the measurements no longer holds.
The sensor speed determines the maximum region a sensor can measure in a fixed time interval.A faster sensor can reach and observe at a farther location in the state space to achieve better observability.More importantly, when the data contains local structures, it is essential to plan the sensor trajectories to capture those structures.Faster sensors can achieve this purpose when local structures are well separated in the state space, without the need to assign additional sensors.
The effect of these timescales will be further discussed in the numerical experiments in Section 4.

Computing Mobile Sensor Trajectories
With the problem formulation 3, and the discussion in Section 2.3, we consider the objective to minimize the condition number of the observability matrix under the schedule σ: The observability matrix with respect to trajectory σ of length l is written as where O Ψ is the projected observability matrix with full measurements, and C σ is a block-diagonal selection matrix determined by σ.Problem 4 becomes a submatrix selection problem minimizing the condition number.In the special case when the length of the periodic trajectory is 1, the objective becomes max σ:|σ|=k κ(C σ Ψ), which is identical to that of a stationary sensor placement problem under the DMD basis.Solving such a problem is in general NP-hard, but just as in the stationary sensor placement problem, we can leverage greedy algorithms and utilize the same idea as QRcp/Q-DEIM for under-sampling and GappyPOD+E or over-sampling to solve it approximately.
We introduce our greedy time-forwarding algorithm in Section 3.1 and illustrate on a synthetic example of sparse linear dynamics on a torus in Section 3.2 before presenting the main results in Section 4.

Algorithm
The projected full observability matrix O Ψ is by definition segmented into blocks, so for the purpose of efficient computation, we propose a greedy algorithm that finds sensor locations σ 1 , σ 2 , ..., σ l by sequentially focusing on each block Ψ, ΨΛ, ..., ΨΛ l−1 .

The algorithms are as follows:
Algorithm 1: Greedy Time-forwarding Observability-based Path Planning Algorithm Input: low-dimension dynamics matrix Λ, projection basis Ψ, number of sensors k, trajectory period l.Output: trajectory σ.
:] 2 ; end Sort r in descending order and select s be the first valid ordered index in S; In the selection step, we want to find a row in X from the candidate set S to append to the current observability matrix in order to minimize κ(O σ ).We use the same selecting rules as in QRcp and Gappy-POD+E.The candidate set S is critical when sensor movement constraints are present.When the sensor is unrestricted to move in time, we can simply set S = [n]\σ i .However, in practice, the sensors have a limit on their speed so there is a restricted region in which the sensors can move between time steps.Additionally, the state space can have special multiply-connected topological structure such that not all locations are accessible from every other.Future work will incorporate the background flow field in this estimated restriction region, although this is neglected for simplicity in the present work.
Under a sensor speed constraint, we only consider the locations where • a sensor can move to within a time step from its current location; • it can go back to its initial location at the end of the period to form a cycle.
These requirements guide the selection of the candidate set S in the algorithm.When the topology of the state space is regularly shaped, a simple Euclidean distance can be used; while it is irregular with obstructions or complex network structures, we can resolve to other types of distance functions.

Illustrative Example: Sparse Linear Dynamics on a Torus
To show the effectiveness of mobile sensors, we demonstrate the algorithm with a random simulation of sparse linear dynamical system on a torus [64].We design the system to contain two types of structures: the 2D discrete inverse Fourier transform function and the Gaussian basis function.A Fourier mode is a global feature present across the state space, while a Gaussian mode is a local feature that only concentrates in a small neighbor around a center.On a 128x128 spatial grid, we build the sparse system with 2 conjugate Fourier modes and 3 conjugate Gaussian modes by generating randomly their oscillation frequencies and damping rates.This is a system of size n = 128 2 = 16384 with a low-dimensional linear representation of rank m = 10, where the projection basis Ψ contains the modes, and the low-dimension linear dynamics matrix Λ is diagonal with the oscillation and damping information.We generate the data by adding system disturbances and measurement noise.Since all parameters in the model are known in the synthetic example, we use the trace of the error covariance matrix as an accurate representation of the expected squared error to evaluate the estimation.
First, we estimate the system with sensors at fixed locations selected by applying QRcp on the basis Ψ, a common sensor placement strategy.We see from Figure 3 there is significant performance improvement as we increase the number of fixed sensors up to 3. At least 3 sensors are needed to obtain a good estimation of the system so that they can be placed to observe the local regions of the Gaussian modes.
We then show that equivalent performance can be achieve using only one mobile sensor with the same sampling rate and fast enough speed.We choose a trajectory period such that the cycle is complete within the Nyquist rate of the system.When the sensor is slow, there is no significant improvement since the sensor cannot move to other local features within a cycle.But, with fast enough speed, our algorithm is able to direct the sensor to reach the localization of all three Gaussian modes and make better estimation (Figure 4).Under the same sampling rate, three sensors collect three times as many measurements as only one sensor within any time interval.This fundamental differences in measurement size due to number of sensors contributes to the difference between three stationary sensors and one mobile sensors.We can narrow this performance gap by increasing the sampling rate.At a sampling rate of 0.001, the difference in estimation error is minimal.
Through this synthetic experiment we see that a mobile sensor can indeed improve Kalman filter estimation, and the trajectory planned by our greedy method is effective to pinpoint local structures and achieve good observability.

Numerical Experiments
In practice, it is often the case that for spatiotemporal data and systems, the underlying low-rank dynamic model, the disturbance, and the noise are not known.In this case, we would fit an estimated model representation from data via DMD, and approximate disturbance and noise covariances either from data or through hyper-parameter tuning.Here, we look at two examples: (i) the Kuramoto Sivashinsky (KS) system, and  (ii) the Sea Surface Temperature (SST) dataset from NOAA [65], [66].We study the performance when we use a DMD approximated model for Kalman filter estimation and sensor path planning by our greedy algorithm.In both examples, we fit a linear DMD model with a chosen low rank.The dynamics matrix Λ is diagonal with DMD eigenvalues and the basis Ψ consists of the DMD modes.We further add a white measurement noise with variance R = I to the data, and we set the system disturbance to be Q = qI where the uniform variance q is a hyperparameter tuned by experiment.
These two examples are representative in different aspects.The KS system is known for its chaotic behavior.Therefore, a linear representation of the system is extremely difficult.Additionally, the modes from the linear approximated model are mostly local since linear correlation between locations is small, so that full observability is hard to achieve with few fixed sensors.We show in Section 4.1 that mobile sensors can perform particularly well comparing to fixed sensors by reaching more locations and capturing more local structures.
The SST dataset from NOAA contains weekly mean optimum interpolated sea surface temperature measurements from global satellite data.The dataset can be well approximated by a linear model and most modes in the approximated system are global so that observability is easily achieved with even just one stationary sensor.We show then in Section 4.2 mobile sensors further accelerate the convergence of error.

Kuramoto Sivashinsky System
The KS system is given by the equation u t + uu x + u xx + u xxxx = 0. We consider the numerical solution of the system on a spatial grid of size 2048 over x ∈ [0, 2π].The initial condition is randomly generated over a standard normal distribution.With a random initial condition, we numerically solve the KS equation and collect data on the time interval t ∈ [0, 10] with a time step of dt = 10 −4 .We first perform singular value decomposition (SVD) to find a low rank representation of the data.The first 100 singular values capture 99.99% of the energy, so we estimate a low-dimensional linear representation of the system by fitting a standard dynamic mode decomposition (DMD) model [52] with an SVD rank of 100.
Because of the chaotic nature of the system, accurate estimation is not possible with a limited number of 10 sparse fixed sensors.Indeed, we need 100 fixed sensors, equivalent to the full rank of our approximated linear system, to effectively estimate the system (Figure 5).Additionally, we see that there is no significant improvement in performance with increasing sampling rate using fixed sensors since more frequent measurements at the same locations add little information of the unobserved states.
On the other hand, mobile sensors can move to measure different locations and gain more information of the entire state space.With fast enough sensor speed, 10 moving sensors can achieve a significantly improved estimation comparing to 10 fixed sensors by increasing the sampling rate (Figure 5).The improvement in performance is limited by the sensor speed.We set the minimum speed in this example to be v min = 2π 2048 * 10 4 so that the sensor is able to move to its left and right neighbor at a discrete sampling time step of 10 −4 .With higher sensor speed, sensors can make observations over a wider spatial range, thus giving better estimation.As v → ∞, the performance of 10 moving sensors approaches that of 100 fixed sensors with fast enough sampling rate.
Due to the greedy nature of our algorithm, it selects based on the immediate reward at the next time step and cannot look ahead.When the sampling is sufficiently fast, the greedy algorithm makes a decision based on the closest neighbors of the current location.Such a decision is not informative, and the trajectory planned fails to have a good performance.One way to reduce the greediness of the algorithm is to perform a multiscale path completion.We start by finding a trajectory at a slower sampling rate.Then, we gradually decrease the time step and apply the same path planning algorithm, except using the previously found trajectory as guidance and filling in the gap to construct a more complete trajectory at the faster sampling rate.We apply this multiscale expansion procedure on the KS example, initiated at the sampling rate with the smallest error, and expand to faster sampling rate based on that path.We see the performance is no longer worse with a fast sampling rate in the KS experiments, but it flattens and reaches a limit determined by the sensor speed (Figure 6).

Sea Surface Temperature
The Sea Surface Temperature (SST) dataset contains weekly collection from satellite data of sea surface temperature measurements on the 1 degree latitude by 1 degree longitude (180 by 360) global grid from 1990 to the present.We fit a standard DMD model with a low dimension of 10.
First, we see that Kalman filter estimation using one stationary sensor is comparable to that of ten sensors (Figure 7).This verifies that the approximated linear model contains mostly global features that can be observed well with very few sensors.However, with a bad initial estimation, KF with 10 fixed sensors converges to low error very quickly (below 0.1 within one year), while it takes 1 sensor almost 28 years to reach a comparable error.
Next, we consider estimating the SST data with one mobile sensor.Since the limiting estimation is already good with one stationary sensor, we show the improvements in the KF convergence speed for mobile sensors.The DMD eigenvalues show the max frequency of the system to be around a half-year, so we choose the period of mobile sensor trajectory to be about a quarter-year (14 weeks).Since the globe has continental land as obstructions, we need to ensure the planned trajectory does not cross any land as the sensor moves in water.We build a connectivity graph and adjacency matrix for the candidate selection step in our algorithm instead of a simple Euclidean distance function.
Figure 7 shows the results of one mobile sensor with different sensor speed limits.Mobile sensor estimation indeed produces much faster convergence compared to a stationary sensor.As the speed limit increases, the sensor can move to farther locations with better observability, further improving the convergence of estimation.Figure 8 shows the paths of the sensor.When v is small, the initial location plays an important role since the trajectory does not move far from it.The first location picked by the algorithm is close to Alaska, so the first two trajectories with low speed center around the North Pacific and the Arctic Ocean.As v increases, the sensor explores the equator and south hemisphere regions, especially the El Nino regions around the equatorial Pacific, which is an important local feature.Figure 9 shows the planned trajectory using 2 mobile sensors.
KF estimation convergence also matters when the underlying dynamics is nonstationary and changes over time.If the estimation does not reach a meaningful error in time, the shifting dynamics will further slow down the convergence and increase the limiting error.To show this, we instead fit a DMD model using only the first half of the SST data and use it as the approximated linear model for KF estimation.In this case, the fitted linear model is representative and relevant only in the first training half, and does not reflect any possible changes in the data dynamics afterwards.Then, one stationary sensor performs significantly worse due to slow convergence (Figure 10).On the other hand, the error from one mobile sensor converges fast enough within the training period so that in the second half the error is still relatively low.Therefore, fast convergence with mobile sensors ensures a fast adjustment in estimation when the dynamics change in time.

Conclusion and Future Work
In this work, we developed a mathematical strategy for planning a periodic trajectory for limited mobile sensors to estimate a spatiotemporal system using Kalman filter estimation.We examine the system observability as a metric that influences the estimation performance in terms of the limiting squared error as well as the convergence rate.We consider an objective to minimize the condition number of the discrete observability matrix along the trajectory and formulate it as a submatrix selection problem.We then propose a time-forwarding greedy algorithm that selects sensor locations along the trajectory using the same rules as QRcp and GappyPOD+E from a carefully chosen candidate subset.The experiments show that the method is able to plan a trajectory that locates the local features and improves the estimation performance.In these experiments, we explore Kalman filter design factors and their impact on estimation as they relate to the three important timescales: the Nyquist rate of the underlying dynamics, the rate of sampling, and the velocity of the sensors.We find that mobile sensors are especially beneficial for a complex, non-linear system to capture local features in an approximated linear model, without deploying a large amount of sensors.We also see an improvement in estimation convergence rate using mobile sensors, which more rapidly reaches an accurate estimation.
In future work, a weighted cost function can be added to the objective to better incorporate different costs to the path planning.Sensor speed can be formulated as a cost instead of a hard constraint imposed in the selection process.Furthermore, energy consumption caused by sensor movement can also be included to plan a trajectory that is also more energy efficient.For example, as in the flow field applications, we can consider the flow field information and the energy cost associated with it as the sensor moves with or against the flow.Incorporating the background flow velocity in the set of possible next locations is an important future extension of this work.We can also refer to many different sensor control laws as cost constraints for incorporating other tasks such as simultaneous structure tracking.
For multi-sensor planning, it will be interesting to consider different asynchronous periodic trajectory for each individual sensor instead of all having the same period.This will be particularly useful for multiscale systems so that each sensor can be responsible for estimating features of different timescales.As shown in the numerical experiments, the performance of Kalman filter estimation fundamentally depends on the accurate modeling of the linear system and the correct choices of hyper-parameters.Better data-driven linear system identification can be explored.An alternating model fitting and estimation approach can be explored to update both the model and the sensor trajectory continuously to achieve even better performance.

Figure 2 :
Figure 2: A snapshot of the random system in 2D (left) and on a 3D torus (right).

Figure 3 :
Figure 3: Expected squared error of the KF estimation in time, (a) stationary sensor placement by number of sensors; (b) one mobile sensor by velocity constraints.

Figure 4 :
Figure 4: Sensor locations along trajectory with speed of 5 (left) and 37 (right) units per time step of 0.01.

Figure 5 :
Figure 5: (a) True spatiotemporal dynamics of the KS system in T ∈ [9, 10]; (b) Bode plot of estimation error against sampling rate; (c) Estimated x-t plot by 10 mobile sensors with sampling rate dt = 0.001, 0.005, 0.01, 0.02 (corresponding with the dashed vertial lines on the bode plot).

Figure 6 :
Figure 6: Estimation error against sampling rate with multiscale expansion plotted in full transparency connected by line.

Figure 8 :
Figure 8: Planned sensor trajectory with a cycle period of 14 weeks, where the movement speed is limited to 5, 15, 25, 50 spatial units (1 degree of latitude or longitude).Zoomed in map on the right.

Figure 9 :Figure 10 :
Figure 9: Planned sensor trajectory for 2 sensors with a cycle period of 14 weeks, with sensor speed limit at 15 spatial units (1 degree of latitude or longitude).
update candidate rows S in X; Calculate scores and select the next valid row s ∈ S by Algorithm 2; Assign the closest and previously unassigned sensor s ∈ σ i−1 to move to s, and put s to corresponding position in σ i ; Append the selected row X[s, :] to O σ ; end Add σ i to σ; Input: target matrix X, current observability matrix O σ , candidates S. Output: row index s.p ← row size of O σ ; if p < m then // under-sampling, use QRcp rule [