Fátima Oliva-Palomo

Anand Sanchez-Orta

Hussain Alazki
email: halazki@pampano.unacar.mx

Pedro Castillo
email: castillo@hds.utc.fr

Aldo-Jonathan Muñoz-Vázquez
email: aldo.munoz-vazquez@tamu.edu

Robust Global Observer Position-Yaw Control based on Ellipsoid Method for Quadrotors

Keywords: Position Control, Unmanned Aerial Vehicle, Real-time Validation, Quaternion Control

An observed-based control strategy for position-yaw tracking of quadrotors is proposed. From a virtual controller for the position dynamics, attitude and angular velocity references are estimated with the aim of allowing the position-control components to enforce the tracking of a position reference. In addition, no measurement of the linear velocity is assumed, and the measurement of the position signal is noisy, thus, a dynamic observer for the position dynamics is employed. The ellipsoid method is considered to estimate the feedback controller and observer gains, based on an optimization procedure, in order to produce the ultimately boundedness of attitude and position tracking errors inside a compact convex vicinity of the origin, whose Lebesgue measure is minimized via LMIs. The advantages of the proposed scheme are illustrated through numerical simulations that include aerodynamic disturbances and noisy measurements. Finally, an experimental assessment is carried out to illustrate the feasibility and performance of the proposed controller in a real-time application.

Introduction

The quadrotor is a highly maneuverable vehicle that is capable of performing aggressive trajectories. This characteristics make it attractive for many different applications. Unfortunately, the fast dynamics and disturbances make it prone to oscillations, therefore, this small vehicle requires a fast and robust controller. The controller of the full actuated attitude dynamics has to guarantee the stability of the entire system, whilst being fast enough, considering that the position dynamics depend on the attitude [START_REF] Oliva-Palomo | A fractional nonlinear pi-structure control for robust attitude tracking of quadrotors[END_REF][START_REF] Izaguirre-Espinosa | Attitude control of quadrotors based on fractional sliding modes: theory and experiments[END_REF]. For this purpose, several contributions have been reported, some of them using the Euler angles as attitude representation because are more intuitive, but such representation has singularities [START_REF] Bouabdallah | Full control of a quadrotor[END_REF][START_REF] Ferrin | Differential flatness based control of a rotorcraft for aggressive maneuvers[END_REF][START_REF] Cabecinhas | A globally stabilizing path following controller for rotorcraft with wind disturbance rejection[END_REF][START_REF] Mellinger | Trajectory generation and control for precise aggressive maneuvers with quadrotors[END_REF][START_REF] Shao | Robust back-stepping output feedback trajectory tracking for quadrotors via extended state observer and sigmoid tracking differentiator[END_REF][START_REF] Cai | Quadrotor trajectory tracking and obstacle avoidance by chaotic grey wolf optimization-based active disturbance rejection control[END_REF][START_REF] Zhao | Attitude control for quadrotors subjected to wind disturbances via active disturbance rejection control and integral sliding mode control[END_REF][START_REF] Wang | Disturbance observer-based adaptive fault-tolerant control for aquadrotor helicopter subject to parametric uncertainties andexternal disturbances[END_REF][START_REF] Cai | Quadrotor trajectory tracking and obstacle avoidance by chaotic grey wolf optimization-based active disturbance rejection control[END_REF], which can restrict the attitude space [START_REF] Ferrin | Differential flatness based control of a rotorcraft for aggressive maneuvers[END_REF], discarding the tracking of highly aggressive maneuvers. Additionally, the use of quaternion algebra has been employed for using Euler parameters for attitude representation, avoiding singularities to produce a global controller [START_REF] Tayebi | Attitude stabilization of a vtol quadrotor aircraft[END_REF], even for aggressive maneuvers [START_REF] Wang | Switching control of attitude tracking on a quadrotor uav for large-angle rotational maneuvers[END_REF][START_REF] El-Badawy | Quadcopter aggressive maneuvers along singular configurations: An energy-quaternion based approach[END_REF][START_REF] Hua | Auto-tuning nonlinear pidtype controller for rotorcraft-based aggressive transportation[END_REF]. Taking in to account aerodynamic disturbances, there some control approaches consider sliding mode [START_REF] Zhao | Attitude control for quadrotors subjected to wind disturbances via active disturbance rejection control and integral sliding mode control[END_REF][START_REF] Xiong | Position and attitude tracking control for a quadrotor uav[END_REF][START_REF] Zhao | Attitude control for quadrotors subjected to wind disturbances via active disturbance rejection control and integral sliding mode control[END_REF] or adaptive [START_REF] Roberts | Adaptive position tracking of vtol uavs[END_REF] controllers, which can induce high frequencies in the controller signals, increasing the oscillation of the system. Recently, a neuroadaptive integral robust controller for visual quadrotor tracking is proposed in [START_REF] Shao | Neuroadaptive integral robust control of visual quadrotor for tracking a moving object[END_REF]. Asymptotic stability of the closed-loop system is guaranteed using Barbalat's lemma, however, conservative bounds of uncertainties are assumed to compute the controller gains and only simulations are provided.

Additionally, given the underactuation of the vehicle, the attitude has to be slaved to the position to allow tracking. How to modify the attitude dynamics makes the position control a nontrivial problem, for which, different strategies have been developed to solve underactuation. One approach is based on a position controller that is accounted for a disturbance to the attitude dynamics [START_REF] Sanchez | Continuous reactive-based positionattitude control of quadrotors[END_REF], which causes the vehicle to tilt to the desired direction in order to track the position reference, however, it has the drawback that the attitude controller is capable to compensate at least a part of the position controller, leading to a large position error. Another method consists in computing the desired attitude reference from the desired position trajectory [START_REF] Mellinger | Trajectory generation and control for precise aggressive maneuvers with quadrotors[END_REF], but in that case, all dynamic parameters have to be known in order to compute an accurate trajectory for the attitude dynamics, and even so, any of the disturbance causes a poor performance. A further more strategy is to compute the desired attitude from the desired acceleration that can be in turns computed from a virtual controller [START_REF] Roberts | Adaptive position tracking of vtol uavs[END_REF], nevertheless, in this method the rotation about the z axis is restrained at 0 in order to allow the tracking, leavening one degree of freedom unused. In addition, a virtual controller (desired force on the vehicle) can be used in order to compute the desired references for the attitude [START_REF] Cabecinhas | A globally stabilizing path following controller for rotorcraft with wind disturbance rejection[END_REF], leaving only the design of a controller that does not induce high frequencies in the system dynamics.

By virtue of the aforementioned discussion, it is evident that the quadrotor needs a robust controller to compensate for dynamic disturbances, whilst avoiding high frequencies and kinematic singularities. Along this direction, the attractive ellipsoid method allows one to design a continuous controller that is robust against a broad class of disturbances and uncertainties, assuring practical stability [START_REF] Poznyak | Attractive Ellipsoids in Robust Control[END_REF]. It is assumed that there is not complete information of the plant, such as the external disturbances, but at least, an upper bound is known. In this methodology, the feedback gains are calculated from an optimization problem, where Linear Matrix inequalities (LMIs) appear, and by solving them numerically, the region at which the tracking errors converge is reduced. The solution of convex inequalities has allowed the development of important contributions to the control of aerial vehicles, such as the computation of observer gains in fault estimation schemes [START_REF] Zhu | Robust fault estimation for a 3-dof helicopter considering actuator saturation[END_REF].

The contribution is then, an observer based on position feedback by using the Ellipsoid Method. The virtual controller is defined in the position dynamics, which represents the desired force in the vehicle for position tracking. From the virtual controller, the attitude references are estimated, which are represented by unit quaternions in order to avoid possible singularities. Additionally, a desired rotation about the z axis is performed while tracking the position. The controller is designed based on the ellipsoid method, assuring a small region at which the errors converge, depending on the solution of LMIs to estimate the optimal gains. The validity of proposed scheme is shown throughout a simulation study, subject to a time-varying wind gust and noisy measurements, and experimental results, which confirm the robustness of the proposed controller to wide range of uncertainties and disturbances. where ξ = [ξ x , ξ y , ξ z ] T ∈ I is the position of the vehicle, e z = [0, 0, 1] T , m ∈ R is the mass, g ∈ R represents the gravitational acceleration magnitude, f d (t) ∈ I denotes the bounded disturbance forces, q = [q 0 , q 1 , q 2 , q 3 ] T ∈ S 3 is the unit quaternion that represents the rotation of the vehicle, whose scalar and vector parts are q 0 ∈ R and q = [q 1 , q 2 , q 3 ] T , respectively, this leads to

Background and Problem statement

m ξ = -T R(q)e z + mge z + f d (t) (1)
q = 1 2 q ⊗ Ω (2) J Ω = -Ω × JΩ + τ + τ d (t) (3)
q = [ q 0 q ] = [ cos(θ) λ sin(θ). ] (4)
The quaternion product ⊗ is given by [START_REF] Kuipers | Quaternions and rotation sequences: a primer with applications to orbits, aerospace and virtual reality[END_REF] p

⊗ q = [ p 0 -p T p I 3×3 p 0 + [p×] ] [ q 0 q ] (5)
and the conjugate of a quaternion is

q * = [q 0 , -q 1 , -q 2 , -q 3 ] T (6)
Note that, in the case unit quaternions, q * = q -1 . Ω = [0, Ω T ] T denotes a pure quaternion with scalar part of zero and vector part Ω, where

Ω = [Ω 1 , Ω 2 , Ω 3 ]
T ∈ B defines the angular velocity of the airframe. In addition the rotation matrix R(q) is expresed in terms of the unit quaternion, as is equivalent to

