Multi-Objective Optimization of Gas Pipeline Networks

: The main goal of this paper is to prove that bi-objective optimization of high-pressure gas networks ensures grater system e ﬃ ciency than scalar optimization. The proposed algorithm searches for a trade-o ﬀ between minimization of the running costs of compressors and maximization of gas networks capacity (security of gas supply to customers). The bi-criteria algorithm was developed using a gradient projection method to solve the nonlinear constrained optimization problem, and a hierarchical vector optimization method. To prove the correctness of the algorithm, three existing networks have been solved. A comparison between the scalar optimization and bi-criteria optimization results conﬁrmed the advantages of the bi-criteria optimization approach.


Introduction
The introduction of a Third Party Access (TPA) regime has been a central element of liberalization of the gas industry.The objective is to foster competition in the gas market, improve supply efficiency and bolster infrastructural investments, thereby strengthening energy security.However, TPA provisions have a game-changing impact on the gas supply market by introducing competition, breaking the incumbent's monopoly, lowering prices and adding complexities to service provision by the system operator.The main responsibility of the gas network operator is to maintain safe and efficient operation of the gas supply system in real time and to ensure security of supplies in the short, medium and long term.Operation optimization of natural gas pipelines has received increasing attentions due to such advantages as maximizing the operating economic benefit and the gas delivery amount.
This paper is concerned with the high-pressure gas networks under steady-state conditions.If we assume that the goal is to minimize the operational costs, then the steady-state optimization algorithm determines the pressure values in the network nodes, flow values in the pipes, flows in the compressor stations and the compression ratio of compressor stations, minimizing the objective function.
Different methods are proposed in the literature to solve the problem of this nature.In particular, non-linear models are applied because non-linearities can be represented in their original form.In [1] mixed integer non-linear programming is used to tackle the steady-state case; In [2] sequential linear programming is employed to study the behaviour of an integrated gas and power system.In [3] mixed integer linear programming is applied for the steady-state optimization of gas flows.Nowadays, it is possible to state that mixed integer linear programming techniques are mature because they are fast, robust, and are able to solve problems with up to hundreds of thousands of variables [4].
Several locally mixed integer linear programming formulations used to piecewise linearize a non-linear function have been proposed in the literature.A unified framework can be found in [5].
It should be noted that the gas network control is a complex process which can be operated on the basis of various criteria.The criteria are formulated by the network operator and depend primarily on the operating conditions of the network, structure of the network, its hydraulic properties and type of consumers.Under certain conditions, scalar optimization is sufficient to satisfy the operator requirement while in other cases one should look for compromise solutions, which are optimization problems with an objective function in the form of vector.
Solving the single objective optimization problems is far more common than solving the multi-objective problems, since there appears to be no generally effective and efficient method available for solving the multi-objective problems directly as they are.Typically, a multi-objective problem is to be effectively converted to a single objective problem before applying an optimization algorithm.This conversion can be done easily by first deciding the relative importance for each objective a priori.Then, for example, the decision-maker may combine the individual objective functions into a scalar cost function (linear or non-linear combination), which effectively converts a multi-objective problem into a single objective one.
Nonetheless, there are more and more papers in which the task of gas transmission system optimization is formulated as a multi-criteria task.
In [6] a bi-criteria optimization algorithm is formulated (running cost and maximum capacity) to solve a simple structure of the network, assuming a steady-state flow.
In [7] an algorithm for the optimization of design and operation of a gas transmission system is employed.The objective function consists of two components: one describes the costs of construction and operation of pipelines while the other specifies the costs of construction and operation of the compressor station.The problem was solved using mixed non-linear programming In [8] an algorithm for investment and operational costs minimization of a high-pressure gas network is applied.The mixed integer non-linear programming was used to find optimal suction and discharge pressures of a fixed number of compressor stations, as well as the length and the diameter of the pipeline segments.The optimization goal was to minimize the annual costs, assuming 10-year investment depreciation.The main disadvantage of this algorithm is the assumption that pipeline diameters are treated as continuous variables.
In the paper by [9], a multi-objective ant colony optimization technique for pipeline optimization has been developed to solve the multi-objective gas pipeline transportation problem.The multi-objective problem considered is about minimizing fuel consumption in the compressors and maximizing the throughput.
In [10] a multi-objective optimization method to trade-off reliability and power demand in the decision-making process is developed.In the optimization, the steady-state behaviour of the natural gas pipeline networks is considered, but the uncertainties of the supply conditions and the customer consumptions are accounted for.The multi-objective optimization is about finding the operational strategies which minimize power demand and the risk of gas supply shortage.
In [11] three objective functions to simultaneously optimize: minimizing the fuel consumption in the compressor stations, maximizing the network throughput and maximizing the percentage of added hydrogen at the network entrance is employed.The NSGA-IIb genetic algorithm coupled with a Newton-Raphson procedure of the MATLAB toolbox (MathWorks, Inc., Natick, MA, USA) is implemented.
In [12] the short-and medium-term planning problems of the regimes of multi-line technical gas pipeline corridors (MLGP) of the Russian gas supply system is considered.The fall in gas production caused by depleted gas fields leads to a decrease in the load of some operating MLGPs.
The purpose of this article is to develop and test mathematical models and a computer program in order to support the adoption of dispatch solutions for managing modes of large MLGPs under conditions of incomplete loading.
The paper by [13] addresses the Line Pack Management of the "GZ1 Hassi R'mell-Arzew" gas pipeline.For a gas pipeline system, the decision-making on the gas line pack management scenarios Energies 2020, 13, 5141 3 of 19 usually involves a delicate balance between minimizing fuel consumption in the compression stations and maximizing the gas line pack.In order to select an acceptable Line Pack Management of Gas Pipeline scenario from these two angles for "GZ1 Hassi R'mell-Arzew" gas pipeline, the idea of multi-objective decision-making has been introduced.
Pipeline failures due to natural disasters, corrosion and manufacturing defects have generated a tremendous economic loss.To recover efficiently from the disruption with a minimum loss, [14] propose a multi-objective optimization model that minimizes both the loss from disruption and recovery time.Recovery on disrupted pipelines may incur an additional recovery cost but, on the other hand, it can alleviate the cost of unmet demand.Considering the gas production and transmission costs, we use multi-objective optimization to study the recovery decisions in a pipeline network.
The paper by [15] proposes the multi-objective optimization of the design of natural gas transmission networks to support the decision of regulatory authorities.Problem formulation involves two objective functions: minimization of the transportation fare and maximization of the transported gas volume.These design parameters of the pipeline project must be previously established by the regulatory agency, considering an attractive return on investment for the entrepreneurs and the demands of current and future consumers.
The first part of this function, which expresses the running cost of the compressors, should be minimized (minimization of the operating costs), and the second part, which expresses the line pack, should be maximized (ensuring delivery security).In order to establish a relationship between the two functions, weighting factors were used.In this case study, it was assumed that both functions are equally important for the user, therefore the weighting coefficients are equal to 0.5.Weights must be constant and must be functions of the original objectives.The values of the weighting factors are arbitrarily chosen by the user (transmission system operator).Depending on the hydraulic properties of the system, the capacity, the structure, the number of compressor stations and on the amount of the transported gas, the operator decides what is more important: the running cost of the compressor stations or security of gas supply.
The elements of vector x are discharge pressure values of the working compressor stations.Determining components of the x vector, which will minimize the formulated objective function is the solution to the bi-criteria optimization problem.This involves finding such discharge pressure values for each compressor station which will ensure that the operational costs of the transmission system will be a compromise between maximizing line pack costs and minimizing running costs.
Minimizing the running costs of the transmission system is equivalent to minimizing power consumption in the compressor stations.The total power consumption in the system is expressed by Equation (2): where: Energies 2020, 13, 5141 The second part of the objective function relates to gas line pack in the network.Total gas volume in the network as a function of geometrical dimensions of the pipelines, and gas parameters can be expressed by the following Equation ( 3): where: V In order to compare the power consumption as a result of running cost minimization with power consumption as a result of line pack maximization, the second term of objective function was expressed by Equation ( 5): where: N s -power loss corresponding to given flow parameters [W], Q v,i -volumetric gas flow through the i-th pipeline at real conditions [m 3 /s], ∆p i -pressure drop for the i-th pipeline [Pa].
The results of calculations (see Table 1) have shown that the greater line pack in the network, the higher operational costs of the compressor station.Calculations were made on the pipeline with an internal diameter of 0.7 m and a length of 50 km.This means that the first part of the objective function is to be minimized and the second part is to be maximized.We are searching for an appropriate discharge pressure values for working compressor stations which will ensure a compromise between these two components.

