A Multi-Objective Finite-Element Method Optimization That Reduces Computation Resources Through Subdomain Model Assistance, for Surface-Mounted Permanent-Magnet Machines Used in Motion Systems

Surface-mounted permanent-magnet (PM) motors are widely used in motion control systems because of their high peak torque-to-inertia ratio. These motors exhibit high magnetic saturation to produce high peak torque. A precise finite-element model (FEM) is needed to optimize these motors. It requires high computation power, especially for multi-objective optimizations. On the contrary, subdomain models require low computing power but are inaccurate when magnetic saturation occurs. To solve this problem, we use a subdomain-assisted FEM optimization method (SAFEMOM) that combines subdomain models for preoptimization and FEM for refined optimization. We provide the quantitative measurement methods to compare SAFEMOM with FEM-only optimization. If the constraints in the magnetic fluxes in the subdomain part of SAFEMOM are based on the actual saturation value of the materials, then the advantage of SAFEMOM is not significant. The machines in this study use non-linear material with a knee point at 1.3 T and hence show heavy magnetic saturation above 1.5 T. In that case, contrary to what could be intuitively thought, we need to increase the magnetic flux density limit to 2.6 T - 3.0 T in SAFEMOM to have a significant advantage. SAFEMOM reduces about 80% of computing time to obtain a slightly better convergence than the one using FEM only. Also, if limited computing resource is allowed, SAFEMOM gives an error reduced by a factor of nearly eight compared to the optimization error using FEM only. Those results are validated on a family of surface-mounted permanent-magnet machines with a wide range of design parameters.


I. INTRODUCTION
Surface-mounted permanent-magnet (SMPM) motors have the advantage of a high ratio between torque and inertia, as well as simple manufacture. They are used for many The associate editor coordinating the review of this manuscript and approving it for publication was Philip Pong . applications, including many motion control systems, which is the background of the work we present here.
Each motion system needs to be optimized individually depending on its application. Much work has been published on optimizing electrical machines using single-objective optimization algorithms, such as genetic algorithms [1] and space-mapping algorithms [2], [3]. However, a Pareto front VOLUME 11, 2023 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ that shows the optimal designs for a few key features can often be of great use for those who make complex designs using those motors. Therefore, a multi-objective optimization of the motor is needed, which can bring a Pareto front that gives the designer an intuitive understanding of the relationship between the different performance parameters [4]. This Pareto front can also be used for system optimization [5]. For this purpose, previous research has used different types of multi-objective optimization algorithms: differential evolution optimization algorithms [4], [6], particle swarm optimization algorithms [7], [8], the strength Pareto evolutionary algorithm 2 [9], and the nondominated sorting in genetic algorithm II (NSGA-II) [10]. Whatever optimization algorithm is used, a considerable computation burden is needed to do precise multi-objective optimization of complex non-linear systems such as SMPM motors. Much work has been done using FEM to design and optimize motor structures [11]. FEM is a technique that is used for a wide range of topologies such as SMPM machines [12], [13], [14], [15], [16], permanent-magnet (PM) machines with inserted PMs [17], [18], [19], [20], [21], machines combining inserted PMs and PM-assisted synchronous reluctance motors [22]. Also, FEM has been used widely for axialflux permanent-magnet motors [23]. Because of the computation load related to FEM, a field reconstruction method was designed to reduce the computational load [24], [25], [26]. Nevertheless, this approach brings a high error when magnetic saturation occurs.
The problem is that many calculations are needed to create those Pareto fronts, which is very costly in terms of computing power. If multiple objectives need to be optimized simultaneously, it is problematic. It is even more problematic when the non-linear finite-element method (FEM) needs to be used to achieve a high-precision multi-objective optimization process. Indeed, because of the non-linearity, its computation burden is prohibitive. Hence, a reduction of the calculation burden of the optimization procedure is needed.
Previous research provided many kinds of analytical models to reduce the computation burden, including the linear relative permeance model [27], the linear complex permeance model [28], or the linear subdomain model [29], [30], [31], [32]. The linear subdomain model has been developed for many configurations. Many structures and phenomena, including dynamic effects, can be considered in those models [33]. Subdomain models have a wide application range [34], including structures involving dual rotors [35]. Recently, subdomain models have been improved to deal with structures containing some non-orthogonal boundaries [36].
It has been shown that the linear subdomain model is generally superior to the linear relative permeance model and the linear complex permeance model in terms of accuracy [37] and is also very fast. The accuracy is excellent for a system with linear magnetic features, but when magnetic saturation is present, the accuracy of the linear subdomain method decreases. There have been attempts to include magnetic saturation directly in the subdomain models, but the computational burden is high [38]. Much work has also been published using linear subdomain models to model and optimize structures [39], [40]. Nevertheless, as those models are linear, they lead to incorrect structures. Indeed, motors used in motion control systems need a high peak torque, which leads to high magnetic saturation.
It is natural to use a linear subdomain model for preoptimization first, and then a non-linear FEM to further increase the accuracy of the results. This approach can reduce the calculation resource for optimizing SMPM machines in motion systems. This technique is named the subdomain assisted FEM optimization process (SAFEMOM). If no constraint is imposed on the subdomain model, the designs resulting from the preoptimization using the subdomain model are far from the optimum. This is because magnetic saturation is absent from the linear subdomain model. In the research presented below, we found out that if constraints are imposed to forbid the magnetic flux density from being greater than the saturation point, the subdomain model does not provide any helpful assistance in the optimization of machines used in motion control systems.
Therefore, this article discusses a combination of the subdomain model and FEM to solve this problem properly. We compare optimization done through SAFEMOM and through FEM only. This comparison is made through quantitative measurement techniques. To the best of our knowledge, we have not seen any discussion of this comparison, nor quantitative measurements of convergence of such methods. We have not seen in-depth discussions with quantitative results on how linear models should be properly constrained to be of use for a non-linear optimization process. Also, we have not seen any article showing how this combination of a subdomain model and a FEM can achieve multi-objective optimization, which shows high precision and reduced computation burden.
In this article, we first define the optimization problem in section II. In section III, we discuss the framework and models for optimization. In section IV, we give a quantitative technique to measure optimization convergence, and we also give the SAFEMOM results. The SAFEMOM given in this article is compared to FEM-only optimization techniques. Finally, in section VI, we summarize the key contributions.

