International Journal of Automation and Computing  2018, Vol. 15 Issue (4): 402-416 PDF
Unmanned Aerial Vehicle Formation Inspired by Bird Flocking and Foraging Behavior
Tian-Jie Zhang
Beijing University of Aeronautics and Astronautics, Beijing 100191, China
Abstract: This paper considers a multiple unmanned aerial vehicles (UAV) formation problem and proposes a new method inspired by bird flocking and foraging behavior. A bidirectional communication network, a navigator based on bird foraging behavior, a controller based on bird interaction and a movement switch are developed for multi-UAV formation. Lyapunov′s second method and mechanical energy method are adopted for stability analysis. Parameters of the controller are optimized by Levy-flight based pigeon inspired optimization (Levy-PIO). Patrol missions along a square and an S shaped trajectory are designed to test this formation method. Simulations prove that the bird flocking and foraging strategy can accomplish the mission and obtain satisfying performance.
Key words: Unmanned aerial vehicle (UAV)     formation     bird flocking     foraging     pigeon-inspired optimization     Levy-flight
1 Introduction

Unmanned aerial vehicles (UAVs) have been developed at a fast speed in recent years. There appears a popular tendency, simple UAVs swarm for sophisticated missions instead of a single complex UAV. This situation partly attributes to the limit of UAV technology, still cooperation always works better than solo. Natural phenomena have proved the success of animal cooperation such as mammals herd, fish school, insects swarm and birds flock. Hence, it is a wise choice to adopt a bio-inspired method for missions with multiple UAVs.

As for a specific task, UAVs should operate in ordered formation. There are four main methods for UAV formation at present, leader-wingman, artificial potential field (APF), virtual structure and behavior-based method. The leader-wingman strategy acts as a developed method for multilevel UAV systems. This method was adopted without communication between the leader and its wingman, instead, vision sensors were employed to collect necessary information[1]. Although the ignorance of communication can simplify the whole system, it cannot deal with the problem of occluded target, especially in a large UAV flock. APF was adopted to achieve desirable UAV formation[2, 3]. APF is easily adapted to vehicles with different dynamics but parameter optimization is always a difficult problem. Virtual structure also works effectively for multi-UAV formation[4], however when turning, UAVs at different positions fly at different speed and trajectories are asymmetric and redundant. Behavior-based method solves the problem of multi-UAV formation, a UAV possesses a limited number of actions to choose, every team member makes a decision and coordinates together[5], the main problem is that the stability and convergence of the system cannot be guaranteed.

The idea of bio-inspired multiple UAV formation originates from behavior of small gregarious birds, such as pigeons and starlings, which are quite familiar to human. The underlying rules of bird flocking and foraging are quite fascinating. When gathering and forming a flock with numerous members, how do pigeons make a common decision and act as one unit? When travelling in close formation, how do starlings with sight partly blocked at the back detect their front space and choose the direction? When flying over a long distance, what is the communication method to avoid flock separation? When attacked by predators, what is the reaction strategy to make a sharp turn?

Much time and effort has been devoted into exploring bird behavior. Masure and Allee[6] studied social order in pigeon flocks and drew a conclusion that sex and body strength influenced dominance. Nagy et al.[7] studied hierarchical group dynamics in pigeon flocks and proposed a multi-level leader-follower strategy. Chen et al.[8] analyzed the same global positioning system (GPS) data as Nagy′s but simplified the flock model and extracted a two-level leader-follower structure. Pettit et al.[9] proposed a different flock mechanism where speed determined leadership, slow pigeons followed and obeyed the fast ones. Cavagna et al.[10] recorded experiment data and put forward scale-free behavioral correlations among starlings, it was considered that a single starling could receive and transfer information with remote partners over the whole flock. Lukeman et al.[11] built a pigeon model with strong short-range repulsion, intermediate-range alignment, and long-range attraction within circular zones as well as weak but significant front-sector interaction to describe flocking behavior. Under this assumption, birds communicated with each other and acquired useful surrounding information to avoid collision and keep stable formation. Cavagna et al.[12] also investigated ordered structure of starling flocks and found boundary information inflow enhanced correlation in flock. Ballerini et al.[13] observed European starlings and adopted topology distance instead of metric distance to model their interaction. It is found that a bird only communicated with six or seven fixed neighbors. Topology distance could effectively condense the flock when disturbed by a predator. Pettit et al.[14] tracked homing pigeons and showed that interaction between pigeons stabilized a side-by-side configuration, promoting bidirectional information exchange and reducing the risk of separation. Pomeroy and Heppner[15] photographed turning movement of rock doves and found doves did not maintain fixed relative positions during a turn. Yomosa et al.[16] quantitatively investigated the internal motion of a flock based on pairwise statistics and showed pigeons adopted a mixture of two idealized and fundamentally different turning types, namely, parallel-path and equal-radius type turning. The latter type was more flexible and often used during predatory raid. All these results have provided background knowledge and experiment data for this work.

The outline of this paper is as follows. Section 1 describes background knowledge and gives a brief introduction of innovation. Section 2 gives description of formation patrol missions and a UAV model in details. It is demonstrated in Section 3 how a controller and a navigator for a single UAV are devised based on bird behavior. In Section 4, the stability of the bio-inspired formation plan is analyzed utilizing the stability theory of ordinary differential equations with limited simplification of the model. Next comes Section 5, Levy-flight based pigeon inspired optimization (Levy-PIO) is adopted for controller parameter optimization. In Section 6, numerical simulations are implemented and results of multi-UAV formation inspired by bird flocking and foraging behavior are shown, following analysis reflects both advantages and disadvantages of the method. The final part is the conclusion, which summarizes the whole content and points out our future work.

2 Description of the patrol mission and the UAV model 2.1 Formation patrol mission

This paper considers two patrol missions along different predesigned trajectories with n UAVs. One UAV acts as the leader which is able to acquire the coordinates of a series of waypoints, while others are followers and only know their aim points in the relative moving coordinate system fixed to the leader. There is a fixed communication topology simulating scale-free interaction[10], while separation, alignment, attraction and front interaction force come into effect based on metric distance[13]. Both topology and metric distances are combined to design a UAV′s controller for multi-UAV formation[17]. Topology distance communication makes formation more robust while metric distance communication works well in near distance. Both distances provide a balance for performance and cost.

The first desired patrol trajectory surrounds a rectangle area, the trajectory simulates the outline of an essential building. For the sake of simplicity, the size of the rectangle is given in body length of a UAV, rather than metric length. In most biology papers, interaction range is measured by body length because different birds have different sizes. This tradition is inherited in our work, the distance based interaction, such as separation, alignment, attraction and front interaction, works according to the body length of a UAV. To avoid the inconvenience of unit conversion, the size of the rectangle is 300×300 in body length.

The second patrol trajectory is an S shaped route, it simulates complex terrain where a UAV formation should fly along a curve to avoid obstacles. The trajectory is in symmetry to its midpoint and the circumcircle radius of the half trajectory is 200 body length.

 Download: larger image Fig. 1. Communication topology and formation shape

Fig. 1 provides an example of communication topology, n UAVs are supposed to form the shape of an equilateral triangle, one apex angle shows the direction of the formation movement. The leader UAV flies at the front and stays at the vertex of the moving triangle. Follower UAVs stay at the symmetry positions and along the triangle edges.

In Fig. 1, dashed lines show UAVs are located on edges of an equilateral triangle, while solid lines are adopted to measure L in (8) and have no concern with the predesigned formation structure. Black arrows demonstrate there is a directional communication topology link between the two neighbors.

The patrol trajectories are shown in Fig. 2, respectively.

 Download: larger image Fig. 2. Patrol trajectory of the mission

In Fig. 2(a), n UAVs patrol counterclockwise along a square trajectory shown by arrows. Their flight direction is collinear with the aircraft longitudinal axis pointing from the tail to the nose. In Fig. 2(b), n UAVs fly along the S shape trajectory from (0, 0) to (0, 800) directed by arrows.

These missions order UAV formation to fly along designed trajectories while covering a large area. It is assumed that each UAV carries a camera to photograph the ground and detect abnormal issues. Under the assumption that the UAVs execute the task in a safe environment without any obstacle or armed resistance on the way, the formation would better fly horizontally, which means UAVs should keep level flight. In the patrolling process, n UAVs should maintain the configuration of the equilateral triangle as long as possible, except turning. The two missions appear quite similar to pigeons′ normal circular flight, so it is reasonable to refer to bird behavior.

2.2 UAV model and performance constraint

Because the model originates from small gregarious birds, the UAVs described by the model are supposed to be small, light, flexible with limited load. These features mean every UAV is inexpensive and performance defects can be easily compensated by swarm. For the patrol mission, a single UAV can only monitor a small area, while n UAVs can consist of a sensor network and cover a larger place. An orderly and effective formation is then necessary.

