Global optimization of mixed-integer nonlinear - Springer Link

Computing (2012) 94:325–343 DOI 10.1007/s00607-011-0175-7

Global optimization of mixed-integer nonlinear (polynomial) programming problems: the Bernstein polynomial approach Bhagyesh V. Patil · P. S. V. Nataraj · Sharad Bhartiya

Received: 28 March 2011 / Accepted: 12 November 2011 / Published online: 23 November 2011 © Springer-Verlag 2011

Abstract In this paper, we propose an algorithm for constrained global optimization of mixed-integer nonlinear programming (MINLP) problems. The proposed algorithm uses the Bernstein polynomial form in a branch-and-bound framework. Ingredients such as continuous relaxation, branching for integer decision variables, and fathoming for each subproblem in the branch-and-bound tree are used. The performance of the proposed algorithm is tested and compared with several state-of-the-art MINLP solvers, on two sets of test problems. The results of the tests show the superiority of the proposed algorithm over the state-of-the-art solvers in terms of the chosen performance metrics. Keywords Bernstein polynomials · Branch-and-bound · Global optimization · Mixed-integer nonlinear programming Mathematics Subject Classification (2000)

90C11 · 90C26 · 65G30

The authors have presented the results of this paper during the SCAN 2010 conference in Lyon, September 2010. B. V. Patil (B) · P. S. V. Nataraj Systems and Control Engineering Group, Indian Institute of Technology Bombay, Powai, Mumbai 400 076, India e-mail: [email protected] S. Bhartiya Chemical Engineering Department, Indian Institute of Technology Bombay, Powai, Mumbai 400 076, India

123

326

B. V. Patil et al.

1 Introduction Mixed-integer nonlinear programming (MINLP) refers to the class of optimization problems that involve continuous as well as integer decision variables. MINLP problems in their general form can be represented as [1]: min x

f (x)

subject to gi (x) ≤ 0,

i = 1, 2, . . . , m

h j (x) = 0, xk ∈ x ⊆ R,

j = 1, 2, . . . , n k = 1, 2, . . . , ld

xk ∈ {0, 1, . . . , q} ⊂ Z,

(1)

k = ld + 1, . . . , l

where f : Rl → R is the (possibly nonlinear) objective function, gi : Rl → R (i = 1, 2, . . . , m), and h j : Rl → R ( j = 1, 2, . . . , n) are the (possibly nonlinear) inequality and equality constraint functions. Further, x := [x, x] is an interval in R, xk (k = 1, 2, . . . , ld ) are continuous decision variables, and the rest of xk (k = ld + 1, . . . , l) are integer decision variables with values 0 to q, q ≥ 0, q ∈ Z. We call the set of all x satisfying inequality and equality constraints, gi and h j , respectively as the feasible set D of points. We shall assume that D is nonempty. MINLP problems are known to be challenging from theoretical, algorithmic, and computational points of view [1]. Several techniques exist in literature to solve MINLP problems, such as Outer Approximation (OA) [2,3], Generalized Benders Decomposition (GBD) [4], Branch-and-Bound (BB) [5,6], and Extended Cutting Plane methods [7]. All these techniques more or less rely on obtaining successive solutions of closely related nonlinear programming (NLP) problems, and only guarantee global optimality under (generalized) convexity. Several commercial software packages to solve convex MINLPs such as GAMS/ DICOPT [8], SBB [8], MINLP_BB [9], Bonmin [10], SCICONIC [11], LaGo[12], and LOGMIP [13] have been published. Moreover, commercial solvers such as LINDOGlobal [14], COUENNE [15], BARON [16] based on global optimization methods to solve nonconvex MINLPs have also been published. Most of the above commercial solvers provide the successive outer-approximation of the initial relaxation till global optimality, or they exploit the structure of the problem and perform decomposition and reformulation to solve the reformulated problem in the branch-and-bound framework. If decomposition is not possible, then each NLP subproblem is attempted to solve to global optimality, and the success depends heavily upon the type of NLP solver used. In this work, we propose a global optimization algorithm for solving the MINLP problem, where the objective function ( f ) and constraints (gi and h j ) are restricted to be polynomial functions in x. The proposed algorithm is based on the Bernstein form of polynomials, and extends the work recently proposed for solving NLPs in [17–21]. The proposed algorithm uses tools such as continuous relaxation, branching for the integer decision variables, and fathoming for each subproblem in a branch-and-bound framework. These tools are new in the context of the Bernstein polynomial based approach to global optimization.

123

Global optimization of mixed-integer nonlinear (polynomial) programming problems

327

