Technical and economic optimization of natural gas transmission network operation to balance node delivery flow rate and operation cost

Many researches have been conducted in optimization of gas transmission network operation. Some published review literatures describe developments in this area. Rios-Mercado and Borraz-Sanchez [1] provides a review on the most relevant research works conducted to solve natural gas transportation problems via pipeline systems. Wu et al. [8] provides a review on the progress related to the steady-state operation optimization models of natural gas pipelines as well as corresponding solution methods based on stochastic optimization algorithms. Next, this article will introduce the related research progress of gas transmission network operation optimization from the aspects of mathematical model types, solving methods and objective functions.

Long before, researchers have recognized operation optimization models of gas transmission network as a nonlinear programing (NLP) model because of the nonlinear pipeline and compressor constraints [9]. Besides, if consider the RCs number in each station, the model then turns into a mixed-integer NLP (MINLP) model. Hence, the models can be classified into NLP and MINLP model depending on the decision variable. NLP models assume that the RCs number in the network is known in advance. Besides, discrete variables that represent the number of compressors to be turned on are introduced into MINLP models.

The solving method of gas transmission network operation optimization have been studied by many researchers. However, there is still existing a tremendous potential. This is because different methods have different advantages, disadvantages, and range of applications [10]. Dynamic programming (DP) is the first method applied to the operating optimization of gas network. In 1968, Wong and Larson [11] firstly applied DP to the minimum fuel cost problem. However, as the dimension of the model increases, the computation quantity of the DP algorithm increases exponentially. So, it is difficult to solve the problem of large pipeline optimization. Percell and Ryan [12] studied the generalized reduced gradient (GRG) for non-cyclic structures. The GRG method conquered the dimension problem of DP algorithm and is suitable for complex structural networks. Futhermore, the linear programming (LP) was applied to deal with this problem [13], by transforming the non-linear model into a linear model. The advantage of LP is that the calculation speed is usually faster. However, if the initial value is not chosen properly, the solution of the algorithm is easy to fall into a inferior local optimal solution.

Different from above traditional deterministic methods, some emergent metaheuristic optimization algorithms have shown many advantages in terms of dealing with the large-scale pipeline and MINLP problem [8]. In recent years, some studies applied metaheuristic algorithms in this area effectively, for instance, ant colony optimization [14], Genetic Algorithm [15, 16], Simulated Annealing [17], Particle Swarm Optimization [18], and their extensions.

Gas transmission networks are usually characterized as a complex and huge system. We are supposed to fully consider the economic, safety and social benefits of the plan in the process of developing the network operation plan. Therefore, operation optimization of the gas transmission network is essentially a multi-objective optimization problem. However, because of the difficulty and computational complexity of multi-objective optimization techniques, most previous scholars only conducted single-objective optimization studies. The main objective functions are: delivery flow rate maximization [16], line pack maximization [19, 20], power consumption minimization [21, 22], profits maximization [23], etc. Some common optimization objectives in natural gas network operation optimization are summarized in Fig. 2.

OPEN IN VIEWER

Fig. 2

Common objectives of gas transmission network operation optimization.

With the development of multi-objective optimization algorithms in recent years, some researchers have begun to apply multi-objective optimization techniques to gas transmission network operation optimization problems. Wu et al. [24] utilized weighted sum method and inertia-adaptive PSO algorithm to solve the bi-objective model considered maximum operation benefit and maximum flow rate. Alinia Kashani and Molaei [25] built a MINLP model to minimizing operating cost, maximizing throughput and maximizing line pack. The developed model was solved by NSGA-II. Demissie et al. [26] presented a multi-objective model for different pipeline configurations of linear, tree and looped topologies. They took energy consumption minimization and transport flow rate maximization as objective functions. Su et al. [27] developed a multi-objective optimization method to trade-off reliability and power demand of natural gas pipeline network. Nevertheless, there is still a big gap in the study of multi-objective optimization compared with single objective optimization.

The main contributions of this paper are summarized below. Firstly, previous studies generally applied a large number of assumptions and simplified pipeline and compressor station models, leading to the reduction of accuracy and practicability of the model. For instance, Wu et al. [24] simplified the compressor model without considering the working rule and working interval of compressors. Demissie et al. [26] supposed that the gas Compressibility factor is an average constant and did not consider its variety in different pipelines. In this paper, a relatively precise and realistic optimization model of compressor stations and pipelines is established. This model considers the working rule and working feasible area of the compressor. The gas density, friction factor and compressibility factor are accurately calculated, thus replacing the method of picking empirical constants directly. Besides, the previous optimization of natural gas network operations mostly belonged to single-objective optimization problems. Even for multi-objective optimization problems, most of them focused on maximizing the total pipeline output [24–26]. Unfortunately, few researches have considered maximizing the delivery flow rate for a specific consumer. Hence, this paper focuses on the delivery flow rate at a specific gas delivery node GDN. And an optimization model considering maximum delivery flow rate of specified GDN and minimum CS power consumption cost is built. Finally, given that previous research usually selected small pipeline examples with simple structure and few users [25, 26], few studies were conducted on real-world large pipelines, and failed to effectively combine multi-objective optimization with engineering practice. In view of this, this study selected a practical large-scale gas transmission network as the optimization object. This pipeline is 2,229 km in length, and contains 7 CSs, 2 GINs and 20 GDNs.

