Multi-Objective Optimization of the Halbach Array Permanent Magnet Spherical Motor Based on Support Vector Machine

: The fundamental harmonic amplitude and waveform distortion rate of the air-gap ﬂux density directly a ﬀ ect the performance of a permanent magnet spherical motor (PMSM). Therefore, in the paper, the axial air-gap magnetic ﬁeld including the end leakage of the Halbach array PMSM is analyzed and optimized. In order to reduce the calculation time of the objective function, the air gap magnetic ﬁeld model adopts a non-linear regression model based on support vector machine (SVM). At the same time, the improved grid search (GS) algorithm is used to optimize the parameters of SVM model, which improves the e ﬃ ciency and accuracy of parameter optimization. Considering the inﬂuence of moment of inertia on the dynamic response of the motor, the moment of inertia of the PMSM is calculated. This paper takes the air gap magnetic density fundamental wave amplitude, waveform distortion rate and rotor moment of inertia as the optimization objectives. The particle swarm optimization (PSO) algorithm is used to optimize the motor structure with multiple objectives. The optimal structure design of the PMSM is selected from all of non-dominated solutions by the technique for order preference by similarity to an ideal solution (TOPSIS). The performance of the motor before and after the optimization is analyzed by the method of ﬁnite element (FEM) and experimental veriﬁcation. The results verify the e ﬀ ectiveness and e ﬃ ciency of the optimization method for the optimal structure designing of the complex PMSM.


Introduction
The structure optimization design for a motor is to further improve its performance, such as back-electromotive force (EMF) performance, and torque performance, which is largely influenced by the air-gap magnetic field. A permanent magnet spherical motor (PMSM) has a wide application in multi-degree-of-freedom motion systems such as spacecraft, manipulators, and panoramic cameras because of its small size, light weight, and high torque output [1,2]. In [3], the Halbach magnet array is applied to spherical motors, making the air-gap magnetic field closer to the sinusoidal distribution. However, due to the short axial length of PMSM, the axial air-gap magnetic density shows a rapid attenuation from the equator to the two ends of the permanent magnet (PM), which results in the distortion of the main air-gap magnetic field [4]. Moreover, the amplitude and waveform of the air-gap magnetic field affect the magnitude and fluctuation of electromagnetic torque, which directly influences the performance of the motor. Therefore, it is necessary to analyze and optimize the air-gap magnetic field of the PMSM.
The modeling of the magnetic field is the premise of the optimization design of the motor. Various methods have been proposed to formulate the magnetic field distribution. The finite element method (FEM) is employed to optimize the rotor structure of the PMSM in [5]. A rough optimization result is obtained based on a given set of specifications. In [6,7], an analytical model of the air-gap magnetic field of the PMSM is derived by solving the Laplace equation. However, in the modeling process, the approximate equivalent of the shape of the permanent magnet pole is adopted, which leads to low modeling accuracy. A three-dimensional (3D) magnetic equivalent circuit model of PMSM is established in [8]. However, the approximated magnetic circuit and neglected end effect influences accuracy to some extent. In [6], a novel magnetic field method is proposed for spherical actuator with cylindrical shaped PM poles based on geometrical equivalence principle. On the basis of optimal equivalent principle, an accurate harmonic model is set up, and the magnetic field is formulated. However, the shape of PM pole is transformed into an approximated dihedral cone. The accuracy of the model is then reduced with the simplified analysis. Considering time cost and generality, the conventional modeling methods are not suitable for the parameters optimization which involves large-scale iterative calculation. Compared with conventional modeling methods, the support vector machine (SVM) modeling method does not require the users to have a deep understanding of the intrinsic mechanism of the motor; instead, it only needs to establish the mapping relationship between input and output [9]. Moreover, only a limited amount of sample data is required to obtain a more precise model, which greatly improves the efficiency of optimization [10]. In [11,12], a non-parametric model of the claw-pole alternator is established, and the flux leakage is optimized by SVM. The results of FEM verify the effectiveness and efficiency of SVM. In [13,14], an SVM model for the electromagnetic torque of the PMSM is established. Based on this model, the stator and rotor structures of the motor are optimized.
Multi-objective optimization of the motor can be achieved by using optimization algorithm. Among the developed optimization algorithms, particle swarm optimization (PSO) is widely applied in the optimization of the motor relying on its high convergence speed and simple algorithm [15]. In [16], the multi-objective PSO optimization algorithm is adopted to optimize the amplitude and waveform distortion of the permanent magnet synchronous motor air-gap flux density. In [17], based on the derived analytical expression of the air-gap magnetic field for the PMSM, the waveform distortion of air-gap magnetic field is optimized by PSO.
In this paper, the axial air-gap magnetic field of the Halbach array PMSM is analyzed. The non-parametric model of axial air-gap magnetic field is established based on SVM. An improved Grid Search (GS) algorithm is proposed to optimize the parameters of SVM model. Since the PMSM is mostly applied in the servo system that requires high control accuracy, the moment of inertia has a great influence on the response speed of the motor. Therefore, the moment of inertia is considered in this paper. Finally, a multi-objective optimization problem with three design variables and three objectives is designed by combining SVM and PSO. The optimal design of the PMSM is selected by the technique for order preference by similarity to an ideal solution (TOPSIS).

