Stochastic H∞ Filtering of the Attitude Quaternion

1. Introduction

The attitude quaternion [1][p. 758] is a very popular spacecraft attitude parametrization, whose mathematical modeling and filtering have been ongoing topics of research for more than four decades [2]. Numerous successful quaternion stochastic estimators have been developed in the realm of optimal filtering, e.g. [3,4] and [5, Ch. 6] for an in-depth survey. An inherent drawback of the optimal filtering approach consists in its sensitivity to the model noise variance. Although adapting the filter gain computations might provide satisfactory performances in some cases, like adaptive process noise estimation in [6] and uncompensated random biases in [7,8], the designer might prefer a less sensitive approach: rather than trying to estimate the unknown parameters, the filter shall attenuate their impact on the estimation error for a transmission level that should be as small as possible. This general approach was followed in [9,10,11] for spacecraft attitude determination. In [9] an $H_{\infty }$ estimator of the spacecraft quaternion and gyro biases was developed. In [10,11] extended $H_{\infty }$ filters were applied to spacecraft attitude determination and gyro calibration by processing space-borne telemetry of the CBERS-2 satellite, including outputs from rate gyros, two sun sensors, and an earth’s sensor. The filters produced estimates of Euler angles and the quaternion, respectively, along with estimates of the gyro biases. Their performance compared favorably with those of standard extended Kalman filters. Reference [12] presented an invariant extended $H_{\infty }$ filter for attitude and position estimation using Lie groups algebra, which conveniently produces state-independent Jacobians at any linearization point. All the above works are rooted in the deterministic $H_{\infty }$ estimation theory. In that realm, the measurement and process noises, and the initial errors, are modeled as energy-finite disturbances, a.k.a. ${\mathcal{L}}_2$ or $l_2$ signals.

In this work, stochastic $H_{\infty }$ quaternion filters are developed for models with process and measurement white noise rather than finite-energy disturbances. We assume that a time-varying line-of-sight (LOS) signal is continuously acquired, that a triad of body-mounted gyros provides a measurement of the spacecraft’s angular velocity, and that both measurement processes are corrupted by additive white noise. The case of a gyro bias is addressed through Brownian motion modeling. All white noise variances are assumed to be unknown. A distinctive feature of our modeling and filtering approach is that the plant equations admit white noise as inputs and that their variances are assumed to be finite-energy perturbations. The stochastic filter aims at estimating the quaternion while attenuating these perturbations. This article is a revised and expanded version of a paper [16] which was presented at the AIAA Guidance Navigation and Control Conference, Toronto, Ontario, Canada, 2010. It is here extended to include the development of an $H_{\infty }$ filter for quaternion and gyro bias estimation, along with extensive Monte-Carlo simulations for validation. The proposed approach is inspired by works on stochastic $H_{\infty }$ filtering and control for nonlinear systems, which are grounded in dissipativity theory [13,14].

Several filters are developed: first, the gyro noise variance is the sole unknown, then both the gyro and the LOS noise variances are unknown, and finally the case of biased gyro measurements is considered. The state multiplicative nature of the errors in both the process and measurement provides a useful structure to the model. The proposed approach avoids linearization by exploiting a pseudo-linear quaternion plant model, which was introduced in [6] and extended to continuous-time stochastic settings in [15]. The filter implementation requires solving differential linear matrix inequalities that do not depend on the state estimate. Extensive Monte-Carlo simulations were run in order to evaluate the novel filter’s performances as an attitude estimator, and to compare them with those of a quaternion multiplicative extended Kalman filter.

Section 2 presents the quaternion $H_{\infty }$ filters for the case of drift-free rate gyros. Section 3 is concerned with the non-zero drift case. Section 4 presents the results of Monte-Carlo simulations. Conclusions are proposed in the last section.

2. Quaternion Stochastic $H_{\infty }$ Filtering

2.1. Problem statement

For the sake of simplicity and clarity, the case of drift-free gyroscopes is treated first. Consider the following Itô stochastic differential equation (SDE) for the attitude quaternion, and the associated measurement equation [15]:

$$\begin{array}{cc}\hfill & \mathbf{d}{\mathbf{q}}_{\mathbf{t}}={\displaystyle \frac{1}{2}}({\mathsf{\Omega }}_t-{\displaystyle \frac{3{\sigma }_{ϵ}^2}{4}}{I}_4){\mathbf{q}}_{t}dt-{\displaystyle \frac{\mathtt{1}}{\mathtt{2}}}\mathsf{\Xi }\left({\mathbf{q}}_t\right){\sigma }_{ϵ}{\mathbf{d}\beta }_t;{\mathbf{q}}_t\left(0\right)\stackrel{a.e.}{=}{\mathbf{q}}_0;t\in [0,T]\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & \mathbf{d}{\mathbf{y}}_{\mathbf{t}}=H_t{\mathbf{q}}_{t}dt-{\displaystyle \frac{1}{2}}\mathsf{\Xi }\left({\mathbf{q}}_t\right){\sigma }_b{\mathbf{d}\mathit{\eta }}_t\hfill \end{array}$$

(2)

where ${\mathbf{q}}_t$ denotes the attitude quaternion, ${\mathsf{\Omega }}_t$ is the following matrix function of the measured angular velocity ${\omega }_t$,

$$\begin{array}{cc}\hfill & {\mathsf{\Omega }}_t=\left[\begin{array}{cc}-\left[{\omega }_t\times \right]& {\omega }_t\\ -{\omega }_t^T& 0\end{array}\right]\hfill \end{array}$$

(3)

where ${\omega }_t$, which is acquired by a triad of body-mounted gyroscopes, is corrupted by an additive standard Brownian motion, ${\beta }_t$, with infinitesimal independent increments ${\mathbf{d}\beta }_t$ such that $E\left\{{\mathbf{d}\beta }_t{\mathbf{d}\beta }_t^T\right\}=I_{3}dt$. The parameter ${\sigma }_{ϵ}$ denotes the standard deviation of the gyro noise ${\beta }_t$. The drift term in Eq. (1) features a damping term in ${\sigma }_{ϵ}$ that ensures mean-square stability of the process ${\mathbf{q}}_t$. Equation (2) is a continuous-time quaternion measurement equation where the measurement value is identically zero [6]. Hence,

$$\begin{array}{c}\hfill \mathbf{d}{\mathbf{y}}_{\mathbf{t}}=\mathbf{0}\end{array}$$

(4)

and the measurement matrix $H_t$ is constructed from LOS measurements. Let ${\mathbf{b}}_t$ and ${\mathbf{r}}_t$ denote the projections of a measured LOS in the spacecraft body frame axes and a reference inertial frame, respectively, then $H_t$ is computed as follows

$$\begin{array}{cc}\hfill & {\mathbf{s}}_t={\displaystyle \frac{1}{2}}({\mathbf{b}}_t+{\mathbf{r}}_t)\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & {\mathbf{d}}_t={\displaystyle \frac{1}{2}}({\mathbf{b}}_t-{\mathbf{r}}_t)\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & H_t=\left[\begin{array}{cc}-\left[{\mathbf{s}}_t\times \right]& {\mathbf{d}}_t\\ -{\mathbf{d}}_t^T& 0\end{array}\right]\hfill \end{array}$$

(7)

The matrix, $\mathsf{\Xi }\left({\mathbf{q}}_t\right)$, which appears both in the process and measurement multiplicative noises, is the following linear matrix function of the quaternion ${\mathbf{q}}_t={\left[{\mathbf{e}}_t^T{q}_t\right]}^T$:

$$\mathsf{\Xi }\left({\mathbf{q}}_t\right)=\left[\begin{array}{c}{q}_t{I}_3+\left[{\mathbf{e}}_t\times \right]\\ -{\mathbf{e}}_t^T\end{array}\right]$$

(8)

The measurement noise is modeled as a standard Brownian motion, ${\mathit{\eta }}_t$, multiplied by the parameter ${\sigma }_b$. Assume that ${\sigma }_{ϵ}$ and ${\sigma }_b$ are unknown, possibly time-varying random parameters, and that there is no prior knowledge of their possible values or bounds. The filtering problem consists of estimating the quaternion ${\mathbf{q}}_t$ in the presence of unknown and random${\sigma }_{ϵ}$ and ${\sigma }_b$ and is formulated as a disturbance attenuation problem. The following estimator is considered:

