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

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 ismodeled 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 centralizedKF 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.


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 Section 2, 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.

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 Appendix A.

Dynamic Model.
The state of each spacecraft is defined as the Keplerian orbital elements.Let x   be the Keplerian orbital elements of satellite  at time step , x = (, ,   , Ω, , )  , where  is the semimajor axis,  is the eccentricity,   is the inclination angle, Ω is the right ascension of ascending node,  is the argument of perigee, and  is the mean anomaly at reference time.Without considering system inputs such as orbit control, the dynamic model of satellite  is given by where w   is the noise and w   ∼ N(0, Q   ), N(0, Q   ) is a Gaussian distribution of mean 0 and variance Q   .Linearize (1) at x  using first order Tyler expansion, and we have where F   is the state transition matrix and (3)

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.

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.Let r  +1 be the position of observed GPS satellite given by the broadcast ephemeris, and let r  +1 be the position of spacecraft  at time step  + 1.The single platform measurement based on GPS pseudorange,   +1 , is given by where |r| means the length of vector r and V  +1 is the measurement noise, V  +1 ∼ N(0,   +1 ).In general, the single platform measurement function may be written as Linearize (5) at x +1 using first order Tyler expansion, and we have where H  1,+1 is the measurement matrix and where ).In general, the single platform measurement function may be written as Linearize (9) at x +1 and x +1 using first order Tyler expansion, and we have

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].

Centralized EKF.
Consider a formation of  spacecrafts.Let X  be the state of the whole formation at time step , . Without considering the system input, the linearized model of cooperative navigation could be rewritten as follows.
Dynamic Model.
Measurement Model.
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.
Measurement Update.
where mean X+  and covariance P +  are the best estimates of system state at time step .

FOEKF.
Considering a local FOEKF estimator on spacecraft , it only processes the measurements related to spacecraft , which means that the measurement Z  +1 , measurement matrix H  +1 , and measurement noise V  +1 are subsets of Z +1 , H +1 , and V +1 in EKF, respectively.The measurement update of FOEKF at spacecraft  is as follows: 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.

Original ROEKF.
ROEKF is a further simplification of FOEKF.Taking spacecraft  for example, suppose the states of other spacecrafts are known, and the measurement model between spacecrafts  and  could be simplified as ; the ROEKF measurement model could be rewrite as where 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, the R matrix, to absorb such errors.This leads to the increased R extended Kalman filter (IREKF) [5].

IREKF. Considering the measurement 𝑧
+1 , the measurement noise caused by using estimates of the state of spacecraft  could be estimated by the following equation: Notice that only estimates of the states of other spacecrafts could be obtained, which means that we could only use P − +1 to approximate P + +1 .Thus the increase to the measurement noise of   is given by H

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.Figure 1 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.

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.Figure 2 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.

