diff --git a/docs/src/man/reach_zonotopes_hybrid.md b/docs/src/man/reach_zonotopes_hybrid.md
index de63d51281..939a302af3 100644
--- a/docs/src/man/reach_zonotopes_hybrid.md
+++ b/docs/src/man/reach_zonotopes_hybrid.md
@@ -49,7 +49,7 @@ function reach_hybrid(As, Ts, init, δ, μ, T, max_order, instant_transitions)
         println("currently in location $loc at time $t")
         R = reach_continuous(As[loc], init, δ, μ, T-t, max_order)
         found_transition = false
-        for i in 1:length(R)-1
+        for i in 1:(length(R)-1)
             S = R[i]
             push!(res, (S, loc))
             for (guard, tgt_loc) in Ts[loc]
@@ -130,7 +130,7 @@ The following function does that for 2-mode models.
 ```@example example_reach_zonotopes_hybrid
 function plot_res(res)
     p = plot()
-    for i in 1:length(res)
+    for i in eachindex(res)
         if res[i][2] == 1
             c = "blue"
         elseif res[i][2] == 2
diff --git a/src/ConcreteOperations/convex_hull.jl b/src/ConcreteOperations/convex_hull.jl
index f38f62e869..4cb06c0488 100644
--- a/src/ConcreteOperations/convex_hull.jl
+++ b/src/ConcreteOperations/convex_hull.jl
@@ -448,7 +448,7 @@ function monotone_chain!(points::Vector{VN}; sort::Bool=true) where {N,VN<:Abstr
 
     # build lower hull
     lower = Vector{VN}()
-    iterator = 1:length(points)
+    iterator = eachindex(points)
     build_hull!(lower, iterator, points)
 
     # build upper hull
diff --git a/src/ConcreteOperations/intersection.jl b/src/ConcreteOperations/intersection.jl
index a3284ffec0..a2a2f92e6f 100644
--- a/src/ConcreteOperations/intersection.jl
+++ b/src/ConcreteOperations/intersection.jl
@@ -553,7 +553,7 @@ function intersection(P1::AbstractHPolygon, P2::AbstractHPolygon;
     i1 = 1
     i2 = 1
     duplicates = 0
-    for i in 1:length(c)
+    for i in eachindex(c)
         if c1[i1].a <= c2[i2].a
             if c2[i2].a <= c1[i1].a
                 duplicates += 1
diff --git a/src/ConcreteOperations/isdisjoint.jl b/src/ConcreteOperations/isdisjoint.jl
index dfe64de08b..8c7e323150 100644
--- a/src/ConcreteOperations/isdisjoint.jl
+++ b/src/ConcreteOperations/isdisjoint.jl
@@ -742,7 +742,7 @@ function isdisjoint(H1::HalfSpace, H2::HalfSpace, witness::Bool=false)
     end
     # compute witness
     v = zeros(N, length(a1))
-    for i in 1:length(a1)
+    for i in eachindex(a1)
         a_sum_i = a1[i] + a2[i]
         if a_sum[i] != 0
             v[i] = (H1.b + H2.b) / a_sum_i
@@ -1100,7 +1100,7 @@ function isdisjoint(X::CartesianProductArray, Y::CartesianProductArray,
     @assert same_block_structure(array(X), array(Y)) "block structure has to " *
                                                      "be identical"
 
-    for i in 1:length(X.array)
+    for i in eachindex(X.array)
         if isdisjoint(X.array[i], Y.array[i])
             return _witness_result_empty(witness, true, X, Y)
         end
diff --git a/src/ConcreteOperations/issubset.jl b/src/ConcreteOperations/issubset.jl
index d540f7c0fc..5c5955ef4e 100644
--- a/src/ConcreteOperations/issubset.jl
+++ b/src/ConcreteOperations/issubset.jl
@@ -1210,7 +1210,7 @@ function ⊆(X::CartesianProductArray, Y::CartesianProductArray,
     end
 
     N = promote_type(eltype(X), eltype(Y))
-    for i in 1:length(aX)
+    for i in eachindex(aX)
         result = ⊆(aX[i], aY[i], witness)
         if !witness && !result
             return false
@@ -1218,7 +1218,7 @@ function ⊆(X::CartesianProductArray, Y::CartesianProductArray,
             # construct a witness
             w = Vector{N}(undef, dim(X))
             k = 1
-            for j in 1:length(aX)
+            for j in eachindex(aX)
                 Xj = aX[j]
                 l = k + dim(Xj)
                 w[k:(l - 1)] = j == i ? result[2] : an_element(Xj)
diff --git a/src/Interfaces/AbstractHyperrectangle.jl b/src/Interfaces/AbstractHyperrectangle.jl
index 30d3ef7f01..af37ddfc84 100644
--- a/src/Interfaces/AbstractHyperrectangle.jl
+++ b/src/Interfaces/AbstractHyperrectangle.jl
@@ -235,7 +235,7 @@ function vertices_list(H::AbstractHyperrectangle; kwargs...)
     vlist = Vector{typeof(c)}(undef, m)
     vlist[1] = copy(v)
     @inbounds for i in 2:m
-        for j in 1:length(v)
+        for j in eachindex(v)
             if trivector[j] == Int8(-1)
                 trivector[j] = Int8(1)
                 v[j] = c[j] + radius_hyperrectangle(H, j)
@@ -750,7 +750,7 @@ A scalar representing the distance between point `x` and hyperrectangle `H`.
     # compute closest point
     y = similar(x)
     outside = false
-    @inbounds for i in 1:length(x)
+    @inbounds for i in eachindex(x)
         ci = center(H, i)
         ri = radius_hyperrectangle(H, i)
         d = x[i] - ci
diff --git a/src/Interfaces/AbstractPolyhedron_functions.jl b/src/Interfaces/AbstractPolyhedron_functions.jl
index 49d285fbe7..ea228d2817 100644
--- a/src/Interfaces/AbstractPolyhedron_functions.jl
+++ b/src/Interfaces/AbstractPolyhedron_functions.jl
@@ -739,7 +739,7 @@ function _affine_map_inverse_hrep(A::AbstractMatrix, P::LazySet,
     end
     if has_undefs  # there were redundant constraints, so remove them
         constraints_res = [constraints_res[i]
-                           for i in 1:length(constraints_res)
+                           for i in eachindex(constraints_res)
                            if isassigned(constraints_res, i)]
     end
     return constraints_res
diff --git a/src/LazyOperations/CartesianProductArray.jl b/src/LazyOperations/CartesianProductArray.jl
index 1298d42a50..dd6b0e4bb9 100644
--- a/src/LazyOperations/CartesianProductArray.jl
+++ b/src/LazyOperations/CartesianProductArray.jl
@@ -478,7 +478,7 @@ function same_block_structure(x::AbstractVector{S1},
     if length(x) != length(y)
         return false
     end
-    for i in 1:length(x)
+    for i in eachindex(x)
         if dim(x[i]) != dim(y[i])
             return false
         end
@@ -577,7 +577,7 @@ function block_to_dimension_indices(cpa::CartesianProductArray{N},
     constrained_blocks = 0
     start_index, end_index = 1, 0
     v_i = 1
-    @inbounds for i in 1:length(cpa.array)
+    @inbounds for i in eachindex(cpa.array)
         end_index += dim(cpa.array[i])
         if v_i <= length(vars) && vars[v_i] <= end_index
             ranges[i] = (start_index, end_index)
@@ -599,7 +599,7 @@ function block_to_dimension_indices(cpa::CartesianProductArray{N}) where {N}
     ranges = Vector{Tuple{Int,Int}}(undef, length(cpa.array))
 
     start_index, end_index = 1, 0
-    @inbounds for i in 1:length(cpa.array)
+    @inbounds for i in eachindex(cpa.array)
         end_index += dim(cpa.array[i])
         ranges[i] = (start_index, end_index)
         start_index = end_index + 1
@@ -633,8 +633,8 @@ function substitute_blocks(low_dim_cpa::CartesianProductArray{N},
                            blocks::Vector{Tuple{Int,Int}}) where {N}
     array = Vector{LazySet{N}}(undef, length(orig_cpa.array))
     index = 1
-    for bi in 1:length(orig_cpa.array)
-        start_ind, end_index = blocks[bi]
+    for bi in eachindex(orig_cpa.array)
+        start_ind, _ = blocks[bi]
         if start_ind == -1
             array[bi] = orig_cpa.array[bi]
         else
diff --git a/src/LazyOperations/Intersection.jl b/src/LazyOperations/Intersection.jl
index 5849170d7b..a083c83e76 100644
--- a/src/LazyOperations/Intersection.jl
+++ b/src/LazyOperations/Intersection.jl
@@ -984,7 +984,7 @@ function get_constrained_lowdimset(cpa::CartesianProductArray{N,S},
     sizehint!(block_structure, non_empty_length)
 
     last_var = 1
-    for i in 1:length(blocks)
+    for i in eachindex(blocks)
         start_index, end_index = blocks[i]
         block_end = last_var + end_index - start_index
         if start_index != -1
diff --git a/src/Plotting/plot_recipes.jl b/src/Plotting/plot_recipes.jl
index d35230a146..23bec24a14 100644
--- a/src/Plotting/plot_recipes.jl
+++ b/src/Plotting/plot_recipes.jl
@@ -161,7 +161,7 @@ julia> plot(Bs, 1e-3)  # default accuracy value (explicitly given for clarity)
 julia> plot(Bs, 1e-2)  # faster but less accurate than the previous call
 ```
 """
-@recipe function plot_list(list::AbstractVector{VN}, ε::Real=N(PLOT_PRECISION),
+@recipe function plot_list(list::AbstractVector{VN}, ε::Real=N(PLOT_PRECISION),  # COV_EXCL_LINE
                            Nφ::Int=PLOT_POLAR_DIRECTIONS;
                            same_recipe=false) where {N,VN<:LazySet{N}}
     if same_recipe
@@ -227,7 +227,7 @@ function _plot_list_same_recipe(list::AbstractVector{VN}, ε::Real=N(PLOT_PRECIS
 end
 
 # recipe for vector of singletons
-@recipe function plot_list(list::AbstractVector{SN}) where {N,SN<:AbstractSingleton{N}}
+@recipe function plot_list(list::AbstractVector{SN}) where {N,SN<:AbstractSingleton{N}}  # COV_EXCL_LINE
     label --> DEFAULT_LABEL
     grid --> DEFAULT_GRID
     if DEFAULT_ASPECT_RATIO != :none
@@ -241,7 +241,7 @@ end
 end
 
 # plot recipe for the union of singletons
-@recipe function plot_list(X::UnionSetArray{N,SN}) where {N,SN<:AbstractSingleton{N}}
+@recipe function plot_list(X::UnionSetArray{N,SN}) where {N,SN<:AbstractSingleton{N}}  # COV_EXCL_LINE
     label --> DEFAULT_LABEL
     grid --> DEFAULT_GRID
     if DEFAULT_ASPECT_RATIO != :none
@@ -321,7 +321,7 @@ julia> plot(B, 1e-3)  # default accuracy value (explicitly given for clarity her
 julia> plot(B, 1e-2)  # faster but less accurate than the previous call
 ```
 """
