Merge pull request #537 from JuliaReach/mforets/TMJets

Add TMJets algorithm

.travis.yml

@@ -21,6 +21,8 @@ script:
              Pkg.develop("LazySets");
              Pkg.develop("MathematicalSystems");
              Pkg.develop("HybridSystems");
+             Pkg.add(PackageSpec(name="TaylorModels", rev="validatedODEs"));
+             Pkg.develop("TaylorSeries");
              Pkg.clone(pwd());
              Pkg.build("Reachability");
              Pkg.test("Reachability"; coverage=true)'

REQUIRE

@@ -3,11 +3,14 @@ Compat
 Expokit
 HybridSystems
 LazySets
-MathematicalSystems 0.6.3
+MathematicalSystems 0.6.4
 Memento
 Optim
 Plots
 ProgressMeter
 RecipesBase
 Reexport
 Suppressor
+IntervalArithmetic
+TaylorModels
+TaylorSeries

src/Options/default_options.jl

@@ -92,7 +92,7 @@ function validate_solver_options_and_add_default_values!(options::Options)::Opti
                 error("No property has been defined.")
             end
         elseif key == :property
-            expected_type = Union{Nothing, Property, Dict{Int, <:Property}}
+            expected_type = Union{Nothing, Property, Dict{Int, <:Property}, Function}
         elseif key == :T
             expected_type = Float64
             domain_constraints = (v::Float64  ->  v > 0.)

src/ReachSets/ContinuousPost/TMJets/TMJets.jl

@@ -0,0 +1,21 @@
+export TMJets
+
+using IntervalArithmetic: IntervalBox
+
+struct TMJets <: ContinuousPost
+    options::TwoLayerOptions
+
+    function TMJets(𝑂::Options)
+        𝑂new = validate_and_wrap_options(𝑂, options_TMJets())
+        return new(𝑂new)
+    end
+end
+
+# convenience constructor from pairs of symbols
+TMJets(𝑂::Pair{Symbol, <:Any}...) = TMJets(Options(Dict{Symbol, Any}(𝑂)))
+
+# default options (they are added in the function validate_and_wrap_options)
+TMJets() = TMJets(Options())
+
+include("init.jl")
+include("post.jl")

src/ReachSets/ContinuousPost/TMJets/init.jl

@@ -0,0 +1,29 @@
+# out-of-place initialization
+init(𝒫::TMJets, 𝑆::AbstractSystem, 𝑂::Options) = init!(𝒫, 𝑆, copy(𝑂))
+
+function options_TMJets()
+
+    𝑂spec = Vector{OptionSpec}()
+
+    # step size and termination criteria
+    push!(𝑂spec, OptionSpec(:abs_tol, 1e-10, domain=Float64, info="absolute tolerance"))
+    push!(𝑂spec, OptionSpec(:max_steps, 500, domain=Int, info="use this maximum number of steps in the validated integration"))
+
+    # approximation options
+    push!(𝑂spec, OptionSpec(:orderT, 10, domain=Int, info="order of the Taylor model in t"))
+    push!(𝑂spec, OptionSpec(:orderQ, 2, domain=Int, info="order of the Taylor model for Jet transport variables"))
+
+    return 𝑂spec
+end
+
+# in-place initialization
+function init!(𝒫::TMJets, 𝑆::AbstractSystem, 𝑂::Options)
+
+    # state dimension
+    𝑂[:n] = statedim(𝑆)
+
+    # adds default values for unspecified options
+    𝑂init = validate_solver_options_and_add_default_values!(𝑂)
+
+    return 𝑂init
+end

src/ReachSets/ContinuousPost/TMJets/post.jl

