Simulation for hybrid systems

[schillic]

This PR is based on #539.

• better cut-off for invariants (either create more small steps in the ODE solver or use an event-based detection)
• nondeterministic transitions
• return useful struct (currently it is a vector of vectors of solutions)
• test
  • initial set
  • multiple initial locations
  • small invariant-guard intersections
  • guards that get activated-deactivated-activated

Possible future additions:

• urgent transitions
• interpolation between time steps (currently we need to enforce a small time step for reasonable results)

Example

using ReachabilityAnalysis
using HybridSystems, MathematicalSystems, LazySets, LinearAlgebra
using LazySets: HalfSpace  # resolve name-space conflicts with Polyhedra

function thermostat()
    c_a = 0.1

    # automaton structure
    automaton = LightAutomaton(2)

    # mode on
    A = hcat(c_a)
    inv = HalfSpace([1.0], 22.0)  # x ≤ 22
    # @vars x v
    # @set x ≤ 22
    m_on = ConstrainedLinearContinuousSystem(A, inv)

    # mode off
    A = hcat(-c_a)
    inv = HalfSpace([-1.0], -18.0)  # x ≥ 18
    m_off = ConstrainedLinearContinuousSystem(A, inv)

    # modes
    modes = [m_on, m_off]

    # transition from on to off
    add_transition!(automaton, 1, 2, 1)
    guard = HalfSpace([-1.0], -21.0)  # x ≥ 21
    t_on2off = ConstrainedIdentityMap(2, guard)

    # transition from off to on
    add_transition!(automaton, 2, 1, 2)
    guard = HalfSpace([1.0], 19.0)  # x ≤ 19
    t_off2on = ConstrainedIdentityMap(2, guard)

    # transition annotations
    resetmaps = [t_on2off, t_off2on]

    # switching
    switchings = [AutonomousSwitching()]

    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)

    # initial condition in mode on
    X0 = Singleton([18.0])
    initial_condition = [(1, X0)]

    return InitialValueProblem(ℋ, initial_condition)
end;

S = thermostat();
T = 5.0;

# ---

sol = solve(S, T=T);

using Plots

plot(sol, vars=(0, 1))

# ---

import DifferentialEquations

# will have many intermediate steps
solsd = solve(S, T=T, ensemble=true, trajectories=10, use_discrete_callback=true);
# will handle invariants continuously but not have many intermediate steps
solsc = solve(S, T=T, ensemble=true, trajectories=10);
# will handle invariants continuously with additionally many intermediate steps
solscstep = solve(S, T=T, ensemble=true, trajectories=10, dtmax=0.1);

for piece in ensemble(solsd)
    plot!(piece, vars=(0, 1), c=:red, lw=2, alpha=1, lab="")
end
for piece in ensemble(solsc)
    plot!(piece, vars=(0, 1), c=:black, lw=2, alpha=1, lab="")
end
for piece in ensemble(solscstep)
    plot!(piece, vars=(0, 1), c=:blue, lw=2, alpha=1, lab="")
end
plot!()

[schillic]

The plot shows that the set-based results are incorrect: the last two jumps are missing!

[schillic]

Tests:

# transition from off to on
guard1 = HalfSpace([1.0], 19.0)  # original guard
guard2 = HalfSpace([1.0], 18.01)  # small invariant-guard intersection
# guard is a union (does not work with reachability analysis; I needed to manually deactivate it)
guard3 = UnionSet(Interval(18., 18.3), Interval(19., 19.3))

# initial condition
X0 = Singleton([18.0])
X1 = BallInf([18.0], 0.3)
X2 = BallInf([21.0], 0.3)
X3 = BallInf([21.0], 0.1)
initial_condition1 = [(1, X0)]  # original singleton condition
initial_condition2 = [(1, X1)]  # whole set of points
initial_condition3 = [(1, X1), (2, X2)]  # several initial locations
initial_condition4 = X3  # set intersecting with both location invariants

1. initial_condition1 and guard1
2. initial_condition2 and guard1
3. initial_condition3 and guard1
4. initial_condition4 and guard1
5. initial_condition1 and guard2
6. initial_condition1 and guard3

As can be seen, there are many transitions missing in the set-based analysis 😱

[mforets]

The confusion is that the fixpoint check is true by default.

[mforets]

i think its good; the docs build is unrelated