Oil & Gas Science and Technology - Rev.

IFP Energies nouvelles 74, 47 (2019) Available online at:

Ó J. Shi et al., published by IFP Energies nouvelles, 2019 ogst.ifpenergiesnouvelles.fr
https://doi.org/10.2516/ogst/2019016

REGULAR ARTICLE

Application of Particle Swarm Optimization (PSO) algorithm for

Black Powder (BP) source identification in gas pipeline network
based on 1-D model
Jing Shi*, Ahmed Al-Durra, Imad Matraji, Khaled Al-Wahedi, and Mohamed Abou-Khousa

Khalifa University, Sas Al-Nakhl Campus, Abu Dhabi, United Arab Emirates

Received: 20 August 2018 / Accepted: 18 March 2019

Abstract. Black Powder (BP) is a worldwide challenge that spans all stages of the natural gas industry from
the producing wells to the consuming points. It can endanger the pipeline operations, damage instruments and
contaminate customer supplies. The formation of BP inside natural gas pipeline mainly results from the corro-
sion of internal walls of the pipeline, which is a complex chemical reaction. This work aims to develop a novel
algorithm for BP source identification within gas pipelines network based on a 1-D model of BP transport and
deposition. The optimization algorithm for BP source identification is developed based on the well-known Par-
ticle Swarm Optimization (PSO) algorithm, which can solve constrained optimization problems. By applying
this optimization algorithm on the gas transmission pipeline network, the BP source at different junctions could
be identified and quantified simultaneously. Extensive simulation studies are conducted to validate the effec-
tivity of the optimization algorithm.

1 Introduction certain amount of water, providing a necessary condition

for corrosion. The presence of BP in gas pipeline network
With the wide use of natural gas, it is challenging to main- threatens the safety of gas industry and also has lots of
tain operational efficiency and safety for gas transmission undesirable influence on the operating companies all over
pipeline network, which is a complex system with pipeline the world (literature reviews of Khan and Alshehhi,
length varying from hundreds to thousands of kilometers 2015). First of all, the gas quality will be greatly affected
(Banda et al., 2006; Tobin and Shambaugh, 2006). Black by intolerable concentration of solid particles. Secondly,
Powder (BP), occurring in both liquid and gas pipelines, gas transmission cost will gradually increase, resulting from
is the name given to the black particulates and sediment raised pressure drop due to changes of internal wall rough-
found in pipelines. It is mainly generated due to the chem- ness and diameter over a long time of operation. Thirdly,
ical reaction of H2S, water and iron, resulting in a mixture the contaminations of BP could cause compressor failures,
of fine particle corrosion product and other solids, such as erosion of control valve and instrument clogging, etc.
sands, clays, metal or construction debris, and liquid hydro- (Baldwin, 1998). Finally, the cleaning/removal of BP and
carbons chemically incorporated with any quantity of iron frequent replacement of customers’ cartridge filter elements
sulfide, iron carbonate and iron oxide contamination (Khan could increase the expenses each year.
and Alshehhi, 2015; Khan et al., 2015; Sherik et al., 2008; To reduce the influence of BP on gas industry, many
Sherik, 2008; Trifilieff and Wines, 2009). Although there companies attempted to manage and control the BP in
exist various compositions of BP, they possess some com- gas transmission pipeline network. In general, the existing
mon characteristics, for example, adsorption, high specific methods could be classified into two different aspects (Al-
gravity and difficult to clean. In terms of BP generation, Qabandi et al., 2015; Cattanach et al., 2011; Khan et al.,
they could come from the gas source in the gas field and also 2015; Trifilieff and Wines, 2009; Tsochatzidis and Maroulis,
result from movement of the upstream point with gas flow, 2007): (1) BP removal and (2) prevention methods. The
which originates from the corrosion of internal wall of any approach of removing BP from gas transmission pipeline is
pipeline. Specifically, corrosion could only occur in the pipe- successful in several respects, including increasing opera-
line with the flow of “wet” gas (Beavers and Thompson, tional safety, altering operational parameters, increasing
2006; Baldwin, 1998), which is natural gas mixed with a operational efficiency and facilitating effective corrosion
inspection. If BP problem is not serious (small quantities),
* Corresponding author: jing.shi5@student.unsw.edu.au mechanical cleaning is widely used for most of the pipelines.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0),
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
2 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019)

