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.

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 [2325]. 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.

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].

Let G = (V, E, A) be a weighted graph, where V = {1, 2, …, n} is the set of vertexes, is the set of edges, and A = [aij] is the weighted adjacency matrix. The adjacency elements associated with the edges are positive, that is, . Moreover, the elements aii = 0 for all . The set of neighbors of vertex is denoted by . The Laplacian of graph G is denoted by L = [lij] with , 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 . The L2 norm is .

2.2. Rigid Graph

As shown in [30], let qi(t) be the trajectory of vertex 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 . Or else, it is called flexible.

For a graph , the rigidity matrix M is defined as where each row corresponds to an edge and the columns corresponds to the coordinates of the vertexes.

Lemma 1 [30, 31]. Let M be the rigidity matrix of a framework of n vertexes in R c. A framework G = (V, E, A) with n (n ≥ 2) vertexes in R c is infinitesimally rigid if and only if rank (M) = c n − c(c + 1)/2.

Lemma 2 [30, 31]. An infinitesimally rigid framework G = (V, E, A) with n (n ≥ 2) vertexes in R c is minimally rigid if and only of it has c n − c(c + 1)/2 edges.

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.

2.3. 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 is the joint position, velocity, and acceleration. is the positive definite inertia matrix. is the Coriolis and centripetal matrix. For the sake of simplicity, the enclosed parameter(s) of are attached in Appendix A. is the human force applied on the device with the compensated gravitational force. 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., ), except for those which are time delayed (e.g., 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 endpoints of the X-shaped frame. Let ϖi = [xi, yi, zi, φi, θi, ψi] be the generalized coordinates where qi = [xi, yi, zi]TR3 denotes the absolute position of the ith UAV in the inertial frame and Θi = [φi, θi, ψi]TR3 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 mi is the mass. is the moment of inertia matrix. 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 is the Coriolis matrix.

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, . Here, is the translational force applied on the ith UAV with the compensated gravitational force. is the roll, pitch, and yaw torques applied on the ith UAV. is the force/torque caused by the interaction of the ith UAV with the environment. Finally, the dynamics of the ith UAV (iV) can be defined as

Besides, let vi ∈ R3 be the linear velocity and the variable ρi ∈ R3 is defined as where γ → 0+, 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 .

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. is skew symmetric.

Property 2. such that .

Property 3. such that .

Assumption 1. The human operator and the environment are passive for i ∈ V, that is,

Assumption 2. The time-varying delays Tm1(t), T1m(t), Tij(t), and Tji(t) are bounded, that is, 0 ≤ Tm1(t) ≤ , 0 ≤ T1m(t) ≤ , 0 ≤ Tij(t) ≤ , and 0 ≤ Tji(t) ≤ , where is the maximum time-varying delay from the master to the leader UAV, is the maximum time-varying delay from the leader UAV to the master, is the maximum time-varying delay from the ith UAV to the jth UAV, and 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 time-varying delay , there exists a positive-definite matrix Γ such that the following inequality holds:

Lemma 6. For any vector function x(∙) and y(∙) and any time-varying delay , there exists a constant α such that the following inequality holds: where 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 Figure 3.

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 where kn is a positive constant.

The commanded velocity for the leader UAV is 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 coupling , 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 . 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 [14].

First, let be the local control loop with a storage function ; the coupling is passive with , that is,

Further, by introducing a scaling into this strategy, let . By properly choosing the parameters μ and λ, it is possible to neglect the contribution of and make the commanded velocity rm proportional to the master scaled position with .

The new energy storage is chosen as

The time derivative of is

Therefore, the master device is passive with respect to the coupling .

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 Figure 4. Therefore, the consensus formation controller is designed as 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 (iV) 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.

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 is the moment of inertia matrix in yaw coordinates and is the control torque applied on the ith UAV in yaw coordinates. 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 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 such that the following inequalities hold: then the overall teleoperation system is stable and the variables , , , , , , and are bounded.

Proof 1. The total control force 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 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 is divided into two parts, that is, where Firstly, the relationship between the master and the leader UAV is investigated. The total forces 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 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, can be rewritten as Noticing that , 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 .
Similarly, the force-velocity coupling is passive, Integrating equation , one has Let . When the parameters μ and λ are properly selected, variable σ is proportional to the variable ρ1 with .
Further, the new energy storage is selected as The time derivative of is Therefore, the fourth term of (31) is passive with respect to coupling .
In conclusion, to ensure the stability of the master and the leader, 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 and are the momentum and the inertia matrix, respectively, of the ith UAV. With Assumption 1, is greater than zero.
With (23), the time derivation of is Considering the following total scaled energy function, where , , R is the relative distance matrix, and is the standard Kronecker product. Noticing that , one has By integrating the above equation from zero to t with Lemma 5, we obtain where δi, εi > 0. Using the fact that H > 0, Defining and , Then H(0) can be rewritten as . There exists for such that and . Therefore, setting the damping injection βi, for any arbitrary . Each slave UAV ensures that and , which in turn implies that , for all .
Combining (21), (36), and (45), one has , , , . 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. Let Eiobs denote the force field and C denote any closed curves around the obstacle in the 3D space. The flux is defined as 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 Figure 5, 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 where qi ∈ R3 is the position of the ith UAV and 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.

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.

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

Initialize the complete graph G = (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 matrix M of the complete graph according to the sequence of edges;
Initialize matrix Mc as the first row of M;
for
if rank
  Add the next row of M to Mc to form a new matrix M’;
  if M’ is full rank
  Record the edge set in M 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 neighbors Ni in the consensus formation controller.

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 and . 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).

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.

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 = [aij] and the connectivity of each vertex d(i);
Establish the subspaces according to the contour map and the number of UAVs in subspaces is ni;
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;
forspace = 1 : subspaces
fori = 1 : ni
  ifd(j) = 2, d(k) = 2, and aij = aik = 1,
   Establish the directed edge Vij, 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 And to receive this message, the energy dissipation is where the electronics energy depends on many factors, such as digital coding, modulation, and spreading of the signal. The amplifier energy 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 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-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.

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 are qm = (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,  = 0.7 s,  =  =  = 1.2 s. The desired relative distance matrix of the consensus formation is . The turning procedures of the control parameters are as follows: (1)Setting kn 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)Setting ks 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)Setting km 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)Setting kv 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: Combined with the previously selected coefficients kc and kv, the damping coefficient βi in different modes can be obtained by Theorem 1 (45).(6)Setting the factor of obstacle repulsive force kobs 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 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 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. 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). 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.

