NAG CL Interface
g13ekc (kalman_​unscented_​state)

Settings help

CL Name Style:


1 Purpose

g13ekc applies the Unscented Kalman Filter (UKF) to a nonlinear state space model, with additive noise.
g13ekc uses direct communication for evaluating the nonlinear functionals of the state space model.

2 Specification

#include <nag.h>
void  g13ekc (Integer mx, Integer my, const double y[], const double lx[], const double ly[],
void (*f)(Integer mx, Integer n, const double xt[], double fxt[], Nag_Comm *comm, Integer *info),
void (*h)(Integer mx, Integer my, Integer n, const double yt[], double hyt[], Nag_Comm *comm, Integer *info),
double x[], double st[], Nag_Comm *comm, NagError *fail)
The function may be called by the names: g13ekc, nag_tsa_kalman_unscented_state or nag_kalman_unscented_state.

3 Description

g13ekc 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:
xt+1 =F(xt)+vt yt =H(xt)+ut  
where xt represents the unobserved state vector of length mx and yt the observed measurement vector of length my. The process noise is denoted vt, which is assumed to have mean zero and covariance structure Σx, and the measurement noise by ut, which is assumed to have mean zero and covariance structure Σy.

3.1 Unscented Kalman Filter Algorithm

Given x^0, an initial estimate of the state and P0 and initial estimate of the state covariance matrix, the UKF can be described as follows:
  1. (a)Generate a set of sigma points (see Section 3.2):
    Xt = [ x ^ t-1     x ^ t-1 +γPt-1    x ^ t-1 -γPt-1] (1)
  2. (b)Evaluate the known model function F:
    Ft =F(Xt) (2)
    The function F is assumed to accept the mx×n matrix, Xt and return an mx×n matrix, Ft. The columns of both Xt and Ft correspond to different possible states. The notation Ft,i is used to denote the ith column of Ft, hence the result of applying F to the ith possible state.
  3. (c)Time Update:
    x^t¯ = i=1 n Wim Ft,i (3)
    Pt¯ = i=1 n Wic (Ft,i- x^ t^ ) (Ft,i- x^ t^ ) T + Σx (4)
  4. (d)Redraw another set of sigma points (see Section 3.2):
    Yt = [ x^ t^     x^ t^ +γPt¯    x^ t^ -γPt¯] (5)
  5. (e)Evaluate the known model function H:
    Ht =H(Yt) (6)
    The function H is assumed to accept the mx×n matrix, Yt and return an my×n matrix, Ht. The columns of both Yt and Ht correspond to different possible states. As above Ht,i is used to denote the ith column of Ht.
  6. (f)Measurement Update:
    y ^ t = i=1 n Wim Ht,i (7)
    Pyyt = i=1 n Wic (Ht,i- y ^ t ) (Ht,i- y ^ t ) T + Σy (8)
    P xyt = i=1 n Wic (Ft,i- x^ t^ ) (Ht,i- y^ t ) T (9)
    Kt = P xyt Pyyt-1 (10)
    x^t = x^ t^ + Kt (yt- y^ t ) (11)
    Pt = Pt¯ - Kt Pyyt KtT (12)
Here Kt is the Kalman gain matrix, x^t is the estimated state vector at time t and Pt the corresponding covariance matrix. Rather than implementing the standard UKF as stated above g13ekc uses the square-root form described in the Haykin (2001).

