Multilateral teleoperation under asymmetric time delays: L 2 stability and robustness

Abstract

Multilateral teleoperation is the general name given to haptic teleoperation with multiple robots and operators. It is the generalization of two robot bilateral teleoperation systems to N robots. One of the applications of multilateral teleoperation is haptic training, where multiple operators can move multiple master robots, and in consensus, they can control a single slave robot. This enables the operators to influence each other’s motion with their force inputs and a more experienced operator can instruct/train a novice operator for teleoperation. In multilateral systems, existence of multiple robots and multiple communication channels between robots poses difficulties in the analysis of stability and control design. Because of the distributed nature of the problem, and especially in the presence of time delays in the communication links between robots, stability of the networked system becomes harder to guarantee. This article presents a method for checking the L ₂ stability of multilateral teleoperation systems with N robots and symmetric or asymmetric time delays between them. Following this, a scalable four-channel-based distributed control law is proposed under the light of the proposed stability criterion which can guarantee delay-independent stability and enable high-performance haptic teleoperation. Furthermore, robust stability of the proposed control architecture is shown using μ analysis. Finally, theoretical results are validated with experiments that include comparisons of the proposed method with two-channel-based passive and absolutely stable systems.

Introduction

Haptic or bilateral teleoperation is the remote operation of robotic systems with haptic/force feedback to the operator from the remote side. Exchange of information takes place between two robots, namely a master and a slave robot. The master robot is moved by an operator with the application of force, and the motion of the master is copied or scaled by the slave robot and any forces acting on the slave are reflected to the operator. Position and force data are transmitted between the master and the slave on a communication network. The problem of communication time delay in bilateral teleoperation has been a subject of ongoing research for a long time. For the first time in the study of Anderson and Spong, ¹ it was shown that a bilateral teleoperation system could be made stable under time delay through passivation. The work of Niemeyer and Slotine ² proposed a simple way of utilizing wave variables to passivate the closed-loop teleoperation system. Lawrence³ showed that stability and transparency, representing the performance of the teleoperation system and which can be measured through hybrid parameters, ⁴ cannot be achieved at the same time. Four-channel architecture, with the transmission of both position and force signals from the master and slave robots to one another, was also proposed in the study by Lawrence ³ to maximize the transparency of a bilateral teleoperation system. However, it is known that, four-channel architectures can optimize transparency when there is no time delay but cannot achieve stability or transparency when there is time delay between the robots. Methods such as three-channel controllers ⁵ or communication disturbance observers ⁶ have been proposed to improve the transparency of teleoperation systems under time delay but still cannot guarantee stability under a wide range of environmental conditions. The presence of multiple position and force communication channels make stability more difficult to guarantee. To solve this problem, damping injection was proposed as an effective tool and as a trade-off between stability and transparency. ⁷ With damping injection, stability and transparency can be achieved with the caveat that the system becomes relatively slow and damped.

With the recent advances in networking technologies, multirobot systems, multirobot coordination, flocking and teleoperation have emerged as promising research fields. Dynamic consensus algorithms making use of graph theory have provided a breakthrough in the analysis and design of distributed and networked control systems. ⁸ Similarly, in the field of haptic teleoperation, Lee and Spong ⁹ and Sirouspour and Setoodeh ¹⁰ have proposed bilateral teleoperation architectures where one master operator can cooperatively control N robots working on a single environment. Also, Sirouspour ¹¹ has proposed a robust stability analysis and design method for cooperative teleoperation. Polushin et al. ¹² have proposed a cooperative teleoperation scheme with nonlinear small gain theorem for the stability of the networked system. Katsura and Ohnishi ¹³ have proposed a multilateral control system, as an extension of the four-channel architecture in bilateral teleoperation, with N master robots/operators to control a single slave robot for the purpose of haptic training. Katsura et al. ¹⁴ have extended that method to cooperative motion control. Multilateral teleoperation systems have drawn attention especially because of their application in haptic training. Khademian and Hashtrudi-Zaad ¹⁵ have proposed a robust multilateral control architecture, Tumerdem and Ohnishi ¹⁶ have utilized consensus algorithms for multilateral teleoperation on changing network topologies, and Tumerdem and Ohnishi ¹⁷ have proposed a solution for multilateral teleoperation under time delay with damping injection. However, Tumerdem and Ohnishi ¹⁷ do not provide an analytic stability proof for choosing the controller parameters, and furthermore, the proposed analysis method cannot accommodate asymmetric time delays, delays that are not identical on the communication channels. When the robots are remotely placed and there are different physical communication lines between them, it is almost impossible to guarantee symmetric delays between the robots, and such an assumption is rather limiting in real applications. Other architectures for stable teleoperation under time delay, mostly excluding the use of explicit force feedback channels, and methods for checking stability have also been proposed. For instance, an analytic method for checking the stability of multilateral teleoperation systems have been provided in studies of Mendez ¹⁸ and Shahbazi ¹⁹ by using n-port passivity. Following these papers, Huang and Lee, ²⁰ Kanno and Yokokohji ²¹ and Van Quang and Ryu ²² have used passivity to guarantee stability of multilateral teleoperation systems. Huang and Lee ²⁰ have proposed a consensus-based algorithm which does not make use of force channels and can provide force reflection only in steady state, which is not desirable in terms of performance. Kanno and Yokokohji ²¹ have made use of wave variables for stability, but the system has low force bandwidth and steady-state errors which actually is typical for wave variable methods. Van Quang and Ryu ²² have utilized a time domain passivity observer and have achieved better results, but again suffer from steady-state error and restrictions in performance due to forcing the teleoperation system to be passive. Wang et al. ²³ also employ a four-channel architecture but passivate the controller according to passivity criterion. As a result, the performance of the system deteriorates and there are large errors in velocity (used instead of position) tracking and force reflection. Sun et al. ²⁴ also use wave variables to guarantee the nonlinear stability of a multilateral system. However, the system in that article is different from the haptic training system which is the focus of this article as it assumes multiple slaves and a single master. Furthermore, the experiment results show large steady-state errors in contact. Chen et al. ²⁵ have proposed an adaptice controller and show the nonlinear stability of an n master m slave cooperative teleoperation system. Therefore, it is very hard to gauge the effectiveness and performance of that control system in a haptic training scenario. Furthermore, an adaptive control strategy can be risky to implement in real teleoperation experiments due to the fact that environment and human parameters can change quite fast and destabilize the system. To get better performance and still guarantee stability, for trilateral teleoperation (two masters and one slave), Razi and Hashtrudi-Zaad ²⁶ and Li et al. ²⁷ have made use of absolute stability criterion which is a less conservative stability criterion compared to passivity. In such a set-up, the environment and operator are still assumed to be passive but the closed-loop system does not have to be passive. However, the absolute stability criterion is also conservative, because it does not take the parameters of the environment and operator into consideration and thus limits the performance of the system. It can be further relaxed by taking environmental/operator parameters into consideration. ²⁸ For this purpose, Tumerdem and Ohnishi ¹⁷ have proposed an L ₂ small gain-based numerical analysis method for checking the stability of multilateral system with symmetric time delay, and Shahbazi et al. ²⁹ have presented an L ₁ small gain method for dual master single slave multilateral teleoperation. Furthermore, no paper in the literature has performed a stability analysis by making use of robust control theory for multilateral teleoperation systems under time delay. Sirouspour ¹¹ has proposed a μ analysis and synthesis method for cooperative teleoperation without time delay, and Khademian and Hashtrudi-Zaad ¹⁵ have proposed a robust stability analysis and H_∞ robust control method for dual user single master teleoperation without time delay.