Consequently, the main purpose of this paper is to establish a multi-objective operation optimization model for large gas transmission network with considering two conflict objectives of minimum CS energy consumption cost and maximum delivery flow rate of a specific GDN. Then ɛ-constraint method and GAMS/DICOPT solver were adopted to solve the multi-objective optimization problem.

The rest of the article is structured as follows: Section 2 defines the system assumptions, bi-objective functions, constraints, and decision variables. Section 3 studies optimization method and solver. Section 4 discusses bi-objective optimization results and sensitivity analyses. Finally, Section 5 summarizes the full text.

2.1 System assumptions

(1) Steady state flow

When natural gas flow in pipelines, it generally in an unstable state, in which the hydraulic and thermal parameters of gas vary over time. However, when it comes to the annual or monthly operational planning of the gas transmission network system, due to the long time, the system instability has little impact on the research results and can be ignored. Therefore, the optimization problem of long-term natural gas network operation is usually simplified to a steady state.

(2) Isothermal flow

Generally, natural gas flow temperature changes along pipelines, and finally approximate the environment temperature where the pipeline is constructed. However, practical applications show that the gas temperature changing value is much smaller than the Kelvin temperature reference value (273.15 K), and the variability within a certain range has little impact on the fluid pressure and flow rate. Hence, this paper adopted the assumption of gas flow in pipelines isothermally, and the temperature is 293.15 K.

2.2 Multi-objective functions

(1) Maximize delivery flow rate of a specified GDN

Various gas supply ways may exist for a certain gas user. With the continuous increase of natural gas consumption in the future, gas transmission network, as the most important and economical means of transportation, should increase its gas transmission capacity as much as possible to meet natural gas demand of users. Hence, the delivery flow rate at the specified GDN maximization might be considered as an objective function, which can be modeled as Equation (1). max f 1 = Q id out id ∈ N id(1)

Where Q id out is the delivery flow rate at a specified GDN, m³/s.

(2) Minimize CS power consumption cost

Minimum operating cost is a major interests of the network company. Gas transmission network operating cost mainly includes CS power consumption cost, CS equipment maintenance and repair costs, pipeline’s maintenance cost, and so on. This paper mainly studies the variable operation cost of network without considering the maintenance and repair cost of equipment. The variable operation cost is mainly made up of the CS power consumption cost. Therefore, this paper chooses to minimize the CS power consumption cost as another optimization goal. The objective function is represented as Equation (2). min f 2 = ∑ j = 1 M j E j δ j W j t j j ∈ M j(2) Where E _j is the power price of the jth CS, ¥/(kW·h); δ _j is the RCs number of the jth CS; W _j is the single compressor power of the jth CS, kW; t _j is the compressor running time of the jth CS, h.

(1) Node flow balance equation

The node flow balance equation is derived from the law of mass conservation. At any node, that is for all GINs, GDNs and GTNs, the inlet flow and outlet flow should be equal. This constraint can be expressed as Equation (3). Q i in - Q i out + ∑ β ij Q j CS + ∑ β ik Q k p = 0 i ∈ N(3)

Where Q i in is the injection flow rate of the ith node, m³/s; Q i out is the injection flow rate of the ith node, m³/s; Q j CS is the flow rate of the jth CS connected to the ith node, m³/s; β _ij is the fluid direction coefficient, if the gas flows into ith node, β _ik  = 1, otherwise, β _ik  = - 1; Q k p denotes the flow passes through the kth pipeline connected to the ith node, m³/s; and β _ik has the same meaning as β _ij .

(2) Gas flow equation

When natural gas flows in the pipeline, part of the energy will be lost due to the roughness of the inner wall of the pipe, causing the pressure gradual loss along the pipeline. The general flow equation describes the relationship between the gas flow rate Q k p and inlet pressure P k a , outlet pressure P k b of the pipe, which is derived from one-dimensional momentum conservation equation [28]. The general flow equation is expressed as Equation (4). Q k p = C [ ( P k a ) 2 - ( P k b ) 2 ] D k 5 λ k Z k avg Δ TL k k ∈ M k(4)

Where Q k p is the flow rate of the kth pipeline, m³/s; P k a is the inlet pressure of the kth pipeline, Pa; P k b is the outlet pressure of the kth pipeline, Pa; C is the constant, 0.0384 (m²·s · K /kg); D _k is the diameter of the kth pipeline, m; λ _k is the friction coefficient of the kth pipeline; Z k avg is the average compressibility factor of the kth pipeline; Δ is the relative density of gas; T is gas flow temperature, K; L _k is the length of the kth pipeline, m.

Based on the inlet pressure P k a and outlet pressure P k b of pipeline, the average pressure P k avg can be calculated by Equation (5). P k avg = 2 3 ( P k a + ( P k b ) 2 P k a + P k b )(5)

