Decentralized state estimation: An approach using pseudomeasurements and preintegration

Abstract

This paper addresses the problem of decentralized, collaborative state estimation in robotic teams. In particular, this paper considers problems where individual robots estimate similar physical quantities, such as each other’s position relative to themselves. The use of pseudomeasurements is introduced as a means of modeling such relationships between robots’ state estimates and is shown to be a tractable way to approach the decentralized state estimation problem. Moreover, this formulation easily leads to a general-purpose observability test that simultaneously accounts for measurements that robots collect from their own sensors, as well as the communication structure within the team. Finally, input preintegration is proposed as a communication-efficient way of sharing odometry information between robots, and the entire theory is appropriate for both vector-space and Lie-group state definitions. To overcome the need for communicating preintegrated covariance information, a deep autoencoder is proposed that reconstructs the covariance information from the inputs, hence further reducing the communication requirements. The proposed framework is evaluated on three different simulated problems, and one experiment involving three quadcopters.

1. Introduction

Decentralized state estimation is a fundamental requirement for real-world multi-robot deployments. Whether the task is collaborative mapping, relative localization, or collaborative dead-reckoning, the multi-robot estimation problem seeks to estimate the state of each robot given all the measurements that each robot obtains locally. This problem is made difficult by the fact that not all robots can communicate with each other and, furthermore, that high-frequency sensor measurements would require substantial communication bandwidth to simply share across the team. A robot might even have insufficient sensors to observe their own state, and hence is dependent on its neighbor’s sensors to have a stable estimate. Hypothetically, an estimator that could somehow collect all these sensor measurements on each robot, and fuse them all to jointly estimate the states of every robot in one large system, would have the lowest possible estimation error variance. This is called the centralized estimator, but is often infeasible to implement in practice.

A common approach is for robots to share their current state and associated covariance rather than of a history of measurement values (Julier and Uhlmann (1997); Julier (2001); Carrillo-Arce et al. (2013)). This approach has the benefit of simplicity, low communication cost, and fixed message size, but suffers from a well-known issue of not being able to compute cross-correlations between the robots’ state estimates (Shalaby et al. (2021b)). Furthermore, in certain problems, robots may be estimating the same physical quantities. As an example, consider two robots estimating each other’s position, in addition to their own positions, as shown in Figure 1 (left). Their state vectors are both robots’ positions, and therefore, both seek to estimate the same physical quantities, a situation referred to here as full state overlap. When robot states have similar, if not identical state definitions, it is straightforward to compute the error between their state estimates using simple subtraction. However, for more complicated problems, especially those with state definitions belonging to arbitrary Lie groups, a generalized measure of error between different robots’ state estimates must be introduced.

Figure 1.

Three examples of decentralized estimation problems within the scope of this paper. Left: A toy problem with 1D robots, each estimating both of their positions. Middle: A problem with an incomplete communication graph. Robots observe landmarks, have range measurements to each other, and estimate their own and neighbor absolute poses. Right: A more complicated experimentally tested problem, where robots equipped with ultra-wideband radios estimate both their own absolute pose and relative poses of neighbors, in addition to IMU biases.

As a concrete example, which is featured experimentally in this paper, consider the case where quadcopters each possess GPS sensors, IMUs, and inter-robot range measurements using ultra-wideband radio. Suppose one quadcopter loses GPS functionality due to traveling under a bridge or a sensor fault, thus losing absolute positioning information. The challenge is to design an algorithm where this faulty robot maintains accurate absolute positioning estimation by appropriately sharing information with its neighbors.

1.1. Contributions

This paper has three main contributions:

1. a framework for decentralized state estimation that uses pseudomeasurements to allow for generic nonlinear relationships between robot states;

2. a preintegration-based method for constant-time, constant-memory, and constant-communication odometry sharing;

3. a theory compatible with Lie-group state definitions, including the familiar vector space.

The pseudomeasurements are shown to be a tractable and effective way to model any generic nonlinear relationship between robot state definitions, including full or partial state overlap as special cases. For example, a nonlinear state relationship is present when robots estimate each other’s poses in their own body frames. A pseudomeasurement is introduced for each edge in the communication graph, and the proposed framework also naturally leads to an observability test that takes into account both the local measurements obtained by each robot and the communication structure between them. Usage of pseudomeasurements requires robots to communicate their states and corresponding covariances. Furthermore, the common states between robots must be at the same time step, which potentially requires odometry information to also be shared, so that states can be propagated forward to a common time.

The proposed use of preintegration allows sharing odometry information over an arbitrary duration of time, in a lossless matter. The naive alternative is for robots to share a history of odometry measurements since the last time they communicated, which has processing, memory, and communication requirements that grow linearly with the time interval between communications. Preintegration provides a constant-time, constant-memory, and constant-communication alternative that is algebraically identical to simply sharing the input measurements themselves. This makes preintegration a natural choice for multi-robot estimation problems. Moreover, preintegration preserves statistical independence assumptions that typical Kalman filtering prediction steps rely on. Preintegration is best known from the visual-inertial odometry literature (Lupton and Sukkarieh (2012); Forster et al. (2017)), where the same concept, adapted for relative pose estimation is introduced by Shalaby et al. (2023). This paper generalizes the multi-robot-preintegration concept to other common process models in robotics, presents a solution for simultaneous input bias estimation, and also further proposes a deep autoencoder to compress the associated covariance information.

Finally, the proposed solution is general to any state definition, process model, and measurement model subject to typical Gaussian noise assumptions. The complexity of the proposed estimation algorithm is identical to a standard extended Kalman filter, and the communicated messages are lightweight and of fixed length. In the experimental test demonstrated in this paper, each robot transmits information at a rate of only 53 kB/s.

This paper does not focus on the treatment of cross-correlations, and hence employs the simple, well-known covariance intersection (CI) (Julier and Uhlmann (1997); Julier (2001)) method. This allows for an arbitrary communication graph within the robot team, while remaining lightweight and avoiding cumbersome bookkeeping. The main drawback is that CI is proven to yield suboptimal estimation. However, in practice, the performance can remain adequate, and sometimes very comparable to centralized estimation (Shalaby et al. (2021b); Julier (2001)), as shown here from simulated and experimental results.

The remainder of this paper is as follows. Related work is discussed in Section 2 and mathematical preliminaries and notation are shown in Section 3. The paper then starts with a simplified “toy” problem showcasing the proposed method in Section 4, and the theory is generalized in Section 5. Preintegration sharing is presented in Section 6. Finally, Section 7 contains an application to a ground robot simulation and Section 8 applies the method to a more complicated experimental quadcopter problem.

2. Related work

There are many sub-problems associated with the overall decentralized state estimation problem. Even if communication links between robots are assumed to be lossless and have infinite bandwidth, there is still the issue of propagating information over general, incomplete communication graphs. For two robots, it is straightforward to compute centralized-equivalent estimators on each robot as done by Grime and Durrant-Whyte (1994), which use information filters to accumulate the information from a series of measurements, and then communicate this quantity. Grime and Durrant-Whyte (1994) also derived centralized-equivalent solutions for fully connected graphs, as well as tree-shaped graphs. However, they show that for generic graphs, it is impossible to obtain a centralized-equivalent estimate with only neighboring knowledge, and that more knowledge of the graph topology is required.

Leung et al. (2010) present a general centralized-equivalent algorithm for arbitrary time-varying graphs, and is formulated over distributions directly, hence allowing inference using any algorithm such as an extended or sigma-point Kalman filter. Robots must still share raw measurements with each other, therefore requiring substantial bookkeeping. The approach is extended to the SLAM problem by Leung et al. (2012). Roumeliotis and Bekey (2002) decompose the centralized Kalman filter equations using a singular value decomposition to generate independent equations that each robot can compute. Provided that robots have broadcasting ability, and obtain direct pose measurements of their neighbors, a centralized-equivalent solution can be obtained. The consensus Kalman filter (Olfati-Saber and Murray (2004); Olfati-Saber (2005)) aims to asymptotically send the state of n arbitrary nodes to a common value, which is effectively a problem with full state overlap. Battistelli et al. (2015) proposes alternate consensus approaches such as “consensus on information” or “consensus on measurements.” However, the problem does not consider the fact that robots collect their own, separate odometry measurements that are not necessarily available to neighboring robots.

As previously mentioned, one of the simplest solutions to the decentralized estimation problem is to use covariance intersection (Julier and Uhlmann (1997); Julier (2001)). CI conservatively assumes maximum correlation between robot estimates. Although the performance is theoretically suboptimal, the implementation is extremely simple, and imposes no constraints whatsoever on the communication frequency or graph topology. Carrillo-Arce et al. (2013) apply CI to a collaborative localization problem, where each robot estimates their own absolute state given direct relative pose measurements to other robots. Meanwhile, Arambel et al. (2001) present a decentralized state estimation algorithm for multiple spacecraft, where each spacecraft estimates the full state of all vehicles and then utilizes CI to fuse the neighbors’ full state. Recently, split-CI has been introduced to separate states into groups of correlated and independent substates (Li and Nashashibi (2013)), while Li and Yang (2021) exploit CI for the fusion of poses on Lie groups.

When employing CI, a user-defined weighting parameter has to be chosen, which affects the level of inflation of the block-diagonal components of the covariance matrix. Zhu and Kia (2019) formulate an optimization problem where the logarithm of the determinant of the posterior covariance matrix is minimized as a function of the CI weighting parameter, alongside an alternative linear-matrix-inequality approach that estimates the most conservative posterior covariance matrix. Meanwhile, Luft et al. (2018) use an EKF-like filter for decentralized estimation where cross-correlations are also explicitly tracked for both the prediction step and the fusion of local measurements. When relative measurements are encountered, an improved approximation to the joint covariance matrix is developed, which outperforms CI. The approach of Luft et al. (2018) assumes that process model inputs between robots are uncorrelated, which is not applicable in some of the problems in this paper. The work by Jung et al. (2020) builds off of Luft et al. (2018) to solve a full 3D collaborative state estimation problem where each robot has a camera and an IMU.

Another approach using scattering theory has recently been presented for two robots (Allak et al. (2019; 2022)), with the objective of reducing the communication cost associated with high-rate sensor measurements. Also making reference to the IMU preintegration technique (Lupton and Sukkarieh (2012); Forster et al. (2017)), covariance pre-computations are derived by Allak et al. (2019) and later extended to also include the mean (Allak et al. (2022)). It is shown that by sharing pre-computed matrices with twice the size as the state vector, a centralized-equivalent state estimate can be directly obtained with no measurement reprocessing. However, the generalization to more robots does not seem straightforward.

A variety of optimization-based approaches can be seen in the literature, especially when applied to multi-robot simultaneous localization and mapping (SLAM). Tian et al. (2022) have released Kimera-Multi, which uses a distributed pose-graph optimization algorithm to perform metric-semantic SLAM. Lajoie and Beltrame (2023) propose Swarm-SLAM, which performs multi-robot SLAM with an emphasis on using sparsity to minimize the number of data exchanges. However, these distributed SLAM methods are appropriate for situations where each robot has sufficiently rich sensor information via cameras or LIDARs and can perform individual SLAM in the first place. The method of this paper does not impose such a requirement.

3. Preliminaries

This paper will address problems where an individual robot’s process model f(⋅) and measurement model g(⋅) are modeled in the standard form of

$\begin{array}{ll}{\mathcal{X}}_{i_k}& =\mathbf{f}\left({\mathcal{X}}_{i_{k-1}},{\mathbf{u}}_{i_{k-1}},{\mathbf{w}}_{i_{k-1}}\right),\\ {\mathbf{y}}_{i_k}& =\mathbf{g}\left({\mathcal{X}}_{i_k}\right)+{\mathbf{v}}_{i_k},\end{array}\hspace{1em}\begin{array}{ll}{\mathbf{w}}_{i_{k-1}}& \sim \mathcal{N}\left(\mathbf{0},{\mathbf{Q}}_{i_{k-1}}\right),\\ {\mathbf{v}}_{i_k}& \sim \mathcal{N}\left(\mathbf{0},{\mathbf{R}}_{i_k}\right),\end{array}$

(1)

for Robot i, where ${\mathbf{u}}_{i_k}\in {\mathbb{R}}^{n_u}$ is the process input at time step k, ${\mathbf{y}}_{i_k}\in {\mathbb{R}}^{n_y}$ are the measurements, and ${\mathcal{X}}_{i_k}\in G$ denotes the robot state belonging to any Lie group G. As a notational convenience, the shorthand ${\mathcal{X}}_{i:j}=\left\{{\mathcal{X}}_i\dots {\mathcal{X}}_j\right\}$ will refer to a collection of arbitrary objects with indices in the range [i, j].

3.1. Lie groups

