Added JKS16 submethod for order reduction.

Co-authored-by: Christian Schilling <[email protected]>

Author: crisiumnih

docs/Project.toml

@@ -38,4 +38,4 @@ RecipesBase = "1"
 StaticArrays = "1"
 SymEngine = "0.7 - 0.12"
 Symbolics = "6.1"
-TaylorModels = "0.6 - 0.7"
+TaylorModels = "0.6 - 0.8"

docs/src/lib/interfaces/AbstractZonotope.md

@@ -243,6 +243,7 @@ LazySets.AbstractReductionMethod
 LazySets.ASB10
 LazySets.COMB03
 LazySets.GIR05
+LazySets.JKS16
 ```
 
 ## Implementations

docs/src/refs.bib

@@ -215,6 +215,20 @@ @article{AlthoffSB10
   doi={10.1016/j.nahs.2009.03.009}
 }
 
+@article{SCOTT2016126,
+  title = {Constrained zonotopes: A new tool for set-based estimation and fault detection},
+  journal = {Automatica},
+  volume = {69},
+  pages = {126-136},
+  year = {2016},
+  issn = {0005-1098},
+  doi = {https://doi.org/10.1016/j.automatica.2016.02.036},
+  url = {https://www.sciencedirect.com/science/article/pii/S0005109816300772},
+  author = {Joseph K. Scott and Davide M. Raimondo and Giuseppe Roberto Marseglia and Richard D. Braatz},
+  keywords = {State estimation, Fault detection, Set-based computing, Zonotopes, Reachability analysis},
+  abstract = {This article introduces a new class of sets, called constrained zonotopes, that can be used to enclose sets of interest for estimation and control. The numerical representation of these sets is sufficient to describe arbitrary convex polytopes when the complexity of the representation is not limited. At the same time, this representation permits the computation of exact projections, intersections, and Minkowski sums using very simple identities. Efficient and accurate methods for computing an enclosure of one constrained zonotope by another of lower complexity are provided. The advantages and disadvantages of these sets are discussed in comparison to ellipsoids, parallelotopes, zonotopes, and convex polytopes in halfspace and vertex representations. Moreover, extensive numerical comparisons demonstrate significant advantages over other classes of sets in the context of set-based state estimation and fault detection.}
+}
+
 @inproceedings{FrehseGDCRLRGDM11,
   author       = {Goran Frehse and
                   Colas Le Guernic and

src/Approximations/overapproximate.jl

@@ -603,6 +603,84 @@ function overapproximate(P::VPolygon, ::Type{<:LinearMap{N,<:Hyperrectangle}}) w
     return LinearMap(R, min_rectangle)
 end
 
+"""
+    overapproximate(P::SparsePolynomialZonotope{N}, ::Type{<:VPolytope}) where {N}
+
+
+Overapproximate a sparse polynomial zonotope with a polytope in vertex representation.
+
+### Input
+
+- `P`         -- sparse polynomial zonotope
+- `VPolytope` -- target type
+
+### Output
+
+A `VPolytope` that overapproximates the sparse polynomial zonotope.
+
+### Algorithm
+
+This method implements [Kochdumper21a; Proposition 3.1.15](@citet).
+The idea is to split `P` into a linear and nonlinear part (such that `P = P₁ ⊕ P₂`). 
+The nonlinear part is enclosed by a zonotope. Then we combine the vertices
+of both sets and finally apply a convex-hull algorithm.
+"""
+function overapproximate(P::SparsePolynomialZonotope{N}, ::Type{<:VPolytope}) where {N}
+    c = center(P)
+    G = genmat_dep(P)
+    GI = genmat_indep(P)
+    E = expmat(P)
+    idx = P.idx
+
+    H = [j for j in 1:size(E, 2) if any(E[:, j] .> 1)]
+    K = setdiff(1:size(E, 2), H)
+
+    if !isempty(H)
+        SPZ₂ = SparsePolynomialZonotope(c, G[:, H], zeros(N, length(c), 0), E[:, H], idx)
+        Z = overapproximate(SPZ₂, Zonotope)
+        c_z = center(Z)
+        GI_mod = hcat(GI, genmat(Z))
+    else
+        c_z = c
+        GI_mod = GI
+    end
+
+    G_mod = G[:, K]
+    E_mod = E[:, K]
+    
+    # P̄ = SparsePolynomialZonotope(c_z, G_mod, GI_mod, E_mod, idx)
+    # Compute vertices of a Z-representation
+    p = size(E, 1)
+    dep_params = Iterators.product(fill([-one(N), one(N)], p)...)
+    indep_params = Iterators.product(fill([-one(N), one(N)], size(GI_mod, 2))...)
+
+    V = Vector{Vector{N}}()
+    for α in dep_params
+        dep_term = zeros(N, size(c))
+        for j in axes(G_mod, 2)
+            prod = one(N)
+            for k in 1:p
+                if E_mod[k, j] == 1
+                    prod *= α[k]
+                end
+            end
+            dep_term += prod * G_mod[:, j]
+        end
+        for β in indep_params
+            indep_term = zeros(N, size(c))
+            for j in axes(GI_mod, 2)
+                indep_term += β[j] * GI_mod[:, j]
+            end
+            point = c_z + dep_term + indep_term
+            push!(V, point)
+        end
+    end
+
+    convex_hull!(V)
+
+    return VPolytope(V)
+end
+
 # function to be loaded by Requires
 function load_paving_overapproximation()
     return quote

src/Interfaces/AbstractZonotope.jl

@@ -808,6 +808,25 @@ Zonotope order-reduction method from [AlthoffSB10](@citet).
 """
 struct ASB10 <: AbstractReductionMethod end
 
+"""
+    JKS16 <: AbstractReductionMethod
+
+`JKS16` uses a greedy factorization and iterative generator reduction
+
+### Algorithm
+
+- The JKS16 method reorders the generator matrix using reduced row echelon form (rref) to form a `T`, 
+    then iteratively removes excess generators from `V` while updating `T`.
+- The default `JKS16()` uses `ϵ = 1e-6` (pivot threshold) and `δ = 1e-3` (volume threshold).
+Referenced from [SCOTT2016126](@citet).
+"""
+struct JKS16{N<:Number} <: AbstractReductionMethod 
+    ϵ::N
+    δ::N
+
+    JKS16(ϵ::N=1e-6, δ::N=1e-3) where {N<:Number} = new{N}(ϵ, δ)
+end
+
 """
     reduce_order(Z::AbstractZonotope, r::Real,
                  [method]::AbstractReductionMethod=GIR05())
@@ -831,16 +850,18 @@ The available algorithms are:
 
 ```jldoctest; setup = :(using LazySets: subtypes, AbstractReductionMethod)
 julia> subtypes(AbstractReductionMethod)
