Annals of Operations Research 104, 213–230, 2001  2001 Kluwer Academic Publishers. Manufactured in The Netherlands.

Globally Optimized Spherical Point Arrangements: Model Variants and Illustrative Results JÁNOS D. PINTÉR ∗

[email protected]

Pintér Consulting Services, Inc. (PCS) & Dalhousie University, PCS: 129 Glenforest Drive, Halifax, NS, Canada B3M 1J2

Abstract. The “well-balanced” distribution of points over the surface of a sphere is of significant interest in various fields of science. The quality of point configurations is typically expressed by criterion functions that have many local optima. A general global optimization framework is suggested to solve such problems. To illustrate the viability of this approach, the model development and solver system LGO is applied to four different model versions. Numerical results – including the visual representation of criterion functions in these models – are presented. The global optimization approach can be tailored to specific problem settings, and it is also applicable to a large variety of other model forms. Keywords: spherical point arrangements, multiextremal criterion functions, global optimization, LGO model development and solver system, illustrative results AMS subject classification: 65K30, 90C05, 90C31

1.

Introduction

We shall consider the problem of “optimal” point distributions on the surface of the unit sphere in R3 . The quality of such point arrangements will be expressed by a given, context-specific criterion function. This verbal problem statement covers many specific cases, and it can also be extended into various directions. For instance, one could consider arbitrary dimensionality of the embedding space Rd (d
3), other feasible sets (surfaces as well as fulldimensional bodies in Rd ), added constraints regarding feasible point configurations, and various measures of configuration quality. The basic problem statement and its versions referred to above are not only of pure theoretical interest: they are also closely related to important problems in applied mathematics, physics, chemistry, biology, and engineering. A few examples of related fields of application are computational complexity, numerical approximation, packing problems, multi-body potential energy models, celestial mechanics, electrostatics, crystallography, molecular structure modeling, and viral morphology. For related expositions consult, for instance, Conway and Sloan [3], Greengard [6], Pain [17], Dean [4], Pardalos, Shalloway and Xue [20], Neumaier [14,15], Hartke [11], Wales and Scheraga [32]. In the ∗ http://is.dal.ca/∼jdpinter/.

214

PINTÉR

context of spherical point arrangements, Saff and Kuijlaars [30] present a concise review of the theoretical background for several model versions. Optimal point arrangement problems are – as a rule – very difficult to solve in the strict mathematical sense, or even numerically. The essential reason for this is the massive multiextremality of the criterion functions. In general, one can expect the number of local optima to be some exponential function of the number of points: see, e.g., the review of Pardalos, Shalloway and Xue [19] which also provides an extensive list of related references. Apparently, the issue of finding suitable strategies to produce globally optimized configurations has been a subject of intensive research for decades. For instance, Pardalos [18] states one of the model variants as an open problem. Smale [31] lists 18 prominent mathematical problems for the next century: two of these – problems 5 and 6 – are closely related to the present subject. (In fact, problem 6 is one of the models discussed later on.) The purpose of our work is to illustrate the relevance and applicability of global optimization (GO) to this important research area. Four spherical arrangement models, reviewed by Saff and Kuijlaars, will be considered. We shall briefly describe a general GO based solution approach. Illustrative numerical results – obtained by using the LGO modeling and solver system – are presented for relatively small, yet non-trivial model instances. The visualization capability of LGO clearly indicates the numerical challenge inherent in point arrangement problems. The paper is concluded by pointing towards diverse problem specifications and model extensions, still within the scope of the general GO methodology proposed. 2.

Problem statement and model versions

The general model of optimal spherical point distributions is the following. Given the unit sphere B in R3 , and a positive integer n, we want to find an n-tuple of unit vectors x(n) = {x1 , . . . , xn }, xi = (xi1 , xi2 , xi3 ), xi  = 1 ( ·  denotes the Euclidean norm)

(1)

distributed on the surface S of B such that x(n) optimizes the value of a given, contextspecific criterion function. Let us denote by dj k = d(xj , xk ) = xj − xk  the distance between points xj and xk. . We shall consider all such pair-wise distances dj k , 1  j < k  n, when defining the quality of point arrangements x(n) by the subsequent formulae (2)–(5). Obviously, the optimization is related to all possible point configurations over S. Four frequently considered model versions will be studied: these are distinguished by their respective objective functions given below. Elliptic Fekete problem.
max dj k subject to

x(n) ∈ S (n)

(i.e., xi ∈ S, i = 1, . . . , n).

(2)

GLOBALLY OPTIMIZED ARRANGEMENTS

215

For reasons of better numerical tractability, one can apply a logarithmic transformation. Hence, the objective in (2) is replaced by the equivalent logarithmic potential function  max log(dj k ) subject to x(n) ∈ S (n) (here log denotes 10-based logarithm). (2 ) Fekete problem. min



