Optimal Design of IPMSM for EV Using Subdivided Kriging Multi-Objective Optimization

: In this paper, subdivided kriging multi-objective optimization (SKMOO) is proposed for the optimal design of interior permanent magnet synchronous motor (IPMSM). The SKMOO with surrogate kriging model can obtain a uniform and accurate pareto front set with a reduced computation cost compared to conventional algorithms which directly adds the solution in the objective function area. In other words, the proposed algorithm uses a kriging surrogate model, so it is possible to know which design variables have the value of the objective function on the blank space. Therefore, the solution can be added directly in the objective function area. In the SKMOO algorithm, a non-dominated sorting method is used to ﬁnd the pareto front set and the ﬁll blank method is applied to prevent premature convergence. In addition, the subdivided kriging grid is proposed to make a well-distributed and more precise pareto front set. Superior performance of the SKMOO is conﬁrmed by compared conventional multi objective optimization (MOO) algorithms with test functions and are applied to the optimal design of IPMSM for electric vehicle.


Introduction
Recently, the problem of environment pollution related to vehicles has been a constant issue around the world. In response to these issues, electric vehicles (EVs) are receiving a lot of attention [1,2]. The use of a permanent magnet synchronous motor (PMSM) in EV driving motors is increasing because the PMSMs have high efficiency, high torque power density, and a high ratio of output torque to motor weight [3][4][5].
The PMSMs are classified as the interior permanent magnet synchronous motor (IPMSM) with permanent magnet inserted inside the rotor and the surface mounted permanent magnet synchronous motor (SPMSM) with a permanent magnet attached to the surface of the rotor. The effective air gap in the d-axis and q-axis inductance are not the same in IPMSM [6]. Therefore, the IPMSM has the torque caused by the magnet and reluctance torque. The IPMSM torque greater than surface mounted permanent magnet synchronous motor [7]. However, such a permanent magnet type motor has a disadvantage of a high cogging torque, a high torque ripple, and a high total harmonic distortion (THD) of back electromotive force (EMF) which causes a reduction in the motor's noise, vibration, and efficiency. This is highly important to EV for comfortable driving [8]. Therefore, we aim to design the IPMSM with reduced cogging torque and THD of back EMF.
The IPMSM can improve the cogging torque characteristics according to the rotor shape that shows a complex shape in which the objective function region has multiple peaks [8]. However, the IPMSM has nonlinear characteristic. Therefore, finite element method (FEM) analysis is required because the accurate analysis is difficult by using simple mathematical models due to the nonlinear characteristic of the IPMSM [9]. The actual motor design using FEM analysis should meet several requirements. It takes a lot of time and effort to simultaneously consider a multi-objective function in motor design that requires the FEM interpretation [10,11]. To solve this problem, a lot of studies of multiobjective optimization (MOO) are in progress [12][13][14][15][16]. The most widely used algorithms are multi-objective particle swarm optimization (MOPSO) and non-dominated sorting genetic algorithm-II (NSGA-II). However, the MOPSO and NSGA-II still need many function calls to get the accurate and well distributed pareto front set. Therefore, subdivided kriging multi-objective optimization (SKMOO) is proposed which applies the surrogate model to reduce function calls.
Outstanding performance of the proposed algorithm is verified through comparison of test results with the MOPSO and NSGA-II. The optimal design was derived by applying the SKMOO to design of the IPMSM with reduced cogging torque and THD of back EMF.

Proposed Algorithm
MOO performance can be determined based on the accuracy, uniform distribution, and number of function calls. The SKMOO uses the following methods to improve algorithm performance.

Latin Hypercube Sampling
Initial samples are created using a Latin hypercube sampling (LHS) method. This method ensures that samples are generated evenly in the whole problem areas while maintaining randomness to generate samples efficiently. The problem area is divided to the number of samples multiply the number of samples. It is a method of sampling so that each row and column of a divided grid contains only one sample. To distribute the sample more evenly, the entire area is divided into several layers and into several blocks. LHS is implemented on each later to generate samples. In addition to these methods, there are many improved LHS methods. However, this paper used the existing LHS method because initial sample generation is not a main idea [17][18][19][20][21].

