# Stable Bilateral Teleoperation Control Method for Biped Robots with Time-Varying Delays

Authors: Viviana Moya; Emanuel Slawiñski; Vicente Mut; Danilo Chávez; Bernardo Wagner
Journal: Journal of Robotics
 (2023)

---

## Abstract

This document proposes a control scheme applied to delayed bilateral teleoperation of the forward and turn speed of a biped robot against asymmetric and time-varying delays. This biped robot is modeled as a hybrid dynamic system because it behaves as a continuous system when the leg moves forward and discrete when the foot touches the ground generating an impulsive response. It is proposed to vary online the damping according to the time delay present in the communication channel, and the walking cycle time using an optimization criterion, to decrease the teleoperation system errors. To accomplish this, a three-phase cascade calibration process is used, and their benefits are evidenced in a comparative simulation study. The first phase is an offline calibration of the inverse dynamic compensation and also the parameters of the bilateral controller. The second phase guarantees the bilateral coordination of the delayed teleoperation system, using the Lyapunov–Krasovskii stability theory, by changing the leader damping and the equivalent follower damping together. The third phase assures a stable walk of the hybrid dynamics by controlling the walking cycle time and the real damping to move the eigenvalues of the Poincaré map, numerically computed, to stable limit cycles and link this result with an equivalent continuous system to join both phases. Additionally, a fictitious force was implemented to detect and avoid possible collisions with obstacles. Finally, an intercontinental teleoperation experiment of an NAO robot via the Internet including force and visual feedback is shown.


---

## Body

## 1. Introduction

The teleoperation of robots allows the capacities and skills of a human operator to be extended and transferred to remote working environments. In these systems, a human operator situated at a local site sends control commands from a leader device to a robot named a follower located in a remote environment, such as humanoid robots, which are used to perform dangerous work applied to mining, nuclear plant inspection, explosives deactivation, defense, rehabilitation, healthcare, and among others. The leader and follower are connected via a communication channel that adds varying time delays, which cause instability or generally poor performance, resulting in higher coordination errors during the execution of a task as well as inadequate transparency [1].Stability is one of the most important properties of bipedal locomotion. “Stable” walking can be defined as any walking that does not result in a fall. This implies that the definition of stability should relate to the set of all states that a bipedal walker can experience and still avoid falling [2]. However, this definition is difficult to handle in practice, due to the size of this set and the absence of systematic tools to synthesize controllers; therefore, another more heuristic but the more focused definition of “stable” walking is that of persisting also in the presence of perturbations, i.e., the ability of a person, animal or robot to continue the motion that was planned in the face of the perturbation. Moreover, that these recovery actions sometimes completely modify the movement or even the type of movement, e.g., recovering from a walking perturbation by coming to a complete stop. Another way to fundamentally study stability is to look at the entire motion or at least one whole walking step or double step at a time and evaluate how small perturbations would affect this particular motion. This approach follows Lyapunov’s mathematical stability theory and has been frequently used in robotics in the field of walking (passive [3]-dynamic [4]). The concept of Lyapunov stability can be applied especially to motions, also called limit cycles, and is therefore very interesting for the study of walking and running. In this case, the propagation of perturbations over a walking cycle is typically one step. The study uses the Poincaré Map that maps the states at the beginning and at the end of the cycle. A periodic solution of the dynamic equations corresponds to a fixed point of the Poincaré map, as shown in Figure 1.Figure 1
Poincaré maps.To realize the combination of bilateral teleoperation with humanoid robots, the control problem has become more challenging in the following aspects: the simultaneous stability of walking, the error between the commands (operator-robot), the motion of the follower taking into account time-varying delays, the analysis of how many degrees of freedom (DOF) should include the haptic device, how the force feedback should be established, and the hybrid dynamics of humanoid robots with nonlinearities and uncertainties. The autonomous control of humanoid robots has been investigated for decades, which mainly includes the methods based on ZMP (zero moment point) and approaches based on dynamic walking.On the other hand, most papers addressed to the teleoperation of biped or humanoid robots are based on using motion capture devices in order to map the movements of a human operator to the robot, as in reference [5], or employing two traditional joysticks [6]. The main disadvantage of both articles is that the user does not receive any kind of tactile feedback. For the case of bilateral teleoperation, which includes force feedback, there are proposals that use haptic devices with low DOF with respect to the follower DOF and others that use a number of DOF in the leader as similar as possible to the DOF of the follower. In general, the last research is focused on the leader mechanical custom design, as in [7], where body teleoperation is used to close the loop by applying force feedback to the operator’s waist, and [8], whose leader is called TABLIS full-body exoskeleton cockpit that allows back feed to the human operator, a sum of many forces, although the inclusion of more tactile information not necessarily improves the performance based on metrics for task-oriented human-robot interaction (HRI). In the strategies using low DOF haptic devices, we remark [9] where the force feedback is based on the zero moment point (ZMP) method. Such force is computed depending on the stability margin, relative to the support polygon limits. Getting close to the edges of the feet will decrease stability and the system is more likely to lose balance. Quantifying this margin represents the basic principle involved in force feedback generation. On the other hand [10, 11], we consider dynamic walking with a force cue based on the coordination error but both use equivalents followers and controllers that are not formally defined, lacking thus the explicit link between the stability of the limit cycle using Poincaré map about the walking and the stability of the errors of delayed teleoperation system.The main weakness found in the state of the art is that the stability of walking and teleoperation errors considering hybrid dynamic modeling is not analyzed objectively, therefore is partially solved as the controller parameters should be online changed to hold stability. A comparative table of the contribution of different articles and this work is presented in Figure2, where it is remarked that the pink boxes are the approaches that we are going to deal with in the following article.Figure 2
Comparative analysis.This research is focused on the bilateral teleoperation of biped robots, where it is not clear whether, as in the case of teleoperated manipulators and mobile robots, if the stability of teleoperated bipedal robots can be ensured by increasing only the damping of leader and follower. In this document, a control scheme aimed at the bilateral teleoperation of the forward and turning speeds of a biped robot represented by a hybrid dynamic model in the presence of variable time delays is presented. The main contribution is the following statements:(1)
A control scheme designed to assure simultaneous stability of walking and teleoperation errors considering the hybrid dynamics of the biped robot as well as the presence of time-varying delays(2)
Explicit definition of equivalent follower and equivalent follower controller, which adds a follower equivalent damping depending on the time delay and allows to join the Lyapunov–Krasovskii analysis applied delayed nonlinear systems with the stability of a Poincaré map about walking of biped represented with a hybrid dynamic model(3)
Use of optimization criterion because there are multiple combinations of real follower damping and walking cycle time that hold the conditions of proposed equivalent damping(4)
Calibration process of controller parameters involving how the damping and walking cycle time should be changed onlineThis article is organized as follows: Section2 presents the leader and follower models and the assumptions used in this document. Section 3 describes the controller for bilateral teleoperation control and its stability analysis also is demonstrated. Next, Section 4 presents the experimental results achieved where a human operator located in San Juan-Argentina drives a humanoid NAO robot located in Hannover, Germany. Finally, the conclusions are presented in Section 5.


## 2. Models

### 2.1. Leader Model

The leader device is represented in Cartesian coordinates as follows:(1)Mmxmx..m+Cmxm,x.mx.m+gmxm=fm+fh,where xm=xmvxmδT∈Rnm×1 and x.m are the position and velocity of the leader, Mmxm∈Rnm×nm is the inertia matrix, Cmxm,x.m∈Rnm×nm is the matrix representing the centripetal and coriolis forces, gmxm∈Rnm is the gravity vector, fh is the force of the human operator, and fm is the control action applied to the leader.


### 2.2. Follower Model

A bipedal robot can be modeled in a simplified way in a sagittal plane for straight-line walking references. The turn is performed using the hip yaw angles of the 3D bipedal robot, taking into account the invariance property against yaw angle rotation [12]. The presence of impact of the legs on the environment during the walking cycle is represented with hybrid models. These systems, also called systems with impulsive effects, are used to represent a bipedal robot [13, 14], that includes two behaviors: continuous when the leg moves forward, and discrete when the foot touches the ground generating impulsive responses [15].


#### 2.2.1. Continuous Dynamic

The continuous dynamics of the bipedal robot between successive impacts presents a total number of physical actuators equal tonb. The state of the bipedal robot is characterized by x=qq.T, where q∈ℝnb is the position and q.∈ℝnb is the velocity. Using the Euler–Lagrange equations, the state between impacts is represented by x.=fx+gxu, as follows:(2)fx=q.−M−1qCq,q.q.+Gq,gx=0M−1qBq,where Mq∈ℝnb×nb represents the inertial matrix, Cq,q.∈ℝnb×nb is the matrix representing the centripetal and Coriolis forces, G∈ℝnb is the gravity vector, B∈ℝnb×nb is the input matrix, and u∈ℝnb is the matrix of control actions. In addition, the angles are defined as follows: q=qp_nshqy_nshqp_nskqp_nsaqp_shqy_shqp_skqp_saT, where qp_nsh is pitch non-stance hip, qp_nsk is pitch non-stance knee, qp_nsa is pitch non-stance ankle, qp_sh is pitch stance hip, qp_sk is pitch stance knee, qp_sa is pitch stance ankle, qy_nsh is yaw non-stance hip, and qy_sh is yaw stance hip, as shown in Figure 3.Figure 3
Joint representation.


