diff --git a/docs/src/lib/utils.md b/docs/src/lib/utils.md
index 7dfbb090f5..58dac114d5 100644
--- a/docs/src/lib/utils.md
+++ b/docs/src/lib/utils.md
@@ -39,6 +39,7 @@ LazySets.Arrays.allequal
 LazySets.Arrays.distance(::AbstractVector, ::AbstractVector, ::Real=2.0)
 LazySets.Arrays.same_sign
 LazySets.Arrays.to_matrix
+LazySets.Arrays._rationalize
 ```
 
 ## Functions and Macros
diff --git a/src/Initialization/init_Symbolics.jl b/src/Initialization/init_Symbolics.jl
index 5edd2e6dd1..c784807f75 100644
--- a/src/Initialization/init_Symbolics.jl
+++ b/src/Initialization/init_Symbolics.jl
@@ -41,7 +41,11 @@ eval(quote
     # case with a single variable
     _vec(vars::Tuple{Num}) = [vars[1]]
 
-    _get_variables(expr::Num) = Symbolics.get_variables(expr)
+    # reduce for several variables e.g. when vars = @variables x[1:3] t
+    _vec(vars::Vector{Any}) = reduce(vcat, vars)
+    _vec(vars::Vector{Num}) = vars
+
+    _get_variables(expr::Num) = convert(Vector{Num}, Symbolics.get_variables(expr))
     _get_variables(expr::Vector{<:Num}) = unique(reduce(vcat, _get_variables(ex) for ex in expr))
 end)
 
diff --git a/src/Sets/HPolyhedron.jl b/src/Sets/HPolyhedron.jl
index 85ac0f80f3..680f59a0a5 100644
--- a/src/Sets/HPolyhedron.jl
+++ b/src/Sets/HPolyhedron.jl
@@ -705,7 +705,7 @@ HPolyhedron{Float64,Array{Float64,1}}(HalfSpace{Float64,Array{Float64,1}}[HalfSp
 ace{Float64,Array{Float64,1}}([-1.0, -0.0], 0.0), HalfSpace{Float64,Array{Float64,1}}([0.0, 1.0], -0.0)])
 ```
 """
-function HPolyhedron(expr::Vector{<:Num}, vars=_get_variables(expr); N::Type{<:Real}=Float64)
+function HPolyhedron(expr::Vector{<:Num}, vars::Vector{Num}; N::Type{<:Real}=Float64)
     clist = Vector{HalfSpace{N, Vector{N}}}()
     sizehint!(clist, length(expr))
     got_hyperplane = false
@@ -734,4 +734,7 @@ function HPolyhedron(expr::Vector{<:Num}, vars=_get_variables(expr); N::Type{<:R
     return HPolyhedron(clist)
 end
 
+HPolyhedron(expr::Vector{<:Num}; N::Type{<:Real}=Float64) = HPolyhedron(expr, _get_variables(expr); N=N)
+HPolyhedron(expr::Vector{<:Num}, vars; N::Type{<:Real}=Float64) = HPolyhedron(expr, _vec(vars); N=N)
+
 end end  # quote / load_modeling_toolkit_hpolyhedron()
diff --git a/src/Sets/HPolytope.jl b/src/Sets/HPolytope.jl
index 99c6de2520..62321ad087 100644
--- a/src/Sets/HPolytope.jl
+++ b/src/Sets/HPolytope.jl
@@ -282,9 +282,11 @@ HalfSpace{Float64,Array{Float64,1}}([-1.0, 0.0], 0.0), HalfSpace{Float64,Array{F
 HalfSpace{Float64,Array{Float64,1}}([0.0, -1.0], 0.0)])
 ```
 """
-function HPolytope(expr::Vector{<:Num}, vars=_get_variables(expr);
-                   N::Type{<:Real}=Float64, check_boundedness::Bool=false)
+function HPolytope(expr::Vector{<:Num}, vars::Vector{Num}; N::Type{<:Real}=Float64, check_boundedness::Bool=false)
     return HPolytope([HalfSpace(ex, vars; N=N) for ex in expr], check_boundedness=check_boundedness)
 end
 
