# Energy-Optimized Consensus Formation Control for the Time-Delayed Bilateral Teleoperation System of UAVs

Authors: Hui-yu Sun; Guang-ming Song; Zhong Wei; Ying Zhang
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

This paper proposes an energy-optimized consensus formation scheme for the time-delayed bilateral teleoperation system of multiple unmanned aerial vehicles (UAVs) in the obstructed environment. To deal with the asymmetric time-varying delays in aerial teleoperation, the local damping is independently distributed on both sides to enforce consensus formation and force tracking of the master haptic device and the slave UAVs. The stability of the time-delayed aerial teleoperation system is analyzed by the Lyapunov function. In addition, a flux-conserved force field is incorporated into the aerial teleoperation system to guarantee a collision-free consensus formation in the obstructed environment. Moreover, to reduce the communication complexity and energy dissipation of the formation, a top-down strategy of 3D optimal persistent graph is first proposed to optimize the formation topology. Under the optimized topology with environmental constraints, communication complexity and energy dissipation can be minimized while the rigid formation can be maintained and transformed persistently in the obstructed environment. Finally, the human-in-the-loop simulations are performed to validate the effectiveness of the proposed scheme.


---

## Body

## 1. Introduction

Groups of unmanned aerial vehicles (UAVs) rather than a single UAV have proven to be very effective in solving complex tasks like surveillance, exploration, and search and rescue [1, 2]. Nevertheless, when the task is extremely complex and high-level decisions are required online, complete autonomy of the agents is far from being reached and the human intervention is still necessary. Therefore, the bilateral teleoperation of a semiautonomous group seems a very promising solution [3]. In this way, the operator can manipulate a group of UAVs by the master device and receive the feedback force on the status of the slave group.The research of the bilateral teleoperation system has been studied for several decades [4, 5]. One critical problem is the time delays that existed in communication. Due to the long distance, limited bandwidth, and packet loss, the time delays are unavoidable and they will degrade the system performance or even cause instability [6]. Without proper control strategy, small time delay (e.g., tens of milliseconds) can destabilize the bilateral teleoperation system [7], thus leading to high risk to the human operator and the environment. A large number of researchers have proposed methods to make a trade-off between stability and transparency as reported in literature [8, 9]. For the constant or time-varying delays, the single-master-single-slave configuration of the teleoperation system has been studied for several decades. Anderson and Spong [7] firstly proposed the scattering theory to passify the system with constant time delays. Niemeyer and Slotine [10] developed the concept of the wave variables introduced from the scattering theory. Both notions have been virtually the only way to enforce passivity of the delayed bilateral teleoperation system.However, the scattering method brings about the position drifting problem that will degrade the stability of the system [11]. Therefore, Nuño et al. [12, 13] put forward the damping intervention method to overcome the position drifting problem. The large damping injections can be independently distributed on both sides. They not only guarantee the passivity of the teleoperation system but deal with the position drifting problem in the scattering system. Compared with the single slave configuration, there are fewer literatures about the single-master-multi-UAV configuration due to its complexities of model and coupling. Franchi et al. [14] proposed the large damping injections that were distributed on the master and each slave UAV to ensure the passivity of the teleoperation system, and the stability conditions of the multi-UAV configuration were given, but the time delays were not considered. Rodríguez-Seda et al. [15] addressed the task of remotely controlling a formation of n-degree-of-freedom (DOF) slave agents coupled bilaterally through constant time-delayed communication channels to a single n-DOF master robot. The proportional-derivative controller was designed to enforce formation control, master-to-slave position coverage, and static force reflection. However, the time delays were constant. Also, the master-slave robots were isomorphic and had the same kinematics and dynamics. Only the position convergence could be achieved under the proposed controller. The mismatches of master-slave kinematics and dynamics were not considered. Chawda and O’Malley [16] studied the time-varying delays in the multislave system while the consensus formation and the interaction with the environment were not considered. Slawiñski and Mut [17] put emphasis on the consensus problem for the teleoperation system over communication networks with time-varying delays by adding local damping at both sites, but the master and slave robots have the same dynamical models. It is not suitable for UAVs with underactuated dynamics and unbounded workspace. The local damping action on the master was not discussed. To our best knowledge, the consensus formation control for the multi-UAV bilateral teleoperation system in the obstructed environment while considering the asymmetric time-varying delays has as yet little researches.For the cooperative control of the multi-UAV bilateral teleoperation system, another critical problem is the consensus formation control in the obstructed environment. In the obstructed environment, the researches about collision avoidance have been active for many years. Usually, the collision avoidance approaches can be classified into behavior-based [18, 19] and potential field approaches [20, 21]. In the aerial teleoperation system, a flux-conserved force field that is inspired by the electronic field is incorporated to avoid the obstacles and escape the local minima. As pointed by Olfati-Saber et al. [22], besides collision avoidance, the connectivity maintenance is the fundamental rule for the coordination of multiagent systems in practice. In the multiagent system, moving in formation is a cooperative task and requires collaboration of every agent in the formation. The fundamental of the consensus formation control is that the relative position and velocity of neighbor vehicles are all in agreement. In these literatures, the relative-position-based approach is applied to coordinate the neighbor vehicles, in which each vehicle is required to communicate with its neighbors [23–25]. However, some interactions between the UAVs are not necessary and they will make the communication much more complex. If consensus can be reached with less neighbor information, the communication links and energy dissipation of the slave formation can be decreased. Jakovetic et al. [26] designed the weights in consensus algorithms for spatially correlated random topologies to decrease the topology complexity of graphs while maintaining their shapes. Yan et al. [27] firstly used the min-weighted rigid graph to describe the topology relation of slave robots in the teleoperation system. Correspondingly, the communication complexity and energy consumption in slave sire can be decreased. Lin et al. [28] concentrated on the fundamental coordination problem that requires a network of agents to achieve a specific but arbitrary formation shape. The complex Laplacian was introduced to address the problems of which formation shapes specified by interagent relative positions can be formed and how they can be achieved with distributed control ensuring global stability. However, these emphases were only placed on formation in 2D space or static formation in 3D space.Inspired by these issues, an energy-optimized consensus formation scheme for the asymmetric time-varying bilateral teleoperation system is studied based on UAV dynamics. The main contributions of this paper mainly lie in the following: (1) coordination control. A passive consensus formation controller is proposed to keep the system be closed-loop stable with the asymmetric time-varying delays while enforcing formation control, master-to-slave velocity convergence, and force reflection; (2) topology optimization. A top-down strategy of 3D optimal persistent graph is first used to describe the topology relationships of slave robots in the obstructed environment. Under the optimized topology with environmental constraints, the communication complexity and energy dissipation can be minimized while the rigid formation can be maintained and transformed persistently in the obstructed environment; and (3) collision avoidance control. Inspired by electromagnetics, each obstacle is regarded as the electronic objects. A new flux-conserved force field is proposed to keep the consensus formation away from the obstacles and escape the local minima.The rest of the paper is organized as follows. Section2 introduces a brief summary of the relevant results on graph theory, rigid graph, and system dynamics. Section 3 proposes the time-delayed bilateral aerial teleoperation system with consensus formation in the obstructed environment. Section 4 proposes a top-down strategy of 3D optimal persistent graph to minimize the communication complexity and energy dissipation of the formation. The human-in-the-loop simulation results and discussions that evaluate the effectiveness of the proposed framework are performed in Section 5. Finally, the conclusions and future work are given in Section 6.


## 2. Preliminaries

### 2.1. Graph Theory

Each UAV in the formation can be regarded as a vertex, and the topology of the slave UAVs is conveniently described as an undirected graph. Some basic concepts of graph theory are introduced in [29].LetG = (V, E, A) be a weighted graph, where V = {1, 2, …, n} is the set of vertexes, E⊆V×V is the set of edges, and A = [aij] is the weighted adjacency matrix. The adjacency elements associated with the edges are positive, that is, i,j∈V⇔aij>0. Moreover, the elements aii = 0 for all i∈V. The set of neighbors of vertex i is denoted by Ni=j∈V:i,j∈E. The Laplacian of graph G is denoted by L = [lij] with lii=∑j∈Niwij, lij=−wij where wij > 0 if j ∈ Ni. In is an identity matrix of size n, 1n is a column vector of size n with all elements equal to one. The L∞ norm is f∞=supt≥0ft. The L2 norm is f22=∫0tft2dt.


### 2.2. Rigid Graph

As shown in [30], let qi(t) be the trajectory of vertex i in the formation. A graph is said to be rigid if and only if the distances between every pair ||qi(t) − qj(t)|| are constant for all i,j∈E. Or else, it is called flexible.For a graphG⊆Rcc∈2,3, the rigidity matrix M is defined as
(1)Mn×dn=⋮⋱⋮⋮⋮⋯⋮⋮⋮⋱⋮0⋯0qiT−qjT0⋯0qjT−qiT0⋯0⋮⋱⋮⋮⋮⋯⋮⋮⋮⋱⋮,where each row 0⋯0qiT−qjT0⋯0qjT−qiT0⋯0 corresponds to an edge i,j∈E and the columns corresponds to the coordinates of the vertexes.Lemma 1 [30, 31].
LetM be the rigidity matrix of a framework of n vertexes in Rc. A framework G = (V, E, A) with n (n ≥ 2) vertexes in Rc is infinitesimally rigid if and only if rank (M) = cn − c(c + 1)/2.Lemma 2 [30, 31].
An infinitesimally rigid frameworkG = (V, E, A) with n (n ≥ 2) vertexes in Rc is minimally rigid if and only of it has cn − c(c + 1)/2 edges.Lemma 3 [30, 31].
If every edge of the frameworkG = (V, E, A) is weighted by its length, a min-weighted rigid graph is the minimally rigid that has the minimally weighted sum in all infinitesimally rigid graphs.Lemma 4 [30, 31].
If the frameworkG = (V, E, A) is the optimal persistent graph if and only if the out-degrees of any vertex in the min-weighted rigid graph are no more than 2.Remark 1.
The optimal persistent graph mainly has two important features. First, it is a directed rigid graph with the smallest communication complexity and the links. Then the weighted sum in all minimally rigid graphs is also the smallest. Based on these features, the optimal persistent graph can be adopted to optimize the communication of the slave UAVs. Examples of the flexible graph, rigid graph, min-weighted graph, and optimal persistent graph in 2D space are shown in Figure1.Figure 1
Examples of graphs in 2D space. (a) Flexible graph. (b) Rigid graph. (c) Min-weighted rigid graph. (d) Optimal persistence graph. The four agents in the plane are (0, 0), (3, 0), (4, 4), and (0, 5).
(a)
(b)
(c)
(d)


### 2.3. System Dynamics