Harmonic Analysis Based on FEM
The 3D FEM model of the Halbach array PMSM to be optimized is shown in Figure 1a. A total of 12 (four poles and three blocks for each pole) NdFeB PMs with approximated dihedral cone structure are embedded in the rotor equator. The Halbach array PM magnetizing direction on the central plane is shown in Figure 1b. The stator coils are evenly embedded in the stator shell in three layers, with 12 in each layer. Table 1 shows the basic parameters of the motor.  Based on the FEM model shown in Figure 1a, adaptive mesh subdivision is employed. The boundaries between the PM, air gap, and winding meet the natural boundary condition. The outermost boundary region satisfies the first kind of homogeneous boundary conditions. Figure 2 shows the variation of radial component of the axial air-gap magnetic flux density (B 1r ) with θ, when ϕ = 0 • , ϕ = 15 • , and ϕ = 30 • , respectively. Figure 3 shows that B 1r has the same trend with θ under different paths, which is rapidly attenuated from the equator to the end of the PM, thus leading to the distortion of the main magnetic field inside the motor and the reduced the peak value of back-EMF [18]. Therefore, the axial air-gap magnetic field should be optimized to improve the performance of the motor. After the materials of the PM and iron core are determined, the main factors influencing the air-gap magnetic field of the motor are the structural parameters of the PM and the air-gap length of the motor. Therefore, in this paper, the thickness of PM (h), the latitude angle of PM (θ 1 ), and the air-gap length (g) are selected as the optimization variables. The influence of the optimization variables on the fundamental harmonic amplitude and waveform distortion of the axial air-gap magnetic field is investigated. In this paper, the waveform distortion of air-gap magnetic field is defined as the ratio of the root mean square of the higher harmonic amplitude to the fundamental harmonic amplitude.
1r is the fundamental harmonic amplitude and B ν 1r is the higher harmonic amplitude. The harmonics amplitudes are obtained by Fourier transform.
Keeping the other parameters unchanged, h is changed from 11.1 mm to 25.1 mm. Figure 3a shows the changing trend of the fundamental wave amplitude and waveform distortion rate of the axial air gap magnetic density with h. It is shown that the fundamental harmonic amplitude increases with the increase of h, but the increasing amplitude gradually gets smaller. The waveform distortion decreases with the increase of h. Considering the performance of the motor and the difficulty of magnet processing, h is selected within the range of 17.1-25.1 mm.
Keeping the other parameters unchanged, θ 1 is changed from 40 • to 160 • . Figure 3b shows the changing trend of the fundamental wave amplitude and waveform distortion rate of the air gap magnetic density with θ 1 . In Figure 3b, with the increase of θ 1 , the fundamental harmonic amplitude first increases and then decreases, while the waveform distortion first decreases and then increases. Considering the performance of the motor and the fabrication difficulty of magnet, θ 1 is selected within the range of 80 • -120 • .
Keeping the other parameters unchanged, g is changed from 0.2 mm to 2.0 mm. Figure 3c shows the changing trend of the fundamental wave amplitude and waveform distortion rate of the air gap magnetic density with g. Figure 3c shows that both the fundamental harmonic amplitude and waveform distortion decrease as g increases. Considering the performance of the motor, g is selected with the range of 0.4-1.6 mm.
The above harmonic analysis of the air gap magnetic field determines the reasonable selection range of design variables and lays a foundation for further optimization in the future.

