Robust Three-Dimensional Collision Avoidance for ...

Robust Three-Dimensional Collision Avoidance for Fixed-Wing Unmanned Aerial Systems Thomas J. Stastny∗ Gonzalo A. Garcia† and Shawn S. Keshmiri‡ University of Kansas, Lawrence, KS, 66044, USA This work considers a new ”morphing” potential field navigation technique, where potential fields morph to distribute potential cost to areas necessary for collision avoidance in the particular case of fixed-wing unmanned aerial systems, vehicles with high speed and high inertia. The original algorithm, developed for the two-dimensional lateral case, is extended to the third dimension through a layered potential plane approach in both longitudinal and lateral axes. The developed online navigation method is used in series with an integrated guidance and nonlinear model predictive controller, in which robustness is included through the augmentation of frequency dependent weighting matrices in the control structure. The combined approach is demonstrated in full nonlinear six-degrees-of-freedom simulation of a large fixed-wing unmanned aerial system showcasing successful online three-dimensional collision and obstacle avoidance and simultaneous fixed altitude and inclined trajectory following in a congested urban setting while subjected to high crosswind.

I.

Introduction

UTONOMOUS unmanned aerial systems (UASs) are becoming increasingly more relevant, making emersion into civil airspace imminent. Consideration of UASs in air traffic control scenarios is well documented with works such as [1, 2, 3]; however, domestic use of UASs entails low altitude operations possibly in proximity to urban areas where fixed obstacles (e.g. buildings, radio towers, trees, etc.) pose safety concerns. Payload and power constraints of UASs also require miniaturized detection systems with limited perception, e.g. miniaturized radar systems [4], further exacerbating the problem. UASs will require the capability to generate collision and obstacle free flight paths which a controller must tightly track while remaining robust to environmental disturbances (e.g. wind gusts) endemic to low altitude operations. A widely exploited approach to the collision avoidance problem in path planning is that of artificial potential field (APF) navigation methods and similar boid-esque algorithms [5, 6, 7, 8, 9, 10], where navigation in evolving environments is achieved by applying repulsive and attractive potentials to obstacles and goals, respectively, and typically commanding the agent to follow the negative gradient of the summed potential field. The approach is elegant in its simplicity utilizing simple relations to achieve desired emergent behavior in multi-agent systems. Its computational expense is also low. APF methods have been applied with much success for collision and obstacle free path planning and control in robotics [7, 11, 12], and more recently applied to multi-agent UAS guidance [1, 13, 14, 15]. The bulk of APF or boid-esque research has been conducted on relatively slow moving vehicles/robots [7, 11, 12, 16] and helicopters [17, 18, 19], or simply considers point-mass models where commanded accelerations are instantaneously realized by agents [2, 3, 8, 9, 10], i.e. higher order vehicle dynamic models are often not considered in the formulation. While adequate for many applications, as well as most air traffic conflict resolutions (where separation distances are very large), effectiveness wanes when considering high speed vehicles with high inertia operating in proximal settings. Beyond navigational concerns, robust collision avoidance in proximal conditions, and within the presence of environmental hazards, requires tight spatial and temporal tracking of highly dynamic and time-varying paths. Typical segment-based waypoint logic and linear control methods (controllers optimized about some defined trim point relying on linearized vehicle dynamics) found in commercial off-the-shelf (COTS) autopilot systems often fall short of such demands, suggesting the need for more advanced guidance and control schemes [20, 21]. In [22], a new morphing potential field (MPF) navigation approach was developed where APFs morph with respect to six-degrees-of-freedom (6DoF) physical characteristics specific to fixed-wing aircraft. The approach is used in a

A

∗ Graduate

Research Assistant, Department of Aerospace Engineering, [email protected], AIAA Student Member. Associate, Department of Aerospace Engineering, [email protected], AIAA Member. ‡ Associate Professor, Department of Aerospace Engineering, [email protected], AIAA Senior Member.

† Post-Doctoral

1 of 18 American Institute of Aeronautics and Astronautics

predictive way with an integrated guidance and nonlinear model predictive control (NMPC) scheme, demonstrating superior performance to classic APF methods in real-time collision and obstacle avoidance for the particular case of fixed-wing UAS and expanding the range of operations by the capability of considering a full nonlinear model of a large UAS in the controller. However, the navigation scheme is limited to fixed altitude path planning and lacks robustness to external disturbances in the controller. Though many UASs operations remain at a desired fixed altitude, for operations in or near urban areas, the ability to plan avoidance paths in the vertical dimension is a must. In this work, we extend the MPF navigation algorithm to the third dimension and also include a means of providing robustness in the NMPC. The resulting algorithm is capable of autonomous collision avoidance in congested conditions and maintenance of a fixed or inclined straight line time-varying trajectory through a full three dimensions. The approach is demonstrated through nonlinear 6DoF simulation of a large fixed-wing UAS.

II. A.

Nonlinear Modeling, Guidance, and Robust Control

Augmented Aircraft Model and Moving Waypoint Logic

