Global optimization of bilinear programs with a multiparametric disaggregation technique

In this paper, we present the derivation of the multiparametric disaggregation technique (MDT) by Teles et al. (J. Glob. Optim., 2011) for solving nonconvex bilinear programs. Both upper and lower bounding formulations corresponding to mixed-integer linear programs are derived using disjunctive programming and exact linearizations, and incorporated into two global optimization algorithms that are used to solve bilinear programming problems. The relaxation derived using the MDT is shown to scale much more favorably than the relaxation that relies on piecewise McCormick envelopes, yielding smaller mixed-integer problems and faster solution times for similar optimality gaps. The proposed relaxation also compares well with general global optimization solvers on large problems.


Introduction
Bilinear programs, for the purpose of this paper, can be written as the following nonconvex nonlinear programming problem: Min z = f 0 = (i, j)∈B L 0 Subject to f q = (i, j)∈B L q a i jq x i x j + h q (x) ≤ 0 q ∈ Q\{0} (P) x ∈ S ∩ ⊂ R n where h q (x) is convex and twice differentiable, a i jq is a scalar with i ∈ I, j ∈ J , and q ∈ Q represents the set of all functions f q , including the objective function f 0 and all constraints.B L q is an (i, j)-index set which defines the bilinear terms present in the problem.Although i = j for strictly bilinear problems, i = j can be allowed to accommodate quadratic problems.The set S ⊂ R n is convex, and ⊂ R n is an n-dimensional hyperrectangle defined in terms of the initial variable bounds x L and x U : The global optimization of bilinear programs of the form of (P) is important in such areas as water networks and petroleum blending operations [1][2][3][4].Nonconvex, bilinear constraints are required to model the mixing of various streams in these systems, and are in some cases the only nonlinearities in the models.The pooling problem, stemming from the original Haverly paper [5], contains these bilinear constraints and has received much attention in the literature [1,2,[6][7][8][9].Recently, Misener and Floudas have demonstrated a novel logarithmic relaxation for modeling bilinear terms with piecewise McCormick envelopes while addressing various classes of pooling problems [1].
Water network optimization problems containing bilinear terms have also received much attention in the literature [3,4,[10][11][12].The same blending constraints present in the pooling problem are present in water network problems, and thus numerous advances in solving bilinear programs have been made addressing these problems.
The global optimization of general nonconvex bilinear programs has received significant attention in the literature [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27].The convex McCormick envelopes [1] coupled with spatial branch and bound search frameworks have been the basis for many of these global optimization techniques, with piecewise McCormick envelopes being a more recent development.Variations of this approach have been suggested, generalizing the convex envelopes to piecewise over-and under-estimators [1,17].Novel ways of representing bilinear terms through reformulation have been another approach reported in the literature [16].Misener et al. [1], building on the work of Vielma and Nemhauser [28,29], have shown that a relaxation of the bilinear terms can be achieved with a logarithmic number of binary variables.Teles et al. [30] have introduced a technique to approximate polynomial constraints that exhibits a similar property, where the main element is the discretization of a subset of the variables.
In this paper, we show that the mixed-integer constraints of the multiparametric disaggregation technique (MDT) presented in [30], when applied to the bilinear terms in program (P), can be derived from disjunctive programming and convex hull reformulation.This approximation technique can be used under some conditions to obtain an upper bound.As the main novelty, we propose a rigorous lower bounding formulation derived from the upper bounding constraints.After proving that the solution from the upper and lower bounding formulations converge to that of the original nonlinear formulation in the limit of an infinite number of discretization intervals, two global optimization algorithms are proposed based on such bounds.Finally, we conclude with a comparison of this new relaxation approach and the one based on piecewise McCormick envelopes, and report computational results on small and large problems.

Discretization with multiparametric disaggregation
Given a nonconvex bilinear term w i j = x i • x j , the MDT described by Teles et al. [30] can be used to obtain an upper bound on problem (P).The complete formulation can be derived from a generalized disjunctive programming (GDP) model [31] followed by a convex hull reformulation [33] and exact linearization [32].For simplicity in the notation, we first rewrite the bilinear product w i j = x i • x j as a single bilinear term w = u • v.This product can be represented exactly with the following constraints and disjunction: where v is discretized through the disjunction in (3) that selects one digit k ∈ K = {0, 1, . .., 9} for each power ∈ Z.Here we assume a basis of 10, although other bases can be selected [11].Note that since (3) is defined over the domain of all the integer numbers, this implies an infinite number of disjunctions.Furthermore, v can represent any positive real number.First, we consider the convex hull reformulation [33] of the disjunction in (3) after which we introduce the disaggregated variables, Substituting ( 5) into (4) and then into (2) leads to the fully discretized (but still exact representation) of v: Considering the product w = u • v by substituting ( 8) into (1) leads to (9), which involves nonlinear terms u • z k, .
Performing an exact linearization [32], we introduce new continuous variables ûk, = u • z k, so that: and ûk, is non-negative, we introduce the following lower and upper bounding constraints, where u U and u L are the non-negative upper and lower bounds on u.
Finally, multiplying equation ( 6) by u and replacing the bilinear terms by the new continuous variables, results in (12).The full set of mixed integer linear constraints for the exact representation of bilinear product w = u • v is thus given by Eqs.(6)(7)(8) and (10)(11)(12).
We should note that the same model formulation can be obtained by the special structured RLT reformulation as described in Sherali et al. [48].Constraints (6)(7)(8) are directly linked to the choice of using a base-10 representation for variable v.Then, the bounds of variable u are known: u L ≤ u ≤ u U .Multiplying the lower and upper bounding constraints u − u L ≥ 0 and u U − u ≥ 0 by z k, ≥ 0 and replacing the bilinear term u • z k, by ûk, , leads to (11).Performing the same variable replacement in the constraint defining variable w, results in (10) while ( 12) can be obtained as previously described.

Upper bounding formulation
Since in practice it is infeasible to compute the infinite sums over all integers, we represent v to a finite level of precision, v , leading to a continuous but approximate representation of the bilinear term, w .The constraints in (8) and (10) are modified in (13)(14) to allow for a maximum power of 10 (P) and a minimum power of 10 ( p).For the remaining constraints, (6-7) and (11)(12), it suffices to replace ∈ Z with ∈ L = {p, p + 1, . . ., P}.These sets of constraints correspond to the equations proposed by Teles et.al. [30]: When we incorporate them into problem (P) by redefining w i j = x i • x j , and selecting x j as the variable on which discretization is performed, the resulting problem (P ) will represent a mixed-integer approximation to the original problem: Subject to f q (x) = (i, j)∈B L q a i jq w i j + h q (x) ≤ 0 q ∈ Q\{0} (P ) where x j and w i j represent discrete and continuous approximations to the variables, respectively.Note that if the convex functions h q (x) are linear, problem (P ) represents a mixed integer linear program (MILP).Otherwise it corresponds to an MINLP.
Further, note that problem (P ) is a restricted version of problem (P), or equivalently problem (P) is a relaxation of problem (P ).It then follows that the solution of problem (P ) either yields an upper bound on problem (P) such that z = z U ≥ z, or else problem (P ) is infeasible.This restricted feasible region can be seen in Figs. 1 and 2.

Infeasibilities in the discretized problem
The mixed-integer approximation problem (P ) can be infeasible even if the original problem (P) is feasible.For example, if bounds such as 10 p−1 ≤ x j ≤ 2 × 10 p−1 are enforced, (P ) will be infeasible, while such a constraint is completely valid and will not necessarily result in an infeasible problem (P).Thus, the parameters p and P must be chosen appropriately to avoid such infeasibilities.
Some general guidelines for ensuring precision-based feasibility can be established.For example, the largest power of 10 (P) must be large enough such that 10 P is of the same order of magnitude as the upper bound on x j .More precisely: P = log 10 x U j .Additionally, p must be small enough to ensure that at least one (and preferably many) discretization points lie between the lower and upper bounds for x j .Thus, p ≤ P is the absolute minimum requirement, but feasibility is more likely as p is decreased.Note that these guidelines do not guarantee feasibility of (P ) in all cases, but represent the minimum level of precision needed given reasonable bounds on x j .

Lower bounding with multiparametric disaggregation
To obtain a lower bounding problem using multiparametric disaggregation, we first note that in the discretized problem, there always exists a gap between discretization points for a finite p.Thus, we can introduce a slack variable x j such that x R j = x j + x j , where x j is the discretized representation of x j from Sect. 2, x R j is the continuous representation of x j and the slack variable x j is bounded between 0 and 10 p .
Again switching to the notation w = u • v for the bilinear term, we have for the continuous representation of v, denoted as v R : For the continuous representation of the bilinear term, w R , note that: where w and v are given by (13)(14).The slack variable w replaces the bilinear term u • v that can be relaxed using the McCormick envelope, (18)(19), which in this case coincides with the RLT bound factor products u L ≤ u ≤ u U and 0 ≤ v ≤ 10 p : Introducing these constraints into Problem (P), and expressing the variables in terms of the original variables w i j = x i • x j , we obtain the new optimization problem, (PR): While (PR) does not exactly represent the product w i j = x i • x j and is feasible for values of w i j , x i , and x j that do not satisfy w i j = x i • x j , the bilinear term is feasible in (PR).Thus, (PR) is a relaxation of (P).The relaxed feasible region resulting from (PR) can be seen in Figs. 3 and 4.
The following property can be readily established from the above discussion: Property 1 The solution of problem (PR) yields a lower bound for problem (P), i.e. z R ≤ z.

Discussion of global optimization algorithms
The upper and lower bounding schemes described can be combined into a global optimization algorithm.First, the following property can be established: Thus, since the variables x j and w i j are eliminated, this yields problem (P ), and hence z approaches z R .
From Property 2, we can establish that as precision is increased (i.e.p approaches−∞), both (P ) and (PR) converge to the same value.Assuming P is large enough such that 10 P ≥ x U j , we can further state that (P ) and (PR) converge such that z = z R = z.

Algorithm 1
The first global optimization algorithm that can be established from the aforementioned upper and lower bounding is as follows.First, we start at some coarse level of discretization such that P ≥ p, and solve both (PR) and (P ).If the difference in solutions to the upper and lower bounding problems is sufficiently small, then the algorithm terminates; otherwise, precision is increased and the problems are resolved.The algorithm is then as follows: Algorithm 1 Step 0. Choose p = P = log 10 x U j Step 1. Solve (PR) to obtain the lower bound z R .
Step 3. If (z −z R )/z R ≤ ε, STOP, the solution is globally optimal.Otherwise, set p = p−1, and return to step 1.

Algorithm 2
While Algorithm 1 follows most naturally from the problems (P ) and (PR), it has several shortcomings.Notably, because (P ) and (PR) are fairly similar and increase similarly in problem size as precision is added, solving (P ) and (PR) repeatedly becomes increasingly expensive.Thus, an alternative method for obtaining an upper bound is to use a local NLP algorithm in place of solving the problem (P ).

Algorithm 2
Step 0. Choose p = P = log 10 x U j Step 1. Solve (PR) to obtain the lower bound z R .
Step 2. Solve (P) using a local NLP algorithm to obtain some upper bound z using the solution to (PR) as a starting point.Step 3. If (z − z R )/z R ≤ ε, STOP, the solution is globally optimal.Otherwise, set p = p −1 and return to step 1.
Algorithm 2 is generally more computationally efficient than Algorithm 1, as the solution of (P) using a local NLP algorithm is obtained much faster than the increasingly large MILP that (P ) becomes as P − p grows.

Extensions to MINLP
The algorithms in Sects.4.1 and 4.2 can be readily extended for solving MINLPs with the following general form: where h q (x, y) is jointly convex in x and y.For most cases in practice the variables y appear in linear form in these terms [34].Analogous problems (P -MINLP) and (PR-MINLP) can be derived, and Algorithm 1 can be used without modification.However, Algorithm 2 requires some modification, as (P-MINLP) cannot be solved using a local NLP algorithm because it is an MINLP.A heuristic can be used to compute the upper bound as in Algorithm 2:

123
Step 2. Fix the binary variables y in (P-MINLP) to the values found by the solution of (PR-MINLP) in Step 1, reducing it to an NLP.Solve (P-MINLP) with these fixed binary variables using a local NLP algorithm to obtain some z using the solution to (PR-MINLP) as a starting point.Step 3. If (z − z R )/z R ≤ ε, STOP, the solution is globally optimal.Otherwise, set p = p −1 and return to step 1.
A disadvantage of this algorithm is that it could take, in the worst case, an infinite number of iterations to converge.By fixing the binary variables to that of the solution to (PR-MINLP), (P-MINLP) can be rendered infeasible, and algorithm would continue until (PR-MINLP) has enough discretization points to exactly represent (P-MINLP).If this occurs, heuristics or other approaches may be utilized to obtain an upper bound in place of solving (P-MINLP) with fixed binary variables.However, in practice, this is unlikely to occur.

Comparison with piecewise McCormick envelopes
A common approach to solving bilinear programs of the form (P) is to bound the bilinear terms using McCormick convex envelopes [35].This formulation results in a lower bounding LP if the only nonlinearities are bilinear, or a lower bounding convex NLP if there are remaining convex nonlinearities.For each bilinearity w i j = x i • x j , we introduce instead the following constraints: This relaxation yields a lower bound on problem (P).However, this lower bound can be weak depending on the bounds on x i and x j .To improve the quality of the lower bound, these convex envelopes can be used on discretized portions of the variable range.Thus, piecewise McCormick envelopes can be introduced in (P) to obtain a tighter lower bound at the cost of becoming an MILP or convex MINLP [10,17,36].The envelopes, defined over a set of N points on the domain of variable x j , can be represented by the following disjunctive constraints: Applying the convex hull reformulation [33] for the above disjunctive constraints yields By adding these piecewise McCormick envelope constraints (22) to problem (P), we can define a new relaxed MILP or convex MINLP, (PR-PCM).Furthermore, we can compare the performance of (PR-PCM) to the lower bounding problem derived from multiparametric disaggregation, (PR).

Lower bounding problems
In order to compare the lower bounds predicted by multiparametric disaggregation and piecewise McCormick envelopes, we solve the relaxation problems (PR) and (PR-PCM) resulting from multiparametric disaggregation and piecewise McCormick, respectively.For the specific example Problem (P1), we can derive analogous relaxed problems (P1R) and (P1R-PCM) using the MDT and piecewise McCormick envelopes, respectively.Note that x 1 is the variable being discretized.
Problem size and computational results are shown in Table 1 for the lower bounds predicted by (P1R) at P = 0 and p = {0, −1, . .., −6} and (P1R-PCM) at N = {1, 15, 150, 1,500, 15,000}.Solving (P1R-PCM) at various levels of discretization, it becomes clear that problem size increases approximately exponentially with each order of magnitude decrease in the optimality gap.However, this is not the case when solving (P1R) since as precision is added, i.e. p is decreased, a linear relationship holds.Note that the discretization level of MDT at p = −1 matches that of PCM at N = 15( p = −2 is equivalent to N = 150 and so on) and that the lower bounds are exactly the same.Upper bounds are also reported by using as a starting point the solution of the lower bounding relaxation problem as in Algorithm 2 in Sect.4.2.Thus, the MDT columns are equivalent to the iterations of the algorithm with the reported CPU time being the total time required for the lower and upper bounding problems.These times would need to be accumulated to compare Algorithm 2 directly to BARON, unless the discretization level ( p) is chosen a priori.Results for solving the original NLP (P1) with the global optimization solver BARON are reported in the second column.
As seen in Table 1, the relaxed problems (P1R) and (P1R-PCM) are considerably larger than the original NLP (P1) since the addition of both continuous and binary variables increases problem size.However, note that the multiparametric disaggregation problem (P1R) requires fewer additional continuous variables and fewer constraints than when utilizing piecewise McCormick envelopes in (P1R-PCM).For a single McCormick envelope, a relatively weak lower bound of −1.5 is obtained leading to an upper bound that is in fact a suboptimal solution.Multiparametric disaggregation with p = −4 approaches an optimality gap of 0.003 %, while piecewise McCormick envelopes, even with 1,500 partitions, only approach an optimality gap of 0.034 %.Further discretization and refinement of the solution is quickly solved using multiparametric disaggregation, as the p = −5 and p = −6 problems are solved in <1.2 s and fewer than 130 variables.In contrast, increasing the number of partitions to 15,000 in (P1R-PCM), leads to a MILP that cannot be solved in 1-h of computational time.

Numerical experiments
The performance of the underestimating problems from multiparametric disaggregation (PR) and piecewise McCormick (PR-PCM) is further evaluated through the solution of three    additional small test problems from the literature for different accuracy levels.Since in all problems the functions h q (x) in (P) are linear, the resulting bounding MILP problems were solved in GAMS 23.8.2 [38] using CPLEX 12.4 (1 thread) [39].Default options were used except for the relative optimality tolerance, equal to 10 −6 and a maximum computational effort equal to 3,600 CPU seconds.Similarly as in problem (P1), the original nonlinear programs were solved by CONOPT 3 [40], following initialization with the values from the MILP, and by BARON 10.2 [41].We also report results for larger problems in the areas of water networks and multiperiod blending.Such problems were also solved by GloMIQO 1.0.0 [49].The computational experiments were performed on an Intel i7 950 processor running at 3.07 GHz, with 8 GB of RAM, running Windows 7. 7

P3
Test problem (P3) corresponds to Problem 106 in Hock and Schittowski [43].While the objective function and the first three constraints are linear, there are three bilinear inequalities.To study how the choice of parameterized/partitioning variables affects computational performance, three versions are considered: (a) Test problem (P4) is taken from Shen and Zhang [44] and after a few transformations the following bilinear problem results.It features a global optimal solution f = 460212.25 at (43.165468, 45, 70, 1.228374, 27.749229, 0.814083).Variables x 2 , x 5 and x 6 are chosen as the parameterized variables.
22.85714 ≤ x 5 ≤ 33, 0.714286 ≤ x 6 ≤ 10 7.2 Computational statistics Tables 1, 2, 3 and 4 give the computational results for problems (P1-P4) as a function of the discretization level, and the columns correspond to the iterations of Algorithm 2, where the relaxation problem is either generated by multiparametric disaggregation (PR), MDT columns, or by the piecewise McCormick envelopes (PR-PCM), PCM columns.We provide the problem size, lower bound from the relaxation problem, and upper bound from a local NLP solver (the values of the latter remain the same independent of the accuracy level).The optimality gap and total computational effort (CPLEX plus CONOPT, the latter being almost negligible) are also reported.Note that the triplets in Table 2 are related to the number of discrete points for the three variables selected.The n-tuples in Tables 3 and 4 have a similar meaning.In problems (P2-P3), the discretization points of PCM match exactly those being used by MDT (e.g.second column of PCM should be compared with second column of MDT).An interesting observation is that for the same accuracy level, the relaxation from piecewise McCormick (PCM) is tighter than MDT in (P2) but identical in (P1) and (P3).If one looks at the original NLP problems, the latter problems have strictly bilinear terms whereas (P2) features a quadratic term in x 3 , which is one of the variables being discretized.Applying the piecewise relaxation to a quadratic term w ii = x i • x i makes it possible to use tighter bounds for both parts of x i in each partition n, see (23), whereas with MDT just one part is discretized, meaning the use of the overall bounds for the other part.This provides the    Italic values indicate that the maximum computational time was reached explanation for multiparametric disaggregation being less tight than piecewise McCormick, which can be observed in Fig. 5.
For (P4), the lower bounds of x 5 and x 6 are not integer values and so, while the number of discrete points is roughly the same, their location is not, making it possible for the relaxation from MDT to be better than the one from PCM.For instance, p = (−2, −2, −2) leads to a slightly better lower bound than N = (500, 1, 000, 1, 000).Due to the better performance of the new method, there are more columns for MDT than for PCM, meaning that higher accuracy levels, i.e. lower optimality gaps, can be achieved by the former for a given computational time.
While MDT is not as tight as PCM, the latter generates considerably larger MILP problems for the same accuracy level.More specifically, with PCM we get exactly an exponential increase in the number of binary variables with the change of power and roughly the same behavior with respect to the number of total variables and constraints.In other words, the number of 0-1 variables grows linearly with the number of partitions.In contrast, the number of 0-1 variables for MDT grows logarithmically with the number of discrete points and linearly with the change of power, keeping problems tractable for a wider accuracy range.Notice that with MDT all problems can be solved with an optimality gap of less than 0.01 % in less than one hour, while PCM closes only to a gap of 6.63 % for (P3).
Since problem size is related to the number of parameterized/partitioned variables, one might be tempted to keep this number as low as possible.However, the results for (P3) give opposite results, since the worst performance is obtained for 3 parameterized variables (gap = 0.009 %), followed by 4 (gap = 0.0067 % in 1898 CPUs) and then 5 (gap = 0.0071 % in 520 CPUs, which improves to just 0.0020 % for p = (−3, . .., −3), calculated using the best possible solution at time of termination).

Evaluation of Algorithm 1
Compared to Algorithm 2 in Tables 2, 3 and 4, Algorithm 1 generally leads to larger optimality gaps, particularly in the first iterations at coarse discretization levels, and it also requires more CPU time since two MILPs are solved at each iteration.Comparatively, solving (P) with a local solver can be done almost instantaneously.The lower bounding problem is the same as in Algorithm 2, while the upper bounding problem (P ) yields considerably worse bounds than the ones returned from the local NLP solver (recall that these were always global optimal solutions).As is illustrated in Fig. 6 for test problem (P3a), the solutions from both (P ) and (PR) typically become closer to the optimum with an increase in accuracy.However, it is clear that using a local NLP solver, as in Algorithm 2, is a superior approach.While it is possible for the solution of (P ) to be closer to the global optimum than the solution of (P) using a local NLP algorithm, in practice, this generally does not occur, especially when using the solution of (PR) as a starting point.

Results for larger problems
As seen in Tables 1, 2, 3 and 4, when compared to the commercial solver BARON, multiparametric disaggregation is competitive in (P1), (P2) and (P4) but is orders of magnitude slower in (P3).This is to be expected given the very small size of problems and the reduced number of bilinear terms, which facilitates the spatial branch-and-bound procedure in BARON.In order to evaluate how both methods scale with problem size, we solved the five most difficult (higher optimality gap at termination when solved with BARON) water-using network design problems in Teles et al. [11], which also correspond to bilinear programs.In such problems, the discretized variables are concentrations featuring P = 3 in Ex17 and P = 2 in the other examples.The relative optimality tolerance was set to 0.01 % and the total computational time to 1-h.Table 5 lists the results for global optimization solvers BARON and GloMIQO together with the implementation of algorithm 2 for MDT.While BARON fails to prove optimality in all five cases, GloMIQO and MDT are successful in two, which interestingly are not the same.GloMIQO takes roughly half an hour to tackle Ex14 and Ex15, while it takes just a few seconds for MDT to solve EX17 and Ex18.In the full set of problems from [11], MDT performs better than BARON 84 % of the time and slightly better than GloMIQO (53 %) when considering the CPUs performance metric.The average optimality gaps are however considerably lower, 0.088 % for MDT, 0.466 % for GloMIQO and 0.711 % for BARON.
In Table 6, several blending problems [45,46] are solved using BARON, GloMIQO, and Algorithm 3 at a single level of discretization.These are multiperiod blending problems with varying numbers of tanks, time periods, and product qualities.Each problem is identified such that 6T-3P-2Q-029 is a 6 tank, 3 time period, 2 quality problem, with a unique 3-digit identifier to distinguish it from other problems of the same size.The results for MDT are reported for P = 0 and p = −3 and were solved using Gurobi [47] using 12 threads.As these problems are originally MINLPs, larger than those in Table 5, and 12 threads are utilized, the computational difference is much more significant.GloMIQO and BARON both use MILP solvers in their algorithms which can utilize multiple threads, but this effect is not significant as the subproblems being solved are generally very small.Because of the use of multiple threads and multiple CPUs in this comparison, wall times are reported for all three methods.While GloMIQO and BARON are unable to converge to an optimality gap of 0.1 % within the time limit of 2 h, the MDT is able to close the gap in all but two cases.In these two cases an additional level of discretization would close the gap, but these results show that if a reasonably fine discretization is chosen a priori, the optimality gaps can be significantly less than those of commercial global optimization solvers.For the one problem that BARON and GloMIQO were able to close the gap within the time limit, MDT outperformed them by approximately 17 and 1,400 times, respectively, although the problem is small enough and is solved fast enough that meaningful conclusions should not be drawn from this result alone.

Table 6
Comparison for multiperiod blending problems [45,46] Problem This paper has presented a derivation based on disjunctive programming and exact linearization of the multiparametric disaggregation technique (MDT) for solving bilinear programming problems.The derivation, which was shown to be equivalent to applying the reformulation linearization technique (RLT), gives rise to an MILP approximation that yields an upper bound if feasible with respect to the original problem.A lower bounding MILP relaxation problem has also been derived, which can then be used as basis for rigorous global optimization algorithms.As has been shown with the smaller test problems, the relaxation from the MDT was shown to be as tight as the relaxation based on piecewise McCormick envelopes (PCM) for strictly bilinear terms but weaker for quadratic terms.The real advantage of MDT comes from the more favorable scaling of problem size as discretization is increased, which is translated into an ability to reach considerably lower optimality gaps and a clearly better computational performance.For large problems it was shown that multiparametric disaggregation can outperform global solvers BARON and GloMIQO.

Fig. 1 Fig. 2
Fig. 1 Feasible region for a bilinear curve x i • x j = 0.1, p = −2, P = −1.The solid curve represents the exact bilinear curve, while the dots represent the approximated (and incomplete) feasible region resulting from the upper bounding formulation

Fig. 3 Fig. 4
Fig. 3 Plot of the lower bounding feasible region for w i j = 0.1 and p = P = −1.The solid line is the true curve, while the dotted lines represent the boundaries of the relaxed feasible region of (PR).Note the similarities to piecewise McCormick envelopes

Fig. 5 1 Fig. 6
Fig. 5 Feasible region from conventional McCormick, piecewise McCormick (PCM) and multiparametric disaggregation relaxation (MDT) for a quadratic term over a specific region.Notice that the underestimation from PCM is strictly tighter than MDT in the region x i ∈]25, 35[\{30} time (s) Gap (%) Wall time (s) Gap (%Italic values indicate that the maximum computational time was reached 8 Conclusions