From 0f27f9b98624b65802186ae3ffbda34bb1732b90 Mon Sep 17 00:00:00 2001 From: AlCap23 Date: Tue, 19 Nov 2019 14:55:54 +0100 Subject: [PATCH 1/4] Added DMDc method --- src/DataDrivenDiffEq.jl | 4 +++ src/dmdc.jl | 72 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+) create mode 100644 src/dmdc.jl diff --git a/src/DataDrivenDiffEq.jl b/src/DataDrivenDiffEq.jl index 8a106c67f..7e54e7fa6 100644 --- a/src/DataDrivenDiffEq.jl +++ b/src/DataDrivenDiffEq.jl @@ -23,6 +23,10 @@ export ExtendedDMD export dynamics, linear_dynamics export reduce_basis, update! +include("./dmdc.jl") +export DMDc +export get_dynamics, get_input_map + include("./sindy.jl") export SInDy diff --git a/src/dmdc.jl b/src/dmdc.jl new file mode 100644 index 000000000..27ff6ab3e --- /dev/null +++ b/src/dmdc.jl @@ -0,0 +1,72 @@ +using DifferentialEquations +using Plots +using DataDrivenDiffEq +using LinearAlgebra + +mutable struct DMDc{K, B, Q, P} + koopman::K + forcing::B + + Qₖ::Q + Pₖ::P +end + + +function DMDc(X::AbstractArray, Y::AbstractArray, Γ::AbstractArray; Δt::Float64 = 0.0) + @assert all(size(X) .== size(Y)) + @assert size(X)[2] .== size(Γ)[2] + + + nₓ = size(X)[1] + nᵤ = size(Γ)[1] + + Ω = vcat(X, Γ) + G = X * pinv(Ω) + + Ã = G[:, 1:nₓ] + B̃ = G[:, nₓ+1:end] + + # Eigen Decomposition for solution + Λ, W = eigen(Ã) + + if Δt > 0.0 + # Casting Complex enforces results + ω = log.(Complex.(Λ)) / Δt + else + ω = [] + end + + koopman = ExactDMD(Ã, Λ, ω, W, nothing, nothing) + + + return DMDc(koopman, B̃, nothing, nothing) +end + +function DMDc(X::AbstractArray, Γ::AbstractArray; Δt::Float64 = 0.0) + @assert size(X)[2]-1 == size(Γ)[2] "Provide consistent input data." + return DMDc(X[:, 1:end-1], X[:, 2:end], Γ, Δt = Δt) +end + +# Some nice functions +LinearAlgebra.eigen(m::DMDc) = eigen(m.koopman) +LinearAlgebra.eigvals(m::DMDc) = eigvals(m.koopman) +LinearAlgebra.eigvecs(m::DMDc) = eigvecs(m.koopman) + +get_dynamics(m::DMDc) = m.koopman +get_input_map(m::DMDc) = m.forcing + +# TODO this can be done better, maybe use macros +function DataDrivenDiffEq.dynamics(m::DMDc; discrete::Bool = true) + if discrete + nᵢ = size(m.forcing)[2] + function zero_callback(u, p, t) + return zeros(eltype(m.forcing), nᵢ) + end + @inline function dudt_(du, u, p, t; y = zero_callback) + du .= m.koopman.Ã * u + m.forcing * y(u, p, t) + end + return dudt_ + else + throw(ErrorException("Continouos dynamics are not implemented right now.")) + end +end From 0ed83a2265a3bb2345a4c5daa6424c5fa7be04d8 Mon Sep 17 00:00:00 2001 From: AlCap23 Date: Mon, 25 Nov 2019 16:08:24 +0100 Subject: [PATCH 2/4] Better interface --- src/dmdc.jl | 57 ++++++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 45 insertions(+), 12 deletions(-) diff --git a/src/dmdc.jl b/src/dmdc.jl index 27ff6ab3e..37667db36 100644 --- a/src/dmdc.jl +++ b/src/dmdc.jl @@ -1,7 +1,7 @@ -using DifferentialEquations -using Plots +using OrdinaryDiffEq using DataDrivenDiffEq using LinearAlgebra +using Plots mutable struct DMDc{K, B, Q, P} koopman::K @@ -12,7 +12,7 @@ mutable struct DMDc{K, B, Q, P} end -function DMDc(X::AbstractArray, Y::AbstractArray, Γ::AbstractArray; Δt::Float64 = 0.0) +function DMDc(X::AbstractArray, Y::AbstractArray, Γ::AbstractArray; B::AbstractArray = [], Δt::Float64 = 0.0) @assert all(size(X) .== size(Y)) @assert size(X)[2] .== size(Γ)[2] @@ -20,11 +20,16 @@ function DMDc(X::AbstractArray, Y::AbstractArray, Γ::AbstractArray; Δt::Float6 nₓ = size(X)[1] nᵤ = size(Γ)[1] - Ω = vcat(X, Γ) - G = X * pinv(Ω) + if isempty(B) + Ω = vcat(X, Γ) + G = Y * pinv(Ω) - Ã = G[:, 1:nₓ] - B̃ = G[:, nₓ+1:end] + Ã = G[:, 1:nₓ] + B̃ = G[:, nₓ+1:end] + else + Ã = (Y - B*Γ)*pinv(X) + B̃ = B + end # Eigen Decomposition for solution Λ, W = eigen(Ã) @@ -42,9 +47,10 @@ function DMDc(X::AbstractArray, Y::AbstractArray, Γ::AbstractArray; Δt::Float6 return DMDc(koopman, B̃, nothing, nothing) end -function DMDc(X::AbstractArray, Γ::AbstractArray; Δt::Float64 = 0.0) + +function DMDc(X::AbstractArray, Γ::AbstractArray; B::AbstractArray = [], Δt::Float64 = 0.0) @assert size(X)[2]-1 == size(Γ)[2] "Provide consistent input data." - return DMDc(X[:, 1:end-1], X[:, 2:end], Γ, Δt = Δt) + return DMDc(X[:, 1:end-1], X[:, 2:end], Γ, B = B, Δt = Δt) end # Some nice functions @@ -52,21 +58,48 @@ LinearAlgebra.eigen(m::DMDc) = eigen(m.koopman) LinearAlgebra.eigvals(m::DMDc) = eigvals(m.koopman) LinearAlgebra.eigvecs(m::DMDc) = eigvecs(m.koopman) +DataDrivenDiffEq.isstable(m::DMDc) = isstable(m.koopman) + get_dynamics(m::DMDc) = m.koopman get_input_map(m::DMDc) = m.forcing # TODO this can be done better, maybe use macros function DataDrivenDiffEq.dynamics(m::DMDc; discrete::Bool = true) if discrete - nᵢ = size(m.forcing)[2] + dims = size(m.forcing) + nᵢ = length(dims)<=1 ? 1 : dims[2] + println(zeros(eltype(m.forcing), nᵢ)) function zero_callback(u, p, t) return zeros(eltype(m.forcing), nᵢ) end - @inline function dudt_(du, u, p, t; y = zero_callback) - du .= m.koopman.Ã * u + m.forcing * y(u, p, t) + @inline function dudt_(u, p, t; y = zero_callback) + m.koopman.Ã * u + m.forcing .* y(u, p, t) end return dudt_ else throw(ErrorException("Continouos dynamics are not implemented right now.")) end end + +X = [4 2 1 0.5 0.25; 7 0.7 0.07 0.007 0.0007] +U = [-4 -2 -1 -0.5] +B = Float32[1; 0] + +sys = DMDc(X, U) +sys = DMDc(X, U, B = B) + +isstable(sys) + +get_input_map(sys) + +get_dynamics(sys) + +eigen(sys) + +dudt_ = dynamics(sys) +dudt_(X[:, 1], [], 0.0) + +prob = DiscreteProblem(dudt_, X[:, 1], (0.0, 10.0)) +sol = solve(prob) + +plot(sol) From 3d491079bf894bd1d2ec3941dd57a537209d1f39 Mon Sep 17 00:00:00 2001 From: Julius Martensen Date: Fri, 29 Nov 2019 21:55:13 +0100 Subject: [PATCH 3/4] Works now. Add more nice function and think about control input --- examples/DMDc_Examples.jl | 31 +++++++++++++++++++++++ src/DataDrivenDiffEq.jl | 3 ++- src/dmdc.jl | 52 +++++++-------------------------------- 3 files changed, 42 insertions(+), 44 deletions(-) create mode 100644 examples/DMDc_Examples.jl diff --git a/examples/DMDc_Examples.jl b/examples/DMDc_Examples.jl new file mode 100644 index 000000000..efcc4c7fe --- /dev/null +++ b/examples/DMDc_Examples.jl @@ -0,0 +1,31 @@ +using DataDrivenDiffEq +using DifferentialEquations +using Plots +gr() + +# Define measurements from unstable system with known control input +X = [4 2 1 0.5 0.25; 7 0.7 0.07 0.007 0.0007] +U = [-4 -2 -1 -0.5] +B = Float32[1; 0] + +# See fail with unknown input +sys = DMDc(X, U) +# But with a little more knowledge +sys = DMDc(X, U, B = B) + +# Extract the DMD from inside DMDc +get_dynamics(sys) +# Acess all the other stuff +eigen(sys) .≈ eigen(get_dynamics(sys)) +# Get unforced dynamics +dudt_ = dynamics(sys) +prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.)) +sol_unforced = solve(prob) +plot(sol_unforced) + +# Create a system with cos control input to stabilize +dudt_ = dynamics(sys, control = (u, p, t) -> -0.5u[1]) +prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.)) +sol = solve(prob) + +plot!(sol) diff --git a/src/DataDrivenDiffEq.jl b/src/DataDrivenDiffEq.jl index 2f9e2b96d..4797d6e0b 100644 --- a/src/DataDrivenDiffEq.jl +++ b/src/DataDrivenDiffEq.jl @@ -25,7 +25,8 @@ export reduce_basis, update! include("./dmdc.jl") export DMDc -export get_dynamics, get_input_map +export eigen, eigvals, eigvecs +export get_dynamics, get_input_map, dynamics include("./sindy.jl") export SInDy diff --git a/src/dmdc.jl b/src/dmdc.jl index 37667db36..9c6d1a2b0 100644 --- a/src/dmdc.jl +++ b/src/dmdc.jl @@ -1,9 +1,4 @@ -using OrdinaryDiffEq -using DataDrivenDiffEq -using LinearAlgebra -using Plots - -mutable struct DMDc{K, B, Q, P} +mutable struct DMDc{K, B, Q, P} <: abstractKoopmanOperator koopman::K forcing::B @@ -43,7 +38,6 @@ function DMDc(X::AbstractArray, Y::AbstractArray, Γ::AbstractArray; B::Abstract koopman = ExactDMD(Ã, Λ, ω, W, nothing, nothing) - return DMDc(koopman, B̃, nothing, nothing) end @@ -64,42 +58,14 @@ get_dynamics(m::DMDc) = m.koopman get_input_map(m::DMDc) = m.forcing # TODO this can be done better, maybe use macros -function DataDrivenDiffEq.dynamics(m::DMDc; discrete::Bool = true) - if discrete - dims = size(m.forcing) - nᵢ = length(dims)<=1 ? 1 : dims[2] - println(zeros(eltype(m.forcing), nᵢ)) - function zero_callback(u, p, t) - return zeros(eltype(m.forcing), nᵢ) - end - @inline function dudt_(u, p, t; y = zero_callback) - m.koopman.Ã * u + m.forcing .* y(u, p, t) - end - return dudt_ - else - throw(ErrorException("Continouos dynamics are not implemented right now.")) +function dynamics(m::DMDc; control = nothing) + control = isnothing(control) ? zero_callback(m) : control + @inline function dudt_(u, p, t; y = control) + m.koopman.Ã * u + m.forcing * y(u, p, t) end + return dudt_ end -X = [4 2 1 0.5 0.25; 7 0.7 0.07 0.007 0.0007] -U = [-4 -2 -1 -0.5] -B = Float32[1; 0] - -sys = DMDc(X, U) -sys = DMDc(X, U, B = B) - -isstable(sys) - -get_input_map(sys) - -get_dynamics(sys) - -eigen(sys) - -dudt_ = dynamics(sys) -dudt_(X[:, 1], [], 0.0) - -prob = DiscreteProblem(dudt_, X[:, 1], (0.0, 10.0)) -sol = solve(prob) - -plot(sol) +function zero_callback(m::DMDc) + return length(size(m.forcing)) <= 1 ? (u, p, t) -> zero(eltype(m.forcing)) : (u, p, t) -> zeros(eltype(m.forcing), size(m.forcing)[2]) +end From 37027d36aa3c1934a70ac0b9bd57286c3b056dd9 Mon Sep 17 00:00:00 2001 From: Julius Martensen Date: Thu, 5 Dec 2019 20:55:18 +0100 Subject: [PATCH 4/4] Added tests, ready for merge --- examples/DMDc_Examples.jl | 13 ++++++++----- src/dmdc.jl | 27 +++++++++++++++------------ test/runtests.jl | 30 +++++++++++++++++++++++++++++- 3 files changed, 52 insertions(+), 18 deletions(-) diff --git a/examples/DMDc_Examples.jl b/examples/DMDc_Examples.jl index efcc4c7fe..df666aa62 100644 --- a/examples/DMDc_Examples.jl +++ b/examples/DMDc_Examples.jl @@ -1,5 +1,5 @@ using DataDrivenDiffEq -using DifferentialEquations +using OrdinaryDiffEq using Plots gr() @@ -16,16 +16,19 @@ sys = DMDc(X, U, B = B) # Extract the DMD from inside DMDc get_dynamics(sys) # Acess all the other stuff -eigen(sys) .≈ eigen(get_dynamics(sys)) +eigen(sys) +eigvals(sys) +eigvecs(sys) +isstable(sys) # Get unforced dynamics dudt_ = dynamics(sys) prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.)) -sol_unforced = solve(prob) +sol_unforced = solve(prob, FunctionMap()) plot(sol_unforced) - +sol_unforced[:,:] # Create a system with cos control input to stabilize dudt_ = dynamics(sys, control = (u, p, t) -> -0.5u[1]) prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.)) -sol = solve(prob) +sol = solve(prob, FunctionMap()) plot!(sol) diff --git a/src/dmdc.jl b/src/dmdc.jl index 9c6d1a2b0..91304b00d 100644 --- a/src/dmdc.jl +++ b/src/dmdc.jl @@ -2,49 +2,51 @@ mutable struct DMDc{K, B, Q, P} <: abstractKoopmanOperator koopman::K forcing::B - Qₖ::Q + Qₖ::Q # TODO implement for update Pₖ::P end -function DMDc(X::AbstractArray, Y::AbstractArray, Γ::AbstractArray; B::AbstractArray = [], Δt::Float64 = 0.0) +function DMDc(X::AbstractArray, Y::AbstractArray, U::AbstractArray; B::AbstractArray = [], dt::Float64 = 0.0) @assert all(size(X) .== size(Y)) - @assert size(X)[2] .== size(Γ)[2] + @assert size(X)[2] .== size(U)[2] nₓ = size(X)[1] - nᵤ = size(Γ)[1] + nᵤ = size(U)[1] if isempty(B) - Ω = vcat(X, Γ) + Ω = vcat(X, U) G = Y * pinv(Ω) Ã = G[:, 1:nₓ] B̃ = G[:, nₓ+1:end] else - Ã = (Y - B*Γ)*pinv(X) + Ã = (Y - B*U)*pinv(X) B̃ = B end # Eigen Decomposition for solution Λ, W = eigen(Ã) - if Δt > 0.0 + # Right now now really useful, since we do not know + # how to deal with B + if dt > 0.0 # Casting Complex enforces results - ω = log.(Complex.(Λ)) / Δt + ω = log.(Complex.(Λ)) / dt else ω = [] end koopman = ExactDMD(Ã, Λ, ω, W, nothing, nothing) - + # TODO implement update method return DMDc(koopman, B̃, nothing, nothing) end -function DMDc(X::AbstractArray, Γ::AbstractArray; B::AbstractArray = [], Δt::Float64 = 0.0) - @assert size(X)[2]-1 == size(Γ)[2] "Provide consistent input data." - return DMDc(X[:, 1:end-1], X[:, 2:end], Γ, B = B, Δt = Δt) +function DMDc(X::AbstractArray, U::AbstractArray; B::AbstractArray = [], dt::Float64 = 0.0) + @assert size(X)[2]-1 == size(U)[2] "Provide consistent input data." + return DMDc(X[:, 1:end-1], X[:, 2:end], U, B = B, dt = dt) end # Some nice functions @@ -66,6 +68,7 @@ function dynamics(m::DMDc; control = nothing) return dudt_ end + function zero_callback(m::DMDc) return length(size(m.forcing)) <= 1 ? (u, p, t) -> zero(eltype(m.forcing)) : (u, p, t) -> zeros(eltype(m.forcing), size(m.forcing)[2]) end diff --git a/test/runtests.jl b/test/runtests.jl index 5cf7e119b..09c0caf97 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,7 +49,6 @@ end @testset "EDMD" begin - # Test for linear system function linear_sys(u, p, t) x = -0.9*u[1] @@ -97,6 +96,35 @@ end @test sol[:,:] ≈ s4[:,:] end + +@testset "DMDc" begin + # Define measurements from unstable system with known control input + X = [4 2 1 0.5 0.25; 7 0.7 0.07 0.007 0.0007] + U = [-4 -2 -1 -0.5] + B = Float32[1; 0] + + # But with a little more knowledge + sys = DMDc(X, U, B = B) + @test isa(get_dynamics(sys), ExactDMD) + @test sys.koopman.Ã ≈[1.5 0; 0 0.1] + @test get_input_map(sys) ≈ [1.0; 0.0] + @test !isstable(sys) + @test_nowarn eigen(sys) + + # Check the solution of an unforced and forced system against each other + dudt_ = dynamics(sys) + prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.)) + sol_unforced = solve(prob, FunctionMap()) + + dudt_ = dynamics(sys, control = (u, p, t) -> -0.5u[1]) + prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.)) + sol = solve(prob, FunctionMap()) + + @test all(abs.(diff(sol[1,:])) .< 1e-5) + @test sol[2,:] ≈ sol_unforced[2,:] +end + + @testset "SInDy" begin # Test the pendulum function pendulum(u, p, t)