# Polyhedron of JuMP model with Polyhedra or Lazyset

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

``````

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

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?

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.

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.

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.

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.

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)
``````

2 Likes