1 Motivation
We consider dimensional continuous constrained optimization problems over a hyperrectangular domain :
(1)  
When , and are nonconvex, the problem may have multiple local minima. Such difficult problems are generally solved using generic exhaustive branch and bound (BB) methods. The objective function and constraints are bounded on disjoint subspaces using enclosure methods. By keeping track of the best known upper bound of the global minimum , subspaces that cannot contain a global minimizer are discarded (pruned).
Several authors proposed hybrid approaches in which a BB algorithm cooperates with another technique to enhance the pruning of the searchspace. Hybrid algorithms may be classified into two categories
[24]: integrative approaches, in which one of the two methods replaces a particular operator of the other method, and cooperative methods, in which the methods are independent and are run sequentially or in parallel. Previous works include
integrative approaches: [34]
integrates a stochastic genetic algorithm (GA) within an interval BB. The GA provides the direction along which a box is partitioned, and an individual is generated within each subbox. At each generation, the best evaluation updates the best known upper bound of the global minimum. In
[9], the crossover operator is replaced by a BB that determines the best individual among the offspring. 
cooperative approaches: [28] sequentially combines an interval BB and a GA. The interval BB generates a list of remaining small boxes. The GA’s population is initialized by generating a single individual within each box of . [12] (BB and memetic algorithm) and [7] (beam search and memetic algorithm) describe similar parallel strategies: the BB identifies promising regions that are then explored by the metaheuristic. [1] hybridizes a GA and an interval BB. The two independent algorithms exchange upper bounds and solutions through shared memory. New optimal results are presented for the rotated Michalewicz () and Griewank functions ().
In this communication, we build upon the cooperative scheme of [1]. The efficiency and reliability of their solver remain very limited; it is not competitive against stateoftheart solvers. Their interval techniques are naive and may lose solutions, while the GA may send evaluations subject to roundoff errors. We propose to hybridize a stochastic differential evolution algorithm (close to a GA), described in Section 2, and a deterministic interval branch and contract algorithm, described in Section 3. Our hybrid solver Charibde is presented in Section 4. Experimental results (Section 5) show that Charibde is highly competitive against stateoftheart solvers.
2 Differential Evolution
Differential evolution (DE) [29]
is among the simplest and most efficient metaheuristics for continuous problems. It combines the coordinates of existing individuals (candidate solutions) with a given probability to generate new individuals. Initially devised for continuous unconstrained problems, DE was extended to mixed problems and constrained problems
[23].Let denote the size of the population, the amplitude factor and the crossover rate. At each generation (iteration), new individuals are generated: for each individual , three other individuals (called base individual), and , all different and different from , are randomly picked in the population. The coordinates of the new individual are computed according to
(2) 
where is picked in with uniform probability. The index , picked in with uniform probability for each , ensures that at least a coordinate of differs from that of . replaces in the population if it is “better” than (e.g. in unconstrained optimization, is better than if it improves the objective function).
Figure 1 depicts a twodimensional crossover between individuals , (base individual), and . The contour lines of the objective function are shown in grey. The difference , scaled by , yields the direction (an approximation of the direction opposite the gradient) along which is translated to yield .
3 Reliable Computations
Reliable (or rigorous) methods provide bounds on the global minimum, even in the presence of roundoff errors. The only reliable approaches for achieving a numerical proof of optimality in global optimization are intervalbased methods that interleave branching of the searchspace and pruning of the subdomains that cannot contain an optimal solution.
Section 3.1 introduces interval arithmetic, an extension of real arithmetic. Reliable global optimization is detailed in Section 3.2, and interval contractors are mentioned in Section 3.3.
3.1 Interval Arithmetic
An interval with floatingpoint bounds defines the set . denotes the set of all intervals. The width of is . is the midpoint of . A box is a Cartesian product of intervals. The width of a box is the maximum width of its components. The convex hull of and is the smallest interval enclosing and .
Interval arithmetic [19] extends real arithmetic to intervals. Interval arithmetic implemented on a machine must be rounded outward (the left bound is rounded toward , the right bound toward ) to guarantee conservative properties. The interval counterparts of binary operations and elementary functions produce the smallest interval containing the image. Thanks to the conservative properties of interval arithmetic, we define interval extensions (Definition 1) of functions that may be expressed as a finite composition of elementary functions.
Definition 1 (Interval extension)
Let . is an interval extension (or inclusion function) of iff
(3)  
Interval extensions with various sharpnesses may be defined (Example 1). The natural interval extension replaces the variables with their domains and the elementary functions with their interval counterparts. The Taylor interval extension is based on the Taylor expansion at point .
Example 1 (Interval extensions)
Let , and . The exact range is . Then