There is still a need for guaranteeing the stability of N robot multilateral teleoperation systems with asymmetric time delays without conservative assumptions. For this purpose, in this article, we are proposing an L ₂ input–output stability criterion for checking the stability of multilateral teleoperation systems and a robust stability analysis method for dealing with environmental uncertainties. For the first time, we are providing a stability proof and relaxed conditions for designing four-channel-based multilateral controllers guaranteeing delay independent stability of multilateral systems by taking operator and environmental parameters into consideration as mentioned in the literature. ²⁸ This relaxed analysis method helps us employ four-channel-based controllers with high position and force gains and obtain better position tracking/force reflection results in haptic training experiments compared with the methods in the literature, which are limited due to absolute stability or passivity criteria and neglect environmental parameters. In the experiments section, to demonstrate this claim, we compare the performance of the proposed controller with controllers proposed according to absolute stability and passivity criteria. ²⁷ Furthermore, we show that the controller we propose can deal with asymmetric time delays between the robots. Because we make use of graph theory and ideas from dynamic consensus algorithms, stability analysis and control laws are scalable to teleoperation systems with arbitrary number of master and slave robots. Finally, in this article, we also propose, for the first time, a robust stability analysis technique for multilateral teleoperation under time delay (including asymmetric time delay), using parametric uncertainty representation and μ analysis and show that the proposed control architecture is capable of achieving robust stability. Thus, the novelty of this article is threefold: first, we propose a more relaxed criterion for checking the stability of multilateral teleoperation systems under asymmetric time delays; second, we propose force reflecting four-channel-based controllers using damping injection whose parameters can be determined by the proposed stability criterion; third, we propose a robust stability analysis method for multilateral teleoperation systems with N robots and show that the control architecture proposed here can guarantee robust stability.

The article is divided into following sections: in the second section, acceleration control for robotic manipulators with DC motors using disturbance observers is discussed. Third section provides the basic graph theoretical notation that will be used in the article. In the fourth section, multilateral teleoperation is introduced as an extension of four-channel bilateral teleoperation architecture and local damping injection is proposed to stabilize the multilateral teleoperation system under time delays. In the fifth section, a method based on small gain theorem is proposed to check L ₂ stability of the multilateral teleoperation systems. Using this analysis, controller parameters, including damping injection terms, are determined to guarantee that the proposed control system is L ₂ stable. In the sixth section, a robust stability analysis method for multillateral teleoperation is introduced and robustness of the proposed control architecture is shown. In the seventh section, experiment results validating the efficacy of the proposed method and discussions are provided. Eighth section concludes the article and discusses future directions.

Disturbance observer and acceleration control

Acceleration controllers with disturbance observers are effective motion controllers that are used as low-level controllers in this article. Disturbance observer has a simple structure: it estimates the disturbance forces acting on a DC motor by comparing the responses of the actual plant and a model of the plant. Figure 1 shows the structure of a disturbance observer implemented on a brushless DC linear motor. For such a system, where g_dis is the filter cut-off frequency, $I_a^{ref}$ is the current supplied, K_tn is the motor constant, M_n is the nominal mass of the motor, x is the encoder-based position measurement, $\dot{x}$ is the motor velocity obtained by numerical derivation and low-pass filtering, and the disturbance force estimate ${\widehat{F}}_{dis}$ on a motor can be found by ³⁰

Figure 1.

Disturbance observer.

${\widehat{F}}_{dis}=\frac{g_{dis}}{s+g_{dis}}(I_a^{ref}{K}_{tn}+g_{dis}{M}_n\dot{x})-g_{dis}{M}_n\dot{x}$

1

${\widehat{F}}_{dis}=\frac{g_{dis}}{s+g_{dis}}{F}_{dis}$

2

This disturbance estimate can be fed back to the system in addition to acceleration control law

${\ddot{x}}_{ref}={\ddot{x}}^{cmd}+K_p(x^{cmd}-x)+K_v({\dot{x}}^{cmd}-\dot{x})$

3

where ${\ddot{x}}_{ref}$ is the acceleration reference for the motor control system and x_cmd is the desired trajectory. Position response of a system trying to track a reference trajectory with this controller becomes

$x^{res}=\frac{1}{M_n{s}^2}\left(M_n{\ddot{x}}_{ref}-\frac{s}{s+g_{dis}}{F}_{dis}\right)$

4

Ideal disturbance compensation and trajectory tracking is achieved if g_dis → ∞ which implies $\frac{s}{s+g_{dis}}\to 0$ . This means that effect of the disturbances on the system response is eliminated within the bandwidth of the disturbance observer. Disturbance observers can also be used for the measurement of external forces. When they are modified for external force estimation, they are called reaction force observers. The reaction force observer equations are given as ³¹

${\widehat{F}}_{reac}=\frac{g_{reac}}{s+g_{reac}}(I_a^{ref}{K}_{tn}+g_{reac}{M}_n\dot{x}-F_{int}-f-D\dot{x})-\hspace{0.17em}\hspace{0.17em}{g}_{reac}{M}_n\dot{x}$

5

where D is the damping term due to viscous friction, F_int is the dynamic coupling force acting on the motor that is part of a multi–Degree of Freedom (DOF) robotic system, and f is the kinetic friction term and g_reac is the cut-off frequency for the reaction force observer. In effect, this estimator subtracts the known dynamic forces and disturbances acting on the system from the total disturbance measurement to extract the external forces.

Graph theory notation

Recently, especially following Olfati-Saber and Murray, ⁸ graph theory has become an effective tool in the analysis and design of controllers for multi-agent systems that are connected to each other with communication networks. Although our results could also be explained without any reference to graph theory, the use of it makes it easier to scale our results to systems with N robots and explain different control/communication architectures with a unified/simple representation. For designing a multilateral control system, we have to analyse the information flow among multiple robots. The information flow between these robots can be represented by graphs. And by using tools from graph theory such as adjacency, degree and Laplacian matrices, it becomes possible to design and analyse scalable control laws for multirobot systems. Figure 2 shows examples of such information graphs, where Figure 2(a) represents a strongly connected architecture with all robots communicating with each other. Figure 2(b) represents a ring architecture: each robot is able to communicate with one other robot adjacent to it and the robots complete a ring among themselves. Figure 2(c) shows a hub architecture with one robot acting as a communication hub between other robots. In this article, we will assume a strong connection among robots, but it is possible to design multilateral controllers with different topologies. Now let us assume that the graph $\mathcal{G}$ = $(\mathcal{V},\mathcal{E},\mathcal{A})$ is a directed graph with n nodes: $\mathcal{V}$ ={v ₁, …, v_n }. The graph also has directional edges connecting the node: $\mathcal{E}\subseteq \mathcal{V}\times \mathcal{V}$ . The definition of a single edge between nodes i and j is given as follows: e_ij = $\left(v_i,v_j\right)\subseteq \mathcal{E}$ . Adjacency matrix of a graph $\mathcal{G}$ is $\mathcal{A}$ which is composed of the following elements a_ij

