---
abstract: 'In this paper, we propose a robust Kalman filtering framework for systems with probabilistic uncertainty in system parameters.  We consider two cases, namely discrete time systems, and continuous time systems with discrete measurements.  The uncertainty, characterized by mean and variance of the states, is propagated using conditional expectations and polynomial chaos expansion framework.  The results obtained using the proposed filter are compared with existing robust filters in the literature.  The proposed filter demonstrates better performance in terms of root mean squared error and rate of convergence.' 
author:
- 'Sunsoo Kim$^{1}$, Vedang M. Deshpande$^{2}$ and Raktim Bhattacharya$^{3}$ [^1][^2] [^3]'
bibliography:
- 'robustEstimation.bib'
title: '**Kalman Filtering with Probabilistic Uncertainty in System Parameters [^4]** '
---

Robust Kalman filter, estimation of uncertain systems, probabilistic uncertainty, polynomial chaos. 

Introduction
============

Robust filtering algorithms such as $\mathcal{H}_{2}/\mathcal{H}_{\infty}$ filters and robust Kalman filters, have been developed to address uncertainty in system models.  In the $\mathcal{H}_{2}$/$\mathcal{H}_{\infty}$ framework, filters are designed to minimize the impact of exogenous signals, i.e. process and sensor noise, on the estimation error [@geromel2002robust; @lacerda2011robust; @duan2013lmis; @green2012linear; @lewis2017optimal].  A robust Kalman filter is an extension of the well known Kalman filter, which can handle uncertainties in the system [@xie1994robustDT].  In this framework, the filter is designed to minimize an upper bound on the estimation error variance [@zhu2002design; @yang2002robust; @liu2018robust; @abolhasani2018robust; @xie1994robustCT; @shi1999robust], or the worst-case error variance [@sayed2001framework; @zorzi2016robust; @zorzi2018robust].  Our work falls in the category of robust Kalman framework. 

Existing robust Kalman filter algorithms can be categorized based on how system uncertainty is represented, which is assumed to be parametric.  The uncertainty is either represented as norm bounded parameter uncertainty [@xie1994robustDT; @zhu2002design; @yang2002robust; @liu2018robust; @abolhasani2018robust; @xie1994robustCT; @shi1999robust], or polytopic parametric uncertainty [@shaked2001new; @lewis2017optimal].  In this work, we model parametric uncertainty as random variables with known probability density function.  To the best of our knowledge, this is the first work on robust Kalman filtering with probabilistic system uncertainty. 

We present two robust Kalman filtering algorithms with probabilistic uncertainty in system parameters.  The first algorithm is for discrete-time (DT) system where the dynamics and measurements are both in discrete time.  The second algorithm is for continuous-time (CT) dynamical systems with discrete-time measurements.  In both these cases, mean and variance of uncertain states are calculated using a formulation based on conditional expectation.  For the CT system, we apply polynomial chaos (PC) framework to propagate the uncertainty. 

The rest of the paper is organized as follows.  We first present the problem formulation with uncertainty in CT and DT domain in §\[sec:Problem Formulation\] followed by a discussion on polynomial chaos framework in §\[sec:Polynomial Chaos Theory\]. § \[sec:Robust Kalman filter\] presents the proposed robust filter.  Simulation results obtained using the proposed filter are presented, and compared with the existing methods in §\[sec:result\].  Concluding remarks and future research directions are provided in §\[sec:Conclusion\]. 

