Note: This post was written before Gurobi supported nonlinear optimization. It has been updated to work with Python 3.
A common problem in handling geometric data is determining the center of a given polygon. This is not quite so easy as it sounds as there is not a single definition of center that makes sense in all cases. For instance, sometimes computing the center of a polygon’s bounding box may be sufficient. In some instances this may give a point on an edge (consider a right triangle). If the given polygon is non-convex, that point may not even be inside or on its boundary.
This post looks at computing Chebyshev centers for arbitrary convex polygons. We employ essentially the same model as in Boyd & Vandenberghe’s Convex Optimization text, but using Gurobi instead of CVXOPT.
Consider a polygon defined by the intersection of a finite number of half-spaces, $Au \le b$. We assume we are given the set of vertices, $V$, in clockwise order around the polygon. $E$ is the set of edges connecting these vertices. Each edge in $E$ defines a boundary of the half-space $a_i^\intercal u \le b_i$
$$ V = {(1,1), (2,5), (5,4), (6,2), (4,1)}\\ E = {((1,1),(2,5)), ((2,5),(5,4)), ((5,4),(6,2)), ((6,2),(4,1)), ((4,1),(1,1))} $$
The Chebyshev center of this polygon is the center point $(x, y)$ of the maximum radius inscribed circle. That is, if we can find the largest circle that will fit inside our polygon without going outside its boundary, its center is the point we are looking for. Our decision variables are the center $(x, y)$ and the maximum inscribed radius, $r$.
In order to do this, we consider the edges independently. The long line segment below shows an arbitrary edge, $a_i^\intercal u \le b_i$. The short line connected to it is orthogonal in the direction $a$. $(x, y)$ satisfies the inequality.
The shortest distance from $(x, y)$ will be in the direction of $a$. We’ll call this distance $r$. If we were to move the edge so it had the same slope but went through $(x, y)$, its distance from $a_i^\intercal u = b_i$ would be $r||a_i||_2$. Thus we can add a constraint of the form $a_i’u + r||a_i||_2 \le b_i$ for each edge and maximize the value of $r$ as our objective function.
$$ \begin{align*} & \text{max} && r \\ & \text{s.t.} && (y_i-y_j)x + (x_j-x_i)y + r\sqrt{(x_j-x_i)^2 + (y_j-y_i)^2} \le (y_i-y_j)x_i + (x_j-x_i)y_i \\ & && \quad \forall \quad ((x_i,y_i), (x_j,y_j)) \in E \\ \end{align*} $$
As this is linear, we can solve it using any LP solver. The following code does so with Gurobi.
#!/usr/bin/env python3
from gurobipy import Model, GRB
from math import sqrt
vertices = [(1,1), (2,5), (5,4), (6,2), (4,1)]
edges = zip(vertices, vertices[1:] + [vertices[0]])
m = Model()
r = m.addVar()
x = m.addVar(lb=-GRB.INFINITY)
y = m.addVar(lb=-GRB.INFINITY)
m.update()
for (x1, y1), (x2, y2) in edges:
dx = x2 - x1
dy = y2 - y1
m.addConstr((dx*y - dy*x) + (r * sqrt(dx**2 + dy**2)) <= dx*y1 - dy*x1)
m.setObjective(r, GRB.MAXIMIZE)
m.optimize()
print('r = %.04f' % r.x)
print('(x, y) = (%.04f, %.04f)' % (x.x, y.x))
The model output shows our center and its maximum inscribed radius.
$$ r = 1.7466\\ (x, y) = (3.2370, 2.7466) $$
Question for the reader: in certain circumstances, such as rectangles, the Chebyshev center is ambiguous. How might one get around this ambiguity?