Dynamic optimization of natural gas networks under customer demand uncertainties

Accepted Manuscript Dynamic optimization of natural gas networks under customer demand uncertainties Hesam Ahmadian Behrooz, R.Bozorgmehry Boozarjomehry PII:

S0360-5442(17)31081-2

DOI:

10.1016/j.energy.2017.06.087

Reference:

EGY 11095

To appear in:

Energy

Please cite this article as: Hesam Ahmadian Behrooz, R.Bozorgmehry Boozarjomehry, Dynamic optimization of natural gas networks under customer demand uncertainties, Energy (2017), doi: 10.1016/ j.energy.2017.06.087 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

EP

C2

AC C

P1, P2, P3 & P4 N1

TE D

S1

P5

D1

D2

N8

P6 N9

ACCEPTED MANUSCRIPT

RI PT

Dynamic optimization of natural gas networks under customer demand uncertainties

Chemical Engineering Faculty, Sahand University of Technology, Tabriz, Iran Department of Chemical and Petroleum Engineering, Sharif University of Technology, Tehran, Iran

AC C

EP

TE D

b

M AN U

a

SC

Hesam Ahmadian Behrooz a,1, R.Bozorgmehry Boozarjomehryb,2

Abstract

1 2

Corresponding author. Tel.: (+9841)33459150, E-mail: [email protected] Tel.: (+9821)66166445, E-mail: [email protected]

1

ACCEPTED MANUSCRIPT

In natural gas transmission networks, the efficiency of daily operation is strongly dependent on our knowledge about the customer future demands. Unavailability of

and quantify the corresponding uncertainty.

RI PT

accurate demand forecasts makes it more important to be able to characterize the loads

A planning strategy for natural gas transmission networks under future demand uncertainty is addressed, which involves coupling of a gas transmission network

SC

dynamic simulator with the stochastic optimization framework. Loads from a gas-fired

M AN U

power plant are studied where the loads are characterized by a number of uncertain parameters, and unscented transform is utilized for uncertainty propagation. This algorithm is compared to deterministic approaches in terms of energy consumption, supply flow rate and line-pack manipulations by means of an illustrative example where 1% increase in energy requirement is observed for stochastic formulation compared to

TE D

expected condition while 10% extra energy is required for the worst case scenario. It is also shown that preparation for uncertain future loads can be done with appropriate management of the available compression units, without further natural gas supplies.

EP

Stochastic programming provides an optimal way to determine the minimum pressure

AC C

buffer required at delivery points for safe operation of the network under customer demand uncertainties.

Keywords: High-pressure gas networks; uncertainty; stochastic optimization; unscented transformation; planning;

2

ACCEPTED MANUSCRIPT

1. Introduction Natural gas travels thousands of miles from the production sites through the

RI PT

transmission lines to the local distribution networks [1]. Large industrial customers such as gas-fired power plants are also supplied from the transmission lines. It is the responsibility of the pipeline monitoring and control center operators (i.e., dispatchers) to ensure that all customers receive adequate amounts of gas at (or above) their

SC

required pressure level. Accurate anticipation of customer needs can provide managers

M AN U

with valuable information that can be helpful in satisfying customer expectations. Hence, these companies attempt to gain insight into their customers’ future needs by the solution of gas demand forecast problem as accurate as possible [2] which can enable them to operate the system smoothly, efficiently and cost effectively. The uncertainties resulting from the estimated demands often play an important role in achieving these

TE D

goals in natural gas industry. Also, it is worth mentioning that there are many uncertainty sources to cope with: modeling uncertainties, measurement uncertainties

EP

and uncertainties resulting from external disturbances [3]. Accurate estimation of the customers future demands coupled with a dynamic gas

AC C

network simulator [4] can be used in an optimization framework to provide the dispatchers with operational guidance [2]. The resulting optimization framework can play an essential role in a successful and safe operation of a gas transmission system due to the significant effect of the future demands on the planning strategy. This optimization framework is composed of three major parts: dynamic simulator, demand forecasting tool and the optimizer [5]. 3

ACCEPTED MANUSCRIPT

Nomenclature velocity along the axis of the pipe [m/s] isentropic sound wave speed [m/s] spatial coordinate [m]

v

cross section area of the pipe[m2]

Cp D

heat capacity at constant Vw pressure [J/kg K] x inner diameter of the pipe [m]

fr

friction factor [-]

Z

g

standard gravity [m/s2]

γ =

m&

mass flow rate [kg/s]

η id

p

pressure [Pa]

R

specific gas constant [J/kg K]

H

isentropic head [kJ/kg]

Nr

Compressor speed (rpm)

Nr,max

Maximum [rpm]

compressor

Nr,min

Minimum [rpm]

compressor

Qac

Actual flow rate [m3/s]

Subscripts

a1 − a6

coefficients for centrifugal compressor map [-]

s

suction side of a compressor

d

discharge side of a compressor

Cp Cv

M AN U

τw

compressibility factor [-]

ρ

speed θ l

TE D

speed

ratio of constant pressure and constant volume heat capacities thermodynamic efficiency of compressors heat flow into the pipe per unit length of pipe and per unit time [J/m s] shear stress at the wall [kg/m s2]

SC

Ω

RI PT

Ac

Ν( µ ,
 )

density of the gas [kg/m3] angle of inclination of pipe to the horizontal [radian] Normal distribution with mean µ and standard deviation


EP

Dynamic simulation of natural gas transmission networks (GTNs) includes several aspects such as transient analysis [6], pack manipulation [7] and discontinuities [4], and

AC C

a summary of the issue is discussed in Section 2. Natural gas demand forecasting is an aiding tool in future development, planning, operation, optimization and control of gas networks. This can ensure us in meeting the contractual requirements of customer demand flow and minimum pressure level. Also, there are several other reasons that motivate natural gas (transmission or distribution) companies toward developing more 4

ACCEPTED MANUSCRIPT

efficient and accurate methodologies of forecasting natural gas consumption. Among these are: Improving profitability, flexibility and reliability

•

Managing capacity and planning for peak-shaving and gas storage facilities

•

Reducing the possibility of sudden switches between different modes of operation Enhancing the efficiency of daily operation

•

Maintenance planning and scheduling

M AN U

•

SC

RI PT

•

Natural gas networks [8] are studied in the literature from various points of view including modeling [9], (steady and dynamic) simulation [10], centralized [11] and decentralized [12] state estimation, steady optimization [13] and dynamic optimization

TE D

[14], however, uncertainty was always ignored. The counterpart of the optimization problems encountered in GTNs can be found in electrical engineering as power optimization problems [15].

EP

One of the early works in the area of dynamic optimization of gas networks is the work of Marqués and Morari [16] in which a two-layer hierarchical optimization algorithm was

AC C

proposed. The lower layer sets the optimal discharge pressure of compressor stations near to off-takes in a pre-specified time intervals during a day. On the other hand, the upper layer calculates the daily set points of the rest of the compressor stations. An estimator is utilized to compensate for the demand forecast errors. Aalto et al. [17] used data obtained from a pipeline simulator to identify linear models of Finnish natural gas 5

ACCEPTED MANUSCRIPT

pipeline. Then, on the basis of these identified linear models, a constrained quadratic programming optimization problem was established using linearized constraints and quadratic approximation of energy consumption of the compressor stations as the

RI PT

objective function. The proposed framework was evaluated for the Finnish natural gas pipeline, where 5 to 8% savings in compressor energy consumption was achieved. Abbaspour et al. [5] proposed a transient optimization framework for the minimization of

SC

the total fuel consumption in the daily operation of gas networks. The line-pack profile to be maintained and the gas withdrawal profile of the customers are assumed to be

M AN U

known. Then, the rotational speed of the compressors is determined in such a way that the desired line-pack profile is maintained while operational constraints are respected. The effect of the number of subintervals on discretization of simulation time was also studied. Baumrucker and Biegler [18] formulated the optimal control of gas pipelines

TE D

systems as a mathematical programming problem with equilibrium constraints to handle non-smooth issues encountered such as flow reversals. The work by Carter and Rachford [19] is the earliest work to consider load uncertainties

EP

in gas network operation planning. In their work, different scenarios were assumed for

AC C

the future loads, and also a decision point representing the time at which a power plant (as a customer) will actually come online is determined. Finally, different strategies are provided by the optimizer for each scenario, which are common up to the decision point, and after this point, an appropriate strategy corresponding to each scenario is suggested. They do not provide the mathematical background of their work; however, the conceptual analysis of the results especially line-pack inventory variations in the 6

