Gaussian Filtering for Simultaneously Occurring Delayed and Missing Measurements

Approximate filtering algorithms in nonlinear systems assume Gaussian prior and predictive density and remain popular due to ease of implementation as well as acceptable performance. However, these algorithms are restricted by two major assumptions: they assume no missing or delayed measurements. However, practical measurements are frequently delayed and intermittently missing. In this paper, we introduce a new extension of the Gaussian filtering to handle the simultaneous occurrence of the delay in measurements and intermittently missing measurements. Our proposed algorithm uses a novel modified measurement model to incorporate the possibility of the delayed and intermittently missing measurements. Subsequently, it redesigns the traditional Gaussian filtering for the modified measurement model. Our algorithm is a generalized extension of the Gaussian filtering, which applies to any of the traditional Gaussian filters, such as the extended Kalman filter (EKF), unscented Kalman filter (UKF), and cubature Kalman filter (CKF). A further contribution of this paper is that we study the stochastic stability of the proposed method for its EKF-based formulation. We compared the performance of the proposed filtering method with the traditional Gaussian filtering (particularly the CKF) and three extensions of the traditional Gaussian filtering that are designed to handle the delayed and missing measurements individually or simultaneously.

The literature [18], [32], [33], [34] on filtering with simul-100 taneously occurring delayed and missing measurements wit-101 nesses some developments for different classes of systems. 102 For example, [32] considered this problem in coupled neu-103 ral networks and developed an estimation method. Simi-104 larly, [18] develops a filtering algorithm for the delayed and 105 missing measurements. In this regard, it introduces sepa-106 rate stochastic models for incorporating the delay and miss-107 ing measurements possibilities. Subsequently, it introduces a 108 Ricatti-like equation for designing Kalman filtering method. 109 However, [18] restricts the delay up to one sampling interval, 110 wherein the practical delays can often be larger. Further-111 more, [33] introduced unbiased finite impulse response-based 112 filtering approach for finite-horizon case, considering the 113 presence of delayed and missing measurements. However, 114 it assumes that the delay is time-stamped, while delays 115 without time-stamping are observed in many practical sys-116 tems [35], [36]. Finally, [34] introduces three new designs 117 of robust linear Kalman filtering for handling the simulta-118 neously occurring delayed and missing measurements. How-119 ever, in every design, it restricts the delay as one sampling 120 interval, which can lead to poor accuracy if the real mea-121 surement delay is higher. Moreover, [18], [34] rely on the 122 augmented state-space approaches, which may increase the 123 computational complexity significantly. 124 Although the above-discussed filtering methods 125 [18], [33], [34] can handle the simultaneously occurring 126 delayed and missing measurements, they have two broad 127 limitations: i) they are designed for linear dynamical systems, 128 ii) they fail to handle delays of more than one sampling 129 interval without time-stamping. For the nonlinear dynami-130 cal systems, [37] is probably the only existing method to 131 handle the simultaneously occurring delayed and missing 132 measurements. It integrates a likelihood-based technique 133 with the nonlinear Gaussian filtering framework for handling 134 the simultaneously occurring delayed and missing measure-135 ments. However, [37] uses an augmented state transition 136 equation, with the size of the system growing by a factor 137 of d, where d the maximum number of delays. The authors 138 also use Gaussian mixture with d components as a likeli-139 hood function, and mention computational issues with large 140 number of components in the likelihood function. Indeed, the 141 increased size of the covariance matrices itself can present 142 significant computational issues in multidimensional integra-143 tion involved. 144 In this paper, we develop a new extension of Gaussian 145 filtering to handle the simultaneously occurring delayed 146 and missing measurements for nonlinear dynamical sys-147 tems. In this regard, we reformulate the measurement model 148 using a set of Bernoulli random variables to incorporate 149 the possibilities of delayed and missing measurements. 150 Subsequently, we re-derive the Gaussian filtering method 151 for the modified measurement model. It should be men-152 tioned that we modify only the measurement model.  sequently, in the redesigned Gaussian filtering, only the 154 filtering parameters related to the measurements, such as 155 VOLUME 10, 2022 the measurement estimate, the measurement covariance, and 156 the state-measurement cross-covariance, are re-derived. It is In view of the above discussion, we highlight the main

