Dynamic optimization of natural gas pipeline networks with demand and composition uncertainty

Journal Pre-proofs Dynamic optimization of natural gas pipeline networks with demand and composition uncertainty Kai Liu, Lorenz T. Biegler, Bingjian Zhang, Qinglin Chen PII: DOI: Reference:

S0009-2509(19)30939-X https://doi.org/10.1016/j.ces.2019.115449 CES 115449

To appear in:

Chemical Engineering Science

Received Date: Revised Date: Accepted Date:

24 December 2018 29 November 2019 18 December 2019

Please cite this article as: K. Liu, L.T. Biegler, B. Zhang, Q. Chen, Dynamic optimization of natural gas pipeline networks with demand and composition uncertainty, Chemical Engineering Science (2019), doi: https://doi.org/ 10.1016/j.ces.2019.115449

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier Ltd.

Dynamic optimization of natural gas pipeline networks with demand and composition uncertainty Kai Liua , Lorenz T. Bieglerb,∗, Bingjian Zhanga,∗, Qinglin Chena a

School of Materials Science and Engineering, Guangdong Engineering Center for Petrochemical Energy Conservation, Sun Yat-Sen University, No.135, Xingang West Road, Guangzhou, 510275, China b Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, USA

Abstract Pipelines are one of the most efficient methods to transport large quantities of natural gas from gas reserves to main markets. However, because of the unsteady nature of gas flow, the operation of pipelines is always dynamic. Furthermore, the uncertainties in demand and the fluctuation of gas composition also make the efficient operation of these dynamic processes challenging. This study addresses the problem of determining the optimal operation to minimize compression costs, while considering demand and gas composition uncertainties. A dynamic pipeline network model is developed with rigorous thermodynamic equations, allowing accurate calculation of gas compressibility factor at any temporal and spatial point. The supply gas composition and demand nodes flow rates are assumed to be uncertain. To deal with these uncertainties, a robust optimization algorithm is applied using back-off constraints calculated from Monte Carlo simulation. Through successive iterations, the algorithm terminates at a solution that is robust to a specified level of process variability with minimal cost. We show from the case studies that the formulated model and the algorithm can successfully address the problem with acceptable computational cost. Keywords: Dynamic optimization, Uncertainty, Pipeline network, Back-off constraints, Monte Carlo simulation

∗

Corresponding author Email addresses: [email protected] (Kai Liu), [email protected] (Lorenz T. Biegler), [email protected] (Bingjian Zhang), [email protected] (Qinglin Chen)

Preprint submitted to Chemical Engineering Science

November 28, 2019

1. Introduction As an energy source, natural gas is expected to grow substantially due to its sustainable supply and reduced adverse impacts on the climate and the environment [1]. The recent shale gas revolution makes natural gas even more attractive[2]. Gas processing plants, however, are often far away from the main markets, and thus pipeline networks are often used to transport large volumes of gas from processing plants to different receivers, including power generation plants and chemical plants. Pipeline networks typically consist of supply nodes, demand nodes, compressor stations, and pipelines with different diameters and lengths. Compressor stations are built to compensate for the pressure loss due to friction, which is a major operating cost of gas pipeline networks [3]. One of the key issues in the operation of pipeline networks is the presence of uncertainty arising from supply, demand, equipment failure and etc. The increasing demand for natural gas is accompanied by greater stresses on gas pipeline networks, and unless the uncertainties are anticipated and dealt with accordingly, deterioration of performance and production losses may occur. Optimization of pipeline networks has been extensively studied. These can be classified into two primary groups, steady-state pipeline problems and dynamic pipeline problems. For long-term planning problems, the goal is to design costefficient network configuration or obtain optimal operation strategy for years or decades, and thus is often assumed to be steady-state [4]. Because of the existence of discrete variables, such as flow directions and compressor station locations, long-term problems are typically formulated as nonconvex mixed-integer nonlinear programming models (MINLP) [5] [6], which are difficult to solve directly and are usually reformulated to make the model tractable. Geißler et al. used tailored graphbased model reformulation techniques that yield block-separable MINLP structures, which scale well with the size of the problem [7]. Mikolajkov´a et al. linearized the MINLP model of pipeline networks to an MILP model to achieve a good trade-off between model accuracy and computing efficiency [8]. On the other hand, the reformulation of MINLP models of gas pipelines as continuous NLP models also gains popularity, since NLP models can be solved much faster with standard NLP solvers. Koch et al. reformulated discrete aspects with complementarity constraints and a mathematical program with equilibrium constraints (MPEC) model was solved instead [9]. Rose et al. compared MINLP models with continuous reformulations and concluded that the continuous reformulation could be solved faster than the MINLP models and this tendency is more evident for larger instances [10]. Another popular topic for long-term pipeline network problems is multi-objective optimization, where typical conflicting objectives are the minimization of compressor power consumption and the maximization of transported gas volume [11] [12] [13]. For real life planning 2

tasks, the optimization model can be arbitrarily complicated, and therefore models with different fidelity are developed. Schmidt et al. proposed detailed stationary optimization models for gas networks, which considers highly accurate physical models for pipeline network components [14]. This model can be used either as a validation model for the final-stage design or a benchmark model to compare different designs. However, steady-state models cannot deal with transient behaviors of gas for short-term operations, such as daily operation of compressors, transitional stages, and process control, as the state of gas in pipelines changes frequently over time under these circumstances. Therefore, dynamic models of pipeline networks are used to handle it. Osiadacz proposed a partial differential equation (PDE) to simulate the isothermal transient gas flows in networks [15], which has been extensively used for transient pipeline network problems [16] [17] [18]. When formulated as an optimization problem, PDE-constrained models often require a lot of effort to solve. Baumrucker and Biegler developed an MPEC model for the operation of pipeline networks that allows flow reversal, and a multi-period approach was proposed to simplify the original PDE model [19]. Gopalakrishnan and Biegler proposed an economic Nonlinear Predictive Control (NMPC) formulation to bridge the gap between real-time optimization and advanced control [3]. Chiang and Zavala presented an optimization model that considers the spatio-temporal interactions between gas and electric transmission networks, improving the economic performance and process flexibility significantly by coordinated dispatch [20]. While most transient models of pipeline networks are formulated as NLPs, MINLP or Mixed-Integer PDE-constrained models can also be found in [21] [22]. Herr´an-Gonz´alez et al. presented two numerical schemes within MATLAB-Simulink library to model and simulate gas distribution pipeline networks, which shows that the inclination term has an impact on the pressure drop on pipelines [23]. A comprehensive review on the optimization of gas pipeline network can be found in [4]. Moreover, many approaches for pipeline network optimization apply the ideal gas law and assume a constant gas compressibility factor (Z-factor). This assumption is acceptable when the variation of pressure and gas composition is moderate. However, in transient or dynamic flows, especially under circumstances of start-up, shutdown, and flow-rate changes, the state of gas in pipeline networks may undergo significant changes that the constant Z-factor assumption may lead to inaccuracy. Behrooz and Boozarjomehry used Hall-Yarborough implicit equation and the Papay explicit equation in dynamic simulation and estimation for pipeline network model [24]; this does not consider the influence of gas composition on gas transportation. Guandalini et al. developed a dynamic pipeline model to investigate the impact of hydrogen injection on the transport of gas in pipelines, which is applied to 3

a single pipeline and limited to the scope of simulation instead of optimization [25]. In this study, we focus on transient processes of natural gas pipeline networks under uncertainty, and a PDE-constrained optimization model is developed to represent the dynamic pipeline networks. In addition, an empirical model is incorporated into the dynamic model to calculate the Z-factor. Although optimization under nominal condition has demonstrated its effectiveness in dealing with pipeline problems, the optimal solution obtained under nominal condition can be suboptimal and even infeasible in practice. Especially in dynamic processes, parametric uncertainty and process disturbances, combined with complex nonlinearity of the models, can lead to unexpected performance. The main feature that distinguishes parametric uncertainty from disturbances is that the former is often time-invariant or slow time-varying while disturbances can vary quickly with time. Nevertheless, the approaches and tools used for parametric uncertainty can be easily extended to processes with disturbances [26]. Therefore, only parametric uncertainty is taken into account in the following sections. Uncertainty in pipeline networks comes from several sources, including measurement noise, weather change, sales gas processing fluctuation, emergency events, and market uncertainty. Behrooz used an unscented transform as the uncertainty propagation tool to manage uncertainty in a stationary gas transmission network model [27]. This approach was extended to the application of dynamic optimization of gas pipeline networks with 9 nodes and two compressor stations [28]. In these studies, a derivative free optimization strategy is applied and optimization under demand uncertainty was completed within a CPU hour of computation. More generally, strategies used to solve optimization problems with uncertainty are classified as robust optimization and stochastic programming [29]. In stochastic programming, a subset of decisions is made in the first stage and the remaining decisions can be made in subsequent stages, whenever the uncertainty is realized. In robust optimization, the approach guarantees feasibility over a specified uncertainty. Because of the lack of recourse action, robust optimization often generates more conservative solutions in comparison with stochastic programming, but with far less computational cost. Therefore, one can argue that robust optimization is more appropriate for short-term planning problems, since the feasibility over a specified set of uncertain parameters is the major concern and recourse decisions are not compulsorily needed. In robust optimization, the most crucial step is to identify the worst-case scenario and optimize it [30]. A less expensive optimization approach is to use back-off terms to tighten the constraints, in order to guarantee constraints satisfaction in the presence of uncertainty to a certain level of confidence [26]. Diehl et al. proposed an approximate 4

