LaGO - An object oriented
library for solving MINLPs†
Ivo Nowak‡ , Hernán Alperin, and Stefan Vigerske
Humboldt-Universität zu Berlin
Institut für Mathematik
Rudower Chaussee 25,
D-12489 Berlin, Germany
http://www-iam.mathematik.hu-berlin.de/˜eopt/

Abstract. The paper describes a software package called LaGO for
solving nonconvex mixed integer nonlinear programs (MINLPs). The
main component of LaGO is a convex relaxation which is used for generating solution candidates and computing lower bounds of the optimalvalue. The relaxation is generated by reformulating the given MINLP asa block-separable problem, and replacing nonconvex functions by convex
underestimators. Results on medium size MINLPs are presented.
Keywords. mixed integer nonlinear programming, convex relaxation,
heuristics, decomposition, software
AMS classifications. 90C22, 90C20, 90C27, 90C26, 90C59

1

Introduction

Several strategies for solving nonconvex MINLPs have been proposed [BP03].
Exact solution approaches include branch-and-reduce [TS02], branch-and-bound[ADFN98,SP99], interval analysis [VEH96] and outer approximation [KAGB01].
On the other side, heuristic solution algorithms explicitely addressing MINLP,
such as the multistart scatter search heuristic [ULP+ 02], appeared rather rarelyin the literature. More information on global optimization algorithms can befound in [Flo00,HP95,BSV+ 01].
LaGO (LagrangianGlobalOptimizer) is an object oriented library for solving
nonconvex mixed integer nonlinear programs (MINLPs) written in C++. The
basic component of the solver is a relaxation, which is used for generating solution
candidates and computing lower bounds of the optimal value. It is very importantthat the quality of the relaxation is good enough, in the sense that the amountof work to retrieve a good solution from the relaxation is not too high. For
improving a given relaxation, a series of operations, which can be performed†

The work was supported by the German Research Foundation (DFG) under grant
NO 421/2-1.
‡
ivo@mathematik.hu-berlin.de, corresponding author

relatively fast and take advantage of a partially separable structure of the given
optimization problem, is used.
The software can be used as a general purpose MINLP solver or as a toolbox
for developing specialized solution algorithms. The object oriented design ofthe library allows easy replacement of modules, such as linear algebra routinesand various solvers. Optimization problems can be either provided in AMPLformat [FGK93] or coded in C++ using LaGO ’s class MinlpProblem. The solver
supports two types of structural properties: partial separability, i.e. the Hessians
have almost block-diagonal structure, and sparse and low rank Hessians, i.e.
there exist fast methods for multiplying a vector with a Hessian.
The basic components of the current version of LaGO are: preprocessing,convex relaxation and relaxation-based heuristic, which are described in Sections
2, 3 and 4. A branch-and-cut algorithm is currently under development and will
be added in the future. In Section 5 we report preliminary results, and giveconclusions in Section 6.
Notation. Let J ⊂ {1, . . . , n} be an index set. We denote by |J| the numberof elements of J. A subvector xJ ∈ IR|J| of x ∈ IRn is defined by (xj )j∈J . Similarly, a function fJ (x) is defined by (fj (x))j∈J . The space of k times differentiable
functions f : IRn 7→ IRm is denoted by C k (IRn , IRm ).

2

Preprocessing

We assume that the given MINLP problem has the form:
min h0 (x)s.t. hI (x) ≤ 0
hE (x) = 0
x∈Y

(P)

where
Y = {x ∈ [x, x] | xj ∈ {xj , xj } for j ∈ B},nx, x ∈ IR , B ⊆ {1, . . . , n}, I ∪ E = {1, . . . , m} with I ∩ E = ∅, and hi ∈C 2 (IRn , IR) for i = 0, . . . , m. Note that MINLPs with piecewise twice differentiable functions and integrality constraints, xi ∈ [xi , xi ]∩ZZ, can be transformed
into the form (P) by introducing additional binary variables.
To facilitate the construction of a relaxation, problem (P) is transformed into
the following block-separable extended reformulation.