Surrogate Model Creation
A surrogate model is created using the kriging method with initial samples. Kriging is a method of predicting unknown values as a linear combination. It uses the value of the objective function already known to estimate the shape and objective function values of the entire problem area [22][23][24]. The number of samples is very important because the kriging method is predicted using the initial samples. If the samples are too small, the problem area will not be predicted properly, and if there are too many samples, unnecessary sample generation will increase the number of function calls. Therefore, it is essential to select the appropriate number of samples. For this paper, appropriate initial samples were selected by referring to the [25].

Non-Dominated Sorting
The non-dominated sorting compares all samples. The problem of minimizing the objective functions f 1 , f 2 is that if "A" solution has a smaller value in both f 1 , f 2 objective functions than "B" solution, the "B" solution is dominated by "A" solution. On the other hand, when comparing "A" solution and "C" solution, the "A" solution is small in f 1 and "C" is small in f 2 . When the dominant relation of samples for each objective function value is indistinguishable, such as the "A" and "C" solution, they are called non-dominant relations and the samples are defined as non-dominated solution [26][27][28]. Furthermore, this set of non-dominated solutions is called the pareto front set. Therefore, the initial samples are classified through a non-dominated sorting and the pareto front set is extracted.

Euclidean Distance Calculation
The Euclidean distance means the distance of adjacent solutions. It is calculated to obtain a uniform pareto front set. ED max is defined as the maximum Euclidean distance. The ED max is calculated as follows where Norm i is range of ith objective function to normalize, N n is the number of solutions, and N o is the number of objective functions. In Figure 1, the dotted line means the Euclidean distance. The predicted solution is added directly between the two solutions with the ED max . In addition, the spline interpolation is used to increase the accuracy of the solution and has a more continuous pareto front set. Therefore, the predicted solution is filled the blank of pareto set, similar to a green point in Figure 1.

Euclidean Distance Calculation
The Euclidean distance means the distance of adjacent solutions. It is calculated to obtain a uniform pareto front set. EDmax is defined as the maximum Euclidean distance. The EDmax is calculated as follows where Normi is range of ith objective function to normalize, Nn is the number of solutions, and No is the number of objective functions. In Figure 1, the dotted line means the Euclidean distance. The predicted solution is added directly between the two solutions with the EDmax. In addition, the spline interpolation is used to increase the accuracy of the solution and has a more continuous pareto front set. Therefore, the predicted solution is filled the blank of pareto set, similar to a green point in Figure 1.

Inverse Searching
The combination of design variables should be explored to generate the predicted solution. Inverse searching is possible because it uses a kriging surrogate model. The spatial distribution is obtained through the kriging on a predefined grid area from calculated samples without additional function calls [29]. In Figure 2, the predicted solution of objective function area is placed at the same position in the f1 and f2 of the design variable area. The most similar combination of design variables is found by comparing the objective function value of the surrogate kriging model with the predicted value. The formula for inverse searching is as follows where x1, x2 are the design variable, i, j are values of uniformly divided design variables, respectively ∈ 1,2, ⋯ , , ∈ 1,2, ⋯ , where Nx1, Nx2 are number of each design variable, fk,pre is the value of the kth objective function of the predicted solution.

Inverse Searching
The combination of design variables should be explored to generate the predicted solution. Inverse searching is possible because it uses a kriging surrogate model. The spatial distribution is obtained through the kriging on a predefined grid area from calculated samples without additional function calls [29]. In Figure 2, the predicted solution of objective function area is placed at the same position in the f 1 and f 2 of the design variable area. The most similar combination of design variables is found by comparing the objective function value of the surrogate kriging model with the predicted value. The formula for inverse searching is as follows where x 1 , x 2 are the design variable, i, j are values of uniformly divided design variables, where N x1 , N x2 are number of each design variable, f k,pre is the value of the kth objective function of the predicted solution.

