Quadrupole Magnet Design based on Genetic Multi-Objective Optimization

This work suggests to optimize the geometry of a quadrupole magnet by means of a genetic algorithm adapted to solve multi-objective optimization problems. To that end, a non-domination sorting genetic algorithm known as NSGA-III is used. The optimization objectives are chosen such that a high magnetic field quality in the aperture of the magnet is guaranteed, while simultaneously the magnet design remains cost-efficient. The field quality is computed using a magnetostatic finite element model of the quadrupole, the results of which are post-processed and integrated into the optimization algo-rithm. An extensive analysis of the optimization results is performed, including Pareto front movements and identification of best designs.


Introduction
Multipole magnets are essential components of particle accelerators and of paramount importance for the success of particle physics experiments that are performed in accelerator facilities [44].Quadrupole magnets, in particular, are crucial for keeping particle beams focused on their desired trajectories [9].The very nature of particle physics experiments places very high demands on the operation and performance of accelerator magnets.Therefore, their production must be based on meticulous design procedures.
Focusing on the case of quadrupole magnet design, common practice dictates to formulate the design objectives into an optimization problem with respect to the geometric parameters of the quadrupole, while certain operational parameters must also be taken into account during the optimization process.Therein, the electromagnetic phenomena taking place inside the magnet are typically simulated using a digital magnet model and are included into the optimization as performance measures [1,28,31,43,47].However, such geometry optimization problems are notoriously hard to solve, in fact, NP-hard [2], due to the fact that their computational complexity increases exponentially with the number of solution candidates.In many cases, the only viable option is to employ search algorithms that are capable of comprehensively exploring the parameter space and finding an adequate configuration at a reasonable computational cost.Popular methods of choice are Monte Carlo algorithms [4,37,38], simulated annealing [16,40] or genetic algorithms (GAs) [26,32,45,50], all of which have been successfully applied to a number of geometry optimization problems [3,46,49].
This work focuses on the use of GAs for the purpose of quadrupole design.Originally, GAs have been developed for unconstrained, single-objective optimization problems, which constitute their natural domain of application [36].However, most geometry optimization applications, also including quadrupole magnet design, feature multiple and often conflicting objectives.In such cases, a multi-objective optimization (MOO) [23,30] must be solved instead.Therein, the goal is to find so-called Pareto optimal solutions, which cannot be further improved with respect to one objective without worsening another [35].Such MOO problems are faced and tackled increasingly more often in various engineering applications that concern geometry optimization [5,15,19,51].To address these problems, traditional GAs have been extended, e.g., in the form of so-called non-dominated sorting genetic algorithms (NSGAs) [11,13,29,48].
Regarding the quadrupole magnet design problem considered in this work, the objectives of the MOO are: (a) the maximization of magnetic field quality within the aperture of the quadrupole magnet, and (b) the minimization of the magnet's radius.The former objective is aligned with the desired magnet operation with respect to beam focusing.The latter objective aims at designs that use the minimum necessary amount of magnet material, so that production costs remain acceptable.The magnetic field quality is first computed using a magnetostatic finite element (FE) magnet model and then expressed in the form of relative multipole coefficients on a reference radius.The evaluation of the relative multipole coefficients is integrated into the fitness function of an NSGA, such that only solutions which maximize the magnetic field quality are allowed to propagate through the evolution process, as dictated by the first optimization objective.In particular, the algorithm known as NSGA-III [11] is employed.The solutions are further constrained by the second objective, i.e., the need for a minimal magnet radius.The Pareto-optimal solutions, i.e., those that provide an acceptable balance with respect to both objectives, are then analyzed such that the best magnet design candidates are identified.
While works featuring GA-based MOO algorithms for the optimization of accelerator systems do appear in the literature [18,25,27,33,42,52], the application of such optimization algorithms to accelerator magnet design seems to have been so far neglected.In fact, the only relevant work that the authors are aware of is [44], which mentions the option of optimizing accelerator magnets by means of MOO based on GAs, however, without providing any verification in the form of numerical experiments.The present work aims to fill this gap.
The rest of this paper is organized as follows.Section 2 presents the problem formulation with respect to computing the magnetic field of the magnet for a given design, the magnet model considered in this work, and the computation of the magnetic field quality in the aperture of the magnet.Section 3 presents the general formulation of MOO problems and discusses GAs suitable for the solution of such problems, in particular the NSGA-III.In Section 4, the MOO problem with respect to the quadrupole magnet design is formulated.The numerical results obtained by means of the NSGA-III algorithm are presented and extensively discussed in Section 5. Finally, conclusions are drawn in Section 6.