A Lie group G is a smooth manifold whose elements, given a group operation ◦: G × G → G, satisfy the group axioms (Solà et al. (2018)). The application of this operation to two arbitrary group elements $\mathcal{X},\mathcal{Y}\in G$ is written as $\mathcal{X}\mathcal{\circ }\mathcal{Y}\in G$ . For any G, there exists an associated Lie algebra $\mathfrak{g}$ , a vector space identifiable with elements of ${\mathbb{R}}^m$ , where m is referred to as the degrees of freedom of G. Lie algebra elements are related to group elements through the exponential and logarithmic maps, denoted $\exp :\mathfrak{g}\to G$ and $\log :G\to \mathfrak{g}$ . The “vee” and “wedge” operators are denoted ${(\cdot )}^{\vee }:\mathfrak{g}\to {\mathbb{R}}^m$ and ${(\cdot )}^{\wedge }:{\mathbb{R}}^m\to \mathfrak{g}$ , which can be used to associate Lie algebra elements with vectors. Composing these operators, group elements can be associated with vectors using

$\mathcal{X}=\exp \left({\mathit{\xi }}^{\wedge }\right)\triangleq \mathrm{Exp}\left(\mathit{\xi }\right),\mathit{\xi }=\log {\left(\mathcal{X}\right)}^{\vee }\triangleq \mathrm{Log}\left(\mathcal{X}\right),$

where $\mathcal{X}\in G,\mathit{\xi }\in {\mathbb{R}}^m$ , and the shorthand notation $\mathrm{Exp}:{\mathbb{R}}^m\to G$ and $\mathrm{Log}:G\to {\mathbb{R}}^m$ has been defined. Following Solà et al. (2018), the adjoint matrix representation of an element $\mathcal{X}\in G$ is denoted $\mathbf{Ad}:G\to {\mathbb{R}}^{m\times m}$ and defined such that

$\mathbf{Ad}\left(\mathcal{X}\right)\mathit{\xi }={\left(\mathcal{X}{\mathit{\xi }}^{\wedge }{\mathcal{X}}^{-1}\right)}^{\vee }.$

The most common Lie groups appearing in robotics are SO(n), representing rotations in n-dimensional space, SE(n), representing poses, and SE₂(3) representing “extended” poses that also contain velocity information. In these cases, the elements $\mathcal{X}$ are invertible matrices and the group operation ◦ is regular matrix multiplication.

3.1.1. ⊕ and ⊖ operators

Estimation theory for vector-space states and Lie groups can be elegantly aggregated into a single mathematical treatment by defining generalized “addition” $\oplus :G\times {\mathbb{R}}^m\to G$ and “subtraction” $\ominus :G\times G\to {\mathbb{R}}^m$ operators, whose precise definitions will depend on the problem at hand. For example, possible implementations include

$$\begin{array}{ll}\mathcal{X}\oplus \delta \mathbf{x}=\mathcal{X}\circ \mathrm{Exp}\left(\delta \mathbf{x}\right)& \left(\text{Lie} \text{group} \text{right}\right),\\ \mathcal{X}\oplus \delta \mathbf{x}=\mathrm{Exp}\left(\delta \mathbf{x}\right)\mathcal{\circ }\mathcal{X}& \left(\text{Lie} \text{group} \text{left}\right),\\ \mathbf{x}\oplus \delta \mathbf{x}=\mathbf{x}+\delta \mathbf{x}& \left(\text{vector} \text{space}\right),\end{array}$$

for addition and, correspondingly,

$$\begin{array}{ll}\mathcal{X}\ominus \mathcal{Y}=\mathrm{Log}\left({\mathcal{Y}}^{-1}\mathcal{\circ }\mathcal{X}\right)& \left(\text{Lie} \text{group} \text{right}\right),\\ \mathcal{X}\ominus \mathcal{Y}=\mathrm{Log}\left(\mathcal{X}{\mathcal{\circ }\mathcal{Y}}^{-1}\right)& \left(\text{Lie} \text{group} \text{left}\right),\\ \mathbf{x}\ominus \mathbf{y}=\mathbf{x}-\mathbf{y}& \left(\text{vector} \text{space}\right),\end{array}$$

for subtraction. This abstraction is natural since a vector space technically qualifies as a Lie group with regular addition + as the group operation.

3.1.2. Gaussian distributions on Lie groups

As an example use of this abstraction, consider defining a normally distributed Lie group element with mean $\overline{\mathcal{X}}$ and covariance Σ, as done by Barfoot and Furgale (2014), with

$\mathcal{X}=\overline{\mathcal{X}}\circ \mathrm{Exp}\left(\delta \mathbf{x}\right),\delta \mathbf{x}\sim \mathcal{N}\left(\mathbf{0},\mathbf{\Sigma }\right),$

when using a right parameterization, or a similar definition for left parameterizations. This can alternatively be written in an abstract way, applicable to any group or vector space, with

$\mathcal{X}=\overline{\mathcal{X}}\oplus \delta \mathbf{x},\delta \mathbf{x}\sim \mathcal{N}\left(\mathbf{0},\mathbf{\Sigma }\right).$

(2)

Moreover, given that $\delta \mathbf{x}=\mathcal{X}\ominus \overline{\mathcal{X}}$ , it follows from (2) that

$p\left(\mathcal{X}\right)=\eta \hspace{0.17em}\exp \left(-\frac{1}{2}{\left(\mathcal{X}\ominus \overline{\mathcal{X}}\right)}^{\mathrm{T}}{\mathbf{\Sigma }}^{-1}\left(\mathcal{X}\ominus \overline{\mathcal{X}}\right)\right)\triangleq {\mathcal{N}}_L\left(\overline{\mathcal{X}},\mathbf{\Sigma }\right),$

where the reader should note the definition of the generalized Gaussian ${\mathcal{N}}_L\left(\overline{\mathcal{X}},\mathbf{\Sigma }\right)$ (Bourmaud et al. (2016)).

3.1.3. Derivatives on Lie groups

Again following Solà et al. (2018), the Jacobian of a function f: G → G, taken with respect to $\mathcal{X}$ can be defined as

${\frac{Df\left(\mathcal{X}\right)}{D\mathcal{X}}|}_{\overline{\mathcal{X}}}\triangleq {\frac{\partial f\left(\overline{\mathcal{X}}\oplus \delta \mathbf{x}\right)\ominus f\left(\overline{\mathcal{X}}\right)}{\partial \delta \mathbf{x}}|}_{\delta \mathbf{x}=\mathbf{0}},$

(3)

where it should be noted that the function $f\left(\overline{\mathcal{X}}\oplus \delta \mathbf{x}\right)\ominus f\left(\overline{\mathcal{X}}\right)$ of δx has ${\mathbb{R}}^m$ as both its domain and codomain, and can thus be differentiated using any standard technique. With the above general definition of a derivative, it is easy to define the so-called Jacobian of G as J = DExp(x)/Dx, where left/right group Jacobians are obtained with left/right definitions of ⊕ and ⊖.

3.1.4. Composite groups

A composite groups is simply the concatenation of N other Lie groups G₁, …, G _N (Solà et al. (2018)), with elements of the form

$\mathcal{X}=\left({\mathcal{X}}_1,\dots {\mathcal{X}}_N\right)\in G_1\times \cdots \times G_N.$

The group operation, inverse, and identity are defined elementwise. For example,

$\mathcal{X}\circ \mathcal{Y}=\left({\mathcal{X}}_1\circ {\mathcal{Y}}_1,\dots ,{\mathcal{X}}_N\circ {\mathcal{Y}}_N\right).$

Furthermore, defining $\delta \mathbf{x}={\left[\delta {\mathbf{x}}_1^{\mathrm{T}}\dots \delta {\mathbf{x}}_N^{\mathrm{T}}\right]}^{\mathrm{T}}$ the ⊕, operator is given by

$\mathcal{X}\oplus \delta \mathbf{x}=\left({\mathcal{X}}_1\oplus \delta {\mathbf{x}}_1,\dots ,{\mathcal{X}}_N\oplus \delta {\mathbf{x}}_N\right)$

and a similar definition applies to ⊖.

3.2. Maximum a posteriori

Maximum a posteriori (MAP) is the standard approach taken in the robotics literature. Popular algorithms such as the extended Kalman filter (EKF), iterated EKF, sliding-window filter, and batch estimator can all be derived from a MAP approach, thus unifying them under a common theory. Given a statistically independent measurement y with a standard model as in (1), as well as a prior distribution $p\left(\mathcal{X}\right)={\mathcal{N}}_L\left(\check{\mathcal{X}},\check{\mathbf{P}}\right)$ , the estimate $\widehat{\mathcal{X}}$ produced by the MAP approach is

$\begin{array}{ll}\widehat{\mathcal{X}}& =\underset{\mathcal{X}}{\mathrm{argmax}}\hspace{0.17em}p\left(\mathcal{X}|\mathbf{y}\right)\\ & =\underset{\mathcal{X}}{\mathrm{argmax}}\hspace{0.17em}\eta p\left(\mathbf{y}|\mathcal{X}\right)p\left(\mathcal{X}\right)\\ & =\underset{\mathcal{X}}{\mathrm{argmax}}\hspace{0.17em}\eta \mathcal{N}\left(\mathbf{g}\left(\mathcal{X}\right),\mathbf{R}\right){\mathcal{N}}_L\left(\check{\mathcal{X}},\check{\mathbf{P}}\right)\end{array}$

where η is a normalization constant. Equivalently, minimizing the negative logarithm yields a nonlinear least-squares problem of the form

$\begin{array}{ll}\widehat{\mathcal{X}}& =\underset{\mathcal{X}}{\mathrm{argmax}}\frac{1}{2}\mathbf{e}{\left(\mathcal{X}\right)}^{\mathrm{T}}\mathbf{We}\left(\mathcal{X}\right),\end{array}$

(4)

$\begin{array}{ll}\mathbf{e}\left(\mathcal{X}\right)& =\left[\begin{array}{l}\mathcal{X}\ominus \check{\mathcal{X}}\\ \mathbf{y}-\mathbf{g}\left(\mathcal{X}\right)\end{array}\right],\end{array}$

where $\mathbf{W}=\text{diag}\left({\check{\mathbf{P}}}^{-1},{\mathbf{R}}^{-1}\right)$ . Using an on-manifold optimization approach, (4) can be solved by first parameterizing the state with $\mathcal{X}=\widehat{\mathcal{X}}\oplus \delta \mathbf{x}$ and solving the problem

$\delta \widehat{\mathbf{x}}=\underset{\delta \mathbf{x}}{\mathrm{argmax}}\frac{1}{2}\mathbf{e}{\left(\widehat{\mathcal{X}}\oplus \delta \mathbf{x}\right)}^{\mathrm{T}}\mathbf{We}\left(\widehat{\mathcal{X}}\oplus \delta \mathbf{x}\right).$

(5)

Using (3), the Jacobian of the error vector is given by

$$\mathbf{H}\triangleq {\frac{D\mathbf{e}\left(\mathcal{X}\right)}{D\mathcal{X}}|}_{\overline{\mathcal{X}}}={\frac{\partial \mathbf{e}\left(\widehat{\mathcal{X}}\oplus \delta \mathbf{x}\right)}{\partial \delta \mathbf{x}}|}_{\delta \mathbf{x}=\mathbf{0}},$$

and an approximate solution to (5) can be obtained by solving the Gauss–Newton system

$\left({\mathbf{H}}^{\mathrm{T}}\mathbf{WH}\right)\delta \widehat{\mathbf{x}}={\mathbf{H}}^{\mathrm{T}}\mathbf{We}\left(\widehat{\mathcal{X}}\right).$

The above is iterated with $\widehat{\mathcal{X}}\leftarrow \check{\mathcal{X}}\oplus \delta \widehat{\mathbf{x}}$ and initialized with $\widehat{\mathcal{X}}\leftarrow \check{\mathcal{X}}$ . A common approximation for the posterior covariance $\widehat{\mathbf{P}}$ where $p\left(\mathcal{X}|\mathbf{y}\right)\approx {\mathcal{N}}_L\left(\widehat{\mathcal{X}},\widehat{\mathbf{P}}\right)$ is given by $\widehat{\mathbf{P}}={\left({\mathbf{H}}^{\mathrm{T}}\mathbf{WH}\right)}^{-1}$ with H evaluated at $\widehat{\mathcal{X}}$ .

3.3. Covariance intersection

Covariance intersection (CI) is a tool introduced by Julier and Uhlmann (1997) for the purposes of decentralized data fusion under unknown cross-correlations, and can be summarized with the following lemma.

Lemma 1

Consistency of Covariance Intersection. The inequality

$\left[\begin{array}{cc}\frac{1}{w}{\mathbf{\Sigma }}_{xx}& \mathbf{0}\\ \mathbf{0}& \frac{1}{1-w}{\mathbf{\Sigma }}_{yy}\end{array}\right]\ge \left[\begin{array}{ll}{\mathbf{\Sigma }}_{xx}& {\mathbf{\Sigma }}_{xy}\\ {\mathbf{\Sigma }}_{xy}^{\mathrm{T}}& {\mathbf{\Sigma }}_{yy}\end{array}\right],$

(6)

which applies in the positive definite sense, holds for all w ∈ (0, 1), where Σ _xx , Σ _yy , and the right-hand-side of (6) are positive definite.

There are several known strategies for choosing w (Julier (2001)). Following Shalaby et al. (2021b), a fixed value of w = 0.99 is chosen for all the results shown in this paper, as it is a simple approach that yields acceptable results.

4. A toy problem