𝑝 (X
Measurement Update. There are three basic mathematic operations in BE: multiplication, division, and marginalization.Interested readers may refer to Appendix B for details.Notice that (X  | Z 1: ) = N( X+  , P +  ); from (12) we have (X +1 | X  ) = N(Φ  X  , Q  ), and from (13) we have (Z +1 | X +1 ) = N(H +1 X +1 , R +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.

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.
Figure 3 gives the DBN model of the example in Section 3.4, where the random variable x   ( = 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  be the full set of all the random variables in DBN.For the example in Figure 3,  = {x 1  , x 2  , x 3  , x 1 +1 , x 2 +1 , x 3 +1 }.DBN holds the joint distribution of all the variables in .Inference on DBN means to find the joint distribution of interested variables, which is (x 1 +1 , x 2 +1 , x 3 +1 | Z 1:+1 ) for this example.To perform the inference, we need to remove the barren variables in DBN, which are x 1   , x 2  , and x 3  .This procedure is often referred to as barren node removal in the probabilistic graphical model.Moreover, an important structural constraint named the running intersection property [16] must be satisfied.
If a variable is contained in cliques   and   , then it must also be contained in the (unique) path between   and   .
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.
Figure 4 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.[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.

Receiving Message.
A message to a clique is defined as a distribution (or potential) passed from an adjacent clique.
If clique   have received a message from clique   , the potential of   is updated by the following equation through multiplication: where (  ) is the original potential of clique   , (Msg) is the message, and  * (  ) is the new potential assigned to   .In fact, the operation of receiving message is a distributed implementation of the time update in centralized BE (or EKF).

Collecting Evidence.
Take clique   for example.All evidences deployed to   are in the form of ( |   ).To collect the evidence of  means to use the measurement of  to update the estimate of local state, which is the distribution of state conditioned on .The potential of   after collecting evidence is updated by The operation of collecting evidence may be seen as a distributed implementation of the measurement update in centralized BE (or EKF).

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 clique   generating a message to clique   , let    be the complementary set of the separator   in set   .The message from   to   is defined as the distribution of   , (  ).To get (  ), we need to remove the variables in    from   through marginalization.The message is given by Example 1.Consider the example of JT in Figure 4. Select clique  3 as the root of JT.The operations at each spacecraft return the following results.
Collecting Evidence.
Generating Message.
Collecting Evidence.
Generating Message.

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 Table 1.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),   −       , is based on satellite 1, where the  axis is from the center of Earth to the center of satellite 1, the  axis is along the direction of velocity in the orbital plane, and the  axis is established using the righthand rule.The relationship between   −       and the Earth-Centered Inertial (ECI) coordinate system,   −      , is given in Figure 5.Under the two-body assumption, the relative motion of the formation is given in Figure 6.

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 Section 2. Figure 7 illustrates the error of estimation of the absolute state of satellite 3 in the ECI coordinate system using centralized 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.Figure 8 gives result of absolute state estimation of satellite 3 in the ECI coordinate system using FOEKF.The root-meansquare error (RMSE) tends to fluctuate around 0.05 m, which is inferior to centralized EKF.
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. Figure 9 shows that a notable precision loss is caused by this simplification compared with Figure 8.
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.Figure 10 gives an example of satellite 3 using IREKF.The results are greatly improved in IREKF compared with ROEKF.However, because the covariance matrixes  10 with Figure 8, it can be seen that IREKF tends to approximate FOEKF.
A common method to compensate the precision loss caused by simplification or linearization is to use iterations.ICEKF is based on the same intuition.Figure 11 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 12 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.
In order to provide an intuitive comparison of different algorithms, Table 2 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.
The computational complexity of decentralized algorithms could be estimated using similar methods.For FOEKF, suppose Z  +1 ∈ R   ×1 .The computational complexity of FOEKF time update and measurement update at spacecraft  is ( 3  3 ) and ( 3  3 )+(   2  2 )+( 2  )+ ( 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  is (  ) spacecraft .Thus the computational complexity of ROEKF time update and measurement update at spacecraft  is ( 3 ) and ( 3 ) + (   2 ) + ( 2  ) + ( 3  ), respectively.Because the calculation of the increase to the measurement noise of R is (   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 spacecraft , suppose that the number of evidences deployed to clique  is    .Thus the computational complexity of JT at spacecraft  is ( 3  3 ) + (    2  2 ) + ( 2  ) + ( 3  ).Notice that, in JT, the measurement update of centralized BE is split and implemented in different cliques of JT.Thus    ≤   .Hence we could sort the local computational complexity of decentralized algorithms as ROEKF < IREKF < JT ⩽ FOEKF.

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 Section 3 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.

Figure 7 :
Figure 7: Error of absolute position estimation for satellite 3 using EKF.

Figure 8 :
Figure 8: Error of absolute position estimation of satellite 3 using FOEKF.

Table 1 :
Keplerian orbital elements of a satellite formation.Figure 5: Relative coordinate system and ECI coordinate system.The relative coordinate system (or the Hill frame) is represented by   −       , and the ECI coordinate system is represented by   −       .

Table 2 :
RMSE of all algorithms at the end of simulation.