II. DEFINITION OF THE OPTIMIZATION PROBLEM
As the background of the article is motion control systems, the speed of the motor is relatively low. Such applications use SMPM machines (structure shown in Fig. 1) as their speed is relatively limited. These applications need a significant acceleration and deceleration to allow a displacement from one point to another to be done quickly. For the displacement to be repeated between two points at a high frequency, the SMPM motor is optimized using the following three objectives. 1) The peak torque is maximized. Therefore, minus peak torque T is minimized. 2) The inertia of the rotor I  is minimized, to allow high accelerations and decelerations. As the speed is low, the power losses in copper represent more than 98% of the losses. Therefore, for thermal and efficiency reasons 3) the power losses L in the copper are minimized. Each design p is identified by its three objectives p = (L, I , T ) and hence a point in the objective space.
In the optimization, the three following constraints are used. 1) For cost reasons, the volume of PM is limited. 2) The PM should not demagnetize. 3) Joule losses are limited to a maximum value to avoid overheated designs. Joule losses are also limited to a minimum value, as designs with lower values do not provide enough torque. The values are given in Table 1 and Figs 2 and 3.
The following eight continuous variables are optimized: 1) the current density, 2) the thickness of the PMs, 3) the PM opening angle ratio, 4) the tooth-tip opening ratio, 5) the tooth length, 6) the stator yoke thickness, 7) the thickness of toothtip 8) and the width of teeth. The ranges are shown in Table 2.
In the optimization, the PM (N48SH) at 80℃ is modeled according to the magnetization curve shown in Fig. 2. The ferromagnetic material is M19 steel. Its properties are shown in Figs 3 and 4. The knee point of its B-H curve is at 1.3 T. The knee point is calculated, based on IEC 61869, as the point at which an increase of 50% of the magnetic field intensity leads to an increase of 10% of the magnetic flux density. In Fig. 3, the magnetic flux density of the M19 steel is  represented, and its value for low magnetic field intensity can be clearly seen. In Fig. 3, we can see the magnetic saturation of the M19 steel beyond 1.5 T more explicitly. The outer radius and active length are fixed to keep the same cooling surface area. The inner radius is fixed because the mechanical design of the bearing is fixed. The air gap is assumed to be 0.8 mm. The pole-pairs and slots combination are also fixed during the optimization. The optimization is done using two-dimensional modeling for the magnetic modeling. Magnetic end effects are not considered for the torque calculation, and only the active length is taken for the inertia calculation. Therefore, torque and inertia are proportional to the active length. For the copper loss, we consider the end-winding influence. The values of the motor fixed parameters are given in Table 1.
The optimization is highly dependent on magnetic saturation. Hence, there is a need for quick and precise optimization that considers magnetic saturation. This optimization is done through the framework shown in the next section.

