Machine based optimization using genetic algorithms in a storage ring

The genetic algorithm (GA) has been a popular technique in optimizing the design of particle accelerators. As a population based algorithm, GA requires a large number of evaluations of the objective functions, which can be time consuming. One can benefit from parallel computing with significantly reduced computing time when fulfilling the function evaluation by a numerical machine model in simulation codes. Indeed, this is the most common approach in GA applications. In this paper, instead of applying GA in the conventional numerical calculations as described above, we present a successful experimental demonstration of implementing GA in real machine based optimization. We conduct the minimization of the average vertical beam size of the SPEAR3 storage ring using GA. Beam loss rate is chosen as the sole objective function because it is inversely proportional to the vertical beam size and can be measured instantaneously in SPEAR3. The decision variables are the strengths of SPEAR3 skew quadrupoles, by varying which we can change both the betatron coupling and the vertical dispersion while searching for the minimum beam size. The results in this paper can shed light on new applications of GAs in the particle accelerator community, for example, optimizing the luminosity of a high energy collider or the injection efficiency of a diffraction limited storage ring in real time.


I. INTRODUCTION
The interest in optimizing the operation and design of particle accelerators has long been driven by various beam requirements for different scientific applications.Since modern accelerators are always complex, optimization typically deals with multi-interrelated variables and physical quantities, and a global solution sometimes cannot be reached using classical techniques such as exhaustive search or gradient based search methods.The concept of population based search allows evolutionary algorithms (EAs) to be used as very attractive tools for global optimization.Consequently, there are growing interests in applying EAs, especially genetic algorithms, for accelerator design and optimization in recent years.
Genetic algorithm (GA) is one of three major independent implementation instances of EAs along with evolution strategies [1] and evolutionary programming [2].These algorithms are broadly similar and all based on the concept of solving complex problems by mimicking the processes of Darwinian evolution.Individual solutions compete with each other continuously to discover the optimum in the whole search space.However, significant differences between genetic operators and population selections [3] make them suitable for different applications.Ever since GAs were developed by Holland in the 1970s [4], they have become popular and their applications have been extended from the classical single-objective optimization to multiobjective optimization [5].GAs were originally introduced to particle accelerator related research for a wiggler design in the early 1990s [6].Thereafter, they have been successfully applied in more areas such as designing injector systems [7], diagnosing and designing accelerating cavities [8][9][10][11][12], and optimizing beam optics design in storage rings [13][14][15].Some of these applications of GAs in the particle accelerator community are well summarized in Ref. [16].One core process of the GA is to evaluate objective functions from a given set of decision variables, i.e. knobs, to be adjusted for optimum search.In spite of various approaches used in previous efforts, all of them evaluate the objective functions based on numerical simulations using particle tracking codes such as ELEGANT [17], which has been used for optimizing the lattice of the SPEAR3 storage ring [18], or analytical models, such as an equivalent circuit of accelerating cavities [12].
In principle, when the optimization targets or their correlating parameters are measurable experimentally, it is possible to use the real machine as the function evaluator to directly measure the objective functions, rather than using a computer model.Compared to the simulation based optimization, the machine based optimization has the most accurate representation of the machine condition that includes all lattice errors.Moreover, the machine based optimization can excel in speed when the underlying physics quantities can be measured quickly in experiment.One example of such quantity is the luminosity of the high energy colliders.It can take hours to calculate the luminosity in a particle accelerator even using parallel computing.However, its measured value can normally be obtained in real time during the machine operation.In SPEAR3, we have identified a similar parameter of interest that can be measured nearly instantaneously, which in turn enables us to carry out machine based GA optimization.In this paper, we report an innovative application of GA to minimize the average vertical beam size of the SPEAR3 storage ring in real time.Smaller vertical beam sizes in a storage ring can produce higher photon brilliance.Therefore, there have been several results reported for achieving ultralow vertical beam size or emittance recently [19].To construct the GA, we use 13 skew quadrupoles in the SPEAR3 storage ring as decision variables and the average vertical beam size as the objective function.The results can serve as a proof of principle for using GA in machine based optimization.However, the average beam size cannot be directly measured in SPEAR3; instead, we optimize the beam loss rates due to their negative correlation.Thus, to minimize the vertical beam size, we need to find the maximum beam loss rate.In SPEAR3, the lattice with magnetic errors can be derived from a computer code called LOCO (Linear Optics from Close Orbits) [20] by minimizing the deviation between the model and measured orbit response matrix data.It has been proven that LOCO can give accurate predictions of parameters such as emittances, dispersions, and beta functions.Using LOCO data, optics correction can also be made, although typically several iterations are required.Since the ideal model lattice is free of transverse coupling and vertical dispersion, LOCO correction has been the regular method to minimize vertical beam size in SPEAR3.We have reduced the emittance ratio to 0.05% in the past [21].In this paper, we will compare the results using LOCO correction with those from the GA method.
The structure of this paper is as follows.In Sec.II, we derive the relationship between the beam loss rate, emittance ratio, and the average vertical beam size in SPEAR3.In Sec.III, we briefly describe the algorithm details.The optimization results using GA are presented in Sec IV.Finally, in Sec.V, we discuss the work presented in the paper and the direction of future research.

