I want to use jump and guard conditions for my hybrid automata. But it is not working. Can anyone help?
using ReachabilityAnalysis, Symbolics
# Define the variables: height (x), velocity (v), and time (t)
@variables x v t
const vars = [x, v]
# Define the constants
const g = 9.81 # gravity (m/s^2)
const c = 0.8 # restitution coefficient (energy loss during bounce)
# Mode 1: Falling dynamics
function falling!(du, u, p, t)
x, v = u
du[1] = v # x' = v (position derivative)
du[2] = -g # v' = -g (velocity derivative due to gravity)
return du
end
# Mode 2: Bouncing dynamics (after hitting the ground)
function bouncing!(du, u, p, t)
x, v = u
du[1] = v # x' = v
du[2] = -g # v' = -g (same as falling, but with new initial velocity)
return du
end
# Define the hybrid automaton
function bouncing_ball()
n = 2 # number of dimensions (x, v)
automaton = GraphAutomaton(2) # 2 modes: falling and bouncing
# Mode 1: Falling
invariant_falling = HalfSpace(x >= 0, vars) # Ball cannot go below ground level
falling_system = @system(x' = falling!(x), dim:2, x ∈ invariant_falling)
# Mode 2: Bouncing
invariant_bouncing = HalfSpace(x >= 0, vars) # Also, ball cannot go below ground
bouncing_system = @system(x' = bouncing!(x), dim:2, x ∈ invariant_bouncing)
# Transition from falling to bouncing when ball hits the ground (x = 0)
guard = HPolyhedron([x <=0,v<=0], vars) # Trigger mode switch when x == 0
transition_map = Dict(v=>-c*v,x=>x) # Bounce with reduced velocity
add_transition!(automaton, 1, 2, 1) # Mode 1 (falling) -> Mode 2 (bouncing)
# Transition from bouncing back to falling after the bounce
add_transition!(automaton, 2, 1, 2) # Mode 2 (bouncing) -> Mode 1 (falling)
# Create the hybrid system
H = HybridSystem(; automaton=automaton, modes=[falling_system, bouncing_system], resetmaps=[transition_map])
# Initial conditions: ball starts from height 10 meters, with 0 initial velocity
X0 = Hyperrectangle([1.0, 0.0], [1, 0.0])
init = [(1, X0)] # Start in mode 1 (falling) at initial condition
println(1)
return InitialValueProblem(H, init)
end
# Solve the hybrid system
prob = bouncing_ball()
sol = solve(prob, tspan=(0.0, 2.0),
alg=TMJets21a(; abstol=1e-5, maxsteps=10_000, orderT=5, orderQ=1))
# Overapproximate the solution using zonotopes
solz = overapproximate(sol, Zonotope)