ACCEPTED MANUSCRIPT

system is carried out in detail. Zavala [20] proposed a scenario-based stochastic optimal control model for the optimization of gas networks isothermal operation under load uncertainties, where a weighted risk-mean objective function is considered.

RI PT

Increased robustness in operation of the system is demonstrated in a case study. The applications of expert systems in optimization of natural gas pipeline operations is studied in [21]. Nguyen and Chan [22] discuss the applications of evolutionary

SC

scheduling in compressor scheduling.

M AN U

By reviewing the previous studies, it is noted that there are a few studies about the load uncertainties in gas networks planning, where all are on the basis of characterizing the uncertainties using multiple scenarios.

The present paper addresses the daily planning problem for natural gas transmission systems under demand uncertainties. The appearance of these uncertain disturbances

TE D

in the dynamic model of the plant makes the objective function and the non-equality constraints uncertain as well. If the probability density function (PDF) of the uncertain

EP

parameters are known, accurate calculations of the PDF corresponding to the objective function and other dependent variables require the calculation of multidimensional

AC C

integrals [23] mostly using numerical methods where an efficient method such as orthogonal collocation (on finite elements) method can be used. However, as the number of

the uncertain parameters is increased, these methods become

computationally intractable and curse of the dimensionality is the main drawback of these methods. Long run times of the simulations in the case of the dynamic problems, make the aforementioned methods less efficient. To overcome these issues, mean7

ACCEPTED MANUSCRIPT

variance optimization models [24] are commonly adopted where the calculation of the expectation of the objective function and feasibility of the non-equality constraints are

RI PT

the main concerns. Considering PDF-based characterization of random effects [25] for general nonlinear models and considering the fact that normal or Gaussian distribution is believed in process engineering practice to be a sufficient assumption [23], some uncertain

SC

parameters with known mean and variance are used to characterize the user demand.

M AN U

However, the PDF of the dependent variables of the process model might not be normal in general. In these cases, Monte Carlo methods, power series expansion (PSE) based uncertainty propagation methods [26] and polynomial chaos expansion (PCE) method [27] can be used. PCE is a probabilistic method that uses orthogonal stochastic polynomials in the random inputs in order to project model outputs.

TE D

Chaffart et al. [28] applied the PSE method in the analysis of the parametric uncertainty of a heterogeneous multi-scale catalytic reactor model. Pintarič et al. [29] used

EP

sensitivity analyses to evaluate the influence of uncertain parameters on the first-stage design variables and the objective function in order to reduce the scenarios for

AC C

deterministic equivalent problem. Ostrovsky et al. [30] developed an approximate transformation to convert chance constraints into deterministic ones to avoid calculation of multiple integrals. Rasoulian and Ricardez-Sandoval [31] applied PSE method to quantify the effect of parametric uncertainty in an epitaxial thin film growth process. Nagy and Braatz [32] evaluated the efficiency and accuracy of PSE and PCE in a simulated batch cooling crystallization process. 8

ACCEPTED MANUSCRIPT

Unscented transform is an alternative to these methods which requires a small number of sample points compared to its computationally expensive counterpart. This makes the unscented transform a powerful tool that can be used to propagate a random

RI PT

variable through a nonlinear function [33]. The application of this transform is demonstrated in the closed-loop set point optimization of extractive distillation processes [34] and pressure swing distillation processes [35] under feed composition

SC

disturbances.

M AN U

Some features of this algorithm make it more suitable for dynamic optimization of natural gas transmission network under demand uncertainties: 1. This algorithm is suitable for the cases with black-box models as the calculation of the outputs are required for a finite number of sigma points for any given transformation.

TE D

2. Considerable reduction of computational load can be achieved. For example, if n is the number of uncertain parameters, uncertainty propagation on the basis of

EP

the evaluation of the multidimensional integrals using Gaussian quadrature [36], a third order PSE approximation [26] with forth-order accuracy and unscented

AC C

transform requires 3n, 6n and 2n+1 grid points, respectively. 3. The algorithm can be used with discontinuous transformations. The linearized estimates of a discontinuous function cannot capture the information. Discontinuities can occur in GTNs in situations such as valve position changes from open to fully closed, compressor startup or shut-down and flow reversals in the network as well [4]. 9

ACCEPTED MANUSCRIPT

4. Derivative information are not required. However, as shown in the case study, the objective function and the non-equality

RI PT

constraints of the problem under study can be well approximated using Gaussian PDFs which motivated the application of the unscented transform as the uncertainty propagation tool.

SC

Accordingly, the unscented transform (UT) [37] is used as a means for the propagation of the mean and covariance of the uncertain parameters through nonlinear

M AN U

transformations (i.e., objective function and constraints). Then, an optimization problem on the basis of the non-isothermal model of GTN is formulated to which its solution provides the daily planning, mainly discharge pressure set-points of the compressor stations. These set-points enables us to manage the line-pack of the pipes in order to

TE D

achieve some pre-specified goals with minimum cost while the operational constraints are respected.

Section 2 describes the dynamic model of a GTN. Section 3 discussed various aspects

EP

of natural gas demand of a customer, and uncertainty related to the forecasted demand is the topic of Section 4. Formulation of the optimization algorithm is described in

AC C

Section5. A benchmark is studied in Section Error! Reference source not found. to evaluate the performances of the stochastic optimization approach, and finally conclusions are drawn. 2.

Model

10

ACCEPTED MANUSCRIPT

Various equipment such as compressors, valves, and pressure regulators are put together to form a natural gas transmission network. For the purpose of dynamic simulation in order to obtain the temporal behavior of the states of such a network (i.e.,

RI PT

pressure, mass flow rate and temperature), its mathematical model can be derived on the basis of the model of each equipment together with the appropriate boundary conditions. Accordingly, the dynamic model of a GTN includes:

a set of nonlinear partial differential equations (PDEs) including mass,

SC

•

M AN U

momentum and energy balance equations [38] corresponding to the pipes (Section 2.1) •

equations required to evaluate thermo-physical properties of natural gas

•

pressure-flow and energy balance equations of a compressor (Section 2.2)

•

equations corresponding to the conservation of mass and energy at junction

•

TE D

nodes

the conditions specified at the network boundaries

EP

The governing equations corresponding to a pipe segment, a compressor and a valve will be discussed in the subsequent sections. Furthermore, various aspects of solution

AC C

method of this system of ODEs can be found in [4]. 2.1

Governing equation for gas pipeline

The model representing the one dimensional transient behavior (i.e., pressure, mass flow rate and temperature) of the gas flowing in a pipeline consists of the following governing equations [4]: 11

ACCEPTED MANUSCRIPT

∂p ZRT ∂m& =− ∂t Ac ∂x

(1)

&m & Ac pg & & 2 ZRT  ∂m ∂ m ∂p  fr ZRT  m =−  − − − A sin (θl )  c   ∂t ∂x  pAc  ∂x  2DAc  p ZRT

(2)

where

τw

2  Vw +   Cp

2 p ∂Z  Ω + τ wπ Dv  T ∂Z   ∂v  Vw  + − 1 ( ) 1− ( )T p    Z ∂T   ∂x  C p p  Z ∂p  Ac 

(3)

RI PT

∂T  ∂T = v ∂t  ∂x

is the shear stress at the wall and can be obtained by Darcy-Weisbach

SC

equation (Eq.(4)) on the basis of friction factor [5] and the parameter Vw can be

Vw =