A full nonlinear 6DoF rigid body description of the Meridian UAS [23] will be used (see [24, 25]). The Meridian UAS is a 1,100 pound (499 kg), long range, high endurance autonomous aircraft designed, built, and flight tested by the University of Kansas Aerospace Engineering Department for the National Science Foundation (NSF) Center for Remote Sensing of Ice Sheets (CReSIS) to provide an aerial platform for ice-penetrating radar developed for research in Polar Regions [23], see Fig. 1. Following the approach in [24, 26, 27, 28, 29, 30], the model output is augmented to include trajectory information in the form of time-varying waypoints, or ”moving points” characterized by a position and velocFigure 1: Meridian UAS flight test at NEEM, Greenland. ity. Using the current and previous moving point, o [k] and o [k − 1], a trajectory segment is generated at each time step by placing waypoints~a [k] and ~b [k] an arbitrarily large distance in both extremes of the segment defined between the moving points. It is the extended trajectory segment and waypoints~a [k] and ~b [k] the guidance logic (see Fig. 2) considers at each time step. n o o [k] : ~po [k] = [poN [k] , poE [k] , poH [k]]T ; vo [k] = [voN [k] , voE [k] , voH [k]]T (1) In Fig. 2, Iˆ is the local (flat Earth) inertial coordinate frame with components North N, East E, and Down D (Height H being the negative ⃗ + 1] �[� ⃗ + 2] �[� of Down). The aim of navigation scheme in Section III is the generation ⃗ [�] � ⃗ + 3] of these moving points throughout the NMPC horizon. �[� Abstractly (for brevity), the discretized augmented UAS model used in this work can be expressed as xk+1 = fd (xk , uk ) ∈ Rn with output �[� �[� + 2] �[� + 3] ⃗ + 3] r ek = he (xk ) ∈ R . The chosen discretization is Euler, for computational �[� + 1] efficiency, noting the NMPC is a numerical control technique. �-axis ⃗ + 2] �[�

B.

Robust Nonlinear Model Predictive Control

Embedding the augmented model and guidance logic described in Section II A into a NMPC scheme has shown promise for excellent aircraft control characteristics [22, 24, 26] and exhibits predictive qualities desirable for the present application, i.e. premonition in planning of avoidance maneuvers. The NMPC formulation also provides a means to incorporate a full nonlinear description of the aircraft system and its actuator constraints, allowing generation of control sequences adequate for a an expanded range of flight regimes (only possible with extensive

⃗ + 1] �[�

�[�]

�̂

�[� − 1]

�-axis

�-axis �-axis

�[�] ⃗

Figure 2: Waypoint generation logic.


gain scheduling in linear control techniques). The task of the NMPC will be to minimize the output errors e (from the augmented model output) within the prediction horizon, generating a robust, feasible control sequence while adhering to state output and control input constraints and keeping a tight trajectory following. Work completed in [24, 27] demonstrates the capability of simultaneous robustness to low frequency disturbances (e.g. wind, gusts), high frequency measurement noise, unmodeled aircraft dynamics, and over actuation of control surfaces through the use of frequency dependent weighting matrices integrated into the NMPC formulation. The approach is similar to mixed sensitivity techniques commonly used in H2 and H∞ control, see [31, 32, 33], relying on the characterization of some expected frequency range for different classes of disturbances and designing weighting matrices to penalize only these ranges for similar output signal. This allows concurrent disturbance rejection and noise mitigation without degradation of performance during the minimization of weighted model outputs z. An additional benefit of the robustness scheme is the low computational burden of the weighting matrices. In this work, only wind disturbances will be considered, though control actuation bandwidth is still desirable to regulate. Two weighting matrices W1 and W2 are incorporated into the NMPC formulation. Wind disturbances, characterized by low frequencies (0.1 − 10 Hz), are seen in the error outputs ek , and undesirable high frequency control actuation in uk . Design of dynamic weighting matrices for each set of signals is accomplished through formation of discrete state space models considering the targeted frequency range. The formulation is shown in Eq. 2 and Fig. 3.  xW1 = AW1 xW1 + BW1 e k k k+1 W1 : zW1 = CW1 xW1 + DW1 e k k  k (2) xW2 = AW2 xW2 + BW2 u k k k+1 W2 : zW2 = CW2 xW2 + DW2 u k

k

k

r T r T 1 1 W1 W2 r , xW2 = 1 2 xW , . . . , x ∈ R x , . . . , xW ∈ Rm are intermediate weighting matrix k k k k k h T iT W2 T 1 states, and new system output zk = zW ∈ Rr+m . The cost function is augmented for the new minimization k , zk objectives as follows: 1 where xW k =

N

J (z, u) =

∑ zTk+1 Qk+1 zk+1 + uTk Rk uk

(3)

k=1

� � Model + Guidance

�

�

�1 �2

��1 ��2

��

Optimization Algorithm

�

NMPC �̂ Figure 3: Robust NMPC formulation (left) and robustness weighting matrices singular values (right). where Qk+1 and Rk are positive definite and semi-positive definite weighting matrices, respectively. Figure 3 shows the singular weighting matrix design for each set of signals, e and u. Minimization of the cost function is completed using Sequential Quadratic Programming [34].


III. A.

Three-dimensional MPF Navigation Extension

Morphing Potential Field Formulation

In [22], the classical APF formulation was augmented to ”morph” for high speed, high inertia vehicles by inclusion of a morphing factor Γ into a Gaussian shaped potential function (see Eq. (4)) in a 2-D lateral case (i.e. NE inertial plane). The morphing is based on the angle of approach, magnitude of the moving point velocity relative to the object velocity~vor , where~vor =~vo − ~Vob j , and 6DoF kinematic aircraft constraints. The function extends the potential shape in appropriate directions distributing potential cost to areas necessary for adequate avoidance without the undesirable effects of amplitude or avoidance radius enlargement seen in the generic formulation. An additional reference shifting term ~S is also included in the distance norm as a means of further shaping the potential to avoid unnecessary levels of cost beyond the avoided obstacle. This term shifts the potential function origin away from the centroid of the obstacle.  !2   k~pob j −~po +~Sk  (4) mp f = exp −Γ   σ where σ is the desired avoidance radius, and ~pob j and ~po are defined as the 2-D lateral inertial position coordinates (considering only North and East components) of the object to avoid and the moving point, respectively. The vectors making up the distance norm in Eq. (4) can be represented as a single distance vector ~dc defined from the moving point position ~po to the shifted potential centroid ~c (i.e. ~dc = ~pob j −~po +~S), shown in Fig. 4 (a). The reference shifting term ~S is defined from the obstacle position ~pob j in the direction opposite the relative moving point velocity ~vor some distance ds = f (Rmin ), where ds ≤ Rmin and Rmin is the minimum turning radius of the aircraft (see Eq. (5)) as a function of the speed and maximum allowed bank angle φmax . In this work, ds is simply expressed as a fraction of Rmin , though more complex relations could be devised. Shifting the potential centroid acts to redistribute some of the girth of the potential cost toward an approaching moving point, also leaning out the cost behind an obstacle. Rmin =

k~vor k2 g tan φmax

(5)

The morphing factor is expressed in Eq (6) Γ = (1 − Γ0 ) sin2 ηv + Γ0 ;

Γ0 ∈ (0, 1] , ηv ∈ [−π, π]

(6)

where ηv is the angle between the moving point relative velocity vector ~vor and the distance vector ~dc and Γ0 is the extension term of the morphing factor. Note that the magnitude of the extension term determines the extent of morphing while ηv determines the direction of morphing and contour shape about the shifted center. Γ0 is defined as a smooth piecewise function such that the potential contour line corresponding to the avoidance radius σ extends to σ + Rmin when ηv = 0 (i.e. head-on approach), remains at the nominal radius σ when |ηv | = 90◦ , and extends back to the initial contour line when |ηv | = 180◦ (see Eq. (7) and Fig. 4 (a) and (b)). As k~vor k increases, Γ0 decreases and the morphing distance increases. When ~vor = 0, Γ0 = 1 and ds = 0; no morphing or shifting occurs and the function reduces to the general potential formulation.  2  σ  , |ηv | < π/2 σ +Rmin −ds 2 Γ0 = (7)   σ , else σ +ds Taking the gradient of mp f with respect to distance vector components dcN , dcE (see Eq. 8), negating, and normalizing (i.e. forming a unit vector, except where a gradient of zero gives a zero vector) will produce a vector field usable in the moving point planning algorithm, see [22]. Figure 4 (b) shows an example of a morphed vector field. ( ) 2 ~ (1 − Γ0 ) cos ηv k~dc k o ~vr ∇mp f = mp f 2 dc − (8) σ k~vor k To avoid uncoordinated avoidance in cooperative flight, Chang et. al. in [35] inspired the utilization of a gyroscopic vector field (see Eq. (9), where dcN , dcE are the North and East inertial components of the distance vector ~dc ). gyro = [−dcE , dcN ]


(9)

𝑝𝑓 𝑐𝑜𝑛𝑡𝑜𝑢𝑟 𝑚𝑝𝑓 𝑐𝑜𝑛𝑡𝑜𝑢𝑟

𝜎

⃗ 𝑜𝑏𝑗 𝒑

𝜎 + 𝑅𝑚𝑖𝑛

⃗𝑺 ⃗ 𝒄 𝜂𝑣 ⃗𝒅𝑐

⃗ 𝒐𝑟 𝒗 ⃗𝒐 𝒑

(a)

(b)

Figure 4: Contour lines and vector definitions for a regular potential function and morphed potential function about a fixed obstacle (a). Vector field generated by the gradient of a morphed potential field (b). Parameters used: ~pob j = {750, 500} f t, ~Vob j = {0, 0} f t/s,~vor =~vo = {169, 0} f t, σ = 150 f t, φmax = 60◦ , ds = Rmin /3.

Adding a minimally weighted gyroscopic field to the morphing potential gradient in some obstacle avoidance and uncooperative vehicle-to-vehicle avoidance scenarios also produced some desirable effects in [22]. However, in these scenarios, the gyroscopic spin is directed clockwise or counter-clockwise based on the current moving point heading with respect to the distance vector ~dc (similar to the split/rejoin logic seen in [35]). Negative and positive values of the approach angle (±ηv ) instigate clockwise and counter-clockwise spin, respectively. During normalization of the morphed vector field, a weighted sum of the two fields, produces a combined field ∇mp f = %gyro ·gyro−%mp f ·∇mp f , where %gyro and %mp f sum to 100%. B.

Layered Potential Planes

In this section, the morphing potential formulation is extended to the third spacial dimension by ”layering” poten0 , see Eq. (10). A layering term tial planes vertically (in the H-axis) through manipulation of the avoidance radius σLAT j `LAT is introduced as a function of the height difference ∆pH = |poH − pob H | and a longitudinal avoidance radius σLON . Note the designation LAT identifies the lateral plane, in which the formulation still resides; and LON, the longitudinal plane, wherein the layers are ”stacked”. This distinction will be important when considering potentials in both lateral and longitudinal planes, to be developed throughout this section. 0 σLAT = `LAT σLAT

(10)

The morphing potential function remains as it was defined in Section III A; however, each σ in the formulation is 0 . Recall the extension term Γ in Γ, Eq. (6), is also a function of σ , and thus in the new definition, replaced with σLAT 0 0 . ΓLAT is a function of σLAT 

2  

~

  

dcLAT   (11) (mp f )LAT = exp −ΓLAT  0   σLAT   where ~dcLAT takes on the same definition from Section III A as a vector of only N and E components. The layering term `LAT will determine what lateral avoidance radius is considered at each height with respect to the object to avoid. Choice of a function for `LAT is mostly determined by any mission constraints and/or obstacle/vehicle geometry, allowing for a fair amount of flexibility in the definition of a volume to be avoided. Here, a spherical definition will be used.


s `LAT =

∆pH 1− σLON

2 (12)

where (mp f )LAT must be redefined as a piecewise function in order to describe a finite volume (and avoid complex values of `LAT !).  ( 2 )  ~dc  k k LAT exp −ΓLAT , |∆pH | < σLON 0 σLAT (mp f )LAT = (13)   0, else Terming the present form of `LAT as ”spherical” suggests a representation of a spherical volume; however, bear in mind that defining σLAT and σLON at different radii provides a means to create ellipsoidal-esque volumes as well. Another variation could be the choice of a piecewise `LAT function to extend the infinite cylinder (seen in the original 2D formulation) for some range (of course, entailing a bound adjustment in the current piecewise (mp f )LAT definition), then tailor off the volume to a point, perhaps to describe a skyscraper or radio tower. Figure 5 (left) shows several layers of potential planes for an arbitrary obstacle. At each plane, potential cost is calculated for the given avoidance radius until the radius reaches zero, i.e. the piecewise bound of (mp f )LAT is reached. Figure 5 (right) shows the same obstacle; however, as a complementary visualization, a ”contour shell” representation displays the variation in potential cost throughout the entire range of the finite volume in the H-axis.