$$\begin{array}{cc}\hfill & \mathbf{d}{\widehat{\mathbf{q}}}_{\mathbf{t}}={\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t{\widehat{\mathbf{q}}}_{t}dt+K\left({\widehat{\mathbf{q}}}_t\right)\left(\mathbf{d}{\mathbf{y}}_{\mathbf{t}}-H_t{\widehat{\mathbf{q}}}_{t}dt\right)\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & \widehat{\mathbf{q}}\left(0\right)={\widehat{\mathbf{q}}}_0\hfill \end{array}$$

(10)

The Itô correction term, ${\displaystyle \frac{3{\sigma }_{ϵ}^2}{4}}$, is dropped for simplicity. It is not required for the filter development since the dynamics of the filter estimate and estimation error, $\{{\widehat{\mathbf{q}}}_t,{\tilde{\mathbf{q}}}_t\}$, are designed to be stable in the mean-square sense independently from that correction [13]. Let ${\tilde{\mathbf{q}}}_t$ denote the additive quaternion estimation error, i.e.,

$$\begin{array}{cc}\hfill & {\tilde{\mathbf{q}}}_t={\mathbf{q}}_t-{\widehat{\mathbf{q}}}_t\hfill \end{array}$$

(11)

Given a scalar $\gamma >0$, a gain process $K\left({\widehat{\mathbf{q}}}_t\right)$ is sought such that the following $H_{\infty }$ criterion is satisfied:

$$\begin{array}{c}\hfill E\{{\int }_0^T\parallel {\tilde{\mathbf{q}}}_t{\parallel }^{2}dt\}\le {\gamma }^{2}E\{\parallel {\tilde{\mathbf{q}}}_0{\parallel }^2+{\int }_0^T\parallel {\mathbf{v}}_t{\parallel }^{2}dt\}\end{array}$$

(12)

under the constraints (1)(2), and where ${\mathbf{v}}_t$ denotes the augmented process of admissible disturbance functions. Whenever Eq. (12) is true, the so-called $L_2$-gain property from $\{{\tilde{\mathbf{q}}}_0,{\mathbf{v}}_t\}$ to ${\tilde{\mathbf{q}}}_t$, is satisfied for $0\le t\le T$. Two cases for $\mathbf{v}$ will be considered in the following: $\mathbf{v}={\sigma }_{ϵ}$ and $\mathbf{v}=\{{\sigma }_{ϵ},{\sigma }_b\}$ and significant differences will be highlighted.

2.2. Augmented stochastic process

Following standard techniques, the following augmented process is defined:

$$\begin{array}{cc}\hfill & {\mathbf{q}}_t^a\stackrel{▵}{=}\left[\begin{array}{c}{\widehat{\mathbf{q}}}_t\\ {\tilde{\mathbf{q}}}_t\end{array}\right]\hfill \end{array}$$

(13)

The design model for the process ${\mathbf{q}}_t$ obeys the following equation

$$\begin{array}{cc}\hfill & \mathbf{d}{\mathbf{q}}_{\mathbf{t}}={\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t{\mathbf{q}}_{t}dt-{\displaystyle \frac{\mathtt{1}}{\mathtt{2}}}\mathsf{\Xi }\left({\mathbf{q}}_t\right){\sigma }_{ϵ}{\mathbf{d}\beta }_t\hfill \end{array}$$

(14)

The estimator’s equation is given by Eq. (9), i.e.

$$\begin{array}{cc}\hfill \mathbf{d}{\widehat{\mathbf{q}}}_{\mathbf{t}}& ={\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t{\widehat{\mathbf{q}}}_{t}dt+{\widehat{K}}_t\left(\mathbf{d}{\mathbf{y}}_{\mathbf{t}}-H_t{\widehat{\mathbf{q}}}_{t}dt\right)\hfill \\ \hfill & =\left[{\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t-{\widehat{K}}_t{H}_t\right]{\widehat{\mathbf{q}}}_{t}dt\hfill \end{array}$$

(15)

where ${\widehat{K}}_t$ denotes $K\left({\widehat{\mathbf{q}}}_t\right)$. Using Eqs. (14), (15), and (2) yields the SDE for the estimation error:

$$\begin{array}{cc}\hfill \mathbf{d}{\tilde{\mathbf{q}}}_{\mathbf{t}}& =\left[{\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t-{\widehat{K}}_t{H}_t\right]{\tilde{\mathbf{q}}}_{t}dt-{\displaystyle \frac{1}{2}}\mathsf{\Xi }{\sigma }_{ϵ}{\mathbf{d}\beta }_t+{\displaystyle \frac{1}{2}}{\widehat{K}}_t\mathsf{\Xi }{\sigma }_b{\mathbf{d}\mathit{\eta }}_t\hfill \\ \hfill & =\left[{\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t-{\widehat{K}}_t{H}_t\right]{\tilde{\mathbf{q}}}_{t}dt-{\displaystyle \frac{1}{2}}\sum _{i=1}^3{\mathsf{\Xi }}_{Ci}{\sigma }_{ϵ}d{\beta }_i+{\displaystyle \frac{1}{2}}{\widehat{K}}_t\mathsf{\Xi }{\sigma }_b{\mathbf{d}\mathit{\eta }}_t\hfill \end{array}$$

(16)

where $\mathsf{\Xi }$ denotes the matrix $\mathsf{\Xi }\left(\mathbf{q}\right)$, ${\mathsf{\Xi }}_{Ci}$, $i=1,2,3$, denote the columns of the matrix $\mathsf{\Xi }$, and the scalar processes ${\beta }_i$, $i=1,2,3$, are the components of the vector Brownian motion ${\mathit{\beta }}_i$. Notice that $\mathsf{\Xi }$ is a function of the augmented process $\{{\widehat{\mathbf{q}}}_t,{\tilde{\mathbf{q}}}_t\}$ since it is a function of the quaternion $\mathbf{q}$. ${\widehat{\mathsf{\Xi }}}_t$ and ${\tilde{\mathsf{\Xi }}}_t$ denote the matrices $\mathsf{\Xi }\left({\widehat{\mathbf{q}}}_t\right)$ and $\mathsf{\Xi }\left({\tilde{\mathbf{q}}}_t\right)$, respectively, while ${\widehat{\mathsf{\Xi }}}_{Ci}$ and ${\tilde{\mathsf{\Xi }}}_{Ci}$, $i=1,2,3$, denote their columns. Appending Eqs. (15) and (16) yields the following augmented SDE:

$$\begin{array}{cc}\hfill & \left[\begin{array}{c}\mathbf{d}{\widehat{\mathbf{q}}}_{\mathbf{t}}\\ \mathbf{d}{\tilde{\mathbf{q}}}_{\mathbf{t}}\end{array}\right]=\left[\begin{array}{c}({\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t-{\widehat{K}}_t{H}_t){\widehat{\mathbf{q}}}_t\\ ({\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t-{\widehat{K}}_t{H}_t){\tilde{\mathbf{q}}}_t\end{array}\right]dt+\sum _{i=1}^3\left[\begin{array}{c}{O}_{4\times 1}\\ -{\displaystyle \frac{1}{2}}{\mathsf{\Xi }}_{Ci}\end{array}\right]{\sigma }_{ϵ}d{\beta }_i+\left[\begin{array}{c}{O}_{4\times 3}\\ {\displaystyle \frac{1}{2}}{\widehat{K}}_t\mathsf{\Xi }\end{array}\right]{\sigma }_b{\mathbf{d}\mathit{\eta }}_t\hfill \end{array}$$

(17)

where $O_{m\times n}$ denotes an $m\times n$ matrix of zeros. Equation (17) can be re-written in the following compact form:

$$\begin{array}{c}\hfill \mathbf{d}{\mathbf{q}}_t^a=F^a{\mathbf{q}}_t^{a}dt+\sum _{i=1}^3{\mathbf{g}}_2^i\left({\mathbf{q}}_t^a\right){\sigma }_{ϵ}d{\beta }_i+G\left({\mathbf{q}}_t^a\right){\sigma }_b{\mathbf{d}\mathit{\eta }}_t\end{array}$$

(18)

where $F^a$ , ${\mathbf{g}}_2^i\left({\mathbf{q}}_t^a\right)$ and $G\left({\mathbf{q}}_t^a\right)$ are effectively defined from Eq. (17).

2.3. Hamilton-Jacobi-Bellman inequality

The desired $L_2$-gain property will be satisfied if and only if the augmented system (18) is dissipative with respect to the supply rate $S(\mathbf{v}\left(t\right),{\mathbf{q}}_t^a)={\gamma }^2{\parallel \mathbf{v}\left(t\right)\parallel }^2-{\parallel {\tilde{\mathbf{q}}}_t\parallel }^2$, for a given positive scalar $\gamma$ [13]. Thus a non-negative scalar-valued function, $V({\mathbf{q}}^a,t)$, is sought that satisfies the fundamental property [14]:

$$\begin{array}{cc}\hfill & E\left\{V({\mathbf{q}}_t^a,t)\right\}\le E\{V({\mathbf{q}}_s^a,s)+{\int }_s^t\{{\gamma }^2\parallel \mathbf{v}\left(\tau \right){\parallel }^2-\parallel {\tilde{\mathbf{q}}}_{\tau }{\parallel }^{2}d\tau \}\forall 0\le s\le t\le T\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & E\left\{V({\mathbf{q}}_0^a,0)\right\}\le {\gamma }^{2}E\{\parallel {\tilde{\mathbf{q}}}_0{\parallel }^2\}\hfill \end{array}$$

(20)

for all ${\mathbf{q}}^a$ and for all admissible $\mathbf{v}\left(t\right)$. When Eq. (19) is satisfied, the function V is called a storage function for the supply rate S. A sufficient condition for Eq. (19) is:

$$\begin{array}{c}\hfill E\left\{dV({\mathbf{q}}^a,t)\right\}\le E\{{\gamma }^2\parallel \mathbf{v}\left(t\right){\parallel }^2-\parallel {\tilde{\mathbf{q}}}_t{\parallel }^2\}\forall 0\le t\le T\end{array}$$

(21)

for all ${\mathbf{q}}^a$ and for all admissible $\mathbf{v}\left(t\right)$, where $dV$ is the Itô differential of the function V. Notice that the processes $\left\{{\mathsf{\Omega }}_t\right\}$ and $\left\{H_t\right\}$ are observed and thus belong to the information pattern. Let ${\mathcal{I}}^t$ denote the information pattern at time t, which includes the observations of the angular velocity and of the LOS until t, i.e. ${\mathcal{I}}^t=\{{\left\{\omega \right\}}_0^t,{\left\{\mathbf{b}\right\}}_0^t\}$. In the development of the sufficiency conditions, one can use the following property of the expectation operator to write (21) as follows:

$$\begin{array}{cc}\hfill & E\left\{E\right\{dV({\mathbf{q}}^a,t)-{\gamma }^2\parallel {\mathbf{v}}_t{\parallel }^2+\parallel {\tilde{\mathbf{q}}}_t{\parallel }^2|{\mathcal{I}}^t\left\}\right\}\le 0\hfill \end{array}$$

(22)

for all ${\mathbf{q}}^a$ and t. Then a sufficient condition for (22) is

$$\begin{array}{cc}\hfill & E\{dV({\mathbf{q}}^a,t)-{\gamma }^2\parallel {\mathbf{v}}_t{\parallel }^2+\parallel {\tilde{\mathbf{q}}}_t{\parallel }^2|{\mathcal{I}}^t\}\le 0\hfill \end{array}$$

(23)

for all ${\mathbf{q}}^a$ and t. Given the conditioning on ${\mathcal{I}}^t$ in (23), the functions ${\mathsf{\Omega }}_t$ and $H_t$ can be considered deterministic functions. Using the Itô differentiation rule [17, p. 112] and dropping the expectation operator on both sides of Eq. (21) yields the following sufficient condition, i.e. the Hamilton-Jacobi-Bellman (HJB) inequality for $V({\mathbf{q}}^a,t)$:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial V}{\partial t}}+{\displaystyle \frac{\partial V}{\partial {\mathbf{q}}^a}}{F}^a{\mathbf{q}}^a+{\displaystyle \frac{1}{2}}{\sigma }_b^2\mathrm{tr}\left\{G{G}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}\right\}+{\displaystyle \frac{1}{2}}{\sigma }_{ϵ}^2\sum _{i=1}^3{{\mathbf{g}}_2^i}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}{\mathbf{g}}_2^i\le {\gamma }^2{\parallel \mathbf{v}\parallel }^2-{\parallel \tilde{\mathbf{q}}\parallel }^2\end{array}$$

(24)

for all $0\le t\le T$, ${\mathbf{q}}^a$, and for all admissible $\mathbf{v}\left(t\right)$, where G and ${\mathbf{g}}_2^i$ are functions of ${\mathbf{q}}^a$. Dropping the integral and the expectation operators is a standard procedure in stochastic $H_{\infty }$ estimation and control [13] or in stochastic optimal control [18, p. 321].

2.3.1. Case where $\mathbf{v}=\{{\sigma }_{ϵ},{\sigma }_b\}$

Assume both ${\sigma }_{ϵ}$ and ${\sigma }_b$ are perturbations. Bringing all terms to the left-hand-side (LHS) of Eq. (24) yields

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial V}{\partial t}}+{\displaystyle \frac{\partial V}{\partial {\mathbf{q}}^a}}{F}^a{\mathbf{q}}^a+{{\mathbf{q}}^a}^{T}L{\mathbf{q}}^a+\left[{\displaystyle \frac{1}{2}}\mathrm{tr}\left\{G{G}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}\right\}-{\gamma }^2\right]{\sigma }_b^2+\left[{\displaystyle \frac{1}{2}}\sum _{i=1}^3{{\mathbf{g}}_2^i}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}{\mathbf{g}}_2^i-{\gamma }^2\right]{\sigma }_{ϵ}^2\le 0\end{array}$$

(25)

where L denotes the following $8\times 8$ matrix

$$\begin{array}{c}\hfill L=\left[\begin{array}{cc}{O}_4& O_4\\ O_4& I_4\end{array}\right]\end{array}$$

(26)

and $I_4$ denotes the four-dimensional identity matrix. For a solution to exist for all ${\mathbf{q}}^a$ and for all admissible $\{{\sigma }_{ϵ},{\sigma }_b\}$ the coefficients multiplying the arbitrary disturbance functions, ${\sigma }_{ϵ}$ and ${\sigma }_b$, in Eq. (25) must be negative, which yields the following conditions:

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2}}\mathrm{tr}\left\{G{G}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}\right\}-{\gamma }^2\le 0\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2}}\sum _{i=1}^3{{\mathbf{g}}_2^i}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}{\mathbf{g}}_2^i-{\gamma }^2\le 0\hfill \end{array}$$

