Abstract
A linear magnetization model is built to replace the Jiles–Atherton model in order to describe the relationship between the magnetic field intensity and the magnetization intensity of the giant magnetostrictive vibrator (GMV). The systematic modeling of the GMV is composed of three aspects, i.e., the structural mechanic model, the magnetostrictive model, and the Jiles–Atherton model. The Jiles–Atherton model has five parameters to be defined; hence, its solution is so complex that it is not convenient in application. Therefore, the immune genetic algorithm (IGA) is applied in the identification of the five parameters of the Jiles–Atherton model and it showed a higher stability compared with the identification result of the differential evolution algorithm (DEA). The identification parameters of the two algorithms were employed, respectively, to calculate the excitation force and it was found that the relative error of IGA was evidently smaller than that of DEA, indicating that the former was more reliable than the latter. According to the identification results of IGA and based on the least square method (LSM), curve-fittings to the magnetic field intensity and magnetization intensity were conducted by using the linear function. And the linear magnetization model was built to replace the Jiles–Atherton model. Research results show that the linear model of the GMV can be established by combining the linear magnetization model with the structural mechanic model as well as the giant magnetostrictive model. The linear magnetization model, which has great engineering application value, can be applied in the open-loop control of the vibrator.
1. Introduction
Giant magnetostrictive material (GMM) is a functional material with superior performance and outstanding advantages, such as large range of magnetostrictive strain, high electrical motor conversion efficiency, quick response, high intensity of power, weak driving magnetic field, good frequency performance, and high Curie point [1]. The GMV, developed with the application of GMM, can make up for the deficiencies of the mechanical vibrator. For instance, it is difficult to meet the inherent frequency of the components with high rigidity as a consequence of its vibration frequency below 200 Hz [2]. The excitation force cannot be regulated smoothly once the vibration frequency is fixed. The reliability is low and the service life is short with the application of DC motor as the driving force. However, for a rather long period of time, the mechanical vibrator dominates in the field of vibration stress relief (VSR). Therefore, the research and development of the GMV is of great value to the further technological promotion of VSR and the acceleration of the updating of both the VSR equipment and the industry as well.
In recent years, in the application of GMM, the majority of researchers focused their study on giant magnetostrictive microdisplacement actuator [3], and little attention was given to the study of the GMV. Meng et al. [4] and Wang et al. [5] optimized the vibrator structure and the magnetic circuit design. They analyzed the magnetic field homogeneity of different magnetic circuits but did not mention the giant magnetostrictive model yet. The modeling of GMV involves three aspects: structural dynamics model, magnetization model, and magnetostrictive model which describes the relation between magnetization intensity and magnetostrictive coefficient and can be built with quadratic domain transformation model [6]. The magnetization model, which describes the relation between magnetic field intensity and magnetization intensity, can be built with Jiles–Atherton model [7]. Although the Jiles–Atherton model has the advantage of clear physical meanings, it has five unidentified parameters. As a consequence of its complex solution process, it is inconvenient for engineering application. In order to increase the local research ability, Cao et al. [8] adopted the strategy of combining the genetic algorithm and trust region algorithm and took the simulation method to identify the parameters of the Jiles–Atherton model. The identification results are found to be close to the partial known parameters after comparison. Knypinski et al. [9] and Rong et al. [10] studied the identification of parameters of the Jiles–Atherton model using the particle algorithm. Unfortunately, these abovementioned algorithms adopted the simulating data to replace the testing data, compared the identification result with the known parameters; however, for the given GMM, the undetermined parameters are related to the magnetic field; hence, the above identification results cannot be used directly.
In this paper, IGA was applied to identify the five parameters of the Jiles–Atherton model. The identification results of IGA and DEA were compared to verify the reliability of the parameters. The linear magnetization model was built to replace the Jiles–Atherton model based on the identification results of IGA and LSM and combining the fitting curve of the magnetic field intensity and magnetization intensity with linear function. The linear model of the GMV was built with the combination of linear magnetization model, structural dynamic model, and magnetostrictive model.
2. Overall Structure of the GMV
The GMV mainly consists of the drive coil, coil skeleton, GMM rod, magnetic conductive cover, magnetic conductive outer wall, output rod, pretightening spring, adjustment bolt, radiator, support plate, upper cover, lower end cover, tension bolt, housing and base, etc., and its structure is shown in Figure 1. Among them, GMM rod, magnetic conductive cover, and magnetic conductive outer wall form a close magnetic circuit. The pretightening spring exerts prepressure stress on the GMM rod, and the stress can be adjusted by the bolt so that it can give the most appropriate stress on the GMM rod to obtain the largest magnetostrictive coefficient [11]. To utilize the second harmonic generation (SHG) of the GMM [12], the bias magnetic field is not set and the coil and magnetic circuit structure is simplified. The drive coil adopts a scheme of reducing the number of turns, increasing the wire diameter, and driving with high electric current to improve its dynamic performance.