Constraints
For high pressure gas networks, operating above 7.0 Bar-gauge, the Panhandle equation is used.For each pipe k [having node l (k) at its left, and r (k) at its right] the pressure drop equation has the form: where: p-pressure; Q-flow through pipe. where: L-pipe length; D-pipe diameter; E-pipe efficiency stated as a constant (~0.9).
For the whole of the network, the equality constraints are the following (more details in [16,17]): where: A-the nodal-branch incidence matrix (dim A = n × m); n-the number of nodes; m-the number of branches; P-the vector of squared nodal pressures (dim P = n × 1); K-the unit -nodal incidence matrix (dim K = n × r); r-the number of units; Energies 2020, 13, 5141 6 of 19 F-the vector of flows through units (dim F = r × 1); L-the vector of loads; ∆P-the vector of squared drop pressures (dim ∆P = m × 1); K f -the vector of pipe constants (dim K f = m × m); Q-the vector of flows through pipes (dim Q = m × 1).

Operational Constraints
The gas transmission system contains three types of controllable units, referred to as "non-pipe elements".The upper and lower bounds are placed on all non-pipe element flows, as well as inlet and outlet pressures.These apply only when the non-pipe element is switched on.
In addition, bounds may be placed on any pressure in the network (including outlets and inlets).Sources and regulators are modelled using linear constraints, which are bounds on the pressure and flow.Compressors have both linear and non-linear constraints.

