- Email: [email protected]
Piecewise polyhedral formulations for a multilinear term
Journal Pre-proof Piecewise polyhedral formulations for a multilinear term Kaarthik Sundar, Harsha Nagarajan, Jeff Linderoth, Site Wang, Russell Bent
PII: DOI: Reference:
S0167-6377(20)30190-5 https://doi.org/10.1016/j.orl.2020.12.002 OPERES 6677
To appear in:
Operations Research Letters
Received date : 6 August 2020 Revised date : 2 December 2020 Accepted date : 3 December 2020 Please cite this article as: K. Sundar, H. Nagarajan, J. Linderoth et al., Piecewise polyhedral formulations for a multilinear term, Operations Research Letters (2020), doi: https://doi.org/10.1016/j.orl.2020.12.002. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2020 Published by Elsevier B.V.
Journal Pre-proof
Piecewise Polyhedral Formulations for a Multilinear Term
b Department
lP repro of
Kaarthik Sundara , Harsha Nagarajana,˚, Jeff Linderothb , Site Wangc , Russell Benta a Los Alamos National Laboratory, Los Alamos, NM, USA of Industrial and Systems Engineering, Department of Computer Sciences, University of Wisconsin-Madison, Madison, WI, USA c Department of Industrial Engineering, Clemson University, SC, USA
Abstract
In this paper, we present a mixed-integer linear programming (MILP) formulation of a piecewise, polyhedral relaxation (PPR) of a multilinear term using its convex-hull representation. Based on the PPR’s solution, we also present a MILP formulation whose solutions are feasible for nonconvex, multilinear equations. We then present computational results showing the effectiveness of proposed formulations on standard benchmark nonlinear programs (NLPs) with multilinear terms and compare with a traditional formulation that is built using recursive bilinear groupings of multilinear terms. Keywords: multilinear, piecewise polyhedral relaxations, convex hull, global optimization 1. Introduction
Jou
rna
In the global optimization literature, many solution techniques used to solve nonconvex Mixed-Integer Nonlinear Programs (MINLPs) or Nonlinear Programs (NLPs) to optimality rely on convex relaxations [1] and piecewise convex relaxations of nonconvex functions [2, 3, 4, 5, 6] to obtain tight bounds on the optimal objective values of MINLPs/NLPs. The three primary topics which this letter focuses are as follows. First, on developing piecewise, polyhedral relaxations (PPR) of a multilinear term using the convex hull representation and second, on comparing this relaxation with traditional, recursive bilinear piecewise relaxations. Our work generalizes formulations that approximate functions with piecewise, linear functions [7, 8, 9, 10] to formulations that relax functions with piecewise, convex polyhedra. Third, since recovering feasible solutions for MINLPs/NLPs is of great interest for quantifying the quality of the solutions of PPR, this letter focuses on developing a MILP formulation which requires solutions to satisfy multilinear terms, thereby producing a locally optimal solution that is feasible for all nonconvex, multilinear equations. Throughout the rest of the letter, boldface is used to denote vectors and we formally define a multilinear term as φpxq : r`, us Ñ R, where ź φpxq “ xi . (1) iPI
Here, I is an index set for the set of variables, x, and r`, us “ tx P R|I| : ` ď x ď uu. For every i P I, we associate a spatial disjunction with variable xi defined by discretization points Si “ tsi,1 , si,2 , . . . , si,ni u ˚ Corresponding
author Email addresses: [email protected] (Kaarthik Sundar), [email protected] (Harsha Nagarajan), [email protected] (Jeff Linderoth), [email protected] (Site Wang), [email protected] (Russell Bent) Preprint submitted to Elsevier
December 3, 2020
Journal Pre-proof
lP repro of
and associated intervals Pi “ tδi1 , δi2 , ¨ ¨ ¨ , δini ´1 u where the intervals form a partition of the domain of xi given by r`i , ui s, i.e. δik “ rsi,k , si,k`1 s and `i “ si,1 ă si,2 ă ¨ ¨ ¨ ă si,ni “ ui . Given this notation, each xi , i P I is constrained to satisfy xi P rsi,1 , si,2 s _ rsi,2 , si,3 s _ ¨ ¨ ¨ _ rsi,ni ´1 , si,ni s.
(2)
Given Pi and Si for every i P I, we let pni ´ 1q and ni be the cardinality of the respective sets. For any i P I and k P t1, ¨ ¨ ¨ , ni ´ 1u, the interval δik is defined to be active when the value of xi lies in the interval δik . Eq. (2) then enforces one interval to be active per variable. Finally, the full set of discretization points and the partition sets defined by these points and intervals in the space r`, us are Ś Ś given by the sets S “ iPI Si and P “ iPI Pi , respectively. Given the multilinear term φpxq and the associated sets S and P, this work presents piecewise, polyhedral relaxations (PPR) for the graph of φpxq, given by X “ tpx, wq P r`, us ˆ R : w “ φpxqu.
(3)
An important special case occurs when φpxq is bilinear with a single partition on each variable, i.e., |I| “ 2 and |P| “ 1. Here, the McCormick relaxation, given by the following equations w ě u2 x1 ` u1 x2 ´ u1 u2 ,
w ě `2 x1 ` `1 x2 ´ `1 `2 ,
w ď u2 x1 ` `1 x2 ´ `1 u2 ,
w ď `2 x1 ` u1 x2 ´ u1 `2 ,
(4a) (4b)
rna
defines the convex hull [11, 12] of the set X, i.e., convpXq. For general multilinear terms i.e., |I| ě 3 with |P| “ 1, McCormick described a procedure that recursively applies the McCormick relaxation to products of variables. We refer to this procedure as the “recursive McCormick relaxation”. On general multilinear terms with asymmetric bounds on the variables, despite not capturing the convex hull [13], the recursive McCormick relaxation for a multilinear term is the basis for relaxations used in many global optimization solvers such as BARON [14], SCIP [15] and Couenne [16]. In contrast, a relaxation based on a vertex representation, for |I| ě 2 and |P| “ 1 does characterize the convex hull of X [17, 18] and is given by |I| |I| 2ÿ 2ÿ ! ) ˆs, w “ ˆsq , convpXq “ Proj px, w, λq P r`, us ˆ R ˆ ∆2|I| : x “ λs x λs φpx
x,w
s“1
(5)
s“1
Jou
ˆ1, x ˆ2, . . . , x ˆ 2|I| is the collection of all the vertices of the hyper-rectangle r`, us given by S and where x ∆2|I| is the 2|I| -dimensional 0-1 simplex. Note that for the case where |I| “ 2 and |P| “ 1 (bilinear term without any partitions), the set of feasible points are identical for formulations in Eq. (5) and Eq. (4), and thus are equivalent. When |I| ě 2 and |P| ě 2, the resulting MILP relaxation is expressed as a disjunctive union of |P| polytopes where each polytope is obtained by applying Eq. (5) to the domain defined by the corresponding partition set. The convex hull of this disjunctive union is completely characterized by utilizing the theory of disjunctive programming [19, 20] and the introduction of a sufficient number of auxiliary variables; in particular, by utilizing the formulations in [20]. The formulation in [20] is an extended formulation that generates constraints for each partition separately and then aggregating them. In contrast, the approach presented in this letter works directly with the combinatorial structure underlying the shared extreme points and presents non-extended PPR for a multilinear term. Furthermore, logarithmically-sized MILP 2
Journal Pre-proof
lP repro of
relaxations for disjunctive constraints in the special case of a bilinear term with spatial disjunction on one variable is addressed in [21, 22], and it’s generalization to multilinear terms in [23], albeit without much computational studies. To the best of our knowledge, the work in [21, 22, 23] is the state-of-the-art for PPR of a multilinear term. The contributions of this letter include (a) a special-order-set (SOS)-2-based PPR that uses a vertex representation for a single multilinear term and a spatial disjunction on each variable, and (b) given the PPR’s solution, a piecewise formulation that uses the vertex representation of the convex hull polytope to recover the feasible solution for the nonconvex multilinear equation. This PPR and the feasible solutions obtained are compared with a relaxation and feasible solutions recovered based on recursively relaxing bilinear groupings of the multilinear term akin to the recursive McCormick relaxation. We refer to this relaxation as the recursive piecewise polyhedral relaxation (R-PPR). To the best of our knowledge, this is the first work that presents a systematic comparison between an SOS-2-based PPR and a recursive piecewise relaxation, along with the feasible solution recovery for a general multilinear term with spatial disjunctions on every variable. 2. Piecewise Polyhedral Relaxation for a Multilinear Term
C
ź
xi
iPI
GPPR
rna
In this section, we present the PPR for X given I (the index set of variables), P (the partition set), and S (the set of discretization points). For notational simplicity, for any i P I and k P t1, . . . , ni ´ 1u, we use δ ki “ si,k and δsik “ si,k`1 to denote the lower and upper limits for this interval, respectively. ˆ s P S. We also define yi,k as a binary We associate a non-negative multiplier variable, λs , with each x indicator for the k th interval of variable xi , i P I. This variable takes a value of 1 when interval δik is active for variable xi and takes a value of 0 otherwise. We let λ and y denote the vector of continuous multiplier variables and indicator variables, respectively. Next, for each i P I, we define a set-valued ˆ s P S : eTi x ˆ s “ ru where ei is a unit vector whose ith coordinate is equal to 1 and function, µi prq “ tx is 0 everywhere else. This function defines the subset of points in S whose ith component is equal to r. ˆ “ř ˆ Finally, we let λpSq ˆ s PSˆ λs for any S Ď S. x Given these notations, the PPR for a multilinear term is given by ! ) ř “ px, w, λ, yq P r`, us ˆ R ˆ ∆|S| ˆ t0, 1u iPI pni ´1q | px, w, λ, yq satisfies Eq. (7)
(6)
where, Eq. (7) is given by ÿ ÿ ˆs, w “ ˆ s q, x“ λs x λs φpx k“1
ˆ s PS x
yi,k “ 1,
@i P I,
λpµi pδ ki qq ď yi,k´1 ` yi,k ,
λpµi pδ ni i ´1 qq
ď yi,ni ´1 ,
ˆ s P S, and, @x
yi,k P t0, 1u,
(7b)
@i P I,
λpµi pδ 1i qq ď yi,1 ,
λs ě 0,
(7a)
Jou
nÿ i ´1
ˆ s PS x
(7c)
@i P I, k P 2, . . . , ni ´ 2,
@i P I,
(7d) (7e) (7f) (7g)
@i P I, k P 1, 2, . . . , ni ´ 1. 3
Journal Pre-proof
x2 s2,4 x ˆ 13 ˆ9 x
s2,2
ˆ5 x
s2,1
ˆ1 x
ˆ 10 x
ˆ 11 x
ˆ6 x
ˆ7 x
ˆ 16 x ˆ 12 x
lP repro of
s2,3
ˆ 14 x ˆ 15 x
ˆ2 x
ˆ3 x
ˆ8 x ˆ4 x
s1,1 s1,2 s1,3 s1,4
x1
Figure 1: φpxq “ x1 x2 with |P1 | “ |P2 | “ 3
rna
Constraints in Eq. (7) use the binary variables of each variable partition to activate and deactivate the non-negative multiplier variables λs , s P S. In particular, Eq. (7b) ensures that only one partition per variable is active. The constraints in Eq. (7c) – (7e) enforce adjacency conditions on the λs variables akin to SOS-2 constraints [24, 25]. In the context of mathematical models for formulating general disjunctive union of polytopes, this model is also referred to as the “Convex-Combination (CC)" model or the “lambda" formulation in the literature [8]. The constraints in Eq. (7c) – (7e) ensure that the relaxation is the convex combination of the vertices of the active partition in P. The traditional SOS-2 is an ordered set of non-negative variables, of which at most two can be non-zero and they must be consecutive in their ordering. In this case, instead of an ordered set of non-negative variables, we have an ordered set that is a sum of non-negative variables i.e., for any i P I, λpµi pδ ki qqq, with k “ t1, . . . , ni ´ 1u, forms an ordered set of non-negative variable sums where SOS-2 constraints are imposed. It is not difficult to establish that the PPR in Eq. (6) has the property of being locally sharp, i.e., the projection of the LP relaxation of Eq. (6) to px, wq-space is exactly the convpXq. This is a property of the SOS-2 formulation for piecewise, linear approximations of functions [26]. In the literature, there are other vertex-based formulations which use this variable space and other formulations which use a logarithmic number of binary variables [26, 10, 27], however, we focus on providing a detailed comparison study between the SOS-2-based PPR in (7) and its recursive counterparts (in Sec. 3).
x“
16 ÿ
Jou
PPR Example. The PPR for a bilinear term with |P| “ 9 (3 partitions on each variable), shown in Fig. 1, is given in Eq. (8).
s“1
ˆs, w “ λs x
16 ÿ
s“1
ˆ s q, λs φpx
16 ÿ
s“1
(8a)
λs “ 1,
y1,1 ` y1,2 ` y1,3 “ 1, y2,1 ` y2,2 ` y2,3 “ 1,
(8b)
λ2 ` λ6 ` λ10 ` λ14 ď y1,1 ` y1,2 ,
λ5 ` λ6 ` λ7 ` λ8 ď y2,1 ` y2,2 ,
(8d)
λ4 ` λ8 ` λ12 ` λ16 ď y1,3 ,
λ13 ` λ14 ` λ15 ` λ16 ď y2,3 ,
(8f)
λ1 ` λ5 ` λ9 ` λ13 ď y1,1 ,
λ3 ` λ7 ` λ11 ` λ15 ď y1,2 ` y1,3 , λs ě 0 @s “ 1, . . . , 16
and,
λ1 ` λ2 ` λ3 ` λ4 ď y2,1 ,
λ9 ` λ10 ` λ11 ` λ12 ď y2,2 ` y2,3 ,
yi,k P t0, 1u
@i “ 1, 2, k “ 1, 2, 3. 4
(8c) (8e) (8g)
Journal Pre-proof
3. Recursive Piecewise Polyhedral Relaxation for a Multilinear Term with |I| ě 3
C
ź iPI
xi
GR-PPR
lP repro of
For a general multilinear term with |I| ě 3 and no partitions i.e., |P| “ 1, the recursive McCormick relaxation successively uses the McCormick relaxation in Eq. (4) on bilinear groupings of the terms. We extend this approach to build a R-PPR by recursively grouping bilinear terms and applying the PPR, given in Eq. (6), to each bilinear term. For a particular grouping of bilinear terms with |I| ě 3, the R-PPR is succinctly expressed as CB GPPR FPPR A EPPR PPR xx1 ¨ x2 y “ ... ¨ x|I|´1 ¨ x|I|
(9)
Notice that a different R-PPR can be created by choosing a different order of recursive grouping of variables, leading to different relaxation qualities [28, 29]. In computational studies, we present the performance of different groupings of bilinear terms in addition to the grouping presented in (9). Also, while applying R-PPR, only the original space of variables in the multilinear term are partitioned, and not the lifted/auxiliary variables resulting from recursive relaxations. 4. MILP-based Piecewise Formulation for Recovering Feasible Solutions
Jou
rna
Thus far, the focus of the paper has been on developing piecewise polyhedral relaxations for multilinear terms. The optimal solution to this relaxation provides an active partition for each variable. We let P ˚ and S ˚ denote the optimal active partition and the extreme points of the active partition chosen by the the PPR of the multilinear term, respectively. Let r`˚ , u˚ s denote the variable bound vector defined by the active partition P ˚ . Applying the convex hull formulation given in Eq. (5) for x P r`˚ , u˚ s yields the same solution produced by the PPR. Additionally, enforcing the solution to lie on the one-dimensional faces of the polytope in Eq. (5) for variable bounds r`˚ , u˚ s yields a feasible solution to the multilinear term, which is straight-forward to observe. This provides a simple MILP-based technique to compute a feasible solution of the original nonlinear problem of interest. Hence, in this section, we present an alternate MILP that recovers a feasible solution by requiring the solution to lie on an edge of the polytope which represents the relaxation of the multilinear term’s active partition. The formulation presented below is easily generalized to the case of requiring a solution to lie on the edges of the convex hull relaxation for any partition, but we omit this for the sake of simplicity of exposition. To enforce these additional constraints, we introduce some additional notation. Let V and E denote the vertices and the edges (one-dimensional faces) of the polytope defined by Eq. (5) for x P r`˚ , u˚ s. We know that |V| “ 2|I| and |E| “ |I| ¨ 2|I|´1 . Then, for every v P V, we define a set function γpvq “ te P E : e is incident on vu. Additionally, we introduce a binary variable ze for e P E that takes a value 1 if the solution lies on the edge e. Given these notations, the piecewise formulation that enforces the solution to lie on one of the edges of the polytope defined by Eq. (5) for x P r`˚ , u˚ s is defined by ! ) px, w, λ, zq P r`˚ , u˚ s ˆ R ˆ ∆|V| ˆ t0, 1u|E| | px, w, λ, zq satisfies Eq. (11) (10) where, Eq. (11) is given by ÿ ÿ ˆs, w “ ˆ s q, x“ λs x λs φpx ˆ s PS ˚ x
(11a)
ˆ s PS ˚ x
5
Journal Pre-proof
ePE
ze “ 1,
λs ě 0,
λv ď
ÿ
ePγpvq
(11b)
ze @v P V,
ˆ s P S ˚ , and, @x
ze P t0, 1u @e P E,
(11c)
lP repro of
ÿ
(11d)
where Eqs. (11a) and (11b) ensure that any solution to the formulation (11) is a convex combination of the vertices that constitute a one-dimensional face of the polytope in Eq. (5) and that only one face is active. Note that the above described method is also applicable if the original problem is an MINLP, as recovering feasible solutions to MINLPs is by itself an arduous task using off-the-shelf solvers. In the next section, we present extensive computational results that compare PPR with R-PPR, and the quality of recovered feasible solutions. 5. Computational Results
In this section, we present three sets of computational results. The first two sets of results use the NLP instance shown in Eq. (12) and the third set of results uses 60 nonlinear optimization problems that contain only multilinear terms taken from [30]. F “ max px1 x2 x3 x4 ` x3 x4 x5 x6 ` x5 x6 x7 x8 q
subject to: 100x1 ´ x2 ´ x3 ` 833x4 ` 95x5 ` x6 ´ x7 ` 100x8 ď 50000,
100 ď x1 ď 500, 1000 ď x2 , x3 ď 2000, 10 ď xi ď 100 @i P t4, . . . , 8u.
(12a) (12b) (12c)
Jou
rna
The first set of results shows the strength of the various relaxations (PPR and R-PPRs) and the effectiveness of feasible solution recovery performed on the active partition obtained from the respective relaxations. The second set of results is aimed at comparing the different ways in which feasible solution recovery can be performed by using the template formulation in Sec. 4 in terms of computation time and quality of the feasible solutions. The final set of results compares the strength of the proposed PPRs on instances used in the literature for which the optimal solutions are known a-priori. For both the first and the third set of computational experiments, we partition the variables of the multilinear terms with uniformly located discretization points. All models were implemented in Julia using JuMP v0.19.5 [31] and the experimental results used Gurobi v8.0 [32] with Gurobi’s presolver and heuristics turned off. All nonconvex models were solved using BARON v19.12.7 with CPLEX 12.10 and Ipopt 3.12.8 as BARON’s sub-solvers. The computational experiments were performed on an Intel Xeon CPU L5420 @2.50GHz with 64 GB RAM. 5.1. Effectiveness of the formulations on the NLP in (12) Table 1 compares the objective value of the PPR with three R-PPRs (denoted by upper bound, U B), each with a different variable grouping, on the NLP in Eq. (12). Note that the choice of grouping changes the quality of the relaxation. The active partition of each variable in the relaxed solution is used to compute a feasible solution (whose objective is denoted by lower bound, LB) through formulation Eq. T T ´LB (11). In this table, U B and LB gaps are given by U B´OP ¨ 100 and OP LB ¨ 100, respectively; OP T UB corresponds to the global optimum value of (12) which is equal to 3.2642E10. It is clear from these results that, in both the relaxation and the feasible solution recovery formulations, PPR outperforms all the R-PPRs. For example, with only 4 partitions, PPR finds a relaxed solution (U B) that is better 6
Journal Pre-proof
# of partitions per variable
2
% gaps xa xb xc xd xxa xxb xc yyxd xxxa xb yxc yxd xa xxb xxc xd yy
LB, U B 2.33, 23.99 2.33, 65.47 2.33, 65.47 2.33, 47.37
lP repro of
than any R-PPR (up to 12 partitions). It is also clear that the quality of the feasible solutions (LB) recovered based on the active partition of the PPR formulation is consistently better than those of R-PPR’s for every chosen number of partitions. Note that, the LB gaps do not necessarily decrease monotonically as the number of partitions increase; the new partitions for PPR/R-PPR formulations yield different active partitions, thus leading to different local optimal solutions. In these experiments, the run-times of both PPR and R-PPR formulations were less than a minute and hence we do not report them explicitly. As is the case with any piecewise formulation, the computational performance of PPR degrades as the number of partitions increases on larger NLPs, however the strength of PPR’s relaxation remains superior compared to any other recursive variation. 4
6
8
10
12
LB, U B 0.15, 3.20 36.74, 25.37 36.74, 25.37 5.73, 25.27
LB, U B 1.11, 2.98 2.33, 21.73 2.33, 21.73 1.11, 16.78
LB, U B 0.15, 0.83 0.15, 14.72 0.15, 14.72 0.15, 12.29
LB, U B 0.00, 0.69 0.00, 12.98 0.00, 12.98 0.00, 9.59
LB, U B 0.05, 0.43 2.33, 10.34 2.33, 10.34 1.11, 8.19
Table 1: Relative optimality gaps of upper and lower bounds in percent for the PPR and the R-PPRs for varying number of partitions per variable. The x¨y represents the grouping of the different variables into bilinear terms in the R-PPRs.
rna
BARON’s performance on the NLP in (12): In order to compare the relative optimality gaps obtained by the PPR formulation in Table 1, we also solved instance (12) using BARON solver. The best optimality gap BARON reports after a run time of one hour is 20.36% with 153,555 open spatial branch-and-bound nodes, where the best upper bound reported by BARON is 4.0986 ˆ 1010 and the best lower bound by Ipopt within BARON is 3.2642 ˆ 1010 . In contrast, the PPR reports 3.34% with only four partitions per variable (see table 1). The best gap PPR reports is 0.47% (12 partitions) which is close to global optimum, thus reinforcing the value of the convex hull formulations proposed in this letter. 5.2. Feasible solution recovery for the NLP in (12)
Jou
Here, we discuss the quality of the feasible solutions obtained by using different variants of the MILP formulation presented in Sec. 4 for recovering solutions to the NLP described in (12). To this end, we introduce auxiliary variables t1 , t2 and t3 , such that the objective function of (12) will be to maximize pt1 ` t2 ` t3 q subject to constraints t1 “ x1 x2 x3 x4 , t2 “ x3 x4 x5 x6 , t3 “ x5 x6 x7 x8 , (12b) and (12c). The MILP formulation, as presented in Sec. 4, when applied to each multilinear term in Eq. (12), results in a feasible solution to the NLP. This formulation computes the best solution that lies on the one-dimensional face of the active partition of each partitioned variable as dictated by the solution of the relaxation. From here on, we refer to this formulation as F a . Instead, the formulation can be extended to compute the best solution that lies on the one-dimensional face of any combination of the variable partition (the formulation is a trivial extension of the formulation in Eq. (11), and is not presented). We denote this formulation as F1f . Though formulation F1f may sometimes lead to better quality feasible solutions, the number of one-dimensional faces grows exponentially with the number of partitions on each variable and thus can result in a MILP with a large number of binary variables. Note 7
Journal Pre-proof
# of partitions
xa xb xc xd
xxa xxb xc yyxd
xxxa xb yxc yxd
xa xxb xxc xd yy
Fa
F1f
Fa
F1f
F2f
Fa
F1f
F2f
Fa
F1f
F2f
2 4 6 8 10 12
2.33 0.15 1.11 0.15 0.0 0.05
2.33 0.15 1.11 0.15 0.0 0.05
86.25 0.15 4.22 0.15 1.12 0.15
2.33 36.74 2.33 0.15 0.0 2.33
0.17 0.05 0.03 5E-3 5E-4 6E-5
86.25 0.15 4.22 0.15 1.12 0.15
2.33 36.74 2.33 0.15 0.0 2.33
0.17 0.01 0.03 0.05 5E-3 6E-5
2.33 2.33 1.11 0.17 0.0 0.05
2.33 5.73 1.11 0.15 0.0 1.11
0.01 2E-4 4E-3 2E-5 4E-5 2E-5
lP repro of
per variable
Table 2: Quality (percent gaps) of feasible solutions obtained using different variants of formulation (11) on the NLP in Eq. (12).
rna
that while applying F1f on various recursive grouping of bilinear terms, we do not partition the auxiliary lifted variables. The third formulation, F2f , is applicable only when using recursive groupings of bilinear terms. Here, we enforce partitions on the auxiliary lifted variables introduced while performing the recursion. We also remark that comparing the solutions of F a and F1f with F2f is not a fair comparison. F2f naturally provides better quality solutions since more variables are partitioned. Nevertheless, we compare all the formulations in Table 2. In this table, all percent gaps are computed relative to the T ´LB global optimum value (OP T ), that is, OP LB ¨ 100. As expected and consistent with our previous observations, even in the feasible solution recovery process, the convex hull representation of the multilinear terms outperforms the recursive approaches when comparing formulations F a and F1f . As for formulation F2f , it outperforms even the convex hull formulation as it partitions the auxiliary lifted variables. As a side note, it is not difficult to prove that for a fixed partition count, F2f on any recursive bilinear grouping will always produce the best feasible solution among all the formulations i.e., F1f and F a . 5.3. Strength of the relaxations on instances from [30]
Jou
Our final set of results are based on a collection of 60 multilinear problems taken from [30], which are available at: https://minlp.com/downloads/testlibs/barlibs/mult3.zip. The bounds on all the variables for every problem instance in [30] are r0, 1s. To demonstrate the strength of the relaxations discussed in this paper, we created a modified set of instances. For each variable, we uniformly distributed the lower and upper bounds in the interval r0.1, 0.2s and r0.9, 1s, respectively. This modification is motivated by the following known result in [13]–the recursive McCormick relaxation for a multilinear term is the convex hull when all variable bounds are symmetric about the origin or are r0, 1s. Though the bounds were modified, the global optimum values for all the instances remained the same as those in [30] when solved using BARON. For each problem instance we relaxed the problems using PPR and a R-PPR with a lexicographic ordering. A time limit of one hour was imposed on every run of the instance. The results of all the runs are summarized using the box-plot shown in Fig. 2. In Fig. 2, the optimality gap for the PPR in most instances is 0% when the number of partitions per variable is three. This is not the case for R-PPR. However, as expected, Fig. 2(a) suggests that when the bound on every variable in the problem is r0, 1s, R-PPR is generally the better alternative. This is not surprising since similar recursive relaxations are known to capture the convex hull of a multilinear term with symmetric variable bounds [13]. Once the variable bounds are randomized 8
40
PPR
lP repro of
Journal Pre-proof
3500
R-PPR
Computation time (sec)
Optimality gap (%)
30
3000
20
PPR
R-PPR
2500
2000
1500
1000
10
500
0
0
2
3 4 Number of partitions
2
3 4 Number of partitions
(a) Percent gaps and run times for instances with variable bounds set to r0, 1s. The optimality gaps for the PPRs increase from 3 partitions to 4 partitions because of the computation time limit of 1 hour imposed on every run.
35
PPR
R-PPR
15 10 5
2
PPR
R-PPR
2500
2000
1500
1000
500
0
Jou
0
Computation time (sec)
25
3000
rna
Optimality gap (%)
30
20
3500
3 4 Number of partitions
2
3 4 Number of partitions
(b) Percent gaps and run times for instances with variable bounds randomly chosen as described earlier. In general, the R-PRR has better computation times because it has fewer multiplier variables (λ). Nevertheless, the PPR closes the optimality gap with fewer partitions despite taking a little more computation time. Figure 2: The optimality gap in % refers to the relative gap between the objective value of the relaxation of the problem and the global optimum to the problem instance. The computation time is the time taken to solve the PPR and the R-PPR to optimality.
9
Journal Pre-proof
6. Concluding Remarks
lP repro of
and made asymmetric, the PPR formulation outperforms the recursive relaxation (R-PPR) in both optimality gap and computation time (Fig. 2(b)). Though the (bottom right) figure suggests that the computation times of PPR are larger than that of R-PPR, it is worth noting that the number of partitions and the run times necessary for the PPR to achieve 0% optimality gap on all 60 instances is much superior when compared to the R-PPR.
This letter presents piecewise polyhedral relaxation formulations for a multilinear term with spatial disjunctions on every variable. An SOS-2-based formulation that generalizes similar formulations for piecewise linear approximations of functions was developed to formulate the convex hull of a single multilinear term. Extensive computational experiments demonstrated that our proposed formulations have significant computational advantages over the well-known recursive grouping of bilinear terms. A MILP-based piecewise formulation was developed to recover high-quality feasible solutions for problems with multilinear terms. Given exponential growth in the number of variables in the PPR formulation, future work will focus on developing formulations with subsets of variables for those rare cases when the number of variables in a multilinear term is large and further add valid inequalities, by generalizing the ideas proposed in [30]. Acknowledgements
7. References
rna
The work was funded by the Center for Nonlinear Studies (CNLS) at LANL and LANL’s Directed Research and Development (LDRD) projects, “20170201ER: POD: A Polyhedral Outer-approximation, Dynamic-discretization optimization solver" and “20190590ECR: Discrete Optimization Algorithms for Provably Optimal Quantum Circuit Design". This work was carried out under the U.S. DOE Contract No. DE-AC52-06NA25396.
[1] M. Tawarmalani, N. V. Sahinidis, A polyhedral branch-and-cut approach to global optimization, Mathematical Programming 103 (2) (2005) 225–249. [2] M. Hasan, I. Karimi, Piecewise linear relaxation of bilinear programs using bivariate partitioning, AIChE journal 56 (7) (2010) 1880–1893.
Jou
[3] L. S. de Assis, E. Camponogara, B. Zimberg, E. Ferreira, I. E. Grossmann, A piecewise McCormick relaxation-based strategy for scheduling operations in a crude oil terminal, Computers & Chemical Engineering 106 (2017) 309–321. [4] P. M. Castro, Tightening piecewise McCormick relaxations for bilinear problems, Computers & Chemical Engineering 72 (2015) 300–311. [5] H. Nagarajan, M. Lu, E. Yamangil, R. Bent, Tightening McCormick relaxations for nonlinear programs via dynamic multivariate partitioning, in: International Conference on Principles and Practice of Constraint Programming, Springer, 2016, pp. 369–387.
10
Journal Pre-proof
lP repro of
[6] H. Nagarajan, M. Lu, S. Wang, R. Bent, K. Sundar, An adaptive, multivariate partitioning algorithm for global optimization of nonconvex programs, Journal of Global Optimization 74 (4) (2019) 639–675. [7] S. Sridhar, J. Linderoth, J. Luedtke, Locally ideal formulations for piecewise linear functions with indicator variables, Operations Research Letters 41 (6) (2013) 627–632. [8] J. P. Vielma, S. Ahmed, G. Nemhauser, Mixed-integer models for nonseparable piecewise-linear optimization: Unifying framework and extensions, Operations research 58 (2) (2010) 303–315. [9] J. Lee, D. Wilson, Polyhedral methods for piecewise-linear functions i: the lambda method, Discrete applied mathematics 108 (3) (2001) 269–285. [10] M. Padberg, Approximating separable nonlinear functions via mixed zero-one programs, Operations Research Letters 27 (1) (2000) 1 – 5. [11] G. P. McCormick, Computability of global solutions to factorable nonconvex programs: Part I–convex underestimating problems, Mathematical programming 10 (1) (1976) 147–175. [12] F. A. Al-Khayyal, J. E. Falk, Jointly constrained biconvex programming, Mathematics of Operations Research 8 (2) (1983) 273–286. [13] J. Luedtke, M. Namazifar, J. Linderoth, Some results on the strength of relaxations of multilinear functions, Mathematical programming 136 (2) (2012) 325–351. [14] N. V. Sahinidis, Baron: A general purpose global optimization software package, Journal of global optimization 8 (2) (1996) 201–205. [15] S. Vigerske, A. Gleixner, SCIP: Global optimization of mixed-integer nonlinear programs in a branch-and-cut framework, Optimization Methods and Software 33 (3) (2018) 563–593.
rna
[16] P. Belotti, Couenne: a user’s manual, Tech. rep., Lehigh University (2009).
[17] C. A. Floudas, Deterministic global optimization: theory, methods and applications, Vol. 37, Springer Science & Business Media, 2013. [18] A. D. Rikun, A convex envelope formula for multilinear functions, Journal of Global Optimization 10 (4) (1997) 425–437.
Jou
[19] E. Balas, Disjunctive programming, Annals of discrete mathematics 5 (1979) 3–51. [20] R. G. Jeroslow, J. K. Lowe, Modelling with integer variables, Mathematical Programming at Oberwolfach II (1984) 167–184. [21] R. Misener, C. A. Floudas, Global optimization of mixed-integer quadratically-constrained quadratic programs (MIQCQP) through piecewise-linear and edge-concave relaxations, Mathematical Programming 136 (1) (2012) 155–182. [22] R. Misener, J. P. Thompson, C. A. Floudas, APOGEE: global optimization of standard, generalized, and extended pooling problems via linear and logarithmic partitioning schemes, Computers & Chemical Engineering 35 (5) (2011) 876–892. 11
Journal Pre-proof
[23] J. Huchette, J. P. Vielma, A combinatorial approach for small and strong formulations of disjunctive constraints, Mathematics of Operations Research 44 (3) (2019) 793–820.
lP repro of
[24] E. Beale, J. J. Forrest, Global optimization using special ordered sets, Mathematical Programming 10 (1) (1976) 52–69. [25] E. M. L. Beale, J. A. Tomlin, Special facilities in a general mathematical programming system for non-convex problems using ordered sets of variables, OR 69 (447-454) (1970) 99. [26] J. Huchette, J. P. Vielma, Nonconvex piecewise linear functions: Advanced formulations and simple modeling tools, arXiv preprint:1708.00050. [27] R. Misener, C. A. Floudas, Global optimization of mixed-integer quadratically-constrained quadratic programs (miqcqp) through piecewise-linear and edge-concave relaxations, Mathematical Programming 136 (1) (2012) 155–182. [28] P. Belotti, S. Cafieri, J. Lee, L. Liberti, A. J. Miller, On the composition of convex envelopes for quadrilinear terms, in: Optimization, Simulation, and Control, Springer, 2013, pp. 1–16. [29] J. Lee, D. Skipper, E. Speakman, Algorithmic and modeling insights via volumetric comparison of polyhedral relaxations, Mathematical Programming 170 (1) (2018) 121–140. [30] X. Bao, A. Khajavirad, N. V. Sahinidis, M. Tawarmalani, Global optimization of nonconvex problems with multilinear intermediates, Mathematical Programming Computation 7 (1) (2015) 1–37. [31] I. Dunning, J. Huchette, M. Lubin, JuMP: A modeling language for mathematical optimization, SIAM Review 59 (2) (2017) 295–320.
Jou
rna
[32] Gurobi Optimizer 7.5, http://www.gurobi.com.
12
PII: DOI: Reference:
S0167-6377(20)30190-5 https://doi.org/10.1016/j.orl.2020.12.002 OPERES 6677
To appear in:
Operations Research Letters
Received date : 6 August 2020 Revised date : 2 December 2020 Accepted date : 3 December 2020 Please cite this article as: K. Sundar, H. Nagarajan, J. Linderoth et al., Piecewise polyhedral formulations for a multilinear term, Operations Research Letters (2020), doi: https://doi.org/10.1016/j.orl.2020.12.002. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2020 Published by Elsevier B.V.
Journal Pre-proof
Piecewise Polyhedral Formulations for a Multilinear Term
b Department
lP repro of
Kaarthik Sundara , Harsha Nagarajana,˚, Jeff Linderothb , Site Wangc , Russell Benta a Los Alamos National Laboratory, Los Alamos, NM, USA of Industrial and Systems Engineering, Department of Computer Sciences, University of Wisconsin-Madison, Madison, WI, USA c Department of Industrial Engineering, Clemson University, SC, USA
Abstract
In this paper, we present a mixed-integer linear programming (MILP) formulation of a piecewise, polyhedral relaxation (PPR) of a multilinear term using its convex-hull representation. Based on the PPR’s solution, we also present a MILP formulation whose solutions are feasible for nonconvex, multilinear equations. We then present computational results showing the effectiveness of proposed formulations on standard benchmark nonlinear programs (NLPs) with multilinear terms and compare with a traditional formulation that is built using recursive bilinear groupings of multilinear terms. Keywords: multilinear, piecewise polyhedral relaxations, convex hull, global optimization 1. Introduction
Jou
rna
In the global optimization literature, many solution techniques used to solve nonconvex Mixed-Integer Nonlinear Programs (MINLPs) or Nonlinear Programs (NLPs) to optimality rely on convex relaxations [1] and piecewise convex relaxations of nonconvex functions [2, 3, 4, 5, 6] to obtain tight bounds on the optimal objective values of MINLPs/NLPs. The three primary topics which this letter focuses are as follows. First, on developing piecewise, polyhedral relaxations (PPR) of a multilinear term using the convex hull representation and second, on comparing this relaxation with traditional, recursive bilinear piecewise relaxations. Our work generalizes formulations that approximate functions with piecewise, linear functions [7, 8, 9, 10] to formulations that relax functions with piecewise, convex polyhedra. Third, since recovering feasible solutions for MINLPs/NLPs is of great interest for quantifying the quality of the solutions of PPR, this letter focuses on developing a MILP formulation which requires solutions to satisfy multilinear terms, thereby producing a locally optimal solution that is feasible for all nonconvex, multilinear equations. Throughout the rest of the letter, boldface is used to denote vectors and we formally define a multilinear term as φpxq : r`, us Ñ R, where ź φpxq “ xi . (1) iPI
Here, I is an index set for the set of variables, x, and r`, us “ tx P R|I| : ` ď x ď uu. For every i P I, we associate a spatial disjunction with variable xi defined by discretization points Si “ tsi,1 , si,2 , . . . , si,ni u ˚ Corresponding
author Email addresses: [email protected] (Kaarthik Sundar), [email protected] (Harsha Nagarajan), [email protected] (Jeff Linderoth), [email protected] (Site Wang), [email protected] (Russell Bent) Preprint submitted to Elsevier
December 3, 2020
Journal Pre-proof
lP repro of
and associated intervals Pi “ tδi1 , δi2 , ¨ ¨ ¨ , δini ´1 u where the intervals form a partition of the domain of xi given by r`i , ui s, i.e. δik “ rsi,k , si,k`1 s and `i “ si,1 ă si,2 ă ¨ ¨ ¨ ă si,ni “ ui . Given this notation, each xi , i P I is constrained to satisfy xi P rsi,1 , si,2 s _ rsi,2 , si,3 s _ ¨ ¨ ¨ _ rsi,ni ´1 , si,ni s.
(2)
Given Pi and Si for every i P I, we let pni ´ 1q and ni be the cardinality of the respective sets. For any i P I and k P t1, ¨ ¨ ¨ , ni ´ 1u, the interval δik is defined to be active when the value of xi lies in the interval δik . Eq. (2) then enforces one interval to be active per variable. Finally, the full set of discretization points and the partition sets defined by these points and intervals in the space r`, us are Ś Ś given by the sets S “ iPI Si and P “ iPI Pi , respectively. Given the multilinear term φpxq and the associated sets S and P, this work presents piecewise, polyhedral relaxations (PPR) for the graph of φpxq, given by X “ tpx, wq P r`, us ˆ R : w “ φpxqu.
(3)
An important special case occurs when φpxq is bilinear with a single partition on each variable, i.e., |I| “ 2 and |P| “ 1. Here, the McCormick relaxation, given by the following equations w ě u2 x1 ` u1 x2 ´ u1 u2 ,
w ě `2 x1 ` `1 x2 ´ `1 `2 ,
w ď u2 x1 ` `1 x2 ´ `1 u2 ,
w ď `2 x1 ` u1 x2 ´ u1 `2 ,
(4a) (4b)
rna
defines the convex hull [11, 12] of the set X, i.e., convpXq. For general multilinear terms i.e., |I| ě 3 with |P| “ 1, McCormick described a procedure that recursively applies the McCormick relaxation to products of variables. We refer to this procedure as the “recursive McCormick relaxation”. On general multilinear terms with asymmetric bounds on the variables, despite not capturing the convex hull [13], the recursive McCormick relaxation for a multilinear term is the basis for relaxations used in many global optimization solvers such as BARON [14], SCIP [15] and Couenne [16]. In contrast, a relaxation based on a vertex representation, for |I| ě 2 and |P| “ 1 does characterize the convex hull of X [17, 18] and is given by |I| |I| 2ÿ 2ÿ ! ) ˆs, w “ ˆsq , convpXq “ Proj px, w, λq P r`, us ˆ R ˆ ∆2|I| : x “ λs x λs φpx
x,w
s“1
(5)
s“1
Jou
ˆ1, x ˆ2, . . . , x ˆ 2|I| is the collection of all the vertices of the hyper-rectangle r`, us given by S and where x ∆2|I| is the 2|I| -dimensional 0-1 simplex. Note that for the case where |I| “ 2 and |P| “ 1 (bilinear term without any partitions), the set of feasible points are identical for formulations in Eq. (5) and Eq. (4), and thus are equivalent. When |I| ě 2 and |P| ě 2, the resulting MILP relaxation is expressed as a disjunctive union of |P| polytopes where each polytope is obtained by applying Eq. (5) to the domain defined by the corresponding partition set. The convex hull of this disjunctive union is completely characterized by utilizing the theory of disjunctive programming [19, 20] and the introduction of a sufficient number of auxiliary variables; in particular, by utilizing the formulations in [20]. The formulation in [20] is an extended formulation that generates constraints for each partition separately and then aggregating them. In contrast, the approach presented in this letter works directly with the combinatorial structure underlying the shared extreme points and presents non-extended PPR for a multilinear term. Furthermore, logarithmically-sized MILP 2
Journal Pre-proof
lP repro of
relaxations for disjunctive constraints in the special case of a bilinear term with spatial disjunction on one variable is addressed in [21, 22], and it’s generalization to multilinear terms in [23], albeit without much computational studies. To the best of our knowledge, the work in [21, 22, 23] is the state-of-the-art for PPR of a multilinear term. The contributions of this letter include (a) a special-order-set (SOS)-2-based PPR that uses a vertex representation for a single multilinear term and a spatial disjunction on each variable, and (b) given the PPR’s solution, a piecewise formulation that uses the vertex representation of the convex hull polytope to recover the feasible solution for the nonconvex multilinear equation. This PPR and the feasible solutions obtained are compared with a relaxation and feasible solutions recovered based on recursively relaxing bilinear groupings of the multilinear term akin to the recursive McCormick relaxation. We refer to this relaxation as the recursive piecewise polyhedral relaxation (R-PPR). To the best of our knowledge, this is the first work that presents a systematic comparison between an SOS-2-based PPR and a recursive piecewise relaxation, along with the feasible solution recovery for a general multilinear term with spatial disjunctions on every variable. 2. Piecewise Polyhedral Relaxation for a Multilinear Term
C
ź
xi
iPI
GPPR
rna
In this section, we present the PPR for X given I (the index set of variables), P (the partition set), and S (the set of discretization points). For notational simplicity, for any i P I and k P t1, . . . , ni ´ 1u, we use δ ki “ si,k and δsik “ si,k`1 to denote the lower and upper limits for this interval, respectively. ˆ s P S. We also define yi,k as a binary We associate a non-negative multiplier variable, λs , with each x indicator for the k th interval of variable xi , i P I. This variable takes a value of 1 when interval δik is active for variable xi and takes a value of 0 otherwise. We let λ and y denote the vector of continuous multiplier variables and indicator variables, respectively. Next, for each i P I, we define a set-valued ˆ s P S : eTi x ˆ s “ ru where ei is a unit vector whose ith coordinate is equal to 1 and function, µi prq “ tx is 0 everywhere else. This function defines the subset of points in S whose ith component is equal to r. ˆ “ř ˆ Finally, we let λpSq ˆ s PSˆ λs for any S Ď S. x Given these notations, the PPR for a multilinear term is given by ! ) ř “ px, w, λ, yq P r`, us ˆ R ˆ ∆|S| ˆ t0, 1u iPI pni ´1q | px, w, λ, yq satisfies Eq. (7)
(6)
where, Eq. (7) is given by ÿ ÿ ˆs, w “ ˆ s q, x“ λs x λs φpx k“1
ˆ s PS x
yi,k “ 1,
@i P I,
λpµi pδ ki qq ď yi,k´1 ` yi,k ,
λpµi pδ ni i ´1 qq
ď yi,ni ´1 ,
ˆ s P S, and, @x
yi,k P t0, 1u,
(7b)
@i P I,
λpµi pδ 1i qq ď yi,1 ,
λs ě 0,
(7a)
Jou
nÿ i ´1
ˆ s PS x
(7c)
@i P I, k P 2, . . . , ni ´ 2,
@i P I,
(7d) (7e) (7f) (7g)
@i P I, k P 1, 2, . . . , ni ´ 1. 3
Journal Pre-proof
x2 s2,4 x ˆ 13 ˆ9 x
s2,2
ˆ5 x
s2,1
ˆ1 x
ˆ 10 x
ˆ 11 x
ˆ6 x
ˆ7 x
ˆ 16 x ˆ 12 x
lP repro of
s2,3
ˆ 14 x ˆ 15 x
ˆ2 x
ˆ3 x
ˆ8 x ˆ4 x
s1,1 s1,2 s1,3 s1,4
x1
Figure 1: φpxq “ x1 x2 with |P1 | “ |P2 | “ 3
rna
Constraints in Eq. (7) use the binary variables of each variable partition to activate and deactivate the non-negative multiplier variables λs , s P S. In particular, Eq. (7b) ensures that only one partition per variable is active. The constraints in Eq. (7c) – (7e) enforce adjacency conditions on the λs variables akin to SOS-2 constraints [24, 25]. In the context of mathematical models for formulating general disjunctive union of polytopes, this model is also referred to as the “Convex-Combination (CC)" model or the “lambda" formulation in the literature [8]. The constraints in Eq. (7c) – (7e) ensure that the relaxation is the convex combination of the vertices of the active partition in P. The traditional SOS-2 is an ordered set of non-negative variables, of which at most two can be non-zero and they must be consecutive in their ordering. In this case, instead of an ordered set of non-negative variables, we have an ordered set that is a sum of non-negative variables i.e., for any i P I, λpµi pδ ki qqq, with k “ t1, . . . , ni ´ 1u, forms an ordered set of non-negative variable sums where SOS-2 constraints are imposed. It is not difficult to establish that the PPR in Eq. (6) has the property of being locally sharp, i.e., the projection of the LP relaxation of Eq. (6) to px, wq-space is exactly the convpXq. This is a property of the SOS-2 formulation for piecewise, linear approximations of functions [26]. In the literature, there are other vertex-based formulations which use this variable space and other formulations which use a logarithmic number of binary variables [26, 10, 27], however, we focus on providing a detailed comparison study between the SOS-2-based PPR in (7) and its recursive counterparts (in Sec. 3).
x“
16 ÿ
Jou
PPR Example. The PPR for a bilinear term with |P| “ 9 (3 partitions on each variable), shown in Fig. 1, is given in Eq. (8).
s“1
ˆs, w “ λs x
16 ÿ
s“1
ˆ s q, λs φpx
16 ÿ
s“1
(8a)
λs “ 1,
y1,1 ` y1,2 ` y1,3 “ 1, y2,1 ` y2,2 ` y2,3 “ 1,
(8b)
λ2 ` λ6 ` λ10 ` λ14 ď y1,1 ` y1,2 ,
λ5 ` λ6 ` λ7 ` λ8 ď y2,1 ` y2,2 ,
(8d)
λ4 ` λ8 ` λ12 ` λ16 ď y1,3 ,
λ13 ` λ14 ` λ15 ` λ16 ď y2,3 ,
(8f)
λ1 ` λ5 ` λ9 ` λ13 ď y1,1 ,
λ3 ` λ7 ` λ11 ` λ15 ď y1,2 ` y1,3 , λs ě 0 @s “ 1, . . . , 16
and,
λ1 ` λ2 ` λ3 ` λ4 ď y2,1 ,
λ9 ` λ10 ` λ11 ` λ12 ď y2,2 ` y2,3 ,
yi,k P t0, 1u
@i “ 1, 2, k “ 1, 2, 3. 4
(8c) (8e) (8g)
Journal Pre-proof
3. Recursive Piecewise Polyhedral Relaxation for a Multilinear Term with |I| ě 3
C
ź iPI
xi
GR-PPR
lP repro of
For a general multilinear term with |I| ě 3 and no partitions i.e., |P| “ 1, the recursive McCormick relaxation successively uses the McCormick relaxation in Eq. (4) on bilinear groupings of the terms. We extend this approach to build a R-PPR by recursively grouping bilinear terms and applying the PPR, given in Eq. (6), to each bilinear term. For a particular grouping of bilinear terms with |I| ě 3, the R-PPR is succinctly expressed as CB GPPR FPPR A EPPR PPR xx1 ¨ x2 y “ ... ¨ x|I|´1 ¨ x|I|
(9)
Notice that a different R-PPR can be created by choosing a different order of recursive grouping of variables, leading to different relaxation qualities [28, 29]. In computational studies, we present the performance of different groupings of bilinear terms in addition to the grouping presented in (9). Also, while applying R-PPR, only the original space of variables in the multilinear term are partitioned, and not the lifted/auxiliary variables resulting from recursive relaxations. 4. MILP-based Piecewise Formulation for Recovering Feasible Solutions
Jou
rna
Thus far, the focus of the paper has been on developing piecewise polyhedral relaxations for multilinear terms. The optimal solution to this relaxation provides an active partition for each variable. We let P ˚ and S ˚ denote the optimal active partition and the extreme points of the active partition chosen by the the PPR of the multilinear term, respectively. Let r`˚ , u˚ s denote the variable bound vector defined by the active partition P ˚ . Applying the convex hull formulation given in Eq. (5) for x P r`˚ , u˚ s yields the same solution produced by the PPR. Additionally, enforcing the solution to lie on the one-dimensional faces of the polytope in Eq. (5) for variable bounds r`˚ , u˚ s yields a feasible solution to the multilinear term, which is straight-forward to observe. This provides a simple MILP-based technique to compute a feasible solution of the original nonlinear problem of interest. Hence, in this section, we present an alternate MILP that recovers a feasible solution by requiring the solution to lie on an edge of the polytope which represents the relaxation of the multilinear term’s active partition. The formulation presented below is easily generalized to the case of requiring a solution to lie on the edges of the convex hull relaxation for any partition, but we omit this for the sake of simplicity of exposition. To enforce these additional constraints, we introduce some additional notation. Let V and E denote the vertices and the edges (one-dimensional faces) of the polytope defined by Eq. (5) for x P r`˚ , u˚ s. We know that |V| “ 2|I| and |E| “ |I| ¨ 2|I|´1 . Then, for every v P V, we define a set function γpvq “ te P E : e is incident on vu. Additionally, we introduce a binary variable ze for e P E that takes a value 1 if the solution lies on the edge e. Given these notations, the piecewise formulation that enforces the solution to lie on one of the edges of the polytope defined by Eq. (5) for x P r`˚ , u˚ s is defined by ! ) px, w, λ, zq P r`˚ , u˚ s ˆ R ˆ ∆|V| ˆ t0, 1u|E| | px, w, λ, zq satisfies Eq. (11) (10) where, Eq. (11) is given by ÿ ÿ ˆs, w “ ˆ s q, x“ λs x λs φpx ˆ s PS ˚ x
(11a)
ˆ s PS ˚ x
5
Journal Pre-proof
ePE
ze “ 1,
λs ě 0,
λv ď
ÿ
ePγpvq
(11b)
ze @v P V,
ˆ s P S ˚ , and, @x
ze P t0, 1u @e P E,
(11c)
lP repro of
ÿ
(11d)
where Eqs. (11a) and (11b) ensure that any solution to the formulation (11) is a convex combination of the vertices that constitute a one-dimensional face of the polytope in Eq. (5) and that only one face is active. Note that the above described method is also applicable if the original problem is an MINLP, as recovering feasible solutions to MINLPs is by itself an arduous task using off-the-shelf solvers. In the next section, we present extensive computational results that compare PPR with R-PPR, and the quality of recovered feasible solutions. 5. Computational Results
In this section, we present three sets of computational results. The first two sets of results use the NLP instance shown in Eq. (12) and the third set of results uses 60 nonlinear optimization problems that contain only multilinear terms taken from [30]. F “ max px1 x2 x3 x4 ` x3 x4 x5 x6 ` x5 x6 x7 x8 q
subject to: 100x1 ´ x2 ´ x3 ` 833x4 ` 95x5 ` x6 ´ x7 ` 100x8 ď 50000,
100 ď x1 ď 500, 1000 ď x2 , x3 ď 2000, 10 ď xi ď 100 @i P t4, . . . , 8u.
(12a) (12b) (12c)
Jou
rna
The first set of results shows the strength of the various relaxations (PPR and R-PPRs) and the effectiveness of feasible solution recovery performed on the active partition obtained from the respective relaxations. The second set of results is aimed at comparing the different ways in which feasible solution recovery can be performed by using the template formulation in Sec. 4 in terms of computation time and quality of the feasible solutions. The final set of results compares the strength of the proposed PPRs on instances used in the literature for which the optimal solutions are known a-priori. For both the first and the third set of computational experiments, we partition the variables of the multilinear terms with uniformly located discretization points. All models were implemented in Julia using JuMP v0.19.5 [31] and the experimental results used Gurobi v8.0 [32] with Gurobi’s presolver and heuristics turned off. All nonconvex models were solved using BARON v19.12.7 with CPLEX 12.10 and Ipopt 3.12.8 as BARON’s sub-solvers. The computational experiments were performed on an Intel Xeon CPU L5420 @2.50GHz with 64 GB RAM. 5.1. Effectiveness of the formulations on the NLP in (12) Table 1 compares the objective value of the PPR with three R-PPRs (denoted by upper bound, U B), each with a different variable grouping, on the NLP in Eq. (12). Note that the choice of grouping changes the quality of the relaxation. The active partition of each variable in the relaxed solution is used to compute a feasible solution (whose objective is denoted by lower bound, LB) through formulation Eq. T T ´LB (11). In this table, U B and LB gaps are given by U B´OP ¨ 100 and OP LB ¨ 100, respectively; OP T UB corresponds to the global optimum value of (12) which is equal to 3.2642E10. It is clear from these results that, in both the relaxation and the feasible solution recovery formulations, PPR outperforms all the R-PPRs. For example, with only 4 partitions, PPR finds a relaxed solution (U B) that is better 6
Journal Pre-proof
# of partitions per variable
2
% gaps xa xb xc xd xxa xxb xc yyxd xxxa xb yxc yxd xa xxb xxc xd yy
LB, U B 2.33, 23.99 2.33, 65.47 2.33, 65.47 2.33, 47.37
lP repro of
than any R-PPR (up to 12 partitions). It is also clear that the quality of the feasible solutions (LB) recovered based on the active partition of the PPR formulation is consistently better than those of R-PPR’s for every chosen number of partitions. Note that, the LB gaps do not necessarily decrease monotonically as the number of partitions increase; the new partitions for PPR/R-PPR formulations yield different active partitions, thus leading to different local optimal solutions. In these experiments, the run-times of both PPR and R-PPR formulations were less than a minute and hence we do not report them explicitly. As is the case with any piecewise formulation, the computational performance of PPR degrades as the number of partitions increases on larger NLPs, however the strength of PPR’s relaxation remains superior compared to any other recursive variation. 4
6
8
10
12
LB, U B 0.15, 3.20 36.74, 25.37 36.74, 25.37 5.73, 25.27
LB, U B 1.11, 2.98 2.33, 21.73 2.33, 21.73 1.11, 16.78
LB, U B 0.15, 0.83 0.15, 14.72 0.15, 14.72 0.15, 12.29
LB, U B 0.00, 0.69 0.00, 12.98 0.00, 12.98 0.00, 9.59
LB, U B 0.05, 0.43 2.33, 10.34 2.33, 10.34 1.11, 8.19
Table 1: Relative optimality gaps of upper and lower bounds in percent for the PPR and the R-PPRs for varying number of partitions per variable. The x¨y represents the grouping of the different variables into bilinear terms in the R-PPRs.
rna
BARON’s performance on the NLP in (12): In order to compare the relative optimality gaps obtained by the PPR formulation in Table 1, we also solved instance (12) using BARON solver. The best optimality gap BARON reports after a run time of one hour is 20.36% with 153,555 open spatial branch-and-bound nodes, where the best upper bound reported by BARON is 4.0986 ˆ 1010 and the best lower bound by Ipopt within BARON is 3.2642 ˆ 1010 . In contrast, the PPR reports 3.34% with only four partitions per variable (see table 1). The best gap PPR reports is 0.47% (12 partitions) which is close to global optimum, thus reinforcing the value of the convex hull formulations proposed in this letter. 5.2. Feasible solution recovery for the NLP in (12)
Jou
Here, we discuss the quality of the feasible solutions obtained by using different variants of the MILP formulation presented in Sec. 4 for recovering solutions to the NLP described in (12). To this end, we introduce auxiliary variables t1 , t2 and t3 , such that the objective function of (12) will be to maximize pt1 ` t2 ` t3 q subject to constraints t1 “ x1 x2 x3 x4 , t2 “ x3 x4 x5 x6 , t3 “ x5 x6 x7 x8 , (12b) and (12c). The MILP formulation, as presented in Sec. 4, when applied to each multilinear term in Eq. (12), results in a feasible solution to the NLP. This formulation computes the best solution that lies on the one-dimensional face of the active partition of each partitioned variable as dictated by the solution of the relaxation. From here on, we refer to this formulation as F a . Instead, the formulation can be extended to compute the best solution that lies on the one-dimensional face of any combination of the variable partition (the formulation is a trivial extension of the formulation in Eq. (11), and is not presented). We denote this formulation as F1f . Though formulation F1f may sometimes lead to better quality feasible solutions, the number of one-dimensional faces grows exponentially with the number of partitions on each variable and thus can result in a MILP with a large number of binary variables. Note 7
Journal Pre-proof
# of partitions
xa xb xc xd
xxa xxb xc yyxd
xxxa xb yxc yxd
xa xxb xxc xd yy
Fa
F1f
Fa
F1f
F2f
Fa
F1f
F2f
Fa
F1f
F2f
2 4 6 8 10 12
2.33 0.15 1.11 0.15 0.0 0.05
2.33 0.15 1.11 0.15 0.0 0.05
86.25 0.15 4.22 0.15 1.12 0.15
2.33 36.74 2.33 0.15 0.0 2.33
0.17 0.05 0.03 5E-3 5E-4 6E-5
86.25 0.15 4.22 0.15 1.12 0.15
2.33 36.74 2.33 0.15 0.0 2.33
0.17 0.01 0.03 0.05 5E-3 6E-5
2.33 2.33 1.11 0.17 0.0 0.05
2.33 5.73 1.11 0.15 0.0 1.11
0.01 2E-4 4E-3 2E-5 4E-5 2E-5
lP repro of
per variable
Table 2: Quality (percent gaps) of feasible solutions obtained using different variants of formulation (11) on the NLP in Eq. (12).
rna
that while applying F1f on various recursive grouping of bilinear terms, we do not partition the auxiliary lifted variables. The third formulation, F2f , is applicable only when using recursive groupings of bilinear terms. Here, we enforce partitions on the auxiliary lifted variables introduced while performing the recursion. We also remark that comparing the solutions of F a and F1f with F2f is not a fair comparison. F2f naturally provides better quality solutions since more variables are partitioned. Nevertheless, we compare all the formulations in Table 2. In this table, all percent gaps are computed relative to the T ´LB global optimum value (OP T ), that is, OP LB ¨ 100. As expected and consistent with our previous observations, even in the feasible solution recovery process, the convex hull representation of the multilinear terms outperforms the recursive approaches when comparing formulations F a and F1f . As for formulation F2f , it outperforms even the convex hull formulation as it partitions the auxiliary lifted variables. As a side note, it is not difficult to prove that for a fixed partition count, F2f on any recursive bilinear grouping will always produce the best feasible solution among all the formulations i.e., F1f and F a . 5.3. Strength of the relaxations on instances from [30]
Jou
Our final set of results are based on a collection of 60 multilinear problems taken from [30], which are available at: https://minlp.com/downloads/testlibs/barlibs/mult3.zip. The bounds on all the variables for every problem instance in [30] are r0, 1s. To demonstrate the strength of the relaxations discussed in this paper, we created a modified set of instances. For each variable, we uniformly distributed the lower and upper bounds in the interval r0.1, 0.2s and r0.9, 1s, respectively. This modification is motivated by the following known result in [13]–the recursive McCormick relaxation for a multilinear term is the convex hull when all variable bounds are symmetric about the origin or are r0, 1s. Though the bounds were modified, the global optimum values for all the instances remained the same as those in [30] when solved using BARON. For each problem instance we relaxed the problems using PPR and a R-PPR with a lexicographic ordering. A time limit of one hour was imposed on every run of the instance. The results of all the runs are summarized using the box-plot shown in Fig. 2. In Fig. 2, the optimality gap for the PPR in most instances is 0% when the number of partitions per variable is three. This is not the case for R-PPR. However, as expected, Fig. 2(a) suggests that when the bound on every variable in the problem is r0, 1s, R-PPR is generally the better alternative. This is not surprising since similar recursive relaxations are known to capture the convex hull of a multilinear term with symmetric variable bounds [13]. Once the variable bounds are randomized 8
40
PPR
lP repro of
Journal Pre-proof
3500
R-PPR
Computation time (sec)
Optimality gap (%)
30
3000
20
PPR
R-PPR
2500
2000
1500
1000
10
500
0
0
2
3 4 Number of partitions
2
3 4 Number of partitions
(a) Percent gaps and run times for instances with variable bounds set to r0, 1s. The optimality gaps for the PPRs increase from 3 partitions to 4 partitions because of the computation time limit of 1 hour imposed on every run.
35
PPR
R-PPR
15 10 5
2
PPR
R-PPR
2500
2000
1500
1000
500
0
Jou
0
Computation time (sec)
25
3000
rna
Optimality gap (%)
30
20
3500
3 4 Number of partitions
2
3 4 Number of partitions
(b) Percent gaps and run times for instances with variable bounds randomly chosen as described earlier. In general, the R-PRR has better computation times because it has fewer multiplier variables (λ). Nevertheless, the PPR closes the optimality gap with fewer partitions despite taking a little more computation time. Figure 2: The optimality gap in % refers to the relative gap between the objective value of the relaxation of the problem and the global optimum to the problem instance. The computation time is the time taken to solve the PPR and the R-PPR to optimality.
9
Journal Pre-proof
6. Concluding Remarks
lP repro of
and made asymmetric, the PPR formulation outperforms the recursive relaxation (R-PPR) in both optimality gap and computation time (Fig. 2(b)). Though the (bottom right) figure suggests that the computation times of PPR are larger than that of R-PPR, it is worth noting that the number of partitions and the run times necessary for the PPR to achieve 0% optimality gap on all 60 instances is much superior when compared to the R-PPR.
This letter presents piecewise polyhedral relaxation formulations for a multilinear term with spatial disjunctions on every variable. An SOS-2-based formulation that generalizes similar formulations for piecewise linear approximations of functions was developed to formulate the convex hull of a single multilinear term. Extensive computational experiments demonstrated that our proposed formulations have significant computational advantages over the well-known recursive grouping of bilinear terms. A MILP-based piecewise formulation was developed to recover high-quality feasible solutions for problems with multilinear terms. Given exponential growth in the number of variables in the PPR formulation, future work will focus on developing formulations with subsets of variables for those rare cases when the number of variables in a multilinear term is large and further add valid inequalities, by generalizing the ideas proposed in [30]. Acknowledgements
7. References
rna
The work was funded by the Center for Nonlinear Studies (CNLS) at LANL and LANL’s Directed Research and Development (LDRD) projects, “20170201ER: POD: A Polyhedral Outer-approximation, Dynamic-discretization optimization solver" and “20190590ECR: Discrete Optimization Algorithms for Provably Optimal Quantum Circuit Design". This work was carried out under the U.S. DOE Contract No. DE-AC52-06NA25396.
[1] M. Tawarmalani, N. V. Sahinidis, A polyhedral branch-and-cut approach to global optimization, Mathematical Programming 103 (2) (2005) 225–249. [2] M. Hasan, I. Karimi, Piecewise linear relaxation of bilinear programs using bivariate partitioning, AIChE journal 56 (7) (2010) 1880–1893.
Jou
[3] L. S. de Assis, E. Camponogara, B. Zimberg, E. Ferreira, I. E. Grossmann, A piecewise McCormick relaxation-based strategy for scheduling operations in a crude oil terminal, Computers & Chemical Engineering 106 (2017) 309–321. [4] P. M. Castro, Tightening piecewise McCormick relaxations for bilinear problems, Computers & Chemical Engineering 72 (2015) 300–311. [5] H. Nagarajan, M. Lu, E. Yamangil, R. Bent, Tightening McCormick relaxations for nonlinear programs via dynamic multivariate partitioning, in: International Conference on Principles and Practice of Constraint Programming, Springer, 2016, pp. 369–387.
10
Journal Pre-proof
lP repro of
[6] H. Nagarajan, M. Lu, S. Wang, R. Bent, K. Sundar, An adaptive, multivariate partitioning algorithm for global optimization of nonconvex programs, Journal of Global Optimization 74 (4) (2019) 639–675. [7] S. Sridhar, J. Linderoth, J. Luedtke, Locally ideal formulations for piecewise linear functions with indicator variables, Operations Research Letters 41 (6) (2013) 627–632. [8] J. P. Vielma, S. Ahmed, G. Nemhauser, Mixed-integer models for nonseparable piecewise-linear optimization: Unifying framework and extensions, Operations research 58 (2) (2010) 303–315. [9] J. Lee, D. Wilson, Polyhedral methods for piecewise-linear functions i: the lambda method, Discrete applied mathematics 108 (3) (2001) 269–285. [10] M. Padberg, Approximating separable nonlinear functions via mixed zero-one programs, Operations Research Letters 27 (1) (2000) 1 – 5. [11] G. P. McCormick, Computability of global solutions to factorable nonconvex programs: Part I–convex underestimating problems, Mathematical programming 10 (1) (1976) 147–175. [12] F. A. Al-Khayyal, J. E. Falk, Jointly constrained biconvex programming, Mathematics of Operations Research 8 (2) (1983) 273–286. [13] J. Luedtke, M. Namazifar, J. Linderoth, Some results on the strength of relaxations of multilinear functions, Mathematical programming 136 (2) (2012) 325–351. [14] N. V. Sahinidis, Baron: A general purpose global optimization software package, Journal of global optimization 8 (2) (1996) 201–205. [15] S. Vigerske, A. Gleixner, SCIP: Global optimization of mixed-integer nonlinear programs in a branch-and-cut framework, Optimization Methods and Software 33 (3) (2018) 563–593.
rna
[16] P. Belotti, Couenne: a user’s manual, Tech. rep., Lehigh University (2009).
[17] C. A. Floudas, Deterministic global optimization: theory, methods and applications, Vol. 37, Springer Science & Business Media, 2013. [18] A. D. Rikun, A convex envelope formula for multilinear functions, Journal of Global Optimization 10 (4) (1997) 425–437.
Jou
[19] E. Balas, Disjunctive programming, Annals of discrete mathematics 5 (1979) 3–51. [20] R. G. Jeroslow, J. K. Lowe, Modelling with integer variables, Mathematical Programming at Oberwolfach II (1984) 167–184. [21] R. Misener, C. A. Floudas, Global optimization of mixed-integer quadratically-constrained quadratic programs (MIQCQP) through piecewise-linear and edge-concave relaxations, Mathematical Programming 136 (1) (2012) 155–182. [22] R. Misener, J. P. Thompson, C. A. Floudas, APOGEE: global optimization of standard, generalized, and extended pooling problems via linear and logarithmic partitioning schemes, Computers & Chemical Engineering 35 (5) (2011) 876–892. 11
Journal Pre-proof
[23] J. Huchette, J. P. Vielma, A combinatorial approach for small and strong formulations of disjunctive constraints, Mathematics of Operations Research 44 (3) (2019) 793–820.
lP repro of
[24] E. Beale, J. J. Forrest, Global optimization using special ordered sets, Mathematical Programming 10 (1) (1976) 52–69. [25] E. M. L. Beale, J. A. Tomlin, Special facilities in a general mathematical programming system for non-convex problems using ordered sets of variables, OR 69 (447-454) (1970) 99. [26] J. Huchette, J. P. Vielma, Nonconvex piecewise linear functions: Advanced formulations and simple modeling tools, arXiv preprint:1708.00050. [27] R. Misener, C. A. Floudas, Global optimization of mixed-integer quadratically-constrained quadratic programs (miqcqp) through piecewise-linear and edge-concave relaxations, Mathematical Programming 136 (1) (2012) 155–182. [28] P. Belotti, S. Cafieri, J. Lee, L. Liberti, A. J. Miller, On the composition of convex envelopes for quadrilinear terms, in: Optimization, Simulation, and Control, Springer, 2013, pp. 1–16. [29] J. Lee, D. Skipper, E. Speakman, Algorithmic and modeling insights via volumetric comparison of polyhedral relaxations, Mathematical Programming 170 (1) (2018) 121–140. [30] X. Bao, A. Khajavirad, N. V. Sahinidis, M. Tawarmalani, Global optimization of nonconvex problems with multilinear intermediates, Mathematical Programming Computation 7 (1) (2015) 1–37. [31] I. Dunning, J. Huchette, M. Lubin, JuMP: A modeling language for mathematical optimization, SIAM Review 59 (2) (2017) 295–320.
Jou
rna
[32] Gurobi Optimizer 7.5, http://www.gurobi.com.
12