Improved Method for Distributed Parameter Model of Solenoid Valve Based on Kriging Basis Function Predictive Identification Program

In this paper, a method for the improvement of the calculation accuracy of the distributed parameter model (DPM) of electromagnetic devices is proposed based on the kriging basis function predictive identification program (PIP). Kriging is mainly an optimal interpolation method which uses spatial self-covariance, and takes a polynomial as the basis function. The accuracy of the krigingbased surrogate model can be improved by adjusting the related functions and hyperparameters. Based on the DPM of a solenoid valve, there exist certain errors in the estimation. They can be summarized as follows: Firstly, the estimation error of magnetic flux leakage (MFL) permeance is caused directly by the deviation of the magnetic flux tube due to the segmented magnetic field line. Secondly, the estimation error of soft magnetic resistance because of the nonlinearity of the permeability of soft magnetic material leads to the change of soft magnetic resistance alongside the magnetic flux. In this paper, an improved kriging error correction method is applied to modify the leak permeance and soft magnetic resistance calculation. The kriging basis function is adjusted to adapt to the data curve of the MFL permeance error data. The calculated MFL permeance data are compared with the error variation data to select the appropriate basis function. To improve the computational efficiency, the PIP is proposed to select the appropriate basis function. The modified MFL permeance data and soft magnetic resistance are substituted into the DPM for improving the computational accuracy and efficiency of the solenoid valve.


Introduction
The optimization design of modern solenoid valves faces more and more challenges with increasing performance requirements to realize higher latching force, higher output force, and more extensive robustness [1][2][3][4][5][6][7][8][9]. As a means of simulation, the FEM method has the absolute advantage of accuracy, and the disadvantage of efficiency. The distributed parameter model (DPM) method can effectively speed up the calculation, and kriging models can decrease error. The DPM and kriging models have become important factors affecting the design efficiency and accuracy of solenoid valves. Therefore, evaluating the applicability of the kriging model has an important theoretical significance and a practical engineering value for product design optimization [10][11][12][13][14][15].
With the increasing levels of computational technology and theoretical analysis for the DPM, the kriging model, whose inputs affect its outputs, plays an important role in engineering design and error correction [16][17][18][19][20][21]. The calculation accuracy of air leakage permeance and the soft magnetic resistance has an effect on the electromagnetic force calculation accuracy of solenoid valves. To reduce the relative error and improve the accuracy of target prediction based on sample data, various data fitting methods based on the kriging model have been proposed. In the design and optimization of products, kriging has a good performance in error correction [22][23][24][25][26][27][28]. An adaptive finite difference weighted essentially non-oscillatory (WENO) method with Gauss-kriging reconstruction is proposed to reduce dissipation in smooth regions of flow, while preserving high resolution around discontinuities for hyperbolic systems of conservation laws [29]. The work in [30] proposed a kriging model-based method for the reliability design and optimization of planetary gear, using a genetic algorithm. The kriging model was used to establish the gear reliability model to simplify the reliability calculation in this method [30]. Shin et al. performed an optimization study of response time to improve the dynamic performance of a direct-acting solenoid valve based on a kriging model [31]. A kriging metamodel-based multi-objective optimization strategy has been employed to optimize the valve-plate shape of the axial piston pump [32]. The kriging method was also used to build a surrogate model, which presents the relationship between dynamic responses and dynamic simulation of the valve train [33]. Many of the works where the kriging model has been adopted for error calculation mainly focus on the calculation errors associated with replacing the original analytical model with the kriging-based surrogate model.
In this paper, a method for the improvement of the calculation efficiency and accuracy of the DPM of electromagnetic devices is proposed based on the kriging basis function predictive identification program (PIP). To obtain an improved DPM performance, an appropriate basis function is selected by contrasting various basis functions with error item curves of air leakage permeance. To further improve the calculation efficiency, the PIP is introduced to prejudge the error item curves against standard functions. The modified leakage permeance and the soft magnetic resistance data are considered in the DPM of the electromagnetic device to calculate the electromagnetic force, thereby verifying the method proposed in this paper.

