diff --git a/src/Sets/HParallelotope.jl b/src/Sets/HParallelotope.jl
index 82a2c198d3..9d48acefba 100644
--- a/src/Sets/HParallelotope.jl
+++ b/src/Sets/HParallelotope.jl
@@ -344,6 +344,15 @@ A random parallelotope.
 ### Notes
 
 All numbers are normally distributed with mean 0 and standard deviation 1.
+
+### Algorithm
+
+The directions matrix and offset vector are created randomly. On average there
+is a good chance that this resulting set is empty. We then modify the offset to
+ensure non-emptiness.
+
+There is a chance that the resulting set represents an unbounded set. This
+implementation checks for that case and then samples a new set.
 """
 function rand(::Type{HParallelotope};
               N::Type{<:Real}=Float64,
@@ -351,9 +360,26 @@ function rand(::Type{HParallelotope};
               rng::AbstractRNG=GLOBAL_RNG,
               seed::Union{Int, Nothing}=nothing)
     rng = reseed(rng, seed)
-    D = randn(rng, N, dim, dim)
-    offset = randn(rng, N, 2 * dim)
-    return HParallelotope(D, offset)
+
+    while true
+        D = randn(rng, N, dim, dim)
+        offset = randn(rng, N, 2 * dim)
+
+        # make sure that the set is not empty:
+        # `offset[i] >= -offset[i-dim]` for all `i ∈ dim+1:2*dim`
+        @inbounds for i in dim+1:2*dim
+            offset[i] = -offset[i-dim] + abs(offset[i])
+        end
+
+        P = HParallelotope(D, offset)
+
+        # convert to polyhedron to check boundedness
+        Q = HPolyhedron(vcat(P.directions, -P.directions), P.offset)
+        if isbounded(Q)
+            return P
+        end
+        # set is unbounded; sample a new set in the next iteration
+    end
 end
 
 """