(28)

for all $0\le t\le T$, ${\mathbf{q}}^a$. Thus, according to Eq. (25), any non-zero disturbance will only add a negative term to the LHS, increasing the system’s dissipativity for the chosen supply rate. Henceforth, the worst-case scenario consists in vanishing disturbances e.g.,

$$\begin{array}{cc}\hfill & {\sigma }_{ϵ}^{*}=0\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & {\sigma }_b^{*}=0\hfill \end{array}$$

(30)

Some interpretation of Eqs. (29),(30) is required. It might appear at first sight counter-intuitive that the worst-case scenario for an estimator consists of vanishing noise variances. However, the proposed solution follows the attenuation $H_{\infty }$ criterion (12), not a minimum error condition. Henceforth performance is measured via the ratio between the (${\mathcal{L}}_2$ norms of the) estimation error and the noise variance, not via the estimation error itself. It is thus not surprising that the proposed estimator will perform better, i.e. attenuate better in the presence of high variances than in the presence of low variances.

Summary:

Using Eqs. (29) and (30) in Eq. (25) yields the following sufficient conditions. For all ${\mathbf{q}}^a$ and $0\le t\le T$,

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial V}{\partial t}}+{\displaystyle \frac{\partial V}{\partial {\mathbf{q}}^a}}{F}^a{\mathbf{q}}^a+{{\mathbf{q}}^a}^{T}L{\mathbf{q}}^a\le 0\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2}}\mathrm{tr}\left\{G{G}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}\right\}-{\gamma }^2\le 0\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2}}\sum _{i=1}^3{{\mathbf{g}}_2^i}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}{\mathbf{g}}_2^i-{\gamma }^2\le 0\hfill \end{array}$$

(33)

2.3.2. Case where $\mathbf{v}={\sigma }_{ϵ}$

Consider the case where ${\sigma }_{ϵ}$ is the sole perturbation. In this case, the HJB inequality evaluated at the worst-case yields the following sufficient conditions: for all ${\mathbf{q}}^a$ and $0\le t\le T$,

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial V}{\partial t}}+{\displaystyle \frac{\partial V}{\partial {\mathbf{q}}^a}}F{\mathbf{q}}^a+{{\mathbf{q}}^a}^{T}L{\mathbf{q}}^a+{\displaystyle \frac{1}{2}}{\sigma }_b^2\mathrm{tr}\left\{G{G}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}\right\}\le 0\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2}}\sum _{i=1}^3{{\mathbf{g}}_2^i}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}{\mathbf{g}}_2^i-{\gamma }^2\le 0\hfill \end{array}$$

(35)

In Eq. (34), the parameter ${\sigma }_b$ is a known deterministic parameter.

2.4. Storage function V

A standard approach to circumvent the formidable task of solving the partial differential inequality for V [Eq. (24)] consists of guessing the solution in a parameterized form and developing sufficient conditions for its parameters. The classical quadratic form for V will be used here:

$$\begin{array}{c}\hfill V({\mathbf{q}}^a,t)={{\mathbf{q}}^a}^T{P}_t{\mathbf{q}}^a\end{array}$$

(36)

where $P_t$ is assumed to be symmetric, positive, and block diagonal, i.e.,

$$\begin{array}{c}\hfill P_t=\left[\begin{array}{cc}{\widehat{P}}_t& O_4\\ O_4& {\tilde{P}}_t\end{array}\right]\end{array}$$

(37)

2.5. Sufficient conditions on the matrices $\widehat{K}$, $\tilde{P}$, $\widehat{P}$