The main merit of the proposed approach is that it guarantees that the global minimum of the MINLP is found to a user-specified accuracy. Further merits of the proposed approach are that (i) it does not require any initial guess for starting the optimization, but only an initial box bounding the region of interest; (ii) in case of multiple global solutions in this box, it guarantees that all solutions are found; (iii) it does not require any convexification, decomposition, or linearization of the primal MINLP problem; and (iv) it does not require any prior knowledge of the stationary points. The performance of the proposed algorithm is tested and compared with some of the well-known MINLP solvers. For this purpose, two sets of test problems are used. The first set has 18 standard MINLP problems taken from literature [1,22,23], while the second set comprises 4 MINLP problems constructed by the present authors, on the lines given in [24]. The rest of the paper is organized as follows: In Sect. 2, we present a brief review of the Bernstein form, and outline a basic Bernstein global optimization algorithm for solving NLPs. In Sect. 3, we first introduce the ingredients of the proposed algorithm for MINLPs, and then present the main algorithm. In Sect. 4, we first compare the performance of the proposed algorithm for three different branching rules on integer variables. Then, we compare the performance of the proposed algorithm with those of several state-of-the-art MINLP solvers. Lastly, we give the conclusions of the work in Sect. 5. 2 The Bernstein polynomial approach 2.1 Background Let l ∈ N be the number of variables and x = (x1 , x2 , . . . , xl ) ∈ Rl . A multi-index I is defined as I = (i 1 , i 2 , . . . , il ) ∈ Nl and the multi-power x I is defined as x I = (x1i1 , x2i2 , . . . , xlil ). A multi-index N is defined as N = (n 1 , n 2 , . . . , nl ). Inequalities I ≤ N for multi-indices are meant component-wise. With I = (i 1 , . . . , ir −1 , ir , ir +1 , . . . , il ) we associate the index Ir,k given by Ir,k = (i 1 , . . . , ir −1 , ir + k, ir +1 , . . . , il ), where 0 ≤ ir + k ≤ nr . Also we write NI nl n 1 for i1 · · · il and (N /I ) for (n 1 /i 1 , n 2 /i 2 , . . . , nl /il ) provided that 0 < i k , k = 1, 2, . . . , l. A real bounded and closed interval x is defined as x = [x, x] := [inf x, sup x] ∈ IR, where IR denotes the set of compact intervals. Let w(x) denote the width of x, that is w(x) := x − x, and m(x) denote the midpoint of x, that is m(x) := (x + x)/2. Similarly, for an l-dimensional interval vector or box x = (x1 , x2 , . . . , xl ) ∈ IRl , the width of x is w(x) := max(w(x1 ), w(x2 ), . . . , w(xl )). We can write an l-variate polynomial p in the form p(x) =

aI x I ,

x ∈ Rl ,

I ≤N

123

328

B. V. Patil et al.

with N being the degree of p. We expand a given multivariate polynomial into Bernstein polynomials to obtain bounds for its range over an l -dimensional box x. The I th Bernstein basis polynomial of degree N is defined as B IN (x) = Bin11 (x1 ) · · · Binl l (xl ),

x ∈ Rl ,

where, for i j = 0, 1, . . . , n j , j = 1, 2, . . . , l i n −i nj n j (x j − x j ) j (x j − x j ) j j Bi j (x j ) = . ij (x j − x j )n j The Bernstein coefficients b I (x) of p over the box x are given by I J K w(x) J (inf x) K −J a K , b I (x) = I ≤ N. J N J ≤I K ≤J J The Bernstein form of a multivariate polynomial p is defined by p (x) = b I (x) B IN (x) .

(2)

(3)

(4)

(5)

I ≤N

The Bernstein coefficients are collected in an array (b I (x)) I ∈S , where S = {I : I ≤ N }. We denote S0 as a special subset of the index set S comprising indices of the vertices of this array, that is S0 := {0, n 1 } × {0, n 2 } × · · · × {0, nl }. Equation (4) is the basis for the conventional method of computing the Bernstein coefficients from the coefficients of the power form. Other methods for computation of the Bernstein coefficients are: Ray’s matrix method [21,25], Garloff’s method [26], the scaled coefficients method [27], the parallel scheme method [28], and Smith’s method [29]. An attractive feature of the matrix method is that the Bernstein coefficients are computed using only matrix operations such as multiplication, inverse, transpose and reshape. We shall use Ray’s matrix method for computation of the Bernstein coefficients. Theorem 1 (Range enclosure property) Let p be a polynomial of degree N , and let p(x) denote the range of p on a given box x ∈ IRl . Then, (6) p(x) ⊆ B(x) := min (b I (x)) I ∈S , max (b I (x)) I ∈S . Proof See [30]. Remark 1 The above theorem says that the minimum and maximum coefficients of the array (b I (x)) I ∈S provide lower and upper bounds for the range. This forms the Bernstein range enclosure, defined by B(x) in equation (6). Lemma 2 (Vertex property) [30] Consider the Bernstein form in equation (5) for a polynomial p of degree N , and let the range p(x) = [a, b]. Then

123


329

a = min (b I (x)) if and only if min (b I (x)) = min (b I (x)) 0≤I ≤N

0≤I ≤N

I ∈S0

b = max (b I (x)) if and only if max (b I (x)) = max (b I (x)) 0≤I ≤N

0≤I ≤N

I ∈S0

Remark 2 The above Lemma says that the lower bound (respectively upper bound) is sharp if and only if min (b I (x)) I ∈S (respectively max(b I (x)) I ∈S ) is attained at a Bernstein coefficients of the array (b I (x)) with I ∈ S0 . This condition is known as the vertex property. The following properties follow immediately from Theorem 1. Lemma 3 Let B(x) be the Bernstein range enclosure for a polynomial p(x) on a given box x ∈ IRl . Then, the following properties hold 1. 2. 3. 4.

B(x) ≤ 0 ⇒ p(x) ≤ 0 for all x ∈ x. B(x) > 0 ⇒ p(x) > 0 for all x ∈ x. 0∈ / B(x) ⇒ p(x) = 0 for all x ∈ x. B(x) ⊆ [−, ] ⇒ p(x) ∈ [−, ] for all x ∈ x, where > 0.

2.2 Basic Bernstein algorithm for constrained global optimization In this subsection, we outline the basic optimization algorithm using the Bernstein form for solving polynomial NLPs with only continuous decision variables. The algorithm is similar to the one described in [31] and uses the Bernstein form as a inclusion function (refer to Appendix B) for the optimization. The aim of the Bernstein algorithm described below is to determine numerically the global minimum f ∗ , and the set of global minimizer(s) x ∗ ∈ x of the MINLP problem (1), where f, gi , h j are all polynomial functions in x. However, on a computer with finite precision it is difficult to numerically verify the equality constraint h j (x) = 0. Hence, we replace this constraint in (1) with the relaxed constraint h j (x) ∈ [−zer o , zer o ], j = 1, 2, . . . , n where zer o > 0 is a very small number that represents the ‘numerical zero’ of the computing machine. In the sequel, we shall consider this relaxed MINLP problem. The algorithm called BernsteinGlobal generates • a set U enclosing the set of feasible points D, • the value y, a guaranteed underestimate of f ∗ , • the value p , an upper bound on f ∗ . Briefly, the algorithm works as follows. We start the algorithm by setting the working box y to the given initial box x and compute the Bernstein coefficients of the objective, inequality and equality constraint polynomials on y. We store them in arrays p to (bo (y)), bgi (y) , and bh j (y) , respectively. We also initialize the upper bound infinity, and the current minimum y to the minimum Bernstein coefficient of objective function over y. We further initialize each component of a flag vector R to zero, the working list L with the item {(y, y, R)}, and the solution list Lsol to the empty list. We then pick the first item (y, y, R) from L, and delete its entry from L. For this item we then subdivide the box y along the longest width direction creating two subboxes v1 , v2 . We compute the Bernstein coefficient arrays

