add constraint to CartesianProduct

[yupbank]

hey, I'm new to Julia, but I'm trying to add a constraint to cartesian product and it's not very easy.

What I was trying to do is to express this:

A zonotope Z \subset R^2, (z_1, z_2) \in BZ = Z \times Z \subset R^4 s.t. z_1 + z_2 \in Z.

I can easily get Z \times Z, but to add the z_1 + z_2 \in Z isn't easy for me, can anyone point me to the right direction?

[mforets]

I guess the "easy" part is:

julia> Z = rand(Zonotope, dim=2)
Zonotope{Float64,Array{Float64,1},Array{Float64,2}}([-0.5237041738184478, -0.6673040904143759], [0.28818313743592594 1.4749992558137126 0.19340492845344384; -1.5847339350618719 -0.5202837139329893 -1.0559788345007433])

julia> Z2 = Z × Z
CartesianProduct{Float64,Zonotope{Float64,Array{Float64,1},Array{Float64,2}},Zonotope{Float64,Array{Float64,1},Array{Float64,2}}}(Zonotope{Float64,Array{Float64,1},Array{Float64,2}}([-0.5237041738184478, -0.6673040904143759], [0.28818313743592594 1.4749992558137126 0.19340492845344384; -1.5847339350618719 -0.5202837139329893 -1.0559788345007433]), Zonotope{Float64,Array{Float64,1},Array{Float64,2}}([-0.5237041738184478, -0.6673040904143759], [0.28818313743592594 1.4749992558137126 0.19340492845344384; -1.5847339350618719 -0.5202837139329893 -1.0559788345007433]))

For the rest of the question, i assume the notation (z_1, z_2) refer to the concatenation of the 2-dimensional vectors z_1 and z_2. Let me try to reformulate (correct me if this is wrong) you'd like to find a set Y such that BZ := Z2 ∩ Y satisfies that all points z ∈ BZ are of the form z = (z_1, z_2) with z_1 + z_2 ∈ Z.

[yupbank]

yes.

[mforets]

Then if A, B are subsets of R² that satisfy that their Minkowski sum is Z, i.e. A ⊕ B = Z, one can take Y := A x B. However, finding such sets A and B does't seem an easy problem to me -- in fact, i think that there are many (infinite?) solutions.

Let A := 0 and B = Z, and let Y be the cartesian product Y = 0 x Z. Then BZ = Z2 ∩ Y is such that all elements of B2 are of the form (z_1, z_2) with z_1 + z_2 ∈ Z.

[yupbank]

hmm, let me do a 1-d zonotope illustration. for Z = [0, 1], then, Z2 = Z x Z would be a square.

▢

And the with the constraint applied (a, b) \in Z2, a+b \in [0, 1], the result BZ would be

◺

and sorry if I didn't explain clearly, I have some concrete 2d zonotope of interest, and I was trying to plot if possible the BZ using some projection in 3d. since BZ is in 4d.

[mforets]

To summarize a discussion on gitter, i think the following code gives an over-approximation of B2. In fact, if B2 is a (convex) polytope -- which i think it's the case..? -- this gives an exact representation so long as one chooses a sufficiently high number of support directions to evaluate. The idea is to overload the support function with that of B2 and use LazySets machinery to lazily represent the set, obtain a representation in constraint form, and project/plot the result.

using LazySets, JuMP, GLPK, Distributions, Plots
using LazySets: center, dim

struct Foo{N, VN, MN} <: LazySet{N}
    Z::Zonotope{N, VN, MN}
end
Foo(Z::AbstractZonotope) = Foo(convert(Zonotope, Z))

function LazySets.ρ(d::AbstractVector{N}, F::Foo{N}) where {N<:Real}

    Z = F.Z
    c = center(Z)
    G = genmat(Z)
    n, m = size(G)

    model = Model(GLPK.Optimizer)

    @variable(model, z1[1:n])
    @variable(model, z2[1:n])
    @variable(model, -1 <= ξ[1:m] <= 1)
    @variable(model, -1 <= α[1:m] <= 1)
    @variable(model, -1 <= β[1:m] <= 1)

    @constraint(model, z1 .+ z2 .- c .== G * ξ)
    @constraint(model, z1 .- c .== G * α)
    @constraint(model, z2 .- c .== G * β)
    @objective(model, Max, sum(d[i] * z1[i] for i in 1:n) + sum(d[n+i] * z2[i] for i in 1:n))
    optimize!(model)

    return objective_value(model)
end

Results

In one dimension:

using Plots

Z = Interval(0, 1)
F = Foo(Z)
BZ = overapproximate(F, OctDirections(2))

plot(Z × Z, lab="Z x Z")
plot!(BZ, lab="BZ")

In two dimensions:

Z = rand(Zonotope, dim=2);
F = Foo(Z)
BZ = overapproximate(F, OctDirections(4)) # this is 4-dimensional

plot(Projection(BZ, [1, 2]))