Force-Sensor-Less Bilateral Teleoperation Control of Dissimilar Master-Slave System with Arbitrary Scaling

This study designs a high-precision bilateral teleoperation control for a dissimilar master-slave system. The proposed nonlinear control design takes advantage of a novel subsystem-dynamics-based control method that allows designing of individual (decentralized) model-based controllers for the manipulators locally at the subsystem level. Very importantly, a dynamic model of the human operator is incorporated into the control of the master manipulator. The individual controllers for the dissimilar master and slave manipulators are connected in a specific communication channel for the bilateral teleoperation to function. Stability of the overall control design is rigorously guaranteed with arbitrary time delays. Novel features of this study include the completely force-sensor-less design for the teleoperation system with a solution for a uniquely introduced computational algebraic loop, a method of estimating the exogenous operating force of an operator, and the use of a commercial haptic manipulator. Most importantly, we conduct experiments on a dissimilar system in 2 degrees of freedom (DOF). As an illustration of the performance of the proposed system, a force scaling factor of up to 800 and position scaling factor of up to 4 was used in the experiments. The experimental results show an exceptional tracking performance, verifying the real-world performance of the proposed concept.


Introduction
Bilaterally teleoperated robotic systems can bring the perception and precision of direct manipulation into challenging and risk-intensive tasks in environments that may be hazardous or hostile for humans. In contrast to unilateral teleoperation where the command flow goes only from the master to the slave, bilateral teleoperation provides the operator with information about the slave manipulator in the form of force feedback, to assist in the coordination and decision-making processes. To broaden the application scope of teleoperation, arbitrary motion and forces scaling has been pursued by many researchers, but no rigorously stability guaranteed method have been shown to work in a multi-DOF system.
Currently, one of the most interesting applications for teleoperation lies in Learning from Demonstrations (LdD) applications with heavy-duty manipulators. LfD is an established technique in robotics, where a robot is taught to perform tasks by demonstrations from a human teacher. The robot can then repeat these tasks in even slightly varying conditions [1]. The key enabler for LfD applications with heavy-duty manipulators is teleoperation of asymmetric systems with motion and force scaling. Conventional kinesthetic teaching methods, an established method for providing teaching samples, cannot be applied for such manipulators due to the size and force limitations (workspace over 2 m and payload over 500 kg) [2]. Instead, teaching samples can be captured using teleoperation with motion and force scaling between the manipulators. Teleoperation has the advantage of an intuitive and efficient communication and operation strategy between humans and robots. Teleoperated demonstrations have been successfully used for LfD applications with promising results using 1:800 force scaling in [2] and using 1:1 scaling in [3], and [4].
In applications where heavy objects are handled or a great amount of force is required, hydraulic actuation has remained the most attractive solution due to its great power-to-weight ratio. Hydraulic actuators further have the benefits of simplicity, robustness, and low cost. However, control of such actuators is significantly challenged by their complex nonlinear dynamic behavior. When the actuators are used in articulated systems, the control design is further complicated by the associated nonlinear multi-body dynamics, and the overall dynamics can be described by coupled nonlinear third-order differential equations. Consequently, the constrained motion control of multiple degrees-of-freedom (n-DOF) hydraulic robotic manipulators has been a well-recognized challenge [5].
As an additional challenge to the above, contact force measurements are often required for contact control. In conventional applications, a 6-DOF force/torque sensor is often attached to the tip of the manipulator for this purpose. However, these force/torque sensors are expensive and prone to overloading and shocks, a situation frequently occurring with hydraulic heavy-duty manipulators [6]. Therefore, methods avoiding direct contact force measurements have become desirable.
Needless to say, teleoperation of hydraulic manipulators has been an extremely difficult problem due to unresolved challenges in their high-precision control [7,8,9,10]. However, due to recent advances in hydraulic manipulators' highprecision control and leaps in the state-of-the-art (see [5,11,6,12]), teleoperation of hydraulic manipulators is suddenly becoming a feasible and interesting field of study again. Moreover, time delay, a focused research topic, especially in the teleoperation of extraterrestrial systems [13], can be alternatively addressed for terrestrial applications in the advent of 5G cellular networks with ultra-low latencies [14]. Terrestrial applications are within the author's main scope.
In this paper, we target an asymmetric bilateral teleoperation system comprised of a commercial haptic master manipulator and a hydraulic slave manipulator. The system has notable asymmetry between the manipulators due to substantial differences between the dynamics of the master and slave manipulators. Due to this asymmetry, handling motion and force scaling in the teleoperation architecture becomes necessary. The current state-of-the-art in teleoperation control has been focusing on purely electrical manipulators in symmetrical configurations, including multi-master or multi-slave setups, shared control, or dealing with time delays [15,16,17,18]. The existing methods for teleoperation of hydraulic manipulators have mainly relied on linear control theory and system linearization [8,9,10,19]. However, these methods have limitations in teleoperation of complex, highly nonlinear, and asymmetric systems. In contrast, the adaptive teleoperation scheme proposed by Zhu and Salcudean in [20] was reported to be capable of addressing nonlinear dynamics of asymmetric master and slave manipulators with arbitrary motion and force scaling. However, experiments with only a 1-DOF symmetrical system were presented. Moreover, both manipulators were equipped with force sensors.
In the present study, the results of [20] and [21] are used as the foundation for designing high-precision bilateral teleoperation control for significantly asymmetric systems. In [21], preliminary attempts for full-dynamics-based (and high-precision) bilateral teleoperation for an asymmetric hydraulic/electric system were demonstrated, while in [22] artificial constraints in the task space were implemented. However, sufficient stability analysis and theoretical discussions were not included.
To improve the preliminary theory and control performance reported in [21], the following distinguishable contributions are demonstrated in the present study. 1) We propose a master manipulator contact force estimation by using joint control torques and estimated manipulator dynamics. A solution to a computational algebraic loop formed around the actuation and force estimation is proposed. 2) We propose a novel method for estimating the exogenous force of the human operator. 3) Stability of the overall control design is rigorously guaranteed with robustness against an arbitrary time delay.
With the control theoretical developments described above, the experiments demonstrate significant improvements in relation to our preliminary study [21]. The experiments with a 2-DOF system with a force scaling ratio of up to 800 and a position scaling ratio of up to 4, in lieu of the 1-DOF experiments in [20], serve a critical step toward practical 6-DOF applications.
The rest of this paper is organized as follows: Section 2 presents the mathematical preliminaries. Section 3 discusses control of the master manipulator, while Section 4 discusses control of the slave manipulator. Section 5 presents the teleoperation scheme and discusses properties of the teleoperation method. Section 6 presents the experimental system and results. Finally, conclusions are drawn in Section 7.

