Ngoc Thinh Nguyen

Ionela Prodan

Laurent Lefèvre

ON CONTROL SYSTEMS TECHNOLOGY 1 Stability guarantees for translational thrust-propelled vehicles dynamics through NMPC designs

Keywords: NMPC, terminal constraint, unconstrained NMPC, feedback linearization, ellipsoidal invariant set, asymptotic stability, VTOL (Vertical Take-Off and Landing) vehicles

We propose two different NMPC (Nonlinear Model Predictive Control) schemes (with and without terminal cost and constraints) for stabilizing the translational dynamics of thrustpropelled vehicles. Both approaches make use of an elaborated nonlinear feedback linearization controller and its associated ellipsoidal invariant set under restrictive input constraints, hence, guaranteeing the closed-loop asymptotic stability. The terminal constraint set of the corresponding NMPC design is easy to tune due to its clear formulation expressed directly in terms of the tuning variables while for the NMPC scheme without terminal constraint, the design allows to stabilize the system with a significantly shorter prediction horizon in comparison with existing method in the literature. Simulation and experimental tests over a nano-drone platform validate the proposed approaches.

I. INTRODUCTION

N OWADAYS, there is a tremendous interest in the re- search community on UAVs (Unmanned Aerial Vehicles) and thrust-propelled vehicles (like quadcopters, hexacopters), in particular [START_REF] Zanelli | Nonlinear model predictive control of a human-sized quadrotor[END_REF], [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF]. Usually, the drone platforms already have their built-in controllers which control the rotors to track the four inputs which in turn govern the thrust level and the three Euler angles [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF]- [START_REF] Nguyen | Effective angular constrained trajectory generation for thrust-propelled vehicles[END_REF]. Therefore, more recently, controlling the thrust propelled vehicles requires only the maneuvering of their translation dynamics [START_REF] Nguyen | Effective angular constrained trajectory generation for thrust-propelled vehicles[END_REF]. This is already a challenging task as the dynamics are not only strongly nonlinear but also subject to many operating constraints [START_REF] Zanelli | Nonlinear model predictive control of a human-sized quadrotor[END_REF]- [START_REF] Budaciu | On the evaluation of the crazyflie modular quadcopter system[END_REF]. Designing a controller respecting the system dynamics and constraints can be fulfilled by considering an Model Predictive Control (MPC) [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF] approach. This control strategy provides the input obtained by solving a finite horizon optimal control problem where the dynamics and constraints are explicitly taken into account. It is well known that the computational burden of the MPC-based optimization problem may pose serious issues in many time-critical applications. However, with recent technological advances, it becomes possible to employ MPC for many applications requiring fast response [START_REF] Badgwell | Model-predictive control in practice[END_REF] as is the case in the control and coordination of aerial vehicles [START_REF] Zanelli | Nonlinear model predictive control of a human-sized quadrotor[END_REF], [START_REF] Mueller | A model predictive controller for quadrocopter state interception[END_REF]. E.g., [START_REF] Zanelli | Nonlinear model predictive control of a human-sized quadrotor[END_REF] approximates the NMPC (Nonliner MPC) of stabilizing a human-sized quadcopter vehicle with a N.T. Nguyen is with the Univ. of Lübeck, Institute for Robotics and Cognitive Systems, 23562 Lübeck, Germany (email: nguyen@rob.uni-luebeck.de).

I. Prodan and L. Lefèvre are with Univ. Grenoble Alpes, Grenoble INP Institute of Engineering and Management, LCIS, 26000 Valence, France (email: {ionela.prodan,laurent.lefevre} @ lcis.grenoble-inp.fr). This work was done while N.T. Nguyen was pursuing his PhD degree at laboratory LCIS. quadratic problem and succeeds in solving it by using a Cortex A9 800-MHz microcontroller. [START_REF] Mueller | A model predictive controller for quadrocopter state interception[END_REF] employs an MPC scheme to generate quadcopter trajectories in real time by formulating the problem using the system's jerks. Also, it has been proved that a linear MPC-based optimization can be even executed over a 8-bit microcontroller by using a convex lifting approach [START_REF] Gulan | Embedded linear model predictive control for 8bit microcontrollers via convex lifting[END_REF]. Another important issue of MPC is that it does not provide an explicit controller equation [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF], thus, causing difficulty in analyzing the stability and the feasibility of the closed-loop scheme (except explicit MPC whose numerical limitations make it difficult for many real-time implementations [START_REF] Bemporad | The explicit linear quadratic regulator for constrained systems[END_REF], [START_REF] Munir | Low complexity constrained control using higher degree lyapunov functions[END_REF]). The use of MPC laws without stability guarantees may prove shortsighted in certain applications. In particular, the authority loss for attitude and position control in UAVs will not only risk the mission objectives but lead to critical behavior. In the literature, the stability of an MPC controller can be achieved either by adding a terminal constraint [START_REF] Munir | Low complexity constrained control using higher degree lyapunov functions[END_REF]- [START_REF] Cannon | Nonlinear model predictive control with polytopic invariant sets[END_REF] or by enlarging the prediction horizon [START_REF] Grüne | NMPC without terminal constraints[END_REF]. In this paper, we employ both approaches to design two different NMPC schemes for stabilizing the translation dynamics of thrustpropelled vehicles. These designs make use of an elaborated feedback linearization controller which was introduced in the previous work of the authors [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF] and was successfully applied to a nano-drone. In here, we employ this controller as the local controller (i.e., locally constraint admissible) to design the two NMPC schemes with guaranteed stability. Note that the role of a local controller is to ensure that once inside the terminal region, the system's trajectories remain inside. Thus, the local controller only acts as a guarantee of feasibility and asymptotic stability and is not actually used throughout the simulation. Next we enumerate this paper's contributions which, to the best of our knowledge, are new to the state of the art. We: i) propose an NMPC design with terminal ellipsoidal region which is invariant and constraint admissible under a nonlinear feedback linearization controller. ii) propose an NMPC design without terminal constraint and the corresponding tuning procedure for choosing the optimal prediction horizon which guarantees the stability. iii) validate the NMPC controllers through simulation and experimental testing over the Crazyflie 2.0 nano-quadcopter platform [START_REF] Budaciu | On the evaluation of the crazyflie modular quadcopter system[END_REF].

The paper structure is as follows. Section II presents first the translational thrust-propelled systems, then, the principles of the two NMPC design schemes, i.e., NMPC with and without terminal constraint. Section III details several properties of the elaborated feedback linearization controller which are employed in Section IV to ensure the stability of the two NMPC designs. Section VI presents the simulation results and Section VII shows the corresponding experimental validation. Section VIII presents the conclusions and the future work.

Notations: Denote by I n the identity matrix of size n and by 0 m×n the zero matrix of size m × n. 0 represents a zero matrix of appropriate dimension inferred from the context. Let λ min (A) and λ max (A) denote the minimum and maximum eigenvalues of the square matrix A having only real eigenvalues. For a vector x ∈ R n and a positive semidefinite matrix P ∈ R n×n , x2 P denotes the weighted norm x P x while x represents the Euclidean norm of the vector x. Furthermore, denote by x(i|k) and u(i|k) the values of the state x and the input u at time instant i predicted upon the information available at time k ∈ N.

II. MODEL DYNAMICS AND NMPC SETUPS A. System modeling

This section recapitulates the discrete translation dynamics of a standard thrust-propelled aircraft system. The dynamics are obtained by applying the Runge-Kutta fourth-order discretization method [START_REF] Butcher | Runge-kutta methods: some historical notes[END_REF] to the continuous thrust-propelled system found in [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF] and are expressed in their state-space representation as follows:

x k+1 = f (x k , u k ) = Ax k + h ψ (u k ), (1)
where the state x k ∈ R 6 gathers the position and the velocity at the time step k,