fr ρv v 8 ZRT  p ∂Z p )T − 1 − ( ρ C pT  Z ∂p

2  T ∂Z   1 + Z ( ∂T )p   

(4)

(5)

TE D

τw =

M AN U

calculated by Eq. (5):

These PDEs are converted to ODEs by applying the method of lines technique along with the orthogonal collocation method used for spatial discretization [39]. Centrifugal compressor model

EP

2.2

Integration of a centrifugal compressor model and its performance map enables us to

AC C

build a more realistic simulation model of a GTN. Multiple compressors might be grouped together (as in series or parallel connections) to form compressor stations (CSs). In practice, compressor stations are operated under discharge pressure control [19], where setting the discharge pressure, completely determine other parameters including suction pressure, gas flow, rotational speed and adiabatic efficiency. Choosing station discharge pressure as the control input has the added benefit of separating 12

ACCEPTED MANUSCRIPT

compressor internal model from the network model [40]. The behavior of a centrifugal compressor is expressed by two different types of characteristic curves: Curves of head against flow rate for various fixed values of compressor rotational

RI PT

i.

speed (rpm). ii.

Curves of head against flow rate for various constant adiabatic efficiency lines.

These two curves together with the surge and stonewall limits completely describe a

SC

centrifugal compressor (i.e., compressor wheel map [41]). For simulation purposes, a 2

M AN U

wheel map can be described by three normalized parameters, where H Nr and

ηid

are

considered to be a quadratic polynomial function of Qac Nr (Eqs. (6) and (7)), and constant coefficients ai can be obtained by the method of least squares [5]. 2

Q  Q  ηid = a4 + a5  ac  + a5  ac   Nr   Nr 

2

TE D

Q  Q  H = a1 + a2  ac  + a3  ac  2 Nr  Nr   Nr 

(6) (7)

EP

Eqs. (8) and (9) (i.e., adiabatic power and discharge temperature) together with Eqs. (6) and (7) can be used to model the quasi-steady behavior of a centrifugal compressor as

AC C

the dynamics of a CS is much faster than the pipeline dynamics [17].  γ −1      γ Ts  pd  γ   ϒ = 47.4 − 1   γ − 1   ps   

(8)

 γ −1      Ts  pd  γ   Td = Ts + − 1   ηid  ps   

(9)

13

ACCEPTED MANUSCRIPT

A typical centrifugal compressor wheel map is shown in Fig. 1. The rotational speed of the compressor and constant efficiency lines parameterizes the compressor wheel map where constant efficiency lines are drawn with dashed lines. Considering a compressor

RI PT

model in the simulation of a gas transmission system, imposes five constraints in the optimization studies. Surge line, choke line, maximum allowable speed, minimum required speed are the four constraints imposed by the compressor wheel map, where

SC

the operating point of a compressor should remain inside this envelope without any violation of the limits (Fig. 1). Minimum suction pressure is the other constraint required

M AN U

from the operational point of view. Inclusion of a compressor wheel map in the simulation model of a GTN introduces further nonlinearity into the model, which adds another dimension of complexity. Also, the small amount of natural gas used as fuel in compressor stations is ignored in our simulations. Pressure buffer

TE D

2.3

Gas transmission companies must provide the minimum pressure required by the customers. The pressure buffer is defined as the difference between the pressure of the

EP

gas delivered to the customer and the contractually guaranteed delivery pressure. Fig. 2 shows a typical display of actual pressure at a delivery point and the minimum pressure

AC C

obligation of the customer during the operation of the transmission system. According to Fig. 2, pressure buffer has a transient behavior called transient pressure buffer which can be regarded as a safety factor used by the dispatchers in the operation of the system. This parameter is important as its higher values guarantee the satisfaction of the customer at the expense of higher operating cost of the compression units. On the 14

ACCEPTED MANUSCRIPT

other hand, a smaller pressure buffer means less energy consumption in compression units, causing an increased possibility of failure in the satisfaction of the customer demands considering unpredicted disturbances. The application of the concept of

RI PT

transient pressure buffer and its optimal calculation will be discussed further in Section 6. A similar concept called comfort zone is defined and applied by Nguyen et al. [42].

SC

3. Demand forecasting

Customers’ consumption behavior (also known as consumption pattern) should be

M AN U

taken into account in order to be able to provide the customers with required amount and minimum pressure while maintaining the operation of the transmission system reliable and economical. The pattern of consumption is also important in the development of a forecasting procedure [43] as there are different methods for different

TE D

types of consumers. For example, domestic customers have large variations in their demands in cold weathers during a day as they are mostly influenced by outdoor temperature, so they are called temperature-sensitive customers [1]. Alternatively, gas-

EP

fired power plants might have unnoticed large loads as they might become operational to compensate a sudden increase in electricity demand in summer [44]. Types of natural gas consumers

AC C

3.1

As stated by Vitullo [1] and American Gas Association [45], natural gas consumers can be categorized into five main types: residential customers, commercial customers, industrial customers, power generators and natural gas vehicles. A brief summary of key features of the demand by these customers can be found in Table 1, where each 15

ACCEPTED MANUSCRIPT

customer has its unique demand characteristics [46]. Despite these vast differences in characteristics, in this paper, the study is carried out only for the gas-fired power plants due to their unnoticed nature [44] to illustrate the application of the proposed method for

RI PT

considering demand uncertainty in natural gas network planning.

Apart from the type of the customer, the duration of prediction is an important issue in the process of demand forecasting as well, which might have a wide range including

SC

next hour, next month, next year or even more [47]. Forecasting problem is a

M AN U

challenging issue as many uncertainties can arise, and the duration of prediction has a significant effect on the accuracy of the forecasting. An extremely broad range of timescales is one of the aspects of the forecasting problem, where a short, medium or long range might be considered in the development of the forecasting tool. Most of the developed models in this regards are for medium or long term periods. Planning for in new

infrastructural

development

TE D

investments

or

other

strategic

planning,

generally requires a long term future demand forecast. On the other hand, in the development of daily strategies for the optimal operation of the system, a day ahead

EP

prediction horizon is sufficient. Here, a model capable of making hourly gas load

AC C

forecasts is required as we are interested in the optimal operation of a GTN on a daily basis (called a gas day). 3.2

Forecasting Methods

As the development of a model for forecasting is strongly depended on the pattern of the consumption by the customer and the prediction horizon, different models have been introduced. Different approaches have been applied in the load forecasting 16

ACCEPTED MANUSCRIPT

problem in natural gas industry, among which the most common mathematical modeling techniques are the ones with the statistical background such as autoregressive moving average models and the ones with artificial intelligence background especially neural

RI PT

and fuzzy neural network models [47].

One of the main key issues in the successful application of statistical methods is the availability of a rich collection of historical data. On the other hand, black box nature is

SC

the main disadvantage of artificial intelligence based models due to their non-intuitive

M AN U

nature which makes their explanation to non-experts a challenging task [48]. Also, a review of various methods concerning prediction of natural gas consumption [46] and relevant discussions such as short-term load forecasting [43], industrial end-users [49] and artificial intelligence-based methods [48] can be found in [50] and references therein.

Customer daily load model characterization

3.3.1 Power plants

TE D

3.3

Compensation of large loads form large end users, is a challenging task for distribution

EP

companies [2]. Loads from a gas-fired power plants or combined cycle plants fueled

AC C

with natural gas, as generally occur with little or no advance warning [19], can introduce uncertainties in the operation of gas distribution systems. Due to their unnoticed nature, the onset time, load duration and even the amount of the required gas for this type of customers have some levels of uncertainty. These parameters and the uncertainties associated with them can be estimated on the basis of the past experiences. Planning on the basis of the average demand condition for the future might result in catastrophic 17

ACCEPTED MANUSCRIPT

results. As a power plant comes online, a large load is introduced to the network and commonly a prior preparation is necessary to meet this unnoticed large load, and since the onset time is uncertain, the available time for preparation is uncertain as well. Here,

characterized as follow. The assumptions are:

RI PT

the load pattern of gas-fired power plant and the corresponding uncertainties are

1. Whether a particular power plant will become online or not is known in advance.

SC

2. The load of a power plant which will become online is predictable, where it is

M AN U

possible to characterize its load with a rectangular pulse-shaped signal as shown in Fig. 3 with three parameters, where (ts), (td≠0) and (Μ) are the onset time, duration and load amount, respectively.

3. Parameters (ts), (td) and (Μ) are independent Gaussian random variables with

TE D

known mean and covariance.

Accordingly, for a deterministic scenario we know that at time t = (ts) hr., the power plant will become online requiring (Μ) MMSCFD for (td) hours. On the other hand, for

EP

uncertain scenarios, (ts), (td) and (Μ) are assumed to be Gaussian random variables, whose mean and covariance are determined for a specific plant on the basis of the past

AC C

experiences and also time of the year, weather condition and etc. [19]. 4. Uncertainty

Due to several reasons, the most accurate predictions achieved will have some level of uncertainty [2], and appropriate approach to deal with uncertainty in operational planning must be adopted. For example, electricity demand is affected by the day of the 18

ACCEPTED MANUSCRIPT

week, seasonal variations, holiday periods, feast days and the weather conditions [51]. Whether conditions can also be predicted with uncertainty. Accordingly, planning of a

RI PT

gas-fired power plant will suffer from the weather forecast uncertainties as well. Due to the uncertainty in the forecasted demands, a method that provides hedging against demand uncertainty is essential. If the expected natural gas consumption is considered as the basis for the planning of a gas day, there might be under or over-

SC

forecasting. Under-forecast might result in catastrophic outcomes such as inefficiency in

M AN U

meeting peak demands where in the case of industrial customer might results in a complete shutdown of the plant or in the case of a power generation unit could result in black-outs. Over-forecast will have significant economic consequences. Additionally, beside the development of accurate methods for the estimation of the best forecast, the uncertainty levels corresponding to these forecasts should be quantified.

TE D

Then, this uncertainty can be used in the calculation of the risk of the decisions made on the basis of these forecasts.

Dealing with uncertainty

EP

4.1

The uncertain parameters (i.e., customers’ demands) can be described by possible

AC C

scenarios [52] or probability distributions. These parameters appear in the objective function and constraints defined in the optimization problem arising in daily planning of a GTN, as well. Accordingly, stochastic programming is applied to deal with uncertain parameters in the formulation of the optimization problem.

19

ACCEPTED MANUSCRIPT

In the scenario-based approach, several simultaneous future load scenarios are considered, however there is an uncertainty about which of several situations will actually occur. This approach has been studied for the planning of GTNs previously by

RI PT

Carter and Rachford [19] and Zavala [20], and one drawback of scenario-based approach is that a decision time is required as input to the optimization program. This parameter shows the time when we expect to know which of the assumed scenarios will

SC

actually occur. However, in practice, this parameter might be uncertain, as well.

M AN U

Here, the load by a customer is characterized by some parameters with a known probability distribution. Then, the method of the chance-constrained programming is applied in which each constraint is satisfied at a predetermined confidence level specified by the decision maker. In this approach, on the basis of the previously observed behavior of the customer, we can increase or decrease the level of

TE D

uncertainty in its characterizing parameters. For example, a power plant onset time might be somewhat invariable, however, its requested load shows more fluctuations. Applications of stochastic programming approach

EP

4.2

Short-term storage of natural gas in transmission pipelines (called line-packing) and

AC C

facilities designed to provide additional natural gas on very short notice (called peakshaving facilities) are the two viable alternatives available to a gas distribution company to meet the customer loads beyond the normal requirements. Packing natural gas is the easiest way to enhance the storage capacity of the pipeline network which provides the gas supplier with an effective method to meet customer demands in peak conditions. However, this strategy is applicable up to a point when a decreased accumulation 20

ACCEPTED MANUSCRIPT

capability of the network is encountered. Beyond that point, it might still be possible to enhance the capacity of the transmission system by means of secondary supply

RI PT

sources (i.e., peak shaving facilities). The main objective is to operate economically while simultaneously maintaining the existing level of service. In order to achieve these aims, the most important prerequisite is the prediction of customer consumption pattern as well as the peak demands. Then,

SC

this information can be utilized in the development of the appropriate strategies aimed

M AN U

at reducing the risk of meeting the demands during peak periods. 5. Stochastic dynamic optimization problem formulation

The general formulation of a stochastic dynamic optimization problem for GTN planning has the form: J ( x ( t ) ,u ( t ) ,θ ) min ()

Subject to:

f ( xad ( t ) ,u ( t ) ,θ ) = 0

TE D

ad

u t

EP

g ( xad ( t ) ,u ( t ) ,θ ) ≤ 0

(10)

(11) (12)

AC C

In this formulation, xad represents differential and algebraic state variables, u ( t ) and

y( t ) are the vectors of manipulated and output variables, and θ is the vector of uncertain parameters used in the representation of disturbances. The process dynamics is described by a set of differential-algebraic equations (Eq. (11)), and Eq. (12) is used to impose path inequality constraints on the state and control variables. Due to the large 21

ACCEPTED MANUSCRIPT

number of state variables and special attentions required in the solution of the initial value problems in GTN dynamic simulation [4], the control parameterization method is

RI PT

applied [53]. The chance constrained programming technique can be utilized to convert the stochastic problem stated in Eqs. (10)-(12) into an equivalent deterministic nonlinear

SC

programming problem as follows.

The deterministic equivalent of the stochastic objective function (10) is formulated using

M AN U

a weighted objective function as: J ( x ( t ) , u ( t ) ,θ ) = min µ E [ J ] + ( 1 − µ ) min () () ad

u t

u t

Var [ J ]

(13)

Where µ ≥ 0 indicate the relative importance of the mean and standard deviation of stochastic objective function (10) in minimization. Inequality path-constraints are

tf

hi = ∫ ( ri ( t ) .ri ( t ) dt ) −ε ≤ 0 where

(14)

EP

t0

TE D

handled within the control parameterization framework as follows:

ε is a small positive number and ri ( t ) corresponds to the ith member of the ng

and measures the degree of the constraint violation:

AC C

g (t ) ∈ R

ri ( t ) = max ( gi ( t ) ,0 )

(15)

As random parameters are present in constraints, theses constraints are probabilistic, and has to be satisfied with a pre-specified probability of pj , in other words:

22

ACCEPTED MANUSCRIPT

P [ hi ≤ 0 ] ≥ p j

(16)

which means that the probability of realizing hi less than or equal to zero must be

RI PT

greater than or equal to pj . Knowing the mean and covariance of hi , constraint (16) can be restated as: E [ hi ] + φ ( p j ) Var [ hi ] ≤ 0

φ

SC

where

(17)

is the distribution function of the standard normal variable. This formulation

satisfied with a certain probability. 5.1

The unscented transform

M AN U

guarantees that for example, the customer contracted minimum pressure will be

The transformation of uncertainty is implemented using the unscented transform. Let

σ θ2 ∈ R nθ ×nθ

TE D

θ ∈ R nθ be a vector of Gaussian random variables whose mean  ∈   and covariance are assumed to be known. The random vector

y ∈R y is obtained from the n

transformation y = f (θ ) , where f : R θ → R y is a nonlinear and (in our case also n

EP

n

discontinuous) map. The advantages of using UT for discontinuous functions is

AC C

discussed in [54]. UT can be used to propagate the pdf of the system random parameters (i.e., θ ), and if a Gaussian approximation   ,
  is assumed to represent the stochastic variable y , it is accurate up to second order in estimating mean and covariance [37]. Now, to approximate the mean µy = E [ y] and covariance
 = [] of y , the following steps are involved: 23

ACCEPTED MANUSCRIPT

1) 2nθ + 1 sigma-points are calculated on the basis of the mean and covariance of