123

330

B. V. Patil et al.

(bo (vk )), bgi (vk ) , bh j (vk ) , k = 1, 2 for v1 , v2 and the Bernstein range enclosures Bo (vk ), Bgi (vk ), Bh j (vk ) of the objective, inequality and equality constraint polynomials respectively. We next check the feasibility of the inequality and equality constraints for v1 , v2 by doing the following tests (cf. Lemma 3): • If Bgi (vk ) ≤ 0, 0 ∈ Bh j (vk ), Bh j (vk ) ⊆ [−zer o , zer o ], for all i = 1, . . . , m, j = 1, . . . , n, then vk is a feasible box, • If Bgi (vk ) > 0 for some i, then vk is infeasible, and can be discarded. • If 0 ∈ / Bh j (vk ), Bh j (vk ) [−zer o , zer o ] for some j, then vk is infeasible, and can be discarded. If vk survives these tests, and if each component of the flag vector R k (see Remark 3 below) is equal to unity, we update the upper bound p and add the item (vk , dk , R k ) to the end of list L such that second members of all items of the list do not decrease (here dk := min Bo (vk )). We next discard item(s) (z, z, R) from list L whose minimum value z is greater than the upper bound p , and denote the first item of the list L by (y, y, R). Further, if item (y, y, R) satisfies the termination criteria, then we deposit it in the solution list Lsol . If the termination criteria is not satisfied then we continue the algorithm until the list L becomes empty. The algorithm can be made more efficient by using a flag vector R, based on the following ideas. Remark 3 Suppose gi (x) ≤ 0 for x ∈ x, meaning that the ith inequality constraint is satisfied in x. So, for any subbox x0 ⊆ x, there is no need to again check gi (x) ≤ 0 for x ∈ x0 . The same holds true for h j (x). We handle such information with a flag vector R = (R1 , . . . , Rm , Rm+1 , . . . , Rm+n ), where the components Ri takes the values 0 or 1, as follows: • Ri = 1 if the ith constraint (inequality or equality) is satisfied for the box. • Ri = 0 if the ith constraint satisfaction has not yet been verified for the box. Algorithm Bernstein [ y, p , U ] = BernsteinGlobal(N , a I , x, p , x , zer o ) Inputs Degree N of the variables occurring in the objective and constraint polynomials, the coefficients a I of the objective and constraint polynomials in the power form, the initial search domain x, the tolerance parameters p and x on the global minimum and global minimizer(s), and the tolerance parameter zer o to which the equality constraints are to be satisfied. Outputs A lower bound y and an upper bound p on the global minimum f ∗ , along with a set U containing all the global minimizer(s) x(i) . BEGIN Algorithm 1. Set y := x. 2. From a I , compute the Bernstein coefficient arrays of the objective and constraint polynomials on the box y respectively as (bo (y)), bgi (y) , bh j (y) , i = 1, 2, . . . , m, j = 1, 2, . . . , n.

123


331

Set p := ∞ and y := min (bo (y)). Set R = (R1 , . . . , Rm , Rm+1 , . . . , Rm+n ) := (0, . . . , 0). Initialize list L := {(y, y, R)}, Lsol := {}. If L is empty then go to step 14. Otherwise, pick the first item (y, y, R) from L, and delete its entry from L. 7. Choose a coordinate direction λ parallel to which y1 × · · · × yl has an edge of maximum length, that is λ ∈ {i : w(y) := w(yi ), i = 1, 2, . . . , l}. 8. Bisect y normal to direction λ, getting boxes v1 , v2 such that y = v1 ∪ v2 . 9. for k = 1, 2 k k , Rk (a) Set R k = (R1k , . . . , Rm m+1 , . . . , Rm+n ) := R. (b) Find the Bernstein coefficient array and the corresponding Bernstein range enclosure of the objective function ( f ) over vk as (b0 (vk )) and B0 (vk ), respectively. (c) Set dk := min Bo (vk ). (d) If p < dk then go to substep (j). (e) for i = 1, 2, . . . , m if Ri = 0 then (i) Find the Bernstein coefficient array and the corresponding Bernstein range enclosure of the inequality constraint polynomial (gi ) over vk as (bgi (vk )) and Bgi (vk ), respectively. (ii) If Bgi (vk ) > 0 then go to substep (j). (iii) If Bgi (vk ) ≤ 0 then set Rik := 1. (f) for j = 1, 2, . . . , n if Rm+ j = 0 then (i) Find the Bernstein coefficient array and the corresponding Bernstein range enclosure of the equality constraint polynomial (h j ) over vk as (bh j (vk )) and Bh j (vk ), respectively. (ii) If 0 ∈ / Bh j (vk ) then go to substep (j). k (iii) If Bh j (vk ) ⊆ [−zer o , zer o ] then set Rm+ j := 1. k (g) If R = (1, . . . , 1) then set p := min( p , max Bo (vk )). (h) Enter (vk , dk , R k ) into the list L such that the second members of all items of the list do not decrease. (j) end (of k−loop). 10. Discard all items (z, z, R) in the list L that satisfy p < z. 11. Denote the first item of the list L by (y, y, R). 12. If (w(y) < x ) & (max Bo (y) − min Bo (y)) < p then remove the item from the list L and enter it into the solution list Lsol . 13. Go to step 6. 14. {Compute the global minimum} Set the global minimum y to the minimum of the second entries over all the items in Lsol . 15. {Compute the global minimizers} Find all those items in Lsol for which the second entries are equal to y. The first entries of these items contain the global minimizer(s) x(i) . 16. Return the lower bound y and upper bound p on the global minimum f ∗ , along with the set U containing all the global minimizer(s) x(i) . 3. 4. 5. 6.