Consider first one of the simplest multi-robot estimation problems, shown on the left of Figure 1. Two robots are located at positions r₁ and r₂, respectively, and both robots seek to estimate both robots’ positions. By design, each robot carries distinct, conceptually independent estimates, even though their states represent the same true physical variables. This mimics exactly what will occur in implementation, as each robot’s processor will have a live estimate of both robots’ positions. Their state vectors can therefore be defined as

${\mathbf{x}}_1={\left[\begin{array}{ll}{r}_1^{\left[1\right]}& r_2^{\left[1\right]}\end{array}\right]}^{\mathrm{T}},{\mathbf{x}}_2={\left[\begin{array}{ll}{r}_1^{\left[2\right]}& r_2^{\left[2\right]}\end{array}\right]}^{\mathrm{T}},$

where the square bracket superscript (⋅)^[i] is used when necessary to denote Robot i’s estimate or “instance” of a common physical variable. Each robot also collects local measurements from its sensors. Robot 1 is capable of measuring its own position,

$y_1={\mathbf{G}}_1{\mathbf{x}}_1+v_1,v_1\sim \mathcal{N}\left(0,R_1\right),{\mathbf{G}}_1=\left[\begin{array}{ll}1& 0\end{array}\right],$

while Robot 2 is only capable of measuring its position relative to Robot 1,

$y_2={\mathbf{G}}_2{\mathbf{x}}_2+v_2,v_2\sim \mathcal{N}\left(0,R_2\right),{\mathbf{G}}_2=\left[\begin{array}{ll}-1& 1\end{array}\right].$

To keep things simple for this demonstrative problem, robots are assumed to have access to each other’s input measurements, such as wheel odometry. This allows them to predict their state forward in time using a conventional Kalman filter. However, a more communication-efficient solution will be proposed in Section 6. Neither robot is capable of estimating their full state vector from local measurements only, meaning that some form of communication will be required.

To reflect the knowledge that the two robots’ state vectors are physically the same, a key design choice of this paper is to incorporate a pseudomeasurement of the form

${\mathbf{y}}_{12}={\mathbf{x}}_1-{\mathbf{x}}_2+{\mathbf{v}}_{12},{\mathbf{v}}_{12}\sim \mathcal{N}\left(\mathbf{0},\mathbf{\Psi }\right),$

whose “measured” value is always exactly zero. This pseudomeasurement can be viewed as a soft constraint on the problem, inversely weighted by the arbitrary pseudomeasurement covariance Ψ. The estimation problem is now to compute, as accurately as possible, the posterior distribution

$p\left({\mathbf{x}}_1,{\mathbf{x}}_2|y_1,y_2,{\mathbf{y}}_{12}\right).$

4.1. Solution via MAP

Applying MAP to this simplified problem is to say that

${\widehat{\mathbf{x}}}_1,{\widehat{\mathbf{x}}}_2=\underset{{\mathbf{x}}_1,{\mathbf{x}}_2}{\mathrm{argma}}xp\left({\mathbf{x}}_1,{\mathbf{x}}_2|y_1,y_2,{\mathbf{y}}_{12}\right).$

(7)

Assuming that v₁, v₂, v₁₂ are all independent random variables, that is, that p(v₁, v₂, v₁₂) = p(v₁)p(v₂)p(v₁₂), allows the use of Bayes’ rule to write

$\begin{array}{l}p\left({\mathbf{x}}_1,{\mathbf{x}}_2|y_1,y_2,{\mathbf{y}}_{12}\right)=\eta p\left(y_1|{\mathbf{x}}_1\right)p\left(y_2|{\mathbf{x}}_2\right)\\ \times p\left({\mathbf{y}}_{12}|{\mathbf{x}}_1,{\mathbf{x}}_2\right)p\left({\mathbf{x}}_1,{\mathbf{x}}_2\right)\end{array},$

(8)

where η is a normalization constant that does not depend on x₁ or x₂. Next, assume that the prior distributions of the robots are independent and Gaussian, possibly as a result of using CI,

$p\left({\mathbf{x}}_1,{\mathbf{x}}_2\right)=p\left({\mathbf{x}}_1\right)p\left({\mathbf{x}}_2\right)=\mathcal{N}\left({\check{\mathbf{x}}}_1,{\check{\mathbf{P}}}_1\right)\mathcal{N}\left({\check{\mathbf{x}}}_2,{\check{\mathbf{P}}}_2\right).$

(9)

Substituting (9) into (8) and grouping terms into those available to each robot yields

$\begin{array}{l}p\left({\mathbf{x}}_1,{\mathbf{x}}_2|y_1,y_2,{\mathbf{y}}_{12}\right)=\\ \eta p\left({\mathbf{y}}_{12}|{\mathbf{x}}_1,{\mathbf{x}}_2\right)\left(p\left(y_1|{\mathbf{x}}_1\right)p\left({\mathbf{x}}_1\right)\right)\left(p\left(y_2|{\mathbf{x}}_2\right)p\left({\mathbf{x}}_2\right)\right).\end{array}$

(10)

Since the local measurement models are linear, it is straightforward to exactly compute the terms

$p\left(y_i|{\mathbf{x}}_i\right)p\left({\mathbf{x}}_i\right)={\eta }_{i}p\left({\mathbf{x}}_i|y_i\right)={\eta }_i\mathcal{N}\left({\stackrel{\sim }{\mathbf{x}}}_i,{\stackrel{\sim }{\mathbf{P}}}_i\right),i=1,2,$

(11)

using the regular Kalman filter equations. The means and covariances ${\stackrel{\sim }{\mathbf{x}}}_i,{\stackrel{\sim }{\mathbf{P}}}_i$ (with tildes) represent the distribution of each robot’s state conditioned on only local measurements, without the information that the robots’ states are physically the same. Substituting (11) into (10) yields a simplified expression for the posterior, and the optimization problem (7) now leads to the least-squares problem

${\widehat{\mathbf{x}}}_1,{\widehat{\mathbf{x}}}_2=\underset{{\mathbf{x}}_1,{\mathbf{x}}_2}{\mathrm{argmin}}\frac{1}{2}\mathbf{e}{\left({\mathbf{x}}_1,{\mathbf{x}}_2\right)}^{\mathrm{T}}\mathbf{We}\left({\mathbf{x}}_1,{\mathbf{x}}_2\right),$

where

$\begin{array}{c}\mathbf{e}\left({\mathbf{x}}_1,{\mathbf{x}}_2\right)=\left[\begin{array}{cc}\mathbf{1}& \mathbf{0}\\ \mathbf{0}& \mathbf{1}\\ \mathbf{1}& -\mathbf{1}\end{array}\right]\left[\begin{array}{l}{\mathbf{x}}_1\\ {\mathbf{x}}_2\end{array}\right]-\left[\begin{array}{l}{\stackrel{\sim }{\mathbf{x}}}_1\\ {\stackrel{\sim }{\mathbf{x}}}_2\\ \mathbf{0}\end{array}\right]\triangleq \mathbf{Hx}-\mathbf{z},\\ \mathbf{W}=\text{diag}\left({\stackrel{\sim }{\mathbf{P}}}_1^{-1},{\stackrel{\sim }{\mathbf{P}}}_2^{-1},{\mathbf{\Psi }}^{-1}\right).\end{array}$

In this linear case the unique solution $\widehat{\mathbf{x}}$ is given by

$\left[\begin{array}{l}{\widehat{\mathbf{x}}}_1\\ {\widehat{\mathbf{x}}}_2\end{array}\right]={\left({\mathbf{H}}^{\mathrm{T}}\mathbf{WH}\right)}^{-1}{\mathbf{H}}^{\mathrm{T}}\mathbf{Wz},$

(12)

which is also known to be the mean. The relevant matrices expand to

$\begin{array}{ll}{\mathbf{H}}^{\mathrm{T}}\mathbf{WH}& =\left[\begin{array}{ll}{\stackrel{\sim }{\mathbf{P}}}_1^{-1}+{\mathbf{\Psi }}^{-1}& \hspace{1em}-{\mathbf{\Psi }}^{-1}\\ -{\mathbf{\Psi }}^{-1}& {\stackrel{\sim }{\mathbf{P}}}_2^{-1}+{\mathbf{\Psi }}^{-1}\end{array}\right],\\ {\mathbf{H}}^{\mathrm{T}}\mathbf{Wz}& =\left[\begin{array}{l}{\stackrel{\sim }{\mathbf{P}}}_1^{-1}{\stackrel{\sim }{\mathbf{x}}}_1\\ {\stackrel{\sim }{\mathbf{P}}}_2^{-1}{\stackrel{\sim }{\mathbf{x}}}_2\end{array}\right].\end{array}$

The next steps involve various applications of the Sherman-Morrison-Woodbury (SMW) identities to analytically invert the inverse covariance matrix ${\mathbf{H}}^{\mathrm{T}}\mathbf{WH}$ , as well as solve for the solution using (12). The derivation details are omitted for brevity but follow the same steps as (Barfoot 2023, Ch. 3.3.2). The eventual result is

$\begin{array}{ll}\widehat{\mathbf{x}}& \triangleq \left[\begin{array}{l}{\widehat{\mathbf{x}}}_1\\ {\widehat{\mathbf{x}}}_2\end{array}\right]=\left[\begin{array}{l}{\stackrel{\sim }{\mathbf{x}}}_1+{\mathbf{K}}_1\left({\stackrel{\sim }{\mathbf{x}}}_2-{\stackrel{\sim }{\mathbf{x}}}_1\right)\\ {\stackrel{\sim }{\mathbf{x}}}_2+{\mathbf{K}}_2\left({\stackrel{\sim }{\mathbf{x}}}_1-{\stackrel{\sim }{\mathbf{x}}}_2\right)\end{array}\right],\\ \widehat{\mathbf{P}}& \triangleq {\left({\mathbf{H}}^{\mathrm{T}}\mathbf{WH}\right)}^{-1}\\ & =\left[\begin{array}{ll}\left(\mathbf{1}-{\mathbf{K}}_1\right){\stackrel{\sim }{\mathbf{P}}}_1& \hspace{1em}-{\mathbf{K}}_1{\stackrel{\sim }{\mathbf{P}}}_2\\ -{\mathbf{K}}_2{\stackrel{\sim }{\mathbf{P}}}_1& \left(\mathbf{1}-{\mathbf{K}}_2\right){\stackrel{\sim }{\mathbf{P}}}_2\end{array}\right],\\ {\mathbf{K}}_1& \triangleq {\stackrel{\sim }{\mathbf{P}}}_1{\left(\mathbf{\Psi }+{\stackrel{\sim }{\mathbf{P}}}_2+{\stackrel{\sim }{\mathbf{P}}}_1\right)}^{-1},\\ {\mathbf{K}}_2& \triangleq {\stackrel{\sim }{\mathbf{P}}}_2{\left(\mathbf{\Psi }+{\stackrel{\sim }{\mathbf{P}}}_2+{\stackrel{\sim }{\mathbf{P}}}_1\right)}^{-1},\end{array}$

(13)

and furthermore $p\left({\mathbf{x}}_1,{\mathbf{x}}_2|{\mathbf{y}}_1,{\mathbf{y}}_2,{\mathbf{y}}_{12}\right)=\mathcal{N}\left(\widehat{\mathbf{x}},\widehat{\mathbf{P}}\right)$ , which is a standard result from MAP approaches (Barfoot (2023)). The final individual estimates are obtained by marginalizing out the other robots’ states, which is trivial to do in covariance form by simply extracting the corresponding blocks from $\widehat{\mathbf{x}},\widehat{\mathbf{P}}$ , yielding

$\begin{array}{ll}p\left({\mathbf{x}}_1|{\mathbf{y}}_1,{\mathbf{y}}_2,{\mathbf{y}}_{12}\right)& =\mathcal{N}\left({\widehat{\mathbf{x}}}_1,\left(\mathbf{1}-{\mathbf{K}}_1\right){\stackrel{\sim }{\mathbf{P}}}_1\right),\\ p\left({\mathbf{x}}_2|{\mathbf{y}}_1,{\mathbf{y}}_2,{\mathbf{y}}_{12}\right)& =\mathcal{N}\left({\widehat{\mathbf{x}}}_2,\left(\mathbf{1}-{\mathbf{K}}_2\right){\stackrel{\sim }{\mathbf{P}}}_2\right),\end{array}$

The equations (13) has a form similar to a situation where robots simply treated the other robot’s state estimate as a “measurement” of their own state. This is a result that is specific to this simple toy problem, where robots have full state overlap.

Conditioning on the pseudomeasurement y₁₂ has introduced cross-correlation terms in (13), which are feasible to keep track of for this two-robot scenario, but introduce substantial complexity for an arbitrary multi-robot scenario. Therefore, this paper simply employs the CI approximation as required, including for this toy problem for the sake of consistency. Specifically, before each state fusion using (13), inflate the covariance matrices with

${\stackrel{\sim }{\mathbf{P}}}_1\leftarrow \frac{1}{w}{\stackrel{\sim }{\mathbf{P}}}_1,{\stackrel{\sim }{\mathbf{P}}}_2\leftarrow \frac{1}{1-w}{\stackrel{\sim }{\mathbf{P}}}_{2,}$

(14)

where w = 0.99 is used.