Using Eqs. (36),(37) in Eq. (31), straightforward computations yield the following identity:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial V}{\partial t}}+{\displaystyle \frac{\partial V}{\partial {\mathbf{q}}^a}}{F}^a{\mathbf{q}}^a+{{\mathbf{q}}^a}^{T}L{\mathbf{q}}^a=\left[\begin{array}{cc}{{\widehat{\mathbf{q}}}_t}^T& {{\tilde{\mathbf{q}}}_t}^T\end{array}\right]\left[\begin{array}{cc}{\displaystyle \frac{d{\widehat{P}}_t}{dt}}+{{\widehat{F}}_t}^T{\widehat{P}}_t+{\widehat{P}}_t{\widehat{F}}_t& O_4\\ O_4& {\displaystyle \frac{d{\tilde{P}}_t}{dt}}+{{\widehat{F}}_t}^T{\tilde{P}}_t+{\tilde{P}}_t{\widehat{F}}_t+I_4\end{array}\right]\left[\begin{array}{c}{\widehat{\mathbf{q}}}_t\\ {\tilde{\mathbf{q}}}_t\end{array}\right]\end{array}$$

(38)

where ${\widehat{F}}_t$ is defined as follows:

$$\begin{array}{c}\hfill {\widehat{F}}_t={\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t-{\widehat{K}}_t{H}_t\end{array}$$

(39)

Exploiting the following property of the matrix $\mathsf{\Xi }$ [5, Eq. (8.20b)],

$$\begin{array}{cc}\hfill & \mathsf{\Xi }{\mathsf{\Xi }}^T={\mathbf{q}}_t^T{\mathbf{q}}_t{I}_4-{\mathbf{q}}_t{\mathbf{q}}_t^T\hfill \end{array}$$

(40)

yields the following identities

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2}}\sum _{i=1}^3{{\mathbf{g}}_2^i}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}{\mathbf{g}}_2^i={\mathbf{q}}_t^T\left[{\displaystyle \frac{1}{4}}(\mathrm{tr}{\tilde{P}}_t{I}_4-{\tilde{P}}_t)\right]{\mathbf{q}}_t\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{2}}\mathrm{tr}\left\{G{G}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}\right\}={\mathbf{q}}_t^T\left\{{\displaystyle \frac{1}{4}}\left[\mathrm{tr}\left({\widehat{K}}_t^T{\tilde{P}}_t{\widehat{K}}_t\right)I_4-{\widehat{K}}_t^T{\tilde{P}}_t{\widehat{K}}_t\right]\right\}{\mathbf{q}}_t\hfill \end{array}$$

(42)

2.5.1. Convexity condition with respect to ${\sigma }_{ϵ}$

Using Eq. (41) in Eq. (33) [Eq. (35)], and noting that the inequality must be verified for all ${\mathbf{q}}_t$ yields the following condition on ${\tilde{P}}_t$:

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{4}}(\mathrm{tr}{\tilde{P}}_t{I}_4-{\tilde{P}}_t)-{\gamma }^2{I}_4\le 0\end{array}$$

(43)

This condition must be satisfied in both cases, whether ${\sigma }_b$ is known or is an unknown perturbation. The latter condition can be easily reformulated in terms of the characteristic values of the symmetric positive matrix $\tilde{P}$. Let ${\tilde{\lambda }}_i$, $i=1,2,3,4$, denote the four positive eigenvalues of $\tilde{P}$, then Eq. (43) is equivalent to the following condition:

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{4}}max\left\{{\tilde{\lambda }}_2+{\tilde{\lambda }}_3+{\tilde{\lambda }}_4,{\tilde{\lambda }}_1+{\tilde{\lambda }}_3+{\tilde{\lambda }}_4,{\tilde{\lambda }}_1+{\tilde{\lambda }}_2+{\tilde{\lambda }}_4,{\tilde{\lambda }}_1+{\tilde{\lambda }}_2+{\tilde{\lambda }}_3\right\}\le {\gamma }^2\hfill \end{array}$$

(44)

2.5.2. Case where $\mathbf{v}=\{{\sigma }_{ϵ},{\sigma }_b\}$

Using Eq. (42) in Eq. (32), and noting that the inequality must be verified for all ${\mathbf{q}}_t$, yields the following condition on ${\tilde{P}}_t$ and $\widehat{K}$:

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{4}}\left[\mathrm{tr}\left({\widehat{K}}_t^T{\tilde{P}}_t{\widehat{K}}_t\right)I_4-{\widehat{K}}_t^T{\tilde{P}}_t{\widehat{K}}_t\right]-{\gamma }^2{I}_4\le 0\hfill \end{array}$$

(45)

Using Eq. (38) in Eq. (31) yields two uncoupled differential matrix inequalities for ${\widehat{P}}_t$ and ${\tilde{P}}_t$, respectively:

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d{\widehat{P}}_t}{dt}}+F_t^T{\widehat{P}}_t+{\widehat{P}}_t{F}_t\le 0\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d{\tilde{P}}_t}{dt}}+F_t^T{\tilde{P}}_t+{\tilde{P}}_t{F}_t+I_4\le 0\hfill \end{array}$$

(47)

Assuming that the matrices ${\widehat{P}}_t$ and ${\tilde{P}}_t$ are identical for the sake of simplicity allows dropping Eq. (46). Combining Eqs. (43), (45), and (47) yields the inequalities to be solved for $\widehat{K}$ and $\tilde{P}$.

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d{\tilde{P}}_t}{dt}}+F_t^T{\tilde{P}}_t+{\tilde{P}}_t{F}_t+I_4\le 0\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{4}}(\mathrm{tr}{\tilde{P}}_t{I}_4-{\tilde{P}}_t)-{\gamma }^2{I}_4\le 0\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{4}}(\mathrm{tr}{M}_t{I}_4-M_t)-{\gamma }^2{I}_4\le 0\hfill \end{array}$$

(50)

where

$$\begin{array}{cc}\hfill & F_t={\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t-\widehat{K}{H}_t\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill & M_t={\widehat{K}}^T{\tilde{P}}_t\widehat{K}\hfill \end{array}$$

(52)

Notice that the left-hand sides (LHS) of the above inequalities are independent of the quaternion estimate; thus, the filter gain $\widehat{K}$ is also independent of ${\widehat{\mathbf{q}}}_t$. It will be denoted by $K_t$ in the following.

Sufficient conditions in the form of Linear Matrix Inequalities:

Since the above inequalities are not linear with respect to $\tilde{P}$ and K, some manipulations are required in order to bring them to a Linear Matrix Inequality (LMI) structure. The bilinear dependence with respect to $\tilde{P}$ and K is readily coped with via a standard parametrization approach. Let ${\tilde{Y}}_t$ denote the following four-dimensional matrix:

$$\begin{array}{c}\hfill {\tilde{Y}}_t={\tilde{P}}_t{K}_t\end{array}$$

(53)

then, using Eq. (53) in Eq. (48) yields

$$\begin{array}{c}\hfill {\displaystyle \frac{d{\tilde{P}}_t}{dt}}+{\displaystyle \frac{1}{2}}({\mathsf{\Omega }}_t^T{\tilde{P}}_t+{\tilde{P}}_t{\mathsf{\Omega }}_t)-(H_t^T{\tilde{Y}}_t^T+{\tilde{Y}}_t{H}_t)+I_4\le 0\end{array}$$

(54)

To circumvent the difficulty arising from the quadratic structure of $M_t$ with respect to ${\tilde{P}}_t$ and K, a symmetric positive definite matrix $W_t$ is sought such that

$$\begin{array}{c}\hfill M_t-W_t={\tilde{Y}}_t^T{\tilde{P}}_t^{-1}{\tilde{Y}}_t-W_t\le 0\end{array}$$

(55)

Notice that ${\tilde{P}}_t^{-1}$ exists since ${\tilde{P}}_t$ is assumed to be positive definite. Then, the following bounds on the LHS of Eq. (50) are used:

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{4}}(\mathrm{tr}{M}_t{I}_4-M_t)-{\gamma }^2{I}_4\le ({\displaystyle \frac{1}{4}}\mathrm{tr}{M}_t-{\gamma }^2)I_4\le ({\displaystyle \frac{1}{4}}\mathrm{tr}{W}_t-{\gamma }^2)I_4\le 0\end{array}$$

(56)

and Eq. (50) is replaced with the following sufficient condition on W:

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{4}}\mathrm{tr}{W}_t-{\gamma }^2\le 0\hfill \end{array}$$

(57)

where W, $\tilde{Y}$, and $\tilde{P}$ satisfy Eq. (55), which by the Schur complement can be written as the following LMI :

$$\begin{array}{cc}\hfill & \left[\begin{array}{cc}-W_t& -{\tilde{Y}}_t\\ -{\tilde{Y}}_t^T& -{\tilde{P}}_t\end{array}\right]\le 0\hfill \end{array}$$

(58)

