electronics

Article
Electric Field Evaluation Using the Finite Element
Method and Proxy Models for the Design of Stator
Slots in a Permanent Magnet Synchronous Motor
Stéfano Frizzo Stefenon 1,2, * , Laio Oriel Seman 3 , Clodoaldo Schutel Furtado Neto 1 ,
Ademir Nied 1 , Darlan Mateus Seganfredo 2 , Felipe Garcia da Luz 2 , Pablo Henrique Sabino 2 ,
José Torreblanca González 4 and Valderi Reis Quietinho Leithardt 5,6,7
1 Electrical Engineering Graduate Program, Electrical Engineering Department, Santa Catarina State
University (UDESC), R. Paulo Malschitzki 200 (North Industrial Zone), Joinville 89219-710, Brazil;
[email protected] (C.S.F.N.); [email protected] (A.N.)
2 Electrical Engineering Department, University of Planalto Catarinense (UNIPLAC),
Av. Mal. Castelo Branco 170 (Universitário), Lages 88509-900, Brazil;
[email protected] (D.M.S.); [email protected] (F.G.d.L.);
[email protected] (P.H.S.)
3 Graduate Program in Applied Computer Science, University of Vale do Itajaí (UNIVALI),
R. Uruguai 458 (Centro), Itajaí 88302-202, Brazil; [email protected]
4 Expert Systems and Applications Lab, Faculty of Science, University of Salamanca, Plaza de los Caídos s/n,
37008 Salamanca, Spain; [email protected]
5 COPELABS, Universidade Lusófona de Humanidades e Tecnologias, 1749-024 Lisboa, Portugal;
[email protected]
6 Departamento de Informática da Universidade da Beira Interior, 6200-001 Covilhã, Portugal
7 VALORIZA, Research Center for Endogenous Resources Valorization, Instituto Politécnico de Portalegre,
7300-555 Portalegre, Portugal
* Correspondence: [email protected]




Received: 30 September 2020; Accepted: 19 November 2020; Published: 22 November 2020


Abstract: The efficiency 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
field 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.

Keywords: finite element method; evaluation of design; proxy model; polynomial regression

1. Introduction
Efficiency in electrical equipment has been improved over the years, and materials and methods
that reduce energy consumption are increasingly being studied [1–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

Electronics 2020, 9, 1975; doi:10.3390/electronics9111975 www.mdpi.com/journal/electronics

Electronics 2020, 9, 1975 2 of 20

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–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
Electronics 2020, 9, 1975 3 of 20

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.

2. Machine Geometry
The PMSM can be analyzed in parts, i.e., divided into regions as can be seen in Figure 1,
wherein dsb denotes the stator depth, sst denotes the depth of the stator teeth, g denotes the gap,
dm denotes the depth of the permanent magnet, drb denotes the depth of the rotor, di denotes the depth
of the inert region, and rrs 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.
Electronics 2020, 9, 1975 4 of 20

r rg
r sb r st
r rb

r ri
r ss Slot
r rs
di Shaft g Tooth
Rotor
d rb
dm
Magnetically
d st Inert Region
Air
G ma
nen
Per agnet
t
d sb
ap M

Stator

Figure 1. Cross-section of the Permanent Magnet Synchronous Machine (PMSM).

From the values of these variables, it is possible to calculate the rri radius of the inert region of the
rotor, the rrb radius of the bottom of the rotor, the rrg radius of the rotor air gap, and the rst radius of
the stator tooth, according to the equations:

rri = rrs + di (1)

rrb = rri + drb (2)
rrg = rrb + dm (3)
rst = rrg + g. (4)

The slot dimensions are as follows: the θtt angle to the tip of the tooth in relation to the radius
rst ; the θst angle to the slot in relation to the radius rst ; the rsi internal tooth radius; the θti angle to the
radius rsi ; the θtb angle to the radius rsb ; the wtb width of the base of the tooth; the dtb depth of the base
of the tooth; the dtte depth to the edge of the tooth; and the dttc depth of the center of the tip of the
teeth. On the basis of the ratio of the number of slots Ss and the angles of the motor dimensions, the rsb
radius of the inert region of the stator and the rss radius of the stator housing were obtained through
the following equations:

2 · π · αt
θt = (5)
Ss
2 · π · αtt
θtt = (6)
Ss
2·π
θst = − θtt (7)
Ss
 
θt
wtb = 2 · r st · sin (8)
2
 
θtt
wtt = 2 · r st · sin (9)
2
s
 w 2   2
tb θt
rsb = + (rst · cos + dtb + dttc ) (10)
2 2
Electronics 2020, 9, 1975 5 of 20

 
wtb
θtb = 2 · a · sin (11)
2 · rsb
dst = rsb − rst (12)
rss = rsb + dsb . (13)

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.

w stR
w ttR d ttR

d wR
d siR w siR w tbR

Figure 2. Rectangular machine stator slot approximation.

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].

