Application of Particle Swarm Optimization (PSO) algorithm for Black Powder (BP) source identi ﬁ cation in gas pipeline network based on 1-D model

. 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 corrosion of internal walls of the pipeline, which is a complex chemical reaction. This work aims to develop a novel algorithm for BP source identi ﬁ cation within gas pipelines network based on a 1-D model of BP transport and deposition. The optimization algorithm for BP source identi ﬁ cation is developed based on the well-known Particle 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 identi ﬁ ed and quanti ﬁ ed simultaneously. Extensive simulation studies are conducted to validate the effec-tivity of the optimization algorithm.


Introduction
With the wide use of natural gas, it is challenging to maintain operational efficiency and safety for gas transmission pipeline network, which is a complex system with pipeline length varying from hundreds to thousands of kilometers (Banda et al., 2006;Tobin and Shambaugh, 2006). Black Powder (BP), occurring in both liquid and gas pipelines, is the name given to the black particulates and sediment found in pipelines. It is mainly generated due to the chemical reaction of H 2 S, water and iron, resulting in a mixture of fine particle corrosion product and other solids, such as sands, clays, metal or construction debris, and liquid hydrocarbons chemically incorporated with any quantity of iron sulfide, iron carbonate and iron oxide contamination Sherik et al., 2008;Sherik, 2008;Trifilieff and Wines, 2009). Although there exist various compositions of BP, they possess some common characteristics, for example, adsorption, high specific gravity and difficult to clean. In terms of BP generation, they could come from the gas source in the gas field and also result from movement of the upstream point with gas flow, which originates from the corrosion of internal wall of any pipeline. Specifically, corrosion could only occur in the pipeline with the flow of "wet" gas (Beavers and Thompson, 2006;Baldwin, 1998), which is natural gas mixed with a certain amount of water, providing a necessary condition for corrosion. The presence of BP in gas pipeline network threatens the safety of gas industry and also has lots of undesirable influence on the operating companies all over the world (literature reviews of . First of all, the gas quality will be greatly affected by intolerable concentration of solid particles. Secondly, gas transmission cost will gradually increase, resulting from raised pressure drop due to changes of internal wall roughness and diameter over a long time of operation. Thirdly, the contaminations of BP could cause compressor failures, erosion of control valve and instrument clogging, etc. (Baldwin, 1998). Finally, the cleaning/removal of BP and frequent replacement of customers' cartridge filter elements could increase the expenses each year.
To reduce the influence of BP on gas industry, many companies attempted to manage and control the BP in gas transmission pipeline network. In general, the existing methods could be classified into two different aspects (Al-Qabandi et al., 2015;Cattanach et al., 2011;Trifilieff and Wines, 2009;Tsochatzidis and Maroulis, 2007): (1) BP removal and (2) prevention methods. The approach of removing BP from gas transmission pipeline is successful in several respects, including increasing operational safety, altering operational parameters, increasing operational efficiency and facilitating effective corrosion inspection. If BP problem is not serious (small quantities), mechanical cleaning is widely used for most of the pipelines.
Normally, mechanical scrapers are used to clean pipeline wall, such as mechanical pigs are deployed into a pipeline to scrape debris from pipeline wall and remove BP. Although mechanical scrap is efficient to keep the pipeline in fairly clean condition, frequent cleaning could damage the pipeline wall by exposing fresh steel surface under "wet" gas, resulting in excessive corrosion and BP generation. In addition, mechanical cleaning is not effective with complex BP formation, which calls for a combination of chemical cleaning. It has been practiced that chemical cleaning, acting as liquid soap, could significantly improve the efficiency of BP removal combined with mechanical cleaning, such as pigging operations, which acts as brush (Trifilieff and Wines, 2009). Separators and cyclo-filters are also widely installed to reduce the BP concentration, where these devices could physically knock out the BP particles in the gas stream. Then, the BP particles are collected at the bottom in a collection hub. However, this method only applies to gas stream with a high concentration of solid particles and also relatively large particle size (more than 10 microns). The removal approach is practiced as a good option to protect the downstream operations, but it has several disadvantages (Trifilieff and Wines, 2009). These methods are not single trial, but frequent application, which account for a large amount of expense. Also, these solutions cannot address the location of BP formation, that is to say, they are afterward remedy. In addition, chemical cleaning needs subsequent handling procedures, which could be costly and challenging if the chemical disposal is toxic or harmful to the environment. To compensate the above drawbacks of removal approach, gas operators have an alternative solution to prevent the occurring of BP generation, which is considered as a consequence of corrosions at the internal wall of gas transmission pipeline (Cattanach et al., 2011;Trifilieff and Wines, 2009). To prevent the occurring of corrosions, the inner wall of gas transmission pipeline is normally coated with high solids solvent-based epoxy polyamine films, which is used to protect the inner surface of pipelines. This approach is practiced to be cost effective; however, it is difficult to be applied to buried pipelines. Reducing water contamination is another approach to prevent the occurring of pipeline corrosions, which is based on the philosophy of internal corrosions, are largely related to the case of "wet" gas. Namely, this approach is moisture control.
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 problem of BP source until now. The main reason is the unknown information about the BP generation and where it could be generated. The developed dynamic models for BP transport and deposition are limited for gas transmission pipeline network, which in turn restricts the study of model-based methods for BP identification. It is well known that full three-dimension (3-D) Computational Fluid Dynamics (CFD) simulation software is specialized for multiphase flow modelling; however, it seems impossible to simulate the 3-D dynamics of large pipeline networks (e.g. 100 km), which is under studied in these works (Kharoua et al., 2015(Kharoua et al., , 2017. Filali et al. (2016) proposed a 1-D approach for modelling transport and deposition of BP particles in gas transmission network, where the authors discussed behaviors of particles with different diameters. Moreover, two different deposition models (Fan and Ahmadi, 1993;Wood, 1981) were compared to calculate the bed height, which was validated by Discrete Phase Model (DPM) based on CFD software.
The main contribution of this work is to develop a tool to identify the BP source in gas pipeline network, which is modeled as tree-shaped gas transmission network with BP dynamics of motion and distribution along each pipe. BP source identification is formulated as a constrained optimization problem, which is solved by the Particle Swarm Optimization (PSO) techniques (Delice et al., 2017;Eberhart and Shi, 1998;Kennedy and Eberhart, 1995;Min et al., 2017). PSO is a popular stochastic optimization technique with some features and advantages compared to other optimization algorithms such as Ant Colony Optimization (ATO) and Genetic Algorithm (GA) (Wiak et al., 2008;Saravanan, 2006). These features and advantages include (1) taking real numbers as particles; (2) few parameters need to be tuned and (3) simple implementation and effective global search capability, hence it could be a good candidate. In fact, PSO technique has been widely used in oil and gas industry, for example, Wu et al. (2014) optimized the operation of trunk natural gas pipelines via PSO based algorithm, and Madoliat et al. (2017) also successfully applied PSO to the transient analysis of natural gas pipeline. The PSO algorithm is used to solve a formulated optimization problem for similar applications. In particular, the basic PSO algorithm has been improved by incorporating the inertia-adaption technique to solve a constrained nonlinear optimization problem in Wu et al. (2014). Then, it was used to solve a more complex problem in Madoliat et al. (2017), where the solution of nonlinear PDE flow equations could be obtained simultaneously.
The paper is structured as follows. Section 2 will briefly discuss the 1-D approach for modelling BP particles in gas networks, which is detailed in Filali et al. (2016). A general structure of tree-shaped model of gas transmission pipeline network is given. In Section 3, PSO-based optimization algorithm for BP identification is presented. In Section 4, extensive simulation results and discussions of the optimization algorithm for BP source identification are presented. Finally, conclusions are given in Section 5.

Modelling for dynamics of BP particles in gas transmission pipeline
In this section, a schematic of tree-shaped gas transmission pipeline network will be built based on a set of pipeline connection rules. In addition, the methodology of a 1-D approach for modelling the dynamics of transport and deposition of BP particles in the gas pipeline network will be explained (Filali et al., 2016).

A tree-shaped model of gas network
In gas industry, the natural gas pipeline network is a highly integrated transmission and distribution grid that could transport natural gas from its origin to any position of high gas consumption demand. In many cases, the natural gas produced from a well has to be transported on a very long distance to the point of use. In order to maintain effective distribution of natural gas, the gas transmission system is extensive, consisting of complex pipeline topology. However, the entire gas transmission network could be separated into several sub-network with single gas well. This simplification also contributes to the investigation of BP identification, which will be included in following sections of this paper. In this case, the gas transmission pipeline network is properly generalized as tree-shaped model (Babonneau et al., 2012;Shiono and Suzuki, 2016), which could be built based on following rules for pipe connection: 1. Gas transmission pipeline network without loop; 2. Junction connected with three pipes; 3. One-way flow of gas stream.
Based on these connection rules, a schematic of treeshaped gas pipeline transmission network is given in Figure 1. There is only one source for gas supply, and each junction has three pipe connection, including two main pipes and one branch pipe. The arrow on the pipe shows the flow direction of gas stream. The customer is located at the end of each branch pipe. The sources, main pipes, branch pipes and customers are denoted as S, P, B, and M respectively.

One-D approach of modelling BP particle dynamics
In this section, the 1-D approach for modelling dynamics of BP particles in gas pipeline (Filali et al., 2016) is described. It is well known that the CFD software is popular to simulate 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 kilometers. In this case, Filali et al. (2016) proposed a simplified 1-D approach of modelling, where the flow of gas is continuously mixed with BP particles, and the behavior of particle movement is modelled based on the dusty gas assumption and the usage of analytical solutions of steady scalar advection/reaction equation.
The governing equation for gas stream mixed with solid particles is given by the following one-dimensional advection-diffusion-reaction equation: where C is the particle concentration in the gas flow, U is the average velocity of gas flow inside the pipeline, D diff is the diffusion coefficient, S is a term to describe the deposition, pickup or generation of solid particles. For fully developed turbulent flow, the axial diffusion term D diff @ 2 C @x 2 is negligible compared to the dominant advection term U @C @x and the following expression is valid, For steady state solution, equation (1) could be simplified as: where, b dep is the deposition rate and it could be calculated based on two methods, for example, Wood (1981), Fan and Ahmadi (1993), and n gen is the generation rate of the BP inside the pipeline and it is considered as an unknown parameter.

PSO-based optimization algorithm for BP identification
In this section, the objective of BP identification is presented as a constrained optimization problem, which will be solved by a PSO-based optimization algorithm. The PSO algorithm is a population-based evolutionary search algorithm inspired by social behavior of animals such as bird flocking. Basically, each particle in the swarm has a position and velocity, with its position representing a candidate solution in the multi-dimensional solution space and velocity indicating moves from one position to another. A fitness function is defined to evaluate each particle until certain convergence criteria is satisfied. During the searching process, the particle with fitness value will be selected as local/global best particle. The updated equations of position and velocity are given in the literature (Marini and Walczak, 2015) as follows: where, x ij (t) denotes the position of particle i in dimension j at time t, with j = 1, 2, 3, . . ., n x , and n x is the dimension of the solution space. The updating of position is calculated by adding the velocity v ij (t) to current position. y ij (t) is the personal best solution in dimension j for particle i, whileŷ j t ð Þ represents the best solution found so far. c 1 and c 2 are positive acceleration factors used to scale the contribution of cognitive and social component. r 1j ; r 2j 2 0; 1 ½ are random variables with normal distribution.
The structure of this section is outlined as follows: firstly, the objective of BP identification is formulated, which is developed based on the 1-D model of BP particles transport and deposition. Secondly, the BP source could be identified and quantified by sole running of the PSO-based optimization algorithm. Finally, a proof will be provided to support the uniqueness of the formulated optimization problem.

Formulation of BP identification problem
It is noted that additional unknown sources of BP generation are considered in this study. Along with that, the objective of this paper is to identify the BP sources and quantify the generated BP concentration using several measurements of BP concentration at the end node of each branch pipe, which is a point of natural gas use (client). As most of the gas transmission pipeline network is buried underground, therefore the sensor used for BP concentration measurement is only installed at each client (Abou-Khousa et al., 2015). The presence of an optimization algorithm for BP source identification is developed based on a 1-D model, which has several assumptions (Filali et al., 2016), such as, one-way flow of gas stream, constant velocity of gas flow at each pipe, average mixture of solid particles and gas stream, fixed size/diameter of solid particle, etc., but some additional assumptions should be also taken into consideration. The accuracy of BP source identification is largely restricted because of limited measurement of BP concentration along the pipe, which mainly results from two aspects: (1) lack of reliable measuring device that can be installed underground; (2) the "additional" BP generation from the internal walls is unpredictable, that is, the occurrence of chemical corrosion is highly random and the location is difficult to be identified. Therefore, the gas supply station is assumed in this paper as the main BP generation source, while the junctions stand for additional source, which essentially represent the amount of BP generated along its downstream. Herein, two assumptions are given for the problem formulation of BP identification: (a) BP source is located in junction only. (b) Additional BP particle generation occurs in the main pipe only.
According to Figure 1, there is one main BP particle source and several unknown additional sources located at each junction, which are the connection points for every two main pipes and one branch pipe. Let us consider the gas transmission pipeline network as the system. Measurements of BP concentration at each client are the system inputs and estimations of BP concentration at each source are the system outputs. Therefore, the following constrained optimization problem (Aguirre et al., 2007;Cagnina et al., 2008;Hu and Eberhart, 2002;Liu, 2008;Parsopoulos and Vrahatis, 2002) can be formulated with weighted sum of errors (e.g. absolute error) as its cost function. s:t: where, S 1 , S 2 , S 3 , . . ., S i represent the estimation of BP concentration at sources (S 1 is the main source, and S 2 , S 3 , . . ., S i are the additional source). m is the number of clients and BP measurement points. M k is the measurement of BP concentration at the client k, which is measured by BP concentration sensor located at the end of each branch pipe and it is simulated using the 1-D model presented in the previous section.M k is the estimated value of BP concentration, k k is the weighting parameters for each system input, which could be selected as inverse of each measurement of BP concentration at corresponding client. The optimization algorithm for BP source identification is presented as a flowchart in Figure 2.

Simulation studies
In this section, the optimization algorithm for BP source identification is applied to a real gas transmission pipeline network, including 15 pipes, 16 nodes, and 8 junctions.
The geometric parameters of this gas pipeline network are practical data from an existing network. The topology of this network is given in Figure 1, as well as its geometric parameters in are divided into two aspects, including investigations on disturbed measurements and model parameter mismatch. Finally, some discussions are presented.

Ideal situations
In this section, the optimization algorithm for BP identification is validated on the gas network given in Table 1. It is noted that, initially, the measurements of BP concentration at the end node of each branch pipe (client) are simulated based a sequence of BP sources given as follows. By sole running of the 1-D model, the measurements of BP concentration at each client could be generated according to the BP concentration at each source given in Table 2. Then these measurements shown in Table 3, will be saved and assumed as known parameters, i.e. BP concentration measurements. This is the first step for following simulation studies. The second step is the estimation of BP concentration at each source using the measurements of BP concentration. It is noticeable that the aim of BP optimization problem is the reverse process, where the BP concentration at each source will be estimated by application of the optimization algorithm.
Simulation results of identification for BP concentration at sources are given in Figure 4. In this simulation study, the parameters of PSO are selected as: k k ¼ 1 M k , n = 300, m = 50, c 1 = 2, c 2 = 2, where, k k is the weighting parameter in equation (6), n is the quantity of particles, m is the maximum number of iteration (termination condition), and c 1 , c 2 are the acceleration constants. It is noticeable that the PSO parameters are tuned by trial and error, which can ensure satisfactory performance with acceptable computational load. PSO is a stochastic optimization technique, therefore the algorithm is repeated three times (blue, red  9.5 6 7.5 5 4.5 8 5.5 0.5 start Set the BP Source at each junction, as variables to be optimized Initialization of the position and velocity of each particle Evaluation of the fitness of each particle For each particle, set: Local best fitness = current fitness Local best position = current position Global best fitness = min(local best fitness) Update velocity and position for each particle Evaluation of the fitness of each particle   Table 3. In addition, the statistic results of each test are summarized in Table 4. The estimation value and absolute error are given for each test and source, and an overall average with the Standard Derivation (SD) are given for all sources. It could be concluded that the optimization algorithm for BP source identification is satisfactory and able to achieve small estimation errors with perfect measurements of BP concentration at clients.

Sensitivity studies
As stated in Section 4.1, the BP concentration at sources can be well estimated/reproduced by the PSO algorithm when the model is completely known and there are no uncertainties in measurements. However, the efficacy of this approach is still questionable in practice, because all the models are at best only approximations for reality due to various model errors or parametric uncertainties. This problem can be formulated as follows:  with h the vector of geometric and flow conditions parameters, the measurement error, M k the real measurements, b M k the 1-D model approximation and S* the vector of BP concentration at sources which are unknown. The aim is to effectively recover S* via PSO algorithm relying on M k , regardless of the presence of uncertainties in h and . A good approach to solve this problem is associated to robust model calibration and inverse problem theory (Tarantola, 1987;Kaipio, 2008). As a preliminary study, three case studies are simulated in this section. Firstly, the measurement error on M 3 and M 5 are simulated. Secondly, the disturbances in flow measurements are considered. Thirdly, parametric uncertainties 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.

Black Powder measurement uncertainty
In this section, the optimization algorithm for BP source identification is tested under the situations of measurement uncertainty, where M 3 and M 5 are increased by 5% respectively. Basically, e is regarded as an additive nonzero noise and no variations in h in this scenario. The tuned parameters for the algorithm are the same with ideal situations. Some figures and statistic data are presented as preliminary results. It is noted that the problem formulation in equation (7) is a deterministic case, which will be generalized in the next stage of experimental validations and tackled by the likelihood approaches (Tarantola, 1987), where measurement errors can be modelled with a probabilistic distribution and confidence intervals will be computed for estimated BP concentration at sources within Bayesian framework.

Uncertainty in M 3 measurement
In this case, the BP concentration measurement M 3 is increased by 5%. The performance of identification is shown in Figure 5, and also the corresponding statistic data is presented in Table 5. It can be seen from Figures 5c and 5d, 10% and 5% estimation error for sources 3 and 4 (e) (f )       respectively, which are expected and acceptable. Although the inaccurate measurement occurring at M 3 affects the BP concentration estimation at sources 3 and 4, and it has negligible influence on the BP concentration estimation at sources 5, 6, and the downward sources.

Uncertainty in M 5 measurement
In this case study, M 5 is increased by 5% and the performance of the optimization algorithm is presented in Table 6. It is noted that the estimation errors of sources 5 and 6 are 15% and 10% respectively, which are high relatively compared with their downward sources. The estimation error of sources 7 and 8 could be negligible.

Mass flow rate measurement uncertainty
Similarly, variations in h are simulated in this section, where the mass flow rate in pipe 3 (Q 3 ) and pipe 5 (Q 5 ) are increased by 10% and 5% respectively. The tuned parameters for PSO are same with the ideal situations. One objective of this paper is to study the robustness of the model fidelity to measurement data given parametric uncertainties. Although these simulation tests are not representative in practice where noise/uncertainties always exist in a stochastic way, they can provide some preliminary results at this stage.

Uncertainty in Q 3 measurement
In this case study, the mass flow rate measurement of Q 3 is increased by 10%. The simulation results are shown in Figure 6, along with its statistic data in Table 7. It can be seen in Figure 6, the variation of flow rate Q 3 has a very small influence on the optimization, which proves the robustness of the optimization algorithm under uncertainty.

Uncertainty in Q 5 measurement
In this case study, the mass flow rate measurement of Q 5 is increased by 5%. The simulation results are shown in Table 8. The estimation error is less than 2% for each source, which is acceptance and negligible for this work.

Deposition rate uncertainty
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 concentration for each client is always higher than the expected value, and the estimation error is shown in Table 9. This is reasonable because the increase of roughness will result in more BP deposition, in other words, more BP should be generated in order to achieve the same measured BP concentration as expected. However, the average estimation error is nearly 2%, which is still acceptable in the practical situation.

Conclusion
In this paper, the problem of BP source identification is studied by applying a PSO-based optimization algorithm, which is developed upon a 1-D model of BP transport and deposition. The 1-D model is a simplified approach for modelling the dynamics of BP particles in gas transmission pipeline network. A schematic of tree-shaped gas network is proposed, which is generalized with a set of connection rules. Through some preliminary simulation studies, the PSO-based algorithm is validated to be a useful technique for BP source identification. At this stage, the model errors and uncertainties are assumed deterministic though they are not representative in practice. The main purpose of this paper is to discuss the application of PSO techniques on BP source identification. More work will be done together with the experimental studies in the next stage, including (1) analysis of likelihood approach which has specified probabilistic distribution to model error; (2) improvement of the 1-D model by considering pickup case and various particle sizes. determine the concentration at source S 2 . Consequently, the remaining sources S k can be determined step by step.
2) Necessity: If the number of client node measurements is less than the number of unknown BP source, for example, þ 1) is not available, i.e. is not measured due to fault.
Let us consider equation M j+1 : where, P jÀ1 i¼1 f jþ1;i S i ð Þ represents the influence of source S 1 to S jÀ1 on M j+1 . f j+1,j (S j ) and f j+1,j+1 (S j+1 ) represent the influence of S j and S j+1 on M j+1 respectively. Therefore, S j is a preliminary condition to identify S j+1 . However, S j cannot be identified because M j is not available. In other word, unique solution of S j and S j+1 cannot be achieved in this situation.  Fig. 3. Generalized tree-shaped gas transmission pipeline network with additional BP source.