Jiles-Atherton-Based Hysteresis Identification of Joint Resistant Torque in Active Spacesuit Using SA-PSO Algorithm

To eliminate the influence of spacesuits’ joint resistant torque on the operation of astronauts, an active spacesuit scheme based on the joint-assisted exoskeleton technology is proposed. Firstly, we develop a prototype of the upper limb exoskeleton robot and theoretically analyse the prototype to match astronauts’ motion behavior. Then, the Jiles-Atherton model is adopted to describe the hysteretic characteristic of joint resistant torque. Considering the parameter identification effects in the Jiles-Atherton model and the local optimum problem of the basic PSO (particle swarm optimization) algorithm, a SA(simulated annealing-) PSO algorithm is proposed to identify the Jiles-Atherton model parameters. Compared with the modified PSO algorithm, the convergence rate of the designed SA-PSO algorithm is advanced by 6.25% and 20.29%, and the fitting accuracy is improved by 14.45% and 46.5% for upper limb joint model. Simulation results show that the identified J-A model can show good agreements with the measured experimental data and well predict the unknown joint resistance torque.


Introduction
In 2018, the ISECG (International Space Exploration Coordination Group), composed of 14 space agencies, released the Global Exploration Roadmap 3 rd edition (GER III) [1], which emphasized that human beings should not only return to the moon and establish a long-term manned lunar base but also complete the manned mars and further deep space exploration. Subsequently, the United States took the lead in announcing the Artemis Program [2], which was expected to achieve the manned lunar landing project in 2024 and the sustainable survival or even long-term residence on the moon before 2028. In this manner, the future maneuvering missions of planetary walking, equipment installation, sample collection, material handling, and base construction will put forward higher requirements for the flexibility and mobility of the EVA (extravehicular activity) spacesuit. However, due to the pressure protection, the multilayer EVA spacesuit would produce an obvious resistance torque, which increases energy expenditure, limits joint mobility, and reduces ergonomic performance [3]; some simple daily operations would not be easily completed after wearing the spacesuit. In order to eliminate the influence of joint resistant torque, some structure schemes have been proposed for spacesuits, such as the MCP (mechanical counter pressure) spacesuit [4,5], the power assist elbow [6], the X1 exoskeleton [7], the elbow-assisted exoskeleton [8], and the hard thigh-hip joint [9]. On the basis of those schemes, we have proposed the concept of active spacesuit [10], that is to say, a joint-assisted exoskeleton is directly worn outside spacesuits to enhance the operational capability, as well as assist astronauts to complete various maneuvering missions. Furthermore, it is well known that astronauts mainly rely on the upper limbs to complete the orbital missions, so active spacesuit designed in this paper would provide assistance to astronauts' upper limbs.
Joint resistant torque is not only the main basis for astronauts' operation intensity and fatigue estimation but also an important parameter for the extravehicular mission planning [11]. Therefore, many efforts should be made into obtaining the accurate resistance torque. Considering the joint resistant torque, one of the important features is the hysteretic characteristics due to the energy loss caused by the joint friction and material elastic distortion [12], so it is a highly challenging task for modeling and predicting the resistant torque. Currently, polynomial fitting based on the measured data can not reasonably describe the hysteresis. In [13], the Preisach model was firstly employed to describe the resistant torque of the EMU (extravehicular mobility unit) spacesuit. In [14], the Preisach model was further used to simulate the hysteretic characteristics and predict the joint resistant torque. Nonetheless, the selection of the Preisach model weight function and the experimental measurement of related parameters are complex, which brings a lot of trouble to numerical calculation and computer programming. In view of the above defects, the J-A (Jiles-Atherton) model [15], which has fewer parameters and less calculation cost, is adopted to describe the hysteretic characteristic of joint resistant torque in this paper. In addition, for the conventional J-A model, the model parameters are not easy to be identified directly by using the common identification methods. Therefore, the parameter identification methods based on intelligent algorithm are addressed as the primary goal.
PSO (particle swarm optimization), which was introduced by Kennedy and Eberhart [16], starts from the random solution in n-dimensional search space, the optimal solution is found through updating generations, so it is suitable for the parameter identification of the J-A model [17,18]. For example, in [19], the authors compared the identification results of PSO, DSM (direct search method), and GA (genetic algorithm) used in the J-A model parameters for two magnetic materials and found that PSO was proved to be better in calculation accuracy and convergence performance. In [20], a modified PSO was proposed to identify the parameters of the J-A model and obtained better parametric solutions than the basic PSO. However, similar to other EAs (evolutionary algorithms), the application of PSO also faces challenges in terms of diversity and convergence. Motivated by some existing works [21][22][23][24], combining with other methods has been extensively investigated, such as C-(coevolutionary-) PSO [25], AGLD-(adaptive granularity learning distributed-) PSO [26], TA-(triple archives-) PSO [27], and CL-(comprehensive learning-) PSO [28]. Considering engineering applications, a SA-(simulated annealing-) PSO algorithm is proposed in this paper. Therein, the SA can accept the worse solutions with a certain probability according to Metropolis criterion and ensure the ability to reduce the chance of getting trapped in a local optimum [29]. By using the random value to update of the global optimal solution, the parameters of the J-A model identified using the SA-PSO algorithm can describe the hysteretic characteristic of joint resistant torque in active spacesuit more accurately. To verify the effectiveness of the proposed SA-PSO algorithm, the simulation comparison test with the modified PSO has been performed in terms of convergence rate and fitting accuracy.
The main content of the paper is as follows: Section 2 introduced the upper limb exoskeleton robot designed in active spacesuit. The Jiles-Atherton hysteresis model adopted to describe the hysteretic characteristic of joint resistant torque is introduced in Section 3. In Section 4, the identification procedure of the proposed SA-PSO algorithm is presented. The comparison of model fitting accuracy and prediction result is given in Section 5 to verify the effectiveness of the identified J-A model. Finally, the conclusions about the identified J-A model for joint resistant torque are discussed in Section 6. The contributions are summarized as follows: (1) A prototype of the upper limb exoskeleton robot is developed and theoretically analyse the prototype to match the astronaut's kinematic behavior

