Domain reduction techniques for global NLP and MINLP optimization

TL;DR

This paper surveys domain reduction techniques used in NLP and MINLP optimization, analyzing their impact on solver performance through computational experiments on a large test set.

Contribution

It provides a comprehensive review of domain reduction methods and evaluates their effectiveness across multiple global solvers on extensive test problems.

Findings

Domain reduction techniques significantly improve solver convergence.

Performance gains vary depending on problem characteristics.

The survey highlights promising methods for future research.

Abstract

Optimization solvers routinely utilize presolve techniques, including model simplification, reformulation and domain reduction techniques. Domain reduction techniques are especially important in speeding up convergence to the global optimum for challenging nonconvex nonlinear programming (NLP) and mixed-integer nonlinear programming (MINLP) optimization problems. In this work, we survey the various techniques used for domain reduction of NLP and MINLP optimization problems. We also present a computational analysis of the impact of these techniques on the performance of various widely available global solvers on a collection of 1740 test problems.

Tables

Table 1. Table 1 : Test set statistics

Test Library	Global	Princeton	CMU-IBM	MINLP
Number of problems	369	980	142	249
Avg. no. of continuous variables	1092	1355	369	346
Avg. no. of binary variables	0	0	139	235
Avg. no. of integer variables	0	0	0	24
Avg. no. of constraints	785	836	956	534

Table 2. Table 2 : BARON options used for tests

Option	Value
LBTTDo	0
MDo	0
OBTTDo	0
PDo	0
TDo	0

Table 3. Table 3 : Performance deterioration for BARON when domain reduction is turned off

Test Library	Global	Princeton	CMU-IBM	MINLP
% increase in number of nodes	1180	261	546	802
% increase in computational time	67	70	75	47

Table 4. Table 4 : Couenne options used for tests

Option	Value
aggressive_fbbt	no
branch_fbbt	no
feasibility_bt	no
max_fbbt_iter	0
optimality_bt	no
redcost_bt	no
fixpoint_bt	no
two_implied_bt	no

Table 5. Table 5 : Performance deterioration for Couenne when domain reduction is turned off

Test Library	Global	Princeton	CMU-IBM	MINLP
% increase in number of nodes	129	21	186	171
% increase in computational time	26	19	12	16

Table 6. Table 6 : SCIP options used for tests

Option	Value
presolving/maxrounds	0
presolving/components/maxrounds	0
presolving/convertinttobin/maxrounds	0
presolving/domcol/maxrounds	0
presolving/dualagg/maxrounds	0
presolving/dualinfer/maxrounds	0
presolving/gateextraction/maxrounds	0
presolving/implfree/maxrounds	0
presolving/implics/maxrounds	0
presolving/inttobinary/maxrounds	0
presolving/redvub/maxrounds	0
presolving/stuffing/maxrounds	0
presolving/trivial/maxrounds	0
presolving/tworowbnd/maxrounds	0
propagating/maxrounds	0
propagating/dualfix/maxprerounds	0
propagating/genvbounds/maxprerounds	0
propagating/obbt/freq	-1
propagating/obbt/maxprerounds	0
propagating/obbt/tightintboundsprobing	False
propagating/probing/maxprerounds	0
propagating/pseudoobj/maxprerounds	0
propagating/redcost/maxprerounds	0
propagating/rootredcost/maxprerounds	0
propagating/vbounds/maxprerounds	0
conflict/useprop	False
conflict/useinflp	False
conflict/usepseudo	False
heuristics/bound/maxproprounds	0
heuristics/clique/maxproprounds	0
heuristics/randrounding/maxproprounds	0
heuristics/shiftandpropagate/freq	-1
heuristics/shifting/freq	-1
heuristics/vbounds/maxproprounds	0
misc/allowdualreds	False
misc/allowobjprop	False

Table 7. Table 7 : Performance deterioration for SCIP when domain reduction is turned off

Test Library	Global	Princeton	CMU-IBM	MINLP
% increase in number of nodes	174	152	417	56
% increase in computational time	24	33	141	18

Equations

$$\left.\begin{array}[]{lll}\min&f(\vec{x})\\ \mathop{\rm s.t.}&\vec{g(\vec{x})}\leq 0\\ &\vec{x}_{l}\leq\vec{x}\leq\vec{x}_{u}\\ &\vec{x}\in\mathbb{R}^{n-m}\times\mathbb{Z}^{m}\end{array}\right\}$$

$$\left.\begin{array}[]{lll}\min&3x+4y\\ &x-y\leq 4\\ &xy\leq 3\\ &x^{2}+y^{2}\geq 1\\ &1\leq x\leq 5\\ &1\leq y\leq 5\\ \end{array}\right\}$$

$$\min$$

$$x-y\le 4$$

$$z_3\le 3$$

$$z_1+z_2\ge 1$$

$$z_1=x^2$$

$$z_2=y^2$$

$$z_3=xy$$

$$1\le x\le 5$$

$$1\le y\le 5$$

$$\vec{x}+\vec{y}=[{\vec{x}}^l+{\vec{y}}^l,{\vec{x}}^u+{\vec{y}}^u]$$

$$\vec{x}+\vec{y}=\{x+y|x\in \vec{x}\text{ and }y\in \vec{y}\}$$

$$\vec{x}-\vec{y}$$

$$\vec{x}\times \vec{y}$$

$$\frac{1}{\vec{x}}$$

$$\frac{\vec{x}}{\vec{y}}$$

$$\left.\begin{array}[]{lll}\min&\pm x_{k}\\ \mathop{\rm s.t.}&\vec{g(\vec{x})}\leq 0\\ &\vec{x}^{l}\leq\vec{x}\leq\vec{x}^{u}\\ &\vec{x}\in\mathbb{R}^{n-m}\times\mathbb{Z}^{m}\end{array}\right\}$$

$$\left.\begin{array}[]{lll}\min&\pm x_{k}\\ \mathop{\rm s.t.}&\vec{g}_{\rm conv}(\vec{x})\leq 0\\ &\vec{x}^{l}\leq\vec{x}\leq\vec{x}^{u}\\ &\vec{x}\in\mathbb{R}^{n}\end{array}\right\}$$

$$\left.\begin{array}[]{lll}\min&\pm x_{k}\\ \mathop{\rm s.t.}&\vec{g}_{\rm lin}(\vec{x})\leq 0\\ &\vec{x}^{l}\leq\vec{x}\leq\vec{x}^{u}\\ &\vec{x}\in\mathbb{R}^{n}\end{array}\right\}$$

$$\left.\begin{array}[]{lll}\min&\pm\vec{a}_{i}^{T}\vec{x}\\ \mathop{\rm s.t.}&\vec{a}_{j}^{T}\vec{x}\leq b_{j}\qquad j=1,\ldots,i-1,i+1,\ldots,n\\ &\vec{x}^{l}\leq\vec{x}\leq\vec{x}^{u}\\ &\vec{x}\in\mathbb{R}^{n}\end{array}\right\}$$

$$\sum _{j=1}^n{a}_{ij}{x}_j\le$$

$${\vec{x}}^l\le \vec{x}\le$$

$$x_h\le \frac{1}{a_{ih}}(b_i-\sum _{j\ne h}\min \left(a_{ij}{x}_j^U,a_{ij}{x}_j^L\right)),a_{ih}>0$$

$$x_h\ge \frac{1}{a_{ih}}(b_i-\sum _{j\ne h}\min \left(a_{ij}{x}_j^U,a_{ij}{x}_j^L\right)),a_{ih}<0$$

$$x_1+x_2$$

$$x_2+x_3$$