+HPolytope(expr::Vector{<:Num}; N::Type{<:Real}=Float64, check_boundedness::Bool=false) = HPolytope(expr, _get_variables(expr); N=N, check_boundedness=check_boundedness)
+HPolytope(expr::Vector{<:Num}, vars; N::Type{<:Real}=Float64, check_boundedness::Bool=false) = HPolytope(expr, _vec(vars); N=N, check_boundedness=check_boundedness)
+
 end end  # quote / load_modeling_toolkit_hpolytope()
diff --git a/src/Sets/HalfSpace.jl b/src/Sets/HalfSpace.jl
index 40c22ed408..c10e070bf4 100644
--- a/src/Sets/HalfSpace.jl
+++ b/src/Sets/HalfSpace.jl
@@ -639,7 +639,7 @@ Note in particular that strict inequalities are relaxed as being smaller-or-equa
 Finally, the returned set is the half-space with normal vector `[a1, …, an]` and
 displacement `b`.
 """
-function HalfSpace(expr::Num, vars=_get_variables(expr); N::Type{<:Real}=Float64)
+function HalfSpace(expr::Num, vars::Vector{Num}; N::Type{<:Real}=Float64)
     valid, sexpr = _is_halfspace(Symbolics.value(expr))
     if !valid
         throw(ArgumentError("expected an expression describing a half-space, got $expr"))
@@ -655,10 +655,8 @@ function HalfSpace(expr::Num, vars=_get_variables(expr); N::Type{<:Real}=Float64
     return HalfSpace(coeffs, β)
 end
 
-function HalfSpace(expr::Num, vars::NTuple{L, Union{<:Num, <:Vector{Num}}}; N::Type{<:Real}=Float64) where {L}
-    vars = _vec(vars)
-    return HalfSpace(expr, vars, N=N)
-end
+HalfSpace(expr::Num; N::Type{<:Real}=Float64) = HalfSpace(expr, _get_variables(expr); N=N)
+HalfSpace(expr::Num, vars; N::Type{<:Real}=Float64) = HalfSpace(expr, _vec(vars), N=N)
 
 end end  # quote / load_modeling_toolkit_halfspace()
 
diff --git a/src/Sets/Hyperplane.jl b/src/Sets/Hyperplane.jl
index 2e777f00b3..f0d37bfa6a 100644
--- a/src/Sets/Hyperplane.jl
+++ b/src/Sets/Hyperplane.jl
@@ -556,7 +556,7 @@ Therefore, the order in which the variables appear in `vars` affects the final r
 Finally, the returned set is the hyperplane with normal vector `[a1, …, an]` and
 displacement `b`.
 """
-function Hyperplane(expr::Num, vars=_get_variables(expr); N::Type{<:Real}=Float64)
+function Hyperplane(expr::Num, vars::Vector{Num}; N::Type{<:Real}=Float64)
     valid, sexpr = _is_hyperplane(Symbolics.value(expr))
     if !valid
         throw(ArgumentError("expected an expression of the form `ax == b`, got $expr"))
@@ -572,9 +572,7 @@ function Hyperplane(expr::Num, vars=_get_variables(expr); N::Type{<:Real}=Float6
     return Hyperplane(coeffs, β)
 end
 
-function Hyperplane(expr::Num, vars::NTuple{L, Union{<:Num, <:Vector{Num}}}; N::Type{<:Real}=Float64) where {L}
-    vars = _vec(vars)
-    return Hyperplane(expr, vars, N=N)
-end
+Hyperplane(expr::Num; N::Type{<:Real}=Float64) = Hyperplane(expr, _get_variables(expr); N=N)
+Hyperplane(expr::Num, vars; N::Type{<:Real}=Float64) = Hyperplane(expr, _vec(vars), N=N)
 
 end end  # quote / load_symbolics_hyperplane()