An exact algorithm for separable non-convex MINLPs Claudia DAmbrosio - - PowerPoint PPT Presentation

An exact algorithm for separable non-convex MINLPs

Claudia D’Ambrosio ¶ Jon Lee § Andreas W¨ achter §

¶University of Bologna §IBM T.J. Watson Research Center, NY. Part of this research was carried out while the first author was intern at the IBM T.J. Watson Research Center whose support is gratefully acknowledged.

Aussois 2010, 7 January, 2010

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 1 / 19

Introduction

Motivation: Real world problems;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 2 / 19

Introduction

Motivation: Real world problems; Non-convex MINLP problems are challenging.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 2 / 19

Introduction

Motivation: Real world problems; Non-convex MINLP problems are challenging. Choose between... Global solution:

Convex under-estimators and spatial branch-and-bound;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 2 / 19

Introduction

Motivation: Real world problems; Non-convex MINLP problems are challenging. Choose between... Global solution:

Convex under-estimators and spatial branch-and-bound;

“Limited” running time:

Local search, linearization and approximation.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 2 / 19

The class of problems

min

j∈N Cjxj

subject to f (x) ≤ 0 ; ri(x) +

k∈H(i) gik(xk) ≤ 0 , ∀i ∈ M ;

Lj ≤ xj ≤ Uj , ∀j ∈ N ; xj integer, ∀j ∈ I , (P) where: f : Rn → Rp and ri : Rn → R ∀i ∈ M , are convex functions, H(i) ⊆ N ∀i ∈ M , the gik : R → R are non-convex univariate function ∀i ∈ M, and I ⊆ N. For j ∈ H we assume that the bounds are finite.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 3 / 19

The SC-MINLP (Sequential Convex MINLP) framework

Init Lower bounding relaxation Q Upper bounding restriction R Refinement MATLAB MINLP solver NLP solver

AMPL

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 4 / 19

The Lower Bounding problem Q: step 1

For simplicity, let us consider a term of the form g(xk) := gik(xk): g : R → R is a univariate non-convex function of xk , for some k (1 ≤ k ≤ n).

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 5 / 19

The Lower Bounding problem Q: step 1

For simplicity, let us consider a term of the form g(xk) := gik(xk): g : R → R is a univariate non-convex function of xk , for some k (1 ≤ k ≤ n).

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

Automatically detect the concavity/convexity intervals or piecewise definition: [Pp−1, Pp] := the p-th subinterval of the domain of g (p ∈ {1 . . . p});

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 5 / 19

The Lower Bounding problem Q: step 1

For simplicity, let us consider a term of the form g(xk) := gik(xk): g : R → R is a univariate non-convex function of xk , for some k (1 ≤ k ≤ n).

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

Automatically detect the concavity/convexity intervals or piecewise definition: [Pp−1, Pp] := the p-th subinterval of the domain of g (p ∈ {1 . . . p}); ˇ H := the set of indices of subintervals on which g is convex;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 5 / 19

The Lower Bounding problem Q: step 1

For simplicity, let us consider a term of the form g(xk) := gik(xk): g : R → R is a univariate non-convex function of xk , for some k (1 ≤ k ≤ n).

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

Automatically detect the concavity/convexity intervals or piecewise definition: [Pp−1, Pp] := the p-th subinterval of the domain of g (p ∈ {1 . . . p}); ˇ H := the set of indices of subintervals on which g is convex; ˆ H := the set of indices of subintervals on which g is concave.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 5 / 19

The Lower Bounding problem Q: step 2

Replace the term g(xk) with: p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) ,

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 6 / 19

The Lower Bounding problem Q: step 2

Replace the term g(xk) with: p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) ,

and we include the following set of new constraints: P0 + p

p=1 δp − xk = 0 ;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 6 / 19

The Lower Bounding problem Q: step 2

Replace the term g(xk) with: p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) ,

and we include the following set of new constraints: P0 + p

p=1 δp − xk = 0 ;

δp − (Pp − Pp−1)zp ≥ 0 , ∀p ∈ ˇ H ∪ ˆ H ; δp − (Pp − Pp−1)zp−1 ≤ 0 , ∀p ∈ ˇ H ∪ ˆ H ; 0 ≤ δp ≤ Pp − Pp−1, ∀p ∈ {1, . . . , ¯ p}; with two dummy variables z0 := 1 and zp := 0 and two new sets of variables zp (binary) and δp (continuous).

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 6 / 19

The Lower Bounding problem Q: step 2

Replace the term g(xk) with: p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) ,

and we include the following set of new constraints: P0 + p

p=1 δp − xk = 0 ;