3.2 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 μ and covariance matrix Σ, the sigma points are usually constructed as:
Xt = [μ   μ+γΣ   μ-γΣ]  
When calculating the weighted sample mean and covariance matrix two sets of weights are required, one used when calculating the weighted sample mean, denoted Wm and one used when calculating the weighted sample covariance matrix, denoted Wc. The weights and multiplier, γ, are constructed as follows:
λ =α2(L+κ)-L γ =L+λ Wim = { λL+λ      i=1 12(L+λ)      i=2,3,,2L+1 Wic = { λL+λ +1-α2+β i=1 12(L+λ) i=2,3,,2L+1  
where, usually L=m and α,β and κ are constants. The total number of sigma points, n, is given by 2L+1. The constant α is usually set to somewhere in the range 10−4α1 and for a Gaussian distribution, the optimal values of κ and β are 3-L and 2 respectively.
The constants, κ, α and β are given by κ=3-mx, α=1.0 and β=2. If more control is required over the construction of the sigma points then the reverse communication function, g13ejc, can be used instead.

4 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 (1997a) 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 (1997b) A new extension of the Kalman Filter to nonlinear systems International Symposium for Aerospace/Defense, Sensing, Simulation and Controls (Volume 3) 26

5 Arguments

1: mx Integer Input
On entry: mx, the number of state variables.
Constraint: mx1.
2: my Integer Input
On entry: my, the number of observed variables.
Constraint: my1.
3: y[my] const double Input
On entry: yt, the observed data at the current time point.
4: lx[mx×mx] const double Input
Note: the (i,j)th element of the matrix is stored in lx[(j-1)×mx+i-1].
On entry: Lx, such that LxLxT=Σx, i.e., the lower triangular part of a Cholesky decomposition of the process noise covariance structure. Only the lower triangular part of the matrix stored in lx is referenced.
If Σx is time dependent, the value supplied should be for time t.
5: ly[my×my] const double Input
Note: the (i,j)th element of the matrix is stored in ly[(j-1)×my+i-1].
On entry: Ly, such that LyLyT=Σy, i.e., the lower triangular part of a Cholesky decomposition of the observation noise covariance structure. Only the lower triangular part of the matrix stored in ly is referenced.
If Σy is time dependent, the value supplied should be for time t.
6: f function, supplied by the user External Function
The state function, F as described in (b).
The specification of f is:
void  f (Integer mx, Integer n, const double xt[], double fxt[], Nag_Comm *comm, Integer *info)
1: mx Integer Input
On entry: mx, the number of state variables.
2: n Integer Input
On entry: n, the number of sigma points.
3: xt[mx×n] const double Input
On entry: Xt, the sigma points generated in (a). For the jth sigma point, the value for the ith parameter is held in xt[(j-1)×mx+i-1], for i=1,2,,mx and j=1,2,,n.
4: fxt[mx×n] double Output
On exit: F(Xt) .
For the jth sigma point the value for the ith parameter should be held in fxt[(j-1)×mx+i-1], for i=1,2,,mx and j=1,2,,n.
5: comm Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to f.
userdouble *
iuserInteger *
pPointer 
The type Pointer will be void *. Before calling g13ekc you may allocate memory and initialize these pointers with various quantities for use by f when called from g13ekc (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
6: info Integer * Input/Output
On entry: info=0.
On exit: set info to a nonzero value if you wish g13ekc to terminate with fail.code= NE_USER_STOP.
Note: f should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by g13ekc. If your code inadvertently does return any NaNs or infinities, g13ekc is likely to produce unexpected results.
7: h function, supplied by the user External Function
The measurement function, H as described in (e).
The specification of h is:
void  h (Integer mx, Integer my, Integer n, const double yt[], double hyt[], Nag_Comm *comm, Integer *info)
1: mx Integer Input
On entry: mx, the number of state variables.
2: my Integer Input
On entry: my, the number of observed variables.
3: n Integer Input
On entry: n, the number of sigma points.
4: yt[mx×n] const double Input
On entry: Yt, the sigma points generated in (d). For the jth sigma point, the value for the ith parameter is held in yt[(j-1)×mx+i-1], for i=1,2,,mx and j=1,2,,n, where mx is the number of state variables and n is the number of sigma points.
5: hyt[my×n] double Output
On exit: H(Yt) .
For the jth sigma point the value for the ith parameter should be held in hyt[(j-1)×my+i-1], for i=1,2,,my and j=1,2,,n.
6: comm Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to h.
userdouble *
iuserInteger *
pPointer 
The type Pointer will be void *. Before calling g13ekc you may allocate memory and initialize these pointers with various quantities for use by h when called from g13ekc (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
7: info Integer * Input/Output
On entry: info=0.
On exit: set info to a nonzero value if you wish g13ekc to terminate with fail.code= NE_USER_STOP.
Note: h should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by g13ekc. If your code inadvertently does return any NaNs or infinities, g13ekc is likely to produce unexpected results.
8: x[mx] double Input/Output
On entry: x^t-1 the state vector for the previous time point.
On exit: x^t the updated state vector.
9: st[mx×mx] double Input/Output
Note: the (i,j)th element of the matrix is stored in st[(j-1)×mx+i-1].
On entry: St, such that St-1St-1T=Pt-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 the matrix stored in st is referenced.
On exit: St, the lower triangular part of a Cholesky factorization of the updated state covariance matrix.
10: comm Nag_Comm *
The NAG communication argument (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
11: fail NagError * Input/Output
The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

6 Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
See Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
NE_BAD_PARAM
On entry, argument value had an illegal value.
NE_INT
On entry, mx=value.
Constraint: mx1.
On entry, my=value.
Constraint: my1.
NE_INTERNAL_ERROR
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
See Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
NE_MAT_NOT_POS_DEF
A weight was negative and it was not possible to downdate the Cholesky factorization.
Unable to calculate the Cholesky factorization of the updated state covariance matrix.
Unable to calculate the Kalman gain matrix.
NE_NO_LICENCE
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library CL Interface for further information.
NE_USER_STOP
User requested termination in f.
User requested termination in h.

7 Accuracy

Not applicable.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.
g13ekc makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this function. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

None.

10 Example

This example implements the following nonlinear state space model, with the state vector x and state update function F given by:
mx =3 xt+1 = ( ξt+1 ηt+1 θt+1 ) T =F(xt)+vt = xt+ ( cosθt -sinθt 0 sinθt cosθt 0 001 ) ( 0.5r 0.5r 00 r/d -r/d ) ( ϕRt ϕLt ) +vt  
where r and d are known constants and ϕRt and ϕLt are time-dependent knowns. The measurement vector y and measurement function H is given by:
my =2 yt =(δt,αt)T =H(xt)+ut = ( Δ-ξtcosA-ηtsinA θt-A ) +ut  
where A and Δ are known constants. The initial values, x0 and P0, are given by
x0 = ( 0 0 0 ) , P0 = ( 0.100 00.10 000.1 )  
and the Cholesky factorizations of the error covariance matrices, Lx and Lx by
Lx = ( 0.100 00.10 000.1 ) , Ly = ( 0.010 00.01 )  

10.1 Program Text

Program Text (g13ekce.c)

10.2 Program Data

Program Data (g13ekce.d)

10.3 Program Results

Program Results (g13ekce.r)
The example described above can be thought of relating to the movement of a hypothetical robot. The unknown state, x, is the position of the robot (with respect to a reference frame) and facing, with (ξ,η) giving the x and y coordinates and θ the angle (with respect to the x-axis) that the robot is facing. The robot has two drive wheels, of radius r on an axle of length d. During time period t the right wheel is believed to rotate at a velocity of ϕRt and the left at a velocity of ϕLt. In this example, these velocities are fixed with ϕRt=0.4 and ϕLt=0.1. The state update function, F, calculates where the robot should be at each time point, given its previous position. However, in reality, there is some random fluctuation in the velocity of the wheels, for example, due to slippage. Therefore, the actual position of the robot and the position given by equation F will differ.
In the area that the robot is moving there is a single wall. The position of the wall is known and defined by its distance, Δ, from the origin and its angle, A, from the x-axis. The robot has a sensor that is able to measure y, with δ being the distance to the wall and α the angle to the wall. The measurement function H gives the expected distance and angle to the wall if the robot's position is given by xt. Therefore, the state space model allows the robot to incorporate the sensor information to update the estimate of its position.
GnuplotProduced by GNUPLOT 4.6 patchlevel 3 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 Example Program Illustration of Position and Orientation of a Hypothetical Robot Wall Position gnuplot_plot_1 Initial gnuplot_plot_2 Actual gnuplot_plot_3 Updated gnuplot_plot_4 gnuplot_plot_5 gnuplot_plot_6 gnuplot_plot_7