1 Introduction↩︎

Nonlinear state estimation for discrete-time stochastic systems is a classical problem with broad applications in control, signal processing, and autonomous systems. A standard approach is to linearize the model, leading to the extended Kalman filter (EKF) [1]. However, EKF is effective only for mildly nonlinear systems, since it relies on local linearization based on a first-order Taylor expansion. To achieve a more accurate treatment of nonlinearities while retaining moderate computational complexity, sigma-point Kalman filters, including the unscented Kalman filter (UKF) [2], the cubature Kalman filter (CKF) [3], and the Gauss–Hermite Kalman filter [4], [5], have become standard tools for nonlinear estimation [6].

In many practical applications, however, the estimation problem is characterized not only by nonlinearities, but also by model uncertainty [7], [8]. Such uncertainty may arise from imprecisely known model parameters, unmodeled dynamics, or unknown sensor drifts [9], [10]. Under these conditions, standard sigma-point filters may suffer substantial performance degradation, since they are built upon a nominal model. Existing robust sigma-point filtering methods are often specialized to particular types of uncertainty, such as outliers or heavy-tailed non-Gaussian noise [11], [12], rather than to more general perturbations affecting the underlying probabilistic model. A natural framework for handling this broader form of uncertainty is provided by minimax robust estimation paradigms, which have been extensively studied in the linear setting [13]–[19]. In that framework, the estimator is designed according to the least favorable model selected within a prescribed ambiguity set, centered at the nominal model, with radius quantifying the level of uncertainty. The few existing works on robust minimax state estimation for nonlinear models are based on EKF-type approximations [20]–[22]. Moreover, the latter mainly focus on the construction of the robust estimator, while they do not characterize the associated least favorable model. It is worth noting that characterizing the least favorable model is equally important, since it provides the basis for generating adversarial data and assessing whether such a type of uncertainty is realistic.

Extending a minimax approach in [13], [23] to the nonlinear case is highly nontrivial. Indeed, the challenge is not merely to counteract uncertainty by enlarging the covariance matrices of the noise processes in a standard nonlinear algorithm, but rather to provide a clear interpretation of a robust estimator, as a minimizer of a minimax game, obtained via a suitable approximation. Moreover, unlike in the linear case, it is difficult to derive a state-space representation of the least favorable model, i.e. the maximizer, which in turn makes the generation of adversarial data challenging.

In this paper, we consider a minimax game for nonlinear state estimation under model uncertainty in the measurement equation which extends the linear framework recently proposed in [23]. The characterization of the solution to this problem is not tractable. We show that replacing the center of the ambiguity set with a sigma-point approximation allows us to characterize a minimizer, thereby yielding a robust sigma-point-like filter.

Another contribution of the paper concerns the characterization and assessment of the least favorable model, i.e. the maximizer of the corresponding problem. The latter does not admit a state-space representation, which makes direct simulation highly nontrivial. To address this challenge, we propose a Markov chain Monte Carlo (MCMC) scheme for generating adversarial data from the least favorable model. In particular, we design a proposal state-space model that remains close to the target density while being easy to sample from. We further establish the convergence of the resulting Markov chain.

The numerical experiments confirm the convergence of the proposed MCMC algorithm and show that the robust sigma-point filter exhibits the expected optimal behavior relative to the standard nonlinear filter under adversarial data. In addition, we consider a more pragmatic setting with parametric model mismatch and different sensor non-idealities scenarios, where it is not guaranteed that the actual model belongs to the prescribed ambiguity set. Even in this case, the proposed approach compares favorably with the existing state-of-the-art robust nonlinear filters, such as the robust EKF [20], [21], the maximum correntropy UKF [12], and the maximum entropy particle filter [24], in terms of both estimation performance and computational efficiency. Moreover, the simulation results show that the robust EKF, which relies on a simpler approximation of the nonlinear dynamics, yields poor estimation performance. This illustrates that, in practice, it is important not only to guarantee robustness against uncertainty, but also to retain a good approximation of the nonlinear dynamics. This paper builds upon some preliminary results presented in the conference paper [25].

The remainder of the paper is organized as follows. Section 2 introduces the problem formulation. Section 3 introduces the minimax problem and derives a saddle point solution. Section 4 presents the MCMC scheme. Section 5 reports the numerical examples. Finally, Section 6 concludes the paper.

Notation. \([0,N]_{\mathbb{Z}}\) denotes the interval of integers between \(0\) and \(N\). \(z\sim f(z)\) means that the random vector \(z\) is distributed according to the probability density \(f(z)\); \(f(z)= \mathcal{N}(\mu ,P)\) means that the probability density \(f(z)\) is Gaussian with mean \(\mu\) and covariance matrix \(P\). Given a sequence of random variables \(a_r\), with \(r\in\mathbb{N}\), \(a_r\overset{p}{\longrightarrow} a\) means that the sequence \(a_r\) converges in probability to the random variable \(a\). Given a symmetric matrix \(P\), \(P>0\) and \(P\geq 0\) mean that \(P\) is positive definite and semi-definite, respectively. Moreover, \(|P|\) and \(\operatorname{tr}(P)\) denote the determinant and the trace of \(P\); \(I_n\) denotes the identity matrix of dimension \(n\); \(A^\top\) denotes the transpose of matrix \(A\). Finally, throughout the paper, we omit the arguments of probability density functions, i.e., \(\psi_t\) instead of \(\psi_t(y_t|x_t)\), when the context is clear in order to keep the notation concise.

2 Problem Formulation↩︎

