Globally Optimal Distributed Kalman Filtering for Multisensor Systems with Unknown Inputs

In this paper, the state estimation for dynamic system with unknown inputs modeled as an autoregressive AR (1) process is considered. We propose an optimal algorithm in mean square error sense by using difference method to eliminate the unknown inputs. Moreover, we consider the state estimation for multisensor dynamic systems with unknown inputs. It is proved that the distributed fused state estimate is equivalent to the centralized Kalman filtering using all sensor measurement; therefore, it achieves the best performance. The computation complexity of the traditional augmented state algorithm increases with the augmented state dimension. While, the new algorithm shows good performance with much less computations compared to that of the traditional augmented state algorithms. Moreover, numerical examples show that the performances of the traditional algorithms greatly depend on the initial value of the unknown inputs, if the estimation of initial value of the unknown input is largely biased, the performances of the traditional algorithms become quite worse. However, the new algorithm still works well because it is independent of the initial value of the unknown input.


Introduction
The classic Kalman filtering (KF) [1] requires the model of the dynamic system is accurate. However, in many realistic situations, the model may contain unknown inputs in process or measurement equations. The issue concerning estimating the state of a linear time-varying discrete time system with unknown inputs is widely studied by researchers. One widely adopted approach is to consider the unknown inputs as part of the system state and then estimate both of them. This leads to an augmented state Kalman filtering (ASKF). Its computational cost increases due to the augmented state dimension. It is proposed by Friedland [2] in 1969 a two-stage Kalman filtering (TSKF) to reduce the computation complexity of the ASKF, which is optimal for the situation of a constant unknown input. On the basis of the work in [2], it is proposed by Hsieh et al. an optimal two-stage algorithm (OTSKF) for the dynamic system with random bias and a robust two-stage algorithm for the dynamic system with unknown inputs in 1999 [3] and 2000 [4] respectively. It is assumed in [3][4][5] that the unknown inputs were an autoregressive AR (1) process, with the two-stage algorithms being optimal in the mean square error (MSE) sense. However, the optimality of the ASKF and OTSKF depends on the premise that the initial value of the unknown measurement can be estimated correctly. Under the condition of incorrect initial value of the unknown measurement, the ASKF and OTSKF will have poor performance (see Examples 1 and 2 in Section 5), especially, when the unknown measurement is not stationary as regarded in [4,5]. Due to the difficulty of knowing the exact initial value of the

Problem Formulation
Consider a discrete time dynamic system: where x k ∈ R m is the system state, ν k ∈ R m is the measurement vector, the process noise and measurement noise ω k ∈ R n are zero-mean white noise sequences with the following covariances: where: d k ∈ R p is the unknown input. Matrices F k , H k and A k are of appropriate dimensions by assuming that A k ∈ R n×p is of full column rank, i.e., rank(A k ) = p. Therefore, we have (A k ) † A k = I, where the superscript " †" denotes Moore-Penrose pseudo inverse. It is assumed d k follows an autoregressive AR (1), i.e.,: where B k is nonsingular and ω d k is a zero-mean white noise sequences with covariance: This model is widely considered in [2][3][4][5]. For example, in radar systems, the measurement often contains a fixed unknown deviation or an unknown deviation that gradually increases as the distance becomes longer. Such deviations can be described by Equation (6).
ASKF and OTSKF are two classic algorithms to handle this problem. The ASKF regards x k and d k as an augmented state and estimates them together, while the OTSKF estimates x k and d k respectively at first and then fusions them to achieve the optimal estimation. As a matter of fact, the unknown inputs can be eliminated easily by difference method. Denote: Equations (1) and (2) can be represented as: where: From Equation (12), it is not difficult to find out that the new measurement noise u k is one-step correlated and correlates with the process noise, i.e.,:

