Optimal Area Polygonisation Problems: Mixed Integer Linear Programming Models

Optimal Area Polygonisation Problems: Mixed Integer Linear Programming Models Hip´olito Hern´andez-P´erez, Jorge Riera-Ledesma, Inmaculada Rodr´ıguez-Mart´ın, Juan-Jos´e Salazar-Gonz´alez IMAULL, Universidad de La Laguna, 38200 Tenerife, Spain Abstract This paper describes approaches to finding polygons with maximum and minimum area through a given set of points on the Euclidean plane. Both problems are NP-hard and of interest in Computational Geometry. They have been extensively studied, although only the recent literature presents mathematical formulations to consistently find optimal solutions to instances from a benchmark collection with up to 25 points. We present and compare eight new Mixed Integer Linear Programming models for the two problems. Four represent a polygon as a solution to the Asymmetric Travelling Salesman Problem, while the other four use the Symmetric Travelling Salesman Problem representation. Six of the models select triangles to compute the area of a polygon, while the other two models select subpolygons. A computational analysis on the benchmark collection shows that one of our new formulations can consistently solve most of the benchmark instances with up to 50 points. The success of our proposal is a linear system to select a triangulation of the point set, a partition of the triangulation to cover the internal and external area of a polygon, and a representation of the polygon as the solution of an Asymmetric Travelling Salesman Problem. Keywords: Computational Geometry, Triangulations, Polygonisations, Mixed Integer Linear Programming, Travelling Salesman Problem. ⋆This work has been partially supported by the Spanish Government through the research grants PID2023-148599NB-I00 and PCI2024-155092-2. Email addresses: [email protected] (Hip´olito Hern´andez-P´erez), [email protected] (Jorge Riera-Ledesma), [email protected] (Inmaculada Rodr´ıguez-Mart´ın), [email protected] (Juan-Jos´e Salazar-Gonz´alez) Preprint submitted to European Journal of Operational Research November 11, 2025 1. Introduction Given a finite set of points and a cost matrix representing the travel cost between each pair of points, the well-known Travelling Salesman Problem (TSP) looks for a sequence of the points such that a closed route (tour) following that sequence has the minimum cost. When the cost matrix is a metric between the points, the problem is called Metric TSP. A particular case of the Metric TSP occurs when the metric is the Euclidean norm, happening when the points are in the Cartesian plane, each one iis associated with coordinates (αi, βi), and the travel cost between point iand point jis p(αi−αj)2+ (βi−βj)2. The problem is then called Euclidean TSP and aims at computing a polygon whose vertices are the given points and with the shortest perimeter. The computational complexity of the Euclidean TSP is interesting since it is an N P-hard problem (see Papadimitriou (1977)), but it is not in NP because the square root of a number may be irrational and, therefore, there is no efficient procedure to compute the total travel cost of a tour. Another interesting property of the Euclidean TSP is that optimal solutions are simple polygons, meaning their line segments may intersect only at the given points. This property has been used in the literature to describe Polynomial-Time Approximation Schemes (i.e. algorithms that compute, in time polynomial in the size of the problem instance, a solution within a factor of (1 + ϵ) of the optimal solution for any fixed ϵ > 0). The property does not hold when the objective function sense is reversed and the aim is to compute a polygon with the longest perimeter (known as Maximum TSP). Our manuscript concerns two problems closely related to the Euclidean TSP and the Maximum TSP. They consist of finding sequences of the given points such that the internal area is minimal in one problem and maximal in the other problem. The point sequence of points, seen above as a tour, is now seen as a polygon. The optimisation problems are named MinArea and MaxArea, respectively, and both together Optimal Area polygonisation Problem (OAP), of interest in Computational Geometry (see e.g. Mark Berg and Schwarzkopf (1997)). Computational Geometry delves into various challenges surrounding the construction of geometric objects from a finite set of points on a plane. Tasks like triangulations of point sets, Voronoi diagrams, arrangements, etc., can be considered within this field of research. As our paper remarks, particularly close to the OAP are the Triangulation Problems, where one must look for triangles with vertices spanning a given point set (see e.g. Dantzig et al. (1985), Loera et al. (2010)). Some problems are optimisation problems, like the Minimum Weight Triangula2 tion Problem (see e.g. Kyoda et al. (1997)), where the weight of a triangle is the perimeter of the triangle, and the aim is to find a triangulation with minimal total weight. Other problems are feasibility questions, like checking whether there exists a triangulation using only line segments in a given set. The mentioned optimisation problem is NP-hard (see Mulzer and Rote (2008)), and the mentioned feasibility problem is N P-complete (see Lloyd (1977)). Fekete (2000) shows that MinArea and MaxArea (and generalisations in higher-dimensional spaces) are NP-complete. Both problems garnered significant attention during a global competition called “Computational Geometry: Solving Hard Optimisation Problems” (CGSHOP 2019)1. Held in 2019, the contest tasked participants with developing effective algorithms to solve (exactly and/or heuristically) a diverse range of benchmark instances. Seven papers were published in a special issue of the ACM Journal of Experimental Algorithmics devoted to this challenge. Demaine et al. (2022) surveyed publications related to the MinArea and MaxArea problems until 2019. Five papers were devoted to heuristic algorithms, and two papers were devoted to exact algorithms. Heuristic procedures are described by Crombez et al. (2022), Goren et al. (2022), Eder et al. (2022), and Ramos et al. (2022a). Fekete et al. (2022) and Ramos et al. (2022b) solve exactly the MinArea and MaxArea problems. Fekete et al. (2022) present exact algorithms based on two Mixed Integer Linear Programming (MILP) models. The first formulation models the polygon as a counterclockwise-oriented circuit visiting all the given points, and it needs an additional point (reference point) to define the cost of an arc equal to the area between this arc and the reference point, with the sign depending on whether the arc is clockwise oriented or anticlockwise oriented. The second formulation aims at selecting triangles to cover the interior of the optimal polygon, and the major drawback is the connectivity constraints between the selected triangles. Ramos et al. (2022b) present an alternative approach to force the connectivity constraints, consisting of finding a tree spanning all the selected triangles in a graph. While computationally, the triangle-based model showed better results than the other formulations, it is based on O(n6) variables. The formulation is the best approach in the literature. On a benchmark collection of instances for the CGSHOP 2019 challenge, this formulation allowed finding all optimal solutions for instances with 10 and 15 points, some with 20, a few with 25, and none with 30. 1https://cgshop.ibr.cs.tu-bs.de/competition/cg-shop-2019/ 3 Our paper describes eight new MILP models to solve the MinArea and MaxArea problems. Four represent the polygon as an oriented tour, and the other four as a non-oriented tour. Six models compute the objective function based on selecting triangles, either inside the polygon, outside the polygon, or both. The other two models replace triangles with subpolygons. After having implemented all variants, there is one that succeeded in finding optimal solutions to all the benchmark instances with up to 30 points. Even more, we apply this new model to larger benchmark instances, and it finds the optimal solutions to all instances with up to 45 points, and most of the instances with 50 points. Our major contribution is a set of linear systems to model the selection of triangles and subpolygons. The polygon is modelled with TSP constraints, and the non-crossing edges requirement is implicit in the linear systems. The remainder of this paper is organised as follows. Section 2 gives some definitions and notations of MinArea and MaxArea. Section 3 describes two MILP models taken from the literature. Section 4 introduces our new MILP formulations. Section 5 describes valid inequalities and other implementation details. Section 6 analyses computational experiments on benchmark instances. The paper ends with conclusions. 2. Definitions and notations This section introduces the notation needed in the rest of the paper. We assume to be working on the Cartesian plane R2with the Euclidean norm, and we are given the Cartesian coordinates of the points in N. Without loss of generality, and in coherence with benchmark instances in the literature, the coordinates are assumed to be integer numbers. The set of all convex combinations of points in N, known as the convex hull of N, is denoted by CH(N). Let N′be the subset of Non the boundary of CH(N). Given two points iand jon the plane, CH({i, j}) is the line segment between these points and, for simplicity, we denote it by [i, j] or [j, i] indistinctly. A polygon Pis a region of the plane bounded by a finite collection of line segments forming a simple closed curve which is the homeomorphic image of a disk, i.e. it is a certain deformation of a disk. Each line segment on the polygon’s border is called an edge, and the extreme point of an edge is called a vertex. Given N, a polygonisation of Nis a polygon Pwhose vertex set is N. A polygon is simple if its edges’ intersections occur only at its vertices. Thus, an edge of a simple polygon cannot include a third point of N. For the rest of this manuscript, we refer to simple polygons only. 4 0 1 2 34 5 6 78 9 (a) MinArea. Area = 9096, Perimeter = 854.50355 0 1 2 34 5 6 78 9 (b) MaxArea. Area = 17532, Perimeter = 659.91797 0 1 2 34 5 6 78 9 (c) Euclidean TSP. Area = 16096, Perimeter = 589.25625 0 1 2 34 5 6 78 9 (d) Maximum TSP. Area = 11974, Perimeter = 928.97496 Figure 1: Optimal solutions for the instance uniform-0000010-1. The Jordan Curve Theorem states that a polygon Pdivides CH(N) into two disjoint regions: the interior and the exterior. A line segment connecting two vertices of Pand not crossing an edge of Pis called an internal diagonal when the whole segment lies in the interior of Pand an external diagonal when it lies in the exterior. In this paper, the area of the polygon Pmeans the area c(P) of its interior region. We assume N= N′, thus |N|>|N′|≥ 3, any polygonisation Pof Nis non-convex, and c(CH(N)) > c(P)>0. MinArea and MaxArea look for a polygonisation Pof Nwith minimum and maximum area c(P), respectively. Figure 1 shows optimal solutions for the MinArea and MaxArea problems with ten points from Ramos et al. (2022b). Optimal solutions for the Euclidean TSP and the Maximum TSP with non-crossing edges are also displayed. The four pictures confirm that, while the problems are related, they may all have different optimal objective values (and therefore different optimal solutions). 5 Triangles are useful objects for computing the area of a polygon. Let V be the set of all triangles whose three vertices are in N, discarding triangles containing a fourth point in Nand triangles with the three vertices in a line. We assume that each triangle is counterclockwise oriented, meaning that a triangle with vertices {i, j, k} ⊂ Ncan be uniquely represented by three (oriented) arcs {(i, j),(j, k),(k, i)}. Let Vij be the subset of triangles in V containing the arc (i, j). Graphically, Vij and Vji represent the triangles on the left and right, respectively, of a line segment [i, j] with i, j ∈N. Trivially, |Vij ∪Vji|≤ |N|−2. For some formulations in our paper, a polygon is represented as a Hamiltonian circuit with non-crossing arcs on a directed graph G(N, A) where N is the node set, and the above-introduced arcs define A. The previous assumption on the triangles’ counterclockwise orientation helps us to have a unique representation of a polygon: a tour is a counterclockwise Hamiltonian circuit in G(N, A). Following the points along a tour, the interior region of the polygon is on the left, and the exterior region is on the right. Let A′be the subset of Acorresponding to line segments on the boundary of CH(N); these arcs determine a counterclockwise circuit spanning N′. Let A′′ =A\A′represent the set of all possible (internal or external) diagonals. (|A′′|≥ 1.) For other formulations in our paper, a polygon is represented as a nonoriented Hamiltonian tour with non-crossing edges in the undirected graph G(N, E) where Nis the node set and E={{i, j}: (i, j)∈A}. Let E′= {{i, j} ∈ E: (i, j)∈A′}be the edges on the boundary of CH(N), and E′′ =E\E′be the set of all possible (internal or external) diagonals. (|E′′|≥ 1.) The constraint that the polygon must be simple (called non-crossing edges requirement) implies an incompatibility between some pairs of line segments. Precisely, a pair of line segments [i, j] and [k, l] are incompatible when there exists a common point to both, not in {i, j, k, l}. We denote by I2the set of all pairs of incompatible line segments with vertices in N. Observe that no pair in I2includes a line segment on the border of CH(N), meaning that {(i, j),(j, i),(k, l),(l, k)} ⊂ A′′ when {[i, j],[k, l]} ∈ I2. Two triangles sand tare incompatible when the interior regions of sand tare not disjoint, and they are adjacent when there is an edge {i, j} ∈ E′′ such that s∈Vij and t∈Vji. We denote by I3the set of all pairs {s, t} being sand tincompatible triangles in V, and Fthe set of all pairs {s, t} being sand tadjacent triangles in V. We denote by G(V, F) the undirected graph that has a node for each triangle t∈Vand an edge e={s, t} ∈ Ffor each pair of adjacent triangles sand t. Observe that G(V, F) is a connected 6 graph where a node’s degree is at most three. This paper uses the following standard graph notation. For S⊂N, let A(S) = {(i, j)∈A:i, j ∈S},δ+(S) = {(i, j)∈A:i∈S, j ∈N\S}and δ−(S) = {(i, j)∈A:i∈N\S, j ∈S}. For S⊂N, let E(S) = {{i, j} ∈ E: i, j ∈S}and δ(S) = {{i, j} ∈ E:i∈S, j ∈N\Sor i∈V\S, j ∈S}. A triangulation Tof a polygon Pis a pair-disjoint subset of triangles such that the union is the interior area of P. In this paper, we call this triangulation internal because we also consider another triangulation T′covering the external area of the polygon within the convex hull, which we call the external triangulation of P. A triangulation of the set of points Nis defined in the literature as a pair-disjoint subset of triangles such that the union is CH(N). Thus, the union of the internal and external triangulations of a polygon Pis a triangulation of its set of vertices N. Any internal triangulation Tof Pincludes |N|−2 triangles and |N|−3 non-crossing internal diagonals. Any external triangulation T′of Pincludes |N|−|N′|triangles and |N|−|N′|non-crossing external diagonals. Thus, a triangulation of N contains 2|N|−|N′|−2 triangles and involves 3|N|−|N′|−3 non-crossing line segments. It is also folklore (see e.g. Loera et al. (2010)) that the dual graph of a triangulation Tof P, which is the subgraph of G(V, F ) induced by T, is a tree (with each node of degree at most three). This tree can be extended to have the following result: Lemma 1. A polygon triangulation corresponds to a binary tree where each internal node (with degree three) corresponds to a triangle and each leaf node (with degree one) corresponds to an edge of the polygon. The above result refers to internal triangulations, but it also applies to external triangulations, replacing the word “tree” with “forest”. Figure 2 shows two examples, with Tin grey colour, T′in white colour, and the binary trees from Tin red colour. Triangulations are of interest when solving OAP for two reasons. One reason is that the area of a polygon reduces to the summation of the triangles’ areas. Given a triangle t∈V, let ctbe the area of t. Then, the area of a polygon Pcan be computed as c(P) = Pt∈Tctfor any internal triangulation Tof P. Alternatively, if T′is an external triangulation of P, then c(P) = c(CH(N)) −Pt∈T′ct. Another reason, and more important, is that triangulations suggest a nice alternative to model the non-crossing edges requirement, as Section 4 shows. Before, Section 3 summarises how this requirement is formulated in previous articles. 7 0 1 2 34 5 6 78 9 (a) MinArea 0 1 2 34 5 6 78 9 (b) MaxArea Figure 2: Examples of triangulations of the polygons in Figures 1a and 1b. 3. Known models This section presents mathematical formulations taken from the recent literature on the MinArea and the MaxArea. One formulation is based on the directed graph G(N, A), determining the optimal polygon by selecting arcs in the graph to have a Hamiltonian circuit. Other formulations are based on the larger undirected graph G(V, E), modelling the optimal polygon as a selection of triangles. For this reason, the first one is called the arc-based formulation, and the others are called the triangle-based formulations. 3.1. Arc-Based Formulation The first mathematical formulations for the MinArea and the MaxArea were introduced by Fekete et al. (2022) with the name “Edge-Based Formulation”. It is renamed here as “Arc-Based Formulation” because it makes use of binary variables related to the arcs of a directed graph, and we introduce a different formulation in Section 4.5 using binary variables related to edges of an undirected graph. We now present the formulations in Fekete et al. (2022). Consider a given finite set Nof points and one additional (arbitrary but fixed) reference point r∈ N. For any pair i, j in N, let cij be the area of the triangle if {(i, j),(j, r),(r, i)}is counterclockwise oriented, and minus this number otherwise. It is known in the literature as the signed area. Then, the area of any (simple) polygonisation of Ncan be computed by considering the polygon as a circuit and adding cij for each arc (i, j) along this circuit. The result is a positive or negative number depending on whether the circuit orientation is counterclockwise or clockwise, respectively. Hence, MaxArea 8 can be mathematically formulated as the Asymmetric TSP with no crossing arcs. More precisely, for each a∈A, let xabe a binary variable to determine whether ais in the circuit representing the polygonisation. A formulation for MaxArea is: min/max X a∈A caxa(1) subject to: X a∈δ+(i) xa=X a∈δ−(i) xa= 1 for all i∈N(2) X a∈A(S) xa≤ |S|−1 for all S⊂N(3) xa∈ {0,1}for all a∈A(4) xij +xji +xkl +xlk ≤1 for all {[i, j],[k, l]} ∈ I2.(5) Constraints (2)–(4) force a Hamiltonian circuit on the directed graph G(N, A) and Inequalities (5) guarantee the non-crossing edge requirement. Both maximising and minimising in (1) end with polygons with maximum area, and the difference is the circuit’s orientation. By our definition of Ain Section 2, where the reverse arcs in A′are excluded from A, the two optimisation problems may find different optimal polygons. But still, the above model cannot solve the MinArea without additional constraints to forbid clockwise-oriented circuits. To this end, Fekete et al. (2022) introduce inequalities called Slap Constraints. These inequalities use a straight line with points of Non both sides and no point of Non the line. This line intersects some arcs in any Hamiltonian circuit, and these arcs alternate in direction when crossing the line. The inequalities force the first selected arc to be counterclockwise, the second to be clockwise, etc. To be more precise, let us consider a vertical straight line. Going from bottom to top along this line, the first intersected edge must correspond to an arc in the circuit oriented from left to right, the next from right to left, and so on until we reach the last one, which must be oriented from right to left. To formally write the inequalities, let (i1, j1),...,(ik, jk) be the arcs in Aintersecting the vertical line from left to right. Then the slab inequalities for the fixed straight line are: 0≤ m X l=1 (xiljl−xjlil)≤1 for all m= 1, . . . , k. (6) 9 0 1 2 34 5 6 78 9 Figure 3: subpolygons after drawing all possible edges of instance uniform-0000010-1. for partitioning these triangles into Tand T′. The new model, referred by MT3D, is (1)–(4) and (19)–(26). Note that this combined model does not need (5) to guarantee the non-crossing edges requirement. 4.4. Subpolygon-based Formulation The previous formulations determine a polygon as the union of triangles. Alternatively, a polygon can be seen as the union of other smaller polygons (called subpolygons), not necessarily triangles. A natural candidate set of subpolygons to form a polygonisation of Nis suggested by the graph G(N, A), drawing each arc (i, j)∈Aas a straight line between iand j. Some lines intersect, creating other points not in N, and the result is a mosaic of tiles, each one being a candidate subpolygon to be part of the optimal polygon. The result is a tessellation of the convex hull CH(N) (the subpolygons are the disjoint tiles). Figure 3 shows the 59 subpolygons on instance uniform-0000010-1 with |N|= 10. This section describes a new model based on selecting these subpolygons, and we start with notation. Let Wbe the set of subpolygons generated with the above-described procedure on N. Without abuse of notation, each subpolygon s∈Wis the subset of arcs in A, each arc related to each edge of the subpolygon. As done in the previous section with the triangles, each subpolygon s∈W is counterclockwise oriented, which determines the orientation of the arc related to each edge. Differently to the previous section, now the arc set s may not represent a circuit in G(N, A) since the vertices of the subpolygon smay not be in N. Given s, t ∈W, we say that sand tare adjacent if sand tshare a line segment on the plane, meaning that there exist (i, j)∈A′′ such that 16 (i, j)∈sand (j, i)∈t. In such a case, we say that sis on the left of (i, j) and tis on the right of (i, j). The set of all subpolygons on the left of an arc (i, j)∈Ais denoted by Wij. Each edge {i, j} ∈ E′′ is associated with a subset Hij of pairs {s, t}of adjacent subpolygons with s∈Wij and t∈Wji. For example, the line segment [2,5] in Figure 3 separates three pairs of adjacent subpolygons. Finally, we denote by csthe area of each subpolygon s∈W. We use here again the variables xaused in the previous models to determine the polygon and introduce the binary variable ηsfor each s∈W, assuming value 1 if sis selected to cover the interior area of the polygon, and value 0 otherwise. Now, models for the MinArea and MaxArea are: (1)–(4) and ηs=xij for all (i, j)∈A′and s∈Wij (27) ηs−ηt=xij −xji for all {i, j} ∈ E′′ and {s, t} ∈ Hij : (28) s∈Wij, t ∈Wji xij ≤ηs≤1−xji for all (i, j)∈A′′ and s∈Wij.(29) Equation (27) selects all subpolygons adjacent to a boundary line segment when the tour traverses this line segment. For each arc (i, j) separating the subpolygons sand t, with s∈Wij (i.e. left of (i, j)) and t∈Wji (i.e. right of (i, j)): Equations (28) force ηs=ηtwhen [i, j] is not traversed by the tour, and Constraints (29) imply ηs= 1 and ηt= 0 if xij = 1, and ηs= 0 and ηt= 1 if xji = 1. When the variables xij represent a Hamiltonian circuit, the linear system (27)–(29) has only one solution; thus, the integrability of the ηsvariables is unnecessary. This linear system is related to the one proposed by Dantzig et al. (1985) for the triangulation of a set of points, aiming at covering each point in the convex hull with exactly one triangle. Model (1)–(4) and (27)–(29) is referred by MSD since it uses subpolygonbased variables and the directed graph G(N, A). As in MT1D, Formulation MSD selects objects to cover the internal area of the polygon. One could expect another model, similar to MT2D, by selecting subpolygons to cover the external area. However, it is the same formulation. It is worth observing that Model MSD is an aggregated version of MT3D. Theorem 4. The LP bound from MSD is not better than the one from MT3D. Proof. The result arises after duplicating (27)–(29) and replacing ηswith Pt∈Vsytin one copy and with 1 −Pt∈Vsy′ tin the other copy, where Vs represents the triangles including the subpolygon s. 17 A potential advantage of MSD could be that |W|could be expected to be smaller than |V|, thus a model with a smaller number of variables, and therefore an MILP could be faster in practice. Section 6 shows that this does not occur on our instances in part because Theorem 10 reduces V. 4.5. Edge-based Formulations The previous new formulations replace the connectivity constraints (9) in the seminar triangle-based formulation with the Asymmetric TSP constraints (2)–(4). However, MinArea and MaxArea concern geometric objects in the Euclidean space. Thus, one could expect other formulations with Symmetric TSP constraints instead. This section shows how to adapt the formulations in the previous subsections to use edge-based variables xe instead of arc-based variables xa. The key element is to represent the polygon as a non-oriented tour in G(N, E) rather than a Hamiltonian circuit in G(N, A), and to compute an internal triangulation Tof the polygon. To this end, we make use of the set Vof counterclockwise-oriented triangles and the subsets Vij and Vji for each {i, j} ∈ E. Note that, while Vij ∪Vji are all the triangles with the line segment [i, j] being an edge, each triangle in Vij ∪Vji has a third vertex kthat is on one side of the line segment [i, j]. Precisely, all triangles in Vij have the third vertex on the same side of {i, j}, and all triangles in Vji have the third vertex on the other side. Note also that one of the two subsets is empty when {i, j} ∈ E′. The new model uses the already-introduced binary variable ytfor each t∈V, thus defining the selected triangulation by T={t∈V:yt= 1}. Let us now introduce a new binary variable x′ efor each e∈E, assuming value 1 when eis in the tour and 0 otherwise. We write indistinctly x′ eor x′ ij or x′ ji when e={i, j}. All the formulations in this section use the x′ e variables to determine the edges of the polygon. A Symmetric TSP solution is determined by (for example) the following integer linear system: X e∈δ(i) x′ e= 2 for all i∈N(30) X e∈E(S) x′ e≤ |S|−1 for all S⊂N(31) x′ e∈ {0,1}for all e∈E. (32) The non-crossing edges requirement can be guaranteed with x′ ij +x′ kl ≤1 for all {[i, j],[k, l]} ∈ I2.(33) 18 Now, Model MT1D suggests formulating MinArea and MaxArea by the objective function (7) subject to (30)–(33) and X t∈Vij yt+X t∈Vji yt=x′ ij for all {i, j} ∈ E′(34) −x′ ij ≤X t∈Vij yt−X t∈Vji yt≤x′ ij for all {i, j} ∈ E′′ (35) x′ ij ≤X t∈Vij yt+X t∈Vji yt≤2−x′ ij for all {i, j} ∈ E′′ (36) yt≥0 for all t∈V. (37) Equations (34) force selecting one adjacent triangle to each edge on the boundary of the convex hull of N, and no one otherwise. These equations ensure that the selected triangles are in the internal region of the polygon. Inequalities (35) force that either zero or two triangles are selected if the edge is not in the tour. Inequalities (36) force one adjacent triangle to be selected if the edge is in the tour, and at most two adjacent triangles may be selected otherwise. A proof similar to that of Theorem 2 produces the following result. Theorem 5. When variables x′ erepresent a simple polygonization of N, all the extreme solutions of (34)–(37) are integers. The new formulation (7) and (30)–(37) is referred to by MT1U. In contrast to previous models, the objective function (1) is not an option since variables xaare not present. Model MT2D suggests selecting an external triangulation T′using y′ tbinary variables, so T′={t∈V:y′ t= 1}. The new formulations for MinArea and MaxArea are defined by the objective function min/max c(CH(N)) −X t∈V cty′ t subject to (30)–(33) and X t∈Vij y′ t+X t∈Vji y′ t= 1 −x′ ij for all {i, j} ∈ E′(38) −x′ ij ≤X t∈Vij y′ t−X t∈Vji y′ t≤x′ ij for all {i, j} ∈ E′′ (39) x′ ij ≤X t∈Vij y′ t+X t∈Vji y′ t≤2−x′ ij for all {i, j} ∈ E′′ (40) y′ t≥0 for all t∈V. (41) 19 Theorem 6. When variables x′ erepresent a simple polygonization of N, all the extreme solutions of (38)–(41) are integers. The new formulation is referred to by MT2U. As already observed on models MT1D and MT2D, the system (38)–(40) becomes the system (34)–(36) by replacing Pt∈Vij y′ twith 1 −Pt∈Vij ytfor all (i, j)∈A. Model MT3D suggests formulating MinArea and MaxArea with the constraints (30)–(32) and (34)–(41). This model does not need (33) to force the non-crossing edges requirement. As for the objective function, we can use the one in MT1D or MT2D, or even a convex combination like: min/max 1 2c(CH(N)) + 1 2X t∈V ct(yt−y′ t). The formulation with this objective function is referred to by MT3U. Using the linear equation x′ ij =xij +xji for each {i, j} ∈ Eone gets: Lemma 7. Each solution of (19)–(26) is also a solution of (34)–(41). Theorem 8. The LP bound from MT3U is not better than the one from MT3D. Given a Hamiltonian circuit with non-crossing arcs described by the oriented arc-based variables xafor a∈A, one can check the existence of the internal and external triangulations with the system (34)–(41) instead of the system (19)–(26) (since both systems define the same integer polytope). Empirically, we have observed that MILP solvers manage the system (34)–(41) with more difficulties because it has fewer equations and more inequalities than the system (19)–(26). In addition, when the variables xadetermine a fractional TSP solution, we also observed that the polytope defined by (34)–(41) may include other non-integer extreme points than the ones in (19)–(26). For that reason, there is no gain in using (34)–(41) with oriented variables xa. Instead, the smaller number of non-oriented variables x′ ecould perhaps compensate for the weak LP bound of MT3U. Unfortunately, Section 6 shows that this is not the case in our benchmark instances. Given a Hamiltonian tour with non-crossing edges described by the nonoriented variables x′ e, one could also think of a counterclockwise oriented Hamiltonian circuit to later check the existence of a triangulation with the system (19)–(26). However, in this case, we do not have a linear system defining xij and xji from x′ ij. 20 Finally, Model MSD suggests another edge-based formulation. Let W and Hdefined as in Section 4.4. Now, the MinArea and MaxArea can be modelled by the objective function min/max X s∈W csηs subject to (30)–(32) and ηs=x′ ij for all {i, j} ∈ E′and s∈Wij ∪Wji −x′ ij ≤ηs−ηt≤x′ ij for all {i, j} ∈ E′′ and {s, t} ∈ Hij x′ ij ≤ηs+ηt≤2−x′ ij for all {i, j} ∈ E′′ and {s, t} ∈ Hij. The new formulation is referred to by MSU. Theorem 9. The LP bound from MSU is not better than the one from MSD. 4.6. Reformulation of Model (1)–(6) We close this section with a reformulation of Model (1)–(6) to have it in the same framework as the other models and perform a fair computational comparison in Section 6. Model MT1D includes ytvariables to triangulate inside the circuit. Model MT2D includes y′ tvariables to triangulate outside the circuit. Model MT3D includes both ytand y′ t. Now Model (1)–(6) by Fekete et al. (2022) can be seen as the remaining case, i.e., it is an arc-based model without ytand y′ t. For that reason, we refer to it as Model M0D. Since there is no triangulation in Model M0D, the area of the polygon is forced to be computed with the objective function (1). This objective function needs additional constraints on the xavariables to avoid a counterclockwise-oriented circuit (and make the model valid for MinArea). An alternative to the slab constraints (6) when using the TSP representation (2)–(4) is to exploit the precedence relations between the points on the convex-hull boundary N′. See e.g. Balas et al. (1995) for valid linear inequalities to force a point ito be before a point jin a circuit starting and ending at point r, with i, j, r ∈N′. Section 5.2 shows a simpler procedure when the TSP constraints are modelled with flow variables. 5. Implementation details This section describes the main ingredients of a method to solve the formulations in this paper with a general-purpose MILP solver. 21 5.1. Variable fixing The triangle-based models suffer the inconvenience of a (simple) polygon admitting many triangulations. The number of triangulations of a convex polygon with nvertices is the so-called (n−2)-Catalan number (see e.g. Loera et al. (2010)), which is large. The number of triangulations is smaller for non-convex polygons, but it may still be very large, meaning that trianglebased models may include many variables ytand y′ t(and many alternative solutions representing the same polygon). The following result helps to eliminate some variables. Let A′′ I={(i, j)∈A: (i, j)∈A′′ :i, j ∈N′}and E′′ I={{i, j}: (i, j)∈A′′ I}. Fekete et al. (2022) observe that a polygonisation of Ncannot have an edge on an arc (i, j)∈A′′ Ibecause there are other points on both sides of the line segment determined by iand j. Thus, they fix xa= 0 when a∈A′′ I. The following result extends the fixing also to the ytand y′ tvariables. Theorem 10. A non-convex polygon on the plane admits a triangulation where triangles do not include edges in E′′ I. Proof. For each {i, j} ∈ E′′ I, a triangulation with a triangle in Vij has another triangle in Vji. Let k, l be the other vertices of these triangles. If k∈ N′or l∈ N′, replace the two triangles with one in Vkl and another in Vlk. Otherwise, select the triangle from the triangulation adjacent to one of the previous triangles. The new triangle brings in another vertex q. The procedure repeats until q∈ N′, and then one must replace the considered triangles with ones from Vq. The replacement is possible because all the considered vertices determine a convex subpolygon. As a consequence, the LP relaxations of the models are stronger with: •xa= 0 for all a∈A′′ I, •x′ e= 0 for all e∈E′′ I, •yt=y′ t= 0 for all t∈Vaand a∈A′′ I, •ηs= 1 for all s∈Waand a∈A′′ I. 5.2. TSP representations The new models for MinArea and MaxArea force the connectivity constraints between selected objects using TSP variables and constraints. TSP has been deeply investigated, and there are several alternative MILP formulations. We have used (2)–(4) for the oriented variant and (30)–(32) 22 for the non-oriented variant. Modern TSP solvers deal with large instances managing (3) and (31) through branch-and-cut algorithms. However, since MinArea and MaxArea seem to be more complex problems (30 points is challenging), it makes sense to replace these exponential-size formulations with compact ones. The price to pay for having a polynomial number of constraints is to add a polynomial number of continuous variables. This section shows an alternative. Let rbe the arbitrary point taken from N′. The so-called “multi-commodity flow formulation” for the TSP uses a variable fk afor each k∈V\ {r}and each a∈A. It takes a positive value if a path along the tour from rto kgoes through arc a. Then, an oriented tour in G(N, A) can be modelled with (2), (4) and X a∈δ−(i) fk a−X a∈δ+(i) fk a=     1 if i=k −1 if i=r 0 otherwise for all k∈N\ {r}(42) 0≤fk a≤xafor all k∈N\ {r}and a∈A provides LP bounds equal to the TSP formulation using (3). In addition, the flow variables also allow modelling the counterclockwise orientation with X a∈δ+(i) fj a= 1 for all (i, j)∈A′:i, j ∈N′\ {r} instead of the Slab Constraints (6). Similarly, a non-oriented tour can be modelled with (30), (32), (42) and fk ij ≥0, fk ji ≥0, fk ij +fk ji =x′ ij/2 for all k∈N\ {r}and {i, j} ∈ E instead of the so-called Subtour Elimination Constraints (31). 5.3. Non-crossing edges requirement Covering the internal and external areas of a polygon by triangles or subpolygons guarantees the non-crossing edges requirements, i.e. the polygon generated from models MT3D,MT3U and MSD is simple. Hence, Inequalities (5) and (33) are unnecessary to have valid formulations. Still, adding them may strengthen the LP relaxation, especially when they are lifted as we now describe. Given {[i, j],[k, l]} ∈ I2, Inequality (5) is dominated by X t∈Vij yt+xji +X t∈Vkl yt+xlk ≤1 23 and X t∈Vij y′ t+xij +X t∈Vkl y′ t+xkl ≤1. Based on our computational experiments, they strengthen the LP relaxation of MT1D and MT2D, respectively. Similarly, Inequality (33) suggests X t∈Vij (yt+y′ t) + X t∈Vkl (yt+y′ t)≤1, which may strengthen the LP relaxation of MT3U. It is known in the MILP literature that the above inequalities create incompatibilities between variables, and the LP relaxations can then be strengthened with further linear inequalities. In this direction, Fekete et al. (2022) refer to maximal clique inequalities and halfspace constraints, which are heuristically generated within a branch-and-cut framework. Modern MILP solvers already implement internal functions to generate similar inequalities when incompatibilities are detected in the given formulation. Thus, we do not extend this article in this direction. 5.4. Other valid linear constraints Equalities (11), Pt∈Ty′ t=|N|−|N′|and Pt∈Tct(yt+y′ t) = c(CH(N)), and the inequalities X t∈Vij (yt+y′ t)≤1 for all (i, j)∈A strengthen LP relaxations, according to our computational experiments. 5.5. Alternative objective-function coefficients When the Cartesian coordinates of iand jare (αi, βi) and (αj, βj), respectively, the coefficient cij in the objective function (1) may be fixed to (αiβj−αjβi)/2. Other alternative definitions of cij to compute the polygon’s area are possible, like (αi+αj)(βj−βi)/2 or (βi+βj)(αi−αj)/2. In all cases, cij =−cji, which makes the coefficient matrix highly asymmetric. As observed in Section 3, when using a formulation with only the arcbased variables xa, additional constraints are necessary to force the counterclockwise orientation of the circuit. These constraints are unnecessary for the new arc-based models in Section 4 because the counterclockwise orientation is in each triangle or subpolygon, and then the linear systems extend this orientation to the circuit. This simple observation implies that, for 24 example, Models MT1D and MT3D can be reformulated using the objective function (7) instead of (1). Computational experiments showed better performance when using (7). This choice does not apply to Model M0D because variables ytand y′ tare absent, hence (1) is mandatory. There is no choice on the edge-based models either because variables xaare absent. All formulations compute the area of a polygon by adding the areas of pieces (triangles or subpolygons). The literature on Computational Geometry has different alternatives to compute these areas. We now show one of these procedures, known as Pick’s Theorem (see e.g. Mark Berg and Schwarzkopf (1997)), which we found particularly interesting for the subpolygon-based models. The procedure assumes that each vertex in Nis a grid point (i.e. a point with integer coordinates) on the Cartesian plane, and it determines the area of a polygonisation Pof Nby counting the grid points in the interior of Pand on the boundary of P, as follows: Theorem 11. Let Pbe a polygon whose vertices have integer coordinates on the Cartesian plane. Then the area of Pis the sum of the number of grid points in the interior of P, plus half the number of grid points on the boundary of P, minus one. The above result gives an alternative evaluation of the area ctof a triangle t∈Tand csof a subpolygon s∈W. Recall that we assume that the coordinates of Nare integers. Then, the area of a triangle is already a half-integer. However, the vertices of a subpolygon may likely have noninteger coordinates and, therefore, the area may be a fractional number. Instead, with the above result, this number is replaced by a half-integer in the objective function. The replacement reduces numerical issues and benefits modern solvers. Computational experiments showed better performance of MSD and MSU when using this approach. 5.6. Integrality constraints Theorem 2 guarantees that the variables ytassume integer values when the variables xaare forced to be binary. Conversely, the variables xaassume integer values when ytis required to be a binary variable. Computational experiments showed us that the performance is better in the former setting. 5.7. Initial heuristic An initial feasible solution may help an MILP solver on a model. We propose a simple constructive procedure starting from the convex hull CH(N), improved with a 2-edge-exchange search. The constructive heuristic starts 25