# JuMP model Polyhedron

**URL:** <https://discourse.julialang.org/t/jump-model-polyhedron/50108>\
**Category:** Optimization (Mathematical)\
**Tags:** plotting, jump, polyhedra\
**Created:** [November 13, 2020, 6:53pm UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108 "2020-11-13T18:53:08Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![JohnZ](https://avatars.discourse-cdn.com/v4/letter/j/ba8739/32.png) [@JohnZ](https://discourse.julialang.org/u/JohnZ)\
**Post date:** [November 13, 2020, 6:53pm UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/1 "2020-11-13T18:53:08Z")

</div>

Hi, I am trying to plot the feasible set of a JuMP model as a 3D figure. I have used the syntax provided in Polyhedra.jl documentation to create a polyhedron: [Optimization · Polyhedra](https://juliapolyhedra.github.io/Polyhedra.jl/stable/optimization/)

A simple example is given below:

```julia
   using JuMP, Clp, Polyhedra, CDDLib, Plots, Makie
   Data1 =rand(1.0:20,12)     
   Data2 = rand(1.0:15,12)
    m = Model(Clp.Optimizer)
    @variable(m,0<=x[1:12]<=1)
    @variable(m,0<=y[1:12]<=1)
    @objective(m,Max,sum(Data1[i]*x[i] for i in 1:12))
    @constraint(m,con2[i=1:12], sum(Data1[i]*x[i] -Data2[i]*y[i] for i in 1:12) <=0.3)
    #optimize!(m)
    #JuMP.value.(x)
    poly = polyhedron(m, CDDLib.Library(:exact))

```

I am not sure how to plot it as `poly` represents several half spaces and I do not know what I need to pass on to the `Plots` or `Makie` package. For example, the code below would not work:

```julia
 Makie.mesh(poly, color=:blue)

```

Any ideas?

---

<div class="post-metadata">

**Author:** ![blegat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blegat/32/217090_2.png) [@blegat](https://discourse.julialang.org/u/blegat)\
**Post date:** [November 14, 2020, 6:56pm UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/2 "2020-11-14T18:56:14Z")

</div>

You need to wrap it in a `Polyhedra.Mesh` as it is a subtype of the abstract type that the 3D libraries expect, so do

```julia
Makie.mesh(Polyhedra.Mesh(poly), color=:blue)

```

See [the doc](https://juliapolyhedra.github.io/Polyhedra.jl/stable/plot/#D-plotting-with-Plots-2) for more details.

---

<div class="post-metadata">

**Author:** ![JohnZ](https://avatars.discourse-cdn.com/v4/letter/j/ba8739/32.png) [@JohnZ](https://discourse.julialang.org/u/JohnZ)\
**Post date:** [November 14, 2020, 10:42pm UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/3 "2020-11-14T22:42:08Z")

</div>

> [@blegat](#):
>
> `Makie.mesh(Polyhedra.Mesh(poly), color=:blue)`

This does not work for me. I get this error:

```julia
ERROR: MethodError: no method matching fulldecompose(::Polyhedra.Mesh{24,Rational{BigInt},CDDLib.Polyhedron{Rational{BigInt}}}, ::Type{Rational{BigInt}})
Closest candidates are:
  fulldecompose(::Polyhedra.Mesh{3,T,PT} where PT<:Polyhedron{T} where T, ::Type{T}) where 
T at C:\Users\A\.julia\packages\Polyhedra\FADlK\src\decompose.jl:79
  fulldecompose(::Polyhedra.Mesh{N,T,PT} where PT<:Polyhedron{T}) where {N, T} at C:\Users\A\.julia\packages\Polyhedra\FADlK\src\decompose.jl:203
Stacktrace:
 [1] fulldecompose(::Polyhedra.Mesh{24,Rational{BigInt},CDDLib.Polyhedron{Rational{BigInt}}}) at C:\Users\A\.julia\packages\Polyhedra\FADlK\src\decompose.jl:203
 [2] fulldecompose! at C:\Users\A\.julia\packages\Polyhedra\FADlK\src\decompose.jl:38 [inlined]
 [3] coordinates at C:\Users\A\.julia\packages\Polyhedra\FADlK\src\decompose.jl:205 [inlined]
 [4] decompose(::Type{Point{24,Float32}}, ::Polyhedra.Mesh{24,Rational{BigInt},CDDLib.Polyhedron{Rational{BigInt}}}) at C:\Users\A\.julia\packages\GeometryBasics\csguK\src\interfaces.jl:109
 [5] mesh(::Polyhedra.Mesh{24,Rational{BigInt},CDDLib.Polyhedron{Rational{BigInt}}}; pointtype::Type{T} where T, facetype::Type{T} where T, uv::Type{T} where T, normaltype::Type{T} where T) at C:\Users\A\.julia\packages\GeometryBasics\csguK\src\meshes.jl:104
 [6] uv_normal_mesh at C:\Users\A\.julia\packages\GeometryBasics\csguK\src\meshes.jl:177 [inlined]
 [7] convert_arguments at C:\Users\A\.julia\packages\AbstractPlotting\JCbJs\src\conversions.jl:480 [inlined]
 [8] plot!(::Scene, ::Type{Combined{AbstractPlotting.mesh,ArgType} where ArgType}, ::Attributes, ::Polyhedra.Mesh{24,Rational{BigInt},CDDLib.Polyhedron{Rational{BigInt}}}; kw_attributes::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\A\.julia\packages\AbstractPlotting\JCbJs\src\interfaces.jl:610
 [9] plot! at C:\Users\A\.julia\packages\AbstractPlotting\JCbJs\src\interfaces.jl:597 [inlined]
 [10] plot(::Type{Combined{AbstractPlotting.mesh,ArgType} where ArgType}, ::Attributes, ::Polyhedra.Mesh{24,Rational{BigInt},CDDLib.Polyhedron{Rational{BigInt}}}; kw_attributes::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\A\.julia\packages\AbstractPlotting\JCbJs\src\interfaces.jl:558
 [11] plot at C:\Users\A\.julia\packages\AbstractPlotting\JCbJs\src\interfaces.jl:555 [inlined]
 [12] #plot#192 at C:\Users\A\.julia\packages\AbstractPlotting\JCbJs\src\interfaces.jl:549 
[inlined]
 [13] #mesh#148 at C:\Users\A\.julia\packages\AbstractPlotting\JCbJs\src\recipes.jl:12 [inlined]
 [14] top-level scope at none:1

```

---

<div class="post-metadata">

**Author:** ![blegat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blegat/32/217090_2.png) [@blegat](https://discourse.julialang.org/u/blegat)\
**Post date:** [November 14, 2020, 11:56pm UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/4 "2020-11-14T23:56:40Z")

</div>

It seems your polyhedron has dimension 24, you can only plot 3D polyhedra. Maybe you need to project it first; see [Projection/Elimination · Polyhedra](https://juliapolyhedra.github.io/Polyhedra.jl/stable/projection/).

---

<div class="post-metadata">

**Author:** ![JohnZ](https://avatars.discourse-cdn.com/v4/letter/j/ba8739/32.png) [@JohnZ](https://discourse.julialang.org/u/JohnZ)\
**Post date:** [November 15, 2020, 12:07pm UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/5 "2020-11-15T12:07:35Z")

</div>

Thanks. I have looked at the documentation for projection and I think I need to use  
`project(p::Polyhedron, pset, algo)` or `fixandeliminate(p::HRep{T}, I, v)`

Are there any examples of projecting to a 3D space? I am not sure what input arguments are required for the above functions, particularly for `pset` in the first function, and `I `and `v` in the second.

---

<div class="post-metadata">

**Author:** ![blegat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blegat/32/217090_2.png) [@blegat](https://discourse.julialang.org/u/blegat)\
**Post date:** [November 15, 2020, 1:08pm UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/6 "2020-11-15T13:08:40Z")

</div>

There is an example [here](https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/master/examples/3D%20Plotting%20a%20projection%20of%20the%204D%20permutahedron.ipynb) which projects on an orthogonal 3D basis which is not aligned with the original axis.

To make it clearer, you have a 24-dimensional H-represented polyhedron. You would like to have a 3-dimensional view of it. You have 3 operations you can do

1. Fix some dimension to a fixed value to reduce the dimension with `fixandeliminate`, this is **very cheap** to do.
2. Project the polyhedron along some of the 24 axis. I would recommend using CDDLib with `CDDLib.BlockElimination()` as algorithm, this is **expensive** to do.
3. Project it along combinations of the 24 axis as in the example mentioned above, this requires computing the V-representaiton of the 24 axis first, this is **very expensive** to do.

You can of course combine different operations, I would recommend starting with 1) to cheaply reduce the dimensionality as lower dimensional operations are usually cheaper.

The `pset` and `I` are vector if indices from 1 to 24. For `project`, they are the dimension you project to and for `fixandeliminate`, they are the dimension that are fixed to the corresponding values in `v` and then eliminated.

To see which dimension corresponds to which JuMP variable, use `dimension_names` as shown in [this example](https://juliapolyhedra.github.io/Polyhedra.jl/stable/generated/Extended%20Formulation/).  
Of course, this suppose you set names to your JuMP variables.

---

<div class="post-metadata">

**Author:** ![JohnZ](https://avatars.discourse-cdn.com/v4/letter/j/ba8739/32.png) [@JohnZ](https://discourse.julialang.org/u/JohnZ)\
**Post date:** [November 15, 2020, 5:36pm UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/7 "2020-11-15T17:36:07Z")

</div>

Thank you for the detailed explanation. I will start with 1) as I was looking to plot the polyhedron of problem of much higher dimension than the MWE.

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [November 16, 2020, 3:09pm UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/8 "2020-11-16T15:09:21Z")

</div>

1. View \hat{X} \supseteq X obtained by taking the lower-dim projection lazily (and exactly) and overapproximating with support functions to the desired accuracy. This approach is both **efficient** and **scalable** , but is somehow more advanced as it requires combining some operations.

---

<div class="post-metadata">

**Author:** ![JohnZ](https://avatars.discourse-cdn.com/v4/letter/j/ba8739/32.png) [@JohnZ](https://discourse.julialang.org/u/JohnZ)\
**Post date:** [November 16, 2020, 6:48pm UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/9 "2020-11-16T18:48:39Z")

</div>

> [@mforets](#):
>
> View ^X⊇X\hat{X} \supseteq X obtained by taking the lower-dim projection lazily (and exactly) and overapproximating with support functions to the desired accuracy.

Sorry, I do not understand what does this mean.

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [November 16, 2020, 8:16pm UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/10 "2020-11-16T20:16:31Z")

</div>

@JohnZ, to extend on my comment:

There is a classic method for approximating polyhedra [1], useful for applications that require operating with high-dimensional subsets of \mathbb{R}^n [1] (“high” meaning n \> 10, i.e. when options 2 and 3 above are not a viable alternative). This is the kind of approach taken in our library [LazySets.jl](https://github.com/JuliaReach/LazySets.jl).

[1] Kamenev, G. K. (1996). An algorithm for approximating polyhedra. Computational Mathematics and Mathematical Physics, 36(4), 533–544,

Let me comment in passing that `LazySets.jl` does heavily rely on `Polyhedra.jl`: whenever your computation requires to work with H-rep and V-rep interchangeably, quantifier elimination and so on. The tools interface with each other nicely. On the other hand, there are many other things that can be done “on-the-fly” very efficiently! This is where the mathematical methods implemented in `LazySets.jl` shine.

* * *

The following code implements the example from the OP. Note that `m` had 12 times the same constraint, so I only kept one as the half-space `C`. I also changed the sign in front of `Data2`; otherwise plots were too boring (the result was always a square). This result shows that the feasible set is certainly contained in the triangle. Moreover, this triangle is tight by construction. But we can say more than that: the result is actually exact, because the underapproximation, coincides with the overapproximation for eps small enough (you can also specify the directions in the overapproximation too, but Kamenev’s method picks them automatically).

Overapproximation in 3D is also possible, though it seems like a good idea that you first get familiar with the 2D code. See the documentation of LazySets.jl for details or ask back if you need more help!

About runtime, `X̂(0.01)` takes \< 1ms and the call to `underapproximate` is 3 times as fast. This approach scales well with the dimension, as it is internally just solving LPs where appropriate.

```julia
using LazySets, Plots

Data1 = rand(1.0:20,12)     
Data2 = rand(1.0:15,12)

B = Hyperrectangle(low=zeros(24), high=ones(24))
C = HalfSpace(vcat(Data1, Data2), 0.3)
P = intersection(B, C)

X = Projection(P, [23, 24]) # lazy projection
ρ(vcat(Data1, zeros(12)), P) # compute optimum value max sum(Data1[i]*x[i])
X̂(ε) = overapproximate(X, HPolygon, ε); # epsilon-close approximation

Xunder = plot(underapproximate(X, BoxDirections(2)))

pover = plot()
plot!(pover, X̂(1.0), alpha=.5, title="Overapproximation", lab="ε = 1.0", lw=3.0)
plot!(pover, X̂(0.01), alpha=.5, lab="ε = 0.01")

punder = plot(Xunder, alpha=0.1, title="Underapproximation", lab="boxdirs")

plot(pover, punder, layout=(1, 2))

```

 ![Screenshot from 2020-11-16 17-08-31](https://global.discourse-cdn.com/julialang/original/3X/0/e/0e860125f0c1bd9b57d2fad892ac83b3c7589bd7.png)

---

<div class="post-metadata">

**Author:** ![JohnZ](https://avatars.discourse-cdn.com/v4/letter/j/ba8739/32.png) [@JohnZ](https://discourse.julialang.org/u/JohnZ)\
**Post date:** [November 18, 2020, 12:33am UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/11 "2020-11-18T00:33:00Z")

</div>

Thanks for the detailed example. I was after generic code that I can simply use to do a 3D plot for any LP defined in JuMP. But it seems like, that is not possible as the projection is problem dependent.

> [@mforets](#):
>
> ```julia
> X = Projection(P, [23, 24]) # lazy projection
> 
> ```

What does 23 and 24 represent in this line of code? I believe 24 corresponds to the total number of variables x and y bounds, but I am not sure about 23.

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [November 18, 2020, 1:33am UTC](https://discourse.julialang.org/t/jump-model-polyhedron/50108/12 "2020-11-18T01:33:34Z")

</div>

> [@JohnZ](#):
>
> What does 23 and 24 represent in this line of code?

`y[11]` and `y[12]` respectively.

> I was after generic code that I can simply use to do a 3D plot for any LP defined in JuMP.

Yes, I saw. My comment gives a quick way to project in 2D, which is usually fine, since you can pick different combinations of variables. For 3D, if the exact projection methods are too expensive (the LazySets command to take concrete projections is `project`), you can try overapproximating with octagonal or spherical directions.