Electric Field Evaluation Using the Finite Element Method and Proxy Models for the Design of Stator Slots in a Permanent Magnet Synchronous Motor

: The efﬁciency of electric motors is being improved every day and projects with design variations can improve their performance. Among electric motors, the Permanent Magnet Synchronous Machine (PMSM) is being increasingly used, because of its growing use in electric vehicles. Simulating design variations using the Finite Element Method (FEM) can improve PMSM design, and by optimizing the parameters based on the FEM, even better results can be achieved. The design of the PMSM stator slots must be evaluated, as conductors are accommodated and an electrical potential is applied at this location. The FEM parameters are varied, and the results can be used to build an approximate model, known as a proxy model. The proxy model can then be used in a mathematical programming problem to optimize the design of stators that have less electric ﬁeld in certain regions, thus reducing the chance of developing a failure. The results of the proposed methodology show that its application is promising for machine design and can also be used for the design of other systems.


Introduction
Efficiency in electrical equipment has been improved over the years, and materials and methods that reduce energy consumption are increasingly being studied [1][2][3]. Currently, designers need to develop robust products for greater durability and more efficiency [4]. To improve the performance evaluation of projects before developing prototypes, it is possible to use software that simulates the conditions that the equipment will be exposed to, such as mechanical or electrical stress [5]. The Finite Element Method (FEM) is being increasingly explored to simulate adverse conditions in electrical equipment. For electric field simulations using FEM, COMSOL Multiphysics software is a complete platform for simulating variations in the machine design [6].
The application of FEM to assess the condition of electric motors has gained notoriety in scientific research. Fonteyn et al. [7] developed a study on the effect of magnetostriction and electromagnetic stress on iron in induction machines through FEM. Using the same method, Parreira et al. [8] published an evaluation of a switched reluctance machine. The Permanent Magnet Synchronous Machine (PMSM) is widely used, as it has an application for electric motorization of hybrid urban vehicles [9]. Meessen et al. [10] successfully used the FEM to calculate the inductance through the flow on the d and q axes in a PMSM motor. However, to have an effective machine design, in addition to applying the FEM, it is necessary to perform the optimization of the parameters.
In order to solve the nonlinear analysis of the optimization parameters, Luukko and Pyrhonen [11] showed that the definition of the problem constraints must be considered. One of the restrictions in the analysis performed is the constant flow control. According to the results presented, the permanent magnet flux connection parameters and PMSM axis inductances can be calculated. From these parameters, a procedure was performed to generate a proposal for the design of a more efficient machine. To optimize the PMSM efficiency using the FEM, the changeable flow's weakening angle can be used in the armature current control strategy; thus, the ideal angle can be obtained [12]. Another application of the FEM using optimization for the design of the PMSM is presented by Jeong, Kim, and Hur [13], in a hybrid-type permanent magnet. In this work, an equivalent circuit was evaluated in relation to the motor's magnetic characteristics to improve the motor's performance and reliability.
The FEM can be applied in other ways to improve PMSM performance. The method can be used to demonstrate that the resulting harmonic content in the phase currents, when operating at speeds above the bandwidth of a controller, result in increased magnet losses [14]. In addition, evaluating the harmonic content, Zhang et al. [15] proposed a machine with windings connected circularly, making it suitable for high-power and high-performance marine propulsion applications. In this research, the best design was also validated by simulation using the FEM. Prototypes can be simulated with the FEM before they are manufactured, and thus, the built model has a better performance [16]. Even the use of the FEM for PMSM design can be extended to robotic applications to determine the design through an analytical and experimental evaluation [17].
A comparison between several PMSM rotor topologies for electric vehicles in relation to geometric parameters was presented by Liu et al. [18]. The paper showed that saturation effects, flow on the d and q axes, and other factors must be considered in analyses using the FEM to obtain accurate evaluated motor parameters. FEM analysis can also be used to improve the assessment of interference to reduce the breakdown of rotary axes in PMSM, by evaluating the stress to which the engine will be subjected through simulation [19]. In [20], the variation in the slots' design through an optimization process was also presented using the FEM for the design of a motorcycle engine equivalent.
The use of nonexplicit functions in optimization (i.e., only available through simulators) is given by means of "simulated constraints" that model the system behavior. Among the approaches found in the literature to deal with simulated constraints, two are mentioned here. In the first approach, derivative-free optimization (DFO) algorithms are used to communicate directly with the simulator, as in [21], where the authors used Mesh Adaptive Direct Search (MADS) to tune oil well parameters. Although this approach can generate good results, it creates a dependency on the simulator and results in an increased computational time.
Another approach is presented in [22][23][24], in which Particle Swarm Optimization (PSO) was applied. The swarm topology defines which particles a subset can exchange information with, allowing all particles to communicate, and so, the entire examination tends to converge to the same position [25]. The PSO can be applied to reduce vibration and noise and increase the operational efficiency of the motor [26]. In addition, the use of this algorithm can be applied to maximize the torque output in which the method is applied to optimize the dimensions of a coaxial magnetic gear [27].
In the second approach, proxy models can be used to create approximations of the simulators' models, as in [28], where the authors used polynomial proxy models to approximate pressure drops in pipelines. In this approach, the simulator is only used during the data collection and proxy model generation phase. The optimization process can then be performed based on the proxy model, with a reduced computational cost. The second approach is considered in this paper.
There is also a line of research that makes use of the Taguchi method [29] in combination with the FEM to optimize parameters. For instance, [30] employed the Taguchi method and finite element analysis in the optimization of automobile roof technical parameters. The authors in [31] used the Taguchi Method combined with the FEM to analyze a robust new design for titanium alloy prick hole extrusion.
On the basis of the demands of highly efficient machine designs, this study aimed to evaluate the influence of variations in the PMSM stator design. For this, we employed the approach to create a proxy model to approximate an optimized stator machine design in order to integrate the simulated constraint to state-of-the-art solvers to derive a stator slot design. For this evaluation, possible design variations and their influence on the electric field distribution were compared.
The greater intensity of the electric field in specific regions of a piece of electrical equipment can generate partial discharges (PDs). PDs are, in general, the result of electrical stresses in the insulation or surroundings and generally appear as short pulses. This phenomenon, when continuous, may give rise to small disruptive discharges, which over time may generate failures in electrical equipment [5].
Because the stator grooves accommodate the turns, there is a greater intensity of electric and magnetic fields in this region. With the focus on decreasing the intensity of the electric field applied in a specific location, calculating the intensity of the electric field can improve the field distribution capacity and thereby reduce the chance of failure in the PMSM over time.
With the methodology presented in this paper, after creating the model, the designer will no longer need to use the simulator to run different scenarios (i.e., changing the restrictions involved in the limits of the variables of interest). In addition, after the development of a proxy model that fits the results obtained by the FEM, the model can be integrated into different state-of-the-art optimization algorithms, enabling the search for the maximum or minimum global optimum involved in the design.
The scope of this research is limited to presenting a simplified methodology to optimize the parameters of electric fields based on the FEM. For this purpose, an illustrative example is used based on a small dataset, in order to be didactic and easily reproducible. It is understood that what is presented here can serve as a guide for the development of other works. This paper is organized as follows: In Section 2, the necessary calculations and parameters used for the design of a stator for an electric motor of the PMSM type are presented. Section 3 discusses the proposed method. In Section 4, the variations in the machine design are evaluated, and then the design and results are discussed. At the end of the paper, a conclusion is presented based on the analysis performed.