robust formulation that linearizes the uncertainty set, and then used the dual norm to represent the back-off terms [31]. Such a linearization strategy, though accurate in linearized worst-case feasibility, may result in inaccurate estimation of back-off terms when the constraints are highly nonlinear and the range of uncertain parameters is large. Shi et al. used Monte Carlo Simulation to calculate the magnitudes of the back-off terms for the robust optimization of grade transitions in a polyethylene solution polymerization process [32], and demonstrated its performance when additional uncertainty parameters are considered. Koller et al. extended this strategy for simultaneous design, control and scheduling of multiproduct systems under uncertainty in a systematic way[33]. Maußner and Freund proposed a general algorithm for the optimization of maleric anhydride synthesis under uncertainty, in which cubature rules are used to speed up the the calculation of back-off terms[34]. Demand flow rates and gas composition are assumed to be uncertain. To simplify the problem, we divide the gas components into methane and the other gases, assuming the relative ratios of the components in the other gases remain in the same proportion. To deal with these uncertainties, a robust optimization algorithm is applied using back-off constraints, and Monte Carlo simulation is used to calculate the back-off terms. The advantages of the algorithm are the easy implementation and adjustable confidence level for the solution. Through successive iterations, the algorithm terminates when no constraint violations are found and the change of back-off terms is within tolerance. Both constant and time-dependent back-off terms are compared in this study. The two primary contributions of this study include: (1) development of a dynamic model of gas pipeline networks that takes into account non-ideal gas behavior and demand changes, which has been discretized into a standard NLP problem that can be solved in near real-time, and (2) application of a robust optimization strategy to deal with the dynamic natural gas pipeline network models with demand and gas composition uncertainty. The paper is organized as follows. In section 2, we develop a dynamic pipeline network model with an empirical equation of state. Section 3 explains the optimization strategy under uncertainty, followed by four case studies in section 4 to verify the model and algorithm, in which a large-scale network is used to check the scalability of the model. Finally, concluding remarks are made in section 5. 2. Model development In this section, a detailed PDE-constrained model is presented to describe natural gas transmission pipeline networks. These pipelines often extend thousands of miles and gas pressure decreases along the pipelines due to friction. Therefore, compressor 5

stations are built to compensate for the pressure loss. In this model, our goal is to minimize the compressor power consumption, and satisfy customer demand and operational safety requirement at the same time. 2.1. Pipeline model development Gas dynamics in a pipe is governed by a hypobolic systems of PDEs [35], We assume that a set L of pipes is included in the model. For a given pipe l ∈ L, the conservation of mass, momentum, and energy is described as [14]: ∂ρl (τ, z) ∂(ρl (τ, z)vl (τ, z)) + = 0, (1a) ∂τ ∂z ∂(ρl (τ, z)vl (τ, z)) ∂pl (τ, z) ∂(ρl (τ, z)vl2 (τ, z)) ∂hl (z) + + + gρl (τ, z) ∂τ ∂z ∂z ∂z |vl (τ, z)|vl (τ, z) + λl (τ, z) ρl (τ, z) = 0, (1b) 2Dl ∂Tl (τ, z) ∂Tl (τ, z) Tl (τ, z) ∂Zl (τ, z) ∂pl (τ, z) Al ρl (τ, z)cp ( + vl (τ, z) ) − Al (1 + ) ∂τ ∂z Zl (τ, z) ∂τ ∂τ Tl (τ, z) ∂Zl (τ, z)(τ, z) ∂pl (τ, z) ∂hl (z) − Al vl (τ, z) + Al (τ, z)ρl (τ, z)vl (τ, z)g Zl (τ, z) ∂Tl (τ, z) ∂z ∂z + πDl U (Tl (τ, z) − Tsoil ) = 0, (1c) where τ ∈ T := [0, T ] is the time dimension, T is the final time, z ∈ Xl := [0, Ll ] is the axial dimension, Ll is the length, ρl (τ, z) is the gas density, vl (τ, z) is the gas velocity, pl (τ, z) is the gas pressure, Zl (τ, z) is the gas compressibility factor, g is the gravitational acceleration, hl (z) is the pipe height, λl (τ, z) is the pipe friction coefficient, Dl is the pipe diameter, Al is the constant cross-sectional area, cp is the heat capacity at constant pressure, Tl (τ, z) is the gas temperature, U is the heat transfer coefficient, and Tsoil is the soil temperature. In the momentum conservation equation (1b), p is much greater than ρv 2 , therefore, the third term in Eq. (1b) is normally ignored without significant impact. The large-scale pipeline networks can extend to regions with different geological conditions, leading to pipes with elevation changes and temperature fluctuations. However, this requires detailed geological information on each segment of the pipelines. Besides, the energy conservative equation (1c) describing the temperature profiles is highly nonlinear PDE, which can be difficult to solve in the scope of dynamic optimization. Recently, simulation based optimization methods are proposed to include energy conservative equation to optimize the dynamic optimization of pipeline 6

networks [28], which can provide accurate gas temperature profiles. Moreover, since temperature changes only occur in very long pipelines or due to downstream pressure changes, such as compression units and lamination, the impact of temperature is commonly neglected during design and operation optimization [3] [19] [36]. Therefore, we assume an isothermal gas flow model with horizontal pipes, and the resulting pipeline model is simplified as follows: ∂ρl (τ, z) ∂fl (τ, z) + = 0, ∂τ Al ∂z ∂fl (τ, z) ∂Al pl (τ, z) −λl (τ, z) fl (τ, z)|fl (τ, z)| + = , ∂τ ∂z 2Dl Al ρl (τ, z)

(2a) (2b)

where fl (τ, z) is the mass flow rate. Depending on the Reynolds number Rel (τ, z), gases can be classified either laminar flow or turbulent flow, and different types of gases applies to different empirical equations for the calculation of friction factor. To simplify the model, we don’t consider gas reversal or peculiar conditions where gas flow rate approaches zero due to demand changes. Our experiments show that Rel (τ, k) is always present in the turbulent region under the operating conditions of interest; therefore, the Colebrook-White equation is applied for turbulent friction coefficient calculation, which are given as

Rel (τ, z) =

fl (τ, z)Dl , Al ν

(3a)

1 2.51 l p p = −2log10 ( + ), λl (τ, z) Rel (τ, z) λl (τ, z) 3.71Dl

(3b)

where ν is dynamic viscosity, and l is the pipe rugosity. Our numerical trials show 2.51 √ the Reynolds number to be very high (above 2000000), which makes Rel (τ,z)

λl (τ,z)

l only a small fraction of 3.71D . As a result, a simplified explicit friction coefficient l equation is used, which retains to dominant term in the Colebrook-White equation:

1 l p = −2log10 ( ). 3.71Dl λl (τ, z)

(4)

The gas law is given by ρl (τ, z) =

pl (τ, z)M , Zl (τ, z)RTgas

7

(5)

where M is molecular weight, R is gas constant and Tgas is the operating temperature. The compressibility factor Zl (τ, z) represents the deviation of real gas from ideal gas, the calculation of which is described in the following sections. PDEs of each pipe are discretized by using finite difference scheme in space with Nz points and equal length ∆z, and implicit Euler in time with Nt points and equal length ∆τ , which has proved to be effective in [37]. We define the sets X¯ := {1, ..., Nz } and T¯ := {1, ..., Nt }, and the discretized version of pipe Eq. (2) is given by: ρl,t+1,k − ρl,t,k fl,t+1,k+1 − fl,t+1,k + = 0, l ∈ L, t ∈ T¯ /{Nt }, k ∈ X¯ /{Nz } ∆τ Al ∆z fl,t+1,k − fl,t,k Al (pl,t+1,k+1 − pl,t+1,k ) −λl fl,t+1,k |fl,t+1,k | + = , ∆τ ∆z 2Dl Al ρl,t+1,k l ∈ L, t ∈ T¯ /{Nt }, k ∈ X¯ /{Nz } pl,t,k M ρl,t,k = , l ∈ L, t ∈ T¯ , k ∈ X¯ Zl,t,k RTgas 1 l p = −2log10 ( ). 3.71Dl λl,t,k

(6a)

(6b)

(6c) (6d)

The junctions in the real networks are represented by a set N of nodes in the model. We assume that there is no loss of mass and energy at the nodes and the physical junctions have no appreciable volume. Therefore, for each node n ∈ N , the sum of flows into the node equals the sum of flows out of them as shown in Fig. 1, which is given by: X X supply X X customer fl,t,Nz + fn,t = fl,t,1 + fn,t , t ∈ T¯ (7) l∈Lin n

n∈Ns

l∈Lout n

n∈Nd