Because the time constant of aircraft attitude loop always has a smaller magnitude than its position loop, a UAV′s attitude should be simplified when discussing a formation patrol problem, so a mass point model is suitable to describe the movement of a UAV. In what follows, a UAV flies in the inertial coordinate frame (x, y, h), the mass center of UAV is fixed to its airframe, the gravitational acceleration g is a constant, moment acting on a UAV is ignored. Under these assumptions, the movement can be described as[18, 19]

 $\begin{split}& {v_x} = v\cos \gamma \cos \chi \\& {v_y} = v\cos \gamma \sin \chi \\& \dot h = v\sin \gamma \end{split}$ (1)

where vx is the speed along x-axis and vy is y-axis speed. v is the airspeed, h is the altitude, γ is the flight-path angle and $\chi$ is the heading angle. UAV dynamics is given in (2):

 $\begin{split}& \dot v= F\\& \dot \gamma = \frac{g}{v}(n\cos \phi - \cos \gamma )\\& \dot \chi = \frac{{gn\sin \phi }}{{v\cos \gamma }}\end{split}$ (2)

where F is the sum of forces acting on a UAV and is described in (3). $\phi$ is the bank angle and n is the normal overload. The control variables are F, $\phi$ and n. The model outputs are v, γ and χ. n and $\phi$ are adopted to reduce the absolute value of γ to 0 so as to maintain level flight, F drives each UAV to its destination.

Flexibility is a key factor for movement of a UAV. Considering a pigeon or starling flock in the wild, each member must avoid collision in close formation with its sight blocked in most directions and its body situated in complicated flow fields, still for the sake of acquiring the movement information of neighbors and adjusting motion state in time, each UAV should possess a large range of velocity and a long communication distance. The velocity falls in the range of [0, 5] in body length per second and the communication range is [0, 17] in body length. The design process will be fully illustrated in the controller part in Section 3.

3 Navigator and controller design based on bird behavior 3.1 Communication topology with individual difference

Communication ability is crucial for multi-UAV formation, only with adequate redundancy, large range, channel availability, high signal-noise ratio and power, can UAVs exchange information in time and react to the varying environment[20]. Small UAVs have sacrificed their load capacity and can not accomplish perfect environment perception by themselves. There appears a popular tendency of combination of contributed processing and distributed perception to solve the problem[21], meanwhile, distributed control has attracted much attention, especially in robot swarm field. Distributed control abandons the traditional concept of “central processing unit” (not referring to central processing unit (CPU) in a computer, but a central point in communication and processing topology), and proposes that all the members obey a specific protocol. Although each UAV decides what to do next only with several neighbors, the whole swarm reaches consensus and acts as one unit.

It is a wise idea to combine the concepts above, a new structure with distributed perception, global or local contributed processing and distributed control must possess a broad application prospect in the future. The following design has been learnt from this wonderful idea.

N UAVs are numbered from 1 to n counterclockwise, UAV1 stays at the vertex of the triangle at the front of the formation and serves as the leader with the superior authority, it can directly receive the coordinate of next waypoint and lead the whole team towards the destination, this information can be programmed before or transmitted by a ground control station. UAV1 also receives position and velocity information from UAV2 and UAVn by the communication topology. The job of leader UAV simulates foraging behavior of a wild bird flock, it is supposed that some experienced or visually keen individuals usually spot food at first and fly towards the feeding area, others just follow them. The flock is naturally divided into leaders and followers[9]. Meanwhile in the formation, the left n–1 UAVs become followers without knowing their final destination, they can only fly around the leader and acquire the relative coordinates in the moving coordinate system fixed to UAV1, the leader. This strategy reflects individual difference among team members.

Traditional bird behavior studies utilize the same model to describe every flock member for the sake of simplicity, however it is unreasonable to treat everyone in the same way for the fact that flocking birds are distinguished by size, age, strength and social status[22]. Different models should be employed for different UAVs because a future UAV swarm will be heterogeneous and composed of different kinds of UAVs. For example, in a patrol mission, a UAV team is possibly consisted of guides, photographers, relays and a commander. As for a future battle, an armed UAV squadron should be composed of fighters, bombers, scouters, early-warning UAVs, etc. That is the reason why individual difference should be considered in the UAV formation patrol mission.

UAV2 and UAVn are the second level leaders, both of them receive position and velocity information from UAV1, the top leader, and from their own followers. UAV2 acquires information from UAV3, while UAVn from UAVn–1. The second level leaders implicate the feature of contributed processing in the way they gather necessary information from their own leader and followers. Contributing all their positions and velocities together certainly helps decision because abundant information can reflect system status clearly. Communication relationship between UAV2 and UAVn is bidirectional, both of them work as transmitters and receivers. Under the desirable condition, UAV2 serves as a relay which passes the leader information from UAV1 to UAV3, also monitors the condition of UAV3. The relationship between UAVn and UAVn–1 is quite the same.

All n UAVs collect information of their own positions and velocities, and then exchange it with their neighbors. After considering conditions of the neighbors, each UAV can take next step. The information is first collected distributedly, then processed concentratedly and applied distributedly finally. This communication topology implicates the advanced control strategy of distributed perception, contributed processing and distributed control. Although the quite simple topology can hardly be a typical example of the advanced control strategy, we still want to bring this idea into our fundamental research and that is why the topology is devised as it is.

3.2 Navigator based on bird foraging

There are many similarities between bird foraging and UAV patrol. Birds, especially starlings and pigeons, are always busy looking for food, the flocks can search their surroundings over and over again during day time. Once a smart one finds a food resource, it would fly towards there directly and send a signal to its mates[6, 10]. Other birds do not have to look for the resource by themselves, instead they only need to follow the smart bird[9].

Both birds and UAVs can explore the local environment and collect the information they are interested in, once the destination located, they would fly there as fast as possible. On the way, a large number of starlings crowd in the sky without any block or collision[7], a UAV swarm should learn from this behavior.

Double-integrator dynamics is employed to describe UAVs′ behavior. The higher the integrator order is, the smoother the UAV trajectory is, so navigator based on bird foraging can directly generate a desirable trajectory with a smooth curve. Force acting on a UAV has a common form in (3) with its mass omitted.

 $\frac{{{\rm{d}}\vec v}}{{{\rm{d}}t}} = {\vec f_{individual}} + {\vec f_{social}}.$ (3)

As for the UAV1, it is driven by a navigation force ${\vec f_{individual,\,1}}$ produced from its velocity ${\vec v_1}$ and the distance ${\vec a_1}$ between its present position and destination. ${\vec f_{social,\,1}}$ is an interaction force demonstrated in (11). For the sake of simplicity, noise is not introduced into the whole system. ${\vec f_{individual,\,1}}$ is consisted of two parts[11]:

 ${\vec f_{individual,1}} = {\vec a_1} - {\gamma _1}{\vec v_1}$ (4)
 ${\vec a_1} = {\vec x_{des,1}} - {\vec x_1}\quad\quad\quad\;\;$ (5)

where the distance ${\vec a_1}$ works as intrinsic driving force towards the waypoints on the patrol trajectory and ${\gamma _1}{\vec v_1}$ works as drag. γ1 is a drag regulator influencing the speed of UAV1, the larger γ1 is, the slower UAV1 moves and it takes much more time to accomplish the patrol mission. However, the team leader cannot fly too fast because its followers may not keep pace with it and the formation breaks. The drag ${\gamma _1}{\vec v_1}$ is a rate feedback factor which can increase damping ratio and stabilize the double integrator system. The details would be illustrated in Section 4.

Follower UAVs employ the same form of force as UAV1 in (4) and (5). There is difference in destination ${\vec x_{des,flo}}$ and the drag regulator γflo. ${\vec x_{des,flo}}$ is the relative destination of followers UAV2 – UAV n in the moving coordinate system fixed to UAV1 and γflo is smaller than γ1 because $\left\| {{{\vec x}_{des,flo}} - {{\vec x}_{flo}}} \right\| \ll \left\| {{{\vec x}_{des,1}} - {{\vec x}_1}} \right\|$ . To avoid UAV1 flying too fast, the drag force should increase, another regulator kslow is introduced and (5) is transformed into the following equation:

 $\left[ {\begin{array}{*{20}{c}}{{{\vec a}_1}}\\{{{\vec a}_{flo}}}\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{k_{slow}}({{\vec x}_{des,1}} - {{\vec x}_1})}\\{\displaystyle\frac{1}{{{k_{slow}}}}({{\vec x}_{des,flo}} - {{\vec x}_{flo}})}\end{array}} \right]\quad\quad\;\;$ (6)
 ${k_{slow}} = {{\rm{e}}^{ - di{s^2}}}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$ (7)
 $dis = \left\| {{{\vec x}_1} - \left[ {\frac{{({{\vec x}_2} + {{\vec x}_n})}}{2} + \frac{{{{\vec v}_2} + {{\vec v}_n}}}{{\left\| {{{\vec v}_2} + {{\vec v}_n}} \right\|}} \times L} \right]} \right\|$ (8)