$$x_1\in$$

$$x_2\in$$

$$x_3\in$$

$$x_1+x_2$$

$$-x_1+x_2$$

$$x_1\in$$

$$x_2\in$$

$$x_i\le \max \{x_j^L{x}_k^L,x_j^L{x}_k^U,x_j^U{x}_k^L,x_j^U{x}_k^U\}$$

$$x_i\ge \min \{x_j^L{x}_k^L,x_j^L{x}_k^U,x_j^U{x}_k^L,x_j^U{x}_k^U\}$$

$$x_j\le v+(x_j^u-v)x_i$$

$$x_j\ge v-(v-x_j^l)x_i$$

$$-230x_{10}-200x_{16}-400x_{17}\le -5$$

Full text

Keywords: Constraint propagation; feasibility-based bounds tightening; optimality-based bounds tightening; domain reduction

Introduction

We consider the following mixed-integer nonlinear programming problem (MINLP):

[TABLE]

MINLP is a very general representation for optimization problems and includes linear programming (LP), mixed-integer linear programming (MIP) and nonlinear programming (NLP) in its subclasses. A variety of applications in diverse fields are routinely formulated using this framework including water network design [94, 62], hydro energy systems management [44], protein folding [143], robust control [18], trim loss [84], heat exchanger network synthesis [69], gas networks [128], transportation [68], chemical process operations [76], chemical process synthesis [77], crystallographic imaging [178], and seizure predictions [148]. Modelling via nonconvex objective functions or constraints is necessitated for many of these practical applications. Aside from the combinatorial complexity introduced by the integer variables, nonconvexities in objective function or the feasible region lead to multiple local minima and provide a challenge to the optimization of such problems.

Branch-and-bound [138] based methods can be exploited to solve Problem 1 to global optimality. Inspired by branch-and-bound for discrete programming problems [109], branch-and-bound was adapted for continuous problems by Falk and Soland [58]. The algorithm proceeds by bounding the global optimum by a valid lower and upper bound throughout the search process. Whenever these bounds are within an acceptable tolerance, the algorithm terminates with the upper bound as a (near-)global optimum. The algorithm exhaustively searches the space by branching on variables to divide the space into subdomains. The lower and upper bounding procedures are recursively applied to the new subdomains in a tree search until the bounds converge. Branch-and-bound algorithms are known to terminate within $\epsilon$-accuracy ($\epsilon>0$) to a global optimum as long as branching, node selection and lower bounding are designed to satisfy certain conditions [89]. For special cases, the true global optimum can be finitely achieved with branch-and-bound [173, 8, 38]. The success of branch-and-bound methods for global optimization is evident from the numerous software implementations available, including ANTIGONE [135], BARON [185], Couenne [25], LindoGlobal [117] and SCIP [3].

In this work, we survey the various domain reduction techniques that are employed within branch-and-bound algorithms. While these techniques are not necessary to ensure convergence to the global optimum, they typically speed up convergence. These techniques often exploit feasibility analysis to eliminate infeasible parts of the search space. Alternatively, the methods can also utilize optimality arguments to shrink the search space while ensuring at least one optimal solution is retained. Domain reduction techniques constitute the major component of the solution methods for satisfiability problems through unit propagation [126] and for constraint programming (CP) through various filtering algorithms that achieve differing levels of consistencies [29]. They are also exploited in artificial intelligence (AI) [54] and interval analysis [142, 98]. Some of the other names used in the literature for these methods include bound propagation, bounds tightening, bound strengthening, domain filtering, bound reduction and range reduction.

Mathematical-programming-based methods for solving nonconvex MINLPs often rely on the solution of a relaxation problem for finding a valid lower bound. The strength of the relaxations employed for lower bounding depends on the diameter of the feasible region. Smaller domains lead to tighter relaxations. For example, Figure 1 shows the convex relaxation for a simple univariate concave function. A convex relaxation for a given function defined on a nonempty convex set is a convex function that underestimates the given function on its domain. The convex relaxation over the reduced domain provides a better approximation for the univariate concave function thereby providing better lower bounds. Domain reduction techniques not only reduce the diameter of the search space, but they also improve the tightness of convex relaxations.

Domain reduction techniques have been extensively studied in various communities including AI, CP, mathematical programming, interval analysis and computer science where they can be viewed as chaotic iterations [14]. The primary objective of this paper is to review key domain reduction techniques applicable to (nonconvex) NLPs and MINLPs and point to connections between methods from constraint programming and interval arithmetic communities wherever applicable. We also present results on standard test libraries with different global solvers to demonstrate the impact of domain reduction strategies on their performance.

The remainder of the paper is organized as follows. Representation of general MINLPs through factorable reformulations and directed acyclic graphs is described in Section 2. Methods for constraint propagation and bounds tightening in global optimization of MINLPs often rely on interval arithmetic. A brief introduction to this topic is provided in Section 3. We introduce presolving for optimization in Section 4. Domain reduction techniques that rely on eliminating infeasible regions for the problem are described in Section 5. Techniques that utilize optimality arguments for carrying out domain reduction are described in Section 6. Computational tradeoffs in the implementation of some of the advanced techniques are discussed in Section 7. The computational impact of many of these techniques on the performance of widely available solvers is investigated in Section 8. Finally, we conclude in Section 9.

Representation

A crucial step in the branch-and-bound algorithm is the construction of relaxations. One of the most widely used methods for this purpose is the idea of factorable reformulations. It involves splitting a problem into basic atomic functions that are utilized for computing the function, an idea exploited by McCormick [131], who developed a technique that constructed non-differentiable relaxations of optimization problems. Ryoo and Sahinidis [159] gain differentiability by introducing new variables and equations for each of the intermediate functional forms. These functional forms are simple in nature like the bilinear term. A similar idea was proposed by Kearfott [97] where new variables and equations are introduced to decompose nonlinearities to allow for more accurate computations of interval Jacobian matrices.

Consider the following example:

[TABLE]

A factorable reformulation can proceed by introducing a new variable for every nonlinearity occurring in the model. We replace $z_{1}$ for $x^{2}$, $z_{2}$ for $y^{2}$ and $z_{3}$ for $xy$.

A factorable reformulation of the model is thus given by:

[TABLE]

It is trivial to outer approximate the univariate terms of this model (cf. Figure 1), while the bilinear term in 4 may be outer approximated by its convex and concave envelopes [131]. The combination of these outer approximators provides a convex relaxation of the original problem. In general, factorable reformulations decompose the problem into simpler functional forms for which convex relaxations are known. Thus, a convex relaxation can be obtained by reformulating a problem into its factorable form and relaxing each of the simple functional forms present in the model. For example, see Algorithm Relax f(x) in Tawarmalani and Sahinidis [185].

Factorable reformulations can be conveniently represented through a directed acyclic graph [179, 167]. In this graph representation, variables ($x,y,z$) and constants are leaf nodes, vertices are elementary operations ($+,-,*,/$, $\log$, $\exp$, etc.) and the functions to be represented are the root nodes. Variable and constraint bounds are represented through suitable intervals for the root and leaf nodes. Common subexpressions are combined to reduce the size of the graph as doing so is known to tighten the resulting relaxations [185, Theorem 2]. Different mathematical formulations can be generated from the same DAG depending on the needs of the solver. Expressions are evaluated by propagating values from the leaves to the root node through the edges in a forward mode. Backward propagation is utilized for the computation of slopes and derivatives. Slopes can also be utilized to construct linear relaxations for the problem. Figure 2 represents the DAG for Problem 2.

Interval arithmetic