wttR = rst · θtt . (14)

The tooth tip’s rectangular depth is adjusted so that it has the same area as the cross-section.
att
dttR = . (15)
wttR

The width between the tips of the teeth is approximate to the circumferential distance between
the teeth.
wstR = rst · θst . (16)

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.

   
π θ π θ
wsiR = rsi · sin − ti + rsb · sin − tb . (17)
Ss 2 Ss 2

Finally, maintaining the slot area and the base of the teeth, the depth and width of the teeth must
obey the following equations:
a
dsiR = slt , (18)
wsiR
a
wtbR = tb . (19)
dsiR
For analysis purposes, a complete filling of the grooves with the windings was considered,
i.e., dwR = dsiR . From the Pareto-optimal design [34], Table 1 shows the dimensions of the PMSM to
be analyzed.
Electronics 2020, 9, 1975 6 of 20

Table 1. Geometric characteristics of the machine [34].

Parameter Unit Value

External diameter (2· rss ) cm 16.6
Stator slot depth (dsiR ) mm 23.8
Rotor radius (rrs ) mm 34.8
Gap (g) mm 2.28
Stator depth (dsb ) mm 13.2
Rotor depth (drb ) mm 14.8
Number of poles (P) - 4
Number of stator slots (Ss ) - 24

3. 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.

3.1. 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 Ve inside a structure. A general solution for the contour
structure is given by
N
V ( x, y) ≈ ∑ a + bx + cy, (20)
e =1

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 Ve1 , Ve2 , and Ve3 and are presented in following multiplication:
     
Ve1 1 x1 y1 a
V =
 e2  1 x2 y2  ·  b  . (21)
     
Ve3 1 x3 y3 c

By solving this multiplication, it is possible to determine a, b, and c. When replaced in the Ve

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.

3
Ve = ∑ ∝i (x, y)Vei (22)
i =1
Electronics 2020, 9, 1975 7 of 20

with

∝1 = 1
2λ [( x2 y3 − x3 y2 ) + (y2 − y3 ) x + ( x3 − x2 ) y]
∝2 = 1
2λ [( x3 y1 − x1 y3 ) + (y3 − y1 ) x + ( x1 − x3 ) y] (23)
∝3 = 1
2λ [( x1 y2 − x2 y1 ) + (y1 − y2 ) x + ( x2 − x1 ) y]

where λ is the area of the structure and will be positive if the intersections are numbered clockwise.
It is obtained by the following equation:

1
λ= [( x2 − x1 ) (y3 − y1 ) − ( x3 − x1 ) (y2 − y1 )] (24)
2
in which ∝i is defined as the shape function of the elements and obeys the following properties:
(
1, i=j
∝i ( x, y) = (25a)
0, i 6= j
3
∑ ∝i (x, y) = 1. (25b)
i =1

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.

Figure 3. Normal mesh for finite element analysis.

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].
Electronics 2020, 9, 1975 8 of 20

3.2. 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):

e(k ) = y(k ) − ŷ(k). (26)

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 R2 was used. The R2 is a measure of the
system response variation that is explained by the model and is given by

Explained Variation
R2 = . (27)
Total Variation
Consider y to be the mean of the training data (simulation data) as follows:

K
1
y=
K ∑ y ( k ), (28)
k =1

then, the sum of the squares total (SST) is

N
SST = ∑ ( y ( k ) − y )2 , (29)
i =1

the sum of squared residuals (SSR) is