196
The remaining part of the paper is organized as follows.

197
In Section II, we mathematically formulate the problem where x k ∈ R n and z k ∈ R q are state and measure-214 ment variables, respectively, at k th sampling instant, k ∈ 215 {1, 2, . . . , T s } with T s representing the number of sampling 216 intervals. Moreover, f k : x k−1 → x k and h k : x k → z k are 217 general nonlinear functions. Finally, η k and ν k are zero-mean 218 Gaussian noises representing the process and measurement 219 noises, respectively. The covariances of η k and ν k are denoted 220 as Q k and R k , respectively.

221
Following our problem statement, we need to reformulate 222 the measurement model (Eq. (2)) to address the simultaneous 223 occurrence of randomly delayed and randomly missing mea-224 surements. Our reformulation of the measurement model is 225 based on two sets of Bernoulli random variables, denoted by 226 α and : α corresponds to the missing measurements and 227 corresponds to the delayed measurements.

228
The measurements are generally received from multiple 229 sources, and they all may not be missing at the same time. 230 Thus, we consider that the measurement at any sampling 231 instant may be partly missing, i.e., specific elements of 232 the measurement may be missing at any particular instant. 233 Thus, we define a matrix of Bernoulli random variables, 234 It should be mentioned that α i k is either 237 0 or 1, with α i k = 0 representing that the i th element of the 238 received measurement y k , denoted as y k (i), is missing.

239
For modeling the delay portion, we restrict the maximum 240 delay to N max . N max is the practitioner's choice and it can 241 be assigned with a fairly large value if the expected delay 242 is large. Therefore, our model and the proposed filtering 243 technique should not be deemed to be restricted to small 244 delays. We define N max + 1 equiprobable Bernoulli random 245 variables: one for each of the current and the N max possible 246 delayed instants. At k th instant, we denote them as Note that j+1 k corresponds to j th delayed instant. We assign 249 0 k = 0, and model the actual measurement as The coefficients of z k−m ∀m ∈ {1, 2, . . . , N max } gov-254 ern the delay extent. For example, if the measurement 255 is one time-step delayed, i.e., y k = z k−1 , then coeffi-256 cient of z k−1 , i.e., (1 − 0 k )(1 − 1 k ) 2 k takes the value 257 one, while the random variables associated with z k−m 258 ∀m = 1 remain zero. At the same time, λ k regulates the 259 missing measurement possibility. The diagonal elements of 260 FIGURE 1. Pictorial diagram representing the sequence to be followed to obtain the received measurement y k from the ideal z 1 , z 2 , . . . , z k that would have been received in the lenient environment.  To this end, let us simplify the notation for the coefficients Thus, the measurement model can be finally given as

287
The above discussions emphasize the importance of the 288 measurement model (6) for developing the proposed fil-289 tering algorithm. As mentioned in the previous section, 290 some of the existing filters, such as [19], [20] The traditional Gaussian filtering is designed with respect 313 to the measurement z, modeled in Eq. (2). In this section, 314 we derive the necessary modifications to the algorithm to deal 315 with the modified measurement y, modeled in Eq. (6). As the 316 measurement model is changed, we re-derive all the related 317 expressions in the traditional Gaussian filtering to propose the 318 advanced Gaussian filtering for y. The traditional Gaussian 319 filtering uses only three such expressions, namely the mea-320 surement estimateẑ, measurement error covariance P zz , and 321 the cross covariance P xz , derived for z. We re-derive all the 322 measurement related expressions in the Gaussian filtering 323 algorithm for the modified measurement model above. On a 324 different note, it should be mentioned that the state dynam-325 ics remains unaffected from the simultaneous occurrence of 326 the delayed and missing measurements. Therefore, the time 327 update step of the proposed filtering technique remains the 328 same as the traditional Gaussian filtering [9], [10].

