Nonlinear Systems

Estimation examples for nonlinear systems.

Nonlinear System with Extended Dynamic Mode Decomposition

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)

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)

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}

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

@variables u[1:2]
Ψ = Basis([u; u[1]^2], u)
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, 73132.14114913782]
R² : [1.0, 1.0]

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

The underlying dynamics have been recovered correctly by the algorithm!

The eigendecomposition of the (generator of the) Koopman operator can be accessed via generator.

LinearAlgebra.Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
3-element Vector{Float64}:
3×3 Matrix{Float64}:
 -1.56356e-15  -1.0          -3.2209e-14
  0.613941     -1.56725e-17   1.0
 -0.789352     -8.74311e-15  -3.59306e-15

Nonlinear Systems - Sparse Identification of Nonlinear Dynamics

To find the underlying system without any Algorithms related to Koopman operator theory, we can use Sparse Identification of Nonlinear Dynamics - SINDy for short. As the name suggests, it finds the sparsest basis of functions which build the observed trajectory. Again, we will start with a nonlinear system:

using DataDrivenDiffEq
using LinearAlgebra
using ModelingToolkit
using OrdinaryDiffEq
using Plots
using Random
using Symbolics: scalarize

Random.seed!(1111) # For the noise

# Create a nonlinear pendulum
function pendulum(u, p, t)
    x = u[2]
    y = -9.81sin(u[1]) - 0.3u[2]^3 -3.0*cos(u[1]) - 10.0*exp(-((t-5.0)/5.0)^2)
    return [x;y]

u0 = [0.99π; -1.0]
tspan = (0.0, 15.0)
prob = ODEProblem(pendulum, u0, tspan)
sol = solve(prob, Tsit5(), saveat = 0.01)

# Create the data with additional noise
X = sol[:,:] + 0.1 .* randn(size(sol))
DX = similar(sol[:,:])

for (i, xi) in enumerate(eachcol(sol[:,:]))
    DX[:,i] = pendulum(xi, [], sol.t[i])

ts = sol.t

To estimate the system, we first create a DataDrivenProblem via feeding in the measurement data. Using a Collocation method, it automatically provides the derivative. Control signals can be passed in as a function (u,p,t)->control or an array of measurements.

prob = ContinuousDataDrivenProblem(X, ts, GaussianKernel() ,
    U = (u,p,t)->[exp(-((t-5.0)/5.0)^2)], p = ones(2))

Now we infer the system structure. First we define a Basis which collects all possible candidate terms. Since we want to use SINDy, we call solve with an Optimizer, in this case STLSQ which iterates different sparsity thresholds and returns a pareto optimal solution of the underlying sparse_regression!. Note that we include the control signal in the basis as an additional variable c.

@variables u[1:2] c[1:1]
@parameters w[1:2]
u = scalarize(u)
c = scalarize(c)
w = scalarize(w)

h = Num[sin.(w[1].*u[1]);cos.(w[2].*u[1]); polynomial_basis(u, 5); c]

basis = Basis(h, u, parameters = w, controls = c)

λs = exp10.(-10:0.1:-1)
opt = STLSQ(λs)
res = solve(prob, basis, opt, progress = false, denoise = false, normalize = false, maxiter = 5000)
Linear Solution with 2 equations and 7 parameters.
Returncode: solved
L₂ Norm error : [24.098805956305238, 53.763768116208304]
AIC : [12417.437062312865, 10008.51971337545]
R² : [0.9959910957032931, 0.996294784415434]

A more detailed description of the result can be printed via print(res, Val{true}), which also includes the discovered equations and parameter values.

Where the resulting DataDrivenSolution stores information about the inferred model and the parameters:

system = result(res);
params = parameters(res);
Model ##Basis#393 with 2 equations
States : u[1] u[2]
Parameters : 7
Independent variable: t
Differential(t)(u[1]) = p₃*u[2]
Differential(t)(u[2]) = p₆*(u[2]^3) + p₄*sin(u[1]*w[1]) + p₅*cos(u[1]*w[2]) + p₇*c[1]
[1.0, 1.0, 1.01, -9.58, -3.0, -0.31, -9.85]

Since any system obtained via a solve command is a Basis and hence a subtype of an AbstractSystem defined in ModelingToolkit, we can simply simulate the result via:

infered_prob = ODEProblem(system, u0, tspan, parameters(res))
infered_solution = solve(infered_prob, Tsit5(), saveat = ts)

As of now, the control input is dropped in the simulation of a system. We are working on this and pull requests are welcome!

As we can see above, the estimated system matches the ground truth reasonably well.