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 aerophotography 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:
(1) 
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 mixedinteger nonlinear programming problem (MINLP) with binary variables and nonconvex quadratic constraints. This approach is not a new one and has been widely used in studies of packing problems, e.g. see
[3]. 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, “maxmin”, distance. The main advantage of global optimization – ability to use generalpurpose branchandbound solver including parallel implementations for highperformance clusters.Paper Organisation.
This paper is organised as follows. Formulation of Flat Torus Packing problem as mixed integer nonlinear (bilinear nonconvex) mathematical programming problem is done in Section 2. The Section 3 presents results of computing experiments for , including brief description of SCIP solver [7] that was used for global optimization. The next Section 4 shows results for obtained by ParaSCIP solver [8] 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:
(2) 
Formulation (2) has nondifferentiable functions in constraints. Let’s avoid nonsmoothness at the expense of introducing auxiliary continuous and binary variables:^{2}^{2}2There are a number of literature on various (similar to each other) “tricks” to perform such conversion of nonsmooth problems to MILP or MINLP, but, it seems that the article [9] was one of the first.:
(3) 
Take the first equation of (3). It is equivalent to the following system of inequations ():
(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:
(5) 
Note that definition of is equivalent to the following system of linear inequations:
(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 mixedinteger nonlinear (and nonconvex) problem with quadratic constraints (, , ):
(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 mixedinteger constraints with continuous and binary variables and quadratic nonconvex 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.:
(9) 
Simplest additional constraint, in the last line of (9), reduces computing time almost twice (due to twice less volume of domain in multidimension space of continuous variables).
3 Computing experiments with SCIP, N8
Computing experiments were performed with SCIP (Solving Constrained Integer Programming), [7]. This is a fairly popular an opensource 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 branchcutandprice. 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 NLformat (representing an instance of mathematical programming problem) from AMPL (A Mathematical Programming Language) [10], which has almost 35 years long history (since 1985). SCIP supports NLformat by special scipampl build. Originally, AMPL required usage of special commercially licensed translator: to create NLfile that might be passed to any AMPLcompatible solvers, including free ones; to parse solution SOLfiles returned from solvers. But in 2005 AMPL developers disclosed internal formats of NLfiles [11]. Thanks to that, Pyomo (PYthon Optimization Modeling Objects, pyomo.org, [12], opensource and free optimization modeling tool) now supports creation of NLfiles. Thus, very popular in scientific researches Python programming language may be used to generate optimization problems, which may be solved by a proper AMPLcompatible solver.
Important feature of SCIP is its capability to solve optimization problems having polynomials in constraints (other nonlinearities are admitted also). In our experiments, the Thomson problem, which has been formulated as NLP with polynomials of the 4th degree in equalityconstraints, was solved by SCIP [6]. Details of implementations of branchandbound algorithm in SCIP for the case of bilinear and nonconvex polynomial constraints may be found in the article [13]. For brevity we give the following citation: “SCIP uses convex envelopes for wellknown 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 NLfile from MINLP presented in (8) and (9). Then NLfile is processed by scipampl application, which returns SOLfile with solution and some LOGfile with auxiliary information about solving process. Finally, SOLfile 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 i76700 @ 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).
N  4  5  6  7  8 

Solving time, sec  3  30  118  2552  27240 (454 min) 
Results for N=7.
The case deserves more attention as, see [1, 2], there are three different (up to isometric transformation) optimal arrangements, see [2], 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 SCIPconsole. So, after successful completion of solving user can compare a number of solutions. Results are presented in the figures 4–6. They have been selected manually (“drawandcompare”) from 11 optimal solutions found by standalone SCIP process. Because branchandbound 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).
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 [2] presents only conjecture about optimal configuration for . Our attempts to solve problem for by “single–threaded” SCIP have failed. Computing efficiency of branchandbound algorithm implemented in SCIP may be substantially increased by its parallel implementation as ParaSCIP solver does.
ParaSCIP — finegrained parallelization implementation of B&B.
ParaSCIP, ug.zib.de, [8], 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 POSIXthreads. SCIP, CPLEX and Xpress solvers are supported with MPI and only SCIP with both mechanisms. Only SCIPbased specializations of the framework (ParaSCIP for MPIbased and FiberSCIP for POSIXthreads 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 [14].
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 toolchain installed which does not support C++11 extensions required for SCIP build. We used another machine with a similar CentOS version but with devtoolset7 (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 nondefault settings. We used the following ParaSCIP settings:
Quiet = FALSE LogSolvingStatusFilePath = "./logs/" LogNodesTransferFilePath = "./logs/" CheckpointFilePath = "./logs/" SolutionFilePath = "./logs/" NotificationInterval = 10 LogSolvingStatus = TRUE Checkpoint = TRUE CheckpointInterval = 1800
Also settings for depthfirst 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 CIPformat from NLformat; dfs.set
is the settings file for depthfirst 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).
HPC4/HPC5, CPU Xeon E52680v3 12C @ 2.5GHz 
cores 
solving time, min 

SCIP  1  780 
ParaSCIP  8 (7 solvers)  126 
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 branchandbound 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.
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 semiempirical rule (see [6]) that simple auxiliary constraints that reduces the volume of feasible domain (i.e. last inequation of (9)) may substantially reduce solving time for branchandbound 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. 161110352).
This work has been carried out using computing resources of the
federal collective usage centre Complex for Simulation and Data
Processing for Megascience Facilities at NRC ”Kurchatov Institute”
(ministry subvention under agreement RFMEFI62117X0016), ckp.nrcki.ru.
References
 [1] O. R. Musin and A. V. Nikitenko. Optimal packings of congruent circles on a square flat torus. ArXiv eprints, dec 2012.
 [2] 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.
 [3] 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.
 [4] V. Voloshinov, S. Smirnov, and O. Sukhoroslov. Implementation and use of coarsegrained parallel branchandbound in Everest distributed environment. Procedia Computer Science, 108:1532–1541, 2017.
 [5] S. Smirnov and V. Voloshinov. Implementation of concurrent parallelization of branchandbound algorithm in Everest distributed environment. Procedia Computer Science, 119:83–89, 2017.
 [6] S. Smirnov and V. Voloshinov. On domain decomposition strategies to parallelize branchandbound method for global optimization in Everest distributed environment. Procedia Computer Science, 136:128–135, 2018.
 [7] 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 1826, ZIB, Takustr. 7, 14195 Berlin, 2018.
 [8] 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.

[9]
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.  [10] 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. https://ampl.com/resources/theamplbook.
 [11] D. M. Gay. Writing .nl files. Technical report, No. 20057907P. Sandia National Laboratories, 2005.
 [12] 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.
 [13] Stefan Vigerske and Ambros Gleixner. SCIP: Global optimization of mixedinteger nonlinear programs in a branchandcut framework. Optimization Methods and Software, pages 1–31, 2017.
 [14] 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.
Comments
There are no comments yet.