END Algorithm

123

332

B. V. Patil et al.

3 Proposed algorithm for polynomial MINLPs In this section, we present an algorithm for global optimization using the Bernstein form for solving polynomial MINLPs with continuous as well as integer decision variables. The proposed algorithm extends the above described Bernstein algorithm for polynomial NLPs with continuous decision variables. To cope with integer decision variables, tools from MINLP literature such as relaxation, bounding, branching, and node selection are incorporated into the algorithm. We first briefly describe these tools and then present our proposed algorithm. 3.1 Relaxation An optimization problem P has a relaxation RP, if the set of feasible solutions F S of P is a subset of the set of feasible solutions of RP [1]: F S(P) ⊆ F S(RP). A simple relaxation of the problem P is obtained by relaxing the integrality condition on each integer variable xk as x k ≤ xk ≤ x k (x k and x k indicate lower and upper bounds on xk , respectively). 3.2 Bounding Bounding is used to compute bounds on the optimal objective function value for the subproblems to be solved. Let CS denote the candidate subproblem currently under consideration1 , RCS a relaxation of CS, z RC S the optimal value of RCS, and z inc the current incumbent (best solution obtained so far). Then, CS is said to be fathomed if any of the following two conditions are satisfied [1]: (I) No feasible solution to the RCS is found. Then, CS has no feasible solution, and hence CS is fathomed. (II) A feasible solution to the RCS is found. Then, one of the following cases arises: (i) The solution to the RCS is greater than or equal to current (incumbent) solution z inc , that is z RC S ≥ z inc . Then, z inc is the best solution obtained so far. Hence, CS is fathomed. (ii) The solution to the RCS is less than the current (incumbent) solution z inc , that is z RC S < z inc . In this case if all integer variables xk take integer values from their respective integer domains, then z inc can be updated with z RC S and CS is fathomed.

1 CS refers to an optimization problem which includes the objective and the constraint polynomials, and

the initial search domain x.

123


333

3.3 Branching Different rules exist for branching the integer variables in MINLPs [5,32–34]. We outline three of the widely used branching rules. We shall use these rules in the proposed algorithm (refer to Sect. 3.5) for branching the integer variables. (I) Maximal fractional branching [5,32] In this rule, we select that integer variable for branching for which the maximal integer violation occurs. That is, we select the branching variable that maximizes the value of the argument defined as: arg max{min(xk − xk , xk + 1 − xk )}, k

where xk are integer variables to be branched and x is f loor (x). (II) Strong branching [32,34] This rule evaluates both child nodes for all integer variables with highest priority and selects the branching variable that degrades the objective value the most. This is achieved by selecting the branching variable that maximizes the value of the scor e function defined as: + − + − + scor e(Q − k , Q k ) := (1 − μ). min{Q k , Q k } + μ. max{Q k , Q k }, + where Q − k and Q k are the two subproblems obtained by adding the inequalities xk ≤ x k and xk ≥ x k (x k is the fractional value of x for the kth integer variable, and the score factor μ is some number between 0 to 1). Though the rule is computationally expensive, it helps to obtain better branching decisions and can be useful for finding early good solutions. (III) Most infeasible branching [34] In this rule, we select that integer variable whose fractional part is closest to 0.5. That is, we select the branching variable that maximizes the value of the scor e function defined as:

s(k) := 0.5 − |x k − x k − 0.5|. 3.4 Node selection After branching the integer variable, children nodes with corresponding candidate subproblems are created. For selection of a candidate subproblem, different criteria exist [1,5,32]. Basically, there are three criteria, namely depth-first search with backtracking, breadth-first search, and the best bound search. In the first criterion, we consider any one of the children nodes and solve the corresponding candidate subproblem. Once a node is fathomed, we backtrack from this node towards the root node, until we identify a node that has a child node not considered. We shall use this criterion in the proposed algorithm.

123

334

B. V. Patil et al.

3.5 Proposed algorithm Briefly, the proposed algorithm called Bernstein Mixed-Integer Optimizer (BMIO) works as follows. At the outset, we initialize a working list Lcs with the original problem CS (the root problem), the current estimate z inc of the global minimum to infinity, and global minimizer xinc to the initial domain x. From the list Lcs , we pick the last subproblem, delete its entry from Lcs , and denote this subproblem as CS0 . We next obtain a relaxation (see Sect. 3.1) RCS0 of the candidate subproblem CS0 . We apply the algorithm BernsteinGlobal to solve RCS0 and compute the global minimum z RC S0 and global minimizers x RC S0 for the relaxed root problem. If z RC S0 meets the bounding criteria (see Sect. 3.2), and the integer decision variables take integer values from their respective domains, then we update z inc and xinc to z RC S0 and x RC S0 , respectively. If the problem is found to be feasible, but some of the integer decision variables take noninteger values, then we choose a branching variable based on the branching rule (see Sect. 3.3). We now generate two new subproblems CS1 and CS2 from CS0 , and add them at the end of list Lcs . We continue the entire process until Lcs becomes empty. At termination, the best incumbent solution z inc and xinc are returned as the global minimum z ∗ and global minimizers x∗ , respectively. Algorithm [z ∗ , x∗ ] = BMIO(CS) Inputs The original problem CS−this includes the objective and the constraint polynomials, and the initial search domain x. Outputs The global minimum z ∗ and global minimizers x∗ in the search domain x. BEGIN Algorithm 1. {Initialization} Set Lcs ← {CS}, z inc ← ∞, xinc ← x. 2. {Start new iteration} If Lcs is empty, then terminate with z ∗ ← z inc and x∗ ← xinc . Otherwise, pick the last subproblem from Lcs , delete its entry from Lcs , and denote it as CS0 . 3. {Relaxation} Obtain a relaxation RCS0 of the candidate subproblem CS0 (refer to Sect. 3.1). Solve RCS0 using the algorithm BernsteinGlobal of Sect. 2.2. Denote the obtained global minimum of the relaxed problem as z RC S0 and the global minimizers as x RC S0 . 4. {Bounding} (a) If RCS0 is found to be infeasible (z RC S0 = ∞), then CS0 has no feasible solution. Go to step 2. (b) If z RC S0 ≥ z inc , then CS0 has no feasible solution better than z inc . Go to step 2. (c) If z RC S0 < z inc and all integer variables take integer values from their respective domains, then update z inc ← z RC S0 and xinc ← x RC S0 . Go to step 2. (d) If some of the integer variables takes a noninteger value, then go to step 5. 5. {Branching} Apply a branching rule to CS0 so as to choose an integer variable for branching (refer to Sect. 3.3). Generate two new branched subproblems CS1 and

