Rate-Dependent Hysteresis Model of a Giant Magnetostrictive Actuator Based on an Improved Cuckoo Algorithm

: A rate-dependent asymmetric Prandtl–Ishilinskii (RAPI) model was proposed to tackle the serious rate-dependent hysteresis nonlinearity of the giant magnetostrictive actuator(GMA) output. First, a polynomial function was introduced based on the PI model, and hysteresis factors were introduced to the Play operator, which accurately described the asymmetrical characteristic of the actuator output. On this basis, rate-dependent parameters were added to establish a rate-dependent RAPI model. Second, an improved cuckoo search (ICS) algorithm was proposed to solve the difﬁculty in the parameter identiﬁcation of the RAPI model. For the ICS algorithm, the algorithm stability and optimization accuracy were improved using the adaptive step (AS) strategy and bird’s nest disturbance strategy. Then, the effectiveness of the ICS algorithm was tested by comparing it with other parameter identiﬁcation algorithms. Finally, the rate-dependent RAPI model was veriﬁed by combining the output data of the giant magnetostrictive actuator under different frequencies. The results show that the rate-dependent RAPI model exhibits a higher accuracy than the PI model, thus verifying the effectiveness of the rate-dependent RAPI model.


Introduction
A giant magnetostrictive actuator, a new-type actuator with giant magnetostrictive materials as the drive element, is a power output device with the magnetostrictive effect, accompanied by the advantages of a high output frequency, strong stability, and large output force [1].In recent years, the increasing speed of high-speed trains has increased the vibration of the train due to the irregularities of the track.The application of a giant magnetostrictive actuator to a train's active suspension allows for different damping forces to be output based on the disturbance excitation of distinct lines, thereby lessening train vibration and enhancing the stability and comfort of the train operation.Therefore, the application of a super magnetostrictive actuator in train active suspension has a good prospect.
Due to the serious hysteresis nonlinearity of the giant magnetostrictive actuator output, however, it is difficult to control the precise output, which greatly limits the practical application of giant magnetostrictive actuators [2].The hysteresis characteristic of a giant magnetostrictive actuator can be expressed by a hysteresis model, and establishing an accurate hysteresis model has become a significant foundation for solving hysteresis nonlinearity.At present, hysteresis models are generally divided into two types: physical models and phenomenological models, the former of which mainly include the J-A model [3] and the free energy model [4] and the latter of which include the Preisach model [5], the Prandtl-Ishilinskii (PI) model [6], and neural network models [7].In application, the J-A model and the free energy model require high-precision experimental data and more model parameters.The Preisach model is characterized by dual integrals and a complicated solving process, thus restricting its application.The solution precision of neural network models excessively relies upon the accuracy of experimental data.Over the other methods, the PI model integrates the advantages of a small calculated quantity, high flexibility, and easy identification and solution.In this study, therefore, a giant magnetostrictive actuator was subjected to hysteresis modeling based on the PI model.
The PI model was put forward by Prandtl et al. to model magnetic materials.In the process of popularization, the basic PI model has been improved by research scholars and applied in various fields.Zhong, Y. [8] introduced a cubic polynomial into the primary function of the PI model to construct an enhanced PI model, endowing the model with asymmetry so that it could describe the magnetic hysteresis phenomenon within a wider range and identify model parameters using the moth-flame algorithm.The results showed that the fitting accuracy of the improved PI model is significantly higher than that of the classical PI model under high frequencies.Zhao, X.H. [9] divided hysteresis nonlinearity into a linear module and a nonlinear model, described the nonlinear model with the PI model, replaced the linear model with an electromechanical transfer function, and finally connected the PI model with the transfer function in series to form a dynamic hysteresis model and solve the inverse model.The results manifested that the designed dynamic PI model displays good performance within 0-50 Hz.Tang, H.B. [10] improved the PI model by adding a dead-zone operator and introducing a frequency function.It was then experimentally verified that the improved model is much better than the classical PI model and can fit the GMM output curve well within 80 Hz.Pan, Y. [11] optimized the parameters of the classical PI model using the sequential quadratic programming algorithm, added a shift operator to make the model asymmetric, and improved the rate dependence of the model through the series ARX model given the shortcomings of the asymmetric PI model.Finally, it was verified that the model has a good fitting effect on hysteresis curves above 200 Hz, but the improving effect on low-frequency hysteresis curves is not evident.Subhash Rakheja et al. [12] established a rate-dependent PI model based on the memoryless function of the dead-zone operator and compared the comprehensive PI model data and the actual output.The comparison results revealed that the improved model can describe the nonlinear hysteresis characteristics in a wide range of excitation amplitudes and frequencies.Dargahi et al. [13] described the hysteresis curve of the magnetorheological elastomer using the improved generalized PI model and replaced the original equal interval distribution with the threshold distribution function.The results showed that the improved generalized model can accurately characterize the hysteresis characteristics of magnetorheological materials under different loading conditions and magnetic field ranges.Peng et al. [14] applied the PI model to the hysteresis nonlinear modeling of a giant magnetostrictive actuator, proposed a rate-dependent modified PI model, and then identified the parameters of the model by using the constrained least squares method.The results showed that, within a certain frequency range, the rate-dependent modified PI model could effectively describe the hysteresis characteristics of the giant magnetostrictive actuator.Xiao et al. [15] modified the classical PI model and used the least square method to identify parameters of the modified PI model.The results show that the modified PI model can reduce the modeling error, but the model accuracy is still not guaranteed under a high frequency.Nie, L.L. et al. [16] established a rate-dependent asymmetric hysteresis RDAPI model by introducing a dynamic envelope function into the play operator of the PI model and then combined with other improved PI models to characterize the rate-dependent and asymmetric hysteresis behavior of the piezoelectric micro-positioning level.Experiments show that the accuracy of the RDAPI model is improved significantly.On this basis, a robust adaptive control method is proposed.Wang, W. et al. [17] improved the Play operator in the PI model to an asymmetric operator, established an asymmetric API model, and optimized the parameters of the API model, which reduced the number of parameters and maintained a high accuracy.The experimental results show that the compensator based on this model can effectively suppress the hysteresis of the piezoelectric actuator.
In summary, the parameter identification and hysteresis modeling techniques for the magnetostrictive actuator based on the PI model are still immature.In this paper, aiming at the rate-dependent asymmetric hysteresis characteristics of the output curve of the active suspension magnetostrictive actuator, a rate-dependent asymmetric PI (RAPI) model was constructed by adding auxiliary functions based on the PI model, which can accurately describe the hysteresis curve of the output of the magnetostrictive actuator.In order to solve the problem of the RAPI model having multiple parameters that are difficult to identify, an improved cuckoo search (ICS) algorithm was proposed, which has higher optimization accuracy and stability compared to traditional algorithms.Finally, through comparative experiments with the PI model, it is verified that the RAPI model effectively improves the fitting accuracy of the actuator's hysteresis curve.Based on the RAPI model, precise dynamic control of the giant magnetostrictive actuator can be further achieved, which is of great significance for the practical application of the magnetostrictive actuator in the active suspension of trains.