out where Lin n represents the set of pipes that send gas to node n, Ln represents the supply set of pipes that receive gas from node n, Ns is the set of supply nodes, fn,t customer represents flows from supply nodes, Nd is the set of demand nodes and fn,t represents customer flows. According to the zero energy loss assumption, the pressures of pipes connected to the node n are expected to be equal to the node pressure pn,t . On the other hand, compressor stations raise the pressures of the pipes connected to them. Therefore, two sets of pipes exist in our model, active pipes with compressors Lact and passive pipes without compressors Lpas . For l ∈ Lact , the inlet pressure is the sum of the starting node pressure pn,t and the boost pressure of the compressor ∆pl,t located

8

2344"5

! )1*+

"#.0

! )"*+*./

! )"*+*,

!

"#$-% %

"#$&'( %

732+89:; ! )1*+

"#.6

Figure 1: Pipeline network node illustration

at the beginning of the pipes; for l ∈ Lpas , the inlet pressure is the starting node pressure pn,t . This leads to: ¯ pl,t,1 = pn,t + ∆pl,t , n ∈ N , l ∈ Lact ∩ Lout n ,t ∈ T ¯ pl,t,1 = pn,t , n ∈ N , l ∈ Lpas ∩ Lout n ,t ∈ T ¯ pl,t,Nz = pn,t , n ∈ N , l ∈ Lin n ,t ∈ T

(8a) (8b) (8c)

For boundary conditions, we fix the pressures at supply nodes and require that the customer flows are greater than customer demands. Moreover, pressures at demand nodes are also required to be greater than minimal contract pressures, pn,t = ps ,

n ∈ Ns , t ∈ T¯

customer fn,t ≥ dn,t , pn,t ≥ pdmin , n

(9a)

n ∈ Nd , l ∈ ∈ Nd , t ∈ T¯

Lin n ,t

∈ T¯

(9b) (9c)

where ps is supply node pressure, dn,t is customer demand flows, and pdmin is minimal contract pressures. Compressor power consumption for active pipes is calculated by:   γ−1 cp Tgas  pn,t + ∆pl,t  γ powerl,t = fl,t,1 −1 , ηe pn,t

¯ n ∈ N , l ∈ Lact ∩ Lout n , t ∈ T (10)

where powerl,t is compressor power, ηe is compressor energy efficiency, γ is the ratio of heat capacities. In this isothermal model, we assume that the temperature is constant at Tgas . While heat capacities cp and cv are also influenced by the gas composition, we observed less than a 5% variation over the entire gas composition range. As a result, only nominal heat capacity values were used in this study. 9

2.2. Thermodynamics for compressibility factor During the dynamic process, the gas pressure changes over time in every segment. When gas pressure ranges from 37 bar to 75 bar at 293.15 K, the Z-factor, used to describe the deviation of real gas from idea gas, would change from roughly 0.95 to 0.85. Considering the nonlinear nature of the model, this should have a significant impact on our final results. Equations of state (EOS), such as the Peng-Robinson (PR) method [38] and Redlich-Kwong (RK) method [39], have been commonly used to describe fluid properties including Z-factor. Yet this comes with two drawbacks. First, the equations are often in a form of cubic function with three roots. When multiple real roots exist, extra efforts are required to identify the proper root depending on its phase state [40]. Besides, the methods are general-purpose tools designed for all the gases; therefore they could be either inaccurate or inefficient for specified gases and mixtures. Rather than apply general-purpose models designed for all gases, we employ simpler empirical models, which are specific to our system and range of temperature and pressure. Popular empirical models for natural gas include the AGA, Papay and AGA8-DC92 [14]. The AGA model is valid for pressures up to 70 bar, while the other mentioned models have a wider ranges of validity. Given that natural gas exists only in the gaseous phase and methane remains to be the main component in this study, we selected the AGA model due to its range of validity for our problem of interest. Because natural gas is a mixture of different gases, pseudo-critical pressure and temperature are used, which is the weighted sum of critical properties of all the components. The resulting thermodynamic model is given as follows: prl,t,k − 0.533 r , T l ∈ L, t ∈ T , k ∈ X

0.257prl,t,k

Zl,t,k = 1 + pl,t,k prl,t,k = c , pmix Tgas Tr = c , Tmix X pcmix = xi pci ,

l ∈ L, t ∈ T , k ∈ X

(11a) (11b) (11c) (11d)

i∈Nc c Tmix =

X

xi Tic ,

(11e)

i∈Nc

where prl,t,k is pseudo reduced pressure for the mixture, T r is the pseudo reduced temperature for the mixture, pcmix is the pseudo critical pressure for the mixture, 10

c Tmix is the pseudo critical temperature for the mixture, xi is the gas mole fraction for component i, pci is critical pressure for component i, and Tic is critical temperature for component i.

2.3. Bounds and objective function In pipeline operations, the pipeline pressure is constrained by the pipeline property and customer demand. If the pipeline pressure is too low, the pipeline will fail to transport the gas from suppliers to customers as demanded. Otherwise, if the pipeline pressure is too high, the operating cost will rise, along with the risk of pipeline burst and gas leakage. We impose pressure constraints on all the nodes and pipeline segments at the same level except the demand nodes, which are constrained in Eq. 9. On the other hand, although supply pressures are fixed, the supply flow rates are important manipulated variables constrained by contracts and system capacity. Compressors also are required to operate in certain ranges with the corresponding bounds as follows: l ∈ L, t ∈ T¯ , k ∈ X¯ f L ≤ fl,t,k ≤ f U , l ∈ L, t ∈ T¯ , k ∈ X¯ supply sL ≤ fn,t ≤ sU , n ∈ Ns , t ∈ T¯

(12b)

power ≤ powerl,t ≤ power . l ∈ Lact , t ∈ T¯

(12d)

pL ≤ pl,t,k ≤ pU ,

L

U

(12a) (12c)

The first constraint sets the lower and upper bounds of pressures for all the pipelines over time. The second constraint imposes bounds on nodal flow rates, with the third constraint setting the bounds for supply flow rates. The last constraint is used to bound the compressor power for all the active pipes. We assume that the compressor is driven by electricity, and the object of this study is to minimize the electricity cost. If electricity cost is the only term in the objective function, the model will deplete the gas in the pipeline to minimize the cost. In order to maintain stable operation, some researchers have penalized the deviation of the total mass of gas in the networks. Since the initial state changes all the time with only total mass constraints, the control profiles should be calculated for each time horizon. Moreover, the continuously changing initial state may contribute to convergence failure at some point. On the other hand, if we add cyclic state constraints to the model, the degree of freedom will be negative. Thus, penalty terms of pressures and flow rates are added to the objective function in order to penalize the deviation of the final pipeline state from the initial state. The resulting objective function is as follows: 11

φ = Ce

X X l∈Lact t∈T¯

powerl,t ∆τ + Cp

X X pl,1,k − pl,N ,k X X fl,1,k − fl,N ,k t t ( ( )2 + Cf )2 . ¯ p¯ f ¯ ¯ l∈L l∈L k∈X

k∈X

(13) where Ce is the unit cost of electricity, Cp is the penalty parameter for terminal pressure constraint and Cf is the penalty parameter for terminal flow constraint, p¯ and f¯ are average pressure and mass flow rate, respectively, which are used to make the tracking terms unitless. For the case studies we selected penalty parameter values to 10000. Penalty parameters should be large enough so that the penalty terms will approach zero, which indicates that the final state is close to the initial state. On the other hand, if penalty terms are too large, the objective function will put too much emphasis on the penalty terms, leading to poor economic performance. The initial conditions can be obtained by solving the steady-state version of this model, and the entire pipeline network optimization model is thus formulated as an NLP with the following structure: min Objective function, Eq. (13) s.t. pipeline equations, Eq. (6) Node equations, Eq. (7), (8) Boundary conditions, Eq. (9) Compressor equation, Eq. (10) Equation of state, Eq. (11) Bounds, Eq. (12)

(14a) (14b) (14c) (14d) (14e) (14f) (14g)

3. Optimization strategy under uncertainty As described in the previous sections, there are two types of strategies for optimization under uncertainty, stochastic optimization and robust optimization. The stochastic optimization approach is computationally prohibitive when a large number of scenarios are considered. Therefore, the robust optimization approach based on back-off constraints is applied here to handle the uncertainty.

12

3.1. Concept of back-off constraints Consider an optimization problem: min φ(z, y, u, t) ¯ s.t z˙ = f (z, y, u, t, θ), ¯ =0 h(z, y, u, t, θ) ¯ ≤0 g(z, y, u, t, θ)

z(0) = z0

(15)

where z represents differential state variables, z0 is the initial condition, y represents algebraic state variables, u is the manipulated variable vector, t is time and θ¯ is the nominal value of uncertain parameters. Now consider an equivalent problem with a more concise formulation: minu F (u, θ) s.t gi (u, θ) ≤ 0,

i∈I

(16)