Optimal Estimation for the Remodeled System
It is assumed in the classic Kalman filtering that the process noises and measurement noises are uncorrelated temporally; in the meantime, both noises are mutually uncorrelated except at the same time instant. The noises in Equations (13)-(15) apparently violate these assumptions. Using the latest research achievements about Kalman filtering with correlated noises in [15][16][17][18][19][20], we can give an optimal estimation for the remodeled system (9) and (10) in the MSE sense. The recursive state estimate of the new system is presented in the following theorem. Theorem 1. The globally optimal estimate for the remodeled system (9) and (10) is given by: where: Remark 1. From Theorem 1, the new algorithm presented in this section is optimal in the MSE sense. In theory, the ASKF and OTSKF are also optimal in the MSE sense (see [2,3]). Nevertheless, the optimality of the ASKF and OTSKF depends on the assumption that the initial condition of the unknown measurement d 0|0 = E(d 0 ) , which is difficult to meet in real situations. It will be demonstrated by numerical examples in the Section 5 that if the initial value of the unknown input is wrong, their performances will be greatly influenced. By contrast, the new algorithm will continue its good performance as it does not rely on the initial value of the unknown input.

Remark 2.
A flop is defined as one addition, subtraction and multiplication. To estimate the complexity of an algorithm, the total number of flops is counted, expressing it as a polynomial of the dimensions of the matrices and vectors involved, and simplifying the expression by ignoring all terms except the leading terms. Then the complexities of the ASKF, OTSKF and the new algorithm are equivalent to O(m 3 + n 3 + p 3 + m 2 n + mn 2 + m 2 p + mp 2 + n 2 p + np 2 + mnp). The evaluation complexities of the three algorithms are the same order polynomials. We will compare their complexities more precisely by numerical examples in Section 5.

Multisensor Fusion
The l-sensor dynamic system is given by: where x k ∈ R m is the system state, y i k ∈ R n i is the measurement vector in the i-th sensor, ν k ∈ R m is the process noise and ω i k ∈ R n i is measurement noise, d i k ∈ R p i is the unknown input in i-th sensor. Matrices F k , H i k and A i k are of appropriate dimensions. We assume the system has the following statistical properties: (1) Every single sensor satisfies the assumption in Section 2.
Similarly to Equations (9) and (10), Equation (21) could be converted to: where: The stacked measurement equation is written as: According to Theorem 1, the local Kalman filtering at the i-th sensor is: with covariances of filtering error given by: where: According to Theorem 1, the centralized Kalman filtering with all sensor data is given by: The covariance of filtering error given by: where: diag is the diagonalization of a matrix.

Remark 3.
There are two key points to express the centralized filtering Equations (27) and (28) in terms of the local filtering: (1) Taking into consideration the measurement noise of single sensor in new system Equations (22) and (23), it can be observed that the sensor noises of the converted system are cross-correlated even if the original sensor noises are mutually independent. (2) ∆z k in Equation (27) is not stacked by local ∆z i k in Equation (26) directly and includes ∆z k−1 in its expression, which makes our problem more complicated than the previous distributed problem in [9][10][11][12][13][14]21].
Next, we will solve these two problems to express the centralized filtering Equation (28) in terms of the local filtering. We assume that H T k is of full column rank. Thus, we have (H T k ) † H T k = I. Using (28), we can get: Substituting (29) and (30) into (27), we have: Using (26), we have: We assume that J i k ∈ R m×n i is of full column rank, i.e., rankJ i k = n i . Thus, we have (J i k ) † J i k = I. Then, using (24), we can get: To express the centralized filtering x k|k in terms of the local filtering, by (25), (32) and (33), we have: where L † k ( * i) is the i-th column block of L † k . Thus, substituting (34) into (31) yields: Similarly to Equation (35), using Equations (24), (26), (29) and (32), we have: That means the centralized filtering is expressed in terms of the local filtering. Therefore, the distributed fused state estimate is equal to the centralized Kalman filtering adopting all sensor measurements, which means the distributed fused state estimate achieves the best performance.

Remark 4. From this new algorithm, it is easy to see that local sensors should transmit x i
k|k , x i k|k−1 , P i k|k and P i k|k−1 to the fusion center to get global fusion result. The augmented methods greatly increase the state dimension and computation complexity for the multisensor system. Since the difference method does not increase the state dimension, the computation complexity is lower than that of the augmented method for the multisensory system.

