Adaptive piecewise linear relaxations for enclosure computations for nonconvex multiobjective mixed-integer quadratically constrained programs
Abstract
EconStor is a publication server for scholarly economic literature, provided as a non-commercial public service by the ZBW.
Full text
Link, Moritz; Volkwein, Stefan Article — Published Version Adaptive piecewise linear relaxations for enclosure computations for nonconvex multiobjective mixed-integer quadratically constrained programs Journal of Global Optimization Provided in Cooperation with: Springer Nature Suggested Citation: Link, Moritz; Volkwein, Stefan (2023) : Adaptive piecewise linear relaxations for enclosure computations for nonconvex multiobjective mixed-integer quadratically constrained programs, Journal of Global Optimization, ISSN 1573-2916, Springer US, New York, NY, Vol. 87, Iss. 1, pp. 97-132, https://doi.org/10.1007/s10898-023-01309-5 This Version is available at: https://hdl.handle.net/10419/307034 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/
Journal of Global Optimization (2023) 87:97–132 https://doi.org/10.1007/s10898-023-01309-5 Adaptive piecewise linear relaxations for enclosure computations for nonconvex multiobjective mixed-integer quadratically constrained programs Moritz Link1 ·Stefan Volkwein1 Received: 29 August 2022 / Accepted: 20 June 2023 / Published online: 5 July 2023 © The Author(s) 2023 Abstract In this paper, a new method for computing an enclosure of the nondominated set of multiobjective mixed-integer quadratically constrained programs without any convexity requirements is presented. In fact, our criterion space method makes use of piecewise linear relaxations in order to bypass the nonconvexity of the original problem. The method chooses adaptively which level of relaxation is needed in which parts of the image space. Furthermore, it is guaranteed that after finitely many iterations, an enclosure of the nondominated set of prescribed quality is returned. We demonstrate the advantages of this approach by applying it to multiobjective energy supply network problems. Keywords Mixed-integer nonlinear programming ·Multiobjective optimization ·Box enclosure ·Adaptive piecewise linear relaxation ·Energy supply networks Mathematics Subject Classification 65K05 ·90C11 ·90C26 ·90C29 ·90C35 1 Introduction Multiobjective optimization problems arise if there is more than one objective of interest and the objectives are in conflict with each other, i.e. in general it is not possible to obtain a solution that is optimal w.r.t. each of the objective functions simultaneously. Therefore, the nondominated set (cf. Definition 2.1) consists of so-called optimal compromises.Tohavean overview of (at least some of) these optimal compromises is an advantage as one can observe the interplay between the objective functions in a more direct and easy way. Especially in real-world applications, it occurs rather rarely that only one objective is of interest. This is also the case in our application, where we are looking for energy supply network plans which BMoritz Link [email protected] Stefan Volkwein stefan.v[email protected] 1Department of Mathematics and Statistics, University of Konstanz, Universitätsstraße 10, 78464 Constance, Baden-Württemberg, Germany 123
98 Journal of Global Optimization (2023) 87:97–132 are optimal w.r.t. the occurring costs as well as its accompanying CO2-emissions. Naturally, a network plan being low in emissions is more costly and vice versa. Thus, having a list of optimal compromises between these two extreme solutions helps the political decisionmakers to get a broader picture of the situation. In most cases, solving a multiobjective optimization problem (MOP) numerically means that one computes either an enclosure (cf. [23–26]) or an approximation (cf. [5,10,17,30, 33]) of the nondominated set. This is due to the fact that the nondominated set is in general infinite. In the linear case, it is also possible to compute the entire nondominated set (cf., e.g., [47]). In this paper, we introduce a deterministic approach to compute an enclosure of the nondominated set with a prescribed quality ε>0. Although the focus of our algorithm is to compute such an enclosure, it also returns an approximation of the nondominated set, namely a finite set of ε-nondominated points of (MOP) (see Definition 2.3). In this work, we present two new methods to compute such enclosures. The first method (see Algorithm 2) tackles the problem more directly without being too complicated. In fact, it is a slight modification of the method presented in [23] in the way that we use the Restricted Weighted Sum Method (cf. [31,32]) to obtain nondominated points of the problem. While being very effective for problems where single objective subproblems can be solved very quickly, it is not applicable (in terms of efficiency) to more complex problems—which can happen even with a comparably small number of variables in the nonconvex mixedinteger setting. Furthermore, this approach depends heavily on the availability of solvers for the resulting single objective problems. These drawbacks are addressed with the second method, where we combine the direct approach with piecewise linear relaxations to bypass the nonlinearities of the original problems—and therefore only single objective mixed-integer linear programming (MILP) problems have to be solved, where powerful solvers are available (cf. [14,29]). The accuracy of these relaxations is chosen adaptively by the method itself during the solution process. Consequently, it is possible to use different levels of accuracy in different parts of the criterion space, which can be of importance if also the relaxed problems are difficult to solve. This approach is the first of its kind for a (MOP). In [24] a method is presented to solve multiobjective convex mixed-integer optimization problems via relaxations. It makes use of cuts in the criterion space in order to discard infeasible integer assignments as well as to separate parts that do not belong to the image of the feasible set. In [25] the authors present the first deterministic method being able to deal with multiobjective mixed-integer nonconvex problems. However, the method presented there is based on a branch-and-bound scheme in the decision space which is a fundamental difference from the methods presented in this paper. Note further that due to clarity and conciseness, we restrict ourselves to nonconvex quadratically constrained problems in this paper as this is of relevance to our application. However, the presented method can be straightforwardly generalized to polynomially constrained programs using reduction techniques to decompose polynomials into quadratic and bilinear terms (cf., e.g., [45] and Remark 4.3). The manuscript is organized as follows: After introducing preliminary notions in Sect. 2 we define enclosures and related notions as well as provide a basic approach for computing such an enclosure in Sect.3. In Sect.4we introduce and relate piecewise linear relaxations of (MOP) to so-called lower bound sets of the nondominated set, which results in a new procedure for computing an enclosure using adaptively chosen relaxations presented in Sect.5.We further prove correct and finite termination of this procedure. In Sect. 6we demonstrate the advantages of using this novel approach by applying it to different instances of an energy supply network optimization problem, which are in fact nonconvex multiobjective mixed-integer quadratically constrained problems. Finally, a conclusion is drawn in Sect.7. 123
Journal of Global Optimization (2023) 87:97–132 99 2 Notations and definitions Throughout we write x≤yfor x,y∈Rkif for any i∈[k]:={1,...,k}we have that xi≤yi. In the same manner we write x<yif xi<yifor any i∈[k]. The ingredients of a general multiobjective optimization problem are the continuous components fi:Rn+m→R,i∈[k], of the objective function f=(f1,..., fk):Rn+m→Rk, the continuous components hi:Rn+m→R,i∈[p], of the equality constraint function h=(h1,...,hp):Rn+m→Rpand the continuous components gi:Rn+m→R,i∈[q], of the inequality constraint function g=(g1,...,gq):Rn+m→Rq. For all of these functions, we only allow quadratic and bilinear terms. The corresponding multiobjective optimization problem is then given by min f(x)s.t. h(x)=0,g(x)≤0,x∈X:= XC×XI,(MOP) where XC=[C,uC]R nand XI=[I,uI]∩ZmR mare non-empty boxes1with C,uC∈Rn, C<uC,and I,uI∈Zm, I<uI, respectively. These boxes represent the box constraints of the continuous and integer variables, respectively. If m=0, i.e. if there are no integrality constraints for any variable, we call (MOP)acontinuous multiobjective optimization problem, and if n=0 we call it a pure integer multiobjective optimization problem. Note that the components of the occurring functions f,g,andhare neither assumed to be linear nor convex. For a survey on solution methods for general versions of (MOP)we refer the reader to [22]. Single objective versions of (MOP) we denote by MIQCP (or more general MINLP) problems. Interested readers are referred to [1,6,9,36,43,52]. As we may require integrality of some variables and due to the possible non-convexity of the occurring functions of (MOP) we obtain its possibly nonconvex feasible set S={x∈X|h(x)=0,g(x)≤0}, which we assume to be nonempty, as well as its possibly nonconvex image set f(S)=f(x)∈Rk|x∈S. We define the projection of the feasible set on Rmby SI=xI∈Zm|∃xC∈Rn:(xC,xI)∈S, which we call the set of feasible integer assignments of (MOP). For a given integer assignment ˆxI∈Zmwe define the sets SˆxI=x∈S|xI=ˆxIand XˆxI=x∈X|xI=ˆxI, describing the feasible and box-feasible solutions corresponding to ˆxI, respectively. Notice that SˆxI=∅holds, provided that ˆxI∈Zm\SI. Since we assume all constraint functions to be continuous and by the definition of X, the feasible set Sis compact. Furthermore, by the continuity of the objective functions, we know that f(S)is also compact and we can therefore find a box B=[z,zu]⊆Rkwith f(S)⊆int(B)for some z,zu∈Rk,z≤zu,z= zu. In general, the kobjective functions of (MOP) are in conflict with each other, i.e. we cannot assume that there is a solution ¯x∈Sminimizing all objective functions simultaneously. This motivates the concepts of (non-)dominance and efficiency (cf. [20]). 1A closed n-dimensional box is defined as [, u]:=({}+Rn +)∩({u}−Rn +)={x∈Rn|≤x≤u}, where , u∈Rn,l<u. The open box is defined as (, u):= ({}+int(Rn +)) ∩({u}−int(Rn +)) ={x∈Rn| <x<u}. 123
100 Journal of Global Optimization (2023) 87:97–132 Definition 2.1 (1) Let y1,y2∈Rk.Theny1dominates y2(and y2is dominated by y1)if y1 i≤y2 ifor all i∈[k]and y1 j<y2 jfor some j∈[k]. We write y1≤y2,y1= y2.Ifeveny1 i<y2 ifor all i∈[k],theny1strictly dominates y2(and y2is strictly dominated by y1) and we write y1<y2. (2) Let y∈Rkand N⊂Rk. Then the vector yis (strictly) dominated by the set Nif there exists a vector ˆy∈N(strictly) dominating y. Similarly, if yis not (strictly) dominated by Nwe call y (weakly) nondominated by N. (3) A point ¯x∈Sis called an efficient solution of (MOP)ifthereexistsnox∈Ssuch that f(x)dominates f(¯x). If there exists no x∈Ssuch that f(x)strictly dominates f(¯x),then ¯xis called a weakly efficient solution of (MOP). The set of (weakly) efficient solutions of (MOP) is called its (weakly) efficient set. (4) A point ¯y∈f(S)is called a (weakly) nondominated point of (MOP) if there exists no y∈f(S)(strictly) dominating ¯y. The set of (weakly) nondominated points of (MOP)is called its (weakly) nondominated set and is denoted by N(and Nw). Remark 2.2 (1) Note that N=f(x)∈Rk|x∈Sis efficient. (2) Note that even if we assume all objective and constraint functions of (MOP)tobeconvex (or linear), the nondominated set of (MOP) can be disconnected and nonconvex due to integrality constraints (e.g., c.f. [16, Figure 1]). The goal of multiobjective optimization is to determine the nondominated set Nof a given (MOP). As the set Ncan be infinite or of a very complicated structure it can be impossible to compute it exhaustively using numerical methods. Therefore, numerical methods for multiobjective optimization mostly focus on computing a finite approximation Aof Nsubsequent to the following three goals (cf. [49]): •Coverage all parts of the nondominated set Nshould be represented in the approximation A. •Uniformity the points of the approximation Ashould be distributed uniformly along the nondominated set N. •Cardinality the set Ashould contain an appropriate number of points naturally depending on the number and extent of the cost functions. Another approach of numerically solving an (MOP) is—instead of computing a sole approximation—to compute a coverage of the nondominated set Nas presented in, e.g., [23–26]. Besides the introduction of enclosures of the nondominated set, the authors present a branch-and-bound framework for computing such an enclosure of the nondominated set with a prescribed quality. As we are aiming at the same goal in this paper, we present the enclosure-related notions in more detail in Sect.3. However, speaking of enclosures of a certain quality we come to the issue of obtaining convergence of numerical methods mostly only in the limit, i.e. after possibly infinitely many iterations. In general, this is overcome by using termination criteria together with specific tolerances—in fact, one terminates a method if it computed a solution satisfying an a-priori chosen criterion up to a specified tolerance. In single objective optimization, one way to use such tolerances is embodied in branch-and-bound methods which iteratively compute lower and upper bounds of the optimal value while checking if the gap between these two becomes smaller than the prescribed tolerance. In particular, these methods terminate with so-called ε-optimal solutions. This concept can also be applied to nondominance in multiobjective optimization. 123
Journal of Global Optimization (2023) 87:97–132 101 Definition 2.3 Let ε>0. A point ¯y∈f(S)is called a (weakly) ε-non-do-mi-na-ted point of (MOP)ifthereexistsnoy∈f(S)such that y+εe(strictly) dominates ¯y,wheree∈Rk denotes the all-one vector. Similarly, a solution ¯x∈Sis called (weakly) ε-efficient for (MOP) if there is no x∈Ssuch that f(x)+εe(strictly) dominates f(¯x).Thesetofε-nondominated points is denoted by Nε. Remark 2.4 Depending on the considered problem (MOP) one could also use another vector instead of e. E.g., this could be useful if the attained magnitudes of the objective functions differ largely. For an example see the explanations of the numerical examples in Sects.3 and6. 3 Enclosures and local bounds In [23,24,26] the authors present methods for solving problems like (continuous or convex versions of) (MOP). The aim of these methods is to find an enclosure of the nondominated set N. Following the sandwiching idea of single objective branch-and-bound methods, an extension to the multiobjective setting is needed. Definition 3.1 Let L,U⊆Rkbe two finite sets satisfying N⊆L+Rk +and N⊆U−Rk +. Then Lis called a lower bound set,Uis called an upper bound set,andthesetE:= E(L,U) defined as N⊆E(L,U):= L+Rk +∩U−Rk += ∈L u∈U, ≤u [, u] is called enclosure of the nondominated set Nof (MOP). For transferring the gap termination criterion to the multiobjective setting we need some sort of quality measure for the set E(L,U).In[23–26] the authors use the so-called width w(E)which is defined as the optimal value of the problem max ,us(, u)s.t. ∈L,u∈U,≤u,(W(E)) where s(, u):= min{ui−i|i∈[k]} represents the shortest edge length of a given box [, u]. The justification for the choice of this width measure is given by the following lemma. Lemma 3.2 (cf. [26, Lemma 3.1]) For sets L ,U⊆Rkwith N⊆L+Rk +and some ε>0 let w(E)<ε. Then the relation E(L,U)∩f(S)⊆Nε(1) holds, i.e. for any feasible point x ∈S with f (x)∈E we know that f (x)is an ε-nondominated point of (MOP). Due to Lemma 3.2, finding lower and upper bound sets L,Uand thus an enclosure Ewith w(E)<εbecomes important. In fact, computing such an enclosure is the aim of the methods from [23,24], whereas in the branch-and-bound based method from [26] the authors use the width measure in the preimage space to declare a box sufficiently branched. As we propose a criterion space method, we focus on computing an enclosure Esuch that w(E)<ε.In order to do so, we make use (as done in [23–25]) of the concept of Local Lower and Local 123
102 Journal of Global Optimization (2023) 87:97–132 Upper Bounds (LLBs and LUBs, respectively) which is an extension of the work in [15,34] andalso[21]. In [15] the authors call a set Y⊆Rkstable if there are no elements in Ydominating each other, i.e. for any two distinct y1,y2∈Ythere are i,j∈[k]such that y1 i<y2 iand y2 j<y1 j. One can immediately see that for any (MOP) the corresponding nondominated set Nis stable. We start with the definition of local upper bounds as introduced in [34]. Definition 3.3 Let N⊆f(S)be a finite and stable set. Then the lower search region for N is s(N):= y∈int (B)|yyfor every y∈N, and the lower search zone for some u∈Rkis given by c(u):= {y∈int (B)|y<u}. AsetU=U(N)is called local upper bound set given Nif (i) s(N)=u∈U(N)c(u), (ii) cu1cu2for any u1,u2∈U(N)with u1= u2. Each point u∈U(N)is called a local upper bound (LUB). In [23,24] the authors extended the concept of LUBs to so-called Local Lower Bounds. Definition 3.4 Let N⊆int (B)be a finite and stable set. Then the upper search region for Nis S(N):= y∈int (B)|yyfor every y∈N, and the upper search zone for some ∈Rkis given by C() := {y∈int (B)|<y}. AsetL=L(N)is called a local lower bound set given Nif (i) S(N)=∈L(N)C(), (ii) C1C2for any 1, 2∈L(N)with 1= 2. Each point ∈L(N)is called a local lower bound (LLB). In [26] the authors show that given a stable set N, the corresponding local upper bound set is uniquely determined. In the same manner, one can show that also the corresponding local lower bound set is unique. The following result from [23] relates the concept of local lower/upper bounds to lower/upper bounds as introduced in Definition 3.1. Lemma 3.5 (cf. [23, Corollary 3.6]) Suppose that the sets N1⊆f(S)and N2⊆ int(B)\f(S)+int(Rk +)are finite and stable. Then U(N1)is an upper bound set and L(N2)is a lower bound set in the sense of Definition 3.1. As the local lower and upper bound sets depend on the stable set N, we have to update both if we add a new point yto the set N. These procedures come from [34, Algorithm 3], [23, Algorithm 2]. In Algorithm 1we provide the scheme for updating a local upper bound set. Similar to [34] we use the following notation: for y∈Rk,c∈Rand i∈[k]we define y−i:= (y1,...,yi−1,yi+1,...,yk)and 123
Journal of Global Optimization (2023) 87:97–132 103 Algorithm 1 Updating a Local Upper Bound Set Require: Local upper bound set U(N)and update point y∈f(S) 1: A={u∈U(N)|y<u} 2: for i∈[k]do 3: Bi=u∈U(N)|yi=uiand y−i<u−i 4: Pi=∅ 5: end for 6: for i∈[k]do 7: for u∈Ado 8: Pi=Pi∪yi,u−i 9: end for 10: end for 11: for i∈[k]do 12: Pi=u∈Pi|uufor all u∈Pi∪Bi,u= u 13: end for 14: U(N∪{y})=(U(N)\A)∪i∈[k]Pi return Updated local upper bound set U(N∪{y}) (c,y−i):= (y1,...,yi−1,c,yi+1,...,yk). In order to obtain an algorithm for updating a local lower bound set L(N)w.r.t. an update point y∈int(B), one can modify Algorithm 1in the following way. We replace U(N)by L(N), and change the following: •Step 1 to A={∈L(N)|y> }, •Step 3 to Bi={∈L(N)|yi=iand y−i> −i}, •Step 12 to Pi=∈Pi|for all ∈Pi∪Bi, = . Furthermore, by analyzing the update procedure of U(N)w.r.t. a point y∈f(S), one can relate the local upper bounds from U(N)to the ones from U(N∪{y})in the following way: for a local upper bound u∈U(N∪{y})we have either u∈U(N)or u=yi,u −i for some i∈[k]and u∈U(N). For the latter case, we call uthe parent of u. Otherwise, u is its own parent. Lemma 3.6 [23, Lemma 3.8] Let u ∈U(N∪{y})be a local upper bound. Then its parent u∈U(N)is unique. Similarly, we define the parents of a given local lower bound . Furthermore, the analogue of Lemma 3.6 holds also for any local lower bound ∈L(N∪{y}). In [15,34] the authors present a general scheme for computing an approximation of the nondominated set of a multiobjective optimization problem (cf., e.g., [15, Algorithm 2]). In [23] the authors extend this scheme for computing an enclosure Eof Nwhich, in fact, consists of boxes [, u],whereis a local lower and ua local upper bound, and satisfies w(E)<εfor given ε>0. Both methods make use of the search zones determined by a given local upper bound u. For this paper, we use a slightly modified approach of these as a basic scheme which is written in Algorithm 2. In Algorithm 2, we loop through the set Uloop, i.e. the set of local upper bounds at the beginning of the iteration. If the current local upper bound uis part of a not sufficiently small box (see Step 5), we try to improve this local upper bound, i.e. find a nondominated point in the search region determined by u. In order to do so, we solve the weighted-sum problem min αf(x)s.t. x∈Sand f(x)<u−δe,(WSP (α;u)) 123
104 Journal of Global Optimization (2023) 87:97–132 Algorithm 2 General scheme for computing an enclosure of the nondominated set relying on a scalarization technique. Require: box B=z,zuwith f(S)⊆int(B), termination tolerance ε>0 and off-set factor δ>0 1: Initialize nondominated set N=∅, set of local lower bounds L=zand set of local upper bounds U=zu 2: while w(E)≥εdo 3: Uloop =U 4: for u∈Uloop do 5: if there exists ∈Lwith ≤uand s(, u)≥εthen 6: if there exists y∈Nwith y<u−δethen 7: Update Uw.r.t. yusing Algorithm 1 8: Update Lw.r.t. yusing modified Algorithm 1 9: else 10: Update Lw.r.t. u−δeusing modified Algorithm 1 11: end if 12: end if 13: end for 14: end while return Enclosure E(L,U)satisfying w(E)<ε for some weight vector α∈int(Rk +). Note that if the weight vector αhas only strictly positive entries we know that any global solution ¯xof (WSP (α;u)) is efficient for (MOP) (cf. [4, Lemma 1.5.2],[20, p. 214f]), i.e. ¯y=f(¯x)∈Nand ¯y<u−δe. Furthermore, we know that there exists y∈Nwith y<u−δeif and only if (WSP (α;u)) has a solution. In sum, this means that we can decide whether to enter and execute the if-statement in Step 6 after solving (WSP (α;u)) for a strictly positive weight vector α. This weight vector could be chosen always the same, e.g., α=(1,...,1)/k. But one could also compute it individually for any u, as done for the numerical examples presented later in this work. For example, one could compute the weight vector αfor the problem (WSP (α;u)) using a computation technique proposed in [48]. Given the current local upper bound of interest u∈U(N),where Nis the current approximation of N,foranyi∈[k]let yi∈Nbe a defining point of the i–th component of u, i.e. for any i∈[k]we have yi i=uiand yi −i<u−i. We then solve the system ⎛ ⎜ ⎝ y1 1... y1 k . . .. . . yk 1... yk k ⎞ ⎟ ⎠⎛ ⎜ ⎝ ˜α1 . . . ˜αk ⎞ ⎟ ⎠=⎛ ⎜ ⎝ 1 . . . 1 ⎞ ⎟ ⎠,(2) and set α:= |˜α|/˜α2.Notethatαis the normal vector of the hyperplane determined by the points y1,...,yk.Using(WSP (α;u)), Algorithm 2can be seen in the context of the so-called Adaptive Weighted Sum Method as presented in [31,32]. Remark 3.7 For a clear definition of defining points we refer to [34]. Furthermore, in [34] the authors present methods for updating local upper bound sets together with the sets of defining points, i.e. an algorithm doing the same as Algorithm 1, but also updating the corresponding defining points (see [34, Algorithm 5]). In fact, [34, Algorithm 5] is used in any implementation of the present paper instead of Algorithm 1.Wewaivedtopresent[34, Algorithm 5] in detail to keep things more concise. If (WSP (α;u)) is feasible, it is solved to optimality and its solution ¯y∈Nis used to update the local upper and lower bound sets. Afterward, it proceeds to the next iteration. If 123
Journal of Global Optimization (2023) 87:97–132 111 4 Piecewise linear relaxations and lower bounding As we have seen, using Algorithm 2we are able to compute an enclosure of the nondominated set of a given (MOP). However, as all the weighted-sum problems that have to be solved during Algorithm 2are still MINLP problems in general, there may arise inconvenience while tackling, e.g., larger instances (see Sect.6). One idea to cope with this issue is to bypass the non-linearity of the problems, e.g., trying to solve just MILP problems. We have seen before that the idea of computing enclosures instead of a finite approximation of the nondominated set is somehow lifted from the single objective setting to the multiobjective one. In the same spirit, we are going to use an idea coming from single objective optimization in regard to solving MINLP problems. In [13,40,45] several methods for solving MINLP problems, e.g., a (MOP) with k=1, by iteratively solving adaptively refined convex or linear MIP relaxations of the original problem are presented. Note that there is plenty of literature and ongoing research in the area of efficiently computing linear or convex relaxations (or approximations) of nonlinear problems (cf., e.g., [6,12,13,28, 38,44,53]). Although the basic idea—but not the realization—of consequently refining the relaxations and therefore guaranteeing convergence of the sequence of optimal values of the relaxed problems to the optimal value of the original problem is the same, the respective termination criteria of the algorithms from [45] on the one hand and [13,40] on the other hand differ. In the latter, both algorithms terminate if the maximal possible constraint violation of the relaxed problem in comparison to the original problem is below an a-priori given, but arbitrarily small error bound. Morally speaking, we say that the original MINLP problem is solved to optimality if we computed an optimal solution of a relaxation violating the nonlinear constraints only in an acceptable manner. The authors prove that under certain assumptions on the chosen relaxation technique, the proposed methods terminate after a finite number of steps. In [45] the authors terminate their algorithm by the gap criterion which is known from classical branch-and-bound methods. The algorithm uses parts of the optimal solutions of the relaxed problems to fix variable values in the original problem, e.g., one can fix all integer variables in the MINLP problem to take the corresponding values of the relaxed solution and obtain an NLP problem which—if it is feasible—might be easier to solve and then gives an upper bound on the optimal value. In case of fixing the integer variables of the relaxed solution (˜xC,˜xI)in the MINLP problem the resulting NLP problem has the following form min f(x)s.t. x∈S˜xI.(redMOP(˜xI)) Furthermore, as before, the lower bound on the optimal value is consequently improved by refining the relaxations. In every step, the smallest available upper bound and the best lower bound is used to sandwich the optimal value of the original MINLP problem and terminate the algorithm if the gap between these bounds is small enough. However, it is not obvious that the resulting NLP problems become feasible at any time during the procedure, and therefore there is no guarantee that any upper bound is available for the gap criterion at all. Nevertheless, in practice, it is rather unlikely to produce exclusively infeasible integer assignments with the solution of the relaxed problems, in particular with increasing accuracy. In this paper, we want to merge both ideas as the use of upper bounds together with the gap criterion may cause termination even if the maximal possible constraint violation of the relaxation is not below the tolerance and therefore accelerate the computations. In fact, we use the scheme presented in [13] to preserve the theoretically ensured convergence of the method, but also add the upper bounding part as well as the gap termination criterion from 123
112 Journal of Global Optimization (2023) 87:97–132 [45] to avoid unnecessary computations. Furthermore, the relaxation techniques used in [13] are able to handle general nonlinear functions, whereas the techniques in [45]arebasedon the so-called McCormick relaxations for bilinear and multilinear terms (cf. [42]) and are therefore only able to handle general polynomial terms. For this paper, we restrict ourselves to the explicit handling of quadratic and bilinear terms as these are the only ones appearing in our application (see Sect.6). However, the method presented in the remainder of this paper is capable of handling general nonlinearities if an adequate relaxation technique is available. We start by introducing the necessary notions. Definition 4.1 Let (MOP) be given. Further let ˜ S⊆Rn+mbe a set and ˜ f:Rn+m→Rka k-vector-valued function such that S⊆˜ Sand ˜ f(x)≤f(x)for any x∈S.Thenwecall min ˜ f(x)s.t. x∈˜ S,(RMOP ˜ S) arelaxation of (MOP). We denote the nondominated set of (RMOP ˜ S)byN˜ S. Note that, if k=1, we have that the optimal value ˜ f(˜x∗)of (RMOP ˜ S) is a lower bound on the optimal value f(x∗)of (MOP). For brevity we assume without loss of generality that the objective fin (MOP) is linear. This is also motivated by our energy supply network considered in Sect.6Consequently, fneeds no further relaxation and for the remainder of that paper a relaxation (RMOP ˜ S) is characterized by its feasible set ˜ S. Given two different relaxations ˜ Sand ˆ S, we call the relaxation ˜ S finer than ˆ Sif ˜ S⊆ˆ S.Again,ifk=1, for two relaxations ˜ S⊆ˆ Swe have that the optimal value f(ˆx∗)of the relaxation ˆ Sis a lower bound on the optimal value f(˜x∗)of the relaxation ˜ S. There is a plenty of possibilities to obtain relaxations of a given (MOP). One way is to use the concept of piecewise linear underand overestimators (cf., e.g., [11,13,28]). Definition 4.2 Let h(x1,...,xn)be a nonlinear real-valued function with compact domain Dh⊂Rnappearing in (MOP), e.g., as a constraint function. •We call a continuous piecewise linear function hu:Dh→Rapiecewise linear underestimator of hif hu(x)≤h(x)for all x∈Dh. •We call a continuous piecewise linear function ho:Dh→Rapiecewise linear overestimator of hif h(x)≤ho(x)for all x∈Dh. •We call Rh=(hu,ho),wherehuis a piecewise linear underestimator and hoa piecewise linear overestimator of h,apiecewise linear relaxation of h. Given a piecewise linear relaxation Rhof a term h, we obtain a relaxation in the sense of Definition 4.1 by replacing any appearance of hby an additional variable ˆ hand adding the constraints hu(x)≤ˆ hand ˆ h≤ho(x). Indeed, this yields a relaxation in the sense of Definition 4.1 by the fact that (h(x), x)∈R1+n|x∈Dh ⊆ˆ h,x∈R1+n|hu(x)≤ˆ h≤ho(x)for all x∈Dh. Thus, given (MOP) with nonlinear terms hi,i∈I, for some finite index set Itogether with piecewise linear relaxations Rhi,i∈I, and proceeding as above, we obtain a relaxation of (MOP) with feasible set denoted by ˜ SRor simply by R,whereRis the collection of all Rhi. For the remainder of that paper we only consider relaxations coming from piecewise linear relaxations of all appearing nonlinear terms of (MOP). One can see immediately that given two relaxations Rand Rwe have that R⊆Rif and only if for any nonlinear term 123
Journal of Global Optimization (2023) 87:97–132 113 happearing in (MOP)wehaveR h⊆Rh, i.e. hu(x)≤h u(x)and h o(x)≤ho(x)for all x∈Dh. Consequently, by improving the piecewise linear underand overestimators we refine the corresponding relaxation. Accordingly with [11], we measure the quality of a given piecewise linear relaxation Rhof a nonlinear term hby •the overestimation error ε(h,Rh) o:= max x∈Dh ho(x)−h(x), •the underestimation error ε(h,Rh) u:= max x∈Dh h(x)−hu(x), •and the overall relaxation error ε(h,Rh) rel := max ε(h,Rh) u,ε (h,Rh) o. Clearly, the relaxation error decreases while refining relaxations, i.e. for two relaxations R h⊆Rhwe have that ε(h,R h) rel ≤ε(h,Rh) rel . Naturally, we can also extend this concept to measure the quality of a relaxation Rof a given (MOP). In fact, the overestimation error of Ris defined by εR o:= max εhi,Rhi o|hiappears in (MOP)and Rhi∈Rfor i∈I. Similar for the underestimation error.Theoverall relaxation error is then given by εR rel := max{εR o,ε R u}. Similar as before, if R⊆Rfor two relaxations Rand Rwe have that εR rel ≤εR rel. Let us now briefly introduce the relaxation technique of interest for that paper, namely the so-called McCormick piecewise linear relaxations (cf. [42]). Remark 4.3 As already mentioned, we put the focus in this paper on the ingredients necessary for our application. As there appear only bilinear and quadratic constraints, we restrict ourselves to the description of exactly this case. However, by decomposing general polynomials into bilinear and quadratic terms one can apply the method from this work to general polynomially constrained problems. For example, the polynomial x3ycan be decomposed into x1=x2, x2=xy, x3=x1x2, where xi,i=1,2,3, are additional variables. We start with the bilinear case, i.e. h(x,y)=xy, with compact domain Dh=(x,y)∈R2|x≤x≤xu,y≤y≤yu. The McCormick piecewise linear relaxation depends on a so-called partition (or triangulation) of the domain Dh. For our purpose it is enough to consider simple partitions of the intervals [x,xu]and [y,yu]in requidistant intervals. Consequently, we obtain gridpoints xi= x+i(yu−x)/rand yi=y+i(yu−y)/rfor i∈{0,...,r}with corresponding intervals [xi−1,xi]and [yi−1,yi]for i∈[r]. 123
114 Journal of Global Optimization (2023) 87:97–132 Fig. 4 McCormick relaxations of the bilinear term xyusing one partition per variable (green) and two partitions per variable (red). (Color figure online) For any i,j∈[r]and any (x,y)∈[xi−1,xi]×[yj−1,yj]the linear underestimator function on Dij h:= [xi−1,xi]×[yj−1,yj]is given by hij u(x,y):= max xi−1y+yj−1x−xi−1yj−1,xiy+yjx−xiyj, and the linear overestimator is given by hij o(x,y):= min xi−1y+yjx−xi−1yj,xiy+yj−1x−xiyj−1. The piecewise linear underand overestimator functions huand hoare then given by concatenation as shown in Fig. 4. Note that huand hoare continuous piecewise linear functions and hu(x,y)≤h(x,y)≤ ho(x,y)for all (x,y)∈Dh, i.e. huis a piecewise linear underestimator and hois a piecewise linear overestimator of hin the sense of Definition 4.2. We denote the piecewise linear relaxation based on requidistant partitions for each variable by Rr xy. The resulting relaxation error is given by ε(h,Rr h) rel =xu−xyu−y 4r2. Although the quadratic case is just a special case of the bilinear case we give the explicit formulation. Therefore, let h(x)=x2with Dh={x∈R|x≤x≤xu}. We consider again the partition in requidistant intervals, i.e. Di h=[xi−1,xi]for i∈[r].Foranyi∈[r] and any x∈Di hthe linear underestimator is given by hi u(x):= max 2xi−1x−x2 i−1,2xix−x2 i, and linear overestimator is given by hi o(x):= (xi−1+xi)x−xi−1xi. 123
Journal of Global Optimization (2023) 87:97–132 115 The resulting relaxation error corresponding to requidistant partitions is ε(h,Rr h) rel =xu−x2 4r2. Clearly, for the McCormick piecewise linear relaxations we have that ε(h,Rr h) rel →0for r→∞, i.e. using the McCormick piecewise linear relaxation technique, for any ε>0we find r∈Nsuch that ε(h,Rr h) rel <ε. Remark 4.4 In [45] the authors provide an adaptive partitioning scheme, i.e. there is no requirement for equidistant intervals. New breakpoints are added in parts of the variable domains close to the value of the previous solution in order to refine the relaxation. The idea is to only partition on regions of the variable domain that appear to influence optimality the most. In contrast, using uniformly distributed break points, one runs the risk of refining in uninteresting regions of the variable domain and therefore unnecessarily blowing up the problem. The scheme presented in [45] could also be incorporated2into the methods of the present paper. But for simplicity of presentation and since the equidistant procedure works very well for our application, we restrict ourselves to the equidistant partitioning scheme presented in this paper. After we have introduced the relaxations of interest for the present paper, we now relate relaxations (RMOP ˜ S) and lower bounds of the nondominated set Nof (MOP). For the ease of notation, we subsequently assume that for a given relaxation ˜ Swe have that f(˜ S)⊆int(B). We obtain the following lemma. Lemma 4.5 Let f (˜x)∈N˜ Sfor some ˜x∈˜ S and some relaxation ˜ Sof(MOP).Then f(˜x)is nondominated by N,i.e. f(˜x)∈int(B)\(f(S)+int(Rk +)). Proof Assume for a contradiction that there is some ¯y∈Ssuch that f(¯y)dominates f(˜x). This contradicts ¯y∈S⊆˜ Stogether with the nondominance of f(˜x)w.r.t. (RMOP ˜ S). The above lemma tells us that any stable set ˜ Nconsisting of nondominated points of possibly different relaxations of a given (MOP) is nondominated w.r.t. N. Therefore, by employing Lemma 3.5, we obtain a lower bound set L(˜ N)of the nondominated set N.If we have furthermore any potentially nondominated—and therefore particularly stable—set N⊆f(S)available, we obtain an enclosure E(L(˜ N), U(N)) of the nondominated set in the sense of Definition 3.1 again by Lemma 3.5. That is the core idea of the method presented in the next section. We finalize this section with some aspects concerning the sets Nand ˜ N. We want to make use of the idea from [45], where the optimal value is sandwiched between lower bounds coming from relaxed problems and upper bounds coming from solutions of the reduced problems. In the multiobjective setting, the lower bound is not a single value anymore, but a stable set consisting of optimal solutions of possibly different relaxations, namely ˜ N.Inthe same spirit, also the upper bound is not a single value but a stable set consisting of potentially nondominated points, namely N.Aswewanttosomehowdecrease this set Ntowards N,it is updated w.r.t. nondominance, i.e. if we compute a new potentially nondominated point y, we add it to Nif yis nondominated w.r.t. N. Furthermore, we discard any points in Nwhich are dominated by y. This is written in Algorithm 3. 2Note that the argument showing that the refinement technique presented in this paper satisfies Assumption 5.1 heavily relies on the equidistant partitioning scheme. 123
116 Journal of Global Optimization (2023) 87:97–132 Algorithm 3 Updating the set Nw.r.t. a point y∈f(S) Require: Stable set Nand update point y∈f(S) 1: Set N=N\y∈N|y≤y 2: if yis nondominated w.r.t. Nthen 3: Set N=N∪{y} 4: end if return Updated set N By doing so, we ensure that Nstays stable and consequently improves towards Nas shown in the next lemma. Note that since N⊆f(S)we have that Nis nondominated w.r.t. N, i.e. Nactually approximates Nfrom above. Lemma 4.6 Let N1be a stable input set and let y ∈f(S). Then Algorithm 3returns a stable set N2. Furthermore, either N1⊆N2or there exists y∈N1such that y dominates y. Proof We have to distinguish two base cases. Firstly, assume there exists no y∈N1with y≤y. Then either yis not nondominated w.r.t. N1andwehavethatN1=N2,oryis nondominated w.r.t. N1andwehavethatN1⊂N2=N1∪{y}. In both cases, we have that N2is stable. Secondly, assume that there exist yi∈N1with y≤yi,i∈[s]for some s∈N.Due to the stability of N1we have that yis nondominated w.r.t. N1\{yi|i∈[s]}. Hence, N2=(N1\{yi|i∈[s]})∪{y}and N2is stable. Further, if y=y1we have that N1=N2. If otherwise, we have that ydominates y1. Let us now turn to the set ˜ Nwhich is meant to approach Nfrom below, i.e. we only want to use the best relaxed solutions available. In that setting, we call a nondominated point ˜yof a relaxation ˜ Sbetter than a nondominated point ˆyof a relaxation ˆ S,if ˆydominates ˜y.Notethat if ydominates ˜ywe know that ˜ SSholds for the corresponding relaxations. In a similar but reversed way as for N, we update the set ˜ Nas written in Algorithm 4. Algorithm 4 Updating the set ˜ Nw.r.t. a point y∈int(B)\f(S)+int Rk + Require: Stable set ˜ Nand update point y∈int(B)\f(S)+int Rk + 1: Set ˜ N=˜ N\˜y∈˜ N|˜y≤y 2: if ˜ Nis nondominated w.r.t. ythen 3: Set ˜ N=˜ N∪{y} 4: end if return Updated set ˜ N We obtain the analogue to Lemma 4.6. Note that by Lemma 4.5 we know that ˜ Nis nondominated w.r.t. N, i.e. ˜ Nactually approximates Nfrom below. Lemma 4.7 Let ˜ N1be a stable input set and let y ∈int(B)\(f(S)+int(Rk +)). Then Algorithm 4returns a stable set ˜ N2. Furthermore, either ˜ N1⊆˜ N2or there exists ˜y∈˜ N1such that ˜y dominates y. Proof We have to distinguish two base cases. Firstly, assume there exists no ˜y∈˜ N1with ˜y≤y. Then either ˜ N1is not nondominated w.r.t. yand we have that ˜ N1=˜ N2,or ˜ N1is nondominated w.r.t. yand we have that ˜ N1⊂˜ N2=˜ N1∪{y}. In both cases, we have that ˜ N2is stable. 123
Journal of Global Optimization (2023) 87:97–132 117 Secondly, assume that there exist ˜yi∈˜ N1with ˜yi≤y,i∈[s]for some s∈N.Due to the stability of ˜ N1we have that ˜ N1\{˜yi|i∈[s]} is nondominated w.r.t. y. Hence, ˜ N2=(˜ N1\{˜yi|i∈[s]})∪{y}and ˜ N2is stable. Further, if y=˜y1we have that ˜ N1=˜ N2. If otherwise, we have that ˜y1dominates y. 5 General scheme We have seen at the end of Sect.3that Algorithm 2is able to compute an enclosure as well as an approximation of the nondominated set of a given (MOP). However, if nonlinear constraint functions are present in (MOP), the scalarized problems arising during Algorithm 2 are MINLP problems. At the end of Sect.3we have also seen that if the used solver, like, e.g., SCIP, is capable of handling the occurring nonlinear constraints one could solve the scalarized single objective problems directly. Nevertheless, if the complexity of (MOP)and therefore the one of the resulting MINLP problems increases, the run time of solvers like SCIP for computing a solution to such an MINLP problem may increase, too. Note that for any computed nondominated point in our approximation, we have to solve at least one such MINLP problem, so even a small increase of computational time per problem may cause a tremendous upturn of run time of the whole procedure. In Sect.4we have presented ideas and concepts from single objective optimization of MINLP problems which are meant to reduce complexity and therefore facilitate the computations while solving a scalar MINLP problem. Now one could think of choosing a specific relaxation Rand then using one of the present algorithms for computing a representation (or enclosure) of the nondominated set of the relaxed mixed-integer linear (or convex) problem, see, e.g., [47,51] for the biobjective linear case and [24] for the multiobjective convex case, or even Algorithm 2or [23]. One could argue that if the relaxation Rsatisfies some quality criterion, e.g., a small enough estimation error εR rel, the approximation (or enclosure) of NR can be considered to be an approximation (or enclosure) of N, similar as proposed in [13] for the single objective case. Note that convergence then only relies on the theory of the used multiobjective method. However, computing such a relaxation and solving the arising problem using an available solution method may be very time-consuming as the complexity, even of the relaxed problems, may increase with ongoing refinement—in particular, as the number of integer variables increases while tightening the relaxations. One strategy for avoiding this is trying to use cheap relaxations whenever possible and refining them only when necessary, e.g., only in specific parts of the image space. This idea of adaptively refining the relaxations while computing an enclosure of the nondominated set of (MOP) is the core of this work. Note that one could additionally incorporate adaptivity in the variable domains into the refinement procedure (see Remark 4.4). In the following, we present an algorithm similar to Algorithm 2which makes use of these ideas in order to compute an enclosure of the nondominated set without solving scalarized MINLP problems, but only MILP and NLP problems (see Algorithm 5). 123
118 Journal of Global Optimization (2023) 87:97–132 Algorithm 5 General scheme for computing an enclosure of the nondominated set relying on relaxation and scalarization techniques. Require: box B=z,zuwith f(S)⊆int(B), termination tolerance εencl >0, off-set factor εencl >δ>0, initial relaxation RI, estimation error tolerance εrel >0 1: Initialize potentially nondominated set N=∅and set of local upper bounds U={zu} 2: Initialize set of best relaxed solutions ˜ N=∅and set of local lower bounds L={z} 3: Initialize the sets E=E(L,U),D(U)=zu,RIand D(˜ N)=∅ 4: while w(E)≥εencl do 5: Uloop =U 6: for u∈Uloop do 7: if there exists ∈Lwith ≤uand s(, u)≥εencl then 8: done =false 9: while done =false do 10: Set Rcurrent =min Ru,R,whereu,Ru∈D(U)and R min R˜y|˜y,R˜y∈D˜ Nand ˜y<u−εencle 11: Set relaxation ˜ S=SRcurrent 12: if there exists ˜y=f(˜x)∈N˜ Swith ˜y<u−δethen 13: if ˜ Nis nondominated w.r.t. ˜ythen 14: Update ˜ Nand Lw.r.t. ˜yusing Alg. 4and 1for local lower bounds and set R˜y=Rcurrent 15: if εRcurrent rel <ε rel then 16: Update Nand Uw.r.t. ˜yusing Alg. 3and 1and set Ru=Rcurrent for any new local upper bound u 17: done =true 18: else 19: if there exists solution yto (redMOP(˜xI)) with y<u−δethen 20: Update Nand Uw.r.t. yusing Alg. 3and 1and set Ru=Rcurrent for any new local upper bound u 21: done =true 22: else 23: Choose Rwith Rcurrent Rand set Ru=R 24: end if 25: end if 26: else 27: Choose Rwith Rcurrent Rand set Ru=R 28: end if 29: else 30: Update Lw.r.t. u−δeusing Algorithm ?? 31: done =true 32: end if 33: end while 34: end if 35: end for 36: end while return Enclosure E(L,U)satisfying w(E)<ε encl and approximation Nof Nεencl Before starting the procedure we fix a relaxation technique guaranteeing that the relaxation error quality criterion, namely εR rel <ε rel, is satisfied after a finite number of refinement steps. We introduce the set consisting of all such relaxations := R|εR rel <ε rel, and assume the following for the remainder of the paper. 123
Journal of Global Optimization (2023) 87:97–132 119 Assumption 5.1 For any relaxation technique and any initial relaxation RI.LetRIR1 R2...be a chain of strictly decreasing relaxations. Then for any εrel >0 there exists an s∈Nwith εRs rel <ε rel, i.e. Rl∈for all l≥s. Remark 5.2 We should mention that there is a wide range of alternative relaxation techniques in the literature, including the use of convex underestimators (cf. [2,3,37,50]) instead of piecewise linear ones. One could then either solve the resulting convex MINLP problems or combine them with outer approximation techniques (cf. [18,27,35,39,54]). However, here we restrict ourselves to the case of piecewise linear relaxations since the relaxation error computation and especially convergence in the sense of Assumption 5.1 is straightforward. For the above-mentioned approaches, one has to ensure that both of these requirements are fulfilled. Furthermore, if R∈we consider any feasible point ˜x∈˜ Sof the corresponding relaxed problem (RMOP ˜ S) based on the feasible set ˜ S=SRas a feasible point of (MOP), i.e. ˜x∈S. Consequently, by Lemma 4.5 for any efficient point ˜x∈˜ Sof (RMOP ˜ S)wehavethat f(˜x)∈N. This means, that for a relaxation R∈we consider any nondominated point of (RMOP ˜ S) as a nondominated point of (MOP). For ease of notation, we write ˜x∈Sif ˜x∈SRfor some R∈for the remainder of that paper. Note that Assumption 5.1 holds for the McCormick piecewise linear relaxations introduced in Sect. 4. In Step 3, we initialize the set D(U):= {(u,R)|u∈U,Rcaused the computation of u}, consisting of any present local upper bound together with the relaxation Rwhich was needed to obtain this specific local upper bound. We say that the relaxation Rcaused the computation of a local upper bound uif uentered the set of local upper bounds after it was updated w.r.t. a point ywhose computation relied on R.Thiscouldbeeitherthecaseifyis the solution of the relaxed problem corresponding to Rand R∈or if yis a nondominated point of (redMOP(˜xI)), where ˜xIwas computed using the relaxation R. Furthermore, we initialize the set D˜ N:= (˜y,R)|˜y∈˜ N,Rcaused the computation of ˜y, consisting of solutions of relaxed problems together with their corresponding relaxation R. Suppose now we are at the beginning of the l-th call of the outer while-loop in Step 4 and we have that w(E)≥εencl.WefixthesetUloop to be the current assignment of the set of local lower bounds Uand start the for-loop in Step 6. In that for-loop, let ˆu∈Uloop such that there exists ∈Lwith ≤ˆuand s(, ˆu)≥εencl, i.e. the search zone determined by and ˆuis not yet well enough explored and we set done =false. We use the indicator done to determine whether we achieved an improvement w.r.t. ˆu, i.e. the inner while-loop ensures that we concentrate on ˆuuntil we made some improvement. We say that we improved ˆuif we entered one of the if-statements in the Steps 15 and 19 or the else-statement in Step 29 as in these steps either a potentially nondominated point ywith y<ˆu−δeis found or the search region c(ˆu)is declared to be well enough explored. However, given the current local upper bound ˆuwe have to choose an appropriate relaxation Rcurrent for executing our computations. This is realized in Step 10. We choose a relaxation at least as fine as the relaxation which led to ˆu, i.e. Rˆu⊇Rcurrent. Furthermore, if there exists some relaxed solution ˜y∈˜ Nwith ˜y<ˆu−εenclewe choose a strictly finer relaxation than R˜y, i.e. we ensure that R˜yRcurrent. Note that the strictness of the inclusion is not necessary for the convergence of Algorithm 5since the method also refines the relaxations if 123
120 Journal of Global Optimization (2023) 87:97–132 Fig. 5 A strict refinement of the relaxation is needed in order to obtain solutions closer to the desired area as it is dominated by ˜y the incumbent relaxed solution did not lead to an improvement of the lower bound set. However, for some problems, it may be of advantage to refine the relaxations more aggressively instead of solving cheaper problems that do not have a high chance of leading to a significant improvement of the lower bound set. In fact, not forcing the inclusion Rcurrent ⊆R˜yto be strict makes it impossible to find a relaxed solution ydominating ˜y, and therefore satisfying ˆu−εencle≤y, as can be seen in Fig.5. The possible negative effect of not forcing strictness can be seen in the comparison of Figs.7and 8, where we can observe an increase in runtime as well as in the number of problems to be solved. However, refining too aggressively may also be a problem as it may result in solving harder problems than necessary. From our first observations, it is a good strategy to take the coarsest possible relaxation without losing the strictness of the inclusion—at least with using our basic refinement strategy. After that we initialize the relaxation ˜ S=SRcurrent ,solve(WSP (α;u)) with u=ˆufor some α∈int(Rk +)and feasible set ˜ Sand then decide whether we are able to enter the ifstatement in Step 12. If (WSP (α;u)) is infeasible, we declare the current search region c(ˆu) as well enough explored by the same arguments as in Algorithm 2,andsetdone =true. If, otherwise, there exists a solution ˜yto (WSP (α;u)) we check if ˜ Nis nondominated w.r.t. ˜y, i.e. if ˜yimproves the set ˜ N. If that is not the case, i.e. if there exists y∈˜ Nwith ˜y≤y,we have to restart the inner while-loop with a finer relaxation as the current one did not lead to any improvement. If otherwise, ˜ Nis nondominated w.r.t. ˜y, it is reasonable to move on as ˜ysuggests an improvement of ˆu. Consequently, we update the sets ˜ Nand Lw.r.t. ˜yand save the corresponding relaxation. If now the relaxation is fine enough in the sense of [13], i.e. Rcurrent ∈, we consider ˜yas a nondominated point of (MOP), update the sets Nand Uw.r.t. ˜yand set done =true. If otherwise, the relaxation is not yet fine enough we try to make use of the idea from [45], i.e. using parts of the relaxed solution to set up the reduced problem (redMOP(˜xI)). We solve the corresponding (WSP (α;u)) and if it has a solution we obtain a potentially nondominated point, i.e. update the sets Nand Uand set done =true. If it is infeasible, we have to restart with a finer relaxation, since Rcurrent suggested wrongly that we would find a potentially nondominated point. 123
Journal of Global Optimization (2023) 87:97–132 127 of the problems one might have an advantage by only considering relaxations instead of the original problem, as can be seen in the next section. 6 Application to the multiobjective optimization of decentralized energy supply networks In this section, we present numerical results of the described method on some network optimization problem. In fact, aiming to model a decentralized energy supply network we obtain a MIQCP problem. The general network structure is a graph, where the nodes represent individual consumers and the edges connect the consumer nodes with the so-called source node, where energy is supplied. The mixed-integer character is coming from certain decision options available in the optimization process, e.g., if a gas pipe is laid at some edge or not. Furthermore, we take stationary models of energy flow into account, namely an equation based on the Ohmic law for the electricity flow as well as the Darcy–Weisbach equation for gas flow. As both of them contain bilinear or quadratic terms the resulting optimization problem has the mentioned MIQCP structure. As objective functions, we use the overall costs for realizing a given network plan on the one hand and the carbon emissions of that network plan on the other hand. Naturally, a cheap network plan results in high carbon emissions, and a low carbon emission can be obtained by, e.g., investing in energy-efficient house renovation which results in higher costs. Thus, we have a classical (MOP) with two conflicting objective functions. Details regarding the modeling aspects can be found in [38] and more recently in [19]. For the present paper, we consider three network instances of such decentralized energy supply networks, namely: If e.g., we set up a single objective optimization problem with network 1 ⊂network 2 ⊂network 3 #Nodes122039 # Binaries 108 189 360 # Variables 484 829 1570 # Constraints 620 1064 2014 cost minimization as objective function and put the carbon emissions to the constraints, we obtain the following computational times •network 1:0.57s •network 2:3.35s •network 3:>3h, using the SCIP solver with the standard settings3from the pyscipopt-package (cf. [41]). For testing the new method we use a relative width tolerance εencl =0.03 as well as two off-set factors ˜ δ=0.95εencl and δ=0.8εencl. In the network models the present nonlinearities are of the following form: •For modeling the low-voltage energy flow we use for instance Re i,jfe in,i,j=aeuj¯ui,j,(13) 3Note that the computational time needed for solving network 3 is drastically reduced if one increases the tolerance for termination. 123
128 Journal of Global Optimization (2023) 87:97–132 where Re i,j>0 denotes the resistance of the underlying cable at arc (i,j),ae>0the calorific multiplier of three-phase electric power flows, fe in,i,jis a variable representing the electric power flow on arc (i,j)into j.Thevariableuidenotes the electrical voltage at node iand the variable ¯ui,j=ui−ujthe voltage drop on arc (i,j). Consequently, we have a quadratic term u2 iand a bilinear term uiujappearing in (13). For the computation of the corresponding relative relaxation errors the box constraints 360 ≤ui,uj≤440 are relevant. Thus, if we partition the corresponding intervals into requidistant intervals, i.e. use the relaxation Rr,weobtain εRr rel =802 4r2 1 max{x2|360 ≤x≤440}, and therefore the number of partitions of each interval to fall below a given tolerance εrel is given by r=40 440 1 √εrel . Thus, if we require a relative relaxation error εrel =0.01 we have to partition the corresponding intervals into at least r=1 partitions, i.e. we do not have to partition at all. •For modeling low-pressure gas supply we use a reformulation of the Darcy–Weisbach equation avoiding the use of the sign-function as proposed in [8]. By doing so, we obtain for instance Rg i,j¯q2 ij ≤¯pmax yi,j,(14) where Rg i,j>0 denotes the resistance constant of the underlying gas pipeline on arc (i,j),¯qij the gas flow on arc (i,j),¯pmax >0 the maximal pressure loss allowed in the network as well as a binary decision variable yi,jindicating if a gas pipe is laid at arc (i,j). The relevant box constraints are −150 ≤¯qij ≤150 and partitioning into r intervals, i.e. using relaxation Rr, we obtain εRr rel =3002 4r2 1 max{x2|−150 ≤x≤150}, and therefore the number of partitions of each interval to fall below a given tolerance εrel is given by r=1 √εrel . Thus, if we require a relative relaxation error εrel =0.01 we have to partition the corresponding intervals into at least r=10 partitions. Note that even if we just require a relative relaxation error εrel =0.03 we still have to partition into at least r=6intervals. In sum, this yields that—if we use requidistant partitions for any variable appearing in any nonlinear term of our problem—we fall below a relaxation error tolerance of εrel =0.01 as soon as we use a relaxation Rrwith r≥10. Note that if we did not use the adaptive approach given in Algorithm 5, but chose a relaxation with r≥10 and then used a method for solving multiobjective linear mixed-integer problems we would have to solve a problem with at least 10 ·|Edges|additional integer variables, i.e. in the case of network 3 about 400 if we just use the ones for the quadratic terms. 123
Journal of Global Optimization (2023) 87:97–132 129 Fig. 9 Computational results on network 3. Left: enclosure given by LLBs and LUBs of (Circles) obtained by Algorithm 5with WSM. Right: counter of used degrees of relaxations Looking at the results for network 3 (see Fig.9; the results for network 1 and network 2 are similar) we can see that the method uses only the relaxation Rrwith r=1, i.e. the coarsest relaxation possible using the McCormick relaxations. This shows the power of the proposed method dealing with the considered large network instances. 7 Conclusion In the present work, a general MIQCP problem is considered and two novel methods for computing an enclosure of the nondominated set are presented. For both of them, we proved correct and finite termination as well as demonstrated the respective advantages and disadvantages. The implementation of the second approach is currently only able to deal with bilinear and quadratic terms. However, one could handle general polynomial terms still relying on McCormick relaxations. For general nonlinear terms, one has to go for a more elaborate relaxation technique as presented in, e.g., [11]. However, these are only implementation issues. As long as the relaxation technique satisfies Assumption 5.1 the theoretical results presented in this paper still apply. Acknowledgements This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) via the German Excellence Strategy. We are also very thankful for the support by Piet Hensel and Dirk König from Rechenzentrum für Versorgungsnetze Wehr GmbH (http://www.rzvn.de) for their support for setting up the numerical model in Sect.6. Funding Open Access funding enabled and organized by Projekt DEAL. Data availability The datasets used for computations in Sect.6are available from the corresponding author on reasonable request. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 123
130 Journal of Global Optimization (2023) 87:97–132 References 1. Achterberg, T., Wunderling, R.: Mixed integer programming: analyzing 12 years of progress. In: Facets of Combinatorial Optimization, pp. 449–481. Springer, Heidelberg (2013). https://doi.org/10.1007/9783-642-38189-8_18 2. Adjiman, C.S., Dallwig, S., Floudas, C.A., Neumaier, A.: A global optimization method, αBB, for general twice-differentiable constrained NLPs—I. Theoretical advances. Comput. Chem. Eng. 22(9), 1137–1158 (1998). https://doi.org/10.1016/S0098-1354(98)00027-1 3. Androulakis, I.P., Maranas, C.D., Floudas, C.A.: αBB: a global optimization method for general constrained nonconvex problems, vol. 7, pp. 337–363 (1995). https://doi.org/10.1007/BF01099647. State of the art in global optimization: computational methods and applications (Princeton, NJ, 1995). https://doi. org/10.1007/BF01099647 4. Banholzer, S.: Rom-Based Multiobjective Optimization with PDE Constraints. Ph.D. thesis, Universität Konstanz, Konstanz (2021) 5. Banholzer, S., Gebken, B., Dellnitz, M., Peitz, S., Volkwein, S.: ROM-based multiobjective optimization of elliptic PDEs via numerical continuation. In: Non-smooth and Complementarity-based Distributed Parameter Systems—Simulation and Hierarchical Optimization. Int. Ser. Numer. Math., vol. 172, pp. 43–76. Birkhäuser/Springer, Cham (2022). https://doi.org/10.1007/978-3-030-79393-7_3 6. Belotti, P., Kirches, C., Leyffer, S., Linderoth, J., Luedtke, J., Mahajan, A.: Mixed-integer nonlinear optimization. Acta Numer. 22, 1–131 (2013). https://doi.org/10.1017/S0962492913000032 7. Bestuzheva, K., Besançon, M., Chen, W.-K., Chmiela, A., Donkiewicz, T., van Doornmalen, J., Eifler, L., Gaul, O., Gamrath, G., Gleixner, A., Gottwald, L., Graczyk, C., Halbig, K., Hoen, A., Hojny, C., van der Hulst, R., Koch, T., Lübbecke, M., Maher, S.J., Matter, F., Mühmer, E., Müller, B., Pfetsch, M.E., Rehfeldt, D., Schlein, S., Schlösser, F., Serrano, F., Shinano, Y., Sofranac, B., Turner, M., Vigerske, S., Wegscheider, F., Wellner, P., Weninger, D., Witzig, J.: The SCIP Optimization Suite 8.0. ZIB-Report 21-41, Zuse Institute Berlin (2021). http://nbn-resolving.de/urn:nbn:de:0297-zib-85309 8. Borraz-Sánchez, C., Bent, R., Backhaus, S., Hijazi, H., Van Hentenryck, P.: Convex relaxations for gas expansion planning. INFORMS J. Comput. 28(4), 645–656 (2016). https://doi.org/10.1287/ijoc.2016. 0697 9. Boukouvala, F., Misener, R., Floudas, C.A.: Global optimization advances in mixed-integer nonlinear programming, MINLP, and constrained derivative-free optimization, CDFO. Eur. J. Oper. Res. 252(3), 701–727 (2016). https://doi.org/10.1016/j.ejor.2015.12.018 10. Burachik, R.S., Kaya, C.Y., Rizvi, M.M.: Algorithms for generating pareto fronts of multi-objective integer and mixed-integer programming problems. Eng. Optim. (2021). https://doi.org/10.1080/0305215X.2021. 1939695 11. Burlacu, R.: Adaptive mixed-integer refinements for solving nonlinear problems with discrete decisions. Ph.D. thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg (2019) 12. Burlacu, R.: On refinement strategies for solving MINLPs by piecewise linear relaxations: a general red refinement. Optim. Lett. (2021). https://doi.org/10.1007/s11590-021-01740-1 13. Burlacu, R., Geißler, B., Schewe, L.: Solving mixed-integer nonlinear programmes using adaptively refined mixed-integer linear programmes. Optim. Methods Softw. 35(1), 37–64 (2020). https://doi.org/ 10.1080/10556788.2018.1556661 14. Cplex, I.I.: V12. 1: user’s manual for CPLEX. Int. Bus. Mach. Corp. 46(53), 157 (2009) 15. Dächert, K., Klamroth, K., Lacour, R., Vanderpooten, D.: Efficient computation of the search region in multi-objective optimization. Eur. J. Oper. Res. 260(3), 841–855 (2017). https://doi.org/10.1016/j.ejor. 2016.05.029 16. De Santis, M., Eichfelder, G., Niebling, J., Rocktäschel, S.: Solving multiobjective mixed integer convex optimization problems. SIAM J. Optim. 30(4), 3122–3145 (2020). https://doi.org/10.1137/19M1264709 17. Diessel, E.: An adaptive patch approximation algorithm for bicriteria convex mixed-integer problems. Optimization 0(0), 1–46 (2021). https://doi.org/10.1080/02331934.2021.1939699 18. Duran, M.A., Grossmann, I.E.: An outer-approximation algorithm for a class of mixed-integer nonlinear programs. Math. Program. 36(3), 307–339 (1986). https://doi.org/10.1007/BF02592064 19. Eggen, C., Huynh, T.-V., Link, M., Stephan, P., Volkwein, S.: An MINLP model for designing decentralized energy supply networks. Technical report. arXiv:2212.06527 (2022) 20. Ehrgott, M.: Multicriteria Optimization, 2nd edn., p. 323. Springer, Berlin (2005) 21. Ehrgott, M., Gandibleux, X.: Bound sets for biobjective combinatorial optimization problems. Comput. Oper. Res. 34(9), 2674–2694 (2007). https://doi.org/10.1016/j.cor.2005.10.003 22. Eichfelder, G.: Twenty years of continuous multiobjective optimization in the twenty-first century. EURO J. Comput. Optim. 9, 100014 (2021). https://doi.org/10.1016/j.ejco.2021.100014 123
Journal of Global Optimization (2023) 87:97–132 131 23. Eichfelder, G., Warnow, L.: An approximation algorithm for multi-objective optimization problems using a box-coverage. J. Global Optim. (2021) 24. Eichfelder, G., Warnow, L.: A hybrid patch decomposition approach to compute an enclosure for multiobjective mixed-integer convex optimization problems (2021) 25. Eichfelder, G., Stein, O., Warnow, L.: A deterministic solver for multiobjective mixed-integer convex and nonconvex optimization (2022) 26. Eichfelder, G., Kirst, P., Meng, L., Stein, O.: A general branch-and-bound framework for continuous global multiobjective optimization. J. Global Optim. 80(1), 195–227 (2021). https://doi.org/10.1007/ s10898-020-00984-y 27. Fletcher, R., Leyffer, S.: Solving mixed integer nonlinear programs by outer approximation. Math. Program. 66(3, Ser. A), 327–349 (1994). https://doi.org/10.1007/BF01581153 28. Geißler, B., Martin, A., Morsi, A., Schewe, L.: Using piecewise linear functions for solving MINLPs. In: Mixed Integer Nonlinear Programming. IMA Vol. Math. Appl., vol. 154, pp. 287–314. Springer, New York (2012) 29. Gurobi Optimization, LLC: Gurobi Optimizer Reference Manual (2022). https://www.gurobi.com 30. Iapichino, L., Trenz, S., Volkwein, S.: Reduced-order multiobjective optimal control of semilinear parabolic problems. In: Numerical Mathematics and Advanced Applications—ENUMATH 2015. Lect. Notes Comput. Sci. Eng., vol. 112, pp. 389–397. Springer, Cham (2016) 31. Kim, I.Y., de Weck, O.: Adaptive weighted sum method for bi-objective optimization: Pareto front generation. Struct. Multidiscip. Optim. 29, 149–158 (2005). https://doi.org/10.1007/s00158-004-0465-1 32. Kim, I.Y., de Weck, O.: Adaptive weighted sum method for multiobjective optimization: a new method for Pareto front generation. Struct. Multidiscip. Optim. 31(2), 105–116 (2006). https://doi.org/10.1007/ s00158-005-0557-6 33. Kirlik, G., Sayın, S.: Bilevel programming for generating discrete representations in multiobjective optimization. Math. Program. 169(2), 585–604 (2018). https://doi.org/10.1007/s10107-017-1149-0 34. Klamroth, K., Lacour, R., Vanderpooten, D.: On the representation of the search region in multi-objective optimization. Eur. J. Oper. Res. 245(3), 767–778 (2015). https://doi.org/10.1016/j.ejor.2015.03.031 35. Kronqvist, J., Lundell, A., Westerlund, T.: The extended supporting hyperplane algorithm for convex mixed-integer nonlinear programming. J. Global Optim. 64(2), 249–272 (2016). https://doi.org/10.1007/ s10898-015-0322-3 36. Lee, J., Leyffer, S. (eds.): Mixed Integer Nonlinear Programming. The IMA Volumes in Mathematics and its Applications, vol. 154, p. 690. Springer, New York (2012). https://doi.org/10.1007/978-1-4614-19273. Selected papers based on the IMA Hot Topics Workshop “Mixed-Integer Nonlinear Optimization: Algorithmic Advances and Applications” held in Minneapolis, MN, November 17–21, 2008. https://doi. org/10.1007/978-1-4614-1927-3 37. Liberti, L.: Reformulation and convex relaxation techniques for global optimization. Q. J. Belg. Fr. Ital. Oper. Res. Soc. 2, 255–258 (2004). https://doi.org/10.1007/s10288-004-0038-6 38. Lu, J.: Mixed-Integer Nonlinear Modeling and Optimization of Designing Decentralized Energy Supply Networks. Ph.D. thesis, Universität Konstanz, Konstanz (2023) 39. Lundell, A., Kronqvist, J.: Polyhedral approximation strategies for nonconvex mixed-integer nonlinear programming in SHOT. J. Global Optim. 82(4), 863–896 (2022). https://doi.org/10.1007/s10898-02101006-1 40. Lundell, A., Skjäl, A., Westerlund, T.: A reformulation framework for global optimization. J. Global Optim. 57(1), 115–141 (2013). https://doi.org/10.1007/s10898-012-9877-4 41. Maher, S., Miltenberger, M., Pedroso, J.P., Rehfeldt, D., Schwarz, R., Serrano, F.: PySCIPOpt: mathematical programming in python with the SCIP optimization suite. In: Mathematical Software—ICMS 2016, pp. 301–307. Springer, Cham (2016). https://doi.org/10.1007/978-3-319-42432-3_37 42. McCormick, G.P.: Computability of global solutions to factorable nonconvex programs. I. Convex underestimating problems. Math. Program. 10(2), 147–175 (1976). https://doi.org/10.1007/BF01580665 43. Misener, R., Floudas, C.A.: Global optimization of mixed-integer quadratically-constrained quadratic programs (MIQCQP) through piecewise-linear and edge-concave relaxations. Math. Program. 136(1, Ser. B), 155–182 (2012). https://doi.org/10.1007/s10107-012-0555-6 44. Morsi, A., Geißler, B., Martin, A.: Mixed integer optimization of water supply networks. In: Mathematical Optimization of Water Networks. Internat. Ser. Numer. Math., vol. 162, pp. 35–54. Birkhäuser/Springer Basel AG, Basel (2012). https://doi.org/10.1007/978-3-0348-0436-3_3 45. Nagarajan, H., Lu, M., Wang, S., Bent, R., Sundar, K.: An adaptive, multivariate partitioning algorithm for global optimization of nonconvex programs. J. Global Optim. 74(4), 639–675 (2019). https://doi.org/ 10.1007/s10898-018-00734-1 46. Pascoletti, A., Serafini, P.: Scalarizing vector optimization problems. J. Optim. Theory Appl. 42(4), 499– 524 (1984). https://doi.org/10.1007/BF00934564 123
132 Journal of Global Optimization (2023) 87:97–132 47. Perini, T., Boland, N., Pecin, D., Savelsbergh, M.: A criterion space method for biobjective mixed integer programming: the boxed line method. INFORMS J. Comput. 32(1), 16–39 (2020). https://doi.org/10. 1287/ijoc.2019.0887 48. Ryu, N., Min, S.: Multiobjective optimization with an adaptive weight determination scheme using the concept of hyperplane. Int. J. Numer. Methods Eng. 118(6), 303–319 (2019). https://doi.org/10.1002/ nme.6013 49. Sayın, S.: Measuring the quality of discrete representations of efficient sets in multiple objective mathematical programming. Math. Program. 87(3, Ser. A), 543–560 (2000). https://doi.org/10.1007/ s101070050128 50. Skjäl, A.: On the use of convex under estimators in global optimization. Ph.D. thesis, Abo Akademi University (2014) 51. Stidsen, T., Andersen, K.A.: A hybrid approach for biobjective optimization. Discrete Optim. 28, 89–114 (2018). https://doi.org/10.1016/j.disopt.2018.02.001 52. Tawarmalani, M., Sahinidis, N.V.: Convexification and Global Optimization in Continuous and Mixedinteger Nonlinear Programming. Nonconvex Optimization and its Applications, vol. 65, p. 475. Kluwer Academic Publishers, Dordrecht (2002). https://doi.org/10.1007/978-1-4757-3532-1. Theory, algorithms, software, and applications 53. Vielma, J.P., Ahmed, S., Nemhauser, G.: Mixed-integer models for nonseparable piecewise-linear optimization: unifying framework and extensions. Oper. Res. 58(2), 303–315 (2010). https://doi.org/10.1287/ opre.1090.0721 54. Westerlund, T., Pettersson, F.: An extended cutting plane method for solving convex MINLP problems. Comput. Chem. Eng. 19, 131–136 (1995). https://doi.org/10.1016/0098-1354(95)87027-X. European Symposium on Computer Aided Process Engineering 3–5 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123