Modeling of SVM and Calculation of Moment of Inertia
The main idea of SVM regression is to map the training sample to high-dimensional space through a nonlinear relationship. The key to establishing the SVM model is how to construct the relationship between input and output. In this paper, the radial basis function is chosen to construct the SVM model.

Support Vector Machine Regression Modeling Method
Let the training sample be (x i , y i ) | i = 1, 2 . . . , n , where n is the number of training samples, x i ∈ R N is the input of the i-th sample, and y i ∈ R is the output of the i-th sample. The sample space R N is mapped to a high-dimensional feature space through the nonlinear function ϕ(x i ), and then linear regression is performed. The regression function is where w is the weight and b is the offset. It is assumed that the linear regression function can accurately fit the sample space data under the accuracy ε. Then, Then, the regression problem can be reduced to the following optimal one.
Then, the following formula can be obtained: where w is a term related to the complexity of the function f, ε is the insensitive loss function, ξ i , ξ * i is the relaxation factor, and C is the penalty factor.
Lagrange multipliers are used to solve the convex quadratic naturalization problem, the results are shown as follows: and the constraint conditions are satisfied.
where α i , α * i represents the Lagrange multiplier. Then, the available weight w is SV stands for support vector. Then, the regression function of the support vector machine can be expressed as For nonlinear regression problems, the sample space can be transformed into a high-dimensional space through nonlinear transformation, and linear regression analysis can be performed in the high-dimensional space. In a high-dimensional space, the kernel function K(x i , x) is used to replace the dot product (x i , x) in the above formula, so that

The Choice of Kernel Function and the Influence of its Parameters on Regression
The form of the kernel function and the value of its parameters affect the complexity and accuracy of the SVM regression model. SVM nonlinear regression modeling is to map a nonlinear sample space to a high-dimensional space through a kernel function, and then perform linear regression. Only with an appropriate choice of kernel function type and corresponding kernel parameters can an SVM regression model with good generalization ability be constructed. This paper chooses the RBF kernel function to construct the SVM regression model, and its expression is The kernel function parameter σ and penalty factor C in the SVM model have an important impact on the establishment of an accurate SVM model. The penalty factor C is used to measure the tolerance to training errors, and the kernel function parameter σ determines the distribution of the sample space mapped to the high-dimensional space. When σ is small, the training error is small. At this time, almost all of the sample data are support vectors, and the generalization ability of the SVM is poor. When σ gets larger, the support vector decreases, and the corresponding training error and test error increase. The generalization ability of the support vector machine is still relatively poor. When the value of the penalty factor C is smaller, the penalty for experience error is smaller. The less complex the learning machine, the greater the experience risk, and vice versa.

Establishment of Sample Space
Two sets of sample dots by FEM are built for the training and testing of SVM model. Table 2 shows the data distribution of the training sample space. From Table 2, 5 3 sample dots need to be established. It is time consuming to build 125 models by FEM. Here, the Taguchi method based on orthogonal experiment is adopted [19]. The orthogonal experiment mainly relies on the number and level of the studied parameters. The advantage is that fewer experiments need to be conducted to obtain the best combination of design parameters. According to Table 2, the orthogonal experiment L 25 (5 3 ) with three parameters and five levels is established. The corresponding parameter design and the finite element results of two objective functions are shown in Table 3. Twenty-five sample dots are used to build the SVM model, thus reducing the computational burden.

