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 interval-based 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 state-of-the art solvers including mathematical programming methods, population-based 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. State-of-the-art 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 search-space that cannot contain a global minimizer. We introduce interval computations and interval-based 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 floating-point computations. Real numbers that are not representable on a computer are enclosed within intervals with floating-point 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 noteits midpoint. In the following, capital letters represent interval quantities (interval ) and bold letters represent vectors (box , vector ).
Interval arithmetic defines the interval counterparts of real-valued operators () and elementary functions (, , ). For example, and , where (resp. ) denotes downward (resp. upward) rounding.
Let be a real-valued 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 Interval-based 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 search-space 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 floating-point 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 search-space 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 fixed-point 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 interval-based 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:
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 search-space, the bounce-back method [PriceStornLampinen2006] replaces with a component that lies between and the admissible bound:
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
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
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 non-differentiable functions – at can be computed. Charibde handles an interval extension proposed by Kearfott [Kearfott1996], based on the values of its subderivative:
This expression is used to compute an enclosure of the partial derivatives of through automatic differentiation.
4 State-of-the-art Solvers
The certified global minima obtained by Charibde are compared with state-of-the-art optimization softwares available on the NEOS (Network-Enabled Optimization System) server111 http://www.neos-server.org/neos/, a free web service for solving optimization problems in AMPL/GAMS formats. These solvers include local methods, population-based metaheuristics and spatial branch and bound: Ipopt implements a primal-dual interior point algorithm, which uses a filter line search method. LOQO is based on an infeasible primal-dual interior-point method. MINOS employs a reduced-gradient method. PGAPack
is a Parallel Genetic Algorithm library.PSwarm combines pattern search and particle swarm. Couenne uses a reformulation-based branch-and-bound algorithm for a globally optimum solution. It is a complete solver, i.e. it performs an exhaustive exploration of the search-space.
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.
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|
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|
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 population-based 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 search-space. 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|
|(167 eval)||(24 eval)||(15 eval)||(8 eval)||(16 eval)||(8 eval)|
|(88 eval)||(17 eval)||(19 eval)||(5924 eval)||(138 eval)||(50 eval)|
|(3 eval)||(38 eval)||(3 eval)||(3 eval)||(1 eval)||(3 eval)|
|(9582 eval)||(9615 eval)||(9602 eval)||(9626 eval)||(9622 eval)||-|
|(2035 eval)||(2049 eval)||(2008 eval)||(2040 eval)||(2046 eval)||-|
|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 x-axis 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 search-space 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 real-valued 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.
We provided a comparison between Charibde, a cooperative solver that combines an EA and interval-based methods, and state-of-the-art solvers (stemming from mathematical programming and population-based 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.