High-Precision Displacement and Force Hybrid Modeling of Pneumatic Artiﬁcial Muscle Using 3D PI-NARMAX Model

: Pneumatic artiﬁcial muscle (PAM) is attractive in rehabilitation and biomimetic robots due to its ﬂexibility. However, there exists a strong hysteretic nonlinearity in PAMs and strong coupling between the output displacement and the output force. At present, most commonly used hysteresis models can be treated as two-dimensional models, which only consider the nonlinearity between the input and the output displacement of the PAM without considering the coupling of the output force. As a result, high-precision modeling and estimation of the PAM’s behavior is difﬁcult, especially when the external load of the system varies signiﬁcantly. In this paper, the inﬂuence of the output force on the displacement is experimentally investigated. A three-dimensional model based on the modiﬁed Prandtl–Ishlinskii (MPI) model and the Nonlinear AutoRegressive Moving Average with eXogenous inputs (NARMAX) model is proposed to describe the relationship and couplings among the input, the output displacement, and the output force of the PAM. Experiments are conducted to verify the modeling accuracy of the proposed model when the external load of the PAM varies across a wide range. The experimental results show that the proposed model captures well the hysteresis and couplings of the PAM and can precisely predict the PAM’s behavior.


Introduction
Pneumatic artificial muscle (PAM), a flexible actuator, mainly consists of an inner elastomeric tube and an outer braided mesh [1]. When the inner pressure increases, the PAM will produce radial expansion and axial contraction, similar to muscles. Compared with other popular actuators, e.g., motors, PAM features a simple structure, easy installation, light weight, and large output force/weight ratio [2]. Further, PAM shows good flexibility because it generates an output force and displacement based on its elastic deformations [3]. This helps to guarantee the safety in human-machine interactions. Therefore, it has been widely used in rehabilitation and biomimetic robots [4]. For example, Festo developed different bionic manipulators using PAMs. A PAM-driven intrinsically soft ankle rehabilitation robot was developed by Meng et al. [5]. A novel bionic elbow joint system actuated by PAM was developed by Yang et al. [6].
It must be noted that the precise control of the PAM is challenging. PAM exhibits a complex hysteresis nonlinearity due to the internal friction, the deformation of the elastomeric tube during the gas charging and discharging, etc. [6][7][8]. Further, the hysteresis is sensitive to many factors. For instance, the rate dependence refers to the fact that the hysteresis loop will become wider if the input rate increases, especially at higher frequencies. Consequently, the motion accuracy of the PAM is significantly affected, which increases the difficulty in the motion control of PAM-based systems. Therefore, modeling and compensation of the PAM's hysteresis have become some of the hot topics. In the research work of Sofla et al. [9], a PAM-driven multi-section-compliant robotic manipulator was developed, where the Bouc-Wen model was modified based on the constant curvature method and the concentrated masses to compensate the asymmetric hysteresis. For the miniature PAM-based catheter control, a feedforward controller based on a generalized Prandtl-Ishlinskii model was proposed by Shakiba et al. [10]. Liang et al. proposed an energy-based nonlinear control method for a two-link PAM-driven robot to realize accurate positioning control [11].
The hysteresis models can be briefly classified into three categories: the integral type, the differential type, and the system identification type [12]. The integral model, also called the operator-based model, mainly includes the Maxwell model, Preisach model, and Prandtl-Ishlinskii (PI) model. The differential model mainly includes the Duhem model and Bouc-Wen model. In the system identification model, the popular algorithms include the least mean squares, the recursive least squares, the genetic algorithm, the particle swarm optimization, and the neural network. Based on the above basic structure, a variety of hysteresis modeling methods can be obtained. For instance, a rate-dependent PI model was used to improve the control performance of a compliant mechanism [13]. A novel black-box approach was proposed to build the pressure-contraction hysteresis model of the PAM across multiple operating ranges using an adaptive-network-based fuzzy inference system (ANFIS) [14].
The above hysteresis models can be labeled as two-dimensional (2D) models reflecting the mapping from the input to the output displacement of the PAM. For PAMs, due to the flexibility, there exist strong couplings between the output displacement and output force. In other words, the behavior of the PAM will vary significantly at different loads, i.e., the so-called load dependence [15]. As a result, the input and outputs relationship of the PAM is three-dimensional (3D) in nature. However, the effect of the output force is not considered in the above 2D models. Consequently, they cannot precisely predict the PAM's behavior if confronted with varying loads or strong disturbances.
Some solutions have been proposed to account for the influence of the output force (or the external load). In Sofla's research [16], a modified Bouc-Wen model was used to describe the relationship among the inner pressure, the hysteretic restoring force, and the length of the PAM. The influence of the PAM's length on the hysteretic restoring force was also considered. This model was then adopted in the design of a PAM-driven manipulator [9]. In Zhang's work, a comprehensive dynamic model was proposed for predicting the dynamic hysteresis behaviors with rate-dependent and load-dependent effects [17]. In addition, in Konda's work [18], the influence of the load on the PAM, shape memory polymer, and super-coiled polymer was investigated, and the corresponding modeling and compensation methods were developed.
In this paper, the modified Prandtl-Ishlinskii (MPI) model is integrated with the Nonlinear AutoRegressive Moving Average with eXogenous inputs (NARMAX) model to describe the relationship among the input-displacement-force of the PAM. The proposed model is composed of a hysteresis modeling unit and a coupling unit. The hysteresis modeling unit characterizes the hysteresis between the input and output displacement of the PAM. The coupling unit is used to account for the couplings between the output displacement and output force. The effectiveness of the proposed model is tested on an in-house-built testbench for PAMs. The experimental results show that the proposed MPI-NARMAX hysteresis model can describe the complex relationship among the inputdisplacement-force well. Compared to the standalone MPI model individually identified at different loads, it does not need to repeat the identification process for the MPI-NARMAX model if the external load varies significantly. Therefore, the proposed MPI-NARMAX model achieves higher modeling accuracy and applicability.

