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.
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 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.
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 Appendix A.
2.1. Dynamic Model
The state of each spacecraft is defined as the Keplerian orbital elements. Let be the Keplerian orbital elements of satellite at time step , , 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 bywhere is the noise and , is a Gaussian distribution of mean and variance .
Linearize (1) at using first order Tyler expansion, and we havewhere is the state transition matrix and
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. Let be the position of observed GPS satellite given by the broadcast ephemeris, and let be the position of spacecraft at time step . The single platform measurement based on GPS pseudorange, , is given bywhere means the length of vector and is the measurement noise, .
In general, the single platform measurement function may be written as
Linearize (5) at using first order Tyler expansion, and we havewhere is the measurement matrix and
2.2.2. Interformation Measurement
Similar to single platform measurement, let and be the position of spacecraft and at time step , respectively. The interformation measurement between spacecrafts and , , is given bywhere is the measurement noise, .
In general, the single platform measurement function may be written as
Linearize (9) at and using first order Tyler expansion, and we havewhere and are the measurement matrix and
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 of spacecrafts. Let 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, is the noise of the dynamic model, is the vector of measurements, is the measurement matrix, and 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 and covariance are the best estimates of system state at time step .
3.2. FOEKF
Considering a local FOEKF estimator on spacecraft , it only processes the measurements related to spacecraft , which means that the measurement , measurement matrix , and measurement noise are subsets of , , and 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.
3.3. ROEKF
3.3.1. 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
Notice that is known. Let ; the ROEKF measurement model could be rewrite aswhere is the set of all the measurements related with spacecraft and is the measurement matrix related with .
Thus the local estimator at spacecraft only needs to estimate its own state, and . The measurement update of FOEKF at spacecraft is as follows:
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 matrix, to absorb such errors. This leads to the increased extended Kalman filter (IREKF) [5].
3.3.2. IREKF
Considering the measurement , 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 to approximate . Thus the increase to the measurement noise of is given by .
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. 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.
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. 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.
4.1. Inference Engine
BE is the inference engine on JT. Let be the history measurements from initial time step to current time step , and let be the best estimation of current state under history measurements. If the conditional probability between new state and current state, , and the conditional probability between new measurement and new state, , 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.
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 ; from (12) we have , and from (13) we have . 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.
Figure 3 gives the DBN model of the example in Section 3.4, where the random variable 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, .
DBN holds the joint distribution of all the variables in . Inference on DBN means to find the joint distribution of interested variables, which is for this example. To perform the inference, we need to remove the barren variables in DBN, which are , , and . 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 is a subset of the full set . Each clique is assigned with a potential function, which is a probabilistic distribution of the variables in . The set of variables shared by two adjacent cliques and is called a separator, denoted by . 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 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.
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 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 , 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).
4.4.2. 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).
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 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 as the root of JT. The operations at each spacecraft return the following results.(i)Spacecraft 1:
Receiving Message.
Generating Message.(ii)Spacecraft 2:
Receiving Message.
Collecting Evidence.
Generating Message.(iii)Spacecraft 3:
Receiving Message.
Collecting Evidence.
Generating Message.
The message generated by spacecraft 3, , is the final result at time step , 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 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 right-hand 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.
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 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-mean-square 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 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.
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.
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 Section 3.1. Suppose and . 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 , , , , and , respectively. Thus the total computational complexity of EKF time update and measurement update is and , respectively.
The computational complexity of decentralized algorithms could be estimated using similar methods. For FOEKF, suppose . The computational complexity of FOEKF time update and measurement update at spacecraft is and , 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 and , respectively. Because the calculation of the increase to the measurement noise of is , 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 . 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.
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 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.
Appendices
A. State Transition Matrix and Measurement Matrix under Two-body Assumptions
A.1. State Transition Matrix
Recall the dynamic model in Section 2. Under two-body assumptions, the state transition matrix is given bywhere is the time step and is the gravity constant.
A.2. Measurement Matrix
Recall the measurement model in Section 2. While computing the measurement matrices, we need to consider the Keplerian-to-Cartesian partial derivatives, which are partials of positions to Keplerian orbital elements. Under two-body assumptions, the Keplerian-to-Cartesian partial derivatives are given bywhere is solved by the Kepler’s equation as follows:
B. Basic Operations on CG Potentials
B.1. Multiplication
Consider two random vectors, and , where and
Equation (B.1) defines the distribution of conditioned on , which is
The joint distribution of and is given by where
B.2. Division
If and are Gaussian,and the distribution of conditioned on is given bywhere and .
B.3. Marginalization
If the joint distribution of and is Gaussian, The distribution of each random vector is given by
Conflict of Interests
The authors declare no conflict of interests.
Authors’ Contribution
All authors contributed equally to this work.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant no. 61203200).