1 Introduction
Numerical solvers usually embed advanced methods to tackle nonlinear optimization problems. Stochastic methods, in particular evolutionary algorithms, handle large problems and provide sufficiently good solutions with limited computation capacity, but may easily get trapped in local minima. On the other hand, local and global (exhaustive) deterministic methods may guarantee local or global optimality, but are often limited by the size or the nonlinearity of the problems and may suffer from numerical approximations.
In [Vanaret2013], a new reliable hybrid solver named Charibde has been introduced to reconcile stochastic methods and numerical analysis methods. An evolutionary algorithm and an intervalbased algorithm are combined in a cooperative framework: the two methods run in parallel and cooperate by exchanging the best known upper bound of the global minimum and the best current solution. The contribution of this paper is the achievement of new certified optimality results by Charibde for six highly multimodal nonlinear test functions, for which no or few results were available. We also compare Charibde with stateofthe art solvers including mathematical programming methods, populationbased metaheuristics and spatial branch and bound. Charibde proves to be highly competitive with the best solvers on the given test functions, while being fully reliable.
Charibde is presented in Section 2. It is evaluated on a benchmark of difficult test functions given in Section 3. Stateoftheart solvers to which Charibde is compared are described in Section 4. Numerical results, including proofs of optimality, values of global minima and corresponding solutions, are provided and discussed in Section 5.
2 Charibde: a Rigorous Solver
The rigorous nonlinear solver Charibde was introduced by Vanaret et al. [Vanaret2013], building on an original idea by Alliot et al. [Alliot2012]. It combines the efficiency of a Differential Evolution algorithm and the reliability of interval computations to discard more efficiently subspaces of the searchspace that cannot contain a global minimizer. We introduce interval computations and intervalbased methods in Sections 2.1 and 2.2. Details on the implementation of the Differential Evolution algorithm are given in Section 2.3. Finally, the cooperation scheme of Charibde is explained in Section 2.4.
2.1 Interval Analysis
Interval Analysis (IA) is a method of numerical analysis introduced by Moore [Moore1966] to bound rounding errors in floatingpoint computations. Real numbers that are not representable on a computer are enclosed within intervals with floatingpoint bounds. Each numerical computation is safely carried out by using outward rounding.
An interval is the set . We denote by its midpoint. is the set of all intervals. A box
is an interval vector. We note
its midpoint. In the following, capital letters represent interval quantities (interval ) and bold letters represent vectors (box , vector ).Interval arithmetic defines the interval counterparts of realvalued operators () and elementary functions (, , ). For example, and , where (resp. ) denotes downward (resp. upward) rounding.
Let be a realvalued function. is an interval extension of if
The natural interval extension is obtained by replacing elementary operations in with their interval extensions.
Dependency is the main source of overestimation when using interval computations: multiple occurrences of a same variable are considered as different variables. For example, the interval evaluation of over the interval yields , which crudely overestimates the exact range . However, an appropriate rewriting of the syntactic expression of may reduce or overcome dependency: if is continuous inside a box, yields the optimal range when each variable occurs only once in its expression. Completing the square in the expression of provides the optimal syntactic expression . Then .
2.2 Intervalbased Techniques
2.2.1 Interval Branch and Bound Algorithms
Interval Branch and Bound algorithms (IB&B) exploit the conservative properties of interval extensions to rigorously bound global optima of numerical optimization problems [Hansen1992]. The method consists in splitting the initial searchspace into subspaces (branching) on which an interval extension is evaluated (bounding). By keeping track of the best upper bound of the global minimum , boxes that certainly do not contain a global minimizer are discarded (example 2.2.1). Remaining boxes are stored to be processed at a later stage until the desired precision is reached. The process is repeated until all boxes have been processed. Convergence certifies that , even in the presence of rounding errors. However, the exponential complexity of IB&B hinders the speed of convergence on large problems.
Consider the problem over . Then . The floatingpoint evaluation provides an upper bound of . Evaluating on the subinterval reduces the overestimation induced by dependency: . Because , the interval cannot contain a global minimizer and can be safely discarded.
2.2.2 Interval Contraction
Propagating the (in)equality constraints of the problem, as well as the constraints and , may narrow the domains of the variables or prove that a subdomain of the searchspace cannot contain a global minimizer.
Stemming from the IA and Interval Constraint Programming communities, filtering/contraction algorithms [ChabertJaulin2009] narrow the bounds of the variables without loss of solutions. Standard contraction algorithms generally integrate a filtering procedure into a fixedpoint algorithm. HC4 [Benhamou1999] handles one constraint after the other and performs the optimal contraction w.r.t. to a constraint if variables occur only once in its expression. Box [VanHentenryck1997] narrows one variable after the other w.r.t. all constraints, using an interval version of Newton’s method. Mohc [Araya2010] exploits the monotonicity of the constraints to enhance contraction of HC4 and interval Newton.
The intervalbased algorithm embedded in Charibde follows an Interval Branch and Contract (IB&C) scheme (algorithm 1) that interleaves steps of bisection and filtering. We note the priority queue in which the remaining boxes are stored, the desired precision and the best known solution, such that .
2.3 Differential Evolution
Differential Evolution (DE) is an Evolutionary Algorithm that combines the coordinates of existing individuals with a particular probability to generate new potential solutions
[StornPrice1997]. It was embedded within Charibde for its ability to solve extremely difficult optimization problems, while having few control parameters.We denote by the population size, the weighting factor and the crossover rate. For each individual of the population, three other individuals , and , all different and different from , are randomly picked in the population. The newly generated individual is computed as follows:
(1) 
where is a random index in ensuring that at least one component of differs from that of , and
is a random number uniformly distributed in
, picked for each . replaces in the population if .The following advanced rules have been implemented in Charibde:
 Boundary constraints:

