Tight Convex Relaxations for the Expansion Planning Problem

Secure energy transport is considered as highly relevant for the basic infrastructure of nowadays society and economy. To satisfy increasing demands and to handle more diverse transport situations, operators of energy networks regularly expand the capacity of their network by building new network elements, known as the expansion planning problem. A key constraint function in expansion planning problems is a nonlinear and nonconvex potential loss function. In order to improve the algorithmic performance of state-of-the-art MINLP solvers, this paper presents an algebraic description for the convex envelope of this function. Through a thorough computational study, we show that this tighter relaxation tremendously improves the performance of the MINLP solver SCIP on a large test set of practically relevant instances for the expansion planning problem. In particular, the results show that our achievements lead to an improvement of the solver performance for a development version by up to 58%.


Introduction
Expansion planning problems of energy networks (see Sect. 2) are becoming more relevant as well as more challenging due to increasing demands of modern society. The goal of this paper is to improve the algorithmic performance of an LP-based branch-and-bound MINLP solver for the expansion planning problem. In particular, we focus on building a tighter LP relaxation by exploiting the convex envelope of the main nonlinearity present in this problem. To this end, we compute the convex envelope and implement it in the MINLP solver SCIP.
The state of the art for solving nonconvex MINLPs to global optimality is spatial branch-and-bound. Within spatial branch-and-bound one of the most important ingredients is the construction of tight convex relaxations. For a nonconvex constraint {x ∈ D | f (x) ≤ 0}, where D ⊆ R n is convex, the most common method for constructing a convex relaxation is to compute a convex underestimator, that is, a function g : D → R that satisfies g(x) ≤ f (x) for all x ∈ D. Then, g(x) ≤ 0 yields a convex relaxation of the constraint f (x) ≤ 0.
The best possible convex underestimator of a nonconvex function f over the domain D is given by the convex envelope vex D ( f ). For general functions, computing the convex envelope is a difficult task but for specific functions it might be tractable. Note that constraints can also be of the form f (x) ≥ 0, in which case concave overestimators are required to build a convex relaxation. In particular, for equality constraints, both the convex under and concave overestimators are needed.
In this paper, we analytically derive the convex and concave envelope of the potential loss function f (x, y) = y sgn(x) |x| α with α > 1. This is the only source of nonlinearity and nonconvexity of the expansion planning problems that we consider, see Sect. 2. Thus, the convex relaxations depend on the convexification of these constraints. As we will see in Sect. 4.1, the concave envelope of f can be deduced from the convex envelope of f , and thus, we concentrate on the convex envelope in the following. To the best of our knowledge, the convex envelope of this function is unknown in the literature. Contributions (i) We derive an explicit description for the convex and concave envelope of f over a rectangular domain D = [x, x] × [y, y] with x < 0 < x and 0 < y < y. (ii) We perform an exhaustive computational study, we show that these tighter relaxations tremendously improve the performance of the MINLP solver SCIP on a large test set of practically relevant instances that stem from real-world applications. Our results show that the performance improves up to 58% for hard instances.
Outline The rest of the paper is structured as follows: Sect. 2 introduces the expansion planning problem. In Sect. 3, we introduce preliminary theoretical results from convex analysis that we need to derive the convex envelope of y sgn(x) |x| α . In Sect. 4, we present our main contribution, which is the derivation of a closed-form expression for the convex envelope of f . In Sect. 5, we compare the convex envelope with the standard factorable relaxation for the constraint function f followed by an extensive computational study in Sect. 6. We test the impact of the convex envelope over real-world expansion planning problem instances using the MINLP solver SCIP. To this end, we first study the impact on the dual bound in the root node and secondly the overall performance within the spatial branch-and-bound search. Finally, Sect. 7 provides conclusions.

