Fast source term estimation using the PGA-NM hybrid method

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

Contents lists available at ScienceDirect

Engineering Applications of Artificial Intelligence journal homepage: www.elsevier.com/locate/engappai

Fast source term estimation using the PGA-NM hybrid method a,⁎

MARK

b

Hui Li , Jianwen Zhang a b

College of Information Science and Technology, Beijing University of Chemical Technology, Beijing 100029, China Fluid Mechanics and Heat Transfer Laboratory, Beijing University of Chemical Technology, Beijing 100029, China

A R T I C L E I N F O

A BS T RAC T

Keywords: Gaussian dispersion model Pasquill-Gifford dispersion model Nelder-Mead algorithm Parallel genetic algorithm Inversion problem

There are significant challenges related to estimating the source term of the atmospheric release. Urged on by robots in performing emergency responding tasks, a fast and accurate algorithm for this inversion problem is indispensable. Sometimes the NM simplex algorithm is efficient in the optimization problem, but sometimes the quality of convergence is unacceptable as a numerical breakdown, even for smooth and well-behaved functions. In contrast, full convergence might be seen in parallel genetic algorithms with a comparative slower convergence. In this paper we combine the PGA and the NM simplex algorithm by initializing simplex from the final individual of PGA results and obtaining the best vertex through simplex algorithm thereafter. A numerical simulation of the proposed algorithm shows noteworthy improvement of efficiency and robustness, compared with the PGA or the NM algorithm only.

1. Introduction Accidental releases may result in fires, spills or explosions that involve hazardous materials, such as chemicals or radionuclides. The atmospheric dispersion model, together with the source term, are used to determine the consequences of accidental releases of hazardous or toxic materials. Computational fluid dynamic (CFD) models simulate the interaction of liquids, gases or solid particles with surfaces defined by boundary conditions. The CFD model provides complex analysis of single-phase fluid flow based on conservation of mass and momentum by resolving the Navier-Stokes equation using finite difference and finite volume methods in three dimensions (Bonifacio et al., 2015). CFD models are the most accurate model for situations in complex scenarios such as transonic or turbulent flows (Stern and Yamartino, 2001). Lagrangian models define a region of air as a box containing an initial concentration of pollutants, then mathematically follow the trajectory of the parcels as it moves downwind (Chow et al., 2008). Lagrangian models allow the model to track the wind fields and temporal changes as it moves along its trajectory (Oettl et al., 2001). Eulerian dispersion models are similar to Lagrangian models in that they also track the movement of a large number of plume parcels moving from their initial location. The difference between Lagrangian models and Eulerian models is that Eulerian models use a fixed 3D Cartesian grid as a frame of reference rather than a moving frame of reference (Guha, 1997). The technical literature on Gaussian dispersion model is quite ⁎

extensive and dates back to the 1930s or earlier. It is perhaps the oldest and perhaps the most commonly used model type. Sir Graham Sutton derived an air pollutant plume dispersion equation in 1947 which included the assumption of Gaussian distributions for the vertical and crosswind dispersion of the plume and also included the effect of ground reflection of the plume (Sutton, 1947). The width of the plume is determined by σy and σz, which are defined either by stability classes (Pasquill, 1961 and Gifford, 1976) or distance from the source. One severe limitation of the Gaussian model is that the model is not designed to model the dispersion at sites too close to the source. Also Gaussian models have been shown to consistently overpredict concentrations in low wind conditions (Sokhi et al., 1998). A rapid and accurate estimation of the source term allows faster and more efficient action for emergency response teams, in addition to estimating the impacted areas and ambient concentrations. The proper emergency action depends partly on knowing the source locations, the source emission rate, and the likely future dispersion routes of the releases. The methods for estimating source term parameters usually are divided into two categories: probability or statistics methods and optimization methods. The former methods use atmospheric dispersion models and some statistical approaches to establish mathematical models for quantitative estimation. The feature of this category is that it demands the priori distribution of parameters and requires a large amount of observed data for statistics (Schmidpeter et al., 2015). The inverse problem can also be formulated as a non-linear optimization problem, whose cost function is given by the least-square

Corresponding author. E-mail address: [email protected] (H. Li).

http://dx.doi.org/10.1016/j.engappai.2017.03.010 Received 16 August 2016; Received in revised form 20 February 2017; Accepted 31 March 2017 Available online 11 April 2017 0952-1976/ © 2017 Elsevier Ltd. All rights reserved.

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

H. Li, J. Zhang

2014). The hybrid algorithm has been effectively applied in the optimization problem of the Muskingum parameter estimation (Haddad et al., 2015). Also Ma has mentioned using the GA-NM hybrid algorithm to identify gas emission source in Ma et al. (2013) and compared it with other eight algorithms without going into great detail. However, to the best of our knowledge, the combination of PGA and NM has not emerged yet. In this paper, we attempted to aid a mobile robot in finding the gas emission location. We implemented the PGA-NM hybrid algorithm to estimate the location and the source emission rate as fast and accurately as possible. Based on these results, the mobile robot carried out several successive tasks including applying other new techniques such as the terahertz time-domain spectroscopy (THz-TDS) to locate the leakage source more accurately, and rehearsing preset emergency procedures thereafter. The numerical simulation validated the capability of the proposed algorithm to achieve satisfying goal of converging fast and being tolerant to measurement errors. The paper is organized as follows. The Gaussian dispersion model and its inversion problem are described in Section 2. The GA, the PGA, the NM simplex algorithm, the hybrid algorithm and the solutions for the source term estimation problem are presented in Section 3. Section 4 provides a numerical simulation and discuss the results of the hybrid algorithm in the ideal condition and the condition with measurement errors. Section 5 concludes the paper finally.