III. THE FRAMEWORK AND MODELS FOR THE OPTIMIZATION A. FRAMEWORK OF SAFEMOM
The SAFEMOM we propose is a hybrid optimization with two steps: the preoptimization using the subdomain method and the refined optimization using FEM. More details about the subdomain method and FEM are given in section III-B. The optimization process is shown in Fig. 5. We chose a   population of 28, which means that 28 designs are computed synchronously using an evolutionary optimization algorithm.
In the first step, the linear subdomain model is used to approach optimal designs quickly. As the calculation time of the subdomain method is extremely short, we take 300 generations in the preoptimization process. As the subdomain model is based on iron with infinite permeability, we need to give additional constraints to allow the algorithm to converge closer to the optimal point. In Fig. 6, the different lines indicate where the fluxes are limited due to saturation. As the magnetic vector potential A is known in the different points shown in the figures, the magnetic flux between these two points is known. The flux going through the rotor yoke between points 1 and 2 is known. Similarly, the flux in the tooth tip between points 3 and 4, the flux in the tooth going between points 5 and 6, and the flux in the stator yoke between points 7 and 8 are all known. For each of those four sections, the average magnetic flux density is constrained  below a given density limit. If the subdomain model gives a magnetic flux higher than the limit, the calculated design is removed from the design space. A more detailed discussion of the magnetic flux density limit value is given in section IV. In this article, every time we mention a magnetic flux density limit, it is a limit imposed only in the subdomain and not in the FEM part of SAFEMOM. Indeed, the FEM part is non-linear. Therefore, it takes magnetic saturation into account directly.
In the second step, based on the preoptimal results of the first step, a non-linear FEM is used to refine the optimization. As the non-linear FEM consumes much more computing power, the number of generations is restrained to 50. This restriction is reasonable, as the starting point of the non-linear FEM optimization is not too far from the optimum.
The calculation time required for the optimization is shown in Table 3, using AMD 3700X CPU and using an eight-core parallel computer. The subdomain computation time is less than 0.6% of the total computation time.
As we will see later, the fact that the subdomain model is quicker than FEM is not a sufficient condition to allow the reduction of the computation time in an optimization process, nor is it a sufficient condition to increase the accuracy in a given time. Indeed, a key question that needs to be solved is how to constrain the analytical model in a way that helps to increase the optimization speed and to allow it to get closer to the Pareto front. Before discussing those constraints, we give a short overview of the magnetic models used in the optimization.