The master device of the bilateral aerial teleoperation system is shown in Figure2(a). The operator is incorporated into the overall framework and remotely manipulates the formation by the master device. The device is a fully actuated mechanical system described by the following Euler–Lagrange equation [32]:
(2)Mmqmq¨m+Cmqm,q̇mq̇m=fm+fh,where qm,q̇m,q¨m∈R3 is the joint position, velocity, and acceleration. Mmqm∈R3×3 is the positive definite inertia matrix. Cmqm,q̇m∈R3×3 is the Coriolis and centripetal matrix. For the sake of simplicity, the enclosed parameter(s) of Mmqm,Cmqm,q̇m are attached in Appendix A. fh∈R3 is the human force applied on the device with the compensated gravitational force. fm∈R3 is the haptic force feedback to the device. In order to keep equations as clear as possible, the sign (t) of all variables at time t are omitted (e.g., qm≡qmt), except for those which are time delayed (e.g., qmt−Tm1≡qmt−Tm1t for a variable time delay).Figure 2
The devices in the bilateral teleoperation system. (a) The physical map and kinematic model diagram of the master device. (b) The slave UAV.
(a)
(b)As shown in Figure2(b), the slave ith UAV is a rigid body that is actuated by a thrust force generated by four rotors on the endpoints of the X-shaped frame. Let ϖi = [xi, yi, zi, φi, θi, ψi] be the generalized coordinates where qi = [xi, yi, zi]T ∈ R3 denotes the absolute position of the ith UAV in the inertial frame and Θi = [φi, θi, ψi]T ∈ R3 denotes the three Euler angles representing the roll, pitch, and yaw, respectively. The dynamic model of the ith UAV can be derived using the Euler–Lagrange equations that partitioned into the translational and rotational coordinates [33].
(3)miq¨i=fi−00mig,(4)JiΘ¨i+ddtJiΘ̇i−12∂∂tΘ̇iJiΘ̇i=τi,where mi is the mass. Ji∈R3 is the moment of inertia matrix. fi,τi∈R6 is the external generalized force applied to the ith UAV. For the sake of simplicity, the derivations of the Euler–Lagrange equations are attached in Appendix B.The Coriolis-centripetal vector is defined by(5)V¯Θi,Θ̇i=J̇iΘ̇i−12∂∂tΘ̇iJiΘ̇i=J̇i−12∂∂tΘ̇iJiΘ̇i=CΘi,Θ̇iΘ̇i,where CΘi,Θ̇i∈R3×3 is the Coriolis matrix.Then (4) can be rewritten as
(6)JiΘ¨i+CΘi,Θ̇iΘ̇i=τi.The dynamics of theith UAV can be yielded as
(7)miq¨i=fi−mige3,JiΘ¨i+CΘi,Θ̇iΘ̇i=τi.In the bilateral teleoperation of UAVs, UAVs suffered the external generalized force caused by the interaction of neighbor UAVs and UAVs with the environment, that is,fi,τi=fis+fie,τis+τie. Here, fis is the translational force applied on the ith UAV with the compensated gravitational force. τis is the roll, pitch, and yaw torques applied on the ith UAV. fie,τie is the force/torque caused by the interaction of the ith UAV with the environment. Finally, the dynamics of the ith UAV (i ∈ V) can be defined as
(8)miq¨i=fis+fie,JiΘ¨i+CΘi,Θ̇iΘ̇i=τis+τie.Besides, letvi ∈ R3 be the linear velocity and the variable ρi ∈ R3 is defined as
(9)v̇i=ρi+γρ̇i,where γ → 0+, v̇i∈R3 is the acceleration of the ith UAV, and ρi ∈ R3 represents the acceleration at an infinitesimal time instant before t so it can be assumed that ρi≈v̇i.To simplify the stability analysis of the time-delayed teleoperation system, the following properties, assumptions, and lemmas will be used in this paper [34, 35].Property 1.
Mmqm−2Cmqm,q̇m is skew symmetric.Property 2.
∃λα,λβ∈R+ such that λαI≤Mm≤λβI.Property 3.
∃kα∈R+ such that Cmqm,q̇mq̇m≤kαq̇m.Assumption 1.
The human operator and the environment are passive fori ∈ V, that is,
(10)∫0Tq̇mTtfhtdt≥0,−∫0TviTtfietdt≥0.Assumption 2.
The time-varying delaysTm1(t), T1m(t), Tij(t), and Tji(t) are bounded, that is, 0 ≤ Tm1(t) ≤ T¯m1, 0 ≤ T1m(t) ≤ T¯1m, 0 ≤ Tij(t) ≤ T¯ij, and 0 ≤ Tji(t) ≤ T¯ji, where T¯m1 is the maximum time-varying delay from the master to the leader UAV, T¯1m is the maximum time-varying delay from the leader UAV to the master, T¯ij is the maximum time-varying delay from the ith UAV to the jth UAV, and T¯ji is the maximum time-varying delay from the jth UAV to the ith UAV.Lemma 5.
For any vector functiona(∙) and b(∙) and any time-varying delay 0≤Tt≤T¯, there exists a positive-definite matrix Γ such that the following inequality holds:
(11)−2aTt∫t−Tttbξdξ−∫t−TttbTξΓbξdξ≤TaTtΓ−1at≤T¯aTtΓ−1at.Lemma 6.
For any vector functionx(∙) and y(∙) and any time-varying delay 0≤Tt≤∞, there exists a constant α such that the following inequality holds:
(12)−2∫0txTσ∫0−T0yσ+θdθdσ≤αx22+T2αy22,where ·22 is the L2 norm of the signal.


## 2.1. Graph Theory

Each UAV in the formation can be regarded as a vertex, and the topology of the slave UAVs is conveniently described as an undirected graph. Some basic concepts of graph theory are introduced in [29].LetG = (V, E, A) be a weighted graph, where V = {1, 2, …, n} is the set of vertexes, E⊆V×V is the set of edges, and A = [aij] is the weighted adjacency matrix. The adjacency elements associated with the edges are positive, that is, i,j∈V⇔aij>0. Moreover, the elements aii = 0 for all i∈V. The set of neighbors of vertex i is denoted by Ni=j∈V:i,j∈E. The Laplacian of graph G is denoted by L = [lij] with lii=∑j∈Niwij, lij=−wij where wij > 0 if j ∈ Ni. In is an identity matrix of size n, 1n is a column vector of size n with all elements equal to one. The L∞ norm is f∞=supt≥0ft. The L2 norm is f22=∫0tft2dt.


## 2.2. Rigid Graph

As shown in [30], let qi(t) be the trajectory of vertex i in the formation. A graph is said to be rigid if and only if the distances between every pair ||qi(t) − qj(t)|| are constant for all i,j∈E. Or else, it is called flexible.For a graphG⊆Rcc∈2,3, the rigidity matrix M is defined as
(1)Mn×dn=⋮⋱⋮⋮⋮⋯⋮⋮⋮⋱⋮0⋯0qiT−qjT0⋯0qjT−qiT0⋯0⋮⋱⋮⋮⋮⋯⋮⋮⋮⋱⋮,where each row 0⋯0qiT−qjT0⋯0qjT−qiT0⋯0 corresponds to an edge i,j∈E and the columns corresponds to the coordinates of the vertexes.Lemma 1 [30, 31].
LetM be the rigidity matrix of a framework of n vertexes in Rc. A framework G = (V, E, A) with n (n ≥ 2) vertexes in Rc is infinitesimally rigid if and only if rank (M) = cn − c(c + 1)/2.Lemma 2 [30, 31].
An infinitesimally rigid frameworkG = (V, E, A) with n (n ≥ 2) vertexes in Rc is minimally rigid if and only of it has cn − c(c + 1)/2 edges.Lemma 3 [30, 31].
If every edge of the frameworkG = (V, E, A) is weighted by its length, a min-weighted rigid graph is the minimally rigid that has the minimally weighted sum in all infinitesimally rigid graphs.Lemma 4 [30, 31].
If the frameworkG = (V, E, A) is the optimal persistent graph if and only if the out-degrees of any vertex in the min-weighted rigid graph are no more than 2.Remark 1.
The optimal persistent graph mainly has two important features. First, it is a directed rigid graph with the smallest communication complexity and the links. Then the weighted sum in all minimally rigid graphs is also the smallest. Based on these features, the optimal persistent graph can be adopted to optimize the communication of the slave UAVs. Examples of the flexible graph, rigid graph, min-weighted graph, and optimal persistent graph in 2D space are shown in Figure1.Figure 1
Examples of graphs in 2D space. (a) Flexible graph. (b) Rigid graph. (c) Min-weighted rigid graph. (d) Optimal persistence graph. The four agents in the plane are (0, 0), (3, 0), (4, 4), and (0, 5).
(a)
(b)
(c)
(d)


## 2.3. System Dynamics

The master device of the bilateral aerial teleoperation system is shown in Figure2(a). The operator is incorporated into the overall framework and remotely manipulates the formation by the master device. The device is a fully actuated mechanical system described by the following Euler–Lagrange equation [32]:
(2)Mmqmq¨m+Cmqm,q̇mq̇m=fm+fh,where qm,q̇m,q¨m∈R3 is the joint position, velocity, and acceleration. Mmqm∈R3×3 is the positive definite inertia matrix. Cmqm,q̇m∈R3×3 is the Coriolis and centripetal matrix. For the sake of simplicity, the enclosed parameter(s) of Mmqm,Cmqm,q̇m are attached in Appendix A. fh∈R3 is the human force applied on the device with the compensated gravitational force. fm∈R3 is the haptic force feedback to the device. In order to keep equations as clear as possible, the sign (t) of all variables at time t are omitted (e.g., qm≡qmt), except for those which are time delayed (e.g., qmt−Tm1≡qmt−Tm1t for a variable time delay).Figure 2
The devices in the bilateral teleoperation system. (a) The physical map and kinematic model diagram of the master device. (b) The slave UAV.
(a)
(b)As shown in Figure2(b), the slave ith UAV is a rigid body that is actuated by a thrust force generated by four rotors on the endpoints of the X-shaped frame. Let ϖi = [xi, yi, zi, φi, θi, ψi] be the generalized coordinates where qi = [xi, yi, zi]T ∈ R3 denotes the absolute position of the ith UAV in the inertial frame and Θi = [φi, θi, ψi]T ∈ R3 denotes the three Euler angles representing the roll, pitch, and yaw, respectively. The dynamic model of the ith UAV can be derived using the Euler–Lagrange equations that partitioned into the translational and rotational coordinates [33].
(3)miq¨i=fi−00mig,(4)JiΘ¨i+ddtJiΘ̇i−12∂∂tΘ̇iJiΘ̇i=τi,where mi is the mass. Ji∈R3 is the moment of inertia matrix. fi,τi∈R6 is the external generalized force applied to the ith UAV. For the sake of simplicity, the derivations of the Euler–Lagrange equations are attached in Appendix B.The Coriolis-centripetal vector is defined by(5)V¯Θi,Θ̇i=J̇iΘ̇i−12∂∂tΘ̇iJiΘ̇i=J̇i−12∂∂tΘ̇iJiΘ̇i=CΘi,Θ̇iΘ̇i,where CΘi,Θ̇i∈R3×3 is the Coriolis matrix.Then (4) can be rewritten as
(6)JiΘ¨i+CΘi,Θ̇iΘ̇i=τi.The dynamics of theith UAV can be yielded as
(7)miq¨i=fi−mige3,JiΘ¨i+CΘi,Θ̇iΘ̇i=τi.In the bilateral teleoperation of UAVs, UAVs suffered the external generalized force caused by the interaction of neighbor UAVs and UAVs with the environment, that is,fi,τi=fis+fie,τis+τie. Here, fis is the translational force applied on the ith UAV with the compensated gravitational force. τis is the roll, pitch, and yaw torques applied on the ith UAV. fie,τie is the force/torque caused by the interaction of the ith UAV with the environment. Finally, the dynamics of the ith UAV (i ∈ V) can be defined as
(8)miq¨i=fis+fie,JiΘ¨i+CΘi,Θ̇iΘ̇i=τis+τie.Besides, letvi ∈ R3 be the linear velocity and the variable ρi ∈ R3 is defined as
(9)v̇i=ρi+γρ̇i,where γ → 0+, v̇i∈R3 is the acceleration of the ith UAV, and ρi ∈ R3 represents the acceleration at an infinitesimal time instant before t so it can be assumed that ρi≈v̇i.To simplify the stability analysis of the time-delayed teleoperation system, the following properties, assumptions, and lemmas will be used in this paper [34, 35].Property 1.
Mmqm−2Cmqm,q̇m is skew symmetric.Property 2.
∃λα,λβ∈R+ such that λαI≤Mm≤λβI.Property 3.
∃kα∈R+ such that Cmqm,q̇mq̇m≤kαq̇m.Assumption 1.
The human operator and the environment are passive fori ∈ V, that is,
(10)∫0Tq̇mTtfhtdt≥0,−∫0TviTtfietdt≥0.Assumption 2.
The time-varying delaysTm1(t), T1m(t), Tij(t), and Tji(t) are bounded, that is, 0 ≤ Tm1(t) ≤ T¯m1, 0 ≤ T1m(t) ≤ T¯1m, 0 ≤ Tij(t) ≤ T¯ij, and 0 ≤ Tji(t) ≤ T¯ji, where T¯m1 is the maximum time-varying delay from the master to the leader UAV, T¯1m is the maximum time-varying delay from the leader UAV to the master, T¯ij is the maximum time-varying delay from the ith UAV to the jth UAV, and T¯ji is the maximum time-varying delay from the jth UAV to the ith UAV.Lemma 5.
For any vector functiona(∙) and b(∙) and any time-varying delay 0≤Tt≤T¯, there exists a positive-definite matrix Γ such that the following inequality holds:
(11)−2aTt∫t−Tttbξdξ−∫t−TttbTξΓbξdξ≤TaTtΓ−1at≤T¯aTtΓ−1at.Lemma 6.
For any vector functionx(∙) and y(∙) and any time-varying delay 0≤Tt≤∞, there exists a constant α such that the following inequality holds:
(12)−2∫0txTσ∫0−T0yσ+θdθdσ≤αx22+T2αy22,where ·22 is the L2 norm of the signal.


## 3. Bilateral Aerial Teleoperation with Consensus Formation

In the bilateral aerial teleoperation with consensus formation, UAV1 is selected as the leader and the others are the followers. The leader is controlled by the master haptic device and the followers are required to communicate with the neighbors to form a rigid formation. The movement of each UAV depends on the surrounding UAVs and obstacles. The teleoperation system is shown in Figure3.Figure 3
The overall bilateral teleoperation system.


### 3.1. Workspace Mapping