x k [x k y k z k v x k v y k v z k ] , the input u k ∈ R 3 , u k [T k φ k θ k ]
consists of the normalized thrust T (the normal thrust divided by the system's mass), the roll, φ, and pitch, θ, angles. The yaw angle ψ ∈ [-π, π] is an assumed known constant influencing the system (usually considered to be zero as in [START_REF] Lu | Constrained anti-disturbance control for a quadrotor based on differential flatness[END_REF] but here we aim to exploit a general fix value of ψ). Explicitly, the matrix A and the input terms h ψ (u k ) are given by:

A = I 3 ∆ t I 3 0 3×3 I 3 , (2)
h ψ (u k ) =   ∆ 2 t 2 I 3 ∆ t I 3     T k (c φ k s θ k c ψ + s φ k s ψ) T k (c φ k s θ k s ψ -s φ k c ψ) -g + T k c φ k c θ k   , (3)
where ∆ t is the discretization time step and g stands for the gravitational acceleration 1 . Note that, the discrete system (1) is linear in the state since the original continuous system taken from [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF] already possesses this property. The vehicle has to respect thrust limits, with T limit > g, and roll and pitch constraints as follows:

u(k) ∈ U = {0 ≤ T k ≤ T limit , |φ k | ≤ c , |θ k | ≤ c }, (4)
with c ∈ (0, π/2), the desired maximum value of the angles. Then, without loss of generality, we aim to stabilize the system (1) around the equilibrium point:

x e = 0, u e = [g 0 0] . (5)
Note that, the thrust-propelled system (1) can be stabilized at any position x e ∈ R 3 with zero velocity by simply choosing the origin of the coordinate system at x e .

B. NMPC with and without terminal constraint

Model Predictive Control (MPC) is a control strategy in which at each time step, an open-loop optimal control problem of prediction horizon N , subject to the system dynamics (1) and constraints on states (including the current state as the initial condition) and inputs (4), is solved to obtain an optimal control sequence. From the sequence, only the first element is applied to the system. At the next sampling instant, the state is measured again and introduced to the optimization problem. The process is iteratively executed to establish the closed-loop controlled system [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF]. Hereinafter, we introduce two different NMPC setups, i.e., with and without terminal constraint for the considered dynamics (1) which are adopted from their general designs given in [START_REF] Chen | A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability[END_REF], [START_REF] Grüne | NMPC without terminal constraints[END_REF]. We also recapitulate the corresponding design principles for guaranteeing the stability of the two NMPC schemes [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF], [START_REF] Grüne | NMPC without terminal constraints[END_REF].

1) NMPC with terminal constraint (NMPC-t): The NMPC design detailed hereby employs both the terminal cost and the terminal constraint 2 . For more details, the open-loop optimization problem V T,N T (x k ) at time step k with the measured state x k is solved over the prediction horizon N T :

V T,N T (x k ) := min u(•|k) J T,N T (x k , u(•|k)), (6)
subject to

x(i + 1|k) = f (x(i|k), u(i|k)), (7a) u
(i|k) ∈ U, i ∈ {0, . . . , N T -1}, (7b) x(0|k) = x k , x(N T |k) ∈ X f , (7c)
with the cost function J T,N T (x k , u(•|k)) defined as:

J T,N T (x k , u(•|k)) = N T -1 i=0 (x(i|k), u(i|k)) + F (x(N T |k)).
The predicted state and input sequences employed for the optimization problem at time step k are denoted by x(•|k) [x(0|k) . . . x(N T |k)] and u(•|k) [u(0|k) . . . u(N T -1|k)], respectively. The terminal region X f ⊂ R 6×6 is designed according to the stability conditions introduced hereinafter. Furthermore, the stage cost, (x, u), and terminal cost, F (x), are defined as:

(x, u) = x -x e 2 Q + u -u e 2 R , (8)
F (x) = x -x e 2 P , (9)
in which, the symmetric matrices Q ∈ R 6×6 (positive definite) and R ∈ R 3×3 (positive semi-definite) are defined by the users, while the positive definite matrix P ∈ R 6×6 is a solution of the Lyapunov equation introduced later (see (32)). Note that, the formulation of (•) in ( 8) allows R = 0 as similar to [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF], [START_REF] Limón | On the stability of constrained mpc without terminal constraint[END_REF]. Assuming that the optimization problem ( 6)-( 7) is successfully solved, it provides the optimal state and input trajectories (x * (•|k), u * (•|k)). Then, the MPC control action applied to the system (1) at time step k is defined as:

µ T (x k ) = u * T (0|k). (10)
Finally, the nominal closed-loop system is given by:

x k+1 = f (x k , µ T (x k )). (11)
Hereinafter, we summarize the NMPC principles given in [START_REF] Chen | A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability[END_REF], [START_REF] Grüne | Nonlinear Model Predictive Control: Theory and Algorithms[END_REF] to establish the sufficient conditions for the stability of the nominal closed-loop system [START_REF] Chen | A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability[END_REF]. Assumption 1: Consider the closed terminal set X f associated to [START_REF] Badgwell | Model-predictive control in practice[END_REF] and the terminal cost F from [START_REF] Bemporad | The explicit linear quadratic regulator for constrained systems[END_REF]. We assume that for each x k ∈ X f at any time step k, there exists an admissible control value3 u loc (x k ) ∈ U such that the two following conditions hold:

x k+1 = f (x k , u loc (x k )) ∈ X f , (12)
F (x k+1 ) + (x k , u loc (x k )) ≤ F (x k ). ( 13
)
Theorem 1: Let Assumption 1 hold. Then, the nominal closed-loop system [START_REF] Chen | A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability[END_REF] under the feedback law µ T (10) from the NMPC-t scheme ( 6)-( 7) achieves recursive feasibility 4 and asymptotic stability.

Proof: See Chapter 5 in [START_REF] Grüne | Nonlinear Model Predictive Control: Theory and Algorithms[END_REF] and Theorem 1 in [START_REF] Chen | A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability[END_REF]. 2) NMPC without terminal constraint (u-NMPC): The NMPC scheme without terminal cost and terminal stabilizing constraint is also referred as unconstrained NMPC as in [START_REF] Grüne | Nonlinear Model Predictive Control: Theory and Algorithms[END_REF], [START_REF] Reble | Unconstrained model predictive control and suboptimality estimates for nonlinear continuous-time systems[END_REF], hence, being denoted by u-NMPC within the paper. The following optimization control problem V U,N U (x k ) at time k is solved over the prediction horizon N U :

V U,N U (x k ) := min u(•|k) J U,N U (x k , u(•|k)), (14)
subject to

x(i + 1|k) = f (x(i|k), u(i|k)), (15a) u(i|k) ∈ U, i ∈ {0, . . . , N U -1}, (15b) x(0|k) = x k , (15c)
with the cost function J U,N U (x k , u(•|k)) defined in terms of the stage cost from (8) as:

J U,N U (x k , u(•|k)) = N U -1 i=0 (x(i|k), u(i|k). (16)
Then, the u-NMPC control action and the nominal closed-loop system at time k are given by:

µ U (x k ) = u * U (0|k), (17)
x k+1 = f (x k , µ U (x k )), (18)
with x * U (•|k) the optimal input trajectories. The stability of system [START_REF] Chen | A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability[END_REF] is well studied and presented in [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF], [START_REF] Grüne | Nonlinear Model Predictive Control: Theory and Algorithms[END_REF] and summarized hereinafter.

Assumption 2: There exist constants γ, c ∈ R >0 such that for any N U ≥ 2 and for all the initial state x 0 satisfying x 0x e 2 Q ≤ c, we have:

V U,N U (x k ) ≤ γ x 0 -x e 2 Q . (19)
Theorem 2: Let Assumption 2 hold. Then, there exists an N 0 ∈ N, such that for the prediction horizon length N U ≥ N 0 , the equilibrium point x e is uniformly exponentially stable under the nominal closed-loop dynamics [START_REF] Chen | A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability[END_REF] for any initial state x 0 satisfying V U,N U (x 0 ) ≤ cγ.

Proof: See Theorem 1 in [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF] and Theorem 3.6 in [START_REF] Grüne | NMPC without terminal constraints[END_REF]. At first, in [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF], the authors show that

V N,U (x k ) ≤ cγ implies that V U,N U (x k ) ≤ γ x k -x e 2
Q with a case dictinction based on whether x kx e 2 Q ≤ c or not [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF]. Secondly, whenever V U,N U (x k ) ≤ γ x k -x e 2 Q holds, in [START_REF] Grüne | NMPC without terminal constraints[END_REF], the authors show that V U,N U (x k ) decreases for all N U ≥ N 0 with N 0 given by (see Variant 3 in [START_REF] Grüne | NMPC without terminal constraints[END_REF] for more details):

N 0 = 2 + ln(γ -1) ln γ -ln(γ -1) . ( 20
)
Lastly, the recursive feasibility and exponential stability are obtained with the initial condition V U,N U (x 0 ) ≤ cγ. Remark 1: One can impose state constraints into the two NMPC controllers by adding the following state condition into (7b) and (15b):

x(i|k) ∈ X ,

with X ⊆ R 6 , the non-empty state constraint set required to be convex and to contain the equilibrium point x e . Then, Theorems 1 and 2 accordingly require X f ⊆ X and

x k ∈ X , ∀ x k -x e 2 Q ≤ c.
Within this paper, we do not consider state constraints for the thrust-propelled dynamics (1) even though the proposed approaches provide the means to scale the terminal region X f (by tuning its radius as will be discussed in Section III) and the set of x 0x e 2 Q ≤ c (by increasing the value of Q). This is due to the fact that we want to analyze the maximum performance (convergence speed -shown by the velocity) of the NMPC controllers. Furthermore, the initial condition has to be chosen close to the equilibrium point for guaranteeing the stability (c.f. Theorem 2) due to hardware limitations (with largest possible prediction horizon of only 10 steps), hence, realistic constraints on position are never activated under our experimental tests.

C. Motivation

Assumptions 1 and 2, which are required to design the two NMPC constructions provided in Sections IV-A and IV-B are hard to validate for a generic nonlinear system given as in (7a) and (15a). In the literature, the standard approach is to linearize the dynamics around the equilibrium point, from which it follows the construction of a linear controller and its associated invariant set (of ellipsoidal form as in [START_REF] Chen | A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability[END_REF], [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF] or polyhedral form as in [START_REF] Cannon | Nonlinear model predictive control with polytopic invariant sets[END_REF]). However, employing a linear controller for a general nonlinear system (7a),(15a) obviously restricts the corresponding invariant set (also serving as the terminal region X f in (7)), hence, arguably reducing the efficiency of the NMPC-t controller ( 6)- [START_REF] Munir | Low complexity constrained control using higher degree lyapunov functions[END_REF] and leading to an impractically large prediction horizon for the u-NMPC controller ( 14)- [START_REF] Limón | On the stability of constrained mpc without terminal constraint[END_REF]. Thus, it is worthwhile to ask whether, for particular dynamics (as those shown in ( 1)), we may dispense with the linearized dynamics/linear controller construction and, instead, check Assumptions 1 and 2 by applying a feedback linearization controller. In our opinion this helps not only to account for the nonlinearities of the system (1) but also to apply various tools available for linear dynamics when considering the resulting closed-loop linear system (i.e., invariant set construction becomes much easier, as the difficulties are "pushed" in the nonlinear input). As is to be expected, the simplified closed-loop dynamics come at the price of a complex input formulation which, in turn, makes the input constraint and stage cost (x, u) in ( 4),( 8) formulations difficult to handle. These difficulties will be tackled by employing previously obtained results of the authors [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF], [START_REF] Nguyen | A stablizing nmpc design for thrust-propelled vehicles dynamics via feedback linearization[END_REF] and which will be recapitulated in the forthcoming sections.

III. FEEDBACK LINEARIZATION LAW AND INPUT CONSTRAINT SATISFACTION

This section recapitulates the results of a feedback linearization controller for quadcopter trajectory tracking [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF], [START_REF] Nguyen | Effective angular constrained trajectory generation for thrust-propelled vehicles[END_REF] which will be employed as the local controller for guaranteeing the stability of the two NMPC controllers detailed in Section II-B. In general, feedback linearization control approach reduces the controlled system to an equivalent linear system (assuming no mismatches on state feedback and system dynamics, as similar to nominal NMPC applications in ( 10) and ( 17)). Thus, the method does not require any approximation and does not restrict the system to operate only in the neighborhood of the linearization point [START_REF] Craig | Introduction to robotics: mechanics and control[END_REF]. Let us start by introducing the standard feedback linearization law

u loc (x k ) T loc (x k ) φ loc (x k ) θ loc (x k
) defined as follows (more details can be found in [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF], [START_REF] Nguyen | Effective angular constrained trajectory generation for thrust-propelled vehicles[END_REF]):

T loc (x k ) = u 2 x (x k ) + u 2 y (x k ) + (u z (x k ) + g) 2 , (22a) φ loc (x k ) = arcsin u x (x k ) s ψ -u y (x k ) c ψ T loc (x k ) , (22b)
θ loc (x k ) = arctan u x (x k ) c ψ + u y (x k ) s ψ u z (x k ) + g , (22c)
where the virtual control inputs u q (x k ) with q ∈ {x, y, z} are calculated in terms of the positions and velocities (i.e.,

x k [x k y k z k v x k v y k v z k ] from (1)): u q (x k ) = K 1q q k + K 2q v q k , (23)
with K 1q , K 2q ∈ R the control gains. If the feedback linearization control action u loc (x k ) as in ( 22) is not affected by the saturation limit, i.e., u loc (x k ) ∈ U with U in (4), then u loc (x k ) transforms the dynamics (1) into the following linear system:

x k+1 = A cl x k , (24)
with the matrix A cl ∈ R 6×6 given by:

A cl = I 3 + ∆ 2 t 2 K 1 ∆ t I 3 + ∆ 2 t 2 K 2 ∆ t K 1 I 3 + ∆ t K 2 , (25)
and with

K 1 = diag(K 1x , K 1y , K 1z ) and K 2 = diag(K 2x , K 2y , K 2z ).
Furthermore, the control gains K 1q , K 2q are required to satisfy the following conditions:

- 2 ∆ t < K 2q < ∆ t 2 K 1q < 0, (26)
which serve to establish the asymptotic stability of the resulted linear system (24) according to the Routh-Hurwitz criterion.

Then, A cl is a Schur matrix having all its eigenvalues strictly inside the unit circle. However, the control action u loc (x k ) has to respect the input constraints (i.e., u loc (x k ) ∈ U (4)) for the linearization to be valid. Thus, hereinafter, we will introduce some conditions on the virtual control inputs u x (x k ), u y (x k ) and u z (x k ) in [START_REF] Wächter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF] to ensure the input constraint satisfaction. Proposition 1: By choosing the values of the three positive saturation limits U x , U y and U z such that:

U z < g, ( 27
) U 2 x + U 2 y ≤ (-U z + g) 2 tan 2 c , (28)
U 2 x + U 2 y + (U z + g) 2 ≤ T limit , (29)
we have that, if the virtual inputs (23

) satisfy |u x (x k )| ≤ U x , |u y (x k )| ≤ U y and |u z (x k )| ≤ U z ,
u loc (x k ) -u e 2 ≤ x k (K K + 2Γ)x k , (30)
where

K = [K 1 K 2 ]
(24) and Γ ∈ R 6×6 is given by:

Γ = 1 (-U z + g) 2 K xy K xy , (31)
with