Compressors
The centrifugal compressors used on the transmission system have specific characteristics.The operating regime can be expressed by what is known as an "envelope" (Figure 1).If the increase in pressure produced by a machine is plotted graphically against the flow, four limits which enclose an area in which the compressor can properly run emerge.The limits are defined as: • "surge": this is the point at which the flow through the compressor becomes so low that a reversal of flow can occur, which can be damaging to the compressor by causing high stress in the bearings or in the impeller; • "choke": at the opposite end of the diagram, a compressor can reach choke.When the pressure ratio is low, there comes a point at which no further increase in the flow through the compressor is possible; • "maximum and minimum speed": obviously a compressor can run up to a given maximum speed consistent with machine safety, and equally there is a minimum speed line.
Energies 2020, 13, x FOR PEER REVIEW 6 of 18 In addition, bounds may be placed on any pressure in the network (including outlets and inlets).Sources and regulators are modelled using linear constraints, which are bounds on the pressure and flow.Compressors have both linear and non-linear constraints.

Compressors
The centrifugal compressors used on the transmission system have specific characteristics.The operating regime can be expressed by what is known as an "envelope" (Figure 1).If the increase in pressure produced by a machine is plotted graphically against the flow, four limits which enclose an area in which the compressor can properly run emerge.The limits are defined as: • "surge": this is the point at which the flow through the compressor becomes so low that a reversal of flow can occur, which can be damaging to the compressor by causing high stress in the bearings or in the impeller; • "choke": at the opposite end of the diagram, a compressor can reach choke.When the pressure ratio is low, there comes a point at which no further increase in the flow through the compressor is possible; • "maximum and minimum speed": obviously a compressor can run up to a given maximum speed consistent with machine safety, and equally there is a minimum speed line.The surge line is formulated by the inequality: where: Q-the flow through compressor (m 3 /h) and a1, b1 are specified coefficients.
The choking line is formulated by the inequality: The surge line is formulated by the inequality: Energies 2020, 13,5141 where: Q-the flow through compressor (m 3 /h) and a 1 , b 1 are specified coefficients.
The choking line is formulated by the inequality: By construction, compressors have a maximum and a minimum speed limit.Operating beyond this upper limit may result in damaging the compressor.On the other hand, an unacceptable efficiency occurs while operating below the lower limit.These constraints are formulated by the inequality: where: RPMAX-maximum speed (R.P.M.); RPMIN-minimum speed (R.P.M.); a, b, c, d, e-coefficients.
For the purpose of the developed algorithm, non-linear constraints of the envelope have been linearized [18,19].Linearization is obtained in the form of a convex quadrangle describing the operative envelope from the outside.The latter is transformed to a set of linear inequalities in terms of pressure and volumetric flow.
When a converged optimum solution is obtained, the real envelopes are inspected to ensure that no violations are present.
Since all the compressors do not usually have the same characteristics, problems may occur when running different types of compressors at the same time.For instance, in a series arrangement the maximum flow limit through one of the compressors may be lower than for the rest of the compressors, therefore, a flow limit is imposed on those compressors with a higher flow constraint.This is stated by the inequality: where: Q max -maximum compressor flow; Q-compressor inlet flow.
For the parallel arrangement, all the pressure ratios have to be equal.Therefore, the compressor with the lowest pressure ratio imposes a limit on the rest of the compressors with a higher pressure ratio.This is stated by the inequality, CR max ≥ p d p s (15) where: CR max -maximum compressor ratio.
Finally, the compressors have a maximum power, which should not be exceeded.Developing more power than indicated by the limit may destroy the compressor.The power required by the prime mover is stated by Equation (1).

