# Polyhedron of JuMP model with Polyhedra or Lazyset

**URL:** <https://discourse.julialang.org/t/polyhedron-of-jump-model-with-polyhedra-or-lazyset/89218>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, optimization, geometry, polygons, polyhedra\
**Created:** [October 24, 2022, 10:42pm UTC](https://discourse.julialang.org/t/polyhedron-of-jump-model-with-polyhedra-or-lazyset/89218 "2022-10-24T22:42:53Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![papar](https://avatars.discourse-cdn.com/v4/letter/p/c5a1d2/32.png) [@papar](https://discourse.julialang.org/u/papar)\
**Post date:** [October 24, 2022, 10:42pm UTC](https://discourse.julialang.org/t/polyhedron-of-jump-model-with-polyhedra-or-lazyset/89218/1 "2022-10-24T22:42:53Z")

</div>

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?

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

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 25, 2022, 12:46am UTC](https://discourse.julialang.org/t/polyhedron-of-jump-model-with-polyhedra-or-lazyset/89218/2 "2022-10-25T00:46:51Z")

</div>

Polyhedra has what you want:

```julia
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:

- [Optimization · Polyhedra](https://juliapolyhedra.github.io/Polyhedra.jl/stable/optimization/#Creating-a-polyhedron-from-the-feasible-set-of-a-JuMP-model)
- [Projection/Elimination · Polyhedra](https://juliapolyhedra.github.io/Polyhedra.jl/stable/projection/)
- [Plot · Polyhedra](https://juliapolyhedra.github.io/Polyhedra.jl/stable/plot/)

---

<div class="post-metadata">

**Author:** ![papar](https://avatars.discourse-cdn.com/v4/letter/p/c5a1d2/32.png) [@papar](https://discourse.julialang.org/u/papar)\
**Post date:** [October 25, 2022, 6:31am UTC](https://discourse.julialang.org/t/polyhedron-of-jump-model-with-polyhedra-or-lazyset/89218/3 "2022-10-25T06:31:02Z")

</div>

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

> [@odow](#):
>
> `poly = Polyhedra.eliminate(poly, [4, 5])`

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:

```julia
@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?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 25, 2022, 8:05pm UTC](https://discourse.julialang.org/t/polyhedron-of-jump-model-with-polyhedra-or-lazyset/89218/4 "2022-10-25T20:05:51Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![papar](https://avatars.discourse-cdn.com/v4/letter/p/c5a1d2/32.png) [@papar](https://discourse.julialang.org/u/papar)\
**Post date:** [October 26, 2022, 9:35pm UTC](https://discourse.julialang.org/t/polyhedron-of-jump-model-with-polyhedra-or-lazyset/89218/5 "2022-10-26T21:35:47Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 26, 2022, 10:06pm UTC](https://discourse.julialang.org/t/polyhedron-of-jump-model-with-polyhedra-or-lazyset/89218/6 "2022-10-26T22:06:39Z")

</div>

> [@papar](#):
>
> I believe if a model has multiple variables say x[1:500] and y[1:10] then it will be

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.

---

<div class="post-metadata">

**Author:** ![papar](https://avatars.discourse-cdn.com/v4/letter/p/c5a1d2/32.png) [@papar](https://discourse.julialang.org/u/papar)\
**Post date:** [October 26, 2022, 10:12pm UTC](https://discourse.julialang.org/t/polyhedron-of-jump-model-with-polyhedra-or-lazyset/89218/7 "2022-10-26T22:12:03Z")

</div>

I agree. I wanted to draw a shape like this: [Polyhedron -- from Wolfram MathWorld](https://mathworld.wolfram.com/Polyhedron.html) or anything spherical! But I am not sure what type of linear program would produce a similar shape.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 26, 2022, 10:53pm UTC](https://discourse.julialang.org/t/polyhedron-of-jump-model-with-polyhedra-or-lazyset/89218/8 "2022-10-26T22:53:48Z")

</div>

A simple example would be tangent planes on a sphere:

```plaintext
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](https://global.discourse-cdn.com/julialang/original/3X/0/1/0127122e8c4dbb2268fe6bdaebe462891439d425.png)