Ali Srour
email: ali.srour@irisa.fr

Antonio Franchi
email: a.franchi@utwente.nl

Paolo Robuffo Giordano

Controller and Trajectory Optimization for a Quadrotor UAV with Parametric Uncertainty

In this work, we exploit the recent notion of closed-loop state sensitivity to critically compare three typical controllers for a quadrotor UAV with the goal of evaluating the impact of controller choice, gain tuning and shape of the reference trajectory in minimizing the sensitivity of the closedloop system against uncertainties in the model parameters. To this end, we propose a novel optimization problem that takes into account both the shape of the reference trajectory and the controller gains. We then run a large statistical campaign for comparing the performance of the three controllers which provides some interesting insight for the goal of increasing closed-loop robustness against parametric uncertainties.

I. INTRODUCTION

The requirement for autonomous/intelligent systems such as robots to operate in unpredictable and/or uncertain realworld conditions is one of the biggest challenges in our community. In fact, whether directly or indirectly, robot decisions and control actions are based on a model of the "world", which is oftentimes only a rough representation of the physical reality. When considering motion tasks, one of the common sources of uncertainty is in the robot model parameters that can be hard to measure or also vary at runtime. Adaptive [START_REF] Astrom | Adaptive control 2nd edition[END_REF] or robust control [START_REF] Zhou | Essentials of robust control[END_REF] techniques are typical ways to deal with parametric uncertainty by either trying to estimate the parameters online, or by trading off performance for robustness vs. parametric uncertainty. Also, Model Predictive Control [START_REF] Sun | A comparative study of nonlinear mpc and differential-flatness-based control for quadrotor agile flight[END_REF], [START_REF] Hanover | Performance, precision, and payloads: Adaptive nonlinear mpc for quadrotors[END_REF] predicts system behavior and optimizes a cost function over a finite time horizon using a dynamic model of the system and its environment. However, real-time rescheduling for optimization at each time step can result in significant computational overhead. Another recent approach to deal with these issues has been proposed in [START_REF] Giordano | Trajectory generation for minimum closed-loop state sensitivity[END_REF]- [START_REF] De Pierrepont | Model-in-the-loop testing of control systems and path planner algorithms for quadrotor uavs[END_REF] where several metrics based on the notion of closed-loop state/input sensitivity have been introduced. These metrics are able to quantify how uncertainties in the parameters of a robot model can affect the evolution of the robot states and inputs in closed-loop and as a function of the reference trajectory being tracked (i.e., the 'motion task'). One can then attempt to increase robustness against uncertain parameters by planning a suitable feedforward (or desired) trajectory that minimizes the sensitivity metrics, thereby generating an intrinsically-robust and control-aware motion plan. These ideas have also been extended to the problem of combining state observabilty with parametric robustness in [START_REF] Bohm | Cop: Control & observability-aware planning[END_REF], and in [START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF] a recent extension is proposed for obtaining the 'tubes' of perturbed trajectories given a known range of parametric deviation.

All these previous works have shown very promising results by focusing on the optimization of the reference trajectory to be tracked by the robot under suitable state and input constraints. However, one may also consider the concurrent optimization of the gains of the chosen control strategy (besides the reference trajectory) for the purpose of minimizing a sensitivity metric. In this paper we then consider this latter possibility applied, as case study, to a tracking task for a quadrotor UAV under three popular control strategies: a near-hovering (or linear) controller, a nonlinear geometric controller, and a sliding-mode controller. We present and discuss the results of a thorough statistical analysis aimed at comparing the three controllers in different conditions (non-optimized, by only optimizing for the reference trajectory, by optimizing for the control gains only and by optimizing for both the trajectory and the control gains). The aim is twofold: (i) to show how gain optimization (or tuning) can be done while taking into account a concrete robustness metric (the state sensitivity vs. parametric uncertainties) and (ii) to get statistical insights into the behavior, performance and robustness of the three controllers in different conditions.

The rest of the paper is structured as follows: in Sect. II we recall the main notions related to the closed-loop state sensitivity. Then in Sect. III we present the optimization problems considered in this work, and we detail the quadrotor model and the three tracking controller. Section IV discusses the results of the statistical analysis on the controller performance, and Sect. V concludes the paper and draws some future perspectives.

II. PRELIMINARIES