Due to the mismatches of the workspaces between the master haptic device and the leader UAV, the position of the master device is regarded as a velocity setpoint for the leader. In order to solve the dissimilarity, a position-velocity workspace mapping is defined as(13)rm=knqm,kn∈R+,where kn is a positive constant.The commanded velocity for the leader UAV is(14)rs=knqmt−Tm1,kn∈R+,where Tm1 is the time-varying delay from the master device to the leader UAV. rs ∈ R3 is the commanded velocity for the leader UAV.It is well known that the master device is passive with respect to the force-velocity couplingfm,q̇m, where the velocity of the master device is synchronized with the velocity of the slave device. However, in our setting, in order to consider the difference between the workspace of the master and that of the robots at the slave site, it is necessary to synchronize the position of the master device with the velocity of the leader UAV. Unfortunately, the master device is not passive with respect to the force-position coupling fm,rm. In order to couple the position of the master device to the velocity of the leader, it is possible to implement a local control loop on the master that makes it passive with respect to the coupling fm,rm [14].First, letr=q̇m+λqm,λ∈R+ be the local control loop with a storage function κm=1/2q̇mTMmq̇m; the coupling fm,r is passive with λ→0+, that is,
(15)rTfm≥κ̇m.Further, by introducing a scaling into this strategy, letrm=μr=μq̇m+μλqm,μ,λ∈R+. By properly choosing the parameters μ and λ, it is possible to neglect the contribution of q̇m and make the commanded velocity rm proportional to the master scaled position with rm≈knqm,kn=μλ.The new energy storage is chosen as(16)κ¯m=μκm.The time derivative ofκ¯m is
(17)κ¯̇m=μκ̇m≤μrTfm=rmTfm.Therefore, the master device is passive with respect to the couplingfm,rm.


### 3.2. Master-Slave Coupling with Local Damping

It is well known that the stability of the bilateral teleoperation system with time-varying delays can be guaranteed by adding large damping injections at both sides. The damping distributions are shown in Figure4. Therefore, the consensus formation controller is designed as
(18)fm=−kmknqm−v1t−T1m−αmq̇m,(19)fis=−biksvi−knqmt−Tm1−αsv̇i−kv∑j∈Nivi−vjt−Tij−kc∑j∈Niqi−qjt−Tij−dij−βiq̇i,where T1m is the time-varying delay from the leader UAV to the master, Tm1 is the time-varying delay from the master to the leader UAV, and Tij is the time-varying delay from the ith UAV (i ∈ V) to the jth UAV. km, ks, kv, and kc are the positive constants. αm, αs, and βi are the damping coefficients. bi = 1 when the 1th UAV is connected to the master, or else, bi = 0. dij is the desired relative distance between the ith UAV and the jth UAV. Ni is the set of neighbors of the ith UAV. Obviously, this controller is a proportional-derivative controller. The proportional terms are aimed to reduce the tracking error in velocity and relative distance, and the derivative terms are used to maintain the system stable with less overshoot. km and ks are the positive constants which represent the proportional gains to achieve velocity coordination between the master device and the leader UAV. kv and kc are the positive coefficients which represent the proportional gains to achieve velocity and relative position coordination between slave UAVs. αm, αs, and βi are the damping coefficients which can make the system stable with less oscillatory. The error decreases with increasing proportional gains km, ks, kv, and kc, and the system may become more oscillatory. However, the system becomes damped when the damping coefficients αm, αs, and βi are increased. The constraint relationships can be seen in Theorem 1.Figure 4
The distributions of local damping in the aerial teleoperation system. (a) The distribution of the damping between the master and the leader UAV. (b) The distribution of the damping among the slave UAVs.
(a)
(b)Due to the consensus formation controller in (18), the roll and pitch angles in the rotational coordinates will gradually converge to the constant values. Therefore, the yaw coordinates of the ith UAV can be yielded as
(20)Jiψψ¨i=τisψ+τieψ,where Jiψ is the moment of inertia matrix in yaw coordinates and τisψ is the control torque applied on the ith UAV in yaw coordinates. τieψ is the torque caused by the interaction of the ith UAV with the environment in yaw coordinates. In the bilateral teleoperation system of UAVs with consensus formation control, the yaw angle of the ith UAV only depends on the repulsive forces which are produced by the obstacles in the obstructed environment, that is,
(21)ψ¨i=Jiψ−1Rifiez,for each slave UAV ensures that ψi−ψj,ψi∈L∞ for all i ∈ V, j ∈ Ni.Theorem 1.
For the single-master-multi-UAV teleoperation system with the consensus formation controller in (18), if there exist positive-definite matrices P and Q and any positive constants δi,εi such that the following inequalities hold:
(22)−αmI+T¯m1Q+T¯1mkm24P−1<0,−m+αsknksI+T¯1mP+T¯m1km24Q−1<0,βi>δi+εikc2lii+12δi∑j∈NiwijTij2−kv,∀i∈V,j∈Ni,then the overall teleoperation system is stable and the variables knqm−v1t−T1m, qi−qjt−Tij−dij, vi−vjt−Tij, ψi−ψj, q̇m, vi, and v̇i are bounded.Proof 1.
The total control forcefis applied on each UAV mainly contains two parts. In the first part, the velocity of the leader is required to converge to the velocity of the master. In this way, the operator can manipulate the movement of the slave formation by the master haptic device. In the second part, the UAVs are required to keep a relative position consensus with the neighbors. The control force fis applied on each UAV is linearly mixed. The first part is to track the master and the latter part is to keep consensus with the neighbors.
Thus, the control forcefis is divided into two parts, that is,
(23)fis=fis#+fis∗,where
(24)fis#=−biksvi−knqmt−Tm1−αsv̇i,fis∗=−kv∑j∈Nivi−vjt−Tij−kc∑j∈Niqi−qjt−Tij−dij−βiq̇i.
Firstly, the relationship between the master and the leader UAV is investigated. The total forces∑i=1nfis# which are applied on the leader UAV are
(25)fv1−knqmt−Tm11s#=−ks−αsv̇1.
Because the haptic device only has three translational output DOFs of force feedback, the rotational part of the bilateral teleoperation is unilateral. Therefore, only the stability of the translational part of the system is analyzed. The following Lyapunov function is defined as(26)H#=12q̇TmMmq̇m+m1γkm2knksρ1Tρ1+km2knknqm−v1T⋅knqm−v1+∫−T¯m10∫t+θtq̇mξTQq̇mξdξdθ+∫−T¯1m0∫t+θtv̇1ξTPv1ξdξdθ,where Q and P are positive-definite matrices. Obviously, function H# is always greater than zero.
With (9), Property 1, and Lemma 4, the time derivative of H# is
(27)Ḣ#=q̇mTfm+q̇mT⋅fh+kmknksρ1Tf1s#+f1e−m1kmknksρ1Tρ1+kmknknqm−v1Tknq̇m−v̇1+T¯m1q̇mTQq̇m−∫t−T¯m1tq̇mξTQq̇mξdξ+T¯1mv̇1TPv̇1−∫t−T¯1mtv̇1ξTPv̇1ξdξ.
With (18), (23), and Assumption 2, Ḣ# can be rewritten as
(28)Ḣ#≤q̇mTkmv1t−T1mt−knqm−αmq̇m+q̇mT⋅fh+kmknksρ1Tksknqmt−Tm1t−v1−αsv̇1+f1e−m1kmknksρ1Tρ1+kmknknqm−v1Tknq̇m−v̇1+T¯m1q̇mTQq̇m−∫t−T¯m1tq̇mξTQq̇mξdξ+T¯1mv̇1TPv̇1−∫t−T¯1mtv̇1ξTPv̇1ξdξ.
Noticing thatv1t−T1m−v1=−∫t−T1mtv̇1ξdξ,qmt−Tm1−qm=−∫t−T¯m1tq̇mξdξ, one has
(29)Ḣ#≤−kmq̇mT∫t−T1mtv̇1ξdξ−kmq̇mTknqm−v1−αmq̇mTq̇m+q̇mT⋅fh−kmρ1T∫t−T¯m1tq̇mξdξ+kmknρ1Tknqm−v1−m1+αskmknksρ1Tρ1+kmknksρ1Tfe+kmknqm−v1Tq̇m−kmknknqm−v1Tρ1+T¯m1q̇mTQq̇m−∫t−T¯m1tq̇mξTQq̇mξdξ+T¯1mv̇1TPv̇1−∫t−T¯1mtv̇1ξTPv̇1ξdξ.
With Lemma4, the following inequalities are obtained:
(30)−kmq̇mT∫t−T1mtv̇1ξdξ−∫t−T1mtv̇1ξTPv̇1ξdξ≤T¯1mkm24q̇mTP−1q̇m,−kmρ1T∫t−Tm1tq̇mξdξ−∫t−Tm1tq̇mξTQq̇mξdξ≤T¯m1km24ρ1TQ−1ρ1.
With (26), (27), (28), (29), and (30), one has
(31)Ḣ#≤−αmq̇mTq̇m+q̇mTfh−m1+αsknksρ1Tρ1+kmknksρ1Tf1e+T¯m1km24ρ1TQ−1ρ1+T¯m1q̇mTQq̇m+T¯1mkm24q̇mTP−1q̇m+T¯1mv̇1TPv̇1≤q̇mT−αmI+T¯m1Q+T¯1mkm24P−1q̇m+ρ1T−m1+αsknksI+T¯1mP+T¯m1km24Q−1ρ1+q̇Tmfh+kmknksρ1Tf1e.
With Assumption1, the third term of (31) is passive with respect to the force-velocity coupling fh,q̇m.
Similarly, the force-velocity couplingf1e,v1 is passive,
(32)v1Tf1e≥κ̇s,κs=m12v1Tv1.
Integrating equationv̇1=ρ1+γρ̇1,γ∈R+, one has
(33)v1=∫0Tρ1dt+γρ1,γ∈R+.
Letσ=μv1=μ∫0Tρ1dt+μγρ1,μ,γ∈R+. When the parameters μ and λ are properly selected, variable σ is proportional to the variable ρ1 with σ≈μγρ1.
Further, the new energy storage is selected as(34)κ¯s=μκs.
The time derivative ofκ¯s is
(35)κ¯̇s=μκ̇s≤μv1Tf1e=σTf1e≈μγρ1Tf1e.
Therefore, the fourth term of (31) is passive with respect to coupling f1e,ρ1.
In conclusion, to ensure the stability of the master and the leader,Ḣ# should be negative which is equivalent to the following inequalities:
(36)−αmI+T¯m1Q+T¯1mkm24P−1<0−m+αsknksI+T¯1mP+T¯m1km24Q−1<0.Setting αm, αs in (36) for the master and the leader UAVs implies that knqm−v1t−T1m,q̇m,v̇1∈L∞.
For the follower UAVs, the following Lyapunov function for theith UAV (i ∈ V) is constructed as
(37)Hi∗=12ṗiTMi−1pi−∫0tvifiedt,where pi∈R3 and Mi∈R3×3 are the momentum and the inertia matrix, respectively, of the ith UAV. With Assumption 1, Hi∗ is greater than zero.
With (23), the time derivation of Hi∗ is
(38)Ḣi∗=viTfis∗+fie−viTfie=−βivi2−kv∑j∈Nivi−vjt−Tij−kcviT∑j∈Niqi−qjt−Tij−dij=−βivi2−kvviT∑j∈Nivi−vj−kvviT∑j∈Nivj−vjt−Tij−kcviT∑j∈Niqi−qj−dij−kcviT∑j∈Niqj−qjt−Tij.
Considering the following total scaled energy function,(39)H∗=12qTL⊗In−diagARq+∑i=1nHi∗kc,where q=q1,…,qn, v=v1,…,vn, R is the relative distance matrix, and ⊗ is the standard Kronecker product.
(40)Ḣ∗=q̇TL⊗In−diagARq+∑i=1n−βi+kvkcvi2−kvkcvj−vjt−Tij−viT∑j∈Niqi−qj−dij−viT∑j∈Niqj−qjt−Tij=−∑i=1nβi+kvkcvi2−kvkc∑i=1n∑j∈NiviTvj−vjt−Tij−∑i=1n∑j∈NiviTqj−qjt−Tij.
Noticing thatqj−qjt−Tij=∫−Tij0q̇jξdξ, one has
(41)Ḣ∗≤−∑i=1nβi+kvkcvi2−∑i=1nviT∑j∈Ni∫−Tij0q̇jξdξ−kvkc∑i=1nviT∑j∈Ni∫−Tij0v̇jξdξ.
By integrating the above equation from zero tot with Lemma 5, we obtain
(42)H∗−H0≤−∑i=1nβi+kvkcvi22+∑i=1n∑j∈Niδi+εi2vi22+Tij22δivj22+Tij22εiv̇j22,where δi, εi > 0. Using the fact that H∗ > 0,
(43)H0≥∑i=1n∑j∈Niβi+kvkclii−δi+εi2vi22−Tij22δivj22.
DefiningΓ=v122,…,vn22∈Rn and Ω∈Rn×n,
(44)Ωn×n=β1kcl11−δ12−T1222δ1⋯−T1n22δ1−T2122δ2β2kcl22−δ22⋯−T2n22δ2⋮⋮⋱⋮−Tn122δn−Tn222δn⋯βnkclnn−δn2.
ThenH(0) can be rewritten as H0≥1nTΩΓ. There exists ℓl≥0 for l∈V such that H0≥ℓ1v122+⋯+ℓnvn22 and vi,v̇i∈L2. Therefore, setting the damping injection βi,
(45)βi>δi+εikc2lii+12δi∑j∈NiwijTij2−kv,for any arbitrary δi,εi>0. Each slave UAV ensures that vi∈L2 and H∗∈L∞, which in turn implies that qi−qjt−Tij−dij,vi−vjt−Tij,vi∈L∞, for all i∈V,j∈Ni.
Combining (21), (36), and (45), one has knqm→v1t−T1m, qi−qjt−Tij→dij, vi→vjt−Tij, ψi→ψj. The proof is completed.