330
In this part, we re-derive the measurement parameters, such as 331 y, P yy , and P xy . Before proceeding to the derivation, it should 332 VOLUME 10, 2022 be mentioned that only one of k (m, j) ∀m is one and the 333 others are zero at any instant t k to ensure that only one mea-334 surement is received. Although this consideration violates 335 the independence of k (m, j) for different m (i.e., different 336 delay), they will be assumed to be statistically independent 337 in our derivation. We now derive the expressions ofŷ, P yy , 338 and P xy . For y k given in Eq. (6), the measurement estimate is 342 As the missing and delay occurrences are mutually indepen-343 dent events, λ k and k (m, j) are statistically independent.

344
Moreover, λ k and k (m, j) are independent of the measure-345 ment value z k also. Thus, we simplify the above equation as

347
Following our previous notations,

348
Recalling the previous discussion, we get The measurement error covariance is

357
From Eqs. (6) and (10), we get We can rewrite this expression as

364
From Eqs. (11) and (13), we can write We can now compute every expectation term individually 368 for A 1 and A 2 defined in Eq. (13) and add them to obtain 369 P yy k|k−1 . In this regard, for A 1 given in Eq. (13), we get Similarly, for A 1 and A 2 given in Eq. (13), we can write After further simplification, we get After substituting all the expectation terms from previous 389 discussions and considering that This also leads to Finally, for A 2 given in Eq. (13), we have Applying binomial expansion and simplifying further, we get The cross-covariance between the state and measurement is

413
Substituting y k −ŷ k|k−1 from Eq. (13), we get For A 1 given in Eq. (13), we get Similarly, for A 2 given in Eq. (13), we get Substituting this into Eq. (28), we get from Eqs. (27) and (30), respectively, into Eq. (25), we get As discussed at the beginning of this section, the pro-435 posed filtering method modifies the traditional Gaussian 436 filtering by re-deriving the expressions of measurement 437 estimate, measurement covariance, and state-measurement 438 cross-covariance (Eqs. (10), (23), and (31), respectively). 439 Please follow [9], [10] for a detailed discussion on the 440 traditional Gaussian filtering. The proposed filtering algo-441 rithm also follows the same filtering strategy by replac-442 ing the expressions ofẑ, P zz , and P xz with the re-derived 443 expressions ofŷ, P yy , and P xy , respectively. We provide the 444 pseudo-code for implementing the proposed filtering method 445 in Algorithm 1.

446
In advancing the traditional Gaussian filtering for handling 447 various measurement irregularities, such as the delayed and 448 missing measurements, the major difficulty appears in incor-449 porating those irregularities through mathematical models. 450 The problem becomes yet more challenging if the irregu-451 larities are uncertain to appear at any particular sampling 452 instant, as considered in this paper. We handled this problem 453 by mathematically characterizing such irregularities, particu-454 larly the delayed and missing measurements, by formulating 455 a stochastic model, as in Eq. (6).

456
Remark 2: Our measurement model utilizes a sequence of 457 Bernoulli random variables to characterize the possibility of 458 a measurement coming from various possible past instants. 459 A future research problem may be to introduce a more con-460 venient model by reducing the required number of random 461 variables.

462
Remark 3: Our filter design strategy concludes that han-463 dling the measurement irregularities becomes convenient if 464 an efficient mathematical model for characterizing the con-465 cerned irregularities is formulated.

564
To this end, it should be mentioned that the literature 565 contains several notions of stability for nonlinear systems.

566
The readers may please refer to [39] for a detailed discussion.

567
In this paper, we particularly use exponential stability, where 568 the stability is ensured if the system's convergence is bounded 569 with an exponential envelope.

571
Before proceeding forward to prove that the error dynamics 572 presented in Eq. (38) is exponentially bounded, we mathemat-573 ically define the exponentially bounded process, as follows. 574 Definition 1: Let us consider that e k denotes a stochastic 575 process and κ > 0, ξ > 0, and 0 < β < 1 are real numbers. 576 Then, e k is said to be exponentially bounded in mean square 577 if it satisfies where · represents the spectral norm for matrices and 580 Euclidean norm for vectors.