Problem Formulation and Magnet Model
The physical behavior of a quadrupole magnet can be fully described by the Maxwell equations, which provide the mathematical foundation for classical electrodynamics.However, in most cases, it suffices to provide an approximate description of the electromagnetic phenomena appearing in a given problem setting, thus simplifying the underlying equations and the resulting simulation model.For our particular application of a quadrupole magnet with nonlinear materials, it is sufficient to consider a magnetostatic representation of the underlying magnetic field quantities.Therefore, we employ the magnetic vector potential formulation in the steady state regime.
Considering the 3D case, the magnetic vector potential formulation is given as where ν is the magnetic reluctivity tensor, b the magnetic flux density, a the magnetic vector potential, and j the current density.Equation ( 1) is discretized by employing a projection on (1) with a test function v on a computational domain Ω with boundary ∂Ω = ∂Ω N ∪ ∂Ω D , where ∂Ω N and ∂Ω D respectively denote the boundary parts where Neumann and Dirichlet boundary conditions (BCs) are imposed.The corresponding weak formulation reads: [7] find a ∈ L 0 (curl, Ω) such that for all test functions v in the space L 0 (curl, Ω), defined as In (3), H(curl, Ω) is the space of square-integrable functions with squareintegrable weak curl.Note that the Neumann boundary integral arising in the deduction of (2) vanishes by applying the Neumann BC n • b = 0, where n is the outer normal unit vector.The Dirichlet BC n × a = 0 is enforced in L 0 (curl, Ω).
The magnetic vector potential is approximated within a finite element space as where âj are the degrees of freedom (dofs), E is the number of dofs, and w j denotes Nédélec basis functions of the first kind and the first order [41].We apply the Ritz-Galerkin procedure, such that the set of test functions is the same as the set of shape (basis) functions.

2D Magnet Model
The magnet model is further simplified taking into account the translation invariance of the magnet along the z-axis.Therefore, we may consider a magnetic vetor potential perpendicular to a 2D cross section of the 3D domain, in which case a = [0, 0, a z = a z (x, y)] and j = [0, 0, j z ] .We denote the computational domain of the 2D cross section of the quadrupole magnet with Ω 2D := {x ∈ R 2 | x 2 ≤ R}, which corresponds to a circle with radius R ∈ R around the origin.Figure 1 depicts the full 2D quadrupole magnet model contained in Ω 2D , where the corresponding geometry descriptors and numerical identifiers are given in Table 1.The Dirichlet BC n × a = 0 is imposed on the boundary ∂Ω D , also shown in Figure 1.
Table 1: Geometrical parameters of the quadrupole magnet.

Description Identifier Notation Value Units
Pole width Fig. 1: 2D cross section of the quadrupole magnet.The computational domain, its boundaries, and the different material-based subdomains are presented as in the legend.The identifiers 1 − 6 and the letters A, B, denote the geometrical parameters of the magnet, as described in Table 1.
Table 1 also presents the intervals [a i , b i ], within which the geometrical parameters are allowed to vary during the optimization procedure, i.e., a i ≤ x i ≤ b i , i = 1, . . ., 6. Parameters with constant values throughout the optimization are denoted with the identifiers A and B. Shims are included in the quadrupole model as important pole adjustments, which can lead to significant improvements in the homogeneity of the magnetic field, thus, to field quality improvements as well [24].The non-linear material of the yoke and the poles is modeled with a Brauer curve approximation [8] upon 1010-Steel and implemented in a FE solver with the Newton method [17].The model is implemented using three open-source tools, namely, the mesh generator Gmsh [20], the GetDP FE solver [17], and the ONELAB interface [21].

