# Decentralized State Estimation Algorithm of Centralized Equivalent Precision for Formation Flying Spacecrafts Based on Junction Tree

Authors: Mengyuan Dai; Hua Mu; Meiping Wu; Zhiwen Xian
Journal: International Journal of Aerospace Engineering
 (2015)

---

## Abstract

As centralized state estimation algorithms for formation flying spacecraft would suffer from high computational burdens when the scale of the formation increases, it is necessary to develop decentralized algorithms. To the state of the art, most decentralized algorithms for formation flying are derived from centralized EKF by simplification and decoupling, rendering suboptimal estimations. In this paper, typical decentralized state estimation algorithms are reviewed, and a new scheme for decentralized algorithms is proposed. In the new solution, the system is modeled as a dynamic Bayesian network (DBN). A probabilistic graphical method named junction tree (JT) is used to analyze the hidden distributed structure of the DBNs. Inference on JT is a decentralized form of centralized Bayesian estimation (BE), which is a modularized three-step procedure of receiving messages, collecting evidences, and generating messages. As KF is a special case of BE, the new solution based on JT is equivalent in precision to centralized KF in theory. A cooperative navigation example of a three-satellite formation is used to test the decentralized algorithms. Simulation results indicate that JT has the best precision among all current decentralized algorithms.


---

## Body

## 1. Introduction

Formation flying has been recognized as an enabling technology for many future space missions [1], for example, the synthetic aperture radar (SAR) [2], and the stellar interferometry [3]. Compared with conventional space missions based on a single spacecraft, formation flying has several advantages: lower price, robustness after spacecraft failures, modularity, and good flexibility for reconfigurations [4]. The concept of formation flying is driven by the idea that a fleet of small spacecrafts could form a virtual spacecraft of unlimited scale which may serve as a substitution of a monolithic spacecraft [4]. As the size of the formation becomes larger, the guidance, navigation, and control (GNC) task becomes more and more onerous due to the increased size of state and measurement data [5]. Thus it becomes a necessity to develop decentralized algorithms to balance the computational loads among spacecrafts.This paper takes the navigation task as an example to study the decentralized state estimation algorithms. There are generally two kinds of sensors used for formation navigation: the proprioceptive sensors and the exteroceptive sensors [6]. The first kind of sensors is used for single platform measurements to estimate the state of local spacecraft, for example, GPS receivers [4, 5, 7]. The second kind of sensors is used for interformation measurements to estimate the states of local spacecraft and its neighbors, for example, optical devices [8, 9] and radio frequency (RF) devices [10, 11]. With interformation measurements and communications, the navigation resources at different spacecrafts may be shared among the formation. This kind of navigation is also known as cooperative navigation [12]. As a typical navigation application for formation flying spacecrafts, this paper studies formation flying missions using RF range augmented GPS sensors, which are introduced by Park [4] and Ferguson and How [5]. Though the sensors are specially appointed in this paper, the technique used in designing decentralized algorithms may inspire other cooperative navigation missions.To the state of the art, almost all the decentralized estimation algorithms for formation flying are based on distributing centralized EKF through simplifying the measurement models or decoupling the interformation measurements through approximations. In fact, as interformation measurements make the states of spacecrafts strongly coupled, the simplification and decoupling will result in suboptimal estimations, which means that the precision of these algorithms is inferior to EKF. From the perspective of probabilistic graphical models in artificial intelligence, Lauritzen introduced a method for modeling and local computing of exact mean and covariance in dynamic systems based on JT [13]. The most remarkable advantage of the JT algorithm to other decentralized algorithms is that the JT algorithm has equivalent precision to centralized EKF [14]. Thus it is possible to develop a new decentralized cooperative navigation algorithm for formation flying spacecrafts.This paper is arranged in 5 sections. In Section2, the model for cooperative navigation is introduced. A brief introduction of existing decentralized solutions is provided in Section 3. And Section 4 introduces the new decentralized solution as a series of local computations and messages passing on JT. Simulation results in Section 5 show that the new solution has centralized equivalent precision which is better than the decentralized algorithms reviewed.


## 2. Model of Cooperative Navigation for Formation Flying

This paper takes the formation flying mission based on RF range augmented GPS as an example. The GPS pseudorange is used for position estimate at local spacecraft, while the RF range devices are used for relative motion measurements between spacecrafts. As the dynamic model and the measurement model are both nonlinear, some partial derivatives are needed for linearization. The partials used in the numerical simulation of this research are listed in AppendixA.


### 2.1. Dynamic Model

The state of each spacecraft is defined as the Keplerian orbital elements. Letx
k
i be the Keplerian orbital elements of satellite i at time step k, x
=
a
,
e
,
i
′
,
Ω
,
ω
,
M
T, where a is the semimajor axis, e is the eccentricity, i
′ is the inclination angle, Ω is the right ascension of ascending node, ω is the argument of perigee, and M is the mean anomaly at reference time. Without considering system inputs such as orbit control, the dynamic model of satellite i is given by(1)
x
k
+
1
i
=
f
x
k
i
+
w
k
i
,where w
k
i is the noise and w
k
i
~
N
0
,
Q
k
i, N
0
,
Q
k
i is a Gaussian distribution of mean 0 and variance Q
k
i.Linearize (1) at x
^
k
i using first order Tyler expansion, and we have(2)
x
k
+
1
i
=
f
x
^
k
i
+
F
k
i
·
x
k
i
-
x
^
k
i
+
w
k
i
,where F
k
i is the state transition matrix and(3)
F
k
i
=
∂
f
∂
x
k
i
x
k
i
=
x
^
k
i
.


### 2.2. Measurement Model

There are two kinds of measurements in the formation flying example: single platform measurement based on GPS pseudorange and interformation measurement based on RF range.


#### 2.2.1. Single Platform Measurement

For simplicity, the time delay for signal transmission and errors such as hardware delay, ionosphere delay, multipath, relativistic effect, and clock error are not considered. Letr
k
+
1
s be the position of observed GPS satellite given by the broadcast ephemeris, and let r
k
+
1
i be the position of spacecraft i at time step k
+
1. The single platform measurement based on GPS pseudorange, z
k
+
1
i, is given by(4)
z
k
+
1
i
=
r
k
+
1
i
-
r
k
+
1
s
+
v
k
+
1
i
,where r means the length of vector r and v
k
+
1
i is the measurement noise, v
k
+
1
i
~
N
0
,
R
k
+
1
i.In general, the single platform measurement function may be written as(5)
z
k
+
1
i
=
h
1
x
k
+
1
i
+
v
k
+
1
i
.Linearize (5) at x
^
k
+
1
i using first order Tyler expansion, and we have(6)
z
k
+
1
i
=
h
1
x
^
k
+
1
i
+
H
1
,
k
+
1
i
·
x
k
+
1
i
-
x
^
k
+
1
i
+
v
k
+
1
i
,where H
1
,
k
+
1
i is the measurement matrix and (7)
H
1
,
k
+
1
i
=
∂
h
1
∂
x
k
+
1
i
x
k
+
1
i
=
x
^
k
+
1
i
.


#### 2.2.2. Interformation Measurement

Similar to single platform measurement, letr
k
+
1
i and r
k
+
1
j be the position of spacecraft i and j at time step k
+
1, respectively. The interformation measurement between spacecrafts i and j, z
k
+
1
i
j, is given by(8)
z
k
+
1
i
j
=
r
k
+
1
i
-
r
k
+
1
j
+
v
k
+
1
i
j
,where v
k
+
1
i
j is the measurement noise, v
k
+
1
i
j
~
N
0
,
R
k
+
1
i
j.In general, the single platform measurement function may be written as(9)
z
k
+
1
i
j
=
h
2
x
k
+
1
i
,
x
k
+
1
j
+
v
k
+
1
i
j
.Linearize (9) at x
^
k
+
1
i and x
^
k
+
1
j using first order Tyler expansion, and we have(10)
z
k
+
1
i
j
=
h
2
x
^
k
+
1
i
,
x
^
k
+
1
j
+
H
2
,
k
+
1
i
·
x
k
+
1
i
-
x
^
k
+
1
i
+
H
2
,
k
+
1
j
·
x
k
+
1
j
-
x
^
k
+
1
j
+
v
k
+
1
i
j
,where H
2
,
k
+
1
i and H
2
,
k
+
1
j are the measurement matrix and (11)
H
2
,
k
+
1
i
=
∂
h
2
∂
x
k
+
1
i
x
k
+
1
i
=
x
^
k
+
1
i
,
H
2
,
k
+
1
j
=
∂
h
2
∂
x
k
+
1
j
x
k
+
1
j
=
x
^
k
+
1
j
.


## 2.1. Dynamic Model

