Symbolic Regression of a Nonlinear System via Lifting

To infer more complex examples, EQSearch also can be called with a Basis to use predefined features. Let's look at the well-known pendulum model

using DataDrivenDiffEq
using LinearAlgebra
using OrdinaryDiffEq
using DataDrivenSR
using Plots

function pendulum!(du, u, p, t)
    du[1] = u[2]
    du[2] = -9.81 * sin(u[1])
end

u0 = [0.1, π / 2]
tspan = (0.0, 10.0)
sys = ODEProblem{true, SciMLBase.NoSpecialize}(pendulum!, u0, tspan)
sol = solve(sys, Tsit5());

We will use the data provided by our problem, but add the control signal U = sin(0.5*t) to it. Instead of using a function, like in another example

prob = DataDrivenProblem(sol)
Continuous DataDrivenProblem{Float64} ##DDProblem#405 in 2 dimensions and 43 samples

And plot the problems data.

plot(prob)
Example block output

To solve our problem, we will use EQSearch, which provides a wrapper for the symbolic regression interface. We will stick to simple operations, use a L1DistLoss, and limit the verbosity of the algorithm. Note that we do not include sin, but rather lift the search space of variables.

@variables u[1:2]
u = collect(u)

basis = Basis([polynomial_basis(u, 2); sin.(u)], u)

eqsearch_options = SymbolicRegression.Options(binary_operators = [+, *],
    loss = L1DistLoss(),
    verbosity = -1, progress = false, npop = 30,
    timeout_in_seconds = 60.0)

alg = EQSearch(eq_options = eqsearch_options)
EQSearch
  weights: Nothing nothing
  numprocs: Nothing nothing
  procs: Nothing nothing
  addprocs_function: Nothing nothing
  parallelism: Symbol serial
  runtests: Bool true
  eq_options: Options{SymbolicRegression.CoreModule.OptionsStructModule.ComplexityMapping{Int64, Int64}, OperatorEnum, Node, Expression, @NamedTuple{}, MutationWeights, false, false, nothing, Nothing, 5}

Again, we solve the problem to obtain a DataDrivenSolution with similar options as the previous example but provide a Basis along the arguments.

res = solve(prob, basis, alg, options = DataDrivenCommonOptions(maxiters = 100))
println(res)
"DataDrivenSolution{Float64}" with 2 equations and 0 parameters.
Returncode: Success
Residual sum of squares: 412.43777329559555
Note

Currently, the parameters of the result found by EQSearch are not turned into symbolic parameters. This affects some functions like dof, aicc, bic.

system = get_basis(res)
Model ##Basis#411 with 2 equations
States : u[1] u[2]
Independent variable: t
Equations
Differential(t)(u[1]) = u[2]
Differential(t)(u[2]) = -sin(u[1])

Copy-Pasteable Code

using DataDrivenDiffEq
using LinearAlgebra
using OrdinaryDiffEq
using DataDrivenSR

function pendulum!(du, u, p, t)
    du[1] = u[2]
    du[2] = -9.81 * sin(u[1])
end

u0 = [0.1, π / 2]
tspan = (0.0, 10.0)
sys = ODEProblem{true, SciMLBase.NoSpecialize}(pendulum!, u0, tspan)
sol = solve(sys, Tsit5());

prob = DataDrivenProblem(sol)

@variables u[1:2]
u = collect(u)

basis = Basis([polynomial_basis(u, 2); sin.(u)], u)

eqsearch_options = SymbolicRegression.Options(binary_operators = [+, *],
    loss = L1DistLoss(),
    verbosity = -1, progress = false, npop = 30,
    timeout_in_seconds = 60.0)

alg = EQSearch(eq_options = eqsearch_options)

res = solve(prob, basis, alg, options = DataDrivenCommonOptions(maxiters = 100))

system = get_basis(res)

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

This page was generated using Literate.jl.