This notebook shows how to plot and get the nodes of a convex polyhedral function.
Consider the polyhedral function $$ f(x) = \max(-4x - 1, 2x - 1, -x/4, x/2) $$ in the interval $x \in [-1, 1]$.
As the function is convex (as it is the maximum of linear functions), its epigraph is convex. It is defined as the set of $x, y$ satisfying \begin{align*} -1 & \le x\\ x & \le 1\\ -4x - 1 & \le y\\ 2x - 1 & \le y\\ -x/4 & \le y\\ x/2 & \le y \end{align*} We can create this H-representation in Polyhedra as follows:
using Polyhedra
h = HalfSpace([-1, 0], 1) ∩ HalfSpace([1, 0], 1) ∩ HalfSpace([-4, -1], 1) ∩ HalfSpace([2, -1], 1) ∩
HalfSpace([-1/4, -1], 0) ∩ HalfSpace([1/2, -1], 0)
H-representation Polyhedra.Intersection{Float64,Array{Float64,1},Int64}: 6-element iterator of HalfSpace{Float64,Array{Float64,1}}: HalfSpace([-1.0, 0.0], 1.0) HalfSpace([1.0, 0.0], 1.0) HalfSpace([-4.0, -1.0], 1.0) HalfSpace([2.0, -1.0], 1.0) HalfSpace([-0.25, -1.0], 0.0) HalfSpace([0.5, -1.0], 0.0)
We now transform it to a polyhedron using CDDLib, the library that is chosen to compute the V-representation.
using CDDLib
p = polyhedron(h, CDDLib.Library())
Polyhedron CDDLib.Polyhedron{Float64}: 6-element iterator of HalfSpace{Float64,Array{Float64,1}}: HalfSpace([-1.0, 0.0], 1.0) HalfSpace([1.0, 0.0], 1.0) HalfSpace([-4.0, -1.0], 1.0) HalfSpace([2.0, -1.0], 1.0) HalfSpace([-0.25, -1.0], 0.0) HalfSpace([0.5, -1.0], 0.0)
Note that creating the polyhedron does not trigger the computation of the V-representation by itself.
This is done when the V-representation is queried, e.g. by vrep
.
We can see that 5 nodes of the polyhedral function with the ray $(0, 1)$ which is expected for an epigraph.
vrep(p)
V-representation CDDGeneratorMatrix{Float64,Float64}: 5-element iterator of Array{Float64,1}: [1.0, 1.0] [0.0, 0.0] [0.666667, 0.333333] [-0.266667, 0.0666667] [-1.0, 3.0], 1-element iterator of Ray{Float64,Array{Float64,1}}: Ray([-1.77636e-16, 1.0])
We see now that the V-representation is known by the polyhedron. Now both the H- and V-representation are available.
p
Polyhedron CDDLib.Polyhedron{Float64}: 6-element iterator of HalfSpace{Float64,Array{Float64,1}}: HalfSpace([-1.0, 0.0], 1.0) HalfSpace([1.0, 0.0], 1.0) HalfSpace([-4.0, -1.0], 1.0) HalfSpace([2.0, -1.0], 1.0) HalfSpace([-0.25, -1.0], 0.0) HalfSpace([0.5, -1.0], 0.0): 5-element iterator of Array{Float64,1}: [1.0, 1.0] [0.0, 0.0] [0.666667, 0.333333] [-0.266667, 0.0666667] [-1.0, 3.0], 1-element iterator of Ray{Float64,Array{Float64,1}}: Ray([-1.77636e-16, 1.0])
We can plot the epigraph as follows. Note that rays are not supported for 2D plotting, the polyhedron needs to be bounded so we add the halfspace $y \le 4$.
using Plots
plot(p ∩ HalfSpace([0, 1], 4))
scatter!([x[1] for x in points(p)], [x[2] for x in points(p)])
JuMP variables can be constrained to be in the epigraph.
Note that if the optimization model simply involves minimizing or maximizing a linear objective function over the polyhedron, the problem can simply the solved by iterating over the V-representation.
One can use Polyhedra.linear_objective_solver
to get VRepOptimizer
when the V-representation is computed.
This optimizer does exactly what we just described.
using JuMP
model = Model(Polyhedra.linear_objective_solver(p))
@variable(model, x)
@variable(model, y)
@constraint(model, [x, y] in p)
@objective(model, Min, x + y)
optimize!(model)
The optimal solution of this problem is one of the extreme points of p
:
termination_status(model)
OPTIMAL::TerminationStatusCode = 1
value(x), value(y)
(-0.26666666666666666, 0.06666666666666665)
plot!([-0.2, -1.2], [0.0, 1.0], linewidth=2) # Objective function
scatter!([value(x)], [value(y)], markershape=:star5, markersize=8) # Optimal solution
@objective(model, Max, x + y)
optimize!(model)
As we can see, maximizing $x + y$ is an unbounded problem, the infeasibility certificate is given by the extreme ray $(0, 1)$.
termination_status(model)
DUAL_INFEASIBLE::TerminationStatusCode = 3
value(x), value(y)
(-1.7763568394002506e-16, 1.0)
The epigraph can also be used in a more complex optimization problem involving other variables and constraints. Here Polyhedra.linear_objective_solver
should not be used because the feasible set of the new model is not exactly p
so the V-representation will have to be recomputed which is less efficient than solving the problem using a linear programming solver such as GLPK.
using JuMP
import GLPK
model = Model(GLPK.Optimizer)
@variable(model, x <= -2)
@variable(model, y)
@variable(model, 0 <= u <= 3)
@constraint(model, [2x + u, y] in p)
@objective(model, Min, u + y)
optimize!(model)
termination_status(model)
OPTIMAL::TerminationStatusCode = 1
value(x), value(y), value(u)
(-2.0, 3.0, 3.0)
Consider the polyhedral function $$ f(x) = \max(-4x - 2y - 1, 2x - y - 1, -x/4 + y/2, x/2 + y) $$ in the interval $(x, y) \in [-1, 1]^2$.
As the function is convex (as it is the maximum of linear functions), its epigraph is convex. It is defined as the set of $x, y, z$ satisfying \begin{align*} -1 & \le x\\ x & \le 1\\ -1 & \le y\\ y & \le 1\\ -4x - 2y - 1 & \le z\\ 2x - y - 1 & \le z\\ -x/4 + y/2 & \le z\\ x/2 + y & \le z \end{align*} We can create this H-representation in Polyhedra as follows:
using Polyhedra
h = HalfSpace([-1, 0, 0], 1) ∩ HalfSpace([1, 0, 0], 1) ∩
HalfSpace([0, -1, 0], 1) ∩ HalfSpace([0, 1, 0], 1) ∩
HalfSpace([-4, -2, -1], 1) ∩ HalfSpace([2, -1, -1], 1) ∩
HalfSpace([-1/4, 1/2, -1], 0) ∩ HalfSpace([1/2, 1, -1], 0)
H-representation Polyhedra.Intersection{Float64,Array{Float64,1},Int64}: 8-element iterator of HalfSpace{Float64,Array{Float64,1}}: HalfSpace([-1.0, 0.0, 0.0], 1.0) HalfSpace([1.0, 0.0, 0.0], 1.0) HalfSpace([0.0, -1.0, 0.0], 1.0) HalfSpace([0.0, 1.0, 0.0], 1.0) HalfSpace([-4.0, -2.0, -1.0], 1.0) HalfSpace([2.0, -1.0, -1.0], 1.0) HalfSpace([-0.25, 0.5, -1.0], 0.0) HalfSpace([0.5, 1.0, -1.0], 0.0)
We create the polyhedron the same way as for 1-dimensional polyhedral functions.
using CDDLib
p = polyhedron(h, CDDLib.Library())
Polyhedron CDDLib.Polyhedron{Float64}: 8-element iterator of HalfSpace{Float64,Array{Float64,1}}: HalfSpace([-1.0, 0.0, 0.0], 1.0) HalfSpace([1.0, 0.0, 0.0], 1.0) HalfSpace([0.0, -1.0, 0.0], 1.0) HalfSpace([0.0, 1.0, 0.0], 1.0) HalfSpace([-4.0, -2.0, -1.0], 1.0) HalfSpace([2.0, -1.0, -1.0], 1.0) HalfSpace([-0.25, 0.5, -1.0], 0.0) HalfSpace([0.5, 1.0, -1.0], 0.0)
We obtain the 10 nodes as follows along with the $(0, 0, 1)$ ray that is expected as the polyhedron is an epigraph.
vrep(p)
V-representation CDDGeneratorMatrix{Float64,Float64}: 10-element iterator of Array{Float64,1}: [0.222222, -0.333333, -0.222222] [1.0, 0.25, 0.75] [1.0, 1.0, 1.5] [-0.666667, 1.0, 0.666667] [1.0, -1.0, 2.0] [0.0888889, -0.533333, -0.288889] [0.166667, -1.0, 0.333333] [-0.933333, 1.0, 0.733333] [-1.0, 1.0, 1.0] [-1.0, -1.0, 5.0], 1-element iterator of Ray{Float64,Array{Float64,1}}: Ray([-6.66134e-16, -7.40149e-17, 1.0])
m = Polyhedra.Mesh(p)
using MeshCat
vis = Visualizer()
setobject!(vis, m)
IJuliaCell(vis)
using Makie
Makie.mesh(m, color=:blue)
using Makie
Makie.wireframe(m)