$a_{ij}=\{\begin{array}{cc}1\hspace{0.17em}\hspace{0.17em}& if\hspace{0.17em}\hspace{0.17em}(v_i,v_j)\in \mathcal{E}\\ 0\hspace{0.17em}\hspace{0.17em}& else\end{array}$

6

Figure 2.

Information graphs: (a) strongly connected; (b) ring; (c) hub.

The degree matrix D = [δ_ij ] is a matrix containing the degree (number of connections to other nodes) of every node on the diagonals and is defined as follows

$D=diag\left({\displaystyle \sum _{j=1}^N}{a}_{ij}\right)$

7

Finally, Laplacian matrix is used in our control design and is the difference between the degree matrix and the adjacency matrix: $\mathcal{L}=D-\mathcal{A}$ . When time delays are present on the informational links between the nodes, we can represent this using a delayed adjacency matrix in the form used by Münz et al. ^32,33

${\mathcal{A}}_{\tau }=[a_{ij}{e}^{-s{T}_{ij}}]$

8

where T_ij is the delay in the information channel from node i to node j in seconds, and e ^−s_Tij is the delay element written in frequency domain. For the delayed system, the degree matrix is the same as the no delay case, and the Laplacian matrix can be written as ${\mathcal{L}}_T=D-{\mathcal{A}}_T$ .

Controller design for multilateral teleoperation systems with time delay

In conventional four-channel-based multilateral teleoperation, ¹³ position and force measurements by each robot must be transmitted to the others. In the special case of haptic training, there are multiple masters and a single slave. In this case, master robots and the slave should be doing the same motion at all times and should communicate their position data to one another so that they can synchronize their motion. Furthermore, since this is a haptic teleoperation scheme, the forces acting on the slave during contact should be reflected to the master operators. Hence, all robots should also share their force measurements at all times. The distribution of the force reflection among the master robots is dynamically realized by the operators. The operators can influence the motion of other master robots as well as the slave through force input, thus enabling haptic training. These goals can be encapsulated by the following equations, which mean that during contact, the sum of the force estimates (from disturbance observers) should be 0, implying collaborative force reflection, and the positions of the robots at all times should be equal, implying synchronization ¹³

$F_1+F_2+\cdots +F_N=0$

9

$x_1=x_2=\dots =x_N$

10

Here, F_i are scalar force measurements and x_i are the position measurements of each robot. For these goals, N acceleration control laws can be written for each robot in a distributed manner

${\ddot{x}}_1^{ref}={\ddot{x}}_1^{cmd}+C_p(x_1^{cmd}-{\delta }_1{x}_1)$

11

${\ddot{x}}_2^{ref}={\ddot{x}}_2^{cmd}+C_p(x_2^{cmd}-{\delta }_2{x}_2)$

12

$\begin{array}{c}⋮\\ {\ddot{x}}_N^{ref}={\ddot{x}}_N^{cmd}+C_p(x_N^{cmd}-{\delta }_n{x}_N)\end{array}$

13

Here, C_p (s) = K_p + K_vs is a PD compensator transfer function, δi/δii are the diagonal elements of D. The position references of each robot is determined by the position of other robots and the connection matrix $\mathcal{A}$

$x^{cmd}=\mathcal{A}x$

14

where x is the combined position vector. Acceleration feedforwards for the robots are identical and can be written as

${\ddot{x}}_1^{com}={\ddot{x}}_2^{com}=\dots ={\ddot{x}}_N^{cmd}=C_f(F_1+F_2+\cdots +F_N)$

15

Here, C_f is a gain for force feedforward and we use it as 1 in this article, although it can be increased or decreased to change the effective inertia of the overall control system.

In our control design and analysis, we require strong informational connection among the robots. A strongly connected graph can be seen in Figure 2(a). Application of this structure to multilateral teleoperation implies that all robots pass their position and force measurements to all the other robots. A strongly connected graph with N nodes has the adjacency matrix

$\mathcal{A}=\left[\begin{array}{cccc}0& 1& \cdots & 1\\ 1& 0& & ⋮\\ ⋮& & \ddots & 1\\ 1& \cdots & 1& 0\end{array}\right]$

16

And the degree matrix is diagonal with N − 1 in the diagonals. With such a graph, we can write a distributed controller of the form

${\ddot{x}}^{ref}=-C_p\mathcal{L}x-C_f(\mathcal{I}+\mathcal{A}))F$

17

where F is the force vector containing all the force measurements and $\mathcal{I}$ is the identity matrix. This control law is distributed, and all the robots have their individual controllers implemented locally. In Figure 3, the local controller for each robot can be seen as a block diagram along with the communication channels to other robots.

Figure 3.

Block diagram of the local controllers in the distributed multilateral control scheme.

The various controller transfer functions in each robot i are chosen as

$C{1}_i=C_p(s)=K_p+K_{v}s$

18

$C{2}_i=C_f$

19

$C_i=K_p+K_{v}s$

20

Here, all the position measurements transmitted across the communication channels pass through a transfer function C1 _i which is chosen to be the same as the PD transfer function C_p . Similarly, the force measurements pass through C _2i which is chosen as C2 _i = C_f = 1. C_i is the local compensator transfer function for each robot. It can be seen that the outputs y of each robot contain both the position and force information and these are fed as inputs to the other robots according to the network topology represented by the interconnection matrix $\mathcal{A}$ . In the presence of time delays among the robots, we can represent the time delay from robot i to robot j with T_ij in seconds. So, in our analysis, there are no constraints on the time delay between the robots except that they are constant. The delays between robots do not have to be equal or even symmetric, which leads to a powerful, general and scalable stability result. If there are delays between robots, then the controller equation will be revised as

${\ddot{x}}^{ref}=-C_p{\mathcal{L}}_{\tau }x-C_f(\mathcal{I}+{\mathcal{A}}_{\tau })F$

21

where

${\mathcal{A}}_{\tau }=\left[\begin{array}{cccc}0& e^{-s{T}_{12}}& \cdots & e^{-s{T}_{1N}}\\ e^{-s{T}_{21}}& 0& \cdots & ⋮\\ ⋮& & & e^{-s{T}_{(N-1)N}}\\ e^{-s{T}_{N1}}& \cdots & e^{-s{T}_{N(N-1)}}& 0\end{array}\right]$

22

and

${\mathcal{L}}_{\tau }=\left[\begin{array}{cccc}N-1& -e^{-s{T}_{12}}& \cdots & -e^{-s{T}_{1N}}\\ -e^{-s{T}_{21}}& N-1& \cdots & ⋮\\ ⋮& & & -e^{-s{T}_{(N-1)N}}\\ -e^{-s{T}_{N1}}& \cdots & -e^{-s{T}_{N(N-1)}}& N-1\end{array}\right]$

23

Stability and performance of the distributed control system are negatively affected by the time delay elements in the Laplacian and adjacency matrices, and actually, the conventional system becomes unstable as will be seen in the following section. To overcome this problem, we compensate for the delay elements in the interconnection by using local damping in the robots. These terms are not affected by time delay and dissipate any excess energy created by the delay elements in the communication channels. We will prove that this approach will stabilize our system in the next section. The local damping injection–based distributed multilateral control law can be realized by choosing the control parameters as

$C{1}_i=C_p(s)=K_p+K_{v}s$

24

$C{2}_i=C_f$

25

$C_i=K_p+K_{v}s+b_{i}s$

26

where b_i are the local damping terms added to the robot controllers. As a result, the distributed control law becomes

${\ddot{x}}^{ref}=-C_p{\mathcal{L}}_{\tau }x-B\dot{x}-C_f(\mathcal{I}+{\mathcal{A}}_{\tau })F$

27

where B = diag(b_i ).

In the next section, we will introduce the small gain theorem and input–output stability concepts to prove the L ₂ stability of the proposed method. Here, it is important to note that the controller parameters C_p (s), C_f and b_i are identical for each robot. As will be shown in the next section, these values will be chosen based on the robot/environment transfer function with the largest H_∞ gain.

Delay-independent L ₂ stability analysis for multilateral teleoperation

Input–output L ₂ stability

In our analysis, we start with the assumption that the signals in our system are all in the L ₂ signal space, which means that they satisfy

${\displaystyle {\int }_0^{\infty }}{\left|y(t)\right|}^{2}dt<\infty$

28

And their L ₂ norm is

$\left|\right|y{\left|\right|}_2={\left({\displaystyle {\int }_0^{\infty }}{\left|y(t)\right|}^{2}dt\right)}^{1/2}$

29

The extended L ₂ space L _2e is a space of functions satisfying

${\displaystyle {\int }_0^T}{\left|y(t)\right|}^{2}dt<\infty$

30

where T is finite. Assume that L _2e (U) is an input signal space, L _2e (Y) is an output signal space and S is a mapping y = S(u) from the input signal space to the output signal space. Here, S can represent our system. S is said to be finite gain L ₂ stable if its gain γ ₂ satisfies the following condition for T ≥ 0

$\left|\right|S(u{)||}_2\le {\gamma }_2||u{\left|\right|}_2+b_2$

31

Here, b ₂ is an arbitrary constant. This implies that the output is bounded if the input is bounded.

We can transform systems connected to each other with feedback connection to input–output form. Figure 4 shows a system connected in input–output form and equations representing this system can be written in vector form

Figure 4.

Input–output connection.

$U=E-HY$

32

$Y=GU$

33

or in scalar form

$u_i=e_i-{\displaystyle \sum _{j=1}^N}{H}_{ij}{y}_j$

34

$y_i=G_i(u_i)$

35

Here, G and H are matrices and G is diagonal with G_i as its diagonal elements. Small gain theorem states that any input–output connected system will be L ₂ input–output stable if both G and H are stable meaning that the gains γ ₂(G) and γ ₂(H) are bounded, and if together they satisfy the small gain condition

${\gamma }_2(G)\cdot {\gamma }_2(H)<1$

36

L ₂ stability analysis for multilateral teleoperation systems

Stability analysis by checking the poles of a multivariable system like the multilateral control system we are interested in can be increasingly difficult as the number of robots and the amount of information exchange among them increase, and this becomes even more difficult when time delay is present in the communication channels. Therefore, using the input–output approach to check the stability for a multivariable system is a common approach. In our stability analysis, we represent the multilateral system in the input–output framework. Here, we also assume that the environment and the operators are a part of the closed loop teleoperation system. For this purpose, we model the environment and human forces acting on the robots using the following equation, which has been commonly employed in stability analysis for bilateral teleoperation systems ^3,5,28

