# 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)
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`

.

`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
Equations
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`

.

`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 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]
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[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
Equations
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.