Figure 2 shows the estimation error of each robot as multiple pseudomeasurements are fused in succession. The two robots’ estimates not only converge to zero error, but also to a common value, which is the main effect of the pseudomeasurement. A pseudomeasurement covariance of Ψ = 10 ⋅1 was chosen for this simulation to show its effect, but smaller values can be used. Since the expressions in (13) are in covariance form, it is even possible to use Ψ = 0, in which case the two estimates will converge together after the first pseudomeasurement. For the prior distributions, arbitrary Gaussian distributions were chosen, with the initial true states drawn from these distributions.

Figure 2.

Estimation convergence for a single trial of the two-robot toy problem with Ψ = 10 ⋅1. Due to pseudomeasurements, the robot states successfully converge to a common value.

In Figure 3, 100 Monte Carlo trials are performed on a simulation of this toy problem, but extended to four robots using the methods from the next section. The root-mean-squared error (RMSE) and normalized estimation error squared (NEES), calculated as per (Bar-Shalom et al., 2001, Ch. 5.4), are plotted through time. The lines marked “Proposed” fuse pseudomeasurements as described, and use CI before each state fusion. The naive solution is identical, but does not perform the covariance intersection step in (14) before state fusion, thus completely neglects cross-correlations. The centralized solution is simply a Kalman filter with state $\mathbf{x}={\left[r_1{r}_2\right]}^T$ fusing both the measurements y₁ and y₂ using the standard equations. Although the use of CI does introduce error compared to the centralized solution, it is still vastly better than the naive approach.

Figure 3.

Results of 100 Monte Carlo trials for a four-robot version of the toy problem. The top two plots consist of a NEES plot, which is a measure of consistency. The bottom plot is the RMSE of the state. The proposed solution, which performs CI, remains statistically consistent and has reasonably low error in many cases.

5. General problem

Consider now N robots, which communicate in correspondence with an arbitrary undirected graph $\mathcal{G}=\left(\mathcal{V},\mathcal{E}\right)$ where $\mathcal{V}=\left\{1,\dots ,N\right\}$ is the set of nodes or robot IDs, and $\mathcal{E}$ is the set of edges. The robots have states ${\mathcal{X}}_i\in G_i,i\in \mathcal{V}$ belonging to possibly different groups. Pseudomeasurement models ${\mathbf{c}}_{ij}:G_i\times G_j\to {\mathbb{R}}^c$ are now defined in a generic way

${\mathbf{y}}_{i{j}_k}={\mathbf{c}}_{ij}\left({\mathcal{X}}_{i_k},{\mathcal{X}}_{j_k}\right)+{\mathbf{v}}_{i{j}_k},{\mathbf{v}}_{i{j}_k}\sim \mathcal{N}\left(\mathbf{0},\mathbf{\Psi }\right),$

for each pair $\left(i,j\right)\in \mathcal{E}$ . These models are designed by the user, and should correspond with redundant, common, or over-parameterizations of states appearing on different robots, such as when two robots are estimating the same physical quantities. However, these quantities can also differ in the way they are represented from robot to robot, such as each robot resolving the same physical pose in their own frame, and hence, the pseudomeasurement function is left as a general nonlinear function. Examples shall be given in Section 7 and 8.

Using MAP, the general state fusion problem is now

${\widehat{\mathcal{X}}}_{1:N}=\underset{{\mathcal{X}}_{1:N}}{\mathrm{argma}}xp\left({\mathcal{X}}_{1:N}|{\mathbf{y}}_{ij}\right),\text{for} \text{all} \left(i,j\right)\in \mathcal{E}.$

It is easier to instead consider a single pseudomeasurement at a time, such as, without loss of generality, y₁₂. The posterior distribution given only y₁₂ is

$\begin{array}{ll}p\left({\mathcal{X}}_{1:N}|{\mathbf{y}}_{12}\right)& =\eta p\left({\mathbf{y}}_{12}|{\mathcal{X}}_1,{\mathcal{X}}_2\right)p\left({\mathcal{X}}_{1:N}\right)\\ & =\eta \mathcal{N}\left({\mathbf{c}}_{12}\left({\mathcal{X}}_1,{\mathcal{X}}_2\right),\mathbf{\Psi }\right){\displaystyle \prod }_{i=1}^N{\mathcal{N}}_L\left({\stackrel{\sim }{\mathcal{X}}}_i,{\stackrel{\sim }{\mathbf{P}}}_i\right)\end{array}$

(15)

where robot state priors have again been assumed to be independent, as a result of using covariance intersection. Due to this independence, the variables ${\mathcal{X}}_{3:N}$ can be removed from the optimization problem since their optimal values are simply ${\widehat{\mathcal{X}}}_{3:N}={\stackrel{\sim }{\mathcal{X}}}_{3:N}$ and have no effect on ${\mathcal{X}}_1,{\mathcal{X}}_2$ . Minimizing the negative logarithm of (15), omitting terms involving ${\mathcal{X}}_{3:N}$ , leads to a least-squares problem with error vector given by

$\mathbf{e}\left({\mathcal{X}}_1,{\mathcal{X}}_2\right)=\left[\begin{array}{c}{\mathcal{X}}_1\ominus {\stackrel{\sim }{\mathcal{X}}}_1\\ {\mathcal{X}}_2\ominus {\stackrel{\sim }{\mathcal{X}}}_2\\ -{\mathbf{c}}_{12}\left({\mathcal{X}}_1,{\mathcal{X}}_2\right)\end{array}\right],$

and weight $\mathbf{W}=\text{diag}\left({\stackrel{\sim }{\mathbf{P}}}_1^{-1},{\stackrel{\sim }{\mathbf{P}}}_2^{-1},{\mathbf{\Psi }}^{-1}\right)$ . Defining J _i as the group Jacobian associated with G _i , as well as S _i , S _j being the Jacobians of c _ij with respect to ${\mathcal{X}}_i,{\mathcal{X}}_j$ , respectively, the Jacobian of the error vector is

$\mathbf{H}=\left[\begin{array}{cc}{\mathbf{J}}_1^{-1}& \mathbf{0}\\ \mathbf{0}& {\mathbf{J}}_2^{-1}\\ -{\mathbf{S}}_1& -{\mathbf{S}}_2\end{array}\right],$

which is written without arguments $\left({\mathcal{X}}_1,{\mathcal{X}}_2\right)$ for brevity. The relevant terms of the Gauss–Newton system are

$\begin{array}{l}{\mathbf{H}}^{\mathrm{T}}\mathbf{WH}=\\ \left[\begin{array}{cc}{\mathbf{J}}_1^{-\mathrm{T}}{\mathbf{P}}_1^{-1}{\mathbf{J}}_1^{-1}+{\mathbf{S}}_1^{\mathrm{T}}{\mathbf{\Psi }}^{-1}{\mathbf{S}}_1& {\mathbf{S}}_1^{\mathrm{T}}{\mathbf{\Psi }}^{-1}{\mathbf{S}}_2\\ {\mathbf{S}}_2^{\mathrm{T}}{\mathbf{\Psi }}^{-1}{\mathbf{S}}_1& {\mathbf{J}}_2^{-\mathrm{T}}{\mathbf{P}}_2^{-1}{\mathbf{J}}_2^{-1}+{\mathbf{S}}_2^{\mathrm{T}}{\mathbf{\Psi }}^{-1}{\mathbf{S}}_2\end{array}\right],\end{array}$

$\begin{array}{l}{\mathbf{H}}^{\mathrm{T}}\mathbf{We}\left({\mathcal{X}}_1,{\mathcal{X}}_2\right)=\\ \left[\begin{array}{l}{\mathbf{J}}_1^{-\mathrm{T}}{\mathbf{P}}_1^{-1}\left({\mathcal{X}}_1\ominus {\stackrel{\sim }{\mathcal{X}}}_1\right)-{\mathbf{S}}_1^{\mathrm{T}}{\mathbf{\Psi }}^{-1}{\mathbf{c}}_{12}\left({\mathcal{X}}_1,{\mathcal{X}}_2\right)\\ {\mathbf{J}}_2^{-\mathrm{T}}{\mathbf{P}}_2^{-1}\left({\mathcal{X}}_2\ominus {\stackrel{\sim }{\mathcal{X}}}_2\right)-{\mathbf{S}}_2^{\mathrm{T}}{\mathbf{\Psi }}^{-1}{\mathbf{c}}_{12}\left({\mathcal{X}}_1,{\mathcal{X}}_2\right)\end{array}\right],\end{array}$

which, by substantial manipulation with the SMW identities, can be used to analytically compute $\delta \widehat{\mathbf{x}}={\left({\mathbf{H}}^{\mathrm{T}}\mathbf{WH}\right)}^{-1}{\mathbf{H}}^{\mathrm{T}}\mathbf{We}\left({\widehat{\mathcal{X}}}_1,{\widehat{\mathcal{X}}}_2\right)$ , producing on-manifold iterated-EKF-like expressions. The result is

$\begin{array}{ll}\delta {\widehat{\mathbf{x}}}_1& =-{\mathbf{J}}_1\left({\widehat{\mathcal{X}}}_1\ominus {\stackrel{\sim }{\mathcal{X}}}_1\right)+{\mathbf{K}}_1\mathbf{z},\\ \delta {\widehat{\mathbf{x}}}_2& =-{\mathbf{J}}_2\left({\widehat{\mathcal{X}}}_2\ominus {\stackrel{\sim }{\mathcal{X}}}_2\right)+{\mathbf{K}}_2\mathbf{z},\\ {\mathbf{K}}_1& ={\mathbf{J}}_1{\stackrel{\sim }{\mathbf{P}}}_1{\mathbf{J}}_1^{\mathrm{T}}{\mathbf{S}}_1^{\mathrm{T}}{\mathbf{V}}^{-1},\\ {\mathbf{K}}_2& ={\mathbf{J}}_2{\stackrel{\sim }{\mathbf{P}}}_2{\mathbf{J}}_2^{\mathrm{T}}{\mathbf{S}}_2^{\mathrm{T}}{\mathbf{V}}^{-1},\\ \mathbf{z}& =-{\mathbf{c}}_{12}\left({\stackrel{\sim }{\mathcal{X}}}_1,{\stackrel{\sim }{\mathcal{X}}}_2\right)+{\mathbf{S}}_1{\mathbf{J}}_1\left({\widehat{\mathcal{X}}}_1\ominus {\stackrel{\sim }{\mathcal{X}}}_1\right)\\ & +{\mathbf{S}}_2{\mathbf{J}}_2\left({\widehat{\mathcal{X}}}_2\ominus {\stackrel{\sim }{\mathcal{X}}}_2\right),\\ \mathbf{V}& =\mathbf{\Psi }+{\mathbf{S}}_1{\mathbf{J}}_1{\stackrel{\sim }{\mathbf{P}}}_1{\mathbf{J}}_1^{\mathrm{T}}{\mathbf{S}}_1^{\mathrm{T}}+{\mathbf{S}}_2{\mathbf{J}}_2{\stackrel{\sim }{\mathbf{P}}}_2{\mathbf{J}}_2^{\mathrm{T}}{\mathbf{S}}_2^{\mathrm{T}},\end{array}$

(16)

where iteration is done with ${\widehat{\mathcal{X}}}_i\leftarrow {\widehat{\mathcal{X}}}_i\oplus \delta {\widehat{\mathbf{x}}}_i$ , after initialization with ${\widehat{\mathcal{X}}}_i\leftarrow {\stackrel{\sim }{\mathcal{X}}}_i$ , until a convergence condition is met, such as $\delta {\widehat{\mathbf{x}}}_i$ being sufficiently small. The marginal posterior covariances of Robots 1 and 2 are obtained from the corresponding diagonal blocks of ${\left({\mathbf{H}}^{\mathrm{T}}\mathbf{WH}\right)}^{-1}$ , and can be shown to be

$\begin{array}{ll}{\widehat{\mathbf{P}}}_1& =\left(\mathbf{1}-{\mathbf{K}}_1{\mathbf{S}}_1\right){\mathbf{J}}_1{\stackrel{\sim }{\mathbf{P}}}_1{\mathbf{J}}_1^{\mathrm{T}},\\ {\widehat{\mathbf{P}}}_2& =\left(\mathbf{1}-{\mathbf{K}}_2{\mathbf{S}}_2\right){\mathbf{J}}_2{\stackrel{\sim }{\mathbf{P}}}_2{\mathbf{J}}_2^{\mathrm{T}}.\end{array}$

(17)

The above fusion step introduces cross-correlations between the state estimates of ${\mathcal{X}}_1$ and ${\mathcal{X}}_2$ , and hence require a covariance intersection step

${\widehat{\mathbf{P}}}_1\leftarrow \frac{1}{w}{\widehat{\mathbf{P}}}_1,{\widehat{\mathbf{P}}}_2\leftarrow \frac{1}{1-w}{\widehat{\mathbf{P}}}_2,w\in \left(0,1\right).$