@@ -0,0 +1,76 @@
+using TaylorModels: validated_integ
+using TaylorSeries: set_variables
+using LazySets.Approximations: box_approximation
+
+function post(𝒜::TMJets,
+              𝑃::InitialValueProblem{<:Union{BBCS, CBBCS, CBBCCS}, <:LazySet},
+              𝑂_global::Options)
+
+    # ==================================
+    # Initialization
+    # ==================================
+
+    𝑂 = merge(𝒜.options.defaults, 𝑂_global, 𝒜.options.specified)
+
+    # system of ODEs
+    f! = 𝑃.s.f
+    n = statedim(𝑃)
+
+    # initial time and time horizon
+    t0 = 0.0
+    T = 𝑂[:T]
+
+    # maximum allowed number of steps
+    max_steps = 𝑂[:max_steps]
+
+    # unrap algorithm-specific options
+    abs_tol, orderQ, orderT = 𝑂[:abs_tol], 𝑂[:orderQ], 𝑂[:orderT]
+
+    # initial sets
+    box_x0 = box_approximation(𝑃.x0)
+    q0 = center(box_x0)
+    δq0 = IntervalBox(low(box_x0)-q0, high(box_x0)-q0)
+
+    # fix the working variables and maximum order in the global
+    # parameters struct (_params_TaylorN_)
+    set_variables("x", numvars=length(q0), order=2*orderQ)
+
+    # define the property
+    if 𝑂[:mode] == "check"
+        property = 𝑂[:property]
+    elseif 𝑂[:mode] == "reach"
+        property = (t, x) -> true
+    end
+
+    # =====================
+    # Flowpipe computation
+    # =====================
+
+    info("Reachable States Computation...")
+    @timing begin
+        tTM, xTM = validated_integ(f!, q0, δq0, t0, T, orderQ, orderT, abs_tol,
+                                   maxsteps=max_steps, check_property=property)
+    end
+
+    # convert to hyperrectangle and wrap around the reach solution
+    N = length(xTM)
+    Rsets = Vector{ReachSet{Hyperrectangle, Float64}}(undef, N-1)
+    @inbounds for i in 1:N-1
+        Hi = convert(Hyperrectangle, xTM[i])
+        t0 = tTM[i]; t1 = tTM[i+1]
+        Rsets[i] = ReachSet{Hyperrectangle, Float64}(Hi, t0, t1)
+    end
+
+    Rsol = ReachSolution(Rsets, 𝑂)
+
+    # ===========
+    # Projection
+    # ===========
+
+    if 𝑂[:project_reachset]
+        info("Projection...")
+        Rsol = @timing project(Rsol)
+    end
+
+    return Rsol
+end

src/ReachSets/ReachSets.jl

@@ -100,6 +100,8 @@ include("ContinuousPost/BFFPSV18/reach_blocks_wrapping_effect.jl")
 
 include("ContinuousPost/GLGM06/GLGM06.jl")
 
+include("ContinuousPost/TMJets/TMJets.jl")
+
 # ========================
 # Reachability Algorithms
 # ========================

src/Utils/normalization.jl

@@ -61,6 +61,9 @@ function normalize(system::AbstractSystem)
     throw(ArgumentError("the system type $(typeof(system)) is currently not supported"))
 end
 
+# "black-box" like systems are not normalized; algorithms should handle this
+normalize(S::BBCS) = S
+
 # x' = Ax, in the continuous case
 # x+ = Ax, in the discrete case
 for (L_S, CL_S) in ((:LCS, :CLCS), (:LDS, :CLDS))

src/Utils/systems.jl

@@ -14,8 +14,11 @@ const CACCS = ConstrainedAffineControlContinuousSystem
 const CACDS = ConstrainedAffineControlDiscreteSystem
 const CACS = ConstrainedAffineContinuousSystem
 const CADS = ConstrainedAffineDiscreteSystem
+const BBCS = BlackBoxContinuousSystem
+const CBBCS = ConstrainedBlackBoxContinuousSystem
+const CBBCCS = ConstrainedBlackBoxContinuousSystem
 
-export LCS, LDS, CLCS, CLDS, CLCCS, CLCDS, CACCS, CACDS, IVP
+export LCS, LDS, CLCS, CLDS, CLCCS, CLCDS, CACCS, CACDS, IVP, BBCS, CBBCS, CBBCCS
 
 import Base: *
 import LazySets.constrained_dimensions

test/Reachability/solve_continuous.jl

@@ -192,3 +192,30 @@ X = HalfSpace([-1.0, 0.0, 0.0, 0.0], 0.0)
 X0 = BallInf(zeros(4), 0.5)
 p = IVP(ConstrainedAffineContinuousSystem(A, b, X), X0)
 sol = solve(p, :T=>0.1)
+
+# ==================================
+# Test TMJets with Van der Pol model
+# ==================================
+using TaylorIntegration
+using TaylorModels: @taylorize
+using Reachability: solve
+
+@taylorize function vanderPol!(t, x, dx)
+    local μ = 1.0
+    dx[1] = x[2]
+    dx[2] = (μ * x[2]) * (1 - x[1]^2) - x[1]
+    return dx
+end
+
+𝑆 = BlackBoxContinuousSystem(vanderPol!, 2)
+X0 = Hyperrectangle(low=[1.25, 2.35], high=[1.55, 2.45])
+𝑃 = InitialValueProblem(𝑆, X0)
+
+# reach mode
+𝑂 = Options(:T=>7.0, :mode=>"reach")
+solve(𝑃, 𝑂, op=TMJets(:abs_tol=>1e-10, :orderT=>10, :orderQ=>2));
+
+# check mode
+property=(t,x)->x[2] <= 2.75
+𝑂 = Options(:T=>7.0, :mode=>"check", :property=>property)
+solve(𝑃, 𝑂, op=TMJets(:abs_tol=>1e-10, :orderT=>10, :orderQ=>2))