where gi (u, θ) ≤ 0, i ∈ I are the inequality constraints, and all the equality constraints and state variables are eliminated. In the absence of uncertainty, the problem ¯ the nominal value of θ. With active inequality constraint set i ∈ A, is solved with θ, the KKT conditions are given by X ¯ + ¯ = 0, ∇u F (u, θ) λi ∇u gi (u, θ) (17a) i∈A

gi (u, θ) ≤ 0, i ∈ I ¯ ≤ 0. i ∈ A 0 ≤ λi ⊥ gi (u, θ)

(17b) (17c)

In the presence of uncertainty, we consider the linearization at the nominal value ¯ θ: ¯ + ∇θ F (u, θ) ¯ T (θ − θ), ¯ F (u, θ) ≈ F (u, θ) ¯ + ∇θ gi (u, θ) ¯ T (θ − θ). ¯ g(u, θ) ≈ g(u, θ) i∈I

(18a) (18b)

¯ and ∇θ g(u, θ) ¯ are determined by sensitivity methods. Even though, where ∇θ F (u, θ) ¯ is large. Here, we use back-off the approximation may be inaccurate when ||θ − θ|| ¯ T (θ − θ), ¯ and g˜(u; θ) to replace ∇θ g(u, θ) ¯ T (θ − θ). ¯ term F˜ (u; θ) to replace ∇θ F (u, θ) Shi et al. [32] showed that the the following conditions are necessary for convergence of the back-off procedure to the optimum: ¯ >> ||∇u F˜ (u; θ)|| ≈ 0, ||∇u F (u, θ)|| ¯ >> ||∇u g˜i (u; θ)|| ≈ 0. i ∈ I ||∇u gi (u, θ)|| 13

(19a) (19b)

After replacing the back-off terms with their approximations, the optimization problem becomes: ¯ + F˜ (u; θ), minu F (u, θ) (20) ¯ + g˜i (u; θ) ≤ 0. i ∈ I s.t gi (u, θ) With the assumption (19) satisfied, we have X ¯ + ¯ = 0, ∇u F (u, θ) λi ∇u gi (u, θ)

(21a)

i∈A

gi (u, θ) + g˜i (u; θ) ≤ 0, i ∈ I ¯ + g˜(u; θ) ≤ 0. i ∈ A 0 ≤ λi ⊥ g(x, θ)

(21b) (21c)

Consider an optimization model with KKT conditions as (21): ¯ min φ(u, θ) ¯ + bc ≤ 0 s.t g(u, θ)

(22)

where bc is the back-off term, and the subscript c is the confidence level, indicating the probability of constraints gi (u∗ , θ) satisfied under uncertainty (P [gi (u∗ , θ) + bc ≤ 0] = c). The back-off terms are determined by the variance of gi (u∗ , θ). Srinivasan et al. used Monte Carlo simulation to calculate the variance of the constraints and approximate the back-off terms [26], which is given as: Pm

g(u∗ , θi ) g¯(u , θ) = , m Pm (g(u∗ , θi ) − g¯(u∗ , θ)2 S 2 = i=1 , m−1 bc = ηS. ∗

i=1

(23a) (23b) (23c)

where m is the number of sampling points, parameter η corresponds to the confidence level parameter c, and S 2 is the variance of g(u∗ , θ). If g(u, θ) is normally distributed, η can be selected as the quantiles of the normal distribution, e.g., when η equals 3, ¯ − 3S ≤ g(u∗ , θ) ≤ g(u∗ , θ) ¯ + 3S] ≈ 99.7%. P [g(u∗ , θ) 3.2. Back-off constraints optimization strategy To determine values of the back-off terms, we combine Monte Carlo simulation along with an iterative strategy in order to reduce the number of samples and guarantee feasibility of constraints under uncertainty. The strategy to solve the back-off 14

constraints problem is illustrated in Fig. 2 and is adapted from [32]. Several modifications are made in order to better address the problem of interest. We combine the differential variables z and algebraic variables y into state variables x, and the strategy can be described as follows: ¯ and obtain 1: Solve the optimization problem at nominal conditions with θ = θ, ∗ ∗ optimal control variables u and state variables x . 2: Fix the control variables u∗ , generate values of θ from an uncertainty distribution, initialize the model with state variables x∗ from Step 1, remove variable bounds, perform simulations, and calculate the back-off terms using Eq. (23). 3: Solve the optimization with back-off terms and obtain new control variables u∗ and state variables x∗ . 4: Repeat Step 2 with updated control variables u∗ and state variables x∗ . 5: If there are no constraint violations in Step 4 and back-off terms don’t change (within a tolerance), stop; Else, go back to Step 3 with updated back-off terms.

Figure 2: Back-off constraints optimization strategy framework

We consider two major modifications to the strategy in [32]. First, in Step 3. the Monte Carlo Simulation model is initialized with state variables x∗ from the previous step and variable bounds are removed. The purpose of doing this is twofold: 1) to guarantee convergence of the simulation model when the degrees of freedom are zero; 2) to make sure the model converges to local solutions that are close enough, 15

such that the calculated back-off terms don’t fluctuate wildly and fewer iterations are needed. Second, we modify the terminal conditions with the more strict rule that no constraint violations are allowed in the Monte Carlo Simulation step. These modifications prove to be effective from our numerical tests, as well as from the case studies in the next section. This approach guarantees a robust solution without excessive sampling. Here, both constant and time-varying back-off terms are determined. Constant back-off terms are calculated by choosing the largest variance of the constraints in the time horizon, bc = ηmax(S(t)), while time-varying back-off terms are calculated at each time segment, bc (t) = ηS(t). Because of the additional variables, the time-varying back-off method is supposed to be more time-consuming but leads to less conservative solutions. The following section will discuss their performance differences. 4. Natural gas pipeline network case studies Gas transmission pipeline networks are used to move massive amounts of highpressure sales gas from processing plants to distribution systems and large industrial users [1]. In this study, we mainly consider uncertainties arise from two sources. One source is the uncertainty of demand and the other is the composition of supply gas. Demand uncertainty of pipeline network is extensively studied, which stems from several factors such as weather, emergency or economic concerns. The uncertainty in gas composition, on the other hand, is less studied. Several factors can contribute to the variation of gas composition in pipeline networks, including gas biogenic origins and treatment processes. Moreover, gas refineries often produce natural gas with specified heating value rather than compositions, therefore, pipeline networks may receive gases with the same heating value but different compositions. In this section, we consider several case studies to demonstrate the proposed model and algorithm. The first three case studies consider a tree-shaped small network. Case Study 4 is used to verify the scalability of the model for a large network. All the models are formulated as NLPs in Pyomo 5.3 [41] [42] and solved by using the interior-point solver IPOPT 12.9 [43] on an Intel Core i7 CPU @ 2.9 GHz Macbook Pro running MacOS 10.14. HSL MA97 is used as the linear solver [44], and Pyomo.DAE is used to discretize the PDE [45]. 4.1. Small network case studies As shown in Fig. 3, a small network with 1 supply node, 5 demand nodes, 6 pipes, and 3 compressor stations is used to demonstrate the model and algorithm. Table 1 shows the parameters and bounds of the model. The demands are 25, 37, 16

37, 25, 25 kg/s for nodes 3-7, respectively, and the demand pressures are required to be greater than 40 bar.

Figure 3: Small pipeline network topology

Table 1: Physical parameters and variable bounds

Parameters Ll Tgas cv Cp

100 km 293.15 K 1.85 kJ/kg/K 10000

Dl l ηl Cf

0.9 m R −5 5 × 10 m cp 85% Ce 10000 η

Bounds pL sU

34 bar 200 kg/s

pU fL

70 bar 0 kg/s

sL fU

8.314 J/mol/K 2.34 kJ/kg/K 0.1 $/kW h 3

0 kg/s 200 kg/s

As described in the previous section, each pipe is discretized by using finite difference scheme in space and implicit Euler in time. The time horizon is 24 hours and 48 time segments are used. As for the spatial discretization resolution, the number of segments heavily influences the model accuracy and the computational efficiency. Generally, increasing discretization points leads to more accurate models but longer computation time. Therefore, a trade-off should be made to strike a good balance between accuracy and efficiency. To resolve this, we conducted a simulation experiment of the impact of spatial discretization. The model is simulated with up to 100 discretization points for each pipe. Fig. 4 shows the pressure profiles at node 4 with varying discretization, and other nodes have similar trends. As the increase of the number of discretization points, the pressure profiles have a tendency to be 17

lower. However, the pressure profile with 10 discretization points is very close to the pressure profiles with more discretization points. Given the large number of simulations required, the model is thus discretized with 10 spatial segments for each pipe. Nevertheless, this model can be easily extended to higher-resolution with larger computational budget or using parallelization [36].

Nx = 2 Nx = 3 Nx = 4 Nx = 5 Nx = 10 Nx = 20 Nx = 30 Nx = 40 Nx = 50 Nx = 100

44

pressure/[bar]

43 42 41 40 0

10

20 30 time/[h]

40

50

Figure 4: Impact of spatial discretization on pressure at node 4

4.1.1. Case Study 1: small network optimization under demand uncertainty In this case study, we assume that the demand is uncertain and the supply gas has a constant composition as shown in Table 2 1 . The demand profile of each demand node is divided into three stages, the first stage and the last stage are known and the second stage is uncertain, which is as follows:

dn,t =

 

0 ≤ t < 10.5 10.5 ≤ t ≤ 16.5 16.5 < t ≤ 24

dn θ d × dn  dn

Table 2: Gas compositions

Component N2 CH4 C2 H6 C3 H8 n-C4 H10

Mole fraction % 0.56 87.98 9.00 1.99 0.47 18

(24)

where dn is the nominal demand for node n, θd is the uncertainty parameter for demand, and assumed to behave according to Gaussian distribution θd ∼ N (1.2, 0.03). As illustrated in Fig. 5, the shadowed areas represent the uncertain demand regions within three times the variance. In the Monte Carlo simulation step, 200 points are sampled as shown in Fig. 6. The demand pressures are required to be over 40 bar to satisfy the contract.

nodes 3,6,7 nodes 4,5 3 σ deviation

demand/[kg/s]

45 40 35 30 25 0

5

10 15 time/[h]

20

25

Figure 5: Demand profile for small network

After discretization, there are 26916 variables 26769 equality constraints and 245 inequality constraints. As described in Section 2, the NLP model is first solved under nominal conditions to obtain optimal control variables u∗ and state variables x∗ . In the simulation step, we fix the control variables at u∗ , and remove the inequality constraints to obtain a square matrix. In order to guarantee convergence of every realization of process uncertainty, the model is initialized with x∗ , and variable bounds are removed. Fig. 7 shows the simulation results, where 82 out of 200 samples fail to satisfy the pressure constraints for at least one demand node. According to Fig. 7, three trends can be observed. First, for each demand node, the pressure profiles overlap before entering the uncertain region, to an extent that the pressure profiles are almost the same for different realizations. After entering the uncertain regions, the widths of pressure profile band expand quickly and peak at the end of the uncertain region. Second, the gas pressures are higher in the first stage and decrease in the later stages. Finally, the constraint violations are mostly 1

We assume the gas is from Algeria - Arzew, and its composition information is published in

[46]

19

found in nodes 4, 5, 7, and nodes 3 and 6 have very few violations. The first trend can be explained by the demand profile shown in Fig. 5, in which the demands are the same for all the simulations in the first and third stages. During the first stage, all simulations share the same initial conditions and boundary conditions. After entering the second stage, different realizations of the demand profile lead to different demand node pressure profiles. At the beginning of the third stage, the demand profiles start to be the same again, which explains why the pressure profile band widths stop growing after the second stage. As for the second trend, the demand change between the first and the second stage is so fast that the pipeline inventory must be built up before entering the second stage. Therefore, higher pressures are required in the first stage. Because the operator is required to restore the system to the original state at the end of every time horizon, the inventory decreases slowly after the first stage. The third trend can be explained by the relative locations of demand nodes. As illustrated in Fig. 3, node 3 is followed by nodes 4, 5, and node 6 is followed by node 7. The terminal nodes, such as nodes 4, 5 and 7, tend to operate at the lowest possible pressures, and therefore are more vulnerable to process uncertainty.

14

Probability(%)

12 10 8 6 4 2 0

1.125

1.150

1.175 1.200 1.225 Demand uncertaity level

1.250

1.275

Figure 6: Monte Carlo simulation distribution of demand uncertainty level for small network, m = 200

20

50

48

48 pressure/[bar]

pressure/[bar]

50

46 44

46 44

42

42

40

40

38

0

5

10 15 time/[h]

20

38

25

0

5

50

50

48

48

46 44

40

40 5

20

25

44 42

0

25

46

42

38

20

(b) node 4

pressure/[bar]

pressure/[bar]

(a) node 3

10 15 time/[h]

10 15 time/[h]

20

38

25

0

5

(c) node 5

10 15 time/[h]

(d) node 6

50

pressure/[bar]

48 46 44 42 40 38

0

5

10 15 time/[h]

20

25

(e) node 7 Figure 7: Case Study 1: Monte Carlo simulation results under demand uncertainty without back-off terms

After Monte Carlo simulation, the constraint variance S 2 and back-off terms bc can be obtained via Eq.(23), and bc are added to the model constraints. The updated NLP model is solved again to generate new u∗ and x∗ , and the algorithm continues 21

until terminal conditions are satisfied. For this case, the constant back-off method is applied and only 5 back-off terms are used to deal with 5 demand nodes. The algorithm takes 1216 CPU seconds and 2 iterations to converge. Compared with the model under nominal condition, the consideration of demand uncertainty raises the objective function value (including the penalty terms) from 73566.5 to 78629.8, and the power consumption from 292.6 M W h to 295.2 M W h. Fig. 8 shows the Monte Carlo simulation results at the last iteration with back-off terms. All constraints are now satisfied and all pressure profiles of demand nodes maintain some distance from the lower bound at the cost of a generally higher operating pressures.

22

52

50

50

48

48

pressure/[bar]

pressure/[bar]

52

46 44

46 44

42

42

40

40

38

0

5

10 15 time/[h]

20

38

25

0

5

52

52

50

50

48

48

46 44

42 40 5

20

25

44

40 0

25

46

42

38

20

(b) node 4

pressure/[bar]

pressure/[bar]

(a) node 3

10 15 time/[h]

10 15 time/[h]

20

38

25

0

5

(c) node 5

10 15 time/[h]

(d) node 6

52

pressure/[bar]

50 48 46 44 42 40 38

0

5

10 15 time/[h]

20

25

(e) node 7 Figure 8: Case Study 1: Monte Carlo simulation results under uncertainty with back-off terms

23

4.1.2. Case Study 2: small network optimization under demand and composition uncertainty In practice, gas supply comes from a combination of different sources, which makes gas composition uncertain. For this case study we adopt the same demand uncertainty as in Case Study 1 is adopted, and we derive a simple, but comprehensive uncertain description for gas composition. Table 3 is the worldwide average LNG compositions in 2011 [46], from which we can find that the Z-factor is largely proportional to the fraction of methane. To simplify our uncertainty description, we classify the feed gases into two parts: methane and the remaining gases. We assume that the relative ratios of chemical components, other than methane, are constant as shown in Table 2, so that only the fraction of methane, xCH4 , is required to represent gas composition uncertainty. The thermodynamic model Eq. 11 can be modified as prl,t,k Zl,t,k = 1 + 0.257prl,t,k − 0.533 r , l ∈ L, t ∈ T , k ∈ X T pl,t,k r pl,t,k = c , l ∈ L, t ∈ T , k ∈ X pmix Tgas Tr = c , Tmix c pmix = xCH4 pcCH4 + (1 − xCH4 )pcothers , c c c Tmix = xCH4 TCH + (1 − xCH4 )Tothers , 4

(25a) (25b) (25c) (25d) (25e)

c where xCH4 is the fraction of methane, pcCH4 is critical pressure of methane, TCH 4 c c is critical temperature of methane, pothers and Tothers is pseudo critical pressure and pseudo critical temperature for the mixture of components shown in Table 2 without considering methane, respectively. According to Table 3, the mean and variance of methane fractions are 0.9084 and 0.0451. We assume that the methane fractions of supply gas is consistent with Gaussian distribution θc ∼ N (0.9084, 0.0451). Fig. 9 shows uncertain description with 300 sampling points. Because of the symmetric nature of Gaussian distribution, the methane fraction of some points in Fig. 9 are greater than one. While this is physically impossible to have natural gas with methane fraction greater than one, the calculation of Eq.25 is still valid at those points when the natural gas consists of light gases such nitrogen and hydrogen. By retaining points greater than one, we can maintain the assumed mean and variance of Gaussian distribution and extend the variation of gas composition to include gases with decent amount of light components.

24

Based on the simplified composition model and the assumed Gaussian distribution, we calculate the distribution of Z-factors and compare them with the real-world gases in Table 3. Fig. 10 shows the range of Z-factor values using the simplified composition model. The scatter points are Z-factor values of gases from Table 3 at various pressures, and the solid lines represent minimal and maximal values of Z-factor at corresponding pressures calculated with the gas compositions sampled from the assumed Gaussian distribution. This plot shows that the simplified composition model can well represent the real world gas; all Z-factor points from real world gases lie within this range. Table 3: Worldwide average LNG compositions

N2 (%) Algeria – Arzew 0.56 Algeria – Bethioua 1 1.20 Algeria – Bethioua 2 0.92 Algeria – Skikda 1.02 Egypt – Damietta 0.08 Egypt – Idku 0.00 Libya 0.69 Nigeria 0.08 Abu Dhabi 0.29 Oman 0.35 Qatar 0.36 Trinidad 0.03 USA – Alaska 0.17 Australia – NWS 0.09 Brunei 0.05 Indonesia – Arun 0.06 Indonesia – Badak 0.02 Malaysia 0.16

2

CH4 (%) 87.98 87.59 91.39 91.19 97.70 97.20 81.57 91.28 84.77 87.89 90.10 96.82 99.73 87.39 90.61 91.16 89.76 91.15

Calculated by AGA method at 293.15 K, 37 bar.