#### 2.2.2. Domain and Guard

The domainDH specifies the allowed system settings based on the height of the non-stance foot l, measured from a rigid and plain floor, and indicates that it must be above the ground. Additionally, the guard SH represents the impact surface. The domain and guard are defined in reference [13] as follows:(3)DH=x∈R2nb:lq≥0,SH=x∈R2nb:lq=0and∂lq∂qq.<0,where x∈R2nb is the total number of control states given by q and q. and ∂lq/∂q is the partial derivative of l with respect to q.


#### 2.2.3. Discrete Dynamic

The discrete dynamics of the robot determine how the robot’s speeds changes when the foot hits the ground, while at the same time switching legs. In particular, the reset mapΔ:SH⟶DH is given by: Δq which represents the switching between the legs when going from “stance” to “no stance” and vice versa after the end of the double support phase, whereas Δq. determines the variation in the speed of each joint due to impact, and therefore, Δ is the vector defined in [13] that links both Δq and Δq.. It is important to mention that a stable continuous dynamic between impacts does not guarantee a stable hybrid system. Hence, it is assumed that there is a dynamic hybrid system represented as follows [14, 15]:(4)W=x.=fx+gxux,ifx∈DHandSH,x+=Δx−,ifx−∈SH.


### 2.3. Assumptions, Properties, and Nomenclature

This document uses the assumptions and properties commonly used in teleoperation systems as in [16, 17]. The main assumptions are presented, while Table 1 shows the nomenclature used in this article.Table 1
Nomenclature.
VariableDefinitionkgxmt−h1=Forward and turn speedsvhiprefvδrefTSent by the operatorzr=v^hipv^δTEstimated real follower speedsused in error dynamicsze=vhipfevδfeTEquivalent smoothed speedsused for force feedbackη=vhipvδTIt is a vector that represents the linearand turning speeds of the equivalent followery˜˙1tLinear speed errory˜2tJoint errorsy˜˙3tTurn speed errorTWalking cycle timeαmLeader dampingσeEquivalent follower dampingσrReal damping of the robotAssumption 1.
The communication channel adds a forward time delayh1 and a backward time delay h2, such that there exist positive values h¯1 and h¯2 such that 0≤h1t≤h¯1 and 0≤h2t≤h¯2 for all t [16].Assumption 2.
The human operator injects a finite amount of energy, which is described as follows:(5)Eh=ϕ−∫0tfhTx˙mdt>0,where ϕ>0 is a positive value [16].Assumption 3.
In the movement of the robot when the swing leg touches the ground, considering a flat horizontal surface, there is no bounce or slip off the swing leg [15, 18].


## 2.1. Leader Model

The leader device is represented in Cartesian coordinates as follows:(1)Mmxmx..m+Cmxm,x.mx.m+gmxm=fm+fh,where xm=xmvxmδT∈Rnm×1 and x.m are the position and velocity of the leader, Mmxm∈Rnm×nm is the inertia matrix, Cmxm,x.m∈Rnm×nm is the matrix representing the centripetal and coriolis forces, gmxm∈Rnm is the gravity vector, fh is the force of the human operator, and fm is the control action applied to the leader.


## 2.2. Follower Model

A bipedal robot can be modeled in a simplified way in a sagittal plane for straight-line walking references. The turn is performed using the hip yaw angles of the 3D bipedal robot, taking into account the invariance property against yaw angle rotation [12]. The presence of impact of the legs on the environment during the walking cycle is represented with hybrid models. These systems, also called systems with impulsive effects, are used to represent a bipedal robot [13, 14], that includes two behaviors: continuous when the leg moves forward, and discrete when the foot touches the ground generating impulsive responses [15].


### 2.2.1. Continuous Dynamic

The continuous dynamics of the bipedal robot between successive impacts presents a total number of physical actuators equal tonb. The state of the bipedal robot is characterized by x=qq.T, where q∈ℝnb is the position and q.∈ℝnb is the velocity. Using the Euler–Lagrange equations, the state between impacts is represented by x.=fx+gxu, as follows:(2)fx=q.−M−1qCq,q.q.+Gq,gx=0M−1qBq,where Mq∈ℝnb×nb represents the inertial matrix, Cq,q.∈ℝnb×nb is the matrix representing the centripetal and Coriolis forces, G∈ℝnb is the gravity vector, B∈ℝnb×nb is the input matrix, and u∈ℝnb is the matrix of control actions. In addition, the angles are defined as follows: q=qp_nshqy_nshqp_nskqp_nsaqp_shqy_shqp_skqp_saT, where qp_nsh is pitch non-stance hip, qp_nsk is pitch non-stance knee, qp_nsa is pitch non-stance ankle, qp_sh is pitch stance hip, qp_sk is pitch stance knee, qp_sa is pitch stance ankle, qy_nsh is yaw non-stance hip, and qy_sh is yaw stance hip, as shown in Figure 3.Figure 3
Joint representation.


### 2.2.2. Domain and Guard

The domainDH specifies the allowed system settings based on the height of the non-stance foot l, measured from a rigid and plain floor, and indicates that it must be above the ground. Additionally, the guard SH represents the impact surface. The domain and guard are defined in reference [13] as follows:(3)DH=x∈R2nb:lq≥0,SH=x∈R2nb:lq=0and∂lq∂qq.<0,where x∈R2nb is the total number of control states given by q and q. and ∂lq/∂q is the partial derivative of l with respect to q.


### 2.2.3. Discrete Dynamic

The discrete dynamics of the robot determine how the robot’s speeds changes when the foot hits the ground, while at the same time switching legs. In particular, the reset mapΔ:SH⟶DH is given by: Δq which represents the switching between the legs when going from “stance” to “no stance” and vice versa after the end of the double support phase, whereas Δq. determines the variation in the speed of each joint due to impact, and therefore, Δ is the vector defined in [13] that links both Δq and Δq.. It is important to mention that a stable continuous dynamic between impacts does not guarantee a stable hybrid system. Hence, it is assumed that there is a dynamic hybrid system represented as follows [14, 15]:(4)W=x.=fx+gxux,ifx∈DHandSH,x+=Δx−,ifx−∈SH.


## 2.2.1. Continuous Dynamic

The continuous dynamics of the bipedal robot between successive impacts presents a total number of physical actuators equal tonb. The state of the bipedal robot is characterized by x=qq.T, where q∈ℝnb is the position and q.∈ℝnb is the velocity. Using the Euler–Lagrange equations, the state between impacts is represented by x.=fx+gxu, as follows:(2)fx=q.−M−1qCq,q.q.+Gq,gx=0M−1qBq,where Mq∈ℝnb×nb represents the inertial matrix, Cq,q.∈ℝnb×nb is the matrix representing the centripetal and Coriolis forces, G∈ℝnb is the gravity vector, B∈ℝnb×nb is the input matrix, and u∈ℝnb is the matrix of control actions. In addition, the angles are defined as follows: q=qp_nshqy_nshqp_nskqp_nsaqp_shqy_shqp_skqp_saT, where qp_nsh is pitch non-stance hip, qp_nsk is pitch non-stance knee, qp_nsa is pitch non-stance ankle, qp_sh is pitch stance hip, qp_sk is pitch stance knee, qp_sa is pitch stance ankle, qy_nsh is yaw non-stance hip, and qy_sh is yaw stance hip, as shown in Figure 3.Figure 3
Joint representation.


## 2.2.2. Domain and Guard

The domainDH specifies the allowed system settings based on the height of the non-stance foot l, measured from a rigid and plain floor, and indicates that it must be above the ground. Additionally, the guard SH represents the impact surface. The domain and guard are defined in reference [13] as follows:(3)DH=x∈R2nb:lq≥0,SH=x∈R2nb:lq=0and∂lq∂qq.<0,where x∈R2nb is the total number of control states given by q and q. and ∂lq/∂q is the partial derivative of l with respect to q.


## 2.2.3. Discrete Dynamic

The discrete dynamics of the robot determine how the robot’s speeds changes when the foot hits the ground, while at the same time switching legs. In particular, the reset mapΔ:SH⟶DH is given by: Δq which represents the switching between the legs when going from “stance” to “no stance” and vice versa after the end of the double support phase, whereas Δq. determines the variation in the speed of each joint due to impact, and therefore, Δ is the vector defined in [13] that links both Δq and Δq.. It is important to mention that a stable continuous dynamic between impacts does not guarantee a stable hybrid system. Hence, it is assumed that there is a dynamic hybrid system represented as follows [14, 15]:(4)W=x.=fx+gxux,ifx∈DHandSH,x+=Δx−,ifx−∈SH.


## 2.3. Assumptions, Properties, and Nomenclature