Interval arithmetic is a system of arithmetic based on intervals of real numbers [137]. An interval variable is defined using a variable’s lower and upper bounds; the variable itself is restricted to lie between the bounds. Consider the interval variables $\vec{x}=[\vec{x}^{l},\vec{x}^{u}]$ and $\vec{y}=[\vec{y}^{l},\vec{y}^{u}]$. Addition on two intervals can be defined as:

[TABLE]

The addition operator defined by equation 5 has the following property:

[TABLE]

Extensions to the definitions of other classic operators can be similarly defined:

[TABLE]

These operators can be suitably defined to account for infinities within the intervals and for the case when 0 lies in the interval of the denominator in a division operation [107]. Natural extensions of factorable functions can be computed by replacing all the elementary operations involved in the computation of a factorable function with their interval counterparts. A natural extension provides valid lower and upper bounds for the value of a function.

Floating point arithmetic is susceptible to rounding errors, which can lead to solutions being lost even for simple problems. Neumaier and Shcherbina [145] provide an example of a simple MIP problem that leads to incorrect solutions by major commercial solvers due to roundoff errors. Interval arithmetic methods utilize outward rounding to ensure no points are lost due to roundoff errors. These methods are utilized to design rigorous branch-and-bound-based optimization and root finding methods that ensure no solutions are lost due to roundoff errors [193, 15]. The methods bound function values through interval arithmetic and carry out domain reduction through branching and fathoming tests based on monotonicity, convexity, infeasibility as well as interval Newton type methods which provide conditions for existence and uniqueness of solutions within a box [81, 142]. Interval-based methods have been used for practical applications like solvent blend design [177]. Interval-based solvers include ICOS [110], Globsol [99] and Numerica [193]. Interval-based solvers are often slower than their nonrigorous counterparts. More recently, aspects of safe computations are being introduced with different parts of nonrigorous optimization algorithms. For example, Neumaier and Shcherbina [145] describe methods to guarantee safety in the solution of MIPs through suitable preprocessing of the LP relaxations and postprocessing of their solutions. Their methods lead to valid results even when the LP solver itself does not use rigorous methods for obtaining a solution. Borradaile and van Hentenryck [34] provide ways of constructing numerically safe linear underestimators for univariate and multivariate functions. Numerically safe methods for computation have also been developed by the CP community, and are referenced in the remainder of the paper including in Sections 5 and 6.

Presolving optimization models

The idea of analyzing and converting an optimization model into a form more amenable to fast solution is old. Starting with the work of Brearley et al. [35], a number of techniques have been used for analyzing models in the operations research community. We use the term presolve to denote all the techniques used for simplification of optimization models. These techniques have been developed extensively for linear programming models [189, 190, 93, 74] and for mixed-integer linear programming models [157, 87, 141, 164, 121]. Bixby and Rothberg [31] observe that turning off root node presolve at the start of a branch-and-bound search degrades performance of CPLEX 8.0 by a factor of 10.8, while turning off presolve at every node other than the root node degrades the performance by a factor of 1.3 on certain MIP models, demonstrating the importance of presolve. See Achterberg and Wunderling [6] for an extensive computational analysis of the impact of various components of presolve algorithms for MIPs. Some of the ideas for simplification of models include:

• •

Elimination of redundant constraints

• •

Identification and elimination of dominated constraints (dominated constraints are constraints with a feasible region that is a superset of the feasible region of other constraints in the model)

• •

Elimination of redundant variables

• •

Assimilating singleton rows into bounds on variables

• •

Tightening bounds on dual variables

• •

Fixing variables at their bounds

• •

Increasing sparsity in the model

• •

Rearrangement of variables and constraints to induce structure

Similar preprocessing operations can be derived for nonlinear problems [10, 119, 63, 133, 75]. Some of the general guidelines include:

• •

Avoid potentially undefined functions

• •

Reduce nonlinearity in the model

• •

Improve scaling of model

• •

Increase convexity in the model through reformulations

Amarger et al. [10] provide a software implementation REFORM to carry out many of these reformulations. Presolve techniques are usually implemented at the solver level by developers of various software for optimization. The modelling systems AIMMS [7] and AMPL [11] also provide dedicated presolve systems that are applied to all optimization models [67, 91]. The success of presolve in mathematical programming has led to efforts for its extension to general constraint programming such as automated reformulation of CP models [115, 147]. These methods can lead to considerable simplifications in the model, and can also lead to a reduction in the memory required to solve the problem. Another advantage of these presolve methods is that they are often able to detect infeasibility in optimization models. If a presolved model is infeasible, then the original model is also infeasible. Presolve methods can also be used to detect and correct the causes of infeasibilities for infeasible optimization models. Chinneck [46] provides examples of simple cases where infeasibilities can be correctly diagnosed by analyzing the sequence of reductions obtained with presolve. More recently, Puranik and Sahinidis [152] provide an automated infeasibility diagnosis methodology through the isolation of irreducible inconsistent sets (IISs). An IIS is defined as an infeasible set of constraints that has every subset feasible. Isolating an IIS can help accelerate the process of model correction by allowing the model expert to focus onto a smaller problem area within the model. The authors of [152] propose a deletion presolve procedure that exploits feasibility-based bounds tightening techniques in order to accelerate the isolation of an IIS.

While simplified models are often easier to solve than their original counterparts, once an optimal solution has been obtained, the modelling system must return a solution that can be interpreted by the user with respect to the original model form that was specified. Restoration procedures to obtain primal and dual solutions to the original problem are described by Andersen and Andersen [12] and Fourer and Gay [67]. Such procedures are not discussed here. In the subsequent sections, we review domain reduction techniques that are a major component of all presolve algorithms.

Reduction of infeasible domains

The methods described in this section eliminate regions from the search space where no feasible points of Problem 1 can exist. Global optimization algorithms maintain continuous and discrete variable domains through upper and lower bounds. Domain reduction is achieved by making these bounds tighter. Therefore, domain reduction is also commonly referred to as bounds tightening.

The tightest possible bounds based on feasibility of the constraints of Problem 1 can be obtained by solving the following problems for each of the $n$ variables ($k=1,\ldots,n$):

[TABLE]

$\pm x_{k}$ in problem 11 denotes two optimization problems, where $x_{k}$ and $-x_{k}$ are individually minimized. Solution of these optimization problems returns the tightest possible bounds on the feasible region. If these bounds are tighter than the user-specified bounds in the model, we can achieve domain reduction by using them. However, since the constraints of Problem 11 are potentially nonconvex and due to the presence of integer variables, these are expensive global optimization problems, which might be as hard to solve as Problem 1. A computationally inexpensive way as compared to problem 11 to obtain valid bounds for each of the variables is by solving the following problems for each of the $n$ variables ($k=1,\ldots,n$):

[TABLE]

where $\vec{g}_{\rm conv}(\vec{x})\leq 0$ refers to a convex relaxation of the constraints. The convex relaxation can be nonlinear. The integrality restrictions on the variables are relaxed as well. While bounds obtained from Problem 12 in general are weaker than the bounds obtained from Problem 11, they require the solution of convex optimization problems and are therefore obtained more efficiently. Convex relaxations are utilized in branch-and-bound methods for obtaining valid lower bounds to the optimum, and are thus already available for bounds tightening. To exploit the efficiency and robustness of LP solvers, these convex relaxations are linearized through outer approximation to obtain linear relaxations [186]. Linear relaxations may provide weaker bounds than nonlinear ones, but can be solved very efficiently. This linearization can also be utilized for obtaining tighter bounds on variables through the solution of the following problems for each of the $n$ variables ($k=1,\ldots,n$):

[TABLE]