We consider the nominal discrete-time nonlinear state space model: \[\label{nomi95mod} \left\{\begin{align} x_{t+1}&= f(x_t) + Bv_t \\ y_t &= h( x_t) + Dv_t.\\ \end{align}\right.\tag{1}\] where \(x_t\in\mathbb{R}^n\) is the state, \(y_t\in\mathbb{R}^m\) is the observation, \(v_t \in \mathbb{R}^{n+m}\) is white Gaussian noise (WGN) with covariance matrix equal to \(I_{n+m}\), \(x_0\) is Gaussian distributed. Matrices \(B \in \mathbb{R}^{ n \times (n+m)}\) and \(D \in \mathbb{R}^{ m \times (n+m)}\) are full row rank and such that \(BD^\top=0\). In plain words, the noise processes \(Bv_t\) and \(Dv_t\) are assumed independent. Moreover, we assume that \(v_t\) is independent from the initial state \(x_0\). Consider the nominal state space model (1 ). The latter over the finite time interval \([0,N]_{\mathbb{Z}}\) is characterized by the nominal density of \(Z_N :=\left[\begin{array}{lllllll}{x^\top_{0}} & {\ldots} & {x^\top_{N+1}}, & {y^\top_{0}} & {\ldots} & {y^\top_{N}}\end{array}\right]^\top\) which is \[\label{nomi1} p\left(Z_{N}\right)= p_0\left(x_{0}\right) \prod_{t=0}^{N} {p_{t}}\left(x_{t+1} | x_{t}\right) \psi_{t}\left(y_{t} | x_{t}\right) ,\tag{2}\] where \(p_0(x_0) = \mathcal{N}(\hat{x}_0,\tilde{P}_0)\) denotes the density of \(x_0\), and \[\begin{align} \label{pxt1t} p_t(x_{t+1}|x_t) & = \mathcal{N} (f(x_t), ~ BB^\top),\\ \psi_t(y_t|x_t) & = \mathcal{N} (h(x_t), ~ DD^\top). \nonumber \end{align}\tag{3}\] We assume that the actual density of \(Z_N\) takes the form \[\tilde{p}\left(Z_{N}\right)=p_0\left(x_{0}\right) \prod_{t=0}^{N} {p_{t}}\left(x_{t+1} | x_{t}\right) \tilde{\psi}_{t}\left(y_{t} | x_{t}\right),\nonumber\] i.e. the measurement model – described by the conditional density \(\tilde{\psi}_t\) – is different from the nominal one. We can measure the mismatch between the actual and the nominal measurement model at time \(t\) by the conditional Kullback–Leibler (KL) divergence \[\begin{align} \mathcal{D}(\tilde{\psi}_{t},& {\psi}_{t}) : =\iint \tilde{\psi}_t(y_t | x_t) \tilde{p}_t(x_t|Y_{t-1})\ln\left(\frac{\tilde{\psi}_t}{\psi_t}\right)d y_t d x_t \end{align}\] where \(\tilde{p}_t(x_t|Y_{t-1})\) is the actual a priori density of \(x_t\) given \(Y_{t-1}:=[\,y_0^\top\,\ldots\, y_{t-1}^\top\,]^\top\). We assume that \(\tilde{\psi}_t\) belongs to the ambiguity set \[\label{ambi95st} {\mathcal{B}}_t:=\left\{\, \tilde{\psi}_{t}s.t.\mathcal{D}(\tilde{\psi}_{t}, {\psi}_{t}) \leq c_{t}\right\},\tag{4}\] where \(c_t\) is the so-called tolerance. In plain words, we place an upper bound on model uncertainty at each time step, ensuring that uncertainty is not concentrated on specific time steps where the estimator is more vulnerable. This type of ambiguity set naturally arises when the nominal model 1 is estimated from data, see [23] and [26]. Notice that, \(\mathcal{B}_t\) is a convex set, moreover it contains conditional densities which are not necessarily Gaussian. Our aim is to develop a general approach for sigma point Kalman filtering, which provides a robust a posteriori state estimator of \(x_{t}\) given \(Y_t\), in the case that the actual measurement model is unknown and different from the nominal one in (1 ). When the functions \(f(x)\) and \(h(x)\) are linear in \(x\), the problem has been addressed in [23], where a robust filter can be derived and interpreted within a dynamic minimax game framework. In that setting, the robust estimator can be explicitly characterized as a Kalman-like filter with an enlarged estimation error covariance, and the associated least-favorable model can also be derived in closed form as a state space model, which allows the systematic generation of data from such a model. However, when the system is nonlinear, although one may heuristically introduce robustness in a given estimation algorithm by enlarging the estimation error covariance, the interpretation of this modification is far from trivial. In particular, it is generally unclear which underlying dynamic game leads to such a solution, and how the corresponding least-favorable model should be constructed in order to generate adversarial data consistent with that robustness interpretation. In the following we address these challenges by providing a dynamic game interpretation and characterizing the associated least-favorable model in the nonlinear setting. Throughout the paper, to ease the exposition, we will consider the case in which the nonlinear model (1 ) is time-invariant and autonomous, i.e. there is no exogenous input \(u_t\). However, the results we present can be extended to the time-varying case where \(f_t(x_t,u_t)\) and \(h_t(x_t,u_t)\) depend on an input \(u_t\), see Remark [rem95gen] for more details.

3 Robust sigma point Kalman filtering↩︎

Our framework relies on a dynamic minimax game composed by two players: one player, i.e. the state estimator, minimizes the variance of the state estimation error, while the other one, i.e. the nature, selects the least favorable model belonging to a set of possible models about the nominal model (1 ). The latter is the prescribed ambiguity set (4 ). In this setting, the robust estimation problem characterizing the estimator of \(x_t\) given \(Y_t\) is defined as \[\label{minimaxrsu2} (\tilde{\psi}^\star_t, g^\star_t ) =\arg \underset{g_t \in \mathcal{G}_{t}}{\mathrm{min}}\max_{\tilde{\psi}_{t} \in \mathcal{B}_{t}} J_t(\tilde{\psi}_t,g_t)\tag{5}\] where \[\label{JtU} \begin{align} J_t(\tilde{\psi}_t,g_t)=&\frac{1}{2} \iint \left\| x_{t}-g_{t}\left(y_{t}\right) \right\|^{2} \tilde{\psi}_{t}(y_t | x_t) \tilde{ p}_{t}\left(x_{t} | Y_{t-1}\right) d y_{t} d x_{t} \end{align}\tag{6}\] is the variance of the estimation error under the actual model; \(\mathcal{G}_t\) is the convex set of estimators with finite second order moments which are square integrable with respect to densities \(\tilde{\psi}_t \tilde{ p}_{t}(x_{t} | Y_{t-1})\) and \(\tilde{\psi}_t\) satisfies the condition \[\iint \tilde{\psi}_{t}(y_t | x_t)\tilde{ p}_{t}\left(x_{t} | Y_{t-1}\right) d y_{t} d x_{t}=1. \label{Iuukf}\tag{7}\] Our aim is to characterize a saddle point solution to 5 . Solutions of this type of problems are typically non-unique, see [13], [26] for instance. In view of the Von Neumann’s minimax theorem [27], there exists a saddle point \((\tilde{\psi}^{\star}_t,~ g^{\star}_t)\) such that \[J_t(\tilde{\psi}_t,~ g^{\star}_t) \leq J_t(\tilde{\psi}^{\star}_t,~ g^{\star}_t) \leq J_t(\tilde{\psi}^{\star}_t,~ g_t),\nonumber\] since the corresponding sets \(\mathcal{B}_t\) and \(\mathcal{G}_t\) are convex and compact. The next result characterizes the structure of the maximizer of (5 ), i.e. the least favorable model.

Let \[\begin{align} \bar p_t(w_t|Y_{t-1})&= \psi_t(y_t|x_t)\tilde{p}_t(x_t|Y_{t-1})\nonumber\\ \tilde{p}_t(w_t|Y_{t-1})&= \tilde{\psi}_t(y_t|x_t)\tilde{p}_t(x_t|Y_{t-1})\nonumber \end{align}\] denote the pseudo-nominal and the actual conditional density of \(w_t:=[\,x_t^\top\; y_t^\top\,]^\top\) given \(Y_{t-1}\), respectively. The idea is to characterize a saddle point solution by mapping Problem 5 into a minimax game whose players are \(g_t\) and \(\tilde{p}_t(w_t|Y_{t-1})\) \[\begin{align} \label{J95eqpb95minimax} \underset{g_t \in \mathcal{G}_{t}}{\mathrm{min}}\max_{\tilde{p}_{t} \in \tilde{\mathcal{B}}_{t}}\bar J_t(\tilde{p}_t,g_t) \end{align}\tag{10}\] where \[\begin{align} \bar J_t(\tilde{p}_t,g_t):=&\frac{1}{2} \iint \left\| x_{t}-g_{t}\left(y_{t}\right) \right\|^{2} \tilde{ p}_{t}\left(w_{t} | Y_{t-1}\right) d w_{t}; \end{align}\] the ambiguity set \(\tilde{\mathcal{B}}_t\) is obtained using the correspondence between \(\tilde{\psi}_t(y_t|x_t)\) and \(\tilde{p}_t(w_t|Y_{t-1})\) given by the optimality condition 8 : \[\label{ball95uukf2} \tilde{\mathcal{B}}_t:=\left\{\, \tilde{p}_{t}(w_t|Y_{t-1})s.t.\mathcal{D}_{KL}(\tilde{p}_{t},\bar {p}_{t}) \leq c_{t}\right\}\tag{11}\] where \[\begin{align} \mathcal{D}_{KL}(\tilde{p}_{t},\bar {p}_{t})=\int\ln\left(\frac{\tilde{p}_t}{\bar p_t}\right)\tilde{p}_t dw_t. \nonumber \end{align}\] Unfortunately, in this case it is not possible to characterize a saddle point solution to (10 ) because the pseudo-nominal density \(\bar p_t(w_t|Y_{t-1})\) is not Gaussian. Thus, we construct an approximation of \(\bar p_t(w_t|Y_{t-1})\) using the same approximation mechanism exploited in the traditional sigma point Kalman filter, [6]. First, we introduce the definition of sigma points, which is used to compute an approximation of the expectation of a nonlinear transformation of a Gaussian random vector.

There are many ways to define the sigma points and the weights (see [6]). Different choices result in various nonlinear Kalman filters, among which the most popular are:

• Unscented Kalman filter (UKF). The sigma points are obtained through the unscented transformation: \[\nonumber \begin{align} \sigma_{i}(\hat{x},P)=\left\{\begin{array}{ll} \hat{x} + \sqrt{\lambda +n} (\sqrt{{P}} )_i, & if1\leq i\leq n \\ \hat{x} - \sqrt{\lambda +n} (\sqrt{{P}} )_{i-n}, & ifn+1\leq i\leq 2n\\ \hat{x}, & ifi=2n + 1 \end{array}\right. \end{align}\]where \((\sqrt{ P })_i\in\mathbb{R}^n\) is the \(i\)-th column of \(\sqrt{ P }\) which is a square root matrix of \(P\). The corresponding weights are \[\begin{align} \label{W95ukf} W_m^{i}&=\lambda / (n+\lambda), ~ W^i_c = W^i_m +1- a^2 +b, \quad ifi=2n + 1 \nonumber\\ W_m^i&= W_c^i = 1 / (2(n+\lambda)), \quad if1\leq i \leq 2n \end{align}\tag{12}\] where \(\lambda = a^2 (\kappa +n) - n\), and the parameters \(a\), \(b\) and \(\kappa\) can be chosen as suggested in [28].

• Cubature Kalman filter (CKF). The sigma points are generated using the spherical cubature transformation: \[\nonumber \begin{align} \sigma_{i}(\hat{x},P)=\left\{\begin{array}{ll} \hat{x} + (\sqrt{{n P}} )_i, & if1\leq i\leq n \\ \hat{x} - (\sqrt{{n P}} )_{i-n}, & ifn+1\leq i\leq 2n \end{array}\right. \end{align}\] and the corresponding weights are \(W^i_c = W_m^i = {1}/{2n}\).

• Gauss-Hermite Kalman filter. The sigma points are generated by the Gauss-Hermite moment transformation \[\nonumber \begin{align} \sigma_{i}(\hat{x},P)=\hat{x} +\sqrt{{P}} \lambda^i, \quad i= 1\ldots q^n \end{align}\] where \(\lambda^i \in \mathbb{R}^{n}\) is the \(i-\)th vector of the set formed by the \(n-\)dimensional Cartesian products of the roots, say \(\nu_{k}\) with \(k=1, \ldots, q\), of the Hermite polynomial of order \(q\), say \(H_q(\nu)\). The corresponding weights \(W_m^i = W_c^i\) are formed as the products of the \(n\) terms \[\begin{align} \frac{q!}{q^2\left(H_{q-1}\left(\lambda_{k}^i\right)\right)^2} \end{align}\] where \(\lambda_{k}^i\), with \(k = 1, \ldots, n\), denotes the \(k\)-th element of \(\lambda^i\).

The approximation \(\bar {p}_{t}^a(w_{t} | Y_{t-1})=\mathcal{N}(m_t,K_t)\) of \(\bar {p}_{t}(w_{t} | Y_{t-1})\) in standard sigma point Kalman filters is obtained from the sigma points corresponding to the approximation \(\tilde{p}_{t}(x_{t} | Y_{t-1})\simeq\mathcal{N}(\hat{x}_t,P_t)\): \[\begin{align} \label{defmK}m_t&= \left[\begin{array}{c} \hat{x}_t \\ m_{y_t}\end{array}\right],\; K_t= \left[\begin{array}{cc} P_t & K_{x_t y_t} \\ K_{y_t x_t} & K_{y_t}\end{array}\right], \end{align}\tag{13}\] \[\label{mmy} \begin{align} m_{y_t} &= \sum_{i=1}^{p} W_m^i h(\mathcal{X}^i_t)\\ K_{y_t} &= \sum_{i=1}^{p} W_c^i(h(\mathcal{X}^i_t) - m_{y_t})(h(\mathcal{X}^i_t) - m_{y_t})^\top + D D^\top \\ K_{x_t y_t} &= \sum_{i=1}^{p} W_c^i(\mathcal{X}^i_t - \hat{x}_t)(h(\mathcal{X}^i_t) - m_{y_t})^\top\\ \mathcal{X}^i_t&=\sigma_i(\hat{x}_t, P_t), \quad i=1\ldots p. \end{align}\tag{14}\] Applying such approximation in 10 we obtain the approximate problem \[\label{minimaxa95up} \hat{x}_{t|t} = \underset{g_t \in \mathcal{G}_{t}}{\mathrm{argmin}}\max_{\tilde{p}_{t} \in {\tilde{\mathcal{B}}}^a_{t}} {\bar J_t}(\tilde{p}_{t} ,g_t)\tag{15}\] where \(\tilde{\mathcal{B}}^a_t\) is obtained from (11 ) replacing \(\bar{p}_{t}(w_t|Y_{t-1})\) with \(\bar{p}_{t}^a(w_t|Y_{t-1})\). In plain words, the approximation regards the center of (the ball describing) the ambiguity set. Before characterizing a saddle point solution to Problem 15 , we must specify how \(\hat{x}_t\) and \(P_t\), i.e. the mean and covariance matrix of the a priori conditional density at time \(t-1\) used in the prediction stage, are obtained. These quantities are computed from: i) the mean and covariance matrix, \(\hat{x}_{t-1|t-1}\) and \(\tilde{P}_{t-1|t-1}\), of the posterior density \(\tilde{p}_{t-1}(x_{t-1}|Y_{t-1})\) at time \(t-1\), which is assumed to be reasonably approximated as a Gaussian density; and ii) the state transition model \(p_t(x_{t+1}| x_t)\). Since \(p_t(x_{t+1}| x_t)\) is not affected by uncertainty, i.e. it represents both the nominal and the actual state process model, we can employ the same approximation used in the standard sigma-point Kalman filter to compute \(x_t\) and \(P_t\), [6]: \[\begin{align} {\mathcal{X}}_{t-1|t-1}^i&=\sigma_i(\hat{x}_{t-1|t-1},\tilde{P}_{t-1|t-1}),\quad i=1\ldots p\nonumber\\ \hat{\mathcal{X}}_{t}^i&=f( {\mathcal{X}}_{t-1|t-1}^i),\quad i=1\ldots p\nonumber\\ \hat{x}_t&=\sum_{i=1}^{p} W_m^i\hat{\mathcal{X}}_t^i,\; \nonumber\\ \label{defPP}{ P}_{t}&=\sum_{i=1}^{p} W_c^i(\hat{\mathcal{X}}_{t}^i-\hat{x}_{t})(\hat{\mathcal{X}}_{t}^i- \hat{x}_{t})^\top + B B^\top. \end{align}\tag{16}\] We are ready to characterize a saddle point solution to Problem 15 .

The theorem above not only characterizes the robust estimator, but also shows that the least favorable posterior density, i.e. the one corresponding to 10 , can be approximated by the Gaussian density \(\tilde{p}_t^a(x_t| Y_t)\) in (19 ). The latter allows us to conclude that the sigma-point approximation is valid over the entire time horizon, and this can be shown by an inductive argument: if the posterior \(\tilde{p}_{t-1}(x_{t-1}|Y_{t-1})\) at time \(t-1\) can be approximated as Gaussian, then the prior density \(\tilde{p}_t(x_t| Y_{t-1})\) can be approximated by means of the sigma point transformation, and the theorem ensures that the updated posterior at time \(t\) remains approximately Gaussian. Consequently, at the next time step, the approximation of \(\tilde{p}_{t+1}(x_{t+1} \mid Y_t)\) can again be constructed using sigma points.

The resulting robust sigma point Kalman filter is outlined in Algorithm 1. Note that, the computation of \(\theta_{t|t}\) in Step [AL95theta] can be efficiently computed through a bisection method, see [29]. Finally, in the limit case \(c_t = 0\), i.e. the absence of uncertainty in the measurement model at time \(t\), we have that \(\theta_{t|t}=0\) and \(\tilde{P}_{t|t}=P_{t|t}\), thus the robust sigma point Kalman filter coincides with the standard sigma point Kalman filter.

In view of Theorem [th2], the approximate least favorable density takes the following form: \[\begin{align} \label{psi095t} &\tilde{\psi}_{t}^{0}(y_t|x_t)=\frac{1}{M_{t}} \exp \left(\frac{ \theta_{t|t}}{2} \left\| x_{t}- \hat{x}_{t|t}\right\|^{2}\right) \psi_{t}(y_t|x_t) \end{align}\tag{20}\] where \(\hat{x}_{t|t}\) and \(\theta_{t|t}\) are the state filter and the risk sensitivity parameter obtained by the robust sigma point Kalman filter; \[\begin{align} \label{M95T3} & M_{t}=\iint \exp \left( \frac{\theta_{t|t}}{2} \left\| x_{t}- \hat{x}_{t|t} \right\|^{2}\right) \psi_{t} \tilde{p}^a_{t}\left(x_{t} | Y_{t-1}\right) d y_{t} d x_{t} \end{align}\tag{21}\] where \[\tilde{p}^a_t(x_t|Y_{t-1})= \int \bar p^a_t(w_t|Y_{t-1}) d y_t =\mathcal{N}(\hat{x}_t,P_t)\] is the marginal obtained from the approximation of the pseudo-nominal density \(\bar p_t(w_t|Y_{t-1})\). It is worth noting that the optimality of \(\tilde{\psi}^0_t\) is guaranteed by Proposition [l95lfm2] since by direct substitution we have \[\begin{align} \mathcal{D}( \tilde{\psi}_t^0, \psi_t) =\mathcal{D}_{KL}(\tilde{p}_{t},\bar {p}^a_{t})=c_t \nonumber \end{align}\] where \(\tilde{p}_t\) denotes the maximizer of 15 . We now investigate the conditions under which the approximate maximizer is well defined.

We conclude that \(( \hat{x}_{t|t},\tilde{\psi}_{t}^{0})\) represents an approximate solution to Problem (5 ) in the sense that: \[\begin{align} J_t(\tilde{\psi}_t,~ \hat{x}_{t|t}) \leq J_t(\tilde{\psi}^{0}_t, \hat{x}_{t|t}) \approx J_t(\tilde{\psi}^{\star}_t, g^{\star}_{t}) \leq J_t(\tilde{\psi}^{\star}_t, g_t), \end{align}\nonumber\] for any \(\tilde{\psi}_t \in \mathcal{B}_t\) and \(g_t \in {\mathcal{G}}_t\). Note that, the first of the above inequalities follows from Proposition [l95lfm2]. Thus, the approximate least favorable model over the finite time interval \([0,N]_{\mathbb{Z}}\) takes the form: \[\label{tilde95p0ZN} \begin{align} \tilde{p}^0(Z_N) \propto p_0(x_0) \prod_{t=0}^N p_t(x_{t+1} | x_t) \tilde{\psi}^0_t(y_t | x_t) \end{align}\tag{31}\] where the symbol \(\propto\) means that the two sides are equal up to a constant factor. In what follows, we will simply refer to (31 ) as least favorable model without specifying that it is an approximation.

4 Adversarial data generation↩︎

In order to assess the performance of the robust sigma point filter in the least favorable scenario, we need to develop a simulator for generating random samples from the least favorable density (31 ), i.e adversarial data. In the linear setup, it was shown that the least favorable model admits a state space representation over a finite time interval, [23]. Such result, however, cannot be exploited in this nonlinear setting. Indeed, a fundamental aspect in the linear setup is that the least favorable density is Gaussian, but the least favorable density in (31 ) is not Gaussian in general. Thus, it is not straightforward to draw samples from (31 ). Markov chain Monte Carlo (MCMC) algorithms are effective techniques for approximately sampling from complex probability densities in high-dimensional spaces. Thus, we use the Metropolis-Hastings (MH) algorithm [31] in order to tackle the problem. Our target density is \(\tilde{p}^0(Z_N)\), defined in (31 ). Suppose it is easy to generate a random sample from a proposal probability density \(\bar q( Z_N | Y_N^k)\) where \(Z_N^k=\{Y_N^k,X_{N+1}^k\}\), \(Y_N^k\) and \(X_{N+1}^k\) are the subvectors of \(Z_N^k\) containing the observation and state components, respectively. We consider the MH scheme outlined in Algorithm 2 which provides samples of \(Z_N\), following the target \(\tilde{p}^0\) generated by the proposal \(\bar q\).

Next, we introduce the proposal density \(\bar q( Z_N | Y_N^k )\), as well as how to evaluate \(\pi (Z_N)\) and \(q( Z_N | Y_N^k)\). In what follows, we assume that \(f:\, \mathbb{R}^{n}\rightarrow \mathbb{R}^n\) and \(h:\, \mathbb{R}^{n}\rightarrow \mathbb{R}^m\) are continuous functions.

4.1 Proposal density↩︎

We construct the proposal density relying on the least favorable state space model in [23] derived for the linear case. More precisely, given the observations \(Y_{N}^k\), we can compute \(\hat{x}_{t|t}\), \(\hat{x}_{t+1}\), \(L_t\) and \(\theta_{t|t}\), with \(t\in[0,N]_{\mathbb{Z}}\), using the robust sigma point Kalman filter outlined in Algorithm 1. Then, we linearize the nominal nonlinear model (1 ) along the estimated state trajectory by the robust sigma point Kalman filter: \[\label{linearized95model} \begin{array}{rl} x_{t+1} &= A_{t} x_t - A_{t} \hat{x}_{t|t} + f(\hat{x}_{t|t}) + B v_t\\ y_t &= C_{t} x_t - C_{t} \hat{x}_{t} + h(\hat{x}_t)+ D v_t\\ \end{array}\tag{33}\] where \(A_t\) and \(C_t\) are the Jacobian matrices:\[\label{AC}A_{t} := \left. \frac{\partial f(x)}{\partial x}\right |_{x=\hat{ x}_{t|t}},\; \; C_{t}: = \left. \frac{\partial h(x)}{\partial x}\right |_{x=\hat{ x}_{t}}.\tag{34}\] The least favorable model over the time interval \([0,N]_{\mathbb{Z}}\) corresponding to the ambiguity set about the linearized model (33 ) and with tolerance sequence \(\{c_t,\; t\in[0,N]_{\mathbb{Z}}\}\) is given by using the result in [23]: \[\label{LF95x2}\begin{align} \eta_{t+1} &=\bar {A}_t \eta_{t} + \bar a_t + \bar B_t \epsilon_t \\ y_{t} &=\bar{C}_t \eta_{t}+ b_t + \bar D_t \epsilon_t \end{align}\tag{35}\] where \[\label{barABCD}\begin{align} & \eta_t = \left[\begin{array}{c} x_t \\ e_{t-1} \\ Bv_{t-1} \end{array} \right],~ \bar a_t= \left[\begin{array}{c} a_t \\ 0 \\ 0 \end{array}\right], ~\epsilon_t = \left[\begin{array}{c} B v_t \\ \upsilon_{t} \end{array} \right], \\ &a_t = - A_{t} \hat{{x}}_{t|t} + f(\hat{{x}}_{t|t}), \\ &\bar{A}_{t}:=\left[\begin{array}{ccc} A_t & 0 & 0\\ 0 & \Delta_t - L_t F_t A_{t-1} & \Lambda_t \\ 0 &0 &0\\ \end{array}\right],\\ &\bar {B}_{t}:=\left[\begin{array}{cc} I_n & 0 \\ 0 & -L_t \Upsilon_t\\ I_n & 0\\ \end{array}\right],~~ b_t = - C_{t} \hat{{x}}_{t} + h(\hat{{x}}_t),\\ &\bar {C}_{t}:=\left[\begin{array}{llll} C_t & F_t A_{t-1} & F_{t} \end{array}\right], \quad \bar {D}_{t}:=\left[\begin{array}{cc} 0 & \Upsilon_t \end{array}\right],\\ &\Delta_t : = A_{t-1} - L_{t}C_t A_{t-1}, ~ \Lambda_t := I_n - L_{t} F_{t} - L_t C_t; \end{align}\tag{36}\] \(\varepsilon_{t} \in \mathbb{R}^{m}\) is normalized WGN; \[\begin{align}e_{t} &:= \left( \Delta_t - L_t F_t A_{t-1} \right) e_{t-1} + \Lambda_t B v_{t-1} - L_t \Upsilon_t \varepsilon_{t}; \\ O_t&:=\left((DD^\top)^{-1} -L_t^\top(\Omega_{t+1}^{-1}+\theta_{t|t}I_n) L_t \right)^{-1},\\ F_{t}& := - O_t L^\top_t (\Omega_{t+1}^{-1}+\theta_{t|t}I_n) (I_{ n} - L_{t} C_{t}); \end{align}\] \(\Upsilon_t\) is the Cholesky factor of \({O}_{t}\), i.e. \(\Upsilon_t \Upsilon^\top_t = O_t\); \(L_t\) is the Kalman filtering gain obtained by Algorithm 1. Moreover, \(\Omega^{-1}_{t}\) is computed by the following backward recursion: \[\label{omgea}\begin{align}\Omega^{-1}_{t} = A^\top_{t-1} F^\top_{t} O_t F_{t} A_{t-1} + \Delta_t^\top (\theta_{t|t} I_n + \Omega^{-1}_{t+1}) \Delta_t \end{align}\tag{37}\] with \(\Omega^{-1}_{N+1} = 0\). Notice that we can rewrite (35 ) as \[\label{linearized95lfm2} \begin{align} x_{t+1} &= {A_{t} x_t }+ a_t + B v_t \\ y_t &= {C_{t} x_t } + b_t + F_t A_{t-1} e_{t-1} +F_t B v_{t-1} +\Upsilon_t\varepsilon_t \end{align}\tag{38}\] where \(e_{t-1}\) is independent from \(v_{t-1}\) and it is Gaussian distributed with zero mean and covariance matrix \(\Pi_{e_t}\). The latter is the \(n \times n\) submatrix of \(\Pi_t\) obtained selecting rows and columns from \(n+1\) to \(2n\); \(\Pi_t\) is the solution to the Lyapunov equation (see [23]): \[\begin{align} \Pi_{t+1}=\Gamma_t \Pi_t \Gamma_t^\top +X_t\Xi X_t^\top \nonumber \end{align}\] where \[\begin{align} \label{a95lfm1} & \Gamma_t:=\left[\begin{array}{ccc} \Delta_t & -L_{t} F_{t} A_{t-1} &\Lambda_t \\ 0 & \Delta_t - L_t F_{t} A_{t-1} & \Lambda_t \\ 0 & 0& 0 \end{array}\right], \; X_t=\left[\begin{array}{cc} 0 & -L_t \Upsilon_t\\ 0 & -L_t \Upsilon_t\\ I_n & 0 \end{array}\right], \end{align}\tag{39}\] and \(\Xi\) is the block diagonal matrix with main blocks \(\{ B B^\top,I_m\}\). Therefore, the measurement model in (38 ) is characterized by the transition density of \(y_t\) given \(x_t\) and \(Y_N^k\) which is Gaussian with mean \(C_{t} x_t + b_t\) and covariance matrix \(Q^L_t\) where \[\begin{align} Q^L_{t} &:= F_{t} A_{t-1} \Pi_{e_{t-1}} A^\top_{t-1} F^{\top}_{t} + F_{t} B B^\top F^\top_{t} + \Upsilon_t \Upsilon^\top_t. \end{align}\] However, our aim is to develop a proposal density that effectively captures the essential characteristics of the target density. Notice that, \(A_t x_t +a_t\) and \(C_{t} x_t + b_t\) represent the first order Taylor expansion of \(f(x_t)\) and \(h(x_t)\) around \(\hat{ x}_{t|t}\) and \(\hat{ x}_{t}\), respectively. Then, a refined version of (38 ) is the one in which \(A_t x_t +a_t\) and \(C_{t} x_t + b_t\) are replaced by \(f(x_t)\) and \(h(x_t)\) , respectively. The corresponding “proposal” state space model is \[\label{linearized95lfm3} \begin{array}{rl} x_{t+1} &= f(x_t) + B v_t \\ y_t &= h(x_t) + F_t A_{t-1} e_{t-1} +F_t B v_{t-1} +\Upsilon_t\upsilon_t \\ \end{array}\tag{40}\] and the corresponding proposal density takes the form \[\label{defPROP} \bar q( Z_N|Y_N^k)= p_0( x_0)\prod_{t=0}^N p_t( x_{t+1}|x_t)\psi^L_t\left({ y_t | x_t ,Y_{t-1}^k}\right)\tag{41}\] where \[\label{psiL}{\psi^L_t( y_t| x_t, Y_{t-1}^k)} = \mathcal{N} (h(x_t), Q^L_{t}).\tag{42}\] To sum up, the state space model (40 ) is constructed as follows: first, given the observations \(Y_N^k\), we perform a forward sweep to compute \(\hat{x}_{t|t}\), \(\hat{x}_{t+1}\), \(L_t\) and \(\theta_{t|t}\) by means of the proposed robust sigma point Kalman filter, then we perform a backward sweep to compute \(\Omega^{-1}_t\) and finally we perform a forward sweep to compute the matrices characterizing the state space model. In this respect, it is easy to draw a proposal sample \(Z_N\) from the state space model (40 ). Finally, in view of Remark [remk1], we only need to evaluate \[\label{bar95p2} q({Z_N | Y_N^k})= \prod_{t=0}^N \psi^L_t\left({ y_t |x_t ,Y_t^k}\right).\tag{43}\] Accordingly, in view of (42 ) and (43 ), we obtain \[\begin{align} q( Z_N | Y_N^k ) = \prod_{t=0}^{N}\frac{1}{ \bar M_t} \exp\left( - \frac{1}{2}\|y_t - h(x_t) \|^2_{( Q^L_t)^{-1}} \right) \end{align}\nonumber\] where \(\bar M_t = \sqrt{(2 \pi)^{m}\left| Q^L_t\right|}.\)

4.2 Target density↩︎

In regard to the evaluation of target density \(\pi(Z_N)\), by (20 ) and (32 ), we have \[\begin{align} \pi(Z_N) = \prod_{t=0}^N \frac{1}{M_{t}} \exp \left(\frac{ \theta_{t|t}}{2} \left\| x_{t}- \hat{x}_{t|t}\right\|^{2}\right) \psi_{t}(y_t|x_t),\nonumber \end{align}\] so the challenging aspect regards the computation of \(M_t\) defined in (21 ). It is worth noting that it is not possible to find a closed form expression for the latter because of the presence of \(\psi_t(y_t|x_t)\), i.e. its mean is a nonlinear function of \(x_t\). Then, an approximation of \(M_t\) can be obtained through Monte Carlo integration: \[\label{t95M2} \hat{M}_{t,r}:= \frac{1}{r} \sum^{r}_{j=1} N_t (x_t^j),\tag{44}\] where \(N_t(x_t)\) has been defined in 24 ; \(x^{1}_t \ldots x^r_t\) are independent samples drawn from \(\tilde{p}_t^a(x_t| Y_{t-1}) = \mathcal{N}(\hat{x}_t, P_t)\). Note that, \(N_t(x_t)\) is a random variable obtained from a nonlinear transformation of \(x_t\) given \(Y_{t-1}\). The next theorem shows under which conditions \(N_t(x_t)\) has finite variance.

By Proposition [propos2minor], we have that the following condition approximately holds for \(r\) large \[\hat{M}_{t,r} \sim \mathcal{N}\left(M_{t},~ \frac{\hat{\sigma}^2_{t,r}}{ r} \right)\] and \[\hat{\sigma}^2_{t,r} := \frac{1}{r-1} \sum_{j=1}^{r}(N_t(x^j_t) - \hat{M}_{t,r})^2\nonumber\] is the sample variance estimator of \(\sigma^2_t\). Thus, it remains to design the number of samples \(r\) in such a way that \(\hat{M}_{t,r}\) approximates \(M_t\) with a certain accuracy. With a \(95\%\) confidence level, it holds that \[\frac{\sqrt{r}|\hat{M}_{t,r} - M_t|}{\hat{\sigma}_{t,r}} \leq 1.96,\] and thus the following inequality approximately holds \[\frac{|\hat{M}_{t,r} - M_t| }{M_t}\leq \frac{1.96\,\hat{\sigma}_{t,r}}{\sqrt{r}\hat{M}_{t,r}}:=\tau_r.\] Therefore, given the desired relative accuracy \(\tau^\star\) on the computation of \(M_t\) and an initial value for \(r\), we check whether \(\tau_r\leq \tau^\star.\) If not, we increase \(r\) until \(\tau_r\leq \tau^\star.\) Notice that, if we need to increase \(r\), then \(\hat{\sigma}_{t,r}\) and \(\hat{M}_{t,r}\) can be updated recursively using well known formulas.

In the special case \(h(x)\) is a linear function, i.e. \(h(x) = C x\), the normalization constant \(M_t\) can be computed in closed form.

Finally, the next results establish the convergence properties of the proposed algorithm.

5 Numerical Experiments↩︎

5.1 Worst-case performance↩︎

We assess the performance of the proposed filters in the worst-case scenario (i.e. using adversarial data). We consider the nominal nonlinear state space model 1 with \[\left\{ \begin{array}{rcl} f(x) & = & \begin{bmatrix} 0.2x_1+x_2+0.5\cos(2x_2)-1\\ 0.95x_2 \end{bmatrix} \\[1.2em] h(x) & = & \begin{bmatrix} x_1+0.1\left(\sqrt{1+x_2^2}-1\right)\\ 0.1x_1-0.1x_2 \end{bmatrix}. \end{array} \right.\] Moreover, \(x=[\, {x}_{1}^{\top} \; {x}_{2}^\top \,]^\top\), \(x_0\sim \mathcal{N} ( 0, I_2)\) and \(B=[\, I_2 \; 0\,]\), \(D=[\,0 \; I_2 \,]\). Notice that, the chosen nominal model satisfies the hypothesis on \(h\) in Theorem 2. In what follows, we consider the following nonlinear filters:

• UKF denotes the unscented Kalman filter [2] where the parameters in (12 ) are set as \(a = 0.5\), \(b = 2\), and \(\kappa = 1\);

• CKF denotes the cubature Kalman filter [3];

• RUKF denotes the robust sigma point Kalman filter of Section 3 with tolerance \(c =2 \cdot 10^{-4}\); the sigma points are obtained using the unscented transformation whose parameter setting is the same as the one of UKF;

• RCKF denotes the robust sigma point Kalman filter of Section 3 with tolerance \(c = 2 \cdot 10^{-4}\); the sigma points are obtained using the spherical cubature rule.

We evaluate their performance under the least favorable models corresponding to RUKF and RCKF introduced above. More precisely, we have generated \(M=500\) adversarial samples \(Z^k_N\) of length \(N=50\) from each least favorable model by means of the MH scheme proposed in Section 4. We set the initial value of \(r\) equal to \(100\) and the desired relative accuracy \(\tau^\star = 2 \cdot 10^{-3}\) for the computation of \(\hat{M}_{r,t}\). In this way, the cumulative relative error over the time interval \([0,50]_{\mathbb{Z}}\) is approximately equal to 10%. An upper bound equal to 1000 is imposed on \(r\) to control the computational time in edge cases. Fig. 3 shows the boxplot of the value of \(r\) required to compute \(\hat{M}_{r,t}\) (both rejected and accepted samples have been considered) by the two MH algorithms. We can see that the typical range of \(r\) is 250\(\div\)​400 and values greater than 600 are considered as outliers; thus, the chosen upper bound for \(r\) is adequate. Fig. 4 shows the corresponding acceptance rate for the two different simulators. Regarding the least favorable model of RUKF, i.e. the first subfigure in Fig. 4, we have considered the trajectories generated by two different samples of \(Z^0_N\). Notably, the algorithm exhibits a similar burn-in period and acceptance rate for the two trajectories. Regarding the least favorable model of RCKF, as in the second subfigure in Fig. 4, we have depicted only one trajectory, as we observed a similar behavior. Overall, the acceptance rate in both cases is approximately equal to \(80\%\), indicating that the proposal density has been properly designed and provides a good approximation to the target density.

Next, we assess the filters using the aforementioned dataset. More precisely, for each sample \(Z_N^k\) we extract the state trajectory \(X_N^k\) and the corresponding measurement trajectory \(Y_N^k\). Then, we evaluate the performance of the filters computing the mean squared error at time \(t\): \[\label{def95perf} \overline{\text{MSE}}_t = \frac{1}{M} \sum^M_{k=1} \| {{x}}_t^k-\hat{x}_{t|t}^k \|^2\tag{48}\] where \(x_t^k\) is the state at time \(t\) extracted from \(X_N^k\), while \(\hat{x}^k_{t|t}\) is the filtered state estimate obtained by the observations \(Y^{k}_N\). Fig. 5 and Fig. [FIG_RCKF] show the mean squared error for the four filters using the adversarial data generated by the least favorable model of RUKF and RCKF, respectively. The results highlight that, in each case, the robust filters outperform the standard filters. Moreover, the robust filters converge to a stable regime more rapidly, while also avoiding the pronounced fluctuations observed during the initial stage.

Finally, we consider a nominal model that is identical to the previous one, except that the state transition function \(f(x)\), which is now chosen as \[f(x)= \begin{bmatrix} 0.2x_1+x_2+0.5\cos(x_2)-1\\ 1.2x_2 \end{bmatrix}.\] Under adversarial data, the state trajectory sometimes has a drift. In this situation, we observed that the standard filters completely lose the ability to effectively track the true state, whereas robust filters can still maintain satisfactory tracking performance. More precisely, we evaluate the filter’s performance using the average mean squared error: \[\overline{\mathrm{MSE}} = \frac{1}{N}\sum_{t=1}^{N}\overline{\mathrm{MSE}}_t.\] The results are reported in Table 1. This finding is particularly meaningful in practical applications, since the presence of a drift may arise in scenarios such as tracking and navigation, where model uncertainty is often unavoidable.

5.2 Mass-spring system↩︎

In what follows, we assess the effectiveness of the proposed filters in a more challenging and realistic setting: the true model is not guaranteed to belong to the ambiguity set. We consider the mass-spring system depicted in Fig. 7. The reference position of the mass is chosen such that the restoring force is equal to zero.

Actual state space model. Let \(x=[p~s]^\top\) be the state vector, where \(p~[\texttt{m}]\) is the displacement and \(s~[\texttt{m/s}]\) is the velocity of the mass. The true nonlinear process model is described by the dynamics of a mass-spring system subject to friction, spring restoring force, and an external force. Such model can be written as (1 ) with \[\label{state95space95mass95spr95m} \begin{align} f(x) = x +\frac{T_s}{m} \left[\begin{array}{c} m s \\ -F_{f} - F_{s} \end{array}\right], ~ B = T_s \left[\begin{array}{ccc} \epsilon & 0 & 0 \\ 0 & \sqrt q / m & 0 \end{array}\right] \end{align}\tag{49}\] where \(m=1~\texttt{[kg]}\) is the wight of the mass, \(T_s=0.1~\texttt{[s]}\) is the sampling time, \(\epsilon = 0\) so that the relationship between displacement and velocity is exact. Specifically, \(F_f~[\texttt{N}]\) denotes the resistive force due to friction, \(F_s~[\texttt{N}]\) represents the restoring force of the spring; \(F_e~[\texttt{N}]\) is the external force, which is modeled as WGN with zero mean and variance \(q=0.25\), thus the second component of \(Bv_t\) represents \(T_sF_e/m\). The resistive force consists of static, Coulomb, and viscous friction components [32], i.e., \[\label{F95f} F_f = \alpha s + \eta(x),\tag{50}\] where \(\alpha=0.5~\texttt{[Ns/m]}\) is the viscous friction coefficient, and \(\eta(x)\) is defined by the following piecewise function: \[\label{friction95f} \eta(x)= \left\{ \begin{array}{ll} \mu_k m g \operatorname{sign}(s), & \text{for } |s|>0, \\[1mm] -k p, & \text{for } s=0 \text{ and } |p|\leq \mu_s m g/k, \\[1mm] -\mu_s m g \operatorname{sign}(p), & \text{for } s=0 \text{ and } |p|>\mu_s m g/k, \end{array} \right.\tag{51}\] where \(g=9.81~[\texttt{m/s}^2]\) is the gravitational acceleration. The parameters \(\mu_k\) and \(\mu_s\) denote the kinetic and static friction coefficients, respectively, which are given by \[\mu_k \sim \mathcal{U}(0.3,0.8), ~~~ \mu_s \sim \mathcal{U}(0.3,0.8).\] The restoring force is modeled by a hardening spring [32], where a small displacement increment beyond a certain threshold leads to a large increase in force, i.e., \[\label{spring95f} F_s = kp + k a^2 p^3,\tag{52}\] where \(k=10~\texttt{[N/m]}\) is the spring constant, and \(a \sim \mathcal{U}(0.01,0.05)\) is the nonlinear stiffness parameter characterizing the hardening effect. We assume the displacement \(p\) is measured by a sensor. More precisely, four different measurement non-idealities are considered:
\(\bullet\) sensor drift \[y_t = p_t + \tilde{\varepsilon}_t, \qquad \tilde{\varepsilon}_t \sim \mathcal{N}(0.5,1);\] \(\bullet\) intermittent observations \[y_t = \gamma_k p_t + \tilde{\varepsilon}_t, \qquad \tilde{\varepsilon}_t \sim \mathcal{N}(0,1)\] where \(\gamma_k\) is a binary random variable taking values in \(\{0,~1 \}\) with success probability \(\mathbb{P}[\gamma_{i,k}=1] = 0.8\);
\(\bullet\) outlier-contaminated noise \[y_t = p_t + \tilde{\varepsilon}_t, \qquad \tilde{\varepsilon}_t \sim \begin{cases} \mathcal{N}(0, 1), & \text{with probability } 0.9, \\ \mathcal{N}(0, 5), & \text{with probability } 0.1; \end{cases}\] \(\bullet\) nonlinearity dead zone \[y_t = \begin{cases} p_t + \tilde{\varepsilon}_t, & \text{if } |p_t + \tilde{\varepsilon}_t| \ge 0.1, \\ 0, & \text{if } |p_t + \tilde{\varepsilon}_t| < 0.1, \end{cases}\] where \(\tilde{\varepsilon}_t \sim \mathcal{N}(0,1)\).

Monte Carlo experiments. We consider Monte Carlo experiments with \(M=1000\) independent trials. In each trial, the actual process model is specified by independently sampling the uncertain parameters \(a\), \(\mu_k\), and \(\mu_s\) from the aforementioned uniform distributions. In addition, the corresponding measurement data are generated according to one of the four actual measurement models described above. Each trajectory is simulated over \(N=100\) time steps, with the initial state given by \(x_0 \sim \mathcal{N}\left([3\;0]^\top,\,0.1I\right)\). Since the sampling time is \(T_s=0.1\,\mathrm{s}\), this corresponds to a total simulation duration of \(10\) seconds. For each generated observation, we estimate the state using as nominal model (49 ) with nominal parameters \[a=0.03, \qquad \mu_k=0.6, \qquad \mu_s=0.5, \qquad \epsilon=10^{-8},\] where a small positive value of \(\epsilon\) is introduced to ensure the invertibility of the matrix \(BB^\top\), and nominal measurement model \[y_t = p_t + D v_t,\] where \(D=[\, 0 \; 0 \; 1\,].\) To this end,

we consider the following nonlinear filtering algorithms: the standard EKF, UKF, CKF, the bootstrap particle filter (PF) with \(N_p\) particles [6], and the proposed robust sigma point filters, i.e. RUKF and RCKF. Moreover, to further highlight the effectiveness of the proposed approach, we conduct a systematic comparison with several state-of-the-art robust nonlinear filtering methods available in the literature.

• REKF: the robust extended Kalman filter [20] with the following changes. The ambiguity set is chosen as in 4 . In addition, since this method is not directly applicable in the present setting due to the non-differentiability of \(f(x)\), we compute the numerical Jacobian matrix via a Newton approximation method.

• MC-UKF: the robust maximum correntropy UKF [12], where \(\sigma\) denotes the kernel width and the parameters in (12 ) are set identically to those of UKF and RUKF. The implementation follows the code publicly made available online by the authors.

• ME-PF: the robust maximum entropy particle filter proposed in [24] with \(1000\) particles. For each particle, the corresponding maximum entropy distribution is assumed to lie within a Kullback–Leibler divergence ball centered at the nominal distribution, with radius denoted by \(\tau\). The implementation follows the code publicly made available online by the authors.

All these filters are initialized with \(\hat{x}_0=[3\; 0]^\top\) and \(P_0=0.1 I\).

Results. The performance of these filters is assessed using the average of the mean squared error over the entire time horizon \(\overline{\text{MSE}}_t\) defined in (48 ). Table ¿tbl:tab:rmse95comparison? summarizes the performance of various filters under different measurement non-idealities and filter parameter settings. Note that when \(c = 0\), RUKF, RCKF, and REKF reduce to the standard UKF, CKF, and EKF, respectively. It is clear that all the proposed robust sigma point Kalman filters outperform the standard filters (i.e. UKF, CKF, EKF, and PF), regardless of the specific value of \(c\). For sigma point filters, when \(c\) is very small, their performance is comparable to that of the standard filters. As \(c\) increases, the performance improves; however, for large values of \(c\), a slight performance degradation is observed, as the robust filters become overly conservative. With regard to REKF, its performance is worse than that of the proposed estimator. It is worth noting that REKF, RUKF and RCKF are based on the same robust paradigm and therefore adopt the same strategy for coping with uncertainty, but differ in the way they approximate the nonlinear dynamics. This indicates that, in practice, it is important not only to achieve robustness against uncertainty, but also to preserve an accurate approximation of the nonlinear dynamics. Moreover, unlike the other robust filters, the performance of REKF fluctuates irregularly as \(c\) varies across almost all Monte Carlo experiments. This behavior is primarily attributable to a small number of realizations in which the estimated velocity enters a neighborhood of zero. Since \(\eta(x)\) exhibits a discontinuity on the hyperplane \(s=0\), the first-order Taylor approximation employed by REKF becomes highly inaccurate in these realizations, resulting in catastrophic estimation errors that significantly affect the overall performance metrics. Then, PF fails to achieve a good performance due to its high sensitivity to model uncertainties, even when a large number of particles is employed. The unsatisfactory performance of ME-PF can be attributed to the fact that its notion of robustness is achieved by maximizing the entropy of the particles ignoring the potential weakness of the robust estimator. Finally, the proposed robust filters also outperform MC-UKF. One possible reason is that the latter does not incorporate a time-varying strategy for adjusting the kernel width \(\sigma\) at each time step. To further investigate the impact of the parameter \(\alpha\) on the performance of RUKF, different values of \(\alpha\) are selected under the sensor drift scenario, as summarized in Table 2: changing parameter setting does not alter the fact that the RUKF performs better than the UKF. Finally, we compare the computational time of the robust filters. All simulations were implemented in MATLAB and executed on a Mechanical Revolution notebook equipped with an AMD R9-7845HX CPU and a GeForce RTX 4070 GPU. To ensure a fair comparison, fixed parameter values (i.e., \(c\), \(\tau\), or \(\sigma\)) were used for all the robust filters, chosen as the mean of the candidate values listed in Table [tbl:tab1]. As shown in Table 3, the proposed robust sigma-point Kalman filters have running time similar to MC-UKF and REKF, yet are significantly more computationally efficient than ME-PF.

6 Conclusion↩︎

We proposed a robust sigma-point filter and interpreted it as the minimizer of a minimax game whose ambiguity set is centered on a sigma-point approximation. We have shown that such an interpretation does not always hold, but it depends on whether function \(h\) satisfies certain properties. Moreover, we proposed a MCMC scheme for generating adversarial data from the least favorable model, thus allowing uncertainty to be assessed by examining the resulting realizations. The numerical results confirmed the convergence of the proposed MCMC scheme and validated the expected behavior of the robust sigma-point filter under worst-case conditions. Finally, they also showed the effectiveness of the proposed filtering approach compared to the existing robust nonlinear filters.

References↩︎

---

1. This paper was not presented at any IFAC meeting. Corresponding author: S. Yi.↩︎