### 3.3. Flux-Based Collision Avoidance

In general, the artificial potential fields [15] are usually used to avoid the obstacles. The gradient of the potential around the obstacles is used as the repulsive force exerted on the vehicle and the potential disappears outside of the influence region.Here, a flux-conserved force field is proposed to keep the consensus formation away from the obstacles and escape the local minima. Inspired by electromagnetics, each obstacle can be regarded as the electronic objects. LetEiobs denote the force field and C denote any closed curves around the obstacle in the 3D space. The flux is defined as
(46)Wiobs=∮Cn⋅Eiobsds,where Wiobs is the flux of obstacle i for i ∈ {1, …, nobs} and n is the normal vector of the curve.Similar to the conservation law of electric flux in the electromagnetic field, the most important feature of the flux-conserved force field is that the flux of the force field is conserved for any curves. In Figure5, for any curve C surrounding obstacle i, the value of Wiobs remains the same. The repulsive force of the ith UAV from the obstacles is
(47)fie=kobs∑i=1nobs∂Wiobs∂qi−xoi,where qi ∈ R3 is the position of the ith UAV and xoi∈R3 is the position of obstacle i. kobs is a positive coefficient. For any curve surrounding an electric point, the value of Wiobs is always same.Figure 5
The flux of the force field is conserved for any curves.With the flux-conserved force field, it can be proved that there are no local minima in this force field. Proving can be demonstrated as follows. As shown in Figure6, assuming a local minimum, point P appears around obstacle i. The local minimum is a point where the sum of the force field is close to zero. In this figure, the value of the flux generated by the obstacle is defined as positive and the value of the flux generated by the target is defined as negative. Firstly, considering an imaginary little closed curve around the local minimum, point P and the force on the curve point toward the inside. Then because the force filed is flux conserved in the closed curve, the curve must have a target with a negative flux value. Therefore, if there is not a target in the closed curve, the local minimum around the obstacle must certainly not occur.Figure 6
No local minima in flux-conserved force field.


## 3.1. Workspace Mapping

Due to the mismatches of the workspaces between the master haptic device and the leader UAV, the position of the master device is regarded as a velocity setpoint for the leader. In order to solve the dissimilarity, a position-velocity workspace mapping is defined as(13)rm=knqm,kn∈R+,where kn is a positive constant.The commanded velocity for the leader UAV is(14)rs=knqmt−Tm1,kn∈R+,where Tm1 is the time-varying delay from the master device to the leader UAV. rs ∈ R3 is the commanded velocity for the leader UAV.It is well known that the master device is passive with respect to the force-velocity couplingfm,q̇m, where the velocity of the master device is synchronized with the velocity of the slave device. However, in our setting, in order to consider the difference between the workspace of the master and that of the robots at the slave site, it is necessary to synchronize the position of the master device with the velocity of the leader UAV. Unfortunately, the master device is not passive with respect to the force-position coupling fm,rm. In order to couple the position of the master device to the velocity of the leader, it is possible to implement a local control loop on the master that makes it passive with respect to the coupling fm,rm [14].First, letr=q̇m+λqm,λ∈R+ be the local control loop with a storage function κm=1/2q̇mTMmq̇m; the coupling fm,r is passive with λ→0+, that is,
(15)rTfm≥κ̇m.Further, by introducing a scaling into this strategy, letrm=μr=μq̇m+μλqm,μ,λ∈R+. By properly choosing the parameters μ and λ, it is possible to neglect the contribution of q̇m and make the commanded velocity rm proportional to the master scaled position with rm≈knqm,kn=μλ.The new energy storage is chosen as(16)κ¯m=μκm.The time derivative ofκ¯m is
(17)κ¯̇m=μκ̇m≤μrTfm=rmTfm.Therefore, the master device is passive with respect to the couplingfm,rm.


## 3.2. Master-Slave Coupling with Local Damping

It is well known that the stability of the bilateral teleoperation system with time-varying delays can be guaranteed by adding large damping injections at both sides. The damping distributions are shown in Figure4. Therefore, the consensus formation controller is designed as
(18)fm=−kmknqm−v1t−T1m−αmq̇m,(19)fis=−biksvi−knqmt−Tm1−αsv̇i−kv∑j∈Nivi−vjt−Tij−kc∑j∈Niqi−qjt−Tij−dij−βiq̇i,where T1m is the time-varying delay from the leader UAV to the master, Tm1 is the time-varying delay from the master to the leader UAV, and Tij is the time-varying delay from the ith UAV (i ∈ V) to the jth UAV. km, ks, kv, and kc are the positive constants. αm, αs, and βi are the damping coefficients. bi = 1 when the 1th UAV is connected to the master, or else, bi = 0. dij is the desired relative distance between the ith UAV and the jth UAV. Ni is the set of neighbors of the ith UAV. Obviously, this controller is a proportional-derivative controller. The proportional terms are aimed to reduce the tracking error in velocity and relative distance, and the derivative terms are used to maintain the system stable with less overshoot. km and ks are the positive constants which represent the proportional gains to achieve velocity coordination between the master device and the leader UAV. kv and kc are the positive coefficients which represent the proportional gains to achieve velocity and relative position coordination between slave UAVs. αm, αs, and βi are the damping coefficients which can make the system stable with less oscillatory. The error decreases with increasing proportional gains km, ks, kv, and kc, and the system may become more oscillatory. However, the system becomes damped when the damping coefficients αm, αs, and βi are increased. The constraint relationships can be seen in Theorem 1.Figure 4
The distributions of local damping in the aerial teleoperation system. (a) The distribution of the damping between the master and the leader UAV. (b) The distribution of the damping among the slave UAVs.
(a)
(b)Due to the consensus formation controller in (18), the roll and pitch angles in the rotational coordinates will gradually converge to the constant values. Therefore, the yaw coordinates of the ith UAV can be yielded as
(20)Jiψψ¨i=τisψ+τieψ,where Jiψ is the moment of inertia matrix in yaw coordinates and τisψ is the control torque applied on the ith UAV in yaw coordinates. τieψ is the torque caused by the interaction of the ith UAV with the environment in yaw coordinates. In the bilateral teleoperation system of UAVs with consensus formation control, the yaw angle of the ith UAV only depends on the repulsive forces which are produced by the obstacles in the obstructed environment, that is,
(21)ψ¨i=Jiψ−1Rifiez,for each slave UAV ensures that ψi−ψj,ψi∈L∞ for all i ∈ V, j ∈ Ni.Theorem 1.
For the single-master-multi-UAV teleoperation system with the consensus formation controller in (18), if there exist positive-definite matrices P and Q and any positive constants δi,εi such that the following inequalities hold:
(22)−αmI+T¯m1Q+T¯1mkm24P−1<0,−m+αsknksI+T¯1mP+T¯m1km24Q−1<0,βi>δi+εikc2lii+12δi∑j∈NiwijTij2−kv,∀i∈V,j∈Ni,then the overall teleoperation system is stable and the variables knqm−v1t−T1m, qi−qjt−Tij−dij, vi−vjt−Tij, ψi−ψj, q̇m, vi, and v̇i are bounded.Proof 1.
The total control forcefis applied on each UAV mainly contains two parts. In the first part, the velocity of the leader is required to converge to the velocity of the master. In this way, the operator can manipulate the movement of the slave formation by the master haptic device. In the second part, the UAVs are required to keep a relative position consensus with the neighbors. The control force fis applied on each UAV is linearly mixed. The first part is to track the master and the latter part is to keep consensus with the neighbors.
Thus, the control forcefis is divided into two parts, that is,
(23)fis=fis#+fis∗,where
(24)fis#=−biksvi−knqmt−Tm1−αsv̇i,fis∗=−kv∑j∈Nivi−vjt−Tij−kc∑j∈Niqi−qjt−Tij−dij−βiq̇i.
Firstly, the relationship between the master and the leader UAV is investigated. The total forces∑i=1nfis# which are applied on the leader UAV are
(25)fv1−knqmt−Tm11s#=−ks−αsv̇1.
Because the haptic device only has three translational output DOFs of force feedback, the rotational part of the bilateral teleoperation is unilateral. Therefore, only the stability of the translational part of the system is analyzed. The following Lyapunov function is defined as(26)H#=12q̇TmMmq̇m+m1γkm2knksρ1Tρ1+km2knknqm−v1T⋅knqm−v1+∫−T¯m10∫t+θtq̇mξTQq̇mξdξdθ+∫−T¯1m0∫t+θtv̇1ξTPv1ξdξdθ,where Q and P are positive-definite matrices. Obviously, function H# is always greater than zero.
With (9), Property 1, and Lemma 4, the time derivative of H# is
(27)Ḣ#=q̇mTfm+q̇mT⋅fh+kmknksρ1Tf1s#+f1e−m1kmknksρ1Tρ1+kmknknqm−v1Tknq̇m−v̇1+T¯m1q̇mTQq̇m−∫t−T¯m1tq̇mξTQq̇mξdξ+T¯1mv̇1TPv̇1−∫t−T¯1mtv̇1ξTPv̇1ξdξ.
With (18), (23), and Assumption 2, Ḣ# can be rewritten as
(28)Ḣ#≤q̇mTkmv1t−T1mt−knqm−αmq̇m+q̇mT⋅fh+kmknksρ1Tksknqmt−Tm1t−v1−αsv̇1+f1e−m1kmknksρ1Tρ1+kmknknqm−v1Tknq̇m−v̇1+T¯m1q̇mTQq̇m−∫t−T¯m1tq̇mξTQq̇mξdξ+T¯1mv̇1TPv̇1−∫t−T¯1mtv̇1ξTPv̇1ξdξ.
Noticing thatv1t−T1m−v1=−∫t−T1mtv̇1ξdξ,qmt−Tm1−qm=−∫t−T¯m1tq̇mξdξ, one has
(29)Ḣ#≤−kmq̇mT∫t−T1mtv̇1ξdξ−kmq̇mTknqm−v1−αmq̇mTq̇m+q̇mT⋅fh−kmρ1T∫t−T¯m1tq̇mξdξ+kmknρ1Tknqm−v1−m1+αskmknksρ1Tρ1+kmknksρ1Tfe+kmknqm−v1Tq̇m−kmknknqm−v1Tρ1+T¯m1q̇mTQq̇m−∫t−T¯m1tq̇mξTQq̇mξdξ+T¯1mv̇1TPv̇1−∫t−T¯1mtv̇1ξTPv̇1ξdξ.
With Lemma4, the following inequalities are obtained:
(30)−kmq̇mT∫t−T1mtv̇1ξdξ−∫t−T1mtv̇1ξTPv̇1ξdξ≤T¯1mkm24q̇mTP−1q̇m,−kmρ1T∫t−Tm1tq̇mξdξ−∫t−Tm1tq̇mξTQq̇mξdξ≤T¯m1km24ρ1TQ−1ρ1.
With (26), (27), (28), (29), and (30), one has
(31)Ḣ#≤−αmq̇mTq̇m+q̇mTfh−m1+αsknksρ1Tρ1+kmknksρ1Tf1e+T¯m1km24ρ1TQ−1ρ1+T¯m1q̇mTQq̇m+T¯1mkm24q̇mTP−1q̇m+T¯1mv̇1TPv̇1≤q̇mT−αmI+T¯m1Q+T¯1mkm24P−1q̇m+ρ1T−m1+αsknksI+T¯1mP+T¯m1km24Q−1ρ1+q̇Tmfh+kmknksρ1Tf1e.
With Assumption1, the third term of (31) is passive with respect to the force-velocity coupling fh,q̇m.
Similarly, the force-velocity couplingf1e,v1 is passive,
(32)v1Tf1e≥κ̇s,κs=m12v1Tv1.
Integrating equationv̇1=ρ1+γρ̇1,γ∈R+, one has
(33)v1=∫0Tρ1dt+γρ1,γ∈R+.
Letσ=μv1=μ∫0Tρ1dt+μγρ1,μ,γ∈R+. When the parameters μ and λ are properly selected, variable σ is proportional to the variable ρ1 with σ≈μγρ1.
Further, the new energy storage is selected as(34)κ¯s=μκs.
The time derivative ofκ¯s is
(35)κ¯̇s=μκ̇s≤μv1Tf1e=σTf1e≈μγρ1Tf1e.
Therefore, the fourth term of (31) is passive with respect to coupling f1e,ρ1.
In conclusion, to ensure the stability of the master and the leader,Ḣ# should be negative which is equivalent to the following inequalities:
(36)−αmI+T¯m1Q+T¯1mkm24P−1<0−m+αsknksI+T¯1mP+T¯m1km24Q−1<0.Setting αm, αs in (36) for the master and the leader UAVs implies that knqm−v1t−T1m,q̇m,v̇1∈L∞.
For the follower UAVs, the following Lyapunov function for theith UAV (i ∈ V) is constructed as
(37)Hi∗=12ṗiTMi−1pi−∫0tvifiedt,where pi∈R3 and Mi∈R3×3 are the momentum and the inertia matrix, respectively, of the ith UAV. With Assumption 1, Hi∗ is greater than zero.
With (23), the time derivation of Hi∗ is
(38)Ḣi∗=viTfis∗+fie−viTfie=−βivi2−kv∑j∈Nivi−vjt−Tij−kcviT∑j∈Niqi−qjt−Tij−dij=−βivi2−kvviT∑j∈Nivi−vj−kvviT∑j∈Nivj−vjt−Tij−kcviT∑j∈Niqi−qj−dij−kcviT∑j∈Niqj−qjt−Tij.
Considering the following total scaled energy function,(39)H∗=12qTL⊗In−diagARq+∑i=1nHi∗kc,where q=q1,…,qn, v=v1,…,vn, R is the relative distance matrix, and ⊗ is the standard Kronecker product.
(40)Ḣ∗=q̇TL⊗In−diagARq+∑i=1n−βi+kvkcvi2−kvkcvj−vjt−Tij−viT∑j∈Niqi−qj−dij−viT∑j∈Niqj−qjt−Tij=−∑i=1nβi+kvkcvi2−kvkc∑i=1n∑j∈NiviTvj−vjt−Tij−∑i=1n∑j∈NiviTqj−qjt−Tij.
Noticing thatqj−qjt−Tij=∫−Tij0q̇jξdξ, one has
(41)Ḣ∗≤−∑i=1nβi+kvkcvi2−∑i=1nviT∑j∈Ni∫−Tij0q̇jξdξ−kvkc∑i=1nviT∑j∈Ni∫−Tij0v̇jξdξ.
By integrating the above equation from zero tot with Lemma 5, we obtain
(42)H∗−H0≤−∑i=1nβi+kvkcvi22+∑i=1n∑j∈Niδi+εi2vi22+Tij22δivj22+Tij22εiv̇j22,where δi, εi > 0. Using the fact that H∗ > 0,
(43)H0≥∑i=1n∑j∈Niβi+kvkclii−δi+εi2vi22−Tij22δivj22.
DefiningΓ=v122,…,vn22∈Rn and Ω∈Rn×n,
(44)Ωn×n=β1kcl11−δ12−T1222δ1⋯−T1n22δ1−T2122δ2β2kcl22−δ22⋯−T2n22δ2⋮⋮⋱⋮−Tn122δn−Tn222δn⋯βnkclnn−δn2.
ThenH(0) can be rewritten as H0≥1nTΩΓ. There exists ℓl≥0 for l∈V such that H0≥ℓ1v122+⋯+ℓnvn22 and vi,v̇i∈L2. Therefore, setting the damping injection βi,
(45)βi>δi+εikc2lii+12δi∑j∈NiwijTij2−kv,for any arbitrary δi,εi>0. Each slave UAV ensures that vi∈L2 and H∗∈L∞, which in turn implies that qi−qjt−Tij−dij,vi−vjt−Tij,vi∈L∞, for all i∈V,j∈Ni.
Combining (21), (36), and (45), one has knqm→v1t−T1m, qi−qjt−Tij→dij, vi→vjt−Tij, ψi→ψj. The proof is completed.