;

.
Example 1 shows that interval arithmetic often overestimates the range of a realvalued function. This is due to the dependency problem, an inherent behavior of interval arithmetic. Dependency decorrelates multiple occurrences of the same variable in an analytical expression (Example 2).
Example 2 (Dependency)
Let . Then
(4)  
Interval extensions have different convergence orders, that is the overestimation decreases at different speeds with the width of the interval.
3.2 Global Optimization
Interval arithmetic computes a rigorous enclosure of the range of a function over a box. The first branch and bound algorithms for continuous optimization based on interval arithmetic were devised in the 1970s [20][27], then refined during the following years [14]: the searchspace is partitioned into subboxes. The objective function and constraints are evaluated on each subbox using interval arithmetic. The subspaces that cannot contain a global minimizer are discarded and are not further explored. The algorithm terminates when .
To overcome the pessimistic enclosures of interval arithmetic, interval branch and bound algorithms have recently been endowed with filtering algorithms (Section 3.3) that narrow the bounds of the boxes without loss of solutions. Stemming from the Interval Analysis and Interval Constraint Programming communities, filtering (or contraction) algorithms discard values from the domains by enforcing local (each constraint individually) or global (all constraints simultaneously) consistencies. The resulting methods, called interval branch and contract (IBC) algorithms, interleave steps of contraction and steps of bisection.
3.3 Interval Contractors
Stateoftheart contractors (contraction algorithms) include HC4 [6], Box [32], Mohc [2], 3B [17], CID [31] and XNewton [3]. Only HC4 and XNewton are used in this communication.
HC4Revise is a twophase algorithm that exploits the syntax tree of a constraint to contract each occurrence of the variables. The first phase (evaluation) evaluates each node (elementary function) using interval arithmetic. The second phase (propagation) uses projection functions to inverse each elementary function. HC4Revise is generally invoked as the revise procedure (subcontractor) of HC4, an AC3like propagation loop.
XNewton computes an outer linear relaxation of the objective function and the constraints, then computes a lower bound of the initial problem using LP techniques (e.g. the simplex algorithm). additional calls may contract the domains of the variables.
4 Charibde: a Cooperative Approach
4.1 Hybridization of Stochastic and Deterministic Techniques
Our hybrid algorithm Charibde (Cooperative Hybrid Algorithm using Reliable IntervalBased methods and Differential Evolution), written in OCaml [16], combines a stochastic DE and a deterministic IBC for nonconvex constrained optimization. Although it embeds a stochastic component, Charibde is a fully rigorous solver.
4.1.1 Previous Work
Preliminary results of a basic version of Charibde were published in 2013 [33]
on classical multimodal problems (7 boundconstrained and 4 inequalityconstrained problems) widely used in the Evolutionary Computation community. We showed that Charibde benefited from the start of convergence of the DE algorithm, and completed the proof of convergence faster than a standard IBC algorithm. We provided new optimal results for 3 problems (Rana, Eggholder and Michalewicz).
4.1.2 Contributions
In this communication, we present two contributions:

we devised a new cooperative exploration strategy MaxDist that

selects boxes to be explored in a novel manner;

periodically reduces DE’s domain;

restarts the population within the new (smaller) domain.
An example illustrates the evolution of the domain without loss of solutions;


we assess the performance of Charibde against stateoftheart rigorous (GlobSol, IBBA, Ibex) and nonrigorous (Couenne, BARON) solvers on a benchmark of difficult problems.
4.1.3 Cooperative Scheme
Two independent parallel processes exchange bounds, solutions and searchspace via MPI message passing (Figure 2).
The cooperation boils down to three main steps:

whenever the DE improves its best evaluation, the best individual and its evaluation are sent to the IBC to update the best known upper bound ;

whenever the IBC finds a better punctual solution (e.g. the center of a box), it is injected into DE’s population;