θ :  µθ , i = 0  θ( i ) =  µθ + γ Si , i = 1,K , nθ  µ − γ S , i = n + 1,K ,2 n  θ i θ θ

RI PT

(18)

Where Si is the i -th column of the matrix square root of
 and γ is a scaling parameter [55],

SC

γ = nθ + λ , λ = α 2 ( nθ + κ ) − nθ ,

(19)

M AN U

2) Each of the 2nθ + 1 sigma points is propagated through the function f to form

y(i ) :

( )

y(i) = f θ(i) , y(i) ∈R y n

(20)

3) A weighted average of the transformed points is used to predict the mean: 2 nθ

i =0

(21)

TE D

E [ y ] = ∑ wim y(i )

Moreover, the covariance is approximated using the weighted outer product of the

2 nθ

(

EP

transformed points:

)(

Cov [ y] = ∑ wic y(i ) − E [ y] y(i ) − E [ y]

T

(22)

AC C

i =0

)

4) Each sigma point is associated with a weight, where the weights are defined as:

λ λ  m c 2  wi = n + λ , wi = n + λ + ( 1 − α + β ) , i = 0  θ θ  1  wim = wic = , i = 1,K ,2nθ  2 ( nθ + λ )

24

(23)

ACCEPTED MANUSCRIPT

α, β and κ are the design parameters of UT. The semi-positive definiteness of the covariance matrix can be guaranteed choosing κ ≥ 0, which is usually set to 0 [56]. The parameter α (0 ≤ α ≤ 1) determines the spread of sigma points around  [55]. β is a

RI PT

non-negative parameter which is used to incorporate our prior knowledge about the distribution of µy and for a Gaussian distribution the optimal choice is β = 2 [55]. For

SC

UT, α, β and κ are set to the following values: α=1, β=2 and κ=0. In this work, nθ is the total number of the parameters used in the representation of the uncertain customer

M AN U

loads. It is also worth mentioning that the tuning parameters of the unscented transform (i.e., α, β and κ) can be used to handle the cases with non-normal distribution if the skew of the distribution is nonzero. 6. Case study Problem statement

TE D

6.1

A transmission line including one supply (S1), one end delivery point (D2) and one side delivery (D1) is considered as the case study. Pipes P1, P3 and P5 are each 100 km

EP

long and Pipes P2, P4 and P6 are each 40 km long. Each of six pipes has a diameter of 1.422 m and roughness of 17.78 micrometers. Compressor stations are located every

