Leo Liberti
email: liberti@lix.polytechnique.fr

Sonia Cafieri
email: cafieri@lix.polytechnique.fr

Fabien Tarissan
email: tarissan@lix.polytechnique.fr

Reformulations in Mathematical Programming: A Computational Approach

Mathematical programming is a language for describing optimization problems; it is based on parameters, decision variables, objective function(s) subject to various types of constraints. The present treatment is concerned with the case when objective(s) and constraints are algebraic mathematical expressions of the parameters and decision variables, and therefore excludes optimization of black-box functions. A reformulation of a mathematical program P is a mathematical program Q obtained from P via symbolic transformations applied to the sets of variables, objectives and constraints. We present a survey of existing reformulations interpreted along these lines, some example applications, and describe the implementation of a software framework for reformulation and optimization.

Introduction

Optimization and decision problems are usually defined by their input and a mathematical description of the required output: a mathematical entity with an associated value, or whether a given entity has a specified mathematical property or not. Mathematical programming is a language designed to express almost all practically interesting optimization and decision problems.

Mathematical programming formulations can be categorized according to various properties, and rather efficient solution algorithms exist for many of the categories. As in most languages, the same semantics can be conveyed by many different syntactical expressions. In other words, there are many equivalent formulations for each given problem (what the term "equivalent" means in this context will be defined later). Furthermore, solution algorithms for mathematical programming formulations often rely on solving a sequence of different problems (often termed auxiliary problems) related to the original one: although these are usually not fully equivalent to the original problem, they may be relaxations, projections, liftings, decompositions (among others). Auxiliary problems are reformulations of the original problem.

Consider for example the Kissing Number Problem (KNP) in D dimensions [START_REF] Kucherenko | New formulations for the kissing number problem[END_REF], i.e. the determination of the maximum number of unit D-dimensional spheres that can be arranged around a central unit D-dimensional sphere. As all optimization problems, this can be cast (by using a bisection argument) as a sequence of decision problems on the cardinality of the current spheres configuration. Namely, given the positive integers D (dimension of Euclidean space) and N , is there a configuration of N unit spheres around the central one? For any fixed D, the answer will be affirmative or negative depending on the value of N . The highest N such that the answer is affirmative is the kissing number for dimension D. The decision version of the KNP can be cast as a nonconvex Nonlinear Programming (NLP) feasibility problem as follows. For all i ≤ N , let x i = (x i1 , . . . , x iD ) ∈ R D be the center of the i-th sphere. We look for a set of vectors {x i | i ≤ N } satisfying the following constraints:

∀ i ≤ N ||x i || = 2 ∀ i < j ≤ N ||x i -x j || ≥ 2 ∀ i ≤ N -2 ≤ x i ≤ 2.
It turns out that this problem is numerically quite difficult to solve, as it is very unlikely that the local NLP solution algorithm will be able to compute a feasible starting solution. Failing to find an initial feasible solution means that the solver will immediately abort without having made any progress. Most researchers with some practical experience in NLP solvers (such as e.g. SNOPT [START_REF] Gill | User's Guide for SNOPT 5.3[END_REF]), however, will immediately reformulate this problem into a more computationally amenable form by squaring the norms to get rid of a potentially problematic square root, and treating the reverse convex constraints ||x ix j || ≥ 2 as soft constraints by multiplying the right hand sides by a non-negative scaling variable α, which is then maximized:

max α (1) ∀ i ≤ N ||x i || 2 = 4 (2) ∀ i < j ≤ N ||x i -x j || 2 ≥ 4α. ( 3
) ∀ i ≤ N -2 ≤ x i ≤ 2 (4) α ≥ 0. (5)
In this form, finding an initial feasible solution is trivial; for example, x i = (2, 0, . . . , 0) for all i ≤ N will do. Subsequent solver iterations will likely be able to provide a better solution. Should the computed value of α be ≥ 1, the solution would be feasible in the hard constraints, too. Currently, we are aware of no optimization language environment that is able to perform the described reformulation automatically. Whilst this is not a huge limitation for NLP experts, people who simply wish to model a problem and get its solution will fail to obtain one, and may even be led into thinking that the formulation itself is infeasible.

Another insightful example of the types of limitations we refer to can be drawn from the KNP. We might wish to impose ordering constraints on some of the spheres to reduce the number of symmetric solutions. Ordering spheres packed on a spherical surface is hard to do in Euclidean coordinates, but it can be done rather easily in spherical coordinates, by simply stating that the value of a spherical coordinate of the i-th sphere must be smaller than the corresponding value in the j-th sphere. We can transform a Euclidean coordinate vector x = (x 1 , . . . , x D ) in D-spherical coordinates (ρ, ϑ 1 , . . . , ϑ D-1 ) such that ρ = ||x|| and ϑ ∈ [0, 2π] D-1 by means of the following equations:

ρ = ||x|| (6) ∀k ≤ D x k = ρ sin ϑ k-1 D-1 h=k cos ϑ h (7)
(this yields another NLP formulation of the KNP). Applying the D-spherical transformation is simply a matter of term rewriting and algebraic simplification, and yet no currently existing optimization language environment offers such capabilities. By pushing things further, we might wish to devise an algorithm that dynamically inserts or removes constraints expressed in either Euclidean or spherical coordinates depending on the status of the current solution, and re-solves the (automatically) reformulated problem at each iteration. This may currently be done (up to a point) by optimization language environments such as AMPL [START_REF] Fourer | The AMPL Book[END_REF], provided all constraints are part of a pre-specified family of parametric constraints. Creating new constraints by term rewriting, however, is not a task currently addressed by current mathematical programming implementations. The limitations emphasized in the KNP example illustrate a practical need for very sophisticated software including numerical as well as symbolic algorithms, both applied to the unique goal of solving optimization problems cast as mathematical programming formulations. The current state of affairs is that there are many numerical optimization solvers and many Computer Algebra Systems (CAS) -such as Maple or Mathematica -whose efficiency is severely hampered by the full generality of their capabilities. In short, we would ideally need (small) parts of the symbolic kernels driving the existing CASes to be combined with the existing optimization algorithms, plus a number of super-algorithms capable of making automated, dynamic decisions on the type of reformulations that are needed to improve the current search process.

Although the above paradigm might seem far-fetched, it does in fact already exist in the form of the hugely successful CPLEX [START_REF]ILOG: ILOG CPLEX 11.0 User's Manual[END_REF] solver targeted at solving Mixed-Integer Linear Programming (MILP) problems. The initial formulation provided by the user is automatically simplified and improved with a sizable variety of pre-processing steps which attempt to reduce the number of variables and constraints. Thereafter, at each node of the Branch-and-Bound algorithm, the formulation may be tightened as needed by inserting and removing additional valid constraints, in the hope that the current relaxed solution of the (automatically obtained) linear relaxation is improved. Advanced users may of course decide to tune many parameters controlling this process, but practitioners needing a practical answer can simply use default parameters and to let CPLEX decide what is best. Naturally, the task carried out by CPLEX is greatly simplified by the assumption that both objective function and constraints are linear forms, which is obviously not the case in a general nonlinear setting.

In this chapter we attempt to move some steps in the direction of endowing general mathematical programming with the same degree of algorithmic automation enjoyed by linear programming. We propose: (a) a theoretical framework in which mathematical programming reformulations can be formalized in a unified way, and (b) a literature review of the most successful existing reformulation and relaxation techniques in mathematical programming. Since an all-comprehensive literature review in reformulation techniques would extend this chapter to possibly several hundreds (thousands?) pages, only a partial review has been provided. In this sense, this should be seen as "work in progress" towards laying the foundations to a computer software which is capable of reformulating mathematical programming formulations automatically. Note also that for this reason, the usual mathematical notations have been translated to a data structure framework that is designed to facilitate computer implementation. Most importantly, "functions" -which as mathematical entities are interpreted as maps between sets -are represented by expression trees: what is meant by the expression x + y, for example, is really a directed binary tree on the vertices {+, x, y} with arcs {(+, x), (+, y)}. For clarity purposes, however, we also provide the usual mathematical languages.

One last (but not least) remark is that reformulations can be seen as a new way of expressing a known problem. Reformulations are syntactical operations that may add or remove variables or constraints, whilst keeping the fundamental structure of the problem optima invariant. When some new variables are added and some of the old ones are removed, we can usually try to re-interpret the reformulated problem and assign a meaning to the new variables, thus gaining new insights to the problem. One example of this is given in Sect. 3.5. One other area in mathematical programming that provides a similarly clear relationship between mathematical syntax and semantics is LP duality with the interpretation of reduced costs. This is important insofar as it offers alternative interpretations to known problems, which gains new and useful insights.

The rest of this chapter is organized as follows. In Section 2 we propose a general theoretical framework of definitions allowing a unified formalization of mathematical programming reformulations. The definitions allow a consistent treatment of the most common variable and constraint manipulations in mathematical programming formulations. In Section 3 we present a systematic study of a set of well known reformulations. Most reformulations are listed as symbolic algorithms acting on the problem structure, although the equivalent transformation in mathematical terms is given for clarity purposes. In Section 4 we present a systematic study of a set of well known relaxations. Again, relaxations are listed as symbolic algorithms acting on the problem structure whenever possible, the equivalent mathematical transformation being given for clarity. Section 5 describes the implementation of ROSE, a Reformulation/Optimization Software Engine.

General framework

In Sect. 2.1 we formally define what a mathematical programming formulation is. In Sect. 2.2 we discuss the expression tree function representation. Sect. 2.3 lists the most common standard forms in mathematical programming.

A data structure for mathematical programming formulations

In this chapter we give a formal definition of a mathematical programming formulation in such terms that can be easily implemented on a computer. We then give several examples to illustrate the generality of our definition. We refer to a mathematical programming problem in the most general form:

min f (x) g(x) ⋚ b x ∈ X,    (8)
where f, g are function sequences of various sizes, b is an appropriately-sized real vector, and X is a cartesian product of continuous and discrete intervals.

The precise definition of a mathematical programming formulation lists the different formulation elements: parameters, variables having types and bounds, expressions depending on the parameters and variables, objective functions and constraints depending on the expressions. We let P be the set of all mathematical programming formulations, and M be the set of all matrices. This is used in Defn. 1 to define leaf nodes in mathematical expression trees, so that the concept of a formulation can also accommodate multilevel and semidefinite programming problems. Notationwise, in a digraph (V, A) for all v ∈ V we indicate by δ + (v) the set of vertices u for which (v, u) ∈ A and by δ -(v) the set of vertices u for which (u, v) ∈ A.

Definition 1. Given an alphabet L consisting of countably many alphanumeric names N L and operator symbols O L , a mathematical programming formulation P is a 7-tuple (P, V, E, O, C, B, T ), where:

• P ⊆ N L is the sequence of parameter symbols: each element p ∈ P is a parameter name; • V ⊆ N L is the sequence of variable symbols: each element v ∈ V is a variable name; • E is the set of expressions: each element e ∈ E is a Directed Acyclic Graph (DAG) e = (V e , A e ) such that: (a) V e ⊆ L is a finite set (b) there is a unique vertex r e ∈ V e such that δ -(r e ) = ∅ (such a vertex is called the root vertex) (c) vertices v ∈ V e such that δ + (v) = ∅ are called leaf vertices and their set is denoted by λ(e); all leaf vertices v are such that v ∈ P ∪ V ∪ R ∪ P ∪ M (d) for all v ∈ V e such that δ + (v) = ∅, v ∈ O L (e) two weight functions χ, ξ : V e → R are defined on V e : χ(v) is the node coefficient and ξ(v)

is the node exponent of the node v; for any vertex v ∈ V e , we let τ (v) be the symbolic term of v: namely, v = χ(v)τ (v) ξ (v) . elements of E are sometimes called expression trees; nodes v ∈ O L represent an operation on the nodes in δ + (v), denoted by v(δ + (v)), with output in R; |V| is the sequence of variable types: for all v ∈ V, v is called a continuous variable if T (v) = 0, an integer variable if T (v) = 1 and a binary variable if T (v) = 2.

• O ⊆ {-1, 1} × E
c ≡    e c ≤ b c if s c = -1 e c = b c if s c = 0 e c ≥ b c if s c = 1; • B ⊆ R |V| × R |V| is the sequence of variable bounds: for all v ∈ V let B(v) = [L v , U v ] with L v , U v ∈ R; • T ⊆ {0, 1, 2}
We remark that for a sequence of variables z ⊆ V we write T (z) and respectively B(z) to mean the corresponding sequences of types and respectively bound intervals of the variables in z. Given a formulation P = (P, V, E, O, C, B, T ), the cardinality of P is |P | = |V|. We sometimes refer to a formulation by calling it an optimization problem or simply a problem.

Definition 2. Any formulation Q that can be obtained by P by a finite sequence of symbolic operations carried out on the data structure is called a problem transformation.

Examples

In this section we provide some explicitly worked out examples that illustrate Defn. 1.

A quadratic problem

Consider the problem of minimizing the quadratic form 3x 2 1 + 2x 2 2 + 2x 2 3 + 3x 2 4 + 2x 2 5 + 2x 2 6 -2x 1 x 2 -2x 1 x 3 -2x 1 x 4 -2x 2 x 3 -2x 4 x 5 -2x 4 x 6 -2x 5 x 6 subject to x 1 + x 2 + x 3 + x 4 + x 5 + x 6 = 0 and x i ∈ {-1, 1} for all i ≤ 6. For this problem,

• P = ∅; • V = (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ); • E = (e 1 , e 2 )
where e 1 , e 2 are the graphs shown in Fig. 1;

• O = (-1, e 1 ); • C = ((e 2 , 0, 0)); • B = ([-1, 1], [-1, 1], [-1, 1], [-1, 1], [-1, 1], [-1, 1]); • T = (2, 2, 2, 2, 2, 2). ^^^^^+ × × × × × × × × × × × × × × × × × × × × 3 3 2 2 2 2 2 2 2 2 2 2 -2 -2 -2 -2 -2 -2 -2 x1 x1 x1 x1 x2 x2

Balanced graph bisection

Example 2.1 is a (scaled) mathematical programming formulation of a balanced graph bisection problem instance. This problem is defined as follows.

Balanced Graph Bisection Problem (BGBP). Given an undirected graph G = (V, E) without loops or parallel edges such that |V | is even, find a subset

U ⊂ V such that |U | = |V | 2 and the set of edges C = {{u, v} ∈ E | u ∈ U, v ∈ U } is as small as possible.
The problem instance considered in Example 2.1 is shown in Fig. 2. To all vertices i ∈ V we associate variables xi =

 1 i ∈ U -1 i ∈ U . The number of edges in C is counted by 1 4 {i,j}∈E (x i -x j ) 2 .
The fact that |U | = |V | 2 is expressed by requiring an equal number of variables at 1 and -1, i.e. We can also express the problem in Example 2.1 as a particular case of the more general optimization problem:

min x x ⊤ Lx s.t. x1 = 0 x ∈ {-1, 1} 6 ,    where L =         3 -1 -1 -1 0 0 -1 2 -1 0 0 0 -1 -1 2 0 0 0 -1 0 0 3 -1 -1 0 0 0 -1 2 -1 0 0 0 -1 -1 2         and 1 = (1, 1, 1, 1, 1, 1)
⊤ . We represent this class of problems by the following mathematical programming formulation:

• P = (L ij | 1 ≤ i, j ≤ 6); • V = (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ); • E = (e ′
1 , e 2 ) where e ′ 1 is shown in Fig. 3 and e 2 is shown in Fig. 1 (below);

• O = (-1, e ′ 1 ); • C = ((e 2 , 0, 0)); • B = ([-1, 1], [-1, 1], [-1, 1], [-1, 1], [-1, 1], [-1, 1]); • T = (2, 2, 2, 2, 2, 2).

The Kissing Number Problem

The kissing number problem formulation (1)-( 5) is as follows:

• P = (N, D); • V = (x ik | 1 ≤ i ≤ N ∧ 1 ≤ k ≤ D) ∪ {α}; • E = (f, h j , g ij | 1 ≤ i < j ≤ N )
, where f is the expression tree for α, h j is the expression tree for ||x j || 2 for all j ≤ N , and g ij is the expression tree for ||x ix j || 2 -4α for all i < j ≤ N ;

• O = (1, f ); • C = ((h i , 0, 4) | i ≤ N ) ∪ ((g ij , 1, 0) | i < j ≤ N ); • B = [-2, 2] N D ; • T = {0} N D .
As mentioned in Section 1, the kissing number problem is defined as follows.

2 2 2 2 2 2 replacemen + × × × × × × × × × × × × × × × × × × × × L 11 L 22 L 33 L 44 L 55 L 66 L ′ 12 L ′ 13 L ′ 14 L ′ 23 L ′ 45 L ′ 46 L ′ 56
L ′ ij = Lij + Lji for all i, j.
Kissing Number Problem (KNP). Find the largest number N of non-overlapping unit spheres in R D that are adjacent to a given unit sphere.

The formulation of Example 2.1 refers to the decision version of the problem: given integers N and D, is there an arrangement of N non-overlapping unit spheres in R D adjacent to a given unit sphere?

A data structure for mathematical expressions

Given an expression tree DAG e = (V, A) with root node r(e) and whose leaf nodes are elements of R or of M (the set of all matrices), the evaluation of e is the (numerical) output of the operation represented by the operator in node r applied to all the subnodes of r (i.e. the nodes adjacent to r); in symbols, we denote the output of this operation by r(δ + (r)), where the symbol r denotes both a function and a node. Naturally, the arguments of the operator must be consistent with the operator meaning. We remark that for leaf nodes belonging to P (the set of all formulations), the evaluation is not defined; the problem in the leaf node must first be solved and a relevant optimal value (e.g. an optimal variable value, as is the case with multilevel programming problems) must replace the leaf node.

For any e ∈ E, the evaluation tree of e is a DAG ē = ( V , A) where

V = {v ∈ V | |δ + (v)| > 0 ∨ v ∈ R} ∪ {x(v) | |δ + (v)| = 0 ∧ v ∈ V} (in
short, the same as V with every variable leaf node replaced by the corresponding value x(v)). Evaluation trees are evaluated by Alg. 1. We can now naturally extend the definition of evaluation of e at a point x to expression trees whose leaf nodes are either in V or R. Definition 3. Given an expression e ∈ E with root node r and a point x, the evaluation e(x) of e at x is the evaluation r(δ + (r)) of the evaluation tree ē.

We consider a sufficiently rich operator set O L including at least +, ×, power, exponential, logarithm, and trigonometric functions (for real arguments) and inner product (for matrix arguments). Note that since any term t is weighted by a multiplier coefficient χ(t) there is no need to employ aoperator, for it suffices to multiply χ(t) by -1 = ξ(v) in the appropriate term(s) t; a division u/v is expressed by multiplying u by v raised to the power -1. Depending on the problem Algorithm 1 The evaluation algorithm for expression trees.

double Eval(node v) { double ρ; if (v ∈ OL) { // v is an operator array α[|δ + (v)|]; ∀ u ∈ δ + (v) { α(u) =Eval(u); } ρ = χ(v)v(α) ξ(v) ; } else { // v is a constant value ρ = χ(v)v ξ(v) ; } return ρ; }
form, it may sometimes be useful to enrich O L with other (more complex) terms. In general, we view an operator in O L as an atomic operation on a set of variables with cardinality at least 1.

A standard form for expressions

Since in general there is more than one way to write a mathematical expression, it is useful to adopt a standard form; whilst this does not resolve all ambiguities, it nonetheless facilitates the task of writing symbolic computation algorithms acting on the expression trees. For any expression node t in an expression tree e = (V, A):

• if t is a sum: 1. |δ + (t)| ≥ 2 2.
no subnode of t may be a sum (sum associativity); 3. no pair of subnodes u, v ∈ δ + (t) must be such that τ (u) = τ (v) (i.e. like terms must be collected); as a consequence, each sum only has one monomial term for each monomial type 4. a natural (partial) order is defined on δ + (t): for u, v ∈ δ + (t), if u, v are monomials, u, v are ordered by degree and lexicographically

• if t is a product: 1. |δ + (t)| ≥ 2 2.
no subnode of t may be a product (product associativity); 3. no pair of subnodes u, v ∈ δ + (t) must be such that τ (u) = τ (v) (i.e. like terms must be collected and expressed as a power)

• if t is a power: 1. |δ + (t)| = 2
2. the exponent may not be a constant (constant exponents are expressed by setting the exponent coefficient ξ(t) of a term t) 3. the natural order on δ + (t) lists the base first and the exponent later.

The usual mathematical nomenclature (linear forms, polynomials, and so on) applies to expression trees.

Standard forms in mathematical programming

Solution algorithms for mathematical programming problems read a formulation as input and attempt to compute an optimal feasible solution as output. Naturally, algorithms which exploit problem structure are usually more efficient than those that do not. In order to be able to exploit the structure of the problem, solution algorithms solve problems that are cast in a standard form that emphasizes the useful structure. In this section we review the most common standard forms.

Linear Programming

A mathematical programming problem P is a Linear Programming (LP) problem if (a) |O| = 1 (i.e. the problem only has a single objective function); (b) e is a linear form for all e ∈ E; and (c) T (v) = 0 (i.e. v is a continuous variable) for all v ∈ V.

An LP is in standard form if (a) s c = 0 for all constraints c ∈ C (i.e. all constraints are equality constraints) and (b) B(v) = [0, +∞] for all v ∈ V. LPs are expressed in standard form whenever a solution is computed by means of the simplex method [START_REF] Dantzig | Linear Programming and Extensions[END_REF]. By contrast, if all constraints are inequality constraints, the LP is known to be in canonical form.

Mixed Integer Linear Programming

A mathematical programming problem P is a Mixed Integer Linear Programming (MILP) problem if (a) |O| = 1; and (b) e is a linear form for all e ∈ E.

A MILP is in standard form if s c = 0 for all constraints c ∈ C and if B(v) = [0, +∞] for all v ∈ V. The most common solution algorithms employed for solving MILPs are Branch-and-Bound (BB) type algorithms [START_REF]ILOG: ILOG CPLEX 11.0 User's Manual[END_REF]. These algorithms rely on recursively partitioning the search domain in a tree-like fashion, and evaluating lower and upper bounds at each search tree node to attempt to implicitly exclude some subdomains from consideration. BB algorithms usually employ the simplex method as a sub-algorithm acting on an auxiliary problem, so they enforce the same standard form on MILPs as for LPs. As for LPs, a MILP where all constraints are inequalities is in canonical form.

Nonlinear Programming

A mathematical programming problem P is a Nonlinear Programming (NLP) problem if (a)

|O| = 1 and (b) T (v) = 0 for all v ∈ V.
Many fundamentally different solution algorithms are available for locally solving NLPs, and most of them require different standard forms. One of the most widely used is Sequential Quadratic Programming (SQP) [START_REF] Gill | User's Guide for SNOPT 5.3[END_REF], which requires problem constraints c ∈ C to be expressed in the form

l c ≤ c ≤ u c with l c , u c ∈ R ∪ {-∞,

Mixed Integer Nonlinear Programming

A mathematical programming problem P is a Mixed Integer Nonlinear Programming (MINLP) problem if |O| = 1. The situation as regards MINLP standard forms is generally the same as for NLPs, save that a few more works have appeared in the literature about standard forms for MINLPs [113,114,[START_REF] Pantelides | Mixed integer linear/nonlinear programming interface specification[END_REF][START_REF] Liberti | Writing global optimization software[END_REF]. In particular, the Smith standard form [114] is purposefully constructed so as to make symbolic manipulation algorithms easy to carry out on the formulation. A MINLP is in Smith standard form if: Essentially, the Smith standard form consists of a linear part comprising objective functions and a set of constraints; the rest of the constraints have a special form ⊗(x 1 , . . . , x p )-v = 0 for some p ∈ N, with v, x 1 , . . . , x p ∈ V(P ) and ⊗ a nonlinear operator in O L . By grouping all nonlinearities in a set of equality constraints of the form "variable = operator(variables)" (called defining constraints) the Smith standard form makes it easy to construct auxiliary problems. The Smith standard form can be constructed by recursing on the expression trees of a given MINLP [112] and is an optreformulation. Solution algorithms for solving MINLPs are usually extensions of BB type algorithms [114, [START_REF] Liberti | Writing global optimization software[END_REF][START_REF] Liberti | Reformulation and convex relaxation techniques for global optimization[END_REF]124,[START_REF] Nowak | Relaxation and Decomposition Methods for Mixed Integer Nonlinear Programming[END_REF].

• O = {d o , e o }

Separable problems

A problem P is in separable form if (a) O(P ) = {(d o , e o )}, (b) C(P ) = ∅ and (c) e o is such that: • r(e o ) is the sum operator • for all distinct u, v ∈ δ + (r(e o )), λ(u) ∩ λ(v) = ∅,
where by slight abuse of notation λ(u) is the set of leaf nodes of the subgraph of e o whose root is u. The separable form is a standard form by itself. It is useful because it allows a very easy problem decomposition: for all u ∈ δ + (r(e o )) it suffices to solve the smaller problems

Q u with V(Q u ) = λ(v) ∩ V(P ), O(Q u ) = {(d o , u)} and B(Q u ) = {B(P )(v) | v ∈ V(Q u )}. Then u∈δ + (r(eo))
x(V(Q u )) is a solution for P .

Factorable problems

A problem P is in factorable form [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF]130,111,124]

if: 1. O = {(d o , e o )} 2. r(e o ) ∈ V (consequently, the vertex set of e o is simply {r(e o )}) 3. for all c ∈ C: • s c = 0 • r(e c ) is the sum operator • for all t ∈ δ + (r(e c ))
, either (a) t is a unary operator and δ + (t) ⊆ λ(e c ) (i.e. the only subnode of t is a leaf node) or (b) t is a product operator such that for all v ∈ δ + (t), v is a unary operator with only one leaf subnode.

The factorable form is a standard form by itself. Factorable forms are useful because it is easy to construct many auxiliary problems (including convex relaxations, [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF][START_REF] Adjiman | A global optimization method, αBB, for general twice-differentiable constrained NLPs: I. Theoretical advances[END_REF]111]) from problems cast in this form. In particular, factorable problems can be reformulated to emphasize separability [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF]124,[START_REF] Nowak | Relaxation and Decomposition Methods for Mixed Integer Nonlinear Programming[END_REF].

D.C. problems

The acronym "d.c." stands for "difference of convex". Given a set

Ω ⊆ R n , a function f : Ω → R is a d.c.
• O = {(d o , e o )}; • e o is a d.c. function;
• c is a linear form for all c ∈ C.

D.C. programming problems have two fundamental properties.

The first is that the space of all d.c. functions is dense in the space of all continuous functions. This implies that any continuous optimization problem can be approximated as closely as desired, in the uniform convergence topology, by a d.c. optimization problem [129,[START_REF] Horst | Global Optimization: Deterministic Approaches[END_REF]. The second property is that it is possible to give explicit necessary and sufficient global optimality conditions for certain types of d.c. problems [129,116]. Some formulations of these global optimality conditions [115] also exhibit a very useful algorithmic property: if at a feasible point x the optimality conditions do not hold, then the optimality conditions themselves can be used to construct an improved feasible point x ′ .

Linear Complementarity problems

Linear complementarity problems (LCP) are nonlinear feasibility problems with only one nonlinear constraint. An LCP is defined as follows [START_REF] Di Giacomo | Mathematical programming methods in dynamical nonlinear stochastic supply chain management[END_REF], p. 50:

• O = ∅;
• there is a constraint c ′ = (e, 0, 0) ∈ C such that (a) t = r(e) is a sum operator; (b) for all u ∈ δ + (t), u is a product of two terms v, f such that v ∈ V and (f, 1, 0) ∈ C; • for all c ∈ C {c ′ }, e c is a linear form.

Essentially, an LCP is a feasibility problem of the form:

Ax ≥ b x ≥ 0 x ⊤ (Ax -b) = 0,    where x ∈ R n , A is an m × n matrix and b ∈ R m .
Many types of mathematical programming problems (including MILPs with binary variables [START_REF] Di Giacomo | Mathematical programming methods in dynamical nonlinear stochastic supply chain management[END_REF][START_REF] Judice | Reformulation of mathematical programming problems as linear complementarity problems and investigation of their solution methods[END_REF]) can be recast as LCPs or extensions of LCP problems [START_REF] Judice | Reformulation of mathematical programming problems as linear complementarity problems and investigation of their solution methods[END_REF]. Furthermore, some types of LCPs can be reformulated to LPs [START_REF] Mangasarian | Linear complementarity problems solvable by a single linear program[END_REF] and as separable bilinear programs [START_REF] Mangasarian | The linear complementarity problem as a separable bilinear program[END_REF]. Certain types of LCPs can be solved by an interior point method [START_REF] Kojima | An interior point potential reduction algorithm for the linear complementarity problem[END_REF][START_REF] Di Giacomo | Mathematical programming methods in dynamical nonlinear stochastic supply chain management[END_REF].

Bilevel Programming problems

The bilevel programming (BLP) problem consists of two nested mathematical programming problems named the leader and the follower problem.

A mathematical programming problem P is a bilevel programming problem if there exist two programming problems L, F (the leader and follower problem) and a subset ℓ = ∅ of all leaf nodes of E(L) such that any leaf node v ∈ ℓ has the form (v, F ) where v ∈ V(F ).

The usual mathematical notation is as follows [START_REF] Falk | On bilevel programming, part I: General nonlinear cases[END_REF][START_REF] Audet | Links between linear bilevel and mixed 0-1 programming problems[END_REF]:

min y F (x(y), y) min x f (x, y) s.t. x ∈ X, y ∈ Y,    (9)
where X, Y are arbitrary sets. This type of problem arises in economic applications. The leader knows the cost function of the follower, who may or may not know that of the leader; but the follower knows the optimal strategy selected by the leader (i.e. the optimal values of the decision variables of L) and takes this into account to compute his/her own optimal strategy. BLPs can be reformulated exactly to MILPs with binary variables and vice-versa [START_REF] Audet | Links between linear bilevel and mixed 0-1 programming problems[END_REF], where the reformulation is as in Defn. 6. Furthermore, two typical Branch-and-Bound (BB) algorithms for the considered MILPs and BLPs have the property that the the MILP BB can be "embedded" in the BLP BB (this roughly means that the BB tree of the MILP is a subtree of the BB tree of the BLP); however, the contrary does not hold. This seems to hint at a practical solution difficulty ranking in problems with the same degree of worst-case complexity (both MILPs and BLPs are NP-hard).

Semidefinite Programming problems

Consider known symmetric n×n matrices C, A k for k ≤ m, a vector b ∈ R m and a symmetric n×n matrix X = (x ij ) where x ij is a problem variable for all i, j ≤ n. The following is a semidefinite programming problem (SDP) in primal form:

min X C • X ∀k ≤ m A k • X = b k X 0,    ( 10
)
where X 0 is a constraint that indicates that X should be symmetric positive semidefinite, and

C • X = tr(C ⊤ X) = i,j c ij x ij .
We also consider the SDP in dual form:

max y,S b ⊤ y k≤m y k A k + S = C S 0,    ( 11
)
where S is a symmetric n × n matrix and y ∈ R m . Both forms of the SDP problem are convex NLPs, so the duality gap is zero. Both forms can be solved by a particular type of polynomial-time interior point method (IPM), which means that solving SDPs is practically efficient [START_REF] Alizadeh | Interior point methods in semidefinite programming with applications to combinatorial optimization[END_REF]125]. SDPs are important because they provide tight relaxations to (nonconvex) quadratically constrained quadratic programming problems (QCQP), i.e. problems with a quadratic objective and quadratic constraints (see Sect. 4.3). SDPs can be easily modelled with the data structure described in Defn. 1, for their expression trees are linear forms where each leaf node contains a symmetric matrix. There is no need to explicitly write the semidefinite constraints X 0, S 0 because the solution IPM algorithms will automatically find optimal X, S matrices that are semidefinite.

Reformulations

In this section we define some types of reformulations and establish some links between them (Sect. 3.1) and we give a systematic study of various types of elementary reformulations (Sect. 3.2) and exact linearizations (Sect. 3.3). Sect. 3.5 provides a few worked out examples. In this summary, we tried to focus on two types of reformulations: those that are in the literature, but may not be known to every optimization practitioner, and those that represent the "tricks of the trade" of most optimization researchers but have never, or rarely, been formalized explicitly; so the main contributions of this section are systematic and didactic. Since the final aim of automatic reformulations is let the computer arrive at an alternative formulation which is easier to solve, we concentrated on those reformulations which simplified nonlinear terms into linear terms, or which reduced integer variables to continuous variables. By contrast, we did not cite important reformulations (such as the LP duality) which are fundamental in solution algorithms and alternative problem interpretation, but which do not significantly alter solution difficulty.

Reformulation definitions

Consider a mathematical programming formulation P = (P, V, E, O, C, B, T ) and a function x : V → R |V| (called point) which assigns values to the variables. Definition 4. A point x is type feasible if:

x(v) ∈    R if T (v) = 0 Z if T (v) = 1 {L v , U v } if T (v) = 2 for all v ∈ V; x is bound feasible if x(v) ∈ B(v) for all v ∈ V; x is constraint feasible if for all c ∈ C we have: e c (x) ≤ b c if s c = -1, e c (x) = b c if s c = 0, and e c (x) ≥ b c if s c = 1. A point x is feasible in P if it
is type, bound and constraint feasible.

A point x feasible in P is also called a feasible solution of P . A point which is not feasible is called infeasible. Denote by F (P ) the feasible points of P . Definition 5. A feasible point x is a local optimum of P with respect to the objective o ∈ O if there is a non-empty neighbourhood N of x such that for all feasible points y = x in N we have

d o f o (x) ≥ d o f o (y). A local optimum is strict if d o f o (x) > d o f o (y). A feasible point x is a global optimum of P with respect to the objective o ∈ O if d o f o (x) ≥ d o f o (y) for all feasible points y = x. A global optimum is strict if d o f o (x) > d o f o (y).
Denote the set of local optima of P by L(P ) and the set of global optima of P by G(P ). If O(P ) = ∅, we define L(P ) = G(P ) = F (P ).

Example 1. The point x = (-1, -1, -1, 1, 1, 1) is a strict global minimum of the problem in Example 2.1 and |G| = 1 as U = {1, 2, 3} and V U = {4, 5, 6} is the only balanced partition of V leading to a cutset size of 1.

It appears from the existing literature that the term "reformulation" is almost never formally defined in the context of mathematical programming. The general consensus seems to be that given a formulation of an optimization problem, a reformulation is a different formulation having the same set of optima. Various authors make use of this definition without actually making it explicit, among which [107,114,101,[START_REF] Liberti | An exact reformulation algorithm for large nonconvex NLPs involving bilinear terms[END_REF][START_REF] Fischer | New constrained optimization reformulation of complementarity problems[END_REF][START_REF] Grant | Disciplined convex programming[END_REF][START_REF] Bjorkqvist | Automated reformulation of disjunctive constraints in MINLP optimization[END_REF][START_REF] Plateau | Reformulations quadratiques convexes pour la programmation quadratique en variables 0-1[END_REF][START_REF] Judice | Reformulation of mathematical programming problems as linear complementarity problems and investigation of their solution methods[END_REF][START_REF] Di Giacomo | Mathematical programming methods in dynamical nonlinear stochastic supply chain management[END_REF]. Many of the proposed reformulations, however, stretch this implicit definition somewhat. Liftings, for example (which consist in adding variables to the problem formulation), usually yield reformulations where an optimum in the original problem is mapped to a set of optima in the reformulated problem (see Sect. 3.2). Furthermore, it is sometimes noted how a reformulation in this sense is overkill because the reformulation only needs to hold at global optimality [START_REF] Adams | Comparisons and enhancement strategies for linearizing mixed 0-1 quadratic programs[END_REF]. Furthermore, reformulations sometimes really refer to a change of variables, as is the case in [START_REF] Mladenović | Reformulation descent applied to circle packing problems[END_REF]. Throughout the rest of this section we give various definitions for the concept of reformulation, and we explore the relations between them. We consider two problems P = (P(P ), V(P ), E(P ), O(P ), C(P ), B(P ), T (P ))

Q = (P(Q), V(Q), E(Q), O(Q), C(Q), B(Q), T (Q)).
Reformulations have been formally defined in the context of optimization problems (which are defined as decision problems with an added objective function). As was noted in Sect. 1, we see mathematical programming as a language used to describe and eventually solve optimization problems, so the difference is slim. The following definition is found in [START_REF] Audet | Links between linear bilevel and mixed 0-1 programming problems[END_REF]. Definition 6. Let P A and P B be two optimization problems. A reformulation B(•) of P A as P B is a mapping from P A to P B such that, given any instance A of P A and an optimal solution of B(A), an optimal solution of A can be obtained within a polynomial amount of time.

This definition is directly inspired to complexity theory and NP-completeness proofs. In the more practical and implementation oriented context of this chapter, Defn. 6 has one weak point, namely that of polynomial time. In practice, depending on the problem and on the instance, a polynomial time reformulation may just be too slow; on the other hand, Defn. 6 may bar a non-polynomial time reformulation which might be actually carried out within a practically reasonable amount of time. Furthermore, a reformulation in the sense of Defn. 6 does not necessarily preserve local optimality or the number of global optima, which might in some cases be a desirable reformulation feature. It should be mentioned that Defn. 6 was proposed in a paper that was more theoretical in nature, using an algorithmic equivalence between problems in order to attempt to rank equivalent NP-hard problems by their Branch-and-Bound solution difficulty.

The following definition was proposed by H. Sherali [105].

Definition 7. A problem Q is a reformulation of P if:
• there is a bijection σ :

F (P ) → F (Q); • |O(P )| = |O(Q)|;
• for all p = (e p , d p ) ∈ O(P ), there is a q = (e q , d q ) ∈ O(Q) such that e q = f (e p ) where f is a monotonic univariate function.

Defn. 7 imposes a very strict condition, namely the bijection between feasible regions of the original and reformulated problems. Although this is too strict for many useful transformations to be classified as reformulations, under some regularity conditions on σ it presents some added benefits, such as e.g. allowing easy correspondences between partitioned subspaces of the feasible regions and mapping sensitivity analysis results from reformulated to original problem.

In the rest of the section we discuss alternative definitions which only make use of the concept of optimum (also see [START_REF] Liberti | Reformulations in mathematical programming: Definitions[END_REF][START_REF] Liberti | Reformulations in mathematical programming: Definitions and systematics[END_REF]). These encompass a larger range of transformations as they do not require a bijection between the feasible regions, the way Defn. 7 does. Definition 8. Q is a local reformulation of P if there is a function ϕ : F (Q) → F (P ) such that (a) ϕ(y) ∈ L(P ) for all y ∈ L(Q), (b) ϕ restricted to L(Q) is surjective. This relation is denoted by P ≺ ϕ Q.

Informally, a local reformulation transforms all (local) optima of the original problem into optima of the reformulated problem, although more than one reformulated optimum may correspond to the same original optimum. A local reformulation does not lose any local optimality information and makes it possible to map reformulated optima back to the original ones; on the other hand, a local reformulation does not keep track of globality: some global optima in the original problem may be mapped to local optima in the reformulated problem, or vice-versa (see Example 2). x + sin(x) and Q ≡ min x∈ [-2π,2π] sin(x). It is easy to verify that there is a bijection between the local optima of P and those of Q. However, although P has a unique global optimum, every local optimum in Q is global.

Definition 9. Q is a global reformulation of P if there is a function ϕ : F (Q) → F (P ) such that (a) ϕ(y) ∈ G(P ) for all y ∈ G(Q), (b) ϕ restricted to G(Q) is surjective. This relation is denoted by P ¡ ϕ Q.
Informally, a global reformulation transforms all global optima of the original problem into global optima of the reformulated problem, although more than one reformulated global optimum may correspond to the same original global optimum. Global reformulations are desirable, in the sense that they make it possible to retain the useful information about the global optima whilst ignoring local optimality. At best, given a difficult problem P with many local minima, we would like to find a global reformulation Q where

L(Q) = G(Q).
Example 3. Consider a problem P with O(P ) = {f }. Let Q be a problem such that O(Q) = { f } and F (Q) = conv(F (P )), where conv(F (P )) is the convex hull of the points of F (P ) and f is the convex envelope of f over the convex hull of F (P ) (in other words, f is the greatest convex function underestimating f on F (P )). Since the set of global optima of P is contained in the set of global optima of Q [START_REF] Horst | On the convexification of nonlinear programming problems: an applications-oriented approach[END_REF], the convex envelope is a global reformulation.

Unfortunately, finding convex envelopes in explicit form is not easy. A considerable amount of work exists in this area: e.g. for bilinear terms [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF][START_REF] Al-Khayyal | Jointly constrained biconvex programming[END_REF], trilinear terms [START_REF] Meyer | Trilinear monomials with mixed sign domains: Facets of the convex and concave envelopes[END_REF], fractional terms [122], monomials of odd degree [START_REF] Liberti | Convex envelopes of monomials of odd degree[END_REF][START_REF] Liberti | Comparison of convex relaxations for monomials of odd degree[END_REF] the envelope is known in explicit form (this list is not exhaustive). See [119] for recent theoretical results and a rich bibliography. Definition 10. Q is an opt-reformulation (or exact reformulation) of P (denoted by P < Q) if there is a function ϕ : F (Q) → F (P ) such that P ≺ ϕ Q and P ¡ ϕ Q. This type of reformulation preserves both local and global optimality information, which makes it very attractive. Even so, Defn. 10 fails to encompass those problem transformations that eliminate some global optima whilst ensuring that at least one global optimum is left. Such transformations are specially useful in Integer Programming problems having several symmetric optimal solutions: restricting the set of global optima in such cases may be beneficial. One such example is the pruning of Branch-and-Bound regions based on the symmetry group of the problem presented in [START_REF] Margot | Pruning by isomorphism in branch-and-cut[END_REF]: the set of cuts generated by the procedure fails in general to be a global reformulation in the sense of Defn. 9 because the number of global optima in the reformulated problem is smaller than that of the original problem.

Lemma 1. The relations ≺, ¡, < are reflexive and transitive, but in general not symmetric.

Proof. For reflexivity, simply take ϕ as the identity. For transitivity, let P ≺ Q ≺ R with functions ϕ : F (Q) → F (P ) and ψ : F (R) → F (Q). Then ϑ = ϕ • ψ has the desired properties. In order to show that ≺ is not symmetric, consider a problem P with variables x and a unique minimum x * and a problem Q which is exactly like P but has one added variable w ∈ [0, 1]. It is easy to show that P ≺ Q (take ϕ as the projection of (x, w) on x). However, since for all w ∈ [0, 1] (x * , w) is an optimum of Q, there is no function of a singleton to a continuously infinite set that is surjective, so Q ≺ P . Given a pair of problems P, Q where ≺, ¡, < are symmetric on the pair, we call Q a symmetric reformulation of P . We remark also that by Lemma (1) we can compose elementary reformulations together to create chained reformulations (see Sect. 3.5 for examples).

Definition 11. Any problem Q that is related to a given problem P by a formula f (Q, P ) = 0 where f is a computable function is called an auxiliary problem with respect to P .

Deriving the formulation of an auxiliary problem may be a hard task, depending on f . The most useful auxiliary problems are those whose formulation can be derived algorithmically in time polynomial in |P |.

We remark that casting a problem in a standard form is an opt-reformulation. A good reformulation framework should be aware of the available solution algorithms and attempt to reformulate given problems into the most appropriate standard form.

Elementary reformulations

In this section we introduce some elementary reformulations in the proposed framework.

Objective function direction

Given an optimization problem P , the optimization direction d o of any objective function o ∈ O(P ) can be changed by simply setting d o ← -d o . This is an opt-reformulation where ϕ is the identity, and it rests on the identity min f (x) =max -f (x). We denote the effect of this reformulation carried out on each objective of a set O ⊆ O by ObjDir(P, O).

Constraint sense

Changing constraint sense simply means to write a constraint c expressed as e c ≤ b c as -e c ≥ -b c , or e c ≥ b c as -e c ≤ -b c . This is sometimes useful to convert the problem formulation to a given standard form. This is an opt-reformulation where ϕ is the identity. It can be carried out on the formulation by setting χ(r(e c )) ← -χ(r(e c )), s c ← -s c and b c = -b c . We denote the effect of this reformulation carried out for all constraints in a given set C ⊆ C by ConSense(P, C).

Liftings, restrictions and projections

We define here three important classes of auxiliary problems: liftings, restrictions and projections. Essentially, a lifting is the same as the original problem but with more variables. A restriction is the same as the original problem but with some of the variables replaced by either parameters or constants. A projection is the same as the original problem projected onto fewer variables. Whereas it is possible to give definitions of liftings and restrictions in terms of symbolic manipulations to the data structure given in Defn. 1, such a definition is in general not possible for projections. Projections and restrictions are in general not opt-reformulations nor reformulations in the sense of Defn. 7.

Lifting

A lifting Q of a problem P is a problem such that:

P(Q) P(P ), V(Q) V(P ), O(Q) = O(P ), E(Q) E(P ), C(Q) = C(P ), B(Q) B(P ), T (Q) T (P )
. This is an opt-reformulation where ϕ is a projection operator from V(Q) onto V(P ): for y ∈ F(Q), let ϕ(y) = (y(v) | v ∈ V(P )). We denote the lifting with respect to a new set of variables V by Lift(P, V ).

Essentially, a lifting is obtained by adding new variables to an optimization problem.

Restriction

A restriction Q of a problem P is such that:

• P(Q) ⊇ P(P ) • V(Q) V(P ) • |O(Q)| = |O(P )| • |C(Q)| = |C(P )|
• for each e ∈ E(P ) there is e ′ ∈ E(Q) such that e ′ is the same as e with any leaf node v ∈ V(P ) V(Q) replaced by an element of P(Q) ∪ R.

We denote the restriction with respect to a sequence of variable V with a corresponding sequence of values R by Restrict(P, V, R).

Essentially, a restriction is obtained by fixing some variables at corresponding given values.

Projection

A projection Q of a problem P is such that:

• P(Q) ⊇ P(P ) • V(Q) V(P ) • E, O, C, B, T (Q) are so that for all y ∈ F(Q) there is x ∈ F(P ) such that x(v) = y(v) for all v ∈ V(Q).
In general, symbolic algorithms to derive projections depend largely on the structure of the expression trees in E. If E consists entirely of linear forms, this is not difficult (see e.g. [START_REF] Balas | Projection, lifting and extended formulation in integer and combinatorial optimization[END_REF], Thm. 1.1). We denote the projection onto a set of variables V = V(Q) as Proj(P, V ). Essentially, F (Q) = {y | ∃x (x, y) ∈ F(P )}.

Equations to inequalities

Converting equality constraints to inequalities may be useful to conform to a given standard form. Suppose P has an equality constraint c = (e c , 0, b c ). This can be reformulated to a problem Q as follows:

• add two constraints c 1 = (e c , -1, b c ) and c 2 = (e c , 1, b c ) to C; • remove c from C.
This is an opt-reformulation denoted by Eq2Ineq(P, c).

Essentially, we replace the constraint e c = b c by the two constraints e c ≤ b c , e c ≥ b c .

Inequalities to equations

Converting inequalities to equality constraints is useful to convert problems to a given standard form: a very well known case is the standard form of a Linear Programming problem for use with the simplex method. Given a constraint c expressed as e c ≤ b c , we can transform it into an equality constraint by means of a lifting operation and a simple symbolic manipulation on the expression tree e c , namely:

• add a variable v c to V(P ) with interval bounds B(v c ) = [0, +∞] (added to B(P )) and type T (v c ) = 0 (added to T (P )); • add a new root node r 0 corresponding to the operator + (sum) to e c = (V, A), two arcs (r 0 , r(e c )), (r 0 , v) to A, and we then set r(e c ) ← r 0 ; • set s c ← 0.

We denote this transformation carried out on the set of constraints C by Slack(P, C). Naturally, for original equality constraints, this transformation is defined as the identity.

Performing this transformation on any number of inequality constraints results into an optreformulation.

Proposition 1. Given a set of constraints C ⊆ C(P ), the problem Q = Slack(P, C) is an optreformulation of P .

Proof. We first remark that V(P ) ⊆ V(Q). Consider ϕ defined as follows:

for each y ∈ F(Q) let ϕ(y) = x = (y(v) | v ∈ V(P )).
It is then easy to show that ϕ satisfies Defn. 10.

Absolute value terms

Consider a problem P involving a term e = (V, A) ∈ E where r(e) is the absolute value operator | • | (which is continuous but not differentiable everywhere); since this operator is unary, there is a single expression node f such that (r(e), f ) ∈ A. This term can be reformulated so that it is differentiable, as follows:

• add two continuous variables t + , t -with bounds [0, +∞];

• replace e by t + + t -;

• add constraints (f -t + -t -, 0, 0) and (t + t -, 0, 0) to C.

This is an opt-reformulation denoted by AbsDiff(P, e).

Essentially, we replace all terms |f | in the problem by a sum t + + t -, and then add the constraints f = t + + t -and t + t -= 0 to the problem.

Product of exponential terms

Consider a problem P involving a product g = i≤k h i of exponential terms, where h i = e fi for all i ≤ k. This term can be reformulated as follows:

• add a continuous variable w to V with T (w) = 0 and bounds B(w) = [0, +∞];

• add a constraint c = (e c , 0, 0) where e c = i≤k f ilog(w) to C; • replace g with w.

This is an opt-reformulation denoted by ProdExp(P, g). It is useful because many nonlinear terms (product and exponentials) have been the reduced to only one (the logarithm).

Essentially, we replace the product i e fi by an added nonnegative continuous variable w and then add the constraint log(w) = i f i to the problem.

Binary to continuous variables

Consider a problem P involving a binary variable x ∈ V with (T (x) = 2). This can be reformulated as follows:

• add a constraint c = (e c , 0, 0) to C where e c = x 2x;

• set T (x) = 0.
This is an opt-reformulation denoted by Bin2Cont(P, x). Since a binary variable x ∈ V can only take values in {0, 1}, any univariate equation in x that has exactly x = 0 and x = 1 as solutions can replace the binary constraint x ∈ {0, 1}. The most commonly used is the quadratic constraint x 2 = x, sometimes also written as

x(x -1) ≥ 0 ∧ x ≤ 1 [118].
In principle, this would reduce all binary problems to nonconvex quadratically constrained problems, which can be solved by a global optimization (GO) solver for nonconvex NLPs. In practice, GO solvers rely on an NLP subsolver to do most of the computationally intensive work, and NLP solvers are generally not very good in handling nonconvex/nonlinear equality constraints such as x 2 = x. This reformulation, however, is often used in conjunction with the relaxation of binary linear and quadratic problems (see Sect. 4.4).

Integer to binary variables

It is sometimes useful, for different reasons, to convert general integer variables to binary (0-1) variables. One example where this yields a crucial step into a complex linearization is given in Sect. 3.5. There are two established ways of doing this: one entails introducing binary assignment variables for each integer values that the variable can take; the other involves the binary representation of the integer variable value. Supposing the integer variable value is n, the first way employs O(n) added binary variables, whereas the second way only employs O(log 2 (n)). The first way is sometimes used to linearize complex nonlinear expressions of integer variables by transforming them into a set of constants to choose from (see example in Sect. 3.5). The second is often used in an indirect way to try and break symmetries in 0-1 problems: by computing the integer values of the binary representation of two 0-1 vectors x 1 , x 2 as integer problem variables v 1 , v 2 , we can impose ordering constraints such as v 1 ≤ v 2 + 1 to exclude permutations of x 1 , x 2 from the feasible solutions.

Assignment variables

Consider a problem P involving an integer variable v ∈ V with type T (v) = 1 and bounds

B(v) = [L v , U v ] such that U v -L v > 1. Let V = {L v , .
. . , U v } be the variable domain. Then P can be reformulated as follows:

• for all j ∈ V add a binary variable w j to V with T (w j ) = 2 and B(w j ) = [0, 1];

• add a constraint c = (e c , 0, 1) where e c = j∈V w j to C; • add an expression e = j∈V jw j to E; • replace all occurrences of v in the leaf nodes of expressions in E with e. This is an opt-reformulation denoted by Int2Bin(P, v).

Essentially, we add assignment variables w j = 1 if v = j and 0 otherwise. We then add an assignment constraint j∈V w j = 1 and replace v with j∈V jw j throughout the problem.

Binary representation

Consider a problem P involving an integer variable v ∈ V with type T (v) = 1 and bounds

B(v) = [L v , U v ] such that U v -L v > 1. Let V = {L v , .
. . , U v } be the variable domain. Then P can be reformulated as follows:

• let b be the minimum exponent such that |V | ≤ 2 b ; • add b binary variables w 1 , . . . , w b to V such that T (w j ) = 2 and B(w j ) = [0, 1] for all j ≤ b; • add an expression e = L v + j≤b w j 2 j
• replace all occurrences of v in the leaf nodes of expressions in E with e.

This is an opt-reformulation denoted by BinaryRep(P, v).

Essentially, we write the binary representation of v as L v + j≤b w j 2 j .

Feasibility to optimization problems

The difference between decision and optimization problems in computer science reflects in mathematical programming on the number of objective functions in the formulation. A formulation without objective functions models a feasibility problem; a formulation with one or more objective models an optimization problem. As was pointed out by the example in the introduction (see Sect. 1, p. 2), for computational reasons it is sometimes convenient to reformulate a feasibility problem in an optimization problem by introducing constraint tolerances. Given a feasibility problem P with O = ∅, we can reformulate it to an optimization problem Q as follows:

• add a large enough constant

M to P(Q); • add a continuous nonnegative variable ε to V(Q) with T (ǫ) = 0 and B(ǫ) = [0, M ]; • for each equality constraint c = (e c , 0, b c ) ∈ C, apply Eq2Ineq(P, c); • add the expression ε to E(Q); • add the objective function o = (ε, -1) to O(Q); • for each constraint c = (e c , s c , b c ) ∈ C (we now have s c = 0), let e ′ c = e c +s c ε and c ′ = (e ′ c , s c , b c ); add c ′ to C(Q).
As the original problem has no objective function, the usual definitions of local and global optima do not hold. Instead, we define any point in F (P ) to be both a local and a global optimum (see paragraph under Defn. 5). Provided the original problem is feasible, this is an opt-reformulation denoted by Feas2Opt(P ).

Proposition 2. Provided F (P ) = ∅, the reformulation Feas2Opt(P ) is an opt-reformulation.

Proof. Let F be the projection of F (Q) on the space spanned by the variables of P (i.e. all variables of Q but ε, see Sect. 3.2), and let π be the projection map. We then have F (P ) ⊆ F (this is because the constraints of Q essentially define a constraint relaxation of P , see Sect. 4.1 and Defn. 14). Let x ′ ∈ F(P ). We define ψ : F → F (P ) to be the identity on F (P ) and trivially extend it to F (Q) F by setting ψ(z) = x ′ for all z ∈ F(Q) F . The function φ = ψ • π maps F (Q) to F (P ), and preserves local minimality by construction, as per Defn. 8. Since ε is bounded below by zero, and the restriction (see Sect. 3.2) of Q to ε = 0 is exactly P , any x ′ ∈ G(Q) is also in F (P ). Moreover, by definition G(P ) = F (P ) as O(P ) = ∅, showing that the identity (projected on F ) preserves global minimality in the sense of Defn. 9.

Exact linearizations

Definition 12. An exact linearization of a problem P is an opt-reformulation Q of P where all expressions e ∈ E(P ) are linear forms.

Different nonlinear terms are linearized in different ways, so we sometimes speak of a linearization of a particular nonlinear term instead of a linearization of a given problem.

Piecewise linear objective functions

• add a continuous variable t to V bounded in [-∞, +∞]; • for all k ∈ K, add the constraint c k = (e k -t, d o , 0) to C.

This is an opt-reformulation denoted by MinMax(P ).

Essentially, we can reformulate an objective function min max k∈K e k as min t by adding a continuous variable t and the constraints ∀k ∈ K t ≥ e k to the problem.

Product of binary variables

Consider a problem P where one of the expressions e ∈ E(P ) is

k∈ K v k , where v k ∈ V(P ), B(v k ) = [0, 1] and T (v k ) = 2
for all k ∈ K (i.e. v k are binary 0-1 variables). This product can be linearized as follows:

• add a continuous variable w to V bounded in [0, 1]; • add the constraint ( k∈ K v k -w, -1, | K| -1) to C; • for all k ∈ K add the constraint (w -v k , -1, 0) to C.
This is an opt-reformulation denoted by ProdBin(P, K).

Essentially, a product of binary variables k∈ K v k can be replaced by an added continuous variable w ∈ [0, 1] and added constraints ∀k

∈ K w ≤ v k and w ≥ k∈ K v k -| K| + 1.
Proposition 3. Given a problem P and a set K ⊂ N, the problem Q = ProdBin(P, K) is an opt-reformulation of P .

Proof. Suppose first that ∀k ∈ K, v k = 1. We have to prove that w = 1 in that case. It comes from the hypothesis that k∈ K v k -| K| + 1 = 1 which implies by the last constraint that w = 1. The other constraints are all reduced to w ≤ 1 which are all verified. Suppose now that at least one of the binary variable is equal to zero and call i the index of this variable. Since ∀k ∈ K w ≤ v k , we have in particular the constraint for k = i. This leads to w = 0 which is the expected value. Besides, it comes that k∈ K v k -| K| ≤ -1. We deduced from this inequality that the last constraint is verified by the value of w.

As products of binary variables model the very common AND operation, linearizations of binary products are used very often. Hammer and Rudeanu [START_REF] Hammer | Boolean Methods in Operations Research and Related Areas[END_REF] cite [START_REF] Fortet | Applications de l'algèbre de Boole en recherche opérationelle[END_REF] as the first published appearance of the above linearization for cases where | K| = 2. For problems P with products v i v j for a given set of pairs {i, j} ∈ K where v i , v j are all binary variables, the linearization consists of |Q| applications of Prodbin(P, {i, j}) for each {i, j} ∈ K. Furthermore, we replace each squared binary variable v 2 i by simply v i (as v 2 i = v i for binary variables v i ). We denote this linearization by ProdSet(P, K).

Product of binary and continuous variables

Consider a problem P involving products v i v j for a given set K of ordered variable index pairs (i, j) where v i is a binary 0-1 variable and v j is a continuous variable with B(v j ) = [L j , U j ]. The problem can be linearized as follows:

• for all (i, j) ∈ K add a continuous variable w ij bounded by [L j , U j ] to V; • for all (i, j) ∈ K replace the product terms v i v j by the variable

w ij ; • for all (i, j) ∈ K add the constraints (w ij , -1, U j v i ), (w ij , 1, L j v i ), (w ij , -1, v j -(1 -v i )L j ), (w ij , 1, v j -(1 -v i )U j ) to C.
This is an opt-reformulation denoted by ProdBinCont(P, K).

Essentially, a product of a binary variable v i and a continuous variable v j bounded by [L j , U j ] can be replaced by an added variable w ij and added constraints:

       w ij ≤ U j v i w ij ≥ L j v i w ij ≤ v j -(1 -v i )L j w ij ≥ v j -(1 -v i )U j
Proposition 4. Given a problem P and a set K of ordered variable index pairs (i, j), the problem Q = ProdBinCont(P, K) is an opt-reformulation of P .

Proof. We have to prove that the reformulation ensures w ij = v i v j for all possible values for v i and v j . We do it by cases on the binary variable v i . Suppose first that v i = 0. Then the two first constraints implies that w ij = 0 which corresponds indeed to the product v i v j . It remains to see that the two other constraints don't interfere with this equality. In that case, the third constraint becomes w ij ≤ v j -L j . Since v j ≥ L j by definition, we have v j -L j ≥ 0 implying that w ij is less or equal to a positive term. With a similar reasoning, it comes from the fourth constraint that w ij is greater or equal to a negative term. Thus, for the case v i = 0, the constraints lead to w ij = 0. Suppose now that v i = 1. The two first inequalities lead to L i ≤ w ij ≤ U j which corresponds indeed to the range of the variable. The two last constraints become w ij ≥ v j and w ij ≤ v j . This implies w ij = v j which is the correct result.

Complementarity constraints

Consider a problem P involving constraints of the form c = (e c , 0, 0) where (a) r(e c ) is the sum operator, (b) for each node e outgoing from e c , e is a product operator, (c) each of these product nodes e has two outgoing nodes f, g. We can linearize such a constraint as follows:

• for each product operator node e outgoing from r(e c ) and with outgoing nodes f, g:

1. add a (suitably large) constant parameter M > 0 to P; 2. add a binary variable w to V with T (v) = 2 and B = [0, 1] 3. add the constraints (f -M w, -1, 0) and (g + M w, -1, M ) to C • delete the constraint c.

Provided we set M as an upper bound to the maximum values attainable by f and g, this is an opt-reformulation which is also a linearization. We denote it by CCLin(P ).

Essentially, we linearize complementarity constraints k∈K f k g k = 0 by eliminating the constraint, adding 0-1 variables w k for all k ∈ K and the linearization constraints f k ≤ M w k and g k ≤ M (1w k ). This reformulation, together with AbsDiff (see Sect. 3.2), provides an exact linearization (provided a suitably large but finite M exists) of absolute value terms.

Minimization of absolute values

• add continuous variables t + k , t - k with bounds [0, +∞]; • replace e k by t + k + t - k ; • add constraints (f k -t + k -t - k , 0, 0) to C.

This is an opt-reformulation denoted by MinAbs(P, K).

Essentially, we can reformulate an objective function min

k∈ K |f k | as min k∈ K (t + k + t - k ) whilst adding constraints ∀k ∈ K f k = t + k + t - k to the problem.
This reformulation is related to AbsDiff(P, e) (see Sect. 3.2), however the complementarity constraints t + k t - k = 0 are not needed because of the objective function direction: at a global optimum, because of the minimization of t + k + t - k , at least one of the variables will have value zero, thus implying the complementarity.

Linear fractional terms

Consider a problem P where an expression in E has a sub-expression e with a product operator and two subnodes e 1 , e 2 where ξ(e 1 ) = 1, ξ(e 2 ) = -1, and e 1 , e 2 are affine forms such that e 1 = i∈V a i v i + b and e 2 = i∈V c i v i + d, where v ⊆ V and T (v i ) = 0 for all i ∈ V (in other words e is a linear fractional term a ⊤ v+b c ⊤ v+d on continuous variables v). Assume also that the variables v only appear in some linear constraints of the problem Av = q (A is a matrix and q is a vector in P). Then the problem can be linearized as follows:

• add continuous variables α i , β to V (for i ∈ V ) with T (α i ) = T (β) = 0;
• replace e by i∈V a i α i + bβ;

• replace the constraints in Av = q by Aαqβ = 0;

• add the constraint i∈V c i α i + dβ = 1; • remove the variables v from V.
This is an opt-reformulation denoted by LinFract(P, e).

Essentially, α i plays the role of vi c ⊤ v+d , and β that of

1 c ⊤ v+d .
It is then easy to show that e can be re-written in terms of α, β as a ⊤ α + bβ, Av = q can be re-written as Aα = qβ, and that c ⊤ α + dβ = 1. Although the original variables v are removed from the problem, their values can be obtained by α, β after the problem solution, by computing v i = αi β for all i ∈ V .

Advanced reformulations

In this section we review a few advanced reformulations in the literature.

Hansen's Fixing Criterion

This method applies to unconstrained quadratic 0-1 problems of the form min x∈{0,1} n x ⊤ Qx where Q is an n × n matrix [START_REF] Hansen | Method of non-linear 0-1 programming[END_REF], and relies on fixing some of the variables to values guaranteed to provide a global optimum.

Let P be a problem with

P = {n ∈ N, {q ij ∈ R | 1 ≤ i, j ≤ n}}, V = {x i | 1 ≤ i ≤ n}, E = {f = i,j≤n q ij x i x j }, O = {(f, -1)}, C = ∅, B = [0, 1] n , T = 2.
This can be restricted (see Sect. 3.2) as follows:

• initialize two sequences V = ∅, A = ∅;
• for all i ≤ n:

1. if q ii + j<i min(0, q ij ) + j>i min(0, q ij ) > 0 then append x i to V and 0 to A; 2. (else) if q ii + j<i max(0, q ij ) + j>i max(0, q ij ) < 0 then append x i to V and 1 to A; • apply Restrict(P, V, A).

This opt-reformulation is denoted by FixQB(P ).

Essentially, any time a binary variable consistently decreases the objective function value when fixed, independently of the values of other variables, it is fixed.

Compact linearization of binary quadratic problems

This reformulation concerns a problem P with the following properties:

• there is a subset of binary variables x ⊆ V with |x| = n, T (x) = 2, B(x) = [0, 1] n ; • there is a set E = {(i, j) | 1 ≤ i ≤ j ≤ n} in P such that the terms x i x j appear as sub- expressions in the expressions E for all (i, j) ∈ E; • there is an integer K ≤ n in P and a covering {I k | k ≤ K} of {1, . . . , n} such that ( i∈I k x i , 0, 1) is in C for all k ≤ K; • there is a covering {J k | k ≤ K} of {1, . . . , n} such that I k ⊆ J k for all k ≤ K such that, letting F = {(i, j) | ∃k ≤ K((i, j) ∈ I k × J k ∨ (i, j) ∈ J k × I k )}, we have E ⊆ F .
Under these conditions, the problem P can be exactly linearized as follows:

• for all (i, j) ∈ F add continuous variables w ij with T (w ij ) = 0 and B(w ij ) = [0, 1];

• for all (i, j) ∈ E replace sub-expression x i x j with w ij in the expressions E;

• for all k ≤ K, j ∈ J k add the constraint ( i∈I k w ijx j , 0, 0) to C.

• for all (i, j) ∈ F add the constraint

w ij = w ji to C.
This opt-reformulation is denoted by RCLin(P, E). It was shown in [START_REF] Liberti | Compact linearization of binary quadratic problems[END_REF] that this linearization is exact and has other desirable tightness properties. See [START_REF] Liberti | Compact linearization of binary quadratic problems[END_REF] for examples.

Reduced RLT Constraints

This reformulation concerns a problem P with the following properties:

• there is a subset x ⊆ V with |x| = n and a set E = {(i, j) | 1 ≤ i ≤ j ≤ n} in P such that the terms x i x j appear as sub-expressions in the expressions E for all (i, j) ∈ E; • there is a number m ≤ n, an m × n matrix A = (a ij ) and an m-vector b in P such that ( j≤n a ij x j , 0, b i ) ∈ C for all i ≤ m.

Let F = {(i, j) | (i, j) ∈ E ∨ ∃k ≤ m(a kj = 0}.
Under these conditions, P can be reformulated as follows:

• for all (i, j) ∈ F add continuous variables

w ij with T (w ij ) = 0 and B(w ij ) = [-∞, +∞];
• for all (i, j) ∈ E replace sub-expression x i x j with w ij in the expressions E;

• for all i ≤ n, k ≤ m add the constraints ( j≤n a kj w ij -b k x i , 0, 0) to C: we call this linear system the Reduced RLT Constraint System (RCS) and ( j≤n a kj w ij , 0, 0) the companion system;

• let B = {(i, j) ∈ F | w ij is basic in the companion}; • let N = {(i, j) ∈ F | w ij is non-basic in the companion}; • add the constraints (w ij -x i x j , 0, 0) for all (i, j) ∈ N .
This opt-reformulation is denoted by RedCon(P ), and its validity was shown in [START_REF] Liberti | Linearity embedded in nonconvex programs[END_REF]. It is important because it effectively reduces the number of quadratic terms in the problem (only those corresponding to the set N are added). This reformulation can be extended to work with sparse sets E [START_REF] Liberti | An exact reformulation algorithm for large nonconvex NLPs involving bilinear terms[END_REF], namely sets E whose cardinality is small with respect to 1 2 n(n + 1). Essentially, the constraints w ij = x i x j for (i, j) ∈ B are replaced by the RCS ∀i ≤ n (Aw i = x i ), where w i = (w i1 , . . . , w in ).

Advanced examples

We give in this section a few advanced examples that illustrate the power of the elementary reformulations given above.

The Hyperplane Clustering Problem

As an example of what can be attained by combining these simple reformulations presented in this chapter, we give a MINLP formulation to the Hyperplane Clustering Problem (HCP) [START_REF] Dhyani | [END_REF][START_REF] Caporossi | Using the computer to study the dyamics of the handwriting processes[END_REF]. Given a set of points p = {p i | 1 ≤ i ≤ m} in R d we want to find a set of n hyperplanes w = {w j1 x 1 + . . . + w jd = w 0 j | 1 ≤ j ≤ n} in R d and an assignment of points to hyperplanes such that the distances from the hyperplanes to their assigned points are minimized.

We then derive a MILP reformulation. For clarity, we employ the usual mathematical notation instead of the notation given Defn. 1.

The problem P can be modelled as follows:

• Parameters. The set of parameters is given by p ∈ R m×d , m, n, d ∈ N.

• Variables. We consider the hyperplane coefficient variables w ∈ R n×d , the hyperplane constants w 0 ∈ R n , and the 0-1 assignment variables x ∈ {0, 1} m×n . • Objective function. We minimize the total distance, weighted by the assignment variable:

min i≤m j≤n |w j p i -w 0 j |x ij .
• Constraints. We consider assignment constraints: each point must be assigned to exactly one hyperplane:

∀i ≤ m j≤n x ij = 1,
and the hyperplanes must be nontrivial:

∀j ≤ n k≤d |w jk | = 1,
for otherwise the trivial solution with w = 0, w 0 = 0 would be optimal. This is a MINLP formulation because of the presence of the nonlinear terms (absolute values and products in the objective function) and of the binary assignment variables. We shall now apply several of the elementary reformulations presented in this chapter to obtain a MILP reformulation

Q of P . Let K = {(i, j) | i ≤ m, j ≤ n}.
1. Because x is nonnegative and because we are going to solve the reformulated MILP to global optimality, we can apply an reformulation similar to MinAbs(P, K) (see Sect. 3.3) to obtain an opt-reformulation P 1 as follows:

min i,j (t + ij x ij + t - ij x ij ) s.t. ∀i j x ij = 1 ∀j |w j |1 = 1 ∀i, j t + ij -t - ij = w j p i -w 0 j ,
where t + ij , t - ij ∈ [0, M ] are continuous added variables bounded above by a (large and arbitrary) constant M which we add to the parameter set P. We remark that this upper bound is enforced without loss of generality because w, w 0 can be scaled arbitrarily. 2. Apply ProdBinCont(P 1 , K) (see Sect. 3.3) to the products t + ij x ij and t - ij x ij to obtain a opt-reformulation P 2 as follows:

min i,j (y + ij + y - ij ) s.t. ∀i j x ij = 1 ∀j |w j |1 = 1 ∀i, j t + ij -t - ij = w j p i -w 0 j ∀i, j y + ij ≤ min(M x ij , t + ij ) ∀i, j y + ij ≥ M x ij + t + ij -M ∀i, j y - ij ≤ min(M x ij , t - ij ) ∀i, j y - ij ≥ M x ij + t - ij -M,
where y + ij , y - ij ∈ [0, M ] are continuous added variables. 3. For each term e jk = |w jk | apply AbsDiff(P 2 , e jk ) to obtain an opt-reformulation P 3 as follows:

min i,j (y + ij + y - ij ) s.t. ∀i j x ij = 1 ∀i, j t + ij -t - ij = w j p i -w 0 j ∀i, j y + ij ≤ min(M x ij , t + ij ) ∀i, j y + ij ≥ M x ij + t + ij -M ∀i, j y - ij ≤ min(M x ij , t - ij ) ∀i, j y - ij ≥ M x ij + t - ij -M ∀j k≤d (u + jk + u - jk ) = 1 ∀j, k u + jk -u - jk = w jk ∀j, k u + jk u - jk = 0, where u + jk , u - jk ∈ [0, M ]
are continuous variables for all j, k. Again, the upper bound does not enforce loss of generality. P 3 is an opt-reformulation of P : whereas P was not everywhere differentiable because of the absolute values, P 3 only involves differentiable terms. 4. We remark that the last constraints of P 3 are in fact complementarity constraints. We apply CCLin(P 3 ) to obtain the reformulated problem Q:

min i,j (y + ij + y - ij ) s.t. ∀i j x ij = 1 ∀i, j t + ij -t - ij = w j p i -w 0 j ∀i, j y + ij ≤ min(M x ij , t + ij ) ∀i, j y + ij ≥ M x ij + t + ij -M ∀i, j y - ij ≤ min(M x ij , t - ij ) ∀i, j y - ij ≥ M x ij + t - ij -M ∀j k≤d (u + jk + u - jk ) = 1 ∀j, k u + jk -u - jk = w jk ∀j, k u + jk ≤ M z jk ∀j, k u - jk ≤ M (1 -z jk )
, where z jk ∈ {0, 1} are binary variables for all j, k. Q is a MILP reformulation of P (see Sect. 2.3). This reformulation allows us to solve P by using a MILP solver -these have desirable properties with respect to MINLP solvers, such as numerical stability and robustness, as well as scalability and an optimality guarantee. A small instance consisting of 8 points and 2 planes in R 2 , with p = {(1, 7), (1, 1), (2, 2), (4, 3), (4, 5), (8, 3), (10, 1), [START_REF] Androulakis | alphaBB: A global optimization method for general constrained nonconvex problems[END_REF][START_REF] Adjiman | A global optimization method, αBB, for general twice-differentiable constrained NLPs: II. Implementation and computational results[END_REF]} is solved to optimality by the ILOG CPLEX solver [START_REF]ILOG: ILOG CPLEX 11.0 User's Manual[END_REF] to produce the following output:

Selection of software components

Large software systems consist of a complex architecture of interdependent, modular software components. These may either be built or bought off-the-shelf. The decision of whether to build or buy software components influences the cost, delivery time and reliability of the whole system, and should therefore be taken in an optimal way [START_REF] Cortellessa | Automated selection of software components based on cost/reliability tradeoff[END_REF].

Consider a software architecture with n component slots. Let I i be the set of off-the-shelf components and J i the set of purpose-built components that can be plugged in the i-th component slot, and assume I i ∩J i = ∅. Let T be the maximum assembly time and R be the minimum reliability level. We want to select a sequence of n off-the-shelf or purpose-built components compatible with the software architecture requirements that minimize the total cost whilst satisfying delivery time and reliability constraints. This problem can be modelled as follows.

• Parameters:

1. Let N ∈ N;

2. for all i ≤ n, s i is the expected number of invocations; 3. for all i ≤ n, j ∈ I i , c ij is the cost, d ij is the delivery time, and µ ij the probability of failure on demand of the j-th off-the-shelf component for slot i; 4. for all i ≤ n, j ∈ J i , cij is the cost, t ij is the estimated development time, τ ij the average time required to perform a test case, p ij is the probability that the instance is faulty, and b ij the testability of the j-th purpose-built component for slot i. • Variables:

1. Let x ij = 1 if component j ∈ I j ∪ J i is chosen for slot i ≤ n, and 0 otherwise; 2. Let N ij ∈ Z be the (non-negative) number of tests to be performed on the purpose-built component j ∈ J i for i ≤ n: we assume N ij ∈ {0, . . . , N }. • Objective function. We minimize the total cost, i.e. the cost of the off-the-shelf components c ij and the cost of the purpose-built components cij (

t ij + τ ij N ij ): min i≤n   j∈Ii c ij x ij + jinJi cij (t ij + τ ij N ij )x ij   .
• Constraints:

1. assignment constraints: each component slot in the architecture must be filled by exactly one software component ∀i ≤ n

j∈Ii∪Ji x ij = 1;
2. delivery time constraints: the delivery time for an off-the-shelf component is simply d ij , whereas for purpose-built components it is

t ij + τ ij N ij ∀i ≤ n j∈Ii d ij x ij + j∈Ji (t ij + τ ij N ij )x ij ≤ T ;
3. reliability constraints: the probability of failure on demand of off-the shelf components is µ ij , whereas for purpose-built components it is given by

ϑ ij = p ij b ij (1 -b ij ) (1-bij )Nij (1 -p ij ) + p ij (1 -b ij ) (1-bij )Nij ,
so the probability that no failure occurs during the execution of the i-th component is

ϕ i = e si P j∈I i µij xij + P j∈J i ϑij xij ! , whence the constraint is i≤n ϕ i ≥ R;
notice we have three classes of reliability constraints involving two sets of added variables ϑ, ϕ.

This problem is a MINLP with no continuous variables. We shall now apply several reformulations to this problem (call it P ).

1. Consider the term g = i≤n ϕ i and apply ProdExp(P, g) to P to obtain P 1 as follows:

min i≤n   j∈Ii c ij x ij + j∈Ji cij (t ij + τ ij N ij )x ij   ∀i ≤ n j∈Ii ∪Ji x ij = 1 ∀i ≤ n j∈Ii d ij x ij + j∈Ji (t ij + τ ij N ij )x ij ≤ T p ij b ij (1 -b ij ) (1-bij )Nij (1 -p ij ) + p ij (1 -b ij ) (1-bij )Nij = ϑ ij w ≥ R i≤n s i   j∈Ii µ ij x ij + j∈Ji ϑ ij x ij   = log(w),
and observe that w ≥ R implies log(w) ≥ log(R) because the log function is monotonically increasing, so the last two constraints can be grouped into a simpler one not involving logarithms of problem variables:

i≤n s i   j∈Ii µ ij x ij + j∈Ji ϑ ij x ij   ≥ log(R).
2. We now make use of the fact that N ij is an integer variable for all i ≤ n, j ∈ J i , and apply Int2Bin(P, N ij ). For k ∈ {0, . . . , N } we add assignment variables ν k ij so that ν k ij = 1 if N ij = k and 0 otherwise. Now for all k ∈ {0, . . . , N } we compute the constants

ϑ k = pij bij (1-bij ) (1-b ij )k (1-pij )+pij (1-bij ) (1-b ij )k
, which we add to the problem parameters. We remove the constraints defining ϑ ij in function of N ij : since the following constraints are valid:

∀i ≤ n, j ∈ J i k≤N ν k ij = 1 (12) ∀i ≤ n, j ∈ J i k≤N kν k ij = N ij ( 13
) ∀i ≤ n, j ∈ J i k≤N ϑ k ν k ij = ϑ ij , (14)
the second constraints are used to replace N ij and the third to replace ϑ ij . The first constraints are added to the formulation. We obtain:

min i≤n   j∈Ii c ij x ij + j∈Ji cij (t ij + τ ij k≤N kν k ij )x ij   ∀i ≤ n j∈Ii∪Ji x ij = 1 ∀i ≤ n j∈Ii d ij x ij + j∈Ji (t ij + τ ij k≤N kν k ij )x ij ≤ T i≤n s i   j∈Ii µ ij x ij + j∈Ji x ij k≤N ϑ k ν k ij   ≥ log(R) ∀i ≤ n, j ∈ J i k≤N ν k ij = 1.
3. We distribute products over sums in the formulation to obtain the binary product sets

{x ij ν k ij | k ≤ N }
for all i ≤ n, j ∈ J i : by repeatedly applying the ProdBin reformulation to all binary products of binary variables, we get a MILP opt-reformulation Q of P where all the variables are binary.

We remark that the MILP opt-reformulation Q derived above has a considerably higher cardinality than |P |. More compact reformulations are applicable in step 3 because of the presence of the assignment constraints (see Sect. 3.4).

Reformulation Q essentially rests on linearization variables w k ij which replace the quadratic terms x ij ν k ij throughout the formulation. A semantic interpretation of step 3 is as follows. We notice that for i ≤ n, j ∈ J i , if x ij = 1, then x ij = k ν k ij (because only one value k will be selected), and if x ij = 0, then x ij = k ν k ij (because no value k will be selected). This means that ∀i ≤ n, j

∈ J i x ij = k≤N ν k ij ( 15
)
is a valid problem constraint. We use it to replace x ij everywhere in the formulation where it appears with j ∈ I i , obtaining a opt-reformulation with x ij for j ∈ I i and quadratic terms ν k ij ν h lp . Now, because of (12), these are zero when (i, j) = (l, p) or k = h and are equal to ν k ij when (i, j) = (l, p) and k = h, so they can be linearized exactly by replacing them by either 0 or ν k ij according to their indices. What this really means is that the reformulation Q, obtained through a series of automatic reformulation steps, is a semantically different formulation defined in terms of the following decision variables:

∀i ≤ n, j ∈ I i x ij = 1 if j ∈ I i is assigned to i 0 otherwise. ∀i ≤ n, j ∈ J i , k ≤ N ν k ij = 1 if j ∈ J i
is assigned to i and there are k tests to be performed 0 otherwise. This is an important hint to the importance of automatic reformulation in problem analysis: it is a syntactical operation, the result of which, when interpreted, can suggest a new meaning.

Relaxations

Loosely speaking, a relaxation of a problem P is an auxiliary problem of P whose feasible region is larger; often, relaxations are obtained by simply removing constraints from the formulation.

Relaxations are useful because they often yield problems which are simpler to solve yet they provide a bound on the objective function value at the optimum.

Such bounds are mainly used in Branch-and-Bound type algorithms, which are the most common exact or ε-approximate (for a given ε > 0) solution algorithms for MILPs, nonconvex NLPs and MINLPs. Although the variants for solving MILPs, NLPs and MINLPs are rather different, they all conform to the same implicit enumeration search type. Lower and upper bounds are computed for the problem over the current variable domains. If the bounds are sufficiently close, a global optimum was found in the current domain: store it if it improves the incumbent (i.e. the current best optimum). Otherwise, partition the domain and recurse over each subdomain in the partition. Should a bound be worse off than the current incumbent during the search, discard the domain immediately without recursing on it. Under some regularity conditions, the recursion terminates. The Branch-and-Bound algorithm has been used on combinatorial optimization problems since the 1950s [START_REF] Aho | Data Structures and Algorithms[END_REF]. Its first application to nonconvex NLPs is [START_REF] Falk | An algorithm for separable nonconvex programming problems[END_REF]. More recently, Branch-and-Bound has evolved into Branch-and-Cut and Branch-and-Price for MILPs [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF]133,[START_REF]ILOG: ILOG CPLEX 11.0 User's Manual[END_REF], which have been used to solve some practically difficult problems such as the Travelling Salesman Problem (TSP) [START_REF] Applegate | The Travelling Salesman Problem: a Computational Study[END_REF]. Some recent MINLP-specific Branch-and-Bound approaches are [102, [START_REF] Androulakis | alphaBB: A global optimization method for general constrained nonconvex problems[END_REF][START_REF] Adjiman | A global optimization method, αBB, for general twice-differentiable constrained NLPs: I. Theoretical advances[END_REF][START_REF] Adjiman | A global optimization method, αBB, for general twice-differentiable constrained NLPs: II. Implementation and computational results[END_REF]114,124,[START_REF] Liberti | Writing global optimization software[END_REF].

A further use of bounds provided by mathematical programming formulations is to evaluate the performance of heuristic algorithms without an approximation guarantee [START_REF] Davidović | Towards the optimal solution of the multiprocessor scheduling problem with communication delays[END_REF]. Bounds are sometimes also used to guide heuristics [START_REF] Puchinger | Relaxation guided variable neighbourhood search[END_REF].

In this section we define relaxations and review the most useful ones. In Sect. 4.1 we give some basic definitions. We then list elementary relaxations in Sect. 4.2 and more advanced ones in Sect. 4.3. We discuss relaxation strengthening in Sect. 4.4.

Definitions

Consider an optimization problem P = (P, V, E, O, C, B, T ) and let Q be such that:

P(Q) ⊇ P(P ), V(Q) = V(P ), E(Q) ⊇ E(P ) and O(Q) = O(P ).
We first define relaxations in full generality.

Definition 13. Q is a relaxation of P if (a) F (P ) ⊆ F(Q); (b) for all (f, d) ∈ O(P ), ( f , d) ∈ O(Q) and x ∈ F(P ), d f (x) ≥ df (x).
Defn. 13 is not used very often in practice because it does not say anything on how to construct Q. The following elementary relaxations are more useful. Definition 14. Q is a:

• constraint relaxation of P if C(P ) C(Q); • bound relaxation of P if B(P ) B(Q); • a continuous relaxation of P if ∃v ∈ V(P ) (T (v) > 0) and T (v) = 0 for all v ∈ V(Q).

Elementary relaxations

We shall consider two types of elementary relaxations: the continuous relaxation and the convex relaxation. The former is applicable to MILPs and MINLPs, and the latter to (nonconvex) NLPs and MINLPs. They are both based on the fact that whereas solving MILPs and MINLPs is considered difficult, there are efficient algorithms for solving LPs and convex NLPs. Since the continuous relaxation was already defined in Defn. 14 and trivially consists in considering integer/discrete variables as continuous ones, in the rest of this section we focus on convex relaxations.

Formally (and somewhat obviously), Q is a convex relaxation of a given problem P if Q is a relaxation of P and Q is convex. Associated to all sBB in the literature there is a (nonconvex) NLP or MINLP in standard form, which is then used as a starting point for the convex relaxation.

Outer approximation

Outer approximation (OA) is a technique for defining a polyhedral approximation of a convex nonlinear feasible region, based on computing tangents to the convex feasible set at suitable boundary points [START_REF] Duran | An outer-approximation algorithm for a class of mixed-integer nonlinear programs[END_REF][START_REF] Fletcher | Solving mixed integer nonlinear programs by outer approximation[END_REF][START_REF] Kesavan | Outer-approximation algorithms for nonconvex mixed-integer nonlinear programs[END_REF]]. An outer approximation relaxation relaxes a convex NLP to an LP, (or a MINLP to a MILP) and is really a "relaxation scheme" rather than a relaxation: since the tangents to all boundary points of a convex set define the convex set itself, any choice of (finite) set of boundary points of the convex can be used to define a different outer approximation. OA-based optimization algorithms identify sets of boundary points that eventually guarantee that the outer approximation will be exact near the optimum. In [START_REF] Kesavan | Outer-approximation algorithms for nonconvex mixed-integer nonlinear programs[END_REF], the following convex MINLP is considered:

min L 0 (x) + cy s.t. L(x) + By ≤ 0 x L ≤ x ≤ x U y ∈ {0, 1} q ,        (16)
where L 0 : R n → R, L : R n → R m are convex once-differentiable functions, c ∈ R q , B is an m × q matrix. For a given y ′ ∈ {0, 1} q , let P (y ′ ) be ( 16) with y fixed at y ′ . Let {y j } be a sequence of binary q-vectors. Let T = {j | P (y j ) is feasible with solution x j }. Then the following is a MILP outer approximation for [START_REF] Balas | Canonical cuts on the unit hypercube[END_REF]:

min x,y,η η ∀j ∈ T L 0 (x j ) + ∇L 0 (x j )(x -x j ) + cy ≤ η ∀j L(x j ) + ∇L(x j )(x -x j ) + By ≤ 0 x L ≤ x ≤ x U y ∈ {0, 1} q ,           
where x j is the solution to F (y j ) (defined in [START_REF] Fletcher | Solving mixed integer nonlinear programs by outer approximation[END_REF]) whenever P (y j ) is infeasible. This relaxation is denoted by OuterApprox(P, T ).

αBB convex relaxation

The αBB algorithm [START_REF] Androulakis | alphaBB: A global optimization method for general constrained nonconvex problems[END_REF][START_REF] Adjiman | A global optimization method, αBB, for general twice-differentiable constrained NLPs: I. Theoretical advances[END_REF][START_REF] Adjiman | A global optimization method, αBB, for general twice-differentiable constrained NLPs: II. Implementation and computational results[END_REF][START_REF] Floudas | Deterministic Global Optimization[END_REF] targets single-objective NLPs where the expressions in the objective and constraints are twice-differentiable. The convex relaxation of the problem P :

min x f (x) s.t. g(x) ≤ 0 h(x) = 0 x L ≤ x ≤ x U        (17)
is obtained as follows.

1. Apply the Eq2Ineq reformulation (see Sect. a) if the root node r of the expression tree e c is a sum operator, for every subnode s ∈ δ + (r) replace s with a specialized convex underestimator if s is a bilinear, trilinear, linear fractional, fractional trilinear, univariate concave term. Otherwise replace with α-underestimator; b) otherwise, replace r with a specialized if s is a bilinear, trilinear, linear fractional, fractional trilinear, univariate concave term. Otherwise replace with α-underestimator.

The specialized underestimators are as follows: McCormick's envelopes for bilinear terms [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF][START_REF] Al-Khayyal | Jointly constrained biconvex programming[END_REF], the second-level RLT bound factor linearized products [108, 107, 104] for trilinear terms, and a secant underestimator for univariate concave terms. Fractional terms are dealt with by extending the bilinear/trilinear underestimators to bilinear/trilinear products of univariate functions and then noting that x/y = φ 1 (x)φ 2 (y) where φ 1 is the identity and φ 2 (y) = 1/y [START_REF] Maranas | Finding all solutions to nonlinearly constrained systems of equations[END_REF]. Recently, the convex underestimator for trilinear terms have been replaced with the convex envelopes [START_REF] Meyer | Trilinear monomials with mixed sign domains: Facets of the convex and concave envelopes[END_REF].

The general-purpose α-underestimator:

α(x L -x) ⊤ (x U -x) (18)
is a quadratic convex function that for suitable values of α is "convex enough" to overpower the generic nonconvex term. This occurs for α ≥ max{0, -1 2 min

x L ≤x≤x U λ(x)},
where min λ(x) is the minimum eigenvalue of the Hessian of the generic nonconvex term in function of the problem variables. The resulting αBB relaxation Q of P is a convex NLP. This relaxation is denoted by αBBRelax(P ).

Branch-and-Contract convex relaxation

The convex relaxation is used in the Branch-and-Contract algorithm [134], targeting nonconvex NLPs with twice-differentiable objective function and constraints. This relaxation is derived essentially in the same way as for the αBB convex relaxation. The differences are:

• the problem is assumed to only have inequality constraints of the form c = (e c , -1, 0); • each function (in the objective and constraints) consists of a sum of nonlinear terms including: bilinear, linear fractional, univariate concave, and generic convex.

The convex relaxation is then constructed by replacing each nonconvex nonlinear term in the objective and constraints by a corresponding envelope or relaxation. The convex relaxation for linear fractional term had not appeared in the literature before [134].

Symbolic reformulation based convex relaxation

This relaxation is used in the symbolic reformulation spatial Branch-and-Bound algorithm proposed in [113,114]. It can be applied to all NLPs and MINLPs for which a convex underestimator and a concave overestimator are available. It consists in reformulating P to the Smith standard form (see Sect. 2.3) and then replacing every defining constraint with the convex and concave under/over-estimators. In his Ph.D. thesis [112], Smith had tried both NLP and LP convex relaxations, finding that LP relaxations were more reliable and faster to compute, although of course with slacker bounds. The second implementation of the sBB algorithm he proposed is described in [START_REF] Liberti | Reformulation and convex relaxation techniques for global optimization[END_REF][START_REF] Liberti | Writing global optimization software[END_REF] and implemented in the ooOPS software framework [START_REF] Liberti | ooOPS[END_REF]. Both versions of this algorithm consider under/overestimators for the following terms: bilinear, univariate concave, univariate convex (linear fractional being reformulated to bilinear). The second version also included estimators for piecewise convex/concave terms. One notable feature of this relaxation is that it can be adapted to deal with more terms. Some recent work in polyhedral envelopes, for example [119], gives conditions under which the sum of the envelopes is the envelope of the sum: this would yield a convex envelope for a sum of terms. It would then suffice to provide for a defining constraint in the Smith standard form linearizing the corresponding sum. The Smith relaxation is optionally strengthened via LP-based optimality and feasibility based range reduction techniques. After every range reduction step, the convex relaxation is updated with the new variable ranges in an iterative fashion until no further range tightening occurs [112, [START_REF] Liberti | Reformulation and convex relaxation techniques for global optimization[END_REF][START_REF] Liberti | Writing global optimization software[END_REF]. This relaxation, denoted by SmithRelax(P ) is at the basis of the sBB solver [START_REF] Liberti | Writing global optimization software[END_REF] in the ooOPS software framework [START_REF] Liberti | ooOPS[END_REF], which was used to obtain solutions of many different problem classes: pooling and blending problems [START_REF] Haverly | Studies of the behaviour of recursion for the pooling problem[END_REF][START_REF] Liberti | An exact reformulation algorithm for large nonconvex NLPs involving bilinear terms[END_REF], distance geometry problems [START_REF] Kucherenko | New formulations for the kissing number problem[END_REF][START_REF] Lavor | Computational experience with the molecular distance geometry problem[END_REF], and a quantum chemistry problem [START_REF] Lavor | Solving Hartree-Fock systems with global optimization metohds[END_REF][START_REF] Liberti | Reformulation in mathematical programming: an application to quantum chemistry[END_REF].

BARON's convex relaxation

BARON (Branch And Reduce Optimization Navigator) is a commercial Branch-and-Bound based global optimization solver (packaged within the GAMS [START_REF] Brook | GAMS, a user's guide[END_REF] modelling environment) which is often quoted as being the de facto standard solver for MINLPs [124,123]. Its convex relaxation is derived essentially in the same way as for the symbolic reformulation based convex relaxation. The differences are:

• better handling of fractional terms [120, 121] • advanced range reduction techniques (optimality, feasibility and duality based, plus a learning reduction heuristic) • optionally, an LP relaxation is derived via outer approximation.

Advanced relaxations

In this section we shall describe some more advanced relaxations, namely the Lagrangian relaxation, the semidefinite relaxation, the reformulation-linearization technique and the signomial relaxation.

Lagrangian relaxation

Consider a MINLP f * = min x f (x) s.t. g(x) ≤ 0 x ∈ X ⊆ R n ,    ( 19
)
where f : R n → R and g : R n → R m are continuous functions and X is an arbitrary set. The Lagrangian relaxation consists in "moving" the weighted constraints to the objective function, namely:

L(µ) = inf x f (x) + µ ⊤ g(x) x ∈ X ⊆ R n ,
for some nonnegative µ ∈ R m + . For all x ∈ X with g(x) ≤ 0, we have µ ⊤ g(x) ≤ 0, which implies L(µ) ≤ f * for all µ ≥ 0. In other words, L(µ) provides a lower bound to (19) for all µ ≥ 0. Thus, we can improve the tightness of the relaxation by solving the Lagrangian problem

max µ≥0 L(µ), (20)
(namely, we attempt to find the largest possible lower bound). If ( 19) is an LP problem, it is easy to show that the Lagrangian problem [START_REF] Bjorkqvist | Automated reformulation of disjunctive constraints in MINLP optimization[END_REF] is the dual LP problem. In general, solving [START_REF] Bjorkqvist | Automated reformulation of disjunctive constraints in MINLP optimization[END_REF] is not a computationally easy task [START_REF] Nowak | Relaxation and Decomposition Methods for Mixed Integer Nonlinear Programming[END_REF]. However, one of the nice features of Lagrangian relaxations is that they provide a lower bound for each value of µ ≥ 0, so [START_REF] Bjorkqvist | Automated reformulation of disjunctive constraints in MINLP optimization[END_REF] does not need to be solved at optimality. Another useful feature is that any subset of problem constraints can be relaxed, for X can be defined arbitrarily. This is useful for problems that are almost block-separable, i.e. those problems that can be decomposed in some independent subproblems bar a few constraints involving all the problem variables (also called complicating constraints). In these cases, one considers a Lagrangian relaxation of the complicating constraints and then solves a block-separable Lagrangian problem. This approach is called Lagrangian decomposition.

The Lagrangian relaxation has some interesting theoretical properties: (a) for convex NLPs it is a global reformulation [START_REF] Boyd | Convex Optimization[END_REF]; (b) for MILPs, it is at least as tight as the continuous relaxation [133]; (c) for MINLPs, under some conditions (i.e. some constraint qualification and no equality constraints) it is at least as tight as any convex relaxation obtained by relaxing each nonconvex term or each constraint one by one [START_REF] Horst | Duality bound methods in global optimization[END_REF], such as all those given in Sect. 4.2. Further material on the use of Lagrangian relaxation in NLPs and MINLPs can be found in [START_REF] Nowak | Relaxation and Decomposition Methods for Mixed Integer Nonlinear Programming[END_REF][START_REF] Horst | Duality bound methods in global optimization[END_REF].

Consider a problem P such that O(P ) = {(e o , d o )} and a subset of constraints C ⊆ C(P ). A Lagrangian relaxation of C in P (denoted by LagRel(P, C)) is a problem Q defined as follows.

• V(Q) = V(P ), B(Q) = B(P ), T (Q) = T (P ), • P(Q) = P(P ) ∪ {µ c | c ∈ C}, • C(Q) = C(P ) C, • O(Q) = {(e ′ o , d ′ o )}
, where e ′ o = e o + c∈C µ c c. The Lagrangian problem cannot itself be defined in the data structure of Defn. 1, for the max operator is only part of O L as long as it has a finite number of arguments.

Semidefinite relaxation

As was pointed out in Sect. 2.3, SDPs provide very tight relaxations for quadratically constrained quadratic MINLPs (QCQP). A QCQP in general form is as follows [START_REF] Anstreicher | SDP versus RLT for nonconvex QCQPs[END_REF]:

min x x ⊤ Q 0 x + a ⊤ 0 x ∀i ∈ I x ⊤ Q i x + a ⊤ i x ≤ b i ∀i ∈ E x ⊤ Q i x + a ⊤ i x = b i x L ≤ x ≤ x U ∀j ∈ J x i ∈ Z,            (21)
where

I ∪ E = {1, . . . , m}, J ⊆ {1, . . . , n}, x ∈ R n , Q i is an n × n symmetric matrix for all i ≤ m.
For general matrices Q i and J = ∅, the QCQP is nonconvex. Optionally, the integer variables can be reformulated exactly to binary (see Int2Bin, Sect. 3.2) and subsequently to continuous (see Bin2Cont, Sect. 3.2) via the introduction of the constraints x 2 ix i = 0 for all i ∈ J: since these constraints are quadratic, they can be accommodated in formulation ( 21) by suitably modifying the Q i matrices. Many important applications can be modelled as QCQPs, including graph bisection (see Sect. 2.1) and graph partitioning [START_REF] Liberti | Compact linearization of binary quadratic problems[END_REF], scheduling with communication delays [START_REF] Davidović | Towards the optimal solution of the multiprocessor scheduling problem with communication delays[END_REF], distance geometry problems such as the KNP (see Sect. 2.1) [START_REF] Kucherenko | New formulations for the kissing number problem[END_REF] and the Molecular Distance Geometry Problem (MDGP) [START_REF] Lavor | Computational experience with the molecular distance geometry problem[END_REF][START_REF] Liberti | Double VNS for the molecular distance geometry problem[END_REF], pooling and blending problems from the oil industry [START_REF] Haverly | Studies of the behaviour of recursion for the pooling problem[END_REF][START_REF] Liberti | An exact reformulation algorithm for large nonconvex NLPs involving bilinear terms[END_REF] and so on.

The SDP relaxation of the QCQP, denoted by SDPRelax(P ) is constructed as follows:

• replace all quadratic products x i x j in ( 21) with an added linearization variable X ij • form the matrix X = (X ij ) and the variable matrix

X = 1 x ⊤ x X • for all 0 ≤ i ≤ m form the matrices Qi = -b i a ⊤ i /2 a i /2 Q i
• the following is an SDP relaxation for QCQP:

min X Q0 • X ∀i ∈ I Qi • X ≤ 0 ∀i ∈ E Qi • X = 0 x L ≤ x ≤ x U X 0.            (22)
As for the SDP standard form of Sect. 2.3, the SDP relaxation can be easily represented by the data structure described in Defn. 1.

Reformulation-Linearization Technique

The Reformulation-Linearization Technique (RLT) is a relaxation method for mathematical programming problems with quadratic terms. The RLT linearizes all quadratic terms in the problem and generates valid linear equation and inequality constraints by considering multiplications of bound factors (terms like x ix L i and x U ix i ) and constraint factors (the left hand side of a constraint such as n j=1 a j x jb ≥ 0 or n j=1 a j x jb = 0). Since bound and constraint factors are always non-negative, so are their products: this way one can generate sets of valid problem constraints. In a sequence of papers published from the 1980s onwards (see e.g. [START_REF] Adams | A tight linearization and an algorithm for 0-1 quadratic programming problems[END_REF]108,110,107,103,111, 109]), RLT-based relaxations were derived for many different classes of problems, including IPs, NLPs, MINLPs in general formulation, and several real-life applications. It was shown that the RLT can be used in a lift-and-project fashion to generate the convex envelope of binary and general discrete problems [106,[START_REF] Adams | A hierarchy of relaxations leading to the convex hull representation for general discrete optimization problems[END_REF].

Basic RLT

The RLT consists of two symbolic manipulation steps: reformulation and linearization. The reformulation step is a reformulation in the sense of Defn. 10. Given a problem P , the reformulation step produces a reformulation Q ′ where:

• P(Q ′ ) = P(P ); • V(Q ′ ) = V(P ); • E(Q ′ ) ⊇ E(P ); • C(Q ′ ) ⊇ C(P ); • O(Q ′ ) = O(P ); • B(Q ′ ) = B(P ); • T (Q ′ ) = T (P ); • ∀x, y ∈ V(P ), add the following constraints to C(Q ′ ): (x -L x )(y -L y ) ≥ 0 (23) (x -L x )(U y -y) ≥ 0 (24) (U x -x)(y -L y ) ≥ 0 ( 25
) (U x -x)(U y -y) ≥ 0; (26)
• ∀x ∈ V(P ), c = (e c , s c , b c ) ∈ C(P ) such that e c is an affine form, s c = 1 and b c = 0 (we remark that all linear inequality constraints can be easily reformulated to this form, see Sect. 3.2), add the following constraints to C(Q ′ ):

e c (x -L x ) ≥ 0 (27) e c (U x -x) ≥ 0; (28)
• ∀x ∈ V(P ), c = (e c , s c , b c ) ∈ C(P ) such that e c is an affine form, s c = 0 and b c = 0 (we remark that all linear equality constraints can be trivially reformulated to this form), add the following constraint to C(Q ′ ):

e c x = 0. ( 29
)
Having obtained Q ′ , we proceed to linearize all the quadratic products engendered by ( 23)- [START_REF] Dhyani | [END_REF]. We derive the auxiliary problem Q from Q ′ by reformulating Q ′ to Smith's standard form (see Sect. 2.3) and then performing a constraint relaxation with respect to all defining constraints. Smith's standard form is a reformulation of the lifting type, and the obtained constraint relaxation Q is a MILP whose optimal objective function value f is a bound to the optimal objective function value f * of the original problem P . The bound obtained in this way is shown to dominate, or be equivalent to, several other bounds in the literature [START_REF] Adams | A hierarchy of relaxations leading to the convex hull representation for general discrete optimization problems[END_REF]. This relaxation is denoted by RLTRelax(P ).

We remark in passing that ( 23)-( 26), when linearized by replacing the bilinear term xy with an added variable w, are also known in the literature as McCormick relaxation, as they were first proposed as a convex relaxation of the nonconvex constraint w = xy [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF], shown to be the convex envelope [START_REF] Al-Khayyal | Jointly constrained biconvex programming[END_REF], and widely used in spatial Branch-and-Bound (sBB) algorithms for global optimization [114,[START_REF] Adjiman | A global optimization method, αBB, for general twice-differentiable constrained NLPs: I. Theoretical advances[END_REF][START_REF] Adjiman | A global optimization method, αBB, for general twice-differentiable constrained NLPs: II. Implementation and computational results[END_REF]124,[START_REF] Liberti | Writing global optimization software[END_REF]. RLT constraints of type [START_REF] Dhyani | [END_REF] have been the object of further research showing their reformulating power [START_REF] Liberti | Reduction constraints for the global optimization of NLPs[END_REF][START_REF] Liberti | Reformulation and convex relaxation techniques for global optimization[END_REF][START_REF] Liberti | Linearity embedded in nonconvex programs[END_REF][START_REF] Liberti | An exact reformulation algorithm for large nonconvex NLPs involving bilinear terms[END_REF][START_REF] Liberti | Compact linearization of binary quadratic problems[END_REF] (also see Sect 3.4, where we discuss compact linearization of binary quadratic problems and reduced RLT constraints).

RLT Hierarchy

The basic RLT method can be extended to provide a hierarchy of relaxations, by noticing that we can form valid RLT constraints by multiplying sets of bound and constraint factors of cardinality higher than 2, and then projecting the obtained constraints back to the original variable space. In [106,[START_REF] Adams | A hierarchy of relaxations leading to the convex hull representation for general discrete optimization problems[END_REF] it is shown that this fact can be used to construct the convex hull of an arbitrary MILP P . For simplicity, we only report the procedure for MILP in standard canonical form (see Sect. 2.3) where all discrete variables are binary, i.e. T (v) = 2 for all v ∈ V(P ). Let |V(P )| = n. For all integer d ≤ n, let P d be the relaxation of P obtained as follows:

• for all linear constraint c = (e c , 1, 0) ∈ C(P ), subset V ⊆ V(P ) and finite binary sequence B with |V | = |B| = d such that B x is the x-th term of the sequence for x ∈ V , add the valid constraint:

e c    x∈V Bx=0 x       x∈V Bx=1 (1 -x)    ≥ 0; (30)
we remark that ( 30) is a multivariate polynomial inequality; • for all monomials of the form a x∈J⊆V(P )

x with a ∈ R in a constraint [START_REF] Di Giacomo | Mathematical programming methods in dynamical nonlinear stochastic supply chain management[END_REF]

(P )) ⊆ F(X n ) ⊆ F(X n-1 ) . . . ⊆ F(X 1 ) ⊆ F(P ).
We recall that for a set Y ⊆ R n , conv(Y ) is defined as the smallest convex subset of R n containing Y .

A natural practical application of the RLT hierarchy is to generate relaxations for polynomial programming problems [103], where the various multivariate monomials generated by the RLT hierarchy might already be present in the problem formulation.

Signomial programming relaxations

A signomial programming problem is an optimization problem where every objective function is a signomial function and every constraint is of the form c = (g, s, 0) where g is a signomial function of the problem variables, and s = 0 (so signomial equality constraints must be reformulated to pairs of inequality constraints as per the Eq2Ineq reformulation of Sect. 3.2). A signomial is a term of the form:

a K k=1 x r k k , (31)
where a, r k ∈ R for all k ∈ K, and the r k exponents are assumed ordered so that r k > 0 for all k ≤ m and r k < 0 for m ≤ k ≤ K. Because the exponents of the variables are real constants, this is a generalization of a multivariate monomial term. A signomial function is a sum of signomial terms. In [19], a set of transformations of the form x k = f k (z k ) are proposed, where x k is a problem variable, z k is a variable in the reformulated problem and f k is suitable function that can be either exponential or power. This yields an opt-reformulation where all the inequality constraints are convex, and the variables z and the associated (inverse) defining constraints x k = f k (z k ) are added to the reformulation for all k ∈ K (over each signomial term of each signomial constraint). We distinguish the following cases:

• If a > 0, the transformation functions f k are exponential univariate, i.e. x k = e z k . This reformulates [START_REF] Duran | An outer-approximation algorithm for a class of mixed-integer nonlinear programs[END_REF] as follows:

a e P k≤m r k z k Q K k=m+1 x |r k | k ∀k ≤ K x k = e z k .
• If a < 0, the transformation functions are power univariate, i.e.

x k = z 1 R k for k ≤ m and x k = z -1 R k for k > m, where R = k≤K |r k |.
This is also called a potential transformation. This reformulates [START_REF] Duran | An outer-approximation algorithm for a class of mixed-integer nonlinear programs[END_REF] as follows:

a k≤K z |r k | R k ∀k ≤ m x k = z 1 R k ∀k > m x k = z -1 R k R = k≤K |r k |.           
This opt-reformulation isolates all nonconvexities in the inverse defining constraints. These are transformed as follows:

∀k ≤ K x k = e z k → ∀k ≤ K z k = log x k ∀k ≤ m z k = x R k ∀k > m z k = x -R k ,
and then relaxed using a piecewise linear approximation as per Fig. 4. This requires the introduction of binary variables (one per turning point). The signomial relaxation is a convex MINLP; it can be further relaxed to a MILP by outer approximation of the convex terms, or to a convex NLP by continuous relaxation of the discrete variables. This relaxation is denoted by SignomialRelax(P ).

Valid cuts

Once a relaxation has been derived, it should be strengthened (i.e. it should be modified so that the deriving bound becomes tighter). This is usually done by tightening the relaxation, i.e. by adding inequalities. These inequalities have the property that they are redundant with respect to the original (or reformulated) problem but they are not redundant with respect to the relaxation. Thus, they tighten the relaxation but do not change the original problem. In this section we discuss such inequalities for MILPs, NLPs and MINLPs.

′ obtained by Q from adding c to C(Q) has F (P ) ⊆ F(Q ′ ).
Naturally, because Q can be seen as a constraint relaxation of Q ′ , we also have F (Q ′ ) ⊆ F(Q). Linear valid inequalities are very important as adding a linear inequality to an optimization problem usually does not significantly alter the solution time.

For any problem P and any c ∈ C(P ), let F c be the set of points in R n that satisfy c. Let Q be a relaxation of P . Definition 16. A linear valid inequality c is a valid cut if there exists y ∈ Q such that y ∈ F c .

Valid cuts are linear valid inequalities that "cut away" a part of the feasible region of the relaxation. They are used in two types of algorithms: cutting plane algorithms and Branch-and-Bound algorithms. The typical iteration of a cutting plane algorithm solves a problem relaxation Q (say with solution x ′ ), derives a valid cut that cuts away x ′ ; the cut is then added to the relaxation and the iteration is repeated. Convergence is attained when x ′ ∈ F(P ). Cutting plane algorithms were proposed for MILPs [START_REF] Gomory | Essentials of an algorithm for integer solutions to linear programs[END_REF] but then deemed to be too slow for practical purposes, and replaced by Branch-and-Bound. Cutting plane algorithms were also proposed for convex [START_REF] Kelley | The cutting plane method for solving convex programs[END_REF] and bilinear [START_REF] Konno | A cutting plane algorithm for solving bilinear programs[END_REF] NLPs, and pseudoconvex MINLPs [132,131].

Valid cuts for MILPs

This is possibly the area of integer programming where the highest number of papers is published annually. It would be outside the scope of this chapter to relate on all valid cuts for MILPs, so we limit this section to a brief summary. The most effective cutting techniques usually rely on problem structure. See [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF], Ch. II.2 for a good technical discussion on the most standard techniques, and [START_REF] Margot | Pruning by isomorphism in branch-and-cut[END_REF][START_REF] Margot | Exploiting orbits in symmetric ILP[END_REF][START_REF] Kaibel | Packing and partitioning orbitopes[END_REF] for recent interesting group-theoretical approaches which are applicable to large subclasses of IPs. Valid inequalities are generated by all relaxation hierarchies (like e.g. Chvátal- Gomory [133] or ). The best known general-purpose valid cuts are the Gomory cuts [START_REF] Gomory | Essentials of an algorithm for integer solutions to linear programs[END_REF], for they are simple to define and can be written in a form suitable for straightforward insertion in a simplex tableau; many strengthenings of Gomory cuts have been proposed (see e.g. [START_REF] Letchford | Strengthening Chvátal-Gomory cuts and Gomory fractional cuts[END_REF]). Lift-and-project techniques are used to generate new cuts from existing inequalities [START_REF] Balas | Projection, lifting and extended formulation in integer and combinatorial optimization[END_REF]. Families of valid cuts for general Binary Integer Programming (BIP) problems have been derived, for example, in [START_REF] Balas | Canonical cuts on the unit hypercube[END_REF][START_REF] Maculan | Geometrical cuts for 0-1 integer programming[END_REF], based on geometrical properties of the definition hypercube {0, 1} n . In [START_REF] Balas | Canonical cuts on the unit hypercube[END_REF], inequalities defining the various faces of the unit hypercube are derived. The cuts proposed in [START_REF] Maculan | Geometrical cuts for 0-1 integer programming[END_REF] are defined by finding a suitable hyperplane separating a unit hypercube vertex x from its adjacent vertices. Intersection cuts [START_REF] Balas | Intersection cuts -a new type of cutting planes for integer programming[END_REF] are defined as the hyperplane passing through the intersection points between the smallest hypersphere containing the unit hypercube and n halflines of a cone rooted at the current relaxed solution of Q. Spherical cuts are similar to intersection cuts, but the considered sphere is centered at the current relaxed solution, with radius equal to the distance to the nearest integral point [START_REF] Liberti | Spherical cuts for integer programming problems[END_REF]. In [START_REF] Boyd | Fenchel cutting planes for integer programs[END_REF], Fenchel duality arguments are used to find the maximum distance between the solution of Q and the convex hull of the F (P ); this gives rise to provably deep cuts called Fenchel cuts. See [START_REF] Cornuéjols | Valid inequalities for mixed integer linear programs[END_REF] for a survey touching on the most important general-purpose MILP cuts, including Gomory cuts, Lift-and-project techniques, Mixed Integer Rounding (MIR) cuts, Intersection cuts and Reduce-and-split cuts.

Valid cuts for NLPs

Valid cuts for NLPs with a single objective function f subject to linear constraints are described in [START_REF] Horst | Global Optimization: Deterministic Approaches[END_REF] (Ch. III) when an incumbent x * with f (x * ) = γ is known, in order to cut away feasible points x ′ with f (x ′ ) > γ. Such cuts are called γ-valid cuts. Given a nondegenerate vertex x ′ of the feasible polyhedron for which f (x ′ ) > γ, we consider the n polyhedron edges emanating from x ′ . For each i ≤ n we consider a point x i on the i-th edge from x ′ such that f (x i ) ≥ γ. The hyperplane passing through the intersection of the x i is a γ-valid cut (see Fig. 5). More precisely, let Q be the matrix whose i-th column is x ix ′ and e the unit n-vector. Then by [START_REF] Horst | Global Optimization: Deterministic Approaches[END_REF] Thm. III.1 eQ -1 (xx ′ ) ≥ 1 defines a γ-valid cut. Under some conditions, we can find x i such that f (x) = x i and define the strongest possible γ-valid cut, also called concavity cut.

The idea for defining γ-valid cuts was first proposed in [128]; this was applied to 0-1 linear programs by means of a simple reformulation in [100]. It is likely that this work influenced the inception of intersection cuts [START_REF] Balas | Intersection cuts -a new type of cutting planes for integer programming[END_REF] (see Sect. 4.4), which was then used as the basis for current work on Reduce-and-Split cuts [START_REF] Andersen | Reduce-and-split cuts: Improving the performance of mixedinteger Gomory cuts[END_REF]. Some valid cuts for pseudoconvex optimization problems are proposed in [132]. An optimization problem is pseudoconvex if the objective function is a linear form and the constraints are in the form c = (g, -1, 0) where g(x) is a pseudoconvex function of the problem variable vector x. A function g : S ⊆ R n → R is pseudoconvex if for all x, y ∈ S, g(x) < g(y) implies ∇g(y)(xy) < 0. So it follows that for each x, y ∈ S with g(y) > 0, there is a constant α ≥ 1 such that

x ′ x 1 x 2 f (x) = γ
g(y) + α(∇g(y))(x -y) ≤ g(x) (32)
is a (linear) outer approximation to the feasible region of the problem. If g is convex, α = 1 suffices.

In [START_REF] Nowak | Relaxation and Decomposition Methods for Mixed Integer Nonlinear Programming[END_REF], Ch. 7 presents a non-exhaustive list of NLP cuts, applicable to a MINLP standard form ( [START_REF] Nowak | Relaxation and Decomposition Methods for Mixed Integer Nonlinear Programming[END_REF] Eq. ( 7.1): minimization of a linear objective subject to linear inequality constraints and nonlinear inequality constraints): linearization cuts (outer approximation, see Sect. 4.2), knapsack cuts (used for improving loose convex relaxations of given constraints), interval-gradient cuts (a linearization carried out on an interval where the gradient of a given constraint is defined), Lagrangian cuts (derived by solving Lagrangian subproblems), level cuts (defined for a given objective function upper bound), deeper cuts (used to tighten loose Lagrangian relaxation; they involve the solution of separation problems involving several variable blocks).

Another NLP cut based on the Lagrangian relaxation is proposed in [124]: consider a MINLP in the canonical form min g(x)≤0 f (x) and let L(•, µ) = f (x)+µ ⊤ g(x) be its Lagrangian relaxation. Let f be a lower bound obtained by solving L and f be an upper bound computed by evaluating f at a feasible point x ′ . From f ≤ f (x)+µ ⊤ g(x) ≤ f +µ ⊤ g(x) one derives the valid cut g i (x) ≥ -1 µi ( f -f ) for all i ≤ m (where g : R n → R m ).

Valid cuts for MINLPs

Naturally, both MILP and NLP cuts may apply to MINLPs. Some more specific MINLP cuts can be derived by reformulating integer variables to binary (see Sect. 3.2) and successively to continuous (see Sect. 3.2). The added quadratic constraints may then be relaxed in a Lagrangian (see Sect. 4.3) or SDP fashion (see Sect. 4.3) [START_REF] Plateau | Reformulations quadratiques convexes pour la programmation quadratique en variables 0-1[END_REF]: any of the NLP cuts described in Sect. 4.4 applied to such a reformulation is essentially a specific MINLP valid cut.

Reformulation/Optimization Software Engine

Although specific reformulations are carried out by most LP/MILP preprocessors [START_REF]ILOG: ILOG CPLEX 11.0 User's Manual[END_REF][START_REF] Guéret | Applications of optimization with Xpress-MP. Dash Optimization[END_REF], and a few very simple reformulations are carried out by some mathematical programming language environments [START_REF] Fourer | The AMPL Book[END_REF][START_REF] Brook | GAMS, a user's guide[END_REF], there is no software optimization framework that is able to carry out reformulations in a systematic way. In this section we describe the Reformulation/Optimization Software Engine (ROSE), a C++ software framework for optimization that can reformulate and solve mathematical programs of various types. ROSE is work in progress; currently, it is more focused on reformulation than optimization, but it has nonetheless a few native solvers (e.g. a Variable Neighbourhood Search (VNS) based algorithm for nonconvex NLPs [START_REF] Liberti | Variable neighbourhood search for the global optimization of constrained NLPs[END_REF]) and wrappers to various other external solvers (e.g. the LP solver GLPK [START_REF] Makhorin | GNU Linear Programming Kit[END_REF] and the local NLP solver SNOPT [START_REF] Gill | User's Guide for SNOPT 5.3[END_REF]). In our research, we currently use ROSE's reformulation capabilities with AMPL's considerable set of numerical solvers in order to obtain solutions of complex MINLPs. ROSE consists of a set of interfaces with external clients (currently, it has a direct user interface and an AMPL [START_REF] Fourer | The AMPL Book[END_REF] interface), a problem class, a virtual solver class with many implementations, and an expression tree manipulation library called Ev3 (see Sect. 5.3). Reformulations may occur within the problem class, within the solvers, or within Ev3. Solvers may embed either a numerical solution algorithm or a symbolic reformulation algorithm, or both. The problem class builds the problem and simplifies it as much as possible; solvers are generally passed one or more problem together with a set of configuration parameters, and provide either a numerical solution or a reformulation. Reformulation solvers usually change the structure of their input problems; there is a special dedicated reformulation solver that makes an identical copy of the input problem. Most reformulation solvers acting on the mathematical expressions call specific methods within Ev3.

Development history

ROSE has a long history. Its "founding father" is the GLOP software ([71] Sect. 3.3), conceived and used by E. Smith to write his Ph.D. thesis [112] at CPSE, Imperial College, under the supervision of Prof. Pantelides. GLOP was never publically released, although test versions were used by CPSE students and faculty over a number of years. GLOP, however, was not so much a software framework rather than an implementation of the reformulation-based sBB algorithm described in [114]. The same algorithm (in a completely new implementation) as well as some other global optimization algorithms were put together in the ooOPS (object-oriented OPtimization System) software framework ([71] Sect. 3.4), coded by the first author of this chapter during his Ph.D. thesis [START_REF] Liberti | Reformulation and convex relaxation techniques for global optimization[END_REF] at CPSE, Imperial College, and drawing a few software architecture ideas from its MILP predecessor, ooMILP [127,126]. The ooOPS software framework [START_REF] Liberti | ooOPS[END_REF] includes an sBB algorithm for MINLPs (which has a few glitches but works in a lot of instances), a VNS algorithm for nonconvex NLPs, a wrapper to the GO solver SobolOpt [START_REF] Kucherenko | Application of deterministic low-discrepancy sequences in global optimization[END_REF], and a wrapper to SNOPT. ooOPS was used to compile the results of several research papers, but unfortunately Imperial College never granted the rights to distribute its source publically. Besides, ooOPS used binary expression trees, which made it much more difficult to reformulate sums and products with more than two operands. The MINLP Object-oriented Reformulation/Optimization Navigator (MORON) was conceived to address these two limitations. MORON has an extensive API for dealing with both reformulation and optimization and includes: a prototypical Smith reformulator and convexifier ( [START_REF] Liberti | Writing global optimization software[END_REF], Sect. 2.3 and 7.3); a preliminary version of the sBB algorithm; a wrapper to SNOPT. A lot of work was put into the development of Ev3, a separate expression tree library with reformulating capabilities [START_REF] Liberti | Framework for symbolic computation in C++ using n-ary trees[END_REF]. Unfortunately, due to lack of time, development of MORON was discontinued. ROSE is MORON's direct descendant: it has a leaner API, almost the same software architecture (the main classes being Problem and Solver), and it uses Ev3 to handle expression trees. We expect to be able to publically distribute ROSE within the end of 2008; for using and/or contributing to its development, please contact the first author. We also remark that many of the ideas on which ooOPS's and MORON's sBB solvers are built are also found in Couenne [START_REF] Belotti | Branching and bound reduction techniques for non-convex MINLP. Optimization Methods and Software[END_REF], a modern sBB implementation coded by P. Belotti within a CMU/IBM project, targeted at general MINLPs, and publically distributed within COIN-OR [START_REF] Lougee-Heimer | The common optimization interface for operations research: Promoting opensource software in the operations research community[END_REF].

Software architecture

The ROSE software relies on two main classes: Problem and Solver. The standard expected usage sequence is the following. The client (either the user or AMPL) constructs and configures a problem, selects and configures a solver, then solves a problem using the selected solver, and finally collects the output from the problem.

The Problem class has methods for reading in a problem, access/modify the problem description, perform various reformulations to do with adding/deleting variables and constraints, evaluate the problem expressions and their first and second derivatives at a given point, and test for feasibility of a given point in the problem. The Solver class is a virtual class that serves as interface for various solvers. Implementations of the solver class are passed a pointer to a Problem object and a set of user-defined configuration parameters. Solver implementations may either find numerical solutions and/or change the problem structure. Numerical solvers normally re-insert the numerical solution found within the Problem object. The output of a reformulation solver is simply the change carried out on the problem structure. Every action carried out on a mathematical expression, be it a function evaluation or a symbolic transformation, is delegated to the Ev3 library (see Sect. 5.3).

The Problem class ROSE represents optimization problems in their flat form representation; i.e. variables, objective functions and constraints are arranged in simple linear lists rather than in jagged arrays of various dimensions. The reason for this choice is that languages such as AMPL and GAMS already do an excellent job of translating structured form problem formulations to their flat counterparts. Problems are defined in problem.h and problem.cxx.

This class rests on three structs defining variables, objectives and constraints.

• struct Variable, storing the following information concerning decision variables.

-ID, an integer (int) storing an ID associated to the variable. This ID does not change across reformulations, except in case of reformulations which delete variables. In this case, when a variable is deleted the IDs of the successive variables are shifted. The lists storing variable objects do not make any guarantee on the ordering of IDs across the list. -Name, a string (std::string) storing the variable name. This is only used for printing purposes. -LB, a floating point number (double) storing the variable lower bound.

-UB, a floating point number (double) storing the variable upper bound.

-IsIntegral, a flag (bool) set to 1 if the variable is integer and 0 otherwise. Binary variables occur when IsIntegral is set to 1, LB to 0, and UB to 1. -Persistent, a flag (bool) set to 1 if the variable cannot be deleted by reformulation algorithms, and 0 otherwise. -Optimum, a floating point number (double) storing a value for the variable. Notwithstanding the name, this is not always the optimum value. • struct Objective, storing the following information concerning objective functions.

-ID, an integer (int) storing an ID associated to the objective. This ID does not change across reformulations, except in case of reformulations which delete objectives. In this case, when an objective is deleted the IDs of the successive objectives are shifted. The lists storing objective objects do not make any guarantee on the ordering of IDs across the list. -Name, a string(std::string) storing the objective name. This is not currently used.

-Function, the expression tree (Expression) of the objective function.

-FunctionFET, a fast evaluation expression tree (see Sect. 5.3 on p. 46) pointer (FastEval-Tree*) corresponding to Function. -NonlinearPart, the expression tree (Expression) of the nonlinear part of Function. This may contain copies of subtrees of Function. The nonlinear part of an expression includes all subexpressions involving variables that appear nonlinearly at least once in the expression. For example, the nonlinear part of x + y + z + yz is y + z + yz. -NonlinearPartFET, a fast evaluation expression tree pointer (FastEvalTree*) corresponding to NonlinearPart. -OptDir, a label (int) which is 0 if the objective is to be minimized and 1 if it is to be maximized. -Diff, the first-order partial derivatives (std::vector<Expression>) of Function.

-DiffFET, the fast evaluation tree pointers (std::vector<FastEvalTree*>) corresponding to the first-order partial derivatives. -Diff2, the second-order partial derivatives (std::vector<std::vector<Expression> >) of Function. -Diff2FET, the fast evaluation tree pointers (std::vector<std::vector<FastEvalTree*> >) corresponding to second-order partial derivatives. • struct Constraint, storing the following information concerning constraints.

-ID, an integer (int) storing an ID associated to the constraint. This ID does not change across reformulations, except in case of reformulations which delete constraints. In this case, when a constraint is deleted the IDs of the successive constraints are shifted. The lists storing constraint objects do not make any guarantee on the ordering of IDs across the list. -Name, a string(std::string) storing the constraint name. This is not currently used.

-Function, the expression tree (Expression) of the constraint function.

-FunctionFET, a fast evaluation expression tree (see Sect. 5.3 on p. 46) pointer (FastEval-Tree*) corresponding to Function. -NonlinearPart, the expression tree (Expression) of the nonlinear part of Function. This may contain copies of subtrees of Function. The nonlinear part of an expression includes all subexpressions involving variables that appear nonlinearly at least once in the expression. For example, the nonlinear part of x + y + z + yz is y + z + yz. -NonlinearPartFET, a fast evaluation expression tree pointer (FastEvalTree*) corresponding to NonlinearPart. -LB, a floating point number (double) storing the constraint lower bound.

-UB, a floating point number (double) storing the constraint upper bound.

-Diff, the first-order partial derivatives (std::vector<Expression>) of Function.

-DiffFET, the fast evaluation tree pointers (std::vector<FastEvalTree*>) corresponding to the first-order partial derivatives. -Diff2, the second-order partial derivatives (std::vector<std::vector<Expression> >) of Function. -Diff2FET, the fast evaluation tree pointers (std::vector<std::vector<FastEvalTree*> >) corresponding to second-order partial derivatives. We remark that Constraint objects are expressed in the form LB ≤ Function ≤ UB; in order to deactivate one constraint side, use the defined constant MORONINFINITY (1 × 10 30 ).

Indexing of problem entities

Pointers to all variable, objective and constraint objects (also called entities) in the problem are stored in STL vectors. Thus, on top of the entity indexing given by the ID property, we also have the natural indexing associated to these vectors, referred to as local indexing. Whereas ID-based indices are constant throughout any sequence of reformulations, local indices refer to the current problem structure. Direct and inverse mappings between indices and local indices are given by the following Problem methods:

• int GetVariableID(int localindex) • int GetVarLocalIndex(int varID) • int GetObjectiveID(int localindex) • int GetObjLocalIndex(int objID) • int GetConstraintID(int localindex) • int GetConstrLocalIndex(int constrID).
Individual problem entities can be accessed/modified by their ID; a subset of the methods also exist in the "local index" version -such methods have the suffix -LI appended to their names. All indices in the API start from 1.

Parameters

The parameter passing mechanism is based on a Parameters class with the following methods.

Problem API

The API of the Problem class is given in the tables on pages 43-44. Within the reformulation methods, the Add-methods automatically call a corresponding New-method to produce the next available ID. The DeleteVariable methods does not eliminate all occurrences of the variable from the problem (i.e. this is not a projection). The AMPL-based construction methods were made possible by an undocumented AMPL solver library feature that allows clients to access AMPL's internal binary trees [38,[START_REF] Galli | Parsing AMPL internal format for linear and non-linear expressions[END_REF].

The Solver virtual class

Solver is a virtual class whose default implementation is an inactive (empty) solver. This is not a pure virtual class because it represent the union of all possible solver implementations, rather than the intersection; in other words, not all methods in Solver are implemented across all solvers (check the source files solver*.h, solver*.cxx to make sure).

Existing Solver implementations

Each solver implementation consists of a header and an implementation file. Currently, ROSE has three functional numerical solvers: VNS solver for nonconvex NLPs [START_REF] Liberti | Variable neighbourhood search for the global optimization of constrained NLPs[END_REF], a wrapper to SNOPT [START_REF] Gill | User's guide for SNOPT version 7[END_REF], a wrapper to GLPK [START_REF] Makhorin | GNU Linear Programming Kit[END_REF]; and various reformulator solvers, among which: a problem analyser that returns problem information to AMPL, a problem copier that simply makes an identical copy of the current problem (for later reformulations), an outer approximation reformulator, a

Ev3

Ev3 is a library providing expression tree functionality and symbolic transformations thereof (see Sect. 2.2). This library may also be used stand-alone, and the rest of this section actually refers to the stand-alone version. The only adaptation that was implemented for usage within ROSE was to provide additional structures for Fast Evaluation Trees (FETs). Ev3's native trees are very easy to change for reformulation needs, but unfortunately turn out to be slow to evaluate by Alg. 1. Since in most numerical algorithms for optimization the same expressions are evaluated many times, a specific data structure fevaltree with relative source files (fastexpression.h, fastexpression.cxx) have been added to Ev3. FETs are C-like n-ary (as opposed to binary) trees that have none of the reformulating facilities of their Ev3 counterparts, but which are very fast to evaluate. Construction and evaluation of FETs is automatic and transparent to the user.

Architecture

The Ev3 software architecture is mainly based on 5 classes. Two of them, Tree and Pointer, are generic templates that provide the basic tree structure and a no-frills garbage collection based on reference count. Each object has a reference counter which increases every time a reference of that object is taken; the object destructor decreases the counter while it is positive, only actually deleting the object when the counter reaches zero. This type of garbage collecting is due to Collins, 1960 (see [START_REF] Kaltofen | Challenges of symbolic computation: My favorite open problems[END_REF]). Other two classes, Operand and BasicExpression, implement the actual semantics of an algebraic expression. The last class, ExpressionParser, implements a simple parser (based on the ideas given in [117]) which reads in a string containing a valid mathematical expression and produces the corresponding n-ary tree. The constructor of this class allocates a new integer for the reference counter ncount and a new NodeType object, and the copy constructor increases the counter. The destructor deletes the reference counter and invokes the delete method on the NodeType object. In order to access the data and methods of the NodeType object pointed to by node, the -> operator in the Pointer class is overloaded to return node.

A mathematical expression, in Ev3, is defined as a pointer to a BasicExpression object (see below for the definition of a BasicExpression object):

typedef Pointer<BasicExpression> Expression;
The Tree class This is a template class defined as template<class NodeType> class Tree { vector<Pointer<NodeType> > nodes; // methods }; This is the class implementing the n-ary tree (subnodes are contained in the nodes vector). Notice that, being a template, the whole implementation is kept independent of the semantics of a NodeType. Notice also that because pointers to objects are pushed on the vector, algebraic substitution is very easy: just replace one pointer with another one. This differs from the implementation of GiNaC [START_REF] Bauer | Introduction to the ginac framework for symbolic computation within the C++ programming language[END_REF] where it appears that algebraic substitution is a more convoluted operation.

The Operand class

This class holds the information relative to each expression term, be they constants, variables or operators. This makes it easy to perform symbolic manipulation on like terms (e.g. x + 2x = 3x). • every leaf term (variables and constants) can be raised to a numeric power. This makes it easy to perform symbolic manipulation of polynomials.

Introducing numeric coefficients and exponents is a choice that has advantages as well as disadvantages. GiNaC, for example, does not explicitely account for numeric coefficients. The advantages are obvious: it makes symbolic manipulation very efficient for certain classes of basic operations (operations on like terms). The disadvantage is that the programmer has to explicitely account for the case where terms are assigned coefficients: whereas with a pure tree structure recursive algorithms can be formulated as "for each node, do something", this becomes more complex when numeric coefficients are introduced. Checks for non-zero or non-identity have to be performed prior to carrying out certain operations, as well as having to manually account for cases where coefficients have to be used. However, by setting both multiplicative and exponent coefficients to 1, the mechanism can to a certain extent be ignored and a pure tree structure can be recovered.

The BasicExpression class

This class is defined as follows:

class BasicExpression : public Operand, public Tree<BasicExpression> { // methods };

It includes no data of its own, but it inherits its semantic data from class Operand and its tree structure from template class Tree with itself (BasicExpression) as a base type. This gives BasicExpression an n-ary tree structure. Note that an object of class BasicExpression is not a Pointer, only its subnodes (if any) are stored as Pointers to other BasicExpressions. This is the reason why the client code should never explicitely use BasicExpression; instead, it should use objects Expression, which are defined as Pointer<BasicExpression>. This allows the automatic garbage collector embedded in Pointer to work.

The ExpressionParser class

This parser originates from the example parser found in [117]. The original code has been extensively modified to support exponentiation, unary functions in the form f (x), and creation of n-ary trees of type Expression. For an example of usage, see Section 5.3 below.

Application Programming Interface

The Ev3 API consists in a number of internal methods (i.e., methods belonging to classes) and external methods (functions whose declaration is outside the classes). Objects of type class Expression can be built from strings containing infix-format expressions (like, e.g. "log(2*x*y)+ sin(z)") by using the built-in parser. However, they may also be built from scratch using the supplied construction methods (see Section 5.3 for examples). Since the fundamental type Expression is an alias for Pointer<BasicExpression>, and BasicExpression is in turn a mix of different classes (including a Tree with itself as a template type), calling internal methods of an Expression object may be confusing. Thus, for each class name involved in the definition of Expression, we have listed the calling procedure explicitly in the tables on pages [START_REF] Horst | On the convexification of nonlinear programming problems: an applications-oriented approach[END_REF] -they combine references of their arguments; -they may change their arguments; -they may change the order of the subnodes where the operations are commutative; -they may return one of the arguments. Thus, it is advisable to perform the operations on copies of the arguments when the expression being built is required to be independent of its subnodes. In particular, all the expression building functions (e.g. operator+(), . . . , Log(), . . . ) do not change their arguments, whereas their -Link counterparts do. • The built-in parser (ExpressionParser) uses linking and not copying (also see Section 5.3) of nodes when building up the expression. method. It is important that variable indices should start from 1 and increase monotonically by 1, as variable indices are used to index the array of values passed to the Eval() method.

Copying vs. Linking

One thing that is immediately noticeable is that this architecture gives a very fine-grained control over the construction of expressions. Subnodes can be copied or "linked" (i.e., a reference to the object is put in place, instead of a copy of the object -this automatically uses the garbage collection mechanism, so the client code does not need to worry about these details). Copying an expression tree entails a set of advantages/disadvantages compared to linking. When an expression is constructed by means of a copy to some other existing expression tree, the two expressions are thereafter completely independent. Manipulation one expression does not change the other. This is the required behaviour in many cases. The symbolic differentiation routine has been designed using copies because a derivative, in general, exists independently of its integral. Linking, however, allows for facilities such as "propagated simplification", where some symbolic manipulation on an expression changes all the expressions having the manipulated expression tree as a subnode. This may be useful but calls for extra care. The built-in parser has been designed using linking because the "building blocks" of a parsed expression (i.e. its subnodes of all ranks) will not be used independently outside the parser.

Simplification Strategy

The routine for simplifying an expression repeatedly calls a set of simplification rules acting on the expression. These rules are applied to the expression as long as at least one of them manages to further simplify it.

Simplifications can be horizontal, meaning that they are carried out on the same list of subnodes (like e.g. x + y + y = x + 2y), or vertical, meaning that the simplification involves changing of node level (like e.g. application of associativity: ((x + y) + z) = (x + y + z)).

The order of the simplification rules applied to an object Expression e is the following:

1. e->ConsolidateProductCoeffs(): in a product having n subnodes, collect all multiplicative coefficients, multiply them together, and set the result as the multiplicative coefficient of the whole product:

n i=1 (c i f i ) = ( n i=1 c i )( n i=1 f i ).
2. e->DistributeCoeffOverSum(): in a sum with n subnodes and a non-unit multiplicative coefficient, distribute this coefficient over all subnodes in the sum: 7. SimplifyRecursive(e): deals with the most common simplification rules, i.e.:

c n i=1 f i = n i=1 cf i . 3
• try to simplify like terms in fractions where numerator and denominator are both products;

• x ± 0 = 0 + x = x; • x × 1 = 1 × x = x; • x × 0 = 0 × x = 0; • x 0 = 1; • x 1 = x; • 0 x = 0; • 1 x = 1.

Differentiation

Derivative rules are the usual ones; the rule for multiplication is expressed in a way that allows for n-ary trees to be derived correctly:

∂ ∂x n i=1 f i = n i=1   ∂f i ∂x j =i f j   .

Algorithms on n-ary Trees

We store mathematical expressions in a tree structure so that we can apply recursive algorithms to them. Most of these algorithms are based on the following model.

if expression is a leaf node do something else recurse on all subnodes do something else end if

In particular, when using Ev3, the most common methods used in the design of recursive algorithms are the following:

• IsLeaf(): is the node a leaf node (variable or constant)?

• GetSize(): find the number of subnodes of any given node.

• GetOpType(): return the type of operator node.

• GetNode(int i): return the i-th subnode of this node (nodes are numbered starting from 0). • DeleteNode(int i): delete the i-th subnode of this node (care must be taken to deal with cases where all the subnodes have been deleted -Ev3 allows the creation of operators with 0 subnodes, although this is very likely to lead to subsequent errors, as it has no mathematical meaning). • Use of the operators for manipulation of nodes: supposing Expression e, f contain valid mathematical expressions, the following are all valid expressions (the new expressions are created using copies of the old ones).

Expression e1 = e + f; Expression e2 = e * Log(Sqrt(e^2 -f^2)); Expression e3 = e + f -f; // this is automatically by simplified to e

Ev3 usage example

The example in this section explains the usage of the methods which represent the core, high-level functionality of Ev3: fast evaluation, symbolic simplification and differentiation of mathematical expressions.

The following C++ code is a simple driver program that uses the Ev3 library. Its instructions should be self-explanatory. First, we create a "parser object" of type ExpressionParser. We then set the mapping variable names / variable indices, and we parse a string containing the mathematical expression log(2xy) + sin(z). We print the expression, evaluate it at the point (2, 3, 1), and finally calculate its symbolic derivatives w.r.t. x, y, z, and print them.

#include "expression.h" #include "parser.h" int main(int argc, char** argv) { ExpressionParser p; // create the parser object p.SetVariableID("x", 1) // map between symbols and variable indices p.SetVariableID("y", 2) // x --> 0, y --> 1, z --> 2 p.SetVariableID("z", 3) int parsererrors = 0; // number of parser errors /* call the parser's Parse method, which returns an Expression which is then used to initialize Expression e */ Expression e(p.Parse("log(2*x*y)+sin(z)", parsererrors)); cout << "parsing errors: " << parsererrors << endl; cout << "f = " << e->ToString() << endl; // print the expression double val[3] = {2, 3, 1}; cout << "eval(2,3,1):

" << e->Eval(val, 3) << endl; // evaluate the expr. cout << "numeric check: " << ::log(2*2*3)+::sin(1) << endl; // check result // test diff Expression de1 = Diff(e, 1); // calculate derivative w.r.t. x cout << "df/dx = " << de1->ToString() << endl; // print derivative Expression de2 = Diff(e, 2); // calculate derivative w.r.t. y cout << "df/dy = " << de2->ToString() << endl; // print derivative Expression de3 = Diff(e, 3); // calculate derivative w.r.t. z cout << "df/dz = " << de3->ToString() << endl; // print derivative return 0; } The corresponding output is parsing errors: 0 f = (log((2*x)*(y)))+(sin(z)) eval(2,3,1):

3.32638 numeric check: 3.32638 df/dx = (1)/(x) df/dy = (1)/(y) df/dz = cos(z)

Notes

• In order to evaluate a mathematical expression f (x 1 , x 2 , . . . , x n ), where x i are the variables and i are the variable indices (starting from 1 and increasing by 1), we use the Eval() internal method, whose complete declaration is as follows:

double Expression::Eval(double* varvalues, int size) const;

The array of doubles varvalues contains size real constants, where size >= n. The variable indices are used to address this array (the value assigned to x i during the evaluation is varvalues[i-1]), so it is important that the order of the constants in varvalues reflects the order of the variables. This method does not change the expression object being evaluated. • The core simplification method is an external method with declaration bool Simplify(Expression* e); It consists of a number of different simplifications, as explained in Section 5.3. It takes a pointer to Expression as an argument, and it returns true if some simplification has taken place, and false otherwise. This method changes its input argument.

• The symbolic differentiation procedure is an external method:

Expression Diff(const Expression& e, int varindex); It returns a simplified expression which is the derivative of the expression in the argument with respect to variable varindex. This method does not change its input arguments.

• External class methods take Expressions as their arguments. According as to whether they need to change their input argument or not, the Expression is passed by value, by reference, or as a pointer. This may be a little confusing at first, especially when using the overloaded -> operator on Expression objects. Call a method in the Pointer<> class. In particular, care must be taken between the two forms e->MethodName() and ePtr->Method-Name() as they are syntactically very similar but semantically very different.

Validation examples

As validation examples, we show ROSE's output on simple input problems by using two kind of reformulations. In order to ease the reading of the examples, we use an intuitive description format for MINLPs problems [71, pages 237-239]. It is worth noticing that the symbol '<' stands here for '≤' and that we use an explicit boundary (1e 30 ) for dealing with infinity.

The first example performs the reformulation of products between continuous and binary variables.

Original Problem As presented in Section 3.3, ROSE identifies all the terms involving a continuous and a binary variable (respectively x1 and x2 in the example) and add exactly one variable (w3 here) and four constraints. The reader might now check that both the objective function and the constraints are linear terms and that the computed values are similar in the two formulations of the problem.

The second example is an optimization problem whose objective function contains four nonlinear terms. We show how ROSE is able to find a convex relaxation for the problem using the convexifier reformulator (see Section 4.2). The reformulation process is performed in various steps. In order to explain how the reformulator/convexifier works, we show in the following how the original problem is modified during the main steps.

Original Problem

The first step consists in reformulating the problem to the Smith standard form. Each nonconvex term in the objective function is replaced by an added variable w and defining constraints of the form w = nonconvex term are added to the problem. The objective function of the reformulated problem is one linearizing variable only, that is the sum of all the added variables, and a constraint for this equation is also added to the problem. We remark that the obtained reformulation is a lifting reformulation, since a new variable is added for each nonconvex term. This first-stage reformulation is the following: Then, each defining constraint is replaced by a convex under-estimator and concave overestimator of the corresponding nonlinear term. In particular, the term 2*x1^2 is treated as a convex univariate function f(x) and a linear under-estimator is obtained by considering five tangents to f at various given points, an over-estimator is obtained by considering the secant through the points (x1^L,f(x1^L)),(x1^U,f(x1^U)), where x1^L,and x1^U are the bounds on x1.For the term y2^3, where the range of y2 includes zero, the linear relaxation given in [START_REF] Liberti | Convex envelopes of monomials of odd degree[END_REF] Finally, the Smith defining constraints are removed, obtaining the final reformulation (of the relaxation type).

Conclusion

This chapter contains a study of mathematical programming reformulation and relaxation techniques. Section 1 presents some motivations towards such a study, the main being that Mixed Integer Nonlinear Programming solvers need to be endowed with automatic reformulation capabilities before they can be as reliable, functional and efficient as their industrial-strength Mixed Integer Linear Programming solvers are. Section 2 presents a general framework for representing and manipulating mathematical programming formulations, as well as some definitions of the concept of reformulation together with some theoretical results; the section is concluded by listing some of the most common standard forms in mathematical programming. In Section 3 we present a partial systematic study of existing reformulations. Each reformulation is presented both in symbolic algorithmic terms (i.e. a prototype for carrying out the reformulation automatically in terms of the provided data structures is always supplied) and in the more usual mathematical terms. This should be seen as the starting point for a more exhaustive study: eventually, all known useful reformulations might find their place in an automatic reformulation preprocessing software for Mixed Integer Nonlinear Programming. In Section 4, we attempt a similar work with respect to relaxations. Section 5 describes the implementation of ROSE, a reformulation/optimization software engine.

Index

  is the sequence of objective functions; each objective function o ∈ O has the form (d o , f o ) where d o ∈ {-1, 1} is the optimization direction (-1 stands for minimization, +1 for maximization) and f o ∈ E; • C ⊆ E × S × R (where S = {-1, 0, 1}) is the sequence of constraints c of the form (e c , s c , b c ) with e c ∈ E, s c ∈ S, b c ∈ R:

Fig. 1 .

 1 Fig. 1. The graphs e1 (above) and e2 (below) from Example 2.1.

6 i=1

 6 x i = 0.

Fig. 2 .

 2 Fig. 2. The BGBP instance in Example 2.1.

Fig. 3 .

 3 Fig. 3. The graph e ′ 1 from Example 2.1. L ′ ij = Lij + Lji for all i, j.

  +∞}. More precisely, an NLP is in SQP standard form if for all c ∈ C (a) s c = 0 and (b) there is c ′ ∈ C such that e c = e c ′ and s c = -s c ′ .

  where e o is a linear form; • C can be partitioned into two sets of constraints C 1 , C 2 such that c is a linear form for all c ∈ C 1 and c = (e c , 0, 0) for c ∈ C 2 where e c is as follows: 1. r(e c ) is the sum operator 2. δ + (r(e c )) = {⊗, v} where (a) ⊗ is a nonlinear operator where all subnodes are leaf nodes, (b) χ(v) = -1 and (c) τ (v) ∈ V.

Example 2 .

 2 Consider the problem P ≡ min x∈[-2π,2π]

  Consider a problem P having an objective function o = (d o , e o ) ∈ O(P ) and a finite set of expressions e k for k ∈ K such that e o = d o min k∈K d o e k (this is a piecewise linear objective function of the form min max k e k or max min k e k depending on d o ). This can be linearized by adding one variable and |K| constraints to the problem as follows:

  Consider a problem P with a single objective function o = (d o , e o ) ∈ O where e o = (-d o ) k∈ K e k where the operator represented by the root node r(e k ) of e k is the absolute value | • | for all k ∈ K ⊆ K. Since the absolute value operator is unary, δ + (r(e k )) consists of the single element f k . Provided f k are linear forms, this problem can be linearized as follows. For each k ∈ K:

  0.452055) x_1 + (-1.20548) x_2 + (1.50685) = 0 2: (0.769231) x_1 + (1.15385) x_2 + (-8.84615) = 0 Assignment of points to hyperplanar clusters: hyp_cluster 1 = { 2 3 4 8 } hyp_cluster 2 = { 1 5 6 7 }.

3 . 2 )

 32 to each nonlinear equality constraint in C, obtaining an opt-reformulation P 1 of P . 2. For every nonconvex inequality constraint c = (e c , s c , b c ) ∈ C(P 1 ):

Fig. 4 .

 4 Fig. 4. Piecewise linear underestimating approximations for concave (left) and convex (right) univariate functions.

Definition 15 .

 15 Given an optimization problem P and a relaxation Q, a valid inequality is a constraint c = (e c , s c , b c ) such that the problem Q

Fig. 5 .

 5 Fig. 5. A γ-valid cut.

.

  DifferenceToSum(e): replace all differences and unary minus with sums, multiplying the coefficient of the operands by -1. 4. SimplifyConstant(e): simplify operations on constant terms by replacing the value of the node with the result of the operation. 5. CompactProducts(e): associate products; e.g. ((xy)z) = (xyz). 6. CompactLinearPart(e): this is a composite simplification consisting of the following routines: a) CompactLinearPartRecursive(e): recursively search all sums in the expression and perform horizontal and vertical simplifications on the coefficients of like terms. b) ReorderNodes(e): puts each list of subnodes in an expression in standard form (also see Sect. 2.2): constant + monomials in rising degree + complicated operands (where complicated operands are sublists of subnodes).

  constraints = [ -1e+30 < (-30*x2)+(w3) < 0 ], [ -1e+30 < (15*x2)+(-1*w3) < 0 ], [ -1e+30 < (15)+(-1*x1)+(-15*x2)+(w3) < 0 ], [ -1e+30 < (-30)+(x1)+(30*x2)+(-1*w3) < 0 ];

#

  ROSE problem: convexifier # Problem has 3 variables and 1 constraints # Variables: variables = -1 < x1 < 1 / Continuous, -2 < y2 < 3 / Continuous, 1 < t3 < 2 / Continuous; # Objective Function: objfun = min [ (2*x1^2)+(y2^3) +((x1)*(y2))+((x1)/(t3)) ]; constraints = [ 2 < (x1)+(y2) < 1e+30 ]; # Constraints: constraints = [ 2 < (x1)+(y2) < 1e+30 ], [ 0 < (w4)+(w5)+(w6)+(w7)+(-1*w8) < 0 ], [ -2 < (4*x1)+(w4) < 1e+30 ], [ -2 < (-4*x1)+(w4) < 1e+30 ], [ -0.5 < (2*x1)+(w4) < 1e+30 ], [ -0.5 < (-2*x1)+(w4) < 1e+30 ], [ -2 < (-3*y2)+(w5) < 1e+30 ], [ -54 < (-27*y2)+(w5) < 1e+30 ], [ -1e+30 < (-6.75*y2)+(w5) < 6.75 ], [ -1e+30 < (-12*y2)+(w5) < 16 ], [ -2 < (2*x1)+(y2)+(w6) < 1e+30 ], [ -3 < (-3*x1)+(-1*y2)+(w6) < 1e+30 ], [ -1e+30 < (-3*x1)+(y2)+(w6) < 3 ], [ -1e+30 < (2*x1)+(-1*y2)+(w6) < 2 ], [ 0.5 < (-0.5*x1)+(w7)+(z9) < 1e+30 ], [ -1 < (-1*x1)+(w7)+(-1*z9) < 1e+30 ], [ -1e+30 < (-1*x1)+(w7)+(z9) < 1 ], [ -1e+30 < (-0.5*x1)+(w7)+(-1*z9) < -0.5 ];

#

  ROSE problem: convexifier # Problem has 8 variables and 6 constraints # Variables: variables = -1 < x1 < 1 / Continuous, -2 < y2 < 3 / Continuous, 1 < t3 < 2 / Continuous, 0 < w4 < 2 / Continuous, -8 < w5 < 27 / Continuous, -3 < w6 < 3 / Continuous, -1 < w7 < 1 / Continuous, -12 < w8 < 33 / Continuous; # Objective Function: objfun = min [ w8 ]; # Constraints: constraints = [ 2 < (x1)+(y2) < 1e+30 ], [ 0 < (-1*w4)+(2*x1^2) < 0 ], [ 0 < (-1*w5)+(y2^3) < 0 ], [ 0 < (-1*w6)+((x1)*(y2)) < 0 ], [ 0 < (-1*w7)+((x1)/(t3)) < 0 ], [ 0 < (w4)+(w5)+(w6)+(w7)+(-1*w8) < 0 ];

  is used. McCormick's envelopes are considered for the bilinear term x1*y2.The fractional term is reformulated as bilinear by considering z=1/t3 and McCormick's envelopes are exploited again. We obtain the following relaxation: # ROSE problem: convexifier # Problem has 9 variables and 22 constraints # Variables: variables = -1 < x1 < 1 / Continuous, -2 < y2 < 3 / Continuous, 1 < t3 < 2 / Continuous, 0 < w4 < 2 / Continuous, -8 < w5 < 27 / Continuous, -3 < w6 < 3 / Continuous, -1 < w7 < 1 / Continuous, -12 < w8 < 33 / Continuous, 0.5 < z9 < 1 / Continuous; # Objective Function: objfun = min [ w8 ]; # Constraints: constraints = [ 2 < (x1)+(y2) < 1e+30 ], [ 0 < (-1*w4)+(2*x1^2) < 0 ], [ 0 < (-1*w5)+(y2^3) < 0 ], [ 0 < (-1*w6)+((x1)*(y2)) < 0 ], [ 0 < (-1*w7)+((x1)/(t3)) < 0 ], [ 0 < (w4)+(w5)+(w6)+(w7)+(-1*w8) < 0 ], [ -2 < (4*x1)+(w4) < 1e+30 ], [ -2 < (-4*x1)+(w4) < 1e+30 ], [ -0.5 < (2*x1)+(w4) < 1e+30 ], [ -0.5 < (-2*x1)+(w4) < 1e+30 ], [ -2 < (-3*y2)+(w5) < 1e+30 ], [ -54 < (-27*y2)+(w55) < 1e+30 ], [ -1e+30 < (-6.75*y2)+(w5) < 6.75 ], [ -1e+30 < (-12*y2)+(w5) < 16 ], [ -2 < (2*x1)+(y2)+(w6) < 1e+30 ], [ -3 < (-3*x1)+(-1*y2)+(w6) < 1e+30 ], [ -1e+30 < (-3*x1)+(y2)+(w6) < 3 ], [ -1e+30 < (2*x1)+(-1*y2)+(w6) < 2 ], [ 0.5 < (-0.5*x1)+(w7)+(z9) < 1e+30 ], [ -1 < (-1*x1)+(w7)+(-1*z9) < 1e+30 ], [ -1e+30 < (-1*x1)+(w7)+(z9) < 1 ], [ -1e+30 < (-0.5*x1)+(w7)+(-1*z9) < -0.5 ];

  Implementations of this class may be numerical solvers, working towards finding a solution, or reformulation solvers, working towards analysing or changing the problem structure. Normally, Problem bounds (both variable and constraint) can be changed dynamically by a solver without the original problem bounds being modified. Numerical solvers can add both linear and nonlinear cuts (see Sect. 4.4) to the formulation before solving it. Cuts are dealt with via two auxiliary classes UserLinearCut and UserCut.Because of their simplicity, UserLinearCut and UserCut do not offer a full set/get interface, and all their properties are public. Cuts can only be added, never deleted; however, they can be enabled/disabled as needed.

	class Problem. Basic methods. class Problem. Evaluation methods.
	Method name Method name	Purpose Purpose
	Problem(bool nosimplify) double EvalObj(int oID)	constructor with optional nosimplify evaluate an objective
	void SetName(string) double EvalNLObj(int oID)	set the problem name evaluate the nonlinear part of an objective
	string GetName(void) double EvalObjDiff(int oID, int vID)	get the problem name evaluate the derivative of an obj.
	bool IsProblemContinuous(void) double EvalObjDiffNoConstant(int oID, int vID) class Solver. Basic and cut-related methods.	true if no integer variables eval. non-const. part of a deriv.
	bool IsProblemLinear(void) double EvalObjDiff2(int oID, int vID1, int vID2) Method name	true if no nonlinear expressions evaluate 2nd derivative of an obj. Purpose
	Problem* GetParent(void) double EvalConstr(int cID) string GetName(void)	get parent problem in a tree of problems evaluate a constraint solver name
	int GetNumberOfChildren(void) double EvalNLConstr(int cID) void SetProblem(Problem* p)	get number of children problems evaluate nonlinear part of a constraint set the problem for the solver
	Problem* GetChild(int pID) double EvalConstrDiff(int cID, int vID) Problem* GetProblem(void)	get pID-th child in list of children problems evaluate a constr. derivative get the problem from the solver
	string GetFormulationName(void) double EvalConstrDiffNoConstant(int cID, int vID) eval. non-const. of constr. deriv. name of reform. assigned to this prob. bool CanSolve(int probtype) can this solve a certain problem type?
	void SetOptimizationDirection(int oID, int minmax) double EvalConstrDiff2(int cID, int vID1, int vID2) evaluate 2nd constr. derivative set opt. dir. of oID-th objective (0=LP,1=MILP,2=NLP,3=MINLP)
	void SetOptimizationDirectionLI(int li, int minmax) bool IsObjConstant(int oID) void Initialize(bool force)	local index version is the objective a constant? initialize solver
	int GetOptimizationDirection(int oID) bool IsObjDiffConstant(int oID, int vID) bool IsProblemInitialized(void)	get opt. dir. of oID-th objective is the obj. derivative a constant? is solver initialized?
	int GetOptimizationDirectionLI(int li) bool IsObjDiff2Constant(int oID, int vID1, int vID2) is the 2nd obj. deriv. a const.? local index version int Solve(void) solve/reformulate the problem
	bool HasDeleted(void) bool IsConstrConstant(int cID) Parameters GetParams(void)	true if simplification deleted some entity is the constraint a constant? get the parameter set
	void SetSolved(bool s) bool IsConstrDiffConstant(int cID, int vID) Parameters& GetParamsRef(void)	mark problem as solved/unsolved is the constr. deriv. a constant? get a reference to the parameters
	bool IsSolved(void) bool IsConstrDiff2(int cID, int vID1, int vID2) void ReplaceParams(Parameters& p)	return solved/unsolved mark is the 2nd constr. deriv. a const.? replace the parameters
	void SetFeasible(int feas) int IsFeasible(void) bool IsConstrActive(int cID, double tol, int& LU) void SetOptimizationDirection(int maxmin) int GetOptimizationDirection(void)	mark problem as feasible/infeasible is the constraint active L/U bound? set 1st objective opt. dir. return feasible/infeasible mark get 1st objective opt. dir.
	Parameters GetParams(void) Parameters& GetParamsRef(void) void ReplaceParams(Parameters& prm) int GetNumberOfVariables(void) int GetNumberOfIntegerVariables(void) void GetSolution(map<int,double>& ofval, class Problem. Construction methods. Method name map<int,double>& soln) get solution returns a copy of the set of parameters returns a reference to a set of parameters replace the current set of parameters return number of variables return number of integer variables Purpose void Parse(char* file) void GetSolutionLI(vector<double>& ofval, parse a ROSE-formatted file vector<double>& soln) get solution Ampl::ASL* ParseAMPL(char** argv, int argc) parse AMPL-formatted .nl file void SetMaxNumberOfCuts(int) set max number of cuts
	int GetNumberOfObjectives(void) int GetMaxNumberOfCuts(void)	return number of objectives get max number of cuts
	int GetNumberOfConstraints(void) Variable* GetVariable(int vID) Variable* GetVariableLI(int li) Variable* GetObjective(int vID) Variable* GetObjectiveLI(int li) Variable* GetConstraint(int vID) Variable* GetConstraintLI(int li) void SetOptimalVariableValue(int vID, double val) void SetOptimalVariableValueLI(int li, double val) double GetOptimalVariableValue(int vID) double GetOptimalVariableValueLI(int li) void SetCurrentVariableValue(int vID, double val) void SetCurrentVariableValueLI(int li, double val) int GetNumberOfCuts(void) class Problem. Reformulation methods. Method name int AddCut(Expression e, double L, double U) int NewVariableID(void) int AddCut(vector<pair<int,double> >&, double L, double U) return number of constraints return pointer to variable entity local index version return pointer to objective entity local index version return pointer to constraint entity local index version set optimal variable value local index version get optimal variable value local index version set optimal variable value local index version get number of cuts added till now add a nonlienar cut Purpose add a linear cut returns next available variable ID int NewObjectiveID(void) double EvalCut(int cutID, double* xval) evaluate a cut returns next available variable ID int NewConstraintID(void) double EvalNLCut(int cutID, double* xval) evaluate the nonlinear part returns next available variable ID double EvalCutDiff(int cutID, int vID, double* xval) evaluate derivatives void AddVariable(string& n, bool i, bool pers, double EvalCutDiffNoConstant(int cutID, int vID, double* xval) as above, without constants double L, double U, double v) adds a new variable void AddObjective(string& n, Expression e, int dir, double v) double EvalCutDiff2(int cutID, int vID1, int vID2, double* xval) evaluated 2nd derivatives adds a new objective bool IsCutLinear(int cutID) is this cut linear? void AddConstraint(string& n, Expression e, double L, double U) adds a new constraint void DeleteVariable(int vID) void EnableCut(int cutID) enables a cut deletes a variable void DeleteObjective(int oID) void DisableCut(int cutID) disables a cut deletes an objective void DeleteConstraint(int cID) void SetCutLB(int cutID, double L) set lower bound deletes a constraint double GetCutLB(int cutID) get lower bound
	double GetCurrentVariableValue(int vID) void SetCutUB(int cutID, double U) class UserCut. Method name UserCut(Expression e, double L, double U) constructor get optimal variable value set upper bound Expression Function cut's body FastEvalTree* FunctionFET fast eval tree of Function Expression NonlinearPart nonlinear part of Function FastEvalTree* NonlinearPartFET fast eval tree of NonlinearPart double LB lower bound double UB upper bound bool IsLinear marks a linear cut vector<Expression> Diff derivatives vector<FastEvalTree*> DiffFET corresponding fast eval trees vector<vector<Expression> > Diff2 2nd derivatives vector<vector<FastEvalTree*> > Diff2FET corresponding fast eval trees double& disc) test feas. of current point w.r.t. a constraint bool TestConstraintsFeasibility(double tol, double& disc) test feasibility of current point in problem bool TestVariablesFeasibility(double tol, double& disc) test feasibility of current point in bounds double GetStartingPoint(int vID) get starting point embedded in the problem double GetStartingPointLI(int localindex) local index version void SetOptimalObjectiveValue(int oID, double val) double GetOptimalObjectiveValue(int oID) double GetObj1AdditiveConstant(void) int Nonzeroes number of nonzeroes in linear form get additive constant of 1st objective double UB upper bound vector<double>& soln) local index version double LB lower bound void GetSolutionLI(vector<double>& ofval, double L, double U) C-style constructor map<int,double>& soln) get solution UserLinearCut(int* varIDs, double* coeffs, int size, void GetSolution(map<int,double>& ofval, UserLinearCut(vector<pair<int,double> >&, double L, double U) C++-style constructor get optimal obj. fun. value Method name Purpose set optimal obj. fun. value solvers are initialized and then activated. class UserLinearCut. bool TestConstraintsFeasibility(int cID, double tol, Purpose double GetCurrentVariableValueLI(int li) local index version double GetCutUB(int cutID) get upper bound
	int* Varindices	variable indices in row
	double* Coeffs	coefficients of row

  • value, the value of a constant numeric term, only has meaning if oplabel is CONST;• varindex, the variable index, only has meaning if oplabel is VAR; • every term, (variables, constants and operators), can be multiplied by a numeric coefficient.

	Operand {
	int oplabel;	// operator label
	double value;	// if constant, value of constant
	long varindex;	// if variable, the variable index
	string varname;	// if variable, the variable name
	double coefficient; // terms can be multiplied by a number
	double exponent;	// leaf terms can be raised to a number
	// methods
	};
	• oplabel can be one of the following labels (the meaning of which should be clear):
	enum OperatorType {
	SUM, DIFFERENCE, PRODUCT, FRACTION, POWER, PLUS, MINUS, LOG,
	EXP, SIN, COS, TAN, COT, SINH, COSH, TANH, COTH, SQRT, VAR,
	CONST, ERROR
	};

  -51. Notes Class Operand. Call: ret = (Expression e)->MethodName(args). Expression Parse(char* buf, int& errors) parse buf and return an Expression errors is the number of parsing errors occcurred • There exist a considerable number of different constructors for Expression. See their purpose and syntax in files expression.h, tree.cxx. See examples of their usage in file expression.cxx. • Internal class methods usually return or set atomic information inside the object, or perform limited symbolic manipulation. Construction and extended manipulation of symbolic expressions have been confined to external methods. Furthermore, external methods may have any of the following characteristics:

	Class BasicExpression (inherits from Operand, Tree<BasicExpression>).
	Method name Call: ret = (Expression e)->MethodName(args). Purpose
	int GetOpType(void) Method name	returns the operator label Purpose
	double GetValue(void) string ToString(void)	returns the value of the constant leaf (takes returns infix notation expression in a string
	void Zero(void)	multiplicative coefficient and exponent sets this to zero
	void One(void)	into account) sets this to one
	double GetSimpleValue(void) bool IsEqualTo(Expression&)	returns the value (takes no notice of is this equal to the argument?
	bool IsEqualToNoCoeff(Expression&)	coefficient and exponent) [like above, ignoring multiplicative coefficient]
	long GetVarIndex(void) int NumberOfVariables(void)	returns the variable index of the variable leaf number of variables in the expression
	string GetVarName(void) double Eval(double* v, long vsize)	returns the name of the variable leaf evaluate; v[i] contains the value for variable
	double GetCoeff(void)	returns the value of the multiplicative coefficient with index i, v has length vsize
	double GetExponent(void) bool DependsOnVariable(long i)	returns the value of the exponent (for leaves) does this depend on variable i?
	void SetOpType(int) int DependsLinearlyOnVariable(long i)	sets the operator label does this depend linearly on variable i?
	void SetValue(double)	sets the numeric value of the constant leaf (0=nonlinearly, 1=linearly, 2=no dep.)
	void SetVarIndex(long) void ConsolidateProductCoeffs(void)	sets the variable index of the variable leaf if node is a product, move product of
	void SetVarName(string)	sets the name of the variable leaf all coefficients as coefficient of node
	void SetExponent(double) void DistributeCoeffOverSum(void)	sets the exponent (for leaves) if coeff. of a sum operand is not 1,
	void SetCoeff(double)	sets the multiplicative coefficient distribute it over the summands
	bool IsConstant(void) void VariableToConstant(long varindex, double c) substitute a variable with a constant c is the node a constant?
	bool IsVariable(void) void ReplaceVariable(long vi1, long vi2, string vn2) replace occurrences of variable vi1 is the node a variable?
	bool IsLeaf(void)	is the node a leaf? with variable vi2 having name vn2
	bool HasValue(double v) string FindVariableName(long vi)	is the node a constant with value v? find name of variable vi
	bool IsLessThan(double v) bool IsLinear(void)	is the node a constant with value ≤ v? is this expression linear?
	void ConsolidateValue(void) bool GetLinearInfo(...)	set value to coeff*value*exponent returns info about the linear part
	Expression Get[Pure]LinearPart(void)	and set coeff to 1 and exponent to 1 returns the linear part
	void SubstituteVariableWithConstant Expression Get[Pure]NonlinearPart(void)	returns the nonlinear part
	double RemoveAdditiveConstant(void) void Interval(...) (long int varindexMethod name Class ExpressionParser.	returns any additive constant and removes it performs interval arithmetics on the expression Purpose
	Pointer<NodeType> Copy(void) Method name	returns a copy of this node Purpose
	void SetTo(Pointer<NodeType>& t) void SetVariableID(string x, long i)	this is a reference of t assign index i to variable x;
	void SetToCopyOf(Pointer<NodeType>& t) this is a copy of t var. indices start from 1 and increase by 1
	Pointer<NodeType> long GetVariableID(string x)	return index of variable x
	operator=(Pointer<NodeType> t)	assigns a reference of t to this
	void Destroy(void)	destroys the node (collects garbage)
	Template class Tree<NodeType>. Call: ret = (Expression e)->MethodName(args).
	Method name	Purpose
	void AddNode(Pointer<NodeType>)	pushes a node at the end of the node vector
	void AddCopyOfNode(Pointer<NodeType> n) pushes a copy of node n at the end of the node vector
	bool DeleteNode(long i)	deletes the i-th node,
		returns true if successful
	void DeleteAllNodes(void)	empties the node vector
	Pointer<NodeType> GetNode(long i)	returns a reference to the
		i-th subnode
	Pointer<NodeType> * GetNodeRef(long i)	returns a pointer to the
		i-th subnode
	Pointer<NodeType> GetCopyOfNode(long i) returns a copy of the i-th subnode
	long GetSize(void)	returns the length of the
		node vector

, double c) substitute a variable with a constant c Template class Pointer<NodeType>. Call: ret = (Expression e).MethodName(args).

• The lists given above only include the most important methods. For the complete lists, see the files expression.h, tree.cxx, parser.h in the source code distribution.

  Consider an Expression e object and a pointer Expression* ePtr = &e. The following calls are possibile: e->MethodName(args); (*ePtr)->MethodName(args); Call a method in the BasicExpression, Operand or Tree<> classes. e.MethodName(args); (*ePtr).MethodName(args); ePtr->MethodName(args);

sBB, 30, 32, 35 scalability, 26 scheduling communication delays, 33 SDP, 12, 33 dual, 12 primal, 12 semantic interpretation, 28 separable, 10 sequential quadratic prorgramming, 10 signomial, 36 signomial function, 36 simplex tableau, 38 simplex method, 9, 10, 17 Slack, 17 SNOPT, 2 software

Acknowledgements

Financial support by ANR grant 07-JCJC-0151 and by the EU NEST "Morphex" project grant is gratefully acknowledged. We also wish to thank: Claudia D'Ambrosio and David Savourey for help on the ROSE implementation; Pierre Hansen, Nenad Mladenović, Frank Plastria, Hanif Sherali and Tapio Westerlund for many useful discussions and ideas; Kanika Dhyani and Fabrizio Marinelli for providing interesting application examples.

• The symbolic derivative routine Diff() uses copying and not linking of nodes when building up the derivative. • The method BasicExpression::IsEqualToNoCoeff() returns true if two expressions are equal apart from the multiplicative coefficient of the root node only. I.e., 2(x + y) would be deemed "equal" to x + y (if 2 is a multiplicative coefficient, not an operand in a product) but x + 2y would not be deemed "equal" to x + y. • The Simplify() method applies all simplification rules known to Ev3 to the expression and puts it in standard form. • The methods GetLinearInfo(), GetLinearPart(), GetPureLinearPart(), GetNonlinear-Part(), GetPureNonlinearPart() return various types of linear and nonlinear information from the expression. Details concerning these methods can be found in the Ev3 source code files expression.h, expression.cxx. • The method Interval() performs interval arithmetic on the expression. Details concerning this method can be found in the Ev3 source code files expression.h, expression.cxx. • Variables are identified by a variable index, but they also know their variable name. Variable indices are usually assigned within the ExpressionParser object, with the SetVariableID()