δp − (Pp − Pp−1)zp ≥ 0 , ∀p ∈ ˇ H ∪ ˆ H ; δp − (Pp − Pp−1)zp−1 ≤ 0 , ∀p ∈ ˇ H ∪ ˆ H ; 0 ≤ δp ≤ Pp − Pp−1, ∀p ∈ {1, . . . , ¯ p}; with two dummy variables z0 := 1 and zp := 0 and two new sets of variables zp (binary) and δp (continuous). If xk ∈ [Pp∗−1, Pp∗]

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 6 / 19

The Lower Bounding problem Q: step 2

Replace the term g(xk) with: p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) ,

and we include the following set of new constraints: P0 + p

p=1 δp − xk = 0 ;

δp − (Pp − Pp−1)zp ≥ 0 , ∀p ∈ ˇ H ∪ ˆ H ; δp − (Pp − Pp−1)zp−1 ≤ 0 , ∀p ∈ ˇ H ∪ ˆ H ; 0 ≤ δp ≤ Pp − Pp−1, ∀p ∈ {1, . . . , ¯ p}; with two dummy variables z0 := 1 and zp := 0 and two new sets of variables zp (binary) and δp (continuous). If xk ∈ [Pp∗−1, Pp∗] then zp = 1 for p = 1, . . . , p∗ − 1 and

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 6 / 19

The Lower Bounding problem Q: step 2

Replace the term g(xk) with: p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) ,

and we include the following set of new constraints: P0 + p

p=1 δp − xk = 0 ;

δp − (Pp − Pp−1)zp ≥ 0 , ∀p ∈ ˇ H ∪ ˆ H ; δp − (Pp − Pp−1)zp−1 ≤ 0 , ∀p ∈ ˇ H ∪ ˆ H ; 0 ≤ δp ≤ Pp − Pp−1, ∀p ∈ {1, . . . , ¯ p}; with two dummy variables z0 := 1 and zp := 0 and two new sets of variables zp (binary) and δp (continuous). If xk ∈ [Pp∗−1, Pp∗] then zp = 1 for p = 1, . . . , p∗ − 1 and zp = 0 for p = p∗, . . . , ¯ p,

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 6 / 19

The Lower Bounding problem Q: step 2

Replace the term g(xk) with: p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) ,

and we include the following set of new constraints: P0 + p

p=1 δp − xk = 0 ;

δp − (Pp − Pp−1)zp ≥ 0 , ∀p ∈ ˇ H ∪ ˆ H ; δp − (Pp − Pp−1)zp−1 ≤ 0 , ∀p ∈ ˇ H ∪ ˆ H ; 0 ≤ δp ≤ Pp − Pp−1, ∀p ∈ {1, . . . , ¯ p}; with two dummy variables z0 := 1 and zp := 0 and two new sets of variables zp (binary) and δp (continuous). If xk ∈ [Pp∗−1, Pp∗] then zp = 1 for p = 1, . . . , p∗ − 1 and zp = 0 for p = p∗, . . . , ¯ p, δp = Pp − Pp−1 for p = 1, . . . , p∗ − 1,

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 6 / 19

The Lower Bounding problem Q: step 2

Replace the term g(xk) with: p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) ,

and we include the following set of new constraints: P0 + p

p=1 δp − xk = 0 ;

δp − (Pp − Pp−1)zp ≥ 0 , ∀p ∈ ˇ H ∪ ˆ H ; δp − (Pp − Pp−1)zp−1 ≤ 0 , ∀p ∈ ˇ H ∪ ˆ H ; 0 ≤ δp ≤ Pp − Pp−1, ∀p ∈ {1, . . . , ¯ p}; with two dummy variables z0 := 1 and zp := 0 and two new sets of variables zp (binary) and δp (continuous). If xk ∈ [Pp∗−1, Pp∗] then zp = 1 for p = 1, . . . , p∗ − 1 and zp = 0 for p = p∗, . . . , ¯ p, δp = Pp − Pp−1 for p = 1, . . . , p∗ − 1, 0 ≤ δp∗ ≤ Pp∗ − Pp∗−1, and

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 6 / 19

The Lower Bounding problem Q: step 2

Replace the term g(xk) with: p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) ,

and we include the following set of new constraints: P0 + p

p=1 δp − xk = 0 ;

δp − (Pp − Pp−1)zp ≥ 0 , ∀p ∈ ˇ H ∪ ˆ H ; δp − (Pp − Pp−1)zp−1 ≤ 0 , ∀p ∈ ˇ H ∪ ˆ H ; 0 ≤ δp ≤ Pp − Pp−1, ∀p ∈ {1, . . . , ¯ p}; with two dummy variables z0 := 1 and zp := 0 and two new sets of variables zp (binary) and δp (continuous). If xk ∈ [Pp∗−1, Pp∗] then zp = 1 for p = 1, . . . , p∗ − 1 and zp = 0 for p = p∗, . . . , ¯ p, δp = Pp − Pp−1 for p = 1, . . . , p∗ − 1, 0 ≤ δp∗ ≤ Pp∗ − Pp∗−1, and δp = 0 for p = p∗ + 1, . . . , ¯ p.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 6 / 19