where L shown in Fig. 1 is the distance between UAV1 and the midpoint of UAV2 and UAVn in the standard formation. dis reflects the distance difference between the actual and desirable positions of UAV1, $\left\| {\; \cdot \;} \right\|$ refers to Euclidean norm. If UAV1, UAV2 and UAVn are located on the vertexes of an equilateral triangle, dis = 0, kslow = 1, (7) and (8) do not work. Otherwise dis > 0, k slow < 1, (6) slows down UAV 1, but speeds up the left followers.

3.3 Controller based on neighborhood interaction

First of all, attitude problem should be solved. When patrolling along a rectangle side in Fig. 2, UAVs are supposed to maintain level flight, control variables are φ = 0 and n = cosγ – 10γ, (2) can be substituted by

 $\begin{split}& \dot \gamma = - \frac{{10g}}{v}\gamma \\& \dot \chi = 0\end{split}$ (9)

where the gravitational acceleration g is a constant, v is limited to [0, 5] body length per second. The flight path angle γ will decay to 0 quickly and the heading angle $\chi=\chi$ (0).

When turning, coordinated turn is adopted

 $\dot \chi = \frac{g}{v}\tan \phi$ (10)

then n = cosγ – 10γ → γ = acos(ncosφ), the control variable is $\phi = \arctan \left( {\displaystyle\frac{{v{{\dot \chi }_{idl}}}}{g}} \right)$ and $n = \displaystyle\frac{1}{{\cos \phi }}$ , ${\dot \chi _{idl}}$ is the ideal yaw rate.

The navigator has been designed to lead the UAV swarm, the communication relationship based on topology distance assures cohesion. Even if the formation configuration is broken by exterior disturbance and an individual leaves its partners far away, the communication link is not broken and can guide the individual back to the team.

However, flight trajectories sometimes have oscillations if only with navigation force because a navigator merely provides desirable aim points and position relationship, varying motion states are not fully considered. Positions are directly given as a command, in order to arrive at the desirable positions, velocities and accelerations are then calculated, last are driving forces. This procedure seems like a pattern of passive execution. It implicates a control style from positions to forces.

On the contrary, if there is another control mode from forces to positions, formation result would become better. Neighborhood interaction in starling or pigeon flocks complies with the pattern of active perception, a bird can collect position information of its neighbors autonomously. If two individuals stay far from each other, they tend to close[11]. If they stay too near, both would like to separate. When keeping a suitable interval distance, they may maintain the same velocity to reserve the current state[23].

The controller is based on metric distance without fixed neighbors[17, 23]. Once an object UAV enters the range of interaction, the subject can measure the distance and a corresponding force is adopted. The interaction zone of a UAV is a series of concentric circles added with a front sector of θ ± 30°, θ is the angle between the line of action and the longitudinal axis of the UAV, as is shown in Fig. 3. Following bird studies[6, 24], the social interaction is modelled in terms of separation ${\vec f_{sep}}$ , alignment ${\vec f_{ali}}$ , attraction ${\vec f_{att}}$ and frontal interaction ${\vec f_{front}}$ . To determine the relative contribution of four different kinds of forces, normalized forces weighted by coefficients are employed. The sum force is[11]

 ${\vec f_{social}} = {\omega _{sep}}{\vec f_{sep}} + {\omega _{ali}}{\vec f_{ali}} + {\omega _{att}}{\vec f_{att}} + {\omega _{front}}{\vec f_{front}}$ (11)

where ωsep, ωali, ωatt and ωfront are the weighting coefficients to be optimized in Section 5.

As to separation, UAVi is led by force ${\vec f_{sep,i}}$ to fly in the direction of the vector pointing to itself from an object UAVj which enters its perception zone in (12). Inside the separation zone, ${\vec f_{sep,i}}$ increases with the distance following a halved Gaussian, ${{\rm{e}}^{\frac{{ - {{\left\| {{{\vec x}_i} - {{\vec x}_j}} \right\|}^2}}}{{{r_{sep}}^2}}}}$ , $\left\| {{{\vec x}_i} - {{\vec x}_j}} \right\| \ge 0$ and rsep is the standard deviation of the Gaussian set. When the distance between UAVi and UAVj $\left\| {{{\vec x}_i} - {{\vec x}_j}} \right\| = 0$ , $\left\| {{{\vec f}_{sep,i}}} \right\| = 1$ reaches its maximum, when $\left\| {{{\vec x}_i} - {{\vec x}_j}} \right\| = {r_{sep}}$ reaches the maximum, $\left\| {{{\vec f}_{sep,i}}} \right\| = {{\rm{e}}^{ - 1}}$ is the smallest value. $\frac{{{{\vec x}_i} - \vec x{}_j}}{{\left\| {{{\vec x}_i} - {{\vec x}_j}} \right\|}}$ controls the direction of the separation force

 $\begin{split}& {{\vec f}_{sep,i}} = {{\rm{e}}^{\frac{{ - {{\left\| {{{\vec x}_i} - {{\vec x}_j}} \right\|}^2}}}{{{r_{sep}}^2}}}} \times \frac{{{{\vec x}_i} - {{\vec x}_j}}}{{\left\| {{{\vec x}_i} - {{\vec x}_j}} \right\|}}\\& 0 < \left\| {{{\vec x}_i} - {{\vec x}_j}} \right\| \le {r_{sep}}\\&\qquad \!\left| \theta \right| \ge {30^\circ }.\end{split}$ (12)

As to alignment, UAVi is subject to a force ${\vec f_{ali,i}}$ to align with any UAVj that enters its notched ring zone. The radius of the inner ring is rsep, while the outer radius is rali. Alignment doesn′t employ the form of Gaussian, instead it is

 $\begin{split}& {{\vec f}_{ali,i}} = \frac{{{{\vec v}_{norm,j}} - {{\vec v}_{norm,i}}}}{{\left\| {{{\vec v}_{norm,j}} - {{\vec v}_{norm,i}}} \right\|}}\\& {r_{sep}} \le \left\| {{{\vec x}_i} - {{\vec x}_j}} \right\| < {r_{ali}}\\&\qquad\!\left| \theta \right| \ge {30^\circ }\end{split}$ (13)

where ${\vec v_{norm,i}}$ is the normalized velocity of ${\vec v_i}$ . ${\vec f_{ali,i}}$ aims at changing the direction of UAVi once the two UAVs do not coordinate their velocities. If ${\vec v_{norm,i}}$ deviate to left compared to ${\vec v_{norm,j}}$ , ${\vec f_{ali,i}}$ points to the right side so as to drive UAVi to turn right.

As regards its cohesion behavior, UAVi is attracted by the force pointing to UAVj located in its metric neighborhood. To avoid breaking the formation configuration, UAVs are set to cohere more strongly at the border of the swarm than its interior by adopting a Gaussian force like the separation part

 $\begin{split}& {{\vec f}_{att,i}} = {{\rm{e}}^{\frac{{ - {{\left( {{r_{att}} - \left\| {{{\vec x}_j} - {{\vec x}_i}} \right\|} \right)}^2}}}{{{{\left( {{r_{att}} - {r_{ali}}} \right)}^2}}}}} \times \frac{{{{\vec x}_j} - {{\vec x}_i}}}{{\left\| {{{\vec x}_j} - {{\vec x}_i}} \right\|}}\\& {r_{ali}} \le \left\| {{{\vec x}_i} - {{\vec x}_j}} \right\| \le {r_{att}}\\&\left| \theta \right| \ge {30^ \circ }\end{split}$ (14)

where ratt is the attraction radius, $\displaystyle\frac{{{{\vec x}_j} - {{\vec x}_i}}}{{\left\| {{{\vec x}_j} - {{\vec x}_i}} \right\|}}$ represents the direction of attraction which pulls UAVi towards its neighbor UAVj. The Gaussian ${{\rm e}^{\frac{{ - {{\left( {{r_{att}} - \left\| {{{\vec x}_j} - {{\vec x}_i}} \right\|} \right)}^2}}}{{{{\left( {{r_{att}} - {r_{ali}}} \right)}^2}}}}}$ translated by $\left\| {{{\vec x}_j} - {{\vec x}_i}} \right\|$ is monotonously increasing when $\left\| {{{\vec x}_j} - {{\vec x}_i}} \right\| \in [{r_{ali}},{r_{att}}]$ . When $\left\| {{{\vec x}_j} - {{\vec x}_i}} \right\| = {r_{ali}}$ , $\left\| {{{\vec f}_{att,i}}} \right\|$ reaches the minimum of e–1. When $\left\| {{{\vec x}_j} - {{\vec x}_i}} \right\| = {r_{att}}$ , $\left\| {{{\vec f}_{att,i}}} \right\|$ reaches the maximum of 1 on the border of the interaction zone.

