 # Packing of Circles on Square Flat Torus as Global Optimization of Mixed Integer Nonlinear problem

The article demonstrates rather general approach to problems of discrete geometry: treat them as global optimization problems to be solved by one of general purpose solver implementing branch-and-bound algorithm. This approach may be used for various types of problems, i.e. Tammes problems, Thomson problems, search of minimal potential energy of micro-clusters, etc. Here we consider a problem of densest packing of equal circles in special geometrical object, so called square flat torus R^2/Z^2 with the induced metric. It is formulated as Mixed-Integer Nonlinear Problem with linear and non-convex quadratic constraints. The open-source solver SCIP, http://scip.zib.de, and its parallel implementation ParaSCIP, http://ug.zib.de, had been used in computing experiments. The main result is "computer aided" proof of the conjecture on optimal packing for N=9 that was published in 2012. To do that, ParaSCIP took about 2059 CPU*hours (16 hours x 128 CPUs) of cluster HPC4/HPC5, National Research Centre "Kurchatov Institute".

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Densest packing problems appear in many areas of discrete geometry. Hereinafter one of these problems is considered: a kind of Tammes problem for so called square flat torus. Flat torus in a nutshell is a factor–space with metric induced by “ordinary” Euclidean metric ( - is an integer lattice in ). The problem of optimal packings of congruent circles into flat torus has been studied in [1, 2]

. As to practical reasons this problem relates to the problem of “super resolution of images” for aero-photography and space imagery.

The problem has very simple formulation: find arrangement of points in flat torus to maximize minimal distance between any pair of these points. The feature of the problem is an “unusual” distance between points of a torus, which are treated as equivalence classes of factor–space . Figure 1 illustrates definition of distance between points and on Square Flat Torus - the length of red segment, i.e. we have the following formula:

 d(x,y)≐√(min{|x1−y1|,1−|x1−y1|})2+(min{|x2−y2|,1−|x2−y2|})2. (1) Figure 1: Distance between x and y on Square Flat Torus

Oleg Musin and Anton Nikitenko have studied this problem for and pool results in the articles [1, 2] The proof of optimality presented there is based on a computer enumeration of irreducible contact graphs corresponding to potentially optimal arrangement of points on the flat torus. They have proved optimal packings for up to 8 and presented a conjecture of optimal arrangement for .

Numerical proof presented hereinafter is based on more general global optimization approach. The packing problem is formulated as mixed-integer nonlinear programming problem (MINLP) with binary variables and non-convex quadratic constraints. This approach is not a new one and has been widely used in studies of packing problems, e.g. see

. We have paid attention to flat torus packing problem in our experiments on global and discrete optimization in distributed computing environment [4, 5, 6]. The main advantage of “geometric” approach is an analytical form of optimal, “max-min”, distance. The main advantage of global optimization – ability to use general-purpose branch-and-bound solver including parallel implementations for high-performance clusters.

#### Paper Organisation.

This paper is organised as follows. Formulation of Flat Torus Packing problem as mixed integer nonlinear (bilinear non-convex) mathematical programming problem is done in Section 2. The Section 3 presents results of computing experiments for , including brief description of SCIP solver  that was used for global optimization. The next Section 4 shows results for obtained by ParaSCIP solver  that is parallel implementation of SCIP based on MPI. Some details of ParaSCIP compilation on HPC4/HPC5 cluster of NRC “Kurchatov Institute” are provided also. Finally, we present solution found by ParaSCIP for , which confirmed conjecture made in [1, 2]. Conclusion Section 4 is followed by Acknowledgements.

## 2 Formulation as global optimization MINLP

Let be a set of unordered pairs of indices . The problem may be formulated as follows:

 D→maxxik:D⩽∑k=1:2(min{∣∣xik−xjk∣∣,1−∣∣xik−xjk∣∣})2  ((i,j)∈IJ),0⩽xik⩽1(k=1:2, i=1:N). (2)

Formulation (2) has non-differentiable functions in constraints. Let’s avoid non-smoothness at the expense of introducing auxiliary continuous and binary variables:222There are a number of literature on various (similar to each other) “tricks” to perform such conversion of non-smooth problems to MILP or MINLP, but, it seems that the article  was one of the first.:

 (3)