Expansion Planning Problem
Transmission system operators of energy networks such as natural gas, hydrogen and water networks regularly face both increasing demand and more diverse transport situations. To deal with these challenges, they must expand the capacity of the existing networks without having to resort to the expensive options of setting up new pipeline corridors or replacing pipelines by larger ones. In practice, operators preferably build new pipelines parallel to existing ones. This process is referred to as "looping" in the terminology of the industry. It is significantly cheaper than, for example, building completely new pipelines, since existing rights of way can be used simplifying the regulatory and the planning processes and, moreover, the additional encroachment on the environment is rather moderate. In the literature, two different strands have been considered for modeling the looping problem and related expansion planning problems: the split-pipe and the discrete approach. In the first approach, the entire length of the pipeline is split into several segments each of which may have its own length and diameter (see, e.g., [9,18,41]), and the discrete approach, where the optimal diameters are chosen from a set of commercially available diameters (see, e.g., [4,7,17]).
In this paper, we consider expansion planning based on the split-pipe approach. More precisely, we use the formulation of Model E for the looping problem from [21], since the authors show that it performed best among all studied model formulations. In the following, we briefly explain the split-pipe approach and state its model formulation for the loop expansion planning problem.
Split-pipe approach The starting point is the concept of the so-called equivalent diameter. It allows to replace multiple pipelines in parallel by a single pipeline without changing any physical properties of the network, see, for example, [4,7,19] for two parallel pipelines in gas and water networks. For a generalization of the equivalent diameter to more than two pipelines in parallel, we refer to [20]. A type of equivalent diameter also exists for the case of multiple serial pipelines with no intermediate source or sink. When modeling split-pipes, it suffices to consider only those equivalent diameters resulting from parallel pipelines (multiple loops) and serial pipelines (splitpipe setting) that are not dominated by others with lower costs. These nondominated points together with their nondominated convex combinations form the so-called efficient frontier. This means that all points on the efficient frontier represent the cost minimal way of equipping a pipeline with a certain equivalent diameter. For more detailed information, we refer to Section 2 in [21].
Model formulation for the split-pipe approach Let G = (V, A) be a directed and connected graph with node set V and pipe set A. A balanced demand vector b ∈ R |V| specifies nodal inflow and outflow values of the network. Here, we deal with a stationary model, this means that the overall network demand satisfies v∈V b v = 0. The physical state of the network is described by flow variables x a ∈ [x a , x a ] for all arcs a = (v, w) ∈ A and nonnegative potential variables π v ∈ [π v , π v ] at each node v ∈ V. The arc flows arise from a friction-induced potential difference between the arcs' adjacent nodes. In stationary models, it is a common approach to represent this interdependency by virtue of the Weymouth equation, see [38], with α = 2 for gas network operations and by virtue of Darcy-Weisbach with α = 2 or Hazen-Williams with α = 1.852 for water network operations, see [37], in any case taking the shape of Eq. (1a). Note that sgn(x a ) indicates the direction of the arc flow. Given that, in general, bi-directional arc flow is possible, we assume that x < 0 and x > 0. Here, for each pipe a ∈ A, positive variables y a ∈ [y a , y a ] and c a ≥ 0 model the efficient frontier of the split-pipe approach described above, where y a represents the physical impact of the expansion on the potential loss along arc a with corresponding costs c a that can be assumed as proportional to the pipe length L a . Lenz [20] shows that the best expansion candidates in the split-pipe setting are determined by the efficient frontier [Eq. (1b)], which is given by k a -many nondominated points in the (y a , c a )-system. In the case that pipe a is not looped, then by construction it follows y a = y a with associated costs c a = 0. Otherwise, y a > y and c a > 0. Finally, as in classical network flow problems, Eq. (1c) models the flow conservation at every network node. Then, the model reads, In this paper, we compute the convex hull of the constraint (1a) by studying the convex and concave envelope of the nonlinear function (x, y) → y sgn(x) |x| α and then investigating its impact on the performance of expansion planning problem instances. Finally, let us mention that the nonlinear nature of the expansion planning problem stimulated a number of different model formulations and solution approaches. For example, [3] introduced the linear programming gradient (LPG) method, a two-stage heuristic that alternates between solving a linear program with fixed flows in order to obtain pipeline diameters and modifying the flow values on the basis of a sensitivity analysis. André et al. [4] tackle the discrete version of the problem heuristically in two stages, where a first step solves the continuous relaxation in order to identify pipes to be looped and a second step determines discrete-valued diameters for these selected pipes based on branch-and-bound. Bragalli et al. [7] utilize a continuous reformulation of the cost function to solve the model directly with an MINLP solver. Raghunathan [27] proposes a disjunctive program together with a convex relaxation and then solves it to global optimality. Humpola [15] formulates a model and solves it to global optimality by using convex reformulation techniques and special tailored cuts. Borraz-Sánchez et al. [6] explicitly model forward and backward flow directions by using binary variables and present a solution approach based on a second-order cone relaxation.

