Application of Particle Swarm Optimization to ... - ACS Publications

ARTICLE pubs.acs.org/IECR

Application of Particle Swarm Optimization to Fourier Series Regression of Non-Periodic Data Eldin Wee Chuan Lim* Department of Chemical and Biomolecular Engineering, National University of Singapore, 4 Engineering Drive 4, Singapore 117576 ABSTRACT: A new methodology for solving black-box optimization problems by the continuous approach has been developed in this study. A discrete Fourier series method was derived and used for reformulation of black-box objective functions as continuous functions. Particle swarm optimization (PSO) was then applied to locate the global optimal solutions of the continuous functions derived. The continuous functions generated by the proposed discrete Fourier series method correlated almost exactly with their original black-box counterparts. The PSO algorithm was observed to be highly successful in achieving global optimization of all such objective functions considered in this study. Case studies were also carried out in which the methodology developed here was applied to the optimization of the formulation for culture media used for cultivating a freshwater microalga, Haematococcus pluvialis, as well as for a one-dimensional regression problem. Subsequently, the methodology was extended to solving a two-dimensional black-box optimization problem that was simulated based on the Rosenbrock function. The possibility of applying the same approach for solving the chemical plume tracing problem in the area of chemical defense was also successfully demonstrated. The results obtained indicate that the discrete Fourier series method coupled with the PSO algorithm is indeed a promising methodology for solving black-box optimization problems by a continuous approach.

’ INTRODUCTION Black-box optimization refers to the maximization or minimization of an objective function over a set of feasible parameter values where the objective function cannot be evaluated analytically. Such optimization problems can arise in situations where a complex system whose behavior is not well understood and not formally described needs to be optimized.1 In such problems, values of the objective function might have to be measured, estimated, or evaluated from simulations, and the domain(s) of the objective function and/or feasible region is usually noncontinuous in nature. Such optimization problems are ubiquitous in engineering, operations research, and computer science. However, research on black-box optimization techniques seems to be limited largely to the domain of theoretical computer science.2 In particular, the concepts and techniques developed to date have not been applied to solving problems that are of interest to chemical engineers. Audet and Orban3 devised a framework for identifying locally optimal algorithmic parameters by minimizing some measure of performance of the algorithm being optimized. The measure of performance was treated as a black box, and a derivative-free method known as mesh adaptive direct search was applied. It was shown that the optimization strategy allowed the use of a surrogate function, which is a simplification of the real objective function that exhibits similar behavior but is less costly to evaluate, to guide the search process. Driessen et al.4 incorporated the D-optimality criterion into various sequential derivative-free optimization algorithms for solving black-box optimization problems. The criterion was used to define a trust region that adapts its shape to the locations of points in which the objective function has been evaluated and to determine an optimal geometryimproving point. More recently, Ramani et al.5 considered the problem of optimizing the parameters of a given denoising r 2010 American Chemical Society

algorithm for restoration of a signal corrupted by white Gaussian noise. A black-box optimization approach that was based on a Monte Carlo technique was applied to minimize a quantity referred to as Stein's unbiased risk estimate. The applicability of the method to denoising algorithms in the wavelet and variational settings was demonstrated. Here, a new approach for solving black-box optimization problems that draws analogy from concepts rooted in the realm of discrete optimization is presented. In what follows, a brief review of the subject of discrete optimization and techniques used to transform such problems into equivalent continuous ones is first provided, and then the proposed analogous approach for solving black-box optimization problems is described. One approach to solving stochastic optimization problems with continuous decision variables is to apply a stochastic approximation algorithm to locate a root of the gradient of the objective function. Estimates of the gradient can be obtained by using finite differences, the likelihood ratio method, or perturbation analysis. Andradottir6 presented a method for stochastic discrete optimization that resembled a stochastic approximation algorithm that generated a sequence of estimates of the solution where each new estimate was obtained from the previous estimate by taking a small step in the direction of the negative of the gradient of the objective function. Gong et al.7 proposed a stochastic comparison algorithm that was designed for discrete optimization problems where only noisy estimates of the objective function were available and the search space was very large. Received: December 1, 2009 Accepted: December 9, 2010 Revised: December 1, 2010 Published: December 29, 2010 2307

dx.doi.org/10.1021/ie101399r | Ind. Eng. Chem. Res. 2011, 50, 2307–2322

Industrial & Engineering Chemistry Research The main advantage of this algorithm was that it was essentially parameter-free and required no knowledge about the structure of the search space. Martinelli8 presented a stochastic comparison algorithm for the optimization of time-varying objective functions defined on a discrete set whose values were known only through estimates. The algorithm was derived from a stochastic comparison algorithm by introducing a noise filter that made bad movements more difficult, without compromising good movements. Kleywegt et al.9 studied a Monte Carlo simulation-based approach to stochastic discrete optimization problems where a random sample was generated and the expected value function was approximated by the corresponding sample average function. The sample average optimization problem was then solved, and the procedure was repeated until a stopping criterion was satisfied. Alrefaei and Andradottir10 proposed a modified stochastic ruler method that showed better overall performance than the original stochastic ruler method. It involved generating a Markov chain sequence using values in the feasible set of the discrete optimization problem and used the number of visits of this sequence to the different states to estimate the optimal solution. This method was subsequently further improved by the same authors11 by considering the state that had the best average estimated objective function value obtained from all previous observations of the objective function values as the estimate of the optimal solution. Solution methods for discrete optimization problems can generally be classified into combinatorial and continuous approaches. In the former, a sequence of states is generated from a discrete finite set to represent a partial solution, whereas in continuous approaches, the discrete optimization problem is characterized using equivalent continuous formulations or continuous relaxations.12 There are many methods to formulate discrete optimization problems as equivalent continuous ones. For example, Goyal and Ierapetritou13 recently proposed combining a simplicial-based approach and the sample average approximation approach for solving convex stochastic mixed-integer nonlinear programming (MINLP) problems. This was based on the idea of closely describing the feasible region by a set of linear constraints representing an approximation of its convex hull. The objective function was linearized at the simplicial boundary points, and the global optimal solution was obtained using the linear representation of the feasible space and the linear approximation of the objective function. Aytug and Sayin14 used support vector machines (SVMs) to learn the efficient set of a multiple-objective discrete optimization problem. They conjectured that a surface generated by SVMs could provide a good approximation to the efficient set. To evaluate their idea, the authors incorporated the SVM-learned efficient set into a multiobjective genetic algorithm through the fitness function to generate feasible solutions that were representative of the true efficient set. In this article, a new continuous approach to the solution of black-box optimization problems is presented. As is typical in such problems, the objective function is composed of a set of points that can be assumed to be generated from simulations or experiments. It can be assumed that an analytical function representation of the objective function is unavailable, regardless of whether such a function is derivable by other theoretical or modeling means. In principle, the full objective function can be derived empirically by conducting an infinite number of such simulations or experiments to obtain the complete set of data points. However, such an approach is almost always impractical, as these simulations or experiments are usually expensive to conduct. The goal of the optimization process is to obtain the value(s) of the independent variable(s) that will maximize