dj−1 k

subject to x(n) ∈ S (n) .

(3)

Note that this objective function corresponds to the classical Coulomb potential model. Power sum problem.  max djak subject to x(n) ∈ S (n)

(a > 0 is a given model parameter).

(4)

Saff and Kuijlaars in their cited paper remark that (4) is of interest only in the case, when a < 2; for even n and a
2, optimal configurations consist of two equal groups of coinciding points located at the antithetic points of an arbitrary main circle of S. Tammes (Hard-spheres) problem. max min dj k

subject to x(n) ∈ S (n) .

(5)

Notice the relation of this model-type (and its possible extensions) also to more general “hard object” packing problems. One can observe that the criterion functions in (2)–(5) pose various levels of numerical difficulty. Models (2 ) and (3) need the exclusion of points that are “too close” to each other, to avoid poles and implied numerical problems; (4) seems to have the analytically “easiest” objective function; and model (5) seems to be the “hardest”, especially since its objective is non-differentiable. Let us note that numerous other (scientifically meaningful) objective function forms could also be considered. As a rule, such functions possess a high degree of symmetry, induced by equivalent permutations of point arrangements. The ambiguity of solutions can be excluded by imposing, e.g., the lexicographic arrangement of the points in a complete configuration x(n). The problem of multiextremality is another inherent characteristics of such problems. Among others, Saff and Kuijlaars comment on the fact that the criterion functions in (2)–(5) have a very large number of local optima, and that some of these are fairly close to the “true” (most typically unknown, but conjectured) global optimum. The circumstances outlined above imply the genuine relevance of global scope search strategies.

216

3.

PINTÉR

A global optimization approach to spherical point distribution problems

3.1. The continuous GO model The objective of global optimization is to find the “absolutely best” solution of nonlinear decision models, in the presence of multiple local optima. A detailed exposition of the most prominent GO model types and solution approaches can be found, for instance, in Horst and Pardalos [12]. Several dozens of other textbooks and quite a few informative WWW sites are also devoted to this emerging subject; among others, Pintér [24] presents a concise, Web-enabled review. We shall consider a fairly general GO model form defined by the following ingredients: • • • • •

z f (z) D l, u g(z)

decision vector, an element of the real m-space Rm ; continuous objective function; non-empty set of admissible decisions, D is described by explicit, finite bounds of z which define an embedding “box” in Rm ; k-vector of continuous constraint functions (k
0).

Applying the notation above, the continuous global optimization model is stated as min f (z) subject to z ∈ D ⊂ Rm ,   D := x: l  z  u, g(z)  0 .

(6) (7)

For numerical reasons, the following additional requirements are often postulated: • D is a full-dimensional subset (“body”) in Rm ; • the set of globally optimal solutions to (6), (7) is at most countable; • f and g (the latter component-wise) are Lipschitz-continuous functions on [l, u]. Evidently, the point arrangement problems presented above belong to the scope of model (6), (7). The arrangement x(n) corresponds to z; f is one of the criterion functions (2)– (5); the bounds l and u will naturally follow by a suitable model parameterization of the points in x(n); and – in the point arrangement models stated above – the additional constraints g are absent. Consequently, the model variants (2)–(5) are well-posed global optimization problems that can be handled by suitable GO tools. 3.2. GO formulation of spherical point arrangement problems The configurations x(n) defined by (1) will be most naturally represented by using spherical coordinates. Let us denote by e1 , e2 , and e3 the three unit vectors in the Cartesian coordinate system that has its origin at the center of B. For xi ∈ S, let βi denote the angle between xi and its projection onto the plane defined by e1 and e2 . Further, let αi

GLOBALLY OPTIMIZED ARRANGEMENTS

217

denote the angle between this projected segment and e1 . Then the component vectors xi of x(n) are expressed by −π π  βi  . 2 2 Introducing γi = sin βi as part of an equivalent parameterization, one can write   1/2  1/2  xi = cos αi 1 − γi2 , sin αi 1 − γi2 , γi , 0  αi  2π, −1  γi  1. xi = (cos αi cos βi , sin αi cos βi , sin βi ),

0  αi  2π,

(8)

(8 )

The latter form has the advantage over (8) that if both αi and γi are selected uniformly from their respective interval domains, then the points xi will have a uniform distribution over S. This fact will be numerically important, since randomized sample points are (also) used by the solver procedure. The collection of inequality relations in (8 ) for i = 1, . . . , n defines the explicit lower and upper bounds l, u appearing in the GO model (6), (7). The above parameterization implies that each point arrangement x(n) can be described by 2n variables. In order to eliminate the rotational symmetries of otherwise completely equivalent solutions, one can fix three angles in the spherical representation of x(n). We choose α1 = β1 = β2 = 0,

i.e.,

α1 = γ1 = γ2 = 0.

(9)