Problem Formulation {#sec:Problem Formulation}
===================

The objective of filtering is to estimate the state-trajectory $\x(t)$ or $\x_k$ of a physical process in CT or DT, given noisy measurements.  The uncertainty in the system parameters, in the external excitation (process noise), and in the measurement errors (sensor noise), are all treated as probabilistic.  The model for the evolution of the state is assumed to be the following linear-time-invariant *stochastic* system,

$$\begin{aligned}
\text{CT:} &\  \xdot(t) = \Ap \x(t) + \Bp\w(t), \elab{sys_ct0} \\
\text{DT:} & \ \x_{k} = \Ap \x_{k-1} + \Bp \w_{k-1},  \elab{sys_dt0}\end{aligned}$$

where $\x\in\mathcal{R}^n$ is the model’s state vector, $\w\in\real^m$ is white noise with zero mean, $\Ap: \real^{d}\mapsto\real^{n\times n}$ and $\Bp: \real^{d}\mapsto\real^{n\times m}$ are system matrices, and $\param \in \real^{d}$ represents uncertain parameters in the system matrix with known probability density function.  We also assume, the initial condition for is a random variable with a given distribution.  Measurement from sensors is modeled as $$\begin{aligned}
 \y_k = \C \x_k + \n_k,\elab{y0}\end{aligned}$$ which maps the state $\x$ to the output space $\y$ and is corrupted by noise $\n$. In the output model, $\C\in\real^{m\times n}$ is deterministic and $\n$ is white noise with zero mean.  The variances of $\w$ and $\n$ are assumed to be $\vo{Q}$ and $\vo{R}$ respectively. 

The objective here is to determine the *unbiased* estimate of $\x(t)$ or $\x_k$ with *minimum error-variance*, using the model defined by and .  This is achieved by extending the formulation for standard Kalman filtering, to LTI systems with probabilistic uncertainty in system parameters, which is discussed in §\[sec:Robust Kalman filter\].  However before that, we briefly discuss the polynomial chaos framework that is used for propagation of uncertainty in CT systems. 

Polynomial Chaos Theory {#sec:Polynomial Chaos Theory}
=======================

Polynomial chaos is a deterministic framework to determine the evolution of a stochastic process $\vo{\xi}(t,\param)$, where $\param\in\mathcal{D}_{\param} \subset\real^d$ represents the parameter space with probability density function $\pdfp$. Differential equations with probabilistic parameters e.g. $$\vo{\dot{\xi}}(t,\param) = \vo{F}(t,\vo{\xi}(t),\param), \eqnlabel{autoDyn}$$ are examples of such stochastic processes that are amenable for analysis using polynomial chaos theory.  Assuming $\vo{\xi}(t,\param)$ to be a second-order process, it can be expanded, with $\set{L}_2$ convergence [@CameronMartin], as $$\vo{\xi}(t,\param) = \sum_{i=0}^\infty \vo{\xi}_i(t)\phi_i(\param),$$ where $\vo{\xi}_i(t)$ are time varying coefficients, and $\phi_i(\param)$ are known basis polynomials.  For exponential convergence, $\phi_i(\param)$ are chosen to be orthogonal with respect to the probability density function $p(\param)$, i.e.

$$\Exp{\basis{i}{\param}\basis{j}{\param}}:= \int_{\mathcal{D}_{\param}}{\basis{i}{\param}\basis{j}{\param} \pdfp
\,d\param}  = h_i^2\delta_{ij}, \eqnlabel{basisFcn}$$

where $\delta_{ij}$ is the Kronecker delta, $h_i:=\int_{\mathcal{D}_{\param}}{\phi_i^2\pdfp\,d\param}$.

For computational purposes, we truncate the expansion to a finite number of terms, i.e. the solution of is approximated by the polynomial chaos expansion as $$\begin{aligned}
    \vo{\xi}(t,\param) \approx \vo{\hat{\xi}}(t,\param) =  \sum_{i=0}^N \vo{\xi}_i(t)\basis{i}{\param}.  \elab{approxSoln}\end{aligned}$$ For a more compact representation of the ensuing expressions, we define $\mo{\Phi}(\param)$ to be $$\begin{aligned}
\mo{\Phi}(\param) &:= \begin{bmatrix}\basis{0}{\param}, & \cdots, & \basis{N}{\param}\end{bmatrix}^T, \text{ and } \\
\Phin{n} &:= \mo{\Phi}(\param)\otimes \I{n},\end{aligned}$$ where $\I{n}\in\real^{n\times n}$ is identity matrix.  We also define matrix $\vo{\Xi}\in\real^{n\times(N+1)}$, with polynomial chaos coefficients $\x_i$, as $$\begin{aligned}
\vo{\Xi} = \begin{bmatrix} \vo{\xi}_0, & \cdots, & \vo{\xi}_N \end{bmatrix}.  \elab{defX}\end{aligned}$$ Therefore, $\vo{\hat{\xi}}(t,\param)$ can be written as $$\begin{aligned}
\vo{\hat{\xi}}(t,\param) := \vo{\Xi}(t)\mo{\Phi}(\param) \eqnlabel{compactX}.\end{aligned}$$ Noting that $\vo{\hat{\xi}} \equiv \vec{\vo{\hat{\xi}}}$, we obtain an alternate form for , $$\begin{aligned}
\vo{\hat{\xi}} & \equiv  \vec{\vo{\hat{\xi}}} = \vec{\vo{{\Xi}}\mo{\Phi}(\param)} = \vec{\I{n}\vo{{\Xi}}\mo{\Phi}(\param)} \nonumber \\
 & = (\mo{\Phi}^T(\param)\otimes \I{n})\vec{\vo{{\xi}}}  =  \mo{\Phi}_n^T(\param)\vo{\xi}_\text{pc}, \eqnlabel{compactxpc}\end{aligned}$$ where $\vo{\xi}_\text{pc} := \vec{\vo{\Xi}}$, and $\vec{\cdot}$ is the vectorization operator. 

The unknown coefficients $\vo{\xi}_\text{pc}$ are determined using one of many methods including Galerkin projection[@Ghanem1; @ghanem:91], stochastic collocation [@Marzouk2009a; @NASA-CR-2003-212153], and least-square minimization [@walters2003towards; @hosder2007efficient].  In this work, we pursue the Galerkin projection approach to determine the coefficients $\vo{\xi}_\text{pc}$ by first defining error $\vo{e}(t,\param):=\vo{\xi}(t,\param) - \mo{\Phi}_n^T(\param)\vo{\xi}_\text{pc}(t)$, and determining $\vo{\xi}_\text{pc}(t)$ from the following algebraic equations $$\int_\domain{\param} \vo{e}(t,\param) \phi_i(\param) p(\param) d\param = 0,$$ for $i=0,\cdots,N$. If $\vo{\xi}(t,\param)$ is solution of a differential equation , then the error is defined in terms of the equation error, as shown in . 

In general, polynomial chaos does not scale well with state-space and parameter dimension.  The number of basis functions for a given order $r$ with $d$ independent random variables is $\frac{(d+r)!}{d!r!}$. With large number of parameters (increasing $d$), the number of basis functions, for a given order of approximation, will increase factorially and the computational cost will be prohibitive.  This limits how large both $d$ and $r$ can be.  Recent development in sparse polynomial chaos may scale better [@constantine2012sparse; @conrad2013adaptive].  However, usually we can get quite good performance with low order approximations [@bhattacharyarobust; @Fisher2009; @bhattacharya2012linear; @Dutta2010].  Unfortunately, the order of approximation, for which acceptable accuracy is achieved, has to be determined empirically. 

For $d>1$, the polyvariate basis functions are determined from tensor-products of univariate polynomials, with limit on the total order of the product using Pascal’s triangle, the univariate polynomials can be determined from different distributions.  For a given distribution, using polynomials that are orthogonal with respect to the distribution, is usually chosen for exponential convergence [@Xiu].  Poor scalability of polynomial chaos is due to the tensor product of the basis functions.  However, anisotropic tensor products [@babuska2004galerkin; @babuvska2007stochastic] or anisotropic Smolyak cubature methods result in improved scaling [@nobile2008sparse]. 

In this paper, we consider elements of $\param$ to be independent.  However, in several applications this assumption may not valid.  For such applications, suitable transformation such as Rosenblatt [@rosenblatt1952remarks], Nataf [@der1986structural] and Box-Cox [@box1964analysis] transformation can be applied to arrive at a set of independent parameters.  An overview of such techniques is described in the work by Elred [@eldred2009comparison]. 

Robust Kalman filter {#sec:Robust Kalman filter}
====================

In Kalman filtering, state estimation involves two steps: a) model-based uncertainty propagation to obtain the *prior* state uncertainty, and b) incorporation of measurements to update the prior to *posterior* state uncertainty by minimizing the error variance.  With probabilistic uncertainty in the system parameters, along with process noise, the propagation step becomes complicated.  In this paper, we solve this by computing the mean and variance of the propagated states using conditional expectations. 

The new robust Kalman filtering algorithms, for uncertain DT and CT systems, are presented next. 

Discrete Robust Kalman Filter {#sec:DRF}
-----------------------------

Let us consider the DT model given by and as

$$\begin{aligned}
\x_{k} &= \Ap \x_{k-1} + \Bp \w_{k-1}, \elab{sys}  \\
\y_k &= \C\x_{k} + \n_k .  \elab{y_dt}\end{aligned}$$

### **Uncertainty propagation** {#sec:ddkf_prop}

The uncertainty in $\x_k(\param,\w)$, the solution of , is due to uncertainty in the initial condition $\x_0$, the uncertainty in the system parameters $\param$, and the process noise $\w_k$. We will assume $\x_0$ and $\param$ are random variables with given probability density functions.  For simplicity, we assume $\x_0$, $\param$, and $\w_k$ are statistically independent.  It is noteworthy, that due to the uncertainty in the system matrices, the state probability density function will not be Gaussian, even if $\x_0$ is Gaussian.  However, we restrict ourselves to characterizing the first two moments of $\x_k(\param,\w)$ as defined below, since in this paper we are focusing on Kalman filtering.  Let us define

$$\begin{aligned}
\mub_k &:= \Exp{\x_k(\param,\w)}, \text{ and }\\
\sigb_k &:= \Exp{\left(\x_k(\param,\w)-\mub_k\right)\left(\x_k(\param,\w)-\mub_k\right)^T}.\end{aligned}$$

Consequently, the propagation equation for $\mub_k$ is given by $\mub_k = \E{\A(\param)\x(\param)_{k-1}} + \E{\B(\param)\w_{k-1}}$. We use the conditional expectation with respect to $\param$ to calculate the quantities in the previous equation.  For a given $\param$, the propagation equations are similar to those in standard Kalman filtering.  Assuming that the posteriors $\mub^{+}_{{k-1}}$ and $\sigb^{+}_{{k-1}}$ have *no uncertainty*, we can write the propagation equation for conditional mean and variance as

$$\begin{aligned}
\mub^{-}_{k}(\param) &= \Ap\mub^{+}_{{k-1}}, \elab{cond:mu}\\
\sigb^{-}_{k}(\param) &= \Ap\sigb^{+}_{{k-1}}\Apt + \Bp \Q \Bpt, \elab{cond:var}\end{aligned}$$

where $\mub^{+}_{k}(\param)$ and $\sigb^{+}_{k}(\param)$ are *stochastic* since they depend on $\param$. The *total* mean and variance of $\x_k(\param)$ can be computed from the conditional mean and variance as

$$\begin{aligned}
\mub^{-}_{k} &:= \E{\mub^{-}_{k}(\param)}, \\
\sigb^{-}_{k}  &:= \E{\sigb^{-}_{k}(\param)} + \textbf{Var}\left(\mub^{-}_{k}(\param)\right).   \eqnlabel{totalStatVar}
%  \Exp{\left(\mub^{-}_{k}(\param)-\mub^{-}_{k}\right)\left(\mub^{-}_{k}(\param)-\mub^{-}_{k}\right)^T}.\end{aligned}$$

With slight abuse of notation, we represent the conditional mean and variance as functions of $\param$, i.e. $\mub^{-}_{k}(\param)$ and $\sigb^{-}_{k}(\param)$. Whereas the total mean and variance are represented without the functional dependence, i.e. $\mub^{-}_{k}$ and $\sigb^{-}_{k}$.

Since the posterior $\mub^{+}_{{k-1}}$ is independent of $\param$, the total prior mean is calculated as $$\begin{aligned}
\mub^{-}_{k} &:= \E{\mub^{-}_{k}(\param)} = \E{\Ap\mub^{+}_{k-1}} =\Bar{\A} \mub^{+}_{{k-1}}, \eqnlabel{mean_prior} %= \E{\Ap} \mub^{+}_{{k-1}} \nonumber\\ &=\Bar{\A} \mub^{+}_{{k-1}},\end{aligned}$$ where $\Bar{\A}:=\E{\Ap}$. The variance of conditional mean, $\textbf{Var}\left(\mub^{-}_{k}(\param)\right)$, can be determined as $$\begin{aligned}
&\textbf{Var}\left(\mub^{-}_{k}(\param)\right) = \Exp{\left(\mub^{-}_{k}(\param)-\mub^{-}_{k}\right)\left(\mub^{-}_{k}(\param)-\mub^{-}_{k}\right)^T}\nonumber\\
&= \Exp{\left(\Ap-\bar{\A}\right) \left(\mub^{+}_{k-1}  \mub^{+}_{k-1}{}^T \right) \left(\Ap-\bar{\A} \right)^T}.  \eqnlabel{var_mu}\end{aligned}$$ Therefore, the total prior variance follows from , , and as $$\begin{aligned}
& \sigb^{-}_{k}  := \E{\sigb^{-}_{k}(\param)} + \textbf{Var}\left(\mub^{-}_{k}(\param)\right) \nonumber \\
&=\E{ \Ap\sigb^{+}_{k-1}\Apt} + \E{\Bp \Q \Bpt }\nonumber \\
&+ \Exp{\left(\Ap-\bar{\A}\right) \left(\mub^{+}_{{k-1}}  \mub^{+}_{k-1}{}^T \right) \left(\Ap-\bar{\A}\right)^T}.  \eqnlabel{var_prior}
% &+ \Exp{\Ap\mub^{+}_{{k-1}} ( \mub^{+}_{{k-1}})^T \Ap^T} - \bar{\A} \mub^{+}_{{k-1}}( \mub^{+}_{{k-1}})^T \bar{\A}^T\end{aligned}$$

### **Update** {#sec:update}

Since we have assumed the matrix $\C$ in the measurement model to be independent of $\param$, we can simply follow the standard Kalman update equations.  For the brevity of discussion, we omit the step by step derivation of the well known Kalman gain and update equations, which can be found in many textbooks, e.g. [@bryson2018applied].  Once we have the propagated priors from equations and , the posteriors are given by

$$\begin{aligned}
 \mub^{+}_{k} &= \mub^{-}_{k} + \K_k\left(\y_k - \C\mub^{-}_{k}\right), \eqnlabel{mean:update1} \\
\sigb^{+}_{k} &= \left(\vo{I}-\K_k \C \right)\sigb^{-}_{k} + \K_k\R\K_k^T, \eqnlabel{var:update1}\end{aligned}$$

where $\y_k$ is is the sensor measurement, and $\K_k :=\sigb^{-}_{k}\C^T$ $[\C\sigb^{-}_{k}\C^T+\R]^{-1}$ is the optimal Kalman gain to achieve an unbiased estimate with the minimum state error variance. 

Continuous-Discrete Robust Kalman Filter {#sec:CRF}
----------------------------------------

The continuous-discrete filter, also known as the hybrid Kalman filter, is more practical than other filters as it is suitable for most physical dynamical systems that are governed by continuous time ODEs, and sensor measurements are available only at discrete time instants.  The system and sensor equations follow from and as

$$\begin{aligned}
 \xdot(t) &= \Ap \x(t) + \Bp\w(t), \elab{sys_ct} \\
\y_k &= \C\x_{k} + \n_k .\elab{y_ct}\end{aligned}$$

### **Uncertainty Propagation** {#sec:UP}

Determining the moments of $\x(t,\param,\w)$ is nontrivial in this case, particularly due to $\param$. This can be shown by first defining mean and covariance as

$$\begin{aligned}
\mub(t) &:= \Exp{\x(t,\param,\w)}, \text{ and }\\
\sigb(t) &:= \Exp{\left(\x(t,\param,\w)-\mub(t)\right)\left(\x(t,\param,\w)-\mub(t)\right)^T}.\end{aligned}$$

The propagation equation for $\mub(t)$ is given by $$\dot{\mub}(t) = \E{\A(\param)\x(\param,t)} + \E{\B(\param)\w(t)},$$ which presents a challenge in solving the differential equation due to uncertain matrices $\Ap$ and $\Bp$. Similar difficulty is faced in the propagation equation for $\sigb(t)$. We next present an approach based on the polynomial chaos theory to determine the first two moments of $\x(t,\param,\w)$.

As in the previous section, we adopt the formulation based on the conditional expectation with respect to $\param$. For a given $\param$, we can write the propagation equation for conditional mean and variance as

$$\begin{aligned}
\dot{\mub}(t,\param) &= \Ap\mub(t,\param), \eqnlabel{cond:muCT}\\ %\elab{cond:mu}\\
\dot{\sigb}(t,\param) &= \Ap\sigb(t,\param)  + \sigb(t,\param)\Apt \nonumber \\ & \quad \quad + \Bp\Q\Bpt, \eqnlabel{cond:varCT}\end{aligned}$$

Equations and are integrated using $\mub_0$ and $\sigb_0$ as the initial condition, which are *deterministic* and define the uncertainty in the initial condition of the state.  It is assumed to be independent of $\param$.

The total mean and variance of $\x(t,\param)$ can be computed as

$$\begin{aligned}
\mub(t) &:= \E{\mub(t,\param)},\\
\sigb(t) &:= \E{\sigb(t,\param)}  + \textbf{Var}(\mub(t,\param)). %  \Exp{\left(\mub(t,\param)-\mub(t)\right)\left(\mub(t,\param)-\mub(t)\right)^T}.\end{aligned}$$

Stochastic processes $\mub(t,\param)$ and $\sigb(t,\param)$ are expanded with polynomial chaos basis functions as follows.\
The expansion for $\mub(t,\param)$ follows from as $$\begin{aligned}
\hat{\mub}(t,\param) &= \sum_{i=0}^{N} \mub_i(t)\phi_i(\param)
=\begin{bmatrix} \mub_0(t) & \cdots & \mub_N(t)\end{bmatrix}\mo{\Phi}(\param) \nonumber\\
&=  \vo{{\tilde{\mu}_{\text{pc}}}}\Phin{n} = \Phint{n}\mub_{\text{pc}} , \eqnlabel{mu_pc}
%\vo{\widetilde{\mu_{\text{pc}}}}\Phin{n} \\
% &= (\mo{\Phi}^{T}(\param)\otimes \I{n}) \vo{\mub}_{\text{pc}}\end{aligned}$$ where, $$\begin{aligned}
 \vo{{\tilde{\mu}_{\text{pc}}}} &:= \begin{bmatrix} \mub_0(t) & \cdots & \mub_N(t)\end{bmatrix} \in \real^{n(N+1)}, \nonumber\\
 \mub_{\text{pc}} &:= \vec{\vo{{\tilde{\mu}_{\text{pc}}}}} \in \real^{n(N+1)}.  \nonumber\end{aligned}$$

Since $\sigb(t,\param)\geq \vo{0}$, the stochastic process $\sigb(t,\param)$ is expanded using quadratic basis functions constructed from $\phi_i$. We adopt the expansion presented in [@bhattacharyarobust], i.e. $$\begin{aligned}
\hat{\sigb}(t,\param) &= \Phint{n} \begin{bmatrix} \sigb_{00} & \cdots & \sigb_{0N} \\ \vdots &  & \vdots \\ \sigb_{N0}  & \cdots & \sigb_{NN}
\end{bmatrix} \Phin{n}.  \nonumber\end{aligned}$$ Since $\hat{\sigb}(t,\param)$ is symmetric and $\hat{\sigb}(t,\param)\geq0$ , it follows that $\sigb_{ij} = \sigb^T_{ij} = \sigb_{ji} \geq 0$. Therefore, $\hat{\sigb}(t,\param)$ can be expanded as $$\begin{aligned}
\hat{\sigb}(t,\param) &=  \sum_{ij} \sigb_{ij}(t)\phi_i(\param)\phi_j(\param).  \nonumber\end{aligned}$$ Moreover, we note that the quadratic basis functions, $\{\phi_i(\param)\phi_j(\param)\}$, are not linearly independent.  Therefore, the PC expansion for $\hat{\sigb}(t,\param)$ can be effectively written as $$\begin{aligned}
\hat{\sigb}(t,\param) &=  \sum_{i=0}^M \sigb_{i}(t)\theta_i(\param) = \big[\sigb_{0}(t),\cdots, \sigb_{M}(t) \big] \Thp  \nonumber \\ &=  (\Thpt \otimes \vo{I}_n)\sigb_{\text{pc}} =: \Thnpt{n}\sigb_{\text{pc}}, \eqnlabel{var_pc}\end{aligned}$$ where, $M:=2(N-1)$, $0\leq\sigb_{i}(t) \in \real^{n\times n}$ and $$\begin{aligned}
\sigb_{\text{pc}} &:=
\begin{bmatrix} \sigb_0(t) & \cdots & \sigb_M(t)\end{bmatrix}^T \in \real^{n(M+1)\times n}.\nonumber\end{aligned}$$ The basis functions $\theta_i(\param)$ are linearly independent polynomials chosen from quadratic terms resulting from the expansion of $\big(\phi_0(\param)+\phi_1(\param)+\cdots+\phi_N(\param)\big)^2$, i.e. $\theta_i(\param)$ are linearly independent basis functions selected from the following set $$\begin{aligned}
\begin{Bmatrix} \phi_0(\param)\phi_0(\param) \\ 2\phi_0(\param)\phi_1(\param) \\ \vdots\\ 2\phi_{N-1}(\param)\phi_{N}(\param) \\ \phi_{N}(\param)\phi_{N}(\param)\end{Bmatrix},\end{aligned}$$ and, $$\begin{aligned}
\Thp :=
\left[ \theta_0(\param), \theta_1(\param), \cdots, \theta_M(\param)\right]^T \in \real^{M+1}.  \nonumber \end{aligned}$$ With this mean and variance approximation, the error equations in and are

$$\begin{aligned}
&\vo{e}_{\mub}(t,\param) := \Phint{n}\dot{\mub}_{\text{pc}} - \Ap\Phint{n}\mub_{\text{pc}}, \text{ and }\\
&\vo{e}_{\sigb}(t,\param) := \Thnpt{n}\dot{\sigb}_{\text{pc}} - \Ap\Thnpt{n}\sigb_{\text{pc}} \nonumber \\
 & \quad \quad \quad -\Thnpt{n}\sigb_{\text{pc}} \Apt - \Bp\Q\Bpt.\end{aligned}$$

The differential equations for $\dot{\mub}_i(t)$ and $\dot{\sigb}_i(t)$ are obtained by setting $$\begin{aligned}
\E{\vo{e}_{\mub}(t,\param)\phi_i(\param)} &= 0, \;\text{for } i=0,\cdots,N;\\
\E{\vo{e}_{\sigb}(t,\param)\theta_j(\param)} &= 0,\; \text{for } j=0,\cdots,M;
\end{aligned}$$ resulting in $$\begin{aligned}
 \dot{\mub}_{\text{pc}}  = \A_{\mub}  \mub_{\text{pc}}\, , \, \text{ and } \,
 \dot{\sigb}_{\text{pc}} = \F_{\sigb}  + \B_{\sigb}\, ,
\eqnlabel{propagation}\end{aligned}$$ where, $$\begin{aligned}
&\A_{\mub} :=\E{\Phin{n}\Phint{n}}^{-1}\E{\Phin{n}\Ap\Phint{n}},  \\
&\F_{\sigb} := \E{\Thnp{n}\Thnpt{n}}^{-1}  \times \\ & \E{\Thnp{n}\Ap\Thnpt{n}\sigb_{\text{pc}}
+ \Thnp{n}\Thnpt{n}\sigb_{\text{pc}}\Apt},  \\
&\B_{\sigb} := \E{\Thnp{n}\Thnpt{n}} ^{-1} \E{\Thnp{n}\Bp\Q\Bpt}.\end{aligned}$$

Given the posteriors $\mub^+(t_{k-1})$ and $\sigb^+(t_{k-1})$, at time instant $t_{k-1}$, the evolution of the state uncertainty is determined by integrating over $[t_{k-1},t_k]$ to arrive at $\mub^-(t_k,\param)$ and $\sigb^-(t_k,\param)$, the conditional prior mean and the conditional prior variance of the state.  The total mean and covariance priors, i.e. $\mub^-(t_k)$ and $\sigb^-(t_k)$, are then determined from . 

Integration of requires initial conditions $\mub_\text{pc}^+(t_{k-1})$ and $\sigb_\text{pc}^+(t_{k-1})$, which are determined by projecting $\mub^+(t_{k-1})$ and $\sigb^+(t_{k-1})$ on the basis functions $\{\phi_i(\param)\}_{i=0}^{N}$, and $\{\theta_i(\param)\}_{i=0}^{M}$ respectively.  Noting that $\mub^+(t_{k-1})$ and $\sigb^+(t_{k-1})$ are $\param$ independent, initial conditions $\mub_\text{pc}(t_{k-1})$ and $\sigb_\text{pc}(t_{k-1})$ are given by $$\begin{aligned}
\mub_\text{pc}^+(t_{k-1}) := \begin{bmatrix}
\mub^+(t_{k-1}) \\
\vo{0}_{nN}
\end{bmatrix},
\sigb_\text{pc}^+(t_{k-1}) := \begin{bmatrix}
\sigb^+(t_{k-1}) \\
\vo{0}_{nM\times n}
\end{bmatrix}.  \nonumber\end{aligned}$$ With these initial conditions, linear ODEs can be integrated to calculate $\mub_\text{pc}^-(t_{k})$ and $\sigb_\text{pc}^-(t_{k})$ at time $t_k$. Therefore, conditional mean and covariance priors at $t_k$ follow from and as $$\begin{aligned}
\mub^-(t_k,\param) &= \Phint{n}\mub_\text{pc}^-(t_{k}), \\ \sigb^-(t_k,\param) &= \Thnpt{n}\sigb_\text{pc}^-(t_{k}).\end{aligned}$$ The total mean and covariance priors $\mub^-(t_k)$ and $\sigb^-(t_k)$ are calculated using as follows.  $$\begin{aligned}
\mub^-(t_k) &= \E{\Phint{n}} \mub_\text{pc}^-(t_{k}), \nonumber \\
\sigb^-(t_k) &= \E{\Thnpt{n}}\sigb_\text{pc}^-(t_{k}) \nonumber  \\
&+ \vo{\tilde{\mu}^-_{\text{pc}}}
(t_{k}) \Big(\textbf{Var}\left(\mo{\Phi}(\param)\right) \Big) \big( \vo{\tilde{\mu}^-_{\text{pc}}}(t_{k})\big)^T, \nonumber\end{aligned}$$ where, $\textbf{Var}\left(\mo{\Phi}(\param)\right):= \mathbb{E}\Big[\big(\mo{\Phi}(\param) - \overline{\vo{\Phi}} \big) \big(\mo{\Phi}(\param)- \overline{\vo{\Phi}} \big)^T  \Big]$, and $\overline{\vo{\Phi}}:=\E{\mo{\Phi}(\param)}$.

### **Update** {#update}

Since the measurements are obtained at discrete time instants, we can use the Kalman update equations from §\[sec:update\].  The updated posteriors are given by

$$\begin{aligned}
 &\mub^{+}(t_k) = \mub^{-}(t_k) + \K_k\left(\y(t_k) - \C\mub^{-}(t_k)\right), \nonumber \\
&\sigb^{+}(t_k) = \big(\vo{I}-\K(t_k) \C \big)\sigb^{-}(t_k) + \K(t_k)\R\K^T(t_k), \nonumber\end{aligned}$$

where, $\y(t_k)$ is the sensor measurement, and $$\K(t_k):=\sigb^{-}(t_k)\C^T[\C\sigb^{-}(t_k)\C^T+\R]^{-1}.$$

Result {#sec:result}
======

Performance of the proposed robust Kalman filter is tested with two cases of simulation: 1) Case I: Initial mean, $\mu_0 = [0 \ 0]^T $, for checking steady state error, 2) Case II: Initial mean, $\mu_0 \neq [0 \ 0]^T $, for checking convergence rate with initial uncertainty.  We compare the performance of the filter in terms of the estimation accuracy characterized by the mean and standard deviation (SD) of root mean squared (RMS) error, and the rate of convergence. 

Discrete Robust Kalman filter {#discrete-robust-kalman-filter}
-----------------------------

The proposed discrete robust Kalman filter discussed in §\[sec:DRF\] is applied to the example (\[eqn:example1\]) that was previously considered as a test problem in [@xie1994robustDT; @zhu2002design]. 

$$\begin{aligned}
\x_{k} &= \begin{bmatrix}0 & -0.5 \\ 1 &1 + \delta\end{bmatrix} \x_{k-1} +\begin{bmatrix} -6 \\ 1 \end{bmatrix}\w_{k-1}, \\
\y_k &= \begin{bmatrix} -100 & 10\end{bmatrix}\x_{k} + \n_k .\end{aligned}$$

\[eqn:example1\]

where $\delta$ is a uniformly distributed random parameter in $[-0.3 \ 0.3]$, and the variance of process and measurement noise is assumed to be unity, i.e. $Q=1$, $R=1$.

We choose uniformly spaced $10$ points in $[-0.3 \ 0.3]$ as samples for $\delta$. Then, mean and standard deviation of RMS error obtained for different realizations of the plant corresponding to different values of $\delta$, are considered as metrics for the estimation accuracy.  We compare the performance of the proposed filter with standard Kalman filter with nominal plant realization corresponding to $\delta = 0$. As claimed by the authors of [@zhu2002design], and verified by us, the filter presented in [@zhu2002design] performs better than the one discussed in [@xie1994robustDT].  Therefore, herein, we compare the performance of the proposed filter only with [@zhu2002design].  The simulation results for the proposed discrete robust Kalman filter, the nominal Kalman filter, and the filter from [@zhu2002design], are shown in Fig.  \[fig:DT\_2\] and TABLE \[Table:RMS\_DT\].  In both simulation Cases I and II, the proposed robust Kalman filter has the least RMS error than the other filters, as shown in TABLE \[Table:RMS\_DT\]. 

Moreover, for Case II as shown in Fig.  \[fig:DT\_2\], the proposed filter converges faster than the nominal KF, and its convergence rate is comparable to the filter from [@zhu2002design].  We also note that the computational time required for the proposed filter is comparable to that of nominal KF and the filter from [@zhu2002design]. 

\[Table:RMS\_DT\]

  -- ------------- ------------------------- ------------------------- --
        Case I              Case II                                    
      Mean  /  SD         Mean  /  SD                                  
         $x_1$          2.7325 / 1.7518           3.4675 / 2.0216      
         $x_2$          4.3049 / 2.9640         6.0837 / **3.5487**    
         $x_1$          0.4438 / 0.4136           2.4085 / 1.0595      
         $x_2$          4.4418 / 4.1355          24.0821 / 10.5982     
         $x_1$      **0.3182** / **0.2914**   **0.5666** / **0.4314**  
         $x_2$      **3.1846** / **2.9114**     **5.6669** / 4.3099    
  -- ------------- ------------------------- ------------------------- --

  : Comparison of RMS error in Discrete time filters. 

Continuous-Discrete Robust Kalman Filter {#continuous-discrete-robust-kalman-filter}
----------------------------------------

The proposed hybrid robust Kalman filter in §\[sec:CRF\] is applied to the example (\[eqn:example2\]) and its performance is compared with the nominal Kalman filter. 

$$\begin{aligned}
 \xdot(t) &= \begin{bmatrix} 0 & -1 + \delta \\ 1 & -0.5 \end{bmatrix} \x(t) +\begin{bmatrix} -2 \\ 1 \end{bmatrix}\w(t), \\
\y_k &= \begin{bmatrix} -100 & -100\end{bmatrix}\x(t_k) + \n(t_k) .\end{aligned}$$

\[eqn:example2\]

where $\delta$ is uniformly distributed in the interval $[-0.95 \ 0.95]$, and the variances of process and measurement noise are $Q=1$ and $R=1$. We use the similar performance metrics discussed in the previous subsection. 

The proposed robust Kalman filter is 2 times more accurate than the nominal KF in steady state as shown in TABLE \[Table:RMS\_Hybrid\].  Moreover, it shows faster convergence than the nominal KF as shown in Fig.\[fig:CT2\].  Again, we note that the computational time required for the proposed filter is comparable to the nominal KF. 

  -- ------------- ------------------------- ------------------------- --
        Case I              Case II                                    
      Mean  /  SD         Mean  /  SD                                  
         $x_1$          0.0223 / 0.0204           0.2052 / 0.1694      
         $x_2$          0.0195 / 0.0211           0.2038 / 0.1695      
         $x_1$      **0.0155** / **0.0092**   **0.1833** / **0.0783**  
         $x_2$      **0.0137** / **0.0077**   **0.1822** / **0.0782**  
  -- ------------- ------------------------- ------------------------- --

  : Comparison of RMS error in hybrid filters.[]{data-label="Table:RMS_Hybrid"}

Conclusion {#sec:Conclusion}
==========

In this paper, we proposed robust Kalman filter with probabilistic uncertainty in system parameters.  Mean and variance of the uncertain system are propagated using conditional probability and the polynomial chaos (PC) expansion framework.  The proposed approach demonstrates better performance in terms of estimation accuracy and convergence rate as compared with existing methods, which makes it very attractive for state estimation of uncertain systems. 

[^1]: $^{1}$Sunsoo Kim is a Ph. D student in the Department of Electrical and Computer Engineering, Texas A&M University, College Station, TX 77840, USA.  Email: [kimsunsoo@tamu.edu]{}

[^2]: $^{2}$Vedang M. Deshpande is a Ph. D student in the Department of Aerospace Engineering, Texas A&M University, College Station, TX 77843, USA.  Email: [vedang.deshpande@tamu.edu]{}

[^3]: $^{3}$Raktim Bhattacharya is with the Faculty of Aerospace Engineering, Texas A&M University, College Station, TX 77840, USA.  Email: [raktim@tamu.edu]{}

[^4]: This work was supported by NSF Grant \#1762825.