Parameter Optimization of SVM Model
The kernel function parameter σ and the penalty factor C in the SVM model have an important influence on the accuracy of SVM model [20]. The penalty factor C is used to measure the tolerance to training errors, and the kernel function parameter σ determines the distribution of sample space mapped to high dimensional space. In this paper, a GS algorithm is used to obtain the optimal SVM model parameters. The accuracy of the SVM regression model is measured by the root Mean Square Error (MSE), and its expression is as follows.
where f FEM and f SVM are the values from FEM and SVM, respectively. Figure 4 shows the 3D diagram of the parameter optimization of the fundamental harmonic amplitude and waveform distortion of air-gap magnetic flux density. C and σ are within the range of (0,100), and the step is 0.1.
Energies 2020, 13, x FOR PEER REVIEW 8 of 20 According to Table 2, the orthogonal experiment L 25 (5 ) with three parameters and five levels is established. The corresponding parameter design and the finite element results of two objective functions are shown in Table 3. Twenty-five sample dots are used to build the SVM model, thus reducing the computational burden.

Parameter Optimization of SVM Model
The kernel function parameter σ and the penalty factor C in the SVM model have an important influence on the accuracy of SVM model [20]. The penalty factor C is used to measure the tolerance to training errors, and the kernel function parameter σ determines the distribution of sample space mapped to high dimensional space. In this paper, a GS algorithm is used to obtain the optimal SVM model parameters. The accuracy of the SVM regression model is measured by the root Mean Square Error (MSE), and its expression is as follows.
where and are the values from FEM and SVM, respectively. Figure 4 shows the 3D diagram of the parameter optimization of the fundamental harmonic amplitude and waveform distortion of air-gap magnetic flux density. C and σ are within the range of (0,100), and the step is 0.1.  It can be seen from Figure 4 that, within a certain range, the MSE corresponding to C and σ is very small, while in most cases, the MSE is very large. Under normal circumstances, the optimal solution can be found using a large-scale and small-step grid, which, however, is often time-consuming. Therefore, an improved GS algorithm is proposed in this paper. First, the rough search is performed with large steps within a large range, and then the search space and step are gradually narrowed for detailed search. The detailed process is as follows: 1. The initial search space is (0,100) and the initial step is 1. The rough search is performed to obtain an optimal parameter C and σ. 2. Based on the obtained optimal parameter C and σ, the search space is halved, and a more accurate search is performed with a small step. 3.
Step 2 is repeated until the step size is reduced to 0.1, and with the search stopped, the optimal solution is then output.
In order to verify the effectiveness of the improved grid optimization algorithm, Table 4 and Table 5 show the comparison results of the optimization time and prediction error of the air gap magnetic density fundamental wave amplitude and waveform distortion rate obtained using the improved GS algorithm, the traditional GS algorithm, the PSO algorithm, and the genetic algorithm (GA), respectively.  It can be seen from Figure 4 that, within a certain range, the MSE corresponding to C and σ is very small, while in most cases, the MSE is very large. Under normal circumstances, the optimal solution can be found using a large-scale and small-step grid, which, however, is often time-consuming. Therefore, an improved GS algorithm is proposed in this paper. First, the rough search is performed with large steps within a large range, and then the search space and step are gradually narrowed for detailed search. The detailed process is as follows: 1.
The initial search space is (0,100) and the initial step is 1. The rough search is performed to obtain an optimal parameter C and σ.

2.
Based on the obtained optimal parameter C and σ, the search space is halved, and a more accurate search is performed with a small step. 3.
Step 2 is repeated until the step size is reduced to 0.1, and with the search stopped, the optimal solution is then output.
In order to verify the effectiveness of the improved grid optimization algorithm, Tables 4 and 5 show the comparison results of the optimization time and prediction error of the air gap magnetic density fundamental wave amplitude and waveform distortion rate obtained using the improved GS algorithm, the traditional GS algorithm, the PSO algorithm, and the genetic algorithm (GA), respectively.  According to Tables 4 and 5, compared with traditional GS algorithms, the improved GS algorithm shows advantages including less optimization time and lower prediction error. Moreover, compared with the optimization results of PSO algorithm and GA algorithm, the improved GS algorithm is also better in optimization time and prediction error. The result shows that the improved GS algorithm improves the efficiency and accuracy of parameter optimization.
Based on the improve GS algorithm, the prediction errors of the fundamental harmonic amplitude and waveform distortion are 1.6382% and 1.3278%, respectively. This shows that the SVM model has good generalization capabilities. In addition, for the same sample data, much less time is used to calculate the output of SVM model than that of the FEM. Therefore, the SVM model can be used for the online calculation of parameter optimization of the PMSM [14].