Notice that the successive bounds in Eq. (56) yield a sufficient condition for the attenuation filtering problem.

2.5.3. Case where $\mathbf{v}={\sigma }_{ϵ}$

Using Eq. (42) and the definition of ${\mathbf{q}}^a$ yields the following identity:

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{2}}{\sigma }_b^2\mathrm{tr}\left\{G{G}^T{\displaystyle \frac{{\partial }^{2}V}{\partial {{\mathbf{q}}^a}^2}}\right\}& ={\mathbf{q}}_t^T\left\{{\displaystyle \frac{{\sigma }_b^2}{4}}\left[\begin{array}{c}\mathrm{tr}\left({\widehat{K}}_t^T{\tilde{P}}_t{\widehat{K}}_t\right)I_4-{\widehat{K}}_t^T{\tilde{P}}_t{\widehat{K}}_t\end{array}\right]\right\}{\mathbf{q}}_t\hfill \\ \hfill & =\left[\begin{array}{cc}{{\widehat{\mathbf{q}}}_t}^T& {{\tilde{\mathbf{q}}}_t}^T\end{array}\right]\left[\begin{array}{c}{I}_4\\ I_4\end{array}\right]\left\{{\displaystyle \frac{{\sigma }_b^2}{4}}\left[\begin{array}{c}\mathrm{tr}\left({\widehat{K}}_t^T{\tilde{P}}_t{\widehat{K}}_t\right)I_4-{\widehat{K}}_t^T{\tilde{P}}_t{\widehat{K}}_t\end{array}\right]\right\}\left[\begin{array}{cc}{I}_4& I_4\end{array}\right]\left[\begin{array}{c}{\widehat{\mathbf{q}}}_t\\ {\tilde{\mathbf{q}}}_t\end{array}\right]\hfill \end{array}$$

(59)

Using Eqs. (38) and (59) in Eq. (34) yields

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{{\widehat{\mathbf{q}}}_t}^T& {{\tilde{\mathbf{q}}}_t}^T\end{array}\right]\left[\begin{array}{cc}{\displaystyle \frac{d{\widehat{P}}_t}{dt}}+{{\widehat{F}}_t}^T{\widehat{P}}_t+{\widehat{P}}_t{\widehat{F}}_t+{\displaystyle \frac{{\sigma }_b^2}{4}}\left[\begin{array}{c}\mathrm{tr}{\widehat{M}}_t{I}_4-{\widehat{M}}_t\end{array}\right]& {\displaystyle \frac{{\sigma }_b^2}{4}}\left[\begin{array}{c}\mathrm{tr}{\widehat{M}}_t{I}_4-{\widehat{M}}_t\end{array}\right]\\ {\displaystyle \frac{{\sigma }_b^2}{4}}\left[\begin{array}{c}\mathrm{tr}{\widehat{M}}_t{I}_4-{\widehat{M}}_t\end{array}\right]& {\displaystyle \frac{d{\tilde{P}}_t}{dt}}+{{\widehat{F}}_t}^T{\tilde{P}}_t+{\tilde{P}}_t{\widehat{F}}_t+{\displaystyle \frac{{\sigma }_b^2}{4}}\left[\begin{array}{c}\mathrm{tr}{\widehat{M}}_t{I}_4-{\widehat{M}}_t\end{array}\right]+I_4\end{array}\right]\left[\begin{array}{c}{\widehat{\mathbf{q}}}_t\\ {\tilde{\mathbf{q}}}_t\end{array}\right]\le 0\end{array}$$

(60)

for all $({\widehat{\mathbf{q}}}_t,{\tilde{\mathbf{q}}}_t,t)$, where ${\widehat{M}}_t={\widehat{K}}_t^T{\tilde{P}}_t{\widehat{K}}_t$. Assuming that $\widehat{P}=\tilde{P}$ as in the previous case yields the following matrix differential inequality

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\displaystyle \frac{d{\tilde{P}}_t}{dt}}+F^T{\tilde{P}}_t+{\tilde{P}}_{t}F+{\displaystyle \frac{{\sigma }_b^2}{4}}\left[\begin{array}{c}\left(\mathrm{tr}M\right)I_4-M\end{array}\right]& {\displaystyle \frac{{\sigma }_b^2}{4}}\left[\begin{array}{c}\left(\mathrm{tr}M\right)I_4-M\end{array}\right]\\ {\displaystyle \frac{{\sigma }_b^2}{4}}\left[\begin{array}{c}\left(\mathrm{tr}M\right)I_4-M\end{array}\right]& {\displaystyle \frac{d{\tilde{P}}_t}{dt}}+F^T{\tilde{P}}_t+{\tilde{P}}_{t}F+I_4+{\displaystyle \frac{{\sigma }_b^2}{4}}\left[\begin{array}{c}\left(\mathrm{tr}M\right)I_4-M\end{array}\right]\end{array}\right]\le 0\end{array}$$

(61)

for all $0\le t\le T$, where

$$\begin{array}{cc}\hfill & F={\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t-K_t{H}_t\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill & M=K_t^T{\tilde{P}}_t{K}_t\hfill \end{array}$$

(63)

Thus, when ${\sigma }_b$ is a known parameter, Eqs. (61)-(63), and Eq. (43) are the sufficient conditions for $\tilde{P}$ and K.

Sufficient conditions in the form of LMI:

Introducing a matrix variable $W_t$ that satisfies Eqs. (53), (55), and using the same upper bounds sequence as in Eq. (56) yields the following differential LMI for this case:

$$\begin{array}{cc}\hfill & \left[\begin{array}{cc}{\displaystyle \frac{d{\tilde{P}}_t}{dt}}+{\displaystyle \frac{1}{2}}({\mathsf{\Omega }}_t^T{\tilde{P}}_t+{\tilde{P}}_t{\mathsf{\Omega }}_t)-(H_t^T{\tilde{Y}}_t^T+{\tilde{Y}}_t{H}_t)+{\displaystyle \frac{{\sigma }_b^2}{4}}\left(\mathrm{tr}{W}_t\right)I_4& {\displaystyle \frac{{\sigma }_b^2}{4}}\left(\mathrm{tr}{W}_t\right)I_4\\ {\displaystyle \frac{{\sigma }_b^2}{4}}\left(\mathrm{tr}{W}_t\right)I_4& {\displaystyle \frac{d{\tilde{P}}_t}{dt}}+{\displaystyle \frac{1}{2}}({\mathsf{\Omega }}_t^T{\tilde{P}}_t+{\tilde{P}}_t{\mathsf{\Omega }}_t)-(H_t^T{\tilde{Y}}_t^T+{\tilde{Y}}_t{H}_t)+I_4+{\displaystyle \frac{{\sigma }_b^2}{4}}\left(\mathrm{tr}{W}_t\right)I_4\end{array}\right]\le 0\hfill \end{array}$$

(64)

2.6. Quaternion stochastic $H_{\infty }$ filters summary

Given ${\widehat{\mathbf{q}}}_0$, choose $\tilde{P}\left(0\right)$ such that Eq. (20) is satisfied. Solve the following set of (differential) LMIs for ${\tilde{P}}_t={\tilde{P}}_t^T>0\in {\mathbb{R}}^4$, ${\tilde{Y}}_t\in {\mathbb{R}}^4$, and $W_t=W_t^T>0\in {\mathbb{R}}^4$:

2.6.1. Case where $\mathbf{v}=\{{\sigma }_{ϵ},{\sigma }_b\}$

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d{\tilde{P}}_t}{dt}}+{\displaystyle \frac{1}{2}}({\mathsf{\Omega }}_t^T{\tilde{P}}_t+{\tilde{P}}_t{\mathsf{\Omega }}_t)-(H_t^T{\tilde{Y}}_t^T+{\tilde{Y}}_t{H}_t)+I_4\le 0\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill & \left[\begin{array}{cc}-W_t& -{\tilde{Y}}_t\\ -{\tilde{Y}}_t^T& -{\tilde{P}}_t\end{array}\right]\le 0\hfill \end{array}$$

(66)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{4}}(\mathrm{tr}{\tilde{P}}_t{I}_4-{\tilde{P}}_t)-{\gamma }^2{I}_4\le 0\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{4}}\mathrm{tr}{W}_t-{\gamma }^2\le 0\hfill \end{array}$$

(68)

2.6.2. Case where $\mathbf{v}={\sigma }_{ϵ}$

