A New Stochastic Optimization Algorithm based on a ... - CiteSeerX

6th World Congresses of Structural and Multidisciplinary Optimization Rio de Janeiro, 30 May - 03 June 2005, Brazil

A New Stochastic Optimization Algorithm based on a Particle Collision Metaheuristic Wagner F. Sacco1, Cassiano R.E. de Oliveira2 1

Nuclear & Radiological Engineering Program, Georgia Institute of Technology, Atlanta, GA, USA, [email protected] 2 Nuclear & Radiological Engineering Program, Georgia Institute of Technology, Atlanta, GA, USA, [email protected]

In many complex optimization problems, the so-called stochastic optimization methods as genetic algorithms and simulated annealing have been applied successfully. However, these methods have drawbacks: simulated annealing is too sensitive to the choice of parameters and the canonical genetic algorithm is liable to suboptimal convergence. With the purpose of overcoming these drawbacks, a new stochastic optimization algorithm is introduced in this article. It is loosely inspired by nuclear collision reactions, particularly scattering (where an incident neutron is scattered by a collision with a target nucleus) and absorption (where the incident neutron is absorbed by the target nucleus). This algorithm was named “Particle Collision Algorithm” (PCA). PCA resembles in its structure simulated annealing. A basic difference is that, besides the number of iterations, it does not rely on user-defined parameters. First, the new algorithm is validated using three well known test functions. Then, PCA is applied to an engineering design optimization problem: a nuclear core design optimization. This problem consists in adjusting several reactor cell parameters, such as dimensions, enrichment and materials, in order to minimize the average peak-factor in a 3-enrichment-zone reactor, considering operational restrictions. Minimizing the average peak-factor is important for safety reasons, avoiding the risk of a fuel rod meltdown, and for economical reasons, so that the cycle length is maximized. The results are compared to those obtained in previous efforts using the canonical genetic algorithm and the niching genetic algorithm. PCA performs better than both methods, showing its potential for other applications. Keywords: Stochastic Optimization, Metaheuristics, Design Optimization. 1. Introduction Stochastic optimization methods based on the simulated annealing paradigm [1] have been actively developed in the last 30 years and successfully applied to many complex optimization problems. These methods though very powerful are not free from practical drawbacks, the main one being that performance is too sensitive to the choice of free parameters, such as, for example, the annealing schedule and initial temperature [2]. Ideally, an optimization algorithm should not rely on user-defined or problem-dependent parameters and should not converge to a suboptimal solution. Given the appropriate annealing schedule, simulated annealing is guaranteed to converge to the global optimum [3], as it is a Metropolis-based algorithm [4], where a worse solution can be accepted with a certain probability. However, its rate of convergence is strongly dependent on the user-specified initial parameters. So, the best outcome would be the development of an algorithm similar to simulated annealing, but without the burden of parameter specification. One such algorithm, which does not rely on user-supplied parameters to perform the optimality search, is presented here. Named the “Particle Collision Algorithm” (PCA), the algorithm is loosely inspired by the physics of nuclear particle collision reactions [5], particularly scattering (where an incident particle is scattered by a target nucleus) and absorption (where the incident particle is absorbed by the target nucleus). Thus, a particle that hits a high-fitness “nucleus” would be “absorbed” and would explore the boundaries. On the other hand, a particle that hits a low-fitness region would be scattered to another region. This permits us to simulate the exploration of the search space and the exploitation of the most promising areas of the fitness landscape through successive scattering and absorption collision events. The PCA is applied to a nuclear reactor design optimization problem that was previously attacked using the standard genetic algorithm [6] and the niching genetic algorithm [7]. The three algorithms are submitted to the same computational effort and their results are compared. The remainder of the paper is organized as follows. In the next section the details of the implementation of the method are presented. In section 3, the new algorithm is validated. In section 4 the optimization problem is outlined. The next section presents the system implementation and the comparative tests. Finally, in Section 6 conclusions are presented along with suggestions of future improvements to the algorithm. 2. The PCA The PCA resembles in its structure that of simulated annealing: first an initial configuration is chosen; then there is a modification of the old configuration into a new one. The qualities of the two configurations are compared. A decision then is made on whether the new configuration is “acceptable”. If it is, it serves as the old configuration for the next step. If it is not acceptable, the algorithm proceeds with a new change of the old configuration. PCA can also be considered a Metropolis algorithm, as a trial solution can be accepted with a certain probability. This acceptance may avoid the convergence to local optima. The pseudo code description of the PCA is shown in Figure 1 on its default version for maximization problems.

