From 017ab014e8a6a1cf310d2d0bfe3928b1c5381c90 Mon Sep 17 00:00:00 2001 From: lucaferranti Date: Fri, 8 Jul 2022 11:29:12 +0300 Subject: [PATCH 1/5] add quadratic map for sspz --- src/Sets/SimpleSparsePolynomialZonotope.jl | 44 +++++++++++++++++++++- 1 file changed, 43 insertions(+), 1 deletion(-) diff --git a/src/Sets/SimpleSparsePolynomialZonotope.jl b/src/Sets/SimpleSparsePolynomialZonotope.jl index 19f055e121..0f8cabc2d7 100644 --- a/src/Sets/SimpleSparsePolynomialZonotope.jl +++ b/src/Sets/SimpleSparsePolynomialZonotope.jl @@ -1,5 +1,5 @@ export SimpleSparsePolynomialZonotope, PolynomialZonotope, expmat, nparams, - linear_map + linear_map, quadratic_map """ SimpleSparsePolynomialZonotope{N, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}, ME<:AbstractMatrix{<:Integer}} <: AbstractPolynomialZonotope{N} @@ -152,3 +152,45 @@ The set resulting from applying the linear map `M` to `P`. function linear_map(M::AbstractMatrix, P::SSPZ) return SimpleSparsePolynomialZonotope(M * center(P), M * genmat(P), expmat(P)) end + +""" + quadratic_map(Q::Vector{MT}, S::SimpleSparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}} + +Return an overapproximation of the quadratic map of the given zonotope. + +### Input + +- `S` -- simple sparse polynomial zonotope +- `Q` -- array of square matrices + +### Output + +An overapproximation of the quadratic map of the given zonotope. +""" +function quadratic_map(Q::Vector{MT}, S::SimpleSparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}} + m = length(Q) + c = center(S) + h = ngens(S) + G = genmat(S) + E = expmat(S) + + cnew = similar(c, m) + Gnew = similar(G, m, h^2 + h) + tmp = similar(Q) + @inbounds for (i, q) in enumerate(Q) + cnew[i] = dot(c, q, c) + Gnew[i, 1:h] = c' * (q + q')*G + tmp[i] = q * G + end + + Enew = repeat(E, 1, h + 1) + @inbounds for i in 1:h + idxstart = h * i + 1 + idxend = (i + 1) * h + Enew[:, idxstart:idxend] .+= E[:, i] + for j in eachindex(tmp) + Gnew[j, idxstart:idxend] = G[:, i]' * tmp[j] + end + end + return SimpleSparsePolynomialZonotope(cnew, Gnew, Enew) +end From 67db7d63cd72eb297d5200194ed2e8904009a580 Mon Sep 17 00:00:00 2001 From: lucaferranti Date: Fri, 8 Jul 2022 11:39:15 +0300 Subject: [PATCH 2/5] added docstring and tests --- docs/src/lib/sets/SimpleSparsePolynomialZonotope.md | 1 + test/Sets/SimpleSparsePolynomialZonotope.jl | 12 ++++++++++++ 2 files changed, 13 insertions(+) diff --git a/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md b/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md index bd21404b08..4bcdce9091 100644 --- a/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md +++ b/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md @@ -20,4 +20,5 @@ cartesian_product(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZono linear_combination(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) convex_hull(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) convex_hull(::SimpleSparsePolynomialZonotope) +quadratic_map(::Vector{MT}, ::SimpleSparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}} ``` diff --git a/test/Sets/SimpleSparsePolynomialZonotope.jl b/test/Sets/SimpleSparsePolynomialZonotope.jl index 73e1b863bc..0962d38017 100644 --- a/test/Sets/SimpleSparsePolynomialZonotope.jl +++ b/test/Sets/SimpleSparsePolynomialZonotope.jl @@ -53,4 +53,16 @@ for N in [Float64, Float32, Rational{Int}] @test center(CH1) == _c @test genmat(CH1) == _g @test expmat(CH1) == _e + + S2 = SimpleSparsePolynomialZonotope(N[0.2, -0.6], N[1 0;0 0.4], [1 0;0 1]) + Q = [[0. 1;1 0], [-3.2 1.2;1.2 0]] + + q1 = quadratic_map(Q, S2) + + @test center(q1) ≈ [-0.24, -0.416] + @test genmat(q1) ≈ [-1.2 0.16 0.0 0.4 0.4 0.0; + -2.72 0.192 -3.2 0.48 0.48 0.0] + + @test expmat(q1) == [ 1 0 2 1 1 0; + 0 1 0 1 1 2] end From 2d121fa4043b76508b82fddb713dfffc3c9dde2e Mon Sep 17 00:00:00 2001 From: Luca Ferranti <49938764+lucaferranti@users.noreply.github.com> Date: Sat, 9 Jul 2022 10:42:59 +0300 Subject: [PATCH 3/5] Apply suggestions from code review Co-authored-by: Christian Schilling --- src/Sets/SimpleSparsePolynomialZonotope.jl | 26 +++++++++++++++------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/Sets/SimpleSparsePolynomialZonotope.jl b/src/Sets/SimpleSparsePolynomialZonotope.jl index 0f8cabc2d7..f91759bc1f 100644 --- a/src/Sets/SimpleSparsePolynomialZonotope.jl +++ b/src/Sets/SimpleSparsePolynomialZonotope.jl @@ -160,12 +160,22 @@ Return an overapproximation of the quadratic map of the given zonotope. ### Input +- `Q` -- vector of square matrices - `S` -- simple sparse polynomial zonotope -- `Q` -- array of square matrices ### Output An overapproximation of the quadratic map of the given zonotope. + +### Algorithm + +This method implements Proposition 12 in [1]. +See also Proposition 3.1.30 and Proposition 3.1.31 in [2]. + +[1] Kochdumper, Althoff. *Sparse polynomial zonotopes - a novel set +representation for reachability analysis*. 2021 +[2] Kochdumper. *Extensions of polynomial zonotopes and their application to +verification of cyber-physical systems*. 2021. """ function quadratic_map(Q::Vector{MT}, S::SimpleSparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}} m = length(Q) @@ -176,11 +186,11 @@ function quadratic_map(Q::Vector{MT}, S::SimpleSparsePolynomialZonotope) where { cnew = similar(c, m) Gnew = similar(G, m, h^2 + h) - tmp = similar(Q) - @inbounds for (i, q) in enumerate(Q) - cnew[i] = dot(c, q, c) - Gnew[i, 1:h] = c' * (q + q')*G - tmp[i] = q * G + QiG = similar(Q) + @inbounds for (i, Qi) in enumerate(Q) + cnew[i] = dot(c, Qi, c) + Gnew[i, 1:h] = c' * (Qi + Qi') * G + QiG[i] = Qi * G end Enew = repeat(E, 1, h + 1) @@ -188,8 +198,8 @@ function quadratic_map(Q::Vector{MT}, S::SimpleSparsePolynomialZonotope) where { idxstart = h * i + 1 idxend = (i + 1) * h Enew[:, idxstart:idxend] .+= E[:, i] - for j in eachindex(tmp) - Gnew[j, idxstart:idxend] = G[:, i]' * tmp[j] + for j in eachindex(QiG) + Gnew[j, idxstart:idxend] = G[:, i]' * QiG[j] end end return SimpleSparsePolynomialZonotope(cnew, Gnew, Enew) From cb4517ca62aa34e173ecc6ed7faec8270d94bbc7 Mon Sep 17 00:00:00 2001 From: Luca Ferranti <49938764+lucaferranti@users.noreply.github.com> Date: Sat, 9 Jul 2022 12:03:23 +0300 Subject: [PATCH 4/5] Update src/Sets/SimpleSparsePolynomialZonotope.jl Co-authored-by: Christian Schilling --- src/Sets/SimpleSparsePolynomialZonotope.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Sets/SimpleSparsePolynomialZonotope.jl b/src/Sets/SimpleSparsePolynomialZonotope.jl index f91759bc1f..297290c528 100644 --- a/src/Sets/SimpleSparsePolynomialZonotope.jl +++ b/src/Sets/SimpleSparsePolynomialZonotope.jl @@ -170,11 +170,11 @@ An overapproximation of the quadratic map of the given zonotope. ### Algorithm This method implements Proposition 12 in [1]. -See also Proposition 3.1.30 and Proposition 3.1.31 in [2]. +See also Proposition 3.1.30 in [2]. -[1] Kochdumper, Althoff. *Sparse polynomial zonotopes - a novel set +[1] N. Kochdumper, M. Althoff. *Sparse polynomial zonotopes: A novel set representation for reachability analysis*. 2021 -[2] Kochdumper. *Extensions of polynomial zonotopes and their application to +[2] N. Kochdumper. *Extensions of polynomial zonotopes and their application to verification of cyber-physical systems*. 2021. """ function quadratic_map(Q::Vector{MT}, S::SimpleSparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}} From df694de33f450df59f66fd57e471908135d61ff4 Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Fri, 15 Jul 2022 10:48:24 -0300 Subject: [PATCH 5/5] Update src/Sets/SimpleSparsePolynomialZonotope.jl --- src/Sets/SimpleSparsePolynomialZonotope.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Sets/SimpleSparsePolynomialZonotope.jl b/src/Sets/SimpleSparsePolynomialZonotope.jl index 297290c528..1d008cc8ac 100644 --- a/src/Sets/SimpleSparsePolynomialZonotope.jl +++ b/src/Sets/SimpleSparsePolynomialZonotope.jl @@ -165,7 +165,7 @@ Return an overapproximation of the quadratic map of the given zonotope. ### Output -An overapproximation of the quadratic map of the given zonotope. +The quadratic map of the given zonotope represented as a polynomial zonotope. ### Algorithm