# LowLevelParticleFilters

This is a library for *state estimation*, that is, given measurements $y(t)$ from a dynamical system, estimate the state vector $x(t)$. Throughout, we assume dynamics on the form

\[x(t+1) = f(x(t), u(t), p, t, w(t))\\ y(t) = g(x(t), u(t), p, t, e(t))\]

where $x$ is the state vector, $u$ an input, $p$ some form of parameters, $t$ is the time and $w,e$ are disturbances (noise). Throughout the documentation, we often call the function $f$ `dynamics`

and the function $g$ `measurement`

.

The dynamics above describe a *discrete-time* system, i.e., the function $f$ takes the current state and produces the *next state*. This is in contrast to a *continuous-time* system, where $f$ takes the current state but produces the *time derivative* of the state. A continuous-time system can be *discretized*, described in detail in Discretization.

The parameters $p$ can be anything, or left out. You may write the dynamics functions such that they depend on $p$ and include parameters when you create a filter object. You may also override the parameters stored in the filter object when you call any function on the filter object. This behavior is modeled after the SciML ecosystem.

Depending on the nature of $f$ and $g$, the best method of estimating the state may vary. If $f,g$ are linear and the disturbances are additive and Gaussian, the `KalmanFilter`

is an optimal state estimator. If any of the above assumptions fail to hold, we may need to resort to more advanced estimators. This package provides several filter types, outlined below.

## Estimator types

We provide a number of filter types

`KalmanFilter`

. A standard Kalman filter. Is restricted to linear dynamics (possibly time varying) and Gaussian noise.`SqKalmanFilter`

. A standard Kalman filter on square-root form (slightly slower but more numerically stable with ill-conditioned covariance).`ExtendedKalmanFilter`

: For nonlinear systems, the EKF runs a regular Kalman filter on linearized dynamics. Uses ForwardDiff.jl for linearization. The noise model must be Gaussian.`UnscentedKalmanFilter`

: The Unscented Kalman filter often performs slightly better than the Extended Kalman filter but may be slightly more computationally expensive. The UKF handles nonlinear dynamics and measurement model, but still requires an additive Gaussian noise model and assumes all posterior distributions are Gaussian, i.e., can not handle multi-modal posteriors.`ParticleFilter`

: The particle filter is a nonlinear estimator. This version of the particle filter is simple to use and assumes that both dynamics noise and measurement noise are additive. Particle filters handle multi-modal posteriors.`AdvancedParticleFilter`

: This filter gives you more flexibility, at the expense of having to define a few more functions. This filter does not require the noise to be additive and is thus the most flexible filter type.`AuxiliaryParticleFilter`

: This filter is identical to`ParticleFilter`

, but uses a slightly different proposal mechanism for new particles.`DAEUnscentedKalmanFilter`

: An Unscented Kalman filter for differential-algebraic systems (DAE).

## Functionality

This package provides

- Filtering, estimating $x(t)$ given measurements up to and including time $t$. We call the filtered estimate $x(t|t)$ (read as $x$ at $t$ given $t$).
- Smoothing, estimating $x(t)$ given data up to $T > t$, i.e., $x(t|T)$.
- Parameter estimation.

All filters work in two distinct steps.

- The
*prediction*step (`predict!`

). During prediction, we use the dynamics model to form $x(t|t-1) = f(x(t-1), ...)$ - The
*correction*step (`correct!`

). In this step, we adjust the predicted state $x(t|t-1)$ using the measurement $y(t)$ to form $x(t|t)$.

In general, all filters represent not only a point estimate of $x(t)$, but a representation of the complete posterior probability distribution over $x$ given all the data available up to time $t$. One major difference between different filter types is how they represent these probability distributions.

# Particle filter

A particle filter represents the probability distribution over the state as a collection of samples, each sample is propagated through the dynamics function $f$ individually. When a measurement becomes available, the samples, called *particles*, are given a weight based on how likely the particle is given the measurement. Each particle can thus be seen as representing a hypothesis about the current state of the system. After a few time steps, most weights are inevitably going to be extremely small, a manifestation of the curse of dimensionality, and a resampling step is incorporated to refresh the particle distribution and focus the particles on areas of the state space with high posterior probability.