The Error Correction Method of DPM Based on Kriging Model
The kriging model is a global surrogate model derived from geostatistics, and it can mathematically provide an optimal linear unbiased estimation. This makes it an appropriate method to achieve the primary goal of this paper. To estimate the linear, unbiased and minimum estimated variance using kriging, the geometric features of information sample shape, size and space on samples are taken into account. Nonetheless, kriging is still a smooth interpolation method [34,35]. The electromagnetic calculation software is Matlab in this paper.
The kriging model can be summarized in two parts, which are represented by where f j (x) is the basis, β j is the regression coefficient, Z(x) is the stochastic distribution function.
A solenoid valve with a fan-shaped permanent magnet (PM) is selected as the example in this study, and its electromagnetic system is shown in Figure 1. The armature is a moving part, the angle of the fan-shaped PM is 20 • , its inner radius is 8.5 mm, its outer radius is 17.5 mm and its height is 8 mm. The solenoid valve rated coil ampere-turns equal 1800 AT. The PM is divided into seven segments. The magnetic flux leakage (MFL) permeances of the flux tube from the outer radius side to the inner radius side are G 1 , G 2 , G 3 , G 4 , G 5 , G 6 and G 7 [36][37][38][39][40].
The solenoid valve is excited in 0 AT, and its finite element method (FEM) responses are compared with the unmodified DPM, as shown in Figure 2. There are 20 sample points for every curve.   Figure 1a).
The solenoid valve is excited in 0 AT, and its finite element method (FEM) responses are compared with the unmodified DPM, as shown in Figure 2. There are 20 sample points for every curve.
Armature stroke(mm) MFL Permeance (H) The FEM and DPM data are subtracted to obtain the MFL permeance error terms (i.e., G1, G2, G3, G4, G5, G6 and G7), which are then stored in the error matrix Er. The MFL permeance data are shown in Figure 3 which reveal a non-monotonic variation.
The solenoid valve is excited in 0 AT, and its finite element method (FEM) responses are compared with the unmodified DPM, as shown in Figure 2. There are 20 sample points for every curve. The FEM and DPM data are subtracted to obtain the MFL permeance error terms (i.e., G1, G2, G3, G4, G5, G6 and G7), which are then stored in the error matrix Er. The MFL permeance data are shown in Figure 3 which reveal a non-monotonic variation.   The FEM and DPM data are subtracted to obtain the MFL permeance error terms (i.e., ∆G 1 , ∆G 2 , ∆G 3 , ∆G 4 , ∆G 5 , ∆G 6 and ∆G 7 ), which are then stored in the error matrix Er. The MFL permeance data are shown in Figure 3 which reveal a non-monotonic variation.   Figure 1a).
The solenoid valve is excited in 0 AT, and its finite element method (FEM) responses are compared with the unmodified DPM, as shown in Figure 2. There are 20 sample points for every curve.
Armature stroke(mm) MFL Permeance (H) The FEM and DPM data are subtracted to obtain the MFL permeance error terms (i.e., G1, G2, G3, G4, G5, G6 and G7), which are then stored in the error matrix Er. The MFL permeance data are shown in Figure 3 which reveal a non-monotonic variation.
Armature stroke(mm) MFL Permeance error(H) Polynomial Function Exponential Function Gaussian Function where |d i | is the distance between the known and predicted quantities. The basis function output characteristics can be determined using the key nodes and the output of the model. Given a function k(x 1 , x 2 . . . x n ) and an input parameter x 1 , x 2 . . . x n , the domain of the function can be expressed as follows: where ∆F is the function output, ∆ 2 (∆x) is the higher derivative of the input parameter fluctuations. Given that a function relation exists between k(x 1 , x 2 . . . x n ) and ∆x 1 , ∆x 2 . . . ∆x n , k(x 1 , x 2 . . . x n ), this function relation can be described as follows: where w is the weight coefficient. For every node in the count, w = 1. However, if the sample points are not key nodes, then 0 < w < 1. Given that, the boundary condition of the input parameter is as follows: where x il , x ih are the top and bottom limitations of the ith input value. Combining Equations (7)-(9), ∆F i can be defined as: Actuators 2021, 10, 10 (7) and (8), the following equations are obtained: where d i = |x − x i |, x is the dependent variable, θ i is the PM segment positions with respect to the armature displacement. The armature stroke is 2 mm, and θ i can be obtained by particle swarm optimization (PSO) based on 20 nodes. Based on the MFL permeance error data, and adopting a periodic function (i.e., Fourier series) for ∆G 3 fitting, it can be written as follows: where x is the armature stroke, f (x) is the dependent variable of the MFL permeance.  (13), df (x)/dx has four solutions defined in the range 0 to 2. By adopting a polynomial function for ∆G 3 fitting, it can be written as follows: where Taking the derivative of Equation (14), i.e., df (x)/dx, yields four solutions. By adopting an exponential function for ∆G 3 fitting, it can be written as follows: where a = −7.591 × 10 −7 , b = 4.862, c = 7.583 × 10 −7 , d = 4.863. Taking the derivative of Equation (15), i.e., df (x)/dx, has only one solution. It is obvious that the exponent is unfit for ∆G 3 fitting. By adopting a Gaussian function for ∆G 3 fitting, it can be written as follows: where Taking the derivative of Equation (16), i.e., df (x)/dx, yields four solutions. Figure 4 shows the contrasting condition between basis functions and MFL permeance error ∆G 3 . It is obvious that there is a large error between the exponential function-based fitted curve and ∆G 3 . The mean error for this particular case is over 100%. The polynomial function-based fitted curve, Fourier function-based fitted curve and Gaussian function-  To further improve computational efficiency, the reliability index of samples is calculated by four algorithms to achieve optimal basis function. The algorithms are: The relative error of the basis functions can be obtained by taking MC as the base. Firstly, the failure probability of the samples is calculated, then the reliability index can be obtained using the calculation results. The reliability index is calculated by the MC method and it is used as the standard with a termination condition of 10 −7 . Table 1 shows the reliability index of ΔG3.  Table 1 shows that Gaussian function with PSO has the smallest count for ΔG3, while the Fourier function with PSO has the smallest count for ΔG1.
ΔG3 is presented as an example. The polynomial function is adopted for the armature stroke, Gaussian error data and ΔG3. This function can be written as follows: ΔG p p x p y p x p xy To further improve computational efficiency, the reliability index of samples is calculated by four algorithms to achieve optimal basis function. The algorithms are:
Monte Carlo (MC) The relative error of the basis functions can be obtained by taking MC as the base. Firstly, the failure probability of the samples is calculated, then the reliability index can be obtained using the calculation results. The reliability index is calculated by the MC method and it is used as the standard with a termination condition of 10 −7 . Table 1 shows the reliability index of ∆G 3 .  Table 1 shows that Gaussian function with PSO has the smallest count for ∆G 3 , while the Fourier function with PSO has the smallest count for ∆G 1 .
∆G 3 is presented as an example. The polynomial function is adopted for the armature stroke, Gaussian error data and ∆G 3 . This function can be written as follows: where p 00 = −2.057 × 10 −7 , p 10 = 7.069 × 10 −7 , p 01 = −0.02261, p 20 = −4.141 × 10 −7 , p 11 = 0.3128. x is the armature stroke, y is the Gaussian error data. Figure 5 shows Gaussian error data fitting. where p00 = −2.057 × 10 −7 , p10 = 7.069 × 10 −7 , p01 = −0.02261, p20 = −4.141 × 10 −7 , p11 = 0.3128. x is the armature stroke, y is the Gaussian error data. Figure 5 shows Gaussian error data fitting. The MFL permeance of the solenoid valve G3 can be improved. Now, assuming that P is the FEM data matrix and Q is the DPM data matrix, Er can be written as follows: If ΔG depicts the permeance error of the sampling point, Er can be derived as follows: where n is the PM segment number, m is the sampling point position. i = 1, 2, …n. j = 1, 2, …m. The basis function can be defined as fk(x1, x2…xq), k = 1, 2, 3, …n, where q is the dimension. fk(x1, x2…xq). The function can be written as: Based on Equation (20), the polynomial function is adopted for the armature stroke, l, ΔG and f(x1, x2…xq). This function is written as follows: Equation (21) can also be written as follows: , , ,..., , , , ,..., where fer(l, x1, x2…xq) is the penalty function for the error data of the sampling point with armature displacement. The ith MFL permeance or soft magnetic resistance can be defined as qi(x), where x is the unknown quantity. Therefore, wi(x) can be obtained as follows: where wi(x) is the improved DPM data, and it is very close to the FEM data. The MFL permeance of the solenoid valve G 3 can be improved. Now, assuming that P is the FEM data matrix and Q is the DPM data matrix, Er can be written as follows: If ∆G depicts the permeance error of the sampling point, Er can be derived as follows: where n is the PM segment number, m is the sampling point position. i = 1, 2, . . . n. j = 1, 2, . . . m. The basis function can be defined as f k (x 1 , x 2 . . . x q ), k = 1, 2, 3, . . . n, where q is the dimension. f k (x 1 , x 2 . . . x q ). The function can be written as: Based on Equation (20), the polynomial function is adopted for the armature stroke, l, ∆G and f (x 1 , x 2 . . . x q ). This function is written as follows: Equation (21) can also be written as follows: f er l, x 1 , x 2 , . . . , x q = X ∆G, l, f x 1 , x 2 , . . . , x q where f er (l, x 1 , x 2 . . . x q ) is the penalty function for the error data of the sampling point with armature displacement. The ith MFL permeance or soft magnetic resistance can be defined as q i (x), where x is the unknown quantity. Therefore, w i (x) can be obtained as follows: where w i (x) is the improved DPM data, and it is very close to the FEM data.