v i = q ⊗ v b ⊗ q * = [ 0 R(q)v b ] , (7)
where v i and v b are pure quaternions,

v i = [v i1 , v i2 , v i3 ] T ∈ I, v b = [v b1 , v b2 , v b3 ]
T ∈ B, and R(q) is given by

R(q) =   1 -2q 2 2 -2q 2 3 2q1q2 -2q0q3 2q1q3 + 2q0q2 2q1q2 + 2q0q3 1 -2q 2 1 -2q 2 3 2q2q3 -2q0q1 2q1q3 -2q0q2 2q2q3 + 2q0q1 1 -2q 2 1 -2q 2 2   (8)
J ∈ R 3×3 stands for the constant inertia matrix around the center of mass and is expressed in the body fixed frame B. The matrix [Ω×] ∈ R 3×3 denotes the skew-symmetric matrix of the vector

Ω [Ω×] =   0 -Ω 3 Ω 2 Ω 3 0 -Ω 1 -Ω 2 Ω 1 0   , (9)
the external bounded torque disturbance is τ d (t), and τ ∈ B is the control torque.

The Class of quasi-Lipschitz Functions

The quadrotor model is nonlinear, then the following definition will be useful in the stability analysis.

Definition 1 [START_REF] Alazki | Robust output stabilization for a class of nonlinear uncertain stochastic systems under multiplicative and additive noises: The attractive ellipsoid method[END_REF]: A vector function g : R n → R k is said to be from the class C(C, δ 0 , δ 1 ) of quasi-Lipschitz functions if there exists a matrix C ∈ R k×n and nonnegative constants δ 0 and δ 1 such that for any x ∈ R n the following inequality holds:

∥g(x) -Cx∥ 2 ≤ δ 0 + δ 1 ∥x∥ 2 (10)
This implies that the growth rates of g(x) as ∥x∥ → ∞ are not faster than linear, [START_REF] Poznyak | Attractive Ellipsoids in Robust Control[END_REF]. Notice that if C = 0, δ 0 = 0 and g(0) = 0 the inequality (10) characterizes the Lipschitz continuity property of the function g(x) with the Lipschitz constant L = √ δ 1 .

Attractive Ellipsoids

In order to determine the convergence region and to apply the ellipsoid method, the following definition is needed.

Definition 2 [START_REF] Poznyak | Attractive Ellipsoids in Robust Control[END_REF]: An ellipsoid, represented by ε ⊂ R n with center in x c is given by:

ε = { x ∈ R n | (x -x c ) T P (x -x c ) ≤ 1 } (11)
is said to be attractive for the trajectories x :

lim sup t→∞ (x -x c ) T P (x -x c ) ≤ 1
where the ellipsoidal matrix P is a symmetric positive definite matrix 0 < P = P T ∈ R n×n ,

Problem formulation

The problem is to design a closed-loop controller allowing robust position tracking for system (3), assuring that the error trajectories are confined to an invariant attractive set. The controller is continuous and independent of model parameters, but its design consider the system model, ensuring robustness to parametric uncertainties, aerodynamical effects and exogenous disturbances. In addition, considering underactuation, the desired attitude consigns that guarantee position tracking are also needed.

Mapping based Position Control for a Quadrotor

The quadrotor is an underactuated vehicle, then, in order to reach a desired position, the vehicle needs to rotate to change the direction of the thrust to generate force components in a desired direction. In consequence, the vehicle gives up some degrees of freedom in the attitude to gain degrees of freedom in position. In this sense, a virtual control u d is defined, where the components of this vector corresponds to the input forces needed to track the desired position. Then, the virtual controller u d corresponds to the desired thrust T d in a desired direction, which is defined by rotation represented by the unit quaternion q d , see Fig. 2, consequently the virtual control can be defined as

u d = -T d R(q d )e z ( 12
)
Therefore, using (4), ( 8), [START_REF] Tayebi | Attitude stabilization of a vtol quadrotor aircraft[END_REF], and considering that the magnitude of the virtual controller is the desired thrust T d , then the vector that define the desired direction ûd =

u d ∥u d ∥ = -R(q d
)e z and is given by

ûd =   -2q 1d q 3d -2q 0d q 2d -2q 2d q 3d + 2q 0d q 1d -1 + 2q 2 1d + 2q 2 2d   =   ûd1 ûd2 ûd3   (13)
Now using the unit quaternion restriction q 2 0d + q 2 1d + q 2 2d + q 2 3d = 1 and from [START_REF] Wang | Switching control of attitude tracking on a quadrotor uav for large-angle rotational maneuvers[END_REF], the scalar part of the quaternion, expressed in terms of q 3d and ûd3 , equals to

q 0d = 1 2 √ -2û d3 + 2 -4q 2 3d (14)
then, from ( 13) and ( 14), the solution for q 1d and q 2d is given by

[ q 1d q 2d ] = 1 2(q 2 0d +q 2 3d ) [ -q 3d q 0d -q 0d -q 3d ] [ ûd1 ûd2 ] =   q 3d ûd 3 -1 -1 2 √ -2 ûd 3 +2-4 q 3d 2 ûd 3 -1 1 2 √ -2 ûd 3 +2-4 q 3d 2 ûd 3 -1 q 3d ûd 3 -1   [ ûd1 ûd2 ] (15)
The scalar part q 0d and the two components of the vector part q 1d and q 2d are in terms of q 3d and the direction defined by the vector ûd . Given that there is no restriction of which value should be q 3d , let define the final rotation in two steps. First, assuming that the desired rotation axis is in the xy plane, ie q 3d = 0, the desired quaternion equals to

q dxy =       1 2 √ -2û d3 + 2 ûd 2 √ -2û d 3 +2 -û d 1 √ -2û d 3 +2 0       (16)
After the vehicle is rotated to generate the desired force components, we can perform another rotation of a desired ψ d angle about the z axis on the body of the quadrotor, as can be seen in Fig. 3, this rotation is defined by

q dz =       cos ( ψ d 2 ) 0 0 sin ( ψ d 2 )       , ( 17
)
This last rotation does not change the direction of the desired force, which implies that the tracking of the position it is not modified. Therefore, the desired quaternion for the attitude is given by the consecutive rotations of q dxy and q dz as follows

q d = q dxy ⊗ q dz =           1 2 √ -2û d3 + 2 cos ( ψ d 2 ) -û d 1 sin ( ψ d 2 ) +û d 2 cos ( ψ d 2 ) √ -2û d 3 +2 -û d 1 cos ( ψ d 2 ) -û d 2 sin ( ψ d 2 ) √ -2û d 3 +2 1 2 √ -2û d3 + 2 sin ( ψ d 2 )           (18)
Note that in (15) q 1d and q 2d are not well defined if ûd3 = 1, that means that the direction of the desired force is e z . It is assumed that this scenario is highly improbable due to the gravitational force compensation in the control law, which means that the component ûd3 is always pointing upwards or ûd3 < 0. Also, note that the desired quaternion in ( 18) is equivalent to the equations ( 14) and ( 15) using

q 3d = 1 2 √ -2û d3 + 2 sin ( ψ d 2
) . Now, to obtain the desired angular velocity, from (2), there exists a correspondence between the desired quaternion and the desired angular velocity, and its given by

qd = 1 2 q d ⊗ Ω d (19)
then, the desired angular velocity is equal to

Ω d = 2 q * d ⊗ qd (20)
with the time derivative of (18), ûT d ûd = 1, its derivative ûT d ud = 0, we obtain that ( 20) is a pure quaternion and the angular velocity equals to

Ω d =      -sin(ψ d ) ud 1 +cos(ψ d ) ud 2 + ud 3 ( sin(ψ d )û d 1 -cos(ψ d )û d 2 ) ûd 3 -1 -cos(ψ d ) ud 1 -sin(ψ d ) ud 2 + ud 3 ( cos(ψ d )û d 1 +sin(ψ d )û d 2 ) ûd 3 -1 ψd + ûd 1 ud 2 -ûd 2 ud 1 ûd 3 -1      (21)
Note that for the computation of the desired angular velocity the time derivative of the virtual position control ud is needed, therefore this controller needs to be at least one time differentiable.

Attitude Control Design

With the purpose of allow to generate the desired input forces defined by the virtual control, an attitude controller has to assure that the systems tracks the desired references computed from the virtual controller. The control has to compensate the disturbances and aerodynamical effects present in the attitude dynamics and the design has to consider the nonlinearities of the model.

In order to design the attitude controller, first let define the attitude and angular velocity errors.

Attitude Error Equations

The rotation error quaternion is defined by

q e ≜ q * d ⊗ q (22)
which represent the rotation difference between the actual and the desired attitude. Now to obtain the relation of the quaternion error and a angular velocity error, we obtain the derivative of ( 22), using ( 2), ( 6), ( 5), ( 19) and ( 22), is equal to

qe = 1 2 q e ⊗ Ω e (23)
where the error of angular velocity is

Ω e = Ω -R T (q e )Ω d ( 24
)
The attitude dynamics has 3 degrees of freedom, the same that the unit quaternion has (4 different parameters and one restriction ∥q e ∥ = 1). Assuring the 3 components of the vector part q e → 0 implies that the scalar part q 0e → ±1, then we only need to assure q e → 0 in order to have q → q d . Also, considering small errors, with q 0e ≈ 1, it can be assumed that the scalar part of the quaternion error doesn't change sign, if the initial condition is q 0e (t 0 ) ≈ 1 that means that q e0 ≥ 0 and with the unit quaternion restriction, then