This document uses the assumptions and properties commonly used in teleoperation systems as in [16, 17]. The main assumptions are presented, while Table 1 shows the nomenclature used in this article.Table 1
Nomenclature.
VariableDefinitionkgxmt−h1=Forward and turn speedsvhiprefvδrefTSent by the operatorzr=v^hipv^δTEstimated real follower speedsused in error dynamicsze=vhipfevδfeTEquivalent smoothed speedsused for force feedbackη=vhipvδTIt is a vector that represents the linearand turning speeds of the equivalent followery˜˙1tLinear speed errory˜2tJoint errorsy˜˙3tTurn speed errorTWalking cycle timeαmLeader dampingσeEquivalent follower dampingσrReal damping of the robotAssumption 1.
The communication channel adds a forward time delayh1 and a backward time delay h2, such that there exist positive values h¯1 and h¯2 such that 0≤h1t≤h¯1 and 0≤h2t≤h¯2 for all t [16].Assumption 2.
The human operator injects a finite amount of energy, which is described as follows:(5)Eh=ϕ−∫0tfhTx˙mdt>0,where ϕ>0 is a positive value [16].Assumption 3.
In the movement of the robot when the swing leg touches the ground, considering a flat horizontal surface, there is no bounce or slip off the swing leg [15, 18].


## 3. Stable Bilateral Controller

The proposed control scheme is shown in Figure4 and is addressed to the bilateral teleoperation of the forward speed and turning of a biped robot. This schematic illustrates the parts of the proposed control scheme for the delayed bilateral teleoperation system. The controller design is oriented to achieve stable hybrid dynamics using Lyapunov-Krasovskii stability theory and Poincaré Map. Therefore, it is defined as an equivalent follower based on the stability of hybrid dynamics and then uses an equivalent follower controller which puts an equivalent damping. The injected leader damping and the equivalent follower damping depend on the time delay. However, the last one cannot be applied directly but in an indirect way through walking cycle time and real follower damping.Figure 4
Complete control scheme: leader controller, time delays, and remote site.


### 3.1. Errors

To teleoperate the bipedal robot, an error vector (for controlling forward speed, gait length, and turning speed) is defined asex:=y˜˙1,y˜2,y˜2.,y˜˙3,y˜4v,y˜4δT, where their components will be described next.First, the forward speed error is defined as follows:(6)y˜˙1t=y˙1t−vhipt,where y˙1 is the forward speed reference set as follows:(7)y˙1t=kgvxmvt−h1=vhipref,where kgv linearly maps the position of the leader to the forward velocity and xmv is the position of the leader. Furthermore, vhip is the real speed of the humanoid in the X-axis obtained by linearization. The position of the hip depends on the pitch angle of the posture ankle qp_sa and the pitch angle of the posture knee qp_sk, and its time derivative is calculated as in reference [13] as follows:(8)phipq=a1−qp_sa+a2−qp_sa−qp_sk,vhip=ddtphipq,where phip is the linearization of the X-position of the hip.To handle the walking cycle to obtain shorter or larger gaits, it is necessary to control the error between the reference values of the joints angles and the real values obtained by the sensors of the humanoid as a gait length error represented as follows:(9)y˜2t=y2t−yr,2t,where y2t=qp_nsaqp_saqp_nskqp_skqp_nshqp_shT are the joint angles that define the gait following the reference. yr,2t=qp_nsarqp_sarqp_nskrqp_skrqp_nshrqp_shrT are the real joint angles. To obtain the references of stance knee angle qp_sk and stance hip angle qp_sh, during a walking cycle time, a procedure based on the use of Bezier curves is applied as in reference [10]. Taking into account the percentages of the single support and double support phases during a walking cycle, and considering that the legs are symmetrical, the angles qp_nsh and qp_nsk are calculated using a phase shift of the angles qp_sh and qpsk. The ankle angles for the stance and non-stance legs are obtained considering that the torso always stays perpendicular to the ground during the walk. Considering this, we use formula (15) [13] to find the ankle angles as: qp_nsa=−qp_nsh+qp_nsk and qp_sa=−qp_sh+qp_sk. Besides, the turning speed error is described as follows:(10)y˜˙3q=y.3t−vδt,where y.3 is the turning speed reference and it is given as follows:(11)y˙3t=kgδxmδt−h1=vδref,where kgδ linearly maps the position of the leader to the turning velocity and xmδ is the position of the leader. In addition, vδt=fqy_sh,qy_nsh is the actual turn speed of the humanoid robot, and it is obtained by differentiating the turn angle of the robot δt. This yaw rotation reference does not affect the control in the sagittal plane due to the invariance property as stated in [12], where the within-stride feedback does not depend on the yaw orientation of the robot. The references for the angles qy_sh,qy_nsh, depend on the turn speed reference and are obtained through a sequence illustrated in Figure 5 based on [19]. The variable φ is the percentage of the walking cycle that both feet spend in the ground at the same time (double support) and is set in practice to get a trade-off between a greater average speed (lower φ) and a better balance (higher φ). For the humanoid to turn a given angle, we determine the phases of a single and double support. Half of the total turn angle is sent to the hip’s yaw angle of the support leg during the single support phase. In order to achieve the same turn direction in each leg, the angles are sent with opposite signs due to the orientation of the motors of each joint. In order to perform a successful turn, each hip joint performs the turn when its corresponding leg is used as support. Once the single support phase ends, the hip joint stops moving so the robot will not fall. Once a total walking cycle has been performed, the process starts again until the robot has turned the desired angle.Figure 5
Turning sequence.Finally, the error componenty˜4 is defined by: y˜4=η−zr=y˜4vy˜4δ=vhip−v^hipvδ−v^δ. Where zr=v^hipv^δT is the estimated forward and turn speed and it is found by using the following linear observer:(12)z˙r=βr−1η−zr,where βr−1>0 is the observer gain.


### 3.2. Stability of Follower between Impacts

#### 3.2.1. Closed Loop Follower Control

The continuous part of the model (4) can be used to represent the auxiliary variable yx through the notation of Lie’s derivative as follows:(13)y..x=Lf2yx+LgLfyxux,where yx is designed to work simultaneously with speed and position errors as yx=∫y˜˙1y˜2∫y˜˙3y˜4T.The feedback linearization controller based on reference [14] is set as follows:(14)ux=LgLfyx+−Lf2yx−υx,where LgLfyx+ is the corresponding pseudo-inverse matrix. A disadvantage of formula (14) is that it is based on a perfect compensation of the robot dynamics. If uncertainties in the dynamics are considered and it is assumed that the functions fx,gx of formula (2) are estimated, the controller must be designed based on the estimated functions f^x,g^x. Thus, the law of control (14) is rewritten as follows [20]:(15)ux=Lg^Lf^yx+−Lf^2yx−υx.Hence, replacing formula (15) in formula (13), it is obtained as follows:(16)y..x=−υx−W,where(17)W=Δ1+Δ2υ,Δ1=Lf2yx−LgLfyxLg^Lf^yx+Lf^2yx,Δ2=−LgLfyxLg^Lf^yx+Lf^2yx−I.To compensate for the uncertainties, a combined controller is applied withυ=υ1+υ2, where the first part υ1 allows to follow the reference model considering a perfect knowledge of inverse dynamics, and the second part υ2 compensates for nonlinear uncertainties W.From formula (16), the closed-loop dynamics of e is represented in state space as follows:(18)e.=Fe+υ1+W,in which we set υ1=Ke+Kvfv and υ2=−W^. Therefore, W˜=W−W^=w˜v0w˜gw˜δw˜4vw˜4δ. A nonlinear dynamics compensation technique such as robust adaptive control [20–22], neural networks [23], and deep learning, among others, must be used to estimate W^, where we assumed that W˜ is bounded. As the compensator gets better, such errors will be lower.Describing the terms of (18) as follows:(19)F=00000000I000000000000000000000000000,Kv=kzv00000000000000000000kzδ00000000000000,K=−M−1kv000M−1σβr00000000−1ε2−1ε000000−I−1kδ0I−1τβrM−1kv000−1βrM−1σ+10000I−1kδ0−1βrI−1τ+1,where σ,τ are the real injected damping σr=σ00τ, kv,kδ are the proportional gains ks=kv00kδ, ε>0 is a control gain, M is the mass of the reference model, I is the inertia of the reference model D=M00I, and kz=kzv00kzδ represents the impedance applied to the virtual force fv which is computed using the distance and orientation measured between the robot’s ultrasonic sensors and the objects close to it. Then, fv is used to reduce the robot’s advance speed as it gets near to the obstacle.The robot has a fictitious force to avoid independently obstacles. To detect obstacles, NAO is equipped with two ultrasonic sensors (or sonars), which allow it to estimate the distance to obstacles in its environment. The fictitious forcefv is computed with the distance and orientation measured between the robot ultrasonic sensors and the objects close to it (obstacles). Then, fv is defined according to its tangential and normal components, where only the tangential component is going to be used to reduce the robot’s forward speed, providing the operator the opportunity to apply a turning speed to evade obstacles.


#### 3.2.2. Stability of Hybrid Dynamics Follower