These relations imply that the unit vector e1 = (1, 0, 0) is always a component of the solution; further, that at least another – namely, the second – vector in the set x(n) belongs to the (e1 , e2 )-plane. This convention reduces the number of unknown parameters in x(n) to m = 2n − 3, corresponding to the dimensionality in the resulting box-constrained GO problem (6), (7). As mentioned above, a lexicographic arrangement of the points could also be postulated, but we shall not consider such additional constraints here. In the next two sections we shall briefly describe the LGO program system which will be applied to solve illustrative examples of models (2)–(5), under the assumptions (8 ), (9). 4.

The LGO integrated model development and solver system

LGO – abbreviating Lipschitz(-Continuous) Global Optimizer – serves to analyze and solve optimization problems of the form (6), (7), under the general structural assumptions listed in section 3. Therefore it is particularly suitable to handle GO problems related to models that are supported by limited – or difficult to use – analytical information. For instance, the model (5) has a non-differentiable objective, while it still belongs to the Lipschitz problem-class. First, we present a summary of the principal LGO features; this will be followed in section 5 by a short description of the application development steps. For the necessary theoretical background, including also some key implementation details, consult Pintér [21]. Additional aspects – related to more recent developments, current LGO features, and detailed user guidance – are discussed by Pintér [22,23].

218

PINTÉR

4.1. Solver options The core of LGO is a suite of robust and efficient nonlinear optimization methods. Currently, the following component algorithms are offered: • • • •

global adaptive partition and search (continuous branch-and-bound method), global adaptive random search, unconstrained local search, applying an exact penalty function approach, constrained local search, based upon sequentially generated local approximations of the constraint functions g and the objective function f .

The built-in global scope algorithms are also equipped with proper statistical tools. The latter generate a lower bound estimate, based upon the search points and corresponding function values sampled in the global phase of LGO. Both global scope algorithms are gradient-free procedures: this fact makes their application easy in many practical modeling situations in which higher order information is difficult (or computationally expensive) to obtain. The global solvers – activated with suitable parameterization – enable a robust and theoretically established approximation of the solution. In most practical cases, numerical efficiency requires the additional use of local search strategies: this is the main reason for using also built-in local solvers. Note at the same time that the exclusive use of local solvers per se may miss the global optimum, unless the problem is known to have a convex structure. In both local solver options, gradient information is numerically approximated. 4.2. Automatic and interactive operational modes The LGO solvers can be activated automatically or interactively. The automatic option only asks the LGO user to select one of the global search modes; this will be then automatically followed by both local searches. The interactive option makes possible to select solvers one by one, followed by the explicit allocation of computational effort – the maximal number of objective and constraint function evaluations – to each solver mode invoked. This approach makes possible to provide globally established optimum estimates, even in case of a limited search effort. The latter may be necessitated, e.g., by computationally intensive function evaluations. Note the relevance of this point with respect to large instances of the models discussed here, since – at least theoretically – the solver needs to consider all pair-wise distances during the evaluation of any of the objective functions (2)–(5). Finally, LGO also offers an automatic local search procedure (that includes both types of local solvers) launched from a user-supplied initial solution. This feature evidently supports the (correct) solution of convex models, and it also enables the fast use of pre-specified expert estimates of the global solution. (Such estimates are often postulated, including the models discussed here.)

GLOBALLY OPTIMIZED ARRANGEMENTS

219

4.3. LGO versions, development environments, and connectivity issues The currently used standard version of LGO has been developed in Fortran 90 (LF90, by Lahey Computer Systems [13]) for personal computers. In addition to immediate LF90 connectivity, Dynamic Link Library (DLL) connection is also supported by LF90 – at the time of writing – with respect to the following development environments: Borland C/C++ and Delphi; Microsoft Visual Basic and Visual C/C++. Connection to other applications via DLLs (and in certain cases, via static link options) is also possible. LGO versions can be run in “plain” DOS sessions, in a DOS “box” of any of the currently used Windows (95/98/2000/NT) versions, as a stand-alone Windows application, or directly from other Windows applications that support external program calls. LGO can also invoke external application programs, in order to evaluate functions, execute subroutines and auxiliary calculations that may contribute to modules of the GO problem analyzed. Let us note here that the “DOS command-line” LGO implementations (available for workstation and mainframe platforms, or for other professional PC Fortran platforms) do not have the interactive graphics capabilities to be discussed later; these provide plain text output to the screen and to files. The description presented in the next section refers to the current menu-driven Windows version, but model development using any other LGO version follows essentially the same guidelines. (Please contact the author for information regarding currently available implementations.) 5.

Model development using LGO

5.1. LGO menu interface In the Windows 95/98/2000/NT implementation, the solver system is embedded under a menu-driven interface, to assist LGO users. Upon launching the program, its front page – including also a logo and license information – appears. From this screen the following menu options can be activated (the corresponding submenu items are also shown): Model Formulation

