Merge pull request #783 from JuliaReach/schillic/examples

Revise models

docs/generate.jl

@@ -19,12 +19,12 @@ MODELS = [
           "Lorenz",
           "LotkaVolterra",
           "OpAmp",
-          "SquareWaveOscillator",
           "Platoon",
           "ProductionDestruction",
           "Quadrotor",
           "SEIR",
           "Spacecraft",
+          "SquareWaveOscillator",
           "TransmissionLine",
           "VanDerPol"
           #
@@ -54,8 +54,10 @@ for model in MODELS
             # Literate.notebook(input, target_dir_md; execute=true, credit=false)
             # if used, add the following to the top of the script files (where `MODELNAME` is the model name):
             #md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated_examples/MODELNAME.ipynb)
-        elseif any(endswith.(file, [".jld2", ".png"]))
-            # ignore *.jld2 and *.png files without warning
+        elseif any(endswith.(file, [".png", ".jpg", ".gif"]))
+            cp(joinpath(model_path, file), joinpath(target_dir_md, file); force=true)
+        elseif any(endswith.(file, [".jld2"]))
+            # ignore *.jld2 files without warning
         else
             @warn "ignoring $file in $model_path"
         end

docs/make.jl

@@ -44,7 +44,7 @@ LINEAR_METHODS = ["Linear methods" => [
                                        #"Propagating hyperrectangles"   => "tutorials/linear_methods/hyperrectangles_setprop.md",
                                        #"Propagating support functions" => "tutorials/linear_methods/supfunc_setprop.md",
                                        #"Helicopter model"              => "tutorials/linear_methods/helicopter.md",
-                                       "Structural model" => "generated_examples/ISS.md"
+                                       "International Space Station" => "generated_examples/ISS.md"
                                        #
                                        ]]
 
@@ -86,8 +86,8 @@ HYBRID_SYSTEMS = ["Hybrid systems" => [
 CLOCKED_SYSTEMS = ["Clocked systems" => [
                                          #
                                          #"Introduction"           => "tutorials/clocked_systems/introduction.md",
-                                         "Platoon"                => "generated_examples/Platoon.md",
-                                         "Electro-mechanic break" => "generated_examples/EMBrake.md"
+                                         "Platoon"                 => "generated_examples/Platoon.md",
+                                         "Electromechanical brake" => "generated_examples/EMBrake.md"
                                          #
                                          ]]
 
@@ -209,7 +209,7 @@ makedocs(;
                                        "Production-Destruction" => "generated_examples/ProductionDestruction.md",
                                        "Brusselator"            => "generated_examples/Brusselator.md",
                                        "Quadrotor"              => "generated_examples/Quadrotor.md",
-                                       "Epidemic model"         => "generated_examples/SEIR.md",
+                                       "Epidemic (SEIR) model"  => "generated_examples/SEIR.md",
                                        "Spacecraft"             => "generated_examples/Spacecraft.md"
                                        #
                                        ],

examples/Brusselator/Brusselator.jl

@@ -1,126 +1,136 @@
 # # Brusselator
-#
+
 #md # !!! note "Overview"
-#md #     System type: polynomial continuous system\
+#md #     System type: Polynomial continuous system\
 #md #     State dimension: 2\
 #md #     Application domain: Chemical kinetics
-#
+
 # ## Model description
-#
-# A chemical reaction is said to be *autocatalytic* if one of the reaction products is
-# also a catalyst for the same or a coupled reaction, and such a reaction is called an autocatalytic reaction.
-# We refer to the wikipedia article [Autocatalysis](https://en.wikipedia.org/wiki/Autocatalysis) for details.
 
+# A chemical reaction is said to be *autocatalytic* if one of the reaction
+# products is also a catalyst for the same or a coupled reaction. We refer to
+# [Wikipedia](https://en.wikipedia.org/wiki/Autocatalysis) for details.
+#
 # The Brusselator is a mathematical model for a class of autocatalytic reactions.
-# The dynamics of the Brusselator is given by the two-dimensional ODE
+# The dynamics of the Brusselator is given by the following two-dimensional
+# differential equations.
 #
 # ```math
-#   \left\{ \begin{array}{lcl} \dot{x} & = & A + x^2\cdot y - B\cdot x - x \\
-#    \dot{y} & = & B\cdot x - x^2\cdot y \end{array} \right.
+# \begin{aligned}
+#     \dot{x} &= A + x^2 ⋅ y - B ⋅ x - x \\
+#     \dot{y} &= B ⋅ x - x^2 ⋅ y
+# \end{aligned}
 # ```
-
-# The numerical values for the model's constants (in their respective units) are
-# given in the following table.
 #
-# |Quantity|Value|
-# |-----|-------|
-# |A    | 1     |
-# |B    | 1.5   |
-
-using ReachabilityAnalysis #!jl
+# The numerical values for the model's constants (in their respective units) are
+# ``A = 1, B = 1.5``.
 
-const A = 1.0
-const B = 1.5
-const B1 = B + 1
+using ReachabilityAnalysis  #!jl
 
 @taylorize function brusselator!(du, u, p, t)
+    local A = 1.0
+    local B = 1.5
+    local B1 = B + 1
+
     x, y = u
-    x² = x * x
-    aux = x² * y
-    du[1] = A + aux - B1 * x
-    du[2] = B * x - aux
+
+    x²y = x^2 * y
+    du[1] = A + x²y - B1 * x
+    du[2] = B * x - x²y
     return du
 end
 
 #md # !!! tip "Performance tip"
-#md #     The auxiliary variables `B1`, `x²` and `aux` have been defined to make
-#md #     better use of `@taylorize` and help to reduce allocations.
+#md #     The auxiliary variables `B1` and `x²y` improve the performance of
+#md #     `@taylorize`.
 
-# ## Reachability settings
-#
-# The initial set is defined by ``x \in [0.8, 1]``, ``y \in [0, 0.2]``.
-# These settings are taken from [^CAS13].
+# ## Specification
 
-U₀ = (0.8 .. 1.0) × (0.0 .. 0.2); #!jl
-prob = @ivp(u' = brusselator!(U), u(0) ∈ U₀, dim:2); #!jl
+# The initial set is ``x ∈ [0.8, 1]``, ``y ∈ [0, 0.2]``. The time horizon is
+# ``18``. These settings are taken from [^CAS13].
 
-# ## Results
+U₀ = (0.8 .. 1.0) × (0.0 .. 0.2)  #!jl
+prob = @ivp(u' = brusselator!(u), u(0) ∈ U₀, dim:2)  #!jl
 
-# We use `TMJets` algorithm with sixth-order expansion in time and second order expansion
-# in the spatial variables.
+T = 18.0;  #!jl
 
-sol = solve(prob; T=18.0, alg=TMJets20(; orderT=6, orderQ=2)); #!jl
+# ## Analysis
 
-#-
+# We use the `TMJets` algorithm with sixth-order expansion in time and
+# second-order expansion in the spatial variables. The version `TMJets20` is
+# more efficient here.
 
-using Plots #!jl
+alg = TMJets20(; orderT=6, orderQ=2)  #!jl
+sol = solve(prob; T=T, alg=alg);  #!jl
 
-fig = plot(sol; vars=(1, 2), xlab="x", ylab="y", lw=0.2, color=:blue, lab="Flowpipe", #!jl
-           legend=:bottomright) #!jl
-plot!(U₀; color=:orange, lab="Uo") #!jl
+# ## Results
 
+using Plots  #!jl
 #!jl import DisplayAs  #hide
+
+fig = plot(sol; vars=(1, 2), xlab="x", ylab="y", lw=0.2, color=:blue,  #!jl
+           lab="Flowpipe", legend=:bottomright)  #!jl
+plot!(fig, U₀; color=:orange, lab="U₀")  #!jl
+
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
-# We observe that the system converges to the equilibrium point `(1.0, 1.5)`.
+# We observe that the system converges to the equilibrium point ``(1.0, 1.5)``.
 
-# Below we plot the flowpipes projected into the time domain.
+# Below we plot the projected flowpipes over time.
 
 fig = plot(sol; vars=(0, 1), xlab="t", lw=0.2, color=:blue, lab="x(t)", legend=:bottomright)  #!jl
-plot!(sol; vars=(0, 2), xlab="t", lw=0.2, color=:red, lab="y(t)")  #!jl
+plot!(fig, sol; vars=(0, 2), lw=0.2, color=:red, lab="y(t)")  #!jl
 
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
 # ## Changing the initial volume
 
-# The model was considered in [^GCLASG20] but using a different set of initial conditions. Let us parametrize
-# the initial states as a ball centered at ``x = y = 1`` and radius ``r > 0``:
-
-U0(r) = Singleton([1.0, 1.0]) ⊕ BallInf(zeros(2), r)
+# The model was considered in [^GCLISG20], but using a different set of initial
+# states and a time horizon ``30``. Let us parametrize the initial states as a
+# square centered at ``x = y = 1`` and radius ``r > 0``:
 
-#-
+U0(r) = BallInf([1.0, 1.0], r);
 
 # The parametric initial-value problem is defined accordingly.
 
-bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim:2)
+bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim:2);
 
 # First we solve for ``r = 0.01``:
 
-sol_01 = solve(bruss(0.01); T=30.0, alg=TMJets20(; orderT=6, orderQ=2))  #!jl
+sol_01 = solve(bruss(0.01); T=30.0, alg=alg);  #!jl
+
+#-
 
-LazySets.set_ztol(Float64, 1e-15)  #!jl
+old_ztol = LazySets._ztol(Float64)  #!jl
+LazySets.set_ztol(Float64, 1e-15)  # use higher precision for the plots #!jl
 
 fig = plot(sol_01; vars=(1, 2), xlab="x", ylab="y", lw=0.2, color=:blue, lab="Flowpipe (r = 0.01)",  #!jl
            legend=:bottomright)  #!jl
-
-plot!(U0(0.01); color=:orange, lab="Uo", xlims=(0.6, 1.3))  #!jl
+plot!(fig, U0(0.01); color=:orange, lab="Uo", xlims=(0.6, 1.3))  #!jl
 
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
-# We observe that the wrapping effect is controlled and the flowpipe doesn't blow up even for the large time horizon ``T = 30.0``.
-# Next we plot the flowpipe zoomed to the last portion and compare ``r = 0.01`` with a set of larger initial states, ``r = 0.1``.
+# We observe that the wrapping effect is controlled and the flowpipe does not
+# blow up even for the large time horizon ``T = 30.0``. Next we plot the
+# flowpipe zoomed to the last portion and compare ``r = 0.01`` with a larger set
+# of initial states, ``r = 0.1``.
 
-sol_1 = solve(bruss(0.1); T=30.0, alg=TMJets20(; orderT=6, orderQ=2)) #!jl
+sol_1 = solve(bruss(0.1); T=30.0, alg=alg);  #!jl
+
+#-
 
 fig = plot(; xlab="x", ylab="y", xlims=(0.9, 1.05), ylims=(1.43, 1.57), legend=:bottomright)  #!jl
+plot!(fig, sol_1; vars=(1, 2), lw=0.2, color=:red, lab="r = 0.1", alpha=0.4)  #!jl
+plot!(fig, sol_01; vars=(1, 2), lw=0.2, color=:blue, lab="r = 0.01")  #!jl
 
-plot!(sol_1; vars=(1, 2), lw=0.2, color=:red, lab="r = 0.1", alpha=0.4)  #!jl
+#!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
-plot!(sol_01; vars=(1, 2), lw=0.2, color=:blue, lab="r = 0.01")  #!jl
+#-
 
-#!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
+LazySets.set_ztol(Float64, old_ztol);  # reset precision #!jl
 
-# The volume at time ``T = 9.0`` can be obtained by evaluating the flowpipe and computing the volume of the hyperrectangular overapproximation:
+# The volume at time ``T = 9.0`` can be (over)estimated by evaluating the
+# flowpipe and computing the volume of a hyperrectangular overapproximation:
 
 vol_01 = volume(set(overapproximate(sol_01(9.0), Hyperrectangle)))  #!jl
 
@@ -130,6 +140,11 @@ vol_1 = volume(set(overapproximate(sol_1(9.0), Hyperrectangle)))  #!jl
 
 # ## References
 
-# [^CAS13]: X. Chen, E. Abraham, S. Sankaranarayanan. *Flow*: An Analyzer for Non-Linear Hybrid Systems.* In Proceedings of the 25th International Conference on Computer Aided Verification (CAV’13). Volume 8044 of LNCS, pages 258-263, Springer, 2013.
-
-# [^GCLASG20]: S. Gruenbacher, J. Cyranka, M. Lechner, Md. Ariful Islam,, Scott A. Smolka, R. Grosu *Lagrangian Reachtubes: The Next Generation*. arXiv: 2012.07458
+# [^CAS13]: X. Chen, E. Abraham, S. Sankaranarayanan. *Flow*\*: *An Analyzer for
+#           Non-Linear Hybrid Systems*. In Proceedings of the 25th International
+#           Conference on Computer Aided Verification (CAV'13). Volume 8044 of
+#           LNCS, pages 258-263, Springer, 2013.
+#
+# [^GCLISG20]: S. Gruenbacher, J. Cyranka, M. Lechner, Md. Ariful Islam, Scott
+#              A. Smolka, R. Grosu *Lagrangian Reachtubes: The Next Generation*.
+#              CDC 2020.