Multiobjective Optimization Design of Spherical Axial Split-Phase Permanent Maglev Flywheel Machine Based on Kriging Model

In order to improve the torque and suspension performance and reduce the torque ripple of spherical axial split-phase permanent maglev flywheel machine (S-AP-MF), a multi-objective optimization design method based on kriging model was proposed. A joint simulation model of S-AP-MF is constructed using Maxwell and Isight, and the sample data set is obtained. Then, latin hypercube method is used to select samples evenly and hierarchically to conduct sensitivity analysis on S-AP-MF structural parameters and obtain the variables to be optimized. Based on this, to improve the model accuracy and optimization efficiency, the kriging model and neighborhood cultivation genetic algorithm (NCGA) are used to obtain the optimal parameter set of S-AP-MF. Finally, the torque performance and levitation force characteristics before and after optimization are analyzed and compared by finite element method. The results show that the torque and average levitation force of the optimized S-AP-MF are increased by 9.54% and 5.51% respectively, and the torque ripple is reduced by 27.95%, which verifies the effectiveness of the proposed method.


I. INTRODUCTION
Flywheel energy storage system (FESS) has the advantages of zero pollution, high energy storage density and rapid charge and discharge. It can be used in the fields of uninterruptible power supply, electric vehicle, rail transit, wind power system frequency modulation and so on [1]- [4]. The machine used in flywheel energy storage as the core component of FESS directly affects the energy storage density and operation reliability of the system. Considering the integration and reliable and efficient operation of FESS, a new type of axial split-phase permanent maglev flywheel machine (AP-MF) is proposed in [5], and the characteristics of natural decoupling of electromagnetic torque and levitation force of AP-MF are verified by experiments. AP-MF sets a magnet separator The associate editor coordinating the review of this manuscript and approving it for publication was Shihong Ding . between the suspension pole and the torque pole to weaken the coupling effect of the levitation magnetic circuit and the torque magnetic circuit. The permanent magnet arranged between the two phases is also beneficial to improve the suspension stability of the machine. However, the magnetic pole of the machine adopts cylindrical salient pole structure, which leads to the uneven distribution of magnetic flux density in space. Changing all the magnetic poles of the machine into a spherical structure can effectively improve the global uniformity of air gap magnetic density and obtain better torque and suspension performance. Therefore, based on the characteristics of the improved spherical machine, this paper studies and optimizes the torque and suspension performance of the spherical AP-MF (S-AP-MF).
The finite element model (FEM) analysis method is often used in the optimization design of the machine. The high-order nonlinear machine model FEM with complex VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ structure usually needs to be iterated many times, which greatly weakens the optimization efficiency of the complex structure. An agent model is a model that uses mathematical methods to approximate the relationship between input and output parameters. After meeting the accuracy, it can replace FEM, reduce the workload of calculation and improve the optimization rate. Therefore, to improve the optimization efficiency of machine, different proxy models are applied to the field of optimization, such as response surface model [6], [7], neural network model [8], [9], Kriging model [10], [11], support vector machines [12], [13] and polynomial chaos expansion [14], [15] and so on. Response surface method is based on experimental design, which uses the determined polynomial function to approximate the real function relation, and has the shortcoming of low fitting accuracy in high-order complex system [16]. The generalization ability of the model established by the neural network method is poor. The support vector machine method is a complex machine learning modeling process, while the polynomial chaos expansion method is abstruse and suitable for uncertainty analysis [17]. Compared with the surrogate model above, the Kriging model has higher fitting accuracy for highly nonlinear and high-order processes, simple modeling, and less sample points required for the optimization of the machine [10].
The proposal of intelligent optimization algorithm broadens the train of thought in the field of optimization design. The intelligent solar tracking system proposed in reference [18] integrates the two-axis solar tracking system with the maximum power point tracker. The particle swarm optimization algorithm is used to optimize the parameters of the PI controller to improve the output power of solar cells. The solar powered of Stirling engine is optimized by genetic algorithm in reference [19]. Further, some researchers combine intelligent algorithm with agent model to carry out optimization design.
Literature [20] combines kriging proxy model with genetic algorithm and selects spherical function to construct Kriging surface interpolation to optimize the torque ripple of switched reluctance motor. In [21], Kriging model and genetic algorithm are combined to improve the optimization efficiency for torque ripple optimization of permanent magnet motors. However, the traditional genetic algorithm has the problems of poor local search ability. The internal magnetic field of S-AP-MF has the characteristics of highorder nonlinearity, which puts forward higher requirements for the accuracy and reliability of the optimization algorithm. The neighborhood cultivation genetic algorithm (NCGA) can solve the multi-objective optimization problem. And the neighborhood crossover mechanism in it helps to improve the local search ability and improve the accuracy of the optimal solution.
Considering the present situation of the research above, in order to improve the accuracy and speed of the optimization algorithm, taking improving the torque and suspension performance of the S-AP-MF as the optimization goal, this paper proposes a NCGA optimization method based on Kriging proxy model to realize the optimization design. In order to improve the accuracy of the proxy model, the multi-disciplinary analysis software Isight and Maxwell were used for joint simulation. The latin hypercube method was used to select samples uniformly and stratified, and the sensitivity analysis of S-AP-MF structural parameters was carried out to obtain the variables to be optimized. Combined with the kriging model, the high-precision kriging model of S-AP-MF was established, and the multi-objective optimization model was established. NCGA was used to optimize the parameters to get the best combination. Finally, the torque and suspension performance before and after optimization were compared by finite element analysis.
The core contributions of this study are summarized as follows: 1) The integration of Isight and Maxwell realizes the efficient and accurate acquisition of sample point data sets by the data interaction between the two software.
2) The sensitivity analysis of the parameters affecting the optimization goal is carried out to obtain the main variables affecting the optimization goal.
3) The kriging model with few sample points and high fitting accuracy is applied to the machine optimization problem instead of the traditional FEM, which greatly improves the optimization efficiency. 4) The optimization method of the combination of kriging model and NCGA is used to optimize the three optimization objectives globally. The neighborhood crossover mechanism in NCGA helps to elevate the local search ability and improve the accuracy of the optimization solution. The remainder of this article is organized as follows. In Section II, we provide the basic structure of S-AP-MF. In Section III, we propose an optimization method based on kriging model applicable to S-AP-MF. In Section IV, we describe the multi-objective optimization and thoroughly analyze the optimization results by finite element method in Section V. Finally, we draw conclusions in Section VI.