The state of each spacecraft is defined as the Keplerian orbital elements. Letx
k
i be the Keplerian orbital elements of satellite i at time step k, x
=
a
,
e
,
i
′
,
Ω
,
ω
,
M
T, where a is the semimajor axis, e is the eccentricity, i
′ is the inclination angle, Ω is the right ascension of ascending node, ω is the argument of perigee, and M is the mean anomaly at reference time. Without considering system inputs such as orbit control, the dynamic model of satellite i is given by(1)
x
k
+
1
i
=
f
x
k
i
+
w
k
i
,where w
k
i is the noise and w
k
i
~
N
0
,
Q
k
i, N
0
,
Q
k
i is a Gaussian distribution of mean 0 and variance Q
k
i.Linearize (1) at x
^
k
i using first order Tyler expansion, and we have(2)
x
k
+
1
i
=
f
x
^
k
i
+
F
k
i
·
x
k
i
-
x
^
k
i
+
w
k
i
,where F
k
i is the state transition matrix and(3)
F
k
i
=
∂
f
∂
x
k
i
x
k
i
=
x
^
k
i
.


## 2.2. Measurement Model

There are two kinds of measurements in the formation flying example: single platform measurement based on GPS pseudorange and interformation measurement based on RF range.


### 2.2.1. Single Platform Measurement

For simplicity, the time delay for signal transmission and errors such as hardware delay, ionosphere delay, multipath, relativistic effect, and clock error are not considered. Letr
k
+
1
s be the position of observed GPS satellite given by the broadcast ephemeris, and let r
k
+
1
i be the position of spacecraft i at time step k
+
1. The single platform measurement based on GPS pseudorange, z
k
+
1
i, is given by(4)
z
k
+
1
i
=
r
k
+
1
i
-
r
k
+
1
s
+
v
k
+
1
i
,where r means the length of vector r and v
k
+
1
i is the measurement noise, v
k
+
1
i
~
N
0
,
R
k
+
1
i.In general, the single platform measurement function may be written as(5)
z
k
+
1
i
=
h
1
x
k
+
1
i
+
v
k
+
1
i
.Linearize (5) at x
^
k
+
1
i using first order Tyler expansion, and we have(6)
z
k
+
1
i
=
h
1
x
^
k
+
1
i
+
H
1
,
k
+
1
i
·
x
k
+
1
i
-
x
^
k
+
1
i
+
v
k
+
1
i
,where H
1
,
k
+
1
i is the measurement matrix and (7)
H
1
,
k
+
1
i
=
∂
h
1
∂
x
k
+
1
i
x
k
+
1
i
=
x
^
k
+
1
i
.


### 2.2.2. Interformation Measurement

Similar to single platform measurement, letr
k
+
1
i and r
k
+
1
j be the position of spacecraft i and j at time step k
+
1, respectively. The interformation measurement between spacecrafts i and j, z
k
+
1
i
j, is given by(8)
z
k
+
1
i
j
=
r
k
+
1
i
-
r
k
+
1
j
+
v
k
+
1
i
j
,where v
k
+
1
i
j is the measurement noise, v
k
+
1
i
j
~
N
0
,
R
k
+
1
i
j.In general, the single platform measurement function may be written as(9)
z
k
+
1
i
j
=
h
2
x
k
+
1
i
,
x
k
+
1
j
+
v
k
+
1
i
j
.Linearize (9) at x
^
k
+
1
i and x
^
k
+
1
j using first order Tyler expansion, and we have(10)
z
k
+
1
i
j
=
h
2
x
^
k
+
1
i
,
x
^
k
+
1
j
+
H
2
,
k
+
1
i
·
x
k
+
1
i
-
x
^
k
+
1
i
+
H
2
,
k
+
1
j
·
x
k
+
1
j
-
x
^
k
+
1
j
+
v
k
+
1
i
j
,where H
2
,
k
+
1
i and H
2
,
k
+
1
j are the measurement matrix and (11)
H
2
,
k
+
1
i
=
∂
h
2
∂
x
k
+
1
i
x
k
+
1
i
=
x
^
k
+
1
i
,
H
2
,
k
+
1
j
=
∂
h
2
∂
x
k
+
1
j
x
k
+
1
j
=
x
^
k
+
1
j
.


## 2.2.1. Single Platform Measurement

For simplicity, the time delay for signal transmission and errors such as hardware delay, ionosphere delay, multipath, relativistic effect, and clock error are not considered. Letr
k
+
1
s be the position of observed GPS satellite given by the broadcast ephemeris, and let r
k
+
1
i be the position of spacecraft i at time step k
+
1. The single platform measurement based on GPS pseudorange, z
k
+
1
i, is given by(4)
z
k
+
1
i
=
r
k
+
1
i
-
r
k
+
1
s
+
v
k
+
1
i
,where r means the length of vector r and v
k
+
1
i is the measurement noise, v
k
+
1
i
~
N
0
,
R
k
+
1
i.In general, the single platform measurement function may be written as(5)
z
k
+
1
i
=
h
1
x
k
+
1
i
+
v
k
+
1
i
.Linearize (5) at x
^
k
+
1
i using first order Tyler expansion, and we have(6)
z
k
+
1
i
=
h
1
x
^
k
+
1
i
+
H
1
,
k
+
1
i
·
x
k
+
1
i
-
x
^
k
+
1
i
+
v
k
+
1
i
,where H
1
,
k
+
1
i is the measurement matrix and (7)
H
1
,
k
+
1
i
=
∂
h
1
∂
x
k
+
1
i
x
k
+
1
i
=
x
^
k
+
1
i
.


## 2.2.2. Interformation Measurement

Similar to single platform measurement, letr
k
+
1
i and r
k
+
1
j be the position of spacecraft i and j at time step k
+
1, respectively. The interformation measurement between spacecrafts i and j, z
k
+
1
i
j, is given by(8)
z
k
+
1
i
j
=
r
k
+
1
i
-
r
k
+
1
j
+
v
k
+
1
i
j
,where v
k
+
1
i
j is the measurement noise, v
k
+
1
i
j
~
N
0
,
R
k
+
1
i
j.In general, the single platform measurement function may be written as(9)
z
k
+
1
i
j
=
h
2
x
k
+
1
i
,
x
k
+
1
j
+
v
k
+
1
i
j
.Linearize (9) at x
^
k
+
1
i and x
^
k
+
1
j using first order Tyler expansion, and we have(10)
z
k
+
1
i
j
=
h
2
x
^
k
+
1
i
,
x
^
k
+
1
j
+
H
2
,
k
+
1
i
·
x
k
+
1
i
-
x
^
k
+
1
i
+
H
2
,
k
+
1
j
·
x
k
+
1
j
-
x
^
k
+
1
j
+
v
k
+
1
i
j
,where H
2
,
k
+
1
i and H
2
,
k
+
1
j are the measurement matrix and (11)
H
2
,
k
+
1
i
=
∂
h
2
∂
x
k
+
1
i
x
k
+
1
i
=
x
^
k
+
1
i
,
H
2
,
k
+
1
j
=
∂
h
2
∂
x
k
+
1
j
x
k
+
1
j
=
x
^
k
+
1
j
.


## 3. Review of Decentralized Algorithms

To the best knowledge of the authors, most current decentralize solutions are based on centralized EKF. Generally, there are two types of decentralized algorithms for formation flying spacecraft: the full order extended Kalman filter (FOEKF) and the reduced order extended Kalman filter (ROEKF) [5]. In FOEKF, each spacecraft runs a full order estimator which estimates the full states of the system, but only based on local measurements related to the very spacecraft. In ROEKF, each spacecraft runs a reduced order estimator which only processes the state of itself. The reduction of the order of the filter is based on decoupling the measurement models. As FOEKF and ROEKF all require simplifications and approximations of the EKF model, they all suffer from suboptimal precision which is inferior to centralized EKF. To compensate the precision loss, a direct method is to use iterations, which leads to the iterated cascade extended Kalman filter (ICEKF) [5].


### 3.1. Centralized EKF