Preliminaries on Convex Envelopes
In general, the derivation of convex and concave envelopes is a challenging task depending on the nonlinear constraint function and on the domain of the variables involved, see, for example, [14]. It can be shown that for a lower semicontinuous function f over a compact and convex domain D, min x∈D f (x) = min x∈D vex D [ f ](x) and hence computing the convex envelope is in general N P-hard. Unless P = N P, no efficient method exists to compute convex envelopes of general nonlinear functions.
In order to build tight convex relaxations, convex envelopes have been widely studied for particular functions or classes of functions. For example, [25] constructs a linear relaxation over a box domain for the bilinear term x y, which is proven to be the convex envelope by [2]. In the last decade, several approaches arose to convexify the bilinear term x y over special domains, see, e.g., [13,23,24,26,30]. The convex envelope of the fractional term y/x has been studied, e.g., in [33,40].
Moreover, [33] computes the convex envelope of the function f (x, y) = g(y) x 2 over a rectangular domain D ⊆ R 2 , where f : D → R is convex in x and concave in y. Even though this function is close to our constraint function f (x, y) = y sgn(x) |x| α , this approach is not applicable in our setting, given that our function involves the nonconvex term sgn(x) |x| α instead of x 2 .
Note that f (x, y) = y sgn(x) |x| α is twice differentiable only for α > 2 and thus it is not possible to compute the convex envelope using techniques based on the Hessian, as done, for example, for (n − 1)-convex functions in [16]. In principle, when f ∈ C 2 (i.e., α > 2), we could generate valid tangential hyperplanes of function f to separate violated LP solutions without explicitly calculating the convex envelope, as done in [5] for twice continuously differentiable functions with fixed convexity behavior over a box. However, our analysis is valid for any α > 1.
In the remainder of this section, we introduce some mathematical notation and necessary theoretical results that are required to derive vex D ( f ). Even though our function is continuous, many of the following results hold for lower semicontinuous functions, which we define as follows.
We continue by formally defining the convex envelope of a function.

Definition 3.3
The convex envelope of a function f : D → R over a convex set D ⊆ R n is given by the tightest convex underestimator η : The concave envelope of a function f : D → R over D ⊆ R n is given by the tightest concave overestimator η : D → R of f , defined pointwise as A geometrical characterization of the convex envelope is as follows.
Theorem 3.1 [28,Theorem 5.3] Let E be any convex set in R n+1 , and let Then, g : R n → R is a convex function.
The convex envelope of a l.s.c. function f corresponds to the function g, obtained by applying Theorem 3.1 to the convex set conv(epi D f ). Therefore, it follows see also [28]. In order to determine the convex envelope of a l.s.c. function f at point x ∈ D, its representation in (2) leads to the following minimization problem Note that the n + 1 is justified by Caratheodory's theorem [28,Theorem 17.1]. Problem (3) is nonconvex and nonlinear as the multipliers λ k and the variables x k are unknown. However, it is possible to reduce the degree of nonlinearity and nonconvexity of Problem (3) by studying the convex set conv(epi D f ) in more detail.
In particular, we reduce the domain D by determining a subset A ⊆ D, such that conv(epi D f ) = conv(epi A f ) holds.

Definition 3.4
Let D be compact and convex subset of R n , and let f : D → R be a l.s.c. function.
The generating set of the convex envelope of f is defined by Given that a closed convex set containing no lines is characterized by the set of its extreme points and extreme directions [28,Theorem 18.5], we have the following result, see also [34].

Theorem 3.2 Let D be compact and convex subset of
Even though there exists no general formula to compute G D ( f ) of a nonconvex function f , [34] provides a condition that determines when a point For this, we need the relative interior of a convex set D, ri(D)

Computing the Convex Envelope of f (x, y) = y sgn(x) |x|Į
n this section, we show how to compute the convex and concave envelope of f : In the remainder of this paper, we assume y > 0, x < 0 and x > 0, as it typically occurs in practical applications, see Sect. 2. For a visualization of the nonconvex function f with α = 2, we refer to Fig. 5.

Preliminaries about the Envelopes of f(x, y) = y sgn(x) |x|P
the concave envelope of f y Note that the above problem has nine variables (x k , y k , λ k for 1 ≤ k ≤ 3). In the next section, we show how to reduce this problem to a one-dimensional problem.
Here, the structure of the function f (x, y) = y sgn(x) |x| α allows us to deduce both envelopes from each other, see Fig. 1. The function f (x, y) = y sgn(x) |x| α is odd in x when y is fixed. This allows us to retrieve the concave envelope from the convex envelope, and thus, from now on, we restrict ourselves to the derivation of the convex envelope. Thus,