$\text{Operator/Environment}\hspace{0.17em}\hspace{0.17em}i\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}:F_i=F_i^{*}-Z_i{V}_i$

37

Here, $F_i^{*}$ are exogenous inputs independent of the system states, V_i are robot velocities in contact with the human/operator, and Z_i are the human/operator impedances mapping velocity to force

$Z_i=\frac{k_{zi}}{s}+b_{zi}+m_{zi}s$

38

Here, k_zi is the environmental/operator stiffness parameter, b_zi is the damping parameter, and m_zi is the mass parameter. It is assumed in bilateral teleoperation literature that exogenous forces are L _2e signals and bounded. ³ By modelling the operators and environments as mentioned, it becomes straightforward to transform the multilateral teleoperation system to the input–output form as in Figure 4. The system block G is a transfer matrix that contains the transfer functions G_i of each robot and the operator/environment the robot is in contact with. System block H is the interconnection matrix between these robot transfer function outputs. As mentioned before, in multilateral teleoperation, feedback interconnection between the robots is represented by an adjacency matrix for a strongly connected graph. Time delays in the communication channels are an integral part of this framework in the form of the delayed adjacency matrix which makes up the H block. The equations representing this system are

$e_i={\displaystyle \sum _{j=1,j\ne i}^N}{e}^{-T_{ij}s}{F}_j^{*}+F_i^{*}$

39

$u_i=e_i-{\displaystyle \sum _{j=1}^N}{H}_{ij}{y}_j$

40

$y_i=G_i(u_i)$

41

$G_i=\frac{(C{1}_i(s)-C{2}_i{Z}_{i}s)}{m_i{s}^2+(b_i+Z_i(s))s+{\delta }_i{C}_p(s)}$

42

$H={\mathcal{A}}_{\tau }$

43

And G is of the form

$G=\left[\begin{array}{cccc}{G}_1& 0& \cdots & 0\\ 0& \ddots & 0& 0\\ ⋮& \cdots & G_{N-1}& 0\\ 0& \cdots & 0& G_N\end{array}\right]$

44

where G_i is the transfer function for a robot which is controlled with the four-channel-based control law explained in the previous section and which is in contact with an environment that has an impedance Z_i .

The L ₂-induced gain of a matrix operator S can be obtained by the H_∞ norm, which is the supremum of the maximum singular value over the frequency domain

${\gamma }_2(S)=\underset{\omega >0}{\text{sup}}\hspace{0.17em}{\sigma }_{\text{max}}(S(j\omega ))$

45

${\sigma }_{\text{max}}(S(j\omega ))=\sqrt{{\lambda }_{\text{max}}\{S^{*}(j\omega )S(j\omega )\}}$

46

Where S ^*(jω) is the conjugate transpose of S(jω) and λ are the eigenvalues of a matrix. So, according to small gain theorem, if γ ₂(G(ω)) ⋅ γ ₂(H(ω)) < 1, ω > 0 can be guaranteed, then the multilateral teleoperation system becomes L ₂ stable. In order to obtain the sufficient conditions to guarantee stability, we will try to guarantee locally that each subsystem is stable. So we decouple the overall stability problem into two parts: making sure that individual system plants satisfy a certain gain condition and that the feedback interconnection satisfies a certain gain condition. Since the feedback connection between the robots can be represented by the delayed interconnection matrix ${\mathcal{A}}_{\tau }$ , gain of the interconnection is determined by

${\sigma }_{\text{max}}(H(j\omega ))=\sqrt{{\lambda }_{\text{max}}\{{\mathcal{A}}_{\tau }^{*}{\mathcal{A}}_{\tau }}\}$

47

It can be shown that ${\mathcal{A}}_{\tau }^{*}{\mathcal{A}}_{\tau }$ is Hermitian. The eigenvalues of a Hermitian matrix are always real which in this particular case mean that they are independent of time delay. By the Gershgorin circle theorem, we have that the eigenvalues are always constrained to be in a circle with the radius determined by the rows of the interconnection matrix. Let $\mathcal{H}={\mathcal{A}}_{\tau }^{*}{\mathcal{A}}_{\tau }$ , then the eigenvalues of $\mathcal{H}$ satisfy

$|\lambda -h_{ii}|\hspace{0.17em}\le {\displaystyle \sum _{j\ne i}}\left|h_{ij}\right|$

48

$|\lambda -{\delta }_i|\hspace{0.17em}\le (N-1)(N-2)$

49

$|\lambda -(N-1)|\hspace{0.17em}\le (N-1)(N-2)$

50

$\left|\lambda \right|\hspace{0.17em}\le {(N-1)}^2$

51

from which we obtain a bound on the maximum singular value of the interconnection matrix H(jω)

${\sigma }_{\text{max}}(H(j\omega ))\le (N-1)$

52

Thus, we obtain the result that the interconnection matrix of a strongly connected graph with time delay has an H_∞ norm/L ₂ gain equal to N − 1. This result is important in that the gain of the interconnection is independent of time delay, and furthermore, the delay elements in the matrix can be asymmetric. We will use this result now to design subsystem transfer functions G_i such that the overall connected system is stable by the small gain theorem. The small gain theorem can be used to obtain the constraints we will impose on G based on the interconnection gain

${\gamma }_2(G(j\omega ))\cdot (N-1)<1\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\forall \omega >0$

53

${\gamma }_2(G(j\omega ))<\frac{1}{N-1}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\forall \omega >0$

54

Furthermore, because G is diagonal, the maximum singular value of the matrix will belong to the subsystem G_i (s) with the maximum gain in frequency domain. So the control gains of all robots will be determined based on the robot transfer function with maximum gain (the worst case condition). This condition can be written as follows

${\gamma }_2(G(j\omega ))=\underset{\omega >0}{\text{sup}}\left\{\sqrt{{\lambda }_{\text{max}}\{G^{*}(j\omega )G(j\omega )\}}\right\}$

55

$=\underset{\omega >0}{\text{sup}}\hspace{0.17em}\left\{\underset{i}{\text{max}}\right\{\sqrt{G_i(-j\omega )G_i(j\omega )}\left\}\right\}$

56

$=\underset{\omega >0}{\text{sup}}\hspace{0.17em}\left\{\underset{i}{\text{max}}\right\{|G_i(j\omega )|\left\}\right\}$

57

Hence, we can guarantee the stability of the overall system by locally guaranteeing that

$|G_i(j\omega )|\hspace{0.17em}<\frac{1}{N-1}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\forall \omega \hspace{0.17em},i$

58

where G_i was defined for each robot as follows

$G_i(j\omega )=\frac{(C{1}_i(j\omega )-C{2}_i{Z}_i(j\omega ))}{m_i{(j\omega )}^2+(b_i+Z_i(j\omega ))j\omega +{\delta }_i{C}_p(j\omega )}$

59

We will next show that damping injection method will help us satisfy the above condition given that there is adequate damping, whereas the conventional four-channel-based architecture cannot. If the conventional four-channel-based controllers are selected as in equations (18) to (20), these gains are functions of frequency ω, the controller parameters, the impedance parameters and are independent of delay. For the sake of simplicity, we will assume that C _2i = C_f is chosen as 1 and then the gains become

$\left|G_i\right|=((K_p-k_{zi}{)}^2+2(K_p-k_{zi})(m_{zi}){\omega }^2+{(K_v-b_{zi})}^2{\omega }^2+{(m_{zi}{\omega }^2)}^2)/({({\delta }_i{K}_p+k_{zi})}^2-2({\delta }_i{K}_p+k_{zi})(m_{zi}+m_i){\omega }^2{+{(b_{zi}+{\delta }_i{K}_v)}^2{\omega }^2+{(m_{zi}+m_i)}^2{\omega }^4))}^{1/2}$

60

For the purposes of demonstration, we first examine the simplest case of free motion, where we neglect the effects of the environmental and operator impedances, and assuming m_i = 1 without loss of generality, we have the identical L ₂ gains for each robot