Machine Geometry
The PMSM can be analyzed in parts, i.e., divided into regions as can be seen in Figure 1, wherein d sb denotes the stator depth, s st denotes the depth of the stator teeth, g denotes the gap, d m denotes the depth of the permanent magnet, d rb denotes the depth of the rotor, d i denotes the depth of the inert region, and r rs denotes the radius of the rotor axis. All of these variables, except for the rotor radius, can be adjusted to change the machine's design and improve its efficiency. From the values of these variables, it is possible to calculate the r ri radius of the inert region of the rotor, the r rb radius of the bottom of the rotor, the r rg radius of the rotor air gap, and the r st radius of the stator tooth, according to the equations: The slot dimensions are as follows: the θ tt angle to the tip of the tooth in relation to the radius r st ; the θ st angle to the slot in relation to the radius r st ; the r si internal tooth radius; the θ ti angle to the radius r si ; the θ tb angle to the radius r sb ; the w tb width of the base of the tooth; the d tb depth of the base of the tooth; the d tte depth to the edge of the tooth; and the d ttc depth of the center of the tip of the teeth. On the basis of the ratio of the number of slots S s and the angles of the motor dimensions, the r sb radius of the inert region of the stator and the r ss radius of the stator housing were obtained through the following equations: In the stator, the flow enters and leaves the slots in a tangential direction, and in the stator slots, the conductors are accommodated where the electric potential is applied [32]. Figure 2 shows a cut rectangular approximation of the slot, in which the teeth of the stator slot can also be seen. For the dispersion inductance calculation, the slot geometry is approximate to a rectangle and the tooth tip's width is approximated as the circumferential length [33].
The tooth tip's rectangular depth is adjusted so that it has the same area as the cross-section.
The width between the tips of the teeth is approximate to the circumferential distance between the teeth.
The width of the base of the stator slots, where the windings are accommodated, is the average of the distance of the rope length from the inner corners.
Finally, maintaining the slot area and the base of the teeth, the depth and width of the teeth must obey the following equations: For analysis purposes, a complete filling of the grooves with the windings was considered, i.e., d wR = d siR . From the Pareto-optimal design [34], Table 1 shows the dimensions of the PMSM to be analyzed.