B. THE MAGNETIC MODEL USED IN THE OPTIMIZATION 1) SUBDOMAIN MODEL
The subdomain model relies on the following assumptions. 1) The permeability of iron is infinite. 2) The magnetization of the PMs is parallel.
3) The space between the PMs is assumed to have the same permeability as the one of the PMs. 4) No end effect is taken into account. 5) The teeth are assumed to have both sides that are radial, whereas in reality, both sides are parallel, as shown in Fig. 7. As the area of the slot in the model is reduced compared to the reality, the current is increased accordingly. The model is based on the resolution of Laplace's equation expressed in the cylindrical coordinate system [41], in the air (layers 2 and 3 of Fig. 1) Poisson's equation in the slots (layer 4 of Fig. 1) and Poisson's equation in the PMs (layer 1 of Fig. 1) where A is the magnetic vector potential along z-axis, B R is the remanent flux density, J ex is the excitation current density, ϕ is the angle in the cylindrical coordinate system, and r is the radius in the cylindrical coordinate system. A Fourier series is used to apply these differential equations to each harmonic of the magnetic vector potential. As those are second-order differential equations, two constants must be solved for each harmonic of the magnetic vector potential in each layer or sector. The values of those constants are found using the different boundary conditions. The detailed resolution procedure is given in [33] and [34]. In this form, the subdomain model does not include magnetic saturation. Therefore, saturation is considered by removing designs with a magnetic flux density higher than the imposed limit in the different areas mentioned in Fig. 6.

2) FEM MODEL
We use Matlab to control a FEM software called FEMM, which allows us to obtain 2D FEM magnetic field results. In a typical calculation, the mesh has 7083 nodes and 12858 elements (Fig. 8). The mesh size has been decreased by using symmetry and further reduced where it does not impact the results significantly. Further calculation time has been decreased through periodic air-gap boundary conditions that remove the need to re-mesh the machine between two simulations with different rotor angles. This model considers magnetic saturation without any need for additional constraints, but this is at the cost of the computing time.

C. THE SELECTION OF THE OPTIMIZATION ALGORITHM
We can divide the computational time the evolutionary optimization algorithm needs into two parts. Part A is the computation time the model requires, in our case, non-linear FEM or subdomain, to give values for the objective function. Part B is the time used by the algorithm to define which designs need to be calculated in the next generation.
For the optimization of SMPM motors, more than 99% of the computation resources are used to calculate part A. For this reason, the question of the complexity of the optimization algorithm itself is not critical in our choice. We chose the Bi-criterion evolution optimization (BCE-IBEA [42]) as the optimization algorithm, which is known to have a relatively high computing time required for part B, compared to other evolutionary optimization algorithms, like the NSGA-II or NSGA-III [43], [44]. It combines the indicator-based and non-dominated optimization algorithm advantages. It is based on ideas coming from the strength Pareto evolutionary algorithm 2 [9] and the indicator-based evolutionary algorithm (IBEA) [45], the latest being inspired from the nondominated sorting in genetic algorithm II (NSGA-II) [10]. We chose it because it gives a better diversity of the results in the Pareto front [46], which is a highly desirable feature, as the resulting data are used in a Gaussian process regression (GPR) [47].