AC C

140 km along the transmission line, where each station has three identical parallel centrifugal compressors. A constant ambient temperature of 10°C and overall heattransfer coefficient of 4.4 W/(m.K) are specified at all points in the network. A schematic of this transmission line can be seen in Fig. 4.

25

ACCEPTED MANUSCRIPT

The supply pressure at node N1 is maintained at 1080 psia (7446.3 kPa) and natural gas at the constant temperature of 50 °C with a spe cific gravity of 0.6 is provided at this node. The pipeline delivers natural gas to the end point, where the delivery to the end

RI PT

user at node N9 must be above a certain minimum pressure (i.e., 700 psia = 4826.3 kPa). Along the pipeline, a location is considered at which gas can be delivered to a natural gas-fired power plant. The demand of this customer is not known a priori and

SC

uncertain scenario is assumed for its requested demand. This side delivery is included to investigate the effect of side load on the operational planning of the transmission line

6.2

M AN U

under study.

Power plant side customer

It is assumed that the side customer is a power plant. For this case, it is also assumed that the flow must be delivered continuously at the rate of 3100 MMSCFD (87.8

TE D

MMSCMD) to the end user. Initially, side delivery is assumed to be zero, and steady state optimization is applied to obtain the optimal steady state operating condition of the compressor stations, which is shown in Fig. 4. This is considered as the initial condition

EP

for the simulations. As can be seen in this figure, the pressure at node N9 is at its

AC C

minimum allowed value.

6.2.1 Power plant side customer at D1 The following scenarios for power plant load are imposed on the side delivery D1 to investigate the effectiveness of the proposed stochastic approach. The request time by the power plant is uncertain, which means that the time available for preparation is uncertain as well. In this example, it is assumed that the power plant will definitely come 26

ACCEPTED MANUSCRIPT

online and the uncertain parameters characterizing its load (i.e., load onset time, duration and rate) are shown in Table 2. These parameters can be determined on the basis of past experiences and operational data of the plant. Accordingly, the power

RI PT

plant comes online at time Ν (9.00, 0.3872) hours and it draws Ν(500.0, 502) MMSCFD for Ν(4.00, 0.3872) hours. Also, four different values for constraint satisfaction

SC

probability ( pj ) are assumed.

In order to investigate the effectiveness of the stochastic approach, other than the

M AN U

uncertain load scenarios with stochastic parameters, three deterministic scenarios for load pattern are defined. The first one is an average (or expected) condition scenario called “Expected” in which power plant comes online at time 9 hours for 4 hours requiring 500 MMSCFD of natural gas. Also, two other scenarios on the basis of Monte

TE D

Carlo random simulations are considered.

1e5 random simulations are run, and the load parameters corresponding to the minimum and maximum total amount of natural gas request realizations, called “Min”

EP

and ”Max” respectively, are chosen. Accordingly, as shown in Table 2, it is known that at 7.95 hours and 10.14 hours the power plant will come online requiring 660.5 and

AC C

385.6 MMSCFD of natural gas for 5.69 and 2.18 hours for maximum and minimum load condition, respectively.

If any of the loads in Table 2 is imposed on the system and the compressor discharge pressures are kept constant at their steady values, the end point delivery pressure will fall below 700 psia. Accordingly, the requirement of advanced preparation is the 27

ACCEPTED MANUSCRIPT

common feature of these different scenarios, however, preparation for each scenario requires unique set point trajectories for compressors discharge pressure.

RI PT

The main objectives are the satisfaction of the end and side customers under uncertain peaking load without any interruption while the power consumption of the compressors is minimized. The objective function of the optimization problem can be expressed as

J =∑ i

t f =24

∫

SC

the total power consumption of all compression units:

ϒi ( t ) dt

t0 =0

(24)

M AN U

where ϒi ( t ) denotes the power consumption of compression unit i at time t . The discharge pressure set point of the stations discretized on ½-hour intervals during a day are chosen as the optimization decision variables. The simulation model of the GTN is implemented in C++ environment [4] and all the optimization problems are solved using

TE D

the optimization package NOMAD [57] that implements the “Mesh Adaptive Direct Search” algorithm for optimization under general nonlinear constraints. This algorithm has extended the “Generalized Pattern Search” algorithm [58] using local exploration.

EP

The main feature of this algorithm that makes it more appealing for solving the

AC C

stochastic dynamic optimization problems is its capability for solving black-box optimization problems where the simulation model has the following characteristics: 1. Simulations are time-consuming 2. Simulation might fail to return a valid output even for feasible points 3. No analytic derivative information is available about the simulations.

28

ACCEPTED MANUSCRIPT

The constraints to be respected are: 1. The supply pressure is maintained at the constant value of 1080 psia.

RI PT

2. End customer satisfaction which means the delivered pressure at N9 must exceed 700 psia and flow rate must be 3100 MMSCFD.

3. The side customer (i.e., power plant) satisfaction which means the delivered pressure at N8 must exceed 600 psia and the requested flow rate must be

SC

provided without any interruption.

M AN U

4. Operating constraints of compressors (maximum available power, minimum suction pressure, maximum discharge temperature and wheel map) must be respected.

5. Maximum operating pressure of the line is set to 1200 psia. It is also assumed that preparation for upcoming scenarios is possible as soon as the

TE D

gas day starts. It is also necessary to prevent the depletion of the system at the end of the gas day as it might introduce difficulties in the operation of the next day. This goal is accomplished by imposing a non-equality constraint stated as the total line-pack of the

EP

system at the end of the gas day must be greater than or equal to the initial line-pack of the system. Although different end constraints can be considered in the optimization

AC C

problem, such as the initial condition of the current day or the initial condition of the next day, the aforementioned constraint significantly reduces the number of end constraints. No constraint is imposed on the inlet flow rate at supply node during the optimization runs and it is implicitly controlled by the downstream conditions. It is worth mentioning that, if the obtained inlet flow trajectory violates the production capacity of the gas 29

ACCEPTED MANUSCRIPT

production facility, peak-shaving facilities might be planned to be utilized at the appropriate moments.

RI PT

6.2.2 Optimization results and discussion The optimization problem was solved considering different scenarios as shown in Table 2. The run time of each deterministic and stochastic optimization problem on an Intel Core i7-4770 PC, 3.40 GHz and 8 GB RAM is about 15 minutes and one hour,

SC

respectively. First, the assumption of the normal or Gaussian distribution of the

M AN U

objective function and constraints should be validated.

As the values of the decision variables are changing during the optimization iterations, the number of active constraints are changing as well. For example, considering the wheel map of a compressor, the operating points of the compressors which are determined by their actual flow rates and isentropic head characteristics, should remain

TE D

inside the envelope which can be stated using the following three constraints. If surge and stonewall limits which are the minimum and maximum stable compressor

EP

flow points, respectively, are approximated with first order polynomials correlating the actual flow rate with isentropic head at the surge and chocking points of the compressor

AC C

as ( H surge = b1 + b2Qac ) and ( Hstonewall = b3 + b4Qac ) , respectively, the wheel map constraint can be mathematically stated as [33]: 2   Qac   Qac   N a1 + a2   + a3    ≤ b1 + b2Qac 424 3  Nr   Nr   1  HSurge  1444442444443 2 r

(25)

H From Eq.( 6)

30

ACCEPTED MANUSCRIPT

2   Qac   Qac   N a1 + a2   + a3    ≥ b3 + b4Qac 424 3  Nr   Nr   1  HStonewall  1444442444443

(26)

Nr,min ≤ Nr ≤ Nr,max

(27)

2 r

RI PT

H From Eq.( 6)

It is clear that all the constraints (25) to (27) cannot simultaneously be active. On the

SC

basis of the formulation of the stochastic dynamic optimization, for example, the

tf

(

)

(

)

h1 = ∫ max( Nr ( t ) − Nr ,max ,0 ) dt −ε ≤ 0 2

t0 tf

M AN U

constraint (27) is decomposed into two different constraints as:

h2 = ∫ max( Nr ,min − Nr ( t ) ,0 ) dt −ε ≤ 0 2

t0

(28) (29)

TE D

The distributions of the objective function and constraint corresponding to the delivery pressure at node N9, due to the imposed demand uncertainties are obtained through 1e5 Monte-Carlo simulations and are depicted as a histogram in Fig. 5 and Fig. 6 for the

EP

first optimization iteration which can be considered as normally distributed. It can been seen that the objective function is a normal random variable with the

AC C

