From ce6008088eac9026e6e3ca5d5fd2df8b2aea26e8 Mon Sep 17 00:00:00 2001 From: mforets Date: Fri, 8 Mar 2019 12:39:55 -0300 Subject: [PATCH 1/3] add new option to discretize --- src/ReachSets/discretize.jl | 86 ++++++++++++++++++++++++++++++++++--- 1 file changed, 80 insertions(+), 6 deletions(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index dbdc0256..76bf874e 100644 --- a/src/ReachSets/discretize.jl +++ b/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,7 +645,7 @@ 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 @@ -631,6 +653,19 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous U0 = next_set(U, 1) Eψ0 = sih(Phi2Aabs * sih(A * U0)) + Ω0, Ud = _discretize_interpolation_inhomog(δ, U0, U, X0, ϕ, Einit, Eψ0, 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, Phi2Aabs, set_operations::Val{:lazy}) Ω0 = ConvexHull(X0, ϕ * X0 ⊕ δ*U0 ⊕ Eψ0 ⊕ Einit) if U isa ConstantInput @@ -647,6 +682,45 @@ 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(X0, 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, Phi2Aabs, set_operations::Val{:zonotope}) + Einit = convert(Zonotope, Einit) + Eψ0 = convert(Zonotope, Eψ0) + Z1 = X0 + ϕX0 = linear_map(ϕ, X0) + U0 = convert(Zonotope, U0) + δI = Matrix(δ*I, size(U0)) + δ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 From 268bef00a3fde7eed84362331e127c59fa886618 Mon Sep 17 00:00:00 2001 From: mforets Date: Fri, 8 Mar 2019 12:57:32 -0300 Subject: [PATCH 2/3] fix arguments --- src/ReachSets/discretize.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 76bf874e..229f6a0d 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -653,7 +653,7 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous U0 = next_set(U, 1) Eψ0 = sih(Phi2Aabs * sih(A * U0)) - Ω0, Ud = _discretize_interpolation_inhomog(δ, U0, U, X0, ϕ, Einit, Eψ0, Phi2Aabs, Val(set_operations)) + Ω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 @@ -665,7 +665,7 @@ function _discretize_interpolation_homog(X0, ϕ, Einit, set_operations::Val{:laz end # version using lazy sets and operations -function _discretize_interpolation_inhomog(δ, U0, U, X0, ϕ, Einit, Eψ0, Phi2Aabs, set_operations::Val{:lazy}) +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 @@ -695,7 +695,7 @@ function _discretize_interpolation_homog(X0, ϕ, Einit, set_operations::Val{:zon end # version using concrete operations with zonotopes -function _discretize_interpolation_inhomog(δ, U0, U, X0, ϕ, Einit, Eψ0, Phi2Aabs, set_operations::Val{:zonotope}) +function _discretize_interpolation_inhomog(δ, U0, U, X0, ϕ, Einit, Eψ0, A, sih, Phi2Aabs, set_operations::Val{:zonotope}) Einit = convert(Zonotope, Einit) Eψ0 = convert(Zonotope, Eψ0) Z1 = X0 From 1170f92abe0b3af4e7417ddd1fef354d3247bae5 Mon Sep 17 00:00:00 2001 From: mforets Date: Fri, 8 Mar 2019 16:14:26 -0300 Subject: [PATCH 3/3] fix usage of msum with zonotopes --- src/ReachSets/discretize.jl | 9 +++++---- test/ReachSets/unit_discretization.jl | 18 ++++++++++++++++++ 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 229f6a0d..4b83fc54 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -689,19 +689,20 @@ end function _discretize_interpolation_homog(X0, ϕ, Einit, set_operations::Val{:zonotope}) Einit = convert(Zonotope, Einit) Z1 = convert(Zonotope, X0) - Z2 = linear_map(ϕ, minkowski_sum(X0, Einit)) + 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 = X0 - ϕX0 = linear_map(ϕ, X0) + Z1 = convert(Zonotope, X0) + ϕX0 = linear_map(ϕ, Z1) U0 = convert(Zonotope, U0) - δI = Matrix(δ*I, size(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) diff --git a/test/ReachSets/unit_discretization.jl b/test/ReachSets/unit_discretization.jl index 560bbae1..7d33e356 100644 --- a/test/ReachSets/unit_discretization.jl +++ b/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 # ===================================================================