Merge pull request #528 from JuliaReach/mforets/512

#512 - Add new option to discretize

src/ReachSets/discretize.jl

@@ -4,7 +4,7 @@ const CLCDS = ConstrainedLinearControlDiscreteSystem
 @inline Id(n) = Matrix(1.0I, n, n)
 
 """
-    discretize(𝑆, δ; [algorithm], [exp_method], [sih_method])
+    discretize(𝑆, δ; [algorithm], [exp_method], [sih_method], [set_operations])
 
 Apply an approximation model to `S` obtaining a discrete initial value problem.
 
@@ -42,6 +42,15 @@ Apply an approximation model to `S` obtaining a discrete initial value problem.
     - `"lazy"`     -- compute a wrapper set type around symmetric interval hull
                       in a lazy way using `SymmetricIntervalHull`
 
+- `set_operations`  -- (optional, default: `"lazy"`) set operations used for the
+                       discretized initial states and transformed inputs, choose among:
+
+    - `"lazy"`     -- use lazy convex hull for the initial states and lazy linear
+                      map for the inputs
+    - `"zonotope"` -- use concrete zonotope operations (linear map and Minkowski sum)
+                      and return zonotopes for both the initial states and the
+                      inputs of the discretized system
+
 ### Output
 
 The initial value problem of a discrete system.
@@ -84,14 +93,23 @@ function discretize(𝑆::InitialValueProblem{<:AbstractContinuousSystem},
                     δ::Float64;
                     algorithm::String="forward",
                     exp_method::String="base",
-                    sih_method::String="concrete")
+                    sih_method::String="concrete",
+                    set_operations::String="lazy")
 
     if algorithm in ["forward", "backward"]
         return discretize_interpolation(𝑆, δ, algorithm=algorithm,
-                    exp_method=exp_method, sih_method=sih_method)
+                    exp_method=exp_method, sih_method=sih_method, set_operations=set_operations)
+
     elseif algorithm == "firstorder"
+        if set_operations != "lazy"
+            throw(ArgumentError("the algorithm $algorithm with set operations=$set_operations is not implemented"))
+        end
         return discretize_firstorder(𝑆, δ, exp_method=exp_method)
+
     elseif algorithm == "nobloating"
+        if set_operations != "lazy"
+            throw(ArgumentError("the algorithm $algorithm with set operations=$set_operations is not implemented"))
+        end
         return discretize_nobloating(𝑆, δ, exp_method=exp_method)
     else
         throw(ArgumentError("the algorithm $algorithm is unknown"))
@@ -594,7 +612,11 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous
                                    δ::Float64;
                                    algorithm::String="forward",
                                    exp_method::String="base",
-                                   sih_method::String="concrete")
+                                   sih_method::String="concrete",
+                                   set_operations::String="lazy")
+
+    # used to dispatch on the value of the set operation
+    set_operations = Symbol(set_operations)
 
     if sih_method == "concrete"
         sih = symmetric_interval_hull
@@ -623,14 +645,27 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous
 
     # early return for homogeneous systems
     if inputdim(𝑆) == 0
-        Ω0 = ConvexHull(X0, ϕ * X0 ⊕ Einit)
+        Ω0 = _discretize_interpolation_homog(X0, ϕ, Einit, Val(set_operations))
         return IVP(LDS(ϕ), Ω0)
     end
 
     U = inputset(𝑆)
     U0 = next_set(U, 1)
 
     Eψ0 = sih(Phi2Aabs * sih(A * U0))
+    Ω0, Ud = _discretize_interpolation_inhomog(δ, U0, U, X0, ϕ, Einit, Eψ0, A, sih, Phi2Aabs, Val(set_operations))
+
+    return IVP(CLCDS(ϕ, Id(size(A, 1)), nothing, Ud), Ω0)
+end
+
+# version using lazy sets and operations
+function _discretize_interpolation_homog(X0, ϕ, Einit, set_operations::Val{:lazy})
+    Ω0 = ConvexHull(X0, ϕ * X0 ⊕ Einit)
+    return Ω0
+end
+
+# version using lazy sets and operations
+function _discretize_interpolation_inhomog(δ, U0, U, X0, ϕ, Einit, Eψ0, A, sih, Phi2Aabs, set_operations::Val{:lazy})
     Ω0 = ConvexHull(X0, ϕ * X0 ⊕ δ*U0 ⊕ Eψ0 ⊕ Einit)
 
     if U isa ConstantInput
@@ -647,6 +682,46 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous
     else
         throw(ArgumentError("input of type $(typeof(U)) is not allowed"))
     end
+    return Ω0, Ud
+end
 
-    return IVP(CLCDS(ϕ, Id(size(A, 1)), nothing, Ud), Ω0)
+# version using concrete operations with zonotopes
+function _discretize_interpolation_homog(X0, ϕ, Einit, set_operations::Val{:zonotope})
+    Einit = convert(Zonotope, Einit)
+    Z1 = convert(Zonotope, X0)
+    Z2 = linear_map(ϕ, minkowski_sum(Z1, Einit))
+    Ω0 = overapproximate(ConvexHull(Z1, Z2), Zonotope)
+    return Ω0
+end
+
+# version using concrete operations with zonotopes
+function _discretize_interpolation_inhomog(δ, U0, U, X0, ϕ, Einit, Eψ0, A, sih, Phi2Aabs, set_operations::Val{:zonotope})
+    n = size(A, 1)
+    Einit = convert(Zonotope, Einit)
+    Eψ0 = convert(Zonotope, Eψ0)
+    Z1 = convert(Zonotope, X0)
+    ϕX0 = linear_map(ϕ, Z1)
+    U0 = convert(Zonotope, U0)
+    δI = Matrix(δ*I, n, n)
+    δU0 = linear_map(δI, U0)
+    Z2 = reduce(minkowski_sum, [ϕX0, δU0, Eψ0, Einit])
+    Ω0 = overapproximate(ConvexHull(Z1, Z2), Zonotope)
+
+    if U isa ConstantInput
+        Ud = ConstantInput(minkowski_sum(δU0, Eψ0))
+
+    elseif U isa VaryingInput
+        Ud = Vector{LazySet}(undef, length(U))
+        for (k, Uk) in enumerate(U)
+            Eψk = convert(Zonotope, sih(Phi2Aabs * sih(A * Uk)))
+            Uk = convert(Zonotope, Uk)
+            δUk = linear_map(δI, Uk)
+            Ud[k] = minkowski_sum(δUk, Eψk)
+        end
+        Ud = VaryingInput(Ud)
+
+    else
+        throw(ArgumentError("input of type $(typeof(U)) is not allwed"))
+    end
+    return Ω0, Ud
 end

test/ReachSets/unit_discretization.jl

@@ -31,6 +31,16 @@ let # preventing global scope
     discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="backward")
     @test inputdim(discr_sys_homog) == 0
 
+    # bloating, forward interpolation using lazy operations (default)
+    discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="forward", set_operations="lazy")
+    @test discr_sys_homog.x0 isa ConvexHull
+    @test inputdim(discr_sys_homog) == 0
+
+    # using concrete zonotopes
+    discr_sys_homog = discretize(cont_sys_homog, δ, set_operations="zonotope")
+    @test discr_sys_homog.x0 isa Zonotope
+    @test inputdim(discr_sys_homog) == 0
+
     # ===============================================================
     # Discretization of a continuous-time system with constant input
     # ===============================================================
@@ -66,6 +76,14 @@ let # preventing global scope
     discr_sys = discretize(cont_sys, δ, algorithm="firstorder")
     @test inputdim(discr_sys) == 4
 
+    # using concrete zonotopes
+    # since a Ball2 cannot be converted to a Zonotope, we use a box instead
+    U = ConstantInput(BallInf(ones(4), 0.5))
+    cont_sys = IVP(CLCCS(A, B, nothing, U), X0)
+    discr_sys = discretize(cont_sys, δ, set_operations="zonotope")
+    @test discr_sys.x0 isa Zonotope
+    @test inputdim(discr_sys) == 4
+
     # ===================================================================
     # Discretization of a continuous-time system with time-varying input
     # ===================================================================