123


335

CS2 by partitioning at the branch point and add them to the end of list Lcs . Go to step 2. END Algorithm 4 Numerical tests We consider two sets of test problems. The first set has 18 test problems taken from [1,22,23]. For these test problems, we first compare the performance of the three branching rules described in Sect. 3.3. We next compare the performance of the proposed algorithm BMIO (using the best branching rule) with those of the four state-of-the-art MINLP solvers, viz AlphaECP, BARON, Bonmin, DICOPT, whose GAMS interface is available through the NEOS server [35], and one MATLAB based open-source solver BNB20 [36]. The second set comprises 4 specially constructed test problems, built on the lines suggested in [24]. For this set, we compare the performance of the proposed algorithm BMIO with the above mentioned MINLP solvers too. For all computations, we use a desktop PC with Pentium IV 2.40 GHz processor, and 2 GB RAM. The proposed algorithm BMIO is implemented in MATLAB [37]. We specify an accuracy = 10−6 for computing the global minimum and global minimizers, and allow a maximum number of 500 subdivisions in algorithm BernsteinGlobal (see Sect. 2.2). Table 1 describes the list headers for Tables 2, 3 and 5. Table 2 reports the test problems along with the references from which they were taken, problem dimensions, total number of continuous, binary and integer variables, and total number of constraints. Table 2 also gives the performance comparison of the proposed algorithm BMIO for the three different branching rules. We find that for two test problems (st_testph4, st_test1) strong branching is worse in terms of nodes searched. On the other hand, we find for one test problem (nvs02) strong branching performs better in terms of nodes searched. Similarly, we find for one test problem (zhu2) the most infeasible branching performs far better compared to the other branching rules. We find from Table 2 the average number of nodes searched by algorithm BMIO as 13, 13, and 12 for maximal fractional, strong, and most infeasible branchings respectively. The Table 1 Description of headers for Tables 2, 3, and 5

Header Example

Description Name of test problem along with reference from where it is taken

l

Number of variables (binary, integer and continuous)

lb

Number of binary variables

li

Number of integer variables

M

Total number of constraints (m + n)

f*

Global minimum obtained

Nodes

Number of nodes searched for the global minimum

Boxes

Number of boxes processed after subdivision in the Bernstein algorithm

123

336

B. V. Patil et al.

Table 2 Test problems with their characteristics and performance comparison of the proposed algorithm BMIO for different branching rules on integer variables Example

l

lb

li

M

Maximal fractional branching

Strong branching

Most infeasible branching

Nodes Boxes Time (s)

Nodes Boxes Time Nodes Boxes Time (s) (s)

floudas1[1]

3

0

1

1

1

1,003

0.45

1

1,003

0.46

1

1,003

0.45

st_e13 [22]

3

1

0

3

1

297

0.16

1

297

0.16

1

297

0.16

0.05

0.05

nvs03 [22]

3

0

2

2

1

37

0.05

1

37

1

37

zhu1 [23]

3

0

2

2 54

1,166

1.07

54

1,166

1.10 54

1,166

1.05

st_testph4 [22]

4

0

3 10 25

1,835

2.02

33

5,101

4.77 26

1,870

2.21

nvs21 [22]

4

0

2

2 18

1,149

0.82

18

1,149

0.85 18

1,149

0.81

gbd [22]

5

3

0

4

5

2,201

1.51

5

2,201

1.70

5

2,201

1.40

st_e27 [22]

5

2

0

6

1

572

0.40

1

572

0.41

1

572

0.40

0.08

0.08

ex1226 [22]

6

3

0

5

1

117

0.07

1

117

1

117

st_test1 [22]

6

0

5

1 15

2,145

1.21

17

2,870

1.45 15

2,145

1.14

zhu2 [23]

6

0

5

4 55

3,230

3.79

55

3,230

3.90 23

2,571

2.71

st_miqp4 [22]

7

0

3

4

1

1,003

0.78

1

1,003

0.79

1

1,003

0.77

alan [22]

9

4

0

7

1

1

0.05

1

1

0.04

1

1

0.04

nvs14 [22]

9

0

5

4

1

1,043

1.19

1

1,043

1.20

1

1,043

1.15

ex1225 [22]

9

6

0 10 11

6,869

6.79

11

6,869

7.02 11

6,869

6.60

gear3 [22]

9

0

4

4

8

2,024

2.94

8

2,024

2.98

2,024

2.80

9

0

5

3 43

6,304

8.44

35

5,145

11.07 43

6,304

8.41

4,016

53.45

7

4,012

50.72

3,960

48.50

nvs02 [22] st_test3 [22]

