add discrete post for VREP

JuliaReach · Jan 9, 2025 · bf3ccd2 · bf3ccd2
1 parent 0619171
commit bf3ccd2
Show file tree

Hide file tree

Showing 2 changed files with 44 additions and 26 deletions.
diff --git a/src/Algorithms/VREP/post.jl b/src/Algorithms/VREP/post.jl
@@ -1,17 +1,11 @@
-# ===========================
-# Continuous post interface
-# ===========================
-
-# this algorithm uses polygons (two dimensions) or polytopes (any dimension) in vertex representation
-function post(alg::VREP{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
-    @unpack δ, approx_model, static, dim = alg
+# continuous post
+function post(alg::VREP, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
     # normalize system to canonical form
-    # x' = Ax, x in X, or
-    # x' = Ax + u, x in X, u in U
     ivp_norm = _normalize(ivp)
 
     # homogenize the initial-value problem
@@ -20,20 +14,37 @@ function post(alg::VREP{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::VREP{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+    @unpack δ, approx_model, static, dim = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
 
     # the initial state should be expressed in v-rep
-    Ω0 = initial_state(ivp_discr)
+    Ω0 = initial_state(ivp)
     if !(Ω0 isa VPolygon || Ω0 isa VPolytope)
         n = size(Φ, 1)
         Ω0 = n == 2 ? convert(VPolygon, Ω0) : convert(VPolytope, Ω0)
     end
 
-    X = stateset(ivp_discr)
+    X = stateset(ivp)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # reconvert the set of initial states and state matrix, if needed
     Ω0 = _reconvert(Ω0, static, dim)

diff --git a/test/algorithms/VREP.jl b/test/algorithms/VREP.jl
@@ -41,12 +41,19 @@ end
     A = [0 1.0
          -1 0]
     X0 = VPolygon([[1.0, 0.0], [1.2, 0.0], [1.1, 0.2]])
-    prob = @ivp(x' = A * x, x(0) ∈ X0)
-    sol = solve(prob; tspan=(0.0, 2.0), alg=VREP(; δ=1e-2))
+    ivp = @ivp(x' = A * x, x(0) ∈ X0)
+    alg = VREP(; δ=1e-2)
+    # continuous algorithm
+    sol = solve(ivp; tspan=(0.0, 2.0), alg=alg)
     @test setrep(sol) == VPolygon{Float64,Vector{Float64}}
+    # discrete algorithm
+    ivp_norm = ReachabilityAnalysis._normalize(ivp)
+    ivp_discr = discretize(ivp_norm, alg.δ, alg.approx_model)
+    NSTEPS = 500
+    fp_d = ReachabilityAnalysis.post(alg, ivp_discr, NSTEPS)
 
     # statically sized array
-    sol = solve(prob; tspan=(0.0, 2.0), alg=VREP(; δ=1e-2, static=true, dim=2))
+    sol = solve(ivp; tspan=(0.0, 2.0), alg=VREP(; δ=1e-2, static=true, dim=2))
     @test setrep(sol) == VPolygon{Float64,SVector{2,Float64}}
 end
 
@@ -62,21 +69,21 @@ end
 
     # default Polyhedra backend works, but shows warning messages of this sort:
     # glp_simplex: unable to recover undefined or non-optimal solution
-    #prob = @ivp(x' = A*x, x(0) ∈ X0)
-    #sol = solve(prob, tspan=(0.0, 1.0), alg=VREP(δ=1e-3, static=true, dim=4))
+    #ivp = @ivp(x' = A*x, x(0) ∈ X0)
+    #sol = solve(ivp, tspan=(0.0, 1.0), alg=VREP(δ=1e-3, static=true, dim=4))
 
     # change the backend
-    prob = @ivp(x' = A * x, x(0) ∈ X0)
-    sol = solve(prob; tspan=(0.0, 1.0),
+    ivp = @ivp(x' = A * x, x(0) ∈ X0)
+    sol = solve(ivp; tspan=(0.0, 1.0),
                 alg=VREP(; δ=1e-3, static=true, dim=4, backend=CDDLib.Library()))
 end
 
 @testset "VREP for sdof with distinct approximations" begin
-    prob, _ = harmonic_oscillator()
+    ivp, _ = harmonic_oscillator()
     tmax = 2.0
 
-    a = solve(prob; tspan=(0.0, tmax), alg=VREP(; δ=0.1, approx_model=Forward(; setops=:concrete)))
-    b = solve(prob; tspan=(0.0, tmax),
+    a = solve(ivp; tspan=(0.0, tmax), alg=VREP(; δ=0.1, approx_model=Forward(; setops=:concrete)))
+    b = solve(ivp; tspan=(0.0, tmax),
               alg=VREP(; δ=0.1, approx_model=StepIntersect(; setops=:concrete)))
-    c = solve(prob; tspan=(0.0, tmax), alg=VREP(; δ=0.1, approx_model=CorrectionHull(; exp=:base)))
+    c = solve(ivp; tspan=(0.0, tmax), alg=VREP(; δ=0.1, approx_model=CorrectionHull(; exp=:base)))
 end