diff --git a/src/solve.jl b/src/solve.jl
index 909475ad..701cd6bf 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -1,4 +1,5 @@
 export solve
+import LazySets.constrained_dimensions
 
 function default_operator(system::InitialValueProblem{S}) where
         {S<:Union{AbstractContinuousSystem, AbstractDiscreteSystem}}
@@ -117,6 +118,32 @@ function solve(system::InitialValueProblem{<:HybridSystem, <:LazySet{N}},
     return solve!(sys_new, copy(options), opC, opD)
 end
 
+
+"""
+    constrained_dimensions(HS::HybridSystem)::Dict{Int,Vector{Int}}
+
+Return all coordinates which appear in any guard or invariant constraint for each location.
+
+### Input
+
+- `HS`  -- hybrid system
+"""
+function constrained_dimensions(HS::HybridSystem)::Dict{Int,Vector{Int}}
+    result = Dict{Int,Vector{Int}}()
+    sizehint!(result, nstates(HS))
+    for mode in states(HS)
+        vars = Vector{Int}()
+        append!(vars, constrained_dimensions(stateset(HS, mode)))
+        for transition in out_transitions(HS, mode)
+            append!(vars, constrained_dimensions(stateset(HS, transition)))
+        end
+        result[mode] = unique(vars)
+    end
+
+    return result
+end
+
+
 function init_states_sys_from_init_set_sys(
         system::InitialValueProblem{<:HybridSystem, <:LazySet{N}}) where N<:Real
     HS = system.s
diff --git a/test/Reachability/solve_hybrid_bouncing_ball.jl b/test/Reachability/solve_hybrid_bouncing_ball.jl
index 8d494468..f782b198 100644
--- a/test/Reachability/solve_hybrid_bouncing_ball.jl
+++ b/test/Reachability/solve_hybrid_bouncing_ball.jl
@@ -4,6 +4,7 @@
 
 using Reachability, HybridSystems, MathematicalSystems, LazySets, Polyhedra
 import LazySets.HalfSpace
+import Reachability.constrained_dimensions
 
 # Transition graph (automaton)
 a = LightAutomaton(1);
@@ -38,7 +39,11 @@ system = InitialValueProblem(HS, inits);
 options = Options(:mode=>"reach", :vars=>[1, 2], :T=>5.0, :δ=>0.1,
                   :plot_vars=>[1, 2], :max_jumps=>1, :verbosity=>1);
 
-# # default algorithm
+
+# test constrained_dimensions
+@test constrained_dimensions(HS) == Dict{Int,Vector{Int}}(1=>[1, 2])
+
+# default algorithm
 sol = solve(system, options);
 
 # specify lazy discrete post operator
diff --git a/test/Reachability/solve_hybrid_thermostat.jl b/test/Reachability/solve_hybrid_thermostat.jl
index 543f0a09..f05a7b79 100644
--- a/test/Reachability/solve_hybrid_thermostat.jl
+++ b/test/Reachability/solve_hybrid_thermostat.jl
@@ -8,6 +8,7 @@ using Reachability, HybridSystems, MathematicalSystems, LazySets, Polyhedra
 # fix namamespace conflicts with Polyhedra
 import LazySets.HalfSpace
 import Reachability.ReachSet
+import Reachability.constrained_dimensions
 
 c_a = 0.1;
 
@@ -58,6 +59,9 @@ system = InitialValueProblem(HS, X0);
 options = Options(:mode=>"reach", :vars=>[1], :T=>5.0, :δ=>0.1,
                   :max_jumps=>1, :verbosity=>1);
 
+# test constrained_dimensions
+@test constrained_dimensions(HS) == Dict{Int,Vector{Int}}(1=>[1], 2=>[1])
+
 # default algorithm
 sol = solve(system, options);