Aperture Field Quality
One of the most important quantities of interest to be taken into account during the design of a quadrupole magnet, is the field quality represented by the harmonic distortion factor Q, which can be computed upon the multipole coefficients of the calculated field in the magnet's aperture.To that end, the FE solution a z is evaluated at a reference circle with radius r ref .The result is then represented by a Fourier series with Fourier coefficients a p and b p , such that using the polar coordinate system (r, ϕ).Then, the Fourier characterization of the magnetic vector potential and the magnetic flux density in the beam aperture are given as [44] The evaluation at r = r ref for the radial magnetic flux density yields where B p and A p are called normal and skew harmonic coefficients, respectively.The harmonic distortion factor Q(r ref ) in the aperture of a 2P -pole magnet considered at a reference radius r ref can be obtained from the harmonic Fourier coefficients using the formula [44] Q For a quadrupole magnet (P = 2), with a quadrupole field gradient the magnetic field should be close to a pure quadrupole field, i.e., the field gradient g should be high in relation to the other multipole coefficients B p , with p = 1, . . ., C, and p = 2, where C refers to the truncation coefficient.Accordingly, the corresponding harmonic distortion factor for a quadrupole magnet, given by should be in the order of 10 −4 [44].

Multi-Objective Optimization
In MOO, we seek to minimize (or maximize, depending on the problem at hand) a set of -possibly conflicting -objective functions f i (x) : X → R, i = 1, ..., k, k ≥ 2, where x = (x 1 , . . ., x n ) ∈ X denotes the vector of decision variables, equivalently, optimization parameters, and the domain X ⊂ R n is referred to as the feasible decision region [23,30].For a feasible parameter vector x ∈ X , a corresponding feasible objective vector z = (f 1 (x), . . ., f k (x)) is obtained, where z ∈ Z and Z ⊂ R k is called the feasible objective region.
The structure of X is induced by a set of constraints applied to the decision variables.In the specific case of geometry optimization, these constraints are often given in the form of bounding box intervals, similar to the ones shown in Table 1, such that Given this set of constraints, the MOO problem reads In most cases, a parameter vector that minimizes all objective functions simultaneously does not exist.It is therefore necessary to have a method of comparing a set of solutions while taking into account the satisfaction of all objectives.This issue is resolved using the concept of dominating solutions.Assuming two parameter vectors x, x ∈ X arising in a minimization procedure, then x dominates x if A parameter vector that cannot be dominated is called Pareto-optimal [35].In essence, Pareto-optimality means that the current solution cannot be further improved with respect to one of the objectives, without simultaneously deteriorating another objective.The set of Pareto-optimal solutions is referred to as the Pareto front.A Pareto front is bounded by the ideal and nadir objective vectors, respectively denoted with z * and z nad .The former is obtained by individually optimizing the objective functions and the latter by approximating the worst objective values of the Pareto front.

Genetic Algorithms
GAs belong to a class of population based, stochastic optimization algorithms which solve optimization problems by only allowing candidate solutions with a promising "gene pool" to reproduce and propagate through generations t = 1, 2, . . ., T , where T is the final generation [45,50].This generation-based evolution of candidates is realized by creating a sequence of subsets P in the feasible decision space called populations, which eventually converge to a set of minimizers.In that way, a sequence (P t ) t∈N ⊂ X with is generated, such that each x ∈ P * is a minimizer of (13).To deal with the limitations of realistic and thus finite calculations, Equation ( 15) needs to be truncated by a final generation T , so that with a set of generation-related minimizers x ∈ P T is considered.The sequence of populations (P t ) t=1,2,...,T is dictated by the genetic operators crossover C σ , mutation M r , and fitness selection F , such that where • expresses the concatenation of the sequentially applied genetic operators on the current population P t .The indices σ and m respectively denote the crossover and mutation rates of the corresponding operators, where σ, m ∈ [0, 1].The crossover and mutation rates determine the probability that the given operator is applied to a given population.The crossover operator C σ describes how sample solutions are recombined to generate new solutions for the next population.The mutation operator M m describes random distortions to the elements of a population and is particularly significant for the convergence of the GA, as it ensures that the objective space X is searched comprehensively and that the limit P * is initialization-independent.Last, the fitness selection operator F allocates a fitness value to the population members and selects those with the highest values to progress to the next generation.In this work, simulated binary crossover [12,14] and polynomial mutation [22,53] are employed as crossover and mutation operators, respectively.For fitness evaluation, we use the FE model of the quadrupole to compute the magnetic field distribution in the magnet and evaluate the aperture field quality, which in turn determines the fitness of a given population.As a selection operator, we use binary tournament selection [39], due to its ease of implementation and robustness against stochastic noise.