PI Model
In the PI model, the hysteresis of the output is generally described by the weighted superposition of the Play operator [18].The basic Play operator is defined as follows: C m [0,t E ] is assumed to be the set of all continuous piecewise functions in the definitional domain of [0,t E ], in which 0 = t 0 < t 1 < t 2 . . .< t E , u(t)∈C m , and the output of the Play operator is: Here, F r [u(0)] is the output of a single Play operator at the initial moment; F r [u(t)] denotes the output of a single Play operator at time t, in which u(t) stands for the input signal at time t; F r [u(t − 1)] represents the output of the Play operator at the time before t.r is the threshold of the Play operator.
The output of the PI model is obtained through the weighted superposition of multiple Play operators with different thresholds and the summation with the input linear function, with the following output formula: Here, PI(t) is the output of the PI model; u(t) denotes the input of the PI model; α is the weight of the input; n represents the number of Play operators in the PI model; P ri stands for the weight of each Play operator.

Asymmetric PI (API) Model
The Play operator of the PI model is centrosymmetric, and the output of giant magnetostrictive actuators is characterized by asymmetric hysteresis, leading to a big error in the description of the asymmetric hysteresis curve by the PI model.Therefore, the PI model should be improved to make up for this deficiency.If introduced into the PI model, high-order polynomials, which are generally asymmetric, can endow the model with the ability to fit asymmetric curves.In this study, the linear function in the PI model was substituted by a polynomial function as follows: Despite the asymmetry of the model output thanks to the polynomial function, the model details were not described, so the hysteresis factor β was introduced into the model to optimize the operator Play.The improved expression is as follows: Here, F r , β [u(t)] is the output of the improved Play operator at the current moment; F r , β [u(t − 1)] stands for the output of the operator at the previous moment.The value of β could change the output amplitude and output gap of the Play operator, thus affecting the output of the asymmetric PI model.With the above improvements combined, the asymmetric PI model is expressed as shown in Formula (7).