Subdivided Kriging Grid
The real solution is obtained using inverse searching. However, there is a difference between the real solution and the predicted solution because the real solution is created only on a kriging grid line, similar to red point in the design variable area of Figure 3. Therefore, the SKMOO is proposed to reduce difference. The right of design variable area in Figure 3 shows the subdivided kriging grid. The red dotted line makes the pareto front set more uniform and the solution more accurate. The more kriging grid, the higher accuracy of the solution; however, the reason why the grid cannot be infinitely finely divided is because all the values of all grid points must be compared, which can take a long time. The proper kriging grid setting is required considering the accuracy of the solution and interpretation time [25]. Therefore, in this paper, the kriging grid is not finely divided from the beginning of algorithm, it is divided later in the algorithm to reduce the interpretation time and increase the accuracy of the solution. In addition, if the two solutions are too close to each other, the solutions are deleted, and new solution is added between the two solutions.

Fill Blank Method
Conventional MOO algorithms used mutation operations to prevent premature convergence [30]. However, existing mutation operations generate samples randomly throughout the problem area, which can lead to redundant exploration of the area already explored. To solve this problem, a fill blank method is used [31,32]. This is a method to add samples where the samples are least located for efficient exploration. When the samples are placed irregularly in the design variable area, a sample is added to the area to increase the variety of solutions. This method is greatly efficient because it does not overlap the area already explored.

Subdivided Kriging Grid
The real solution is obtained using inverse searching. However, there is a difference between the real solution and the predicted solution because the real solution is created only on a kriging grid line, similar to red point in the design variable area of Figure 3. Therefore, the SKMOO is proposed to reduce difference. The right of design variable area in Figure 3 shows the subdivided kriging grid. The red dotted line makes the pareto front set more uniform and the solution more accurate. The more kriging grid, the higher accuracy of the solution; however, the reason why the grid cannot be infinitely finely divided is because all the values of all grid points must be compared, which can take a long time. The proper kriging grid setting is required considering the accuracy of the solution and interpretation time [25]. Therefore, in this paper, the kriging grid is not finely divided from the beginning of algorithm, it is divided later in the algorithm to reduce the interpretation time and increase the accuracy of the solution. In addition, if the two solutions are too close to each other, the solutions are deleted, and new solution is added between the two solutions.

Subdivided Kriging Grid
The real solution is obtained using inverse searching. However, there is a difference between the real solution and the predicted solution because the real solution is created only on a kriging grid line, similar to red point in the design variable area of Figure 3. Therefore, the SKMOO is proposed to reduce difference. The right of design variable area in Figure 3 shows the subdivided kriging grid. The red dotted line makes the pareto front set more uniform and the solution more accurate. The more kriging grid, the higher accuracy of the solution; however, the reason why the grid cannot be infinitely finely divided is because all the values of all grid points must be compared, which can take a long time. The proper kriging grid setting is required considering the accuracy of the solution and interpretation time [25]. Therefore, in this paper, the kriging grid is not finely divided from the beginning of algorithm, it is divided later in the algorithm to reduce the interpretation time and increase the accuracy of the solution. In addition, if the two solutions are too close to each other, the solutions are deleted, and new solution is added between the two solutions.

Fill Blank Method
Conventional MOO algorithms used mutation operations to prevent premature convergence [30]. However, existing mutation operations generate samples randomly throughout the problem area, which can lead to redundant exploration of the area already explored. To solve this problem, a fill blank method is used [31,32]. This is a method to add samples where the samples are least located for efficient exploration. When the samples are placed irregularly in the design variable area, a sample is added to the area to increase the variety of solutions. This method is greatly efficient because it does not overlap the area already explored.

