Linear Time Continuous System with Controls

Now we will extend the previous example by adding some exegeneous control signals. As always, we will generate some data via OrdinaryDiffEq.jl

using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq
using Plots

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

f(u,p,t) = A*u .+ B .* sin(0.5*t)

sys = ODEProblem(f, u0, tspan)
sol = solve(sys, Tsit5(), saveat = 0.05);

We will use the data provided by our problem, but add the control signal U = sin(0.5*t) to it.

X = Array(sol)
t = sol.t
control(u,p,t) = [sin(0.5*t)]
prob = ContinuousDataDrivenProblem(X, t, U = control)

Continuous DataDrivenProblem{Float64} ##DDProblem#447 in 2 dimensions and 201 samples with controls

And plot the problems data.

plot(prob)

Again, we will use gDMD to estimate the systems dynamics. Since we have a control signal defined in the problem, the algorithm will detect it automatically and use gDMDc:

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

Linear Solution with 2 equations and 4 parameters.
Returncode: solved
L₂ Norm error : [0.7901973541577341, 0.015529506581979992]
AIC : [-39.32998241537054, -829.1676962648277]
R² : [0.9992996011479618, 0.9999315221716872]

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#459 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₄*u₁

And also plot the prediction of the recovered dynamics

plot(res)

Again, we can have a look at the generator of the system, which is independent from the inputs.

generator(system)

LinearAlgebra.Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
 -0.9126905858281008
 -0.20190983250690914
vectors:
2×2 Matrix{Float64}:
 -0.912685    0.0563158
 -0.00315083  0.193897

Sticking to the same procedure as earlier, we now use a linear sparse regression to solve the problem

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

basis = Basis([x; u], x, controls = u, independent_variable = t, name = :LinearBasis)

┌ Warning: The variable syntax (x[1:2])(t) is deprecated. Use (x(t))[1:2] instead.
│                   The former creates an array of functions, while the latter creates an array valued function.
│                   The deprecated syntax will cause an error in the next major release of Symbolics.
│                   This change will facilitate better implementation of various features of Symbolics.
└ @ Symbolics ~/.julia/packages/Symbolics/J8IHJ/src/variable.jl:129
┌ Warning: The variable syntax (u[1:1])(t) is deprecated. Use (u(t))[1:1] instead.
│                   The former creates an array of functions, while the latter creates an array valued function.
│                   The deprecated syntax will cause an error in the next major release of Symbolics.
│                   This change will facilitate better implementation of various features of Symbolics.
└ @ Symbolics ~/.julia/packages/Symbolics/J8IHJ/src/variable.jl:129
Model LinearBasis with 3 equations
States : x[1](t) x[2](t)
Independent variable: t
Equations
φ₁ = x[1](t)
φ₂ = x[2](t)
φ₃ = u[1](t)

Note that we added a new variable u[1](t) as a control to both the equations and the basis constructor. 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 4 parameters.
Returncode: solved
L₂ Norm error : [0.32342057352979164, 0.012784387671186405]
AIC : [-218.88914521865723, -868.2656437058345]
R² : [0.9997133331347746, 0.9999436268564356]

Which holds the same equations

sparse_system = result(sparse_res)
println(sparse_system)

Model ##Basis#466 with 2 equations
States : x[1](t) x[2](t)
Parameters : p₁ 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₄*u[1](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. To include the control signal, we simply substitute the control variables in the corresponding equations.

subs_control = (u[1] => sin(0.5*t))

eqs = map(equations(sparse_system)) do eq
    eq.lhs ~ substitute(eq.rhs, subs_control)
end

@named sys = ODESystem(
    eqs,
    get_iv(sparse_system),
    states(sparse_system),
    parameters(sparse_system)
    );
nothing #hide

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);
nothing #hide

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]
B = [0.0; 1.0]
u0 = [10.0; -10.0]
tspan = (0.0, 10.0)

f(u,p,t) = A*u .+ B .* sin(0.5*t)

sys = ODEProblem(f, u0, tspan)
sol = solve(sys, Tsit5(), saveat = 0.05);

X = Array(sol)
t = sol.t
control(u,p,t) = [sin(0.5*t)]
prob = ContinuousDataDrivenProblem(X, t, U = control)

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

system = result(res)

generator(system)

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

basis = Basis([x; u], x, controls = u, independent_variable = t, name = :LinearBasis)

sparse_res = solve(prob, basis, STLSQ(1e-1))

sparse_system = result(sparse_res)

subs_control = (u[1] => sin(0.5*t))

eqs = map(equations(sparse_system)) do eq
    eq.lhs ~ substitute(eq.rhs, subs_control)
end

@named sys = ODESystem(
    eqs,
    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.