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 bioinspired 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, leaderwingman, artificial potential field (APF), virtual structure and behaviorbased method. The leaderwingman 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 multiUAV formation^{[4]}, however when turning, UAVs at different positions fly at different speed and trajectories are asymmetric and redundant. Behaviorbased method solves the problem of multiUAV 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 bioinspired 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 multilevel leaderfollower strategy. Chen et al.^{[8]} analyzed the same global positioning system (GPS) data as Nagy′s but simplified the flock model and extracted a twolevel leaderfollower 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 scalefree 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 shortrange repulsion, intermediaterange alignment, and longrange attraction within circular zones as well as weak but significant frontsector 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 sidebyside 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, parallelpath and equalradius 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 bioinspired formation plan is analyzed utilizing the stability theory of ordinary differential equations with limited simplification of the model. Next comes Section 5, Levyflight based pigeon inspired optimization (LevyPIO) is adopted for controller parameter optimization. In Section 6, numerical simulations are implemented and results of multiUAV 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 missionThis 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 scalefree 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 multiUAV 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:


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:


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 constraintBecause 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 v_{x} is the speed along xaxis and v_{y} is yaxis speed. v is the airspeed, h is the altitude, γ is the flightpath angle and
$\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).
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 differenceCommunication ability is crucial for multiUAV formation, only with adequate redundancy, large range, channel availability, high signalnoise 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, UAV_{1} 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. UAV_{1} also receives position and velocity information from UAV_{2} and UAV_{n} 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 UAV_{1}, 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, earlywarning UAVs, etc. That is the reason why individual difference should be considered in the UAV formation patrol mission.
UAV_{2} and UAV_{n} are the second level leaders, both of them receive position and velocity information from UAV_{1}, the top leader, and from their own followers. UAV_{2} acquires information from UAV_{3}, while UAV_{n} from UAV_{n–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 UAV_{2} and UAV_{n} is bidirectional, both of them work as transmitters and receivers. Under the desirable condition, UAV_{2} serves as a relay which passes the leader information from UAV_{1} to UAV_{3}, also monitors the condition of UAV_{3}. The relationship between UAV_{n} and UAV_{n–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 foragingThere 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.
Doubleintegrator 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 UAV_{1}, it is driven by a navigation force
${\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
Follower UAVs employ the same form of force as UAV_{1} in (4) and (5). There is difference in destination
$\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 UAV_{1} and the midpoint of UAV_{2} and UAV_{n} in the standard formation. dis reflects the distance difference between the actual and desirable positions of UAV_{1},
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
When turning, coordinated turn is adopted
$\dot \chi = \frac{g}{v}\tan \phi $  (10) 
then n = cosγ – 10γ → γ = acos(ncosφ), the control variable is
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_{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, UAV_{i} is led by 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, UAV_{i} is subject to a force
$\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
As regards its cohesion behavior, UAV_{i} is attracted by the force pointing to UAV_{j} 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 r_{att} is the attraction radius,
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
Download:


Fig. 3. Different interaction zones 
Download:


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 flockingDuring 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, UAV_{1} flies towards the first waypoint, because the navigation force of UAV_{1} is stronger than other forces, UAV_{1} would overtake other UAVs and occupy a leader position in the front.
The desirable position of UAV_{2} is generated according to the present positions and velocities of UAV_{1}, UAV_{2} and UAV_{3}.
${\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
The desirable position of UAV_{3} is produced by its present position and a predicted position by UAV_{2}. UAV_{n} and UAV_{2} have symmetric positions, UAV_{n–1} and UAV_{3} are symmetric too. Attracted by navigation forces towards desirable positions, n UAVs can soon flock, coordinate velocities and form a stable formation.
Download:


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, UAV_{3} is located on the left side of UAV_{1}, after several steps, UAV_{1} has already finished turning with UAV_{3} located at the back. UAV_{3} 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 UAV_{1} 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
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, UAV_{1} overtakes both UAV_{2} and UAV_{n}, adjustment ends and level flight pattern begins to work. Adjustment criterion is the angles
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
The dynamics equation of UAV_{1} 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 x_{des} is reducible to a fixed value k and
$\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
V (x) = 0 if and only if x = 0.
V (x) > 0 if and only if 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,
As to follower UAVs, more simplification is needed because coupling among followers is much more complex, for example, the dynamics equation of UAV_{2} 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),
$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,
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 UAV_{i},
${\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,
$ \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
Considering UAV_{i}, 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 UAV_{i} in zone A, so do B, C, D and E. The sum of energy of UAV_{i} 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 r_{sep}, r_{ali}, r_{att} and r_{mid} are the interaction ranges in (11)–(15), △v_{k} is the relative difference of velocity of UAV_{i} and its neighbors in zone B. j, k, l, m and n represent neighbors of UAV_{i}. When V_{i} reaches its minimum, it is considered that UAV_{i} is flying under its most stable condition.
When x_{j} = r_{sep}, x_{l} = r_{ali}, x_{m} = x_{n} = r_{mid}, all UAVs are flying on the border of different interactions,
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, Levyflight strategy is adopted and two new operators are employed for improvement^{[31]}. The operators in LevyPIO are Levyflight based search operator and landmark operator with variable step length.
Levyflight, 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 X_{i} (t) is the pigeon i′s present location, X_{leader} is the global best location of the flock leader selected by the cost function in (33). X_{itemp} is the output of Levyflight based PIO. s and step are vectors with the same dimension as X_{i} (t) → s and step are vectors and they have the same dimension as X_{i} (t). N (0, σ^{2}) represents 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 X_{itemp} 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
Levyflight 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.
Download:


Fig. 6. Fitness curve of Levyflight based PIO 
According to Fig. 6, final fitness is f_{cost} = 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 analysisNumerical simulations of multiUAV 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 squareThe 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 RungeKutta method is employed for numerical simulation and time step is h = 0.01, 80 000 steps in total.
The formation result is shown in Fig. 7.
Download:


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 xaxis while compressed along yaxis, this process can be observed in the zone
When the leader UAV_{1} 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 UAV_{1} 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
Download:


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:


Fig. 9. Heading angle

In Fig. 9, heading angles of five UAVs are clearly shown. The heading angle of UAV_{1} is representative, its initial value approaches 0, that is consistent with the heading angle in Table 2 and means UAV_{1} flies towards the right side. It increases to
Download:


Fig. 10. Flight path angle γ of UAV_{1} 
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 UAV_{1} 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:


Fig. 11. Average distance among UAVs 
where
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:


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 UAV_{2}, 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 routeThe 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 RungeKutta method is employed and time step is h = 0.01, 150 000 steps in total.
The formation result is shown in Fig. 13.
Download:


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). UAV_{1} works as the leader and UAV_{2}, UAV_{3} 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:


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,
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
Download:


Fig. 15. Resolution of frontal attraction force 
If UAV_{j} is located in the frontal attraction zone of UAV_{i}, it is attracted towards UAV_{i} by
Bird flocking and foraging has enlightened people in many fields including multiUAV 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 multiUAV 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.
[1] 
M. A. Dehghani, M. B. Menhaj. Communication free leaderfollower formation control of unmanned aircraft systems. Robotics and Autonomous Systems, vol.80, pp.6975, 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.14531462, 2008. DOI:10.1016/j.simpat.2008.08.005 
[3] 
Y. Fu, X. K. Wang, H. Y. Zhu, Y. G. Yu. Parameters optimization of multiUAV 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 multirobot formation control. Robotics and Autonomous Systems, vol.56, no.4, pp.324333, 2008. DOI:10.1016/j.robot.2007.08.005 
[5] 
T. Balch, R. C. Arkin. Behaviorbased formation control for multirobot teams. IEEE Transactions on Robotics and Automation, vol.14, no.6, pp.926939, 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.306327, 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.890893, 2010. DOI:10.1038/nature08891 
[8] 
Z. Y. Chen, H. T. Zhang, X. Chen, D. X. Chen, T. Zhou. Twolevel leaderfollower organization in pigeon flocks. Europhysics Letters, vol.112, no.2, pp.16, 2015. DOI:10.1209/02955075/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.31323137, 2015. DOI:10.1016/j.cub.2015.10.044 
[10] 
A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefaninib, M. Viale. Scalefree correlations in starling flocks. Proceedings of the National Academy of Sciences of the United States of America, vol.107, no.26, pp.1186511870, 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.1257612580, 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.12321237, 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.256267, 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. Metrictopological interaction model of collective behavior. Ecological Modelling, vol.222, no.17, pp.30413049, 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.19731988, 2014. DOI:10.1007/s1107101415797 
[19] 
S. Kim, Y. Kim. Three dimensional optimum controller for multiple UAV formation flight using behaviorbased 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 multiagent systems with time delay under fixed and switching topologies. International Journal of Automation and Computing, vol.11, no.3, pp.340346, 2014. DOI:10.1007/s1163301407981 
[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.254274, 2017. DOI:10.1007/s1163301710702 
[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.6366, 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 Lyapunovbased approach. IEEE Transactions on Control Systems Technology, vol.23, no.4, pp.14401450, 2015. DOI:10.1109/TCST.2014.2375812 
[27] 
X. Y. Zhang, H. B. Duan. Altitude consensus based 3D flocking control for fixedwing unmanned aerial vehicle swarm trajectory tracking. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, vol.230, no.14, pp.26282638, 2016. DOI:10.1177/0954410016629692 
[28] 
M. Favata. Conservative selfforce correction to the innermost stable circular orbit: Comparison with multiple postNewtonianbased methods. Physical Review D, vol.83, no.2, pp.126, 2011. DOI:10.1103/PhysRevD.83.024027 
[29] 
Y. M. Deng, H. B. Duan. Control parameter design for automatic carrier landing system via pigeoninspired optimization. Nonlinear Dynamics, vol.85, no.1, pp.97106, 2016. DOI:10.1007/s110710162670z 
[30] 
H. B. Duan, X. H. Wang. Echo state networks with orthogonal pigeoninspired optimization for image restoration. IEEE Transactions on Neural Networks and Learning Systems, vol.27, no.11, pp.24132425, 2016. DOI:10.1109/TNNLS.2015.2479117 
[31] 
R. Dou, H. B. Duan. Levy flight based pigeoninspired optimization for control parameters optimization in automatic carrier landing system. Aerospace Science and Technology, vol.61, pp.1120, 2017. DOI:10.1016/j.ast.2016.11.012 