Consider a formation ofn spacecrafts. Let X
k be the state of the whole formation at time step k, X
k
=
x
k
1
T
x
k
2
T
⋯
x
k
n
T
T. Without considering the system input, the linearized model of cooperative navigation could be rewritten as follows.Dynamic Model.
(12)
X
k
+
1
=
Φ
k
X
k
+
W
k
,
W
k
~
N
0
,
Q
k
.Measurement Model.
(13)
Z
k
+
1
=
H
k
+
1
X
k
+
1
+
V
k
+
1
,
V
k
+
1
~
N
0
,
R
k
+
1
,where Φ is the state transition matrix, W is the noise of the dynamic model, Z is the vector of measurements, H is the measurement matrix, and V is the noise of measurements. The procedure of EKF could be described as a sequence of time update and measurement update.Time Update.
(14)
X
^
k
+
1
-
=
Φ
k
X
^
k
+
,
(15)
P
k
+
1
-
=
Φ
k
P
k
+
Φ
k
T
+
Q
k
.Measurement Update.
(16)
K
k
+
1
=
P
k
+
1
-
H
k
+
1
T
H
k
+
1
P
k
+
1
-
H
k
+
1
T
+
R
k
+
1
-
1
,
(17)
X
^
k
+
1
+
=
X
^
k
+
1
-
+
K
k
+
1
Z
k
+
1
-
H
k
+
1
X
^
k
+
1
-
,
(18)
P
k
+
1
+
=
I
-
K
k
+
1
H
k
+
1
P
k
+
1
-
,where mean X
^
k
+ and covariance P
k
+ are the best estimates of system state at time step k.


### 3.2. FOEKF

Considering a local FOEKF estimator on spacecrafti, it only processes the measurements related to spacecraft i, which means that the measurement Z
k
+
1
i, measurement matrix H
k
+
1
i, and measurement noise V
k
+
1
i are subsets of Z
k
+
1, H
k
+
1, and V
k
+
1 in EKF, respectively. The measurement update of FOEKF at spacecraft i is as follows:(19)
K
k
+
1
=
P
k
+
1
-
H
k
+
1
i
T
H
k
+
1
i
P
k
+
1
-
H
k
+
1
i
T
+
R
k
+
1
i
-
1
,
X
^
k
+
1
+
=
X
^
k
+
1
-
+
K
k
+
1
Z
k
+
1
i
-
H
k
+
1
i
X
^
k
+
1
-
,
P
k
+
1
+
=
I
-
K
k
+
1
H
k
+
1
i
P
k
+
1
-
.Given that the RF range and Doppler measurements are both nonlinear (the linearization requires the states of other spacecrafts), spacecrafts doing relative measurements must communicate with each other exchanging information of mean and covariance of its own state. Because only part of the measurements is considered at each spacecraft, FOEKF is suboptimal. Meanwhile, as FOEKF runs a full order estimator at each spacecraft, the computational load may not be well balanced.


### 3.3. ROEKF

#### 3.3.1. Original ROEKF

ROEKF is a further simplification of FOEKF. Taking spacecrafti for example, suppose the states of other spacecrafts are known, and the measurement model between spacecrafts i and j could be simplified as(20)
z
k
+
1
i
j
=
H
2
,
k
+
1
i
x
k
+
1
i
+
H
2
,
k
+
1
j
x
k
+
1
j
+
v
k
+
1
i
j
,
v
k
+
1
i
j
~
N
0
,
R
k
+
1
i
j
.Notice thatx
k
+
1
j is known. Let z
-
k
+
1
i
j
=
z
k
+
1
i
j
-
H
2
,
k
+
1
j
x
k
+
1
j; the ROEKF measurement model could be rewrite as(21)
Z
-
k
+
1
i
=
H
-
k
+
1
i
x
k
+
1
i
+
V
k
+
1
i
,
V
k
+
1
i
~
N
0
,
R
k
+
1
i
,where Z
-
k
+
1
i is the set of all the measurements related with spacecraft i and H
-
k
+
1
i is the measurement matrix related with Z
-
k
+
1
i.Thus the local estimator at spacecrafti only needs to estimate its own state, x
^
k
i and P
k
i. The measurement update of FOEKF at spacecraft i is as follows:(22)
K
k
+
1
=
P
k
+
1
i
-
H
-
k
+
1
i
T
H
-
k
+
1
i
P
k
+
1
i
-
H
-
k
+
1
i
T
+
R
k
+
1
i
-
1
,
x
^
k
+
1
i
+
=
x
^
k
+
1
i
-
+
K
k
+
1
Z
-
k
+
1
i
-
H
-
k
+
1
i
x
^
k
+
1
i
-
,
P
k
+
1
i
+
=
I
-
K
k
+
1
H
-
k
+
1
i
P
k
+
1
i
-
.Compared with FOEKF, the computational load is reduced remarkably in ROEKF, because ROEKF only updates the state of local spacecraft. However, as only estimates instead of real states of other spacecrafts could be obtained, the errors of state estimation of other spacecrafts are not considered in ROEKF, rendering a further precision loss than FOEKF. We could use a proper increase to the measurement noise, theR matrix, to absorb such errors. This leads to the increased R extended Kalman filter (IREKF) [5].


#### 3.3.2. IREKF

Considering the measurementz
k
+
1
i
j, the measurement noise caused by using estimates of the state of spacecraft j could be estimated by the following equation:(23)
E
z
k
+
1
i
j
-
H
2
,
k
+
1
i
x
k
+
1
i
z
k
+
1
i
j
-
H
2
,
k
+
1
i
x
k
+
1
i
T
=
E
H
2
,
k
+
1
j
x
k
+
1
j
+
v
k
+
1
i
j
H
2
,
k
+
1
j
x
k
+
1
j
+
v
k
+
1
i
j
T
=
E
H
2
,
k
+
1
j
x
k
+
1
j
H
2
,
k
+
1
j
x
k
+
1
j
T
+
H
2
,
k
+
1
j
x
k
+
1
j
v
k
+
1
i
j
T
+
v
k
+
1
i
j
H
2
,
k
+
1
j
x
k
+
1
j
T
+
v
k
+
1
i
j
v
k
+
1
i
j
T
=
E
H
2
,
k
+
1
j
x
k
+
1
j
H
2
,
k
+
1
j
x
k
+
1
j
T
+
E
v
k
+
1
i
j
v
k
+
1
i
j
T
=
H
2
,
k
+
1
j
P
k
+
1
j
+
H
2
,
k
+
1
j
T
+
R
k
+
1
i
j
.Notice that only estimates of the states of other spacecrafts could be obtained, which means that we could only useP
k
+
1
j
- to approximate P
k
+
1
j
+. Thus the increase to the measurement noise of R
i
j is given by H
2
,
k
+
1
j
P
k
+
1
j
+
H
2
,
k
+
1
j
T.


### 3.4. ICEKF

Iteration is a straight forward method to compensate the precision loss caused by approximations. ICEKF could serve as a general scheme to improve FOEKF and ROEKF. In ICEKF, the local estimation at each spacecraft is passed to other spacecrafts for the linearization of measurement functions through a cascade communication chain. Figure1 gives an example of ICEKF for a three-spacecraft formation compared with EKF. Although ICEKF reduces the computational load while obtaining a better precision than FOEKF or ROEKF, the introduction of iteration will result in an increase in the communicational load.Figure 1
ICEKF compared with centralized EKF [15].


## 3.1. Centralized EKF

Consider a formation ofn spacecrafts. Let X
k be the state of the whole formation at time step k, X
k
=
x
k
1
T
x
k
2
T
⋯
x
k
n
T
T. Without considering the system input, the linearized model of cooperative navigation could be rewritten as follows.Dynamic Model.
(12)
X
k
+
1
=
Φ
k
X
k
+
W
k
,
W
k
~
N
0
,
Q
k
.Measurement Model.
(13)
Z
k
+
1
=
H
k
+
1
X
k
+
1
+
V
k
+
1
,
V
k
+
1
~
N
0
,
R
k
+
1
,where Φ is the state transition matrix, W is the noise of the dynamic model, Z is the vector of measurements, H is the measurement matrix, and V is the noise of measurements. The procedure of EKF could be described as a sequence of time update and measurement update.Time Update.
(14)
X
^
k
+
1
-
=
Φ
k
X
^
k
+
,
(15)
P
k
+
1
-
=
Φ
k
P
k
+
Φ
k
T
+
Q
k
.Measurement Update.
(16)
K
k
+
1
=
P
k
+
1
-
H
k
+
1
T
H
k
+
1
P
k
+
1
-
H
k
+
1
T
+
R
k
+
1
-
1
,
(17)
X
^
k
+
1
+
=
X
^
k
+
1
-
+
K
k
+
1
Z
k
+
1
-
H
k
+
1
X
^
k
+
1
-
,
(18)
P
k
+
1
+
=
I
-
K
k
+
1
H
k
+
1
P
k
+
1
-
,where mean X
^
k
+ and covariance P
k
+ are the best estimates of system state at time step k.


## 3.2. FOEKF