difference between the measured and simulated concentrations by the mathematical model (Haupt, 2005). CFD models, Lagrangian models and Eulerian models are not suitable for fast source term estimation because they would consume a lot more time than the Gaussian models to finish a forward simulation process. Therefore, Gaussian models are often used for predicting the dispersion of continuous plume releases and solving the optimization problem. The approaches to solve inverse problem based on optimization theory could be divided into three categories: indirect optimization, direct optimization and intelligent optimization. The indirect optimization computes the first and second derivative of the cost function to determine the direction of reducing the cost function, thus the cost function is minimized gradually to hit an optimum. Typical methods of this category are the Newton method (Venkatakrishnan, 2015), the gradient-type search (Yu, 2015) etc. The direct optimization does not need any derivative of the cost function, but rely on the initial value. The Nelder-Mead(NM) simplex method, proposed by John Nelder and Roger Mead in 1965 (Nelder and Mead, 1965), is a representative method of this category. The NM method is suitable for the multidimensional search space, and it has been successfully applied to the identification of multiple-point releases in paper (Baghmisheh et al., 2012). The main drawback of this method is that the search may be trapped in a local optimum, even for smooth and well-behaved functions. The intelligent optimization method applies artificial intelligence to repeatedly tune calculation model or parameters and compare cost function values, ultimately achieves an optimum. Genetic algorithms (GA), particle swarm optimization (PSO) and simulated annealing (SA) are all heuristic intelligent optimization algorithms. Inspired by natural evolution, such as selection, mutation and crossover, GAs are very often used for large sized optimization problems including the inversion problem as mentioned previously (Haupt et al., 2009). In many problems, GAs may have a tendency to converge towards the global optimum rather than some local optimum of the problem, so it is more suitable for global optimization problems (Singh et al., 2013). But slow convergence is a drawback of GAs. To make GA faster, one of the most promising choices is to use parallel implementations. The basic idea in parallel genetic algorithm (PGA) is to divide a task into chunks and perform the chunks (divideand-conquer) simultaneously using multiple processors. The first method to parallelize GAs is the global PGA or master-slave PGA. There is only one population in the global PGA as in the serial GA, and the evaluation of individuals are parallelized. Abramson and Abela (1991) has implemented a global GA on a shared memory computer to search for efficient timetables for schools. Fine grained PGAs and Coarse-grained PGAs are more sophisticated PGAs. Fine grained PGAs partition the population into single individuals. Selection and crossover are restricted to a small neighborhood, but neighborhoods overlap permitting some interaction among all individuals (Manderick and Spiessens, 1989). Coarse-grained PGAs use a few relatively large subpopulations (deme) and exchange individuals occasionally. The migration and communication topologies are two important factors in the performance of the Coarse-grained PGAs (Pettey et al., 1995). Although PGAs have better performance on the computation time than serial GAs, the improvement is not huge because the time required to send individuals back and forth will offset any performance gains. A good strategy for avoiding slow convergence toward the optimum is to combine NM simplex algorithm with PGA. The advantage of this hybrid algorithm is that it neither requires information of priori knowledge, nor a large number of observational data, nor the gradient of the cost function. Moreover, the NM method is protected from being trapped in a local suboptimum because of global optimization provided by PGA. The idea of combining GA with NM has sparsely emerged in some papers with little detail, such as (Haddad et al., 2015) and (Garca et al.,

2. Problem formulation 2.1. Point source dispersion model Airborne particulate dispersion can be described mathematically using the advection-diffusion equation: (1)

Ct + u ·(∇C ) − ∇·(K ∇C ) = Q

which is a partial differential equation (Hosseini and Stockie, 2016). However, a direct application of such equation is often too expensive. Methods based on the Gaussian approximation are widespread in commercial software used for regulatory and monitoring purposes (Agency,). For demonstration purposes we choose a basic Gaussian plume model for a continuous point source, which could easily be substituted with a more refined model. Suppose the source is located at the coordinates (x0, y0, z0), the wind blows along the x axis direction, and the wind speed is constant. Based on the Pasquill-Gifford dispersion model (Crowl,), Gaussian plume dispersion could be formulated as:

C (x , y , z ) =

⎡ (z − z 0 ) 2 (z + z 0 ) 2 ⎤ − 1 y − y0 2 − Q0 e− 2 ( σy ) ⎢e 2σz2 + e 2σz2 ⎥ , ⎢ ⎥ 2πσy σz u ⎣ ⎦

(2)

g / m3 ,

where C (x, y, z ) is the concentration of emissions, in at any sensor located at (x, y, z ); x, in meters, denotes the distance downwind from the emission source point; y, in meters, denotes the distance crosswind from the emission plume centerline; z, in meters, denotes the height above ground level; Q0 is the source emission rate, in g/s; u denotes the horizontal wind velocity along the plume centerline, in m/s; σy, σz, in m, are horizontal and vertical standard deviation of the emission distribution respectively. The stability conditions have a strong influence on dispersion. The dispersion coefficients σy and σz are functions of atmospheric conditions and the distance downwind from the release. The atmospheric conditions are classified according to six different stability classes, shown in Table 1 for rural conditions and Table 2 for urban conditions (Crowl). The graphic representation of the Pasquill-Gifford dispersion model is depicted in Fig. 1. Detailed parameters will be elaborated in the numerical simulation of Section 4.1. 69

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

H. Li, J. Zhang

At each iteration, the transformation is determined by computing one or more test points, together with their cost function values, and by comparison of these cost function values with those at these vertices. This process is terminated when the working simplex S becomes sufficiently small in some sense. The general algorithm is given in Algorithm 1.

Table 1 Recommended Equations for Pasquill-Gifford Dispersion Coefficients for Plume Dispersion (Rural conditions) (Crowl,). PG stability

σ y (m )

A

0.22x (1 + 0.0001x )−0.5

0.2x

B

0.16x (1 + 0.0001x )−0.5

0.12x

C

0.11x (1 + 0.0001x )−0.5

0.08x (1 + 0.0002x )−0.5

D

0.08x (1 +

0.0001x )−0.5

0.0015x )−0.5

E

0.06x (1 + 0.0001x )−0.5

F

0.0001x )−0.5