Normally, mechanical scrapers are used to clean pipeline particles in gas transmission network, where the authors
wall, such as mechanical pigs are deployed into a pipeline discussed behaviors of particles with different diameters.
to scrape debris from pipeline wall and remove BP. Moreover, two different deposition models (Fan and
Although mechanical scrap is efficient to keep the pipeline Ahmadi, 1993; Wood, 1981) were compared to calculate
in fairly clean condition, frequent cleaning could damage the bed height, which was validated by Discrete Phase
the pipeline wall by exposing fresh steel surface under Model (DPM) based on CFD software.
“wet” gas, resulting in excessive corrosion and BP genera- The main contribution of this work is to develop a tool
tion. In addition, mechanical cleaning is not effective with to identify the BP source in gas pipeline network, which is
complex BP formation, which calls for a combination of modeled as tree-shaped gas transmission network with BP
chemical cleaning. It has been practiced that chemical clean- dynamics of motion and distribution along each pipe. BP
ing, acting as liquid soap, could significantly improve the source identification is formulated as a constrained opti-
efficiency of BP removal combined with mechanical clean- mization problem, which is solved by the Particle Swarm
ing, such as pigging operations, which acts as brush (Trifil- Optimization (PSO) techniques (Delice et al., 2017; Eber-
ieff and Wines, 2009). Separators and cyclo-filters are also hart and Shi, 1998; Kennedy and Eberhart, 1995; Min
widely installed to reduce the BP concentration, where these et al., 2017). PSO is a popular stochastic optimization tech-
devices could physically knock out the BP particles in the nique with some features and advantages compared to
gas stream. Then, the BP particles are collected at the bot- other optimization algorithms such as Ant Colony Opti-
tom in a collection hub. However, this method only applies mization (ATO) and Genetic Algorithm (GA) (Wiak
to gas stream with a high concentration of solid particles and et al., 2008; Saravanan, 2006). These features and advan-
also relatively large particle size (more than 10 microns). tages include (1) taking real numbers as particles; (2) few
The removal approach is practiced as a good option to pro- parameters need to be tuned and (3) simple implementation
tect the downstream operations, but it has several disadvan- and effective global search capability, hence it could be a
tages (Trifilieff and Wines, 2009). These methods are not good candidate. In fact, PSO technique has been widely
single trial, but frequent application, which account for a used in oil and gas industry, for example, Wu et al.
large amount of expense. Also, these solutions cannot (2014) optimized the operation of trunk natural gas pipeli-
address the location of BP formation, that is to say, they nes via PSO based algorithm, and Madoliat et al. (2017)
are afterward remedy. In addition, chemical cleaning needs also successfully applied PSO to the transient analysis of
subsequent handling procedures, which could be costly natural gas pipeline. The PSO algorithm is used to solve
and challenging if the chemical disposal is toxic or harmful a formulated optimization problem for similar applications.
to the environment. To compensate the above drawbacks In particular, the basic PSO algorithm has been improved
of removal approach, gas operators have an alternative solu- by incorporating the inertia-adaption technique to solve a
tion to prevent the occurring of BP generation, which is con- constrained nonlinear optimization problem in Wu et al.
sidered as a consequence of corrosions at the internal wall of (2014). Then, it was used to solve a more complex problem
gas transmission pipeline (Cattanach et al., 2011; Trifilieff in Madoliat et al. (2017), where the solution of nonlinear
and Wines, 2009). To prevent the occurring of corrosions, PDE flow equations could be obtained simultaneously.
the inner wall of gas transmission pipeline is normally coated The paper is structured as follows. Section 2 will briefly
with high solids solvent-based epoxy polyamine films, which discuss the 1-D approach for modelling BP particles in gas
is used to protect the inner surface of pipelines. This networks, which is detailed in Filali et al. (2016). A general
approach is practiced to be cost effective; however, it is dif- structure of tree-shaped model of gas transmission pipeline
ficult to be applied to buried pipelines. Reducing water con- network is given. In Section 3, PSO-based optimization
tamination is another approach to prevent the occurring of algorithm for BP identification is presented. In Section 4,
pipeline corrosions, which is based on the philosophy of extensive simulation results and discussions of the optimiza-
internal corrosions, are largely related to the case of “wet” tion algorithm for BP source identification are presented.
gas. Namely, this approach is moisture control. Finally, conclusions are given in Section 5.
Although several removal and prevention methods have
been put into practice, the location of BP generation is still
unpredictable and there is no research on the identification 2 Modelling for dynamics of BP particles in
problem of BP source until now. The main reason is the gas transmission pipeline
unknown information about the BP generation and where
it could be generated. The developed dynamic models for In this section, a schematic of tree-shaped gas transmission
BP transport and deposition are limited for gas transmis- pipeline network will be built based on a set of pipeline con-
sion pipeline network, which in turn restricts the study of nection rules. In addition, the methodology of a 1-D
model-based methods for BP identification. It is well known approach for modelling the dynamics of transport and
that full three-dimension (3-D) Computational Fluid deposition of BP particles in the gas pipeline network will
Dynamics (CFD) simulation software is specialized for mul- be explained (Filali et al., 2016).
tiphase flow modelling; however, it seems impossible to sim-
ulate the 3-D dynamics of large pipeline networks (e.g. 2.1 A tree-shaped model of gas network
100 km), which is under studied in these works (Kharoua
et al., 2015, 2017). Filali et al. (2016) proposed a 1-D In gas industry, the natural gas pipeline network is a highly
approach for modelling transport and deposition of BP integrated transmission and distribution grid that could
 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019) 3

transport natural gas from its origin to any position of high @2C @C
gas consumption demand. In many cases, the natural gas Ddiff U : ð2Þ
@x 2 @x
produced from a well has to be transported on a very long
distance to the point of use. In order to maintain effective For steady state solution, equation (1) could be simpli-
distribution of natural gas, the gas transmission system is fied as:
extensive, consisting of complex pipeline topology. How-
@C
ever, the entire gas transmission network could be sepa- U ¼ bdep C þ ngen ; ð3Þ
rated into several sub-network with single gas well. This @x
simplification also contributes to the investigation of BP where, bdep is the deposition rate and it could be calcu-
identification, which will be included in following sections lated based on two methods, for example, Wood (1981),
of this paper. In this case, the gas transmission pipeline net- Fan and Ahmadi (1993), and ngen is the generation rate
work is properly generalized as tree-shaped model (Babon- of the BP inside the pipeline and it is considered as an
neau et al., 2012; Shiono and Suzuki, 2016), which could be unknown parameter.
built based on following rules for pipe connection:

1. Gas transmission pipeline network without loop;

2. Junction connected with three pipes;
3 PSO-based optimization algorithm
3. One-way flow of gas stream. for BP identification
In this section, the objective of BP identification is pre-
Based on these connection rules, a schematic of tree- sented as a constrained optimization problem, which will
shaped gas pipeline transmission network is given in Fig- be solved by a PSO-based optimization algorithm. The
ure 1. There is only one source for gas supply, and each PSO algorithm is a population-based evolutionary search
junction has three pipe connection, including two main algorithm inspired by social behavior of animals such as
pipes and one branch pipe. The arrow on the pipe shows bird flocking. Basically, each particle in the swarm has a
the flow direction of gas stream. The customer is located position and velocity, with its position representing a candi-
at the end of each branch pipe. The sources, main pipes, date solution in the multi-dimensional solution space and
branch pipes and customers are denoted as S, P, B, and velocity indicating moves from one position to another.
M respectively. A fitness function is defined to evaluate each particle until
certain convergence criteria is satisfied. During the search-
ing process, the particle with fitness value will be selected
2.2 One-D approach of modelling BP particle dynamics as local/global best particle. The updated equations of
position and velocity are given in the literature (Marini
In this section, the 1-D approach for modelling dynamics of and Walczak, 2015) as follows:
BP particles in gas pipeline (Filali et al., 2016) is described.
It is well known that the CFD software is popular to simu-
x ij ðt þ 1Þ ¼ x ij ðt Þ þ v ij ðt þ 1Þ; ð4Þ
late the dynamics of multi-phase flow, which could give a
solution to the simulation of BP particles in the gas stream.
However, the computation load of CFD is extremely heavy