IV. THE MEASUREMENT OF THE OPTIMIZATION RESULTS
As we have suggested a new framework for the optimization of surface SMPM machines, we need to have ways to evaluate the quality of the results as a function of the different magnetic flux density limits used in the subdomain part of the optimization process. The results need to be evaluated quantitatively and fairly. We discuss here mainly two indicators.
1) The spread of the results describes how the obtained Pareto surface covers a wide range of parameters (objectives). It is calculated through the hypervolume indicator. It tells us how much hypervolume there is between the Pareto points and the reference point, as shown in Fig. 9. Hence a more extensive range of objectives increases the volume. Also, as we see in the figure, a better optimization results in a bigger hypervolume, as the goal is to minimize the objectives.
2) The convergence of the results describes the distance between the calculated non-dominated solutions of every generation and the approximate Pareto surface, which is defined in section IV-A1. It is calculated through the modified generational distance (GD+).
The GD+ tells us how far a set of points is from the approximate Pareto surface. If the point is more optimal than the approximate Pareto surface (objectives that are lower than those of the Pareto surface), the value of the distance between the point and the approximate Pareto front is taken to be zero for the calculation of the GD+, as shown in Fig. 10. This can happen because of the computation noise of the models, or if the approximate Pareto surface is slightly suboptimal. The spread should be maximized, which means that the hypervolume indicator should be maximized. Also, a good convergence should be obtained, which means that the average GD+ should be minimized.
The traditional way to assess both spread and convergence is through a visual check of the different optimization points. In Fig. 11, we see the result of our optimization process. In many papers, such figures are used to visually assess the convergence and spread of the optimization. This approach has clear limits, as it is not quantitative. This approach is even more problematic if the Pareto front is embedded in a four-dimensional space. Therefore, our quantitative approach, which also works in higher dimensions, is of critical importance.
Through the genetic algorithm, multiple designs are obtained. Each design p is identified by the three values (L, I , T ), which are the objectives of the multi-objective optimization. At the same time, each design's parameters (dimensions, material properties, current) are stored and associated with the design p. The different designs are grouped to form a dataset.  The dataset is processed to measure the hypervolume and GD+ indicators correctly. Two filters are used for this purpose. They are defined in more detail in sections IV-C1 and IV-C2. Their key contribution and role are given here: Firstly, the Niche-radius-based filter (NRBF) does mainly two things. It keeps only non-dominated designs, and then it reduces the size of the dataset of the non-dominated designs to obtain a better distribution of the points. As those points are used for a GPR, this distribution is essential for the confidence interval calculated by the GPR. More details on this filter are given in section IV-C1.
Secondly, the Robust filter removes unreliable designs from the dataset. It removes noise on the indicators and allows us to obtain better measurements. More details on this filter are given in section IV-C2.

A. SPREAD CALCULATION 1) CALCULATION PROCESS
A general diagram of the whole computation that includes this pre-processing for calculating the hypervolume indicator is shown in Fig. 12. An intuitive representation of the hypervolume is represented in Fig. 9. First, we start with the complete dataset. This dataset includes different generation results taken during the optimization process. In this calculation, some designs are kept in the next generation, leading to an accumulation of points in the same area and a bad distribution of the results in the objective space around the Pareto front. Some areas have many points, and some others do not. This distribution is problematic, as these data are used for a GPR. This problem is solved by using the NRBF defined in section IV-C1.
As evolutionary optimization involves a random process, it induces some results far from the other results, which mislead the calculation of the hypervolume. These results are typically in an area of low distribution of points. This low distribution of points has an impact on the GPR. Indeed, the GPR has the great advantage of including both prediction values and also a confidence interval for the predicted values. This confidence interval relies on the distribution. Because of this, the points of low distribution have low confidence, as calculated with the GPR. These points are hence removed for the calculation of the hypervolume. The robust filter defined in section IV-C2 is used to remove those points.
From the filtered data of the complete dataset, we can obtain the minimum and maximum for each objective, which allow scaling the whole dataset and getting the reference point for the hypervolume computation. The dataset for each magnetic flux density limit is then filtered in the same way, and the hypervolume indicator is calculated.
The calculation of the hypervolume indicator is given in the Appendix.

2) COMPARISON OF THE SPREAD OF OPTIMIZATION RESULTS
As in the evolution algorithm, the generations are calculated using a process that involves randomness, the whole process has been done 4 times to obtain average results. The computation results of the hypervolume indicator as a function of the magnetic flux density limits are shown in Fig. 13. From the figure, we see that the hypervolume indicator, and hence the spread of the points indicating the Pareto front, stays relatively constant, with a slight increase as a function of the magnetic flux density limit. There is no substantial influence of the magnetic flux density limit on the spread. Also, using FEM only, the hypervolume was calculated as 0.66. The optimization using FEM only does not give better results in terms of the spread of the results. Although not significantly better, SAFEMOM has a slight advantage if high magnetic flux density limits are used.