Optimized Finite Element Method
In this section, the proposed method is presented. In this study, the FEM was used to evaluate the electrical potential distribution in the PMSM stator slots. From the results of the variation of the parameters, the analyses were performed, the values generated the cost function to be optimized, and the design was obtained.
For the FEM analysis, an object is divided into equal parts, and based on the equations, each element's behavior is estimated [35]. Through this method, all individual behaviors are analyzed and added, thus obtaining the object's total behavior resulting from the sum of each element [36]. With the sum of the elements, a mesh is obtained, which can have variations that result in greater precision and computational cost, depending on the configuration used.

FEM Elements Calculation
The analysis of the finite element method is divided into steps, which describe how the analysis is produced, so that it is possible to achieve the expected results. These steps are from 1 to 4 and are classified in order of application. There are several ways to calculate each element, which are subdivided for mathematical evaluation [33]. In this paper, we used triangular shapes, which are explained in this subsection.
In the first stage, the analyzed structure should be subdivided into triangular shapes, with the objective of approximating the potential V e inside a structure. A general solution for the contour structure is given by where N is the number of elements. The second step is to analyze the triangular figure which is composed of three ends, which in turn have powers called V e1 , V e2 , and V e3 and are presented in following multiplication: By solving this multiplication, it is possible to determine a, b, and c. When replaced in the V e equation, such a resolution generates a new equation responsible for demonstrating the potential at any (x, y) points within the structure, as long as its potentials are known. with where λ is the area of the structure and will be positive if the intersections are numbered clockwise. It is obtained by the following equation: in which ∝ i is defined as the shape function of the elements and obeys the following properties: After defining and calculating each element, the algorithm performs a refinement process from an initial mesh, evaluating finer meshes [23]. After evaluating the meshes, an asymptotic behavior is found, and in this way, the variations become smaller until the algorithm converges [37]. The elements are connected through nodal points, and these form a mesh [38], as shown in Figure 3 for analysis with triangular elements. The evaluation through the finer meshes achieves the necessary accuracy of the electric field for this model, considering that a more detailed evaluation results in a greater computational effort, without significant improvement in the results. The design of electromagnetic devices for PMSM based on the FEM can be accomplished by evaluating significant parameters using the gradient of difference between the analytical model's outputs and those obtained with the FEM. The parameters are used to build a response in a small space searching for an optimal project [39].