Upper Limb Exoskeleton Robot
Since the real spacesuits are high-value equipment, a mock spacesuit, which is composed of the airtight dry diving suit and the simulated spacesuit, is necessary to be developed. And then, the upper limb exoskeleton robot is directly worn outside the mock spacesuit to achieve the integration or shedding of the existing equipment, as well as enhance the operational capability. In addition, the design parameters of the exoskeleton robot are determined by the size and motion range of the mock spacesuit, including the connecting rod size and the DoFs (degrees of freedom) of joint rotation. Therefore, one side of the exoskeleton robot arm is simplified as a 4-DoF system (2 DoFs-shoulder joint, 1 DoF-elbow joint, and 1 DoF-wrist joint). Due to flexion/ extension in the sagittal plane being the main form of upper limbs' motion, and the additional metabolic rate penalty increases dramatically when mass is added distally during the astronauts' movement [30], we only choose to install the driving motor on the flexion/extension direction of the shoulder and elbow joint in the sagittal plane, so that realize the joints assistance. The structural design of the upper limb exoskeleton robot is shown in Figure 1. The overall structure is composed of base frame structure, arm structure, and some arched fixed protective gears. In the base frame structure, the DoFs of the curved shoulder fixing rod and the curved back rod is convenient for astronauts to put on and remove the exoskeleton robot, and the DoFs of the curved back rod and the curved waist fixing rod is used to adjust the position of the base frame. Moreover, the arm structure is connected with the shoulder support of the base frame structure, and the DoFs can ensure the abduction/adduction of the spacesuit's shoulder joint. In the arm structure, the U-shaped ring is designed to lift heavy objects and fix spacesuit gloves, and arched protective gears can reduce the stress exerted on the upper limbs by the exoskeleton robot. It is worth mentioning that the rods are 2 International Journal of Aerospace Engineering designed to be replaceable and curved, which are convenient to adapt to different EVA spacesuits.