In this section we briefly summarize for the reader's convenience the main notions related to the closed-loop state sensitivity introduced in [START_REF] Giordano | Trajectory generation for minimum closed-loop state sensitivity[END_REF], [START_REF] Brault | Robust trajectory planning with parametric uncertainties[END_REF] and recently extended in [START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF]. We consider a robot model

ẋ = f (x, u, p) (1)
where x ∈ R nq is the state, u ∈ R nu the control inputs and p ∈ R np a vector of (possibly uncertain) model parameters.

We assume some quantity of interest s(x) ∈ R ns (e.g., the position of the robot center) needs to track a desired trajectory s d (t, a) defined for t ∈ [t 0 , t f ] and function of a finite set of parameters a ∈ R na . For this task a tracking controller is used

ξ = g (ξ, x, a, p c , k c , t) u = h (ξ, x, a, p c , k c , t) , (2)
which is evaluated at a nominal value p c for the uncertain parameters p. The vector ξ ∈ R n ξ represents the possible internal controller states (e.g., an integral action) and k c ∈ R n k is the vector of control gains. The state sensitivity for the closed-loop system (1-2) is defined as

Π(t) = ∂x(t) ∂p p=pc (3)
and the input sensitivity

Θ(t) = ∂u(t) ∂p p=pc . (4)
Matrix Π(t) quantifies how variations in the parameters p around a nominal value p c will affect the evolution of the states x (in closed-loop). Analogously, matrix Θ(t) relates variations in p to variations in the inputs u. A closedform expression for matrixes Π(t) and Θ(t) is in general not available, but in [START_REF] Giordano | Trajectory generation for minimum closed-loop state sensitivity[END_REF], [START_REF] Brault | Robust trajectory planning with parametric uncertainties[END_REF] it is shown how these two quantities can be easily evaluated via forward integration of a differential equation along the system trajectories. Matrixes Π(t) and Θ(t) can be used for several purposes such as optimization of the reference trajectory s d (t, a) for producing minimally sensitive motion plans (by minimizing some combined norm of Π(t) and Θ(t) as done in [START_REF] Giordano | Trajectory generation for minimum closed-loop state sensitivity[END_REF], [START_REF] Brault | Robust trajectory planning with parametric uncertainties[END_REF]). An extension of these ideas has been recently proposed in [START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF] where a suitable weighted norm for the state sensitivity Π(t) is introduced. Assume that each parameter p i can vary in a given range δp i centered at a nominal p ci

p i ∈ [p ci -δp i , p ci + δp i ] (5)
and define the diagonal weight matrix W = diag(δp 2 i ). Letting ∆p = pp c , an ellipsoid in parameter space centered at p c and with semi-axes δp i has equation

∆p T W -1 ∆p = 1. (6)
From ( 6) and [START_REF] Sun | A comparative study of nonlinear mpc and differential-flatness-based control for quadrotor agile flight[END_REF][START_REF] Hanover | Performance, precision, and payloads: Adaptive nonlinear mpc for quadrotors[END_REF] one can obtain the corresponding ellipsoids in state space ∆x T (ΠW Π T ) -1 ∆x = 1 [START_REF] Candido | Minimum Uncertainty Robot Path Planning using a POMDP Approach[END_REF] and in input space

∆u T (ΘW Θ T ) -1 ∆u = 1, (8)
see [START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF]. Here ∆x stands for ∆x = x-x nom where x nom is the state evolution of (1-2) in the unperturbed case p = p c , and analogously for ∆u = uu nom .

The state and input space ellipsoids can be used for different purposes. First of all, one can define a 'weighted sensitivity norm' by considering the eigenvalues λ i of the kernel matrix ΠW Π T in [START_REF] Candido | Minimum Uncertainty Robot Path Planning using a POMDP Approach[END_REF]. In particular, in [START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF] and in this work we consider the following matrix norm

∥Π∥ W = max(λ i (ΠW Π T )) (9)
which represents the largest (worst-case) deviation of the states x assuming a parametric uncertainty as in [START_REF] Giordano | Trajectory generation for minimum closed-loop state sensitivity[END_REF]. Furthermore, one can also exploit [START_REF] Candido | Minimum Uncertainty Robot Path Planning using a POMDP Approach[END_REF][START_REF] Ansari | Minimum Sensitivity Control for Planning with Parametric and Hybrid Uncertainty[END_REF] for obtaining the tubes of perturbed trajectories for the individual components of the states and the inputs. By again referring to [START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF] for all details, for each direction of interest in the state space one can obtain the 'tube radius' r i (t) such that q nom,i (t) -r i (t) ≤ q i (t) ≤ q nom,i (t) + r i (t).