distribution of Ν(1955.6, 48.552) and Ν(1955.5, 48.132) MWh for Monte Carlo simulations and unscented transform, respectively. Also, the constraint corresponding to the end point delivery pressure is a normal random variable with the distribution of Ν(2.185e5, 2.21e42) psi2.hr and Ν(2.18e5, 2.18e42) psi2.hr for Monte Carlo simulations and unscented transform, respectively. Comparison of the mean and variance of the 31

ACCEPTED MANUSCRIPT

distributions for the unscented transform and Monte-Carlo simulations shows satisfactory agreements. While the mean and variance of the all active constraints and objective function change during the optimization runs, the assumption of the normality

RI PT

can be validated through Monte-Carlo simulations at each iteration and final optimization solution as well.

The results of stochastic scenario “St. 1” are compared to those with deterministic

SC

scenarios. In order to validate the results of the stochastic formulation, for the optimum

M AN U

values of the compressor discharge pressures obtained using stochastic formulation, 1e5 Monte-Carlo simulations are run and 1015 out of the all realized simulations (i.e., approximately 1% of realizations) violated the imposed constraint corresponding to the delivery pressure at node N9.

However, the trajectory of the operating points of the compressors have enough

TE D

distance from the boundaries of the wheel map and all the constraints corresponding to the wheel map boundaries are inactive for all the Monte Carlo realizations. This is due

EP

to the fact that while the optimizer is trying to keep the operating points away from the boundaries of the wheel map to be feasible, at the same time it tries to optimize the total

AC C

energy consumption of the compressors which forces the operating point toward the high efficiency regions which are located inside the wheel map. Figs. 7-9 show the optimal line-pack trajectory relative to the initial steady state condition for pipes 1-6, respectively, in preparation for different scenarios. The line-pack trajectories of pipes 1-4 roughly follow the same trend, where initially their gas content is depleted, then gas is packed to restore them to their initial condition. On the other hand, 32

ACCEPTED MANUSCRIPT

the line-pack trajectories of pipes 5-6 show the opposite trend of the rest of the pipes, where initially gas is packed to make them ready for the upcoming load, and then their packed gas is utilized to handle the demand of the customers. It is also shown that

RI PT

the difference between the line-pack of the pipes for the stochastic and average scenario are notable. It is also clear from the results that the single end condition for the total line-pack of the system can restore the system to its initial condition. The

SC

trajectories are plotted for 30 hours to show the restoration of the system to its initial

is shown in Fig. 10 and Fig. 11.

M AN U

condition at the end of the gas day. The computed rotational speed for the compressors

It can be concluded that the use of appropriate set-points for each scenario will manipulate the line-pack of each pipe to meet all the loads imposed on the transmission line while all operational constraints are respected. The difference between

TE D

deterministic cases was not unexpected due their different onsets, loads and durations. However, the stochastic case is a representative of a range of expected possibilities which are unknown apriori and following the appropriate compressor set-points make

EP

the system ready to respond, whichever possibility occurs. The comparison of the

AC C

compressors operating conditions for the uncertain future scenario and the expected future condition clearly shows that in order to be able to cope with the customer load uncertainty, a different strategy other than a strategy on the basis of the expected condition, should be employed. The time integral of the power used in compression units over a 24-hour duration are compared in Table 3. The total daily energy used for the stochastic case is 9.58% 33

ACCEPTED MANUSCRIPT

higher compared to the realization of the “Min“, which is the highest price we pay for a robust preparation schedule. Preparation for uncertain future also needs 0.84% more energy compared to the expected condition. This extra energy is the price we pay to

RI PT

gain greater operational flexibility in handling customer load uncertainties.

The optimal value for the end point pressure as obtained using steady-state optimization is 700 psia that equals to its minimum contracted value. However, from the

SC

practical point of view, as the transmission line is under various disturbances, the

M AN U

minimum allowable pressure considered for the delivery point must be somewhat higher than the contracted threshold. This issue can be discussed in the deterministic cases as follows (Fig. 12). Initially, the end pressure of 700 psia raises to some values higher due to the requirement of packing in the pipe P6; however, it meets the threshold line in its transitions to prevent energy losses. In other words, the optimizer tries to keep the end

the system.

TE D

point pressure as close as possible to its minimum value, unless packing is required in

EP

As uncertainties might result in the violation of customer contracted pressure, a minimum pressure buffer as discussed in Section 2.3 should be considered in

AC C

optimization runs. The decision about the value of the minimum pressure buffer can be made on the basis of the experiences of the dispatchers. For example, a fixed value of 20 psi might be used. However, according to the temporal behavior of the pressure at the end delivery for the stochastic case, the pressure falls to the minimum value of the 715 psia at the time its deterministic counterparts reach 700 psia. This provides approximately 15 psi pressure buffer. Accordingly, stochastic programming provides an 34

ACCEPTED MANUSCRIPT

optimal way to determine the minimum pressure buffer and under the assumed conditions for our case study, the delivery pressure will remain above 700 psia with a

RI PT

pre-specified probability. We find from the temporal behavior of the supply mass flow rate for different cases as shown in Fig. 13 that during the preparation period (i.e., before the time when the load is imposed), the input of the pipeline is increased in order to prepare the system for the

SC

upcoming load. However, depending on the size and duration of the load, different

M AN U

packing rates are required. This packed gas makes the system flexible enough to respond to the imposed load by the power plant while all the constraints are satisfied. After that, the input gas is used to restore the system to its initial condition. Investigation of the supply mass flow rates for “Expected” and “St. 1“ as shown in Fig. 13, reveals that the difference between the temporal profiles of these two scenarios is noticeably

TE D

less than the difference between other operating parameters such as compressors rotational speed or line-pack. This issue might be interpreted as indicating that preparation for upcoming uncertain load conditions does not require extra natural gas

EP

supplies, and it is just the matter of the management of the available compression units

AC C

for appropriate line-pack positioning. 6.2.3 Effect of pj

The effect of different values of pj is also studied. The total energy for each stochastic scenario is shown in Table 4. The outlet pressures are also compared in Fig. 14. It can be seen that as the value of pj is increased, the pressure buffer, i.e., the distance 35

ACCEPTED MANUSCRIPT

between the outlet pressure profile and constant pressure line of 700 psia is increased in order to decrease the possibility of constraint violation. This results in a slight

RI PT

increase in the energy consumption as well. Also, for pj =0.999 and pj =.9999 the pressure buffers are 21 psia and 26 psia, respectively, corresponding to total energy requirements of 2160.1 and 2172.3 MWh. The small increase in the energy

SC

consumption originates from the fact that even for the case with pj = 0.99, the operating points of the compressors have enough distance from the boundaries of the wheel map

M AN U

that all the constraints imposed by the wheel map are inactive, indicating that the optimal operating points tend to be located at high efficiency regions rather than the wheel map boundaries which prevents the movements of the operating point with pj variations. The effects of the different values of the pj for the steady state counterpart

TE D

of the problem under study can be found in [33]. 7. Conclusions and future directions

The extremely large gas holdup in gas transmission networks makes taking preventive

EP

actions by gas network dispatchers essential in order to be able to overcome the

AC C

uncertain upcoming loads. In this work, a stochastic optimization framework for daily planning of natural gas transmission networks is proposed. It is found that the strategy provided by the solution of the stochastic optimization problem compared with its deterministic counterparts can lead to the operation of the system in a more efficient and reliable manner, regardless of the actual realized load.

36

ACCEPTED MANUSCRIPT

The requirement of the decision time as a drawback of the scenario-based approaches where different scenarios are assumed for the future loads is resolved. The historical data can be analyzed to quantify the uncertainty level and this information along with

RI PT

accurate methods for the estimation of the future loads can be incorporated into the optimization framework.

SC

For nθ uncertain parameters, 2 nθ + 1 and 3 nθ grid points are required for uncertainty propagation with the aid of the unscented transform and multidimensional integrals

M AN U

using Gaussian quadrature, respectively.

In a computational study, it is shown that the line-pack manipulation strategy plays the main role to handle uncertain future demands while a 0.84% increase in the total energy consumption is observed. This is a reasonable price to pay to provide our customers

TE D

with acceptable levels of service.

The solution of the stochastic optimization problem can also determine the optimal pressure buffer required in each delivery point instead of a pre-specified pressure buffer

EP

to assure the safe operation of the network. This makes stochastic programming a valuable tool for gas network (transmission and distribution) companies in planning for

AC C

future operations.

