naginterfaces.library.tsa.kalman_​unscented_​state_​revcom

naginterfaces.library.tsa.kalman_unscented_state_revcom(irevcm, y, ly, x, st, xt, fxt, comm, lx=None, ropt=None)

kalman_unscented_state_revcom applies the Unscented Kalman Filter to a nonlinear state space model, with additive noise.

kalman_unscented_state_revcom uses reverse communication for evaluating the nonlinear functionals of the state space model.

For full information please refer to the NAG Library document for g13ej

https://support.nag.com/numeric/nl/nagdoc_30/flhtml/g13/g13ejf.html

Parameters

irevcmint

On initial entry: must be set to $0$ or $3$.

If $irevcm=0$, it is assumed that $t=0$, otherwise it is assumed that $t\ne 0$ and that kalman_unscented_state_revcom has been called at least once before at an earlier time step.

On intermediate entry: $irevcm$ must remain unchanged.

yfloat, array-like, shape $\left(\text{my}\right)$

$y_t$, the observed data at the current time point.

lyfloat, array-like, shape $(\text{my},\text{my})$

$L_y$, such that $L_y{L}_y^T={\Sigma }_y$, i.e., the lower triangular part of a Cholesky decomposition of the observation noise covariance structure. Only the lower triangular part of $ly$ is referenced.

If ${\Sigma }_y$ is time dependent, the value supplied should be for time $t$.

xfloat, ndarray, shape $\left(\text{mx}\right)$, modified in place

On initial entry: ${\hat{x}}_{t-1}$ the state vector for the previous time point.

On intermediate exit: when

$irevcm=1$

$x$ is unchanged.

$irevcm=2$

${\hat{x}}_{\overline{t}}$.

On intermediate entry: $x$ must remain unchanged.

On final exit: ${\hat{x}}_t$ the updated state vector.

stfloat, ndarray, shape $(\text{mx},\text{mx})$, modified in place

On initial entry: $S_t$, such that $S_{t-1}{S}_{t-1}^T=P_{t-1}$, i.e., the lower triangular part of a Cholesky decomposition of the state covariance matrix at the previous time point. Only the lower triangular part of $st$ is referenced.

On intermediate exit: when

$irevcm=1$

$st$ is unchanged.

$irevcm=2$

$S_{\overline{t}}$, the lower triangular part of a Cholesky factorization of $P_{\overline{t}}$.

On intermediate entry: $st$ must remain unchanged.

On final exit: $S_t$, the lower triangular part of a Cholesky factorization of the updated state covariance matrix.

xtfloat, ndarray, shape $(\text{mx},:)$, modified in place

Note: the required extent for this argument in dimension 2 is determined as follows: if $irevcm=0$: $max(\text{my},n)$; if $irevcm\text{in}(1,2,3)\text{and}\text{licomm}\ge 10$: $max(\text{my},comm{[}^{\prime }icom{m}^{\prime }]\left[9\right])$; otherwise: $0$.

On initial entry: need not be set.

On intermediate exit: $X_t$ when $irevcm=1$, otherwise $Y_t$.

For the $j$th sigma point, the value for the $i$th parameter is held in $xt[\text{i}-1,\text{j}-1]$, for $\text{j}=1,2,\dots ,n$, for $\text{i}=1,2,\dots ,\text{mx}$.

On intermediate entry: $xt$ must remain unchanged.

On final exit: the contents of $xt$ are undefined.

fxtfloat, ndarray, shape $(:,:)$, modified in place

Note: the required extent for this argument in dimension 1 is determined as follows: if $irevcm\text{in}(1,2)$: $max(\text{mx},\text{my})$; otherwise: $0$.

Note: the required extent for this argument in dimension 2 is determined as follows: if $irevcm=1$: $n+\text{mx}$; if $irevcm=2$: $n+\text{my}$; otherwise: $0$.

On initial entry: need not be set.

On intermediate exit: the contents of $fxt$ are undefined.

On intermediate entry: $F\left(X_t\right)$ when $irevcm=1$, otherwise $H\left(Y_t\right)$ for the values of $X_t$ and $Y_t$ held in $xt$.

For the $j$th sigma point the value for the $i$th parameter should be held in $fxt[i-1,\text{j}-1]$, for $\text{j}=1,2,\dots ,n$.

When $irevcm=1$, $i=1,2,\dots ,\text{mx}$ and when $irevcm=2$, $i=1,2,\dots ,\text{my}$.