Model Solution

Result Analysis

Project File Names User Main File User Function File Compile LGO Model External System Call

Input Parameter File Run LGO

Summary Output File Detailed Output File External System Call

Help

Exit

LGO Information Summary Quit Model Development Summary About LGO

Below we shall briefly describe the functionality of these menu options.

220

PINTÉR

5.2. Model formulation In order to apply LGO, the user first defines the optimization problem, according to standard specifications. The Project File Names option prompts the user to give names to the following program items (their functionality is shown in brackets, and will be subsequently discussed): • user main file (provides main application program frame), • user function file (provides model description), • project executable file (combines the user files named above with the LGO system, and subsequently generates result files), • input parameter file (provides runtime information for the project executable file), • summary output file (contains principal information regarding run results), • detailed output file (contains more detailed information regarding run results). Upon launching this menu option, default file names are provided. These names correspond to a test problem file system available to LGO users, as an example to set up their own applications. The defaults can – and, of course, should – be overwritten following the actual application needs, while keeping the sample user files. The sample files are commented in sufficient detail, to assist their usage. Selecting next the User Main File option invokes the application frame program file, under the name chosen in the previous menu option. In standard LGO shipments the simple, but generally available Notepad editor (a Windows accessory) is activated to display and manipulate this file, as well as all other LGO text files discussed later on. (Note that if LGO users prefer to use a different text editor, then this can also be simply arranged upon request.) The main user file opens the input/output files, and calls the LGO solver system. Optional application-dependent user actions can be also included here, before and/or after executing the LGO run. The User Function File option serves to open a file, in order to describe the optimization model that will be submitted to LGO. The objective function and the constraints are defined here, following again a commented sample file. The Compile LGO Model option invokes the generation of object files on the basis of the user source files discussed above; then links these to the LGO object file system and the auxiliary files (modules, libraries, and resource file) needed. The actual compilation obviously depends on the (Lahey Fortran 90 or 95, Borland C/C++ or Delphi, Microsoft C/C++ or Visual Basic) development environment used. If all compilation and linking operations are successful, then control is returned to the main menu level, and the user may proceed to the model solution stage. In case of compilation and/or linking errors, corresponding messages appear on the screen, in a DOS “box” under Windows. The LGO system also indicates the errors: a message prompts the user to return to editing mode. (Obviously, all errors need to be corrected as this point, before proceeding further.) The External System Call option can be used to launch direct calls – by simply typing the name of an executable program (available from the given path), or by invoking

GLOBALLY OPTIMIZED ARRANGEMENTS

221

Windows Explorer – to arbitrary program systems available to the user. This option (found also under the Result Analysis menu item) may be helpful in describing and analysing models, or in importing/exporting information. 5.3. Model solution Following successful model compilation and linking, the LGO executable program is ready to run, but it still needs certain input parameters. These parameters are set by selecting the Input Parameter File option. Again, a commented sample of this file is provided to assist users, thus its preparation is straightforward. In particular, the number of variables and constraints, explicit lower and upper bounds, nominal values (to support stand-alone local search) and a few key optimization parameters are defined in this file. This structure allows for repeated runs on different feasible interval sets, launched from different starting poins, and/or using modified optimization parameters, using the same executable program. The Run LGO menu option launches the application. As mentioned earlier, LGO executable programs can be run in automatic and interactive operational modes. All user actions are prompted by the system, with added brief on-screen explanations. During program execution, the LGO system generates troubleshooting messages (as needed) and terse progress reports – currently applied solver mode, incumbent optimum estimate, and the number of function calls – to the screen. In interactive mode, the user is prompted after completing each solver option to select another one, or to terminate the optimization procedure. Upon successfully terminating the run, two automatically generated output files are available for inspection and further work. (In case of runtime errors caught by LGO or by the operating system – caused typically by modeling flaws – the user needs to return to the model formulation stage.) 5.4. Visualization LGO has built-in visualization capabilities, to assist the model development procedure. Specifically, upon completing the optimization process the user can view projections of the aggregated objective plus constraint-induced penalties (as needed), in interactively selected variable subspaces. In these projections, all other coordinates are held at their optimized estimate. The search progress is also summarized in the figures, by displaying the projected scatterplot of improving search points; all other coordinates correspond to the actually found sequence of improving solution vectors. Let us note here that feasible/infeasible set projections with respect to selected constraint functions – not relevant in the models discussed in this paper – can also be generated. Although the use of this visualization option can be simply skipped, in most cases it is a salient idea to check the graphical information – at least in subspaces of special interest. By visual analysis, the user can verify his/her preliminary concept regarding the shape of the objective and constraints. This can effectively help in finding modified problem formulations and/or LGO input parameters to be applied in subsequent LGO runs. Such pictures may also be used to suggest the reduction of the search region, in