Problem Reduction by Using the Generating Set
To simplify Problem (3), we use Corollary 3.1 to find a set First, notice that after fixing x to a value in [x, x], the function f (x, y) is linear and thus concave over the line segment l = {(x, y) | y ∈ [y, y]}. Then, Corollary 3.1 implies that all points in ri(l) The convex hulls conv(epi D y f ) and conv(epi D y f ) correspond to the epigraphs of the convex envelopes of the one-dimensional functions f y (x) := f (x, y) and f y (x) := f (x, y). The following proposition introduces these one-dimensional convex envelopes, see Fig. 2 for an illustration.
where γ = min{x, βx} and β is the unique negative root of (α − 1) is not the convex envelope of f (x, y) as can be seen from Fig. 6 (middle). This function is not even convex! Given that conv(epi D y f ) = epi D y ϕ y and conv(epi D y f ) = epi D y ϕ y , we can further rewrite (4) to x Fig. 2 The convex envelope ϕ y (red) of the one-dimensional function f y (blue) depends on the relation of the values βx and x. In the case of βx < x, the convex envelope is given by the secant between the points Hence, we just need to determine the convex envelope of f for y ∈ (y, y). In the following, we assume that y ∈ (y, y).

Reduction to a One-Dimensional Optimization Problem
We can further reduce the optimization Problem (6) to a one-dimensional problem.
Firstly, (6b) enables us to fix the multiplier λ Given that λ is fixed, we can use (6a) to define x 2 in terms of x 1 . For this, we define the functions t x,y , T x,y : [x, x] → R given by These functions satisfy x 2 = t x,y (x 1 ) and x 1 = T x,y (x 2 ) for every feasible solution x 1 , x 2 of Problem (6). Both functions are affine linear, strictly decreasing and inverse to each other. For an interpretation of T x,y , see Fig. 3. Note that t x,y and T x,y are well defined, since we assume λ y ∈ (0, 1), as described above.
Projecting out x 2 in Problem (6) yields a one-dimensional problem, which we can express with the help of t x,y and T x,y as The expression for the bounds come from combining both bounds of x 1 , namely The one-dimensional functions in blue f y (x) = y sgn(x) |x| α and f y (x) = y sgn(x) |x| α are obtained after fixing y to y and y in f (x, y) = y sgn(x) |x| α . Their corresponding convex envelopes ϕ y (x) and ϕ y (x) given by (5) are shown in red, and their epigraphs epi D y ϕ y and epi D y ϕ y in gray. Red crosses in the (x, y)-space correspond to the points (x 1 , y) and (t x,y (x 1 ), y) with x 2 = t x,y (x 1 ), and red crosses in the (x, y, ϕ y (x 1 ))-space correspond to the points (x 1 , y, ϕ y (x 1 )) and t x,y (x 1 ), y, ϕ y (t x,y (x 1 )) Then, Problem (9) is equivalent to (11) allows for a geometrical interpretation as illustrated in Fig. 4: For a given point (x, y) ∈ D, Problem (11) corresponds to finding the points (x 1 , y) and (t x,y (x 1 ), y) such that the line segment between these points contains (x, y), and the line segment with endpoints (x 1 , y, ϕ y (x 1 )) and t x,y (x 1 ), y, ϕ y (t x,y (x 1 )) has the lowest height at (x, y). Note that the compact formulation of Problem (11) only contains x 1 as a variable, where constraint (7) allows us to recover x 2 = t x,y (x 1 ).

Remark 4.2 Problem
At this point, we can already compute the convex envelope for a simple case.
Proof When x < βx holds, then ϕ y (x) = αx α−1 x y + (1 − α)x α y. Therefore, F x,y is a linear function and (11) reduces to minimizing a linear function over an interval. A simple computation shows that the coefficient of x 1 in F x,y (x 1 ) is negative, and so the minimum is achieved at the upper bound, i.e., vex D [ f ](x, y) = F x,y (min x, T x,y (x) ).
In light of the previous proposition, we assume in the remainder of this section that x ≥ βx. In particular, γ = βx (see Proposition 4.2).

Properties of the Objective Function
The objective function of Problem (11), F x,y , is convex as it is a convex combination of convex functions. Note that the function ϕ y (t x,y (·)) is convex since it is the composition of a convex (ϕ y ) with an affine linear function (t x,y ).
The following proposition characterizes the optimal solution for a large class of problems of a form similar to Problem (11).