RAPI Model
A study [19] has shown that when the excitation frequency is greater than 5 Hz, the hysteresis characteristics of giant magnetostrictive actuators will be affected by the frequency, showing rate-dependent hysteresis characteristics.Since the PI model is rateindependent, and neither the added hysteresis factor nor the auxiliary function contain frequency parameters, rate dependence processing is required for the asymmetric PI model.
The current rate-dependent hysteresis modeling methods can be divided into two types: separated rate-dependent hysteresis models and overall rate-dependent hysteresis models.The separated rate-dependent hysteresis models divide the hysteresis nonlinearity of the system into rate-independent hysteresis links of giant magnetostrictive materials and linear time invariant links consisting of other components, which are then respectively modeled and connected in series so as to establish an integral system model.Overall rate-dependent hysteresis models do not differentiate the system and establish a holistic model.Given that the model accuracy of separated modeling will be reduced under the high input frequency, the asymmetric PI model was improved in terms of the aspect of rate dependence through the overall modeling method.Moreover, the coefficient c of the weighted Play operator sum was added, and a frequency function specific to the auxiliary function G[u(t)] and the hysteresis factor β was constructed: Here, c(f ), β(f ), a 0 (f ). . .a n (f ) denote frequency-related parameters.The parameter values were identified by combining the output under different frequencies, and the unknown parameters in Formula (9) were acquired through curve fitting.Then, each frequency function was substituted into Formula (8) so as to obtain the rate-dependent RAPI model.

CS Algorithm
With the development of meta-heuristic optimization algorithms in recent years, such new-type optimization algorithms as the artificial fish swarm algorithm [20], CS algorithm [21], and grey wolf algorithm [22] have been developed by researchers and gradually applied to the field of parameter identification [23,24].Therein, the CS algorithm has been widely used by virtue of a few control parameters and a strong global optimization ability.The CS algorithm is a meta-heuristic optimization algorithm proposed by Yang et al., imitating the predation (Levy flight) and brooding behavior (discovery probability (Pa)) of cuckoos.The algorithm assumes the following conditions for simulation: • There was only one bird egg in each nest, and the eggs were randomly distributed.

•
In each evolution, the bird egg with the highest fitness was reserved.

•
Under a fixed number of bird's nests, the host bird abandoned the cuckoo eggs at Pa, and the eliminated nests would be replaced by new nests.
Based on the above assumptions, the iterative formula for the CS algorithm is as follows: Here, X is the position of the i-th bird's nest in the t + 1 iteration; X t i stands for the position of the i-th bird's nest in the t-th iteration; X t best represents the best bird's nest for the t-th iteration; α is the factor that controls the step of cuckoo movement, usually taken as 0.01; ⊗ is the product of vectors; Levy(λ) is a random number that obeys the Levy distribution, which is solved by the following formula: Here, µ~N(0,1), ν~N(0,1), β was taken as 1.5, Φ through the following formula: Here, τ is the standard Gamma function.Based on the above formula, the formula of the CS algorithm is as follows: When the bird's nests were updated according to Levy flight, the host bird would still abandon some nests as per Pa if it found the cuckoo's eggs, and a random number distributed in [0, 1] was generated through the rand function and compared with Pa.Stochastically biased walk was adopted for the nests whose random number was greater than Pa, while the nests satisfying < Pa were kept unchanged.The random preference formula is as follows: X t j and X t k are two random nests of the t-th generation, and r is a random number following a normal distribution.
The Cuckoo Search algorithm has a simple iterative process and a small number of control parameters.However, it also has limitations: the step size and discovery probability during the optimization process remain constant, which leads to the inability of the Cuckoo Search algorithm to balance global search and local exploration.For complex problems, it cannot rely on Levy flights to escape from local optima.