ARTICLE

(or minimize) the objective function, keeping in mind that the optimal solution might not correspond to any of the available data points representing the objective function. Currently available methods that make use of linearization, smoothing, or averaging techniques to derive an equivalent continuous formulation of the black-box objective function are approximations to the true analytical model at best. On the other hand, it is also wellestablished that a function (continuous or otherwise) can be represented as a Fourier series consisting of a large number of trigonometric terms. In the limit that the number of such terms approaches infinity, the Fourier series will converge exactly to the original function it represents, thus making it effectively a universal mathematical regression technique. Based on this premise, it was conjectured that a Fourier series representation of a black-box objective function can be derived that will correspond exactly or almost exactly to the true analytical model and that also effectively converts the optimization problem involving black-box functions into a continuous problem. A stochastic global optimization method can then be applied to solve the continuous optimization problem. Such an approach that exploits the superior regression capabilities of the Fourier series method to solving black-box optimization problems has not been attempted in the research literature to date. In the next section, the formulation of the Fourier series method as applied to black-box functions is derived, and the basic principles of a stochastic method, particle swarm optimization, that will be applied to solving the resulting continuous optimization problem are introduced.

’ COMPUTATIONAL METHODS Discrete Fourier Analysis. A periodic function, f(x), with period L can be expressed as a Fourier series given by the expression ¥ X nπx nπx f ðxÞ ¼ a0 þ an cos þ bn sin ð1Þ L L n¼1