Regulators
Flow through regulators is always unidirectional, which is stated by: p o > p i (17) where: p o -outlet nodal pressure; p i -inlet nodal pressure.
There is an equation which relates inlet-outlet pressure difference with the flow.If, for instance, the maximum drop pressure characteristic is known for a particular regulator, the maximum permissible flow can be obtained through this equation, and vice versa.In our problem, we do not use that equation; instead, we use the maximum flow obtained from that equation.The flow inequality constraint is where: Q max -maximum flow through the regulator; Q-actual flow through the regulator.
Finally, the outlet pressure remains constant for any flow rate.This is stated by:

Valves
There are several types of valves, for our purposes only the isolating valves are of interest.These are used to interrupt the flow and to shut off sections of the network.In our problem valves are represented by two inequalities.The first one states that pressure difference is always greater upstream, i.e., p o > p i (20) where: p o -outlet pressure; p i -inlet pressure.
There is an equation for valves which relates their flow within the pressure difference throughout the valve.Since the maximum flow through the valve is given by the constructors, we can obtain the maximum pressure difference from the valve equation.Thus, we only need to define one constraint.In our problem, we have chosen the flow inequality: where: Q max -maximum flow through the valve; Q-actual flow through the valve.

Compressor Stations
The inequality constraints imposed on each compressor station are as follows: Energies 2020, 13, 5141 N ≤ N max (24) Q f ≤ Q max (25) where:

Nodes and Pipelines
Constraints for nodes and pipelines are as follows: w ≤ w max (27) where:

Bi-Criteria Algorithm
To prove the advantages of the bi-criteria method over the scalar method to control the gas transmission system, a comparative analysis using both methods was performed for the same gas structures [20].An arithmetic average of the optimization results of running cost minimization and line pack maximization calculated using scalar optimization were compared with vector optimization results.The Rosen method was chosen for scalar optimization The Rosen method was also used, together with hierarchical method to develop the bi-criteria algorithm.The choice of methods was based on the literature analysis and our own research.