On final exit: the contents of $fxt$ are undefined.

commdict, communication object, modified in place

Communication structure.

On initial entry: need not be set.

lxNone or float, array-like, shape $(:,\text{mx})$, optional

$L_x$, such that $L_x{L}_x^T={\Sigma }_x$, i.e., the lower triangular part of a Cholesky decomposition of the process noise covariance structure. Only the lower triangular part of $lx$ is referenced.

If $\text{ldlx}=0$, there is no process noise ($v_t=0$ for all $t$) and $lx$ is not referenced.

If ${\Sigma }_x$ is time dependent, the value supplied should be for time $t$.

roptNone or float, array-like, shape $\left(\text{lropt}\right)$, optional

Options. The default value will be used for $ropt[i-1]$ if $\text{lropt}<i$. Setting $\text{lropt}=0$ will use the default values for all options and $ropt$ need not be set.

$ropt\left[0\right]$

If set to $1$ then the second set of sigma points are redrawn, as given by equation [equation]. If set to $2$ then the second set of sigma points are generated via augmentation, as given by equation [equation].

Default is for the sigma points to be redrawn (i.e., $ropt\left[0\right]=1$)

$ropt\left[1\right]$

${\kappa }_x$, value of $\kappa$ used when constructing the first set of sigma points, $X_t$.

Defaults to $3-\text{mx}$.

$ropt\left[2\right]$

${\alpha }_x$, value of $\alpha$ used when constructing the first set of sigma points, $X_t$.

Defaults to $1$.

$ropt\left[3\right]$

${\beta }_x$, value of $\beta$ used when constructing the first set of sigma points, $X_t$.

Defaults to $2$.

$ropt\left[4\right]$

Value of $\kappa$ used when constructing the second set of sigma points, $Y_t$.

Defaults to $3-2\times \text{mx}$ when $\text{ldlx}\ne 0$ and the second set of sigma points are augmented and ${\kappa }_x$ otherwise.

$ropt\left[5\right]$

Value of $\alpha$ used when constructing the second set of sigma points, $Y_t$.

Defaults to ${\alpha }_x$.

$ropt\left[6\right]$

Value of $\beta$ used when constructing the second set of sigma points, $Y_t$.

Defaults to ${\beta }_x$.

Returns

irevcmint

On intermediate exit: $irevcm=1$ or $2$. The value of $irevcm$ specifies what intermediate values are returned by this function and what values the calling program must assign to arguments of kalman_unscented_state_revcom before re-entering the routine. Details of the output and required input are given in the individual argument descriptions.

On final exit: $irevcm=3$

Raises

NagValueError

(errno $11$)

On entry, $irevcm=⟨value⟩$.

Constraint: $irevcm=0$, $1$, $2$ or $3$.

(errno $21$)

On entry, $\text{mx}=⟨value⟩$.

Constraint: $\text{mx}\ge 1$.

(errno $22$)

$\text{mx}$ has changed between calls.

On intermediate entry, $\text{mx}=⟨value⟩$.

On initial entry, $\text{mx}=⟨value⟩$.

(errno $31$)

On entry, $\text{my}=⟨value⟩$.

Constraint: $\text{my}\ge 1$.

(errno $32$)

$\text{my}$ has changed between calls.

On intermediate entry, $\text{my}=⟨value⟩$.

On initial entry, $\text{my}=⟨value⟩$.

(errno $121$)

On entry, augmented sigma points requested, $n=⟨value⟩$ and $\text{mx}=⟨value⟩$.

Constraint: $n\ge ⟨value⟩$.

(errno $122$)

On entry, redrawn sigma points requested, $n=⟨value⟩$ and $\text{mx}=⟨value⟩$.

Constraint: $n\ge ⟨value⟩$.

(errno $123$)

$\text{n}$ has changed between calls.

On intermediate entry, $n=⟨value⟩$.

On intermediate exit, $n=⟨value⟩$.

(errno $171$)

On entry, $ropt\left[0\right]=⟨value⟩$.

Constraint: $ropt\left[0\right]=1$ or $2$.

(errno $172$)

On entry, $ropt[⟨value⟩]=⟨value⟩$.

Constraint: $\kappa >⟨value⟩$.

(errno $173$)

On entry, $ropt[⟨value⟩]=⟨value⟩$.

Constraint: $\alpha >0$.

(errno $181$)

On entry, $\text{lropt}=⟨value⟩$.