where $\vec{g}_{\rm lin}(\vec{x})\leq 0$ refers to a linearized outer approximation of the feasible region. The bounds obtained from Problem 13 are in general weaker than the bounds obtained from Problem 12, but require the solution of linear instead of nonlinear programming problems and are therefore solved more efficiently.

Feasibility-based arguments for reduction can also be utilized for tightening constraints. Consider a linear set of constraints $\vec{Ax}\leq\vec{b}$. Tight bounds on constraint $i$ can be obtained with the solution of following optimization problems [164]:

[TABLE]

Problem 14 can help identify redundancy if the maximum value of $\vec{a}_{i}^{T}\vec{x}$ is strictly less than $b_{i}$. Conversely, if the minimum value of $\vec{a}_{i}^{T}\vec{x}$ is strictly greater than $b_{i}$, this also identifies infeasibility.

In general, full solution of LP or NLP problems for bounds tightening can be expensive. To balance the computational effort involved with the reduction in bounds obtained, these techniques are usually carried out only at the root node and/or for a subset of variables. They are utilized only sparingly through the rest of the search or not at all. In some cases, however, solving optimization problems for tightening methods throughout the search has been shown to provide significant computational benefits [43].

Note that reduction methods exploiting Problems 11, 12 or 13 are often referred to in the literature as optimality-based bounds tightening. While these methods utilize the solution of optimization problems, they only carry out reduction of infeasible regions from the search space. For this reason, we prefer to classify them under feasibility-based reduction techniques.

Bounds propagation techniques

Propagation-based bounds tightening techniques will be simply referred to as propagation in the remainder of the paper. These techniques find their roots in several works in the literatures of mathematical programming [35], constraint logic programming [47], interval arithmetic [50], and AI [54]. These methods are often referred to as bounds propagation techniques in CP and are a specialization of constraint propagation. Davis [54] refers to a constraint network with nodes (variables) which can take possible labels (domains) and are connected by constraints. He further refers to six different categories of constraint propagation based on the type of information which is updated:

• •

Constraint inference: New constraints are inferred and added.

• •

Label inference: Constraints are utilized to restrict the sets of possible values for nodes.

• •

Value inference: Nodes are partially initialized and constraints are utilized to complete assignments for all nodes.

• •

Expression inference: Nodes are labelled with values expressed over other nodes.

• •

Relaxation: Nodes are assigned exact values which may violate certain constraints.

• •

Relaxation labelling: Nodes are assigned labels using probabilities. Updates in the network involve update of the probabilities.

Propagation is equivalent to label inference in the AI community. Davis utilizes the Waltz algorithm [199] to describe bounds propagation in a constraint network. In the CP community, these methods are often utilized for solving discrete constraint satisfaction problems (CSP) [88] with backtracking-based search methods [192]. However, they are also utilized for continuous CSPs [28, 54].

Hager [79] discusses a reduce operator to eliminate regions outside the solution set for solving systems of constraints. These methods have also been developed extensively in the mathematical programming community, first for LPs, and later for nonlinear programming problems [187, 82, 80, 108, 132]. Mclinden and Mangasarian [122] demonstrate inference of bounds for simple monotonic complementarity problems and convex problems. Lodwick [118] analyzed the relationship between the constraint propagation methods from AI and bound tightening methods from mathematical programming. These methods typically operate by systematically analyzing one constraint at a time to infer valid bounds on variables.

Propagation techniques are computationally inexpensive. For example, consider a set of linear constraints:

[TABLE]

The following inequalities are implied by every linear constraint:

[TABLE]

[TABLE]

If these inequalities imply tighter bounds on $x_{h}$ than the ones specified by the model, the bounds can be updated. For example, consider the set of inequalities:

[TABLE]

From inequality 17, we can infer $x_{1}\geq 4-\max(0,4)=0$. Thus, the lower bound of $x_{1}$ is updated to 0. From inequality 18, we can infer that $x_{2}\leq 1-\min(-1,1)=2$. Thus, the upper bound of $x_{2}$ is updated to 2. By analyzing inequality 17 again, we can update the lower bound of $x_{1}$ to 2. The domains of the variables after these bounds tightening steps are $x_{1}\in[2,4],x_{2}\in[0,2]\text{ and }x_{3}\in[-1,1]$. Harvey and Schimpf [85] describe how bounds can be iteratively tightened in sublinear time for long linear constraints with many variables.

Inequalities 15 and 16 indicate that only one of the bounds for a variable can be updated from an inequality constraint based on the sign of its coefficient. Achterberg [2] formalizes this through the concept of variable locks. Thus, inequality 17 down locks variables $x_{1}$ and $x_{2}$, since the lower bounds for $x_{1}$ and $x_{2}$ cannot be moved arbitrarily without violating inequality 17. The concept of variable locks allows for efficient duality fixing of variables. Note that there is no up lock for variable $x_{1}$ and no down lock for variable $x_{3}$ based on the inequalities 15 and 16. Thus, if the coefficient of $x_{1}$ in the objective function is negative, $x_{1}$ can be set to its upper bound. Similarly, if the coefficient of $x_{3}$ in the objective function is positive, $x_{3}$ can be set to its lower bound. The concept of variables locks was extended for CP to develop presolve procedures for cumulative constraints by Heinz et al. [86]. Duality fixing is extended by Gamrath et al. [70] to allow for fixing of a singleton column. The authors also define dominance between two variables for a MIP and show how dominance information can be used for fixing variables at their bounds.

Figure 3 considers three cases of propagation. In Case 1, bounds for $x$ and $y$ are successfully tightened through iterative application of propagation, with no further tightening possible via Problem 11. In Case 2, the bound obtained for $x$ by propagation is the tightest, however the bound for $y$ can be tightened further by Problem 11. In Case 3, bounds for neither $x$ nor $y$ can be tightened by propagation, whereas significant reduction is possible via Problem 11.

In general, reduction in domain through propagation is not guaranteed. Consider the following example:

[TABLE]

The tightest possible domain for $x_{1}$ and $x_{2}$ based on feasibility arguments is $[-1,3]$ and can be obtained by solving Problem 11. However, analyzing constraints one-at-a-time leads to no reduction of bounds for this problem. The reason for this drawback is that propagation only analyzes one constraint at a time, whereas Problem 11 utilizes information from all the constraints in the model simultaneously.

Belotti [22] considers bound reduction by generating a convex combination of two linear inequalities and utilizing this combination for bound reduction. The author proposes a univariate optimization problem to determine the optimal value for the convex multipliers. However, the computational complexity of considering all possible combinations of $m$ constraints in $n$ variables is $O(m^{2}n^{3})$. The author proposes a heuristic scheme which, when implemented only on a subset of the nodes of the branch-and-bound tree, shows computational benefits of using this method. For other presolve-based methods for mixed-integer programs analyzing more than one constraint at a time, see Achterberg et al. [4]. Domes and Neumaier [55] propose the generation of a new constraint by taking the linear combination or aggregation of all constraints for quadratic problems. This new constraint can uncover new relationships between the variables and can lead to domain reduction. The authors propose the use of the dual multipliers of a local solution in the case of a feasible subproblem and the use of a constraint violation measure in the case of an infeasible subproblem to aggregate constraints. Their method shows computational benefits in reducing the cluster effect (Section 6.4). The notion of global constraints in CP [154] is related. A global constraint provides a concise representation for a set of individual constraints. Filtering algorithms on a global constraint effectively carry out domain reduction by considering multiple individual constraints simultaneously. Lebbah et al. [111] present a global constraint for a set of quadratic constraints and describe filtering algorithms. Similar ideas have been extended to polynomial constraints [113].

