# Kalman.jl

The [Kalman.jl](https://github.com/wkearn/Kalman.jl) package aims to provide a Julian interface for various kinds of Kalman filter.


```
using Kalman, Winston
```

The Kalman filter is a way to estimate the state of a system given some knowledge of the (possibly stochastic) dynamics of that system and some noisy data collected from the system. Basically you take your state estimate, run it through the system model to advance it in time (along with the covariance matrix of your estimate), then you update your state estimate and state covariance matrix with the information gained from the data you have collected.

We'll try out the Kalman filter using an example from [Welch and Bishop, 2006](http://www.cs.unc.edu/~welch/media/pdf/kalman_intro.pdf).

Imagine we have a voltmeter taking measurements of a noisy but, on average, constant circuit. We first need a linear state-space model of this system that looks like $ x_{k+1} = Ax_k+w_k$. Because the voltage of the circuit doesn't change from time to time, we set $A=1$ and find our model $ x_{k+1} = x_k + w_k$ where $w_k$ is a vector of white Gaussian noise with variance $Q$.

The second half of the Kalman filter model is an observation model which relates the measurements you take from the system to the system state. It takes the form $ z_{k} = Hx_k + v_k$ where $v_k$ is white Gaussian noise with variance $R$. For our example, we just take direct measurements of the voltage, so $H=1$.


## Types

Now we look at how we can represent those models in Julia.


```
type LinearModel <: Model
    a::Matrix
    g::Matrix
    q::Matrix
end
```
A `LinearModel` takes a matrix `a` representing the state-transition matrix, `g` which explains how the process noise affects the system state, and `q` which is the process noise covariance matrix.

We can set one of those up for our model now. Our two transition matrices are just 1 as above. We'll give our process some small (but non-zero!) noise. The circuit isn't going to be perfectly constant.


```
f = LinearModel([1]',[1]',[1e-10]')
```




    LinearModel(1x1 Array{Int64,2}:
     1,1x1 Array{Int64,2}:
     1,1x1 Array{Float64,2}:
     1.0e-10)



The observation model is represented with a `LinearObservationModel` type

```
type LinearObservationModel <: ObservationModel
    h::Matrix
    r::Matrix
end
```

with definitions as above. For our example h is just `[1]'`, and we'll give our voltmeter $0.1\ \mathrm{V}$ RMS noise.


```
z = LinearObservationModel([1]',[0.01]')
```




    LinearObservationModel(1x1 Array{Int64,2}:
     1,1x1 Array{Float64,2}:
     0.01)



Finally, we need a way to represent the state of our system. To completely describe the state, we need both a state vector and a state covariance matrix

```
type State
    x::Vector
    p::Matrix
end
```

We have to seed our Kalman filter with an initial state and covariance estimate. It turns out that the covariance will converge regardless of the starting value. For our example as well, the state will converge, but there are systems with dynamics which might lead the state into a local minimum, so it's important to make a reasonable estimate for the state.

We'll start with $\hat{x}_0 = 0$ and $P_0 = 0.5$


```
x0 = State([0.0],[1.0]')
```




    State{Float64}([0.0],1x1 Array{Float64,2}:
     1.0)



Now we have everything to make our Kalman Filter. We just stick the two models and the state into a `BasicKalmanFilter` type

```
type BasicKalmanFilter <: LinearKalmanFilter
    x::State
    f::LinearModel
    z::LinearObservationModel
    adv::Bool
end
```

The `adv` field is a bad way of ensuring that we first advance our filter in time and then update it with some measurements. That may go away.


```
kf = BasicKalmanFilter(x0,f,z,false)
```




    BasicKalmanFilter(State{Float64}([0.0],1x1 Array{Float64,2}:
     1.0),LinearModel(1x1 Array{Int64,2}:
     1,1x1 Array{Int64,2}:
     1,1x1 Array{Float64,2}:
     1.0e-10),LinearObservationModel(1x1 Array{Int64,2}:
     1,1x1 Array{Float64,2}:
     0.01),false)



Now we need some data. Data to be fed to the Kalman filter is wrapped in an `Observation` type (which should probably just be a typealias for a Vector). Let's assume that our true voltage $x = -0.37727\ \mathrm{V}$.


```
x = -0.37727
y = map(y->Observation([y]),x+0.1*randn(50))
```




    50-element Array{Observation{Float64},1}:
     Observation{Float64}([-0.556689])
     Observation{Float64}([-0.28849]) 
     Observation{Float64}([-0.429203])
     Observation{Float64}([-0.27673]) 
     Observation{Float64}([-0.461203])
     Observation{Float64}([-0.348198])
     Observation{Float64}([-0.587995])
     Observation{Float64}([-0.333402])
     Observation{Float64}([-0.413611])
     Observation{Float64}([-0.502472])
     Observation{Float64}([-0.43491]) 
     Observation{Float64}([-0.465556])
     Observation{Float64}([-0.275459])
     ⋮                                
     Observation{Float64}([-0.4012])  
     Observation{Float64}([-0.494106])
     Observation{Float64}([-0.516698])
     Observation{Float64}([-0.389034])
     Observation{Float64}([-0.283263])
     Observation{Float64}([-0.198191])
     Observation{Float64}([-0.318942])
     Observation{Float64}([-0.229493])
     Observation{Float64}([-0.260345])
     Observation{Float64}([-0.485391])
     Observation{Float64}([-0.297279])
     Observation{Float64}([-0.203137])



Now we can run our filter by sequentially calling `predict(kf::BasicKalmanFilter)` and `update(kf::BasicKalmanFilter,y::Observation)` or by simply calling `predictupdate(kf,y)`. We'll also store the state and variance information and plot it up with Winston.


```
xs = x0.x[1]*ones(50)
ps = x0.p[1]*ones(50)

for i = 1:49
    predictupdate!(kf,y[i])
    xs[i+1] = kf.x.x[1]
    ps[i+1] = kf.x.p[1]
end
```


```
plot([0,50],[x,x],"b")
oplot(map(y->y.y[1],y),"r.")
oplot(xs,"g")
xlim(0,50)
```


```
plot(ps[2:50],"g")
```

## Unscented Kalman filters.

Unscented Kalman filters let you apply the prediction-update method of the Kalman filter to nonlinear models. They do this by deterministically sampling a set of "sigma points" which are then run through the filter equations to reconstruct the state and covariance matrices after a nonlinear transformation.

Unscented Kalman filters aren't pushed to Github yet, and there are some restrictions on the kind of model you can use for it (additive Gaussian noise for now; you might be able to get away with non-Gaussian noise, but multiplicative noise definitely won't work).