Solving formula (18) and considering that the reference forward and turning speeds are kept constant between impacts and considering a feedforward of the reference, the following dynamic equations are obtained:(20)Mv˙hip=kvy˜˙1−σv˙hipf−kzvfv−Mw˜v,(21)y˜2..=−1εy˜2.−1ε2y˜2+w˜g,(22)Iv˙δ=kδy˜˙3−τv˙δf−kzδfv−Iw˜δ,(23)y˜4.=y˜˙4vy˜˙4δ=M−1kvy˜˙1−1βrM−1σ+1y˜4vI−1kδy˜˙3−1βrI−1τ+1y˜4δ+w˜4vw˜4δ.Therefore, from formulas (20)–(23), we infer that ex⟶0 considering null both W˜ and fv, or in general the error tends to a ball centered in the origin. In order to analyze the stability of the hybrid system, these disturbances and the fictitious force will not be considered, but in reference [24], the authors studied how bounded perturbations affect the hybrid stability. Hence, to analyze the stability of ex considering the hybrid dynamics, a Poincaré Map is used which describes the relationship of the states between impacts. The stability of the map determines the stability of the periodic orbit O. Besides, the exponential stability of a periodic orbit in a nonlinear system with impulse effects can be studied by linearizing the Poincaré return map around a fixed point and evaluating its eigenvalues. The Poincaré map is represented as follows:(24)ek+1=Pek.So, the fixed pointe∗ is locally exponentially stable if and only if the periodic orbit O is locally exponentially stable, and this is achieved if the eigenvalue λi module of (24) lie inside the unit circle. Following the procedure of linearization described in reference [25], it is applied to formula (24) around the fixed point e∗∈S, we get the following equation:(25)ϕek+1=Aϕek,where ϕek=ek−e∗ and A=A0A1A2…Antnt×nt is the Jacobian of the Poincaré Map, with i=1,…,nt, nt=2ny2+np, where the states for position and speed of the vector y˜2 are ny2, the number of states to be disturbed is np and is corresponding to the forward speed error y˜˙1, turning speed error y˜˙3, and the states y˜4v and y˜4δ.Therefore, the Jacobian is calculated as follows:(26)Ai=Pe∗+Δei−Pe∗−Δei2Δei,where the value of Δe depends on the states: Δx=ΔqiΔq˙i and another vector is defined as Δp=Δpvhipiq˙iΔpvδiq˙iΔpy˜4viq˙iΔpy˜4δiq˙inp×1. The values of Δqi,Δq˙i,Δpvhipiq˙i, Δpvδiq˙i, Δpy˜4viq˙i, and Δpy˜4δiq˙i are small perturbations introduced to calculate the linearized model (26), performing nt evaluations of Pe∗±Δei [18]. Linking all parameters that affect the dynamics of the error within a single vector is obtained as follows:(27)ci=ΔxΔpDksβr,where we set D=4.5001, ks=1001, βr=0.5000.5, ε=0.5, vhipref=5cm/s, vδref=4°/s, Δqi=1, Δq˙i=0.375, and Δpvhipi=0.7 to numerically evaluate the Poincaré Map whose linearization is calculated using the fixed point e∗ and the parameters shown previously. This fixed point corresponds to the point where the periodic orbit of the proposed system intersects the reset surface. Using this fixed point and taking different values of Δp (Δpvhip depends on the type of terrain, the mass of the robot, and among others), the walking cycle time, and the real damping and constant keeping the rest of parameters of ci, the Jacobian’s eigenvalues of the Poincaré Map are calculated.It is possible to get a similarλhybrid by applying different combinations of real follower damping and walking cycle time, as shown in Figure 6. As the walking cycle time increases for a given real damping value, the eigenvalues decrease in magnitude and generate a more stable behavior. Also, as the walking cycle time decreases and the value Δpvhip moves away from 1, the real damping should decrease. Therefore, it is possible to apply a wider range of walking cycle time and real damping when the discontinuity caused by the reset map is smaller (Δpvhip closer to one). Additionally, it should be noted that more surfaces can be obtained by changing the other parameters of ci like D, βr, among others.Figure 6
Stability surfaces, is visualized as the real damping and walking cycle time change, taking differentΔp: (a) Δpvhip=0.25, (b) Δpvhip=0.5, and (c) Δpvhip=0.7.
(a)(b)(c)


### 3.3. Equivalent Follower and Equivalent Follower Controller

From the equations of the hybrid real follower between impacts, shown in equations (20) and (22), the following expression is derived:(28)Dη˙=kskgxmt−h1−zr−σrz˙r−kzfv−w˜l,where w˜l=w˜vw˜δT is the nonlinear uncertainty that sustains w˜l∈ℒ∞ and kg=kgv00kgδ. We define a linear continuous time follower named equivalent follower for all t that has the following constraints: Δp=Δpmin,σr=σrmin,ε=εmin,σe=σemin, such that the values of βe,Tmin are defined as follows:


#### 3.3.1. Definition

Definition 1.
Tmin:=minT/maxdi,Δp=Δpmin,σr=σrmin,Tmin,ε=εminλhybrid=1Definition 2.
σrmin:=minσr/maxdi,Δp=Δpmin,Tmin,ε=εminλhybrid=1Definition 3.
βe:=minβe/maxdi,σe=σemin,βeλequivalent=1 for the nonperturbed discrete-time version of formula (28):(29)Dη˙=kskgxmt−h1−ze−σez˙e,ze=η−βez˙e,where λhybrid are the eigenvalues of the hybrid system computed from the Jacobian’s eigenvalues of the Poincaré map. λequivaalent are the eigenvalues of the discretized version of (29), di=Dksβr, σe is the linear damping coefficient injected by the equivalent follower controller, and ze is the velocity of the equivalent follower. The values σrmin,εmin,σemin are positive values near to zero while 0<Δpmin<1 establish the range of the reset map from Δp=Δpmin,1.
Due to the effect caused by the increasing of equivalent damping (equivalent follower controller),λequivalent decreases. It is possible to get a similar λhybrid by applying different combinations of real follower damping and walking cycle time, as shown in Figure 7. Now, we propose to online set σr and T through the following optimization criteria, visually illustrated in Figure 7:(30)σr,T,ε=minkσσr−σe+kwtT+kεεsuch that:maxci,σr,T,ελhybrid=maxdi,σeh1,h2λequivalent,where kσ,kwt,kε are the weights of the functional. This optimization cost function aims to find a trade-off between a walking velocity as near as possible to the velocity of the non-delay case and get the lowest error between the follower equivalent velocity ze and the follower real velocity zr. Finally, from Figure 7, it can be seen that σr approaches σe to ensure stable behavior. Part (a) of Figure 7 shows that if h↑⇒σe↑ therefore λequivalent decreases. However, without a proper control λhybrid will be greater than λequivalent, which does not guarantee the stability of the delayed bilateral teleoperation system. Part (b) of Figure 7 states that if h↑⇒σe↑ therefore λequivalent decreases and if h↑⇒T↑∧σr↑ therefore, λhybrid decreases in the same way as the equivalent, which allows to hold a stable behavior assuring bounded errors.Figure 7
Evolution of the behavior of the equivalent follower robot (continuous system) and the humanoid robot (hybrid system).


### 3.4. Stable Bilateral Controller

In this work, aP + d such as the controller is used to analyze what conditions must be held to keep the errors ec=kgxm−ze bounded. The proposed control strategy adds the following force feedback law for fm:(31)fm=−kmkgxm−zet−h2+gmxm−αmx˙m,where km is a gain matrix, kg linearly maps the leader’s position to a reference forward and turn speeds, and αm is the linear damping coefficient injected into the leader. The only difference with a traditional P + d scheme is that the equivalent follower velocity is used over the real one. Besides, the control actions applied to the real follower are given by formula (15).


### 3.5. Stability Analysis

Here, the stability of the proposed teleoperation system is analyzed and the process of calibrating the control parameters is explained. Iffm (31) is included in the leader (1) and with the equivalent robot (29), the following closed-loop dynamics are obtained:(32)x¨m=Mm−1−Cmx˙m−kmkgxm−zet−h2−αmx˙m+fh,(33)η˙=D−1kskgxmt−h1−ze−σez˙e−kzfv−w˜l.According to formulas (32) and (33), and taking also into account (5), it can be seen that x¨m=0,η˙=0, if the vector xt:=x˙mzez˙eecT remains zero for a sufficiently high time interval, considering fh, fv, and w˜l null and whose analysis is associated with Lyapunov–Krasovskii’s concept of stability, where there is a continuous of states that we represent using integrals with finite limits.A positive definite functionVxt=V1+V2+V3+V4 is given as follows:(34)V1=12l1z˙eTDz˙e,(35)V2=12l2x˙mTMmxmx˙m+l2Eh,(36)V3=12l3ecTec+12l4zeTze,(37)V4=l5∫−h¯20∫t+θtz˙eξTz˙eξdξdθ+l6∫−h¯10∫t+θtx˙mξTx˙mξdξdθ.Following the procedure described in reference [11] along the trajectories of the closed-loop system, considering the follower robot, leader dynamics, time delay, and human operator and environment forces, is bounded by:(38)V˙≤−λmx˙mTx˙m−λseze.Tz˙e−ρez˙e,where(39)ρe=w˜l+kzf¯v.Ifαm and σe are sufficiently high to comply as follows:(40)λm=αmkm−h1−h24>0,λse=σe+Dks−14h1−h2>0,then, the variables x˙m,z˙e∈ℒ∞. If formula (38) is integrated in time, the following is verified:(41)Vt−V0≤−λmx˙m22−λsez˙e22−∫0tz˙eϵTρeϵdϵ.