25

C2 H6 (%) 9.00 8.39 7.17 7.02 1.80 2.30 13.38 4.62 13.22 7.27 6.23 2.74 0.08 8.33 4.97 6.01 5.06 4.96

C3 H8 (%) 1.99 2.12 0.52 0.66 0.22 0.30 3.67 2.62 1.63 2.92 2.32 0.31 0.01 3.35 2.89 1.84 3.54 2.79

C4+ (%) 0.47 0.70 0.00 0.11 0.20 0.20 0.69 1.40 0.09 1.57 0.99 0.10 0.00 0.84 1.48 0.93 1.62 0.94

Z-factor2 0.9075 0.9077 0.9163 0.9158 0.9245 0.9234 0.8952 0.9093 0.9027 0.9037 0.9091 0.9231 0.9287 0.9033 0.9077 0.9107 0.9052 0.9099

1.275 demand parameter

1.250 1.225 1.200 1.175 1.150 1.125 1.100 0.80

0.85 0.90 0.95 fraction of methane

1.00

Figure 9: Case Study 2: Monte Carlo simulation distribution of demand and composition uncertainty level, m=300

min max

0.92

Z-factor

0.90 0.88 0.86 0.84 0.82 0.80 40

45

50 55 Pressure [bar]

60

65

70

Figure 10: The range of gas Z-factor values for simplified gas composition model

We first use Monte Carlo simulation to check the feasibility of the new problem without back-off terms.Fig. 11 shows the pressure violations of demand nodes under gas composition and demand uncertainty, where 171 out of 300 simulations violate the constraints. Fig. 12 shows the minimal pressure over time, which has similar constraint violations as Case Study 1.

26

node 3

120 100 80 60 40 20 0

node 4

140 number of violations

number of violations

140

120 100 80 60 40 20

0

5

10 15 time/[h]

20

0

25

0

5

(a) node 3

number of violations

number of violations

25

node 6

140

120 100 80 60 40 20 0

20

(b) node 4

node 5

140

10 15 time/[h]

120 100 80 60 40 20

0

5

10 15 time/[h]

20

0

25

0

5

(c) node 5

20

25

(d) node 6

node 7

140 number of violations

10 15 time/[h]

120 100 80 60 40 20 0

0

5

10 15 time/[h]

20

25

(e) node 7 Figure 11: Case Study 2: Number of constraint violations over time for small network without back-off terms

This case with back-off terms requires 3685 CPU seconds and 3 iterations to converge with all constraints satisfied. Compared with the solution at nominal conditions, the objective function value increases from 76073 to 83527, and the com27

pressor power consumption increases from 292.7 to 297.1 M W h. Compared with the result in Case Study 1, the consideration of gas composition uncertainty raises the objective function value by 6.3%, and the compressor power consumption by 0.60%, which indicates that the consideration of gas composition uncertainty is mainly at the cost of the penalty terms.

28

node 3

node 4

55 minimal pressure/[bar]

minimal pressure/[bar]

55

50

45

40

35

50

45

40

35 0

5

10 15 time/[h]

20

25

0

5

(a) node 3

20

25

(b) node 4

node 5

node 6

55 minimal pressure/[bar]

55 minimal pressure/[bar]

10 15 time/[h]

50

45

40

35

50

45

40

35 0

5

10 15 time/[h]

20

25

0

5

(c) node 5

20

25

(d) node 6

node 7

55 minimal pressure/[bar]

10 15 time/[h]

50

45

40

35 0

5

10 15 time/[h]

20

25

(e) node 7 Figure 12: Case Study 2: Minimal pressures over time with demand and composition uncertainty without back-off terms

29

4.1.3. Case Study 3: small network optimization under demand and composition uncertainty with time varying back-off constants The constant back-off method is easy to implement and has fewer variables added to the model. On the other hand, pressures at different times have different profiles. As shown in Fig. 12, the minimal pressure is lower at the second and third stages than the first stage. If back-off terms are time-varying, tighter constraints are expected. In this case study, the model in Case Study 2 is used but optimized with time-varying back-off terms, where the back-off terms for different demand nodes are set to be time-dependent. Compared with Case Study 2, 245 back-off terms are used for 5 demand nodes and 49 time steps. The back-off approach from Fig. 2 requires 3 iterations and 3672 CPU seconds for the algorithm to converge, which is similar to Case Study 2. However, its objective function value with back-off terms decreases from 83527.3 to 82876.4 M W h, and the compressor boost cost decreases from 297.1 to 296.1 M W h. Fig. 13 shows the comparison of back-off term values of the two methods.

30

1.2

0.8 0.6 0.4 0.2 0.0

time-varying back-off terms constant back-off terms

1.0 back-off term value

1.0 back-off term value

1.2

time-varying back-off terms constant back-off terms

0.8 0.6 0.4 0.2

0

5

10 15 time/[h]

20

0.0

25

0

5

(a) node 3

1.2

back-off term value

back-off term value

20

25

1.0

0.6 0.4

0.8 0.6 0.4

0.2

0.2

0.0

0.0

5

25

1.2

0.8

0

20

(b) node 4

time-varying back-off terms constant back-off terms

1.0

10 15 time/[h]

10 15 time/[h]

20

25

time-varying back-off terms constant back-off terms 0

5

(c) node 5

10 15 time/[h]

(d) node 6

1.2

back-off term value

1.0 0.8 0.6 0.4 0.2 0.0

time-varying back-off terms constant back-off terms 0

5

10 15 time/[h]

20

25

(e) node 7 Figure 13: Comparison of back-off terms between constant and time-varying back-off methods

Table 4 compares the three case studies and shows that the consideration of uncertainty raises the computational time significantly. Particularly, the composition uncertainty makes the problem much more complicated, tripling the CPU time of 31

Case Study 1. However, the robustness of the model, with satisfaction of all constraints, is achieved with only 1.2% increase in economic cost and 12.6% of the overall objective function value. We note that the time-varying back-off method can obtain a slightly better result than constant back-off method with similar computing time. Table 4: Comparison of results for small network

Case Nominal case Case Study 1 Case Study 2 Case Study 3

OBJ

Compressor power consumption (M W h) 73566.5 292.6 78629.8 295.2 83527.3 297.1 82876.4 296.1

CPU time (s) 11 1216 3685 3672

4.2. Case Study 4: Large network optimization under demand and composition uncertainty In order to test the scalability of the model, a large pipeline network adapted from [3] is studied. As illustrated in Fig. 14, the model consists of 3 supplier nodes 1, 15, and 16; 12 demand nodes with different demand patterns; 5 compressors on pipes 1, 8, 14, 20, and 27. We assume that the gas is supplied at a constant pressure of 30 bar. For demand nodes both constant and periodic demand patterns are considered: nodes 4 and 25 have a constant demand of 10 kg/s, nodes 30, 32, 34, 35 have a constant demand of 5 kg/s, and nodes 6, 12, 13, 19, 21, 23 have a sinusoidal demand pattern with a 24 h period and an amplitude of 5% around the nominal value at 25 kg/s, which is illustrated in Fig. 15. The demand node pressures are required to be higher than 40 bar. To simulate the real demand, an uncertain demand disturbance is added to the sinusoidal demand patterns (nodes 6, 12, 13, 19, 21) through Gaussian noise θd ∼ N (0, 0.15). In addition, composition uncertainty of methane θc ∼ N (0.9084, 0.0451) is also added to the model. Otherwise, the same model parameters are used as the previous case studies, and only the time-varying back-off method is applied in this case. Therefore, 588 back-off terms are used for 12 demand nodes and 49 time steps, leading to a model with 149388 variables, 149143 equality constraints and 588 inequality constraints.

32

Figure 14: Large pipeline network topology

33

30

demand/[kg/s]

25 20 n30,32,34,35 n4,25 n6,12,13,19,21,23 3 σ deviation

15 10 5 0

5

10 15 time /[h]

20

25

Figure 15: Large pipeline network demand profile

4

demand parameter

3 2 1 0 −1 −2 −3 −4 0.80

0.85 0.90 0.95 fraction of methane

1.00

Figure 16: Case Study 4: Monte Carlo simulation distribution of uncertain parameters for large network, m=300

Fig. 17 shows the Monte Carlo simulation results without back-off terms. According to the results, 160 out of 300 simulations violate the constraints at least once. Compared with the small network case studies, the pressure profiles have smoother changes and the widths of pressure bands are identical over the whole time horizon. This can be explained by the demand profile features of the model. First, the demand disturbance here is not as volatile as the previous cases, which can be easily observed from the demand profiles. Second, the demand uncertainty applies over the

34

whole time horizon, resulting in similar pressure variance at each time. The algorithm in Fig. 2 converges after 3 iterations with 20162 CPU seconds, to a point that satisfies all constraints. The objective function value increases 14.6% from 26627.0 to 30509.5, and the compressor power consumption increases 14.4% from 244.1 to 279.2 M W h. Fig. 18 shows the back-off term values over time.

35

50

48

48

46 44 42 40 38

pressure/[bar]

50

48

pressure/[bar]

pressure/[bar]

