Optimization of capacity expansion in potential-driven networks including multiple looping: a comparison of modelling approaches

Lenz, Ralf; Becker, Kai Helge Article — Published Version Optimization of capacity expansion in potential-driven networks including multiple looping: a comparison of modelling approaches OR Spectrum Provided in Cooperation with: Springer Nature Suggested Citation: Lenz, Ralf; Becker, Kai Helge (2021) : Optimization of capacity expansion in potential-driven networks including multiple looping: a comparison of modelling approaches, OR Spectrum, ISSN 1436-6304, Springer, Berlin, Heidelberg, Vol. 44, Iss. 1, pp. 179-224, https://doi.org/10.1007/s00291-021-00648-7 This Version is available at: https://hdl.handle.net/10419/286855 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/ Vol.:(0123456789) OR Spectrum (2022) 44:179–224 https://doi.org/10.1007/s00291-021-00648-7 1 3 ORIGINAL ARTICLE Optimization ofcapacity expansion inpotential‑driven networks includingmultiple looping: acomparison ofmodelling approaches RalfLenz1 · KaiHelgeBecker1,2,3 Received: 18 September 2019 / Accepted: 26 July 2021 / Published online: 19 November 2021 © The Author(s) 2021 Abstract In commodity transport networks such as natural gas, hydrogen and water networks, flows arise from nonlinear potential differences between the nodes, which can be represented by so-called potential-driven network models. When operators of these networks face increasing demand or the need to handle more diverse transport situations, they regularly seek to expand the capacity of their network by building new pipelines parallel to existing ones (“looping”). The paper introduces a new mixedinteger nonlinear programming model and a new nonlinear programming model and compares these with existing models for the looping problem and related problems in the literature, both theoretically and experimentally. On this basis, we give recommendations to practitioners about the circumstances under which a certain model should be used. In particular, it turns out that one of our novel models outperforms the existing models with respect to computational time, the number of solutions found, the number of instances solved and cost savings. Moreover, the paper extends the models for optimizing over multiple demand scenarios and is the first to include the practically relevant option that a particular pipeline may be looped several times. Keywords Nonlinear programming· MINLP· Potential-driven networks· Network expansion· OR in energy * Ralf Lenz [email protected] Kai Helge Becker [email protected] 1 Zuse Institute Berlin, Takustr. 7, 14195Berlin, Germany 2 Berlin School ofEconomics andLaw, Alt-Friedrichsfelde 60, 10315Berlin, Germany 3 FOM Hochschule Berlin, Bismarckstrasse 107, 10625Berlin, Germany 180 R.Lenz, K.H.Becker 1 3 1 Introduction Operators of commodity transport networks such as natural gas, hydrogen and water networks regularly have to face both increasing demand and the need to handle more diverse transport situations. To deal with these challenges without having to resort to the expensive options of setting up new pipeline corridors or of demolishing pipelines and replacing them by larger ones, they must expand the capacity of the existing network. Without using compressor stations to increase the pressure in the pipelines, this can only be achieved by building pipelines parallel to existing ones, a process that is referred to as “looping” in the terminology of the industry. Due to the high costs involved, there has long been research about finding the cost minimal way of looping an existing network, i.e. of determining both the pipelines that are to be looped and the diameters of the pipelines to be built. Commodity flows such as flows in gas and water networks arise from frictioninduced potential differences between the nodes in the network, which leads to a type of nonlinear network model that is referred to as potential-driven network in the literature (Birkhoff and Diaz 1956; Raghunathan 2013; Robinius etal. 2018). Apart from the nonlinear character of potential-driven networks, a specific difficulty in finding optimal loops arises from the fact that the diameters of the pipelines typically have to be selected from a discrete set of commercially available diameters (André etal. 2009; Hansen etal. 1991; Fasold 1999), which additionally imparts to the problem a combinatorial flavour. For these reasons, the problem of capacity expansion in potential-driven networks belongs to the family of mixed-integer nonlinear programming (MINLP) problems. As a consequence, a large body of the literature on the topic is concerned with developing sophisticated special-purpose algorithms geared at finding locally optimal solutions or approximating a global optimum, partly based on novel ways of modelling the capacity expansion problem. While this research focus has certainly advanced our ability to solve this problem, it has also led to a variety of mathematical models for the looping problem that remain unconnected in the literature. In view of the recent advances in the development of general-purpose solvers for NLP and MINLP problems that allow for an efficient implementation of models that is significantly faster than implementing one of the specialpurpose algorithms developed, it seems to be useful for practitioners involved in the design of networks to shift the research focus to the modelling stage. In fact, despite the algorithmic advances in the previous decades, practitioners tackling the capacity expansion problem often find it more convenient to resort to simulation studies instead of implementing the complex special-purpose optimization models proposed in the literature, thereby not realizing the potential that an optimization-based approach has to offer. For this reason, the present paper discusses various models for finding optimal loops in potential-driven networks that can be implemented easily by using state-of-the-art MINLP solvers. In particular, this paper brings together, for the first time, the diversity of existing models in the literature and compares these, both theoretically and 181 1 3 Optimization ofcapacity expansion inpotential‑driven… experimentally. For this comparison, we choose a comprehensive approach that includes two different, rather unconnected strands of research on the problem that are referred to as the discrete approach and the split-pipe approach (see below), and also takes into account literature that does not explicitly focus on the looping problem, but addresses closely related problems. Moreover, the present paper also introduces two novel models for the looping problem, one of which has turned out to outperform the existing models in the literature with respect to computational time, the number of solutions found, the number of instances solved and cost savings. Additionally, we extend our models to include the case of optimizing capacity expansion across multiple demand scenarios. Finally, in discussing these models, the paper goes beyond the existing literature by including, throughout the paper, the practically highly relevant case that a particular pipeline may be looped several times. Altogether, by focussing on the modelling aspect in conjunction with state-ofthe-art MINLP solvers, the present paper should provide a useful guide for practitioners who look for an optimization-based alternative to simulation studies, but do not have the option to implement the special-purpose algorithms suggested in the literature. Additionally, our overview and comparison of modelling approaches in the literature, as well as our new models, may be interesting for researchers who seek to develop new special-purpose algorithms for the capacity expansion problem. In line with our aim to bring together the diversity of existing approaches in the literature, the remainder of this section discusses the literature on the topic rather comprehensively before presenting an outline of the structure of the paper. During the past 40 years of research on capacity problems in potential-driven networks, three different modelling approaches have been considered for the looping problem (cf. Shiono and Suzuki 2016): (1) a direct approach where the optimal diameters for the loops are chosen from the set of commercially available diameters (discrete approach), (2) a continuous approach where continuous diameters are typically used to approximate the problem and (3) an extended approach where the entire length of the pipeline to be looped is split into several segments of variable lengths each of which may have its own diameter from the discrete set of available diameters (split-pipe approach). In the following discussion of the literature about these three approaches, we will not only look at papers addressing the looping problem (in fact, to the best of our knowledge, there are so far, apart from the present paper, only two papers, namely André etal. (2009) and Pietrasz etal. (2008), that solely focus on the looping problem), but will also consider papers about two closely related problems: the network design problem, which looks for the optimal diameters of pipelines between given unconnected nodes or determines both the location of nodes and the optimal pipelines between these, and the network expansion problem, which is about the optimal placement of new network elements of different types (such as pipelines, compressors and valves) at pre-defined (previously connected or unconnected) locations in the network. (1) Discrete diameters One of the first solution approaches using mathematical optimization techniques for this NP-hard problem (Yates et al. 1984) has been proposed by Jacoby (1968). The author solves a nonlinear model using a 182 R.Lenz, K.H.Becker 1 3 gradient approximation method where the resulting continuous diameters are finally rounded to the nearest discrete-valued pipe diameters. The approach was tested on a small water network that contains seven pipelines and two cycles. Other early work where the discrete approach is used for networks with a simple structure include Liang (1971), who solved a gunbarrel system using dynamic programming, and Rothfarb etal. (1970), who used a serial and parallel merge algorithm to design tree shaped networks. Later, Gessler (1985) applied selective enumeration techniques to tackle a small sized network with two cycles. In the 1990s, a class of approaches was developed that relies on meta-heuristics, such as genetic algorithms, see, e.g. Simpson etal. (1994) and Savic and Walters (1997) for water networks, and Boyd etal. (1994) and Castillo and González (1998) for gas networks. In the past decade, different papers applied a MINLP formulation to determine discrete pipe sizes, e.g. André etal. (2009), Bragalli etal. (2012) and Robinius etal. (2018). While André etal. (2009) solve the problem heuristically in two stages, where a first step identifies pipes to be looped by solving the continuous relaxation and a second step determines discrete-valued diameters for these selected pipes using Branch and Bound, Bragalli etal. (2012) solve the model directly with an MINLP solver using a continuous reformulation of the cost function. Robinius etal. (2018) determine discrete arc sizes in the context of tree shaped networks, which allows flows on arcs to be fixed and thus simplifying the MINLP model to an Mixed Integer Programming model. The nonlinear nature of the problem has led to a number of different MINLP models and approaches: Raghunathan (2013) presents a disjunctive program together with a convex relaxation, which is then solved to global optimality and Borraz-Sánchez etal. (2016) propose a new solution approach by presenting a new model together with a second-order cone relaxation. Humpola (2014) formulates a model that is solved to global optimality using convex reformulation techniques, special tailored cuts. In this paper, we investigate different existing MINLP approaches and propose a new model for the discrete looping problem. (2) Continuous diameters have been considered e.g. by Bhaskaran and Salzborn (1979a), Rowell and Barnes (1982), Bhave (1985), De Wolf and Smeers (1996), DeWolf and Bakhouya (2012), Babonneau etal. (2012) and André etal. (2013). Continuous diameters typically have been used in the literature to approximate the discrete diameters that are commercially available. Hansen etal. (1991), for example, use successive linear programming with a trust region strategy, where the algorithm adjusts the continuous diameters in each iteration to elements in the set of available discrete diameters. Osiadacz and Górecki (1995) apply sequential quadratic programming to the continuous relaxation of a gas network design problem and round the solution to the closest available diameter size. Shiono and Suzuki (2016) introduce an analytic approach that calculates the optimal diameter costs for the pipe-sizing problem of a tree-shaped gas network with continuous diameters, which are then heuristically converted to discrete pipe diameters. As the present paper is concerned with models that lead to exact solutions, we will not further consider this line of research here. 183 1 3 Optimization ofcapacity expansion inpotential‑driven… (3) Split-pipe approach This approach, where the pipes can be split into several segments with different diameters each, combines features of both the discrete and the continuous approaches: while the diameters are chosen from the discrete set of available diameters, the option to split pipes at arbitrary points into sections of different diameters leads, as we will see in the next section, to a situation that is equivalent to choosing diameters from a continuous set. This concept was first used in designing networks with a tree structure, which allows the flows in the pipelines to be treated as constants and leads to a linear programming model (e.g. Karmeli etal. (1968) and Gupta etal. (1972) for water networks, and, independently, Bhaskaran and Salzborn (1979b) for gas networks). The first paper to attempt the split-pipe version of the capacity expansion problem as a nonlinear problem for general networks was Alperovits and Shamir (1977), who introduced the linear programming gradient (LPG) method, a two-stage heuristic that alternates between solving the linear program with fixed flows from Karmeli etal. (1968) for obtaining the pipeline diameters and modifying these flows on the basis of a sensitivity analysis of the solution of the first stage. This idea stimulated a number of subsequent papers (e.g. Quindry etal. 1981; Fujiwara and Khang 1990; Kessler and Shamir 1989) that improved on the method. Starting with Eiger etal. (1994), a strand of genuinely NLP-based global solution methods emerged, with further contributions by Zhang and Zhu (1996) and Sherali etal. (2001), for example. Surprisingly, despite 50 years of research on the split-pipe approach, all these contributions proceed from the very same basic model that goes back to the linear model by Karmeli etal. (1968). For this reason, the present paper presents a novel, alternative model for the split-pipe approach and compares it to the split-pipe model in the literature. We have now sketched the history of research on the three different approaches to our looping problem and closely related problems. In the light of the variety of models and optimization methods for these problems, it is astonishing that there is nearly no research that brings together the discrete and the split-pipe approaches and compares the models proposed in the literature. In fact, to the best of our knowledge, the very early paper by Bhaskaran and Salzborn (1979b) with its linear model for a tree-shaped network is the only paper to compare the split-pipe approach with the discrete approach, albeit for very small trees with 9 and 14 nodes. The present paper addresses this gap and presents comparisons of all models discussed in the literature. Another important aspect of the looping problem that has not been addressed in the literature is multiple looping, where each pipeline may be reinforced with several parallel pipelines of different diameters. This is a practically highly relevant problem as multiple looping can first replace large diameters that are commercially not available; may second lead to cost savings by substituting several parallel pipelines with smaller diameters for one pipeline with a large diameter; may third, as we will see in Sect.2.3, allow for pipe characteristics that cannot be realized with single loops; and can finally provide a tool for strategic planning where several stages of successively looping a given network are to be considered. For these reasons, we will take into account multiple looping throughout the paper. The remainder of the paper is organized as follows: In the next section, we formalize the looping problem in potential-driven networks, show some of its basic properties and explain our approach of dealing with multiple loops. In Sect.3, we present a 184 R.Lenz, K.H.Becker 1 3 new model for the discrete looping problem and contrast it with the existing models in the literature. The subsequent Sect.4 turns to the split-pipe approach of the looping problem. Again we present a new model and address its relationship with the splitpipe model that can be found in the literature. Moreover, we discuss the way in which the feasible regions of all models presented so far are related to each other. Section5 extends the models to the case of multi-scenarios. In Sect.6, we carry out extensive computational experiments on instances of both natural gas and water networks that allows us to compare the performance of all models and give recommendations regarding their use. The paper ends with some concluding remarks in Sect.7. 2 The expansion planning problem Let us begin by formally defining the planning problem that this paper discusses. 2.1 Problem statement Let G=(V,A,r) be a directed multigraph with node set V , arc set A and a function r∶A→V×V that maps each arc to its end points. The nodes can be partitioned into supply nodes (sources), consumption nodes (sinks) and transshipment nodes. In this paper, we restrict to passive and connected networks, where the only arc type are pipes that are needed to transport the commodities. We are given a demand vector b∈ℝ|V| where bv>0 denotes injection into the network at source v∈V and bv<0 withdrawal at sink v∈V . Since we work with stationary and isothermal models, the demand vector is balanced, i.e. ∑v∈V b v = 0 . With each arc a∈A , flow variables xa∈[x a,xa] are associated. Positive values of xa indicate flow along the arc a=(v,w) from v to w, whereas negative values indicate flow in the reversed direction. As in classical network flow problems, flow conservation is required at every node i.e. ∑a∈𝛿 + (v) x a − ∑a∈𝛿 − (v) x a =b v with 𝛿+(v) ∶= {(v,w)∈A} and 𝛿−(v) ∶= {(w,v)∈A} . In potential-based networks, the physical state is additionally described by nonnegative potential variables 𝜋v∈[𝜋 v,𝜋v] at each node v∈V . In applications such as water transport problems, the nodal potentials 𝜋v correspond to hydraulic heads, i.e. the sum of the elevation head, velocity head and pressure head (Walski etal. 2001), while they represent squared pressure variables in gas transport problems (Koch etal. 2015). The flow along a pipe a=(v,w) depends on the potential difference at its adjacent nodes, the pipe length La>0 , a physical parameter Ra>0 representing phenomena such as friction and density and the diameter da>0 . The potential difference is given by a function Φ(da,.)∶ℝ → ℝ that is strictly increasing, Φ(da,.)∈ C 1 and anti-symmetrical, i.e. Φ(da,−xa) = −Φ(da,xa) , and typically takes the form (1) 𝜋v−𝜋w= Φ(da,xa)∀a∈Aand (v,w)=r(a), 185 1 3 Optimization ofcapacity expansion inpotential‑driven… with 𝛼,𝛽>0 . In line with most authors in the literature, we model Ra as a constant and assume the pipes to have zero slope, i.e. there is no influence of gravity on the potential drop 𝜋v −𝜋w (e.g. Zhang and Zhu 1996; André etal. 2009; Babonneau etal. 2012). The value of 𝛼 depends on the commodity and the type of approximation used for modelling the commodity flow. In water transport problems, the potential loss function is typically given by the equations of Darcy–Weisbach with 𝛼=2 or Hazen–Williams with 𝛼=1.852 (Walski etal. 2001), whereas in gas transport problems Eq. (2) takes the shape of the Weymouth equation with 𝛼=2 ( Weymouth (1912)). In the approximations proposed by Darcy–Weisbach and Weymouth, the exponent of the diameter is 𝛽=5 , while it is 𝛽=4.87 in the case of the approximation by Hazen–Williams. The solution of the capacity expansion problem involves two decisions: (a) which pipelines a∈A should be looped? and (b) what pipeline diameters da should be used for the loops? As we wish to solve the problem to global optimality and any pre-selection of certain pipes are looping candidates would be a heuristic procedure, we will allow all pipes to be looped. For the diameters da , we have in the discrete case da∈Da∶= {Da,0,Da,1, ..., D a , k a} , where da,0 refers to the diameter of the already existing pipe and d a , k a denotes the maximal possible diameter when looping. In the continuous case the domain of the diameter is given by d a ∈D a ∶= [D a ,0,D a , k a] . While the present paper is not concerned with the continuous looping problem due to its approximative nature, we will see in Sect.2.3 that a continuous interval of diameters can also be interpreted as representing diameters in the split-pipe problem. The general looping problem then reads Note that even when no particular bounds x a and xa on the flow variables xa are given, flow bounds are implied by the bounds on 𝜋v and da by virtue of Eq.(1). (2) Φ( da,xa)= L a R a d 𝛽 a sgn (xa) | xa | 𝛼 . (3a) min x a,da,𝜋v ∑ a∈A c(da)L a (3b) s.t. 𝜋v−𝜋w= Φ(da,xa)∀a∈Aand (v,w)=r(a) (3c) ∑ a∈𝛿 + (v) xa− ∑ a∈𝛿 − (v) xa=bv∀v∈V (3d) 𝜋 v≤𝜋v≤𝜋v∀v∈V (3e) x a≤xa≤xa∀a∈A (3f) da∈Da∀a∈A. 186 R.Lenz, K.H.Becker 1 3 To provide the reader with some intuition for the problem statement of the expansion planning problem, let us briefly demonstrate that we can always find a solution to the problem provided that there is no flow bound (3e), that the intersection of the bounds of the potential variables is non-empty, i.e. ∩v∈ V [𝜋 v,𝜋v]=[𝜋,𝜋] with 𝜋≤𝜋 and that we can choose sufficiently large pipeline diameters. Let b the overall network inflow, i.e. the sum of all flows that enter the network at the entry nodes, which is clearly an upper bound on any flow along any arc in the network. We now select for each arc a∈A a diameter  da such that La R a ∕  d 𝛽 a sgn (b) | b |𝛼 ≤(𝜋−𝜋)∕ | A | . Let (𝜋v,xa)v∈V,a∈A be a corresponding solution of Problem(3a–c) where the values for  da are fixed. The solution has unique flow values and the potential values are uniquely determined up to a constant shift (cf. Collins etal. 1978; Humpola 2014 for the existence of such a solution), which obviously satisfies Summing up these inequalities along any path ������ v0vk between two connected nodes v0 and vk yields: We now choose a node v0 with the highest potential in our solution and shift this node’s potential by a constant r∈ℝ such that it is equal to the upper bound 𝜋 , i.e. such that we have 𝜋 v0 ∶= 𝜋v0 +r=𝜋 for the new potential of the node v0 . If we shift the potentials 𝜋v of all other nodes by the same constant r, all our potentials are guaranteed to satisfy (3d). 2.2 Convexity analysis In the previous section, we have seen that we can model the expansion planning problem as a Mixed-Integer Nonlinear Program (MINLP) for discrete looping decisions and we will see in Sect.2.3 that the split-pipe approach leads to a Nonlinear Program (NLP). We will now show that the feasible regions of the discrete and continuous (or split-pipe) capacity expansion planning problem are non-convex. When neglecting loop expansions, i.e. da=da,0 for all pipes a∈A , the continuous and the discrete problems (3) reduce to the same existence problem of validating a given demand scenario for feasibility. The feasible region of the resulting problem is convex (Maugis 1977; Collins etal. 1978), even though it comprises nonlinear nonconvex constraints of type(3b). However, the following proposition shows that this property does not hold for the feasible region of the expansion planning problem. 𝜋 v−𝜋w= LaRa  d𝛽 a sgn (xa)|xa|𝛼≤ LaRa  d𝛽 a b𝛼 ≤ 𝜋−𝜋 |A|∀a∈Aand (v,w)=r(a) . | 𝜋v0−𝜋k|=| ∑ k−1 i=0𝜋vi−𝜋vi+1|=| ∑ k−1 i=0 L i,i+1 R i,i+1  d𝛽 i,i+1 sgn (xi,i+1)|xi,i+1|𝛼 | ≤ ∑ k−1 i=0 Li,i+1Ri,i+1  d𝛽 i,i+1 b𝛼 ≤ ∑ k−1 i=0 𝜋−𝜋 | A | k≤|A| ≤𝜋−𝜋. 193 1 3 Optimization ofcapacity expansion inpotential‑driven… a one-dimensional function 𝜙(xa) . This is achieved by shifting the loop diameters from the term that represents the flow towards the potential loss term (see Eq.9b). To this end, we introduce variables Δa,i to model the potential loss that corresponds to the loop diameter chosen (constraint 9c) and use the binary variables 𝜆a,i not to choose the diameters, but to select the potential loss. As a trade-off we have to introduce big-M constraints (9d) to determine the potential loss variables Δa,i . where M(1) a,i ∶= Φ(Da,i,x a ),M (2) a,i ∶= Φ(Da,i,xa ) . 3.4 Further approaches fordiscrete looping intheliterature There are two further models in the literature for the network design problem, both of which are not intended as “stand-alone models”, but as starting points for solution algorithms. While we will not go into the details of the models here because preliminary computational experiments have demonstrated that they do not perform particularly well (see Sect.6.2), the approaches still deserve mentioning. Borraz-Sánchez etal. (2016) present an MINLP formulation that is similar to the model developed by Raghunathan (2013) (Sect.3.2) in the sense that it distinguishes (9a) min 𝜆,x,𝜋,Δ ∑a∈A La ∑k a i=0 𝜆a,ica, i (9b) s.t. sgn (xa) | xa | 𝛼− ∑ ka i=0 D 𝛽 a,i LaRa Δa,i=0∀a∈A , (9c) ∑k a i=0 Δa,i=𝜋v−𝜋w∀a∈A,(v,w)=r(a) , (9d) M(1) a,i 𝜆a,i≤Δa,i≤M (2) a,i 𝜆a,i∀a∈A,∀i∈[ka] , (9e) ∑k a i=0 𝜆a,i=1∀a∈A , (9f) ∑ a∈𝛿 + (v) xa− ∑ a∈𝛿 − (v) xa=bv∀v∈V , (9g) 𝜋 v ≤ 𝜋v ≤ 𝜋v∀v∈V, (9h) 𝜋v−𝜋w ≤ Δa,i ≤ 𝜋v−𝜋w∀a∈A,∀i∈[ka], (9i) x a ≤x a ≤x a∀ a ∈A, (9j) 𝜆a,i∈{0, 1}∀ a ∈A,∀ i ∈[ k a], 194 R.Lenz, K.H.Becker 1 3 between forward and backward flow directions. Here, however, this is achieved by binary variables z+ a,z− a that are then used to model the potential loss function in product form 𝜆 a,i(z+ a −z− a )(𝜋v−𝜋w)=LaRa∕D 𝛽 a,i x 𝛼 . Humpola (2014) uses indicator constraints to select arcs a in the network design problem (or diameters Da,i in the context of our expansion problem), i.e. constraints of the form 𝜆a,i=1 ⇒ 𝜋v−𝜋w= Φ(da,i,xa) and 𝜆a,i=0 ⇒ xa,i=0 , which may be represented by big-M constraints in a MINLP model. 3.5 Equivalence ofthediscrete models To finish our section on discrete models for the capacity expansion problem, we show that the three models presented so far are equivalent. More precisely, we show that for each solution to one of the models, there exist solutions to the other two models such that all of them represent the (1) same physical network state in terms of flow and pressure, (2) same decisions on loop extensions, and (3) same objective function values. To this end, we use the standard projection to map the feasible regions of the discrete models onto the space of their common variables, i.e. onto the ( (𝜋 v ,𝜆 a,i ) v∈V,a∈A,i∈[ka]) -space: In this context we denote the feasible regions of the discrete models as: XA for ModelA, XB for ModelB and XC for ModelC. Proposition 2 The discrete ModelsA, B and C are equivalent, i.e. Proof “ X A ⊆X C ”: Let ( x a ,𝜋 v ,  𝜆 a,i ) a∈ A ,v∈ V ,i∈[ka] be a solution of Model A, i.e. ( 𝜋 v ,  𝜆 a,i ) a∈ A ,v∈ V ,i∈[ka] ∈X A . Then (9e–g) and (9i,j) follow right away with ( x a ,𝜋 v ,  𝜆 a,i ) a∈ A ,v∈ V ,i∈[ka] . For pipe a, let  ia∈[ka] be such that  𝜆a,  ia = 1 and define Then Eq.(9b) follows by construction, (9c) equals (7b) and (9d) holds by virtue of M(1) a,i = Φ(Da,  i a ,x a ),M (2) a,i = Φ(Da,  i a ,xa ) and (9h) follows from (9c) and (9g). Thus, ( x a ,𝜋 v ,  𝜆 a,i ,Δ a,i ) a∈ A ,v∈ V ,i∈[ka] is a solution of ModelC, and ( 𝜋 v ,  𝜆 a,i ) v∈ V ,a∈ A ,i∈[ka] ∈X C . “ X C ⊆XA ”: Let ( xa,𝜋v,𝜆a,i,Δa,i ) a∈A,v∈V,i∈[k a] be a solution of Model C, i.e. ( 𝜋v,𝜆a,i )v∈ V ,a∈ A ,i∈[ka] ∈X C . With ( xa,𝜋v,𝜆a,i,Δa,i )a∈ A ,v∈ V ,i∈[ka] , Eq. (7b) follows from (9b–d) and (7c–g) equals (9e–g), (9j–i). Hence, ( 𝜋v,𝜆a,i )v∈ V ,a∈ A ,i∈[ka] ∈X A . XA=Proj 𝜋v,𝜆a,i {( 𝜋v,𝜆a,i,xa ) | ( 𝜋v,𝜆a,i,xa ) satisfy Eqs. (7b)−(7g) } , XB=Proj 𝜋v,𝜆a,i{(𝜋v,𝜆a,i,x+ a,i,x− a,i,Δ+ a,i,Δ− a,i,za)|(𝜋v,𝜆a,i,x+ a,i,x− a,i,Δ+ a,i,Δ− a,i,za) satisfy Eqs. (8b)−(8k)}, X C=Proj 𝜋v,𝜆a,i {( 𝜋v,𝜆a,i,xa,Δa,i )|( 𝜋v,𝜆a,i,xa,Δa,i ) satisfy Eqs. (9b)−(9j) }. XA=XC=XB. Δ a,i∶= { 0∀i∈[ka]⧵{  ia } LaRa D𝛽 a,i sgn (xa) | xa | 𝛼i= ia. 195 1 3 Optimization ofcapacity expansion inpotential‑driven… “ X A ⊆X B ”: Let ( 𝜋 v ,  𝜆 a,i ,x a ) v∈ V ,a∈ A ,i∈[ka] be a solution of Model A, that is, ( 𝜋 v ,  𝜆 a,i ) v∈ V ,a∈ A ,i∈[ka] ∈X A . For pipe a let  ia∈[ka] be such that  𝜆a,  ia = 1 . We set x+ a,i ,x − a,i ,Δ + a,i ,Δ − a,i ∶= 0 for all i∈[ka]⧵{  ia} . Case 1: if xa ≥ 0 , define x+ a,  ia ∶= xa,x − a,  ia = 0 and Δ+ a,  ia ∶= 𝜋 v−𝜋 w and za∶= 1 . Case 2: if xa<0 , define x + a,  ia =0, x − a,  ia ∶= −  x a and Δ− a,  ia ∶= 𝜋 w −𝜋 v and za∶= 0 . Then in both cases (8b)–(8k) hold. Hence, ( 𝜋 v ,  𝜆 a,i ) v∈ V ,a∈ A ,i∈[ka] ∈X B . “ X B ⊆X A”: Let ( x+ a,i,x− a,i,za,𝜋v,𝜆a,i,Δ+ a,i,Δ− a,i ) a∈A,v∈V,i∈[k a] be a solution of Model B, i.e. (𝜋 v ,𝜆 a,i ) v∈V,a∈A,i∈[k a ] ∈XB . Set  x a∶= ∑i∈ka (x + a,i −x − a,i) for all a∈A , then the equations of ModelA are satisfied. Hence, (𝜋 v ,𝜆 a,i ) v∈V,a∈A,i∈[k a ] ∈XA . ◻ 4 Split‑pipe loop expansions In this section, we present two equivalent modelling approaches for continuous loop expansions. The first one is a model common in the literature (e.g. Alperovits and Shamir 1977; Zhang and Zhu 1996), while the second one is a novel approach. 4.1 Split‑pipe looping withlength variables (Model D) This split-pipe model is identical to the first discrete Model A, except that the indicator variables 𝜆a,i∈{0, 1} have been relaxed to continuous length variables 𝜆a,i∈[0, 1] here, turning the MINLP of Sect. 3.1 into a NLP. As explained in Sect.2.3, these length variables denote the proportion of pipe a that is looped with diameter Da,i , i.e. for each arc a, the pipeline consists of segments of those equivalent diameters Da,i for which 𝜆a,i>0 . 4.2 Split‑pipe looping withefficient frontier constraints (Model E) In the previous section, we modelled the efficient frontier of equivalent diameters using length variables 𝜆a,i to express convex combinations of the equivalent (10a) minimize 𝜆,x,𝜋 ∑a∈A La ∑k a i=0 𝜆a,ica, i (10b–f) s.t. (7b)−(7f) (10g) 𝜆a,i∈[0, 1]∀a∈A,∀i∈[ka] 196 R.Lenz, K.H.Becker 1 3 diameters Da,i that are the extreme points of the frontier. In the present section, the efficient frontier is modelled explicitly by using linear constraints. As the efficient frontier consists of points (D−𝛽 ,c) (cf. Sect.2.3), we introduce new continuous variables ya to model the exponentiated diameter (see constraints(11b, g) and variables ca that represent the costs per unit length of pipe a for the equivalent diameters (see (11h) and the objective function). On this basis, the efficient frontier can be represented by linear constraints(11d), each of which models the frontier between a pair of adjacent extreme points (cf. the right side of Fig.1). The parameters sa,i and ta,i can be calculated in advance as s a,i=(ca,i−ca,i+1)∕(D −𝛽 a,i −D −𝛽 a,i+1) and t a,i=−sa,iD −𝛽 a,i+1 +ca,i+ 1 . Note that the bounds in (11g) can be calculated as ya =D −𝛽 a,ka and ya=D −𝛽 a,0 . 4.3 Equivalence ofthesplit‑pipe models To compare the feasible regions of the split-pipe ModelsDandE, we use the standard projection to map their feasible regions onto the space of their common variables, i.e. onto the ( (𝜋 v ,x a ) v∈V,a∈A) -space: (11a) min y,c,x,𝜋 ∑a∈A Lac a (11b) s.t. 𝜋v−𝜋w=LaRayasgn (xa)|xa|𝛼∀a∈A,(v,w)=r(a) (11c) ∑ a∈𝛿 + (v) xa− ∑ a∈𝛿 − (v) xa=bv∀v∈ V (11d) ca ≥s a,i y a +t a,i ∀a∈ A ,∀i∈[k a −1 ] (11e) 𝜋 v ≤𝜋v≤𝜋v ∀ v ∈V (11f) x a ≤xa≤xa∀a∈ A (11g) y a ≤ya≤ya∀a∈ A (11h) ca ≥ 0 ∀a∈ A X D=Proj 𝜋v,xa {( 𝜋v,xa,𝜆a,i )|( 𝜋v,xa,𝜆a,i ) satisfy Eqs. (10b)−(10g) } , X E=Proj 𝜋v,xa {( 𝜋v,xa,ya,ca )|( 𝜋v,xa,ya,ca ) satisfy Eqs. (11b)−(11h) }. 197 1 3 Optimization ofcapacity expansion inpotential‑driven… The feasible region of ModelD is a subset of the feasible region of ModelE due to the inequality constraints(11d), i.e. XD⊆XE . But since the objective is to minimize the cost for building loops, solutions of both models are forced to be on the efficient frontier. The proof of the following proposition formalizes this argument. Proposition 3 Each optimal solution to ModelD is an optimal solution to ModelE, and vice versa. Proof (i) Let (x∗ a,𝜋∗ v,𝜆∗ a,i) a∈A,v∈V,i∈[k a] be an optimal solution of Model D, i.e. ( x ∗ a,𝜋 ∗ v ) a ∈ A,v ∈ V,i ∈[ k a] ∈X D . According to Fujiwara and Dey (1987), an optimal solution has the property that each arc consists of at most two pipe segments, where the corresponding equivalent diameters correspond to adjacent extreme points of the efficient-frontier {( D −𝛽 a,0 ,ca,0),…,(D −𝛽 a , k a ,ca,k a)} . Hence, for all a∈A there exists a unique i0∈[ka−1] such that 0 ≤ 𝜆∗ a,i0 ,𝜆∗ a,i0+1 ≤ 1 with 𝜆∗ a,i0 +𝜆∗ a,i0+1=1 and 𝜆∗ a,i=0∀ i ∈[ ka ]⧵{ i0 , i0 +1} . Define Then, (x∗ a,𝜋∗ v,ya,ca)a∈ A ,v∈ V ,i∈[ka] obviously fulfils the variable bounds(11e–g) and Eq.(11c). From (10b) and the definition of y follows (11b). We now show that (11d) holds: From our definition of ya , we obtain Since the objective is to minimize and the slopes sa,i in (11d) are decreasing for increasing i, we have ca=sa,i0 ya+ta,i0=maxi∈[ka−1]sa,i ya+ta,i which implies ca≥sa,i ya+ta,i for all a∈A and for all i∈[ka−1] , and hence, the solution of ModelD satisfies (11d). Therefore, all optimal solutions of ModelD are feasible for Model E, i.e. ( x ∗ a,𝜋 ∗ v )a∈ A ,v∈ V ,i∈[ka] ∈X E . The optimal objective function value of ModelD is ∑ a∈ALa( 𝜆∗ a,i0 ca,i 0 + 𝜆∗ a,i0+1 ca,i 0 +1 ) and equals ∑a∈A L a  c a by construction. (ii) Let ( x ∗ a ,𝜋 ∗ v ,y ∗ a ,c ∗ a)a∈A,v∈V be an optimal solution of Model E, that is, ( x ∗ a ,𝜋 ∗ v)a∈A,v∈V ∈X E . Consider i0∈[ka−1] such that D−𝛽 a,i0+1 ≤y∗ a≤D −𝛽 a,i0 . Then, there exist unique 0≤𝜆 a,i 0,𝜆 a,i 0+ 1 ≤1 with 𝜆 a,i 0 +𝜆 a,i 0+ 1 =1 such that  y a∶= 𝜆∗ a,i 0 D −𝛽 a,i 0 +𝜆∗ a,i 0 +1D −𝛽 a,i 0 +1and ca∶= 𝜆∗ a,i 0 ca,i0 +𝜆∗ a,i 0 +1ca,i0+1 . (12) ya=𝜆a,i 0 (D −𝛽 a,i 0 −D −𝛽 a,i 0 +1)+D −𝛽 a,i 0 + 1 (13) ⇔ 𝜆a,i0(ca,i0−ca,i0+1)= c a,i0 −c a,i0+1 D−𝛽 a,i0 −D−𝛽 a,i0+1 (ya−D−𝛽 a,i0+1 ) (14) ⇔ 𝜆a,i0ca,i0+𝜆a,i0+1ca,i0+1= c a,i0 −c a,i0+1 D−𝛽 a,i 0 −D−𝛽 a,i 0 +1 (ya−D−𝛽 a,i0+1)+ca,i0+ 1 (15) ⇔c a =s a,i 0 y a +t a,i 0. 198 R.Lenz, K.H.Becker 1 3 y ∗ a=𝜆a,i0D −𝛽 a,i 0 +𝜆a,i0+1D −𝛽 a,i 0 + 1 for all a∈A . We set 𝜆a,i∶= 0 for all i∈[ka]⧵{i0,i0+1} ; then, (x∗ a,𝜋∗ v,𝜆a)a∈ A ,v∈ V ,i∈[ka] fulfils (10b–g). Hence, ( x∗ a,𝜋∗ v )a∈ A ,v∈ V ,i∈[ka] ∈X D . Therefore, all optimal solutions of ModelE are feasible for ModelD, i.e. ( x ∗ a,𝜋 ∗ v )a∈ A ,v∈ V ,i∈[ka] ∈X D . From (11d), we have c∗ a=sa,i0y∗ a+ta,i0 . However, rearranging y ∗ a =𝜆a,i 0 D −𝛽 a,i0 +𝜆a,i 0 +1D −𝛽 a,i0+1 as in (12) to (14) yields 𝜆a,i0ca,i0 +𝜆a,i0+1ca,i0+1=sa,i0y∗ a+ta,i0 ; hence, ∑a∈A L a ( 𝜆a,i0 c a,i0 + 𝜆a,i0+1 c a,i0+1 )= ∑a∈A L a c ∗ a . ◻ 4.4 Comparison ofrelaxations In this section, we consider the continuous relaxations of the discrete models. As solvers relax the combinatorial part during the solution procedure, the tightness of the continuous relaxation plays a major role for the performance of the models. To this end, we use the standard projection to map their feasible regions onto the space of their common variables, i.e. onto the ( (𝜋 v ,x a ) v∈V,a∈A) -space: In the following, we denote the feasible regions of the continuous relaxations of the discrete ModelsA, B and C by XA, rel , XB, rel and XC, rel , respectively (Fig.2). XA,rel =Proj 𝜋v,𝜆a,i {( 𝜋v,𝜆a,i,xa ) | ( 𝜋v,𝜆a,i,xa ) satisfy Eqs. (7b)−(7f) and 𝜆a,i∈[0, 1]∀a∈A,∀i∈[ka]}, XB,rel =Proj 𝜋v,𝜆a,i{(𝜋v,𝜆a,i,x+ a,i,x− a,i,Δ+ a,i,Δ− a,i,za)|(𝜋v,𝜆a,i,x+ a,i,x− a,i,Δ+ a,i,Δ− a,i,za) satisfy Eqs. (8b)−(8h),(8j)and za,𝜆a,i∈[0, 1]∀a∈A,∀i∈[ka]} , XC,rel =Proj 𝜋v,𝜆a,i{(𝜋v,𝜆a,i,xa,Δa,i)|(𝜋v,𝜆a,i,xa,Δa,i)satisfy Eqs. (9b)−(9i) and 𝜆 a,i ∈[0, 1]∀a∈A,∀i∈[k a ] } . Fig. 2 Relation of feasible regions. The highlighted relations in grey are proved as propositions. The other relations X A ⊆X A ,rel , X C ⊆XC, rel and X B ⊆XB, rel are the canonical continuous relaxations of the corresponding discrete models. Proposition2 also implies X C ,X B ⊆X A ,rel . Since X D =X A ,rel , the splitpipe model can be seen as a relaxation of ModelsA, B and C 199 1 3 Optimization ofcapacity expansion inpotential‑driven… Proposition 4 Let 𝛼≠1 . Then, the following relationships hold for the continuous relaxation of ModelsA, B and C: 1. XA,rel ⊆XC, rel but XC, rel ≠XA,rel , 2. XA,rel ⊈ X B, rel and XB, rel ⊈ X A,rel , 3. XC,rel ⊈ X B, rel and XB, rel ⊈ X C,rel . Proof “ X A,rel ⊆XC, rel ”: Let ( x a ,𝜋 v ,  𝜆 a,i ) a∈ A ,v∈ V ,i∈[ka] be a solution to the continuous relaxation of ModelA, i.e. ( 𝜋 v ,  𝜆 a,i ) a∈ A ,v∈ V ,i∈[ka] ∈X A,rel . We show that it can be transformed to a point in XC, rel . First, ( x a ,𝜋 v ,  𝜆 a,i ) a∈ A ,v∈ V ,i∈[ka] satisfy Eqs. (9e–g) and (9i). By defining the potential loss constraint (9c) and the bounds of Δa,i (9h) are satisfied. With M(1) a,i = Φ(Da,i,x a) , M(2) a,i = Φ(Da,i,xa ) , the big M-formulation (9d) holds by construction. Equation(9b) also holds, since hence, ( x a ,𝜋 v ,  𝜆 a,i ,Δ a,i ) a∈ A ,v∈ V ,i∈[ka] is a solution of the continuous relaxation of ModelC, i.e. XA,rel ⊆XC, rel . The remaining relations between the continuous relaxations of the discrete models are shown by counter examples. Throughout we consider a single pipe a=(v,w) with La R a =1, D 𝛽 a,0 =1, D 𝛽 a,1 = 2 for given 𝛼>0 and 𝛽>0 . “XC, rel ⊈ XA,rel and XB, rel ⊈ XA,rel”: Let bv= 𝛼 √100 , bw=− 𝛼 √100 and 𝜋v,𝜋w∈[0, 3600] be given. Then, (𝜆a,0,𝜆a,1 , 𝜋v,𝜋w)∶=(0.5, 0.5, 3600, 3500)∈ X B,rel, X C,rel , where the the remaining variables of ModelB are given by e.g. x + a,0 =𝛼 √ 100, x+ a,1 =0, x− a,0 =0, x− a,1 =0, za= 1 , Δ+ a,0 =100, Δ+ a,1 =0, Δ− a,0 =0, Δ− a,1 =0 , and of Model C by xa = 𝛼 √100 and Δa,0 =100, Δa,1 =0 . However, (𝜆a,0,𝜆a,1,𝜋v,𝜋w)∶=(0.5, 0.5, 3600, 3500)∉ X D , where the remaining variable x a = 𝛼 √100 is uniquely determined by Eq.(7d), which contradicts Eq.(7b), i.e. 𝜋 v−𝜋w=100 ≠0.75 � 𝛼 √ 100 �𝛼 = ∑ ka i=0(𝜆a,i∕D𝛽 a,i)sgn (xa) � xa �𝛼 . “XA,rel ⊈ XB, rel and XC, rel ⊈ X B, rel ”: Now slightly change the above example to bv=1000 =−bw . Moreover, we set xa =1000 . Then (𝜆a,0,𝜆a,1, 𝜋v,𝜋w)∶=(0.5, 0.5, 3∕4×1000𝛼,0)∈ X A,rel, X C,rel , (16) Δ a,i∶= 𝜆a,isgn (xa)|xa| 𝛼 LaRaD −𝛽 a,i∀a∈A,∀i∈[ka] ⇒∑ka i=0 Δa,i (16) =sgn (xa) | xa | 𝛼∑ka i=0 𝜆a,iLaRaD−𝛽 a , i (10b) =𝜋v−𝜋w∀a∈A , (16) ⇒𝜆a,isgn (xa)|xa|𝛼= D 𝛽 a,i LaRa Δa,i∀a∈A,∀i∈[ka], ⇒ ka ∑ i=0 𝜆a,isgn (xa) | xa | 𝛼= ka ∑ i=0 D𝛽 a,i L a R a Δa,i (9e) ⇒sgn (xa) | xa | 𝛼= ka ∑ i=0 D𝛽 a,i L a R a Δa,i∀a∈A , 200 R.Lenz, K.H.Becker 1 3 where the remaining variables of the continuous relaxation of Model A are given by xa =1000 and of the continuous relaxation of Model C by Δa,0 =1∕2×1000𝛼 , Δa,1 =1∕4×1000𝛼 and xa =1000 . However, (𝜆a,0,𝜆a,1,𝜋v,𝜋w)=(0.5, 0.5, 0.75 ×1000𝛼,0)∉XB, rel , because it follows from Eqs. (8e–f) that x+ a,0 =x + a,1 =500 =0.5 x a and x − a,0 =x − a,1 =0 , and from (8c, g) it follows that za =0.5 and Δ+ a,0 =500 𝛼 and Δ+ a,1 =0.5 ×500 𝛼 and Δ− a,0,Δ− a,1 =0 . With these values, (8b) implies ∑ i ∈[ k a] (Δ + a,i −Δ − a,i )=3∕2×500 𝛼 ≠3∕4×1000 𝛼 =𝜋v−𝜋wfor 𝛼≠ 1. “ X B, rel ⊈XC, rel ”: We now choose bv= 𝛼 √75 =−bw and set x a >0 . Then (𝜆a,0,𝜆a,1,𝜋v,𝜋w, x+ a,0 ,x + a,1 ,x − a,0 ,x − a,1 ,z a, Δ + a,0 ,Δ+ a,1 ,Δ− a,0 ,Δ− a,1 )∶=(0.75, 0.25, 3600, 3525, 𝛼 √75, 0, 0, 0, 1, 75, 0, 0, 0)∈XB, rel . However, (𝜆a,0,𝜆a,1,𝜋v,𝜋w)=(0.75, 0.25, 3600, 3525)∉XC, rel , since xa = 𝛼 √75 follows uniquely from (9f), which together with Eq.(9b) leads to Δa,0 +2Δa,1 =75 . Moreover, from Eq. (9d) it follows that Δa,0,Δa,1 >0 , since 𝜆a,0,𝜆a,1 >0 and M(1) a,i = Φ(D a,i ,x a )> 0 holds by virtue of x a >0 . But then, Eqs. (9b) and (9c) cannot be satisfied simultaneously, because (9c) yields Δa,0 +Δ a,1 =75 . 5 Multi‑scenario modelling Expansion planning is a long-lasting process that involves high investment costs. Hence, practitioners seek network extensions that enable feasible operations even for different possible future demand scenarios. In practice, this process is typically approached as follows: Firstly operators typically identify multiple “bottleneck scenarios” that stress the network. Afterwards, they run simulations in order to decide on network expansions that ideally resolve all bottlenecks at once. For the problem of finding optimal gas network extensions for multiple demand scenarios, a model formulation and an algorithmic solution approach based on scenario decomposition has already been developed in the literature, see e.g. Schweiger (2016); Schweiger and Liers (2018). While these authors consider general new pipelines as extension candidates, we focus on loops in this paper. More precisely, we extend the formulations of the ModelsA, D and E from Sects.3 and 4 to the case of multiple scenarios, which we simply denote as Models A-MS, D-MS and E-MS. These three models have been chosen on the basis of our computational results in the single-scenario case (see Sects.6.2–6.4). It is well known that network expansions might worsen the flow situation and might allow for less network throughput or might render previously feasible instances as infeasible. This phenomenon is known as Braess’s paradox (Braess 1968) and has been shown to appear in nonlinear potential-driven networks, too, see Calvert and Keady (1993). To prevent the Braess phenomenon, most authors in the literature equip pipeline extensions with valves to allow for the flow and pressure patterns that were possible prior to the extension. Similarly, our modelling approaches in the following are equipped with binary variables that enable the model to switch off loops in individual scenarios. 201 1 3 Optimization ofcapacity expansion inpotential‑driven… 5.1 Discrete looping withlength variables (model A‑MS) The expansion planning problem for multiple demand scenarios 𝜔∈Ω concerns the cost-minimal selection of pipes to be looped such that the corresponding network enables feasible operation for each single scenario involved. For this purpose, we present a two-stage stochastic program, where the first-stage variables 𝜆a,i denote whether extensions are built. All remaining variables of the original model A have become scenariodependent recourse variables. As we would like to allow the model to switch off loops in individual scenarios, we modify the potential loss function to 𝜋 v,𝜔−𝜋w,𝜔=LaRa � 𝜆a,𝜔 ∑ ka i=0𝜆a,iD−𝛽 a,i+(1−𝜆a,𝜔)D−𝛽 a,0 � sgn (xa,𝜔) � xa,𝜔 �𝛼 , where the newly introduced binary variable 𝜆a,𝜔 indicates whether a loop extension of pipe a∈A is used in scenario 𝜔∈Ω . This happens if 𝜆a,𝜔=1 and 𝜆a,i>0 for some i∈[ka]⧵{0} is satisfied. In the case that an extension is not used, i.e. 𝜆a,𝜔=0 , the potential loss function reduces to its original form as given in Eq.(7b). Moreover, to avoid the resulting nonlinear product of binary variables, we introduce a continuous variable 𝛾 a,𝜔=𝜆a,𝜔 ∑k a i=0 𝜆a,iD −𝛽 a,i and linearize the product by virtue of Constraints(17–f). (17a) min 𝜆,x,𝜋,𝛾∑ a∈A La k a ∑ i=0 𝜆a,ica,i , (17b) s.t. 𝜋 v,𝜔 −𝜋 w,𝜔 = L aRa(𝛾a,𝜔+1−𝜆a,𝜔 D𝛽 a,0 )sgn (xa,𝜔) | xa,𝜔 | 𝛼∀a∈A,(v,w)=r(a),∀𝜔 ∈Ω, (17c) k a ∑ i=0 𝜆a,i=1∀a∈A , (17d) 𝛾 a,𝜔≤D −𝛽 a,0 𝜆a,𝜔∀a∈A,∀𝜔 ∈Ω, (17e) 𝛾 a,𝜔≤ ∑ i∈[k a ] 𝜆 a,i D𝛽 a,i ∀a∈A,∀𝜔 ∈Ω, (17f) ∑ i ∈[k a ] 𝜆 a,i D𝛽 a,i ≤ 1−𝜆 a,𝜔 D𝛽 a,0 +𝛾a,𝜔∀a∈A,∀𝜔 ∈Ω, (17g) ∑ a∈𝛿 + (v) xa,𝜔− ∑ a∈𝛿 − (v) xa,𝜔=bv,𝜔∀v∈V,∀𝜔 ∈Ω, (17h) 𝜋 v ≤ 𝜋v,𝜔 ≤ 𝜋v∀ v ∈V , ∀𝜔∈Ω, 202 R.Lenz, K.H.Becker 1 3 5.2 Split‑pipe looping withlength variables (model D‑MS) As it is the case with ModelsA andD, the following split-pipe looping modelD-MS is identical to the discrete ModelA-MS with the only difference being that the variables 𝜆a,i∈{0, 1} are relaxed to continuous variables 𝜆a,i∈[0, 1] . 5.3 Split‑pipe looping withefficient‑frontier constraints (model E‑MS) In our extension of ModelE to a two-stage stochastic program, the first-stage variables are given by the ya -variables, whereas the remaining variables are scenario-dependent second-stage variables. To indicate whether a loop extension of pipe a∈A is used in scenario 𝜔∈Ω and to represent the corresponding impact of the loop extension on the pressure loss, we introduce a binary variable 𝜆a,𝜔 and a continuous variable ya,𝜔 . Then, Constraints(19j–m) guarantee that ya,𝜔 =ya holds if the loop extension of pipea is used, i.e. if 𝜆a,𝜔=1 . Accordingly, if the loop extension is not used, i.e. 𝜆a,𝜔=0 , then ya,𝜔 =ya is satisfied, which guarantees that the pressure loss Constraint(19b) reduces to the original form as given in(11b). The model then reads as follows: (17i) x a ≤ xa,𝜔 ≤ xa∀a∈A,∀𝜔∈Ω, (17j) 𝛾a,𝜔∈ ℝ ≥0∀a∈A,∀𝜔∈Ω, (17k) 𝜆a,𝜔∈{0, 1}∀a∈A,∀𝜔∈Ω, (17l) 𝜆a,i∈{ 0, 1 }∀ a ∈A , ∀ i ∈[ k a]. (18a) min 𝜆,x,𝜋,𝛾∑ a∈A La k a ∑ i=0 𝜆a,ica,i , (18b–18k) s.t. (17b)−(17k) (18l) 𝜆a,i∈[0, 1]∀a∈A,∀i∈[ka] (19a) min y,c,x,𝜋,𝜆 ∑a∈A Lac a (19b) s.t. 𝜋v,𝜔−𝜋w,𝜔=LaRaya,𝜔sgn (xa,𝜔)|xa,𝜔|𝛼∀a∈A,(v,w)=r(a),∀𝜔∈Ω, (19c) ∑ a∈𝛿 + (v) xa,𝜔− ∑ a∈𝛿 − (v) xa,𝜔=bv,𝜔∀v∈V,∀𝜔 ∈Ω, 209 1 3 Optimization ofcapacity expansion inpotential‑driven… section All instances, we provide a summary of the results of all 500 instances and indicate the number of solved instances ( ∑ # solved), which includes the optimal solved instances (# opt) and the instances that were detected as infeasible (# infeas). Moreover, we state the number of instances, where at least one feasible solution was found within the time limit (# sol found), no matter whether optimality was shown for this solution. For all remaining instances, no feasible solution was found within the time limit nor was infeasibility proved. (2) The section All opt reports data for all instances for which all models under comparison have found an optimal solution. (3) In the section Only opt by, we compare the models with respect to the instances that they alone were able to solve to global optimality. In all three cases, we report computational time by means of the shifted geometric mean (sgm) and the arithmetic mean (amm). For instances that are solved to optimality by all models, we additionally present the sgm of the number of Branchand-Bound nodes. 6.2 Comparison ofdiscrete models (Models A, B andC) We begin our computational experiments with the Belgian gas network, which has frequently been used as a tool in research about network optimization (see e.g. Table 4 Comparison of discrete models New York network Test set Demand 100 Demand 200 Demand 500 Models A B C A B C A B C Allinstances # opt 0 142 0 9 412 58 0 0 0 # infeas 0 0 0 36 79 73 500 500 500 ∑ # solved 0 142 0 45 491 131 500 500 500 # sol found 500 500 500 384 421 348 0 0 0 Time (amm) 14400.0 12781.4 14400.0 13227.5 1768.0 11639.1 0.0 0.3 0.0 Time (sgm) 14400.0 12183.0 14400.0 8759.6 1000.0 7465.9 0.0 0.3 0.0 All opt # opt 0 8 0 Time (amm) – – – 2917.1 734.3 1470.1 – – – Time (sgm) – – – 563.3 537.6 881.1 – – – Nodes (sgm) – – – 126220 8918 109148 – – – Only opt by # opt 0 142 0 0 355 0 0 0 0 Time (amm) 14400.0 8903.2 14400.0 14400.0 2190.2 14400.0 – – – Time (sgm) 14400.0 8185.9 14400.0 14400.0 1570.5 14400.0 – – – 210 R.Lenz, K.H.Becker 1 3 Table 5 Comparison of discrete models Increasing the circuit rank of the Belgian network Test set + 0 Arcs + 2 Arcs + 4 Arcs + 6 Arcs + 8 Arcs + 10 Arcs Models A B C A B C A B C A B C A B C A B C All instances # opt 499 500 493 393 416 227 251 364 127 307 297 214 326 262 245 372 248 299 # infeas 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ∑ # solved 499 500 493 393 416 227 251 364 127 307 297 214 326 262 245 372 248 299 # sol found 500 500 500 500 500 500 498 479 500 500 381 500 499 357 500 500 316 500 Time (amm) 80.4 69.7 446.9 4688.8 4726.2 9094.5 8489.5 7306.5 11352.2 7018.2 9003.2 9463.2 6719.1 9817.5 8613.4 5123.6 9969.2 7151.8 Time (sgm) 9.1 39.1 20.3 1368.6 2069.0 4766.3 3492.2 4329.8 6399.7 2327.2 5812.4 4357.7 2489.2 7343.8 3805.4 1443.4 6824.2 2262.1 All opt # opt 493 221 117 183 192 207 Time (amm) 29.8 64.4 248.7 545.8 859.6 2632.2 753.0 1440.6 2148.4 809.8 3632.4 2276.1 1381.0 4934.4 2254.0 969.5 4870.5 1899.4 Time (sgm) 7.8 37.6 17.7 249.5 585.5 1213.8 258.0 787.3 495.3 278.2 1879.5 663.9 489.2 3291.6 814.1 326.4 2712.0 510.2 Nodes (sgm) 294 305 564 33817 8509 146470 20679 7216 25610 18061 16784 33299 30406 27272 37800 19120 20619 15523 Only opt by # opt 0 1 0 22 39 0 24 129 0 48 36 1 62 25 6 67 3 8 Time (amm) 14400.0 363.7 14400.0 11918.7 9747.6 14400.0 13283.2 8441.5 14400.0 9907.2 12030.4 14358.8 8811.9 12824.1 13859.5 5925.5 14154.0 13321.2 Time (sgm) 14400.0 363.7 14400.0 10717.3 8217.9 14400.0 12488.9 7165.3 14400.0 8101.0 11101.1 14352.9 6849.8 12245.5 12989.9 3866.6 13979.7 11744.4 211 1 3 Optimization ofcapacity expansion inpotential‑driven… De Wolf and Smeers (1996), DeWolf and Smeers (2000) and Babonneau et al. (2012)). Our computational results for the network are found in Table2. Demands are given in 106m3∕day and computational time in seconds. We first consider the performance of our models as depending on the demand. For the lowest demand B=100 about all instances are solved to optimality within a couple of minutes. With increasing demand, computational time increases and the number of instances solved to optimality decreases. However, at a demand of B=1000 , we have reached a point where an optimum has been found for less than a quarter of theinstances and more of 60 % of the instances are detected as infeasible. As a consequence of the high number of infeasibilities, the number of solved instances has increased and computational time decreased. This pattern suggests that the range of demand used here is sufficient to obtain an overview of the general behaviour of the models and that it is not necessary to test the models for a demand of less than B=100 and greater than B=1000 . Due to this typical pattern, which was frequently observed also during our preliminary computational experiments, we will generally use a test range from a low demand to a high demand with a large number of infeasibilities. We observe that, with respect to the sgm, Model A is the fastest model throughout, up to a factor 4.3 faster than Model B for a demand of B=100 , while B is the fastest with respect to the amm. This is due to the fact that Model B has a high number of instances that were solved fast (see e.g. those instances that were only solved by B). Model B also solves the largest number of instances. We note though that here Model B does not stand out with respect to the infeasible instances discovered and instances for which a feasible solution has been found. Let us remark in passing that during our preliminary tests of the two models mentioned in Sect.3.4, which were carried out on the rather simple Belgian network, the model in Humpola (2014) ran into the time limit at a demand of B=200 with about 2/3 of all instances, while the model in Borraz-Sánchez etal. (2016) exceeded the time limit in all cases for the same demand. (In contrast with this, the three models considered here ran into the time limit only with 0.8% to 6.6% of the same instances.) We now turn our attention to GasLib-40 (see Table3). Again we can see that Model A is the fastest model, with B only having an advantage at a demand of B=1000 with respect to the amm. This time Model A also tends to be better with respect to the number of instances solved. Finally note, that a demand of B=500 is too difficult for all three models, i.e. for unfavourable demand situations we have reached the limits of what is computational tractable. (But obviously, this does not impply that expansion problems on more complex networks cannot be solved for favourable demand situations.) The famous New York network is of particular interest for the capacity expansion problem as it is the only network in the water library where the arcs are already laid out with given diameters. Introduced by Schaake and Lai (1969), this instances has been used extensively (see e.g. Quindry etal. (1981), Bhave (1985), Bragalli etal. (2012)). On this network, Model B clearly performs best, both regarding the number of instances solved and the runtime (Table4). We note that we have instances here where Model C performs better than Model A in terms of both solved instances and runtime ( B=200 m3∕s ). 212 R.Lenz, K.H.Becker 1 3 It is well known that optimization problems on potential-driven networks get more difficult to solve the more cycles they contain because the existence of cycles leads to more complex patterns of flow directions (cf. the literature review in Shiono and Suzuki (2016)). To systematically test our models also with respect to this type of difficulty, we successively increase the circuit rank by adding up to ten new pipelines to the Belgian network, which results in five new networks Belgium + (2, 4, 6, 8, 10) arcs, see Fig.3. Here, the intention was to add arcs that connect different regions of the network and have a high impact on the network topology. For the new pipelines we use in our data set an original diameter Da,0 that is equal to the average diameter of the existing pipelines in the network and use moderate demand scenarios with B=200 . We observe (see Table5) that for increasing circuit rank the number of instances solved to optimality at first decreases and then increases again for Models A and C, while it continues to decrease for Model B. A similar pattern applies with respect to computational time, where the runtime at first increases for Models A and C and then decreases for the instances with 6 or more additional arcs, while the runtime of Model B increases nearly throughout. This leads to a situation where Model A turns out to have the best sgm runtime for all circuit ranks, B solves the most instances for smaller circuit ranks, and A solves the most instances for larger circuit ranks. Finally, while Model C shows poor results for small circuit ranks, it gains ground Table 6 Comparison of split-pipe models Belgian network with different demand scenarios Test set Demand 100 Demand 200 Demand 500 Demand 1000 Models D E D E D E D E All instances # opt 500 500 494 496 318 315 77 82 # infeas 0 0 0 0 5 5 310 310 ∑ # solved 500 500 494 496 323 320 387 392 # sol found 500 500 500 500 495 495 189 189 Time (amm) 1.0 0.7 256.2 169.9 5850.9 5891.4 3655.4 3533.2 Time (sgm) 1.0 0.6 7.0 4.8 676.1 640.9 80.3 79.1 All opt # opt 500 492 269 60 Time (amm) 1.0 0.7 84.0 42.0 841.4 951.2 1903.4 1657.4 Time (sgm) 1.0 0.6 5.5 3.8 80.0 78.3 205.0 196.8 Nodes (sgm) 24 21 240 197 16381 15153 50809 53529 Only opt by # opt 0 0 2 4 49 46 17 22 Time (amm) – – 9662.2 5916.3 8554.7 8457.3.1 10335.6 9148.2 Time (sgm) – – 3106.3 415.1 4037.7 3234.6 6748.7 6067.2 213 1 3 Optimization ofcapacity expansion inpotential‑driven… Table 7 Comparison of the split-pipe models GasLib-40 and New York networks GasLib-40 New York Test set Demand 50 Demand 100 Demand 500 Demand 1000 Demand 100 Demand 200 Demand 500 Models D E D E D E D E D E D E D E All instances # opt 499 500 371 413 237 244 206 243 500 500 421 421 0 0 # infeas 0 0 0 0 0 0 160 160 0 0 79 79 500 500 ∑ # solved 499 500 371 413 237 244 366 403 500 500 500 500 500 500 # sol found 500 500 500 500 500 500 340 340 500 500 421 421 0 0 Time (amm) 32.4 3.2 4045.7 3071.0 9273.6 9002.8 6154.6 5477.9 3.1 2.9 1.9 2.9 0.0 0.0 Time (sgm) 3.5 1.4 215.7 164.3 4832.4 4452.5 756.7 687.8 3.0 2.8 1.8 2.8 0.0 0.0 All opt # opt 499 354 214 186 490 417 0 Time (amm) 3.6 3.2 375.3 341.7 3291.6 2738.2 5434.6 4247.7 3.1 2.9 1.3 2.4 – – Time (sgm) 3.3 1.4 40.4 37.1 1300.6 1050.9 2842.7 2226.8 3.0 2.8 1.3 2.4 – – Nodes (sgm) 42 16 4270 4040 157062 33984 338695 285031 698 677 182 176 – – Only opt by # opt 0 1 17 59 23 30 20 57 0 0 0 0 0 0 Time (amm) 14400.0 4.5 11605.1 5349.1 10890.4 10570.0 12434.0 10905.4 – – – – – – Time (sgm) 14400.0 4.5 5266.3 1317.9 7837.6 8512.7 11240.0 10117.4 – – - – – – 214 R.Lenz, K.H.Becker 1 3 Table 8 Comparison of split-pipe models Increasing the circuit rank of the Belgian network Test set +0 Arcs +2 Arcs +4 Arcs +6 Arcs +8 Arcs +10 Arcs Models D E D E D E D E D E D E All instances # opt 500 500 155 235 58 136 26 170 12 194 30 308 # infeas 0 0 0 0 0 0 0 0 0 0 0 0 ∑ # solved 500 500 155 235 58 136 26 170 12 194 30 308 # sol found 500 500 500 500 500 500 499 500 500 500 500 500 Time (amm) 1.0 0.7 10188.8 7963.9 12928.2 10819.0 13766.4 9781.5 14087.2 9346.4 13634.8 6224.3 Time (sgm) 1.0 0.6 3192.4 1173.1 9611.2 3952.4 11843.6 2790.4 12944.1 2675.5 11243.6 844.2 All opt # opt 500 101 35 22 12 29 Time (amm) 1.0 0.7 761.7 477.2 1087.5 1015.0 2047.6 246.8 1364.9 172.1 1275.8 23.4 Time (sgm) 1.0 0.6 91.2 34.4 258.5 37.2 237.8 14.5 160.3 32.0 194.0 10.9 Nodes (sgm) 24 21 21956 4642 71439 4495 29888 1023 19818 2046 17299 582 Only opt by # opt 0 0 54 134 23 101 4 148 0 182 1 279 Time (amm) – – 10526.8 4762.4 12222.9 3738.5 14103.6 1256.1 14400.0 1454.7 14392.9 1289.6 Time (sgm) – – 3777.5 407.2 8689.6 387.3 13637.5 155.7 14400.0 199.6 14392.4 172.6 215 1 3 Optimization ofcapacity expansion inpotential‑driven… with increasing circuit ranks and performs second best for the instances with Belgium + 10 arcs. In sum, we can state that overall ModelsA and B perform better than ModelC, with Model A clearly outperforming Model B on problems with a higher circuit rank (including GasLib-40). The performance of Model C may partly be a consequence of the fact that Model A, as we have shown in Sect.4.4, has a tighter relaxation than Model C, even though there are a smaller number of entire data sets and a larger number of individual instances where ModelC outperforms both ModelA and ModelB, in particular instances with a higher circuit rank. ModelC is therefore a valuable alternative in cases where ModelsA and B require too much computational time to prove infeasibility or solve an instance to optimality. 6.3 Comparison ofsplit‑pipe models (models D andE) Again we start with the Belgian network (Table6) and proceed from easily solved instances to scenarios with a large number of infeasible instances. It turns out that both models solve a comparable number of instances, with the new ModelE being somewhat faster. It is worth noting that for a demand of B=500 , i.e. for more difficult instances, 10% of the instances are solved to optimality by only one of the two models, i.e. in this case both models complement each other well. To compare the performance on gas networks of larger sizes with a more complex cyclic structure, we turn to the GasLib-40 network (Table7 on the left). Here, Model E performs clearly better: it solves up to 11 % more instances (demand B=100, 1000 ) and requires less runtime (up to a factor of 2.5 in the case of B=50 ). To further test the performance of the split-pipe models, we consider the New York water network again (Table7 on the right). Here the results are balanced: both models solve the majority of instances in at most 3 s (sgm). Finally, we investigate the performance of the models for increasing circuit rank. We observe (Table8) that our novel Model E consistently outperforms D, with the ratio of the runtimes of the model becoming the more favourable for E the higher the circuit rank, up to a factor of 13.9 for 10 additional arcs. Similarly, Model E performs significantly better with respect to the number of solved instances, with the best ratios of solved instances reached at 16.2 for Belgium + 8 arcs and 10.0 for Belgium + 10 arcs. On the basis of our computational experiments for the split-pipe case, we can clearly recommend the use of our novel model instead of the model from the literature, particularly for networks with a more complex cyclic structure, even though one will, of course, encounter instances where the other model is more successful. 6.4 Comparison ofdiscrete andsplit‑pipe models Let us now compare the discrete with the split-pipe approaches. As mentioned in the introduction, this is the first time that such a comparison is carried out for networks of a practically relevant size and complexity. Of course, one may expect that the split-pipe model will perform better on average, but as our computational 216 R.Lenz, K.H.Becker 1 3 experiments show (compare the results for Belgium +4, +6, and +8 arcs in Tables5and8, for example), there are data sets where the discrete approach yields more optimal solutions within the time limit of 4 hours than the split-pipe approach. For this reason it is unclear what size the performance gap may actually take. In the following, we will look at two criteria: (a) potential cost savings and (b) computational time. Concerning cost savings, Table9 shows the average gain of the best known solutions of the split-pipe models over the discrete models. Here, we consider only those instances where an expansion actually takes place and where either at least one splitpipe and at least one discrete model have been solved to optimality (column Optimal) or where both approaches provide at least one solution (column Feasible). For these instances, we then consider the average gain of the best known split-pipe over the best-known discrete one. The column Feasible additionally shows the percentage of instances where the best-known solution is provided by the split-pipe models. As Table 9 demonstrates, the split-pipe approach yields cost savings on all data sets. While the results for the comparably simple networks of Belgium and New York are rather low, the benefit of realizing a split-pipe solution can be considerable for networks with a complex cyclical structure. Moreover, Fig.4 shows that for all feasible instances with non-zero objective function value, the best solutions were practically always found by the split-pipe models, with our novel Model E delivering the best results. Even more, further analysis of the data presented in Tables2, 3, 4, 5, 6, 7 and 8 reveals that the split-pipe models optimally solve or detect (in)feasibility in all cases except 1 (Model E) or 2 (Model D) instances, whereas the number of instances with unknown status is much higher for the discrete models, namely 318 (Model A), 578 (Model B) and 348 (Model C) out of the 8000 instances we have calculated per model. Therefore, in view of the much larger number of feasible solutions by the split-pipe models, the economic benefit of these models goes well beyond the cost savings and the percentage of instances depicted in Table9 and Fig.4. To compare the overall runtime performance of all discrete and split-pipe models, we use a performance profile(Dolan and Moré (2002)). It is based on the performance ratio, i.e. the runtime of a particular data set with the model under consideration divided by the best runtime for that data set with any of the five models. The performance profile describes on the y-axis the fraction of instances among all solved (i.e. optimal or infeasible) instances that the model could solve with a performance ratio of up to the corresponding number on the x-axis. Clearly, models are to be preferred when their profile shows higher y-values for fixed x-values and lower x-values for fixed y-values. To exclude trivial cases, we disregard instances that were not solved by any model and those that were solved by all models in less than 1 s CPU time. As we can see, the runtime of Model E dominates the runtime of all other models across the spectrum of performance ratios (Fig.5). Further insights about the model performances for different data sets are gathered from Fig.6, which depicts the sgm of the runtime of all instances for different network types and confirms the dominance of the split-pipe models. Again ModelE turns out to perform best. 217 1 3 Optimization ofcapacity expansion inpotential‑driven… 6.5 Outlook: comparison ofmulti‑scenario models To conclude our computational experiments and provide an outlook on potential future applications, we study the behaviour of ModelsA, D and E in the multiple-scenario case (ModelsA-MS, D-MS and E-MS). ModelsB and C have been left out here due to their comparably weak performance as seen in Sect.6.2, 6.4. In particular, we present the test results of the multiple-scenario models on the (original) Belgian network and on its extended versions with additional arcs (circuit rank). These networks have been chosen because we have seen in the previous sections that on these networks all three models were able to find feasible solutions for a similar number of instances and the networks are sufficiently small to suggest that the models may be able to deal with the additional complexity resulting from the multiple-scenario setting. When evaluating potential pipeline expansions, practitioners typically look at about 5–10 scenarios. For this reason, all instances considered here will consist of 10 scenarios each. Each instance was constructed by combining 10 of the demand scenarios previously used in the single-scenario case in Sects.6.2, 6.3. More precisely, for the Belgian network each of the 125 instances used here has been composed of two scenarios with a total demand of B=100 , three scenarios with B=200 , three scenarios with B=500 , and two scenarios with B=1000 . In the circuit rank experiment, the previously 500 demand scenarios with B=200 in Sects.6.2–6.3 were grouped into 50 multiple-scenario instances with 10 demand scenarios each. As preliminary computational tests showed that the Table 9 Gain of split-pipe over discrete problem Network Test set Optimal Feasible Gain (%) # instances Gain (%) # instances SP better on # instances [%] Belgium Demand 100 2.2 500 2.2 500 100.0 Belgium Demand: 200 1.6 498 1.6 500 100.0 Belgium Demand: 500 1.1 332 1.2 495 100.0 Belgium Demand: 1000 1.3 74 1.4 189 100.0 GasLib-40 Demand: 50 321.5 147 321.5 147 100.0 GasLib-40 Demand 100 2.1 326 3.7 500 100.0 GasLib-40 Demand: 500 – 0 10.7 498 100.0 GasLib-40 Demand: 1000 0.7 2 8.0 175 100.0 New York Demand: 100 1.0 142 1.1 500 100.0 New York Demand: 200 0.8 413 0.8 415 100.0 New York Demand: 500 – 0 – 0 – Belgium+2 arcs Cycle rank 1.6 268 1.5 500 99.8 Belgium+4 arcs Cycle rank 8.2 151 5.3 499 96.4 Belgium+6 arcs Cycle rank 14.1 169 12.1 498 98.8 Belgium+8 arcs Cycle rank 17.4 185 14.5 498 99.6 Belgium+10 arcs Cycle rank 71.0 284 50.7 490 100.0 218 R.Lenz, K.H.Becker 1 3 additional complexity of the multiple-scenario setting led to a rather small number of instances for which the models could find feasible solutions, we will provide in the following the results of computational experiments where the demand data from Sects. 6.2–6.3was scaled down by a factor of 0.4 for all instances. (Note though that this does not mean that our demand sets are just scaled versions of each other. The data sets used in the previous sections and here were generated independently for each given total network demand B .) Let us first look at the instances on the Belgian network. As shown in Table10 none of the models was able to solve any of the 125 instances to optimality or prove infeasibility. Clearly, the additional complexity of the multiple-scenario setting is taking its toll and we can only evaluate the models by looking at the instances for which feasible solutions have been found. Here ModelE-MS takes the lead by finding solutions for 119 out of 125 instances, with 68 of these being better than the solutions found by the other two models. ModelD-MS has found the best solution for another 60 instances, while the discrete ModelA-MS can only offer best solutions for 3 instances. A similar pattern can be observed in the case of the circuit rank experiment (Table11) where ModelE-MS was able to provide solutions for all but 2 instances of the 300 = 6 × 50 instances, while the next best ModelD did so for all but 22 instances. Moreover, again ModelE turns out to have found the largest number of best solutions among the three models, namely for 245 out of the 300 instances. This time, however, in contrast with the situation in the previous experiment (Table10), optimal solutions were found, and interestingly so mainly by ModelA. It should be noted though that there are only 3 out of 300 instances for which ModelA has found a better solution than the two split-pipe models, which again demonstrates the superiority of the split-pipe approach. Even though the split-pipe models cannot be solved to proven optimality, the authors of the present paper strongly assume on the basis of their experience with solving capacity expansion models that a considerable number of instances is actually solved to optimality by both D-MS and E-MS. This assumption rests on the observations that (a) we know from the discrete model in the case of several Fig. 4 Percentage of instances for which the best solution was found Belgium GasLib-40 New York Cycle Rank 0 20 40 60 80 100 No instancesin% Model A Model B Model C Model D Model E