MFL Permeance Prejudge the Error Data Based on Kriging PIP
To improve the overall calculation efficiency of the program, the predictive identification program is applied to prejudge the scatter data before employing the kriging base function.
The basis function can be defined as f (x), a = (x 1 , x 2 , . . . x n ) ∈ R n , where f (x) is continuous and differentiable in the field of definition.
where f (x) is the maximum with where f (x) is the minimum with x = x k . The extreme function can be analyzed according to Equations (24) and (25), and the gradient of f (x) can be written as follows: Based on the characteristic of the first derivative, the adjacent points x i and x i+1 are the local maxima or minima. Therefore, the adjacent points are used as the criteria to judge the curve of the function. If the adjacent points of the function data are monotonically increasing, and the second derivative product is positive, this shows that the function is unimodal and it has a maximum. Therefore, it can be obtained as follows: Similarly, From Equation (28), the function is unimodal and has a minimum. Similarly, If the adjacent points of the function data are monotonically increasing or decreasing, then the function has a multi-crest. As such, the derivative can identify the extremum of the scatter data, and the shape of the multi-crest based on the number of extremum points.
To identify the more complicated scatter data curve, the standard function for numerical matching of the scatter data is applied. When the trend of the curve for the standard function is the same as the trend of the distribution for the scatter data, it can be said that the corresponding kriging basis function is practical for this type of standard function.
Based on the trend of the ∆G i curve, the Schwefel function has similar characteristics. Figure 6 shows the standard Schwefel function.
To identify the more complicated scatter data curve, the standard function for numerical matching of the scatter data is applied. When the trend of the curve for the standard function is the same as the trend of the distribution for the scatter data, it can be said that the corresponding kriging basis function is practical for this type of standard function.
Based on the trend of the ΔGi curve, the Schwefel function has similar characteristics. Figure 6 shows the standard Schwefel function. The Schwefel function can be written as follows: From Equation (30), the parameters can be modified to have the following: The Schwefel functions are compared with ΔG3 and ΔG4, as shown in Figure 7. The similarity between the Schwefel function, ΔG3 and ΔG4 can be obtained using the Tanimoto method. This method is stated as follows: The Schwefel function can be written as follows: From Equation (30), the parameters can be modified to have the following: The Schwefel functions are compared with ∆G 3 and ∆G 4 , as shown in Figure 7.
To identify the more complicated scatter data curve, the standard function for numerical matching of the scatter data is applied. When the trend of the curve for the standard function is the same as the trend of the distribution for the scatter data, it can be said that the corresponding kriging basis function is practical for this type of standard function.
Based on the trend of the ΔGi curve, the Schwefel function has similar characteristics. Figure 6 shows the standard Schwefel function. The Schwefel function can be written as follows: From Equation (30), the parameters can be modified to have the following: The Schwefel functions are compared with ΔG3 and ΔG4, as shown in Figure 7. The similarity between the Schwefel function, ΔG3 and ΔG4 can be obtained using the Tanimoto method. This method is stated as follows: The similarity between the Schwefel function, ∆G 3 and ∆G 4 can be obtained using the Tanimoto method. This method is stated as follows: By calculation, the similarity of ∆G 3 with the Schwefel function is 0.7654, and the similarity of ∆G 4 with the Schwefel function is 0.7775. The similarities of ∆G with the Schwefel function are further detailed in Table 2. In Table 2, the similarities of ∆G 3 and ∆G 4 with the Schwefel function are over 0.7, while the similarities of other ∆G 1 , ∆G 2 , ∆G 5 , ∆G 6 and ∆G 7 are the same as the trigonometric function. The similarities of ∆G are compared with the trigonometric function, as shown in Table 3. It can be concluded that the similarities which are over 0.7 conform to the corresponding standard function. Hence, the kriging basis function can be selected. When the similarity of the ∆G curve and the Schwefel function is over 0.7, ∆G can be modified by the Gaussian function. When the similarity of the ∆G curve and the trigonometric function is over 0.7, ∆G can be modified by the Fourier function.