Rosen Method
The Rosen's gradient projection method [21][22][23][24] is an alternative to the Zoutendijk's method [25].Both methods are based on the concept of feasible directions.The Zountendijk's method, however, requires a solution of an auxiliary linear optimization problem to find a usable direction, while the Rosen's method uses the projection of the negative of the objective function gradient onto the active constraints.The Rosen's method requires linear constraints.In other words, it is assumed that non-linear constraints have first been linearized at some point in the design space.However, in general, the objective function remains non-linear.The Rosen method is the most common in the group of feasible directions methods due to its effectiveness.The key factor of this method is the direction gradient projection on the surface tangent to the constraints.At the current point, the constraints given by the equation g j (x) ≤ 0, j = 1, 2, . . ., m (28) can be split into two groups: active constraints g j (x) = 0, j = j 1 , j 2 , . . ., j p ; -non-active constraints g j (x) < 0, j j 1 , j 2 , . . ., j p .
The gradients of the p active constraints are given by: ∇g j (x) = a 1j , a 2j , . . ., a nj T , j = j 1 , j 2 , . . ., j p ( 29) Matrix N(n × p) is given by: N = ∇g j1 , ∇g j2 , . . ., ∇g jp (30) Matrix: is called the projection matrix.It projects the vector-∇f(x) onto the intersection of all the hyperplanes perpendicular to the vectors ∇g(x).
It is assumed that the constraints which are active are independent so that the columns of the matrix N will be linearly independent and hence (N T N) will be non-singular.
The new iteration point is given by where: α-is a step length along the search direction S k .
Solution search is along such a direction S which is periodically updated.This method eliminates the constraints which do not contribute to improving the search direction from the projection matrix.