Considering a local FOEKF estimator on spacecrafti, it only processes the measurements related to spacecraft i, which means that the measurement Z
k
+
1
i, measurement matrix H
k
+
1
i, and measurement noise V
k
+
1
i are subsets of Z
k
+
1, H
k
+
1, and V
k
+
1 in EKF, respectively. The measurement update of FOEKF at spacecraft i is as follows:(19)
K
k
+
1
=
P
k
+
1
-
H
k
+
1
i
T
H
k
+
1
i
P
k
+
1
-
H
k
+
1
i
T
+
R
k
+
1
i
-
1
,
X
^
k
+
1
+
=
X
^
k
+
1
-
+
K
k
+
1
Z
k
+
1
i
-
H
k
+
1
i
X
^
k
+
1
-
,
P
k
+
1
+
=
I
-
K
k
+
1
H
k
+
1
i
P
k
+
1
-
.Given that the RF range and Doppler measurements are both nonlinear (the linearization requires the states of other spacecrafts), spacecrafts doing relative measurements must communicate with each other exchanging information of mean and covariance of its own state. Because only part of the measurements is considered at each spacecraft, FOEKF is suboptimal. Meanwhile, as FOEKF runs a full order estimator at each spacecraft, the computational load may not be well balanced.


## 3.3. ROEKF

### 3.3.1. Original ROEKF

ROEKF is a further simplification of FOEKF. Taking spacecrafti for example, suppose the states of other spacecrafts are known, and the measurement model between spacecrafts i and j could be simplified as(20)
z
k
+
1
i
j
=
H
2
,
k
+
1
i
x
k
+
1
i
+
H
2
,
k
+
1
j
x
k
+
1
j
+
v
k
+
1
i
j
,
v
k
+
1
i
j
~
N
0
,
R
k
+
1
i
j
.Notice thatx
k
+
1
j is known. Let z
-
k
+
1
i
j
=
z
k
+
1
i
j
-
H
2
,
k
+
1
j
x
k
+
1
j; the ROEKF measurement model could be rewrite as(21)
Z
-
k
+
1
i
=
H
-
k
+
1
i
x
k
+
1
i
+
V
k
+
1
i
,
V
k
+
1
i
~
N
0
,
R
k
+
1
i
,where Z
-
k
+
1
i is the set of all the measurements related with spacecraft i and H
-
k
+
1
i is the measurement matrix related with Z
-
k
+
1
i.Thus the local estimator at spacecrafti only needs to estimate its own state, x
^
k
i and P
k
i. The measurement update of FOEKF at spacecraft i is as follows:(22)
K
k
+
1
=
P
k
+
1
i
-
H
-
k
+
1
i
T
H
-
k
+
1
i
P
k
+
1
i
-
H
-
k
+
1
i
T
+
R
k
+
1
i
-
1
,
x
^
k
+
1
i
+
=
x
^
k
+
1
i
-
+
K
k
+
1
Z
-
k
+
1
i
-
H
-
k
+
1
i
x
^
k
+
1
i
-
,
P
k
+
1
i
+
=
I
-
K
k
+
1
H
-
k
+
1
i
P
k
+
1
i
-
.Compared with FOEKF, the computational load is reduced remarkably in ROEKF, because ROEKF only updates the state of local spacecraft. However, as only estimates instead of real states of other spacecrafts could be obtained, the errors of state estimation of other spacecrafts are not considered in ROEKF, rendering a further precision loss than FOEKF. We could use a proper increase to the measurement noise, theR matrix, to absorb such errors. This leads to the increased R extended Kalman filter (IREKF) [5].


### 3.3.2. IREKF

Considering the measurementz
k
+
1
i
j, the measurement noise caused by using estimates of the state of spacecraft j could be estimated by the following equation:(23)
E
z
k
+
1
i
j
-
H
2
,
k
+
1
i
x
k
+
1
i
z
k
+
1
i
j
-
H
2
,
k
+
1
i
x
k
+
1
i
T
=
E
H
2
,
k
+
1
j
x
k
+
1
j
+
v
k
+
1
i
j
H
2
,
k
+
1
j
x
k
+
1
j
+
v
k
+
1
i
j
T
=
E
H
2
,
k
+
1
j
x
k
+
1
j
H
2
,
k
+
1
j
x
k
+
1
j
T
+
H
2
,
k
+
1
j
x
k
+
1
j
v
k
+
1
i
j
T
+
v
k
+
1
i
j
H
2
,
k
+
1
j
x
k
+
1
j
T
+
v
k
+
1
i
j
v
k
+
1
i
j
T
=
E
H
2
,
k
+
1
j
x
k
+
1
j
H
2
,
k
+
1
j
x
k
+
1
j
T
+
E
v
k
+
1
i
j
v
k
+
1
i
j
T
=
H
2
,
k
+
1
j
P
k
+
1
j
+
H
2
,
k
+
1
j
T
+
R
k
+
1
i
j
.Notice that only estimates of the states of other spacecrafts could be obtained, which means that we could only useP
k
+
1
j
- to approximate P
k
+
1
j
+. Thus the increase to the measurement noise of R
i
j is given by H
2
,
k
+
1
j
P
k
+
1
j
+
H
2
,
k
+
1
j
T.


## 3.3.1. Original ROEKF

ROEKF is a further simplification of FOEKF. Taking spacecrafti for example, suppose the states of other spacecrafts are known, and the measurement model between spacecrafts i and j could be simplified as(20)
z
k
+
1
i
j
=
H
2
,
k
+
1
i
x
k
+
1
i
+
H
2
,
k
+
1
j
x
k
+
1
j
+
v
k
+
1
i
j
,
v
k
+
1
i
j
~
N
0
,
R
k
+
1
i
j
.Notice thatx
k
+
1
j is known. Let z
-
k
+
1
i
j
=
z
k
+
1
i
j
-
H
2
,
k
+
1
j
x
k
+
1
j; the ROEKF measurement model could be rewrite as(21)
Z
-
k
+
1
i
=
H
-
k
+
1
i
x
k
+
1
i
+
V
k
+
1
i
,
V
k
+
1
i
~
N
0
,
R
k
+
1
i
,where Z
-
k
+
1
i is the set of all the measurements related with spacecraft i and H
-
k
+
1
i is the measurement matrix related with Z
-
k
+
1
i.Thus the local estimator at spacecrafti only needs to estimate its own state, x
^
k
i and P
k
i. The measurement update of FOEKF at spacecraft i is as follows:(22)
K
k
+
1
=
P
k
+
1
i
-
H
-
k
+
1
i
T
H
-
k
+
1
i
P
k
+
1
i
-
H
-
k
+
1
i
T
+
R
k
+
1
i
-
1
,
x
^
k
+
1
i
+
=
x
^
k
+
1
i
-
+
K
k
+
1
Z
-
k
+
1
i
-
H
-
k
+
1
i
x
^
k
+
1
i
-
,
P
k
+
1
i
+
=
I
-
K
k
+
1
H
-
k
+
1
i
P
k
+
1
i
-
.Compared with FOEKF, the computational load is reduced remarkably in ROEKF, because ROEKF only updates the state of local spacecraft. However, as only estimates instead of real states of other spacecrafts could be obtained, the errors of state estimation of other spacecrafts are not considered in ROEKF, rendering a further precision loss than FOEKF. We could use a proper increase to the measurement noise, theR matrix, to absorb such errors. This leads to the increased R extended Kalman filter (IREKF) [5].


## 3.3.2. IREKF

Considering the measurementz
k
+
1
i
j, the measurement noise caused by using estimates of the state of spacecraft j could be estimated by the following equation:(23)
E
z
k
+
1
i
j
-
H
2
,
k
+
1
i
x
k
+
1
i
z
k
+
1
i
j
-
H
2
,
k
+
1
i
x
k
+
1
i
T
=
E
H
2
,
k
+
1
j
x
k
+
1
j
+
v
k
+
1
i
j
H
2
,
k
+
1
j
x
k
+
1
j
+
v
k
+
1
i
j
T
=
E
H
2
,
k
+
1
j
x
k
+
1
j
H
2
,
k
+
1
j
x
k
+
1
j
T
+
H
2
,
k
+
1
j
x
k
+
1
j
v
k
+
1
i
j
T
+
v
k
+
1
i
j
H
2
,
k
+
1
j
x
k
+
1
j
T
+
v
k
+
1
i
j
v
k
+
1
i
j
T
=
E
H
2
,
k
+
1
j
x
k
+
1
j
H
2
,
k
+
1
j
x
k
+
1
j
T
+
E
v
k
+
1
i
j
v
k
+
1
i
j
T
=
H
2
,
k
+
1
j
P
k
+
1
j
+
H
2
,
k
+
1
j
T
+
R
k
+
1
i
j
.Notice that only estimates of the states of other spacecrafts could be obtained, which means that we could only useP
k
+
1
j
- to approximate P
k
+
1
j
+. Thus the increase to the measurement noise of R
i
j is given by H
2
,
k
+
1
j
P
k
+
1
j
+
H
2
,
k
+
1
j
T.


## 3.4. ICEKF