Finally frontal interaction zone is a sector area symmetric about the sight line of a bird, as to a UAV, the sight line is replaced by the longitudinal axis. In this sector, only separation and attraction work because alignment mainly brings improvement to lateral individuals. Like forces in (12) and (14), the frontal interaction also adopts a Gaussian force

 {\vec f_{front,i}} = \left\{ \begin{aligned}& {{\rm{e}}^{\frac{{ - {{\left( {{r_{att}} - \left\| {{{\vec x}_j} - {{\vec x}_i}} \right\|} \right)}^2}}}{{{{\left( {{r_{att}} - {r_{mid}}} \right)}^2}}}}} \times \frac{{{{\vec x}_j} - {{\vec x}_i}}}{{\left\| {{{\vec x}_j} - {{\vec x}_i}} \right\|}}\\& \;{r_{mid}} \le \left\| {{{\vec x}_j} - {{\vec x}_i}} \right\| < {r_{att}}\\&\qquad \left| \theta \right| < {30^ \circ }\\& {{\rm{e}}^{\frac{{ - {{\left\| {{{\vec x}_i} - {{\vec x}_j}} \right\|}^2}}}{{{r_{mid}}^2}}}} \times \frac{{{{\vec x}_i} - {{\vec x}_j}}}{{\left\| {{{\vec x}_i} - {{\vec x}_j}} \right\|}}\\& \left\| {{{\vec x}_i} - {{\vec x}_j}} \right\| < {r_{mid}}\\&\qquad \left| \theta \right| < {30^ \circ }\end{aligned} \right. (15)

where ${r_{mid}} = \displaystyle\frac{{({r_{sep}} + {r_{att}})}}{2}$ is the border of separation and attraction. ${\vec f_{front,i}}$ is monotonously decreasing in [0, rmid) but increasing in [rmid, ratt]. ${\vec f_{front,i}}$ is smallest at $\left\| {{{\vec x}_j} - {{\vec x}_i}} \right\| = {r_{mid}}$ , $\;\;\min({\vec f_{front,i}}) = {{\rm e}^{ - 1}}$ and $\;\;\max({\vec f_{front,i}}) = 1$ at $\left\| {{{\vec x}_j} - {{\vec x}_i}} \right\| = 0\;\;{\rm or}\;\; {r_{att}}$ . Concave function of ${\vec f_{front,i}}$ assures stability of the frontal interaction and $\left\| {{{\vec x}_j} - {{\vec x}_i}} \right\| = {r_{mid}}$ is the stable point. Once UAVi enters the frontal separation zone of UAVj, that is $\left\| {{{\vec x}_i} - {{\vec x}_j}} \right\| < {r_{mid}}$ , it is repelled out of the district. Otherwise, UAVi would be pulled back when it flies in the frontal attraction zone with ${r_{mid}} \le \left\| {{{\vec x}_j} - {{\vec x}_i}} \right\| < {r_{att}}.$

 Download: larger image Fig. 3. Different interaction zones

 Download: larger image Fig. 4. Interaction zone in 3D perspective

Different neighborhood interaction zone is shown in Fig. 3, separation effects in the central notched circle, attraction works in the outer notched ring and frontal interaction is employed in the sector. In Fig. 4, force magnitude distribution is shown, the altitude of peaks and gray levels are utilized to measure magnitude. At the central point and on the peripheral border, altitude reaches 1, the peak and the margin are colored white in Fig. 4. Magnitude decreases in the middle range with low height and dark color. Alignment is not displayed in Fig. 4 because it doesn′t change according to distance. Still it should be demonstrated that the magnitude in the picture is relative value, the true force needs to be modified by the weighting coefficients in (11).

3.4 Flight pattern switch based on flocking

During patrol along the predesigned trajectory, UAVs are supposed to accomplish a series of tasks including flocking, formation, level flight, turning and adjustment. Many tasks consist of the entire mission just like many organs consist of the whole human body.

At the beginning of the patrol, n UAVs take off from different sites whose average distance is larger than the metric communication distance. At that moment, only topology communication comes into operation. The leader, UAV1 flies towards the first waypoint, because the navigation force of UAV1 is stronger than other forces, UAV1 would overtake other UAVs and occupy a leader position in the front.

The desirable position of UAV2 is generated according to the present positions and velocities of UAV1, UAV2 and UAV3.

 ${\vec v_{temp,2}} = \frac{{{{\vec v}_1} + {{\vec v}_3}}}{{\left\| {{{\vec v}_1} + {{\vec v}_3}} \right\|}}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\;$ (16)
 $\begin{split}{{\vec x}_{des,2}}= & \frac{1}{3} \times \left({{\vec x}_2} + {{\vec x}_1} + \left[ {\begin{array}{*{20}{c}}{\cos \beta }&{ - \sin \beta }\\{\sin \beta }&{\cos \beta }\end{array}} \right] \times {{\vec v}_{temp,2}} +\right. & \left.{{\vec x}_3} + \left[ {\begin{array}{*{20}{c}}{\cos \;(\pi - \beta )}&{\sin \;(\pi - \beta )}\\{ - \sin \; (\pi - \beta )}&{\cos \;(\pi - \beta )}\end{array}} \right] \times {{\vec v}_{temp,2}}\right)\end{split}$ (17)

where ${\vec v_{temp,2}}$ , a temporal variable, is a normalized vector of average velocities of UAV1 and UAV3, the relative destination of UAV2 is the average of its present position and the predicted position by UAV1 and UAV3. In the desirable formation configuration, UAV2 is located at an angle of β counterclockwise to flight direction of the team, on left back side of UAV1 and at an angle of (π – β) clockwise on the right front side of UAV 3, as shown in Fig. 1. Limited by the topology communication, UAV2 can only acquire information from its neighbors, so the formation flight direction is substituted by ${\vec v_{temp,2}}$ in (17). $\left[ {\begin{array}{*{20}{c}}{\cos \beta }&{ - \sin \beta }\\{\sin \beta }&{\cos \beta }\end{array}} \right]$ aims at rotating the direction.

The desirable position of UAV3 is produced by its present position and a predicted position by UAV2. UAVn and UAV2 have symmetric positions, UAVn–1 and UAV3 are symmetric too. Attracted by navigation forces towards desirable positions, n UAVs can soon flock, coordinate velocities and form a stable formation.

 Download: larger image Fig. 5. Turning problem of the fixed triangle formation

Near the turning point, n UAVs stop the level flight pattern, instead they switch into equal radius turning mode. In order to clarify the shortcoming of the fixed triangle formation turning, an example with five UAVs is adopted in Fig. 5. The problem is that inner loop UAVs would turn around. It′s unreasonable that one UAV turns 90° while another turns 270° in the same turning. The root of the problem is that before turning, UAV3 is located on the left side of UAV1, after several steps, UAV1 has already finished turning with UAV3 located at the back. UAV3 generates a trajectory of an irregular circle. To solve the problem, a new turning pattern should be introduced.

Parallel path turning and equal radius turning are two common turning types observed in pigeon flocks[15, 16]. Parallel path turning urges all pigeons to keep a fixed formation configuration while turning, equal radius turning allows them to select different centers but the same radius and turn freely according to their present positions. The configuration is sacrificed for flexibility. To achieve pleasant performance, the equal radius turning is adopted. When the team approaches a waypoint, the leader UAV1 receives the coordinate of next waypoint and informs its followers of switching into turning. All the members decide a turning radius R and calculate their centers, if to turn left, the direction of velocity is rotated 90° counterclockwise and the center lies on the new direction with the distance of R between the center and its present position. The turning angle is divided into 100 parts, for each time step, a UAV turns one part[15]. During the equal radius turning, topology communication doesn′t work, while metric communication is still on duty to avoid collision. The configuration is completely ignored during this period, however for convenience of reconfiguration as follows, turning speeds are supposed to be negotiated. A simple method is that inner UAVs fly more slowly, while outer ones move faster. A way to judge whether a UAV is an inner or outer one is the third component of cross product, employing the velocity ${\vec v_i}$ of UAVi and relative position $({\vec x_j} - {\vec x_i})$ for cross product, if the third component is negative, UAVi is the inner one referring to UAVj, otherwise the outer one.

After turning, all five UAVs fly towards next waypoint and begin to reform the standard formation, however inner UAVs perhaps turn around because its leader drops back behind it and its desirable position is further behind. Adjustment is designed for solution. All n UAVs maintain their directions but adjust speeds according to leadership. The leader tends to fly faster and the follower tends to move more slowly. When the team leader, UAV1 overtakes both UAV2 and UAVn, adjustment ends and level flight pattern begins to work. Adjustment criterion is the angles $\left\langle {{{\vec v}_1},\;{{\vec x}_2} - {{\vec x}_1}} \right\rangle$ and $\left\langle {{{\vec v}_1},{{\vec x}_n} - {{\vec x}_1}} \right\rangle$ , if both angles are larger than 120°, adjustment reaches its goal.

4 Stability analysis of UAV formation based on bird flocking and foraging 4.1 Stability analysis based on ordinary differential equations

Newton second law is employed to describe formation movement and the math model is based on second order differential equations. Stability problems of the bird foraging strategy can be easily transformed into stability of second order differential equations. Individual difference still needs to be considered in this section, the leader and its followers are analyzed differently. The navigation force ${\vec f_{individual}}$ is mainly considered in this part.

The dynamics equation of UAV1 from (3) to (5) can be summarized as

 $\frac{{{{\rm{d}}^2}{x_1}}}{{{\rm{d}}{t^2}}} = {{\rm e}^{ - {{[{x_1} - (\frac{{{x_2} + {x_n}}}{2} + L)]}^2}}}({x_{des}} - {x_1}) - \gamma \frac{{{\rm{d}}{x_1}}}{{{\rm{d}}t}}.$ (18)

For the sake of simplicity, a more abstract form is adopted to substitute (18):

 $\frac{{{{\rm{d}}^2}{x_1}}}{{{\rm{d}}{t^2}}} + \gamma \frac{{{\rm{d}}{x_1}}}{{{\rm{d}}t}} + \alpha {{\rm e}^{ - {{({x_1} - A)}^2}}}({x_1} - k) = 0$ (19)

where xdes is reducible to a fixed value k and $\displaystyle\frac{{\left( {{x_2} + {x_n}} \right)}}{{2 + L}}$ in exponent can be substituted by another fixed value A, because x1 does not appear in $\displaystyle\frac{{\left( {{x_2} + {x_n}} \right)}}{{2 + L}}$ directly. This simplification eliminates coupling in (18), which is very important for the following proof. Apparently, x1 = k is a stable solution of (19), in order to study stability of zero solution, y = x1–k is utilized to change the variable and (19) is transformed into the differential equations as

 $\frac{{\rm{d}}}{{{\rm{d}}t}}\left[ {\begin{array}{*{20}{c}}y\\v\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}v\\{ - \alpha {{\rm{e}}^{ - {{(y + k - A)}^2}}}y - \gamma v}\end{array}} \right].$ (20)