Figure 2 shows the magnetic field intensity-magnetostrictive coefficient curve under different prepressure stresses [13]. The optimal pretightening force is the basis of GMM rod design.

Figure 3 shows the relationship among the pretightening stress of the GMM bar, the magnetic field intensity, and the electromechanical coupling coefficient [14]. As can be seen from Figure 3, when the pretightening stress is , the maximum magnetostrictive coefficient is 1000 ppm and the electromechanical coupling coefficient is 0.6.

The main structural parameters of the GMV are shown in Table 1.
3. System Model of the GMV
3.1. Magnetostrictive Model and Magnetization Model
In real GMM, the electromechanical coupling is very complicated under the effect of alternating magnetic field. However, when applying a certain pretightening stress on the direction of giant magnetostrictive rod and considering that the material crystals prefer growing in the axial direction, the GMM magnetic domain will distribute along the tangential and preferred direction of magnetization; at this time, the relation between magnetostrictive coefficient of the GMM rod and magnetization intensity is approximately the quadratic domain transformation model based on the energy [15]:where is the saturated magnetic field intensity, is the saturated magnetostrictive coefficients, and is the magnetization intensity.
The micromagnetic theory holds that the domain of ferromagnetic material varies in structure with form and size. In addition, the stress does not distribute evenly, which makes the movement of the domain wall and the process of magnetization irreversible and leads to the containing of the pressure-stabilized state in the material energy free moving equation, so that it will cause energy loss in the magnetization process. Based on Weise molecular field theory and micromagnetics, Jiles and Atherton built theoretically the Jiles–Atherton model of ferromagnetic material [16]:where is related to the magnetic field intensity , when increases, , when decreases, ; is the prestress of the GMM rod. In the Jiles–Atherton model, the undetermined parameters are as follows: domain and wall interaction coefficient, irreversible loss coefficient, reversible coefficient, anhysteretic magnetization form factor coefficient, and saturated magnetization intensity.
3.2. Structural Dynamic Model of the GMV
Without taking into consideration the magnetic stagnant effect, the linear constitutive piezomagnetic equation of GMM is [17]where is the magnetostrictive strain, is the smooth coefficient, is the stress, is the piezomagnetic coefficient, and is the magnetic field intensity.
When GMM works under dynamic magnetic field, hysteresis nonlinearity will exist. The second item of the linear piezomagnetic equation is presented by quadratic domain transformation model, and according to all these, the nonlinear piezomagnetic equation can be deduced [18]:
When considering the mass and damper, formula (4) can be written as nonlinear piezomagnetic equation [19], including inertia item and damper item:where is the expansion and contraction quantity of the GMM bar.
The upper end of the GMM rod bears the inertia of rear mass, and the bottom bears the output rod reaction force. Since the forces are dynamic, the mass of the output rod, the base, the damper, and the rigidity should all be taken into consideration. The stress of the GMM rod consists of both the rear mass inertia and the output force of itself, namely,where is the output force of the GMM rod, is the inertial force of the rear mass, is the magnetostrictive stress, is the inertial stress of the rear mass, and is the cross sectional area of the GMM rod.
According to formula (4), the magnetostrictive stress can be written as
Based on formulas (3) and (4), the output force of the GMM rod can be written as
According to mechanic theory, the equilibrium between the GMM rod and the output rod is
By adjusting formula (9) and using Laplace transformation, we obtain the vibrator displacement equation as follows:where is the Laplace operator, is the total mass, is the total damper, and is the total stiffness.
The excitation force and the reaction force of the output rod are equal; however, they are in opposite direction, so
By using Laplace transformation for formula (8), we obtain
Substituting formulas (7) and (10) into (9), we obtainwhere is the equivalent mass and is the equivalent stiffness.
The inherent frequency of the vibrator and the damping ratio, respectively, are
From formulas (14) and (15), we obtain the inherent frequency of the vibrator , and the damping ratio .
With quadratic domain transformation model (1), Jiles–Atherton model (2), displacement model (10), and excitation force model (13), we can build a vibrator system model, and its structure is shown in Figure 4.