Fill Blank Method
Conventional MOO algorithms used mutation operations to prevent premature convergence [30]. However, existing mutation operations generate samples randomly throughout the problem area, which can lead to redundant exploration of the area already explored. To solve this problem, a fill blank method is used [31,32]. This is a method to add samples where the samples are least located for efficient exploration. When the samples are placed irregularly in the design variable area, a sample is added to the area to increase the variety of solutions. This method is greatly efficient because it does not overlap the area already explored.  Figure 4 represent the flow chart of the proposed algorithm. First, the parameters, such as the grid size, number of populations, number of mutations, constraints, and so on, are defined. Initial samples allow the surrogate model to have uniform and irregularity using the Latin hypercube sampling to ensure even accuracy overall. Pareto front set is obtained using non-dominated sorting where the non-dominated solution is not dominated by another solution for any one objective function. In the Euclidean distance calculation, the Euclidean distance is normalized so that it does not have biased values for one objective function. Predicted solution is added using calculated ED max and design variable is explored through the inverse searching. The Fill blank method is used to diversify the pareto front set. After the set iteration is executed, the kriging grid is subdivided and then the solution is added and deleted to improve the characteristic of algorithm. The inverse searching and the Fill blank method are proceeded. Finally, if the specified convergence condition is satisfied, the algorithm is finished. Figure 4 represent the flow chart of the proposed algorithm. Fi such as the grid size, number of populations, number of mutations, co are defined. Initial samples allow the surrogate model to have unifo using the Latin hypercube sampling to ensure even accuracy overall obtained using non-dominated sorting where the non-dominated so nated by another solution for any one objective function. In the Euclid lation, the Euclidean distance is normalized so that it does not have b objective function. Predicted solution is added using calculated EDmax is explored through the inverse searching. The Fill blank method is u pareto front set. After the set iteration is executed, the kriging grid is s the solution is added and deleted to improve the characteristic of alg searching and the Fill blank method are proceeded. Finally, if the sp condition is satisfied, the algorithm is finished.

Verification of Proposed Algorithm
The SKMOO, the MOPSO, and NSGA-II are applied to the test fu tive functions are defined as follows

Verification of Proposed Algorithm
The SKMOO, the MOPSO, and NSGA-II are applied to the test functions. The objective functions are defined as follows Maximize f = ( f 1 (x), f 2 (x, y)) (7) where the range of x is from 0.1 to 1 and y is from 0 to 5, respectively. Figures 5 and 6 show the results of test function 1 and test function 2. The solution of SKMOO is more similar to the real pareto set and uniform distribution than conventional algorithms.
where the range of x 1 and x 2 are from −4 to 4 in Equations (9) and (10). The range of x 1 , x 2 , and x 3 are from 0 to 1 in Equations (11) and (12).
where d i is the shortest Euclidean distance in the objective function area between measured pareto front set and real pareto front set. d is average value of all distance between adjacent solutions. In a multi-objective problem, it is better to distribute throughout the problem area so that users can choose a solution from a variety of solutions. GD and SP values mean better performance if there are values close to 0. Table 1 shows numerical indicators. It shows that the proposed algorithm is closer to zero than the other two algorithms. In addition, since the number of function call means the number of FEM analyses, the smaller the function call, the fewer FEM there are to be solved. Therefore, the proposed algorithm is superior in terms of time. The tests were repeated 100 times each. Table 1 summarizes the results of the test and shows that the SKMOO has superior performance on GD, SP, and the number of function calls than NSGA-II and MOPSO, respectively. values mean better performance if there are values close to 0. Table 1 shows numeric indicators. It shows that the proposed algorithm is closer to zero than the other two alg rithms. In addition, since the number of function call means the number of FEM analyse the smaller the function call, the fewer FEM there are to be solved. Therefore, the propose algorithm is superior in terms of time. The tests were repeated 100 times each. Table  summarizes the results of the test and shows that the SKMOO has superior performan on GD, SP, and the number of function calls than NSGA-II and MOPSO, respectively.