When a component lies outside the bounds of the searchspace, the bounceback method [PriceStornLampinen2006] replaces with a component that lies between and the admissible bound:
(2)  Evaluation:

Given inequality constraints , the evaluation of an individual is computed as a triplet (, , ), where is the objective value of , the number of violated constraints and . If at least one of the constraints is violated, the objective value is not computed
 Selection:

Given the evaluation triplets () and () of two candidate solutions and , the best individual to be kept for the next generation is computed as follows:

if or ( and ) or ( and ) then is kept

otherwise, replaces

2.4 Charibde: a Cooperative Algorithm
Charibde combines an Interval Branch and Contract algorithm and a Differential Evolution algorithm in a cooperative way: neither of the algorithms is embedded within the other, but they run in parallel and exchange bounds and solutions using an MPI implementation (Figure 1).
The cooperation scheme boils down to 3 main steps:

Whenever the best known DE evaluation is improved, the best individual is evaluated using IA. The upper bound of the image – an upper bound of the global minimum – is sent to the IB&C thread

In the IB&C algorithm, is compared to the current best upper bound . An improvement of the latter leads to a more efficiently pruning of the subspaces that cannot contain a (feasible) global minimizer

Whenever the evaluation of the center of a box improves , the individual replaces the worst individual of DE, thus preventing premature convergence
3 Benchmark of Test Functions
The highly multimodal nonlinear test functions considered in this study can be found in Table 1. Contrary to standard test functions that have a global minimum 0 at (Griewank function) or have a global minimizer with identical components (Schwefel function), we have selected six functions with nontrivial global minima: Michalewicz function, Sine Envelope Sine Wave function (shortened to Sine Envelope), Shekel’s Foxholes function [PriceStornLampinen2006], Egg Holder function [Whitley1996], Rana’s function [Whitley1996] and Keane’s function [Keane1994]. Except for the Michalewicz function, all are nonseparable. Their surfaces and contour lines can be observed on Figure 2 for . It illustrates the numerous local minima and the ruggedness of the functions.
The first inequality constraint of Keane’s function describes a hyperbola in two dimensions and is active at the global minimizer, which hinders the efficiency of solvers. The second inequality constraint is linear and is not active at the global minimizer. The Egg Holder (resp. Rana) function is strongly subject to dependency: and occur three (resp. five) times in its expression, and occur six (resp. ten) times. Its natural interval extension therefore produces a large overestimation of the actual range.
The last three functions (Egg Holder, Rana and Keane) contain absolute values. is differentiable everywhere except for , however its subderivative – generalizing the derivative to nondifferentiable functions – at can be computed. Charibde handles an interval extension proposed by Kearfott [Kearfott1996], based on the values of its subderivative:
(3) 
This expression is used to compute an enclosure of the partial derivatives of through automatic differentiation.
Function  Expression  Domain 