Resistant Torque Modeling
3.1. Spacesuit Resistant Torque. Spacesuit resistant torque, which needs to be overcome in the flexion or extension of the EVA spacesuits, is produced by the soft material structure effect, volume effect, and pressure effect during the bending process [31]. Among these, the structure effect is caused by the stretching, squeezing, and friction of the spacesuits' soft material, and the volume effect and pressure effect are caused by changes in the internal volume and pressure of the spacesuits [32]. Based on the previous work [10], we have measured the spacesuit resistant torque of the mock spacesuit by a single-joint measuring device, and the angle range of joints is needed be determined in advance. When the motion range of the shoulder joint and the elbow joint is set as 0.35-1.22 rad and 0.52-1.57 rad, the first-order transition curves, which are defined as bending from the minimum angle to different maximum angles (0.52 rad, 0.7 rad, 0.87 rad, and 1.05 rad and 0.87 rad, 1.05 rad, 1.22 rad, and 1.4 rad) and then stretching back to the same minimum angle (0.35 rad and 0.52 rad), are shown in Figure 2. It can be seen from Figure 2 that the joint resistant torque has hysteretic characteristics. When the joint angle is the same, the joint resistant torque is different in different directions of the joint rotation. Notably, it is necessary to obtain the motion range of the upper limb joint in advance, before manually measuring the resistant torque corresponding to the mock spacesuit.

J-A Hysteresis
Model. The J-A hysteresis model is derived from ferromagnetic magnetization theory to show the relationship between the applied magnetic field intensity H and magnetization M [33]. Usually, the magnetization M is decomposed into the reversible magnetization M rev and the irreversible magnetization M irr , M irr satisfies the differential equation: where k is the irreversible loss parameter; δis the direction coefficient, when dH/dt > 0, δ = 1; when dH/dt < 0, δ = −1, and H e is the effective magnetic field, which is defined as H e = H + αM, α is the domain wall interaction parameter, and M an is the anhysteretic magnetization, which can be provided by the Langevin function: where a is the shape parameter, M s is the saturation magnetization.
The reversible magnetization M rev satisfies the differential equation: where c is the reversible coefficient. According to equations (1) and (3), the differential equation of magnetization M with magnetic field intensity H can be expressed as where the parameters α, a, c, k, and M s have a clear physical definitions, only one first-order differential equation requires less memory storage, which makes it suitable to describe the hysteretic characteristics of the joint resistant torque.  3 International Journal of Aerospace Engineering which depends on the optimization function. Corresponding to the identified parameters α, a, c, k, and M s of the J-A model, a 5-dimensional search space is selected in this paper. Accordingly, the initial position, velocity and fitness value of the ith particle can be represented by

SA-PSO Algorithm
where P max , P min is set as the maximum and minimum value of the particle position and V max , V min is set as the maximum and minimum value of the particle velocity, T c,j is the jth calculated resistant torque, T m,j is the jth measured resistant torque, T m,max is the maximum measured value, and N is the number of raw data to be compared. Each particle updates itself by following two values in k -iteration: (1) the personal best particle P best , which is defined as the optimal solution found by the personal best fitness F Pbest and (2) the global best particle G best , which is the overall optimal solution found by the global best fitness F Gbest . The velocity and position of particles can follow the motion equation as where ω is the inertia weight; c 1 and c 2 are the learning fac-tors, c 1 = c 2 = 2; and r 1 and r 2 are two random numbers being generated in ½0, 1.
It is well known that the parameters of PSO significantly affect its computational behavior, so the selection of its parameters is crucial [34]. In this paper, a negative tangent curve is selected to control the change of the inertia weight [29]. Specifically, in the initial stage, a large value with a gradual decline rate is set to give sufficient time for the particles to conduct a global search and reduce the probability of falling into the local optimum. In the mid-term stage, the local search ability is gradually strengthened through the approximate linear descent of the inertia weight. When reduced to a small value with a gradual decline rate in the later stages, the algorithm focuses on a detailed local search to obtain the global optimal solution. The inertia weight function can be expressed as where ω max and ω min are the maximum and minimum values, and ω max = 0:95, ω min = 0:4. K is the maximum iteration number, K = 100 in this paper.

SA Algorithm.
For the basic PSO algorithm, it is easy to get trapped in the local optimum because particles tend to be homogenized. To solve the problem, the Metropolis criterion of the SA algorithm is introduced into each iteration of PSO. In other words, the worse solution can be accepted at a certain probability in the process of temperature drop, so as to reduce the probability of falling into the local optimal solution. When the temperature is low, the probability of accepting the worse solution becomes lower, which makes the candidate solution more optimized [35].   International Journal of Aerospace Engineering Start Initialize particles' position P i and velocity V i using (5), (6) Calculate the fitness value of each particle using (7), and record P best , G best , F Pbest , F Gbest .
Update the inertia weight using (10), and the learning factors c 1 and c 2 using (11), (12) Set the temperature drop function T(k) using (11) Calculate the probability p i (k) using (14) p i (k) ≥ r 3 Replace G best Y Update the particles' velocity V i,k and position P i,k , and check their boundary conditions N Evaluate the fitness value and decide whether to updateP best , G best , F Pbest , and F Gbest .