the exploration strategy MaxDist periodically reduces the searchspace of DE, then regenerates the population in the new searchspace.
4.2 Differential Evolution
Populationbased metaheuristics, in particular DE, are endowed with mechanisms that help escape local minima. They are quite naturally recommended to solve difficult multimodal problems for which traditional methods struggle to converge. They are also capable of generating feasible solutions without any a priori knowledge of the topology. DE has proven greatly beneficial for improving the best known upper bound , a task for which standard branch and bound algorithms are not intrinsically intended.
4.2.1 Base Individual
In the standard DE strategy, all the current individuals have the same probability to be selected as the base individual . We opted for an alternative strategy [23] that guarantees that all individuals of the population play this role once and only once at each generation: the index of the base individual is obtained by summing the index of the individual and an offset in , drawn with uniform probability.
4.2.2 Bound Constraints
When a coordinate of (computed during the crossover) exceeds the bounds of the component of the domain , the bounceback method [23] replaces with a valid coordinate that lies between the base coordinate and the violated bound:
(5) 
where is drawn in with uniform probability.
4.2.3 Constraint Handling
The extension of evolutionary algorithms to constrained optimization has been addressed by numerous authors. We implemented the direct constraint handling [23]
that assigns to each individual a vector of evaluations (objective function and constraints), and selects the new individual
(see Section 2) based upon simple rules:
and are feasible and has a lower or equal objective value than ;

is feasible and is not;

and are infeasible, and does not violate any constraint more than .
4.2.4 Rigorous Feasibility
Numerous NLP solvers tolerate a slight violation (relaxation) of the inequality constraints (e.g. instead of ). The evaluation of a “pseudofeasible” solution (that satisfies such relaxed constraints) is not a rigorous upper bound of the global minimum; roundoff errors may even lead to absurd conclusions: may be lower than the global minimum, and (or) may be very distant from actual feasible solutions in the searchspace.
To ensure that an individual is numerically feasible (i.e. that the evaluations of the constraints are nonpositive), we evaluate the constraints using interval arithmetic. is considered as feasible when the interval evaluations are nonpositive, that is .
4.2.5 Rigorous Objective Function
When is a feasible point, the evaluation may be subject to roundoff errors; the only reliable upper bound of the global minimum available is
(the right bound of the interval evaluation). However, evaluating the individuals using only interval arithmetic is much costlier than cheap floatingpoint arithmetic.
An efficient inbetween solution consists in systematically computing the floatingpoint evaluations , and computing the interval evaluation only when the best known approximate evaluation is improved. is then compared to the best known reliable upper bound : if is improved, is sent to the IBC. This implementation greatly reduces the cost of evaluations, while ensuring that all the values sent to the IBC are rigorous.
4.3 Interval Branch and Contract
Branching aims at refining the computation of lower bounds of the functions using interval arithmetic. Two strategies may be found in the early literature:

the variable with the largest domain is bisected;

the variables are bisected one after the other in a roundrobin scheme.
More recently, the Smear heuristic
[10] has emerged as a competitive alternative to the two standard strategies. The variable for which the interval quantity is the largest is bisected.Charibde’s main contractor is detailed in Algorithm 3. We exploit the contracted nodes of HC4Revise to compute partial derivatives via automatic differentiation [26]. HC4Revise is a revise procedure within a quasifixed point algorithm with tolerance : the propagation loop stops when the box is not sufficiently contracted, i.e. when the size of becomes larger than a fraction of the initial size . Most contractors include an evaluation phase that yields a lower bound of the problem on the current box. Charibde thus computes several lower bounds (natural, Taylor, LP) as long as the box is not discarded. Charibde calls ocamlglpk [18]
, an OCaml binding for GLPK (GNU Linear Programming Kit). Since the solution of the linear program is computed using floatingpoint arithmetic, it may be subject to roundoff errors. A cheap postprocessing step
[21] computes a rigorous bound on the optimal solution of the linear program, thus providing a rigorous lower bound of the initial problem.4.4 MaxDist: a New Exploration Strategy
The boxes that cannot be discarded are stored in a priority queue to be processed at a later stage. The order in which the boxes are extracted determines the exploration strategy of the searchspace (“bestfirst”, “largest first”, “depthfirst”). Numerical tests suggest that