581
It should be mentioned that the above definition is general 582 and does not limit e k to be the estimation error only. However, 583 our notations and discussions will be focused on the estima-584 tion error.

585
In our stability analysis, we approach Eq. (39) in a different 586 way. In this regard, please refer to the subsequent discussion. 587 Remark 6: Let us consider that τ 1 > 0, τ 2 > 0, γ > 0, 588 and 0 < φ < 1 denote real numbers, and V (e k ) represents a 589 scalar-valued stochastic process, which satisfies Then, the stochastic process e k satisfies Please refer to [38] and [40], for a detailed discussion.

596
In the subsequent discussion, we conclude that this remark 597 is another way of defining the exponential bound in mean 598 square.
VOLUME 10, 2022 It should be mentioned that with τ 2 /τ 1 = κ , 1 − φ = β, and process V (e k ), then the stochastic process e k is exponentially 611 bounded. 612 We will use Remark 8 as our stability criterion. Alter-613 natively, we will conclude the stability of the MDEKF by

617
• F k is a non-singular matrix.

618
• The system, noise, and filter parameters satisfy the fol-619 lowing bounds: It is worth mentioning that H k P k−m H T k ≥ 0 because 649 P k−m|k−1 is a positive definite matrix. Subsequently, applying 650 the norm property, we can write For any invertible matrix M, please note that σ + (M −1 ) = 655 (σ − (M)) −1 , where σ + (·) and σ − (·) represent the largest 656 and smallest singular values, respectively. Thus, substituting 657 µ k = µ and using the inequalities presented in Eqs. (47) 658 and (48), we get The bound presented in Eq. (52) assumes that R k is a posi-663 tive definite matrix. Moreover, h(x k−m|k−1 )h(x k−m|k−1 ) T is a 664 positive semidefinite matrix. Thus, the matrix in the second 665 factor of the above equation is also positive definite. Subse-666 quently, it follows that σ − (·) = λ − (·), with λ − (·) represent-667 ing the smallest eigenvalue. We now calculate the smallest 668 eigenvalue of the second factor by using the Rayleigh-Ritz 669 characterization [41] in Eq. (55), as shown at the bottom of 670 the page.
second term on the right side of Eq. (55) is zero. Subse-673 quently, using the bound given in Eq. (52), we obtain It is worth mentioning thatẑ k−m|k−1 , P zz k−m|k−1 , and P xz can be determined using Eq. (33). Substituting these param-  Please note that R k is a positive definite matrix in Eq. (58). 701 Similarly, the positive definiteness of P k−m|k−1 ensures that 702 T is a posi-704 tive semidefinite matrix. Since the sum of the positive definite 705 and positive semidefinite matrices is positive definite, we fur-706 ther conclude Consequently, the following inequality can be deduced from 712 Eq. (58): Substituting Eq. (32), we can rearrange the above inequality 716 as Applying the bounds of Eqs. (47), (48), and (52), and taking 722 the inverse, the above inequality can be rearranged as We choose q, f , and ρ 2 for 726 which 0 Proof: Throughout the theorem, we adopt the fol-

750
It should be mentioned that the above equation is scalar. Thus, 751 applying norm property to the second expression on the right 752 side, we get We will now calculate the bound of each term in the right side 756 individually. In this regard, let us recall B k defined in Eq. (38) 757 and use norm property to get Similarly, for A k given in Eq. (38), we get Applying the bound given in Eq. (48), we have 774 P −1 k|k ≤ 1/ρ 1 . Thus, substituting Eqs. (49), (64), and (65) 775 into Eq. (63), it is simplified as Let us now consider the fourth expression on the right side of 782 Eq. (62). For C k defined in Eq. (38), applying C T k = C k 783 gives 787 We now substitute µ k = µ, λ k ≤ 1, N max m=0 k (m, j) = 788 1, and P −1 k|k ≤ 1/ρ 1 . Then, using the inequalities presented 789 in Eqs. (44) and (47), and substituting the bound of K from 790 Lemma 1, we finally get