II. MINIMIZING THE VERTICAL BEAM SIZE IN SPEAR3
SPEAR3 is a third-generation, double-bend achromat electron storage ring with a circumference of 234.1 m.With the present low emittance lattice, the nominal horizontal emittance is 10 nm rad at 3 GeV with a typical emittance ratio of 0.1%.The transverse beam sizes vary along the ring due to the effect of local coupling and dispersion.By optimizing the strength of the 13 skew quadrupoles distributed around the ring, we can minimize the coupling and vertical dispersion and hence reduce the average vertical beam size.Beam size reduction leads to the shrinkage of the beam volume, which increases the Touschek beam loss in SPEAR3.In the following, we briefly describe the relationship of Touschek beam loss rates, vertical beam sizes, and emittance ratio.
In a storage ring, circulating electrons are lost due to collisions with gas molecules and electron-electron scattering inside the bunch, where the latter is known as the Touschek effect.As with most modern electron storage rings, the loss of stored beam current in SPEAR3 is mainly due to the Touschek scattering [22].In the Touschek effect, the collisions of particles with transverse oscillations can lead to the loss of particles by transferring transverse momenta to longitudinal momenta.The Touschek lifetime, τ T , can be used to characterize the electron loss rate due to Touschek scattering [23]: where N is the number of electrons in the bunch and is proportional to the beam current I; r 0 is the classical electron radius; γ is the Lorentz factor; σx;y;z represents average sizes around the ring for all directions; c Δp represents the rf bucket height; a typical range from 0.001 to 1 in a storage ring; DðξÞ is the Touschek integral, a slowly varying function.Normally σx ≫ σy , so when varying the skew quadrupoles of the storage ring, the vertical beam size changes much faster than the horizontal beam size.Neglecting the beam lengthening effect with small change in a stored beam current, σz can be treated as a constant as well.Furthermore, assuming a fixed momentum aperture and a nearly constant DðξÞ, one can derive the simple scaling law between beam loss rates dI=dt and vertical beam size as shown below: In the literature, the emittance ratio is also used to characterize the vertical beam size in the storage ring.However, depending on the definition, transverse emittance and beam size could refer to different quantities.The rms apparent emittance E is the observable quantity that varies around the ring and can be directly derived from beam size σ, dispersion D, energy spread δ, and beta function β: However, in the presence of transverse betatron coupling, the rms apparent emittance is different from the rms projected emittance ε (the second-order statistical moment represents the beam phase space) and the eigenemittance ϵ (the invariant characterizing the lattice).The influence of K. TIAN, J. SAFRANEK, AND Y. YAN Phys.Rev. ST Accel.Beams 17, 020703 (2014) 020703-2 betatron coupling on these three transverse emittances has been described using the resonance driving term formalism in Ref. [19].For a coupling free lattice, these three quantities coincide, E ¼ ε ¼ ϵ.When coupling sources exist, the apparent emittance oscillates around the ring and its amplitude generally increases with the coupling.This effect is especially important for the vertical emittance.The projected emittance stays constant in coupling free regions and exhibits abrupt jumps at locations of coupling sources.However, with enough observation points, the mean projected emittance is equal to the mean apparent emittance.Therefore, the emittance ratio can be calculated from both of them, i.e. εy =ε x ≈ hE y i=hE x i.When the average vertical beam size is minimized in a storage ring, the vertical dispersion is also small.Hence, one can claim that minimizing vertical beam size is equivalent to minimizing the emittance ratio.As a matter of fact, the emittance ratio has been widely used interchangeably with vertical beam size in the literature on reducing vertical beam size in storage rings.In SPEAR3, direct validation of Eq. ( 2) can be demonstrated by experimentally measured data.During the experiment, 100 mA is filled in SPEAR3 with the normal filling pattern of 280 bunches.We achieve 20 conditions in the ring by varying all skew quads in a fixed step size.The current decay during the experiment is about 4%, hence is negligible.For each condition, we measure the vertical beam size at 1 Hz with an x-ray pinhole camera [21].The beam loss rate is recorded at 1 Hz by a beam loss monitor (BLM), a 2-inch NaI scintillator and photomultiplier tube [24], installed downstream of the SPEAR3 horizontal beam scraper, which is inserted to make it the minimum physical aperture in the ring.Both beam size and beam loss data are averaged for 2 min to reduce noise.The global beam loss is also acquired by recording the current decay during this 2-min period using a dc current transformer (DCCT).The current droop during the entire period of the experiment is less than 4%; thus, the stored current can be considered as a constant for the analysis in Eq. (2).Touschek scattered electrons tend to be lost at locations with a large dispersion area or smaller horizontal aperture.In SPEAR3, the horizontal scraper is parked beyond −33 mm during operation so that it does not interfere with the stored beam.A large amount of the Touschek beam loss during operation occurs at a high dispersion section that is designed intentionally as a beam loss point.Another BLM was installed there to monitor the beam loss rate.By inserting the horizontal scraper closer to the beam, more and more electrons are lost at the scraper.As shown in Fig. 1, at −15 mm, the scraper replaces the high dispersion section as the dominant place for the Touschek beam loss.In order to capture as much loss as possible, we insert the scraper to −6 mm for a loss rate above 16; 000 counts= sec.We believe that this insures that we have most of the Touschek beam loss captured at the scraper and recorded by the scraper BLM.
Considering the demagnification of the beam from the optics setup and pixel resolution of the CCD sensor in the digital camera, the effective resolution for imaging the beam with the pinhole camera is about 6 μm=pixel.This resolution is fine for the beam with relatively larger vertical beam size, but can be less accurate for vertical beam size measurement when the beam size is small.In addition, the measurement of beam size at one location of the ring is not enough to determine the Touschek life time.Therefore it can only be used as verification to the change of the lifetime.The measured results are shown in Fig. 2. Agreeing with the scaling law in Eq. ( 2), the beam loss rate is inversely proportional to the vertical beam size.The global beam loss, measured from the stored DCCT, has more pronounced noise, but it still serves as a reasonable cross-check to our assumption that the BLM near the scraper captures most of the beam loss.As the data shows, the vertical beam size measured at the x-ray pinhole camera tends to overestimate the average beam size in SPEAR3.
As previously discussed, the emittance ratio and vertical beam size can be calculated from LOCO analysis.Sampling 5 out of the 20 different cases, we measured the response matrix of the ring for each, and then fit the model to the measured data using LOCO.Emittance ratio and average vertical beam size then can be calculated from the fitted model.The results are shown in Fig. 2, where the quadraticlike relationship between the emittance ratio and beam loss rates agrees with Eqs. ( 2) and (3).In addition to reduction of the average vertical beam size, decreased off energy dynamic aperture through resonance excitation could also be the cause for increasing beam loss rate while adjusting the skew quadrupoles.In such a case, we could fail to minimize the vertical beam size by maximizing the beam loss rate.However, the monotonic relationship of the measured beam loss monitor data, the vertical beam size, and the emittance ratio in Fig. 2 helps to ease this concern.Furthermore, the small beam size we achieved by maximizing the measured beam loss further eliminates concerns about resonant effects.One can also verify the absence of resonance excitation by repeating the measurement with smaller rf acceptance by reducing rf cavity gap voltage.But as it has been well proved that the resonance excitation is not present while varying the skew quadrupoles in SPEAR3, this is not necessary in our study.Therefore, based on the above discussion, we can conclude that, by setting the horizontal scraper to −6 mm, finding the maximum of the beam loss rates measured by the BLM near the horizontal scraper is equivalent to minimizing the vertical beam size or the emittance ratio.