the “bestfirst” strategy is rarely relevant because of the overestimated range (due to the dependency problem);

the “largest first” strategy does not give advantage to promising regions;

the “depthfirst” strategy tends to quickly explore the neighborhood of local minima, but struggles to escape from them.
We propose a new exploration strategy called MaxDist. It consists in extracting from the box that is the farthest from the current solution . The underlying ideas are to

explore the neighborhood of the global minimizer (a tedious task when the objective function is flat in this neighborhood) only when the best possible upper bound is available;

explore regions of the searchspace that are hardly accessible by the DE algorithm.
The distance between a point and a box is the sum of the distances between each coordinate and the closest bound of . Note that MaxDist is an adaptive heuristic: whenever the best known solution is updated, is reordered according to the new priorities of the boxes.
Preliminary results (not presented here) suggest that MaxDist is competitive with standard strategies. However, the most interesting observation lies in the behavior of : when using MaxDist, the maximum size of (the maximum number of boxes simultaneously stored in ) remains remarkably low (a few dozens compared to several thousands for standard strategies). This offers promising perspectives for the cooperation between DE and IBC: the remaining boxes of the IBC may be exploited in the DE to avoid exploring regions of the searchspace that have already been proven infeasible or suboptimal.
The following numerical example illustrates how the remaining boxes are exploited to reduce DE’s domain through the generations. Let
(7)  
s.t.  
be a constrained optimization problem defined on the box (Figure (a)a). The dotted curves represent the frontiers of the two inequality constraints, and the contour lines of the objective function are shown in solid dark. The feasible region is the bananashaped set, and the global minimizer is located in its lower right corner.
The initial domain of DE (which corresponds to the initial box in the IBC) is first contracted with respect to the constraints of the problem. The initial population of DE is then generated within this contracted domain, thus avoiding obvious infeasible regions of the searchspace. This approach is similar to that of [11]. Figure (b)b depicts the contraction (the black rectangle) of the initial domain with respect to the constraints (sequentially handled by HC4Revise): .
Periodically, we compute the convex hull of the remaining boxes of and replace DE’s domain with . Note that

the convex hull (linear complexity) may be computed at low cost, because the size of remains small when using MaxDist;

by construction, MaxDist handles boxes on the rim of the remaining domain (the boxes of ), which boosts the reduction of the convex hull.
Figures (c)c and (d)d represent the convex hull of the remaining subboxes in the IBC, respectively after 10 and 20 DE generations. The population is then randomly regenerated within the new contracted domain . The convex hull operation progressively eliminates local minima and infeasible regions. The global minimum eventually found by Charibde with precision is ; both constraints are active.
5 Experimental Results
Currently, GlobSol [15], IBBA [22] and Ibex [8] are among the most efficient solvers in rigorous constrained optimization. They share a common skeleton of interval branch and bound algorithm, but differ in the acceleration techniques. GlobSol uses the reformulationlinearization technique (RLT), that introduces new auxiliary variables for each intermediary operation. IBBA calls a contractor similar to HC4Revise, and computes a relaxation of the system of constraints using affine arithmetic. Ibex is dedicated to both numerical CSPs and constrained optimization; it embeds most of the aforementioned contractors (HC4, 3B, Mohc, CID, XNewton). Couenne [5] and BARON [25] are stateoftheart NLP solvers. They are based on a nonrigorous spatial branch and bound algorithm, in which the objective function and the constraints are over and underestimated by convex relaxations. Although they perform an exhaustive exploration of the searchspace, they cannot guarantee a given precision on the value of the optimum.
All five solvers and Charibde are compared on a subset of 11 COCONUT constrained problems (Table 1), extracted by Araya [3] for their difficulty: ex2_1_7, ex2_1_9, ex6_2_6, ex6_2_8, ex6_2_9, ex6_2_11, ex6_2_12, ex7_2_3, ex7_3_5, ex14_1_7 and ex14_2_7. Because of numerical instabilities of the ocamlglpk LP library (“assert failure”), the results of the problems ex6_1_1, ex6_1_3 and ex_6_2_10 are not presented. The second and third columns give respectively the number of variables and the number of constraints . The fourth (resp. fifth) column specifies the type of the objective function (resp. the constraints): L is linear, Q is quadratic and NL is nonlinear. The logsize of the domain (sixth column) is .
The comparison of CPU times (in seconds) for solvers GlobSol, IBBA, Ibex, Couenne, BARON and Charibde on the benchmark of 11 problems is detailed in Table 2
. Mean times and standard deviations (in brackets) are given for Charibde over 100 runs. The numerical precision on the objective function
and the tolerance for equality constaints were identical for all solvers. TO (timeout) indicates that a solver could not solve a problem within one hour. The results of GlobSol (proprietary piece of software) were not available for all problems; only those mentioned in [22] are presented. The results of IBBA were also taken from [22]. The results of Ibex were taken from [3]: only the best strategy (simplex, XNewIter or XNewton) for each benchmark problem is presented. Couenne and BARON (only the commercial version of the code is available) were run on the NEOS server [13].Type  