Lyapunov′s second method for stability is employed[25, 26], this classical theory demonstrates that a system $\displaystyle\frac{{{\rm d}x}}{{{\rm d}{t}}} = f(x)$ has a point of equilibrium at x = 0. Consider a Lyapunov function V (x): R2R such that

V (x) = 0 if and only if x = 0.

V (x) > 0 if and only if x ≠ 0.

$\displaystyle\frac{{{\rm d}V(x)}}{{{\rm d}t}} ={\Large\sum}_{i = 1}^2 {\frac{{\partial V}}{{\partial {x_i}}}{f_i}(x) \le 0}$ if and only if x ≠ 0. This equilibrium is stable in the sense of Lyapunov, if for any ε > 0, there exists δ > 0 such that when $\left\| {x(0) - {x_e}} \right\| < \delta$ , then for any t ≥ 0, $\left\| {x(t) - {x_e}} \right\| < \varepsilon$ (Note that f has an equilibrium at xε, f (xε) = 0).

For (20), Lyapunov function is

 $V(y,v) = \frac{1}{2}{v^2} + \int_0^y {\alpha s{{\rm{e}}^{ - {{(s + k - A)}^2}}}} {\rm{d}}s$ (21)

where (21) satisfies the conditions above, V (y, v) = 0 if and only if y = 0, v = 0, V (y, v) > 0 if and only if y ≠ 0, v ≠ 0, $\displaystyle\frac{{{\rm{d}}V(y,v)}}{{{\rm{d}}t}} \!\!=\!\! \displaystyle\frac{{\partial V}}{{\partial y}}\displaystyle\frac{{{\rm{d}}y}}{{{\rm{d}}t}} \!\!+\!\! \displaystyle\frac{{\partial V}}{{\partial v}}\displaystyle\frac{{{\rm{d}}v}}{{{\rm{d}}t}} \!\!=\! \alpha {{\rm{e}}^{ - {{(y \!+\! k \!-\! A)}^2}}}\!\!yv \!-\! v\alpha {{\rm{e}}^{ - {{(y + k \!-\! A)}^2}}}y-$ $\gamma {v^2} = \gamma {v^2} \le 0$ , then the solution of (19), x1 = k = xdes, is stable and UAV1 would fly to its destination by the navigation force.

As to follower UAVs, more simplification is needed because coupling among followers is much more complex, for example, the dynamics equation of UAV2 is

 $\frac{{{{\rm{d}}^2}{x_2}}}{{{\rm{d}}{t^2}}} = {{\rm{e}}^{{{[{x_1} - (\frac{{{x_2} + {x_n}}}{2} + L)]}^2}}}({x_1} + {k_2} - {x_2}) - \gamma \frac{{{\rm{d}}{x_2}}}{{{\rm{d}}t}}$ (22)
 $\frac{{\rm{d}}}{{{\rm{d}}t}}\left[ {\begin{array}{*{20}{c}}y\\v\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}v\\{ - \alpha x - \gamma v}\end{array}} \right].\quad\quad\quad\quad\quad\quad\quad\quad\;\;$ (23)

Equation (22) is reducible to (23), ${{\rm{e}}^{{{[{x_1} - (\frac{{{x_2} + {x_n}}}{2} + L)]}^2}}}$ is replaced by a fixed number α because in the following numerical simulations, ${{\rm{e}}^{{{[{x_1} - (\frac{{{x_2} + {x_n}}}{2} + L)]}^2}}}$ is always a constant in [0.2, 0.4], the stability of follower UAVs is deeply influenced by UAV1, once UAV1 flies stably, others can produce suitable trajectories. Equation (23) has been translated to meet the demand. Lyapunov function is

 $V(y,v) = \frac{\alpha }{2}{y^2} + \frac{1}{2}{v^2}$ (24)

where V (y, v) = 0 if and only if y = 0, v = 0, V (y, v) > 0 if and only if y ≠ 0, v ≠ 0, $\displaystyle\frac{{{\rm{d}}V(y,v)}}{{{\rm{d}}t}} \!=\! \displaystyle\frac{{\partial V}}{{\partial y}}\frac{{{\rm{d}}y}}{{{\rm{d}}t}} \!+\! \displaystyle\frac{{\partial V}}{{\partial v}}\frac{{{\rm{d}}v}}{{{\rm{d}}t}} \!=\! \alpha vy +$ $v( - \gamma v + \alpha y) = - \gamma {v^2} \le 0$ , then movement of UAV2 is stable in the sense of Lyapunov.

4.2 Stability analysis based on artificial potential field

Neighborhood interaction resembles traditional artificial potential field in the way forces can be classified into separation, attraction and alignment, and they are based on distance[27]. Work is the product of force and the displacement in the direction of the force applied. When a UAV is propelled by the forces in (11), it moves in a conservation field and work is done[28]. For separation of UAVi,

 ${\vec f_{sep,i}} = {f_{sep}}(r)\frac{{\vec r}}{r} = {f_{sep}}(r)\frac{x}{r}\vec i + {f_{sep}}(r)\frac{y}{r}\vec j$ (25)

where r is distance variable, $\vec r = (x\vec i,y\vec j)$ shows the position of a UAV, x, y are distance coordinates, $\vec i$ , $\vec j$ are directional vectors.

 $\begin{split}\nabla \times {{\vec f}_{sep}} = & \frac{{\partial {{\vec f}_{sep}}}}{{\partial x}} - \frac{{\partial {{\vec f}_{sep}}}}{{\partial y}} = y\frac{{\partial f(r)}}{{\partial r}}\frac{{\partial r}}{{\partial x}} - x\frac{{\partial f(r)}}{r}\frac{{\partial r}}{{\partial y}} = \\ & \frac{{\partial f(r)}}{r}\frac{{yx - xy}}{{\sqrt {{x^2} + {y^2}} }} = 0 \end{split}$ (26)

where curl $\nabla \times {\vec f_{sep}} = 0$ proves ${\vec f_{sep,i}}$ is a conservation force, then there is a conservation field corresponding to ${\vec f_{sep,i}}$ . Other forces, ${\vec f_{att,i}}$ , ${\vec f_{front,i}}$ have the same result. Energy intensity in the field can be employed to analyze stability of movement. UAVs with a high energy level tend to change their condition in a short period of time, while another with a low level has a tendency of moving slowly within its vicinity. This is a simple method, however, it is an adoptable way to release potential energy of a UAV for stability.

Considering UAVi, its interaction district is divided into five zones as shown in Fig. 3, neighbors in A, C, D, E zones are influenced by conservation force and have potential energy, while kinetic energy is utilized as substitute in zone B. We also define A as a set including all the neighbors of UAVi in zone A, so do B, C, D and E. The sum of energy of UAVi is

 $\begin{split}{V_i} = & \sum_{j \in A} {\int_{{x_j}}^{{r_{sep}}} {{\rm{e}}^{\frac{{ - {s^2}}}{{{r^2}_{sep}}}}} {\rm{d}}s + \frac{1}{2}\sum\limits_{k \in B} {\Delta v_k^2} + } \\& \sum\limits_{l \in C} {\int_{{r_{ali}}}^{{x_l}} {{{\rm{e}}^{\frac{{ - {{({r_{att}} - s)}^2}}}{{{r^2}_{att}}}}}{\rm{d}}s} } + \\& \sum\limits_{m \in D} {\int_{{r_{mid}}}^{{x_m}} {{{\rm{e}}^{\frac{{ - {{({r_{att}} - s)}^2}}}{{{{({r_{att}} - {r_{mid}})}^2}}}}}{\rm{d}}s} + } \\& \sum\limits_{n \in E} {\int_{{x_n}}^{{r_{mid}}} {{{\rm{e}}^{\frac{{ - {s^2}}}{{r_{mid}^2}}}}{\rm{d}}s} } \end{split}$ (27)

where rsep, rali, ratt and rmid are the interaction ranges in (11)–(15), △vk is the relative difference of velocity of UAVi and its neighbors in zone B. j, k, l, m and n represent neighbors of UAVi. When Vi reaches its minimum, it is considered that UAVi is flying under its most stable condition.

When xj = rsep, xl = rali, xm = xn = rmid, all UAVs are flying on the border of different interactions, ${V_i} =\displaystyle\frac{1}{2}\sum\limits_{k \in B} {\Delta v_k^2}$ , only with kinetic energy left in zone B. Then, alignment would coordinate their velocities, the ideal sum of energy reduces to zero at that moment. Although the ideal result is not achieved in the following simulation, this method still leads a way for stable formation. Borders of different forces are stable rings of potential energy for the specific configuration. Once a UAV enters the separation zone, it is driven out by the repulsive force and once it enters the attraction zone, the UAV is again pulled by attraction and moves back. By choosing optimal coefficients, UAV configuration can achieve its best performance.

5 Parameter optimization with Levy-flight based pigeon inspired optimization

Pigeon inspired optimization (PIO) was proposed in 2014 from pigeons′ behavior of seeking routes back to their loft, it has been widely adopted in parameter optimization and image processing[29, 30]. Parameters for optimization form the location of a pigeon. In this paper, Levy-flight strategy is adopted and two new operators are employed for improvement[31]. The operators in Levy-PIO are Levy-flight based search operator and landmark operator with variable step length.

Levy-flight, one of random walk models, simulates a predation process. When a predator observes clues to its prey, it tends to search the area thoroughly with short step, on the other hand, when the predator finds no food there, it will move to another area with long step. This strategy can combine diversity with convergence speed and help a lot to avoid premature problem which is a main disadvantage of basic PIO.

 $\begin{split}& s = \frac{\mu }{{{{\left| v \right|}^{\frac{1}{\delta }}}}},\;\delta = 1.5\\& \mu \sim {{N}}\left( {0,\sigma _\mu ^2} \right),\;v \sim {{N}}\left( {0,\sigma _v^2} \right)\end{split}$ (28)
 ${\sigma _\mu } = {\left\{ {\left. {\frac{{\Gamma \left( {1 + \delta } \right)\sin \left( {\displaystyle\frac{{\pi \delta }}{2}} \right)}}{{\Gamma \left[ {\displaystyle\frac{{\left( {1 + \delta } \right)}}{2}} \right]\delta \times {2^{\frac{{\left( {\delta - 1} \right)}}{2}}}}}} \right\}} \right.^{\frac{1}{\delta }}},\;{\sigma _v} = 1$ (29)
 $\begin{split}& step = s \circ \left( {{X_i}\left( t \right) - {X_{leader}}} \right)\\& {X_{itemp}} = {X_i}\left( t \right) + step \circ randn\end{split}$ (30)