The methodology proposed in this work can be applied in larger networks as the one discussed in [12] to explore its benefits in risk management. Also, it can be used in combined gas and electricity network optimization for the optimal dispatching of power and gas, in order to minimize the integrated gas-electricity system operational risks [59]. 37

ACCEPTED MANUSCRIPT

References

AC C

EP

TE D

M AN U

SC

RI PT

[1] Vitullo SR. Disaggregating time series data for energy consumption by aggregate and individual customer. MARQUETTE UNIVERSITY; 2011. [2] Lundberg W, Westhoff M. Fundamentals Of Load Forecasting And Uncertainty WN, Pipeline Simulation Interest Group; 1997. [3] Kandepu R, Foss B, Imsland L. Applying the unscented Kalman filter for nonlinear state estimation. J Process Control 2008;18:753–768. [4] Behrooz HA, Boozarjomehry RB. Modeling and state estimation for gas transmission networks. J Nat Gas Sci Eng 2015;22:551–570. [5] Abbaspour M, Krishnaswami P, Chapman K. Transient optimization in natural gas compressor stations for linepack operation. J Energy Resour Technol 2007;129:314–324. [6] Gato LMC, Henriques JCC. Dynamic behaviour of high-pressure natural-gas flow in pipelines. Int J Heat Fluid Flow 2005;26:817–25. doi:10.1016/j.ijheatfluidflow.2005.03.011. [7] Abbaspour M, Krishnaswami P, Chapman K. Transient optimization in natural gas compressor stations for linepack operation. J Energy Resour Technol 2007;129:314–324. [8] Liu E, Li C, Yang Y. Optimal energy consumption analysis of natural gas pipeline. Sci World J 2014;2014. [9] Herrán-González A, De La Cruz JM, De Andrés-Toro B, Risco-Martín JL. Modeling and simulation of a gas distribution pipeline network. Appl Math Model 2009;33:1584–600. doi:10.1016/j.apm.2008.02.012. [10] Dorao C, Fernandino M. Simulation of transients in natural gas pipelines. J Nat Gas Sci Eng 2011;3:349–355. [11] Reddy HP, Narasimhan S, Bhallamudi SM. Simulation and state estimation of transient flow in gas pipeline networks using a transfer function model. Ind Eng Chem Res 2006;45:3853–3863. [12] Behrooz HA, Boozarjomehry RB. Distributed and decentralized state estimation in gas networks as distributed parameter systems. ISA Trans 2015;58:552–566. [13] Ríos-Mercado RZ, Kim S, Boyd EA. Efficient operation of natural gas transmission systems: A network-based heuristic for cyclic structures. Comput Oper Res 2006;33:2323–2351. [14] Mahlke D, Martin A, Moritz S. A simulated annealing algorithm for transient optimization in gas networks. Math Methods Oper Res 2007;66:99–115. [15] Ferruzzi G, Cervone G, Delle Monache L, Graditi G, Jacobone F. Optimal bidding in a Day-Ahead energy market for Micro Grid under uncertainty in renewable energy production. Energy 2016;106:194–202. [16] Marqués D, Morari M. On-line optimization of gas pipeline networks. Automatica 1988;24:455–69. doi:10.1016/0005-1098(88)90091-X. [17] Hans A. Optimal Control of Natural Gas Pipeline Networks: A Real-Time, ModelBased, Receding Horizon Optimisation Approach. VDM Publishing; 2008. 38

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

[18] Baumrucker B, Biegler LT. MPEC strategies for cost optimization of pipeline operations. Comput Chem Eng 2010;34:900–913. [19] Carter RG, Rachford Jr HH. Optimizing line-pack management to hedge against future load uncertainty, Pipeline Simulation Interest Group; 2003. [20] Zavala VM. Stochastic optimal control model for natural gas networks. Comput Chem Eng 2014;64:103–113. [21] Uraikul V, Chan CW, Tontiwachwuthikul P. Development of an expert system for optimizing natural gas pipeline operations. Expert Syst Appl 2000;18:271–82. doi:10.1016/S0957-4174(00)00009-9. [22] Nguyen HH, Chan CW. Applications of artificial intelligence for optimization of compressor scheduling. Eng Appl Artif Intell 2006;19:113–126. [23] Li P, Arellano-Garcia H, Wozny G. Chance constrained programming approach to process optimization under uncertainty. Comput Chem Eng 2008;32:25–45. [24] Arellano-Garcia H, Wozny G. Chance constrained optimization of process systems under uncertainty: I. Strict monotonicity. Comput Chem Eng 2009;33:1568– 1583. [25] Darlington J, Pantelides CC, Rustem B, Tanyi B. An algorithm for constrained nonlinear optimization under uncertainty. Automatica 1999;35:217–228. [26] Bahakim SS, Rasoulian S, Ricardez-Sandoval LA. Optimal design of large-scale chemical processes under uncertainty: A ranking-based approach. AIChE J 2014;60:3243–3257. [27] Crestaux T, Le Maıˆtre O, Martinez J-M. Polynomial chaos expansion for sensitivity analysis. Reliab Eng Syst Saf 2009;94:1161–1172. [28] Chaffart D, Rasoulian S, Ricardez-Sandoval LA. Distributional uncertainty analysis and robust optimization in spatially heterogeneous multiscale process systems. AIChE J 2016. [29] Pintarič ZN, Kasaš M, Kravanja Z. Sensitivity analyses for scenario reduction in flexible flow sheet design with a large number of uncertain parameters. AIChE J 2013;59:2862–2871. [30] Ostrovsky GM, Ziyatdinov NN, Lapteva TV. Optimization problem with normally distributed uncertain parameters. AIChE J 2013;59:2471–2484. [31] Rasoulian S, Ricardez-Sandoval LA. Robust multivariable estimation and control in an epitaxial thin film growth process under uncertainty. J Process Control 2015;34:70–81. [32] Nagy Z, Braatz R. Distributional uncertainty analysis using power series and polynomial chaos expansions. J Process Control 2007;17:229–240. [33] Behrooz HA. Managing demand uncertainty in natural gas transmission networks. J Nat Gas Sci Eng 2016;34:100–111. [34] Ahmadian Behrooz H. Robust design and control of extractive distillation processes under feed disturbances. Ind Eng Chem Res 2017;56:4446–4462. [35] Behrooz HA. Robust synthesis of the pressure-swing distillation process under azeotropic feed composition disturbance—Study of the Tetrahydrofuran/Methanol system. Comput Chem Eng 2017;104:211–230.

39

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

[36] Paramasivan G, Kienle A. Decentralized Control System Design under Uncertainty Using Mixed-Integer Optimization. Chem Eng Technol 2012;35:261– 271. [37] Julier SJ, Uhlmann JK. Unscented filtering and nonlinear estimation. Proc IEEE 2004;92:401–422. [38] Abbaspour M, Krishnaswami P, Chapman K. Transient optimization in natural gas compressor stations for linepack operation. J Energy Resour Technol 2007;129:314–324. [39] Solsvik J, Jakobsen HA. Effects of Jacobi polynomials on the numerical solution of the pellet equation using the orthogonal collocation, Galerkin, tau and least squares methods. Comput Chem Eng 2012;39:1–21. [40] Letniowski FW. Compressor station modeling in networks, Pipeline Simulation Interest Group; 1993. [41] Mohitpour M, Golshan H, Murray MA. Pipeline Design & Construction: A Practical Approach. ASME Press; 2003. [42] Nguyen HH, Uraikul V, Chan CW, Tontiwachwuthikul P. A comparison of automation techniques for optimization of compressor scheduling. Adv Eng Softw 2008;39:178–188. [43] Chen Q, Shi Y, Xu X. Combination Model for Short-Term Load Forecasting. Open Autom Control Syst J 2013;5:124–132. [44] Wallooppillai R, Laud M-L. Planning for power: Pipeline design to meet gas-fired power plant needs, Pipeline Simulation Interest Group; 2003. [45] American Gas Association 2010. http://www.aga.org. [46] Vitullo SR, Brown RH, Corliss GF, Marx BM. Mathematical models for natural gas forecasting. Can Appl Math Q 2009;17:807–827. [47] Sharma SK, Sharma V. ANN Approach for Daily Prediction of Gas Load Forecasting. 5th Natl. Conf Comput. Nation Dev. INDIACom, 2011. [48] Azadeh A, Asadzadeh S, Ghanbari A. An adaptive network-based fuzzy inference system for short-term natural gas demand estimation: uncertain and complex environments. Energy Policy 2010;38:1529–1536. [49] Sánchez-Úbeda EF, Berzosa A. Modeling and forecasting industrial end-use natural gas consumption. Energy Econ 2007;29:710–742. [50] De Gooijer JG, Hyndman RJ. 25 years of time series forecasting. Int J Forecast 2006;22:443–473. [51] McSharry PE, Bouwman S, Bloemhof G. Probabilistic forecasts of the magnitude and timing of peak electricity demand. Power Syst IEEE Trans On 2005;20:1166– 1172. [52] Li W, Hui C-W, Li P, Li A-X. Refinery planning under uncertainty. Ind Eng Chem Res 2004;43:6742–6755. [53] Feehery WF, Barton PI. Dynamic optimization with state variable path constraints. Comput Chem Eng 1998;22:1241–1256. [54] Julier SJ, Uhlmann JK. A general method for approximating nonlinear transformations of probability distributions. Technical report, Robotics Research Group, Department of Engineering Science, University of Oxford; 1996. 40

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