Calculation of Moment of Inertia
PMSM is mostly applied in servo systems, which have a high requirement on the dynamic response and control accuracy of the motor. Therefore, in the optimization design of the PMSM, the moment of inertia of the rotor should be reduced. The expression of moment of inertia is where J is the moment of inertia of the rigid body, m i is the mass of the particle, and r i is the distance from the particle to the axis of rotation.
The rotor of the PMSM contains the PMs and the rotor yoke. According to (3), the expression of the moment of inertia of the PMSM can be obtained as where J z is the moment of inertia with the z axis as the axis of rotation, R m is the outer radius of the PM, R r is the inner radius of the PM, ρ 1 is the mass density of the rotor yoke (ρ 1 = 1.02 g/cm 3 ), and ρ 2 is the mass density of the PM (ρ 2 = 7.5 g/cm 3 ).

Optimization Theory
In this paper, there are three optimization goals, which belong to the multi-objective optimization problem. In multi-objective optimization, it is difficult to find a true optimal solution due to the conflicts between objectives. There are a series of solutions, which are characterized by the existence of at least one solution whose objective is superior to other objectives. Such solutions are called non-dominated solutions.
Based on the related mathematical model of the Halbach array permanent magnet spherical motor established in the previous chapter, the linear weighted summation method is used to transform the multi-objective optimization problem into a single-objective one. By assigning different weights to different objective functions, a set of non-inferior solutions are obtained. The ordinal preference method that approximates the ideal solution is adopted to find the most satisfactory result from the non-inferior solutions, which was developed by Hwang and Yoon [21]. The basic theory is based on the smallest difference between the selected plan and the ideal plan, and the largest difference between the selected plan and the negative ideal plan.
Supposing that there are n attributes and m solutions in an optimization problem, which correspond to a geometric system of m points in a n-dimensional space, the coordinates of each point are determined by the standardized weighted attribute values of each solution. Within comparison with the positive ideal solution and the negative ideal solution, the advantages and disadvantages of each scheme are judged. The specific process is as follows: 1.
The normalized decision matrix is obtained by normalization method. Suppose that the decision matrix of the multi-attribute decision-making problem is Y = y ij , and the standardized decision matrix is Z = z ij , then 2.
Construct a weighted normalization matrix. According to the weight coefficient given by the designer, multiply it by the standardized decision matrix x ij = w j z ij (16) where w j is the weight of the j-th objective function.

3.
Determine positive and negative ideal solutions. Positive ideal solution: Negative ideal solution: x − j = max x ij j j 2 min x ij j j 1 (18) where j 1 is benefit type attribute, and j 2 is cost type attribute. x + j is the j-th attribute value of positive ideal solution x + , and x − j is the j-th attribute value of negative ideal solution x − . 4.
Calculate the distance from each scheme to the positive ideal and negative ideal solution under different weight coefficients. The distance from the objective function to the positive ideal solution x + j is d + j , and the distance to the negative ideal solution

5.
Calculate the closeness of the solution.
where 0≤ C * i ≤ 1, When C * i = 0, it means that the scheme is the worst scheme; when C * i = 1, it means that the scheme is the best one. 6.
Sort the according to the calculated size of C * i . The larger the C * i , the better the solution.