where, as usual, q nom,i (t) is the behavior of the state q i (t) in the unperturbed case p = p c . Equation (10) bounds from above/below the envelope of perturbed states when the parameter uncertainty is bounded as in [START_REF] Giordano | Trajectory generation for minimum closed-loop state sensitivity[END_REF], and an analogous upper/lower bound can also be obtained for the inputs components u i (t).

III. OPTIMIZATION PROBLEM

As explained in Sect. I, in this work we are interested in comparing the performance of three typical trajectory tracking controllers for a quadrotor UAV in the context of minimally sensitive trajectory optimization. To this end, and considering the various quantities introduced in Sect. II, we consider a first trajectory optimization problem

a * = arg min a ∥Π(t f )∥ W s.t. M a = b U min,i ≤ u nom,i (t) -r i (t) ∀i ∀t ∈ [t 0 , t f ] u nom,i (t) + r i (t) ≤ U max,i ∀i ∀t ∈ [t 0 , t f ]. (11)
We seek to find the optimal value a * for the shape parameter a of the reference trajectory s d (a, , t) for minimizing the weighted norm (9) at the final time t f . The constraints consist of given initial/final conditions for s d (a, , t), represented by the linear constraints M a = b, and constraints that bound the envelope of perturbed inputs within actuation limits U min,i ≤ U max,i , ensuring that the tracking of the optimized reference trajectory will be feasible for any value of the uncertain parameters p in the range [START_REF] Giordano | Trajectory generation for minimum closed-loop state sensitivity[END_REF]. These constraints leverage the 'tubes' as in [START_REF] De Pierrepont | Model-in-the-loop testing of control systems and path planner algorithms for quadrotor uavs[END_REF] but are evaluated for the inputs.

Variants of this problem have already been considered in [START_REF] Giordano | Trajectory generation for minimum closed-loop state sensitivity[END_REF], [START_REF] Brault | Robust trajectory planning with parametric uncertainties[END_REF], [START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF] but by always focusing on a single control strategy. We also consider a second optimization problem

k * c = arg min kc ∥Π(t f )∥ W s.t. M a = b K min,i ≤ k ci ≤ K max,i U min,i ≤ u nom,i (t) -r i (t) ∀i ∀t ∈ [t 0 , t f ] u nom,i (t) + r i (t) ≤ U max,i ∀i ∀t ∈ [t 0 , t f ] (12)
where we optimize w.r.t. the control gains k c inside suitable bounds 0 < K min,i ≤ K max,i . We then consider a third optimization problem that optimizes both the shape parameter a and the control gains k c subject to the same constraints as in equation ( 12).

(a * , k * c ) = arg min a, kc ∥Π(t f )∥ W s.t. M a = b K min,i ≤ k ci ≤ K max,i U min,i ≤ u nom,i (t) -r i (t) ∀i ∀t ∈ [t 0 , t f ] u nom,i (t) + r i (t) ≤ U max,i ∀i ∀t ∈ [t 0 , t f ] (13)
We note that in the previous problems we focus on the minimization of the state sensitivity norm at the final time t f . In other words, we look for the combination of reference trajectory and/or control gains that will result in the largest tracking accuracy at the final time (for, e.g., reaching at best a target location). Of course, different choices are also possible such as minimization of (the integral) of the sensitivity norm along the whole trajectory (for increasing the average tracking accuracy during motion). Note also that additional constraints, such as obstacle avoidance, could be easily added to the optimization problem by exploiting the tubes on the states, see [START_REF] Wasiela | A Sensitivity-Aware Motion Planner (SAMP) to Generate Intrinsically-Robust Trajectories[END_REF] for an example in this sense. One can also optimize the sensitivity at specific waypoints along the trajectory where, for instance, the effects of uncertainties in high speed regimes can be minimized.

A. Quadrotor Model

Let x = (r, v, q, ω) ∈ R 6 × S 3 × R 3 be the quadrotor state vector consisting of the world-frame position r = (r x , r y , r z ) and velocity v = (v x , v y , v z ), the rotation from body to world frame represented by a unit quaternion q = (q w , q x , q y , q z ) and the body-frame angular velocity ω = (ω x , ω y , ω z ). The quadrotor model is

                     ṙ = v v = f m   2(qwqy + qxqz) 2(qyqz -qwqx) 1 -2(q 2 x + q 2 y )   -g q = 1 2 0 w ⊗ q ẇ = I -1 (τ -w × Iw) (14)
