#1015 - Add split method for zonotopes

docs/src/lib/representations.md

@@ -592,7 +592,9 @@ linear_map(::AbstractMatrix{N}, ::Zonotope{N}) where {N<:Real}
 scale(::Real, ::Zonotope)
 ngens(::Zonotope)
 reduce_order(::Zonotope, r)
+split(::Zonotope, ::Int)
 ```
+
 Inherited from [`LazySet`](@ref):
 * [`norm`](@ref norm(::LazySet, ::Real))
 * [`radius`](@ref radius(::LazySet, ::Real))

src/Zonotope.jl

@@ -1,5 +1,6 @@
 import Base: rand,
-             ∈
+             ∈,
+             split
 
 export Zonotope,
        order,
@@ -438,6 +439,58 @@ function reduce_order(Z::Zonotope{N}, r)::Zonotope{N} where {N<:Real}
     return Zonotope(c, Gred)
 end
 
+"""
+    split(Z::Zonotope, j::Int)
+
+Return two zonotopes obtained by splitting the given zonotope.
+
+### Input
+
+- `Z` -- zonotope
+- `j` -- index of the generator to be split
+
+### Output
+
+The zonotope obtained by splitting `Z` into two zonotopes such that
+their union is `Z` and their intersection is possibly non-empty.
+
+### Algorithm
+
+This function implements [Prop. 3, 1], that we state next. The zonotope
+``Z = ⟨c, g^{(1, …, p)}⟩`` is split into:
+
+```math
+Z₁ = ⟨c - \\frac{1}{2}g^{(j)}, (g^{(1, …,j-1)}, \\frac{1}{2}g^{(j)}, g^{(j+1, …, p)})⟩ \\\\
+Z₂ = ⟨c + \\frac{1}{2}g^{(j)}, (g^{(1, …,j-1)}, \\frac{1}{2}g^{(j)}, g^{(j+1, …, p)})⟩,
+```
+such that ``Z₁ ∪ Z₂ = Z`` and ``Z₁ ∩ Z₂ = Z^*``, where
+
+```math
+Z^* = ⟨c, (g^{(1,…,j-1)}, g^{(j+1,…, p)})⟩.
+```
+
+[1] *Althoff, M., Stursberg, O., & Buss, M. (2008). Reachability analysis of
+nonlinear systems with uncertain parameters using conservative linearization.
+In Proc. of the 47th IEEE Conference on Decision and Control.*
+"""
+function split(Z::Zonotope, j::Int)
+    @assert 1 <= j <= ngens(Z) "cannot split a zonotope with $(ngens(Z)) generators along index $j"
+    c, G = Z.center, Z.generators
+    Gj = G[:, j]
+    Gj_half = Gj / 2
+
+    c₁ = c - Gj_half
+    c₂ = c + Gj_half
+
+    G₁ = copy(G)
+    G₁[:, j] = Gj_half
+    G₂ = copy(G₁)
+
+    Z₁ = Zonotope(c₁, G₁)
+    Z₂ = Zonotope(c₂, G₂)
+    return Z₁, Z₂
+end
+
 """
     constraints_list(P::Zonotope{N}
                     )::Vector{LinearConstraint{N}} where {N<:Real}

test/unit_Zonotope.jl

@@ -96,6 +96,11 @@ for N in [Float64, Rational{Int}, Float32]
     Z = convert(Zonotope, Hyperrectangle(N[2, 3], N[4, 5]))
     @test Z.center == N[2, 3] && diag(Z.generators) == N[4, 5]
     convert(Zonotope, BallInf(N[5, 3], N(2)))
+
+    # split a zonotope
+    Z = Zonotope(N[0, 0], N[1 1; -1 1])
+    Z1, Z2 = split(Z, 1) # in this case the splitting is exact
+    @test Z1 ⊆ Z && Z2 ⊆ Z
 end
 
 for N in [Float64, Rational{Int}]