II. THE BASIC STRUCTURE OF THE MACHINE
The structure of S-AP-MF is shown in Figure 1. The basic structure of the machine mainly includes flywheel, rotor, stator core, magnetic isolation ring, permanent magnet and so on.
S-AP-MF is divided into A and B phases in the axial direction, and the axial length is short, which is beneficial to the integration of the machine. There are 8 torque poles and 4 suspension poles in each phase, and a magnetic isolation ring is arranged between the torque pole and the suspension pole, which weakens the electromagnetic coupling between the torque magnetic circuit and the suspension magnetic circuit. The torque pole control coils are connected in series to form torque windings, and the radial relative suspension control coils on the two suspension poles are connected in series to form suspension windings; An axially magnetized  permanent magnet is arranged between the two phases to provide permanent magnet bias magnetic flux to generate levitation force. The stator and rotor cores adopt spherical structure with a common spherical center, and when the rotor is deflected, the electromagnetic force always points to the spherical center of the rotor, thus avoiding the generation of disturbing torque. The initial structural parameters of the S-AP-MF are listed in Table 1.

III. OPTIMIZATION METHOD BASED ON KRIGING MODEL
Kriging surrogate model was first put forward by krige in the field of geostatistics, Matheron developed the concept of geostatistics and used kriging proxy model to deal with deposit preservation and error estimation [22]. In recent years, among many agent models, kriging agent model has attracted wide attention because it can accurately and effectively deal with high-dimensional nonlinear problems. The kriging agent model can effectively improve the optimization speed of the machine instead of the FEM, and its construction process is shown in Figure 2.

A. OPTIMIZATION OBJECTIVES AND SELECTION OF OPTIMIZATION VARIABLES
S-AP-MF adds a permanent magnet from the point of view of improving the suspension control precision, but its operation principle is still similar to that of the switched reluctance motor and follows the ''magnetoresistance minimum principle,'' which leads to the inherent defect of the switched reluctance motor-large torque ripple. Therefore, in order to improve the torque and suspension output and reduce the torque ripple, the torque, suspension force and torque ripple are selected as the target variables, in which the torque ripple m can be theoretically calculated by formula (1).
where, T max , T min and T avg refer to the maximum, minimum and average of the measured torque, respectively. The change of the parameters between the stator and rotor poles will affect the flux conduction path, and then affect the electromagnetic performance of the machine [23]. Therefore, this paper selects 8 design variables: stator ball radius R, rotor tooth width W 1 , suspension tooth width W 2 , torque pole arc P, torque tooth width W 3 , torque yoke high H 1 , rotor tooth height H 2 , rotor yoke height H 3 to study the influence on torque performance and suspension force.

