Concrete Pontryagin Difference of Polytopes

[ueliwechsler]

Similar to #586 the Pontryagin (or Minkwoski) Difference of two HPolytopes from Theorem 2.2 in Book on RPI sets.

I am not sure about it, but it seems straight forward to implement (...too straightforward 🤔) and if think it should work.

I have done a naive implementation and I'd like to hear your opinion on it.

""" Pontryagin Set Difference for HPolyhedron and HPolygon
    From Theory and Computation of Disturbance Invariant Sets for discrete time
    linear systems https://www.hindawi.com/journals/mpe/1998/934097/abs/"""
Base.:-(X::HPolytope, Y::HPolytope) = begin
    ax, bx, n = get_constraints(X)
    b_xminusy = bx .- ρ.(ax,Ref(Y))
    XminusY = HPolytope(vcat(ax'...), b_xminusy)
    return XminusY
end

function get_constraints(𝒫::AbstractPolyhedron)
    a = getfield.(𝒫.constraints,:a) # normal vector of halfspace
    b = getfield.(𝒫.constraints,:b) # support value of halfspace
    n = length(𝒫.constraints)
    return a,b, n
end

[schillic]

I have not yet looked at the algorithm, so here are just two high-level comments.

1. Not all subtypes of AbstractPolyhedron have a field 𝒫.constraints, but you can obtain the list of constraints via constraints_list(𝒫).
2. What is the effect of Ref(Y), compared to just Y?

[ueliwechsler]

Thanks!

2.What is the effect of Ref(Y), compared to justY?

It tells the broadcast operator, that Yis a single element and, therefore, it iterates only over the vectors in ax.

[mforets]

Similar to #586 the Pontryagin (or Minkwoski) Difference of two HPolytopes from Theorem 2.2 in Book on RPI sets.

Very interesting reference, thanks. Now, given a direction d and two polytopes X and Y and if X ⊖ Y denotes the Minkowski difference between X and Y, Theorem 2.2 on Kolmanovsky and Gilbert claims that ρ(d, X ⊖ Y) ≤ ρ(d, X) - ρ(d, Y). Note that equality doesn't hold -- they show a counterexample in Remark 2.1, namely

julia> X = BallInf(zeros(2), 1.0)
BallInf{Float64}([0.0, 0.0], 1.0)

julia> Y = Ball2(zeros(2), 1/2)
Ball2{Float64}([0.0, 0.0], 0.5)

Then, the approach that i suggest is the following one:

1- define a new wrapper type for the Minkowski difference operator, the lazy Minkowski difference MinkowskiDifference : https://github.com/JuliaReach/LazySets.jl/issues/1509 2- define a new overapproximation function that dispatches on the lazy Minkowski difference of two sets whose constraints_list is available (the implementation of 2 is essentially the one in the OP)

[mforets]

The reason to write it this way is that the operation is now expressed exactly, which can later be approximated in different ways, or nested with other lazy sets in some formula. And we can also add algorithms that compute the minkowski difference without approximation by dispatching on a new function minkowski_difference.

The case of zonotopes is treated in Althoff's On computing the Minkowski difference of zonotopes.

[ueliwechsler]

For the implementation of (2) the fact that ρ(d, X ⊖ Y) ≤ ρ(d, X) - ρ(d, Y) has no consequences, it just results in the fact that X ⊖ Y ⊕ Y ⊆ X. But how would this affect the implementation (1). I assume, for (1), we could use most of the Code from MinkowskiSum. E.g. How would you then calculate ρ(x, md::MinkowskiDifference)?

For (2), I have an adjustment to the OP. I think, we could use tosimplehrep. Do think this would work for X::LazySet in general?

Base.:-(X::AbstractPolytope, Y::LazySet) = begin
    # ax, bx = get_constraints(X)
    F, g = tosimplehrep(X)
    g_xminusy = [g[i] - ρ(F[i,:],Y) for i in eachindex(g)]
    XminusY = HPolytope(F, vcat(g_xminusy...))
    return XminusY
end

[mforets]

Ok! Looking more carefully in the following page (page 326), i see now that Theorem 2.3 says that if X is a polyhedron then X ⊖ Y can be computed by a finite intersection of halfspaces, so i agree with the behavior of your proposed function.

(Approximating a lazy MinkowskiDifference still can be added for a general LazySet, although it will only be an overapproximation).

I leave some comments below, feel free to submit a PR ! 👍

In LazySets we would rather name this function in lower case minkowski_difference, leaving - as an alias to the lazy constructor MinkowskiDifference to be added in #1509 (compare between minkowski_sum and +, or ⊕, the latter which are just aliases for MinkowskiSum).

Do think this would work for X::LazySet in general?

X::AbstractPolytope LGTM, stating clearly in the docstring that both X and Y are assumed to be bounded (or add assertions like @assert isbounded(Y)). We can then generalize to possibly unbounded X and Y in a second stage. (*)

You can add the new function in AbstractPolytope.jl.

Regarding the implementation, it works without the splatting too and getting rid of it gives a slight performance improvement, compare with:

function minkowski_difference(X::AbstractPolytope, Y::LazySet)
    F, g = tosimplehrep(X)
    g_xminusy = [g[i] - ρ(F[i, :], Y) for i in eachindex(g)]
    XminusY = HPolytope(F, g_xminusy)
    return XminusY
end

---

(*) Remark 2.4 (page 327) states that neither X nor Y need to be bounded for the theorem to hold. There is also a discussion about the case of Y being unbounded: if Y is unbounded in one of the directions normal to the faces of X then the difference is the empty set. We can dispatch on unbounded sets in LazySets using the AbstractPolyhedron interface, and unboundedness of the support function can be checked too, see e.g. this function.

[schillic]

We can dispatch on unbounded sets in LazySets using the AbstractPolyhedron interface

Y<:AbstractPolyhedron does not imply that Y is unbounded.

[ueliwechsler]

When talking to different people about the set difference topic, I found out, that the definition of the Minkowski difference is ambiguous.

Therefore, I propose the distinction between the following two differences.

1. Pontryagin Difference (as described above and accordingly implemented) as A⊖B={x|x+b∈A|∀b∈B}
2. Minkowski difference (which is, e.g. used in collision detection) as A⊖B={a−b|a∈A,b∈B} = A ⊕ (-B) And can be computed by just taking the Minkowski sum of the first set with the "negative" second set.

[mforets]

Thanks for sharing the thought. Indeed, in the previous issue we saw that some people mean different things for Minkowski difference. My opinion would be that we reserve the name MinkowskiDifference for the Pontryagin (aka geometric difference), since it is the meaning used eg. in control theory applications which is closer to the application domain of LazySets.

However, it's good to handle the second meaning too; for the lazy reflection of a set wrt the origin i've opened https://github.com/JuliaReach/LazySets.jl/issues/1558, so the operation can be built lazily. We can also add concrete counterparts reflect(...) (see https://github.com/JuliaReach/LazySets.jl/issues/1559) to be used as f(A, B) = minkowski_sum(A, reflect(B)), such that a downstream package can name this new function to its convenience.

[ueliwechsler]

I agree, it creates some confusion and there are ambiguities.

I would prefer the name pontryagin difference 😉, but I am totally fine with your proposition as long as ⊖ respectively - does the right operation.

[mforets]

i like the idea of adding an alias