Document Type: Original Article
Authors
Chemical Engineering Department, Engineering Faculty, Ferdowsi University of Mashhad, Mashhad, Iran
Abstract
Keywords
Main Subjects
1. Introduction
Gas transport pipelines is one of the most important methods in transporting gas from gas wells to consumers. Since gas pressure is dropping in long distances, compressor stations should be installed in appropriate places. In view of consumption variation of consumers, the producer should inject an appropriate amount of gas in pipeline and control the output pressure of compressor stations. Gas injection flow and compressor stations’ output pressures are two important parameters in gas pipeline operation. For better operation, the consumption rate of all nodes, pipeline dynamics and delay times should be identified. Current pipeline control involves an offline manual control by considering supply and demand all over the gas network. The main objective of researchers is to design a control scheme capable of minimizing the energy consumption while keeping the system in the safe region. The safe region conditions are involved in are critical points such as minimum delivery pressure, minimum and maximum compressor flow rate. (Marques & Morari, 1988).
The first step in designing a model based control scheme is the modeling step. After that, different control schemes could be surveyed to choose the best one. The mathematical model of an unsteady state gas flow includes several partial differential equations that depend on spatial coordinates and time (Ke & Ti, 2000).
An isothermal gas flow was modeled without ignoring any terms in momentum equation (Osiadacz, 1987). Additionally, an unsteady state gas flow in pipeline was modeled assuming a steady state heat transfer term and the constant compressibility factor (Tentis et al., 2003). Their work showed a vital difference between pressure profiles of isothermal and nonisothermal models. In another work in 2009, Gonzales et al. modeled an unsteady state flow in pipeline network with the help of MATLAB^{®/}SIMULINK^{®} software (HerranGonzales et al., 2009). They have developed two simplified models and solved equations by using CrankNicolson algorithm and the method of characteristics.
The Energy Information Administration of USA reported that nearly half of the natural gas price for residential customers came from transportation costs (EIA, 2007). This created a great motive for researchers for reducing high operational cost of gas transport pipelines. One of the early attempts to develop a rational control policy is represented by (Batey et al., 1961), in which some rules of thumbs are presented for operating the system with low energy consumption. In 1988, Marqués and Morari presented a moving horizon optimization formulation, in which a dynamic simulator computed optimal operating profiles for a single source pipeline network (Marques & Morari, 1988). Moreover, compressor performance was optimized by a dynamic optimization formulation in order to follow a desired line pack profile (Abbaspour et al., 2007).
Since 2001, predictive control has been used in control and optimization of pipeline operation. In this year, a linear predictive control for industrial oxygen distribution pipeline networks was developed (Zhu et al., 2001). In a recent study, a nonlinear predictive control was used for the operation costs of the optimization of transport pipeline network (Gopalakrishnan & Biegler, 2013). In contrast with other set point tracking predictive controllers, their controller had an operation cost function and tried to reduce it. They reported that an economic predictive control could greatly reduce the operational costs of pipeline.
In this article, the pipeline network was over 300 Km long from Khangiran refinery to Farooj compressor station in Iran. This pipeline network, after necessary simplifications, includes 11 consumers, one compressor station in the middle of the line and one refinery at the beginning of the line as a producer. The authors developed a precise dynamic gas pipeline network model based on continuity, momentum and energy balances. Afterwards, three control schemes were designed for this pipeline and their performances were evaluated.
2. Dynamic Modeling of Gas Pipeline Networks
Figure 1 is a schematic of gas pipeline considered in this paper. In this figure nodes 1 and 8 represent refinery and compressor station, respectively. Other nodes are consumers. The pipeline network is divided into two sections: section one consists of node 2 to 7 and their pipelines. The rest of the pipeline network, after node 8, is called section two. The diameter, length and elevation change of each pipe segment is shown in Table 1. The conditions of the inlet gas and the flow rate of each consumer are shown in Table 2.
Figure 1. Pipeline Network
Table 1. Pipe Segments Characteristics
Elevation Change (m) 
Length (km) 
Pipe Thickness (in) 
Outside Diameter (in) 
Pipe Segment 
627 
49.5 
0.406 
36 
12 
224 
21.5 
0.406 
36 
23 
182 
68.8 
0.406 
36 
34 
13 
10 
0.406 
36 
45 
15 
7 
0.406 
36 
56 
19 
6 
0.375 
30 
67 
92 
19 
0.375 
30 
78 
89 
15.6 
0.375 
30 
89 
31 
42.5 
0.375 
30 
910 
101 
20.5 
0.375 
30 
1011 
88 
35.7 
0.375 
30 
1112 
21 
24.5 
0.375 
30 
1213 
Table 2. Inlet and Outlet Conditions of Pipeline Network
68.3 
Pipeline injection pressure (bar) 
38 
Pipeline injection temperature (^{o}C) 
0.014 
Consumption flow of node 2 (MMSCMD) 
7.200 
Consumption flow of node 3 (MMSCMD) 
2.768 
Consumption flow of node 4 (MMSCMD) 
8.953 
Consumption flow of node 5 (MMSCMD) 
7.790 
Consumption flow of node 6 (MMSCMD) 
2.886 
Consumption flow of node 7 (MMSCMD) 
0.395 
Consumption flow of node 9 (MMSCMD) 
0.504 
Consumption flow of node 10 (MMSCMD) 
0.019 
Consumption flow of node 11 (MMSCMD) 
0.603 
Consumption flow of node 12 (MMSCMD) 
10.589 
Outlet flow of pipeline, node 13 (MMSCMD) 


