*Editorial* **Robust Design Optimization of Electrical Machines and Devices**

**Tamás Orosz 1,\*, David Pánek 2, Anton Rassõlkin <sup>3</sup> and Miklós Kuczmann <sup>1</sup>**


This article introduces a Special Issue (SI) that contains fourteen chosen articles from robust design optimization of electrical machines and devices. Optimization is essential for the research and design of electromechanical devices, especially electrical machines. Finding the optimal solutions may lead to cheaper, more economical products, faster and more efficient production, or more sustainable solutions. However, optimizing such a complex system as an electrical machine is a computationally expensive optimization problem, where many physical domains should be considered together. However, a good, practical design needs to consider the electrical device's design parameters; it should be insensitive to parameter changes or manufacturing tolerances. This Special Issue focused on papers showing how modern artificial intelligence (AI) tools can be used for robust design optimization of electric machines and electrical devices, how these tools can be benchmarked, or the correctness of the result validated.

The articles which are published in this special issue present the latest results of current research fields. Hopefully, the presented models and various application fields will provide useful information for researchers and professionals interested in these techniques themselves or who have other problems from different fields.

Testing and benchmarking the numerical tools for electromagnetic analysis is an important task. The Compumag Society provides openly accessible, challenging benchmark problems (TEAM problems) for testing novel numerical solvers. In [1], the authors deal with a solution of a robust design of a solenoid, and the test problem aims to search for the optimal shape of a coil, which ensures a uniform field 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 paper points out that if we are looking for designs with acceptable tolerances, not only symmetrical designs can be favoured. The paper points out 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.

A variety of electromechanical systems requires special techniques for optimization; each optimization is unique and focused on specific parameters aimed at performance improvement. In [2], a fast and accurate optimization tool is presented for optimal exploitation of permanent magnet synchronous machine with hairpin winding intended for transport applications. The focus of the optimization is maximizing power density and efficiency. As a benchmark case study, a surface-mounted permanent magnet synchronous motor designed for a student racing competition vehicle was considered. Several optimization steps are presented in the paper, and as a result, the main indexes, such as efficiency, volume power density, and power losses, were improved by 0.15%, 10.55%, and 3.4%, respectively.

**Citation:** Orosz, T.; Pánek, D.; Rassõlkin, A.; Kuczmann, M. Robust Design Optimization of Electrical Machines and Devices. *Electronics* **2022**, *11*, 1427. https://doi.org/ 10.3390/electronics11091427

Received: 20 April 2022 Accepted: 25 April 2022 Published: 29 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Another study on permanent magnet synchronous machines [3] is focused on 8-pole 9-slot and 8-pole 12-slot machines and considers rotor eccentricity. Authors have performed a magnetic field analysis using an analytical method and the torque characteristics for benchmark machines. The optimization is based on perturbation and electromagnetic theories. In both cases, two models (slotted and simplified without slots) were analyzed using the FEM. The research work confirms that the slotted model can obtain similar results, even if the magnetic flux density is predicted without slots and the back-EMF is derived using it. The results obtained using the analytical method are compared with the FEM and experimental results.

The next paper in SI [4] aims to investigate the reconfigurations of rotor flux barriers for a five-phase permanent magnet assisted synchronous reluctance machines. That type of electrical machine is relatively new on the market. However, they are gaining popularity in industrial applications due to their electromagnetic characteristics (robustness, torque/power density, performance, etc.). In this research work, a Lumped Parameter Model conducted to a 2D FEM was applied to the proposed permanent magnet assisted synchronous reluctance motor models under the steady-state condition. Based on the FEM results, the maximum torque, minimum cogging torque, and minimum torque ripples were achieved. As a result, the optimal model of the electrical machine operates at high-performance values with desirable values of line-to-line back-EMF and air-gap flux density.

Wind generators are integrated with electrical machines that require a reliable operation. However, the increasing use of non-linear loads introduces undesired disturbances that may compromise the integrity of the electrical machines inside the wind generator. Ref. [5] proposed a five-step methodology for power quality disturbance detection in grids with the injection of wind farm energy. The proposed method is validated using a set of synthetic signals and is then tested using two different sets of real signals from an IEEE workgroup and from a wind park located in Spain.