222

PINTÉR

order to enable a more detailed “exploration” around the optimum estimate found. As mentioned earlier, runs on modified box regions, or from various starting points can be launched without recompiling the LGO executable program. Illustrative examples to visualize point arrangement problems will be presented later on. After the LGO optimization run is terminated, control is returned again to the main menu level. Accordingly, repeated model formulation, compilation and linking, input file changes, or result analysis may follow. 5.5. Result analysis The solution procedure is typically followed by an inspection of the text result files, automatically generated by LGO. This includes two options. The Summary Output File option displays a basic report of the results. Specifically, it presents the optimum estimate, (global search phase based) lower bound estimate, solution vector estimate, constraint function values, and the total runtime. The Detailed Output File option displays a more thorough description of the program run. Namely, it echoes the variable ranges, nominal values, and the input parameters of all sequentially invoked solver options. Furthermore, it also keeps track of the improving solutions found in the global search phase(s), and the stopping criteria met during execution, in addition to the information presented in the summary result file. Obviously, these two result files may have different importance at different stages of model development. The summary file is most useful in “routine” LGO runs, while the additional information provided by the detailed output file can assist both model development and the choice of LGO solver parameterizations. The results obtained may be directly communicated to other programs. These external programs are typically available via the user main and function files of LGO, but can also be set up independently; they can also be attained by the External System Call option discussed above. Note that returning to LGO immediately after a completed run will overwrite the currently generated output files. Therefore in sequentially executed runs with changing model forms, input parameters, or runtime operations it is advisable to save the results. (This is simply done by making use of a standard Windows dialog that can be enabled from the text editor used.) After the analysis of results is completed, control is returned again to the main menu level. In a typical application development process, the model formulation, run and analysis stages are applied in an iterative fashion. 5.6. Help In addition to the detailed User’s Guide [23], concise on-line help is available. The LGO Information Summary menu option serves to introduce novice users to the program system. The general mathematical model form encompassed by LGO, main system features, hardware and operating system requirements, notes on installation and usage, connectivity to other development platforms, available LGO versions, and a list of existing applications are summarized in this help file.

GLOBALLY OPTIMIZED ARRANGEMENTS

223

The Model Development Summary menu option provides step-by-step instructions to users, regarding the essential stages of application development using LGO. The About LGO option displays the LGO front screen, including copyright and licensee information. 5.7. Exit Upon selecting the Exit | Quit option, all current project files are closed and saved; all temporary files are deleted; and the LGO run is terminated. Summing up the review of LGO menu options, it can be seen that essential background information – related to the subject of global optimization, LGO, and to the model structure used – is permanently available. Moreover, the principal stages of user application (code) development – editing, compilation, linking, program execution, visual and text result analysis – can be directly activated. The LGO menu effectively supports rapid prototyping, making the application development process faster and easier. 6.

Illustrative numerical results and discussion

To show the viability of an LGO-based solution approach, spherical point arrangement problems for n = 25 will be considered below, for the model versions (2 )–(5). Applying the parameterization (8 ) and fixing three variables by (9), each of these problems has 47 globally optimized decision variables. Although such problem sizes are not very large, in the general “black box” GO context they are considered to be far from trivial. In the input parameter file, the maximal number of global search phase steps – i.e., objective function evaluations – is set to 50000. This setting is admittedly heuristic, but based upon extensive numerical experience a value of Km (in which m is the number of variables, and the integer K satisfies the relations max(m, 100)  K  max(m, 1000)) seems to be sufficient in most practical problems solved by using LGO in the past few years. Note that frequently a much smaller number of steps would be acceptable to generate a globally established estimate of the optimal solution. At the same time, it is obviously possible to “fabricate” GO problems of extreme numerical difficulty, even for small values m. Theoretically established formulae for the number of necessary global search steps, in order to guarantee a given accuracy, lead to exponential functions of the model dimension. (These functions typically contain some unknown constant factors, depending on the function-class considered; consult, for instance the related discussion in Pintér [21].) The global solution estimate can be automatically or interactively refined by a subsequent local search: in the numerical examples below the automatic mode is used (starting with the global adaptive random search mode). Without going into details, we note that the local search stopping rules applied in LGO are based upon usual numerical optimality and convergence criteria. (The actual settings in the input parameter file are in line with the expected accuracy, to enable comparison with known benchmark results.)

224

PINTÉR