14

0 13 10

8

8 6

average computational time in seconds for algorithm BMIO are 4.73, 4.93, and 4.37 for maximal fractional, strong, and most infeasible branchings respectively. Overall, we find the performance of most infeasible branching with BMIO more satisfactory compared to the other branching rules. Table 3 reports for the 18 test problems the global minimum obtained, total number of nodes searched, and the computational time taken in seconds to find the global minimum. The bold values in the table indicate the local minimum value. Similarly, Table 4 reports the average number of nodes searched to obtain the global minimum. We find from Table 4 that BMIO visits significantly fewer nodes (on the average) than existing solvers, except for DICOPT. However, DICOPT gives a local minimum in several test problems (see floudas1, zhu2, nvs02), and is unable to solve four test problems (see st_e13, zhu1, st_miqp4, nvs14). This is perhaps because DICOPT is based upon extension of outer-approximation algorithms (as the reader is aware, in outer-approximation algorithms, we underestimate the objective function with linear support functions, and this underestimation may not be appropriate if the objective function is nonlinear [1]). We next present the performance results for the four specially constructed test problems, P1 to P4. These test problems are described in Appendix A. Table 5 summarizes

123


337

Table 3 Performance comparison of the proposed algorithm BMIO with state-of-the-art MINLP solvers Example Statistics Solver/Algorithm

f* floudas1 Nodes

st_e13

nvs03

AlphaECP

BARON

Bonmin

BNB20

DICOPT

BMIO

−8.5

−8.5

−8.5

−5

−4

−8.5

3

1

1

1

1

1

Time (s) 1.046

0.25

0.14

0.01

0.21

f*

2

2

2

2

Nodes

1

1

1

1

Time (s) 1.11

0.23

0.14

0.37

f*

16

16

16

16

16

16

Nodes

3

2

7

1

3

1

0.26

0.25

0.02

0.44

0.05

Time (s) 5.64 zhu1

f*

−3.9374E+10 −3.9374E+10 −3.9374E+10 −3.9374E+10

Nodes

8

gbd

st_e27

ex1226

st_test1

1 0.16

−3.9374E+10

−

1

Time (s) 1.36

0.25

0.16

0.07

f*

−80.5

−80.5

−80.5

−80.5

−80.5

−80.5

1

3

4

25

3

26

*

54 1.05

Time (s) 0.89

0.26

0.26

0.22

0.47

2.21

f*

−5.68

−5.68

−5.68

−5.68

−5.68

−5.68

Nodes

1

1

−

25

1

18

Time (s) 15.54

1.06

0.16

0.29

0.23

0.81

f*

2.2

2.2

2.2

2.2

2.2

2.2

Nodes

1

1

1

7

1

5

Time (s) 0.5

0.25

0.22

0.03

0.22

1.40

f*

2

2

2

2

2

2

Nodes

2

1

1

1

1

1

Time (s) 0.71

0.26

0.13

0.01

0.22

0.40

f*

−17

−17

−17

−17

−17

−17

Nodes

3

1

1

1

1

1

Time (s) 0.70

0.34

0.21

0.15

0.22

0.08

f*

0

0

0

0

0

0

Nodes

1

Time (s) 3.19 zhu2

*

1

st_testph4 Nodes

nvs21

0.45 2

1

14

17

3

15

0.38

0.40

0.17

0.69

1.14

f*

0

−50040.63

0

−42585

0

−51568

Nodes

4

1

−

23

1

23

Time (s) 1.94

0.25

0.16

1.38

0.23

2.71

f*

−4596

−4596

−4595.99

−4566

4

2

−

3

Time (s) 0.95

0.27

11.02

0.017

f*

2.92

2.92

2.92

st_miqp4 Nodes

2.92

−4596 *

1 0.77

2.92

2.92

123

338

B. V. Patil et al.

Table 3 continued Example

alan

nvs14

ex1225

gear3

nvs02

st_test3

Statistics

Solver/Algorithm AlphaECP

BARON

Bonmin

BNB20

DICOPT

BMIO 1

Nodes

1

1

1

11

2

Time (s)

0.61

0.23

0.20

0.14

1.02

f*

−30885.46

−40358.15

−37791.63

−40292

Nodes

192

7

70

35

Time (s)

42.01

0.36

0.76

0.48

0.04 −40358.05

*

1 1.15

f*

31

31

31

31

31

31

Nodes

2

1

6

25

3

11

Time (s)

0.72

0.26

0.28

0.28

0.47

6.60

f*

0

0

0

0

0

0

Nodes

4

199

130

93

12

8

Time (s)

1.61

0.83

1.094

0.69

1.812

2.80

f*

5.92

5.92

5.92

5.92

5.96

5.92

Nodes

29

7

20

1078

49

43

Time (s)

128.02

0.36

0.37

14.13

0.86

8.41

f*

−7

−7

−7

−7

−7

Nodes

3

1

6

3

6

Time (s)

0.94

0.27

0.61

0.98

48.50

**

* indicates that the solver failed giving message “relaxed NLP is infeasible” ** indicates that the solver did not give the result even after an hour and is therefore terminated − indicates that the nodes were not available for the comparison Table 4 Comparison of nodes searched with the proposed algorithm BMIO and state-of-the-art MINLP solvers Performance metric

Average number of nodes


BARON

Bonmin

BNB20

DICOPT

BMIO

15

13

15

75

5

12

the global minimum obtained, total number of nodes searched, and the computational time taken in seconds to find the global minimum. • Test problem P1 is taken from [24] wherein the correct global minimum is reported. We modified this problem restricting x1 to be an integer decision variable. While Bonmin and DICOPT are unable to solve this problem, all the other existing solvers find only a local minimum. However, the proposed algorithm BMIO is able to find the correct global minimum at the root node. • Test problem P2 is a problem constructed based on nvs03 [22]. For P2, BARON is unable to capture the correct global minimum, and BNB20 failed to solve this problem. On the other hand, AlphaECP, Bonmin, and DICOPT find the correct global minimum value, but visited more nodes than the proposed algorithm BMIO.