The Lower Bounding problem Q: step 2

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

0 ≤ δ1 ≤ (1 − 0) = 1, 0 ≤ δ2 ≤ (2 − 1) = 1, 0 ≤ δ3 ≤ (4 − 2) = 2.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 7 / 19

The Lower Bounding problem Q: step 2

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

0 ≤ δ1 ≤ (1 − 0) = 1, 0 ≤ δ2 ≤ (2 − 1) = 1, 0 ≤ δ3 ≤ (4 − 2) = 2. If x = 1.5 then δ1 = 1, δ2 = 0.5, δ3 = 0.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 7 / 19

The Lower Bounding problem Q: step 2

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

0 ≤ δ1 ≤ (1 − 0) = 1, 0 ≤ δ2 ≤ (2 − 1) = 1, 0 ≤ δ3 ≤ (4 − 2) = 2. If x = 1.5 then δ1 = 1, δ2 = 0.5, δ3 = 0.Then g(x) =

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 7 / 19

The Lower Bounding problem Q: step 2

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

0 ≤ δ1 ≤ (1 − 0) = 1, 0 ≤ δ2 ≤ (2 − 1) = 1, 0 ≤ δ3 ≤ (4 − 2) = 2. If x = 1.5 then δ1 = 1, δ2 = 0.5, δ3 = 0.Then g(x) = p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp)

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 7 / 19

The Lower Bounding problem Q: step 2

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

0 ≤ δ1 ≤ (1 − 0) = 1, 0 ≤ δ2 ≤ (2 − 1) = 1, 0 ≤ δ3 ≤ (4 − 2) = 2. If x = 1.5 then δ1 = 1, δ2 = 0.5, δ3 = 0.Then g(x) = p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) = g(0 + δ1) + g(1 + δ2) +

g(2 + δ3) − g(1) − g(2)

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 7 / 19

The Lower Bounding problem Q: step 2

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

0 ≤ δ1 ≤ (1 − 0) = 1, 0 ≤ δ2 ≤ (2 − 1) = 1, 0 ≤ δ3 ≤ (4 − 2) = 2. If x = 1.5 then δ1 = 1, δ2 = 0.5, δ3 = 0.Then g(x) = p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) = g(0 + δ1) + g(1 + δ2) +

g(2 + δ3) − g(1) − g(2) = g(1) + g(1.5) + g(2) − g(1) − g(2)

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 7 / 19

The Lower Bounding problem Q: step 2

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

0 ≤ δ1 ≤ (1 − 0) = 1, 0 ≤ δ2 ≤ (2 − 1) = 1, 0 ≤ δ3 ≤ (4 − 2) = 2. If x = 1.5 then δ1 = 1, δ2 = 0.5, δ3 = 0.Then g(x) = p

p=1 g(Pp−1 + δp) − p−1 p=1 g(Pp) = g(0 + δ1) + g(1 + δ2) +

g(2 + δ3) − g(1) − g(2) = g(1) + g(1.5) + g(2) − g(1) − g(2) = g(1.5)

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 7 / 19

The Lower Bounding problem Q: step 3

Still nonconvex;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 8 / 19

The Lower Bounding problem Q: step 3

Still nonconvex; Use piece-wise linear approximation for the concave intervals:

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 8 / 19

The Lower Bounding problem Q: step 3

Still nonconvex; Use piece-wise linear approximation for the concave intervals:

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 8 / 19

The Lower Bounding problem Q: step 3

Still nonconvex; Use piece-wise linear approximation for the concave intervals:

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 8 / 19

The Lower Bounding problem Q: the convex MINLP model

Replace the term g(xk) with:

• p∈ ˇ

H g(Pp−1 + δp) + p∈ ˆ H

• b∈Bp g(Xp,b) αp,b − p−1

p=1 g(Pp) ,

and we include the following set of new constraints: P0 + p

p=1 δp − xk = 0 ;

δp − (Pp − Pp−1)zp ≥ 0 , ∀p ∈ ˇ H ∪ ˆ H ; δp − (Pp − Pp−1)zp−1 ≤ 0 , ∀p ∈ ˇ H ∪ ˆ H ;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 9 / 19

The Lower Bounding problem Q: the convex MINLP model

Replace the term g(xk) with:

• p∈ ˇ

H g(Pp−1 + δp) + p∈ ˆ H

• b∈Bp g(Xp,b) αp,b − p−1

p=1 g(Pp) ,

and we include the following set of new constraints: P0 + p

p=1 δp − xk = 0 ;

δp − (Pp − Pp−1)zp ≥ 0 , ∀p ∈ ˇ H ∪ ˆ H ; δp − (Pp − Pp−1)zp−1 ≤ 0 , ∀p ∈ ˇ H ∪ ˆ H ; Pp−1 + δp −

b∈Bp Xp,b αp,b = 0 , ∀p ∈ ˆ

H ;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 9 / 19

The Lower Bounding problem Q: the convex MINLP model

Replace the term g(xk) with:

• p∈ ˇ

H g(Pp−1 + δp) + p∈ ˆ H

• b∈Bp g(Xp,b) αp,b − p−1

p=1 g(Pp) ,

and we include the following set of new constraints: P0 + p

p=1 δp − xk = 0 ;

δp − (Pp − Pp−1)zp ≥ 0 , ∀p ∈ ˇ H ∪ ˆ H ; δp − (Pp − Pp−1)zp−1 ≤ 0 , ∀p ∈ ˇ H ∪ ˆ H ; Pp−1 + δp −

b∈Bp Xp,b αp,b = 0 , ∀p ∈ ˆ

H ;

• b∈Bp αp,b = 1 , ∀p ∈ ˆ

H ; {αp,b : b ∈ Bp} := SOS2 , ∀p ∈ ˆ H . with two dummy variables z0 := 1, zp := 0 and the new set of variables αp,b.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 9 / 19

The Upper Bounding problem R

Upper Bound of the original problem:

1 The integer variables are fixed; Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 10 / 19

The Upper Bounding problem R

Upper Bound of the original problem:

1 The integer variables are fixed; 2 We solve the resulting non-convex NLP problem to local optimality;

min

j∈N Cjxj

subject to f (x) ≤ 0 ; ri(x) +

k∈H(i) gik(xk) ≤ 0 , ∀i ∈ M ;

Lj ≤ xj ≤ Uj , ∀j ∈ N ; xj = xj, ∀j ∈ I . (R)

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 10 / 19

SC-MINLP (Sequential Convex MINLP) Algorithm

Choose tolerances ε, εfeas > 0; initialize LB := −∞; UB := +∞; Find Pi

p , ˆ

Hi, ˇ Hi, X i

pb (∀i ∈ M, p ∈ {1 . . . pi}, b ∈ Bi p). Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 11 / 19

SC-MINLP (Sequential Convex MINLP) Algorithm

Choose tolerances ε, εfeas > 0; initialize LB := −∞; UB := +∞; Find Pi

p , ˆ

Hi, ˇ Hi, X i

pb (∀i ∈ M, p ∈ {1 . . . pi}, b ∈ Bi p).

repeat Solve the convex MINLP relaxation Q of the original problem P to obtain x; if (val(Q) > LB) then LB := val(Q); if (x is feasible for the original problem P (within tolerance εfeas)) then return x end if end if

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 11 / 19

SC-MINLP (Sequential Convex MINLP) Algorithm

Choose tolerances ε, εfeas > 0; initialize LB := −∞; UB := +∞; Find Pi

p , ˆ

Hi, ˇ Hi, X i

pb (∀i ∈ M, p ∈ {1 . . . pi}, b ∈ Bi p).

repeat Solve the convex MINLP relaxation Q of the original problem P to obtain x; if (val(Q) > LB) then LB := val(Q); if (x is feasible for the original problem P (within tolerance εfeas)) then return x end if end if Solve the non-convex NLP restriction R of the original problem P to obtain x; if (solution x could be computed and val(R) < UB) then UB := val(R); xUB := x end if

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 11 / 19

SC-MINLP (Sequential Convex MINLP) Algorithm

Choose tolerances ε, εfeas > 0; initialize LB := −∞; UB := +∞; Find Pi

p , ˆ

Hi, ˇ Hi, X i

pb (∀i ∈ M, p ∈ {1 . . . pi}, b ∈ Bi p).

repeat Solve the convex MINLP relaxation Q of the original problem P to obtain x; if (val(Q) > LB) then LB := val(Q); if (x is feasible for the original problem P (within tolerance εfeas)) then return x end if end if Solve the non-convex NLP restriction R of the original problem P to obtain x; if (solution x could be computed and val(R) < UB) then UB := val(R); xUB := x end if if (UB − LB > ε) then Update Bi

p , X i pb ;

end if

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 11 / 19

SC-MINLP (Sequential Convex MINLP) Algorithm