In this model, f and τ denote the total thrust and bodyframe moment, g = (0, 0, g) is the gravity vector in world frame, and m and I represent the mass and bodyframe inertia. Although quaternions are used to represent the quadrotor orientation in model ( 14), the conversion of quaternions to the roll ϕ, pitch θ, and yaw ψ angles is also considered in two of the control strategies.

Let w i be the squared velocity of the i-th propeller and define the quadrotor control input u = (w 1 , , . . . , , w 4 ). An allocation matrix is used to relate the inputs u (i.e., the squared propeller speeds) to the thrust/torques (f, , τ ).

f τ = k f    1 1 1 1 -gy L -gy -gy -(L + gy) -(L -gx) gx L + gx gx kt -kt kt -kt    u = = T u (15)
The aerodynamic coefficients of the propellers are denoted by k f and k t , the x and y coordinates of the barycenter in the body frame are denoted by g x and g y , and the length of the propeller arms is denoted by L (see [START_REF] Antonelli | Adaptive trajectory tracking for quadrotor mavs in presence of parameter uncertainties and external disturbances[END_REF] for a more detailed derivation). These coefficients are difficult to measure reliably and may change at runtime depending on the flight regime. Therefore, the vector of uncertain parameters considered in this work is p = (k f , k t , g x , g y ).

B. Near-Hovering Controller (NH)