K xy = diag(K 1x , K 1y , 0) diag(K 2x , K 2y , 0) . Proof:
The proof of point (i) is presented in previous works of the authors [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF], [START_REF] Nguyen | A stablizing nmpc design for thrust-propelled vehicles dynamics via feedback linearization[END_REF] and the proof of point (ii) is given in Appendix A. The goal of the condition ( 27) is to constrain the thrust-propelled dynamics (1) not to free fall. Also, (28) and (29) are to guarantee the constraints on the angle (|φ loc |, |θ loc | ≤ c ) and the thrust (T loc ≤ T limit ), respectively, with c , T limit from (4). Further discussions on the origin and the advantages/ disadvantages of ( 27)-(29) are given in [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF]. By Proposition 1, we guarantee the input constraint satisfaction for the feedback linearization controller u loc (x k ) [START_REF] Hart | Pyomo: modeling and solving mathematical programs in python[END_REF] for the whole region of |u q (x k )| ≤ U q , ∀q ∈ {x, y, z} [START_REF] Wächter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF]. Hereinafter, we will construct an ellipsoidal invariant set within this region, thus, being able to employ a simple analysis resulted from the linear system (24). Let us consider a symmetric positive definite matrix M ∈ R 6×6 . We obtain the matrix P ∈ R 6×6 , also being symmetric positive definite, as the unique solution of the following Lyapunov equation:

A cl P A cl = P -M. ( 32
)
with the Schur matrix A cl from (24). Proposition 2: Let us consider the set X f defined as follows:

X f = {x ∈ R 6 | x 2 P ≤ δ}, (33)
with δ given by: δ = λ min (P )r 2 , with r 2 = min q∈{x,y,z}

U 2 q K 2 1q + K 2 2q , (34)
and U q from (27)-(29). Then, for all x k ∈ X f , we have that:

(i) the feedback linearization controller u loc (x k ) ∈ U (22),
for all ψ ∈ [-π, π];

(ii) X f is forward invariant 5 under u loc (x k ). Furthermore, x k+1 satisfies:

x k+1 2 P ≤ ρ x k 2 P , s.t. x k+1 = f (x k , u loc (x k )), (35)
with ρ ∈ (0, 1) given by: ρ = 1λ min (M )/λ max (P ).

(36)

Proof: First, since the matrix P obtained from (32) is symmetric positive definite, for all x k such that x k 2 P ≤ δ with δ as in (34), we have that:

λ min (P ) x k 2 ≤ x k 2 P ≤ δ, (37)
thus, x k 2 ≤ r 2 , which implies that q 2 k + v 2 q k ≤ r 2 for all q ∈ {x, y, z} with r as in (34). Hence, by using Cauchy-Schwarz inequality, the virtual inputs u q (x k ) (23) satisfy:

|u q (x k )| ≤ r 2 (K 2 1q + K 2 2q ) ≤ U q . ( 38
)
Then, Proposition 1 provides that u loc (x k ) ∈ U and the system

x k+1 = f (x k , u loc (x k )) is linearized into the linear system x k+1 = A cl x k (24)
. Hence, we arrive to:

x k+1 2 P = x k A cl P A cl x k = x k 2 P -x k 2 M (39) ≤ x k 2 P -λ min (M )/λ max (P ) x k 2 P ρ x k 2 P < δ.
Thus, X f is forward invariant and (35) is validated.

IV. STABILITY OF DIFFERENT NMPC SCHEMES

In the following, we gather all the results obtained from the feedback linearization controller u loc (x k ) [START_REF] Hart | Pyomo: modeling and solving mathematical programs in python[END_REF] to design two NMPC schemes with and without terminal constraint. We check that Assumptions 1 and 2 given in Section II-A hold, hence, guaranteeing the asymptotic stability property for stabilizing the thrust-propelled system (1).

A. NMPC with terminal constraint (NMPC-t)

Proposition 3: Consider the setup of the NMPC-t scheme ( 8)-(10) using the terminal region X f in (33) with the matrix P from the Lyapunov equation (32). If the symmetric matrix M in (32) is chosen such that M Q * with the matrix Q * verifying:

Q * = Q + λ max (R)(K K + 2Γ), (40)
with Q and R from (8), K and Γ from (31), then, Assumption 1 is satisfied with the feedback linearization controller u loc (x k ) from [START_REF] Hart | Pyomo: modeling and solving mathematical programs in python[END_REF]. Thus, the NMPC-t scheme ( 8)-( 10) guarantees the recursive feasibility and the (nominal) asymptotic stability for the closed-loop dynamics [START_REF] Chen | A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability[END_REF]. Proof: At first, since M Q * with Q * as in (40), the matrices M and P obtained from (32) are both symmetric positive definite. Then, Proposition 2 holds and hence, X f is forward invariant under u loc (x k ) (35). Furthermore, since the virtual inputs u q (x k ) ≤ U q , ∀q ∈ {x, y, z} (38) for all x k ∈ X f , 5 A sub-level set of a Lyapunov function is a forward invariant set.

Proposition 1 provides u loc (x k )-u e 2 ≤ (K K +2Γ) with K and Γ from (30)-(31). Thus, the stage cost (x k , u loc (x k )) is bounded by:

(x k , u loc (x k )) ≤ x k Q * x k (41)
with Q * from (40). Then, by employing the terminal cost [START_REF] Bemporad | The explicit linear quadratic regulator for constrained systems[END_REF], the condition ( 13) is verified as follows:

F (x k ) = x k -x e 2 P in
F (f (x k , u loc (x k ))) -F (x k ) + (x k , u loc (x k )) (42) ≤ x k (A cl P A cl -P )x k + x k Q * x k = x k (-M + Q * )x k ≤ 0,
in which P satisfies (32) and M Q * from (40). Hence, completing the proof. Hereinafter, let us summarize the design procedure of the NMPC scheme with terminal constraint given in ( 8)- [START_REF] Munir | Low complexity constrained control using higher degree lyapunov functions[END_REF] for stabilizing the dynamics [START_REF] Zanelli | Nonlinear model predictive control of a human-sized quadrotor[END_REF].

Procedure 1: NMPC-t design 1) Choose the positive definite symmetric matrices Q ∈ R 6×6 and R ∈ R 3×3 to formulate the stage cost in (8). 2) Choose the saturation limits U x , U y and U z satisfying ( 27)-(29). 3) Choose the control gains K 1q , K 2q with q ∈ {x, y, z} satisfying (26). 4) Find the matrix Q * in (40), then, define the symmetric matrix M Q * and solve the Lyapunov equation (32) for P . 5) Find δ in (34) to obtain the terminal region X f in (33). 6) Define the prediction horizon N T such that the first NMPC iteration w.r.t. the initial state x 0 is feasible. Remark 2: For a predefined initial state x 0 , Step 6 in Procedure 1 requires the users to run the NMPC-t algorithm once, while progressively increasing prediction horizon until the optimization problem becomes feasible. However, note that, the computation time is greatly affected by any increase of the prediction length. Thus, in order to continue increasing the region of attraction when the prediction horizon N T is already large, one can increase the size of the terminal constraint set X f given in (33) instead. The parameters which affect most the size of the terminal region X f (33) are the control gains K 1q , K 2q with q ∈ {x, y, z} in (26). In general, decreasing the magnitudes of the control gains increase the size of the set X f and vice versa. Therefore, when state constraints are considered as mentioned in Remark 1, the control gains K 1q , K 2q with q ∈ {x, y, z} are required to be appropriately chosen such that X f ⊆ X with the state constraint set X from [START_REF] Craig | Introduction to robotics: mechanics and control[END_REF].

B. NMPC without terminal constraint (u-NMPC)

Proposition 4: Consider the translational thrust-propelled dynamics (1) and the corresponding u-NMPC setup ( 14)- [START_REF] Grüne | Nonlinear Model Predictive Control: Theory and Algorithms[END_REF]. By choosing the matrix M ∈ R 6×6 used in (32) being symmetric positive definite and setting the constants c, γ ∈ R >0 from [START_REF] Reble | Unconstrained model predictive control and suboptimality estimates for nonlinear continuous-time systems[END_REF] as:

c = δ λ min (Q) λ max (P ) , γ = λ max (Q * )λ max (P ) λ min (Q)λ min (P )(1 -ρ) , (43)
with δ in (34), Q in (8), Q * in (40), ρ in (36) and P satisfying (32), then, Assumption 2 is satisfied. Furthermore, the equilibrium point x e in ( 5) is uniformly exponentially stable under the nominal closed-loop dynamics [START_REF] Grüne | Nonlinear Model Predictive Control: Theory and Algorithms[END_REF] for all the initial states x 0 satisfying V U,N U (x 0 ) ≤ cγ.