(Pext )

min cT x + c0
s.t. AI x + bI ≤ 0
AE x + b E = 0
k = 1, . . . , p
g k (xJk ) ≤ 0,x∈Y

where {J1 , . . . , Jp } is a partition of {1, . . . , n}, i.e.

Sp

k=1 Jk= {1, . . . , n} andJi ∩ Jk = ∅ for i 6= k, g k ∈ C 2 (IRnk , IRmk ), k = 1, . . . , p, with nk = |Jk |, c ∈ IRn ,b ∈ IRm , AI ∈ IR(|I|,n) and AE ∈ IR(|E|,n) . To obtain the block functions g k (xJk )used in (Pext ) we construct the sparsity graph Gs = (V, Es ) with the nodes
V = {1, . . . , n} and the edgesno∂ 2 hk (x)Es = (i, j) ∈ V 2
6= 0 for some k ∈ {0, . . . , m} and x ∈ [x, x] .
∂xi ∂xjSp
Let {J˜1 , . . . , J˜p } be a partition of V . We define the set of nodes of l=k+1 J˜lSpconnected to J˜k by Rk = {i ∈ l=k+1 J˜l | (i, j) ∈ Es , j ∈ Jk }, for k = 1, . . . , p.
The set Rk can be interpreted as the set of flows of a nonlinear network problem
connecting a component J˜k with components J˜l , where k < l. After introducing
new variables y k ∈ IR|Rk | for every set Rk , and using the copy-constraints xRk =y k , problem (P) is formulated as a block-separable program with respect to theblocks Jk = (J˜k , Rk ), with k = 1, . . . , p. Finally, nonlinear equality constraintshi (x) = 0 for i ∈ E are replaced by the two inequality constraints hi (x) ≤ 0 and−hi (x) ≤ 0, and block-separable constraintspX

hki (xJk ) ≤ 0

k=1

are replaced bypX

tik ≤ 0,hki (xJk ) ≤ tik ,k = 1, . . . , p.

k=1

The resulting program has the form (Pext ), with gik (xJk , tik ) = hki (xJk )−tik . For
the sake of simplicity we keep the notation n, m, B, I, E, Y as in (P) including
the new variables tik into the block xJk .
Furthermore, the type (linear, convex, concave) of all functions is determinedby evaluating the minimum and maximum eigenvalue of each Hessian at sample
points. Since the given MINLP is coded in AMPL [FGK93], all functions are
given in a black-box representation, i.e. there exist only procedures for evaluating
functions, gradients and Hessians. Therefore, the generation of the sparsity graphand the determination of function types is performed by sampling techniques.
Problem (P) assumes the existence of bounds xj and xj for all j = 1, . . . , n.
This is not the usual case in the instances collected from MINLPLib [BDM03]
for instance. Since we need the box [x, x] to compute convex relaxations andfor the above sampling technique, we generate missing bounds by determining
the maximum and minimum value of a variable over a domain defined by theconvex constraints of problem (P). If there are not enough convex constraints todetermine a variable bound, we use a simple guessing algorithm.

3

Convex relaxation

A convex relaxation of problem (Pext ) is defined bymin cT x + c0
s.t. AI x + bI ≤ 0
AE x + b E = 0ğ k (xJk ) ≤ 0,
k = 1, . . . , p
x ∈ [x, x](C)

where ğ k are convex underestimators of g k over Xk = [xJk , xJk ], i.e. ğ k (x) ≤g k (x) for all x ∈ Xk and ğ k is convex over Xk . The best convex underestimators are convex envelopes. Since computing convex envelopes of nonconvex
functions can be computationally very expensive, we use the following procedurefor generating convex relaxations.
We construct at first a polynomial relaxation of (Pext ) by replacing nonquadratic functions by polynomial underestimators. We assume that the nonquadratic functions of (Pext ) are of the form g(x) = bT x + f (xN ) where thenumber |N | of nonlinear variables is small. We use as an underestimator for f amultivariate polynomial defined byp(x) = aT ϕ(x)(1)