Defining a particle filter is straightforward, one must define the distribution of the noise `df`

in the dynamics function, `dynamics(x,u,p,t)`

and the noise distribution `dg`

in the measurement function `measurement(x,u,p,t)`

. Both of these noise sources are assumed to be additive, but can have any distribution. The distribution of the initial state `d0`

must also be provided. An example for a linear Gaussian system is given below.

`using LowLevelParticleFilters, LinearAlgebra, StaticArrays, Distributions, Plots`

Define problem

```
n = 2 # Dimension of state
m = 1 # Dimension of input
p = 1 # Dimension of measurements
N = 500 # Number of particles
const dg = MvNormal(p,1.0) # Measurement noise Distribution
const df = MvNormal(n,1.0) # Dynamics noise Distribution
const d0 = MvNormal(randn(n),2.0) # Initial state Distribution
```

Define linear state-space system (using StaticArrays for maximum performance)

```
const A = SA[0.97043 -0.097368
0.09736 0.970437]
const B = SA[0.1; 0;;]
const C = SA[0 1.0]
```

Next, we define the dynamics and measurement equations, they both take the signature `(x,u,p,t) = (state, input, parameters, time)`

```
dynamics(x,u,p,t) = A*x .+ B*u
measurement(x,u,p,t) = C*x
vecvec_to_mat(x) = copy(reduce(hcat, x)') # Helper function
```

the parameter `p`

can be anything, and is often optional. If `p`

is not provided when performing operations on filters, any `p`

stored in the filter objects (if supported) is used. The default if none is provided and none is stored in the filter is `p = SciMLBase.NullParameters()`

.

We are now ready to define and use a filter

`pf = ParticleFilter(N, dynamics, measurement, df, dg, d0)`

```
ParticleFilter{PFstate{StaticArraysCore.SVector{2, Float64}, Float64}, typeof(Main.dynamics), typeof(Main.measurement), Distributions.ZeroMeanIsoNormal{Tuple{Base.OneTo{Int64}}}, Distributions.ZeroMeanIsoNormal{Tuple{Base.OneTo{Int64}}}, Distributions.IsoNormal, DataType, Random.Xoshiro, SciMLBase.NullParameters}
state: PFstate{StaticArraysCore.SVector{2, Float64}, Float64}
dynamics: dynamics (function of type typeof(Main.dynamics))
measurement: measurement (function of type typeof(Main.measurement))
dynamics_density: Distributions.ZeroMeanIsoNormal{Tuple{Base.OneTo{Int64}}}
measurement_density: Distributions.ZeroMeanIsoNormal{Tuple{Base.OneTo{Int64}}}
initial_density: Distributions.IsoNormal
resample_threshold: Float64 0.1
resampling_strategy: LowLevelParticleFilters.ResampleSystematic <: LowLevelParticleFilters.ResamplingStrategy
rng: Random.Xoshiro
p: SciMLBase.NullParameters SciMLBase.NullParameters()
threads: Bool false
```

With the filter in hand, we can simulate from its dynamics and query some properties

```
du = MvNormal(m,1.0) # Random input distribution for simulation
xs,u,y = simulate(pf,200,du) # We can simulate the model that the pf represents
pf(u[1], y[1]) # Perform one filtering step using input u and measurement y
particles(pf) # Query the filter for particles, try weights(pf) or expweights(pf) as well
x̂ = weighted_mean(pf) # using the current state
```

```
2-element Vector{Float64}:
-0.9677269094523568
-1.3127763339779588
```

If you want to perform filtering using vectors of inputs and measurements, try any of the functions

```
sol = forward_trajectory(pf, u, y) # Filter whole vectors of signals
x̂,ll = mean_trajectory(pf, u, y)
plot(sol, xreal=xs, markersize=2)
```

`u`

ad `y`

are then assumed to be vectors of vectors. StaticArrays is recommended for maximum performance.

If MonteCarloMeasurements.jl is loaded, you may transform the output particles to `Matrix{MonteCarloMeasurements.Particles}`

with the layout `T × n_state`

using `Particles(x,we)`

. Internally, the particles are then resampled such that they all have unit weight. This is conventient for making use of the plotting facilities of MonteCarloMeasurements.jl.