Choose tolerances ε, εfeas > 0; initialize LB := −∞; UB := +∞; Find Pi

p , ˆ

Hi, ˇ Hi, X i

pb (∀i ∈ M, p ∈ {1 . . . pi}, b ∈ Bi p).

repeat Solve the convex MINLP relaxation Q of the original problem P to obtain x; if (val(Q) > LB) then LB := val(Q); if (x is feasible for the original problem P (within tolerance εfeas)) then return x end if end if Solve the non-convex NLP restriction R of the original problem P to obtain x; if (solution x could be computed and val(R) < UB) then UB := val(R); xUB := x end if if (UB − LB > ε) then Update Bi

p , X i pb ;

end if until ((UB − LB ≤ ε) or (time or iteration limited exceeded))

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 11 / 19

SC-MINLP (Sequential Convex MINLP) Algorithm

Choose tolerances ε, εfeas > 0; initialize LB := −∞; UB := +∞; Find Pi

p , ˆ

Hi, ˇ Hi, X i

pb (∀i ∈ M, p ∈ {1 . . . pi}, b ∈ Bi p).

repeat Solve the convex MINLP relaxation Q of the original problem P to obtain x; if (val(Q) > LB) then LB := val(Q); if (x is feasible for the original problem P (within tolerance εfeas)) then return x end if end if Solve the non-convex NLP restriction R of the original problem P to obtain x; if (solution x could be computed and val(R) < UB) then UB := val(R); xUB := x end if if (UB − LB > ε) then Update Bi

p , X i pb ;

end if until ((UB − LB ≤ ε) or (time or iteration limited exceeded)) return the current best solution xUB

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 11 / 19

Refining the Lower Bounding problem

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

Add a breakpoint where the solution of problem Q of the previous iteration lies (global convergence);

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 12 / 19

Refining the Lower Bounding problem

0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5

Add a breakpoint where the solution of problem Q of the previous iteration lies (global convergence); Add a breakpoint where the solution of problem R of the previous iteration lies (speed up the convergence).

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 12 / 19

Convergence Analysis

Assumptions (l the iteration counter for the repeat loop):

• A1. f (x) and ri(x) continuous; gik in (P) uniformly Lipschitz-continuous

with a bounded Lipschitz constant Lg.

• A2. P has a feasible point. For each l, globally solve Ql.
• A3. Adds a breakpoint for every Q problem solution xl.
• A4. εfeas = 0, ε = 0 and no time/iteration limit is set.

Theorem

Under assumptions A1-A4, Algorithm SC-MINLP either terminates at a global solution of the original problem P, or each limit point of the sequence {xl}∞

l=1 is a global solution of P.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 13 / 19

Convergence Analysis

Assumptions (l the iteration counter for the repeat loop):

• A1. f (x) and ri(x) continuous; gik in (P) uniformly Lipschitz-continuous

with a bounded Lipschitz constant Lg.

• A2. P has a feasible point. For each l, globally solve Ql.
• A3. Adds a breakpoint for every Q problem solution xl.
• A4. εfeas = 0, ε = 0 and no time/iteration limit is set.

Theorem

Under assumptions A1-A4, Algorithm SC-MINLP either terminates at a global solution of the original problem P, or each limit point of the sequence {xl}∞

l=1 is a global solution of P.

The basic idea: at each iteration, the “error” of problem Q is shrinked because of the first refinement rule.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 13 / 19

Computational Results: Intro

Framework implemented under the AMPL environment;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 14 / 19

Computational Results: Intro

Framework implemented under the AMPL environment; Automatically detect the intervals of concavity/convexity using Matlab;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 14 / 19

Computational Results: Intro

Framework implemented under the AMPL environment; Automatically detect the intervals of concavity/convexity using Matlab; A solver for general convex MINLPs and a solver for NLPs used as black-box;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 14 / 19

Computational Results: Intro

Framework implemented under the AMPL environment; Automatically detect the intervals of concavity/convexity using Matlab; A solver for general convex MINLPs and a solver for NLPs used as black-box; Three classes of MINLPs of the form P (real-world applications);

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 14 / 19

Computational Results: Intro

Framework implemented under the AMPL environment; Automatically detect the intervals of concavity/convexity using Matlab; A solver for general convex MINLPs and a solver for NLPs used as black-box; Three classes of MINLPs of the form P (real-world applications); Selected and reformulated instances of GlobalLib and MinlpLib.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 14 / 19

Computational Results

Intel Core2 CPU 6600, 2.40 GHz, 1.94 GB of RAM. AMPL, Matlab R2007a, Ipopt, Bonmin. εfeas = 10−4, ε = 10−5, time limit of 2 hours.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 15 / 19

Computational Results

