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

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.