Constraint: $0\le \text{lropt}\le 7$.

(errno $191$)

$comm$[‘icomm’] has been corrupted between calls.

(errno $211$)

$comm$[‘rcomm’] has been corrupted between calls.

(errno $301$)

A weight was negative and it was not possible to downdate the Cholesky factorization.

(errno $302$)

Unable to calculate the Kalman gain matrix.

(errno $303$)

Unable to calculate the Cholesky factorization of the updated state covariance matrix.

Warns

NagAlgorithmicWarning

(errno $202$)

On entry, $\text{licomm}=⟨value⟩$ and $\text{lrcomm}=⟨value⟩$.

Constraint: $\text{licomm}\ge 30$ and $\text{lrcomm}\ge 30+\text{my}+\text{mx}\times \text{my}+2\times max(\text{mx},\text{my})$.

The minimum required values for $\text{licomm}$ and $\text{lrcomm}$ are returned in $comm{[}^{\prime }icom{m}^{\prime }\left]\right[0]$ and $comm{[}^{\prime }icom{m}^{\prime }\left]\right[1]$ respectively.

Notes

kalman_unscented_state_revcom applies the Unscented Kalman Filter (UKF), as described in Julier and Uhlmann (1997b) to a nonlinear state space model, with additive noise, which, at time $t$, can be described by:

$$\begin{array}{c}\begin{array}{cc}{x}_{t+1}& =F\left(x_t\right)+v_t\\ y_t& =H\left(x_t\right)+u_t\end{array}\end{array}$$

where $x_t$ represents the unobserved state vector of length $m_x$ and $y_t$ the observed measurement vector of length $m_y$. The process noise is denoted $v_t$, which is assumed to have mean zero and covariance structure ${\Sigma }_x$, and the measurement noise by $u_t$, which is assumed to have mean zero and covariance structure ${\Sigma }_y$.

Unscented Kalman Filter Algorithm

Given ${\hat{x}}_0$, an initial estimate of the state and $P_0$ and initial estimate of the state covariance matrix, the UKF can be described as follows:

1. Generate a set of sigma points (see Sigma Points):

   $$\begin{array}{cc}{X}_t& =[{\hat{x}}_{t-1}{\hat{x}}_{t-1}+\gamma \sqrt{P_{t-1}}{\hat{x}}_{t-1}-\gamma \sqrt{P_{t-1}}]\end{array}$$

2. Evaluate the known model function $F$:

   $$\begin{array}{cc}{F}_t& =F\left(X_t\right)\end{array}$$

   The function $F$ is assumed to accept the $m_x\times n$ matrix, $X_t$ and return an $m_x\times n$ matrix, $F_t$. The columns of both $X_t$ and $F_t$ correspond to different possible states. The notation $F_{t,i}$ is used to denote the $i$th column of $F_t$, hence the result of applying $F$ to the $i$th possible state.

3. Time Update:

   $$\begin{array}{c}\begin{array}{cc}{\hat{x}}_{\overline{t}}& ={\sum }_{i=1}^n{W}_i^m{F}_{t,i}\\ P_{\overline{t}}& ={\sum }_{i=1}^n{W}_i^c(F_{t,i}-{\hat{x}}_{\hat{t}}){(F_{t,i}-{\hat{x}}_{\hat{t}})}^T+{\Sigma }_x\end{array}\end{array}$$

4. Redraw another set of sigma points (see Sigma Points):

   $$\begin{array}{cc}{Y}_t& =[{\hat{x}}_{\hat{t}}{\hat{x}}_{\hat{t}}+\gamma \sqrt{P_{\overline{t}}}{\hat{x}}_{\hat{t}}-\gamma \sqrt{P_{\overline{t}}}]\end{array}$$

5. Evaluate the known model function $H$:

   $$\begin{array}{cc}{H}_t& =H\left(Y_t\right)\end{array}$$

   The function $H$ is assumed to accept the $m_x\times n$ matrix, $Y_t$ and return an $m_y\times n$ matrix, $H_t$. The columns of both $Y_t$ and $H_t$ correspond to different possible states. As above $H_{t,i}$ is used to denote the $i$th column of $H_t$.