Bounds can also be inferred for nonlinear constraints via propagation. Consider a bilinear term of the form $x_{i}=x_{j}x_{k}$. Then,

[TABLE]

and

[TABLE]

represent valid bounds for variable $x_{i}$ and can be used to potentially tighten any user-specified bounds for this variable.

In the context of factorable reformulations, an optimization problem is already decomposed into simple functional forms, which can be readily utilized for domain propagation. Consider, for instance, the factorable reformulation for Problem 2. From constraint 3, we can infer that $0\leq z_{3}\leq 3$. From constraint 4, we can infer that $1\leq y\leq 3$.

The iterative application of simple constraints such as 15–16 and 19–20 is referred to as poor man’s LPs and poor man’s NLPs [173, 160] because this iterative application is an inexpensive way to approximate the solution of linear and nonlinear optimization problems that aim to reduce domains of variables. Propagation based on these techniques can be applied not only to the original problem constraints but also to any cutting planes and other valid inequalities that might be derived by the branch-and-bound solution algorithm. Moreover, all these techniques can be implemented efficiently via the DAG representation. If a bound on a variable $x$ has changed, it can be propagated to other variables that depend on $x$ through a simple forward propagation on the DAG based on interval arithmetic. Similarly, the variables on which $x$ depends can be updated through a backward propagation.

Convergence of propagation

Propagation methods can be carried out iteratively as long as there is an improvement in variable bounds. However, these methods can fail to reach a fixed point finitely. Consider as an example [144]: $x_{1}+x_{2}=0$, $x_{1}-qx_{2}=0$ with $q\in(0,1)$. Propagation will converge to $(0,0)$ at the limit, i.e., as the number of iterations approaches infinity. On the contrary, any linear programming based method will terminate at $(0,0)$ quickly as the only feasible point. A necessary condition for nonconvergence of propagation iterations is the presence of cycles in expression graphs [59]. The problem of achieving a fixed point with propagation iterations in the presence of only integer variables has been shown to be NP-complete [32, 33]. The fixed point obtained is independent of the order in which the variables are considered [41]. In practice, the iterations are usually terminated when the improvement in variable domains is insignificant or when an upper limit on the number of iterations is reached.

For linear inequalities in the presence of continuous variables alone, the fixed point can be obtained in polynomial time by the solution of a large linear program [23]. For nonlinear problems and for problems with integer variables, a linear relaxation can be solved to obtain reduced bounds [24].

Consistency

Constraint programming algorithms for constraint satisfaction problems rely on the formalized notion of consistency to propagate domains along constraint networks and remove values from variable domains that cannot be a part of any solution. Differing levels of consistency exist, including arc consistency and k-consistency, while various algorithms have been proposed to achieve these consistencies [29]. These concepts were originally proposed for discrete problems but were also extended to continuous constraint satisfaction problems [116, 27, 162, 26, 49, 30].

The domain reduction algorithms typically employed for global optimization of MINLPs are usually not iterated until they reach a certain level of consistency because attempting to establish consistency can lead to a prohibitively large computational effort. However, domain reduction techniques are not even necessary to prove convergence of branch-and-bound based methods for global optimization. In contrast to CP methods, the availability of relaxation methods for generating bounds allows the mathematical-programming-based branch-and-bound algorithms to converge without establishing any levels of consistency.

Techniques for mixed-integer linear programs

Many techniques for presolving LPs can also be directly utilized for MIPs. However, specialized reduction methods can be developed for MIPs. For example, a number of methods have been developed for the analysis of 0-1 binary programs [78, 51, 164]. Probing and shaving techniques are equivalent to the singleton consistency techniques from CP. Singleton consistency techniques [150] proceed by assigning a fixed value to a variable and carrying out propagation. If this assignment leads to infeasibility of the constraint program, the value assigned to the variable cannot be a part of any solution and can be eliminated. Probing techniques, that can be used for binary programs, involve the fixing of a binary variable $x_{i}$ to say 0. If basic preprocessing and other domain reduction methods are able to prove infeasibility for this subproblem, then the variable $x_{i}$ can be fixed to 1. If fixing the binary variable $x_{i}$ to both 0 and 1 lead to infeasibility, the problem can be declared as infeasible. For binary programs, if probing is carried on all binary variables, it leads to singleton consistency. Probing has also been utilized for satisfiability problems [120]. A similar technique for continuous problems is called shaving [129, 144]. The method proceeds by removing a fraction of the domain of a variable $[x_{i}^{l},x_{i}^{u}]$ as either $[x_{i}^{l}+\epsilon,x_{i}^{u}]$ or $[x_{i}^{l},x_{i}^{u}-\epsilon]$ and testing whether the reduced domain leads to provable infeasibility of the model. If infeasibility is indeed proved, the domain of the variable can be reduced to the complement of the box used for testing. The author suggests using $\epsilon=0.1\times(x_{i}^{u}-x_{i}^{l})$ for this method. Shaving is widely used in the sub-field of CP based scheduling [191]. Belotti et al. [25] call this technique aggressive bounds tightening (ABT) and implement it in the solver COUENNE. Once a fraction of the domain of a variable has been eliminated, propagation is invoked in the hopes of proving infeasibility in the subproblem. However, ABT and shaving techniques in general are expensive and reduction in the domain is not guaranteed. Faria and Bagajewicz [62] propose several variants of this shaving strategy for bilinear terms and utilize it in a branch-and-bound algorithm for water management and pooling problems [61, 60]. Nannicini et al. [140] propose a similar strategy which they refer to as aggressive probing. In contrast to ABT technique, the nonconvex restrictions created by restricting a variable domain in their method are solved to global optimality with branch-and-bound in order to carry out the maximum possible domain reduction.

Implication-based reductions utilize relations between the values of different variables that must be satisfied at an optimal solution. For instance, if fixing two binary variables $x_{i}$ and $x_{j}$ to 0 leads to infeasibility, then we have an implication requiring that one of the variables must be nonzero, which can be represented by the inequality $x_{i}+x_{j}\geq 1$. Such relations can be efficiently represented through the use of conflict graphs [17]. Conflict graphs are utilized to generate valid inequalities that strengthen the MIP formulation. In general, if fixing a binary variable $x_{i}$ to 0 implies that variable $x_{j}$ must take value $v$, then the following inequalities are valid for the problem [121]:

[TABLE]

Implications can be derived by carrying out probing by fixing more than one variable at a time and/or through the analysis of problem structure. Implications can also be utilized for identifying and eliminating redundant constraints. Inequalities derived from implications lead to automatic disaggregation for some constraints [164]. Disaggregated constraints, while redundant for the MIP formulation, lead to tighter LP relaxations for the problem. Achterberg et al. [5] show how implications can be derived from conflicts (See Section 5.5) and from knapsack covers.

Tighter formulations can also be obtained for a problem through coefficient reduction. Consider the example of a linear constraint on binary variables from [127]:

[TABLE]

The coefficients of this constraint can be reduced to:

[TABLE]

While the set of binary values satisfying both of the above constraints are the same, the reduced constraint has a tighter LP relaxation. Approaches for coefficient reduction are provided by [51, 164]. Andersen and Pochet [13] prove that, if no coefficients for an MIP system can be strengthened, then there does not exist a dominating constraint that can be used to replace an existing constraint in the MIP system to tighten its relaxation. They also describe an optimization formulation and an algorithmic solution to the problem of strengthening a coefficient in a constraint as much as possible.

Conflict analysis