50

46 44 42 40 38

0

5

10

time/[h]

15

20

25

40

5

10

time/[h]

15

20

25

0

50

48

48

46 44 42 40

pressure/[bar]

50

48

46 44 42 40 38

0

5

10

time/[h]

15

20

25

5

10

time/[h]

15

20

25

0

pressure/[bar]

48

pressure/[bar]

50

46 44 42 40 38

5

10

time/[h]

15

20

25

5

10

time/[h]

15

20

25

0

pressure/[bar]

48

pressure/[bar]

50

46 44 42 40 38

5

10

time/[h]

15

(j) node 32

5

20

25

10

time/[h]

15

20

25

20

25

(i) node 30

48

38

25

40

50

40

20

42

48

42

15

44

50

44

time/[h]

46

(h) node 25

46

10

38 0

(g) node 23

0

5

(f) node 21

48

38

25

40

50

40

20

42

48

42

15

44

50

44

time/[h]

46

(e) node 19

46

10

38 0

(d) node 13

0

5

(c) node 12

50

pressure/[bar]

pressure/[bar]

42

(b) node 6

38

pressure/[bar]

44

38 0

(a) node 4

pressure/[bar]

46

46 44 42 40 38

0

5

10

time/[h]

15

(k) node 34

20

25

0

5

10

time/[h]

15

(l) node 35

Figure 17: Case Study 4: Monte Carlo simulation results under demand and composition uncertainty for large network without back-off terms

36

3.3

3.2

3.2

3.1 3.0 2.9 2.8 2.7

back-off term value

3.3

3.2

back-off term value

back-off term value

3.3

3.1 3.0 2.9 2.8 2.7

0

5

10

time/[h]

15

20

25

2.8

5

10

time/[h]

15

20

25

0

3.3

3.2

3.2

3.1 3.0 2.9 2.8

back-off term value

3.3

3.2

3.1 3.0 2.9 2.8 2.7

0

5

10

time/[h]

15

20

25

5

10

time/[h]

15

20

25

0

back-off term value

3.2

back-off term value

3.3

3.1 3.0 2.9 2.8 2.7

5

10

time/[h]

15

20

25

5

10

time/[h]

15

20

25

0

back-off term value

3.2

back-off term value

3.3

3.1 3.0 2.9 2.8 2.7

5

10

time/[h]

15

(j) node 32

5

20

25

10

time/[h]

15

20

25

20

25

(i) node 30

3.2

2.7

25

2.8

3.3

2.8

20

2.9

3.2

2.9

15

3.0

3.3

3.0

time/[h]

3.1

(h) node 25

3.1

10

2.7 0

(g) node 23

0

5

(f) node 21

3.2

2.7

25

2.8

3.3

2.8

20

2.9

3.2

2.9

15

3.0

3.3

3.0

time/[h]

3.1

(e) node 19

3.1

10

2.7 0

(d) node 13

0

5

(c) node 12

3.3

back-off term value

back-off term value

2.9

(b) node 6

2.7

back-off term value

3.0

2.7 0

(a) node 4

back-off term value

3.1

3.1 3.0 2.9 2.8 2.7

0

5

10

time/[h]

15

(k) node 34

20

25

0

5

10

time/[h]

15

(l) node 35

Figure 18: Case Study 4: Time-varying back-off terms for large network with demand and composition uncertainty

37

5. Conclusions We develop a dynamic model for pipeline network optimization with empirical equations of state, which allow accurate calculation of gas compressibility factor. A robust optimization algorithm based on back-off constraints is implemented to deal with the demand and composition uncertainty in natural gas pipeline networks, with minimum compression costs. Several case studies are used to demonstrate the model and algorithm. In Case Study 1, demand uncertainty is considered and a robust solution is obtained with an increase of objective function value and power consumption by 6.9% and 0.9%, respectively. In Case Study 2, we test the model in the presence of inlet gas composition and demand uncertainty, and a constant backoff method is applied to solve the problem. Compared with the result in Case Study 1, the consideration of gas composition uncertainty raises the objective function value by 6.3%, and the compressor power consumption by 0.60%. From Case Study 2, we observe that large offset areas exist between the minimal pressures and the constraints in the pressure profiles. Therefore, a time-varying back-off method is applied to resolve the problem, which generates a solution with tighter offsets and lower operating cost with similar computing time. To verify the scalability of the model, a larger network is considered in Case Study 4, and the time-varying back-off method obtains a robust optimal solution. Our case studies show that the algorithm can generate a robust solution at the cost of economic performance, compressor power consumption, with 1.2% for the small network, and 14.4% for the large network when considering both demand and inlet gas composition uncertainty. This work has dealt with robust optimization of dynamic natural gas pipeline networks with demand and inlet gas composition uncertainty. The developed model can be easily extended to even larger pipeline networks due to its flexibility and the platform it is built on. With rigorous thermodynamics, the model can also be used to address pipelines with multiple supply sources, where gases with different compositions mix in the networks. Finally, parallel computing and model decomposition methods can be applied to study models with higher resolutions and additional uncertain parameters. 6. Acknowledgements The authors gratefully acknowledge financial support from China Scholarship Council, the National Natural Science Foundation of China (No. 21776323), as well as the support from Center for Advanced Process Decision-making at Carnegie Mellon University.

38

Nomenclature Greek symbols θ¯ nominal value of uncertain parameters ∆τ discretization segment length in time, h ∆pl,t pressure boost, bar ∆z discretization segment length in space, km l pipe rugosity, m η tuning parameter for conservatism ηe compressor efficiency γ ratio of heat capacities cp and cv λl (τ, z) friction coefficient in continuous form λl,t,k friction coefficient in discretized form ν dynamic viscosity of gas, kg/(ms) ρl (τ, z) gas density in continuous form, kg/m3 ρl,t,k gas density in discretized form, kg/m3 τ time dimension, h θ uncertain parameters θc uncertain parameter for gas composition θd uncertain parameter for demand Cf penalty parameter for terminal flow constraint, M W h/(kg/s) Cp penalty parameter for terminal pressure constraint, M W h/bar vl (τ, z) gas speed in continuous form, m/s Latin symbols T¯ set of discretization points in time X¯ set of discretization points in space f¯ average gas flow rate, kg/s p¯ average pressure, bar L set of pipes in the pipeline network in Ln set of pipes that send gas to node n out Ln set of pipes that receive gas from node n Lact set of pipes with compressor stations Lpas set of pipes without compressor stations N set of nodes in pipeline network 39

Nd Ns T Xl Al bc Ce cp Dl dn,t dn fL fU customer fn,t

set of demand nodes set of supply nodes time horizon, h axial dimension, m cross-sectional area, m2 back-off terms with confidence level of c electricity unit price, $/kW heat capacity at constant pressure, kJ/(kgK) pipe diameter,m customer demand flow rate, kg/s nominal demand gas flow for node n, kg/s lower bound for flow rate, kg/s upper bound for flow rate, kg/s gas flow to customers, kg/s

supply fn,t fl (τ, z) fl,t,k g hl (z) Ll M Nt Nz pcCH4 Pic c Pmix pcothers pL r Pl,t,k U p pl (τ, z) pl,t,k pdmin pn,t ps

gas flow from supply nodes, kg/s gas flow rate in continuous form, kg/s gas flow rate in discretized form, kg/s gravitational acceleration, m/s2 pipe height, m pipe length, km molecular weight discretization number in time discretization number in space critical pressure of methane, bar pseudo critical pressure for component i, bar pseudo critical pressure for mixture, bar pseudo critical pressure of gas components other than methane, bar lower bound for pressure, bar pseudo reduced pressure, bar upper bound for pressure, bar gas pressure in continuous form, bar gas pressure in discretized form, bar minimal pressure allowed for demand nodes, bar pressure at node n, bar pressure at supply node, bar 40

powerL powerU powerl,t R Rel (τ, z) S sL sU T Tic c Tmix c Tothers Tr Tl (τ, z) Tgas Tsoil U u u∗ x x∗ xCH4 xi y z

lower bound for compressor power, M W upper bound for compressor power, M W compressor power, M W gas constant, J/(molK) Reynolds number square root of constraint variance lower bound for supply flow rate, kg/s upper bound for supply flow rate, kg/s final time, h pseudo critical temperature for component i, K pseudo critical temperature for mixture, K pseudo critical temperature of gas components other than methane, K pseudo reduced temperature, K gas temperature, K temperature, K soil temperature, K heat transfer coefficient, W/(m2 K) control variables optimal control variables state variables optimal state variables methane fraction gas mole fraction for component i algebraic state variables differential variables in Section 3, and axial dimension in the rest of the paper Zl (τ, z) compressibility factor in continuous form Zl,t,k compressibility factor in discretized form

41