Mathematical Preliminaries
Let an orthogonal coordinate system (i.e., a frame) {A} be attached to a rigid body. Then, the linear/angular velocity vector A V ∈ R 6 and the force/moment vector A F ∈ R 6 of the rigid body, expressed in frame {A}, can be expressed as [23]: are the linear and angular velocity vectors of frame {A}, expressed in frame {A}, and A f ∈ R 3 and A m ∈ R 3 are the force and moment vectors that are being measured and expressed in frame {A}. Transformation of linear/angular and force/moment vectors between two frames, attached to a common rigid body, namely {A} and {B}, can be expressed as [23]: where A U B ∈ R 6×6 is a force/moment transformation matrix that also transform velocities between frames {A} and {B}. The dynamics of a freely moving rigid body, expressed in the fixed rigid body frame {A}, can be defined as where M A ∈ R 6×6 is the mass matrix, C A ∈ R 6×6 the Coriolis and centrifugal terms, and G A ∈ R 6 the gravity vector of the rigid body. For a detailed formulation of M A , C A and G A , readers are referred to [23].

The Master Manipulator
Phantom premium 3.0/6DOF, a commercial haptic manipulator without any modifications to the hardware, has been chosen to act as the master manipulator in this study. It possesses 6-DOF manipulability and force-feedback along each individual DOF, with a workspace mimicking human arm motion pivoting from the shoulder. For this study, we developed a new control system for the manipulator to rigorously address the dynamics of the lightweight manipulator. As a challenge to the control design, the manipulator lacks force/torque sensors. Therefore, human operator contact force estimation is required. Without loss of generality, we consider manipulation and force perception only within a specific 2-DOF plane (by using joint 2 and joint 3), while the rest of the DOFs are locked with software. Benefits of using a commercial haptic device as the master manipulator is the ease of implementation for wide range of applications. Frame {S tcp } is assigned to the tool center point (TCP) of the master manipulator (see Fig. 1) and the external forces resulting from the dynamics of the human operator as well as the exogenous operating hand force are estimated and expressed in this frame. Orientation of frame {S tcp } is aligned with handle of the master manipulator and the handle is held horizontal as shown in Fig. 1 by position control at the spherical wrist.

Concept of the Virtual Decomposition Control
To design the intended high-precision teleoperation for complex asymmetric system, the study takes advantage of a novel virtual decomposition control (VDC) approach (see [23,24]). The method is developed especially for controlling complex robotic systems, with a number of significant state-of-the-art control performance improvements with robotic systems (see, e.g., [20,25,6,5,11,21,12]). As a key feature, VDC enables to virtually break down complexity of the original system to a set of manageable modular subsystems [23,25] such that the control design and stability analysis can be performed locally at the subsystem level without imposing additional approximations. This allows, e.g., that changing the control (or dynamics) of a subsystem does not affect the control equations of the rest of the system [23].
The subsystem-dynamics-based control design philosophy in VDC originates from two unique concepts, namely virtual stability and virtual power flows (VPFs); see Appendix B. The VPFs uniquely define the dynamic interactions among the subsystems such that the virtual stability of every subsystem ensures that a positive VPF is connected to its corresponding negative VPF in the adjacent subsystem (and vice versa). Thus, when every subsystem qualify as virtually stable, all the VPFs cancel each other out, eventually, leading to the stability of the entire system in the sense of Lebesgue integrable functions (see Appendix A). For more detailed information and additional benefits of VDC, see [23,5].

Kinematics
Relevant coordinate frames in terms of control of the master manipulator are shown in Fig. 1 Fig. 1.
Remark 1. Link 4 and 5 are virtually cut from the same rigid body, using design principles of the VDC approach, which allows separate computations later.
The independent joint velocity coordinates are denoted aṡ  Figure 1: Frame assignment of the haptic manipulator.
The respective joint angles q 2 and q 3 are shown in Fig. 2. Then, velocities of all links of the master manipulator can be determined using the geometrical transformation matrices between each frame with the independent joint velocity coordinates as Let the independent velocity coordinates at master manipulator's handle be where J m ∈ R 2×2 is the invertible Jacobian matrix of the master manipulator. Then, another mapping matrix can be defined as

Dynamics
Dynamics of each rigid body of the master manipulator can be determined using (3) and (5). According to [26], and under Assumption 1, dynamic model of the master manipulator can be expressed as where f m is the net reaction force from the master manipulator toward the human operator and will be defined later in more detail, τ m ∈ R 2 denotes the applied torques of the manipulator, and Assumption 1. Bearing friction of all the revolute joints of the master manipulator are zero.

Human Operator
Based on literature [27]; [20], and [28], sufficient accuracy for modeling the human operator can be achieved using a simple second-order linear time-invariant model. The following model is used here where M h ∈ R 2×2 , D h ∈ R 2×2 and K h ∈ R 2×2 are symmetric positive-definite matrices approximating the inertia, damping and stiffness of the arm of the human operator, respectively; while f m ∈ R 2 , appeared first in (10), denotes the net force vector, exerted by the master manipulator toward the operator, and f * h ∈ R 2 denotes the exogenous force vector actively generated by the operator. The position of the arm of the operator is denoted by x h ∈ R 2 , whileẋ h ∈ R 2 andẍ h ∈ R 2 denote the first and second time-derivatives of the position vector, respectively, subject toẋ In [16] and [29], it was suggested that the exogenous force of the operator could be estimated using a fast parameter adaptation function. This differs from the approach in [20], where a switching term with a constant force, instead of an estimate, was used to ensure stability. The precise expression of the exogenous force, denoted f * h ∈ R 2 in (14), would necessarily involve research on complex human motor neuron actions. In this paper, we describe this exogenous force as a general linear-in-parameter form as where Ψ(t) is a time-variant matrix and p is a parameter vector. We treat vector p as constant by moving all time-variant properties into Ψ(t). (16) is quite general. It covers the expressions used in [16] and [17], in which Ψ(t) = 1 are used. Most importantly, this expression takes the same form commonly used in neural networks, allowing flexible incorporation of basis radial functions into machine learning mechanisms.

Remark 2. Note that expression
With more elegant design of Ψ(t), for example, muscle activation measured by electromyography could be used for the intent force modeling.

Control of the Master Manipulator with a Human Operator
For accurate control of the master manipulator, dynamics of both the manipulator itself and the human operator need to be addressed together. The required control law must therefore define required contact force towards the human operator.
The estimated human operator exogenous force is written aŝ wherep is an estimate of the parameter vector. The time-invariant parameters are estimated using the following parameter adaptation law aṡ T as well as its upper and lower bounds, respectively, ρ i is the adaptation gain of the ith element of p = [p 1 ,p 2 , ...,p i , ...] T , Ψ i (t) denotes the ith column of the time-variant matrix Ψ(t), and V mr ∈ R 2 denotes the required velocity at the tip of the master manipulator, expressed in frame {S tcp }.
In Section 3.3, dynamics were calculated using the measured independent joint velocity vectorq m . However, since the proposed control method is velocity based, we need to define the required velocities in Cartesian space. Let V md ∈ R 2 be the desired velocity of the tip of the master manipulator, to be defined later in Section 5. Then, the required velocity vector, V mr ∈ R 2 , is designed as where A ∈ R 2×2 is a diagonal positive-definite gain matrix, andf m denotes a filtered estimate of the forces of the master manipulator to be determined later in this section. Computeq whereq mr ∈ R 2 denotes the required counterpart ofq m and V mr ∈ R 30 denotes the required counterpart of V m .
Remark 3. The second term in right hand side of (19) acts as a local force feedback term within the control design.
The linear parametrization of the required rigid body dynamics can be written according to [23] as Interested reader is referred to the formulation of the regressor matrix Y A ∈ R 6×13 and the parameter vector θ A ∈ R 13 in [23]. Using A ∈ {B 11 , B 12 , B 21 , B 22 , O 1 }, dynamics of each rigid body can be calculated with (22) as Furthermore, dynamics of the human operator are calculated with a similar linear parametrization form as Then, control equations for the master manipulator can be defined as where K m ∈ R 2×2 is a positive-definite gain matrix. The last term in (25) is a velocity feedback term used to ensure the control stability.

Force Estimation
The net reaction force from the master manipulator toward the operator can be estimated using the known dynamics of the master manipulator as a base. This method is similar to the inverse dynamics based estimation methods, described in [30]. The main difference here is that the estimated actuator torque and applied torque are calculated based on the inverse dynamics, yielding that the external force can be estimated in addition to the mere collision detection, possible with the simpler method. The force estimate can be expressed aŝ where τ m is the master robot control input defined in (25), and τ mm is the estimated master robot dynamics, defined as whereq m andq m are estimates ofq m andq m , respectively, obtained by differentiation from the measured joint angles q m . The filtered estimate of the master manipulator force vector is obtained usinġ where C ∈ R 2×2 is a diagonal positive definite matrix.

Computation Algorithms
Differentiating (19) and expressingV mr as an affine function ofḟ m yieldṡ where A 1 (t) ∈ R 2×2 is a known matrix and B 1 (t) =V md ∈ R 2 is a known vector (which will be given in Section 5), andḟ m ∈ R 2 is a vector to be specified later in this subsection. Using (20) and (21), it follows from (22)- (26), that where A 2 (t) ∈ R 2×2 is a known matrix and B 2 (t) ∈ R 2 is a known vector. Then, it follows from (28)˙f where B 3 (t) ∈ R 2 is a known vector. The existence of a computational algebraic loop can be clearly seen in (31).
To ensure numerical stability, we must have This means both C and A must be restricted. Finally,ḟ m can be computed from (31) aṡ Onceḟ m is obtained,f m in (19) can be computed using integration withḟ m (0) = 0.

Stability
Substituting (25) and (14) into (10) yields Then the non-negative function for the master manipulator is chosen as The time-derivative of the non-negative function in (35) is obtained using (34), (18) and the skew-symmetric properties of C * m aṡ Theorem 1. Analyzing the master manipulator (10) with the human operator (14) subject to control (25) with estimated exogenous operator force using adaptation law (18), it yields The proof directly follows (35) and (36). For the concept of L 2 and L ∞ stability (having similarities to Lyapunov functions method ), see Appendix A.

Slave Manipulator
A commercial HIAB-031 hydraulic manipulator is chosen to act as the slave manipulator of the teleoperation system. The manipulator is retrofitted with fast hydraulic servo valves, pressure transducers to measure cylinder chamber pressures and high accuracy incremental encoders to measure joint angles. Although the manipulator is retrofitted, it does not have force/ torque sensor at the TCP. Consequently, a force-sensor-less control method with external force estimation is used for the slave manipulator as was the case with the master manipulator.
In the experiments, manipulation and force perception is considered within the same 2-DOF plane as with the master manipulator. The extension cylinder and rotation of the boom was mechanically locked. Fig. 3 a illustrates the slave manipulator and shows several important frames of the manipulator. Frame {B s } is fixed to the base of the slave manipulator, frame {O 2 } is attached to the last link of the slave manipulator and frame {G} is attached to the tip of the slave manipulator and has the same orientation as frame {O 2 }. Frame {C} has the same origin as frame {G}, but is aligned with frame {B s }.
As discussed in Section 3.1, VDC enables modularity in the control design. Consequently, the slave manipulator can be considered as a subsystem (with its own local subsystems) of the overall system. Stability-guaranteed constrained motion control of the manipulator is described in [6]. To incorporate the control system designed in [6], control equations of the last object need to be adjusted while rest of the control system is kept identical to that of [6]. Fig. 3 presents the decomposed structure of the slave manipulator, with the re-used control design circled by a dashed line.

Object 2 -Kinematics and Dynamics
Let the linear/angular velocity vector TO2 V ∈ R 6 at the driven VCP of Object 2 be known from the kinematic chain through the previous subsystems (see [6]). Then, kinematic transformations among the frames in Object 2 (see Fig. 4) can be written as Next, dynamics of the environment are defined. In this work, we assume flexible environment with dynamics described by second-order linear time-invariant model [20] as where M e ∈ R 2×2 , D e ∈ R 2×2 and K e ∈ R 2×2 are symmetric positive-definite matrices approximating the inertia, damping and stiffness of the environment, respectively, and x s ∈ R 2 denotes the tip position of the slave manipulator, expressed in frame {B s }, subject toẋ s = V s with Then, dynamics of the environment can be included on the slave manipulator as where The net force/moment vector (rigid body dynamics) O2 F * of Object 2 can be written in view of (3) as  Figure 3: (a) The slave manipulator, (b) Virtual decomposition of the slave manipulator (c) simple oriented graph (SOG) of the slave manipulator. The circled area in the SOG represents subsystem of the slave manipulator, for which the control has been designed in [6]. and, eventually, the force balance (i.e. force resultant) equation of Object 2 can be written as

Object 2 -Control
Let the required velocity of the slave manipulator be designed as where V sd ∈ R 2 is to be defined in Section 5, andf s ∈ R 2 is obtained fromf s using a first order filter as˙f andf s is obtained using (15) in [6].
Required piston velocities of the slave manipulator are then redesigned from (87) in [6] into where J −1 x ∈ R 2×2 is the invertible Jacobian matrix of the slave manipulator, defined in [6].
Then, in view of (38), (39) and (41), the required linear/angular velocity vectors in Object 2 can be written as The required contact force of the slave manipulator is designed as Finally, using (22) and (42)-(45) the required control laws for Object 2 dynamics can be written as In line with (22), Y O2 θ θ θ ∈ R 6 in (54) is the model-based feedforward compensation term for the rigid body dynamics and K O2 ∈ R 6×6 is a positive-definite velocity feedback matrix to ensure the control stability. By defining the estimated parameter vector θ θ θ O2 ∈ R 13 in (54) is updated aṡ where θ O2i is the i th element of θ θ θ O2 ; s O2i is the i th element of s O2 ; ρ i > 0 is the update gain; θ − O2i is the lower bound of θ O2i ; and θ + O2i is the upper bound of θ O2i .

Stability
The remaining system, for which the control was designed first in [6], qualifies as virtually stable according to Theorem 2 Theorem 2. Consider the system encircled by a dashed line in Fig. 3. The subsystem qualifies virtually stable with its affiliated vector A V r − A V , ∀A ∈ Ψ r and its affiliated scalar variables (f pir − f pi ) for the hydraulic cylinder i, ∀i ∈ 1, 2, where Ψ r contains rigid body frames of each rigid link and object of the remaining subsystem. A non-negative accompanying function for this system can be found as Proof. The proof for Theorem 2 can be obtained from the results of [6]. . This subsystem is virtually stable with its affiliated vector O2 V r − O2 V being a virtual function in both L 2 and L ∞ in the sense of Definition 3. This is because a non-negative accompanying function can be found such thaṫ holds with 0 γ s < ∞. Note that p TO2 is the virtual power flow by Definition 2 in the driven VCP of Object 2, and p G characterizes the virtual power flow between the end-effector and the environment while in constrained motion (i.e., σ f = 1).

Proof. See Appendix C.
Theorem 4. Considering (106) and Definition 3, the contact with the environment qualifies virtually stable. The non-negative accompanying function for the entire slave manipulator can be written by summing the individual functions from (59) and (61) as Then stability analysis for the remaining subsystems follows exactly as shown in [6], ultimately yielding stability of the entire slave robot. Then it follows that

Teleoperation
After individual velocity-based controllers for both master and slave manipulators (see (25) and [6]) of the teleoperation system have been designed, a scheme for connecting the manipulators can be designed. Connection between the two manipulators is made with a communication channel that virtually connects the manipulators together. This section designs bilateral teleoperation and specifies two design vectors V md and V sd , used in (19) and (46), respectively. Using position control δ = 1 in [20], V md and V sd can be designed as where κ p > 0 and κ f > 0 are position and force scaling factors for arbitrary motion/force scaling between the manipulators, Λ ∈ R 2×2 is a diagonal positivedefinite matrix, and P m ∈ R 2 and P s ∈ R 2 denote the position/orientation of the master and slave manipulator, respectively, subject toṖ m = V m andṖ s = V s . Furthermore,Ṽ m ,Ṽ s ,P m ,P s , are filtered values of V m , V s , P m , P s , respectively, obtained using the following first order filteṙ where X ∈ R 2 is the input signal andX ∈ R 2 is the filtered signal. The use of filtered variables in the two design vectors, (69) and (68), makes the required accelerationsV sd andV md , functions of V m , V s , f m and f s .

Tracking
Substituting (68) and (69) into (37) and (67), then subtracting the resulting error terms from each other yields It can be easily seen that (71) can be further written as andZ is obtained from Z using (70). Following Lemma 2.4 in [20] yields Eventually, it follows from (73), (74) and Lemma 1 in [20] that hold, which guarantees the L 2 and L ∞ stability of the velocity and position tracking of the teleoperation system.
Transparency of the teleoperation system can be clearly seen from (79). Within a limited frequency range, the filtered signals can be assumed to be approximately equal to their non-filtered counterparts. The last term on the right hand side of (79) is bounded to converge to zero. Then, transparency error can be described by the second term on the right hand side of (79). It comprises of velocity and acceleration dependent terms. The acceleration related term κ −1 f κ p A −1 C −1 acts as a virtual mass on the teleoperation system, while the velocity dependent term κ −1 f κ p A −1 C −1 Λ determines the damping of the teleoperation system.

Stability under time delay
time delay under teleoperation is a much investigated issue especially in space teleoperation. Although the focus of this study is in terrestrial applications, robustness against arbitrary time delay of the proposed method is discussed briefly. In [13] similar approach was used for longer and varying delays.
Without loss of generality, we consider a one-dimensional system in the stability analysis. The extension to multiple-dimensional systems can be proceeded accordingly. Due to the fact that both master and slave manipulators have independent stability-guaranteed controllers linked only by the communication channel, the stability under time delay can be analyzed by modifying (68) and (69) as where the communication channel is represented as pure time delay of T ; see Fig. 5. The stability under arbitrary time delay can be analyzed similarly to the method presented in [20]. Fig. 5 represents the teleoperation system under arbitrary time delay based on (82) and (83). In the figure, Z h is the operator dynamics defined in (14) (now considered one-dimensional) and Z e is environment dynamics; here approximated with second-order linear dynamics.
To analyze the effect of time delay on the system stability, transfer functions for both sides of the communication channel need to be defined. The transfer function for master side, from input D to output A (see Fig. 5) can be formed as where Following the same procedure, transfer function from the input B of the slave side to output C can be formed as Assume flexible environment with dynamics as where M e , D e and K e define the inertia, damping and stiffness of the environment, respectively. Then, the transfer function from B to C can be written as To guarantee stability under arbitrary time delay, the gain of each manipulator together with their respective local controllers must remain equal or smaller than one across the entire frequency spectrum. Thus, to ensure stability of the entire teleoperation system with arbitrary time delay (see Fig. 5), the following conditions need to be satisfied To satisfy the stability conditions in (91) and (92), the following relation must be satisfied for both the slave and master manipulators by substituting M , D and K with M e , D e and K e (the slave side), or M * h , D * h and K * h (the master side), respectively. Furthermore, to satisfy (93), It follows directly from the positive-definite properties of M , D, K, C and Λ that a ≥ 0 and c > 0 hold indefinitely. Consequently, it follows from (93)-(95) that b + 4C ΛAK (2ACM + 1) ≥ 0 (96) must hold to fulfill the stability conditions in (91) and (92).

Slave side
Master side

Experiments
This section evaluates the performance of the proposed teleoperation system. First, Section 6.1 addresses the system implementation issues. Then, Section 6.2 provides the experiments without time delay, followed by the experiment with time delay in Section 6.3.

Experiment description
The experimental implementation comprises four main components, visualised in Fig. 6, which are the electric master manipulator (Phantom Premium 6DOF/3.0L haptic device), the host computer for the master manipulator, the real-time computer and the hydraulic slave manipulator (HIAB-031 manipulator). The two-DOF hydraulic manipulator (in Fig. 7) has a maximum reach of approximately 3.2 m, and a payload of 475 kg is attached to its tip. For the real-time control system, the following components were used: a DS1005 processor board, a DS3001 incremental encoder board, a DS2103 DAC board, a DS2003 ADC board, and a DS4504 100 Mb/s ethernet interface. The remaining hardware implementations can be found in [6], [11] or [21]. Control computations have been run with 500 Hz frequency. The communication channel parameters of each experiment are shown in Table 1. Teleoperation control between the master and slave manipulators was engaged by pressing a pushbutton and disengaged by releasing the button. At first, the slave manipulator was driven from free space to contact with the environment along the Cartesian y-axis; see Fig. 7 for the directions of the Cartesian coordinate system. After contact with the environment was established, the slave manipulator was driven along the surface of the pallets, while maintaining constant force against the environment. Finally, the slave manipulator was driven back to free space after approximately 0.5 m of sliding against the wooden pallets. The above described task was repeated without any time delay with two different sets of motion/force scaling parameters (κ p and κ f ) of the communication system; see Fig. 8. Then, an experiment with 80 ms one-way time delay in the communication channel was performed; see Fig. 9. In Figs. 8 and 9, the master manipulator data is shown in blue, and the slave manipulator data in red.

The experiments without time delay
In the first experiment (Fig. 8a), 1:1 position mapping was used between the manipulators (κ p = 1), while forces of the master manipulator were scaled up by a factor of κ f = 300; see Table 1. In the second experiment (Fig. 8b), 4:1 position scaling between the master and slave manipulator was used (κ p = 4), yielding 4 times larger movement of the slave manipulator compared to the movement of the master manipulator, along with force scaling by a factor of κ f = 800; see Table 1.
As Fig. 8 shows, accurate position (see the first and second rows) and force tracking (the third and fourth rows) between the master and slave ma-

Communication Channel
The control addressed in Section III The control addressed in Section IV The control addressed in Section V nipulators is achieved with different scaling parameters, as predicted by the theory, despite the inherent challenges of force control of hydraulic manipulator in the teleoperation system (see the discussion in Section 1). The forces in third and fourth rows presents the estimated contact forces of the slave-and master manipulator along the x-and y-axes, respectively. In the results, the master manipulator forces are scaled up by the respective scaling factor. Note that the slave manipulator contact forces are estimated from the cylinders' chamber pressures. Thus, some inaccuracies can exist in the measured contact forces (see [11] for more details). However, the operator is still able to effectively sense the contact forces between the slave and the environment, and excessive contact forces can be prevented. It is valid to mention that the proposed force-sensorless approach provides a practical solution for teleoperation of extremely powerful hydraulic manipulators, as conventional six-DOF force/torque sensors are fragile and prone to overloading [6].
Remark 4. Note that in the second experiment the transition from free space motion to contact motion, was done rapidly with a velocity of approximately 0.2 m/s. This is to demonstrate the stable behavior of the control system even with high velocity and rapid changes of system states.

Conclusions
This study presented force-reflected bilateral teleoperation using asymmetric electrical master and hydraulic slave manipulators. The control design takes advantage of the VDC approach, which allows us to design local subsystemdynamics-based controllers for the master and slave manipulators independently. Then, the teleoperation system was completed by designing the communication channel between the master and slave controllers as motivated by [20]. The teleoperation scheme provided unique features, such as arbitrary motion/force scaling between the manipulators, effectively enabling the connection of two very dissimilar manipulators.
The experimental results demonstrated the performance of the proposed method and showed excellent motion and force tracking between the manipulators.

Time [s]
Position tracking, y-direction Furthermore, robustness against an arbitrary time delay was demonstrated theoretically and experimentally. Similar to our previous studies [6,11], and [31], tracking performance improvements can be expected after rigorous application of a full parameter adaptation implementation and with tuning of the system parameters.
This paper advances force-reflected bilateral teleoperation control one more step toward practical applications and implements novel features to make it applicable to a large class of manipulators remotely operated over 5G cellular network. The theoretical and experimental studies in this paper concluded that the use of a force-sensor-less design for bilateral teleoperation is feasible.
Future work will focus on maximizing the system performances of position tracking and transparency, as well as forming basis functions for human operator exogenous forces through machine learning. Moreover, we intend to focus on expanding the experimental system to possess more degrees of freedom with the goal of achieving 6-DOF manipulation.
A L 2 and L ∞ Stability Definition 1 provides a definition for the Lebesgue space. Lemma 1 (a simplified version of Lemma 2.3 in [23]) provides that a system is L 2 and L ∞ stable with its affiliated vector x(t), being a function in L ∞ and its affiliated vector y(t), being a function in L 2 .
Lemma 2 provides an alternative to Barbalat's lemma.
Remark 5. As a distinction to Lyapunov approaches, Lemma 1 allows different appearances of variables in the non-negative function itself and in its timederivative. When all error signals are proven to belong to L 2 and L ∞ in the sense of Lemma 1, then asymptotic stability can be proven with Lemma 2, if the time-derivatives of all error signals belong to L ∞ . Note that s(t) = 0 is a special case that satisfies (100) in Lemma 1.

B Virtual Stability
The unique feature of the VDC approach is the introduction of a scalar term, namely the virtual power flow (VPF) [23]; see Definition 2. The VPFs uniquely define the dynamic interactions among the subsystems and play an important role in the definition of virtual stability [23], which is defined in a simplified form in Definition 3.

Definition 2 ([23]
). The virtual power flow with respect to frame {A} is the inner product of the linear/angular velocity vector error and the force/moment vector error as where A V r ∈ R 6 and A F r ∈ R 6 represent the required vectors of A V ∈ R 6 and A F ∈ R 6 , respectively.

Definition 3 ( [23]).
A subsystem with a driven VCP to which frame {A} is attached and a driving VCP to which frame {C} is attached is said to be virtually stable with its affiliated vector x(t) being a virtual function in L ∞ and its affiliated vector y(t) being a virtual function in L 2 , if and only if there exists a non-negative accompanying function where 0 γ 0 < ∞, P and Q are two block-diagonal positive-definite matrices, and p A and p C denote the virtual power flows (by Definition 2) at frames {A} and {C}, respectively.
Remark 6. In view of Theorem 2.1 in [23], when all subsystems qualify as virtually stable (in the sense of Definition 3), the L 2 and L ∞ stability of the entire system can be guaranteed in the sense Lemma 1.