diff --git a/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md b/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md
index 08e911c465..2088356d9b 100644
--- a/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md
+++ b/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md
@@ -7,6 +7,7 @@ CurrentModule = LazySets
 ```@docs
 SimpleSparsePolynomialZonotope
 PolynomialZonotope
+rand(::Type{SimpleSparsePolynomialZonotope})
 center(::SimpleSparsePolynomialZonotope)
 genmat(::SimpleSparsePolynomialZonotope)
 expmat(::SimpleSparsePolynomialZonotope)
diff --git a/src/Sets/SimpleSparsePolynomialZonotope.jl b/src/Sets/SimpleSparsePolynomialZonotope.jl
index 33ec4fa54e..28f5093e6b 100644
--- a/src/Sets/SimpleSparsePolynomialZonotope.jl
+++ b/src/Sets/SimpleSparsePolynomialZonotope.jl
@@ -311,3 +311,56 @@ function remove_redundant_generators(S::SimpleSparsePolynomialZonotope)
 
     return SimpleSparsePolynomialZonotope(cnew, Gnew, Enew)
 end
+
+
+"""
+    rand(::Type{SimpleSparsePolynomialZonotope};
+         [N]::Type{<:Real}=Float64, [dim]::Int=2, [nparams]::Int=2, [maxdeg]::Int=3,
+         [num_generators]::Int=-1,
+         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)
+
+Create a random zonotope.
+
+### Input
+
+- `Zonotope`       -- type for dispatch
+- `N`              -- (optional, default: `Float64`) numeric type
+- `dim`            -- (optional, default: 2) dimension
+- `nparams`        -- (optional, default: 2) number of parameters
+- `maxdeg`         -- (optinal, default: 3) maximum degree for each parameter
+- `rng`            -- (optional, default: `GLOBAL_RNG`) random number generator
+- `seed`           -- (optional, default: `nothing`) seed for reseeding
+- `num_generators` -- (optional, default: `-1`) number of generators of the
+                      zonotope (see comment below)
+
+### Output
+
+A random zonotope.
+
+### Algorithm
+
+All numbers are normally distributed with mean 0 and standard deviation 1.
+
+The number of generators can be controlled with the argument `num_generators`.
+For a negative value we choose a random number in the range `dim:2*dim` (except
+if `dim == 1`, in which case we only create a single generator). Note that the final number of
+generators may be lower if redundant monomials are generated.
+"""
+function rand(::Type{SimpleSparsePolynomialZonotope};
+              N::Type{<:Real}=Float64,
+              dim::Int=2,
+              nparams::Int=2,
+              maxdeg::Int=3,
+              rng::AbstractRNG=GLOBAL_RNG,
+              seed::Union{Int, Nothing}=nothing,
+              num_generators::Int=-1)
+    rng = reseed(rng, seed)
+    center = randn(rng, N, dim)
+    if num_generators < 0
+        num_generators = (dim == 1) ? 1 : rand(dim:2*dim)
+    end
+    generators = randn(rng, N, dim, num_generators)
+    expmat = rand(0:maxdeg, nparams, num_generators)
+    SSPZ = SimpleSparsePolynomialZonotope(center, generators, expmat)
+    return remove_redundant_generators(SSPZ)
+end
diff --git a/test/Sets/SimpleSparsePolynomialZonotope.jl b/test/Sets/SimpleSparsePolynomialZonotope.jl
index 91cd54854a..a70822b8ed 100644
--- a/test/Sets/SimpleSparsePolynomialZonotope.jl
+++ b/test/Sets/SimpleSparsePolynomialZonotope.jl
@@ -4,6 +4,8 @@ for N in [Float64, Float32, Rational{Int}]
     G = N[1 2;2 2.]
     E = [1 4;1 2]
 
+    @test rand(SimpleSparsePolynomialZonotope) isa SimpleSparsePolynomialZonotope
+
     S = SimpleSparsePolynomialZonotope(c, G, E)
 
     @test genmat(S) == G