- Email: [email protected]
Minimizing fuel consumption of a gas pipeline in transient states by dynamic programming
Journal of Natural Gas Science and Engineering 28 (2016) 193e203
Contents lists available at ScienceDirect
Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse
Minimizing fuel consumption of a gas pipeline in transient states by dynamic programming Xiaorui Zhang, Changchun Wu*, Lili Zuo National Engineering Laboratory for Pipeline Safety/Beijing Key Laboratory of Urban Oil and Gas Distribution Technology, China University of Petroleum (Beijing), 18 Fuxue Road, Changping, Beijing, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 12 June 2015 Received in revised form 15 November 2015 Accepted 17 November 2015 Available online 22 November 2015
Gas pipelines are often operated in transient states because of time-varying consumer demands for natural gas. For a certain operating period of a gas pipeline, given the delivery rates at the pipeline's delivery points and the flow rates or pressures of its gas sources as functions of time, pipeline dispatchers face a transient optimization problem with the objective of minimizing pipeline fuel consumption for the period. To address this problem, a two-level model is constructed. The model includes a pipeline level that optimizes the discharge pressure of each compressor station and a station level that optimizes the operating scheme of each compressor station. A key assumption, that the discharge pressures of each compressor station on the pipeline are time-independent, is made to simplify the problem. Based on this assumption, a dynamic programming (DP)-based algorithm is proposed for the pipeline level in which the discharge pressures are chosen as the state variables. The algorithm is essentially a recursion from the pipeline destination to the origin. The i-th stage of the recursion covers the i-th and the (i þ 1)-th compressor stations together with the pipeline segment between them, and the state transition of the stage reflects the relationship between the discharge pressures of the two compressor stations. Unlike in steady-state operation, determining this relationship involves the partial differential equations governing transient gas flow in the pipeline segment between two adjacent stations. Because these equations must be solved by numerical simulation, it is much more difficult to determine the relationship for transient operation than for steady-state operation. To compress the gas going through a compressor station to a compression ratio assigned at the pipeline level, the operating scheme of the station is solved at the station level. A preprocessing method is devised to reduce the computational time of the DP approach. The method approximates the transient optimization problem as a series of steady-state optimization problems. The results provide information that can narrow the discharge pressure region of each compressor station. A simple example shows that the approximate minimum fuel consumption given by the proposed algorithm is 6.2% higher than the true optimum. Typical cases for different operating periods are tested on a gas pipeline. The proposed approach is also tested on a network studied in the literature. For most of the case studies, optimal solutions can be computed by the proposed approach within 15 min. In addition, the computational complexity grows modestly as the period becomes longer. © 2015 Elsevier B.V. All rights reserved.
Keywords: Gas pipeline Optimization Transient state Dynamic programming
1. Introduction The energy costs of gas pipelines usually constitute the majority of pipeline operating costs. For many gas pipelines, each compressor station is equipped with centrifugal compressors
* Corresponding author. E-mail addresses: [email protected] (X. Zhang), [email protected] (C. Wu), [email protected] (L. Zuo). http://dx.doi.org/10.1016/j.jngse.2015.11.035 1875-5100/© 2015 Elsevier B.V. All rights reserved.
driven by gas turbines, and the fuel consumed by the drivers is natural gas taken from the pipelines. Thus, the energy costs associated with operating gas pipelines depend mainly on the fuel consumption of the drivers. Many studies have been conducted to minimize pipeline energy consumption in recent decades (Wong and Larson, 1968; Carter, 1996; Wright et al., 1998; Ferber et al., 1999; Casoetto et al., 2011; Thanh et al., 2013; Fasihizadeh et al., 2014; Schmidt et al., 2014), most of which have been based on steady-state operation. The most
194
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
widely used method for the steady-state optimization of gas pipelines is dynamic programming (DP), which was first used for the optimization of gas pipeline operation in the late 1960s (Wong and Larson, 1968). The DP method offers advantages in searching for global optima and in robustness (Peretti and Toth, 1982; Anglard, 1988; Liu et al., 2014). In addition to single pipelines, the method has been applied to gas pipeline networks (Lall and Percell, 1990; Carter, 1998). Evolutionary algorithms, such as genetic algorithms (MohamadiBaghmolaei et al., 2014), ant colony optimization algorithms (Chebouba et al., 2009), and particle swarm optimization algorithms (Zheng and Wu, 2012; Wu et al., 2014), have also been used for steady-state optimization. However, gas pipelines are often operated in transient states because of time-varying consumer demand for natural gas. Hence, a transient optimization problem must be solved. This problem results in a nonlinear mixed-integer programming model with partial differential equations (PDEs) as constraints. In the model, the PDEs describe the transient dynamics of gas flow in a pipeline, the integer variables denote compressor states, and the control parameters are time-dependent. A limited number of studies addressing this problem have been reported in the literature (Dupont and Rachford, 1987; Krishnaswami et al., 2004; Rossi et al., 2011, 2012, 2014). The transient optimization problem is usually formulated as a nonlinear mixed-integer model by discretizing the PDE constraints. Mantri et al. (1985) reported a hybrid algorithm consisting of a general reduced gradient method and dynamic programming to solve the nonlinear model. A gradient-based approach was also proposed by Rachford and Carter (2000, 2004). Multiple numerical experiments have been carried out to demonstrate the feasibility of the method (Carter et al., 2001; Pietsch et al., 2001; Carter et al., 2010). Furthermore, this approach has enhanced the ability to cope with load forecast uncertainty (Carter and Rachford, 2003) and to optimize the operation of branch networks in transient states (Carter et al., 2006; Rachford et al., 2009). An industrial application of the improved approach was reported recently (Alfred et al., 2013). In other approaches, linear models are used to solve the transient optimization problem. For example, by formulating a linear model near the solution of each previous iteration, Sekirnjak (1996) and Kelling et al. (2000) solved this problem by sequential linear programming. Geissler et al. (2011) and Mahlke et al. (2007) linearized the model over its entire feasible region and devised a mixed-integer linear programming (MIP) algorithm to solve the problem. A simulated annealing method was used to provide a good initial solution for the MIP approach (Moritz, 2007). For these linear approaches, a compromise between algorithm speed and accuracy is inevitable. This paper addresses a specific type of fuel minimization problem in transient states. The consumer demands are given together with the source flow rates or pressures. At least one of these sources should offer gas at a prescribed pressure trend with respect to time in the case of a severe flow-rate imbalance and nonexistence of a feasible solution. The objective is to minimize the total fuel consumption of the gas pipeline during the operating period. The results include the optimal discharge pressures and the optimal compressor combinations for each station, in addition to the online compressor speeds. A DP-based approach is proposed in this paper to solve this type of problem. A two-level mathematical model is first formulated to describe the fuel minimization problem addressed in this paper. The model is then reformulated to fit the dynamic programming method, and the algorithm is described in detail. The effect of a key assumption made in developing the proposed algorithm is examined in Section 4. The proposed approach is then tested on an in-service gas
pipeline. Two case studies involving different operating periods, namely, 24 h and 168 h, are then described. The DP approach is also adopted to solve a transient optimization problem studied by Geißler et al. (2011) and Domschke et al. (2011), which involves a complex gas pipeline network. Finally, concluding remarks are presented. 2. Mathematical model In this section, a mathematical model describing the fuel minimization problem of a gas pipeline operated in transient states is formulated. The model features two levels, a pipeline level and a station level. At the pipeline level, the discharge pressures of each station at each time step are determined. Transient simulations are carried out to calculate the station suction pressures and flow rates. The station level verifies whether each compressor station can operate at the given flow rate, suction pressure, and discharge pressure imposed by the pipeline level. If a station can operate under these conditions, its minimum fuel consumption is computed and returned to the pipeline level. If it cannot, the pipeline level is informed of the failure, and the discharge pressure of the station is adjusted. Two simplifications are made in the model to make the problem easier: 1. It is assumed that gas flows isothermally in a pipeline. This assumption is made because taking the thermal calculations into account makes the fuel consumption minimization problem much more complex. Moreover, this assumption is widely used in addressing this problem (Rachford and Carter, 2000, 2004; Moritz, 2007; Geißler et al., 2011; Mahlke et al., 2007). 2. The fuel consumption of each compressor station is ignored in the flow-rate balance of the pipeline system. This simplification facilitates the construction of the DP approach. Furthermore, because the fuel consumption of a compressor station is very small compared with the total flow rate for a gas pipeline, this assumption creates only a small error in the isothermal simulation.
2.1. Objective function The objective of the model is to minimize the fuel consumption of a gas pipeline during an operating period Tend:
E¼
Tend XNCS Z i¼1
Fi ðQstd ðtÞ; Ps ðtÞ; Ts ðtÞ; Pd ðtÞÞdt;
(1)
t¼0
where Ncs is the number of compressor stations and Fi the fuel consumption of the i-th compressor station, which is related to its flow rate Qstd, suction pressure Ps, suction temperature Ts, and discharge pressure Pd. According to the trapezoidal rule of numerical integration, Eq. (1) can be approximated as follows:
E¼
XNcs XNtime i¼1
n¼1
0:5fFi ½ðn 1ÞDt þ Fi ðnDtÞgDt;
Fi ðnDtÞ ¼ Fi ½Qstd ðnDtÞ; Ps ðnDtÞ; Ts ðnDtÞ; Pd ðnDtÞ;
(2) (3)
where Dt is the time step and Ntime is the number of time steps. Hence, Ntime*Dt is equal to the operating period, Tend. As previously stated, the fuel consumption of a station F is computed at the station level. Although it is assumed that gas flows
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
isothermally in a pipeline, thermal calculations are carried out at the station level. Before gas leaves a station, it is cooled to a prescribed temperature. Moreover, only steady-state operation is considered at the station level because a compressor reacts much faster than a pipeline does if its operating condition is changed. Because identical compressors are arranged in parallel in each compressor station for the pipeline addressed in this paper, only the number of online compressors is adjusted to minimize the fuel consumption at the station. This approach is taken because in this type of station configuration, distributing the load among the online compressors equally is optimal once the number of compressors to run is determined (Wright, 1998; Carter, 1996; Mahmoudimehr and Sanaye, 2013; Zhang et al., 2014).
2.2. Constraints The model constraints arise mainly from simulating a compressor unit and from the transient dynamics of gas flow in a pipeline segment. The limits of the compressor speed, driver power, and operating pressure are also taken into consideration.
2.2.1. Simulating a compressor unit The operating envelope of a centrifugal compressor is defined by its maximum speed, minimum speed, surge line, and stone line, as illustrated in Fig. 1, Eq. (4), and Eq. (5). In Eq. (4), S is the compressor speed, and Smin and Smax are the minimum and maximum speeds, respectively. In Eq. (5), Qac represents the compressor flow rate under its suction pressure and suction temperature, and the coefficients a0, a1, a2, b0, b1, and b2 are compressor-specific. The coefficients are calculated by regressing the surge line and the stone line of the compressor.
Smin S Smax
(4)
a0 þ a1 S þ a2 S2 Qac b0 þ b1 S þ b2 S2
(5)
The performance map of the compressor is regressed using Eqs. (6) and (7) (Odom and Muster, 2009; Haijun et al., 2015). In Eq. (6), H is the polytropic head of the compressor, and in Eq. (7), hc denotes the compressor's efficiency. The coefficients c0,0, c0,1, …, c2,2, d0, d1, d2, and d3 are also compressor-specific and are computed by regressing the performance map of the compressor:
H S2 ¼ c0;0 þ c0;1 S þ c0;2 S2 þ c1;0 þ c1;1 S þ c1;2 S2 ðQac =SÞ þ c2;0 þ c2;1 S þ c2;2 S2 ðQac =SÞ2 hc ¼ d0 þ d1 ðQac =SÞ þ d2 ðQac =SÞ2 þ d3 ðQac =SÞ3
195
(6)
(7)
Compressing gas with a compressor is regarded as a polytropic process and expressed as shown in Eq. (8), where mV is the polytropic exponent. According to the compressor suction conditions and the compressor's discharge conditions, Eq. (9) calculates the polytropic head, where Zs is the gas compressibility factor under the compressor suction conditions, R is the universal gas constant, and MGAS is the molar mass of the gas. The polytropic efficiency of the compressor is computed using Eq. (10), where kave is the average specific capacity ratio, which is calculated using Eq. (11) according to the specific capacity ratio at the compressor suction, ks, and at the compressor discharge, kd:
mV 1 Td Zs Pd mV ¼ Ts Zd Ps H¼
mV 1 mV R Zs Ts ðPd =Ps Þ mV 1 mV 1 MGAS
hc ¼
kave 1 P T lg d lg d kave Ps Ts
kave ¼ ðks þ kd Þ=2
(8)
(9)
(10)
(11)
To simulate a compressor, the system of equations composed of Eqs. (6)e(11) and an equation of state are solved. The compressor input power is calculated using Eq. (12), where m denotes the compressor mass flow rate and hm the mechanical efficiency, which is assumed to be constant in most cases. The fuel consumption of the compressor unit is computed using Eq. (13), where hGT is the driver efficiency and LHV the lower heating value of the fuel. Note that only compressors driven by gas turbines are considered in this paper.
Pshaft ¼ ðHmÞ=ðhc hm Þ
(12)
. f ¼ Pshaft LHV=hGT
(13)
The efficiency performance map of a gas turbine is regressed using Eqs. (14)e(18). In these equations, Ta and Pa are the atmospheric temperature and pressure onsite, respectively, and Ta,0 and Pa,0 are the design atmospheric temperature and pressure, respectively. In addition, the maximum available power of the gas turbine can be computed using Eq. (19). This value constrains the input power of the compressor, as illustrated by Eq. (20). The coefficients e0, e1, … e5, f0, f1, … f5 in these equations are gas-turbinespecific and are computed by regressing the performance map of the turbine.
hGT ¼ e0 þ e1 nPT þ e2 NPT þ e3 ðnPT Þ2
2 þe4 NPT þ e5 nPT NPT nPT ¼ S Fig. 1. Operating envelope of a compressor.
.pffiffiffi q
(14)
(15)
196
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
q ¼ Ta Ta;0
(16)
first reformulated. The proposed approach is then described in detail.
. pffiffiffi N PT ¼ Pshaft d q
(17)
3.1. Model reformulation
d ¼ Pa Pa;0
(18)
Pmax ¼ f0 þ f1 S þ f2 Ta þ f3 STa þ f4 S2 þ f5 Ta2
(19)
Pshaft Pmax
(20)
2.2.2. Governing equations of transient gas flow in a pipe Eqs. (21) and (22) describe the transient dynamics of gas flow in a pipeline segment. In Eq. (21), r is the gas density and v the gas velocity. In Eq. (22), g is the acceleration due to gravity, P is the gas pressure, h is the pipe height, and d is the pipe's inner diameter. In addition, l is the friction coefficient, which is estimated using Eq. (23) according to the absolute roughness of the internal wall, k, and the Reynolds number, Re:
vr vðrvÞ þ ¼0 vt vx
(21)
vðrvÞ v rv2 dh vP l v2 þ r ¼ gr vt dx vx d 2 vx
(22)
1 k 2:51 pffiffiffi ¼ 2lg þ pffiffiffi 3:7d Re l l
(23)
2.2.3. Pressure limits The discharge pressure of a compressor station is constrained by the maximum allowable operating pressure (MAOP) of the pipeline:
Pd PMAOP ;
(24)
and the inlet pressure of a compressor station is constrained by its minimum allowable pressure:
Pi Pi;min ; i ¼ 1; 2; …; NCS
(25)
In addition, the pressure at each delivery point should not fall below the contract pressure:
P Pcontract
(26)
2.2.4. Flow-rate balance At each node, Kirchhoff's first law is applied:
X
in ¼ Qiso
X
out Qiso
(27)
As previously mentioned, the fuel consumption of each compressor station is ignored in the equation. 3. Dynamic programming approach A dynamic programming approach is described in this section to solve the model formulated in Section 2. Because dynamic programming models have a special structure, the previous model is
The model described in the previous section is reformulated to fit the dynamic programming approach. Elements of the reformulated model described in this section include the definition of a stage, the choice of state variables and decision variables, and the index function. 3.1.1. Stage The dynamic programming recursion process starts from the pipeline destination and progresses to its origin. In this procedure, each stage addresses the zone from the discharge of a downstream station to that of an adjacent upstream station. The gas pipeline terminal is regarded as a station without any compressor, and its origin is considered to be a station with only one feasible discharge pressure, i.e., the source pressure. Hence, there are a total of NCS þ 1 stages in the dynamic programming approach. 3.1.2. State variables In the DP model, the station discharge pressures are chosen as the state variables. These pressures should vary over time because transient states are considered in this formulation. However, the pressures are assumed to be constant because multiple simulations have revealed that a gas pipeline can operate with constant station discharge pressures if its consumer demands vary modestly over time. This assumption makes it possible to use the dynamic programming method, which is fast, easy, and robust, to minimize the fuel consumption of a gas pipeline in transient states. Indeed, the assumption may make the DP approach fail in severe transient situations. However, these situations are rare in the operation of trunk gas pipelines. Therefore, failure under these conditions should be considered a minor defect of the approach. The contract pressure at the pipeline terminal is regarded as the initial state. Note that the initial state can be arbitrary and that the choice does not affect the algorithm results. This result derives from the fact that the pipeline terminal is considered a station without any compressor, as previously mentioned. Consequently, only the suction pressure is of concern. The possible discharge pressure of the i-th station is [Piþ1,min, PMAOP], which is continuous. However, only finite states can be considered in a DP approach. Hence, the region is discretized using a fixed step size, 0.1 MPa. The resulting discretized set is further refined by the procedure described in the following subsection. Note that the lower bound for the last compressor station is the contract pressure of the terminal consumer, and only the source pressure is considered feasible for the pipeline origin, as previously described. 3.1.3. Decision variables and state transition For simplicity, in each stage, the upstream station is referred to as CSA and the adjacent downstream station is referred to as CSB. At the beginning of each stage, a set of possible CSB discharge pressures is determined. Corresponding to each discharge pressure, there is a flow-rate trend record of CSB. For the initial state, there is only one record, which is the terminal consumer demand. In each stage, the system state transitions from the CSA discharge pressure to that of CSB by selecting a discharge pressure value in the set determined at the beginning of the stage. Hence, the decision variables are also the station discharge pressures. 3.1.4. Stage effect Once the CSB discharge pressure has been determined, the
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
transient gas flow in the pipeline segment connecting CSA and CSB is simulated. In the simulation, the upstream boundary condition is the CSA discharge pressure and the downstream boundary condition is the flow-rate trend record corresponding to the CSB discharge pressure. In this paper, the transient simulation is carried out using Stoner Pipeline Simulator (SPS) 9.5 because of its widespread use in the gas pipeline industry, good agreement with field data, flexibility and robustness. SPS is also the only available simulator we could invoke from our optimization codes. The transient simulation yields the CSA flow-rate trend and the CSB suction pressure trend. The latter trend is examined to determine whether it falls below the minimum allowable CSB suction pressure at any time. If it does, another CSB discharge pressure is chosen, and the simulation process is repeated. If all possible CSB discharge pressures have been tested and the minimum allowable CSB suction pressure cannot be fulfilled in any of the cases, the CSA discharge pressure is considered infeasible, as is any other lower pressure value. Hence, the computational effort is reduced. If the CSB suction pressure trend stays above its minimum allowable suction pressure, the station level is invoked. This level examines whether CSB can operate at the given suction pressure, flow rate, and discharge pressure. If it can, its minimum fuel consumption during the operating period is computed and returned to the pipeline level. This fuel consumption is the stage effect. In addition, the CSA flow-rate trend is also recorded with its discharge pressure for future use. If CSB cannot operate under the given conditions, another CSB discharge pressure is chosen, and the process is repeated. Simulating transient gas-flow dynamics in a pipeline segment is computationally demanding, and numerous transient simulations are conducted in the DP approach. This process makes the computational time of this approach very long. However, these simulations are essential for obtaining a satisfactory and reliable solution. Coarse approximation may speed up the approach but at the risk of providing impractical results. Moreover, the DP approach is accelerated by detecting infeasible station discharge pressures, as previously described. The approach also benefits from selecting the best way to transition to a target state and considering previous decisions irrelevant once the transition is complete. 3.1.5. Index function The index function of the DP model is
F Pd;CSA ¼
min f Pd;CSA ; Pd;CSB þ F Pd;CSB ; Pd;CSB 2SðPd;CSB Þ
(28)
where S(Pd,CSB) is the set of possible CSB discharge pressures and f(Pd,CSA, Pd,CSB) is the CSB fuel consumption. If an infeasible situation occurs, f is set to a very large number. For the initial state, F(Pd,CSB) ¼ 0. 3.1.6. Markov property One important quality of a DP model is its Markov property. This property is related to the choice of state variables, which are the station discharge pressures in the current approach. In each stage, a transient simulation is conducted once the CSB discharge pressure is determined, and the simulation results are not related to previous decision procedures. Thus, the Markov property of the DP model is guaranteed. 3.2. Algorithm procedures Based on the reformulated model, the proposed algorithm is described in detail in the following steps:
197
0. Set k ¼ Ncs, Sstate ¼ {Pcontract}, Sobjfun ¼ {0}, and Scsdisflow ¼ {qenddelivery}. 1. Discretize the discharge pressure region [Pkþ1,min, PMAOP] with step size 0.1 MPa, and generate a set of discrete points Sdecision ¼ {P1, P2 … PNdecision}. 2. Sstate,pre ¼ Sstate, Sobjfun,pre ¼ Sobjfun, Scsdisflow,pre ¼ Scsdisflow, Sstate ¼ F, Sobjfun ¼ F, and Scsdisflow ¼ F. 3. Set i ¼ 1. 4. Set j ¼ 1. 5. Calculate the flow rate into the (k þ 1)-th station from the upstream pipeline by a flow-rate balance. Its discharge flow rate is Qj2Scsdisflow,pre. 6. Simulate the transient gas-flow dynamics in the pipeline segment that connects the k-th and the (k þ 1)-th stations. In the simulation, the upstream boundary condition is the constant suction pressure Pi2Sdecision, and the downstream boundary condition is the flow rate determined in Step 5. This transient simulation yields the discharge flow-rate trend of the k-th station and the suction pressure trend of the (k þ 1)-th station. 7. If the suction pressure of the (k þ 1)-th station stays above its minimum allowable suction pressure, go to 8. Otherwise, go to 11. 8. Examine whether the (k þ 1)-th station can operate at the given suction pressure, flow rate, and discharge pressure. If it can, calculate its fuel consumption, Fkþ1, and go to 9. If it cannot, go to 11. 9. If j ¼ 1, add Pi to Sstate, and record objfunj þ Fkþ1 as objfuni in Sobjfun, where objfunj is the fuel consumption of the previous stages corresponding to Qj and is stored in Sobjfun,pre. Record the discharge flow-rate trend of the k-th station in Scsdisflow. If j > 1, go to 10. 10. If objfunj þ Fkþ1 < objfuni, update objfuni and the corresponding discharge flow-rate trend stored in Scsdisflow. If objfunj þ Fkþ1 >¼ objfuni, go to 11. 11. j ¼ jþ1; if j > jSstate,prej, go to 12. Otherwise, go to 5. 12. If the suction pressure of the (k þ 1)-th station falls below its minimum allowable suction pressure for all of the flow-rate trends stored in Scsdislfow,pre, then Pi is considered infeasible and is eliminated from Sdecision together with all other lower pressure values. Otherwise, go to 13. 13. i ¼ iþ1; if i > jSdecisionj, go to 14. Otherwise, go to 4. 14. k ¼ k1; if k 0, go to 1. Otherwise, go to 15. 15. The minimum value in Sobjfun corresponds to the optimal solution.
3.3. Algorithm acceleration The forward method often requires more than 1 h to compute an optimal solution according to algorithm tests. We devise a preprocessing method to cut the computational time. By assuming that the gas pipeline is operated in a steady state at each time step, we approximate a transient optimization problem as a series of steadystate optimization problems. Then, each of them is solved. The minimum discharge pressure Pk,min,ss of each compressor station is searched among the solutions, as is the maximum pressure Pk,max,ss. Finally, the discharge pressure region [Pk,min,ss, Pk,max,ss] replaces [Pkþ1,min, PMAOP] in the DP algorithm. The former is typically much smaller than the latter. Thus, the algorithm can be sped up. To approximate a transient optimization problem as a series of steady-state optimization problems, the time step should be small. In the following tests, a time step of 1 h is used. In addition, for the last compressor station, Pmax,ss is set to min(Pmax,ss þ 500, PMAOP) in case the pressure at the pipeline end falls below the contract
198
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
pressure. 4. Algorithm tests The proposed DP approach is first tested on the Central AsiaChina gas pipeline, illustrated in Fig. 2. This pipeline extends across more than 1800 km and transports 30 billion cubic meters of natural gas per year from Turkmenistan and Kazakhstan to China. The pipeline is composed of two identical parallel pipelines, eight compressor stations, two sources, and four delivery points. Identical compressors driven by gas turbines are installed at each compressor station. Table 1 presents the number of compressors at each station and the speed limits of these compressors. Note that compressors in distinct stations may be different. The maximum discharge pressure of each station is 9810 kPa, and the minimum suction pressure is 6000 kPa. The contract pressure at C4 is 7500 kPa. To examine the effect of the key assumption made in Section 3, the proposed approach is used to minimize the fuel consumption of the last section of the gas pipeline. The result is compared with the true optimum to assess the effect of the assumption. Multiple tests are carried out on the pipeline for the original DP approach and the accelerated one. Two of these tests are described in detail below. The operating periods of these two tests are 24 h and 168 h. Finally, the accelerated DP approach is adopted to solve the problem studied by Geißler et al. (2011) and Domschke et al. (2011). This test involves a more complex network.
Table 1 Compressor station details. Station
Number of compressors
Minimum speed
Maximum speed
CS CS CS CS CS CS CS CS
3 3 3 3 3 4 4 4
3965 3965 3965 3965 3120 3965 3120 3120
6405 6405 6405 6405 5040 6405 5040 5040
1 2 3 4 5 6 7 8
4.1. Influence of the time-independent assumption To examine the effect of the assumption that the discharge pressure of each station is time-independent, the proposed approach is used to optimize the operation of the last section of the gas pipeline in transient states. In the problem, the suction pressure of the last compressor station and the terminal consumer demand are given. The fuel consumption of the last station is minimized over a certain period. This problem is chosen because its true optimum is available (Wong and Larson, 1968; Goldberg, 1983), and the method proposed by Wong and Larson (1968) is used to compute this true optimum. However, to make the comparison fair, the compressor station discharge pressure is chosen as the decision variable instead of its flow rate, as in the original paper, and SPS is used to perform transient simulations. A typical demand at the pipeline destination is plotted in Fig. 3. For the test problem, the optimal fuel consumption is 23.89 104 Sm3/d, whereas that calculated by the proposed method is 25.37 104 Sm3/d, or 6.2% higher. The optimum solution is depicted in Fig. 4 together with the solution provided by the proposed algorithm, which is close to the average of the optimal solution. The pressure at C4 is plotted in Fig. 5. For the true optimum, the pressure at C4 remains slightly above its lower bound except during the first few time steps. This result shows that the pressure energy is fully used to transport the gas. In contrast, for the proposed algorithm, the pressure at C4 is somewhat higher than
Fig. 3. Demand at C4.
Fig. 4. Optimal solution.
the contract pressure, which means that part of the pressure energy is wasted. Although a definitive conclusion cannot be drawn from a single case study, the study can at least shed some light on this issue. Because the state of a gas pipeline rarely fluctuates severely, it is believed that the solution computed by the proposed algorithm is not too far from the true optimum that in most cases arises during daily operations. Fig. 2. Topological structure of the Central Asia-China gas pipeline.
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
199
Fig. 7. Demand at C4. Fig. 5. Pressure at C4.
4.2. 24-hour case In the current case, S1 supplies gas at a constant pressure (7000 kPa) and S2 supplies gas at a constant flow rate (41.667 104 Sm3/h). Fig. 6 shows the demands of intermediate consumers, and the demand at C4 is plotted in Fig. 7. Most of the gas is delivered to the pipeline terminal, and the demand at C2 fluctuates severely within the operating period. Note that the flow rates considered in this paper are flow-rate values under standard conditions of 101.325 kPa and 20 C, as used in the Chinese petroleum industry. Table 2 shows the optimal solution computed by the original DP method. The optimal discharge pressures of most stations are equal to the maximum allowable operating pressure. Note that the design flow rate of the pipeline is approximately 357 104 Sm3/h, whereas the total flow rate in this case study is 378.04 104 Sm3/h. The pipeline system operates under a heavy load. Hence, the discharge pressure of most compressor stations in the solution is equal to the maximum allowable pressure. Fig. 8 shows the discharge flow rate of each compressor station, and their suction pressures are plotted in Fig. 9. The figures show that CS 5 and CS 6, which are located adjacent to C2, are most
Table 2 Optimum discharge pressure. Station
Discharge pressure
CS CS CS CS CS CS CS CS
9810 9810 9810 9810 9810 9810 9310 9110
1 2 3 4 5 6 7 8
Fig. 8. Discharge flow rate of each compressor station.
Fig. 6. Demands of intermediate consumers.
strongly affected by their demand fluctuation. With increasing distance, the effect diminishes gradually. In addition, the suction pressure of each compressor station remains far above the minimum allowable suction pressure (6000 kPa), as illustrated in Fig. 9. The pressure at C4 also remains above the contract pressure, 7500 kPa, during the operating period, as shown in Fig. 10. However, the pressure is not far from the lower bound. This situation is the result of efforts to minimize fuel consumption because higher pressure indicates extra compression power. The station level is also examined. CS 5 and CS 6 are investigated carefully because their operating conditions fluctuate severely. Two compressors are online in CS 5 during the operating period. Because identical compressors are installed in CS 5, the speed and input power of only one compressor are plotted in Fig. 11. The figure
200
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
Fig. 9. Suction pressure of each compressor station.
Fig. 12. Speed and fuel consumption of a compressor in CS 6.
plotted in Fig. 12. The figure shows that the compressor speed is flexible and that the input power is no greater than 16,600 kW, whereas the maximum available power of its driver is never less than 27,050 kW. Table 3 shows the solution computed by the accelerated algorithm. The solution is nearly the same as that computed by the original approach. After acceleration, the objective function value increases only by 0.72%, and the computational time is reduced by 92.29% e464 s. Twenty-four more tests are carried out. For twenty tests, the accelerated algorithm can compute an optimal solution within 15 min. The average computational time is 6.6 min. For the four remaining tests, the computational time is set to 1 h. The average computational time for each test is 15.2 min. In addition, the objective function value of the optimal solutions increases by 1.22% on average after acceleration. Fig. 10. Pressure at C4.
4.3. 168-hour case In the current case, all of the gas is delivered to C4, and Fig. 13 shows the corresponding demand, which fluctuates dramatically during the operating period. The sources supply gas under the same conditions as those in the 24-h case. The optimal solution yielded by the original DP approach is presented in Table 4. The discharge pressure of each station is equal to the maximum allowable operating pressure, except for CS 8. This result indicates that the pipeline capacity is nearly reached. Fig. 14 displays the discharge flow rate of each station, and the suction pressure trends are plotted in Fig. 15. CS 8 and CS 7 are markedly affected by the demand fluctuation at C4. With increasing distance, the influence diminishes gradually. Moreover, the suction pressure of each station remains well above the minimum allowable suction pressure, 6000 kPa, as shown in Fig. 15. The pressure at C4 is illustrated in Fig. 16. Although the pressure
Fig. 11. Speed and fuel consumption of a compressor in CS 5.
shows that the compressor speed varies within the speed limits and that the input power remains below 14,500 kW. Note that the maximum available power of a gas turbine varies as its speed changes. However, this available power for the current case is never less than 24,303 kW, which is far more power than required. Two compressors are also online in CS 6 during the operating period. The compressor speed and input power of one compressor are
Table 3 Optimum discharge pressure. Station
Discharge pressure
CS CS CS CS CS CS CS CS
9810 9810 9810 9810 9810 9718 9526 9334
1 2 3 4 5 6 7 8
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
201
Fig. 15. Suction pressure of the compressor stations.
Fig. 13. Demand at C4.
Table 4 Optimum discharge pressure. Station
Discharge pressure
CS CS CS CS CS CS CS CS
9810 9810 9810 9810 9810 9810 9810 9710
1 2 3 4 5 6 7 8
Fig. 16. Pressure at C4.
Fig. 14. Discharge flow rates of the compressor stations.
fluctuates significantly, it still remains above the contract pressure, 7500 kPa. The station level is also examined carefully, and no constraint is found to be violated. However, because the minimum up time and the minimum down time of a compressor unit are not taken into consideration in the current model, some compressors may run for a short period and then be turned off. This problem is found in CS 8, as shown in Fig. 17, and should be addressed in future research. Table 5 presents the optimal solution computed by the accelerated algorithm. The solution is nearly the same as that offered by the original approach. The accelerated algorithm computes this solution in approximately 32 min e 10.2% of the computational time of the original algorithm. However, the objective function
Fig. 17. Number of operating compressors in CS 8.
value increases by only 0.42%. Although the computational time is relatively long, it is considered acceptable considering the long operating period of the problem and the details of the optimization model. In addition, the computational time increases modestly as the operating period
202
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203 Table 5 Optimum discharge pressure.
Table 6 Pipe details.
Station
Discharge pressure
Pipe
Length (km)
Outer diameter (mm)
Thickness (mm)
CS CS CS CS CS CS CS CS
9810 9810 9810 9810 9810 9810 9716 9616
PL1 PL2 PL3 PL4 PL5 PL6 PL7 PL8 PL9 PL10 PL11 PL12
100 13 182 159 70 136 110 165 101 45 200 220
1016 1016 1016 1016 457 762 660 762 660 660 660 660
21 21 21 21 17.5 19 19 19 19 19 19 19
1 2 3 4 5 6 7 8
becomes longer. Two factors contribute to this effect. First, the assumption that the discharge pressure of each compressor station is time-independent makes problems with a longer operating period more constrained, leading to fewer feasible solutions. Second, the computational time for a transient simulation does not grow linearly as the operating period becomes longer.
4.4. Test on a more complex network In this section, we use the accelerated algorithm to solve the problem studied by Geißler et al. (2011) and Domschke et al. (2011). Fig. 18 shows the pipeline network. The transient process lasts for 4 h. Source1 supplies gas at a constant pressure of 7000 kPa, and source2 supplies gas at a constant flow rate of 36 104 Sm3/h. The flow rate of all of the deliveries except for sink A04 increases linearly from 36 104 Sm3/h to 42.235 104 Sm3/h. Sink A04 delivers gas at a constant flow rate of 18.5 104 Sm3/h. In their respective studies, Geißler et al. (2011) and Domschke et al. (2011) offered only part of the network details. To conduct the test, assumptions are made about the missing information. Network elevation change is ignored. Table 6 shows the pipe details. The roughness of all of the pipes is 0.01 mm. The maximum discharge pressure of each compressor station is 10,000 kPa, the minimum suction pressure of each station is 6000 kPa, and the contract pressure of each delivery is 4000 kPa. With the preceding assumptions, the accelerated algorithm solves this problem within 10 min. The optimum discharge pressures of CS 1 and CS 2 are 10,000 kPa and 9500 kPa, respectively, and CS 3 is bypassed.
5. Conclusions Based on the preceding analyses and tests, we draw the following conclusions: 1. The proposed approach was used to minimize the fuel consumption of the last section of a gas pipeline. Its result was only 6.2% higher than the true optimal. 2. The proposed approach was also tested on various problems arising in daily gas pipeline operations. For most case studies, the approach can yield a solution within 15 min.
Fig. 18. Network topology.
3. The computational effort of the proposed algorithm grows modestly as the operating period becomes longer. To increase the applicability of the proposed algorithm, we suggest carrying out the following work in the future: 1. Consider the minimum up time and the minimum down time of a compressor unit at the station level. 2. The SPS simulator we implemented can only use a portion of the power of multi-core CPUs. These CPUs are common today, and full use of all of the processors can greatly reduce the computational time. Therefore, we suggest replacing the simulator with one that is capable of conducting parallel computations. Acknowledgments This research was supported by the Australian and Western Australian Governments and the North West Shelf Joint Venture Partners, as well as the Western Australian Energy Research Alliance through the Australia-China Natural Gas Technology Partnership Fund Top Up Scholarship. This study was also supported by the Beijing Higher Education Young Elite Teacher Project (No. YETP0677). We also thank the anonymous reviewer for his/her insightful comments, which helped us improve our manuscript. Nomenclature d Fe fe ge He he ke ks, kd e
inner diameter of a pipe, m fuel consumption of a compressor station, Nm3/h fuel consumption of a compressor, Nm3/h gravity, 9.81 m/s2 polytropic head of a compressor, kJ/kg height, m absolute roughness of an internal pipe wall, m heat capacity ratio of natural gas at the suction and discharge of a compressor LHV e lower heating value of the fuel, kJ/Nm3 mV e polytropic exponent MGAS e molecular mass, kg/kmol Ncs e number of compressor stations Ntime e number of time steps P0 e ISO pressure, 101.325 kPa Pa e ambient pressure at a site, kPa Pa,0 e ambient pressure of a design, 101.325 kPa Pd e discharge pressure of a compressor or a station, kPa Pmax e maximum output power of a gas turbine, kW Ps e suction pressure of a compressor or a station, kPa Pshaft e input power of a compressor, kW Qac e volumetric flow rate under specific conditions, m3/h
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
Qiso e volumetric flow rate under standard conditions, Nm3/h Re e reynolds number Re gas constant, 8.314 kJ/(kg*K) Se compressor speed, RPM Scsdisflow set of possible discharge flow-rate patterns of a station Sdecision e set of possible decisions Sendflow e flow rate at the end of a pipeline Sobjfun e set of fuel consumption records of previous stages te time, s Te optimization horizon, min T0 e ISO temperature, 20 C Ta e ambient temperature at a site, K Ta,0 e ambient design temperature, 288.15 K Td e discharge temperature of a compressor or a station, C Ts e suction temperature of a compressor or a station, C we gas velocity, m/s xe length, m Zs e compressibility factor of natural gas at the station suction Dt e time step, min hC e compressor efficiency, [0,1] hM e mechanical efficiency, 0.98 hGT e gas turbine efficiency, [0,1] re gas density, kg/m3 le friction coefficient of a pipe References Alfred, S., Fasullo, J., Pfister, J., Daniels, A., 2013, April. Capacity determination using state finding and gas transient optimization. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Anglard, P., 1988, January. Hierarchical steady state optimization of very large gas pipelines. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Carter, R.G., 1996, January. Compressor station optimization: computational accuracy and speed. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Carter, R.G., 1998, October. Pipeline optimization: dynamic programming after 30 years. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Carter, R., Goodreau, M., Rachford, H., 2001. Optimizing pipeline operations through mathematical advances. Pipeline Gas J. 228 (10), 51e53. Carter, R.G., Rachford Jr., H.H., 2003, January. Optimizing line-pack management to hedge against future load uncertainty. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Carter, R.G., Rachford, H.H., Dupont, T.F., 2006, January. A new approach to gas control simulation. In: 2006 International Pipeline Conference. American Society of Mechanical Engineers, pp. 865e875. Carter, R., Reisner, M., Sekirnjak, E., 2010, January. Transient optimization-examples and directions. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Casoetto, B., Flottes, E., Ardeois, J., Thanh, N., Joliot, J.B., Nait-Abdallah, R., 2011. How to commercialize reliable capacities on a complex transmission network? J. Nat. Gas Sci. Eng. 3 (5), 657e663. Chebouba, A., Yalaoui, F., Smati, A., Amodeo, L., Younsi, K., Tairi, A., 2009. Optimization of natural gas pipeline transportation using ant colony optimization. Comput. Oper. Res. 36 (6), 1916e1923. Domschke, P., Geißler, B., Kolb, O., Lang, J., Martin, A., Morsi, A., 2011. Combination of nonlinear and linear optimization of transient gas networks. Inf. J. Comput. 23 (4), 605e617. Dupont, T.F., Rachford Jr., H.H., 1987. Optimization of power usage in transient gas transmission lines. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Fasihizadeh, M., Sefti, M.V., Torbati, H.M., 2014. Improving gas transmission networks operation using simulation algorithms: case study of the national Iranian gas network. J. Nat. Gas Sci. Eng. 20, 319e327. Ferber, P., Basu, U., Venkataramanan, G., Goodreau, M., Linden, P., 1999, January. Gas pipeline optimization. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Geißler, B., Kolb, O., Lang, J., Leugering, G., Martin, A., Morsi, A., 2011. Mixed integer linear models for the optimization of dynamical transport networks. Math. Methods Oper. Res. 73 (3), 339e362.
203
Goldberg, D.E., 1983. Computer-Aided Gas Pipeline Operation Using Genetic Algorithms and Rule Learning. Doctoral Dissertation. University of Michigan. Haijun, C., Changchun, W., Xiaorui, Z., 2015, January. Corrected maps narrow compressor performance prediction range. Oil Gas J. 113 (1), 72e79. Kelling, C., Reith, K., Sekirnjak, E., 2000, October. A practical approach to transient optimization for gas networks. In: PSIG Meeting, vol. 5. Krishnaswami, P., Chapman, K.S., Abbaspour, M., 2004, January. Compressor station optimization for linepack maintenance. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Liu, E., Li, C., Yang, Y., 2014. Optimal energy consumption analysis of natural gas pipeline. Sci. World J. 2014. Lall, H.S., Percell, P.B., 1990. A dynamic programming based gas pipeline optimizer. In: Analysis and Optimization of Systes. Springer Berlin Heidelberg, pp. 123e132. Mahlke, D., Martin, A., Moritz, S., 2007. A simulated annealing algorithm for transient optimization in gas networks. Math. Methods Oper. Res. 66 (1), 99e115. Mahmoudimehr, J., Sanaye, S., 2013. Minimization of fuel consumption of natural gas compressor stations with similar and dissimilar turbo-compressor units. J. Energy Eng. 140 (1). Mantri, V.B., Prestion, L.B., Pringle, C.S., 1985. Transient optimization of a natural gas pipeline system. In: PSIG Annual Meeting. MohamadiBaghmolaei, M., Mahmoudy, M., Jafari, D., MohamadiBaghmolaei, R., Tabkhi, F., 2014. Assessing and optimization of pipeline system performance using intelligent systems. J. Nat. Gas Sci. Eng. 18, 64e76. Moritz, S., 2007. A Mixed Integer Approach for the Transient Case of Gas Network Optimization. Doctoral Dissertation. TU Darmstadt. Odom, F.M., Muster, G.L., 2009, January. Tutorial on modeling of gas turbine driven centrifugal compressors. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Peretti, A., Toth, P., 1982. Optimization of a pipe-line for the natural gas transportation. Eur. J. Oper. Res. 11 (3), 247e254. Pietsch, U., Rachford Jr., H.H., Carter, R.G., 2001, January. Investigating real-world applications of transient optimization. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Rachford Jr., H.H., Carter, R.G., 2000, January. Optimizing pipeline control in transient gas flow. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Rachford Jr., H.R., Carter, R.G., 2004. U.S. Patent No. 6,701,223. U.S. Patent and Trademark Office, Washington, DC. Rachford Jr., H.H., Carter, R.G., Dupont, T.F., 2009, January. Using optimization in transient gas transmission. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Rossi, B., Casoetto, B., Le Maitre, A., Jacquiau, A., Moumini, N., 2012, January. Computing and optimizing the available linepack in a gas network to address the issue of flexibility. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Rossi, B., Dalphin, J., Sinegre, L., Jacquiau, A., Renaudie, T., 2011, January. Dealing with CCGT: an explicit dynamic model. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Rossi, B., Djerourou, F., Suez-Crigen, G.D.F., 2014, May. Transient optimization in gas transmission networks, a new approach on GRTgaz network. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Sekirnjak, E., 1996, October. Practical experiences with various optimization techniques for gas transmission and distribution systems. In: Proceedings of the 28th Annual Meeting of the Pipeline Simulation Interest Group, pp. 1e16 (PSIG'96). Schmidt, M., Steinbach, M.C., Willert, B.M., 2014. High Detail Stationary Optimization Models for Gas Networks: Validation and Results. Friedrich-AlexanderUniversit€ at Erlangen-Nürnberg, Department Mathematik. Tech. rep. Submitted. Thanh, N., le Maitre, A., Ardeois, J., Joliot, J.B., 2013, April. An optimization tool for network capacities. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Wong, P., Larson, R., 1968. Optimization of natural-gas pipeline systems via dynamic programming. Automatic control. IEEE Trans. 13 (5), 475e481. Wright, S., Ditzel, C., Somani, M., 1998, January. Compressor station optimization. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Wu, X., Li, C., Jia, W., He, Y., 2014. Optimal operation of trunk natural gas pipelines via an inertia-adaptive particle swarm optimization algorithm. J. Nat. Gas Sci. Eng. 21, 10e18. Zhang, X., Wu, C., Zuo, L., Meng, X., 2014, November. Dynamic programming based algorithm for compressor station optimization. In: ASME 2014 International Mechanical Engineering Congress and Exposition. American Society of Mechanical Engineers (Pp. V011T14A028-V011T14A028). Zheng, Z., Wu, C., 2012. Power optimization of gas pipelines via an improved particle swarm optimization algorithm. Petroleum Sci. 9 (1), 89e92.
Contents lists available at ScienceDirect
Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse
Minimizing fuel consumption of a gas pipeline in transient states by dynamic programming Xiaorui Zhang, Changchun Wu*, Lili Zuo National Engineering Laboratory for Pipeline Safety/Beijing Key Laboratory of Urban Oil and Gas Distribution Technology, China University of Petroleum (Beijing), 18 Fuxue Road, Changping, Beijing, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 12 June 2015 Received in revised form 15 November 2015 Accepted 17 November 2015 Available online 22 November 2015
Gas pipelines are often operated in transient states because of time-varying consumer demands for natural gas. For a certain operating period of a gas pipeline, given the delivery rates at the pipeline's delivery points and the flow rates or pressures of its gas sources as functions of time, pipeline dispatchers face a transient optimization problem with the objective of minimizing pipeline fuel consumption for the period. To address this problem, a two-level model is constructed. The model includes a pipeline level that optimizes the discharge pressure of each compressor station and a station level that optimizes the operating scheme of each compressor station. A key assumption, that the discharge pressures of each compressor station on the pipeline are time-independent, is made to simplify the problem. Based on this assumption, a dynamic programming (DP)-based algorithm is proposed for the pipeline level in which the discharge pressures are chosen as the state variables. The algorithm is essentially a recursion from the pipeline destination to the origin. The i-th stage of the recursion covers the i-th and the (i þ 1)-th compressor stations together with the pipeline segment between them, and the state transition of the stage reflects the relationship between the discharge pressures of the two compressor stations. Unlike in steady-state operation, determining this relationship involves the partial differential equations governing transient gas flow in the pipeline segment between two adjacent stations. Because these equations must be solved by numerical simulation, it is much more difficult to determine the relationship for transient operation than for steady-state operation. To compress the gas going through a compressor station to a compression ratio assigned at the pipeline level, the operating scheme of the station is solved at the station level. A preprocessing method is devised to reduce the computational time of the DP approach. The method approximates the transient optimization problem as a series of steady-state optimization problems. The results provide information that can narrow the discharge pressure region of each compressor station. A simple example shows that the approximate minimum fuel consumption given by the proposed algorithm is 6.2% higher than the true optimum. Typical cases for different operating periods are tested on a gas pipeline. The proposed approach is also tested on a network studied in the literature. For most of the case studies, optimal solutions can be computed by the proposed approach within 15 min. In addition, the computational complexity grows modestly as the period becomes longer. © 2015 Elsevier B.V. All rights reserved.
Keywords: Gas pipeline Optimization Transient state Dynamic programming
1. Introduction The energy costs of gas pipelines usually constitute the majority of pipeline operating costs. For many gas pipelines, each compressor station is equipped with centrifugal compressors
* Corresponding author. E-mail addresses: [email protected] (X. Zhang), [email protected] (C. Wu), [email protected] (L. Zuo). http://dx.doi.org/10.1016/j.jngse.2015.11.035 1875-5100/© 2015 Elsevier B.V. All rights reserved.
driven by gas turbines, and the fuel consumed by the drivers is natural gas taken from the pipelines. Thus, the energy costs associated with operating gas pipelines depend mainly on the fuel consumption of the drivers. Many studies have been conducted to minimize pipeline energy consumption in recent decades (Wong and Larson, 1968; Carter, 1996; Wright et al., 1998; Ferber et al., 1999; Casoetto et al., 2011; Thanh et al., 2013; Fasihizadeh et al., 2014; Schmidt et al., 2014), most of which have been based on steady-state operation. The most
194
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
widely used method for the steady-state optimization of gas pipelines is dynamic programming (DP), which was first used for the optimization of gas pipeline operation in the late 1960s (Wong and Larson, 1968). The DP method offers advantages in searching for global optima and in robustness (Peretti and Toth, 1982; Anglard, 1988; Liu et al., 2014). In addition to single pipelines, the method has been applied to gas pipeline networks (Lall and Percell, 1990; Carter, 1998). Evolutionary algorithms, such as genetic algorithms (MohamadiBaghmolaei et al., 2014), ant colony optimization algorithms (Chebouba et al., 2009), and particle swarm optimization algorithms (Zheng and Wu, 2012; Wu et al., 2014), have also been used for steady-state optimization. However, gas pipelines are often operated in transient states because of time-varying consumer demand for natural gas. Hence, a transient optimization problem must be solved. This problem results in a nonlinear mixed-integer programming model with partial differential equations (PDEs) as constraints. In the model, the PDEs describe the transient dynamics of gas flow in a pipeline, the integer variables denote compressor states, and the control parameters are time-dependent. A limited number of studies addressing this problem have been reported in the literature (Dupont and Rachford, 1987; Krishnaswami et al., 2004; Rossi et al., 2011, 2012, 2014). The transient optimization problem is usually formulated as a nonlinear mixed-integer model by discretizing the PDE constraints. Mantri et al. (1985) reported a hybrid algorithm consisting of a general reduced gradient method and dynamic programming to solve the nonlinear model. A gradient-based approach was also proposed by Rachford and Carter (2000, 2004). Multiple numerical experiments have been carried out to demonstrate the feasibility of the method (Carter et al., 2001; Pietsch et al., 2001; Carter et al., 2010). Furthermore, this approach has enhanced the ability to cope with load forecast uncertainty (Carter and Rachford, 2003) and to optimize the operation of branch networks in transient states (Carter et al., 2006; Rachford et al., 2009). An industrial application of the improved approach was reported recently (Alfred et al., 2013). In other approaches, linear models are used to solve the transient optimization problem. For example, by formulating a linear model near the solution of each previous iteration, Sekirnjak (1996) and Kelling et al. (2000) solved this problem by sequential linear programming. Geissler et al. (2011) and Mahlke et al. (2007) linearized the model over its entire feasible region and devised a mixed-integer linear programming (MIP) algorithm to solve the problem. A simulated annealing method was used to provide a good initial solution for the MIP approach (Moritz, 2007). For these linear approaches, a compromise between algorithm speed and accuracy is inevitable. This paper addresses a specific type of fuel minimization problem in transient states. The consumer demands are given together with the source flow rates or pressures. At least one of these sources should offer gas at a prescribed pressure trend with respect to time in the case of a severe flow-rate imbalance and nonexistence of a feasible solution. The objective is to minimize the total fuel consumption of the gas pipeline during the operating period. The results include the optimal discharge pressures and the optimal compressor combinations for each station, in addition to the online compressor speeds. A DP-based approach is proposed in this paper to solve this type of problem. A two-level mathematical model is first formulated to describe the fuel minimization problem addressed in this paper. The model is then reformulated to fit the dynamic programming method, and the algorithm is described in detail. The effect of a key assumption made in developing the proposed algorithm is examined in Section 4. The proposed approach is then tested on an in-service gas
pipeline. Two case studies involving different operating periods, namely, 24 h and 168 h, are then described. The DP approach is also adopted to solve a transient optimization problem studied by Geißler et al. (2011) and Domschke et al. (2011), which involves a complex gas pipeline network. Finally, concluding remarks are presented. 2. Mathematical model In this section, a mathematical model describing the fuel minimization problem of a gas pipeline operated in transient states is formulated. The model features two levels, a pipeline level and a station level. At the pipeline level, the discharge pressures of each station at each time step are determined. Transient simulations are carried out to calculate the station suction pressures and flow rates. The station level verifies whether each compressor station can operate at the given flow rate, suction pressure, and discharge pressure imposed by the pipeline level. If a station can operate under these conditions, its minimum fuel consumption is computed and returned to the pipeline level. If it cannot, the pipeline level is informed of the failure, and the discharge pressure of the station is adjusted. Two simplifications are made in the model to make the problem easier: 1. It is assumed that gas flows isothermally in a pipeline. This assumption is made because taking the thermal calculations into account makes the fuel consumption minimization problem much more complex. Moreover, this assumption is widely used in addressing this problem (Rachford and Carter, 2000, 2004; Moritz, 2007; Geißler et al., 2011; Mahlke et al., 2007). 2. The fuel consumption of each compressor station is ignored in the flow-rate balance of the pipeline system. This simplification facilitates the construction of the DP approach. Furthermore, because the fuel consumption of a compressor station is very small compared with the total flow rate for a gas pipeline, this assumption creates only a small error in the isothermal simulation.
2.1. Objective function The objective of the model is to minimize the fuel consumption of a gas pipeline during an operating period Tend:
E¼
Tend XNCS Z i¼1
Fi ðQstd ðtÞ; Ps ðtÞ; Ts ðtÞ; Pd ðtÞÞdt;
(1)
t¼0
where Ncs is the number of compressor stations and Fi the fuel consumption of the i-th compressor station, which is related to its flow rate Qstd, suction pressure Ps, suction temperature Ts, and discharge pressure Pd. According to the trapezoidal rule of numerical integration, Eq. (1) can be approximated as follows:
E¼
XNcs XNtime i¼1
n¼1
0:5fFi ½ðn 1ÞDt þ Fi ðnDtÞgDt;
Fi ðnDtÞ ¼ Fi ½Qstd ðnDtÞ; Ps ðnDtÞ; Ts ðnDtÞ; Pd ðnDtÞ;
(2) (3)
where Dt is the time step and Ntime is the number of time steps. Hence, Ntime*Dt is equal to the operating period, Tend. As previously stated, the fuel consumption of a station F is computed at the station level. Although it is assumed that gas flows
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
isothermally in a pipeline, thermal calculations are carried out at the station level. Before gas leaves a station, it is cooled to a prescribed temperature. Moreover, only steady-state operation is considered at the station level because a compressor reacts much faster than a pipeline does if its operating condition is changed. Because identical compressors are arranged in parallel in each compressor station for the pipeline addressed in this paper, only the number of online compressors is adjusted to minimize the fuel consumption at the station. This approach is taken because in this type of station configuration, distributing the load among the online compressors equally is optimal once the number of compressors to run is determined (Wright, 1998; Carter, 1996; Mahmoudimehr and Sanaye, 2013; Zhang et al., 2014).
2.2. Constraints The model constraints arise mainly from simulating a compressor unit and from the transient dynamics of gas flow in a pipeline segment. The limits of the compressor speed, driver power, and operating pressure are also taken into consideration.
2.2.1. Simulating a compressor unit The operating envelope of a centrifugal compressor is defined by its maximum speed, minimum speed, surge line, and stone line, as illustrated in Fig. 1, Eq. (4), and Eq. (5). In Eq. (4), S is the compressor speed, and Smin and Smax are the minimum and maximum speeds, respectively. In Eq. (5), Qac represents the compressor flow rate under its suction pressure and suction temperature, and the coefficients a0, a1, a2, b0, b1, and b2 are compressor-specific. The coefficients are calculated by regressing the surge line and the stone line of the compressor.
Smin S Smax
(4)
a0 þ a1 S þ a2 S2 Qac b0 þ b1 S þ b2 S2
(5)
The performance map of the compressor is regressed using Eqs. (6) and (7) (Odom and Muster, 2009; Haijun et al., 2015). In Eq. (6), H is the polytropic head of the compressor, and in Eq. (7), hc denotes the compressor's efficiency. The coefficients c0,0, c0,1, …, c2,2, d0, d1, d2, and d3 are also compressor-specific and are computed by regressing the performance map of the compressor:
H S2 ¼ c0;0 þ c0;1 S þ c0;2 S2 þ c1;0 þ c1;1 S þ c1;2 S2 ðQac =SÞ þ c2;0 þ c2;1 S þ c2;2 S2 ðQac =SÞ2 hc ¼ d0 þ d1 ðQac =SÞ þ d2 ðQac =SÞ2 þ d3 ðQac =SÞ3
195
(6)
(7)
Compressing gas with a compressor is regarded as a polytropic process and expressed as shown in Eq. (8), where mV is the polytropic exponent. According to the compressor suction conditions and the compressor's discharge conditions, Eq. (9) calculates the polytropic head, where Zs is the gas compressibility factor under the compressor suction conditions, R is the universal gas constant, and MGAS is the molar mass of the gas. The polytropic efficiency of the compressor is computed using Eq. (10), where kave is the average specific capacity ratio, which is calculated using Eq. (11) according to the specific capacity ratio at the compressor suction, ks, and at the compressor discharge, kd:
mV 1 Td Zs Pd mV ¼ Ts Zd Ps H¼
mV 1 mV R Zs Ts ðPd =Ps Þ mV 1 mV 1 MGAS
hc ¼
kave 1 P T lg d lg d kave Ps Ts
kave ¼ ðks þ kd Þ=2
(8)
(9)
(10)
(11)
To simulate a compressor, the system of equations composed of Eqs. (6)e(11) and an equation of state are solved. The compressor input power is calculated using Eq. (12), where m denotes the compressor mass flow rate and hm the mechanical efficiency, which is assumed to be constant in most cases. The fuel consumption of the compressor unit is computed using Eq. (13), where hGT is the driver efficiency and LHV the lower heating value of the fuel. Note that only compressors driven by gas turbines are considered in this paper.
Pshaft ¼ ðHmÞ=ðhc hm Þ
(12)
. f ¼ Pshaft LHV=hGT
(13)
The efficiency performance map of a gas turbine is regressed using Eqs. (14)e(18). In these equations, Ta and Pa are the atmospheric temperature and pressure onsite, respectively, and Ta,0 and Pa,0 are the design atmospheric temperature and pressure, respectively. In addition, the maximum available power of the gas turbine can be computed using Eq. (19). This value constrains the input power of the compressor, as illustrated by Eq. (20). The coefficients e0, e1, … e5, f0, f1, … f5 in these equations are gas-turbinespecific and are computed by regressing the performance map of the turbine.
hGT ¼ e0 þ e1 nPT þ e2 NPT þ e3 ðnPT Þ2
2 þe4 NPT þ e5 nPT NPT nPT ¼ S Fig. 1. Operating envelope of a compressor.
.pffiffiffi q
(14)
(15)
196
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
q ¼ Ta Ta;0
(16)
first reformulated. The proposed approach is then described in detail.
. pffiffiffi N PT ¼ Pshaft d q
(17)
3.1. Model reformulation
d ¼ Pa Pa;0
(18)
Pmax ¼ f0 þ f1 S þ f2 Ta þ f3 STa þ f4 S2 þ f5 Ta2
(19)
Pshaft Pmax
(20)
2.2.2. Governing equations of transient gas flow in a pipe Eqs. (21) and (22) describe the transient dynamics of gas flow in a pipeline segment. In Eq. (21), r is the gas density and v the gas velocity. In Eq. (22), g is the acceleration due to gravity, P is the gas pressure, h is the pipe height, and d is the pipe's inner diameter. In addition, l is the friction coefficient, which is estimated using Eq. (23) according to the absolute roughness of the internal wall, k, and the Reynolds number, Re:
vr vðrvÞ þ ¼0 vt vx
(21)
vðrvÞ v rv2 dh vP l v2 þ r ¼ gr vt dx vx d 2 vx
(22)
1 k 2:51 pffiffiffi ¼ 2lg þ pffiffiffi 3:7d Re l l
(23)
2.2.3. Pressure limits The discharge pressure of a compressor station is constrained by the maximum allowable operating pressure (MAOP) of the pipeline:
Pd PMAOP ;
(24)
and the inlet pressure of a compressor station is constrained by its minimum allowable pressure:
Pi Pi;min ; i ¼ 1; 2; …; NCS
(25)
In addition, the pressure at each delivery point should not fall below the contract pressure:
P Pcontract
(26)
2.2.4. Flow-rate balance At each node, Kirchhoff's first law is applied:
X
in ¼ Qiso
X
out Qiso
(27)
As previously mentioned, the fuel consumption of each compressor station is ignored in the equation. 3. Dynamic programming approach A dynamic programming approach is described in this section to solve the model formulated in Section 2. Because dynamic programming models have a special structure, the previous model is
The model described in the previous section is reformulated to fit the dynamic programming approach. Elements of the reformulated model described in this section include the definition of a stage, the choice of state variables and decision variables, and the index function. 3.1.1. Stage The dynamic programming recursion process starts from the pipeline destination and progresses to its origin. In this procedure, each stage addresses the zone from the discharge of a downstream station to that of an adjacent upstream station. The gas pipeline terminal is regarded as a station without any compressor, and its origin is considered to be a station with only one feasible discharge pressure, i.e., the source pressure. Hence, there are a total of NCS þ 1 stages in the dynamic programming approach. 3.1.2. State variables In the DP model, the station discharge pressures are chosen as the state variables. These pressures should vary over time because transient states are considered in this formulation. However, the pressures are assumed to be constant because multiple simulations have revealed that a gas pipeline can operate with constant station discharge pressures if its consumer demands vary modestly over time. This assumption makes it possible to use the dynamic programming method, which is fast, easy, and robust, to minimize the fuel consumption of a gas pipeline in transient states. Indeed, the assumption may make the DP approach fail in severe transient situations. However, these situations are rare in the operation of trunk gas pipelines. Therefore, failure under these conditions should be considered a minor defect of the approach. The contract pressure at the pipeline terminal is regarded as the initial state. Note that the initial state can be arbitrary and that the choice does not affect the algorithm results. This result derives from the fact that the pipeline terminal is considered a station without any compressor, as previously mentioned. Consequently, only the suction pressure is of concern. The possible discharge pressure of the i-th station is [Piþ1,min, PMAOP], which is continuous. However, only finite states can be considered in a DP approach. Hence, the region is discretized using a fixed step size, 0.1 MPa. The resulting discretized set is further refined by the procedure described in the following subsection. Note that the lower bound for the last compressor station is the contract pressure of the terminal consumer, and only the source pressure is considered feasible for the pipeline origin, as previously described. 3.1.3. Decision variables and state transition For simplicity, in each stage, the upstream station is referred to as CSA and the adjacent downstream station is referred to as CSB. At the beginning of each stage, a set of possible CSB discharge pressures is determined. Corresponding to each discharge pressure, there is a flow-rate trend record of CSB. For the initial state, there is only one record, which is the terminal consumer demand. In each stage, the system state transitions from the CSA discharge pressure to that of CSB by selecting a discharge pressure value in the set determined at the beginning of the stage. Hence, the decision variables are also the station discharge pressures. 3.1.4. Stage effect Once the CSB discharge pressure has been determined, the
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
transient gas flow in the pipeline segment connecting CSA and CSB is simulated. In the simulation, the upstream boundary condition is the CSA discharge pressure and the downstream boundary condition is the flow-rate trend record corresponding to the CSB discharge pressure. In this paper, the transient simulation is carried out using Stoner Pipeline Simulator (SPS) 9.5 because of its widespread use in the gas pipeline industry, good agreement with field data, flexibility and robustness. SPS is also the only available simulator we could invoke from our optimization codes. The transient simulation yields the CSA flow-rate trend and the CSB suction pressure trend. The latter trend is examined to determine whether it falls below the minimum allowable CSB suction pressure at any time. If it does, another CSB discharge pressure is chosen, and the simulation process is repeated. If all possible CSB discharge pressures have been tested and the minimum allowable CSB suction pressure cannot be fulfilled in any of the cases, the CSA discharge pressure is considered infeasible, as is any other lower pressure value. Hence, the computational effort is reduced. If the CSB suction pressure trend stays above its minimum allowable suction pressure, the station level is invoked. This level examines whether CSB can operate at the given suction pressure, flow rate, and discharge pressure. If it can, its minimum fuel consumption during the operating period is computed and returned to the pipeline level. This fuel consumption is the stage effect. In addition, the CSA flow-rate trend is also recorded with its discharge pressure for future use. If CSB cannot operate under the given conditions, another CSB discharge pressure is chosen, and the process is repeated. Simulating transient gas-flow dynamics in a pipeline segment is computationally demanding, and numerous transient simulations are conducted in the DP approach. This process makes the computational time of this approach very long. However, these simulations are essential for obtaining a satisfactory and reliable solution. Coarse approximation may speed up the approach but at the risk of providing impractical results. Moreover, the DP approach is accelerated by detecting infeasible station discharge pressures, as previously described. The approach also benefits from selecting the best way to transition to a target state and considering previous decisions irrelevant once the transition is complete. 3.1.5. Index function The index function of the DP model is
F Pd;CSA ¼
min f Pd;CSA ; Pd;CSB þ F Pd;CSB ; Pd;CSB 2SðPd;CSB Þ
(28)
where S(Pd,CSB) is the set of possible CSB discharge pressures and f(Pd,CSA, Pd,CSB) is the CSB fuel consumption. If an infeasible situation occurs, f is set to a very large number. For the initial state, F(Pd,CSB) ¼ 0. 3.1.6. Markov property One important quality of a DP model is its Markov property. This property is related to the choice of state variables, which are the station discharge pressures in the current approach. In each stage, a transient simulation is conducted once the CSB discharge pressure is determined, and the simulation results are not related to previous decision procedures. Thus, the Markov property of the DP model is guaranteed. 3.2. Algorithm procedures Based on the reformulated model, the proposed algorithm is described in detail in the following steps:
197
0. Set k ¼ Ncs, Sstate ¼ {Pcontract}, Sobjfun ¼ {0}, and Scsdisflow ¼ {qenddelivery}. 1. Discretize the discharge pressure region [Pkþ1,min, PMAOP] with step size 0.1 MPa, and generate a set of discrete points Sdecision ¼ {P1, P2 … PNdecision}. 2. Sstate,pre ¼ Sstate, Sobjfun,pre ¼ Sobjfun, Scsdisflow,pre ¼ Scsdisflow, Sstate ¼ F, Sobjfun ¼ F, and Scsdisflow ¼ F. 3. Set i ¼ 1. 4. Set j ¼ 1. 5. Calculate the flow rate into the (k þ 1)-th station from the upstream pipeline by a flow-rate balance. Its discharge flow rate is Qj2Scsdisflow,pre. 6. Simulate the transient gas-flow dynamics in the pipeline segment that connects the k-th and the (k þ 1)-th stations. In the simulation, the upstream boundary condition is the constant suction pressure Pi2Sdecision, and the downstream boundary condition is the flow rate determined in Step 5. This transient simulation yields the discharge flow-rate trend of the k-th station and the suction pressure trend of the (k þ 1)-th station. 7. If the suction pressure of the (k þ 1)-th station stays above its minimum allowable suction pressure, go to 8. Otherwise, go to 11. 8. Examine whether the (k þ 1)-th station can operate at the given suction pressure, flow rate, and discharge pressure. If it can, calculate its fuel consumption, Fkþ1, and go to 9. If it cannot, go to 11. 9. If j ¼ 1, add Pi to Sstate, and record objfunj þ Fkþ1 as objfuni in Sobjfun, where objfunj is the fuel consumption of the previous stages corresponding to Qj and is stored in Sobjfun,pre. Record the discharge flow-rate trend of the k-th station in Scsdisflow. If j > 1, go to 10. 10. If objfunj þ Fkþ1 < objfuni, update objfuni and the corresponding discharge flow-rate trend stored in Scsdisflow. If objfunj þ Fkþ1 >¼ objfuni, go to 11. 11. j ¼ jþ1; if j > jSstate,prej, go to 12. Otherwise, go to 5. 12. If the suction pressure of the (k þ 1)-th station falls below its minimum allowable suction pressure for all of the flow-rate trends stored in Scsdislfow,pre, then Pi is considered infeasible and is eliminated from Sdecision together with all other lower pressure values. Otherwise, go to 13. 13. i ¼ iþ1; if i > jSdecisionj, go to 14. Otherwise, go to 4. 14. k ¼ k1; if k 0, go to 1. Otherwise, go to 15. 15. The minimum value in Sobjfun corresponds to the optimal solution.
3.3. Algorithm acceleration The forward method often requires more than 1 h to compute an optimal solution according to algorithm tests. We devise a preprocessing method to cut the computational time. By assuming that the gas pipeline is operated in a steady state at each time step, we approximate a transient optimization problem as a series of steadystate optimization problems. Then, each of them is solved. The minimum discharge pressure Pk,min,ss of each compressor station is searched among the solutions, as is the maximum pressure Pk,max,ss. Finally, the discharge pressure region [Pk,min,ss, Pk,max,ss] replaces [Pkþ1,min, PMAOP] in the DP algorithm. The former is typically much smaller than the latter. Thus, the algorithm can be sped up. To approximate a transient optimization problem as a series of steady-state optimization problems, the time step should be small. In the following tests, a time step of 1 h is used. In addition, for the last compressor station, Pmax,ss is set to min(Pmax,ss þ 500, PMAOP) in case the pressure at the pipeline end falls below the contract
198
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
pressure. 4. Algorithm tests The proposed DP approach is first tested on the Central AsiaChina gas pipeline, illustrated in Fig. 2. This pipeline extends across more than 1800 km and transports 30 billion cubic meters of natural gas per year from Turkmenistan and Kazakhstan to China. The pipeline is composed of two identical parallel pipelines, eight compressor stations, two sources, and four delivery points. Identical compressors driven by gas turbines are installed at each compressor station. Table 1 presents the number of compressors at each station and the speed limits of these compressors. Note that compressors in distinct stations may be different. The maximum discharge pressure of each station is 9810 kPa, and the minimum suction pressure is 6000 kPa. The contract pressure at C4 is 7500 kPa. To examine the effect of the key assumption made in Section 3, the proposed approach is used to minimize the fuel consumption of the last section of the gas pipeline. The result is compared with the true optimum to assess the effect of the assumption. Multiple tests are carried out on the pipeline for the original DP approach and the accelerated one. Two of these tests are described in detail below. The operating periods of these two tests are 24 h and 168 h. Finally, the accelerated DP approach is adopted to solve the problem studied by Geißler et al. (2011) and Domschke et al. (2011). This test involves a more complex network.
Table 1 Compressor station details. Station
Number of compressors
Minimum speed
Maximum speed
CS CS CS CS CS CS CS CS
3 3 3 3 3 4 4 4
3965 3965 3965 3965 3120 3965 3120 3120
6405 6405 6405 6405 5040 6405 5040 5040
1 2 3 4 5 6 7 8
4.1. Influence of the time-independent assumption To examine the effect of the assumption that the discharge pressure of each station is time-independent, the proposed approach is used to optimize the operation of the last section of the gas pipeline in transient states. In the problem, the suction pressure of the last compressor station and the terminal consumer demand are given. The fuel consumption of the last station is minimized over a certain period. This problem is chosen because its true optimum is available (Wong and Larson, 1968; Goldberg, 1983), and the method proposed by Wong and Larson (1968) is used to compute this true optimum. However, to make the comparison fair, the compressor station discharge pressure is chosen as the decision variable instead of its flow rate, as in the original paper, and SPS is used to perform transient simulations. A typical demand at the pipeline destination is plotted in Fig. 3. For the test problem, the optimal fuel consumption is 23.89 104 Sm3/d, whereas that calculated by the proposed method is 25.37 104 Sm3/d, or 6.2% higher. The optimum solution is depicted in Fig. 4 together with the solution provided by the proposed algorithm, which is close to the average of the optimal solution. The pressure at C4 is plotted in Fig. 5. For the true optimum, the pressure at C4 remains slightly above its lower bound except during the first few time steps. This result shows that the pressure energy is fully used to transport the gas. In contrast, for the proposed algorithm, the pressure at C4 is somewhat higher than
Fig. 3. Demand at C4.
Fig. 4. Optimal solution.
the contract pressure, which means that part of the pressure energy is wasted. Although a definitive conclusion cannot be drawn from a single case study, the study can at least shed some light on this issue. Because the state of a gas pipeline rarely fluctuates severely, it is believed that the solution computed by the proposed algorithm is not too far from the true optimum that in most cases arises during daily operations. Fig. 2. Topological structure of the Central Asia-China gas pipeline.
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
199
Fig. 7. Demand at C4. Fig. 5. Pressure at C4.
4.2. 24-hour case In the current case, S1 supplies gas at a constant pressure (7000 kPa) and S2 supplies gas at a constant flow rate (41.667 104 Sm3/h). Fig. 6 shows the demands of intermediate consumers, and the demand at C4 is plotted in Fig. 7. Most of the gas is delivered to the pipeline terminal, and the demand at C2 fluctuates severely within the operating period. Note that the flow rates considered in this paper are flow-rate values under standard conditions of 101.325 kPa and 20 C, as used in the Chinese petroleum industry. Table 2 shows the optimal solution computed by the original DP method. The optimal discharge pressures of most stations are equal to the maximum allowable operating pressure. Note that the design flow rate of the pipeline is approximately 357 104 Sm3/h, whereas the total flow rate in this case study is 378.04 104 Sm3/h. The pipeline system operates under a heavy load. Hence, the discharge pressure of most compressor stations in the solution is equal to the maximum allowable pressure. Fig. 8 shows the discharge flow rate of each compressor station, and their suction pressures are plotted in Fig. 9. The figures show that CS 5 and CS 6, which are located adjacent to C2, are most
Table 2 Optimum discharge pressure. Station
Discharge pressure
CS CS CS CS CS CS CS CS
9810 9810 9810 9810 9810 9810 9310 9110
1 2 3 4 5 6 7 8
Fig. 8. Discharge flow rate of each compressor station.
Fig. 6. Demands of intermediate consumers.
strongly affected by their demand fluctuation. With increasing distance, the effect diminishes gradually. In addition, the suction pressure of each compressor station remains far above the minimum allowable suction pressure (6000 kPa), as illustrated in Fig. 9. The pressure at C4 also remains above the contract pressure, 7500 kPa, during the operating period, as shown in Fig. 10. However, the pressure is not far from the lower bound. This situation is the result of efforts to minimize fuel consumption because higher pressure indicates extra compression power. The station level is also examined. CS 5 and CS 6 are investigated carefully because their operating conditions fluctuate severely. Two compressors are online in CS 5 during the operating period. Because identical compressors are installed in CS 5, the speed and input power of only one compressor are plotted in Fig. 11. The figure
200
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
Fig. 9. Suction pressure of each compressor station.
Fig. 12. Speed and fuel consumption of a compressor in CS 6.
plotted in Fig. 12. The figure shows that the compressor speed is flexible and that the input power is no greater than 16,600 kW, whereas the maximum available power of its driver is never less than 27,050 kW. Table 3 shows the solution computed by the accelerated algorithm. The solution is nearly the same as that computed by the original approach. After acceleration, the objective function value increases only by 0.72%, and the computational time is reduced by 92.29% e464 s. Twenty-four more tests are carried out. For twenty tests, the accelerated algorithm can compute an optimal solution within 15 min. The average computational time is 6.6 min. For the four remaining tests, the computational time is set to 1 h. The average computational time for each test is 15.2 min. In addition, the objective function value of the optimal solutions increases by 1.22% on average after acceleration. Fig. 10. Pressure at C4.
4.3. 168-hour case In the current case, all of the gas is delivered to C4, and Fig. 13 shows the corresponding demand, which fluctuates dramatically during the operating period. The sources supply gas under the same conditions as those in the 24-h case. The optimal solution yielded by the original DP approach is presented in Table 4. The discharge pressure of each station is equal to the maximum allowable operating pressure, except for CS 8. This result indicates that the pipeline capacity is nearly reached. Fig. 14 displays the discharge flow rate of each station, and the suction pressure trends are plotted in Fig. 15. CS 8 and CS 7 are markedly affected by the demand fluctuation at C4. With increasing distance, the influence diminishes gradually. Moreover, the suction pressure of each station remains well above the minimum allowable suction pressure, 6000 kPa, as shown in Fig. 15. The pressure at C4 is illustrated in Fig. 16. Although the pressure
Fig. 11. Speed and fuel consumption of a compressor in CS 5.
shows that the compressor speed varies within the speed limits and that the input power remains below 14,500 kW. Note that the maximum available power of a gas turbine varies as its speed changes. However, this available power for the current case is never less than 24,303 kW, which is far more power than required. Two compressors are also online in CS 6 during the operating period. The compressor speed and input power of one compressor are
Table 3 Optimum discharge pressure. Station
Discharge pressure
CS CS CS CS CS CS CS CS
9810 9810 9810 9810 9810 9718 9526 9334
1 2 3 4 5 6 7 8
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
201
Fig. 15. Suction pressure of the compressor stations.
Fig. 13. Demand at C4.
Table 4 Optimum discharge pressure. Station
Discharge pressure
CS CS CS CS CS CS CS CS
9810 9810 9810 9810 9810 9810 9810 9710
1 2 3 4 5 6 7 8
Fig. 16. Pressure at C4.
Fig. 14. Discharge flow rates of the compressor stations.
fluctuates significantly, it still remains above the contract pressure, 7500 kPa. The station level is also examined carefully, and no constraint is found to be violated. However, because the minimum up time and the minimum down time of a compressor unit are not taken into consideration in the current model, some compressors may run for a short period and then be turned off. This problem is found in CS 8, as shown in Fig. 17, and should be addressed in future research. Table 5 presents the optimal solution computed by the accelerated algorithm. The solution is nearly the same as that offered by the original approach. The accelerated algorithm computes this solution in approximately 32 min e 10.2% of the computational time of the original algorithm. However, the objective function
Fig. 17. Number of operating compressors in CS 8.
value increases by only 0.42%. Although the computational time is relatively long, it is considered acceptable considering the long operating period of the problem and the details of the optimization model. In addition, the computational time increases modestly as the operating period
202
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203 Table 5 Optimum discharge pressure.
Table 6 Pipe details.
Station
Discharge pressure
Pipe
Length (km)
Outer diameter (mm)
Thickness (mm)
CS CS CS CS CS CS CS CS
9810 9810 9810 9810 9810 9810 9716 9616
PL1 PL2 PL3 PL4 PL5 PL6 PL7 PL8 PL9 PL10 PL11 PL12
100 13 182 159 70 136 110 165 101 45 200 220
1016 1016 1016 1016 457 762 660 762 660 660 660 660
21 21 21 21 17.5 19 19 19 19 19 19 19
1 2 3 4 5 6 7 8
becomes longer. Two factors contribute to this effect. First, the assumption that the discharge pressure of each compressor station is time-independent makes problems with a longer operating period more constrained, leading to fewer feasible solutions. Second, the computational time for a transient simulation does not grow linearly as the operating period becomes longer.
4.4. Test on a more complex network In this section, we use the accelerated algorithm to solve the problem studied by Geißler et al. (2011) and Domschke et al. (2011). Fig. 18 shows the pipeline network. The transient process lasts for 4 h. Source1 supplies gas at a constant pressure of 7000 kPa, and source2 supplies gas at a constant flow rate of 36 104 Sm3/h. The flow rate of all of the deliveries except for sink A04 increases linearly from 36 104 Sm3/h to 42.235 104 Sm3/h. Sink A04 delivers gas at a constant flow rate of 18.5 104 Sm3/h. In their respective studies, Geißler et al. (2011) and Domschke et al. (2011) offered only part of the network details. To conduct the test, assumptions are made about the missing information. Network elevation change is ignored. Table 6 shows the pipe details. The roughness of all of the pipes is 0.01 mm. The maximum discharge pressure of each compressor station is 10,000 kPa, the minimum suction pressure of each station is 6000 kPa, and the contract pressure of each delivery is 4000 kPa. With the preceding assumptions, the accelerated algorithm solves this problem within 10 min. The optimum discharge pressures of CS 1 and CS 2 are 10,000 kPa and 9500 kPa, respectively, and CS 3 is bypassed.
5. Conclusions Based on the preceding analyses and tests, we draw the following conclusions: 1. The proposed approach was used to minimize the fuel consumption of the last section of a gas pipeline. Its result was only 6.2% higher than the true optimal. 2. The proposed approach was also tested on various problems arising in daily gas pipeline operations. For most case studies, the approach can yield a solution within 15 min.
Fig. 18. Network topology.
3. The computational effort of the proposed algorithm grows modestly as the operating period becomes longer. To increase the applicability of the proposed algorithm, we suggest carrying out the following work in the future: 1. Consider the minimum up time and the minimum down time of a compressor unit at the station level. 2. The SPS simulator we implemented can only use a portion of the power of multi-core CPUs. These CPUs are common today, and full use of all of the processors can greatly reduce the computational time. Therefore, we suggest replacing the simulator with one that is capable of conducting parallel computations. Acknowledgments This research was supported by the Australian and Western Australian Governments and the North West Shelf Joint Venture Partners, as well as the Western Australian Energy Research Alliance through the Australia-China Natural Gas Technology Partnership Fund Top Up Scholarship. This study was also supported by the Beijing Higher Education Young Elite Teacher Project (No. YETP0677). We also thank the anonymous reviewer for his/her insightful comments, which helped us improve our manuscript. Nomenclature d Fe fe ge He he ke ks, kd e
inner diameter of a pipe, m fuel consumption of a compressor station, Nm3/h fuel consumption of a compressor, Nm3/h gravity, 9.81 m/s2 polytropic head of a compressor, kJ/kg height, m absolute roughness of an internal pipe wall, m heat capacity ratio of natural gas at the suction and discharge of a compressor LHV e lower heating value of the fuel, kJ/Nm3 mV e polytropic exponent MGAS e molecular mass, kg/kmol Ncs e number of compressor stations Ntime e number of time steps P0 e ISO pressure, 101.325 kPa Pa e ambient pressure at a site, kPa Pa,0 e ambient pressure of a design, 101.325 kPa Pd e discharge pressure of a compressor or a station, kPa Pmax e maximum output power of a gas turbine, kW Ps e suction pressure of a compressor or a station, kPa Pshaft e input power of a compressor, kW Qac e volumetric flow rate under specific conditions, m3/h
X. Zhang et al. / Journal of Natural Gas Science and Engineering 28 (2016) 193e203
Qiso e volumetric flow rate under standard conditions, Nm3/h Re e reynolds number Re gas constant, 8.314 kJ/(kg*K) Se compressor speed, RPM Scsdisflow set of possible discharge flow-rate patterns of a station Sdecision e set of possible decisions Sendflow e flow rate at the end of a pipeline Sobjfun e set of fuel consumption records of previous stages te time, s Te optimization horizon, min T0 e ISO temperature, 20 C Ta e ambient temperature at a site, K Ta,0 e ambient design temperature, 288.15 K Td e discharge temperature of a compressor or a station, C Ts e suction temperature of a compressor or a station, C we gas velocity, m/s xe length, m Zs e compressibility factor of natural gas at the station suction Dt e time step, min hC e compressor efficiency, [0,1] hM e mechanical efficiency, 0.98 hGT e gas turbine efficiency, [0,1] re gas density, kg/m3 le friction coefficient of a pipe References Alfred, S., Fasullo, J., Pfister, J., Daniels, A., 2013, April. Capacity determination using state finding and gas transient optimization. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Anglard, P., 1988, January. Hierarchical steady state optimization of very large gas pipelines. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Carter, R.G., 1996, January. Compressor station optimization: computational accuracy and speed. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Carter, R.G., 1998, October. Pipeline optimization: dynamic programming after 30 years. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Carter, R., Goodreau, M., Rachford, H., 2001. Optimizing pipeline operations through mathematical advances. Pipeline Gas J. 228 (10), 51e53. Carter, R.G., Rachford Jr., H.H., 2003, January. Optimizing line-pack management to hedge against future load uncertainty. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Carter, R.G., Rachford, H.H., Dupont, T.F., 2006, January. A new approach to gas control simulation. In: 2006 International Pipeline Conference. American Society of Mechanical Engineers, pp. 865e875. Carter, R., Reisner, M., Sekirnjak, E., 2010, January. Transient optimization-examples and directions. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Casoetto, B., Flottes, E., Ardeois, J., Thanh, N., Joliot, J.B., Nait-Abdallah, R., 2011. How to commercialize reliable capacities on a complex transmission network? J. Nat. Gas Sci. Eng. 3 (5), 657e663. Chebouba, A., Yalaoui, F., Smati, A., Amodeo, L., Younsi, K., Tairi, A., 2009. Optimization of natural gas pipeline transportation using ant colony optimization. Comput. Oper. Res. 36 (6), 1916e1923. Domschke, P., Geißler, B., Kolb, O., Lang, J., Martin, A., Morsi, A., 2011. Combination of nonlinear and linear optimization of transient gas networks. Inf. J. Comput. 23 (4), 605e617. Dupont, T.F., Rachford Jr., H.H., 1987. Optimization of power usage in transient gas transmission lines. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Fasihizadeh, M., Sefti, M.V., Torbati, H.M., 2014. Improving gas transmission networks operation using simulation algorithms: case study of the national Iranian gas network. J. Nat. Gas Sci. Eng. 20, 319e327. Ferber, P., Basu, U., Venkataramanan, G., Goodreau, M., Linden, P., 1999, January. Gas pipeline optimization. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Geißler, B., Kolb, O., Lang, J., Leugering, G., Martin, A., Morsi, A., 2011. Mixed integer linear models for the optimization of dynamical transport networks. Math. Methods Oper. Res. 73 (3), 339e362.
203
Goldberg, D.E., 1983. Computer-Aided Gas Pipeline Operation Using Genetic Algorithms and Rule Learning. Doctoral Dissertation. University of Michigan. Haijun, C., Changchun, W., Xiaorui, Z., 2015, January. Corrected maps narrow compressor performance prediction range. Oil Gas J. 113 (1), 72e79. Kelling, C., Reith, K., Sekirnjak, E., 2000, October. A practical approach to transient optimization for gas networks. In: PSIG Meeting, vol. 5. Krishnaswami, P., Chapman, K.S., Abbaspour, M., 2004, January. Compressor station optimization for linepack maintenance. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Liu, E., Li, C., Yang, Y., 2014. Optimal energy consumption analysis of natural gas pipeline. Sci. World J. 2014. Lall, H.S., Percell, P.B., 1990. A dynamic programming based gas pipeline optimizer. In: Analysis and Optimization of Systes. Springer Berlin Heidelberg, pp. 123e132. Mahlke, D., Martin, A., Moritz, S., 2007. A simulated annealing algorithm for transient optimization in gas networks. Math. Methods Oper. Res. 66 (1), 99e115. Mahmoudimehr, J., Sanaye, S., 2013. Minimization of fuel consumption of natural gas compressor stations with similar and dissimilar turbo-compressor units. J. Energy Eng. 140 (1). Mantri, V.B., Prestion, L.B., Pringle, C.S., 1985. Transient optimization of a natural gas pipeline system. In: PSIG Annual Meeting. MohamadiBaghmolaei, M., Mahmoudy, M., Jafari, D., MohamadiBaghmolaei, R., Tabkhi, F., 2014. Assessing and optimization of pipeline system performance using intelligent systems. J. Nat. Gas Sci. Eng. 18, 64e76. Moritz, S., 2007. A Mixed Integer Approach for the Transient Case of Gas Network Optimization. Doctoral Dissertation. TU Darmstadt. Odom, F.M., Muster, G.L., 2009, January. Tutorial on modeling of gas turbine driven centrifugal compressors. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Peretti, A., Toth, P., 1982. Optimization of a pipe-line for the natural gas transportation. Eur. J. Oper. Res. 11 (3), 247e254. Pietsch, U., Rachford Jr., H.H., Carter, R.G., 2001, January. Investigating real-world applications of transient optimization. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Rachford Jr., H.H., Carter, R.G., 2000, January. Optimizing pipeline control in transient gas flow. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Rachford Jr., H.R., Carter, R.G., 2004. U.S. Patent No. 6,701,223. U.S. Patent and Trademark Office, Washington, DC. Rachford Jr., H.H., Carter, R.G., Dupont, T.F., 2009, January. Using optimization in transient gas transmission. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Rossi, B., Casoetto, B., Le Maitre, A., Jacquiau, A., Moumini, N., 2012, January. Computing and optimizing the available linepack in a gas network to address the issue of flexibility. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Rossi, B., Dalphin, J., Sinegre, L., Jacquiau, A., Renaudie, T., 2011, January. Dealing with CCGT: an explicit dynamic model. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Rossi, B., Djerourou, F., Suez-Crigen, G.D.F., 2014, May. Transient optimization in gas transmission networks, a new approach on GRTgaz network. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Sekirnjak, E., 1996, October. Practical experiences with various optimization techniques for gas transmission and distribution systems. In: Proceedings of the 28th Annual Meeting of the Pipeline Simulation Interest Group, pp. 1e16 (PSIG'96). Schmidt, M., Steinbach, M.C., Willert, B.M., 2014. High Detail Stationary Optimization Models for Gas Networks: Validation and Results. Friedrich-AlexanderUniversit€ at Erlangen-Nürnberg, Department Mathematik. Tech. rep. Submitted. Thanh, N., le Maitre, A., Ardeois, J., Joliot, J.B., 2013, April. An optimization tool for network capacities. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Wong, P., Larson, R., 1968. Optimization of natural-gas pipeline systems via dynamic programming. Automatic control. IEEE Trans. 13 (5), 475e481. Wright, S., Ditzel, C., Somani, M., 1998, January. Compressor station optimization. In: PSIG Annual Meeting. Pipeline Simulation Interest Group. Wu, X., Li, C., Jia, W., He, Y., 2014. Optimal operation of trunk natural gas pipelines via an inertia-adaptive particle swarm optimization algorithm. J. Nat. Gas Sci. Eng. 21, 10e18. Zhang, X., Wu, C., Zuo, L., Meng, X., 2014, November. Dynamic programming based algorithm for compressor station optimization. In: ASME 2014 International Mechanical Engineering Congress and Exposition. American Society of Mechanical Engineers (Pp. V011T14A028-V011T14A028). Zheng, Z., Wu, C., 2012. Power optimization of gas pipelines via an improved particle swarm optimization algorithm. Petroleum Sci. 9 (1), 89e92.