Linear Time Continuous System
Similar to the linear time discrete example
, we will now estimate a linear time continuous system $\partial_t u = A u$. We simulate the correspoding system using OrdinaryDiffEq.jl
and generate a ContinuousDataDrivenProblem
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, 10.0)
f(u,p,t) = A*u
sys = ODEProblem(f, u0, tspan)
sol = solve(sys, Tsit5(), saveat = 0.05);
We could use the DESolution
to define our problem, but here we want to use the data for didactic purposes. For a ContinuousDataDrivenProblem
, we need either the state trajectory and the timepoints or the state trajectory and its derivate.
X = Array(sol)
t = sol.t
prob = ContinuousDataDrivenProblem(X, t)
Continuous DataDrivenProblem{Float64} ##DDProblem#429 in 2 dimensions and 201 samples
And plot the problems data.
plot(prob)
We can see that the derivative has been automatically added via a collocation
method, which defaults to a LinearInterpolation
. We can do a visual check and compare our derivatives with the interpolation of the ODESolution
.
DX = Array(sol(t, Val{1}))
scatter(t, DX', label = ["Solution" nothing], color = :red, legend = :bottomright)
plot!(t, prob.DX', label = ["Linear Interpolation" nothing], color = :black)
Since we have a linear system, we can use gDMD
, which approximates the generator of the dynamics
res = solve(prob, DMDSVD())
println(res)
Linear Solution with 2 equations and 4 parameters.
Returncode: solved
L₂ Norm error : [0.23746883585660647, 0.0003580043828573712]
AIC : [2717.891563428103, 5329.784635305852]
R² : [0.9997232799213416, 1.0002549921804418]
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)
println(system)
Model ##Koopman#441 with 2 equations
States : x₁ x₂
Parameters : p₁ p₂ p₃ p₄
Independent variable: t
Equations
Differential(t)(x₁) = p₁*x₁ + p₂*x₂
Differential(t)(x₂) = p₃*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 => [2717.89, 5329.78]
:R² => [0.999723, 1.00025]
:L₂ => [0.237469, 0.000358004]
And check the parameters of the result
parameters(res)
4-element Vector{Float64}:
-0.9151044984
0.2056320742
-0.0002117625
-0.2010050953
or the generator of the system
Matrix(generator(system))
2×2 Matrix{Float64}:
-0.915104 0.205632
-0.000211762 -0.201005
to see that the operator is slightly off, but within expectations. In a real example, this could have many reasons, e.g. noisy data, insufficient time samples or missing states.
Sticking to the same procedure as earlier, we now use a linear sparse regression to solve the problem
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(1e-1))
println(sparse_res)
Linear Solution with 2 equations and 3 parameters.
Returncode: solved
L₂ Norm error : [0.4193909459034067, 0.004991383862028038]
AIC : [2487.2471761218662, 4268.545489443772]
R² : [0.9969812390823322, 0.9823366808911742]
Which holds the same equations
sparse_system = result(sparse_res)
println(sparse_system)
Model ##Basis#448 with 2 equations
States : x[1](t) x[2](t)
Parameters : p₁ p₂ p₃
Independent variable: t
Equations
Differential(t)(x[1](t)) = p₁*x[1](t) + p₂*x[2](t)
Differential(t)(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 an ODESystem
@named sys = ODESystem(
equations(sparse_system),
get_iv(sparse_system),
states(sparse_system),
parameters(sparse_system)
);
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)
ode_prob = ODEProblem(sys, x0, tspan, ps)
estimate = solve(ode_prob, Tsit5(), saveat = prob.t);
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, 10.0)
f(u,p,t) = A*u
sys = ODEProblem(f, u0, tspan)
sol = solve(sys, Tsit5(), saveat = 0.05);
X = Array(sol)
t = sol.t
prob = ContinuousDataDrivenProblem(X, t)
using ModelingToolkit
@parameters t
@variables x[1:2](t)
basis = Basis(x, x, independent_variable = t, name = :LinearBasis)
sparse_res = solve(prob, basis, STLSQ(1e-1))
sparse_system = result(sparse_res)
@named sys = ODESystem(
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)
ode_prob = ODEProblem(sys, x0, tspan, ps)
estimate = solve(ode_prob, Tsit5(), saveat = prob.t);
# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl
This page was generated using Literate.jl.