0.04x (1 +

σz (m )

0.06x (1 +

Algorithm 1. the NM simplex algorithm Require: xin Construct the initial working simplex S from xin. while the termination criteria is not satisfied do transform the working simplex. calculate the termination test values. end while return the best vertex of the current simplex S and its cost function value.

0.03x (1 + 0.0003x )−1

0.016x (1 +

0.0003x )−1

Table 2 Recommended Equations for Pasquill-Gifford Dispersion Coefficients for Plume Dispersion (Urban conditions) (Crowl,). PG stability

σ y (m )

σz (m )

A

0.32x (1 + 0.0004x )−0.5

0.24x (1 + 0.0001x )0.5

B

0.32x (1 + 0.0004x )−0.5

C

0.22x (1 + 0.0004x )−0.5

0.24x (1 + 0.0001x )0.5 0.2x

D

0.16x (1 + 0.0004x )−0.5

0.14x (1 + 0.0003x )−0.5

E

0.11x (1 + 0.0004x )−0.5

0.08x (1 + 0.0015x )−0.5

F

0.11x (1 + 0.0004x )−0.5

0.08x (1 + 0.0015x )−0.5

To construct the initial working simplex, we choose x 0 = xin , where xin is composed of the predetermined source term parameters to be estimated, such as Q0, x0,y0 and z0 in Eq. (2). The remaining n vertices are then generated to obtain a right-angled shapes of S:

xj = x 0 + hj ej ,

j = 1, …, n

(5) n

where hj is the stepsize in the direction of unit vector ej in R . Each step in the loop of Algorithm 1 consists of the following three processes.

2.2. Reverse problem formulation The reconstruction of the source term using the observed data from a set of sensors is an important and challenging inverse problem. Assuming that there are n sensors in the downwind direction of source emission to measure the concentration of emissions. Suppose the concentration calculated by the dispersion model in the position of sensor i is Cmodeli, and Cmeai is the corresponding concentration measured at sensor i. Then the cost function f is defined as sum of squares of the error between the measured concentration and calculated concentration, namely i i min fc = Σin=1 (Cmea − Cmodel )2

(1) Ordering Determine the indexes h , s, l of the worst, the second worst and the best vertex in the current working simplex S

fh = max f j , j

j≠h

fl = min f j j≠h

(6)

(2) Centroid Calculate the centroid c of the best side. The centroid c is the point opposite to the worst vertex xh

(3)

c≔

For the n concentration measurements obtained from sensors, there are several possible source terms that would fit the observations. The parameters mainly include Q0, x0,y0 and z0. Cmodel is derived from Eq. (2), substitute it for Eq. (3), we obtain:

min fc (Q0 , x 0 , y0 , z 0 ) ⎧ ⎡ (z − z 0 ) 2 (z + z 0 ) 2 ⎤ ⎫ − ⎪ i ⎪ 1 y − y0 2 − Q0 = Σin=1 ⎨Cmea − e− 2 ( σy ) ⎢e 2σz2 + e 2σz2 ⎥ ⎬. ⎢ ⎥ 2 πσ σ u ⎪ y z ⎣ ⎦⎪ ⎩ ⎭

fs = max f j ,

1 n

∑ xj . (7)

j≠h

(3) Transformation Compute the new working simplex from the current one. We replace only the worst vertex xh with a better point by using reflection, expansion or contraction with respect to the best side. All test points lie on the line defined by xh and c, and at most two of them are computed in one iteration. If it succeeds, the accepted point becomes the new vertex of the working simplex. If it fails, shrink the simplex towards the best vertex xl. In this case, n new vertices are computed. Simplex transformations in the Nelder-Mead method are controlled by four parameters α for reflection, β for contraction, γ for expansion and δ for shrinkage. We deliberately set the parameters following the constraints below:

(4)

Thus the inverse problem of atmospheric dispersion can be converted into the optimization problem above. 3. Algorithm 3.1. Nelder mead simplex method The Nelder-Mead simplex algorithm is one of the most popular algorithms for direct multidimensional optimization. It is widely used to solve parameter estimation problems where the function values are subject to noise. We apply the NM method to the source term estimation problem partly because measurements cannot be noise free and be accurately described by the dispersion model in real condition. The NM method maintains at each step a nondegenerate simplex A simplex S in Rn that is the convex hull of n + 1 vertices x 0, …, xn ∈ R n . From an initial simplex obtained in the previous phase, the NM simplex method performs a sequence of transformations of the working simplex, aiming at decreasing the cost function values at the vertices.

α > 0,

0 < β < 1,

γ > 1,

γ > α,

0 < δ < 1.

(8)

The following algorithm (refer to Wright (1997)) describes the working simplex transformations in step (3).

• • 70

Reflect Compute the reflection point xr ≔c + α (c − xh ) and fr ≔f (xr ). If fl ≤ fr < fs , accept xr and terminate the iteration (see Fig. 2(a)). Expand If fr < fl , compute the expansion point xe = c + γ (xr − c ) and fe = f (xe ). If fe < fr , accept xe and terminate the iteration. Otherwise

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

H. Li, J. Zhang

Fig. 1. The graphic representation of the Pasquill-Gifford dispersion model.

•

(if fe ≥ fr ), accept xr and terminate the iteration (see Fig. 2(b)). Contract If fr ≥ fs , compute the contraction point xt by using the better of the two point xh and xr. Outside Contraction: If fs ≤ fr < fh , compute xt = c + β (xr − c ) and ft = f (xt ). If ft ≤ fr , accept xt and terminate the iteration (see Fig. 2(c)). Inside Contraction: If fr ≥ fh , compute xt = c + β (xh − c ) and

•

ft = f (xt ). If ft < fh , accept xt and terminate the iteration. Otherwise, perform a shrink transformation (see Fig. 2(d)). Shrink Compute n new vertexes xj = xl + δ (xj − xl ) and f j = f (xj ), for j = 0, …, n , with j ≠ l (see Fig. 2(e)).

The shrink transformation is requisite because the introduction of it can prevent the algorithm from failing in some special case.

71

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

H. Li, J. Zhang

Fig. 2. Transformations in the NM simplex algorithm.

To ensure termination in a finite amount of time, we establish the termination criterion as following:

1 +1 i i Σin=1 {Cmea − Cmodel (x[k ] )}2 < ϵ, n+1

Selection. Crossover. Mutation. Evaluation (Calculation of the fitness). end while return the best individual and its cost function.

(9)

where x[k ] is in the working simplex, and ϵ is the predetermined precision.

Initial population is generated randomly, allowing the entire range of possible solutions. Selection operators give preference to better individuals, allowing them to pass on their genes to the next generation of the algorithm. We employ roulette wheel selection (Lipowski and Lipowska, 2011) for selecting potentially better solutions for recombination. The fitness function assigns a fitness to each individual, then this fitness level is used to associate a probability of selection with each individual. If fi is the fitness of individual i in the population, its probability of being selected should be

3.2. Genetic algorithm The genetic algorithm is an intelligent search technique that mimics the process of natural selection. GA generates solutions to optimization problems using heuristic inspired by natural evolution, such as selection, crossover and mutation. The bit representation traditionally used in genetic algorithms has some drawbacks when applied to this multidimensional, high-precision inversion problem. Genetic representation, selection, crossover and mutation should be elaborately designed to reflect the difference. In this paper, we represent the solution as a set of floating point parameters including Q0, x0,y0 and z0 in Eq. (2), by referring to Michalewicz (1994). Although these parameters are in different metrics, that does not matter because we regard them as a whole vector (Q0, x0,y0,z0) and force each parameter to be within the desired range. The fitness function has a direct connection to the cost function fc in Eq. (4). And it should monotonously decrease over fc, so we set the fitness function as:

f=

1/ f ckf

Pi = fi / Σ jn=1 f j

(11)

where n is the number of individuals in the population. Next, we take every pairs of two parent individuals and produce children from them. Suppose two selected parents are (Q1,x1,y1,z1) and (Q2,x2,y2,z2), we employ arithmetic crossover operator (Mondal and Maiti, 2003) to linearly combine two parent chromosome vectors, according to the equations:

⎫ ⎧ l1, xl1, yl , zl1) = kc·(Q1, x1, y , z1) + (1−kc )·(Q2 , x2 , y , z2 ) ⎪ ⎪ (Q 1 1 2 ⎨ ⎬ ⎪ l ⎪ ⎩ (Q2 , xl2 , yl2 , zl2 ) = (1−kc )·(Q1, x1, y1, z1) + kc·(Q2 , x2 , y2 , z2 )⎭

(10)

where kf is a positive integer or a reciprocal of positive integer. We then perform the GA following a usual procedure as Algorithm 2.

(12)

l2 , xl2 , yl , zl2 ) are two children, and kc is a l1, xl1, yl , zl1) and (Q where Q 2 1 random real number in the range of {0, 1}. We employ a new mutation scheme prevent the genetic algorithm from premature convergence, which is similar to the scheme in Gaffney et al. (2009). For each individual to mutation, we update a random component of the individual vector:

Algorithm 2. The genetic algorithm Generate initial population of individuals. while the termination criteria are not satisfied do

vli = vi + d , 72

d = ξ·(vhi − vli )

(13)

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

H. Li, J. Zhang

point, the NM algorithm attempts to find the best local optimum in its neighborhood. Iterations continue until the termination criterion is satisfied finally. See Algorithm 4.

where vli and vhi is the lower and upper domain bounds of the variable vi , and ξ is a random real number in the range of {−1/k m, 1/k m}. The result of mutation is only accepted if it is within the search scope. The uniform mutation is an alternative mutation scheme used in this paper. The uniform mutation operator selects a random component vi and replace it with vli , where vli is a uniform random value between the upper and lower bounds for that variable. The algorithm terminates when either a predetermined number of iterations has been finished, or a satisfactory fitness level has been reached for the generation.

Algorithm 4. The basic PGA-NM hybrid algorithm Generate initial population of individuals for the PGA-NM algorithm. Call Algorithm 3 Construct the initial working simplex S from output of Algorithm 3. Call Algorithm 1 return the best vertex and its cost function value.

3.3. Parallel genetic algorithm To make GA faster or allow a more extensive coverage of the search space, the most promising choices is to use parallel implementations. The basic idea behind most parallel programs is to divide a task into chunks and to solve the chunks simultaneously using multiple processors. There are three main types of parallel GAs: global PGA, fine-grained PGA, and coarse-grained PGA. The global PGA does not affect the behavior of the algorithm, while the last two methods change the way the serial GA works (Cant-Paz, 1999). We use global PGA to solve the source estimation problem, because it can be a very efficient method when the evaluation needs considerable computations. We compare computation time of the main GA operators and individual evaluation. The comparison based on the dispersion model with parameters (in the later numeric example) is shown in Table 3. In our global PGA implementation, the populations are stored in shared memory and each processor read an individual assign to it and write the evaluation results back. The global PGA also has the advantage of preserving the search behavior of the serial GA, so we do not have to change the algorithm much (Cant-Paz, 1999). As shown in Algorithm 3, the only difference between the global PGA and the serial GA (Algorithm 2) is the implementation of the evaluation part.

To avoid getting stuck in the suboptimal solution, and avoid the numerical breakdown of the NM algorithm, we may select three or more individuals from PGA final population, then perform the successive NM algorithm successively or parallelly from these points. In practice, most of the results are toward a same point, so we could accept this point as the global optimum. If the results are totally different from each other, we have two choices: compute again or choose the best one from the present results. This algorithm is demonstrated in Algorithm 5 as another choice. Algorithm 5. The GA-NM hybrid algorithm (optional) Generate initial population of individuals for the PGA-NM algorithm. Call Algorithm 3 Choose N best or almost best individuals xi (i = 1, …, N ). for each i ∈ [1, N ] do Construct the initial working simplex S from xi. Call Algorithm 1 Collect the best point yi and its cost function value. end for if most of the points are the same and their cost function values is the lowest then Accept this point p. else Choose the point p with the lowest cost function value from yi (i = 1, …, N ). end if return p and its cost function value.

Algorithm 3. The parallel genetic algorithm Generate initial population of individuals. while the termination criteria are not satisfied do Selection. Crossover. Mutation. Parallel evaluation. end while return the best individual and its cost function.

4. Numerical results and discussion In this section, we tested our hybrid algorithm in a numerical simulation and compare it to the results applying PGA only or the results applying NM simplex algorithm only. Also we attempted to compare with the results in other literatures.

3.4. PGA-NM hybrid algorithm It is a rational idea to combine the PGA with the NM algorithm. The PGA can be used to search behavioral parameters of the NM algorithm, like (Przystaka and Katunin, 2016). But in this paper, we prefer an unsophisticated combination. We narrow the searching area by PGA and then hit the optimum by the NM simplex algorithm fast and accurately. In concrete terms, the best points obtained by PGA are used for creating initial simplex of the NM algorithm. Starting from this

4.1. Sensors placement and model profile As sketched in Fig. 3, the sensors were placed based on the

Table 3 The computation time of main GA operations. Operation

Average time (μs)

Selection Crossover Mutation Evaluation

0.2 1.01 1.5 5.22 Fig. 3. the dispersion scenarios and the sensors placement.

73

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

H. Li, J. Zhang

the input of the NM simplex algorithm. We not only compared the PGA-NM algorithm with the NM simplex algorithm only and the PGA only, but also compared it with the PSO algorithm and the PSO-NM algorithm. After that, we made comparisons between the results of the PGANM algorithm with measurement errors and without measurement errors. Finally, we tested the performance of our PGA-NM algorithm on the embedded system. The hardware and software configurations are ARM Cortex-A7 quad-core processor (a Broadcom BCM2836 SoC) 900 MHz, on-board RAM 1 GB and Raspbian OS.

dispersion scenarios. From the position of 250 m downwind, we placed three sensors every 100 m outward. The three sensors were placed every 10 m crosswind from the assumed centerline. To avoid interference from the ground, we placed the sensor at the height of 3 m. A very high placement may not be a good choice, because if some sensor malfunctions the sensor may not be easily replaced or serviced. Assuming the source emission rate Q0 is 1500 g/s, and the point source is located at {x 0 , y0} = {8(m ), 10(m )} with a height of z 0 = 5. Suppose the atmospheric stability classes in the Pasquill-Gifford dispersion model is F, and the average wind speed is 2 m/s. The model with these parameters has already been shown before in Fig. 1. We placed an embedded system in the position of each sensor to calculate the gas concentration for each sensor from the PasquillGifford model using the GPS data of the source and the sensor. The workstation remotely collected all simulated concentrations from these embedded systems via the wireless channel. After computing the source position and emission rate using our algorithm, the workstation sent these results to the mobile robot. The robot then went to the source position automatically and verified the results. We did not pay much attention on the number optimization of sensors. Nine is big enough than the number of variables to be backcalculated, and small enough to avoid wasting computational resources. The observed data were generated from the previous dispersion scenarios and the Pasquill-Gifford model, as illustrated in Table 4.

4.3. NM simplex algorithm only As shown in Fig. 4, from the initial point {0, 0, 0, 0}, the NM simplex algorithm converges to the optimum soon. The convergence of Q0 is illustrated in Fig. 5. The right optimum has been obtained because the initial point is not very far away from the optimum. We changed the initial point and obtained the results in Table 5. In some case (row 2), the algorithm failed unfortunately by getting stuck in a certain state. With the initial point {15000, −150, 150, 10}, the algorithm converges to a sub-optimum. The final Q0 1494.83 is toward 1500, but the position is not well estimated. Additionally, The NM simplex algorithm shows sensitivity to the initial point in some area. In the case of row 3– 5 in Table 5, the final z0 switch frequently between 5 and −5 though the initial points change very little. Although the NM simplex algorithm shows diversity for different initial values, if the initial point is near the optimum, the algorithm is likely to converge, as shown in Fig. 6 and double justified by Fig. 7.

4.2. Runtime environment and algorithm configuration We implemented our algorithms in the java programming language (64 bit v1.8.0_121) and ran it on a laptop workstation. The hardware and software configurations are Intel Core i7-3740QM CPU @ 2.70 GHz, DDR3 SDRAM 8 GB and Windows 10 Pro. First we applied the NM simplex algorithm only. We set stepsize h in the direction of unit vector to {500, 20, 20, 10}. We set the standard simplex transformation parameters {α , β , γ , δ} to {1, 0.5, 2, 0.5}, which was used in most implementations. Several initial points including {0, 0, 0, 0} were tested and compared. We set a bound 10000 for the total iterations to jump out from stucking. Next we applied GA or PGA only to the numerical example. we set the search domain as: Q0 ∈ {20, 6000}, x 0 ∈ {−1000, 1000}, y0 ∈ {−1000, 1000}, z 0 ∈ {0, 30}. The size of population was initially set to 80, the crossover rate was initially set to 0.9 and the mutation rate was initially set to 0.2. To ensure the PGA terminates even when it cannot obtain a good parameter estimation with a cost function value low enough, we set a default bound 20 for the total generations. For different behavioral parameters, we repeated the algorithm 50 or more times to test the performance based on the statistical analysis. Then we enabled PGA and the NM simplex algorithm successively as a whole hybrid algorithm without changing the parameters much. We did not have to specify an initial point when the process enters the NM phase, since we had already obtained a point at the end of the PGA phase. In other words, the point is the output of the PGA process and

4.4. PGA only As a global optimization algorithm, the GA finds individuals of progressively higher fitness. At the beginning of this section, we describe the default GA configuration that will be used throughout this section, as shown in Table 6. The following comparisons are made by changing one or more configurations. Because the GA is a stochastic search, most results in the form of averages and standard deviations were provided base on 50 or more times repeated computations in this section. It should be noted that the default configuration is not the best configuration for the GA. Tables 7 and 8 illustrate how the crossover rate and the mutation rate impact the GA. As usual, better source term estimation can be obtained with more crossover operations. Unlike other GAs, larger mutation rate brought better results, because the proposed mutation operator changed the chromosome less than the uniform mutation. Increasing the population size means more opportunity to obtain good individual. Table 9 supports the claim. However, it is expensive to tune these parameters, because the increase of population size is propositional to the increase of cost function evaluations. Table 10 shows how the generation bound impacts the GA. We can get better source term estimation slowly with the increase of generations. The mutation operator played an important role in the mechanism of evolution. Table 11 shows the advantage of the proposed mutation operator compared with the uniform mutation. Less change of chromosome is helpful for retaining individual characteristics. When we increased the behavioral parameter km, the estimation could be more accurate for the long-term evolution, though less accurate for the short-term evolution. When we decreased km, the result of the proposed mutation degenerated into the uniform mutation. Table 12 illustrates the impact effects of different selections. In the implementation of roulette game selection, if the parameter kf is very large, the worse individual has very little opportunity to be selected. If the parameter kf is small, the opportunity difference become small too. The roulette game selection degenerated into the truncation selection

Table 4 The observed data. No.

x(m)

y(m)

z(m)

C (g / m 3 )

1 2 3 4 5 6 7 8 9

250.0 250.0 250.0 350.0 350.0 350.0 450.0 450.0 450.0

10.0 20.0 30.0 10.0 20.0 30.0 10.0 20.0 30.0

3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0

3.2618 1.8885 0.3665 2.1362 1.6205 0.7074 1.5359 1.2997 0.7874

74

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

H. Li, J. Zhang

Fig. 4. 3D convergent trajectory and its projections of the NM simplex algorithm.

Table 5 Diverse results of the NM simplex algorithm only. Initial point

Result

source term estimation

{0, 0, 0, 0} {1000, 100, 50, 50} {800, 100, 50, 42} {800, 100, 50, 43} {800, 100, 50, 44} {15000, −150, 150, 10}

optimum stuck optimum optimum optimum sub-optimum

{1500.01, 8.00, 10.00, 5.00} {1000, 100, 50, 50} {1500.22, 7.98, 10.00, 5.00} {1499.86, 7.99, 10.00, −5.00} {1500.17, 7.97, 10.00, 5.00} {1494.83, −47.93, 6.83, 2.14}

4.5. PGA-NM hybrid algorithm Combining the global search heuristic of the PGA and the fastness of the NM simplex algorithm, the PGA-NM hybrid algorithm can obtain the optimum accurately and fast enough (see Table 14). The PGA, though faster than the GA, still shows slower convergence than other algorithms. The random starting NM algorithm ran faster than other algorithms, but diverged frequently. As a special case of PGA-NM, the Monte Carlo/random search + NM hybrid algorithm is the second faster algorithm among these algorithms. The simplification by disabling the crossover operation and applying the uniform operation reduced the computation of the Monte Carlo/random search + NM

Fig. 5. The convergent trajectory of Q0.

when kf tended to zero. Fig. 8 shows the change of the best individual in an entire GA process. By applying paralleling mechanism, the global parallel GA speeded up the searching process without changing its other performance, shown in Table 13. Parallel algorithm was executed on the i7-3740QM processor.

75

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

H. Li, J. Zhang

Fig. 6. Convergence analysis with different initial values. Here we fixed the Q, z and let x, y both change from −400 m to 400 m with step size 50 m. The green point means that the NM algorithm will converge from this initial point, while the red point means the initial point leading to diverge. The blue point {8, 10} is the global optimum.

Fig. 7. Convergence analysis with another source parameter {300, −450, 40, 5}. The green point means that the NM algorithm will converge from this initial point, while the red point means the initial point leading to diverge. The blue point {−450, 40} is the global optimum.

running time and the accuracy in the PSO algorithm. We obtained the satisfied estimation ( fc < 10−6 ) in 96.7 ms. The PSO-NM hybrid algorithm ran faster than PSO, and reached the same level as the PGA-NM. But the disadvantage of the PSO-NM

hybrid algorithm, but at the same time raised the possibility of divergence. PSO is also a stochastic optimization technique for the solution of continuous optimization problem. There is a trade-off between the 76

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

H. Li, J. Zhang

Table 6 Default configurations of the GA, PGA and PGA-NM algorithm. Item

Default configuration

selection operator crossover operator crossover rate mutation operator mutation rate

roulette wheel selection arithmetic crossover, random kc in {0, 1} 0.9 proposed in this paper, km=4 0.2

search scope population size fitness function

{(20, 6000), (−500, 500), (−500, 500), (0, 30)} 100

Table 10 Comparison of different generation bounds. Generation bound

Source estimation {μ , σ} {Q (g /s ), x (m ), y (m ), z (m )}

Cost function value fc{μ,σ}

0

{{3659.5, 1584.7}, {−196.8, 220.2} {8.8, 4.7}, {13.05, 6.59}} {{2872.7, 599.9}, {−128.9, 54.0} {2.4, 3.4}, {5.0, 1.6}} {{1995.3, 293.9}, {−101.6, 36.1} {3.2, 2.5}, {0.9, 1.3}} {{1818.6, 211.4}, {−84.9, 29.4} {4.3, 2.0}, {0.97, 1.3}}

{12.9, 25.2}

20 200 2000

k

generation bound

f = 1/fc f , k f = 1/4 20

h of the NM process

{100, 15, 15, 5}

terminate criteria

fc < 10−6

0.1 0.5 0.9

Mutation

0.05 0.2 0.5 0.8

{6.78 × 10−2, 4.42 × 10−2}

Best fc (200 gens) (μ,σ)

Best fc (2000 gens){μ,σ}

uniform mutation

{1.74 ×

10−1}

{1.19 × 10−1, 8.44 × 10−2}

proposed (km=1)

{1.46 × 10−1, 1.06 × 10−1}

{1.11 × 10−1, 6.29 × 10−2}

proposed (km=4 default)

{1.03 × 10−1, 8.04 × 10−2}

{6.55 × 10−2, 3.25 × 10−2}

10−1,

1.08 ×

Source estimation {μ, σ} {Q (g /s ), x (m ), y (m ), z (m )}

Cost function value fc{μ,σ}

{{3296.5, 1030.7}, {−169.8, 92.6} {2.9, 5.6}, {5.8, 2.2}} {{2871.6, 646.1}, {−121.9, 46.8} {3.5, 3.1}, {5.5, 1.8}} {{2612.3, 503.1}, {−110.3, 40.3} {3.9, 2.5}, {4.9, 1.8}}

{1.26, 1.03}

proposed (km=128)

{2.25 × 10−1, 1.61 × 10−1}

{5.73 × 10−3, 4.69 × 10−3}

{4.85 × 10−1, 2.53 × 10−1}

proposed (km=1024)

{3.97 × 10−1, 2.81 × 10−1}

{1.65 × 10−1, 1.72 × 10−1}

{3.79 × 10−1, 2.24 × 10−1} Table 12 Comparison of different fitness functions.

Table 8 Comparison of different mutation rates. Mutation rate

{1.14 × 10−1, 6.94 × 10−2}

Table 11 Comparison of different mutations.

Table 7 Comparison of different crossover rates. Crossover rate

{4.11 × 10−1, 1.98 × 10−1}

Best fc (20 gens) (μ,σ)

Best fc (2000 gens) (μ,σ)

roulette, kf=64

{9.42 × 10−1, 9.98 × 10−1}

{9.46 × 10−2, 5.78 × 10−2}

roulette, kf=4

{6.01 × 10−1, 4.36 × 10−1}

{9.98 × 10−2, 5.92 × 10−2}

Cost function value fc{μ,σ}

Source estimation {μ , σ} {Q (g /s ), x (m ), y (m ), z (m )}

{{2895.4, 484.4}, {−116.3, {3.3, 2.4}, {5.7, 1.2}} {{2770.7, 509.4}, {−121.0, {2.9, 2.8}, {5.3, 1.6}} {{2503.2, 610.8}, {−119.5, {2.8, 3.6}, {3.7, 1.8}} {{2233.4, 431.4}, {−107.9, {3.2, 3.1}, {2.9, 1.5}}

Selection

38.8}

{4.45 × 10−1, 1.99 × 10−1}

roulette, kf=1

{3.81 × 10−1, 2.91 × 10−1}

{8.14 × 10−2, 4.87 × 10−2}

47.7}

{4.05 × 10−1, 2.11 × 10−1}

{3.69 × 10−1, 1.95 × 10−1}

{6.61 × 10−2, 4.24 × 10−2}

56.6}

{2.91 × 10−1, 2.07 × 10−1}

roulette, k f = 1/4 (default)

50.2}

{1.88 × 10−1, 1.12 × 10−1}

roulette, k f = 1/256

{4.25 × 10−1, 2.15 × 10−1}

{6.79 × 10−2, 3.63 × 10−2}

truncation selection

{3.88 × 10−2, 1.69 × 10−1}

{7.27 × 10−2, 3.07 × 10−2}

algorithm is its lower global searching ability than PGA-NM. When we enlarged the searching scope to x: {−1000,1000} and y: {−1000,1000}, the PSO-NM algorithm had a distinct possibility to diverge, even if we increased the PSO generation bound to 1000. In contrast, the PGA-NM with generation bound 20 did work well in the enlarged searching range. No divergence happened in 200000 repetitions. The configurations of the PSO algorithm and the PSO-NM algorithm in this paper include: Neumann neighborhood topology, the PSO size was 20, inertia weight was 0.7, cognitive coefficient was 1.0, and social coefficient was 1.4.

4.6. Test with measurement error Next, we tested the influence of the measurement error in the PGANM algorithm. Case 1 is under the previous condition with no measurement error, see Table 4. In case 2, we added a random 0– 20% increase to the value of each sensor. In case 3, we introduced a random 0–20% decrease to the value of each sensor. In case 4, we randomly increased or decreased 0–20% the observed concentrations. Table 15 compares the results of the PGA-NM algorithm of case 1,2,3 and 4. It demonstrates that the source term estimations are acceptable in case 2–4, and the computation time are much longer but still acceptable.

Table 9 Comparison of different population sizes. Population size

Source estimation {μ , σ} {Q (g /s ), x (m ), y (m ), z (m )}

Cost function value fc{μ,σ}

10 100

{{3422.5, 1253.1}, {−206.7, 145.4} {10.1, 9.7}, {9.1, 4.6}} {{2727.1, 510.3}, {−113.2, 40.9} {3.4, 2.3}, {5.2, 1.5}}

{5.72, 6.29}

1000

{{2101.6, 286.9}, {−88.2, 34.9} {4.8, 2.11}, {3.9, 1.1}}

{1.62 × 10−1, 5.94 × 10−2}

77

{3.64 × 10−1, 2.06 × 10−1}

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

H. Li, J. Zhang

Fig. 8. {Q0 , x 0 , y0 , z0} versus generations for PGA. Table 13 Comparison between the global PGA and the serial GA.

Table 15 Tests with different measurement errors.

Implementations

Population size

Average time consumed (ms)

Case No.

GA only GA only PGA only PGA only

100 1000 100 1000

15.0 283.5 6.4 121.2

4.7. Tests on the embedded system

Source term estimation

Cost function value

time (ms)

1

{1499.9, 7.9, 9.9, 4.9}

8.31 × 10−7

16.2

2

{1512.3, 1.5, 9.8, 4.4}

2.82 × 10−2

309

3

{1349.6, 7.2, 9.7, 4.6}

2.27 × 10−2

307

4

{1348.9, 7.5, 9.4, 4.0}

1.13 × 10−1

308

comparison was made based on the averages of 50 repetitions too.

Java succeeded at running our program on both the Windows system and many other systems by being platform-independent at both the source and binary levels. The performance of our PGA-NM algorithm on the embedded system are shown in Table 16. The

5. Conclusion In this paper, we addressed the source term estimation problem by optimizing the fit between the modeled dispersion and the monitored

Table 14 Comparison of the PGA-NM algorithm and its competitive algorithms. Algorithm

Average best estimation {Q0 , x 0 , y0 , z0}

Average best cost function value

Time (ms)

PGA only (generation bound=20)

{2690.7, −118.9, 3.2, 4.9}

3.53 × 10−1

6.5

PGA only (generation bound=2000)

{1711.3, −70.7, 5.3, 1.41}

5.81 × 10−2

581.4

Random starting+NM

{1497.6, 7.9, 10.1, 4.8} or divergence

9.94 × 10−7

9.0

Monte Carlo/random search (generation bound=10) +NM

{1499.2, 7.9, 10.0, 4.8} or divergence

8.22 × 10−7

10.9

PGA-NM in this paper (GA generation bound=10)

{1499.9, 8.1, 10.0, 4.9}

8.13 × 10−7

11.3

PGA-NM, enlarged range (GA generation bound=20)

{1499.9, 8.1, 10.1, 5.0}

8.53 × 10−7

13.2

PSO (PSO generation bound=2000)

{1500.1, 7.8, 9.9, 4.9}

7.76 × 10−7

96.7

PSO-NM (PSO generation bound=30)

{1499.9, 7.9, 10.0, 4.9}

8.43 × 10−7

13.4

PSO-NM, enlarged range (PSO generation bound=300)

{1497.6, 7.9, 10.1, 4.9} or divergence

8.23 × 10−7

37

78

Engineering Applications of Artificial Intelligence 62 (2017) 68–79

H. Li, J. Zhang

Cant-Paz, E., 1999. A survey of parallel genetic algorithms. Calc. Paralleles Reseau. Et. Syst. Repartis 10 (4), 429–449. Chow, F.K., Kosovic, B., Chan, S., 2008. Source inversion for contaminant plume dispersion in urban environments using building-resolving simulations. J. Appl. Meteorol. Climatol. 47 (6), 1553–1572. Crowl, D.A., Chemical process safety: Fundamentals with applications, 2/e. Gaffney, J., Green, D.A., Pearce, C.E.M., 2009. Binary versus real coding for genetic algorithms: A false dichotomy? Anziam J. 51, C347–C359. Garca, M.L.L., Garca-Rdenas, R., Gmez, A.G., 2014. Hybrid meta-heuristic optimization algorithms for time-domain-constrained data clustering. Appl. Soft Comput. 23 (5), 319–332. Gifford, F., Jr, 1976. Consequences of effluent releases. Nucl. Saf. 17 (1), 68–86. Guha, A., 1997. A unified eulerian theory of turbulent deposition to smooth and rough surfaces. J. Aerosol Sci. 28 (8), 1517–1537. Haddad, O.B., Hamedi, F., Fallah-Mehdipour, E., Orouji, H., Mario, M.A., 2015. Application of a hybrid optimization method in muskingum parameter estimation. J. Irrig. Drain. Eng. 141 (12), 482–489. Haupt, S.E., Beyer-Lout, A., Long, K.J., Young, G.S., 2009. Assimilating concentration observations for transport and dispersion modeling in a meandering wind field. Atmos. Environ. 43 (6), 1329–1338. Haupt, S.E., 2005. A demonstration of coupled receptor dispersion modeling with a genetic algorithm. Atmos. Environ. 39 (37), 7181–7189. Hosseini, B., Stockie, J.M., 2016. Bayesian estimation of airborne fugitive emissions using a gaussian plume model. Atmos. Environ. 141, 122–138. Lipowski, A., Lipowska, D., 2011. Roulette-wheel selection via stochastic acceptance. Phys. A Stat. Mech. Appl. 391 (6), 2193–2196. Ma, D., Deng, J., Zhang, Z., 2013. Comparison and improvements of optimization methods for gas emission source identification. Atmos. Environ. 81 (2), 188–198. Manderick, B., Spiessens, P., 1989. Fine-grained parallel genetic algorithms. In: International Conference on Genetic Algorithms. Michalewicz, Z., 1994. Genetic algorithms + data structures = evolution programs. Comput. Stat. Data Anal. 24 (3), 372–373. Mondal, S., Maiti, M., 2003. Multi-item fuzzy eoq models using genetic algorithm. Comput. Ind. Eng. 44 (1), 105–117. Nelder, J.A., Mead, R., 1965. A simplex method for function minimization. Comput. J. 7 (4), 308–313. Oettl, D., Kukkonen, J., Almbauer, R.A., 2001. Evaluation of a gaussian and a lagrangian model against a roadside data set, with emphasis on low wind speed conditions. Atmos. Environ. 35 (12), 2123–2132. Pasquill, F., 1961. The estimation of the dispersion of windborne material. Meteorol. Mag. 90 (1063), 33–49. Pettey, C.B., Leuze, M.R., Grefenstette, J.J., 1995. A parallel genetic algorithm. Syst. Eng. 367 (3), 155–161. Przystaka, P., Katunin, A., 2016. Multi-Objective Meta-Evolution Method for Large-Scale Optimization Problems. Schmidpeter, A., Lochschmidt, S., Willhalm, A., 2015. Numerical investigation of concentrations and concentration fluctuations around isolated obstacles of different shapes. comparison with wind tunnel results. Environ. Fluid Mech. 15 (5), 1–36. Singh, S.K., Sharan, M., Issartel, J.P., 2013. Inverse modelling for identification of multiple-point releases from atmospheric concentration measurements. Bound.Layer. Meteorol. 146 (2), 277–295. Sokhi, R., Fisher, B., Lester, A., McCrae, I., Bualert, S., Sootornstit, N., 1998. Modelling of air quality around roads. In: Proceedings of the 5th International Conference On Harmonisation with Atmospheric. Stern, R., Yamartino, R.J., 2001. Development and first evaluation of micro-calgrid: a 3d, urban-canopy-scale photochemical model. Atmos. Environ. 35 (14), (149–S165). Sutton, O., 1947. The problem of diffusion in the lower atmosphere. Q. J. R. Meteorol. Soc. 73 (317–318), 257–281. Venkatakrishnan, V., 2015. Newton solution of inviscid and viscous problems. AIAA J. 27 (7), 885–891. Wright, M.H., 1997. Direct search methods: Once scorned, now respectable. Yu, N., 2015. Universal gradient methods for convex optimization problems. Math. Program. 152 (1–2), 381–404.

Table 16 Tests on the embedded system. Algorithm

Average best cost function value

Average time consumed (ms)

PGA-NM (GA generations=10) PGA-NM (GA generations=20) PGA-NM (GA generations=50)

7.93 × 10−7

85

8.12 × 10−7

110.5

8.06 × 10−7

171.3

sensor data. We tackled this optimization problem by using the PGANM hybrid algorithm, which combines the global search ability of the genetic algorithm, the direct search feature of the NM simplex algorithm, and the parallel mechanism. We implemented the proposed algorithm in java language on the Raspbian OS of an ARM Cortex-A7 quad-core processor. The running time of less than 200 ms are considerably shorter than 14.29 s in Table 1 of paper (Ma et al., 2013). We showed in comparison to other mainstream optimization algorithms that the PGA-NM algorithm is able to return a satisfying source term estimation more fast or more stably. The PGA-NM algorithm ran fast than GA, PGA and PSO owing to its NM process. The PGA-NM algorithm is stable than random starting+NM, random search+NM, PSO-NM owing to its global search efficiency. Although the accuracy of the proposed algorithm still depends on the accuracy of the dispersion model like other competitive algorithms, the algorithm has appreciable tolerance to measurement errors. In our further studies we plan to apply the PGA-NM hybrid algorithm to establishing a new atmospheric dispersion model. Acknowledgment This work was supported by the State Commission of Science Technology of China under grants 2015 The national key technology R & D program (Project No. 2015BAK39B02). And this project was supported by the NSFC general technical foundation research joint fund (Project No. U1636208) too. References Abramson, D., Abela, J., 1991. A parallel genetic algorithm for solving the school timetabling problem, pp. 1–11. Agency, E.P., Guideline on air quality models. Appendix W to Part 51 Title 40 Protection of the Environment. Baghmisheh, M.T.V., Peimani, M., Sadeghi, M.H., Ettefagh, M.M., Tabrizi, A.F., 2012. A hybrid particle swarm nelder mead optimization method for crack detection in cantilever beams. Appl. Soft Comput. 12 (8), 2217–2226. Bonifacio, H.F., Maghirang, R.G., Glasgow, L.A., 2015. Numerical simulation of transport of particles emitted from ground-level area source using aermod and cfd. Eng. Appl. Comput. Fluid Mech. 8 (4), 488–502.

79