The force information is shown in Figure 18(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 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.

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 as  = 50 nJ/bit,  = 10 J/bit/m2, 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, 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.

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.

Appendix

A. Enclosed Parameter(s) for (2)

The physical map and kinematic model diagram are shown in Figure 2(a). The workspace of the PHANTOM Omni haptic interface point (HIP) is −210 mm ≤ xm ≤ 210 mm, −110 mm ≤ ym ≤ 205 mm, −85 mm ≤ zm ≤ 130 mm. The kinematic and dynamic models will be defined in the sequel. The enclosed parameters of , will be given in the dynamic model.

A.1. Kinematic Model

All joints of the master haptic device PHANTOM Omni are rotational. Joints θ1, θ2, and θ3 give the translational movements for the end effector while joints θ4, θ5, and θ6 provide the rotational movements. The device is only actuated in the first three joints, which gives the ability to generate haptic forces in the translational directions. After employing the Denavit Hartenberg (D–H) conversion on the manipulator, the transformation matrix can be found and the kinematics can be expressed as [37] where qm = [xm, ym, zm] is the joint position, θ = [θ1, θ2, θ3]T is the vector of joints, and L1 = L2 = 133.35 mm, which means the length of link 1 and link 2. ky = 23.35 mm and kz = −168.35 mm represent workspace transformation offset between the origin of the haptic device and the first joint.

Further, the time derivative of (51) is where Jm = [Jm11, Jm12, Jm13; Jm21, Jm22, Jm23; Jm31, Jm32, Jm33] is the Jacobian matrix. The parameters of the Jacobian matrix are expressed as

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] where M is the inertia matrix, C is the Coriolis and centrifugal forces, represents the gravity, and τ = JmTFm 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. where

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. where

Gravity forces depend on the joint angles, lengths, and mass of each joint

The variable k1 to k10 is related to inertia, gravity effect, mass, and length of joints. The values are given by

Therefore, the dynamics of the master haptic device PHANTOM Omni can be rewritten as where are the joint position, velocity, and acceleration. is the positive-definite inertia matrix. is the Coriolis and centripetal matrix. Fm = fm + fh, fh is the human force applied on the device with the compensated gravitational force, and fm 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 = [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 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 where mi is the mass of each UAV. The rotational kinetic energy is where matrix Ji is the moment of inertia matrix. The potential energy is Ui = migzi, where 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 is obtained from the Euler–Lagrange equations with external generalized force where fi 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 fi = RiFit where Ri is the rotation matrix and Fit = [0, 0, Ui1]T is the translational force in the inertial frame. Ui1 is the summation of thrust moments where Fi1, Fi2, Fi3, and Fi4 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 wji is the angular velocity of the motor j on the ith UAV.

Then the generalized moments for Eular angles are 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

Conflicts of Interest

The authors declare that there is no conflict of interests regarding the publication of this article.

Acknowledgments

The research reported in this paper was carried out at the State Key Laboratory of Bioelectronics, Jiangsu Key Lab of Remote Measurement and Control, School of Instrument Science and Engineering, Southeast University, Nanjing, Jiangsu, China. This work is supported in part by the National Natural Science Foundation of Jiangsu Province under Grant 61375076 and Research & Innovation Program for Graduate Student in Universities of Jiangsu Province under Grant KYLX16_0192. The authors thank all the members of the Robotic Sensor and Control Lab for their great support.