Adaptive Step (AS) Strategy
In the CS algorithm, α is the step coefficient, which directly determines the flight distance of the bird's nest in each iteration.Since α is taken as 0.01 in the standard CS algorithm, the flight distance of the bird's nest is relatively fixed, resulting in an imbalance between global search and local development.By analyzing the step formula, it can be concluded that under a rather lower α value, the iteration range of the algorithm is small, accompanied by a high search accuracy and a low convergence rate.When the α value is rather higher, the iterative step of the algorithm increases, along with an enlarged search scope and low convergence accuracy [25].Based on the above analysis, an AS function was proposed in this study and substituted into the following formula: According to the experiment, w min = 0.8 is the minimum value of the inertia factor; w max = 2 is the maximum value of the inertia factor; K is the number of iterations; and K max is the maximum number of iterations, which was set to K max = 1000.The change in the w value is displayed in Figure 1: The Cuckoo Search algorithm has a simple iterative process and a small number of control parameters.However, it also has limitations: the step size and discovery probability during the optimization process remain constant, which leads to the inability of the Cuckoo Search algorithm to balance global search and local exploration.For complex problems, it cannot rely on Levy flights to escape from local optima.

Adaptive Step (AS) Strategy
In the CS algorithm, α is the step coefficient, which directly determines the flight distance of the bird's nest in each iteration.Since α is taken as 0.01 in the standard CS algorithm, the flight distance of the bird's nest is relatively fixed, resulting in an imbalance between global search and local development.By analyzing the step formula, it can be concluded that under a rather lower α value, the iteration range of the algorithm is small, accompanied by a high search accuracy and a low convergence rate.When the α value is rather higher, the iterative step of the algorithm increases, along with an enlarged search scope and low convergence accuracy [25].Based on the above analysis, an AS function was proposed in this study and substituted into the following formula: According to the experiment, wmin = 0.8 is the minimum value of the inertia factor; wmax = 2 is the maximum value of the inertia factor; K is the number of iterations; and Kmax is the maximum number of iterations, which was set to Kmax = 1000.The change in the w value is displayed in Figure 1: After the AS function was introduced, the algorithm was able to search globally at a large step in the initial stage of iterations and search locally at a small step in the later stage, showing an initially declining, then rising, and finally declining trend on the whole.This not only conformed to the algorithm optimization law but also contributed to the strong population diversity of the algorithm in the middle and later stages.After the AS function was introduced, the algorithm was able to search globally at a large step in the initial stage of iterations and search locally at a small step in the later stage, showing an initially declining, then rising, and finally declining trend on the whole.This not only conformed to the algorithm optimization law but also contributed to the strong population diversity of the algorithm in the middle and later stages.

Bird's Nest Disturbance Strategy
Although the standard CS algorithm can jump out of the local optimal solution to a certain extent with the help of the Levy flight strategy in the later iteration, a long-step flight is carried out after repeated small-step iterations so that the algorithm cannot get rid of the local optimum in the later iterations.Hence, a bird's nest disturbance strategy was proposed to disturb the population with abnormal monitoring results by monitoring the results of each iteration.
The specific implementation steps are as follows: The bird's nests that kept the global optimal solution unchanged in repeated iterations were recorded.When the global optimal solution was unchanged in the n-th iteration of the bird's nests, the current population was trapped in local optima, so it was disturbed at a disturbance step, which was the linear difference between the optimal solution of the current generation and the position of the bird's nest.In this way, the algorithm was capable of optimization near the local optimal solution, ensuring that the algorithm could still remain convergent after the disturbance.The disturbance step length formula is as follows: ) Here, V i,j is the j-th dimensional disturbance step in the i-th cuckoo nest; ζ 1 is the step factor controlling the disturbance radius, in which I is the current population number; N is the maximum population number; s is a constant; ζ 2 stands for the weight coefficient of the contemporary optimal cuckoo nest, taken as 0.03; X i,j is the j-th dimensional nest value in the i-th nest.The disturbance step V i,j was dynamically adjusted by the cuckoo population number V i,j , and the disturbance step factor of the head population was large, which enabled the algorithm to perform optimization in a larger scope.The disturbance step factor of the rear population was small, so the algorithm searched near the local optimal solution.The bird's nest disturbance strategy imposed different disturbance steps on different populations, which facilitated global search and local development with the local optimal solution as the center and effectively improved the optimization accuracy and search ability of the algorithm.