#### 3.5.1. Remark

The term∫0tz˙eϵTρeϵdϵ is bounded since x˙m,z˙e∈ℒ∞ and using the proof of Section 3.3, it is concluded that the estimation errors are limited by W^∈ℒ∞. Therefore, from formula (41), we infer that Vt will be bounded for all t, and therefore, x˙m,z˙e∈ℒ2 and kgxm−ze,ze,xm,∈ℒ∞. If fv,fh,e˜l remain zero, and following an analysis similar to reference [16] based on Barbalat’s Lemma, it is possible to deduce that x˙m and z˙e will tend to zero as t⟶∞.


### 3.6. Cascade Calibration Process

A three-step calibration process is used to set the controller parameters, linking the stability result obtained by using the equivalent follower with the equivalent follower controller, as well as the relation between the equivalent dynamics and the hybrid dynamics as follows:(i)
V.1 Offline calibration step, the inverse dynamics compensation is calibrated, including their gains that depend on the used method, in particular, we apply an adaptive inverse dynamic compensation [20]. Furthermore, the reference parameters D and the gains of the bilateral controller (ks, km) without delay are adjusted while the observer parameter βr is set empirically. Finally, βe is set based on Definition 1.(ii)
V.2 Online step, the leader dampingαm and equivalent follower damping σe are calculated from formula (40).(iii)
V.3 Online step, the values of walking cycle timeT and real damping σr are set through the expression (30).


## 3.1. Errors

To teleoperate the bipedal robot, an error vector (for controlling forward speed, gait length, and turning speed) is defined asex:=y˜˙1,y˜2,y˜2.,y˜˙3,y˜4v,y˜4δT, where their components will be described next.First, the forward speed error is defined as follows:(6)y˜˙1t=y˙1t−vhipt,where y˙1 is the forward speed reference set as follows:(7)y˙1t=kgvxmvt−h1=vhipref,where kgv linearly maps the position of the leader to the forward velocity and xmv is the position of the leader. Furthermore, vhip is the real speed of the humanoid in the X-axis obtained by linearization. The position of the hip depends on the pitch angle of the posture ankle qp_sa and the pitch angle of the posture knee qp_sk, and its time derivative is calculated as in reference [13] as follows:(8)phipq=a1−qp_sa+a2−qp_sa−qp_sk,vhip=ddtphipq,where phip is the linearization of the X-position of the hip.To handle the walking cycle to obtain shorter or larger gaits, it is necessary to control the error between the reference values of the joints angles and the real values obtained by the sensors of the humanoid as a gait length error represented as follows:(9)y˜2t=y2t−yr,2t,where y2t=qp_nsaqp_saqp_nskqp_skqp_nshqp_shT are the joint angles that define the gait following the reference. yr,2t=qp_nsarqp_sarqp_nskrqp_skrqp_nshrqp_shrT are the real joint angles. To obtain the references of stance knee angle qp_sk and stance hip angle qp_sh, during a walking cycle time, a procedure based on the use of Bezier curves is applied as in reference [10]. Taking into account the percentages of the single support and double support phases during a walking cycle, and considering that the legs are symmetrical, the angles qp_nsh and qp_nsk are calculated using a phase shift of the angles qp_sh and qpsk. The ankle angles for the stance and non-stance legs are obtained considering that the torso always stays perpendicular to the ground during the walk. Considering this, we use formula (15) [13] to find the ankle angles as: qp_nsa=−qp_nsh+qp_nsk and qp_sa=−qp_sh+qp_sk. Besides, the turning speed error is described as follows:(10)y˜˙3q=y.3t−vδt,where y.3 is the turning speed reference and it is given as follows:(11)y˙3t=kgδxmδt−h1=vδref,where kgδ linearly maps the position of the leader to the turning velocity and xmδ is the position of the leader. In addition, vδt=fqy_sh,qy_nsh is the actual turn speed of the humanoid robot, and it is obtained by differentiating the turn angle of the robot δt. This yaw rotation reference does not affect the control in the sagittal plane due to the invariance property as stated in [12], where the within-stride feedback does not depend on the yaw orientation of the robot. The references for the angles qy_sh,qy_nsh, depend on the turn speed reference and are obtained through a sequence illustrated in Figure 5 based on [19]. The variable φ is the percentage of the walking cycle that both feet spend in the ground at the same time (double support) and is set in practice to get a trade-off between a greater average speed (lower φ) and a better balance (higher φ). For the humanoid to turn a given angle, we determine the phases of a single and double support. Half of the total turn angle is sent to the hip’s yaw angle of the support leg during the single support phase. In order to achieve the same turn direction in each leg, the angles are sent with opposite signs due to the orientation of the motors of each joint. In order to perform a successful turn, each hip joint performs the turn when its corresponding leg is used as support. Once the single support phase ends, the hip joint stops moving so the robot will not fall. Once a total walking cycle has been performed, the process starts again until the robot has turned the desired angle.Figure 5
Turning sequence.Finally, the error componenty˜4 is defined by: y˜4=η−zr=y˜4vy˜4δ=vhip−v^hipvδ−v^δ. Where zr=v^hipv^δT is the estimated forward and turn speed and it is found by using the following linear observer:(12)z˙r=βr−1η−zr,where βr−1>0 is the observer gain.


## 3.2. Stability of Follower between Impacts

### 3.2.1. Closed Loop Follower Control

The continuous part of the model (4) can be used to represent the auxiliary variable yx through the notation of Lie’s derivative as follows:(13)y..x=Lf2yx+LgLfyxux,where yx is designed to work simultaneously with speed and position errors as yx=∫y˜˙1y˜2∫y˜˙3y˜4T.The feedback linearization controller based on reference [14] is set as follows:(14)ux=LgLfyx+−Lf2yx−υx,where LgLfyx+ is the corresponding pseudo-inverse matrix. A disadvantage of formula (14) is that it is based on a perfect compensation of the robot dynamics. If uncertainties in the dynamics are considered and it is assumed that the functions fx,gx of formula (2) are estimated, the controller must be designed based on the estimated functions f^x,g^x. Thus, the law of control (14) is rewritten as follows [20]:(15)ux=Lg^Lf^yx+−Lf^2yx−υx.Hence, replacing formula (15) in formula (13), it is obtained as follows:(16)y..x=−υx−W,where(17)W=Δ1+Δ2υ,Δ1=Lf2yx−LgLfyxLg^Lf^yx+Lf^2yx,Δ2=−LgLfyxLg^Lf^yx+Lf^2yx−I.To compensate for the uncertainties, a combined controller is applied withυ=υ1+υ2, where the first part υ1 allows to follow the reference model considering a perfect knowledge of inverse dynamics, and the second part υ2 compensates for nonlinear uncertainties W.From formula (16), the closed-loop dynamics of e is represented in state space as follows:(18)e.=Fe+υ1+W,in which we set υ1=Ke+Kvfv and υ2=−W^. Therefore, W˜=W−W^=w˜v0w˜gw˜δw˜4vw˜4δ. A nonlinear dynamics compensation technique such as robust adaptive control [20–22], neural networks [23], and deep learning, among others, must be used to estimate W^, where we assumed that W˜ is bounded. As the compensator gets better, such errors will be lower.Describing the terms of (18) as follows:(19)F=00000000I000000000000000000000000000,Kv=kzv00000000000000000000kzδ00000000000000,K=−M−1kv000M−1σβr00000000−1ε2−1ε000000−I−1kδ0I−1τβrM−1kv000−1βrM−1σ+10000I−1kδ0−1βrI−1τ+1,where σ,τ are the real injected damping σr=σ00τ, kv,kδ are the proportional gains ks=kv00kδ, ε>0 is a control gain, M is the mass of the reference model, I is the inertia of the reference model D=M00I, and kz=kzv00kzδ represents the impedance applied to the virtual force fv which is computed using the distance and orientation measured between the robot’s ultrasonic sensors and the objects close to it. Then, fv is used to reduce the robot’s advance speed as it gets near to the obstacle.The robot has a fictitious force to avoid independently obstacles. To detect obstacles, NAO is equipped with two ultrasonic sensors (or sonars), which allow it to estimate the distance to obstacles in its environment. The fictitious forcefv is computed with the distance and orientation measured between the robot ultrasonic sensors and the objects close to it (obstacles). Then, fv is defined according to its tangential and normal components, where only the tangential component is going to be used to reduce the robot’s forward speed, providing the operator the opportunity to apply a turning speed to evade obstacles.


### 3.2.2. Stability of Hybrid Dynamics Follower