Numerical Examples
In this section, several simulations will be carried out for dynamic system with unknown inputs. It is assumed that the unknown input d k+1 = B k d k + ω k in this paper. Actually, the unknown measurement d k is a stationary time series if the eigenvalue of B k is less than 1, or else the unknown measurement d k is a non-stationary time series. The performances of the new algorithm (denoted as Difference KF) in these two cases are discussed in Example 1 and 2, respectively: Example 1. A two dimension target tracking problem is considered. The target dynamic models are given as Equations (1)- (7). The state transition matrices: and the measurement matrix is given by: Suppose A k is an identity matrix with appropriate dimensions, B k = 0.9I. In this case, d k is a stationary time series. The targets start at x 0 = (50, 1, 50, 1) T and the initial value d 0 = (5, 5) T . The covariance matrices of the noises are given by: In the following, the computer time and performances of the ASKF, OTSKF and Difference KF will be compared respectively.

•
Computer time The complexities of the three algorithms are analyzed in Remark 2, which shows the complexities of the three algorithms are the equivalent order polynomials. Now let us compare their computer time by this example. Table 1 illustrates the computer time of the three algorithms with 1000 Monte-Carlo runs respectively, through which we can find out that the new algorithm is the fastest algorithm in this example. • Estimation Performances In [3], Hsieh et al. has proved that the OTSKF is equivalent to the ASKF, so the tracking results of the two algorithms are the same. In order to make the figure clearer, we will only compare the performances of the following six algorithms: Algorithm 1: KF without considering unknown input. Algorithm 2: ASKF with accurate initial value of unknown input (d 0 = (5, 5) T ).
Using 100 Monte-Carlo runs, we can evaluate estimation performance of an algorithm by estimating the second moment of the tracking error: It must be noticed that Difference KF uses (y 1 , y 2 , · · · , y k , y k+1 ) to estimate x k at step k. However, the KF, ASKF and OTSKF only use (y 1 , y 2 , · · · , y k ) to estimate x k at step k. To make an equal comparison, x k|k−1 in Difference KF with x k|k in the other five algorithms is compared. As d k+1 = 0.9d k + ω k , d k is almost equal to a random white noise with small covariance after several steps and the influence of the initial value d 0 will be gradually weakened. The tracking errors of the six methods are compared in Figure 1 and Table 2. It can be noticed that no matter whether the initial values of the unknown input in ASKF and OTSKF are accurate or wrong, the tracking results of the six algorithms are almost the same after about 25 steps. However, it should be noticed that the Difference KF performs better than the ASKF with inaccurate initial value of unknown measurement in the first stage, which is important for some practical conditions, for instance, in multi-target tracking problems, due to data association errors and heavy clutters, tracking has to restart very often. Therefore, in order to derive an entirely good tracking, initial work status at each tracking restarting should be as good as possible.
of the six algorithms are almost the same after about 25 steps. However, it should be noticed that the Difference KF performs better than the ASKF with inaccurate initial value of unknown measurement in the first stage, which is important for some practical conditions, for instance, in multi-target tracking problems, due to data association errors and heavy clutters, tracking has to restart very often. Therefore, in order to derive an entirely good tracking, initial work status at each tracking restarting should be as good as possible.   Example 2. The dynamic equations are the same as Example 1. Assume B k = I. This model has been considered in [4,5]. d k is a non-stationary time series here. The non-stationary unknown measurement is common in practice. For instance, for an air target, the unknown radar bias is frequently increasing with distance changing between the target and radar.
The targets start at x 0 = (50, 1, 50, 1) T and the initial value d 0 = (5, 5) T .The performances of the following six algorithms are compared: Algorithm 1: KF without considering unknown input. Algorithm 2: ASKF with accurate initial value of unknown input (d 0 = (5, 5) T ). Algorithm 3: OTSKF with accurate initial value of unknown input (d 0 = (5, 5) T ). Algorithm 4: ASKF with wrong inaccurate initial value of unknown input (d 0 = (0, 0) T ). Algorithm 5: ASKF with inaccurate initial value of unknown input (d 0 = (20, 20) T ). Algorithm 6: Difference KF without any information about initial value of unknown input. Figure 2 and Table 3 compare the tracking errors of the six methods. As the new algorithm, ASKF and OTSKF with accurate initial value of unknown input are optimal in the MSE sense. Their performances are almost of no difference. The KF without considering unknown input is worse because it does not use any information of the unknown input. Numerical examples also demonstrate that once the initial value of the unknown input is inaccurate, the performance of the ASKF becomes poorer. We can also see that if the initial value of the unknown input is largely biased, the performance of ASKF is even poorer than KF ignoring unknown input. This is because d k+1 = d k + ω k in this example, the influence of the incorrect initial value d 0 will always exist. Nevertheless, the new algorithm is independent of the initial value of the unknown input and yet performs well.  From Examples 1 and 2, we can see that the performance of the difference KF is almost the same to that of the ASKF and OTSKF with accurate initial value of unknown input. If the initial value of the unknown measurement is largely biased, performances of the ASKF and OTSKF will be badly influenced. Due to the fact that it is not easy to get the exact initial value of the unknown measurement, Difference KF is a better option than the ASKF and OTSKF in practice.