Non-Dominated Sorting Genetic Algorithms & NSGA-III
As noted before, GAs were originally developed to solve single-objective optimization problems.Therefore, they cannot address a number of issues related to MOO, such as dealing with multiple objective functions and ensuring diversity in the populations.These issues have been addressed with the introduction of NSGAs [11,13,29,48].In this work, we resort to the so-called NSGA-III algorithm [11,29], which is briefly discussed in the following.First, NSGA-III deals with the issue of optimising numerous objective functions by preferentially handling solutions that dominate other members of a population according to definition (14), by using non-dominated sorting.
Given the current population P, non-domination sorting partitions the population This hierarchy is induced according to a domination factor n p = 0, 1, ..., N −1, which indicates how many solutions from an equal or lower hierarchy level dominate the given solution.The hierarchy construction is as follows: The set F 1 includes x ∈ P, which have n p = 0, i.e, they are not dominated by any solution.Then, we consider the set Q = P\F 2 and decrement the domination factor for n p → n p − 1 for all q ∈ Q that are dominated by an element of F 1 .This process is repeated iteratively, thus yielding a hierarchical set sequence.The sets with lower domination factors qualify to the next population, whereas the sets with higher domination factors are discarded.It is thus ensured that that the solutions propagating to future generations are Pareto-optimal with respect to the current population they belong to.
To ensure population diversity, NSGA-III adds another operation to the fitness selection procedure.Therein, the objective vectors z are normalized to the unit cube by using the ideal objective vector z * and the nadir objective vector z nad. .In that way, it is possible to consider objective functions that are scaled differently.Then, reference points on the unit hypercube are chosen, which lie on a simplex [10].The reference points typically have a space-filling property and the objective vectors are projected to the reference points.The population members are then determined by an explicit diversity-preserving mechanism.

Quadrupole Magnet Optimization
The MOO concerns maximizing the absolute value of the field gradient g, as introduced in Section 2.2, while minimizing the outer radius R of the magnet.The optimization parameters are the six geometrical parameters listed in Table 1 and the current density j.The latter takes values within the interval [1.0, 20.0] A/mm 2 and in the following is denoted with x 7 .Then, the MOO problem reads where the feasible decision space is Instead of optimizing the full 2D magnet model shown in Figure 1, we exploit the three mirror symmetries of the magnet model in order to reduce it to the one-eighth segment depicted in Figure 2.This model reduction leads to a significant improvement in terms of the computational cost of the finite element method (FEM).
Besides constraints on the geometrical parameters and the current density, we also introduce constraints on the absolute duodecapole gradient g d , as well as the on the saturation behavior of the iron yoke of the magnet, which read Analogous to the quadrupole field gradient g in Equation ( 11), the duodecapole gradient is given by g d = B 6 .In any case,the additional constraints demand that the duodecapole gradient g d should not exceed a certain percentage of the absolute field gradient g, as well as that the magnetic flux density entering the pole nose bin and the magnetic flux density exiting the yoke bend bout are bounded by a fixed saturation threshold value b th .The magnetic flux densities bin and bout are illustrated in Figure 2.Both are obtained by post-processing the FE solution a z .
Fig. 2: One-eighth of the 2D quadrupole magnet model.
The solutions to the MOO problem defined by equations ( 18)-( 19) are obtained, as previously noted, using the NSGA-III algorithm.In particular, the implementation of the algorithm which is available in the open-source, Python-based optimization software pymoo is used [6].

