- Email: [email protected]
Elite genetic algorithms with adaptive mutations for solving continuous optimization problems – application to modeling of the optical constants of solids
15 May 1998
Optics Communications 151 Ž1998. 147–159
Full length article
Elite genetic algorithms with adaptive mutations for solving continuous optimization problems – application to modeling of the optical constants of solids Aleksandra B. Djurisic ˇ´ Department of EEE, UniÕersity of Hong Kong, Pokfulam Road, Hong Kong, China Received 17 September 1997; revised 6 January 1998; accepted 19 February 1998
Abstract The elite genetic algorithm with adaptive mutations is proposed as a tool for solving continuous optimization problems. The new algorithm and the corresponding classical genetic algorithm were severely tested on three families of multiminima test functions for 20, 50 and 100 variables. All performed tests proved that the introduced adaptive mutations significantly improve the ability of the algorithm to find a global minimum. After verifying the superiority of the proposed algorithm over the classical one on the test functions, an elite genetic algorithm with adaptive mutations was applied for solving the model parameter determination problem for modeling optical constants of the following metals: beryllium, chromium, nickel and palladium. Good agreement between calculated and experimental data was obtained for all four metals. q 1998 Elsevier Science B.V. All rights reserved. PACS: 02.70.–c; 78.20.Ci; 78.66.–w Keywords: Genetic algorithms; Model parameter estimation; Optical constants
1. Introduction Genetic algorithms ŽGAs. w1x are stochastic global search methods that mimic the concept of natural evolution. Due to the nature of the algorithm, their successful application was mostly restricted to optimization problems whose solution can be conveniently represented in binary form. However, there is a rising interest in applying genetic algorithms to continuous optimization problems, especially since there is no need for initial estimates, which is an important advantage of GAs over other stochastic search methods such as, for instance, simulated annealing w2x. For that reason, various modifications of original GAs have been reported w3–7x. GAs search for optimal solution employing mechanisms of natural evolution: selection, mating and mutation, which are applied to the set of possible problem solutions, called population. Each element of the population is represented by a vector of variables, called string or chromosome. Each element of the chromosome is, by analogy,
termed a gene. Chromosomes are characterized by their ‘‘fitness’’ – performance with respect to some objective function. Highly fit individuals have a high probability of being selected for reproduction Žthey may survive or give offsprings in the next generation.. Strings with lower fitness have correspondingly lower probability for transferring their genes to the following generations. This is illustrated in Fig. 1, where shade of gray represents fitness – the lighter shade corresponds to the better gene. Many variations of the simple GA have been developed w8x, differing in chromosome representation, selection, reproduction and mutation, existence of fitness scaling, etc. It was shown that floating-point number representation, instead of binary string representation common for classical GAs, is more convenient for continuous optimization applications w5,6,9x. In such a manner, by representing variable values with real numbers, the length of the chromosome is equal to the number of variables, thus being much smaller than in the case of binary coding. For instance, for coding real numbers p - 2 5 with four deci-
0030-4018r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII S 0 0 3 0 - 4 0 1 8 Ž 9 8 . 0 0 1 2 2 - 9
148
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
Fig. 1. Illustration of GA: shade of gray represents fitness: the lighter shade correspods to a more fitted gene.
mal digits we need 21 genes, i.e. the length of chromosome is about twenty times larger. Also, an important feature of floating-point representation is that the variable values cannot be altered or destroyed during crossover operation, while in the case of binary coding such undesired changes of variable values can lead to loss or deterioration of genetic information. More important, conversion of decimal numbers into binary ones and vice versa is avoided by using floating-point representation. The important advantage of real-coded GAs is the absence of Hamming cliff problem inherent to binary coded GAs w1,9,10x. Such problems could be overcome by the use of Gray coding, but that leads to higher nonlinearities with respect to recombination. Also, using real-coding reduces the dimensionality of the problem, thus reducing the opportunity for deception, but the possibility of new obstacles to convergence is introduced. The fact that realcoded GAs can in certain cases be blocked from further progress has been recognized and discussed w9x. Recently, there has been much work in modification of real-coded GAs in order to make them as successful in solving continuous optimization problems as their binary counterparts are in solving discrete optimization problems w3,4,10,11x. However, in solving certain problems such as, for example, parameter estimation problem, even GAs with floating-point representation and appropriate selection, crossover and mutation mechanisms often fail to give satisfactory results, especially for a large number of variables. These shortcomings are manifested mainly due to the discrete sampling of the solution space, while variables have continuous nature. GAs are capable of roughly locat-
ing a global optimum, but a huge number of chromosomes in the population is required for any refinements. Recently, genetic algorithms with parameter space size adjustment have been proposed for model parameter determination w7x. The main idea is narrowing of the solution space by improving the guess of the global minimum position. This is achieved by introducing one additional loop. The initial population is generated, classical GA is performed, and the boundaries for each parameter are narrowed towards the obtained average value for that parameter in the inner loop final population. The new initial population is generated in the new boundaries and GA is performed again. The outer loop is executed for a specified number of iterations n max , which means that the time required for performing this algorithm is n max times larger than that for the corresponding classical GA. In this paper, an elite genetic algorithm with adaptive mutations ŽEGAAM. is proposed, which can successfully solve the presented problems, without any observable increase of CPU time requirements. In other words, the execution time of EGAAM and GA are approximately the same, while the precision of EGAAM is significantly better. Performance of the EGAAM was compared to that of classical GA on three families of multiminima test functions, which represent very difficult tests for any optimization algorithm because of the large number of local minima. After proving that EGAAM is a reliable tool for solving continuous optimization problems, its ability to determine model parameter values on the basis of known experimental results is demonstrated. Model parameter estimation is performed by minimizing the objective function, usually in the form of a sum of squared absolute or
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
relative differences between calculated and experimental values, which often have multiple minima. In many areas, the difficulty or inability of providing initial estimates of the parameter values arises frequently. Sometimes not even the order of magnitude of the model parameters can be correctly guessed. In such cases, application of the global optimization algorithm is imperative. Let us discuss in more detail the model parameter estimation problem and the model employed. One of the frequently employed models for modeling the optical properties of solids is the Lorentz-Drude ŽLD. oscillator model w7,12–15x. The LD model assigns oscillators to major critical points ŽCPs. in joint density of states to model direct interband transitions, with some additional oscillators to model absorption between CPs. Each oscillator is characterized by oscillator strength, frequency and damping constant. The initial values of the oscillator frequencies can be estimated from band structure calculation, while for the oscillator strengths and damping constants even the order of magnitude is unknown. Classical optimization techniques Žsimplex, Levenberg-Marqardt, etc.. require initial estimates close to the final values to obtain a meaningful solution, so that the fitting procedure is reduced to trial-and-error and involves re-running the fitting routine many times with different initial values before the acceptable fit is obtained, if it could be obtained at all. The main aim of this work was to devise a versatile continuous optimization method suitable for model parameter estimation, so that the optical properties of various materials could be investigated by simply changing the set of experimental data, andror changing the subroutine that calculates the model equations if a different type of material is involved Žfor instance, there is no Drude part that models intraband absorption for insulators.. Applications of global optimization methods such as the Metropolis algorithm and simulated annealing to model parameter estimation have been reported recently w14–19x. This is often a very complex optimization problem, since the objective function has a large number of local minima and the parameter values are required to have a certain physical meaning. Therefore, it is important not only to achieve a sufficiently low value of the objective function, but the location in parameter space also matters. Even for only three parameters, the global optimization technique, namely simulated annealing, was clearly superior to the gradient method w16x. However, the simulated annealing algorithm shows a certain dependence on the initial values in practice, which can be reduced but not completely eliminated by careful choice of cooling schedule and move generation procedure. The main advantage in applying GAs to model parameter estimation is the fact that initial estimates are not required. However, the problem of achieving satisfactory precision was encountered – to calculate accurately the optical properties with the LD model, at least two significant digits in estimated model parameter values are re-
149
quired, and for some of the parameters the only a priori available information is that they can have a value between 10y4 and 10. Independence of the initial values together with arbitrary precision in locating the global optimum can be achieved with combinations of simulated annealing and GAs, but the CPU time requirements of such algorithms are immense w20x. Therefore, improvement of the precision of locating the minimum is achieved in EGAAM by introducing adaptive parameter space size in the mutation process. In this paper it is shown that the EGAAM algorithm obtains parameter values of optical constants that give theoretical data exhibiting good agreement with experimental values for all investigated metals. In the following section the investigated algorithms, GA and EGAAM are described. Section 3 is devoted to the comparison of the performance of these algorithms applied on nine multiminima test functions. In Section 4 a short description of the applied Lorentz-Drude ŽLD. model for the optical constants of metals is given and EGAAM was used to estimate the parameters for beryllium, chromium, nickel and palladium.
2. Description of the algorithm In the implementation of GA we must define the representation of chromosomes, generation of initial population and genetic operators: selection, reproduction and mutation. In the following the chromosome representation and genetic operators employed in the investigated algorithms are presented. 2.1. Representation of chromosomes and population generation Continuous variables can be handled either directly, through real-valued Žfloating-point. representation and appropriate genetic operators, or by standard binary representation schemes and standard genetic operators. In case of binary representation, real values are approximated to the necessary degree with a fixed-point binary scheme, or the logarithm of the variable is encoded, thus reducing the required number of bits. However, since floating point representation w5,9x proved to be more convenient for continuous optimization problems, it was applied in this paper. In floating-point chromosome representation, each gene has the value of the corresponding variable pŽ k ., k s 1, n v , where n v is the number of variables. The values pŽ k . in chromosomes of the initial population are generated according to the formula p Ž k . s p l Ž k . q Ž pu Ž k . y p l Ž k . . r ,
Ž1.
where r is a random number r g w0,1x, and p lŽ k . and puŽ k . are initially set boundaries. In such a way, confinement of variables in the specified domain is achieved insuring that all variables have physical interpretation, so
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
150
that, for instance, in case of model parameter determination we cannot obtain a negative value for the frequency. 2.2. Selection and reproduction The many different existing selection methods can be divided in two categories: random selection methods and selection methods based on the fitness measure. The former methods, like roulette wheel w6x and similar methods w21x introduce stochastic errors. These errors are reduced in the latter ones by taking into account the fitness of an individual, like in the binary tournament method w22x, ‘‘stochastic remainder sampling without replacement’’ w23x, or any other selection method incorporating the following concept – the more fitted the chromosome, the higher the chance of being selected as a parent w3,4,7,24–26x. For reducing the bias of the selection method to highly fit individuals and thus preventing premature convergence, fitness scaling can be employed w23,26x. In this paper, an elitist selection mechanism w24,25,27x has been employed. In elitist selection, Ps percent of the new generation is produced by selection, and Pc percent is produced by crossover. Ns s NPs strings with the best fitness, where N is the number of strings in the population which enter directly the next generation. The Nc s NPc strings in the new population are generated by crossover among the parent strings which were chosen by fitness proportionally between all the strings in the current population. The probability that a string will have an offspring in the next generation is given by N
F Ž i . s f Ž i .r
Ý f Ži. ,
Ž2.
is1
where f Ž i . is the fitness value of the ith string. 2.3. CrossoÕer Operation of crossover exchanges subsets of elements between two parent chromosomes. If the subset consists of adjacent elements, it is an ‘‘ordered combination’’ crossover, while in ‘‘uniform combination’’ crossover each element is randomly chosen w21x. Ordered combination crossover can be one-point or two-point Žpoints of crossover are randomly selected, and elements between them are swapped.. Various crossover operators for continuous variables, like intermediate or line recombination w4x or crossover by different continuous mixing functions w3x, have been introduced. Deb and Agrawal w10x have suggested simulated binary crossover ŽSBX., which has search power similar to single-point crossover used in binary GAs. The main advantage of this crossover was demonstrated when the initial population was generated in bounds that do not bracket an optimum variable value, since the SBX operator enables searching outside the initial boundaries. However, this feature represents an important short-
coming for model parameter estimation, since there is no guarantee that the obtained solution would be physically meaningful Žfor instance, the LD model parameters must be positive.. Eshelman and Schaffer w11x suggested the blend crossover operator BLX-a which randomly picks a point in the range Ž p 1 y a Ž p 2 y p1 ., p1 q a Ž p 2 y p1 .. and the best results were reported for a s 0.5. It should be pointed out that the suggested mechanisms have been tested on functions up to 20 variables w4,10x, while here it was required that the algorithm can find global minimum of multimodal test function up to 100 variables. In the algorithms investigated here, uniform combination crossover was employed, providing that the parameter values are not altered or destroyed during the crossover, while the new values are introduced in the process of mutation. Crossover is performed by generating a random integer N1 g w n min ,n v x, where n v is the number of variables, i.e. the number of elements in the strings, and n min is the minimal number of elements exchanged in the crossover. Best results were obtained for n min s n vr2. When the number of elements to be exchanged is determined, random integers n i g w1,n par x, i s 1, N1 are generated and the elements at positions n i are swapped. 2.4. Mutation, concept of adaptiÕe mutations Mutation is necessary for maintaining a certain diversity in the population, thus preventing quick convergence to a local minimum. Mutation is usually performed by randomly altering the individual genes with probability Pm . Real coded GAs usually perturb the solution a little around the current value, which can be done with a uniformly or normally distributed step, or with a specifically designed manner of choosing the random step w4x, while binary GAs usually adopt a bitwise complement operator. The mutation step can be constant, or vary with the number of generations w28x or with the number of successful mutations w29x. If the mutation probability is too large, mutation no longer improves the performance of the population, because it enables losses of genetic information which could cause poor convergence w30x. Random mutations can also be performed either by generating an aleatory new value of individual gene with probability Pm w7x, or by introducing Pm percent of strings in the new generation whose elements are chosen randomly w27x. Although Pm is much smaller than Pc , and ergo mutation is usually considered a background operator w1,9x that ensures that the probability of searching any particular subspace is never equal to zero, replacing random mutations with adaptive mutations substantially improves the performance of the algorithm. The main problem in applying GAs to continuous optimization problems is discrete sampling of the solution space, which results in a roughly determined location of the global minimum. More precise location of the mini-
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
mum requires larger number of chromosomes in the population, thus demanding extensive computer resources, while the random character of population generation prevents any guarantee that the global minimum will be precisely located even if the algorithm is rerun for a number of times. Various schemes of crossover developed for continuous variables intend to introduce new values of genes, but such concepts give new values which are in certain manner bounded to existing ones. Therefore, in this paper the introduction of new values is accomplished by an adaptive mutation process. In the current generation, the average value m ˆ Ž k . of the parameter pŽ k . is computed, and Pm percent of the chromosomes in the next generation are formed by generating their genes in the same manner as during the creation of the initial population, but in narrowed boundaries. New boundaries for each parameter are determined according to
151
where d is a real number having a value between 0 and 1. The concept of such adaptive mutations is incorporated in EGAAM. To investigate clearly the influence of adaptive mutations in the performance of the algorithm, both EGAAM and classical GA have the same selection and crossover mechanism, they differ only in the mutation procedure. In the classical GA investigated here, uniform distributed mutation is performed by changing the parameter value with probability Pm . A new parameter value is given by pmut Ž k . s p Ž k . q sgn D p Ž k . ,
Ž5.
where pmut Ž k ., pŽ k . are the values of parameter k after and before mutations, respectively, and sgn is a random number in the interval wy1,1x, while D pŽ k . is the step value for the parameter k.
pnew - u Ž k . s pold - u Ž k . y c . Ž pold - u Ž k . y m ˆ Ž k.. ,
Ž3.
3. Test of the GA and EGAAM
pnew - l Ž k . s pold - l Ž k . q c . Ž m ˆ Ž k . y pold - l Ž k . . ,
Ž4.
To compare the performance of EGAAM and GA, experiments on three families of test functions, for 20, 50 and 100 variables were performed. It is well known that the performance of GA depends strongly on the parameters of the algorithm: N, Pc , Ps and Pm , and in the case of EGAAM, two additional parameters c and d. In the results presented here the influence of the number of chromosomes in the population N was investigated, since it is only one of the above parameters that significantly influences the required computational time. Before that, test runs on all test functions were performed to determine the optimal values for other algorithm parameters. The best results for classical GA were obtained for Pc s 0.2, Ps s
where m ˆ Ž k . is the average value of the parameter pŽ k . in the current population, and c is a predetermined positive number 0 - c - 1. In such a manner, a specified number Ž Nm s NPm . of new chromosomes is introduced in every generation. During the evolution, while m ˆ Ž k . changes towards the optimal value, the parameter values in the new chromosomes are more and more concentrated around m ˆ Ž k ., providing finer structure and more significant digits in the obtained variable values. For preventing excessive narrowing of the boundaries, their extreme values are set to pmax - lŽ k . s m ˆ Ž k .Ž1 y d . and pmin - uŽ k . s mˆ Ž k .Ž1 q d .,
Fig. 2. Function g for two variables.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
152
Fig. 3. Comparison of the algorithms for the function g of x i g wy10,10x, i s 1,20, for 2000, 1000, 500 and 250 strings in the population.
0.8 and Pm s 0.01, while the optimal values for EGAAM were Pc s 0.35, Ps s 0.6 and Pm s 0.05. Parameter d was set to value 0.05, while c s ccrn max , where n max is the maximal number of generations. The results given here correspond to n max s 200 and cc s 3. The first family of multiminima functions was employed by Aluffi-Pentini et al. w31x and Dekkers and Aarts w32x. It is given by
p gŽ x. s
n
ny 1
k 1sin2p y 1 q
Ý Ž yi y k 2 . 2 is1 2
= Ž 1 q k 1sin2p yiq1 . q Ž yn y k 2 . ,
Ž6.
where yi s 1 q 0.25Ž x i q 1., k 1 s 10, k 2 s 1. This func-
tion, shown for two variables in Fig. 2, has roughly 5 n local minima for x i g wy10,10x, i s 1,n. In Ref. w32x this function was tested for three variables, while in this paper tests with 20, 50 and 100 variables were performed. In order to determine the influence of the number of chromosomes in the population to the final objective function values, calculations for 2000, 1000, 500 and 250 chromosomes were carried out. The obtained results for 20, 50 and 100 variables are presented in Fig. 3, Fig. 4, and Fig. 5, respectively. It can be observed that EGAAM in all cases obtains values lower by several orders of magnitude, compared to the results of GA, and its performance depends considerably less on the number of chromosomes in the population. In other words, EGAAM gives satisfactory results even for 250 chromosomes in the population when GA fails to locate the minimum.
Fig. 4. Comparison of the algorithms for the function g of x i g wy10,10x, i s 1,50, for 2000, 1000, 500 and 250 strings in the population.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
153
Fig. 5. Comparison of the algorithms for the function g of x i g wy10,10x, i s 1,100, for 2000, 1000, 500 and 250 strings in the population.
A section of the second investigated family of functions f Ž x . along an axis for values of a s 0.2, b s 0.1 and c s 2 is shown in Fig. 6. It is given by n
f Ž x. s
Ý ax i2 q bx i2 sin cx i ,
Ž7.
is1
The obtained final objective function values for 20, 50 and
100 variables where x i g wy10,10x, i s 1,n, are given in Fig. 7. For the sake of comprehensibility of presentation, only results for N s 2000 chromosomes in the population are shown, since the dependence of the final solution on the number of chromosomes in the population is similar for all investigated functions. For this family of test functions, EGAAM is once again superior to GA, obtaining significantly lower Ž4–5 orders of magnitude. objective function values.
Fig. 6. Section of f along one axis for a s 0.2, b s 0.1 and c s 2.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
154
Fig. 7. Comparison of the algorithms for the function f of x i g wy10,10x, i s 1,n; n s 20, 50 and 100, for 2000 strings in the population.
The third family of multiminima functions, also investigated by Aluffi-Pentini et al. w31x and Dekkers and Aarts w32x, is given by ny 1
½
h Ž x . s k 3 sin2 Ž p k 4 x 1 . q
Ý Ž x i y k5 .2 is1
2
= 1 q k 6 sin Ž p k 4 x iq1 . 2
q Ž x n y k 5 . 1 q k 6 sin2 Ž p k 7 x n .
5
,
Ž8.
where k 3 s 0.1, k 4 s 3, k 5 s 1, k 6 s 1 and k 7 s 2. In Ref. w32x this function, whose two-dimensional segment is shown in Fig. 8, was investigated for five variables
x i g wy5,5x, i s 1, . . . ,5, i.e. in the area where there are roughly 155 minima. In this paper, a family of test functions hŽ x . was investigated for 20, 50 and 100 variables x i g wy10,10x, i s 1,n. The obtained results, also for N s 2000 strings in the population, are presented in Fig. 9. For this function EGAAM also achieves lower objective function values, by about three orders of magnitude.
4. Application to beryllium, chromium, nickel and palladium In this section, the applied Lorentz-Drude ŽLD. model for the optical dielectric function, which was often em-
Fig. 8. Function h for two variables.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
155
Fig. 9. Comparison of the algorithms for the function h of x i g wy10,10x, i s 1,n; n s 20, 50 and 100, for 2000 strings in the population.
ployed for modeling the optical constants of metals w15,14,12x is briefly discussed. It was shown w33–36x that the dielectric constant e r Ž v . can be expressed in the following form,
e rˆ Ž v . s e rˆŽ f . Ž v . q e rˆŽ b. Ž v . ,
Ž9.
which separates explicitly the intraband effects Žfrequently referred to as free electron effects. from interband effects Žfrequently referred to as bound electron effects.. The intraband part e rˆŽ f . Ž v . of the dielectric constant is described by the well-known free electron or Drude model w37x
e rˆŽ f . Ž v . s 1 y
V p2 v Ž v q i G0 .
,
Ž 10.
and the interband part of the dielectric constant e rˆŽ b. Ž v . is described by the simple semi-quantum model resembling the Lorentz result for insulators w12x f j v p2
k
e rˆŽ b. Ž v . s y
Ý js1
Žv
2
y v j2
. q i vGj
,
Ž 11.
where v p is the plasma frequency, k is the number of interband transitions with frequency v j , oscillator strength f j and lifetime 1rGj , while V p s f 0 v p is associated with intraband transitions with oscillator strength f 0 and damping constant G 0 . The following objective function was used:
'
isN
EŽ p. s
Ý is1
q
exp e r1 Ž v i . y e r1 Ž vi . exp e r1 Ž vi .
e r 2 Ž v i . y e rexp 2 Ž vi . e rexp 2 Ž vi .
2
.
Ž 12.
The values of the plasma frequencies of these metals were determined according to the definition
vp s
N e2
1r2
ž / me0
,
Ž 13 .
where N is the concentration of the valence electrons. For all metals four oscillators were used, and hence there are 14 parameters to be determined. ŽSee Table 1.. 4.1. Beryllium Beryllium has many unique properties: its density is 3.5 times higher than that of Li, and also greater than the density of Mg and Ca. It has very low compressibility and electronic specific heat, its melting temperature is much higher than for the other alkaline earth elements and its Table 1 Values of the LD model parameters Metal
Be
Cr
Ni
Pd
vp f0 G0 f1 G1 v1 f2 G2 v2 f3 G3 v3 f4 G4 v4
18.51 0.081 0.034 0.050 1.345 0.227 0.101 3.441 1.325 0.388 3.862 2.882 0.260 2.916 4.532
10.75 0.168 0.046 0.199 2.622 0.162 0.257 2.280 0.820 0.952 2.436 1.981 0.848 0.155 8.552
15.92 0.089 0.027 0.058 2.793 0.099 0.150 1.362 0.491 0.100 1.962 1.585 0.714 6.018 5.847
9.72 0.331 0.008 0.763 3.487 0.392 0.098 0.512 0.496 0.523 3.815 1.284 0.479 3.387 5.753
156
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
Debye temperature is among the highest in the periodic table. Also, Be is relatively light compared to the other optically important metals, has low Z number, low thermal expansion at cryogenic temperatures, very high specific stiffness Žratio of Young’s modulus to density., high thermal conductivity and high reflectivity in the infrared. Because of these properties Be is particularly attractive for optical components intended for use in space. Unfortunately, only a few studies of the optical properties of Be have been done, and they have produced widely different results, especially in the visible range of the spectrum. The variability of the results may be due to difficulties in sample preparation w38x. Bulk samples are usually produced by high-pressure sintering of microcrystalline powder, which results in the incorporation of significant amounts of BeO in the sample. For many metals oxidation problems can be overcome by vacuum evaporation, but for Be the optical properties depend critically on the evaporation procedure. In this paper tabulated data given in the study of Arakawa et al. w38x were used. The data consists of the results of Hull w39x in the range 0.06–0.8 eV, while in the ranges 0.8–5.0 eV and 1.5–14 eV results of Hunderi et al. w40x and Seignac et al. w41x were used, respectively. In regions of overlap, higher reflectance values were favored. Fig. 10 shows the real and imaginary parts of the optical dielectric constant as a function of energy. Open circles represent tabulated data, and the solid line the LD
model. The inset shows the structure between 0.6 eV and 2.0 eV. It can be observed that the model describes the tabulated data within the experimental error. 4.2. Chromium Chromium is a body-centered cubic metal, antiferromagnetic below 312 K. The magnetic ordering, with spatial period incommensurate with the crystal lattice period, introduces structure in the spectra of the optical constants in the infrared w42x. Cr oxidizes easily, but the Cr2 O 3 overlayer is impervious to Cr or O after a thickness of ˚ is reached. Hence, there is no further oxidation about 30 A at room temperature. Such an overlayer may affect the optical constants in the visible and near infrared obtained by ellipsometry and by any other technique in the vacuum ultraviolet. There have been few studies for which the sample was both prepared and measured in UHV. Data of Bos and Lynch w43x, tabulated in Ref. w42x, resulting from Kramers-Kronig analysis of low-temperature absorptance data on bulk samples, were employed for LD model parameter determination. The measurements of Bos and Lynch extend to 13 mm in the infrared, although their Kramers-Kronig analysis extends further. A number of measurements extend to longer wavelengths, but they do not agree well with each other, and some of them show weak structures which are probably spurious. In this article the data obtained by Kirillova w44x, measured in the range
Fig. 10. Real and imaginary parts of the dielectric constant of Be versus energy Žsolid line – model, open circles – experimental data.; the inset shows the structure between 0.6 eV and 2.0 eV.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
157
Fig. 11. Real and imaginary parts of the dielectric constant of Cr versus energy Žsolid line – model, open circles – experimental data.; the inset shows the structure around 1 eV.
0.069–4.901 eV, and results of Johnson and Christy w45x, measured in the range 0.64–6.10 eV, were also employed. Fig. 11 shows the real and imaginary parts of the dielectric constant of chromium Žopen circles – experimental data, solid line – LD model.. The inset shows the structure around 1 eV. Good agreement with experimental data can be observed. 4.3. Ni The overall agreement among the published n and k spectra of Ni is good. In this article tabulated data from the
study of Lynch et al. w46x, consisting of data obtained by Vehse and Arakawa w47x in the range 3.0–25 eV, and data obtained by Lynch et al. w48x in the range 0.1–3.0 eV were employed. In Fig. 12, showing the real and imaginary parts of the optical dielectric function versus energy Žopen circles – tabulated data, solid line – LD model. good agreement between tabulated and calculated data can be observed. 4.4. Pd The optical properties of Pd have been studied under a variety of experimental conditions and by several tech-
Fig. 12. Real and imaginary parts of the dielectric constant of Ni versus energy Žsolid line – model, open circles – experimental data..
158
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
Fig. 13. Real and imaginary parts of the dielectric constant of Pd versus energy Žsolid line – model, open circles – experimental data..
niques. Unfortunately, many results disagree with the actual magnitude of the optical structures, although the structure energies agree almost in every case. This depends on sample preparation, accuracy of the optical measurements and the method of analysis of the experimental data. In this article the data tabulated in the study of Borghesi et al. w49x were used. These data were obtained from the following works and references therein: Vehse et al. w50x in the range 6.0–12.5 eV, Johnson and Christy w45x in the range 0.6–6.0 eV and Weaver et al. w51,52x in the range 0.1–0.6 eV. In Fig. 13 excellent agreement between calculated data and tabulated values of the real and imaginary parts of the optical dielectric constant can be observed. Open circles represent experimental data, the solid line the data calculated by the LD model.
5. Conclusion In this paper the problem of applying genetic algorithms to continuous optimization was discussed. To this end, the prime concern was the possibility of applying these algorithms to a very specific and delicate problem: determination of the model parameter values by minimizing the difference between experimental and calculated data. In this case, precise and reliable location of the global minimum is imperative, since slightly different objective function values can be obtained for significantly different parameter values. Elite genetic algorithms with adaptive mutations were proposed for overcoming the difficulties in continuous optimization originating from discrete sampling of the solution space. Performance of the new algorithm was compared to that of the classical GA for three families of multiminima test functions, for up to
100 variables and different number of chromosomes in the population. In all investigated cases, EGAAM obtained extensively lower objective function values than GA. After verifying its ability to locate the global minimum on nine multiminima test functions, EGAAM was applied to model parameter determination of the optical constants of Be, Cr, Ni and Pd. For all four metals excellent agreement between calculated and experimental results is obtained. References w1x D.E. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley, Reading, 1989. w2x S. Kirkpatrick, C.D. Gelatt Jr., M.P. Vecchi, Science 220 Ž1983. 671. w3x M.W. Gutowski, J. Phys. A 27 Ž1994. 7893. w4x H. Muhlenbein, D. Schlierkamp-Voosen, Evolutionary Com¨ putation 1 Ž1993. 25. w5x K.P. Wong, Y.W. Wong, in: Proc. ANZIIS-93, Perth, Western Australia, 1993, pp. 512–516. w6x K.P. Wong, Y.W. Wong, IEE Proc. Gen. Transm. Distrib. 141 Ž1994. 507. w7x A.B. Djurisic, ˇ ´ J.M. Elazar, A.D. Rakic, ´ Optics Comm. 134 Ž1997. 407. w8x A. Chipperfield, R. Fleming, Control Computers 23 Ž1995. 88. w9x D.E. Goldberg, Complex Systems 5 Ž1991. 139. w10x K. Deb, R.B. Agrawal, Complex Systems 9 Ž1995. 115. w11x J.L. Eshelman, J.D. Schaffer, in: Proc. Foundations of GA Workshop, 1992, pp. 187–202. w12x A.D. Rakic, ´ Appl. Optics 34 Ž1995. 4755. w13x C.J. Powell, J. Opt. Soc. Am. 60 Ž1970. 78. w14x A.B. Djurisic, ˇ ´ A.D. Rakic, ´ J.M. Elazar, Phys. Rev. E 55 Ž1997. 4797. w15x A.D. Rakic, ´ J.M. Elazar, A.B. Djurisic, ˇ ´ Phys. Rev. E 52 Ž1995. 6862.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159 w16x T.H. Han, W.W. Chang, Electron. Lett. 32 Ž1996. 1256. w17x A. Franke, A. Stendal, O. Stenzel, C. von Borczyskowski, Pure Appl. Optics 5 Ž1996. 845. w18x A.D. Rakic, ´ M.L. Majewski, J. Appl. Phys. 80 Ž1996. 5909. w19x M.K. Vai, S. Prasad, N.C. Li, F. Kai, IEEE Trans. Electron. Devices 36 Ž1989. 761. w20x A.B. Djurisic, ˇ ´ A.D. Rakic, ´ J.M. Elazar, Appl. Optics 36 Ž1997., to be published. w21x F. Curatelli, Int. J. Electronics 78 Ž1995. 435. w22x S.E. Cienawski, J.W. Ehart, S. Ranjithan, Water Resources Research 31 Ž1995. 399. w23x Q. Li, E.J. Rotwell, K.-M. Chen, D.P. Nyquist, IEEE Trans. Antennas Propag. 44 Ž1996. 198. w24x R. Vemuri, R. Vemuri, Electron. Lett. 30 Ž1994. 1270. w25x S.H. Clearwater, T. Hogg, Artificial Intelligence 81 Ž1996. 327. w26x C.B. Lucasius, M.L.M. Beckers, G. Kateman, Anal. Chim. Acta 286 Ž1994. 135. w27x R.R. Brooks, S.S. Iyengar, J. Chen, Artificial Intelligence 81 Ž1996. 327. w28x U.D. Hanebeck, G.K. Schmidt, Fuzzy Sets Systems 79 Ž1996. 59. w29x T. Back, ¨ H.P. Schwefell, in: Genetic Algorithms in Engineering and Computer Science, Wiley, Chichester, 1995, pp. 111–140. w30x D. Raynolds, J. Gonatann, Artificial Intelligence 82 Ž1996. 303. w31x F. Aluffi-Pentini, V. Parisi, F. Zirilli, J. Optim. Theory Appl. 47 Ž1985. 1. w32x A. Dekkers, E. Aarts, Mathematical Program. 50 Ž1991. 367. w33x N.W. Ashcroft, K. Sturm, Phys. Rev. B 3 Ž1971. 1898. w34x H. Ehrenreich, H.R. Philipp, B. Segall, Phys. Rev. 132 Ž1963. 1918.
w35x w36x w37x w38x
w39x w40x w41x w42x
w43x w44x w45x w46x
w47x w48x w49x
w50x w51x w52x
159
H. Ehrenreich, H.R. Philipp, Phys. Rev. 128 Ž1962. 1622. K. Sturm, N.W. Ashcroft, Phys. Rev. B 10 Ž1974. 1343. P. Drude, The Theory of Optics, Dover, New York, 1959. E.T. Arakawa, T.A. Callcott, Y.-C. Chang, in: E.D. Palik ŽEd.., Handbook of Optical Constants of Solids II, Academic Press, San Diego, CA, 1991. A.B. Hull, Perkin-ElmerrAppl. Optics Operations, Box 3115, Garden Grave, CA 92641, private communication. O. Hunderi, P. Myers, J. Phys. F 4 Ž1974. 1088. A. Seignac, S. Robin, Solid State Commun. 19 Ž1976. 343. D.W. Lynch, R.W. Hunter, in: E.D. Palik ŽEd.., Handbook of Optical Constants of Solids II, Academic Press, San Diego, CA, 1991. L.W. Bos, D.W. Lynch, Phys. Rev. B 2 Ž1970. 4567. M.M. Kirillova, M.M. Noskov, Phys. Met. Metallogr. 26 Ž1968. 189. P.B. Johnson, R.W. Christy, Phys. Rev. B 11 Ž1975. 1315. D.W. Lynch, W.R. Hunter, in: E.D. Palik ŽEd.., Handbook of Optical Constants of Solids, Academic Press, Orlando, FL, 1985, pp. 275–367. R.C. Vehse, E.T. Arakawa, Phys. Rev. 180 Ž1969. 695. D.W. Lynch, R. Rosei, J.H. Weaver, Solid State Commun. 9 Ž1971. 2195. A. Borghesi, A. Piaggi, in: E.D. Palik ŽEd.., Handbook of Optical Constants of Solids II, Academic Press, San Diego, CA, 1991. R.C. Vehse, E.T. Arakawa, M.W. Williams, Phys. Rev. B 1 Ž1970. 517. J.H. Weaver, R.L. Benbow, Phys. Rev. B 12 Ž1975. 3509. J.H. Weaver, C. Krafka, D.W. Lynch, E.E. Koch, in: Physik, Daten, Physics Data, Optical Properties of Metals, vol. 18-1, Fach-information zentrum, Karlsruhe, 1981, p. 179.
Optics Communications 151 Ž1998. 147–159
Full length article
Elite genetic algorithms with adaptive mutations for solving continuous optimization problems – application to modeling of the optical constants of solids Aleksandra B. Djurisic ˇ´ Department of EEE, UniÕersity of Hong Kong, Pokfulam Road, Hong Kong, China Received 17 September 1997; revised 6 January 1998; accepted 19 February 1998
Abstract The elite genetic algorithm with adaptive mutations is proposed as a tool for solving continuous optimization problems. The new algorithm and the corresponding classical genetic algorithm were severely tested on three families of multiminima test functions for 20, 50 and 100 variables. All performed tests proved that the introduced adaptive mutations significantly improve the ability of the algorithm to find a global minimum. After verifying the superiority of the proposed algorithm over the classical one on the test functions, an elite genetic algorithm with adaptive mutations was applied for solving the model parameter determination problem for modeling optical constants of the following metals: beryllium, chromium, nickel and palladium. Good agreement between calculated and experimental data was obtained for all four metals. q 1998 Elsevier Science B.V. All rights reserved. PACS: 02.70.–c; 78.20.Ci; 78.66.–w Keywords: Genetic algorithms; Model parameter estimation; Optical constants
1. Introduction Genetic algorithms ŽGAs. w1x are stochastic global search methods that mimic the concept of natural evolution. Due to the nature of the algorithm, their successful application was mostly restricted to optimization problems whose solution can be conveniently represented in binary form. However, there is a rising interest in applying genetic algorithms to continuous optimization problems, especially since there is no need for initial estimates, which is an important advantage of GAs over other stochastic search methods such as, for instance, simulated annealing w2x. For that reason, various modifications of original GAs have been reported w3–7x. GAs search for optimal solution employing mechanisms of natural evolution: selection, mating and mutation, which are applied to the set of possible problem solutions, called population. Each element of the population is represented by a vector of variables, called string or chromosome. Each element of the chromosome is, by analogy,
termed a gene. Chromosomes are characterized by their ‘‘fitness’’ – performance with respect to some objective function. Highly fit individuals have a high probability of being selected for reproduction Žthey may survive or give offsprings in the next generation.. Strings with lower fitness have correspondingly lower probability for transferring their genes to the following generations. This is illustrated in Fig. 1, where shade of gray represents fitness – the lighter shade corresponds to the better gene. Many variations of the simple GA have been developed w8x, differing in chromosome representation, selection, reproduction and mutation, existence of fitness scaling, etc. It was shown that floating-point number representation, instead of binary string representation common for classical GAs, is more convenient for continuous optimization applications w5,6,9x. In such a manner, by representing variable values with real numbers, the length of the chromosome is equal to the number of variables, thus being much smaller than in the case of binary coding. For instance, for coding real numbers p - 2 5 with four deci-
0030-4018r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII S 0 0 3 0 - 4 0 1 8 Ž 9 8 . 0 0 1 2 2 - 9
148
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
Fig. 1. Illustration of GA: shade of gray represents fitness: the lighter shade correspods to a more fitted gene.
mal digits we need 21 genes, i.e. the length of chromosome is about twenty times larger. Also, an important feature of floating-point representation is that the variable values cannot be altered or destroyed during crossover operation, while in the case of binary coding such undesired changes of variable values can lead to loss or deterioration of genetic information. More important, conversion of decimal numbers into binary ones and vice versa is avoided by using floating-point representation. The important advantage of real-coded GAs is the absence of Hamming cliff problem inherent to binary coded GAs w1,9,10x. Such problems could be overcome by the use of Gray coding, but that leads to higher nonlinearities with respect to recombination. Also, using real-coding reduces the dimensionality of the problem, thus reducing the opportunity for deception, but the possibility of new obstacles to convergence is introduced. The fact that realcoded GAs can in certain cases be blocked from further progress has been recognized and discussed w9x. Recently, there has been much work in modification of real-coded GAs in order to make them as successful in solving continuous optimization problems as their binary counterparts are in solving discrete optimization problems w3,4,10,11x. However, in solving certain problems such as, for example, parameter estimation problem, even GAs with floating-point representation and appropriate selection, crossover and mutation mechanisms often fail to give satisfactory results, especially for a large number of variables. These shortcomings are manifested mainly due to the discrete sampling of the solution space, while variables have continuous nature. GAs are capable of roughly locat-
ing a global optimum, but a huge number of chromosomes in the population is required for any refinements. Recently, genetic algorithms with parameter space size adjustment have been proposed for model parameter determination w7x. The main idea is narrowing of the solution space by improving the guess of the global minimum position. This is achieved by introducing one additional loop. The initial population is generated, classical GA is performed, and the boundaries for each parameter are narrowed towards the obtained average value for that parameter in the inner loop final population. The new initial population is generated in the new boundaries and GA is performed again. The outer loop is executed for a specified number of iterations n max , which means that the time required for performing this algorithm is n max times larger than that for the corresponding classical GA. In this paper, an elite genetic algorithm with adaptive mutations ŽEGAAM. is proposed, which can successfully solve the presented problems, without any observable increase of CPU time requirements. In other words, the execution time of EGAAM and GA are approximately the same, while the precision of EGAAM is significantly better. Performance of the EGAAM was compared to that of classical GA on three families of multiminima test functions, which represent very difficult tests for any optimization algorithm because of the large number of local minima. After proving that EGAAM is a reliable tool for solving continuous optimization problems, its ability to determine model parameter values on the basis of known experimental results is demonstrated. Model parameter estimation is performed by minimizing the objective function, usually in the form of a sum of squared absolute or
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
relative differences between calculated and experimental values, which often have multiple minima. In many areas, the difficulty or inability of providing initial estimates of the parameter values arises frequently. Sometimes not even the order of magnitude of the model parameters can be correctly guessed. In such cases, application of the global optimization algorithm is imperative. Let us discuss in more detail the model parameter estimation problem and the model employed. One of the frequently employed models for modeling the optical properties of solids is the Lorentz-Drude ŽLD. oscillator model w7,12–15x. The LD model assigns oscillators to major critical points ŽCPs. in joint density of states to model direct interband transitions, with some additional oscillators to model absorption between CPs. Each oscillator is characterized by oscillator strength, frequency and damping constant. The initial values of the oscillator frequencies can be estimated from band structure calculation, while for the oscillator strengths and damping constants even the order of magnitude is unknown. Classical optimization techniques Žsimplex, Levenberg-Marqardt, etc.. require initial estimates close to the final values to obtain a meaningful solution, so that the fitting procedure is reduced to trial-and-error and involves re-running the fitting routine many times with different initial values before the acceptable fit is obtained, if it could be obtained at all. The main aim of this work was to devise a versatile continuous optimization method suitable for model parameter estimation, so that the optical properties of various materials could be investigated by simply changing the set of experimental data, andror changing the subroutine that calculates the model equations if a different type of material is involved Žfor instance, there is no Drude part that models intraband absorption for insulators.. Applications of global optimization methods such as the Metropolis algorithm and simulated annealing to model parameter estimation have been reported recently w14–19x. This is often a very complex optimization problem, since the objective function has a large number of local minima and the parameter values are required to have a certain physical meaning. Therefore, it is important not only to achieve a sufficiently low value of the objective function, but the location in parameter space also matters. Even for only three parameters, the global optimization technique, namely simulated annealing, was clearly superior to the gradient method w16x. However, the simulated annealing algorithm shows a certain dependence on the initial values in practice, which can be reduced but not completely eliminated by careful choice of cooling schedule and move generation procedure. The main advantage in applying GAs to model parameter estimation is the fact that initial estimates are not required. However, the problem of achieving satisfactory precision was encountered – to calculate accurately the optical properties with the LD model, at least two significant digits in estimated model parameter values are re-
149
quired, and for some of the parameters the only a priori available information is that they can have a value between 10y4 and 10. Independence of the initial values together with arbitrary precision in locating the global optimum can be achieved with combinations of simulated annealing and GAs, but the CPU time requirements of such algorithms are immense w20x. Therefore, improvement of the precision of locating the minimum is achieved in EGAAM by introducing adaptive parameter space size in the mutation process. In this paper it is shown that the EGAAM algorithm obtains parameter values of optical constants that give theoretical data exhibiting good agreement with experimental values for all investigated metals. In the following section the investigated algorithms, GA and EGAAM are described. Section 3 is devoted to the comparison of the performance of these algorithms applied on nine multiminima test functions. In Section 4 a short description of the applied Lorentz-Drude ŽLD. model for the optical constants of metals is given and EGAAM was used to estimate the parameters for beryllium, chromium, nickel and palladium.
2. Description of the algorithm In the implementation of GA we must define the representation of chromosomes, generation of initial population and genetic operators: selection, reproduction and mutation. In the following the chromosome representation and genetic operators employed in the investigated algorithms are presented. 2.1. Representation of chromosomes and population generation Continuous variables can be handled either directly, through real-valued Žfloating-point. representation and appropriate genetic operators, or by standard binary representation schemes and standard genetic operators. In case of binary representation, real values are approximated to the necessary degree with a fixed-point binary scheme, or the logarithm of the variable is encoded, thus reducing the required number of bits. However, since floating point representation w5,9x proved to be more convenient for continuous optimization problems, it was applied in this paper. In floating-point chromosome representation, each gene has the value of the corresponding variable pŽ k ., k s 1, n v , where n v is the number of variables. The values pŽ k . in chromosomes of the initial population are generated according to the formula p Ž k . s p l Ž k . q Ž pu Ž k . y p l Ž k . . r ,
Ž1.
where r is a random number r g w0,1x, and p lŽ k . and puŽ k . are initially set boundaries. In such a way, confinement of variables in the specified domain is achieved insuring that all variables have physical interpretation, so
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
150
that, for instance, in case of model parameter determination we cannot obtain a negative value for the frequency. 2.2. Selection and reproduction The many different existing selection methods can be divided in two categories: random selection methods and selection methods based on the fitness measure. The former methods, like roulette wheel w6x and similar methods w21x introduce stochastic errors. These errors are reduced in the latter ones by taking into account the fitness of an individual, like in the binary tournament method w22x, ‘‘stochastic remainder sampling without replacement’’ w23x, or any other selection method incorporating the following concept – the more fitted the chromosome, the higher the chance of being selected as a parent w3,4,7,24–26x. For reducing the bias of the selection method to highly fit individuals and thus preventing premature convergence, fitness scaling can be employed w23,26x. In this paper, an elitist selection mechanism w24,25,27x has been employed. In elitist selection, Ps percent of the new generation is produced by selection, and Pc percent is produced by crossover. Ns s NPs strings with the best fitness, where N is the number of strings in the population which enter directly the next generation. The Nc s NPc strings in the new population are generated by crossover among the parent strings which were chosen by fitness proportionally between all the strings in the current population. The probability that a string will have an offspring in the next generation is given by N
F Ž i . s f Ž i .r
Ý f Ži. ,
Ž2.
is1
where f Ž i . is the fitness value of the ith string. 2.3. CrossoÕer Operation of crossover exchanges subsets of elements between two parent chromosomes. If the subset consists of adjacent elements, it is an ‘‘ordered combination’’ crossover, while in ‘‘uniform combination’’ crossover each element is randomly chosen w21x. Ordered combination crossover can be one-point or two-point Žpoints of crossover are randomly selected, and elements between them are swapped.. Various crossover operators for continuous variables, like intermediate or line recombination w4x or crossover by different continuous mixing functions w3x, have been introduced. Deb and Agrawal w10x have suggested simulated binary crossover ŽSBX., which has search power similar to single-point crossover used in binary GAs. The main advantage of this crossover was demonstrated when the initial population was generated in bounds that do not bracket an optimum variable value, since the SBX operator enables searching outside the initial boundaries. However, this feature represents an important short-
coming for model parameter estimation, since there is no guarantee that the obtained solution would be physically meaningful Žfor instance, the LD model parameters must be positive.. Eshelman and Schaffer w11x suggested the blend crossover operator BLX-a which randomly picks a point in the range Ž p 1 y a Ž p 2 y p1 ., p1 q a Ž p 2 y p1 .. and the best results were reported for a s 0.5. It should be pointed out that the suggested mechanisms have been tested on functions up to 20 variables w4,10x, while here it was required that the algorithm can find global minimum of multimodal test function up to 100 variables. In the algorithms investigated here, uniform combination crossover was employed, providing that the parameter values are not altered or destroyed during the crossover, while the new values are introduced in the process of mutation. Crossover is performed by generating a random integer N1 g w n min ,n v x, where n v is the number of variables, i.e. the number of elements in the strings, and n min is the minimal number of elements exchanged in the crossover. Best results were obtained for n min s n vr2. When the number of elements to be exchanged is determined, random integers n i g w1,n par x, i s 1, N1 are generated and the elements at positions n i are swapped. 2.4. Mutation, concept of adaptiÕe mutations Mutation is necessary for maintaining a certain diversity in the population, thus preventing quick convergence to a local minimum. Mutation is usually performed by randomly altering the individual genes with probability Pm . Real coded GAs usually perturb the solution a little around the current value, which can be done with a uniformly or normally distributed step, or with a specifically designed manner of choosing the random step w4x, while binary GAs usually adopt a bitwise complement operator. The mutation step can be constant, or vary with the number of generations w28x or with the number of successful mutations w29x. If the mutation probability is too large, mutation no longer improves the performance of the population, because it enables losses of genetic information which could cause poor convergence w30x. Random mutations can also be performed either by generating an aleatory new value of individual gene with probability Pm w7x, or by introducing Pm percent of strings in the new generation whose elements are chosen randomly w27x. Although Pm is much smaller than Pc , and ergo mutation is usually considered a background operator w1,9x that ensures that the probability of searching any particular subspace is never equal to zero, replacing random mutations with adaptive mutations substantially improves the performance of the algorithm. The main problem in applying GAs to continuous optimization problems is discrete sampling of the solution space, which results in a roughly determined location of the global minimum. More precise location of the mini-
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
mum requires larger number of chromosomes in the population, thus demanding extensive computer resources, while the random character of population generation prevents any guarantee that the global minimum will be precisely located even if the algorithm is rerun for a number of times. Various schemes of crossover developed for continuous variables intend to introduce new values of genes, but such concepts give new values which are in certain manner bounded to existing ones. Therefore, in this paper the introduction of new values is accomplished by an adaptive mutation process. In the current generation, the average value m ˆ Ž k . of the parameter pŽ k . is computed, and Pm percent of the chromosomes in the next generation are formed by generating their genes in the same manner as during the creation of the initial population, but in narrowed boundaries. New boundaries for each parameter are determined according to
151
where d is a real number having a value between 0 and 1. The concept of such adaptive mutations is incorporated in EGAAM. To investigate clearly the influence of adaptive mutations in the performance of the algorithm, both EGAAM and classical GA have the same selection and crossover mechanism, they differ only in the mutation procedure. In the classical GA investigated here, uniform distributed mutation is performed by changing the parameter value with probability Pm . A new parameter value is given by pmut Ž k . s p Ž k . q sgn D p Ž k . ,
Ž5.
where pmut Ž k ., pŽ k . are the values of parameter k after and before mutations, respectively, and sgn is a random number in the interval wy1,1x, while D pŽ k . is the step value for the parameter k.
pnew - u Ž k . s pold - u Ž k . y c . Ž pold - u Ž k . y m ˆ Ž k.. ,
Ž3.
3. Test of the GA and EGAAM
pnew - l Ž k . s pold - l Ž k . q c . Ž m ˆ Ž k . y pold - l Ž k . . ,
Ž4.
To compare the performance of EGAAM and GA, experiments on three families of test functions, for 20, 50 and 100 variables were performed. It is well known that the performance of GA depends strongly on the parameters of the algorithm: N, Pc , Ps and Pm , and in the case of EGAAM, two additional parameters c and d. In the results presented here the influence of the number of chromosomes in the population N was investigated, since it is only one of the above parameters that significantly influences the required computational time. Before that, test runs on all test functions were performed to determine the optimal values for other algorithm parameters. The best results for classical GA were obtained for Pc s 0.2, Ps s
where m ˆ Ž k . is the average value of the parameter pŽ k . in the current population, and c is a predetermined positive number 0 - c - 1. In such a manner, a specified number Ž Nm s NPm . of new chromosomes is introduced in every generation. During the evolution, while m ˆ Ž k . changes towards the optimal value, the parameter values in the new chromosomes are more and more concentrated around m ˆ Ž k ., providing finer structure and more significant digits in the obtained variable values. For preventing excessive narrowing of the boundaries, their extreme values are set to pmax - lŽ k . s m ˆ Ž k .Ž1 y d . and pmin - uŽ k . s mˆ Ž k .Ž1 q d .,
Fig. 2. Function g for two variables.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
152
Fig. 3. Comparison of the algorithms for the function g of x i g wy10,10x, i s 1,20, for 2000, 1000, 500 and 250 strings in the population.
0.8 and Pm s 0.01, while the optimal values for EGAAM were Pc s 0.35, Ps s 0.6 and Pm s 0.05. Parameter d was set to value 0.05, while c s ccrn max , where n max is the maximal number of generations. The results given here correspond to n max s 200 and cc s 3. The first family of multiminima functions was employed by Aluffi-Pentini et al. w31x and Dekkers and Aarts w32x. It is given by
p gŽ x. s
n
ny 1
k 1sin2p y 1 q
Ý Ž yi y k 2 . 2 is1 2
= Ž 1 q k 1sin2p yiq1 . q Ž yn y k 2 . ,
Ž6.
where yi s 1 q 0.25Ž x i q 1., k 1 s 10, k 2 s 1. This func-
tion, shown for two variables in Fig. 2, has roughly 5 n local minima for x i g wy10,10x, i s 1,n. In Ref. w32x this function was tested for three variables, while in this paper tests with 20, 50 and 100 variables were performed. In order to determine the influence of the number of chromosomes in the population to the final objective function values, calculations for 2000, 1000, 500 and 250 chromosomes were carried out. The obtained results for 20, 50 and 100 variables are presented in Fig. 3, Fig. 4, and Fig. 5, respectively. It can be observed that EGAAM in all cases obtains values lower by several orders of magnitude, compared to the results of GA, and its performance depends considerably less on the number of chromosomes in the population. In other words, EGAAM gives satisfactory results even for 250 chromosomes in the population when GA fails to locate the minimum.
Fig. 4. Comparison of the algorithms for the function g of x i g wy10,10x, i s 1,50, for 2000, 1000, 500 and 250 strings in the population.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
153
Fig. 5. Comparison of the algorithms for the function g of x i g wy10,10x, i s 1,100, for 2000, 1000, 500 and 250 strings in the population.
A section of the second investigated family of functions f Ž x . along an axis for values of a s 0.2, b s 0.1 and c s 2 is shown in Fig. 6. It is given by n
f Ž x. s
Ý ax i2 q bx i2 sin cx i ,
Ž7.
is1
The obtained final objective function values for 20, 50 and
100 variables where x i g wy10,10x, i s 1,n, are given in Fig. 7. For the sake of comprehensibility of presentation, only results for N s 2000 chromosomes in the population are shown, since the dependence of the final solution on the number of chromosomes in the population is similar for all investigated functions. For this family of test functions, EGAAM is once again superior to GA, obtaining significantly lower Ž4–5 orders of magnitude. objective function values.
Fig. 6. Section of f along one axis for a s 0.2, b s 0.1 and c s 2.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
154
Fig. 7. Comparison of the algorithms for the function f of x i g wy10,10x, i s 1,n; n s 20, 50 and 100, for 2000 strings in the population.
The third family of multiminima functions, also investigated by Aluffi-Pentini et al. w31x and Dekkers and Aarts w32x, is given by ny 1
½
h Ž x . s k 3 sin2 Ž p k 4 x 1 . q
Ý Ž x i y k5 .2 is1
2
= 1 q k 6 sin Ž p k 4 x iq1 . 2
q Ž x n y k 5 . 1 q k 6 sin2 Ž p k 7 x n .
5
,
Ž8.
where k 3 s 0.1, k 4 s 3, k 5 s 1, k 6 s 1 and k 7 s 2. In Ref. w32x this function, whose two-dimensional segment is shown in Fig. 8, was investigated for five variables
x i g wy5,5x, i s 1, . . . ,5, i.e. in the area where there are roughly 155 minima. In this paper, a family of test functions hŽ x . was investigated for 20, 50 and 100 variables x i g wy10,10x, i s 1,n. The obtained results, also for N s 2000 strings in the population, are presented in Fig. 9. For this function EGAAM also achieves lower objective function values, by about three orders of magnitude.
4. Application to beryllium, chromium, nickel and palladium In this section, the applied Lorentz-Drude ŽLD. model for the optical dielectric function, which was often em-
Fig. 8. Function h for two variables.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
155
Fig. 9. Comparison of the algorithms for the function h of x i g wy10,10x, i s 1,n; n s 20, 50 and 100, for 2000 strings in the population.
ployed for modeling the optical constants of metals w15,14,12x is briefly discussed. It was shown w33–36x that the dielectric constant e r Ž v . can be expressed in the following form,
e rˆ Ž v . s e rˆŽ f . Ž v . q e rˆŽ b. Ž v . ,
Ž9.
which separates explicitly the intraband effects Žfrequently referred to as free electron effects. from interband effects Žfrequently referred to as bound electron effects.. The intraband part e rˆŽ f . Ž v . of the dielectric constant is described by the well-known free electron or Drude model w37x
e rˆŽ f . Ž v . s 1 y
V p2 v Ž v q i G0 .
,
Ž 10.
and the interband part of the dielectric constant e rˆŽ b. Ž v . is described by the simple semi-quantum model resembling the Lorentz result for insulators w12x f j v p2
k
e rˆŽ b. Ž v . s y
Ý js1
Žv
2
y v j2
. q i vGj
,
Ž 11.
where v p is the plasma frequency, k is the number of interband transitions with frequency v j , oscillator strength f j and lifetime 1rGj , while V p s f 0 v p is associated with intraband transitions with oscillator strength f 0 and damping constant G 0 . The following objective function was used:
'
isN
EŽ p. s
Ý is1
q
exp e r1 Ž v i . y e r1 Ž vi . exp e r1 Ž vi .
e r 2 Ž v i . y e rexp 2 Ž vi . e rexp 2 Ž vi .
2
.
Ž 12.
The values of the plasma frequencies of these metals were determined according to the definition
vp s
N e2
1r2
ž / me0
,
Ž 13 .
where N is the concentration of the valence electrons. For all metals four oscillators were used, and hence there are 14 parameters to be determined. ŽSee Table 1.. 4.1. Beryllium Beryllium has many unique properties: its density is 3.5 times higher than that of Li, and also greater than the density of Mg and Ca. It has very low compressibility and electronic specific heat, its melting temperature is much higher than for the other alkaline earth elements and its Table 1 Values of the LD model parameters Metal
Be
Cr
Ni
Pd
vp f0 G0 f1 G1 v1 f2 G2 v2 f3 G3 v3 f4 G4 v4
18.51 0.081 0.034 0.050 1.345 0.227 0.101 3.441 1.325 0.388 3.862 2.882 0.260 2.916 4.532
10.75 0.168 0.046 0.199 2.622 0.162 0.257 2.280 0.820 0.952 2.436 1.981 0.848 0.155 8.552
15.92 0.089 0.027 0.058 2.793 0.099 0.150 1.362 0.491 0.100 1.962 1.585 0.714 6.018 5.847
9.72 0.331 0.008 0.763 3.487 0.392 0.098 0.512 0.496 0.523 3.815 1.284 0.479 3.387 5.753
156
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
Debye temperature is among the highest in the periodic table. Also, Be is relatively light compared to the other optically important metals, has low Z number, low thermal expansion at cryogenic temperatures, very high specific stiffness Žratio of Young’s modulus to density., high thermal conductivity and high reflectivity in the infrared. Because of these properties Be is particularly attractive for optical components intended for use in space. Unfortunately, only a few studies of the optical properties of Be have been done, and they have produced widely different results, especially in the visible range of the spectrum. The variability of the results may be due to difficulties in sample preparation w38x. Bulk samples are usually produced by high-pressure sintering of microcrystalline powder, which results in the incorporation of significant amounts of BeO in the sample. For many metals oxidation problems can be overcome by vacuum evaporation, but for Be the optical properties depend critically on the evaporation procedure. In this paper tabulated data given in the study of Arakawa et al. w38x were used. The data consists of the results of Hull w39x in the range 0.06–0.8 eV, while in the ranges 0.8–5.0 eV and 1.5–14 eV results of Hunderi et al. w40x and Seignac et al. w41x were used, respectively. In regions of overlap, higher reflectance values were favored. Fig. 10 shows the real and imaginary parts of the optical dielectric constant as a function of energy. Open circles represent tabulated data, and the solid line the LD
model. The inset shows the structure between 0.6 eV and 2.0 eV. It can be observed that the model describes the tabulated data within the experimental error. 4.2. Chromium Chromium is a body-centered cubic metal, antiferromagnetic below 312 K. The magnetic ordering, with spatial period incommensurate with the crystal lattice period, introduces structure in the spectra of the optical constants in the infrared w42x. Cr oxidizes easily, but the Cr2 O 3 overlayer is impervious to Cr or O after a thickness of ˚ is reached. Hence, there is no further oxidation about 30 A at room temperature. Such an overlayer may affect the optical constants in the visible and near infrared obtained by ellipsometry and by any other technique in the vacuum ultraviolet. There have been few studies for which the sample was both prepared and measured in UHV. Data of Bos and Lynch w43x, tabulated in Ref. w42x, resulting from Kramers-Kronig analysis of low-temperature absorptance data on bulk samples, were employed for LD model parameter determination. The measurements of Bos and Lynch extend to 13 mm in the infrared, although their Kramers-Kronig analysis extends further. A number of measurements extend to longer wavelengths, but they do not agree well with each other, and some of them show weak structures which are probably spurious. In this article the data obtained by Kirillova w44x, measured in the range
Fig. 10. Real and imaginary parts of the dielectric constant of Be versus energy Žsolid line – model, open circles – experimental data.; the inset shows the structure between 0.6 eV and 2.0 eV.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
157
Fig. 11. Real and imaginary parts of the dielectric constant of Cr versus energy Žsolid line – model, open circles – experimental data.; the inset shows the structure around 1 eV.
0.069–4.901 eV, and results of Johnson and Christy w45x, measured in the range 0.64–6.10 eV, were also employed. Fig. 11 shows the real and imaginary parts of the dielectric constant of chromium Žopen circles – experimental data, solid line – LD model.. The inset shows the structure around 1 eV. Good agreement with experimental data can be observed. 4.3. Ni The overall agreement among the published n and k spectra of Ni is good. In this article tabulated data from the
study of Lynch et al. w46x, consisting of data obtained by Vehse and Arakawa w47x in the range 3.0–25 eV, and data obtained by Lynch et al. w48x in the range 0.1–3.0 eV were employed. In Fig. 12, showing the real and imaginary parts of the optical dielectric function versus energy Žopen circles – tabulated data, solid line – LD model. good agreement between tabulated and calculated data can be observed. 4.4. Pd The optical properties of Pd have been studied under a variety of experimental conditions and by several tech-
Fig. 12. Real and imaginary parts of the dielectric constant of Ni versus energy Žsolid line – model, open circles – experimental data..
158
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159
Fig. 13. Real and imaginary parts of the dielectric constant of Pd versus energy Žsolid line – model, open circles – experimental data..
niques. Unfortunately, many results disagree with the actual magnitude of the optical structures, although the structure energies agree almost in every case. This depends on sample preparation, accuracy of the optical measurements and the method of analysis of the experimental data. In this article the data tabulated in the study of Borghesi et al. w49x were used. These data were obtained from the following works and references therein: Vehse et al. w50x in the range 6.0–12.5 eV, Johnson and Christy w45x in the range 0.6–6.0 eV and Weaver et al. w51,52x in the range 0.1–0.6 eV. In Fig. 13 excellent agreement between calculated data and tabulated values of the real and imaginary parts of the optical dielectric constant can be observed. Open circles represent experimental data, the solid line the data calculated by the LD model.
5. Conclusion In this paper the problem of applying genetic algorithms to continuous optimization was discussed. To this end, the prime concern was the possibility of applying these algorithms to a very specific and delicate problem: determination of the model parameter values by minimizing the difference between experimental and calculated data. In this case, precise and reliable location of the global minimum is imperative, since slightly different objective function values can be obtained for significantly different parameter values. Elite genetic algorithms with adaptive mutations were proposed for overcoming the difficulties in continuous optimization originating from discrete sampling of the solution space. Performance of the new algorithm was compared to that of the classical GA for three families of multiminima test functions, for up to
100 variables and different number of chromosomes in the population. In all investigated cases, EGAAM obtained extensively lower objective function values than GA. After verifying its ability to locate the global minimum on nine multiminima test functions, EGAAM was applied to model parameter determination of the optical constants of Be, Cr, Ni and Pd. For all four metals excellent agreement between calculated and experimental results is obtained. References w1x D.E. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley, Reading, 1989. w2x S. Kirkpatrick, C.D. Gelatt Jr., M.P. Vecchi, Science 220 Ž1983. 671. w3x M.W. Gutowski, J. Phys. A 27 Ž1994. 7893. w4x H. Muhlenbein, D. Schlierkamp-Voosen, Evolutionary Com¨ putation 1 Ž1993. 25. w5x K.P. Wong, Y.W. Wong, in: Proc. ANZIIS-93, Perth, Western Australia, 1993, pp. 512–516. w6x K.P. Wong, Y.W. Wong, IEE Proc. Gen. Transm. Distrib. 141 Ž1994. 507. w7x A.B. Djurisic, ˇ ´ J.M. Elazar, A.D. Rakic, ´ Optics Comm. 134 Ž1997. 407. w8x A. Chipperfield, R. Fleming, Control Computers 23 Ž1995. 88. w9x D.E. Goldberg, Complex Systems 5 Ž1991. 139. w10x K. Deb, R.B. Agrawal, Complex Systems 9 Ž1995. 115. w11x J.L. Eshelman, J.D. Schaffer, in: Proc. Foundations of GA Workshop, 1992, pp. 187–202. w12x A.D. Rakic, ´ Appl. Optics 34 Ž1995. 4755. w13x C.J. Powell, J. Opt. Soc. Am. 60 Ž1970. 78. w14x A.B. Djurisic, ˇ ´ A.D. Rakic, ´ J.M. Elazar, Phys. Rev. E 55 Ž1997. 4797. w15x A.D. Rakic, ´ J.M. Elazar, A.B. Djurisic, ˇ ´ Phys. Rev. E 52 Ž1995. 6862.
A.B. Djurisicr ˇ ´ Optics Communications 151 (1998) 147–159 w16x T.H. Han, W.W. Chang, Electron. Lett. 32 Ž1996. 1256. w17x A. Franke, A. Stendal, O. Stenzel, C. von Borczyskowski, Pure Appl. Optics 5 Ž1996. 845. w18x A.D. Rakic, ´ M.L. Majewski, J. Appl. Phys. 80 Ž1996. 5909. w19x M.K. Vai, S. Prasad, N.C. Li, F. Kai, IEEE Trans. Electron. Devices 36 Ž1989. 761. w20x A.B. Djurisic, ˇ ´ A.D. Rakic, ´ J.M. Elazar, Appl. Optics 36 Ž1997., to be published. w21x F. Curatelli, Int. J. Electronics 78 Ž1995. 435. w22x S.E. Cienawski, J.W. Ehart, S. Ranjithan, Water Resources Research 31 Ž1995. 399. w23x Q. Li, E.J. Rotwell, K.-M. Chen, D.P. Nyquist, IEEE Trans. Antennas Propag. 44 Ž1996. 198. w24x R. Vemuri, R. Vemuri, Electron. Lett. 30 Ž1994. 1270. w25x S.H. Clearwater, T. Hogg, Artificial Intelligence 81 Ž1996. 327. w26x C.B. Lucasius, M.L.M. Beckers, G. Kateman, Anal. Chim. Acta 286 Ž1994. 135. w27x R.R. Brooks, S.S. Iyengar, J. Chen, Artificial Intelligence 81 Ž1996. 327. w28x U.D. Hanebeck, G.K. Schmidt, Fuzzy Sets Systems 79 Ž1996. 59. w29x T. Back, ¨ H.P. Schwefell, in: Genetic Algorithms in Engineering and Computer Science, Wiley, Chichester, 1995, pp. 111–140. w30x D. Raynolds, J. Gonatann, Artificial Intelligence 82 Ž1996. 303. w31x F. Aluffi-Pentini, V. Parisi, F. Zirilli, J. Optim. Theory Appl. 47 Ž1985. 1. w32x A. Dekkers, E. Aarts, Mathematical Program. 50 Ž1991. 367. w33x N.W. Ashcroft, K. Sturm, Phys. Rev. B 3 Ž1971. 1898. w34x H. Ehrenreich, H.R. Philipp, B. Segall, Phys. Rev. 132 Ž1963. 1918.
w35x w36x w37x w38x
w39x w40x w41x w42x
w43x w44x w45x w46x
w47x w48x w49x
w50x w51x w52x
159
H. Ehrenreich, H.R. Philipp, Phys. Rev. 128 Ž1962. 1622. K. Sturm, N.W. Ashcroft, Phys. Rev. B 10 Ž1974. 1343. P. Drude, The Theory of Optics, Dover, New York, 1959. E.T. Arakawa, T.A. Callcott, Y.-C. Chang, in: E.D. Palik ŽEd.., Handbook of Optical Constants of Solids II, Academic Press, San Diego, CA, 1991. A.B. Hull, Perkin-ElmerrAppl. Optics Operations, Box 3115, Garden Grave, CA 92641, private communication. O. Hunderi, P. Myers, J. Phys. F 4 Ž1974. 1088. A. Seignac, S. Robin, Solid State Commun. 19 Ž1976. 343. D.W. Lynch, R.W. Hunter, in: E.D. Palik ŽEd.., Handbook of Optical Constants of Solids II, Academic Press, San Diego, CA, 1991. L.W. Bos, D.W. Lynch, Phys. Rev. B 2 Ž1970. 4567. M.M. Kirillova, M.M. Noskov, Phys. Met. Metallogr. 26 Ž1968. 189. P.B. Johnson, R.W. Christy, Phys. Rev. B 11 Ž1975. 1315. D.W. Lynch, W.R. Hunter, in: E.D. Palik ŽEd.., Handbook of Optical Constants of Solids, Academic Press, Orlando, FL, 1985, pp. 275–367. R.C. Vehse, E.T. Arakawa, Phys. Rev. 180 Ž1969. 695. D.W. Lynch, R. Rosei, J.H. Weaver, Solid State Commun. 9 Ž1971. 2195. A. Borghesi, A. Piaggi, in: E.D. Palik ŽEd.., Handbook of Optical Constants of Solids II, Academic Press, San Diego, CA, 1991. R.C. Vehse, E.T. Arakawa, M.W. Williams, Phys. Rev. B 1 Ž1970. 517. J.H. Weaver, R.L. Benbow, Phys. Rev. B 12 Ž1975. 3509. J.H. Weaver, C. Krafka, D.W. Lynch, E.E. Koch, in: Physik, Daten, Physics Data, Optical Properties of Metals, vol. 18-1, Fach-information zentrum, Karlsruhe, 1981, p. 179.