Solving formula (18) and considering that the reference forward and turning speeds are kept constant between impacts and considering a feedforward of the reference, the following dynamic equations are obtained:(20)Mv˙hip=kvy˜˙1−σv˙hipf−kzvfv−Mw˜v,(21)y˜2..=−1εy˜2.−1ε2y˜2+w˜g,(22)Iv˙δ=kδy˜˙3−τv˙δf−kzδfv−Iw˜δ,(23)y˜4.=y˜˙4vy˜˙4δ=M−1kvy˜˙1−1βrM−1σ+1y˜4vI−1kδy˜˙3−1βrI−1τ+1y˜4δ+w˜4vw˜4δ.Therefore, from formulas (20)–(23), we infer that ex⟶0 considering null both W˜ and fv, or in general the error tends to a ball centered in the origin. In order to analyze the stability of the hybrid system, these disturbances and the fictitious force will not be considered, but in reference [24], the authors studied how bounded perturbations affect the hybrid stability. Hence, to analyze the stability of ex considering the hybrid dynamics, a Poincaré Map is used which describes the relationship of the states between impacts. The stability of the map determines the stability of the periodic orbit O. Besides, the exponential stability of a periodic orbit in a nonlinear system with impulse effects can be studied by linearizing the Poincaré return map around a fixed point and evaluating its eigenvalues. The Poincaré map is represented as follows:(24)ek+1=Pek.So, the fixed pointe∗ is locally exponentially stable if and only if the periodic orbit O is locally exponentially stable, and this is achieved if the eigenvalue λi module of (24) lie inside the unit circle. Following the procedure of linearization described in reference [25], it is applied to formula (24) around the fixed point e∗∈S, we get the following equation:(25)ϕek+1=Aϕek,where ϕek=ek−e∗ and A=A0A1A2…Antnt×nt is the Jacobian of the Poincaré Map, with i=1,…,nt, nt=2ny2+np, where the states for position and speed of the vector y˜2 are ny2, the number of states to be disturbed is np and is corresponding to the forward speed error y˜˙1, turning speed error y˜˙3, and the states y˜4v and y˜4δ.Therefore, the Jacobian is calculated as follows:(26)Ai=Pe∗+Δei−Pe∗−Δei2Δei,where the value of Δe depends on the states: Δx=ΔqiΔq˙i and another vector is defined as Δp=Δpvhipiq˙iΔpvδiq˙iΔpy˜4viq˙iΔpy˜4δiq˙inp×1. The values of Δqi,Δq˙i,Δpvhipiq˙i, Δpvδiq˙i, Δpy˜4viq˙i, and Δpy˜4δiq˙i are small perturbations introduced to calculate the linearized model (26), performing nt evaluations of Pe∗±Δei [18]. Linking all parameters that affect the dynamics of the error within a single vector is obtained as follows:(27)ci=ΔxΔpDksβr,where we set D=4.5001, ks=1001, βr=0.5000.5, ε=0.5, vhipref=5cm/s, vδref=4°/s, Δqi=1, Δq˙i=0.375, and Δpvhipi=0.7 to numerically evaluate the Poincaré Map whose linearization is calculated using the fixed point e∗ and the parameters shown previously. This fixed point corresponds to the point where the periodic orbit of the proposed system intersects the reset surface. Using this fixed point and taking different values of Δp (Δpvhip depends on the type of terrain, the mass of the robot, and among others), the walking cycle time, and the real damping and constant keeping the rest of parameters of ci, the Jacobian’s eigenvalues of the Poincaré Map are calculated.It is possible to get a similarλhybrid by applying different combinations of real follower damping and walking cycle time, as shown in Figure 6. As the walking cycle time increases for a given real damping value, the eigenvalues decrease in magnitude and generate a more stable behavior. Also, as the walking cycle time decreases and the value Δpvhip moves away from 1, the real damping should decrease. Therefore, it is possible to apply a wider range of walking cycle time and real damping when the discontinuity caused by the reset map is smaller (Δpvhip closer to one). Additionally, it should be noted that more surfaces can be obtained by changing the other parameters of ci like D, βr, among others.Figure 6
Stability surfaces, is visualized as the real damping and walking cycle time change, taking differentΔp: (a) Δpvhip=0.25, (b) Δpvhip=0.5, and (c) Δpvhip=0.7.
(a)(b)(c)


## 3.2.1. Closed Loop Follower Control

The continuous part of the model (4) can be used to represent the auxiliary variable yx through the notation of Lie’s derivative as follows:(13)y..x=Lf2yx+LgLfyxux,where yx is designed to work simultaneously with speed and position errors as yx=∫y˜˙1y˜2∫y˜˙3y˜4T.The feedback linearization controller based on reference [14] is set as follows:(14)ux=LgLfyx+−Lf2yx−υx,where LgLfyx+ is the corresponding pseudo-inverse matrix. A disadvantage of formula (14) is that it is based on a perfect compensation of the robot dynamics. If uncertainties in the dynamics are considered and it is assumed that the functions fx,gx of formula (2) are estimated, the controller must be designed based on the estimated functions f^x,g^x. Thus, the law of control (14) is rewritten as follows [20]:(15)ux=Lg^Lf^yx+−Lf^2yx−υx.Hence, replacing formula (15) in formula (13), it is obtained as follows:(16)y..x=−υx−W,where(17)W=Δ1+Δ2υ,Δ1=Lf2yx−LgLfyxLg^Lf^yx+Lf^2yx,Δ2=−LgLfyxLg^Lf^yx+Lf^2yx−I.To compensate for the uncertainties, a combined controller is applied withυ=υ1+υ2, where the first part υ1 allows to follow the reference model considering a perfect knowledge of inverse dynamics, and the second part υ2 compensates for nonlinear uncertainties W.From formula (16), the closed-loop dynamics of e is represented in state space as follows:(18)e.=Fe+υ1+W,in which we set υ1=Ke+Kvfv and υ2=−W^. Therefore, W˜=W−W^=w˜v0w˜gw˜δw˜4vw˜4δ. A nonlinear dynamics compensation technique such as robust adaptive control [20–22], neural networks [23], and deep learning, among others, must be used to estimate W^, where we assumed that W˜ is bounded. As the compensator gets better, such errors will be lower.Describing the terms of (18) as follows:(19)F=00000000I000000000000000000000000000,Kv=kzv00000000000000000000kzδ00000000000000,K=−M−1kv000M−1σβr00000000−1ε2−1ε000000−I−1kδ0I−1τβrM−1kv000−1βrM−1σ+10000I−1kδ0−1βrI−1τ+1,where σ,τ are the real injected damping σr=σ00τ, kv,kδ are the proportional gains ks=kv00kδ, ε>0 is a control gain, M is the mass of the reference model, I is the inertia of the reference model D=M00I, and kz=kzv00kzδ represents the impedance applied to the virtual force fv which is computed using the distance and orientation measured between the robot’s ultrasonic sensors and the objects close to it. Then, fv is used to reduce the robot’s advance speed as it gets near to the obstacle.The robot has a fictitious force to avoid independently obstacles. To detect obstacles, NAO is equipped with two ultrasonic sensors (or sonars), which allow it to estimate the distance to obstacles in its environment. The fictitious forcefv is computed with the distance and orientation measured between the robot ultrasonic sensors and the objects close to it (obstacles). Then, fv is defined according to its tangential and normal components, where only the tangential component is going to be used to reduce the robot’s forward speed, providing the operator the opportunity to apply a turning speed to evade obstacles.


## 3.2.2. Stability of Hybrid Dynamics Follower