Proposition 4.4
Suppose that a function F : R → R is convex and has a global minimum z * ∈ R. Then, the solution of min{F(z) : a ≤ z ≤ b} is given by where mid(a, b, z * ) selects the middle value between a, b, z * .
Proof The claim follows directly by comparing F(a), F(b) and F(z * ). If z * is within the bounds, i.e., a ≤ z * ≤ b, then F(mid(a, b, z * )) = F(z * ). If z * ≤ a, then the convexity of F implies that F is increasing in [z * , ∞) and so F(a) = F(mid(a, b, z * )) is the minimum. The argument is analogous for the case z * ≥ b.
In order to apply Proposition 4.4, we need to extend F x,y from x 1 ∈ [x 1 , x 1 ] to x 1 ∈ R and show that it has a global minimum. The extension is immediate, given that all the functions involved can be evaluated in R. To show that F x,y has a global minimum over R, we show that the sublevel sets of F x,y are bounded. This follows from the following proposition.
Proof We compute the limit when x 1 goes to +∞. The case when x 1 goes to −∞ is similar. By definition, For every large enough x 1 , we have x 1 ≥ βx and t x,y (x 1 ) < βx. Thus, This can be rewritten as ax α 1 + bx 1 + c, where a > 0, b < 0 (since t x,y is strictly decreasing and linear), and c ∈ R. Given that α > 1, the sign of a determines the behavior of F x,y at infinity; hence, lim x 1 →∞ F x,y (x 1 ) = ∞.

Solution to Problem (11)
With these properties of F x,y , we are ready to compute its global minimum over R.
To this end, we use that F x,y is differentiable, since ϕ y ∈ C 1 for α > 1, and solve F x,y (x 1 ) = 0. We have Then, given the piecewise nature of ϕ y , see (5), we get where To solve F x,y (x 1 ) = 0, we restrict x 1 to be in L L, L R, RL or R R. Given that y < y, it follows F x,y (x 1 ) = 0 for all x 1 ∈ L L. Thus, the global minimum is not in L L. Case x 1 ∈ L R: The solution of F x,y (x 1 ) = 0 is given by Note, however, that x lr / ∈ L R, as t x,y (x lr ) < βx. Thus, the global minimum is not in L R. Case x 1 ∈ RL: The solution of F x,y (x 1 ) = 0 is given by The point x rl is in RL if and only if x rl ≥ βx and t x,y (x rl ) ≤ βx. Given that β, x < 0, we have that x rl > βx. Thus, the global minimum is x rl if and only if t x,y (x rl ) ≤ βx. Otherwise, the global minimum must be in the next case. Case x 1 ∈ R R: The solution of F x,y (x 1 ) = 0 is given by From the above case distinction, it follows that the optimal solution is x rl if t x,y (x rl ) ≤ βx, otherwise it is x rr . Therefore, from Proposition 4.4, we conclude that Figure 6 shows a visualization of y sgn(x) |x| α and its convex envelope.

Simplifying the Convex Envelope
We now derive a representation of vex D [ f ](x, y) more compact than (15).

Theorem 4.1 The convex envelope of f (x, y) = y sgn(x) |x| α over D = [x, x]×[y, y] is given by
Proof To prove the claim, we rewrite (15). We first show that where x rl and x rr are given in (13) and (14). Then, we can rewrite (15) as Then, we show that This implies that mid(x 1 , x 1 , max{x rl , x rr }) = min{x 1 , max{x rl , x rr }} = min x, T x,y (x), max{x rl , x rr } , which proves the result.
To prove (17) and (18), define G x,y : R → R as and notice that G x,y (x rr ) = 0 and G x,y is strictly increasing. The equivalence (17) follows from the following computation, To prove (18), it is enough to show x 1 ≤ x rr . Recall that x 1 = max x, T x,y (x) . The inequality x ≤ x rr follows from the definition of x rr and λ y , y, y > 0.
Finally, let us show T x,y (x) ≤ x rr . This follows from the following computation, The condition x rl ≥ x depends only on the bounds, and not on the given point (x, y). Therefore, if the bounds are such that x rl ≥ x, then we can further reduce the description of the convex envelope as above.
In summary, we derived an analytic solution for the convex envelope of the function f (x, y) = y sgn(x) |x| α , with α > 1, which is stated in (16). To evaluate the convex envelope at a point (x, y) only requires the computation of the function F x,y at the minimum between x 1 , x rl and x rr .

