Merge pull request #75 from JuliaReach/mforets/various_fixes

Various fixes

docs/src/man/nonlinear.md

@@ -1,3 +1,75 @@
 # Nonlinear ordinary differential equations
 
-- lotka-volterra
+In this section we illustrate the flowpipe computation for a nonlinear system.
+
+## Model description
+
+Our running example is the [Lotka-Volterra model](https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations). The 2-dimensional Lotka-Volterra system depicts the populations change of a class of predators and a class of
+preys. The growth rate of preys’ population $x$ over time is given by
+
+```math
+\dot{x} = x\cdot (\alpha - \beta \cdot y)
+```
+wherein  $\alpha, \beta$ are constant parameters and $y$ is the population of predators.
+
+It gives that the number of preys grows exponentially without predation.
+
+The population growth of predators is governed by the differential equation
+
+```math
+\dot{y} = -y\cdot (\gamma - \delta\cdot x)
+```
+wherein  $\gamma, \delta$ are constant parameters.
+
+We set those parameters as  $\alpha = 1.5 ,  \beta = 1 ,  \gamma = 3$  and  $\delta = 1$.
+
+```@example lotka_volterra
+using ReachabilityAnalysis
+
+@taylorize function lotka_volterra!(du, u, p, t)
+    local α, β, γ, δ = 1.5, 1.0, 3.0, 1.0
+    du[1] = u[1] * (α - β*u[2])
+    du[2] = -u[2] * (γ - δ*u[1])
+    return du
+end
+```
+
+## Reachability settings
+
+The reachability settings are taken from [this resource](https://ths.rwth-aachen.de/research/projects/hypro/lotka-volterra/).
+
+We consider the initial set  $x\in [4.8,5.2], y \in [1.8,2.2]$.
+
+```@example lotka_volterra
+X₀ = Hyperrectangle(low=[4.8, 1.8], high=[5.2, 2.2])
+
+prob = @ivp(x' = lotka_volterra!(x), dim: 2, x(0) ∈ X₀)
+```
+
+## Results
+
+We compute the flowpipe using the TMJets algorithm for the time horizon $[0,5]$:
+
+```@example lotka_volterra
+sol = solve(prob, T=5.0)
+
+setrep(sol)
+```
+
+We can change to the zonotopic overapproximation of the flowpipe using
+the `overapproximate` function:
+
+```@example lotka_volterra
+sol = overapproximate(sol, Zonotope)
+
+setrep(sol)
+```
+
+Finally we plot the solution in phase-space:
+
+```@example lotka_volterra
+using Plots
+
+plot(sol, vars=(1, 2), xlab="x", ylab="y", lw=0.2, color=:lightblue, lab="Flowpipe")
+plot!(X₀, color=:orange, lab="Xo")
+```

src/Algorithms/ASB07/post.jl

@@ -27,7 +27,7 @@ function post(alg::ASB07,
 
     # this algorithm requires Ω0 to be a zonotope
     Ω0 = _convert_or_overapproximate(Zonotope, Ω0)
-    Ω0 = reduce_order(Ω0, max_order)
+    Ω0 = _reduce_order(Ω0, max_order)
     c0 = center(Ω0)
     G0 = genmat(Ω0)
     N = eltype(Ω0)

src/Algorithms/GLGM06/post.jl

@@ -26,7 +26,7 @@ function post(alg::GLGM06, ivp::IVP{<:AbstractContinuousSystem}, tspan; kwargs..
 
     # this algorithm requires Ω0 to be a zonotope
     Ω0 = _convert_or_overapproximate(Zonotope, Ω0)
-    Ω0 = reduce_order(Ω0, max_order)
+    Ω0 = _reduce_order(Ω0, max_order)
 
     # reconvert the set of initial states, if needed
     force_static = haskey(kwargs, :force_static) ? kwargs[:force_static] : false

src/Algorithms/GLGM06/reach.jl

@@ -232,12 +232,12 @@ function reach_inhomog_GLGM06!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
     k = 2
     while k <= NSTEPS
         Rₖ = minkowski_sum(linear_map(Φ_power_k, Ω0), Wk₊, remove_zero_generators=false)
-        Rₖ = reduce_order(Rₖ, max_order)
+        Rₖ = _reduce_order(Rₖ, max_order)
         Δt += δ
         F[k] = ReachSet(Rₖ, Δt)
 
         Wk₊ = minkowski_sum(Wk₊, linear_map(Φ_power_k, U), remove_zero_generators=false)
-        Wk₊ = reduce_order(Wk₊, max_order)
+        Wk₊ = _reduce_order(Wk₊, max_order)
 
         mul!(Φ_power_k_cache, Φ_power_k, Φ)
         copyto!(Φ_power_k, Φ_power_k_cache)

src/Algorithms/TMJets/TMJets.jl

@@ -8,14 +8,14 @@ by Luis Benet and David Sanders in `TalorModels.jl`.
 
 - `max_steps` -- (optional, default: `1000`) maximum number of steps in the
                  validated integration ``x' = f(x)``
-- `abs_tol`   -- (optional, default: `1e-10`) absolute tolerance
+- `abs_tol`   -- (optional, default: `1e-15`) absolute tolerance
 - `orderT`    -- (optional, default: `8`) order of the Taylor model in time
 - `orderQ`    -- (optional, default: `2`) order of the Taylor models for jet
                  transport variales
 """
 @with_kw struct TMJets{N} <: AbstractContinuousPost
-    max_steps::Int=1000
-    abs_tol::N=1e-10
+    max_steps::Int=2000
+    abs_tol::N=1e-15
     orderT::Int=8
     orderQ::Int=2
     # setrep::ST=Zonotope{Float64, Vector{Float64}, Matrix{Float64}}

src/Flowpipes/flowpipes.jl

@@ -181,6 +181,11 @@ function project(fp::Flowpipe, vars::NTuple{D, T}) where {D, T<:Integer}
     end
 end
 
+# convenience alias to match the usage in the plot recipe
+function project(fp::Flowpipe; vars::NTuple{D, T}) where {D, T<:Integer}
+    return projct(fp, vars)
+end
+
 function Base.similar(fp::Flowpipe{N, RT}) where {N, RT<:AbstractReachSet{N}}
    return Flowpipe(Vector{RT}())
 end
@@ -351,8 +356,8 @@ function Base.similar(fp::HybridFlowpipe{N, RT, FT}) where {N, RT, FT}
     return HybridFlowpipe(Vector{FT}())
 end
 
-# first searches the flowpipe that contains `t` in its time-span, then the
-# corresponding reach-set
+# first searches the flowpipe that contains `t` in its time-span, then searches
+# inside that flowpipe for the corresponding reach-set
 function (fp::HybridFlowpipe)(t::Number)
     Fk = array(fp)
     i = 1

src/Flowpipes/reachsets.jl

@@ -335,7 +335,8 @@ end
 
 # interface functions
 set(R::SparseReachSet) = R.X
-setrep(R::SparseReachSet{N, ST}) where {N, ST<:LazySet{N}} = ST
+setrep(::SparseReachSet{N, ST}) where {N, ST<:LazySet{N}} = ST
+setrep(::Type{SparseReachSet{N, ST}}) where {N, ST<:LazySet{N}} = ST
 tstart(R::SparseReachSet) = inf(R.Δt)
 tend(R::SparseReachSet) = sup(R.Δt)
 tspan(R::SparseReachSet) = R.Δt
@@ -488,7 +489,8 @@ end
 
 # interface functions
 @inline set(R::TaylorModelReachSet) = R.X
-@inline setrep(R::TaylorModelReachSet{N}) where {N} = Vector{TaylorModel1{TaylorN{N}, N}}
+@inline setrep(::TaylorModelReachSet{N}) where {N} = Vector{TaylorModel1{TaylorN{N}, N}}
+@inline setrep(::Type{TaylorModelReachSet{N}}) where {N} = Vector{TaylorModel1{TaylorN{N}, N}}
 @inline tstart(R::TaylorModelReachSet) = inf(R.Δt)
 @inline tend(R::TaylorModelReachSet) = sup(R.Δt)
 @inline tspan(R::TaylorModelReachSet) = R.Δt
@@ -618,7 +620,8 @@ end
 # implement abstract reachset interface
 # TODO: use HPolyhedron or HPolytope if the set is bounded or not
 set(R::TemplateReachSet) = HPolytope([HalfSpace(R.dirs[i], R.sf[i]) for i in eachindex(R.sf)])
-setrep(R::TemplateReachSet{N, VN}) where {N, VN} = HPolyhedron{N, VN}
+setrep(::Type{TemplateReachSet{N, VN}}) where {N, VN} = HPolyhedron{N, VN}
+setrep(::TemplateReachSet{N, VN}) where {N, VN} = HPolyhedron{N, VN}
 tspan(R::TemplateReachSet) = R.Δt
 tstart(R::TemplateReachSet) = inf(R.Δt)
 tend(R::TemplateReachSet) = sup(R.Δt)

src/Flowpipes/solutions.jl

@@ -104,3 +104,8 @@ end
 function project(sol::ReachSolution{<:Flowpipe}, args...)
     return project(sol.F, args...)
 end
+
+# convenience alias to match the usage in the plot recipe
+function project(sol::ReachSolution{<:Flowpipe}; vars::NTuple{D, T}) where {D, T}
+    return project(sol.F, vars)
+end