70 65 60

Height

55 50 45 40 150

35 100

30 20

50 40

60

80

0

North

East

Figure 5: Morphed potential contour layers (left). Potential contour shells (right). Reformulation of the morphing potential cost equation entails redefinition of the gradient function. A piecewise 0 , as definition of (∇mp f )LAT should mirror the bounds of (mp f )LAT and all σ in the formulation replaced with σLAT in the cost function. As the lateral gradient is only a function of N and E position components, inclusion of the height term will have no effect on differentiation. Also, gyroscopic augmentation to the formulation remains identical to that seen in Section III A. The final output must be a 3-D vector in the inertial coordinate frame, the resulting lateral morphing gradient can be expressed as    (∇mp f )LATN       (∇mp f )  , |∆pH | < σLON  LATE  (∇mp f )LAT I = (14)   0    ~ 0, else


C.

Longitudinal Formulation

While layering lateral potential planes provided some improvement to the representation of obstacles in three dimensions, the formulation still only allows for avoidance planning in the occupied 2-D lateral plane. This section extends the morphing potential formulation to a longitudinal plane, allowing generation of avoidance vectors with an H-axis component. Lateral components N and E are projected onto a longitudinal N v -axis, defined along the moving point velocity vector and perpendicular to the H-axis, see Fig. 6. For this formulation, longitudinal 2-D position and velocity vectors are defined as follows. ~poLON = [poN v , poH ]T ~voLON = [voN v , voH ]T h i j ob j ob j T ~pob LON = pN v , pH h iT ~Vob j = V obv j ,V ob j N H LON