Branch-and-bound based algorithms often encounter infeasibility in subproblems during the search for global optimum. Solvers for the satisfiability problem (SAT) also utilize a backtracking-based branching scheme for their solution [139]. The SAT problem consists of binary variables constrained by a set of logical conditions. The SAT problem has a solution if there exists an instantiation of the binary variables satisfying all the logical conditions. SAT-based solvers learn and add conflict clauses from infeasible subproblems [126]. These clauses are created by identifying the corresponding instantiations of a subset of variables leading to infeasibility and prohibiting them. This often leads to a reduction in the search tree. The idea of conflict analysis is similar to the idea of no-good learning from the CP community [180, 114, 96, 100]. No-good is a generalization of a conflict clause from SAT to CP.

The ideas of conflict analysis have been extended for MIP [1, 163]. However, generation of conflict clauses is more complicated for MIP due to the presence of both continuous and general integer variables. Infeasibility in SAT problems and CP problems is detected through a chain of logical deductions caused by fixing some variables. However, in MIP, infeasibility can be identified through either such deductions or an infeasible LP relaxation. Achterberg [1] proposes a generalized conflict graph (termed implication graph in [163]) for representation of bound propagations that can be used to identify cause of infeasibility. Note that the conflict graph and the implication graph are defined differently than in Section 5.4. In the case of an infeasible LP relaxation, Achterberg proposes identifying a minimum bound-cardinality IIS. Representation of conflict causes for MIP requires use of disjunctive constraints which must be reformulated with additional binary variables and inequalities. Limited computational analysis demonstrates a reduction in the number of nodes and time with conflict analysis.

Reduction of suboptimal domains

In contrast to the methods of the previous section, the methods described here can lead to the elimination of feasible points from the domain under the condition that at least one globally optimal solution remains within the search space.

Consider the following convex relaxation of Problem 1:

[TABLE]

As discussed before, convex relaxations are often constructed and linearized in a global branch-and-bound search for obtaining valid lower bounds. Assume that the optimal objective for Problem 21 has value $L$ and, at the optimal solution, a bound $x_{j}\leq x_{j}^{u}$ is active with a Lagrange multiplier $\lambda_{j}>0$. Let $U$ be a known valid upper bound for the optimal objective function value of Problem 1. Then, the following constraint does not exclude any optimal solutions better than $U$ (Theorem 2 in [158]):

[TABLE]

A geometric interpretation of this cut can be observed in Figure 4, where $x_{j}^{*}$ denotes the right-hand-side of equation 22. The constraint excludes values of $x_{j}$ for which the convex relaxation is guaranteed to have its value function to be greater than or equal to $U$. Consequently, the nonconvex problem also has its value function guaranteed to be greater than or equal to $U$ in this domain. A corresponding cut can also be derived if, at the optimal solution $L$ of the convex relaxation, a variable is at its lower bound, i.e., $x_{j}=x_{j}^{l}$, with the corresponding Lagrange multiplier $\lambda^{j}>0$. The upper bound can then be potentially tightened without losing optimal solutions by the following cut:

[TABLE]

For variables that are not at their bounds at the relaxation solution, probing tests can be carried out by temporarily fixing variables at their bounds and solving the restricted problem. For example, Tests 3 and 4 in [158] work as follows. Set $x_{j}=x_{j}^{u}$ and solve Problem 21. If the corresponding multiplier $\lambda^{j}$ is positive, then the following constraint is valid:

[TABLE]

A similar probing test can be developed by fixing a variable at its lower bound [158, 159, 173]. Reduction by 22 and 23 only requires the solution of the relaxation problem, and thus can be implemented at every node of the branch-and-bound tree without much computational overhead. On the other hand, probing with fixing variables at their bounds requires the solution of $2n$ convex problems. Probing is usually carried out only at the root node and then only at some nodes of the tree or only for a subset of the variables. The probing tests are a generalization of probing for the integer programming case. However, they differ from the integer programming case in the sense that probing leads to variable fixing in the case of binary variables, whereas in the continuous case it leads to reduction in the variable domains.

Similarly, if at the optimal solution of the relaxation Problem 21, a constraint $g^{j}_{\rm conv}(\vec{x})\leq 0$ is active with the corresponding multiplier $\mu_{i}\geq 0$, then the following constraint does not violate any points with objective values better than $U$ [159]:

[TABLE]

These inequalities can be appended to the original formulation. However, this process can lead to the accumulation of a large number of constraints. Alternatively, these are often utilized for propagation-based tests locally at the current node and then discarded. Lebbah et al. [112] provide a safe way of implementing duality-based reduction techniques in the context of floating point arithmetic. It must be noted that suboptimal Lagrangian multipliers often provide greater reduction through constraints 22, 23, 24 and 25. A more detailed discussion on the use of other than optimal multipliers is provided in Tawarmalani and Sahinidis [185] and Sellmann [171].

If a valid upper bound U is available, it can also be exploited by reformulating Problem 12 as follows:

[TABLE]

Problem 26 differs from Problem 12 by the addition of a single constraint. In this constraint, $f_{\text{conv}}(\vec{x})$ refers to the convex relaxation of the objective function. It can be replaced by a suitable linearization $f_{\text{lin}}(\vec{x})$ for use in formulation 27.

[TABLE]

Zamora and Grossmann [201] refer to Problem 26 as the contraction subproblem. They define a contraction operation that uses the solution of Problem 26 along with multipliers-based reduction steps described in this section to carry out domain reduction. Their algorithm is therefore referred to as branch-and-contract. Optimality-based arguments can also be used for filtering in CP, see for example [66, 188].

Note that problem 26 can be solved iteratively to obtain further tightening. However, convergence to a fixed point can be slow. Caprara and Locatelli [41] propose a different approach to carry out bounds tightening called nonlinearities removal domain reduction (NRDR). NRDR relies on the solution of parametric univariate optimization problems to find tight bounds. The method reduces nonlinearities due to a single variable at every iteration. The authors further demonstrate that, under certain assumptions, their domain reduction method is equivalent to carrying out optimization-based bounds tightening iteratively until it reaches a fixed point. Caprara et al. [42] show that the NRDR method leads to bounds consistency for a special case of linear multiplicative programming problems with only two variables in the objective function.

A unified theory for feasibility- and optimality-based tightening

Tawarmalani and Sahinidis [185] provide a unified theory of feasibility- and optimality-based bounds tightening techniques. This theory evolves around Lagrangian subproblems created by the dualization of constraints. Consider the problem:

[TABLE]

Define the Lagrangian subproblem as follows:

[TABLE]

The dual variables $y_{0}$ and $y$ are nonpositive. Suppose an upper bound $U$ is available for the optimal solution of Problem 28. A domain reduction master problem can be constructed as follows:

[TABLE]

Using the linear objective function $h(\mu_{0},\vec{\mu})=a_{0}\mu_{0}+\vec{a}^{T}\vec{\mu}$, Problem 29 can be equivalently stated as:

[TABLE]

Solving Problem 30 can be as hard as solving Problem 28. Instead, the lower bound for Problem 30 can be obtained from the solution of the following relaxed master problem:

[TABLE]

Many domain reduction operations can be obtained by suitable choice of the coefficients $(a_{0},\vec{a})$ in the objective function for Problem 31. For example, equation 25 can be derived by setting $(\vec{a},a_{0})$ to $\vec{e}_{j}$ in Problem 31, where $\vec{e}_{j}$ is the $j^{\rm th}$ column of the identity matrix.

Similarly, propagation for linear constraints described by equations 15 and 16 can derived by application of duality theory to the relaxed problem formed by relaxing all but the $i^{\rm th}$ linear constraint under consideration:

[TABLE]