Hierarchical Method
A survey of current continuous non-linear multi-objective optimization concepts and methods is presented in [10,[26][27][28][29][30][31].The authors conclude that the selection of a specific method depends on the type of information provided in the problem, the user's preferences, the solution requirements, and the availability of software.A hierarchical method has been chosen because in this method the objectives are ranked by the user in order of importance.The optimum solution x * is then found by minimizing the objective functions starting with the most important ones and then proceeding according to the order of importance of the objectives.Let the subscripts of the objectives indicate not only the objective function number, but also the priorities of the objectives.
Thus, f 1 (x) and f k (x) denote the most and least important objective functions, respectively.The first problem is formulated as: subject to: and its solution x 1 and f 1 (x 1 ) is obtained.
Then the second problem is formulated as: subject to: This procedure is repeated until all the k objectives have been considered.For the i-th problem we have: subject to: Finally, the x k solution is taken as the desired x * solution of the multi-objective optimization problem.

Results
Optimization was performed for three different gas network structures which are presented in Figures 2-4.
Energies 2020, 13, x FOR PEER REVIEW 11 of 18 Finally, the xk solution is taken as the desired x * solution of the multi-objective optimization problem.

Results
Optimization was performed for three different gas network structures which are presented in Figures 2-4.Pipeline geometry of network 1 is presented in Table 2.    Pipeline geometry of network 3 is presented in Table 8.Energies 2020, 13, 5141 13 of 19

Network 1
Pipeline geometry of network 1 is presented in Table 2. Nodal input data of network 1 is presented in Table 3. Optimization results for network 1 are presented in Table 4.  Nodal input data of network 2 is presented in Table 6.Optimization results are presented in Table 7.

Network 3
Pipeline geometry of network 3 is presented in Table 8.Nodal input data of network 3 is presented in Table 9. Optimization results for network 3 are presented in Table 10.

Results Analysis
Table 11 presents the results of the tested networks.In each scenario, the bi-criteria optimization objective function value was lower than the value of arithmetic average of the objective function of minimizing power and maximizing gas capacity.The difference increases as the dimension of the network increases.For networks 1, 2 and 3, the differences are 10%, 17% and 28%, respectively.Comparison of optimization results is presented in Figure 5.
The costs of compressor stations electric energy consumption were calculated for each scenario assuming that 1 kWh of electricity costs 0.55 Polish zloty (PLN).They are shown in Tables 12 and 13.The difference increases as the dimension of the network increases.For networks 1, 2 and 3, the differences are 10%, 17% and 28%, respectively.Comparison of optimization results is presented in Figure 5.The costs of compressor stations electric energy consumption were calculated for each scenario assuming that 1 kWh of electricity costs 0.55 Polish zloty (PLN).They are shown in Tables 12 and 13.Cost minimization should not be the only criterion to manage the transmission system.It should be noted that the operator is responsible for gas delivery to the customer under specific terms.It may occur that in certain situations (forecast error, sudden change in ambient temperature, inconsistent gas consumption with the contract by one of the recipients), the operator will not be able to ensure the gas supply under specific parameters resulting from the concluded contracts.Security of supply in all conditions can be achieved by increasing the gas line pack in the network, which in turn results in increased system operational costs.This means that gas transmission system optimization should take into account two factors: minimization of fuel consumption in terms of operating costs and maximization of gas line pack in the system in terms of ensuring delivery security.Therefore, the algorithm should find a compromise solution between these two factors.
The presented tables clearly show that using bi-criteria optimization to manage the gas transmission network is cheaper than in the case of scalar optimization treated as an arithmetic average of the objective function of power minimization and line pack maximization.

Conclusions
The results of the study confirm that bi-criteria optimization allows for much cheaper gas transmission system management.It should be noted that the steady gas flow state in the case

4 of 19 N
i -power consumption for the i-th compressor station [W]; n-isentropa exponent; p a -absolute pressure [Pa]; Q s -volumetric gas flow through the i-th compressor station at standard conditions [m 3 /h]; Z s -gas compressibility factor for the suction [-]; T s -temperature for the suction [K]; T a -absolute temperature [K]; η-adiabatic efficiency; p d,i -discharge pressure for the i-th compressor station [Pa]; p s,i -suction pressure in for the i-th compressor station [Pa]; k-number of switch-on compressor stations in the system[-].
-volume of gas in the network [m 3 ]; p av,i -average pressure in the i-th pipeline [Pa]; p x -pressure at the sending node of the pipeline [Pa]; p y -pressure at the receiving node of the pipeline [Pa]; Z i -gas compressibility factor for the i-th pipeline [-]; T-gas temperature in the pipeline [K]; D i -internal diameter of the i-th pipeline [m]; L i -length of the i-th pipeline [m]; m-number of pipelines in the system [-].
p min -minimum nodal pressure [Pa]; p max -maximum nodal pressure [Pa]; p-pressure value at node [Pa]; w max -maximum gas velocity in the pipeline [m/s]; w-gas velocity in the pipeline [m/s].

Figure 2 .
Figure 2. Structure of network 1.The network consists of: number of pipelines-17; number of nodes-20; number of sources-1; number of compressor stations-2; source pressure-5000 kPa.

Figure 2 .
Figure 2. Structure of network 1.The network consists of: number of pipelines-17; number of nodes-20; number of sources-1; number of compressor stations-2; source pressure-5000 kPa.

Table 1 .
Pressure drop, volumetric flow rate for real and standard conditions.

Table 2 .
Pipeline geometry of network 1.

Table 8 .
Pipeline geometry of network 3.

Table 2 .
Pipeline geometry of network 1.

Table 3 .
Nodal input data of the network 1.

Table 4 .
Optimization results for network 1.Pipeline geometry of network 2 is presented in Table5.

Table 5 .
Pipeline geometry of network 2.

Table 6 .
Nodal input data of network 2.

Table 7 .
Optimization results for network 2.

Table 8 .
Pipeline geometry of network 3.

Table 9 .
Nodal input data of network 3.

Table 10 .
Optimization results for network 3.

Table 11 .
Power consumption for each optimization scenario.

Table 12 .
Comparison of optimization results-24 h scenario.

Table 13 .
Comparison of optimization results-1 year scenario.

Table 13 .
Comparison of optimization results-1 year scenario.