Problem  Domain logsize  
ex2_1_7  20  10  Q  L  
ex2_1_9  10  1  Q  L  
ex6_2_6  3  1  NL  L  
ex6_2_8  3  1  NL  L  
ex6_2_9  4  2  NL  L  
ex6_2_11  3  1  NL  L  
ex6_2_12  4  2  NL  L  
ex7_2_3  8  6  L  NL  
ex7_3_5  13  15  L  NL  
ex14_1_7  10  17  L  NL  
ex14_2_7  6  9  L  NL 
Rigorous  Non rigorous  

Problem  GlobSol  IBBA  Ibex  Charibde  Couenne  BARON 
ex2_1_7  16.7  7.74  34.9 (13.3)  476  16.23  
ex2_1_9  154  9.07  35.9 (0.29)  3.01  3.58  
ex6_2_6  306  1575  136  3.3 (0.41)  TO  5.7 
ex6_2_8  204  458  59.3  2.9 (0.37)  TO  TO 
ex6_2_9  463  523  25.2  2.7 (0.03)  TO  TO 
ex6_2_11  273  140  7.51  1.96 (0.06)  TO  TO 
ex6_2_12  196  112  22.2  8.8 (0.17)  TO  TO 
ex7_2_3  TO  544  1.9 (0.30)  TO  TO  
ex7_3_5  TO  28.91  4.5 (0.09)  TO  4.95  
ex14_1_7  TO  406  4.2 (0.13)  13.86  0.56  
ex14_2_7  TO  66.39  0.2 (0.04)  0.01  0.02  
Sum  1442  TO  1312.32  101.26  TO  TO 
Charibde was run on an Intel Xeon E31270 @ 3.40GHz x 8 with 7.8 GB of RAM. BARON and Couenne were run on 2 Intel Xeon X5660 @ 2.8GHz x 12 with 64 GB of RAM. IBBA and Ibex were run on similar processors (Intel x86, 3GHz). The difference in CPU time between computers is about 10% [4], which makes the comparison quite fair.
The hyperparameters of Charibde for the benchmark problems are given in Table
3; is the population size, and is the quasifixed point ratio. The amplitude , the crossover rate and the MaxDist strategy are common to all problems. Tuning the hyperparameters is generally problemdependent, and requires structural knowledge about the problem: the population size may be set according to the dimension and the number of local minima, the crossover rate is related to the separability of the problem, and the techniques based on linear relaxation have little influence for problems with few constraints, but are cheap when the constraints are linear.Problem  Bisections  Fixedpoint ratio ()  LP  XNewton  

