Methods and Instruments for Beam Lines Global Optimization

PHYSCON 2011, León, Spain, September, 5–September, 8 2011 METHODS AND INSTRUMENTS FOR BEAM LINES GLOBAL OPTIMIZATION ∗ Serge Andrianov Eugeny Podzyvalov Andrew Ivanov Applied Mathematics and Control Processes, Saint Petersburg State University, Russia. [email protected] Applied Mathematics and Control Processes, Saint Petersburg State University, Russia. podzyvalov [email protected] Applied Mathematics and Control Processes, Saint Petersburg State University, Russia. [email protected] Abstract In the paper some methods and instruments for beam lines global optimization are discussed. It is known that modern beam machines (from short beam lines to circular accelerators) are based on very complex control systems, which ensure desired beam characteristics. The necessary demands formulated in the form of a functional set. These problems are referred to the so-called NP-complete problems and require searching for a possible solutions manifold in multi-objective parameters spaces. In the paper we discuss a optimization problem for a special class of optimal beam lines — the class of ion-optical systems (such as micro- and nanoprobe systems, matching channels as so on). For successful decision of this problem we suggest to use methods and instruments for global optimization problem solution. This approach is based on genetic algorithms and neural networks tools. Key words Beam control, nonlinear dynamics, applications, nonlinear control, global optimization. 1 Introduction It is well known that optimization tools are needed in every step of any accelerator project: from the design to corresponding beam line start-up. In different steps the project require different optimization methods and technologies, for example, the initial (design) step of any project is based on a scanning of appropriate parameters set. The next step demands some local optimizers giving us more adequate during operations to choose the next best control parameters set. It is obviously that the optimization procedure effectiveness depends on used mathematical methods and ∗ THE WORK IS SUPPORTED BY FEDERAL TARGETED PROGRAMME “SCIENTIFIC AND SCIENTIFICPEDAGOGICAL PERSONNEL OF THE INNOVATIVE RUSSIA IN 2009-2013” (GOVERNMENTAL CONTRACT NO. P 793) computer modeling instruments. The modern electrophysical facilities (from enough small machines such as ion-optical systems (IOS) up to large linear and circular accelerators) are considered as multi-parametric control systems with complex organized set of requirements (see, for example, [Andrianov, 2004]). The modern beam line represent complex dynamical systems, described by nonlinear motion equations for beam particles. The choice of appropriate mathematical tools for both dynamics simulations and optimization procedures plays a significant role. At present all problems of optimal solution searching processes can be separated in two follows stages: systems design with the given properties; control system optimization using corresponding set of functionals. In the paper main attention is focused on the first stage. Modern control systems of beam lines even in the case of small control elements number (for instance, limited ten) present itself complex adjustments. This feature brings the necessary to create the special methods and algorithm providing efficient decision of the problem. The analysis of the constructive features of the modern beam lines allows to formulate the following main positions: a nonlinear nature of particles motion equations; a collective character of control object (the beam is considered as particles ensemble); a multi-objective character of optimization problems; the presence of some representative set of optimal solutions; it is necessary to consider the feasibility of optimal solutions (the tolerances problem). The own experience of paper authors (see e. g. [Andrianov, 2009]) and a review of existing publications [Yang, 2008], urge us to use for solution of beam line design as methods of global optimization based on corresponding mathematical methods and information and computational technologies. In the next this approach we name as a global optimization concept. As an additional demand we also mention the requirement fitting parameters tolerances. The main advantage of the presented approach is the possibility of finding a global extremum in the multi-objective parameters spaces [Michalewicz, 1996]. The presented concept also includes the analysis of local extremums, which may be more resistant to parameters values change, and therefore easier and cheaper realizable in practice. Considering the global optimization problems in accelerator theory, we should mention the papers [Gao, 2010; Yang, 2008], dedicated to the application of the genetic algorithms ideology for the solving problems of beam physics. In the article [Gao, 2010], for example, consider the problem of finding the optimal parameters of a linear system subject to obtain low beam emittance. However, there is no considered the impact that additional functionals may have included on the constructed set of Pareto optimal solutions (eg., restrictions on the beam features on the target). This will require the researcher to conduct a further analysis of their mutual influence, which is not always trivial. It is worth mentioning the need for a mechanism of global optimization in conjunction with modern computer simulation. Thus, it is important to properly assess the possibility of a development environment and existing methodologies to create appropriate models and corresponding program complexes. 2 Mathematical model of beam dynamics The beam physics is based on four basic points: first of all we have to consider an evolution of ensemble of particles (one has to take account of a collective character of beams); secondly — the motion of particles are described by nonlinear equations (the Newton–Lorentz equations) [Chao, 1999]; in the third place the control system in beam physics is designed as an optimal structure of control system (generated by electromagnetical steering elements) with very many control parameters; fourthly — the corresponding optimal control problems have two features: 1) multi-parameter and nonlinear character of control system; 2) conflicting character of a set of objective functions. Such problems solutions always begin with the construction of a linear model, consisting of quadrupole lenses, dipole magnets and free spaces between them. However, this approach could be considered taking into account the nonlinear effects and the effect of self-charge, because the mathematical model of IOS is based on the matrix formalism, which has a certain degree of universality. In addition to accounting for various restrictions on the control parameters of the system is important to detect parameter values sets that are optimal in terms of system implementing. Feature of control systems of beam lines is a set of particles, which can be regarded with equations of motion for each particle individually or use the approaches describe the entire ensemble of particles in general. 2.1 Beam Particle Motion Equations It is known that each particle beam motion is described using the well-known Newton–Lorentz equation [Chao, 1999]: dX = F(X, s), ds (1) where F(X, s) = F(X, E, B, s) is a Newton–Lorentz force, and E — a electrical field intensity vector, B — a magnetic induction vector. These vectors are described in an appropriate curvilinear system of coordinates, associated with some reference orbit [Chao, 1999]. The equations of motion (1) in general form can be written as: ∞ dX ∑ 1k P (E, B, s)X[k] , = ds (2) k=0 where X[k] = X ⊗ . . . ⊗ X. The symbol “⊗” means Kronecker product and X = {x, x′ , y, y ′ }∗ in two dimensional case of transversal phase coordinates or in ∗ the general case we have X = (x1 , x2 , . . . , x2n ) , x, y — transverse coordinates to the reference orbit, the symbol “*” means transparent and x′ = dx/ds. Similarly, we can describe the equation for the six-phase vector which also includes the longitudinal phase coordinates. Usually, on the first step a researcher considers so called “linear approximation” of motion equations: dX = P11 (s)X, X(s) = R11 (s|s0 )X0 , ds this case presented, for example, in [Andrianov, 2004]. The solutions of equations (2) can be found in the following form: X(s) = ∞ ∑ [k] R1k (E, B; s|s0 )X0 , k=0 where X0 = X(s0 ) — “the initial value” of the phase vector. The program implements a similar problem has the opportunity to represent the system structure in graphical form. 2.2 Representation of the Particle Beam There are several ways to write the equations of motion of a particle beam as an ensemble. It is worth mentioning three of them: 1) Single-Particle Description. Set ( of particles is rep) resented as a matrix MN (s) = X1 (s), . . . , XN (s) , where Xi (s) — current phase coordinates vector (at the moment s), N — number of particles, dim X = 2n, where 2n — dimension of phase space and MN is a matrix describing a particle ensemble M(s) including N particles. 2) Envelope Formalism Description. Another way to view a description of the beam is in the terms of envelope S(s) — matrix dimension n×n. The most popular form of this matrix is described in terms of rms (root– mean–square) values [Chao, 1999]: S(s) = ∫ ⃗ s)X(s)X∗ (s)dX, f (X, M(s) where M(s) — a current phase manifold occupied by ⃗ s) — a distribution function. beam particles, f (X, The transverse manifold occupied by beam particles can described by S(s)–matrix (i. e. rms-envelope matrix). Besides the rms-envelope matrix we can introduce another envelope matrices (see, for example, [Andrianov, 2004]). For all forms of envelope matrices one can write the following form of the matrix dynamical equation: S(s) = R(s|s0 )S0 R∗ (s|s0 ). (3) 3) Phase Space Distribution Function. Another way of representing the beam line is to describe an ensemble of particles in the shape of the distribution function f0 (X) in phase space, satisfying the Vlasov system equations: ∂f (X, s) ∂f (X, s) + FL (X, E, B, s) = 0, ∂s ∂X∗ where x1 , . . . , x2n are phase coordinates of a particle and Tik1 ...k2n — are tensor components representing partial derivations in the Taylor expansion [Chao, 1999] and k = (k1 , . . . , k2n ). From authors point of view the matrix formalism based on two-dimensional matrices Rik (s|s0 ), which combine corresponding partial derivations, is more comfortable. 3 Formalization of a Global Optimization Problem In the previous section we consider some basic mathematical approach for beam, dynamics presentation. Before formalization of problems for optimization of beam line structure, one has to define the vector of con⃗ and functions U ⃗ (s). The control patrol parameters W rameters describe some characteristics of a beam line, which can not be changed during a experiment session. For example, lengths of drifts, lenses and other geometrical parameters can play the role of the control parameters. On the other hand characteristics of electromagnetical fields, generated by control elements, can play the role of control functions. 3.1 Control Function Formalization The motion equations (2) and the evolution equation of the matrix MN , envelope matrix S(s) and distribution function f0 (X, s) allow to follow the beam evolution with given external fields (in this paper we neglect the space charge). Usually in the theory and practice of beam physics electro-magnetic field induced by the control elements, such as dipoles, quadrupoles, sextupoles and solenoids. In this case, these fields can be represented in the form of Taylor’s series in the transverse coordinates, whose factors in general depend on the longitudinal coordinate along the reference orbit s (see, for example, [Chao, 1999]). For example, a scalar magnetic potential ψ can be represented in the following form: ψ(x, y, s) = ∞ ∑ k, i=0 where FL — the Lorentz force acting on the particle ensemble, E, B — electric and magnetic components of electromagnetic field satisfying Maxwell’s equations [Chao, 1999]. All these approaches to particle beam description are well agreement with so called matrix formalism (see, for example, [Andrianov, 2004]). The matrix formalism allows a researcher to use an unified tools for all problems of beam particle evolution. It should be note that there is a difference between this formalism and the widely used formalism based on tensor description of particle coordinates evolution xi (s) = ∞ ∑ |k|=1 ∑ k1 ,...,k2n k1 +...+k2n =|k| k2n Tik1 ...kn xk11 xk22 · · · x2n , aik (s) xi y k , i! k! where aik (s) satisfy algebraic-differential equations (see, i. e. [Andrianov, 2004]). For all control elements can be found aik to any order (using some reference coefficients aik (s)). For example, in the linear approximation for the quadruple magnet: Bx = gy, By = gx, Bs = 0, where g(s) = ∂By (0, 0, s)/∂x = ∂Bx (0, 0, s)/∂y. In these representations the functional factors facing the coordinates serve as the control functions. In the case of quadrupole lens is controlled by g(s). In the particles beam lines control elements follow each after another. All control functions for ⃗ (s) = elements can be grouped into one vector U ∗ (u1 (s), u2 (s), . . . , um (s)) , where ui (s), s = 1, m — describe appropriate control actions generated by control elements. For example, we can consider the linear approximation for beam dynamics in a system with dipole magnets, quadrupole lenses (and may be solenoids): x′′ + a1 x + a2 y = a3 + a4 s, y ′′ + b1 x + b2 y = b3 + b4 s, (4) where the functions ai (s), bj (s) describe fields in the corresponding control elements [Podzyvalov, 2008]. Using the piece-wise approximation of the control field these functions can be presented by a set of control parameters. Some of these parameters have a geometrical nature and other describe control field characteristics (see, for example, [Andrianov, 2004]). From the practical point of view the usage of control functions for control elements description it is not convenient. In this paper, following [Andrianov, 2004; Andrianov, 2009], we use a parametrization for function description in some appropriate function classes (see e. g. [Andrianov, 2008]). For example, for above mentioned case of piece-wise approximation for magnetic field gradient g(s) we receive a parameter gm , describing the value of field gradient and parameter Lm for the corresponding effective length. So in our problems we change vectors of control parameters and functions by one vector of parameters, in which one part have geometrical interpretation and other connected with control field characteristics. 3.2 Functional Requirements In general, the functional requirements of the particles beam can be represented in the form of the integral functional: J[A] = ∫sT ∫ s0 g1 (A, X, τ )dXdτ + M(s) + ∫ g2 (A, X, T )dX, MT where M(s) — a current phase set, occupied by beam particles. The function g1 describes the functional criteria distribution within the system (if necessary) and function g2 — the terminal beam requirements. The choice of the functions g1 , g2 is determined by the specific task and requirements to the system. Often, the first integral is represented as a finite sum, and the functional J[A] can be written in the form: J[A] = p ∑ i=1 αi Ji [A], where Ji [A] — partial functional responsibility for certain characteristics of the beam, αi — weight coefficients determining the contribution of a functional, i. e. its importance. Even in the case of short systems (with a small parameters number) value p can be quite large (see e. g. equation (4)), which leads to multi-parametric optimization problems. It should be note that Ji [A] often, to a certain extent, contradict each other, so we can speak about the contribution antagonism of the Ji [A] in J[A]. The meaning of these functionals and their contribution to the decision of choosing one or another solution will be described below in the global optimization terms. The approach is only support tool for the optimal solutions choice by researchers, witch have sufficient experience in the described class of problems and able to build a functional based on their experience. 4 Global Optimization Methods And Instruments Methods of global optimization can be divided into two categories: deterministic and heuristic. The first of these relates to the methods with a proof of optimality and allow to obtain a solution that differs from the optimum by no more than the specified amount. These methods allow to obtain solutions only for very small dimensions. In practice, often using heuristic algorithms [Weise, 2009] allowing to find a approximate minimum without proof of it’s optimality. This permits us to solve the problem of high dimensionality in a reasonable time with a fairly acceptable result. As already mentioned, even simple control systems of beam particles can contain a large number of parameters, therefore in the paper considered the heuristic search algorithms application. 4.1 A Classical Global Optimization Problem We can formulate a classical global optimization problem as a problem of detection the best possible elements x∗ from a set X in according to a set of criteria F = {f1 , f2 , . . . , fn }. These criteria are expressed as mathematical functions, the so-called objective functions. The exact mutual influences between these objective functions can apparently become complicated and are not always obvious. The simplest method accounting the generated criten ∑ ωi f i ria set is to define some weighted sum g(x) = i=1 of all functions f ∈ F. In some practical problems it is convenient to consider the following form of the n ∑ ωi fi2pi , where pi are optifunction g(x): g(x) = i=1 mization parameters too. Each objective function fi is multiplied with a weight wi representing its importance [Deb, 2001]. However, weighted sums are only suitable to optimize functions that at least share the same degree. Usually global optimization is founded on evolution- ary (genetic) algorithms realizing probabilistic searching for the parameter spaces by means of such genetic operator: crossover, mutation and inversion. The form of these evolution operators and their presence characterizes the properties of the optimal solutions research (see e. g. [Michalewicz, 1996; Weise, 2009]). Fig. 1.). Here, Li — the lenses length with their gradients k1 , k2 , and free field space a, g, λi are served as system control parameters. Focusing lenses are shown as elements of length L1 , defocusing as L2 . In this case, for simplicity of problem statement considered the similar system so called “russian quadruplet” [Dymnikov, 1997], often used like test system type. 4.2 Genetic Algorithms Genetic algorithms described in the works [Michalewicz, 1996; Weise, 2009] and many others are based on rules similar to most of them. These rules are described in the following sequence of steps: the initial population step – to create an initial population of random individuals; the evaluation step – to compute the objective values of the solution candidates; the fitness assignment step – to use the objective values to determine fitness values; the selection step – to select the fittest individuals for reproduction; the reproduction step – to create new individuals from the mating pool by crossover mutation and inversion; the condition breakpoint step – to continue the previous steps until a condition breakpoint. Each of these steps is a kind of subproblem which has an impact on the general problem solutions. Choosing the heuristic search methods is necessary to speak about “conditional global” and local extremums. This convention arises from a number of objective reasons. Indeed, in the process of modeling of such systems of accelerator physics, first consider the linear model, further complication and taking into account nonlinear effects occur only when the linear model satisfies certain requirements. This implies that the architecture of any system is only optimized with the specific task. And every time a researcher will attempt to build such a system and optimize it, he will receive only conditional optimal solutions. Sometimes, in the process of control system modeling developers are trapped in some design features. For example, the system must satisfy the features of the matching channel and therefore must primarily be considered with the number of controls and their layout (matching problem) [Ovsyannikov, 2008]. It follows that the process of modeling should be primarily organized in such a way that the developer can independently choose the major and minor functional and limitations. In the next section we consider the effectiveness of these methods and instruments of global optimization for a simple control system of beam particles. 5 An Ion-Optical System Example In the paper we demonstrate the suggested approach for a simple system consisting of two quadrupole doublets with non-zero distance between the lenses (see Figure 1. The control system of beam particles and envelopes in x-plane — the solid line, in y — the dashed line. The structure considered in Fig. 1 is designated as FODO–structure and widely used in various systems design in accelerator physics. This system can be optimized for its control parameters. For simplicity, in this paper the structure optimization is not considered. The usage of the matrix formalism allows us exclude some parameters from the optimization problem by introducing some analytical dependence. It is necessary to note that the similar reduction can be applied for nonlinear motion equation [Andrianov, 2004]. The others should impose some physical restrictions. 5.1 Symbolic Forms of Objective Functions Consider some of the character entry imposed on the optimization parameters in the modeling process. The transverse manifold described by S(s)-matrix, and occupied by beam particles as shown in the equation (3). Assume that M is a matriciant of focusing part of the system, then the total matriciant can be represented in following form: Rtot = Rg · Mobj · Ra , (5) where Mobj is the matrix for a objective (consisting from four quadrupoles, see Fig. 1). Taking into account the matrices Rg and Ra and using the necessary condition for beam focusing from point to point r12 = r34 = 0 (r12 corresponds to the coordinate x, r34 – to the coordinate y), we can write the equations for the parameters a and g: am11 + m12 , am21 + m22 am33 + m34 am11 + m12 = , am21 + m22 am43 + m44 g=− where mij – elements of the matrix M. It is necessary to note that elements r11 and r33 describe the increment (for r11 > 1, r33 > 1) or the decrement (for r11 < 1, r33 < 1) for beam size on the target [Andrianov, 2004]. Solving the equation (5) and substituting into the (3) we obtain the matrix S(s) on the target. Thus, the parameter g is excluded from the optimization process and due to the ideology of the matrix formalism [Andrianov, 2004], we recorded the focusing condition from point to point in an analytical form. Further, we can introduce some restrictions on the other parameters in order to simplify the optimization process. It is useful to note that similar simplification can be applied due to usage of the matrix formalism for beam dynamics description. 5.2 Conditional objective functions Using the above described linear approximation we can express almost all the physically significant beam characteristics √ in terms of the elements S(s)-matrix (for example, S11 is the maximum value √ of the coordinate x in envelope line, then element S33 is the maximum value of the coordinate y). Let us suppose that Sopt (s) is an optimal envelope matrix of the beam on the target. Then the deviation of “output” S(s)–matrix from optimal can be expressed in terms of some metric: ∥ S(s) − Sopt (s) ∥6 ε, where ε is determined by the nature of the specific physical problem. Assuming that the beam particles can move only in the working channel we can included the searching for parameters values that satisfy in the whole plane s. Resulting set of feasible solutions will require from the system designer to select the optimum solution (or a solutions subset) in terms of simplicity implementation, cost and other operational constraints. The corresponding analysis the set of optimal control parameters one can choose the values that correspond to the best tolerances (see, for example, [Andrianov, 2004; Podzyvalov, 2008]). This ensures that the adjustment will be stable enough for a small drift of the parameters in the manufacturing and assembling processes of the system and in its exploitation. 5.3 Graphical Interpretation Results An important aspect of any optimization process is a graphical representation of the optimization variables. However, in the case of a large number of parameters, the user of a special optimization tools must be able to choose what parameters should be reflected graphically. In the system described above, we are interested in a graphic study of the magnetic field gradient k1 , k2 of the lenses. Selection of optimal solutions of these parameters (or some subset) depends on the functional 2 2 , which minimizes the size of the beam + R33 F = R11 on the target. Beam characteristics on target incorporated in the matrix R(s) from equation (5). Shown in Fig. 2 the dependence of magnetic field gradients k1 , k2 suggests that there are wide and narrow areas, satisfying the functional F ≤ 1 (this means that the beam particles do not deteriorate during the channel). The magnetic field lens gradient values (for the functional F ≤ 1). Figure 2. To choose the optimal solution we use the data reflected in Fig. 3, 4, which displays the dependence of the each lenses gradient from the functional F . It is seen that most of the solutions for the gradient k1 , which have the minimum values of functional F are located at intervals [1.2,. . . ,1.8; 2.4, . . . ,2.8]. For the gradient k2 the minimum values of functional F are located at interval [0.6,. . . ,1.2]. Figure 3. F ≤ 1). The magnetic field gradient k1 (for the functional gain experience. The problem described the ion-optical systems optimization in term of the matrix formalism. It is shown practicability of the using the global optimization approaches in such considered problems. This ideology implies to use technologies of the parallel programming and can be solved by them. Figure 4. The magnetic field gradient k2 (for the functional F ≤ 1). Further analysis of these intervals can be made by specifying these intervals as the initial for the global optimization algorithm. This is possible thanks to a flexible approach in the program complex creating “Global Optimization Approach” (see Fig. 5.), which allows you to create different structures, configure a global optimization algorithm, graphically display the results and further investigate the extreme suspicious areas. Program complex screenshot “The global optimization approach”. Figure 5. The corresponding applied software include described global optimization concept ideology and implies to use distributed and parallel technologies for necessary computing. This allows the researcher to immerse themselves in the solution of the problem and make the system modeling process flexible and manageable. 6 Conclusion The approach and its implementation presented in this paper allow to an experienced researcher quickly find the optimal sets in parameters spaces with acceptable features of the beam, as a novice researcher can learn on test problems to investigate such problems and to References Andrianov, S. N. (2004). Dynamical Modeling of Beam Particle Control Systems, (St. Petersburg State University, St. Petersburg) (in Russian). Andrianov, S. N., Edamenko, N. S. and Podzyvalov, E. A. (2009). Some Problems of Global Optimization for Beam Lines, in J. of Mathematical Physics, Catania, Italy, http://lib.physcon.ru/download/p1998.pdf. Yang, L., Robin, D. (2008). Global Optimization of the Magnetic Lattice Using Genetic Algorithms, Proc. European Particle Accelerator Conference EPAC’08, Italy, Genova, pp. 3050–3052. Michalewicz, Z. (1996). Genetic Algorithms+Data Structures=Evolution Programs, 3-rd Edition. Springer. Gao, W. W., Li, W., Wang, L. (2010). Low Emittance Lattice Optimization Using Multiobjective Genetic Algorithm, in Proc. International Particle Accelerator Conference IPAC10, Kyoto, Japan, pp. 4515–4517. Chao, A., Tigner, M. (1999). Handbook of Accelerator Physics and Engineering, World Scientific. Podzyvalov, E. A. (2008). Modeling of adjustment errors in beam lines in XXXIX Int. Conf. on Control Processes and Stability — CPS’08, SPb, Russia, pp. 164–170 (in Russian). Tereshonkov, Yu., Andrianov, S. (2008). Load Curves Distortion Induced by Fringe Fields Effects in the Ion Nanoprobe, in Proc. European Particle Accelerator Conference EPAC’08, Italy, Genova, pp. 1514–1516. Weise, T. (2009). Global Optimization Algorithms: Theory and Application, 2-nd Edition. Deb, K. (2001). Multi-Objective Optimization using Evolutionary Algorithms, Chichester. Ovsyannikov, A. D, Ovsyannikov, D. A., Sheng-Luen (2008). Optimization of radial matching section, Proc. of BDO08, St.Petersburg, Florida, pp. 952–958. Dymnikov, A. D., Marti’nez, G. (1997). Optimal magnetic and electrostatic Russian Quadruplet microprobe lens system with high demagnification, in J. Nuclear Instruments and Methods in Physics Research, Volume 130, Issues 1-4, pp. 64–69.