diff --git a/src/AbstractHPolygon.jl b/src/AbstractHPolygon.jl
index bae0304ab7..8f828916ce 100644
--- a/src/AbstractHPolygon.jl
+++ b/src/AbstractHPolygon.jl
@@ -185,7 +185,8 @@ end
 
 """
     rand(::Type{HPOLYGON}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
-         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
+         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing,
+         [num_constraints]::Int=-1
         )::HPOLYGON{N} where {HPOLYGON<:AbstractHPolygon}
 
 Create a random polygon in constraint representation.
@@ -198,7 +199,7 @@ Create a random polygon in constraint representation.
 - `rng`             -- (optional, default: `GLOBAL_RNG`) random number generator
 - `seed`            -- (optional, default: `nothing`) seed for reseeding
 - `num_constraints` -- (optional, default: `-1`) number of constraints of the
-                       polygon (see comment below)
+                       polygon (must be 3 or bigger; see comment below)
 
 ### Output
 
@@ -209,6 +210,8 @@ A random polygon in constraint representation.
 We create a random polygon in vertex representation and convert it to constraint
 representation.
 See [`rand(::Type{VPolygon})`](@ref).
+For non-flat polygons the number of vertices and the number of constraints are
+identical.
 """
 function rand(::Type{HPOLYGON};
               N::Type{<:Real}=Float64,
@@ -217,7 +220,9 @@ function rand(::Type{HPOLYGON};
               seed::Union{Int, Nothing}=nothing,
               num_constraints::Int=-1
              )::HPOLYGON{N} where {HPOLYGON<:AbstractHPolygon}
-    @assert dim == 2 "a polygon must have dimension 2"
+    @assert dim == 2 "cannot create a random $HPOLYGON of dimension $dim"
+    @assert num_constraints >= 3 "cannot construct a random $HPOLYGON with " *
+        "only $num_constraints constraints"
     rng = reseed(rng, seed)
     vpolygon = rand(VPolygon; N=N, dim=dim, rng=rng, seed=seed,
                     num_vertices=num_constraints)
diff --git a/src/Ellipsoid.jl b/src/Ellipsoid.jl
index 92a2b247df..620d901a0f 100644
--- a/src/Ellipsoid.jl
+++ b/src/Ellipsoid.jl
@@ -199,6 +199,12 @@ deviation 1.
 The idea for the shape matrix comes from
 [here](https://math.stackexchange.com/a/358092).
 The matrix is symmetric positive definite, but also diagonally dominant.
+
+```math
+Q =  \frac{1}{2}(S + S^T) + nI,
+```
+where ``n`` = `dim` (defaults to 2), and ``S`` is a ``n \\times n`` random
+matrix whose coefficients are uniformly distributed in the interval ``[-1, 1]``.
 """
 function rand(::Type{Ellipsoid};
               N::Type{<:Real}=Float64,
@@ -209,7 +215,18 @@ function rand(::Type{Ellipsoid};
     rng = reseed(rng, seed)
     center = randn(rng, N, dim)
     # random entries in [-1, 1]
-    shape_matrix = rand(N(-1):N(.0001):N(1), dim, dim)
+    # this needs a bit of code because 'rand' only samples from [0, 1]
+    shape_matrix = Matrix{N}(undef, dim, dim)
+    for j in 1:dim
+        for i in 1:dim
+            entry = rand(N)
+            if rand(Bool)
+                entry = -entry
+            end
+            shape_matrix[i, j] = entry
+        end
+    end
+    # make diagonally dominant
     shape_matrix = N(0.5) * (shape_matrix + shape_matrix') +
                    Matrix{N}(dim*I, dim, dim)
     return Ellipsoid(center, shape_matrix)
diff --git a/src/EmptySet.jl b/src/EmptySet.jl
index 6b45ee493d..17487d6671 100644
--- a/src/EmptySet.jl
+++ b/src/EmptySet.jl
@@ -103,7 +103,7 @@ function an_element(∅::EmptySet)
 end
 
 """
-    rand(::Type{EmptySet}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
+    rand(::Type{EmptySet}; [N]::Type{<:Real}=Float64, [dim]::Int=0,
          [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
         )::EmptySet{N}
 
diff --git a/src/Interval.jl b/src/Interval.jl
index 3e8fa1774e..2d48cb3e00 100644
--- a/src/Interval.jl
+++ b/src/Interval.jl
@@ -317,7 +317,7 @@ function rand(::Type{Interval};
               rng::AbstractRNG=GLOBAL_RNG,
               seed::Union{Int, Nothing}=nothing
              )::Interval{N}
-    @assert dim == 1 "an Interval must have dimension 1"
+    @assert dim == 1 "cannot create a random Interval of dimension $dim"
     rng = reseed(rng, seed)
     x = randn(rng, N)
     y = randn(rng, N)
diff --git a/src/Line.jl b/src/Line.jl
index 44075ba8c3..9763e24f45 100644
--- a/src/Line.jl
+++ b/src/Line.jl
@@ -207,7 +207,7 @@ function rand(::Type{Line};
               rng::AbstractRNG=GLOBAL_RNG,
               seed::Union{Int, Nothing}=nothing
              )::Line{N}
-    @assert dim == 2 "a Line must have dimension 2"
+    @assert dim == 2 "cannot create a random Line of dimension $dim"
     rng = reseed(rng, seed)
     a = randn(rng, N, dim)
     while iszero(a)
diff --git a/src/LineSegment.jl b/src/LineSegment.jl
index ad507d1b09..d9616058f2 100644
--- a/src/LineSegment.jl
+++ b/src/LineSegment.jl
@@ -245,7 +245,7 @@ function rand(::Type{LineSegment};
               rng::AbstractRNG=GLOBAL_RNG,
               seed::Union{Int, Nothing}=nothing
              )::LineSegment{N}
-    @assert dim == 2 "a LineSegment must have dimension 2"
+    @assert dim == 2 "cannot create a random LineSegment of dimension $dim"
     rng = reseed(rng, seed)
     p = randn(rng, N, dim)
     q = randn(rng, N, dim)
diff --git a/src/VPolygon.jl b/src/VPolygon.jl
index 5fc633361e..3de7b5e836 100644
--- a/src/VPolygon.jl
+++ b/src/VPolygon.jl
@@ -277,6 +277,60 @@ function ∈(x::AbstractVector{N}, P::VPolygon{N})::Bool where {N<:Real}
     return true
 end
 
+"""
+    _random_zero_sum_vector(rng::AbstractRNG, N::Type{<:Real}, n::Int)
+
+Create a random vector with entries whose sum is zero.
+
+### Input
+
+- `rng` -- random number generator
+- `N`   -- numeric type
+- `n`   -- length of vector
+
+### Output
+
+A random vector of random numbers such that all positive entries come first and
+all negative entries come last, and such that the total sum is zero.
+
+### Algorithm
+
+This is a preprocessing step of the algorithm
+[here](https://stackoverflow.com/a/47358689) based on
+[P. Valtr. Probability that n random points are in convex
+position](https://link.springer.com/article/10.1007%2FBF01271274).
+"""
+function _random_zero_sum_vector(rng::AbstractRNG, N::Type{<:Real}, n::Int)
+    # generate a sorted list of random x and y coordinates
+    list = sort!(randn(rng, N, n))
+    while (length(remove_duplicates_sorted!(list)) < n)
+        # make sure that no duplicates exist
+        list = sort!(append!(list, randn(rng, N, length(list) - n)))
+    end
+    # lists of consecutive points
+    in_l1 = rand(rng, Bool, n-2)
+    l1 = Vector{N}() # normal
+    l2 = Vector{N}() # inverted
+    push!(l1, list[1])
+    push!(l2, list[1])
+    for i in 2:n-1
+        push!(in_l1[i-1] ? l1 : l2, list[i])
+    end
+    push!(l1, list[end])
+    push!(l2, list[end])
+    # convert to vectors representing the distance (order does not matter)
+    dist = Vector{N}()
+    sizehint!(dist, n)
+    for i in 1:length(l1)-1
+        push!(dist, l1[i+1] - l1[i])
+    end
+    for i in 1:length(l2)-1
+        push!(dist, l2[i] - l2[i+1])
+    end
+    @assert isapprox(sum(dist), zero(N), atol=1e-6)
+    return dist
+end
+
 """
     rand(::Type{VPolygon}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
          [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
@@ -303,9 +357,11 @@ A random polygon in vertex representation.
 We follow the idea [here](https://stackoverflow.com/a/47358689) based on
 [P. Valtr. Probability that n random points are in convex
 position](https://link.springer.com/article/10.1007%2FBF01271274).
+There is also a nice video available
+[here](http://cglab.ca/~sander/misc/ConvexGeneration/convex.html).
 
 The number of vertices can be controlled with the argument `num_vertices`.
-For a negative value we choose a random number in the range `1:10`.
+For a negative value we choose a random number in the range `3:10`.
 """
 function rand(::Type{VPolygon};
               N::Type{<:Real}=Float64,
@@ -314,63 +370,24 @@ function rand(::Type{VPolygon};
               seed::Union{Int, Nothing}=nothing,
               num_vertices::Int=-1
              )::VPolygon{N}
-    @assert dim == 2 "a polygon must have dimension 2"
+    @assert dim == 2 "cannot create a random VPolygon of dimension $dim"
     rng = reseed(rng, seed)
     if num_vertices < 0
-        num_vertices = rand(1:10)
+        num_vertices = rand(3:10)
     end
 
     # special cases, 0 or 1 vertex
-    if num_vertices <= 1
-        vertices = [randn(rng, N, 2) for i in 1:num_vertices]
-        return VPolygon(vertices; apply_convex_hull=false)
+    if num_vertices == 0
+        return VPolygon{N}()
+    elseif num_vertices == 1
+        return VPolygon([randn(rng, N, 2)])
     end
 
     # general case, >= 2 vertices
 
-    function remove_duplicates_sorted!(v::AbstractVector)
-        for i in length(v)-1:-1:1
-            if v[i] == v[i+1]
-                splice!(v, i+1)
-            end
-        end
-        return v
-    end
-
-    function random_axis_aligned_vectors(rng, N, n)
-        # generate a sorted list of random x and y coordinates
-        list = sort!(randn(rng, N, n))
-        while (length(remove_duplicates_sorted!(list)) < n)
-            # make sure that no duplicates exist
-            list = sort!(append!(list, randn(rng, N, length(list) - n)))
-        end
-        # lists of consecutive points
-        in_l1 = rand(rng, Bool, n-2)
-        l1 = Vector{N}() # normal
-        l2 = Vector{N}() # inverted
-        push!(l1, list[1])
-        push!(l2, list[1])
-        for i in 2:n-1
-            push!(in_l1[i-1] ? l1 : l2, list[i])
-        end
-        push!(l1, list[end])
-        push!(l2, list[end])
-        # convert to vectors representing the distance (order does not matter)
-        dist = Vector{N}()
-        sizehint!(dist, n)
-        for i in 1:length(l1)-1
-            push!(dist, l1[i+1] - l1[i])
-        end
-        for i in 1:length(l2)-1
-            push!(dist, l2[i] - l2[i+1])
-        end
-        @assert isapprox(sum(dist), zero(N), atol=1e-6)
-        return dist
-    end
-
     # get random horizontal and vertical vectors
-    horiz = random_axis_aligned_vectors(rng, N, num_vertices)
-    vert = random_axis_aligned_vectors(rng, N, num_vertices)
+    horiz = _random_zero_sum_vector(rng, N, num_vertices)
+    vert = _random_zero_sum_vector(rng, N, num_vertices)
 
     # randomly combine horizontal and vertical vectors
     m = num_vertices
diff --git a/src/helper_functions.jl b/src/helper_functions.jl
index eb4704d7f6..35287d286a 100644
--- a/src/helper_functions.jl
+++ b/src/helper_functions.jl
@@ -86,6 +86,28 @@ function ispermutation(u::AbstractVector{T}, v::AbstractVector{T})::Bool where T
     return true
 end
 
+"""
+    remove_duplicates_sorted!(v::AbstractVector)
+
+Remove duplicate entries in a sorted vector.
+
+### Input
+
+- `v` -- sorted vector
+
+### Output
+
+The input vector without duplicates.
+"""
+function remove_duplicates_sorted!(v::AbstractVector)
+    for i in length(v)-1:-1:1
+        if v[i] == v[i+1]
+            splice!(v, i+1)
+        end
+    end
+    return v
+end
+
 """
     samedir(u::AbstractVector{N},
             v::AbstractVector{N})::Tuple{Bool, Real} where N<:Real