where ϕ(x) = (xβ1 , . . . , xβr )T , βj ∈ IN0n , for 1 ≤ j ≤ r, is a vector of monomials,
and a ∈ IRr is the vector of coefficients for each monomial.
P The degree of thepolynomial p is the number d = maxrj=1 |βj | with |β| = ni=1 βi . For our published results, we use d = 2, but in our implementation larger degree polynomial
can be used. Let Dk be a differential operator defined by k∂ f (x)k,
D (f (x)) =∂xβl 1≤l≤rk
with |βl | = k for l = 1, . . . , rk ≤ kr can be at most all the possible monomialsof degree k, but in fact the sparsity pattern is used and fewer monomial thanrr
k are needed. In order to determine the coefficients a ∈ IR of the polynomial
underestimator (1), we solve the linear program

(U)

minr

a∈IR

s.t.

2
X

k=0

δk

X

kDk (f (x) − aT ϕ(x))k1x∈SkaT ϕ(x) ≤ f (x), x ∈ S0 .

The coefficients δ0 > δ1 , δ2 give the relative importance to the information thatcomes from the evaluation of the function, the gradient and the Hessian respectively. The finite sample sets Sk , k = 0, 1, 2 contain sample points for thecomputation of the function, gradient and Hessian respectively. In general, thisapproach is only rigorous for certain kind of functions such as convex and concave.
The polynomial relaxation of (Pext ) can be reformulated by adding someextra variables and constraints as the following mixed-integer quadratic program(MIQQP)(Q)min cT x + c0
s.t. AI x + bI ≤ 0
AE x + b E = 0
q k (xJk ) ≤ 0,
k = 1, . . . , px∈Y

where q k , k = 1, . . . , p, are quadratic forms. We use two methods for convexifying(Q).
The first method is based on replacing all nonconvex quadratic forms by
the so-called α-underestimators introduced by Adjiman and Floudas [AF97]. An
α-underestimator of a function f ∈ C 2 (IRn , IR) is the function
˘f(x)= f (x) + αT r(x)where
r(x) = Diag(x − x)(x − x)(2)nand Diag(·) denotes a diagonal matrix. The parameter α ∈ IR is computed
according to1
α = max{0, −ρ} Diag(w)−2 e2nwhere e ∈ IR is the vector of ones, w = x − x is the diameter vector of theinterval, andρ ≤ λ1 (Diag(w)∇2 f (x) Diag(w)), x ∈ [x, x],
is a bound on the the minimum eigenvalue of the transformed Hessian of f
over [x, x]. If f is a quadratic form, i.e. ∇2 f is constant, ρ can be determinedby eigenvalue computation. In the general case it can be computed by interval
Hessians [AF97] or using a sampling technique. It is clear that f˘(x) ≤ f (x) for
all x ∈ [x, x], and f˘ is convex. Note that f = f˘ if f is convex.
Applying the α-underestimator technique directly to the original functions
would also give a convex relaxation. However, our method is often tighter because
the α-convexification depends only on the curvature of the function and noton the function behaviour. For more clarification see the example in Figure 1
where f is the original function, f˘ the α-convexification of f , q the polynomial
underestimator, and q̆ the α-convexification of q.
The second method for convexifying problem (Q) is based on replacing thebox-constraints xC ∈ [xC , xC ] and the binary constraints xB ∈ {xB , xB } by
the quadratic constraints rC (x) ≤ 0, and rB (x) = 0 respectively, where thequadratic form r is defined as above, and C = {1, . . . , n} \ B is the index set ofcontinuous variables. The resulting reformulation of (Q) is all-quadratic, and its
dual is a semidefinite program. There exist several methods for solving semidefinite programs [WSV00]. We use a method based on formulating the dual as an
eigenvalue optimization problem [Now02]. This approach allows fast computation of near optimal dual solutions and supports decomposition. For each dual
point µ produced by this algorithm, the Lagrangian is convex, and defines thefollowing convex Lagrangian relaxation

f
q
q
f

Fig. 1. α-convex underestimator versus the convexification of the polynomial underestimator.