Iteration is a straight forward method to compensate the precision loss caused by approximations. ICEKF could serve as a general scheme to improve FOEKF and ROEKF. In ICEKF, the local estimation at each spacecraft is passed to other spacecrafts for the linearization of measurement functions through a cascade communication chain. Figure1 gives an example of ICEKF for a three-spacecraft formation compared with EKF. Although ICEKF reduces the computational load while obtaining a better precision than FOEKF or ROEKF, the introduction of iteration will result in an increase in the communicational load.Figure 1
ICEKF compared with centralized EKF [15].


## 4. JT Based Decentralized Solution

The decentralized solutions reviewed in the previous section are all based on distributing centralized EKF through simplifications or decoupling, rendering suboptimal estimations with inferior precisions to EKF. This section proposes a general scheme for the design of decentralized algorithms from another point of view: analyzing the hidden decentralized structure of EKF with the aid of probabilistic graphical models.From the perspective of BE, the dynamic model and the measurement model used in EKF are all conditional Gaussian (CG), which could be interpreted by the probabilistic graphical model of DBN. Previous researches have demonstrated that inference on DBN could be distributed by message passing on a decentralized data structure named JT [16]. By properly splitting and implementing the time update and measurement update of BE on different cliques of JT, it is possible to design a decentralized solution which is centralized-equivalent.The algorithm of JT contains two parts: symbolic calculation and mathematical calculation. Figure2 gives the flowchart of JT. Symbolic calculation returns the decentralized structure of JT based on CG distributions defined by the model of cooperative navigation, while mathematical calculation returns the result of inference through message passing on JT. The key elements of the JT algorithm are introduced in the following part of this section.Figure 2
Flowchart of JT.


### 4.1. Inference Engine

BE is the inference engine on JT. LetZ
1
:
k be the history measurements from initial time step to current time step k, and let p
X
k
∣
Z
1
:
k be the best estimation of current state under history measurements. If the conditional probability between new state and current state, p
X
k
+
1
∣
X
k, and the conditional probability between new measurement and new state, p
Z
k
+
1
∣
X
k
+
1, are given, we could use BE to estimate the most likely state of the next time step. The procedure of BE is as follows.Time Update.
(24)
p
X
k
+
1
∣
Z
1
:
k
=
∫
p
X
k
+
1
∣
X
k
·
p
X
k
∣
Z
1
:
k
d
X
k
.Measurement Update.
(25)
p
X
k
+
1
∣
Z
1
:
k
+
1
=
p
X
k
+
1
∣
Z
1
:
k
·
p
Z
k
+
1
∣
X
k
+
1
∫
p
X
k
+
1
∣
Z
1
:
k
·
p
Z
k
+
1
∣
X
k
+
1
d
X
k
+
1
.There are three basic mathematic operations in BE: multiplication, division, and marginalization. Interested readers may refer to AppendixB for details. Notice that p
X
k
∣
Z
1
:
k
=
N
X
^
k
+
,
P
k
+; from (12) we have p
X
k
+
1
∣
X
k
=
N
Φ
k
X
k
,
Q
k, and from (13) we have p
Z
k
+
1
∣
X
k
+
1
=
N
H
k
+
1
X
k
+
1
,
R
k
+
1. It can be verified that, by applying the basic operations on CG potentials to BE, BE returns the same result as EKF. The proof is omitted here for simplicity.


### 4.2. DBN

In the probabilistic graph theory, DBN is used to visualize the relationship of conditional distribution between random variables. Since the dynamic model and the measurement model of EKF each define a CG distribution, the time update and the measurement update of EKF (or BE) may be represented by DBN.Figure3 gives the DBN model of the example in Section 3.4, where the random variable x
k
i
i
=
1,2
,
3 represents the state of each spacecraft. Notice that measurements are hidden behind the solid lines in DBN, while dynamic models are represented by dashed lines. Let U be the full set of all the random variables in DBN. For the example in Figure 3, U
=
x
k
1
,
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3.Figure 3
Example of DBN.DBN holds the joint distribution of all the variables inU. Inference on DBN means to find the joint distribution of interested variables, which is p
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1 for this example. To perform the inference, we need to remove the barren variables in DBN, which are x
k
1, x
k
2, and x
k
3. This procedure is often referred to as barren node removal in the probabilistic graphical model.


### 4.3. Decentralize DBN Based on JT

JT is first introduced by Lauritzen [17] while studying the inference algorithm of the Bayesian network. JT is a tree structure of cliques. A clique C is a subset of the full set U. Each clique is assigned with a potential function, which is a probabilistic distribution of the variables in C. The set of variables shared by two adjacent cliques C
i and C
j is called a separator, denoted by S
i
j
=
C
i
∩
C
j. Measurements are treated as evidences deployed to different cliques. To generate JT from DBN, the sequence of barren node removal must be defined. Moreover, an important structural constraint named the running intersection property [16] must be satisfied.If a variable is contained in cliques
C
i
and
C
j
, then it must also be contained in the (unique) path between
C
i
and
C
j
.The Kruskal algorithm and the Prim algorithm are the most common algorithms to generate JT. Interested readers may refer to [18] for details. When a JT is constructed, it must be initialized with potentials and evidences. Here we use two criteria for the initialization of JT.Criterion 1.
The probabilistic distribution defined by a set of random variables is assigned as a potential to the first clique which contains the set. If a clique is assigned with more than one potential, the multiplication of these potentials is treated as the potential for the clique; if there is no distribution assigned to a clique, the potential of this clique is initialized as standard normal.Criterion 2.
A measurement is deployed as an evidence to the first clique which contains all the variables related to the very measurement.Figure4 gives a JT based on the DBN model in the previous subsection, where the shaded ellipsoids represent the cliques, the hollow ellipsoids represent the evidences, and the dashed lines indicate the allocation of cliques to spacecrafts.Figure 4
Example of JT.


### 4.4. Inference on JT

Inference on JT is driven by a series of local computations on cliques and messages passing between adjacent cliques [19]. As the message passing is directional, an arbitrary clique is picked up as the root of JT. The message passing starts from leaves to root, and then it rolls back from root to leaves. If all cliques have received messages from both its parents and children, the inference at this time step is finished. A remarkable advantage of JT is its good modularity, rendering flexibility for reconfigurations and robustness after spacecraft failures. The operations at each clique are in a uniformed manor of three steps: receiving message, collecting evidence, and generating message.


#### 4.4.1. Receiving Message

A message to a clique is defined as a distribution (or potential) passed from an adjacent clique. If cliqueC
i have received a message from clique C
j, the potential of C
i is updated by the following equation through multiplication:(26)
p
∗
C
i
=
p
C
i
·
p
M
s
g
,where p
C
j is the original potential of clique C
i, p
M
s
g is the message, and p
∗
C
i is the new potential assigned to C
i.In fact, the operation of receiving message is a distributed implementation of the time update in centralized BE (or EKF).


#### 4.4.2. Collecting Evidence

Take cliqueC
i for example. All evidences deployed to C
i are in the form of p
Z
∣
C
i. To collect the evidence of Z means to use the measurement of Z to update the estimate of local state, which is the distribution of state conditioned on Z. The potential of C
i after collecting evidence is updated by(27)
p
∗
C
i
=
p
C
i
·
p
Z
∣
C
i
∫
C
i
·
p
Z
∣
C
i
d
C
i
.The operation of collecting evidence may be seen as a distributed implementation of the measurement update in centralized BE (or EKF).


#### 4.4.3. Generating Message

After updating the potential of a clique based on messages and evidences, the next operation is to generate the message passed to the next clique. The message generated is based on the separator between adjacent cliques. Considering cliqueC
i generating a message to clique C
j, let S
-
i
j
i be the complementary set of the separator S
i
j in set C
i. The message from C
i to C
j is defined as the distribution of S
i
j, p
S
i
j. To get p
S
i
j, we need to remove the variables in S
-
i
j
i from C
i through marginalization. The message is given by(28)
p
M
s
g
=
∫
p
C
i
d
S
-
i
j
i
.Example 1.
Consider the example of JT in Figure4. Select clique C
3 as the root of JT. The operations at each spacecraft return the following results.(i)
Spacecraft 1:Receiving Message.
(29)
p
x
k
1
,
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
=
p
x
k
+
1
1
∣
x
k
1
·
p
x
k
1
,
x
k
3
,
x
k
3
∣
Z
1
:
k
.Generating Message.
(30)
p
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
=
∫
p
x
k
1
,
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
d
x
k
1
.
(ii)
Spacecraft 2:Receiving Message.
(31)
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
=
p
x
k
+
1
2
∣
x
k
2
·
p
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
.Collecting Evidence.
(32)
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
=
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
·
p
z
k
+
1
12
∣
x
k
+
1
1
,
x
k
+
1
2
p
z
k
+
1
12
.Generating Message.
(33)
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
=
∫
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
d
x
k
2
.
(iii)
Spacecraft 3:Receiving Message.
(34)
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
,
z
k
+
1
12
=
p
x
k
+
1
3
∣
x
k
3
·
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
.Collecting Evidence.
(35)
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1
=
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
,
z
k
+
1
12
p
z
k
+
1
13
·
p
z
k
+
1
23
·
p
z
k
+
1
13
∣
x
k
+
1
1
,
x
k
+
1
3
·
p
z
k
+
1
23
∣
x
k
+
1
2
,
x
k
+
1
3
.Generating Message.
(36)
p
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1
=
∫
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1
d
x
k
3
.The message generated by spacecraft 3,p
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1, is the final result at time step k
+
1, which is the same as centralized BE.