$\left|G_i\right|=\sqrt{\frac{{(K_p)}^2+{(K_v)}^2{\omega }^2}{{({\delta }_i{K}_p)}^2-2({\delta }_i{K}_p){\omega }^2+{({\delta }_i{K}_v)}^2{\omega }^2+{\omega }^4}}$

61

It can easily be seen that for ω ² < 2(δ_iK_p ) the gain will be greater than $\frac{1}{{\delta }_i}=\frac{1}{N-1}$ , rendering the system unstable in the sense of L ₂. On the other hand, if we use the local damping terms b_i in our proposed controller, these terms will be present in the denominator only, decreasing the gain

$\left|G_i\right|=\sqrt{\frac{{(K_p)}^2+{(K_v)}^2{\omega }^2}{{({\delta }_i{K}_p)}^2-2({\delta }_i{K}_p){\omega }^2+{({\delta }_i{K}_v+b_i)}^2{\omega }^2+{\omega }^4}}$

62

If they are chosen as

$b_i\ge \sqrt{2{\delta }_i{K}_p}$

63

This will force the system L ₂ gain to satisfy the small gain condition

$\begin{array}{l}\left|G_i\right|\hspace{0.17em}=\sqrt{\frac{{(K_p)}^2+{(K_v)}^2{\omega }^2}{{({\delta }_i{K}_p)}^2+{({\delta }_i{K}_v)}^2{\omega }^2+2({\delta }_i{K}_v{b}_i{)}^2{\omega }^2+{\omega }^4}}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}<\frac{1}{{\delta }_i}=\frac{1}{N-1}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\forall \omega >0\end{array}$

64

and guarantee the L ₂ stability of the multilateral system. Next, we proceed to show similar results for the general form of the multilateral teleoperation system where human and operator impedances are also included in the local gains. The L ₂ gain for each robot in contact with an environment/operator can be written explicitly and should satisfy the following small gain condition

$\left|G_i\right|\hspace{0.17em}=((K_p-k_{zi}{)}^2+2(K_p-k_{zi})(m_{zi}){\omega }^2+{(K_v-b_{zi})}^2{\omega }^2+{(m_{zi}{\omega }^2)}^2)/({({\delta }_i{K}_p+k_{zi})}^2-2({\delta }_i{K}_p+k_{zi})(m_{zi}+1){\omega }^2+{(b_i+b_{zi}+{\delta }_i{K}_v)}^2{\omega }^2{+{(m_{zi}+m_i)}^2{\omega }^4))}^{1/2}<\frac{1}{{\delta }_i}=\frac{1}{N-1}$

65

This can be rewritten as

${({\delta }_i{K}_p-{\delta }_i{k}_{zi})}^2+{(2{\delta }_i^2(K_p-k_{zi})m_{zi}+({\delta }_i{K}_v-{\delta }_i{b}_{zi}))}^2{\omega }^2+{({\delta }_i{m}_{zi}{\omega }^2)}^2<({\delta }_i{K}_p+k_{zi}{)}^2-2({\delta }_i{K}_p+k_{zi})(m_{zi}+m_i){\omega }^2+{(b_i+b_{zi}+{\delta }_i{K}_v)}^2{\omega }^2+{(m_{zi}+m_i)}^2{\omega }^4$

66

which would be satisfied if the controller gains are selected according to the impedance parameters

$K_p\ge \frac{k_{zi}({\delta }_i^2-1)}{2({\delta }_i^2+{\delta }_i)}$

67

$m_i\ge ({\delta }_i-1)m_{zi}$

68

$b_i\ge (2({\delta }_i{K}_p+k_{zi})(m_{zi}+m_i){+2{\delta }_i^2(K_p-k_{zi})m_{zi}+({\delta }_i^2-1)b_{zi}^2)}^{1/2}$

69

If these conditions are satisfied for each robot, the small gain condition is also satisfied, and therefore, closed-loop multilateral system is L ₂ stable. Equation (67) provides a bound for the position feedback gains, depending on environmental stiffness and the number of robots. Equation (68) shows that the ratio of the environmental masses with respect to the robot masses plays an important role in this analysis. This condition relating the robot masses and environmental/operator masses can be achieved by changing the nominal mass M_n = m_i in the disturbance observers, which also becomes a control design parameter. Equation (69) provides the minimum values that can be used as damping injection parameters. Again, it should be noted that we choose all the individual robot controller parameters identical and based on the worst case gain |G_i |, and this simplifies the control design. These analytical solutions provide us with an insight in designing the controller parameters with respect to the environment/operators and the number of robots, independent of time delay.

L ₂ stability analysis for dual-user teleoperation with dominance ratios

Some articles have also utilized methods for determining dominance among masters in dual-master teleoperation for various training purposes. ³⁴ These methods are modified versions of the multilateral teleoperation architecture presented here. In this section, we show that the proposed stability analysis and control design method can also be utilized to guarantee stability with dominance factors in dual master teleoperation. The method proposed in the study by Khademian and Hashtrudi-Zaad ³⁴ and used in some other multilateral teleoperation papers can also be written with our notation with time delay as

${\ddot{x}}^{ref}=-C_p(\mathcal{I}-{\mathcal{A}}_{\tau }^{\alpha })x-C_f(\mathcal{I}+{\mathcal{A}}_{\tau }^{\alpha })F$

70

where

${\mathcal{A}}_{\tau }^{\alpha }=\left[\begin{array}{ccc}0& (1-\alpha )e^{-s{T}_{12}}& \alpha e^{-s{T}_{13}}\\ \alpha e^{-s{T}_{21}}& 0& (1-\alpha )e^{-s{T}_{23}}\\ \alpha e^{-s{T}_{31}}& (1-\alpha )e^{-s{T}_{32}}& 0\end{array}\right]$

71

and α term here is the dominance factor which determines the command ratio of one master over the other and can be set between 0 ≤ α ≤ 1. For instance with the above law, when α = 0, the second robot becomes the sole master and the first robot can only follow the motion of the master which bilaterally operates the third robot. When α = 1, the first robot becomes the sole master and the second robot follows the master; when α = 0.5, the system becomes equivalent to the multilateral teleoperation without dominance. For values of alpha between 0, 0.5 and 1, the dominance is skewed towards one side. Stability of this system can be proven in a similar manner to the previous section. First, we introduce the damping injection–based control law

${\ddot{x}}^{ref}=-C_p(\mathcal{I}-{\mathcal{A}}_{\tau }^{\alpha })x-B\dot{x}-C_f(\mathcal{I}+{\mathcal{A}}_{\tau }^{\alpha })F$

72

In the stability analysis, the main difference is in the G_i and H functions. New system transfer functions become

$u_i=e_i-{\displaystyle \sum _{j=1}^N}{H}_{ij}{y}_j$

73

$y_i=G_i(u_i)$

74

$G_i=\frac{(C{1}_i(s)-C{2}_i{Z}_{i}s)}{m_i{s}^2+(b_i+Z_i(s))s+C_p(s)}$

75

$H={\mathcal{A}}_{\tau }^{\alpha }$

76

The difference in G_i is the presence of 1 instead of δ_i term, and ${\mathcal{A}}_{\tau }^{\alpha }$ is used instead of ${\mathcal{A}}_{\tau }$ for H. The L ₂ stability condition is the same as before and can be written as

${\gamma }_2(G(j\omega ))\cdot {\gamma }_2({\mathcal{A}}_{\tau }^{\alpha }(j\omega ))<1,\omega >0$

77

The gain ${\gamma }_2({\mathcal{A}}_{\tau }^{\alpha }(j\omega ))$ can be obtained in a similar vein. It can be shown that ${\mathcal{A}}_{\tau }^{\alpha *}{\mathcal{A}}_{\tau }^{\alpha }$ is also Hermitian. The eigenvalues are therefore real and again utilizing the Gershgorin circle theorem we have that

$1\le \lambda ({\mathcal{A}}_{\tau }^{\alpha *}{\mathcal{A}}_{\tau }^{\alpha })\le 2$

78

H_∞ gain of ${\mathcal{A}}_{\tau }^{\alpha }(j\omega ))$ is

$\underset{\omega >0}{\text{sup}}\left\{{\sigma }_{\text{max}}\right({\mathcal{A}}_{\tau }^{\alpha }(j\omega )\left)\right\}\le \sqrt{2}$

79

Then from small gain theorem, local robot gain condition is

$\left|G_i\right|<\frac{1}{\sqrt{2}}$