For a full usage example, see the benchmark section below or example_lineargaussian.jl

### Resampling

The particle filter will perform a resampling step whenever the distribution of the weights has become degenerate. The resampling is triggered when the *effective number of samples* is smaller than `pf.resample_threshold`

$\in [0, 1]$, this value can be set when constructing the filter. How the resampling is done is governed by `pf.resampling_strategy`

, we currently provide `ResampleSystematic <: ResamplingStrategy`

as the only implemented strategy. See https://en.wikipedia.org/wiki/Particle_filter for more info.

# Particle Smoothing

Smoothing is the process of finding the best state estimate given both past and future data. Smoothing is thus only possible in an offline setting. This package provides a particle smoother, based on forward filtering, backward simulation (FFBS), example usage follows:

```
N = 2000 # Number of particles
T = 80 # Number of time steps
M = 100 # Number of smoothed backwards trajectories
pf = ParticleFilter(N, dynamics, measurement, df, dg, d0)
du = MvNormal(m,1) # Control input distribution
x,u,y = simulate(pf,T,du) # Simulate trajectory using the model in the filter
tosvec(y) = reinterpret(SVector{length(y[1]),Float64}, reduce(hcat,y))[:] |> copy
x,u,y = tosvec.((x,u,y))
xb,ll = smooth(pf, M, u, y) # Sample smooting particles
xbm = smoothed_mean(xb) # Calculate the mean of smoothing trajectories
xbc = smoothed_cov(xb) # And covariance
xbt = smoothed_trajs(xb) # Get smoothing trajectories
xbs = [diag(xbc) for xbc in xbc] |> vecvec_to_mat .|> sqrt
plot(xbm', ribbon=2xbs, lab="PF smooth")
plot!(vecvec_to_mat(x), l=:dash, lab="True")
```

We can plot the particles themselves as well

```
downsample = 5
plot(vecvec_to_mat(x), l=(4,), layout=(2,1), show=false)
scatter!(xbt[1, 1:downsample:end, :]', subplot=1, show=false, m=(1,:black, 0.5), lab="")
scatter!(xbt[2, 1:downsample:end, :]', subplot=2, m=(1,:black, 0.5), lab="")
```

# Kalman filter

The `KalmanFilter`

(wiki) assumes that $f$ and $g$ are linear functions, i.e., that they can be written on the form

\[x(t+1) = Ax(t) + Bu(t) + w(t)\\ y(t) = Cx(t) + Du(t) + e(t)\]

for some matrices $A,B,C,D$ where $w \sim N(0, R_1)$ and $e \sim N(0, R_2)$ are zero mean and Gaussian. The Kalman filter represents the posterior distributions over $x$ by the mean and a covariance matrix. The magic behind the Kalman filter is that linear transformations of Gaussian distributions remain Gaussian, and we thus have a very efficient way of representing them.

A Kalman filter is easily created using the constructor. Many of the functions defined for particle filters, are defined also for Kalman filters, e.g.:

```
R1 = cov(df)
R2 = cov(dg)
kf = KalmanFilter(A, B, C, 0, R1, R2, d0)
sol = forward_trajectory(kf, u, y) # filtered, prediction, pred cov, filter cov, loglik
```

It can also be called in a loop like the `pf`

above

```
for t = 1:T
kf(u,y) # Performs both correct and predict!!
# alternatively
ll, e = correct!(kf, y, nothing, t) # Returns loglikelihood and prediction error
x = state(kf)
R = covariance(kf)
predict!(kf, u, nothing, t)
end
```

The matrices in the Kalman filter may be *time varying*, such that `A[:, :, t]`

is $A(t)$. They may also be provided as functions on the form $A(t) = A(x, u, p, t)$. This works for both dynamics and covariance matrices.

The numeric type used in the Kalman filter is determined from the mean of the initial state distribution, so make sure that this has the correct type if you intend to use, e.g., `Float32`

or `ForwardDiff.Dual`

for automatic differentiation.

## Smoothing using KF

Kalman filters can also be used for smoothing

```
kf = KalmanFilter(A, B, C, 0, cov(df), cov(dg), d0)
xT,R,lls = smooth(kf, u, y, p) # Smoothed state, smoothed cov, loglik
```