794
(68) 795 For D k defined in Eq. (38), applying the matrix norm 796 property, the fifth expression on the right side of Eq. (62) 797 can be expressed as We now substitute µ k = µ, λ k ≤ 1, and Let us now consider the sixth term in the summation on 816 the right side of Eq. (62). Applying the norm property,

820
where χ 4 is a positive real number given as 821 Then, using the 832 definition of V (e k ) and substituting χ 1 3 + χ 2 + χ 3 + χ 4 = 833 γ , Eq. (74) can be expressed as We select the parameters γ and φ such that V (e k−1|k−1 ) ≥ 837 γ /φ satisfies. Then, the above equation is the same as 838 Eq. (41). 839 We now consider the inequality presented in Eq. (48). 840 Taking inverse, then multiplying e T k|k and e k|k from left and 841 right side, respectively, we get In this section, we proved the stochastic stability of the 851 proposed method for its EKF-based formulation. The essen-852 tial requirements of the proof involve a set of bounds on the 853 system, noise, and filter parameters. Moreover, the stability 854 analysis requires the initial estimation error to be bounded. 855 Note that boundedness of noise and system parameters does 856 not automatically imply exponential stability. The stability 857 analysis of the proposed MDEKF reduces to the stability 858 analysis of the traditional EKF for the delay and missing 859 measurements probabilities being zero (µ i = δ = 1) and 860 N max = 0, considering the convention 0 0 = 1.

862
In real-life problems, the measurement systems (including the 863 measuring devices and the supplementary units) are usually 864 designed to efficiently capture the measurements. Therefore, 865 they may not be expected to miss many measurements. Sub-866 sequently, the missing measurement probability is usually 867 small. Thus, we consider the missing measurement prob-868 ability up to 0.2 for the simulation, which means around 869 20% of the measurements are missing. On the other hand, 870 the delay inherently appears in the measurements. Therefore, 871 we consider a sufficiently large range of the delay probability 872 (0.1 ≤ 1 − δ d ≤ 0.9). Moreover, it should be mentioned 873 that the practical measurement systems are designed for small 874 delays. Hence, we restrict the maximum possible delay to one 875 or two time-steps (denoted as 1-delay and 2-delay scenarios), 876 such as the delay up to one or two sampling intervals. 877 VOLUME 10, 2022  For the performance analysis, we considered three pop-878 ular and advanced Gaussian filters, namely, the CKF [13],

881
MDCQKF, and MDGHF, respectively. We use the root mean 882 square error (RMSE) as our performance metrics. Please note 883 that we will frequently use the notation µ m = 1−µ to denote 884 the missing measurements probability. 885 We compare the MDCKF with the following filters: 886 i) traditional CKF [13], ii) the CKF-based formulation 887 of [22] In the first problem, we consider a two-dimensional nonlinear 900 dynamical system with the state-space model given as x k = 901 2 cos(x k−1 ) + η k and z k = 1 + x T k x k + ν k [22].

902
The true data of the state and measurement are gener-903 ated by considering the initial state as The filter is initialized with the initial estimatex 0|0 = 0.9x 0 905 and P 0|0 = 7I 2 . The noise covariances are assigned as 906 Q k = 0.1I 2 and R k = 0. In the second simulation problem, we consider an identifica-927 tion problem of individual sinusoids from the measurements 928 of multiple superimposed sinusoids [9], [22]. We consider 929 that the superimposed signal consists of three sinusoids. Iden-930 tification of individual sinusoids is equivalent to estimating 931      (Table 3). This research did not use any experimentally generated data 1004 or data from any publicly available dataset. Model def-1005 initions (including the specified probability distributions) 1006 and parameter values (including the initialization parame-1007 ters) provided in the paper are adequate for reproducing the 1008 exact qualitative behavior of the algorithms illustrated in the 1009 paper.