q e0 = √ 1 -q 2 e1 -q 2 e2 -q 2 e3 ( 25
)
with this the differential equation of the vectorial part of the quaternion is given by

qe = Q(q e )Ω e , Q(q e ) = 1 2
[ (1 -q T e q e ) 1/2 I 3×3 + [q e ×] ] (26
) Now the state space can be defined by error variables that must converge to a small invariant set as follows

x a = [ q e Ω e ] , ẋa = B a τ + ζ a (x a , Ω, Ω d , Ωd , t) (27)
where ζ a (x a , Ω d , Ωd , t) corresponds to the disturbances and non-linearities of the attitude dynamics

ζ(x a , Ω d , Ωd , t) = [ Q(qe)Ωe J -1 (-Ω×JΩ+τ d )-R T e Ωd -R T e Ω d ×Ωe
]

and 27), the quasi-linear model can be represented as

B a = [0 3×3 J -1 T ] T From (
ẋa = A a x a + B a τ + φ a (x a , Ω d , Ωd , t) (28)
where

φ(x a , Ω, Ω d , Ωd , t) = ζ(x a , Ω, Ω d , Ωd , t)-A a x a , A a = [ -aI 3×3 1 2 I 3×3 0 3×3 0 3×3 ] , with a > 0 ∈ R is a positive scalar.
The form of A a comes from the linearization around the origin of Q(q) and the constant a for and additional parameter for the solution of the ellipsoid method.

Assumptions of the attitude dynamics Assumption 1:

The desired Ω d and Ωd are defined by the user, and are bounded by

∥Ω d ∥ 2 ≤ Ω + d ∥ Ωd ∥ 2 ≤ Ω+ d ( 29
)
Assumption 2: The norm of the disturbances vector τ d has an upper bound, then there exists a positive constant d + such that

∥τ d ∥ 2 < d + (30)
Assumption 3: The nonlinear term in (3) can be bounded using two constants c 0 ≥ 0 and c 1 > 0 such that

∥J -1 (-Ω × JΩ)∥ 2 ≤ c 0 + c 1 ∥Ω∥ 2 ∀ Ω ∈ R 3 (31)

Bound of the nonlinearities and disturbances

Now with the purpose of proceed to the stability analysis, the bound of φ is needed. With the assumptions 1-3, the norm of φ(x, Ω d , Ωd , t) is calculated as

φ(x a , Ω d , Ωd , t) T φ(x, Ω d , Ωd , t) = φ T 1 φ 1 + φ T 2 φ 2 (32) φ T 1 φ 1 ≤ 3 2 ∥Ω e ∥ 2 + 3a 2 ∥q e ∥ 2 φ T 2 φ 2 ≤ c 3 + 4c 1 ∥Ω∥ 2 + c 4 ∥Ω e ∥ 2 where c 3 = 4c 0 + 4J + d + + 4 Ω+ d y c 4 = 4Ω + d . Then, with b 0 = c 3 + 3a 2 , b 1 = 4c 1 , b 2 = 3 2 + c 4 , the norm of ∥φ∥ 2 is given by ∥φ∥ 2 = b 0 + b 1 ∥Ω∥ 2 + b 2 ∥Ω e ∥ 2 . ( 33
)
The bound of the nonlinearities and disturbances depends on the angular velocity of the vehicle. In order to express (33) only in terms of the error and desired angular velocity, an extended vector with the system state and the desired error velocity is defined

xa =   q e Ω e R T e Ω d   = [ x a R T e Ω d ] , (34)
using [START_REF] Alazki | Robust output stabilization for a class of nonlinear uncertain stochastic systems under multiplicative and additive noises: The attractive ellipsoid method[END_REF] we can calculate the angular velocity of the vehicle Ω = Ω e + R T e Ω d , that can be expressed as:

Ω = G 1 xa , G 1 = [ 0 3×3 I 3×3 I 3×3 ] ,
the norm of the angular velocity ∥Ω∥ 2 can be expressed in terms of the extended error vector (34) as

∥Ω∥ 2 = xT a G T 1 G 1 xa = xT a G 2 xa , G 2 =   0 3×3 0 3×3 0 3×3 0 3×3 I 3×3 I 3×3 0 3×3 I 3×3 I 3×3 ,  
then the norm of ∥φ∥ 2 from (33) can be expressed by

∥φ∥ 2 = b 0 + xT a Gx a , G = [ 03×3 03×3 03×3 03×3 (b1+b2)I3×3 b1I3×3 03×3 b1I3×3 b1I3×3 ] . ( 35
)
and the state x a can be defined as

x a = H xa , H = [ I 6×6 0 6×3 ] (36)

Attitude Practical Stability

Consider the following control law

τ = Kx a (37) Theorem 1. Suppose that 1.
-The assumptions 1-3 are fulfilled 2.-For the system (27) with (37), for some matrices K and P a = P T a > 0 (solution of the optimization problem (46), when it exists), also some positive constants α a , ε 1 , ε 2 that satisfy the inequality W a (P a , K, α a , ε 1 , ε 2 ) < 0 where

W a =   G a ε 1 G 12 P a ε 1 G T 12 (ε 1 b 1 -ε 2 )I 3×3 0 3×6 P a 0 6×3 -ε 1 I 6×6   (38) G a = (A a + B a K + αa 2 I 6×6 ) T P a + P a (A a + B a K + αa 2 I 6×6 ) + ε 1 G 11 , G 12 = [ 0 3×3 b 1 I 3×3 ] , G 11 = [ 0 3×3 0 3×3 0 3×3 (b 1 + b 2 )I 3×3 ] ,
there exists an energy function V a (P a , x a ) that has a differential inequality defined as

Va ≤ -α a V a (P a , x a ) + β a (39)
where, β a is a positive constant.

Then, there exists an invariant ellipsoid for the closedloop system [START_REF]Optitrack for robotics[END_REF] and (37).

Proof. In order to analyze the behavior for the system, the energy function is proposed:

V a = x T a P a x a , (40)
Now, adding and subtracting ε 1 ∥φ∥ 2 and αV in the derivative of (40) follows that 35) and (36) in (41) it follows

Va = (A a x a + B a Kx a + φ) T P a x a + x T a P a (A a x a + B a Kx a + φ) ± ε 1 ∥φ∥ 2 ± α a V a (41) Using (
Va ≤ [ xa R T e Ω d φ ] T [ G3 H T Pa PaH -ε1I6×6 ] [ xa R T e Ω d φ
]

+ε 1 b 0 -α a V a ± ε 2 ∥R T e Ω d ∥ 2 ≤ [ xa R T e Ω d φ ] T [ Ga ε1G12 Pa ε1G T 12 (ε1b1-ε2)I3×3 03×6 Pa 06×3 -ε1I6×6 ] [ x R T e Ω d φ
]

+ε 1 b 0 + ε 2 Ω + d -α a V a (42) where G 3 = H T ((A a + B a K + αa 2 I 6×6 ) T P a + P a (A a + B a K + αa 2 I 6×6 ))H + ε 1 G, Then, the equation (41) can be expressed as Va ≤ [ xa φ ] T W a [ xa φ ] -α a V a + β a ( 43
)
where

β a = ε 1 b 0 + ε 2 Ω + d . Assuring that the matrix W a ≤ 0, the following in- equality is guaranteed Va ≤ -α a V a + β a ( 44
)
and assures that the stability region is defined by the ellipsoid given by

E a = { x a ∈ R n | x T a P a x a ≤ β a α a } (45)
Now, the minimization of the stability region can be found using the following optimization problem tr

{ βa αa P -1 a } → min αa,βa,ε1,ε2,K,Pa
subject to the restrictions

α a > 0, ε 1 > 0, ε 2 > 0, 0 < P a , W a = W a (K, P a , α a , ε 1 , ε 2 ) ≤ 0 (46)
which can be solved by using the MATLAB toolboxes Se-DuMi and YALMIP through the Interior Point Method [START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF].

Position Control Design

Now with an attitude control, we proceed to the design of the control law for the position dynamics, and how should be computed the virtual control in order to assure the tracking of the desired position.

Observer equations

A way of estimate the velocity of the vehicle is needed because only position measurement is available. Using a dirty derivative, although easy to implement, increases the noise present in the signals of position. In order to reduce the noise an observer is proposed to estimate the velocity of the vehicle.

A space state of the position is defined as

x ξ = [ x ξ1 x ξ2 ] = [ ξ ξ] (47)
with the position dynamics

ξ = 1 m (u o + f d (t)), u o = -T R(q)e z + mge z (48)
ẋξ and the output of the system are

ẋξ = A p x ξ + B p u o + B p f d (t), y = Cx ξ + η ( 49
)
where η ∈ R 3 represents the noise present in the measurements,

A p = [ 0 3×3 I 3×3 0 3×3 0 3×3 ] , B p = [ 0 3×3 1 m I 3×3 ] , C = [ I 3×3 0 3×3 ]
and the observer equation is given by

ẋξ = A p xξ + B p u o + L(y -C xξ ). (50)
where L ∈ R 6×3 .

Position Error Equations

In order to analyze how the virtual control affects the position stability, let introduce this term in the position dynamics. From (1) adding and substracting the virtual control u d defined in [START_REF] Tayebi | Attitude stabilization of a vtol quadrotor aircraft[END_REF], and T d R(q)e z , the acceleration is given by

m ξ = u d + f d (t) + ζ ξ (T, T d , q, q d ) + mge z ( 51
)
where

ζ ξ (T, T d , q, q d ) = (-T e R(q) + T d R(q d )(I 3×3 -R(q e )
))e z (52) Note that ζ ξ (T, T d , q, q d ) = 0 if T = T d and q = q d . This means that the term ζ ξ can be seen as a disturbance produced by the difference of the desired and the actual input forces in the position dynamics. Since the desired attitude gives the direction of the input force vector, then the norm of ζ ξ increases as the attitude error grows.

The desired position state is defined by Then, using (47), ( 54) and (51) the error dynamics is expressed as

x d = [ ξ d ξd ] (53
ẋe = [ ξ -ξd 1 m (u d + f d (t) + mge z + ζ ξ (T, T d , q, q d )) -ξd ]
(55) Given that the gravity, mass and desired acceleration are well known, in the virtual control can be compensated the gravitational force and the desired acceleration, as a result the virtual controller is given by

u d = u p -mge z + m ξd (56)
where u p is the feedback of the system signals. With (56), we can define a state space for the error dynamics as

ẋe = A p x e + B p u p + B p ζ p , ( 57
)
where

ζ p = [0 1×3 , (f d (t) + ζ ξ (T, T d , q, q d )) T ] T .
Defining the observer feedback controller as

u p = K p xe , (58)
from ( 49), ( 50) and ( 54), the equation of the observer error is defined by

ẋe = (A p -LC)x e -B p f d (t) + Lη. ( 59
)
also, using (54) it follows that xe = x e + xe , then from (55), the position error dynamics is expressed as

ẋe = (A p + B p K p )x e + B p K p xe + B p ζ p ( 60
)
Defining an extended error for the position error variables that we seek to stabilize x p = [x T e , xT e ] T , which time derivative is using (60) and ( 59)

ẋp = Ãx p + B ζp (61)
where

ζT p = [ζ ξ T f d (t) T η T ] T , Ã = [ Ap + BpKp BpKp 06×6 Ap -LC ] , B = [ Bp Bp 06×3 06×3 -Bp L ]

Assumptions of the Position Dynamics Assumption 4:

The rotor dynamics are fast, meaning that the time delay between the instant that desired thrust (T d ) is commanded and the instant that the actuators actually give that force is considered small enough to be neglected, and it can be assumed that the desired thrust is equal to the system thrust T d = T .

Assumption 5: Part of the disturbances in the acceleration error in (55) are caused by the attitude error represented in the term ζ ξ (T d , q, q d ). That means that the disturbances grows as the attitude error q e increases. Assuming that this term has an upper bound implies that the attitude error vector q e has an upper bound, that indicates that the attitude controller assures that the norm of the vectorial part of the error quaternion ∥q e ∥ ≤ q + e < 1, where q + e is a positive constant.

Assumption 6: The thrusters of the vehicle have an upper limit in the generated force, then there exists a positive constant T + > T .

Assumption 7: The external forces in the position dynamics f d (t) has an upper bound, then, there exists a positive constant f + d such that ∥f d (t)∥ 2 < f + d . Assumption 8: The position noise measurement has a bound, there exists a positive constant η + such that ∥η∥ 2 < η + , where η ∈ R 3 represents the noise present in the measurements.

The disturbance due to the attitude error can be bounded if the assumptions 4-6 are fulfilled as follows:

ζ T ξ ζ ξ = T 2 d ((I 3×3 -R(q e )))e z ) T (I 3×3 -R(q e )))e z ≤ 2T + 2 q + e 2 (62)