Generate an initial solution Old_Config For n = 0 to # of iterations Generate a stochastic perturbation of the solution If Fitness(New_Config) > Fitness(Old_Config) Old_Config := New_Config Exploration () Else Scattering () End If End For Exploration () For n = 0 to # of iterations Generate a small stochastic perturbation of the solution If Fitness(New_Config) > Fitness(Old_Config) Old_Config := New_Config End If End For return Scattering () Fitness ( New _ Config ) Best Fitness If pscattering > random (0, 1) Old_Config := random solution Else Exploration (); End if return pscattering = 1 −

Figure 1. PCA’s pseudo code. The “stochastic perturbation” mentioned in the beginning of the loop are random variations in each variable’s values within their ranges. If the quality or fitness of the new configuration is better than the fitness of the old configuration, then the “particle” is “absorbed”, there is an exploration of the boundaries searching for an even better solution. Function “Exploration()” performs this local search, generating a small stochastic perturbation of the solution inside a loop. In PCA’s current version, it is a one-hundred-iteration loop. The “small stochastic perturbation” is similar to the previous stochastic perturbation, but each variable’s new value is kept within the boundaries of the original value. Otherwise, if the quality of the new configuration is worse than the old configuration’s, the “particle” is “scattered”. The scattering probability (pscattering) is inversely proportional to its quality. A low-fitness particle will have a greater scattering probability. In a process similar to Monte Carlo’s “Russian Roulette” [8], the configuration is “scattered” (replaced by a random configuration) or, following Metropolis, survives, with its boundaries explored (“else” branch of the function). Note once more that PCA does not need user-defined parameters as the annealing schedule in SA. 3. Validation As the search space of the reactor design optimization problem is not known beforehand, prior to applying PCA to this problem it was necessary to validate it by application to commonly employed test functions and comparison to the results achieved by the widely employed genetic algorithm, which is briefly outlined in the next subsection. 3.1. Genetic Algorithms Genetic Algorithms [9] are searching methods based upon the biological principles of natural selection and survival of the fittest introduced by Charles Darwin in his seminal work “The Origin of Species” (1859). They were rigorously introduced by Holland [10]. GAs consist of a population of individuals that are possible solutions and each one of these individuals receives a reward, known as “fitness”, that quantifies its suitability to solve the problem. In ordinary applications, fitness is simply the objective function. Individuals with better than average fitnesses receive greater opportunities to cross. On the other hand, low fitness individuals will have less chance to reproduce until they are extinguished. Consequently, the good features of the best individuals are disseminated over the generations. In other words, the most promising areas of the search space are explored, making the GA converge to the optimal or near optimal solution. Genetic Algorithms have proven to be efficient in a great variety of domains, as the population of candidate solutions converge to a single optimum, in a phenomenon known as genetic drift [9]. Many populational diversity mechanisms, called niching methods [11], have been proposed to force the GA to maintain a heterogeneous population throughout the evolutionary process. The niching method which will be referred herein is Fuzzy Clearing [7], where the individuals are clustered using an algorithm called

“Fuzzy Clustering Means” [12] and the individual with best fitness (dominant) is determined for each cluster or class, using Fuzzy Logic [13] terminology. Following that, the dominant’s fitness is preserved and all the others individuals have their fitnesses zeroed (in the case of a maximization problem). Another way of preserving populational diversity that will be soon referred is using the Island Genetic Algorithm [14]. In this model, the population is divided into multiple populations that evolve isolated from each other most of the time, but exchange individuals occasionally with their neighbors (migration). 3.2. Test Functions Easom’s function [15] is a unimodal test function where the global minimum region has a small area relatively to the search space, making it a great challenge for optimization algorithms. Easom is defined by Eq.(1).

z(x, y)=−cos(x)⋅cos(y)⋅ exp{−1⋅[(x −π)2 +(y −π)2]}, −100 ≤ x, y ≤100,

(1)