(15)

poN v = poN cos ξv − poE sin ξv voN v = voN cos ξv − voE sin ξv ob j ob j j pob N v = pN cos ξv − pE sin ξv ob j ob j VN v = VN cos ξv −VEob j sin ξv ξv = arctan (voE /voN )

(16)

where

𝐻

𝐻

𝐿𝑂𝑁𝐺𝐼𝑇𝑈𝐷𝐼𝑁𝐴𝐿 𝑃𝐿𝐴𝑁𝐸 𝑁

⃗ 𝑜𝑏𝑗 𝒑

𝑰̂

⃗ 𝑜𝑏𝑗 𝒑 𝐿𝑂𝑁

⃗𝒐 𝒗 ⃗𝒅𝑐 𝐿𝑂𝑁

⃗𝒐 𝒑 𝜉𝑣

⃗ 𝒄 𝑁𝑣

𝐸

Figure 6: Longitudinal potential plane. From this point, progression through the morphing potential field development, for the most part, mirrors the original formulation and layered extension in the lateral plane. This, with the distinction of longitudinal components, as oppose to their lateral counterparts. The potential cost is defined as follows  ( 2 )  ~dc  k k LON exp −ΓLON , |∆p⊥ | < σLAT 0 σLON (mp f )LON = (17)   0, else j 0 o ~ where ~dcLON =~pob LON −~pLON − SLON and σLON = `LON σLON . The longitudinal layering term will be chosen similar to the lateral case, i.e. s ∆p⊥ 2 `LON = 1 − (18) σLAT


where ∆p⊥ is defined as the magnitude of the distance between the object to avoid and the longitudinal plane. ~vo ×H ¯ · ~pob j −~po ∆p⊥ = ¯ k~vo ×Hk

(19)

¯ is the unit vector corresponding to the H-axis. Note the position and velocity vectors used here are the full where H 3-D expressions, unless otherwise specified with a LON tag. ~vor The longitudinal reference shifting term is defined ~SLON = −ds LON , where ~vor =~vo − ~Vob j is the lon o

~v

LON

LON

LON

rLON gitudinal relative velocity and ds = f RminLON , again chosen as a constant fraction of RminLON for the present work. RminLON is the minimum distance in the N v -axis required to climb or descend (solely) to avoid an object with longitudinal avoidance radius σLON . Here, RminLON will be defined as a function of a fixed conservative flight path angle γmax > 0.

RminLON =

σLON tan γmax

(20)

The morphing factor ΓLON retains the same structure as in the lateral case with longitudinal parameter substitutions. ΓLON = (1 − Γ0 ) sin2 ηvLON + Γ0 ;

Γ0 ∈ (0, 1] , ηvLON ∈ [−π, π]

where ηvLON is defined between ~dcLON and~vrLON and  2 0 σLON   0 , |ηvLON | < π/2 σLON +RminLON −ds Γ0 = 0 2   0σLON , else σ +d LON

(21)

(22)

s

With potential cost defined in the longitudinal plane, the gradient function can be derived with respect to N v and H position components.

 

~

  (1 − Γ ) cos η d

v c 0 LON LON 2 o ~dc −

~ v (23) (∇mp f )LON = (mp f )LON 0 r LON LON

~vor  σLON  LON As in the lateral case, gyroscopic augmentation to the longitudinal mp f formulation is conducted inthe same man ner as originally developed. However, the vector is defined with longitudinal components, i.e. gyro = −dcH , dcN v . After calculation of (∇mp f )LON , a direction vector in the inertial reference frame must be extracted. The longitudinal morphing potential gradient is expressed in the inertial coordinate frame as    cos ξ (∇mp f )  v LONN v       (∇mp f )  , |∆p⊥ | < σLAT  LONN v sin ξv  (24) (∇mp f )LON I =   (∇mp f )LONH    ~ 0, else D.

Combined Approach

With both lateral and longitudinal layered formulations capable of generating avoidance vectors in their respective planes, the thought naturally occurs to combine the results. However, considering the complex vehicle dynamics, aerodynamics, and performance characteristics of UASs on top of specific mission objectives, payload constraints, and varying priorities, the ”optimal” manner of such a merger is not immediately evident. For example, collision avoidance based purely on aircraft performance might weigh the lateral and longitudinal capabilities of the given platform as metrics for a ”path of least resistance” decision logic to blend each vector plane. On the other hand, specific payload constraints (e.g. remote sensing with an altitude sensitive instrument) may dictate avoidance to first consider all possible paths in the lateral plane before a last resort to diving or climbing. The possibilities are wide ranging and extremely scenario specific, and it will not be the aim of this work to determine the optimal settings. However, to explore, briefly, the notion, a simple weighting logic will be introduced to show, as a proof of concept, that combined


avoidance maneuvers can be achieved using the developed algorithm. Both lateral and longitudinal potential costs and normalized gradients are calculated separately, then combined in a weighted sum, and re-normalized. (mp f )3D = %LON (mp f )LON + %LAT (mp f )LAT

(25)

(∇mp f )3D = %LON (∇mp f )LON I + %LAT (∇mp f )LAT I

(26)

where %LON and %LAT sum to 100%. E.

Trajectory Following Vector Field

As a base, the navigation unit generates moving points to approach a straight trajectory defined by a desired heading ξtra j and coordinate point ~ptra j (located at any point on the desired trajectory path). An arctangent function is used to generate a vector field steering moving points asymptotically to the path similar to the approach seen in [36]. Though many UAS operations remain at a desired fixed altitude, the ability to plan moving points through a full three-dimensions to achieve some trajectory is, nevertheless, valuable. This section will extend the trajectory following vector field to a 3-D formulation capable of steering moving points toward a straight trajectory, whether at some incline or fixed altitude, in a combination of both lateral and longitudinal planes. For the lateral case, the vector field is calculated using the following equation. ξapp =

n h io 2 j tra j o ξappmax arctan Kξ poN − ptra sin ξ − p − p cos ξ tra j tra j E N E π

