diff --git a/docs/src/lib/sets/SparsePolynomialZonotope.md b/docs/src/lib/sets/SparsePolynomialZonotope.md
index 022d1ceffb..4dd9360364 100644
--- a/docs/src/lib/sets/SparsePolynomialZonotope.md
+++ b/docs/src/lib/sets/SparsePolynomialZonotope.md
@@ -19,4 +19,5 @@ nindependentgens(::SparsePolynomialZonotope)
 nparams(::SparsePolynomialZonotope)
 order(::SparsePolynomialZonotope)
 linear_map(::AbstractMatrix, ::SparsePolynomialZonotope)
+exact_sum(::SparsePolynomialZonotope, ::SparsePolynomialZonotope)
 ```
diff --git a/src/Sets/SparsePolynomialZonotope.jl b/src/Sets/SparsePolynomialZonotope.jl
index be303eefdd..f2cc9941e1 100644
--- a/src/Sets/SparsePolynomialZonotope.jl
+++ b/src/Sets/SparsePolynomialZonotope.jl
@@ -1,5 +1,5 @@
 export SparsePolynomialZonotope, expmat, nparams, ndependentgens, nindependentgens,
-       dependent_genmat, independent_genmat, indexvector,
+       dependent_genmat, independent_genmat, indexvector, exact_sum, ⊞,
        linear_map, quadratic_map, remove_redundant_generators
 
 """
@@ -288,6 +288,37 @@ function rand(::Type{SparsePolynomialZonotope};
     return SparsePolynomialZonotope(center(SSPZ), genmat(SSPZ), GI, expmat(SSPZ))
 end
 
+
+"""
+    exact_sum(P1::SparsePolynomialZonotope, P2::SparsePolynomialZonotope)
+
+Compute the exact sum of sparse polyomial zonotopes ``P₁`` and ``P₂``.
+
+### Input
+
+- `P1` -- sparse polynomial zonotope
+- `P2` -- sparse polynomial zonotope
+
+### Output
+
+exact sum ``P₁ ⊞ P₂``
+
+"""
+function exact_sum(P1::SparsePolynomialZonotope, P2::SparsePolynomialZonotope)
+
+    indexvector(P1) == indexvector(P2) || throw(ArgumentError("Exact sum currently only implemented for Sparse Polynomial Zonotopes with the same index vector"))
+
+    c = center(P1) + center(P2)
+    G = hcat(dependent_genmat(P1), dependent_genmat(P2))
+    GI = hcat(independent_genmat(P1), independent_genmat(P2))
+    E = hcat(expmat(P1), expmat(P2))
+    idx = indexvector(P1)
+
+    return SparsePolynomialZonotope(c, G, GI, E, idx)
+end
+
+const ⊞ = exact_sum
+
 # function quadratic_map(Q::Vector{MT}, S::SparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}}
 #     m = length(Q)
 #     c = center(S)
diff --git a/test/Sets/SparsePolynomialZonotope.jl b/test/Sets/SparsePolynomialZonotope.jl
index 8ce5256671..649f73b633 100644
--- a/test/Sets/SparsePolynomialZonotope.jl
+++ b/test/Sets/SparsePolynomialZonotope.jl
@@ -22,13 +22,20 @@ for N in [Float64, Float32, Rational{Int}]
     @test order(PZ) == 2//1
 
     M = N[-0.5 0.2;-0.1 0.6]
-    PZ = SparsePolynomialZonotope(zeros(N, 2), N[2 0 1;1 2 1], zeros(N, 2, 0), [1 0 1;0 1 3])
-    LMPZ = linear_map(M, PZ)
-    @test center(PZ) == zeros(N, 2)
+    PZ2 = SparsePolynomialZonotope(zeros(N, 2), N[2 0 1;1 2 1], zeros(N, 2, 0), [1 0 1;0 1 3])
+    LMPZ = linear_map(M, PZ2)
+    @test center(LMPZ) == zeros(N, 2)
     @test dependent_genmat(LMPZ) ≈ N[-0.8 0.4 -0.3;0.4 1.2 0.5] atol=1e-7
     @test isempty(independent_genmat(LMPZ))
     @test expmat(LMPZ) == [1 0 1;0 1 3]
     @test indexvector(LMPZ) == indexvector(PZ)
+
+    ESPZ = PZ ⊞ PZ2
+    @test center(ESPZ) == [4, 4]
+    @test dependent_genmat(ESPZ) == [2 1 2 2 0 1;0 2 2 1 2 1]
+    @test independent_genmat(ESPZ) == hcat([1, 0])
+    @test expmat(ESPZ) == [1 0 3 1 0 1;0 1 1 0 1 3]
+    @test indexvector(ESPZ) == indexvector(PZ)
 end
 
 SSPZ = SimpleSparsePolynomialZonotope([0.2, -0.6], [1 0;0 0.4], [1 0;0 1])