Solving formula (18) and considering that the reference forward and turning speeds are kept constant between impacts and considering a feedforward of the reference, the following dynamic equations are obtained:(20)Mv˙hip=kvy˜˙1−σv˙hipf−kzvfv−Mw˜v,(21)y˜2..=−1εy˜2.−1ε2y˜2+w˜g,(22)Iv˙δ=kδy˜˙3−τv˙δf−kzδfv−Iw˜δ,(23)y˜4.=y˜˙4vy˜˙4δ=M−1kvy˜˙1−1βrM−1σ+1y˜4vI−1kδy˜˙3−1βrI−1τ+1y˜4δ+w˜4vw˜4δ.Therefore, from formulas (20)–(23), we infer that ex⟶0 considering null both W˜ and fv, or in general the error tends to a ball centered in the origin. In order to analyze the stability of the hybrid system, these disturbances and the fictitious force will not be considered, but in reference [24], the authors studied how bounded perturbations affect the hybrid stability. Hence, to analyze the stability of ex considering the hybrid dynamics, a Poincaré Map is used which describes the relationship of the states between impacts. The stability of the map determines the stability of the periodic orbit O. Besides, the exponential stability of a periodic orbit in a nonlinear system with impulse effects can be studied by linearizing the Poincaré return map around a fixed point and evaluating its eigenvalues. The Poincaré map is represented as follows:(24)ek+1=Pek.So, the fixed pointe∗ is locally exponentially stable if and only if the periodic orbit O is locally exponentially stable, and this is achieved if the eigenvalue λi module of (24) lie inside the unit circle. Following the procedure of linearization described in reference [25], it is applied to formula (24) around the fixed point e∗∈S, we get the following equation:(25)ϕek+1=Aϕek,where ϕek=ek−e∗ and A=A0A1A2…Antnt×nt is the Jacobian of the Poincaré Map, with i=1,…,nt, nt=2ny2+np, where the states for position and speed of the vector y˜2 are ny2, the number of states to be disturbed is np and is corresponding to the forward speed error y˜˙1, turning speed error y˜˙3, and the states y˜4v and y˜4δ.Therefore, the Jacobian is calculated as follows:(26)Ai=Pe∗+Δei−Pe∗−Δei2Δei,where the value of Δe depends on the states: Δx=ΔqiΔq˙i and another vector is defined as Δp=Δpvhipiq˙iΔpvδiq˙iΔpy˜4viq˙iΔpy˜4δiq˙inp×1. The values of Δqi,Δq˙i,Δpvhipiq˙i, Δpvδiq˙i, Δpy˜4viq˙i, and Δpy˜4δiq˙i are small perturbations introduced to calculate the linearized model (26), performing nt evaluations of Pe∗±Δei [18]. Linking all parameters that affect the dynamics of the error within a single vector is obtained as follows:(27)ci=ΔxΔpDksβr,where we set D=4.5001, ks=1001, βr=0.5000.5, ε=0.5, vhipref=5cm/s, vδref=4°/s, Δqi=1, Δq˙i=0.375, and Δpvhipi=0.7 to numerically evaluate the Poincaré Map whose linearization is calculated using the fixed point e∗ and the parameters shown previously. This fixed point corresponds to the point where the periodic orbit of the proposed system intersects the reset surface. Using this fixed point and taking different values of Δp (Δpvhip depends on the type of terrain, the mass of the robot, and among others), the walking cycle time, and the real damping and constant keeping the rest of parameters of ci, the Jacobian’s eigenvalues of the Poincaré Map are calculated.It is possible to get a similarλhybrid by applying different combinations of real follower damping and walking cycle time, as shown in Figure 6. As the walking cycle time increases for a given real damping value, the eigenvalues decrease in magnitude and generate a more stable behavior. Also, as the walking cycle time decreases and the value Δpvhip moves away from 1, the real damping should decrease. Therefore, it is possible to apply a wider range of walking cycle time and real damping when the discontinuity caused by the reset map is smaller (Δpvhip closer to one). Additionally, it should be noted that more surfaces can be obtained by changing the other parameters of ci like D, βr, among others.Figure 6
Stability surfaces, is visualized as the real damping and walking cycle time change, taking differentΔp: (a) Δpvhip=0.25, (b) Δpvhip=0.5, and (c) Δpvhip=0.7.
(a)(b)(c)


## 3.3. Equivalent Follower and Equivalent Follower Controller

From the equations of the hybrid real follower between impacts, shown in equations (20) and (22), the following expression is derived:(28)Dη˙=kskgxmt−h1−zr−σrz˙r−kzfv−w˜l,where w˜l=w˜vw˜δT is the nonlinear uncertainty that sustains w˜l∈ℒ∞ and kg=kgv00kgδ. We define a linear continuous time follower named equivalent follower for all t that has the following constraints: Δp=Δpmin,σr=σrmin,ε=εmin,σe=σemin, such that the values of βe,Tmin are defined as follows:


### 3.3.1. Definition

Definition 1.
Tmin:=minT/maxdi,Δp=Δpmin,σr=σrmin,Tmin,ε=εminλhybrid=1Definition 2.
σrmin:=minσr/maxdi,Δp=Δpmin,Tmin,ε=εminλhybrid=1Definition 3.
βe:=minβe/maxdi,σe=σemin,βeλequivalent=1 for the nonperturbed discrete-time version of formula (28):(29)Dη˙=kskgxmt−h1−ze−σez˙e,ze=η−βez˙e,where λhybrid are the eigenvalues of the hybrid system computed from the Jacobian’s eigenvalues of the Poincaré map. λequivaalent are the eigenvalues of the discretized version of (29), di=Dksβr, σe is the linear damping coefficient injected by the equivalent follower controller, and ze is the velocity of the equivalent follower. The values σrmin,εmin,σemin are positive values near to zero while 0<Δpmin<1 establish the range of the reset map from Δp=Δpmin,1.
Due to the effect caused by the increasing of equivalent damping (equivalent follower controller),λequivalent decreases. It is possible to get a similar λhybrid by applying different combinations of real follower damping and walking cycle time, as shown in Figure 7. Now, we propose to online set σr and T through the following optimization criteria, visually illustrated in Figure 7:(30)σr,T,ε=minkσσr−σe+kwtT+kεεsuch that:maxci,σr,T,ελhybrid=maxdi,σeh1,h2λequivalent,where kσ,kwt,kε are the weights of the functional. This optimization cost function aims to find a trade-off between a walking velocity as near as possible to the velocity of the non-delay case and get the lowest error between the follower equivalent velocity ze and the follower real velocity zr. Finally, from Figure 7, it can be seen that σr approaches σe to ensure stable behavior. Part (a) of Figure 7 shows that if h↑⇒σe↑ therefore λequivalent decreases. However, without a proper control λhybrid will be greater than λequivalent, which does not guarantee the stability of the delayed bilateral teleoperation system. Part (b) of Figure 7 states that if h↑⇒σe↑ therefore λequivalent decreases and if h↑⇒T↑∧σr↑ therefore, λhybrid decreases in the same way as the equivalent, which allows to hold a stable behavior assuring bounded errors.Figure 7
Evolution of the behavior of the equivalent follower robot (continuous system) and the humanoid robot (hybrid system).


## 3.3.1. Definition

Definition 1.
Tmin:=minT/maxdi,Δp=Δpmin,σr=σrmin,Tmin,ε=εminλhybrid=1Definition 2.
σrmin:=minσr/maxdi,Δp=Δpmin,Tmin,ε=εminλhybrid=1Definition 3.
βe:=minβe/maxdi,σe=σemin,βeλequivalent=1 for the nonperturbed discrete-time version of formula (28):(29)Dη˙=kskgxmt−h1−ze−σez˙e,ze=η−βez˙e,where λhybrid are the eigenvalues of the hybrid system computed from the Jacobian’s eigenvalues of the Poincaré map. λequivaalent are the eigenvalues of the discretized version of (29), di=Dksβr, σe is the linear damping coefficient injected by the equivalent follower controller, and ze is the velocity of the equivalent follower. The values σrmin,εmin,σemin are positive values near to zero while 0<Δpmin<1 establish the range of the reset map from Δp=Δpmin,1.
Due to the effect caused by the increasing of equivalent damping (equivalent follower controller),λequivalent decreases. It is possible to get a similar λhybrid by applying different combinations of real follower damping and walking cycle time, as shown in Figure 7. Now, we propose to online set σr and T through the following optimization criteria, visually illustrated in Figure 7:(30)σr,T,ε=minkσσr−σe+kwtT+kεεsuch that:maxci,σr,T,ελhybrid=maxdi,σeh1,h2λequivalent,where kσ,kwt,kε are the weights of the functional. This optimization cost function aims to find a trade-off between a walking velocity as near as possible to the velocity of the non-delay case and get the lowest error between the follower equivalent velocity ze and the follower real velocity zr. Finally, from Figure 7, it can be seen that σr approaches σe to ensure stable behavior. Part (a) of Figure 7 shows that if h↑⇒σe↑ therefore λequivalent decreases. However, without a proper control λhybrid will be greater than λequivalent, which does not guarantee the stability of the delayed bilateral teleoperation system. Part (b) of Figure 7 states that if h↑⇒σe↑ therefore λequivalent decreases and if h↑⇒T↑∧σr↑ therefore, λhybrid decreases in the same way as the equivalent, which allows to hold a stable behavior assuring bounded errors.Figure 7
Evolution of the behavior of the equivalent follower robot (continuous system) and the humanoid robot (hybrid system).


## 3.4. Stable Bilateral Controller

In this work, aP + d such as the controller is used to analyze what conditions must be held to keep the errors ec=kgxm−ze bounded. The proposed control strategy adds the following force feedback law for fm:(31)fm=−kmkgxm−zet−h2+gmxm−αmx˙m,where km is a gain matrix, kg linearly maps the leader’s position to a reference forward and turn speeds, and αm is the linear damping coefficient injected into the leader. The only difference with a traditional P + d scheme is that the equivalent follower velocity is used over the real one. Besides, the control actions applied to the real follower are given by formula (15).


## 3.5. Stability Analysis

