Design Optimization for Superconducting Bending Magnets using Pareto Front Curve

A novel limit design method for superconducting magnets is presented. It is particularly suitable for ion core magnets such as those used in accelerator magnets. In general, a stochastic optimization whose objective functions consist of values, e.g., the magnetic field, experience field of superconducting coils, current density, and multipole field integral, is often used. However, it is well known that the obtained solution strongly depends on the initial one. Furthermore, once the calculation model is fixed, the range of solutions is also fixed, i.e., there are times when it may be impossible to find the global optimum solution even with a lot of parameter sweeps. In this study, we draw the Pareto front curve to obtain the range and infer whether the solution is an optimum one. In addition, the Pareto front curve indicates the neighborhood solution that is substituted for the initial one. After this process a stochastic optimization is implemented with its initial design parameters. To confirm the validity, we designed a superconducting bending magnet, and it showed that this method works well.


Introduction
Stochastic optimization methods have often been used in engineering design [1][2][3]. These methods have many merits such as easy handling. The optimization algorithm minimizes the single objective function or multiple objective functions. The single objective optimization requires breaking down several constraints into a single objective function while considering the differences in sensitivity of each function against the parameters. As a countermeasure, we have already proposed a method to generate an objective function with automated weighting [4]. In cases where a lot of constraints exist, the weighting is even more critical. In this technique, the ratio of the evaluated function value over each penalty function value is used to determine each weight. Each function value is smoothed through spline fitting. These weights should be updated automatically at predefined times. If they are constant the relationship among their function values changes with times. This automated weighting technique can lower the values of each function impartially, thereby avoiding obtaining local optimum solutions. The multiple objective optimization can estimate each constraint partially and obtain many non-inferior solutions. However, this needs a huge amount of calculation because there are a lot of design cases.
Furthermore, the optimization of design parameters is often conducted with the calculation model fixed. Once the model is fixed, the obtained solution could be limited to the area where the global optimum solutions does not exist, i.e., the global optimum solution may not be included in the solution set even with the use of an advanced variable optimization.
For a practical approach to dealing with optimization problems, the concept of Pareto optimization [5] The Pareto front solutions are not dominated by any other  feasible solution, i.e., no other solution exists that is better for at least one objective function value and equal or superior with respect to the other objective function values. This fact has been widely used in engineering to aid designers in their decision-making processes. In this study, we proposed a method: at the first step determine the possibility of satisfying the demanded specifications from the Pareto front curve and then implement the optimization with the neighborhood solution as the initial solution.
To draw the Pareto front curve we used one of the multi-objective optimizations, PSO (particle swarm optimization). PSO enables effective plotting on the solution space without wasteful computation. A feature of this flow leads to avoiding obtaining the local optimum solutions. The proposed method was applied to a bending magnet design, which showed good performance in solving this problem.

Drawing Pareto front curve and determining computational model
A solution space typically exists under inequality constraints, and its boundary is composed of a set of Pareto optimum solutions. In this study, PSO was used to generate the solution space. This design set of PSO parameters is shown in Table 1. The PSO concept was first introduced by Kennedy and Everhart in 1995 [6], and they showed that problems are optimized by having a population of candidate solutions, using dubbed particles, and moving these particles around in the search space according to simple mathematical formulae that describe the particle's position and velocity. As shown in Fig. 1, several sets of initial parameters generate the Pareto front curve through successive alternation of generations. This curve indicates the boundary of the solution space, i.e., we can judge whether the current model is a solution candidate through comparative evaluation with the required specifications. The curve also gives the neighborhood solution. With this solution as the initial solution, we can optimize the parameters as described below.

Number in population 3000
Number in generation 25

Number in parameters 18
Maximum archive number 5000 Sharing Crowding distance

Single objective stochastic optimization
In the previous section, the computational model was determined using the Pareto front curve. To optimize the rest of the parameters, an adaptive simulated annealing (ASA) method [7] was chosen as the single objective optimization algorithm. ASA departs from traditional simulated annealing (SA) by considering a generating probability density function with fatter tails than the typical Boltzmann distribution, allowing ASA to possibly escape local minimum solutions by considering farther points. The objective function should be created by summing each constraint function with weights. ASA has some characteristic features, for example, the parameter sensitivity, that is, the gradient of the objective function, can be estimated, separate annealing temperature parameters are assigned for each variable, and they are updated as the optimization progresses. The entire computational flow in this study is shown in Fig. 2.