The next step is to make the approximation that $p\left({\mathcal{X}}_{1:N}|{\mathbf{y}}_{12}\right)\approx {\prod }_{i=1}^N{\mathcal{N}}_L\left({\widehat{\mathcal{X}}}_i,{\widehat{\mathbf{P}}}_i\right)$ , and proceed with the fusion of a second pseudomeasurement, using $p\left({\mathcal{X}}_{1:N}|{\mathbf{y}}_{12},{\mathbf{y}}_{13}\right)=\eta p\left({\mathbf{y}}_{13}|{\mathcal{X}}_1,{\mathcal{X}}_3\right)p\left({\mathcal{X}}_{1:N}|{\mathbf{y}}_{12}\right)$ . This new posterior can again be approximated as Gaussian using the expressions in (16) and (17), and the process is repeated until all pseudomeasurements are incorporated. The order in which the pseudomeasurements are fused is arbitrary and will correspond to the times at which information is shared between robots.

If only one iteration is performed, equations (16) and (17) simplify to the on-manifold EKF equations since ${\mathbf{J}}_i\left({\stackrel{\sim }{\mathcal{X}}}_i\ominus {\stackrel{\sim }{\mathcal{X}}}_i\right)={\mathbf{J}}_i\left(\mathbf{0}\right)=\mathbf{1}$ . Much like the EKF is often sufficient compared to the iterated EKF, performing a single iteration in this multi-robot case is also sufficient for some problems. Algorithm 1 summarizes the general-purpose decentralized estimation algorithm from the point of view of an arbitrary robot. The algorithm is presented in a callback format, describing how the current state estimate is updated when various events occur.

A feature of Algorithm 1 is that the robots can share their state information at anytime, with performance improving the more often sharing occurs, at the cost of increased communication bandwidth. After states are shared and pseudomeasurements are fused, the common states between robots will naturally drift due to sensor noise, until the next pseudomeasurement re-synchronizes the common states. This therefore becomes a tunable trade-off between estimation accuracy and communication bandwidth, which needs to be evaluated for a specific problem. In this paper’s experiments, sharing is done at a set frequency of 10 Hz. Regarding computational complexity, Algorithm 1 will share a nearly identical runtime to a Kalman filter.

5.1. An observability test

In a decentralized state estimation context, observability refers to the ability for each robot to uniquely determine their state trajectory, given the inputs and measurements obtained by all robots. Determining observability for a decentralized estimator is non-trivial due to the need to capture the communication topology within the test itself. To illustrate this, consider again the linear toy problem with two robots from Section 4. A naive approach to determining observability would be to construct one “total state” $\mathbf{x}={\left[{\mathbf{x}}_1{\mathbf{x}}_2\right]}^{\mathrm{T}}$ , and to perform a standard observability test on this augmented system using the collected sensor measurements from both robots y₁, y₂. However, such a test would falsely conclude that the system is unobservable, whereas the results from Section 4 clearly show successful estimation. The pseudomeasurement must additionally be incorporated into the test to give the correct result, as the robots are reliant on communication to attain observability of their individual states. To the best of the authors’ knowledge, existing observability tests do not take into the full problem scope that this paper is concerned with. A decentralized observability test is presented in Pilloni et al. (2013) for linear-time-invariant systems with unknown input, but assumes each node is observable. The observability analysis in Huang et al. (2011) has a similar approach to this paper, but is specific to their system where robots have identical state definitions.

An advantage of the proposed approach is that the effects of communication of observability can be accurately captured by incorporating the pseudomeasurements themselves into a standard observability test. Concretely, for nonlinear systems, a local observability test can be formed by considering the MAP problem on an entire trajectory simultaneously, but without prior information on the initial state (Psiaki (2013)). Applying this to the multi-robot system, let the bolded ${\mathcal{X}}_k=\left({\mathcal{X}}_{1_k},\dots ,{\mathcal{X}}_{N_k}\right)$ denote the “total state” of all robots at time step k, and ${\mathbf{y}}_k={\left[{\mathbf{y}}_{1_k}^{\mathrm{T}}\dots {\mathbf{y}}_{N_k}^{\mathrm{T}}\right]}^{\mathrm{T}}$ denote a stacked vector containing all the robots’ local measurements at time step k. Let ${\mathbf{\psi }}_k={\left[\dots {\mathbf{y}}_{ij}^{\mathrm{T}}\dots \right]}^{\mathrm{T}},\left(i,j\right)\in \mathcal{E}$ denote the stacked pseudomeasurements between all robots at time step k. The MAP problem is

${\widehat{\mathcal{X}}}_{0:K}=\underset{ {\mathcal{X}}_{0:K}}{\mathrm{argma}}xp\left( {\mathcal{X}}_{0:K}|{\mathbf{y}}_{0:K},{\mathbf{\psi }}_{0:K}\right)$

with

$\begin{array}{l}p\left( {\mathcal{X}}_{0:K}|{\mathbf{y}}_{0:K},{\mathbf{\psi }}_{0:K}\right)=\\ \eta {\displaystyle \prod }_{k=0}^{K}p\left({\mathbf{y}}_0| {\mathcal{X}}_0\right)p\left({\mathbf{\psi }}_0| {\mathcal{X}}_0\right){\displaystyle \prod }_{k=1}^{K}p\left( {\mathcal{X}}_k| {\mathcal{X}}_{k-1}\right),\end{array}$

leading to a nonlinear least-squares problem with weight W and error vector

$\begin{array}{ll}\mathbf{e}\left( {\mathcal{X}}_{0:K}\right)& ={\left[\dots {\mathbf{e}}_{u,k}\dots {\mathbf{e}}_{y,k}\dots {\mathbf{e}}_{\psi ,k}\dots \right]}^{\mathrm{T}},\\ {\mathbf{e}}_{u,k}& =\left[\dots {\left({\mathcal{X}}_{i_k}\ominus \mathbf{f}\left({\mathcal{X}}_{i_{k-1}},{\mathbf{u}}_{i_{k-1}}\right)\right)}^{\mathrm{T}}\dots \right],\\ {\mathbf{e}}_{y,k}& =\left[\dots {\left({\mathbf{y}}_{i_k}-\mathbf{g}\left({\mathcal{X}}_{i_k}\right)\right)}^{\mathrm{T}}\dots \right],i=1,\dots ,N,\\ {\mathbf{e}}_{\psi ,k}& =\left[\dots -{\mathbf{c}}_{ij}{\left({\mathcal{X}}_{i_k},{\mathcal{X}}_{j_k}\right)}^{\mathrm{T}}\dots \right],\left(i,j\right)\in \mathcal{E}.\end{array}$

The error Jacobian is

$\mathbf{H}=\left[\begin{array}{cccc}-{\mathbf{F}}_0& \mathbf{1}& & \\ & \ddots & \ddots & \\ & & -{\mathbf{F}}_{K-1}& \mathbf{1}\\ -{\mathbf{G}}_0& & & \\ & -{\mathbf{G}}_1& & \\ & & \ddots & \\ & & & -{\mathbf{G}}_K\\ -{\mathbf{\Phi }}_0& & & \\ & -{\mathbf{\Phi }}_1& & \\ & & \ddots & \\ & & & -{\mathbf{\Phi }}_K\end{array}\right],$

$\begin{array}{ll}{\mathbf{F}}_k& =\text{diag}\left(\dots ,\frac{D\mathbf{f}\left({\mathcal{X}}_{i_k},{\mathbf{u}}_{i_k}\right)}{D{\mathcal{X}}_{i_k}},\dots \right),\\ {\mathbf{G}}_k& =\text{diag}\left(\dots ,\frac{D\mathbf{g}\left({\mathcal{X}}_{i_k}\right)}{D{\mathcal{X}}_{i_k}},\dots \right),i=1,\dots ,N,\\ {\mathbf{\Phi }}_k& =-\frac{D{\mathbf{e}}_{\psi ,k}\left({\mathcal{X}}_k\right)}{D {\mathcal{X}}_k},\end{array}$

with all undisplayed entries in H equal to zero. For the solution to the MAP problem to be unique, then to first order, $\left({\mathbf{H}}^{\mathrm{T}}\mathbf{WH}\right)$ must be invertible, and thus full rank. Fortunately, W is always positive definite regardless of any cross-correlations that would add off-diagonal entries. Hence,

$\text{rank}\left({\mathbf{H}}^{\mathrm{T}}\mathbf{WH}\right)=\text{rank}\left(\mathbf{H}\right),$

and it is thus required that H be full column rank. This implies that the proposed observability test is unaffected by the approximation induced by CI, or any cross-correlation terms that may or may not be successfully tracked. In a similar way to (Barfoot 2023, Ch 3.1.4), it can be shown that by performing a variety of elementary row/column operations, the rank of H is equivalent to the rank of

$\mathcal{O}=\left[\begin{array}{c}{\mathbf{M}}_0\\ {\mathbf{M}}_1{\mathbf{F}}_0\\ ⋮\\ {\mathbf{M}}_K{\mathbf{F}}_{K-1}\dots {\mathbf{F}}_0\end{array}\right],{\mathbf{M}}_k=\left[\begin{array}{l}{\mathbf{G}}_k\\ {\mathbf{\Phi }}_k\end{array}\right].$

Hence, if $\mathcal{O}$ has maximum rank, the solution to the MAP problem is locally unique, and the system is said to be observable. Note that this test easily allows for time-varying graphs, which would yield a different Φ _k for each time step k.

6. Efficient odometry sharing using preintegration

Many problems, especially those where robots estimate their neighbors’ positions, will require robots to have access to their neighbors’ process model input values u. Previously in this paper, it has been assumed that all robots have unrestricted access to each other’s inputs. In robot state estimation applications, the input is often the odometry measurements, such as wheel encoder or IMU measurements. These can occur at frequencies of 100–1000 Hz, and can therefore be infeasible to share in real time, especially if multiple robots are to simultaneously share measurements at high frequency. This could quickly reach a bandwidth limit on the common communication channel, such as ultra-wideband radio. While Xu et al. (2020) solve this problem by directly sharing pose changes between two points in time, this violates statistical independence assumptions and leads to inconsistent estimates.

The proposed solution to this problem is to use preintegration. That is, robots will instead share preintegrated input measurements over an arbitrary duration of time instead of individual input measurements. Specifically, consider the following generic process model

${\mathcal{X}}_k=\mathbf{f}\left({\mathcal{X}}_{k-1},{\mathbf{u}}_{k-1},{\mathbf{w}}_{k-1}\right).$

The action of preintegration is to directly iterate this process model by repeated compositions in order to, after algebraic manipulation, generate a new preintegrated process model f _pq that relates two states at arbitrary time steps k = p and k = q. That is,

$\begin{array}{ll}{\mathcal{X}}_{p+1}& =\mathbf{f}\left({\mathcal{X}}_p,{\mathbf{u}}_p,{\mathbf{w}}_p\right),\\ {\mathcal{X}}_{p+2}& =\mathbf{f}\left(\mathbf{f}\left({\mathcal{X}}_p,{\mathbf{u}}_p,{\mathbf{w}}_p\right),{\mathbf{u}}_{p+1},{\mathbf{w}}_{p+1}\right),\\ & ⋮\\ {\mathcal{X}}_q& =\mathbf{f}\left(\mathbf{f}\left(\dots \mathbf{f}\left({\mathcal{X}}_p,{\mathbf{u}}_p,{\mathbf{w}}_p\right)\dots \right),{\mathbf{u}}_{q-1},{\mathbf{w}}_{q-1}\right)\\ & \triangleq {\mathbf{f}}_{pq}\left({\mathcal{X}}_i,\mathrm{\Delta }{\mathcal{X}}_{pq}\right)\oplus {\mathbf{w}}_{pq,}\end{array}$

(18)

where $\mathrm{\Delta }{\mathcal{X}}_{pq}$ is the relative motion increment (RMI), which in general may also belong to a Lie group, and ${\mathbf{w}}_{pq}\sim \mathcal{N}\left(\mathbf{0},{\mathbf{Q}}_{pq}\right)$ is the preintegrated noise. The advantages of preintegration will stem from the careful choice of RMI definition, which is ideally done such that the RMI has the following properties.

• 1. The RMI is determined from the input measurements exclusively, and is independent of the state estimate:

$\mathrm{\Delta }{\mathcal{X}}_{pq}=\mathrm{\Delta }{\mathcal{X}}_{pq}\left({\mathbf{u}}_{p:q-1}\right).$

• 2. Far fewer numbers are required to represent the RMI than the (q − p) raw measurements that occurred during the preintegration interval:

$\dim \left(\mathrm{\Delta }{\mathcal{X}}_{pq}\right)\ll \left(q-p\right)\dim \left({\mathbf{u}}_k\right).$

If the above points are true, communicating $\mathrm{\Delta }{\mathcal{X}}_{pq},{\mathbf{Q}}_{pq}$ instead of u_p:q−1 will not only reduce the communication cost, but will also result in a fixed message size and ability to directly predict the state forward over a long duration of time, instead of sequentially processing the measurements. Note that it is not always possible to define an RMI such that (18) holds exactly. However, it turns out that many common process models in robotics are amenable to preintegration (Barrau and Bonnabel (2019); Eckenhoff et al. (2019)), and furthermore are typically extremely fast to preintegrate incrementally as input measurements are obtained (Forster et al. (2017)). In other words, there exists a function Increment(⋅) defined such that it satisfies