Assuming $a_{ih}<0$ and $\vec{x}_{h}$ is not at its lower bound, the optimal dual solution for the linear program can be obtained as:

[TABLE]

Here, $\mu$ is the dual multiplier corresponding to $\vec{a}_{i}^{T}\vec{x}\leq b_{i}$, $\vec{\lambda}_{j}$ is the dual multiplier corresponding to $\vec{x}_{j}\leq\vec{x}_{j}^{u}$, and $\vec{\sigma}_{j}$ is the dual multiplier corresponding to $\vec{x}_{j}\geq\vec{x}_{j}^{l}$. The domain reduction master problem can be constructed as:

[TABLE]

Equation 16 follows from Problem 32. Equation 15 can be similarly derived.

Tawarmalani and Sahinidis [185] also present a duality-based reduction scheme that utilizes dual feasible solutions and a learning reduction heuristic. In branch-and-bound algorithms, branching is typically carried out by various heuristics for the selection of the branching variable and the branching point. If one of the nodes created due to a branching decision is proven to be inferior, the learning reduction heuristic attempts to expand on the region defined by this node that is proven to be inferior by the construction of dual solutions for the other node. The authors remark that all dual feasible solutions can be utilized to carry out reductions. This idea is instantiated in the Lagrangian variable bounds propagation by Gleixner and coworkers [71, 72]. To avoid solving $2n$ optimization problems at every node with Problem 27, valid Lagrangian variable inequalities can be generated by aggregating the linear relaxation constraints with the dual solution of Problem 27. While redundant for the linear relaxation, these inequalities approximate the effect of solving Problem 27 locally at every node and can be used to infer stronger bounds on variables when variable bounds are updated through branching or otherwise if a better upper bound is obtained. Gleixner et al. [71] provides heuristics for ordering the generated Lagrangian variable inequalities to achieve maximum tightening. An aggressive filtering strategy is proposed that involves the solution of LPs to determine variables for which bounds cannot be tightened. The optimization steps for these variables in Problem 27 can therefore be skipped, leading to computational savings. The authors also provide heuristics for ordering the LP solves in Problem 27 to allow for more efficient warm starting and reduce the number of simplex iterations.

Pruning

If the lower bound obtained at a node of the branch-and-bound search is worse than the best known upper bound $U$, the current node can be fathomed. We are guaranteed that the globally optimal solution cannot lie at this node, since all feasible solutions at this node have an objective value greater than $U$. This process is also referred to as pruning and is a special case of a general concept called dominance; node $n_{1}$ dominates node $n_{2}$ if, for every feasible solution $s$ in $n_{2}$, there is a complete solution in $n_{1}$ that is as good or better than $s$. The idea of dominance is old, first introduced by Kohler and Steiglitz [105] and developed in more detail by Ibaraki [92]. Dominance relations can be utilized to speed up the search. For MIP problems, Fischetti and Toth [65] propose solution of an auxiliary MIP involving only fixed variables at a node to determine whether the current node is dominated by another node which may or may not have been explored yet. However, the overhead of solving the auxiliary MIP is fairly large. Fischetti and Salvagnin [64] propose improvements to this scheme and demonstrate computational benefits for network loading problems arising in telecommunication. Sewell et al. [172] propose a memory-based dominance rule for a scheduling application. Memory-based dominance rules require storage of the entire search tree, and their performance is dependent on the memory available. However, since information from the entire tree is available, considerably stronger pruning rules can be determined. Memory-based search algorithms are related to the heuristic search algorithms from AI [57, 182].

Problems involving integer variables often have a high degree of symmetry. Symmetry is detrimental for branch-and-bound algorithms as it can lead to repetitive work for the solver. Most symmetry breaking approaches rely on exploiting information about the specific problem being considered. A more general technique called isomorphic pruning has been proposed by Margot [123, 124]. Isomorphic pruning relies on lexicographic tests to determine if a node can be pruned. Orbital branching [146] is another method that can be utilized to tackle symmetry. The branching method identifies orbits of equivalent variables. Orbital branching proceeds by fixing a variable in the orbit to one at a node, and fixes all the variables in the node to zero in another node. Thus, orbital branching implicitly prunes all the other nodes which involve each of the other variables in the orbit fixed to one.

Exploiting optimality conditions

The Karush-Kuhn-Tucker [95, 106] conditions are necessary for a point to be locally optimal for a nonlinear programming problem under certain constraint qualifications. See Schichl and Neumaier [168] for a derivation and a general discussion of these conditions. The conditions can be used to reduce the search space for a nonconvex NLP, since globally optimal points also satisfy them. Vandenbussche and Nemhauser generate valid inequalities for quadratic programs with box constraints through the analysis of optimality conditions [196] and also utilize them in a branch-and-cut scheme [195]. Optimality conditions have been extensively utilized in branch-and-bound algorithms for quadratic programs [83, 38, 39, 37, 45, 90]. Optimality conditions can also be utilized for pruning of nodes. For example, for an unconstrained optimization problem, if 0 is not contained within the interval inclusion function for the partial derivatives of the objective function at a node, the corresponding node can be pruned [125]. This test is often referred to as the monotonicity test.

Sahinidis and Tawarmalani [161] have added a modelling language construct for BARON which allows for the specification of certain constraints as relaxation-only. Relaxation-only constraints are utilized for the construction of convex relaxations and for domain reduction, but are not utilized in local search for obtaining upper bounds. The authors use first-order optimality conditions explicitly as relaxation-only constraints and observe improved convergence for some univariate optimization problems. Amaran and Sahinidis [9] use a similar strategy and show significant computational benefits for parameter estimation problems. They analyze their results and demonstrate that explicit use of optimality conditions aids in domain reduction steps of BARON leading to computational speedups. Puranik and Sahinidis [151] propose a strategy for carrying out implicit bounds tightening on optimality conditions for bound-constrained optimization problems. Their method does not require the generation of optimality conditions which can be time consuming and lead to increase in memory requirements. For a large collection of test problems, this strategy leads to computational speedups and reduction in the number of nodes.

Cluster effect

Branch-and-bounds methods often undergo repeated branching in the neighbourhood of the global solution before converging. This problem is referred to as the cluster effect and was first studied by Du and Kearfott [56] in the context of interval-based branch-and-bound methods. Subsequent analysis [144, 169, 200] also demonstrates the importance of convergence order of the bounding operation (see Definition 1 in [200]). These results indicate that at least second-order convergence is required to overcome the cluster effect. Thus, tighter relaxations can indeed help mitigate the cluster effect and their development is the subject of extensive research in the area [174, 175, 155, 176, 184, 134, 20, 102, 103, 183, 202, 101, 19]. Neumaier and coworkers provide methods to construct exclusion regions for the solution of systems of equations [166] and for the solution of optimization problems [165]. These exclusion regions guarantee that no other solution can lie within the exclusion box around a local minimizer or a solution for systems of equations. These boxes can then be eliminated from the search space. These boxes are constructed through existence and uniqueness tests based on Krawczyk operator or the Kantorovich Theorem (see Chapter 1, [98]). Other methods to construct exclusion regions include back boxing [194] and $\epsilon$-inflation [130].

Implementation of domain reduction techniques