In [6], the dielectric properties of a Bi3−xNdxTi1.5W0.5O9 material is investigated. Many recent studies showed that replacements of atoms in A and also in B-positions of an AP's crystal lattice led to a change in the structure, the dielectric properties and significantly influenced the polarization processes in this compounds.The dependences of the relative permittivity *-*/*-*<sup>0</sup> and the tangent of loss *tanδ* at different frequencies on temperature were examined in this paper, together with the piezoelectric properties of the material.

Nanotechnology provides an effective way to upgrade the thermophysical characteristics of dielectric oils and creates optimal transformer design. The properties of insulation materials have a significant effect on the optimal transformer design. Ester-based nanofluids (NF) are introduced as an energy-efficient alternative to conventional mineral oils, prepared by dispersing nanoparticles in the base oil. Ref. [7] presents the effect of graphene oxide and TiO2 nanoparticles on the thermophysical properties of pure natural ester and synthetic ester oils with temperature varied from ambient temperature up to 80 °C. A range of concentrations of graphene oxide (GO) and TiO2 nanoparticles were used in the study to upgrade the thermophysical properties of ester-based oils. The experimental results show that nanoparticles have a positive effect on the thermal conductivity and viscosity of oils which reduces with an increase in temperature

Low voltage cables are widely used in nuclear power plants and photovoltaic generators. In the case of nuclear power plants the low voltage cables link the system components with the controlling and monitoring instrumentation and control equipment and supplying power to the devices. During the service period these cables are exposed to a wide range of stresses: high-temperature, radiation, mechanical stresses, etc. Since the proper function of the low voltage cables is essential for these power plants' continuous and reliable operation. In [8], the authors examine the effect of mechanical stresses during the aging procedure of these cables, it shows that the Shore D hardness was also higher on the thermo-mechanically aged samples. These findings show the combined aging has a higher impact on the insulation properties. Hence, involving the mechanical stress in

the aging procedure of cable qualification enables the design of more robust cables in a harsh environment.

New materials and manufacturing technologies have influence also on design optimization process of electromechanical devices. In [9], the authors present an optimization of a additively manufactured permanent magnet coupling. Two approches are introduced time-consuming Genetic Algorithm method and faster Taguchi method. The research work analyze the abilities of compared methods within the optimization of studied coupling with minimization of volume and maximization of transmitted torque as objectives. Taking into account that resulting optimal geometry (the clutch volume is reduced by 17%) and characteristics (magnetic torque density is enhanced by more than 20%) achieved by compared methods are nearly identical, the Taguchi method is found to be more time-efficient and effective within the considered optimization problem. The permanent magnet coupling was manufactured and simulation results were validated using an experimental setup.

Model predictive current control has recently become a powerful advanced control technology in industrial drives. In [10], the authors proposed a computationally efficient calcualtion of the current prediction control for synchronous reluctance motors. The porposed methodology can reduce the computational cost by a merging the predictive current control model with a simple hysteresis current control. Therefore, only four voltage vectors should used to predict the current and evaluate the cost function. The proposed methodology can reduce the computation cost of a classical predictve current model by about 20%.

Linear motors are a special type of electrical machine that requires special attention due to nonlinearities caused by side effects. The authors of the paper [11] propose a modified dynamic equivalent circuit model for a linear induction motor. A proposed model considering both longitudinal (speed-dependent) end effect based on conventional Duncan's approach and transverse edge effect investigated by using additional correction factors. In addition, the field-analysis method is used to include the typical linear motors iron saturation effect, the skin effect, and the air-gap leakage effect. Model simulation results show a good agreement between field analysis and FEM estimation of the electrical parameters. Moreover, to validate the proposed paper method, 3-D FEM was employed. Thrust-velocity characteristic of studied linear induction machines shows that the proposed method provides more precision as compared to Duncan's model.

