NAG FL Interface
g13ejf (kalman_​unscented_​state_​revcom)

Settings help

FL Name Style:


FL Specification Language:


1 Purpose

g13ejf applies the Unscented Kalman Filter to a nonlinear state space model, with additive noise.
g13ejf uses reverse communication for evaluating the nonlinear functionals of the state space model.

2 Specification

Fortran Interface
Subroutine g13ejf ( irevcm, mx, my, y, lx, ldlx, ly, ldly, x, st, ldst, n, xt, ldxt, fxt, ldfxt, ropt, lropt, icomm, licomm, rcomm, lrcomm, ifail)
Integer, Intent (In) :: mx, my, ldlx, ldly, ldst, ldxt, ldfxt, lropt, licomm, lrcomm
Integer, Intent (Inout) :: irevcm, n, icomm(licomm), ifail
Real (Kind=nag_wp), Intent (In) :: y(my), lx(ldlx,*), ly(ldly,*), ropt(lropt)
Real (Kind=nag_wp), Intent (Inout) :: x(mx), st(ldst,*), xt(ldxt,*), fxt(ldfxt,*), rcomm(lrcomm)
C Header Interface
#include <nag.h>
void  g13ejf_ (Integer *irevcm, const Integer *mx, const Integer *my, const double y[], const double lx[], const Integer *ldlx, const double ly[], const Integer *ldly, double x[], double st[], const Integer *ldst, Integer *n, double xt[], const Integer *ldxt, double fxt[], const Integer *ldfxt, const double ropt[], const Integer *lropt, Integer icomm[], const Integer *licomm, double rcomm[], const Integer *lrcomm, Integer *ifail)
The routine may be called by the names g13ejf or nagf_tsa_kalman_unscented_state_revcom.

3 Description

g13ejf 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 g13ejf 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.
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 (5) with:
Yt = [Xt   Xt,1+γΣx   Xt,1-γΣx] (13)
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.

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

Note: this routine uses reverse communication. Its use involves an initial entry, intermediate exits and re-entries, and a final exit, as indicated by the argument irevcm. Between intermediate exits and re-entries, all arguments other than fxt must remain unchanged.
1: irevcm Integer Input/Output
On initial entry: must be set to 0 or 3.
If irevcm=0, it is assumed that t=0, otherwise it is assumed that t0 and that g13ejf has been called at least once before at an earlier time step.
On intermediate exit: irevcm=1 or 2. The value of irevcm specifies what intermediate values are returned by this routine and what values the calling program must assign to arguments of g13ejf before re-entering the routine. Details of the output and required input are given in the individual argument descriptions.
On intermediate re-entry: irevcm must remain unchanged.
On final exit: irevcm=3
Constraint: irevcm=0, 1, 2 or 3.
Note: any values you return to g13ejf as part of the reverse communication procedure should not include floating-point NaN (Not a Number) or infinity values, since these are not handled by g13ejf. If your code does inadvertently return any NaNs or infinities, g13ejf is likely to produce unexpected results.
2: mx Integer Input
On entry: mx, the number of state variables.
Constraint: mx1.
3: my Integer Input
On entry: my, the number of observed variables.
Constraint: my1.
4: y(my) Real (Kind=nag_wp) array Input
On entry: yt, the observed data at the current time point.
5: lx(ldlx,*) Real (Kind=nag_wp) array Input
Note: the second dimension of the array lx must be at least mx.
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 lx is referenced.
If ldlx=0, there is no process noise (vt=0 for all t) and lx is not referenced.
If Σx is time dependent, the value supplied should be for time t.
6: ldlx Integer Input
On entry: the first dimension of the array lx as declared in the (sub)program from which g13ejf is called.
Constraint: ldlx=0 or ldlxmx.
7: ly(ldly,*) Real (Kind=nag_wp) array Input
Note: the second dimension of the array ly must be at least my.
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 ly is referenced.
If Σy is time dependent, the value supplied should be for time t.
8: ldly Integer Input
On entry: the first dimension of the array ly as declared in the (sub)program from which g13ejf is called.
Constraint: ldlymy.
9: x(mx) Real (Kind=nag_wp) array Input/Output
On initial entry: x^t-1 the state vector for the previous time point.
On intermediate exit: when
irevcm=1
x is unchanged.
irevcm=2
x^t¯.
On intermediate re-entry: x must remain unchanged.
On final exit: x^t the updated state vector.
10: st(ldst,*) Real (Kind=nag_wp) array Input/Output
Note: the second dimension of the array st must be at least mx.
On initial 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 st is referenced.
On intermediate exit: when
irevcm=1
st is unchanged.
irevcm=2
St¯, the lower triangular part of a Cholesky factorization of Pt¯.
On intermediate re-entry: st must remain unchanged.
On final exit: St, the lower triangular part of a Cholesky factorization of the updated state covariance matrix.
11: ldst Integer Input
On entry: the first dimension of the array st as declared in the (sub)program from which g13ejf is called.
Constraint: ldstmx.
12: n Integer Input/Output
On initial entry: the value used in the sizing of the fxt and xt arrays. The value of n supplied must be at least as big as the maximum number of sigma points that the algorithm will use. g13ejf allows sigma points to be calculated in two ways during the measurement update; they can be redrawn or augmented. Which is used is controlled by ropt.
If redrawn sigma points are used, then the maximum number of sigma points will be 2mx+1, otherwise the maximum number of sigma points will be 4mx+1.
On intermediate exit: the number of sigma points actually being used.
On intermediate re-entry: n must remain unchanged.
On final exit: reset to its value on initial entry.
Constraints:
if irevcm=0 or 3,
  • if redrawn sigma points are used, n2×mx+1;
  • otherwise n4×mx+1.
