update docs

docs/make.jl

@@ -38,8 +38,7 @@ makedocs(
                           "Platoon" => "models/platoon.md",
                           "Quadrotor" => "models/quadrotor.md",
                           "SEIR model" => "models/seir.md",
-                          "Spacecraft" => "models/spacecraft.md",
-                          "Platoon" => "models/platoon.md"],
+                          "Spacecraft" => "models/spacecraft.md"],
                           # "Electromechanic break" => "man/applications/embrake.md"
         "Algorithms" => Any["ASB07" => "lib/algorithms/ASB07.md",
                             "BFFPSV18" => "lib/algorithms/BFFPSV18.md",

examples/ISS.jl

@@ -2,18 +2,29 @@
 #md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/ISS.ipynb)
 #
 #md # !!! note "Overview"
-#md #     System type: linear continuous system\
+#md #     System type: linear time-invariant continuous system\
 #md #     State dimension: 270\
-#md #     Application domain: aerospace engineering
+#md #     Application domain: Aerospace Engineering
 #
 # ## Model description
 #
 # The International Space Station (ISS) is a continuous linear time-invariant
-# system ``\dot{x}(t) = Ax(t) + Bu(t)`` proposed as a benchmark in ARCH
-# 2016 [TLT16]. It has 270 state variables.
+# system
+#
+# ```math
+#   \begin{array}{lcl}
+#   \dot{x}(t) &=& Ax(t) + Bu(t) \\
+#   y(t) & = & C x(t)
+#   \end{array}
+# ```
+# proposed as a benchmark in ARCH 2016 [TLT16]. It has 270 state variables,
+# ``A ∈ \mathbb{R}^{270\times 270}``, ``B ∈ \mathbb{R}^{270\times 3}``,
+# ``C ∈ \mathbb{R}^{3\times 270}``.
 
 # The A, B, and C matrices are available in MATLAB format
-# ![](slicot.org/objects/software/shared/bench-data/iss.zip)
+# ![](slicot.org/objects/software/shared/bench-data/iss.zip) (`iss.mat`). For convenience
+# such `.mat` has been converted to the [JLD2 format](https://github.com/JuliaIO/JLD2.jl)
+# and stored in the file `iss.jld2`.
 
 using ReachabilityAnalysis, JLD2
 using ReachabilityAnalysis: add_dimension
@@ -27,58 +38,103 @@ const C3_ext = vcat(C3, fill(0.0, 3));
 
 # ## Reachability settings
 #
-# Initially all the variables are in the range ``[-0.0001, 0.0001]``, ``u_1``
-# is in ``[0, 0.1]``, ``u_2`` is in ``[0.8, 1]``, and ``u_3`` is in ``[0.9, 1]``.
+# Initially all the variables are in the range ``[-0.0001, 0.0001]``, and the inputs
+# are bounded: ``u_1(t) ∈ [0, 0.1]``, ``u_2(t) ∈ [0.8, 1]``, and ``u_3`` is in ``[0.9, 1]``.
 # The time bound is 20.
 
-# There are two versions of this benchmark:
+# There are two versions of this benchmark, ISSF01 (time-varying inputs) and
+# ISSC01 (constant inputs).
 
-# - ISSF01: The inputs can change arbitrarily over time: ``\forall t: u(t)\in \mathcal{U}``.
+# - **ISSF01** (time-varying inputs): In this setting, the inputs can change
+#   arbitrarily over time: ``\forall t: u(t) \in \mathcal{U}``.
 
 function ISSF01()
     @load ISS_path A B
 
-    U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.]);  # input set
-    X0 = BallInf(zeros(size(A, 1)), 0.0001)  # -0.0001 <= xi <= 0.0001 for all i
-    prob_ISSF01 = @ivp(x' = A*x + B*u, x(0) ∈ X0, u ∈ U, x ∈ Universe(270))
+    U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.])
+    X0 = BallInf(zeros(size(A, 1)), 0.0001)
+    return @ivp(x' = A*x + B*u, x(0) ∈ X0, u ∈ U, x ∈ Universe(270))
 end
 
-# - ISSC01 (constant inputs): The inputs are uncertain only in their initial
-# value, and constant over time: ``u(0)\in \mathcal{U}``, ``\dot u(t)= 0``.
+# - **ISSC01** (constant inputs): The inputs are uncertain only in their initial
+#   value, and constant over time: ``u(0)\in \mathcal{U}``, ``\dot u(t)= 0``.
 
 function ISSC01()
     @load ISS_path A B
 
-    Aext = add_dimension(A, 3)
-    Aext[1:270, 271:273] = B
-    S = LinearContinuousSystem(Aext)
+    A_ext = add_dimension(A, 3)
+    A_ext[1:270, 271:273] = B
 
-    U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.])  # input set
-    X0 = BallInf(zeros(size(A, 1)), 0.0001)  # -0.0001 <= xi <= 0.0001 for all i
+    U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.])
+    X0 = BallInf(zeros(size(A, 1)), 0.0001)
     X0 = X0 * U