for gas pipeline network, which could extend to hundreds of v ij ðt þ 1Þ ¼ v ij ðt Þ þ c1 r 1j ðt Þ y ij ðt Þ  x ij ðt Þ
kilometers. In this case, Filali et al. (2016) proposed a sim-

plified 1-D approach of modelling, where the flow of gas is þ c2 r 2j ðt Þ y^j ðt Þ  x ij ðt Þ ; ð5Þ
continuously mixed with BP particles, and the behavior
of particle movement is modelled based on the dusty gas where, xij(t) denotes the position of particle i in dimension
assumption and the usage of analytical solutions of steady j at time t, with j = 1, 2, 3, . . ., nx, and nx is the dimension
scalar advection/reaction equation. of the solution space. The updating of position is
The governing equation for gas stream mixed with solid calculated by adding the velocity vij(t) to current position.
particles is given by the following one-dimensional advec- yij(t) is the personal best solution in dimension j for
tion–diffusion-reaction equation: particle i, while y^j ðt Þ represents the best solution
found so far. c1 and c2 are positive acceleration factors
used to scale the contribution of cognitive and social com-
@C @C @2C ponent. r 1j ; r 2j 2 ½0; 1 are random variables with normal
þU ¼ Ddiff 2 þ S; ð1Þ
@t @x @x distribution.
The structure of this section is outlined as follows:
where C is the particle concentration in the gas flow, U is firstly, the objective of BP identification is formulated,
the average velocity of gas flow inside the pipeline, Ddiff is which is developed based on the 1-D model of BP particles
the diffusion coefficient, S is a term to describe the depo- transport and deposition. Secondly, the BP source could be
sition, pickup or generation of solid particles. For fully identified and quantified by sole running of the PSO-based
developed turbulent flow, the axial diffusion term optimization algorithm. Finally, a proof will be provided to
Ddiff @@xC2 is negligible compared to the dominant advection
2
support the uniqueness of the formulated optimization
term U @C @x
and the following expression is valid, problem.
4 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019)

M2 M5 M7

B2 B5 B7
S1 P1 S2 P2 P3 S4 S5 P5 P6 S7 P7 B8
M8
B1 B3 P4 B6
S3 S6 S8
B4
M1 M3 M6

M4

Fig. 1. Schematic of gas transmission pipeline network.

3.1 Formulation of BP identification problem weighted sum of errors (e.g. absolute error) as its cost
function.
It is noted that additional unknown sources of BP genera- Xm  
 b 
tion are considered in this study. Along with that, the min k k M k  M k ðS 1 ; S 2 ; S 3 ; . . . ; S i Þ; ð6Þ
objective of this paper is to identify the BP sources and S 1 ; S 2 ; S 3 ; ...; S i
k¼1

quantify the generated BP concentration using several mea- 

surements of BP concentration at the end node of each S n > 0; n ¼ 1; 2; 3; . . . ; i
s:t: ;
branch pipe, which is a point of natural gas use (client). kk > 0
As most of the gas transmission pipeline network is buried
underground, therefore the sensor used for BP concentra- where, S1, S2, S3, . . ., Si represent the estimation of BP
tion measurement is only installed at each client (Abou- concentration at sources (S1 is the main source, and S2,
Khousa et al., 2015). The presence of an optimization algo- S3, . . ., Si are the additional source). m is the number of
rithm for BP source identification is developed based on a 1- clients and BP measurement points. Mk is the measure-
D model, which has several assumptions (Filali et al., 2016), ment of BP concentration at the client k, which is mea-
such as, one-way flow of gas stream, constant velocity of gas sured by BP concentration sensor located at the end of
flow at each pipe, average mixture of solid particles and gas each branch pipe and it is simulated using the 1-D model
stream, fixed size/diameter of solid particle, etc., but some presented in the previous section. M ^ k is the estimated
additional assumptions should be also taken into considera- value of BP concentration, kk is the weighting parame-
tion. The accuracy of BP source identification is largely ters for each system input, which could be selected as
restricted because of limited measurement of BP concentra- inverse of each measurement of BP concentration at cor-
tion along the pipe, which mainly results from two aspects: responding client. The optimization algorithm for BP
(1) lack of reliable measuring device that can be installed source identification is presented as a flowchart in
underground; (2) the “additional” BP generation from the Figure 2.
internal walls is unpredictable, that is, the occurrence of
chemical corrosion is highly random and the location is dif-
ficult to be identified. Therefore, the gas supply station is 3.2 Discussion on unique solution of
assumed in this paper as the main BP generation source, BP optimization problem
while the junctions stand for additional source, which essen-
tially represent the amount of BP generated along its down- Proposition 3.1 Based on the given assumptions in
stream. Herein, two assumptions are given for the problem Section 3.1, the solution of the formulated optimization
formulation of BP identification: problem is unique if and only if the number of measurement
is equal to the number of unknown BP source.
(a) BP source is located in junction only. The proof is referred to the Appendix A.
(b) Additional BP particle generation occurs in the main
pipe only.
4 Simulation studies
According to Figure 1, there is one main BP particle
source and several unknown additional sources located at In this section, the optimization algorithm for BP source
each junction, which are the connection points for every identification is applied to a real gas transmission pipeline
two main pipes and one branch pipe. Let us consider the network, including 15 pipes, 16 nodes, and 8 junctions.
gas transmission pipeline network as the system. Measure- The geometric parameters of this gas pipeline network are
ments of BP concentration at each client are the system practical data from an existing network. The topology of
inputs and estimations of BP concentration at each source this network is given in Figure 1, as well as its geometric
are the system outputs. Therefore, the following parameters in Table 1. The remaining parts are extensive
constrained optimization problem (Aguirre et al., 2007; simulation studies, including BP identification algorithm
Cagnina et al., 2008; Hu and Eberhart, 2002; Liu, 2008; applied on gas network with perfect measurements, and
Parsopoulos and Vrahatis, 2002) can be formulated with sensitivity studies on this algorithm. The sensitivity studies
 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019) 5

