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#456 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 : [2234.5885378704893, 3814.263965569404]
R² : [1.0045382370541798, 0.994496551140509]
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#468 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.9126905858281018
-0.20190983250690941
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)
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 : [2593.706863477063, 3892.459860451417]
R² : [0.9989304796373335, 0.995114599626362]
Which holds the same equations
sparse_system = result(sparse_res)
println(sparse_system)
Model ##Basis#475 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)
);
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]
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.