As shown in Figure 4, the system model consists of the hysteresis nonlinearity model and the second-order linear system model in series, the input is the driving current, and the output is the displacement and the force. There are five unidentified parameters to be identified and obtained in the Jiles–Atherton model.
3.3. Discrete Transformation of the Continuous Model
The displacement and the excitation force of the vibrator can be obtained by applying the quadratic domain transformation model, the Jiles–Atherton model, and the displacement and excitation force model. By using inverse Laplace transformation of equation (10), we obtain the differential equation of the oscillator displacement model as follows:
In formula (16), defining as the output vector, we obtain the continuous-time state-space equation of the displacement:
Transforming the excitation force formula (13), we write
Take as the state variable and as the output vector, from formulas (19) and (20), and we obtain the continuous-time state-space expression of the excitation force as follows:
Formulas (17) and (21) are the expressions of the vibrator displacement and excitation force, and by discrete transforming of these two formulas, respectively, we obtain the corresponding discrete-time state-space expression for displacement and excitation force as follows:where represents the sampling period.
4. Identification Principle and IGA
4.1. Principle of Parameter Identification
In the Jiles–Atherton model, the 5 parameters to be identified, i.e., , , , and , can be presented as
The identification of the parameters of the Jiles–Atherton model can be regarded as a minimum optimization and its objective function is as follows:where is the sampling time, is the total number of sample, is the measured value, is the calculated value, is the error of the excitation force, is vector th of parameter , and , is the top and bottom limit of , respectively.
4.2. Immune Genetic Algorithm
The evaluation function is judged by the adaptive value using the floating-point encoding scheme. The larger the adaptive value is, the better the individual will be. The adaptive value is the key information of the genetic operation. According to the problem to be solved, the estimation function [20] can be defined by the compound function in formula (25):where is the population algebra and is the individual sequence number.
Cross probability and mutation probability , which are key factors affecting the genetic behavior and algorithm performance, play an important role in the convergence of algorithm. The larger the is, the faster the new individual is reproduced and the higher the probability of damaging the genetic pattern. On contrary, the smaller is, the slower the search is, even close to a standstill [21]. If is too small, it is not apt to reproduce new individual; however, if is too large, it will become a complete random searching method. To overcome the abovementioned deficiencies, self-adaptive strategy is adopted, the population similarity is represented by information entropy, and cross probability and mutation probability are regulated dynamically. The information entropy [22] is calculated at first:where is the number of the total gene, is the sequence number, is the population quantity, and is the information entropy of the th gene:where is the quantity of gene and is the probability that the value of gen appears.
The similarity of the population is as follows:
represents the overall degree of the population similarity, satisfying , the larger the value is, the higher the similarity of the population will be, and when it equals to 1, it means that the antibodies of the two populations are totally the same [23].
The adaptive strategy dynamically adjusts the cross probability and the mutation probability according to the population diversity. When the individual adaptive values tend to converge, and will increase appropriately; and when the adaptive values disperse, and will decrease appropriately, as shown in the following formula:
4.3. First Identification Results
In the process of identification, the excitation current increased in the range of 1~30 A at a step size of 0.3 A, and the driving frequency increased in the range of 1~400 Hz at a step size of 2 Hz. The output pulse of each frequency point was repeated for 50 cycles.
The largest excitation force of each cycle is collected, the average excitation force of the 50 cycles is calculated, and 100 sets of data are collected as sample each time. The initial setting of the population size is , and the maximum number of iteration is . Taking 100 sets of data as samples, the noise interference is temporarily ignored, that is, the noise variance . DEA and IGA were adopted, respectively, to identify the parameters. The adaptive values are calculated uniformly with formula (26). The cross probability of DEA is , and the transformation operator of DEA is . The cross probability of IGA is , and the transformation operator of IGA is .To evaluate the identification effect of the algorithm, in addition to the objective function, the error of the excitation force was also taken into consideration. The measured excitation force was taken as the true value, and the calculated excitation force was taken as the theoretical value; hence, the relative error of the excitation force is expressed as
In which, is the relative error of the excitation force, is the measured excitation force, and is the calculated excitation force.
The first identification data of DEA and IGA are shown in Table 2, where is the convergence generation. By comparison and analysis of the data in Table 2, the objective function values of DEA and IGA are 0.873 and 0.635, respectively. Their difference is 0.238, and the relative value is 37% lower. As for the relative error of excitation force, it is 4.938% for IGA, which is lower than 5%; however, it is as high as 13.625% for DEA. When comparing the number of iteration, only the convergence generation of DEA is smaller than that of IGA. It needs further observation to comprehensively assess the advantages and disadvantages of the algorithms and the parameter identification results.
The objective function and identification parameters of the iteration process are shown in Figure 5. Figure 5(a) shows the objective iterating process. In the first 150 generations, DEA had the largest descending gradient; then, it started to slow down. From about the 500th generation, it changed little. For IGA, however, its gradient generally kept a descending trend from the beginning to the end of iteration and gradually approached the mimimal value, and its result was much smaller than that of DEA. To further compare the iterative process of the five parameters, the changing range was very large for the first 200 generations, indicating that the algorithm had a comparatively strong global search ability, the evolution of the population was very active in the early period, and after about 240 generations, the change of each parameter tended to stagnate, after the 300th generation, except for parameter a, the other parameters did not change anymore. Combining with the iterative process of the objective function, it can be interpreted that the premature phenomenon appeared in the population, and the result was trapped into regional optimal solution. However, for IGA, except for parameter , it started to converge stably from the 150th generation, no matter ascending or descending, the gradient was relatively moderate, and it approached gradually to the stable value, and till about the 350th generation; the current value basically reached about 90% of the final value. However, the evolution did not stagnate, the iteration changed from global search to region search, and till the 550th generation, the iteration process got the optimal solution.

