Update models

models/Airplane/Airplane.jl

@@ -176,12 +176,12 @@ U₀ = ZeroSet(6);
 ivp = @ivp(x' = Airplane!(x), dim: 18, x(0) ∈ X₀ × U₀)
 prob = ControlledPlant(ivp, controller, vars_idx, period);
 
-# The safety specification is that ``(x_2, x_7, x_8, x_9) ∈ ±[0.5, 1, 1, 1]``
+# The safety specification is that ``(x_2, x_7, x_8, x_9) ∈ ±[1, 1, 1, 1]``
 # for 20 control periods. A sufficient condition for guaranteed violation is to
 # overapproximate the result with hyperrectangles.
 
 safe_states = concretize(CartesianProductArray([
-    Universe(1), Interval(-0.5, 0.5), Universe(4),
+    Universe(1), Interval(-1.0, 1.0), Universe(4),
     BallInf(zeros(3), 1.0), Universe(9)]))
 
 predicate_set(R) = isdisjoint(overapproximate(R, Hyperrectangle), safe_states)
@@ -197,7 +197,7 @@ function predicate(sol; silent::Bool=false)
 end
 
 if falsification
-    k = 4  # falsification can run for a shorter time horizon
+    k = 7  # falsification can run for a shorter time horizon
 else
     k = 20
 end
@@ -274,8 +274,8 @@ vars = (2, 7)
 fig = plot_helper(vars)
 plot!(fig; xlab="s_y", ylab="ϕ", leg=:bottomleft)
 if falsification
-    xlims!(-0.01, 0.65)
-    ylims!(0.85, 1.01)
+    xlims!(-0.01, 1.15)
+    ylims!(0.5, 1.01)
 else
     xlims!(-0.55, 0.55)
     ylims!(-1.05, 1.05)

models/InvertedPendulum/InvertedPendulum.jl

@@ -1,22 +1,17 @@
 # # Inverted Pendulum
 #
-# The Inverted Pendulum benchmark is a classical model of motion.
+# The Inverted Pendulum benchmark is a classical model of motion. We consider
+# two scenarios, which we refer to as the verification and the falsification
+# scenario.
 
 module InvertedPendulum  #jl
 
 using ClosedLoopReachability
 import OrdinaryDiffEq, Plots, DisplayAs
 using ReachabilityBase.CurrentPath: @current_path
 using ReachabilityBase.Timing: print_timed
-using ClosedLoopReachability: SingleEntryVector
-using Plots: plot, plot!, xlims!, ylims!
-
-# The following option determines whether the falsification settings should be
-# used. The falsification settings are sufficient to show that the safety
-# property is violated. Concretely, we start from an initial point and use a
-# smaller time horizon.
-
-const falsification = true;
+using ClosedLoopReachability: SingleEntryVector, Specification, NoSplitter
+using Plots: plot, plot!, xlims!, ylims!, lens!, bbox, savefig
 
 # ## Model
 
@@ -67,136 +62,230 @@ period = 0.05;
 
 # ## Specification
 
-# The uncertain initial condition is ``θ \in [1, 1.2], \dot{θ} \in [0, 0.2]``.
+# The following script creates a different problem instance for the two
+# scenarios, respectively.
 
-X₀ = BallInf([1.1, 0.1], 0.1)
-if falsification
-    ## Choose a single point in the initial states (here: the top-most one):
-    X₀ = Singleton(high(X₀))
-end
-U₀ = ZeroSet(1);
+function InvertedPendulum_spec(verification::Bool)
+    ## The uncertain initial condition is ``\dot{θ} \in [0, 0.2]``, and ``θ``
+    ## depends on the scenario.
+    if verification
+        ## ``θ \in [1, 1.175]``.
+        X₀ = Hyperrectangle(low=[1.0, 0], high=[1.175, 0.2])
+    else
+        ## ``θ \in [1, 1.2]``. We choose a single point (here: the top-most one):
+        X₀ = Singleton(high(BallInf([1.1, 0.1], 0.1)))
+    end
+    U₀ = ZeroSet(1);
 
-# The control problem is:
+    ## The control problem is:
 
-ivp = @ivp(x' = InvertedPendulum!(x), dim: 3, x(0) ∈ X₀ × U₀)
-prob = ControlledPlant(ivp, controller, vars_idx, period);
+    ivp = @ivp(x' = InvertedPendulum!(x), dim: 3, x(0) ∈ X₀ × U₀)
+    prob = ControlledPlant(ivp, controller, vars_idx, period);
 
-# The safety specification is that ``θ ∈ [0, 1]`` for ``t ∈ [0.5, 1]`` (i.e.,
-# the control periods ``10 ≤ k ≤ 20``). A sufficient condition for guaranteed
-# violation is to overapproximate the result with hyperrectangles.
+    ## The safety specification is that ``θ ∈ [0, 1]`` for ``t ∈ [0.5, 1]``
+    ## (i.e., the control periods ``10 ≤ k ≤ 20``). A sufficient condition for a
+    ## guaranteed verdict is to overapproximate the result with hyperrectangles.
 
-unsafe_states = HalfSpace(SingleEntryVector(1, 3, -1.0), -1.0)
+    if verification
+        unsafe_states = UnionSet(HalfSpace(SingleEntryVector(1, 3, -1.0), -1.0),
+                                 HalfSpace(SingleEntryVector(1, 3, 1.0), 0.0))
+    else
+        unsafe_states = HalfSpace(SingleEntryVector(1, 3, -1.0), -1.0)
+    end
 
-function predicate_set(R)
-    t = tspan(R)
-    return t.lo >= 0.5 && t.hi <= 1.0 &&
-           overapproximate(R, Hyperrectangle) ⊆ unsafe_states
-end
+    function predicate_set_safe(R)
+        t = tspan(R)
+        return t.hi <= 0.5 ||
+            isdisjoint(overapproximate(R, Hyperrectangle), unsafe_states)
+    end
 
-function predicate(sol; silent::Bool=false)
-    for F in sol
-        t = tspan(F)
-        if t.hi < 0.5 || t.lo > 1.0
-            continue
+    function predicate_safe(sol; silent::Bool=false)
+        for F in sol
+            t = tspan(F)
+            if t.hi <= 0.5
+                continue
+            end
+            for R in F
+                if !predicate_set_safe(R)
+                    silent || println("  Potential violation for time range $(tspan(R)).")
+                    return false
+                end
+            end
         end
-        for R in F
-            if predicate_set(R)
-                silent || println("  Violation for time range $(tspan(R)).")
-                return true
+        return true
+    end
+
+    function predicate_set_unsafe(R)
+        t = tspan(R)
+        return t.lo >= 0.5 && t.hi <= 1.0 &&
+            overapproximate(R, Hyperrectangle) ⊆ unsafe_states
+    end
+
+    function predicate_unsafe(sol; silent::Bool=false)
+        for F in sol
+            t = tspan(F)
+            if t.hi < 0.5
+                continue
+            end
+            for R in F
+                if predicate_set_unsafe(R)
+                    silent || println("  Violation for time range $(tspan(R)).")
+                    return true
+                end
             end
         end
+        return false
     end
-    return false
-end
 
-if falsification
-    k = 11  # falsification can run for a shorter time horizon
-else
-    k = 20
+    if verification
+        predicate = predicate_safe
+    else
+        predicate = predicate_unsafe
+    end
+
+    if verification
+        T = 1.0
+    else
+        T = 11 * period  # falsification can run for a shorter time horizon
+    end
+
+    spec = Specification(T, predicate, unsafe_states)
+
+    return prob, spec
 end
-T = k * period
+
 T_warmup = 2 * period;  # shorter time horizon for warm-up run
 
 # ## Analysis
 
-# To enclose the continuous dynamics, we use a Taylor-model-based algorithm:
-
-algorithm_plant = TMJets(abstol=1e-7, orderT=4, orderQ=1);
-
-# To propagate sets through the neural network, we use the `DeepZ` algorithm:
+# To enclose the continuous dynamics, we use a Taylor-model-based algorithm. We
+# also use an additional splitting strategy to increase the precision. These
+# algorithms are defined later for each scenario. To propagate sets through the
+# neural network, we use the `DeepZ` algorithm:
 
 algorithm_controller = DeepZ();
 
 # The falsification benchmark is given below:
 
-function benchmark(; T=T, silent::Bool=false)
+function benchmark(prob, spec; T, algorithm_plant, splitter, verification, silent::Bool=false)
     ## Solve the controlled system:
     silent || println("Flowpipe construction:")
     res = @timed solve(prob; T=T, algorithm_controller=algorithm_controller,
-                       algorithm_plant=algorithm_plant)
+                       algorithm_plant=algorithm_plant, splitter=splitter)
     sol = res.value
     silent || print_timed(res)
 
     ## Check the property:
     silent || println("Property checking:")
-    res = @timed predicate(sol; silent=silent)
-    silent || print_timed(res)
-    if res.value
-        silent || println("  The property is violated.")
-        result = "falsified"
+    if verification
+        res = @timed spec.predicate(sol; silent=silent)
+        silent || print_timed(res)
+        if res.value
+            silent || println("  The property is verified.")
+            result = "verified"
+        else
+            silent || println("  The property may be violated.")
+            result = "not verified"
+        end
     else
-        silent || println("  The property may be satisfied.")
-        result = "not falsified"
+        res = @timed spec.predicate(sol; silent=silent)
+        silent || print_timed(res)
+        if res.value
+            silent || println("  The property is violated.")
+            result = "falsified"
+        else
+            silent || println("  The property may be satisfied.")
+            result = "not falsified"
+        end
     end
 
     return sol, result
+end
+
+function run(; verification::Bool)
+    if verification
+        println("# Running analysis with verification scenario")
+        algorithm_plant = TMJets(abstol=1e-9, orderT=5, orderQ=1)
+        splitter = BoxSplitter([[1.1, 1.16], [0.09, 0.145, 0.18]])
+    else
+        println("# Running analysis with falsification scenario")
+        algorithm_plant = TMJets(abstol=1e-7, orderT=4, orderQ=1)
+        splitter = NoSplitter()
+    end
+    prob, spec = InvertedPendulum_spec(verification)
+
+    ## Run the verification/falsification benchmark:
+    benchmark(prob, spec; T=T_warmup, algorithm_plant=algorithm_plant, splitter=splitter,
+              verification=verification, silent=true)  # warm-up
+    res = @timed benchmark(prob, spec; T=spec.T, algorithm_plant=algorithm_plant,  # benchmark
+                           splitter=splitter, verification=verification)
+    sol, result = res.value
+    if verification
+        @assert (result == "verified") "verification failed"
+    else
+        @assert (result == "falsified") "falsification failed"
+    end
+    println("Total analysis time:")
+    print_timed(res)
+
+    ## Compute some simulations:
+    println("Simulation:")
+    trajectories = verification ? 10 : 1
+    res = @timed simulate(prob; T=spec.T, trajectories=trajectories,
+                          include_vertices=verification)
+    sim = res.value
+    print_timed(res)
+
+    return sol, sim, prob, spec
 end;
 
-# Run the falsification benchmark and compute some simulations:
+# Run the analysis script for the verification scenario:
 
-benchmark(T=T_warmup, silent=true)  # warm-up
-res = @timed benchmark(T=T)  # benchmark
-sol, result = res.value
-@assert (result == "falsified") "falsification failed"
-println("Total analysis time:")
-print_timed(res)
+sol_v, sim_v, prob_v, spec_v = run(verification=true);
 
-println("Simulation:")
-res = @timed simulate(prob; T=T, trajectories=falsification ? 1 : 10,
-                      include_vertices=!falsification)
-sim = res.value
-print_timed(res);
+# Run the analysis script for the falsification scenario:
+
+sol_f, sim_f, prob_f, spec_f = run(verification=false);
 
 # ## Results
 
 # Script to plot the results:
 
-function plot_helper()
+function plot_helper(sol, sim, prob, spec, verification)
     vars = (0, 1)
-    fig = plot()
-    unsafe_states_projected = cartesian_product(Interval(0.5, 1.0),
-                                                project(unsafe_states, [vars[2]]))
-    plot!(fig, unsafe_states_projected; color=:red, alpha=:0.2, lab="unsafe")
-    plot!(fig, sol; vars=vars, color=:yellow, lw=0, alpha=1, lab="")
-    initial_states_projected = cartesian_product(Singleton([0.0]), project(X₀, [vars[2]]))
-    plot!(fig, initial_states_projected; c=:cornflowerblue, alpha=1, lab="X₀")
-    if falsification
-        xlims!(0, T)
-        ylims!(0.95, 1.22)
-    else
-        xlims!(0, T)
-        ylims!(0.55, 1.3)
+    fig = plot(leg=:topright)
+    lab = "unsafe"
+    unsafe_states = spec.ext isa UnionSet ? spec.ext : [spec.ext]
+    for B in unsafe_states
+        unsafe_states_projected = cartesian_product(Interval(0.5, 1.0), project(B, [vars[2]]))
+        plot!(fig, unsafe_states_projected; color=:red, alpha=:0.2, lab=lab)
+        lab = ""
     end
-    lab_sim = falsification ? "simulation" : ""
+    plot!(fig, sol; vars=vars, color=:yellow, lw=0, alpha=1, lab="")
+    initial_states_projected =
+        cartesian_product(Singleton([0.0]), project(initial_state(prob), [vars[2]]))
+    plot!(fig, initial_states_projected; c=:cornflowerblue, alpha=1, m=:none, lw=7, lab="X₀")
+    lab_sim = verification ? "" : "simulation"
     plot_simulation!(fig, sim; vars=vars, color=:black, lab=lab_sim)
+    xlims!(0, spec.T)
     plot!(fig; xlab="t", ylab="θ")
     return fig
 end;
 
 # Plot the results:
 
-fig = plot_helper()
-## Plots.savefig(fig, "InvertedPendulum.png")  # command to save the plot to a file
+fig = plot_helper(sol_v, sim_v, prob_v, spec_v, true)
+ylims!(fig, 0.5, 1.2)
+lens!(fig, [0.49, 0.52], [0.99, 1.01]; inset=(1, bbox(0.1, 0.6, 0.3, 0.3)),
+      lc=:black, xticks=[0.5], yticks=[1.0], subplot=3)
+## Plots.savefig(fig, "InvertedPendulum_verification.png")  # command to save the plot to a file
+fig = DisplayAs.Text(DisplayAs.PNG(fig))
+
+#-
+
+fig = plot_helper(sol_f, sim_f, prob_f, spec_f, false)
+ylims!(fig, 0.95, 1.22)
+## Plots.savefig(fig, "InvertedPendulum_falsification.png")  # command to save the plot to a file
 fig = DisplayAs.Text(DisplayAs.PNG(fig))
 
 end  #jl