## 4.1. Inference Engine

BE is the inference engine on JT. LetZ
1
:
k be the history measurements from initial time step to current time step k, and let p
X
k
∣
Z
1
:
k be the best estimation of current state under history measurements. If the conditional probability between new state and current state, p
X
k
+
1
∣
X
k, and the conditional probability between new measurement and new state, p
Z
k
+
1
∣
X
k
+
1, are given, we could use BE to estimate the most likely state of the next time step. The procedure of BE is as follows.Time Update.
(24)
p
X
k
+
1
∣
Z
1
:
k
=
∫
p
X
k
+
1
∣
X
k
·
p
X
k
∣
Z
1
:
k
d
X
k
.Measurement Update.
(25)
p
X
k
+
1
∣
Z
1
:
k
+
1
=
p
X
k
+
1
∣
Z
1
:
k
·
p
Z
k
+
1
∣
X
k
+
1
∫
p
X
k
+
1
∣
Z
1
:
k
·
p
Z
k
+
1
∣
X
k
+
1
d
X
k
+
1
.There are three basic mathematic operations in BE: multiplication, division, and marginalization. Interested readers may refer to AppendixB for details. Notice that p
X
k
∣
Z
1
:
k
=
N
X
^
k
+
,
P
k
+; from (12) we have p
X
k
+
1
∣
X
k
=
N
Φ
k
X
k
,
Q
k, and from (13) we have p
Z
k
+
1
∣
X
k
+
1
=
N
H
k
+
1
X
k
+
1
,
R
k
+
1. It can be verified that, by applying the basic operations on CG potentials to BE, BE returns the same result as EKF. The proof is omitted here for simplicity.


## 4.2. DBN

In the probabilistic graph theory, DBN is used to visualize the relationship of conditional distribution between random variables. Since the dynamic model and the measurement model of EKF each define a CG distribution, the time update and the measurement update of EKF (or BE) may be represented by DBN.Figure3 gives the DBN model of the example in Section 3.4, where the random variable x
k
i
i
=
1,2
,
3 represents the state of each spacecraft. Notice that measurements are hidden behind the solid lines in DBN, while dynamic models are represented by dashed lines. Let U be the full set of all the random variables in DBN. For the example in Figure 3, U
=
x
k
1
,
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3.Figure 3
Example of DBN.DBN holds the joint distribution of all the variables inU. Inference on DBN means to find the joint distribution of interested variables, which is p
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1 for this example. To perform the inference, we need to remove the barren variables in DBN, which are x
k
1, x
k
2, and x
k
3. This procedure is often referred to as barren node removal in the probabilistic graphical model.


## 4.3. Decentralize DBN Based on JT

JT is first introduced by Lauritzen [17] while studying the inference algorithm of the Bayesian network. JT is a tree structure of cliques. A clique C is a subset of the full set U. Each clique is assigned with a potential function, which is a probabilistic distribution of the variables in C. The set of variables shared by two adjacent cliques C
i and C
j is called a separator, denoted by S
i
j
=
C
i
∩
C
j. Measurements are treated as evidences deployed to different cliques. To generate JT from DBN, the sequence of barren node removal must be defined. Moreover, an important structural constraint named the running intersection property [16] must be satisfied.If a variable is contained in cliques
C
i
and
C
j
, then it must also be contained in the (unique) path between
C
i
and
C
j
.The Kruskal algorithm and the Prim algorithm are the most common algorithms to generate JT. Interested readers may refer to [18] for details. When a JT is constructed, it must be initialized with potentials and evidences. Here we use two criteria for the initialization of JT.Criterion 1.
The probabilistic distribution defined by a set of random variables is assigned as a potential to the first clique which contains the set. If a clique is assigned with more than one potential, the multiplication of these potentials is treated as the potential for the clique; if there is no distribution assigned to a clique, the potential of this clique is initialized as standard normal.Criterion 2.
A measurement is deployed as an evidence to the first clique which contains all the variables related to the very measurement.Figure4 gives a JT based on the DBN model in the previous subsection, where the shaded ellipsoids represent the cliques, the hollow ellipsoids represent the evidences, and the dashed lines indicate the allocation of cliques to spacecrafts.Figure 4
Example of JT.


## 4.4. Inference on JT

Inference on JT is driven by a series of local computations on cliques and messages passing between adjacent cliques [19]. As the message passing is directional, an arbitrary clique is picked up as the root of JT. The message passing starts from leaves to root, and then it rolls back from root to leaves. If all cliques have received messages from both its parents and children, the inference at this time step is finished. A remarkable advantage of JT is its good modularity, rendering flexibility for reconfigurations and robustness after spacecraft failures. The operations at each clique are in a uniformed manor of three steps: receiving message, collecting evidence, and generating message.


### 4.4.1. Receiving Message

A message to a clique is defined as a distribution (or potential) passed from an adjacent clique. If cliqueC
i have received a message from clique C
j, the potential of C
i is updated by the following equation through multiplication:(26)
p
∗
C
i
=
p
C
i
·
p
M
s
g
,where p
C
j is the original potential of clique C
i, p
M
s
g is the message, and p
∗
C
i is the new potential assigned to C
i.In fact, the operation of receiving message is a distributed implementation of the time update in centralized BE (or EKF).


### 4.4.2. Collecting Evidence

Take cliqueC
i for example. All evidences deployed to C
i are in the form of p
Z
∣
C
i. To collect the evidence of Z means to use the measurement of Z to update the estimate of local state, which is the distribution of state conditioned on Z. The potential of C
i after collecting evidence is updated by(27)
p
∗
C
i
=
p
C
i
·
p
Z
∣
C
i
∫
C
i
·
p
Z
∣
C
i
d
C
i
.The operation of collecting evidence may be seen as a distributed implementation of the measurement update in centralized BE (or EKF).


### 4.4.3. Generating Message

After updating the potential of a clique based on messages and evidences, the next operation is to generate the message passed to the next clique. The message generated is based on the separator between adjacent cliques. Considering cliqueC
i generating a message to clique C
j, let S
-
i
j
i be the complementary set of the separator S
i
j in set C
i. The message from C
i to C
j is defined as the distribution of S
i
j, p
S
i
j. To get p
S
i
j, we need to remove the variables in S
-
i
j
i from C
i through marginalization. The message is given by(28)
p
M
s
g
=
∫
p
C
i
d
S
-
i
j
i
.Example 1.
Consider the example of JT in Figure4. Select clique C
3 as the root of JT. The operations at each spacecraft return the following results.(i)
Spacecraft 1:Receiving Message.
(29)
p
x
k
1
,
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
=
p
x
k
+
1
1
∣
x
k
1
·
p
x
k
1
,
x
k
3
,
x
k
3
∣
Z
1
:
k
.Generating Message.
(30)
p
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
=
∫
p
x
k
1
,
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
d
x
k
1
.
(ii)
Spacecraft 2:Receiving Message.
(31)
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
=
p
x
k
+
1
2
∣
x
k
2
·
p
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
.Collecting Evidence.
(32)
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
=
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
·
p
z
k
+
1
12
∣
x
k
+
1
1
,
x
k
+
1
2
p
z
k
+
1
12
.Generating Message.
(33)
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
=
∫
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
d
x
k
2
.
(iii)
Spacecraft 3:Receiving Message.
(34)
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
,
z
k
+
1
12
=
p
x
k
+
1
3
∣
x
k
3
·
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
.Collecting Evidence.
(35)
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1
=
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
,
z
k
+
1
12
p
z
k
+
1
13
·
p
z
k
+
1
23
·
p
z
k
+
1
13
∣
x
k
+
1
1
,
x
k
+
1
3
·
p
z
k
+
1
23
∣
x
k
+
1
2
,
x
k
+
1
3
.Generating Message.
(36)
p
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1
=
∫
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1
d
x
k
3
.The message generated by spacecraft 3,p
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1, is the final result at time step k
+
1, which is the same as centralized BE.


