A small polygon is a unit diameter polygon. Finding the polygon with maximum area is known as the “Largest Small Polygon” (LSP) problem. We are proposing some non-linear formulations to find the small polygon that maximizes the area for a given number of sides to solve the problem computationally. This is an open geometry problem, so we are providing a computational and approximate way to tackle it.
Problem statement
Given a number n, find the n-sided polygon that has diameter one and that has the largest area among all diameter-one n-gons. The diameter of a polygon is the maximum distance between any two vertices in the polygon, and distances between points are assumed to be Euclidean distances.
Some theoretical results
The LSP problem has been solved for odd number of sides. In 1922, Reinhardt showed that the optimal solution for $n = 2k+1$ is the regular n-gon. The area of the regular unit n-gon is:
$$\frac{n}{2}(sin\frac{\pi}{n}-tan\frac{\pi}{2n}) \text{ , when n is odd}$$
$$\frac{n}{8}sin\frac{2\pi}{n} \text{ , when n is even}$$
These are lower bounds for our problem, but the even one might be improved from the regular polygon of $n-1$ sides:
Take the odd regular $n-1$-gon. Add a vertex (at distance 1) along the bisectriz of any angle.
The obtained polygon will be larger than the even regular $n$-gon. This fact proves that the regular polygon is not optimal for $n \geq 6$.
Although LSP has not been solved for general even number of sides, there are some conjectures related to the solutions and variants of the problem that assume some kind of symmetry in terms of the diameter graph of the polygon (the diameter graph of a small polygon is the graph with the vertices of the polygon, where vertices are connected if and only if the distance between them is equal to 1).
Model with polar coordinates
The problem can be stated and solved by describing the problem in an straightforward way, formulating a non linear problem in polar coordinates. We are proposing a similar model to the one in the netlib repository, an AMPL model (by David M. Gay) that started as a GAMS model by Francisco J. Prieto. This model solves the general n instance of the problem.
Variables
If $(x_i, y_i)$ are the cartesian coordinates (real numbers) of the points, they will be the vertices of a small polygon if the distances between them are less or equal to 1.
$$\sqrt{(x_i-x_j)^2+(y_i-y_j)^2} \leq 1$$
That is equivalent to:
$$(x_i-x_j)^2+(y_i-y_j)^2 \leq 1$$
Compute area
Given the n vertices of a polygon, we can compute the area of the polygon by using the Shoelace formula, like a 2D determinant:
Example:
The area of the polygon defined with the vertices
{(-1,0), (0.2,1), (2,0.7), (2,-1)} is
$$\frac{1}{2}|-1\times1+0.2\times0.7+2\times(-1)+2\times0-0\times0.2-1\times2-0.7\times2-(-1)\times(-1)| = 3.63$$
The formula works when the points are sorted counter clock-wise or clock-wise, then, in order to use the formula is is necessary to ensure the order of the points along the perimeter. How could we impose an order among the vertices, so first vertex (x1, y1) is connected to the second (x2, y2), which at the same time is connected to (x3, y3)… and so on until closing the polygon?
$$\frac{1}{2}\sum_{i=1}^{n-1}(x_i \cdot y_{i+1} - y_i \cdot x_{i+1})+\frac{1}{2}(x_n \cdot y_1 - y_n \cdot x_1)$$
Polar coordinates to sort vertices
We need to sort vertices somehow. Intuitively, we could sort them by angle respect to any point in the plane.
We are using polar coordinates to sort vertices by angle in a straightforward way. We are sorting respect to (0,0). Each cartesian point $(x,y)$ distinct to (0,0) is matched to a polar point $(r, t)$, where $r$ is the radius of the point respect to the origin, and $t$ the angle respect to the $X$ axis, with $0 < r$, and $0 \leq t < 2\pi$.
$$x = r \cdot cos(t)$$
$$y = r \cdot sin(t)$$
We can substitute $x$ and $y$ as a function of $r$ and $t$, adding sorting constraints as:
Two vertices having the same angle would cause a suboptimal solution as there would be three points in the same line (origin vertex and two others).
Complete model
Substitute $x$ and $y$ by their polar transformation coordinates. Notice that is is assumed vertex (x1,y1) = (0,0) = (r1, t1), as we can always fix one of the vertices.
With the $\pi$ bound over $t$ we force the polygon to happen over the first semiplane (y >= 0).
$$\frac{1}{2}\sum_{i=1}^{n-1}(r_i \cdot cos(t_i) \cdot r_{i+1} \cdot sin(t_{i+1}))-(r_i \cdot sin(t_i) \cdot r_{i+1} \cdot cos(t_{i+1}))+\frac{1}{2}(r_n \cdot cos(t_n) \cdot r_1 \cdot sin(t_1))-(r_n \cdot sin(t_n) \cdot r_1 \cdot cos(t_1))$$
We can simplify the previous expression. Rearrange each term:
$$r_i \cdot cos(t_i) \cdot r_{i+1} \cdot sin(t_{i+1})-r_i \cdot sin(t_i) \cdot r_{i+1} \cdot cos(t_{i+1}) = $$
$$ = r_i r_{i+1} (cos(t_i) \cdot sin(t_{i+1})-sin(t_i) \cdot cos(t_{i+1})$$
Sine of the sum of $t_{i+1}$, $-t_{i}$:
$$r_i r_{i+1} sin(t_{i+1}-t_i)$$
The objective can be rewritten as:
$$\frac{1}{2}\sum_{i=1}^{n-1}(r_i r_{i+1} sin(t_{i+1}-t_i))+\frac{1}{2}(r_1 r_n \cdot sin(t_1-t_n))$$
Distance between vertices. For all $i = 1..n, j = i+1..n$,
$$(r_i \cdot cos(t_i)-r_j \cdot cos(t_j))^2+(r_i \cdot sin(t_i)^2-r_j \cdot sin(t_j))^2 \leq 1$$
$$(r_i \cdot cos(t_i)-r_j \cdot cos(t_j))^2 = r_i^2 cos^2(t_i) + r_j^2 cos^2(t_j) - 2r_ir_jcos(t_i)cos(t_j)$$
$$(r_i \cdot sin(t_i)-r_j \cdot sin(t_j))^2 = r_i^2 sin^2(t_i) + r_j^2 sin^2(t_j) - 2r_ir_jsin(t_i)sin(t_j)$$
Sum both:
$$r_i^2 + r_j^2 - 2r_ir_j(cos(t_i)cos(t_j) + sin(t_i)sin(t_j)) = r_i^2 + r_j^2 - 2r_ir_j cos(t_i-t_j)$$
So the constraints could be rewritten as:
$$r_i^2 + r_j^2 - 2r_ir_j cos(t_i-t_j) \leq 1$$
This way, non-linear expressions in the model are simpler, there are less sines and cosines, so it will be easier to find better solutions more efficiently.
Angle order. For all $i = 2..n$,
$$t_i >= t_{i-1}$$
Finally, we can fix one of the vertices:
$$r_1 = 0, \ t_1 = 0$$
AMPL model
An area_lb parameter is introduced so we can add better lower bounds for the area. This way it is easier to deal with suboptimal and degenerated solutions. By default, the lower bound is 0.
Overwriting polar_polygon.mod
We use ipopt (version 3.12.13) open source solver for cloud usage. The output is very verbose, so we set “outlev=0”.
Ipopt 3.12.13: outlev=0
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
Ipopt 3.12.13: Optimal Solution Found
suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
Now that we have the solution, we are going to plot it with matplotlib.
We define the solve_polar() method to call AMPL and solve the problem for $N$ sides. polar2cart() transforms from polar to cartesian coordinates. plot_polygon() draws the result.
In the below loop, solve the problem for all the sides between min_sides and max_sides.
Problem size: 10
Ipopt 3.12.13: outlev=0
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
Ipopt 3.12.13: Optimal Solution Found
suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
Problem size: 11
Ipopt 3.12.13: outlev=0
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
Ipopt 3.12.13: Optimal Solution Found
suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
Problem size: 12
Ipopt 3.12.13: outlev=0
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
Ipopt 3.12.13: Optimal Solution Found
suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
Sides: 10
Area: 0.7491373270801958
Sides: 11
Area: 0.7587484213086267
Sides: 12
Area: 0.760729846642839
Second model with cartesian coordinates
There is a way to deal with orientation without using polar coordinates. It is easy to check if three points within the plane are clockwise or counter-clockwise oriented. We are going to undo the polar transformation and use new angle constraints, so that the objective and the distance constraints get simplified.
Three points orientation
Given three points, $A=(A_x, A_y)$, $B=(B_x,B_y)$, $C=(C_x,C_y)$, we can compute the slopes of the vectors $AB$ and $BC$.
$$AB_s = \frac{B_y-A_y}{B_x-A_x}$$
$$BC_s = \frac{C_y-B_y}{C_x-B_x}$$
If both slopes are positive, BC’s slope being less (greater) than AB’s implies that the three points will be (counter-)clockwise oriented. The comparison could be flipped if some of the slopes is negative, but in the end, the following expression let’s us know the orientation:
$$AB_s-BC_s = \frac{B_y-A_y}{B_x-A_x}-\frac{C_y-B_y}{C_x-B_x} = \frac{(B_y-A_y)(C_x-B_x)-(C_y-B_y)(B_x-A_x)}{(B_x-A_x)(C_x-B_x)}$$
If the numerator is positive, points are clockwise.
If the numerator is negative, points are counter-clockwise.
If the numerator is 0, points belong to the same line.
Orientation in the model
As we are fixing one of the vertices (0,0), we can force all the other points to have the same orientation. For $i = 2..N-1$, we can force the same orientation for $A=(0,0)$, $B=(x_i,y_i)$, and $C=(x_{i+1},y_{i+1})$. This way, vertices are going to end up ordered. By using the previous formula, this is possible if we have the same sign in the following expression for each $i$:
$$(B_y-A_y)(C_x-B_x)-(C_y-B_y)(B_x-A_x) = y_i(x_{i+1}-x_i)-(y_{i+1}-y_i)x_i = $$
$$ = y_i x_{i+1}- y_i x_i - x_i y_{i+1} + x_i y_i = y_i x_{i+1} - x_i y_{i+1}$$
Then, we can simply add the following constraint to order our vertices counter-clockwise. For each $i = 2..N-1$:
$$y_i x_{i+1} - x_i y_{i+1} < 0$$
$$x_i y_{i+1} - y_i x_{i+1} > 0$$
By simplicity, in the model we are using the non-strict inequality. Three vertices should not be places in the same segment as that is a suboptimal solution. It would be suboptimal as the point in the middle is pointless, but the polygon could be bigger by using it (just place it in the mediatrix of two other vertices).
For further information about this formula and orientation of polygons, take a look to Graham’s scan Algorithm, where the formula is used to compute convex hulls.
Area computation with triangles
Now, we can use cartesian coordinates as variables, and simplify constraints and objective (by removing sines and cosines).
If we project a ray from the fixed vertex of the polygon to each other vertex, we are generating triangles whose area is (half of the area of the parallelogram computed via scalar product):
$$u_i = \frac{1}{2}(x_i y_{i+1} - y_i x_{i+1})$$
The area of the polygon is equal to the sum of areas of these triangles, so we could add continuous variables $u_i$ for $i = 2..N-1$, where $u_i$ is the area of the triangle that contains vertices $(0,0)$, $(x_i,y_i)$, and $(x_{i+1},y_{i+1})$. These variables are non-negative, and at optimality, they should be equal to the area of their triangle.
Then, the objective function is rewritten as:
$$maximize \quad \sum_{i=2..N-1} u_i$$
$u_i$ should not exceed the area:
$$u_i \leq \frac{1}{2}(x_i y_{i+1} - y_i x_{i+1})$$
$$2u_i \leq x_i y_{i+1} - y_i x_{i+1}$$
So we are replacing orientation constraints with:
$$0 \leq 2u_i \leq x_i y_{i+1} - y_i x_{i+1}$$
This idea is used in state-of-the-art models such as the one published by Christian Bingane (2022), Largest small polygons: A sequential convex optimization approach.