-    prob_ISSC01 = InitialValueProblem(S, X0)
+    return @ivp(x' = A_ext*x, x(0) ∈ X0)
 end
 
+# ## Specifications
+
+# The verification goal is to check the ranges reachable by the output ``y_3(t)``,
+# which is a linear combination of the state variables. In addition to the safety specification,
+# for each version there is an UNSAT instance that serves as sanity checks to ensure
+# that the model and the tool work as intended. But there is a caveat: In principle,
+# verifying an UNSAT instance only makes sense formally if a witness is provided
+# (counter-example, under-approximation, etc.).
+
+# Since `ReachabilityAnalysis.jl` currently doesn't have the capability of verifying
+# UNSAT instances, we run the tool with the same accuracy settings on an SAT-UNSAT pair of instances.
+# The SAT instance demonstrates that the over-approximation is not too coarse,
+# and the UNSAT instance demonstrates that the over-approximation is indeed conservative,
+# at least in the narrow sense of the specification.
+
+# - ISS01: Bounded time, safe property: For all ``t ∈ [0, 20]``, ``y_3(t) ∈ [−0.0007, 0.0007]``.
+#          This property is used with the uncertain input case (ISSF01) and assumed to be satisfied.
+#
+# - ISS02: Bounded time, safe property: For all ``t ∈ [0, 20]``, ``y_3(t) ∈ [−0.0005, 0.0005]``.
+#          This property is used with the constant input case (ISSC01) and assumed to be satisfied.
+#
+# - ISU01: Bounded time, unsafe property: For all ``t ∈ [0, 20]``, ``y_3 (t) ∈ [−0.0005, 0.0005]``.
+#          This property is used with the uncertain input case (ISSF01) and assumed to be unsatisfied.
+#
+# - ISU02: Bounded time, unsafe property: For all ``t ∈ [0, 20]``, ``y_3 (t) ∈ [−0.00017, 0.00017]``.
+#          This property is used with the constant input case (ISSC01) and assumed to be unsatisfied.
 
 # ## Results
 
-# We use `LGG09` algorithm (based on support functions [LGG09])
+# The specification involves only the output ``y(t) := C_3 x(t)``, where ``C_3`` denotes
+# the third row of the output matrix ``C ∈ \mathbb{R}^{3\times 270}``. Hence,
+# it is sufficient to compute the flowpipe associated to ``y(t)`` directly,
+# without the need of actually computing the full 270-dimensional flowpipe associated to
+# all state variables. The flowpipe associated to a linear combination of state
+# variables can be computed efficiently using the support-function based algorithm
+# [[LGG09]](@ref)). The idea is to define a template polyhedron with only two
+# supporting directions, namely ``C_3`` and ``-C_3``.
 
-# - ISSF01
+# The step sizes chosen are ``6×10^{-4}`` for ISSF01 and ``1×10^{-2}`` for ISSC01.
+
+# ### ISSF01
 
 dirs = CustomDirections([C3, -C3]);
 prob_ISSF01 = ISSF01();
 sol_ISSF01 = solve(prob_ISSF01, T=20.0, alg=LGG09(δ=6e-4, template=dirs, sparse=true, cache=false));
 
-#-
+# We can transform the flowpipe on the output ``y_3(t)`` by taking the projection
+# with the vector `C3`.
 
-using Plots, Plots.PlotMeasures, LaTeXStrings
+πsol_ISSF01 = project(sol_ISSF01, C3);
 
-out = project(sol_ISSF01[1:10:length(sol_ISSF01)], C3, vars=(0,1));
-# - out = [Interval(tspan(sol_ISSF01[i])) × Interval(-ρ(-C3, sol_ISSF01[i]), ρ(C3, sol_ISSF01[i])) for i in 1:10:length(sol_ISSF01)];
+#md # !!! note "Technical note"
+#         Note that projecting the solution along direction ``C_3`` corresponds to computing
+#         the min and max bounds for each reach-set `X`, `Interval(-ρ(-C3, X), ρ(C3, X)`.
+#         However, the method `project(sol_ISSF01, C3)` is more efficient than manually
+#         evaluating the support functions because it doesn't recompute them alongf
+#         directions ``C_3`` and ``-C_3``, because they are already known to the flowpipe
+#         structure, or more precisely, `project` checks and finds that the given
+#         directions belong to those represented by each `TemplatReachSet`.
+
+using Plots, Plots.PlotMeasures, LaTeXStrings
 
 fig = Plots.plot();
-Plots.plot!(fig, out, linecolor=:blue, color=:blue, alpha=0.8,
+Plots.plot!(fig, πsol_ISSF01[1:10:end], vars=(0, 1), linecolor=:blue, color=:blue, alpha=0.8,
 #!nb     #tickfont=font(30, "Times"), guidefontsize=45, #!jl
      xlab=L"t",
      ylab=L"y_{3}",
@@ -88,19 +144,21 @@ Plots.plot!(fig, out, linecolor=:blue, color=:blue, alpha=0.8,
 #!nb     #size=(1000, 1000)); #!jl
 fig
 
-# - ISSC01
+# ### ISSC01
 
 dirs = CustomDirections([C3_ext, -C3_ext]);
 prob_ISSC01 = ISSC01();
 sol_ISSC01 = solve(prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=dirs, sparse=true, cache=false));
 
 #-
 
-out = project(sol_ISSC01, C3_ext, vars=(0,1));
-# - out = [Interval(tspan(sol_ISSC01[i])) × Interval(-ρ(-C3_ext, sol_ISSC01[i]), ρ(C3_ext, sol_ISSC01[i])) for i in 1:1:length(sol_ISSC01)];
+# We can transform the flowpipe on the output ``y_3(t)`` by taking the projection
+# with the vector `C3`.
+
+πsol_ISSC01 = project(sol_ISSC01, C3_ext);
 
 fig = Plots.plot();
-Plots.plot!(fig, out, linecolor=:blue, color=:blue, alpha=0.8, lw=1.0,
+Plots.plot!(fig, πsol_ISSC01[1:10:end], vars=(0, 1), linecolor=:blue, color=:blue, alpha=0.8, lw=1.0,
 #!nb    #tickfont=font(30, "Times"), guidefontsize=45, !jl
     xlab=L"t",
     ylab=L"y_{3}",
@@ -112,5 +170,6 @@ fig
 
 # ## References
 
-# - [TLT16] Tran, Hoang-Dung, Luan Viet Nguyen, and Taylor T. Johnson. "Large-scale linear systems from order-reduction (benchmark proposal)." 3rd Applied Verification for Continuous and Hybrid Systems Workshop (ARCH), Vienna, Austria. 2016.
-# - [LGG09] Le Guernic, Colas, and Antoine Girard. "Reachability analysis of linear systems using support functions." Nonlinear Analysis: Hybrid Systems 4.2 (2010): 250-262.
+# - [TLT16] Tran, Hoang-Dung, Luan Viet Nguyen, and Taylor T. Johnson.
+#   *Large-scale linear systems from order-reduction (benchmark proposal).*
+#   3rd Applied Verification for Continuous and Hybrid Systems Workshop (ARCH), Vienna, Austria. 2016.

examples/building.jl

@@ -8,16 +8,21 @@
 #
 # ## Model description
 #
-# This benchmark is quite straightforward: The system is described by
-# ``\dot{x}(t) = Ax(t) + Bu(t)``, ``u(t) \in \mathcal{U}``, ``y(t) = Cx(t)``
+# The system is described by the linear differential equations:
+#
+# ```math
+#   \begin{array}{lcl}
+#   \dot{x}(t) &=& Ax(t) + Bu(t),\qquad u(t) \in \mathcal{U} \\
+#   y(t) &=& C x(t)
+#   \end{array}
+# ```
 #
-# Discrete-time analysis for the building system should use a step size of ``0.01``.
 # There are two versions of this benchmark:
-# - The inputs can change arbitrarily over time: $\forall t: u(t)\in \mathcal{U}$.
-# - (constant inputs) The inputs are uncertain only in their initial value, and
-#    constant over time: ``u(0)\in \mathcal{U}``, ``\dot u (t)= 0``. The purpose
-#    of this model instance is to accommodate tools that cannot handle
-#    time-varying inputs.
+#
+# - *(time-varying inputs):* The inputs can change arbitrarily over time: $\forall t: u(t)\in \mathcal{U}$.
+#
+# - *(constant inputs):* The inputs are uncertain only in their initial value, and
+#    constant over time: ``u(0)\in \mathcal{U}``, ``\dot u (t)= 0``.
 
 using ReachabilityAnalysis, SparseArrays, JLD2
 
@@ -94,9 +99,10 @@ end
 #    account. A tool should be run with the same accuracy settings on BLDF01-BDU02
 #    and BLDC01-BDU02, returning UNSAT on the former and SAT on the latter.
 
-
 # ## Results
 
+# For the discrete-time analysis we use a step size of ``0.01``.
+
 using Plots
 
 # ### BLDF01

examples/lorenz.jl

@@ -1,5 +1,5 @@
-# # SEIR
-#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/seir.ipynb)
+# # Lorenz equations
+#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/lorenz.ipynb)
 #
 #md # !!! note "Overview"
 #md #     System type: Continuous blackbox system\
@@ -13,9 +13,9 @@
 #
 # ```math
 #   \begin{array}{lcl}
-#   \frac{dx}{dt} & = & \sigma (y-x), \\
-#   \frac{dy}{dt} & = & x(\rho - z) - y, \\
-#   \frac{dz}{dt} & = & xy-\beta z
+#   \dfrac{dx}{dt} & = & \sigma (y-x), \\ \\
+#   \dfrac{dy}{dt} & = & x(\rho - z) - y, \\ \\
+#   \dfrac{dz}{dt} & = & xy-\beta z
 #   \end{array}
 # ```
 #
@@ -28,8 +28,6 @@
 # the Prandtl number, Rayleigh number, and certain physical dimensions of the
 # layer itself.
 
-
-
 using ReachabilityAnalysis, Plots
 
 @taylorize function lorenz!(dx, x, params, t)
@@ -42,24 +40,24 @@ using ReachabilityAnalysis, Plots
     return dx
 end
 
-X0 = Hyperrectangle(low=[0.9, 0.0, 0.0], high=[1.1, 0.0, 0.0])
-prob = @ivp(x' = lorenz!(x), dim=3, x(0) ∈ X0);
-
 # ## Reachability settings
 #
-# The initial values used were ``X₀ \in [0.9, 1.1] \times [0, 0] \times [0, 0]``,
-# for a time span of 200.
-# The algorithm used was `TMJets` with ``n_T=10`` and ``n_Q=2``
+# The initial values considered are ``X_0 \in [0.9, 1.1] \times [0, 0] \times [0, 0]``,
+# for a time span of `10`.
 
+X0 = Hyperrectangle(low=[0.9, 0.0, 0.0], high=[1.1, 0.0, 0.0])
+prob = @ivp(x' = lorenz!(x), dim=3, x(0) ∈ X0);
 
 # ## Results
 
-alg = TMJets(abs_tol=1e-15, orderT=10, orderQ=2,
-                         max_steps=50_000);
-# ### Compute flowpipe
+# We compute the flowpipe using the TMJets algorithm with ``n_T=10`` and ``n_Q=2``.
+
+alg = TMJets(abs_tol=1e-15, orderT=10, orderQ=2, max_steps=50_000);
+
 sol = solve(prob, T=10.0, alg=alg)
-solz = overapproximate(sol, Zonotope);
 
-# ### Plot the result
+# Below we plot the flowpipe projected on the `(1, 3)` plane.
+
+solz = overapproximate(sol, Zonotope);
 
 plot(solz, vars=(1, 3))

src/Flowpipes/flowpipes.jl

@@ -341,7 +341,7 @@ function project(fp::Flowpipe, vars::NTuple{D, T}) where {D, T<:Integer}
     end
 end
 
-project(fp::Flowpipe, vars::AbstractVector) = project(fp, Tuple(vars))
+project(fp::Flowpipe, vars::AbstractVector{<:Int}) = project(fp, Tuple(vars))
 project(fp::Flowpipe; vars) = project(fp, Tuple(vars))
 project(fp::Flowpipe, i::Int, vars) = project(fp[i], vars)
 
@@ -357,9 +357,9 @@ function project(fp::Flowpipe, M::AbstractMatrix; vars=nothing)
 end
 
 # concrete projection of a flowpipe for a given direction
-function project(fp::Flowpipe, M::AbstractVector; vars=nothing)
+function project(fp::Flowpipe, dir::AbstractVector{<:AbstractFloat}; vars=nothing)
     Xk = array(fp)
-    πfp = Flowpipe(map(X -> project(X, M, vars=vars), Xk))
+    return Flowpipe(map(X -> project(X, dir, vars=vars), Xk))
 end
 
 # concrete linear map of a flowpipe for a given matrix