SPZ overapproximation of CH(::SPZ,::SPZ)

src/Approximations/overapproximate.jl

@@ -661,3 +661,42 @@ function load_paving_overapproximation()
         end
     end
 end  # quote / load_paving_overapproximation
+
+"""
+# Extended help
+
+    overapproximate(CH::ConvexHull{N,<:SparsePolynomialZonotope,<:SparsePolynomialZonotope},
+                    ::Type{<:SparsePolynomialZonotope}) where {N}
+
+### Algorithm
+
+This method implements Proposition 3.1.28 in [1].
+
+[1] Kochdumper, Niklas. *Extensions of polynomial zonotopes and their application
+    to verification of cyber-physical systems.* PhD diss., Technische Universität
+    München, 2022.
+"""
+function overapproximate(CH::ConvexHull{N,<:AbstractSparsePolynomialZonotope,<:AbstractSparsePolynomialZonotope},
+                         ::Type{<:SparsePolynomialZonotope}) where {N}
+    PZ₁ = first(CH)
+    c₁ = center(PZ₁)
+    G₁ = genmat_dep(PZ₁)
+    GI₁ = genmat_indep(PZ₁)
+    E₁ = expmat(PZ₁)
+    PZ₂ = second(CH)
+    c₂ = center(PZ₂)
+    G₂ = genmat_dep(PZ₂)
+    GI₂ = genmat_indep(PZ₂)
+    E₂ = expmat(PZ₂)
+
+    # zonotope overapproximation of convex hull of zonotopes
+    cZ = zeros(N, dim(CH))
+    Z = overapproximate(ConvexHull(Zonotope(cZ, GI₁), Zonotope(cZ, GI₂)))
+
+    # exact convex hull of simple polynomial zonotopes
+    PZ₁_bar = SimpleSparsePolynomialZonotope(c₁, G₁, E₁)
+    PZ₂_bar = SimpleSparsePolynomialZonotope(c₂, G₂, E₂)
+    PZ_bar = convex_hull(PZ₁_bar, PZ₂_bar)
+
+    return SparsePolynomialZonotope(center(PZ_bar), genmat_dep(PZ_bar), genmat(Z), expmat(PZ_bar))
+end

test/Approximations/overapproximate.jl

@@ -312,6 +312,12 @@ for N in [Float64, Float32]
     # overapproximation with HPolyhedron
     H = HalfSpace(N[1, 0], N(1))
     @test overapproximate(H, HPolyhedron, BoxDirections{N}(2)) == HPolyhedron([H])
+
+    # Example 3.1.29 from thesis
+    PZ1 = SparsePolynomialZonotope(N[-5, 0], N[2 0 2; 0 2 2], Matrix{N}(undef, 2, 0), [1 0 1; 0 1 1])
+    PZ2 = SparsePolynomialZonotope(N[3, 3], N[1 -2 2; 2 3 1], hcat(N[1//2; 0]), [1 0 2; 0 1 1])
+    PZ = overapproximate(CH(PZ1, PZ2), SparsePolynomialZonotope)
+    @test center(PZ) == N[-1, 3//2]   # no reasonable tests available here
 end
 
 for N in [Float64]