where n = 1, 2, 3, ..., and the coefficients are given by the integrals Z 1 L a0 ¼ f ðxÞ dx ð2Þ L 0 Z 2 L nπx f ðxÞ cos dx ð3Þ an ¼ L 0 L Z 2 L nπx bn ¼ f ðxÞ sin dx ð4Þ L 0 L For mathematical analysis, the Fourier series can also be expressed in exponential form with complex coefficients as ¥ X 2πnx f ðxÞ ¼ Rn exp i ð5Þ L n¼ -¥ where Rn = 1/2(an - ibn), n = 0, ( 1, ( 2, ..., and Rn is given by Z 1 L=2 2πnx Rn ¼ f ðxÞ exp i dx ð6Þ L - L=2 L By the convolution theorem, it can be shown that the coefficients Rn of the Fourier series expansion of f(x) are equal to 2308

dx.doi.org/10.1021/ie101399r |Ind. Eng. Chem. Res. 2011, 50, 2307–2322

Industrial & Engineering Chemistry Research

ARTICLE

the values of the corresponding Fourier transform of f(x), denoted as F(f), evaluated at n/L. In other words, the Fourier transform of a periodic function is then an infinite set of sinusoids (i.e., an infinite sequence of equidistant impulses) with amplitudes F(n/L).15 Consider a continuous function f(x) and its Fourier transform F(f) that are discretized by sampling and truncation so that only a finite number of points, say, N, are considered. If it is assumed that the N samples of the original function f(x) are one period of a periodic waveform, then the Fourier transform of this periodic function is given by the N samples as computed by the expression N -1 X n 2πnk f ðkl Þ exp - i ¼ F ð7Þ Nl N k¼0 where l = L/N and n = 0, 1, ..., N - 1. This is the well-established discrete Fourier transform.15 The sampling theorem states that, if the Fourier transform of a function f(x) is zero for all frequencies greater than a certain frequency fc, then the continuous function f(x) can be uniquely determined from knowledge of its sampled values. Consider a black-box function containing N data points for which the corresponding continuous function is to be determined. Without loss of generality, the smallest interval between any two adjacent data points can be used to define the highest-frequency component of the Fourier transform of the corresponding continuous function. Further, the domain of the black-box function, say, 0 e x e L, can be assumed to be one period of a periodic waveform of the corresponding continuous function. This allows eq 7 to be applied for calculating values of the discrete Fourier transform of the black-box function, and by the convolution theorem mentioned earlier, these values will be equal to the coefficients of the Fourier series expansion of the corresponding continuous function. By this argument, the following expressions for the Fourier cosine half-series of a continuous function can be derived from its black-box counterpart N -1 X nπx ð8Þ f ðxÞ ¼ a0 þ an cos L n¼1 " # -1 1 f ð0Þ þ f ðLÞ NX k þ f a0 ¼ N 2 N k¼1

ð9Þ

" # -1 2 f ð0Þ þ f ðLÞ cosðnπÞ NX k nπk þ an ¼ f cos N 2 N N k¼1 ð10Þ Particle Swarm Optimization. The most challenging global optimization problems are those that do not contain any known structure that can be exploited. However, stochastic methods, which generally require few assumptions about the optimization problem, are particularly suited for such problems. The particle swarm optimization (PSO) algorithm is a population-based search algorithm based on the simulation of the social behavior of birds within a flock.16,17 This is a

well-established stochastic global optimization algorithm, and only a brief description is presented here for the sake of completeness. The position and velocity of each particle, Pi, in the swarm are updated at each time step according to the equations xi ðtÞ ¼ xi ðt - 1Þ þ vi ðtÞΔt

ð11Þ

vi ðtÞ ¼ φvi ðt - 1Þ þ F1 ½xpbest, i - xi ðt - 1Þ þ F2 ½xgbest - xi ðt - 1Þ

ð12Þ where xi(t) and vi(t) are the position and velocity, respectively, of particle Pi at time t; Δt is the time step, which is assigned an arbitrary value of 1.0; φ is an inertia weight; F1 and F2 are random variables; xpbest,i is the position giving the best performance of Pi up to time t; and xgbest is the position giving the globally best performance of the entire swarm up to the current time. In the current implementation, particles that attempt to travel out of the domain of the function to be optimized are constrained to remain stationary at the edge of the domain until they reverse their direction of motion. The second and last terms above are commonly referred to as the cognitive and social components, respectively. The farther a particle is from its own best solution and the global best solution, the larger the driving force provided by these components to move the particle toward the best solutions. The random variables F1 and F2 are defined as F1 = r1c1 and F2 = r2c2, where r1 and r2 are randomly selected from a uniform distribution U(0,1) and c1 and c2 are positive acceleration constants. Kennedy18 has shown that c1 þ c2 e 4 is a necessary condition in this formulation. Otherwise, velocities and positions tend to diverge to infinity. The inertia weight φ provides improved performance of the algorithm by controlling the influence of previous velocities on the new velocity. Physically, larger inertia weights cause larger exploration of the search space and vice versa. However, studies have also shown that the PSO algorithm does not converge for all combinations of the inertia weight and the acceleration constants c1 and c2. To ensure convergence, the following relation must hold19 1 φ > ðc1 þ c2 Þ - 1 2

ð13Þ

with φ e 1. The PSO algorithm tends to exhibit cyclic or divergent behavior if the above relation is not satisfied. In the present study, the values used for the above parameters were c1 = 1.0, c2 = 1.0, and φ = 0.5. The PSO algorithm has been applied to a variety of problems in both chemical engineering and other research areas. For example, Li et al.20 proposed a new multiobjective particle swarm optimization (MOPSO) procedure based on the Pareto dominance hybrid algorithm and applied it to a naphtha industrial cracking process. The proposed algorithm was observed to have better convergence and diversity Pareto solutions than the modified version of the nondominated sorting genetic algorithm (NSGAII) method. Tang and Yan21 investigated a typical campaign planning problem (CPP) and proposed a hybrid approach comprising heuristic and PSO algorithms whereby PSO was applied to solve one subproblem with binary variables whereas the heuristic algorithm was applied to another subproblem with the remaining variables. The computational results showed that the proposed hybrid 2309



ARTICLE

algorithm is superior in terms of solution speed and quality of solutions obtained.

’ RESEARCH APPROACH The approach undertaken in this study to validate the effectiveness of the proposed methodology for deriving continuous formulations of black-box optimization problems and subsequently solving such continuous problems with the PSO algorithm can be considered to consist of three parts. First, general black-box functions were constructed by generating N random data points with a single independent variable. Without loss of generality, the domain and range of all such functions generated were 0.0 e x e 1.0 and 0.0 e f(x) e 1.0, respectively. The discrete Fourier series method, consisting of eqs 8-10 presented previously, was applied to derive the continuous equivalent of the black-box function. The agreement between the discrete Fourier series and the original black-box function was validated either visually through a graphical approach or by computation of the mean absolute error. This validation process was repeated for different numbers of random data points used to represent each black-box function, N, that span over 5 orders of magnitude: N = 10, 102, 103, 104, and 105. In the second part of this study, the effectiveness of the PSO algorithm in locating the maximum point of the continuous function representation of the original black-box function was investigated systematically. It is well-known that the number of particles, P, used in such a search process has strong implications for the number of iterations required to locate the optimal solution. In correspondence with the five types of black-box functions used in the previous part defined based on the order of magnitude of N, five values of P spanning the same 5 orders of magnitude, P = 10, 102, 103, 104, and 105, were applied for each type of black-box function, giving a total of 25 simulation cases. In each case, particles were first randomly distributed over the entire domain of the objective function and initialized with random velocities. Subsequently, the positions and velocities of all particles were updated automatically by the PSO algorithm described in the previous section. The algorithm was allowed to undergo 1000 iterations, which was more than sufficient for the optimal solution to be located in all cases considered. It can also be noted at this point that, despite the noisy functions considered in this study, the initial distribution of particles was found to have insignificant influence on the performance of the PSO algorithm in terms of the number of iterations required to locate the optimal solution. The algorithm was tested with different uniform and random initial distributions of particles for each case during the initial phase of this study, and the performance of the optimization algorithm was observed to be insensitive to its initial condition. The methodology described above was applied to a case study based on an experimental work reported by Fabregas et al.22 In this experimental work, the researchers implemented a single-variable optimization strategy to maximize the productivity of vegetative cells of the microalga Haematococcus pluvialis. The optimization strategy was applied to 12 types of nutrients found in the culture media, and the steady-state cell density obtained was shown to be higher than that obtained with traditional media. In the present study, the discrete sets of experimental data reported by Fabregas et al.22

Figure 1. Comparisons between black-box functions and their corresponding discrete Fourier series for (a) N = 10 and (b) N = 102. (c) Comparison of random data values generated with corresponding function values of the discrete Fourier series for N = 103.

were used to derive the continuous equivalent of the optimization problem, and a new optimal formulation for the culture medium was obtained and compared with that reported by the 2310



ARTICLE

Figure 2. Stochastic search process to locate the global maximum of the discrete Fourier series representation of the black-box objective function by the PSO algorithm. The number of data points in the black-box function was N = 10, and the number of particles used in the PSO algorithm was P = 10. Each panel represents two iterations of the PSO algorithm.

previous researchers, as well as traditional formulations of culture media used for cultivating this microalga. Next, a case study that involved a one-dimensional regression problem

that was used recently by Coleman and Block23 to illustrate the ability of a Bayesian regularized feed-forward neural network to predict global optima was considered. The simulated 2311



ARTICLE

Figure 3. Stochastic search process to locate the global maximum of the discrete Fourier series representation of the black-box objective function by the PSO algorithm. The number of data points in the black-box function was N = 102, and the number of particles used in the PSO algorithm was P = 10. Panels are presented at intervals of 20 iterations of the PSO algorithm.

data points generated by the researchers were used to derive a continuous function representation of the problem, and the global maximum was then obtained.

With the completion of the above case studies, it was deemed prudent to extend the current study to consideration of higher-dimensional black-box optimization problems. 2312



ARTICLE

Figure 4. Performance chart of the PSO algorithm in solving optimization problems involving black-box functions of sizes N = 103, 104, and 105 that had been formulated as continuous problems using the discrete Fourier series methodology. The numbers of particles used in the search for the global maxima of the continuous functions were P = 10, 102, 103, 104, and 105, as indicated on the horizontal axis.

To this end, the two-dimensional Rosenbrock function was used as a test function to simulate a black-box optimization problem f ðx, yÞ ¼ ð1 - xÞ2 þ 100ðy - x2 Þ2

ð14Þ

Without loss of generality, the part of the function in the domains 0 e x e 2 and 0 e y e 2, where the global minimum is known to be located at (x, y) = (1, 1), was considered. For simulating a black-box optimization problem, it was assumed that the above function was unknown analytically, and only function values at specific points within the domains were available. The two-dimensional continuous equivalent of this black-box function was then derived by applying a discrete double Fourier sine series ! M X N X mπx nπy amn sin sin ð15Þ f ðx, yÞ ¼ Lx Ly m¼1 n¼1 amn ¼

M X N 4 X h k mπh nπk , f sin sin ð16Þ MN h ¼ 1 k ¼ 1 M N M N

Subsequently, the global minimal point located by applying the PSO algorithm to this discrete double Fourier sine series representation of the Rosenbrock function was then compared with the true global minimum of the function. In the final part of this study, the possibility of applying the black-box optimization methodology presented here to solving a long-standing problem in the area of chemical defense, referred to as the chemical plume tracing problem, is illustrated. During a toxic chemical release and dispersion incident, the imperative need of first responders is to determine the physical location of the source of chemical release in the shortest possible time, so that appropriate countermeasures can be carried out in a timely fashion so as to minimize the extent of human casualties. The first

Figure 5. Continuous equivalents of the optimization problem investigated experimentally by Fabregas et al.22 Nutrient concentrations were considered in multiples of the most common culture medium used.24 Colored symbols indicate global optimal solutions found by the PSO algorithm. The optimization procedure was carried out for 12 types of nutrients typically found in culture media used for cultivating the microalga Haematococcus pluvialis and included (a) Mg, Ca, Fe, and Co; (b) Cu, Cr, Mn, and Mo; and (c) Se, thiamine, B12, and biotin. 2313



ARTICLE

Table 1. Comparisons of Formulations of Optimized Medium for Haematococcus pluvialis Derived Experimentally22 and by BlackBox Optimization (Present Work) with Three Culture Media Commonly Used for Freshwater Microalgae: Most Common Culture Medium,24 Bold Basal Medium (BBM), and Modified Chu Medium (CHU) experiment22

nutrient -1

MgSO4 3 7H2O (mg L )

present work

Nichols24

BBM

CHU

246.5

493.0

4.93

75.0

36.9

CaCl2 3 2H2O (mg L-1)

110.9

222.0

2.22

25.0

36.7

2.62

37.1

5.25

33.5

CoCl2 3 6H2O (mg L-1)

0.011

0.0078

0.024

0.02

CuSO4 3 5H2O (mg L-1) Cr2O3 (mg L-1)

0.012 0.076

0.0067 0.24

0.025 0.015

1.57

MnCl2 3 4H2O (mg L-1)

0.989

1.52

0.198

1.44

0.12

0.075

0.242

FeC6H5O7 3 5H2O (mg L-1)

Na2MoO4 3 2H2O (mg L-1)

SeO2 (mg L-1)

0.005

1.1

0.011

thiamine (μg mL-1)

17.5

80.4

35.0

B12 (μg mL-1)

15.0

231.0

30.0

biotin (μg mL-1)

25.0

349.0

50.0

responders to such an incident might consist of ground troops from the defense forces or armed forces that are equipped with the necessary chemical suits and detectors. A massive search operation to locate the source of chemical release is required before containment, and subsequently, decontamination procedures can be carried out. Alternatively, to minimize human exposure to such potentially lethal chemicals, advanced technologies such as robots and sensors can be exploited to carry out the task of locating the position of the source. It is envisaged that such technologies, when coupled with the necessary artificial intelligence afforded by appropriate numerical algorithms, will be able to accomplish such tasks with minimal or even no human intervention. However, the requirement for an effective numerical algorithm to complement a sensor network or fleet of autonomous robots installed with chemical sensors for carrying out chemical plume tracing has received little attention to date. The development of such a capability would provide all nations throughout the world with a new technological edge to enable fast response to toxic chemical release scenarios.

’ RESULTS AND DISCUSSION Figure 1a-c shows that good agreements between the discrete Fourier series and the respective black-box functions were achieved for the cases N = 10, 102, and 103. It can also be observed visually that the discrete Fourier series was successful in representing these black-box functions in the simplest possible continuous fashion. The mean absolute errors of the regression analyses for these three cases were 16%, 5%, and 0.4%, respectively. Here, the mean absolute error was obtained by first calculating the absolute percentage difference between the values of f(x) for the original black-box function and that for the discrete Fourier series at the same values of x and averaging such absolute percentage differences over the entire set of data points. This is used as an indication of the goodness of fit between the discrete Fourier series function derived and the original black-box function. For N = 104 and 105, because of the large numbers of data points, the spatial extensions of the corresponding figures are too large for them to be adequately presented here. However, it can be stated here that the corresponding mean absolute errors for these two cases were 3% and 0.08%, respectively, indicating that the methodology was equally successful for black-box functions containing large numbers of points. In fact, it can be assumed that

0.02 0.0126 0.0126

no restriction exists on the size of the black-box function to be formulated as a continuous function using such an approach as the discrete Fourier series method, which was derived from the conventional Fourier series and principles associated with the discrete Fourier transform. As mentioned earlier, to the extent that a certain number of data points are available to carry out the regression analysis, the discrete Fourier series so derived is expected to provide a good representation of the true analytical model of the system under consideration. More importantly, the results obtained here support the claim mentioned earlier that a discrete Fourier series method can indeed be used to transform a black-box optimization problem into a continuous one. As mentioned in the previous section, a total of 25 cases of discrete Fourier series regression followed by PSO global optimization of the continuous representation of the black-box objective function were investigated in this study. For brevity and clarity of presentation, only two cases (N = 10, P = 10) and (N = 102, P = 10) are presented in detail graphically to illustrate the stochastic search process for the global maximum of the respective continuous functions. The remaining cases are subsequently presented in the form of a performance chart. Figure 2 shows that particles in the first case, initially randomly distributed over the entire domain, were able to converge fairly quickly to the highest point of the discrete Fourier series. Each panel presented represents two iterations of the PSO algorithm, indicating that the desired global maximal solution was located within 10 iterations per particle. As with other discrete optimization techniques that have been reviewed previously, this simulation shows that the present methodology is also capable of locating global optimal solutions that might not correspond to any of the original points in the black-box objective function. Thus, this approach can find applications in design of experiments (DOE) studies, for example, where the location of such a solution might help to point the user toward the likely combination of optimal conditions to conduct an experiment so that further experimental efforts can be targeted at locations close to this optimal point for actual verifications. Similarly, Figure 3 shows that the PSO algorithm was able to locate the global maximum of a blackbox objective function with N = 102. However, with the larger objective function and using the same number of particles, P = 10, a larger number of iterations was required, as was certainly expected. The global solution was found after about 30 iterations per particle of the PSO algorithm, and it took approximately 2314



ARTICLE

Figure 6. Continuous equivalent of the one-dimensional regression and optimization problem investigated by Coleman and Block.23 The black symbol indicates the global optimal solution found by the PSO algorithm.

another 40 iterations per particle for almost all particles in the swarm to converge at this global maximal point. Figure 4 shows a typical performance chart of the PSO algorithm in solving black-box optimization problems of various sizes that had been formulated as continuous problems using the discrete Fourier series methodology. Here, problem size can be defined based on the number of data points in the black-box function, the number of particles used in the PSO algorithm, the number of iterations needed to locate the global optimal solution, or various combinations of these factors. As might have been expected, the number of iterations of the algorithm required to locate the global maximum of each continuous function increases with problem size and decreases with number of particles used for the search. It can be noted that, because the PSO algorithm is stochastic in nature, such performances are expected to differ quantitatively between simulation runs, and those presented here are illustrative of their general qualitative trends only. The continuous functions generated in this study can be described as fairly noisy in nature because of the fact that the constituents of the corresponding black-box functions were random data points generated by a random number generator. Consequently, the structure of each continuous function was unknown a priori and always nonconvex in nature, criteria that will almost certainly render conventional gradient-based optimization techniques inapplicable. In contrast and as mentioned previously, the stochastic PSO algorithm was observed to be highly successful in achieving global optimization of all objective functions considered here. It can therefore be concluded that the discrete Fourier series methodology coupled with a stochastic optimization technique such as the PSO algorithm is indeed a promising approach for solving black-box optimization problems. Furthermore, although only unconstrained optimization problems were considered in this study, it is anticipated that this methodology can also be extended easily to solving constrained black-box optimization problems. This will require modification of the PSO algorithm to restrict the movements of particles to regions within the domain that satisfy the constraints defined for the optimization problem. Figure 5 shows results of applying the optimization methodology to the derivation of an optimal formulation for the culture

Figure 7. (a) Rosenbrock function in the domains 0 e x e 2 and 0 e y e 2. The color contour scales linearly from blue (0) to red (1000). The function is assumed to be unknown analytically, and white dots indicate positions at which function values are available for black-box optimization. (b) Continuous equivalent of this black-box function derived by applying a discrete double Fourier sine series to this set of about 100 discrete data points.

medium for cultivating the microalga Haematococcus pluvialis. The previous researchers, Fabregas et al.,22 formulated their culture medium for the microalga based on that reported by Nichols24 and defined that as their reference or control. They then repeated their experiments using culture media in which the concentrations of 12 types of nutrients were present at levels 0, 10, or 100 times those in the control. The optimization objective was to find the concentration of each type of nutrient that maximized the cell concentration of the microalga. It can be seen that the discrete sets of experimental data reported by Fabregas et al.22 for the 12 types of nutrients considered in their optimization strategy could be converted into equivalent continuous objective functions by the discrete Fourier series approach, and global optimal solutions could then be located by the PSO algorithm. However, the agreement between the discrete Fourier 2315



ARTICLE

Figure 8. Global optimization by the PSO algorithm on the discrete double Fourier sine series representation of the Rosenbrock function generated from about 100 discrete data points. (a) Initial random distribution of 1000 PSO particles over the domain. Snapshots of the optimization process showing positions of the particles after (b) 200, (c) 400, (d) 600, (e) 800, and (f) 1000 iterations of the PSO algorithm.

series and the corresponding experimental data was observed to be poorer than in the previous cases, and this is most likely due to the fact that the experimental data were collected at irregular intervals of the independent variable. This represents a possible limitation of the present methodology that makes use of a standard discrete Fourier series for transformation of black-box functions into continuous equivalents that inherently assumes that data points in the former exist at regular intervals. With this limitation in mind, comparisons between the optimal formulations derived experimentally by Fabregas et al.22 and those derived numerically in the present study with three other commonly used culture media for cultivating this microalga are presented in Table 1. The optimal concentration value represented in terms of the multiplicative factor relative to the amount present in the culture medium reported by Nichols24 was obtained from Figure 5 for each type of nutrient. The actual

nutrient concentrations as presented in Table 1 were then obtained by multiplying the reported values by the nutrient concentration values of Nichols.24 For a few types of nutrients (MgSO4 3 7H2O, CaCl2 3 2H2O, FeC6H5O7 3 5H2O, Cr2O3, and MnCl2 3 4H2O), Fabregas et al.22 derived optimal concentrations that were higher than those of commonly used culture media such as that due to Nichols,24 namely, BBM and CHU media. However, the results of the present study seem to indicate that the global optimal concentrations of these nutrients might be even higher than those derived experimentally. On the other hand, the converse is true for a second group of nutrients (CoCl2 3 6H2O, CuSO4 3 5H2O, and Na2MoO4 3 2H2O), whereby optimal concentrations were found experimentally to be lower than those in commonly used culture media but the global optimal concentrations obtained numerically were even lower. These results might be illustrative of the inherent restrictions 2316



Figure 9. (a) Rosenbrock function in the domains 0 e x e 2 and 0 e y e 2. The color contour scales linearly from blue (0) to red (1000). The function is assumed to be unknown analytically, and white dots indicate positions at which function values are available for black-box optimization. (b) Continuous equivalent of this black-box function derived by applying a discrete double Fourier sine series to this set of about 400 discrete data points.

associated with discrete experimental optimization strategies in which only a finite number of experiments can be conducted but global optimal conditions do not correspond to any of those actually used. The reformulation of such black-box optimization problems into their equivalent continuous counterparts and subsequent location of the corresponding global optimal solutions will certainly be helpful for experimentalists in planning further experiments to obtain global optimal conditions. In the absence of a priori knowledge of the experimental system, an iterative approach whereby the current optimization technique is used to derive a possible global optimum of the system and the optimal conditions so derived are verified experimentally and in turn used to refine the continuous representation of the experimental system can be employed to locate the actual global optimal point. Such an approach that alternates between experimentation and computation might allow

ARTICLE

experimentalists who are involved in DOE studies to determine global optimal conditions of their experimental systems efficiently, especially in situations where the systems concerned are black boxes. Coleman and Block23 considered an arbitrary one-dimensional regression problem that consisted of seven simulated data points to illustrate the ability of a Bayesian regularized feedforward neural network to predict global optima. They applied a trained neural network model containing one hidden layer with 10 nodes and one output node to approximate the underlying function of the simulated data. Figure 6 shows that the discrete Fourier series method can also be used to obtain an approximate of such an underlying function of the data and the PSO algorithm is able to locate the global maximum of the function. In contrast with the previous case study, the agreement between the discrete Fourier series and the simulated data is much better, as the data in this case are regularly spaced. Figure 7a shows the two-dimensional Rosenbrock function in the domains 0 e x e 2 and 0 e y e 2, where the global minimum point is known to exist at the position (x, y) = (1, 1). Values of the function at about 100 points indicated by the white dots were assumed to be known, whereas the rest of the function was assumed to be unavailable. This set of discrete data points was then treated as part of a black-box function whose global minimum point was to be located. The discrete double Fourier sine series was first applied to this set of discrete data points, and Figure 7b shows that the discrete Fourier series representation obtained was qualitatively similar to the original Rosenbrock function, although quantitative differences can be discerned. The PSO algorithm with 100 particles was then applied to this discrete Fourier series representation of the original function, and Figure 8 shows that the particles were able to converge gradually over about 1000 iterations. The position of the global minimal point reported at the end of 1000 iterations was (x, y) = (1.10, 1.19), which is in close agreement with the actual global minimum of the Rosenbrock function. It is anticipated that closer agreement can be achieved by having a more accurate reconstruction of the Rosenbrock function using the discrete double Fourier sine series and that this can, in turn, be achieved by knowing the function values at a larger number of points on the original function. This is simulated in Figure 9a, where values of the function at about 400 points are assumed to be known. Figure 9b shows that a better representation of the original function is obtained by applying the discrete double Fourier sine series to this set of 400 data points. As with the previous case, Figure 10 shows that the PSO algorithm with 100 particles was able to converge gradually and that the position of the global minimal point reported at the end of 1000 iterations was (x, y) = (0.951, 0.900). As anticipated, a closer agreement with the actual global minimal solution was achieved with a more accurate reconstruction of the original function. In terms of application of such a methodology to solving real black-box optimization problems such as in DOE studies, this implies that more accurate global optimal solutions can be obtained with knowledge of the black-box objective function at larger number of points. There is thus a tradeoff between the amount of resources to be invested to collect information, through either experiments or simulations, about a black-box objective function to be optimized and the accuracy of the global optimal solution that can be achieved through this black-box optimization strategy. The present case study also illustrates the applicability of the methodology for solving higher-dimensional black-box optimization problems. In fact, it can be assumed that there is no limitation to the 2317



ARTICLE

Figure 10. Global optimization by the PSO algorithm on the discrete double Fourier sine series representation of the Rosenbrock function generated from about 400 discrete data points. (a) Initial random distribution of 1000 PSO particles over the domain. Snapshots of the optimization process showing positions of the particles after (b) 200, (c) 400, (d) 600, (e) 800, and (f) 1000 iterations of the PSO algorithm.

number of independent variables of a function that can be regressed with the use of this methodology, as the discrete Fourier series method was derived from the conventional Fourier series and principles associated with the discrete Fourier transform. However, it is also anticipated that the computational time required to evaluate the Fourier coefficients in a high-dimensional discrete Fourier series might become significant as the number of independent variables in the black-box function becomes very large. The application of the black-box optimization methodology developed here to solving the chemical plume tracing problem is considered in the final part of this study. The localization of a source of chemical release by dynamic measurements of the plume it forms is a complex real-world problem because of the stochastic and dynamic nature of the physical environment. It

contains a high degree of freedom and is a highly nonlinear problem for which no straightforward mathematical solution exists. Although there have been several attempts to address the problem of chemical plume tracing by the use of autonomous vehicles equipped with chemical sensors, the equally promising and alternative approach of utilizing a network of stationary sensors has received less attention. We first suppose a scenario whereby a region measuring 1 km 1 km is contaminated with a toxic chemical or pollutant that is released from an unknown location within the region. The chemical forms a gaseous plume that is invisible and odorless. Because of turbulent atmospheric conditions, geographical landscapes, and natural and man-made structures within the region, a complex distribution of the chemical results. Figure 11a shows a randomly generated twodimensional function that was used to simulate such a scenario. 2318



Figure 11. Concentration distribution of a chemical released from an unknown location within a 1 km 1 km region. The color contours represent concentrations of the chemical, with red and blue indicating high and low concentrations, respectively. (a) Network of 10 10 chemical sensors represented by white dots is assumed to have been preinstalled before the chemical release incident. Each sensor is capable of making concentration measurements in real time at its own location and transmitting those measurements to a central computer system for analyses. (b) Concentration distribution of the chemical reconstructed from the sensor measurements using the discrete Fourier series method. Major features of the original contour map have been reconstructed although quantitative differences can also be discerned.

The color contours represent concentrations of the chemical, with red and blue indicating high and low concentrations, respectively. In comparison with Gaussian plume or puff types of dispersion patterns, such a concentration distribution is more likely to arise in a chemical release incident that occurs in an urban setting because of the complex boundary conditions imposed by the buildings and structures, as well as other factors such as street canyon effects. Pockets of chemicals can form at various positions within the region as a result of the chemical being trapped in dead zones around buildings and structures. However, in an actual chemical release incident, such a contour map is almost certainly unavailable for analysis. If it were, the chemical plume tracing problem would be solved instantaneously simply by a visual

ARTICLE

Figure 12. Concentration distribution for a more complex chemical release scenario. (a) Correspondingly, a larger number of sensors represented by white dots is assumed to have been preinstalled. As with the previous scenario, the sensors report point measurements of the continuous concentration distribution. (b) Treating the set of discrete values from the sensor network as a discrete function whose continuous equivalent is to be derived, the discrete Fourier series method is applied, and a much more accurate representation of the original concentration distribution is obtained.

search for the location of highest chemical concentration within the region. In principle, the entire contour map can be constructed by performing a very large number of measurements at different positions in the region, but this is almost never possible in practice in an emergency situation. Instead, we assume that a network of stationary sensors had been preinstalled within the region prior to the chemical release incident. Each sensor, represented by a white dot in Figure 11a, is capable of making concentration measurements in real time at its own location and transmitting those measurements to a central computer system for analysis. With a network containing 10 10 sensors, this is equivalent to taking 100 point measurements of the continuous concentration distribution, and we then ask whether it is possible to achieve an accurate reconstruction of the entire contour map with these discrete data points. Treating this set of discrete 2319



ARTICLE

Figure 13. Chemical plume tracing by the PSO algorithm performed on a contour map that was reconstructed by the discrete Fourier series method. (a) Initial random distribution of 10 PSO particles over the domain of the reconstructed contour map. Snapshots of the search process showing positions of the particles after (b) 200, (c) 400, (d) 600, (e) 800, and (f) 1000 iterations of the PSO algorithm. All particles converged to a common location on the contour map after 1000 iterations. This final positions of all particles represent the global maximal point of both the original and reconstructed contour maps. This indicates that the PSO algorithm has been successful in achieving global optimization of these functions and that localization of the source of chemical release has been successfully accomplished.

values as a discrete function whose continuous equivalent was to be derived, the discrete double Fourier sine series was applied to this discrete function. Figure 11b shows the resulting contour map that was reconstructed from this set of discrete values. On comparison between the two panels, it can be seen that the major features of the original contour map in terms of approximate regions where high and low concentrations of chemical occur was reconstructed, although quantitative differences in the actual concentration distributions can be discerned. The reconstructed contour map might be adequate for chemical plume tracing if it allows a good approximate of the global maximum point of the original concentration distribution to be located. Otherwise, a larger number of discrete values of the

concentration distribution might be required for a more accurate reconstruction. Figure 12a shows another arbitrarily generated concentration distribution that contains a larger number of regions of high and low concentrations. Correspondingly, as depicted by the white dots shown on the figure, it was assumed that a larger number (20 20) of sensors were present in the region for concentration measurements. After the reconstruction process had been performed by applying the discrete double Fourier sine series to this set of discrete values, Figure 12b shows that a much more accurate representation of the original concentration distribution was obtained. In fact, the original and reconstructed contour maps are almost indistinguishable from each other. It can then be 2320


Industrial & Engineering Chemistry Research expected that a good approximate of the global maximum of the original concentration distribution can be obtained by applying the PSO algorithm to the reconstructed contour map, as presented next. It can be noted at this point that the reconstruction process applied here did not require any prior assumptions about the characteristics of spread of the chemical plume or a priori knowledge of wind profiles, wind speeds, or wind directions within the region. The sensors were assumed to be capable of measuring the concentration of the chemical or pollutant at their respective positions only, and no other information was utilized. With an accurate reconstruction of the chemical concentration distribution, the chemical plume tracing problem is essentially solved, and what remains is simply determination of the location of highest concentration within the region. This necessarily assumes that the source is continuing to emit the chemical or pollutant during the plume tracing process. The contrasting scenario whereby the source is empty and no longer emitting is considered to be outside the scope of the present study. The chemical plume tracing problem can then be treated as an optimization problem in which the objective function to be maximized is the discrete Fourier series representation of the concentration distribution. The PSO algorithm can then be applied to locate the global optimum of this objective function. Figure 13a shows positions of 10 PSO particles randomly distributed over the domain of the reconstructed contour map from Figure 12b at the start of the optimization process. The positions of these particles were updated iteratively, and after 1000 iterations, Figure 13f shows that all particles converged to a common location on the contour map. It was verified manually that this final position of all particles was the global maximum point of both the original and reconstructed contour maps. This indicates that the PSO algorithm was successful in achieving global optimization of these functions. More importantly, this also means that localization of the source of chemical release was successfully achieved. To our knowledge, this is the first report of solving the chemical plume tracing problem by the approach of global optimization of a reconstructed concentration distribution of the chemical. Several such chemical plume tracing simulations were conducted in the course of this study, and the methodology presented here, comprising the discrete Fourier series method and PSO, was found to be successful in solving the chemical plume tracing problem in each case. The program codes for the algorithm were executed on a typical personal computer, and computation times were always negligible in the sense that results were always generated almost instantaneously. It can also be noted at this point that the only assumption required for the methodology to be successful is that the position of the source of chemical release to be located is also the point of maximum chemical concentration within the entire region. This assumption is expected to hold as long as chemical is still being released from its source during the search process, which allows the chemical plume tracing problem to be transformed into an optimization problem as was done here. The contrasting problem where the source is depleted before the search begins can be referred to as a chemical puff tracing problem. Although not shown here for the sake of brevity, direct application of the methodology to such a scenario would lead to particles in the PSO algorithm “chasing” after the most concentrated puff rather than converging to the position of the depleted source. For both chemical plume and puff scenarios, where a highly lethal but invisible cloud might be sweeping across a country rapidly, the ability to have a visual image of the entire dispersion process can be very beneficial for warning purposes and evacuation planning.

ARTICLE

This is especially so for residents who are situated in the path of the toxic cloud. In the context of such chemical release scenarios, it is believed that the algorithm developed in this study will be able to play an instrumental role in the national defense for any country in the world that is subjected to such threats.

’ CONCLUSIONS A new methodology for solving black-box optimization problems by the continuous approach has been developed in this study. A discrete Fourier series method was derived from the conventional Fourier series formulation and principles associated with the discrete Fourier transform and used for reformulation of black-box objective functions as continuous functions. A stochastic global optimization technique known as particle swarm optimization (PSO) was then applied to locate the global optimal solutions of the continuous functions derived. The structure of each continuous function was unknown a priori and always nonconvex in nature, and the stochastic PSO algorithm was observed to be highly successful in achieving global optimization of all such objective functions considered in this study. The methodology developed here was applied successfully to optimization of the formulation for a culture medium used for cultivating a freshwater microalga, Haematococcus pluvialis, as well as a one-dimensional regression problem. Subsequently, the methodology was extended to solving a two-dimensional black-box optimization problem that was simulated based on the Rosenbrock function. The possibility of applying the same approach for solving the chemical plume tracing problem in the area of chemical defense was also successfully demonstrated. The discrete Fourier series method coupled to the PSO algorithm is thus a promising methodology for solving black-box optimization problems by a continuous approach. Possible avenues deemed worthy of pursuit beyond the current work include extending the above methodology to higher dimensions for solving optimization problems with more than two independent variables. In particular, the possibility of applying a three- or higher-dimensional discrete Fourier series method for generating an accurate response hypersurface for design of experiments (DOE) studies through the response surface methodology (RSM) approach can be explored.

’ AUTHOR INFORMATION Corresponding Author

*Tel.: (65) 6516 2542. Fax: (65) 6779 1936. E-mail: chelwce@ nus.edu.sg.

’ ACKNOWLEDGMENT This study was supported by the Ministry of Defence (MINDEF) of Singapore and the National University of Singapore (NUS) under the MINDEF-NUS Joint Applied R&D Co-Operation Programme (JPP) under Grants R-279-000-302-232 and R-279-000-302-646. The award of a Lee Kuan Yew Postdoctoral Fellowship by the National University of Singapore is gratefully acknowledged. ’ REFERENCES (1) Droste, S.; Jansen, T.; Wegener, I. Optimization with randomized search heuristics: The (A)NFL theorem, realistic scenarios, and difficult functions. Theor. Comput. Sci. 2002, 287, 131. (2) Droste, S.; Jansen, T.; Wegener, I. Upper and lower bounds for randomized search heuristics in black-box optimization. Theory Comput. Syst. 2006, 39, 525. 2321



ARTICLE

(3) Audet, C.; Orban, D. Finding optimal algorithmic parameters using derivative-free optimization. SIAM J. Optim. 2006, 17, 642. (4) Driessen, L.; Brekelmans, R.; Hamers, H.; den Hertog, D. On D-optimality based trust regions for black-box optimization problems. Struct. Multidisc. Optim. 2006, 31, 40. (5) Ramani, S.; Blu, T.; Unser, M. Monte-Carlo Sure: A black-box optimization of regularization parameters for general denoising algorithms. IEEE Trans. Image Process. 2008, 17, 1540. (6) Andradottir, S. A method for discrete stochastic optimization. Manage. Sci. 1995, 41, 1946. (7) Gong, W.; Ho, Y.; Zhai, W. Stochastic comparison algorithm for discrete optimization with estimation. SIAM J. Optim. 1999, 10, 384. (8) Martinelli, F. Stochastic comparison algorithm for discrete optimization with estimation of time-varying objective functions. J, Optim. Theory Appl. 1999, 103, 137. (9) Kleywegt, A. J.; Shapiro, A.; Homem-de Mello, T. The sample average approximation method for stochastic discrete optimization. SIAM J. Optim. 2001, 12, 479. (10) Alrefaei, M. H.; Andradottir, S. A modification of the stochastic ruler method for discrete stochastic optimization. Eur. J. Oper. Res. 2001, 133, 160. (11) Alrefaei, M. H.; Andradottir, S. Discrete stochastic optimization using variants of the stochastic ruler method. Nav. Res. Logist. 2005, 52, 344. (12) Pardalos, P. M.; Romeijin, H. E.; Tuy, H. Recent developments and trends in global optimization. J. Comput. Appl. Math. 2000, 124, 209. (13) Goyal, V.; Ierapetritou, M. G. Stochastic MINLP optimization using simplicial approximation. Comput. Chem. Eng. 2007, 31, 1081. (14) Aytug, H.; Sayin, S. Using support vector machines to learn the efficient set in multiple objective discrete optimization. Eur. J. Oper. Res. 2009, 193, 510. (15) Brigham, E. O. The Fast Fourier Transform; Prentice-Hall: Englewood Cliffs, NJ, 1974. (16) Engelbrecht, A. P. Computational Intelligence: An Introduction; John Wiley & Sons Ltd.: Chichester, U.K., 2002. (17) Kennedy, J.; Eberhart, R. C. Particle swarm optimization. Proc. IEEE Int. Conf. Neural Networks 1995, 4, 1942. (18) Kennedy, J. The behavior of particles. In Proceedings of the 7th International Conference on Evolutionary Programming; Porto, V. W., Saravanan, N., Waagen, D., Eds.; Springer: New York, 1998; pp 581-589. (19) van den Bergh, F. An analysis of particle swarm optimizers. Ph.D. Thesis, University of Pretoria, Pretoria, South Africa, 2002. (20) Li, C.; Zhu, Q.; Geng., Z. Multi-objective particle swarm optimization hybrid algorithm: An application on industrial cracking furnace. Ind. Eng. Chem. Res. 2007, 46, 3602. (21) Tang, L. X.; Yan, P. Particle swarm optimization algorithm for a campaign planning problem in process industries. Ind. Eng. Chem. Res. 2008, 47, 8775. (22) Fabregas, J.; Dominguez, A.; Regueiro, M.; Maseda, A.; Otero, A. Optimization of culture medium for the continuous cultivation of the microalga Haematococcus pluvialis. Appl. Microbiol. Biotechnol. 2000, 53, 530. (23) Coleman, M. C.; Block, D. E. Nonlinear experimental design using Bayesian regularized neural networks. AIChE J. 2007, 53, 1496. (24) Nichols, H. W. Growth media-freshwater. In Handbook of Phycological Methods. Culture Methods and Measurements; Stein, J. R., Ed.; Cambridge University Press: London, 1973; pp 7-24.

2322


Application of Particle Swarm Optimization to ... - ACS Publications

Recommend Documents