Implicit Sparse Identification of Nonlinear Dynamics

While SINDy works well for ODEs, some systems take the form of DAEs. A common form is f(x, p, t) - g(x, p, t)*dx = 0. These can be inferred via ISINDy, which extends SINDy for Implicit problems. In particular, it solves

where $\Xi$ lies in the nullspace of $\Theta$.

Example : Michaelis-Menten Kinetics

Let's try to infer the Michaelis-Menten Kinetics, like in the corresponding paper. We start by generating the corresponding data.

using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra
using Plots
gr()

function michaelis_menten(u, p, t)
    [0.6 - 1.5u[1]/(0.3+u[1])]
end

u0 = [0.5]
tspan = (0.0, 5.0)
problem = ODEProblem(michaelis_menten, u0, tspan)

solution = solve(problem, Tsit5(), saveat = 0.1, atol = 1e-7, rtol = 1e-7)

savefig("iSINDy_example.png")

X = solution[:,:]
DX = similar(X)
for (i, xi) in enumerate(eachcol(X))
    DX[:, i] = michaelis_menten(xi, [], 0.0)
end

@variables u
basis= Basis([u^i for i in 0:4], [u])

##Basis#374 : 5 dimensional basis in ["u"]

The signature of ISINDy is equal to SINDy, but requires an AbstractSubspaceOptimizer. Currently, DataDrivenDiffEq just implements ADM() based on alternating directions. rtol gets passed into the derivation of the nullspace via LinearAlgebra.

opt = ADM(1.1e-1)

ADM{Float64,ProximalOperators.NormL1{Float64}}(0.11, description : weighted L1 norm
domain      : AbstractArray{Real}, AbstractArray{Complex}
expression  : x ↦ λ||x||_1
parameters  : λ = 0.11)

Since ADM() returns sparsified columns of the nullspace we need to find a pareto optimal solution. To achieve this, we provide a sufficient cost function g to ISINDy. This allows us to evaluate each individual column of the sparse matrix on its 0-norm (sparsity) and the 2-norm of the matrix vector product of $\Theta^T \xi$ (nullspace). This is a default setting which can be changed by providing a function f which maps the coefficients and the library onto a feature space. Here, we want to set the focus on the the magnitude of the deviation from the nullspace.

Ψ = ISINDy(X, DX, basis, opt, g = x->norm([1e-1*x[1]; x[2]]), maxiter = 100)

The function call returns a SparseIdentificationResult. As in Sparse Identification of Nonlinear Dynamics, we can transform the SparseIdentificationResult into an ODESystem.

# Transform into ODE System
sys = ODESystem(Ψ)
dudt = ODEFunction(sys)
ps = parameters(Ψ)

estimator = ODEProblem(dudt, u0, tspan, ps)
estimation = solve(estimator, Tsit5(), saveat = 0.1)

The model recovered by ISINDy is correct

print_equations(Ψ)

##Basis#382 : 1 dimensional basis in ["u"]
Parameters : SymbolicUtils.Sym{ModelingToolkit.Parameter{Real}}[p₁, p₂, p₃, p₄]

Independent variable: t
Equations
φ₁ = (-((p₂ * (u ^ 0)) + (p₄ * (u ^ 1)))) / ((p₁ * (u ^ 0)) + (p₃ * (u ^ 1)))

The parameters are off a little, but, as before, we can use DiffEqFlux to tune them.

Example : Cart-Pole with Time-Dependent Control

Implicit dynamics can also be reformulated as an explicit problem as stated in this paper. The algorithm searches the correct equations by trying out all candidate functions as a right hand side and performing a sparse regression onto the remaining set of candidates. Let's start by defining the problem and generate the data:

using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra
using Plots
gr()

function cart_pole(u, p, t)
    du = similar(u)
    F = -0.2 + 0.5*sin(6*t) # the input
    du[1] = u[3]
    du[2] = u[4]
    du[3] = -(19.62*sin(u[1])+sin(u[1])*cos(u[1])*u[3]^2+F*cos(u[1]))/(2-cos(u[1])^2)
    du[4] = -(sin(u[1])*u[3]^2 + 9.81*sin(u[1])*cos(u[1])+F)/(2-cos(u[1])^2)
    return du
end

u0 = [0.3; 0; 1.0; 0]
tspan = (0.0, 16.0)
dt = 0.001
cart_pole_prob = ODEProblem(cart_pole, u0, tspan)
solution = solve(cart_pole_prob, Tsit5(), saveat = dt)

# Create the differential data
X = solution[:,:]
DX = similar(X)
for (i, xi) in enumerate(eachcol(X))
    DX[:, i] = cart_pole(xi, [], solution.t[i])
end

We see that we include a forcing term F inside the model which is depending on t. As before, we will also need a Basis to derive our equations from:

@variables u[1:4] t
polys = Any[]
for i ∈ 0:4
    if i == 0
        push!(polys, u[1]^0)
    else
        if i < 2
            push!(polys, u.^i...)
        else
            push!(polys, u[3:4].^i...)
        end

    end
end
push!(polys, sin.(u[1])...)
push!(polys, cos.(u[1]))
push!(polys, sin.(u[1]).*u[3:4]...)
push!(polys, sin.(u[1]).*u[3:4].^2...)
push!(polys, cos.(u[1]).^2...)
push!(polys, sin.(u[1]).*cos.(u[1])...)
push!(polys, sin.(u[1]).*cos.(u[1]).*u[3:4]...)
push!(polys, sin.(u[1]).*cos.(u[1]).*u[3:4].^2...)
push!(polys, -0.2+0.5*sin(6*t))
push!(polys, (-0.2+0.5*sin(6*t))*cos(u[1]))
push!(polys, (-0.2+0.5*sin(6*t))*sin(u[1]))
basis= Basis(polys, u, iv = t)