As the results presented below will demonstrate, the global search effort leads to reasonable or good quality solutions, within a moderate amount of computational time. We will also illustrate that the test problems discussed are indeed far from trivial, from a general nonlinear optimization viewpoint. All numerical results and corresponding runtimes presented below were obtained in one single illustrative run, for each model version. Since LGO uses also random sample points in the global search phase, a more detailed numerical study could be done, by solving each problem a number of times, and providing the resulting statistics. Here we did not follow this approach for two reasons. First, the results serve merely to illustrate the viability of the GO approach. Second – as a rule – the preset limit of global search steps would be attained in all individual runs, and hence the run-specific total numbers of steps would be fairly close to each other. (The results cited below will illustrate this point.) A Pentium-Pro 200 MHz processor based personal computer, running under the Windows NT Workstation (version 4.0) operating system was used. The numerical values in table 1 are directly taken from the summary output files generated by LGO; most of the details are, however, suppressed for brevity. Additional numerical details – as well as the corresponding LGO executable programs – are available upon request from the author. We emphasize again that the results presented in table 1 are tentative; and that – similarly to other results reported in the literature – their global optimality is not warranted. However, in case of all functions (2 )–(5) benchmark results are known from the literature to which the presented optimum estimates can be compared; this comparison will then indicate the relative quality of the results obtained. All function values cited below are approximate, as being rounded to 10−6 precision. (Note that in some of the works cited below the results are based upon more precise calculations, than the precision shown by table 1; 10−6 was the precision reported by LGO in its standard usage mode, when the illustrative computations reported here were done.) Table 1 Spherical point distribution model versions: summary of illustrative results.

EFP FP PSP TP

FE

EO

ST

57528 57519 57505 63964

35.174804 243.830829 414.623947 0.635677

109 110 115 118

EFP – Elliptic Fekete Problem, the model form (2 ) is used; FP – Fekete Problem, see (3); PSP – Power Sum Problem, see (4), the model parameter value is set to a = 1; TP – Tammes Problem, see (5); FE – Total number of function evaluations; EO – Estimated optimum value; ST – LGO solution time (seconds), including user and system I/O operations.

GLOBALLY OPTIMIZED ARRANGEMENTS

225

For the models (2 )–(4), Rakhmanov, Saff and Zhou [27] present the results of an extensive numerical study. For each model form and corresponding spherical n-point configuration, their results are based upon the best value found using at least 1000 random starting points and several local optimization schemes. Their study also includes a theoretically based ex-post checking and verification of the point configuration structures obtained. Tabulated numerical results are given for n = 2, 3, . . . , 200 points. For the elliptic Fekete problem (2 ) and n = 25, Rakhmanov, Saff and Zhou report the (approximate) optimum value 80.997510 on a natural logarithm scale; this corresponds to the 10-based logarithm value 35.176811. Pintér, Stortelder and de Swart [25] report consistent numerical results for the elliptic Fekete model up to n = 150 points, when applying two entirely different solution approaches. One of these approaches is LGO (in its 1995 implementation), the other one is a differential algebraic equations (DAE) based approach. (For background on the latter method, consult Brenan, Campbell, and Petzold [2], or Hairer and Wanner [8].) In the numerical tests of Pintér, Stortelder and de Swart [25] a high precision SGI Indy workstation with four 194 MHz R1001SC processors was used; the slightly better DAE based result equals 35.176771. Hence, the result cited in table 1 – obtained by LGO in less than 2 minutes on a fairly modest speed PC, and without postulating any prior structure – approximates the best solution value (known to this author) to about 99.9943%. For the Fekete problem (3), Hardin, Sloane and Smith [9] also report putatively optimal arrangements for n = 4, 5, . . . , 132. Similarly to Rakhmanov, Saff and Zhou [27], their investigations make substantial use of a priori postulated symmetries of the point arrangements sought. They report for n = 25 the value 243.812760. This result (as rounded to six decimals) coincides with the value cited by Rakhmanov, Saff and Zhou. Hence, the illustrative LGO result is within 99.9928% of the best known numerical optimum estimate. For the power sum problem (4) using the parameter a = 1, Rakhmanov, Saff and Zhou [27] report the optimum estimate 414.624638. Hence, the value obtained by LGO is within 99.9998% of the best known solution. The fairly small discrepancy between the cited results in the three above models is essentially thought to be a consequence of numerical inaccuracies and stopping parameter settings in the test runs, as well as of using different hardware. Finally, for the – probably by far most difficult – Tammes problem (5), first we cite the classical estimate of Habicht and van der Waerden [7] discussed briefly by Saff and Kuijlaars [30]. According to this result, for the hard-spheres objective function the pair of inequalities (8π )1/2 3−1/4 n−1/2 − Cn−2/3  max min dj k  (8π )1/23−1/4 n−1/2 (C > 0 is some constant)

(10)

is valid. Let us remark that the known theoretical optimum estimates for the models studied are, as a rule, not constructive since – similarly to (10) – they most typically contain some unknown constants. Evidently, the bounds in (10) become more exact, as n increases.

226

PINTÉR