$$\begin{array}{cc}\hfill & \left[\begin{array}{cc}{\displaystyle \frac{d{\tilde{P}}_t}{dt}}+{\displaystyle \frac{1}{2}}({\mathsf{\Omega }}_t^T{\tilde{P}}_t+{\tilde{P}}_t{\mathsf{\Omega }}_t)-(H_t^T{\tilde{Y}}_t^T+{\tilde{Y}}_t{H}_t)+{\displaystyle \frac{{\sigma }_b^2}{4}}\left(\mathrm{tr}{W}_t\right)I_4& {\displaystyle \frac{{\sigma }_b^2}{4}}\left(\mathrm{tr}{W}_t\right)I_4\\ {\displaystyle \frac{{\sigma }_b^2}{4}}\left(\mathrm{tr}{W}_t\right)I_4& {\displaystyle \frac{d{\tilde{P}}_t}{dt}}+{\displaystyle \frac{1}{2}}({\mathsf{\Omega }}_t^T{\tilde{P}}_t+{\tilde{P}}_t{\mathsf{\Omega }}_t)-(H_t^T{\tilde{Y}}_t^T+{\tilde{Y}}_t{H}_t)+I_4+{\displaystyle \frac{{\sigma }_b^2}{4}}\left(\mathrm{tr}{W}_t\right)I_4\end{array}\right]\le 0\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill & \left[\begin{array}{cc}-W_t& -{\tilde{Y}}_t\\ -{\tilde{Y}}_t^T& -{\tilde{P}}_t\end{array}\right]\le 0\hfill \end{array}$$

(70)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{4}}(\mathrm{tr}{\tilde{P}}_t{I}_4-{\tilde{P}}_t)-{\gamma }^2{I}_4\le 0\hfill \end{array}$$

(71)

For any pair of matrices $({\tilde{Y}}_t,{\tilde{P}}_t)$, compute the gain $K_t$ using

$$\begin{array}{cc}\hfill & K_t={\tilde{P}}_t^{-1}{\tilde{Y}}_t\hfill \end{array}$$

(72)

and compute the estimated quaternion via the estimator differential equation