123


339

Table 5 Performance comparison of the proposed algorithm BMIO with state-of-the-art MINLP solvers for the constructed test problems listed in Appendix A Example

Statistics


P1

P2

P3

P4

BARON

f*

0

0

Nodes

8

1

Bonmin

BNB20

*

12

DICOPT

BMIO

***

1

16

16

−

1

0

Time (s)

1.16

0.01

f*

16

17

16

Nodes

15

3

7

3

0.02 **

0.66

Time (s)

0.95

0.06

0.20

0.12

0.12

f*

2

0

2

2

2

2 1

Nodes

1

1

1

1

4

Time (s)

0.27

0.05

0.10

0.32

0.2

f*

2

0

2

2.2056

Nodes

12

1

−

11

Time (s)

0.79

0.06

0.13

0.42

0.07 2

***

1 0.21

− indicates that the nodes were not available for the comparison * indicates that the solver failed giving the message “the LP relaxation is infeasible or too expensive” ** indicates that the solver searched 9081 nodes in 434.57 seconds, still could not find the solution *** indicates that the solver failed giving the message “relaxed NLP is infeasible”

• Test problems P3 and P4 are problems constructed by the authors. For these problems, BARON is unable to capture the correct global minimum, BNB20 and DICOPT find the correct global minimum for problem P3, but fail to capture the global minimum for problem P4. AlphaECP, Bonmin, and the proposed algorithm BMIO find the correct global minimum for both problems. 5 Conclusions We presented an algorithm for finding the global minimum of polynomial MINLPs based on the Bernstein polynomial approach. The performance of the proposed algorithm was tested and compared with some of the well-known MINLP solvers on two set of test problems. The first set comprised of 18 test problems collected from MINLP literature, while the second set comprised 4 specially constructed test problems. The results show the superiority of the proposed algorithm, in that it captures the global minimum while yielding a reduction in the average number of nodes visited. Appendix A Description of the specially constructed test problems (second set): We denote the objective function as f (x), the constraint functions as g(x), and the initial bounds as xk , k = 1, 2, . . . , l.

123

340

B. V. Patil et al.

1. Test problem P1 min f (x) = x1 subject to x2 − x12 ≥ 0 x2 − x12 (x1 − 2) + 10−5 ≤ 0 x1 ∈ Z, x2 ∈ R where x1 = [0, 10] and x2 = [−10, 10]. 2. Test problem P2 min f (x) = (x1 − 8)2 + (x2 − 2)2 subject to −0.1x12 + x2 ≥ 0 −0.333333333333333x1 − x2 ≥ −4.5 0.000002x3 (x3 − 1) ≤ 0 x1 , x2 ∈ Z, x3 ∈ R where xk = [0, 200], k = 1, 2 and x3 = [0, 1]. 3. Test problem P3 min f (x) = x1 + x2 subject to x1 − 106 ≤ 0 x2 − 0.000002(x1 − 2) ≤ 0 x1 x2 − 2 ≤ 0 x1 ∈ R, x2 ∈ Z where x1 = [−10, 10] and x2 = [0, 1]. 4. Test problem P4 min f (x) = x1 subject to x1 x2 ≤ 4 x2 −

− 2) + 10−5 ≤ 0 x1 ∈ R, x2 ∈ Z