F(Pi, k)<10-3 or k>K
Output the current global optimal particle G best

International Journal of Aerospace Engineering
Firstly, the initial temperature of the proposed SA algorithm is set according to the global best fitness F Gbest and attenuates with a certain cooling coefficient after each iteration.
where μ is the cooling coefficient, μ = 0:95 in this paper. Calculate the probability of accepting the new solution.
where I is the number of particle swarms, I = 50.
Compare the sum of the first m terms of the probability p i ðkÞ with the random number r 3 (value range ½0, 1) to determine whether the global best particle G best is replaced by the new particle. When the new particle is the global optimal particle and increase the probability of the basic PSO algorithm jumping out of the local optimal solution, the combination of PSO algorithm and SA algorithm can effectively reduce the local optimum problem.  Figure 3, and the optimization process is as follows.
Step 1. Set the boundary values of search space and search speed and set the particle swarm size and the maximum iteration.
Step 2. Randomly generate the initial position P i and initial velocity V i of the ith particle according to equations (5) and (6).
Step 3. Calculate the fitness value of each particle using equation (7) and record the personal best particle P best , the global best particle G best , the personal best fitness value F Pbest , the global best fitness value F Gbest .
Step 4. Judge whether the fitness value is less than 10 -3 ; if not, continue with Step 5.
Step 6. Set the initial temperature and the temperature drop function TðkÞ in the kth iterative optimization according to equation (11).
Step 7. Calculate the probability p i ðkÞ of accepting the new solution in the kth iterative optimization according to equation (12).
Step 8. Compare the sum of the first m terms of the probability p i ðkÞ with the random number r 3 to determine whether the global best particle G best is replaced by the new particle.
Step 9. Update the position P i,k and velocity V i,k of the moving particle in the kth iterative optimization according to equations (8) and (9) and check their boundary condi- Step 10. Calculate the fitness value of the moving particle using equation (7) and decide whether to update the personal best particle P best , the global best particle G best , the personal best fitness value F Pbest , and the global best fitness value F Gbest .
Step 11. Judge whether the fitness value is less than 10 -3 or the maximum number of iterations K is reached, If not, k = k + 1, and update the temperature Tðk + 1Þ = TðkÞ * μ, then go to Step 5.
Step 12. Output the current optimal particle, that is, the optimization result, and the algorithm terminates.

Modified PSO Algorithm.
In this section, to verify the validity of the SA-PSO algorithm proposed in this paper, the MPSO (modified PSO) algorithm, which is developed for the SMA-based compliant actuator platform based on the J-A model [20], is also applied for the parameter identification in the test. In MPSO, an iteration-varying inertia weight ΛðnÞ is updated by where n is the algorithm iterations; Λ max , Λ min are the selected maximum and minimum values, Λ max = 0:9, Λ min = 0:4; and Ξ is the maximum iteration number, Ξ = 100 in this paper.    The cognitive parameter c 1 and the social parameter c 2 can be expressed as where c 1ini , c 2ini , c 1end , and c 2end are positive constants, c 1ini = 2:5, c 2ini = 0:5, c 1end = 1:3, and c 2end = 1:7.
In addition, we would select the partial resistant torque of the first-order transition curves in Figure 2, the different maximum angles of the shoulder joint and the elbow joint is set as 0.52 rad, 0.87 rad, and 1.22 rad and 0.87 rad, 1.22 rad, and 1.57 rad, respectively. Comparing with the measured resistant torque, the identification results of applying the SA-PSO and MPSO algorithms for the J-A model are given and illustrated. An appropriate search space is obtained according to the literature [11], as shown in Table 1.