ICS Algorithm Flow
The algorithm flow is as follows.
Step1: Input the algorithm iteration number K max , the number of nests N, w max , and w min .
Step2: Initialize the nest positions and calculate the fitness value.
Step3: Update the nest positions using dynamic inertia Levy flight and retain nests with better fitness.
Step4: Select nests based on Pa, leaving nests with a better fitness value.Calculate the global optimal solution and the fitness value of the optimal solution.
Step5: Determine whether to disturb the population based on the number of repetitions of the global optimum in each iteration.
Step6: Determine whether the algorithm reaches the termination condition; if not, go back to Step3; when the termination condition is reached, output the optimal solution and the fitness value.
The flowchart of the algorithm is shown in Figure 2:

Performance Test of the ICS Algorithm
A set of mixed test functions was selected to verify the performance of the ICS algorithm.In the set, functions F1-F3, which are unimodal functions, could test the convergence ability and optimization speed of the algorithm.Functions F4-F6, which are multimodal functions with increasing test difficulties, a large scope of optimization, and multiple local optimal values in the definitional domain, were applied to test the optimization ability of the algorithm for multimodal complex functions.Function F7, a two-dimensional function, was used to test the optimization ability of the algorithm for low-dimensional problems.The expressions of such test functions are listed in Table 1.

Performance Test of the ICS Algorithm
A set of mixed test functions was selected to verify the performance of the ICS algorithm.In the set, functions F1-F3, which are unimodal functions, could test the convergence ability and optimization speed of the algorithm.Functions F4-F6, which are multimodal functions with increasing test difficulties, a large scope of optimization, and multiple local optimal values in the definitional domain, were applied to test the optimization ability of the algorithm for multimodal complex functions.Function F7, a two-dimensional function, was used to test the optimization ability of the algorithm for low-dimensional problems.The expressions of such test functions are listed in Table 1.
The ICS algorithm was comparatively analyzed with particle swarm optimization (PSO), CS, and ASCS to explore its accuracy and convergence.Table 2 gives the specific parameter settings of each algorithm.The number of test populations was set to 50, the number of dimensions was set to 30, and the maximum number of iterations was set to 1000.Each test function was run 30 times in the same environment, and their average function error (Mean), standard deviation (Std), and optimal value (Best) were recorded.The test results are shown in Table 3.It could be seen that in functions F1-F3, the ICS algorithm could accurately converge to the global optimum.F4, a complex multimodal function, exhibited a lot of local optimal values in the definitional domain.As shown in Tables 4 and 5, the ASCS algorithm was subjected to a too-large standard deviation and unstable optimization results despite the ability to search the optimal value of F4, while other algorithms all failed to search the global optimal value.The ICS algorithm could converge to the global optimum even after repeated tests.F5 and F6, which were classical complex functions, were featured with a relatively high optimization difficulty.For F6, the ICS algorithm was stuck in the local optimal value of 8.88 × 10 −16 in the later stage of iterations, but its optimization accuracy remained higher than that of other algorithms.For F7 (a simple low-dimension function), the ICS algorithm and other algorithms converged to the global optimum, indicating its excellent optimization ability on low-dimension functions.According to the Std value, the standard deviation of the ICS algorithm was 0 in many tests under different functions, manifesting its favorable stability in optimization.

Parameter Optimization
The RAPI model was established on the basis of the API model, the parameters of which, therefore, were identified first.Some parameters preset in the API model were optimized.