References [1] S. Mokhatab, W. A. Poe, Handbook of natural gas transmission and processing, Gulf Professional Publishing, 2012. [2] R. S. Middleton, R. Gupta, J. D. Hyman, H. S. Viswanathan, The shale gas revolution: Barriers, sustainability, and emerging opportunities, Applied Energy 199 (2017) 88 – 95, ISSN 0306-2619. [3] A. Gopalakrishnan, L. T. Biegler, Economic Nonlinear Model Predictive Control for periodic optimal operation of gas pipeline networks, Computers & Chemical Engineering 52 (2013) 90 – 99, ISSN 0098-1354. [4] R. Z. R´ıos-Mercado, C. Borraz-S´anchez, Optimization problems in natural gas transportation systems: A state-of-the-art review, Applied Energy 147 (2015) 536–555. [5] Y. Puranik, M. Kilin¸c, N. V. Sahinidis, T. Li, A. Gopalakrishnan, B. Besancon, T. Roba, Global optimization of an industrial gas network operation, AIChE Journal 62 (9) (2016) 3215–3224. ¨ [6] H. Uster, S¸. Dilavero˘glu, Optimization for design and operation of natural gas transmission networks, Applied Energy 133 (2014) 56–69. [7] B. Geißler, A. Morsi, L. Schewe, M. Schmidt, Solving highly detailed gas transport MINLPs: block separability and penalty alternating direction methods, INFORMS Journal on Computing 30 (2) (2018) 309–323. [8] M. Mikolajkov´a, H. Sax´en, F. Pettersson, Linearization of an MINLP model and its application to gas distribution optimization, Energy 146 (2018) 156–168. [9] T. Koch, B. Hiller, M. E. Pfetsch, L. Schewe, Evaluating gas network capacities, SIAM, 2015. [10] D. Rose, M. Schmidt, M. C. Steinbach, B. M. Willert, Computational optimization of gas compressor stations: MINLP models versus continuous reformulations, Mathematical Methods of Operations Research 83 (3) (2016) 409–444. [11] F. da Silva Alves, J. N. M. de Souza, A. L. H. Costa, Multi-objective design optimization of natural gas transmission networks, Computers & Chemical Engineering 93 (2016) 212–220.

42

[12] A. Demissie, W. Zhu, C. T. Belachew, A multi-objective optimization model for gas pipeline operations, Computers & Chemical Engineering 100 (2017) 94–103. [13] R. Song, M. Cui, J. Liu, Single and multiple objective optimization of a natural gas liquefaction process, Energy 124 (2017) 19–28. [14] M. Schmidt, M. C. Steinbach, B. M. Willert, High detail stationary optimization models for gas networks, Optimization and Engineering 16 (1) (2015) 131–164. [15] A. Osiadacz, Simulation and analysis of gas networks, Gulf Publishing Company, 1987. [16] A. J. Osiadacz, M. Chaczykowski, Comparison of isothermal and non-isothermal pipeline gas flow models, Chemical Engineering Journal 81 (1-3) (2001) 41–51. [17] H. P. Reddy, S. Narasimhan, S. M. Bhallamudi, Simulation and state estimation of transient flow in gas pipeline networks using a transfer function model, Industrial & engineering chemistry research 45 (11) (2006) 3853–3863. [18] R. Raoni, A. R. Secchi, E. C. Biscaia Jr, Novel method for looped pipeline network resolution, Computers & Chemical Engineering 96 (2017) 169–182. [19] B. Baumrucker, L. T. Biegler, MPEC strategies for cost optimization of pipeline operations, Computers & Chemical Engineering 34 (6) (2010) 900–913. [20] N.-Y. Chiang, V. M. Zavala, Large-scale optimal control of interconnected natural gas and electrical transmission systems, Applied Energy 168 (2016) 226–235. [21] M. Hahn, S. Leyffer, V. M. Zavala, Mixed-Integer PDE-Constrained Optimal Control of Gas Networks, Tech. Rep., Tech. rep. Aug. 2017. url: https://wiki. mcs. anl. gov/leyffer/images/2/27/GasNetMIP. pdf (visited on 12/04/2017), 2017. [22] M. Gugat, G. Leugering, A. Martin, M. Schmidt, M. Sirvent, D. Wintergerst, MIP-based instantaneous control of mixed-integer PDE-constrained gas transport problems, Computational Optimization and Applications 70 (1) (2018) 267–294. [23] A. Herr´an-Gonz´alez, J. De La Cruz, B. De Andr´es-Toro, J. Risco-Mart´ın, Modeling and simulation of a gas distribution pipeline network, Applied Mathematical Modelling 33 (3) (2009) 1584–1600.

43

[24] H. A. Behrooz, R. B. Boozarjomehry, Modeling and state estimation for gas transmission networks, Journal of Natural Gas Science and Engineering 22 (2015) 551–570. [25] G. Guandalini, P. Colbertaldo, S. Campanari, Dynamic modeling of natural gas quality within transport pipelines in presence of hydrogen injections, Applied energy 185 (2017) 1712–1723. [26] B. Srinivasan, D. Bonvin, E. Visser, S. Palanki, Dynamic optimization of batch processes: II. Role of measurements in handling uncertainty, Computers & chemical engineering 27 (1) (2003) 27–44. [27] H. A. Behrooz, Managing demand uncertainty in natural gas transmission networks, Journal of Natural Gas Science and Engineering 34 (2016) 100–111. [28] H. A. Behrooz, R. B. Boozarjomehry, Dynamic optimization of natural gas networks under customer demand uncertainties, Energy 134 (2017) 968–983. [29] I. E. Grossmann, R. M. Apap, B. A. Calfa, P. Garc´ıa-Herreros, Q. Zhang, Recent advances in mathematical programming techniques for the optimization of process systems under uncertainty, Computers & Chemical Engineering 91 (2016) 3–14. [30] M. Vallerio, D. Telen, L. Cabianca, F. Manenti, J. Van Impe, F. Logist, Robust multi-objective dynamic optimization of chemical processes using the Sigma Point method, Chemical Engineering Science 140 (2016) 201–216. [31] M. Diehl, H. G. Bock, E. Kostina, An approximation technique for robust nonlinear optimization, Mathematical Programming 107 (1-2) (2006) 213–230. [32] J. Shi, L. T. Biegler, I. Hamdan, J. Wassick, Optimization of grade transitions in polyethylene solution polymerization process under uncertainty, Computers & Chemical Engineering 95 (2016) 260–279. [33] R. W. Koller, L. A. Ricardez-Sandoval, L. T. Biegler, Stochastic back-off algorithm for simultaneous design, control, and scheduling of multiproduct systems under uncertainty, AIChE Journal 64 (7) (2018) 2379—-2389. [34] J. Maußner, H. Freund, Efficient calculation of constraint back-offs for optimization under uncertainty: A case study on maleic anhydride synthesis, Chemical Engineering Science 192 (2018) 306–317. 44

[35] M. V. Lurie, Modeling of oil product and gas pipeline transportation, WileyVch, 2008. [36] V. M. Zavala, Stochastic optimal control model for natural gas networks, Computers & Chemical Engineering 64 (2014) 103–113. [37] M. C. Steinbach, On PDE solution in transient optimization of gas networks, Journal of computational and applied mathematics 203 (2) (2007) 345–361. [38] D.-Y. Peng, D. B. Robinson, A new two-constant equation of state, Industrial & Engineering Chemistry Fundamentals 15 (1) (1976) 59–64. [39] O. Redlich, J. N. Kwong, On the thermodynamics of solutions. V. An equation of state. Fugacities of gaseous solutions., Chemical reviews 44 (1) (1949) 233–244. [40] R. S. Kamath, L. T. Biegler, I. E. Grossmann, An equation-oriented approach for handling thermodynamics based on cubic equation of state in process optimization, Computers & Chemical Engineering 34 (12) (2010) 2085–2096. [41] W. E. Hart, C. D. Laird, J.-P. Watson, D. L. Woodruff, G. A. Hackebeil, B. L. Nicholson, J. D. Siirola, Pyomo–optimization modeling in python, vol. 67, Springer Science & Business Media, second edn., 2017. [42] W. E. Hart, J.-P. Watson, D. L. Woodruff, Pyomo: modeling and solving mathematical programs in Python, Mathematical Programming Computation 3 (3) (2011) 219–260. [43] A. W¨achter, L. T. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming, Mathematical programming 106 (1) (2006) 25–57. [44] HSL. A collection of Fortran codes for large scale scientific computation, http: //www.hsl.rl.ac.uk/, accessed: 2018-10-20, 2016. [45] B. Nicholson, J. D. Siirola, J.-P. Watson, V. M. Zavala, L. T. Biegler, pyomo.dae: a modeling and automatic discretization framework for optimization with differential and algebraic equations, Mathematical Programming Computation 10 (2) (2018) 187–223. [46] I. Union, Petroleum B. guidebook to gas interchangeability and gas quality, BP Publication, London, United Kingdom .

45

Highlights

dynamic pipeline network model is developed with rigorous thermodynamic equations optimal dynamic operation of pipelines is determined to minimize compression costs supply gas composition and demand nodes flow rates are assumed to be uncertain. robust optimization is applied using back-off constraints calculated from Monte Carlo simulation.

Declaration of Interest for “Dynamic optimization of natural gas pipeline networks with demand and composition uncertainty” There are no interests to declare

Author Contributions Section

All authors contributed significantly to this paper.