## 3.3. Flux-Based Collision Avoidance

In general, the artificial potential fields [15] are usually used to avoid the obstacles. The gradient of the potential around the obstacles is used as the repulsive force exerted on the vehicle and the potential disappears outside of the influence region.Here, a flux-conserved force field is proposed to keep the consensus formation away from the obstacles and escape the local minima. Inspired by electromagnetics, each obstacle can be regarded as the electronic objects. LetEiobs denote the force field and C denote any closed curves around the obstacle in the 3D space. The flux is defined as
(46)Wiobs=∮Cn⋅Eiobsds,where Wiobs is the flux of obstacle i for i ∈ {1, …, nobs} and n is the normal vector of the curve.Similar to the conservation law of electric flux in the electromagnetic field, the most important feature of the flux-conserved force field is that the flux of the force field is conserved for any curves. In Figure5, for any curve C surrounding obstacle i, the value of Wiobs remains the same. The repulsive force of the ith UAV from the obstacles is
(47)fie=kobs∑i=1nobs∂Wiobs∂qi−xoi,where qi ∈ R3 is the position of the ith UAV and xoi∈R3 is the position of obstacle i. kobs is a positive coefficient. For any curve surrounding an electric point, the value of Wiobs is always same.Figure 5
The flux of the force field is conserved for any curves.With the flux-conserved force field, it can be proved that there are no local minima in this force field. Proving can be demonstrated as follows. As shown in Figure6, assuming a local minimum, point P appears around obstacle i. The local minimum is a point where the sum of the force field is close to zero. In this figure, the value of the flux generated by the obstacle is defined as positive and the value of the flux generated by the target is defined as negative. Firstly, considering an imaginary little closed curve around the local minimum, point P and the force on the curve point toward the inside. Then because the force filed is flux conserved in the closed curve, the curve must have a target with a negative flux value. Therefore, if there is not a target in the closed curve, the local minimum around the obstacle must certainly not occur.Figure 6
No local minima in flux-conserved force field.


## 4. Top-Down 3D Optimal Persistent Graph

In the above controller, each slave UAV is required to communicate with its neighbors to keep a consensus formation. However, some interactions between the UAVs are not necessary and they will make the communication much more complex. If consensus can be reached with less neighbor information, the communication links and energy dissipation of the formation can be decreased. The topology could be optimized in real-time to achieve the energy optimal formation. In the following, a top-down strategy of the 3D optimal persistent graph is proposed to simplify the topology.


### 4.1. Generation of the Min-Weighted Rigid Graph

It is well known that the optimal persistent graph is the persistence of the min-weighted rigid graph. Therefore, the min-weighted rigid graph should be generated firstly.The generation of the min-weighted rigid graph in 3D space mainly contains two steps: (1) generating the infinitesimally rigid graphs and (2) generating the min-weighted rigid graph that has a minimally weighted sum in all infinitesimally rigid graphs. The program of generating the min-weighted rigid graph is shown in Algorithm1. For example, the candidate minimally rigid graphs for the five UAVs are shown in Figure 7. Calculate the weights of the candidate graphs with the length of all edges. The min-weighted rigid graph is a candidate graph with the smallest weight.Algorithm 1: Generation of the min-weighted rigid graph in 3D space.
Initialize the complete graphG = (V, E, A) with n nodes, where n represents the number of UAVs;
Calculate the length of all edges of the complete graph;
Sequence all the edges by increasing;
Build the rigidity matrixM of the complete graph according to the sequence of edges;
Initialize matrixMc as the first row of M;
fori=1:nn−1/2
if rank M≤3n−6
Add the next row ofM to Mc to form a new matrix M’;
if M’ is full rank
Record the edge set inM corresponding to the row;
end
end
end
Draw the min-weighted rigid graph;Remark 2.
If the min-weighted rigid graph is drawn, the adjacency matrix corresponding to the min-weighted rigid graph can be obtained. Therefore, the optimized adjacency matrix can be directly adopted to describe the neighborsNi in the consensus formation controller.Figure 7
The candidate minimally rigid graphs for five UAVs in 3D space.


### 4.2. Generation of 3D Optimal Persistent Graph

In the following, two basic operations of rigid graph extension are introduced [31]. As shown in Figure 8, one of the operations is adding vertex. In this figure, j and k are the original two vertexes of the optimal persistent graph. The operation of adding a vertex is to add vertex i, directed edges ij→ and ik→. If the original minimally rigid graph has n vertexes, it is easy to know that the newly obtained graph containing n + 1 vertexes is also a minimally rigid graph, that is, the operation of adding the vertex has a rigid characteristic. Another operation is splitting edge. The operation is adding vertex i, connecting the associated edges, and removing the original edge. Similarly, the operation of splitting the edge also has a rigid characteristic.Figure 8
Two basic operations of rigid graph extension. (a)Vertex addition and (b) edge splitting.
(a)
(b)In 2D space, the optimal persistent graph can be generated based on the operation of vertex adding. An example of generating the optimal persistent graph of five UAVs in 2D space is shown in Figure9.Figure 9
Formation based on the operation of adding a vertex.However, in 3D space, the generation of the 3D optimal persistent graph is different from that in 2D space. To obtain the optimal persistent graph in 3D space, a top-down strategy that utilizes the different altitudes of UAVs in the formation will be used. The 3D space will be divided into several subspaces according to the altitudes of UAVs. The generation of the 3D optimal persistent graph mainly contains two steps: (1) generation of the min-weighted rigid graph, which can be seen in Figure7 and (2) generation of the optimal persistent graph in subspaces based on the min-weighted rigid graph. In Figure 10(a), a contour map is drawn according the altitude of UAVs. We can see from the figure that the UAVs are distributed on different contour lines and the spaces can be divided into three subspaces. In Figure 10(b), the first subspace consists of UAVs 1, 2, and 3, connecting UAV1 and UAV3 initially and performing the operation of adding vertex UAV2. In Figure 10(c), the second subspace consists of UAVs 2, 3, and 4. The optimal persistent graph on the second subspace is generated by performing the operation of adding vertex UAV4. Similarly, the third subspace consists of UAVs 4 and 5. The optimal persistent graph on the third subspace is generated by performing the operation of adding vertex UAV5. Combined with the optimal persistent graphs in different subspaces, the 3D optimal persistent graph is generated in Figure 10(d).Figure 10
Generation of the 3D minimally persistent graph based on G10 with the top-down strategy. (a) The contour map of UAVs. (b) The formation generation in the first subspace. (c) The formation generation in the second subspace. (d) The formation generation in the third subspace.
(a)
(b)
(c)
(d)


### 4.3. Optimal Persistent Graph with Switching Formation

The goal of persistence is to maintain the formation of multiple UAVs in movement, but the formation needs to be changed according to the environment. In this situation, the team can choose to maintain a rigid formation just to adjust the distance between the relevant UAVs, making the robot able to pass through the obstacles. However, when the distance cannot meet the needs of the case, a new 3D optimal persistent graph needs to be generated. The program of generating the 3D optimal persistent graph in the obstructed environment is shown in Algorithm2.Algorithm 2: Generation of the 3D optimal persistent graph in obstructed environment.
Generation of the min-weighted rigid graph in 3D space;
Detection of the environment;
Screening the expected rigid graphs to meet environmental constraints;
Calculate the adjacency matrixA = [aij] and the connectivity of each vertex d−(i);
Establish the subspaces according to the contour map and the number of UAVs in subspaces isni;
Establish the initial directed edge of the leader and any of its followers in the first subspace;
Set the out-degrees of the vertexes of the initial directed edge to 2;
for space = 1 : subspaces
for i = 1 : ni
if d−(j) = 2, d−(k) = 2, and aij = aik = 1,
Establish the directed edgeVij, Vik and d−(k) = 2;
if for any vertex i, the out-degrees d−(i) = 2
The 3D optimal persistent graph is generated;
end
end
end
end


### 4.4. Energy Dissipation of Switching Formation

To verify whether the energy dissipation of the formation is decreased or not, the communication cost is calculated. A common low-energy radio model which is usually used in the wireless sensor network is introduced (Yan et al., 2014). It is assumed that the UAVs are all equipped with the transmitter and the receiver. In this model, the transmitter runs the radio equipment and the energy amplifier while the receiver only runs the radio equipment. The energy dissipation depends on the distance between the transmitter and the receiver. Therefore, to transmit an l-bit message at distance d, the energy dissipation is
(48)ETxl,d=ETx_elecl+ETx_ampl,d=Eelec+lεfris_ampd2.And to receive this message, the energy dissipation is
(49)ERxl,d=ETx_elecl=lEelec,where the electronics energy Eelec depends on many factors, such as digital coding, modulation, and spreading of the signal. The amplifier energy εfris_amp depends on the acceptable intensity and bit-error rate.


## 4.1. Generation of the Min-Weighted Rigid Graph

It is well known that the optimal persistent graph is the persistence of the min-weighted rigid graph. Therefore, the min-weighted rigid graph should be generated firstly.The generation of the min-weighted rigid graph in 3D space mainly contains two steps: (1) generating the infinitesimally rigid graphs and (2) generating the min-weighted rigid graph that has a minimally weighted sum in all infinitesimally rigid graphs. The program of generating the min-weighted rigid graph is shown in Algorithm1. For example, the candidate minimally rigid graphs for the five UAVs are shown in Figure 7. Calculate the weights of the candidate graphs with the length of all edges. The min-weighted rigid graph is a candidate graph with the smallest weight.Algorithm 1: Generation of the min-weighted rigid graph in 3D space.
Initialize the complete graphG = (V, E, A) with n nodes, where n represents the number of UAVs;
Calculate the length of all edges of the complete graph;
Sequence all the edges by increasing;
Build the rigidity matrixM of the complete graph according to the sequence of edges;
Initialize matrixMc as the first row of M;
fori=1:nn−1/2
if rank M≤3n−6
Add the next row ofM to Mc to form a new matrix M’;
if M’ is full rank
Record the edge set inM corresponding to the row;
end
end
end
Draw the min-weighted rigid graph;Remark 2.
If the min-weighted rigid graph is drawn, the adjacency matrix corresponding to the min-weighted rigid graph can be obtained. Therefore, the optimized adjacency matrix can be directly adopted to describe the neighborsNi in the consensus formation controller.Figure 7
The candidate minimally rigid graphs for five UAVs in 3D space.