Optimization of the Number of Hysteresis Factors
Hysteresis factors could change the shape of the Play operator by influencing the threshold so as to generate an impact on the output, and the model will be affected, to different degrees, by the number of hysteresis factors.Different quantities of hysteresis factors were thereby experimented, and the maximum absolute deviation (MAD), the mean absolute error (MAE), and the root mean square error (RMSE) were used as evaluation indexes for the model fitting accuracy, as shown in Table 4.
Combining the displacement data of the sinusoidal current under the amplitude of 2 A and 10 Hz, the parameters of four schemes were identified, respectively.The number of operators was set to 10, the order of the auxiliary function was set to 3, and the rest of the parameters were obtained through the ICS algorithm.Table 5 presents the evaluation results under four different numbers of hysteresis factors.As the number of hysteresis factors increased, the model fitting error gradually increased, so the number of hysteresis factors was set to 1 in this study.

Order Optimization of the Auxiliary Function
The auxiliary function in the API model is the key factor deciding the model asymmetry, and its order affects the complexity of the model and the difficulty in parameter identification.Therefore, it was necessary to optimize the model order.An experiment was carried out through the control variable method to acquire the optimal number.Therein, the number of operators was set to 10, that of hysteresis factors was set to 1, the polynomial order was set to 1-6, and the other parameters were determined through the relevant parameter identification algorithm.
It could be observed from Figure 3 that the RMSE decreased considerably when the order of the auxiliary function increased from one to three and basically kept unchanged as the order grew from three to four, but it rose again when the order increased to five and six.This is because the increase in the function order aggravates the difficulty in parameter identification and leads to the failure to converge to the global optimum after reaching the given number of iterations.Moreover, the parameter identification time is substantially lengthened as the model order increases.In full consideration of the model accuracy and parameter identification time, the order of the auxiliary function was chosen to be three, expressed in Formula (20).
Actuators 2023, 12, x FOR PEER REVIEW 11 of 16 factors increased, the model fitting error gradually increased, so the number of hysteresis factors was set to 1 in this study.

Order Optimization of the Auxiliary Function
The auxiliary function in the API model is the key factor deciding the model asymmetry, and its order affects the complexity of the model and the difficulty in parameter identification.Therefore, it was necessary to optimize the model order.An experiment was carried out through the control variable method to acquire the optimal number.Therein, the number of operators was set to 10, that of hysteresis factors was set to 1, the polynomial order was set to 1-6, and the other parameters were determined through the relevant parameter identification algorithm.
It could be observed from Figure 3 that the RMSE decreased considerably when the order of the auxiliary function increased from one to three and basically kept unchanged as the order grew from three to four, but it rose again when the order increased to five and six.This is because the increase in the function order aggravates the difficulty in parameter identification and leads to the failure to converge to the global optimum after reaching the given number of iterations.Moreover, the parameter identification time is substantially lengthened as the model order increases.In full consideration of the model accuracy and parameter identification time, the order of the auxiliary function was chosen to be three, expressed in Formula (20).

Optimization of the Number of Play Operators
As the basic operation unit of a hysteresis model, the Play operator plays an important role in the model output.The model fitting accuracy will decline under a small number of operations.A too-large number of operations will intensify the complexity of the calculation though improve the prediction accuracy to some extent.To determine the concrete number of operators, the experiment was carried out under different numbers of operators, with other parameters remaining the same.Figure 4 shows the RMSE curve under different numbers of operators.With the increase in the number of operators, the RMSE gradually decreased, indicating that the fitting accuracy of the model is gradually improved.The subsequent increase in the number of operators >10 generated a very minor influence on the RMSE.Hence, the number of operators was determined to be 10 to