min L(x; µ),(Rµ )x∈[x,x]where L(·; µ) is the Lagrangian to the all-quadratic programming reformulationof (Q).
Besides of the above convexification methods, one further method for constructing a convex relaxation (C) is implemented. It uses quadratic convex global
underestimators proposed by Phillips, Rosen and Walke [PRW95]. Their method
for approximating the convex envelope of a nonlinear function f by a quadratic
function was improved in LaGO to reduce the absolute value of the smallest
eigenvalue of the obtained quadratic underestimator. Let S ⊂ [x, x] be a finite
sample set, and define the quadratic function
q(x; a, b, c) = c + 2bT x + xT Diag(a)x,where a, b ∈ IRn and c ∈ IR. Then q(·; a, b, c) is convex, if and only if a ≥ 0. The
tightest quadratic convex underestimator q(·; a, b, c) over the set S is provided
by the program
min
a,b,c

(CGU)

s.t.

X

f (x) −q(x; a, b, c) + δeT ax∈Sf (x) ≥ q(x; a, b, c),a ≥ 0.

for all x ∈ S

Since q depends linearly on a, b, c, problem (CGU) is a linear program. The term
δeT a reduces the absolute value of the smallest eigenvalue of Diag(a) in thecase where (CGU) is degenerated. The quality of these underestimators depends
strongly on the sample set S. If S contains all minimizers of f over [x, x], thebound is rigorous. Since we cannot guarantee to find all minimizers, this approachprovides a heuristic underestimator.

4

Relaxation-based heuristic

For computing solution candidates of a MINLP a rounding heuristic was developed, which is based on a convex relaxation of the given problem. The heuristic
works by computing a solution point of problem (C), and rounding some binary
components of this point. The rounded variables are fixed and the restrictedconvex relaxation is solved. This procedure is repeated as long, as either the
restricted convex relaxation is infeasible, or all binary variables are fixed. In thelatter case a local search is started from the last solution of the restricted convex relaxation giving a solution candidate. The values of the binary variables arerecursively switched and the whole process is repeated as long, as either all combinations of binary variables are searched, or the number of solution candidates
exceeds a given number smax .
RoundHeu:
Z = ∅;
FixVariables(0, ∅);FixVariables(y, J):ifΩ̆ ∩ Uy,J = ∅, then return;if
J = B,
then Optimize(y);else begin
get a solution x of (Cy,B );get K ⊂ B \ J where min{xj − xj , xj − xj } is small for j ∈ K;
for j ∈ K, round xj ;repeatcall FixVariables(x, J ∪ K);
switch some components of xK to get a new binary combination of xK ;until (|Z| ≥ smax or all combinations of xK are searched);
end
Optimize(y):if|B| = n,then if y ∈ Ω, then add y to Z;
else get a solution candidate of (Py,B ) by using a local search starting from thesolution of (Cy,B ), and add the point to Z;
Fig. 2. Algorithm rounding heuristic

To describe this method more detailed, let us define the original problem andthe relaxed problem with fixed binary variables

(Py.J )min{cT x | x ∈ Ω ∩ Uy,J }and(Cy,J )

min{cT x | x ∈ Ω̆ ∩ Uy,J }where
Uy,J = {x ∈ IRn | xJ = yJ },y ∈ [x, x], J ⊂ {1, . . . , n}, and Ω and Ω̆ are the feasible sets of problems (Pext )and (C) respectively. The recursive algorithm described in Fig. 2 computes a set
Z of binary feasible solution candidates of problem (Pext ).
In order to improve the simple method Optimize for solving the continuous
nonlinear program (Py,B ), a neighbourhood search or a deformation heuristic,
as described in [AN02], can be used.

5

Preliminary numerical results

We tested the described algorithm using a set of instances of MINLP that wereobtained from the GAMS web-site [BDM03,CG02]. For the computation of theminimum eigenvalue and corresponding eigenvector we used the Lanczos method
ARPACK++ [GS97]. The sequential quadratic programming code SNOPT[GMS97] is used for finding local solutions of nonlinear optimization problems.
To analyse its performance, the code was run on a machine with a 1GHz Pentium
III processor and 256 MB RAM. Table 2 reports the performance on a selected
set of problems using Algorithm RoundHeu. The definition of each column is
explained in Table 1. The limit on |Sf |, the maximum number of feasible points
explored was set to 50.