B. CONVERGENCE CALCULATION 1) CALCULATION PROCESS
To evaluate the convergence, we calculate the modified generational distance (GD+) to measure the distance between the optimal points of a generation and the approximate Pareto front. The evolution of this distance as a function of the generation gives an insight into the optimization speed.
The whole computation of the GD+ is shown in Fig. 14. For the computation, the complete dataset is used. Some pre-processing is applied to it. The NRBF and robust filter are used first. They give a dataset that has both a good distribution and can approximate the Pareto front. Using the GPR from the robust filter, we calculate uniform sampling points. According to the GPR, some of those points have low confidence. They are removed from the dataset using the robust filter. A reference setP * is obtained, shown in blue in Fig. 15.
The complete dataset is split into smaller datasets for the measurement of the convergence of each of them. Each of those smaller datasets is of one generation m only, and of one level n of magnetic flux density limit B limn only. Again, the robust filter is used to remove points of high variance. Different datasets S mn are obtained.

2) COMPARISON OF CONVERGENCE OF OPTIMIZATION RESULTS WITH DIFFERENT MAGNETIC FLUX DENSITY LIMITS
We use the GD+ to compare the convergence of the optimization results. The whole process is run four times to obtain the average results shown in Fig. 16. The results show explicitly that the magnetic flux density limit imposed on the subdomain preoptimization greatly impacts the optimization speed using FEM. Preoptimization with a very high magnetic flux density limit, ie 2.8 T, gives designs much closer to the optimum. Those designs converge much more quickly using FEM.
The GD+ shows the average distance between the points that have been found and their closest neighbor on the approximate Pareto front. From Fig. 16, we can see that SAFEMOM obtains excellent results in 30 generations only, if the imposed magnetic flux density limit is high. Designs that are preoptimized using a lower magnetic flux density limit converge slowly.
In Fig. 17, the convergence is compared between FEM-based optimization and the SAFEMOM with different magnetic flux density limits. We see that the SAFEMOM allows a better convergence than optimization based on FEM only. The SAFEMOM using a low magnetic flux density limit does not allow a much quicker convergence than the optimization using FEM only.

C. USEFUL FILTERS FOR THE MEASUREMENT OF THE QUALITY OF THE DATASET
Now that the measurement of the optimization results has been explained, we give here the detailed calculation of the filters used in this process.

1) NICHE-RADIUS-BASED FILTER (NRBF)
A genetic algorithm produces a dataset with many designs. Some designs are optimal, and others are not optimal. The NRBF filter ensures that only non-dominated designs are kept. Therefore, dominated designs, known to be suboptimal, are removed. As defined in the Appendix, dominated designs are designs with at least another design in the dataset that has all three parameters L, I , and T smaller. Evolution of GD+ value in SAFEMOM given as a function of the generations to show the convergence for different magnetic flux density limits of the subdomain model. The convergence using SAFEMOM is compared to the optimization using FEM only. The process has been done 4 times for each magnetic flux density limit to obtain an average value. The low limit is from 1.6 to 1.8 T, which is a typical value for the subdomain model (slightly higher than the knee point). The medium limit is from 2.0 to 2.4 T, the high limit is from 2.6 to 3.0 T.
The second function of the NRBF is to ensure that the distribution of the remaining designs is good. It does so by removing designs where the density of the non-dominated designs is too high. As a result, a better distribution is obtained.
For this purpose, the algorithm contains a loop that reduces the size of the dataset. At each iteration, it looks for the place with the highest density of non-dominated designs and removes a design from that place.
The definition of the NRBF is inspired by [42]. This filter reduces the data and ensures that the remaining data have a good distribution. As the NRBF uses the niche radius, we define it first. The niche radius represents the average distance between the points and the k th nearest neighbor of each of them. It is calculated in this way: wherep is any point representing the normalized objectives of a design in the dataset Q,p k is the k th nearest neighbor ofp, N Q is the number of points in Q, and d(p,p k ) is the distance betweenp andp k . In our study we use k = 6. Also, the crowding degree function is defined as: The computation process of the NRBF uses the notions above and is given in the form of a pseudo-code in Algorithm 1. Fig. 18 shows the impact of the NRBF. The original data are shown in blue, and the data kept by the NRBF is shown Algorithm 1 Niche-Radius-Based Filter (NRBF) Pseudo-Code Input: Q I : Input dataset N s : Number of points in the output dataset Output: Q O : Output dataset Process: k = 6 number of closest neighbors to consider Select the non-dominated points in Q I and store them in Q ND Remove duplicates from Q ND Calculate n niche,k using (6) for the dataset Q ND Calculate C cd (Q, r niche,k , p) ∀p ∈ Q ND using (7) and put in vector v c of length N d .
Update n niche,k using (6) for the dataset Q ND end Update C cd (Q, r niche,k , p) ∀p ∈ Q ND using (7) and put in vector v c of length N d .
in red. The figure shows that the NRBF keeps the key information from the Pareto front and ensures a relatively uniform distribution of points.