Proof: Since the matrices Q as in [START_REF] Gulan | Embedded linear model predictive control for 8bit microcontrollers via convex lifting[END_REF] and P from (32) are symmetric and positive definite, we have the following chain of inequalities:

x k 2 P ≤ λ max (P ) x k 2 ≤ λ max (P ) λ min (Q) x k 2 Q , ( 44
) hence, x k 2 Q ≤ c implies that x k 2 P ≤ δ.
Thus, for all initial states x k such that x k 2 Q ≤ c, Proposition 2 provides a feasible candidate consisting of the state and input trajectories (x loc (•|k), u loc (•|k)), of the prediction horizon length N U , given by:

x loc (i + 1|k) = f (x loc (i|k), u loc (x loc (i|k))), (45)
with i ∈ {0, . . . , N U -1}, the feedback linearization controller u loc from ( 22) and the initial condition x loc (0|k) = x k . Since V U,N U as in ( 14) is the optimal value function, it satisfies:

V U,N U (x k ) ≤ J U,N U (x k , u loc (•|k)). (46)
Furthermore, similar to (41), the stage cost (8) of 16) is bounded as follows:

J U,N U (x k , u loc (•|k)) as in (
(x loc (i|k), u loc (i|k)) ≤ x loc (i|k) 2 Q * , (47)
with Q * in (40). Then, since x k 2 P ≤ δ, x loc (i|k) 2 P exponentially converges as in (35). Thus, introducing (35) to (47) leads to:

(x loc (i|k), u loc (i|k)) ≤ λ max (Q * ) λ min (P ) ρ i x k 2 P ≤ λ max (Q * )λ max (P ) λ min (Q)λ min (P ) ρ i x k 2 Q . ( 48
)
Lastly, combining (46) and (48) implies that:

V U,N U (x k ) ≤ N U -1 i=0 (x loc (i|k), u loc (i|k)) (49)
≤ N U -1 i=0 ρ i λ max (Q * )λ max (P ) λ min (Q)λ min (P ) x k 2 Q ≤ ∞ i=0 ρ i λ max (Q * )λ max (P ) λ min (Q)λ min (P ) x k 2 Q = γ x k 2 Q ,
with γ from (43). Thus, Assumption 2 is validated and the stability of the corresponding u-NMPC design follows Theorem 2. Hence, completing the proof.

In the following, the design procedure of the u-NMPC controller in ( 14)-( 17) is summarized. Procedure 2: u-NMPC design 1) 2) 3) Similar to the first three steps of Procedure 1. 4) Choose the symmetric positive definite matrix M in (32) by following the analysis detailed in Section V and solve the Lyapunov equation (32) for P . 5) Find δ as in (34), ρ as in (35), c and γ as in (43) in order to find N 0 given in [START_REF] Nguyen | A stablizing nmpc design for thrust-propelled vehicles dynamics via feedback linearization[END_REF].

6) Define the prediction horizon N U ≥ N 0 as in [START_REF] Nguyen | A stablizing nmpc design for thrust-propelled vehicles dynamics via feedback linearization[END_REF]. Remark 3: The region of attraction of the u-NMPC design ( 14)-( 17) which consists of all the initial states x 0 satisfying V U,N U (x 0 ) ≤ cγ is difficult to estimate since it does not show an explicit requirement on the initial state x 0 but requires solving the first u-NMPC iteration. Thus, the usual approach in the literature [START_REF] Grüne | NMPC without terminal constraints[END_REF], [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF] is to employ the condition x 0 2

Q ≤ c instead since V U,N U (x 0 ) ≤ γ x 0 2
Q holds as given in (49). However, it is worth noting that the region of x 0 2 Q ≤ c is very conservative w.r.t. the original region of attraction satisfying V U,N U (x 0 ) ≤ cγ.

V. TUNING THE U-NMPC CONTROLLER

It is well-known in the literature that the u-NMPC scheme requires a sophisticated tuning procedure in order to obtain the reasonable values of the required minimum prediction horizon N 0 and also the region of attraction which guarantee the stability [START_REF] Grüne | NMPC without terminal constraints[END_REF], [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF]. However, to the best of our knowledge, the tuning problems of the NMPC without terminal stabilizing constraints have been underestimated in various relating works [START_REF] Grüne | NMPC without terminal constraints[END_REF], [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF], [START_REF] Limón | On the stability of constrained mpc without terminal constraint[END_REF], [START_REF] Reble | Unconstrained model predictive control and suboptimality estimates for nonlinear continuous-time systems[END_REF]. That is to say, people concentrate mostly on the stability proofs of their NMPC designs (e.g. as our contribution in Proposition 4), then, provide one illustrative example with specific parameters [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF]. Theses examples actually aim to show how the results are obtained (e.g. in order for the readers to validate again the calculation process by themselves) but do not give the insight into the actual tuning process. For our particular u-NMPC design ( 14)-( 17), the most influential parameters are the control gains K 1q , K 2q with q ∈ {x, y, z} satisfying (26), the matrix M from (32) and the weighting matrices Q, R from (8) which are not easy to tune due to a large amount of decision variables. In Table I, we show that an arbitrary choice of these parameters can easily lead to an impractically large prediction horizon length while appropriate changes in the tuning parameters allow to reduce the required prediction horizon N 0 and to increase the domain of attraction (represented by p 2 Q ≤ c as in [START_REF] Reble | Unconstrained model predictive control and suboptimality estimates for nonlinear continuous-time systems[END_REF]). In order to analyze how to tune the prediction horizon N 0 given in (20), we will mostly consider the parameter γ as in (43). By introducing ρ from (36) to (43), we can express γ as the multiplication of γ 1 and γ 2 as follows:

γ = λ max (Q * ) λ min (Q) γ1 λ 2 max (P ) λ min (P )λ min (M ) γ2 , (50)
with P obtained from solving (32) with a symmetric positive definite matrix M ∈ R 6×6 , Q in (8) and Q * in (40). At first, from the formulation of N 0 given in [START_REF] Nguyen | A stablizing nmpc design for thrust-propelled vehicles dynamics via feedback linearization[END_REF], decreasing the value of γ as in (50) also reduces the prediction horizon N 0 . Therefore, it is in our interest to minimize the values of γ 1 and γ 2 defined in (50). In order to reduce the decision variables, we follow some reasonable assumptions:

• We give equal importance to the motions along the three axes, i.e.: with q 1 , q 2 > 0 and k 1 , k 2 < 0 satisfying (26). Consequently, the matrix A cl in (24) becomes:

Q = diag{q 1 , q 1 , q 1 , q 2 , q 2 , q 2 }, (51)
K 1x = K 1y = K 1z = k 1 , (52)
K 2x = K 2y = K 2z = k 2 . ( 53
) Description Q R K 1 K 2 M N 0 c First choice I 3 0 3×3 0 3×3 0.
A cl = 1 + ∆ 2 t 2 k 1 I 3 ∆ t 1 + ∆t 2 k 2 I 3 ∆ t k 1 I 3 (1 + ∆ t k 2 )I 3 . (54)
• The matrix M ∈ R 6×6 in (32) is parametrized as follows:

M = m 1 I 3 m 3 I 3 m 3 I 3 m 2 I 3 , (55)
with m 1 = 1 and m 2 , m 3 , m 4 ∈ R satisfying the following conditions:

m 1 m 2 > m 2 3 , (56)
in which ( 56) is to guarantee the positive definiteness of M . With m 1 = 1, (56) is simplified to m 2 > m 2 3 . Remark 4: The reason for fixing m 1 = 1 in (55) is that any scaling on M results on the same scaling on P due to the Lyapunov equation (32). Then, all the eigenvalues of M and P are scaled similarly, hence, γ 1 in (50) does not change and therefore, fixing m 1 = 1 does not cause any loss of generality within the analysis (but reduces the numerical issues). Next, the Lyapunov equation (32) is linear in P which provides a full-rank system of linear equations with A cl and M from (54)-(55). Thus, by solving this, the matrix P is explicitly given as follows:

P = p 1 I 3 p 3 I 3 p 3 I 3 p 2 I 3 , (57)
with p 1 , p 2 , p 3 calculated as:

  p 1 p 2 p 3   =   1 -a 2 1 -2a 1 a 3 -a 2 3 -a 2 2 -2a 2 a 4 1 -a 2 4 -a 1 a 2 1 -a 1 a 4 -a 2 a 3 -a 3 a 4   -1   m 1 m 2 m 3   , (58)
in which, a 1 , . . . , a 4 are taken from A cl a 1 I 3 a 2 I 3 a 3 I 3 a 4 I 3 with A cl in (54). The eigenvalues of P from (57) are given by:

λ max min (P ) = 1 2 p 1 + p 2 ± (p 1 -p 2 ) 2 + 4p 2 3 . (59)
A. Tuning γ 1

By introducing the matrix Q from (51) and the control gains k 1 , k 2 as in ( 52)-( 53) to the formulation of the matrix Q * (40), we obtain the maximum eigenvalue of Q * as follows:

λ max (Q * ) = 1 2 q 1 + q 2 + r(k 2 1 + k 2 2 ) (60) + (q 1 -q 2 + r(k 2 1 -k 2 2 )) 2 + 4r 2 k 2 1 k 2 2 ≤ √ 2(q1-q2) 2 +2r 2 (k 4 1 +k 4 2 )
, where r is directly proportional with the maximum eigenvalue of R:

r = λ max (R) 1 + 2 (-U z + g) 2 , (61)
with g the gravitational acceleration and U z ∈ R + , U z < g from Proposition 1. Then, γ 1 is obtained by introducing (51)-( 60) into (50) and is bounded by:

q 1 + q 2 + r(k 1 + k 2 ) 2 2 min(q 1 , q 2 ) ≤ γ 1 ≤ (62)
q 1 + q 2 + r(k 2 1 + k 2 2 ) + 2(q 1 -q 2 ) 2 + 2r 2 (k 4 1 + k 4 2 ) 2 min(q 1 , q 2 )
, in which, the first inequality is by introducing

q 1 -q 2 + r(k 2 1 -k 2 2 )
2 ≥ 0 to (60) and the latter one is by using the inequality given below (60). From (62), it can be observed that reducing the values of

|q 1 -q 2 |, λ max (R), k 2
1 and k 2 2 probably provide a smaller value of γ 1 (we use "probably" since the reductions of the aforementioned parameters actually make both the upper and lower bounds as in (62) smaller). Especially, in case of employing q 1 = q 2 = q with q ∈ R + , γ 1 from (50) is simplified to γ * 1 as follows:

γ * 1 = 1 + r(k 2 1 + k 2 2 ) q , ( 63
)
which actually allows us to obtain a specific value of γ 1 (i.e., ≥ 1) by tuning only two weighting matrices Q and R from [START_REF] Gulan | Embedded linear model predictive control for 8bit microcontrollers via convex lifting[END_REF] regardless the predefined values of (k 1 , k 2 ) as in ( 52)-(53). Therefore, for tuning γ 1 as in (50), we propose several general directions:

• decrease the ratio λ max (Q)/λ min (Q) as much as possible.

• decrease the value of λ max (R) based on the employed values of (k 1 , k 2 ) but bear in mind that a small value of R causes large input consumption.

• decrease the values of k 2 1 and k 2 2 . However, this is not encouraged since it can cause an unexpected increase in the value of γ 2 defined as in (50) and probably result in a larger prediction horizon N 0 . Remark 5: Formulation (63) explains the effects of the first two tuning steps given in Table I: i) reducing λ max (R) and ii) eliminating the term |q 1 -q 2 |. However, it does not explain for the rest of Table I which requires us to analyze the parameter γ 2 as detailed in the following.

B. Tuning γ 2

Our approach is to construct an explicit formulation of γ 2 from (50) by using the parametrizations of M as in (55) and P as in (59). The function is expressed in terms of the tuning parameters, i.e., the control gains k 1 , k 2 from ( 52)-( 53) and the scalars m 2 , m 3 from (55) as follows:

γ 2 γ 2 (k 1 , k 2 , m 2 , m 3 ). ( 64
)
It will be in our interest to analyze and to find the (local) minimums of the function γ 2 (•) from (64). Since the function is strongly nonlinear, non-convex and contains up to four variables, we have to divide the task into two steps:

1. Find the optimal values of (m 2 , m 3 ) which provide the local minimum value of γ 2 corresponding to a specific values of (k 1 , k 2 ):

γ * 2 (k 1 , k 2 ) = min m2,m3 γ 2 (k 1 , k 2 , m 2 , m 3 ). ( 65
) subject to m 2 > m 2 3 as required in (56), (m 2 , m 3 ) ∈ [m 2min , m 2max ] × [-m 3max , m 3max ],
with (m 2min , m 2max , m 3max ) positive scalars specifying the ranges of m 2 and m 3 as in (55) which are defined accordingly to the condition m 2 > m 2 3 as required in (56) and are chosen by user. The results of this step are illustrated in Figure 1 for (k 1 , k 2 ) = (-0.1, -0.1) and (-1, -1).

Apply

Step 1 for different choices of (k 1 , k 2 ), then, compare the obtained minimum values γ * 2 as in (65) in order to provide the complete analysis as illustrated in Figure 2. Remark 6: The solution of the optimization problem in (65) is obtained by checking a mesh grid of the variables (m 2 , m 3 )

within the specific range [m 2min , m 2max ] × [-m 3max , m 3max ],
hence, the solution's accuracy depends on the resolution of the mesh grid. However, the accuracy problem is not critical due to the fact that we always have to choose the prediction horizon N U ≥ N 0 (c.f. Theorem 2) and N U ∈ N (e.g., it does not matter if we obtain N 0 = 152.4 instead of the precise value assumed to be 152.36 since the minimum prediction horizon needs to be an integer, which leads in both case to N U = 153). Furthermore, using the mesh grid method provides very fast computing times. E.g., it takes only 0.08 seconds to construct the whole data points for Figure 2 which contains 51 × 41 points of (m 2 , m 3 ) and 16 × 16 points of (k 1 , k 2 ). (b) N 0 obtained at the optimal choices of (m 2 , m 3 ), Q = I 6 , R = 0.01I 3 Fig. 2: Illustration of the analysis on N 0 from [START_REF] Nguyen | A stablizing nmpc design for thrust-propelled vehicles dynamics via feedback linearization[END_REF]. ) which provide the smallest N 0 in comparison with the method in [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF] (using Q = I 6 and R = 0.01I 3 as in ( 16)).

k 1 (52) k (53) m 3 (55) m 2 (55) N 0 (20) c (43) -1.3 -1.
We enumerate in Table II all the scenarios in which, we obtain the smallest value of the required prediction horizon N 0 by using the weighting matrices Q = I 6 and R = 0.01I 3 . Table II also shows the corresponding values of c from (43) which helps identifying all the feasible initial states p 0 , i.e., V U,N U (x 0 ) ≤ cγ as stated in Theorem 2. All the choices of (k 1 , k 2 , m 2 , m 3 ) gathered in Table II require the minimum prediction horizon around 153 ∼ 154 but they provide a large range of c from 0.205 to 0.287. Furthermore, using the control gains (k 1 , k 2 ) = (-1.2, -0.8) and M = I 6 provides the shortest prediction horizon N 0 = 153 and the largest value of c = 0.287, hence, being the best choice in our analysis. For comparison, we use the method proposed in [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF] which employs a linear controller as well as the linearized model of the dynamics [START_REF] Zanelli | Nonlinear model predictive control of a human-sized quadrotor[END_REF]. The results are given in the last three lines in which both the value of c and the prediction horizon N 0 are much more conservative than the results of our proposed method. This is due to the restriction of using the linear controller for the nonlinear system (1). This strongly confirms the effectiveness of our NMPC design approach using the local FL controller u loc (x k ) as defined in [START_REF] Hart | Pyomo: modeling and solving mathematical programs in python[END_REF].