Multi Objective Optimization of PMSM
For the multi-objective optimization of the PMSM, the weight coefficient method is employed to transform multi-objective optimization into single objective optimization. In this paper, the aim of optimization is to obtain the relatively maximum fundamental harmonic amplitude, and the relatively minimum waveform distortion and the moment of inertia. So, the fitness function F can be designed as: where F 10 and F 20 are the reference values of fundamental wave amplitude and waveform distortion rate, respectively. When h = 19.1 mm, θ 1 = 60 • , and g = 1 mm, F 10 and F 20 are 0.3170 T, 0.6519, respectively. Q 1 , Q 2 , and Q 3 are the weights corresponding to the fundamental wave amplitude, waveform distortion rate, and moment of inertia. f 1 represents the function of fundamental harmonic amplitude and f 2 represents the function of waveform distortion, both of which are the functions of h, g, θ 1 . f 3 is the normalized expression of moment of inertia, and as the function of h, θ 1 , it can be expressed as where J 0 is the moment of inertia when h = 17.1mm and θ 1 = 80 • , J 1 is the moment of inertia when h = 25.1 mm and θ 1 = 120 • , and J z is the moment of inertia with the z axis as the axis of rotation. From the above analysis, the multi-objective optimization problem in this paper can be expressed as In the optimization design of this article, the population number of the particle swarm algorithm is set as 50, the number of iterations as 1000, and the learning factor as 1.4944. Different weighting coefficient ratios will produce different optimization results, and each weighting ratio is determined according to the researcher's preference or engineering needs. In this paper, 36 sets of weighting ratios are randomly generated in the form of grid division, the sum of weight coefficients is 1, and 36 sets of weighting ratios are generated in the interval 0-1 with a step of 0.1, that is, 36 optimization schemes. Table 6 shows the optimization results obtained by using the particle swarm algorithm under different weight coefficients (the optimization results are the normalized ones). Figure 6 is a three-dimensional graph drawn based on the optimization results in Table 6. Table 6. Optimization results under different weight coefficients.  However, it is difficult to find the optimal solution directly from it. The mathematical evaluation of the non-dominated solutions is implemented based on the ordinal preference method which approximates the ideal solution. Figure 6 shows the basic flow chart of the technique. According to the distribution of feasible solutions in Figure 6 and the evaluation method in Figure 5, Figure 7 gives the evaluation results of the feasible solutions. It can be seen from Figure 7 that the 15th group of feasible solutions has the largest closeness coefficient. Therefore, this solution is the best one. According to the evaluation results, the corresponding optimization result is h = 23.3484 mm, θ 1 = 113.46 • , and g = 0.8951 mm. Considering the fabrication of the PM, h = 23.3 mm, θ 1 = 113 • , and g = 0.9 mm are selected.