The Error Correction of MFL Permeance and Soft Magnetic Resistance of Solenoid Valves
To verify the practicality of the proposed method, the kriging model based on the basis function adjustment is used to modify the error of the MFL permeance and soft magnetic resistance data of the solenoid valve magnetic system. Substituting modified data into the DPM, the electromagnetic force can be obtained.
From Figure 8, it can be observed that the unmodified DPM permeance mean error reaches 13.1% (see Figure 2), and the modified DPM mean error reaches 4.7%.
Actuators 2021, 10, x FOR PEER REVIEW 10 of 14 By calculation, the similarity of ΔG3 with the Schwefel function is 0.7654, and the similarity of ΔG4 with the Schwefel function is 0.7775. The similarities of ΔG with the Schwefel function are further detailed in Table 2.  Table 2, the similarities of ΔG3 and ΔG4 with the Schwefel function are over 0.7, while the similarities of other ΔG1, ΔG2, ΔG5, ΔG6 and ΔG7 are the same as the trigonometric function. The similarities of ΔG are compared with the trigonometric function, as shown in Table 3. It can be concluded that the similarities which are over 0.7 conform to the corresponding standard function. Hence, the kriging basis function can be selected. When the similarity of the ΔG curve and the Schwefel function is over 0.7, ΔG can be modified by the Gaussian function. When the similarity of the ΔG curve and the trigonometric function is over 0.7, ΔG can be modified by the Fourier function.