## 4.4.1. Receiving Message

A message to a clique is defined as a distribution (or potential) passed from an adjacent clique. If cliqueC
i have received a message from clique C
j, the potential of C
i is updated by the following equation through multiplication:(26)
p
∗
C
i
=
p
C
i
·
p
M
s
g
,where p
C
j is the original potential of clique C
i, p
M
s
g is the message, and p
∗
C
i is the new potential assigned to C
i.In fact, the operation of receiving message is a distributed implementation of the time update in centralized BE (or EKF).


## 4.4.2. Collecting Evidence

Take cliqueC
i for example. All evidences deployed to C
i are in the form of p
Z
∣
C
i. To collect the evidence of Z means to use the measurement of Z to update the estimate of local state, which is the distribution of state conditioned on Z. The potential of C
i after collecting evidence is updated by(27)
p
∗
C
i
=
p
C
i
·
p
Z
∣
C
i
∫
C
i
·
p
Z
∣
C
i
d
C
i
.The operation of collecting evidence may be seen as a distributed implementation of the measurement update in centralized BE (or EKF).


## 4.4.3. Generating Message

After updating the potential of a clique based on messages and evidences, the next operation is to generate the message passed to the next clique. The message generated is based on the separator between adjacent cliques. Considering cliqueC
i generating a message to clique C
j, let S
-
i
j
i be the complementary set of the separator S
i
j in set C
i. The message from C
i to C
j is defined as the distribution of S
i
j, p
S
i
j. To get p
S
i
j, we need to remove the variables in S
-
i
j
i from C
i through marginalization. The message is given by(28)
p
M
s
g
=
∫
p
C
i
d
S
-
i
j
i
.Example 1.
Consider the example of JT in Figure4. Select clique C
3 as the root of JT. The operations at each spacecraft return the following results.(i)
Spacecraft 1:Receiving Message.
(29)
p
x
k
1
,
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
=
p
x
k
+
1
1
∣
x
k
1
·
p
x
k
1
,
x
k
3
,
x
k
3
∣
Z
1
:
k
.Generating Message.
(30)
p
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
=
∫
p
x
k
1
,
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
d
x
k
1
.
(ii)
Spacecraft 2:Receiving Message.
(31)
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
=
p
x
k
+
1
2
∣
x
k
2
·
p
x
k
2
,
x
k
3
,
x
k
+
1
1
∣
Z
1
:
k
.Collecting Evidence.
(32)
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
=
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
·
p
z
k
+
1
12
∣
x
k
+
1
1
,
x
k
+
1
2
p
z
k
+
1
12
.Generating Message.
(33)
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
=
∫
p
x
k
2
,
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
d
x
k
2
.
(iii)
Spacecraft 3:Receiving Message.
(34)
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
,
z
k
+
1
12
=
p
x
k
+
1
3
∣
x
k
3
·
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
∣
Z
1
:
k
,
z
k
+
1
12
.Collecting Evidence.
(35)
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1
=
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
,
z
k
+
1
12
p
z
k
+
1
13
·
p
z
k
+
1
23
·
p
z
k
+
1
13
∣
x
k
+
1
1
,
x
k
+
1
3
·
p
z
k
+
1
23
∣
x
k
+
1
2
,
x
k
+
1
3
.Generating Message.
(36)
p
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1
=
∫
p
x
k
3
,
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1
d
x
k
3
.The message generated by spacecraft 3,p
x
k
+
1
1
,
x
k
+
1
2
,
x
k
+
1
3
∣
Z
1
:
k
+
1, is the final result at time step k
+
1, which is the same as centralized BE.


## 5. Simulation Results and Discussion

The effectiveness of the decentralized solution proposed is validated by the simulation of a three-satellite formation. The Keplerian elements of the formation satellites are listed in Table1. For simplicity, only central body gravitation is considered (also known as the two-body problem). In the formation, satellite 1 is selected as the leading spacecraft, with the other two accompanying it. The relative coordinate system (or the Hill frame), O
s
-
X
r
Y
r
Z
r, is based on satellite 1, where the X axis is from the center of Earth to the center of satellite 1, the Y axis is along the direction of velocity in the orbital plane, and the Z axis is established using the right-hand rule. The relationship between O
s
-
X
r
Y
r
Z
r and the Earth-Centered Inertial (ECI) coordinate system, O
e
-
X
i
Y
i
Z
i, is given in Figure 5. Under the two-body assumption, the relative motion of the formation is given in Figure 6.Table 1
Keplerian orbital elements of a satellite formation.
Sat number
a (km)
e
i (°)
Ω (°)
ω (°)
M
0 (°)
1
8400.000000
0.005000
30.000000
120.000000
0.000000
90.000000
2
8400.000000
0.004743
29.996649
119.988151
1.783138
88.250782
3
8400.000000
0.004743
30.003353
120.011847
1.762616
88.250782Figure 5
Relative coordinate system and ECI coordinate system. The relative coordinate system (or the Hill frame) is represented byO
s
-
X
r
Y
r
Z
r, and the ECI coordinate system is represented by O
e
-
X
i
Y
i
Z
i.Figure 6
Relative motion between satellites.


### 5.1. Precision of Decentralized Algorithms

Every spacecraft in the formation is equipped with a GPS receiver as well as RF range and communication devices. All noised are treated as Gaussian. The standard deviation for GPS measurement and RF range are 0.3 m and 0.1 m, respectively. The estimations of the absolute state of satellite 3 using different algorithms are listed as follows.Centralized EKF is setup as a standard to evaluate the precision of other algorithms. As there are no further simplifications or decoupling, the precision of centralize EKF is better than all the decentralized solutions reviewed in Section2. Figure 7 illustrates the error of estimation of the absolute state of satellite 3 in the ECI coordinate system using centralized EKF.Figure 7
Error of absolute position estimation for satellite 3 using EKF.All decentralized algorithms are compared with centralized EKF. For FOEKF, because only local measurements related with each satellite are considered, the simplification of measurement model will result in suboptimal estimations. Figure8 gives result of absolute state estimation of satellite 3 in the ECI coordinate system using FOEKF. The root-mean-square error (RMSE) tends to fluctuate around 0.05 m, which is inferior to centralized EKF.Figure 8
Error of absolute position estimation of satellite 3 using FOEKF.The core idea that differs ROEKF from FOEKF is that the states transferred from other spacecrafts are assumed to be precise, and therefore the interformation measurement functions could be decoupled. After decoupling, the measurement function is only related to the state of the local spacecraft. This assumption will result in a further precision loss of ROEKF compared with FOEKF. Figure9 shows that a notable precision loss is caused by this simplification compared with Figure 8.Figure 9
Error of absolute position estimation for satellite 3 using ROEKF.Notice that the states transmitted from other spacecrafts are not precise, in other words, with an uncertainty represented by the covariance. In IREKF, such uncertainty is considered as a proper bump-up to the noise of measurements. Figure10 gives an example of satellite 3 using IREKF. The results are greatly improved in IREKF compared with ROEKF. However, because the covariance matrixes transferred from other spacecrafts are generated on local measurements, IREKF is suboptimal. The precision of IREKF is still inferior to EKF. Comparing Figure 10 with Figure 8, it can be seen that IREKF tends to approximate FOEKF.Figure 10
Error of absolute position estimation for satellite 3 using IREKF.A common method to compensate the precision loss caused by simplification or linearization is to use iterations. ICEKF is based on the same intuition. Figure11 gives an example of ICEKF based on FOEKF with 5 iterations per epoch. Without iterations, the algorithm is the same as FOEKF. It can be seen that the result has been improved remarkably compared with FOEKF. However, the result is still slightly inferior to centralized EKF.Figure 11
Error of absolute position estimation for satellite 3 using FOEKF-based ICEKF.Figure12 gives an example of the decentralized solution based on JT. As there are no decoupling or further assumptions, the result of JT is equivalent to centralized EKF.Figure 12
Error of absolute position estimation for satellite 3 using JT.In order to provide an intuitive comparison of different algorithms, Table2 summarizes the RMSE of all algorithms. Because some algorithms fail to reach convergence, the end of simulation is picked up as the epoch when the RMSE of different algorithms is evaluated. It can be seen that JT has the best precision among all the decentralized algorithms discussed in this paper.Table 2
RMSE of all algorithms at the end of simulation.
Algorithm
EKF
FOEKF
ROEKF
IREKF
ICEKF-FOEKF
JT
RMSE (m)
0.0248
0.0739
28.1360
0.0690
0.0310
0.0248


### 5.2. Analysis of Computational Complexity