80

To satisfy this condition, the controller parameters can be chosen as

$K_p>0$

81

$b_i\ge (2(K_p+k_{zi})(m_{zi}+m_i){+2\sqrt{2}(K_p-k_{zi})m_{zi}+(\sqrt{2}-1)b_{zi}^2)}^{1/2}$

82

which guarantees the L ₂ stability of the dominance ratio–based multilateral system. It should be noted that these parameters are less involved than the previous results, and actually, the stability depends only on the damping. This is due to the lack of δ_i in the local feedback.

Robust stability analysis

We have shown the conditions necessary for the delay independent L ₂ stability of the closed-loop multilateral teleoperation system. If these conditions are satisfied, internal stability is guaranteed for the nominal teleoperation system with the assumed controller, environment and system parameters. However, in reality, it is known that all models can suffer from inaccuracies. In robust control, the idea is to describe the inaccuracies in the model as uncertainties that perturb a nominal model and to design a controller for the whole class of systems. In our system, the main source of uncertainty is the master and environment impedance parameters. The robot models in our case are rather simple and the model errors can be considered negligible, and the disturbance observers we employ force the system dynamics to that of the nominal model. This can even be extended to more complicated multi-DOF robots. ^30,31 Therefore, in this article, we will assume uncertainties in the environmental and master parameters. A similar analysis has been performed in the study by Lin and Namerikawa ³⁵ for bilateral teleoperation and in the study by Khademian and Hashtrudi-Zaad ¹⁵ for multilateral teleoperation. However, Khademian and Hashtrudi-Zaad ¹⁵ have performed such an analysis for the multilateral teleoperation system without time delay. As we will demonstrate, it is possible to accommodate delays in robust stability analysis of the multilateral teleoperation system. For this purpose, we are going to describe the uncertainty in our system with parametric uncertainty representation and analyse the robustness properties of the controller architecture we propose. We are going to pull out the uncertainty in the multilateral teleoperation system into diagonal structured form which is in feedback configuration with the nominal system and by using the structured singular value μ and the generalized small gain theorem, we will be able to show the robust stability of the proposed multilateral system with time delay. We will first start this analysis by introducing the parametric uncertainty range in our system.

Here, we will be assuming that any of our robots can be in contact with an environment or a master operator, and the bounds for impedance parameters with significant variation percentages of 100% are given in Table 1. The selection of a wide range of master–environment impedances helps us cover a wide range of experiment conditions, including free motion, contact motion, single–dual master teleoperation, master switching, variations in impedances and so on. In our analysis, we make use of notations of Skogestad ³⁶ and Zhou. ³⁷ The uncertain impedances can be written with parametric uncertainty as

${\tilde{k}}_{zi}=k_{zi}(1+{\delta }_{kzi}),\left|{\delta }_{kzi}\right|\hspace{0.17em}\le 1$

83

${\tilde{b}}_{zi}=b_{zi}(1+{\delta }_{bzi}),\left|{\delta }_{bzi}\right|\hspace{0.17em}\le 1$

84

${\tilde{m}}_{zi}=m_{zi}(1+{\delta }_{mzi}),\left|{\delta }_{mzi}\right|\hspace{0.17em}\le 1$

85

Table 1.

Uncertain parameters.

Parameter	Nominal value	Variation range
m z1, m z2, m z3	0.5 kg	0–1 kg
b z1, b z2, b z3	5 N s/m	0–10 N s/m
k z3, k z3, k z3	5000 N/m	0–10,000 N/m

where δ_mzi , δ_kzi , δ_bzi are the scalar parametric uncertainty variations in the environmental impedance parameters with respect to the nominal parameters and ${\tilde{m}}_{zi},{\tilde{b}}_{zi},{\tilde{k}}_{zi}$ are the uncertain family of impedances. With these definitions, the generalized plant form of the robot i in the multilateral control system can be obtained from G_i and be drawn as in Figure 5. In the figure, $\overline{m_i}$ is defined as

$\overline{m_i}=\left[\begin{array}{cc}\frac{1}{m_i+m_{zi}}& \frac{-m_{zi}}{m_i+m_{zi}}\frac{1}{m_i+m_{zi}}\\ 1& \frac{-m_{zi}}{m_i+m_{zi}}\end{array}\right]$

86

Figure 5.

Pulling out the uncertainties and the generalized plant representation of G_i .

Then, $\frac{1}{m_i+m_{zi}+{\delta }_{mzi}{m}_{zi}}$ is represented as a lower LFT of $\overline{m_i}$ with δ_mzi

$\frac{1}{m_i+m_{zi}+{\delta }_{mzi}{m}_{zi}}=F_{\ell }(\overline{m_i},{\delta }_{mzi})={\overline{m_i}}_{11}+{\overline{m_i}}_{12}{\delta }_{mzi}{(\mathcal{I}-{\overline{m_i}}_{22}{\delta }_{mzi})}^{-1}{\overline{m_i}}_{21}$

87

The uncertainty input and output relationships can be written as

$\underset{u_i\Delta }{\underbrace{\left[\begin{array}{c}{u}_i{\Delta }_1\\ u_i{\Delta }_2\\ u_i{\Delta }_3\\ u_i{\Delta }_4\\ u_i{\Delta }_5\\ u_i{\Delta }_6\end{array}\right]}}=\underset{{\Delta }_i}{\underbrace{\left[\begin{array}{cccccc}{\delta }_{mzi}& 0& 0& 0& 0& 0\\ 0& {\delta }_{bzi}& 0& 0& 0& 0\\ 0& 0& {\delta }_{kzi}& 0& 0& 0\\ 0& 0& 0& {\delta }_{kzi}& 0& 0\\ 0& 0& 0& 0& {\delta }_{bzi}& 0\\ 0& 0& 0& 0& 0& {\delta }_{mzi}\end{array}\right]}}\underset{y_i\Delta }{\underbrace{\left[\begin{array}{c}{y}_i{\Delta }_1\\ y_i{\Delta }_2\\ y_i{\Delta }_3\\ y_i{\Delta }_4\\ y_i{\Delta }_5\\ y_i{\Delta }_6\end{array}\right]}}$

88

And the generalized plant representation P_i of G_i can be written as in equation (89). Here, u_Pi and y_Pi are the inputs and outputs to the generalized plant, respectively, and the complete generalized plant P for the multilateral subsystem G can then be obtained as

$\underset{y_{Pi}}{\underbrace{\left[\begin{array}{c}{y}_i{\Delta }_1\\ y_i{\Delta }_2\\ y_i{\Delta }_3\\ y_i{\Delta }_4\\ y_i{\Delta }_5\\ y_i{\Delta }_6\\ {\dot{x}}_i\\ {\ddot{x}}_i\\ y_i\end{array}\right]}}=\underset{P_i}{\underbrace{\left[\begin{array}{cccccccccc}\frac{-m_{zi}}{m_i+m_{zi}}& -1& -1& 0& 0& 0& -K_p-k_{zi}& -b_i-b_{zi}& 1& -1\\ 0& 0& 0& 0& 0& 0& 0& b_{zi}& 0& 0\\ 0& 0& 0& 0& 0& 0& k& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& k& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& b_{zi}& 0& 0\\ \frac{-m_{zi}^2}{{(m_i+m_{zi})}^2}& \frac{-m_{zi}}{m_i+m_{zi}}& \frac{-m_{zi}}{m_i+m_{zi}}& 0& 0& 0& \frac{(-K_p-k_{zi})m_{zi}}{m_i+m_{zi}}& \frac{(-b_i-b_{zi)m_{zi}}}{m_i+m_{zi}}& \frac{m_{zi}}{m_i+m_{zi}}& \frac{-m_{zi}}{m_i+m_{zi}}\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ \frac{-m_{zi}}{{(m_i+m_{zi})}^2}& \frac{-1}{m_i+m_{zi}}& \frac{-1}{m_i+m_{zi}}& 0& 0& 0& \frac{-K_p-k_{zi}}{m_i+m_{zi}}& \frac{-b_i-b_{zi}}{m_i+m_{zi}}& \frac{1}{m_i+m_{zi}}& \frac{-1}{m_i+m_{zi}}\\ \frac{m_{zi}^2}{{(m_i+m_{zi})}^2}& \frac{m_{zi}}{m_i+m_{zi}}& \frac{m_{zi}}{m_i+m_{zi}}& -1& -1& -1& \frac{((K_p-k_{zi})(m_i+m_{zi})+(K_p+K_{zi})m_{zi})}{m_i+m_{zi}}& \frac{((-b_{zi})(m_i+m_{zi})+(b_i+b_{zi})m_{zi})}{m_i+m_{zi}}& \frac{-m_{zi}}{m_i+m_{zi}}& \frac{m_{zi}}{m_i+m_{zi}}\end{array}\right]}}\underset{u_{Pi}}{\underbrace{\left[\begin{array}{c}{u}_i{\Delta }_1\\ u_i{\Delta }_2\\ u_i{\Delta }_3\\ u_i{\Delta }_4\\ u_i{\Delta }_5\\ u_i{\Delta }_6\\ x_i\\ {\dot{x}}_i\\ e_i\\ v_i\end{array}\right]}}$

89

$\underset{Y_P}{\underbrace{\left[\begin{array}{c}{y}_{P1}\\ ⋮\\ y_{P(N-1)}\\ y_{PN}\end{array}\right]}}=\underset{P}{\underbrace{\left[\begin{array}{cccc}{P}_1& 0& \cdots & 0\\ 0& \ddots & 0& 0\\ ⋮& \cdots & P_{N-1}& 0\\ 0& \cdots & 0& P_N\end{array}\right]}}\underset{U_P}{\underbrace{\left[\begin{array}{c}{u}_{P1}\\ ⋮\\ u_{P(N-1)}\\ u_{PN}\end{array}\right]}}$

90

The uncertainty for the multilateral system then becomes

$\underset{u\Delta }{\underbrace{\left[\begin{array}{c}{u}_1\Delta \\ ⋮\\ u_{N-1}\Delta \\ u_N\Delta \end{array}\right]}}=\underset{\Delta }{\underbrace{\left[\begin{array}{cccc}{\Delta }_1& 0& \cdots & 0\\ 0& \ddots & 0& 0\\ ⋮& \cdots & {\Delta }_{N-1}& 0\\ 0& \cdots & 0& {\Delta }_N\end{array}\right]}}\underset{y\Delta }{\underbrace{\left[\begin{array}{c}{y}_1\Delta \\ ⋮\\ y_{N-1}\Delta \\ y_N\Delta \end{array}\right]}}$

91

Once the generalized plant is obtained, the system can be transformed into NΔ structure by taking the feedback structure into account. Here, in our system, the feedback gain K is taken to be the interconnection matrix H which includes time delay. The output of H in equation (32) is denoted as v, such that v = HY. N, which represents the nominal transfer function of the closed loop system as in Figure 6, can be found with the lower LFT

$N=T_{EZ}=F_{\ell }(P,K)=P_{11}+P_{12}K{(\mathcal{I}-P_{22}K)}^{-1}{P}_{21}$

92

$=F_{\ell }(P,H)=P_{11}+P_{12}H{(\mathcal{I}-P_{22}H)}^{-1}{P}_{21}$

93

Figure 6.

Generalized plant representation of the multilateral teleoperation system and the transformation to MΔ structure.

In Figure 6, Z is the performance output of the control system and can be chosen as Z = Y where Y = [y ₁ y ₂⋯y_N ]^T. In H, a second-order Padé approximation of time delay can be utilized for inversion. The first block diagonal element of N, N ₁₁ is particularly important as it determines the MΔ representation in Figure 6 with

$M=N_{11}$

94

Determining the robust stability of the multilateral teleoperation system now becomes equivalent to the internal stability of MΔ configuration, and we can write the generalized small gain condition as

$\mu (M(j\omega ){\sigma }_{\text{max}}(\Delta (j\omega ))<1,\forall \omega$

95

where μ is the structured singular value defined as

$\mu (M)\triangleq \frac{1}{\text{min}({\kappa }_m)|det(I-{\kappa }_{m}M\Delta )=0,\text{for}\hspace{0.17em}\text{structured}\hspace{0.17em}\Delta ,\sigma (\Delta )\le 1}$

96

which is the inverse of the robust stability margin, such that $\mu (M)=\frac{1}{{\kappa }_m}$ . By definition, σ _max(Δ(jω)) ≤ 1 so the robust stability condition becomes

$\mu (M)\triangleq \frac{1}{{\kappa }_m}<1,\forall \omega$

97

By its nature, robust stability analysis is performed numerically. MATLAB’s Robust Control Toolbox has been utilized to check the robust stability with μ/κm analysis. Robust stability margins of various damping injection–based four-channel controllers, with respect to the large uncertainty set in Table 1, are found to be as in Table 2. It can be seen that even with such high levels of uncertainty, by increasing the local damping term and the position gain together in the proposed architecture, it is possible to increase robustness. In Table 2, stability margins greater than 1 imply robust stability over the full uncertainty range. Thus, the robust control design problem reduces to finding a gain and damping term for the assumed uncertainty. In our analyses, we have assumed identical delays in the first six cases and asymmetric delays in the final case. The μ analysis results show that our proposed control architecture is capable of achieving robust stability under parametric uncertainty and asymmetric time delays, without resorting to high-order H_∞ controllers. The proposed analysis method provides a foundation for more involved analyses as well. It is also possible to describe delay as uncertainty, including varying time delays. Nonlinear or unspecified dynamics in the environment can also be taken into account with various Δ assumptions. These are beyond the scope of this article but can be achieved with some additions to the current analysis, for example by utilizing IQC methods.

Table 2.

Robust stability margins.^a

Controller and system parameters	Robust stability margin inf ω>0(min((κm (M(jω)))
Mn = 0.2, Kp = 2500	0.2777
Kv = 0, bi = 100, Tij = 0.5
Mn = 0.2, Kp = 2500	0.6439
Kv = 0, bi = 200, Tij = 0.5
Mn = 0.2, Kp = 3000	0.8249
Kv = 0, bi = 200, Tij = 0.5
Mn = 0.2, Kp = 2500	0.9965
Kv = 0, bi = 300, Tij = 0.5
Mn = 0.2, Kp = 2500	1.000
Kv = 0, bi = 400, Tij = 0.5
Mn = 0.2, Kp = 3000	1.133
Kv = 0, bi = 300, Tij = 0.5
Mn = 0.2, Kp = 3000	1.0222
Kv = 0, bi = 300
T 12 = 1, T 21 = 0.5, T 31 = 0.6,
T 13 = 0.7, T 23 = 0.4, T 32 = 0.8

^aValues greater than 1 imply robust stability.

Experiment results and discussions

For experiments, we present results for a trilateral teleoperation system. We use linear motors as 1 DOF master and slave robots. This enables us to decouple the teleoperation problem from local dynamics and focus on the performance of the teleoperation system. The robots are connected to each other through a single computer, and delay among them is artificially created in the computer. Trilateral teleoperation is a subset of multilateral teleoperation and suffices to confirm our results, as the control law and analysis are scalable and the number of robots can easily be increased. The experiment setup can be seen in Figure 7. In the experiments, to compare the results with different control laws and delays, a similar motion pattern is repeated. Two master robots are simultaneously moved by two masters and a slave robot moves freely and then contacts a rigid wooden surface applying some force. Figure 7 is enumerated for the purposes of clarification, numbers 1–3 denote the linear motors used as the robots. Motors are STA1116 by Dunkermotoren which are equipped with hall sensors/encoders for position measurement, the forcers are fixed to table and the operators move the shafts. Numbers 1 and 2 denote the master robots and number 3 denotes the slave robot. Number 4 denotes the environment that the slave robot contacts, which is solid wood. Number 5 denotes the communication interface between the motors and the computer which is composed of two NI 6321 DAQ cards. Finally, number 6 denotes the desktop computer which runs the real-time control algorithm. Reaction force observers are used instead of force sensors in these experiments. An experiment with the conventional four-channel-based controller was done with identical controller parameters for all robots chosen as follows: M_n = 0.2kg, K_p = 900N/m, K_v = 60Ns/m, g _dis = g _reac = 500. Figure 8 shows that the system becomes unstable with this controller when symmetric time delays of 0.01 s between robots exist. It can be seen from the figure that at the start of the motion, for a very brief period of time, the system is stable and the robots track each other, but when the force input is slightly increased, the system becomes unstable and very high-frequency oscillations start, and then system is shut down for safety reasons. To compare the proposed method with other methods in the literature, passive and absolutely stable controllers as mentioned in the study by Li et al. ²⁷ in the form of two-channel-based PPP (position–position–position) architecture have been implemented. The controller parameters guaranteeing passivity in a trilateral system such as ours have been given in the study by Li et al. ²⁷ according to our notation as

$C{1}_i=C_i=K_{pi}+K_{vi}$

98

$C{2}_i=0$

99

$K_{p1}=K_{p2}=K_{p3}$

100

$K_{v1}=K_{v2}=K_{v3}$

101

Figure 7.

Multilateral teleoperation system.

Figure 8.

Unstable system with 10-ms time delay between robots: (a) position response and (b) force response.

which means that the robots send only their position information to one another, and the controllers for position control are PD with identical gains for each robot. So in the experiments, PD gains were chosen as K_p = 2500, K_v = 100 which are chosen to be the same as the gains of the proposed controller. Note that, unlike the proposed method, this control law does not have any local damping injection and force transmission among the robots. Also, since there is no force transmission, disturbance observers are not used in the control loop and are implemented just to monitor force responses. Experiment results with 0.2-s time delay between each robot are given in Figure 9. It can be seen that the system becomes stable, and position tracking is achieved during free motion, but during contact (between seconds 9 and 16, 29 and 35), a considerable steady-state error in position responses between the slave and the masters can be seen due to the fact that force reflection can only be realized by creating a steady-state error. This is a systemic problem of two-channel/passivity controllers, and although it is possible to decrease the error, it is impossible to eliminate it completely by increasing/decreasing the PD gains. These particular gain values were chosen to enable a comparison with the proposed method. Increasing the gains can decrease the errors but in return make the system much heavier and slower in response.

Figure 9.

Performance of passive two-channel control system with 200-ms time delay between robots: (a) position response and (b) force response.

Another stability analysis method in the literature is absolute stability method. We have mentioned before that the absolute stability conditions are more relaxed compared to passivity conditions. For instance, according to Li et al., ²⁷ the criteria for a control system to be absolutely stable are given as

$C{1}_i=C_i=K_{pi}+K_{vi}$

102

$C{2}_i=0$

103

$\frac{K_{p1}}{K_{v1}}=\frac{K_{p2}}{K_{v2}}$

104

$5-2\sqrt{6}\le \frac{K_{p1}}{K_{v1}}\frac{K_{v3}}{K_{p3}}\le 5+2\sqrt{6}$

105

These conditions again require a two-channel-based PPP controller and imply that, unlike the passive system, gains for each individual robot can be different from other robots’ gains, but the ratio between the proportional and derivative gains should still be related and constrained. It should be noted that any passive system that satisfies the criteria (70–73) is also absolutely stable, but the reverse is not necessarily true. This makes absolute stability a more broad set and provides more options to the control system designer in choosing the control parameters. Absolutely stable but not passive controllers were also tested with experiments by choosing K _p1 = 2500, K _v1 = 100, K _p2 = 2500, K _v2 = 100, K _p3 = 7500, K _v3 = 300. The experiment results with 0.2-s time delay between each robot are given in Figure 10. It can be seen that the system is stable, tracking in free motion is good, but during contact (between seconds 11 and 15, 23 and 32), again a large steady-state error exists although it is smaller in magnitude compared to the passive system. Furthermore, during contact, the sum of the master forces is roughly equal to one-third of the slave force due to K _p3 being 3 times K _p1 and K _p2. In effect, this system becomes a force-scaled teleoperation architecture, but its application is beyond the scope of this article.

Figure 10.

Performance of absolutely stable two-channel control system with 200-ms time delay between robots: (a) position response and (b) force response.

These results show that a controller based on the more relaxed L ₂ stability criteria which enables us to choose gains independently for each robot/environment and utilize the four-channel architecture could really be an advantage and actually provide improvements in performance/transparency. Indeed, experiments show that when the proposed controller with local damping is utilized, the system is stable and position tracking and force reflection are successfully achieved as can be seen from Figure 11. To choose the controller parameters that make the system stable based on the derived criteria, some assumptions about the environmental parameters have been made. Based on the study by Peer and Buss, ³⁸ maximum environment stiffness has been taken as 10,000 N/m and, for the wood surface, the damping and mass are assumed to be zero, because the wood is fixed and cannot move, and based on the study by Niculescu ³⁹ and Speich, ⁴⁰ the human operator impedances have been taken as 48/s + 4.5 + 0.2s. Niculescu et al. ³⁹ neglect operator hand mass but we have used a value equal to motor mass for the sake of robustness. According to the assumed impedance values, controller parameters for all robots have been chosen identical to the maximum one and are M_n = 0.2kg, K_p = 2500N/m, K_v = 0Ns/m, b_i = 100Ns/m, g _dis = g _reac = 500. If K_p is chosen equal to the RHS of (67) and M_n = m_i , the maximum damping of 100 will be necessary in free motion, so choosing the controller parameters is rather straightforward. Note that this controller is the least robust controller in Table 1 and has a robust stability margin of 0.2777. This means, about 27.7% of the specified uncertainty levels in Table 1 can cause instability in the system. However, in our experiments, we have not witnessed such instabilities, implying that our uncertainty assumption may be too conservative (e.g. we are assuming nominal stiffness of 5000 N/m for all robots, which is actually well beyond the real conditions). Figure 11 shows that the system is stable during free motion and contact (between seconds 12 and 18). In free motion, operational forces due to damping are felt by the masters and this is what makes the system slower. This actually is how the system is stabilized; the controller makes the system slow enough such that there is almost no lag between the motions of the master and slave systems. Meanwhile, with the two-channel controllers, the lag between the system position responses is visible. More importantly, it can also be seen that the four-channel-based system moves faster than the two-channel-based systems with the same amount of delay. Compared with the two-channel controllers, the proposed control system takes much less time to complete a similar motion. Since our proposed stabilization method can be applied to four-channel controllers, good performance under time delay can be achieved. In free motion, the robots are perfectly synchronized and when the slave contacts an environment, a small overshoot occurs in the masters but the master robots quickly track the correct position of the slave robot. Furthermore, there is no steady-state error during contact, which is a common problem of the passivity-based approaches to bilateral and multilateral teleoperation systems, and this can be seen from Figures 9 and 10. It can also be seen that the force reflection goal is also achieved during contact and the sum of the forces is equal to zero. The amount of force on a single master robot depends on how much force input the operator applies, but the sum of the forces is always controlled to be zero.

Figure 11.

Performance of the proposed four-channel damping injection control system with 200-ms time delay between each robot: (a) position response and (b) force response.

Figure 12 shows the results for the same controller parameters with a delay of 500 ms in both directions between each robot. This confirms that the stability analysis and controller design are indeed delay independent. Performance is not as good as the 200 ms case, and the system is even slower; however, stability is preserved and position tracking in free motion and during contact (10–18 s) is achieved again without drift. Force reflection is also successfully achieved during contact as desired.

Figure 12.

Performance of the proposed four-channel damping injection control system with 500-ms time delay between each robot: (a) position response and (b) force response.

Finally, Figure 13 confirms another contribution of the article. The same controller parameters used before can be used to accommodate asymmetric time delays in the system. In this experiment, delays between the robots are T ₁₂ = 200 ms, T ₁₃ = 500 ms, T ₂₁ = 300 ms, T ₃₁ = 100 ms, T ₂₃ = 400 ms, T ₃₂ = 600 ms. As can be seen from the figure, the system achieves position tracking and force reflection during contact (8–18 s) successfully. In contrast to identical delay case, a small amount of lag can be seen between the robot responses, caused by the differences in delay.

Figure 13.

Performance of the proposed four-channel damping injection control system with asymmetric time delays between each robot ranging from 100 ms to 600 ms: (a) position response and (b) force response.

Conclusions

This article presented a stability and robustness analysis as well as a four-channel-based control synthesis method for multilateral teleoperation systems under symmetric/asymmetric time delays. The analysis and design method developed in this article can also be employed for multilateral systems with different control architectures not presented in this article. The method can be further developed to include varying time delays. Also, it is especially suitable to designing optimal/suboptimal robust controllers and nonlinear environmental and operator parameters can be inserted in this methodology using uncertainty formulations.