## 4.2. Generation of 3D Optimal Persistent Graph

In the following, two basic operations of rigid graph extension are introduced [31]. As shown in Figure 8, one of the operations is adding vertex. In this figure, j and k are the original two vertexes of the optimal persistent graph. The operation of adding a vertex is to add vertex i, directed edges ij→ and ik→. If the original minimally rigid graph has n vertexes, it is easy to know that the newly obtained graph containing n + 1 vertexes is also a minimally rigid graph, that is, the operation of adding the vertex has a rigid characteristic. Another operation is splitting edge. The operation is adding vertex i, connecting the associated edges, and removing the original edge. Similarly, the operation of splitting the edge also has a rigid characteristic.Figure 8
Two basic operations of rigid graph extension. (a)Vertex addition and (b) edge splitting.
(a)
(b)In 2D space, the optimal persistent graph can be generated based on the operation of vertex adding. An example of generating the optimal persistent graph of five UAVs in 2D space is shown in Figure9.Figure 9
Formation based on the operation of adding a vertex.However, in 3D space, the generation of the 3D optimal persistent graph is different from that in 2D space. To obtain the optimal persistent graph in 3D space, a top-down strategy that utilizes the different altitudes of UAVs in the formation will be used. The 3D space will be divided into several subspaces according to the altitudes of UAVs. The generation of the 3D optimal persistent graph mainly contains two steps: (1) generation of the min-weighted rigid graph, which can be seen in Figure7 and (2) generation of the optimal persistent graph in subspaces based on the min-weighted rigid graph. In Figure 10(a), a contour map is drawn according the altitude of UAVs. We can see from the figure that the UAVs are distributed on different contour lines and the spaces can be divided into three subspaces. In Figure 10(b), the first subspace consists of UAVs 1, 2, and 3, connecting UAV1 and UAV3 initially and performing the operation of adding vertex UAV2. In Figure 10(c), the second subspace consists of UAVs 2, 3, and 4. The optimal persistent graph on the second subspace is generated by performing the operation of adding vertex UAV4. Similarly, the third subspace consists of UAVs 4 and 5. The optimal persistent graph on the third subspace is generated by performing the operation of adding vertex UAV5. Combined with the optimal persistent graphs in different subspaces, the 3D optimal persistent graph is generated in Figure 10(d).Figure 10
Generation of the 3D minimally persistent graph based on G10 with the top-down strategy. (a) The contour map of UAVs. (b) The formation generation in the first subspace. (c) The formation generation in the second subspace. (d) The formation generation in the third subspace.
(a)
(b)
(c)
(d)


## 4.3. Optimal Persistent Graph with Switching Formation

The goal of persistence is to maintain the formation of multiple UAVs in movement, but the formation needs to be changed according to the environment. In this situation, the team can choose to maintain a rigid formation just to adjust the distance between the relevant UAVs, making the robot able to pass through the obstacles. However, when the distance cannot meet the needs of the case, a new 3D optimal persistent graph needs to be generated. The program of generating the 3D optimal persistent graph in the obstructed environment is shown in Algorithm2.Algorithm 2: Generation of the 3D optimal persistent graph in obstructed environment.
Generation of the min-weighted rigid graph in 3D space;
Detection of the environment;
Screening the expected rigid graphs to meet environmental constraints;
Calculate the adjacency matrixA = [aij] and the connectivity of each vertex d−(i);
Establish the subspaces according to the contour map and the number of UAVs in subspaces isni;
Establish the initial directed edge of the leader and any of its followers in the first subspace;
Set the out-degrees of the vertexes of the initial directed edge to 2;
for space = 1 : subspaces
for i = 1 : ni
if d−(j) = 2, d−(k) = 2, and aij = aik = 1,
Establish the directed edgeVij, Vik and d−(k) = 2;
if for any vertex i, the out-degrees d−(i) = 2
The 3D optimal persistent graph is generated;
end
end
end
end


## 4.4. Energy Dissipation of Switching Formation

To verify whether the energy dissipation of the formation is decreased or not, the communication cost is calculated. A common low-energy radio model which is usually used in the wireless sensor network is introduced (Yan et al., 2014). It is assumed that the UAVs are all equipped with the transmitter and the receiver. In this model, the transmitter runs the radio equipment and the energy amplifier while the receiver only runs the radio equipment. The energy dissipation depends on the distance between the transmitter and the receiver. Therefore, to transmit an l-bit message at distance d, the energy dissipation is
(48)ETxl,d=ETx_elecl+ETx_ampl,d=Eelec+lεfris_ampd2.And to receive this message, the energy dissipation is
(49)ERxl,d=ETx_elecl=lEelec,where the electronics energy Eelec depends on many factors, such as digital coding, modulation, and spreading of the signal. The amplifier energy εfris_amp depends on the acceptable intensity and bit-error rate.


## 5. Human-in-the-Loop Simulations

### 5.1. Experimental Setup

The experimental setup for the bilateral teleoperation system is shown in Figure11. The goal of the experiment is to enable the human operator to remotely control the consensus formation of UAVs in the obstructed environment and to provide the operator with some useful haptic and visual feedback cues that convey information about the UAV state and the surroundings. The experiments are carried out in ten runs and the experimental values are all averaged. The framework of the experimental setup is shown in Figure 12. The testbed consists of a haptic device as the master device, five UAVs as the slave robots, and Gazebo as the multirobot simulation platform and ROS operating system which can develop the communication between the master haptic device and the slave UAVs. In each UAV, the microcontroller runs a single process implementing the attitude controller and the inertial measurement unit (IMU) measures the attitude of the UAV. The onboard computer runs a process implementing the position controller, formation controller among the UAVs, and communication with the human operator. The UAVs use the ROS topics to communicate with each other via TCP/IP socket.Figure 11
The experimental setup for the bilateral teleoperation system.Figure 12
The framework of the experimental setup.The commercial Sensable PHANTOM Omni is selected as the master device. The haptic device has six degrees-of-freedom (DOF) translational and rotational inputs and three DOF translational force outputs. The refresh rate of the feedback force is 1000 Hz. In the experiments, the rotational part of the bilateral teleoperation is unilateral. The rotational movements of the device are ignored, leaving the translational coordinates as the actuatable DOFs. The device can be automatically calibrated by placing the stylus in the inkwell and starting the haptically enabled ROS node. The indicator light in the inkwell of the device will be lit a steady blue when the workspace and force are calibrated properly. The device provides a workspace of 160 × 120 × 70 mm in width, height, and depth, respectively, and a maximum exertable force of 3.3 N.Five Parrot AR Drones are selected as the slave robots. The mechanical structure comprises four rotors attached to the four ends of a crossing to which the battery and the Wi-Fi hardware are attached. An ultrasonic module is fixed for measuring the drone’s flight altitude with a maximum range of six meters. A downward camera is used as the optical flow sensor to compute the motion field. Moreover, a forward camera provides the operator several visual feedbacks about the obstructed environment. The net weight is 436 g and the size is 51.5 × 51.5 centimeters.In the Gazebo platform, communications between the human operator and slave UAVs, as well as communications between UAVs, completely simulate the real-time communications. The operator can send commands on UDP port 5556 with 30 times per second. Meanwhile, the information about the UAV, like its status, position, speed, and engine rotation speed, is sent to the master on UDP port 5554 with 200 times per second. A video stream of the frontward or downward camera is sent by the UAV to the master/UAV through TCP socket. The packet format of the control command and navigation and video data is shown in [36]. The time-varying delays are shown in Figure 13. Moreover, the fluctuations on the delays are significant and cannot be ignored.Figure 13
The time-varying delays for the semiexperiments.


### 5.2. Experimental Results and Discussions

The effectiveness of the consensus formation controller for the bilateral teleoperation system in the obstructed environment is evaluated by the human-in-the-loop simulations. The initial conditions for master/slave devices areqm = (0, 0, 0) m, q1(0) = (3, 3, 2.5) m, q2(0) = (0, 0, 0) m, q3(0) = (6, 4, 0) m, q4(0) = (0, 5, 0) m, q5(0) = (5, 0, −1) m, vi(0) = (0, 0, 0) m/s, Θi(0) = (0, 0, 0) rad, mi = 0.436 kg for i = 1, …, 5. There are two obstacles including the corner and circular obstacles in the slave environment. The time-varying delays are Tm1 = 0.3rand() s, T1m = Tij = Tji = 0.8rand() s where rand() is a random number with a variable of 0.01, T¯m1 = 0.7 s, T¯1m = T¯ij = T¯ji = 1.2 s. The desired relative distance matrix of the consensus formation is R=06666;602222;622022;622022;622220. The turning procedures of the control parameters are as follows:
(1)
Settingkn to solve the kinematic/dynamic dissimilarity between the master and slave. The workspace of the haptic device is 160 × 120 × 70 mm. By setting kn = 10 (1/s), the maximum velocities of the UAV are up to 1.6 m/s, 1.2 m/s, and 0.7 m/s, respectively. Thus, the UAV can fly with unbounded workspace and the kinematic/dynamic dissimilarity is solved by the position-velocity mapping.
(2)
Settingks to achieve master-slave velocity tracking. In the unilateral teleoperation case with no time delays, control parameters km, Tm, and Ts are all equal to zero. When ks is set to 1 (Ns/m), the speed of the UAV can completely track that of the master.
(3)
Settingkm to achieve master-slave force tracking. In the bilateral teleoperation case with no time delays, the time delays Tm, Ts are equal to zero. When the control parameter km is set to 1 (Ns/m), the master feedback force fm can completely track the slave control force fs.
(4)
Settingkv and kc to achieve velocity synchronization and formation control of the slave UAVs. In the bilateral teleoperation case with no time delays, the velocities of agents can keep the same in a steady state if kv is set to 0.5 (Ns/m). Simultaneously, the consensus formation control can be achieved while enforcing a safe and minimum distance between any two UAVs at all time if kc is set to 0.006 (Ns/m).
(5)
Setting the damping coefficientsαm, αs, and βi to overcome the impact of time-varying delays. Let the positive-definite matrices P = Q = I, where I is an identity matrix. Thus, combined with the previously selected coefficients kn, ks, km, max{Tm1}, max{T1m}, max{Tij}, and max{Tji}, the damping coefficient αm and αs can be obtained by Theorem 1.
Also, let any positive constantsδi = εi = 1. The Laplacian matrices in the fully connected, G2, G6, and G9 modes are computed as follows:
(50)L=4−1−1−1−1−14−1−1−1−1−14−1−1−1−1−14−1−1−1−1−14,LG2=3−100−1030−1−10−1200−10−14000−1−14,LG6=4000−1−1300−1−1−1400−10−13000−1−14,LG9=4−1−100040−1−10−1300−10−140−100−13.
Combined with the previously selected coefficientskc and kv, the damping coefficient βi in different modes can be obtained by Theorem 1 (45).
(6)
Setting the factor of obstacle repulsive forcekobs to guarantee that the feedback force is not greater than 3.3 N. Therefore, the parameters of the consensus formation in the human-in-the-loop simulations are given as follows: kn = 10, kc = 0.006, kv = 0.5, kobs = 0.0009, P = Q = I, δi = εi = 1, km = ks = 1, αm = 2, αs = 15, βi = 5 for i∈V in a fully connected mode, βi = 2 in G2, G6, and G9 modes, respectively.


#### 5.2.1. Bilateral Aerial Teleoperation with Consensus Formation