j j are the North and East position coordiand ptra where ptra E N nates of any arbitrary point on the straight trajectory path, ξapp is the approach heading relative to the trajectory path, ξappmax is the maximum approach heading (user defined), and Kξ is a gain adjusting the shape of approach vectors. Including the H-axis in the formulation will require the definition of a longitudinal plane in which moving points and trajectory parameters, similar to the lateral case, can be projected. A new N t -axis is defined along the track, see Figure 7. Inclined trajectory capability is included via a user designated flight path angle γtra j . As in the lateral case for heading, a maximum approach flight path angle is defined γappmax as well as a function gain Kγ . Note γappmax need not correspond to a physical flight path angle constraint, more just the shaping of the arctangent function, i.e. how steep of a ”descent” toward the given track the vector field will exhibit. As one indirect constraint on the aircarft longitudinal dynamics, maximum flight path angle constraint for the planned moving points will be described in the following section, limiting any unachievable trajectory vectors generated by the field. The longitudinal approach angle is calculated as

γapp =

(27)

𝐻

𝐻

𝐿𝑂𝑁𝐺𝐼𝑇𝑈𝐷𝐼𝑁𝐴𝐿 𝑃𝐿𝐴𝑁𝐸 𝑁

𝑰̂

𝜉𝑡𝑟𝑎𝑗 ⃗ 𝑡𝑟𝑎𝑗 𝒑

𝛾𝑡𝑟𝑎𝑗 𝑁𝑡

𝐸

Figure 7: Longitudinal trajectory following plane.

n h io 2 j j γappmax arctan Kγ poNt − ptra sin γtra j − poH − ptra cos γtra j t H N π

(28)

j tra j tra j where poNt = poN cos ξtra j − poE sin ξtra j and ptra N t = pN cos ξtra j − pE sin ξtra j . The combined 3-D direction vector

can now be constructed.   cos (γapp + γtra j ) cos (ξapp + ξtra j )   tra j =  cos (γapp + γtra j ) sin (ξapp + ξtra j )  sin (γapp + γtra j )


(29)

F.

Waypoint Planning Constraints

The trajectory vectors and avoidance vectors generated by the potential and arctangent functions in the previous sections are summed in an adaptive, prioritized weighting scheme based on their potential costs, see [22]. In the event a resultant vector is generated with a desired heading change greater than that capable of the aircraft, a maximum heading change ∆ξmax is defined in Eq. (30) based on standard aircraft steady state turning performance (see [37]) and the path planning sample time. Flight path angle constraints must also be imposed to avoid unachievable (or dangerous) climbs and dives. Similar to the maximum heading change defined, the change in flight path angle for a 2-g symmetric pull-up can be expressed discretely as a function of g and vod (desired moving point velocity) in the path planning time scale and used as another indirect constraint on longitudinal aircraft dynamics. ∆ξmax =

g tan φmax ∆t pp vod

(30)

g ∆t pp vod

(31)

∆γmax =

A further constraint is also necessary in the longitudinal plane, as the value of γ should not exceed some practical limits at a given velocity due to the performance of a given aircraft. For simplicity in this work, and due to a fixed speed constraint, a conservative constant flight path limit γ~Rmax (not to be confused with γmax in Section III C) will be set, though more detailed and variable limits could be derived from the aircraft physics. The moving point planner must saturate rise and descent vectors at this maximum angle.

IV.

Simulation and Discussion

This section presents several examples of the 3-D collision avoidance formulation in 6DoF simulation of the Meridian UAS. A fixed desired speed constraint is set to vod = 100 kts. The maximum allowable flight path angle in the moving point planner is set to γ~Rmax = 8◦ . A.

Obstacle Avoidance with Layered Lateral MPF

Seven UASs are initialized at varying heights with respect to a spherical obstacle placed directly in the path of the desired trajectories for several of the aircraft. The UASs must avoid the obstacle only using lateral maneuvers (as these are the only avoidance vectors the lateral MPF will generate). However, using the layered formulation, it can be seen in Fig. 8 that the avoidance path for each UAS differs with respect to the lateral radius of the sphere at each height, deviating only as much as necessary until completely ignoring the obstacle when flying directly overhead. Parameters used by all UASs in this simulation are shown in Table 1. Table 1: Simulation parameters for layered lateral obstacle avoidance. Parameter σLAT ξtra j ξappmax Kξ %gyro %LAT

B.

Value 1.4 · 50 f t 0◦ 45◦ 0.01 0.1 1

Parameter σLON γtra j γappmax Kγ ds %LON

Value 1.2 · 150 f t 0◦ 8◦ 0.005 Rmin · 3/4 0

Obstacle Avoidance with Layered Longitudinal MPF

In this simulation, seven UASs are initialized at varying lateral distances with respect to an ellipsoidal obstacle placed directly in the path of their desired trajectories, all of which are in the lateral plane. Unlike the layered lateral case, the UASs must avoid the obstacle only using longitudinal maneuvers derived from the layered longitudinal MPF formulation. Figure 9 (left) shows similar behavior to the lateral case, only in the longitudinal plane. Each UAS avoids by climbing over the obstacle, though the altitude to which each aircraft must climb varies with respect to the size of the obstacle ”slice” in each particular longitudinal layer. Parameters used by all UASs in this simulation are shown 10 of 18 American Institute of Aeronautics and Astronautics

Obstacle UAS 1 Trajectory

3200

UAS 2 Trajectory UAS 3 Trajectory UAS 4 Trajectory UAS 5 Trajectory

3000

UAS 6 Trajectory UAS 7 Trajectory UAS 1 Flight Path UAS 2 Flight Path

2800

UAS 3 Flight Path UAS 4 Flight Path UAS 5 Flight Path

Obstacle

UAS 6 Flight Path

2600

UAS 1 Trajectory

1000

North (ft)

UAS 7 Flight Path

UAS 2 Trajectory UAS 3 Trajectory

900

UAS 4 Trajectory UAS 5 Trajectory

2400

UAS 6 Trajectory

800

UAS 7 Trajectory UAS 1 Flight Path UAS 2 Flight Path

Height (ft)

700 2200

UAS 3 Flight Path UAS 4 Flight Path

600


500

2000

400 1800

3000

300

2500 −200 −100

0

100 200 East (ft)

300

400

500

2000

200 −500

0

500

North (ft)

East (ft)

Figure 8: Obstacle avoidance using layered lateral MPF navgiation. Top view (left), 3-D view (right).