2) ROBUST FILTER
The second filter we define is called ''robust filter''. The robust filter aims to remove the points related to high uncertainty from the dataset, to increase our confidence in the dataset. For this purpose, the GPR is used with a Matern kernel with parameter 5/2 [47]. The algorithm is shown in Algorithm 2, as a pseudo-code.
The impact of the robust filter applied on the dataset already filtered by the NRBF is shown in Fig. 19. From the figure, we see that the robust filter removes some unreliable points which have low confidence, as calculated by the GPR. If those points had been taken into account, it would have led to a misleading measure of the hypervolume and GD+.

V. DISCUSSION
A comparison summary can be found in Table 4. In this Table, we see that the SAFEMOM with a high magnetic flux density limit has a significant advantage compared to optimization using FEM alone, but this is only the case when a very high magnetic flux density limit, ie 2.8 T, is imposed. In that case, it converges much more quickly than FEM. To reach a GD+ convergence of around 0.025, SAFEMOM needs only 10 generations, whereas optimization with FEM needs 50 generations. This means a reduction of about 80% of the time. If a slightly higher GD+ is accepted, and hence VOLUME 11, 2023    Calculation of the optimal designs used to calculate the Pareto surface: comparison between the three methods. The calculation time given is a calculation time that includes 50 generations and the preoptimization time for the SAFEMOM. ''Low magnetic flux density limit'' is from 1.6 to 1.8 T, which is a typical value for the subdomain model (similar to the knee point). ''High magnetic flux density limit'' is from 2.6 to 3.0 T.
From Fig. 17, we see that when we compare both SAFEMOM and FEM optimization method after four iterations only, the GD+ of SAFEMOM is reduced by a factor of nearly eight. This means that if only limited computation resources are allowed, then SAFEMOM gives much more optimal results. SAFEMOM using a 1.6 T magnetic flux density limit is not advantageous compared to FEM. Indeed, although 1.6 T is close to the knee point of the BH curve of the material, the optimization using this magnetic flux density limit gives designs with dimensions that are far from optimal. However,  Average torque as a function of tooth opening ratio, calculated with the subdomain model (red, yellow, and blue curves) using different magnetic flux density limits, and calculated through non-linear FEM (green curve). The calculation clearly shows the need to use high magnetic flux density limits to obtain designs structurally closer to the optimum calculated using a non-linear method.
an optimization using a high magnetic flux density limit gives designs that have dimensions much closer to the dimensions found on designs on the Pareto surface.
This may seem counterintuitive. Indeed, the torque produced by design with 3 T magnetic flux density limit has more error than that with 1.6 T magnetic flux density limit, as shown in Fig. 20. Moreover, this graph is consistent with the material properties shown in Figs 3 and 4. The error is calculated using non-linear FEM as a reference.
The fact that the subdomain optimization using a high magnetic flux density limit gives machines parameters much closer to the optimum can be explained using both FEM and subdomain model, as seen in Fig. 21. In that calculation, we calculated the torque with FEM and with the subdomain model for different width of teeth, with different magnetic flux density limits. We see clearly on the graph that if the magnetic flux density limit is 2.8 T, although the torque is very different between FEM and subdomain, the optimal width of teeth is similar. Moreover, if the magnetic flux density limit is too high, namely 3.6 T, this allows a more diverse set of parameters in the preoptimization. Therefore, the optimization using a genetic algorithm reaches the optimum more quickly.