Parameter Optimization
To start the design evaluation process, first a range of data is simulated individually and offline via the FEM. Notice that the values to be considered must be within the operator's range of interest. Subsequently, a proxy model is built from the simulations performed. The proxy model can be derived from several mathematical functions, such as linear regressions, polynomial regressions, and neural networks. These mathematical functions have the objective of minimizing the difference between the output of the model and the simulated value (the given input data): Equation (26) is known as the loss function, in which e(·) is the error, given an input parameter y and an estimated valueŷ. The choice of the mathematical function to be used in the representation of the model must be based on the type of process being modeled, which requires the operator's expertise.
The linear and polynomial regression can be easily obtained by solving a least-squares problem; however, other regression methods can be used to adjust the curve, such as lasso, ridge, and fuzzy regression [40,41]. While lasso and ridge regression can be used to better select the features that represent the model, fuzzy regression can be used to consider the uncertainties of the coefficients in a linear regression (i.e., when the relationship between model parameters are vague) [42].
Many models based on artificial intelligence can be used for regression, these models have a wide range of applications in forecasting time series. Combined techniques can further improve the regression model, as shown in [43], wherein the group method of data handling is used with the Wavelet transform. Regarding a deep learning approach, long short-term memory is used in [44] with the same purpose, in which a neuro-fuzzy model [45] evaluates a problem using Wavelet transform for filtering the signal.
For the validation of the proxy model, the parameter R 2 was used. The R 2 is a measure of the system response variation that is explained by the model and is given by Consider y to be the mean of the training data (simulation data) as follows: then, the sum of the squares total (SST) is the sum of squared residuals (SSR) is and the coefficient of determination is given by Furthermore, the mean squared error (MSE) is used to measure the error of the proxy model in predicting the simulator data. The MSE can be calculated as follows: After validating the proxy model by means of R 2 and MSE, it can be used to substitute the simulated constraint in the problem of interest: where H(·) is now given by the approximated function. Typically, the problem given by Equation (33) is a Nonlinear Problem (NLP), which is well studied in the literature. NLP is a class of optimization problems with continuous variables, with nonlinear functions in the objective and/or the constraints. They can be solved efficiently using state-of-the-art solvers and algorithms, such as the interior-point method available in the Interior Point OPTimizer (IPOPT) solver [46] or the Mesh Adaptive Direct Search (MADS) available in NOMAD [47]. Additionally, a bounding constraint can be added to Equation (33) in order to limit the values that x can assume, as follows: It is important to notice that if the approximate function is not convex, algorithms such as the interior-point method can converge to a local optimum. If a global optimum is desired, an alternative is to use algorithms like spatial branch-and-bound combined with convex underestimators such as McCormick envelopes [48]; however, they can incur high computational costs. Thus, the overall design method presented in this paper can be summarized as presented in Figure 4. In this paper, the used methodology is based on a quantitative approach [49].

Simulation Setup:
Perform FEM analysis on the object of interest.

Proxy Model:
Fit proxy model from simulator data. Validate the model using (31) and (32).

Parameters Optimization:
Find the optimal design parameters by means of (33).

