{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook shows how to plot and get the nodes of a convex polyhedral function.\n",
"\n",
"# One dimensional polyhedral function\n",
"\n",
"Consider the polyhedral function\n",
"$$ f(x) = \\max(-4x - 1, 2x - 1, -x/4, x/2) $$\n",
"in the interval $x \\in [-1, 1]$.\n",
"\n",
"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\n",
"\\begin{align*}\n",
" -1 & \\le x\\\\\n",
" x & \\le 1\\\\\n",
" -4x - 1 & \\le y\\\\\n",
" 2x - 1 & \\le y\\\\\n",
" -x/4 & \\le y\\\\\n",
" x/2 & \\le y\n",
"\\end{align*}\n",
"We can create this H-representation in Polyhedra as follows:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"H-representation Polyhedra.Intersection{Float64,Array{Float64,1},Int64}:\n",
"6-element iterator of HalfSpace{Float64,Array{Float64,1}}:\n",
" HalfSpace([-1.0, 0.0], 1.0)\n",
" HalfSpace([1.0, 0.0], 1.0)\n",
" HalfSpace([-4.0, -1.0], 1.0)\n",
" HalfSpace([2.0, -1.0], 1.0)\n",
" HalfSpace([-0.25, -1.0], 0.0)\n",
" HalfSpace([0.5, -1.0], 0.0)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Polyhedra\n",
"h = HalfSpace([-1, 0], 1) ∩ HalfSpace([1, 0], 1) ∩ HalfSpace([-4, -1], 1) ∩ HalfSpace([2, -1], 1) ∩\n",
" HalfSpace([-1/4, -1], 0) ∩ HalfSpace([1/2, -1], 0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now transform it to a polyhedron using [CDDLib](https://github.com/JuliaPolyhedra/CDDLib.jl), the library that is chosen to compute the V-representation."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Polyhedron CDDLib.Polyhedron{Float64}:\n",
"6-element iterator of HalfSpace{Float64,Array{Float64,1}}:\n",
" HalfSpace([-1.0, 0.0], 1.0)\n",
" HalfSpace([1.0, 0.0], 1.0)\n",
" HalfSpace([-4.0, -1.0], 1.0)\n",
" HalfSpace([2.0, -1.0], 1.0)\n",
" HalfSpace([-0.25, -1.0], 0.0)\n",
" HalfSpace([0.5, -1.0], 0.0)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using CDDLib\n",
"p = polyhedron(h, CDDLib.Library())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Getting the extreme points\n",
"\n",
"Note that creating the polyhedron does not trigger the computation of the V-representation by itself.\n",
"This is done when the V-representation is queried, e.g. by `vrep`.\n",
"We can see that 5 nodes of the polyhedral function with the ray $(0, 1)$ which is expected for an epigraph."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"V-representation CDDGeneratorMatrix{Float64,Float64}:\n",
"5-element iterator of Array{Float64,1}:\n",
" [1.0, 1.0]\n",
" [0.0, 0.0]\n",
" [0.666667, 0.333333]\n",
" [-0.266667, 0.0666667]\n",
" [-1.0, 3.0],\n",
"1-element iterator of Ray{Float64,Array{Float64,1}}:\n",
" Ray([-1.77636e-16, 1.0])"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vrep(p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see now that the V-representation is known by the polyhedron. Now both the H- and V-representation are available."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Polyhedron CDDLib.Polyhedron{Float64}:\n",
"6-element iterator of HalfSpace{Float64,Array{Float64,1}}:\n",
" HalfSpace([-1.0, 0.0], 1.0)\n",
" HalfSpace([1.0, 0.0], 1.0)\n",
" HalfSpace([-4.0, -1.0], 1.0)\n",
" HalfSpace([2.0, -1.0], 1.0)\n",
" HalfSpace([-0.25, -1.0], 0.0)\n",
" HalfSpace([0.5, -1.0], 0.0):\n",
"5-element iterator of Array{Float64,1}:\n",
" [1.0, 1.0]\n",
" [0.0, 0.0]\n",
" [0.666667, 0.333333]\n",
" [-0.266667, 0.0666667]\n",
" [-1.0, 3.0],\n",
"1-element iterator of Ray{Float64,Array{Float64,1}}:\n",
" Ray([-1.77636e-16, 1.0])"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plotting the epigraph\n",
"\n",
"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$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"plot(p ∩ HalfSpace([0, 1], 4))\n",
"scatter!([x[1] for x in points(p)], [x[2] for x in points(p)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Using the polyhedral function in an optimization problem\n",
"\n",
"[JuMP](https://github.com/jump-dev/JuMP.jl) variables can be constrained to be in the epigraph.\n",
"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.\n",
"One can use `Polyhedra.linear_objective_solver` to get `VRepOptimizer` when the V-representation is computed.\n",
"This optimizer does exactly what we just described."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"using JuMP\n",
"model = Model(Polyhedra.linear_objective_solver(p))\n",
"@variable(model, x)\n",
"@variable(model, y)\n",
"@constraint(model, [x, y] in p)\n",
"@objective(model, Min, x + y)\n",
"optimize!(model)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The optimal solution of this problem is one of the extreme points of `p`:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"OPTIMAL::TerminationStatusCode = 1"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"termination_status(model)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-0.26666666666666666, 0.06666666666666665)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"value(x), value(y)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot!([-0.2, -1.2], [0.0, 1.0], linewidth=2) # Objective function\n",
"scatter!([value(x)], [value(y)], markershape=:star5, markersize=8) # Optimal solution"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"@objective(model, Max, x + y)\n",
"optimize!(model)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we can see, maximizing $x + y$ is an unbounded problem, the infeasibility certificate is given by the extreme ray $(0, 1)$."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"DUAL_INFEASIBLE::TerminationStatusCode = 3"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"termination_status(model)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-1.7763568394002506e-16, 1.0)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"value(x), value(y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"using JuMP\n",
"import GLPK\n",
"model = Model(GLPK.Optimizer)\n",
"@variable(model, x <= -2)\n",
"@variable(model, y)\n",
"@variable(model, 0 <= u <= 3)\n",
"@constraint(model, [2x + u, y] in p)\n",
"@objective(model, Min, u + y)\n",
"optimize!(model)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"OPTIMAL::TerminationStatusCode = 1"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"termination_status(model)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-2.0, 3.0, 3.0)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"value(x), value(y), value(u)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Two dimensional polyhedral function\n",
"\n",
"Consider the polyhedral function\n",
"$$ f(x) = \\max(-4x - 2y - 1, 2x - y - 1, -x/4 + y/2, x/2 + y) $$\n",
"in the interval $(x, y) \\in [-1, 1]^2$.\n",
"\n",
"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\n",
"\\begin{align*}\n",
" -1 & \\le x\\\\\n",
" x & \\le 1\\\\\n",
" -1 & \\le y\\\\\n",
" y & \\le 1\\\\\n",
" -4x - 2y - 1 & \\le z\\\\\n",
" 2x - y - 1 & \\le z\\\\\n",
" -x/4 + y/2 & \\le z\\\\\n",
" x/2 + y & \\le z\n",
"\\end{align*}\n",
"We can create this H-representation in Polyhedra as follows:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"H-representation Polyhedra.Intersection{Float64,Array{Float64,1},Int64}:\n",
"8-element iterator of HalfSpace{Float64,Array{Float64,1}}:\n",
" HalfSpace([-1.0, 0.0, 0.0], 1.0)\n",
" HalfSpace([1.0, 0.0, 0.0], 1.0)\n",
" HalfSpace([0.0, -1.0, 0.0], 1.0)\n",
" HalfSpace([0.0, 1.0, 0.0], 1.0)\n",
" HalfSpace([-4.0, -2.0, -1.0], 1.0)\n",
" HalfSpace([2.0, -1.0, -1.0], 1.0)\n",
" HalfSpace([-0.25, 0.5, -1.0], 0.0)\n",
" HalfSpace([0.5, 1.0, -1.0], 0.0)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Polyhedra\n",
"h = HalfSpace([-1, 0, 0], 1) ∩ HalfSpace([1, 0, 0], 1) ∩\n",
" HalfSpace([0, -1, 0], 1) ∩ HalfSpace([0, 1, 0], 1) ∩\n",
" HalfSpace([-4, -2, -1], 1) ∩ HalfSpace([2, -1, -1], 1) ∩\n",
" HalfSpace([-1/4, 1/2, -1], 0) ∩ HalfSpace([1/2, 1, -1], 0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We create the polyhedron the same way as for 1-dimensional polyhedral functions."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Polyhedron CDDLib.Polyhedron{Float64}:\n",
"8-element iterator of HalfSpace{Float64,Array{Float64,1}}:\n",
" HalfSpace([-1.0, 0.0, 0.0], 1.0)\n",
" HalfSpace([1.0, 0.0, 0.0], 1.0)\n",
" HalfSpace([0.0, -1.0, 0.0], 1.0)\n",
" HalfSpace([0.0, 1.0, 0.0], 1.0)\n",
" HalfSpace([-4.0, -2.0, -1.0], 1.0)\n",
" HalfSpace([2.0, -1.0, -1.0], 1.0)\n",
" HalfSpace([-0.25, 0.5, -1.0], 0.0)\n",
" HalfSpace([0.5, 1.0, -1.0], 0.0)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using CDDLib\n",
"p = polyhedron(h, CDDLib.Library())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We obtain the 10 nodes as follows along with the $(0, 0, 1)$ ray that is expected as the polyhedron is an epigraph."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"V-representation CDDGeneratorMatrix{Float64,Float64}:\n",
"10-element iterator of Array{Float64,1}:\n",
" [0.222222, -0.333333, -0.222222]\n",
" [1.0, 0.25, 0.75]\n",
" [1.0, 1.0, 1.5]\n",
" [-0.666667, 1.0, 0.666667]\n",
" [1.0, -1.0, 2.0]\n",
" [0.0888889, -0.533333, -0.288889]\n",
" [0.166667, -1.0, 0.333333]\n",
" [-0.933333, 1.0, 0.733333]\n",
" [-1.0, 1.0, 1.0]\n",
" [-1.0, -1.0, 5.0],\n",
"1-element iterator of Ray{Float64,Array{Float64,1}}:\n",
" Ray([-6.66134e-16, -7.40149e-17, 1.0])"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vrep(p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As `p` is a 3-dimensional polyhedron, we use [MeshCat](https://github.com/rdeits/MeshCat.jl) or [Makie](https://github.com/JuliaPlots/Makie.jl).\n",
"Note that plotting unbounded polyhedra is supported by Polyhedra so no need to add, e.g. the hyperplane $z \\le 6$.\n",
"The top of the shape will have no face as it is unbounded in this direction."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m = Polyhedra.Mesh(p)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using MeshCat\n",
"vis = Visualizer()\n",
"setobject!(vis, m)\n",
"IJuliaCell(vis)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using Makie\n",
"Makie.mesh(m, color=:blue)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using Makie\n",
"Makie.wireframe(m)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.0.3",
"language": "julia",
"name": "julia-1.0"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.0.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}