Modeling
A bending magnet design was implemented as an example of the proposed optimization. A schematic drawing of the model is shown in Fig. 3. This magnet is mainly composed of an iron yoke and a pair of superconducting race-track-type coils. The cross-sectional area of the superconducting coil was assumed to fit in a square of side 55 mm based on the result of the investigation of an optimized wiggler design for SAGA-LS [8]. The constraint list for this optimization is shown in Table 2. For an engineering design, a value close to the constraint, that is, a limit design is favorable. The magnetic field was calculated by using a commercial FEM solver (EM Solution, Science Solutions International Laboratory, Inc., Japan). The calculation points were set on the beam axis (z-axis) at 5 mm intervals for the dipole field integral, on the plane of y = 0 at 1 mm intervals along the z-axis and at 0.5mm intervals parallel to the x-axis for the sextupole field integral calculation. The field distribution around the beam axis is expressed in Eq. (1) [9].
To calculate the multipole field integral described later, we found the approximate curve to the distribution parallel to the x-axis as shown in Eq. (2).   Figure 3. Schematic drawing of bending magnet with iron yoke. Race-track-type coil is arranged inside square of side 55 mm.

Drawing Pareto Front Curve
In the beginning, the maximum race-track straight length was assumed to be 126 mm. In consideration of an arbitrary coil cross-section, four coils were arranged in a single layer as shown on the left of Fig. 4. The current density was defined to be the same for each coil. The total number of design parameters is 18 (for 4 sets of coil center positions, each size of cross-section, the race-track straight length, and the current density). Figure 5 shows the results of plotting solutions. The horizontal axis is the sextupole field integral, and the vertical axis is the coil experience field. The circles mean with the race-track straight length of 126 mm and the squares mean with the other lengths. The solutions were plotted only when the value depending on the equality constraint was within a 2% error of the target value. This result shows that the Pareto front curve is made up of only the solutions of the 126 mm case, and a proper length exists in the condition where the sextupole field integral limit is 30 T/m. The sextupole field integral represents the spatial field homogeneity in the x direction, and the curved parts of the race-track affect the homogeneity, i.e., the sextupole field integral can be varied by changing the straight length of the race-track. The current distance from the beam axis (z-axis) to the curved parts is 126 mm/ 2 = 63 mm. As the magnetic field is spatially reduced at an exponential rate, the magnetic field gradient is represented as

Single objective stochastic optimization
Due to the above estimation, the race-track straight length was fixed. In the next step, the rest of the parameters regarding the coil cross-section shape and its current were optimized by using ASA. The single objective function ) (x OF was defined by the following equations. represent the peak field, dipole field integral, sextupole field integral, coil experience field, and current density, respectively. The vector x is a set of the design variables. The dipole field integral, ) (x BL , is calculated with Eq. (9). This value is the integral calculated with magnetic field ) ( 1 z b (T) along the beam axis (z-axis).
The sextupole field integral, , is an indicator of the magnetic field homogeneity. The sextupole field integral is represented with the coefficient At the bottom of Fig. 5, the coil cross-sectional diagrams that are neighborhood solutions are shown. The tendency of the four coil arrangements is that an empty region exists at the upper right corner, and the coil experience field becomes higher than 7 T, which means the current density is high. The spread of the cross-section area is an easy countermeasure to a decrease in the current density. In consideration of the above, the coil cross-section shape that was composed of two adjacent coils was defined as shown on the right of Fig. 4.
In the final step, an optimization with ASA was conducted. The transition of a single objective function against iteration number is shown in Fig. 6. After about 200 iterations, the feasible solution shown in Table 3 was obtained. The sextupole field integral is close to the target value of 30 T/m, which proves that the assumption of a race-track straight length of 106 mm was reasonable.   For validation, solutions with a race-track straight length of 66, 76, 86, 96, 106, 116, and 126 mm were plotted with Pareto front curves as shown in Fig. 7. The gradient of each curve is similar and is simply shifted along the horizontal axis. From Fig. 7, it can be seen that the sextupole field integral with the 106 mm model is about 28 T/m under an experience field of 6.5 T as expected, and the obtained solution is almost on the Pareto front curve. The modified Pareto front curve indicates which model could be an optimum candidate. Incidentally, regarding the experience field, the solutions over 10 T are apparently infeasible. However, these solutions are useful for drawing a Pareto front curve, because they are determined by factors related to not superconductivity but magnetomotive force. The cross-section of this obtained coil was L-shaped. To confirm that this shape contributes to diminishing the sextupole field integral rather than the simple square shape, we drew the magnetic field distributions with both shapes. As shown in Fig. 8, the distribution with the L-shape is flatter than that with the square one, i.e., homogeneity was definitely improved.

Conclusion
A new effective method for an optimal magnet design was proposed. A conventional optimization method has a higher possibility of falling into a local solution because of its high dependency on the initial parameters. The proposed method used a Pareto front curve for design parameter setting. The proposed method indicates the range of obtainable specifications with the current model. Based on this result, we could sophisticate the computational model and its initial parameters. It was shown that the optimization by the proposed method works well, and the solution reaches around the global solution as expected.