VI. CONCLUSION
It is believed that to increase the accuracy of linear subdomain models when used to model PM motors having materials that saturate, we need to limit the field in the ferromagnetic part to stay under the level of magnetic saturation of those materials. This belief seems intuitive. When used in an optimization process, this approach using linear subdomain models produces designs that exhibit torques similar to those of designs optimized using non-linear FEM. However, as shown in this study, the dimensions of the designs optimized with a linear subdomain model, using the actual saturation level of the material for the magnetic flux density limit, are useless. Indeed, non-linear optimization shows that those designs are far from the optimum. It is difficult for the optimization to converge to optimal dimensions starting with the dimensions of those designs. Therefore, this common belief is inefficient for optimization.
The solution to this problem is that the magnetic flux density limit in the linear subdomain model needs to be taken as very high, in our case 2.6 T to 3.0 T, to ensure that the dimension design space is not over-constrained. This high limit leads to designs that exhibit a different torque than their counterpart calculated with non-linear FEM, but those designs have dimensions much closer to the optimum.
When this technique is used as a quick preoptimization, it allows a much quicker and more precise optimization afterward. Indeed, this has been quantitatively assessed with the GD+ measure. Also, it allows the hypervolume to be increased slightly with a high magnetic flux density limit, which means that the space covered is slightly bigger. Therefore, increasing the magnetic flux density limit was not done at the cost of a cover range. SAFEMOM has a similar cover range compared to the optimization using FEM only. SAFEMOM can reduce about 80% of optimization computing time compared to FEM, for a similar GD+ convergence of about 0.025 in this optimization process. If a higher GD+ and hence higher error are accepted, the computing time is even shorter. Also, if only limited computation resources time is allowed, then the GD+ value is divided significantly (about eight times). SAFEMOM hence produces less error when the computation resources are limited.
The optimization results of this study have been obtained using a broad family of SMPM machines, not only a few designs. Also, the steel used in this study is typically used in industry. It exhibits a typical non-linear relationship between magnetic flux density and magnetic field intensity. Therefore, the validity of the conclusion is broad.
As the measurement tools presented in this research are general, and as the subdomain models can be quickly developed for other types of structure, the presented method has broad application for other structure types and applications.
We expect a further decrease in the computation burden to be achieved through more complex constraint functions for the magnetic flux density in the subdomain. Future work will involve using these techniques to optimize motion control systems.

APPENDIX. HYPERVOLUME INDICATOR CALCULATION
The calculation of the hypervolume indicator is based on the definition given by Zitler and Künzli [45]. The definition is reduced to a hypervolume of 3 dimensions.
Each design is identified by its 3 objectives p = (L, I , T ), which are the losses, inertia, and negative peak torque. Taking the maximum and minimum values of (L, I , T ) for our dataset, the losses, inertia, and torque are normalized from 0 to 1. The normalized designs are hence given byp = (L,Ĩ ,T ) ∈ (0, 1) 3 . By definition, the weak Pareto-dominance relation ⪰ between two designsp 1 andp 2 indicates that the designp 1 is at least as good asp 2 , which we write (p 1 ⪰p 2 ), if and only ifL 1 ≤L 2 ,Ĩ 1 ≤Ĩ 2 , andT 1 ≤T 2 .
In the same way, we can write this definition for sets. Let us define two sets of designsP 1 andP 2 , then, by definition, P 1 weakly dominatesP 2 if and only if ∀p 2 ∈P 2 , ∃p 1 ∈P 1 : p 1 ⪰p 2 [45].
This is calculated through the code given in PlatEMO [49]. Currently, he is a Professor with the College of Electrical Engineering, Zhejiang University, Hangzhou, China. His research interests include the application, control, and design of electrical machines. Prof. Fang achieved the Second Prize of National Scientific and Technological Progress Award twice, the First Prize of Provincial Science and Technology Progress Award twice, and the First Prize of Provincial Science and Technology Progress Award three times. VOLUME 11, 2023