where Xi (t) is the pigeon i′s present location, Xleader is the global best location of the flock leader selected by the cost function in (33). Xitemp is the output of Levy-flight based PIO. s and step are vectors with the same dimension as Xi (t) → s and step are vectors and they have the same dimension as Xi (t). N (0, σ2) represents the normal distribution. $\circ$ denotes the Hadamard product, and randn is a vector with random numbers as its elements, obeying the normal distribution.

 \quad\left\{ \begin{aligned}& {X_i}\left( {t + 1} \right) = {X_{itemp}},\;\;{\rm if}\;{f_{cost}}\left( {{X_{itemp}}} \right) < {f_{cost}}\left( {{X_i}\left( t \right)} \right)\\& {X_i}\left( {t + 1} \right) = {X_i}\left( t \right),\;\;\;{\rm if}\;{f_{cost}}\left( {{X_{itemp}}} \right) \ge {f_{cost}}\left( {{X_i}\left( t \right)} \right).\end{aligned} \right. (31)

The result Xitemp still needs to be compared, if it is a better location, the pigeon should fly there, otherwise the pigeon maintains its present location.

When the pigeon flock flies near their loft, some experienced ones would recognize the familiar landmarks and guide other members towards home. The landmark operator can converge all the pigeon to their leader and the variable step length can almost solve the premature problem.

 \begin{aligned}& {X_i}(t + 1) = {X_i}\left( t \right) + Length \times randn \circ \left( {{X_{leader}} - {X_i}\left( t \right)} \right)\\& \quad Length ={\rm logsig}\left( {\frac{{Nc \times \zeta - t}}{k}} \right)\end{aligned} (32)

where ζ and k are adaptive parameters of logsig function which influence the convergence rate. Nc is the maximum iteration number. Length determines how much the leader can attract the followers, in the beginning of this part, Length = 1, all the followers tend to fly towards the location of their leader, in the end, Length = 0, the followers would maintain their own positions.

The cost function is

 ${f_{cost}} = \sum {Va{r_{vx}}} + \sum {Va{r_{vy}} + \sum {Len} + \sum {Di{s_{near}}} }$ (33)

where $\sum {Len}$ is the total length of trajectories of the five UAVs, Disnear is the sum of distances from the different waypoints to the nearest trajectory point for a single UAV. $\sum {Di{s_{near}}}$ includes the five UAVs. $\sum {Va{r_{vx}}}$ and $\sum {Va{r_{vy}}}$ are the sum of variance of x-velocity and y-velocity, respectively. fcost is designed to acquire the shortest trajectory while reaching the waypoints. The variance is adopted to suppress oscillations and smooth trajectories.

Levy-flight based PIO is utilized to optimize the weighting coefficients for the controller in (11). A pigeon′s location is X = [ωsep, ωali, ωatt, ωfront]T, ωsep ∈(0, 15), ωali ∈(0, 5), ωatt ∈(0, 4), ωfront ∈(0, 2). These intervals are referenced to Lukeman′s result[11]. Fifteen pigeons are initialized at random and the maximum iteration number is 15. Finally, a desirable set of coefficients is acquired and shown in Table 1. The fitness curve is shown in Fig. 6 which fully explains the necessity of the optimization.

Table 1
Interaction weighting coefficients for controller

 Download: larger image Fig. 6. Fitness curve of Levy-flight based PIO

According to Fig. 6, final fitness is fcost = 9 164.25, total trajectory length of the n UAVs reduces 23.4% after optimization. From Table 1, separation takes a main part in interaction because ωsep is nearly 5 times larger than ωali and 2.5 times larger than ωatt. ωfront is a very small amount. This result is consistent with Ryan′s work in the way separation occupies high proportion, alignment and attraction own the same magnitude and front interaction serves as a weak force[11]. The outcome reflects in a stable flock, every pigeon possesses a private space and its neighbors do not come in, which may ensure flight safety.

6 Numerical simulation and analysis

Numerical simulations of multi-UAV formation are implemented in this section to persuade the effect of the controller and the navigator designed according to bird flocking and foraging. After the simulations, results are adopted for analysis. Both advantages and disadvantages are clearly shown and the reason behind them should be paid attention to.

6.1 Patrol mission along square

The patrol trajectory is a 300×300 square as is shown in Fig. 2(a). Five UAVs start the patrol mission from scattered positions. Their initial positions and velocities are listed in Table 2. Coefficients in Section 3 are listed in Table 3 below and β = 120° in (17). Fourth order Runge-Kutta method is employed for numerical simulation and time step is h = 0.01, 80 000 steps in total.

Table 2
Initialization of UAV formation patrol mission

Table 3
Interaction distance coefficients for controller

The formation result is shown in Fig. 7.

 Download: larger image Fig. 7. Formation flight trajectories of five UAVs

The five UAVs are located on edges of a nonstandard triangle with different velocity directions at the beginning. Guided by their navigators, these UAVs start to fly towards the first waypoint at (150, 0). On the way, interaction among UAVs changes the formation shape which is stretched along x-axis while compressed along y-axis, this process can be observed in the zone $\left\{ {} \right.$ (x, y)|x ∈ [50, 100], y ∈ [–10, 10] $\left. {} \right\}$ in Fig. 7.

When the leader UAV1 detects that it is close enough to the first waypoint, it broadcasts the coordinate of the second waypoint, (150, 300). This message is transmitted by the hierarchical communication topology and all team members decide to enter equal radius turning mode and choose turning centers and the radius. The five UAVs begin to draw arcs in Fig. 7, when UAV1 has made a quarter turn, it informs its followers of changing flight mode to adjustment, every UAV maintains its velocity direction but differs in speed. As is shown in $\left\{ {} \right.$ (x, y)|x ∈ [100, 200], y ∈ [0, 50] $\left. {} \right\}$ in Fig. 7, UAV1 flies at the front of the formation in the beginning, but moves to the right border of the team because of the equal radius turning, trajectories cross during this period without any collision. Configuration is broken so the five UAVs can fly freely, which satisfies the former demand. After adjustment, the five UAVs adopt level flight pattern again but with residual oscillations in trajectories, because the method how to generate velocity is changed as the movement pattern is changed. If there is no adjustment, the UAVs would turn around. In this way, adjustment would come into effect. However, problems of switchover are not fully solved by adjustment and this needs improvement in our future work. The five UAVs fly towards the second waypoint, then the third one, (–150, 300), the fourth one, (150, 0) and their final destination, (0, 0). A single trajectory is similar to a closed curve, which means the corresponding UAV flies back to its initial position and accomplishes the patrol mission.

 Download: larger image Fig. 8. Velocity variation of five UAVs in patrol

In Fig. 8, variation of velocity is clearly displayed. The tendency of all five UAVs is quite the same in most time. When t∈[0, 120)∪[160, 315)∪[350, 505)∪[530, 690)∪[710, 800) s, v reduces slowly. In this period of time, all five UAVs fly in level flight pattern. They fly straight forward in a level plain. v oscillates violently because UAVs are turning freely which requires fast speed and variable direction. The process above cycles in whole simulation, main difference lies in the equal radius turning pattern, different UAVs own different velocities due to their positions in the formation. Amplitude is limited in [–5, –5] body length per second considering ability of UAVs.

 Download: larger image Fig. 9. Heading angle ${\chi}$ of five UAVs

In Fig. 9, heading angles of five UAVs are clearly shown. The heading angle of UAV1 is representative, its initial value approaches 0, that is consistent with the heading angle in Table 2 and means UAV1 flies towards the right side. It increases to $\displaystyle\frac{{\rm{\pi }}}{2}$ , π, $\displaystyle\frac{{3{\rm{\pi }}}}{2}$ and π, and UAV1 flies along the whole rectangular trajectory and goes back to its initial district. The other four UAVs are quite the same.

 Download: larger image Fig. 10. Flight path angle γ of UAV1

In Fig. 10, the flight path angle γ has a weak magnitude of 10–15 and UAVs can be reasonably considered to fly in a plane without climbing or diving, that is apparently what is needed in Section 2. The attitude controllers (9) and (10) work. Only the curve of UAV1 is drawn because all 5 curves are highly overlapped.

Fig. 11 shows variation of average distance D among UAVs over time, the average distance is defined as

 $D= \frac{{\left( {\left\| {{{\vec x}_1} - {{\vec x}_2}} \right\| + \left\| {{{\vec x}_2} - {{\vec x}_3}} \right\| + \left\| {{{\vec x}_1} - {{\vec x}_5}} \right\| + \left\| {{{\vec x}_5} - {{\vec x}_4}} \right\|} \right)}}{4}$ (34)
 Download: larger image Fig. 11. Average distance among UAVs

where ${\vec x_i}$ is the position of UAVi, $\left\| {{{\vec x}_i} - {{\vec x}_j}} \right\|$ is the distance between UAVi and UAVj. Under ideal conditions, D = 7.5. From Fig. 11, average distance is much larger than 7.5 in the beginning and reduces to the standard value as time goes by, just before turning, D ≈ 7.5, because equal radius turning does not order UAVs to maintain their relative positions, the formation breaks and expands violently. After entering level flight pattern, the formation compresses again.

Fig. 12 shows energy change of five UAVs, the energy includes potential and kinetic energy according to (27). Analysis in Section 4 has proposed that a stable configuration tends to own a low energy level, it can be proved in Fig. 12 because most curves decrease monotonously piecewise. Energy is released for the sake of configuration improvement and the formation becomes more stable in the meantime.

 Download: larger image Fig. 12. Energy variation of five UAVs

When the UAVs are turning, kinetic energy suddenly increases because of acceleration, which can be confirmed in Fig. 8. The only problem happens to UAV2, its energy increases violently. In order to investigate the reason, we replicate the simulation and calculate the different energy separately and find that the UAV is dragged into the attraction zone of its neighbor and potential energy of attraction field begins to work. This process can also be proved in Fig. 11 where average distance suddenly increases.

6.2 Patrol mission along S shaped route

The patrol trajectory is shown in Fig. 2(b). Three UAVs start the patrol mission with initial positions and velocities are listed in Table 4. Coefficients in Section 3 are not changed and listed in Table 3 and β = 120° in (17). Fourth order Runge-Kutta method is employed and time step is h = 0.01, 150 000 steps in total.

Table 4
Initialization of UAV formation patrol mission

The formation result is shown in Fig. 13.

 Download: larger image Fig. 13. Formation flight trajectories of three UAVs

In Fig. 13, three UAVs form a triangular formation and patrol along an S shaped trajectory from (0, 0) to (0, 800). UAV1 works as the leader and UAV2, UAV3 are followers on the back sides. During the flight, the navigator designed in Section 3.2 guides the UAVs to each waypoint and the controller manages to maintain the standard formation. The trajectories are smooth and accurate, every waypoint in Fig. 2(b) is reached.

 Download: larger image Fig. 14. Velocity variation of three UAVs in patrol

Variation of velocity is clearly displayed in Fig. 14 for the three UAVs. The tendency is much similar to the counterpart in Fig. 8. When a UAV reaches one waypoint, it immediately receives the next waypoint and has to change the direction. This process looks like a starling flock accelerates and turns when encountering predation or obstacles.

Other variables are not shown by figures because of length limit and the results own the same tendency as the counterpart in square mission. Initial heading angle χ is 0, increases to 180° when the UAV reaches (0, 400), the midpoint, $\chi$ decreases to 0 in the end. Average flight path angle $\gamma$ is –2.20×10–15°, which means level flight. Total average distance is 8.70 body length and energy curves have the similar impulse trend as in Fig. 12.

6.3 Shortcomings and analysis

Validity of the bird flocking and foraging strategy has been confirmed in the numerical simulations above, trajectories, velocities, neighbor distance and energy have satisfied the predesigned demand, however, there is still another technical specification, the angle between neighbors. The desirable formation shape is an equilateral triangle, which means $\left\langle {{{\vec v}_i},{{\vec x}_j} - {{\vec x}_i}} \right\rangle \; \in$ $\left\{ {} \right.$ 30°, 50° $\left. {} \right\}$ , UAVi and UAVj become neighbors when they have communication topology links in Fig. 1, n = 5. However it is noticed that $\left\langle {\vec v_2},\;{\vec x_1} - {\vec x_2} \!\right\rangle$ ≈ 23°, $\left\langle \!{\vec v_5},{\vec x_1} \!-\! {\vec x_5}\! \right\rangle$ ≈ 23°, $\left\langle {\vec v_3},{\vec x_2} - {\vec x_3} \right\rangle$ ≈ 30°, $\left\langle {\vec v_4},\;{\vec x_5} - {\vec x_4} \right\rangle$ ≈ 30°. The actual formation seems a concave triangle with a very sharp apex angle. After investigating the whole strategy, attraction in the frontal zone should account for this problem.

 Download: larger image Fig. 15. Resolution of frontal attraction force

If UAVj is located in the frontal attraction zone of UAVi, it is attracted towards UAVi by ${\vec f_{att}}$ . ${\vec f_{att}}$ can be resolved as shown in Fig. 15, ${\vec f_{att,y}}$ drags UAVj to its center zone, it is ${\vec f_{att,y}}$ that hinders UAVj from keeping 150° with UAVi. In addition, there is no compensation to improve the configuration. In order to maintain the desirable distance, the formation is stretched out of its standard shape.

7 Conclusions

Bird flocking and foraging has enlightened people in many fields including multi-UAV formation. In this paper, patrol missions along a square and an S shaped trajectory are designed to test the performance of the formation strategy. A bidirectional communication network is employed for information transmission among UAVs. A navigator based on bird foraging is adopted for multi-UAV navigation, individual difference is fully considered and the leader UAV owns superior rights to its followers. The leader can directly receive coordinates of waypoints, while followers merely acquire relative positions in the moving coordinate system fixed to the leader. Learning from interaction in bird flocks, a controller simulates repulsion, alignment and attraction among birds in the meantime the navigator works. The controller assures flight safety and maintains desirable configuration. In order to fly along a predesigned trajectory, different flight patterns are utilized including flocking, level flight, equal radius turning and adjustment. Stability problems of the bird strategy can be easily transformed into stability of second order differential equations. Lyapunov′s second method and mechanical energy method are adopted to confirm stability. Numerical simulations later prove the UAV formation method from bird flocking and foraging can accomplish the formation patrol mission. Trajectories, velocities, neighbor distance and energy have satisfied the predesigned demand. However, there still exist shortcomings, the trajectories have slight oscillations and the angle is smaller than the desirable value. The reasons are found and these problems will be better solved in our future work.

References
 [1] M. A. Dehghani, M. B. Menhaj. Communication free leader-follower formation control of unmanned aircraft systems. Robotics and Autonomous Systems, vol.80, pp.69-75, 2016. DOI:10.1016/j.robot.2016.03.008 [2] T. Paul, T. R. Krogstad, J. T. Gravdahl. Modelling of UAV formation flight using 3D potential field. Simulation Modelling Practice and Theory, vol.16, no.9, pp.1453-1462, 2008. DOI:10.1016/j.simpat.2008.08.005 [3] Y. Fu, X. K. Wang, H. Y. Zhu, Y. G. Yu. Parameters optimization of multi-UAV formation control method based on artificial physics. In Proceedings of the 35th Chinese Control Conference, IEEE, Chengdu, China, pp. 2614–2619, 2016. [4] W. Ren, N. Sorensen. Distributed coordination architecture for multi-robot formation control. Robotics and Autonomous Systems, vol.56, no.4, pp.324-333, 2008. DOI:10.1016/j.robot.2007.08.005 [5] T. Balch, R. C. Arkin. Behavior-based formation control for multirobot teams. IEEE Transactions on Robotics and Automation, vol.14, no.6, pp.926-939, 1998. DOI:10.1109/70.736776 [6] R. H. Masure, W. C. Allee. The social order in flocks of the common chicken and the pigeon. The Auk, vol.51, no.3, pp.306-327, 1934. DOI:10.2307/4077659 [7] M. Nagy, Z. Akos, D. Biro, T. Vicsek. Hierarchical group dynamics in pigeon flocks. Nature, vol.464, no.7290, pp.890-893, 2010. DOI:10.1038/nature08891 [8] Z. Y. Chen, H. T. Zhang, X. Chen, D. X. Chen, T. Zhou. Two-level leader-follower organization in pigeon flocks. Europhysics Letters, vol.112, no.2, pp.1-6, 2015. DOI:10.1209/0295-5075/112/20008 [9] B. Pettit, Z. Akos, T. Vicsek, D. Biro. Speed determines leadership and leadership determines learning during pigeon flocking. Current Biology, vol.25, no.23, pp.3132-3137, 2015. DOI:10.1016/j.cub.2015.10.044 [10] A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefaninib, M. Viale. Scale-free correlations in starling flocks. Proceedings of the National Academy of Sciences of the United States of America, vol.107, no.26, pp.11865-11870, 2010. DOI:10.1073/pnas.1005766107 [11] R. Lukeman, Y. X. Li, E. K. Leah. Inferring individual rules from collective behavior. Proceedings of the National Academy of Sciences of the United States of America, vol.107, no.28, pp.12576-12580, 2010. DOI:10.1073/pnas.1001763107 [12] A. Cavagna, I. Giardina, F. Ginelli. Boundary information inflow enhances correlation in flocking. Physical Review Letters, vol.110, no.16, pp.168107, 2013. DOI:10.1103/PhysRevLett.110.168107 [13] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, V. Zdravkovic. Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study. Proceedings of the National Academy of Sciences of the United States of America, vol.105, no.4, pp.1232-1237, 2008. DOI:10.1073/pnas.0711437105 [14] B. Pettit, A. Perna, D. Biro, D. J. Sumpter. Interaction rules underlying group decisions in homing pigeons. Journal of the Royal Society Interface, vol.10, no.80, pp.20130529, 2013. DOI:10.1098/rsif.2013.0529 [15] H. Pomeroy, F. Heppner. Structure of turning in airborne rock dove (Columba livia) flocks. The Auk, vol.109, no.2, pp.256-267, 1992. DOI:10.2307/4088194 [16] M. Yomosa, T. Mizuguchi, G. Vasarhelyi, M. Nagy. Coordinated behaviour in pigeon flocks. PLoS One, vol. 10, no. 10, Article number e0140558, 2015. [17] T. Niizato, Y. P. Gunji. Metric-topological interaction model of collective behavior. Ecological Modelling, vol.222, no.17, pp.3041-3049, 2011. DOI:10.1016/j.ecolmodel.2011.06.008 [18] H. X. Qiu, H. B. Duan. Receding horizon control for multiple UAV formation flight based on modified brain storm optimization. Nonlinear Dynamics, pp.1973-1988, 2014. DOI:10.1007/s11071-014-1579-7 [19] S. Kim, Y. Kim. Three dimensional optimum controller for multiple UAV formation flight using behavior-based decentralized approach. In Proceedings of International Conference on Control, Automation and Systems, IEEE, Seoul, South Korea, pp. 17–20, 2007. [20] P. P. Dai, C. L. Liu, F. Liu. Consensus problem of heterogeneous multi-agent systems with time delay under fixed and switching topologies. International Journal of Automation and Computing, vol.11, no.3, pp.340-346, 2014. DOI:10.1007/s11633-014-0798-1 [21] M. S. Mahmoud, Y. Q. Xia. The interaction between control and computing theories: New approaches. International Journal of Automation and Computing, vol.14, no.3, pp.254-274, 2017. DOI:10.1007/s11633-017-1070-2 [22] R. Freeman, R. Mann, T. Guilford, D. Biro. Group decisions and individual differences: Route fidelity predicts flight leadership in homing pigeons (Columba livia). Biology Letters, vol.7, no.1, pp.63-66, 2011. DOI:10.1098/rsbl.2010.0627 [23] C. W. Reynolds. Flocks, herds and schools: A distributed behavioral model. In Proceedings of the 14th Conference on Computer Graphics and Interactive Techniques, ACM, Anaheim, USA, pp. 25–34, 1987. [24] C. K. Hemelrijk, H. Hildenbrandt. Some causes of the variable shape of flocks of birds. PLoS One, vol. 6, no. 8, Article number e22479, 2011. [25] R. C. Rafaila, C. F. Caruntu, G. Livint. Centralized model predictive control of autonomous driving vehicles with Lyapunov stability. In Proceedings of the 20th International Conference on System Theory, Control and Computing, IEEE, Sinaia, Romania, pp. 663–668, 2016. [26] Y. N. Wang, Z. Q. Miao, H. Zhong, Q. Pan. Simultaneous stabilization and tracking of nonholonomic mobile robots: A Lyapunov-based approach. IEEE Transactions on Control Systems Technology, vol.23, no.4, pp.1440-1450, 2015. DOI:10.1109/TCST.2014.2375812 [27] X. Y. Zhang, H. B. Duan. Altitude consensus based 3D flocking control for fixed-wing unmanned aerial vehicle swarm trajectory tracking. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, vol.230, no.14, pp.2628-2638, 2016. DOI:10.1177/0954410016629692 [28] M. Favata. Conservative self-force correction to the innermost stable circular orbit: Comparison with multiple post-Newtonian-based methods. Physical Review D, vol.83, no.2, pp.1-26, 2011. DOI:10.1103/PhysRevD.83.024027 [29] Y. M. Deng, H. B. Duan. Control parameter design for automatic carrier landing system via pigeon-inspired optimization. Nonlinear Dynamics, vol.85, no.1, pp.97-106, 2016. DOI:10.1007/s11071-016-2670-z [30] H. B. Duan, X. H. Wang. Echo state networks with orthogonal pigeon-inspired optimization for image restoration. IEEE Transactions on Neural Networks and Learning Systems, vol.27, no.11, pp.2413-2425, 2016. DOI:10.1109/TNNLS.2015.2479117 [31] R. Dou, H. B. Duan. Levy flight based pigeon-inspired optimization for control parameters optimization in automatic carrier landing system. Aerospace Science and Technology, vol.61, pp.11-20, 2017. DOI:10.1016/j.ast.2016.11.012