Practical Stability with Observer Feedback

Now we seek that the observer and position error converges to a small enough stable region around the origin.

Theorem 2. Suppose that:

1.-The assumptions 4-7 are fulfilled.

2.-For the system (59) and (60) with (58), for some matrices K p , L and P = P T > 0, also some positive constants α, ε that satisfy the inequality W p (P, K p , L, α, ε) < 0 where

W p = [ P Ã + ÃT P + αP P B BT P -εI 9×9 ] , ( 63
)
there exists an energy function V (P, x p ) with the following differential equation

V ≤ -αV + β ( 64
)
where β > 0.

As a result, there exists a stability region around the origin of x e and xe defined by the ellipsoid

E = { x p ∈ R n | x T p P x p ≤ β α } ( 65
)
Proof. In order to analyze the system, a energy function is proposed

V p = x T p P x p (66)
which derivative, using (62), adding and subtracting αV and ε∥ ζp ∥ 2 , using assumptions 7 and 8 for the bound of the disturbances and noise as ε∥ ζp

∥ 2 ≤ ε∥ζ ξ ∥ 2 +εf + d +εη + , is equal to V = x T p P ẋp + ẋT p P x p ± αV ± ε∥ ζp ∥ 2 ≤ [ x p ζp ] T W p [ x p ζp ] -αV + β ( 67
)
where β = εT + 2 q + e 2 + εf + d + εn + . As long as it is assured that W p ≤ 0, there exists a stability region around the origin of x e and xe defined by the ellipsoid

E = { x p ∈ R n | x T p P x p ≤ β α } (68)
and the solution to the optimization problem is tr

{ β α P -1 } → min α,β,ε,Kp,Pa subject to the restrictions α > 0, ε > 0, 0 < P, W p = W p (K p , P, α, ε, ) ≤ 0 (69)
As was mentioned in Assumption 5, the disturbances term ζ ξ in the position error dynamics depends on the attitude error q e , meaning that if the attitude controller doesn't assure the tracking of the desired references computed from the virtual controller, the position error will increase. In this sense, the stability of the position dynamics depends on the attitude stability, ie as q e → 0 implies ξ e → 0.

Note that for the derivative of the control ud ( ẋe , ...

ξ d ) = K p ẋe +m ... ξ d = K p ( ẋe -ẋd )+m ... ξ d
, the ẋe can be estimated using (50). Also the term ... ξ d is needed, which implies that the desired position trajectory needs to be at least 4 times differentiable.

Numerical Results

In order to demonstrate the performance of the controller a simulation study subject to a time-varying wind gust and noisy measurements was implemented. The vehicle accomplishes the tracking of a position trajectory while tracking a desired ψ d angle.

Simulator

The simulation was carried out using Matlab ® Simulink ® 9.2. The solver used was Dorman-Price (ode8) with a fixed step of 0.001s. The Fig. 4 shows how the signals were used to calculate the controllers and dynamics of the simulated system. In order to emulate a real scenario, a time-varying wind gust is simulated, where the wind velocity vector is sin 2 (0.1 t) [0.5403, 0.8415, 0] T m/s. Note that the magnitude of such a vector is v(t) = sin 2 (0.1 t) m/s. This wind gust generates state-dependent aerodynamical disturbances in the form of forces and moments that affect both position and orientation dynamics. The aerodynamic coefficients and the equations that describe these effects are [START_REF] Ramirez-Rodriguez | Robust backstepping control based on integral sliding modes for tracking of quadrotors[END_REF]:

f d = 1 2 ρ a v 2 r S R W T d   c L c Y c D   τ d = 1 2 ρ a v 2 r S   b c l c c m b c n  
where the coefficients c and b are the chord and the span of a blade, respectively, c D , c Y , c L , c l , c m and c n are aerodynamical non-dimensional coefficients of drag, sideforce lift, rolling, pitching and yawing moments, respectively.

v r = √ (v sin(α 0 )) 2 + (v cos(α 0 ) + v i ) 2
and ρ a is the air density, α 0 is the freestream angle of attack, v r is the resultant velocity in the propeller slipstream, v i is the induced velocity produced by the rotor,

v i ≈ √ Tm 2ρaAm
, and T m is the thrust of the rotor and A m is the area of the rotor disk. The orthogonal rotation matrix that relates the body fixed frame and the aerodynamic frame is

W d =   c α d c β d s β d s α d c β d -c α d s β d c β d -s α d s β d -s α d 0 c α d  
where c a = cos(a) and s a = sin(a), α d is the angle of attack, α d = arctan The initial conditions of the system are q = [1, 0, 0, 0] T , Ω = [0, 0, 0] T ξ = [0, 0, -1] T , ξ = [0, 0, 0] T and xξ = [0, 0, -1, 0, 0, 0] T . Additionally, white noise is added to the position, attitude quaternion and angular velocity using the white noise block of Matlab ® Simulink, which generates normally distributed random numbers. Each block has different seed and the parameter noise power of 0.1. The values provided by those blocks are then multiplied by 0.0001 in order to obtain values between ±0.003, which correspond approximately to a position error of 3mm, attitude error of 0.4 degrees and angular velocity error of 0.003 rad/s. It is important to point out that a motion capture system Optitrack produces position error less than 0.3mm and rotational error less than 0.05 degrees [START_REF]Optitrack for robotics[END_REF].

( v sin(α0) v cos(α0)+vi
The control parameters are given by: K = [diag [START_REF] Wang | Disturbance observer-based adaptive fault-tolerant control for aquadrotor helicopter subject to parametric uncertainties andexternal disturbances[END_REF][START_REF] Wang | Disturbance observer-based adaptive fault-tolerant control for aquadrotor helicopter subject to parametric uncertainties andexternal disturbances[END_REF][START_REF] Wang | Disturbance observer-based adaptive fault-tolerant control for aquadrotor helicopter subject to parametric uncertainties andexternal disturbances[END_REF], diag(10, 10, 10)], (40,40,40), -diag(80, 80, 80)],