ex2_1_7  20  RR  0.9  ✓  ✓ 
ex2_1_9  100  RR  0.8  ✓  
ex6_2_6  30  Smear  0  ✓  
ex6_2_8  30  Smear  0  ✓  
ex6_2_9  70  Smear  0  
ex6_2_11  35  Smear  0  
ex6_2_12  35  RR  0  ✓  
ex7_2_3  40  Largest  0  ✓  ✓ 
ex7_3_5  30  RR  0  ✓  
ex14_1_7  40  RR  0  ✓  
ex14_2_7  40  RR  0  ✓ 
Charibde outperforms Ibex on 9 out of 11 problems, IBBA on 10 out of 11 problems and GlobSol on all the available problems. The cumulated CPU time shows that Charibde (101.26s) improves the performances of Ibex (1312.32s) by an order of magnitude (ratio: 13) on this benchmark of 11 difficult problems. Charibde also proves highly competitive against nonrigorous solvers Couenne and BARON. The latter are faster or have similar CPU times on some of the 11 problems, however they both timeout on at least five problems (seven for Couenne, five for BARON). Overall, Charibde seems more robust and solves all the problems of the benchmark, while providing a numerical proof of optimality. Surprisingly, the convergence times do not seem directly related to the dimensions of the instances. They may be explained by the nature of the objective function and constraints (in particular, Charibde seems to struggle when the objective function is quadratic) and the dependency induced by the multiple occurrences of the variables.
Table 4 presents the best upper bounds obtained by Charibde, Couenne and BARON at the end of convergence (precision reached or timeout). Truncated digits on the upper bounds are bounded (e.g. denotes and denotes ). The incorrect digits of the global minima obtained by Couenne and BARON are underlined. This demonstrates that nonrigorous solvers may be affected by roundoff errors, and may provide solutions that are infeasible or have an objective value lower than the global minimum (Couenne on ex2_1_9, BARON on ex2_1_7, ex2_1_9, ex6_2_8, ex6_2_12, ex7_2_3 and ex7_3_5). For the most difficult instance ex7_2_3, Couenne is not capable of finding a feasible solution with a satisfactory evaluation within one hour. It would be very informative to compute the ratio between the size of the feasible domain (the set of all feasible points) and the size of the entire domain. On the other hand, the strategy MaxDist within Charibde greatly contributes to finding an excellent upper bound of the global minimum, which significantly accelerates the interval pruning phase.
Problem  Charibde  Couenne  BARON 

ex2_1_7  
ex2_1_9  
ex6_2_6  
ex6_2_8  
ex6_2_9  
ex6_2_11  
ex6_2_12  
ex7_2_3  
ex7_3_5  
ex14_1_7  
ex14_2_7 
6 Conclusion
We proposed a cooperative hybrid solver Charibde, in which a deterministic interval branch and contract cooperates with a stochastic differential evolution algorithm. The two independent algorithms run in parallel and exchange bounds, solutions and searchspace in an advanced manner via message passing. The domain of the populationbased metaheuristic is periodically reduced by removing local minima and infeasible regions detected by the branch and bound.
A comparison of Charibde with stateoftheart intervalbased solvers (GlobSol, IBBA, Ibex) and NLP solvers (Couenne, BARON) on a benchmark of difficult COCONUT problems shows that Charibde is highly competitive against nonrigorous solvers (while bounding the global minimum) and converges faster than rigorous solvers by an order of magnitude.
References

[1]
Alliot, J.M., Durand, N., Gianazza, D., Gotteland, J.B.: Finding and proving the optimum: Cooperative stochastic and deterministic search. In: 20th European Conference on Artificial Intelligence (ECAI 2012), August 2731, 2012, Montpellier, France (2012)
 [2] Araya, I., Trombettoni, G., Neveu, B.: Exploiting monotonicity in interval constraint propagation. In: Proc. AAAI. pp. 9–14 (2010)
 [3] Araya, I., Trombettoni, G., Neveu, B.: A contractor based on convex interval taylor. In: Integration of AI and OR Techniques in Contraint Programming for Combinatorial Optimzation Problems, pp. 1–16. Springer (2012)
 [4] Araya, I., Trombettoni, G., Neveu, B., Chabert, G.: Upper bounding in inner regions for global optimization under inequality constraints. Journal of Global Optimization 60(2), 145–164 (2014)
 [5] Belotti, P., Lee, J., Liberti, L., Margot, F., Wächter, A.: Branching and bounds tighteningtechniques for nonconvex minlp. Optimization Methods & Software 24(45), 597–634 (2009)

[6]
Benhamou, F., Goualard, F., Granvilliers, L., Puget, J.F.: Revising hull and box consistency. In: International Conference on Logic Programming. pp. 230–244. MIT press (1999)