(a) start Table 1. Geometric parameters and flow conditions for

gas network.
Set the BP Source at each junction, as
variables to be optimized
Pipe Pipe Pipe Mass flow Gas Start End
Initialization of the position No. length diameter rate, Q density node node
and velocity of each particle
(m) (inch) (kg/s) (kg/m3)
One-D Model Evaluation of the fitness of P1 10 000 36 130.092 26.3 S1 S2
Running each particle
B1 3000 10 13.364 25 S2 M1
For each particle, set:
Local best fitness = current fitness
P2 13 000 36 116.728 23 S2 S3
Local best position = current position
B2 11 000 24 14.5 26 S3 M2
Global best fitness = min(local best fitness) P3 12 000 36 102.228 25 S3 S4
B3 3000 10 13.364 26.3 S4 M3
Update velocity and position for each particle
P4 3000 36 88.864 25 S4 S5
Evaluation of the One-D Model B4 11 000 24 14.5 23 S5 M4
fitness of each particle Running
P5 11 000 36 74.364 26 S5 S6
NO B5 10 000 24 13.364 25 S6 M5
Current fitness<
Local best fitness P6 10 000 36 61 26.3 S6 S7
YES B6 3000 10 14.5 25 S7 M6
Set local best fitness = current fitness P7 13 000 36 46.5 23 S7 S8
B7 11 000 24 15.5 26 S8 M7
NO
Current fitness< B8 12 000 24 31 25 S8 M8
global best fitness

YES

Set global best fitness = current fitness

Table 2. Real value of BP concentration at each source.

NO
Stopping criteria
met? Source No. S1 S2 S3 S4 S5 S6 S7 S8
YES
BP concentration 9.5 6 7.5 5 4.5 8 5.5 0.5
end
(kg/m3)
(b)
Definition of the
structure of Gas Pipeline
at the end node of each branch pipe (client) are simulated
Initialization of Pipeline based a sequence of BP sources given as follows.
Parameters By sole running of the 1-D model, the measurements of
BP concentration at each client could be generated accord-
Concentration ing to the BP concentration at each source given in Table 2.
of BP Source Then these measurements shown in Table 3, will be saved
Running of One-D Model
and assumed as known parameters, i.e. BP concentration
measurements. This is the first step for following simulation
Concentration Distribution studies. The second step is the estimation of BP concentra-
along each Pipeline tion at each source using the measurements of BP concen-
tration. It is noticeable that the aim of BP optimization
Fig. 2. (a) PSO based BP source optimization algorithm. problem is the reverse process, where the BP concentration
(b) Flowchart of One-D model running. at each source will be estimated by application of the opti-
mization algorithm.
Simulation results of identification for BP concentration
at sources are given in Figure 4. In this simulation study,
are divided into two aspects, including investigations on
the parameters of PSO are selected as: kk ¼ M1k , n = 300,
disturbed measurements and model parameter mismatch.
m = 50, c1 = 2, c2 = 2, where, kk is the weighting parameter
Finally, some discussions are presented.
in equation (6), n is the quantity of particles, m is the max-
imum number of iteration (termination condition), and c1,
4.1 Ideal situations c2 are the acceleration constants. It is noticeable that the
PSO parameters are tuned by trial and error, which can
In this section, the optimization algorithm for BP identifica- ensure satisfactory performance with acceptable computa-
tion is validated on the gas network given in Table 1. It is tional load. PSO is a stochastic optimization technique,
noted that, initially, the measurements of BP concentration therefore the algorithm is repeated three times (blue, red
6 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019)

Table 3. Measurements of BP concentration at each client.

Client No. M1 M2 M3 M4 M5 M6 M7 M8
3
BP concentration (kg/m ) 0.3802 0.1404 0.4831 0.2604 0.1753 1.3014 0.2888 1.3857

(a) (b)

(c) (d)

Fig. 4. Estimations of BP concentration at each source, figures (a)–(h) represent sources 1–8.

and green curves in each figure) to show/compare the 4.2 Sensitivity studies
stability of the optimization results. The dotted line in each
figure represents the real value of BP concentration at each As stated in Section 4.1, the BP concentration at sources
source, which is corresponding to Table 3. In addition, the can be well estimated/reproduced by the PSO algorithm
statistic results of each test are summarized in Table 4. The when the model is completely known and there are no
estimation value and absolute error are given for each test uncertainties in measurements. However, the efficacy of this
and source, and an overall average with the Standard approach is still questionable in practice, because all the
Derivation (SD) are given for all sources. It could be models are at best only approximations for reality due to
concluded that the optimization algorithm for BP source various model errors or parametric uncertainties. This prob-
identification is satisfactory and able to achieve small esti- lem can be formulated as follows:
mation errors with perfect measurements of BP concentra-
tion at clients. b k ðh; S  Þ þ ;
M k ð hÞ ¼ M ð7Þ
 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019) 7

(e) (f )

(g) (h)
Fig. 4. Continued.