2.1. Pipeline Mathematical Model
The main equations of the dynamic model of gas pipeline are equations of mass, momentum and energy balances, which are expressed as follows:
Mass balance equation
(1)
Momentum balance equation
(2)
Energy balance equation of gas
(3)
Energy balance equation of pipe wall
(4)
The main assumptions in the above model are a) neglecting the radial changes of variables (radially lumped and axially differential), b) neglecting the conduction terms in momentum and energy equations and c) constant properties of pipe wall (c_{pw} and r_{w}). By solving equations 1 to 4 simultaneously, one can achieve velocity (flow rate), pressure and temperature profiles along the pipeline according to time. Required initial conditions are steady state values of temperature, pressure, and flow rate along the pipeline which are calculated by steady state assumptions in the model equations. Boundary conditions are pressure, flow rate, and temperature changes by time in the boundaries of pipeline (inlet, outlet and consumers). Each pipeline has one inlet and one outlet that needs three boundary conditions on pressure, flow rate, and temperature. For each consumer, one more boundary condition is required. These boundary conditions can be defined at inlet or outlet. In our simulation, pressure and temperature are defined at inlet, and flow rate is defined at outlet and consumers (see Table 2).
Compressor energy consumption is formulated as
Where BHP is isentropic power (hp), H_{i} is isentropic head (lb_{f} ft/lb_{m}), m is mass flow rate (lb_{m}/s) and ƞ_{c} is isentropic efficiency. H_{i} is calculated continuously with aim of operating conditions and characteristic curve of compressors. The more BHP means the more energy consumed in compressor stations for compressing natural gas. By recording the BHP of compressor stations with time, one could find which control scheme used less energy to deliver naturel gas to consumers and reject disturbances.
2.2. Model Numerical Solution
In order to solve the equations 1 to 3 it was necessary to use numerical solution of partial differential equations methods. To solve equations with method of line (Pregla & Pascher, 1989), they were formed as a set of ordinary differential equations (Kumar, 1987; Chaczykowski, 2010):
for j = 2 to N+1 (6)
for j = 2 to N+1 (7)
for j = 2 to N+1 (8)
where j is the spatial coordinate discretization section index, and P_{j}, T_{j} and m_{j}are pressure (psia), temperature (R) and mass flow rate (lb_{m}/s) of gas at jth discretized section of pipe, respectively. Δ_{x}(P_{j}), Δ_{x}(m_{j}) and Δ_{x}(T_{j}) are approximate differentiation formulas based on backward difference method. A step size of 1640 ft (0.5 km) is used for xcoordinate in these approximate differentiation formula. The value of parameters involved in the model equations (Eq. 48) or the method used for calculating these parameters are shown in Table 3. The model equations are solved in MATLAB/SIMULINK using integration routine ODE15s (Natick, 1997).
Table 4 shows a comparison between actual and simulated node’s pressure at steady state condition. As seen, most of the node’s simulated pressures are in accordance with actual pressures.
Table 3. Model Parameters and Coefficients
Parameter 
Description 
Formula or Value 
Reference 
J 
Energy conversion factor Btu / (psia ft^{3}) 
0.185 
 
g_{c} 
Force conversion factor lb_{m} ft / (lb_{f} s^{2}) 
32.174 
 
g 
Acceleration of gravity (ft/s) 
32.174 
 
Mass universal gas constant psia ft^{3} / (lb_{m} R) 
0.669 
 

c_{p} 
Specific heat of gas Btu / (lb_{m} R) 
Based on the method described in chapter 3 

z 
Compressibilty factor of gas 
Standing and Katz’s chart 

A 
Inside cross sectional area of pipe (ft^{2}) 
 

f 
Moody friction factor 

h_{i} 
Inside heat transfer coefficient, Btu / (s ft^{2} R) 

h_{o} 
Outside heat transfer coefficient, Btu / (s ft^{2} R) 
1/(R_{pipe }+ R_{soil}) where