[7]
Blum, C., Puchinger, J., Raidl, G.R., Roli, A.: Hybrid metaheuristics in combinatorial optimization: A survey. Applied Soft Computing 11(6), 4135–4151 (2011)
 [8] Chabert, G., Jaulin, L.: Contractor programming. Artificial Intelligence 173, 1079–1100 (2009)
 [9] Cotta, C., Troya, J.M.: Embedding branch and bound within evolutionary algorithms. Applied Intelligence 18(2), 137–153 (2003)
 [10] Csendes, T., Ratz, D.: Subdivision direction selection in interval methods for global optimization. SIAM Journal on Numerical Analysis 34(3), 922–938 (1997)
 [11] Focacci, F., Laburthe, F., Lodi, A.: Local search and constraint programming. In: Handbook of metaheuristics, pp. 369–403. Springer (2003)
 [12] Gallardo, J.E., Cotta, C., Fernández, A.J.: On the hybridization of memetic algorithms with branchandbound techniques. Systems, Man, and Cybernetics, Part B: Cybernetics, IEEE Transactions on 37(1), 77–83 (2007)
 [13] Gropp, W., Moré, J.: Optimization environments and the NEOS server. Approximation theory and optimization pp. 167–182 (1997)
 [14] Hansen, E.: Global optimization using interval analysis. Dekker (1992)
 [15] Kearfott, R.B.: Rigorous global search: continuous problems. Springer (1996)
 [16] Leroy, X., Doligez, D., Frisch, A., Garrigue, J., Rémy, D., Vouillon, J.: The objective caml system release 3.12. Documentation and user’s manual. INRIA (2010)
 [17] Lhomme, O.: Consistency techniques for numeric csps. In: IJCAI. vol. 93, pp. 232–238. Citeseer (1993)
 [18] Mimram, S.: ocamlglpk. http://ocamlglpk.sourceforge.net/ (2004)
 [19] Moore, R.E.: Interval Analysis. PrenticeHall (1966)
 [20] Moore, R.E.: On computing the range of a rational function of n variables over a bounded region. Computing 16(1), 1–15 (1976)
 [21] Neumaier, A., Shcherbina, O.: Safe bounds in linear and mixedinteger linear programming. Mathematical Programming 99(2), 283–296 (2004)
 [22] Ninin, J., Hansen, P., Messine, F.: A reliable affine relaxation method for global optimization. Groupe d’études et de recherche en analyse des décisions (2010)
 [23] Price, K., Storn, R., Lampinen, J.: Differential Evolution  A Practical Approach to Global Optimization. Natural Computing, SpringerVerlag (2006)

[24]
Puchinger, J., Raidl, G.R.: Combining metaheuristics and exact algorithms in combinatorial optimization: A survey and classification. In: Artificial intelligence and knowledge engineering applications: a bioinspired approach, pp. 41–53. Springer (2005)
 [25] Sahinidis, N.V.: Baron: A general purpose global optimization software package. Journal of Global Optimization 8(2), 201–205 (1996)
 [26] Schichl, H., Neumaier, A.: Interval analysis on directed acyclic graphs for global optimization. Journal of Global Optimization 33(4), 541–562 (2005)
 [27] Skelboe, S.: Computation of rational interval functions. BIT Numerical Mathematics 14(1), 87–95 (1974)
 [28] Sotiropoulos, D., Stavropoulos, E., Vrahatis, M.: A new hybrid genetic algorithm for global optimization. Nonlinear Analysis: Theory, Methods & Applications 30(7), 4529–4538 (1997)
 [29] Storn, R., Price, K.: Differential evolution  a simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization pp. 341–359 (1997)
 [30] Trombettoni, G., Araya, I., Neveu, B., Chabert, G.: Inner regions and interval linearizations for global optimization. In: AAAI (2011)
 [31] Trombettoni, G., Chabert, G.: Constructive interval disjunction. In: Principles and Practice of Constraint Programming–CP 2007, pp. 635–650. Springer (2007)
 [32] Van Hentenryck, P.: Numerica: a modeling language for global optimization. MIT press (1997)
 [33] Vanaret, C., Gotteland, J.B., Durand, N., Alliot, J.M.: Preventing premature convergence and proving the optimality in evolutionary algorithms. In: International Conference on Artificial Evolution (EA2013). pp. 84–94 (2013)
 [34] Zhang, X., Liu, S.: A new intervalgenetic algorithm. In: Natural Computation, 2007. ICNC 2007. Third International Conference on. vol. 4, pp. 193–197. IEEE (2007)