Merge pull request #931 from JuliaReach/test_rem_gens

schillic · web-flow · commit 4292e7612d41 · 2025-09-01T18:33:46.000+02:00
Add reachability algorithm `HLBS25` for linear parametrix systems
diff --git a/Project.toml b/Project.toml
@@ -31,7 +31,7 @@ IntervalArithmetic = "0.16 - 0.20, =0.20.9"  # new versions require updates and
 IntervalMatrices = "0.11"
 LazySets = "2.14, 3, 4, 5"
 LinearAlgebra = "<0.0.1, 1.6"
-MathematicalSystems = "0.11 - 0.14"
+MathematicalSystems = "0.14.1"
 Parameters = "0.10 - 0.12"
 Random = "<0.0.1, 1.6"
 ReachabilityBase = "0.3.1"
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -24,7 +24,7 @@ JLD2 = "0.4 - 0.6"
 LaTeXStrings = "1"
 LazySets = "2.14, 3, 4, 5"
 Literate = "2"
-MathematicalSystems = "0.11 - 0.14"
+MathematicalSystems = "0.14.1"
 OrdinaryDiffEq = "6"
 Plots = "1"
 ReachabilityBase = "0.2.3 - 0.3"
diff --git a/docs/make.jl b/docs/make.jl
@@ -77,6 +77,12 @@ TAYLOR_METHODS = ["Taylor methods" => [
 #                                                             #
 #                                                             ]]
 
+PARAMETRIC_SYSTEMS = ["Parametric systems" => [
+                                                #
+                                                "Parametric reachability" => "man/parametric.md"
+                                                #
+                                                ]]
+
 HYBRID_SYSTEMS = ["Hybrid systems" => [
                                        #
                                        #"Introduction"           => "tutorials/hybrid_systems/introduction.md",
@@ -153,6 +159,18 @@ NONLINEAR_SOLVERS_API = [
                          #
                          ]
 
+PARAMETRIC_SOLVERS = [
+                      #
+                      "HLBS25" => "man/algorithms/HLBS25.md"
+                      #
+                      ]
+
+PARAMETRIC_SOLVERS_API = [
+                          #
+                          "HLBS25" => "lib/algorithms/HLBS25.md"
+                          #
+                          ]
+
 # ========================
 # Docs contents
 # ========================
@@ -173,7 +191,7 @@ makedocs(; sitename="ReachabilityAnalysis.jl",
                                     LINEAR_METHODS,
                                     UNCERTAIN_INPUTS,
                                     TAYLOR_METHODS,
-                                    #UNCERTAIN_PARAMETERS,
+                                    PARAMETRIC_SYSTEMS,
                                     #LINEAR_PDE,
                                     HYBRID_SYSTEMS,
                                     CLOCKED_SYSTEMS
@@ -184,7 +202,8 @@ makedocs(; sitename="ReachabilityAnalysis.jl",
                              "Systems" => "man/systems.md",
                              #"Reach-sets"        => "man/reachsets.md",
                              #"Flowpipes"         => "man/flowpipes.md",
-                             "Linear solvers" => LINEAR_SOLVERS
+                             "Linear solvers" => LINEAR_SOLVERS,
+                             "Parametric solvers" => PARAMETRIC_SOLVERS
                              #"Nonlinear solvers" => NONLINEAR_SOLVERS,
                              #"Solutions"         => "man/solutions.md",
                              #"Invariants"        => "man/invariants.md",
@@ -221,6 +240,7 @@ makedocs(; sitename="ReachabilityAnalysis.jl",
                                     "Flowpipes"                     => "lib/flowpipes.md",
                                     "Linear solvers"                => LINEAR_SOLVERS_API,
                                     "Nonlinear solvers"             => NONLINEAR_SOLVERS_API,
+                                     "Parametric solvers"            => PARAMETRIC_SOLVERS_API,
                                     "Solutions"                     => "lib/solutions.md",
                                     "Discretization"                => "lib/discretize.md",
                                     "Projections"                   => "lib/projections.md",
diff --git a/docs/src/lib/algorithms/HLBS25.md b/docs/src/lib/algorithms/HLBS25.md
@@ -0,0 +1,4 @@
+# HLBS25
+```@docs
+HLBS25
+```
diff --git a/docs/src/lib/discretize.md b/docs/src/lib/discretize.md
@@ -29,6 +29,7 @@ FirstOrderZonotope
 SecondOrderddt
 CorrectionHull
 StepIntersect
+CorrectionHullMatrixZonotope
 ```
 
 ## Exponentiation
diff --git a/docs/src/man/algorithms/HLBS25.md b/docs/src/man/algorithms/HLBS25.md
@@ -0,0 +1,5 @@
+# Reachability of parametric linear systems (HLBS25)
+
+## Method
+
+This method computes the reachability analysis of linear systems with parametric uncertainty. The method is based on matrix zonotopes to represent the uncertain parameters of the system.
diff --git a/docs/src/man/parametric.md b/docs/src/man/parametric.md
@@ -4,3 +4,15 @@ CurrentModule = ReachabilityAnalysis
 ```
 
 # Parametric reachability
+
+The reachability analysis of parametric systems is performed by the `solve` function.
+The following algorithms are available:
+
+```@contents
+Pages = ["HLBS25.md"]
+Depth = 3
+```
+
+```@meta
+CurrentModule = ReachabilityAnalysis
+```
diff --git a/docs/src/refs.bib b/docs/src/refs.bib
@@ -687,3 +687,18 @@ @inproceedings{ForetsS22
   url          = {https://doi.org/10.1007/978-3-031-07727-2_9},
   doi          = {10.1007/978-3-031-07727-2_9}
 }
+
+@article{HuangLBS2025,
+  author      = {Yushen Huang and
+                 Ertai Luo and
+                 Stanley Bak and
+                 Yifan Sun},
+  title       = {Reachability analysis for linear systems with uncertain parameters using polynomial zonotopes},
+  journal     = {Nonlinear Analysis: Hybrid Systems},
+  volume      = {56},
+  pages       = {101571},
+  year        = {2025},
+  issn        = {1751-570X},
+  url         = {https://www.sciencedirect.com/science/article/pii/S1751570X24001080},
+  doi         = {10.1016/j.nahs.2024.101571}
+}
diff --git a/src/Algorithms/HLBS25/HLBS25.jl b/src/Algorithms/HLBS25/HLBS25.jl
@@ -0,0 +1,58 @@
+import IntervalArithmetic as IA
+
+"""
+    HLBS25{N, AM, RM, R} <: AbstractContinuousPost
+
+Implementation of the reachability algorithm for linear systems with parametric
+uncertainty using matrix zonotopes by [HuangLBS25](@citet).
+
+### Fields
+
+- `δ`                -- step-size of the discretization
+- `approx_model`     -- (optional, default: `CorrectionHullMatrixZonotope`)
+                        model for the discretization of the initial value problem
+- `taylor_order`     -- (optional, default: `5`) order of the Taylor series
+                        expansion of the matrix exponential for each step
+- `max_order`        -- (optional, default: `5`) maximum order of the reach set
+- `reduction_method` -- (optional, default: `LazySets.GIR05()` zonotope
+                        order reduction method used
+- `recursive`        -- (optional, default: `false`) if `true`, compute the
+                        Taylor series expansion of the matrix zonotope
+                        exponential map recursively
+            
+### Notes
+
+The `recursive` option is used to compute the Taylor expansion of the matrix zonotope exponential map.
+If `recursive == true`, each term of the Taylor expansion is computed recursively (e.g., ``A^2 P = A (A P)``).
+
+If `recursive == false`, the Taylor expansion is computed by overapproximating the matrix zonotope exponential 
+map, producing a single matrix that represents the exponential.
+"""
+struct HLBS25{N,AM,RM,R} <: AbstractContinuousPost
+    δ::N
+    approx_model::AM
+    taylor_order::Int
+    max_order::Int
+    reduction_method::RM
+    recursive::R
+end
+
+function HLBS25(; δ::N,
+                approx_model::AM=CorrectionHullMatrixZonotope(),
+                taylor_order::Int=5,
+                max_order::Int=5,
+                reduction_method::RM=LazySets.GIR05(),
+                recursive::Bool=false) where {N,AM,RM}
+    return HLBS25{N,AM,RM,Val{recursive}}(δ, approx_model, taylor_order, max_order,
+                                          reduction_method, Val(recursive))
+end
+
+step_size(alg::HLBS25) = alg.δ
+numtype(::HLBS25{N}) where {N} = N
+
+function rsetrep(::HLBS25{N}) where {N}
+    return ReachSet{N,SparsePolynomialZonotope{N,Matrix{N},Matrix{N},Matrix{Int},Vector{Int}}}
+end
+
+include("post.jl")
+include("reach_homog.jl")
diff --git a/src/Algorithms/HLBS25/post.jl b/src/Algorithms/HLBS25/post.jl
@@ -0,0 +1,44 @@
+# continuous post
+function post(alg::HLBS25, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
+    NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
+
+    # discretize system
+    ivp_discr = discretize(ivp, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::HLBS25{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+    @unpack δ, approx_model, taylor_order, max_order, reduction_method, recursive = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
+    Ω0 = initial_state(ivp)
+
+    # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
+    got_homogeneous = !hasinput(ivp)
+
+    # preallocate output flowpipe
+    ZT = typeof(Ω0)
+    F = Vector{ReachSet{N,ZT}}(undef, NSTEPS)
+
+    if got_homogeneous
+        reach_homog_HLBS25!(F, Ω0, Φ, NSTEPS, δ, taylor_order, recursive, max_order,
+                            reduction_method, Δt0)
+    else
+        error("inhomogeneous algorithm not implemented yet")
+    end
+
+    return Flowpipe(F)
+end
diff --git a/src/Algorithms/HLBS25/reach_homog.jl b/src/Algorithms/HLBS25/reach_homog.jl
@@ -0,0 +1,63 @@
+function reach_homog_HLBS25!(F::Vector{ReachSet{N,S}},
+                             Ω0::S,
+                             Φ::MatrixZonotope{N,MN},
+                             NSTEPS::Integer,
+                             δ::N,
+                             taylor_order::Integer,
+                             recursive::Val{false},
+                             max_order::Integer,
+                             reduction_method::AbstractReductionMethod,
+                             Δt0::IA.Interval) where {N,S<:SparsePolynomialZonotope{N},
+                                                      MN<:AbstractMatrix{N}}
+    # initial reach set
+    Δt = (zero(N) .. δ) + Δt0
+    @inbounds F[1] = ReachSet(Ω0, Δt)
+
+    expΦδ = MatrixZonotopeExp(scale(δ, Φ))
+    expΦδ_approx = overapproximate(expΦδ, MatrixZonotope, taylor_order)
+
+    j = 1
+    @inbounds while j < NSTEPS
+        Zⱼ = set(F[j])
+        Zⱼ₊₁ = overapproximate(expΦδ_approx * Zⱼ, SparsePolynomialZonotope)
+        Zⱼ₊₁ʳ = reduce_order(Zⱼ₊₁, max_order, reduction_method)
+
+        j += 1
+        Δt += δ
+        F[j] = ReachSet(Zⱼ₊₁ʳ, Δt)
+    end
+    return F
+end
+
+function reach_homog_HLBS25!(F::Vector{ReachSet{N,S}},
+                             Ω0::S,
+                             Φ::MatrixZonotope{N,MN},
+                             NSTEPS::Integer,
+                             δ::N,
+                             taylor_order::Integer,
+                             recursive::Val{true},
+                             max_order::Integer,
+                             reduction_method::AbstractReductionMethod,
+                             Δt0::IA.Interval) where {N,S<:SparsePolynomialZonotope{N},
+                                                      MN<:AbstractMatrix{N}}
+    # initial reach set
+    Δt = (zero(N) .. δ) + Δt0
+    @inbounds F[1] = ReachSet(Ω0, Δt)
+
+    expΦδ = MatrixZonotopeExp(scale(δ, Φ))
+
+    j = 1
+    @inbounds while j < NSTEPS
+        Zⱼ = set(F[j])
+
+        em = ExponentialMap(expΦδ, Zⱼ)
+
+        Zⱼ₊₁ = overapproximate(em, SparsePolynomialZonotope, taylor_order)
+        Zⱼ₊₁ʳ = reduce_order(Zⱼ₊₁, max_order, reduction_method)
+
+        j += 1
+        Δt += δ
+        F[j] = ReachSet(Zⱼ₊₁ʳ, Δt)
+    end
+    return F
+end
diff --git a/src/Discretization/CorrectionHullMatrixZonotope.jl b/src/Discretization/CorrectionHullMatrixZonotope.jl
@@ -0,0 +1,81 @@
+module CorrectionHullMatrixZonotopeModule
+
+using LazySets, Reexport
+using MathematicalSystems
+using ..DiscretizationModule
+using LinearAlgebra: I
+
+export CorrectionHullMatrixZonotope
+
+@reexport import ..DiscretizationModule: discretize
+
+const LPCS = LinearParametricContinuousSystem
+const LPDS = LinearParametricDiscreteSystem
+
+"""
+    CorrectionHullMatrixZonotope{R} <: AbstractApproximationModel
+
+Implementation of the discretization algorithm for linear systems with parametric
+uncertainty using matrix zonotopes by [HuangLBS25](@citet).
+
+### Fields
+
+- `taylor_order`     -- (optional, default: `5`) order of the Taylor series
+                        expansion of the matrix exponential for each step
+- `recursive`        -- (optional, default: `false`) if `true`, compute the
+                        Taylor series expansion of the matrix zonotope
+                        exponential map recursively
+
+### Notes
+
+The `recursive` option is used to compute the Taylor expansion of the matrix zonotope exponential map.
+If `recursive == true`, each term of the Taylor expansion is computed recursively (e.g., ``A^2 P = A (A P)``).
+
+If `recursive == false`, the Taylor expansion is computed by overapproximating the matrix zonotope exponential 
+map, producing a single matrix that represents the exponential.
+
+"""
+struct CorrectionHullMatrixZonotope{R} <: AbstractApproximationModel
+    taylor_order::Int
+    recursive::R
+end
+
+function CorrectionHullMatrixZonotope(; taylor_order::Int=5, recursive::Bool=false)
+    return CorrectionHullMatrixZonotope{Val{recursive}}(taylor_order, Val(recursive))
+end
+
+function discretize(ivp::IVP{<:LPCS,<:SparsePolynomialZonotope}, δ,
+                    alg::CorrectionHullMatrixZonotope)
+    taylor_order = alg.taylor_order
+    A = state_matrix(ivp)
+    X0 = initial_state(ivp)
+    n = dim(X0)
+
+    N = eltype(X0) 
+
+    IDₜ = ngens(A) > 0 ? maximum(indexvector(A)) + 1 : 1
+    Tₜ = N(0.5) * δ * Matrix(N(1) * I, n, n)
+    T = MatrixZonotope(Tₜ, [Tₜ], [IDₜ])
+
+    Ω0 = _discretize_CHMZ(A, T, X0, taylor_order, alg.recursive)
+
+    Sdis = LPDS(A)
+    return InitialValueProblem(Sdis, Ω0)
+end
+
+# Recursive case
+function _discretize_CHMZ(A, T, X0, taylor_order, ::Val{true})
+    expAT = MatrixZonotopeExp(A * T)
+    em = ExponentialMap(expAT, X0)
+    return overapproximate(em, SparsePolynomialZonotope, taylor_order)
+end
+
+# Non-recursive case
+function _discretize_CHMZ(A, T, X0, taylor_order, ::Val{false})
+    AT = overapproximate(A * T, MatrixZonotope)
+    expAT = MatrixZonotopeExp(AT)
+    expAT_approx = overapproximate(expAT, MatrixZonotope, taylor_order)
+    return overapproximate(expAT_approx * X0, SparsePolynomialZonotope)
+end
+
+end # module
diff --git a/src/Discretization/DiscretizationModule.jl b/src/Discretization/DiscretizationModule.jl
@@ -70,9 +70,10 @@ in the docstring of each method. Here is a list of all the available approximati
 
 ```jldoctest
 julia> subtypes(ReachabilityAnalysis.DiscretizationModule.AbstractApproximationModel)
-9-element Vector{Any}:
+10-element Vector{Any}:
  Backward
  CorrectionHull
+ CorrectionHullMatrixZonotope
  FirstOrder
  FirstOrderZonotope
  Forward
diff --git a/src/Discretization/discretization.jl b/src/Discretization/discretization.jl
diff --git a/src/Initialization/exports.jl b/src/Initialization/exports.jl
diff --git a/src/Initialization/init.jl b/src/Initialization/init.jl
diff --git a/src/ReachabilityAnalysis.jl b/src/ReachabilityAnalysis.jl
diff --git a/test/algorithms/HLBS25.jl b/test/algorithms/HLBS25.jl
diff --git a/test/runtests.jl b/test/runtests.jl

-Original file line number
+Diff line change
@@ @@ -0,0 +1,4 @@ @@
 +# HLBS25
 +```@docs
 +HLBS25
 +```