##Basis#400 : 26 dimensional basis in ["u₁", "u₂", "u₃", "u₄"]

We added the time dependent input directly into the basis to account for its influence.

NOTE : Including input signals may change in future releases!

Like for a SINDy, we can use any AbstractOptimizer with a pareto front optimization over different thresholds.

λ = exp10.(-4:0.1:-1)
g(x) = norm([1e-3; 10.0] .* x, 2)
Ψ = ISINDy(X[:,:], DX[:, :], basis, λ, STRRidge(), maxiter = 100, normalize = false, t = solution.t, g = g)

# Transform into ODE System
sys = ODESystem(Ψ)
dudt = ODEFunction(sys)
ps = parameters(Ψ)

# Simulate
estimator = ODEProblem(dudt, u0, tspan, ps)
sol_ = solve(estimator, Tsit5(), saveat = dt)

Let's have a look at the equations recovered. They match up.

print_equations(Ψ)

##Basis#414 : 4 dimensional basis in ["u₁", "u₂", "u₃", "u₄"]
Parameters : SymbolicUtils.Sym{ModelingToolkit.Parameter{Real}}[p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈, p₉, p₁₀, p₁₁, p₁₂, p₁₃, p₁₄]

Independent variable: t
Equations
φ₁ = (-(p₂ * (u₃ ^ 1))) / (p₁ * (u₁ ^ 0))
φ₂ = (-(p₄ * (u₄ ^ 1))) / (p₃ * (u₁ ^ 0))
φ₃ = (-(((p₆ * sin(u₁)) + (p₈ * ((sin(u₁) * cos(u₁)) * (u₃ ^ 2)))) + (p₉ * ((-0.2 + (0.5 * sin(6 * t))) * cos(u₁))))) / ((p₅ * (u₁ ^ 0)) + (p₇ * (cos(u₁) ^ 2)))
φ₄ = (-(((p₁₁ * (sin(u₁) * (u₃ ^ 2))) + (p₁₃ * (sin(u₁) * cos(u₁)))) + (p₁₄ * (-0.2 + (0.5 * sin(6 * t)))))) / ((p₁₀ * (u₁ ^ 0)) + (p₁₂ * (cos(u₁) ^ 2)))

Alternatively, we can also use the input as an extended state x.

@variables u[1:4] t x
polys = Any[]
# Lots of basis functions -> sindy pi can handle more than ADM()
for i ∈ 0:4
    if i == 0
        push!(polys, u[1]^0)
    else
        if i < 2
            push!(polys, u.^i...)
        else
            push!(polys, u[3:4].^i...)
        end

    end
end
push!(polys, sin.(u[1])...)
push!(polys, cos.(u[1]))
push!(polys, sin.(u[1]).*u[3:4]...)
push!(polys, sin.(u[1]).*u[3:4].^2...)
push!(polys, cos.(u[1]).^2...)
push!(polys, sin.(u[1]).*cos.(u[1])...)
push!(polys, sin.(u[1]).*cos.(u[1]).*u[3:4]...)
push!(polys, sin.(u[1]).*cos.(u[1]).*u[3:4].^2...)
push!(polys, x)
push!(polys, x*cos(u[1]))
push!(polys, x*sin(u[1]))
basis= Basis(polys, vcat(u, x), iv = t)

##Basis#431 : 26 dimensional basis in ["u₁", "u₂", "u₃", "u₄", "x"]

Now we include the input signal into the extended state array Xᵤ and perform a sparse regression.

U = -0.2 .+ 0.5*sin.(6*solution.t)
Xᵤ = vcat(X, U')

λ = exp10.(-4:0.5:-1)
g(x) = norm([1e-3; 10.0] .* x, 2)
Ψ = ISINDy(Xᵤ[:,:], DX[:, :], basis, λ, STRRidge(), maxiter = 100, normalize = false, t = solution.t, g = g)
print_equations(Ψ, show_parameter = true)

f_1 = (-(-1.0 * (u₃ ^ 1))) / (1.0 * (u₁ ^ 0))
f_2 = (-(-0.9999999999999999 * (u₄ ^ 1))) / (1.0 * (u₁ ^ 0))
f_3 = (-(((0.9999999999999999 * sin(u₁)) + (0.05096839959225282 * ((sin(u₁) * cos(u₁)) * (u₃ ^ 2)))) + (0.0509683995922528 * (x * cos(u₁))))) / ((0.10193679918450561 * (u₁ ^ 0)) + (-0.05096839959225284 * (cos(u₁) ^ 2)))
f_4 = (-(((0.10193679918450552 * (sin(u₁) * (u₃ ^ 2))) + (1.0 * (sin(u₁) * cos(u₁)))) + (0.10193679918450561 * x))) / ((0.2038735983690113 * (u₁ ^ 0)) + (-0.10193679918450563 * (cos(u₁) ^ 2)))

Currently, we can not generate an ODESystem out of the resulting equations, which is a work in progress.

Functions

ISINDy(X, Y, Ψ, opt = ADM(); f, g, maxiter, rtol, p, t, convergence_error)
ISINDy(X, Y, Ψ, opt; f, g, maxiter, rtol, p, t, convergence_error, normalize, denoise)
ISINDy(X, Y, Ψ, lamdas, opt; f, g, maxiter, rtol, p, t, convergence_error, normalize, denoise)