with global minimum in (π, π, -1). Shekel’s Foxholes, introduced by Shekel [16] and adapted for maximization by De Jong [17], is a two-dimensional function with 25 peaks all with different heights, ranging from 476.191 to 499.002. The global optimum is located at (-32, -32). Shekel’s Foxholes is defined by Eq.(2):

z(x, y)=500−

1 24

1 0.002+∑ 6 6 i=0 1+i +(x −a(i)) +(y −b(i))

, −65.536 ≤ x, y ≤ 65.536,

(2)

where a(i) = 16[(imod5) – 2] and b(i) = 16[(i/5) – 2]. The last test function, Rosenbrock's valley, is a classic optimization problem [18]]. The global optimum is inside a long, narrow, parabolic shaped flat valley. To find the valley is trivial, however convergence to the global optimum is difficult and hence this problem has been repeatedly used to test the performance of optimization algorithms. The global optimum of the function, defined by Eq.(3), is (1, 1, 0).

z ( x, y ) =100 ⋅ ( y − x 2 ) 2 + (1 − x) 2 , − 2.048 ≤ x, y ≤ 2.048

(3)

3.3. Test Functions The standard genetic algorithm (SGA) used in the validation tests was GENESIS [19], which employs double-point crossover [9], Stochastic Universal Sampling as selection scheme [20] and elitism [9]. The following parameters were used: 100,000 function evaluations (population size = 100, number of generations = 1000), crossover rate = 0.6, mutation rate = 0.005, generation gap = 0.1. PCA and GDA were set up for 100,000 function evaluations and the rain speed in the latter followed equation (1). This high number of function evaluations is due to the high degree of difficulty of these functions. Table 1 displays the results obtained by each algorithm (the best values are bold-faced).

Experiment #1 #2 #3 #4 #5 #6 #7 #8 #9 #10 Average

Table 1. Results for test functions from the literature Easom’s Function Shekel’s Foxholes SGA PCA SGA PCA 0.000 1.000 499.002 499.002 1.000 1.000 499.002 499.002 0.000 1.000 499.002 499.002 0.000 1.000 498.536 499.002 0.360 1.000 499.002 499.002 0.000 1.000 499.002 499.002 0.000 1.000 499.002 499.002 0.000 1.000 499.002 499.002 0.300 1.000 499.002 499.002 0.000 1.000 498.589 499.002 0.166 1.000 498.914 499.002

Rosenbrock’s Valley SGA PCA 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

The results show that PCA not only reached the optimal values in all test functions but also outperformed the SGA, achieving 100 % effectiveness. This is particularly true in the case of Easom’s function. Because of the extremely narrow valley, the standard genetic algorithm is liable to drift before reaching this region. Even in a relatively easier problem like Shekel’s Foxholes, the SGA got stuck in local optima.

4. Problem Description The nuclear engineering optimization problem used to test PCA’s performance is outlined in this section. Its detailed description is beyond the scope of this article, but the interested reader should refer to [6]. Briefly stated, the problem consists in designing a reactor core which minimizes a parameter called average power peak factor without violating certain operational constraints. The power peak factor is the ratio between the maximum and the average power that occur in the fuel rods. The closer this factor is to 1, the more efficient and the safer is the reactor. A difference of 10-3 in this parameter is significant. The reactor considered here is a cylindrical, three-enrichment- zone Pressurized Water Reactor (PWR), with a typical cell composed by the fuel, the fuel cladding and the moderator (light water). The moderator is used to slow down the particles without absorbing them, as the nuclei from the reactor fuel are particularly liable to fission if bombarded by slow particles. Water is a good choice for moderator, as it also acts as reactor coolant. Besides the dimensions, there is also an array of materials that can be considered in the design. Figure 2 illustrates such reactor and Table 2 shows the design parameters. Fuel

R2

R1

Cladding

R3 Re

∆c

h

Rf Moderator

(a)

(b)

Figure 2. (a) The reactor and (b) its typical cell. Parameter Fuel Radius (cm) Cladding Thickness (cm) Moderator Thickness (cm) Enrichment of Zone 1 (%) Enrichment of Zone 2 (%) Enrichment of Zone 3 (%) Fuel Material Cladding Material

Table 2. Parameters range. Symbol Range Rf 0.508 to 1.270 0.025 to 0.254 ∆c Re 0.025 to 0.762 E1 2.0 to 5.0 E2 2.0 to 5.0 E3 2.0 to 5.0 Mf {U-metal or UO2} Mc {Zircaloy-2, Aluminum or Stainless -304}