6Conclusions

We described the C++ library LaGO for solving general nonconvex MINLPs.
Main features of LaGO are: (i) the objective and constraint functions can begiven in a black-box formulation, i.e. it is only assumed that functions, gradients and Hessians can be evaluated; (ii) the given MINLP is automaticallyreformulated as a block-separable formulation giving the possibility to apply decomposition techniques; (iii) the object oriented design of the software allows
easy extensions and modifications of the current code.
Preliminary results with a relaxation-based rounding heuristic were presented, demonstrating that the current version of LaGO is able to solve medium
size MINLP instances in a reasonable time. The only two instances that were
not solved to known global optimality were those were the limit of |Sf | ≤ 50
feasible solutions was reached.

Table 1. Descriptions of the columns of Table 2.
rel.
error
heu.
time
k
|Sf |
n
m
|E|
|C|
|B|
|U |

|fheu − fbest |/(1 + |fbest |), where fheu , and fbest are the objective
value obtained with the heuristic and the best known objective value
for the problem respectively.
seconds required by Algorithm RoundHeu (without preprocessing) toobtain the heuristic solution xheu , and its corresponding objective
value fheu .
k ≤ 2|B| number of fixed binary combinations explored|Sf | ≤ k, number of feasible solutions obtained from each (Cy,B )description of the MINLP instances.
number of variables,
number of constraints,
number of equality constraints,
number of continuous variables,
number of binary variables,
number of unbounded variables(U = {i | xi = −∞ or xi = +∞}).

The performance of the described solution approach depends mainly on thequality of the relaxation. In order to improve the relaxation, and to search systematically for a global solution, a branch-and-cut algorithm with box-reduction
is currently under development and will be added in the future.

References[ADFN98] C. S. Adjiman, S. Dallwig, C. A. Floudas, and A. Neumaier. A global optimization method, αBB, for general twice-differentiable constrained NLPs
— I. Theoretical advances. Computers and Chemical Engineering, pages(1137–1158), 1998.
[AF97]
C. S. Adjiman and C. A. Floudas. Rigorous convex underestimators forgeneral twice-differentiable problems. J. of Global Opt., 9:23–40, 1997.
[AN02]
H. Alperin and I. Nowak. Lagrangian Smoothing Heuristics for MaxCut.
Technical report, HU–Berlin NR–2002–6, 2002.
[BDM03] M.R. Bussieck, A.S. Drud, and A. Meeraus. MINLPLib - A Collectionof Test Models for Mixed-Iinteger Nonlinear Programming. INFORMS J.
Comput., 15(1), 2003.
[BP03]
M.R. Bussieck and A. Pruessner. Mixed-integer nonlinear programming.
SIAG/OPT Newsletter: Views & News., 2003.
[BSV+ 01] C. Bliek, P. Spellucci, L. Vicente, A. Neumaier, L. Granvilliers, E. Monfroy, F. Benhamou, E. Huens, P. Van Hentenryck, D. Sam-Haroud, andB. Faltings. COCONUT Deliverable D1, Algorithms for Solving Nonlinear Constrained and Optimization Problems: The State of The Art.
http://www.mat.univie.ac.at/~neum/glopt.html, 2001.
[CG02]GAMS Development Corp. and GAMS Software GmbH. MINLP World,
2002. http://www.gamsworld.org/minlp/.

Table 2. Algorithm performance on a selected set of described examples.

example

rel.
error

heu.
time

alan
batch
batchdesex1221
ex1222
ex1223bex1224ex1225ex1226ex1252gbd
sep1
synthes1
synthes3
elf
ex3
ex4
fuel
gkocismeanvarx
oaer
procselsynheatsynthes2