For n = 25, (10) gives the approximate upper bound value 0.761850. Hardin, Sloane and Smith [10] report a 25-point arrangement for which the objective function value of (5) equals 0.710776. Consequently, the value found by LGO – again, based on a fairly limited and entirely “blind” (undirected) global search effort – is within about 89.4342% of the best value known. Especially in the case of this example, it is quite obvious though that good point arrangements need to be symmetric in some way: however, this point is not reflected at all by the formulation considered in this paper (and submitted to LGO). Obviously, a more thorough “black box” search should lead to better results, but – for the purposes of these illustrative calculations – we did not follow this approach. As suggested earlier, probably the PSP is the “easiest” followed by FP, EFP, and TP, in increasing order of complexity. This point is illustrated by figures 1–4. Each picture shows a projection of the corresponding objective function in a subspace that was selected interactively in the corresponding LGO run. Note that – with the exception of function (5) – not all such projections are equally “dramatic”: therefore we selected some of the more interesting pictures generated by LGO. These pictures – reproduced in black and white here – clearly illustrate the non-trivial multiextremality of the various models discussed, and motivate the application of genuine global optimization techniques to such problems. The pictures also show the improving sample point sequence, and its convergence to the optimum estimate generated (the latter is indicated by the larger dark dot on each figure).

Figure 1. A projection of the elliptic Fekete model objective.

GLOBALLY OPTIMIZED ARRANGEMENTS

Figure 2. A projection of the Fekete model objective.

Figure 3. A projection of the power sum model objective.

227

228

PINTÉR

Figure 4. A projection of the Tammes model objective.

7.

Conclusions

“Well-balanced” point distributions on the surface of a sphere are of interest in many areas of scientific modeling. In this work, a general global optimization approach is suggested to analyze such problems, including also their possible far-reaching extensions. The model development and solver system LGO is tentatively applied to four different model versions. Illustrative numerical results and visual representations of the corresponding objective functions are presented. The results show that LGO generates comparable results to those of other – usually far more extensive – numerical studies, in which certain structural symmetries of the configuration sought are postulated a priori. In the context of spherical point arrangements, the review of Saff and Kuijlaars [30] presents a summary of the theoretical background to the model versions (2 )–(5). They postulate that the general pattern for optimal configurations seems to be a regular hexagonal arrangement that is slightly “perturbed” by 12 pentagonal cells, in order to fit on the sphere. This, of course, assumes that the number of points is sufficiently large; furthermore, that it can accommodate such hexagonal-pentagonal formations. Other postulated configurations are also discussed and conjectured by numerous authors. Consult, for instance, Northby [16], Rakhmanov, Saff and Zhou [26,27], Barrón, Gómez and Romero [1], Doye [5], Neumaier [15], Romero, Barrón and Gómez [28], Hartke [11], Rusin [29], and Wales and Scheraga [32], in the context of various point arrangement, crystal and molecular conformation models.

GLOBALLY OPTIMIZED ARRANGEMENTS

229

Evidently, it is both entirely reasonable and possible to draw upon the results of such research, and the underlying scientific expertise. For instance, optimized point configurations can also be produced by starting out directly from – or from an embedding box region around – postulated symmetrical arrangements. Such work can numerically verify the (local) optimality of given configurations, or possibly can lead to improved solutions. In principle, arbitrary conjectured point arrangements – as well as far more general (and quite possibly more difficult) arrangement problems – can also be studied. As mentioned in section 1, model extensions can incorporate arbitrary dimensionality of the embedding space, various feasible sets, weighted point configurations, and other measures of configuration quality. All such model versions, as well as further extensions can be directly considered within the general global optimization framework suggested here. Acknowledgments This work expands upon earlier joint research – related to the elliptic Fekete problem – done with Dr. Walter J.H. Stortelder and Dr. Jacques J.B. de Swart. At the time of that research (1995) both of them were with the National Research Institute for Mathematics and Computer Science (CWI, Amsterdam). The first version of the present paper was written during my visit at Wolfram Research, Inc. (WRI, Champaign, IL). The support received from CWI and WRI during these visits is gratefully acknowledged. Professor E.B. Saff was kind to send several research articles (written by him with co-authors) from which theoretical results, observations and numerical results were cited above. I also wish to thank an anonymous referee for valuable comments that improved the exposition. References [1] C. Barrón, S. Gómez and D. Romero, Lower energy icosahedral atomic clusters with incomplete core, Applied Mathematics Letters 10(5) (1997) 25–28. [2] K.E. Brenan, S.L. Campbell and L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations (North-Holland, Amsterdam, 1989). [3] J.H. Conway and N.J.A. Sloane, Sphere Packings, Lattices and Groups, 2nd ed. (Springer, Berlin, 1988). [4] P.M. Dean, Molecular Similarity in Drug Design (Chapman and Hall, London, 1995). [5] J.P.K. Doye, Lennard–Jones clusters (1999), http://brian.ch.cam.ac.uk/∼jon/structures/LJ.html. [6] L. Greengard, The Rapid Evaluation of Potential Fields in Particle Systems (MIT Press, Cambridge, MA, 1988). [7] W. Habicht and B.L. van der Waerden, Lagerung von Punkten auf der Kugel, Mathematical Annals 123 (1951) 223–234. [8] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed. (Springer, Berlin, 1996). [9] R.H. Hardin, N.J.A. Sloane and W.D. Smith, Minimal energy arrangements of points on a sphere. Arrangements of points on a sphere with minimal 1/r potential (1994), http://www.research.att.com/∼ njas/electrons/index.html. (Last reported web page update: 1997.)