Plot and compare PF and KF

```
plot(vecvec_to_mat(xT), lab="Kalman smooth", layout=2)
plot!(xbm', lab="pf smooth")
plot!(vecvec_to_mat(x), lab="true")
```

# Unscented Kalman Filter

The `UnscentedKalmanFilter`

represents posterior distributions over $x$ as Gaussian distributions, but propagate them through a nonlinear function $f$ by a deterministic sampling of a small number of particles called *sigma points* (this is referred to as the *unscented transform*). This UKF thus handles nonlinear functions $f,g$, but only Gaussian disturbances and unimodal posteriors.

The UKF takes the same arguments as a regular `KalmanFilter`

, but the matrices defining the dynamics are replaced by two functions, `dynamics`

and `measurement`

, working in the same way as for the `ParticleFilter`

above.

`ukf = UnscentedKalmanFilter(dynamics, measurement, cov(df), cov(dg), MvNormal([1.,1.]), nu=m, ny=p)`

```
UnscentedKalmanFilter{typeof(Main.dynamics), typeof(Main.measurement), Matrix{Float64}, Matrix{Float64}, Distributions.ZeroMeanDiagNormal{Tuple{Base.OneTo{Int64}}}, StaticArraysCore.SVector{2, Float64}, Vector{Float64}, Matrix{Float64}, SciMLBase.NullParameters}
dynamics: dynamics (function of type typeof(Main.dynamics))
measurement: measurement (function of type typeof(Main.measurement))
R1: Array{Float64}((2, 2)) [1.0 0.0; 0.0 1.0]
R2: Array{Float64}((1, 1)) [1.0;;]
d0: Distributions.ZeroMeanDiagNormal{Tuple{Base.OneTo{Int64}}}
xs: Array{StaticArraysCore.SVector{2, Float64}}((5,))
x: Array{Float64}((2,)) [0.0, 0.0]
R: Array{Float64}((2, 2)) [1.0 0.0; 0.0 1.0]
t: Base.RefValue{Int64}
ny: Int64 1
nu: Int64 1
p: SciMLBase.NullParameters SciMLBase.NullParameters()
```

If your function `dynamics`

describes a continuous-time ODE, do not forget to **discretize** it before passing it to the UKF. See Discretization for more information.

## UKF for DAE systems

See the docstring for `DAEUnscentedKalmanFilter`

or the test file. This filter is modeled after

"Nonlinear State Estimation of Differential Algebraic Systems" Ravi K. Mandela, Raghunathan Rengaswamy, Shankar Narasimhan

This filter is still considered experimental and subject to change without respecting semantic versioning. Use at your own risk. The `AdvancedParticleFilter`

also supports DAE systems, see this tutorial.

# Extended Kalman Filter

The `ExtendedKalmanFilter`

(EKF) is similar to the UKF, but propagates Gaussian distributions by linearizing the dynamics and using the formulas for linear systems similar to the standard Kalman filter. This can be slightly faster than the UKF (not always), but also less accurate for strongly nonlinear systems. The linearization is performed automatically using ForwardDiff.jl. In general, the UKF is recommended over the EKF unless the EKF is faster and computational performance is the top priority.

The EKF constructor has the following two signatures

```
ExtendedKalmanFilter(dynamics, measurement, R1,R2,d0=MvNormal(Matrix(R1)); nu::Int, p = SciMLBase.NullParameters(), α = 1.0, check = true)
ExtendedKalmanFilter(kf, dynamics, measurement)
```

The first constructor takes all the arguments required to initialize the extended Kalman filter, while the second one takes an already defined standard Kalman filter. using the first constructor, the user must provide the number of inputs to the system, `nu`

.

where `kf`

is a standard `KalmanFilter`

from which the covariance properties are taken.

If your function `dynamics`

describes a continuous-time ODE, do not forget to **discretize** it before passing it to the UKF. See Discretization for more information.

# AdvancedParticleFilter

The `AdvancedParticleFilter`

works very much like the `ParticleFilter`

, but admits more flexibility in its noise models.

The `AdvancedParticleFilter`

type requires you to implement the same functions as the regular `ParticleFilter`

