A convenient and robust procedure for Mössbauer analysis based on Monte Carlo method is described. The method uses simulated annealing approach to find the optimum Mössbauer parameters in the Lorentzian profile as initial values for the Monte Carlo search program. A succession of solutions to the function describing the spectrum is then randomly generated until the solution with the minimum chi-square with respect to the experimental data is reached. The result having the reduced chi-square close to 1 shows the validity of the method.Opisujemo pogodan i snažan postupak za analizu Mössbauerovih spektara zasnovan na Monte Carlo metodi. Primjenjuje se pristup simuliranog opuštanja za nalaženje najpovoljnijih parametara Lorentzovih profila Mossbauerovih spektara koji su početne vrijednosti za Monte Carlo program traženja. Zatim se nasumice stvara niz rješenja za funkciju koja opisuje spektar dok se ne postigne rješenje koje je u najboljem skladu s eksperimentalnim podacima. Postignuti ishod je u dobrom skladu s mjerenim spektrom, što pokazuje vrijednost metode

Ali, Abdelilah Musa

Hsia, Yuanfu

Huang, Hongbo

Li, Jun

Shi, Hairong

Sun, Xuefeng

English

Hrčak - Portal of scientific journals of Croatia

ISSN1330–0016CODEN FIZBE7USE OF MONTE CARLO METHOD WITH SIMULATED ANNEALINGALGORITHM IN THE ANALYSIS OF MÖSSBAUER SPECTRAABDELIILAH MUSA ALI1, YUANFU HSIA, HAIRONG SHI, HONGBO HUANG,XUEFENG SUN and JUN LINDepartment of Physics, Nanjing University, Nanjing, ChinaReceived 29 March 2001;Accepted 24 September 2001 Online 25 May 2002A convenient and robust procedure for Mössbauer analysis based on Monte Carlomethod is described. The method uses simulated annealing approach to find theoptimum Mössbauer parameters in the Lorentzian profile as initial values for theMonte Carlo search program. A succession of solutions to the function describ-ing the spectrum is then randomly generated until the solution with the minimumchi-square with respect to the experimental data is reached. The result having thereduced chi-square close to 1 shows the validity of the method.PACS numbers: 05.10.Ln, 33.45.+x UDC 539.166Keywords: Monte Carlo method, simulated annealing, global optimization, Mössbauerspectra1. IntroductionThe method presented in this paper concerns the Monte Carlo computation ofa given Mössbauer spectrum and compares it with the experimental spectrum tofind the best fit. For the theoretical computation, we express the spectrum in theLorentzian line shape. For each line component, I, we assume [1]I =Aj1 + [(xi − Pj)/Bj ]2 , (1)where Aj , Pj and Bj are the height, the position and the half width at half maxi-mum, respectively, for the jth peak, and xi is the ith channel. Equation (1) leads to1Corresponding author, Tel. & Fax: +86-25-3594234, E-mail: Abdulilah65@hotmail.comFIZIKA B 11 (2002) 1, 19–26 19ali et al.: use of monte carlo method with simulated . . .the expression for the total spectrum of n peaks and a count channel yi such thatyi = E + Fxi + Gx2i +n∑j=1Aj1 + [(xi − Pj)/Bj ]2 , (2)where E, F and G are the parameters for non-linear background. The quadraticbackground is of instrumental origin, mainly due to the periodical change of thedetector solid angle during the movement of the transducer. The summation rep-resents the sum of n Lorentzians.The method proceeds iteratively starting from the solution to Eq. (2) in termsof optimum values of the peaks parameters, which are provided by simulated an-nealing algorithm (SAA), and generates a succession of solutions tending to theone that is consistent with the experimental data. The concepts of annealing ina combinatorial optimization were first introduced, in the early 1980’s, by Kirk-patrick, Gelatt and Vecchi [2], and independently by Cerny [3]. These concepts arebased on a strong analogy between the physical annealing process of solids andthe problem of solving large combinatorial optimization problems. Pursuing thisanalogy, SAA is introduced as a stochastic method to avoid getting stuck in local,non-global minima, when searching for global minima [4]. A strong feature of thisalgorithm is that it finds high-quality solutions, which do not strongly depend onthe choice of the initial solution, i. e., the algorithm finds the global minima by thecurrently tested point rather than the point which has initially been chosen. Thisfeature enables the Monte Carlo search program to start with good initial valuesand thus the user can avoid conducting several trials, which have been necessaryfor setting suitable initial values.In this paper, we use SAA with the modifications made by A. Corana et al.[5] and Goffe et al. [6] to obtain the parameters that will give a “best fit” to theexperimental observations. We define the “best fit” in a least-squared sense, andlook for those parameters which minimize the function χ2, defined byχ2 =m∑i=11yi[y(x) − yi]2 , (3)where y(x) is the experimental value, yi is the calculated (theoretical) data, i isnumber of experimental measurements and m is the number of channels in the datacollection system.In each step of the algorithm search for the global minima, a trial point israndomly generated. The search process follows the path (the direction along whichthe trial points are generated) and moves both uphill and downhill to accept orreject the point. Any downhill step is accepted and the process is repeated fromthis new point. The uphill step is also accepted with a certain probability known as“Metropolis criterion” [7–8]. In this way, the search procedure ‘escapes’ from thelocal minima to the global (or near-global) minima.20 FIZIKA B 11 (2002) 1, 19–26ali et al.: use of monte carlo method with simulated . . .2. Description of the procedure2.1. Simulated annealing algorithmSAA starts at some ‘high’ value of the control parameter, C, and a sequenceof solutions Si, with i = 1, 2, 3, . . ., is generated until their average value of thecost function, χ2, reaches a stable value. The step length is then adjusted and thebest solution reached for this value of C is recorded as Sopt. C is then carefullyreduced until the termination criteria are met. The following steps explain how thealgorithm proceeds:(i) Set the initial value of the control parameter (C0), the initial value of steplength vector (Ls), the number of iterations before the control parameterreduction (Nitr), the number of cycles before the step length adjustment (Ns),the value for termination criterion (εs), the control parameter reduction factor(α) and the starting values of the height (Aj), the peak position (Pj) and thehalf width at half maximum (Bj) of each peak in the spectrum.(ii) Calculate the baseline.(iii) Find the solution, Sj , to Eq. (1) in terms of Aj , Pj and Bj . Incorporatingthe baseline (as in Eq. (2)), find the corresponding χ2 value according to Eq.(3).(iv) Set Sopt = Sj and χ2opt = χ2j .(v) Generate a new solution, S′j , with randomly selected values for the peakparameters by setting A = Aj + RLs, B = Bj + RLs and P = Pj + RLs,where R is a random number generated in the range [-1,1] by a pseudorandomnumber generator, and then calculate χ2 ′j .(vi) If χ2 ′ ≤ χ2j , accept the solution S′j . Set Sj = S′j and χ2j = χ2 ′j . If χ2 ′j > χ2j ,use Metropolis criterion to decide on the acceptance or rejection of S′j . Thatis, accept the solution with the probability p given byp = exp(χ2j − χ2 ′jC)> R .If the solution is accepted, set Sj = S′j and χ2j = χ2 ′j . Otherwise, this indicatesthat the solution lies outside the χ2 domain, reject it and choose another onefor the function evaluation.(vii) If χ2 ′ < χ2opt, set Sopt = S′j and χ2opt = χ2 ′j .(viii) To widely sample the function, adjust Ls so that, approximately, half of theevaluation will be accepted.FIZIKA B 11 (2002) 1, 19–26 21ali et al.: use of monte carlo method with simulated . . .(ix) Check whether the termination criteria are met. That is, if the final χ2j valuesof the last iterations for each control parameter reduction differ from χ2j atthe current control parameter by less than εs, and the final χ2j value at thecurrent control parameter also differs from the current optimum χ2j value byless than εs, report the optimized solution, Sopt, and its corresponding Aopt,Popt and Bopt. If the termination criteria are not met, prepare for anotherloop by staring on the best current optima. Reduce C by an appropriatefactor, α, and go to step (v) until convergence.2.2. Monte Carlo methodIf the values obtained from SAA are not optimal, as they usually are, good near-optimal values will then be obtained [5]. Whatever the case, generation of furthersolutions is required in order to get the best one. The following steps explain thispoint as follows:(i) Use the optimal values Aopt, Popt and Bopt resulting from SAA to calculatethe theoretical spectra (recorded as yk) according to Eq. (2) and calculatethe corresponding χ2k.(ii) Generate a sequence of random solutions, yk+l (with l = 1, 2, 3, . . .), aroundyi. Calculate the χ2k+l value corresponding to each solution.(iii) Sort the theoretical data according to χ2 values. Stop the program at a valueεm, which is low enough that no further useful improvement can be expected,that is, when∣∣∣∣χ2k+l − χ2k+l−1χ2k+l∣∣∣∣ ≤ εm.(iv) Print the theoretical data, which belongs to the best χ2 value (the smallestone).3. Results and discussionThe method was applied to a simple Mössbauer spectrum of α-Fe, where thepeaks are well resolved and their number is known: nine peaks appear in thespectrum presented in this paper. The plot of the computed spectrum and theexperimental one is shown in Fig. 1. For the program execution, the followingvalues of input parameters were used:C = 100, α = 0.85, εs = 10−5, εm = 10−6, Nitr = 5, Ns = 20.The parameter εs was set to ensure that the solution had converged to the globalminima, while εm was chosen to get the difference of successive χ2 values that issufficient to stop the calculation. As the search algorithm proceeds, the searched22 FIZIKA B 11 (2002) 1, 19–26ali et al.: use of monte carlo method with simulated . . .area is restricted to a suitable range by setting upper and lower bound values forFig. 1. Plotted output of Monte Carlo sum to nine Lorentzians (line) and experi-mental spectrum (squares).TABLE 1. Result of simulated annealing.Optimal chi-squired value (reduced) 1.43Number of function evaluations 162001Number of accepted evaluations 81020Final value of C 0.304E-06Step length O(E-4)Number of out of bound evaluations 0the peak parameters. This restriction forces the algorithm to search only the re-gion of interest so that very large values (which may cause floating-point errors) orvery small values (which may not allow a complete search of the function) shouldnot be run through the statistical function χ2. Table 1 shows the result obtainedby simulated annealing. Almost half of the function evaluations is accepted. Thisensures that the step length vector is correctly adjusted and the function is widelyinvestigated. Before the search procedure reaches the optimal values, a large num-ber of intermediate results with different step length values are executed to give aFIZIKA B 11 (2002) 1, 19–26 23ali et al.: use of monte carlo method with simulated . . .valuable information about χ2. Table 2 shows two steps from these intermediateresults (from more than 150); the first step and the last step before using simu-lated annealing subroutine. The step length starts with a large value to give anapproximate estimation of the function χ2, and accordingly, most of the moves areaccepted. As the control parameter C declines, the step length decreases and thecalculation is focused on the most promising area.TABLE 2. The intermediate result from SAA.Intermediate resultsFirst step Last stepCurrent value of C 100.0 0.3419E-6Min function value so far 3.54 1.463Total moves 900 900Downhill 453 363Accepted uphill 396 66Rejected uphill 51 471Step length 4.0 O(E-4)-O(E-5)Trails out of bounds 0 0The optimal solution given by SAA is quite useful for this method since only fewsolutions (21 in our case) need to be performed by Monte Carlo search for the bestcomputed spectrum. Perhaps this is also useful for some other fitting procedureswhere good initial values are recommended. The constraints among the parameterscan be easily applied without affecting smooth running of the program. The abilityof the SAA to identify corner solutions for functions that do not exist in someregions [6] gives the possibility for finding the peak(s), if any, that was (were) notincluded in the computations.The major limitation of this method is that it is costly in terms of functionevaluations, i.e., to reach the global optima, C should be decreased slowly. This isa lengthy, time consuming process. However, the execution time can be minimizedby an appropriate choice of α and Nitr prior to the calculation. Coupled with therapid increase of computer power, this limitation will no longer affect the methodand its applicability with regards to the computation time.4. Summary and conclusionWe have applied the Monte Carlo method to analyze Mössbauer spectra. Theapplication shows that SAA can be effectively used to Mössbauer analysis and offersa viable way of optimizing functions fit to Mössbauer data. The method is conve-nient for its simplicity and general applicability because no specific informationabout the function to be optimized is needed prior to application. Coupled withthe fact that SAA is capable to optimize difficult functions (with large number of24 FIZIKA B 11 (2002) 1, 19–26ali et al.: use of monte carlo method with simulated . . .parameters), this procedure may also be used to analyze other profiles that describ-ing Mössbauer spectra. The source code is written in Fortran and is available fromthe authors.References[1] U. Gonser, Topics in Applied Physics, Vol. V, Mössbauer Spectroscopy, Springer-Verlag,Berlin (1975) p.208.[2] S. Kirkpatrick, C. D. Gelatt and M. P. Vecchi, Science 220 (1983) 671.[3] V. Cerny, J. Optimization Theory and Applications 45 (1985) 41.[4] A. Dekkers and E. Arts, Mathematical Programming 50 (1991) 367-393.[5] A. Corana, M. Marchsi, C. Martini. S. Ridella, ACM Transaction on MathematicalSoftware 13 (1987) 262.[6] W. L. Goffe, G. D. Ferrier and J. Rogers, J. Econometrics 60 (1994) 65.[7] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth and A. H. Teller, J. Chem. Phys.21 (6) (1953) 1087.[8] E. Arts and J. Korst, Simulated Annealing and Boltzmann Mechanics, Anchor Press,Great Britain (1989) p.14.FIZIKA B 11 (2002) 1, 19–26 25ali et al.: use of monte carlo method with simulated . . .PRIMJENA MONTE CARLO METODE S ALGORITMOM SIMULIRANOGOPUŠTANJA U ANALIZI MÖSSBAUEROVIH SPEKTARAOpisujemo pogodan i snažan postupak za analizu Mössbauerovih spektara zas-novan na Monte Carlo metodi. Primjenjuje se pristup simuliranog opuštanja zanalaženje najpovoljnijih parametara Lorentzovih profila Mössbauerovih spektarakoji su početne vrijednosti za Monte Carlo program traženja. Zatim se nasumicestvara niz rješenja za funkciju koja opisuje spektar dok se ne postigne rješenje kojeje u najboljem skladu s eksperimentalnim podacima. Postignuti ishod je u dobromskladu s mjerenim spektrom, što pokazuje vrijednost metode.26 FIZIKA B 11 (2002) 1, 19–26

Primjena Monte Carlo metode s algoritmom simuliranog opuštanja u analizi Mossbauerovih spektara

Abstract

Similar works

Full text

Available Versions

Hrčak - Portal of scientific journals of Croatia