We also notice that even the shortest prediction horizon in our analysis N 0 = 153 steps is still extremely large for a real implementation. However, bear in mind that the obtained results (c.f. Table II) are still more promising than employing the linear controller as considered in [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF]. This also indicates a big gap which still exists between the theory on NMPC design and their practical formulations since through various simulations and experimental tests, the u-NMPC controller ( 16)-( 17) requires a prediction horizon of only 10 steps to stabilize the system (1). This is due to the fact that the u-NMPC controller can fully exploit the inputs of the system while a standard controller (e.g. FL controller u loc (x k ) as in [START_REF] Hart | Pyomo: modeling and solving mathematical programs in python[END_REF] or the linear controller as employed in [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF]) enforces the inputs to follow their explicit formulations, hence, the convergence speed of the standard controller (e.g. ρ as in (35)) can not be fast enough to obtain this short prediction horizon of 10 steps. Also note that, better solutions can be found with more thorough parametrizations of the control gains (52)-(53) and of the matrix M (55). However, the trade-off is between increasing the computation time and the complexity of the optimization problem (65) (which may even cause infeasibility).

VI. SIMULATION RESULTS

To validate the presented results, we consider the simulation model (1) of a Crazyflie 2.0 nano-quadcopter platform [START_REF] Budaciu | On the evaluation of the crazyflie modular quadcopter system[END_REF]. It is assumed that the rotation dynamics are stabilized by a fast torque control mechanism [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF]. As specified in (4), the vehicle has to respect thrust limits, with T limit = 2g, and roll and pitch constraints, i.e., c = 10 • in (4)) which lead to the choice of U x = U y = U z = 1.0875 which satisfy the conditions ( 27)-(28) as given in Tables III and IV. We fix the direction angle ψ = 0 and the initial state at x 0 = 0 0.15 -0.1 0 0 0 . Hereinafter, we introduce first the numerical values of the designs of the two NMPC schemes, i.e., with the terminal constraint (denoted by NMPC-t) and without the terminal constraint (u-NMPC) detailed in Sections IV-A and IV-B, respectively. Their performances are analyzed and compared through simulation and experimental results. The NMPC optimization problems are implemented by using Pyomo [START_REF] Hart | Pyomo: modeling and solving mathematical programs in python[END_REF] and solver IPOPT [START_REF] Wächter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF] in Python 3.0. The sampling time is fixed at ∆ t = 0.1 seconds with ∆ t as in (2)-( 3) which means that the NMPC computing time should be less than ∆ t = 0.1 seconds for making possible the application to real experiments.

A. NMPC setup 1) NMPC-t: The chosen parameters of the NMPC-t control scheme (8)-( 10) are gathered into Table III (c.f. Procedure 1). We choose the control gains K 1q = -5, K 2q = -5 for all q ∈ {x, y, z} (26) so that the corresponding region X f is large enough for employing the prediction horizon N T = 5 and the controlled system also has a good convergence speed.

2) u-NMPC: Table IV presents the parameters of the u-NMPC control scheme ( 14)- [START_REF] Limón | On the stability of constrained mpc without terminal constraint[END_REF]. The values given in Table IV are the best choices according to Section V.

B. Simulation results

We provide the simulation results for four scenarios: Scenario 1: Stabilize the dynamics (1) with ψ = 0 using the NMPC-t controller detailed in ( 8)- [START_REF] Munir | Low complexity constrained control using higher degree lyapunov functions[END_REF] with the parameters gathered in Table III. Using the u-NMPC controller given in ( 14)- [START_REF] Limón | On the stability of constrained mpc without terminal constraint[END_REF] with the parameters in Table IV. Scenario 3: Using the u-NMPC controller given in ( 14)- [START_REF] Limón | On the stability of constrained mpc without terminal constraint[END_REF] with the prediction horizon of 10 steps (enough for stabilizing the system through simulations). Scenario 4: Using the local feedback linearization controller u loc given in [START_REF] Hart | Pyomo: modeling and solving mathematical programs in python[END_REF] with the control gains as in Table IV.

The results are given in Figures 3456where they are plotted in red, green, blue and black corresponding to the four aforementioned scenarios, respectively. Figure 3 presents the convergence of the state trajectories under simulation. Note that, only the results of the y and z axes are plotted since the system does not move along the x axis. The terminal region X f ⊂ R 6 given in (33) of the NMPC-t scheme is illustrated by the yellow ellipsoid (which is obtained by slicing X f along three 2D spaces of the position and the corresponding velocity, then, making a convex union of the three obtained 2D sets). All the closed-loop state trajectories under four scenarios converge to the origin but with different convergence speeds as being illustrated with more details in Figure 4 and also given in Table V. Also, the input results are shown in Figure 5 where all the proposed NMPC controllers respect the input constraints (4). We only provide the input results during the first one second for better clarity since they are all monotonously converging to the equilibrium u e (5). Furthermore, Figure 6 and also Fig. 3: Terminal region X f (approximated illustration) and the trajectories (y, v y ), (z, v z ) under different scenarios. From the presented simulation results, the NMPC-t controller employed under Scenario 1 obtains the shortest transient state (1.8 seconds as given in Table V) due to its appropriate terminal constraint set X f (yellow ellipsoid in Figure 3). This is also better than using the feedback linearization controller u loc defined in ( 22) which provides the second shortest convergence time of 2.1 seconds. Next, when not using the terminal constraints, we can see that the u-NMPC controller under Scenario 2 (using the 153-step prediction horizon) obtains the transition time of 3.1 seconds while employing the u-NMPC controller with the prediction horizon of only 10 steps under Scenario 3 makes the states converging in longer time, i.e., 4.4 seconds. It is also trivial that the computation burden when using the 153-step prediction horizon under Scenario 2 is much heavier than those under other scenarios as can be seen from Figure 6. The average computing time of the NMPC controller with the 153-step prediction horizon is 243.7 milliseconds while it is only 47.9 second for the u-NMPC controller with the 10-step prediction horizon, being the best value in our analysis. We also notice that the u-NMPC controller using the 153step prediction horizon (which is found by our analysis given in Section V) is not ready for real implementation since its minimum computing time is already 209.3 milliseconds while the sampling time is fixed at 100 milliseconds. However, the proposed calculation process based on the feedback linearization controller u loc [START_REF] Hart | Pyomo: modeling and solving mathematical programs in python[END_REF] already reduces significantly the required prediction horizon (153 steps) in comparison with the method using the standard linear controller introduced in [START_REF] Köhler | Nonlinear reference tracking: An economic model predictive control perspective[END_REF] which requires thousands of steps to stabilize the system as mentioned in Table II. In the next section, we show that both the NMPC-t and u-NMPC controllers considered in Scenarios 1 and 3 (using the 10-step prediction horizon) are successfully validated through experiments over a real quadcopter platform.

VII. EXPERIMENTAL RESULTS

This section introduces the experimental validation of the two NMPC controllers, NMPC with and without terminal constraint given in Sections IV-A and IV-B over the laboratory Crazyflie 2.0 (CF) nano-quadcopter. The CF has its inner controller running at 500 Hz to track the set-point of the thrust and the three angles sent from the station computer [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF]. In the station computer, we implement the NMPC optimization problems ( 6) and ( 14) with the chosen sampling time ∆ t = 0.1 seconds (fast enough to stabilize the system) by using Pyomo [START_REF] Hart | Pyomo: modeling and solving mathematical programs in python[END_REF] and the IPOPT solver [START_REF] Wächter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF] in Python 3.0. In the following, we describe the limits of the platform and the mismatches between the theoretically nominal NMPC application and real implementation, then, illustrates the obtained results of the two proposed NMPC controllers.

A. Experimental validation limits and how to overcome them

Firstly, as introduced in Section II-A, the yaw angle ψ is assumed to be a known constant. However, maintaining a constant direction angle for an aerial vehicle is obviously impossible for long running times. Therefore, we try to stabilize the yaw angle of the CF around zero and update the actual yaw value to the NMPC optimization problem at each sampling time. This alternative approach still guarantees the nominal stability of the closed-loop scheme with a less strict assumption, i.e., requiring the yaw angle value to be constant only during the prediction horizon 6 .