(a)

(b)

(c)

(d)

(e)

(f)
4.4. Repeated Identification Results
In order to identify the duplicability and stability of the algorithms, the identification was repeated four times on the basis of the first identification, and the data can be found in Tables 3 and 4. Comparing the data in Tables 3 and 4, the mean value of DEA is as high as 0.917, and the mean value is 252; however, for IGA, the mean value is only 0.615 and the mean value is 332. It means that DEA can be easily trapped in the regional optimal solution and shows an obvious population premature phenomenon. The variances of the five parameters of IGA are all smaller than those of DEA, indicating that the identification data of the former have better duplicability and stability.
The average value of the five parameters was substituted into the model to calculate the excitation force and then was compared with the measured excitation force. The relative error of the excitation force was calculated, respectively, as shown Figure 6. Comparing the measured excitation force and calculated excitation force, the absolute error of IGA is much smaller than that of DEA. The error of the former is relatively large only when the current is less than 10 A, while the latter has a large error in all the range of excitation currents. The maximum relative error of DEA is 13%; however, the maximum relative error of IGA is lower than 5%. Hence, the result of IGA exhibits a better reliability.

(a)

(b)
The above identification results were obtained without considering noise interference. However, in the real test system, noise interference does exist definitely. Therefore, noise interference was added, namely, noise variance , then random identification was made for 20 times, and the data are shown in Table 5. The results show that, even in the case of noise interference, the identification results are consistent with those in Table 4, so the average value of the identification results of IGA is taken as the model parameter.
5. Local Linearization and the Application of the Magnetization Model
5.1. Characteristics of Hysteresis Loops
Based on the parameter identification results of the Jiles–Atherton model, the magnetization strength and magnetic field intensity at different frequencies are calculated, and the magnetic field intensity-magnetization intensity curve of 300 Hz and 400 Hz is plotted, respectively, as shown in Figure 7. When the frequency is 300, the magnetic field intensity is calculated in the case of magnetic field saturation, and the rise and fall of the magnetization process are expressed separately, forming the maximum magnetic hysteresis line. Comparing the 300/400 Hz hysteresis line, when the frequency is 300 Hz, the magnetic field intensity is approximately linear with the magnetization intensity in the range of 0 to 125 kA/m, and when the frequency is 400 Hz, the magnetic field intensity is better in the range of 0 to 100 kA/m. In the linear range, the slope of the curve remains almost constant, and these characteristics can be used for the curve-fitting of the magnetization and description of the approximate linear relationship between the magnetic field intensity and the magnetization intensity in the unsaturated state.

(a)

(b)
Through the analysis of the GMM rod magnetization intensity-magnetic field intensity curve, the shape of the hysteresis return line is found to be related to the magnetization frequency. The lower the frequency is, the higher the linearity will be. It can be inferred that the hysteresis ring in the excitation frequency and excitation current range is included in Figure 7, and the lower the frequency is, the smaller the current and the narrower the hysteresis ring will be.
5.2. Linearization of Magnetization Model
Analyzing the characteristics of the hysteresis curve, the magnetic field intensity in the range of 0 to 100 kA/m is approximately a straight line when the frequency is below 400 Hz. According to the Jiles–Atherton model, the magnetic field intensity and magnetization intensity are calculated (see Table 5). By selecting the linear function and using the least square method to fit the curve, we get the approximate expressions of the magnetic field intensity as well as the magnetization intensity.
Since within a certain range, the magnetization curve approximates a linear relationship; thus, the fitted curve is surely a linear function:where and are the constants to be determined.
Since the problem is boiled down to a polynominal problem, from formula (32), it can be concluded that the corresponding normal equation is
Solve equation (33) with the data in Table 5, and it finds and . Substituting and b into formula (32), we get the linearization model as follows:
In equation (34), the relationship between magnetization intensity and magnetic field intensity is linear. Since the excitation force outputs only when the magnetic field rises, the magnetic field intensity and magnetization intensity are single-value functions, as shown in Figure 8(a). Then, according to equation (1), the relationship curve between magnetic field intensity and magnetostriction coefficient is drawn, as shown in Figure 8(b).