Comparison of Identification Results.
In the SA-PSO algorithm and the MPSO algorithm, the size of the particle swarm is selected as 50 and the maximum iteration number as 100. Following the procedure introduced in Figure 3, the   Table 2, and the fitness values of the SA-PSO and MPSO algorithms are also given in Figure 4.
As for the shoulder joint model, in terms of convergence rate, the SA-PSO algorithm starts to converge at the 48th iteration, while the MPSO algorithm is the 51th iteration, which is advanced by 6.25%. In terms of fitting accuracy, the optimal fitness result of the SA-PSO algorithm is 1.5473, while the MPSO algorithm is 1.7709, which is improved by 14.45%. For the parameter identifi-cation results of elbow joint model, the convergence rate is advanced by 20.29% (the iterative times for convergence of the SA-PSO and the MPSO are 69 and 83, respectively), and the fitting accuracy is improved by 46.5% (the optimal fitness result of the MPSO and the SA-PSO are 1.7243 and 1.1770, respectively). It can be seen that the SA-PSO algorithm is superior to the MPSO algorithm in convergence rate and fitting accuracy, the comparison of the experimental data of measured resistant torque and the identified the J-A model of upper limb joint is shown in Figures 5 and 6.    It can be seen from Figures 5 and 6 that the fitting effect is poor in some areas, mainly due to the manufacturing error of spacesuit joints and the measurement error of resistant torque. As a result, the experimental data of resistant torque does not show good symmetry, which makes the fitting error of the symmetrical J-A model large. However, it can be observed that the optimal fitness results is 1.5473 of shoulder joint and 1.1770 of elbow joint, respectively, which is proved that the simulation results can show good agreements with the measured experimental data, so that the identified J-A hysteresis model can be applied to the research fields of joint motion performance and ergonomic performance.

Prediction
Results. When given a known motion range of upper joints, we can obtain the actual resistant torque through the measuring device. In other words, the resistance torque at a specific angle needs to be tested in advance. However, the joint motion is unknown and uncertain under different operating conditions, so it is desirable to apply the identified J-A model to predict the resistant torque of upper limb joint. In this manner, we would select the partial resistant torque of the first-order transition curves in Figure 2, the different maximum angles of the shoulder joint and the elbow joint is set as 0.7 rad, 1.05 rad and 1.05 rad, 1.4 rad. Comparing with the partial measured resistant torque, the prediction results of applying the SA-PSO and MPSO algorithms for the J-A model are shown in Figure 7 and the modeling errors are shown in Figure 8.
It can be seen from Figure 8 that the SA-PSO algorithm is superior to the MPSO algorithm in prediction effect of unknown resistance torque. For two angles of motion for the shoulder joint, the maximum error of the MPSO algorithm and SA-PSO algorithm is 1.84 Nm, 1.39 Nm and 2.95 Nm, 0.41 Nm, and the prediction error of the SA-PSO algorithm is relatively stable in the motion range of 0.35-1.05 rad. As for the elbow joint, the overall prediction errors of the SA-PSO algorithm is less than that of the MPSO algorithm (the maximum error of the MPSO and SA-PSO algorithm is 0.91 Nm, 1.84 Nm and 1.08 Nm, 0.79 Nm). In addition, the prediction results are proved that the identified J-A hysteresis model can well predict the unknown resistance torque.

Conclusions
In this paper, the parameter identification of spacesuits' joint resistant torque described by the J-A model is addressed. The experimental data of joint resistance torque are limited and need to be measured in advance, so it is desirable to apply the identified J-A model to model and predict the resistant torque of upper limb joint. In addition, the accurate prediction of these joint resistant torques in the upper limb exoskeleton robot is conducive to the application of an active spacesuit. For such purpose, a SA-PSO algorithm is proposed to identify the J-A model parameters. Simulation results show the SA-PSO algorithm is superior to the MPSO algorithm in convergence rate, fitting accuracy, and prediction effect. For this study, some future work to improve some areas with poor fitting effect can be consid-ered by introducing the compensation control algorithm. Besides, we shall focus on the modeling and predicting of multiorder transition curves of joint resistant torque and apply some effective algorithms to identify the Jiles-Atherton model parameters.

Data Availability
The raw data of the joint resistant torque involved in the research process of this paper are obtained by a measuring device independently developed by our laboratory. The corresponding data can be obtained by contacting the first author (lzy_uestc@hotmail.com).

Conflicts of Interest
The authors declare that they have no conflicts of interest.