Intel Core2 CPU 6600, 2.40 GHz, 1.94 GB of RAM. AMPL, Matlab R2007a, Ipopt, Bonmin. εfeas = 10−4, ε = 10−5, time limit of 2 hours.

instance var/int/cons iter # LB UB int change time MINLP # br added Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 15 / 19

Computational Results

Intel Core2 CPU 6600, 2.40 GHz, 1.94 GB of RAM. AMPL, Matlab R2007a, Ipopt, Bonmin. εfeas = 10−4, ε = 10−5, time limit of 2 hours.

instance var/int/cons iter # LB UB int change time MINLP # br added ufl 1 153/39/228 1 4,122.000 4,330.400

• 1.35
• . . .

2 4,324.780 4,330.400 no 11.84 11 . . . 3 4,327.724 4,330.400 no 19.17 5 . . . 4 4,328.993 4,330.400 no 30.75 5 205/65/254 5 4,330.070 4,330.400 no 45.42 5 Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 15 / 19

Computational Results

Intel Core2 CPU 6600, 2.40 GHz, 1.94 GB of RAM. AMPL, Matlab R2007a, Ipopt, Bonmin. εfeas = 10−4, ε = 10−5, time limit of 2 hours.

instance var/int/cons iter # LB UB int change time MINLP # br added ufl 1 153/39/228 1 4,122.000 4,330.400

• 1.35
• . . .

2 4,324.780 4,330.400 no 11.84 11 . . . 3 4,327.724 4,330.400 no 19.17 5 . . . 4 4,328.993 4,330.400 no 30.75 5 205/65/254 5 4,330.070 4,330.400 no 45.42 5 ufl 2 189/57/264 1 27,516.600 27,516.569

• 4.47
• Claudia D’Ambrosio (University of Bologna)

Aussois, 7 January, 2010 15 / 19

Computational Results

Intel Core2 CPU 6600, 2.40 GHz, 1.94 GB of RAM. AMPL, Matlab R2007a, Ipopt, Bonmin. εfeas = 10−4, ε = 10−5, time limit of 2 hours.

instance var/int/cons iter # LB UB int change time MINLP # br added ufl 1 153/39/228 1 4,122.000 4,330.400

• 1.35
• . . .

2 4,324.780 4,330.400 no 11.84 11 . . . 3 4,327.724 4,330.400 no 19.17 5 . . . 4 4,328.993 4,330.400 no 30.75 5 205/65/254 5 4,330.070 4,330.400 no 45.42 5 ufl 2 189/57/264 1 27,516.600 27,516.569

• 4.47
• ufl 3

79/21/101 1 1,947.883 2,756.890

• 2.25
• . . .

2 2,064.267 2,756.890 no 2.75 2 87/25/105 3 2,292.743 2,292.777 no 3.06 2 Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 15 / 19

Computational Results

Intel Core2 CPU 6600, 2.40 GHz, 1.94 GB of RAM. AMPL, Matlab R2007a, Ipopt, Bonmin. εfeas = 10−4, ε = 10−5, time limit of 2 hours.

instance var/int/cons iter # LB UB int change time MINLP # br added ufl 1 153/39/228 1 4,122.000 4,330.400

• 1.35
• . . .

2 4,324.780 4,330.400 no 11.84 11 . . . 3 4,327.724 4,330.400 no 19.17 5 . . . 4 4,328.993 4,330.400 no 30.75 5 205/65/254 5 4,330.070 4,330.400 no 45.42 5 ufl 2 189/57/264 1 27,516.600 27,516.569

• 4.47
• ufl 3

79/21/101 1 1,947.883 2,756.890

• 2.25
• . . .

2 2,064.267 2,756.890 no 2.75 2 87/25/105 3 2,292.743 2,292.777 no 3.06 2 hydro 1 324/142/445 1

• 10,231.039
• 10,140.763
• 18.02
• 332/146/449

2

• 10,140.760
• 10,140.763

no 23.62 4 hydro 2 324/142/445 1

• 3,950.697
• 3,891.224
• 21.73
• . . .

2

• 3,950.583
• 3,891.224

no 21.34 2 . . . 3

• 3,950.583
• 3,891.224

no 27.86 2 336/148/451 4

• 3,932.182
• 3,932.182

no 38.20 2 hydro 3 324/142/445 1

• 4,753.849
• 4,634.409
• 59.33
• . . .

2

• 4,719.927
• 4,660.189

no 96.93 4 336/148/451 3

• 4,710.734
• 4,710.734

yes 101.57 2 Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 15 / 19

Computational Results

Intel Core2 CPU 6600, 2.40 GHz, 1.94 GB of RAM. AMPL, Matlab R2007a, Ipopt, Bonmin. εfeas = 10−4, ε = 10−5, time limit of 2 hours.