00.02
0 57.13
0
0.70
0
0.04
0
0.00
0
0.11
0 13.32
0
0.07
0
0.03
.03 365.73
0
0.01
0
0.11
0
0.10
0
2.04
.40 52.22
01.53
0 31.69
0
0.56
0
0.08
0
0.51
0
0.05
0
0.04
0 15.94
00.26

24 examples,

22 solved[FGK93]

k |Sf |
6
84
12
4
1
16
81
7
6
121
3
3
5
20
118
4
22
4
6
20
6
2
21
12

6
50
8
4
1
16
45
6
5
50
3
3
5
20
50
4
20
4
6
20
6
2
20
12

n m

|E| |C| |B| |U |

9 8
47 74
20 20
6 6
4 4
8 10
12 8
9 11
6 6
40 44
5 5
30 32
7 7
18 24
55 39
33 32
37 31
16 16
12 9
36 45
10 8
11 8
57 65
12 15

3 5 4 13 23 24 7 11 9 3 3 3 1 3 1 1 4 4 3 4 8 3 3 6 2 3 3 23 25 15 1 2 3 23 28 2 1 4 3 3 10 8 7 31 24 30
18 25 8 15
1 12 25 7
7 13 3 6
6 9 3 6
9 22 14 21
4 7 3 65 8 3 621 45 12 32
2 7 5 2Robert Fourer, David M. Gay, and Brian W. Kernighan. AMPL: A Modeling Language for Mathematical Programming. Duxbury Press, Brooks/Cole
Publishing Company, 1993.
[Flo00]C. A. Floudas. Deterministic Global Optimization: Theory, Algorithms andApplications. Kluwer Academic Publishers, 2000.
[GMS97] P. E. Gill, W. Murray, and M. A. Saunders. SNOPT 5.3 user’s guide. Technical report, University of California, San Diego, Mathematics Department
Report NA 97-4, 1997.
[GS97]
F. Gomes and D. Sorensen. ARPACK++: a C++ Implementation of ARPACK
eigenvalue package, 1997. http://www.crpc.rice.edu/software/ARPACK/.
[HP95]
R. Horst and P. Pardalos. Handbook of Global Optimization. Kluwer Academic Publishers, 1995.
[KAGB01] P. Kesavan, R. J. Allgor, E. P. Gatzke, and P. I. Barton. Outer Approximation Algorithms for Separable Nonconvex Mixed-Integer Nonlinear Programs. submitted to Mathematical Programming, 2001.
[Now02] I. Nowak. Lagrangian Decomposition of Mixed-Integer All-Quadratic Programs. Technical report, HU–Berlin NR–2002–7, 2002.

[PRW95] A. Phillips, J. Rosen, and V. Walke. Molecular structure determination
by global optimization. In Dimacs Series in Discrete Mathematics andTheoretical Computer Science, volume 23, pages 181–198. 1995.
[SP99]
E. M. B. Smith and C. C. Pantelides. A Symbolic Reformulation/Spatial
Branch and Bound Algorithm for the Global Optimization of nonconvex
MINLPs. Computers and Chemical Engineering, 23:457–478, 1999.
[TS02]
M. Tawarmalani and N.V. Sahinidis. Convexification and Global Optimization in Continuous and Mixed-Integer Nonlinear Programming: Theory, Algorithms, Software, and Applications. Kluwer Academic Publishers, 2002.
[ULP+ 02] Zsolt Ugray, Leon Lasdon, John Plummer, Fred Glover, Jim Kelly, andRafael Marti. A multistart scatter search heuristic for smooth NLP
and MINLP problems. http://www.utexas.edu/courses/lasdon/papers.htm,
2002.
[VEH96] R. Vaidyanathan and M. EL-Halwagi. Global optimization of nonconvex
MINLP’s by interval analysis. In I. E. Grossmann, editor, Global Optimization in Engineering Design, pages 175–193. Kluwer Academic Publishers,
1996.
[WSV00] H. Wolkowicz, R. Saigal, and L. Vandenberghe. Handbook of Semidefinite
Programming. Kluwer Academic Publishers, 2000.