B. DESIGN OF EXPERIMENT (DOE)
With mathematical statistics as the theoretical support, DOE arranges the experimental scheme scientifically and reasonably, which is helpful to correctly analyze the influence of input parameters on output parameters and improve the accuracy of the agent model to a certain extent. This experimental design adopts latin hypercube sampling and selects eight design variables: stator ball radius, rotor tooth width, suspension tooth width, torque pole arc, torque tooth width, torque yoke height, rotor tooth height and rotor yoke height to study the influence on machine output performance.
The initial values and ranges of design variables are listed in Table 2.  There are many parameters that affect the torque and suspension performance of S-AP-MF. From the point of ensuring the accuracy of proxy model and improving optimization efficiency, sensitivity analysis is carried out on eight design variables: stator ball radius, rotor tooth width, suspension tooth width, torque pole arc, torque tooth width, torque yoke height, rotor tooth height and rotor yoke height. Eight design variables are sampled by latin hypercube and 20 sets of data sets are obtained by simulation, and then the sensitivity analysis of 20 sets of data is shown in Figure 3.
As can be seen from Figure 3, the factors that have great influence on torque are stator ball radius, suspension tooth width, torque tooth width, rotor yoke height and torque yoke height; the factors that have great influence on torque ripple are stator ball radius, suspension tooth width, torque tooth width and torque yoke height; The factors that have great influence on suspension force are stator ball radius, suspension tooth width, torque tooth width and torque yoke height. There is little difference between the absolute values of the sensitivity of rotor yoke and torque yoke to torque, indicating that the effects of rotor yoke and torque yoke on torque are similar. Considering the sensitivity of each design variable to the target variable, the stator ball radius, suspension tooth width, torque tooth width and torque yoke height are selected as the optimization variables.

C. AGENT MODEL
The mathematical expression of the relationship between the variables and responses represented by the Kriging surrogate model is shown in (2): where, y(x) represents the response function, u(x) is the polynomial used to approximate the model, z(x) is a stationary random process related to the deviation of the regression model and the objective function, the mean is zero, and its covariance matrix is where, σ 2 represents the variance of z(x), and R(θ, x i , x j ) takes θ as the parameter to represent the spatial relationship between the two test points i and j. In this paper, the Matern Cubic function is used as the relevant function, and the function is expressed as: In the formula above, θ is the fitting correlation coefficient, and the kriging fitting effect of the three target values is shown in Figure 4.
The analysis of the fitted plots of Figure 4 is as follows: The influence of the stator ball radius on the torque pulsation is more prominent. As the stator ball radius becomes larger, the torque pulsation is effectively suppressed; The torque value shows an increasing trend when the stator ball radius and the torque pole tooth width increase in size. This is due to the fact that the inductance value becomes larger as the square area between the stator torque pole and rotor pole becomes larger; Similarly, the suspension force increases with the grow in the size of the suspension tooth width. It is worth mentioning that the suspension force increases as the stator ball radius decreases. This is probably due to the fact that the overall size of the machine becomes smaller as the radius of the stator ball decreases, and the levitation force generated by the biased permanent magnets becomes relatively larger.

D. ERROR ANALYSIS
In order to ensure the accuracy of the agent model, it is very necessary to analyze the error of the established kriging model. The usual method is to select additional sample points for error analysis, and the error can be calculated approximately instead of the original one within a reasonable range. In this paper, R 2 coefficient (multiple correlation coefficient) is used to evaluate the accuracy of the approximate model [24].
where, n is the number of additional sample points selected for error test; p is a number between 1 and n; y i is the true response value corresponding to the ith sample point; y ip is the predicted value of the proxy model corresponding to the ith sample point;ȳ is the average of the true response values of n sample points. The value of R 2 ranges from 0 to 1. When the value of R 2 is zero, the proxy model has no correlation with the actual model at all. The closer the value of R 2 is to 1, the higher the accuracy of the approximate model is. Generally, if the value is greater than 0.9, it is considered to have a high precision. As shown in Figure 5, the 45 • diagonal line in it represents that the predicted and actual values coincide completely, and the actual and predicted values of torque, torque pulsation and suspension force are basically on this diagonal line, and the R 2 value of them is 0.9865, 0.9851 and 0.917 respectively.
This error result shows that the approximation model selected in this paper has high confidence and can replace the FEM for subsequent optimization calculation.