Example 3. A two-sensor Kalman filtering fusion problem with unknown inputs is
The covariance matrices of the process noises is given by:  Algorithm From Examples 1 and 2, we can see that the performance of the difference KF is almost the same to that of the ASKF and OTSKF with accurate initial value of unknown input. If the initial value of the unknown measurement is largely biased, performances of the ASKF and OTSKF will be badly influenced. Due to the fact that it is not easy to get the exact initial value of the unknown measurement, Difference KF is a better option than the ASKF and OTSKF in practice.

Example 3.
A two-sensor Kalman filtering fusion problem with unknown inputs is considered. The object dynamics and measurement equation are modeled as follows: x k+1 = F k x k + ν k , k = 0, 1, . . . , 100 The state transition matrix F k and the measurement matrices H i k are the same as Example 1, A i k and B i k are identity matrix with appropriate dimensions. The targets start at x 0 = (50, 1, 50, 1) T and the initial value The covariance matrices of the process noises is given by: The covariance matrices of the measurement noises and the unknown inputs are diagonal given by The performances of the following three algorithms are compared as follows: Algorithm 1: Centralized KF without considering unknown input. Algorithm 2: The centralized fusion of the Difference KF. Algorithm 3: The distributed fusion of the Difference KF. The initial states of the three algorithms are set at x 0|0 = x 0 , the initial P i x 0|0 = I. Using 100 Monte-Carlo runs, we can evaluate estimation performance of an algorithm by estimating the second moment of the tracking error.
It is illustrated in Figure 3 and Table 4 that the simulation outcome of distributed fusion and centralized fusion of the new algorithm are exactly the same. Additionally, the new algorithm fusion gives better performance than the KF. Thus, the distributed algorithm has not only the global optimality, but also the good survivability in a poor situation.
The covariance matrices of the measurement noises and the unknown inputs are diagonal given by The performances of the following three algorithms are compared as follows: Algorithm 1: Centralized KF without considering unknown input. Algorithm 2: The centralized fusion of the Difference KF. Algorithm 3: The distributed fusion of the Difference KF.
The initial states of the three algorithms are set at 0|0 0 x x = , the initial 0|0 i x P I = . Using 100 Monte-Carlo runs, we can evaluate estimation performance of an algorithm by estimating the second moment of the tracking error.
It is illustrated in Figure 3 and Table 4 that the simulation outcome of distributed fusion and centralized fusion of the new algorithm are exactly the same. Additionally, the new algorithm fusion gives better performance than the KF. Thus, the distributed algorithm has not only the global optimality, but also the good survivability in a poor situation.

Conclusions
In this paper, the state estimation for dynamic system with unknown inputs modeled as an autoregressive AR (1) process is considered. The main contributions are: (1) A novel optimal algorithm for dynamic system with unknown inputs in the MSE sense is proposed by differential method. The computational burden of the new algorithm is lower than that of ASKF. The performance of the new algorithm is independent of the initial value of the unknown input. (2) A distributed fusion algorithm is proposed for the multisensor dynamic system with unknown inputs, the result of which is equal to the centralized difference Kalman filtering adopting all sensor measurements.
However, it should be noticed that the new algorithm uses y k and y k+1 to estimate x k , which leads the new algorithm to be one-step delayed. The new algorithm can only cope with the unknown inputs in measurement equation, while the ASKF can handle the unknown inputs in both state and measurement equation. Besides, it is assumed throughout the paper that the model of the unknown inputs d k follows an autoregressive AR (1) process. As for the future research, one interesting direction is to extend the difference method to dynamic system with more general unknown inputs.