Convexification of y sgn(x) |x|˛in SCIP
In this section, we compare the convex underestimator induced by SCIP to the convex envelope of y sgn(x) |x| α . To this end, consider the epigraph of y sgn(x) |x| α , By virtue of Theorem 3.1, any convex relaxation of (19) induces a convex underestimator of y sgn(x) |x| α . Therefore, we construct the convex underestimator of y sgn(x) |x| α induced by SCIP's convex relaxation of (19) and compare it to the convex envelope of y sgn(x) |x| α . SCIP, as other MINLP solvers, uses a standard reformulation procedure that allows to generate a convex relaxation for factorable MINLPs. The idea is to decompose a factorable constraint function into simpler functions for which convex underestimators or even convex envelopes are known. For general information about this reformulation method, we refer to [32], and for a description of how it works in SCIP to [35].
The decomposition is achieved by introducing additional variables forming, thus, an extended formulation of the original problem. In our case, SCIP decomposes the constraint θ ≥ y sgn(x) |x| α by introducing an auxiliary variable z into Then, SCIP builds a relaxation by computing the convex envelopes of f 1 (x) = sgn(x) |x| α and f 2 (y, z) = y z. Note that by Proposition 4.2, the convex envelope of f 1 is ϕ y with y = 1, which we denote by ϕ 1 . Moreover, the convex envelope of f 2 is given by the well-known McCormick inequalities [25], The convex relaxation of (19) constructed by SCIP is Hence, the underestimator induced by C is Given that ψ is increasing with respect to z for every y, an optimal solution of Problem (20) satisfies z = ϕ 1 (x). Therefore, the convex underestimator induced by SCIP is φ S (x, y) = ψ(y, ϕ 1 (x)).  (min x, T x,y (x) ). As βx ≥ x, we have that min x, T x,y (x) ≤ βx and t x,y (min x, T x,y (x) ) ≤ βx. Therefore, in F x,y (min x, T x,y (x) ) we have that ϕ y and ϕ y are evaluated on points smaller than βx, and so they are affine linear.
Assume that min x, T x,y (x) = T x,y (x). Notice that T x,y (x) is an affine combination of x and x as T x,y (x) = 1 1−λ y x + −λ y 1−λ y x and 1 1−λ y + −λ y 1−λ y = 1. As ϕ y is affine linear, we have that From ϕ y (x) = yϕ 1 (x) and the definition of λ y , we obtain Notice that z = ϕ 1 (x) (as ϕ 1 is the convex envelope of f 1 ) and so y).
is just a convex underestimator, which implies that φ S (x, y) ≤ F x,y (T x,y (x)). From here, we conclude that φ S (x, y) = F x,y (T x,y (x)) as we wanted to show.
The case of min x, T x,y (x) = x follows from a similar argument.
We now explain the behavior from Fig. 7.

Furthermore, vex D [ f ] is linear in A.
Proof Given that βx < x and T x,y (x) ≤ βx, we have that min x, T x,y (x) = T x,y (x) and t x,y (min x, T x,y (x) ) = x ≤ βx. Thus, we can apply the same reasoning as the one given in the above proof of Proposition 5.1.

Remark 5.1
A similar behavior can be observed between the concave overestimator induced by SCIP and the concave envelope of y sgn(x) |x| α .

Computational Study
In this section, we present a computational study that investigates the impact of the presented convex envelope on the performance of the expansion planning problem (see Sect. 2). For this, we extend SCIP with the convex envelope and compare it against default SCIP. We conduct two experiments to answer the following questions: 1. ROOTGAP: How much gap can be additionally closed in the root node of the branch-and-bound tree when using the envelopes of y sgn(x) |x| α with aggressive separation settings? 2. TREE: To which extent do the presented envelopes affect the performance of SCIP, both in running time and solvability?

Experimental Setup
In order to measure the impact of the convex envelope over the standard relaxation in terms of improving the dual bound, we use an aggressive emphasis setting for the separation in the ROOTGAP experiment. To reduce the impact of side effects, we disable restarts, 1 all primal heuristics, 2 and we give the value of the best known primal solution.
In contrast to the ROOTGAP experiment, the TREE experiment compares SCIP default against SCIP with the additional separation routine that generates gradient cuts to the envelopes of f in every local node during the entire tree search. We remark that these gradient cuts are, in general, only locally valid since we compute them using the local bounds of the node.
In both experiments, the gradient cuts for the constraint y sgn(x) |x| 2 = π v − π w using the convex and concave envelope read, respectively, and cave

