{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Disclaimer\n",
"This notebook is only working under the versions:\n",
"\n",
"- JuMP 0.19 (unreleased, but currently in master)\n",
"\n",
"- MathOptInterface 0.4.1\n",
"\n",
"- GLPK 0.6.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Description**: Shows how to find the Chebyshev center of a polyhedron.\n",
"\n",
"**Author**: Iain Dunning\n",
"\n",
"**License**: ![\"Creative](\"https://i.creativecommons.org/l/by-sa/4.0/88x31.png\")
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Finding the Chebyshev Center using JuMP\n",
"\n",
"*This model is drawn from section 4.3.1 of the book [**Convex Optimization** by Boyd and Vandenberghe](http://stanford.edu/~boyd/cvxbook/)*\n",
"\n",
"Given a polyhedron described as a set of hyperplanes, i.e.\n",
"$$\n",
"P = \\left\\{ x \\mid A x \\leq b \\right\\}\n",
"$$\n",
"find the largest Euclidean ball that lies inside the polyhedron. The center of this ball is known as the **Chebyshev center** of the polyhedron.\n",
"\n",
"This problem can be formulated a linear optimization problem. Let $x$ be the center of the ball, and $r$ be the radius of the ball. Our ball can then be expressed as\n",
"$$\n",
"B = \\left\\{ x + u \\mid \\|u\\|_2 \\leq r \\right\\}\n",
"$$\n",
"and we can express the constraint that the ball lies in a halfspace as $A_i^\\prime x \\leq b_i$, i.e.\n",
"$$\n",
"A_i^\\prime \\left ( x + u \\right) \\leq b_i \\quad \\forall u:\\|u\\|_2 \\leq r\n",
"$$\n",
"\n",
"We can then use the fact that\n",
"$$\n",
"\\sup_{\\|u\\|_2 \\leq r} A_i^\\prime x = r \\| A_i \\|_2\n",
"$$\n",
"to simplify that constraint to the linear equality\n",
"$$\n",
"A_i^\\prime x + r \\| A_i \\|_2 \\leq b_i\n",
"$$\n",
"which allows us to express the problem as\n",
"$$\n",
"\\begin{alignat}{2}\n",
" \\max \\ & r \\\\\n",
"\\text{subject to} \\ & A_i^\\prime x + \\left\\| A_i \\right\\| r \\leq b_i\n",
"\\end{alignat}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"MathOptInterface.Utilities"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Load JuMP\n",
"using JuMP\n",
"using MathOptInterface\n",
"# Load solver package\n",
"using GLPK\n",
"# shortcuts\n",
"const MOI = MathOptInterface\n",
"const MOIU = MathOptInterface.Utilities"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Use a hyperplane representation of the polyhedron\n",
"# i.e. P = { x | Ax ≤ b }\n",
"n = 4\n",
"# Store LHS as vector-of-vectors\n",
"A = Vector{Float64}[\n",
" [ 2, 1], [ 2,-1],\n",
" [-1, 2], [-1,-2]]\n",
"b = ones(n)\n",
"\n",
"# Build JuMP model\n",
"m = Model(optimizer = GLPK.GLPKOptimizerLP())\n",
"@variable(m, r)\n",
"@variable(m, x[1:2])\n",
"@objective(m, Max, 1r)\n",
"@constraint(m, P[i=1:4], sum(A[i][j]*x[j] for j in eachindex(x)) + norm(A[i])*r <= b[i])\n",
"# solve problem\n",
"JuMP.optimize(m)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[-0.0, -0.0]"
]
}
],
"source": [
"# Where is the center?\n",
"print(JuMP.resultvalue.(x))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.4472135954999579"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"JuMP.resultvalue(r)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Use Gadfly to display the solution\n",
"using Gadfly\n",
"Gadfly.set_default_plot_size(8cm, 8cm)\n",
"# Plot lines over [-1.5, 1.5]\n",
"domain = linspace(-1, +1)\n",
"# Plot circle across all angles\n",
"θ = linspace(0,2π)\n",
"plot(\n",
"# a_1 x_1 + a_2 x_2 = b\n",
"# --> x_2 = (b - a_1 x_1)/a_2\n",
"[layer(x=domain,\n",
" y=(b[i]-A[i][1]*domain)/A[i][2],\n",
" Geom.line,\n",
" Theme(line_width=2px,\n",
" default_color=colorant\"blue\")) for i in 1:4]...,\n",
"# x_1' = x_1 + rθ\n",
"# x_2' = x_2 + rθ\n",
"layer(x=JuMP.resultvalue(x[1]) + JuMP.resultvalue(r)*cos.(θ),\n",
" y=JuMP.resultvalue(x[2]) + JuMP.resultvalue(r)*sin.(θ),\n",
" Geom.path,\n",
" Theme(line_width=5px,\n",
" default_color=colorant\"red\")),\n",
"Coord.Cartesian(ymin=-1,ymax=+1)\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.6.0",
"language": "julia",
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.6.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}