d_{i} 
Inside diameter of pipe, ft 
2.932 or 2.437 
 
d_{o} 
Outside diameter of pipe, ft 
3 or 2.5 
 
d_{s} 
Depth from top of soil to pipe centerline, ft 
4.781 or 4.531 
 
k_{s} 
Thermal conductivity of soil, Btu / (s ft R) 
(Perry & Green, 1997), Table 2382 

k_{w} 
Thermal conductivity of pipe (mild steel), Btu / (s ft R) 
 

c_{pw} 
Heat capacity of pipe (mild steel), Btu / (lb_{m} R) 
0.12 
(Perry & Green, 1997), Table 2219 
r_{w} 
Mass density of pipe (mild steel), lbm / ft^{3} 
487 
(Perry & Green, 1997), Table 2118 
k_{g} 
Thermal conductivity of gas Btu / (s ft R) 
 

T_{s} 
Temperature of soil (R) 
519 
 
Table 4. Comparison between Actual and Simulated Node’s Pressure
Relative Error (%) 
Actual Pressure (bar) 
Simulated pressure (bar) 
Nodes 
1.52 
61.07 
60.14 
2 
0.85 
57.15 
56.66 
3 
1.47 
56.1 
55.27 
4 
0.66 
55.41 
55.04 
5 
0.47 
48.44 
48.21 
6 
3.79 
45.62 
47.35 
7 
2.36 
60.03 
58.61 
9 
3.68 
57.32 
55.21 
10 
5.05 
52.87 
50.2 
11 
1.06 
45.87 
45.38 
12 
3. Control Scheme Design
Pipeline pressure control is required for supplying the consumers’ gas needs. In this article, three control schemes, model predictive control, simple PI, and selective PI are surveyed. The control schemes regulate pipeline pressures near their set points for better consumers’ supply. All three schemes try to control pressure and supply consumers in an economic way without violating physical constrains of production and distribution tools.
3.1 Simple PI Controller Scheme
In this scheme, two simple PI controllers were used to regulate pressures of nodes 7 and 12 at the end of each section of pipeline. Changing the consumption rate of a node will result in the pressure change of that node. This pressure change will be sensed in the controlled node with a dynamic that depends on the node distance to the end of pipe section. The controller regulates output pressure of compressor station or refinery to reject this disturbance. This scheme showed in Figure 2.
3.2. Selective PI Controller Scheme
There are only two PI controllers in the simple PI scheme for pressure control. In that scheme, if a disturbance occurs in the middle of pipeline, its effects reach the controlled node by a long dynamic and the controller responds with a huge delay accordingly. The other problem of that scheme is that if the consumption rate of an important node in the middle of the pipeline is disturbed, its pressure will change dramatically before the controller can respond. These reasons establish a motive to put some controllers in the middle of the line to avoid these problems. In the selective controller scheme, nodes 3 and 7 assign refinery injection pressure and nodes 10 and 12 assign the output pressure of compressor station. It should be noticed that consumers in the second section of pipeline indirectly control the refinery injection pressure. When a disturbance occurs in the second section, the controller regulates the output pressure of compressor station which leads to a change in the inlet pressure of compressor station. This change measured is in node 7 and following that, the controller changes the refinery injection pressure. This scheme showed in Figure 3.
Figure 2. Simple PI Controller Scheme
Figure 3. Selective PI Controller Scheme. HS: High Selector
In simple and selective PI schemes SIMC (Skogestad, 2003) method was used for PI controllers design. By considering PI controller in equation 9, controller parameters are shown in Table 5. It should be noticed that controllers for node 7 and 12 are common between simple and selective PI schemes.
Table 5. PI Controller Parameters
Nodes 
K_{C} (%/%) 
K_{I} ((%/%)/s) 
3 
0.9087 
0.0003368 
7 
0.4769 
0.0001135 
10 
0.9722 
0.0006093 
12 
0.6701 
0.0001809 
3.3. Model Predictive Control Scheme
Model predictive control (MPC) refers to a class of controllers that compute a manipulated variable profile with the help of a process model to optimize controller performance and avoid constrains violation over a future time horizon (Muske & Rawlings, 1993).
The main idea of model predictive controllers can be summarized as
In MPC, the vector of control variable is obtained by solving an optimization problem at current time k (MATLAB®, 2010):
w.r.t Du_{j }(k + ik)
s.t
process model (Eqs. 48) is satisfied
55.15 bar < refinery injection pressure <72.4 bar
55.15 bar < compression station outlet pressure <72.4 bar
refinery injection flow rate < 45.32 MMSCMD
Where subscript j denotes the jthcomponent of a vector, (k+ik)denotes the value predicted for time k+i based on the information available at time k and r(k)is the current sampled value of setpoint. w_{i+1}^{y} and w_{i}^{D}^{u}are nonnegative weights for the corresponding variable. The smaller w, the less important is the behavior of the corresponding variable to the overall performance. In our work, the MPC toolbox of MATLAB® (R2011a) was used.
MPC is a centralized controller in which controlled variables, set points, and flow rates of four nodes are inputs and manipulated variables are outputs. MPC controller calculates manipulated variables, which are refinery injection pressure and compressor station output pressure, to reject the node’s pressure disturbances. The MPC scheme is shown in Figure 4 and the value of its parameters are presented in Table 6.
Figure 4. MPC Controller Scheme
Table 6. MPC Designed Parameters
Control Interval (Sampling time) 
50 s 
Prediction Horizon 
600 sampling number 
Control Horizon 
5 sampling number 
Node 3 weight 
1 
Node 7 weight 
2.5 
Node 10 weight 
1 
Node 12 weight 
2.5 
The MPC toolbox in MATLAB software supports linear time invariant model formats such as transfer function and state space models. These models can be imported by user or could be calculated by MATLAB linearization function based on the simulated nonlinear model. In the current work, the second method was used.
4. Results and Discussion
Three designed control schemes were applied to the simulated gas pipeline network. To compare the controller’s performance, three experiments were performed. In these experiments, the controller’s responses to reject known disturbances were monitured. Known disturbances were consumption rate changes in controlled and uncontrolled nodes (nodes 2 to 12). Since the distances between uncontrolled nodes and controlled nodes were not much, the closedloop responses were approximately identical and responses of uncontrolled nodes were not shown. To decide which performance was the best, the Integral of the Square of the Error (ISE) for all nodes and the consumption of energy in compressor station were considered. In addition for safe gas distribution a minimum pressure must be guaranteed at all nodes. This minimum was selected equal to 590 bar. This means that the pressure of all nodes in all times must be greater than 590 bar.
In the first experiment, an increase of 60% was applied to the consumption rate of node 3 by a dynamic with the time constant of 1.38 hour. Figure 5 shows the controllers responses. As can be seen, the MPC controller rejected this disturbance quite effectively. Simple PI and selective PI had 9600% and 5700% more ISE than MPC, respectively. They also consumed more energy in compressor station, which was 105% for simple PI and 62% for selective PI more than the required energy for MPC. In addition, the overshoot, undershoot and the settling time in the case of using MPC is obviously lower than other schemes.
Figure 5. Node 3 Consumption Increased by 60%
In the second experiment, an increase of 60% was applied to the consumption rate of node 7 by a dynamic with the time constant of 1.38 hour. Figure 6 shows the controllers responses. As can be seen, the MPC controller has sligthly better performace than simple PI controller, the performance of selective PI controller was the worst case. Simple PI and selective PI had 37% and 209% more ISE than MPC, respectively. They also consumed more energy in compressor station which was 3% for simple PI and 20% for selective PI more than required the energy for MPC. In addition, the overshoot and undershoot in the case of using MPC was obviously lower than other schemes, but the settling time was almost higher than other schemes.
In the third experiment, an increase of 30% was applied to the outlet flow rate of gas pipeline (node 13) by a dynamic with the time constant of 1.38 hour. This node was at the end of pipeline and its flow rate was relativtely high (see Table 2). Figure 7 shows the controllers responses. As can be seen, in this experiment the selective PI had better perfomaces than other schemes. Unlike the previous two expriments, the overshoot, undershoot, and the settling time in the case of using MPC were obviously higher than other schemes. This means that the selected structure and constraints for MPC were not suitable for this huge disturbance.
Figure 6. Node 7 Consumption Increased by 60%
Figure 7. Node 13 Consumption Increased by 30%
5. Conclusions
Safety and profitability are necessary for the effective control of large scale natural gas transport pipeline. Although pipeline networks were a characteristic feature of the gas transport industry, the application of advanced control to such a system was not common. In this paper, a dynamic model based on continuity, momentum, and energy balances were selected for gas pipeline network from Khangiran refinery to Farooj compressor station (in Iran). A comparison between the simulated and actual pressures at steady state conditions did not show a considerable difference. Next, three controller schemes, MPC, simple PI, and selective PI were applied to the simulated pipeline. All schemes tried to control the pressure and supply consumers in an economic way without violating the physical constrains of production and distribution facilities. The withdrawal rate changes of controlled and uncontrolled nodes showed that the MPC scheme had better performances than other schemes in disturbance rejecting of all nodes except the last node (outlet flow rate of gas pipeline). However, it should be noticed that the MPC controller had a more complicated structure and difficult design procedure than PI controllers.