K p = [-diag

L = [I; 100I]

The solution of LMIs from the optimization problem in the attitude dynamics is expressed as: λ(P a ) i ∈ {0.017587, 0.45995}, ε 1 = 0.0041433, ε 2 = 156.0903 and α a = 0.002973; and for the position controller: λ(P ) i ∈ {0.0014673, 0.74522}, ε = 2.8245, α = 1.012.

Task

The smooth desired trajectory, assuring small initial errors for position and attitude, is given by

ξ d =          r (tanh (δt -4) + 1) 2 cos (f t) r (tanh (δt -4) + 1) 2 sin (f t) -z d          (70)
The function tanh is used to slowly increase the radius of the circle from zero to a desired value r, where its increment depends on δ. This assures small initial errors provided that the initial position of vehicle is in the center of the circle. The trajectory parameters used were δ = 2, r = 1, f = 2 and z d = 1, with ψ = -2 for the desired attitude. The desired trajectory is smooth enough to guarantee existence of ... ξ d , which is a requirement for the solution of underactuation.

In order to calculate the desired trajectories for the attitude controller, the time derivative of the virtual position controller, u d , is required. By using (50), it is computed as

ud = -K p ẋe + m ... ξ d = -K p ( ẋξ -ẋd ) + m ... ξ d
Note also that the term ... ξ d is needed to compute the desired angular velocity in [START_REF] Poznyak | Attractive Ellipsoids in Robust Control[END_REF]. Thus, ... ξ d has to be continuous, which implies that the desired position trajectory needs to be at least 4 times differentiable.

Results

The tracking performance of the system is shown in Fig. 5. In Fig. 6, the 3D trajectory demonstrates how the vehicle begins in the center of the circle and increases the radius of the tracked circle until it reaches target radius of 1. The system follows the desired trajectory with small errors, as can be appreciated in Fig. 7. The position controller is shown in Fig. 8. As seen from this figure, the element u d3 reminds always negative, assuring that there is always a desired attitude at any given time. Fig. 9 shows the successful tracking of the desired attitude, computed from the position controller. However, the noise present in the position control results in a noisy desired attitude trajectory, specially for the desired angular velocity. Even with this noise, the attitude is tracked with small errors, see Fig. 10, and consequently ensures position tracking at the same time. The attitude controller signals are presented in Fig. 11, which has a strong component of noise due to the noise present in the measured and desired signals. Finally, Fig. 12 shows the moment and force aerodynamical disturbances respectively that affect the vehicle.

Experimental Results

The excellent performance of the control strategy is demonstrated when the aerial vehicle executes a similar trajectory to the one used in the numerical results, (70), a circle trajectory while tracking a desired angle ψ d , using the equations ( 18) and [START_REF] Poznyak | Attractive Ellipsoids in Robust Control[END_REF]. The platform used is an AR Drone 2.0 Edition of Parrot shown in Fig. 13. The platform uses a board with an ARM Cortex A8 1GHz processor with 1GB of RAM and the control algorithm is programmed in C++ code using Codeblocks 13 in a Linux Ubuntu 14.04 environment. The embedded system update the IMU measurements at 200Hz (5ms), as well as control outputs to the rotors, interruptions, data fusion and communication to the ground station. The measurements of the position of the vehicle is performed using an Optitrack system with 24 cameras, that send information to the vehicle up to 100Hz (10ms). The communication to the ground station allows to monitor the system measurements and change parameters in the controller and desired position, while all the data is recorded on board.

The parameters of the desired trajectory (70) are δ = 0.7, r = 1, f = 2, z d = 1.5 and ψd with initial condition ψ d (t 0 ) = 0. Using an approximation of the inertia matrix J = diag(0.002237568, 0.002985236, 0.00480374)

and mass m = 0.43kg for the optimization problem, the solution of the LMIs are for the attitude dynamics : λ(P a ) i ∈ {7.5990e -05, 4.4614e -01}, ε 1 = 0.055195, ε 2 = 2.0784 and α a = 0.00208; and for the position controller: λ(P ) i ∈ {1.6117e -04, 3.2167e + 00}, ε = 1.8468, α = 2.99.

The result of the tracking of the desired position trajectory are shown in Fig. 14. The comparison of the measured and desired position is shown in Fig. 14(a), also, the comparison of the observed and desired velocity is displayed in Fig. 14(b), where it can be easily seen the outstanding tracking of the desired trajectory. Usually, given the underactuation of the quadrotor, the tracking of a circle seams that is out of phase [START_REF] Sanchez | Continuous reactive-based positionattitude control of quadrotors[END_REF], which in this case is not appreciable. It is better seen in Fig. 15 the error xe , which magnitude is never greater than 20cm of the tracking of a 2m diameter circle.

In Fig. 16 the comparative of the two different forms to estimate the velocity of the vehicle is presented, one is a dirty derivative and the other one is using the observer. It is perceptible that the dirty derivative amplifies the noise present in the position measurements, while the observer estimate a more regular signal. This is important given that the desired attitude is computed using the error signals of position and velocity. Using more regular signals of the position and velocity signify more uniform desired trajectories for the orientation, making it easier for the attitude controller to perform the tracking, implying a better behavior in the position. The virtual controller presented in Fig. 17, is used to estimate the desired quaternion. Given that the virtual controller contains noise, the desired quaternion and angular velocity also contains noise. The signals of the real quaternion components q 1 , q 3 and the desired q 1d , q 3d (estimated from the virtual controller) are shown in Fig. 18. The main rotation is in the z axis, as can be seen in Fig. 18(b), the rotation about this axis has the same frequency as the circle, meaning that when the quadrotor finish a circle, a full 360º rotation is also completed. Even with noise present in the desired quaternion, the tracking has good performance, with small error in the attitude, as can be seen in the error quaternion vector from Fig. almost the same value plus some terms that come from the virtual controller as [START_REF] Poznyak | Attractive Ellipsoids in Robust Control[END_REF]. The performance of the attitude controller is presented in Fig. 22. Fig. 23 shows the tracking in 3D space, where can be seen that as the diameter of the circle increases, the altitude error increments as well due to the increment of the inclination of the vehicle. Finally, Fig. 24 exhibits how the vehicle has to tilt in order to follow the desired trajectory, almost like looking to the center of the circle as following the desired rotation about the z axis. This trajectories has frequencies greater than usually used in quadrotors [START_REF] Roberts | Adaptive position tracking of vtol uavs[END_REF][START_REF] Cabecinhas | A globally stabilizing path following controller for rotorcraft with wind disturbance rejection[END_REF][START_REF] Cao | Inner-outer loop control for quadrotor uavs with input and state constraints[END_REF], also, for more aggressive trajectories the Euler angles may present singularities that jeopardize the integrity of the vehicle. The video of the experiments can be seen in https://www.youtube.com/watch?v=1ZWpgedxMos.

Implementation considerations

The design of the control is done without considering the latencies present in the vehicle, it is only assumed that the sampling rate is fast enough to be considered continuous. The vehicle has a sampling rate of 5ms for the attitude signals, which is small enough to have a response considered continuous. However solution of the optimization problem could lead to high gains that would induce higher frequencies than the maximum frequency response of the system leading to instability. Also, a more aggressive trajectory in position implies a more aggressive trajectory in orientation, then, the invariant set for the attitude system could grow, leading to a increase in the disturbance in position due to orientation error and therefore a greater error in position.

Given the small errors requirement for the small error initial condition, the vehicle needs to start near from the desired trajectory, or even in the desired trajectory. If the vehicle needs to reach a specific point, the desired path has to be designed in order that gently leads from the starting point(where the vehicle begins) to the end point(the desired position).

Conclusions

A continuous position control for a quadrotor using observer feedback was proposed. The ellipsoid method forced the error to a small enough region around the origin, which allowed the system to track the desired position, even for highly aggressive maneuvers. The solution of underactuation provided the attitude desired references, needed by the attitude controller for the purpose of tracking the position trajectory, through a unit quaternion in order to avoid the singularities. The only restriction on the proposed scheme is that the desired force never points downwards, which implies that there is no error in the x and y directions and the desired point in z is far beneath the ve- hicle, a highly improbable scenario. Finally, the proposed approach allows the tracking in 4 Dof, in this case the 3 Dof of position and the yaw angle.

2. 1 .

 1 Quadrotor Dynamic Model Let I = {e x , e y , e z } be the earth fixed frame, and B = { e b x , e b y , e b z } the body fixed frame, whose origin coincides with the center of mass, see Fig. 1. The orientation of the rigid body is given by an orthogonal rotation matrix R ∈ SO(3) : I → B. The Newton-Euler equations of motion representing the position and orientation dynamics of the quadrotor are

Figure 1 :

 1 Figure 1: Inertial and body frames in the vehicle.

Figure 2 :

 2 Figure 2: The attitude is "sacrificed" in order to control the position.

Figure 3 :

 3 Figure 3: The final rotation about the z axis does not change the direction of the original desired force.

)

  Now the error for the position and observer dynamics are defined as xe = xξ -x d estimated position error xe = xξ -x ξ observed position error x e = x ξ -x d real position error (54)

Figure 4 :

 4 Figure 4: Block diagram of the control algorithm.

)

  and β d is the sideslip angle. The coefficients values used are ρ a = 1.19 kg/m 3 , β d = 0, c l = 0.15, c m = 0.25, c n = 0.8, c L = 0, c Y = 0.2, c D = 0.1, S = 0.5 m 2 , b = 0.05 m, c = 0.02 m. J = diag(0.1241, 0.1241, 0.2483) kg m 2 /rad and mass m = 2kg.

Figure 5 :Figure 6 :

 56 Figure 5: Position and velocity signals in the numerical results.

Figure 7 :

 7 Figure 7: Position error in the numerical results.

Figure 8 :

 8 Figure 8: Position control signals in the numerical results.

Figure 9 :

 9 Figure 9: Attitude signals in the numerical results.

