Robust Stability of Scaled-Four-Channel Teleoperation with Internet Time-Varying Delays

We describe the application of a generic stability framework for a teleoperation system under time-varying delay conditions, as addressed in a previous work, to a scaled-four-channel (γ-4C) control scheme. Described is how varying delays are dealt with by means of dynamic encapsulation, giving rise to mu-test conditions for robust stability and offering an appealing frequency technique to deal with the stability robustness of the architecture. We discuss ideal transparency problems and we adapt classical solutions so that controllers are proper, without single or double differentiators, and thus avoid the negative effects of noise. The control scheme was fine-tuned and tested for complete stability to zero of the whole state, while seeking a practical solution to the trade-off between stability and transparency in the Internet-based teleoperation. These ideas were tested on an Internet-based application with two Omni devices at remote laboratory locations via simulations and real remote experiments that achieved robust stability, while performing well in terms of position synchronization and force transparency.


Introduction
A teleoperation system consists of master and slave mechanical systems where the master is directly manipulated by a human operator and the slave, operating in a remote environment, is designed to track the master closely. The main issues in the analysis and synthesis of these systems are stability and transparency. Transparency refers to the level of achievement of an ideal situation of the operator acting on the remote environment. In practice, there is a compromise between these two goals mainly due to the presence of time delays generated by the communication channel [1].
Given this compromise, many control schemes for teleoperation have been proposed, broadly classified in terms of two main frameworks, namely, intrinsically stable schemes (passivity-based), and delay-dependent stability schemes [2]. In early works on constant delay, stability was addressed by means of frequency, Laplace or passivity techniques applied to linear time-invariant (LTI) master-slave two-port systems [1][2][3][4][5].
In regard to transparency properties, the most successful control scheme in achieving full transparency under ideal conditions is the four-channel (4C) control scheme [1,2,4,6]. The specialized literature contains very detailed analyses of transparency and control architectures based on four, three and two channels, that, in general assume constant delays (see especially [4] and references therein).
However, in recent years a new line of research on Internet-based teleoperation [6][7][8][9][10][11] exposes the system control loop to the time-varying delay of a packet-switched network-the subject of a tutorial presented by Hokayem and Spong [12]. The second goal is to ensure robust stability in Internet-based teleoperation, that is, under time-varying delay. We address the problem of uncertain delays from a continuous-time perspective (not in discrete-time [26,27]) and for time-varying delays (more realistic than uncertain but constant delays [28]). The induced problem is more difficult but, since useful delay encapsulations have been described in the literature [29], we show how to adapt these delay bounds to achieve our second goal of robustness against continuous-time time-varying delay.
The main differences between our previous work and the present work-in other words, the novelty of the present paper-is that here we explicitly discuss ideal transparency issues and rearrange previous solutions [22,25] so that now the resulting controllers are proper, that is, they do not contain single or double differentiators and so avoid problems of noise. We also validated the ideas with simulations and real-life experiments on an Internet teleoperation setup that coordinated two Omni haptic devices between remote laboratory setups in Spain (University of Vigo, Vigo, Spain) and the UK (University of Manchester, Manchester, UK).
The paper is organized as follows: in Section 2 the model for robust stability in 4C-based teleoperation is introduced and in Section 3 a robust stability condition for this control scheme is defined. Section 4 describes the haptic teleoperation case study. Section 5 presents the analysis and simulation results and Section 6 presents the experimental results. Finally, Section 7 concludes the paper with a discussion.

Model for Robust Stability in 4C Based Teleoperation
Consider the generic teleoperation setup depicted in Figure 1, consisting of an exchange of signals through a communication channel between the master side (left) and the slave side (right). The master side is formed of the plant 1/P m and the controller K m . The plant is modelled as a mechanical system with coordinates (positions) given by the nˆ1 vectors x m , x s . In the Laplace domain, if x m p0q " .
x m p0q " 0, we have: P m psqX m psq :" F m psq, P s psqX s psq :" F s psq (1) Note that we use lowercase and uppercase letters for time-domain and Laplace domain signals, respectively.
Sensors 2016, 16, 593 3 of 20 The second goal is to ensure robust stability in Internet-based teleoperation, that is, under time-varying delay. We address the problem of uncertain delays from a continuous-time perspective (not in discrete-time [26,27]) and for time-varying delays (more realistic than uncertain but constant delays [28]). The induced problem is more difficult but, since useful delay encapsulations have been described in the literature [29], we show how to adapt these delay bounds to achieve our second goal of robustness against continuous-time time-varying delay.
The main differences between our previous work and the present work-in other words, the novelty of the present paper-is that here we explicitly discuss ideal transparency issues and rearrange previous solutions [22,25] so that now the resulting controllers are proper, that is, they do not contain single or double differentiators and so avoid problems of noise. We also validated the ideas with simulations and real-life experiments on an Internet teleoperation setup that coordinated two Omni haptic devices between remote laboratory setups in Spain (University of Vigo, Vigo, Spain) and the UK (University of Manchester, Manchester, UK).
The paper is organized as follows: in Section 2 the model for robust stability in 4C-based teleoperation is introduced and in Section 3 a robust stability condition for this control scheme is defined. Section 4 describes the haptic teleoperation case study. Section 5 presents the analysis and simulation results and Section 6 presents the experimental results. Finally, Section 7 concludes the paper with a discussion.

Model for Robust Stability in 4C Based Teleoperation
Consider the generic teleoperation setup depicted in Figure 1, consisting of an exchange of signals through a communication channel between the master side (left) and the slave side (right). The master side is formed of the plant 1/Pm and the controller Km. The plant is modelled as a mechanical system with coordinates (positions) given by the n × 1 vectors xm, xs. In the Laplace domain, if   where fh is the force applied by the human operator and fmc is the force from the controller. The plant has two outputs, namely, the ny signals that send information to the slave (ym) and the nz signals used by the local controller (zm).
On the slave side, fe is the force from the environment. The controller force fsc enters additively in Master and slave exchange signals through the communication channel. Taking into account our theoretical [10], and experimental [11] discussions on delay bounds using the Internet User Datagram Protocol (UDP), the delay τ(t) can be assumed in a non-restrictive way from teleoperation point of view, to be an unknown time-varying function for which the upper bounds on the magnitude and the variation satisfy The input of the master plant is the master force f m " f h´fmc , where f h is the force applied by the human operator and f mc is the force from the controller. The plant has two outputs, namely, the n y signals that send information to the slave (y m ) and the n z signals used by the local controller (z m ).
On the slave side, f e is the force from the environment. The controller force f sc enters additively in f s " f e`fsc .
Master and slave exchange signals through the communication channel. Taking into account our theoretical [10], and experimental [11] discussions on delay bounds using the Internet User Datagram Protocol (UDP), the delay τ(t) can be assumed in a non-restrictive way from teleoperation point of view, to be an unknown time-varying function for which the upper bounds on the magnitude and the variation satisfy @t ě 0: The time-varying delay can be described by the operator D which is a (nonlinear) system used (in the time domain) in any of the equivalent notations: y jd " D j py j q " D j˝yj " D j y j , pj " m, sq . This means that: y md i ptq " D m y m i ptq :" y mi pt´τ m ptqq, y sd i ptq " D s y s i ptq :" y si pt´τ s ptqq where D m , D s are the corresponding master and slave time-delay operators. To simplify matters we also extend the operator D to the Laplace domain in the equivalent forms Y jd " D j pY j q " D j˝Yj " D j Y j , which simply means that: Y jd psq " D j Y j psq ô Y jd psq " L´y jd ptq¯and y jd ptq " D j y j ptq. We now consider a 4C control scheme in the most general and representative case of delay-dependent schemes. Here, the master and slave exchange, through the communication channel, velocities or positions and forces in both directions. The local controllers are defined in terms of the applied external force, their own velocity/position (master/slave) and the delayed velocity/position and force (slave/master) from the other side. When the system exchanges positions and forces (see Figure 2), the controllers K m , K s can be defined as: The time-varying delay can be described by the operatorD which is a (nonlinear) system used (in the time domain) in any of the equivalent notations: . This means that: where D m, D s are the corresponding master and slave time-delay operators. To simplify matters we also extend the operator D to the Laplace domain in the equivalent forms , which simply means that: We now consider a 4C control scheme in the most general and representative case of delay-dependent schemes. Here, the master and slave exchange, through the communication channel, velocities or positions and forces in both directions. The local controllers are defined in terms of the applied external force, their own velocity/position (master/slave) and the delayed velocity/position and force (slave/master) from the other side. When the system exchanges positions and forces (see Figure 2), the controllers Km, Ks can be defined as: Following the procedure described in a previous work [11], we remodeled this whole system to obtain a single feedback loop containing an LTI block and an uncertain time-varying delay block with the same dimension as the number of the signals exchanged between the master and slave. For this purpose, we state the following Proposition 1 to obtain the model for 4C-based teleoperation: Proposition 1. A teleoperation setup as described in Figures 1 and 2, given by Equations (1)-(5), can be formulated as a compact model for robust stability analysis using a generic approach (as developed in [11]) in the form given by equations Equation (7), Equation (8) and the block diagram in Figure 3: with: Following the procedure described in a previous work [11], we remodeled this whole system to obtain a single feedback loop containing an LTI block and an uncertain time-varying delay block with the same dimension as the number of the signals exchanged between the master and slave. For this purpose, we state the following Proposition 1 to obtain the model for 4C-based teleoperation: Proposition 1. A teleoperation setup as described in Figures 1 and 2 given by Equations (1)-(5), can be formulated as a compact model for robust stability analysis using a generic approach (as developed in [11]) in the form given by equations Equation (7), Equation (8) and the block diagram in Figure 3: with:

8)
Sensors 2016, 16, 593 5 of 20  Proof: The proof of this proposition is developed in Appendix A.  Remark 1. The model described by Equation (7) is directly represented by Figure 3, where four time-varying delays channels inherited from Figure 1 are maintained, bearing in mind that they cannot be combined in a single delay because time-varying delays do not commute with LTI systems.

Ideal Transparency Conditions
Note that, from the teleoperation loop in Figure 3, we can also obtain the zero-delay (D = I  Yd = Y) closed loop dynamics for studying ideal behavior and transparency properties, considering: That is: with the following values: The previous admittance matrix definition can be associated with the impedance Z or hybrid H models: The ideal transparency properties are usually defined on the H matrix, regarding which: Remark 1. The model described by Equation (7) is directly represented by Figure 3, where four time-varying delays channels inherited from Figure 1 are maintained, bearing in mind that they cannot be combined in a single delay because time-varying delays do not commute with LTI systems.

Ideal Transparency Conditions
Note that, from the teleoperation loop in Figure 3, we can also obtain the zero-delay (D = I ñ Y d = Y) closed loop dynamics for studying ideal behavior and transparency properties, considering: Xpsq " pΠpsq`CpsqDpsqq´1 ΦpsqFpsq ": ΛpsqFpsq That is:˜x with the following values: Λ 11 " pP s`Cs q¨p1`C 6 q , Λ 12 "´pC 4`C2 P s q¨p1`C 5 q Λ 21 " pC 1`C3 P m q¨p1`C 6 q , Λ 22 " pP m`Cm q¨p1`C 5 q ∆ Λ " pP m`Cm q¨pP s`Cs q`pC 1`C3 P m q pC 4`C2 P s q The previous admittance matrix definition can be associated with the impedance Z or hybrid H models: The ideal transparency properties are usually defined on the H matrix, regarding which: Comparing Equations (13) and (14) with Equation (9), we obtain: Therefore, these properties can also be analyzed through Λ and Z matrix in the following form: Transparency and Condition set "

Robust Stability for 4C Based Teleoperation
Below we apply our previously developed approach to stability under time-varying delay in teleoperation [11,22] to the 4C control scheme. The time-varying delay is an uncertainty that may give rise to instability. To ensure robust stability, of necessity we have to impose some limits on maximum delay and maximum delay variation. Our treatment of time-varying delay is based on an encapsulation described in the literature [29] and adapted by us to a 4C teleoperation setup. This delay encapsulation [29]-based on the concept of the integral quadratic constraint (IQC), which is a powerful framework for robust stability-is combined with input-output stability criteria and µ-analysis and synthesis techniques in order to achieve a final robust stability condition (Theorem 1).

Proposition 2.
A teleoperation setup, as stated in Proposition 1 and depicted in Figure 3, can be modeled as a negative single-feedback loop containing an LTI block G(s) and an uncertain time-varying delay block D , shown in Figure 4a, with: G "¨0 2ˆ2 pP m`Cm q´1 C 4 pP m`Cm q´1 C 2 P m pP m`Cm q´1 C 4 P m pP m`Cm q´1 C 2 pP s`Cs q´1 C 1´p P s`Cs q´1 C 3 P s pP s`Cs q´1 C 1´Ps pP s`Cs q´1 Theorem 1. Consider a 4C-based teleoperation system, modeled by Equations (7) and (8), directly represented by Figure 3 and transformed into Figure 4a by Proposition 2. Given G(s) as defined in Equation (15), let φ(s) [29]: 2kc where c is any positive real number and delay bounds h max and where d is as defined in Equations (2) and (3). If necessary, [29] provides extensions of the delay encapsulations by another multipliers ϕ in Equation (19) that deal with great variability. A sufficient condition for stability of the delayed system as described in Figure 4a, with a complex diagonal structured uncertainty, is given by: ρ rHpjωqs ă 1 (20) where: and where µ p¨q is the structured singular value (SSV) with respect to a repeated complex scalar uncertainty, which is, in turn, equal to ρ p¨q, the spectral radius of a matrix (the maximum of the norms of the eigenvalues).
Proof. This result is obtained in two ways. First, in applying the loop transformation theorem (the loop depicted in Figure 4a is stable if the loop depicted in Figure 4b is stable) and the small gain theorem to Figure 4b, whereby Hpsq L 2¨ ∆ L 2 ă 1, since ∆ L 2 ď 1 by [29], a sufficient condition is Hpsq L 2 ă 1. Second, using the µ-techniques for robust stability [30] we can interpret Figure 4b as a nominal system H(s) connected to a dynamic or complex uncertainty ∆ that is actually diagonal (more details are given in [11]).

Remark 2.
As a prerequisite for robust stability, we must first satisfy nominal stability. Hence, H(s) must be Hurwitz-stable, that is, all the λ i eigenvalues must verify real pλ i q ă 0 @i, i " 1 . . . n . and where    is the structured singular value (SSV) with respect to a repeated complex scalar uncertainty, which is, in turn, equal to   , the spectral radius of a matrix (the maximum of the norms of the eigenvalues).
Proof. This result is obtained in two ways. First, in applying the loop transformation theorem (the loop depicted in Figure 4a is stable if the loop depicted in Figure 4b is stable) and the small gain theorem to Figure 4b, whereby [29], a sufficient condition is 2 ( ) 1 L Hs  . Second, using the μ-techniques for robust stability [30] we can interpret Figure 4b as a nominal system H(s) connected to a dynamic or complex uncertainty Δ that is actually diagonal (more details are given in [11]). 

Remark 2.
As a prerequisite for robust stability, we must first satisfy nominal stability. Hence, H(s) must be Hurwitz-stable, that is, all the λi eigenvalues must verify This frequency approach based on IQCs combined with μ-analysis and synthesis techniques is a powerful framework for robust stability. Therefore, using the stability condition of Equation (20) based on the system in Equation (21), for any delay satisfying Equations (2) and (3) and for some designed closed-loop (non-delayed) dynamic in G, we can ensure robust stability under time-varying delays.

Haptic Teleoperation Case Study
The previous 4C architectural framework was applied to remote teleoperation between master and slave haptic devices located in laboratories in the University of Vigo and the University of Manchester.
Below we describe the details of the setup and of master and slave device identification. We then describe a scaled control scheme γ-4C in which the stability factor γ lets us obtain a practical and internally stable solution to Internet-based teleoperation, where all the controllers are proper (without differentiators), as is usual in classical solutions. Further sections will describe simulations and experimental results.
The master and slave sides consist of two Phantom Omni haptic devices manufactured by SensAble Technologies Inc. (Wilmington, MA, USA). In terms of degree of freedom (DoF), these devices have 6-DoF position sensing and 3-DoF force actuation. To focus on the robust stability issue, we limited movements to 1-DoF-rather than multi-DoF-using the first rotational coordinate (the base rotation) as the free DoF. In this way, the two Omnis rotated around their vertical axes, while the shoulder, elbow and stylus were blocked.
Thus, the master and slave positions were (xm, xs), the sensed Omni base angles in radians (rad). The actuations are (fm, fs), which are the torque inputs to the motors for vertical axis rotation in machine units (we did not have access to physical units). This frequency approach based on IQCs combined with µ-analysis and synthesis techniques is a powerful framework for robust stability. Therefore, using the stability condition of Equation (20) based on the system in Equation (21), for any delay satisfying Equations (2) and (3) and for some designed closed-loop (non-delayed) dynamic in G, we can ensure robust stability under time-varying delays.

Haptic Teleoperation Case Study
The previous 4C architectural framework was applied to remote teleoperation between master and slave haptic devices located in laboratories in the University of Vigo and the University of Manchester.
Below we describe the details of the setup and of master and slave device identification. We then describe a scaled control scheme γ-4C in which the stability factor γ lets us obtain a practical and internally stable solution to Internet-based teleoperation, where all the controllers are proper (without differentiators), as is usual in classical solutions. Further sections will describe simulations and experimental results.
The master and slave sides consist of two Phantom Omni haptic devices manufactured by SensAble Technologies Inc. (Wilmington, MA, USA). In terms of degree of freedom (DoF), these devices have 6-DoF position sensing and 3-DoF force actuation. To focus on the robust stability issue, we limited movements to 1-DoF-rather than multi-DoF-using the first rotational coordinate (the base rotation) as the free DoF. In this way, the two Omnis rotated around their vertical axes, while the shoulder, elbow and stylus were blocked.
Thus, the master and slave positions were (x m , x s ), the sensed Omni base angles in radians (rad). The actuations are (f m , f s ), which are the torque inputs to the motors for vertical axis rotation in machine units (we did not have access to physical units).
We concentrated on 1-DoF movements around a fixed position (the "zero" of the Omnis). Although the Omnis are robotic nonlinear systems, given the linearization principle, it was expected that the local dynamics from forces to positions could be captured as is usual by models in the form G i (s) = 1/(s(Ms + B)) [31]. Note that here, for simplicity sake, we use the terms "forces" and "positions" to denote torques and angles. Note also that, since the forces are in machine units, the inertia-mass M and friction B coefficients were given in machine or virtual units.
Before identifying the models G i (s), our experiments showed that there was a deadzone or static friction at the force inputs. This meant that the effective force f ef was smaller than the applied or written force f i . This effect can be modelled by a deadzone function f ef = DZ(f i ; D pos , D neg ) in the form: Although this static friction slightly damaged linearity, producing small steady-state errors, this was easily mitigated by pre-compensation at input by an inverse function in the form f i = DZ(f ;´D pos ,´D neg ) with negative values´D pos ,´D neg . If the deadzone (asymmetric) widths D pos , D neg were estimated well, the composition of these two functions would give rise to the identity, f ef = f, thus linearizing the system. After linearizing the devices by deadzone compensation, the linear models G i (s) = 1/(s(Ms + B)), were obtained, by a standard procedure [24] applied in closed -loop (proportional controller) square and triangular references, recording force inputs and position outputs and performing standard least-squares identifications of M, B. The results (in machine units, m.u.) were: Therefore, for the teleoperation setup described in Figure 1 formed of two identical OMNI Phantoms as the master and slave systems, for each DoF, the linearized dynamic models were: Seeking a practical solution to Internet-based teleoperation, we fine-tuned a scaled 4C control scheme γ-4C. We also imposed internal stability, that is, the closed-loop system had to be Hurwitz-stable. The operator and environment were modelled as external forces.
These requirements were met in the scheme depicted in Figures 1 and 2 and given by Equations (1)-(5), with master and slave systems as in Equation (24) and with the controllers in Equation (5) defined in the following form: For C m " C s " K ą 0 and γ ą 1 : The Λ and Z matrix in this γ-4C scheme with controllers given by Equation (25) take the values:

26)
Note that the control parameters in Equation (5) with the values selected in Equation (25) were not eight independent controllers. Due to the transparency restriction requirements, there were ultimately only a minimum of two DoFs in the design: the gain K and the coefficient γ (K' can be zero). The gain K was the only parameter that directly influenced the closed-loop poles (∆ Λ = 0) in the non-delayed case, as seen in Equation (26). The value k > 0 was necessary to ensure non-delayed system stability. We could thus tune k > 0 from pole-placement principles, bearing in mind that the feedback loop would be affected by a (time-varying) delay. As for the coefficient γ, this ideally should be set to γ = 1 to ensure perfect transparency, but to avoid singularity in the feedback loop impedance (admittance) matrix, it was set to γ > 1. For greater transparency, γ should be as close as 1 as possible, but has to be greater than 1 to ensure robust stability and loop well-posedness.

Remark 3.
If γ = 1 + ε, with ε Ñ 0, the teleoperation system with the controllers defined by Equation (25) verify the transparency-optimized conditions for zero delay Equation (16). But in the ideal case ε = 0, the impedance matrix Equation (12) loses rank and thus the admittance model Equation (10) is thus not well-defined (ill-posed).

Remark 4.
To avoid the previous issue, the γ factor in Equation (25) can be selected to verify well-posedness and robust stability conditions under Internet time-varying delay, while maintaining the transparency properties of the system as much as possible. This control scheme therefore gives internal stability and a practical solution to the trade-off between performance and stability under time-varying delay.
Remark 5. The proposed control scheme in Figures 1 and 2 as given by Equation (5) allows proper controllers in Equation (25) to be selected based on ideally optimized transparency, that is, they do not contain single or double differentiators like classical 4C solutions and so avoid negative effects with noise. See the explanation in Appendix B. (25) also provide robust stability and performance against parametric uncertainty. The closed-loop poles are now given by (M m s 2 + B m s + K) (M s s 2 + B s s + K) = 0 and, consequently, for zero delay, the system will be stable iff K, M m,s = M˘∆ M , B m,s = B˘∆ B >0.

Remark 6. The controllers defined in Equation
These remarks are further analyzed in the following sections.

Analysis and Simulation Results
We first describe a simulation case for use with Simulink and Matlab, considering linear dynamic models for the master and slave, the rounded identified values given by M = 0.001, B = 0.02 and the controllers as stated in Equation (25), with K = 0.1 and K' = 0.
To justify choosing this control scheme based on performance criteria, using this simulation case and applying zero delay, we compared our scheme with some frequently used two-and three-channel (2C and 3C) control schemes. These schemes can be obtained by canceling certain controller parameters in Equation (25), as shown in Table 1. Table 1. Control Schemes (P m = P s = P = 0.001s 2 + 0.02s, C m = C s = K = 0.1, C 5 = C 6 = K' = 0).
ote that the 4C scheme results were obtained for γ = 1.07 « 1 to avoid singularities. The selected value in K provided an adequate time response by fixing equal and real closed-loop poles, that is, the linear model of the system had critical damping and so did not overshoot.
Performance can be studied by simulations that apply a simulated sine reference force (f h )-like the force imposed by a human operator-to move the master. The slave will follow the master in free motion in the interval t = (0, 10) and (18, 30) s and with an environment step force f e simulated from t = 10 s and t = 18 s Note that the operator and environment are not assumed to be in contact or to be passive.
The best results for tracking properties were obtained for the schemes which approximated the ideal transparency conditions (γ-4C) or, at least, the schemes (PE, 3Ci, 3Cii) in which conditions were achieved in the steady state (s Ñ 0). The simulation results, for brevity, are shown only for these schemes (Figures 5-8). Only the γ-4C scheme provided internal stability to the system. Because the other schemes have a closed-loop pole at the origin, without external forces, some states of the system tended to the same value (in this case x m = x s ) but not to zero, for example, PE in Figure 9a. However, the γ-4C scheme in Figure 9b, without external forces, drove the master and slave to their zero states. Note that the 4C scheme results were obtained for γ = 1.07 ≈ 1 to avoid singularities. The selected value in K provided an adequate time response by fixing equal and real closed-loop poles, that is, the linear model of the system had critical damping and so did not overshoot.
Performance can be studied by simulations that apply a simulated sine reference force (fh)-like the force imposed by a human operator-to move the master. The slave will follow the master in free motion in the interval t = (0, 10) and (18,30) s and with an environment step force fe simulated from t = 10 s and t = 18 s Note that the operator and environment are not assumed to be in contact or to be passive.  Note that the 4C scheme results were obtained for γ = 1.07 ≈ 1 to avoid singularities. The selected value in K provided an adequate time response by fixing equal and real closed-loop poles, that is, the linear model of the system had critical damping and so did not overshoot.
Performance can be studied by simulations that apply a simulated sine reference force (fh)-like the force imposed by a human operator-to move the master. The slave will follow the master in free motion in the interval t = (0, 10) and (18,30) s and with an environment step force fe simulated from t = 10 s and t = 18 s Note that the operator and environment are not assumed to be in contact or to be passive.  Note that the 4C scheme results were obtained for γ = 1.07 ≈ 1 to avoid singularities. The selected value in K provided an adequate time response by fixing equal and real closed-loop poles, that is, the linear model of the system had critical damping and so did not overshoot.
Performance can be studied by simulations that apply a simulated sine reference force (fh)-like the force imposed by a human operator-to move the master. The slave will follow the master in free motion in the interval t = (0, 10) and (18,30) s and with an environment step force fe simulated from t = 10 s and t = 18 s Note that the operator and environment are not assumed to be in contact or to be passive.   The best results for tracking properties were obtained for the schemes which approximated the ideal transparency conditions (γ-4C) or, at least, the schemes (PE, 3Ci, 3Cii) in which conditions were achieved in the steady state (s → 0). The simulation results, for brevity, are shown only for these schemes (Figures 5-8). Only the γ-4C scheme provided internal stability to the system. Because the other schemes have a closed-loop pole at the origin, without external forces, some states of the system tended to the same value (in this case xm = xs) but not to zero, for example, PE in Figure 9a. However, the γ-4C scheme in Figure 9b, without external forces, drove the master and slave to their zero states. Furthermore, it may happen that, for some input forces as used in remote teleoperation, the master and slave positions are not bounded for systems without internal stability, although position tracking is accomplished-for instance, one of the schemes in Figure 10 and the solution given for the 4C scheme in Figure 11.   The best results for tracking properties were obtained for the schemes which approximated the ideal transparency conditions (γ-4C) or, at least, the schemes (PE, 3Ci, 3Cii) in which conditions were achieved in the steady state (s → 0). The simulation results, for brevity, are shown only for these schemes (Figures 5-8). Only the γ-4C scheme provided internal stability to the system. Because the other schemes have a closed-loop pole at the origin, without external forces, some states of the system tended to the same value (in this case xm = xs) but not to zero, for example, PE in Figure 9a. However, the γ-4C scheme in Figure 9b, without external forces, drove the master and slave to their zero states. Furthermore, it may happen that, for some input forces as used in remote teleoperation, the master and slave positions are not bounded for systems without internal stability, although position tracking is accomplished-for instance, one of the schemes in Figure 10 and the solution given for the 4C scheme in Figure 11.  Furthermore, it may happen that, for some input forces as used in remote teleoperation, the master and slave positions are not bounded for systems without internal stability, although position tracking is accomplished-for instance, one of the schemes in Figure 10 and the solution given for the 4C scheme in Figure 11. The best results for tracking properties were obtained for the schemes which approximated the ideal transparency conditions (γ-4C) or, at least, the schemes (PE, 3Ci, 3Cii) in which conditions were achieved in the steady state (s → 0). The simulation results, for brevity, are shown only for these schemes (Figures 5-8). Only the γ-4C scheme provided internal stability to the system. Because the other schemes have a closed-loop pole at the origin, without external forces, some states of the system tended to the same value (in this case xm = xs) but not to zero, for example, PE in Figure 9a. However, the γ-4C scheme in Figure 9b, without external forces, drove the master and slave to their zero states. Furthermore, it may happen that, for some input forces as used in remote teleoperation, the master and slave positions are not bounded for systems without internal stability, although position tracking is accomplished-for instance, one of the schemes in Figure 10 and the solution given for the 4C scheme in Figure 11.   Finally, referring back to Remark 6, we used the γ-4C controllers defined in Equation (25) with master plant as in Equation (24), except that now the slave model was defined as The results can be seen in Figure 12. Due to the fact that the models were different, the master and slave forces-obtained from the external forces and local proportional controllers-were also different in order to accomplish position tracking. Once the performance of the different schemes was analyzed, we studied stability when the time-varying delay was considered. The idea was to maintain as much as possible the ideal tracking properties obtained by the standard controllers in the 4C control scheme by increasing, by means of the γ factor, the stability margin in practical conditions, that is, when there was time-varying delay in signals transmission.
Once the non-delayed part of the master and slave dynamics were fixed (including the gain K), only γ remained to be chosen. All the delay uncertainty was characterized by two measures: the maximum delay hmax and the maximum delay rate d.
The stability condition in Equation (20)  There was an easy way to deal with Equation (20) numerically: the delay uncertainty parameters (hmax, d) could be discretized in the ranges of interest, and for each discretized value, the minimum γ could be computed. These computations can be represented in the form of 3D surfaces Finally, referring back to Remark 6, we used the γ-4C controllers defined in Equation (25) with master plant as in Equation (24), except that now the slave model was defined as P s " M 3 s 2`2 Bs. The results can be seen in Figure 12. Due to the fact that the models were different, the master and slave forces-obtained from the external forces and local proportional controllers-were also different in order to accomplish position tracking. Finally, referring back to Remark 6, we used the γ-4C controllers defined in Equation (25) with master plant as in Equation (24), except that now the slave model was defined as The results can be seen in Figure 12. Due to the fact that the models were different, the master and slave forces-obtained from the external forces and local proportional controllers-were also different in order to accomplish position tracking. Once the performance of the different schemes was analyzed, we studied stability when the time-varying delay was considered. The idea was to maintain as much as possible the ideal tracking properties obtained by the standard controllers in the 4C control scheme by increasing, by means of the γ factor, the stability margin in practical conditions, that is, when there was time-varying delay in signals transmission.
Once the non-delayed part of the master and slave dynamics were fixed (including the gain K), only γ remained to be chosen. All the delay uncertainty was characterized by two measures: the maximum delay hmax and the maximum delay rate d.
The stability condition in Equation (20)  There was an easy way to deal with Equation (20) numerically: the delay uncertainty parameters (hmax, d) could be discretized in the ranges of interest, and for each discretized value, the minimum γ could be computed. These computations can be represented in the form of 3D surfaces Once the performance of the different schemes was analyzed, we studied stability when the time-varying delay was considered. The idea was to maintain as much as possible the ideal tracking properties obtained by the standard controllers in the 4C control scheme by increasing, by means of the γ factor, the stability margin in practical conditions, that is, when there was time-varying delay in signals transmission.
Once the non-delayed part of the master and slave dynamics were fixed (including the gain K), only γ remained to be chosen. All the delay uncertainty was characterized by two measures: the maximum delay h max and the maximum delay rate d.
The stability condition in Equation (20) is based on the Hpsq ": φpsq¨pI`Gpsqq´1 Gpsq system in Equation (21). The closed-loop (non-delayed) system G(s) only depended on γ, and the "multiplier" φ in Equation (19) depended on h max , d. Thus the condition in Equation (20) took the form Γ pγ, h max , dq ă 1. It is important to emphasize that Equation (20) involves the spectral radius of the frequency response H(jω) of some H(s) depending on (γ, h max , d) in a complex way. Thus, Equation (20) had to be treated numerically (not analytically).
There was an easy way to deal with Equation (20) numerically: the delay uncertainty parameters (h max , d) could be discretized in the ranges of interest, and for each discretized value, the minimum γ could be computed. These computations can be represented in the form of 3D surfaces and the designer can use the numerical results to decide, depending on the expected delay uncertainty, the adequate tuning of γ to be included in the controllers.
The stability analysis could be assessed applying Theorem 1 and using Matlab. In all cases it had to be proved, firstly, that the system verified the nominal stability condition given in Equation (22). Here we show, briefly as Figure 13, the effect on stability (from the results on SSV) of different bounds on the delay, fixing the stability factor at γ = 12, which was the value selected for the experiments described in the next section. Note that the delay variation was the most critical aspect of stability. The results depicted in Figure 14 show the minimum value of the stability factor γ (z axis) needed to ensure the stability of the system given the bounds on the delay (h max , d).
Sensors 2016, 16,593 13 of 20 and the designer can use the numerical results to decide, depending on the expected delay uncertainty, the adequate tuning of γ to be included in the controllers. The stability analysis could be assessed applying Theorem 1 and using Matlab. In all cases it had to be proved, firstly, that the system verified the nominal stability condition given in Equation (22). Here we show, briefly as Figure 13, the effect on stability (from the results on SSV) of different bounds on the delay, fixing the stability factor at γ = 12, which was the value selected for the experiments described in the next section. Note that the delay variation was the most critical aspect of stability. The results depicted in Figure 14 show the minimum value of the stability factor γ (z axis) needed to ensure the stability of the system given the bounds on the delay (hmax, d).

Experimental Results
In order to test the above simulation results, real experiments were conducted from the haptic teleoperation case study as described above. The Omni devices were connected to the computer through a FireWire connection. The computers contained the local controllers implemented under Matlab/Simulink and Quarc real-time control software as a real-time machine for the experiments. The Quarc software managed the interface with the haptic device using Quarc Phantom Omni software and also collected and saved all exchanged data. The sampling period in the experiments and the designer can use the numerical results to decide, depending on the expected delay uncertainty, the adequate tuning of γ to be included in the controllers. The stability analysis could be assessed applying Theorem 1 and using Matlab. In all cases it had to be proved, firstly, that the system verified the nominal stability condition given in Equation (22). Here we show, briefly as Figure 13, the effect on stability (from the results on SSV) of different bounds on the delay, fixing the stability factor at γ = 12, which was the value selected for the experiments described in the next section. Note that the delay variation was the most critical aspect of stability. The results depicted in Figure 14 show the minimum value of the stability factor γ (z axis) needed to ensure the stability of the system given the bounds on the delay (hmax, d).

Experimental Results
In order to test the above simulation results, real experiments were conducted from the haptic teleoperation case study as described above. The Omni devices were connected to the computer through a FireWire connection. The computers contained the local controllers implemented under Matlab/Simulink and Quarc real-time control software as a real-time machine for the experiments. The Quarc software managed the interface with the haptic device using Quarc Phantom Omni software and also collected and saved all exchanged data. The sampling period in the experiments

Experimental Results
In order to test the above simulation results, real experiments were conducted from the haptic teleoperation case study as described above. The Omni devices were connected to the computer through a FireWire connection. The computers contained the local controllers implemented under Matlab/Simulink and Quarc real-time control software as a real-time machine for the experiments. The Quarc software managed the interface with the haptic device using Quarc Phantom Omni software and also collected and saved all exchanged data. The sampling period in the experiments was set to 1 millisecond. The two devices were remotely coordinated under programmed force profiles like the human and environment external forces, and the experimental data validated the stability and performance results in the theoretical part.
We describe results for three experiments. The first experiment ( Figure 16) reproduced real local teleoperation of the Omni devices at the University of Vigo, artificially increasing a time-varying delay between master and slave. The simulation case depicted in Figure 15 is compared with this first experiment as depicted in Figure 16.
Sensors 2016, 16,593 14 of 20 was set to 1 millisecond. The two devices were remotely coordinated under programmed force profiles like the human and environment external forces, and the experimental data validated the stability and performance results in the theoretical part. We describe results for three experiments. The first experiment ( Figure 16) reproduced real local teleoperation of the Omni devices at the University of Vigo, artificially increasing a time-varying delay between master and slave. The simulation case depicted in Figure 15 is compared with this first experiment as depicted in Figure 16. The results were very similar to the simulation results, although the difference between the linear identified models used in the simulation case and the OMNI nonlinear dynamic with a pre-compensated deadzone Equation (22) can be appreciated.   was set to 1 millisecond. The two devices were remotely coordinated under programmed force profiles like the human and environment external forces, and the experimental data validated the stability and performance results in the theoretical part. We describe results for three experiments. The first experiment ( Figure 16) reproduced real local teleoperation of the Omni devices at the University of Vigo, artificially increasing a time-varying delay between master and slave. The simulation case depicted in Figure 15 is compared with this first experiment as depicted in Figure 16. The results were very similar to the simulation results, although the difference between the linear identified models used in the simulation case and the OMNI nonlinear dynamic with a pre-compensated deadzone Equation (22) can be appreciated.  The results were very similar to the simulation results, although the difference between the linear identified models used in the simulation case and the OMNI nonlinear dynamic with a pre-compensated deadzone Equation (22) can be appreciated.
The second experiment ( Figure 17) reproduced real remote teleoperation between the University of Vigo and the University of Manchester using linear models of master and slave like the corresponding simulation model described in Section 5. Here we could analyze the impact of time delays separately from model uncertainty and nonlinearities issues. The final experiment ( Figure 18) tested real remote teleoperation between the master haptic device located at the University of Vigo and the slave haptic device located at the University of Manchester. Note that, for implementation over the haptic devices, the deadzone was pre-compensated at the input Equation (22) by an inverse function in the form explained in Section 4. In general, this function caused the system to be more oscillating than the linear model and to have a small steady-state error. The results were very similar to the simulation results, although the difference between the linear identified models used in the simulation case and the OMNI nonlinear dynamic with a pre-compensated deadzone Equation (22) can be appreciated.   The second experiment ( Figure 17) reproduced real remote teleoperation between the University of Vigo and the University of Manchester using linear models of master and slave like the corresponding simulation model described in Section 5. Here we could analyze the impact of time delays separately from model uncertainty and nonlinearities issues. The final experiment ( Figure 18) tested real remote teleoperation between the master haptic device located at the University of Vigo and the slave haptic device located at the University of Manchester. Note that, for implementation over the haptic devices, the deadzone was pre-compensated at the input Equation (22) by an inverse function in the form explained in Section 4. In general, this function caused the system to be more oscillating than the linear model and to have a small steady-state error. The data (forces and positions) were exchanged between the master and the slave computers by a transmission of UDP packets. Each packet was composed of time stamps and the transmitted value with a fixed size of 68 bytes. Figure 19 shows some experimental data for connections via Internet between Vigo and Manchester. The delay was measured in terms of the round trip time (RTT) communication. The maximum delay was near hmax = 0.028 s, with very low variation tested with a transmission rate of 1000 packets per second (Figure 19a). The test over the communication channel also revealed that variations in time delay never exceeds d = 0.45, as can be seen in Figure 19b. It should be noted that the communication network used in the experiment belongs to the GÉANT pan-European research and education network, which has very good quality linkage between facilities.  The data (forces and positions) were exchanged between the master and the slave computers by a transmission of UDP packets. Each packet was composed of time stamps and the transmitted value with a fixed size of 68 bytes. Figure 19 shows some experimental data for connections via Internet between Vigo and Manchester. The delay was measured in terms of the round trip time (RTT) communication. The maximum delay was near h max = 0.028 s, with very low variation tested with a transmission rate of 1000 packets per second (Figure 19a). The test over the communication channel also revealed that variations in time delay never exceeds d = 0.45, as can be seen in Figure 19b. It should be noted that the communication network used in the experiment belongs to the GÉANT pan-European research and education network, which has very good quality linkage between facilities. The second experiment ( Figure 17) reproduced real remote teleoperation between the University of Vigo and the University of Manchester using linear models of master and slave like the corresponding simulation model described in Section 5. Here we could analyze the impact of time delays separately from model uncertainty and nonlinearities issues. The final experiment ( Figure 18) tested real remote teleoperation between the master haptic device located at the University of Vigo and the slave haptic device located at the University of Manchester. Note that, for implementation over the haptic devices, the deadzone was pre-compensated at the input Equation (22) by an inverse function in the form explained in Section 4. In general, this function caused the system to be more oscillating than the linear model and to have a small steady-state error. The data (forces and positions) were exchanged between the master and the slave computers by a transmission of UDP packets. Each packet was composed of time stamps and the transmitted value with a fixed size of 68 bytes. Figure 19 shows some experimental data for connections via Internet between Vigo and Manchester. The delay was measured in terms of the round trip time (RTT) communication. The maximum delay was near hmax = 0.028 s, with very low variation tested with a transmission rate of 1000 packets per second (Figure 19a). The test over the communication channel also revealed that variations in time delay never exceeds d = 0.45, as can be seen in Figure 19b. It should be noted that the communication network used in the experiment belongs to the GÉANT pan-European research and education network, which has very good quality linkage between facilities.

Conclusions
In this paper we have extended our generic framework to stability under time-varying delay [11] to scaled 4C or γ-4C teleoperation. The application of the theoretical framework to the four channel case, with time-varying delays, gave rise to an encapsulated delay uncertainty, size 4ˆ4 and diagonal. Assuming symmetric delays, the related structured singular value or µ-value was equal to the spectral radius of the linear subsystem. The approach represents an appealing frequency technique to deal with the stability robustness of the architecture.
The main novelty of the present research compared to previous related research is that we have explicitly discussed the problems associated with ideal transparency and that we have rearrange previous 4C solutions [22,25] so that the resulting controllers are proper, that is, they do not contain single or double differentiators and so avoid problems of noise. We have also validated this research through new simulations and real-life experiments using an Internet teleoperation setup that coordinates two Omni haptic devices between remote laboratory setups in Spain (University of Vigo) and the UK (University of Manchester).
The experiments achieved remote device coordination under programmed force profiles. This required resolving a large number of practical issues related to hardware, software and communications. Some important practical problems were related to linearized model identification and to deadzone or static friction in the Omni devices, which had to be pre-compensated. Finally, the architecture was fine-tuned and successfully tested, obtaining complete stability to zero of the whole state and robustness against continuous-time time-varying delay for Internet-based teleoperation. The real and simulated experiments give similar results, thus validating the models. Stable remote teleoperation also performed adequately in terms of position synchronization and force transparency.
This research has several possible extensions to future work. First, time-varying delay uncertainty was treated here by means of encapsulation [29], but exploring other novel approaches to uncertain time delays could improve the robustness analysis, for example, Lyapunov-Krasovskii functionals [32]. The generic framework based on IQCs [33] is also a potentially powerful tool in addressing teleoperation robustness. Second, one of the main difficulties encountered was the problem of deadzones as, although the coarse values of the (asymmetrical) deadzone widths can easily be estimated, the related static friction does not repeat well from one experiment to another. Since deadzone is a static monotone bounded uncertainty, it could potentially be addressed by means of the IQC approach that uses Zames-Falb multipliers to deal with robust stability [34].

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

Appendix A
In the proposed scheme the master and slave are modeled as mechanical systems given by the following equation in Laplace domain Equation (A1): P m psqX m psq :" F m psq, P s psqX s psq :" F s psq (A1) The output vector of the plant y m = (x m f m ) T that sends information to the slave in the 4C scheme is a n yˆ1 = 2ˆ1 vector because the system exchanges both positions and forces. It can be expressed in Laplace domain as a function of the scalar X m : Y m psq " pD 1m`D2m s`D 3m P m psqq X m psq ": D m psqX m psq ñ D m "˜1 P m¸, D 1m "˜1 0¸, D 2m "˜0 0¸, D 3m "˜0 1¸(

A2)
The other output of the plant z m , used by the local controller, is a n z = 1 signal which is the position. Hence: Z m psq " pE 1m`E2m sq X m psq ": E m psqX m psq ñ E m "´1 0¯(A3) A generic controller in the Laplace domain that uses the previous notation (see Figure 1) can be described by: F mc psq " C 1m psqZ m psq`C 2m psqY sd psq`C 3m psqF h psq (A4) where the values of C im can be obtained by comparing the master controller in 4C as expressed in (5) with Equation (A4).
C 1m " C m , C 2m "´C 4 C 2¯, C 3m "´C 6 (A5) Combining all the previous equations, and omitting, for simplicity sake, dependence on s, we obtain: P m X m " F m " F h´Fmc " F h´C1m E m X m´C2m Y sd´C3m F h , pP m`C1m E m q X m "´C 2m Y sd`p I´C 3m q F h (A6) On the slave side, the environment force is f e , whereas the controller force is f sc and thus f s " f e`fsc . Hence: pP s´C1s E s q X s " C 2s Y md`p I`C 3s q F e C 1s "´C s , C 2s "´C 1 C 3¯, C 3s " C 5 , E s "´1 0D s "˜1 P s¸, D 1s "˜1 0¸, D 2s "˜0 0¸, D 3s "˜0 1¸(

A7)
In compact form, with 0 rˆc as a zero matrix with r rows and c columns s:

A8)
That is: ΠpsqXpsq "´CpsqY d psq`ΦpsqFpsq (A9) Now Y d = (Y md Y sd ) T can be obtained as a function of Y = (Y m Y s ) T by applying the time-delay expressions in Equations (2)-(4), with I ny as the n y = 2 order identity matrix: Y md Y sd¸"˜I n y¨D m 0 n yˆny 0 n yˆny I n y¨D s¸¨˜Y m Y s¸, or Y d " DY, Considering Equation (A2) and the corresponding slave equation, we obtain: Finally, replacing Equations (A10) and (A11) in Equation (A9), we have a complete description of the system: ΠpsqXpsq "´CpsqDDpsqXpsq`ΦpsqFpsq (A12)

Appendix B
Here we reproduce the classical implementation of the 4C control scheme, but maintaining the sign conventions as in Figure 1, that is, with the same definition of f m " f h´fmc , f s " f e`fsc .
In this case, the master and slave systems exchange external forces f e and f h , instead of f m and f s in our proposal. Accordingly, Equation (5) results as follows: f mc " C m x m`C2 f ed`C4 x sd´C6 f h f sc "´C s x s`C3 f hd`C1 x md`C5 f e (B1) and the addmitance matrix definition for zero delay is now given by: with the following values: Λ 11 " pP s`Cs q¨p1`C 6 q´C 4 C 3 Λ 12 "´pP s`Cs q C 2´C4 p1`C 5 q Λ 21 " C 1 p1`C 6 q`C 3 pP m`Cm q Λ 22 "´C 1 C 2`p P m`Cm q¨p1`C 5 q (B4) ∆ Λ " pP m`Cm q¨pP s`Cs q`C 1 C 4 (B5) The transparency conditions are achieved by selecting: For C m , C s : C 4 psq "´pP m`Cm q , C 1 psq " P s`Cs , C 2 "´p1`C 6 q , C 3 " 1`C 5 (B6) Note that for ∆ Λ Ñ 0 , C 1 and C 4 must be not proper.