$$\begin{array}{cc}\hfill & {\dot{\widehat{\mathbf{q}}}}_t=\left[{\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t-K_t{H}_t\right]{\widehat{\mathbf{q}}}_t\hfill \end{array}$$

(73)

Remark 1: The estimator equation, (73), is not designed to preserve the quaternion unit-norm property. For that purpose, a normalization stage of the estimate is performed along the estimation process [4,15]

$$\begin{array}{cc}\hfill & \widehat{\mathbf{q}}={\displaystyle \frac{\widehat{\mathbf{q}}}{\parallel \widehat{\mathbf{q}}\parallel }}\hfill \end{array}$$

(74)

Remark 2: A key feature of the above filters lies in the fact that the gain computations are independent of the estimated process. As a result, the gain values are insensitive to the initial estimation errors, which are often causes of divergence in linearization-based filtering techniques, like the extended Kalman filter. An additional essential by-product is that the estimate differential equation Eq. (73) can be integrated as an ordinary differential equation.

Remark 3: The above algorithms are solved using standard primal-dual interior-point methods, as implemented in SeDuMi [19,20]. The method formulates a minimization problem over $\gamma$ subject to the constraints described in Eqs. (65)-(68). For the solver SeDuMi, an assessment of the computational complexity is $\mathcal{O}\left(n^4\right)$ for the $2n^2+n$ decision variables, where $n=4$. Compared to the standard computational complexity of a Kalman filter, i.e. $\mathcal{O}\left(n^3\right)$, this yields a ratio of 4.

Remark 4: Inspired by [21], discrete approximations of the differential LMIs are developed via finite-difference formulas. For example, Eq. (65) is implemented as follows:

$$\begin{array}{c}\hfill {\displaystyle \frac{{\tilde{P}}_{k+1}-{\tilde{P}}_k}{\Delta t}}+{\displaystyle \frac{1}{2}}({\mathsf{\Omega }}_{k+1}^T{\tilde{P}}_{k+1}+{\tilde{P}}_{k+1}{\mathsf{\Omega }}_{k+1})-(H_{k+1}^T{\tilde{Y}}_{k+1}^T+{\tilde{Y}}_{k+1}{H}_{k+1})+I_4\le 0\end{array}$$

(75)

where $\Delta t$ denotes the time increment, and $k=0,1,...,N=T/\Delta t$.

3. Quaternion and Gyro Drift Estimation

3.1. Statement of the problem

Assuming that the rate gyro error consists of white noise and a bias, we consider the following stochastic dynamical system in Itô form:

$$\begin{array}{cc}\hfill & \mathbf{d}{\mathbf{q}}_{\mathbf{t}}={\displaystyle \frac{1}{2}}\mathsf{\Omega }({\omega }_t-{\mathbf{c}}_t){\mathbf{q}}_{t}dt-{\displaystyle \frac{\mathtt{1}}{\mathtt{2}}}\mathsf{\Xi }\left({\mathbf{q}}_t\right){\sigma }_{ϵ}\left(t\right){\mathbf{d}\beta }_t;\mathbf{q}\left(0\right)\stackrel{a.e.}{=}{\mathbf{q}}_0;t\in [0,T]\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill & \mathbf{d}{\mathbf{c}}_{\mathbf{t}}={\sigma }_c\left(t\right)\mathbf{d}{\nu }_{\mathbf{t}};\mathbf{c}\left(0\right)\stackrel{a.e.}{=}{\mathbf{c}}_0\hfill \end{array}$$

(77)

$$\begin{array}{cc}\hfill & \mathbf{d}{\mathbf{y}}_{\mathbf{t}}=H_t{\mathbf{q}}_{t}dt-{\displaystyle \frac{1}{2}}\mathsf{\Xi }\left({\mathbf{q}}_t\right){\sigma }_b\left(t\right){\mathbf{d}\mathit{\eta }}_t\hfill \end{array}$$

(78)

where ${\mathbf{c}}_t$ denotes the additive drift, modeled as a random walk process with mean ${\mathbf{c}}_0$ and variance parameter ${\sigma }_c\left(t\right)$. In Eq. (77), ${\nu }_t$ denotes a standard Brownian motion that is independent of ${\beta }_i$ and ${\mathit{\eta }}_t$. Equations (76), (77) stem from a straightforward extension of the quaternion SDE of Section 2.

The filtering problem consists of estimating the quaternion ${\mathbf{q}}_t$ and the gyro drift ${\mathbf{c}}_t$ from the LOS measurements in the presence of unknown noise standard deviations ${\sigma }_{ϵ}\left(t\right)$, ${\sigma }_b\left(t\right)$, and ${\sigma }_c\left(t\right)$. Assuming that ${\sigma }_{ϵ}\left(t\right)$, ${\sigma }_b\left(t\right)$, and ${\sigma }_c\left(t\right)$ are stochastic non-anticipative processes with finite second-order moments, we consider the following estimator:

$$\begin{array}{cc}\hfill & \mathbf{d}{\widehat{\mathbf{q}}}_{\mathbf{t}}={\displaystyle \frac{1}{2}}\mathsf{\Omega }({\omega }_t-{\widehat{\mathbf{c}}}_t){\widehat{\mathbf{q}}}_{t}dt+K_q\left(\mathbf{d}{\mathbf{y}}_{\mathbf{t}}-{\widehat{\mathbf{q}}}_{t}dt\right)\hfill \end{array}$$

(79)

$$\begin{array}{cc}\hfill & \mathbf{d}{\widehat{\mathbf{c}}}_{\mathbf{t}}=K_c\left(\mathbf{d}{\mathbf{y}}_{\mathbf{t}}-{\widehat{\mathbf{q}}}_{t}dt\right)\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill & \widehat{\mathbf{q}}\left(0\right)={\widehat{\mathbf{q}}}_0,\widehat{\mathbf{c}}\left(0\right)={\widehat{\mathbf{c}}}_0\hfill \end{array}$$

(81)

Let ${\tilde{\mathbf{q}}}_t$ and ${\tilde{\mathbf{c}}}_t$ denote the additive quaternion and biases estimation error, i.e.,

$$\begin{array}{cc}\hfill & {\tilde{\mathbf{q}}}_t={\mathbf{q}}_t-{\widehat{\mathbf{q}}}_t\hfill \end{array}$$

(82)

$$\begin{array}{cc}\hfill & {\tilde{\mathbf{c}}}_t={\mathbf{c}}_t-{\widehat{\mathbf{c}}}_t\hfill \end{array}$$

(83)

Given a scalar $\gamma >0$, we seek the gains $K_q,K_c$ such that the following $H_{\infty }$ criterion is satisfied:

$$\begin{array}{c}\hfill E\{{\int }_0^T(\parallel {\tilde{\mathbf{q}}}_t{\parallel }^2+\parallel {\tilde{\mathbf{c}}}_t{\parallel }^2\left)dt\right\}\le {\gamma }^{2}E\{\parallel {\tilde{\mathbf{q}}}_0{\parallel }^2+\parallel {\tilde{\mathbf{c}}}_0{\parallel }^2+{\int }_0^T\parallel {\mathbf{v}}_t{\parallel }^{2}dt\}\end{array}$$

(84)

under the constraints (76)-(78), and where ${\mathbf{v}}_t$ denotes the augmented process of admissible disturbance functions, i.e., ${\mathbf{v}}_t=\{{\sigma }_{ϵ}\left(t\right),{\sigma }_b\left(t\right),{\sigma }_c\left(t\right)\}$. Whenever Eq. (84) is true, it is said that the $L_2$-gain property is satisfied from $\{{\tilde{\mathbf{q}}}_0,{\tilde{\mathbf{c}}}_0,{\mathbf{v}}_t\}$ to $\{{\tilde{\mathbf{q}}}_t,{\tilde{\mathbf{c}}}_t\}$, for $0\le t\le T$.

3.2. Design model development

The SDE of the quaternion-drift system is compactly rewritten as follows:

$$\begin{array}{cc}\hfill & \left[\begin{array}{c}\mathbf{d}{\mathbf{q}}_{\mathbf{t}}\\ \mathbf{d}{\mathbf{c}}_{\mathbf{t}}\end{array}\right]=\left[\begin{array}{c}{\displaystyle \frac{1}{2}}\mathsf{\Omega }({\omega }_t-{\mathbf{c}}_t){\mathbf{q}}_t\\ O_{3\times 1}\end{array}\right]dt+\left[\begin{array}{cc}-{\displaystyle \frac{1}{2}}\mathsf{\Xi }\left({\mathbf{q}}_t\right)& O_{4\times 3}\\ O_3& I_3\end{array}\right]\left[\begin{array}{cc}{\sigma }_{ϵ}\left(t\right)I_3& O_3\\ O_3& {\sigma }_c\left(t\right)I_3\end{array}\right]\left[\begin{array}{c}{\mathbf{d}\beta }_t\\ \mathbf{d}{\nu }_{\mathbf{t}}\end{array}\right]\hfill \end{array}$$

(85)

and the estimator is rewritten as follows:

$$\begin{array}{cc}\hfill & \left[\begin{array}{c}\mathbf{d}{\widehat{\mathbf{q}}}_{\mathbf{t}}\\ \mathbf{d}{\widehat{\mathbf{c}}}_{\mathbf{t}}\end{array}\right]=\left[\begin{array}{c}[{\displaystyle \frac{1}{2}}\mathsf{\Omega }({\omega }_t-{\widehat{\mathbf{c}}}_t)-K_q{H}_t]\\ -K_c{H}_t\end{array}\right]{\widehat{\mathbf{q}}}_{t}dt\hfill \\ \hfill & {\widehat{\mathbf{q}}}_t\left(0\right)={\widehat{\mathbf{q}}}_0;{\widehat{\mathbf{c}}}_t\left(0\right)={\widehat{\mathbf{c}}}_0\hfill \end{array}$$

(86)

The augmented process $\{{\widehat{\mathbf{q}}}_t,{\widehat{\mathbf{c}}}_t,\tilde{\mathbf{q}},\tilde{\mathbf{c}}\}$ is governed by the following SDE

$$\begin{array}{c}\hfill \left[\begin{array}{c}\mathbf{d}{\widehat{\mathbf{q}}}_{\mathbf{t}}\\ \mathbf{d}{\widehat{\mathbf{c}}}_{\mathbf{t}}\\ \mathbf{d}{\tilde{\mathbf{q}}}_{\mathbf{t}}\\ \mathbf{d}{\tilde{\mathbf{c}}}_{\mathbf{t}}\end{array}\right]=\left[\begin{array}{cccc}{\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t-K_q{H}_t& -{\displaystyle \frac{1}{2}}\mathsf{\Xi }\left({\widehat{\mathbf{q}}}_t\right)& O& O\\ -K_c{H}_t& O& O& O\\ O& O& {\displaystyle \frac{1}{2}}{\mathsf{\Omega }}_t-K_q{H}_t& -{\displaystyle \frac{1}{2}}\mathsf{\Xi }\left({\widehat{\mathbf{q}}}_t\right)\\ O& O& -K_c{H}_t& O\end{array}\right]\left[\begin{array}{c}{\widehat{\mathbf{q}}}_t\\ {\widehat{\mathbf{c}}}_t\\ {\tilde{\mathbf{q}}}_t\\ {\tilde{\mathbf{c}}}_t\end{array}\right]dt\\ \hfill +\left[\begin{array}{c}{O}_{4\times 3}\\ O_3\\ -{\displaystyle \frac{1}{2}}\mathsf{\Xi }\left({\widehat{\mathbf{q}}}_t\right)\\ O_3\end{array}\right]{\sigma }_{ϵ}\left(t\right){\mathbf{d}\beta }_t+\left[\begin{array}{c}{O}_{4\times 3}\\ O_3\\ O_{4\times 3}\\ I_3\end{array}\right]{\sigma }_c\left(t\right)\mathbf{d}{\nu }_{\mathbf{t}}+\left[\begin{array}{c}{O}_{4\times 3}\\ O_{4\times 3}\\ K_q{\displaystyle \frac{1}{2}}\mathsf{\Xi }\left({\widehat{\mathbf{q}}}_t\right)\\ K_c{\displaystyle \frac{1}{2}}\mathsf{\Xi }\left({\widehat{\mathbf{q}}}_t\right)\end{array}\right]{\sigma }_b\left(t\right){\mathbf{d}\mathit{\eta }}_t\end{array}$$

(87)

where second-order terms with respect to the noises ${\mathit{\beta }}_i$, ${\mathit{\nu }}_t$, ${\mathit{\eta }}_t$ and to the estimation errors $\tilde{\mathbf{q}}$, $\tilde{\mathbf{c}}$ have been neglected. Equation (87) may be re-written in the following compact form:

$$\begin{array}{cc}\hfill & \mathbf{d}{\mathbf{q}}_t^a=F^a{\mathbf{q}}_t^{a}dt+G_1\left({\mathbf{q}}_t^a\right){\sigma }_{ϵ}\left(t\right){\mathbf{d}\beta }_t+G_2\left({\mathbf{q}}_t^a\right){\sigma }_c\left(t\right)\mathbf{d}{\nu }_{\mathbf{t}}+G\left({\mathbf{q}}_t^a\right){\sigma }_b\left(t\right){\mathbf{d}\mathit{\eta }}_t\hfill \end{array}$$

(88)

The remainder of the filter development is straightforward and is omitted for the sake of brevity.

4. Numerical Simulation

This section is concerned with the numerical validation of the proposed approach in the drift-free case.

4.1. Description

Consider a spacecraft rotating around its center of mass with the following time-varying inertial angular velocity vector, ${\omega }^o\left(t\right)$:

$$\begin{array}{cc}\hfill & {\omega }^o\left(t\right)={[1,-1,1]}^{T}sin(2\pi t/150)[deg/sec]\hfill \end{array}$$

(89)

The measured angular velocity is computed according to

$$\begin{array}{cc}\hfill & \omega \left(t\right)={\omega }^o\left(t\right)+{\sigma }_{ϵ}ϵ\left(t\right)\hfill \end{array}$$

(90)

where $ϵ\left(t\right)$ is a standard zero-mean white Gaussian noise, e.g., $E\left\{ϵ\left(t\right)ϵ{\left(\tau \right)}^T\right\}=I_3\delta (t-\tau )$. Typical values of low-grade gyros are used, i.e., ${\sigma }_{ϵ}\in [{10}^{-4},{10}^0][rad/\sqrt{sec}]$. A time-varying line-of-sight measurement, ${\mathbf{b}}_t$, is assumed to be acquired. It is computed via the classical vector measurement model:

$$\begin{array}{cc}\hfill & {\mathbf{b}}_t=A\left[{\mathbf{q}}_t\right]{\mathbf{r}}_t+{\sigma }_b{\delta \mathbf{b}}_t\hfill \end{array}$$

(91)

where $\mathbf{r}\left(t\right)$ is randomly generated using a zero-mean standard multivariate normal distribution and the attitude matrix $A\left({\mathbf{q}}_t\right)$ is expressed as follows:

$$\begin{array}{cc}\hfill & A\left({\mathbf{q}}_t\right)=(q_t^2-{\mathbf{e}}_t^T{\mathbf{e}}_t)I_3+2{\mathbf{e}}_t{\mathbf{e}}_t^T-2q\left[{\mathbf{e}}_t\times \right]\hfill \end{array}$$

(92)

4.2. Attenuation of gyro ${\sigma }_{ϵ}$

In this section, the numerical study focuses on the impact of the gyro perturbation ${\sigma }_{ϵ}$ on the attenuation performance of the QHF. For that purpose, the parameter ${\sigma }_b$ is set to various known values while ${\sigma }_b$ is kept equal to ${10}^{-6}$ radian. Table 1 presents the values of Monte-Carlo (MC) averages over 500 runs lasting 500 seconds each of ${\delta }_q={max}_t∥{\mathbf{q}}_t-{\widehat{\mathbf{q}}}_t∥$ and of the ratio ${\displaystyle \frac{\delta q}{{\sigma }_{ϵ}\sqrt{\Delta t}}}$, where $\Delta t=0.1s$ is the gyro sampling time. The former is a measure of attitude estimation accuracy while the latter is a measure of attenuation performance. It can be seen that the QHF always converges, that the estimation accuracy is satisfying despite the extreme values of ${\sigma }_{ϵ}$, albeit degraded as ${\sigma }_{ϵ}$ increases, and that the attenuation performance improves with ${\sigma }_{ϵ}$.

Additional MC simulations (500 runs) were performed while varying the parameters ${\sigma }_{ϵ}$ and ${\sigma }_b$. Table 2 depicts the ratios of the time averages (over 6000 seconds) of the angular error, $\delta \varphi$, of the QHF over the QKF. The error $\delta \varphi$ is extracted from the rotational quaternion error’s fourth component. The magnitude of $\delta \varphi$ in the QHF appears in parenthesis (in degree) above the ratios. For a given ${\sigma }_{ϵ}$, the values of $\delta \varphi$ and the ratios increase with ${\sigma }_b$ because the attenuation quality is impaired. It turns out that the ratios are smaller than one in almost all test cases, i.e. the QHF produces a smaller bias than the QKF. The above results suggest that the QHF is advantageous when using low-grade gyros (high ${\sigma }_{ϵ}$) with fine LOS sensors (low ${\sigma }_b$).

4.3. Attenuation of ${\sigma }_{ϵ}$ and ${\sigma }_b$

Next, we test the performances of the QHF when both ${\sigma }_{ϵ}$ and ${\sigma }_b$ are unknown. For that purpose, we evaluate the actual attenuation ratio $AR\left(T\right)$ which is defined as follows:

$$\begin{array}{cc}\hfill & AR\left(T\right)={\displaystyle \frac{E\{{\int }_0^T\parallel \tilde{\mathbf{q}}\left(t\right){\parallel }^{2}dt\}}{E\{\parallel \tilde{\mathbf{q}}\left(0\right){\parallel }^2+{\int }_0^T({\sigma }_{ϵ}^2+{\sigma }_b^2)dt\}}}\hfill \end{array}$$

(93)

where the final time T is 500 sec, the integrals are numerically computed using a time step $\Delta t=0.1$ sec, and the expectations are computed as MC averages over 500 runs. Table 3 shows the values of $AR\left(500\right)$ for various ${\sigma }_{ϵ}$ and ${\sigma }_b$. It also features the steady-state MC means of the best-guaranteed level of attenuation, ${\gamma }_{QHF}^2$, which is calculated within the QHF.

In a nutshell, the performance index $AR\left(500\right)$ is not sensitive to variations in ${\sigma }_{ϵ},{\sigma }_b$ below a threshold of ${10}^{-2}$, above which it decreases rapidly, showing thus improved performance in terms of disturbance attenuation. A similar lack of sensitivity is observed for the parameter ${\gamma }_{QHF}^2$ over the full ranges of ${\sigma }_{ϵ}$ and ${\sigma }_b$, with a small but consistent improvement towards large values. Strikingly, the values of $AR\left(500\right)$ are significantly lower than those of ${\gamma }_{QHF}^2$. In more detail, the gap is about six-fold lower in the case of vanishing variances, and about 30 times lower for very large variances, when the pair $({\sigma }_{ϵ},{\sigma }_b)$ is equal to $(0.1,0.1)$. That is consistent with the $H_{\infty }$ filtering theory, i.e. vanishing disturbances are the worst case in terms of disturbance attenuation. The time variations of the MC averages of $AR\left(t\right)$ and ${\gamma }_{QHF}^2\left(t\right)$, are depicted in Fig. 1, for $0\le t\le 2000$ sec, showing that the gap between them is already large from the start. The properties of the filter are further investigated in Figure 2 and Figure 3 that depict the time variations of $AR$ for various initial estimates of the quaternion, $\widehat{\mathbf{q}}\left(0\right)$, and various initial values of the matrix $\tilde{P}\left(0\right)$, respectively. This is done for the case $({\sigma }_{ϵ},{\sigma }_b)=(0.1,0.1)$, where the disturbance attenuation performance is best. It appears that the transient of $AR$ is strongly shortened when $\widehat{\mathbf{q}}\left(0\right)$ is close to the true quaternion. On the other hand, the steady states are relatively close. Figure 3 further shows the lack of sensitivity of AR to $\tilde{P}\left(0\right)$. These properties stem from the independence of the estimator’s gain computations from the state and are analogous to the convergence properties of covariances in Kalman filters for linear systems.

Figure 4 depicts the MC-means and the MC-standard deviations of the four components of the quaternion estimation error for ${\sigma }_{ϵ}=0.001{\displaystyle \frac{rad}{\sqrt{sec}}}$ and ${\sigma }_b=0.1rad$. The means are close to zero and the standard deviations show satisfying estimation performances, around $3mrad$. Figure 5 presents the time histories of the MC-mean and MC-standard deviation envelop of the angular estimation error, $\delta \varphi$. Albeit oscillating with an amplitude of $0.06\left[deg\right]$ around $0.08\left[deg\right]$, $\delta \varphi$ shows good performances given the measurement noise level ${\sigma }_b$ of about 5 degrees.

Extensive simulations were run to compare the performances of the QHF with those of a quaternion multiplicative EKF (MEKF). In the MEKF, the (quadratic in ${\mathbf{q}}_t$) measurement equation model is linearized, and the filter statistics are matched to the true noise levels. Table 4 shows the time averages, computed on single runs of 2000 seconds, of the quaternion additive estimation error norm in the MEKF (left) and in the QHF (right). These values provide sensible measures of the estimation biases. In addition, the values of the time standard deviations are provided for both filters (in parenthesis). The QHF consistently provides smaller biases than the MEKF. This is explained by the linearization effects impairing the MEKF, whereas the QHF is free of linearization. On the other hand, the MEKF provides smaller standard deviations than the QHF, as expected since the MEKF is a (approximate) minimum variance estimator. Yet, for a given value of ${\sigma }_b$, the gap between them decreases as ${\sigma }_{ϵ}$ increases and becomes negligible for large ${\sigma }_{ϵ}$. Table 4 further shows the low sensitivity of the QHF standard deviations with respect to the parameters ${\sigma }_{ϵ}$ and ${\sigma }_b$. This property is partially explained by the $H_{\infty }$ approach since the gain computations are independent of ${\sigma }_b$ and ${\sigma }_{ϵ}$ per se. Yet they do show some dependence on the level of the noises because the data itself is noisy.

Both filters were tested in cases where the true noise variances were unknown. This might occur as a result of undetected sensor failures or jamming. In Case A, the ${\sigma }_b$ in the MEKF was set to ten times its true value. In Case B, it was lowered to one-tenth of the true value. The results are shown in Figure 6 and Figure 7, respectively. In case A, the MEKF is very slow to converge while, in case B, it converges quicker but to a noisier steady state. The QHF, on the other hand, provides essentially the same performances in both cases, with slight variations due to the data noisiness. In both cases, the QHF outperforms the unmatched MEKF.

5. Conclusion

In this work, stochastic $H_{\infty }$ filtering was applied to the development of novel quaternion attitude estimators from rate gyro and line-of-sight (LOS) measurements. A key assumption is that the variance of the noise affecting the various measurements is unknown and modeled as a disturbance. The estimators compute the quaternion while attenuating the transmission from the noise variance to the estimation error. The $H_{\infty }$ filters involve the solution of a set of differential and algebraic linear matrix inequalities. A remarkable property of the resulting gains computations is that they are independent of the estimated quaternion, in the case of measurement white noise. Extensive Monte-Carlo simulations were run showing that the proposed filter performs well from the standpoint of attitude estimation per se in a wide range of gyro and LOS noise variances. The guaranteed disturbance attenuation level seems to be slightly dependent on these variances since the gain depends on the measurements. The actual disturbance attenuation level seen in the simulations is better than the guaranteed one, by up to one order of magnitude. It improves when the noise level increases and is the worst for (ideal) noise-free sensors. This fact is in agreement with the theory and illustrates the conservative nature of the $H_{\infty }$ approach. When ${\sigma }_{ϵ}$ is the sole unknown the $H_{\infty }$ filter produces lower MC-means than a standard quaternion multiplicative Kalman filter. When both ${\sigma }_{ϵ}$ and ${\sigma }_b$ are unknown, the $H_{\infty }$ filter shows similar MC-means as a multiplicative Kalman filter. When matched the MEKF shows lower MC-standard deviations of the estimation errors than the $H_{\infty }$ filter. The higher the level of the noise, the less obvious the advantage of the Kalman filter. Furthermore, the $H_{\infty }$ filter gain is less sensitive to perturbations than the MEKF gains, in particular to initial estimation errors. This attractive feature is emphasized by comparing the $H_{\infty }$ filter’s performances with those of unmatched Kalman filters. When provided with too high or too low noise variances, the MEKF was outperformed by the $H_{\infty }$ filter, which essentially delivers identical performances within a wide range of noise variances.