Take the first equation of (3). It is equivalent to the following system of inequations ():

 yijk⩽∣∣xik−xjk∣∣,yijk⩽1−∣∣xik−xjk∣∣,yijk⩾∣∣xik−xjk∣∣−ηijkyijk⩾1−∣∣xik−xjk∣∣−1+ηijk=−∣∣xik−xjk∣∣+ηijk. (4)

Equivalence means that satisfies first equation of (3) iff there exists some binary that satisfies (4) with the same . Easy proof may be done considering that 1 is the maximal difference between functions and on the interval (inclusion follows from inclusions and ), see Figure 2.

Continue transformation of (4). Note that the second and the third inequations are equivalent to the following:

 xik−xjk⩽1−yijk, xjk−xik⩽1−yijk,xik−xjk⩽yijk+ηijk, xjk−xik⩽yijk+ηijk.or as two side inequations −yijk−ηijk⩽xik−xjk⩽1−yijk,−1+yijk⩽xik−xjk⩽yijk+ηijk (5)

Let’s use definition of in (3) to transform the first and the last lines in 4:

 zijk⩽−yijk,yijk⩾zijk+ηijk.or as two side inequations zijk+ηijk⩽yijk⩽−zijk (6)

Note that definition of is equivalent to the following system of linear inequations:

 zijk⩽xik−xjk, zijk⩽xjk−xik,zijk⩾xik−xjk−2ζijk,zijk⩾xjk−xik−2(1−ζijk).or as two side inequations zijk⩽xik−xjk⩽zijk+2ζijk,−zijk−2(1−ζijk)⩽xik−xjk⩽−zijk. (7)

The proof of that equivalence is the same as mentioned after system (4) considering that 2 is the maximal difference between functions and on the interval , see Figure 3.

Finally, from definitions (3) and relations (4)–(7) it follows that the problem (2) is equivalent to the following mixed-integer non-linear (and non-convex) problem with quadratic constraints (, , ):

 D→max(with variables xik,yijk,zijk,ηijk,ζijk),s.t.:D⩽∑k=1:2y2ijk,−yijk−ηijk⩽xik−xjk⩽1−yijk,−1+yijk⩽xik−xjk⩽yijk+ηijk,zijk+ηijk⩽yijk⩽−zijk,zijk⩽xik−xjk⩽zijk+2ζijk,−zijk−2(1−ζijk)⩽xik−xjk⩽−zijk,0⩽xik⩽1,yijk∈R,zijk∈R,ηijk∈{0,1},ζijk∈{0,1}. (8)

The problem has continuous variables and binary variables . Then, except the last line, the problem has linear constraints with continuous variables, just as much linear mixed-integer constraints with continuous and binary variables and quadratic non-convex constraints with continuous variables.

In addition to above constraints, some auxiliary constraints may be added to reduce the number of redundant solutions that might be obtained by translation by axis “OX”, “OY”, renumbering of points, mirror image etc.:

 x11=0.5,x12=0 (the first point is fixed to% (0.5,0)),x(i+1)2⩽xi2 (1⩽i⩽N−1) (ascending 2nd coordinates),x21⩽x11. (9)

Simplest additional constraint, in the last line of (9), reduces computing time almost twice (due to twice less volume of domain in multi-dimension space of continuous variables).

## 3 Computing experiments with SCIP, N⩽8

Computing experiments were performed with SCIP (Solving Constrained Integer Programming), . This is a fairly popular an open-source solver, which can be used freely for research and educational purposes. On the home page scip.zib.de one can read: “SCIP is a framework for Constraint Integer Programming oriented towards the needs of mathematical programming experts … as a pure MIP and MINLP solver or as a framework for branch-cut-and-price. SCIP is implemented as C callable library and provides C++ wrapper classes for user plugins. It can also be used as a standalone program to solve mixed integer programs given in various formats such as MPS, LP, flatzinc, CNF, OPB, WBO, PIP, etc.”