Michalewicz  
Sine Envelope  
Shekel  
Egg Holder  
Rana  
Keane  s.t. and 
4 Stateoftheart Solvers
The certified global minima obtained by Charibde are compared with stateoftheart optimization softwares available on the NEOS (NetworkEnabled Optimization System) server^{1}^{1}1 http://www.neosserver.org/neos/, a free web service for solving optimization problems in AMPL/GAMS formats. These solvers include local methods, populationbased metaheuristics and spatial branch and bound: Ipopt implements a primaldual interior point algorithm, which uses a filter line search method. LOQO is based on an infeasible primaldual interiorpoint method. MINOS employs a reducedgradient method. PGAPack
is a Parallel Genetic Algorithm library.
PSwarm combines pattern search and particle swarm. Couenne uses a reformulationbased branchandbound algorithm for a globally optimum solution. It is a complete solver, i.e. it performs an exhaustive exploration of the searchspace.Note: the solver BARON (Branch and Reduce Optimization Navigator), combining constraint propagation, interval analysis and duality, is known as one of the most efficient, although unreliable (see Section 5.3), solvers for solving nonconvex optimization problems to global optimality. However, it is not designed to handle trigonometric functions, and therefore will not be considered in this paper.
5 Numerical Experiments
5.1 Global Minima Certified by Charibde
The global minima to the six test functions are given in Table 6. All solutions are given with a precision . Should the functions have several global minima, only one is provided. A reference given next to a global minimum indicates that a result is available in the literature.
Analyzing the results has brought out identities for finding putative minima for four of the six test functions. They can be found in Table 2.
Function  Putative minimum  

Michalewicz  0.9999981  
Sine Envelope  1  
Egg Holder  0.9999950  
Rana’s  1 
The global minima of the Michalewicz function for up to 75 variables can be found in Table 3. For the sake of conciseness, the corresponding solutions for up to only 10 variables are given in Table 6.
Global minimum  Global minimum  