Figure 10 :

 10 Figure 10: Attitude error signals in the numerical results.

Figure 11 :

 11 Figure 11: Orientation control signals in the numerical results.

Figure 12 :

 12 Figure 12: Torque and force disturbances in the numerical results.

Figure 13 :

 13 Figure 13: Platform used for the experimental validation: Parrot AR Drone 2.0.

19

 19

Figure 14 :Figure 15 :

 1415 Figure 14: Position and observer signals in the experimental validation.

Figure 16 :Figure 17 :

 1617 Figure 16: Comparison between dirty derivative and observer for the velocity estimation in the experimental validation.

Figure 18 :Figure 19 :

 1819 Figure 18: Attitude in the experiments.

Figure 20 :

 20 Figure 20: Angular velocity in the experiments.

Figure 21 :

 21 Figure 21: Angular velocity error in the experiments.

Figure 22 :Figure 23 :Figure 24 :

 222324 Figure 22: Attitude control signals during the experiment.

Acknowledgment

Authors acknowledge partial support from Conacyt-Mexico under Research Grants 264513, 133344 and 133546. This work has been sponsored by the French government research program "Investissements d'avenir" through the Robotex Equipment of Excellence (ANR-10-EQPX-44).


M A Gutierrez-Martinez

L E Cabriales-Ramirez

E U Rojo-Rodriguez

E J Ollervides-Vazquez

P Castillo
email: castillo@hds.utc.fr

O Garcia-Salazar

Genetic algorithm for path planning of UAVs as a maze-solving problem

This research work addresses the problem of path planning for unmanned aerial vehicles (UAVs) in complex environments such as a maze. A genetic algorithm (GA) with variable chromosomes and gene change conditions is designed when the decision point criteria and a collision with a wall are presented; without considering large chromosomes. Our proposed GA obtains the sequence of minimum movements required to solve the maze. Then, the trajectories are generated by high-order polynomials. Finally, numerical simulations and real-time tests are carried out to validate the proposed algorithm.

I. INTRODUCTION

Unmanned aerial vehicles (UAVs) have been used in multiple increasingly complex applications, such as agriculture, parcel delivery, and surveillance [START_REF] Tsouros | A review on uav-based applications for precision agriculture[END_REF], [START_REF] Huang | Optimal control of a hybrid uav/train parcel delivery system[END_REF], and [START_REF] Gu | Multiple moving targets surveillance based on a cooperative network for multi-uav[END_REF], due to their significant development and technological advances in control, intelligence techniques, and algorithms. Due to this, algorithms for optimal navigation have been of interest to researchers in robotics, mechatronics, and aeronautics fields. These algorithms have been widely studied and used in UAVs to achieve optimal navigation to make UAVs autonomous and efficient [START_REF] Rojo-Rodriguez | A fuzzy gain scheduling control algorithm for formation flight of multi-uavs[END_REF]. In [START_REF] Mikhaylov | Algorithm of autonomous uav orientation for applying in complex indoor environment[END_REF], authors worked on the task of exploring unknown environments. They aborded the construction of maps for indoors where the information on the environment was unknown. Then a Lee algorithm was implemented for reduced the memory cost to generate the map, and an A-star algorithm calculated the distance to a particular point. Thus, Lee and A-star algorithms calculated angles and distance to nodes from a point for the orientation of UAVs in a maze.

The tasks such as exploration, mapping, and recognizing unknown environments have been essential for autonomous navigation [START_REF] Manfreda | On the use of unmanned aerial systems for environmental monitoring[END_REF]. Thus, algorithms that solve path planning in complex environments as maze-like have been essential to achieve this navigation. In this sense, path planning algorithms were designed to execute missions considering obstacle penalties or a minimum distance of trajectory criteria [START_REF] Aggarwal | Path planning techniques for unmanned aerial vehicles: A review, solutions, and challenges[END_REF]. In [START_REF] Mac | Improved potential field method for unknown obstacle avoidance using uav in indoor environment[END_REF] and [START_REF] Mac | The development of an autonomous navigation system with optimal control of an uav in partly unknown indoor environment[END_REF], a modified potential field algorithm for path planning of UAV in an indoor environment was designed. The improved potential field algorithm was capable of avoiding unknown obstacles types of barriers and reaching the target of the trajectory. Therefore, the simulations in a virtual environment were performed, and experiments were implemented in real-time in a low-cost UAV to validate the proposed method. In [START_REF] Lupascu | Path planning for autonomous drones using 3d rectangular cuboid decomposition[END_REF], an algorithm for path planning based on cells decomposition was presented, in which authors implemented this algorithm to find trajectories in a 3D static environment for UAVs.

Classical methods were used in path planning in maze environments; however, heuristic and intelligent methods, among which the evolutionary and bio-inspired algorithms are used in cases where the classical methods usually presented problems. In [START_REF] Song | An improved pso algorithm for smooth path planning of mobile robots using continuous high-degree bezier curve[END_REF], a particle distribution algorithm (PSO) was used to find optimal trajectories, while by high degree bezier curves based on polynomial equations, authors smoothed the trajectory and made it continuous for tracking. Hence the proposed method was implemented in environments with barriers as obstacles, avoiding colliding with the walls and generating safe trajectories. In [START_REF] Konatowski | Ant colony optimization algorithm for uav path planning[END_REF] and [START_REF] Gutierrez-Martinez | Collision-free path planning based on a genetic algorithm for quadrotor uavs[END_REF], authors showed that bio-inspired and evolutionary methods such as ant colony and GA solve the path planning problem for UAVs by creating efficient trajectories to avoid collisions. In [START_REF] Choubey | A-mazer with genetic algorithm[END_REF], a GA was designed to solve complex maze-like environments through sequences of movements. However, the authors worked with large chromosomes of 400 genes, which make the function of genetic operators difficult. In [START_REF] Nagata | A genetic algorithm for the picture maze generation problem[END_REF], the authors used GA to find a maze puzzle solution in which used a novel Edge assembly crossover operator (EAX) is used, and it reduces computational time and finds trajectories that solve a maze picture environment. Finally, in [START_REF] Kim | Genetic algorithms for solving shortest path problem in maze-type network with precedence constraints[END_REF], the authors proposed a GA to find the shortest path in a maze-type network. The proposed algorithm solved the classical combinatorial optimization problems and showed that GA is a valuable tool for solving maze-like environments and path planning.

The main contribution of this work is a variable chromosome genetic algorithm to solve mazes in the least number of movements and with the shortest path distance. The motion sequence is converted to waypoints to generate UAV trajectories through high-order polynomials. Finally, these trajectories were simulated and implemented in a UAV to validate the proposed algorithm.

The organization of this work is as follows: section II presents the problem statement to be solved. Section III shows the equations of the mathematical model of UAVs. Section IV describes our genetic algorithms, the genetic operators with their percentages of crossover and mutation, the objective function, and the criteria of gene change. In section V, numerical simulations and real-time validation for the generated high-order polynomials are presented. The conclusions are present in section VI.

II. PROBLEM STATEMENT

This research focuses on the problem of creating a path in complex environments such as a maze. The optimization of the path is considered the criteria of minimum distance and the minimum number of movements since the shortest distance is not necessarily the path with the fewest movements. Thus, having many changes or turns in the path can be a problem for the trajectory generators and the control techniques; then, minimizing the number of movements is essential. A genetic algorithm is proposed to minimize the number of movements and the path distance. Solving a path in a maze becomes a combinatorial optimization problem because a specific order of movements is needed to find a path. For this reason, a genetic algorithm of variable chromosomes is proposed with a criterion of changes of gene only when the conditions of the intersection point, decision point, final point are accomplished; this helps to solve the path in the minimum number of movements and with a short variable of chromosomes.

III. EQUATIONS OF MOTION

The mathematical model of UAV considers the inertial frame I={x I , y I , z I }, the body frame B={x B , y B , z B } that is center of gravity fixed on the vehicle and the wind frame A={x A , y A , z A }; this frame is generated during the flight.

The formulation Newton-Euler for rigid bodies is used to obtain the mathematical model of quadrotor (2)- [START_REF] Rojo-Rodriguez | A fuzzy gain scheduling control algorithm for formation flight of multi-uavs[END_REF].

ξ = V (1) m V = Re 3 (-T T ) + mge 3 (2) Ṙ = R Ω (3) J Ω = -Ω × JΩ + τ a ( 4
)
where the position and velocity to the inertial frame relatives are ξ = (x, y, z) ∈ R 3 and V = ( ẋ, ẏ, ż) ∈ R 3 . η = (φ, θ, ψ) ∈ R 3 represents the rotation coordinates where φ, θ and ψ represent roll, pitch and yaw, respectively. The rotation matrix R ∈ SO(3) : B → I, it satisfies

SO(3) = {R|R ∈ R 3×3 , det[R] = 1, RR = R R = I}
this matrix is Euler angles parameterized φ, θ and ψ. Then the rotation matrix is written as

R =   c θ c ψ s φ s θ c ψ -c φ s ψ c φ s θ c ψ + s φ s ψ c θ s ψ s φ s θ s ψ + c φ c ψ c φ s θ s ψ -s φ c ψ -s θ s φ c θ c φ c θ   Ω = (p
, q, r) ∈ R 3 contains the angular velocity with respect to the body frame, and Ω is a skew-symmetric matrix, that is characterized by