with h the vector of geometric and flow conditions param- uncertainty, where M3 and M5 are increased by 5% respec-
eters,  the measurement error, Mk the real measurements, tively. Basically, e is regarded as an additive nonzero noise
Mb k the 1-D model approximation and S* the vector of BP and no variations in h in this scenario. The tuned parame-
concentration at sources which are unknown. The aim is ters for the algorithm are the same with ideal situations.
to effectively recover S* via PSO algorithm relying on Some figures and statistic data are presented as preliminary
Mk, regardless of the presence of uncertainties in h and results. It is noted that the problem formulation in equation
. A good approach to solve this problem is associated (7) is a deterministic case, which will be generalized in the
to robust model calibration and inverse problem theory next stage of experimental validations and tackled by the
(Tarantola, 1987; Kaipio, 2008). likelihood approaches (Tarantola, 1987), where measure-
As a preliminary study, three case studies are simulated ment errors can be modelled with a probabilistic distribu-
in this section. Firstly, the measurement error on M3 and tion and confidence intervals will be computed for
M5 are simulated. Secondly, the disturbances in flow mea- estimated BP concentration at sources within Bayesian
surements are considered. Thirdly, parametric uncertainties framework.
are simulated for pipeline roughness. It is noted that each
case study is repeated three times with same parameters
to prove the algorithm’s stability and reliability. 4.2.1.1 Uncertainty in M3 measurement
In this case, the BP concentration measurement M3 is
4.2.1 Black Powder measurement uncertainty increased by 5%. The performance of identification is shown
in Figure 5, and also the corresponding statistic data is
In this section, the optimization algorithm for BP source presented in Table 5. It can be seen from Figures 5c and
identification is tested under the situations of measurement 5d, 10% and 5% estimation error for sources 3 and 4
8 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019)

Table 4. Simulation results of optimization algorithm under ideal situations.

BP source Data set
Simulation results Statistic data
Test 1 Test 2 Test 3
Estimation Abs. Error Estimation Abs. Error Estimation Abs. Error Real Avg. SD
S1 9.4994 0.0006 9.4996 0.0004 9.5007 0.0007 9.5000 9.4999 0.0006
S2 6.0004 0.0004 6.0003 0.0003 5.9997 0.0003 6.0000 6.0001 0.0003
S3 7.5000 0.0000 7.4998 0.0002 7.5036 0.0036 7.5000 7.5011 0.0017
S4 5.0005 0.0005 5.0001 0.0001 5.0133 0.0133 5.0000 5.0046 0.0061
S5 4.4994 0.0006 4.4995 0.0005 4.4521 0.0479 4.5000 4.4837 0.0223
S6 8.0001 0.0001 8.0004 0.0004 8.0283 0.0283 8.0000 8.0096 0.0132
S7 5.5000 0.0000 5.5000 0.0000 5.4950 0.0050 5.5000 5.4983 0.0024
S8 0.5000 0.0000 0.5000 0.0000 0.5041 0.0041 0.5000 0.5014 0.0019
Cost function 0.0001 0.0001 0.0037 0.0000 0.0013 0.0017

(a) (b)

(c) (d)

Fig. 5. Estimations of BP concentration at sources for 5% increase in M3 BP measurements; figures (a)–(h) represent sources 1–8.
 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019) 9

(e) (f )

(g) (h)
Fig. 5. Continued.

Table 5. Simulation results of optimization algorithm under situation for 5% increase in M3.
BP source Data set
Simulation results Statistic data
Test 1 Test 2 Test 3
Estimation Abs. Error Estimation Abs. Error Estimation Abs. Error Real Avg. SD
S1 9.4999 0.0001 9.5000 0.0000 9.4963 0.0037 9.5000 9.4987 0.0017
S2 6.0000 0.0000 6.0000 0.0000 6.0000 0.0000 6.0000 6.0000 0.0000
S3 8.2524 0.7524 8.2523 0.7523 8.2519 0.7519 7.5000 8.2522 0.0002
S4 4.4748 0.5252 4.4704 0.5296 4.4768 0.5232 5.0000 4.4740 0.0027
S5 4.5000 0.0000 4.5038 0.0038 4.5001 0.0001 4.5000 4.5013 0.0017
S6 8.0266 0.0266 8.0006 0.0006 7.9997 0.0003 8.0000 8.0090 0.0125
S7 5.4831 0.0169 5.5508 0.0508 5.4971 0.0029 5.5000 5.5103 0.0292
S8 0.5000 0.0000 0.4678 0.0322 0.5100 0.0100 0.5000 0.4926 0.0180
Cost function 0.0012 0.0030 0.0016 0.0000 0.0019 0.0008
10 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019)

Table 6. Simulation results of optimization algorithm under situation for 5% increase in M5.
BP source Data set
Simulation results Statistic data
Test 1 Test 2 Test 3
Estimation Abs. Error Estimation Abs. Error Estimation Abs. Error Real Avg. SD
S1 9.4961 0.0039 9.5000 0.0000 9.5000 0.0000 9.5000 9.4987 0.0018
S2 6.0037 0.0037 5.9952 0.0048 6.0000 0.0000 6.0000 5.9996 0.0035
S3 7.4992 0.0008 7.5030 0.0030 7.5001 0.0001 7.5000 7.5007 0.0016
S4 4.9988 0.0012 5.0000 0.0000 4.9999 0.0001 5.0000 4.9996 0.0006
S5 5.3850 0.8850 5.3828 0.8828 5.3829 0.8829 4.5000 5.3836 0.0010
S6 7.2908 0.7092 7.3189 0.6811 7.3133 0.6867 8.0000 7.3077 0.0121
S7 5.5141 0.0141 5.4965 0.0035 5.4778 0.0222 5.5000 5.4961 0.0148
S8 0.5030 0.0030 0.5002 0.0002 0.5140 0.0140 0.5000 0.5057 0.0060
Cost function 0.0019 0.0007 0.0012 0.0000 0.0012 0.0005

(a) (b)

(c) (d)