N N
SSR = ∑ (y(k) − ŷ(k))2 = ∑ e(k)2 , (30)
i =1 i =1

and the coefficient of determination is given by

SST
R2 = 1 − . (31)
SSR
Electronics 2020, 9, 1975 9 of 20

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:

K
1
MSE =
K ∑ (y(k) − ŷ(k)2 ). (32)
k =1

After validating the proxy model by means of R2 and MSE, it can be used to substitute the
simulated constraint in the problem of interest:

min f ( x ) =
x
∑ f i ( xi ) (33a)
i ∈I
s.t. H ( x ) = 0 (33b)
x∈X (33c)

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:

xmin ≤ x ≤ xmax . (34)

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).

Figure 4. Steps for the synthesis of a scenario.

4. 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.
Electronics 2020, 9, 1975 10 of 20

Figure 5. Locations evaluated in the stator machine.

In “A”, the changes were made in relation to the variation of the dimensions of WstR , 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.

Figure 6. Distribution of the electrical field.

Another parameter that did not result in variation in the distribution of electrical potential was
the change in WttR and dttR . 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 WtbR and dsiR . Consequently, it is not possible to change WsiR and dwR , 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.
Electronics 2020, 9, 1975 11 of 20

Figure 7. Contour electrical potential 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.

Table 2. Electrical voltage in corner of the stator.

Corner Radius of the Stator (mm) Value (V)

00 320.1
01 268.2
02 270.7
03 279.7
04 282.5
05 285.9
06 287.4
07 289.5
08 291.3
09 292.4
10 291.7

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.

4.1. 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.
Electronics 2020, 9, 1975 12 of 20

Figure 8. Electric field lines in the stator machine.

4.1.1. 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 by
n n n n n
ŷ(k ) = β 0 + ∑ β i xi (k) + ∑ ∑ β i1 i2 x i1 ( k ) x i2 ( k ) + · · · + ∑ ··· ∑ β i1 ···i p xi1 (k ) · · · xi p (k). (35)
i =1 i1 =1 i2 = i1 i1 =1 i p = i p −1

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 R2 and MSE values presented in Table 3.






9DOXH>9@






     
5DGLXV>PP@
(a) Proxy model with first order polynomial.

Figure 9. Cont.
Electronics 2020, 9, 1975 13 of 20







9DOXH>9@





     
5DGLXV>PP@
(b) Proxy model with second order polynomial.




9DOXH>9@





     
5DGLXV>PP@
(c) Proxy model with third order polynomial.
Figure 9. Polynomial regression (blue line) of different orders applied to the data presented in Table 2
(red dots).

Table 3. R2 and mean squared error (MSE) values for polynomial regression.

Polynomial Order R2 MSE

1 0.008 168.13
2 0.213 133.37
3 0.591 69.19

Given the nonlinearity of the data presented, it can be seen that although the R2 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:

ŷ(k) = β 0 x (k) + β 1 x (k )2 + β 2 x (k)3 + β 3 ln x (k). (36)

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
Electronics 2020, 9, 1975 14 of 20

in Figure 10. Notice that polynomial transformation is better suited to the data extracted from the
simulator, with an R2 of 0.99 and MSE of 0.79.






9DOXH>9@






     
5DGLXV>PP@
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

V (r ) = 249.58 + 20.25r − 1.97r2 + 0.07r3 − 15.26 ln(r ), (37)

which, after being validated, can be used to recast the problem of interest as

min V (r ) (38a)
r
s.t. V (r ) = 249.58 + 20.25r − 1.97r2 + 0.07r3 − 15.26 ln(r ) (38b)
0 ≤ r ≤ 10 (38c)

which can be easily solved by IPOPT [46], resulting in an optimized radius value of 0.9 mm.

4.1.2. 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

x k = σ (W k x k − 1 + b k ) (39)

in which σ(·) is the activation function, W k ∈ Rnk ×nk −1 is the matrix containing the weights of the
neural network for each layer, along with its biases bk ∈ Rnk . Thus, x k ∈ Rnk is the output of layer k.
When the ReLU function is defined as the activation function, σ (·) can be defined as

σ (y) = max{0, y}. (40)

Electronics 2020, 9, 1975 15 of 20