In our studies on optimization modelling we prefer another, so called NL-format (representing an instance of mathematical programming problem) from AMPL (A Mathematical Programming Language) , which has almost 35 years long history (since 1985). SCIP supports NL-format by special scipampl build. Originally, AMPL required usage of special commercially licensed translator: to create NL-file that might be passed to any AMPL-compatible solvers, including free ones; to parse solution SOL-files returned from solvers. But in 2005 AMPL developers disclosed internal formats of NL-files . Thanks to that, Pyomo (PYthon Optimization Modeling Objects, pyomo.org, , open-source and free optimization modeling tool) now supports creation of NL-files. Thus, very popular in scientific researches Python programming language may be used to generate optimization problems, which may be solved by a proper AMPL-compatible solver.

Important feature of SCIP is its capability to solve optimization problems having polynomials in constraints (other non-linearities are admitted also). In our experiments, the Thomson problem, which has been formulated as NLP with polynomials of the 4th degree in equality-constraints, was solved by SCIP . Details of implementations of branch-and-bound algorithm in SCIP for the case of bilinear and non-convex polynomial constraints may be found in the article . For brevity we give the following citation: “SCIP uses convex envelopes for well-known univariate functions, linearization cuts for convex constraints, and the classical McCormick relaxation of bilinear terms. All of these are dynamically separated, for pure NLPs also at a solution of the NLP relaxation…”.

Returning to Flat Torus Packing Problems, Pyomo is used to create NL-file from MINLP presented in (8) and (9). Then NL-file is processed by scipampl application, which returns SOL-file with solution and some LOG-file with auxiliary information about solving process. Finally, SOL-file is processed by Python code via Pyomo package features to analyse solution obtained (including creation of all illustrations presented below).

Solution times for cases are presented in Table 1. Problems with has been solved on desktop [CPU=Intel Core i7-6700 @ 3.40GHz, Mem=32Gb]. For the problem has been solved on standalone server [CPU=2Xeon5620 @ 2.4Ghz, Mem=32Gb]. SCIP had been run with default settings, except: relative gap had been set to and memory limit - to 28Gb (actually the worst case occupied about 8Gb).

#### Results for N=7.

The case deserves more attention as, see [1, 2], there are three different (up to isometric transformation) optimal arrangements, see , Fig.3b–Fig.3d. A few efforts have been done to find all these configurations in results of SCIP solving.

By default SCIP solver stores optimal solutions in a list that is available by proper commands of SCIP-console. So, after successful completion of solving user can compare a number of solutions. Results are presented in the figures 46. They have been selected manually (“draw-and-compare”) from 11 optimal solutions found by standalone SCIP process. Because branch-and-bound algorithm inherently differs from “exact” combinatorial method (enumeration of irreducible contact graphs) that has been used in [1, 2], SCIP founds redundant solutions (auxiliary constraints (9) are not enough to avoid them).

Every configurations coincides with one of those found in  after some isometric transformation (see captions of the figures 46). Pay attention to a free position of the point number 2 on the Fig.4, its circle can be freely moved within area surrounded by other grey circles. Figure 4: N=7, the 1st configuration (see , Fig.3b, d∗=11+√3) Figure 5: N=7, the 2nd configuration (see , Fig.3c, rotate 90o↻ , d∗=11+√3) Figure 6: N=7, the 3d configuration (see , Fig.3d, flip ↕ and rotate 90o↺ , d∗=11+√3) Figure 7: Optimal configuration for N=8 found by SCIP (see , Fig.3e, d∗=11+√3)

## 4 Computing experiments with ParaSCIP, N=8,9.

The Table 1 demonstrates dramatic growth of solving time (and complexity of Flat Torus Packing Problem) on increasing value of . May be it is the reason for that the article  presents only conjecture about optimal configuration for . Our attempts to solve problem for by “single–threaded” SCIP have failed. Computing efficiency of branch-and-bound algorithm implemented in SCIP may be substantially increased by its parallel implementation as ParaSCIP solver does.

#### ParaSCIP — fine-grained parallelization implementation of B&B.