The ability of the consensus formation to adapt to the obstructed environment is evaluated by the human-in-the-loop simulations. During the movement, the UAVs should always maintain a rigid formation and avoid collisions with the neighbor UAVs and obstacles. In the consensus formation with the unoptimized topology, the communications are fully connected, that is, the numbers of communication edges are always ten regardless of changes in the neighborhood area.In Figure14, the UAVs in the bilateral aerial teleoperation reach the local minima in the corner obstacle and the consensus formation failed. In Figure 15, the consensus formation with flux-based collision avoidance function is shown. The topologies with full connections are marked by the black solid lines. The obstacles are considered as the electronic objects. During the time interval t ∈ [0 s, 150 s], the operator sends the desired velocity to the leader UAV. Under the action of the consensus formation controller, the UAVs are required to form a rigid formation, keep velocity consensus with the master, and avoid collisions with obstacles. When any of the slave UAVs encounters obstacles, the rigid formation breaks up to keep a safe distance from the obstacles. In the flux-conserved force field, the formation can escape the local minima from the corner obstacle. Until all UAVs leave away from the obstacles, the rigid formation is again formed. The relative distances between the slave UAVs are shown in Figure 16. It can be seen that the UAVs can avoid collisions with the neighbors and obstacles while keeping a desired rigid formation in the absence of obstacles. The yaw angles of the slave UAVs are shown in Figure 17.Figure 14
The consensus formation reaches the local minima.Figure 15
The trajectories of the consensus formation.Figure 16
The relative distances.Figure 17
The yaw angles.The velocity information is shown in Figures18(a)–18(c). rsx, rsy, and rsz are the components of the commanded velocity rs sent to the leader UAV. vix, viy, and viz are the components of the velocities of the slave UAVs for i = 1, …, 5. It can be seen that the velocities of the leader UAV v1 are converged to the commanded velocity rs generated by the master haptic device. The velocities of the followers v2, v3, v4, and v5 are converged to the leader’s velocity v1. When the obstacles are encountered, the velocities of the slave UAVs are decreased to prevent collisions with the obstacles. Once all UAVs leave away from the obstacles, the velocities of all of them are increased to track the master velocity.Figure 18
Velocities and force information. (a) The velocities in thex-axis. (b) The velocities in the y-axis. (c) The velocities in the z-axis. (d) The master/slave force.
(a)
(b)
(c)
(d)The force information is shown in Figure18(d). fmx, fmy, and fmz are the components of the feedback force fm applied on the haptic device. fsx, fsy, and fsz are the components of the summed control force fs applied on the slave UAVs. When the rigid formation moves in a stable state, the forces are zero and the operator has no feeling about the environment. Once the obstacles are encountered, the forces increase immediately and the human operator at the master side can perceive the repulsive forces. Therefore, the force reflection makes the operator perceive the movement of the rigid formation better.Moreover, the maneuverability performance of formation and the perceptual sensitivity of the remote environment are shown in Figure19. The cross correlation is to assess the velocity tracking error and smaller values indicate greater velocity tracking. Just noticeable difference (JND) is to assess the force tracking accuracy and smaller values indicate greater perceptual sensitivity to the haptic cue. From the figure, we can see that compared with that of the traditional wave variable method, the performance of velocity tracking with local damping in the consensus formation controller is greater and the force tracking accuracy is asymptotically the same as that of the traditional wave variable method. The time-delayed bilateral aerial teleoperation system with local damping can be in a steady state to guarantee better velocity and force tracking in free motion or obstacle avoidance.Figure 19
Maneuverability performance and perceptual sensitivity.


#### 5.2.2. Bilateral Aerial Teleoperation with Energy-Optimized Consensus Formation

To reduce the communication complexity and energy dissipation of the formation, the 3D top-down optimal persistent graph is used to simplify the formation topology. Moreover, the communication costs are calculated in the unoptimized and optimized mode, respectively. The parameters of the low-energy radio model are set asEelec = 50 nJ/bit, εfris_amp = 10 J/bit/m2, and l = 400 bit.The minimally rigid graphs of the five UAVs are shown in Figure7. There are a total of ten minimally rigid graphs. G1–G10 are the different minimally rigid graphs. In the optimized mode, a min-weighted rigid graph is selected from these minimally rigid graphs according to the weights of the communication edges. The weights of the ten candidate minimally rigid graphs at initial are shown in Figure 20. It can be seen from the figure that the initial min-weighted rigid graph is G2. On the basis of the min-weighted rigid graph, the top-down strategy is used to select the optimal persistent graph. The optimal persistent graphs of the consensus formation used in the optimized mode are shown in Figure 21. It can be seen that the minimally rigid graphs G6 and G9 are selected as the min-weighted rigid graph when the formation reaches the corner obstacle and the circular obstacle, respectively. In the absence of obstacles, the minimally rigid graph G2 is selected. Obviously, the communication links are 3n − 6 = 9, and this value is less than the full connections. Moreover, the weights of the communication are always minimal at different times. The average communication costs of the slave UAVs in fully connected, min-weighted rigid and optimal persistent graphs are shown in Figure 22. Obviously, compared with those with the fully connected graph and min-weighted rigid graph, the communication costs and links with the optimal persistent graph are minimal. The trajectories of the optimized consensus formation are shown in Figure 23. The topologies with optimal persistent graph are marked by the black dashed lines. Moreover, the 3D graphs of the experimental tests with one UAV, three UAVs, and five UAVs which implement the energy-optimized consensus formation scheme for the time-delayed bilateral teleoperation system are shown in Figure 24. Obviously, the optimized consensus formation can be achieved and the advantages are shown in Figure 25. Using the optimal persistent graphs in formation, the trajectories around the obstacles are smoother and the running time is significantly shortened. Due to the shorter running time, the consensus formation with the optimal persistent graph can fly farther.Figure 20
The weights of the initial candidate minimally rigid graphs. G2 is the min-weighted rigid graph at initial.Figure 21
The optimal persistent graphs of the formation used in the optimized mode.Figure 22
Average communication costs and links in different graphs.Figure 23
The consensus formation with 3D optimal rigid graph.Figure 24
3D Graphs of the experimental tests. (a) One UAV in the left view. (b) Three UAVs in the left view. (c) Five UAVs in the top view.
(a)
(b)
(c)Figure 25
The characteristics of consensus formation controller in different graphs. SN denotes the smoothness of the trajectories around the corner or circle obstacle and RT denotes runtime time around the corner or circle obstacle.


## 5.1. Experimental Setup

The experimental setup for the bilateral teleoperation system is shown in Figure11. The goal of the experiment is to enable the human operator to remotely control the consensus formation of UAVs in the obstructed environment and to provide the operator with some useful haptic and visual feedback cues that convey information about the UAV state and the surroundings. The experiments are carried out in ten runs and the experimental values are all averaged. The framework of the experimental setup is shown in Figure 12. The testbed consists of a haptic device as the master device, five UAVs as the slave robots, and Gazebo as the multirobot simulation platform and ROS operating system which can develop the communication between the master haptic device and the slave UAVs. In each UAV, the microcontroller runs a single process implementing the attitude controller and the inertial measurement unit (IMU) measures the attitude of the UAV. The onboard computer runs a process implementing the position controller, formation controller among the UAVs, and communication with the human operator. The UAVs use the ROS topics to communicate with each other via TCP/IP socket.Figure 11
The experimental setup for the bilateral teleoperation system.Figure 12
The framework of the experimental setup.The commercial Sensable PHANTOM Omni is selected as the master device. The haptic device has six degrees-of-freedom (DOF) translational and rotational inputs and three DOF translational force outputs. The refresh rate of the feedback force is 1000 Hz. In the experiments, the rotational part of the bilateral teleoperation is unilateral. The rotational movements of the device are ignored, leaving the translational coordinates as the actuatable DOFs. The device can be automatically calibrated by placing the stylus in the inkwell and starting the haptically enabled ROS node. The indicator light in the inkwell of the device will be lit a steady blue when the workspace and force are calibrated properly. The device provides a workspace of 160 × 120 × 70 mm in width, height, and depth, respectively, and a maximum exertable force of 3.3 N.Five Parrot AR Drones are selected as the slave robots. The mechanical structure comprises four rotors attached to the four ends of a crossing to which the battery and the Wi-Fi hardware are attached. An ultrasonic module is fixed for measuring the drone’s flight altitude with a maximum range of six meters. A downward camera is used as the optical flow sensor to compute the motion field. Moreover, a forward camera provides the operator several visual feedbacks about the obstructed environment. The net weight is 436 g and the size is 51.5 × 51.5 centimeters.In the Gazebo platform, communications between the human operator and slave UAVs, as well as communications between UAVs, completely simulate the real-time communications. The operator can send commands on UDP port 5556 with 30 times per second. Meanwhile, the information about the UAV, like its status, position, speed, and engine rotation speed, is sent to the master on UDP port 5554 with 200 times per second. A video stream of the frontward or downward camera is sent by the UAV to the master/UAV through TCP socket. The packet format of the control command and navigation and video data is shown in [36]. The time-varying delays are shown in Figure 13. Moreover, the fluctuations on the delays are significant and cannot be ignored.Figure 13
The time-varying delays for the semiexperiments.


## 5.2. Experimental Results and Discussions

The effectiveness of the consensus formation controller for the bilateral teleoperation system in the obstructed environment is evaluated by the human-in-the-loop simulations. The initial conditions for master/slave devices areqm = (0, 0, 0) m, q1(0) = (3, 3, 2.5) m, q2(0) = (0, 0, 0) m, q3(0) = (6, 4, 0) m, q4(0) = (0, 5, 0) m, q5(0) = (5, 0, −1) m, vi(0) = (0, 0, 0) m/s, Θi(0) = (0, 0, 0) rad, mi = 0.436 kg for i = 1, …, 5. There are two obstacles including the corner and circular obstacles in the slave environment. The time-varying delays are Tm1 = 0.3rand() s, T1m = Tij = Tji = 0.8rand() s where rand() is a random number with a variable of 0.01, T¯m1 = 0.7 s, T¯1m = T¯ij = T¯ji = 1.2 s. The desired relative distance matrix of the consensus formation is R=06666;602222;622022;622022;622220. The turning procedures of the control parameters are as follows:
(1)
Settingkn to solve the kinematic/dynamic dissimilarity between the master and slave. The workspace of the haptic device is 160 × 120 × 70 mm. By setting kn = 10 (1/s), the maximum velocities of the UAV are up to 1.6 m/s, 1.2 m/s, and 0.7 m/s, respectively. Thus, the UAV can fly with unbounded workspace and the kinematic/dynamic dissimilarity is solved by the position-velocity mapping.
(2)
Settingks to achieve master-slave velocity tracking. In the unilateral teleoperation case with no time delays, control parameters km, Tm, and Ts are all equal to zero. When ks is set to 1 (Ns/m), the speed of the UAV can completely track that of the master.
(3)
Settingkm to achieve master-slave force tracking. In the bilateral teleoperation case with no time delays, the time delays Tm, Ts are equal to zero. When the control parameter km is set to 1 (Ns/m), the master feedback force fm can completely track the slave control force fs.
(4)
Settingkv and kc to achieve velocity synchronization and formation control of the slave UAVs. In the bilateral teleoperation case with no time delays, the velocities of agents can keep the same in a steady state if kv is set to 0.5 (Ns/m). Simultaneously, the consensus formation control can be achieved while enforcing a safe and minimum distance between any two UAVs at all time if kc is set to 0.006 (Ns/m).
(5)
Setting the damping coefficientsαm, αs, and βi to overcome the impact of time-varying delays. Let the positive-definite matrices P = Q = I, where I is an identity matrix. Thus, combined with the previously selected coefficients kn, ks, km, max{Tm1}, max{T1m}, max{Tij}, and max{Tji}, the damping coefficient αm and αs can be obtained by Theorem 1.
Also, let any positive constantsδi = εi = 1. The Laplacian matrices in the fully connected, G2, G6, and G9 modes are computed as follows:
(50)L=4−1−1−1−1−14−1−1−1−1−14−1−1−1−1−14−1−1−1−1−14,LG2=3−100−1030−1−10−1200−10−14000−1−14,LG6=4000−1−1300−1−1−1400−10−13000−1−14,LG9=4−1−100040−1−10−1300−10−140−100−13.
Combined with the previously selected coefficientskc and kv, the damping coefficient βi in different modes can be obtained by Theorem 1 (45).
(6)
Setting the factor of obstacle repulsive forcekobs to guarantee that the feedback force is not greater than 3.3 N. Therefore, the parameters of the consensus formation in the human-in-the-loop simulations are given as follows: kn = 10, kc = 0.006, kv = 0.5, kobs = 0.0009, P = Q = I, δi = εi = 1, km = ks = 1, αm = 2, αs = 15, βi = 5 for i∈V in a fully connected mode, βi = 2 in G2, G6, and G9 modes, respectively.


### 5.2.1. Bilateral Aerial Teleoperation with Consensus Formation