Analysis of the Proposed Methodology
In this section, the analysis of the proposed model is presented. Initially, the FEM application results are discussed, and then the optimizer is presented, discussed, and evaluated. The evaluated locations in this paper are shown in Figure 5. The number of slots and the slots' dimensions were not changed, as these variations would represent changes in the machine's power. In "A", the changes were made in relation to the variation of the dimensions of W stR , which was explained in the previous section. The electric field distribution results were equivalent in all the variations calculated for this parameter, which shows that this variation was not effective, as it reduces the area of use of the stator slot. The distribution of the electric field can be seen in Figure 6. Another parameter that did not result in variation in the distribution of electrical potential was the change in W ttR and d ttR . These results were expected since the electrical potential was applied to the conductors in the stator slots, and the ground was considered in the stator housing.
As the project aimed to improve the motor's performance in relation to the electric field distribution without reducing its electrical capacity, it is not possible to change the stator tooth area, given by W tbR and d siR . Consequently, it is not possible to change W siR and d wR , as these parameters are the area for filling the slot of the conductors in the stator.
The greatest variation in the distribution of electrical potential occurred in "B". This location is the region that has the shortest distance between the applied potential and the ground. This variation in electrical potential can be seen very clearly in Figure 7, which shows a contour distribution. The values of the "B" analysis, for the variation of the radius of the stator corner, are shown in Table 2. Just increasing or decreasing the dimensions of this parameter does not represent a simple answer for a best design, as the response of this variation is not linear. Therefore, it is convenient to optimize the parameters. The justification for "B"'s evaluation is also supported by the field lines' direction and concentration, as shown in Figure 8. Thus, the variation of this parameter is the focus of the application of the optimized finite element method that is presented in the next subsection.

Optimized Finite Element Method
Several methods can be used to create the proxy model with respect to Table 2. In this section, we present two examples, the first using polynomial regression, and the second using a neural network.

Polynomial Regression
Given the parameters presented in Table 2, we started to build the proxy model based on a linear regression and then escalated it to a polynomial regression of degree p. Consider that a n-dimensional polynomial regression of degree p is given bŷ First, we present the tests for polynomial regression with orders varying from 1 to 3, as shown in Figure 9. It can be seen that the regression models do not fit well to the data given by the simulator, i.e., with the R 2 and MSE values presented in Table 3.

5DGLXV>PP@ 9DOXH>9@
(a) Proxy model with first order polynomial.   Table 2 (red dots). Given the nonlinearity of the data presented, it can be seen that although the R 2 increases along with the polynomial degree, the order tends to grow before reaching an adequate fit, resulting in a very complex model. An alternative in the literature for nonlinear data is polynomial transformation, which represents nonlinear relationships between the dependent and the independent variable [50]. Starting from the third-order polynomial previously tested, consider the following regression with the addition of a polynomial transformation: From the results found for the first-, second-, and third-order proxy models, it is not possible to find a regression that defines the optimal value for the evaluated variation. As can be seen in Figure 9, the regression of these orders does not result in an approximate model of the function to be optimized. With the addition of the polynomial regression, the new fitted curve is presented in Figure 10. Notice that polynomial transformation is better suited to the data extracted from the simulator, with an R 2 of 0.99 and MSE of 0.79. 5DGLXV>PP@ 9DOXH>9@ Figure 10. Proxy model with logarithmic transformation (blue line) applied to the data presented in Table 2 (red dots).
The derived regression with the logarithmic transformation was given by which can be easily solved by IPOPT [46], resulting in an optimized radius value of 0.9 mm.

Neural Network
Another possibility for the creation of the proxy model is the use of a neural network. To employ the trained neural network into an optimization problem compatible with a mathematical programming solver such as IPOPT [46] and Gurobi [51], one must carefully choose its activation function. For example, if the ReLU (Rectified Linear Unit) activation function is chosen, techniques such as those presented in [52] can be used to transform the problem. The authors in [52] proposed rewriting the ReLU activation function as a mixed-integer linear problem (MILP) by means of binary variables.
For instance, consider a neural network with K + 1 layers, with 0 being the input layer and K the output layer. Each hidden layer k ∈ {1, . . . , K − 1} has an output vector x k , which is calculated as in which σ(·) is the activation function, W k ∈ R n k ×n k −1 is the matrix containing the weights of the neural network for each layer, along with its biases b k ∈ R n k . Thus, x k ∈ R n k is the output of layer k. When the ReLU function is defined as the activation function, σ(·) can be defined as Following the methodology of [52], we can rewrite the ReLU operator as Notice that following Equation (41), the ReLU is decoupled into a positive part (x ≥ 0) and a negative part (s ≥ 0). To obtain this behavior, we must ensure that at least one of the two terms is zero. This can be achieved by means of big-M constraints, assuming finite values of L and U so that In this sense, using a big-M constraint, with the aid of a binary variable z, it is possible to ascertain the ReLU logic in mixed-integer programming form: Finally, we can rewrite the entire ReLU neural network as a MILP, in which the input layer is given by the hidden layers (∀k = 1, ..., K − 1) by and the output layer (K) by Thus, given a neural network with the ReLU activation function, the Equations (44)-(46) represent its exact formulation as a mixed-integer linear problem. That is, given previous training, considering the weights to be constant, for the same fixed input x 0 , the output will be constant.
Using the data provided by Table 2 to train a ReLU neural network and recasting them to a MILP as suggested by [52], we are also able to find the optimized radius for the illustrative problem. MILPs can be solved to a global optimum by means of efficient algorithms such as branch-and-bound, which is available in solvers like Gurobi [51]. Considering a ReLU neural network with five hidden layers (K = 5) with 10 neurons per layer, we achieve an optimized radius value of 1.1 mm and an MSE of 5.56 × 10 −4 . The resulting neural network is shown in Figure 11. 5DGLXV>PP@ 9DOXH>9@ Figure 11. Proxy model with ReLU neural network (blue line) applied to the data presented in Table 2 (red dots).