The Error Correction of MFL Permeance and Soft Magnetic Resistance of Solenoid Valves
To verify the practicality of the proposed method, the kriging model based on the basis function adjustment is used to modify the error of the MFL permeance and soft magnetic resistance data of the solenoid valve magnetic system. Substituting modified data into the DPM, the electromagnetic force can be obtained.
From Figure 8, it can be observed that the unmodified DPM permeance mean error reaches 13.1% (see Figure 2), and the modified DPM mean error reaches 4.7%.  The soft magnetic resistance in la 1 and la 2 segments is taken as an example, without the presence of an excitation current in the coil and at different armature positions. Figure 9 shows the soft magnetic resistance of FEM and the unmodified and modified DPM. The unmodified DPM permeance mean error reaches 9.94%, while the modified DPM permeance mean error reaches 3.7%.
Actuators 2021, 10, x FOR PEER REVIEW 11 of 14 The soft magnetic resistance in la1 and la2 segments is taken as an example, without the presence of an excitation current in the coil and at different armature positions. Figure  9 shows the soft magnetic resistance of FEM and the unmodified and modified DPM. The unmodified DPM permeance mean error reaches 9.94%, while the modified DPM permeance mean error reaches 3.7%.   The modified MFL permeance and soft magnetic resistance data are substituted into the DPM. The FEM results are compared with the modified and unmodified DPM results of the force on the armature for a coil excitation of 1800 AT and −1800 AT, as shown in Figure 10. The unmodified DPM mean error is 10.2%, and the modified DPM mean error is 3.8%. The FEM and measurement results mean error is 3%. The soft magnetic resistance in la1 and la2 segments is taken as an example, without the presence of an excitation current in the coil and at different armature positions. Figure  9 shows the soft magnetic resistance of FEM and the unmodified and modified DPM. The unmodified DPM permeance mean error reaches 9.94%, while the modified DPM permeance mean error reaches 3.7%. The modified MFL permeance and soft magnetic resistance data are substituted into the DPM. The FEM results are compared with the modified and unmodified DPM results of the force on the armature for a coil excitation of 1800 AT and −1800 AT, as shown in Figure 10. The unmodified DPM mean error is 10.2%, and the modified DPM mean error is 3.8%. The FEM and measurement results mean error is 3%.    Figure 11 shows the flow chart of the improved method for the DPM of solenoid valves based on the kriging basis function PIP. An appropriate function is selected by the PIP for error correction, and by substituting modified data into the DPM, the electromagnetic force can be obtained. Under the same calculation conditions, the output force of one armature displacement point cost 420 s in the FEM, and 56 s in the DPM. MFL Figure 11. Flow chart of method for DPM based on kriging basis function predictive identification program (PIP) error correction.