An investigation of linear induction motors in SI continues with work by Zhang et al. in [12]. An improved equivalent circuit model of double-sided linear induction motors that takes into account the linear motor skin effect and the nonzero leakage reluctance of the secondary, longitudinal, and transverse end effects into consideration is proposed. The proposed equivalent circuit is presented described in detail and highlights the modification in comparison with the traditional equivalent circuit with longitudinal and transverse end effects. 3D FEM is used to verify the proposed equivalent circuit model under varying air gap width and frequency. The results show that the equivalent circuit model that takes into account only the longitudinal end effect considered, and the model considered with both longitudinal and transverse end effect have more than 11% errors with the FEM simulation results in the slip range, while the errors between the value of proposed equivalent circuit and simulation are less than 5%.

There is a great potential in small satellite technology for testing new sensors, processes, and technologies for space applications. The design of their receiving antennas for their ground stations needs a careful design to establish stable communication. Paper [13] shows an interesting solution to the antenna design problem with the antenna array technology. This novel approach can have many advantages over parabolic antennas. From a mechanical point of view, it does not require the design and maintenance of the drive system, which sets the azimuth and the elevation angles. Such systems have a simpler feeding network that cannot be disconnected during the connection time. These tools are insensitive to the moisture and weather conditions during the mission. Moreover, with a pattern reconfigurability algorithm, they can support multi-task missions. This work is

motivated by the design of an antenna array for a future rotatorless base station for the VZLUSAT group of Czech nano-satellites.

The advancement of a device like an insulated core transformer involves the optimization of several parameters. Special attention must be paid to parameters that affect the uniformity of disk output voltage. In the paper, Ref. [14], the accuracy of the FEM model was verified by comparing test data of the insulated core transformer prototype with the simulation results. Particle Swarm Optimization algorithm was implemented for the design parameters (including the number of secondary winding turns and the compensation capacitance) optimization of dummy primary winding. The optimization results presented in the research work show that the maximum non-uniformity of the disk output voltage is reduced from 11.1% to 4.4% from no-load to a full load for a 200 kV/20 mA for an insulated core transformer prototype. The proposed method improves the performance of the insulated core transformer high voltage power supply and cuts down the design time.

**Funding:** The research work by Anton Rassõlkin has been supported by the Estonian Research Council under grant PSG453 "Digital twin for propulsion drive of autonomous electric vehicle".

**Acknowledgments:** For this valuable collection of research works focuses on optimization of electrical machines and devices, the Guest Editors are thankful for all authors who submitted their manuscripts for this SI and congratulate them on publishing their research works with MDPI Electronics. This SI edition would not be possible without the the Academic Editors and all reviewers, our gratitude for their important work. Last but not least, we would like to thank the MDPI team for their support of this SI.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Robust and Multi-Objective Pareto Design of a Solenoid**

**Krisztián Gadó and Tamás Orosz \***

Montana Knowledge Management Ltd., 1111 Budapest, Hungary; gado.krisztian@montana.hu **\*** Correspondence: orosz.tamas@montana.hu

**Abstract:** 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 ensures a uniform field 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 insignificant 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.

**Keywords:** optimization; electrical machines; design optimization; finite element method

**Citation:** Gadó, K.; Orosz, T. Robust and Multi-Objective Pareto Design of a Solenoid. *Electronics* **2021**, *10*, 2139. https://doi.org/10.3390/ electronics10172139

Academic Editor: Davide Astolfi

Received: 2 August 2021 Accepted: 31 August 2021 Published: 2 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. 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–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–15]. A similar design task should be solved during the design and optimization of induction heating or induction brazing processes [16–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)).

#### **2. 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<sup>φ</sup>* = 2 *<sup>A</sup> mm*<sup>2</sup> . The flux density should be *<sup>B</sup>*(*r*, *<sup>z</sup>*)=(0, *<sup>B</sup>*0) with *<sup>B</sup>*<sup>0</sup> = 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.