Fig. 6. Estimations of BP concentration at sources for 10% increase in Q3, figures (a)–(h) represent sources 1–8.
 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019) 11

(e) (f )

(g) (h)
Fig. 6. Continued.

Table 7. Simulation results of optimization algorithm under situation for 10% increase in Q3.
BP source Data set
Simulation results Statistic data
Test 1 Test 2 Test 3
Estimation Abs. Error Estimation Abs. Error Estimation Abs. Error Real Avg. SD
S1 9.5730 0.0730 9.5730 0.0730 9.5728 0.0728 9.5000 9.5730 0.0001
S2 6.1084 0.1084 6.1086 0.1086 6.1237 0.1237 6.0000 6.1136 0.0072
S3 7.5472 0.0472 7.5485 0.0485 7.5379 0.0379 7.5000 7.5445 0.0047
S4 5.0310 0.0310 5.0304 0.0304 5.0310 0.0310 5.0000 5.0308 0.0003
S5 4.5004 0.0004 4.4997 0.0003 4.5003 0.0003 4.5000 4.5001 0.0003
S6 8.0002 0.0002 8.0010 0.0010 8.0895 0.0895 8.0000 8.0302 0.0419
S7 5.5001 0.0001 5.4993 0.0007 5.4441 0.0559 5.5000 5.4812 0.0262
S8 0.5000 0.0000 0.5003 0.0003 0.4994 0.0006 0.5000 0.4999 0.0004
Cost function 0.0001 0.0002 0.0055 0.0000 0.0019 0.0025
12 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019)

Table 8. Simulation results of optimization algorithm under situation for 5% increase in Q5.
BP source Data set
Simulation results Statistic data
Test 1 Test 2 Test 3
Estimation Abs. Error Estimation Abs. Error Estimation Abs. Error Real Avg. SD
S1 9.5364 0.0364 9.5352 0.0352 9.5481 0.0481 9.5000 9.5399 0.0058
S2 6.0551 0.0551 6.0547 0.0547 6.0356 0.0356 6.0000 6.0485 0.0091
S3 7.5437 0.0437 7.5457 0.0457 7.4902 0.0098 7.5000 7.5265 0.0257
S4 5.0098 0.0098 5.0114 0.0114 5.0440 0.0440 5.0000 5.0217 0.0158
S5 4.4390 0.0610 4.4366 0.0634 4.4514 0.0486 4.5000 4.4424 0.0065
S6 8.0688 0.0688 8.0671 0.0671 8.0619 0.0619 8.0000 8.0659 0.0029
S7 5.4989 0.0011 5.4969 0.0031 5.6637 0.1637 5.5000 5.5532 0.0782
S8 0.5009 0.0009 0.5027 0.0027 0.3986 0.1014 0.5000 0.4674 0.0487
Cost function 0.0003 0.0006 0.0157 0.0000 0.0055 0.0072

(a) (b)

(c) (d)

Fig. 7. Estimations of BP concentration at sources for 5% increase in roughness, figures (a)–(h) represent sources 1–8.
 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019) 13

(e) (f )

(g) (h)

Fig. 7. Continued.

respectively, which are expected and acceptable. Although the model fidelity to measurement data given parametric
the inaccurate measurement occurring at M3 affects the uncertainties. Although these simulation tests are not rep-
BP concentration estimation at sources 3 and 4, and it resentative in practice where noise/uncertainties always
has negligible influence on the BP concentration estimation exist in a stochastic way, they can provide some preliminary
at sources 5, 6, and the downward sources. results at this stage.

4.2.2.1 Uncertainty in Q3 measurement

4.2.1.2 Uncertainty in M5 measurement
In this case study, the mass flow rate measurement of
In this case study, M5 is increased by 5% and the perfor-
Q3 is increased by 10%. The simulation results are
mance of the optimization algorithm is presented in Table 6.
shown in Figure 6, along with its statistic data in Table 7.
It is noted that the estimation errors of sources 5 and 6 are
It can be seen in Figure 6, the variation of flow rate Q3
15% and 10% respectively, which are high relatively com-
has a very small influence on the optimization, which proves
pared with their downward sources. The estimation error
the robustness of the optimization algorithm under
of sources 7 and 8 could be negligible.
uncertainty.

4.2.2 Mass flow rate measurement uncertainty

4.2.2.2 Uncertainty in Q5 measurement
Similarly, variations in h are simulated in this section, In this case study, the mass flow rate measurement of
where the mass flow rate in pipe 3 (Q3) and pipe 5 (Q5) Q5 is increased by 5%. The simulation results are
are increased by 10% and 5% respectively. The tuned shown in Table 8. The estimation error is less than 2%
parameters for PSO are same with the ideal situations. for each source, which is acceptance and negligible for this
One objective of this paper is to study the robustness of work.
14 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019)

Table 9. Simulation results of optimization algorithm under situation for 5% increase in roughness.
Data set
BP source Simulation results Statistic data
Test 1 Test 2 Test 3
Estimation Abs. Error Estimation Abs. Error Estimation Abs. Error Real Avg. SD
S1 9.6425 0.1425 9.6367 0.1367 9.6235 0.1235 9.5000 9.6342 0.0097
S2 6.1501 0.1501 6.1489 0.1489 6.1427 0.1427 6.0000 6.1472 0.0040
S3 7.6425 0.1425 7.6349 0.1349 7.6429 0.1429 7.5000 7.6401 0.0045
S4 5.1102 0.1102 5.1045 0.1045 5.1049 0.1049 5.0000 5.1065 0.0032
S5 4.5810 0.0810 4.5727 0.0727 4.5829 0.0829 4.5000 4.5789 0.0054
S6 8.1640 0.1640 8.1594 0.1594 8.1587 0.1587 8.0000 8.1607 0.0029
S7 5.6188 0.1188 5.6037 0.1037 5.6134 0.1134 5.5000 5.6120 0.0077
S8 0.5109 0.0109 0.5107 0.0107 0.5106 0.0106 0.5000 0.5107 1.5e-4
Cost function 0.0023 0.0046 0.0032 0.0000