The ability of the consensus formation to adapt to the obstructed environment is evaluated by the human-in-the-loop simulations. During the movement, the UAVs should always maintain a rigid formation and avoid collisions with the neighbor UAVs and obstacles. In the consensus formation with the unoptimized topology, the communications are fully connected, that is, the numbers of communication edges are always ten regardless of changes in the neighborhood area.In Figure14, the UAVs in the bilateral aerial teleoperation reach the local minima in the corner obstacle and the consensus formation failed. In Figure 15, the consensus formation with flux-based collision avoidance function is shown. The topologies with full connections are marked by the black solid lines. The obstacles are considered as the electronic objects. During the time interval t ∈ [0 s, 150 s], the operator sends the desired velocity to the leader UAV. Under the action of the consensus formation controller, the UAVs are required to form a rigid formation, keep velocity consensus with the master, and avoid collisions with obstacles. When any of the slave UAVs encounters obstacles, the rigid formation breaks up to keep a safe distance from the obstacles. In the flux-conserved force field, the formation can escape the local minima from the corner obstacle. Until all UAVs leave away from the obstacles, the rigid formation is again formed. The relative distances between the slave UAVs are shown in Figure 16. It can be seen that the UAVs can avoid collisions with the neighbors and obstacles while keeping a desired rigid formation in the absence of obstacles. The yaw angles of the slave UAVs are shown in Figure 17.Figure 14
The consensus formation reaches the local minima.Figure 15
The trajectories of the consensus formation.Figure 16
The relative distances.Figure 17
The yaw angles.The velocity information is shown in Figures18(a)–18(c). rsx, rsy, and rsz are the components of the commanded velocity rs sent to the leader UAV. vix, viy, and viz are the components of the velocities of the slave UAVs for i = 1, …, 5. It can be seen that the velocities of the leader UAV v1 are converged to the commanded velocity rs generated by the master haptic device. The velocities of the followers v2, v3, v4, and v5 are converged to the leader’s velocity v1. When the obstacles are encountered, the velocities of the slave UAVs are decreased to prevent collisions with the obstacles. Once all UAVs leave away from the obstacles, the velocities of all of them are increased to track the master velocity.Figure 18
Velocities and force information. (a) The velocities in thex-axis. (b) The velocities in the y-axis. (c) The velocities in the z-axis. (d) The master/slave force.
(a)
(b)
(c)
(d)The force information is shown in Figure18(d). fmx, fmy, and fmz are the components of the feedback force fm applied on the haptic device. fsx, fsy, and fsz are the components of the summed control force fs applied on the slave UAVs. When the rigid formation moves in a stable state, the forces are zero and the operator has no feeling about the environment. Once the obstacles are encountered, the forces increase immediately and the human operator at the master side can perceive the repulsive forces. Therefore, the force reflection makes the operator perceive the movement of the rigid formation better.Moreover, the maneuverability performance of formation and the perceptual sensitivity of the remote environment are shown in Figure19. The cross correlation is to assess the velocity tracking error and smaller values indicate greater velocity tracking. Just noticeable difference (JND) is to assess the force tracking accuracy and smaller values indicate greater perceptual sensitivity to the haptic cue. From the figure, we can see that compared with that of the traditional wave variable method, the performance of velocity tracking with local damping in the consensus formation controller is greater and the force tracking accuracy is asymptotically the same as that of the traditional wave variable method. The time-delayed bilateral aerial teleoperation system with local damping can be in a steady state to guarantee better velocity and force tracking in free motion or obstacle avoidance.Figure 19
Maneuverability performance and perceptual sensitivity.


### 5.2.2. Bilateral Aerial Teleoperation with Energy-Optimized Consensus Formation

To reduce the communication complexity and energy dissipation of the formation, the 3D top-down optimal persistent graph is used to simplify the formation topology. Moreover, the communication costs are calculated in the unoptimized and optimized mode, respectively. The parameters of the low-energy radio model are set asEelec = 50 nJ/bit, εfris_amp = 10 J/bit/m2, and l = 400 bit.The minimally rigid graphs of the five UAVs are shown in Figure7. There are a total of ten minimally rigid graphs. G1–G10 are the different minimally rigid graphs. In the optimized mode, a min-weighted rigid graph is selected from these minimally rigid graphs according to the weights of the communication edges. The weights of the ten candidate minimally rigid graphs at initial are shown in Figure 20. It can be seen from the figure that the initial min-weighted rigid graph is G2. On the basis of the min-weighted rigid graph, the top-down strategy is used to select the optimal persistent graph. The optimal persistent graphs of the consensus formation used in the optimized mode are shown in Figure 21. It can be seen that the minimally rigid graphs G6 and G9 are selected as the min-weighted rigid graph when the formation reaches the corner obstacle and the circular obstacle, respectively. In the absence of obstacles, the minimally rigid graph G2 is selected. Obviously, the communication links are 3n − 6 = 9, and this value is less than the full connections. Moreover, the weights of the communication are always minimal at different times. The average communication costs of the slave UAVs in fully connected, min-weighted rigid and optimal persistent graphs are shown in Figure 22. Obviously, compared with those with the fully connected graph and min-weighted rigid graph, the communication costs and links with the optimal persistent graph are minimal. The trajectories of the optimized consensus formation are shown in Figure 23. The topologies with optimal persistent graph are marked by the black dashed lines. Moreover, the 3D graphs of the experimental tests with one UAV, three UAVs, and five UAVs which implement the energy-optimized consensus formation scheme for the time-delayed bilateral teleoperation system are shown in Figure 24. Obviously, the optimized consensus formation can be achieved and the advantages are shown in Figure 25. Using the optimal persistent graphs in formation, the trajectories around the obstacles are smoother and the running time is significantly shortened. Due to the shorter running time, the consensus formation with the optimal persistent graph can fly farther.Figure 20
The weights of the initial candidate minimally rigid graphs. G2 is the min-weighted rigid graph at initial.Figure 21
The optimal persistent graphs of the formation used in the optimized mode.Figure 22
Average communication costs and links in different graphs.Figure 23
The consensus formation with 3D optimal rigid graph.Figure 24
3D Graphs of the experimental tests. (a) One UAV in the left view. (b) Three UAVs in the left view. (c) Five UAVs in the top view.
(a)
(b)
(c)Figure 25
The characteristics of consensus formation controller in different graphs. SN denotes the smoothness of the trajectories around the corner or circle obstacle and RT denotes runtime time around the corner or circle obstacle.


## 5.2.1. Bilateral Aerial Teleoperation with Consensus Formation

The ability of the consensus formation to adapt to the obstructed environment is evaluated by the human-in-the-loop simulations. During the movement, the UAVs should always maintain a rigid formation and avoid collisions with the neighbor UAVs and obstacles. In the consensus formation with the unoptimized topology, the communications are fully connected, that is, the numbers of communication edges are always ten regardless of changes in the neighborhood area.In Figure14, the UAVs in the bilateral aerial teleoperation reach the local minima in the corner obstacle and the consensus formation failed. In Figure 15, the consensus formation with flux-based collision avoidance function is shown. The topologies with full connections are marked by the black solid lines. The obstacles are considered as the electronic objects. During the time interval t ∈ [0 s, 150 s], the operator sends the desired velocity to the leader UAV. Under the action of the consensus formation controller, the UAVs are required to form a rigid formation, keep velocity consensus with the master, and avoid collisions with obstacles. When any of the slave UAVs encounters obstacles, the rigid formation breaks up to keep a safe distance from the obstacles. In the flux-conserved force field, the formation can escape the local minima from the corner obstacle. Until all UAVs leave away from the obstacles, the rigid formation is again formed. The relative distances between the slave UAVs are shown in Figure 16. It can be seen that the UAVs can avoid collisions with the neighbors and obstacles while keeping a desired rigid formation in the absence of obstacles. The yaw angles of the slave UAVs are shown in Figure 17.Figure 14
The consensus formation reaches the local minima.Figure 15
The trajectories of the consensus formation.Figure 16
The relative distances.Figure 17
The yaw angles.The velocity information is shown in Figures18(a)–18(c). rsx, rsy, and rsz are the components of the commanded velocity rs sent to the leader UAV. vix, viy, and viz are the components of the velocities of the slave UAVs for i = 1, …, 5. It can be seen that the velocities of the leader UAV v1 are converged to the commanded velocity rs generated by the master haptic device. The velocities of the followers v2, v3, v4, and v5 are converged to the leader’s velocity v1. When the obstacles are encountered, the velocities of the slave UAVs are decreased to prevent collisions with the obstacles. Once all UAVs leave away from the obstacles, the velocities of all of them are increased to track the master velocity.Figure 18
Velocities and force information. (a) The velocities in thex-axis. (b) The velocities in the y-axis. (c) The velocities in the z-axis. (d) The master/slave force.
(a)
(b)
(c)
(d)The force information is shown in Figure18(d). fmx, fmy, and fmz are the components of the feedback force fm applied on the haptic device. fsx, fsy, and fsz are the components of the summed control force fs applied on the slave UAVs. When the rigid formation moves in a stable state, the forces are zero and the operator has no feeling about the environment. Once the obstacles are encountered, the forces increase immediately and the human operator at the master side can perceive the repulsive forces. Therefore, the force reflection makes the operator perceive the movement of the rigid formation better.Moreover, the maneuverability performance of formation and the perceptual sensitivity of the remote environment are shown in Figure19. The cross correlation is to assess the velocity tracking error and smaller values indicate greater velocity tracking. Just noticeable difference (JND) is to assess the force tracking accuracy and smaller values indicate greater perceptual sensitivity to the haptic cue. From the figure, we can see that compared with that of the traditional wave variable method, the performance of velocity tracking with local damping in the consensus formation controller is greater and the force tracking accuracy is asymptotically the same as that of the traditional wave variable method. The time-delayed bilateral aerial teleoperation system with local damping can be in a steady state to guarantee better velocity and force tracking in free motion or obstacle avoidance.Figure 19
Maneuverability performance and perceptual sensitivity.


## 5.2.2. Bilateral Aerial Teleoperation with Energy-Optimized Consensus Formation

To reduce the communication complexity and energy dissipation of the formation, the 3D top-down optimal persistent graph is used to simplify the formation topology. Moreover, the communication costs are calculated in the unoptimized and optimized mode, respectively. The parameters of the low-energy radio model are set asEelec = 50 nJ/bit, εfris_amp = 10 J/bit/m2, and l = 400 bit.The minimally rigid graphs of the five UAVs are shown in Figure7. There are a total of ten minimally rigid graphs. G1–G10 are the different minimally rigid graphs. In the optimized mode, a min-weighted rigid graph is selected from these minimally rigid graphs according to the weights of the communication edges. The weights of the ten candidate minimally rigid graphs at initial are shown in Figure 20. It can be seen from the figure that the initial min-weighted rigid graph is G2. On the basis of the min-weighted rigid graph, the top-down strategy is used to select the optimal persistent graph. The optimal persistent graphs of the consensus formation used in the optimized mode are shown in Figure 21. It can be seen that the minimally rigid graphs G6 and G9 are selected as the min-weighted rigid graph when the formation reaches the corner obstacle and the circular obstacle, respectively. In the absence of obstacles, the minimally rigid graph G2 is selected. Obviously, the communication links are 3n − 6 = 9, and this value is less than the full connections. Moreover, the weights of the communication are always minimal at different times. The average communication costs of the slave UAVs in fully connected, min-weighted rigid and optimal persistent graphs are shown in Figure 22. Obviously, compared with those with the fully connected graph and min-weighted rigid graph, the communication costs and links with the optimal persistent graph are minimal. The trajectories of the optimized consensus formation are shown in Figure 23. The topologies with optimal persistent graph are marked by the black dashed lines. Moreover, the 3D graphs of the experimental tests with one UAV, three UAVs, and five UAVs which implement the energy-optimized consensus formation scheme for the time-delayed bilateral teleoperation system are shown in Figure 24. Obviously, the optimized consensus formation can be achieved and the advantages are shown in Figure 25. Using the optimal persistent graphs in formation, the trajectories around the obstacles are smoother and the running time is significantly shortened. Due to the shorter running time, the consensus formation with the optimal persistent graph can fly farther.Figure 20
The weights of the initial candidate minimally rigid graphs. G2 is the min-weighted rigid graph at initial.Figure 21
The optimal persistent graphs of the formation used in the optimized mode.Figure 22
Average communication costs and links in different graphs.Figure 23
The consensus formation with 3D optimal rigid graph.Figure 24
3D Graphs of the experimental tests. (a) One UAV in the left view. (b) Three UAVs in the left view. (c) Five UAVs in the top view.
(a)
(b)
(c)Figure 25
The characteristics of consensus formation controller in different graphs. SN denotes the smoothness of the trajectories around the corner or circle obstacle and RT denotes runtime time around the corner or circle obstacle.


## 6. Conclusions and Future Work

In this paper, we present an energy-optimized consensus formation control scheme for the time-delayed bilateral teleoperation system that enables the operator to manipulate a group of UAVs in the obstructed environment. The damping injections are independently distributed on the master and each slave UAV to deal with the time-varying delays. The stability criteria of the time-delayed teleoperation system are obtained by using the Lyapunov function. The velocity commands generated by the master device are transmitted to the leader UAV to move the formation. The followers are controlled by the consensus formation controller to form a rigid formation, keep velocity consensus with the leader, and avoid obstacles with the flux-conserved force field. Moreover, a top-down strategy of 3D optimal persistent graph is proposed to decrease the communication complexity and energy dissipation of the formation. The proposed framework is evaluated by the human-in-the-loop simulations of UAVs in an obstructed environment. Experimental results show that all UAVs can keep a desired formation in the time-delayed bilateral aerial teleoperation and the operator can perceive the movement of the rigid formation. The good velocity and force tracking are guaranteed in free motion or obstacle avoidance. In addition, with the 3D optimal persistent graph, the communication links and costs of are minimal. The smoothness, runtime, and trajectory of the formation are significantly improved.In the future work, we will research on the bilateral aerial teleoperation of multiple UAVs in the real cooperative tasks, such as cooperative transportation or firefighting.



---

*Source: 2857674-2018-05-22.xml*