Testbench for PAM Characterization
In this paper, a testbench is developed to simultaneously record the control input, the output displacement, and the output force of the PAM to facilitate the 3D hysteresis modeling of PAMs. As shown in Figure 1, the PAM is installed on a linear guiding rail system with one end clamped and the other end sliding along the rail. The inner pressure of the PAM is controlled by a proportional pressure regulator valve (model VPPM-6L-1-G18-0L6H from Festo) with the help of an air compressor. The displacement of the PAM is measured using a displacement sensor (Model 102322/A from novotechnik). A force sensor (Model H3-C3-200 kg-3B from ZEMIC) is used to measure the external force applied on the PAM or the payload of the system. The data acquisition is implemented in the environment of a real-time target machine (Model Mobile from Speedgoat). The overall system runs at a sampling rate of 1 kHz. Different PAMs can be installed on the testbench for the characterization and testing.

The Input-Displacement Hysteresis of the PAM
As previously mentioned, there is an obvious hysteresis nonlinearity between the input and the displacement of the PAM. In this paper, the input voltage applied to the pressure regulator valve is adopted as the input of the overall system. First, the inputdisplacement hysteresis is investigated by removing the external load from the PAM. In this case, the output force can be treated as 0. The following sinusoidal signal is adopted as the input: where u represents the input, f is the frequency of the signal, and A and δ are the maximum allowable voltage and the dead zone of the pressure regulator valve, respectively. For the selected pressure regular valve, A = 10 and δ = 0.1. In human-machine interactions, the PAM generally works at a relatively slow speed. As a result, in the preliminary tests, the frequency of the signal is set to f = 0.05, 0.1, and 0.2 Hz. The displacements of the PAM are measured, and the input-displacement hysteresis loops are shown in Figure 1c. It can be easily found that the PAM's hysteresis loop is not symmetric about its loop center. As the frequencies are low in Figure 1c, the rate dependence of the PAM's hysteresis is observable but not very obvious.