ParaSCIP, ug.zib.de, , is a distributed memory massively parallel MIP and MINLP solver based Ubiquity Generator (UG) framework. This framework is aimed at making MIP and other discrete problem solvers parallel. In the framework the solver and the communication mechanism are abstracted making it easier to parallelize different solvers. Two communication mechanisms are available in the library: distributed memory MPI and shared memory POSIX-threads. SCIP, CPLEX and Xpress solvers are supported with MPI and only SCIP with both mechanisms. Only SCIP-based specializations of the framework (ParaSCIP for MPI-based and FiberSCIP for POSIX-threads based implementations) are publicly available. ParaSCIP can utilize quite big computing resources which allows it to solve problems not solvable even with commercial solvers on a single host: ParaSCIP ran on 80000 cores in parallel while solving open MIP instances from MIPLIB2010 .

Compiling SCIP and ParaSCIP on HPC4/HPC5 (cluster of National Research Center ”Kurchatov Institute”, www.top500.org/site/50615) is not easy because HPC4/HPC5 has CentOS 6 with GCC 4.4 tool-chain installed which does not support C++11 extensions required for SCIP build. We used another machine with a similar CentOS version but with devtoolset-7 (GCC 7.3) installed to build solvers and then copied them to the cluster. Machines had Intel CPUs of different families (Haswell on cluster vs westmere on our machine) so default optimizations during compilation could be not optimal for the cluster. Finally Ipopt, SCIP and ParaSCIP were compiled on our machine and then copied to the cluster where they were running well.

After consultation with ParaSCIP developers we ran parascip with non-default settings. We used the following ParaSCIP settings:

Quiet = FALSE
LogSolvingStatusFilePath = "./logs/"
LogNodesTransferFilePath = "./logs/"
CheckpointFilePath = "./logs/"
SolutionFilePath = "./logs/"
LogSolvingStatus = TRUE
Checkpoint = TRUE
CheckpointInterval = 1800


Also settings for depth-first search were set for node solving to reduce memory consumption in solvers:

nodeselection/dfs/stdpriority = 300000


The following command line was used to run ParaSCIP:

mpirun parascip parascip.set tammesTorus_d2_p9.cip -q -s dfs.set


Here parascip.set is the ParaSCIP settings file listed above; tammesTorus_d2_p9.cip is the problem definition converted to CIP-format from NL-format; dfs.set is the settings file for depth-first search specified above.

#### ParaSCIP vs SCIP for N=8.

This case have been used to compare performance of ParaSCIP and SCIP on HPC4/HPC5. Solving times are presented in the Table 2. One can pay attention that solving with a single SCIP process on the cluster took much more time than that on a standalone server (see last column of the Table 1). The reason is that the problem instance passed to the cluster did not have the last of auxiliary constraints (9).

To evaluate efficiency of parallelization one should remember one feature of ParaSCIP: one of ParaSCIP processes (dedicated for MPI application), plays role of Load Coordinator and, actually, does not work with branch-and-bound algorithm’s search tree. So we have:
efficiency (CPU): 780/126/8 = 0.77
efficiency (solvers): 780/126/7 = 0.88.

#### Results of ParaSCIP for N=9.

ParaSCIP with 128 processes on 8 nodes (127 solvers) of cluster HPC4/HPC5 had solved the problem in 956 minutes.

It is the most promising result presented in this work. Optimal configuration found is shown in the Fig. 8 and coincides with conjecture presented in [1, 2]. Taking into account load balancing between working processes the following evaluation of complexity may be done: . Total complexity, including CPU dedicated for Load Coordination is . We believe that this complexity may be reduced almost twice if the last of auxiliary constraints (9) would be added. Figure 8: Optimal configuration for N=9 found by ParaSCIP (see , Fig.3f, d∗=1√5+2√3)

## 5 Conclusion

Obtained results confirmed the relevance of global optimization approach in solving hard problems of combinatorial geometry. A few tricks in formulation Flat Torus Packing problem as mixed integer nonlinear problem and “general purpose” solver SCIP with its parallel version ParaSCIP let to give “computer aided” proof of optimal arrangement for 9 points. All results obtained before for have been reproduced also.

Once more has been confirmed semi-empirical rule (see ) that simple auxiliary constraints that reduces the volume of feasible domain (i.e. last inequation of (9)) may substantially reduce solving time for branch-and-bound algorithm using McCormik envelopes for global optimization of nonlinear problems.

