JuMP model Polyhedron

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: https://juliapolyhedra.github.io/Polyhedra.jl/stable/optimization/

A simple example is given below:

   using JuMP, Clp, Polyhedra, CDDLib, Plots, Makie
   Data1 =rand(1.0:20,12)     
   Data2 = rand(1.0:15,12)
    m = Model(Clp.Optimizer)
    @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)
    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:

 Makie.mesh(poly, color=:blue)

Any ideas?

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

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

See the doc for more details.

1 Like

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

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
 [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 
 [13] #mesh#148 at C:\Users\A\.julia\packages\AbstractPlotting\JCbJs\src\recipes.jl:12 [inlined]
 [14] top-level scope at none:1

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

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.

There is an example here 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.
Of course, this suppose you set names to your JuMP variables.


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.

  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.

Sorry, I do not understand what does this mean.

@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.

[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.

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


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.

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.

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.