III. ALGORITHM DESCRIPTION
GAs start an optimization process by initializing a random population.Each individual in the population is a chromosome, which at least contains a vector in the hyper dimensional space of decision variables that represent a solution to the problem and the values of corresponding objective functions.Additional information including the rank of the objective functions, fitness, and crowding distance is usually encoded into a chromosome depending on specific problems.Under some selection strategies, some of the chromosomes are chosen as parents to generate the offspring via genetic operations.Values of objective functions are provided by the function evaluator after each offspring is generated.Depending on algorithms, the population size of each generation can be fixed or varied.The parents can be used to form the new generation with offspring, the elitist approach, or always be discarded, the nonelitist approach.After a new generation is created, it repeats the same process for generating the next generation and will not stop until meeting the stopping criterion.In this section, we will detail the formation of a GA based technique for maximizing the beam loss rate of SPEAR3 in real time.

A. Objective functions and decision variables
As discussed in the previous section, the normalized beam loss rate measured by the BLM next to the scraper is the single objective function in our GA formation.The currents, which represent the strengths of the skew quadrupoles in SPEAR3, are the independent decision variables.After varying the currents of the 13 skew quadrupoles, the objective function is evaluated from the direct measurement of the BLM.Distinct from most other applications of GA, the accelerator serves as the function evaluator, instead of a numerical or analytical model.This approach has the obvious advantage of providing the exact evaluation without any approximation, but it is usually difficult for two reasons.First, the experimental evaluation is usually subject to noise, and thus can be inconsistent.Second, it can be too slow to implement a machine based GA.As a global optimizing algorithm, a large number of evaluations are required before reaching the vicinity of the optimum in the solution space.Unlike parallel computing capability in many simulation codes, the machine normally can only evaluate one variable set at a time.Additional time has to be spent if magnet ramping or standardization is involved when changing the decision variables.However, with proper diagnostics in SPEAR3, we can overcome these problems due to the following facts.The SPEAR3 skew quadrupoles are powered by high precision power supplies featured with fast switching and the field variations of the skew quadrupoles are small enough to neglect the hysteresis effects when searching with GA.The data acquisition from the BLM is also nearly instantaneous.
Using the objective function, i.e. the beam loss rate normalized to the square of stored current, decisionvariables, and the rank of the objective function in the current generation, we form each individual in the population, i.e. the chromosomes.A pool of parents is selected in each generation from these individuals to create children through the process of genetic operations.