$\mathrm{\Delta }{\mathcal{X}}_{pq},{\mathbf{Q}}_{pq}=$ Increment $\left(\mathrm{\Delta }{\mathcal{X}}_{p:q-1},{\mathbf{Q}}_{p:q-1},{\mathbf{u}}_{q-1}\right)$ .

A few examples now follow, which describe concrete implementations of $\mathrm{\Delta }{\mathcal{X}}_{pq}$ , f _pq (⋅), and Increment(⋅).

Example 1

Linear preintegration. The linear process model

${\mathbf{x}}_k={\mathbf{F}}_{k-1}{\mathbf{x}}_{k-1}+{\mathbf{L}}_{k-1}{\mathbf{u}}_{k-1}$

can be directly iterated to yield

$\begin{array}{ll}{\mathbf{x}}_q& =\left({\displaystyle \prod }_{k=p}^{q-1}{\mathbf{F}}_k\right){\mathbf{x}}_p+{\displaystyle \sum }_{k=p}^{q-1}\left({\displaystyle \prod }_{\mathcal{\ell }=k+1}^{q-1}{\mathbf{F}}_{\mathcal{\ell }}\right){\mathbf{L}}_k{\mathbf{u}}_k\\ & \triangleq {\mathbf{F}}_{pq}{\mathbf{x}}_p+\mathrm{\Delta }{\mathbf{x}}_{pq},\end{array}$

(19)

where (19) defines f(⋅). Assuming that noise enters the model additively through the input ${\mathbf{u}}_k={\overline{\mathbf{u}}}_k+{\mathbf{w}}_k$ , where ${\mathbf{w}}_k\sim \mathcal{N}\left(\mathbf{0},{\mathbf{Q}}_k\right)$ , (19) becomes ${\mathbf{x}}_j={\mathbf{F}}_{pq}{\mathbf{x}}_i+\mathrm{\Delta }{\overline{\mathbf{x}}}_{pq}+{\mathbf{w}}_{pq}$ where

$\begin{array}{ll}{\mathbf{w}}_{pq}& \triangleq {\displaystyle \sum }_{k=p}^{q-1}\left({\displaystyle \prod }_{\mathcal{\ell }=k+1}^{q-1}{\mathbf{F}}_{\mathcal{\ell }}\right){\mathbf{L}}_k{\mathbf{w}}_k,\\ & ={\mathbf{F}}_{q-1}{\mathbf{w}}_{pq-1}+{\mathbf{L}}_{q-1}{\mathbf{w}}_{q-1}.\end{array}$

The RMI Δx _pq and corresponding covariance are therefore built incrementally with

$\begin{array}{ll}\mathrm{\Delta }{\mathbf{x}}_{pq}& ={\mathbf{F}}_{q-1}\mathrm{\Delta }{\mathbf{x}}_{pq-1}+{\mathbf{L}}_{q-1}{\mathbf{u}}_{q-1},\\ {\mathbf{Q}}_{pq}& ={\mathbf{F}}_{q-1}{\mathbf{Q}}_{pq-1}{\mathbf{F}}_{q-1}^{\mathrm{T}}+{\mathbf{L}}_{q-1}{\mathbf{Q}}_{q-1}{\mathbf{L}}_{q-1,}^{\mathrm{T}}\end{array}$

which together define the Increment(⋅) function.

Example 2

Wheel odometry preintegration on SE(2). Given a robot pose T ∈ SE(2), the wheel odometry process model is given by

${\mathbf{T}}_k={\mathbf{T}}_{k-1}\mathrm{Exp}\left(\mathrm{\Delta }t{\mathbf{u}}_{k-1}\right),$

where $\mathbf{u}={\left[\omega v0\right]}^{\mathrm{T}}$ , ω is the robot’s heading rate-of-change, and v is its forward velocity in its own body frame. Direct iteration yields the preintegrated process model f _pq (⋅) given by

${\mathbf{T}}_q={\mathbf{T}}_p\underset{\triangleq \mathrm{\Delta }{\mathbf{T}}_{pq}}{\underbrace{{\displaystyle \prod }_{k=p}^{q-1}\mathrm{Exp}\left(\mathrm{\Delta }t{\mathbf{u}}_k\right)}}.$

Noise ${\mathbf{w}}_k\sim \mathcal{N}\left(\mathbf{0},{\mathbf{Q}}_k\right)$ is again assumed to enter additively through the input, and a series of first-order approximations lead to

$\begin{array}{ll}\mathrm{\Delta }{\mathbf{T}}_{pq}& ={\displaystyle \prod }_{k=p}^{q-1}\mathrm{Exp}\left(\mathrm{\Delta }t\left({\overline{\mathbf{u}}}_k+{\mathbf{w}}_k\right)\right)\\ & \approx {\displaystyle \prod }_{k=p}^{q-1}\mathrm{Exp}\left(\mathrm{\Delta }t{\overline{\mathbf{u}}}_k\right)\mathrm{Exp}\left(\mathrm{\Delta }t{\mathbf{J}}_k{\mathbf{w}}_k\right)\\ & \approx \mathrm{\Delta }{\overline{\mathbf{T}}}_{pq}\underset{\mathrm{Exp}\left({\mathbf{w}}_{pq}\right)}{\underbrace{{\displaystyle \prod }_{k=p}^{q-1}\mathrm{Exp}\left(\mathrm{\Delta }t\mathbf{Ad}\left(\mathrm{\Delta }{\mathbf{T}}_{k+1j}^{-1}\right){\mathbf{J}}_k{\mathbf{w}}_k\right),}}\end{array}$

where ${\mathbf{J}}_k\triangleq \mathbf{J}\left(\mathrm{\Delta }t{\overline{\mathbf{u}}}_k\right)$ is the right Jacobian of SE(2). Having identified an expression for Exp(w _pq ), under the assumption that w _pq is small,

$\begin{array}{ll}{\mathbf{w}}_{pq}& \approx {\displaystyle \sum }_{k=p}^{q-1}\mathrm{\Delta }t\mathbf{Ad}\left(\mathrm{\Delta }{\mathbf{T}}_{k+1q}^{-1}\right){\mathbf{J}}_k{\mathbf{w}}_k\\ & ={\displaystyle \sum }_{k=p}^{q-2}\mathrm{\Delta }t\mathbf{Ad}\left(\mathrm{\Delta }{\mathbf{T}}_{k+1q}^{-1}\right){\mathbf{J}}_k{\mathbf{w}}_k\\ & +\mathrm{\Delta }t\underset{\mathbf{1}}{\underbrace{\mathbf{Ad}\left(\mathrm{\Delta }{\mathbf{T}}_{qq}^{-1}\right)}}{\mathbf{J}}_{q-1}{\mathbf{w}}_{q-1}\\ & =\underset{\triangleq {\mathbf{F}}_{q-1}}{\underbrace{\mathbf{Ad}\left(\mathrm{\Delta }{\mathbf{T}}_{q-1q}^{-1}\right)}}{\mathbf{w}}_{pq-1}\\ & +\underset{\triangleq {\mathbf{L}}_{q-1}}{\underbrace{\mathrm{\Delta }t{\mathbf{J}}_{q-1}}}{\mathbf{w}}_{q-1},\end{array}$

and the defining operations of the Increment(⋅) function follow,

$\begin{array}{ll}\mathrm{\Delta }{\mathbf{T}}_{pq}& =\mathrm{\Delta }{\mathbf{T}}_{pq-1}\mathrm{Exp}\left(\mathrm{\Delta }t{\mathbf{u}}_k\right),\\ {\mathbf{Q}}_{pq}& ={\mathbf{F}}_{q-1}{\mathbf{Q}}_{pq-1}{\mathbf{F}}_{q-1}^{\mathrm{T}}+{\mathbf{L}}_{q-1}{\mathbf{Q}}_{q-1}{\mathbf{L}}_{q-1}^{\mathrm{T}}.\end{array}$

Example 3

IMU preintegration. Being the most well-known usage of preintegration, a complete reference for IMU preintegration on the $SO\left(3\right)\times {\mathbb{R}}^3\times {\mathbb{R}}^3$ manifold can be obtained from Forster et al. (2017), and alternatively for the SE₂(3) group from Brossard et al. (2022); Barfoot (2023). Either approach can be used with the framework in this paper. However in Section 8 of this paper, T _wi ∈ SE₂(3) matrices are used to represent the extended pose of Robot i relative to an inertial world frame w. Following Shalaby et al. (2023) and Brossard et al. (2022), it can be shown that the discrete-time IMU kinematic equations can be written in the form

${\mathbf{T}}_{w{i}_k}={\mathbf{G}}_{k-1}{\mathbf{T}}_{w{i}_{k-1}}{\mathbf{U}}_{k-1},$

(20)

where

$\begin{array}{ll}{\mathbf{T}}_{w{i}_k}& =\left[\begin{array}{ccc}\mathbf{C}& \mathbf{v}& \mathbf{r}\\ \mathbf{0}& 1& 0\\ \mathbf{0}& 0& 1\end{array}\right]\\ {\mathbf{G}}_{k-1}& =\left[\begin{array}{ccc}\mathbf{1}& \mathrm{\Delta }t\mathbf{g}& \hspace{1em}-\frac{\mathrm{\Delta }{t}^2}{2}\mathbf{g}\\ 0& 1& \hspace{1em}-\mathrm{\Delta }t\\ 0& 0& 1\end{array}\right]\\ {\mathbf{U}}_{k-1}& =\left[\begin{array}{ccc}\exp \left(\mathrm{\Delta }t{\mathbf{\omega }}^{\wedge }\right)& \mathrm{\Delta }t\mathbf{J}\left(\mathrm{\Delta }t\mathbf{\omega }\right)\mathbf{a}& \frac{\mathrm{\Delta }{t}^2}{2}\mathbf{N}\left(\mathrm{\Delta }t\mathbf{\omega }\right)\mathbf{a}\\ 0& 1& \mathrm{\Delta }t\\ 0& 0& 1\end{array}\right]\\ \mathbf{N}\left(\varphi \right)& =\mathbf{z}{\mathbf{z}}^{\mathrm{T}}+2\left(\frac{1}{\varphi }-\frac{\sin \hspace{0.17em}\varphi }{{\varphi }^2}\right){\mathbf{z}}^{\wedge }+2\frac{\cos \hspace{0.17em}\varphi -1}{{\varphi }^2}{\mathbf{z}}^{\wedge }{\mathbf{z}}^{\wedge }\\ \varphi & =\Vert \varphi \Vert ,\mathbf{z}=\varphi /\varphi \end{array}$

and C ∈ SO(3), v, r respectively represent attitude, velocity, position relative to the world frame w, ω is the IMU’s unbiased gyro measurement, a is the IMU’s unbiased accelerometer measurement, g is the gravity vector resolved in frame w, and J( ϕ ) is the left Jacobian of SO(3). Preintegration of these kinematics is easily achieved by direct iteration with

$\begin{array}{ll}{\mathbf{T}}_{w{i}_q}& =\left({\displaystyle \prod }_{k=p}^{q-1}{\mathbf{G}}_{k-1}\right){\mathbf{T}}_{w{i}_p}\left({\displaystyle \prod }_{k=p}^{q-1}{\mathbf{U}}_{k-1}\right)\\ & \triangleq \mathrm{\Delta }{\mathbf{G}}_{pq}{\mathbf{T}}_{w{i}_p}\mathrm{\Delta }{\mathbf{U}}_{pq}\end{array}$

where ΔU _pq is the RMI, and the noise statistics can be propagated through this preintegration process by a standard linearization procedure. Shalaby et al. (2023) present the details of the noise propagation, as well as how to adapt this formulation for relative poses.

6.1. Multi-robot preintegration

In the context of multi-robot estimation problems, an individual robot’s process model may involve the input values of many neighboring robots. To reflect this, rewrite the process model for Robot i as

${\mathcal{X}}_{i_k}=\mathbf{f}\left({\mathcal{X}}_{i_{k-1}},{\mathbf{u}}_{i_{k-1}},{\mathbf{u}}_{j_{k-1}}\right),j\in {\mathcal{N}}_i,$

(21)

where ${\mathbf{u}}_{i_k}$ denotes an input measured by Robot i and ${\mathcal{N}}_i$ denotes the set of neighbor IDs of Robot i. The preintegrated process model would now be written as

${\mathcal{X}}_{i_q}={\mathbf{f}}_{pq}\left({\mathcal{X}}_{i_p},\mathrm{\Delta }{\mathcal{X}}_{i_{pq}},\mathrm{\Delta }{\mathcal{X}}_{j_{pq}}\right),j\in {\mathcal{N}}_i$

(22)

where $\mathrm{\Delta }{\mathcal{X}}_{i_{pq}}$ denotes an RMI calculated from the input measurements of Robot i.

A complication is that the RMIs from neighboring Robots $\mathrm{\Delta }{\mathcal{X}}_{j_{pq}}$ are only available asynchronously, meaning it is not always possible to evaluate (22) directly. To deal with this, assume that the state can be temporarily partially propagated with a value of 0 for the neighbor’s input, and then fully propagated once the RMI is shared. That is, assume the preintegrated process model f _pq is compatible with