230

PINTÉR

[10] R.H. Hardin, N.J.A. Sloane and W.D. Smith, Minimal-energy clusters of hard spheres (1996), http://www.research.att.com/∼njas/packings/index.html. (Last reported web page update: 1997.) [11] B. Hartke, Global cluster geometry optimization by a phenotype algorithm with niches: Location of elusive minima, and low-order scaling with cluster size, Journal of Computational Chemistry 20(16) (1999) 1752–1759. [12] R. Horst and P.M. Pardalos, eds., Handbook of Global Optimization (Kluwer Academic, Dordrecht, 1995). [13] Lahey Computer Systems, Inc., Fortran 90 User’s Guide (Version 4.5), LCS (Incline Village, NV, 1999). [14] A. Neumaier, Molecular modeling of proteins and mathematical prediction of protein structure, SIAM Review 39 (1997) 407–460. [15] A. Neumaier, Molecular modeling of proteins (1999), http://solon.cma.univie.ac.at/∼neum/protein. html. [16] J.A. Northby, Structure and binding of Lennard–Jones clusters: 13  N  147, Journal of Chemical Physics 87 (1987) 6166–6177. [17] R.H. Pain, ed., Mechanisms of Protein Folding (Oxford University Press, Oxford, 1994). [18] P.M. Pardalos, An open global optimization problem on the unit sphere, Journal of Global Optimization 6 (1995) 213. [19] P.M. Pardalos, D. Shalloway and G. Xue, Optimization methods for computing global minima of nonconvex potential energy functions, Journal of Global Optimization 4 (1994) 117–133. [20] P.M. Pardalos, D. Shalloway and G. Xue, Global Minimization of Nonconvex Energy Functions: Molecular Conformation and Protein Folding, DIMACS Series, Vol. 23 (Amer. Math. Soc., Providence, RI, 1996). [21] J.D. Pintér, Global Optimization in Action (Continuous and Lipschitz Optimization: Algorithms, Implementations and Applications) (Kluwer Academic, Dordrecht, 1996). [22] J.D. Pintér, A model development system for global optimization; in: High Performance Software for Nonlinear Optimizaton: Status and Perspectives, eds. R. De Leone, A. Murli, P.M. Pardalos and G. Toraldo (Kluwer Academic, Dordrecht, 1998) pp. 301–314. [23] J.D. Pintér, LGO – A Model Development System for Continuous Global Optimization. User’s Guide (Pintér Consulting Services, Halifax, NS, 2001). [24] J.D. Pintér, Continuous global optimization: An introduction to models, solution approaches, tests and applications, Interactive Transactions in Operations Research and Management Science 2(2) (1999), http://catt.bus.okstate.edu/itorms/pinter/. [25] J.D. Pintér, W.J.H. Stortelder and J.J.B. de Swart, Computation of elliptic Fekete point sets, Research Report MAS-R9705, National Research Institute for Mathematics and Computer Science (CWI), Amsterdam (1997). (Revised version appeared in: CWI Quarterly 12(1) (1999) 63–76.) [26] E.A. Rakhmanov, E.A. Saff and Y.M. Zhou, Minimal discrete energy on the sphere, Mathematical Research Letters 1 (1994) 647–662. [27] E.A. Rakhmanov, E.A. Saff and Y.M. Zhou, Electrons on the sphere, in: Computational Methods and Function Theory, eds. R.M. Ali, St. Ruscheweyh and E.B. Saff (World Scientific, Singapore, 1995) pp. 111–127. [28] D. Romero, C. Barrón and S. Gómez, The optimal geometry of Lennard–Jones clusters (1999) to appear. [29] D. Rusin, Frequently-asked questions about spheres, http://www.math.niu.edu/∼rusin/known-math/ index/spheres.html (1999). [30] E.B. Saff and A.B.J. Kuijlaars, Distributing many points on a sphere, The Mathematical Intelligencer 19(1) (1997) 5–11. [31] S. Smale, Mathematical problems for the next century, The Mathematical Intelligencer 20(1) (1998) 7–15. [32] D.J. Wales and H.A. Scheraga, Global optimization of clusters, crystals, and biomolecules, Science 285 (1999) 1368–1372.