in Table 2. Note the longitudinal trajectory following parameters are more moderate than their lateral counterparts. This is due to the sensitivity of the Meridian dynamics in climbs and dives. A shallow descent back to the desired trajectories is planned due to the strict limits set in the trajectory following parameters. Avoiding too steep of a descent is important as the Meridian has no control surfaces capable of decelerating the aircraft in a dive; reduction in throttle setting, bounded by the NMPC input constraint to a minimum of 10% in the current formulation, is the only form of deceleration possible in this scenario. In such a dive, the aircraft will gain too much speed which can lead to destabilization, making control of the UAS very challenging, or impossible. Table 2: Simulation parameters for layered longitudinal obstacle avoidance. Parameter σLAT ξtra j ξappmax Kξ %gyro %LAT

C.

Value 1.2 · 150 f t 0◦ 45◦ 0.01 0.1 0

Parameter σLON γtra j γappmax Kγ ds %LON

Value 1.4 · 50 f t 0◦ 8◦ 0.005 Rmin · 3/4 1

Obstacle Avoidance with Combined Approach

To demonstrate the capability of blended lateral and longitudinal avoidance planning in 3-D, a spherical obstacle is placed in the path of three UASs, all with identical initial conditions. UAS 1 only considers lateral avoidance, UAS 2 only considers only longitudinal avoidance, and UAS 3 weights each formulation equally. Figure 9 (right) shows each UAS avoiding the obstacle using three different approaches. Parameters used are shown in Table 3.


700

680 800

660 750

640 700

Height (ft)

620 650 Height (ft)

Obstacle UAS 1 Trajectory UAS 2 Trajectory

600

UAS 3 Trajectory

Obstacle UAS 1 Trajectory UAS 2 Trajectory UAS 3 Trajectory UAS 1 Flight Path UAS 2 Flight Path UAS 3 Flight Path

600

580

UAS 4 Trajectory UAS 5 Trajectory UAS 6 Trajectory

550

UAS 7 Trajectory

560


500

540


450

UAS 7 Flight Path

400

2500

520

3000

−300

−200

−100

500 −100

2000

0

100

200

300

−80

−60

−40

−20

0

20

40

60

80

2800 2700 2500 100

North (ft)

North (ft) East (ft)

East (ft)

Figure 9: Obstacle avoidance using layered longitudinal (left) and combined layered longitudinal and lateral (right) MPF navgiation.

Table 3: UAS 3 simulation parameters for combined longitudinal and lateral obstacle avoidance.

UAS 1 UAS 2 UAS 3 All

D.

Parameter σLAT %LAT σLAT %LAT σLAT %LAT ξtra j ξappmax Kξ % gyro

Value 1.4 · 50 f t 1 1.2 · 50 f t 1 1.4 · 50 f t 0.5 0◦ 45◦ 0.01 0.1

Parameter σLON %LON σLON %LON σLON %LON γtra j γappmax Kγ ds

Value 1.2 · 50 f t 0 1.4 · 50 f t 0 1.4 · 50 f t 0.5 0◦ 8◦ 0.005 Rmin · 3/4

Robust Collision and Obstacle Avoidance in a Congested, Unstructured Environment

As a final demonstration, the 3-D planning algorithm along with robust NMPC is put to test in a congested and unstructured environment. The simulation starts with a UAS in a climb (following an inclined desired trajectory), perhaps representative of a climb after take off. However, the climb trajectory is obstructed by the top corner of a building which the planning algorithm must avoid, a scenario not so removed from the realities of human error in flight path planning and autonomous aircraft operation [38]. Shortly after clearing the first building, the desired trajectory is changed to a fixed altitude slightly higher than the top of the incline. Moving points must be planned up to this new trajectory while avoiding a second building and simultaneously experiencing a 10 kts wind step in a spatially and temporally varying wind field. Beyond the second building, a last moving obstacle, and uncooperative airship, must be avoided. The goal of the simulation is to test all elements of the combined navigation, guidance, and robust control algorithms in a 3-D setting. The desired trajectory considered by the navigation algorithm is dependent on the N component of the moving point position at any given time. I.e., if poN < 2600 f t, then ~ptra j = {2600, 0, 600}T , ξtra j = 0◦ , and γtra j = 6◦ ; else ~ptra j = {0, 0, 230}T , ξtra j = 0◦ , and γtra j = 0◦ . A spatially and temporally varying wind field is generated with a 1 − cosine step shape, based on United States Military Specification MIL-F-8785C. An eastward crosswind experienced by the UAS is modeled with the following dynamics.


   wE0 ,

wE1 −wE0 wE =  2

poE < pw E h i 1 − cos lπr (poE − pw ) + wE0 , E

  wE1 ,

poE < pw E + lr

(32)

else w pw E = pE0 + wE1 t

where wE0 = 0 kts and wE1 = 10 kts are the minimum and maximum magnitudes of the wind field, pw E is the position of the wind front (foot of the wind step), pw = 304 f tis the initial position (at t = 0 s) of the wind front, t is the current E0 flight time (assuming initialization at t = 0 s), and lr = 20 f t is the rise time length of the wind step. Building 1 is represented as an infinite cylinder extending in either direction of the longitudinal plane, i.e. no layering effects, positioned to cover the top two corners of the roof in the H-N plane. This choice was made in lieu of introducing new decision logic for lateral vs. longitudinal MPF weighting. It is assumed, for the present scenario, the UAS mission objective is to stay as close to the desired trajectory as possible while avoiding obstacles, thus the longitudinal avoidance path deviates less from the trajectory than a lateral path. Building 2 is similarly represented as an infinite cylinder in the lateral plane, using the original MPF formulation, where again the building is tall enough to disregard any longitudinal avoidance. The airship is modeled after the Goodyear GZ-20 Blimp. An ellipsoidal representation is used in the layered longitudinal MPF formulation based on specifications found in [39]. Figures 10-12 show the simulation. In Figure 11, UAS successfully clears the top of Building 1. Once the desired trajectory switches, the UAS executes a climb while laterally avoiding Building 2 and remaining robust to the step in the wind field. Figure 12 shows the moving points generated during the maneuver and the effect of the wind disturbance on the aircraft. The UAS quickly returns to the moving points and executes a longitudinal dive to avoid the passing airship, all the while remaining robust to the now constant crosswind.

𝐵𝑢𝑖𝑙𝑑𝑖𝑛𝑔 2


Figure 10: Robust collision and obstacle avoidance in a congested and unstructured environment. States and controls from the simulation are shown in Figure 13. Despite the demanding conditions, states are kept within reasonable ranges and control deflections rarely saturate, though the NMPC input constraints successfully limit actuation beyond the maximum value prescribed.