Optimization of the Number of Play Operators
As the basic operation unit of a hysteresis model, the Play operator plays an important role in the model output.The model fitting accuracy will decline under a small number of operations.A too-large number of operations will intensify the complexity of the calculation though improve the prediction accuracy to some extent.To determine the concrete number of operators, the experiment was carried out under different numbers of operators, with other parameters remaining the same.Figure 4 shows the RMSE curve under different numbers of operators.With the increase in the number of operators, the RMSE gradually decreased, indicating that the fitting accuracy of the model is gradually improved.The subsequent increase in the number of operators >10 generated a very minor influence on the RMSE.Hence, the number of operators was determined to be 10 to balance the prediction accuracy and the complexity of the calculation.In this study, the model parameters were identified by combining the experimental data under the current amplitude increasing in the range of 1-5 A and the frequency of 1 Hz.In this process, the number of populations of the ICS algorithm was set to 80, the maximum number of iterations was set to 1000, and the disturbance threshold was set to 3. In the PI model, the weight value was generally positive, but it was experimentally found that a lower RMSE could be acquired when the weight value was negative, so a negative weight value could be assigned to the API model.Since the API model is the sum of the weighted superposition of Play operators, the model characteristics would not be changed by a small number of negative values.The ICS algorithm was repeatedly run ten times, the optimal value and the optimal bird's nest were recorded, and the API model parameters finally obtained are listed in Tables 6 and 7.In this study, the model parameters were identified by combining the experimental data under the current amplitude increasing in the range of 1-5 A and the frequency of 1 Hz.In this process, the number of populations of the ICS algorithm was set to 80, the maximum number of iterations was set to 1000, and the disturbance threshold was set to 3. In the PI model, the weight value was generally positive, but it was experimentally found that a lower RMSE could be acquired when the weight value was negative, so a negative weight value could be assigned to the API model.Since the API model is the sum of the weighted superposition of Play operators, the model characteristics would not be changed by a small number of negative values.The ICS algorithm was repeatedly run ten times, the optimal value and the optimal bird's nest were recorded, and the API model parameters finally obtained are listed in Tables 6 and 7. On the basis of the API model, the rate-related parameters were identified.To reduce the influence of environmental noise on the test displacement as much as possible, the other model parameters were experimentally identified by selecting the sinusoidal current with an amplitude of 5A, the signals with the frequency of [1,25,40,55,70,85,100,120], and the Play operators with a constant weight value.Table 8 shows the values of the model coefficients at different frequencies obtained by the parameter identification algorithm.

Verification of the RAPI Model
In order to obtain the hysteresis input-output data required for the identification model, a giant magnetostrictive actuator test bench is built, as shown in Figure 5.The experimental bench mainly includes: the giant magnetostrictive actuator, YB1635 signal generator, NI-usb6211 acquisition card and display device, Aigtek-ATA309 power amplifier, Panasonic HG-C1030 laser displacement sensor, and PC and Tektronix-MDO3014 oscilloscope.Among them, the signal generator generates a positive rotating AC signal, which is amplified by the power amplifier to drive the vibration of the super magnetostrictive actuator with a suitable excitation current.The laser displacement sensor, oscilloscope, acquisition card, and display device cooperate with the PC to complete the display, acquisition, processing, and post-analysis of data.On the basis of the API model, the rate-related parameters were identified.To reduce the influence of environmental noise on the test displacement as much as possible, the other model parameters were experimentally identified by selecting the sinusoidal current with an amplitude of 5A, the signals with the frequency of [1,25,40,55,70,85,100,120], and the Play operators with a constant weight value.Table 8 shows the values of the model coefficients at different frequencies obtained by the parameter identification algorithm.

Verification of the RAPI Model
In order to obtain the hysteresis input-output data required for the identification model, a giant magnetostrictive actuator test bench is built, as shown in Figure 5.The experimental bench mainly includes: the giant magnetostrictive actuator, YB1635 signal generator, NI-usb6211 acquisition card and display device, Aigtek-ATA309 power amplifier, Panasonic HG-C1030 laser displacement sensor, and PC and Tektronix-MDO3014 oscilloscope.Among them, the signal generator generates a positive rotating AC signal, which is amplified by the power amplifier to drive the vibration of the super magnetostrictive actuator with a suitable excitation current.The laser displacement sensor, oscilloscope, acquisition card, and display device cooperate with the PC to complete the display, acquisition, processing, and post-analysis of data.The experimental curves at 50, 80, and 110 Hz were randomly selected and compared with the curves acquired by the PI model and the RAPI model.The fitting error curves of the two models are shown in Figure 6.It could be observed from the fitting curves of the two models that the RAPI model curves introduced with a frequency function could accurately fit the hysteresis curve of the giant magnetostrictive actuator at variable frequen- The experimental curves at 50, 80, and 110 Hz were randomly selected and compared with the curves acquired by the PI model and the RAPI model.The fitting error curves of the two models are shown in Figure 6.It could be observed from the fitting curves of the two models that the RAPI model curves introduced with a frequency function could accurately fit the hysteresis curve of the giant magnetostrictive actuator at variable frequencies.The error curve diagram shows that the PI model fails to describe the rate-related characteristics.As a result, as the frequency increases, the output error values gradually increase.After the RAPI model is improved with rate correlation, the output error values at different frequencies remain relatively stable, indicating a significant improvement in the fitting accuracy of the RAPI model for the variable frequency output curve of the magnetostrictive actuator.Table 9 exhibits the evaluation results of the RAPI model and the PI model from three evaluation indexes: RMSE, MAD, and MAE.The results revealed that the RMSE of the model was reduced by 3.61 µm, the MAE was reduced by 3.5 µm, and the MAD was reduced by 7.5 µm after the rate-dependence improvement, further verifying the effectiveness of rate-dependence improvement.

