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

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.


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 timevarying 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-offreedom (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][24][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.Section 2 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.

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].
Let G = (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 = [a ij ] is the weighted adjacency matrix.The adjacency elements associated with the edges are positive, that is, i, j ∈ V⇔a ij > 0.Moreover, the elements a ii = 0 for all i ∈ V.The set of neighbors of vertex i is denoted by I n is an identity matrix of size n, 1 n is a column vector of size n with all elements equal to one.The 2.2.Rigid Graph.As shown in [30], let q i (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 ||q i (t) − q j (t)|| are constant for all i, j ∈ E. Or else, it is called flexible.
For a graph G ⊆ R c c ∈ 2, 3 , the rigidity matrix M is defined as where each row 0 ⋯ 0q T i − q T j 0 ⋯ 0q T j − q T i 0 ⋯ 0 corresponds to an edge i, j ∈ E and the columns corresponds to the coordinates of the vertexes.
Lemma 3 [30,31].If every edge of the framework G = (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 framework G = (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 Figure 1.

System Dynamics.
The master device of the bilateral aerial teleoperation system is shown in Figure 2(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]: where q m , q m , q m ∈ R 3 is the joint position, velocity, and acceleration.M m q m ∈ R 3×3 is the positive definite inertia matrix.C m q m , q m ∈ R 3×3 is the Coriolis and centripetal matrix.For the sake of simplicity, the enclosed parameter(s) of M m q m , C m q m , q m are attached in Appendix A. f h ∈ R 3 is the human force applied on the device with the compensated gravitational force.f m ∈ R 3 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., q m ≡ q m t ), except for those which are time delayed (e.g., q m t − T m1 ≡ q m t − T m1 t for a variable time delay).As shown in Figure 2(b), the slave ith UAV is a rigid body that is actuated by a thrust force generated by four rotors on The four agents in the plane are (0, 0), (3, 0), (4,4), and (0, 5).
3 International Journal of Aerospace Engineering the endpoints of the X-shaped frame.Let ϖ i = [x i , y i , z i , φ i , θ i , ψ i ] be the generalized coordinates where q i = [x i , y i , z i ] T ∈ R 3 denotes the absolute position of the ith UAV in the inertial frame and Θ i = [φ i , θ i , ψ i ] T ∈ R 3 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].
where m i is the mass.J i ∈ R 3 is the moment of inertia matrix.f i , τ i ∈ R 6 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 where Then (4) can be rewritten as The dynamics of the ith UAV can be yielded as 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, f i , τ i = f s i + f e i , τ s i + τ e i .Here, f s i is the translational force applied on the ith UAV with the compensated gravitational force.τ s i is the roll, pitch, and yaw torques applied on the ith UAV.f e i , τ e i 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 Besides, let v i ∈ R 3 be the linear velocity and the variable ρ i ∈ R 3 is defined as where γ → 0 + , v i ∈ R 3 is the acceleration of the ith UAV, and ρ i ∈ R 3 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].
Assumption 1.The human operator and the environment are passive for i ∈ V, that is, Assumption 2. The time-varying delays T m1 (t), T 1m (t), T ij (t), and T ji (t) are bounded, that is, 0 ≤ T m1 (t) ≤ T m1 , 0 ≤ T 1m (t) ≤ T 1m , 0 ≤ T ij (t) ≤ T ij , and 0 ≤ T ji (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 function a(•) and b(•) and any timevarying delay 0 ≤ T t ≤ T, there exists a positive-definite matrix Γ such that the following inequality holds: Lemma 6.For any vector function x(•) and y(•) and any timevarying delay 0 ≤ T t ≤ ∞, there exists a constant α such that the following inequality holds: where • 2 2 is the L 2 norm of the signal.

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 Figure 3.

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 where k n is a positive constant.
The commanded velocity for the leader UAV is where T m1 is the time-varying delay from the master device to the leader UAV.r s ∈ R 3 is the commanded velocity for the leader UAV.

International Journal of Aerospace Engineering
It is well known that the master device is passive with respect to the force-velocity coupling f m , 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 f m , r m .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 f m , r m [14].
First, let r = q m + λq m , λ ∈ R + be the local control loop with a storage function Further, by introducing a scaling into this strategy, let r m = μr = μq m + μλq m , μ, λ ∈ R + .By properly choosing the parameters μ and λ, it is possible to neglect the contribution of q m and make the commanded velocity r m proportional to the master scaled position with r m ≈ k n q m , k n = μλ.
The new energy storage is chosen as The time derivative of κ m is Therefore, the master device is passive with respect to the coupling f m , r m .

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 Figure 4. Therefore, the consensus formation controller is designed as where T 1m is the time-varying delay from the leader UAV to the master, T m1 is the time-varying delay from the master to the leader UAV, and T ij is the time-varying delay from the ith UAV (i ∈ V) to the jth UAV.k m , k s , k v , and k c are the positive constants.α m , α s , and β i are the damping coefficients.b i = 1 when the 1th UAV is connected to the master, or else, b i = 0. d ij is the desired relative distance between the ith UAV and the jth UAV.N i 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.k m and k s are the positive constants which represent the proportional gains to achieve velocity coordination between the master device and the leader UAV.k v and k c 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 k m , k s , k v , and k c , 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.
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 where J ψ i is the moment of inertia matrix in yaw coordinates and τ sψ i is the control torque applied on the ith UAV in yaw coordinates.τ eψ i 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, for each slave UAV ensures that 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: then the overall teleoperation system is stable and the variables , and v i are bounded.
Proof 1.The total control force f s i 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 f s i 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 force f s i is divided into two parts, that is, where Firstly, the relationship between the master and the leader UAV is investigated.The total forces ∑ n i=1 f s# i which are applied on the leader UAV are 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
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 With ( 18), (23), and Assumption 2, H # can be rewritten as 7 International Journal of Aerospace Engineering q m ξ dξ, one has With Lemma 4, the following inequalities are obtained: With ( 26), ( 27), ( 28), (29), and ( 30), one has With Assumption 1, the third term of ( 31) is passive with respect to the force-velocity coupling f h , q m .Similarly, the force-velocity coupling f e 1 , v 1 is passive, Integrating equation When the parameters μ and λ are properly selected, variable σ is proportional to the variable ρ 1 with σ ≈ μγρ 1 .
Further, the new energy storage is selected as The time derivative of κ s is Therefore, the fourth term of ( 31) is passive with respect to coupling f e 1 , ρ 1 .In conclusion, to ensure the stability of the master and the leader, H # should be negative which is equivalent to the following inequalities: Setting α m , α s in (36) for the master and the leader UAVs implies that For the follower UAVs, the following Lyapunov function for the ith UAV (i ∈ V) is constructed as where p i ∈ R 3 and M i ∈ R 3×3 are the momentum and the inertia matrix, respectively, of the ith UAV.With Assumption 1, H * i is greater than zero.With (23), the time derivation of H * i is Considering the following total scaled energy function, International Journal of Aerospace Engineering where q = q 1 , … , q n , v = v 1 , … , v n , R is the relative distance matrix, and ⊗ is the standard Kronecker product.
By integrating the above equation from zero to t with Lemma 5, we obtain

42
where δ i , ε i > 0. Using the fact that H * > 0, 43 Then H(0) can be rewritten as H 0 ≥ 1 T n ΩΓ.There exists Therefore, setting the damping injection β i , for any arbitrary δ i , ε i > 0. Each slave UAV ensures that v i ∈ L 2 and H * ∈ L ∞ , which in turn implies that q i − q j t − Combining (21), (36), and (45), one has The proof is completed.

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.Let E i obs denote the force field and C denote any closed curves around the obstacle in the 3D space.The flux is defined as where W i obs is the flux of obstacle i for i ∈ {1, …, n obs } 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 Figure 5, for any curve C surrounding obstacle i, the value of W i obs remains the same.The repulsive force of the ith UAV from the obstacles is where q i ∈ R 3 is the position of the ith UAV and x oi ∈ R 3 is the position of obstacle i. k obs is a positive coefficient.For any curve surrounding an electric point, the value of W i obs is always same.
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 Figure 6, 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.International Journal of Aerospace Engineering

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.

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 minweighted 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 minweighted rigid graph is shown in Algorithm 1.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.
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 neighbors N i in the consensus formation controller.

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.
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 Figure 9.
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 Figure 7 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).International Journal of Aerospace Engineering 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 Algorithm 2.   12 International Journal of Aerospace Engineering 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

Energy Dissipation of
And to receive this message, the energy dissipation is where the electronics energy E elec depends on many factors, such as digital coding, modulation, and spreading of the signal.The amplifier energy ε f ris amp depends on the acceptable intensity and bit-error rate.

Human-in-the-Loop Simulations
5.1.Experimental Setup.The experimental setup for the bilateral teleoperation system is shown in Figure 11.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.The commercial Sensable PHANTOM Omni is selected as the master device.The haptic device has six degrees-offreedom (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 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 matrix A = [a ij ] and the connectivity of each vertex d − (i); Establish the subspaces according to the contour map and the number of UAVs in subspaces is n i ; 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 : , and a ij = a ik = 1, Establish the directed edge V ij , V ik and d − (k) = 2; if for any vertex i, the out-degrees d The 3D optimal persistent graph is generated; end end end end Algorithm 2: Generation of the 3D optimal persistent graph in obstructed environment.13 International Journal of Aerospace Engineering 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.
(2) Setting k s to achieve master-slave velocity tracking.In the unilateral teleoperation case with no time delays, control parameters k m , T m , and T s are all equal to zero.When k s is set to 1 (Ns/m), the speed of the UAV can completely track that of the master.
(3) Setting k m to achieve master-slave force tracking.In the bilateral teleoperation case with no time delays, the time delays T m , T s are equal to zero.When the control parameter k m is set to 1 (Ns/m), the master feedback force f m can completely track the slave control force f s .
(4) Setting k v and k c 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 k v 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 k c 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 k n , k s , k m , max{T m1 }, max{T 1m }, max{T ij }, and max{T ji }, the damping coefficient α m and α s can be obtained by Theorem 1.International Journal of Aerospace Engineering Also, let any positive constants δ i = ε i = 1.The Laplacian matrices in the fully connected, G2, G6, and G9 modes are computed as follows: , , , Combined with the previously selected coefficients k c and k v , the damping coefficient β i in different modes can be obtained by Theorem 1 (45).
(6) Setting the factor of obstacle repulsive force k obs 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: , and G9 modes, respectively.

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 Figure 14, 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.15 International Journal of Aerospace Engineering 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.
The velocity information is shown in Figures 18(a)-18(c).r sx , r sy , and r sz are the components of the commanded velocity r s sent to the leader UAV.v ix , v iy , and v iz 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 v 1 are converged to the commanded velocity r s generated by the master haptic device.The velocities of the followers v 2 , v 3 , v 4 , and v 5 are converged to the leader's velocity v 1 .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.
The force information is shown in Figure 18(d).f mx , f my , and f mz are the components of the feedback force f m applied on the haptic device.f sx , f sy , and f sz are the components of the summed control force f s 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 Figure 19.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.

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 as E elec = 50 nJ/bit, ε f ris amp = 10 J/bit/m 2 , and l = 400 bit.
The minimally rigid graphs of the five UAVs are shown in Figure 7.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, minweighted 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 16 International Journal of Aerospace Engineering 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.

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-   The master/slave force.17 International Journal of Aerospace Engineering 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.A.2. Dynamic Model.Without loss of generality, the master haptic device is a fully actuated system which can be described by the Euler-Lagrange equation [38] M θ θ + C θ, θ θ + g θ, θ = τ, 54 where M is the inertia matrix, C is the Coriolis and centrifugal forces, g represents the gravity, and τ = J m T F m is the vector of torques acting on the joints; θ is The inertia matrix is found by examining the moment of inertia of each link of the manipulator according to the joint angle, length, and mass.The Coriolis and centrifugal effects appear for objects moving with a rotating frame of reference.The Coriolis effect depends on both mass and velocity but is independent of the position.The centrifugal force depends only of the position and the mass and is always oriented away from the axis of rotation.Therefore, the dynamics of the master haptic device PHANTOM Omni can be rewritten as M m q m q m + C m q m , q m q m = F m = f m + f h , 62 where q m , q m , q m ∈ R 3 are the joint position, velocity, and acceleration.M m q m ∈ R 3×3 is the positive-definite inertia matrix.C m q m , q m ∈ R 3×3 is the Coriolis and centripetal 20 International Journal of Aerospace Engineering matrix.F m = f m + f h , f h is the human force applied on the device with the compensated gravitational force, and f m is the feedback force resulting from the interaction between the master device and the slave side.

B. Derivation of the ith UAV's Dynamics
Let ϖ i = [x i , y i , z i , φ i , θ i , ψ i ] be the generalized coordinates where q i = [x i , y i , z i ] T ∈ R 3 denotes the absolute position of the ith UAV in the inertial frame and Θ i = [φ i , θ i , ψ i ] T ∈ R 3 denotes the three Euler angles representing the roll, pitch, and yaw, respectively.The dynamic model of each UAV is derived using the Euler-Lagrange approach.
The dynamic model of each UAV can be separated into the translational and rotational coordinates.The translational kinetic energy of each UAV is T trans i = 1/2m i q T i q i where m i is the mass of each UAV.The rotational kinetic energy is T rot i = 1/2Θ T i J i Θ i where matrix J i is the moment of inertia matrix.The potential energy is U i = m i gz i , where g is the acceleration due to gravity.Then the Lagrange model of the ith UAV can be defined as follows [33]: The dynamics of the ith UAV obtained from the Euler-Lagrange equations with external generalized force where f i is the translational force applied to the ith UAV due to the throttle force in the body frame, τ i is the roll, pitch, and yaw moments.The translational force f i = R i F i t where R i is the rotation matrix and F i t = [0, 0, U i 1 ] T is the translational force in the inertial frame.U i 1 is the summation of thrust moments where F i 1 , F i 2 , F i 3 , and F i 4 are the thrust moments of each propeller.The thrust moment is defined as where b is the thrust factor for each propeller of each UAV and w ji is the angular velocity of the motor j on the ith UAV.
Then the generalized moments for Eular angles are , 67 where l is the distance from the motor to the center of gravity and d is the drag factor.Thus, the control inputs of the ith UAV is given by The Euler-Lagrange equation can be partitioned into the dynamics for the translational coordinates and rotational coordinates.Thus, the final model of the ith UAV can be represented as

Figure 2 :
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.

Figure 4 :
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.

Figure 6 :
Figure 6: No local minima in flux-conserved force field.

Figure 5 :
Figure 5: The flux of the force field is conserved for any curves.

Figure 10 :
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.

Figure 11 :
Figure 11: The experimental setup for the bilateral teleoperation system.

Figure 12 :
Figure 12: The framework of the experimental setup.

Figure 13 :
Figure 13: The time-varying delays for the semiexperiments.

Figure 14 :
Figure 14: The consensus formation reaches the local minima.

Figure 15 :
Figure 15: The trajectories of the consensus formation.

Figure 18 :
Figure 18: Velocities and force information.(a) The velocities in the x-axis.(b) The velocities in the y-axis.(c) The velocities in the z-axis.(d) The master/slave force.

Figure 20 :Figure 21 :Figure 22 :Figure 23 :
Figure 20: The weights of the initial candidate minimally rigid graphs.G2 is the min-weighted rigid graph at initial.
(Yan et al., 2014)n.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