Discussion
On the basis of the results presented to optimize the radius of the stator slot groove, a lower intensity electric field applied to a specific point of the machine is obtained. As the largest electric field can generate partial discharges, and thus failures develop in these locations, the lower the electric field strength, the better the reliability of the machine over time.
Although the illustrative example given in this section is limited to a single dimension problem (i.e., the decision variable is the radius), the design methodology using FEM parameters can be expanded to n-dimensional problems. Creating a proxy model to optimize design parameters increases flexibility and speed for the designer, who is able to search for ideal parameters without accessing the simulator. Thus, it is understood that the methodology proposed here can be used to combine the best of both worlds in the design of stator slots in the PMSM, making use of the FEM to perform precise simulations, combined with exact optimization techniques capable of finding the best parameters for the project.

Conclusions
The FEM application proved to be a promising technique for assessing the electrical potential distribution in the PMSM. On the basis of an evaluation of the applied potential, it is possible to check the locations that are most likely to develop a failure. In addition, the use of a computational model to optimize the parameters of this type of engine is promising, according to the results discussed in the previous section.
The optimization procedure presented in this paper involved three steps: (i) the simulation setup; (ii) the proxy model; and (iii) parameter optimization. In this sense, the PMSM was simulated through the FEM, and its data were collected and a proxy model was built from this.
For the illustrative example of this paper, functions with a low degree of polynomial order did not generate promising results in relation to the accuracy, as represented by the coefficient of determination (R 2 ). To obtain project results that are promising for real applications, it was necessary to use functions that are representative of the project's variations. Other structures may require different regression techniques, depending on their complexity and nonlinearity, such as the example of neural networks, which, when rewritten as mixed-integer linear problems, can be solved to their global optimum.
The technique applied in this paper can be used to improve the design of several electrical system components that are influenced by electrical potential. An optimized profile of electrical components can reduce their vulnerability, improving their robustness to the conditions to which they will be exposed.
The great benefit of applying the method presented in this paper is that, from the variation of the parameters, it is possible to identify an optimized design model, and this assessment can be performed for different components of the electrical power system, in addition to electrical machines. Specifically, it was observed that the variation in the corners of the stator slot can generate a greater intensity electric field in certain places, which can be more susceptible to failures.
In future work, it is understood that the methodology applied here can be applied to other parameters, such as torque, loss, and efficiency. Furthermore, other structures can be analyzed. The assessment of the proxy model's behavior for more decision variables can be investigated, based on the same methodology.

Acknowledgments:
We would like to thank to the Coordination of Superior Level Staff Improvement (CAPES), awarding a doctoral scholarship to one of the authors. The authors also like to acknowledge the collaboration of Computer lab 7 of Instituto de Telecomunicações-IT Branch Covilhã.

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

Abbreviations
The following abbreviations are used in this manuscript: