Admittance-Controlled Teleoperation of a Pneumatic Actuator: Implementation and Stability Analysis

Implementation, experimental evaluation and stability analysis of an admittance-controlled teleoperated pneumatic system is presented. A master manipulator navigates a pneumatic slave actuation, which interacts with a human arm as an environment. Considering the external force in the position control loop in the admittance control, enables the slave to handle the external force independent of the master. The proposed control system is evaluated experimentally using the admittance models with different settings. Stability of the control system is analyzed using the concept of Lyapunov exponents. Parametric stability analysis is conducted to show the effect of changing system parameters on stability.


Introduction
Teleoperated robotic systems have been widely employed in various industrial applications [1,2]. A teleoperated robotic system incorporates a master manipulator (hand-controlled device) operated by an operator, a slave manipulator that emulates the master and, a central controller that coordinates the system through a communication channel [2]. If the slave sends the interaction force with the environment (external force) back to the master, the teleoperation system is called bilateral; otherwise, it is called unilateral [3].
Bilateral teleoperation provides a sense of the remote site to the operator. Thus, the operator can deal with the external force by moving the master manipulator [1,2,4]. In unilateral teleoperation, however, the external force is managed by the slave manipulator [5]. Therefore, control of the slave manipulator can be more challenging. On the other hand, applying unilateral teleoperation eliminates the complexity of rendering the external force on the master and makes the overall teleoperation system more stable [6]. It also helps the operator to focus on navigating the slave by eliminating the burden of dealing with the external force on the master side. The focus of this paper is on unilateral teleoperation of pneumatic actuators. Pneumatic actuators are simple and cost-efficient. However, due to nonlinearities inherent in pneumatic systems caused by compressibility of air and friction, accurate positioning of pneumatics is difficult [7]. Figure 1 shows the experimental setup of the system under study. A PHANTOM joystick is used as the master manipulator. It is capable of serial communication using a high speed IEEE1394 communication protocol (FireWire), which is available on many computers. The slave manipulator is a single degree-of-freedom pneumatic actuator, which interacts with the environment (a human subject) through a handle. A computer with Quanser data acquisition board facilitates the connection of the master and the slave through a local network with negligible delay.

Experimental Setup
Actuators 2020, 9,  The pneumatic setup, shown in Figure 1b, is a double rod (FESTO DNC) with 500 mm stroke. The bore and rod diameter are 40 mm and 16 mm respectively. A proportional directional flow control valve (FESTO MPYE) sends the pressurized air to the cylinder to move it according to the movement of the master. The maximum flow capacity of the control valve is 700 L/min at 7 bar (100 psi) absolute supply pressure. The master, force sensor (ARTECH S type load cell), position encoder (BOURNS incremental rotary encoder) and pressure sensors (by Durham Industries) also send data to the control station. The above components are connected to a local network, which is assumed ideal, meaning the network characteristics such as time delay and packet loss are negligible. The sampling frequency is 500 Hz.
Referring to Figure 1b, the equation of motion of the actuator is defined as: where is the mass of the moving parts, is the actuator position, 1 and 2 are the air pressures at of the actuator chambers, and A represents the piston area. is the viscous friction coefficient and is the externally applied load. The dry friction, , is adopted from LuGre friction model [22] without its viscous friction term as viscous friction is already included as a separate term in (1): where 0 is the equivalent spring constant and 1 is the damping coefficient of bristles. The variable , average bristle deflection, can be achieved by solving the following equation [22]: ,

Control Station Master Side
Slave Side Environment The pneumatic setup, shown in Figure 1b, is a double rod (FESTO DNC) with 500 mm stroke. The bore and rod diameter are 40 mm and 16 mm respectively. A proportional directional flow control valve (FESTO MPYE) sends the pressurized air to the cylinder to move it according to the movement of the master. The maximum flow capacity of the control valve is 700 L/min at 7 bar (100 psi) absolute supply pressure. The master, force sensor (ARTECH S type load cell), position encoder (BOURNS incremental rotary encoder) and pressure sensors (by Durham Industries) also send data to the control station. The above components are connected to a local network, which is assumed ideal, meaning the network characteristics such as time delay and packet loss are negligible. The sampling frequency is 500 Hz.
Referring to Figure 1b, the equation of motion of the actuator is defined as: where m is the mass of the moving parts, x s is the actuator position, P 1 and P 2 are the air pressures at of the actuator chambers, and A represents the piston area. b is the viscous friction coefficient and F ext is the externally applied load. The dry friction, F f , is adopted from LuGre friction model [22] without its viscous friction term as viscous friction is already included as a separate term in (1): where σ 0 is the equivalent spring constant and σ 1 is the damping coefficient of bristles. The variable z, average bristle deflection, can be achieved by solving the following equation [22]: where F c and F s are the Coulomb friction and static friction, respectively, and v sv is the Stribeck velocity.
To move the actuator, air pressure should be charged or discharged into chambers. The air pressure is related to the air mass flows, . m 1 and . m 2 , as follows [21,23]: .
In (4) and (5), R is the ideal gas constant, γ is the ratio of specific heats, α is compressibility correction factor and T is air temperature [24]. V 1 and V 2 are the volumes of each actuator chambers, which are functions of the actuator position as expressed in (6) and (7). L is the length of the actuator stroke and V 0 presents cylinder inactive volume. By defining γ = γ(2/(γ + 1)) (γ+1)/(γ−1) /R, the mass flow rate of air through control valve orifice can be expressed as [25]: where w and . ∅ i (i = 1, 2) are orifice area gradient and mass flow per area unit, respectively. x v represents the displacement of the valve spool, C d is discharge coefficient of the valve, P cr is the valve critical pressure ratio, and P a and P s are the absolute downstream and upstream pressures, respectively. The valve orifice area gradient, w, is related to the control signal by the following equation: where u is the control signal, K v is the valve spool position gain and τ is the valve time constant. The values of the parameters of the test rig are shown in Table 1 that are either provided by the manufacturers or obtained experimentally through previous research [26].

Admittance Control
Admittance control simultaneously maintains the desired force and position in the same direction [13,27]. It provides the system with a function which converts the external force, F ext , to a corresponding position called external position, x ext . Figure 2 shows the general block diagram of admittance control where x ext is added to the primary desired position, x m , determined by the master manipulator. In this approach, the position controller tracks a modified desired trajectory, x d , that is the combination of the primary desired trajectory, x m , and the displacement correspondent to the external force, x ext (x d = x m + x ext ). The relationship between the external force, F ext , and the corresponding displacement, x ext , is usually defined by a second-order linear mass-spring-damper model as follows: where M, B and K are admittance parameters corresponding to inertia, damping and stiffness characteristics and s is the Laplace operator. Response of the manipulator to the external force can be tuned by tuning admittance parameters.

Admittance Control
Admittance control simultaneously maintains the desired force and position in the same direction [13,27]. It provides the system with a function which converts the external force, , to a corresponding position called external position, . Figure 2 shows the general block diagram of admittance control where is added to the primary desired position, , determined by the master manipulator. In this approach, the position controller tracks a modified desired trajectory, , that is the combination of the primary desired trajectory, , and the displacement correspondent to the external force, ( = + ). The relationship between the external force, , and the corresponding displacement, , is usually defined by a second-order linear mass-spring-damper model as follows: where , and are admittance parameters corresponding to inertia, damping and stiffness characteristics and is the Laplace operator. Response of the manipulator to the external force can be tuned by tuning admittance parameters.

Position Controller
To position the pneumatic actuator, Sliding Mode Control (SMC) scheme is employed. SMC is a model-based scheme that provides a robust performance despite model uncertainties [21,28,29]. It is the most popular position controller for pneumatic actuators [7,29]. In order to design an SMC algorithm, an integral sliding surface is defined as: where δ is a positive constant and e is the position error, defined as follows: The SMC control law is obtained by summing up an equivalent control component formulated based on the dynamic model, Av eq , with a robustness control component, Av rb : where K v is the valve spool position gain and w is the valve orifice area. Dynamics of the system on the sliding surface can be expressed as [15]: ..  (14) and applying it to (1), the following expression for the equivalent component is achieved: In (15), Px and Fx are derived as: where Considering the slow dynamics of the system, the rate of changes in dry friction, F f , is slow. Thus, (17) can be neglected. Same assumption was made in previous studies [7,28,29]. The role of the robust part of SMC, Av rb , is to provide robustness despite this assumption. Av rb is formulated as follows: where K rb is a positive gain. Because the discontinuouty of sign function in (19) is not ideal for practical implementation, it is approximated by the continuous tanh function as follows: where a is a large positive number.

Simulation Studies
The performance of the proposed unilateral pneumatic system is first evaluated through simulations. Referring to Figure 1, the teleoperated control system receives the displacement of the master, x m , and the external force imparted to the slave from the environment, F ext . The environment is considered spring-dominant which means F ext is proportional to the slave position: where K ext is the stiffness coefficient of the environment. A simplified admittance model is employed whereby the external force is related to the external position by a virtual spring: In (22), K adm is a positive coefficient corresponding to the stiffness term in the admittance model. The step tracking simulation results for the admittance control are shown in Figure 3. The system step input, x m , is 0.1 m as shown in Figure 3a. Considering K ext = 100 N/m, the external force, F ext , is shown in Figure 3b. The external position, x ext , is derived from (22) where is the stiffness coefficient of the environment. A simplified admittance model is employed whereby the external force is related to the external position by a virtual spring: In (22), is a positive coefficient corresponding to the stiffness term in the admittance model. The step tracking simulation results for the admittance control are shown in Figure 3 Figure 4a. The air pressures in cylinder chambers are shown in Figure 4b. The control signal applied to the slave actuator is shown in Figure 4c and does not saturate. The control signal of this system changes within the range of 0-10 V. Hence, with the 5 V control signal, the valve is fully closed. The position tracking error, defined by (12), is shown in Figure 4d. It is seen that there is a reasonable agreement between motions of slave and master actuators.    Figure 4a. The air pressures in cylinder chambers are shown in Figure 4b. The control signal applied to the slave actuator is shown in Figure 4c and does not saturate. The control signal of this system changes within the range of 0-10 V. Hence, with the 5 V control signal, the valve is fully closed. The position tracking error, defined by (12), is shown in Figure 4d. It is seen that there is a reasonable agreement between motions of slave and master actuators.

Stability Analysis
The concept of Lyapunov exponents (LEs) is a well-established numerical approach to study nonlinear systems [9,20]. LEs indicate stability of a nonlinear system by monitoring the long-term behaviour of adjacent state-space trajectories as time evolves. The total number of LEs is equal to the dimension of the state space and the sign of the LEs is used to deduce the stability property of the system. A positive LE corresponds to chaotic or unstable behavior. When all of the exponents are negative, all of the neighboring trajectories converge as time evolves indicating that the system is exponentially stable with a fixed equilibrium point. If there is one zero exponent and the other exponents are negative, we have a stable system with one-dimensional equilibrium [30]. Presenting the required conditions for employing the Lyapunov exponent method is avoided to keep the focus of the paper. More details can be found in [23].

Calculation of Lyapunov Exponents
Consider the following smooth nonlinear system in an n-dimensional state space: where ∈ ℝ , is the state vector and ( ) is continuous and differentiable. To calculate LEs, a "fiducial" trajectory is found by solving (23) using initial condition (0) = 0 . Orthogonal principal axes, , are then defined on the fiducial trajectory. The asymptotic behaviour of the nonlinear system could be determined by monitoring the length of each principal axis over time. To find the length of the i-th ( = 1, . . . , ) principal axis, a linear equation of motion is formed as: where ( ) is the Jacobian matrix:

Stability Analysis
The concept of Lyapunov exponents (LEs) is a well-established numerical approach to study nonlinear systems [9,20]. LEs indicate stability of a nonlinear system by monitoring the long-term behaviour of adjacent state-space trajectories as time evolves. The total number of LEs is equal to the dimension of the state space and the sign of the LEs is used to deduce the stability property of the system. A positive LE corresponds to chaotic or unstable behavior. When all of the exponents are negative, all of the neighboring trajectories converge as time evolves indicating that the system is exponentially stable with a fixed equilibrium point. If there is one zero exponent and the other exponents are negative, we have a stable system with one-dimensional equilibrium [30]. Presenting the required conditions for employing the Lyapunov exponent method is avoided to keep the focus of the paper. More details can be found in [23].

Calculation of Lyapunov Exponents
Consider the following smooth nonlinear system in an n-dimensional state space: . x = f (x) (23) where x ∈ R n , is the state vector and f (x) is continuous and differentiable. To calculate LEs, a "fiducial" trajectory is found by solving (23) using initial condition x(0) = x 0 . Orthogonal principal axes, δx, are then defined on the fiducial trajectory. The asymptotic behaviour of the nonlinear system could be determined by monitoring the length of each principal axis over time. To find the length of the i-th (i = 1, . . . , n) principal axis, a linear equation of motion is formed as: where F(t) is the Jacobian matrix: Each principal axis is a column of ψ t . Since the resulting principal axes tend to fall along the direction of the fastest growing axis, Gram-Schmidt reorthonormalization scheme is applied to make the axes orthogonal to each other and then normalize their lengths in every iteration [31]. The i-th LE is obtained using the following key equation [31]: where T is the total duration of the observance. [23] details a simple example of calculating LEs.

Stability Analysis of the Proposed Control System
The stability analysis of the entire control system described in Section 3 is conducted. The state space model is formed by defining the state space vector as: where x s and v s are displacement and the velocity of the piston, P 1 and P 2 are air pressures in each actuator chamber, x v is the displacement of spool valve as a result of the control signal, z is average bristle deflection and t 0 e dτ is the integral of the position error as defined by (12). Using the state variables, (27), and Equations (1)-(9), (21) and (22), the state space model is constructed as The control signal in terms of state space variables is obtained by substituting (15) and (19) into (13): The associated variables in (29) are defined as: ...
The dry friction in the state space is as follows: The mass flow rate per area unit in state space is defined as: The equilibrium point of (28) is: To study the stability of the control system, the LEs are calculated and summarized in Table 2. Note that theoretically, Lyapunov exponents are defined when t → ∞ . In practice however, we calculate them over limited time duration and the adequate calculation time can be determined by observing the exponents variation over time. As an example, time evolution of λ 6 is shown in Figure 5.  The signs of Lyapunov exponents determine the stability property of the dynamic system. Negative exponents correspond to exponential stability. The first largest exponent being zero indicates stable system with one-dimensional attractor. To understand the physical meaning of the results in Table 2, (38) and (39) are revisited. According to these equations, 1 , 2 and 5 will eventually reach fixed values where the combination of 3 , 4 and 6 will not hold fixed points but satisfy (39). Having three unknowns ( 3 , 4 , 6 ) and one equation, (39) will have a 2dimentional solution; i.e., the equilibrium of (28) is 2-dimentional. This explains the two zero LEs in Table 2. Therefore, the stability of the proposed control system in presence of external force and model uncertainties is proven, as seen in Table 2.

Parametric Stability Analysis
Parametric stability analysis is the stability analysis of a dynamic system as its parameters change. It is beneficial when some parameters of the system are unknown, inaccurate or varying. The concept of Lyapunov exponents can be used for this reason. Because it provides a quantitative measurement of the stability, one can find the stability region of a system by changing the values of the physical parameters of the system or the controller gains [30]. Not only the concept of LEs can assure the stability of a system when the parameters are varying, it can be used to measure the effect of changing parameters on overall stability.
Parametric stability analysis of the teleoperation system under investigation is conducted using the concept of LEs. The stiffness coefficient of the environment denoted by , and controller gain, , were considered for this purpose. Table 3 shows the numerical values of LEs for admittance unilateral teleoperation as varies. It is evident that changing does not change the values of LEs notably. Table 4 shows the numerical values of LEs while the SMC bandwidth gain, δ, varies. The system is stable for the values of δ ≤ 120. However, for δ = 140, the first LE is positive which denote that the system is chaotic.  The signs of Lyapunov exponents determine the stability property of the dynamic system. Negative exponents correspond to exponential stability. The first largest exponent being zero indicates stable system with one-dimensional attractor. To understand the physical meaning of the results in Table 2, (38) and (39) are revisited. According to these equations, x 1 ss , x 2 ss and x 5 ss will eventually reach fixed values where the combination of x 3 ss , x 4 ss and x 6 ss will not hold fixed points but satisfy (39).
Having three unknowns (x 3 ss , x 4 ss , x 6 ss ) and one equation, (39) will have a 2-dimentional solution; i.e., the equilibrium of (28) is 2-dimentional. This explains the two zero LEs in Table 2. Therefore, the stability of the proposed control system in presence of external force and model uncertainties is proven, as seen in Table 2.

Parametric Stability Analysis
Parametric stability analysis is the stability analysis of a dynamic system as its parameters change. It is beneficial when some parameters of the system are unknown, inaccurate or varying. The concept of Lyapunov exponents can be used for this reason. Because it provides a quantitative measurement of the stability, one can find the stability region of a system by changing the values of the physical parameters of the system or the controller gains [30]. Not only the concept of LEs can assure the stability of a system when the parameters are varying, it can be used to measure the effect of changing parameters on overall stability.
Parametric stability analysis of the teleoperation system under investigation is conducted using the concept of LEs. The stiffness coefficient of the environment denoted by K ext , and controller gain, δ, were considered for this purpose. Table 3 shows the numerical values of LEs for admittance unilateral teleoperation as K ext varies. It is evident that changing K ext does not change the values of LEs notably. Table 4 shows the numerical values of LEs while the SMC bandwidth gain, δ, varies. The system is stable for the values of δ ≤ 120. However, for δ = 140, the first LE is positive which denote that the system is chaotic.

Experimental Results
Experiments were performed to evaluate the performance of the teleoperated pneumatic system under real test scenarios. Different admittance parameter settings were tested. Two experimental scenarios are presented to show the performance of the admittance control in conjunction with SMC. In the first experiment, the admittance parameters are set to study soft reaction to an external force. In the second experiment, an external force is applied to the actuator and the goal of the admittance model is a stiff actuator reaction. Experiment 1: In this experiment, the operator moves the master and, at the same time, a human subject located at the slave side applies a force to the pneumatic actuator. The external force passes through the admittance model (10), having parameters set as M = 10 Kg, B = 50 Ns/m and K = 250 N/m. The desired position trajectory originating from the master, x m , is shown in Figure 6a. Figure 6b shows the external force. The modified desired trajectory is shown in Figure 6c. It is evident that the admittance control module effectively adjusts the primary desired trajectory according to the imposed external force. This figure also compares the modified desired trajectory with the actual position of the actuator and illustrates their reasonable agreement. Figure 6d shows the control signal, which is also reasonable and unsaturated. This experiment shows successful application of admittance control with soft stiffness. Experiment 2: In this experiment, the parameters of the admittance model are set to M = 10 Kg, B = 50 Ns/m and K = 10 4 N/m. The goal is to study the behaviour of the system when the stiffness of admittance model is set high. The actuator is subject to an external force with the magnitude of 100 N applied by a human subject. As Figure 7a shows, the primary desired position is fixed during the experiment. The external force is shown in Figure 7b. The admittance model in (10) converts the external force to small displacement as shown in Figure 7c. For external force with 100 N magnitude, the change in primary desired trajectory is about 0.01 m. By comparing Figures 7c and 6c, the effect of changes in the parameters of the admittance control module can be studied. Figure 7c also shows the actual position of slave actuator, x s . Figure 7d shows the control signal corresponding to the position tracking shown in Figure 7c.
The above experiments show the proposed unilateral control system can successfully control the slave actuator desired position in the presence of the external force. The stability of the actuator motion was evident in the experiments. The SMC position controller worked effectively despite the non-idealities such as friction and air compressibility.
The positioning accuracy of the admittance method presented in this paper is now compared with the accuracy of the impedance method, previously developed by the authors on the same system. The comparison is done by calculating the average position error for the same experiments. The results are summarized in Table 5:     The above experiments show the proposed unilateral control system can successfully control the slave actuator desired position in the presence of the external force. The stability of the actuator motion was evident in the experiments. The SMC position controller worked effectively despite the non-idealities such as friction and air compressibility.  According to Table 5, admittance unilateral teleoperation provides higher tracking accuracy. This is expected since admittance unilateral teleoperation utilizes an SMC position controller. Whereas, in impedance unilateral, an SMC force controller is in charge of tracking the position.

Conclusions
In this paper, a unilateral teleoperation system was analyzed, and successfully implemented on a pneumatic actuator. Admittance control, coupled with a sliding mode position controller, was employed to manage actuator position and external force in the slave side. The proposed control system not only inherited the structural advantages of the unilateral teleoperation system with no need for highly skilled operator, but also displayed satisfactory performance in scenarios involving various levels of environmental stiffness and interactions. Both simulation and experimental results showed the effectiveness of the proposed unilateral teleoperation system in position tracking and handling the external force. Meanwhile, stability of the entire control system was evaluated using the concept of Lyapunov exponents. Parametric stability analysis showed the system remains stable for a wide range of the environmental stiffnesses whereas increasing the controller gain leads to system instability. Future work can focus on extension of the controller presented here to multi-degree of freedom pneumatic manipulators, and the stability analysis of the entire system.