$\begin{array}{l}{\mathfrak{X}}_{i_{k-1}}=\mathbf{f}\left({\mathcal{X}}_{i_{k-2}},{\mathbf{u}}_{i_{k-2}},\mathbf{0}\right),\\ {\mathfrak{X}}_{i_k}=\mathbf{f}\left({\mathfrak{X}}_{i_{k-1}},{\mathbf{u}}_{i_{k-1}},\mathbf{0}\right),\\ {\mathcal{X}}_{i_k}={\mathbf{f}}_{pq}\left({\mathfrak{X}}_{i_k},\mathcal{I},\mathrm{\Delta }{\mathcal{X}}_{j_{pq}}\right),\end{array}$

where $\mathcal{I}$ is the “identity” or “zero” RMI constructed from input values of zero, and f is defined as per (21) with u _j substituted for 0. The variable ${\mathfrak{X}}_{i_k}$ represents an intermediate, non-physical state that is propagated using the process model without input information from neighboring robots. Sometime later, at arbitrary time step k = q, the RMI from a neighboring robot $\mathrm{\Delta }{\mathcal{X}}_{j_{pq}}$ is received and this intermediate state ${\mathfrak{X}}_{i_q}$ is propagated back into a physically meaningful quantity ${\mathcal{X}}_{i_q}$ . A concrete example of this asynchronous intermediate state updating is shown in Section 7, and a summary is shown in Algorithm 2.

6.2 Estimating input biases

For some problems, it may be desired to estimate an input bias b as part of the overall state $\mathcal{X}$ , (a) setup commonly occurring in inertial navigation where accelerometer and rate gyro biases are estimated. The difficulty lies in the frequent inability to express RMIs independently of the bias values, thus leaving RMIs in the form of

$\mathrm{\Delta }{\mathcal{X}}_{pq}\left({\mathbf{u}}_{p:q-1},{\mathbf{b}}_{p:q-1}\right)$

In the context of the multi-robot estimation scheme shown in Algorithm 2, computing RMIs this way causes inconsistency in the filter, since the RMIs are now correlated with the robot states. Accounting for this would require maintaining the cross-correlation between a robot’s state and their neighbors’ biases.

A simpler alternate solution is to have robots estimate their neighbors’ input biases in addition to their own. This requires to exploit the fact that biases are usually modeled to follow a random walk, and therefore have a constant mean in the absence of any correcting information. This motivates the approximation b _p ≈ b_p+1 ≈ … ≈ b _q and hence

$\mathrm{\Delta }{\mathcal{X}}_{pq}\left({\mathbf{u}}_{p:q-1},{\mathbf{b}}_{p:q-1}\right)\approx \mathrm{\Delta }{\mathcal{X}}_{pq}\left({\mathbf{u}}_{p:q-1},{\mathbf{b}}_q\right).$

When robots receive input measurements, they increment their RMIs with raw (biased) inputs to produce biased RMIs $\mathrm{\Delta }{\mathcal{X}}_{pq}\left({\mathbf{u}}_{p:q-1},\mathbf{0}\right)$ . At an appropriate time, they share their current biased RMIs, which is corrected for bias by the receiving robot using the first-order approximation

$\begin{array}{l}\mathrm{\Delta }{\mathcal{X}}_{pq}\left({\mathbf{u}}_{p:q-1},{\mathbf{b}}_q\right)\approx \mathrm{\Delta }{\mathcal{X}}_{pq}\left({\mathbf{u}}_{p:q-1},\mathbf{0}\right)\oplus {\mathbf{B}}_{pq}{\mathbf{b}}_q,\\ {\mathbf{B}}_{pq}\triangleq \frac{D}{D{\mathbf{b}}_q}\left(\mathrm{\Delta }{\mathcal{X}}_{pq}\left({\mathbf{u}}_{p:q-1},{\mathbf{b}}_q\right)\right),\end{array}$

(23)

where B _pq is defined as the bias Jacobian. Equation (23) is an approximation that contributes unmodeled errors to the estimation problem, relying on an assumption that b _q is small. This can be enabled by a proper offline calibration procedure that removes any large bias, and only small deviations from this are estimated online. Such a procedure is used for the quadcopter problem in Section 8, where good, consistent estimation results are still obtained despite the approximation in (23).

6.3. Autoencoding covariance matrices

As shown in Algorithm 2, predicting state estimates that are a function of neighbor inputs requires the RMI $\mathrm{\Delta }{\mathcal{X}}_{pq}$ along with a corresponding covariance Q _pq . As is, these two quantities must be shared between robots. This section proposes an optional method that further reduces communication costs by eliminating the requirement to share the preintegrated covariance Q _pq .

The key insight is that $\mathrm{\Delta }{\mathcal{X}}_{pq}\left({\mathbf{u}}_{p:q-1}\right)$ and Q _pq (u_p:q−1) are both calculated from the same input values u_p:q−1. Hence, if an alternate mapping ${\mathbf{Q}}_{pq}=\mathbf{h}\left(\mathrm{\Delta }{\mathcal{X}}_{pq}\right)$ existed, then it would be sufficient to share $\mathrm{\Delta }{\mathcal{X}}_{pq}$ only, and the receiving robot could infer Q _pq directly from the RMI. In the absence of analytic expressions for h(⋅), this paper approximates the function with a neural network, trained on purely synthetic RMI covariance pairs. An additional complication is that such a function h(⋅) may not always exist since an RMI can correspond to many possible covariances depending on the input values. In this case h(⋅) is not a true function since it is one-to-many, and its definition is modified to also accept a low-dimensional encoding e (Q _pq ), leading to ${\mathbf{Q}}_{pq}=\mathbf{h}\left(\mathrm{\Delta }{\mathcal{X}}_{pq},\mathit{e}\left({\mathbf{Q}}_{pq}\right)\right)$ . This leads to an architecture here referred to as mean-assisted autoencoding, depicted in Figure 4.

Figure 4.

Concept diagram of mean-assisted autoencoder.

The flattened lower-triangular half of Q _pq is given to a simple fully connected encoder network with GELU activation functions and a single hidden layer with 256 neurons. The output of this network is the encoding, which can be as small as one or two numbers. This encoding is then concatenated with a parameterization of the RMI $\mathrm{\Delta }{\mathcal{X}}_{pq}$ and fed to a similar decoder network, again with a single 256-neuron hidden layer. The decoder network outputs the flattened lower-triangular half of a Cholesky decomposition L, which is used to reconstruct the matrix using ${\widehat{\mathbf{Q}}}_{pq}=\mathbf{L}{\mathbf{L}}^{\mathrm{T}}$ . The covariance matrix is guaranteed to be positive definite as long as the diagonal elements of L are non-zero, which is extremely unlikely to occur in practice. For training, the loss function simply uses the Frobenius norm,

$\mathcal{L}\left({\mathbf{Q}}_{pq},{\widehat{\mathbf{Q}}}_{pq}\right)={\Vert {\mathbf{Q}}_{pq}-{\widehat{\mathbf{Q}}}_{pq}\Vert }_{\mathrm{F}}.$

Figure 5 shows the training convergence history for various encoding sizes, applied to IMU preintegration. The Adam optimizer is chosen with an initial learning rate of 10⁻³ that is scheduled to decrease once the loss plateaus. The training process takes just a few minutes on a laptop CPU to achieve less than 1% average reconstruction error on a validation dataset from experimental data. The “No encoding” baseline represents an architecture where the decoder network attempts to infer the covariance matrix directly from the RMI itself, without using the output of the encoder network, and hence eliminating the need to communicate RMI covariance information. However, doing this yields substantially higher error than when an encoding is provided. A visualization of the reconstructed covariance matrices can be seen in Figure 6, using an encoding size of only one number. Section 8 will employ this method “in the loop” for a real quadcopter problem. It will be shown that the reconstruction error is so small that the impact on the estimation results are negligible.

Figure 5.

Mean percentage reconstruction error throughout training for various encoding sizes including no encoding. A single encoding number is sufficient to achieve less than 1% reconstruction error on average.

Figure 6.

Visualization of preintegrated IMU noise covariance matrices along with reconstruction using mean-assisted autoencoding.

The training data is purely synthetic, where RMIs are constructed from a random amount of random IMU measurements, with values covering the realistic range of real IMU measurements. Since the length of the dataset is infinite, the risk of overfitting is completely eliminated, as long as the real IMU measurements lie within the range of randomly generated values. In fact, the networks immediately generalize to any physical sensor of the same type. Physical characteristics such as biases, scale factors, and axis-misalignments are irrelevant since the result after these effects is still a list of values representing the sensor measurements. As long as those values remain within the randomly generated training domain, the network will perform well. In Section 8, the same autoencoder is used on three different quadcopters each with different physical IMUs.

Concretely, using IMU preintegration as an example, the RMI itself must be communicated, which requires 10 floating-point numbers. However, the covariance matrix is 15 × 15, which would require communication of an additional 120 floating-point numbers to represent one of its triangular halves. With the proposed autoencoder, these 120 numbers are replaced with an encoding consisting of one number, thus dramatically reducing the communication cost. Figure 7 shows how the required communication rate varies with the duration between two successive communications between two arbitrary robots. The naive solution without preintegration requires sharing all input measurements that have occurred during that period, whereas preintegration yields a constant message size. The proposed method can be applied to all problems discussed in the paper, where the networks must be trained for each problem.

Figure 7.

RMSE for the ground robot simulation. There are four blue lines for the four robots running the proposed algorithm, and four visibly coincident red lines for the naive algorithm. Left: State fusion occurring at 10 Hz. Right: State fusion occurring at 1 Hz.

7. Simulation with ground robots

The proposed algorithm is tested in a simulation with ground robots, shown in the center of Figure 1. Each robot estimates their own pose and their neighbors’ poses relative to a world frame. Denoting the pose of Robot i relative to the world frame w as T _wi ∈ SE(2), the state of an arbitrary robot is given by

${\mathcal{X}}_{i_k}=\left({\mathbf{T}}_{w{i}_k}^{\left[i\right]},{\mathbf{T}}_{w{j}_k}^{\left[i\right]},\dots \right),\hspace{0.17em}j\in {\mathcal{N}}_i,$

where, again, the (·)^[i] superscript indicates Robot i’s estimate or “instance” of that physical quantity. Each robot collects wheel odometry at 100 Hz, providing ${\mathbf{u}}_{i_k}={\left[{\omega }_{i_k}{v}_{i_k}0\right]}^{\mathrm{T}}$ as input measurements, where ${\omega }_{i_k}$ is Robot i’s angular velocity and $v_{i_k}$ is its forward velocity in its own body frame. The pose kinematics for any single robot along with its preintegration are shown in Example 2. When Robot i receives an input measurement, it updates the part of its state corresponding to its own pose to create

${\mathfrak{X}}_{i_k}=\left({\mathbf{T}}_{i{w}_{k-1}}^{\left[i\right]}\mathrm{Exp}\left(\mathrm{\Delta }t{\mathbf{u}}_{i_{k-1}}\right),{\mathbf{T}}_{w{j}_{k-1}}^{\left[i\right]},\dots \right).$

The neighbor poses ${\mathbf{T}}_{w{j}_{k-1}}$ are now out of date, as neighboring odometry information is not yet accessible to Robot i, and this partially out-of-date state is non-physical and given the symbol ${\mathfrak{X}}_{i_k}$ . Every robot computes their own RMIs from wheel odometry using the equations from Example 2. When a neighbor’s RMI $\mathrm{\Delta }{\mathbf{T}}_{j_{pq}}$ is received at some later time step k = q, the state is updated with

${\mathcal{X}}_{i_q}=\left({\mathbf{T}}_{i{w}_q}^{\left[i\right]},{\mathbf{T}}_{w{j}_p}^{\left[i\right]}\mathrm{\Delta }{\mathbf{T}}_{j_{pq}},\dots \right)$

where p represents the time step index of the last time a neighbor RMI was received.

Each robot also collects range measurements to its neighbors at 10 Hz, with the connectivity graph shown in Figure 1 (middle). Only two robots collect relative position measurements to known landmarks at 10 Hz. At an arbitrary separate frequency, each robot sends its current state and covariance to its neighbors, allowing the neighbors to compute pseudomeasurements of the form

${\mathbf{c}}_{ij}\left({\mathcal{X}}_i,{\mathcal{X}}_j\right)=\left[\begin{array}{c}\mathrm{Log}\left({\mathbf{T}}_{wi}^{{\left[i\right]}^{-1}}{\mathbf{T}}_{wi}^{\left[j\right]}\right)\\ \mathrm{Log}\left({\mathbf{T}}_{wj}^{{\left[i\right]}^{-1}}{\mathbf{T}}_{wj}^{\left[j\right]}\right)\\ \mathrm{Log}\left({\mathbf{T}}_{w\mathcal{\ell }}^{{\left[i\right]}^{-1}}{\mathbf{T}}_{w\mathcal{\ell }}^{\left[j\right]}\right)\\ ⋮\end{array}\right],\mathcal{\ell }\in {\mathcal{N}}_i\cap {\mathcal{N}}_j.$