Optimal Design of IPMSM for EV
In the past several years, environmental pollution caused by many emissions has become increasingly critical. Many countries are trying to reduce emissions. The vehicles are one of the largest causes of pollution [33]. As a result, the eco-friendly EVs have been in the spotlight [34]. EVs must satisfy the power and efficiency in limited space. Therefore, high torque power density and high efficiency characteristics are required. The IPMSM is an appropriate motor to meet these characteristics [35]. The SKMOO is applied to optimal design of IPMSM for EV to verify the practicability. The commercial FEM tool, JMAG, is used to analyze of IPMSM. Load and no-load characteristics are analyzed by using the FEM for accurate interpretation. The objective functions are determined to reduce the cogging torque and THD of back EMF while satisfying the rated torque to reduce noise and vibration. The FEM was carried out using (i5-8600K CPU, 16.0GB RAM) computer. The average number of function calls in the four test functions is 1396.7 for NSGA-II, 1685.8 for MOPSO, and 760.3 for SKMOO. The one interpretation of FEM analysis takes about 1 min 40 s for optimization. One function call means one FEM interpretation. Therefore, the proposed algorithm takes about 21.1 h, the NSGA-II takes about 38.8 h, and the MOPSO takes about 46.8 h. In terms of time, the proposed algorithm improved by 45.6% and 54.9%, respectively. The design specifications are tabulated on Table 2. The many papers use eight-pole type motors to EV [36][37][38]. In the paper based on FEM analysis, however, the model is demonstrated indirectly because the reference model of the proposed paper progressed with reference to the [39]. The accuracy of FEM analysis in existing research cases is very high compared to the simulation [40][41][42]. Furthermore, there are many papers that validate the performance of the algorithm with only FEM interpretation [43][44][45][46].  Figure 9 shows the configuration of reference model. Figure 10 shows the configuration of optimal model and design variables. The design variables, θ 1 and θ 2 , are the angles of the first-and second-layer magnet. The range of the design variables are considered not to overlap the first-and second-layer magnet, and the range is limited according to fix the length and thickness of the magnet, the length of the center posts, and the length of the bridge. The corresponding range of θ 1 is from 47 • to 55.6 • and θ 2 is from 45 • to 70.5 • , respectively. The H_length is the hole length between the bridge of the first-and second-layer magnet. The range of the H_length is from 0.1 to 2.6.  The pareto front set of IPMSM using the SKMOO is shown in Figure 11. The optimal solution is selected assuming that the ratio of the cogging torque and the THD of back EMF are the same. In other words, each objective function value is multiplied by 0.5 and the solution with the smallest sum is selected. At this time, the value is divided the difference between the maximum and minimum values of each objective function to avoid having a biased value for one objective function. Therefore, the 2.70 Nm in the cogging torque and the 3.45% in the THD are selected as the optimal solution. The model compares all conditions in the same except design variables.  The pareto front set of IPMSM using the SKMOO is shown in Figure 11. The optimal solution is selected assuming that the ratio of the cogging torque and the THD of back EMF are the same. In other words, each objective function value is multiplied by 0.5 and the solution with the smallest sum is selected. At this time, the value is divided the difference between the maximum and minimum values of each objective function to avoid having a biased value for one objective function. Therefore, the 2.70 Nm in the cogging torque and the 3.45% in the THD are selected as the optimal solution. The model compares all conditions in the same except design variables. The pareto front set of IPMSM using the SKMOO is shown in Figure 11. The optimal solution is selected assuming that the ratio of the cogging torque and the THD of back EMF are the same. In other words, each objective function value is multiplied by 0.5 and the solution with the smallest sum is selected. At this time, the value is divided the difference between the maximum and minimum values of each objective function to avoid having a biased value for one objective function. Therefore, the 2.70 Nm in the cogging torque and the 3.45% in the THD are selected as the optimal solution. The model compares all conditions in the same except design variables.  Figure 12a shows the comparison of torque between the initial model and the optimal model. The average torque increases 0.47% more compared to the initial model. Figure  12b shows the comparison of back EMF. The THD of back EMF is 4.96% for the initial model and 2.70% for the optimal model, with an optimal model being 45.56% lower. Figure 12c shows a comparison of the cogging torques. The cogging torque is 4.75 Nm for the initial model and 3.45 Nm for the optimal model, with the latter being 27.37% lower.  Figure 12a shows the comparison of torque between the initial model and the optimal model. The average torque increases 0.47% more compared to the initial model. Figure 12b shows the comparison of back EMF. The THD of back EMF is 4.96% for the initial model and 2.70% for the optimal model, with an optimal model being 45.56% lower. Figure 12c shows a comparison of the cogging torques. The cogging torque is 4.75 Nm for the initial model and 3.45 Nm for the optimal model, with the latter being 27.37% lower.  Figure 12a shows the comparison of torque between the initial model and the optimal model. The average torque increases 0.47% more compared to the initial model. Figure  12b shows the comparison of back EMF. The THD of back EMF is 4.96% for the initial model and 2.70% for the optimal model, with an optimal model being 45.56% lower. The results of optimal design are shown in Table 3. The optimal design is satisfied by a rated torque of 285.5 Nm, constraint of THD of 5.00%, and constraint of cogging torque of 4.85 Nm. Figure 13 shows the magnetic flux density of optimal model.  The results of optimal design are shown in Table 3. The optimal design is satisfied by a rated torque of 285.5 Nm, constraint of THD of 5.00%, and constraint of cogging torque of 4.85 Nm. Figure 13 shows the magnetic flux density of optimal model.  The mechanical stress analysis is performed to confirm the force exerted on the steel due to the centrifugal force of the permanent magnet inserted inside when the rotor rotates. The conditions and specifications are tabulated in Table 4. Figure 14 shows the Mises stress at rated and maximum speed. Table 5 shows the analysis result. The Mises stress value is 14.78 MPa in rated speed, the maximum Mises stress value is 262.07 MPa, and the yield stress of the steel is 450 MPa. Therefore, the optimal model is safe from breakage.  The mechanical stress analysis is performed to confirm the force exerted on the steel due to the centrifugal force of the permanent magnet inserted inside when the rotor rotates. The conditions and specifications are tabulated in Table 4. Figure 14 shows the Mises stress at rated and maximum speed. Table 5 shows the analysis result. The Mises stress value is 14.78 MPa in rated speed, the maximum Mises stress value is 262.07 MPa, and the yield stress of the steel is 450 MPa. Therefore, the optimal model is safe from breakage.