10  9.6601517 [Mishra2006]  45  44.6256251 
15  14.6464002  50  49.6248323 [Mishra2006] 
20  19.6370136 [Mishra2006]  55  54.6240533 
25  24.6331947  60  59.6231462 
30  29.6308839 [Mishra2006]  65  64.6226167 
35  34.6288550  70  69.6222202 
40  39.6267489  75  74.6218112 
To illustrate the key role played by the syntactic expression of the function when computing with intervals, we have tested two different – but equivalent – syntaxes of Rana’s function in Charibde. The first syntax is given in Table 1. The second syntax is obtained using the trigonometric identity: . Their impact on the sharpness of the inclusions, therefore on the convergence of Charibde, can be observed in Table 4. The hyphen indicates a computing time greater than one hour.
CPU time (s)  
First syntax  Second syntax  
2  0.25  0.009 
3  6.5  0.12 
4  254  1.45 
5    18.5 
6    244 
7    3300 
5.2 Comparison of Solvers
The comparison of the seven solvers on a particular instance of each test function can be found in Table 5. When available, the number of evaluations of the objective function is given under the found minimum, otherwise we have mentioned the computing time. The hyphens in the last column indicate that PGAPack and PSwarm cannot handle inequality constraints.
Local methods based on mathematical programming (Ipopt, LOQO, MINOS) usually require few iterations to reach a local minimum from a starting point. The quality of the minimum generally depends on both the starting point and the size of the basins of attraction of the function. These solvers turn out to perform poorly on the considered multimodal problems.
Among the populationbased metaheuristics, PGAPack performs consistently better than PSwarm, yet at a higher cost. We have kept the best result ouf of five runs, since it is tedious to run an online solver several times. The default NEOS control parameters have been used; both algorithms would certainly perform better with suitable control parameters. The results of Couenne are discussed in Section 5.3.
Charibde achieves convergence in finite time on the six problems, with a numerical certification of optimality. It benefits from the start of convergence of the DE algorithm that computes a good initial value for . This allows the IB&C algorithm to prune more efficiently subspaces of the searchspace. The number of evaluations of the natural interval extension has the form , where is the number of evaluations in the DE algorithm (whenever the best known solution is improved) and is the number of evaluations in the IB&C algorithm. Note that after converging toward the global minimizer, the DE thread keeps running as long as the certification of optimality has not been obtained.
Michalewicz  Sine Envelope  Shekel  Egg Holder  Rana  Keane  
()  ()  ()  ()  ()  ()  
Ipopt  19.773742  5.8351843  1.8296708  3586.3131827  75.512076  0.2010427 
(167 eval)  (24 eval)  (15 eval)  (8 eval)  (16 eval)  (8 eval)  
LOQO  0.0048572  5.8351843  1.9421143  25.5609238  69.5206  0.0983083 
(88 eval)  (17 eval)  (19 eval)  (5924 eval)  (138 eval)  (50 eval)  
MINOS  0  5.87878  2.5589419  3586.313183  233.592  0.2347459 
(3 eval)  (38 eval)  (3 eval)  (3 eval)  (1 eval)  (3 eval)  
PGAPack  37.60465  5.569544  1.829452  3010.073  2091.068   
(9582 eval)  (9615 eval)  (9602 eval)  (9626 eval)  (9622 eval)    
PSwarm  24.38158  5.835182  1.610072  2840.799  1595.056   
(2035 eval)  (2049 eval)  (2008 eval)  (2040 eval)  (2046 eval)    
Couenne  49.619042  5.9660007  10.4039521  3719.7287498  2046.8320657  0.6222999 
(265s)  (0.4s)  (8s)  (20.6s)  (20.3s)  (2s)  
Charibde  49.6248323  5.9659811  10.4039521  3719.7248363  2046.8320657  0.6222810 
CPU time  8.4s  219s  0.04s  0.8s  17.8s  14860s 
evaluations  717100  227911400  8150  742600  19186150  6351693297 
evaluations  763 + 409769  93 + 21744667  28 + 561  106 + 82751  53 + 1383960  74 + 2036988566 
()  (50, 0.7, 0)  (50, 0.7, 0.9)  (50, 0.7, 0.9)  (50, 0.7, 0.4)  (50, 0.7, 0.5)  (70, 0.7, 0.9) 
5.3 Reliability vs Efficiency
On Figure 3, the evolution of the best known solution in Charibde is compared with that of Couenne for a particular instance of each test function. Intermediate times for other solvers were not available. Note that the scale on the xaxis is logarithmic. These diagrams show that Charibde is highly competitive against Couenne: Charibde achieves convergence faster than Couenne on Michalewicz (ratio 31), Shekel (ratio 200), Egg Holder (ratio 25) and Rana’s function (ratio 1.1). Couenne is faster than Charibde on Sine Envelope (ratio 547) and Keane’s function (ratio 7430).
Is is however crucial to note that Couenne, while being a complete solver (the whole searchspace is exhaustively processed), is not reliable. This stems from the fact that the under and overestimators obtained by linearizing the function may suffer from numerical approximations. Contrary to interval arithmetic that bounds rounding errors, mere realvalued linearizations are not conservative and cannot guarantee the correctness of the result. This problem can easily be observed in Table 5: the global minima obtained by Couenne on Michalewicz, Sine Envelope, Egg Holder and Keane’s function are not correct compared to the certified minima provided by Charibde. The wrong decimal places are underlined.
6 Conclusion
We provided a comparison between Charibde, a cooperative solver that combines an EA and intervalbased methods, and stateoftheart solvers (stemming from mathematical programming and populationbased metaheuristics). They were evaluated on a benchmark of nonlinear multimodal optimization problems among the most challenging: Michalewicz, Sine Envelope, Shekel’s Foxholes, Egg Holder, Rana and Keane. Charibde proved to be highly competitive with the best solvers, including Couenne, a complete but unreliable solver based on spatial branch and bound and linearizations. We provided new certified global minima for the considered test functions, as well as the corresponding solutions. They may be used from now on as references to test stochastic or deterministic optimization methods.
Comments
There are no comments yet.