Test Sets
The computational study is carried out on three test sets taken from [20]: Belgium, GasLib-40 and Circuit rank. These test sets are based upon different gas networks, demand scenarios and expansion candidates. In total, all test sets contain 6500 instances, with Belgium and GasLib-40 having each 2000 instances, and Circuit rank having 2500 instances. As we use gas instances, the exponent α is set to α = 2 in Eq. (1a). Networks The underlying networks of the test sets vary in size and structure. Here, we consider the rather small and almost tree-shaped Belgian gas network from the GAMS model library 3 and the larger GasLib-40 network from [29]. In general, optimization problems get more challenging to solve on potential-driven networks that have a higher circuit rank, which is the number |A| − |V| + 1, because the existence of cycles leads to more complex patterns of flow directions, see [31] and the references therein. To also account for this type of computational difficulty, the Circuit rank test set contains 1 In restarts, SCIP aborts the current search after encountering sufficient variable bound reductions and starts preprocessing the problem again. For more details about restarts in SCIP, we refer to Section 10.9 in [1]. 2 We used the following SCIP settings: limits/totalnodes = 1, separation/emphasis/aggressive = True, limits/restarts = 0 as well as heuristics/emphasis = Off. 3 General Algebraic Modeling System (GAMS) Model Library https://www.gams.com/modlib/modlib. htm. instances where the circuit rank is successively increased by adding up to ten new pipelines to the Belgian network, resulting in five new networks Belgium + (2, 4, 6, 8, 10) arcs. The rationale is to add arcs that connect different regions of the network and have a high impact on the network topology. Details about the different networks, their visualization and information of how to generate the new pipelines are found in Table 1, Figure 6 and Section 6 in [21]. We note that the original Belgian and GasLib-40 networks contain compressor stations. However, their treatment is outside the scope of this paper. In the case of the Belgian network, we modeled them as normal pipelines by using available data from the GAMS model library, while in the case of the GasLib-40 network, compressor station arcs had to be contracted by combining their adjacent nodes due to lack of available data. Demands It is known in practice that increasing network demand typically results in more complex transport situations that stress the networks. Therefore, we generated instances for the Belgium and GasLib-40 test sets that represent a wide range of possible network demand vectors b ∈ R |V| starting from "easier" instances with lower demands up to "harder" instances with higher demands. For the Circuit rank experiment, we utilize demand scenarios with a constant overall network demand. An overview of the total network demand loads (in [10 6 m 3 /day]) used for the particular networks is given in Table 1. Throughout, we randomly generated 500 scenarios for each network and for each given total network demand according to a random uniform distribution of the demand at sources and sinks. For more detailed information about the deployed method, we refer to [21]. Due to the variety of instances, we can assume that the results will provide us with a representative picture for the impact of the convex envelope on gas network related instances. Expansion capacities In all instances of the three test sets, the expansion capacities [y a , y a ] of all arcs a ∈ A result from building at most three pipelines in parallel to the already existing one. The data points for the cost function used follow a quadratic form, see Figure 1 in [21], and were provided to us by practitioners. Implementation We extended a development version 4 of SCIP by a so-called nonlinear handler, where the separation is applied in the separation and enforcement callbacks. With this handler, we add the cuts (21) and (22) in addition to the cuts that SCIP generates per default (see Sect. 5).
Given that the considered gas instances tend to be numerically unstable due to variables in different scales, we observed that many invalid cuts were generated in some preliminary experiments. The following two measures prevented the addition of invalid cuts to SCIP.
-We only add the cuts (21) and (22)  To safeguard against a potential mutual slowdown of parallel processes, we bind the processes to specific cores and run at most one job per node at a time. We used a development version of SCIP 6.0.2.4 [10] with CPLEX 12.7.1.0 as LP solver [8] and Ipopt 3.12.13 as NLP solver [36].