instance var/int/cons iter # LB UB int change time MINLP # br added nck 20 100 144/32/205 1

• 162.444
• 159.444
• 0.49
• 146/33/206

2

• 159.444
• 159.444
• 0.94

1 nck 20 200 144/32/205 1

• 244.015
• 238.053
• 0.67
• . . .

2

• 241.805
• 238.053
• 0.83

1 . . . 3

• 241.348
• 238.053
• 1.16

1 . . . 4

• 240.518
• 238.053
• 1.35

1 . . . 5

• 239.865
• 238.053
• 1.56

1 . . . 6

• 239.744
• 238.053
• 1.68

1 156/38/211 7

• 239.125
• 239.125
• 1.81

1 nck 20 450 144/32/205 1

• 391.499
• 391.337
• 0.79
• 146/32/206

2

• 391.364
• 391.337
• 0.87

1 nck 50 400 356/78/507 1

• 518.121
• 516.947
• 4.51
• . . .

2

• 518.057
• 516.947
• 14.94

2 . . . 3

• 517.837
• 516.947
• 23.75

2 . . . 4

• 517.054
• 516.947
• 25.07

2 372/86/515 5

• 516.947
• 516.947
• 31.73

2 nck 100 35 734/167/1035 1

• 83.580
• 79.060
• 3.72
• . . .

2

• 82.126
• 81.638
• 21.70

2 . . . 3

• 82.077
• 81.638
• 6.45

2 744/172/1040 4

• 81.638
• 81.638
• 11.19

1 nck 100 80 734/167/1035 1

• 174.841
• 171.024
• 6.25
• . . .

2

• 173.586
• 172.631
• 24.71

2 742/171/1039 3

• 172.632
• 172.632
• 12.85

2 Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 16 / 19

Comparison

SC-MINLP COUENNE var/int/cons time time instance

• riginal

(LB) UB (LB) UB Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 17 / 19

Comparison

SC-MINLP COUENNE var/int/cons time time instance

• riginal

(LB) UB (LB) UB ufl 1 45/3/48 116.47 4,330.400 529.49 4,330.400 ufl 2 45/3/48 17.83 27,516.569 232.85 27,516.569 ufl 3 32/2/36 8.44 2,292.777 0.73 2,292.775 Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 17 / 19

Comparison

SC-MINLP COUENNE var/int/cons time time instance

• riginal

(LB) UB (LB) UB ufl 1 45/3/48 116.47 4,330.400 529.49 4,330.400 ufl 2 45/3/48 17.83 27,516.569 232.85 27,516.569 ufl 3 32/2/36 8.44 2,292.777 0.73 2,292.775 hydro 1 124/62/165 107.77 -10,140.763 (-11,229.80) -10,140.763 hydro 2 124/62/165 211.79

• 3,932.182 (-12,104.40)
• 2,910.910

hydro 3 124/62/165 337.77

• 4,710.734 (-12,104.40)
• 3,703.070

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 17 / 19

Comparison

SC-MINLP COUENNE var/int/cons time time instance

• riginal

(LB) UB (LB) UB ufl 1 45/3/48 116.47 4,330.400 529.49 4,330.400 ufl 2 45/3/48 17.83 27,516.569 232.85 27,516.569 ufl 3 32/2/36 8.44 2,292.777 0.73 2,292.775 hydro 1 124/62/165 107.77 -10,140.763 (-11,229.80) -10,140.763 hydro 2 124/62/165 211.79

• 3,932.182 (-12,104.40)
• 2,910.910

hydro 3 124/62/165 337.77

• 4,710.734 (-12,104.40)
• 3,703.070

nck 20 100 40/0/21 15.76

• 159.444

3.29

• 159.444

nck 20 200 40/0/21 23.76

• 239.125

(-352.86)

• 238.053

nck 20 450 40/0/21 15.52

• 391.337

(-474.606)

• 383.149

nck 50 400 100/0/51 134.25

• 516.947

(-1020.73)

• 497.665

nck 100 35 200/0/101 110.25

• 81.638

90.32

• 81.638

nck 100 80 200/0/101 109.22

• 172.632

(-450.779)

• 172.632

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 17 / 19

Comparison

SC-MINLP COUENNE reformulated var/int/cons time time instance

• riginal

(LB) UB (LB) UB ex14 2 1 122/0/124 21.81 0.000 0.18 0.000 ex14 2 2 56/0/57 12.86 0.000 0.10 0.000 ex14 2 6 164/0/166 42.56 0.000 1.03 0.000 ex14 2 7 277/0/280 128.70 0.000 0.63 0.000 ex2 1 1 12/0/8 2.84