The Couplings between the Output Force and Displacement
As the PAM is a flexible actuator, its displacement and force are coupled. The variation in the external load or disturbances will significantly affect the PAM's displacement [8]. In order to intuitively demonstrate the couplings between the output force and displacement, external forces are added to the free end of the PAM when the PAM moves back and forth. First, all the load is removed from the PAM. In this case, the measured hysteresis all comes from the PAM itself. Secondly, a constant load (~91 N) is added to the hook. Finally, the hook is pulled by hand to provide a varying load. Figure 2 shows the PAM's behaviors in the above scenarios at different frequencies. It can be found from Figure 2 that the external load will significantly affect the behavior of the PAM. For the constant load, the stroke of the PAM will decrease, i.e., the hysteresis loop of the PAM will be squeezed along the Y axis. It is noted that the shape of the hysteresis loop is basically unchanged. However, when the external load is varying during the process, there exist several regional distortions on the measured hysteresis loops. As a result, it will be very difficult to model such a distorted hysteresis loop using 2D hysteresis models.
The couplings between the output force and displacement are also experimentally investigated. In this case, the inner pressure is kept constant at different levels, and the relationship between the external load and the displacement of the PAM is measured. The experimental results are plotted in Figure 3. It is observed that hysteresis loops can also be found between the output force and displacement of the PAM. Based on the above results, it is clear that there exist a strong hysteresis and couplings among the input, the output force, and the displacement of the PAM. The 3D hysteresis model is required to obtain an accurate model of such a complex system.

MPI-NARMAX Hysteresis Model
In this paper, a 3D hysteresis model, namely the MPI-NARMAX model, is proposed for modeling the complex relationship among the input-displacement-force of the PAM. The schematic diagram of the model is illustrated in Figure 4, where the current input voltage, the current external force, and the historical output displacement are adopted as inputs. The proposed model can be divided into a hysteresis modeling unit and a coupling unit. In the hysteresis modeling unit, the MPI model is used to characterize the hysteresis between the input voltage and the output displacement of the PAM when no external load is applied. Subsequently, in the coupling unit, the NARMAX model is used to model the couplings between the output force and displacement of the PAM. Compared with the 2D hysteresis models, the proposed 3D hysteresis model is developed to predict the PAM's behavior even when confronted with varying external loads.

The MPI Model
The classical PI model is a widely used hysteresis model consisting of a set of weighted backlash operators. The backlash operator is defined as follows: where u represents the input of the operator, h i represents the output of the ith backlash operator, r i is the threshold, and T represents the sampling period. Then, the classical PI model can be obtained by the weighted summation of a series of backlash operators: where n is defined as the order of the model, ω i is the weight for h i , and b(t) represents the output of the classical PI model. The hysteresis of PAM is rate-dependent, whereas the classical PI model is static. In order to account for the rate dependence, the weights of backlash operators can be linearly adjusted according to the input rate, which can be described as follows: where k i is the slope and d i is the intercept. From Figure 1c, it can be observed that the hysteresis loop of the PAM is asymmetric. However, the classical PI model is a symmetric model, making it difficult to describe the asymmetry of the PAM's hysteresis well. Based on our previous work, a polynomial operator can capture well the asymmetry and improve the modeling accuracy [13]. Therefore, the following polynomial operator is integrated to form the MPI model: where f (t) represents the polynomial output, a i represents the polynomial coefficient, and b i (t) represents the ith power of the output of the classical PI model in (3).

NARMAX Model Based on RFNN
The NARMAX model can be used to deal with the identification problem of nonlinear systems [19]. A nonlinear system can be described as: where y(t) and x(t) represent the output and input of the NARMAX model, respectively, n y and n x represent maximum time lags associated with y(t) and x(t), respectively, and g(•) represents the nonlinear relationship between the outputs and inputs. For the hysteretic system in this paper, the NARMAX model can be described as: whereŷ(t) represents the predictive output of the system, f (t) represents the output of the MPI model, F(t) is the output force of the PAM, and n a , n b , and n c represent the maximum orders inŷ(t), f (t), and F(t), respectively. Based on previous research [20,21], g(•), the nonlinear mapping relationship between y(t), f (t), and F(t), can be identified by neural networks. Therefore, the recurrent fuzzy neural network (RFNN) is selected in this paper to identify the NARMAX model. Generalized from the fuzzy neural network (FNN), the RFNN was shown to possess the same advantages over recurrent neural networks and extend the application domain of the FNN to temporal problems [22]. In this paper, the structure of the RFNN is shown in Figure 5. For the convenience of subsequent expression, N is defined as: There are 5 layers in the RFNN. The first layer includes 3N neurons, i.e., 3 fuzzy operators for each input. It realizes time recursion and the fuzzification process. The output of each node can be described as: where x i represents the ith input element, j represents the ordinal number of the fuzzy operator, m ij and σ ij are parameters of the Gaussian membership function, and θ ij is the parameter of time recursion. The second layer, which contains 3 N neurons, combines the outputs of the first layer to achieve the nonlinearity of the output: The third layer, which has the same number of neurons as the second layer, normalizes the output of each node in the previous layer: The fourth layer multiplies the output of the previous layer as a weight on the input: Finally, the fifth layer sums the outputs of the fourth layer:

Parameter Identification of MPI Model
In the MPI model, the parameters that need to be identified include r i , k i , d i , and a i . The least squares method is used to identify the MPI model from the experimental results of the 0.1 Hz sinusoidal excitation without load. The identified parameters are given in Table 1. The measured displacement and model output are shown in Figure 6. It can be observed that great agreement between the measurement and the model output has been achieved, as shown in Figure 6a,c. In order to quantify the modeling accuracy, the mean absolute error (MAE) and root-mean-square error (RSME) are calculated to be 0.4538 mm and 0.6120 mm, respectively.  The statistics of the modeling error are shown in Figure 6b. It can be observed that the modeling errors center around 0 and lie mainly within the range of ±1 mm. As a result, the MPI model can well predict the PAM's displacement when no load is added to the system. However, the influence of the external load is not considered in MPI model. Typically, if a constant load is added to the hook, the same modeling and identification process has to be repeated to obtain another MPI model. This is straightforward and thus is not be presented herein, to guarantee the conciseness of the paper.

Training of RFNN
In order to capture the influence of the time-varying load on the PAM's displacement, different external loads are applied on the hook when the PAM is moving. The output force and displacement of PAM are measured during the process. The output of the MPI model, the output force, and the displacement of the PAM are used as the training data for the RFNN.
In the training of the RFNN, the standard particle swarm optimization (PSO) method is used to identify the parameters of the RFNN. PSO is initialized by a group of particles, and the optimal solution is found in the subsequent iterative update. The updating rule of particles can be expressed by the following formulas: where q i (t) and v i (t) represent the current position and velocity of the particle, respectively, w v is the inertial weight, c 1 and c 2 represent the learning constants, rand() represents the random number, and pbest and gbest represent the best position for the individual particle and all particles, respectively. The linear decreasing weight is used to update the inertial weight w v . The updating rule of w v can be expressed by the following formula: where m represents the current iteration number. In this paper, the learning constant c 1 and c 2 are set to be 2. w ini and w end , which represent the initial and final value of w v , respectively, are set to 0.9 and 0.4. The overall 3D hysteresis model can be obtained after training the RFNN and cascading it to the MPI model. The experimental results of the training set are shown in Figure 7. The external load with a peak of approximately 100 N is applied and removed from the hook three times. The shaded areas in Figure 7a indicate the durations of the external load. In this manner, the PAM's behaviors with and without external load can be captured in one measurement. It can be found from Figure 7 that when there is no external load, the MPI-NARMAX model can follow the displacement of the PAM well. When there is external load on the PAM, the MPI-NARMAX model can still predict the output displacement of PAM well. There is no obvious change in the magnitude of the modeling error for the PAM with and without external loads. The MAE and RMSE are calculated to be 0.8730 mm and 1.0913 mm, respectively. Due to the influence of the external load, the inputdisplacement hysteresis loops are distorted, as shown in Figure 7b. This will be difficult to model using conventional 2D hysteresis models. The couplings between the force and displacement of the PAM are shown in Figure 7c. It can be observed that the proposed MPI-NARMAX successfully captures the influence of the external load and improves the prediction accuracy on the output displacement of the PAM. As shown in Figure 7b, the input-displacement relationship of the MPI-NARMAX model matches the measurement well.  Similar to the training set, the external load is applied and removed from the hook three times. However, the peak of the external load increases to approximately 150 N. For such a high and varying external force, the proposed model can still precisely predict the output displacement. The MAE and RSME are calculated to be 0.9304 mm and 1.2138 mm, respectively. The modeling error on the testing set is comparable to the modeling error on the training set. The input-displacement hysteresis loops and force-displacement couplings are shown in Figure 8b,c, respectively. Good agreements between the measurements and the model output can be found. This demonstrates that the complex relationship among the input, the output displacement, and the output force of the PAM can be captured well in the proposed 3D hysteresis model.
In order to compare the performance of the proposed model, the performance of the MPI model identified in Section 4.1 is also provided in Figure 8. As this standalone MPI model is identified without any external load, the influence of the external load is totally ignored in this MPI model. As a result, there exist obvious discrepancies between the MPI's output and the measurement. The modeling accuracy of the standalone MPI model is poor. The MAE and RMSE of the MPI model are calculated to be 3.8926 mm and 4.4845 mm, respectively. This demonstrates that the conventional 2D hysteresis models cannot capture the couplings between the output force and displacement of the PAM. The modeling accuracy of the 2D models are poor when confronted with strong and varying external loads.