, but in this case you also need to handle sampling from the noise distributions yourself. The function `dynamics`

must have a method signature like below. It must provide one method that accepts state vector, control vector, parameter, time *and* `noise::Bool`

that indicates whether or not to add noise to the state. If noise should be added, this should be done inside `dynamics`

An example is given below

```
using Random
const rng = Random.Xoshiro()
function dynamics(x, u, p, t, noise=false) # It's important that this defaults to false
x = A*x .+ B*u # A simple linear dynamics model in discrete time
if noise
x += rand(rng, df) # it's faster to supply your own rng
end
x
end
```

The `measurement_likelihood`

function must have a method accepting state, input, measurement, parameter and time, and returning the log-likelihood of the measurement given the state, a simple example below:

```
function measurement_likelihood(x, u, y, p, t)
logpdf(dg, C*x-y) # A simple linear measurement model with normal additive noise
end
```

This gives you very high flexibility. The noise model in either function can, for instance, be a function of the state, something that is not possible for the simple `ParticleFilter`

. To be able to simulate the `AdvancedParticleFilter`

like we did with the simple filter above, the `measurement`

method with the signature `measurement(x,u,p,t,noise=false)`

must be available and return a sample measurement given state (and possibly time). For our example measurement model above, this would look like this

`measurement(x, u, p, t, noise=false) = C*x + noise*rand(rng, dg)`

We now create the `AdvancedParticleFilter`

and use it in the same way as the other filters:

```
apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0)
sol = forward_trajectory(apf, u, y, p)
```

```
LowLevelParticleFilters.ParticleFilteringSolution{AdvancedParticleFilter{PFstate{StaticArraysCore.SVector{2, Float64}, Float64}, typeof(Main.dynamics), typeof(Main.measurement), typeof(Main.measurement_likelihood), Distributions.ZeroMeanIsoNormal{Tuple{Base.OneTo{Int64}}}, Distributions.IsoNormal, DataType, Random.Xoshiro, SciMLBase.NullParameters}, Vector{StaticArraysCore.SVector{1, Float64}}, Vector{StaticArraysCore.SVector{1, Float64}}, Matrix{StaticArraysCore.SVector{2, Float64}}, Matrix{Float64}, Matrix{Float64}, Float64}(AdvancedParticleFilter{PFstate{StaticArraysCore.SVector{2, Float64}, Float64}, typeof(Main.dynamics), typeof(Main.measurement), typeof(Main.measurement_likelihood), Distributions.ZeroMeanIsoNormal{Tuple{Base.OneTo{Int64}}}, Distributions.IsoNormal, DataType, Random.Xoshiro, SciMLBase.NullParameters}
state: PFstate{StaticArraysCore.SVector{2, Float64}, Float64}
dynamics: dynamics (function of type typeof(Main.dynamics))
measurement: measurement (function of type typeof(Main.measurement))
measurement_likelihood: measurement_likelihood (function of type typeof(Main.measurement_likelihood))
dynamics_density: Distributions.ZeroMeanIsoNormal{Tuple{Base.OneTo{Int64}}}
initial_density: Distributions.IsoNormal
resample_threshold: Float64 0.5
resampling_strategy: LowLevelParticleFilters.ResampleSystematic <: LowLevelParticleFilters.ResamplingStrategy
rng: Random.Xoshiro
p: SciMLBase.NullParameters SciMLBase.NullParameters()
threads: Bool false
, StaticArraysCore.SVector{1, Float64}[[0.8795109901273414], [1.3046955972015906], [0.31705874052626987], [1.444574185590547], [-1.1491929791999564], [0.27640665759698513], [-0.5035248207714105], [-1.6006018823293713], [-0.4793351537555876], [-0.5024550001909771] … [-0.34182058078462335], [0.6539941847475343], [-0.8008052990735235], [0.09761753528688945], [-0.6867682494810022], [0.1760239368986566], [-0.3552647409054233], [-1.125615776353175], [0.7346625752534395], [-2.1130777684696396]], StaticArraysCore.SVector{1, Float64}[[0.2714249204411139], [1.2242414292020871], [0.24834069571121886], [-2.3457196118588444], [-3.407214468935609], [-3.383093905244049], [-4.377762937681354], [-2.1168334025813076], [-2.249834799911277], [-0.9937550204436898] … [-2.8017764331371033], [-6.561340944966924], [-5.43425759430829], [-6.140471723865313], [-7.048419880158736], [-6.071152086873574], [-8.16645749639666], [-5.728182642565605], [-3.809027284778431], [-3.1190392495825408]], StaticArraysCore.SVector{2, Float64}[[-0.6145965479863531, -0.7404025168509196] [1.3342329656042855, -0.9427677400054031] … [8.002265302489317, -2.7060597040276013] [6.986038119674148, -3.578004705733763]; [-1.3989836396738076, -1.5417235968188012] [-0.6513995686671716, -0.971501006956388] … [-2.2308717962274605, -6.991295836379237] [8.037581845220743, -0.8013121695395191]; … ; [-1.0783309582777498, -1.1296198455735036] [0.3461236471694079, 0.20169112200337125] … [-0.7382980166266768, -9.015945433727929] [1.8954077341943663, -5.658933422928399]; [-0.2686556045256273, -0.06099057765022367] [-0.9777714617554427, -0.7115122892375569] … [1.5041220350996678, -6.2695327974752875] [-0.09407512410624386, -6.526049807262984]], [-7.331718866542534 -8.964771584273159 … -6.5708585572802605 -7.020812667573451; -8.463575257934966 -10.159306028924501 … -11.705013027283654 -9.601417431385777; … ; -7.801284703302528 -7.609177616694578 … -19.309560128402545 -10.141019228501523; -6.875071516799707 -8.033731094042118 … -9.91605455742218 -12.719348492676403], [0.0006544477115127778 0.0001278348215203822 … 0.0014005944052620871 0.0008930994049285459; 0.00021101628277096852 3.871412532967075e-5 … 8.252346058767788e-6 6.76328036591976e-5; … ; 0.0004092089290805038 0.0004958794938554183 … 4.1111619791732506e-9 3.942859559044348e-5; 0.0010332237431497926 0.0003243358211106567 … 4.937558150863921e-5 2.9926583279006314e-6], -1107.119252397841)
```