Conclusions
In this paper, SKMOO based on a surrogate kriging model is presented to solve MOO problems. The SKMOO can reduce computation cost and show superior performances in the MOO problems. The proposed algorithm is validated by comparing conventional algorithms. When calculating the average amount of reduction for each indicator in four test functions, the proposed algorithm was smaller than NSGA-II about 66.04% in GD, about 59.92% in SP, and about 41.28% in function call, and also smaller than MOPSO by about 52.51% in GD, 54.15% in SP, and 49.31% in function call. Finally, SKMOO is applied to the optimal design of IPMSM and demonstrated feasibility on practical electrical machines by showing improved design results. The results of optimal design model reduced 27.37% in the cogging torque and 45.56% in the THD of back EMF compared to initial model.
Finally, an important contribution of this paper is to consider multi-objective functions simultaneously with only small function calls, and the performance of the algorithm, which was verified through comparison with existing algorithms and applied to the design of electrical appliances. Therefore, the proposed algorithm is expected to be widely used for the MOO problems, which should satisfy multi-objective functions simultaneously.

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

Conclusions
In this paper, SKMOO based on a surrogate kriging model is presented to solve MOO problems. The SKMOO can reduce computation cost and show superior performances in the MOO problems. The proposed algorithm is validated by comparing conventional algorithms. When calculating the average amount of reduction for each indicator in four test functions, the proposed algorithm was smaller than NSGA-II about 66.04% in GD, about 59.92% in SP, and about 41.28% in function call, and also smaller than MOPSO by about 52.51% in GD, 54.15% in SP, and 49.31% in function call. Finally, SKMOO is applied to the optimal design of IPMSM and demonstrated feasibility on practical electrical machines by showing improved design results. The results of optimal design model reduced 27.37% in the cogging torque and 45.56% in the THD of back EMF compared to initial model.
Finally, an important contribution of this paper is to consider multi-objective functions simultaneously with only small function calls, and the performance of the algorithm, which was verified through comparison with existing algorithms and applied to the design of electrical appliances. Therefore, the proposed algorithm is expected to be widely used for the MOO problems, which should satisfy multi-objective functions simultaneously.