Linear Time Discrete System

We will start by estimating the underlying dynamical system of a time discrete process based on some measurements via Dynamic Mode Decomposition on a simple linear system of the form $u(k+1) = A u(k)$.

At first, we simulate the correspoding system using OrdinaryDiffEq.jl and generate a DiscreteDataDrivenProblem from the simulated data.

using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq
using Plots

A = [0.9 -0.2; 0.0 0.2]
u0 = [10.0; -10.0]
tspan = (0.0, 11.0)

f(u,p,t) = A*u

sys = DiscreteProblem(f, u0, tspan)
sol = solve(sys, FunctionMap());

Next we transform our simulated solution into a DataDrivenProblem. Given that the solution knows its a discrete solution, we can simply write

prob = DataDrivenProblem(sol)

Discrete DataDrivenProblem{Float64} ##DDProblem#387 in 2 dimensions and 12 samples

And plot the solution and the problem

plot(sol, label = string.([:x₁ :x₂]))
scatter!(prob)

To estimate the underlying operator in the states $x_1, x_2$, we solve the estimation problem using the DMDSVD algorithm for approximating the operator. First, we will have a look at the DataDrivenSolution

res = solve(prob, DMDSVD(), digits = 1)

Linear Solution with 2 equations and 3 parameters.
Returncode: solved
L₂ Norm error : [0.0, 0.0]
AIC : [-Inf, -Inf]
R² : [1.0, 1.0]

We see that the system has been recovered correctly, indicated by the small error and high AIC score of the result. We can confirm this by looking at the resulting Basis

system = result(res)
using Symbolics

Model ##Koopman#395 with 2 equations
States : x₁ x₂
Parameters : p₁ p₂ p₃
Independent variable: t
Equations
Difference(t; dt=5.5, update=false)(x₁) = p₁*x₁ + p₂*x₂
Difference(t; dt=5.5, update=false)(x₂) = p₃*x₂

And also plot the prediction of the recovered dynamics

plot(res)

Or a have a look at the metrics of the result

metrics(res)

Dict{Any, Any} with 3 entries:
  :AIC => [-Inf, -Inf]
  :R²  => [1.0, 1.0]
  :L₂  => [0.0, 0.0]

To have a look at the representation of the operator as a Matrix, we can simply call

Matrix(system)

2×2 Matrix{Float64}:
  0.9          -0.2
 -1.07938e-17   0.2

to see that the operator is indeed our initial A. Since we have a linear representation, we can gain further insights into the stability of the dynamics via its eigenvalues

eigvals(system)

2-element Vector{Float64}:
 0.19999999999999987
 0.9000000000000004

And plot the stability margin of the discrete System

φ = 0:0.01π:2π
plot(sin.(φ), cos.(φ), xlabel = "Real", ylabel = "Im", label = "Stability margin", color = :red, linestyle = :dash)
scatter!(real(eigvals(system)), imag(eigvals(system)), label = "Eigenvalues", color = :black, marker = :cross)

Similarly, we could use a sparse regression to derive our system from our data. We start by defining a Basis

using ModelingToolkit

@parameters t
@variables x[1:2](t)

basis = Basis(x, x, independent_variable = t, name = :LinearBasis)

Model LinearBasis with 2 equations
States : x[1](t) x[2](t)
Independent variable: t
Equations
φ₁ = x[1](t)
φ₂ = x[2](t)

Afterwards, we simply solve the already defined problem with our Basis and a SparseOptimizer

sparse_res = solve(prob, basis, STLSQ())

Linear Solution with 2 equations and 3 parameters.
Returncode: solved
L₂ Norm error : [0.0, 0.0]
AIC : [-Inf, -Inf]
R² : [1.0, 1.0]

Which holds the same equations

sparse_system = result(sparse_res)

Model ##Basis#403 with 2 equations
States : x[1](t) x[2](t)
Parameters : p₁ p₂ p₃
Independent variable: t
Equations
Difference(t; dt=5.5, update=false)(x[1](t)) = p₁*x[1](t) + p₂*x[2](t)
Difference(t; dt=5.5, update=false)(x[2](t)) = p₃*x[2](t)

Again, we can have a look at the result

plot(
    plot(prob), plot(sparse_res), layout = (1,2)
)

Both results can be converted into a DiscreteProblem

@named sys = DiscreteSystem(equations(sparse_system), get_iv(sparse_system),states(sparse_system), parameters(sparse_system))

ModelingToolkit.DiscreteSystem(Symbolics.Equation[Difference(t; dt=5.5, update=false)(x[1](t)) ~ p₁*x[1](t) + p₂*x[2](t), Difference(t; dt=5.5, update=false)(x[2](t)) ~ p₃*x[2](t)], t, SymbolicUtils.Term{Real, Base.ImmutableDict{DataType, Any}}[x[1](t), x[2](t)], SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}[p₁, p₂, p₃], Dict{Any, Any}(:p₁ => p₁, :p₃ => p₃, :p₂ => p₂, :x => map(Symbolics.CallWith((t,)), x)), Any[], Symbolics.Equation[], :sys, ModelingToolkit.DiscreteSystem[], Dict{Any, Any}(), nothing, nothing, nothing, nothing)

And simulated using OrdinaryDiffEq.jl using the (known) initial conditions and the parameter mapping of the estimation.

x0 = [x[1] => u0[1], x[2] => u0[2]]
ps = parameter_map(sparse_res)

discrete_prob = DiscreteProblem(sys, x0, tspan, ps)
estimate = solve(discrete_prob, FunctionMap());

And look at the result

plot(sol, color = :black)
plot!(estimate, color = :red, linestyle = :dash)

Copy-Pasteable Code

using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq

A = [0.9 -0.2; 0.0 0.2]
u0 = [10.0; -10.0]
tspan = (0.0, 11.0)

f(u,p,t) = A*u

sys = DiscreteProblem(f, u0, tspan)
sol = solve(sys, FunctionMap());

prob = DataDrivenProblem(sol)

res = solve(prob, DMDSVD(), digits = 1)

system = result(res)
using Symbolics

using ModelingToolkit

@parameters t
@variables x[1:2](t)

basis = Basis(x, x, independent_variable = t, name = :LinearBasis)

sparse_res = solve(prob, basis, STLSQ())

sparse_system = result(sparse_res)

@named sys = DiscreteSystem(equations(sparse_system), get_iv(sparse_system),states(sparse_system), parameters(sparse_system))

x0 = [x[1] => u0[1], x[2] => u0[2]]
ps = parameter_map(sparse_res)

discrete_prob = DiscreteProblem(sys, x0, tspan, ps)
estimate = solve(discrete_prob, FunctionMap());

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

This page was generated using Literate.jl.