B. Genetic operations
Except the first generation, all sequential generations are generated via the process of genetic operations.We use two of the most popular genetic operators: real-coded simulated binary crossover [16,25] and polynomial mutation [16,26].In each genetic operation, one of the two operators is chosen randomly but conforms to a predefined ratio.During a crossover, two parents are picked to create two children, while one child is generated from one parent in the case of mutation.Once a child is created, its corresponding objective function is evaluated.Eventually, an offspring population is generated.Crossover and mutation are governed by the user-configurable non-negative tuning parameters η c and η m , respectively.A more detailed K. TIAN, J. SAFRANEK, AND Y. YAN Phys.Rev. ST Accel.Beams 17, 020703 (2014) 020703-4 discussion of these two parameters can be found in the literature such as Ref. [16].In short, these tuning parameters control the probability density function of the likeness between parents and children.For mutation, a smaller η m represents less probability of having a similar child, which in turn provides a global search to the optimum regardless the parent solution.On the other hand, with a bigger η m , it is very likely that the child only varies slightly from the parent and the operation is conducting a local search of optimum around the parent.The behavior of η c is quite similar to η m : the children tend to be close to one of the parents with large η c or be randomly generated with small η c .

C. Replacement, reevaluation, and stopping criteria
To ensure an elitist approach, the current population is replaced by the best solutions chosen from both the offspring and the older generation.To maintain a fixed population size, the remaining solutions are discarded.As the objective functions are measured directly from the SPEAR3 machine, the results may change over time due to variation of machine condition.Therefore, we reevaluate the surviving solutions from all previous generations every 10 generations.Limited by the machine time available for the experiment, we run the algorithm as long as possible, so normally we stop the program manually after a certain amount of time.