`plot(sol, xreal=x)`

We can even use this type as an AuxiliaryParticleFilter

```
apfa = AuxiliaryParticleFilter(apf)
sol = forward_trajectory(apfa, u, y, p)
plot(sol, dim=1, xreal=x) # Same as above, but only plots a single dimension
```

See the tutorials section for more advanced examples, including state estimation for DAE (Differential-Algebraic Equation) systems.

# Troubleshooting and tuning

Tuning a particle filter can be quite the challenge. To assist with this, we provide som visualization tools

`debugplot(pf,u[1:20],y[1:20], runall=true, xreal=x[1:20])`

```
Time Surviving Effective nbr of particles
--------------------------------------------------------------
t: 1 1.000 2000.0
t: 2 1.000 1209.0
t: 3 1.000 761.3
t: 4 1.000 557.1
t: 5 0.227 2000.0
t: 6 1.000 900.7
t: 7 1.000 592.1
t: 8 1.000 355.7
t: 9 1.000 208.3
t: 10 0.184 2000.0
t: 11 1.000 1103.0
t: 12 1.000 709.3
t: 13 1.000 508.3
t: 14 1.000 264.8
t: 15 0.237 2000.0
t: 16 1.000 1549.8
t: 17 1.000 481.2
t: 18 1.000 209.9
t: 19 0.091 2000.0
t: 20 1.000 1529.2
GKS: Ghostscript support not compiled in
```

The plot displays all states and all measurements. The heatmap in the background represents the weighted particle distributions per time step. For the measurement sequences, the heatmap represent the distributions of predicted measurements. The blue dots corresponds to measured values. In this case, we simulated the data and we had access to states as well, if we do not have that, just omit `xreal`

. You can also manually step through the time-series using

`commandplot(pf,u,y; kwargs...)`

For options to the debug plots, see `?pplot`

.

Something seems to be off with this figure as the hottest spot is not really where we would expect it

Optimization of the log likelihood can be done by, e.g., global/black box methods, see BlackBoxOptim.jl, see examples in Parameter Estimation. Standard tricks apply, such as performing the parameter search in log-space etc.