Computational Results
In this section, we present the results for the ROOTGAP and TREE experiments. In the following, Enabled refers to adding the gradient cuts of the convex and concave envelopes, as opposed to Disabled, which refers to SCIP default settings. ROOTGAP Experiment From the 6500 instances, we do not consider the ones that have been detected infeasible (396 instances), no primal solution is known (1113 instances), or none of the versions could increase the dual bound by more than 10 −6 (1847 instances). This restriction results in 3144 instances for the ROOTGAP experiment.
To compare the dual bounds of Enabled and Disabled relative to a given primal bound, we use the following measure. For an instance, let d 1 ∈ R be the dual bound of Enabled, and let d 2 ∈ R be the dual bound of Disabled. Furthermore, let p ∈ R be a reference primal bound, for example, the optimal or best known objective value of that instance. The function GC : measures the improvement of the gap closed when comparing the dual bounds d 1 and d 2 relative to p, see also [26]. Aggregated results for the ROOTGAP experiment are shown in Table 2 and Fig. 8. Table 2 shows that the addition of the cuts closes more gap on 1033 out of the 3144 instances and closes less gap only on 86 instances. (For an explanation, see below.) Additionally, the table presents the average gap closed for instances, where either Enabled improves (better) or deteriorates (worse) by more than 0%, 2%, or 10%, whereas change encompasses the union of those instances that are better or worse by more than 0%, 2%, or 10%.
Among all instances, the average gap closed improvement is 0.87%. However, when one considers all instances for which the gap closed is different, the gap closed improvement increases to 2.44%. Furthermore, if we consider instances for which the gap closed differs by 2% and 10%, then the gap closed improvement increases to 6.75% and 19.69%, respectively.
The presented results show that our convexification of the considered constraint function has a significant impact on the quality of the relaxation. In theory, the dual bounds obtained by Enabled should always be at least as good as the one from Disabled, since we apply the gradient cuts on top of SCIP default.
Nevertheless, as mentioned above, there are 86 where disabling the cuts yields improved dual bounds. The explanation for this behavior is that some of the gradient cuts from the convex envelope increase the computational cost of solving LPs in the OBBT propagator [11]. As a consequence, the internal limit of simplex iterations in OBBT are hit faster, which has the side effect that less genvbounds can be learned. However, less genvbounds reduce the number of variable bound reductions, and as a consequence, lead to worse dual bounds. Indeed, when turning off the propagator genvbounds for Enabled and Disabled on these 86 instances, then better dual bounds are obtained by Enabled for all, but 14 out of these 86 instances. Furthermore, if we disable all propagation methods, then the behavior of Enabled is at least as good as Disabled on these 14 instances. 5 TREE Experiment Aggregated results for the TREE experiment are shown in Tables 3, 4 and 5 for the three test sets under consideration. We compare the number of solved instances and the solver performance when activating the cut generation for the convex envelope (Enabled) with SCIP using default settings (Disabled). The number of solved instances and the computation time is split into different categories: ALL, ALL OPT, [1, tlim], [10, tlim], [100, tlim] and [1000, tlim]. ALL consists of all instances. ALL OPT consists of the instances that are solved to optimality or detected to be infeasible by both settings. [t, tlim] consists of the instances that could be solved by at least one of both settings, Enabled or Disabled, in more than t seconds within the time limit of one hour. Note that the instances in [t 1 , tlim] are contained in [t 2 , tlim] whenever t 1 ≥ t 2 and that one can consider the instances in [t, tlim] to be harder to solve the larger t is.
Most importantly, on all test sets, enabling the cut generation for the convex envelope increases the number of solved instances and decreases the average solving time compared to SCIP default. Concerning the number of solved instances, Enabled solves 39 more instances on Belgium, 19 more instances on GasLib-40 and 68 more instances on Circuit rank. From the tables, we deduce that all the extra instances that Enabled solves belong to the groups [1000, tlim]. This means that our separation routine renders these instances solvable but not trivially solvable. The column "Relative" reports the change of solving time and branch-and-bound nodes relative to Enabled The column "Relative" reports the change of solving time and branch-and-bound nodes relative to Enabled To compare the total running time between Enabled and Disabled, we use the shifted geometric mean (sgm) with a shift value of s = 10 for the average solving time in seconds. For instances that are solved to optimality by both settings (ALL OPT), we additionally present the sgm for the number of explored branch-and-bound nodes with a shift value of s = 100. As can be seen in Tables 3, 4 and 5, there is an improvement in the performance when using Enabled over SCIP default in all test sets and through all time categories. This improvement increases toward the more "difficult" instances to a speed-up of up to 58% on Belgium. Even for ALL OPT, there is throughout a speed-up of 19% (Belgium) and 12% (GasLib-40 and Circuit rank) when using the Enabled separation routine. Likewise, the number of branchand-bound nodes significantly decreases by 22%, 25% and 14% for ALL OPT on these three test sets.
Finally, we remark that the Wilcoxon signed-rank test [39] judges the performance improvements on the Belgium, GasLib-40 and Circuit rank test sets as highly statistical significant being consistently distributed in the shifted geometric mean across all three test sets. 6

Conclusion
In this paper, we derived the convex envelope of the nonconvex and bivariate function f (x, y) = y sgn(x) |x| α for α > 1. This contribution was motivated by the fact that strong relaxations of nonconvex constraints play a key component in solving MINLPs, see [12]. As this function frequently occurs in expansion planning problems, it is also of significant practical relevance. We performed an extensive computational study on real-world instances which showed the benefit of the convex envelope over the standard relaxation used in state-of-the-art solvers. Given that we provide a closedform expression for the envelopes, our implementation in the MINLP solver SCIP allows us to calculate and exploit the envelopes at every node in the branch-andbound tree, which has a crucial impact on the solving process. Apart from yielding improved dual bounds that reduce the gap already in the root node, our procedure drastically boosts the solving process in the branch-and-bound tree. Our procedure significantly reduces the solving time on all test sets, for example, by up to 58% for hard instances on Belgium. Even more important, this tighter relaxation enables to solve more instances on all test sets.
One possible avenue for further research is to consider the convexification of multiple constraints simultaneously. In this fashion, [22] already reported promising results for quadratic functions with absolute value terms for gas network related problems. 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/.