Conclusions
To improve the calculation accuracy and efficiency of the solenoid valve DPM, an improved method for the DPM based on the kriging basis function PIP is proposed in this paper. From the experiments and the results, the following conclusions are drawn: 1. Based on the characteristics of the kriging basis function curve, the relationship between the kriging basis function and the MFL permeance error data can be obtained, and an appropriate function is selected by contrasting various basis functions with error data curves. Then it is applied to gain error compensation between the FEM and DMP data. The PIP is introduced to prejudge the error data by comparing the standard function to the selected basis function. The modified MFL permeance and the soft magnetic resistance data are then substituted into the DPM of the electromagnetic device to calculate the attraction force. 2. The proposed method can effectively improve the calculation accuracy of the solenoid valve electromagnetic system. Compared with the FEM data, the unmodified MFL permeance of the DPM mean error is 13.1%, and the modified MFL permeance of the DPM mean error is 4.7%. The unmodified MFL permeance of the DPM mean error is 9.94%, and the modified MFL permeance of the DPM mean error is 3.7%. 3. The results of the DPM solenoid valve electromagnetic system in the case study showed a significant improvement. Particularly, the calculation accuracy improved by reducing the DPM mean error from 10.2% to 3.8%.

Conclusions
To improve the calculation accuracy and efficiency of the solenoid valve DPM, an improved method for the DPM based on the kriging basis function PIP is proposed in this paper. From the experiments and the results, the following conclusions are drawn:

1.
Based on the characteristics of the kriging basis function curve, the relationship between the kriging basis function and the MFL permeance error data can be obtained, and an appropriate function is selected by contrasting various basis functions with error data curves. Then it is applied to gain error compensation between the FEM and DMP data. The PIP is introduced to prejudge the error data by comparing the standard function to the selected basis function. The modified MFL permeance and the soft magnetic resistance data are then substituted into the DPM of the electromagnetic device to calculate the attraction force.

2.
The proposed method can effectively improve the calculation accuracy of the solenoid valve electromagnetic system. Compared with the FEM data, the unmodified MFL permeance of the DPM mean error is 13.1%, and the modified MFL permeance of the DPM mean error is 4.7%. The unmodified MFL permeance of the DPM mean error is 9.94%, and the modified MFL permeance of the DPM mean error is 3.7%.

3.
The results of the DPM solenoid valve electromagnetic system in the case study showed a significant improvement. Particularly, the calculation accuracy improved by reducing the DPM mean error from 10.2% to 3.8%.

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