Tutorial 2: Duffing oscillator

Install required packages

using Pkg
Pkg.add(["EnKF", "Distributions", "ProgressMeter", "OrdinaryDiffEq", "LinearAlgebra"])

17.9s

Julia

using EnKF
using Distributions
using LinearAlgebra
using ProgressMeter
using OrdinaryDiffEq
using Plots

213.7s

Julia

We are interested in simulating the Duffing oscillator with forcing: with and

First, let us solve the dynamical equation using DifferentialEquations.jl with on with a Runge-Kutta 4 method (denoted RK4 in DifferentialEquations.jl)

We also define an integrator that will can propagate our solution from time to

function duffing(du,u,p,t)
 du[1] = u[2]
 du[2] = u[1] - u[1]^3 -γ*u[2] + d*cos(ω*t)
end

γ = 0.1
d = 0.1
ω = 1.4

u0 = [1.0; -1.0]
tspan = (0.0,50.0)

Δt = 1e-2
T = tspan[1]:Δt:tspan[end]

prob = ODEProblem(duffing, u0, tspan)
sol = solve(prob, RK4(), adaptive = false, dt = Δt)

integrator = init(prob, RK4(), adaptive =false, dt = Δt, save_everystep=false)

1.5s

Julia

t: 0.0 u: 2-element Array{Float64,1}: 1.0 -1.0

Plot solution

plot(sol, vars = (1))

31.8s

Julia

Tools for the ensemble Kalman filter

The ensemble Kalman filter known under the sobriquet EnKF is a Monte-Carlo view of the Kalman filter suited for large and nonlinear systems introduced by Evensen, et al.(1999)

Let us define the propagation function:

function (::PropagationFunction)(t::Float64, ENS::EnsembleState{N, TS}) where {N, TS}
    for (i,s) in enumerate(ENS.S)
        
        set_t!(integrator, deepcopy(t))
        set_u!(integrator, deepcopy(s))
#         for j=1:10
        step!(integrator)
#         end
        ENS.S[i] = deepcopy(integrator.u)

    end
    
    return ENS
end

0.0s

Julia

fprop = PropagationFunction()

0.0s

Julia

PropagationFunction()

Define measurement function m

function (::MeasurementFunction)(t::Float64, s::TS) where TS
    return [s[2]]
end

0.0s

Julia

The measurement matrix is defined by the following method :

function (::MeasurementFunction)(t::Float64) 
    return reshape([0.0, 1.0],(1,2))
end

0.0s

Julia

m = MeasurementFunction()

0.0s

Julia

MeasurementFunction()

Define real measurement function z

function (::RealMeasurementFunction)(t::Float64, ENS::EnsembleState{N, TZ}) where {N, TZ}
    let s = sol(t)
    fill!(ENS, [deepcopy(s[2])])
    end
    return ENS
end

0.0s

Julia

z = RealMeasurementFunction()

0.0s

Julia

RealMeasurementFunction()

Define covariance inflation

A = MultiAdditiveInflation(2, 1.05, MvNormal(zeros(2), 2.0*I))
# A = IdentityInflation()

0.0s

Julia

MultiAdditiveInflation{2}(1.05, IsoNormal( dim: 2 μ: [0.0, 0.0] Σ: [4.0 0.0; 0.0 4.0] ) )

Define noise covariance

ϵ = AdditiveInflation(MvNormal(zeros(1), 1.0*I))

0.0s

Julia

AdditiveInflation{1}(IsoNormal( dim: 1 μ: [0.0] Σ: [1.0] ) )

Define number of ensemble members, size of the measurement vector, and booleans isinflated, isfiltered and isaugmented

N = 50
NZ = 1
isinflated = false
isfiltered = false
isaugmented = false

0.0s

Julia

false

Initialized our ensemble, the mean estimate of the initial condition is

x₀ = [0.5, -0.5]
ens = initialize(N, MvNormal(x₀, 2.0*I))
estimation_state = [deepcopy(ens.S)]

true_state = [deepcopy(x₀)]
covs = []

0.0s

Julia

0-element Array{Any,1}

There is no filtering of the state variable therefore g = FilteringFunction(), i.e no action.

g = FilteringFunction()

0.0s

Julia

FilteringFunction()

Define the ENKF wrapper:

enkf = ENKF(N, NZ, fprop, A, g, m, z, ϵ, isinflated, isfiltered, isaugmented)

0.0s

Julia

ENKF{50,1}(PropagationFunction(), MultiAdditiveInflation{2}(1.05, IsoNormal( dim: 2 μ: [0.0, 0.0] Σ: [4.0 0.0; 0.0 4.0] ) ), FilteringFunction(), MeasurementFunction(), RealMeasurementFunction(), AdditiveInflation{1}(IsoNormal( dim: 1 μ: [0.0] Σ: [1.0] ) ), false, false, false)

Ensemble Kalman filter estimation

Δt = 1e-2
Tsub = 0.0:Δt:50.0-Δt

for (n,t) in enumerate(Tsub)

    global ens
    
    t, ens, cov = enkf(t, Δt, ens)
    
    push!(estimation_state, deepcopy(ens.S))
    push!(covs, deepcopy(cov))
    

end

2.4s

Julia

Plot state estimated from EnKF and true state .

s =  hcat(sol(T).u...)
ŝ =  hcat(mean.(estimation_state)...)

plt = plot(layout = (2, 1), legend = :bottomright)
plot!(plt[1], T, s[1,1:end], linewidth = 3, label = "truth")
plot!(plt[1], Tsub, ŝ[1,1:end-1], linewidth = 3, markersize = 2, label = "EnKF mean", xlabel = "t", ylabel = "x", linestyle =:dash)

plot!(plt[2], T, s[2,1:end], linewidth = 3, label = "truth")
plot!(plt[2], Tsub, ŝ[2,1:end-1], linewidth = 3, markersize = 2, label = "EnKF mean", xlabel = "t", ylabel = "y", linestyle =:dash)

0.9s

Julia

plot(s[1,:], s[2,:], linewidth = 3, label = "truth", legend = true)
plot!(ŝ[1,1:end-1], ŝ[2,1:end-1], linewidth = 3, label = "EnKF mean", xlabel = "x", 
    ylabel = "y", linestyle = :dash, legend = false)#:bottomleft)

0.7s

Julia

Plot diagonal components of the a posteriori covariance

plot(Tsub, map(covs) do P  
    P[1,1]+ eps() end, yscale = :log10, linewidth = 3, label ="sigma_{xx}")
plot!(Tsub, map(covs) do P
    P[2,2]+ eps() end, yscale = :log10, linewidth = 3, xlabel = "t", label = "sigma_{yy}")

0.8s

Julia

Shift+Enter to run

Julia