(3) Friction factor equation

The friction factor λ is a significant gas flow parameter, which is correlated with the gas flow state and pipeline’s roughness. According to the Reynolds number and relative roughness of pipe, the gas flow state is divided into laminar flow and turbulent flow. While the turbulent flow is divided into three areas: hydraulic smooth zone, mixed friction zone and resistance square zone. Different equations have different scope of application. These equations mainly include Weymouth, Pan(A), Pan(B), AGA, Colebrook, etc. Colebrook equation (Equation (6)) is commonly used because of the wide range of application and high accuracy, which is derived from the Prandtl hydraulic smooth tube equation and the Nikolas complete rough tube equation [29]. 1 λ k = - 2 lg ( ɛ k 3.7 D k + 2.51 Re k λ k ) k ∈ M k(6) Where ɛ _k is the inner wall equivalent roughness of the kth pipeline, m; Re _k is the reynolds number of the kth pipeline, which can be calculated by Equation (7). Re k = 1.536 Q k p Δ D k μ(7)

Where Δ is the relative density of gas; μ is the dynamic viscosity of natural gas, Pa·s.

(4) Compressibility factor equation

Compressibility factor Z indicates the volume deviation of the actual gas from the ideal gas after compression by the same pressure, which is a correction factor that must be considered when gas flows in a pipe. BWRS equation is widely adopted to the calculation of compressibility factor due to the character of high computational accuracy and wide applicability. BWRS equation is obtained by modifying the BWR equation. This equation has 11 light hydrocarbon component coefficients, each of which can be obtained by the relevant calculation formula and calculation rules. Since pipelines at different locations have different operating pressure, especially for huge gas pipeline systems, there is a large gap between the starting and ending pressure. Therefore, this model will calculate the average compression factor Z k avg of each pipeline separately, which can be calculated by Equation (8). Z k avg = 1 + ( A 1 + A 2 RT + A 3 RT 3 + A 4 RT 4 + A 5 RT 5 ) × ρ k avg + ( A 6 - A 7 RT + A 8 RT 3 ) × ( ρ k avg ) 2 + A 9 RT ( A 7 + A 8 T ) × ( ρ k avg ) 5 + A 10 ( ρ k avg ) 2 RT 3 ( 1 + A 11 ( ρ k avg ) 2 ) exp ( - A 11 ( ρ k avg ) 2 ) k ∈ M k(8)

Where Z k avg is the average compressibility factor of the kth pipeline; A₁-A₁₁ is the BWRS equation coefficients; R is the constants related to the state of natural gas, 8.314 J/(mol·K); T is natural gas flow temperature, K; ρ k avg is the gas average density of the kth pipeline, kmol/m³, which can be calculated by Equation (9). ρ k avg = P k avg Z k avg RT(9)

It’s worth noting that the calculation equations for the friction factor and the compressibility factor are implicit equations, which cannot be directly solved. Hence, Newton–Raphson method is applied in the modeling process.

(5) Compressor performance equation

The compressor is one of the most complex equipment in a natural gas network. The compressor characteristic curve refers to the relationship between head H _j , isentropic compression efficiency η j is , suction flow rate Q j s and rotational speed ω _j . The maps describe the operating performance of the compressor and estimate its operating parameters. The relationship between pressure head H _j , suction flow rate Q j s and speed ω _j is shown in Fig. 3(a). In addition, the relationship between the isentropic compression efficiency η j is and the suction flow rate Q j s and rotational speed ω _j is also represented in Fig. 3(b). The characteristic curves denote the operating performance of the compressor and estimate its operating parameters.

OPEN IN VIEWER

Fig. 3

Compressor characteristic curves: (a) H _j , Q j s , ω _j ; (b) η j is , Q j s , ω _j .

In the process of modeling, a mathematical language is needed to describe the compressor maps. Hence, the characteristic curves are often fitted to a polynomial equation called compressor characteristic equation. In these equations, H j / ω j 2 , Q j s / ω j , η j is are selected as equation variables, and coefficients A _h , B _h , C _h , D _h , A _e , B _e , C _e , D _e are determined by nonlinear least square method. The head characteristic equation and the efficiency characteristic equation are expressed by Equations (10) and (11) separately [30]. H j ω j 2 = A h + B h ( Q j s ω j ) + C h ( Q j s ω j ) 2 + D h ( Q j s ω j ) 3 j ∈ M j(10) η j is = A e + B e ( Q j s ω j ) + C e ( Q j s ω j ) 2 + D e ( Q j s ω j ) 3(11)

Where H _j is the compressor pressure head of the jth CS, kJ/kg; ω _j is the compressor rotational speed of the jth CS, r/min; Q j s is compressor actual suction flow rate of the jth CS, m³/s; η j is is isentropic compression efficiency of the jth CS, %; A _h , B _h , C _h , D _h is head characteristic equation coefficients; A _e , B _e , C _e , D _e is efficiency characteristic equation coefficients.