13: xt(ldxt,*) Real (Kind=nag_wp) array Input/Output
Note: the second dimension of the array xt must be at least max(my,n).
On initial entry: need not be set.
On intermediate exit: Xt when irevcm=1, otherwise Yt.
For the jth sigma point, the value for the ith parameter is held in xt(i,j), for i=1,2,,mx and j=1,2,,n.
On intermediate re-entry: xt must remain unchanged.
On final exit: the contents of xt are undefined.
14: ldxt Integer Input
On entry: the first dimension of the array xt as declared in the (sub)program from which g13ejf is called.
Constraint: ldxtmx.
15: fxt(ldfxt,*) Real (Kind=nag_wp) array Input/Output
Note: the second dimension of the array fxt must be at least n+max(mx,my).
On initial entry: need not be set.
On intermediate exit: the contents of fxt are undefined.
On intermediate re-entry: F(Xt) when irevcm=1, otherwise H(Yt) for the values of Xt and Yt held in xt.
For the jth sigma point the value for the ith parameter should be held in fxt(i,j), for j=1,2,,n. When irevcm=1, i=1,2,,mx and when irevcm=2, i=1,2,,my.
On final exit: the contents of fxt are undefined.
16: ldfxt Integer Input
On entry: the first dimension of the array fxt as declared in the (sub)program from which g13ejf is called.
Constraint: ldfxtmax(mx,my).
17: ropt(lropt) Real (Kind=nag_wp) array Input
On entry: optional parameters. The default value will be used for ropt(i) if lropt<i. Setting lropt=0 will use the default values for all optional parameters and ropt need not be set.
ropt(1)
If set to 1 then the second set of sigma points are redrawn, as given by equation (5). If set to 2 then the second set of sigma points are generated via augmentation, as given by equation (13).
Default is for the sigma points to be redrawn (i.e., ropt(1)=1)
ropt(2)
κx, value of κ used when constructing the first set of sigma points, Xt.
Defaults to 3-mx.
ropt(3)
αx, value of α used when constructing the first set of sigma points, Xt.
Defaults to 1.
ropt(4)
βx, value of β used when constructing the first set of sigma points, Xt.
Defaults to 2.
ropt(5)
Value of κ used when constructing the second set of sigma points, Yt.
Defaults to 3-2×mx when ldlx0 and the second set of sigma points are augmented and κx otherwise.
ropt(6)
Value of α used when constructing the second set of sigma points, Yt.
Defaults to αx.
ropt(7)
Value of β used when constructing the second set of sigma points, Yt.
Defaults to βx.
Constraints:
  • ropt(1)=1 or 2;
  • ropt(2)>-mx;
  • ropt(5)>-2×mx when ldly0 and the second set of sigma points are augmented, otherwiseropt(5)>-mx;
  • ropt(i)>0, for i=3,4,5,6.
