Optimization

A polyhedron can represents the feasible set of an optimization program. The program is infeasible when the polyhedron is empty.

isempty(p::Rep, solver=Polyhedra.linear_objective_solver(p))

If the V-representation of the polyhedron has been computed, it can be used to solve the linear program.

VRepOptimizer{T} <: AbstractPolyhedraOptimizer{T}

Otherwise, any programming solver implementing the MathOptInterface interface can be used. See here for a list of available solvers.

default_solver(p::Rep)

linear_objective_solver(p::Rep, solver=default_solver(p))

model = Model(solver)
x = @variable(model, [1:fulldim(p)])
@constraint(model, x in p)
@objective(model, c ⋅ x)

Using a polyhedron for in an optimization model

A polyhedron or representation can be used in the constraint of a JuMP model. For instance, consider the 1-simplex:

julia> using Polyhedra

julia> simplex = HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0) ∩ HyperPlane([1, 1], 1)
H-representation Polyhedra.Intersection{Int64, Vector{Int64}, Int64}:
1-element iterator of HyperPlane{Int64, Vector{Int64}}:
 HyperPlane([1, 1], 1),
2-element iterator of HalfSpace{Int64, Vector{Int64}}:
 HalfSpace([-1, 0], 0)
 HalfSpace([0, -1], 0)

and the following JuMP model with two variables

julia> using JuMP

julia> model = Model()
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> @variable(model, λ[1:2])
2-element Vector{VariableRef}:
 λ[1]
 λ[2]

The variables can be constrained to belong to the simplex as follows:

julia> @constraint(model, λ in simplex)
[λ[1], λ[2]] ∈ Polyhedra.PolyhedraOptSet{Int64, Polyhedra.Intersection{Int64, Vector{Int64}, Int64}}(HyperPlane([1, 1], 1) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0))

but a vector of affine or quadratic expressions can also be constrained to belong to the simplex:

julia> A = [1  1
            1 -1]
2×2 Matrix{Int64}:
 1   1
 1  -1

julia> @constraint(model, A * λ in simplex)
[λ[1] + λ[2], λ[1] - λ[2]] ∈ Polyhedra.PolyhedraOptSet{Int64, Polyhedra.Intersection{Int64, Vector{Int64}, Int64}}(HyperPlane([1, 1], 1) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0))

We can verify that the model contains both constraints:

julia> model
A JuMP Model
Feasibility problem with:
Variables: 2
`Array{JuMP.VariableRef,1}`-in-`Polyhedra.PolyhedraOptSet{Int64,Polyhedra.Intersection{Int64,Array{Int64,1},Int64}}`: 1 constraint
`Array{JuMP.GenericAffExpr{Float64,JuMP.VariableRef},1}`-in-`Polyhedra.PolyhedraOptSet{Int64,Polyhedra.Intersection{Int64,Array{Int64,1},Int64}}`: 1 constraint
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
Names registered in the model: λ

When the model is solved, the constraint is automatically transformed into appropriate constraints if the optimizer does not support consraints with the set Polyhedra.PolyhedraOptSet.

julia> import GLPK

julia> set_optimizer(model, GLPK.Optimizer)

julia> optimize!(model)

julia> termination_status(model)
OPTIMAL::TerminationStatusCode = 1

julia> value.(λ)
2-element Array{Float64,1}:
 0.5
 0.5

For instance, GLPK, does not support Polyhedra.PolyhedraOptSet constraints but supports MOI.EqualTo{Float64} and MOI.LessThan{Float64}. The polyhedral constraints are therefore bridged into several MOI.EqualTo{Float64} and MOI.LessThan{Float64} constraints using the following constraint bridge:

PolyhedraToLPBridge{T}

See Polyhedral Function for an example notebook.

Creating a polyhedron from the feasible set of a JuMP model

A typical application of polyhedral computation is the computation of the set of extreme points and rays of the feasible set of an optimization problem. This comes from the fact that given a minimization of a concave function (or maximization of a convex function) on a convex feasible set (e.g. Linear Programming), we are either in the following three situations:

• The feasible set is empty, i.e. the problem is infeasible.
• An extreme ray is optimal, i.e. the problem is unbounded (or it may also be bounded if the objective is constant along the ray).
• An extreme point is optimal.

A JuMP model is treated by polyhedron just like any H-representation. For example, the hypercube of dimension n can be created as follows

using JuMP, Polyhedra
model = Model()
@variable(model, 0 ≤ x[1:2] ≤ 1)
h = hrep(model)

H-representation LPHRep{Float64}:
4-element iterator of HalfSpace{Float64, SparseArrays.SparseVector{Float64, Int64}}:
 HalfSpace(  [1]  =  -1.0, -0.0)
 HalfSpace(  [2]  =  -1.0, -0.0)
 HalfSpace(  [1]  =  1.0, 1.0)
 HalfSpace(  [2]  =  1.0, 1.0)

The name of the variables for each dimension can be recovered as follows

dimension_names(h)

2-element Vector{String}:
 "x[1]"
 "x[2]"

Note that the names of the JuMP variables are lost in the conversion to a polyhedron that does not support names, e.g.,

poly = polyhedron(model, CDDLib.Library(:exact))

However, the ordering of the dimension of the polyhedron is guaranteed to correspond to the order of the JuMP variables as listed by all_variables:

all_variables(model)

2-element Vector{JuMP.VariableRef}:
 x[1]
 x[2]

So all_variables(model)[i] provides the JuMP variable corresponding to the ith dimension. The reverse mapping can be constructed as follows:

Dict(v => i for (i, v) in enumerate(all_variables(model)))

Dict{JuMP.VariableRef, Int64} with 2 entries:
  x[1] => 1
  x[2] => 2