When a compressor is working, if the suction flow rate Q j s is too low, surge might occur; and if the flow rate is too high, stonewall might occur. Both phenomena can affect the performance of the compressor. Therefore, it is necessary to restrict the flow rate and control the compressor within a certain working range. The surge and stonewall limit curve are on the left and right sides of Fig. 4(a) respectively. The limitation of the suction gas flow Q j s (Equation (14)) can be obtained by fitting the surge limit curve and stonewall limit curve to the polynomial equations Equations (12) and (13). Q j , su s = A su + B su ω j + C su ω j 2 + D su ω j 3(12) Q j , st s = A st + B st ω j + C st ω j 2 + D st ω j 3(13) Q j , su s ⩽ Q j s ⩽ Q j , st s(14)

OPEN IN VIEWER

Fig. 4

Procedures of solving the model with the ɛ-constraint method.

Where Q j , su s is the compressor surge limit flow rate of the jth CS, m/s; Q j , st s is the compressor stonewall limit flow rate of the jth CS, m/s; A _su , B _su , C _su , D _su is the surge limit characteristic equation coefficients. A _st , C _st , D _st is the stonewall limit characteristic equation coefficients.

The discharge pressure P j d of the compressor thus can be calculated after calculating the compressor head, which is defined as Equation (15). H j = Z j s RT k ¯ k ¯ - 1 [ ( P j d P j s ) k ¯ - 1 k ¯ - 1 ](15)

Where Z j s is the gas compressibility factor at compressor inlet side of the jth CS; k ¯ is gas adiabatic index; P j s is the compressor suction pressure of the jth CS, MPa; P j d is the compressor discharge pressure of the jth CS, MPa.

Finally, the compressor power W _j can be calculated by the mass flow M _j and head H _j . The compressor power calculation formula is as Equation (16). W j = M j H j η j is η j mech(16) Where W _j is the compressor power of the jth CS, kW; M _j is the compressor mass flow rate of the jth CS, kg/s; η j is is the isentropic compression efficiency of the jth CS, %; η j mech is the compressor mechanical efficiency of the jth CS, %.

(1) Node flow rate constraint

As a resource, the flow of natural gas into the pipeline over a period of time Q ii in is constrained by the maximum production or processing capacity of the producing area. In addition, for different regions and different types of natural gas users, there are different consumption demands. Therefore, at each GDN, the gas distribution volume Q id out should be greater than the minimum volume. The related constraints can be described as Equations (17) and (18). Q ii , min in ⩽ Q ii in ⩽ Q ii , max in ii ∈ N ii(17) Q id , min out ⩽ Q id out ⩽ Q id , max out id ∈ N id(18)

Where Q ii in is the injection flow rate of the mth GIN, m³/s; Q ii , min in and Q ii , max in are the minimum and maximum allowable injection flow rate of the mth GIN, m³/s; Q id out is the delivery flow rate of the ith GDN, m³/s; Q id , min in and Q id , max in are the minimum and maximum allowable delivery flow rate of the mth GDN, m³/s.

(2) Node pressure constraint

Similarly, gas sources or users have certain requirements on the pressure P _i . The pressure at each node of the pipeline can be expressed as Equation (19). P i , min ⩽ P i ⩽ P i , max i ∈ N(19)

Where P _i is the pressure of the ith node, MPa; P_i,min is the minimum allowable pressure of the ith node, MPa; P_i,max is the maximum allowable pressure of the ith node, MPa.

(3) Pipeline pressure constraint

Natural gas flows in the pipe under pressure, and the inner pressure causes stress to the pipe wall. If the stress reaches the yield strength of the pipe, it would cause permanent deformation of the pipe and eventually cause damage to the pipe. Therefore, to ensure the safe transportation of natural gas, maximum allowable operational pressure (MAOP) parameter needs to be introduced to ensure the safe transportation of natural gas. The pipeline pressure constraint is expressed as Equation (20). P k ⩽ MAOP k ∈ M k(20)

Where P _k is the pressure of the kth pipeline, MPa; MAOP is the pipeline maximum allowable operational pressure.

(4) Compressor rotational speed constraint

Compressor rotational speed is related to equipment performance and is within a certain limitation. Maximum and minimum rotational speed of compressors is expressed as Equation (21). ω j , min ⩽ ω j ⩽ ω j , max j ∈ M j(21)

Where ω _j is the compressor rotational speed of the jth CS, r/min; ω_j,min and ω_j,max are the minimum and maximum allowable compressor rotational speed of the jth CS, r/min.

(5) Compressor power constraint

A compressor needs to meet maximum and minimum power constraints. The constraint is represented as Equation (22). W j , min ⩽ W j ⩽ W j , max j ∈ M j(22)

Where W _j is the compressor power of the jth CS, kW; W_j,min and W_j,max are the minimum and maximum allowable compressor power of the jth CS, kW.

3.1 Multi-objective optimization