18: lropt Integer Input
On entry: length of the options array ropt.
Constraint: 0lropt7.
19: icomm(licomm) Integer array Communication Array
On initial entry: icomm need not be set.
On intermediate exit: icomm is used for storage between calls to g13ejf.
On intermediate re-entry: icomm must remain unchanged.
On final exit: icomm is not defined.
20: licomm Integer Input
On entry: the length of the array icomm. If licomm is too small and licomm2 then ifail=201 is returned and the minimum value for licomm and lrcomm are given by icomm(1) and icomm(2) respectively.
Constraint: licomm30.
21: rcomm(lrcomm) Real (Kind=nag_wp) array Communication Array
On initial entry: rcomm need not be set.
On intermediate exit: rcomm is used for storage between calls to g13ejf.
On intermediate re-entry: rcomm must remain unchanged.
On final exit: rcomm is not defined.
22: lrcomm Integer Input
On entry: the length of the array rcomm. If lrcomm is too small and licomm2 then ifail=202 is returned and the minimum value for licomm and lrcomm are given by icomm(1) and icomm(2) respectively.
Suggested value: lrcomm = 30 + my + mx × my + (1+nb) × max(mx,my) , where nb is the optimal block size. In most cases a block size of 128 will be sufficient.
Constraint: lrcomm30+my+mx×my+2×max(mx,my).
23: ifail Integer Input/Output
On initial entry: ifail must be set to 0, −1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of −1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value −1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value −1 is recommended since useful values can be provided in some output arguments even when ifail0 on exit. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On final exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or −1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=11
On entry, irevcm=value.
Constraint: irevcm=0, 1, 2 or 3.
ifail=21
On entry, mx=value.
Constraint: mx1.
ifail=22
mx has changed between calls.
On intermediate entry, mx=value.
On initial entry, mx=value.
ifail=31
On entry, my=value.
Constraint: my1.
ifail=32
my has changed between calls.
On intermediate entry, my=value.
On initial entry, my=value.
ifail=61
On entry, ldlx=value and mx=value.
Constraint: ldlx=0 or ldlxmx.
ifail=81
On entry, ldly=value and my=value.
Constraint: ldlymy.
ifail=111
On entry, ldst=value and mx=value.
Constraint: ldstmx.
ifail=121
On entry, augmented sigma points requested, n=value and mx=value.
Constraint: nvalue.
ifail=122
On entry, redrawn sigma points requested, n=value and mx=value.
Constraint: nvalue.
ifail=123
n has changed between calls.
On intermediate entry, n=value.
On intermediate exit, n=value.
ifail=141
On entry, ldxt=value and mx=value.
Constraint: ldxtmx.
ifail=161
On entry, ldfxt=value and mx=value.
Constraint: if irevcm=1, ldfxtmx.
ifail=162
On entry, ldfxt=value and my=value.
Constraint: if irevcm=2, ldfxtmy.
ifail=171
On entry, ropt(1)=value.
Constraint: ropt(1)=1 or 2.
ifail=172
On entry, ropt(value)=value.
Constraint: κ>value.
ifail=173
On entry, ropt(value)=value.
Constraint: α>0.
ifail=181
On entry, lropt=value.
Constraint: 0lropt7.
ifail=191
icomm has been corrupted between calls.
ifail=201
On entry, licomm=value.
Constraint: licomm2.
icomm is too small to return the required array sizes.
ifail=202
On entry, licomm=value and lrcomm=value.
Constraint: licomm30 and lrcomm30+my+mx×my+2×max(mx,my).
The minimum required values for licomm and lrcomm are returned in icomm(1) and icomm(2) respectively.
ifail=211
rcomm has been corrupted between calls.
ifail=301
A weight was negative and it was not possible to downdate the Cholesky factorization.
ifail=302
Unable to calculate the Kalman gain matrix.
ifail=303
Unable to calculate the Cholesky factorization of the updated state covariance matrix.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

Not applicable.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.
g13ejf 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 routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

As well as implementing the Unscented Kalman Filter, g13ejf can also be used to apply the Unscented Transform (see Julier (2002)) to the function F, by setting ldlx=0 and terminating the calling sequence when irevcm=2 rather than irevcm=3. In this situation, on initial entry, x and st would hold the mean and Cholesky factorization of the covariance matrix of the untransformed sample and on exit (when irevcm=2) they would hold the mean and Cholesky factorization of the covariance matrix of the transformed sample.

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 (g13ejfe.f90)

10.2 Program Data

Program Data (g13ejfe.d)

10.3 Program Results

Program Results (g13ejfe.r)
The example described above can be thought of as 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 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