4.2.3 Deposition rate uncertainty to Abu Dhabi Gas Industries Limited (GASCO) for funding
and supporting this project.
In this case study, the pipeline roughness is increased by
5%. The parameter for the algorithm is same with the ideal
case and some simulation results are shown above:
It can be seen in Figure 7 that the estimation of BP con- References
centration for each client is always higher than the expected
value, and the estimation error is shown in Table 9. This is Abou-Khousa M., Al-Durra A., Al-Wahedi K. (2015) Microwave
reasonable because the increase of roughness will result in sensing system for real-time monitoring of solid contaminants
more BP deposition, in other words, more BP should be in gas flows, IEEE Sens. J. 15, 9, 5296–5302.
generated in order to achieve the same measured BP con- Aguirre A.H., Zavala A.M., Diharce E.V., Rionda S.B. (2007)
centration as expected. However, the average estimation COPSO: constrained Optimization via PSO algorithm, Center
error is nearly 2%, which is still acceptable in the practical for Research in Mathematics (CIMAT), Guanajuato, Mexico,
situation. 30 p. Technical Report No. I-07-04/22-02-2007.
Al-Qabandi S., Khuraibut Y., Al-Inzi M. (2015) Integrated
Solution to Mitigate and Prevent Black Powder Phenomena in
West Kuwait Operation Facilities, Proc. Conference on SPE
5 Conclusion Kuwait Oil and Gas Show, 11–14 October, Mishref, Kuwait.
SPE-175291-MS.
In this paper, the problem of BP source identification is Babonneau F., Nesterov Y., Vial J.P. (2012) Design and
studied by applying a PSO-based optimization algorithm, operations of gas transmission networks, Oper. Res. 60, 1,
which is developed upon a 1-D model of BP transport 34–47.
and deposition. The 1-D model is a simplified approach Baldwin R.M. (1998) Technical assessment: “Black Powder” in
for modelling the dynamics of BP particles in gas transmis- the Gas Industry – sources, characteristics, and treatment,
sion pipeline network. A schematic of tree-shaped gas Gas Machinery Research Council Report TA, pp. 97–104.
network is proposed, which is generalized with a set of Banda M.K., Herty M., Klar A. (2006) Gas flow in pipeline
connection rules. Through some preliminary simulation networks, NHM 1, 1, 41–56.
studies, the PSO-based algorithm is validated to be a useful Beavers J.A., Thompson N.G. (2006) External corrosion of oil
technique for BP source identification. At this stage, the and natural gas pipelines, ASM Handbook 13, 1015–1025.
Cagnina L.C., Esquivel S.C., Coello C.A.C. (2008) Solving
model errors and uncertainties are assumed deterministic
engineering optimization problems with the simple con-
though they are not representative in practice. The main
strained particle swarm optimizer, Informatica 32, 3.
purpose of this paper is to discuss the application of PSO Cattanach K., Jovancicevic V., Ramachandran S. (2011)
techniques on BP source identification. More work will be Development of new corrosion inhibitor to prevent black
done together with the experimental studies in the next powder formation using quartz crystal microbalance tech-
stage, including (1) analysis of likelihood approach which nique, Proc. Conference on NACE international CORRO-
has specified probabilistic distribution to model error; SION. No. 11093.
(2) improvement of the 1-D model by considering pickup Delice Y., Aydoğan E.K., Özcan U., Ilkay _ M.S. (2017) A
case and various particle sizes. modified particle swarm optimization algorithm to mixed-
model two-sided assembly line balancing, J. Intell. Manuf. 28,
Acknowledgments. The authors acknowledge the Gas Process- 1, 23–36.
ing & Materials Science Research Centre (GRC) at Khalifa Eberhart R.C., Shi Y. (1998) Comparison between genetic
University – The Petroleum Institute, Abu-Dhabi, in addition algorithms and particle swarm optimization, Evolutionary
 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019) 15