IV. RESULTS
We choose the population size of 120 in each generation for reasonably big sample size and relatively short time for generating the whole population during the experiment.As a result, it takes less than 3 min to generate one generation.Figure 3 shows the results with 211 generations of GA optimization.To reduce the effect of stored beam current decay, we refilled SPEAR3 twice during this experiment as shown by the normalized stored beam current in Fig. 3(a) (scaled stored current shown as red curve).The optimization was paused during the refill and restarted by loading the dumped data after the fill.The total running time of the algorithm is a bit over 9 hr.Overall, the algorithm behaved well.The normalized beam loss rates in Fig. 3(a) grow steadily for the first 150 generations, and then start to converge.In Fig. 3(b), we also plot the decision variables of select generations with 120 individuals apiece for better understanding of the progress of optimum search.Starting from the first generation (in red), the decision variables or the individuals are generated randomly with a nearly uniform distribution within the boundaries of the hyperspace.The solutions start to cluster at several regions rather than spread out in the whole hyperspace in the 6th and 11th generation.It appears that the final region of the solution is found in the 156th generation.Thereafter, new solutions stop drifting away from this area.This is consistent with the results of the beam loss rate shown in Fig. 3(a).The skew quadrupoles currents of the best solutions in these select generations are shown in Fig. 3(c).In Fig. 3(d), we compare the skew quadrupoles currents of the GA solution and the LOCO solution.The beam loss rates for the best solution found using GA are compared with the solution found using LOCO in Table I.The average vertical beam size and emittance ratio of the ring with these two solutions are also calculated from the fitted model using LOCO.
The normalized beam loss rate measured with the GA solution is increased by 17.9% from that with the LOCO solution.According to the scaling law in Eq. ( 2), this is translated to the reduction of vertical beam size by 15.18%.On the other hand, the calculated values of average vertical beam size using LOCO fitting show a 10.99% reduction of vertical beam size for the GA solution.The discrepancy comes from the accuracy of the fitting in LOCO, loss rates measurement, and the assumption of the solely Touschek effect when deriving the scaling law in Eq. ( 2).
With the GA based optimization, we have found a solution that outperforms the LOCO solution.But the LOCO correction only takes about 15 min compared to the 9-hr run of the GA optimization.We can significantly reduce the time of GA optimization to be within 2 hr by including the LOCO solution in the first generation.The performance and speed of genetic algorithms highly depends on specific problems and can be adjusted with mutation and crossover tuning parameters.Although we have not conducted a thorough study of setting the most appropriate tuning parameters for our problem, we have programmed the code to dynamically adjust these factors according to the diversity of the population in the new generation.It appears that the optimization progresses faster by promoting a more global search with relatively small mutation tuning parameters in early generations.When the population starts to cluster toward the optimum region, it helps to save time by using a large mutation tuning parameter or by shrinking the search space.Without setting the tuning parameters properly, we have observed early convergence of the program, which causes the failure of meeting our goal within a reasonable amount of machine running time.

V. DISCUSSION AND CONCLUSION
Genetic algorithms are believed to be especially suitable for problems with high complexity where traditional gradient based search methods normally fail to optimize.As the storage ring lattice is well designed, the coupling optimization of the ring tends to be a well behaved problem.This is evident from the final solution: out of 13 skew quadrupoles, only 3 are required to be set above 5 A, while most of the others are near zero.This fact weakens the advantage of using GA based optimization.A thorough simulation study has been carried out in SPEAR3 in parallel to comparing different algorithms for real time optimization of the beam loss rate [27].After comparing other techniques including the Nelder-Mead simplex method, Powell's conjugate direction method, and a modified conjugate direction method, it is found that for this particular problem, GA lacks speed for real time optimization, and its performance is sensitive to the random noise of the measured data.Nevertheless, with machine based GA, we are able to find a good solution regardless the time spent.Also, we have more confidence in its global validation.In addition, one should note that machine based GA may show advantage in speed over the traditional gradient based techniques when optimizing problems with more decision variables.As long as the corresponding hardware can be set roughly simultaneously, the time cost by machine based GA is independent of the number of decision variables.However, most traditional algorithms are scaled with the number of decision variables in high order.Thus, machine based GA can be more valuable for large machines.
Genetic algorithms are also usually considered to be good candidates for multiple objective optimizations and not appealing when dealing with single objective optimization such as decision making problems.However, unlike a definite lattice model in simulation, the real machine lattice varies with time, which makes the accurate prediction of objective functions difficult even for a single objective function, for example, the luminosity of a high energy collider.As a result, it is normally impractical to find the exact best solution.Therefore, it is required the one searches for a better solution from time to time with a certain algorithm, where the GA technique can be a powerful candidate.Our experimental results illustrate that with proper hardware, fast data acquisition capability, and solid underlying physics, GA can be successfully used for problems of machine based optimization, even beyond the particle accelerator community.We believe this is the most important contribution of this work.However, the GA techniques we used are far from being refined.In future study, we will focus on improvements in speed and robustness.One possible approach is to create a hybrid algorithm that combines both GA and one of the traditional techniques for a fast local search.When blending the two algorithms, it is challenging to maintain each one's original advantage, which requires thorough study.

FIG. 3 (
FIG. 3 (color online).(a) Normalized beam loss rates at the scraper during the experiment at 1 Hz update rate; (b) population for selected generations (red: 1st generation; green: 6th generation; blue: 11th generation; cyan: 156th generation; black: 211th generation); (c) the best solutions at selected generations; (d) comparison of the GA and LOCO solutions.

TABLE I .
Comparison of optimized solutions from LOCO and GA.