Nonlinear Time Continuous System

Similarly, we can use the Extended Dynamic Mode Decomposition via a nonlinear Basis of observables. Here, we will look at a rather famous example with a finite dimensional solution.

using DataDrivenDiffEq
using LinearAlgebra
using ModelingToolkit
using OrdinaryDiffEq
using Plots

function slow_manifold(du, u, p, t)
    du[1] = p[1] * u[1]
    du[2] = p[2] * (u[2] - u[1]^2)
end

u0 = [3.0; -2.0]
tspan = (0.0, 5.0)
p = [-0.8; -0.7]

problem = ODEProblem(slow_manifold, u0, tspan, p)
solution = solve(problem, Tsit5(), saveat = 0.01);
plot(solution)

Since we are dealing with a continuous system in time, we define the associated DataDrivenProblem accordingly using the measured states X, their derivatives DX and the time t.

prob = ContinuousDataDrivenProblem(solution)

Continuous DataDrivenProblem{Float64} ##DDProblem#466 in 2 dimensions and 501 samples

Additionally, we need to define the Basis for our lifting, before we solve the problem in the lifted space.

@parameters t
@variables u[1:2](t)
Ψ = Basis([u; u[1]^2], u, independent_variable = t)
res = solve(prob, Ψ, DMDPINV(), digits = 1)
system = result(res)

Linear Solution with 2 equations and 3 parameters.
Returncode: solved
L₂ Norm error : [0.0, 1.0114581578656041e-29]
AIC : [-Inf, -33442.5509179254]
R² : [1.0, 1.0]

Model ##Koopman#475 with 2 equations
States : u[1](t) u[2](t)
Parameters : p₁ p₂ p₃
Independent variable: t
Equations
Differential(t)(u[1](t)) = p₁*u[1](t)
Differential(t)(u[2](t)) = p₂*u[2](t) + p₃*(u[1](t)^2)
[-0.8, -0.7, 0.7]

The underlying dynamics have been recovered correctly by the algorithm! Similarly we could use sparse identification to solve the problem

sparse_res = solve(prob, Ψ, STLSQ(), digits = 1)
println(sparse_res)

Linear Solution with 2 equations and 3 parameters.
Returncode: solved
L₂ Norm error : [0.0, 1.0114581578656041e-29]
AIC : [-Inf, -33442.5509179254]
R² : [1.0, 1.0]

And the resulting system

sparse_system = result(sparse_res)
println(sparse_system)

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

We can also directly look at the parameters of each result

parameter_map(res)

3-element Vector{Pair{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}, Float64}}:
 p₁ => -0.8
 p₂ => -0.7
 p₃ => 0.7

Note that we are using parameter_map instead of just parameters, which returns a vector suitable to use with ModelingToolkit.

parameter_map(sparse_res)

3-element Vector{Pair{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}, Float64}}:
 p₁ => -0.8
 p₂ => -0.7
 p₃ => 0.7

To simulate the system, we create an ODESystem from the result

Both results can be converted into an ODESystem

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


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

3-element Vector{Pair{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}, Float64}}:
 p₁ => -0.8
 p₂ => -0.7
 p₃ => 0.7

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

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

And look at the result

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

Copy-Pasteable Code

using DataDrivenDiffEq
using LinearAlgebra
using ModelingToolkit
using OrdinaryDiffEq

function slow_manifold(du, u, p, t)
    du[1] = p[1] * u[1]
    du[2] = p[2] * (u[2] - u[1]^2)
end

u0 = [3.0; -2.0]
tspan = (0.0, 5.0)
p = [-0.8; -0.7]

problem = ODEProblem(slow_manifold, u0, tspan, p)
solution = solve(problem, Tsit5(), saveat = 0.01);

prob = ContinuousDataDrivenProblem(solution)

@parameters t
@variables u[1:2](t)
Ψ = Basis([u; u[1]^2], u, independent_variable = t)
res = solve(prob, Ψ, DMDPINV(), digits = 1)
system = result(res)

sparse_res = solve(prob, Ψ, STLSQ(), digits = 1)

sparse_system = result(sparse_res)

parameter_map(res)

parameter_map(sparse_res)

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


x0 = [u[1] => u0[1], u[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.