Merge 37027d3 into ac9b0e2

examples/DMDc_Examples.jl

@@ -0,0 +1,34 @@
+using DataDrivenDiffEq
+using OrdinaryDiffEq
+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)
+eigvals(sys)
+eigvecs(sys)
+isstable(sys)
+# Get unforced dynamics
+dudt_ = dynamics(sys)
+prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.))
+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, FunctionMap())
+
+plot!(sol)

src/DataDrivenDiffEq.jl

@@ -24,6 +24,11 @@ export ExtendedDMD
 export dynamics, linear_dynamics
 export reduce_basis, update!
 
+include("./dmdc.jl")
+export DMDc
+export eigen, eigvals, eigvecs
+export get_dynamics, get_input_map, dynamics
+
 include("./sindy.jl")
 export SInDy
 

src/dmdc.jl

@@ -0,0 +1,74 @@
+mutable struct DMDc{K, B, Q, P} <: abstractKoopmanOperator
+    koopman::K
+    forcing::B
+
+    Qₖ::Q # TODO implement for update
+    Pₖ::P
+end
+
+
+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(U)[2]
+
+
+    nₓ = size(X)[1]
+    nᵤ = size(U)[1]
+
+    if isempty(B)
+        Ω = vcat(X, U)
+        G = Y * pinv(Ω)
+
+        Ã = G[:, 1:nₓ]
+        B̃ = G[:, nₓ+1:end]
+    else
+        Ã = (Y - B*U)*pinv(X)
+        B̃ = B
+    end
+
+    # Eigen Decomposition for solution
+    Λ, W = eigen(Ã)
+
+    # 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.(Λ)) / dt
+    else
+        ω = []
+    end
+
+    koopman = ExactDMD(Ã, Λ, ω, W, nothing, nothing)
+    # TODO implement update method
+    return DMDc(koopman, B̃, nothing, nothing)
+end
+
+
+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
+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 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
+
+
+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

test/runtests.jl

@@ -50,7 +50,6 @@ end
 
 
 @testset "EDMD" begin
-
     # Test for linear system
     function linear_sys(u, p, t)
         x = -0.9*u[1]
@@ -98,6 +97,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)