Following the methodology of [52], we can rewrite the ReLU operator as

W k x k −1 + b k = x k − s k (41)

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

L k ≤ W k x k −1 + b k = x k − s k ≤ U k . (42)

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:

x k ≤ U k zk (43a)
k k k
s ≤ − L (1 − z ) (43b)
k
zk ∈ {0, 1}n (43c)

Finally, we can rewrite the entire ReLU neural network as a MILP, in which the input layer is
given by
L0 ≤ x 0 ≤ U 0 , (44)

the hidden layers (∀k = 1, ..., K − 1) by

W k x k −1 + b k = x k − s k (45a)
x k , sk ≥ 0 (45b)
k k k
x ≤U z (45c)
s k ≤ − L k (1 − z k ) (45d)
nk
zk ∈ {0, 1} , (45e)

and the output layer (K) by

W K x K −1 + b K = x K (46)
K K K
L ≤x ≤U . (47)

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 x0 , 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.
Electronics 2020, 9, 1975 16 of 20





9DOXH>9@ 







     
5DGLXV>PP@
Figure 11. Proxy model with ReLU neural network (blue line) applied to the data presented in Table 2
(red dots).

4.1.3. 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.

5. 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
(R2 ). 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
Electronics 2020, 9, 1975 17 of 20

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.

Author Contributions: S.F.S.: Conceptualization, Data Acquisition, Formal Analysis, Validation, Writing—Original
Draft; L.O.S.: Conceptualization, Formal analysis, Validation, Writing—Original Draft; C.S.F.N.: Data Acquisition;
A.N.: Review & Editing; D.M.S.: Data Acquisition, Writing—Original Draft; F.G.d.L.: Conceptualization, Writing—
Original Draft; P.H.S.: Conceptualization, Writing—Original Draft; J.T.G.: Review & Editing; V.R.Q.L.: Review &
Editing. All authors have read and agreed to the published version of the manuscript.
Funding: Supported by project Smart following systems, Edge Computing and IoT Consortium, CONSORCIO
TC_TCUE18-20_004, CONVOCATORIA CONSORCIOTC. PLAN TCUE 2018–2020. Project managed by Fundación
General de la Universidad de Salamanca and co-financed with Junta de Castilla y León and FEDER funds.
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:

CAPES Coordination of Superior Level Staff Improvement

DFO Derivative-Free Optimization
FEM Finite Element Method
IPOPT Interior Point OPTimizer
MADS Mesh Adaptive Direct Search
MILP Mixed-integer linear problem
MSE Mean Squared Error
NLP Nonlinear Problem
NOMAD Nonlinear Optimization with the MADS Algorithm
PMSM Permanent Magnet Synchronous Machine
PSO Particle Swarm Optimization
ReLU Rectified Linear Unit
SSR Sum of Squared Residuals
SST Sum of Squares Total

References
1. Stefenon, S.F.; Branco, N.W.; Nied, A.; Bertol, D.W.; Finardi, E.C.; Sartori, A.; Meyer, L.H.; Grebogi, R.B.
Analysis of training techniques of ANN for classification of insulators in electrical power systems. IET Gener.
Transm. Distrib. 2020, 14, 1591–1597. [CrossRef]
2. Ribeiro, M.H.D.M.; Stefenon, S.F.; de Lima, J.D.; Nied, A.; Marini, V.C.; de Coelho, L.S. Electricity Price
Forecasting Based on Self-Adaptive Decomposition and Heterogeneous Ensemble Learning. Energies 2020,
13, 5190. [CrossRef]
3. Stefenon, S.F.; Silva, M.C.; Bertol, D.W.; Meyer, L.H.; Nied, A. Fault diagnosis of insulators from ultrasound
detection using neural networks. J. Intell. Fuzzy Syst. 2019, 37, 6655–6664. [CrossRef]
Electronics 2020, 9, 1975 18 of 20