6. Measurement Update:

   $$\begin{array}{c}\begin{array}{cc}{\hat{y}}_t& ={\sum }_{i=1}^n{W}_i^m{H}_{t,i}\\ P_{y{y}_t}& ={\sum }_{i=1}^n{W}_i^c(H_{t,i}-{\hat{y}}_t){(H_{t,i}-{\hat{y}}_t)}^T+{\Sigma }_y\\ P_{x{y}_t}& ={\sum }_{i=1}^n{W}_i^c(F_{t,i}-{\hat{x}}_{\hat{t}}){(H_{t,i}-{\hat{y}}_t)}^T\\ K_t& =P_{x{y}_t}{P}_{y{y}_t}^{-1}\\ {\hat{x}}_t& ={\hat{x}}_{\hat{t}}+K_t(y_t-{\hat{y}}_t)\\ P_t& =P_{\overline{t}}-K_t{P}_{y{y}_t}{K}_t^T\end{array}\end{array}$$

Here $K_t$ is the Kalman gain matrix, ${\hat{x}}_t$ is the estimated state vector at time $t$ and $P_t$ the corresponding covariance matrix. Rather than implementing the standard UKF as stated above kalman_unscented_state_revcom uses the square-root form described in the Haykin (2001).

Sigma Points

A nonlinear state space model involves propagating a vector of random variables through a nonlinear system and we are interested in what happens to the mean and covariance matrix of those variables. Rather than trying to directly propagate the mean and covariance matrix, the UKF uses a set of carefully chosen sample points, referred to as sigma points, and propagates these through the system of interest. An estimate of the propagated mean and covariance matrix is then obtained via the weighted sample mean and covariance matrix.

For a vector of $m$ random variables, $x$, with mean $\mu$ and covariance matrix $\Sigma$, the sigma points are usually constructed as:

$$X_t=[\mu \mu +\gamma \sqrt{\Sigma }\mu -\gamma \sqrt{\Sigma }]$$

When calculating the weighted sample mean and covariance matrix two sets of weights are required, one used when calculating the weighted sample mean, denoted $W^m$ and one used when calculating the weighted sample covariance matrix, denoted $W^c$. The weights and multiplier, $\gamma$, are constructed as follows:

$$\begin{array}{c}\begin{array}{cc}\lambda & ={\alpha }^2(L+\kappa )-L\\ \gamma & =\sqrt{L+\lambda }\\ W_i^m& =\{\begin{array}{cc}\frac{\lambda }{L+\lambda }& i=1\\ \frac{1}{2(L+\lambda )}& i=2,3,\dots ,2L+1\end{array}\\ W_i^c& =\{\begin{array}{cc}\frac{\lambda }{L+\lambda }+1-{\alpha }^2+\beta & i=1\\ \frac{1}{2(L+\lambda )}& i=2,3,\dots ,2L+1\end{array}\end{array}\end{array}$$

where, usually $L=m$ and $\alpha ,\beta$ and $\kappa$ are constants. The total number of sigma points, $n$, is given by $2L+1$. The constant $\alpha$ is usually set to somewhere in the range ${10}^{-4}\le \alpha \le 1$ and for a Gaussian distribution, the optimal values of $\kappa$ and $\beta$ are $3-L$ and $2$ respectively.

Rather than redrawing another set of sigma points in (d) of the UKF an alternative method can be used where the sigma points used in (a) are augmented to take into account the process noise. This involves replacing equation [equation] with:

$$\begin{array}{cc}{Y}_t& =[X_t{X}_{t,1}+\gamma \sqrt{{\Sigma }_x}{X}_{t,1}-\gamma \sqrt{{\Sigma }_x}]\end{array}$$

Augmenting the sigma points in this manner requires setting $L$ to $2L$ (and hence $n$ to $2n-1$) and recalculating the weights. These new values are then used for the rest of the algorithm. The advantage of augmenting the sigma points is that it keeps any odd-moments information captured by the original propagated sigma points, at the cost of using a larger number of points.

References

Haykin, S, 2001, Kalman Filtering and Neural Networks, John Wiley and Sons

Julier, S J, 2002, The scaled unscented transformation, Proceedings of the 2002 American Control Conference (Volume 6), 4555–4559

Julier, S J and Uhlmann, J K, 1997, A consistent, debiased method for converting between polar and Cartesian coordinate systems, Proceedings of AeroSense ‘97, International Society for Optics and Phonotonics, 110–121

Julier, S J and Uhlmann, J K, 1997, A new extension of the Kalman Filter to nonlinear systems, International Symposium for Aerospace/Defense, Sensing, Simulation and Controls (Volume 3) (26)