**Figure 1.** The examined geometry in an axisymmetric arrangement. The green area shows the controlled region in which the magnetic field is considered. Every single turn of the coil is denoted by yellow rectangles; their radii are optimized turn by turn.

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*0*<sup>z</sup>* = 2 mT) into consideration:

$$F\_1 = \max |B\_{0z} - B\_{iz}|, i = 0, \dots, 99. \tag{1}$$

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*<sup>2</sup> can be defined in the following way:

$$F\_2 = \max\{ \left| \left| \mathbf{B}\_i^+ - \mathbf{B}\_i \right| \right| + \left| \left| \mathbf{B}\_i - \mathbf{B}\_i^- \right| \right|, 0 \le i \le 99 \}, \tag{2}$$

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:

$$F\_3 = \sum R\_j.\tag{3}$$

where *Rj* 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).

#### *2.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 Artap ( ¯ Artap 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 Artap, 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].

**Figure 2.** 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.

#### *2.2. 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 following vector contains all of the design variables for the selected test case (Table 1, [12]):

$$X\_1 = [6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. \tag{4}$$

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*<sup>3</sup> 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.


**Table 1.** The results of the validation run for the *X*<sup>1</sup> input compared to the reference values of *F*<sup>1</sup> and *F*2.

#### *2.3. 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*<sup>1</sup> and *F*<sup>2</sup> 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*<sup>1</sup> and *F*<sup>2</sup> objective functions with the different mesh settings were compared in a selected geometry (Figure 3).

**Table 2.** The settings applied for the FEMM- and Agros2d-based calculations.


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 *F*<sup>1</sup> and *F*<sup>2</sup> functions were calculated from the point values of the magnetic flux density of the control zone; hence, these functions can be sensitive to the numerical errors of the point values of the magnetic flux density calculations. It can be seen from Figure 3c,d that the *Br* and *Bz* components are sensitive to the mesh applied at the top right corner of the control region (r = 5.0 mm, z = 5.0 mm). There is a huge difference in the convergence speed of the two different FEM solvers applied. The smallest FEMM mesh contains about <sup>5</sup> × <sup>10</sup><sup>5</sup> nodes, whereas the Agros2D converges to the result in both *Br* and *Bz*. The difference in the case of the radial component is not significant. However, in the case of the axial component, the difference is significant (about 40%).

**Figure 3.** The pictures (**a**,**b**) show the geometry examined and the resulting flux distribution with the Agros2D (**a**) and FEMM software (**b**). The pictures (**c**,**d**) show the convergence of the radial and the axial component of the flux density in the selected (*r* = 5 mm, *z* = 5 mm) point.

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*<sup>1</sup> 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*<sup>1</sup> and *F*<sup>2</sup> 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*<sup>3</sup> 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 (*Bo* = 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.

**Figure 4.** The pictures (**a**,**b**) show the calculated flux density values in the upper half of the examined region with the two types of software compared: Agros2d (**a**) and FEMM (**b**). The third picture (**c**) shows the dependence of the *Bz* values on the meshing properties.

**Figure 5.** The convergence of the objective functions (*F*<sup>1</sup> and *F*2) in terms of the number of nodes with and without hp-adaptivity. (**a**) The image plots the dependence of *F*1, while (**b**) plots the dependence of the *F*<sup>2</sup> function on the number of nodes.

**Table 3.** The settings used for Agros2D during the optimization.


#### **3. Results and Discussion**

*3.1. 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:


The results of these different optimization tasks are discussed in the following subsections.

#### *3.2. 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*<sup>3</sup> 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*<sup>1</sup> sorting means that the last generation of solutions was sorted based on the *F*<sup>1</sup> 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*<sup>1</sup> were examined (Table 4). This means that these solutions have the best performance for the uniformity objective. As shown in Figure 7a, the *F*<sup>1</sup> 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*<sup>2</sup> 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*<sup>1</sup> sorting in case (a).

**Figure 7.** The image shows the first 100 individuals when they are sorted based on uniformity (*F*1) in case (**a**). The blue dots show the symmetric setups, while the asymmetric solutions are depicted in red. (**a**) The distribution of individuals on the *F*1–*F*<sup>2</sup> axis. (**b**) The distribution on the *F*1–*F*<sup>3</sup> axis. (**c**) The distribution on the *F*2–*F*<sup>3</sup> axis.

**Figure 8.** The first 100 individuals from the three-dimensional Pareto front were sorted based on their robustness (*F*2) in case (**a**). The blue dots show the symmetric setups, while the asymmetric solutions are depicted in red. (**a**) The distribution of individuals on the *F*1–*F*<sup>2</sup> axis. (**b**) The distribution on the *F*1–*F*<sup>3</sup> axis. (**c**) The distribution on the *F*2–*F*<sup>3</sup> axis.

Secondly, the most robust solution was sought, and it was a result of an *F*<sup>2</sup> 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*<sup>2</sup> sorting in case (a).

**Figure 9.** The image shows the distribution of the 100 cheapest individuals (*F*3) in case (**a**). The blue dots show the symmetric setups, while the asymmetric solutions are depicted in red. (**a**) The distribution of individuals on the *F*1–*F*<sup>2</sup> axis. (**b**) The distribution on the *F*1–*F*<sup>3</sup> axis. (**c**) The distribution on the *F*2–*F*<sup>3</sup> axis.

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*<sup>3</sup> 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.

**Figure 10.** The distribution of the radii of the last 100 unsorted individuals. (**a**) Symmetric and asymmetric setups. All coils are allowed to move freely within the logical boundaries. (**b**) Only the symmetric setups where the coils are first allowed to move freely are shown (red); then, they are constrained in the 5.5–20 mm region (blue). (**c**) The distributions of the radii for the symmetric and asymmetric setups using only 12 coils. All coils were constrained to move within the 5.5–20 mm region.

**Figure 11.** The last 100 individuals if the last generation is not sorted based on any of the conditions. These solutions are lying on the Pareto front. The blue dots show the symmetric setups, while the asymmetric solutions are depicted in red. (**a**) The distribution of individuals on the *F*1–*F*<sup>2</sup> axis. (**b**) The distribution on the *F*1–*F*<sup>3</sup> axis. (**c**) The distribution on the *F*2–*F*<sup>3</sup> axis.

#### *3.3. 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*<sup>1</sup> and *F*<sup>2</sup> solutions, the reduced setups perform worse in every way than the 20-coil setups.


**Table 7.** Last individual in case (a).

**Figure 12.** The comparison of the 100 cheapest symmetric solutions. With the blue color, all 20 coils are free to move, and with the red, only 12 coils are used, and they are constrained to move within the 5.5–20 mm region.

After analyzing the solution plot, the search range of the optimized turn radii can be reduced from 5.5–50 to 5.5–20 because the greater radii are not represented in the solution of the three objective problems.

**Figure 13.** The image shows the last 100 symmetric solutions with various setups. The case where all 20 coils were free to move is in red. The blue color shows when all 20 coils are constrained in the 5.5–20 mm region, and finally, the green color depicts the case where only 12 coils are used, and they are constrained to the 5.5–20 mm region. (**a**) The distribution of the various setups on the *F*1–*F*<sup>2</sup> axis, (**b**) on the *F*1–*F*<sup>3</sup> axis, and (**c**) on the *F*2–*F*<sup>3</sup> axis.

The proposed results and the working version of the realized optimization problem can be downloaded from the Artap 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.

#### **4. 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*<sup>1</sup> and *F*<sup>2</sup> 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*<sup>1</sup> and *F*2. However, the difference between those solutions is insignificant for *F*<sup>1</sup> 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*<sup>1</sup> 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.

**Author Contributions:** Conceptualization, T.O.; methodology, T.O., K.G.; software, K.G., T.O.; validation, K.G., T.O.; writing—original draft preparation, T.O., K.G.; writing—review and editing, T.O.; visualization, K.G., T.O.; supervision, T.O.; project administration, T.O.; funding acquisition, T.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

