Robust and Multi-Objective Pareto Design of a Solenoid

: The optimization of the design of a practical electromagnetic device involves many challenging tasks for new algorithms, especially those involving numerical modeling codes in which objective function calls must be minimized for practical design processes. The Compumag Society provides openly accessible, challenging benchmark problems (TEAM problems) for testing novel numerical solvers. This paper deals with a novel solution for the multi-objective TEAM benchmark problem. This solenoid design test problem aims to search for the optimal shape of a coil, which en-sures a uniform ﬁeld distribution in the control region, while the sensitivity and the mass/DC loss of the coil are also considered in the context of robust design. The main differences from the previously published solutions are that the proposed methodology optimizes all three objectives together, not only as two independent two-dimensional sub-problems. We considered the asymmetrical cases in the solution and found that the symmetrical solutions always produced better uniformity and sensitivity measures. However, the difference between the symmetrical and asymmetrical solutions is insigniﬁcant for these objectives. Despite the fact that the cheapest solutions are symmetrical setups, they perform worse than the cheapest asymmetric ones in these uniformity and sensitivity criteria. Therefore, some asymmetric solutions that were previously neglected from the solution space can be competitive and interesting for practical design.


Introduction
The practical design of an electrical device usually leads to a multi-objective optimization task. These problems must involve the resolution of many computationally expensive finite-element-methodology-based numerical field calculations. The goal of the approaches that are applied is to reduce the number of function evaluations or to reduce a numerical model's computational cost without a significant loss of accuracy [1,2]. During industrial design processes, not only the physical parameters of a machine, but also the manufacturing tolerances and the different uncertainties should be considered right from the beginning of the design process [3][4][5][6].
The goal of a design optimization task differs when it is assessed from a mathematical or industrial perspective. Mathematically, the goal of a design optimization task is to find not only a better solution, but also the global optimum of the task, if it exists [7]. However, models that can be used for optimization are usually simplified ones that do not consider many factors, which is important for the easy and robust manufacturability of the product. Moreover, the design of an electrical device usually leads to a general nonlinear optimization problem. The result of these optimization problems is a Pareto-front, i.e., a set of non-dominated solutions [2,8,9]. From the industrial point of view, solutions that are better than the previous ones are usually considered as optimal ones because many factors are neglected during industrial optimization processes.
The TEAM (Testing Electromagnetic Analysis Methods) (https://www.compumag. org/wp/team/, (accessed on 1 June 2021)) benchmark problems offer a wide variety of test problems for benchmarking the accuracy of numerical solvers, and they are openly available from the website of the International Compumag Society [10]. The subject of our analysis is a multi-objective Pareto optimization of a solenoid, which is cataloged as the 35th benchmark problem in the list [11,12]. The goal of this benchmark is to design a coil that generates a uniform magnetic field in a given control region ( Figure 1). The sensitivity to positioning errors of the turns and the power loss or mass of the given design are also considered. This seemingly simple test problem is inspired by a bioelectromagnetic application for Magnetic Fluid Hyperthermia (MFH), where a uniform magnetic field is used to compare the magnetic properties of different nanofluids [5,[13][14][15]. A similar design task should be solved during the design and optimization of induction heating or induction brazing processes [16][17][18][19][20][21][22]. The solution of this test problem requires the resolution of a three-objective optimization problem in which the objective functions depend on finite-element-method-based calculations, and one of the objective functions considers the robustness of the solution during the optimization process.
Many solutions were proposed to resolve this problem in the literature [23]. The first paper proposed the DC problem and computed this optimization task through a gradientbased evolutionary algorithm (EA) search [12]. The following paper resolved the same problem with different EAs [24]. This comparison has great importance because several EAs were used to determine inverse electromagnetic tasks. Due to Wolpert's "No-freelunch" theorem [2,25,26], these metaheuristics must be benchmarked with each other for every kind of optimization task in order to select the most appropriate one. Seo et al. [27] solved this problem with a design sensitivity analysis, which provided almost the same optimization result as the gradient-search-based evolutionary algorithms in a much shorter time. These authors used FEM solvers to calculate the magnetic field. Karban et al. [5] proposed a semi-analytical formula and validated its precision with an hp-adaptive FEM solver [28]. The goal of this semi-analytical formula was to accelerate the solution speed of the magnetostatic problem. This problem was solved by other authors with different evolutionary and genetic algorithms, such as Non-Dominated Sorted Genetic Algorithm (NSGA-II) [29], Wind-Driven Optimization [30], the Micro-Biogeography Inspired Multi-Objective Optimization (µ-BiMo) [31], and the Migration Non-Dominated Sorted Genetic Algorithm (MNSGA-III) [32]. Another paper modified the excitation and solved this problem with the time-harmonic regime as a 2D or a 3D problem, considering the proximity effect of the windings [23]. However, these previously published solutions did not consider the measurement limits and other confounding factors, such as the Earth's magnetic field, which can be greater than the contributions of these effects or the difference between two designs. This paper proposes an approach for the original DC problem that is different from those of the papers that were mentioned earlier [12,24,27,28], in which the original threeobjective problem was resolved as two separate bi-objective problems. The proposed task is handled as a three-objective optimization task because this form of the problem fits better for real-life design tasks. The previous papers excluded asymmetrical solutions from the optimization. However, due to the non-linearity of the optimization problem, some asymmetrical solutions can be competitive with some symmetrical solutions. This 3D solution space is analyzed in this paper, and the analysis makes two other modifications to the parameter space of the problem. Firstly, the boundaries of the radii are changed. Secondly, the number of turns is varied, and the results of these three separate analyses are compared in the paper. The project files of the proposed analysis can be downloaded from the Artap project's homepage (https://github.com/artap-framework/artap/tree/master/ examples, (accessed on 1 June 2021)).

Formulation of the Problem
The goal of this optimization is to create a uniform field distribution in the control zone with a solenoid [12,23,24]. The solenoid is composed of 20 series of connected, singular turns, with a radial position varying from 5 to 50 mm. These turns have exactly the same size. The width of each turn is w = 1.5 mm, the height is h = 1.0 mm, and the prescribed DC current density is j φ = 2 A mm 2 . The flux density should be B(r, z) = (0, B 0 ) with B 0 = 2 mT in the controlled region. The model of the optimized coil is shown in Figure 1, where the green rectangles show the controlled region, and the yellow ones represent the different turns. The main difference from the original description of the TEAM 35 benchmark problem is that the full coil is modeled for asymmetrical calculations, not just the upper half of the solenoid (the description of the examined TEAM benchmark problem can be found at https: //www.compumag.org/wp/wp-content/uploads/2021/07/problem-35.pdf, (accessed on 1 September 2021)). The quality of the uniform magnetic field is assessed by using the point values of the magnetic field, which are evenly spaced among 100 points of the controlled region. Three different objective functions were used to measure the quality of the different solutions. The first one describes the uniformity of the magnetic field in the examined region; the z-component of the flux density is sampled in a 10 × 10 grid. The function takes the maximal difference from the intended flux density (B 0z = 2 mT) into consideration: The second function considers the robustness of a given solution. Let B(R) be the flux density values (2 rows, 100 columns) at a given input R. After that, the B + (X + ∆ξ) and B − (X − ∆ξ) vectors need to be computed, where ∆ξ = ±0.5 mm and B − (X − ∆ξ) represents the magnetic flux density change in the case of a 0.5 mm positioning error in a turn. F 2 can be defined in the following way: where . means the Euclidean norm and i represents the measurement point in the control region. The third function represents the mass or the DC loss of the coils. These quantities are proportional to the input. The summation of R is given by F 3 : where R j represents the radii of the separate turns. There are 20 different turns in the assembled coils, and their distance from the z-axis can be defined separately ( Figure 1).

Modeling and Optimization Frameworks
The coil model described above was modeled with two different FEM tools: Agros2D, an hp-adaptive FEM solver, and a widely used tool, FEMM, which was used to solve this magnetostatic problem [33,34]. The model was defined by the Adze-modeler (Figure 2), which allowed us to switch between the applied solvers with one command and connect the realized model withĀrtap (Ārtap is available for download from: http://www.agros2 d.org/artap/ (accessed on 1 June 2021) [35]). The workflow for the Adze-modeler can be seen in Figure 2. Different pieces of the geometry can be imported from different file types and can be exported to various FEM solvers by using the same description of the physical model. These parametric FEM models are generated by a function call fromĀrtap, which is an optimization framework for robust design optimization. It was developed within the Department of Theoretical Electrical Engineering at the University of West Bohemia in conjunction with a fully hp-adaptive FEM-solver: Agros Suite or Agros2D [28,36,37]. It provides a simple, general interface for facilitating the solution of real-life engineering design problems. The code contains evolutionary and genetic algorithms, wrappers for derivative-free methods, machine learning methods, and an integrated FEM solver. The goal of the realized multi-layered architecture is to separate the problem's definition from those of optimization algorithms and other artificial-intelligencebased methods and to provide automatic parallelization and database connection for the applied algorithms [5,35]. The image (a) shows the functionalities used and the workflow that was realized via the Adze-modeler. (b) The usage of the collision detection function, which automatically replaces the overlapping edges due to the optimization process.

Model Validation
The solenoid was simulated in two different ways in this paper. Firstly, the assumptions of symmetry were used, and only the upper half of the model was calculated with the FEM solvers. This model contains the first ten turns, and a Neuman boundary condition is applied at the z = 0 axis. This model is exactly the same as that used in the reference calculations [12]. However, the second model contains all 20 turns without using any assumptions of symmetry. The goal of this analysis is to let the optimizer select anti-symmetric solutions. These models were prescribed in Adze-modeler, which could export them for the FEMM or Agros2D models. The following test case was selected from Table 1 [12] to validate the correctness of our FEM models and the calculation of the objectives. Only the results of the symmetrical 10-turn model are presented in this paper, as we found the same results with the asymmetrical 20-turn solenoid model.
The simulations were performed with Agros2d and FEMM. In the case of Agros2d, an hp-adaptive mesh was set with 0.001% tolerance, while the mesh size was set to 0.5 mm in FEMM. The results and the settings used are summarized in Table 1. It can be seen that the resulting values are in a relatively good agreement with the reference calculations. The absolute value of the differences is a scale of magnitude smaller than the Earth's magnetic field, which can be a possible measurement limit in practice. The results of the F 3 objective calculations are not shown in Table 1 because the value of this function does not depend on the finite element models; it is simply calculated from the sum of the input vectors.

Sensitivity Analysis of the FEM Models
Another comparative analysis was made for calibration purposes. The analysis aimed to minimize the solution time and the computational demand of the optimization task by selecting the smallest mesh that was large enough to solve the task with the required accuracy. The required accuracy was 1% in the F 1 and F 2 metrics because it meant 20 µT in our problem, which is comparable with the Earth's magnetic field. Therefore, it is smaller than the precision of the measurement or the other neglected modeling details.
During the analysis, the Adze-modeler was used to convert and solve a randomly selected geometry for the different FEM tools. This model was solved with different mesh and solver settings ( Table 2.) in FEMM [33] and Agros2D [28,36]. The sensitivities of the F 1 and F 2 objective functions with the different mesh settings were compared in a selected geometry ( Figure 3). In Agros2D, the polynomial order (p-adaptivity) was set to 2 for all cases. The mesh refinement (h-adaptivity) was considered with different error indicator settings from 5% to 0.005% (Table 2). FEMM can only use first-order polynomials and does not have any adaptivity. It has a "Smart Mesh" feature that is turned on by default. It generates a dense-enough mesh with Triangle [38] to ensure accurate calculations, but it cannot be parameterized. We set up the same mesh size in all regions by turning off this feature, and this mesh size changed during the comparison process. The settings applied are summarized in Table 2.
The  The same solution is plotted in Figure 4 to visualize the differences in the solutions with the following settings: the error indicator was set to 1% in Agros2d and the smartmesh function was used in FEMM ( Figure 4). There is a significant difference in the point values of the magnetic flux density on the right side of the examined region (r = 5 mm). Figure 5 shows that these calculation errors have an effect on the objectives. F 1 considers the maximal difference from the expected value at a single point. These points can cause significant errors during the optimization. The sensitivity of the F 1 and F 2 functions to the mesh selection was examined. It is plotted in Figure 5, where the picture (a) justifies the above-mentioned assumption that the mesh selection has a significant (50%) effect on the values of F 1 .
The F 3 function is independent of the mesh selection, and it has a local minimum when all of the radii have the minimum value because it simply depends on the geometry of the coil. Agros2D was used for further calculations with the following settings (Table 3) because it clearly outperformed FEMM during the analyses. The error indicator was set to 1%. Using a more precise calculation seemed pointless because 1% of our target value (B o = 2 mT) was only 20 µT, which is smaller than the Earth's own magnetic flux density, which would affect the calculation results. If the final application needs more precise results, this effect should also be considered with magnetic and other geometrical simplifications.

Optimization of the Three-Objective Problem
The 35th TEAM benchmark problem was optimized as a three-objective optimization task. This approach is the first difference from the previously proposed solutions, where the whole problem was divided into two multi-objective optimization tasks [12,23]. The 2D segments of the three dimensional Pareto surface are plotted and examined in the following subsections and figures. All three objectives were considered during the optimization process because during the design of a product, all of these aspects should considered together. The symmetrical and asymmetrical solutions were optimized separately. The "symmetrical" model refers to the 10-turn geometry, where only the upper half of the solenoid is calculated during the optimization, while the asymmetrical 20-turn model makes it possible to optimize all of the turns separately. Both the "symmetrical" and "asymmetrical" models were calculated with Agros2D with the setup discussed above ( Table 3).
The optimization was performed with Artap while using the NSGAII algorithm [29]. In all of the cases, the maximum number of generations was 250 with 100 individuals. In the symmetrical cases, the optimization contained 10 independent variables, while in the asymmetrical cases, the problem contained 20 independent variables.
The following analyses were made for the three cases with the three different settings: (a) The radii of the first four and the last four turns varied from 1 to 50 mm, while the radii could be changed from 5.5 to 50 mm in the case of the other turns. (b) The radii of all of the turns could be changed from 5.5 to 50 mm. (c) The number of turns was reduced to 12, and all of the radii could be changed from 5.5 to 50 mm.
The results of these different optimization tasks are discussed in the following subsections.

Optimization of Case (a)
First of all, the three-dimensional Pareto surface after the optimization is plotted in Figure 6. It can be seen that the shape of this function is very spiky, and it is hard to localize one distinct optimum. There are many local optima that are close to each other, but most of them are very sensitive to the parameter changes. Optimizing the given solenoid for only one of the selected goal functions can easily lead to a non-robust solution. Figure 6. The plots (a-c) depict the shape of the optimized F 1 , F 2 , and F 3 objectives, (a-c) images plots the 3D surface of the objective function from different views.
After the optimization was performed, the results were sorted based on the values of the different objective functions, as can be seen in Figures 7-9. The different sortings indicate different priorities, but the optimization was performed with all three functions considered. For example, the F 1 sorting means that the last generation of solutions was sorted based on the F 1 value, and then the first 100 were selected. This sorting aims to show that the asymmetric solutions can be used just as well as the symmetric ones, and for some criteria, they even outperform the symmetric solutions.
Initially, the solutions from the three-dimensional Pareto plot that had the best values for F 1 were examined (Table 4). This means that these solutions have the best performance for the uniformity objective. As shown in Figure 7a, the F 1 values are in the range of 10-40 µT, which is within the measurement's error. This means that they can be considered equal. This holds on the F 2 axis as well. In terms of price, both variant groupings are in the same price range, but the symmetric setups have a slight advantage. Therefore, the best solutions are symmetrical ones, but there is no significant difference between the best symmetrical and asymmetrical solutions if we consider the uniformity of the solutions. Table 4. Solutions based on the F 1 sorting in case (a).  Secondly, the most robust solution was sought, and it was a result of an F 2 sorting ( Table 5). As shown in Figure 8a, the symmetric setups have an advantage, but they are in the region of the measurement error again, so they can be considered equally robust. In terms of uniformity, the difference is negligible. Regarding the price criteria, the symmetric setups can be slightly cheaper. Table 5. Solutions based on the F 2 sorting in case (a).  The goal of the third examination was to sort the results by their price (F 3 ) ( Table 6). It can be seen in Figure 9 that the symmetric setups are cheaper by 12-20%, but choosing one of them would be a sub-optimal solution, since the asymmetric setups perform significantly better in terms of accuracy and robustness. The radii of the Pareto-optimal results are plotted in Figure 10a. This violin plot shows the approximate shape of the symmetrical and asymmetrical solutions. The optimizer can set tiny values of the radii for the coils that are placed above or below the homogenized region in these solutions. Most of the Pareto-optimal solutions have at least one turn that has a radius ≤5.5. There is no big difference between the distributions of the radii for the symmetrical and asymmetrical cases. We can conclude that some asymmetrical results that can be competitive with the symmetrical ones exist. Table 6. Solutions based on the F 3 sorting in case (a). If all of the objective functions are considered, then the symmetrical and asymmetrical variants perform similarly (Table 7) because the difference between them lies in the region of tolerance. The fact that these solutions vary in the values of their radii implies that more than one solution exists for this problem. In Figure 10a, one can see the distribution of the values of the radii for the last 100 individuals with different setups. In Figure 11a, the symmetric and asymmetric setups are compared. In this optimization, all coils are free to move within their logical boundaries. Both setups converge roughly to the same range. The first and last four coils have tiny values of their radii, which makes them hard to manufacture.

Optimization with the Parameter Settings of Cases (b) and (c)
The two other optimization runs were carried out on the symmetric variant to examine the impact on the Pareto surface if we neglected the last 3-3 turn from the optimization (case (c)) or did not allow radii smaller than 5.5 mm to be selected (case (b)).
First, in case (b), the optimization was performed with the constraint 5.5 ≤ radii ≤ 50 for all coils, as given in the original benchmark problem. This solution excludes the best candidates in the solution space, which have small radii. In the second experiment, we examined if the turns mentioned above were cut out, as described in case (c). The main question during this experiment is that of if we can resolve the task with only a twelve-turn coil instead of the original twenty-turn one.
The results of the optimization in case (b), where the minimum radii of the turns had to be greater than 5.5 mm, are plotted in Figures 10b and 12. Figure 10b shows the distribution of the optimal radii in the two examined cases. The red plot shows the distribution of the radii in the previous symmetric case, while the blue plot shows the new constrained solution. As can be seen from the picture, all of the radii are greater in this case. Therefore, the small turns at the two ends of the coil can significantly improve the uniformity of the magnetic field in the homogenized region. In Figure 12, we can see that the goal function values are significantly worse than in the previous case.
In case (c), we solved both the symmetric and the asymmetric problems with 12 coils instead of 20, as shown in Figure 12. The results can be seen in Figures 12 and 13. If the 100 cheapest solutions are considered, then the reduced number of coils produces better results, especially with an asymmetric setup. In terms of cost, contrary to what one would intuitively expect, the reduced setups cost more than their counterparts by 12-20%. The reason for this small difference is that the twelve-turn variant contains turns with generally bigger radii. This difference is not significant if we compare the price of the coil with the solution of the original problem (constrained, Figure 10b). Regarding the best F 1 and F 2 solutions, the reduced setups perform worse in every way than the 20-coil setups.  The proposed results and the working version of the realized optimization problem can be downloaded from theĀrtap project's homepage (the proposed solutions are available for download from: http://www.agros2d.org/artap/, accessed on 1 June 2021) together with a semi-analytical solution, and this can be used to validate the performance of FEM-based multi-objective optimization tools.

Conclusions
A novel three-variable analysis of the multi-objective TEAM benchmark problem is proposed in this paper. The original benchmark problem contains two similar twovariable Pareto optimizations. Firstly, the proposed coil is optimized to produce a uniform magnetic field in the examined region, and the sensitivity of the coil should also be minimized [12]. Secondly, the uniformity of the magnetic field and the mass of the coil are optimized. In practice, these two problems should be handled together. In this paper, the three objectives mentioned above-the uniformity (F 1 ), sensitivity (F 2 ), and price (F 3 )are considered simultaneously. Before optimizing the coil, an FEM simulation was made using two different tools with different mesh settings to find the right FEM setup. It was found that F 1 and F 2 are sensitive to the mesh settings, especially on the right side of the examined region. The FEMM-based calculation of the default (smart mesh function) mesh can produce 100% error in calculating F 1 . It was found that the hp-adaptive solver was significantly faster and gave more accurate results. This solver was used to optimize with 1% in the error indicator, which produced an uncertainty of less than 20 µT in the results. This tolerance is acceptable because it is smaller than the measurement error or the Earth's magnetic field, which perturbs the results. This paper shows that this optimization problem is highly nonlinear, and nothing guarantees that only one optimal solution exists. We considered the asymmetrical cases in the solution and found that the symmetrical solutions always produced better solutions for F 1 and F 2 . However, the difference between those solutions is insignificant for F 1 and F 2 , as it is smaller than 20 µT.
Another interesting discovery is that the cheapest solutions are symmetrical setups, but they perform worse than the cheapest asymmetric ones for F 1 and F 2 . We reduced the number of turns from 20 to 12, and we found that the price of the coil was reduced by only about 12-20%, which was below the expectations. Further studies should be carried out to validate the proposed results by performing measurements on at least a single layout, and these measurement results can be used to benchmark the proposed results.