IV. MULTI-OBJECTIVE OPTIMIZATION
In the multi-objective optimization problem, the optimal solution satisfying one objective will not satisfy other objective functions at the same time, so it is necessary to find the non dominated solution, that is, the pareto solution. The mathematical expression of multi-objective optimization problem is as follows [25]: where: x is the decision variable; D is the decision space; k is the number of target functions; F(x) refers to the objective vector consisting of objective functions f(x), and n is the objective space. Suppose x A is a feasible solution of Eq. (6), and there is no other feasible point x B to make formula (7) hold, then x A is a non-dominated solution (pareto solution).
The Kriging model built above is used to replace FEM, and NCGA is used to optimize the S-AP-MF. NCGA is developed from genetic algorithm, in which all targets have the same importance. After sorting all targets, ''adjacent propagation'' is realized by crossing, thus increasing the probability of cross propagation of solutions close to pareto frontier and accelerating the convergence process [24]. The steps of NCGA optimization algorithm are as follows: x Initialization: t = 0, set the first generation individual P 1 , the population number is N , calculate the fitness function A corresponding to the individual.
y When t = t + 1, P t = A t−1 .
z Sorting: individuals are sorted according to the direction of the target value.
{ Grouping: divided into several groups according to the order of step z, each group has two individuals.
| Crossover and variation: the two parent individuals generated in step { produce two offspring individuals, and at the same time, delete the parent individuals.
} Recombination: all progeny form a new P t . Update: N individuals were selected between A t−1 and P t based on the environment selection mechanism.
Termination: if the conditions are met, end, otherwise return to step y.
Due to the analysis above, the mathematical model of S-AP-MF's torque T , torque pulsatio m and suspended force F avg can be described as follows: (8) According to the fitting analysis of the agent model above, in the multi-objective optimization problem of torque, torque pulsation and suspension force, there is a contradiction in the parameter selection of stator ball radius. At this time, we need a compromise solution to reconcile the two objective function values, so it is necessary to apply the multi-objective optimization algorithm. Figure 6 shows the optimization iteration process. The optimal scheme is obtained after 142 iterations. At this time, the machine optimization scheme is listed in Table 3. Figure 7 shows the pareto front of the feasible solution obtained in the optimization process, and the circled feasible solution is the optimal solution selected by this optimization design. Torque, torque ripple and levitation force can not be optimal at the same time, and there is a contradiction among them, so we can only choose our optimal solution by compromise. It can be seen from Figure 7 that the torque and levitation force of S-AP-MF are ideal and the torque ripple is low in 142 iterations. Therefore, the feasible solution of 142 iterations is selected as the optimal solution of this decision.

V. ANALYSIS OF OPTIMIZATION RESULTS
According to the NCGA optimization method based on the kriging surrogate model, the optimized parameter values are obtained, and the FEM is established. The torque and suspension performance of FEM before and after  optimization are compared to verify the rationality of the proposed optimization method. Figure 8 presents a comparison of torque and suspension performance before and after optimization. The analysis results show that the optimized torque and suspension performance are significantly improved. From Figure 8(a), it can be seen that the torque curve of the optimized machine VOLUME 10, 2022  is smoother than that before optimization. The torque ripple value is dropped from 1.578 N·m, to 1.137 N·m, and the torque ripple is effectively suppressed. While the torque ripple is optimized, the torque output is also improved, and the average torque is rose from 1.699 N·m to 1.861 N·m. The average suspension force of the machine is increased from 269.67 N to 284.52 N, which is 5.51% higher than that before optimization.

VI. CONCLUSION
In this paper, the spherical axial split-phase permanent maglev flywheel machine has been taken as the research object, and the torque, torque pulsation and levitation force of the machine have been optimized based on kriging model and NCGA algorithm. The latin hypercube was used to construct the sample set, and combined with the finite element calculation, an approximate model is fitted to replace the FEM with high calculation cost. For this model, the pareto solution set was derived by NCGA algorithm to obtain the optimal solution. The optimized parameters and initial structural parameters were calculated by FEM. The analysis results show that the torque and average levitation force have been increased by 9.54% and 5.51% respectively, the torque tripple has been reduced by 27.95%, which verifies the effectiveness and reliability of the optimization method. XINYA LI was born in Taizhou, Jiangsu, China, in 1999. She received the B.Sc. degree from the Jiangsu University of Technology, Changzhou, China, in 2020. She is currently pursuing the master's degree with the Nanjing Institute of Technology. Her research interests include flywheel energy storage system for applications in electric power systems, renewable energy generation, and electric vehicle.
BINGYU CONG was born in Weihai, Shandong, China, in 1998. She received the B.Sc. degree from the School of Energy and Power Engineering, Nanjing Institute of Technology, Nanjing, China, in 2020, where she is currently pursuing the master's degree. Her research interests include flywheel energy storage system for applications in electric vehicle and urban rail transit. VOLUME 10, 2022