Polyhedron of JuMP model with Polyhedra or Lazyset

Specific Domains Optimization (Mathematical)
, , , ,
1

Hello,
I want to plot a 3D polyhedron of a JuMP model for education purposes. I have seen these packages Polyhedra and LazySets, but I am unsure how to convert a JuMP model. I believe it requires projection/elimination. A MWE of a JuMP model is given below. Any ideas?

using JuMP, Cbc
function j_poly()
    q = [109,115,102,114,108]
    e = [40 150 2.5 50 40; 40 50 2.5 50 40; 40 50 2.5 50 40; 40 50 2.5 50 40; 40 50 102.5 50 40;40 50 0 105 104;40 105 2 2 2; 104 2 2 2 3]
    l = [16, 18, 10, 13, 31, 16, 18, 10]
    p = [100, 105, 90, 105, 95]
    f  = [0.01, 0.012, 0.016, 0.019, 0.022, 0.027, 0.03, 0.045]

    m = Model(Cbc.Optimizer)
    @variable(m,0<=x[1:5]<=1)
    @constraint(m,c1,sum(x[i]*q[i] for i in 1:5) == 100)
    @constraint(m,c2[k=1:8],sum((l[j] - sum(e[j,i]*x[i] for i in 1:5)*(1+f[j])^(-j)) for j in 1:k) <=30*(1+f[k])^(-k))
    @objective(m,Min,sum(p[i]*x[i] for i in 1:5))

    return m, x
end
2

Polyhedra has what you want:

using JuMP
import CDDLib
import Makie
import GLMakie
import Polyhedra

function j_poly()
    q = [109,115,102,114,108]
    e = [40 150 2.5 50 40; 40 50 2.5 50 40; 40 50 2.5 50 40; 40 50 2.5 50 40; 40 50 102.5 50 40;40 50 0 105 104;40 105 2 2 2; 104 2 2 2 3]
    l = [16, 18, 10, 13, 31, 16, 18, 10]
    p = [100, 105, 90, 105, 95]
    f  = [0.01, 0.012, 0.016, 0.019, 0.022, 0.027, 0.03, 0.045]
    m = Model()
    @variable(m,0<=x[1:5]<=1)
    @constraint(m,c1,sum(x[i]*q[i] for i in 1:5) == 100)
    @constraint(m,c2[k=1:8],sum((l[j] - sum(e[j,i]*x[i] for i in 1:5)*(1+f[j])^(-j)) for j in 1:k) <=30*(1+f[k])^(-k))
    @objective(m,Min,sum(p[i]*x[i] for i in 1:5))
    return m, x
end

model, _ = j_poly()

poly = Polyhedra.polyhedron(model, CDDLib.Library())
poly = Polyhedra.eliminate(poly, [4, 5])

mesh = Polyhedra.Mesh(poly)
Makie.mesh(mesh; color=:blue)

Docs:

1 Like
3

Thank you
I don’t fully understand the following line:

And the documentation is not very clear. What does [4,5] represent here? I can see that the variable x is defined for range [1:5] and c2 generates 8 constraints, so perhaps [4,5] represents half of the number of constraints and total number of variables.
If we take a more generic example and increase number of variables and constrain c2 dimension:

@variable(m,0<=x[1:500]<=1)
@constraint(m,c2[k=1:100],sum((l[j] - sum(e[j,i]*x[i] for i in 1:500)*(1+f[j])^(-j)) for j in 1:k) <=30*(1+f[k])^(-k))

How should this line poly = Polyhedra.eliminate(poly, [4, 5]) be changed?

4

The [4, 5] are the dimensions of the polyhedron to eliminate by projecting them out. Your example has 5 variables, but we can only plot 3-dimensional objects, so I needed to get rid of two dimensions. I choose, arbitrarily, the last two dimensions.

@variable(m,0<=x[1:500]<=1)
How should this line poly = Polyhedra.eliminate(poly, [4, 5]) be changed?

You would need to choose 497 variables to eliminate.

5

That makes sense. I believe if a model has multiple variables say x[1:500] and y[1:10] then it will be difficult to know the order in which a JuMP model contains them, so determining which ones to eliminate could be challenging.

The elimination operation seems quite expensive on a problem of this size, unless I am not doing it correctly.

6

Yes. But I think the underlying problem is that it doesn’t make sense to plot a 510-dimensional space. For educational purposes, stick to two or three variables.

7

I agree. I wanted to draw a shape like this: Polyhedron -- from Wolfram MathWorld or anything spherical! But I am not sure what type of linear program would produce a similar shape.

8

A simple example would be tangent planes on a sphere:

using JuMP
import CDDLib
import Makie
import GLMakie
import Polyhedra

model = Model()
@variable(model, x[1:3])
for _ in 1:200
    θ, ϕ = 2π * rand(), 2π * rand()
    a = [cos(θ) * sin(ϕ), sin(θ) * sin(ϕ), cos(ϕ)]
    @constraint(model, a' * x <= 1)
end
poly = Polyhedra.polyhedron(model, CDDLib.Library())
mesh = Polyhedra.Mesh(poly)
Makie.mesh(mesh; color=:lightblue)

image
image812×662 69.1 KB

2 Likes