programming VII: Proc. Conference on 7th Annual Evolu- Tobin J., Shambaugh P. (2006) The crucial link between natural
tionary Programming, pp. 611–616. gas production and its transportation to market, Proc.
Fan F.G., Ahmadi G. (1993) A sublayer model for turbulent Conference on Stages in the production of pipeline-quality
deposition of particles in vertical ducts with smooth and rough natural gas and NGLs. Energy Information Administration,
surfaces, J. Aerosol Sci. 24, 45–64. Natural Gas Annual. pp. 11.
Filali A., Khezzar L., Alshehhi M., Kharoua N. (2016) A one-D Trifilieff O., Wines T.H. (2009) Black powder removal from
approach for modeling transport and deposition of Black transmission pipelines: Diagnostics and solutions, Proc. Con-
Powder particles in gas network, J. Nat. Gas Sci. Eng. 28, ference on Pipeline Rehabilitation and Maintenance, January
241–253. 19–21, 2009, Bahrain, pp. 19–21.
Hu X., Eberhart R. (2002) Solving constrained nonlinear Tsochatzidis N.A., Maroulis K.E. (2007) Methods help remove
optimization problems with particle swarm optimization, black powder from gas pipelines, Oil Gas J. 105, 10, 52.
Proc. Conference on the Sixth World Multiconference on Wiak S., Krawczyk A., Dolezel I., (eds), (2008) Intelligent
Systemics, Cybernetics and Informatics, 14–18 July, Orlando, computer techniques in applied electromagnetics, Springer,
FL, USA , Vol. 5, pp. 203–206. Berlin, Heidelberg.
Kaipio J. (2008) Modeling of uncertainties in statistical inverse Wood N.B. (1981) A simple method for the calculation of
problems, J. Phys. Confer. Ser. 135, 1, 107–118. turbulent deposition to smooth and rough surfaces, J. Aerosol
Kennedy J., Eberhart R. (1995) Particle Swarm Optimization, Sci. 12, 3, 275–290.
Proc. on IEEE International Conference on Neural Networks. Wu X., Li C., Jia W., He Y. (2014) Optimal operation of trunk
IV, 27 November–1 December, 1995, Perth, WA, Australia, natural gas pipelines via an inertia-adaptive particle swarm
Australia, pp. 1942–1948. optimization algorithm, J. Nat. Gas Sci. Eng. 21, 10–18.
Khan T.S., Alshehhi M.S. (2015) Review of black powder in gas
pipelines–An industrial perspective, J. Nat. Gas Sci. Eng. 25,
66–76. Appendix A
Khan T.S., Alshehhi M.S., Stephen S., Khezzar L. (2015)
Characterization and preliminary root cause identification of Proof for the Proposition 3.1
black powder content in a gas transmission network – a case
study, J. Nat. Gas Sci. Eng. 27, 769–775. 1) Sufficiency: To prove the uniqueness of the optimal solu-
Kharoua N., Alshehhi M., Khezzar L. (2015) Prediction of Black tion, let us consider the general structure of tree-shaped gas
Powder distribution in junctions using the Discrete Phase pipeline network, as shown in Figure 3. The network con-
Model, Powder Technol. 286, 202–211. sists of N pipes, N + 1 nodes and (N  1)/2 junctions.
Kharoua N., Alshehhi M., Khezzar L., Filali A. (2017) CFD According to the 1-D model, the concentration at the
prediction of Black Powder particles’ deposition in vertical and end node of each branch pipe is given as a function of BP
horizontal gas pipelines, J. Petrol. Sci. Eng. 149, 822–833.
Liu C.A. (2008) New dynamic constrained optimization PSO
sources, where we have ðN 21Þ þ 1 unknown BP source to
algorithm, Natural Computation, Proc. Conference on the 4th be identified. The following equations can be easily found
ICNC’08, 18–20 October 2008, Jinan, China, Vol. 7, pp. 650– if we split the influence of each source on each measurement.
653. M 1 ¼ f ðS 1 Þ; ðA1:aÞ
Madoliat R., Khanmirza E., Pourfard A. (2017) Application of
PSO and cultural algorithms for transient analysis of natural M 2 ¼ f ðS 1 ; S 2 Þ ¼ f21 ðS 1 Þ þ f22 ðS 2 Þ; ðA1:bÞ
gas pipeline, J. Petrol. Sci. Eng. 149, 504–514.
Marini F., Walczak B. (2015) Particle Swarm Optimization
(PSO): A tutorial, Chemometr. Intell. Lab. Syst. 149, 153–165. M 3 ¼ f ðS 1 ; S 2 ; S 3 Þ ¼ f31 ðS 1 Þ þ f32 ðS 2 Þ þ f33 ðS 3 Þ;
Min Y., Jiayue Z., Damin Z. (2017) Immunization strategy .. ðA1:cÞ
based on discrete particle swarm optimization algorithm in .
BBV network, Proc. Conference IEEE Intelligent Systems and
Control, 5–6 January 2017, Coimbatore, India, pp. 208–211.
The compact form of above equations is
Parsopoulos K.E., Vrahatis M.N. (2002) Particle swarm opti- XðN 21Þþ1
mization method for constrained optimization problems, Mk ¼ i¼1
fki ðS i Þ: ðA2Þ
Intell. Technol. Theory Appl. New Trends Intell. Technol.
76, 1, 214–220. It is noted that the functions involved must be invert-
Saravanan R. (2006) Manufacturing optimization through intel- ible, which are determined by the physical properties of
ligent techniques, CRC Press, FL, USA. the gas transmission pipelines. Although there is not expli-
Sherik A.M., Zaidi S.R., Tuzan E.V. (2008) Black powder in gas cit mathematical function describing the mapping between
transmission systems, Proc. Conference on CORROSION, S and M, this function must be invertible at a particular
16–20 March, New Orleans, LA, pp. 16–20. time instant because the physical properties of pipeline do
Sherik A.M. (2008) Black powder: study examines sources, not change. Let us start with the first source S1 in equation
makeup in dry gas system, Oil Gas J. 106, 30, 54–59. (A1.a), the concentration could be determined by the first
Shiono N., Suzuki H. (2016) Optimal pipe-sizing problem of tree- measurement M1 only. Based on equation (A1.b), the
shaped gas distribution networks, Eur. J. Oper. Res. 252, 2,
second measurement M2 is represented as a function of
550–560.
(S1) and (S2). Then, f21 (S1) can be obtained according to
Tarantola A. (1987) Method for data fitting and model param-
eter estimation. Inversion Problem Theory, Elsevier Science,
S1, which has been calculated from equation (A1.a). There-
New York, NY. fore, f22 ðS 2 Þ, the remaining term of M2, could be used to
16 J. Shi et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 74, 47 (2019)

M2

S1 S2 S4 S5
S3 Sk
M3
M1

M4
Fig. 3. Generalized tree-shaped gas transmission pipeline network with additional BP source.

Pj1
determine the concentration at source S2. Consequently, where, i¼1 fjþ1;i ðS i Þ represents the influence of source S1
the remaining sources S k can be determined step by step. to Sj1 on Mj+1. fj+1,j(Sj) and fj+1,j+1(Sj+1) represent the
2) Necessity: If the number of client node measurements influence of Sj and Sj+1 on Mj+1 respectively. Therefore,
is less than the number of unknown BP source, for example, Sj is a preliminary condition to identify Sj+1. However,
M j (1 < j < ðN 21Þ þ 1) is not available, i.e. is not measured Sj cannot be identified because Mj is not available. In
other word, unique solution of Sj and Sj+1 cannot be
due to fault.
achieved in this situation.
Let us consider equation Mj+1:
Xj1    
M jþ1 ¼ i¼1
f jþ1; i ðS i Þ þ fjþ1; j S j þ fjþ1; jþ1 S jþ1 ; ðA3Þ