-@recipe function plot_lazyset(X::LazySet{N}, ε::Real=N(PLOT_PRECISION)) where {N}
+@recipe function plot_lazyset(X::LazySet{N}, ε::Real=N(PLOT_PRECISION)) where {N}  # COV_EXCL_LINE
     label --> DEFAULT_LABEL
     grid --> DEFAULT_GRID
     if DEFAULT_ASPECT_RATIO != :none
@@ -416,7 +416,7 @@ Plot a singleton.
 julia> plot(Singleton([0.5, 1.0]))
 ```
 """
-@recipe function plot_singleton(S::AbstractSingleton{N}, ε::Real=zero(N)) where {N}
+@recipe function plot_singleton(S::AbstractSingleton{N}, ε::Real=zero(N)) where {N}  # COV_EXCL_LINE
     label --> DEFAULT_LABEL
     grid --> DEFAULT_GRID
     if DEFAULT_ASPECT_RATIO != :none
@@ -448,7 +448,7 @@ Plot an empty set.
 - `∅` -- empty set
 - `ε` -- (optional, default: `0`) ignored, used for dispatch
 """
-@recipe function plot_emptyset(∅::EmptySet{N}, ε::Real=zero(N)) where {N}
+@recipe function plot_emptyset(∅::EmptySet{N}, ::Real=zero(N)) where {N}  # COV_EXCL_LINE
     label --> DEFAULT_LABEL
     grid --> DEFAULT_GRID
     if DEFAULT_ASPECT_RATIO != :none
@@ -500,9 +500,8 @@ julia> plot(overapproximate(X, PolarDirections(100)))
 julia> plot(X, 0.0, 100)  # equivalent to the above line
 ```
 """
-@recipe function plot_intersection(cap::Intersection{N},
-                                   ε::Real=zero(N),
-                                   Nφ::Int=PLOT_POLAR_DIRECTIONS) where {N}
+@recipe function plot_intersection(cap::Intersection{N}, ε::Real=zero(N),  # COV_EXCL_LINE
+                                   ::Int=PLOT_POLAR_DIRECTIONS) where {N}
     label --> DEFAULT_LABEL
     grid --> DEFAULT_GRID
     if DEFAULT_ASPECT_RATIO != :none
@@ -544,9 +543,8 @@ end
 
 # non-convex sets
 
-@recipe function plot_union(cup::Union{UnionSet{N},UnionSetArray{N}},
-                            ε::Real=N(PLOT_PRECISION); same_recipe=false,
-                            Nφ=PLOT_POLAR_DIRECTIONS) where {N}
+@recipe function plot_union(cup::Union{UnionSet{N},UnionSetArray{N}}, ε::Real=N(PLOT_PRECISION);  # COV_EXCL_LINE
+                            same_recipe=false, Nφ=PLOT_POLAR_DIRECTIONS) where {N}
     n = dim(cup)
     # extract limits and extrema of already plotted sets
     p = plotattributes[:plot_object]
@@ -593,9 +591,8 @@ end
     end
 end
 
-@recipe function plot_polyzono(P::AbstractPolynomialZonotope{N},
-                               ε::Real=N(PLOT_PRECISION); nsdiv=10,
-                               partition=nothing) where {N}
+@recipe function plot_polyzono(P::AbstractPolynomialZonotope{N}, ε::Real=N(PLOT_PRECISION);  # COV_EXCL_LINE
+                               nsdiv=10, partition=nothing) where {N}
     label --> DEFAULT_LABEL
     grid --> DEFAULT_GRID
     if DEFAULT_ASPECT_RATIO != :none
diff --git a/src/Sets/HPolyhedron.jl b/src/Sets/HPolyhedron.jl
index 1ff117f7e1..c1e1f2b1f2 100644
--- a/src/Sets/HPolyhedron.jl
+++ b/src/Sets/HPolyhedron.jl
@@ -238,7 +238,7 @@ function rand(::Type{HPolyhedron};
     P = rand(HPolytope; N=N, dim=dim, rng=rng)
     constraints_P = constraints_list(P)
     constraints_Q = Vector{eltype(constraints_P)}()
-    for i in 1:length(constraints_P)
+    for i in eachindex(constraints_P)
         if rand(rng, Bool)
             push!(constraints_Q, constraints_P[i])
         end
diff --git a/src/Sets/Hyperplane.jl b/src/Sets/Hyperplane.jl
index eab01fda38..22ea5d04fc 100644
--- a/src/Sets/Hyperplane.jl
+++ b/src/Sets/Hyperplane.jl
@@ -369,7 +369,7 @@ to ``d``.
         # not the zero vector, check if it is a normal vector
         N = promote_type(eltype(d), eltype(a))
         factor = zero(N)
-        for i in 1:length(a)
+        for i in eachindex(a)
             if a[i] == 0
                 if d[i] != 0
                     unbounded = true
diff --git a/src/Sets/VPolygon.jl b/src/Sets/VPolygon.jl
index e8a0eca5d6..4386d27036 100644
--- a/src/Sets/VPolygon.jl
+++ b/src/Sets/VPolygon.jl
@@ -94,7 +94,7 @@ function VPolygon(vertices_matrix::MT; apply_convex_hull::Bool=true,
     @assert size(vertices_matrix, 1) == 2 "the number of rows of the matrix " *
                                           "of vertices should be 2, but it is $(size(vertices_matrix, 1))"
 
-    vertices = [vertices_matrix[:, j] for j in 1:size(vertices_matrix, 2)]
+    vertices = [vertices_matrix[:, j] for j in axes(vertices_matrix, 2)]
     return VPolygon(vertices; apply_convex_hull=apply_convex_hull,
                     algorithm=algorithm)
 end
diff --git a/src/Sets/VPolytope.jl b/src/Sets/VPolytope.jl
index eaffefbb8f..215e498095 100644
--- a/src/Sets/VPolytope.jl
+++ b/src/Sets/VPolytope.jl
@@ -67,7 +67,7 @@ VPolytope() = VPolytope{Float64}()
 
 # constructor from rectangular matrix
 function VPolytope(vertices_matrix::MT) where {N,MT<:AbstractMatrix{N}}
-    vertices = [vertices_matrix[:, j] for j in 1:size(vertices_matrix, 2)]
+    vertices = [vertices_matrix[:, j] for j in axes(vertices_matrix, 2)]
     return VPolytope(vertices)
 end
 
diff --git a/test/LazyOperations/UnionSet.jl b/test/LazyOperations/UnionSet.jl
index ba7c9bdd3e..744dead894 100644
--- a/test/LazyOperations/UnionSet.jl
+++ b/test/LazyOperations/UnionSet.jl
@@ -115,7 +115,7 @@ for N in [Float64, Rational{Int}, Float32]
     @test isequivalent(Y.X, linear_map(A, UXY.X))
     @test isequivalent(Y.Y, linear_map(A, UXY.Y))
     Y = linear_map(A, Uarr)
-    for k in 1:length(array(Uarr))
+    for k in eachindex(array(Uarr))
         @test isequivalent(array(Y)[k], linear_map(A, array(Uarr)[k]))
     end