applied sciences Article Cooperative Multi-UAV Collision Avoidance Based on a Complex Network Yang Huang, Jun Tang * and Songyang Lao * College of Systems Engineering, National University of Defence Technology, Changsha 410073, China; huangyang13@nudt.edu.cn * Correspondence: jun.tang@e-campus.uab.cat (J.T.); laosongyang@nudt.edu.cn (S.L.)  Received: 15 August 2019; Accepted: 16 September 2019; Published: 20 September 2019  Abstract: The conflict resolution problem in cooperative unmanned aerial vehicle (UAV) clusters sharing a three-dimensional airspace with increasing air traffic density is very important. This paper innovatively solves this problem by employing the complex network (CN) algorithm. The proposed approach allows a UAV to perform only one maneuver—that of the flight level change. The novel UAV conflict resolution is divided into two steps, corresponding to the key node selection (KS) algorithm based on the node contraction method and the sense selection (SS) algorithm based on an objective function. The efficiency of the cooperative multi-UAV collision avoidance (CA) system improved a lot due to the simple two-step collision avoidance logic. The paper compares the difference between random selection and the use of the node contraction method to select key nodes. Experiments showed that using the node contraction method to select key nodes can make the collision avoidance effect of UAVs better. The CA maneuver was validated with quantitative simulation experiments, demonstrating advantages such as minimal cost when considering the robustness of the global traffic situation, as well as significant real-time and high efficiency. The CN algorithm requires a relatively small computing time that renders the approach highly suitable for solving real-life operational situations. Keywords: Multi-UAV; collision avoidance; complex network; key nodes; robustness; connected component 1. Introduction Collision is an inherent problem in unmanned aerial vehicle (UAV) systems [1–6]. Path planning is important in performing different tasks. Therefore, the demand for the collision avoidance (CA) of UAVs has emerged. With increasing UAV cluster densities, it has become more difficult to perform multi-UAV collaborative tasks, especially if collision avoidance requirement cannot be met. Therefore, the effective planning of flight trajectories and collision avoidance among UAV flight trajectories has become critical in UAV clusters for safe and smooth operation [7]. Various conflict detection and resolution (CDR) theories with different techniques have been proposed to solve this problem. In conventional aviation, several widely used CA systems exist. For example, the traffic alert and CA system (TCAS) is recognized worldwide as an airborne collision avoidance system, and it is the last line of the CDR approach. Technological development has enabled the increase in air traffic density to be mitigated [8,9]. Many studies regarding UAV collision avoidance have been performed in recent years [10–12]. CDR and the optimization of different resolutions have been discussed widely by many researchers and practitioners. These methods can be classified into two primary types of algorithms: Geometric and path planning. Geometric algorithms focus on analyzing the relationship between UAVs and intruders within a specific space to implement passive CA operations. Meanwhile, path planning algorithms focus on obtaining a path for an optimal Appl. Sci. 2019, 9, 3943; doi:10.3390/app9193943 www.mdpi.com/journal/applsci Appl. Sci. 2019, 9, 3943 2 of 18 connection between the start and end points with minimal security constraints [13–15]. Dimensionality, direction selection, the number of participants, and resolutions are four key elements that should be considered for a UAV’s CDR. A detailed summary of publications regarding the different methods of CA is provided in Table 1. Geometric algorithms exploit state information such as the location, velocity, and the heading of UAVs and intruders to geometrically perform the CA operation. Park, Oh, and Tahk (2008) [16] demonstrated a typical example of this approach, in which a single resolution maneuvering logic called ‘vector sharing resolution’ was proposed. This method calculates the subtraction of two UAVs’ movement vectors in a two-dimensional space and adjusts the trajectories of the UAVs based on the shortest distance vector obtained. Another instant of this method can be found in [17] by Chakravarthy et al. In contrast, this paper solves the problem of obstacles being irregular objects. A collision cone approach is employed by which the irregularly shaped moving objects can be modelled through general quadric surfaces and dynamic inversion-based avoidance strategies can be derived. The limitations of this type of method is the need for information from the intruder UAV, as well as sensitivity to noise from the input data of the sensors. Path planning algorithms include sampling-based, artificial heuristic, numerical optimization, and decoupled path planning approaches. Each method exhibits its own strength over others in certain aspects. The philosophy of sampling-based approaches is transforming the path-planning problem to obtain a suitable problem from a limited quantity of candidates that is a continuous state space. Lin and Saripalli (2017) [18] presented a method based on the closed-loop rapidly exploring random tree algorithm and three variations of it. Artificial heuristic approaches are highly advantageous during the optimization process, and typical methods include the genetic algorithm [19], particle swarm optimization [20], and artificial bee colony optimization [21]. Numerical optimization approaches rely on some geometric trajectory calculations and aim to obtain an optimal path with respect to a specific objective in the solution sense. The relevant methods include mixed integer linear programming [22], nonlinear programming [23], Dynamic programming [24], quadratic programming [25], and Pontryagin’s minimum principle [26]. In decoupled path planning approaches, the basic principle involves obtaining an optimal trajectory under the specific constraints in the first phase; subsequently, the trajectory has to be steered out of the dangerous region. The A*(A-Star) search algorithm was applied to obtain a solution in [27]. The Coronoid graph can be used to generate an initial path that can be subsequently refined by curve fitting with curvature and other constrains [28]. The concepts of space–time reachability of aircraft and space–time potential conflict space were proposed in [29] to consider the uncertainty of the pilot intent. The conflict resolution problem was modeled as a constraint optimization problem in [30,31]. Based on the Single European Sky Air Traffic Management Research and Next Generation Air Transportation System initiatives, [32] proposed a new separation management service to shift the centralized air traffic control interventions to more efficient decentralized tactical operations based on a surrounding traffic analysis tool. The rapid development of real-time communication technology has made real-time communication between UAVs possible. Therefore, this paper was able to consider using the real-time information of UAV to prevent collision. It is a great idea to use a complex network to analyze the nature of a system. Networks are everywhere, including the internet, electric power grids, supply chains, urban road networks, and the world trade web [33]. The work already done has mainly been about preventing networks from being destroyed [33]. The collision avoidance network is exactly the opposite of most networks, and it is necessary to destroy the network as soon as possible. The collision avoidance based on a complex network (CACN) model grounds the generation of optimal direction changes of the objective function, which integrates the number of intruder UAVs of designated direction, the robustness of the network, and the connected components of the network. The CACN strategy can reduce the domino effect to a certain extent thus improve computational consuming time. In particular, the generated points by the CACN model are recorded and can be further analyzed. Appl. Sci. 2019, 9, 3943 3 of 18 Table 1. Summary of publications addressing the CDR problem. Dimension Maneuvers Multiple Resolution Geometric Algorithm Path Planning Algorithm Publication 3D T/V/S/C Pairwise Sampling-Based Numerical Decoupled Path Global Artificial Heuristic CPA CCA Path Planning Optimization Planning Approaches Approaches Approaches Approaches Park, Oh, and Tahk (2008) X C(T,S) X X Chakravarthy, Animesh, and Ghose (2012) X C(T,V,S) X X Lin and Saripalli (2017) X C(T,V) X X Pehlivanoglu (2012) X C(T,V) X X Karimi and Pourtakdoust (2013) X C(T,V) X X Xu, Duan and Liu (2010) C(T,V) X X Campbell, Bragg and Neogi (2013) X C(T,V) X X Tisdale, Kim and Hedrick (2009) V X X Jorris and Cobb (2008) C(V,S) X X Grancharova et al. (2014) X C(T,V,S) X X Sridhar, Ng and Chen (2011) V X X Kim, Gu and Postlethwaite (2008) X V X X Dai and Cocharan (2010) V X X Hao, Cheng and Zhang (2018) X X Yang, Yin and Shen. (2018) X X X Yang et al. (2019) X X Radanovic et al. (2018) X X X Note: T = Turns; V = Vertical maneuvers; S = Speed; C() = Combined maneuvers. Appl. Sci. 2019, 9, 3943 4 of 18 The layout of this paper is as follows. Section 2 describes the CA system architecture; Section 3 illustrates the key node selection algorithm and the sense selection algorithm; Section 4 shows the simulation experiments and analyzes the experimental results; finally, conclusions and prospects for future research are summarized in Section 5. 2. The Proposed Collision Avoidance System Architecture Many algorithms in the collision avoidance system support its operation. The core step of the system is to process the information of each trajectory input, which is the four-dimensional (three-dimensional (3D) position + time) information of each UAV, and an analysis on whether any threat is present in the next time interval. If a threat is detected, the system immediately reacts to calculate a suitable CA strategy. The instructions are subsequently communicated to each UAV to synergistically resolve the current threat. The distributed dynamic optimization approach (DDOA) is used to generate suitable direction changes for the UAVs such that the problem of CA can be mitigated. The core of this approach consists of two parts: The conflict detection (CD) module is used to forecast the threats, and the conflict resolution (CR) module is applied to resolve the detected conflicts. The complex network (CN) aims to avoid collisions by choosing the suitable UAV and the corresponding UAV with a local optimization scope. The CN algorithm only allows for one UAV to perform conflict resolution at one time, which is the principle of DDOA. 3. Core Algorithm Analysis This section describes the complex network model for UAV collision avoidance system and the collision avoidance algorithm in detail. The collision avoidance algorithm consists of two parts: The key node selection algorithm and the collision avoidance direction-selection algorithm. 3.1. Representation of Cooperative UAV Flight and Conflict Detection The background of the problem is that when large scale UAV clusters meet, it is very likely that conflicts will occur. For simplicity, weather perturbations are not discussed when analyzing the trajectories. UAVs share a common airspace in the same or different flight levels, and a UAV cluster can be considered safe only if the separation between each UAV is larger than the given safety distance. A sequence of discrete points is used to model the trajectories of the UAV, and speed information can also be obtained from these discrete tracking points. Because UAV clusters frequently experience threats, and some threats are regarded as potential threats, it is more efficient to perform one resolution at a time rather than comprehensively considering each threat. By constructing a complex network of UAV groups and analyzing the nature of the network under threat conditions to determine the appropriate CA strategy, the human–machine group can self-organize CA under the condition of performance constraints, thus significantly improving the real-time performance and applicability of CA in UAV clusters. The key node in the complex network theory exhibits important properties [33–36], and it exhibits the greatest effect on the entire complex network system. The UAV CA algorithm based on the complex network proposed herein considers the safety indicators of the UAV group system at each moment, such that the overall safety of the entire flight process is maximized. The CA system consists of two key algorithms: A key node selection algorithm and a collision direction-selection algorithm. The algorithm expresses the threat between the UAVs using a complex network. The UAV network model is established based on three parameters: Flight speed, flight angle, and security zone. In a real-time communication environment, all UAVs fly toward a given target and must perform CA operations. Each UAV owns a flight trajectory control module, and the CA resolution calculated by the CA system interacts with the surrounding UAVs through real-time data links. The UAV node is used to indicate that the state of the UAV is monitored by the system in real time. Its attributes are expressed as {UAV, velocity, position, angle, t, state, sense, strength, approaching Appl. Sci. 2019, 9, 3943 5 of 18 time, time cost, key-node}. The UAV indicates the UAV number, velocity indicates the UAV speed, position indicates the current position of the UAV, angle indicates UAV flight angle, t indicates the current time of the UAV, state indicates the UAV flight status, ‘sense’ indicates that the UAV collision direction is selected, ‘strength’ indicates the change direction of the UAV, ‘approaching time’ indicates the remaining time of the UAV approaching the nearest point, ‘time cost’ indicates the entire flight time consumption, and ‘key-node’ indicates whether the UAV was selected as a key node. To detect collisions, UAVs are identified in the Cartesian system. The motion of each UAV is driven by the following:  .   →  →   →   v i (t) cos ϕti cos θti    r x,i (t)  .  r. x,i (t)  → , V (t) = →  →     →  R i (t) =  r y,i (t) R i (t) =  →  =  v (t) cos ϕt sin θt (1)  r. y,i (t)   i   i i i   → → v i (t) sin ϕti     r z,i (t)  →  r z,i (t) → In the above formula, R i (t) and Vi (t) are defined as the position vector and speed vector, respectively, of UAVi at time t. The foot mark i represent i-th UAV. The subscript indicates the horizontal axis system and refers to the height. Let θti designate the horizontal angle which is the direction of speed Vi (t) in the x − y plane (counter-clockwise along the x axis). Let ϕti represent the angle of climb, which is the direction of speed Vi (t) in the vertical plane (measured from the x − y plane, with positive being the upwards direction and negative being the downwards direction). The maximum climb (tilt) angle, which is determined by the performance of the UAV, is defined to limit the maximum angle at which the UAV will rise and fall in the vertical plane when performing a CA mission, thus setting the maximum pitch change angle. The constraint can be expressed as follows: ϕi,min ≤ ϕi (t) ≤ ϕi,max (2) The collection of points in the next time step is as follows: → → → → P(i) = {s ∈ V|k R i − R s k = v · ∆t; ϕmax = arg sin(|zi − zs |/k R i − R s k)} (3) Note that zi and zs are defined as the current attitude of UAVi and the candidate attitude of UAVi in the next time step, respectively. The symbols v and ∆t refer to the absolute velocity of UAVs and time interval of each time step, respectively. When the first pair of UAVs involved in the conflict is detected, each UAV maintains its speed before the first pair of UAVs reaches the danger point to simplify the conflict situation. To determine if a threat of conflict exists, the scope and vertical criteria must be met. This time interval is defined as the time closest to the point of approach (CPA) on the horizontal plane: →t →t r h,i − r h,j t →t →t Th,ij = , ( r h,ij · v h,ij < 0) (4) →t →t v h,i − v h,j · cos(αtij − βtij ) →t →t αtij = arctan( r x,ij / r y,ij ) (5) →t →t βtij = arctan( v x,ij / v y,ij ) (6) The formula defines Th,ijt as the time to CPA in the case of a horizontal plane between UAV and i UAV j at time t. The subscript h indicates the horizontal plane, while i and j indicate the serial numbers Appl. Sci. 2019, 9, 3943 6 of 18 of the UAVs. The equation is defined on the condition that the denominator is not equal to zero. The following equation is defined as the time CPA in the case of a vertical plane: →t r z,ij t →t →t Tz,ij = , ( r z,ij · v z,ij < 0) (7) →t v z,ij A traffic alert (TA) event will be triggered if the following conditions are satisfied: t+T t t t+Th,ij t t h,ij (0 < Th,ij < TTA ) ∧ (0 < Tz,ij < TTA ) ∧ (Dh,ij < DMORRA ) ∧ (Dz,ij < ZTHRRA ) (8) DMORRA and ZTHRRA are the range and altitude limit values for safe aviation, respectively, Appl. Sci. 2019, 9, x FOR PEER REVIEW 7 of 19 in case of the appearance of slow-closure-rate encounters which make the time threshold values not t t+Th,ij t t+Th,ij t t + Th ,ij and D t +Th ,ij are the horizontal and vertical distances between UAV and UAV at CPA, t feasible. DDh,ij feasible. h , ij Dz ,ij and z,ij are the horizontal and vertical distances between UAV i i and jUAV j at respectively. We used the time when the first pair of UAVs triggers a TA alert. From that moment, we CPA, respectively. We used the time when the first pair of UAVs triggers a TA alert. From that built a network where each UAV is represented as a node. If a pair of nodes collides, the node will moment, we built a network where each UAV is represented as a node. If a pair of nodes collides, the establish a connection. At that moment, all nodes in the system will verify if they are approaching; if node will establish a connection. At that moment, all nodes in the system will verify if they are they are approaching each other, the two nodes are connected: approaching; if they are approaching each other, the two nodes are connected: →t  t →t t →  tt  t→t ( rh·,ijv⋅h,ij ( r h,ij vh,ij<<00) ) ∧∧ (( rrz ,z,ij ij ⋅ v· zv < 0) ,ijz,ij < 0) (9) (9) In the case of a pair of UAVs, CA is based on the following principles. As shown in Figure 1, In the case of a pair of UAVs, CA is based on the following principles. As shown in Figure 1, UAV1 UAV1 cruises from the right side to the left side, whereas UAV2 and UAV3 cruise from the left side cruises from the right side to the left side, whereas UAV2 and UAV3 cruise from the left side to the right to theAt side. right thisside. time,AtUAV this time, UAV is selected as the key node. According to the scene at this time, 1 is selected1as the key node. According to the scene at this time, the direction thethe of direction of theand CA direction CA the direction and the target point target at the timepoint atpredicted of the the time of CPAthechoosing predictedanCPA choosing upward climb an or upward climb or a downward climb may be selected. Furthermore, UAV a downward climb may be selected. Furthermore, UAV2 and the UAV3 must 2maintain the minimumand the UAV 3 must maintain safety the minimum safety altitude-limitation. altitude-limitation. The specific The specific key node selection methodskey andnode selection direction methods selection and strategies direction are selection discussed later. strategies are discussed later. Figure 1. Figure 1. Traffic alert and Traffic alert and CA CA (collision (collision avoidance) avoidance) system system (TCAS) (TCAS) coordination coordination for for multi-unmanned multi-unmanned aerial vehicles aerial vehicles (UAVs). (UAVs). 3.2. Key Node 3.2. Key Node Selection Selection The The first first problem problem to to solve solve inin an an encounter encounter is is choosing choosing aa key UAV to key UAV to amend amend its its trajectory trajectory toto avoid avoid collision. collision. The Themethod methodproposed proposed herein decides herein the priority decides of theofUAVs the priority via the via the UAVs analysis of the network the analysis of the formed networkinformed the collision in the space. collisionThe mostThe space. threatened one is chosen most threatened one isaschosen the keyasnode. Obviously, the key the UAV node. Obviously, with the higher priority performs collision avoidance operations the UAV with the higher priority performs collision avoidance operations first. first. To To clarify, clarify,ititisisnecessary necessary to to define define the thecollision collision space. space. When When the the system system detects detects aa pair pair of of UAVs UAVs exhibiting a collision risk, the system obtains all the possible UAV numbers exhibiting a collision risk, the system obtains all the possible UAV numbers to be hit in the future to be hit in the future through through global global data data analysis, analysis, and then it and then it obtains obtains the the status status of of the the corresponding corresponding UAV. UAV. TheThe airspace airspace formed formed byby these these UAVs UAVs is is called called the the collision collision space. space. Among Among them, them, aa connection exists between connection exists between UAVsUAVs that are at the risk of collision, and a network is formed for analysis. An adjacency that are at the risk of collision, and a network is formed for analysis. An adjacency matrix is applied matrix is applied to to record record the the status status between between the the UAVs. UAVs. The The node node represents represents the the UAV, UAV, andand the the edge edge represents represents the the relationship between the UAVs. A set is defined to contain the nodes that satisfy the condition. As shown in Figure 2, the nodes represent the corresponding UAVs. the UAV with the higher priority performs collision avoidance operations first. To clarify, it is necessary to define the collision space. When the system detects a pair of UAVs exhibiting a collision risk, the system obtains all the possible UAV numbers to be hit in the future through global data analysis, and then it obtains the status of the corresponding UAV. The airspace formed Appl. by these Sci. 2019, 9, 3943UAVs is called the collision space. Among them, a connection exists between UAVs 7 of 18 that are at the risk of collision, and a network is formed for analysis. An adjacency matrix is applied to record the status between the UAVs. The node represents the UAV, and the edge represents the relationship between the relationship between the UAVs. UAVs. AAsetsetisisdefined definedtoto contain contain thethe nodes nodes thatthat satisfy satisfy the the condition. condition. As As shown in Figure 2, the nodes represent the corresponding shown in Figure 2, the nodes represent the corresponding UAVs. UAVs. Figure 2. Conceptual graph for key node node selection. selection. After the network is completed, a network description is established to describe the internal characteristics of the multi-UAV system. The nodes in the network exhibit different properties, because the position and velocity are different and the degree of threat is different. Thus, the key node can be chosen according to its properties. First, the easiest method to measure the degree of threat of the UAVs is by examining the situation of the connected edges in the network; this is called the key node selecting method based on degree. To render the network closer to the real scene, it is insufficient to only use the number of edges. The degree approach requires more detailed information, and the edge is a good choice. The edges should contain extra properties to monitor the true relationship of the UAVs. The relationship includes the distance and relative speed between UAVs. The closer a pair of UAVs is to each other, the more dangerous they are. Therefore, the weight of each side depends on the relative distance and relative speed between the UAVs, which is defined as the edge weight between the UAVs: → → υ ij · d ij ωij = → → (10) d ij · d ij → → where υ ij is the relative speed of the vector direction, and d ij represent the distance between the UAVs. This implies that the faster the two UAVs approach each other, the greater is the weight of the edge between the two UAVs. The degree of importance of each node can be based on the weight of all edges. The summed value represents the degree of importance of the node. There are many other methods to rank the importance of nodes [34,37]. Another improved method based on node contraction [38] proposed herein aims to excavate the deeper threat of the UAV networks. There are three reasons for using the node contraction method in selecting the key node: First, the principle of key node selection is to fully exploit the importance of the node, for which the network cohesion is a good measure. The choice of key nodes is not only related to the current degree, it is also to the location of the node. If the node is in the fortress position, shrinking the node will result in a much smaller average path length of the network, resulting in greater network cohesion. This method to obtain key nodes is consistent with the key node selection requirements of the UAV group in performing CA operations, i.e., to determine the current and potential threats of the drone to perform CA. Second, the node contraction method has a shorter calculation time than the node deletion method and the median method, and it does not exponentially increase when the number of nodes is small. Finally, the original method of selecting the key nodes has its limitations, and the changes in the structure of the cluster cannot be detected well. Therefore, this method is not considered for some structurally important points. Appl. Sci. 2019, 9, 3943 8 of 18 Notation Sets: U = {u1 , u2 , u3 , . . . , un } sets of node representing UAV flying in the air sector under control. C = {c1 , c2 , c3 , . . . , cn } ⊆ U × U, sets of edge implying pairs of UAV that are potentially in conflict. n→ → → → o V = v 1 , v 2 , v 3 , . . . , v n ⊆ U × U, sets of relative velocity between UAVs. G = (U, C), graph representing the complex network of the UAV group. Matrix: W = [ωij ] the threat matrix implying the degree of threat between pairs of UAVs. H = [hij ] the adjacency matrix of graph G. D = [dij ] matrix that store the relative distance between UAVs. Variables: l(G): The average distance of graph G. ui : Node that is chosen to be removed. First, the UAV group is modelled, the UAV is replaced by a node, and the threatened UAV pairs are connected. Subsequently, the complex network can be represented by G = (U, C) where G is an undirected connected graph. We assume n nodes and m edges exist, and we represent n UAVs and m threats among them, respectively. The threat coefficient ωij = vij /dij shows the degree of threat between UAVi and UAV j , where vij represents the projection value of the relative speed between UAVi and UAV j on its relative position, and dij represents the relative position between UAVi and UAV j . Definition 1. ωij is defined as follows:  → → → →  υ ij · d ij υ ij · d ij >0     → → , → → d ij · d ij d ij · d ij   ωij =   (11)  → → υ ij · d ij  >0  0,     → →   d ij · d ij The adjacency matrix of graph G is represented by H = [hij ]; it contains n rows and n columns. The element hij is defined as follows: 1, 0 < dij < ∞ ( hij = (12) 0, dij = ∞ The network cohesion is defined below. The degree of network cohesion depends on the connectivity of each node in the network, calculated using the arithmetic mean of the shortest distance between all pairs of nodes. The degree of network cohesion is still related to the density of the network; this is because in the UAV network, the more the number of UAVs, the greater is the threat of the entire network and the greater is the network cohesion. The cohesion of graph G is defined as follows: ωij ωij P P (n − 1) · ω i, j∈U i,j∈U ∂[G] = = P = P (13) n·l dij dij i,j∈U n· i, j∈U n(n−1) Appl. Sci. 2019, 9, 3943 9 of 18 The following equation shows the importance of node ui [38]: ωn ∂[G] n·l(G) n · l(G) · ωn−ki − (n − ki ) · l(G − ui ) · ωn IM(ui ) = 1 − = 1− ωn−k = (14) ∂[G − ui ] i n · l(G) · ωn−ki (n−ki )·l(G−ui ) Because l(G − ui ) < l(G), ki > 1, ωn < ωn−ki , the following inequality 0 < IM(ui ) < 1 can be achieved. The algorithm is as follows: [38]: Input: H, W Output: IM Step 1: Calculate the shortest relation distance matrix D = [dij ] for all pairs of UAV nodes; Step 2: Calculate the initial cohesion of the UAV network; Step 3: Assess the node importance of all UAV nodes. Step 4: Calculate the network cohesion ∂[G − ui ] after the node ui shrinks by calculating the shortest distance matrix Di and the threat matrix Wi between the node pairs after contraction, thereby obtaining IM(ui ). As shown in the algorithm steps above, the time complexity of the whole algorithm depends on the calculation of the shortest distance matrix D between all pairs of nodes after the node shrinks, because the calculation of the threat matrix Wi only requires the deletion of the nodes around the current node ui . Because the time complexity of the Floyd algorithm is O(n3 ), the complexity of the algorithm will reach O(n4 ), owing to the need to shrink the operation of each node. The algorithm for calculating Di given below is based on the initial matrix; the calculation can be completed with a small number of operations, and the time complexity of the entire algorithm will reduce to O(n3 ). The direct distance matrix S = [sij ]n×n can be obtained by the adjacency matrix H = [hij ]. At this point, the change law of the shortest distance d0pq relative to dpq between any pair of nodes (up , uq ) ∈ U × U in the UAV network is as follows: (1) If up , ui and uq , ui : O 1 If dpi + diq = dpq , it is implied that node ui is on the shortest path between node up and node uq ; then, d0pq = dpq − 1; O 2 If dpi + diq ≥ dpq + 1, it is implied that node ui is not on the shortest path between node up and node uq ; then d0pq = dpq . (2) If only one node ui exists between node up and node uq , then d0pq = dpq − 1; (3) If up = uq = ui , then d0pq = dpq = 0. 3.3. Direction of Collision Avoidance Most of the CA systems use the upward or downward strategy to avoid collision. The sense selection can reduce the domino effect. For instance, two new encounters will occur if the key UAV descends, and three will occur if it climbs. The implementation of airborne CA systems II (ACAS II) [39] indicates that it is effective and safe to perform resolution in the vertical direction; therefore, the method of this paper only considered the vertical change in direction. As shown in Figure 3, the key node direction choosing method was as follows. To generate a suitable direction for CA, we expect the UAV cluster system to minimize the cost and to choose the fastest method to reach the target, and we expect the domino effect to be as small as possible in the CA process. The objective function for direction choosing should include three elements: (1) The Number of Intruder UAV Designated Directions n 1, zi > zk X ( Zk (m) = sgn(xk (i)), xk (i) = (15) −1, zi ≤ zk i=i Appl. Sci. 2019, 9,Sci. Appl. 39432019, 9, x FOR PEER REVIEW 10 of 19 10 of 18 As shown in the algorithm steps above, the time complexity of the whole algorithm depends on (2) the Network of the shortest distance matrix D between all pairs of nodes after the node shrinks, Robustness calculation because Complex the calculation networks rely on of their the threat matrix Wii.e., robustness, onlytherequires abilitythe ofdeletion a networkof the tonodes around maintain itsthe connectivity 3 current node u . Because the time complexity of the Floyd algorithm is O ( n ) when a fraction of its vertices is damaged, for their function and performance [40]. In the context of i , the complexity of the algorithmmodel, this new conflict will reach n 4 ) , owing totothe it isO(important need to shrink improve the operation reliability as much of each node. The algorithm as possible. The robustness of for calculatingis D the edge connection i given below defined is based on the initial matrix; the calculation can be completed with as follows: a small number of operations, and the time complexity of the entire algorithm will reduce to O(n3 ) . nby P The direct distance matrix S = [ sij ]n×n can be obtainedP n adjacency matrix H = [h ] . the ωij ij At this point, the change law  Mm distance id='pq1 j=relative Imof the shortest  1 υij to d pq between any pair of Rk (m) = R = , M= , ωij = (16) nodes (u p , uq ) ∈ U × U in the UAV Nmnetwork Nmis−asIfollows: 2 dij (1) If u p ≠ ui and uq ≠ ui : (3) Connected Component of the Network ○ 1 If d pi + diq = d pq , it is implied that node ui is on the shortest path between node u p and node u q ; then, d 'pq = d pq − 1 ; Im Cm Ck (m) = C( )= (17) ○ 2 N If d pi + diq ≥ d pq +1 , it is implied that node N m ui is not − Ithe m on m shortest path between node u p ' Nm stand forand node uin N nodes q ; then d pq = d pq . at time m, Im is the number of nodes that are rejected in the the network network at(2)time m, one If only ui exists between Cm represents andnode node u p and the number node u q , then of sub-maps d 'pq =network. of the d pq − 1 ; In summary, the objective 'function at time m consists of three elements and is represented as (3) If u p = uq = ui , then d pq = d pq = 0 . follows: Ok (m) = ω1 Zk (m) + ω2 Rk (m) + ω3 Ck (m) 3.3. Direction of Collision Avoidance (18) Most of the CA systems use the upward or downward strategy to avoid collision. The sense In this formula, selection ω1 , ω2the can reduce , and ω3 are domino the effect. Forcorresponding instance, two newweight coefficients encounters will occur ifofthe the three key UAV indicators. Different weights satisfy descends, and threethe willrequirements occur if it climbs.of different The scenarios implementation of CA. of airborne CAFor instance, systems II (ACASif ωII)2 < ω3 , it is [39] the implied that indicates that it is effective connected component and safe is to perform more resolutionthan important in the the vertical direction; therefore, robustness the of the network. It is method of this paper only considered the vertical change in direction. As shown in Figure 3, the key noteworthy that a normalization process is required in this computation. node direction choosing method was as follows. Figure 3. Conceptual graph for sense selection. Figure 3. Conceptual graph for sense selection. To generate a suitable direction for CA, we expect the UAV cluster system to minimize the cost Assuming and to that choosethe thedirection is amended fastest method to target, to reach the avoidand thewe collision indomino expect the effect[to the period m,bemas+small N − 1] and the corresponding position as possible in theand CA pitch process.angle sequences The objective of UAV function are, for direction k respectively, choosing should include three elements: T = [pk (UAV (1) The Number ofpkIntruder m), pDesignated . . . , pk (m + N − 1)] k (m + 1),Directions ψk = [ψk (m), ψk (m + 1), . . . , ψk (m + N − 1)]T the following constraints of the optimization process should be satisfied: " # " # 1 B q(ϕk (m + ∆t)) = (ϕk (m + ∆t) − ϕk (m + ∆t)) − ≤0 −1 −A  xk (m + ∆t) − [xk (m) + vk cos(ϕk (m)) cos θk · ∆t]     y (m + ∆t) − [ yk (m) + vk cos(ϕk (m)) sin θk · ∆t]   ek (pk (m + ∆t)) =  k   = 0 (19)   zk (m + ∆t) − [zk (m) + vk sin ϕk · ∆t]   ϕk (m + ∆t) − [ϕk (m) + ωk (m) · ∆t]  Appl. Sci. 2019, 9, 3943 11 of 18 In the formula above, q(ϕk (m + ∆t)) implies that the direction change should be in the range of (A,B). ek (pk (m + ∆t)) indicates that the pitch angle at time m + ∆t is determined by the global information at time m and the constrains of the CA system. By repeating the optimization procedure, the complete amending trajectory is generated. The goal of this new model is to decompose the network as quickly as possible. This implies that after the critical node is cleared, the robustness of the connection and the number of connected components should be as low as possible. This implies that the UAVs are leaving each other and the risk of crashes is becoming smaller. Meanwhile, this model also considers the number of intruder UAVs of the corresponding direction, implying that a lesser domino effect will occur because the density of the intruder UAV is small due to the sense selected by the model. In the simulation step, the new model obtains the key nodes and evaluates the robustness of the connections or the number of connected components when different choices are made by the critical nodes. The key nodes will select the course of change and form a new network in the next simulation step. The new model will again select the key node. The model will monitor multiple UAV systems to determine if a TA alert is present; if the model stops searching for critical nodes, multiple UAV systems will be considered safe. It is noteworthy that this method is a solution to the deadlock problem. The reasons for solving the deadlock problem are as follows: During the entire simulation period, the threat detection between real-time UAVs is performed and the network is built in real time; the key nodes are selected, and the direction of CA is selected. Therefore, the deadlock problem can be avoided. The algorithm proposed in this paper has a simple logical structure and strong stability. Decision calculations consume very little time and meet real-time requirements. The algorithm can be applied to engineering practice due to the following reasons. Firstly, the method we proposed is designed from the perspective of scheduling. Our innovation was to abstract the collision avoidance problem of the UAV to the network’s collapse problem, aiming to solve the collision avoidance problem by using network theory. Second, the proposed method allows the UAV to perform only one maneuver—that of the flight level change. This maneuver is used in the traffic collision avoidance system (TCAS), and it is simple but very effective. As such, it can definitely be applied to the UAV collision avoidance system. The algorithm is shown below as Algorithm 1. Algorithm 1 Novel Model Based on Key Node Research Input: Initial position, target position, initial velocity of the UAVs, the number of UAVs N, the simulation step length S Output: Trajectories of all the UAVs R for j = 1:S do for i = 1:N do Calculate whether a TA is issued; while TA is issued do Construct the network of the UAVs involved; Calculate the relative distance and relative velocities of all the involved UAVs; Select the key node; Calculate the weight (relevant to the number of intruders, robustness and connected component); Choose the optimal directions for the key node; Update the current direction of the key node; if no collision then break; else R[i][j] = the coordinate of the step j of UAV i; Calculate the state of all the UAVs involved in the scenario; if all the UAVs have reach the safe area then break; return R; Appl. Sci. 2019, 9, 3943 12 of 18 Appl. Sci. 2019, 9, x FOR PEER REVIEW 13 of 19 4. Simulation and Results distance was set to 35 m. The simulation step size was set to 1 s. The pitch angle ranged from −0.785 We conducted two different simulation experiments on nine UAV scenarios for two different key to 0.785. It was assumed that each UAV can sense the global information and, as such, simulation node selection methods. The cruise speed set by the UAVs was 10 m/s. The minimum separation calculation time can be ignored. The emulation code was run on a T450 laptop(Lenovo Thinkpad distance was set to 35 m. The simulation step size was set to 1 s. The pitch angle ranged from −0.785 T450 Laptop, Lenovo, Beijing, China) with an Intel i7 processor of 2.6 GHz and 8 GB of RAM. to 0.785. It was assumed that each UAV can sense the global information and, as such, simulation calculation time can be ignored. The emulation code was run on a T450 laptop(Lenovo Thinkpad T450 4.1.A Case Scenario Laptop, Lenovo, Beijing, China) with an Intel i7 processor of 2.6 GHz and 8 GB of RAM. Figure 4a shows the UAV trajectory under the first key node selecting method, which is barely 4.1. A Case based ScenarioIn this scenario, the first time a TA appears is when UAV1 and UAV9 collide with on degree. eachFigure other. 4a After the the shows analysis of the CAunder UAV trajectory system, thewe discovered first key node that whenmethod, selecting avoiding the imminent which is barely collision, some potential collisions may occur. Therefore, when the TA issues a warning based on degree. In this scenario, the first time a TA appears is when UAV1 and UAV9 collide at 18:15:22 with for the each first After other. time, the the system of the UAV analysisselects as the key CA1 system, we node and instructs discovered that whenUAVavoiding 1 to descend to avoid a the imminent conflict. As collision, some thepotential simulation step progresses, collisions UAV1 andwhen may occur. Therefore, UAV9theburst into each TA issues other’s a warning collision at 18:15:22 for the firstThis volumes. time, the system time, selects UAV UAV1 becomes as the the 1key node,keyand node and instructs a climb UAV instruction to descend is1 sent to UAV1to toavoid avoida conflict. As the simulation step progresses, collision. The result in Figure 4b demonstratesUAV 1 and UAV burst the conflict 9 into each other’s collision resolution maneuver and that the volumes. This time, UAV becomes the key node, and a climb instruction is UAV is sent to difference in this method arises because the key node in the same scenario 1 1 to avoid different. Whencollision. selecting The result in Figure 4b demonstrates the conflict resolution maneuver and that the key node by the node contracting method, the system selects the UAV9 as the key node the difference in this and method arises because the key node in the same scenario is different. When selecting the key node instructs UAV9 to descend to avoid conflict. These UAVs form a decentralized CA situation. This can by the node contracting method, the system selects the UAV9 as the key node and instructs UAV9 to be explained by the fact that the robustness of the network in this area has been minimized. The descend to avoid conflict. These UAVs form a decentralized CA situation. This can be explained by the intuitive image is that three UAVs have spread out and no cross routes have appeared, thus avoiding fact that the robustness of the network in this area has been minimized. The intuitive image is that the domino effect. three UAVs have spread out and no cross routes have appeared, thus avoiding the domino effect. (a) (b) Figure4. Figure 4. Trajectory Trajectory(a) (a) with with key key node node selecting selectingbased basedon on degree degreeand and(b) (b) with with key key node node selecting selectingbased based onthe on thenode nodecontraction contractionmethod. method. Table Table22summarizes summarizesthethe trajectory of UAV trajectory of and UAV UAV11 and UAV99 (the (thewaypoints waypointsare arerecorded recordedevery every44ss for for simplicity) simplicity) in the same scenario under two different key node selecting methods. The UAV starts to in the same scenario under two different key node selecting methods. The UAV starts to change changeits itsangle angleatat18:15:23 18:15:23and andreturns returnstotothe theoriginal originaltrajectory trajectoryatat18:17:23. 18:17:23. From this simulation, the results indicate that the CA system can facilitate the UAV cluster to effectively avoid collision in bothTable scenarios. The figure 2. Waypoints illustrates of partial the relative distance of each pair of trajectory. UAVs that are involved in a collision scenario involving nine UAVs and four conflicts among them. It Scenario 1 Scenario 2 is noteworthy that the distance between UAVi and UAV j is defined by dij . Evidently, the separation of Amended Trajectory Amended Trajectory each pair of UAVTimeis never ID intruded upon, according to Figure ID 5. For instance, in Scenario 1, UAV1 X(m) Y(m) Z(m) X(m) Y(m) Z(m) encounters UAV18:15:19 UAV9 , and 6 and UAV1 the distances 876.27 6123.72 between 695.67 UAV9them are always 5756.46 1762.27 larger than the separation 688.52 distance even at the closest 18:15:23 UAV1 point904.54such6095.45 that the695.53 safety UAV9 flight is1722.34 guaranteed, and 688.48 5754.09 it is the same in Scenario 2. 18:15:27 UAV1 919.23 6080.77 692.11 UAV9 1695.73 5752.51 685.49 18:15:31 UAV1 933.91 6066.08 688.69 UAV9 1669.11 5750.93 682.51 18:15:35 UAV1 948.60 6051.40 685.27 UAV9 1642.50 5749.35 679.53 18:15:39 UAV1 963.28 6036.71 681.85 UAV9 1615.89 5747.76 676.55 18:15:43 UAV1 977.97 6022.03 678.43 UAV9 1589.28 5746.18 673.57 UAV1 992.65 6007.34 675.02 UAV9 1562.66 5744.60 670.58 18:15:47 UAV9 1482.76 5739.86 688.48 UAV1 1074.25 5925.74 695.53 18:15:51 UAV1 1007.33 5992.66 671.60 UAV9 1536.05 5743.02 667.60 Appl. Sci. 2019, 9, 3943 13 of 18 Table 2. Waypoints of partial trajectory. Scenario 1 Scenario 2 Amended Trajectory Amended Trajectory Time ID ID X(m) Y(m) Z(m) X(m) Y(m) Z(m) 18:15:19 UAV1 876.27 6123.72 695.67 UAV9 1762.27 5756.46 688.52 18:15:23 UAV1 904.54 6095.45 695.53 UAV9 1722.34 5754.09 688.48 18:15:27 UAV1 919.23 6080.77 692.11 UAV9 1695.73 5752.51 685.49 18:15:31 UAV1 933.91 6066.08 688.69 UAV9 1669.11 5750.93 682.51 18:15:35 UAV1 948.60 6051.40 685.27 UAV9 1642.50 5749.35 679.53 18:15:39 UAV1 963.28 6036.71 681.85 UAV9 1615.89 5747.76 676.55 18:15:43 UAV1 977.97 6022.03 678.43 UAV9 1589.28 5746.18 673.57 UAV1 992.65 6007.34 675.02 UAV9 1562.66 5744.60 670.58 18:15:47 UAV9 1482.76 5739.86 688.48 UAV1 1074.25 5925.74 695.53 UAV1 1007.33 5992.66 671.60 UAV9 1536.05 5743.02 667.60 18:15:51 UAV9 1442.83 5737.48 688.48 UAV1 1102.53 5897.46 695.53 UAV1 1022.02 5977.97 668.18 UAV9 1509.44 5741.44 664.62 18:15:55 UAV9 1402.90 5735.11 688.48 UAV1 1130.81 5869.18 695.53 UAV1 1036.70 5963.29 664.76 UAV9 1482.83 5739.86 661.64 18:15:59 UAV9 1362.97 5732.74 688.48 UAV1 1159.10 5840.89 695.53 UAV1 1051.39 5948.60 661.34 UAV9 1456.21 5738.28 658.66 18:16:03 UAV9 1323.04 5730.37 688.48 UAV1 1187.38 5812.61 695.53 UAV1 1066.07 5933.92 657.92 UAV9 1429.60 5736.70 655.67 18:16:07 UAV9 1283.11 5727.99 688.48 UAV1 1215.67 5784.32 695.53 UAV1 1080.76 5919.23 654.50 UAV9 1402.99 5735.12 652.69 18:16:11 UAV9 1243.18 5725.62 688.46 UAV1 1243.95 5756.04 695.53 UAV1 1105.63 5894.36 653.54 UAV9 1366.39 5732.94 651.91 18:16:15 UAV9 1203.25 5723.25 688.48 UAV1 1272.24 5727.75 695.53 18:16:19 UAV1 1133.90 5866.09 653.43 UAV9 1326.46 5730.57 651.86 18:16:23 UAV1 1162.16 5837.83 653.27 UAV9 1286.54 5728.20 651.82 18:16:27 UAV1 1190.43 5809.56 653.13 UAV9 1246.61 5725.83 651.77 18:16:31 UAV1 1218.70 5781.29 652.99 UAV9 1206.68 5723.45 651.72 18:16:35 UAV1 1233.38 5766.61 656.41 UAV9 1180.07 5721.87 654.70 18:16:39 UAV1 1248.07 5751.92 659.82 UAV9 1153.46 5720.29 657.68 18:16:43 UAV1 1262.75 5737.24 663.24 UAV9 1126.84 5718.71 660.67 18:16:47 UAV1 1277.44 5722.55 666.66 UAV9 1100.23 5717.13 663.65 18:16:51 UAV1 1292.12 5707.87 670.08 UAV9 1073.62 5715.55 666.63 18:16:55 UAV1 1306.81 5693.18 673.50 UAV9 1047.01 5713.97 669.61 18:16:59 UAV1 1321.49 5678.50 676.92 UAV9 1020.39 5712.39 672.59 18:17:03 UAV1 1336.18 5663.81 680.34 UAV9 993.78 5710.80 675.58 18:17:07 UAV1 1350.86 5649.13 683.75 UAV9 967.17 5709.22 678.56 18:17:11 UAV1 1365.55 5634.44 687.17 UAV9 940.56 5707.64 681.54 18:17:15 UAV1 1380.23 5619.76 690.59 UAV9 913.94 5706.06 684.52 18:17:19 UAV1 1394.92 5605.07 694.01 UAV9 887.33 5704.48 687.50 18:17:23 UAV1 1416.39 5583.60 695.65 UAV9 850.73 5702.31 688.21 It is noteworthy that the distance between UAVi and UAV j is defined by dij . Evidently, the separation of each pair of UAV is never intruded upon, according to Figure 5. For instance, in Scenario 1, UAV1 encounters UAV6 and UAV9 , and the distances between them are always larger than the separation distance even at the closest point such that the safety flight is guaranteed, and it Appl. Sci. 2019, 9, 3943 14 of 18 is the same in Scenario 2. (Scenario 1) (Scenario 2) Figure Figure 5. 5. Minimal Minimal distance distance of of UAVs UAVs for two scenarios. Because the Because the algorithm algorithmimplements implementsaasingle-machine single-machinechangechangeinin direction, direction,when thethe when UAVUAVdoesdoes not reach Appl. the Sci. pre-set 2019, 9, x CA FOR altitude, PEER REVIEW the relative distance of the UAV will be reduced slightly; not reach the pre-set CA altitude, the relative distance of the UAV will be reduced slightly; however, however, it 15 will of 19 not it collide will because not collide it is atitthe because predicted is at moment the predicted of collision. moment The UAV of collision. has already The UAV flown has already to a safe flown to a safe altitude; altitude; therefore, therefore, although although a closea trend close may trendoccur may during occur during CA, theCA, UAV thecan UAV can maintain maintain a safe a safe relative relative distance distance in the subsequent in the subsequent flight path. flight path. Figure 6 depicts the situation Figure 6 depicts the situation during the during the final final CA CA ofof the the UAV UAV group. group. The The connection connection represents represents the potential threat between the UAVs. The connection method is the same the potential threat between the UAVs. The connection method is the same as the logic of building as the logic of building the the network. It is related to the relative distance between nodes and the relative network. It is related to the relative distance between nodes and the relative speed between the nodes. speed between the nodes. The network node can identify key nodes and subsequently select The network node can identify key nodes and subsequently select the most appropriate CA direction the most appropriate CA direction according according to the network to theproperties. network properties. (Scenario 1) (Scenario 2) Figure 6. Distribution of UAVs UAVs during during aa threat threat in in two two scenarios. scenarios. We Wedefined defineda anewnew amount amount to measure to measurethe effect of UAVs the effect when selecting of UAVs key nodes when selecting keyusing nodesdifferent using methods. We defined different methods. Wethis amount defined this as the residual amount threat factor. as the residual threat The meaning factor. of this of The meaning quantity is the this quantity measure of the sum is the measure of theofsum threats to the UAVs of threats to the in a certain UAVs area centered in a certain on the critical area centered on thenode after critical a conflict node after a avoidance maneuver. That is, it evaluates the safety of the UAV in a certain area conflict avoidance maneuver. That is, it evaluates the safety of the UAV in a certain area centered centered on the point on after performing CA maneuver. Therefore, the threat coefficients of the UAVs the point after performing CA maneuver. Therefore, the threat coefficients of the UAVs in the area in the area are added, and different are added, key and node selections different key nodeare compared. selections The residual are compared. threat The coefficient residual threat under the algorithm coefficient under the compares the advantages and disadvantages of the key node selection algorithm compares the advantages and disadvantages of the key node selection method. method. The residual threat The factor (RTF) can be defined as follows: residual threat factor (RTF) can be defined as follows: N N tk P i,j ij =1 tijk RTFt = i , j =1 (20) (20) RTFt = N To depict the overall threat during the CA procedure, N an average residual threat factor was defined To depict the overall threat during the CA procedure, an average residual threat factor was as follows: trt defined as follows: X RTF = RTFt , TAt