Linear Time Continuous System

Similar to the linear time discrete example, we will now estimate a linear time continuous system $\partial_t u = A u$. We simulate the correspoding system using OrdinaryDiffEq.jl and generate a ContinuousDataDrivenProblem from the simulated data.

using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq
using Plots

A = [-0.9 0.2; 0.0 -0.2]
u0 = [10.0; -10.0]
tspan = (0.0, 10.0)

f(u,p,t) = A*u

sys = ODEProblem(f, u0, tspan)
sol = solve(sys, Tsit5(), saveat = 0.05);

We could use the DESolution to define our problem, but here we want to use the data for didactic purposes. For a ContinuousDataDrivenProblem, we need either the state trajectory and the timepoints or the state trajectory and its derivate.

X = Array(sol)
t = sol.t
prob = ContinuousDataDrivenProblem(X, t)

Continuous DataDrivenProblem{Float64} ##DDProblem#410 in 2 dimensions and 201 samples

And plot the problems data.

plot(prob)

We can see that the derivative has been automatically added via a collocation method, which defaults to a LinearInterpolation. We can do a visual check and compare our derivatives with the interpolation of the ODESolution.

DX = Array(sol(t, Val{1}))
scatter(t, DX', label = ["Solution" nothing], color = :red, legend = :bottomright)
plot!(t, prob.DX', label = ["Linear Interpolation" nothing], color = :black)

Since we have a linear system, we can use gDMD, which approximates the generator of the dynamics

res = solve(prob, DMDSVD())
println(res)

Linear Solution with 2 equations and 4 parameters.
Returncode: solved
L₂ Norm error : [0.23746883585659814, 0.0003580043828573712]
AIC : [-280.9814951941843, -1586.928031133052]
R² : [0.9997919083335011, 0.9999925518753656]

We see that the system has been recovered correctly, indicated by the small error and high AIC score of the result. We can confirm this by looking at the resulting Basis

system = result(res)
println(system)

Model ##Koopman#422 with 2 equations
States : x₁ x₂
Parameters : p₁ p₂ p₃ p₄
Independent variable: t
Equations
Differential(t)(x₁) = p₁*x₁ + p₂*x₂
Differential(t)(x₂) = p₃*x₁ + p₄*x₂

And also plot the prediction of the recovered dynamics

plot(res)

Or a have a look at the metrics of the result

metrics(res)

Dict{Any, Any} with 3 entries:
  :AIC => [-280.981, -1586.93]
  :R²  => [0.999792, 0.999993]
  :L₂  => [0.237469, 0.000358004]

And check the parameters of the result

parameters(res)

4-element Vector{Float64}:
 -0.9151044984
  0.2056320742
 -0.0002117625
 -0.2010050953

or the generator of the system

Matrix(generator(system))

2×2 Matrix{Float64}:
 -0.915104      0.205632
 -0.000211762  -0.201005

to see that the operator is slightly off, but within expectations. In a real example, this could have many reasons, e.g. noisy data, insufficient time samples or missing states.

Sticking to the same procedure as earlier, we now use a linear sparse regression to solve the problem

using ModelingToolkit

@parameters t
@variables x[1:2](t)

basis = Basis(x, x, independent_variable = t, name = :LinearBasis)

Model LinearBasis with 2 equations
States : x[1](t) x[2](t)
Independent variable: t
Equations
φ₁ = x[1](t)
φ₂ = x[2](t)

Afterwards, we simply solve the already defined problem with our Basis and a SparseOptimizer

sparse_res = solve(prob, basis, STLSQ(1e-1))
println(sparse_res)

Linear Solution with 2 equations and 3 parameters.
Returncode: solved
L₂ Norm error : [0.4193909459033963, 0.004991383862028038]
AIC : [-168.65930154106385, -1059.3084582020117]
R² : [0.9996324917308296, 0.9998961564414226]

Which holds the same equations

sparse_system = result(sparse_res)
println(sparse_system)

Model ##Basis#429 with 2 equations
States : x[1](t) x[2](t)
Parameters : p₁ p₂ p₃
Independent variable: t
Equations
Differential(t)(x[1](t)) = p₁*x[1](t) + p₂*x[2](t)
Differential(t)(x[2](t)) = p₃*x[2](t)

Again, we can have a look at the result

plot(
    plot(prob), plot(sparse_res), layout = (1,2)
)

Both results can be converted into an ODESystem

@named sys = ODESystem(
    equations(sparse_system),
    get_iv(sparse_system),
    states(sparse_system),
    parameters(sparse_system)
    );

And simulated using OrdinaryDiffEq.jl using the (known) initial conditions and the parameter mapping of the estimation.

x0 = [x[1] => u0[1], x[2] => u0[2]]
ps = parameter_map(sparse_res)

ode_prob = ODEProblem(sys, x0, tspan, ps)
estimate = solve(ode_prob, Tsit5(), saveat = prob.t);

And look at the result

plot(sol, color = :black)
plot!(estimate, color = :red, linestyle = :dash)

Copy-Pasteable Code

using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq

A = [-0.9 0.2; 0.0 -0.2]
u0 = [10.0; -10.0]
tspan = (0.0, 10.0)

f(u,p,t) = A*u

sys = ODEProblem(f, u0, tspan)
sol = solve(sys, Tsit5(), saveat = 0.05);

X = Array(sol)
t = sol.t
prob = ContinuousDataDrivenProblem(X, t)

using ModelingToolkit

@parameters t
@variables x[1:2](t)

basis = Basis(x, x, independent_variable = t, name = :LinearBasis)

sparse_res = solve(prob, basis, STLSQ(1e-1))

sparse_system = result(sparse_res)

@named sys = ODESystem(
    equations(sparse_system),
    get_iv(sparse_system),
    states(sparse_system),
    parameters(sparse_system)
    );

x0 = [x[1] => u0[1], x[2] => u0[2]]
ps = parameter_map(sparse_res)

ode_prob = ODEProblem(sys, x0, tspan, ps)
estimate = solve(ode_prob, Tsit5(), saveat = prob.t);

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

This page was generated using Literate.jl.