Here, the stability of the proposed teleoperation system is analyzed and the process of calibrating the control parameters is explained. Iffm (31) is included in the leader (1) and with the equivalent robot (29), the following closed-loop dynamics are obtained:(32)x¨m=Mm−1−Cmx˙m−kmkgxm−zet−h2−αmx˙m+fh,(33)η˙=D−1kskgxmt−h1−ze−σez˙e−kzfv−w˜l.According to formulas (32) and (33), and taking also into account (5), it can be seen that x¨m=0,η˙=0, if the vector xt:=x˙mzez˙eecT remains zero for a sufficiently high time interval, considering fh, fv, and w˜l null and whose analysis is associated with Lyapunov–Krasovskii’s concept of stability, where there is a continuous of states that we represent using integrals with finite limits.A positive definite functionVxt=V1+V2+V3+V4 is given as follows:(34)V1=12l1z˙eTDz˙e,(35)V2=12l2x˙mTMmxmx˙m+l2Eh,(36)V3=12l3ecTec+12l4zeTze,(37)V4=l5∫−h¯20∫t+θtz˙eξTz˙eξdξdθ+l6∫−h¯10∫t+θtx˙mξTx˙mξdξdθ.Following the procedure described in reference [11] along the trajectories of the closed-loop system, considering the follower robot, leader dynamics, time delay, and human operator and environment forces, is bounded by:(38)V˙≤−λmx˙mTx˙m−λseze.Tz˙e−ρez˙e,where(39)ρe=w˜l+kzf¯v.Ifαm and σe are sufficiently high to comply as follows:(40)λm=αmkm−h1−h24>0,λse=σe+Dks−14h1−h2>0,then, the variables x˙m,z˙e∈ℒ∞. If formula (38) is integrated in time, the following is verified:(41)Vt−V0≤−λmx˙m22−λsez˙e22−∫0tz˙eϵTρeϵdϵ.


### 3.5.1. Remark

The term∫0tz˙eϵTρeϵdϵ is bounded since x˙m,z˙e∈ℒ∞ and using the proof of Section 3.3, it is concluded that the estimation errors are limited by W^∈ℒ∞. Therefore, from formula (41), we infer that Vt will be bounded for all t, and therefore, x˙m,z˙e∈ℒ2 and kgxm−ze,ze,xm,∈ℒ∞. If fv,fh,e˜l remain zero, and following an analysis similar to reference [16] based on Barbalat’s Lemma, it is possible to deduce that x˙m and z˙e will tend to zero as t⟶∞.


## 3.5.1. Remark

The term∫0tz˙eϵTρeϵdϵ is bounded since x˙m,z˙e∈ℒ∞ and using the proof of Section 3.3, it is concluded that the estimation errors are limited by W^∈ℒ∞. Therefore, from formula (41), we infer that Vt will be bounded for all t, and therefore, x˙m,z˙e∈ℒ2 and kgxm−ze,ze,xm,∈ℒ∞. If fv,fh,e˜l remain zero, and following an analysis similar to reference [16] based on Barbalat’s Lemma, it is possible to deduce that x˙m and z˙e will tend to zero as t⟶∞.


## 3.6. Cascade Calibration Process

A three-step calibration process is used to set the controller parameters, linking the stability result obtained by using the equivalent follower with the equivalent follower controller, as well as the relation between the equivalent dynamics and the hybrid dynamics as follows:(i)
V.1 Offline calibration step, the inverse dynamics compensation is calibrated, including their gains that depend on the used method, in particular, we apply an adaptive inverse dynamic compensation [20]. Furthermore, the reference parameters D and the gains of the bilateral controller (ks, km) without delay are adjusted while the observer parameter βr is set empirically. Finally, βe is set based on Definition 1.(ii)
V.2 Online step, the leader dampingαm and equivalent follower damping σe are calculated from formula (40).(iii)
V.3 Online step, the values of walking cycle timeT and real damping σr are set through the expression (30).


## 4. Simulations Results

The simulation was performed using MATLAB/Simulink, where the equations (20)–(23) are simulated using the parameters set in Section 3.2, and the control parameters are adjusted following the guidelines given in subsection F of Section 3. Using those parameters, two kinds of simulations are performed:First, with constant parameters, we performed 500 simulations using fixed parameters, with real damping from 0 to 10, and a walking cycle time from 0.5 to 5 seconds, taking random values under those intervals for 4 different cases of time delays. Keep fixed the first step and the second is calculated using formula (40).Then, with variable parameters, we performed 500 simulations using real damping and a walking cycle time according to the proposal (30), for the same 4 different cases of time delays.Figures8(a) and 8(b) show that the average coordination errors decrease when using the proposal that varies the damping and the walking cycle time (VP scheme), depending on the time delay, unlike the CP scheme, which uses a constant value of walk cycle time and damping for all delays. Therefore, in the obtained figures, the advantages achieved by the proposed scheme can be evidenced.Figure 8
Coordination errors of the teleoperation system (a) for forward speed and (b) for turn speed.
(a)(b)


## 5. Experiments

This section presents the experimental results for the control scheme proposed in Section3, where the forward and turn speeds of an NAO humanoid robot are remotely teleoperated. At the same time, the operator can perceive force feedback proportional to the coordination error between the leader and the equivalent follower robot and receives video feedback in order to drive the robot through the desired path. A fixed camera was mounted in the laboratory to take the workspace of the remote robot. The image is transmitted to the human operator who receives such visual feedback through a standard display while force feedback is felt simultaneously. The experimental scheme consists of a human operator using a Novint Falcon device with force feedback in Argentina to drive an NAO robot located in Germany. The human operator must evade the obstacle by left and right sides alternately for 10 tests. The parameters are calibrated following the guidelines given in V.1, V.2, and V.3. Also, a moving average filter of 1 second is used to estimate the current time delay as shown in Figure 9 and from then calculate the walking cycle time and real damping. Next, lateral stability control is performed to maintain lateral balance through walking, inducing a lateral oscillation by using only the roll joints in the hips and ankles [26] and merging it with the forward movement. A sequence of images of such a test is presented in Figure 10, while the video of the whole experiment can be seen at https://youtu.be/zeAZTzrCks4.Figure 9
Total internet time delayh1+h2 of the experiment.Figure 10
Test workspace and task to be performed.Figures11(a) and 11(b) show the evolution in time of the NAO’s forward and turn speeds during the execution of the test according to the reference values sent by the human operator. It is important to remark in some time intervals acts the fictitious force due to the near of obstacles, as shown in the rectangle of Figure 11, rising the coordination errors in such moments. Force feedback fm=fmvhipfmδT is also shown, during the time interval (20, 60) seconds, the virtual force acts on the robot slowing down its speed, which increases the coordination error. Thus, the force feedback affects the operator’s decision by letting him know, through a force opposite to the direction of the robot, that a bigger-than-normal synchronism error exists during the experiment, and this pushes the operator’s hand to a lower forward speed in order to decrease the error and get to better control of the robot. The speed tracking curve of the follower and the forward speed according to the values received in the robot, we can see in the following Figure 12. Where it is shown that the actual robot motion follows the human commands with a bounded error throughout the duration of the experiment. The fictitious force acts to avoid collisions at cost of rising the coordination errors. In terms of stability, the coordination error between the follower and leader is bounded by the proposed cascade controller. Finally, Table 2 highlights some important average data results from obstacle evasion for both left and right, where it can be remarked that the average delay of the tests is close to the average transatlantic delay (Buenos Aires-Munich) of 244 ms.Figure 11
Coordination error of the forward and turn speed and force feedback: (a) linear velocity command and linear velocity NAO and (b) turn velocity command and turn velocity NAO.
(a)(b)Figure 12
Forward velocity and Turn angle references generated by the human operator and the signals provided by the humanoid robot: (a) linear velocity command and linear velocity NAO and (b) turning angle command and turning angle NAO.
(a)(b)Table 2
Average data results.
VariableLeft and right side evasionMin delay261.85 msMax delay643.6 msAverage delay360.5 msTask completion time176.05 sForward speed average error1.54 cm/sTurning speed average error1.2 degrees/s


## 6. Conclusions and Discussion

In this article, a control scheme is presented for the teleoperation of the forward and turning speeds in the presence of time delays. The theoretical analysis and the definition proposed for the equivalent follower and the use of an equivalent follower controller allow us to perform a complete stability analysis where in contrast to the mobile robots and manipulators teleoperated cases, as a result, we obtain that it is not sufficient to increase only the damping as the time delay rises for holding stability. For teleoperated biped robots, there are multiple combinations of the real follower damping and the walking cycle time that can be used to assure bounded errors, whose online values are obtained by means of an optimization criterion that weighs the walking velocity and the errors between the real follower velocity and the equivalent follower velocity. It is necessary to emphasize that, due to the practical limitations of the NAO robot used in the present work, it becomes difficult to perceive a clear difference between constant and variable optimization parameters (damping and walking cycle time). The robot has a short walking speed range, where very low speeds produce falls or loss of balance.Adequate use of force feedback allows for decreasing those errors [11, 27]. We point out the work is focused on reaching a controller to ensure bounded errors despite these time delays, and not on optimizing the perception level. A quantitative analysis of the perception of the human operator in front of delays is outside of the scope of this work. However, the use of equivalent follower velocity in the force feedback algorithm instead of the real follower velocity avoids getting worst the level of perception of the human operator about the robot’s walk, feeling in his hand a continuous force (avoiding the discontinuities caused by the impacts) that pushes the human hand towards fewer coordination errors. Finally, it is also important to remark that although the NAO robot has practical limitations, it was used to get a simple test, while we point out that the conceptual contributions here presented can be applied to other bipedal and humanoid robots.



---

*Source: 3197743-2023-02-08.xml*