## 1 Introduction  
In recent years, high-compliant and low-cost pneumatic artificial muscles (PAMs) have been widely implemented in rehabilitation systems [1; 2; 3; 4]. PAMs are shortened in the longitudinal direction and enlarged in the radial direction when being inflated, and they will turn back to their initial form when being completely deflated. PAMs act similar to the human muscle, e.g., the longer muscles produce bigger force and vice versa. Furthermore, these pneumatic muscles are also inherently compliant, which makes them suitable for applying in human-robotic systems. In comparison with the motorized actuators, PAMs are lightweight and have a high power-to-weight ratio. In addition to the aforementioned advantages, the PAM-based applications also have inherent drawbacks, such as very high nonlinearity and uncertainty, and slow response in force generation. These drawbacks make it difficult to model and control PAMs. Using a nonlinear mathematical model to describe the nonlinear characteristic of the PAMs is the most common choice of researchers. In 2003, [PERSON] et al. introduced a three-elements model of PAM, which consists of a contractile (force-generating) element, spring element, and damping element in parallel [5]. Using this type of model, [PERSON] et al. developed the sliding mode control (SMC) based on a nonlinear disturbance observer to improve the tracking performance of a single PAM-mass system [6]. A boundary layer augmented SMC and its modified versions have also been developed for both antagonistic configuration of PAMs and robot orthosis actuated by PAMs [4; 7; 8; 9; 10; 11; 12]. However, the procedure to identify this model's parameters remains complicated with at least two separate experiments: one experiment for determining spring (\(K\)) and contractile (\(F\)) coefficients and another experiment for estimating damping (\(B\)) coefficient. Each experiment must be carried out in three steps [6]. Besides, the parameters of the damping (\(B\)) coefficient must be obtained by measuring the load's acceleration, which is very sensitive to external noise. For this reason, it is difficult to obtain the model's parameters with high accuracy. To deal with hysteresis of PAMs, many hysteresis models have been proposed recently, e.g., Maxwell-slip model [13], Prandtl-Ishlinskii model [14], and Preisach model [15]. In these reports, the dynamic characteristic of PAMs was described by an equivalent pressure/length hysteresis model. The obtained models were used in the feedforward term of the cascade position control scheme for hysteresis compensation. The inner loop of the controllers was designed to regulate the inside pressure of the muscles. The outer loops were designed to deal with the nonlinearity of the PAMs characteristic. Both of the loops use PID-based control strategy. Consequently, some authors continued to develop the modified hysteresis model for both single PAM-mass system and PAMs in antagonistic configuration [16; 17]. However, they mainly focused on modelling of PAMs. Only the trajectory-tracking experiments with low frequency, e.g., up to 0.2 Hz, were conducted in literature. Furthermore, enhanced PID control methods, which were most widely used in these studies, could not deal with hysteresis of PAMs. Another common way to identify the model of PAM-based actuator is the grey-box experiment method [18; 19; 20; 21]. In 2015, to deal with uncertain nonlinearity of PAMs, [PERSON] et al. introduced a grey-box experimental model, which consisted of uncertain, unknown, and nonlinear terms. Based on the built-in model, the authors employed a sliding mode control strategy [18] and an integrated intelligent nonlinear control approach [19] for the tracking purpose. The control performance was significantly improved, and the system could track the \(10^{\circ}\) amplitude sinusoidal signal with 1.5 Hz frequency. The grey-box method was also reported by [PERSON] et al. in 2016 [20] and by [PERSON] et al. in 2018 [21].