#912 - Quadratic map of a zonotope (#2218)

* add quadratic map of a zonotope * Update src/Sets/Zonotope.jl Co-authored-by: Christian Schilling <[email protected]> * Update src/Sets/Zonotope.jl Co-authored-by: Christian Schilling <[email protected]> * apply suggestions * change signature * move c' out of the loop * Update src/Sets/Zonotope.jl * Update docs/src/lib/sets/Zonotope.md * add mathematical description * Update src/Sets/Zonotope.jl Co-authored-by: Marcelo Forets <[email protected]> * Update src/Sets/Zonotope.jl Co-authored-by: Marcelo Forets <[email protected]> * add tests, export function and remove zero generators * Update unit_Zonotope.jl * fix docstring * remove constant calculations from inner loop Co-authored-by: Christian Schilling <[email protected]> Co-authored-by: Marcelo Forets <[email protected]>
JuliaReach · Jul 29, 2020 · bf842a5 · bf842a5
1 parent 2568309
commit bf842a5
Show file tree

Hide file tree

Showing 3 changed files with 78 additions and 1 deletion.
diff --git a/docs/src/lib/sets/Zonotope.md b/docs/src/lib/sets/Zonotope.md
@@ -18,6 +18,7 @@ split(::AbstractZonotope, ::Int)
 split(::AbstractZonotope, ::AbstractVector{Int}, ::AbstractVector{Int})
 remove_zero_generators(::Zonotope)
 linear_map!(::Zonotope, ::AbstractMatrix, ::Zonotope)
+quadratic_map(Q::Vector{MT}, Z::Zonotope{N}) where {N, MT<:AbstractMatrix{N}}
 ```
 
 Inherited from [`LazySet`](@ref):

diff --git a/src/Sets/Zonotope.jl b/src/Sets/Zonotope.jl
@@ -3,7 +3,8 @@ import Base: rand
 export Zonotope,
        scale,
        reduce_order,
-       remove_zero_generators
+       remove_zero_generators,
+       quadratic_map
 
 using LazySets.Arrays: _vector_type, _matrix_type
 
@@ -427,3 +428,67 @@ function linear_map!(Zout::Zonotope, M::AbstractMatrix, Z::Zonotope)
     mul!(Zout.generators, M, Z.generators)
     return Zout
 end
+
+"""
+    quadratic_map(Q::Vector{MT}, Z::Zonotope{N}) where {N, MT<:AbstractMatrix{N}}
+
+Return an overapproximation of the quadratic map of the given zonotope.
+
+### Input
+
+- `Z` -- zonotope
+- `Q` -- array of square matrices
+
+### Output
+
+An overapproximation of the quadratic map of the given zonotope.
+
+### Notes
+
+Mathematically, a quadratic map of a zonotope is defined as:
+
+```math
+Z_Q = \\right\\{ \\lambda | \\lambda_i = x^T Q\\^{(i)} x,~i = 1, \\ldots, n,~x \\in Z \\left\\}
+```
+such that each coordinate ``i`` of the resulting zonotope is influenced by ``Q\\^{(i)}``
+
+### Algorithm
+
+This function implements [Lemma 1, 1].
+
+[1] *Matthias Althoff and Bruce H. Krogh. 2012. Avoiding geometric intersection
+operations in reachability analysis of hybrid systems. In Proceedings of the
+15th ACM international conference on Hybrid Systems: Computation and Control
+(HSCC ’12). Association for Computing Machinery, New York, NY, USA, 45–54.*
+"""
+function quadratic_map(Q::Vector{MT}, Z::Zonotope{N}) where {N, MT<:AbstractMatrix{N}}
+    @assert length(Q) == dim(Z) "the number of matrices needs to match the dimension of the zonotope"
+    G = genmat(Z)
+    c = center(Z)
+    n, p = size(G)
+    h = Matrix{N}(undef, n, binomial(p+2, 2)-1)
+    d = Vector{N}(undef, n)
+    g(x) = view(G, :, x)
+    cᵀ = c'
+    for (i, Qᵢ) in enumerate(Q)
+        cᵀQᵢ = cᵀ * Qᵢ
+        Qᵢc = Qᵢ * c
+        aux = zero(N)
+        for j=1:p
+            aux += g(j)' * Qᵢ * g(j)
+            h[i, j] = cᵀQᵢ * g(j) + g(j)' * Qᵢc
+            h[i, p+j] = 0.5 * g(j)' * Qᵢ * g(j)
+        end
+        d[i] = cᵀQᵢ * c + 0.5 * aux
+        l = 0
+        for j=1:p-1
+            gjᵀQᵢ = g(j)' * Qᵢ
+            Qᵢgj = Qᵢ * g(j)
+            for k=j+1:p
+                l += 1
+                h[i, 2p+l] = gjᵀQᵢ * g(k) + g(k)' * Qᵢgj
+            end
+        end
+    end
+    return Zonotope(d, remove_zero_columns(h))
+end
diff --git a/test/unit_Zonotope.jl b/test/unit_Zonotope.jl
@@ -313,4 +313,15 @@ for N in [Float64]
     H2 = Hyperrectangle(low=N[-2, -2], high=N[2, 0])
     @test issubset(Z, H1)
     @test !issubset(Z, H2)
+
+    # quadratic map
+    Z = Zonotope(N[0, 0], N[1 0; 0 1])
+    Q1 = N[1/2 0; 0 1/2]
+    Q2 = N[0 1/2; 1/2 0]
+    # note that there may be repeated generators (though zero generaors are removed)
+    @test quadratic_map([Q1, Q2], Z) == Zonotope(N[0.5, 0], N[0.25 0.25 0; 0 0 1])
+    Z = Zonotope(N[0, 0], N[1 1; 0 1])
+    Q1 = N[1 1; 1 1]
+    Q2 = N[-1 0; 0 -1]
+    @test quadratic_map([Q1, Q2], Z) == Zonotope(N[2.5, -1.5], N[0.5 2 4; -0.5 -1 -2])
 end