Verification
The structural parameters of the Halbach array permanent magnet spherical motor before and after optimization are shown in Table 7.  Figure 8 shows the comparison of non-optimized and optimized axial air-gap flux density when ϕ = 0 • .  Table 7 and Figure 8, the waveform distortion is obviously reduced, and the fundamental harmonic amplitude is greatly improved for the optimized PMSM.
Under no-load conditions, the rotation speed is 30 rpm. Figure 9a shows the comparison of the influence on the back EMF of the permanent magnet spherical motor with and without the end effect considered before the optimization. Figure 9b shows the comparison of the influence on the back EMF of the permanent magnet spherical motor with and without the end effect considered after the optimization.  It is shown in Figure 9 that the peak value of the spinning back-EMF with the end effect considered is significantly smaller than that without the end effect considered, and the optimized PMSM has a smaller difference between the peak values with and without the end effect considered. Here, a coefficient S is defined to measure the influence of the end effect on the spinning back-EMF.
where P no end−e f f ect is the peak value of the spinning back-EMF without the end effect considered and P end−e f f ect is the peak value of the spinning back-EMF with the end effect considered. The smaller the ratio, the smaller effect of the end effect on the spinning back-EMF. The coefficients corresponding to Figure 9a,b are 0.3499 and 0.2449, respectively. It can be obtained that the end effect has an influence on the peak value of the spinning back-EMF. For the optimized PMSM, the influence of the end effect on the spinning back-EMF is weakened. Figure 10 shows the variation of the spinning torque and tilting torque generated by a single coil. From Figure 10a, the optimized spinning torque has a larger peak value. The Fourier decomposition of the torque is carried out. Before the optimization, the fundamental amplitude and waveform distortion rate of the rotation torque are 113.6380 mT and 0.0989, respectively. After the optimization, the fundamental wave amplitude and waveform distortion rate of the rotation torque are 143.7316 mT and 0.0950, respectively. It is shown in Figure 9 that the peak value of the spinning back-EMF with the end effect considered is significantly smaller than that without the end effect considered, and the optimized PMSM has a smaller difference between the peak values with and without the end effect considered. Here, a coefficient S is defined to measure the influence of the end effect on the spinning back-EMF. where is the peak value of the spinning back-EMF without the end effect considered and is the peak value of the spinning back-EMF with the end effect considered. The smaller the ratio, the smaller effect of the end effect on the spinning back-EMF. The coefficients corresponding to Figure 9a,b are 0.3499 and 0.2449, respectively. It can be obtained that the end effect has an influence on the peak value of the spinning back-EMF. For the optimized PMSM, the influence of the end effect on the spinning back-EMF is weakened. Figure 10 shows the variation of the spinning torque and tilting torque generated by a single coil. From Figure 10a, the optimized spinning torque has a larger peak value. The Fourier decomposition of the torque is carried out. Before the optimization, the fundamental amplitude and waveform distortion rate of the rotation torque are 113.6380 mT and 0.0989, respectively. After the optimization, the fundamental wave amplitude and waveform distortion rate of the rotation torque are 143.7316 mT and 0.0950, respectively.
(a) The result shows that the optimized torque characteristic is improved, which is helpful in improving the torque output and operation stability of the PMSM.

Verification by Experiment
A prototype shown in Figure 11 is manufactured. Figure 11a is the structure of Halbach array PM. Figure 11c is the assembled Halbach array PM of spherical motor rotor. The experimental device for detecting the magnetic field is shown in the Figure 12. The Fourier decomposition of the torque is carried out. Before the optimization, the fundamental amplitude and waveform distortion rate of the tilting torque are 54.6160 mT and 0.9851 respectively. After the optimization, the fundamental wave amplitude and waveform distortion rate of the tilting torque are 73.1590 mT and 0.4810, respectively.
The result shows that the optimized torque characteristic is improved, which is helpful in improving the torque output and operation stability of the PMSM.

Verification by Experiment
A prototype shown in Figure 11 is manufactured. Figure 11a is the structure of Halbach array PM. Figure 11c is the assembled Halbach array PM of spherical motor rotor. The experimental device for detecting the magnetic field is shown in the Figure 12.  Figure 12 shows the optimized experimental device diagram, and the specific dimensions are given in Table 7. The variation of axial air gap flux density with latitude can be detected by making the output shaft yaw in the guide rail. The axial air gap magnetic field is detected by a fixed magnetic field sensor and its value is displayed in a Gauss meter. In this experiment, the reliability of the data needs to be ensured so as to reduce the error. The measurement is conducted every 5 • , and three measurements are taken on the same track and the average value is used. The comparison between the experimental and analytical results of the optimized axial air gap flux density is shown in Figure 13. It can be seen from Figure 13 that there is good agreement between the obtained axial air gap flux density model and the experiment results in waveform distribution and numerical values.
The correctness of the magnetic field model established in this paper and the effectiveness of multi-objective optimization of axial air gap magnetic field are then verified.

Conclusions
In this paper, the axial air-gap magnetic field of the Halbach array PMSM is modeled based on SVM. The reasonable selection range of the design variables is determined through FEM. An improved GS algorithm is proposed to optimize the model parameter. At the same time, the influence of moment of inertia on the response speed of the motor is considered. Based on the established mathematical model, the optimal structure of the PMSM is obtained using PSO. The results show that the axial air-gap magnetic field of the optimized PMSM is improved. The influence of the end effect on the peak value of the back-EMF is weakened, and the performance of the electromagnetic torque is improved.

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