Secondly, the execution time is always significant when considering an NMPC controller, especially for controlling the Crazyflie 2.0 quadcopter system with the required sampling time of 0.1 seconds (as illustrated in Figure 6 for simulation results). Thus, even with a perfect state feedback at time step k, we cannot obtain the MPC control action immediately at the same time step k as assumed in [START_REF] Munir | Low complexity constrained control using higher degree lyapunov functions[END_REF] and [START_REF] Limón | On the stability of constrained mpc without terminal constraint[END_REF]. Hence, we have to relax this assumption by introducing the obtained MPC input of time step k to the CF system at time step k + 1.

Remark 7: If the reader is interested in using embedded NMPC with low-power hardware, a discussion on how to re-formulate the NMPC problems ( 6) and ( 14) into their approximated quadratic formulations, and thereafter, solve them by using a modified interior-point solver is detailed in [START_REF] Zanelli | Nonlinear model predictive control of a human-sized quadrotor[END_REF]. Another discussion on solving a linear MPC problem with a 8-bit microcontroller by using a convex lifting method is presented in [START_REF] Gulan | Embedded linear model predictive control for 8bit microcontrollers via convex lifting[END_REF].

B. Experimental results

For experiments, we consider the two Scenarios 1 and 3 as given in Section VI with the parameters of the controllers detailed in Tables III and IV (the u-NMPC controller has its prediction horizon of 10 steps). Figure 7 presents the results of x (green lines), y (blue lines) and z (red lines) under two scenarios which shows that both the NMPC-t and u-NMPC controller succeed in stabilizing the CF quadcopter. It can be seen again that the NMPC-t provides faster convergence speed than the u-NMPC controller due to its terminal constraint (1.5 seconds vs. 3 seconds, from blue lines in Figure 7). The terminal constraint also results in more input consumption as clearly observed from Figure 8. The maximum values of the thrust and the roll angle of the NMPC-t (dashed lines) are significantly higher than the ones from the u-NMPC controller (solid lines). All the inputs respect their constraints as given in [START_REF] Nguyen | Effective angular constrained trajectory generation for thrust-propelled vehicles[END_REF]. The computing times under experiment of both NMPC controllers are given in Figure 9 which shows that the NMPC-t controller (with N T = 5 steps and the terminal constraint set X f as in (33)) requires more computation effort than the u-NMPC controller (with N U = 10 steps) but both succeeds in running within the sampling time ∆ t = 0.1 seconds.

VIII. CONCLUSION

This paper presented the application of a nonlinear feedback linearization controller for designing two different NMPC schemes, i.e., NMPC with and without terminal constraint, to stabilize the thrust-propelled vehicles dynamics with asymptotic stability guarantee. By using several standard mathematical inequalities, an ellipsoidal invariant set was constructed under the feedback linerization controller ensuring input constraint satisfaction. Within this invariant set, the feedback law was also bounded by the weighted norm of the corresponding state. These ingredients were gathered to ensure the satisfaction of the stability conditions of the two NMPC designs. The proposed control methods were validated through extensive simulations and experiments over a real quadcopter platform. Future works will analyze the mismatches on angle tracking and the delays on sending NMPC input to the controlled system as employed in the experiments.

2 (Fig. 1 :

 21 Fig. 1: Illustration of γ 2 (k 1 , k 2 , m 2 , m 3 ) as in (64) with different values of (k 1 , k 2 ).

  Illustration of the optimal function γ * 2 (k 1 , k 2 ) defined in (65)

Fig. 4 :φ 1 3 Fig. 5 :

 4135 Fig. 4: Convergences of x, y, z under different scenarios.

  u-NMPC with NU = 153 steps Scenario 3: u-NMPC with NU = 10 steps

Fig. 6 :

 6 Fig. 6: Computation time under simulations.

Fig. 7 :

 7 Fig. 7: Convergences of x, y, z under different scenarios.

Fig. 8 :

 8 Fig. 8: Values of the input under different scenarios.

Fig. 9 :

 9 Fig. 9: Computation time under experiments using Pyomo in Python 3.0.

  then, the followings hold: (i) the feedback linearization control action u loc (x

k ) ∈ U (22), for all value of ψ ∈ [-π, π]; (ii) the values of u loc (x k )u e

2 are bounded by:

TABLE I :

 I Prediction horizon length w.r.t different tuning parameters.

TABLE II :

 II Optimal values of (k 1 , k 2 , m 2 , m 3

TABLE III :

 III Parameters of the NMPC-t controller.

			Values
	Q, R in (8)	I 6 , 0.1I 3
	Ux, Uy, Uz in (27)-(29)	1.0875, 1.0875, 1.0875
	K 1q , K 2q , q ∈ {x, y, z} in (26)	-5, -5
	M in (32) (M	Q * (40))	20I 6
	P in (32)	22.02I 3 10.92I 3 10.92I 3 12.68I 3
	δ in (34)		0.1285
	N T in (6)	5

TABLE IV :

 IV Parameters of the u-NMPC controller.

		Values
	Q, R in (8)	I 6 , 0.01I 3
	Ux, Uy, Uz in (27)-(29)	1.0875, 1.0875, 1.0875
	k 1 , k 2 , in (52)-(53)	-1.2, -0.8
	M in (32)	I 6
	P in (32)	18.22I 3 4.02I 3	4.02I 3 12.22I 3
	δ in (34), ρ in (36)	5.78, 0.95
	c, γ in (43)	0.29, 40.96
	N 0 in (20), N U in (14)	N 0 = 152.22, N U = 153

  Table V present the computing time per step of the three NMPC controllers employed under the first three scenarios.

	Controller's	Convergence	Computing time [ms]
	information	time (95%) [s]	Mean	Min	Max
	Scenario 1	1.8	48.5	31.2	84.6
	Scenario 2	3.1	243.7 209.3 298.4
	Scenario 3	4.4	47.9	31.3	84.6
	Scenario 4	2.1	Not applicable

TABLE V :

 V Comparison between different controllers.

		0.15		(v y0 , y 0 )
		0.1
	Position [m]	0 0.05 -0.05	Xf (vy, y)1 (vy, y)2 (vy, y)3 (vy, y)4	(vz, z)1 (vz, z)2 (vz, z)3 (vz, z)4
		-0.4 -0.1	-0.3	-0.2 Velocity [m/s] (v z0 , z 0 ) -0.1	0	0.1

We have used in (3) "c" and "s" to denote the cos(•) and sin(•) functions, respectively.

Usually referred as quasi-infinite horizon NMPC[START_REF] Chen | A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability[END_REF].

The notation u loc stands for local controller, i.e., being locally admissible.

The initial iteration successfully executed implies the feasibility of all the further steps.

The stability is still guaranteed since all the required ingredients detailed in Section III hold with a general constant yaw angle ψ ∈ [-π, π].

APPENDIX A PROOF OF PROPOSITION 1

By applying our previous results presented in [START_REF] Nguyen | Flat trajectory design and tracking with saturation guarantees: a nano-drone application[END_REF], [START_REF] Nguyen | Effective angular constrained trajectory generation for thrust-propelled vehicles[END_REF] to the nonlinear feedback laws (22b)-(22c), we have the following relations:

where is defined as:

Next, as detailed in Appendix B, we have that:

q∈{x,y,z}

Furthermore, by applying arctan(p) ≤ p, ∀p ≥ 0 to (67) and then, combining this result with (66), we arrive to:

in which the denominator is due to |u z (x k )| ≤ U z < g. A similar result is obtained for θ 2 loc (x k ). Finally, from (68)-(69), we have that:

By considering T loc in (22a), we have that:

Next, if u z + g > 0, we have that:

otherwise, if u z + g ≤ 0, it is straightforward to obtain:

Thus, introducing ( 72) and ( 73) to (71) leads to:

This result holds regardless the sign of u z + g.