Discrete-Time Fractional Order Integral Sliding Mode Control of an Antagonistic Actuator Driven by Pneumatic Artiﬁcial Muscles

: Recently, pneumatic artiﬁcial muscles (PAMs), a lightweight and high-compliant actuator, have been increasingly used in assistive rehabilitation robots. PAM-based applications must overcome two inherent drawbacks. The ﬁrst is the nonlinearity due to the compressibility of the air, and the second is the hysteresis due to its geometric construction. Because of these drawbacks, it is difﬁcult to construct not only an accurate mathematical model but also a high-performance control scheme. In this paper, the discrete-time fractional order integral sliding mode control approach is investigated to deal with the drawbacks of PAMs. First, a discrete-time second order plus dead time mathematical model is chosen to approximate the characteristics of PAMs in the antagonistic conﬁguration. Then, the fractional order integral sliding mode control approach is employed together with a disturbance observer to improve the trajectory tracking performance. The effectiveness of the proposed control method is veriﬁed in multi-scenario experiments using a physical actuator.


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, D. B. Reynolds 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, K. Xing 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]. control systems. In addition, the proposed is designed in a discrete-time domain, therefore, it can be easily implemented by any digital control system.

Experiment Platform
A typical configuration of antagonistic configuration of PAMs is shown in Figure 1a, and the proposed experiment platform is demonstrated in Figure 1b. The experimental system consisted of two PAMs which had 1.0 inches of diameter and 22 inches of length. The PAMs were fabricated at our local institute. The pressure inside each PAMs was regulated by two proportional electric control valves series ITV 2030-212S-X26 from SMC company. One potentiometer CP-20H from Midori Precision (Japan) was used to measure the actuator's angle. All the control systems were implemented by using computer-based controller NI cDAQ-9178 from National Instrument (USA). The real-time controller collected the data from the potentiometer via analogue input module and sends the control signal to the electric control valve via analogue output modules. The developed control algorithm was implemented and compiled by LabView software before downloading it to the hardware controller.

System Modeling
Based on the geometry of the typical antagonistic configuration, which is illustrated in Figure 1a, the length of each PAMs can be obtained from the measured joint angle, as in the following equations: where y AN and y PN are the nominal length of the anterior and posterior PAMs when the joint angle θ = 0, and R is the rotation radius of the actuator. Because two similar PAMs are used in the system, we can consider that y AN = y PN = y N . Following that, the relationship between contraction of PAMs and the measured angle can be expressed as where y 0 is the length of PAMs in the complete deflation state. In (2), y 0 and y N are fixed by the deflation and nominal lengths of PAMs. Therefore, the contraction of PAMs can be expressed as the function of the measured joint angle θ. As a result, the dynamic behaviour of an antagonistic muscle can be described by a single input single output (SISO) system, in which the input is the difference pressure of two PAMs (∆P), and the output is the measured angle θ. The input pressure inside the anterior and posterior PAMs can be expressed as where P 0 is the nominal pressure which determines the initial position of antagonistic actuator. The nominal pressure can be chosen so that the joint has the desired compliance for a specific application. It was fixed, so ∆P was chosen as a control variable of trajectory-tracking controller.
All the system parameters P 0 , y 0 , and y N are provided as in Table 1. In this research, the following discrete-time SISO system was chosen to describe the model of antagonistic actuator: where u k represents the control pressure ∆P, y k is the joint angle, d is a positive integer representing the dead time of the system (as a number of the sampling time), p k is the unknown disturbance of the system, a i and b j are the model parameters with b 1 = 0, n and m are integers which satisfy n ≤ m.
The model parameters of the system are obtained by the identification experiment. To verify the mathematical model of PAM, the following experiment procedure was carried out.
Step 1: The initial position of the actuator was set at 0 • by supplying nominal pressure P 0 to each PAMs of the actuator.
Step 2: The actuator angle can be changed by sending different types of control signal to the electrical control valves. Three types of control signals were used in this experiment:

•
Step response: the control signal was a step wave with the final values 0.2, 0.4, 0.5, and 0.8 MPa.