Numerical Results
In this section, the numerical results of the MOO are presented, which are additionally employed to identify the best solutions, i.e., the ones corresponding to the most suitable magnet designs.The notion of the best solution is split into: • Best balanced solution, that is, the Pareto-optimal solution with the lowest Euclidean distance to the ideal objective vector, where the latter is approximated with respect to the final Pareto front.• Best field gradient solution, that is, the Pareto-optimal solution found in the final Pareto front, which results in the highest field gradient.
Using NSGA-III with a crossover rate σ = 1.0 and a polynomial mutation rate m = 0.1, each optimization run is performed for T = 300 generations, with an initial population of 80 individuals, and an offspring population of 56 individuals.These empirically enforced values remain fixed during several optimization runs, where the influence of different saturation behaviors of the iron are investigated.To analyze the saturation behavior, three optimization runs are performed, within which the saturation threshold b th , which is embedded in the optimization's constraints as in Equations (19), is selected as 1.0 T, 1.2 T, or 1.4 T, respectively.Each optimization run is performed upon the same feasible decision space X .The initial point is chosen to be the center point of X and represents a naive but admissible choice of geometry with respect to the optimization constraints, see Figure 3 Figure 3 shows the Pareto front movement and the location of the final Pareto front after t = 300 generations for the three chosen saturation thresholds.Therein, the Pareto front movement is obtained by connecting the mass points of the fronts obtained from each generation of the NSGA-III, where the mass point is defined as the mean of all Pareto optimal solutions.As can be seen, both the final Pareto front and the front movement depend strongly on the chosen saturation threshold.For b th = 1.0 T, the final Pareto front is biased towards low values of the outer radius r.In contrast, the Pareto front for b th = 1.4 T is biased towards high field gradient values.A better balance between the objectives is achieved for b th = 1.2 T. As for the individual evolution of the different objective spaces, Figure A1 in the Appendix shows each of the optimization solutions and their Pareto front movements individually.Last, it is verified that the Pareto optimal solutions of all optimization runs have a harmonic distortion factor Q ≤ 10 −4 .This observation is consistent with the theoretical requirements for the harmonic distortion factor of a sufficiently undistorted quadrupole field, as mentioned in Section 2.2.
Next, focus is shifted to the best solutions of the final Pareto front for each saturation threshold.For a better visual comparison of the best solutions, the parameter combination corresponding to each solution is normalized to the unit cube X norm = [0, 1] 7 .Using this normalization, the parameter distributions and the associated objective values of the best solutions are depicted in Figure 4.As can be observed, the saturation threshold has a strong influence on the magnet's pole width, both for best balanced and for best field gradient solutions.This is attributed to the fact that higher saturation thresholds lead to a larger pole width design.It is further observed that larger pole widths lead to higher field gradients for both cases of best solutions.The saturation threshold has also a slight to moderate influence on the pole bending design.With regard to the current flow, it can be seen that, for both best solution cases, higher saturation thresholds allow higher current densities, as expected.In the case of the best balanced solutions (Figure 4a), it is observed that changes in the saturation threshold have almost no impact on the dimensions of the pole height, yoke height, and the shim geometry.Contrarily, the dimension of the pole height is more important in the case of the best field gradient solution, as can be seen in Figure 4b.Additionally, the choice of the shim geometry seems to be more important for best field gradient solutions than for best balanced solutions.

Conclusion
This work presented a framework for optimizing predominantly the geometrical and secondarily the operational characteristics of a quadrupole magnet.The framework employs a MOO formulation, where two conflicting objectives must be satisfied, namely, high magnetic field quality and acceptable production cost.The MOO problem is solved by means of the so-called NSGA-III algorithm, which is a GA suitably modified to address the issues arising in MOO.Therein, a magnetostatic FE model of the magnet is employed in order to assess the quality of the magnetic field in the aperture of the magnet.Finally, the MOO problem is complemented with additional constraints on the duodecapole gradient and the saturation threshold of iron.
The numerical results indicate that saturation has a major impact on the obtained Pareto-optimal solutions.Further analyzing the connection between optimization parameters and Pareto fronts, it is possible to deduce useful information regarding the impact of the different geometrical characteristics of the magnet onto the optimization objectives.Importantly, all identified optimal magnet designs lead to a sufficiently low harmonic distortion factor, while it is possible to identify designs with an acceptable balance between field quality and production cost.
Future work should consider the utilization of the 3D magnet model, in order to consider fringe field effects during MOO studies.The impact of eddy currents and power losses should also be taken into account.As a result, MOO problems with k > 2 objectives can be formulated and investigated, leading to even more improved magnet designs.To limit the computational burden of MOO, especially for an increased number of objectives, surrogate modeling approaches could be considered [34].

Declarations
Conflict of interest: The authors declared no potential conflict of inter-ests with respect to the research, authorship and/or publication of this article.
. The associated initial objective values are R = 0.115 m and |g| = 11.305T/m.Here, the initial point serves as a reference point to compare the locations of the different final Pareto fronts.

Fig. 3 :
Fig. 3: MOO results for three optimization runs, each corresponding to a different saturation threshold b th .The arrows show the Pareto front movement from the initial point (×) towards the mass point of the final Pareto front ( ).

Fig. 4 :
Fig. 4: Best solutions (normalized) for different saturation threshold values b th .Additional to the best solutions, the corresponding objective values R and |g| are given.Filled black circles: normalized parameter combinations of the best solutions.Transparent blue-colored circles: normalized parameter combinations of all final Pareto-optimal solutions with respect to the different saturation thresholds.