Still, let's see if we can't make them work. We'll take the example from Kandepu et al. (2008), a van der Pol Oscillator.

The equations of motion are

\begin{align}
\dot{x}_1 &= -x_2 \\
\dot{x}_2 &= -\mu (1-x_1^2)x^2 + x_1
\end{align}

where $\mu=0.2$.


```
include(Pkg.dir("Kalman","sandbox","unscented.jl"))
```




    predictupdate! (generic function with 3 methods)




```
const μ = 0.2

x0 = State([0.,5],5.0*eye(2))
Q = 1e-3*eye(2)
R = [0.1 0;
    0 1e-3]

function f2(x::Vector,dt::Float64)
    x1 = zeros(x)
    x1[1] = x[1] + dt*-x[2]
    x1[2] = x[2] + dt*(-μ*(1-x[1]^2)*x[2]+x[1])
    x1
end

h(x) = x

dt = 0.01

xs = fill(zeros(x0.x),4000)
xs[1] = [1.4,0]
x = [1.4,0]
for i = 2:4000
    x = f2(x,dt)
    xs[i] = x
end

ys = xs .+ map(y->sqrt(R)*randn(2),1:4000)
```




    4000-element Array{Array{Float64,1},1}:
     [1.09289,-0.0238201]   
     [1.46761,0.0453131]    
     [1.35845,0.0130256]    
     [1.10434,0.0656155]    
     [1.18755,0.0145539]    
     [1.32602,0.0972942]    
     [1.25023,0.136334]     
     [0.902791,0.112868]    
     [1.22303,0.0854962]    
     [1.11071,0.148996]     
     [1.36694,0.151589]     
     [1.71204,0.152962]     
     [1.61606,0.172073]     
     ⋮                      
     [0.0785064,0.015488]   
     [-0.771996,0.0356125]  
     [0.339114,0.0404495]   
     [0.108062,0.0174942]   
     [0.323293,0.00062621]  
     [-0.0526205,0.041678]  
     [-0.0474085,0.00729233]
     [0.202621,0.0664006]   
     [-0.329399,0.0489221]  
     [-0.157211,0.0354753]  
     [0.131959,0.0343319]   
     [0.109813,0.0808878]   




```
fm = AdditiveUnscentedModel(x->f2(x,dt),Q)
zm = AdditiveUnscentedObservationModel(h,R)
p0 = 5.0*eye(2)
kf = AdditiveUnscentedKalmanFilter(State([0,5.0],p0),fm,zm,0.1,2.0,0.0)

xs1 = 0.0*ones(4000)
xs2 = 5.0*ones(4000)
ps = 5.0*ones(4000)

for i = 2:4000
    predictupdate!(kf,Observation([ys[i]]))
    xs1[i] = kf.x.x[1]
    xs2[i] = kf.x.x[2]
    ps[i] = kf.x.p[1]
end
```


```
plot(xs1,"k")
```

85bf20179627c8dc3455f1e43fa346b6b25df13c
ipynb
Jupyter Notebook
examples/CAJUN_09-25-2014.ipynb
wkearn/Kalman.jl
a080e4e18b2f27af60fa3142be843d09541cbf7d
MIT
2015-02-11T14:27:24.000Z
2021-07-30T20:51:38.000Z
examples/CAJUN_09-25-2014.ipynb
wkearn/Kalman.jl
a080e4e18b2f27af60fa3142be843d09541cbf7d
MIT
2015-10-02T11:47:41.000Z
2018-02-02T09:45:46.000Z
examples/CAJUN_09-25-2014.ipynb
wkearn/Kalman.jl
a080e4e18b2f27af60fa3142be843d09541cbf7d
MIT
2015-03-18T21:46:30.000Z
2021-07-30T20:51:42.000Z
Qwen/Qwen-72B
1. YES
2. YES