A simulation is performed with four robots each executing Algorithm 2, with root-mean-squared error (RMSE) shown in Figure 8. The initial states are initialized to ground truth with some random error with covariance ${\check{\mathbf{P}}}_{i_0}=0.1^2\cdot \mathbf{1}$ . The position and range measurements have Gaussian noise with 0.3 m and 0.1 m of standard deviation, respectively, and the pseudomeasurement covariance is Ψ = 0. The results show that all four robots’ estimation errors successfully stabilize and remain low, despite only two robots having sufficient sensors to make their states observable. Similar to the Toy Problem, the naive solution is implemented which is identical to the proposed solution, but does not do a CI step. A 50-trial Monte Carlo simulation was also performed with the resulting average NEES plotted in Figure 9. The naive solution has significantly higher error than both the centralized or proposed algorithms and is so overconfident that it cannot be plotted within reasonable axis limits in Figure 9. Although, in theory, the proposed algorithm should have lower NEES values than the centralized solution, Figure 9 shows comparable values. This is suspected to be due to the high degree of nonlinearity of the problem.

Figure 8.

50-trial NEES plot for the ground robot simulation for the proposed versus centralized solution, with the multiple blue lines each representing a robot. The naive solution without CI is far outside the plot. The red line represents the expected NEES value.

Figure 9.

Message size in bytes required to share odometry, as a function of the time period between communications between two robots. The red line “without preintegration” naively transmits all input measurements that occurred within the time period. Preintegration maintains a constant message size while providing identical information.

As seen in Figure 8, if state fusion is done at a sufficiently high-frequency, performance is even comparable to the centralized estimator, but this will incur a larger communication and computation requirement as discussed in Section 5.

8. Simulation and experiments with quadcopters

To demonstrate the flexibility of the proposed framework, consider a new problem involving quadcopters. The kinematic state of each quadcopter is modeled using extended pose matrices T ∈ SE₂(3) (Brossard et al. (2022)). Each robot estimates both their absolute pose relative to the world frame T _wi ∈ SE₂(3), their own IMU bias b _i , as well as the relative poses of their neighbors T _ij ∈ SE₂(3) and their IMU bias b _j . The full state of Robot i is then given by

${\mathcal{X}}_i=\left({\mathbf{T}}_{wi},{\mathbf{b}}_i^{\left[i\right]},{\mathbf{T}}_{ij},{\mathbf{b}}_j^{\left[i\right]},\dots \right),\hspace{0.17em}j\in {\mathcal{N}}_i.$

The pose of Robot j relative to Robot i T _ij has kinematics involving the IMU measurements of both robots, and are given in discrete time by

${\mathbf{T}}_{i{j}_k}={\mathbf{U}}_{i_{k-1}}^{-1}{\mathbf{T}}_{i{j}_{k-1}}{\mathbf{U}}_{j_{k-1}},$

where ${\mathbf{U}}_{j_{k-1}}$ has an identical definition as in (20), but computed from Robot j’s IMU measurements. When Robot i receives input measurements from its own IMU u _k , it predicts the part of its own state corresponding to its own pose, and additionally performs a partial prediction on the relative poses with

${\mathfrak{X}}_{i_k}=\left({\mathbf{G}}_{k-1}{\mathbf{T}}_{w{i}_{k-1}}{\mathbf{U}}_{i_{k-1}},{\mathbf{b}}_{i_{k-1}}^{\left[i\right]},{\mathbf{U}}_{i_{k-1}}^{-1}{\mathbf{T}}_{i{j}_{k-1}},{\mathbf{b}}_{j_{k-1}}^{\left[i\right]},\dots \right).$

The terms which are the partially predicted neighbor poses, are a strange, non-physical intermediate state. Only when the neighbor’s RMI $\mathrm{\Delta }{\mathbf{U}}_{j_{pq}}$ is received do the neighbor poses regain meaning with However, since biases are also being estimated in this problem, Robot i must first correct the neighbor’s raw RMIs $\mathrm{\Delta }{\mathbf{U}}_{j_{pq}}\left({\mathbf{u}}_{j_{p:q-1}},\mathbf{0}\right)$ using its estimate of the neighbor’s IMU bias, as described in Section 6.2. That is,

$\mathrm{\Delta }{\mathbf{U}}_{j_{pq}}\approx \mathrm{\Delta }{\mathbf{U}}_{j_{pq}}\left({\mathbf{u}}_{j_{p:q-1}},\mathbf{0}\right)\oplus {\mathbf{B}}_{j_{pq}}{\mathbf{b}}_{j_q}^{\left[i\right]},$

leading to the full state update given by

Finally, the pseudomeasurements chosen for this problem are

${\mathbf{c}}_{ij}\left({\mathcal{X}}_i,{\mathcal{X}}_j\right)=\left[\begin{array}{c}\mathrm{Log}\left({\mathbf{T}}_{wi}{\mathbf{T}}_{ij}{\mathbf{T}}_{wj}^{-1}\right)\\ \mathrm{Log}\left({\mathbf{T}}_{ij}{\mathbf{T}}_{ji}\right)\\ {\mathbf{b}}_i^{\left[i\right]}-{\mathbf{b}}_i^{\left[j\right]}\\ {\mathbf{b}}_j^{\left[i\right]}-{\mathbf{b}}_j^{\left[j\right]}\\ \mathrm{Log}\left({\mathbf{T}}_{ij}{\mathbf{T}}_{j\mathcal{\ell }}{\mathbf{T}}_{i\mathcal{\ell }}^{-1}\right)\\ ⋮\end{array}\right],\mathcal{\ell }\in {\mathcal{N}}_i\hspace{0.17em}\cap \hspace{0.17em}{\mathcal{N}}_j.$

with corresponding covariance Ψ = 0.

8.1. Hardware setup

The hardware setup in these experiments can be seen in Figure 10. Three Uvify IFO-S quadcopters are used that each possess an IMU at 200 Hz, a 1D LIDAR height sensor at 30 Hz, and magnetometers at 30 Hz. Additionally, two ultra-wideband (UWB) transceivers are installed on the quadcopter legs, producing inter-robot distance measurements at 90 Hz for each robot. As shown by Shalaby et al. (2021a), installing multiple UWB tags per robot results in relative position observability. The UWB transceivers are custom-printed modules that use the DW1000 UWB transceiver. The firmware for these modules has been written in C, implementing a double-sided two-way-ranging protocol with details described by Shalaby et al. (2022). Shalaby et al. (2022) also describe the power-based bias calibration and noise characterization procedure used in these experiments. Since all transceivers operate on the same frequency in these experiments, only one can transmit at a time to avoid interference. A decentralized scheduler is therefore implemented that continuously cycles through all transceiver pairs one at a time, obtaining range measurements and potentially transmitting other useful data. In these experiments, the communication graph is complete with all quadcopters capable of communicating with each other. Preintegrated RMIs are shared whenever a UWB measurement occurs, and state sharing occurs at a separate frequency of 10 Hz.

Figure 10.

Simulated estimation error of Robot 3’s estimate of its own kinematic state and IMU biases. The estimate and corresponding bounds with the proposed algorithm are shown in blue, with the centralized estimate overlayed in dark gray. For attitude, the x − y − z components represent roll-pitch-yaw errors, respectively. Note that Robot 3 does not have position measurements, and therefore cannot observe the states shown in this plot without information sharing. The naive solution has rapidly diverging error and is not plotted.

A Vicon motion capture system is used to collect ground truth, from which synthesized absolute position measurements with a standard deviation of 0.3 m are generated for Robots 1 and 2 only. Robot 3 does not receive absolutely position measurements, nor any magnetometer measurements, and therefore has no absolute pose information available without communication with the other two robots. Example trajectories for some of the experimental trials are shown in Figure 11.

Figure 11.

Top: Three quadcopters in flight under a motion capture system. Bottom left: custom UWB module. Bottom right: a close up of the Uvify IFO-S quadcopter, fitted with a UWB module seen on the left leg, as well as on the opposite leg.

8.2. Simulation results

The algorithm is first tested with simulated versions of the described quadcopters, and the estimation results for Robot 3’s absolute pose and bias are shown in Figure 12. Although there are many other states associated with the simulation, these states are the most interesting as they are the ones that are unobservable without incorporation of the pseudomeasurements. Figure 12 shows that Robot 3 is capable of estimating its own absolute pose and bias, using information from sensors located on Robots 1 and 2. Furthermore, the errors remain within the 3-sigma confidence bounds, even with the first-order RMI bias correction, indicating statistical consistency.

Figure 12.

Examples of the various trajectories flown in the experimental trials, where each color represents a different quadcopter.

Figure 13 shows the positioning RMSE for varying frequency at which state information between robots is shared. At lower frequencies, Robot 3’s estimate has more time to drift between communications, and hence, there is higher error. For this problem, roughly the same estimation performance is achieved for state sharing of 20 Hz and above, with 10 Hz being a compromising value providing a trade-off between accuracy and communication cost.

Figure 13.

Average self-positioning RMSE with varying communication rate for the simulated version of the quadcopter problem. Robot 3 does not receive position measurements, and hence is reliant on the other robots to have an observable state.

8.3. Experimental results

Multiple experimental runs are performed on different days, with the absolute positioning results for each robot viewable in Table 1. In some cases, the proposed algorithm even outperforms the centralized solution, which is theoretically optimal. However, the real world contains many unmodelled sources of error, such as frame misalignments, timestamping errors, vibrations, and UWB ranging outliers. These effects may break the assumptions that the optimality of the centralized estimator relies on. Even after tuning covariances to obtain the best centralized performance, it appears that the results benefit from the covariance inflation resulting from CI. For both estimators, the IMU is calibrated to compensate for large biases and scaling factors. The normalized-innovation-squared test (Bar-Shalom et al. (2001)) is also used to reject UWB outliers in both estimators,

Table 1.

Self-Positioning RMSE (m) From Experimental Trials.

Trial #	Centralized			Proposed			Error Reduction
	Robot 1	Robot 2	Robot 3	Robot 1	Robot 2	Robot 3	Robot 1	Robot 2	Robot 3
1	0.43	0.49	0.55	0.22	0.22	0.61	−48%	−54%	10%
2	0.18	0.26	0.34	0.16	0.18	0.40	−13%	−28%	16%
3	0.17	0.24	0.68	0.16	0.17	0.45	−10%	−28%	−33%
4	0.20	0.25	0.31	0.26	0.28	0.48	32%	−13%	55%

A plot of RMSE versus time for Robot 3’s absolute states can be seen in Figure 14, which are states that are unobservable from Robot 3’s own measurements. Again, Figure 14 shows that error magnitudes lie in similar ranges for both the centralized and decentralized estimators. Figure 15 compares two decentralized estimator runs, with one using the mean-assisted autoencoder from Section 6.3. As desired, the lines are identical, and the plot shows that the estimate is unaffected by the autoencoding. This means that the autoencoder is highly effective at compressing the covariance matrix with minimal reconstruction error.

Figure 14.

Position, velocity, and yaw RMSE for Robot three from one of the experimental trials. Since Robot 3 has no position measurements, these quantities are unobservable without the fusion of pseudomeasurements.

Figure 15.

The effect of preintegrated covariance autoencoding, as described in Section 6.3, on the position estimate of Robot 1. The two lines are almost identical, showing that the proposed autoencoder induces minimal error on the estimate. All other states have similar plots.

9. Conclusion

This paper presents a general-purpose algorithm for decentralized state estimation in robotics. The algorithm is the result of a new way to formulate the decentralized state estimation problem, specifically with the assistance of pseudomeasurements that allow the definition of arbitrary nonlinear relationships between robot states. For problems involving relative measurements, a communication-efficient approach is proposed for preintegratable process models, as defined by (18), where state-change information is shared in the form of relative motion increments. The algorithm is tested on three different problems, each involving a variety of state definitions, process models, and measurements. In all of the presented problems, robots only need to share their states, RMIs, and corresponding covariances, which ultimately results in average transmission rates per robot of 0.2 kB/s for the toy problem, 4.5 kB/s for the ground robots, and 53.2 kB/s for the quadcopters.

Thanks to covariance intersection, the algorithm is appropriate for arbitrary graphs, and does not require any bookkeeping, growing memory, buffering of measurements, or reprocessing of data. At the same time, the approximation made by covariance intersection makes the proposed method suboptimal, as it is well-known to be overly conservative. Nevertheless, in the specific problems shown in this paper, the results using CI have been satisfactory provided that the fusion frequency is high enough, and the communication graph is not too sparse. It is also worth mentioning that the proposed algorithm still assumes that process model inputs, whether in raw or preintegrated form, have noise that is uncorrelated with the robot states, just like the sensor measurement noise. These assumptions must hold for a consistent estimator. In this paper, it is only correlations between different robots’ states that are mitigated by covariance intersection.

One limitation of this proposed approach is that the communication cost grows quadratically with the state size, since the state covariance is also shared. While this is not an issue for small state sizes, such as those representing 3D poses, it could become a problem for states involving multiple time steps or a very large number of robots. Furthermore since this paper allows for variable state definitions between robots, the state definition itself for each robot may need to be communicated, or established a priori. The initial value of the state is also assumed to be known through an arbitrary initialization procedure.

Future work can consider improving the approximation made by CI and compressing the covariance matrix associated with the state. Also, using the proposed MAP approach with pseudomeasurements, it should be possible to derive decentralized batch and sliding-window estimators, often termed smoothers, since these algorithms also originate from the MAP problem.