Ω =   0 -r q r 0 -p -q p 0  
In addition, vectors e 1 , e 2 , and e 3 of the canonical basis of R 3 are considered. The inertia moments are considered on J ∈ R 3×3 , and m ∈ R denotes the mass of the quadrotor UAV. T T ∈ R 3 represents the total thrust, and τ a ∈ R 3 is the actuator moments of the quadrotor. The actuator moments τ a of four rotors are described as follows

τ a =   τ φ τ θ τ ψ   =   d[(f 2 + f 4 ) -(f 1 + f 3 )] d[(f 3 + f 4 ) -(f 1 + f 2 )] Q 2 + Q 3 -Q 1 -Q 4  
where d is the distance from the center of mass to the rotor axis, f i is the lift force or thrust force of the propeller for i = 1, 2, 3, 4 and

Q i = ρ a A di ω 2 i r 3 i c
Qi is the reaction moment of the rotors with c Qi that denotes the coefficient moments of the rotor shaft, ρ a is the air density, A di denotes the disk area, ω i denotes the angular velocity of the rotor and r i is the rotor radius for i = 1, 2, 3, 4.

IV. GENETIC ALGORITHM MAZE SOLUTION

A. Initialization of genetic algorithm

A genetic algorithm (GA) is proposed to solve a maze in a static environment. Our genetic algorithm works with a random population of real numbers instead of a binary population to reduce the chromosome length so that the genetic operators calculations are less complicated. Each gene of population is defined as A ∈ R, where A = {1, 2, 3, 4, 5}. In addition, a population of 100 individuals and 20 genes as maximum in the chromosome is implemented. The phenotype of the genetic algorithm is translated as state variables of motion for the resolution of the maze. Then each algorithm gene can take a value between 1 to 4; these states are translatable into the four possible movements, right, up, left, and down, respectively; also, a fifth state is added that represents no movement. In this sense, when the sequence of movements reaches its final point, it remains static. Also, to save processing and avoid using chromosome length, a variable chromosome criterion is used, where a case can be resolved in 3 movements; however, the same case can be resolved in 2 movements, it can be seen in Fig. 1. Our proposed GA aims to optimize the number of movements until reaching the final point; the order of the genes in the chromosome is essential because it translates into the order of movements in the maze. Moreover, in our proposed GA, the algorithm only evaluates the necessary genes until reaching the final point. The chromosome is variable and depends on the number of movements until reaching the final point. Thus, the algorithm works with movements controlled by the index of chromosome gene; in this proposed algorithm, the gene index only changes when the following conditions are achieved.

The first gene change condition to determine subsequent movement is considered if it hits a wall and the same gene movement is maintained; this is shown in Fig 2 . The second condition of gene change is that the sequence of movements is at a decision point, which means the points where different paths can be taken. The third change condition is when the final point is reached; in this case, there is a gene change in the chromosome of movement and a change in individual in the population. Fig. 3 shows the operation of the gene change conditions. It can be seen in Fig. 3 a) when the sequence of movements reaches the first decision point, it can take two accessible paths, which are up, and the right. It can also go down or even backward; however, these are paths penalized by distance and wall collision, which these paths are discarded in future iterations. Then, when the sequence of movements reaches the decision points, there is an increase in the index that controls the movement gene; however, this is the same value. Thus, the algorithm continues moving along the same path to the right. Then, the movements upon reaching the second decision point a change in the chromosome occurs again, this causes increasing the index that controls the movement gene. The next movement in the chromosome is towards up the path is takem up until the final point is reached, as shown in Fig. 3 b). When reaching the final point, it is observed in Fig. 3 c) a number five, is placed in the following gene of the movement chromosome, remaining at that point.

B. Objective function

The objective function considers the following criteria to be able to solve the maze, that the distance of the path is as short as possible by [START_REF] Mikhaylov | Algorithm of autonomous uav orientation for applying in complex indoor environment[END_REF], that it be solved in the least number of movements, and that the distance between the final point of the path and the final point of the maze is the minimum. In addition, penalty pen conditions are added for these sequences of movements that leave the maze and

dist j = (xp j,m -stopx) 2 + (yp j,m -stopy) 2 (5)
In Equation ( 5) dist is the distance from the last movement to the final point, xp and yp are the position of the sequence of movements, and it is updated every cycle. stopx and stopy are the positions of the final point. crash is a counter of times when a sequence of movements collides with a wall, and pen is a high penalty added for leaving the dimensions of the maze. Countm is a counter of movement changes, and Countp is a counter of steps taken to solve the maze until it reaches the desired point. Then, the sequence of movements that reach the final point satisfies the conditions; the sequence of movements that make it in fewer steps satisfies better the established conditions and minimizes the distance of the path. In this sense, the objective function is the following [START_REF] Manfreda | On the use of unmanned aerial systems for environmental monitoring[END_REF], where k i is a constant, with i = {1, 2, 3, 4}

f v j = dist j + k 1 * Countm j + k 2 * Countp j +k 3 * crash j + k 4 * pen j (6)
C. Genetic operators 1) Selection: For the selection of individuals, a deterministic criterion is used. Based on the fitness value of the movement sequences, the individuals are ordered from lowest to highest, and half of the population with the best fitness values are chosen for crossover.

2) Crossover: For the crossover operator, a double crossover point criterion is used. Fig. 4 shows the double crossover point criterion, where a crossover point P C splits the chromosome and sets a boundary between the secondary crossover points P csub1 and P csub2 where the crossover points change in each iteration. Then, because the entire selected population experiment crossover, it has 100 % crossover. 3) Mutation: The high mutation criterion is used with a percentage of 4 % where 80 % of the population experiences the mutation. However, the best individuals are preserved for future iterations to maintain some individuals without mutation and to have convergence. The high mutation provides greater exploration and a combination of movements. The mutation occurs by altering the chain of movements by introducing a new movement, a random number between 1 and 4.

The genetic algorithm, in general, can be observed in algorithm 1 also is shown the criteria of gene change in the chromosome of movements, where the variable i controls the index. Algorithm 2 shows the cases of movements and the counters to detect how many accessible paths there are. If there are two free paths, it is known that it is at the decision point, and this is counted by con1 while con2 is a counter to determine it is within the limits of the maze.

V. VALIDATION

A. Validation genetic algorithm

Two maze-type environments are prepared for six different cases to evaluate the genetic algorithm. Also, numerical simulations and experimental validation are executed to verify that the control algorithms can follow the high-order polynomials.

1) Case 1: In case 1, a maze is prepared as shown in Fig. 5. The green point is the initial point at 0.3m, 0.3m at x, y, and the red point is the final point at 1.2m and 3.9m at x, y. The blue dotted line is the sequence of movements by the genetic algorithm; it can be seen that the algorithm solves this case in 4 movements, this being the minimum number of movements to solve this case. The orange line is the trajectory response of the navigation points generated by the sequence of movements; a high-order polynomial generates this trajectory. This case is solved with seven genes of chromosomes.

2) Case 2: Case 2 is shown in Fig. 6 as in the previous case; it has an initial point at 0.3m and 0.3m at x, y, the final point is at 4.5m and 3.9m at x, y. It can be observed that the sequence of movements generated by the genetic algorithm passes through three decision points and an intersection point; also, the algorithm solves this case in the Algorithm 1: Genetic algorithm to solve maze.

Set variable P population size, N c number of chromosomes. Set P c crosrover rate,P m Mutation rate. Set counter i = 1, j = 1, g = 1 Initializes population M ove P,N c while g <= gene do

Objective function while j <= P do xp = startx; yp = starty; The high-order polynomial generates the orange trajectory, and the points generated by the sequence of movements; are obtained by the GA. This case is solved with six genes of chromosomes.

while i <= N c do n = move j,i Algorithm2(n) if con2 >= 4 then if con1 < 1 then i = i + 1 end end if xp < 0 or yp < 0 or xp > pf x or yp > pf y then f v j = f v j + p i = i + 1 end if xpf = obsx and yf p = obsy then i = i + 1 end if xp = obsx and yp = obsy then crash j = crash j + 1 i = i + 1 end ic = i if xp = stopx
3) Case 3: Case 3 is presented in Fig. 7 in which another maze-like environment is presented with more complicated paths and more decision points. The initial green point is at 0.10m, 0.10m, at x, y respectively, and a final redpoint is at 0.3m y 2.6m in x, y. The path, generated by the sequence of movements, is a blue line dotted. It can be seen that the movements pass through five decision points, and these could have deviated or made unnecessary movements. However, the algorithm finds the appropriate sequence of movements to arrive at the final point with the minimum number of movements, which are six movements. After having the sequence of movements, the sequence of movements is discretized to obtain navigation points; these navigation points obtained are used to generate the continuous trajectory through the high- order polynomial, which can be observed in orange color. This case is solved with nine genes of chromosomes.

4) Case 4:

In case 4, it can be seen that the initial point is at 0.1m, 0.1m in x, y and the final point is at 2.9m, 2.6m, in x, y; the final point is set as far from the inital point to test the algorithm, it can be seen in Fig. 8. The path generated by the sequence of movements is the blue dotted line. It can be seen that the algorithm can go through many decision points, so it is difficult to reach the final point and even more difficult to arrive at the least number of movements. Then, it can be noted that the algorithm finds the trajectory in the least number of movements, which is seven. This case is solved with twelve genes of chromosomes.