(a)

(b)
5.3. Engineering Applications of Magnetization Models
Excitation current is a nonsinusoidal function; it can be unfolded with the Fourier series by taking the first four phases approximation:where is the angle frequency, is the amplitude of the excitation current, and is the time.
Since the response time of the GMM rod can reach 1 ms, the time delay can be ignored. The magnetic field intensity and the magnetostrictive coefficient are considered to be synchronical, and the excitation current and magnetostriction are considered to be frequency-synchronical. And the relation between excitation current and magnetostrictive coefficient is as follows:where is the time constant of the driving coil.
The output force of GMM rod must be applied to the system of mass-damping-spring, and only after vibration transmission can it generate excitation force. The relation between the magnetostrictive coefficient and the excitation force is as follows:
According to equations (35)–(37), the data in Table 6 are applied to calculate the excitation force corresponding to the excitation current and draw the relationship curve of the two, as is shown in Figure 9(a). Obviously, when the excitation current is a discrete point, the excitation force is also discrete, but their variation trend is consistent and their relationship is approximately linear.

(a)

(b)
Applying formulas (35)–(37) to the open-loop control of GMV, we can plot excitation force curves under different excitation frequencies and excitation currents, as shown in Figure 9(b). The general trend is that the excitation force increases with the increase of the excitation current; however, the degree of linearity is different. When the excitation current is below 5 A, the changing rate of the excitation force is relatively small; especially when the current is in the range of 1.0 A to 1.2 A, the excitation forces are 51 N and 130 N, respectively, and the weak excitation force is due to the insensitivity of the GMM to the change of low magnetic field intensity and due to the existence of nonlinear dead zone [24].
When the current is increased to more than 5 A, the changing rate of the excitation force is large, and when the vibration current is below 20 A, a good linearity can be kept. When the excitation current is above 20 A, the changing rate of the vibration force is different from that in the range of 5~20 A, showing the characteristics of local linearization. However, generally, when the excitation current is in the range of 2.5~30 A, the magnetic field intensity is in the range of 1.5~31.5 A/m, the corresponding magnetization intensity is 12.9~216.9 A/m, and the excitation force can be adjusted in the range of 0.343~10 kN, which indicates that the linear magnetized model can accurately describe the relationship between magnetic field intensity and magnetization intensity. Compared with the Jiles–Atherton model, the linear magnetized model is easy to be solved and is convenient for engineering application. To sum up, the linear model of vibrator can be built with the combination of the linear magnetization model, the structural dynamic model, and the magnetostrictive model.
6. Conclusion
In this paper, DEA and IGA are adopted, respectively, to identify the five unidentified parameters of the Jiles–Atherton model. After comparison, it is found that the average value of the objective function of DEA is as high as 0.917, the average value of convergence generation is 252, while the average value of the objective function of IGA is only 0.615, and the average of convergence generation is 332, indicating that the performance of IGA is more duplicable and stable. Substituting the identification parameters of the two algorithms into the Jiles–Atherton model, respectively, and the relative error of the excitation force of DEA is found to be as high as 13%, while the relative error of IGA is less than 5%. It shows that the identification result of IGA has a better reliability than that of DEA. According to the identification result of IGA and based on the least square method, the linear function is used for the curve-fitting of magnetic field intensity and magnetization intensity, and a linear magnetization model is established to replace the Jiles–Atherton model. When the excitation current is in the range of 2.5~30 A, the magnetic field intensity is in the range of 1.5~31.5 A/m, the magnetization intensity is in the range of 12.9~216.9 A/m, and the excitation force can be adjusted in the range of 0.343~10 kN, which shows that the linear model of GMV can be established with the combination of the linear magnetization model, the structural dynamic model, and the magnetostrictive model. The linear magnetization model can be applied to the open-loop control of the vibrator, and it has engineering application value. Based on this, the online identification method [25] can be further explored to expand the application range of parameter identification.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (NSFC) (62063013).