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 LinearAlgebra
using ModelingToolkit
using OrdinaryDiffEq
using Plots

function slow_manifold(du, u, p, t)
du = p * u
du = p * (u - u^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) 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.

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^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 u
Parameters : p₁ p₂ p₃
Independent variable: t
Equations
Differential(t)(u) = p₁*u
Differential(t)(u) = p₂*u + p₃*(u^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.

generator(system)
LinearAlgebra.Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
-1.6000000000000014
-0.7999999999999925
-0.6999999999999995
vectors:
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 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
y = -9.81sin(u) - 0.3u^3 -3.0*cos(u) - 10.0*exp(-((t-5.0)/5.0)^2)
return [x;y]
end

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])
end

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.*u);cos.(w.*u); 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]
Info

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 u
Parameters : 7
Independent variable: t
Equations
Differential(t)(u) = p₃*u
Differential(t)(u) = p₆*(u^3) + p₄*sin(u*w) + p₅*cos(u*w) + p₇*c
[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)
Warning

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.