Based on the established objective functions and constraints, the decision variables include: delivery flow rate at a specified GDN ( Q id out ), RCs number (δ _j ) and their rotational speed (ω _j ). Among them, delivery flow rate and rotational speed are continuous variables, while RCs number is discrete variable. Hence, the established optimization model is a MINLP model. In an actual project, when the decision variables are determined, the network system optimization operation plan is also determined. Operators can directly use the plan to guide the system’s operation. Based on Equations (1)–(22), the Multi-objective optimization model for gas transmission network can be summarized as follows: { max f 1 ( X ) min f 2 ( X ) subjected to h ( X ) = 0 g ( X ) ⩽ 0(23)

Where, X represents an L-dimensional vector composed of some decision variables x _l ; f₁ (X) and f₂ (X) represent the two objective functions; h (X) means an R-dimensional vector consisting of an equality constraint h _r  (X), and g (X) means an S-dimensional vector consisting of an inequality constraint g _s  (X).

The essence of the multi-objective optimization problem is that in most cases, the sub-objectives may conflict with each other, and the improvement of one sub-objective may cause the performance of another sub-objective to decrease. Therefore, it is difficult to optimize multiple sub-objectives simultaneously. Coordination and compromise can only be performed among them to make each sub-objective function as optimal as possible [31]. The economist Pareto first studied the multi-objective optimization problem in the field of economics and proposed Pareto optimal sets, also known as Pareto front. Each solution in the Pareto front cannot be further optimized for a single target without sacrificing other targets. The main purpose of multi-objective optimization technology is to seek one or more satisfactory solutions in the Pareto front [32].

Scholars have continuously adopted ɛ-constraint method to solve multi-objective optimization problems in the oil and gas field. Alves et al. [33] used the ɛ-constraint method to solve the multi-objective design optimization problem of natural gas pipelines. The goals considered were to minimize unit transportation costs and maximize gas transmission volume. Rui et al. [34] thought minimizing the pump running cost and the number of on/off operations of the liquefied light hydrocarbon pipeline system, and also applied the ɛ-constraint method to solve the model. When the ɛ-constraint method is utilized to solve a multi-objective optimization problem, one of the objective functions are selected as the optimization objective of the model, and the others are transformed into constraints. In this paper, the minimization of the energy cost of the compressor station f₂ (x) is regarded as the optimization objective, and the maximum specified GDN delivery flow rate f₁ (x) is converted into a constraint. Therefore, the structure of the bi-objective optimization model established in this paper is thus transformed into a set of single-objective optimization models as showed in Equation (24). { min f 2 ( X ) subjected to f 1 ( X ) ⩾ ɛ u u = 1 , 2 , . . . , U h ( X ) = 0 g ( X ) ⩽ 0(24)

Each ɛ _u represents the lower bound at a specific GDN delivery flow rate in a single-objective model. Before the solution, the variation range of ɛ _u can be adjusted repeatedly by trial calculation. Each time a single-objective model is solved, a Pareto optimal solution of the original bi-objective model thus can be obtained. Finally, by solving a certain number of single-objective models, the Pareto front of the original bi-objective model can be obtained. The steps of solving the model with the ɛ-constraint method are summarized in detail, as shown in Fig. 4.

After the multi-objective optimization model is processed into a set of single-objective optimization models by the ɛ-constraint method, each single-objective optimization model is then programmed in GAMS (General Algebraic Modeling System) version 24.8.2. Figure 5 shows the basic part of modeling along with coding according to GAMS software, which inclouds some important objective functions, equality constraints and inequality constraints. This software provides a concise language system to describe the complex mathematical model and encapsulates the algorithm in the system. After the user describes the mathematical model, the system will automatically call the algorithm to solve it. The solutions are obtained by using the DICOPT solver. Solver decomposes the MINLP model into NLP and MIP subproblems. Two subproblems are solved by generalized reduced gradient method and branch and bound method respectively. The computer utilized has an Intel (R) Pentium (R) CPU G4560 3.50 GHz processor and 8 GB of random access memory (RAM). For the case described in Section 4, a total of 115 single-objective optimization models were solved, taking 34.385 s.

OPEN IN VIEWER

Fig. 5

Basic part of modeling along with coding according to GAMS.

5.1 Optimization results

In the above study, we established a multi-objective optimization model involving two optimization goals (maximizing the specific GDN delivery flow rate and minimizing the energy cost of CSs). The ɛ-constraint method and GAMS/ DICOPT solver were adopted to solve the model. A large practical tree-topology gas transmission network is investigated as an optimization example. During the solution process, the initial value range of ɛ _u is 428.57–950×10⁴m³/d, and the step size is 4.32×10⁴m³/d. However, once ɛ _u exceeds 921.05×10⁴m³/d, no feasible solution can be obtained. This is because the pressure drop increases as the flow increases after delivery flow rate at GDN19 exceeds 921.05×10⁴m³/d. Even if the last CS provides the maximum outbound pressure of 9.85 MPa, it could not meet the requirement that the pressure at GDN15 is not less than 4.8 MPa. Therefore, the value range of ɛ _u is adjusted to 428.57–921.05×10⁴m³/d, and the ɛ-constraint method finally solves 115 single-objective problems.

The Pareto front calculation results of the two-objective optimization problem proposed in this paper are calculated and demonstrated (Fig. 8). Each dot on the Pareto front represents a Pareto optimal solution. It can be seen from the figure that the power consumption of the CS increases with the increase of the delivery flow rate at GDN 19. This tells that the two optimization goals are in conflict with each other and cannot be achieved the optimum value at the same time. As the flow rate increases, the growth rate of compressor energy consumption also gradually accelerates. Points A, B, and C in Fig. 8 represent the Pareto optimal solutions under the conditions of minimum, medium, and maximum GDN19 delivery flow rate respectively. Then this article will use these three points as examples to discuss the differences between different Pareto optimal solutions. Besides, it can be seen from Fig. 8 that the function value of the optimal solution set (Pareto front) of the optimization problem is located at the lower boundary line of the entire objective function space. In the entire feasible region, all feasible solutions except the Pareto front are dominated by the solutions on the Pareto front. Domination means that in the two solutions compared with each other, all the objective functions of one solution are not worse than the other solution, and at least one objective function is better than the other solution. Take the points X1, X2, and X3 in Fig. 2 as an example. Point X1 is located in the feasible region but is not optimal. This point is dominated by points X2 and X3 on the Pareto front. That is, when the flow is the same, the energy consumption of the compressor station at X2 is lower than X1. When the energy consumption is the same, the flow rate delivered by point X3 is higher than that of X1. Therefore, the solution on the Pareto front has the least conflict of goals compared with other solutions, and can provide a better choice for decision-makers.

OPEN IN VIEWER

Fig. 8

Pareto front for bi-objective optimization of the gas transmission network.

5.1.1 Comparison between Pareto optimal solution and actual operation scheme

First, point A in the Pareto front is selected to compare with the actual operation scheme in order to analyze the feasibility and optimization effect of the optimization result. Point A is the point where the power cost of the CS is the smallest in the Pareto front (Fig. 8), while the GDN 19 delivery flow rate at this point is exactly the same as the actual operation scheme. Therefore, the case of point A has a pipeline flow situation that is completely consistent with the actual operation scheme, which can be used to represent the effect of the optimization results. The power consumption cost of CSs in the actual operation scheme is 70.94×10⁴ ¥, while the optimization result is 63.99×10⁴ ¥. By comparing the optimization result of point A with the actual power consumption cost of CSs, the optimization result saves a power consumption cost of 9.79%.

We compared the number of RCs between the optimization results and the actual operation scheme (Fig. 9(a)). It can be seen from the figure that, except for the CS 1, the RCs number in each CS is the same in both schemes. The RCs number in CS 1 is one in the actual operation scheme while two in the optimization result. One more compressor runs in the optimization result than the actual operation scheme. Besides, the comparison between the optimization results and the actual CS rotating speed in each station is shown in Fig. 9(b). It can be seen from the figure that only in CS 3, the rotating speed in the optimization results is higher than the actual operation scheme. As for other CSs, the rotating speed in the optimization results is lower than the actual operation scheme. Besides, the comparison of each CS total power is illustrated, as shown in Fig. 9(c). It is observed that the power of CS 1 and CS 3 in the optimization result is higher than the actual operation scheme. But the power of CS 2, CS 4, CS 5, CS 6 and CS 7 in the optimization result is lower than the actual operation scheme, especially the CS 5. Obviously, optimization strategy of the optimization result is to increase the power of the upstream compressor station and reduce the power of the downstream compressor station. This is also the reason why the RCs number in CS 1 of optimization result is increased. Next, we will further analyze this optimization strategy by pipeline pressure drop.

OPEN IN VIEWER

Fig. 9

Comparison between the optimized results and the actual operation scheme: (a) running compressor number; (b) compressor rotating speed; (c) CS total power.

The comparison of the trunk pipeline pressure drop between the optimization result and the actual operation scheme is shown in Fig. 10. The following conclusions can be drawn through analyzing the figure. In the optimization result, the discharge pressure of CS 1, CS 2, CS 4 and CS 5 are higher than the actual operation scheme, the discharge pressure of CS 6 is almost equal to the actual operation scheme, but the CS 7 discharge pressure is lower than the actual one. It is worth noting that, as the discharge pressure of the upstream 4 CSs in the optimization results increased, the suction pressure of the adjacent CSs also increased, thus reducing the boost degree of this adjacent compressor stations to some extent. This phenomenon is more obvious in CS 5 and CS 6. Combined with Fig. 10(c) above, it can be seen that the reduction of the boost degree of this two CSs save a lot of power consumption for the optimization results. This is also the main reason why the optimization result is better than the actual operation scheme. In addition, the comparison between the optimization results and the actual terminal pressure of each pipeline is shown in Fig. 11.

OPEN IN VIEWER

Fig. 10

Comparison of trunk pipeline pressure drop.

OPEN IN VIEWER

Fig. 11

Comparison of pipeline terminal pressure.

It can be seen from the figure that at the branch 3, optimized end-point pressure of the pipeline is just the minimum allowable pressure of 4.8 MPa, which just meets the minimum pressure constraint of the gas point. However, the end pressure of branch 3 of the actual operation scheme is 5.28 MPa and obviously there is a waste of energy. This is another reason that optimization results are more economical.

5.1.2 Comparison between different pareto optimal solutions

The next stage is to analyze the difference between the optimization results of different Pareto optimal solutions. Three representative points are selected for analysis from the Pareto front, including points A, B, and C (Fig. 8). The objective function values at points A-C are shown in Table 6. As we can be seen from the table: (1) Case of point A has the lowest delivery flow rate at GDN 19 and the power cost of the CS; (2) Case of point C has the highest delivery flow rate at GDN 19 and the power cost of the CS; (3) The two objective function values of point B are between points A and C. Besides, Fig. 12(a) shows the number of RCs at points A to C. The difference mainly occurs at CS 3 and CS 6, of which the number of RCs in case of point A are two, 3 and 2 respectively of CS 6 in case of B, and all 3 in case of C. The compressor rotating speed of each station at points A-C is shown in Fig. 12(b). The differences mainly appear in CS 5, CS 6 and CS 7.

OPEN IN VIEWER

Table 6

Optimization results of objective functions for Points A-C

Objective functions	Optimization results
	Point A	Point B	Point C
Delivery flow rate at GDN 19 (104m3/d)	428.57	679.13	921.05
CS power consumption cost (104¥)	63.99	82.10	109.45

OPEN IN VIEWER

Fig. 12

Points A-C: (a) number of RCs; (b) compressor rotating speed.

5.1.3 Analysis of decision variables variation with pareto front

We analyzed the changes in the number of RCs at each station and in the speed with the Pareto front to further find out the detailed change trend and key change points of the decision variables of the Pareto front. The results are shown in Figs. 13 and 14. Specifically, the dashed lines in the figure represent the upper or lower limits of the variable, respectively. Figure 13(a) and (b) show the changes of δ₃ and δ₆ in the Pareto front, respectively. Since the number of RCs of other CSs has always remained the same, they are not discussed. It can be seen from Fig. 13(a), δ₃ is 2 at the point where delivery flow rate of GDN 19 is 653.21×10⁴m³/d. After reaching this point, δ₃ becomes 3. As can be seen from Fig. 13(b), it is initially used for 2 compressors, and maintains a large flow rate range, and it becomes 3 until the distribution flow rate reaches 813.05×10⁴m³/d.

OPEN IN VIEWER

Fig. 13

Change of running compressors number with the Pareto optimal solution: (a) CS 3; (b) CS 6.

The change of the compressor rotating speed at the Pareto front of each of the seven CSs is shown in Fig. 14. As shown in Fig. 14(a), it grows slowly first, then undergoes a process of rapid decrease and rise, and finally maintains slow growth. It can be seen from Fig. 14(b) that the lower bound is 5900 r/min at all points. So, the changes in flow have little effect on ω₂. As shown in Fig. 14(c), ω₃ is maintained at a high level, and gradually increases first, and then a small decrease occurs at GDN19 at a gas flow rate of 627.29×10⁴m³/d. This is because the power of a single compressor of CS 3 has reached the upper limit of 6600 kW at this point. ω₃ must be reduced in order to meet the power constraint. Then, at the point where the distribution flow is 653.21×10⁴m³/d, ω₃ shows a significant decrease. The reason is that at this point, the number of compressor start-ups of CS 3 became three, and the inlet flow rate of a single compressor decreases, and the speed also decreases accordingly. ω₃ (Fig. 14(d)) keeps increasing after that. It is noted that ω₄ is maintained at the lower limit of the speed at first, followed by a process of rapid rise and then decline occurs, and finally maintains at the lower limit. ω₅ (Fig. 14(e) increases sharply at the beginning, but suddenly stops growing at a point of 679.13×10⁴m³/d, and then keep fluctuating within a small range. This is because that at a point of 679.13×10⁴m³/d, the single compressor power of CS 5 has reached the upper limit of 9800 kW, and no more compressors can be put into use at this station. As shown in Fig. 14(f), ω₆ initially maintains the lower limit, and it does not start to increase significantly until the distribution flow is 679.13×10⁴m³/d. However, at the point where the distribution flow is 813.05×10⁴m³/d, ω₆ suddenly drops, and then continues to increase. This is because the number of start-ups in CS 6 has changed to three at this point, the inlet flow rate of a single compressor has been greatly reduced, and the speed has been greatly reduced. The change of ω₇ is shown in Fig. 14(g). It can be seen from the figure that ω₇ maintains the lower limit of the rotating speed from the beginning. It does not start to rise until the flow is 856.25×10⁴m³/d.

OPEN IN VIEWER

Fig. 14

Change of compressor rotating speed with the Pareto optimal solution: (a) CS 1; (b) CS 2; (c) CS 3; (d) CS 4; (e) CS 5; (f) CS 6; (g) CS 7.

Through the above analysis, it can be found that the key points of each decision variable’s change with the Pareto front point are the same in some CSs. This is because at these key points the number of start-ups has changed, or the power of a single unit has reached the upper limit, resulting in varying of decision variables for CSs along the pipe. Therefore, the various decision variables affect each other and ultimately affect the optimization result in the optimization process. In addition, the solution of this optimization problem includes 15 decision variables, and the decision space where the solution vector is located is a 15-dimensional vector space. According to the previous analysis of the change of decision variables in the Pareto optimal solution, we can see that among all decision variables, the values of 5 decision variables including δ₁, δ₂, δ₅, δ₇, ω₂ are always boundary values. The variable value of the running compressors number δ₁, δ₂, δ₅, δ₇ is always the upper boundary value of 2, and the compressor speed variable ω₂ is always the lower boundary value of 5900 r/min. Therefore, the solution vectors in the Pareto optimal solution set are all on the boundary of the original decision space, and the resulting solution space is the boundary subspace of the original decision space.

In general, every point on the Pareto front is a feasible solution. The final optimal solution can be selected from the Parato front via a process of decision-making. In the natural gas industry, decision makers are one or a group people who have a better understanding of the network operation and the relevant costs. During the decision-making process, operators can use practical information and experience to select the best operational scheme, to satisfy contract gas delivery requirement with the lowest economic cost.

The optimization results are the results of the joint action of multiple operating parameters. Consequently, we performed a sensitivity analysis in order to know the impact of each operating parameter changes on the optimization results. Sensitivity analysis is carried out for pipeline parameters and compressor parameters, respectively. As a result, we can better predict and understand the optimization results.

The section mainly analyzes the effect of inlet flow rate and inlet pressure changes on the pressure drop of the pipeline. A section of natural gas pipeline with an inner diameter of 1016 mm and length of 125.6 km is selected for analysis. The range of values for each parameter for sensitivity analysis is shown in Table 7. The analysis results are obtained through simulation of the previously established pipeline and compressor station models. It can be seen from the figure (Fig. 15(a)) that as the flow increases, the pressure drop also continues to increase, and there is a non-linear relationship between the two. The influence of the inlet pressure changes on the pressure drop is shown in Fig. 15(b). It can be seen from the figure that as the pressure of the pipeline inlet increases, the pressure drop of the pipeline continuously decreases. Therefore, maintaining the pressure of the entire pipeline at a high level can reduce the pressure drop and save system energy consumption.

OPEN IN VIEWER

Fig. 15

Effects of various parameters of the pipeline on the pressure drop: (a) flow rate; (b) inlet pressure.

Coming to the next is the sensitivity analysis of compressor station parameters. This section mainly analyzes the impact of changes in compressor intake gas flow rate, speed and suction pressure on discharge pressure and energy consumption costs. The influence of each parameter, considering the upper and lower bounds of the parameters for different compressor startups, is represented in Fig. 16. A compressor station in the calculation example is selected for analysis in the process analysis. The ranges of the three parameters are shown in Table 8.

OPEN IN VIEWER

Fig. 16

Impact of parameters on compressor performance under different running compressor numbers: (a) inlet flow rate vs. discharge pressure; (b) inlet flow rate vs. energy costs; (c) suction pressure vs. discharge pressure; (d) suction pressure vs. energy costs; (e) rotating speed vs. discharge pressure; (f) rotating speed vs. energy costs.

Figure 16(a) and (b) demonstrate the influence of inlet gas flow rate on discharge pressure and energy cost. As shown in Fig. 16(a), as the inlet gas flow rate increases, the discharge pressure decreases continuously. Besides, the increase of the inlet gas flow rates also leads to an increase in the energy cost of the compressor (Fig. 16(b)), which means that more compressors are needed to meet the increase of inlet flow rate. At the same inlet flow rate, increasing the number of RCs would increase the discharge pressure, causing more energy costs. The relationship between compressor suction pressure, discharge pressure and energy cost is displayed by Fig. 16(c) and (d). As shown in the figure, an increase in suction pressure leads to an increase in discharge pressure and energy costs. Further, at the same suction pressure, an increase in the number of RCs has resulted in an increase in discharge pressure and energy costs. But the increase in energy costs is higher than the discharge pressure. The relationship between rotating speed, discharge pressure and energy costs are shown in Fig. 16(e) and (f). As shown in Fig. 16(e), increasing the compressor’s rotating speed can cause higher discharge pressure. But higher rotational speeds also mean higher energy costs (Fig. 16(f)).

In summary, higher compressor discharge pressure can be obtained by increasing the compressor’s suction pressure, rotating speed, number of RCs or reducing the inlet flow. However, compressor suction pressure and inlet gas flow rate are mainly affected by the inlet natural gas conditions and the operation of upstream units. From the perspective of compressor station variables, it is possible to consider increasing the number of RCs or the rotating speed. In addition, an increase in the discharge pressure of a compressor station would also lead to an increase in the suction pressure of the downstream compressor station, thus increasing the discharge pressure of the downstream compressor station accordingly. But increasing the number of RCs or the rotating speed inevitably leads to an increase in the energy cost of the compressor station. Therefore, an optimization method must be applied to balance the influence between various parameters to achieve the optimal system operating conditions.