•
Sinusoidal signal: The control signal is the 0.2 MPa amplitude sinusoidal signal, where frequency varies from 0.2 to 1.0 Hz. • A sine wave signal with time-varying amplitude and frequency, as in the following equation: where A = 0.05 MPa and f = 0.5 Hz are the basis amplitude and frequency of the control signal, respectively.
All the data, including control signals and measured angles of actuator, were recorded with sampling time T s = 5 ms for further analysis.
Step 3: The discrete-time SOPDT, in which m = n = 2, was chosen as the mathematical model of the actuator with good accuracy. The precise values of the model's parameters are estimated by using the MATLAB software and provided in Table 2.    Table 2. As seen in Table 2, the standard deviations of the model parameters were much smaller than their mean values. Therefore, we can conclude that the model parameters obtained by different methods have similar values. As a result, we can use any aforementioned method for the identification purpose. The model parameters, which were identified by time-varying amplitude and frequency, are chosen to design the controller in Section 3 of this paper.

Control Design
Recently, SMC has been employed for designing the controller for PAMs or systems powered by PAMs [4,[6][7][8][9][10][11]. SMC is able to provide highly accurate tracking performance with a bounded error; however, "chattering" problem is a big challenge that SMC must overcome. SMC is a suitable control approach for PAM-based systems to deal with their uncertain, nonlinear and time varying characteristics. In this research, we addressed a DFISMC to improve the tracking performance of the antagonistic actuator powered by PAMs. The fractional order integral is implemented together with disturbance observer to deal with the "chattering" problem. Figure 3 illustrates the block diagram of the proposed control system. We consider the following fractional integral sliding surface: where e k = y * k − y k is the tracking error with the desired trajectory y * k , and α Ξ e,k is the integral of the tracking error with fractional order α and integral gain K I . α Ξ e,k can be calculated as follows: and α Ξ e,0 = ω N e 0 at the initial state. Please refer to Appendix A for details about fractional integral approximation. We also obtain From (6)- (8), we can obtain Therefore, where e k + 1 is one-step-ahead tracking error, which can be computed from the SISO model of the actuator in (4) as where y * k+1 is one step ahead of the desired trajectory, which is considered to be known when apply the model to a specific application. In (4), disturbance p k is unknown and needs to be estimated. In this study, one-step delayed technique was used to estimate p k . This technique is based on the following assumptions: Assumption 1. Sampling time T s was sufficiently small and system disturbance p k is bounded, so the difference between two consecutive time samples is also bounded, i.e., where O(T s ) is the thickness boundary layer. It means there always exist constants A and B, ∀ k > 0, such that The aforementioned assumption was based on the Taylor expansion described in Appendix B. Estimationp k of disturbance p k can be computed based on (4) aŝ where Hence, the error of estimationp k is Finally, the one-step-ahead tracking error (11) can be expressed by When substituting e k+1 in (11) and p k =p k +p k into (10), we can obtain Disturbance estimation errorp k is unknown in practice; however, it is very small and bounded by assumption 1. Control signal u k can be obtained by solving the reaching law S k+1 = 0 with the absence ofp k as follows: . (21) Adjusting integral gain K I and fractional order integral α may improve performance of the control system.

Experimental Procedure
To verify the effectiveness of the proposed control method, multiple-scenario experiment with different desired trajectories was carried out. In the first scenario of the trajectory-tracking experiment, sinusoidal signals with amplitude 20 • and 0.2, 0.5, 0.8, and 1 Hz frequency were given as desired trajectories. To evaluate the applicability of the proposed control method for a rehabilitation robot, a human-like pattern signal was employed as a desired trajectory in the second scenario of the experiment. The modified knee gait data profile in textbook [29], where the maximum flexion angle is set at 28 • , was used to verify the control performance. In both experimental scenarios, the system was tested under two load conditions: no load and load m = 2.5 kg.
In all experimental scenarios, the sampling time of the discrete-time control system was T s = 5 ms. All the data were recorded for ten cycles from the start time of the experiment. The data were processed by MATLAB software version R2016b. The proposed controller is also compared with the conventional discrete-time sliding mode control (DSMC) method in terms of tracking performance. The parameters of both controllers after being well-tuned are provided in Table 3.

Experiment Result
To quantitatively evaluate the tracking performance, the maximum tracking error (MTE) and root mean square tracking error (RMSTE) are computed. The RMSTE is calculated as follows: where N is the total number of data samples, and e k is the tracking error under k th sample. Figure 4 depicts the experimental results when the actuator tracks the sine wave signals without load. The sinusoidal signal had an amplitude of 20 • and frequency from 0.2 Hz to 1 Hz. In the second scenario, a knee gait pattern was given as a desired trajectory. The proposed controller was evaluated with two different gait cycle (GC) times: 2.5 s and 4 s. The experimental results of this scenario are shown in Figure 5. In both Figures 4 and 5, the upper sub-figure of each image includes the desired trajectory (back line), measured angle controlled by DFISMC (dash-blue line), and measured angle controller by conventional DSMC (dash-dot red line). The lower part of each figure shows the tracking errors of both proposed controller and DSMC controller. In comparison with the traditional DSMC, the DFISMC was able to provide a better performance in both transient and steady states. As demonstrated in Figure 6, in all scenarios of the experiment, both MTE and RMSTE of the proposed control approach are smaller than the ones of the conventional DSMC control method. For example, when tracking the 1.0 Hz frequency sinusoidal signal, the RMSTEs of the DFISMC and DSMC controllers are 1.63 • and 1.43 • , respectively. It means that DFISMC is able to provide a better performance than the conventional DSMC controller. In particular, as seen in the error graphs in Figures 4 and 5, the finite amplitude oscillation of the tracking error in DFISMC is much smaller than in DSMC. It can be concluded that the inherent "chattering" phenomenon of SMC control is reduced with DFISMC. The numerical values of the experimental results in all scenarios are given in Table 4.     Figures 7 and 8 show the control performances of the system when tracking the sinusoidal signals and human-gait pattern with a load of 2.5 kg, respectively. When the antagonistic actuator carries a load m = 2.5 kg, the difference between the DFISMC and DSMC is not significant in terms of MTE. However, the RMSTE of the DFISMC controller are smaller than the ones of the DSMC controller, as shown in Figure 9. For example, when tracking the 2.5 s human-gait trajectory, the RMSTEs of the DFISMC and DSMC controllers are 1.22 • and 1.68 • , respectively. Furthermore, the same conclusion about the "chattering" phenomenon is drawn out in this experiment scenario. All numerical values of MTE and RMSTE in this experimental scenario are also shown in Table 5.
From experimental results with multiple scenarios, we can conclude that the DFISMC controller obtains a better tracking performance than the conventional DSMC controller which used the "sign" function of tracking errors. In addition, the implemented disturbance observer and fractional order integral term are able to deal with the finite-amplitude oscillation of sliding mode implementations. As a result, the "chattering" phenomenon is reduced.

Discussion and Conclusions
This paper proposed an advanced SMC control strategy for PAMs in antagonistic configuration. First, the discrete-time SOPDT is chosen to describe the dynamic behaviour of the antagonistic actuator. The chosen model demonstrated a good approximation of nonlinear characteristics of the actuator: the root mean square errors between estimated and measured values are less than 2.5 • . Based on the built-in model, an DFISMC controller, which employed a fractional order integral of tracking error together with a disturbance observer, was proposed for the tracking purpose. The implemented approximation of FOI and DO was able to reduce the "chattering" phenomenon, which often occurs in SMC implementations. The reduction of the "chattering" phenomenon is very important for applications of the PAMs in rehabilitation robot field. Finally, multi-scenario experiments were carried out to compare the tracking performances between the DFISMC and the conventional DSMC.
In comparison with the three-elements model [5], hysteresis model [13][14][15], and mechanism-based model [22,23], the identification procedure of the proposed method is simplified. Besides, this procedure does not need to measure the load's acceleration, which is very sensitive to noise. Experiments show that, in comparison with DSMC, the DFISMC was able to significantly enhance the tracking performance of 20 • amplitude sinusoidal signals with frequency up to 1.0 Hz. In particular, when the actuator drove a load of m = 2.5 kg, the RMSTEs of DFISMC were about two times less than those of conventional DSMC in most of the desired trajectory frequencies. For example, with a frequency of 0.2 Hz, the RMSTEs are 0.93 • and 1.67 • for DFISMC and DSMC, respectively. The proposed controller achieves a performance comparable to the experimental results with similar configuration and desired trajectory in [23,24]. In [23], when tracking a 0.4 Hz frequency and 5 • sinusoidal signal, the residual error amplitude is 0.5 • equivalent to 10%. When tracking a 0.5 Hz frequency and 20 • amplitude sinusoidal signal, the RMSTE of the DFISMC is 1.47 • , equivalent to 7.35% of amplitude. This result was better than the one in [24], in which a sinusoidal signal with 40 • amplitude and frequency 0.25 Hz is used as a desired trajectory. The experiments also show that the proposed controller can track a human-gait pattern with the MTE of less than 6 • . This result is in accordance with the commercial gait training system LOKOMAT [30], in which the MTE is 15 • . It is shown that the built-in model and proposed controller can be applied in robot gait training system. Furthermore, the proposed DFISMC is designed in discrete-time domain, so that it is convenient for implementing in any digital industrial controller, e.g., the NI instrument in this research.
In summary, this paper presents the control of an antagonistic actuator powered by PAMs. The dynamic behaviour of the antagonistic actuator is described by a discrete-time SOPDT model, which requires a simpler identification procedure. The DFISMC controller based on a DSO and the approximated FOI is used to improve the tracking performance. The implementation of DSO and FOI also helps the system reduce the "chattering" phenomenon. The experimental results illustrate the applicability of the proposed model and controller to a robotic gait training system with a human-gait pattern trackable ability. Future work will involve the impedance control of the antagonistic actuator to increase applicability of PAMs in the field of rehabilitation. The impedance of the actuator can be regulated by manipulating the nominal pressure P 0 of two PAMs. To integrate the impedance controller into the system, the relationship between the actuator compliance and nominal pressure P 0 would be considered and modelled in future work.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Fractional Integral Approximation
Fractional-order calculus is a generalization of the integration and differentiation from integer to non-integer order. This appendix introduces only definitions which are widely used in the area of control systems. First, gamma function Γ(z), which is the extension of the factorial for non-integer number z, is introduced as The most important property of the gamma function is Then, the definition of integral of order α ∈ R is presented. In continuous-time domain, the most often used one is the Riemann-Liouville definition: At this time, the FOI is not supported in any programming language. For this reason, its numerical approximation is required to implement the FOI in any real-time control system. In a digital control system with sampling time T s , interval (0, t) can be approximated by k = t Ts sub-intervals. Therefore, Consider that T s is small enough, so that e is constant in each sub-interval. Therefore, Following that, with the weighting factor ω j as follows: Because of the infinite data in (A7), the approximation of FIO cannot be directly implemented in any digital system. In this research, the recursive approximation of FIO in [31] is employed. Denote Ξ e,k−1 as FIO of the tracking error in the last step, and it can be computed as whereẽ j = e j − e j−1 . We apply the short memory principle to (A10) and we can consider two cases: where The FIO is approximated by Equations (A11) and (A12), which can be easily implemented in any digital control system.

Appendix B. Proof of Assumption 1
Assumption 1 is based on the Taylor expansion and can be explained as follows. For a very small constant T s , we have Then, it can be derived from (A14) that Assume that signal p(t) is smooth, and its differential is bounded. Then there exists a constant A such that and (12) holds. Now, ignore the small term O(T 2 s ) and differentiate both sides of (A15). This gives us dp(t) dt − dp(t − T s ) dt ≈ d 2 p(t) dt 2 T s (A18) By using (A15) on the left side of (A18), Again, assume that the second order differential of p(t) is bounded by a constant B, then it leads to which means that (15) holds.