5. Implementation and Results 5.1. Implementation The optimization algorithm is linked to the Reactor Physics code HAMMER [21] for the cell and diffusion equations calculations. The GA or PCA sends to the code by an interface a solution (the reactor’s dimensions and materials) and receives the fitness. The fitness function was developed in such a way that, if all constraints are satisfied, it has the value of the average peak factor. Otherwise, it is penalized proportionally to the discrepancy on the constraint. In our tests, both GAs had the same setup and number of iterations as in section 3. In the case of the NGA, the number of clusters or classes selected for FCM was ten, the same used in [7]. 5.2. Results Table 3 shows the results obtained by PCA in comparison to those obtained by the standard genetic algorithm and by the niching genetic algorithm.

Table 3. Comparative results for the reactor design optimization problem for 50,000 and 100,000 iterations. SGA NGA PCA Experiment 5·104 iter. 1·105 iter. 5·104 iter. 1·105 iter. 5·104 iter. 1·105 iter. #1 1.3185 1.3185 1.2916 1.2916 1.2849 1.2827 #2 1.3116 1.3116 1.3069 1.3069 1.2876 1.2876 #3 1.3300 1.3300 1.3003 1.3003 1.2964 1.2964 #4 1.3294 1.3294 1.2874 1.2844 1.2953 1.2874 #5 1.3595 1.3595 1.2956 1.2895 1.2829 1.2829 #6 1.3562 1.3562 1.3014 1.3014 1.2791 1.2791 #7 1.3372 1.3372 1.3190 1.2872 1.2975 1.2975 #8 1.3523 1.3523 1.3075 1.3050 1.2865 1.2865 #9 1.3614 1.3614 1.2974 1.2959 1.3010 1.2908 #10 1.3467 1.3467 1.3077 1.3077 1.2852 1.2845 Average 1.3402 1.3402 1.3015 1.2970 1.2896 1.2875 Std. Dev. 0.0175 0.0175 0.0093 0.0085 0.0073 0.0059 We can see that PCA obtained the best overall result, the best average and the lowest standard deviation. Also, PCA’s best result and average at 50,000 iterations are better than NGA’s at 100,000. All SGA’s executions drifted before 50,000 executions. PCA performs well even compared to the parallel Island Genetic Algorithm (IGA), which is the benchmark for this problem [22]. With the same computational effort as PCA’s (4 islands of 50 individuals, 500 generations or 100,000 iterations), IGA obtained 1.2784 as best result with an average of 1.2992, while PCA’s best was 1.2791 with an average of 1.2875. Table 4 shows the best configurations obtained by the SGA with 300 individuals until convergence in [6], by IGA with 400 individuals and 500 generations in [22], by the NGA with 100 individuals and 500 generations in [7] and by PCA in 100,000 iterations, when applied exactly to the same problem. Table 4. Comparison with previously published results. SGA [6] IGA [22] NGA [7] Objectives and Constraints

Fitness

1.310

1.278

1.287

1.279

Minimum average peak factor

1.310

1.278

1.287

1.279

8.02 x 10-5

8.08 x 10-5

8.04 x 10-5

8.06 x 10-5

Average Flux

Parameters

PCA

keff

1.000

0.991

1.000

0.991

Rf (cm)

0.5621

0.7661

0.5441

0.5497

∆r (cm)

0.1770

0.1857

0.1064

0.1450

∆m (cm)

0.6581

0.7569

0.5997

0.6111

E1 (%)

2.756

2.6850

2.5906

2.7953

E2 (%)

4.032

2.8268

2.7559

2.9469

E3 (%)

4.457

4.8819

4.6220

5.0000

Mf

U-metal

U-metal

U-metal

U-metal

Mc

Stainless-304

Stainless-304

Stainless-304

Stainless-304