x12 (x1

where x1 = [−4, 4] and x2 = [0, 8].

123


341

Appendix B Theoretical properties of the Bernstein algorithm Definition 1 (Interval extension) [38] Let f be a real-valued function of real variables x1 , x2 , . . . , xl and F be an intervalvalued function of interval variables x1 , x2 , . . . , xl . Then, the interval function F is said to be an interval extension of f if F(x1 , x2 , . . . , xl ) = f (x1 , x2 , . . . , xl ) for all (x1 , x2 , . . . , xl ) ∈ (x1 , x2 , . . . , xl ). Theorem 4 The Bernstein range enclosure B(x) is an interval extension of p(x). Proof See [39]. Definition 2 (Inclusion monotonicity) [38] An interval extension F is said to be inclusion monotonic if xi ⊆ yi , i = 1, 2, . . . , l ⇒ F(x1 , x2 , . . . , xl ) ⊆ F(y1 , y2 , . . . , yl ). Theorem 5 The Bernstein range enclosure B(x) is an inclusion monotonic interval extension of p(x). Proof See [39]. Definition 3 (Convergence) We say that a sequence of intervals converges to x if both the sequence of left endpoints and the sequence of right endpoints converge to x. Theorem 6 The Bernstein range enclosure B(x) converges quadratically to the range p(x) of the polynomial p on x, that is w(B(x)) − w( p(x)) ≤ λw(x)2 where λ is a constant independent of x. Proof See [39]. Corollary 7 As a polynomial p is a continuous function in x, we have w( p(x)) → 0 as w(x) → 0. Hence, by Theorem 6, w(B(x)) → 0 as w(x) → 0. Remark 4 Algorithm BernsteinGlobal (in Sect. 2.2) initializes a list L consisting of one item (y, y, R), cf. Step 5. Then, at each iteration the list is modified, cf. Steps 6 and 9(h). Let zn1 , . . . , znln be the boxes in the list Ln (the list at the n-th iteration),

n and Un = li=1 zni where ln is the length of Ln . Let yn be the current value of y in p in list Ln . list Ln , and let pn be the current value of By virtue of Theorems 4, 5, 6 and Corollary 7, we have the following result [31]:

123

342

B. V. Patil et al.

Theorem 8 (Convergence property of algorithm BernsteinGlobal) Let algorithm BernsteinGlobal be applied to the optimization problem (1) over a box x ∈ IRl . Assume that at least one global minimizer lies in the interior of D. Then, the∗ yn → f sequence Un is nested and ∞ n=1 Un = D that is, Un → D. Furthermore, with yn ≤ f ∗ and pn f ∗ as n → ∞. References 1. Floudas CA (1995) Nonlinear and mixed-integer optimization: fundamentals and applications. Oxford University Press, New York 2. Duran MA, Grossmann IE (1986) An outer approximation algorithm for a class of mixed-integer nonlinear programs. Math Program 36(3):307–339 3. Fletcher R, Leyffer S (1994) Solving mixed-integer programs by outer approximation. Math Program 66(1-3):327–349 4. Geoffrion AM (1972) A generalized Benders decomposition. J Optim Theory Appl 10(4):237–260 5. Gupta OK, Ravindran A (1985) Branch and bound experiments in convex nonlinear integer programming. Manag Sci 31(12):1533–1546 6. Quesada I, Grossmann IE (1992) An LP/NLP based branch and bound algorithm for convex MINLP optimization problems. Comput Chem Eng 16(10–11):937–947 7. Westerlund T, Pettersson F (1995) A extended cutting plane method for solving convex MINLP problems. Comput Chem Eng 19:131–136 8. GAMS Development Corp (2009) GAMS—the solver manuals. Washington, DC 9. Leyffer S (1999) User manual for MINLP_BB. University of Dundee numerical analysis report NA/XXX 10. Bonami P, Biegler LT, Conn A, Cornuejols ´ G, Grossmann IE, Laird C, Lee J, Lodi A, Margot F, Sawaya N, Wachter ¨ A (2008) An algorithmic framework for convex mixed integer nonlinear programs. Discrete Optim 5(2):186–204 11. SCICON Ltd (1989) SCICONIC user guide version 1.40. Milton Keynes, UK 12. Nowak I (2005) Relaxation and decomposition methods for mixed-integer nonlinear programming. Birkhäuser Verlag, Berlin 13. Vecchietti A, Grossmann IE (1997) LOGMIP: a disjunctive 0-1 nonlinear optimizer for process system models. Comput Chem Eng 21:S427–S432 14. Lindo systems Inc (2009) Lindo API 6.0 15. Belotti P, Lee J, Liberti L, Margot F, Wachter ¨ A (2009) Branching and bounds tightening techniques for non-convex MINLP. Optim Methods Softw 24(4-5):597–634 16. Tawarmalani M, Sahinidis NV (2002) Convexification and global optimization in continuous and mixed-integer nonlinear programming: theory, algorithms, software, and applications (nonconvex optimization and its applications). Kluwer, Dordrecht 17. Nataraj PSV, Kotecha K (2002) An algorithm for global optimization using Taylor-Bernstein form as inclusion function. J Glob Optim 24(4):417–436 18. Nataraj PSV, Kotecha K (2004) Global optimization with higher order inclusion function forms. Part 1: a combined Taylor-Bernstein form. Reliab Comput 10(1):27–44 19. Nataraj PSV, Arounassalame M (2007) A new subdivision algorithm for the Bernstein polynomial approach to global optimization. Int J Autom Comput 4(4):342–352 20. Nataraj PSV, Arounassalame M (2009) An algorithm for constrained global optimization of multivariate polynomials using the Bernstein form and John optimality conditions. Opsearch 46(2):133–152 21. Ray S, Nataraj PSV (2009) An efficient algorithm for range computation of polynomials using the Bernstein form. J Glob Optim 45(3):403–426 22. MINLP Library. http://www.gamsworld.org/minlp/minlplib/minlpstat.htm. Accessed 20 March 2010 23. Zhu W (2005) A provable better branch and bound method for a nonconvex integer quadratic programming problem. J Comput Syst Sci 70(1):107–117 24. Lebbah Y, Michel C, Rueher M (2007) An efficient and safe framework for solving optimization problems. J Comput Appl Math 199(2):372–377 25. Ray S (2007) A new approach to range computation of polynomial problems using the Bernstein form. PhD thesis, Indian Institute of Technology Bombay, India

123


343

26. Garloff J (1985) Convergent bounds for range of multivariate polynomials. In: Nickel K (ed) Interval mathematics. Lecturer notes in computer science, vol 212. Springer, Berlin, pp 37–56 27. Sanchez-Reyes ` J (2003) Algebraic manipulation in the Bernstein form made simple via convolutions. Computer-Aided Des 35(10):959–967 28. Garczarczyk ZA (2002) Parallel schemes of computation for Bernstein coefficients and their application. In: Proceedings of the international conference on parallel computing in electrical engineering, pp 334–337, Warsaw, Poland 29. Smith AP (2009) Fast construction of constant bound functions for sparse polynomials. J Glob Optim 43(2–3):445–458 30. Garloff J (1993) The Bernstein algorithm. Interval Comput 6(2):154–168 31. Ratschek H, Rokne J (1988) New computer methods for global optimization. Ellis Horwood, Chichester 32. Goux JP, Leyffer S (2002) Solving large MINLPs on computational grids. Optim Eng 3(3):327–346 33. Linderoth JT, Savelsbergh MWP (1999) A computational study of search strategies for mixed integer programming. INFORMS J Comput 11(2):173–187 34. Achterberg T, Koch T, Martin A (2005) Branching rules revisited. Oper Res Lett 33(1):42–54 35. NEOS server for optimization. http://www.neos-server.org/neos/solvers/index.html. Accessed 20 March 2010 36. Kuipers K (2009) Branch-and-bound solver for mixed-integer nonlinear optimization problems. MATLAB Central File Exchange. Retrieved 18 Dec 2009 37. The Mathworks Inc (2005) MATLAB version 7.1 (R14) 38. Moore RE (1979) Methods and applications of interval analysis. SIAM, Philadelphia 39. Stahl V (1995) Interval methods for bounding the range of polynomials and solving systems of nonlinear equations. PhD thesis, Johannes Kepler University, Linz

123