4. Feng, G.; Lai, C.; Kar, N.C. Practical testing solutions to optimal stator harmonic current design for PMSM
torque ripple minimization using speed harmonics. IEEE Trans. Power Electron. 2017, 33, 5181–5191.
[CrossRef]
5. Stefenon, S.F.; Americo, J.P.; Meyer, L.H.; Grebogi, R.B.; Nied, A. Analysis of the electric field in porcelain
pin-type insulators via finite elements software. IEEE Lat. Am. Trans. 2018, 16, 2505–2512. [CrossRef]
6. Corso, M.P.; Stefenon, S.F.; Couto, V.F.; Cabral, S.H.L.; Nied, A. Evaluation of methods for electric field
calculation in transmission lines. IEEE Lat. Am. Trans. 2018, 16, 2970–2976. [CrossRef]
7. Fonteyn, K.; Belahcen, A.; Kouhia, R.; Rasilo, P.; Arkkio, A. FEM for directly coupled magneto-mechanical
phenomena in electrical machines. IEEE Trans. Magn. 2010, 46, 2923–2926. [CrossRef]
8. Parreira, B.; Rafael, S.; Pires, A.; Branco, P.C. Obtaining the magnetic characteristics of an 8/6 switched
reluctance machine: From FEM analysis to the experimental tests. IEEE Trans. Ind. Electron. 2005,
52, 1635–1643. [CrossRef]
9. Lee, D.; Park, G.J.; Son, B.; Jung, H.C. Efficiency improvement of IPMSG in the electric power generating
system of a range-extended electric vehicle. IET Electr. Power Appl. 2019, 13, 943–950. [CrossRef]
10. Meessen, K.J.; Thelin, P.; Soulard, J.; Lomonova, E. Inductance calculations of permanent-magnet
synchronous machines including flux change and self-and cross-saturations. IEEE Trans. Magn. 2008,
44, 2324–2331. [CrossRef]
11. Luukko, J.; Pyrhönen, J. Selection of the parameters of a permanent magnet synchronous machine by using
nonlinear optimisation. IET Electr. Power Appl. 2007, 1, 255–263.:20060298. [CrossRef]
12. Wu, J.; Wang, J.; Gan, C.; Sun, Q.; Kong, W. Efficiency optimization of PMSM drives using field-circuit
coupled FEM for EV/HEV applications. IEEE Access 2018, 6, 15192–15201. [CrossRef]
13. Jeong, C.L.; Kim, Y.K.; Hur, J. Optimized design of PMSM with hybrid-type permanent magnet for improving
performance and reliability. IEEE Trans. Ind. Appl. 2019, 55, 4692–4701. [CrossRef]
14. Wallmark, O.; Nikouei, M. Dc-link and machine design considerations for resonant controllers adopted in
automotive PMSM drives. IET Electr. Syst. Transp. 2019, 10, 75–80. [CrossRef]
15. Zhang, Q.; Cheng, S.; Wang, D.; Jia, Z. Multiobjective design optimization of high-power circular winding
brushless DC motor. IEEE Trans. Ind. Electron. 2017, 65, 1740–1750. [CrossRef]
16. Zhang, Y.; Cao, W.; McLoone, S.; Morrow, J. Design and flux-weakening control of an interior permanent
magnet synchronous motor for electric vehicles. IEEE Trans. Appl. Supercond. 2016, 26, 1–6. [CrossRef]
17. Hong, D.K.; Hwang, W.; Lee, J.Y.; Woo, B.C. Design, analysis, and experimental validation of a permanent
magnet synchronous motor for articulated robot applications. IEEE Trans. Magn. 2017, 54, 1–4. [CrossRef]
18. Liu, X.; Chen, H.; Zhao, J.; Belahcen, A. Research on the performances and parameters of interior PMSM
used for electric vehicles. IEEE Trans. Ind. Electron. 2016, 63, 3533–3545. [CrossRef]
19. Ran, X.; Zhao, M.; Shang, J.; He, C. Dynamic analysis and mechanical structure improvement of
submersible rotor. IEEE Access 2019, 7, 51640–51647. [CrossRef]
20. Paula, G.T.; Monteiro, J.R.B.A.; Almeida, T.E.P.; Santana, M.P. Different slot configurations for direct-drive
pm brushless machines. IEEE Lat. Am. Trans. 2015, 13, 634–639. [CrossRef]
21. Seman, L.O.; Miyatake, L.K.; Camponogara, E.; Giuliani, C.M.; Vieira, B.F. Derivative-free parameter tuning
for a well multiphase flow simulator. J. Pet. Sci. Eng. 2020, 192, 107288. [CrossRef]
22. Ding, S.; Hang, J.; Wei, B.; Wang, Q. Modelling of supercapacitors based on SVM and PSO algorithms.
IET Electr. Power Appl. 2018, 12, 502–507. [CrossRef]
23. Nogueira, A.F.L. Calculation of global magnetic forces using analytical and finite element solutions. IEEE Lat.
Am. Trans. 2016, 14, 2365–2371. [CrossRef]
24. Tong, L.; Li, X.; Hu, J.; Ren, L. A PSO optimization scale-transformation stochastic-resonance algorithm with
stability mutation operator. IEEE Access 2017, 6, 1167–1176. [CrossRef]
25. Huynh, D.; Dunnigan, M. Parameter estimation of an induction machine using advanced particle swarm
optimisation algorithms. IET Electr. Power Appl. 2010, 4, 748–760. [CrossRef]
26. Gaing, Z.L.; Lin, C.H.; Tsai, M.H.; Hsieh, M.F.; Tsai, M.C. Rigorous design and optimization of brushless
PM motor using response surface methodology with quantum-behaved PSO operator. IEEE Trans. Magn.
2013, 50, 1–4, doi:10.1109/TMAG.2013.2279672. [CrossRef]
27. Niu, S.; Chen, N.; Ho, S.; Fu, W. Design optimization of magnetic gears using mesh adjustable finite-element
algorithm for improved torque. IEEE Trans. Magn. 2012, 48, 4156–4159. [CrossRef]
Electronics 2020, 9, 1975 19 of 20