Conclusions
A modified RAPI model was proposed based on the hysteresis characteristics of the magnetostrictive actuator.The RAPI model first includes a polynomial function to describe the asymmetry of the hysteresis curve.On this basis, a global modeling approach was used to introduce frequency parameters, enabling the model to describe the raterelated characteristics of the output from the giant magnetostrictive actuator.To accurately identify the parameters of the RAPI model, an adaptive step ICS algorithm with bird-nest disturbance was designed.The optimized ICS algorithm was analyzed and compared with PSO, CS, and ASCS algorithms, showing that the ICS algorithm has higher optimization accuracy and stability.Finally, experimental validation was conducted, and the results show that the RMSE, MAE, and MAD values of the RAPI model were respectively reduced by 3.61 µm, 3.5 µm, and 7.5 µm compared to the PI model, confirming that the accuracy of the RAPI model is far higher than that of the traditional PI model.
The RAPI model accurately describes the hysteresis characteristics of the giant magnetostrictive actuator.On this basis, the controller and control algorithm can be further designed to accurately control the giant magnetostrictive actuator and realize its application in the train active suspension.

Figure 1 .
Figure 1.Variation curve of the inertia factor.

Figure 1 .
Figure 1.Variation curve of the inertia factor.

Figure 2 .
Figure 2. The flowchart of the ICS algorithm.

Figure 2 .
Figure 2. The flowchart of the ICS algorithm.

Figure 3 .
Figure 3. Diagram of RMSE and optimization time under different orders of the auxiliary function.

Figure 3 .
Figure 3. Diagram of RMSE and optimization time under different orders of the auxiliary function.

Figure 4 .
Figure 4. RMSE of the CAPI model under different numbers of operators.

Figure 4 .
Figure 4. RMSE of the CAPI model under different numbers of operators.

Actuators 2023 ,
12,  x FOR PEER REVIEW 14 of 16 increase.After the RAPI model is improved with rate correlation, the output error values at different frequencies remain relatively stable, indicating a significant improvement in the fitting accuracy of the RAPI model for the variable frequency output curve of the magnetostrictive actuator.Table9exhibits the evaluation results of the RAPI model and the PI model from three evaluation indexes: RMSE, MAD, and MAE.The results revealed that the RMSE of the model was reduced by 3.61 μm, the MAE was reduced by 3.5 μm, and the MAD was reduced by 7.5 μm after the rate-dependence improvement, further verifying the effectiveness of rate-dependence improvement.(a) Fitting curves under 50 Hz (b) Error curves under 50 Hz (c) Fitting curves under 80 Hz (d) Error curves under 80 Hz (e) Fitting curves under 110 Hz (f) Error curves under 110 Hz

Figure 6 .
Figure 6.Fitting curves of the RAPI model under different frequencies.

Table 3 .
Experimental result of test functions.

Table 4 .
The expression of the evaluation index.

Table 5 .
Comparison of errors under different numbers of hysteresis factors.

Table 6 .
Parameter identification results of the API model.

Table 7 .
Hysteresis factors and auxiliary coefficients of the API model.

Table 6 .
Parameter identification results of the API model.

Table 7 .
Hysteresis factors and auxiliary coefficients of the API model.

Table 8 .
Parameter values of the RAPI model under different frequencies.

Table 8 .
Parameter values of the RAPI model under different frequencies.

Table 9 .
Evaluation results of the RAPI model under three frequencies.