General purpose solver SCIP and ParaSCIP, its parallel implementation, can be used successfully for global optimization of mixed integer nonlinear problems with bilinear constraints.

### Acknowledgements

Authors are grateful to Alexey Tarasov for recommendation to try Flat Torus Packing problem in our experiments with distributed implementations of B&B. Authors thank Stefan Vigerske and Yuji Shinano for their consultations on using of SCIP and ParaSCIP solvers.
This work is supported by the Russian Science Foundation (project No. 16-11-10352).
This work has been carried out using computing resources of the federal collective usage centre Complex for Simulation and Data Processing for Mega-science Facilities at NRC ”Kurchatov Institute” (ministry subvention under agreement RFMEFI62117X0016), ckp.nrcki.ru.

## References

•  O. R. Musin and A. V. Nikitenko. Optimal packings of congruent circles on a square flat torus. ArXiv e-prints, dec 2012.
•  O. R. Musin and A. V. Nikitenko. Optimal packings of congruent circles on a square flat torus. Discrete & Computational Geometry, 55(1):1–20, 2016.
•  Ignacio Castillo, Frank J Kampas, and János D Pintér. Solving circle packing problems by global optimization: numerical results and industrial applications. European Journal of Operational Research, 191(3):786–802, 2008.
•  V. Voloshinov, S. Smirnov, and O. Sukhoroslov. Implementation and use of coarse-grained parallel branch-and-bound in Everest distributed environment. Procedia Computer Science, 108:1532–1541, 2017.
•  S. Smirnov and V. Voloshinov. Implementation of concurrent parallelization of branch-and-bound algorithm in Everest distributed environment. Procedia Computer Science, 119:83–89, 2017.
•  S. Smirnov and V. Voloshinov. On domain decomposition strategies to parallelize branch-and-bound method for global optimization in Everest distributed environment. Procedia Computer Science, 136:128–135, 2018.
•  A. Gleixner, M. Bastubbe, L. Eifler, T. Gally, G. Gamrath, R. L. Gottwald, G. Hendel, C. Hojny, T. Koch, M. E. Lübbecke, S. J. Maher, M. Miltenberger, B. Müller, M. E. Pfetsch, C. Puchert, D. Rehfeldt, F. Schlösser, C. Schubert, F. Serrano, Y. Shinano, J. M. Viernickel, M. Walter, F. Wegscheider, J. T. Witt, and J. Witzig. The SCIP Optimization Suite 6.0. Technical Report 18-26, ZIB, Takustr. 7, 14195 Berlin, 2018.
•  Yuji Shinano, Tobias Achterberg, Timo Berthold, Stefan Heinz, and Thorsten Koch. ParaSCIP: a parallel extension of SCIP. In Competence in High Performance Computing 2010, pages 135–148. Springer, 2011.
•  George B. Dantzig.

On the significance of solving linear programming problems with some integer variables.

Econometrica, Journal of the Econometric Society, pages 30–44, 1960.
•  R. Fourer, D.M. Gay, and B.W. Kernighan. AMPL: A Modeling Language for Mathematical Programming. Second edition. Duxbury Press/Brooks/Cole Publishing Company, 2003.
•  D. M. Gay. Writing .nl files. Technical report, No. 2005-7907P. Sandia National Laboratories, 2005.
•  W. E. Hart, C. D. Laird, J.-P. Watson, D. L. Woodruff, G. A. Hackebeil, B. L. Nicholson, and J. D. Siirola. Pyomo–optimization modeling in Python. Second edition, volume 67. Springer, 2017.
•  Stefan Vigerske and Ambros Gleixner. SCIP: Global optimization of mixed-integer nonlinear programs in a branch-and-cut framework. Optimization Methods and Software, pages 1–31, 2017.
•  Yuji Shinano, Tobias Achterberg, Timo Berthold, Stefan Heinz, Thorsten Koch, and Michael Winkler. Solving open MIP instances with ParaSCIP on supercomputers using up to 80,000 cores. In Parallel and Distributed Processing Symposium, 2016 IEEE International, pages 770–779. IEEE, 2016.