-3-element Vector{Any}:
+4-element Vector{Any}:
  LazySets.ASB10
  LazySets.COMB03
  LazySets.GIR05
+ LazySets.JKS16
 ```
 
 See the documentation of each algorithm for references. These methods split the
 given zonotopic set `Z` into two zonotopes, `K` and `L`, where `K` contains the
 most "representative" generators and `L` contains the generators that are
-reduced, `Lred`, using a box overapproximation. We follow the notation from
+reduced, `Lred`, using a box overapproximation. This methodology varies slightly 
+for `JKS16` (for details, refer to the method). We follow the notation from
 [YangS18](@citet). See also [KopetzkiSA17](@citet).
 """
 function reduce_order(Z::AbstractZonotope, r::Real,
@@ -852,7 +873,10 @@ function reduce_order(Z::AbstractZonotope, r::Real,
 
     # if r is bigger than the order of Z => do not reduce
     (r * n >= p) && return Z
+    return _reduce_order_zonotope_common(Z, r, n, p, method)
+end
 
+function _reduce_order_zonotope_common(Z, r, n, p, method::Union{ASB10,COMB03,GIR05})
     c = center(Z)
     G = genmat(Z)
 
@@ -878,6 +902,138 @@ function reduce_order(Z::AbstractZonotope, r::Real,
     return Zonotope(c, Gred)
 end
 
+function _reduce_order_zonotope_common(Z, r, n, p, method::JKS16)
+    c = center(Z)
+    G = genmat(Z)
+
+    G, G✶ = _factorG(G, method.ϵ, method.δ)
+
+    T = G[:, 1:n]
+    V = G[:, (n+1):end]
+    R = G✶[:, (n+1):end]
+
+    m = floor(Int, n * r)
+
+    while p > m
+        error_metric = [prod(1 .+ abs.(R[:, j])) - (1 + sum(abs.(R[:, j]))) 
+                  for j in 1:size(R, 2)]
+        
+        j_min = argmin(error_metric)
+
+        r_j = R[:, j_min]
+
+        V = V[:, [1:(j_min-1); (j_min+1):end]]
+        R = R[:, [1:(j_min-1); (j_min+1):end]]
+
+        diag_matrix = Diagonal(1 .+ abs.(r_j))
+        T = T * diag_matrix
+        if !isempty(R)
+            R = inv(diag_matrix) * R
+        end
+
+        p -= 1
+    end
+
+    G_red = hcat(T, V)
+    return Zonotope(c, G_red)
+end
+
+# reorder the generator matrix G
+function _factorG(G::AbstractMatrix, ϵ::N, δ::N) where {N}
+    G✶ = copy(G)
+    n, ng = size(G)
+
+    # rref with column swap
+    for k in 1:n
+        # normalize rows k:n
+        for i in k:n
+            row_norm = norm(G✶[i, :], 1)
+            if row_norm > ϵ
+                G✶[i, :] ./= row_norm
+            end
+        end
+
+        submatrix = G✶[k:n, k:ng]
+        max_val, idx = findmax(abs.(submatrix))
+        i_max = k + idx[1] - 1
+        j_max = k + idx[2] - 1
+
+        if abs(max_val) ≤ ϵ
+            break
+        end
+
+        # swap rows and columns
+        if i_max != k
+            G✶[[k, i_max], :] = G✶[[i_max, k], :]
+        end
+
+        if j_max != k
+            G[:, [k, j_max]] = G[:, [j_max, k]]
+            G✶[:, [k, j_max]] = G✶[:, [j_max, k]]
+        end
+
+        G✶ = _rref(G✶)
+    end
+
+    #  extra column swaps until all |g∗ij| ≤ 1 + δ
+    while true
+        (max_val, max_idx) = findmax(abs.(G✶))
+        (max_val ≤ 1 + δ) && break
+        k, j = Tuple(max_idx)
+
+        G[:,[k,j]] = G[:,[j,k]]
+        G✶[:,[k,j]] = G✶[:,[j,k]]
+
+        pivot = G✶[k,k]
+        for i in 1:n
+            i == k && continue
+            factor = G✶[i,k] / pivot
+            G✶[i,:] -= factor .* G✶[k,:]
+        end
+        G✶[k,:] ./= pivot
+    end
+
+    return G, G✶
+end
+
+# reduced row echelon (rref)
+function _rref(G::AbstractMatrix{N}) where {N}
+    A = copy(G)
+    nr, nc = size(A)
+    pivot = 1
+
+    for r in 1:nr
+        if pivot > nc
+            return A
+        end
+
+        i = r
+        while i ≤ nr && A[i, pivot] == 0
+            i += 1
+            if i > nr
+                i = r
+                pivot += 1
+                pivot > nc && return A
+            end
+        end
+
+        A[[i, r], :] = A[[r, i], :]
+
+        if A[r, pivot] != 0
+            A[r, :] ./= A[r, pivot]
+        end
+
+        for i in 1:nr
+            i == r && continue
+            A[i, :] -= A[i, pivot] * A[r, :]
+        end
+
+        pivot += 1
+    end
+
+    return A
+end
+
 # approximate with a box
 function _approximate_reduce_order(c, G, indices, ::Union{COMB03,GIR05})
     return _interval_hull(G, indices)

src/Plotting/plot_recipes.jl

@@ -373,14 +373,12 @@ julia> plot(B, 1e-2)  # faster but less accurate than the previous call
         else
             x, y = res
             if length(x) == 1 ||
-               (length(x) == 2 && norm([x[1], y[1]] - [x[2], y[2]]) ≈ 0)
+               (length(x) == 2 && isapproxzero(norm([x[1], y[1]] - [x[2], y[2]])))
                 # single point
                 seriestype := :scatter
             elseif length(x) == 2
                 # flat line segment
                 linecolor --> DEFAULT_COLOR
-                markercolor --> DEFAULT_COLOR
-                markershape --> :circle
                 seriestype := :path
             else
                 seriestype := :shape

test/Approximations/overapproximate.jl

@@ -147,6 +147,11 @@ for N in [Float64, Rational{Int}, Float32]
     if N <: AbstractFloat
         @test isequivalent(P, R)
     end
+
+    #overapproximate a SparsePolynomialZonotope with a VPolytope
+    S = SparsePolynomialZonotope(N[-0.5, -0.5], N[1. 1 1 1;1 0 -1 1], zeros(N, 2, 0), [1 0 1 2;0 1 1 0])
+    P = overapproximate(S, VPolytope)
+    @test ispermutation([N[-1.5, -2.5], N[2.5, -0.5], N[3.5, 0.5], N[-0.5, 2.5], N[-1.5, 1.5]], vertices_list(P))
 end
 
 # tests that do not work with Rational{Int}

test/Project.toml

@@ -55,5 +55,5 @@ SetProg = "0.3 - 0.4"
 StaticArrays = "0.12, 1"
 SymEngine = "0.7 - 0.12"
 Symbolics = "1 - 5.30, 6.1"
-TaylorModels = "0.0.1, 0.1 - 0.7"
+TaylorModels = "0.0.1, 0.1 - 0.8"
 WriteVTK = "1"

test/Sets/Zonotope.jl

@@ -141,7 +141,7 @@ for N in [Float64, Rational{Int}, Float32]
     @test ispermutation(collect(generators(Z)), [genmat(Z)[:, j] for j in 1:ngens(Z)])
 
     # test order reduction
-    for method in [LazySets.ASB10(), LazySets.COMB03(), LazySets.GIR05()]
+    for method in [LazySets.ASB10(), LazySets.COMB03(), LazySets.GIR05(), LazySets.JKS16()]
         Z = Zonotope(N[2, 1], N[-1//2 3//2 1//2 1; 1//2 3//2 1 1//2])
         Zred1 = reduce_order(Z, 1, method)
         @test ngens(Zred1) == 2