Cross-Checking with Different Load Statuses
In applications, the influence of the external load can also be modeled by applying different payloads to the PAM and identifying the parameters of the standalone MPI model for each payload. As a result, instead of a single model, multiple MPI models are obtained using this method. In order to compare the performance of the proposed MPI-NARMAX model with these models, the PAM's behaviors in the following load statuses are measured: In status 4, it is very difficult to use a standalone MPI model to fit the measurement, as the measured hysteresis loops contain several regional distortions. As a result, three MPI models are individually identified at status 1, 2, and 3, which are denoted as MPI-1, MPI-2, and MPI-3, respectively.
The identified MPI models, together with the proposed MPI-NARMAX model, are cross-checked using the other measurements. The MAE and RMSE of these models with different load statuses are calculated and listed in Tables 2 and 3, respectively. For standalone MPI models, it is observed that a low modeling error is achieved in identification, as marked in bold. As a result, if the load status remains, a high modeling accuracy is expected. However, if the load status changes, the modeling error increases obviously. For the varying load in status 4, the modeling accuracy of all the standalone MPI models is poor.
On the contrary, the modeling accuracy of the proposed MPI-NARMAX model is consistent, regardless of the change in the load status. For constant loads, the modeling error of the proposed MPI-NARMAX model is comparable to the standalone MPI models individually identified. For the varying load in status 4, the modeling accuracy of the proposed model is the best among the models.
The above experimental results also demonstrate that the proposed MPI-NARMAX model is effective for both constant and varying external loads. It is superior to the standalone MPI models individually identified at different constant loads as it can adapt to different load statuses after only one training, and the parameter identification process does not need to be repeated when the external load changes.

Conclusions
For PAMs, there exists strong hysteresis between the input and outputs. In addition, the coupling between the output displacement and output force is obvious. These make the precise modeling of the PAM's behavior very difficult, especially when varying external load or disturbances occur. This paper proposes a displacement and force hybrid 3D hysteresis model, where the MPI model and RFNN-based NARMAX model are integrated. The MPI model is constructed by cascading a polynomial saturation operator to the classical PI model to account for the asymmetry of the PAM's hysteresis. The RFNN-based NARMAX model is used to capture the couplings between the output force and displacement of the PAM. Similar research work can be found in Zhang's paper [17], where the weight of the constant payload is included in the training data, making the model capable of predicting the behavior of the PAM with different but constant loads. In this paper, with the help of a force sensor, the dynamic variation in the external load can be measured and included as one input to the model. This makes the proposed MPI-NARMAX hysteresis model capable of predicting the PAM's behavior regardless of the change in the external load.
The proposed MPI-NARMAX model only needs to be identified and trained once. The experimental results of the PAM with constant and varying external loads verify that the proposed MPI-NARMAX hysteresis model successfully captures the complex hysteresis and couplings among the input, the output displacement, and the output force of the PAM. It is capable of precisely predicting the output of the PAM when confronted with constant and varying external loads. The modeling accuracy is robust against the variation in the external load. The effectiveness and applicability of the proposed MPI-NARMAX model are verified.