While analyzing the computational complexity of algorithms, the operation of multiplication of two real numbers is defined as a unit element. Recall the procedure of centralized EKF in Section3.1. Suppose x
∈
R
m
×
1 and Z
∈
R
d
×
1. Notice that (16) requires matrix inversio, and take elementary row transformation as the method for matrix inversion. The computational complexity of (14), (15), (16), (17), and (18) is O
m
2
n
2, O
m
3
n
3, O
d
3
+
O
d
2
m
n
+
O
d
m
2
n
2, O
d
m
n, and O
d
m
2
n
2
+
O
m
3
n
3, respectively. Thus the total computational complexity of EKF time update and measurement update is O
m
3
n
3 and O
m
3
n
3
+
O
d
m
2
n
2
+
O
d
2
m
n
+
O
d
3, respectively.The computational complexity of decentralized algorithms could be estimated using similar methods. For FOEKF, supposeZ
k
+
1
i
∈
R
d
i
×
1. The computational complexity of FOEKF time update and measurement update at spacecraft i is O
m
3
n
3 and O
m
3
n
3
+
O
d
i
m
2
n
2
+
O
d
i
2
m
n
+
O
d
i
3, respectively. Notice that, in ROEKF, time update for every spacecraft in the formation is carried out at each local filter, and the computational complexity of decoupling at spacecraft i is O
d
i
m spacecraft i. Thus the computational complexity of ROEKF time update and measurement update at spacecraft i is O
m
3
n and O
m
3
+
O
d
i
m
2
+
O
d
i
2
m
+
O
d
i
3, respectively. Because the calculation of the increase to the measurement noise of R is O
d
i
m
2, the computational complexity of IREKF is at the same level of ROEKF.Because centralized BE is equivalent to centralized EKF, the total computational complexity of JT is the same as centralized EKF. For spacecrafti, suppose that the number of evidences deployed to clique i is d
i
′. Thus the computational complexity of JT at spacecraft i is O
m
3
n
3
+
O
d
i
′
m
2
n
2
+
O
d
i
′
2
m
n
+
O
d
i
′
3. Notice that, in JT, the measurement update of centralized BE is split and implemented in different cliques of JT. Thus d
i
′
≤
d
i. Hence we could sort the local computational complexity of decentralized algorithms as ROEKF < IREKF < JT ⩽ FOEKF.


## 5.1. Precision of Decentralized Algorithms

Every spacecraft in the formation is equipped with a GPS receiver as well as RF range and communication devices. All noised are treated as Gaussian. The standard deviation for GPS measurement and RF range are 0.3 m and 0.1 m, respectively. The estimations of the absolute state of satellite 3 using different algorithms are listed as follows.Centralized EKF is setup as a standard to evaluate the precision of other algorithms. As there are no further simplifications or decoupling, the precision of centralize EKF is better than all the decentralized solutions reviewed in Section2. Figure 7 illustrates the error of estimation of the absolute state of satellite 3 in the ECI coordinate system using centralized EKF.Figure 7
Error of absolute position estimation for satellite 3 using EKF.All decentralized algorithms are compared with centralized EKF. For FOEKF, because only local measurements related with each satellite are considered, the simplification of measurement model will result in suboptimal estimations. Figure8 gives result of absolute state estimation of satellite 3 in the ECI coordinate system using FOEKF. The root-mean-square error (RMSE) tends to fluctuate around 0.05 m, which is inferior to centralized EKF.Figure 8
Error of absolute position estimation of satellite 3 using FOEKF.The core idea that differs ROEKF from FOEKF is that the states transferred from other spacecrafts are assumed to be precise, and therefore the interformation measurement functions could be decoupled. After decoupling, the measurement function is only related to the state of the local spacecraft. This assumption will result in a further precision loss of ROEKF compared with FOEKF. Figure9 shows that a notable precision loss is caused by this simplification compared with Figure 8.Figure 9
Error of absolute position estimation for satellite 3 using ROEKF.Notice that the states transmitted from other spacecrafts are not precise, in other words, with an uncertainty represented by the covariance. In IREKF, such uncertainty is considered as a proper bump-up to the noise of measurements. Figure10 gives an example of satellite 3 using IREKF. The results are greatly improved in IREKF compared with ROEKF. However, because the covariance matrixes transferred from other spacecrafts are generated on local measurements, IREKF is suboptimal. The precision of IREKF is still inferior to EKF. Comparing Figure 10 with Figure 8, it can be seen that IREKF tends to approximate FOEKF.Figure 10
Error of absolute position estimation for satellite 3 using IREKF.A common method to compensate the precision loss caused by simplification or linearization is to use iterations. ICEKF is based on the same intuition. Figure11 gives an example of ICEKF based on FOEKF with 5 iterations per epoch. Without iterations, the algorithm is the same as FOEKF. It can be seen that the result has been improved remarkably compared with FOEKF. However, the result is still slightly inferior to centralized EKF.Figure 11
Error of absolute position estimation for satellite 3 using FOEKF-based ICEKF.Figure12 gives an example of the decentralized solution based on JT. As there are no decoupling or further assumptions, the result of JT is equivalent to centralized EKF.Figure 12
Error of absolute position estimation for satellite 3 using JT.In order to provide an intuitive comparison of different algorithms, Table2 summarizes the RMSE of all algorithms. Because some algorithms fail to reach convergence, the end of simulation is picked up as the epoch when the RMSE of different algorithms is evaluated. It can be seen that JT has the best precision among all the decentralized algorithms discussed in this paper.Table 2
RMSE of all algorithms at the end of simulation.
Algorithm
EKF
FOEKF
ROEKF
IREKF
ICEKF-FOEKF
JT
RMSE (m)
0.0248
0.0739
28.1360
0.0690
0.0310
0.0248


## 5.2. Analysis of Computational Complexity

While analyzing the computational complexity of algorithms, the operation of multiplication of two real numbers is defined as a unit element. Recall the procedure of centralized EKF in Section3.1. Suppose x
∈
R
m
×
1 and Z
∈
R
d
×
1. Notice that (16) requires matrix inversio, and take elementary row transformation as the method for matrix inversion. The computational complexity of (14), (15), (16), (17), and (18) is O
m
2
n
2, O
m
3
n
3, O
d
3
+
O
d
2
m
n
+
O
d
m
2
n
2, O
d
m
n, and O
d
m
2
n
2
+
O
m
3
n
3, respectively. Thus the total computational complexity of EKF time update and measurement update is O
m
3
n
3 and O
m
3
n
3
+
O
d
m
2
n
2
+
O
d
2
m
n
+
O
d
3, respectively.The computational complexity of decentralized algorithms could be estimated using similar methods. For FOEKF, supposeZ
k
+
1
i
∈
R
d
i
×
1. The computational complexity of FOEKF time update and measurement update at spacecraft i is O
m
3
n
3 and O
m
3
n
3
+
O
d
i
m
2
n
2
+
O
d
i
2
m
n
+
O
d
i
3, respectively. Notice that, in ROEKF, time update for every spacecraft in the formation is carried out at each local filter, and the computational complexity of decoupling at spacecraft i is O
d
i
m spacecraft i. Thus the computational complexity of ROEKF time update and measurement update at spacecraft i is O
m
3
n and O
m
3
+
O
d
i
m
2
+
O
d
i
2
m
+
O
d
i
3, respectively. Because the calculation of the increase to the measurement noise of R is O
d
i
m
2, the computational complexity of IREKF is at the same level of ROEKF.Because centralized BE is equivalent to centralized EKF, the total computational complexity of JT is the same as centralized EKF. For spacecrafti, suppose that the number of evidences deployed to clique i is d
i
′. Thus the computational complexity of JT at spacecraft i is O
m
3
n
3
+
O
d
i
′
m
2
n
2
+
O
d
i
′
2
m
n
+
O
d
i
′
3. Notice that, in JT, the measurement update of centralized BE is split and implemented in different cliques of JT. Thus d
i
′
≤
d
i. Hence we could sort the local computational complexity of decentralized algorithms as ROEKF < IREKF < JT ⩽ FOEKF.


## 6. Conclusions

The decentralized state estimation algorithms for formation flying spacecraft studied in this paper could be categorized as two main groups: decomposition through simplification and decoupling and decomposition through structural analysis of models. In the first category, irrelevant measurements are ignored and interformation measurements are decoupled at each spacecraft, rendering suboptimal estimates with inferior precisions to centralized EKF. All algorithms reviewed in Section3 belong to the first category. In the second category, states of different spacecrafts are categorized as cliques of JT and measurements are deployed on different cliques. Since there is no simplification, the new solution based on JT is equivalent in precision to centralized EKF in theory. Simulation results indicate that the new solution of JT holds the best precision among all decentralized algorithms discussed in this paper. For any state estimation applications whose model is in the same form as cooperative navigation, JT may serve as a general scheme for developing decentralized algorithms.



---

*Source: 714302-2015-10-18.xml*