Domain reduction strategies, if successful, typically lead to a reduction in the number of nodes searched in a branch-and-bound tree. Techniques like propagation are computationally inexpensive and can be applied at every node without much overhead. However, they are not as efficient in carrying out domain reduction as the more computationally intensive strategies like Problem 11 or probing. Thus, there exists a need to balance the effort involved in domain reduction in order to reduce the average effort per node. Multiple heuristic or learning strategies are employed for this purpose. These ideas are important in constraint satisfaction problems where multiple filtering algorithms are available achieving varying degrees of consistency. Stergiou [181] experimentally analysed the domain reduction events through filtering techniques. Drawing insights from the clustering of this experimental data, various heuristics are proposed for choosing the propagation algorithm on the fly. Based on insights from Stergiou [181] that propagation events often occur in close clusters, Araya et al. [16] propose an adaptive strategy. Thus, if a constraint propagation mechanism succeeds in carrying out domain reduction at a given node, it should be exploited repetitively in other nodes that are geometrically close to the node until the method fails.

Similar heuristics are also utilized for solution of MINLP problems. For example, the aggressive probing strategy of Nannicini et al. [140] solves probing problems to global optimality and thus has a huge computational overhead. To avoid excessive work, Nannicini et al. propose a strategy based on support vector machines [170] to predict when this aggressive probing strategy is likely to succeed based on the success of propagation operations. Aggressive probing is only carried out when its chances of success are high. Couenne solves linear versions of Problem 11 in the branch-and-bound tree in all nodes up to a depth specified by a parameter $L$. For nodes at a depth $d>L$, the strategy is applied with a probability $2^{L-d}$. A similar strategy is employed by ANTIGONE. The idea of propagation of Lagrangian variable bounds [72, 71] can also be thought of as a means of balancing the computational effort for tightening based on solving variants of Problem 11. Vu et al. [198] show significant computational benefits by carrying out propagation on a subset of the nodes and a partial subgraph of the DAG rather than the entire graph. Vu et al. [197] present multiple strategies for combining the various reduction schemes from constraint programming and mathematical programming for constraint satisfaction problems.

Computational impact of domain reduction

We demonstrate the computational impact of domain reduction techniques on three widely available solvers: BARON [185], Couenne [25] and SCIP [3]. BARON is commercial [21] and also free through the NEOS server [52]. SCIP is free for academics and commercial for all others. Couenne is open-source and free software. All three solvers are available under the GAMS modeling system [36] and provide options that allow us to turn off their domain reduction algorithms. The global MINLP solvers Antigone [135] and LindoGlobal [117] are also available under GAMS but were not used in these experiments since they do not offer facilities that turn off their presolve routines.

The test libraries used in our tests are the Global library [73], Princeton library [149], MINLP library [40] and the CMU-IBM library [48]. Global and Princeton libraries consist of NLP problems, while the MINLP and CMU-IBM libraries consist of MINLP problems. While over 25% of the Princeton library are convex and the NLP relaxations of the CMU-IBM library problems are all convex, the Global and MINLP libraries contain mostly nonconvex problems. In the sequel we present results for each library separately so as to demonstrate that the observed trends are not dominated by any particular library, convexity, or integrality properties. Key statistics on the test sets are summarized in Table 1.

All computational tests were run on a 64-bit Intel Xeon X5650 2.66 Ghz processor running CentOS release 7. The tests were carried out with a time limit of 500 seconds and absolute and relative optimality tolerances set to $10^{-6}$. All solvers were run under two different settings: (1) default options and (2) default options with domain reduction turned off. Our main objective is to compare these two settings. Comparing the relative impact of the individual reduction techniques on the performance of global solvers is beyond the scope of this work. Readers can refer to Kılınç and Sahinidis [104] for experiments showing the effect of different bounds reduction strategies in BARON. In the results presented below, the suffix “nr” is used to denote a solver applied with domain reduction techniques turned off. Comparisons are presented in terms of performance profiles generated through PAVER [136]. In these profiles, a solver is considered to have solved a problem if it obtains the best solution amongst the solvers compared in the profile within a given multiple of the time taken by the fastest solver that solves a problem. CAN_SOLVE denotes the number of problems in the library that can be solved with all the solvers compared in the profile.

BARON

The solver options and their values utilized for tests with BARON are summarized in Table 2. These turn off the various domain reduction techniques utilized in BARON. The remaining options are utilized at their default settings.

Figures 5, 6, 7 and 8 demonstrate the performance of BARON with and without domain reduction techniques employed on the Global, Princeton, CMU-IBM and MINLP libraries. The profiles indicate a huge deterioration in performance across all test libraries when reduction is turned off. Interestingly, for the continuous test libraries, the performance profile of the no-reduction version of BARON “catches up” with that of the reduction-based version at the end of the profiles. This simply means that, without reduction, BARON’s heuristics are still able to come up with (near-)global solutions but the branch-and-bound algorithm is not able to provide sufficiently strong lower bounds in order to prove global optimality. For the MINLP case, the no-reduction-based algorithm is not even able to find good feasible solutions.

Results in Table 3 show that turning off domain reduction techniques leads to a huge increase in the number of nodes required by BARON. The increase in nodes is also accompanied by a significant increase in computational time across all test libraries.

Couenne

The solver options and their values utilized for tests with Couenne are summarized in Table 4. These turn off the various domain reduction techniques utilized in Couenne. The remaining options are utilized at their default settings.

Figures 9, 10, 11 and 12 demonstrate the performance of Couenne with and without domain reduction techniques employed on the Global, Princeton, CMU-IBM and MINLP libraries. Similar to BARON, turning off domain reduction techniques has a huge impact on the performance of Couenne. Contrary to BARON, without reduction, this solver is unable to find good solutions for the NLP test libraries; as a result, the performance profiles for the no-reduction-based algorithm do not catch up with the reduction-based algorithm. Table 5 indicates a significant increase in computational time and the number of nodes explored in the branch-and-bound search when reduction is turned off. It should be pointed out that no comparisons are possible between BARON and Couenne by looking at their respective performance profiles since these profiles depend solely on the solvers included in each figure.

SCIP

The solver options and their values utilized for tests with SCIP are summarized in Table 6. As before, the remaining options are used at their default settings.

Figures 13, 14, 15 and 16 demonstrate the huge impact of turning off domain reductions with SCIP on the Global, Princeton, CMU-IBM and MINLP libraries. Similar to Couenne, this solver is also not able of finding good solutions for continuous problems when reduction is turned off. Results in Table 7 indicate a significant deterioration in performance for SCIP as other solvers.

Relative solver performance

Figure 17 describes the performance of BARON, Couenne and SCIP on all test libraries aggregated together. Even without reduction, BARON dominates over the other two solvers, suggesting that this solver has an edge over the other two solvers in terms of its technology for relaxation construction, branching schemes, and primal feasibility heuristics. However, the impact of domain reduction techniques is clear and rather substantial on all three solvers. Interestingly, the relative order of SCIP and Couenne changes when reduction is turned off, suggesting that SCIP is relying much more on domain reduction than Couenne does. Equivalently, SCIP enjoys a substantial advantage over Couenne thanks to its implementation of a more extensive set of reduction techniques.

Conclusions

We have presented a review of the various domain reduction techniques proposed in literature for the purpose of global NLP and MINLP optimization. These techniques vary in complexity including simple ones like propagation to more computationally intensive ones that involve full solution of optimization subproblems. Application of some of the more complex techniques requires the use of smart heuristics to ensure that they are utilized only when domain reduction is likely. We have also presented computational results with BARON, SCIP and Couenne on publicly available test libraries. The results show that domain reduction techniques have a significant impact on the performance of these solvers. Incorporation of domain reduction within branch-and-bound leads to huge reductions in computational time and number of nodes required for solution. Future research in this area should focus on the development of domain reduction techniques based on sets of constraints (i.e., global filtering methods) for broad classes of structured problems but also for important science and engineering problems, such as pooling problems [153], protein folding [143], and network design problems [53, 156].