• 17.000

0.16

• 17.000

ex2 1 2 14/0/10 1.75

• 213.000

0.06

• 213.000

ex2 1 3 24/0/17 1.62

• 15.000

0.04

• 15.000

ex2 1 4 12/0/10 1.28

• 11.000

0.04

• 11.000

ex2 1 5 29/0/30 2.21

• 268.015

0.13

• 268.015

ex2 1 6 26/0/21 3.31

• 39.000

0.12

• 39.000

ex2 1 7 71/0/61 1,388.00 4,150.410 3.98 4,150.410 ex9 2 2 38/0/37 1,355.47 100.000 0.36 100.000 ex9 2 3 54/0/53 5,755.76 0.000 0.73 0.000 ex9 2 6 57/0/53 (-1.003)

• 1.000

0.78

• 1.000

ex5 2 5 442/0/429 (-23,563.500)

• 3,500.000 (-11,975.600)
• 3,500.000

ex5 3 3 563/0/550 (-16,872.900) 3.056 (-16,895.400) 3.056 Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 18 / 19

Comparison

SC-MINLP COUENNE reformulated var/int/cons time time instance

• riginal

(LB) UB (LB) UB ex14 2 1 122/0/124 21.81 0.000 0.18 0.000 ex14 2 2 56/0/57 12.86 0.000 0.10 0.000 ex14 2 6 164/0/166 42.56 0.000 1.03 0.000 ex14 2 7 277/0/280 128.70 0.000 0.63 0.000 ex2 1 1 12/0/8 2.84

• 17.000

0.16

• 17.000

ex2 1 2 14/0/10 1.75

• 213.000

0.06

• 213.000

ex2 1 3 24/0/17 1.62

• 15.000

0.04

• 15.000

ex2 1 4 12/0/10 1.28

• 11.000

0.04

• 11.000

ex2 1 5 29/0/30 2.21

• 268.015

0.13

• 268.015

ex2 1 6 26/0/21 3.31

• 39.000

0.12

• 39.000

ex2 1 7 71/0/61 1,388.00 4,150.410 3.98 4,150.410 ex9 2 2 38/0/37 1,355.47 100.000 0.36 100.000 ex9 2 3 54/0/53 5,755.76 0.000 0.73 0.000 ex9 2 6 57/0/53 (-1.003)

• 1.000

0.78

• 1.000

ex5 2 5 442/0/429 (-23,563.500)

• 3,500.000 (-11,975.600)
• 3,500.000

ex5 3 3 563/0/550 (-16,872.900) 3.056 (-16,895.400) 3.056 du-opt 242/18/230 13.52 3.556 38.04 3.556 du-opt5 239/15/227 17.56 8.073 37.96 8.073 fo7 338/42/437 (8.759) 22.518 (1.95) 22.833 m6 254/30/327 185.13 82.256 54.13 82.256 no7 ar2 1 394/41/551 (90.583) 127.774 (73.78) 111.141 no7 ar3 1 394/41/551 (81.539) 107.869 (42.19) 113.810 no7 ar4 1 394/41/551 (76.402) 104.534 (54.85) 98.518

• 7 2 338/42/437

(79.365) 124.324 (5.85) 133.988 stockcycle 578/480/195 251.54 119,948.676 84.35 119,948.676 Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 18 / 19

Conclusions

Flexible framework: solves MINLPs where the non-convexity in the

• bjective and constraint functions is manifested as the sum of

non-convex univariate functions;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 19 / 19

Conclusions

Flexible framework: solves MINLPs where the non-convexity in the

• bjective and constraint functions is manifested as the sum of

non-convex univariate functions; No need of additional information (Matlab);

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 19 / 19

Conclusions

Flexible framework: solves MINLPs where the non-convexity in the

• bjective and constraint functions is manifested as the sum of

non-convex univariate functions; No need of additional information (Matlab); No need of complicated structures (expression tree);

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 19 / 19

Conclusions

Flexible framework: solves MINLPs where the non-convexity in the

• bjective and constraint functions is manifested as the sum of

non-convex univariate functions; No need of additional information (Matlab); No need of complicated structures (expression tree); Use any available NLP/MINLP solver interfaced to AMPL;

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 19 / 19

Conclusions

Flexible framework: solves MINLPs where the non-convexity in the

• bjective and constraint functions is manifested as the sum of

non-convex univariate functions; No need of additional information (Matlab); No need of complicated structures (expression tree); Use any available NLP/MINLP solver interfaced to AMPL; Significant success, expecially on instances which naturally belongs to the specified class of problems.

Claudia D’Ambrosio (University of Bologna) Aussois, 7 January, 2010 19 / 19