28. Camponogara, E.; Seman, L.O. Control Optimization of Pump Cycles in Onshore Oilfields with Network
and Electric Power Constraints. J. Energy Resour. Technol. 2021, 143. [CrossRef]
29. Nalbant, M.; Gökkaya, H.; Sur, G. Application of Taguchi method in the optimization of cutting parameters
for surface roughness in turning. Mater. Des. 2007, 28, 1379–1385. [CrossRef]
30. Li, Y.; Liang, D.; Fan, Y.; Xin, J.; Zhuang, W. Application of Taguchi method and finite element analysis in
optimization of automobile roof technical parameters. In Proceedings of the IEEE International Conference
on Smart Manufacturing, Industrial & Logistics Engineering (SMILE), Hangzhou, China, 20–21 April 2019;
pp. 75–79. [CrossRef]
31. Chen, D.C.; You, C.S.; Nian, F.L.; Guo, M.W. Using the Taguchi Method and Finite Element Method to
Analyze a Robust New Design for Titanium Alloy Prick Hole Extrusion. Procedia Eng. 2011, 10, 82–87.
[CrossRef]
32. Hwang, K.Y.; Song, B.K.; Kwon, B.I. Asymmetric dual winding three-phase PMSM for fault tolerance of
overheat in electric braking system of autonomous vehicle. IET Electr. Power Appl. 2019, 13, 1891–1898.
[CrossRef]
33. Stéfano, F.S.; Ademir, N. FEM Applied to Evaluation of the Influence of Electric Field on Design of the Stator
Slots in PMSM. IEEE Lat. Am. Trans. 2019, 17, 590–596. [CrossRef]
34. Krause, P.; Wasynczuk, O.; Sudhoff, S.D.; Pekarek, S. Introduction to the Design of Electric Machinery.
In Analysis of Electric Machinery and Drive Systems; Wiley-IEEE Press: Hoboken, NJ, USA, 2013; Volume 3,
pp. 583–622.
35. Velasco, J.; Frascella, R.; Albarracín, R.; Burgos, J.C.; Dong, M.; Ren, M.; Yang, L. Comparison of positive
streamers in liquid dielectrics with and without nanoparticles simulated with finite-element software.
Energies 2018, 11, 361. [CrossRef]
36. Habibinia, D.; Rostami, N.; Feyzi, M.R.; Soltanipour, H.; Pyrhönen, J. New finite element based method for
thermal analysis of axial flux interior rotor permanent magnet synchronous machine. IET Electr. Power Appl.
2019, 14, 464–470. [CrossRef]
37. Wu, J.; Zhu, X.; Xu, L.; Quan, L.; Fan, D.; Yang, J. Optimisation design of a flux memory motor based on a
new non-linear MC-DRN model. IET Electr. Power Appl. 2019, 13, 2035–2043. [CrossRef]
38. Xiao, Y.; Zhou, L.; Wang, J.; Liu, J. Design and performance analysis of magnetic slot wedge application
in double-fed asynchronous motor-generator by finite-element method. IET Electr. Power Appl. 2018,
12, 1040–1047. [CrossRef]
39. Giurgea, S.; Fodorean, D.; Cirrincione, G.; Miraoui, A.; Cirrincione, M. Multimodel optimization based on
the response surface of the reduced FEM simulation model with application to a PMSM. IEEE Trans. Magn.
2008, 44, 2153–2157. [CrossRef]
40. Shakouri, H.; Nadimi, R.; Ghaderi, S.F. Investigation on objective function and assessment rule in fuzzy
regressions based on equality possibility, fuzzy union and intersection concepts. Comput. Ind. Eng.
2017, 110, 207–215. [CrossRef]
41. Reid, S.; Tibshirani, R.; Friedman, J. A study of error variance estimation in Lasso regression. Stat. Sin.
2016, 26, 35–67. [CrossRef]
42. Seman, L.O.; Rodrigues Machado, V.H.; Koehler, L.A.; Camponogara, E. A framework to estimate dwell
time of BRT systems using fuzzy regression. J. Intell. Fuzzy Syst. 2020, 38, 5279–5293. [CrossRef]
43. Stefenon, S.F.; Ribeiro, M.H.D.M.; Nied, A.; Mariani, V.C.; dos Santos Coelho, L.; da Rocha, D.F.M.;
Grebogi, R.B.; de Barros Ruano, A.E. Wavelet group method of data handling for fault prediction in
electrical power insulators. Int. J. Electr. Power Energy Syst. 2020, 123, 106269. [CrossRef]
44. Kasburg, C.; Stefenon, S.F. Deep learning for photovoltaic generation forecast in active solar trackers.
IEEE Lat. Am. Trans. 2019, 17, 2013–2019. [CrossRef]
45. Stefenon, S.F.; Freire, R.Z.; Coelho, L.d.S.; Meyer, L.H.; Grebogi, R.B.; Buratto, W.G.; Nied, A. Electrical
insulator fault forecasting based on a wavelet neuro-fuzzy system. Energies 2020, 13, 484. [CrossRef]
46. Wächter, A.; Biegler, L.T. On the implementation of an interior-point filter line-search algorithm for large-scale
nonlinear programming. Math. Program. 2006, 106, 25–57. [CrossRef]
47. Le Digabel, S. Algorithm 909: NOMAD: Nonlinear optimization with the MADS algorithm. ACM Trans.
Math. Softw. 2011, 37. [CrossRef]
48. McCormick, G.P. Computability of global solutions to factorable nonconvex programs: Part I—Convex
underestimating problems. Math. Program. 1976, 10, 147–175. [CrossRef]
Electronics 2020, 9, 1975 20 of 20

49. Stefenon, S.F.; Steinheuser, D.F.; da Silva, M.P.; Ferreira, F.C.S.; Klaar, A.C.R.; de Souza, K.E.; Júnior, A.G.;
Venção, A.T.; Branco, R.; Yamaguchi, C.K. Application of Active Methodologies in Engineering Education
through the Integrative Evaluation at the Universidade do Planalto Catarinense, Brazil. Interciencia
2019, 44, 408–413.
50. Benoit, K. Linear Regression Models with Logarithmic Transformations. 2011. Available online:
https://kenbenoit.net/assets/courses/ME104/logmodels2.pdf (accessed on 15 September 2020).
51. Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual. Available online: https://www.gurobi.com/
wp-content/plugins/hd_documentations/documentation/9.1/refman.pdf (accessed on 10 September 2020).
52. Grimstad, B.; Andersson, H. ReLU networks as surrogate models in mixed-integer linear programs.
Comput. Chem. Eng. 2019, 131, 106580. [CrossRef]

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

c 2020 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 (http://creativecommons.org/licenses/by/4.0/).