𝑡 = 15.0 𝑠 𝑡 = 10.0 𝑠 𝑡 = 5.0 𝑠


𝑡 = 0.0 𝑠

𝐵𝑢𝑖𝑙𝑑𝑖𝑛𝑔 2 𝑡 = 25.0 𝑠

𝑡 = 20.0 𝑠

Figure 11: Longitudinal avoidance (side view).

𝐵𝑢𝑖𝑙𝑑𝑖𝑛𝑔 2 𝑡 = 15.0 𝑠

𝑡 = 20.0 𝑠 𝒘 North (ft) Figure 12: Robust lateral avoidance (top view).


𝑡 = 30.0 𝑠

δT (%)

α (deg)

8 6 4 2 0 −2 0

5

10

15 20 Time (s)

25

30

1 0.8 0.6 0.4 0.2

35

2 0 −2 5

30

35

0

5

10

15 20 Time (s)

25

30

35

−2

10

15 20 Time (s)

25

30

35

δA (deg)

φ (deg)

0

φcmd φ 0

5

10 0 −10

10

15 20 Time (s)

25

30

35

δR (deg)

10 5 θcmd

0

θ 0

5

10

15 20 Time (s)

25

30

168

5

10

15 20 Time (s)

25

30

35

0

5

10

15 20 Time (s)

25

30

35

15 20 Time (s)

25

30

35

0 −10

35

170

0

10

˙ (ft/s2 ) W

θ (deg)

25

20

−50

T

15 20 Time (s)

0

50

V (ft/s)

10

−4 0

−5

5

2 δE (deg)

β (deg)

4

0

−20 −40 −60

166 0

5

10

15 20 Time (s)

25

30

35

0

5

10

Figure 13: States and control deflections from robust 3-D collision and obstacle avoidance.


Table 4: Simulation parameters for robust 3-D collision and obstacle avoidance.

Building 1

Building 2

Airship

Tra jectory 1

Tra jectory 2

Parameter σLAT %LAT %gyro σLAT %LAT % gyro σLAT %LAT %gyro ~Vairship ξtra j ξappmax Kξ ξtra j ξappmax Kξ

Value −− 0 0.1 1.2 · 135 f t 1 0.1 1.3 · 97 f t 0 0 {0, 40, 0}T f t/s 0◦ 45◦ 0.01 0◦ 45◦ 0.01

V.

Parameter σLON %LON ds σLON %LON ds σLON %LON ds

Value 1.2 · 170 f t 1 Rmin · 1/3 −− 0 Rmin · 1/2 1.3 · 30 f t 1 Rmin · 1/3

γtra j

6◦ 45◦ 0.005 0◦ 45◦ 0.005

γappmax Kγ γtra j γappmax Kγ

Conclusion

In this work, a new morphing potential field (MPF) navigation technique for fixed-wing unmanned aerial systems (UASs) has been extended to the third dimension for online path planning (in the form of time-varying ”moving” waypoints) including simultaneous three-dimensional (3-D) collision avoidance and fixed altitude or inclined trajectory following. The navigation is conducted in a predictive way planning moving points through a horizon which an integrated guidance and nonlinear model predictive control (NMPC) algorithm is tuned to track. Robustness is included in the NMPC formulation through frequency dependent weighting matrices toward the aim of concurrent disturbance rejection and performance optimization. The combined technique is tested in full nonlinear six-degrees-of-freedom (6DoF) simulation of a large UAS. Fixed obstacle avoidance with both lateral, longitudinal, and blended evasion planning is demonstrated. Successful time-varying trajectory following (both fixed altitude and inclined) is achieved while the aircraft is subjected to high crosswind in a simulation within a congested urban environment containing multiple fixed and moving obstacles.

Acknowledgments This work was completed with funding from the Paul G. Allen Family Foundation (PGAFF) grant KUEA#40956 and the National Science Foundation (NSF) Center for Remote Sensing of Ice Sheet (CReSIS) grant ANT-0424589. The authors would like to thank PGAFF, CReSIS, and NSF for their support.

References 1

Sigurd, K. and How, J., “UAV Trajectory Design Using Total Field Collision Avoidance,” AIAA Guidance, Navigation, and Control Conference and Exhibit, 2003.

2

George, J. and Ghose, D., “A Reactive Inverse PN Algorithm for Collision Avoidance Among Multiple Unmanned Aerial Vehicles,” American Control Conference, 2009, pp. 38903895.

3

Archibald, J. K., Hill, J. C., Jepsen, N. A., Stirling, W. C., and Frost, R. L., “A Satisficing Approach to Aircraft Conflict Resolution,” IEEE Transactions on Systems, Man, and Cybernetics, Part C: Applications and Reviews, Vol. 38, No. 4, 2008, pp. 510–521.


4

Shi, L., Allen, C., Ewing, M., Keshmiri, S., Zakharov, M., Florencio, F., Niakan, N., and Knight, R., “Multichannel Sense-and-Avoid Radar for Small UAVs,” IEEE/AIAA 32nd Digital Avionics Systems Conference (DASC), 2013, pp. 601–610.

5

Khatib, O., “Real-time Obstacle Avoidance for Manipulators and Mobile Robots,” IEEE International Conference on Robotics and Automation. Proceedings, 1985, pp. 500505.

6

Reynolds, C. W., “Flocks, Herds and Schools: A Distributed Behavioral Model,” Proceedings of the 14th annual conference on Computer graphics and interactive techniques, ACM, New York, NY, USA, 1987, pp. 2534.

7

Leonard, N. E. and Fiorelli, E., “Virtual Leaders, Artificial Potentials and Coordinated Control of Groups,” Proceedings of the 40th IEEE Conference on Decision and Control, 2001, pp. 29682973.

8

Olfati-Saber, R., “Flocking for Multi-agent Dynamic Systems: Algorithms and Theory,” IEEE Transactions on Automatic Control, Vol. 51, No. 3, 2006, pp. 401420.

9

Su, H., Wang, X., and Lin, Z., “Flocking of Multi-agents with a Virtual Leader Part I: with a Minority of Informed Agents,” 46th IEEE Conference on Decision and Control, 2007, pp. 29372942.

10

Su, H., Wang, X., and Lin, Z., “Flocking of Multi-agents with a Virtual Leader Part II: with a Virtual Leader of Varying Velocity,” 46th IEEE Conference on Decision and Control, 2007, pp. 14291434.

11

Keeter, M., Moore, D., Muller, R., Nieters, E., Flenner, J., Martonosi, S. E., Bertozzi, A. L., Percus, A. G., and Levy, R., “Cooperative Search with Autonomous Vehicles in a 3D Aquatic Testbed,” American Control Conference (ACC), 2012, pp. 31543160.

12

Nguyen, B. Q., Chuang, Y.-L., Tung, D., Hsieh, C., Jin, Z., Shi, L., Marthaler, D., Bertozzi, A., and Murray, R. M., “Virtual Attractive-Repulsive Potentials for Cooperative Control of Second Order Dynamic Vehicles on the Caltech MVWT,” Proceedings of the American Control Conference (ACC), 2005, pp. 10841089.

13

Bennet, D. J., MacInnes, C. R., Suzuki, M., and Uchiyama, K., “Autonomous Three-Dimensional Formation Flight for a Swarm of Unmanned Aerial Vehicles,” Journal of Guidance, Control, and Dynamics, Vol. 34, No. 6, 2011, pp. 18991908.

14

Gowtham, G. and Kumar, K. S., “Simulation of Multi UAV Flight Formation,” The 24th Digital Avionics Systems Conference, DASC, 2005, pp. 6.

15

Crowther, B., “Flocking of Autonomous Unmanned Air Vehicles,” Aeronautical Journal, Vol. 107, No. 1068, 2003.

16

Wolf, M. T. and Burdick, J. W., “Artificial Potential Functions for Highway Driving with Collision Avoidance,” IEEE International Conference on Robotics and Automation, 2008, pp. 37313736.

17

Shim, D. H., Kim, H. J., and Sastry, S., “Decentralized Nonlinear Model Predictive Control of Multiple Flying Robots,” Proceedings of the 42nd IEEE Conference on Decision and Control, Vol. 4, 2003, pp. 36213626.

18

Shin, J. and Kim, H. J., “Nonlinear Model Predictive Formation Flight,” IEEE Transactions on Systems, Man and Cybernetics, Part A: Systems and Humans, Vol. 39, No. 5, 2009, pp. 11161125.

19

Lam, T. M., Boschloo, H. W., Mulder, M., and Van Paassen, M. M., “Artificial Force Field for Haptic Feedback in UAV Teleoperation,” IEEE Transactions on Systems, Man and Cybernetics, Part A: Systems and Humans, Vol. 39, No. 6, 2009, pp. 13161330.

20

Clark, J. B. and Jacques, D. R., “Flight Test Results for UAVs Using Boid Guidance Algorithms,” Procedia Computer Science, Vol. 8, 2012, pp. 232238.

21

Bamberger, J. R. J., Watson, D. P., Scheidt, D. H., and Moore, K. L., “Flight Demonstrations of Unmanned Aerial Vehicle Swarming Concepts,” Johns Hopkins APL Technical Digest, Vol. 27, No. 1, 2006, pp. 4155.

22

Stastny, T. J., Garcia, G. A., and Keshmiri, S. S., “Collision and Obstacle Avoidance in Unmanned Aerial Systems Using Morphing Potential Field Navigation and Nonlinear Model Predictive Control,” ASME Journal of Dynamic Systems, Measurement, and Control, Vol. 137, No. 1, 2015.


23

University of Kansas, “Center for Remote Sensing of Ice Sheets,” https://www.cresis.ku.edu/, 2014, [Online; accessed 11-Apr-2014].

24

Garcia, G. and Keshmiri, S., “Online Artificial Neural Network Model-Based Nonlinear Model Predictive Controller for the Meridian UAS,” International Journal of Robust and Nonlinear Control, 2013.

25

Royer, D., Sweeten, B. C., and Jones, V., “Modeling and Sensitivity Analysis of the Meridian Unmanned Aircraft,” AIAA Infotech@Aerospace, 2010.

26

Garcia, G. and Keshmiri, S., “Nonlinear Model Predictive Controller for Navigation, Guidance and Control of a Fixed-Wing UAV,” AIAA Guidance, Navigation, and Control Conference, 2011.

27

Garcia, G. A., Decentralized Robust Nonlinear Model Predictive Controller for Unmanned Aerial Systems, Ph.d. thesis, University of Kansas, 2013.

28

Garcia, G., Keshmiri, S., and Colgren, R., “Advanced H-Infinity Trainer Autopilot,” AIAA Modeling and Simulation Technologies Conference, 2010.

29

Park, S., Deyst, J., and How, J., “A New Nonlinear Guidance Logic for Trajectory Tracking,” AIAA Guidance, Navigation, and Control Conference and Exhibit, 2004.

30

Park, S., Deyst, J., and How, J. P., “Performance and Lyapunov Stability of a Nonlinear Path-Following Guidance Method,” AIAA Journal on Guidance, Control, and Dynamics, Vol. 30, 2007.

31

Kwakernaak, H., “H2-OptimizationTheory and Applications to Robust Control Theory,” Annual Reviews in Control, Vol. 26, No. 1, 2002, pp. 69–79.

32

Kwakernaak, H., “Robust Control and H-infinity Optimization Tutorial Paper,” Automatica, Vol. 29, No. 2, 1993, pp. 255–273.

33

Ferreira, H., Rocha, P. H., and Sales, R. M., “Nonlinear H-infinity Control and the Hamilton-Jacobi-Isaac Equation,” 17th World Congress, The international Federation of Automatic Control, Seoul, Korea, 2009, pp. 188–193.

34

Boggs, P. T. and Tolle, J. W., “Sequential Quadratic Programming,” Acta Numerica, Vol. 4, 1995, pp. 151.

35

Chang, D. E., Shadden, S. C., Marsden, J. E., and Olfati-Saber, R., “Collision Avoidance for Multiple Agent Systems,” Proceedings of the 42nd IEEE Conference on Decision and Control, Vol. 1, 2003, pp. 539543.

36

Nelson, D. R., Barber, D. B., McLain, T. W., and Beard, R. W., “Vector Field Path Following for Miniature Air Vehicles,” IEEE Transactions on Robotics, Vol. 23, No. 3, 2007, pp. 519529.

37

Stengel, R. F., Flight Dynamics, Princeton University Press, Princeton, NJ, 2004.

38

The Hindu, “UAV Crashes Into Coconut Grove,” http://www.thehindu.com/news/national/tamil-nadu/ uav-crashes-into-coconut-grove/article5457371.ece, 2013, [Online; accessed 11-Apr-2014].

39

Jane, F. T., Jane’s All the World’s Aircraft, Jane’s Yearbooks, 1976.