[55] Van Der Merwe R. Sigma-point Kalman filters for probabilistic inference in dynamic state-space models. Oregon Health & Science University, 2004. [56] Kandepu R, Imsland L, Foss BA. Constrained state estimation using the unscented Kalman filter. Control Autom. 2008 16th Mediterr. Conf. On, Citeseer; 2008, p. 1453–1458. [57] Audet C, Dennis Jr JE. Mesh adaptive direct search algorithms for constrained optimization. SIAM J Optim 2006;17:188–217. [58] Audet C, Dennis Jr JE. Analysis of generalized pattern searches. SIAM J Optim 2002;13:889–903. [59] Favuzza S, Graditi G, Sanseverino ER. Adaptive and dynamic ant colony search algorithm for optimal distribution systems reinforcement strategy. Appl Intell 2006;24:31–42. [60] Deming MM. A Theory of Weather-Sensitive Energy Demand, Pipeline Simulation Interest Group; 1993.

41

ACCEPTED MANUSCRIPT

Table 1 - Different type of natural gas customers and their characteristics Customer type characteristics

• Commercial

Industrial

• • • •

dependent on both weather and nonweather factors [60] decreasing consumption on weekends less temperature sensitive driven by economic factors [60] decreasing consumption on weekends -

Natural gas vehicles

-

AC C

EP

TE D

M AN U

Power generators

to heat homes, run appliances, and use water heaters;

RI PT

•

Strongly dependent on weather conditions Increasing consumption on weekends.

heating buildings

running boilers and as feedstock for industrial processes

SC

• Residential

Natural gas application

running turbines that drive an electric generator as a fuel

ACCEPTED MANUSCRIPT

Table 2 Power plant load characterizing parameters for different scenarios Required flow (Μ)

Onset time (ts)

Duration (td)

-

500.0 660.5 385.6

9.00 7.95 10.14

4.00 5.69 2.18

0.99 0.90 0.80 0.70

Ν (500.0, 502)

Ν (9.00, 0.3872)

Ν (4.00, 0.3872)

AC C

EP

TE D

M AN U

SC

“Expected” “Max” “Min” “St. 1” “St. 2” “St. 3” “St. 4”

Probability ( pj )

RI PT

Name

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

Table 3 Total daily energy requirement of different scenarios Scenario “Expected” “Max” “Min” “St. 1” Total energy requirement [MWh] 2137.4 2422.8 1967.0 2155.4

ACCEPTED MANUSCRIPT

Table 4 Total daily energy requirement of stochastic scenarios

pj

0.99

0.90

0.80

0.70

AC C

EP

TE D

M AN U

SC

RI PT

Scenario “St. 1” “St. 2” “St. 3” “St. 4” Total energy requirement [MWh-Hours] 2155.4 2155.0 2153.3 2151.7

ACCEPTED MANUSCRIPT Constant speed curves (rpm)

80.00

Constant efficiency curves 70.00

50.00

30.00 20.00 10.00 3.00

4.00

5.00

6.00

7.00

8.00

10.00

SC

Actual Flow [m3/s]

9.00

RI PT

40.00

EP

TE D

M AN U

Fig. 1. A typical centrifugal compressor wheel map

AC C

Head [kJ/kg]

60.00

11.00

12.00

ACCEPTED MANUSCRIPT 5900 5800 5700 5600

RI PT

5400

5100 5000 4900 0

20

40

60

Customer Delivery Pressure

SC

5200

Minimum Pressure Considered By Dispatcher Customer Minimum Contracted Pressure

M AN U

Pressure Buffer

Minimum Pressure Buffer

5300

80

100

120

140

160

Time [hr]

EP

TE D

Fig. 2. Transient behavior of pressure at a delivery point

AC C

Pressure [kPa]

5500

180

200

ACCEPTED MANUSCRIPT

Demand [MMSCFD]

6 5 4 3

1 0 0

2

4

ts

6

8

10

12

ts+t14 d

RI PT

Μ

2

16

18Time [hr] 20

AC C

EP

TE D

M AN U

SC

Fig. 3. A gas fired power plant load characterization

ACCEPTED MANUSCRIPT

7446.4 kPa

5913.1 kPa

5192.8 kPa

5818.8 kPa

6465.9 kPa

7885.5 kPa

7202.1 kPa

N2

P2

N3

C1

N4

P3

N5

N6

C2

AC C

EP

TE D

M AN U

SC

P4

N7

P5

RI PT

P1 N1

5596.7 kPa

D1

S1

N8

4826.3 kPa = 700 psia

D2

P6 N9

M AN U

SC

Frequency

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

Fig. 5. Distribution of the objective function at the first iteration

M AN U

SC

Frequency

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

Fig. 6. Distribution of the constraint corresponding to the end point delivery pressure at the first iteration

ACCEPTED MANUSCRIPT 0 0

5

10

Time [hr]

15

20

25

30

-0.5

%

-1 "Expected" -1.5

"Max"

P1 -2.5 0 0

5

10

Time [hr] 15

20

25

SC

-1

RI PT

"St. 1"

-2

-2 -3

M AN U

-4 -5 -6 -7 -8

"Expected"

P2

EP

TE D

Fig. 7. Pipes P1 and P2 line-pack relative to initial condition

AC C

%

"Min"

"Max" "Min" "St. 1"

30

ACCEPTED MANUSCRIPT 1

Time [hr] 0 -1

0

5

10

15

20

25

30

%

-2 -3

"Expected"

P3

-4

"Max" "Min"

RI PT

-5

-7

2

Time [hr] 0 -2

0

5

10

15

-4

M AN U

%

-6

SC

"St. 1"

-6

-8 -10 -12 -14 -16

20

25

"Expected"

P4

EP

TE D

Fig. 8. Pipes P3 and P4 line-pack relative to initial condition

AC C

30

"Max" "Min" "St. 1"

ACCEPTED MANUSCRIPT 9 8 7 6 "Expected"

%

5

"Max"

P5

4

"Min"

3

"St. 1"

2 0 -1 0

5

10

15

Time [hr] 10

6

25

2 0 0

5

10

-2

M AN U

4

15

20

"Max" "Min" "St. 1"

25

Time [hr]

TE D

Fig. 9. Pipes P5 and P6 line-pack relative to initial condition

EP

30

"Expected"

P6

AC C

%

20

SC

8

RI PT

1

30

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

790 780 "Expected"

770

RI PT

"Max"

760

"Min"

"St. 1"

740

SC

730 720

700 690 0

5

10

M AN U

710

15

20

Time [hr]

EP

TE D

Fig. 12. Pressure at the end delivery point (N9)

AC C

psia

750

25

30

ACCEPTED MANUSCRIPT

510 "Expected"

505

"Max" 500

"Min" "St. 1"

RI PT

490 485 480 475

465 0

5

10

15

Time [hr]

SC

470 20

EP

TE D

M AN U

Fig. 13. Supply mass flow rate (S1)

AC C

kg/s

495

25

30

ACCEPTED MANUSCRIPT

770 760 "St. 1" 750

"St. 2"

RI PT

"St. 3" "St. 4"

700 psia

730 720

SC

710

690 0

5

10

M AN U

700

15

20

Time [hr]

EP

TE D

Fig. 14. Pressure at the end delivery

AC C

psia

740

25

30

ACCEPTED MANUSCRIPT

Optimal operation and daily planning of natural gas networks are addressed

•

An optimization framework for handling customer demand uncertainties is proposed

•

Application of unscented transform in chance constrained optimization is shown

•

Uncertain load of gas-fired power plant is characterized

•

Role of line-pack manipulation strategy in uncertain load handling is illustrated

AC C

EP

TE D

M AN U

SC

RI PT

•