The first control strategy considered is the near-hovering or linear controller, commonly used for its ease of implementation and tuning. However, it only performs well and maintains stability near a hovering state with small roll/pitch angles [START_REF] Mahony | Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor[END_REF]. Tracking is done using a 4D reference trajectory s d (t) = (r d (t), ψ d (t)), which includes the quadrotor 3D position r d (t) = (x d (t), y d (t), z d (t)) and yaw angle ψ d (t). The final control equations are

                                                                 ẍr = rdx + k dx ( ṙdx -ẋ) + kp x (r dx -x) +ki x (r dx -x)dt ÿr = rdy + k dy ( ṙdy -ẏ) + kp y (r dy -y) +ki y (r dy -y)dt ϕ d = -arcsin m f (-sin(ψ)ẍr + cos(ψ)ÿr) θ d = arcsin m f cos(ϕ) (cos(ψ)ẍr + sin(ψ)ÿr) f = m cos(ϕ) cos(θ) (g + rdz + k dz ( ṙdz -ż) + kp z (r dz -z) +ki z (r dz -z)dt τx = -k d ϕ ωx + kp ϕ (ϕ d -ϕ) τy = -k d θ ωy + kp θ (θ d -θ) τz = -k d ψ ωz + kp ψ (ψ d -ψ) . (16)
Finally, the quadrotor inputs are computed as

u = T | -1 pc f τ ( 17
)
where the the allocation matrix T from ( 15) is evaluated at the nominal parameters p c . The vector of control gains k c considered in problem [START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF][START_REF] Wasiela | A Sensitivity-Aware Motion Planner (SAMP) to Generate Intrinsically-Robust Trajectories[END_REF] for this controller consists of all the 15 gains in [START_REF] Lee | Geometric Tracking Control of a Quadrotor UAV on SE(3)[END_REF], that is, the three gains (k p• , k d• , k i• ) for the x, y and z channels, and the two gains (k p• , k d• ) for the ϕ, θ and ψ channels.

C. 3D Geometric Controller (LEE)

The 3D Geometric controller [START_REF] Lee | Geometric Tracking Control of a Quadrotor UAV on SE(3)[END_REF] is a popular control strategy that differs from the near-hovering controller by not relying on approximations of the quadrotor dynamics and providing almost global stability. However, it is more challenging to tune than the near-hovering controller.

Let R(q) ∈ SO(3) be the rotation matrix associated to the quaternion q, and let e 3 = (0, 0, 1). The idea behind this control strategy is to design a desired thrust vector f R d e 3 for stabilizing the translational dynamics (i.e., tracking the reference r d (t)) while the torques τ try to align the quadrotor z axis Re 3 with the desired one and also track the reference signal for the yaw angle ψ d (t). The final controller equations are (see [START_REF] Lee | Geometric Tracking Control of a Quadrotor UAV on SE(3)[END_REF])

                                   er = r d -r, ev = ṙd -v b3 d = Krer + Kvev + Ki erdt + mge3 + mr d ∥Krer + Kvev + Ki erdt + mge3 + mr d ∥ b1 d = (cos ψ d , sin ψ d , 0), b2 d = (b3 d × b1 d )/∥b3 d × b1 d ∥ R d = [b2 d × b3 d b2 d b 3d ] eR = 1 2 (R T d R -R T R d ) ∨ , eω = ω -R T R d ω d f = Krer + Kvev + Ki erdt + mge3 + mr d T Re3 τ = -KReR -Kωeω . ( 18
)
where () ∨ is the usual 'vee' map taking a skew-symmetric matrix into the associated 3D vector, and K r , K v , K i , K R , K ω are 3 × 3 diagonal control gain matrixes. We note that, compared to the original [START_REF] Lee | Geometric Tracking Control of a Quadrotor UAV on SE(3)[END_REF], the moment equation (last row of ( 18)) neglects the compensation for the gyroscopic term as customary in actual implementations of this controller. Note also that, as in the previous nearhovering case, the actual quadrotor inputs u are retrieved by plugging in [START_REF] Bensalah | Full modelling and sliding mode control for a quadrotor uav in visual servoing task[END_REF] the thrust/torques (f, τ ) computed from [START_REF] Labbadi | Robust integral terminal sliding mode control for quadrotor uav with external disturbances[END_REF].

The vector of control gains k c considered in problem (12-13) for this controller consists of all the 15 gains in [START_REF] Labbadi | Robust integral terminal sliding mode control for quadrotor uav with external disturbances[END_REF]

(the diagonal entries of K r , K v , K i , K R , K ω ).

D. Sliding Mode Controller (SMC)

The last controller considered in this work is based on sliding mode, which is a control technique known to generally deal well with model uncertain parameters [START_REF] Bensalah | Full modelling and sliding mode control for a quadrotor uav in visual servoing task[END_REF]- [START_REF] Eltayeb | Sliding mode control design for the attitude and altitude of the quadrotor uav[END_REF]. The controller is built by following the typical steps of the sliding mode design. We start by considering the vertical z dynamics and define the error e z = r dz -z and the sliding surface σ z = ėz + λ z e z where λ z > 0 is a control gain. The Lyapunov candidate V z = 1/2 σ 2 z has time derivative along the system trajectories

Vz = σ z σz = σ z rdz - f m cos(ϕ) cos(θ) + g + λ z ėz . (19)
We then choose the thrust f as

f = m (cos(ϕ) cos(θ)) (g + rdz + k 1z tanh(σ z ) + k 2z σ z ) (20)
where the tanh(•) provides a smooth approximation of the classical sign(•) used in sliding mode. Following analogous arguments, one can define the torque commands as

     τ x = k 1ϕ tanh(σ ϕ ) + k 2ϕ σ ϕ τ y = k 1θ tanh(σ θ ) + k 2θ σ θ τ z = k 1ψ tanh(σ ψ ) + k 2ψ σ ψ (21)
where σ ϕ = ėϕ +λ ϕ e ϕ , σ θ = ėθ +λ θ e θ and σ ψ = ėψ +λ ψ e ψ are the sliding surfaces, and e ϕ = ϕ d -ϕ, e θ = θ d -θ and e ψ = ψ d -ψ. The desired roll/pitch angles are defined as in ( 16) but with

ẍr = rdx + k 1x tanh σ x + k 2x σ x ÿr = rdy + k 1y tanh σ y + k 2y σ y (22)
Here, again σ x = ėx + λ x e x and e x = r dx -r x , and analogously for σ y . The vector of control gains k c considered in problem (12-13) consists of all the 18 gains in [START_REF] Zhou | Bézier curve based smooth path planning for mobile robot[END_REF][START_REF] Mellinger | Minimum Snap Trajectory Generation and Control for Quadrotors[END_REF][START_REF] Loquercio | Learning high-speed flight in the wild[END_REF], that is, the 'sliding' gains λ • and the two gains (k 1• , k 2• ) for the x, y, z, ϕ, θ, and ψ channels.

IV. STATISTICAL ANALYSIS

We now report the results of a statistical analysis for the three controllers of Sect. III when considering the three optimization problems [START_REF] Bohm | Cop: Control & observability-aware planning[END_REF][START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF][START_REF] Wasiela | A Sensitivity-Aware Motion Planner (SAMP) to Generate Intrinsically-Robust Trajectories[END_REF]. For the analysis we generate N traj = 25 initial trajectories s d (a, t) starting at the origin and coping with the initial/final state constraints and input saturations as in [START_REF] Bohm | Cop: Control & observability-aware planning[END_REF][START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF][START_REF] Wasiela | A Sensitivity-Aware Motion Planner (SAMP) to Generate Intrinsically-Robust Trajectories[END_REF]. These initial trajectories are restto-rest motions (from a hovering state to a hovering state) with a final position randomly generated inside a spherical shell of 4 to 6 m centered at the origin, and a final yaw angle randomly generated in the interval [-π/2, π/2]. The trajectories s d (a, t) are implemented as piecewise Bezier curves as already done in [START_REF] Giordano | Trajectory generation for minimum closed-loop state sensitivity[END_REF], [START_REF] Brault | Robust trajectory planning with parametric uncertainties[END_REF], [START_REF] Zhou | Bézier curve based smooth path planning for mobile robot[END_REF] with the goal of achieving a smooth curvature with little snap, as in [START_REF] Mellinger | Minimum Snap Trajectory Generation and Control for Quadrotors[END_REF], [START_REF] Loquercio | Learning high-speed flight in the wild[END_REF]. The framework is implemented in Python and utilizes the COBYLA [START_REF] Conn | On the convergence of derivative-free methods for unconstrained optimization[END_REF] nonlinear optimizer from the nlopt toolbox, along with the symbolic toolbox SymEngine to represent the system symbolically. Additionally, the Jitcode [START_REF] Ansmann | Efficiently and easily integrating differential equations with jitcode, jitcdde, and jitcsde[END_REF] framework is employed for in-time compilation of ordinary differential equations. In the given framework, the NH and LEE controllers have optimization timeframes of 2-3 minutes for each trajectory type, while the SMC controller requires 3-4 minutes per trajectory type.

For the sake of exposition, in the following we will let IN IT represent an initial (non-optimized) trajectory, OP T a the corresponding optimized trajectory solution of ( 11), OP T k the corresponding optimized trajectory solution of ( 12) and OP T ak the corresponding optimized trajectory solution of [START_REF] Wasiela | A Sensitivity-Aware Motion Planner (SAMP) to Generate Intrinsically-Robust Trajectories[END_REF]. We thus obtain a total of 4N traj trajectories: starting from the origin and ending at the same target location, with the parameters p being randomly drawn from [START_REF] Brault | Robust trajectory planning with parametric uncertainties[END_REF]. The final positions r(t f ) of the perturbed trajectories are denoted by small red spheres at t = t f , forming a point cloud that constitutes the empirical variance ellipsoid centered at r d (t f ).

x [m ]
INIT OP T a OP T k OP T ak 0.00 0.03 0.06 0.09 0.12 0.15

r = max(λ i (ΠW Π T )) in [m]
Largest Radius of Uncertainity Ellipsoid at t = t f Fig. 3. Average largest radius of the uncertainty ellipsoids at t = t f , calculated from the combined trajectories of all three controllers, thus reflecting the impact of the optimizations in reducing the ellipsoid size.

the N traj initial ones and the 3N traj optimized ones. Each of these 4N traj trajectories is then provided as reference to the three quadrotor controllers in N sim = 30 simulation runs where in each run the uncertain parameters p are uniformly sampled to lie inside volume of the ellipsoid (6). The uncertainty ranges δp i for building matrix W are set at 15% of the nominal parameter values p c . The control gain bounds in [START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF][START_REF] Wasiela | A Sensitivity-Aware Motion Planner (SAMP) to Generate Intrinsically-Robust Trajectories[END_REF] are set to be 50% and 200% of the initial values used for the INIT case. Therefore problems [START_REF] Bohm | Cop: Control & observability-aware planning[END_REF][START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF][START_REF] Wasiela | A Sensitivity-Aware Motion Planner (SAMP) to Generate Intrinsically-Robust Trajectories[END_REF] will aim at minimizing the sensitivity of the quadrotor position r(t f ) at the final time t f against variations in the parameters k f , k t , g x , g y . For a fist assessment, in Fig. 2 Motor Speed u 0 (t) three corresponding OP T a (2 nd row), OP T k (3 rd row) and OP T ak (4 th row) when considering the three controllers of the previous section on every column. The plots reports different target reach r(t f ) after each perturbed simulation (in green) forming the empirical variance ellipsoid (in red) of the final quadrotor position r(t f ) centered at the desired r d (t f ). By referring to Fig. 2 and Fig. 3, one can note how the ellipsoids are reduced in the OP T a , OP T k cases and even further in the OP T ak cases w.r.t. the INIT cases as expected. The largest radius of the state ellipsoid is decreased on average for the OP T a and OP T k from (12 cm) to (9 and 7.5 cm) respectively but further decreases in the OP T ak to (3 cm). This then confirms that the sensitivity cost minimized in [START_REF] Bohm | Cop: Control & observability-aware planning[END_REF][START_REF] Brault | Tube-based trajectory optimization for robots with parametric uncertainty[END_REF][START_REF] Wasiela | A Sensitivity-Aware Motion Planner (SAMP) to Generate Intrinsically-Robust Trajectories[END_REF] captures well how variations in the parameters affect deviations in the considered states. Furthermore, one may note that in the optimized cases OP T a and OP T ak where the shape parameter a of the reference trajectory is optimized, the optimal trajectory performs an initial maneuver for then arriving at the target location with an almost vertical straight line motion (see Fig. 1 and the attached video for more details). This trend has been empirically verified in almost all the simulation runs and is shared by all the employed controllers. A likely explanation is that a final (almost) vertical motion allows to better minimize the effects of the uncertainties in the considered parameters since the rotational dynamics is not excited during this phase, contrarily to the INIT trajectory where the quadrotor needs instead to perform a non-negligible tilting for reaching r d (t f ) at rest. This insight can then be helpful for designing more robust hover-to-hover motions also in different conditions. This trend can also be deduced from Table . I that shows the percentage contribution of each parameter (N p = 4) in p on the sensitivity of the states r(t f ) defined by C i,j % in [START_REF] Conn | On the convergence of derivative-free methods for unconstrained optimization[END_REF] (see [START_REF] Hafezipour | An uncertainty analysis method for error reduction in end-effector of spatial robots with joint clearances and link dimension deviations[END_REF]). This analysis can effectively assess the impact of each parameter variation on the final target reach error at time t f

C i,j % = Π 2 ij (t f )δ 2 pj Np j=1 Π 2 ij (t f )δ 2 pj × 100 (23)
Table . I shows that changes in the parameters g x and g y have the most significant impact on the sensitivity of states x and y, while state z is predominantly influenced by k f . Consequently, when attempting to reach the final target, a vertical motion is primarily affected only by uncertainties in k f , while the other three parameters affect the states in non-vertical motions. Furthermore, based on Table I, a 15% uncertainty in k t does not affect the states r(t f ) and can then be disregarded. Figure 4 reports the actuator speeds (in revolutions per minute) for the OP T ak case. The red line represents the nominal actuator speed, while the green lines represent the actuation speeds in the N sim perturbed runs. The dashed black lines denote the limit bounds of the actuator, and the blue line represents the input bound tubes, which ensures actuation feasibility when p varies in the ellipsoid [START_REF] Brault | Robust trajectory planning with parametric uncertainties[END_REF].

Fig. 5 reports the so-called violin plots of the norm of the final target reach (defined as ∥r d (t f ) -r(t f )∥) for the complete statistical campaign considering the three controllers and all the 4N traj trajectories in the N sim runs. The interquartile range is indicated by the thick black bar at the center, the median of the samples is shown by the white dot at the center, and the remaining portions of the distribution are indicated by the thin black line. Greater probabilities of the population adopting the specified value are represented by wider areas of the violin plot. In order to obtain a fair comparison of all controllers, we tuned their initial gains so as to have a comparable median in the INIT case, which is around 8 cm as can be seen in Fig. 5.

Observations show that all controllers demonstrate an enhanced average target reach in the OP T a , OP T k , and OP T ak cases. In the OP T a case, the NH and LEE controllers achieve similar results, with an average target reach of around 5.5 to 6 cm, while the SMC controller improvement is comparatively lower. This suggests that the desired trajectory shape is less critical for the SMC performance, but more relevant for the LEE and NH controllers. In the OP T k case the NH controller average target reach is 7 cm, the LEE controller improvement is about 5 cm, and the SMC controller exhibits the most significant improvement of 4 cm. These results imply that the LEE and SMC controllers are more susceptible to controller gain tuning compared to the NH controller. Moreover, in the OP T ak case, all controllers show significant improvement in both the average target reach and its variance. The median of the NH, LEE, and SMC controllers drops from about 8 cm in the IN IT case to around 2 cm for OP T ak case (400% improvement), with a better distribution around the median.

Therefore, from this analysis we can conclude that the SMC controller is not particularly sensitive to the shape of the reference trajectory, but it is quite sensitive to the control gains which, if tuned well, can lead to good performance. The NH and LEE controllers are instead more affected by the trajectory shape in addition to the control gains. Finally, we can note that the NH and LEE controllers perform quite similarly when optimized. Our most likely explanation is that the controllers are evaluated in terms of the target reach ∥r d (t f ) -r(t f )∥ at the final time, and not in terms of how well they are tracking the reference trajectory during motion. The LEE controller should be better at dealing with maneuvers involving, e.g., large tilting or accelerations, but it can have a comparable performance to the NH one when it comes to reach a final hovering state.

V. CONCLUSIONS

In this work we have proposed a novel optimization problem in which the effect of parametric uncertainties in a robot model (quantified by the notion of state sensitivity) can be minimized by acting on the reference trajectory to be tracked and/or on the control gains. A tracking task for a quadrotor UAV has been chosen as case study, by employing three popular tracking controllers widely used in the community. The results of a large statistical analysis have shown the effectiveness of the proposed joint optimization of the reference trajectory and gains in minimizing the effects of uncertainties in selected quadrotor parameters. In addition to identifying the most influential parameters on the system states, we can also gain valuable insights into controller performance and gain tuning robustness concerning parametric uncertainty.

Fig. 1 .

 1 Fig. 1. Drone trajectory tracking under parametric uncertainties (green) to reach a target s d (t f ) (light blue sphere). Non-optimized (IN IT , top left) and optimized (OP T ak , top right) trajectories are compared. The size of the uncertainty ellipses are compared in the bottom left, where the blue and red ellipses correspond to IN IT and OP T ak , respectively. A video demonstration is available at: https://youtu.be/VBRL8XDiJ5c.

Fig. 2 .

 2 Fig.2. Trajectory tracking results for NH, LEE, and SMC controllers are shown in the left, center, and right columns respectively, for the IN IT , OP Ta, OP T k , and OP T ak cases. Each run has N sim trajectories (green), starting from the origin and ending at the same target location, with the parameters p being randomly drawn from[START_REF] Brault | Robust trajectory planning with parametric uncertainties[END_REF]. The final positions r(t f ) of the perturbed trajectories are denoted by small red spheres at t = t f , forming a point cloud that constitutes the empirical variance ellipsoid centered at r d (t f ).

Fig. 4 .

 4 Fig. 4. (Red) of one control input with the nominal parameters (p = pc) in the OP T ak case. The corresponding uncertainty tube is shown in blue, and the green lines represent the perturbed inputs when p ̸ = pc for N sim runs. The allowable input range is depicted by (dashed black) lines.

INITFig. 5 .

 5 Fig. 5. Violin plots displaying the target reach (defined as ∥r d (t f ) -r(t f )∥) in [m] for the different trajectory types using the three controllers NH, LEE and SMC. These results corresponds to the statistical campaign of trajectories (N traj = 25) and perturbed simulations (N sim = 30).

  we report the N sim 'perturbed runs' for one INIT trajectory (1 st row) and the Percentage Contribution of each Parameter p i on the sensitivity of the states q i

	Param\State kf kt gx gy	x 2.46% 0.00 39.28% 58.27%	y 1.03% 0.00 58.91% 40.06%	z 98.36% 0.00 0.69% 0.95%

Average

TABLE I PERCENTAGE

 I CONTRIBUTION OF EACH PARAMETER IN p ON THE FINAL SENSITIVITY OF THE STATES r(t f )

	u 0 (t)	2 6 9 8 7 5 4 3	0 ×10 3	1 bounds	P = Pc	2	OP T ak tube	3	P = Pc	4	5
							t[s]

A. Srour and P. Robuffo Giordano are with CNRS, Univ Rennes, Inria, IRISA, Campus de Beaulieu, 35042 Rennes Cedex, France. email:

We are now working on an experimental validation of these findings and, on a more methodological side, we are also interested in finding more principled ways to impose bounds on the control gains in problems (12-13) besides simple 'box constraints' (e.g., by explicitly considering the stability of the closed-loop via a Lyapunov function).

This work was supported by the project ANR-20-CE33-0003 "CAMP" and Horizon 2020 research and innovation programme [grant agreement No. 871479] AERIAL-CORE