6. Conclusions The comparative performances PCA with the genetic algorithm show that the first algorithm is quite promising and should be applied to other engineering optimization problems. The great advantage of PCA in relation to other optimization algorithms like the genetic algorithm or simulated annealing is that, other than the number of iterations, it does not require any additional parameters. Furthermore, it can be applied to continuous or discrete optimization problems just changing the perturbation function. In genetic algorithms it is necessary to apply special operators for discrete optimization problems [9]. Last but not least, it is extremely easy to implement. The algorithm presented here is in its early stages. As further development, the local search mechanism could be improved, as a more intelligent mechanism could lead to significant gains in computational cost. Another improvement would be a populational algorithm where the particles would interact to obtain better solutions, in a fashion similar to Particle Swarm Optimization [23]. Acknowledgement Wagner F. Sacco is supported by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico, Ministry of Science & Technology, Brazil) under postdoctoral grant 200223/2004-6.

References 1. S. Kirkpatrick, C.D. Gelatt, M.P. Vecchi. Optimization by Simulated Annealing. Science, 1983, 220(4598): 671-680 2. Carter J N. Topology Genetic Algorithms for Incore Fuel Management and other recent developments in optimization. Proc. of APCAM 93, 1993, Sydney, 917-921 Advances in Nuclear Science and Technology 25, Plenum Press, 113-154. 3. Aarts E and Korst J. Simulated Annealing and Boltzmann Machines. New York: John Wiley and Sons, 1989 4. N. Metropolis, A.W. Rosenbluth, M.N. Rosebluth, A.H. Teller, E. Teller. Equations of State Calculations by Fast Computing Machines. J. Chem. Phys., 1953, 21: 1087-1092. 5. Duderstadt J J and Hamilton L J. Nuclear Reactor Analysis. New York: John Wiley and Sons, 1976 6. C.M.N.A. Pereira, R. Schirru, A.S. Martinez. Basic Investigations Related to Genetic Algorithms in Core Designs. Annals of Nuclear Energy, 1999, 26(3): 173-193 7. W.F. Sacco, M.D. Machado, C.M.N.A. Pereira, R. Schirru. The fuzzy clearing approach for a niching genetic algorithm applied to a nuclear reactor core design optimization problem. Annals of Nuclear Energy, 2004, 31(1): 55-69 8. Duderstadt J J and Martin W R. Transport Theory. New York: John Wiley and Sons, 1979 9. Goldberg D E. Genetic Algorithms in Search, Optimization & Machine Learning. Reading: Addison-Wesley, 1989 10.Holland J H. Adaptation in Natural and Artificial Systems. Ann Arbor: University of Michigan Press, 1975 11.Mahfoud S W. Niching Methods for Genetic Algorithms. PhD Thesis, 1995, Illinois Genetic Algorithm Laboratory, University of Illinois at Urbana- Champaign, Urbana, IL. 12.Bezdek J C. Pattern Recognition with Fuzzy Objective Function Algorithms. New York: Plenum Press, 1981 13.L.A. Zadeh. Fuzzy Sets. Information and Control, 1965, 8: 338-352 14.Cantu-Paz E. Efficient and Accurate Parallel Genetic Algorithms. Boston: Kluwer Academic Publishers, 2000 15.Easom E E. A survey of global optimization techniques. M. Eng. Thesis, 1990, University of Louisville, Louisville, KY. 16.Shekel J. Test functions for multimodal search techniques. Proc. of the 5th Princeton Conference on Information Science and Systems, 1971, Princeton, 354-359 17.De Jong K A. An analysis of the behavior of a class of genetic adaptive systems. PhD Thesis, 1975, University of Michigan, Ann Arbor, MI. 18.H.H. Rosenbrock. An Automatic Method for Finding the Greatest or Least Value of a Function. Computer Journal, 1960, 3: 175184. 19.Grefenstette J J. A User's Guide to GENESIS version 5.0. Nashville: Vanderbilt University, 1990 20.Baker J E. Reducing Bias and Inefficiency in the Selection Algorithm. Proc. of the Second International Conference on Genetic Algorithms and their Application, 1987, Hillsdale, NJ, 14-21 21.Suich J E and Honec H C. A User's Guide to GENESIS version 5.0. Aiken: Savannah River Laboratory, 1967 22.C.M.N.A. Pereira, C.M.F. Lapa. Coarse-grained parallel genetic algorithm applied to a nuclear reactor core design optimization problem. Annals of Nuclear Energy, 2003, 30(5): 555-565 23.Kennedy J and Eberhart R. Particle Swarm Optimization. Proc. of the 1995 IEEE International Conference on Neural Networks, 1995, Perth, Western Australia, 1942-1948