5) Case 5: Case 5 is shown in Fig. 9. In this case, the initial point and final point are placed at different extremes, where the initial point is at 2.8m and 0.1m at x, y while the final point is at 0.3m and 2.6m the algorithm solves the maze in the minimum number of movements, which in this case are six. As in the previous cases, the points generated by the sequence of movements are used, and the highorder polynomial is obtained to generate the trajectory, the blue dotted line is the trajectory generated by the sequence of movements. On the other hand, the orange line is the continuous trajectory. This case is solved with nine genes of chromosomes. 6) Case 6: Finally, case 6 is shown in Fig. 10. It can be seen that the initial point is at 1.4m in x and 0.2m in y, the blue dotted line represents the sequence of movements made by the GA. The algorithm could have arrived by the path on the right so that these do not to have to move down. However, the algorithm evaluates and decides to take the path on the left even if it has to move down, this is because the path on the left is solved in fewer steps and fewer movements. The orange path is the one made by the high-order polynomial and the navigation points obtained from the sequence of movements. This case is solved with ten genes of chromosomes. The fitness value graphs are shown in Fig. 11; these graphs show the evolution of the algorithm until the solution is found. The figure shows the fitness value graphs for the case 1,2,3,4,5 and 6 corresponding to letters a), b), c), d), e) and f). It can be seen that cases one and two corresponding to maze 1 (Fig. 11 a) and b) ) are solved in a few less than 20 iterations, because maze 1 has fewer decision points. Also, it can be seen that most of the cases are solved in almost 20 iterations; the only case that took more iterations to solve is case e) corresponding to case 5, where this case is more difficult to solve due to the number of decision points and possible movements to solve the case.

0 1 2 3 X Axis[m] 0 1 2 3 Y Axis[m]

B. Simulation

In this subsection the simulations are carried out to validate the proposed algorithm a numerical validation of the trajectories generated with GA and high-order polynomials. The UAV used to obtain the model is a miniquadrotor that is part of a fleet of aerial vehicles developed in the at CIIIA-FIME-UANL [START_REF] Ollervides-Vazquez | Twolayer saturated pid controller for the trajectory tracking of a quadrotor uav[END_REF]. In the Table I the numerical values of the model parameters are shown.

Parameter Value Parameter Value cT 9 × 10 -6 [N -s 2 ] d 0.2 [m] cQ 3 × 10 -7 [N -m -s 2 ] Jxx 15 × 10 -3 [kg -m 2 ] m 0.3 [kg] Jyy 15 × 10 -3 [kg -m 2 ] g 9.81 [m/s 2 ] Jzz 30 × 10 -3 [kg -m 2 ]
TABLE I: Quadrotor parameters.

The parameters are as follows: c T and c Q are aerodynamic coefficients, and d is the distance of the center of gravity of the quadrotor with respect to the axis of rotation of the propellers.

1) Case 1: The Fig. 12 shows the top view of the trajectory UAV performed in the first maze environment. The trajectory executed for the simulation is the same as the reference by the high-order polynomial. 2) Case 2: The Fig. 13 shows the top view of the trajectory UAV performed and implemented in the first mazetype. Even having four turns in the trajectory, the simulation illustrates that the UAV tracks the reference.

3) Case 3: The Fig. 14 shows the 2D trajectories, the trajectory performed for UAV, and the reference, implemented for the second maze environment. The simulation shows the performance of the trajectory tracked for the UAV. 4) Case 4: The Fig. 15 shows the top view of the trajectory UAV performed. Even having six turns in the trajectory, the simulation shows that the UAV performs the reference signals.

5) Case 5: The Fig. 16 illustrates the 2D of the trajectory UAV performed. As in the previous cases, the simulation demonstrates the tracking of the reference trajectory without any alteration.

6) Case 6: The Fig. 17 shows the 2D trajectories, the trajectory response for UAV and the reference. These 4 cases are considered the second maze-type.

In all above graphics, it can be seen that the trajectories generated from the GA and high-order polynomials can be executed out on a UAV. Reference UAV Fig. 17: Trajectory in the XY -axis for case 6.

2) Experimental platform: Small form-factor quadrotor is selected as experimental platform as well as motion capture system. The testing facilities are located at the Navigation Laboratory, of the Aeronautical Engineering Research and Innovation Center (CIIIA), of the Autonomous University of Nuevo Leon (UANL). The motion capture system consists of 16 VICON T40 cameras, streaming and orientation data to a centralized computer at 100Hz, with an accuracy of 0.1mm and 0.1 • . The quadrotor aircraft used is designed and built for testing purposes in the same facilities, shown at Fig. 18, as basic characteristics it is based on "X" configuration, with 12.5cm arms and 1.9N thrust for each rotor. Individual axis signals are presented for X -axis at Fig. 21 and Fig. 22 for Y -axis. Maze trajectory is initialized approximately at 25s and finalized at 45s, the rest of time is set for take off and landing process. The Fig. 21 and Fig. 22 it is presented the follow-up of the trajectory evolving in time.

4) Case 2:

A simple scenario is presented as second experiment, 3 decision points are performed to accomplish the trajectory. A constant height of 1500mm, X = -2000mm and Y = -900mm as initial point is established, and X = -1200mm and Y = 2700mm for final solution. A trajectory generation is accomplished by the genetic algorithm and implemented using a high-order polynomial. Fig. 23 and Fig. 24 show the ideal and actual trajectory and for 2D perspective in Fig. 24 it is appreciated the maze which is actually the same as experiment 1 but with other final point. The Fig. 25 and Fig. 26 present X, Y -axes respectability for individual signal analysis. Maze trajectory is evolved approximately from 30s to 60s. This is a more aggressive trajectory, showing the effects on X -Y axes, some disturbance is presented due to aerodynamic effects and navigation control tracking. Nonetheless, the performances exhibited is acceptable for the mission. This experimental test is captured on video and is available on the following link https://www.youtube.com/ watch?v=uCr2W3Zdnac

VI. CONCLUSIONS

The proposed algorithm has been efficient for path planning in complex maze-like environments. The criteria of change of gen in the chromosome at the decision and the intersection points have solved the six cases proposed with few genes. The fitness value graphs have showed that the minimization of the objective function has been achieved by solving the cases in a few movements. In addition, the method to obtain a trajectory with high order polynomials has generated smooth trajectories so that the UAV quadrotor tracks these paths in real-time.

Fig. 1 :

 1 Fig. 1: Scheme of movements.

Fig. 2 :

 2 Fig. 2: Change of direction by the intersection point.

Fig. 3 :

 3 Fig. 3: Change of direction by decision point.

Fig. 4 :

 4 Fig. 4: Double crossover point.

  and yp = stopy then move(i + 1) = 5; i = N c end m = m + 1; end countp = m; countm = ic Calculate distance by Equation (5) j = j + 1; Evaluate objective function with (6) end Selection operator Crossover operator Mutation operator Replace population g = g + 1; end minimum number of movements which for this case is four.

Algorithm 2 :

 2 Cases algorithm. switch n do Case 1 Move right xpf = posx + (dp * 2); ypf = yp; xp = xp + dp; yp = yp; Case 2 Move up xpf = xp; ypf = yp + (dp * 2); xp = xp; yp = yp + dp; Case 3 Move left xpf = xp -(dp * 2); ypf = yp; xp = xp -dp; yp = yp; Case 4 Move down xpf = xp; ypf = yp -(dp * 2); xp = xp; yp = yp -dp; Case 5 No move con1 = 0; con2 = 0; x 1 = xp + δp; y 1 = yp + δp; x 2 = xp; y 2 = yp + δp; x 3 = xp -δp; y 3 = yp; x 4 = xp; y 4 = yp -δp; x 5 = xp + δp; y 5 = yp + δp; x 6 = xp + δp; y 6 = yp -δp; x 7 = xp -δp; y 7 = yp + δp; x 8 = xp -δp; y 8 = yp -δp; while n <= 8 do if x n = obsx and y n = obsy then con1 = con1 + 1; if i >= 5 then con2 = con2 + 1;

Fig. 5 :

 5 Fig. 5: Case 1.

Fig. 6 :

 6 Fig. 6: Case 2.

Fig. 7 :

 7 Fig. 7: Case 3.

Fig. 8 :

 8 Fig. 8: Case 4.

Fig. 9 :

 9 Fig. 9: Case 5.

Fig. 10 :

 10 Fig. 10: Case 6.

Fig. 11 :

 11 Fig. 11: Fitness value.

Fig. 12 :

 12 Fig. 12: Trajectory in the XY -axis for case 1.

Fig. 13 :Fig. 14 :

 1314 Fig. 13: Trajectory in the XY -axis for case 2.

C. Implementation 1 )

 1 Real-time experiments: Real-time experiments are performed for path planing demonstration at known space. Two trajectories are determinated offline and commanded to the quadrotor aircraft by a time-dependent continuous polynomial function.

Fig. 15 :Fig. 16 :

 1516 Fig. 15: Trajectory in the XY -axis for case 4.

Fig. 18 :

 18 Fig. 18: Quadrotor used for experimental tests.

3) Case 1 :Fig. 19 :Fig. 20 :

 11920 Fig. 19: 3D trajectory for first experimental scenario.

Fig. 21 :Fig. 22 :

 2122 Fig. 21: Trajectory in the X-axis for first experimental scenario.

Fig. 23 :

 23 Fig. 23: 3D trajectory for second experimental scenario.

Fig. 24 :

 24 Fig. 24: Trajectory in the XY -axis for second experimental scenario.

Fig. 25 :Fig. 26 :

 2526 Fig. 25: Trajectory in the X-axis for second experimental scenario.

Aerospace Engineering Research and Innovation Center, Faculty of Mechanical and Electrical Engineering, Autonomous University of Nuevo Leon, Apodaca NL, Mexico. manuel.gutierrezmrt@uanl.edu.mx;

ACKNOWLEDGMENT

This research work was supported by the Office of Naval Research Global through the grant number N62909-20-1-2030.